水力学实践中经常遇到水流反向演算问题。例如，防洪调度中由河道防洪标准确定水库出库过程^[1]、由下游水文站点实测流量推算区间入流^[2-3]，农田灌溉中由需水量控制渠首引水过程^[4]，以及突发水污染应急处置中由下游考核浓度确定上游的水力调控方案^[5-6]等，都离不开水流反向演算。水流反向演算一般用下游站点所观测的水位流量信息反算上游的水流过程^[7]，是典型的工程水力学反问题。由于水流沿程传播存在坦化的特性，直接对水动力方程进行反向求解存在误差放大、结果不稳定等问题^{[1, 8-11]}。

由于水流反向演算的重要性，国内外关于该问题的研究一直持续不断。关志成等^[8]和钟平安等^[9]分析了采用Muskingum模型进行反向洪水演算存在不稳定的问题，并分别对模型进行了改进。唐文涛等^[12]结合Muskingum演算和区间半分法由坝址处流量推算水库入库过程，将反向演算问题转化为顺时序正演过程。Das^[13]区分正、反向Muskingum演算参数，通过不断校正Muskingum反向参数实现水流过程的反向演算。李兰^[14]采用脉冲谱方法和离散反向演算推导了扩散波的时空反向演算公式，动态反算洪水参数。Eli等^[15]探讨了不同时间、空间反向组合方式下隐格式数值求解Saint Venant方程的稳定性问题，并推荐了空间逆向、时间正向的反向演算方式。Abdulwahid等^[16]论证了水流反向演算中隐格式求解Saint Venant方程相对于显格式具有更好的稳定性。D'Oria等^[17]将Bayes统计法用于河渠水流反向演算，并通过人工案例验证了模型的精度和实用性。此外，一些启发式算法诸如粒子群^[18]、遗传算法^[19]、反向传播(back propagation，BP)神经网络等^[20]也相继被引入河道水流反向演算研究中。

截至目前，虽然关于水流反向演算的研究成果颇多，但大都基于简化的水流方程^[21]，如运动波和扩散波方程，或者仍存在数值计算稳定性不足的问题。为此，本文对水流正向演算模型所在的x-t平面实施逆时针旋转变换，以空间逆向、时间正向的方式对Saint Venant方程进行隐式离散与求解。本研究基于旋转的x-t平面将正向模型的初始、边界条件分别变换为反向模型的边界、初始条件，解决了以往研究中逆向数值迭代不稳定、不收敛的问题。通过理想河道验证后，将模型应用于西江天然河道洪水反向演算实例，进一步验证了模型的实用性。

1 水流计算模型 1.1 正向演算模型

水流在重力、水压力和沿程阻力作用下可由动量守恒推导其运动方程，再联合水流连续方程后得到Saint Venant方程组。具体方程如下：

$ {\frac{{\partial A}}{{\partial t}} + \frac{{\partial Q}}{{\partial x}} = q,} $

(1)

$ {\frac{{\partial Q}}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {\frac{{{Q^2}}}{A}} \right) + gA\frac{{\partial Z}}{{\partial x}} + g\frac{{n_{\rm{s}}^2Q|Q/A|}}{{R_{\rm{w}}^{\frac{4}{3}}}} = 0.} $

(2)

式中：t为时间坐标，s；x为空间坐标，m；A为过流断面面积，m²；B为水面宽度，m；Q为断面流量，m³/s；q为单位河段区间入流量，m³/(s·m)；Z为断面水位，m；n_s为河床糙率系数；R_w为过流断面水力半径，m。

考虑到隐格式稳定性好、对离散时空步长限制要求相对较低，本文采用Preissmann四点偏心隐格式差分方法^[22]离散Saint Venant方程，见图 1。图 1中：j为断面编号；Δx为离散空间步长；Δt为离散时间步长；θ为偏心系数，取值一般为0.5~1.0。

固定计算时间n (n代表离散t后的时间序列编号)，对任意河道断面j，可得到方程(3)和(4)：

$ \begin{array}{*{20}{c}} {a_{2j - 1}^1\Delta {Q_j} + a_{2j - 1}^2\Delta {Z_j} + a_{2j - 1}^3\Delta {Q_{j + 1}} + }\\ {a_{2j - 1}^4\Delta {Z_{j + 1}} = {b_{2j - 1}},} \end{array} $

(3)

$ a_{2j}^1\Delta {Q_j} + a_{2j}^2\Delta {Z_j} + a_{2j}^3\Delta {Q_{j + 1}} + a_{2j}^4\Delta {Z_{j + 1}} = {b_{2j}}. $

