马斯京根河道水流演算方法自1938年提出以来^[1]，进行了大量的相关研究，获得了一系列理论和应用成果，是迄今为止应用最广泛、获得理论研究成果最多的河道水流演算方法^[2-4].马斯京根法的理论研究主要有模型检验研究、模型结构的物理基础分析论证研究、模型结构改进研究、模型参数率定与时变关系研究四类.

模型结构检验研究主要是对实际河道上下断面洪水期的流量观测资料，用马斯京根法模型据上断面流量计算下断面流量，通过与实测下断面流量的比较，分析误差大小和误差特征，进而评判模型结构的合理性与模拟实际洪水的效果^[5-9].

模型结构的物理基础分析论证主要是对蓄泄关系的物理基础和马斯京根法演算方程与简化圣维南方程组的差分模型间的关系，分析模型结构的物理基础，间接证明模型的物理性^[10-13].模型结构改进研究主要是对线性蓄泄关系进行非线性化的改进和结构的进一步简化关系研究，以考虑河段蓄量与出流流量间的非线性关系^[14-15].模型参数率定与时变关系研究主要是参数率定方法研究和考虑流量比重系数X和传播时间K的时间变化影响因素，建立参数随水力特征因素变化的函数关系，应用于模型模拟洪水的计算中^[16-18].

所有的马斯京根法研究和应用成果都限制在有稳定水位流量关系线的河段，对于流量关系线不稳定或不存在水位流量关系的潮汐河段，马斯京根流量演算法模拟实际洪水效果就不好或根本无法使用^[19-22].

本文根据河段蓄量与过水断面和流量蓄泄关系，与水流连续关系，构成过水断面积与流量双变量耦合演算模型，以解决流量关系线不稳定或不存在水位流量关系的潮汐河段洪水演算问题.

1 双变量耦合演算模型结构

河道一维水流，假设没有旁侧入流，由过水断面面积(A)和流量(Q)表达的连续方程:

$\begin{split} \frac{\text{∂}A}{\text{∂}t}+\frac{\text{∂}Q}{\text{∂}x}=0 \end{split}$

(1)

为了使方程闭合，还需要一个描述过水断面面积和流量两者关系的方程，水动力学中以力的平衡为基础构建了两者关系.但由于其中摩阻力与实际差异很大，常导致使用效果不好^[23-26].水文学中提出蓄泄关系，构成马斯京根河道水流演算方法^[1], 但方法只适用于有稳定水位流量关系线的河段^[19-22].本文通过分别建立河段蓄水量与断面面积、河段蓄水量与流量的关系，获得断面面积与流量的关系，构造完备的双变量耦合演算模型(以下简称双变量模型).

对于以体积单位表达的河段蓄量(W)，可表示为河段长(L)和河段平均过水断面面积( $\bar A$ )的乘积:

$\begin{split} W=L·\bar A \end{split}$

(2)

公式（2）中L又可以表达为河段内水流运动平均速度( $\bar u$ )和其传播时间(K)的乘积:

$\begin{split} L=K·\bar u \end{split}$

(3)

由公式(2) 和(3)，蓄量又可以表达为水流传播时间与平均流量( $\bar Q$ )的乘积:

$\begin{split} W=K·\bar Q \end{split}$

(4)

比较公式（2）和（4），得:

$\begin{split} L\bar A=K·\bar Q \end{split}$

(5)

公式(5) 建立了 $\bar Q$ 和 $\bar A$ 间的关系.显然，只要有 $\bar Q$ 和 $\bar A$ 与其上下断面相应要素的关系，就可以与公式（1）构成闭合模型.这里先假设两者都存在加权平均关系:

$\begin{split} \bar A=α·A_j+\left(1－α\right)A_{j+1} \end{split}$

(6)

$\begin{split} \bar Q=χ·Q_j+\left(1－χ\right)Q_{j+1} \end{split}$

(7)

式中, $α$ 和 $χ$ 为权系数，其值一般在0~1之间变化，下标 $j$ 和 $j+1$ 分别表示上断面和下断面位置.

则公式(5)、式(6) 和式(7) 组成的断面面积与流量关系结构：

$\begin{split} L·α·A_j+L\left(1－α\right)A_{j+1}=K·χ·Q_j+K\left(1－χ\right)Q_{j+1} \end{split}$

(8)

天然河道的断面过水断面形状比较复杂且沿程变化^[27]，可以假设平均断面面积等于上、下断面面积的加权平均，公式（7）代入（4）就是马斯京根法的蓄泄关系，其结构的合理性已有大量证明^[1-18].所以公式（5）、（6）和（7）组成的面积与流量关系结构是具有合理性的.

