由于流域降雨径流过程受到降雨和下垫面条件等多重因素影响，实时洪水预报精度难以满足要求.预报模型初始状态值误差^[1]、模型输入误差^[2]、模型结构的不充分性^[3]以及参数非优化^[4]都会影响预报结果的精度.所以经常使用误差修正技术来减少模型在预报时产生的误差^[5-6].

降雨径流过程涉及蒸发、下渗、产流和汇流等多个环节，影响因素复杂.简单的误差修正方法是通过流量误差直接对模型计算结果(流量过程)进行修正^[7].组合卡尔曼滤波(EnKF)或者粒子滤波(PF)方法经常通过数据同化技术对模型初始值或者模型参数进行修正，然而研究结果表明滤波方法在预报修正中的效果过于自信，对于模型预报的不确定性估计不足^[8-9].有学者尝试使用滤波方法对降雨资料误差进行估计和修正，但结果显示仅有在土壤含水量很低且流量资料非常精确的情况下才有明显效果^[10].

影响实时洪水预报系统的误差因素较多，其影响机理非常复杂^[11].产流量是衡量洪水预报水量平衡的重要指标，笔者曾通过直接修正产流量来提高洪水预报精度^[12-14].但引起产流量计算误差的因素较多，通过对这些因素的归纳分析可以发现，面平均雨量是洪水预报模型的最重要输入项，面平均雨量误差直接影响初始值和计算产流量，从而影响模型预报精度^[15]，对于集总式水文模型尤为明显.降雨过程是受到一系列因素(气候、地形、植被以及人类活动)影响的随机过程，此类误差的特性以及误差的大小难以预测和估计^[16].由于遥感降雨资料本身具有误差，卫星和雷达降雨资料的使用同样无法解决降雨资料误差的问题^[17].集总式水文模型虽然在体现和反映真实暴雨过程剧烈的时间和空间变异性方面不如分布式水文模型，但由于在实际洪水预报作业过程中流域水文气象资料(蒸发、降雨)和流域下垫面条件(土壤植被类型、土壤含水量)等无法满足分布式水文模型的要求，所以集总式或半分布式水文模型在实时洪水预报中的应用更广泛^[18]，因此本研究中采用的水文模型为集总式水文模型.

因此，本文提出一种相对简单的降雨系统响应曲线修正方法，通过修正面平均雨量来提高模型预报精度.将此方法与新安江模型结合，通过理想案例验证和实际流域应用，并进一步检验此方法在不同雨量站密度情况下的应用效果.

1 修正方法 1.1 降雨系统响应曲线

将水文模型作为一般系统考虑，图 1中P(t)为面平均降雨系列，E(t)为蒸散发，Q(t)为计算的流量过程，θ为模型参数，X(t)是模型状态变量，X₀为状态初始值. 图 1中的水文模型系统可以表达为由输入、输出、参数和状态变量组成的系统方程：

$ Q\left( t \right) = f\left[ {P\left( t \right),E\left( t \right),X\left( t \right),\theta } \right] $

(1)

式中, Q=[q₁, q₂, q₃, …, q_n]^T为流域出口断面流量过程，P=[p₁, p₂, p₃, …, p_m]^T为面平均雨量向量，E=[e₁, e₂, e₃, …, e_n]^T为蒸发，X=[x₁, x₂, x₃, …, x_n]^T为状态变量向量，θ为参数向量，n为时段流量系列个数，m为时段降雨系列个数.

概念性模型的绝大多数参数都具有明确的物理意义，是流域水文特征的反应，模型参数一旦率定之后一般在短期之内不会改变^[19].模型状态变量取决于模型的输入输出项^[20]，因此，在本研究中认为模型参数不变化，面平均雨量的变化是引起流域出口断面流量变化的主要因素.上述系统方程可以简化为：

$ Q\left( P \right) = f\left( P \right) $

(2)

对于流域水文模型，把降雨输入作为自变量求全微分，可得到式(2)的微分关系：

$ {\rm{d}}Q = \frac{{\partial Q\left( P \right)}}{{\partial P}}\left| {_{P = {P_0}}} \right.{\rm{d}}P $