(4)

式中：ΔQ_j=Q_jⁿ⁺¹-Q_j^*, ΔZ_j=Z_jⁿ⁺¹-Z_j^*, Q_j^*、Z_j^*分别为Q_jⁿ、Z_jⁿ的循环迭代更新值(下文上标中带*号的皆为循环迭代更新值)；系数a_2j-1ⁱ、a_2jⁱ (i=1，2，3，4)与b_2j-1、b_2j由上一时刻变量Q_jⁿ、Z_jⁿ与迭代值Q_j^*、Z_j^*确定，各系数计算表达式如下：

$ {a_{2j - 1}^1 = - a_{2j - 1}^3 = - 1,} $

$ {a_{2j - 1}^2 = a_{2j - 1}^4 = \frac{{B_{j + 1/2}^*\Delta x}}{{2\Delta t{\kern 1pt} \theta }},} $

$ \begin{array}{l} {b_{2j - 1}} = \frac{{\theta (Q_j^{n + 1} + qQ_{j + 1}^{n + 1}) + (1 - \theta )(q_j^n + q_{j + 1}^n)}}{{2\theta }}\Delta x - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{1 - \theta }}{\theta }(Q_{j + 1}^n - Q_j^n) + a_{2j - 1}^2(Z_j^n + Z_{j + 1}^n) - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (a_{2j - 1}^1Q_j^* + a_{2j - 1}^2Z_j^* + a_{2j - 1}^3Q_{j + 1}^* + a_{2j - 1}^4Z_{j + 1}^*), \end{array} $

$ {a_{2j}^1 = \frac{{\Delta x}}{{2\Delta t\theta }} - \left( {\frac{Q}{A}} \right)_j^* + \left( {g\frac{{n_{\rm{s}}^2|Q/A|}}{{2{\kern 1pt} {\kern 1pt} \theta R_{\rm{w}}^{\frac{4}{3}}}}} \right)_j^*\Delta x,} $

$ {a_{2j}^2 = - a_{2j}^4 = - g(A_j^* + A_{j + 1}^*)/2,} $

$ a_{2j}^3 = \frac{{\Delta x}}{{2\Delta t{\kern 1pt} \theta }} + \left( {\frac{Q}{A}} \right)_{j + 1}^* + \left( {g\frac{{n_{\rm{s}}^2|Q/A|}}{{2{\kern 1pt} \theta R_{\rm{w}}^{\frac{4}{3}}}}} \right)_{j + 1}^*\Delta x, $

$ \begin{array}{l} {b_{2j}} = \frac{{\Delta x}}{{2\Delta t\theta }}(Q_{j + 1}^n + Q_j^n) - \frac{{1 - \theta }}{\theta }\left[ {\left( {\frac{{{Q^2}}}{A}} \right)_{j + 1}^n - \left( {\frac{{{Q^2}}}{A}} \right)_j^n} \right] + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{1 - \theta }}{\theta }a_{2j - 1}^2(Z_{j + 1}^n - Z_j^n) - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (a_{2j}^1Q_j^* + a_{2j}^2Z_j^* + a_{2j}^3Q_{j + 1}^* + a_{2j}^4Z_{j + 1}^*). \end{array} $

由于Saint Venant方程的双曲特性，对于河道首末断面，须补充边界条件：

$ {\rm{首断面}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} a_0^3\Delta {Q_1} + a_0^4\Delta {Z_1} = {b_0}, $

(5)

$ {\rm{末断面}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} a_{2L - 1}^1\Delta {Q_L} + a_{2L - 1}^2\Delta {Z_L} = {b_{2L - 1}}{\rm{.}} $

(6)

式中L代表末断面编号。式中系数由外边界条件确定，一般是上游给定入流过程、下游给定流量水位关系曲线^[23]：

$ {{Q_1} = {Q_{{\rm{up}}}}(t),} $

(7)

$ {{Q_L} = \alpha {{({Z_L} + \beta )}^\gamma }.} $

(8)

式中：α、β、γ为经验拟合系数，Q_up(t)为上游实测入流系列。

离散式(7)和(8)，有：

$ {a_0^3 = 1,a_0^4 = 0,{b_0} = {Q_{{\rm{ up }}}}(t) - Q_1^*,} $

$ {a_{2L - 1}^1 = 1,a_{2L - 1}^2 = - \alpha \gamma {{(Z_L^* + \beta )}^{\gamma - 1}},} $

$ {{b_{2L - 1}} = \alpha {{(Z_L^* + \beta )}^\gamma } - Q_L^*.} $

