4. 气压效应对残差序列的初步分离

为什么先分离气压分量

风暴增水残差 R(t) 是实测潮位减去天文潮预报后的剩余量，其中同时包含风应力增水、气压效应、波浪爬高与观测噪声。气压效应在物理上方向明确：海面气压每升高 1 hPa，海面约下降 1 cm，即反气压计效应（inverted barometer, IB）。因此可以先用逐时气压观测构造一个可解释的气压增水分量，把它从原始残差中扣除，再看剩余残差的峰值与持续时间如何变化，从而判断后续风应力标定应使用哪一条序列。

本例采用如下符号约定：气压增水 η_p(t) 以“增水为正”定义，即低压引起增水为正：

η_p(t) = -β · (P(t) - P̄) / 100        # P 单位 hPa，β 单位 cm/hPa，结果单位 m
R_raw(t) = 实测潮位 - 天文潮预报
R_corr(t) = R_raw(t) - η_p(t)

其中 P̄ 为分析窗口内的参考气压（本例取窗口均值），β 为反气压计系数，理论值约 1.0 cm/hPa，实际海洋中因水深、陆架宽度与风场耦合常取 0.7～1.0。除以 100 是把 cm 换算为 m，与潮位单位保持一致。

用一组示例数据走完分离过程

取一次冷锋过境过程，窗口内 6 个代表时刻的气压与原始残差如下（示例数据，用于演示计算逻辑，非实测）。参考气压取窗口均值 P̄ = 1008.0 hPa，反气压计系数先取 β = 0.9 cm/hPa。

以 t5 为例：η_p = −0.9 × (−12.0) / 100 = +0.108 m，即低压贡献约 10.8 cm 增水；R_corr = 0.61 − 0.108 = 0.502 m。峰值从 0.61 m 降到 0.502 m，下降 0.108 m，与 β·|ΔP|/100 完全一致——这是分离计算自洽性的第一道校验。

再看持续时间：若以 0.30 m 为超阈值判据，原始残差在 t3～t6 之间均超阈值（约 4 个时段）；修正后 t3 的 0.183 m 已低于阈值，超阈值时段缩短为 t4～t6（约 3 个时段）。这说明气压修正不仅压低峰值，还会改变“风暴增水持续多久”的判断，对后续风应力标定的样本筛选有直接影响。

关键实现代码

下面函数完成时间对齐、参考气压计算、气压分量构造与残差分离，并输出峰值与超阈值时段的对比指标。输入为两个等长的逐时序列，缺测用 NaN 标记。

import numpy as np

def separate_pressure_surge(t, p_hpa, r_raw, beta_cm_per_hpa=0.9,
                            p_ref=None, threshold=0.30):
    """
    t        : 逐时时间（小时，等间隔）
    p_hpa    : 逐时气压 (hPa)，缺测为 NaN
    r_raw    : 原始风暴增水残差 (m)，缺测为 NaN
    beta     : 反气压计系数 (cm/hPa)
    p_ref    : 参考气压；None 时取有效样本均值
    threshold: 超阈值判据 (m)
    返回 dict：eta_p, r_corr, 峰值与持续时间指标
    """
    t = np.asarray(t, dtype=float)
    p = np.asarray(p_hpa, dtype=float)
    r = np.asarray(r_raw, dtype=float)

    valid = np.isfinite(p) & np.isfinite(r)
    if valid.sum() < 10:
        raise ValueError("有效同步样本不足，无法分离气压分量")

    if p_ref is None:
        p_ref = np.nanmean(p[valid])

    eta_p = np.full_like(p, np.nan)
    eta_p[valid] = -beta_cm_per_hpa * (p[valid] - p_ref) / 100.0

    r_corr = np.full_like(r, np.nan)
    r_corr[valid] = r[valid] - eta_p[valid]

    def peak_and_duration(series):
        s = series[valid]
        i_peak = int(np.argmax(s))
        over = s >= threshold
        # 统计连续超阈值的最长时段（小时）
        longest, cur = 0, 0
        for flag in over:
            cur = cur + 1 if flag else 0
            longest = max(longest, cur)
        return s[i_peak], longest

    peak_raw, dur_raw = peak_and_duration(r)
    peak_corr, dur_corr = peak_and_duration(r_corr)

    return {
        "p_ref": p_ref,
        "eta_p": eta_p,
        "r_corr": r_corr,
        "peak_raw": peak_raw, "peak_corr": peak_corr,
        "peak_drop": peak_raw - peak_corr,
        "dur_raw_h": dur_raw, "dur_corr_h": dur_corr,
        "expected_drop_at_peak": beta_cm_per_hpa
            * abs(p[valid][np.argmax(r[valid])] - p_ref) / 100.0,
    }

调用与校验：

res = separate_pressure_surge(t, p_hpa, r_raw, beta_cm_per_hpa=0.9,
                              threshold=0.30)
print(res["p_ref"], res["peak_raw"], res["peak_corr"])
print(res["peak_drop"], res["expected_drop_at_peak"])
print(res["dur_raw_h"], res["dur_corr_h"])

校验要点：peak_drop 与 expected_drop_at_peak 应在数值精度内相等（两者都等于峰值时刻的 β·|ΔP|/100）。若不等，说明峰值时刻的索引在两条序列中不一致，或 valid 掩码使用有误。

β 取值与适用条件

反气压计系数不是常数。开阔大洋深水区接近 1.0；陆架浅水、半封闭海湾因风场与水位耦合，有效 β 常偏低。可用两种方式检验本例的 β 是否合理：

• 峰值一致性检验：用不同 β 重算，观察 peak_drop 是否随 β 线性变化；若某 β 使修正后残差在气压最低时刻出现负值（即“过度修正”），说明 β 偏大。
• 残差-气压回归检验：对 R_raw 与 ΔP 做线性回归，斜率（换算单位后）可作为 β 的经验估计；若回归斜率显著偏离理论值，应记录该站点的有效 β。

结果判读清单

• 峰值下降量 ≈ β·|ΔP_peak|/100，偏差应在舍入误差内。
• 修正后残差不应在气压极值时刻出现不合理负值；若出现，优先怀疑 β 偏大或时间对齐错位。
• 超阈值持续时间的变化方向应与气压过程一致：低压持续越久，修正后持续时间缩短越明显。
• 若修正前后峰值与持续时间几乎不变，说明该过程气压变幅小，气压效应不是主控因子，可保留原始残差进入风应力分析。