(3)

式中, dP为面平均雨量的变化量，dQ为由面平均雨量变化量引起的流量变化量，P₀为修正之前的面平均雨量系列，${\left. {\frac{{\partial Q}}{{\partial P}}} \right|_{P = {P_{\rm{0}}}}} $为系统在面平均系列为P₀时的系统响应.则观测流量可以表达为：

$ Q\left( P \right) \approx Q\left( {{P_0}} \right) + \frac{{\partial Q\left( P \right)}}{{\partial {p_1}}}\left| {_{P = {P_0}}} \right.\Delta {p_1} + \frac{{\partial Q\left( P \right)}}{{\partial {p_2}}}\left| {_{P = {P_0}}} \right.\Delta {p_2} + L + \frac{{\partial Q\left( P \right)}}{{\partial {p_m}}}\left| {_{P = {P_0}}} \right.\Delta {p_m} $

(4)

式中, Q(P₀)为由降雨系列P₀计算得到的流量系列，${\left. {\frac{{\partial Q\left( P \right)}}{{\partial {p_m}}}} \right|_{P = {P_{\rm{0}}}}}\Delta {p_m} $表示第m时段的降雨变化量Δp_m对出口流量过程Q(P)所造成的系统响应.假设使用观测流量系列长度为L，Q(t)=[Q₁, Q₂, Q₃, …, Q_L]^T.代入式(4)可得：

$ \left\{ \begin{array}{l} {Q_1}\left( P \right) \approx {Q_1}\left( {{P_0}} \right) + \frac{{\partial {Q_1}\left( P \right)}}{{\partial {p_1}}}\left| {_{P = {P_0}}} \right.\Delta {p_1} + \frac{{\partial {Q_1}\left( P \right)}}{{\partial {p_2}}}\left| {_{P = {P_0}}} \right.\Delta {p_2} + L + \frac{{\partial {Q_1}\left( P \right)}}{{\partial {p_m}}}\left| {_{P = {P_0}}} \right.\Delta {p_m}\\ {Q_2}\left( P \right) \approx {Q_2}\left( {{P_0}} \right) + \frac{{\partial {Q_2}\left( P \right)}}{{\partial {p_1}}}\left| {_{P = {P_0}}} \right.\Delta {p_1} + \frac{{\partial {Q_2}\left( P \right)}}{{\partial {p_2}}}\left| {_{P = {P_0}}} \right.\Delta {p_2} + L + \frac{{\partial {Q_2}\left( P \right)}}{{\partial {p_m}}}\left| {_{P = {P_0}}} \right.\Delta {p_m}\\ {Q_L}\left( P \right) \approx {Q_L}\left( {{P_0}} \right) + \frac{{\partial {Q_L}\left( P \right)}}{{\partial {p_1}}}\left| {_{P = {P_0}}} \right.\Delta {p_1} + \frac{{\partial {Q_L}\left( P \right)}}{{\partial {p_2}}}\left| {_{P = {P_0}}} \right.\Delta {p_2} + L + \frac{{\partial {Q_L}\left( P \right)}}{{\partial {p_m}}}\left| {_{P = {P_0}}} \right.\Delta {p_m} \end{array} \right. $

(5)

则式(5)的矩阵形式为：

$ Q\left( P \right) = Q\left( {{P_0}} \right) + S \cdot \Delta P + W $

(6)

式中, ΔP=[Δp₁, Δp₂, Δp₃, L, Δp_m]^T即为面平均雨量误差的估计值；W=[w₁, w₂, w₃, …, w_L]^T为流量观测随机误差项，一般为服从零均值分布的为白噪声向量. S为系统响应矩阵，其表达式为：

$ S = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial {Q_1}\left( P \right)}}{{\partial {p_1}}}}& \cdots &{\frac{{\partial {Q_1}\left( P \right)}}{{\partial {p_m}}}}\\ {\frac{{\partial {Q_2}\left( P \right)}}{{\partial {p_1}}}}& \cdots &{\frac{{\partial {Q_2}\left( P \right)}}{{\partial {p_m}}}}\\ {}& \vdots &{}\\ {\frac{{\partial {Q_L}\left( P \right)}}{{\partial {p_1}}}}& \cdots &{\frac{{\partial {Q_L}\left( P \right)}}{{\partial {p_m}}}} \end{array}} \right] $