在给定上、下游边界条件以及初始条件后，可采用双扫描法求解线性方程组(3)和(4)，实现水流的动态计算^[24]。在求解得到未知量ΔQ_j、ΔZ_j后，判断其是否满足误差控制条件|ΔQ_j| < ε₁和|ΔZ_j| < ε₂(ε₁、ε₂分别为流量水位误差控制限)，若满足，则迭代值Q_j^*、Z_j^*即为当前计算时刻所求Q_jⁿ⁺¹、Z_jⁿ⁺¹，否则令Q_j^*=Q_j^*+ΔL_j和Z_j^*=Z_j^*+ΔZ_j，重新迭代求解方程，直到得到的结果满足误差控制条件。

1.2 反向演算模型

根据时间、空间反向方式的不同组合，水流反向演算存在3种不同形式，如图 2所示(图中T代表正向模型总计算时间)。根据多次数值实验，图 2中3种反算方式都存在不稳定或者不收敛的问题^[15-16]。其中：图 2d的反算方式仍须给定上边界条件，在水流反向演算实践中并不实用；图 2b的双反向方式在稳定性方面明显劣于图 2c的单反向方式^[15]。

特别说明，尽管图 2c采用了时间正向的方式，由历史时刻水流状态确定初始条件，但水动力方程具有双曲特性，水流初始条件影响时域很有限(图 3)，因此可通过恒定流预热的方式给定反向演算所需的初始条件。图 3中，先根据水流传播时间(水流扰动从河道上边界传播至下边界的时间)确定水流预热时长，然后分别以下边界观测的初始和末时刻流量作为恒定流量进行预热，将预热结束的水流状态作为外延后(图 3a中从OT外延至O′T′)的初始条件。

水流反向演算不稳定的主要原因在于^[8-9]：物理上，水流反向演进时会出现陡化；数学上，逆向数值迭代时会出现误差逐渐放大。因此，在图 2c的反算方式下，将x-t平面逆时针旋转90°，交换时、空坐标方向，以空间逆向、时间正向的方式对水动力方程进行离散，x-t平面旋转变换可实现反向模型在时间维度上的正向迭代，进而可解决逆向数值迭代不稳定、不收敛的问题，改善水流反向演算的稳定性，如图 4所示。

图 4中，x-t平面逆时针旋转90°后，正向演算模型的初始、边界条件(图 4中左图)分别变换成反向演算模型的边界、初始条件(图 4中右图)。同样采用图 1中差分格式离散Saint Venant方程以构建反向演算模型。

固定空间断面j，对任意时间n，可得到线性方程(9)和(10)：

$ \begin{array}{l} a_{2n - 1}^1\Delta {Q^n} + a_{2n - 1}^2\Delta {Z^n} + a_{2n - 1}^3\Delta {Q^{n + 1}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} a_{2n - 1}^4\Delta {Z^{n + 1}} = {b_{2n - 1}}, \end{array} $

(9)

$ \begin{array}{l} a_{2n}^1\Delta {Q^n} + a_{2n}^2\Delta {Z^n} + a_{2n}^3\Delta {Q^{n + 1}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} a_{2n}^4\Delta {Z^{n + 1}} = {b_{2n}}. \end{array} $

(10)

式中：ΔQⁿ=Q_jⁿ-Q_*ⁿ, ΔZⁿ=Z_jⁿ-Z_*ⁿ, Q_*ⁿ、Z_*ⁿ分别为Q_j+1ⁿ、Z_j+1ⁿ的循环迭代更新值；系数a_2n-1ⁱ、a_2nⁱ (i=1, 2, 3, 4)与b_2n-1、b_2n由前一断面变量Q_j+1ⁿ、Z_j+1ⁿ和迭代值Q_*ⁿ、Z_*ⁿ确定。各系数计算表达式如下：

$ {a_{2n - 1}^1 = \frac{{1 - \theta }}{\theta },a_{2n - 1}^3 = 1,} $

$ {a_{2n - 1}^2 = - a_{2n - 1}^4 = \frac{{B_*^{n + \theta }\Delta x}}{{2\Delta t\theta }},} $

$ {B_*^{n + \theta } = \theta B_*^{n + 1} + (1 - \theta )B_*^n,} $

