基于洪水預(yù)報(bào)系統(tǒng)誤差反演的多河段聯(lián)合校正方法*

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

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

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

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

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

1 方法原理介紹

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

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

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

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

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

(1)

式中，Qm(t)、Qm(t-Δt)、Qm(t-2Δt)分別為t、t-Δt、t-2Δt時(shí)刻的出口斷面實(shí)測流量；Qc(t)、Qc(t-Δt)、Qc(t-2Δt)分別為t、t-Δt、t-2Δt時(shí)刻的出口斷面預(yù)報(bào)流量；ε(t)、ε(t-Δt)、ε(t-2Δt)分別為t、t-Δt、t-2Δt時(shí)刻的出口斷面洪水預(yù)報(bào)誤差.

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

(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時(shí)刻預(yù)報(bào)誤差ε(t+Δt)的反演方程：

(3)

將估計(jì)的誤差加到原出口斷面預(yù)報(bào)流量Qc(t+Δt)上，則可得到校正后的出口斷面預(yù)報(bào)流量Q′c(t+Δt)：

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

(4)

1.2.2 單河段校正模型 基于動(dòng)力系統(tǒng)反演理論的馬斯京根流量演算的單河段校正模型[14](如圖1c所示)，是考慮河段匯流演算過程中的上斷面入流項(xiàng)、區(qū)間輸入項(xiàng)以及馬斯京根法模型本身誤差，通過建立各項(xiàng)誤差的非線性動(dòng)力系統(tǒng)反演方程(建立原理同單斷面校正模型)，對(duì)各項(xiàng)誤差進(jìn)行校正，以提高河道下斷面洪水預(yù)報(bào)精度.

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

(5)

其中：

(6)

(7)

將估計(jì)的各項(xiàng)誤差代入式(5)中，則可得到校正后的下斷面預(yù)報(bào)值：

(8)

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

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

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

(9)

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

(10)

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

(11)

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

(12)

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

(13)

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

(14)

其中：

(15)

(16)

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

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

(17)

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

2 應(yīng)用與分析

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

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

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

本文采用4個(gè)基本指標(biāo)[21-22]來定量描述洪水預(yù)報(bào)的精度，并借助2個(gè)統(tǒng)計(jì)指標(biāo)來比對(duì)各校正方法的誤差修正能力. 各評(píng)價(jià)指標(biāo)的計(jì)算公式和符號(hào)含義如下：

1)洪峰流量相對(duì)誤差：

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

(18)

式中，Qm,obs為實(shí)測洪峰流量(m3/s),Qm,c為預(yù)報(bào)洪峰流量(m3/s)；δQm為洪峰流量相對(duì)誤差.

2)峰現(xiàn)時(shí)間絕對(duì)誤差：

ΔT=TQm,obs-TQm,c

(19)

式中，TQm,obs為實(shí)測峰現(xiàn)時(shí)間(h),TQm,c為預(yù)報(bào)峰現(xiàn)時(shí)間(h), ΔT為峰現(xiàn)時(shí)間絕對(duì)誤差(h).

3)徑流深相對(duì)誤差：

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

(20)

式中，Robs為實(shí)測次洪徑流深(mm),Rc為預(yù)報(bào)次洪徑流深(mm),δR為徑流深相對(duì)誤差.

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

(21)

5)為了比較不同次洪間多指標(biāo)綜合效果，將(1)～(4)的評(píng)價(jià)指標(biāo)平均值進(jìn)行歸一化處理.

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

(22)

其余指標(biāo)的歸一化算法為：

(23)

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

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

(24)

式中，Qb(t)為t時(shí)刻的基準(zhǔn)預(yù)報(bào)流量(m3/s)；其它變量含義同上.

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

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

圖3 8場次洪統(tǒng)計(jì)指標(biāo)雷達(dá)圖Fig.3 Radar chart of statistical indicators of eight floods

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

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

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

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

3 結(jié)論

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

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

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