洪水预报误差来源多且影响复杂，洪水预报至今仍是水情预报的难点^[1-2]. 实时修正技术是保证实时洪水预报精度的关键^[3]. 现有误差修正方法大多没有考虑流域产汇流过程的物理基础，主要以控制断面误差为基础，将预报对象的误差与模型分离，本质上是“黑箱”分析方法，普遍存在外延性差、预见期损失等问题^[4-5]. 洪水预报误差修正研究中，通常需要考虑模型自变量变化与模型计算结果误差之间的关系. 函数微分中自变量增量与因变量增量之间的微分响应关系，能够较好地描述复杂系统中各影响因素之间的变化关系^[6-7]. 包为民等提出将微分响应关系用于洪水预报误差修正研究，具有物理基础强、不增加参数、不损失预见期等优点^[8]. 现有研究主要利用微分响应方法对流域面平均降雨^[9-10]、面平均产流^[11-13]、面平均状态变量^[14-16]等进行修正，没有考虑洪水预报误差的空间分布差异.

流域降水时空分布和下垫面分布复杂，流域产汇流过程呈高度非线性、时空变化大^[17]. 目前洪水预报误差研究基本上是对各种误差源等权看待，对误差时空变化考虑不足会降低误差修正效果，影响洪水预报精度. 有研究利用雨量站网密度与洪水预报误差分配比例之间的定量关系将洪水预报总误差划分为面雨量输入误差和模型参数误差^[18]，利用各单元计算流量与各单元计算流量和的比例对流域面平均雨量估计值进行各单元分配以减小各单元面雨量误差^[19]. 洪水预报误差修正大多基于流域系统输出信息，虽然该信息能够反映计算过程中各阶段误差的最终结果^[20]，但难以区分反映各阶段误差对洪水预报的影响^[21]. 从流域系统输出的误差序列中辨析误差贡献，考虑误差时空变化进行实时修正有利于提高洪水预报精度.

本文针对洪水预报空间误差修正困难的问题，基于水文系统微分响应关系构建了流域产流分单元修正方法，并比较分析了在不同预见期时流域产流分单元修正、流域面平均产流修正和自回归修正的效果. 该研究对于流域水文要素时空误差修正、洪水预报精度提高具有重要意义和实用价值.

1 方法介绍 1.1 水文系统微分响应修正方法

水文模型的输入、状态变量、参数等任一变量的改变，都会引起流域出口断面流量的变化. 把流域水文模型概化为一个复杂函数，利用微分建立自变量增量与因变量增量之间的响应关系：

$ \text{d}Q=~\frac{\partial f}{\partial X}\text{d}X $

(1)

式中, Q为模型因变量(模型输出)；X为模型自变量(输入变量、状态变量、参数等)；∂f/∂X为变量X的微分响应.

式(1)中自变量增量即为模型影响因素的误差，因变量增量即为模型计算结果的误差. 式(1)的矩阵形式为：

$ \Delta Q=U\Delta X+\xi $

(2)

式中, ΔQ为实测值与计算值的差值序列; U_m×n为雅可比矩阵(微分响应矩阵)，U=[u₁, …, u_j, …, u_n]，列向量$ {{u}_{j}}={{[~\frac{\partial f\left( X, 1 \right)}{\partial {{x}_{j}}}, \frac{\partial f\left( X, 2 \right)}{\partial {{x}_{j}}}, \ldots , \frac{\partial f\left( X, m \right)~}{\partial {{x}_{j}}}]}^{T}} $表示j时刻变量X单位变化所对应的系统微分响应；m为用于误差修正的实测流量时段数; n为变量的修正时段数; ξ为随机误差项，是符合零均值分布特点的白噪声向量.

采用正则化最小二乘法求解式(2)，得变量X的误差估计量为^[16]：

$\Delta X={{({{U}^{T}}U+\beta I)}^{-1}}{{U}^{T}}\Delta Q $

(3)

式中, I为单位矩阵, β为正则化系数.

修正后的变量X′为：

$ {X}'=X+\Delta X $

(4)

流域水文系统是一个非线性系统，考虑非线性问题线性化的影响，采用逐步迫近进行反演求解^[22]：

${{X}^{(k+1)}}={{X}^{(k)}}+\Delta {{X}^{(k)}} $

(5)

式中, 上标k为迫近次数, X^(k)为第k次迭代的变量估计值, ΔX^(k)为第k次迭代的变量误差估计值.

将修正后的变量重新输入模型进行计算可得到修正后的模型结果. 基于水文系统微分响应关系建立误差修正方法，利用了产汇流模型结构的物理基础，概念清晰，计算简便.

1.2 流域产流分单元微分响应修正方法

流域水文模型中每个子过程可以概化为一个系统，多个密切联系的子过程也可以概化为一个系统. 新安江模型包括蒸散发计算、产流计算、分水源计算和汇流计算^[23]，可以将蒸散发计算和产流计算作为一个系统(记为E_R系统)，水源划分和汇流计算作为一个系统(记为R_Q系统). 产流计算是对降雨在数量上的再分配，产流量是影响水量平衡的重要因素，在洪水预报中起着重要作用. 对于产流微分响应修正，可通过R_Q系统利用流域出口流量反演修正产流(图 1).

