黃一昕，梁忠民，胡義明，李彬權(quán)，王 軍

(河海大學水文水資源學院，南京 210098)

洪水預報是一項重要的防洪減災非工程措施[1]. 雖然洪水預報理論和技術得到了長足發(fā)展[2-3]，但仍面臨挑戰(zhàn)：一方面，現(xiàn)行的水文模型理論與方法尚難以對流域水文物理過程進行精準地模擬和預報，另一方面，氣候變化與高強度人類活動的影響導致水文事件發(fā)生的時空格局復雜多變，給洪水預報無形中增加了難度[4-5]. 因此，在實時洪水預報中，實時校正技術作為保障和提升洪水預報精度的最后一道關口，必不可少.

實時校正方法的研究較多[6-9]，從實用角度，一般以對單河段的校正居多，即只利用河段上下游斷面的實時信息進行下斷面預報誤差的修正和更新. 例如：王船海等[10]提出了基于卡爾曼濾波校正技術的單一河道水動力學模型，通過局部校正對全河道預報起帶動作用；常露等[11]建立了基于K-最近鄰非參數(shù)實時校正模型的河道洪水預報系統(tǒng)，可對具有行蓄洪區(qū)的河道流域進行模擬與校正；徐興亞等[12]建立了基于粒子濾波數(shù)據(jù)同化算法的河道洪水實時概率預報模型，不僅可以直接校正水位，同時也可以有效地校正流量和糙率；高益輝等[13]提出了將自適應自回歸模型與河道水流演進基本方程相結(jié)合的多點聯(lián)合校正方法，可對并聯(lián)的單河段匯流系統(tǒng)進行實時校正；梁忠民等[14]提出了基于動力系統(tǒng)反演理論的馬斯京根流量演算誤差校正方法，對單河段洪水演算過程中馬斯京根法的3類誤差源，分別構(gòu)建3種誤差演算方程并進行誤差校正. 上述這些方法的特點是只利用本河段上下游斷面的實測信息，所以計算簡單、應用方便，但也因其忽略了整個河道干支流其它斷面的實測和校正信息，有時并不能取得理想效果，且校正的有效預見期受限. 實際的河流水系結(jié)構(gòu)一般都比較復雜，斷面或河段間存在水力聯(lián)系，下游某一斷面的預報誤差受上游多斷面的共同影響，因此，在實時校正中，可以充分利用河系中河段之間的信息聯(lián)系，考慮上游關聯(lián)斷面對下游校正斷面的影響，降低誤差在匯流過程中的累積傳播效應，以提高校正精度.

為此，本文引入馬斯京根法匯流演算矩陣方程[15]描述河系上下游之間的水力聯(lián)系，并與動力系統(tǒng)誤差反演方程[14]結(jié)合，提出一種適合復雜河流系統(tǒng)的多河段聯(lián)合校正方法，并在大渡河上游流域進行了示例應用.

1 方法原理介紹

1.1 復雜河流系統(tǒng)

圖1a所示是一個自然的河流系統(tǒng)，系統(tǒng)內(nèi)的長河段可劃分為n個子河段(從1級河流至n級河流)和n+1個河段斷面(從第0個斷面至第n個斷面，其中：第0個斷面為河源斷面，第n個斷面為流域出口斷面).由于該系統(tǒng)各級河流的上下游斷面之間存在水力聯(lián)系，各斷面流量Qi(i=1,…,n)的來源包含兩部分：①第i個河段的上斷面入流，即第i-1個斷面的控制面積以上降水徑流Qi-1；②第i個河段的區(qū)間入流qi.相應地，各斷面流量的預報誤差也由兩部分組成：①上斷面入流的預報誤差；②區(qū)間入流的預報誤差. 因此，對于一個如圖1a的復雜河流系統(tǒng)，其出口斷面流量預報誤差的大小，是由上游多個區(qū)間和斷面的預報誤差共同作用與影響決定的. 而對于圖1b和1c所示的概化的單斷面/單河段系統(tǒng)，認為流域下游出口斷面僅受該斷面或上游斷面的影響. 其中，圖1b是將圖1a中的復雜多河段河流系統(tǒng)簡化成一個單斷面系統(tǒng)，出口斷面流量為Q；圖1c是將圖1a中的復雜多河段河流系統(tǒng)簡化成一個單河段系統(tǒng)，下斷面出口斷面流量為QOut，上斷面河源斷面流量為QIn，河段區(qū)間入流為q.

