2. 基于最小二乘的四大主分潮参数初算

为什么选这四个分潮、设计矩阵长什么样

在完整逐时潮位序列上做最小二乘，本质是把“观测潮位”分解为若干余弦分潮的叠加。对 M2、S2、K1、O1 四个主分潮，任意时刻 t 的模型写成：

h(t) = Z₀ + Σᵢ [ Aᵢ·cos(ωᵢt − gᵢ) ] = Z₀ + Σᵢ [ aᵢ·cos(ωᵢt) + bᵢ·sin(ωᵢt) ]

其中 aᵢ = Aᵢ·cos gᵢ，bᵢ = Aᵢ·sin gᵢ，振幅 Aᵢ = √(aᵢ²+bᵢ²)，迟角 gᵢ = atan2(bᵢ, aᵢ)。这样每个分潮贡献两列（cos 列与 sin 列），加上常数项 Z₀ 一列，共 9 个未知量。设计矩阵 X 的行对应每个观测时刻，列对应 9 个基函数；观测向量 y 为去均值前的原始潮位（常数项由 Z₀ 列吸收）。

四个分潮的角速度（度/小时）取标准值：M2 = 28.9841042，S2 = 30.0000000，K1 = 15.0410686，O1 = 13.9430356。时间基准必须与例1统一：以序列起点为 t=0，单位小时，且已做时区归一到 UTC 或统一的地方时，否则迟角会整体平移。

构造与求解的关键实现

下面给出可直接运行的 Python 实现。输入为两列：t_hours（相对起点的逐时小时数，长度 N）与 eta（潮位，米）。要求 N ≥ 9 且无缺测；若有缺测，先按例1的处理日志剔除或插补后再进入本步。

import numpy as np

# 四个主分潮角速度（度/小时）
OMEGA = {
    "M2": 28.9841042,
    "S2": 30.0000000,
    "K1": 15.0410686,
    "O1": 13.9430356,
}
ORDER = ["M2", "S2", "K1", "O1"]

def build_design_matrix(t_hours):
    t = np.asarray(t_hours, dtype=float)
    cols = [np.ones_like(t)]              # Z0 常数项
    for name in ORDER:
        w = np.deg2rad(OMEGA[name])       # 转弧度/小时
        cols.append(np.cos(w * t))
        cols.append(np.sin(w * t))
    return np.column_stack(cols)          # 形状 (N, 9)

def solve_four_constituents(t_hours, eta):
    t = np.asarray(t_hours, dtype=float)
    y = np.asarray(eta, dtype=float)
    X = build_design_matrix(t)
    coef, *_ = np.linalg.lstsq(X, y, rcond=None)
    y_hat = X @ coef
    resid = y - y_hat
    Z0 = coef[0]
    result = {}
    for i, name in enumerate(ORDER):
        a = coef[1 + 2*i]
        b = coef[2 + 2*i]
        A = np.hypot(a, b)
        g = np.rad2deg(np.arctan2(b, a)) % 360.0
        result[name] = {"a": a, "b": b, "A": A, "g": g}
    rms = float(np.sqrt(np.mean(resid**2)))
    return {"Z0": Z0, "constituents": result, "rms": rms,
            "resid": resid, "y_hat": y_hat, "coef": coef}

要点说明：np.linalg.lstsq 用 SVD 求解，对 9 列设计矩阵在 N 远大于 9 时数值稳定；rcond=None 使用机器精度阈值。迟角用 atan2(b, a) 再对 360° 取模，保证落在 [0, 360)。若序列不足约 15 天，M2 与 S2 的频率差（约 1.016 度/小时）在一个记录内只完成约 0.4 个拍周期，两列接近共线，振幅会互相“借能量”，此时应把结果标注为初算值而非最终值。

残差量级怎么判读

四主分潮只能解释潮汐的主要部分，浅水分潮（M4、M6）、长周期分潮与气象扰动都会进入残差。合理判读不是看残差是否为零，而是看：残差 RMS 相对观测振幅的占比、残差是否仍有明显周期性、以及残差是否存在系统性漂移。经验上，若观测潮差约 2～4 m，四主分潮拟合的残差 RMS 常在数厘米到十几厘米量级；若残差 RMS 接近或超过主分潮振幅的 20%，通常提示时间基准错位、单位不一致、缺测被当成 0 填充，或序列过短导致分潮混叠。

一次可复现的验证流程

用合成数据自检，可以确认实现与判读逻辑一致：取 30 天逐时序列，设定 Z₀=1.80 m，M2 振幅 1.10 m、迟角 120°，S2 振幅 0.35 m、迟角 150°，K1 振幅 0.28 m、迟角 60°，O1 振幅 0.22 m、迟角 30°，叠加 0.03 m 的随机噪声，再调用上面的函数。预期恢复的振幅与设定值相差在噪声量级内，迟角相差在几度以内，残差 RMS 约等于噪声标准差。若把时间单位误写成“天”而非“小时”，四个分潮的迟角会整体错乱、残差 RMS 显著抬升，这正是判读环节要抓的典型错误。