基于洪水预报系统误差反演的多河段联合校正方法*

黄一昕，梁忠民，胡义明，李彬权，王 军

(河海大学水文水资源学院，南京 210098)

洪水预报是一项重要的防洪减灾非工程措施[1]. 虽然洪水预报理论和技术得到了长足发展[2-3]，但仍面临挑战：一方面，现行的水文模型理论与方法尚难以对流域水文物理过程进行精准地模拟和预报，另一方面，气候变化与高强度人类活动的影响导致水文事件发生的时空格局复杂多变，给洪水预报无形中增加了难度[4-5]. 因此，在实时洪水预报中，实时校正技术作为保障和提升洪水预报精度的最后一道关口，必不可少.

实时校正方法的研究较多[6-9]，从实用角度，一般以对单河段的校正居多，即只利用河段上下游断面的实时信息进行下断面预报误差的修正和更新. 例如：王船海等[10]提出了基于卡尔曼滤波校正技术的单一河道水动力学模型，通过局部校正对全河道预报起带动作用；常露等[11]建立了基于K-最近邻非参数实时校正模型的河道洪水预报系统，可对具有行蓄洪区的河道流域进行模拟与校正；徐兴亚等[12]建立了基于粒子滤波数据同化算法的河道洪水实时概率预报模型，不仅可以直接校正水位，同时也可以有效地校正流量和糙率；高益辉等[13]提出了将自适应自回归模型与河道水流演进基本方程相结合的多点联合校正方法，可对并联的单河段汇流系统进行实时校正；梁忠民等[14]提出了基于动力系统反演理论的马斯京根流量演算误差校正方法，对单河段洪水演算过程中马斯京根法的3类误差源，分别构建3种误差演算方程并进行误差校正. 上述这些方法的特点是只利用本河段上下游断面的实测信息，所以计算简单、应用方便，但也因其忽略了整个河道干支流其它断面的实测和校正信息，有时并不能取得理想效果，且校正的有效预见期受限. 实际的河流水系结构一般都比较复杂，断面或河段间存在水力联系，下游某一断面的预报误差受上游多断面的共同影响，因此，在实时校正中，可以充分利用河系中河段之间的信息联系，考虑上游关联断面对下游校正断面的影响，降低误差在汇流过程中的累积传播效应，以提高校正精度.

为此，本文引入马斯京根法汇流演算矩阵方程[15]描述河系上下游之间的水力联系，并与动力系统误差反演方程[14]结合，提出一种适合复杂河流系统的多河段联合校正方法，并在大渡河上游流域进行了示例应用.

1 方法原理介绍

1.1 复杂河流系统

图1a所示是一个自然的河流系统，系统内的长河段可划分为n个子河段(从1级河流至n级河流)和n+1个河段断面(从第0个断面至第n个断面，其中：第0个断面为河源断面，第n个断面为流域出口断面).由于该系统各级河流的上下游断面之间存在水力联系，各断面流量Qi(i=1,…,n)的来源包含两部分：①第i个河段的上断面入流，即第i-1个断面的控制面积以上降水径流Qi-1；②第i个河段的区间入流qi.相应地，各断面流量的预报误差也由两部分组成：①上断面入流的预报误差；②区间入流的预报误差. 因此，对于一个如图1a的复杂河流系统，其出口断面流量预报误差的大小，是由上游多个区间和断面的预报误差共同作用与影响决定的. 而对于图1b和1c所示的概化的单断面/单河段系统，认为流域下游出口断面仅受该断面或上游断面的影响. 其中，图1b是将图1a中的复杂多河段河流系统简化成一个单断面系统，出口断面流量为Q；图1c是将图1a中的复杂多河段河流系统简化成一个单河段系统，下断面出口断面流量为QOut，上断面河源断面流量为QIn，河段区间入流为q.

1.2 单断面/单河段校正模型