$ \begin{array}{l} {b_{2n - 1}} = - \frac{{\theta (q_j^{n + 1} + q_{j + 1}^{n + 1}) + (1 - \theta )(q_j^n + q_{j + 1}^n)}}{{2\theta }}\Delta x + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} Q_{j + 1}^{n + 1} + \frac{{1 - \theta }}{\theta }Q_{j + 1}^n + a_{2n - 1}^2(Z_{j + 1}^{n + 1} - Z_{j + 1}^n) - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (a_{2n - 1}^1Q_*^n + a_{2n - 1}^2Z_*^n + a_{2n - 1}^3Q_*^{n + 1} + a_{2n - 1}^4Z_*^{n + 1}), \end{array} $

$ \begin{array}{l} a_{2n}^1 = \frac{{\Delta x}}{{2\Delta t\theta }} + \frac{{1 - \theta }}{\theta }\left( {\frac{Q}{A}} \right)_*^n - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {g\frac{{1 - \theta }}{{2\theta }}\frac{{n_{\rm{s}}^2\left| {Q/A} \right|}}{{R_{\rm{w}}^{\frac{4}{3}}}}} \right)_*^n\Delta x, \end{array} $

$ {a_{2n}^2 = gA_*^{n + \theta }\frac{{1 - \theta }}{\theta },a_{2n}^4 = gA_*^{n + \theta },} $

$ {A_*^{n + \theta } = \theta A_*^{n + 1} + (1 - \theta )A_*^n,} $

$ \begin{array}{*{20}{l}} {a_{2n}^3 = - \frac{{\Delta x}}{{2\Delta t\theta }} + \left( {\frac{Q}{A}} \right)_*^{n + 1} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {g\frac{{n_{\rm{s}}^2|Q/A|}}{{2R_{\rm{w}}^{\frac{4}{3}}}}} \right)_*^{n + 1}\Delta x,} \end{array} $

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{c}} {{b_{2n}} = \frac{{\Delta x}}{{2\Delta t\theta }}(Q_{j + 1}^{n + 1} - Q_{j + 1}^n) + }\\ {\left( {\frac{Q}{A}} \right)_*^{n + 1}Q_{j + 1}^{n + 1} + \frac{{1 - \theta }}{\theta }\left( {\frac{Q}{A}} \right)_*^nQ_{j + 1}^n + } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} gA_*^{n + \theta }\left( {Z_{j + 1}^{n + 1} + \frac{{1 - \theta }}{\theta }Z_{j + 1}^n} \right) + gn_{\rm{s}}^2\frac{{\Delta x}}{2} \cdot \\ \left[ {{{\left( {\frac{{|Q/A|}}{{R_{{\rm{ w }}}^{\frac{4}{3}}}}} \right)}^{n + 1}}Q_{j + 1}^{n + 1} + \frac{{1 - \theta }}{\theta }\left( {\frac{{|Q/A|}}{{R_{{\rm{ w }}}^{\frac{4}{3}}}}} \right)_*^nQ_{j + 1}^n} \right] - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (a_{2n}^1Q_*^j + a_{2n}^2Z_*^n + a_{2n}^3Q_*^{n + 1} + a_{2n}^4Z_*^{n + 1}). \end{array} $

在给定上下游边界条件和初始条件后，同样采用双扫描法迭代求解方程组(9)和(10)。

为便于反向演算模型的理解和实现，进一步比较水流正、反向演算模型，见表 1。

1.3 观测误差模拟

工程实践中所获取的下游实测数据都带有一定的观测误差。在反向演算模型数值实验时，应考虑观测误差的影响。根据前人研究成果^[25]，观测误差一般可作为随机变量模拟，且服从均值为0的Gauss分布，本文采用方程(11)模拟流量观测误差:

$ Q(t) = {Q_0}(t)(1 + {\varepsilon _Q}(t)),{\varepsilon _Q}(t)\backsim N(0,{h^2}). $

(11)

式中h为误差扰动比例系数^[26]，反映观测误差水平。

1.4 模型误差评价指标

本研究采用相关系数(R²)、相对均方根误差$\left( {\frac{{{\rm{RMSE}}}}{{\bar O}}} \right)$和Nash-Sutcliffe效率系数(NSE) 3项误差指标评价模型模拟效果，计算公式如下：

$ {R^2} = {\left( {\frac{{\sum\limits_{t = 1}^T {({O_t} - \bar O)} ({M_t} - \bar M)}}{{\sqrt {\sum\limits_{t = 1}^T {{{({O_t} - \bar O)}^2}} } \sqrt {\sum\limits_{t = 1}^T {{{({M_t} - \bar M)}^2}} } }}} \right)^2}, $

(12)