(7)

式(7)代入式(3)得到由降雨变化量与流量变化量表达的系统响应关系：

$ {\rm{d}}Q = S \cdot {\rm{d}}P $

(8)

式(8)中S矩阵的每一项都可以用式(9)差分近似求解.

$ \frac{{\partial {Q_j}\left( P \right)}}{{\partial {P_i}}} = \frac{{{Q_j}\left( {{p_1}, \cdots ,{p_i} + \Delta {p_i}, \cdots ,{p_m}} \right) - {Q_j}\left( {{p_1}, \cdots ,{p_i}, \cdots ,{p_m}} \right)}}{{\Delta {p_i}}}\left( {i = 1、\cdots 、m,j = 1、\cdots 、L} \right) $

(9)

式中, Q(p₁, …, p_i, …, p_m)为由修正前的面平均雨量系列计算得到的流量系列，Q(p₁, …, p_i+Δp_i, …, p_m)为由面平均雨量系列增加Δp_i之后计算得到的流量系列.

式(8)中当i不变，而j从1~L变化时，该L项差分值就是S矩阵中的一列.而此列正是雨量p_i的单位变化量所对应的系统响应系列，本研究中称之为雨量p_i所对应的系统响应曲线.在实际求解系统响应曲线的过程中，为了计算简单，令式(9)中Δp_i=1，则降雨系统响应曲线就是降雨量的单位变化量引起的系统响应.

计算系统响应曲线的详细步骤如下：(1)面平均雨量系列P中的时段雨量p_i (i=1~m)，在其余时段降雨量p_k(k≠i)均不变的基础上增加1个单位的降雨量Δp_i，得到新的降雨系列用Pⁱ表示. (2)用新的面平均雨量系列Pⁱ通过模型计算后得到流量过程Q_i(t). (3)用Q_i(t)减去用修正前的面平均雨量系列P₀计算得到的流量过程Q(t)所得到的系列值，即为降雨量p_i的系统响应曲线，表示为Sⁱ(t). S矩阵中的每一列均用同样方法求得，示意图如图 2所示.

1.2 降雨误差估计

对于一场洪水，假如有m个时段雨量，其相应的雨量误差依时间顺序为e_P1, e_P2, …, e_Pm，假如有L个时段的流量系列，流量的误差序列为e_Q1, e_Q2, …, e_QL，则降雨误差与流量误差之间的关系可以表达为式(10)：

