geo-toolbox 水文插件算法解析：从 Muskingum 到 Muskingum-Cunge

本文分析 Miku196/geo-toolbox 中 geo-plugin-hydro 的河道洪水演进模块，讨论经典 Muskingum 法与物理参数化的 Muskingum-Cunge 法，并给出可复现的参考实现。代码位于 plugins/geo-plugin-hydro/src/routing/。

代码性质说明：文中 Rust 代码为基于标准文献的参考实现示意，并非仓库逐行摘录。实际实现请以仓库源码为准。

一、Muskingum 法：用槽蓄方程描述洪水波

本节及下一节的参考实现位于 geo-plugin-hydro/src/routing/。

Muskingum 法由 McCarthy 于 1938 年提出，将河段蓄量 \(S\) 表达为入流 \(I\) 与出流 \(O\) 的线性加权组合：

\[S_t = K\left[X I_t + (1 - X) O_t\right] \]

其中 \(K\) 是洪水波通过河段的传播时间（量纲为时间），\(X\) 是入流对槽蓄影响的权重因子，经典 Muskingum 要求 \(0 \le X \le 0.5\)。\(X = 0\) 时蓄量与出流高度相关，退化为线性水库；\(X = 0.5\) 时入流与出流等权，洪水波只平移不衰减。

将槽蓄方程与水量平衡方程联解，得到经典的演算递推式：

\[O_2 = C_0 I_2 + C_1 I_1 + C_2 O_1 \]

三个系数由 \(K\)、\(X\) 和时段 \(\Delta t\) 决定，满足 \(C_0 + C_1 + C_2 = 1\)：

\[C_0 = \frac{\Delta t - 2KX}{2K(1-X) + \Delta t}, \quad C_1 = \frac{\Delta t + 2KX}{2K(1-X) + \Delta t}, \quad C_2 = \frac{2K(1-X) - \Delta t}{2K(1-X) + \Delta t} \]

graph LR I["入流过程 I(t)"] --> M["Muskingum 演算 O2 = C0·I2 + C1·I1 + C2·O1"] K["K：传播时间"] --> M X["X：权重因子"] --> M DT["Δt：演算时段"] --> M M --> O["出流过程 O(t)"] style M fill:#e0f2fe,stroke:#38bdf8,color:#075985 style O fill:#dcfce7,stroke:#4ade80,color:#065f46

1.1 参数 K 与 X 的确定

经典 Muskingum 法的 \(K\) 和 \(X\) 依赖实测入流-出流资料，用试错法、最小二乘法或矩法率定。这限制了它在无资料流域的应用。

1.2 稳定性条件与分段演算

为避免出流出现负值或振荡，演算时段需满足非负系数条件：

\[2KX \le \Delta t \le 2K(1-X) \]

这与 \(C_0 \ge 0\) 且 \(C_2 \ge 0\) 严格等价。经典 Muskingum 的 \(C_1\) 在 \(K>0\)、\(X \ge 0\)、\(\Delta t > 0\) 下恒正，无需额外检查。

当 \(K\) 远大于 \(\Delta t\) 时，将河段分为 \(N\) 个子段，每段用相同的 Muskingum 法逐段演进，\(N \approx K / \Delta t\)。子段参数取 \(K_{\text{sub}} = K / N\)，\(X_{\text{sub}} = X\)——\(X\) 是河道特性，与河段长度无关，不随分段变化。若按原 \(K\)、\(X\) 逐段演进，会得到完全不同的结果。

1.3 Rust 实现

