王家彪，趙建世，雷曉輝，王 浩，，魏雋煜，廖衛(wèi)紅

（1.清華大學(xué)，水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室，北京 100084；2.中國水利水電科學(xué)研究院，北京 100038）

1 研究背景

河道水流計(jì)算模型是洪水管理與水庫調(diào)度的重要工具。隨著流域數(shù)字化與洪水管理需求的不斷提高，水動(dòng)力模型已廣泛應(yīng)用于河道匯流演算以獲取更實(shí)時(shí)精細(xì)的水流動(dòng)態(tài)信息[1-2]。采用水動(dòng)力模型計(jì)算河道水流存在諸多不確定性，其誤差來源主要包括上下游邊界條件，區(qū)間入流以及河道糙率信息。對(duì)于很多河道尤其是山區(qū)小河道而言，其區(qū)間支流雖然可通過現(xiàn)場(chǎng)勘測(cè)或遙感等方式確定具體位置[3]，但入流過程往往缺少實(shí)測(cè)資料，來自于支流入流的誤差嚴(yán)重降低了水流計(jì)算精度。準(zhǔn)確分析出支流入流過程對(duì)提高河道水流計(jì)算精度意義重大。

采用水文學(xué)方法匯流演算時(shí)，通常由水文模型產(chǎn)流模擬或直接由河道入流與出流水量差推算出區(qū)間總?cè)肓鳎肓飨盗袝r(shí)間間隔一般在小時(shí)以上，其精度和時(shí)間尺度都不足以滿足洪水實(shí)時(shí)管理的需求[4-7]。采用水動(dòng)力學(xué)方法進(jìn)行匯流演算時(shí)，區(qū)間入流常作為圣維南方程組的反問題進(jìn)行求解[8]。已有部分研究借助于水文模型等方式獲取初始區(qū)間入流信息，并基于水動(dòng)力模型和數(shù)據(jù)同化方法對(duì)初始入流誤差進(jìn)行動(dòng)態(tài)修正[5，9-10]。如，張薔等[11]與李光熾等[12]基于離散的圣維南方程推導(dǎo)出區(qū)間入流校正量的解析表達(dá)式，然后由主干河道實(shí)測(cè)水位、流量過程動(dòng)態(tài)修正入流過程；吳曉玲等[13]認(rèn)為在河道水流實(shí)時(shí)校正中得到的系統(tǒng)噪聲均值是由區(qū)間入流不精確引起，由此基于卡爾曼濾波方法對(duì)區(qū)間入流初始值動(dòng)態(tài)修正。雖然上述工作取得了一定成效，但由于區(qū)間水量不平衡誤差同時(shí)受到上下游邊界條件、河道地形、模型參數(shù)等多方面影響，該類方法難以確定具體的誤差來源，推算的區(qū)間入流不夠可靠。另外，從總水量平衡進(jìn)行修正，無法直接確定某條支流的單獨(dú)入流過程。為此，有學(xué)者通過直接從下游往上游反算水流的方式推算區(qū)間支流入流，包括馬斯京根反演[4，14-15]和反向求解圣維南方程組等[16-17]。相對(duì)而言，這類水流反算結(jié)果穩(wěn)定性較差[14]，對(duì)下游邊界條件誤差敏感，推算的入流過程仍具有很大不確定性[18]。

本文將數(shù)據(jù)同化應(yīng)用于支流反分析研究中，采用集合卡爾曼濾波方法（EnKF）對(duì)水動(dòng)力模型反算得到的支流入流結(jié)果進(jìn)行校正。EnKF方法已廣泛用于河道水文、水動(dòng)力模型實(shí)時(shí)校正研究中[19-24]，能動(dòng)態(tài)校正流量水位誤差、修正模型糙率系數(shù)，從而改善河道水流計(jì)算精度。本文首先以支流入流斷面為界將河道劃分為兩段，分別采用一維河道水動(dòng)力模型進(jìn)行正、反向水流計(jì)算，根據(jù)入流斷面流量差初步估算支流入流過程。然后，構(gòu)建監(jiān)測(cè)斷面滯時(shí)矩陣，計(jì)算水流擾動(dòng)傳播時(shí)間，確定用于校正支流入流的監(jiān)測(cè)斷面流量觀測(cè)值?？紤]到EnKF校正的結(jié)果仍可能存在較大誤差，研究二次運(yùn)用EnKF對(duì)首次校正結(jié)果進(jìn)一步校正，以提高支流計(jì)算精度。

