基于连续地表水水质自动站的流速与降解系数自动识别算法

摘要

本文提出了一种基于自动站数据的多波峰自动配对与参数反演算法，旨在解决流域水质模型中流速与降解系数难以自动识别的问题。该算法兼顾学术理论与工程应用，能够在复杂水文条件下，自动、批量、稳健地反推出河流区间的流速与降解系数。在重庆 XX 流域的应用中，该算法相比传统方法，模型精度有明显提升。实验结果表明，该方法能够为水环境管理和污染溯源提供有力支撑，并为流域污染物排放控制和应急响应提供决策支持。

1. 研究背景与意义

• 自动站数据价值：地表水自动站可实现高频、连续、时空分布广的水质监测，为流域水质模型，如一维水动力-水质模型，提供了丰富的数据基础。
• 参数识别难点：传统流速和降解系数多依赖人工实验或经验估算，这些方法耗时耗力、主观性强，难以反映动态变化和多事件场景。
• 自动化需求：亟需一种能自动适应不同水文期、自动识别污染事件、批量反演参数的智能算法。

---

2. 算法原理与流程

2.1 一维河流降解模型

污染物在河道中的一阶降解模型：该模型适用于均匀河流，忽略了横向和垂向的浓度差异。

$$

C(x) = C_0 \cdot e^{-\frac{kx}{v}}

$$

• $C(x)$：下游浓度，mg/L;
• $C_0$：上游浓度，mg/L;
• $k$：降解系数 1/d;
• $v$：流速（迁移速度）km/d;
• $x$：区间距离，km。

2.2 多波峰自动配对与参数反演

1. 数据预处理：对自动站数据按站点、因子、时间进行清洗、对齐，并采用线性插值法处理缺失值，对超出合理范围的数值进行剔除。

2. 波峰识别：采用滑动窗口极大值+prominence 算法，自动识别每月的局部波峰，适应不同水文期。其中，prominence 参数表示波峰的显著程度，根据经验，设置为 0.1，以滤除噪声干扰。

3. 相关性配对：对每个上游波峰，在下游自适应时间窗口内，截取波形片段，计算皮尔逊相关系数，筛选相关性大于 0.6 的配对。相关性阈值 0.6 是根据大量实验数据统计分析得到的经验值，可能存在一定的误差，需要根据具体情况进行调整。

4. 参数反演：对每对配对波峰，根据时间差 $\Delta t$ 和浓度变化，自动反推出流速 $v$ 和降解系数 $k$。

   流速 $v$ 的计算公式为：$v = x / \Delta t$

   降解系数 $k$ 的计算公式为：$k = -v \cdot ln(C(x) / C_0) / x$

5. 批量输出：所有站点对、所有因子的参数结果自动拼接汇总，便于后续建模和统计。

---

3. 算法实现与工程细节

• Python 实现：基于 pandas、numpy、scipy、matplotlib 等主流科学计算库，代码结构清晰，易于扩展。
• 健壮性处理：自动跳过无数据或异常数据的站点对与因子组合，保证批量运行稳定。其中，异常数据的判断采用 3σ原则，当数据超出平均值 3 倍标准差范围时，则被认为是异常值。
• 参数自适应：每对站点可单独配置匹配窗口，适应丰/枯水期差异。丰水期采用较短的匹配窗口，枯水期采用较长的匹配窗口，具体窗口长度根据历史数据统计分析确定。
• 可视化辅助：自动生成每月波峰识别与配对结果图，直观展示算法效果，用于人工审核和问题诊断。
• 结果标准化：输出字段包含起点、终点、距离、因子、k、v、时间差、相关性等，便于后续自动化调用。

---

4. 应用案例

以重庆 XX 流域自动站数据为例，选取 XX 站点和 YY 站点之间的数据区间（数据时间范围为 2024 年 1 月 1 日至 2024 年 6 月 30 日），自动识别高锰酸盐指数、总磷等多因子的流速与降解系数。算法可自动适应 180 天历史数据，批量输出多组参数，从输出结果图标看出，峰值关联性识别效果非常好。

自动识别波峰数量

自动批量求 k v 。

5. 优势与创新点

• 全自动化：无需人工干预，适合大规模、长时序数据批量处理。
• 多事件适应：可识别多次污染事件，反映真实水文过程。
• 相关性配对：有效排除偶然波峰，提升参数物理意义，通过排除虚假波峰，避免了将随机扰动误判为污染事件。

---

6. 展望与改进方向

• 更智能的事件识别：可引入机器学习、动态规划等方法，进一步提升波峰配对智能化水平，例如，采用 LSTM 神经网络预测波峰的到达时间，提高配对的准确性。
• 实时在线应用：结合实时数据流，实现参数的动态更新与预警，用于突发水污染事件的预警和应急响应。
• 分析结果，优化滑动区间：结合实际案例，说明如何根据分析结果优化滑动区间。

7. 附录：核心代码片段

# 以站点对和因子为循环，自动批量识别K/V参数
for link in station_links:
    s1, s2, x, default_window = link["起点"], link["终点"], link["距离"], link["默认滑动窗口(h)"]
    for factor in jili_factors:
        # 调用auto_kv_multi_peaks_corr函数，计算流速和降解系数
        kv_df = auto_kv_multi_peaks_corr(
            df, s1, s2, factor, x,
            time_window_dict=None,
            default_time_window=default_window,
            waveform_window=12,
            corr_threshold=0.6
        )
        # ...结果拼接与保存...

本算法已在实际流域水质管理项目中成功应用，欢迎交流与合作！