/// Muskingum 演算系数
///
/// 返回 `(C0, C1, C2)`，满足 C0 + C1 + C2 = 1.0。
///
/// # Panics
/// - `k_hours <= 0`、`delta_t <= 0` 或 `x` 不在 `[0, 0.5]` 内时 panic。
///
/// # Debug 断言
/// 稳定性条件 `2KX <= Δt <= 2K(1-X)`（等价于 `C0 >= 0 && C2 >= 0`）
/// 不满足时，debug 构建下触发 `debug_assert`；release 构建不检查。
/// 注意：此处的 `C1` 在 `K > 0, X >= 0, Δt > 0` 下恒正，无需检查。
pub fn muskingum_coefficients(k_hours: f64, x: f64, delta_t: f64)
    -> (f64, f64, f64)
{
    assert!(k_hours > 0.0, "k_hours 必须为正，当前 {}", k_hours);
    assert!(delta_t > 0.0, "delta_t 必须为正，当前 {}", delta_t);
    assert!((0.0..=0.5).contains(&x), "x 应在 [0, 0.5] 内，当前 {}", x);

    let denom = 2.0 * k_hours * (1.0 - x) + delta_t;
    // k_hours > 0, x <= 0.5, delta_t > 0 => denom > 0，无需额外判断

    let c0 = (delta_t - 2.0 * k_hours * x) / denom;
    let c1 = (delta_t + 2.0 * k_hours * x) / denom;
    let c2 = (2.0 * k_hours * (1.0 - x) - delta_t) / denom;

    debug_assert!(
        c0 >= 0.0 && c2 >= 0.0,
        "稳定性条件不满足：需 2KX <= Δt <= 2K(1-X)，当前 C0={}, C2={}", c0, c2
    );

    (c0, c1, c2)
}

/// 用给定系数逐时段递推
///
/// **初始条件约定**：第一个时段的入流为上一步入流，即
/// `I[-1] = inflow[0]`；出流为上一步出流 `O[-1] = initial_outflow`。
/// 因此第一步实际计算的是
/// `O[0] = (C0 + C1) * I[0] + C2 * O_init`。
/// 若 `initial_outflow == inflow[0]`，则 `O[0] == inflow[0]`，自洽；
/// 否则会产生初始瞬变。
///
/// 空 `inflow` 返回空 `Vec`。
fn route_with_coefficients(
    inflow: &[f64],
    c0: f64, c1: f64, c2: f64,
    initial_outflow: f64,
) -> Vec<f64> {
    let n = inflow.len();
    let mut outflow = Vec::with_capacity(n);
    if n == 0 { return outflow; }

    let mut prev_i = inflow[0];
    let mut prev_o = initial_outflow;

    for &i in inflow {
        let o = c0 * i + c1 * prev_i + c2 * prev_o;
        outflow.push(o);
        prev_i = i;
        prev_o = o;
    }

    outflow
}

/// Muskingum 洪水演进
///
/// `inflow` 为入流序列 (m³/s)，`k_hours` 和 `x` 为 Muskingum 参数，
/// `delta_t` 为时段步长（小时），`initial_outflow` 为初始出流 (m³/s)。
///
/// 初始条件约定见 `route_with_coefficients`。
///
/// # Panics
/// 参数不满足 `muskingum_coefficients` 的约束时 panic。
pub fn muskingum_route(
    inflow: &[f64],
    k_hours: f64,
    x: f64,
    delta_t: f64,
    initial_outflow: f64,
) -> Vec<f64> {
    let (c0, c1, c2) = muskingum_coefficients(k_hours, x, delta_t);
    route_with_coefficients(inflow, c0, c1, c2, initial_outflow)
}

二、Muskingum-Cunge：用河道物理参数替代率定

Muskingum-Cunge 法由 Cunge 于 1969 年提出，核心改进是用河道几何（宽度、坡度）和水力参数（波速、曼宁系数）直接推导 \(K\) 和 \(X\)，摆脱对实测流量资料的依赖。

Cunge 发现：Muskingum 的差分格式在特定参数取值下，其数值扩散可以匹配扩散波方程的物理扩散，从而近似扩散波解。

2.1 物理参数化公式

传播时间 \(K\) 由河段长度 \(\Delta x\) 和洪水波速 \(c\) 决定：

\[K = \frac{\Delta x}{c} \]

权重因子 \(X\) 由 Cell Reynolds 数 \(D\) 导出：

\[X = 0.5 \left(1 - D\right), \quad D = \frac{q_0}{S_0 \cdot c \cdot \Delta x} \]

其中 \(q_0\) 是参考单宽流量，\(S_0\) 是河床坡度。当 \(D = 0\)（\(q_0 = 0\) 或 \(\Delta x \to \infty\)）时 \(X = 0.5\)，无数值扩散；\(D\) 越大，\(X\) 越小，衰减越强。

注意：MC 中 \(X\) 不再受经典 Muskingum 的 \(0 \le X \le 0.5\) 约束。当 \(D > 0.5\) 时 \(X < 0\)，表示更强的数值扩散；当 \(D > 1\) 时 \(X < -0.5\)。这是 MC 与经典 Muskingum 的重要差异。