2 水流計(jì)算

本文通過求解圣維南方程組計(jì)算非恒定水流，水流計(jì)算方程如下：

式中：t為計(jì)算時(shí)間，s；x為縱向坐標(biāo)，m；A為過流斷面面積，m2；B為水面寬度，m；Q為過流流量，m3/s；q為區(qū)間入流量，m3/（s·m）；Z為斷面水位，m；Sf為摩阻坡度，由下式（3）得到：

式中：n為河床糙率系數(shù)；R為過流斷面水力半徑，m。

方程采用四點(diǎn)偏心隱格式差分方法求解（圖1），離散方程（1）和（2）。

圖1 普里斯曼差分格式離散

對(duì)任意離散河段j，可得到線性方程組，

式中：Qj、Zj分別為第j斷面t+1時(shí)刻的流量和水位值；系數(shù)以 及b2j-1，b2j()i=1，2，3，4，j=1，2，…，N由上一時(shí)步流量水位確定，計(jì)算公式如下：

其中對(duì)于支流匯入情形，qj+1/2Δx等于t+1時(shí)刻支流入流QL。

在給定上下游邊界條件和初始條件后，可采用追趕法求解方程組（4）、（5），實(shí)現(xiàn)水流的動(dòng)態(tài)計(jì)算[25]。

3 支流入流計(jì)算模型

3.1 支流入流初步計(jì)算假定河道斷面j處有支流匯入，且在下游斷面G處有水位流量觀測(cè)站點(diǎn)，如圖2所示。

圖2 區(qū)間支流匯入斷面與觀測(cè)斷面

由于無實(shí)測(cè)的支流入流資料，因此需通過干流河道實(shí)測(cè)流量信息對(duì)支流入流進(jìn)行估算。以支流匯入斷面j為界，將河道分成1-j，j-N兩段分別進(jìn)行水流正、反向演算，由斷面j處支流入流前后的流量差推算支流入流。其中水流反向演算：以G斷面為上游、j斷面為下游，逆著水流方向（圖2虛箭頭）計(jì)算經(jīng)支流匯流后j斷面的流量。由于在反算支流時(shí)，j斷面的水位、流量過程未知，因此可通過近似擬合該斷面處的Q-Z關(guān)系曲線作為分段河道的邊界條件，具體步驟如下：

（1）以首斷面1處實(shí)測(cè)流量Q1(t)、末斷面N處實(shí)測(cè)水位ZN(t)（或該處Q-Z關(guān)系曲線）為邊界條件，對(duì)全河段1-N進(jìn)行水流模擬計(jì)算，得到支流匯入斷面j處流量、水位過程。擬合j斷面Q-Z關(guān)系曲線，作為1-j河段水流計(jì)算下邊界條件。其中，Q-Z關(guān)系曲線采用如下常用指數(shù)形式[26]：

式中，α、β和γ為經(jīng)驗(yàn)擬合系數(shù)。

當(dāng)j斷面Q-Z關(guān)系曲線已知時(shí)，此步驟可省略。

（2）以首斷面1處實(shí)測(cè)入流量、下游斷面j處Q-Z關(guān)系曲線為邊界條件，計(jì)算支流匯入前的斷面j處流量Qj(t)和水位過程Zj(t)。

（3）以監(jiān)測(cè)斷面G處實(shí)測(cè)流量過程QG(t)為上邊界條件，步驟（2）計(jì)算的斷面j處水位過程Zj(t)（反向水流計(jì)算對(duì)邊界條件敏感，若將步驟（1）得到的斷面j處Q-Z關(guān)系曲線作為邊界條件，反向結(jié)果不穩(wěn)定）近似為下邊界條件，從下游至上游反向計(jì)算支流匯入后的斷面j處流量過程-Q′j(t)（此時(shí)，G-j河段計(jì)算的流量皆為負(fù)值）。

（4）由步驟（3）計(jì)算的斷面j處流量過程-Q′j(t)減去步驟（2）計(jì)算的斷面j處流量過程Qj(t)，可初步計(jì)算出斷面j處支流入流流量：

3.2 監(jiān)測(cè)斷面滯時(shí)矩陣支流匯入河道后隨干流傳播至監(jiān)測(cè)斷面，其誤差產(chǎn)生的影響需一定傳播時(shí)間才能在斷面G處觀測(cè)到，存在滯時(shí)響應(yīng)的問題[27]。由數(shù)據(jù)同化方法對(duì)支流入流進(jìn)行校正時(shí)，采用的觀測(cè)斷面G處實(shí)測(cè)值應(yīng)盡可能地包含當(dāng)前時(shí)刻支流誤差擾動(dòng)信息。用于校正的監(jiān)測(cè)值QG(tG)與支流入流值QL(tj)存在如下時(shí)間滯后關(guān)系：

