刘可新，李 匡，王紫微，梁犁丽，徐海卿

(中国水利水电科学研究院北京中水科水电科技开发有限公司，北京 100038)

1 研究背景

水文模型主要研究的是时不变的离线系统，大都采用观测到的历史水文资料，先确定好模型参数，然后用于未来的洪水预报。这种方式在实时洪水预报系统中时常得不到满意的结果[1]。这时往往要采用误差修正以取得较好的预报结果，因此误差修正技术一直是洪水预报研究的重点，其中包括传统的自回归模型及其改进方法[2,3]、卡尔曼滤波技术[4]，后来的集合卡尔曼滤波[5]、神经网络技术等[6]。近些年来，诸多学者在这方面做了大量工作，相继提出了抗差分析[7-8]、综合修正方法[9]等多种修正方法，均获得了较好的修正效果。然而，传统的误差修正方法基本上都存在物理基础不强、损失预见期和修正效果不理想等缺陷[10]。鉴于此，文献[10]应用系统响应理论，将单位线作为线性时不变系统，引入实时洪水预报修正中，建立了一种向信息源头追溯的反馈修正模型，继文献[10]之后，文献[11]将动态系统响应曲线引入误差修正中，提出了产流误差动态系统响应曲线(Dynamic System Response Curve, DSRC)修正方法(以下简称DSRC方法)，对产流量进行修正，该方法具有物理基础且不损失预见期，应用于实际流域，效果十分显著，随后系统响应修正方法得以发展与应用[12-17]。

但DSRC方法有时不稳定，没有考虑到流量测量误差对产流修正稳定性的影响，忽视了产流修正系列对流量测量误差的敏感性，因而修正后的产流系列时常出现“震荡”现象，本质上是反演计算的不适定问题。鉴于此，本文深入分析了流量测量误差对产流修正系列的影响程度，通过正则化方法大大降低了产流修正系列对流量测量误差的敏感性，提出了产流动态系统响应正则化(DSRR)修正方法(以下简称DSRR方法) 。

2 修正方法介绍

2.1 DSRC方法讨论

流域出口的实测流量包含信息最全面，计算过程中各阶段的误差最终都将反映在流量上，因此，提取这些误差并追溯其产生源头是误差修正的关键所在[18]。DSRC方法假设误差来源于产流以上部分，以产流系统响应曲线作为误差追溯其源头的桥梁，从而将误差合理分配到各时段产流并进行修正。

如前所述，DSRC方法是应用中的反问题，这类问题往往是不适定的[19]。所谓适定性是指解的存在性，解的唯一性和解的稳定性。应用检验表明DSRC方法往往存在不稳定问题，这一定程度上影响了修正效果。所谓稳定性是指解对观测值的连续依赖性，下面通过简单例子解释该性质。

线性方程组：

现在考虑观测向量的微小误差对解向量的影响，考察线性方程组：

解得x=(1 1)T。

从以上例子看出，观测向量出现微小变化时，解向量却变化很大，这就是解不稳定的现象，其中方程组的系数矩阵A称为病态矩阵，病态矩阵的条件数往往很大(A的条件数为40 002)。反演问题中往往存在病态矩阵，并且观测向量的误差难以避免，因此即使模型合理有效，也很难解得满意的解向量。

DSRC方法同样存在相应稳定性问题，解向量对应于产流修正系列，观测向量对应于流量误差系列，实测流量存在误差是难以避免的，因此流量误差系列并非真值，存在一定的波动，于是产流误差系列往往得不到满意的结果，其不稳定表现为产流修正系列时程上的“振荡”现象。DSRC方法对产流修正系列的稳定性要求较高，因为误差修正是外延过程，产流修正系列不稳定会严重影响修正效果。鉴此，提出了产流动态系统响应正则化修正方法，将不稳定的反演问题转换成邻近的稳定的反演问题，降低了产流修正系列对流量误差系列的敏感程度，提高了产流修正系列的稳定性，从而优化了误差修正效果。

2.2 DSRR方法介绍