2.2 波速的 Kleitz-Seddon 定律

洪水波速 \(c\) 由 Kleitz-Seddon 定律给出（Kleitz 1877 与 Seddon 1900 独立提出）：

\[c = \frac{dQ}{dA}\bigg|_{Q_0} \]

对宽浅矩形河道，曼宁公式 \(Q = \frac{1}{n} A R^{2/3} S_0^{1/2}\) 下，\(c \approx \frac{5}{3} v\)（\(v\) 为断面平均流速）。

2.3 演算系数与 Courant 数

Cunge (1969) 将演算系数表达为 Courant 数 \(C\) 和 Cell Reynolds 数 \(D\) 的函数；Ponce & Yevjevich (1978) 将其推广到变参数情形。

\[C = \frac{c \cdot \Delta t}{\Delta x}, \quad D = \frac{q_0}{S_0 \cdot c \cdot \Delta x} \]

\[C_0 = \frac{-1 + C + D}{1 + C + D}, \quad C_1 = \frac{1 + C - D}{1 + C + D}, \quad C_2 = \frac{1 - C + D}{1 + C + D} \]

分母 \(1 + C + D\) 在 \(C \ge 0, D \ge 0\) 下恒正。

稳定性条件：MC 系数非负要求三个分子均非负：

\(C_0 \ge 0 \iff C + D \ge 1 \iff C \ge 1 - D\)
\(C_1 \ge 0 \iff C - D \ge -1 \iff C \ge D - 1\)
\(C_2 \ge 0 \iff C \le 1 + D\)

综合得：

\[|D - 1| \le C \le 1 + D \]

下界来自 \(C_0\) 与 \(C_1\) 两个条件的较大者，即 \(\max(1-D, D-1) = |D-1|\)；上界来自 \(C_2\)。

分情况：

\(0 \le D \le 1\)：下界为 \(1 - D\)，即 \(1 - D \le C \le 1 + D\)；
\(D > 1\)：下界为 \(D - 1\)，即 \(D - 1 \le C \le 1 + D\)。

\(C_1\) 的约束在 \(D > 1\) 时不可忽略。反例：\(D = 2\)，\(C = 0.5\)，此时 \(C_0 \approx 0.43 \ge 0\)、\(C_2 \approx 0.71 \ge 0\) 均满足，但 \(C_1 \approx -0.14 < 0\)，输出会振荡。工程实践中常取 \(C \approx 1\)（Courant 条件）以最小化数值扩散。

graph LR Geo["河道几何 Δx · S0"] --> MC["Muskingum-Cunge 参数计算"] Hyd["水力参数 n · q0 · c"] --> MC MC --> K["K = Δx / c"] MC --> X["X = 0.5(1 - D)"] K --> Coeff["演算系数 C0, C1, C2"] X --> Coeff Coeff --> Route["逐时段演算"] style MC fill:#e0f2fe,stroke:#38bdf8,color:#075985 style Coeff fill:#fef3c7,stroke:#fcd34d,color:#92400e

上图中 \(D\) 由单宽流量、坡度、波速、河段长度计算；宽度隐含在单宽流量 \(q_0\) 中，不是显式参数。

2.4 Rust 实现

MuskingumCungeParams 只保留路由阶段真正消费的字段。\(K\) 和 \(X\) 是 courant 与 cell_reynolds 的等价表示，二者同时存储会有一致性风险（修改 courant 后 k_hours 不会同步），改为按需求值的只读方法。

/// Muskingum-Cunge 参数
///
/// 只保留路由阶段真正消费的字段。
/// `K` 和 `X` 是 `courant` 与 `cell_reynolds` 的等价表示，
/// 通过方法按需推导，避免结构体内部失去一致性。
#[derive(Debug, Clone)]
pub struct MuskingumCungeParams {
    /// 洪水波速 (m/s)
    pub celerity_ms: f64,
    /// Courant 数 C = c · Δt / Δx
    pub courant: f64,
    /// Cell Reynolds 数 D = q0 / (S0 · c · Δx)
    pub cell_reynolds: f64,
}