式中，Trj-G為斷面j處支流擾動(dòng)傳播至斷面G處時(shí)間（即滯時(shí)），Trj-G的大小可通過對(duì)tj時(shí)刻支流入流值QL(tj)按比例施加正態(tài)分布擾動(dòng)后，由斷面G處出現(xiàn)最大波動(dòng)的時(shí)間來確定：

每一tj時(shí)刻匯入的支流都會(huì)對(duì)干流產(chǎn)生水流擾動(dòng)信息，并且在tG時(shí)刻傳播至監(jiān)測(cè)斷面，由此可得到監(jiān)測(cè)斷面對(duì)支流入流的滯時(shí)矩陣：

式中：Tr(t)為t時(shí)刻的滯時(shí)矩陣；nG為監(jiān)測(cè)斷面數(shù)，本文考慮單個(gè)監(jiān)測(cè)斷面，即nG=1；T為總計(jì)算時(shí)步。

3.3 集合卡爾曼濾波校正集合卡爾曼濾波方法（EnKF）由Monte-Carlo抽樣生成集合樣本模擬校正變量均值和誤差協(xié)方差[28]，并通過不斷引入觀測(cè)數(shù)據(jù)對(duì)模型預(yù)測(cè)值進(jìn)行動(dòng)態(tài)校正。根據(jù)3.1節(jié)支流初步計(jì)算結(jié)果以及3.2節(jié)監(jiān)測(cè)斷面滯時(shí)矩陣，運(yùn)用EnKF對(duì)支流初值進(jìn)行校正，步驟如下（圖3）：

圖3 基于EnKF的支流反分析流程

（1）支流入流初值。由3.1節(jié)方法計(jì)算t時(shí)刻支流入流QL0(t)，并由滯時(shí)矩陣（10）確定用于t時(shí)刻校正的干流觀測(cè)值

（2）集合樣本生成。利用Monte-Carlo抽樣方法生成t時(shí)刻集合樣本：

式中：m為集合樣本數(shù)，取為50；隨機(jī)擾動(dòng)，Rt為模型預(yù)測(cè)誤差協(xié)方差；負(fù)上標(biāo)代表校正前狀態(tài)變量。

（4）支流入流校正。

式中：正上標(biāo)代表校正后狀態(tài)變量；Kt為增益系數(shù)，滿足：

（5）支流入流計(jì)算。計(jì)算t時(shí)刻支流入流量QL()t：

（6）轉(zhuǎn)入下一時(shí)步t+1時(shí)刻。

采用EnKF數(shù)據(jù)同化時(shí)，初值對(duì)校正結(jié)果影響很大，在初始集合均值與真值偏差較大時(shí)，EnKF不能得到最優(yōu)分析[29]。在應(yīng)用上述方法進(jìn)行支流反分析時(shí)，若初始估算結(jié)果誤差較大時(shí)，EnKF校正結(jié)果很可能仍存在很大誤差。此時(shí)可將步驟（1）中由水動(dòng)力模型反算的支流入流初值替換為EnKF首次校正結(jié)果，降低初值誤差，對(duì)支流入流進(jìn)行二次校正，以進(jìn)一步改善支流推算精度。

3.4 評(píng)價(jià)指標(biāo)本文采用相關(guān)系數(shù)R2，相對(duì)均方根誤差和Nash-Sutcliffe效率系數(shù)（NSE）三項(xiàng)誤差指標(biāo)評(píng)價(jià)模型模擬效果，計(jì)算公式如下：

式中：Mt為模擬值；Ot為實(shí)測(cè)值。

4 模型應(yīng)用

4.1 理想案例假設(shè)有一棱柱形梯形斷面河道，底寬10 m，邊坡2，底坡1/10000，糙率系數(shù)0.017。河道全長20 km，在距離上游入口8 km處有一支流匯入，在距離上游入口16 km處設(shè)有水文觀測(cè)站點(diǎn)。干流觀測(cè)斷面和支流入流流量過程如圖4所示（在理想案例中，支流入流不考慮觀測(cè)誤差，監(jiān)測(cè)流量實(shí)際由水動(dòng)力模型正向模擬得到）。

圖4 理想案例干流監(jiān)測(cè)斷面與支流流量過程