近些年来，正则化方法已有较多的研究成果[20,21]。所谓正则化方法，是指用一些与原问题相邻近的适定问题的解去逼近原问题的解，因此如何建立“邻近”的问题是正则化方法的核心。Tikhcmov正则化方法是较成熟的一种正则化方法，其主要思想是通过增加一个罚函数项构建适定的邻近问题。DSRC方法在有些情况下存在稳定性问题，此时流量观测误差对产流修正系列影响较严重，因此本次研究将Tikhcmov正则化方法应用于DSRC方法中，提出了产流动态系统响应正则化修正方法。

为研究方便，以新安江模型产流以下的部分作为系统，产流系列作为系统输入，出口流量作为系统输出，系统结构如图1。

图1 系统示意

图1是新安江模型示意图，其中E0表示实际蒸散发，P表示降雨，R表示产流量，Q表示出口流量，TLOC表示三层蒸发，NSM表示蓄满产流，SRC表示分水源，FT表示汇流，BC表示反演计算，虚线框中的部分表示系统。

根据文献[11]，只考虑产流对出口流量的影响，那么图1所示系统可以表示为：

Q(t)=Q(R,θ,t)

(1)

式中：R=[r1,r2,…,rn]T表示产流系列；Q(R,θ,t)即实测流量系列Q0。

根据动态系统响应理论，建立产流与流量之间的关系，可用式(2)表示：

Q(R,θ,t)=Q(RC,θ,t)+UΔR+E

(2)

式中：RC=[rC1,rC2,…,RCn]T为计算得到的产流量初始系列；ΔR=[Δr1,Δr2,Δr3,…]T为待估计的产流修正系列；Q(RC,θ,t)为计算流量初始系列；E=[e1,e2,…,eL]T为流量观测随机误差项；U为动态系统响应矩阵，其表示各时段产流量与出口流量之间的关系。

DSRC方法根据无约束最小二乘原理，所求ΔR满足式(3)：

Q(RC,θ,t)-UΔR]

(3)

于是有：

ΔR=(UTU)-1UT[Q(R,θ,t)-Q(RC,θ,t)]

(4)

实践发现，UTU矩阵有时十分病态甚至趋于奇异，这时流量观测误差对ΔR的影响很大，甚至导致估计严重失真。

鉴于此，引进了Tikhcmov正则化方法，提出了DSRR方法，在式(3)中加入罚函数项，如式(5)：

(5)

式中:α为权重系数。

经推导得式(6)：

ΔR=(UTU+αI)-1UT[Q(R,θ,t)-Q(RC,θ,t)]

(6)

式中：I为单位矩阵。

DSRR方法较DSRC方法更稳定表现在两个方面，一方面式(5)是约束最小二乘法的目标函数，由两部分组成，误差项和罚函数项，其中误差项表示修正后计算流量与实测流量的接近程度，其值越小表示越接近，罚函数项表示稳定程度，若要满足其值较小，那么ΔR不会出现“振荡”现象；另一方面从式(6)中可以看出，当UTU病态甚至奇异时，加上αI项可解决该问题，降低ΔR对流量观测误差的敏感程度，消除估计严重失真的现象。α作为权重系数，在式(5)中作用是平衡误差项与罚函数项，在式(6)中作用是保证降低UTU病态程度，而又满足与原问题邻近。

修正后计算产流系列R′C为：

R′C=RC+ΔR

(7)

将R′C代入模型重新计算出流过程：

Q′=Q(R′C,θ,t)

(8)

Q′为修正后的计算流量系列。

3 实际流域应用与检验

3.1 评定指标的选取

为了较好地反映修正效果及修正的稳定性，本文采用以下4个指标作为评价指标，其中洪量相对误差、洪峰相对误差和Ns系数反应修正效果，条件数反应修正的稳定性。各指标的计算公式如下：

(1)洪量相对误差：

(9)

(2)洪峰相对误差：

(10)

(3)Nash-Sutcliffe系数：

(11)

(4)条件数：

cond(A)2=‖A‖2‖A-1‖2

(12)

3.2 实际案例检验