采用新安江模型进行流域洪水预报，通常需进行流域划分. 假设流域划分为s个单元流域，流域出口流量对各单元流域产流的响应可以表达为：

$ Q=f({{R}_{i}}, t)\ (i=1\sim s) $

(6)

式中, i为单元流域编码, s为单元流域数目, R_i=[r_{i, 1}, …, r_{i, j}, …, r_{i, n}]^T为第“i”单元流域产流.

考虑误差空间分布，构建流域产流分单元微分响应修正(又称流域产流分块修正或时空产流修正)，记为SDR，如图 2a所示，计算步骤如下：

(1) 计算j时刻各单元流域产流量R_i (i=1~s)，将i单元流域j时刻计算产流量r_{i, j}加上1个单位(其他时刻其他单元流域产流量保持不变)代入R_Q系统计算得到新的流量系列，减去原来的流量系列，得到i单元流域U矩阵j列u_j.

(2) 在不同时刻(j=1~n)重复步骤(1)即可得到i单元流域U矩阵所有列.

(3) 计算产流误差估计值ΔR，根据i单元流域面积占全流域面积的比例将估计的产流误差分配给i单元流域，得到i单元产流误差修正量为ΔR′ _i =ΔR·η_i，其中η_i为i单元流域面积占全流域面积的比例.

(4) 重复步骤(1)~(3)，可得到各单元流域U矩阵及相应的产流修正量.

已有研究采用流域面平均产流修正^[11-13]，记为AMR，如图 2b所示，计算流程如下：

(1) 计算j时刻各单元流域产流量R_i(i=1~s)，将各单元流域j时刻计算产流量加上1个单位(其他时刻产流量保持不变)代入R_Q系统计算得到新的流量系列，减去原来的流量系列，得到全流域U矩阵j列u_j.

(2) 在不同时刻(j=1~n)重复步骤(1)就得到U所有列.

(3) 计算产流量误差估计值ΔR.

(4) 根据单元流域面积占全流域面积的比例将产流误差估计值分配给各单元流域，得到各单元产流误差修正量ΔR′ _i.

流域产流分单元微分响应修正(SDR)和流域面平均产流微分响应修正(AMR)的不同点在于微分响应矩阵的构建和产流估计误差的分配. AMR在求解微分响应矩阵时所有单元流域产流量均增加1个单位，所求得的微分响应反映了流域出口流量对流域平均产流变化的响应. SDR分别计算不同单元流域对应的微分响应，通过各单元流域微分响应提取流域出口断面流量中包含的空间(单元流域)信息，增加了信息利用. AMR微分响应矩阵为2维数组U_m×n，SDR微分响应矩阵为3维数组U_m×n×s. AMR是将产流估计误差按面积比例分配给各个单元流域，而SDR是将利用单元流域微分响应求得的误差估计值分配给相应的单元流域.

2 应用检验 2.1 实例概况

将构建的微分响应误差修正方法用于修正三水源新安江模型计算产流. (1)比较SDR和AMR的修正效果. (2)比较SDR、AMR与二阶自回归(AR)在采用不同预见期时的修正效果.

选择具有不同雨量站点密度和不同洪水特点的淮河大坡岭流域和闽江七里街流域为研究区域. 流域图如图 3，流域信息如表 1. 采用泰森多边形法，将大坡岭流域划分为13个单元流域，将七里街流域划分为41个单元流域，每个单元流域有一个雨量站.

采用三水源新安江模型进行流域洪水模拟. 大坡岭流域采用1987—2009年14场洪水资料，从起涨到洪峰出现的平均历时为36.3 h，其中10场用于新安江模型参数率定，4场用于检验. 七里街流域采用1988—1999年20场洪水资料，从起涨到洪峰出现的平均历时为77.9 h，其中15场用于新安江模型参数率定，5场用于检验.

微分响应修正采用逐步迫近求解，以误差均方差为目标函数判断是否进行下一次迭代. 通过微分响应修正时所使用实测流量资料长度设置了一系列不同预见期，预见期越长所用实测流量时段数越少. 采用径流深相对误差(RRE)、洪峰相对误差(RPE)、误差均方差(RMSE)和Nash-Sutcliffe系数(NSE)^[24]评估计算结果.

2.2 结果分析

淮河大坡岭流域和闽江七里街流域的洪水模拟结果分别见表 2和表 3. 表中$\overline{|RRE|}$为径流深相对误差绝对值的平均值，$\overline{|RPE|}$为洪峰相对误差绝对值的平均值，$\overline{|RMSE|}$为误差均方差平均值，$\overline{|NSE|}$为NSE系数平均值. 部分洪水过程如图 4和图 6，其中XAJ为采用新安江模型未经误差修正的原模型计算结果，SDR为采用流域产流分单元修正的计算结果，AMR为采用流域面平均产流修正的计算结果. 图 5和图 7分别展示了不同预见期时SDR、AMR和AR的修正结果.