impl MuskingumCungeParams {
    /// 传播时间 K (小时)，由 `K = Δx / c / 3600` 推导
    ///
    /// 调用方需提供与构造时一致的 `reach_length_m`。
    /// 本方法与结构体内部的 `courant` 存在恒等关系
    /// `courant = Δt_hours / k_hours`，即 `k_hours = Δt_hours / courant`；
    /// 若需从 `courant` 直接反推 K，调用方也可自行按此式计算。
    /// 保留 `reach_length_m` 参数是为了让本方法可独立于构造过程使用。
    pub fn k_hours(&self, reach_length_m: f64) -> f64 {
        reach_length_m / self.celerity_ms / 3600.0
    }

    /// 权重因子 X = 0.5 · (1 - D)
    ///
    /// 注意：MC 中 X 可为负，不受经典 Muskingum 的 `[0, 0.5]` 约束。
    /// 当 D 很大时 X 可远小于 -1。本实现不做 clamp，返回原始值；
    /// 若业务上需限制，可在上层自行 clamp 到 `[-1, 0.5]` 或 `[-0.5, 0.5]`。
    pub fn x(&self) -> f64 {
        0.5 * (1.0 - self.cell_reynolds)
    }
}

/// 由河道物理参数计算 Muskingum-Cunge 参数
///
/// - `reach_length_m`：河段长度 Δx (m)，须 > 0
/// - `slope`：河床坡度 S0 (m/m)，须 > 0
/// - `manning_n`：曼宁系数，须 > 0
/// - `unit_width_q`：参考单宽流量 q0 (m²/s)，须 > 0
/// - `delta_t_hours`：演算时段 (h)，须 > 0
///
/// 由恒等式 `courant = delta_t_hours / k_hours` 可知，
/// `courant` 与 `K` 不是独立量。
///
/// # Panics
/// 任一参数非正时 panic。
pub fn muskingum_cunge_params(
    reach_length_m: f64,
    slope: f64,
    manning_n: f64,
    unit_width_q: f64,
    delta_t_hours: f64,
) -> MuskingumCungeParams {
    assert!(reach_length_m > 0.0, "reach_length_m 必须为正");
    assert!(slope > 0.0, "slope 必须为正");
    assert!(manning_n > 0.0, "manning_n 必须为正");
    assert!(unit_width_q > 0.0, "unit_width_q 必须为正");
    assert!(delta_t_hours > 0.0, "delta_t_hours 必须为正");

    // 正常水深：由曼宁公式反算 (宽浅矩形简化)
    // q = (1/n) * d^(5/3) * S0^(1/2)  =>  d = (q*n / S0^0.5)^(3/5)
    let d = (unit_width_q * manning_n / slope.sqrt()).powf(0.6);
    // 断面平均流速
    let v = unit_width_q / d;
    // Kleitz-Seddon 波速 (宽浅矩形 c ≈ 5/3 v)
    let c = 5.0 / 3.0 * v;

    let delta_t_s = delta_t_hours * 3600.0;
    let courant = c * delta_t_s / reach_length_m;
    let cell_reynolds = unit_width_q / (slope * c * reach_length_m);

    MuskingumCungeParams { celerity_ms: c, courant, cell_reynolds }
}

/// Muskingum-Cunge 洪水演进
///
/// 系数由 `params.courant` 与 `params.cell_reynolds` 决定，
/// 逐时段用公共递推内核演算。
///
/// **参数时变**：本函数接收固定的 `params`，是定参数版本。
/// 变参数 MC 需在每个时段重算 `muskingum_cunge_params` 并逐时段
/// 更新系数。
///
/// 初始条件约定见 `route_with_coefficients`。
///
/// # Debug 断言
/// 系数非负条件 `|D-1| <= C <= 1+D` 不满足时，
/// debug 构建下触发 `debug_assert`；release 构建不检查。
pub fn muskingum_cunge_route(
    inflow: &[f64],
    params: &MuskingumCungeParams,
    initial_outflow: f64,
) -> Vec<f64> {
    let c = params.courant;
    let d = params.cell_reynolds;
    let denom = 1.0 + c + d;

    let c0 = (-1.0 + c + d) / denom;
    let c1 = (1.0 + c - d) / denom;
    let c2 = (1.0 - c + d) / denom;

    // 三个系数均需非负：C0 约束下界 1-D，C1 约束下界 D-1，
    // C2 约束上界 1+D。合起来是 |D-1| <= C <= 1+D。
    debug_assert!(
        c0 >= 0.0 && c1 >= 0.0 && c2 >= 0.0,
        "MC 稳定性条件不满足：需 |D-1| <= C <= 1+D，当前 C={}, D={}, C0={}, C1={}, C2={}",
        c, d, c0, c1, c2
    );

    route_with_coefficients(inflow, c0, c1, c2, initial_outflow)
}