本次研究将滩坑流域作为实验流域，滩坑流域位于浙江省瓯江流域的小溪支流上，流域面积为3 330 km2，流域内有15个雨量站。滩坑流域受亚热带季风影响，气候温和，雨量丰沛，四季变化明显，日照充足。流域径流以降水补给为主，降水量年际变化大，丰枯水年相互交替出现，流域年平均降水量一般在1 500～2 100 mm，多年平均水面蒸发量为969.9 mm。流域内降水量空间分布不均，年分配有显著差异，全年大致可分为3-4月春雨，5-6月梅雨，7-9月台风雷雨和10-次年2月干季4个时期。流域多年平均流量120 m3/s，多年平均年径流量37.8 亿m3，径流年内分配不均匀，全年以6月份来水量最大，占全年水量的19.7%。流域水系及站网布设见图2。

图2 滩坑流域示意图

本次研究以滩坑流域1985—2004年19场洪水作为检验资料，使用新安江模型对流域的历史洪水进行计算和修正。为全面深入分析DSRR方法的修正效果，选取其中一场洪水(洪号：31980618)作为典型，分别采用DSRC方法和DSRR方法，以不损失预见期为条件，对该场洪水进行修正，其效果见图3和图4。由图3看出，DSRR方法修正后的时段产流平稳变化，没有出现陡涨陡落的情况，并且修正前后产流过程线形状基本保持一致，说明修正效果符合实际情况；而应用DSRC方法，修正后时段产流虽然前期比较平稳，但后期出现了较明显的不稳定现象，而且随着时段的增加，不稳定越来越严重，最后甚至出现了“振荡”现象，表现为相邻时段产流量交替变化，修正效果明显不符合产流规律，脱离实际情形。出现这一状况的原因考虑有以下几点：一方面矩阵UTU条件数(1.2×109)过大，是病态矩阵，使得产流修正系列ΔR对流量观测误差过于敏感，这就意味着微小的观测误差同样会造成ΔR剧烈变化，而流量观测误差是难以避免的，于是采用DSRC方法会导致震荡现象，而采用DSRR方法，矩阵UTU+αI的条件数(780)大大减小了，这就降低了ΔR对流量观测误差的敏感性，增强了方法的稳定性；另一方面DSRR方法中加入了罚函数项，约束产流修正系列，使其在一定范围内变化；此外，以不损失预见期为条件，使得可用的实测流量信息有限，有效信息不足以弥补由于测量误差产生的干扰信息的影响。由4中不难看出，DSRR方法得到的洪水过程线比较光滑，而DSRC方法则出现了局部的陡涨陡落，这与两种方法产流修正效果的稳定程度一致；对洪水的修正效果，DSRR方法整体上优于DSRC方法，这是由于新方法使产流修正更加稳定，修正结果更加符合实际情况，而DSRC方法虽注重洪水过程线拟合程度，但忽视了流量测量误差对产流修正结果的严重影响。

图3 DSRR方法与DSRC方法产流修正效果比较

图4 DSRR方法与DSRC方法典型洪水修正效果对比

滩坑流域19场洪水模拟及修正结果见表1。由表1看出，DSRR方法的修正效果十分显著，各场洪水的预报精度明显提高，平均洪量相对误差、洪峰相对误差和峰现时差分别降低了3.8%、8.2%和1.0h，平均确定性系数提高到0.970；DSRR方法的应用效果明显优于DSRC方法，各项指标均有所改善；同时发现采用DSRC方法，每场洪水条件数均很大，平均达到1.1×109，说明该方法修正结果的稳定性较差，而DSRR方法将条件数降低了7个数量级，平均为773，极大地降低了产流修正系列对流量测量误差的敏感程度，提高了修正结果的稳定性；此外还发现，由于各场洪水的降雨蒸发等水文条件不同，采用DSRR方法使修正结果趋于稳定时，不同洪水其条件数各异。

表1 滩坑流域19场洪水实时修正结果

本文同时对权重系数α做了初步探讨，在不损失预见期条件下，绘制了不同洪水α与产流修正个数的关系图(图5)，同时在产流修正个数逐渐增加条件下，绘制了同一场洪水(洪号：31980618)α的变化趋势图(图6)。由图5看出，不同洪水的α与产流修正个数具有一定的变化一致性，但趋势并不显著，也就是说对于不同的洪水， 的取值与产流修正个数有一定正相关关系，但还要视情况而定，这主要是由于不同洪水水文条件各异所致；由图6看出，同一场水，α随产流修正个数增加而变大，趋势十分显著。