圖5 支流峰值時(shí)刻擾動(dòng)在監(jiān)測(cè)斷面處的響應(yīng)過程

圖6 理想案例支流反分析計(jì)算結(jié)果

4.1.1 支流初步計(jì)算和滯時(shí)矩陣確定 首先，根據(jù)3.1節(jié)所述方法計(jì)算支流入流初值。然后，對(duì)入流初值人為施加不同程度的擾動(dòng)（式（9）中h取不同值），確定監(jiān)測(cè)斷面滯后響應(yīng)時(shí)間。其中，支流峰值時(shí)刻（42 h）產(chǎn)生的擾動(dòng)在監(jiān)測(cè)斷面的響應(yīng)過程如圖5所示。圖5中可看出，雖然絕對(duì)擾動(dòng)程度不同，但支流峰值時(shí)刻產(chǎn)生的擾動(dòng)都在22 min左右對(duì)監(jiān)測(cè)斷面水流產(chǎn)生了最大干擾，此時(shí)監(jiān)測(cè)斷面觀測(cè)的流量值應(yīng)包含了最多的支流入流信息。因此，對(duì)應(yīng)于峰值時(shí)刻的滯時(shí)Trj-G取為22 min。

4.1.2 EnKF校正結(jié)果 根據(jù)支流入流初步計(jì)算結(jié)果以及監(jiān)測(cè)斷面滯時(shí)矩陣，由3.3節(jié)EnKF校正方法對(duì)支流入流初值進(jìn)行校正，校正結(jié)果見圖6。

結(jié)果表明，僅由水動(dòng)力模型初步計(jì)算的支流入流存在很大誤差，峰值和峰現(xiàn)時(shí)間都嚴(yán)重偏離實(shí)測(cè)入流過程；但經(jīng)兩次EnKF校正后，支流入流基本與實(shí)際過程相吻合。其中，第一次校正后，所得結(jié)果相對(duì)均方根誤差仍然偏大，校正的支流入流在38 h左右與實(shí)測(cè)過程出現(xiàn)了較大偏離；經(jīng)過EnKF再次校正后，偏離情況明顯得到改善，最終計(jì)算結(jié)果相對(duì)均方根誤差收斂至5%以內(nèi)?；贓nKF的支流反分析校正結(jié)果誤差指標(biāo)見表1。此外，表1和圖6中還給出了采用馬斯京根方法反算的支流入流結(jié)果。對(duì)比結(jié)果表明，僅經(jīng)一次EnKF校正后所推算的支流入流即比馬斯京根反算的入流精度更高。此外，馬斯京根反算方法對(duì)入流誤差敏感，當(dāng)誤差較大時(shí)，反算結(jié)果出現(xiàn)負(fù)值[14]。

表1 理想案例支流入流反分析結(jié)果誤差評(píng)價(jià)

4.1.3 觀測(cè)誤差影響分析 支流入流反分析過程中，正、反向水動(dòng)力演算精度受觀測(cè)資料誤差影響，從而間接影響EnKF分析結(jié)果。同時(shí)，用于EnKF校正的觀測(cè)斷面流量過程會(huì)直接影響EnKF分析結(jié)果。研究分別考慮入流邊界誤差、觀測(cè)斷面流量誤差以及兩種誤差同時(shí)存在時(shí)的反演結(jié)果，進(jìn)一步分析觀測(cè)誤差對(duì)支流反算結(jié)果的影響，計(jì)算結(jié)果見表1。流量觀測(cè)誤差通常與流量大小正相關(guān)[30]，本文按比例0.1%和1%人工生成正態(tài)分布觀測(cè)誤差（式（19））：

式中a為比例系數(shù)，分別取值0.1%和1%。

由表1可以看出，入流邊界和觀測(cè)斷面誤差都對(duì)支流反演精度有一定負(fù)面影響，其中觀測(cè)斷面誤差較入流邊界誤差影響更大。盡管如此，在不同誤差擾動(dòng)情形下，本文方法推算的支流入流過程各項(xiàng)評(píng)價(jià)指標(biāo)都很接近，本文方法對(duì)入流邊界和觀測(cè)斷面流量誤差不敏感。