公式（1）采用Pressimann四点隐式差分格式，有：

$\begin{split} \frac{θ}{\text{Δ}x}\text{Δ}Q_{j+1}+\frac{1}{2\text{Δ}t}\text{Δ}A_{j+1}=\frac{θ}{\text{Δ}x}\text{Δ}Q_j－\frac{1}{2\text{Δ}t}\text{Δ}A_j+\frac{1}{\text{Δ}x}\left(Q^i_j－Q^i_{j+1}\right) \end{split}$

(9)

公式（8）考虑前后时间相减得：

$\begin{split} －K\left(1－χ\right)\text{Δ}Q_{j+1}+L\left(1－α\right)\text{Δ}A_{j+1}=K·χ·\text{Δ}Q_j－L·α·\text{Δ}A_j \end{split}$

(10)

式中, $θ$ 为Pressimann差分时间加权系数， $\text{Δ}f=f^{i+1}－f^i$ ，上标 $i$ 和 $i+1$ 分别表示时段前和后的时间.公式（10）和（9）中等号右边为已知，等号左边为未知，两个方程，两个时段差变量，构成了双变量模型.两个演算方程系数为：

$\begin{split} \left| \begin{array}{c} \frac{θ}{\text{Δ}x}&\frac{1}{2\text{Δ}t}\\－K\left(1－χ\right)&L\left(1－α\right) \end{array} \right|=\frac{θ}{\text{Δ}x}L\left(1－α\right)+\frac{K\left(1－χ\right)}{2\text{Δ}t}>0 \end{split}$

(11)

公式（11）对于任意在(0, 1) 范围内变化的权系数 $α$ 和 $χ$ 都成立，证明公式（9）和（10）相对于演算变量不线性相关，说明新构建的关系式（10）有效.

2 应用河段选择

为了检验模型结构的合理性、应用于实际河段模拟的有效性，选择河段考虑旁侧入流比例小和河段距离短两个条件，主要考虑如下3方面因素：

（1）两个河段距离短些.因为连续方程差分以河段长为步长，这步长值越大，差分误差就大；

（2）区间面积不大，旁侧入流比例小.因为连续方程式（1）忽略了旁侧入流项，旁侧入流比例大，将导致这忽略带来的误差就大；

（3）不同水力和河道断面特征的河段.使得模型对不同河道和水力特征具有广泛的代表性.

根据这3条原则，本次研究选择了全国涉及长江、黄河、珠江、辽河、松花江五大水系，还选择了塔里木河等内陆河流和入海河流，具有各种水力、各种河道断面特征的代表性，其河段特征与资料年份见表 1.选择的洪水是每年汛期的洪水过程，为了考虑检验的充分性，选择洪水历时尽可能长，一般短则10 d，长的1个月，最长的连续3个月；洪水选择也考虑大、中、小各种洪水的代表性；年份选择主要考虑前后相连；计算时段间隔采用半小时.

3 模型与检验结果分析

为了充分检验模型的效果，同时把马斯京根法模型参数和双变量模型参数进行率定，马斯京根模型把3个汇流参数进行率定（考虑区间来水，3系数之和不为1），这样参数个数与双变量模型相同.马斯京根汇流演算模型为：

$\begin{split} Q·C_{t+1}=C_0·I_{t+1}+C_1·I_t+C_2·Q·C_t \end{split}$

(12)

双变量模型，为了与马斯京根法比较，将式（9）代入（10），消去 $\text{Δ}A_{j+1}$ , 得:

$\begin{split} &ΔQ_{j+1}=B_{1}·ΔQ_{j}+B_{2}·ΔA_{j}+B_{3}(Q^{i}_{j}－Q^{i}_{j+1})\\ &B_{1}=(2θ·L(1－α)Δt－K·χ·Δx)/(2θ·L(1－α)Δt+K(1－χ)·Δx)\\ &B_{2}=L(2α－1)Δx/(2θ·L(1－α)Δt+K(1－χ)·Δx)\\ &B_{3}=L(1－α)2Δt/(2θ·L(1－α)Δt+K(1－χ)·Δx) \end{split}$

(13)

与马斯京根法类似，以I表达上断面流量，以A表达上断面面积，则可得双变量模型的流量演算式为：

$\begin{split} &Q·C_{t+1}=Q·C_{t}+B_{1}(I_{t+1}-I_{t})+B_{2}(A_{t+1}－A_{t})+B_{3}(I_{t}－Q·C_{t}) \end{split}$

(14)