1.2.1 单断面校正模型 基于动力系统误差反演的单断面校正模型，是直接对流域出口断面的预报误差进行校正，以提高洪水预报精度. 动力系统是指状态随时间变化的系统，可由特定的方程(如微分方程)描述. 动力系统反演则是指通过观测资料反求动力系统方程的过程[16-17]. 若以出口断面流量系列为研究对象(如图1b所示)，以一个受3个变量影响的系统为例.

依据出口断面流量的历史预报值和实测值，计算洪水预报误差：

(1)

式中，Qm(t)、Qm(t-Δt)、Qm(t-2Δt)分别为t、t-Δt、t-2Δt时刻的出口断面实测流量；Qc(t)、Qc(t-Δt)、Qc(t-2Δt)分别为t、t-Δt、t-2Δt时刻的出口断面预报流量；ε(t)、ε(t-Δt)、ε(t-2Δt)分别为t、t-Δt、t-2Δt时刻的出口断面洪水预报误差.

根据预报误差时间系列的相依特性，建立误差反演方程：

(2)

式中，参数a1、a2、a3、a4、a5、a6、a7、a8、a9、a10，b1、b2、b3、b4、b5、b6、b7、b8、b9、b10，c1、c2、c3、c4、c5、c6、c7、c8、c9、c10是根据观测资料得到的方程特解.

将微分dε(t)/dt、dε(t-Δt)/dt、dε(t-2Δt)/dt转换成差分，得到t+Δt时刻预报误差ε(t+Δt)的反演方程：

(3)

将估计的误差加到原出口断面预报流量Qc(t+Δt)上，则可得到校正后的出口断面预报流量Q′c(t+Δt)：

Q′c(t+Δt)=Qc(t+Δt)+ε(t+Δt)

(4)

1.2.2 单河段校正模型 基于动力系统反演理论的马斯京根流量演算的单河段校正模型[14](如图1c所示)，是考虑河段汇流演算过程中的上断面入流项、区间输入项以及马斯京根法模型本身误差，通过建立各项误差的非线性动力系统反演方程(建立原理同单断面校正模型)，对各项误差进行校正，以提高河道下断面洪水预报精度.

考虑区间入流的单河段马斯京根法演算方程[18-19]可表示为：

(5)

其中：

(6)

(7)

将估计的各项误差代入式(5)中，则可得到校正后的下断面预报值：

(8)

1.3 多河段联合校正模型

本文提出的多河段联合校正模型，是采用马斯京根法矩阵方程描述多河段汇流过程，基于动力系统反演理论对各河段的区间入流误差进行演算，通过对各河段预报误差联合校正，最终降低河道出口断面洪水预报总误差.

1.3.1 马斯京根法矩阵方程 对于多站点、多断面的河流系统，应考虑河道汇流时的流量演进和误差传递. 如图1a所示的具有n个子河段的河流系统，表示成马斯京根法的矩阵形式(考虑区间入流)[15]：

(9)

以3段为例，如果各子河段的演算参数相等，则有：

(10)

1.3.2 基于动力系统反演理论的区间入流误差演算 考虑河道汇流演进过程中的区间入流误差以及第1个河段的上断面入流误差，基于动力系统反演理论的对预报误差演算. 第1个河段的上断面入流误差可根据历史预报值和实测值计算，而由于区间来水不好测量，区间入流误差系列根据该河段下断面实测流量扣除区间预报流量在下断面的响应来推算. 即：

(11)

参照式(3)建立各区间入流误差(第1个河段的上断面入流误差也可认为是该河段的旁侧入流误差)的动力系统反演方程：

(12)

1.3.3 多河段流量预报误差联合校正 将式(12)代入式(10)，建立基于马斯京根法矩阵方程的流量预报误差演算方程：

(13)

对式(13)左端第一个矩阵求逆矩阵，然后等式两边左乘这个逆矩阵，整理后得到：

(14)

其中：

(15)

(16)

式(14)写成向量形式如下：

Q(t+Δt)=ΗQ(t)+Ρ[q(t+Δt)+ε(t+Δt)]

(17)