4.2 實(shí)例應(yīng)用經(jīng)理想案例驗(yàn)證后，將方法應(yīng)用于西江流域潯江河段的實(shí)例中，以檢驗(yàn)其實(shí)用性。潯江河段為西江主干河道，位于廣西壯族自治區(qū)境內(nèi)。該段河道沿程分布有平南、藤縣、梧州等縣市城鎮(zhèn)，是西江流域洪水高發(fā)河段。本文主要分析大湟江口至長洲庫區(qū)入口（距離大湟江口130.74 km）河段，其中距離大湟江口73.89 km處有支流蒙江匯入主干河道，在距離大湟江口96.19 km處有監(jiān)測(cè)站點(diǎn)藤縣水文站，具體位置如圖7所示。

圖7 潯江河段、支流及監(jiān)測(cè)站點(diǎn)位置

研究分別選取大湟江口不同特征的三場(chǎng)次洪水進(jìn)行分析，其中，蒙江匯入處干流斷面Q-Z關(guān)系曲線以及三場(chǎng)次洪水對(duì)應(yīng)藤縣水文站監(jiān)測(cè)流量過程如圖8所示。

采用第3節(jié)支流反分析方法計(jì)算支流蒙江入流過程，結(jié)果見圖9和表2。

圖8 三實(shí)測(cè)場(chǎng)次洪水過程及潯江斷面Q-Z關(guān)系曲線

圖9 潯江實(shí)例中支流蒙江入流反分析計(jì)算結(jié)果

圖9和表2中結(jié)果表明，本文方法推算的三場(chǎng)次洪水其支流蒙江入流都與實(shí)測(cè)過程十分接近，支流推算精度高。綜合圖9和表2可看出，雖然由水動(dòng)力模型估算的蒙江入流過程已較為接近實(shí)測(cè)入流過程，但仍具有一定誤差；經(jīng)一次EnKF校正后，三場(chǎng)次洪水蒙江支流入流過程均得到明顯改善，計(jì)算的R2和NSE都改善至0.99以上，而相對(duì)均方根誤差縮小至校正前的25%以內(nèi)。與理想案例不同，經(jīng)EnKF再次校正后的蒙江入流過程改善并不明顯，評(píng)價(jià)指標(biāo)反而略微劣化。由此可見，只有當(dāng)首次校正的支流入流過程可靠度較低時(shí)（初始入流估算誤差較大時(shí)，其中初始入流估算誤差與河道地形、擾動(dòng)滯時(shí)有關(guān)）才需二次校正，否則經(jīng)EnKF一次校正即可達(dá)到較好效果?？傮w來說，本文提出的支流入流反分析方法對(duì)于天然河道實(shí)例也具備有效性。

表2 潯江實(shí)例中支流蒙江入流反分析結(jié)果誤差評(píng)價(jià)

5 結(jié)論

本文提出了一種基于EnKF的無實(shí)測(cè)資料區(qū)間支流反分析方法，方法首先由河道正、反向水流計(jì)算初步估算支流入流過程，然后采用EnKF對(duì)初步估算過程進(jìn)行校正。將該方法應(yīng)用于理想案例和西江實(shí)例，主要結(jié)論為：（1）研究以支流匯入斷面為界，基于一維河道水動(dòng)力模型正、反向水流計(jì)算初步估算支流入流，解決了無初始入流資料的問題。構(gòu)建了監(jiān)測(cè)斷面滯時(shí)矩陣，計(jì)算水流擾動(dòng)傳播時(shí)間，以確定用于支流校正的監(jiān)測(cè)斷面觀測(cè)值，解決了EnKF滯時(shí)校正的問題。（2）采用EnKF二次校正有利于改善支流分析精度，但二次校正并非必要。當(dāng)初始估算的支流入流誤差較大、首次校正的支流入流過程可靠度較低時(shí)，運(yùn)用EnKF二次校正可顯著改善首次校正結(jié)果精度。（3）將方法應(yīng)用于理想案例和西江實(shí)例，推算的支流入流過程與實(shí)測(cè)過程十分接近，支流入流R2和NSE皆在0.99以上，相對(duì)均方根誤差也小于0.05。結(jié)果表明，本文提出的方法對(duì)于計(jì)算無實(shí)測(cè)資料的區(qū)間支流有效可行。

本文提出的區(qū)間支流反分析方法對(duì)于其他相關(guān)問題也有借鑒意義。例如，可將EnKF應(yīng)用于沿程均勻旁側(cè)入流校正分析、河道邊界條件校正、區(qū)間取水反分析、污染源溯源等問題中。另外，本文主要討論了單條區(qū)間支流的情況，實(shí)際河道往往存在多條缺乏資料的區(qū)間支流，未來將研究如何對(duì)多條區(qū)間支流進(jìn)行反分析。