也可类似地将式（9）代入（10），消去 $\text{Δ}Q_{j+1}$ 得 $\text{Δ}A_{j+1}$ 的下断面面积演算式，获得了与计算流量类似的结果.本文限于篇幅不讨论断面面积的计算结果.

模型检验，除辽河外的河段选择最后一年的洪水资料为模型检验，其余年份为参数率定.检验期洪水演算模型模拟的效果统计见表 2. DC(Ma)和DC(New)分别为马斯京根和双变量模型率定期确定性系数，VDC(Ma)和VDC(New)分别表示马斯京根和双变量模型检验期确定性系数计算结果.

从率定期和检验期结果看，16个河段中有15个河段双变量模型效果好，河段平均确定性系数也是双变量模型高些；从各河段稳定性看，马斯京根法的确定性系数变幅很大，双变量模型更稳定.

为进一步分析比较模型模拟实际洪水的情况，将资料系列最长的辽河铁岭-马虎山河段检验结果进行详细分析.铁岭-马虎山河段采用的洪水资料为1955-1994年间选择了18年汛期的所有大洪水和几次代表性的中小洪水，其中前期10年的资料用于模型参数率定，后期8年11场洪水用于模型检验. 表 3列出了辽河铁岭-马虎山河段检验期11场洪水的检验结果，QCp(Ma)和QCp(New)分别表示马斯京根模型和双变量模型计算的洪峰；DC(Ma)和DC(New)分别表示马斯京根模型和双变量模型计算的确定性系数.

从表 3检验期各场洪水确定性系数看，双变量模型的平均确定性系数为0.753，略高于马斯京根的0.734，各场洪水变幅也比马斯京根法更小，与16个河段结论一致.

7号是马斯京根模型模拟最差的洪水，8号是双变量模型模拟最差的洪水.对比马斯京根法计算流量与实测流量过程，比较8号洪水（图 1b，1d），8号洪水降雨时空分布不均匀，洪水的区间降雨主要出现在洪峰附近，造成洪峰附近流量叠加使得实测流量大于双变量模型计算流量，双变量模型受区间来水影响；比较7号洪水(图 1a，1c)，马斯京根模型受误差累计影响. “马斯京根模型受误差累计影响”主要体现在储蓄结构中，因为当储蓄关系与实际有出入时，其计算蓄量偏差需要连续几个时段的累积才能弥补其影响.马斯京根法计算的7号洪水，初始计算蓄量偏小，计算出流也偏小，使得涨水段计算流量连续偏小，其误差通过涨水段的累计直到峰后才近似抵消.所以双变量模型改进需要更细地考虑区间来水模拟结构，而马斯京根模型需要改进模型结构，即更细地考虑河段断面面积变化对蓄量的影响结构.

4 结语

本研究根据河段蓄水量等于平均过水断面面积与河段长乘积和蓄水量等于河段平均流量与传播时间的乘积两个物理计算公式，提出了过水断面面积与流量的关系：

$\begin{split} L·α·A_j+L\left(1－α\right)A_{j+1}=K·χ·Q_j+K\left(1－χ\right)Q_{j+1} \end{split}$

(15)

与水流连续方程构成了双变量耦合演算模型：

$\begin{split} \left\{\begin{array}{l}\frac{θ}{\text{Δ}x}\text{Δ}Q_{j+1}+\frac{1}{2\text{Δ}t}\text{Δ}A_{j+1}=\frac{θ}{\text{Δ}x}\text{Δ}Q_j－\frac{1}{2\text{Δ}t}\text{Δ}A_j+\frac{1}{\text{Δ}x}\left(Q^i_j－Q^i_{j+1}\right)\\ －K\left(1－χ\right)\text{Δ}Q_{j+1}+L\left(1－α\right)\text{Δ}A_{j+1}=K·χ·\text{Δ}Q_j－L·α·\text{Δ}A_j\end{array}\right. \end{split}$

(16)

通过四大水系(包括内陆河流和入海河流)的16个河段汛期洪水资料的模拟检验，将双变量耦合通用演算模型与传统的马斯京根法进行了效果比较，结果表明双变量耦合演算模型的模拟精度高于马斯京根法，而且模拟效果更加稳定.

本研究提出的模型在结构合理性方面与现有的马斯京根流量演算模型相近，效果比马斯京根法稍好一些，适用于区间来水影响较小的河段.主要优点是具有通用性，模型结构简单，使用资料要求低，实用性强.与马斯京根模型类似，双变量耦合模型可与新安江模型等结合进行洪水预报，成果具有进一步研究和推广应用价值.