$ \frac{{{\rm{RMSE}}}}{{\bar O}} = \frac{{\sqrt {\frac{1}{T}\sum\limits_{t = 1}^T {{{({O_t} - {M_t})}^2}} } }}{{\frac{1}{T}\sum\limits_{t = 1}^T {{O_t}} }}, $

(13)

$ {\rm{NSE}} = 1 - \frac{{\sum\limits_{t = 1}^T {{{({O_t} - {M_t})}^2}} }}{{\sum\limits_{t = 1}^T {{{({O_t} - \bar O)}^2}} }}. $

(14)

式中: M_t为模拟值，O_t为实测值，M和O分别为模拟、观测系列均值。

在水流反向演算结果评价时存在两种情况：1)上游断面有部分实测数据，此时可直接利用实测数据对反向演算结果进行评价；2)完全缺少上游实测数据，此时可将模型反算的上游入流结果作为上边界，然后再进行水流正向模拟，将模拟得到的观测断面水流过程与实测过程进行对比分析。

2 理想河道案例分析

首先将建立的水流反向演算模型应用于理想河道案例以检验模型的有效性。假设有一棱柱形梯形断面河道，底宽20 m，边坡2.5，底坡0.000 059，糙率系数0.017，河道全长15.1 km。本研究考虑水流缓慢变化(情景1)和快速变化(情景2)两种不同情景，见图 5(不考虑观测误差)。

2.1 结果分析

不考虑观测误差的影响，模型反算的入流过程与实际入流过程对比见图 6。

从图 6中可以看出，在不考虑观测误差影响下，本文构建的水流反向演算模型能很好地还原出实际入流过程。在流量快速变化情景下，反算的入流过程在峰值处存在一定误差，该结果与数值离散误差有关。此外，两种情景下模型都表现出很好的稳定性，反算结果仅在水流退水后期出现轻微的波动。

为进一步评估模型反向演算误差，将反算的入流过程作为上边界条件，由正向模型计算出流过程，并与实际过程进行对比，见图 7。

在缺少入流观测数据时，图 7中拟合结果可用于验证入流反算结果，以判断反向演算模型的合理性。结合图 6和7，在不考虑观测误差影响下，本文构建的水流反向演算模型对于不同水流情景都具有很好的稳定性，且计算精度高，模型有效可行。

2.2 观测误差影响分析

考虑观测误差的影响，分析误差扰动比例系数h=0.1%和h=1%两种观测误差情景下模型反算结果，见图 8。

对比图 6和8可以看出，在观测误差较小时(h=0.1%)，观测误差对反算结果的影响可忽略，而在观测误差相对较大时(h=1%)，观测误差的存在使得反算结果出现了局部波动。在模型实际应用时可考虑采用EnKF等数据同化方法对反算结果进行降噪处理^{[2, 27]}。

为进一步评估反向演算模型的精度，对不同情景下的反算入流过程进行误差分析，结果见表 2。

表 2中，不同情景下的误差指标很接近，R²和NSE皆大于0.95。表 2结果表明，反向演算模型受观测误差影响很小，3个误差指标都表明了模型的精度高，具有实用价值。

3 模型应用

经理想河道案例验证后，将反算模型应用于西江浔江河段的实例中，以检验模型的实用性。浔江河段为西江主干河道，位于中国广西境内。该段河道沿程分布有平南县、藤县、梧州市等，是西江流域洪水高发河段。本研究主要分析大湟江口至长洲库区入口(距离大湟江口130.74 km)河段，在距离大湟江口73.89 km处有支流蒙江汇入，具体位置如图 9所示。

将构建的水流反向演算模型用于该天然河道洪水反向演算，分析大湟江口具有不同特征的3场次洪水过程(考虑支流蒙江入流)。模型反算的入流过程与实测入流过程对比结果见图 10和表 3。

从图 10和表 3中可看出，模型反算的3场次洪水入流过程皆与实测入流过程高度一致，R²和NSE都高达0.99。结果表明，本文构建的反向演算模型可应用于天然河道洪水反向演算，模型精度高、稳定性好。实例应用再次验证了模型有效可行。

4 结论

本文对水流正向演算模型所在的x-t平面进行逆时针旋转变换，以空间逆向、时间正向的方式将正向模型的初始、边界条件分别变换成反向模型的边界、初始条件，构建了河道水流反向演算数值模型。将模型先后应用于理想河道和西江天然河道案例，结果表明：本文提出的水流反向演算模型稳定性好、精度高，反算的入流结果中R²和NSE皆在0.95以上，模型有效可行。未来可考虑将水流反向演算原理推广应用于水质反向演算。