由表 2可知，在淮河大坡岭流域，新安江模型经过SDR和AMR修正后计算精度提高了. SDR修正效果优于AMR，经SDR修正后$\overline{|NSE|}$由0.913提高到0.961，$\overline{|RMSE|}$由102.2 m³/s减小为58.7 m³/s，$\overline{|RRE|}$和$\overline{|RPE|}$也减小. 由图 4可见，采用SDR修正后的洪水计算过程线更好地接近实测过程线，洪峰部分拟合较好. 由图 5可见，随着预见期增大，SDR修正效果总体优于AMR修正，SDR和AMR修正效果优于AR修正. 当预见期大于2 h时，采用SDR和AMR修正的RMSE小于AR修正结果，NSE系数大于AR修正结果. SDR和AMR的修正效果比AR修正稳定，不会随着预见期增长而明显变差.

由表 3可知，在闽江七里街流域，采用SDR和AMR修正后新安江模型计算效果提高了，SDR修正效果总体上优于AMR. 与AMR结果相比，SDR修正后$\overline{|NSE|}$由0.863提高到0.927，$\overline{|RMSE|}$由526.4 m³/s减小到367.1 m³/s. 由图 6可见，经SDR修正后的洪水计算过程线与实测过程线拟合较好. 由图 7可见，与AR修正相比，当预见期大于3 h时采用SDR修正的结果明显较好，RMSE较小，NSE系数较大. 随着预见期增大，AR修正效果明显变差，SDR和AMR的修正效果优于AR，SDR修正效果优于AMR.

2.3 讨论

大坡岭流域和七里街流域所采用的历史洪水从起涨到洪峰出现的平均历时差别较大，能用于误差修正的实测流量个数不同，微分响应修正均获得了较好的修正效果. 在七里街流域部分雨量站雨量资料缺乏，采用邻近站资料替代，雨量资料时空误差较大，通过模型计算引起的产流时空误差也较大. 对于时空误差较大的洪水，AMR修正效果有限，SDR修正效果较好. 实例中微分响应修正以误差均方差作为判断条件进行逐步迫近反演求解，SDR修正和AMR修正的平均迭代次数均小于10次. 经误差修正后，误差均方差减小了，NSE系数提高了，但并不能保证所有洪水的洪峰和径流深模拟精度都提高，可以考虑采用兼顾洪峰和洪水过程的多目标函数作为判断条件.

微分响应修正是基于流域出口流量与被修正量之间的微分响应关系建立的反馈修正方法，理论上可用于修正模型任一自变量，实际应用中通常用于修正主要误差. 在洪水预报中可用于修正降雨、产流、参数、状态变量等，修正效果取决于被修正变量是否是主要误差来源. 产流计算受诸多误差影响，模型输入、状态变量等误差经蒸散发计算和产流计算会传播累积影响计算产流量，产流修正很有必要. 在大坡岭流域预见期小于2 h，在七里街流域预见期小于4 h时，AR修正效果优于SDR，主要是因为AR修正是直接对所有误差影响下的流域出口流量进行修正. 但AR修正效果取决于输出结果误差之间的相关性，随着预见期增大修正效果明显变差.

结合新安江模型的产流微分响应修正是将分水源和汇流部分作为响应系统，利用流域出口流量误差信息反馈修正产流量，将修正后的产流重新输入模型进行后续的分水源、坡地汇流和河道汇流计算，进而修正流域出口断面流量. 流域出口断面流量是所有误差影响的最终结果，对其进行时空误差辨识很难. SDR修正建立了各单元流域产流与流域出口断面流量之间的微分响应关系，利用各单元流域微分响应关系可以从流域出口断面流量中提取有效的空间信息用于产流量修正. 基于微分响应的流域分单元误差修正思路可为分布式流域水文模型误差修正研究提供一种有效途径，进一步可通过水文系统各要素引起的微分响应关系研究水文要素时空变化，成果应用可以有效提高实时洪水预报精度.

3 结论

1) 基于微分响应构建了流域产流分单元修正方法(SDR)，利用汇流系统物理基础建立了各单元流域产流误差与流域出口流量误差之间的响应关系，采用正则化最小二乘法结合逐步迫近求解产流修正估计量，对产流误差修正量进行时空分配实现流域产流分单元修正.

2) 在淮河大坡岭流域和闽江七里街流域结合新安江模型进行方法检验，结果表明：采用产流微分响应修正后的模型效果优于原水文模型效果；SDR效果总体优于AMR；随着预见期增大，SDR和AMR的修正效果优于AR.

3) 微分响应分单元修正能通过水文系统从流域出口信息中提取有效的空间信息用于误差修正，考虑了误差的空间分布，使用简单，不改变原模型结构，可与不同水文模型结合用于修正不同时空误差，对于时空误差研究具有重要意义.