$ \left\{ \begin{array}{l} {e_{{Q_1}}} = {S_{11}} \cdot {e_{{P_1}}} + {S_{12}} \cdot {e_{{P_2}}} + \cdots + {S_{1m}} \cdot {e_{{P_m}}}\\ {e_{{Q_2}}} = {S_{21}} \cdot {e_{{P_1}}} + {S_{22}} \cdot {e_{{P_2}}} + \cdots + {S_{2m}} \cdot {e_{{P_m}}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \vdots \\ {e_{{Q_L}}} = {S_{L1}} \cdot {e_{{P_1}}} + {S_{L2}} \cdot {e_{{P_2}}} + \cdots + {S_{Lm}} \cdot {e_{{P_m}}} \end{array} \right. $

(10)

式(10)表达为向量矩阵形式为：

$ {E_Q} = S \cdot {E_P} $

(11)

式中, $ {E_Q} = \left[{\begin{array}{*{20}{c}} {{e_{{Q_1}}}}\\ \vdots \\ {{e_{{Q_L}}}} \end{array}} \right], S = \left[{\begin{array}{*{20}{c}} {{S_{11}}}& \cdots &{{S_{1m}}}\\ {}& \vdots &{}\\ {{S_{L1}}}& \cdots &{{S_{Lm}}} \end{array}} \right], {E_P} = \left[{\begin{array}{*{20}{c}} {{e_{{P_1}}}}\\ \vdots \\ {{e_{{P_m}}}} \end{array}} \right]$，式(11)表达了降雨误差与流量计算误差的响应关系.

则面平均雨量误差与流量误差之间的系统响应关系为：

$ {E_Q} = S \cdot {E_P} + \zeta $

(12)

面平均雨量误差的最小二乘估计值为：

$ {{E'}_P} = {\left( {{S^T} \cdot S} \right)^{ - 1}} \cdot {S^T} \cdot {E_Q} $

(13)

式中, ζ为残差向量，E′ _P为降雨误差的动态系统响应曲线估计向量.由式(13)计算得到的面平均雨量误差系列加上原来的面平均雨量系列重新使用预报模型计算，便可得到通过降雨系统响应曲线反馈修正以后的流域出口断面流量过程.

2 应用检验

在本研究中将基于面平均雨量误差修正的实时洪水预报修正方法与新安江模型相结合，对此修正方法进行理想案例和实际流域应用检验.本研究之所以选取新安江模型，是因为新安江模型是一个概念性集总式水文模型，其模型输入项为流域蒸散发能力和单元面平均雨量，而本文所提出的修正方法恰恰是通过修正面平均雨量来提高洪水预报精度，且本文所选取的应用流域为淮河上游王家坝流域，此流域为半湿润流域，符合新安江模型的应用范围.

理想模型和实际流域检验应用采用以下评价指标：

(1) NS(Nash-Sutcliffe)系数：

$ NS = 1 - \sum\limits_{j = 1}^L {{{\left( {{M_{\rm{C}}}\left( j \right) - {M_{\rm{O}}}\left( j \right)} \right)}^2}} /\sum\limits_{j = 1}^L {{{\left( {{M_{\rm{O}}}\left( j \right) - \overline {{M_{\rm{O}}}} } \right)}^2}} $

(14)

(2) 相对误差：

$ \Delta M\left( \% \right) = \left[ {\left( {{M_{\rm{C}}} - {M_{\rm{O}}}} \right)/{M_{\rm{O}}}} \right] \times 100\% $

(15)

(3) 修正后相对误差提高幅度：

$ IRM = \left| {R{M_{{\rm{bu}}}} - R{M_{{\rm{au}}}}} \right|/\left| {R{M_{{\rm{bu}}}}} \right| \times 100\% $

(16)

(4) NS系数提高幅度：

$ INS = \left( {N{S_{{\rm{bu}}}} - N{S_{{\rm{au}}}}} \right)/\left( {1 - N{S_{{\rm{bu}}}}} \right) \times 100\% $

(17)

式中, M_O为统计量的实测值，M_C为统计量的计算值，$ \overline {{M_{\rm{O}}}} $为实测系列的平均值. RM_bu表示某一统计指标修正前的值，RM_au表示统计指标修正之后的值，IRM表示统计指标在修正前后提高的百分比(RM可以是径流深相对误差或者洪峰相对误差)，其值越高表示提高的幅度越大，修正的效果更加明显. NS_bu表示修正前的NS系数，NS_au表示修正后的NS系数.

2.1 理想案例应用检验

理想案例中使用的流域面积为2000 km²，流域中雨量站个数为5，马斯京根河道演算的总河段数也为5.以雨量站为依据，将流域按照平均权重划分为5个流域子单元.本理想案例中5个子单元的面平均雨量以及雨量误差均给定，使用本文提出的修正方法估计面平均雨量误差，然后将估计得到的雨量误差与给定的误差系列进行对比，以此来验证此方法理论依据是否可行.在此理想案例中，5个单元的面平均雨量资料以及每个站的时段雨量的误差均不同.详细的降雨资料、降雨误差和降雨误差的估计值见表 1.

表 1中P₁~P₅分别表示5个雨量站单元的面平均雨量，为了保证给定的雨量数据更加接近一场洪水的降雨过程，其具体降雨数据采用湖北省乌溪沟流域2012年8月15日一场历史洪水降雨资料.采用流域内8个雨量站中的5个雨量站数据来作为本理想案例中5个雨量站单元的面平均雨量. ΔP表示给定的雨量误差，误差值为服从零均值分布的随机误差，误差的控制范围为时段面平均雨量值的±20 %以内. ΔP′表示使用降雨系统响应曲线方法估计的面平均雨量误差值.在本理想案例中，由于流域出口断面的实测流量系列只有1个，因此本研究中估计出的面平均雨量误差为整个流域的时段面平均雨量误差值，而非每个雨量单元的面平均雨量误差值.因此，本理想案例检验中只给定全流域的面平均雨量误差系列.

给定的面平均雨量误差系列与通过降雨系统响应曲线方法估计得到的面平均雨量误差系列对比图如图 3所示.

理想案例修正前、后流量过程线对比结果如图 4所示.图中P表示降雨，Q_c表示流量计算系列，Q_u表示修正后的计算流量，Q_o表示实测流量过程.理想案例中没有实测流量，本理想案例中的实测流量系列是由下式计算得到，即：

$ {Q_{\rm{o}}} = Q\left( {P + \Delta P} \right) + \varepsilon $

(18)

式中, Q(P+ΔP)表示使用面平均雨量系列P+ΔP计算得到的流量系列，ε是零均值的随机流量误差，误差范围控制在对应时段流量的±5 %.

通过图 3和图 4可以看出，在各站之间降雨资料以及每个站的时段雨量的误差均不同的情况下，使用降雨系统响应曲线方法估计的面平均雨量误差和给定的误差系列非常接近，2个系列的NS系数高达0.946.给定面平均雨量误差的总和为-8.4 mm，通过动态系统响应曲线方法估计的降雨误差总和为-8.4 mm，相对误差为0，满足水量平衡.洪峰相对误差由修正前的11.3 %减小到0.7 %，径流深相对误差由修正前的3.48 %降低到0.24 %，洪水流量过程确定性系数由修正前的0.834提高到了0.992.在实测流量存在随机误差的前提下，通过此方法能够准确地估计出面平均雨量误差系列，且能够有效地修正流域出口断面流量过程，提高预报精度，说明此方法理论可行，能够有效地提高预报效果，可以将此方法应用于实际流域检验.

2.2 实际流域应用

对于空间集总式水文模型，流域的大小以及流域内雨量站密度会对面平均雨量计算有较大影响.鉴于此，本研究不仅将此方法应用于王家坝流域进行检验，而且对王家坝流域不同雨量站密度情况的应用效果进行研究，分析不同雨量站密度对此方法修正效果的影响.

王家坝流域面积为30500 km², 王家坝水文站是淮河干流上的重要控制站，此流域对于整个流域的洪水控制和管理起着重要作用.王家坝水文站以上90 %的研究区域位于河南境内的东南部，流域主要由丘陵、山区和平原组成，流域内共有41个雨量站.研究区域的气候为亚热带湿润季风气候，多年平均降雨量为940 mm(1953-2000系列)^[21].本研究中将王家坝流域雨量站点分布情况分为3种，站点水系图如图 5所示.

图 5中b图雨量站分布情况为将原有41个雨量站随机去掉8个，c图为在b图的基础上再随机去掉8个.本文选取了王家坝流域1984-2009年共16场历史洪水进行验证.修正过程中所使用的模型参数是经过王家坝流域16场历史洪水率定检验后获得. 16场历史洪水修正前后计算结果，以及不同雨量站分布情况修正效果提高的幅度统计结果见表 2.

在雨量站个数不减少的情况下，16场历史洪水的预报精度均有明显提高(表 2).平均NS系数由修正前的0.748提高到0.905，提高了0.157.径流深相对误差绝对值的平均值由修正前的9.6 %减小到4.5 %，减小了5.1 %.洪峰相对误差绝对值的平均值由修正前的13.4 %减小到8.5 %，减小了4.9 %.结果显示此修正方法对于洪水径流深和洪峰误差均有明显的修正效果.对于径流深误差的修正幅度明显大于对于洪峰误差的修正幅度(图 6)，说明此方法对于控制洪水预报过程中水量平衡具有明显效果.

在不同的雨量站密度条件下，此方法对于径流深、洪峰流量和NS系数修正效果的对比如图 7所示.

降雨系统响应曲线修正方法在3种不同雨量站密度情况下对于径流深和洪峰流量的平均修正效果以及平均NS系数均有明显提高.雨量站密度越低，其修正效果提高的幅度越大(图 7).雨量站密度低，造成模型输入项不确定性增加，面平均雨量误差较大，通过降雨系统响应曲线方法能够有效降低面平均雨量的不确定，从而提高模型预报精度.因此，在没有误差修正时采用预报模型进行洪水预报计算，洪水预报精度相对偏低，修正的效果就比较明显.随着雨量站个数的增加，雨量站密度增大，面平均雨量精度有所提高，预报精度相对增加，在此情况下修正效果的提高幅度就有所减小.

3种雨量站密度情况下，本修正方法对于径流深修正效果的平均提高幅度均比洪峰流量修正效果的提高幅度大，进一步说明此方法能够更好的修正洪水预报过程中水量平衡误差(图 8).

通过此方法在王家坝流域的应用结果以及在不同雨量站密度情况下的结果分析，表明此方法能够有效提高实际流域洪水预报的精度，且对于雨量站点密度较低的流域修正效果更加显著.

3 结论

本文提出了一种基于降雨系统响应曲线的面平均雨量修正方法，通过估计面平均雨量误差来提高水文模型的预报精度.将水文模型概化为系统，以响应函数对面平均降雨系列的偏微分为降雨系统响应曲线的理论基础.建立了面平均雨量误差与流域出口断面流量误差之间的系统响应关系.利用实测流量与计算流量之间的差值系列，通过最小二乘估计面平均雨量的误差系列，使用修正之后的面平均雨量系列重新计算流域出口端面的出流过程，以此提高模型预报精度.虽然此方法对于提高洪水预报精度具有一定的效果，但由于水文系统毕竟是非线性系统，影响因素较为复杂，这两者之间不是绝对的线性关系，且预报流量过程的误差并不是全部由面平均雨量误差所引起的.因此，如何细化面平均雨量误差和计算流量误差之间的响应机制和量化关系，以及面平均雨量误差和其他误差因素的量化影响还有待进一步研究.

现阶段的研究结果表明，通过修正面平均雨量误差来提高洪水预报精度是一种有效的预报误差修正手段，本修正方法通过理想案例和实际流域验证，理想案例中次洪径流深相对误差由修正前的3.48 %减小到0.24 %，洪水流量过程确定性系数由修正前的0.834提高到了0.992；实际流域应用结果为径流深相对误差绝对值的平均值由修正前的9.6 %减小到4.5 %，平均NS系数由修正前的0.748提高到0.905，结果证明此降雨系统响应曲线方法能够有效地修正面平均雨量误差，提高洪水预报精度.此外，还将此方法应用于王家坝流域的不同雨量站密度情况下，应用结果显示对于不同的雨量站分布情况均有明显的修正效果.对于雨量站密度较低的流域，此方法降低面平均雨量的不确定性效果更加显著，洪水预报精度提高幅度更加明显，对此方法的应用范围和适用条件具有一定的指导意义.

本研究提出的修正方法具有结构简单、不增加模型参数、不改变复杂度以及不增加模型中间变量的优点，是一种有效的提高洪水预报精度的修正方法，此修正方法有待更广泛的应用和检验.此外，本方法现阶段是利用具有较好条件的水文测站流域进行研究，水库流域的入库流量过程具有一定的特殊性，由于水位波动造成反推的入库流量过程呈现锯齿状，从实测流量过程中提取面平均雨量误差信息具有一定的难度，如何通过处理入库流量过程将此方法能够应用于水库流域的面平均雨量误差修正值得进一步的探索研究.本文通过分析研究面平均雨量来修正洪水预报误差，面平均雨量的估计也可以通过降雨等值线图来动态计算得到，在后续的研究中可尝试将雨量等值线图方法与本文中提出的方法进行探讨和对比研究，对于研究面平均雨量误差与洪水预报精度之间的量化关系具有一定的参考价值.