1.2 單斷面/單河段校正模型

1.2.1 單斷面校正模型 基于動力系統(tǒng)誤差反演的單斷面校正模型，是直接對流域出口斷面的預報誤差進行校正，以提高洪水預報精度. 動力系統(tǒng)是指狀態(tài)隨時間變化的系統(tǒng)，可由特定的方程(如微分方程)描述. 動力系統(tǒng)反演則是指通過觀測資料反求動力系統(tǒng)方程的過程[16-17]. 若以出口斷面流量系列為研究對象(如圖1b所示)，以一個受3個變量影響的系統(tǒng)為例.

依據(jù)出口斷面流量的歷史預報值和實測值，計算洪水預報誤差：

(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時刻的出口斷面洪水預報誤差.

根據(jù)預報誤差時間系列的相依特性，建立誤差反演方程：

(2)

式中，參數(shù)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是根據(jù)觀測資料得到的方程特解.

將微分dε(t)/dt、dε(t-Δt)/dt、dε(t-2Δt)/dt轉(zhuǎn)換成差分，得到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 單河段校正模型 基于動力系統(tǒng)反演理論的馬斯京根流量演算的單河段校正模型[14](如圖1c所示)，是考慮河段匯流演算過程中的上斷面入流項、區(qū)間輸入項以及馬斯京根法模型本身誤差，通過建立各項誤差的非線性動力系統(tǒng)反演方程(建立原理同單斷面校正模型)，對各項誤差進行校正，以提高河道下斷面洪水預報精度.

考慮區(qū)間入流的單河段馬斯京根法演算方程[18-19]可表示為：

(5)

其中：

(6)

(7)

將估計的各項誤差代入式(5)中，則可得到校正后的下斷面預報值：

(8)

1.3 多河段聯(lián)合校正模型

本文提出的多河段聯(lián)合校正模型，是采用馬斯京根法矩陣方程描述多河段匯流過程，基于動力系統(tǒng)反演理論對各河段的區(qū)間入流誤差進行演算，通過對各河段預報誤差聯(lián)合校正，最終降低河道出口斷面洪水預報總誤差.

1.3.1 馬斯京根法矩陣方程 對于多站點、多斷面的河流系統(tǒng)，應考慮河道匯流時的流量演進和誤差傳遞. 如圖1a所示的具有n個子河段的河流系統(tǒng)，表示成馬斯京根法的矩陣形式(考慮區(qū)間入流)[15]：

(9)

以3段為例，如果各子河段的演算參數(shù)相等，則有：

(10)

1.3.2 基于動力系統(tǒng)反演理論的區(qū)間入流誤差演算 考慮河道匯流演進過程中的區(qū)間入流誤差以及第1個河段的上斷面入流誤差，基于動力系統(tǒng)反演理論的對預報誤差演算. 第1個河段的上斷面入流誤差可根據(jù)歷史預報值和實測值計算，而由于區(qū)間來水不好測量，區(qū)間入流誤差系列根據(jù)該河段下斷面實測流量扣除區(qū)間預報流量在下斷面的響應來推算. 即：

(11)

參照式(3)建立各區(qū)間入流誤差(第1個河段的上斷面入流誤差也可認為是該河段的旁側(cè)入流誤差)的動力系統(tǒng)反演方程：

(12)

1.3.3 多河段流量預報誤差聯(lián)合校正 將式(12)代入式(10)，建立基于馬斯京根法矩陣方程的流量預報誤差演算方程：

(13)

對式(13)左端第一個矩陣求逆矩陣，然后等式兩邊左乘這個逆矩陣，整理后得到：

(14)

其中：

(15)

(16)

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

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

(17)

式中，Η=(ΑΤΑ)-1ΑΤΒ和Ρ=(ΑΤΑ)-1ΑΤ均為系數(shù)矩陣，Q(t+Δt)為t+Δt時刻斷面流量預報值向量，Q(t)為t時刻斷面流量實測值向量，q(t+Δt)為t+Δt時刻區(qū)間流量預報值向量，ε(t+Δt)為t+Δt時刻區(qū)間流量預報誤差預測值向量.

2 應用與分析

2.1 研究區(qū)域及資料

大渡河是中國長江支流岷江的正源，丹巴水文站是大渡河上游的一個重要節(jié)點. 丹巴站以上流域年平均降水量為600～700 mm，年平均徑流量為400～500 mm，流域控制面積為52763 km2，約占大渡河流域總面積的68%. 丹巴站徑流的變化可以直接反映大渡河上游區(qū)的流量變化，也可以決定下游區(qū)的來水變化. 因此，保證丹巴站獲得精準的洪水預報意義重大. 由于丹巴站以上流域的地形地質(zhì)條件復雜，河道迂回曲折，支流較多，且站點布設困難，雨量站網(wǎng)密度稀疏，因而難以準確描繪該流域的下墊面機制，難以測得流域的面降雨量和各河道斷面和區(qū)間的流量，這些都導致在丹巴站的洪水預報不準確. 但丹巴站以上流域內(nèi)有不少水文測站，可基于站點信息，采用本文提出的多河段聯(lián)合校正方法對預報洪水進行校正.

圖2所示為丹巴站以上流域的地理位置、河系概況及測站分布. 對該流域內(nèi)的河道、支流、站點分布等進行合理概化. 概化后以丹巴水文站為研究站點，以丹巴站以上河流系統(tǒng)為研究河道，在干流選取4個代表性水文站(日部、足木足、大金、丹巴)，將長河段劃分為3個子河段(日部-足木足、足木足-大金、大金-丹巴)，對本文方法進行應用檢驗. 根據(jù)該流域2009-2016年的水文氣象觀測數(shù)據(jù)對各站點進行流量預報. 根據(jù)河道洪水傳播特性，預報時間步長Δt取24 h. 根據(jù)已有研究成果[20]，確定采用的馬法參數(shù)為：x=0.4，K=25，C0=0.074，C1=0.815，C2=0.111.選用8場洪水資料對各河段入流誤差的反演方程參數(shù)進行率定和檢驗，其中4場用于率定，4場進行驗證，參數(shù)率定結(jié)果見表1.

圖2 大渡河丹巴站以上河流系統(tǒng)Fig.2 River system of Dadu River above Danba station

表1 各河段入流誤差的反演方程系數(shù)向量的率定結(jié)果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個統(tǒng)計指標來比對各校正方法的誤差修正能力. 各評價指標的計算公式和符號含義如下：

1)洪峰流量相對誤差：

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