来源：Cunge, J.A. (1969). On the subject of a flood propagation computation method (Muskingum method). Journal of Hydraulic Research, 7(2), 205-230.

三、两种方法的工程对比

graph TD subgraph Classic["Muskingum 法"] C_K["K：由实测资料率定"] C_X["X：由实测资料率定"] C_Lim["依赖有资料流域"] end subgraph MC["Muskingum-Cunge 法"] M_K["K = Δx / c c 由曼宁公式推导"] M_X["X = 0.5(1 - D) D 由单宽流量、坡度、波速、河段长度计算"] M_Lim["适用于无资料流域"] end Data["河道几何 + 曼宁 n"] --> MC Gauge["实测入流-出流资料"] --> Classic style Classic fill:#fef3c7,stroke:#fcd34d style MC fill:#dcfce7,stroke:#86efac

维度	Muskingum	Muskingum-Cunge
参数来源	实测流量资料率定	河道几何 + 曼宁系数
适用场景	有资料流域	无资料流域、河网汇流
参数时变性	外部率定常数	可由河道物理参数直接推导；通过逐时段更新 `q0` 可实现变参数版本
数值扩散	经验性	与物理扩散匹配

四、与生态插件的跨插件衔接

MUSLE 的场次产沙估算需要径流总量 \(Q_{surf}\) 和洪峰流量 \(q_{peak}\)。geo-plugin-hydro 的 muskingum_cunge_route 输出出流过程线（m³/s）后：

取峰值即得 \(q_{peak}\)；
对过程线积分即得 \(Q_{surf}\)，积分时需乘以时段步长（秒）：

\[Q_{surf} = \sum_i O_i \cdot \Delta t_{\text{seconds}} \]

读者若直接对流量序列求和，会得到 m³/s 量纲的数值而非 m³，与 MUSLE 的输入单位不匹配。这一步是实际使用中最容易踩的陷阱。

这种“水文演算 → 产沙估算”的流水线正是 geo-toolbox 分层架构的优势所在。MUSLE 侧的因子计算位于 plugins/geo-plugin-ecology/src/musle.rs。

五、发布说明

公式渲染：需开启 MathJax。
Mermaid 流程图：需引入 mermaid.min.js 并初始化。
代码高亮：rust 在博客园默认不支持，可改为 csharp 或 java。
代码性质：文中 Rust 代码为基于标准文献的参考实现示意，实际仓库实现请以仓库源码为准。

参考

代码仓库：Miku196/geo-toolbox
McCarthy, G.T. (1938). The unit hydrograph and flood routing. US Army Corps of Engineers.
Cunge, J.A. (1969). On the subject of a flood propagation computation method. Journal of Hydraulic Research, 7(2), 205-230.
Ponce, V.M. & Yevjevich, V. (1978). Muskingum-Cunge method with variable parameters. Journal of the Hydraulics Division, 104(12), 1663-1667.
Kleitz, F. (1877). Mémoire sur la théorie du mouvement non permanent des eaux. Annales des Ponts et Chaussées.
Seddon, J.A. (1900). River hydraulics. Transactions of the American Society of Civil Engineers, 43, 179-229.
Chow, V.T. (1959). Open-Channel Hydraulics. McGraw-Hill.
包为民. (2006). 《水文预报》. 中国水利水电出版社.

posted @ 2026-09-22 10:21 mikuyyds 阅读(24) 评论(0) 收藏举报

刷新页面返回顶部

mikuyyds

geo-toolbox 水文插件算法解析：从 Muskingum 到 Muskingum-Cunge

geo-toolbox 水文插件算法解析：从 Muskingum 到 Muskingum-Cunge

一、Muskingum 法：用槽蓄方程描述洪水波

1.1 参数 K 与 X 的确定

1.2 稳定性条件与分段演算

1.3 Rust 实现

二、Muskingum-Cunge：用河道物理参数替代率定

2.1 物理参数化公式

2.2 波速的 Kleitz-Seddon 定律

2.3 演算系数与 Courant 数

2.4 Rust 实现

三、两种方法的工程对比

四、与生态插件的跨插件衔接

五、发布说明

参考

公告