4 结 语

本文深入分析了DSRC方法，提出了DSRR方法，对DSRR方法中α的取值做了初步探讨，并对两种方法的修正效果进行了比较，得出以下几点结论：

(1)DSRC方法在有些情况下，产流修正系列对流量测量误差的敏感性较强，难以保证修正结果的稳定性；

(2)产流修正系列对流量测量误差的敏感性可用条件数量化表示，条件数越大敏感性往往越大；

(3)DSRR方法降低了产流修正系列对流量测量误差的敏感性，其修正结果的稳定性及修正效果明显优于DSRC方法；DSRR方法中α的取值视洪水而定，不同洪水产流修正个数对α的取值有一定的参考作用，同一场洪水，随着产流修正时段的增加，α逐渐变大。

□

[1] 田 雨, 雷晓辉, 蒋云钟, 等. 洪水预报实时校正技术研究综述[J]. 人民黄河, 2011,33(3):25-28.

[2] 包为民. 水文预报[M]. 北京: 中国水利水电出版社, 2006.

[3] 周 轶, 菅浩然, 李致家, 等. 基于递推最小二乘改进算法的洪水预报模型研究[J]. 河海大学学报,2007,35(1):77-80.

[4] 葛守西, 程海云, 李玉荣. 水动力学模型卡尔曼滤波实时校正技术[J]. 水利学报, 2005,36(6):1-9.

[5] CLARK M P, RUPP D E,WOODS R A, et al. Hydrological data assimilation with the ensemble Kalman filter: use of streamflow observation to update states in a distributed hydrological model[J]. Advances in Water Rsources, 2008,31:1 309-1 324.

[6] 朱星明, 卢长娜, 王如云, 等.基于人工神经网络的洪水水位预报模型[J].水利学报,2005,36(7):806-811.

[7] 包为民, 嵇海祥, 胡其美, 等. 抗差理论及在水文学中的应用[J]. 水科学进展, 2003,14(4):528- 532.

[8] 赵 超, 洪华生, 包为民, 等. 实时洪水抗差预报系统研究[J]. 水文, 2008,28(2):27-29.

[9] 瞿思敏, 包为民. 实时洪水预报综合修正方法初探[J]. 水科学进展, 2003,14(2):167-171.

[10] 包为民, 司 伟, 沈国华, 等. 基于单位线反演的产流误差修正[J]. 水科学进展, 2012,23(3):315-322.

[11] 司 伟, 包为民, 瞿思敏. 洪水预报产流误差的动态系统响应曲线修正方法[J]. 水科学进展, 2013,24(4):497-503.

[12] Wei Si, Weimin Bao, and Hoshin V. Gupta. Updating real-time flood forecasts via the dynamic system response curve method[J]. Water Resources Research, 2015,(7).

[13] 张小琴, 刘可新, 包为民, 等. 产流误差比例系数的系统响应修正方法[J]. 水科学进展, 2014,25(6):789-796.

[14] 包为民, 阙家骏, 赖善证, 等. 洪水预报自由水蓄量动态系统响应修正方法[J]. 水科学进展, 2015,26(3):365-371.

[15] 杨姗姗, 包为民, 杨小强, 等.微分响应在降雨误差修正中的应用[J]. 中国农村水利水电, 2015,(1):75-79.

[16] 刘可新, 包为民, 李佳佳, 等. 基于系统响应理论的分水源误差修正[J]. 水电能源科学, 2014,32(11):44-47.

[17] 刘可新, 包为民, 赖善证, 等. 动态系统响应曲线修正方法在乌溪沟流域的应用[J]. 中国农村水利水电, 2014,(12):24-26.

[18] 包为民. 洪水预报信息利用问题与讨论[J]. 水文, 2006, 26(2):18-21.

[19] 顾勇为, 归庆明, 边少锋, 等. 地球物理反问题中两种正则化方法的比较[J]. 武汉大学学报·信息科学版, 2005,30(3):238-241.

[20] 董国志. 反问题的正则化方法及其计算[D]. 长沙: 湖南师范大学，2012.

[21] 杨 帆.三类不适定问题的正则化方法研究[D].兰州:兰州大学，2014.