式中，Η=(ΑΤΑ)-1ΑΤΒ和Ρ=(ΑΤΑ)-1ΑΤ均为系数矩阵，Q(t+Δt)为t+Δt时刻断面流量预报值向量，Q(t)为t时刻断面流量实测值向量，q(t+Δt)为t+Δt时刻区间流量预报值向量，ε(t+Δt)为t+Δt时刻区间流量预报误差预测值向量.

2 应用与分析

2.1 研究区域及资料

大渡河是中国长江支流岷江的正源，丹巴水文站是大渡河上游的一个重要节点. 丹巴站以上流域年平均降水量为600～700 mm，年平均径流量为400～500 mm，流域控制面积为52763 km2，约占大渡河流域总面积的68%. 丹巴站径流的变化可以直接反映大渡河上游区的流量变化，也可以决定下游区的来水变化. 因此，保证丹巴站获得精准的洪水预报意义重大. 由于丹巴站以上流域的地形地质条件复杂，河道迂回曲折，支流较多，且站点布设困难，雨量站网密度稀疏，因而难以准确描绘该流域的下垫面机制，难以测得流域的面降雨量和各河道断面和区间的流量，这些都导致在丹巴站的洪水预报不准确. 但丹巴站以上流域内有不少水文测站，可基于站点信息，采用本文提出的多河段联合校正方法对预报洪水进行校正.

图2所示为丹巴站以上流域的地理位置、河系概况及测站分布. 对该流域内的河道、支流、站点分布等进行合理概化. 概化后以丹巴水文站为研究站点，以丹巴站以上河流系统为研究河道，在干流选取4个代表性水文站(日部、足木足、大金、丹巴)，将长河段划分为3个子河段(日部-足木足、足木足-大金、大金-丹巴)，对本文方法进行应用检验. 根据该流域2009-2016年的水文气象观测数据对各站点进行流量预报. 根据河道洪水传播特性，预报时间步长Δt取24 h. 根据已有研究成果[20]，确定采用的马法参数为：x=0.4，K=25，C0=0.074，C1=0.815，C2=0.111.选用8场洪水资料对各河段入流误差的反演方程参数进行率定和检验，其中4场用于率定，4场进行验证，参数率定结果见表1.

图2 大渡河丹巴站以上河流系统Fig.2 River system of Dadu River above Danba station

表1 各河段入流误差的反演方程系数向量的率定结果Tab.1 Calibration results of the coefficient vectors of the inversion equations for the inflow errors of each reach

2.2 评价指标

本文采用4个基本指标[21-22]来定量描述洪水预报的精度，并借助2个统计指标来比对各校正方法的误差修正能力. 各评价指标的计算公式和符号含义如下：

1)洪峰流量相对误差：

δQm=[(Qm,obs-Qm,c)/Qm,obs]×100%

(18)

式中，Qm,obs为实测洪峰流量(m3/s),Qm,c为预报洪峰流量(m3/s)；δQm为洪峰流量相对误差.

2)峰现时间绝对误差：

ΔT=TQm,obs-TQm,c

(19)

式中，TQm,obs为实测峰现时间(h),TQm,c为预报峰现时间(h), ΔT为峰现时间绝对误差(h).

3)径流深相对误差：

δR=[(Robs-Rc)/Robs]×100%

(20)

式中，Robs为实测次洪径流深(mm),Rc为预报次洪径流深(mm),δR为径流深相对误差.

4)确定性系数：

(21)

5)为了比较不同次洪间多指标综合效果，将(1)～(4)的评价指标平均值进行归一化处理.

确定性系数的归一化算法为：

(22)

其余指标的归一化算法为：

(23)

式中，W为归一化处理前的指标值，Wmin为指标的全局最小值，Wmax为指标的全局最大值，W*为归一化处理后的指标值.

不同的评价指标具有不同的量纲和量纲单位，因此不具有可比性. 若要综合比较不同条件下的不同指标，需要对指标进行归一化处理，以消除指标之间的量纲影响. 归一化后的各评价指标处于同一量级，适合进行综合评价(comprehensive evaluation, CE)[23]. 一般采用雷达图的图形数值相结合的方式来综合展示所有归一化指标的对比结果. 在雷达图中，单个归一化指标数值越大，越接近1，表示模拟效果越好；对于不同方法的多个归一化指标，雷达图所围面积越大，表示该方法模拟效果越好.