(18)

式中，Qm,obs為實測洪峰流量(m3/s),Qm,c為預報洪峰流量(m3/s)；δQm為洪峰流量相對誤差.

2)峰現(xiàn)時間絕對誤差：

ΔT=TQm,obs-TQm,c

(19)

式中，TQm,obs為實測峰現(xiàn)時間(h),TQm,c為預報峰現(xiàn)時間(h), ΔT為峰現(xiàn)時間絕對誤差(h).

3)徑流深相對誤差：

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

(20)

式中，Robs為實測次洪徑流深(mm),Rc為預報次洪徑流深(mm),δR為徑流深相對誤差.

4)確定性系數(shù)：

(21)

5)為了比較不同次洪間多指標綜合效果，將(1)～(4)的評價指標平均值進行歸一化處理.

確定性系數(shù)的歸一化算法為：

(22)

其余指標的歸一化算法為：

(23)

式中，W為歸一化處理前的指標值，Wmin為指標的全局最小值，Wmax為指標的全局最大值，W*為歸一化處理后的指標值.

不同的評價指標具有不同的量綱和量綱單位，因此不具有可比性. 若要綜合比較不同條件下的不同指標，需要對指標進行歸一化處理，以消除指標之間的量綱影響. 歸一化后的各評價指標處于同一量級，適合進行綜合評價(comprehensive evaluation, CE)[23]. 一般采用雷達圖的圖形數(shù)值相結(jié)合的方式來綜合展示所有歸一化指標的對比結(jié)果. 在雷達圖中，單個歸一化指標數(shù)值越大，越接近1，表示模擬效果越好；對于不同方法的多個歸一化指標，雷達圖所圍面積越大，表示該方法模擬效果越好.

(24)

式中，Qb(t)為t時刻的基準預報流量(m3/s)；其它變量含義同上.

