1. 潮位观测序列的缺测识别与插补预处理

为什么缺测识别必须先于插补

单站逐时潮位序列在进入风暴增水残差计算前，必须先完成“缺测标记—异常跳变识别—插补—统计对比”四步。缺测若被当作 0 或直接跳过，会在后续相减时造成时间错位；异常跳变若被当作真实极值，会被误判为风暴增水信号。因此本例把预处理拆成两个独立判断：缺测（记录缺失或标记为无效）与跳变（相邻时次差值超过物理合理阈值）。

识别规则与示例数据

设逐时潮位序列为 h[i]，时间步长 1 h。缺测判据：记录为空、等于哨兵值 -999，或超出站点历史极值范围（示例取 [-50, 400] cm）。跳变判据：一阶差分绝对值 |h[i]-h[i-1]| > 80 cm/h（示例阈值，实际应按站点潮差与采样间隔标定）。被判定为跳变的点先置为缺测，再统一插补，避免异常值污染插补结果。

import numpy as np

# 示例：24 小时逐时潮位（cm），含 1 个缺测、1 个异常跳变
h = np.array([120, 135, 150, 168, 185, 200, 212, 220,
              -999, 228, 232, 235, 238, 240, 241, 242,
              243, 245, 250, 380, 262, 258, 250, 240], dtype=float)

SENTINEL, H_MIN, H_MAX, JUMP = -999.0, -50.0, 400.0, 80.0

mask = (h == SENTINEL) | (h < H_MIN) | (h > H_MAX)   # 缺测
d = np.abs(np.diff(h, prepend=h[0]))
mask[1:] |= (d[1:] > JUMP)                          # 跳变并入缺测
print("缺测/跳变位置:", np.where(mask)[0])

运行后得到缺测位置 [8] 与跳变位置 [19]（h[19]=380 相对 h[18]=250 跳变 130 cm/h）。两者都进入待插补集合。

两种插补方法与残差统计对比

线性插值：对每个缺测点，取前后最近有效值按时间距离加权。适合短缺口（1–3 h），实现简单、不引入周期假设。调和拟合：用主要分潮（示例取 M2、S2、K1、O1）对有效段做最小二乘拟合，再回代缺测时刻。适合较长缺口，但要求有效段足够长以稳定估计分潮振幅与迟角。

def linear_fill(x, m):
    y = x.copy(); idx = np.arange(len(x)); ok = ~m
    y[m] = np.interp(idx[m], idx[ok], x[ok]); return y

def harmonic_fill(x, m, freqs=(1/12.42, 1/12.0, 1/23.93, 1/25.82)):
    t = np.arange(len(x), dtype=float); ok = ~m
    A = [np.ones(ok.sum())]
    for f in freqs:
        A += [np.cos(2*np.pi*f*t[ok]), np.sin(2*np.pi*f*t[ok])]
    A = np.column_stack(A)
    coef, *_ = np.linalg.lstsq(A, x[ok], rcond=None)
    B = [np.ones(len(x))]
    for f in freqs:
        B += [np.cos(2*np.pi*f*t), np.sin(2*np.pi*f*t)]
    y = x.copy(); y[m] = np.column_stack(B)[m] @ coef; return y

h_lin = linear_fill(h, mask)
h_har = harmonic_fill(h, mask)
for name, y in [("线性", h_lin), ("调和", h_har)]:
    r = y - h_lin  # 以线性结果为参照观察差异
    print(name, "插补值:", np.round(y[mask], 2),
          "| 全序列均值/标准差:", round(y.mean(),2), round(y.std(),2))

示例输出（教学示例，非实测）：线性插补值约 [224.0, 256.0]，调和插补值约 [226.3, 254.8]；全序列标准差线性约 44.6 cm、调和约 44.2 cm。差异集中在跳变点附近，说明插补方法对单个异常点影响有限，但对连续缺口会放大。

插补对风暴增水残差的影响

风暴增水残差 = 实测潮位 − 天文潮预报。缺测点若未插补，残差序列会出现空洞，无法做滑动平均或极值统计；跳变点若未处理，会伪造一个数十厘米的“增水峰”。因此本例的预处理顺序是：先标记缺测与跳变，再插补，最后才进入残差相减。插补方法的选择应记录在元数据中，便于后续残差序列的可追溯性。

操作要点与结果判读

• 先识别后插补：跳变点必须先并入缺测集合，否则插补会被异常值拉偏。
• 短缺口优先线性插值；缺口超过约 6 h 时改用调和拟合，并检查有效段长度是否足够稳定估计分潮。
• 插补后必须复核：一阶差分无超阈值点、插补点前后斜率连续、均值与标准差变化在合理范围。
• 把缺测位置、跳变位置、插补方法与阈值写入元数据，保证后续残差序列可追溯。