(24)

式中，Qb(t)为t时刻的基准预报流量(m3/s)；其它变量含义同上.

不同校正方法对原始基准预报的提升效果不一样，基准系数可用来评价校正前后的预报精度提升程度和校正方法的相对好坏.BE≤0，表示校正后的预报结果表现较差，校正方法不可取；BE>0，表示校正结果较优，且BE越大，校正效果越明显.

2.3 结果分析及讨论

图3 8场次洪统计指标雷达图Fig.3 Radar chart of statistical indicators of eight floods

图4提供了8场次洪的基准系数BE的对比. 从图4中可以看出：① 多河段联合校正模型和单河段校正模型的所有BE全部大于0，2种校正方法对所有场次洪水过程的预报值都有提升，提升率为100%. ② 多河段联合校正模型比单河段校正模型的提升效果更明显，基准系数能多提升0.1以上. ③ 多河段联合校正模型不仅能提高所有场次洪水的预报效果，还能稳定地提升至一个较高的精度，最低的BE为0.08，最高的BE达0.79. 由此可推断，随着划分的河段数目增加(如：丹巴站以上河道从单段/1段划分成3段)，利用的断面和区间资料信息增多，则根据误差反演措施所带来的校正模型的误差修正作用增强. 相当于每增加一个水文测站，河道上就增加一个误差“传感器”，河道划分越精细，“传感器”越多. “传感器”可以及时地发现误差在什么地方什么时刻出现、误差有多大，校正模型就能据此进行必要的误差修正. 所以，出口断面受上游不确定性影响减弱，最终预报结果的精度提高.

图4 2种校正方法的基准系数对比 Fig.4 Comparison of Benchmark coefficients of two calibration methods

选取8场洪水中的第20110702号和第20161006号洪水作为示例，查看误差修正效果. 图5显示了这2场洪水校正前后的预报流量对比，图5a为丹巴水文站的2011年7月2日发生的洪水，图5b为2016年10月6日发生的洪水. 从图5中可以看出：① 2种校正模型均能有效地校正洪水原始预报结果，将2场不合格的洪水预报变为合格的预报. ② 相较而言，多河段联合校正模型综合表现最佳，不论对洪峰还是洪水过程的校正，效果更为理想，这有利于在实际应用中，准确地模拟流量过程线，尤其是捕捉峰值，以降低洪水风险. ③ 多河段联合校正模型对单峰洪水(如20161006号洪水)和复峰洪水(如20110702号洪水)的预报误差的校正效果相当. 结果在一定程度上说明，本文提出的基于误差反演的多河段联合校正模型，将各河段洪水预报误差演算，并联合进行误差校正的做法更为合理，可进一步提高流域出口断面的洪水预报精度.

图5 两场次洪的预报及校正结果Fig.5 Forecast and correction results of two floods

3 结论

1)本研究将马斯京根矩阵方程和动力系统反演方程相结合，提出洪水预报误差反演的多河段联合校正方法. 该方法考虑了复杂河系中各断面在空间上的水力联系和预报误差在时程上的传递规律，充分利用上游多断面实测和校正信息对下游预报断面的预报误差进行校正. 该方法物理意义明确，理论先进.

2)在大渡河流域的应用结果表明，对不同年月发生的场次洪水均能取得稳定有效的校正效果，显著提升洪水预报精度. 相比于单河段校正方法，多河段联合校正方法在校正能力上整体更优，能保证洪峰、洪量及整个洪水过程的预报精度.

3)本文对大渡河流域构建的洪水预报误差反演方程属于三元三次项的非线性模型，对其它流域未必适合，但文中采用的方法具有普适性，可根据应用流域的实际误差传递规律进行优化调整，以获得更好的校正效果. 因此，本文方法也受限于洪水预报能力和误差多时段外延效果，需进一步拓展研究.