不同校正方法對原始基準預報的提升效果不一樣，基準系數(shù)可用來評價校正前后的預報精度提升程度和校正方法的相對好壞.BE≤0，表示校正后的預報結(jié)果表現(xiàn)較差，校正方法不可?。籅E>0，表示校正結(jié)果較優(yōu)，且BE越大，校正效果越明顯.

2.3 結(jié)果分析及討論

圖3 8場次洪統(tǒng)計指標雷達圖Fig.3 Radar chart of statistical indicators of eight floods

圖4提供了8場次洪的基準系數(shù)BE的對比. 從圖4中可以看出：① 多河段聯(lián)合校正模型和單河段校正模型的所有BE全部大于0，2種校正方法對所有場次洪水過程的預報值都有提升，提升率為100%. ② 多河段聯(lián)合校正模型比單河段校正模型的提升效果更明顯，基準系數(shù)能多提升0.1以上. ③ 多河段聯(lián)合校正模型不僅能提高所有場次洪水的預報效果，還能穩(wěn)定地提升至一個較高的精度，最低的BE為0.08，最高的BE達0.79. 由此可推斷，隨著劃分的河段數(shù)目增加(如：丹巴站以上河道從單段/1段劃分成3段)，利用的斷面和區(qū)間資料信息增多，則根據(jù)誤差反演措施所帶來的校正模型的誤差修正作用增強. 相當于每增加一個水文測站，河道上就增加一個誤差“傳感器”，河道劃分越精細，“傳感器”越多. “傳感器”可以及時地發(fā)現(xiàn)誤差在什么地方什么時刻出現(xiàn)、誤差有多大，校正模型就能據(jù)此進行必要的誤差修正. 所以，出口斷面受上游不確定性影響減弱，最終預報結(jié)果的精度提高.

圖4 2種校正方法的基準系數(shù)對比 Fig.4 Comparison of Benchmark coefficients of two calibration methods

選取8場洪水中的第20110702號和第20161006號洪水作為示例，查看誤差修正效果. 圖5顯示了這2場洪水校正前后的預報流量對比，圖5a為丹巴水文站的2011年7月2日發(fā)生的洪水，圖5b為2016年10月6日發(fā)生的洪水. 從圖5中可以看出：① 2種校正模型均能有效地校正洪水原始預報結(jié)果，將2場不合格的洪水預報變?yōu)楹细竦念A報. ② 相較而言，多河段聯(lián)合校正模型綜合表現(xiàn)最佳，不論對洪峰還是洪水過程的校正，效果更為理想，這有利于在實際應用中，準確地模擬流量過程線，尤其是捕捉峰值，以降低洪水風險. ③ 多河段聯(lián)合校正模型對單峰洪水(如20161006號洪水)和復峰洪水(如20110702號洪水)的預報誤差的校正效果相當. 結(jié)果在一定程度上說明，本文提出的基于誤差反演的多河段聯(lián)合校正模型，將各河段洪水預報誤差演算，并聯(lián)合進行誤差校正的做法更為合理，可進一步提高流域出口斷面的洪水預報精度.

圖5 兩場次洪的預報及校正結(jié)果Fig.5 Forecast and correction results of two floods

3 結(jié)論

1)本研究將馬斯京根矩陣方程和動力系統(tǒng)反演方程相結(jié)合，提出洪水預報誤差反演的多河段聯(lián)合校正方法. 該方法考慮了復雜河系中各斷面在空間上的水力聯(lián)系和預報誤差在時程上的傳遞規(guī)律，充分利用上游多斷面實測和校正信息對下游預報斷面的預報誤差進行校正. 該方法物理意義明確，理論先進.

2)在大渡河流域的應用結(jié)果表明，對不同年月發(fā)生的場次洪水均能取得穩(wěn)定有效的校正效果，顯著提升洪水預報精度. 相比于單河段校正方法，多河段聯(lián)合校正方法在校正能力上整體更優(yōu)，能保證洪峰、洪量及整個洪水過程的預報精度.

3)本文對大渡河流域構(gòu)建的洪水預報誤差反演方程屬于三元三次項的非線性模型，對其它流域未必適合，但文中采用的方法具有普適性，可根據(jù)應用流域的實際誤差傳遞規(guī)律進行優(yōu)化調(diào)整，以獲得更好的校正效果. 因此，本文方法也受限于洪水預報能力和誤差多時段外延效果，需進一步拓展研究.