📚 课程目录第2章基于调和分析的潮汐预报基线构建

2. 基于调和分析的潮汐预报基线构建

为什么需要“调和预报基线”

风暴增水残差 = 实测潮位 − 天文潮预报值。若基线本身偏差 10 cm，残差里就会混入 10 cm 的“假增水”。因此本例的核心不是直接减，而是先用长序列（建议 ≥ 1 年逐时数据）通过最小二乘拟合主要分潮的调和常数，得到一条可外推的天文潮基线，再用不同分潮数目检验基线稳定性。

调和分析的基本模型为：

h(t) = Z0 + Σ_i [ A_i·cos(ω_i·t) + B_i·sin(ω_i·t) ]
A_i = f_i·H_i·cos(g_i)   B_i = f_i·H_i·sin(g_i)
振幅 H_i = sqrt(A_i² + B_i²)   迟角 g_i = atan2(B_i, A_i)

其中 ω_i 为分潮角频率（由天文幅角率给出），Z0 为平均海平面。对每个时刻建立一行设计矩阵，用最小二乘解出 [Z0, A_i, B_i]，即完成调和常数估计。

用最小二乘实现调和常数估计

下面给出可直接运行的 Python 实现。核心是构造设计矩阵：对每个分潮，用其角频率生成 cos、sin 两列，加上常数项，再用 numpy.linalg.lstsq 求解。分潮角频率采用常用近似值（度/小时），实际业务中应替换为天文幅角率表。

import numpy as np

# 主要分潮角频率（度/小时），示例近似值
OMEGA = {
    "M2": 28.9841042, "S2": 30.0000000, "N2": 28.4397295,
    "K2": 30.0821373, "K1": 15.0410686, "O1": 13.9430356,
    "P1": 14.9589314, "Q1": 13.3986609,
    "M4": 57.9682084, "MS4": 58.9841042, "M6": 86.9523127,
}

def build_design(t_hours, constituents):
    cols = [np.ones_like(t_hours)]
    for name in constituents:
        w = np.deg2rad(OMEGA[name])
        cols.append(np.cos(w * t_hours))
        cols.append(np.sin(w * t_hours))
    return np.column_stack(cols)

def fit_constituents(t_hours, h_obs, constituents):
    X = build_design(t_hours, constituents)
    coef, *_ = np.linalg.lstsq(X, h_obs, rcond=None)
    z0 = coef[0]
    amps, phases = {}, {}
    for i, name in enumerate(constituents):
        a, b = coef[1 + 2*i], coef[2 + 2*i]
        amps[name] = np.hypot(a, b)
        phases[name] = np.rad2deg(np.arctan2(b, a)) % 360
    return z0, amps, phases

def predict(t_hours, z0, amps, phases):
    h = np.full_like(t_hours, z0, dtype=float)
    for name, H in amps.items():
        w = np.deg2rad(OMEGA[name])
        g = np.deg2rad(phases[name])
        h += H * np.cos(w * t_hours - g)
    return h

分潮数目对残差量级的影响检验

用同一段长序列分别拟合 4、8、11 分潮，比较拟合残差 RMS 与风暴时段预报的平滑性。示例数据为合成序列（含已知分潮 + 随机噪声），用于演示判读逻辑，非实测结果。

rng = np.random.default_rng(42)
t = np.arange(0, 365*24, 1.0)          # 1 年逐时
true = (1.20
        + 1.10*np.cos(np.deg2rad(28.9841042)*t - np.deg2rad(120))
        + 0.35*np.cos(np.deg2rad(30.0)*t      - np.deg2rad(150))
        + 0.28*np.cos(np.deg2rad(15.0410686)*t - np.deg2rad(60))
        + 0.20*np.cos(np.deg2rad(13.9430356)*t - np.deg2rad(30)))
h_obs = true + rng.normal(0, 0.08, t.size)

plans = {
    "A_4":  ["M2","S2","K1","O1"],
    "B_8":  ["M2","S2","N2","K2","K1","O1","P1","Q1"],
    "C_11": ["M2","S2","N2","K2","K1","O1","P1","Q1","M4","MS4","M6"],
}
for tag, cons in plans.items():
    z0, amps, phs = fit_constituents(t, h_obs, cons)
    res = h_obs - predict(t, z0, amps, phs)
    print(f"{tag:5s} 分潮数={len(cons):2d}  RMS={np.sqrt(np.mean(res**2)):.4f} m")

预期判读：4 分潮时 RMS 较大（漏掉 N2、K2 等造成系统性偏差）；8 分潮 RMS 明显下降；11 分潮相对 8 分潮改善很小（< 5%），说明 8 分潮已基本收敛。此时若继续加潮，只会把噪声拟合进基线，反而在风暴时段产生虚假波动。

选择基线的实用原则：以“增加分潮后 RMS 变化小于 5%”为收敛判据，同时检查风暴时段预报曲线是否平滑。若预报出现与天文潮无关的尖峰，说明分潮过多或序列过短导致过拟合，应回退分潮组合。

调和常数必须来自足够长的序列。少于约 1 年的数据难以分离 K1 与 P1、K2 与 S2 等邻近分潮，容易产生“分潮混叠”，使基线在风暴季节出现系统性偏移。若只有短期数据，应固定使用少量稳健分潮（如 M2、S2、K1、O1）。

基线交付检查表

检查项	通过条件	不通过时的处理
序列长度	≥ 8760 小时（1 年）	缩短分潮列表，仅保留 M2/S2/K1/O1
缺测比例	< 5%，且无长连续缺测	先完成插补，再拟合
拟合残差 RMS	接近观测噪声水平（示例 < 0.15 m）	检查分潮频率、时间基准、单位
分潮收敛性	增加分潮后 RMS 变化 < 5%	停止增加分潮，采用当前组合
风暴时段预报平滑性	无异常尖峰、无负潮位突变	减少分潮，防止过拟合

完成以上步骤后，得到的 z0、amps、phases 即为本站在该时段的调和常数，可用于任意时刻的天文潮预报，作为后续风暴增水残差计算的基线。