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

（1.清華大學，水沙科學與水利水電工程國家重點實驗室，北京 100084；2.中國水利水電科學研究院，北京 100038）

1 研究背景

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

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

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

2 水流計算

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

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

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

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

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

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

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

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

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

3 支流入流計算模型

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

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

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

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

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

當j斷面Q-Z關系曲線已知時，此步驟可省略。

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

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

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

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

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

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

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

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

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

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

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

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

（4）支流入流校正。

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

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

（6）轉入下一時步t+1時刻。

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

3.4 評價指標本文采用相關系數(shù)R2，相對均方根誤差和Nash-Sutcliffe效率系數(shù)（NSE）三項誤差指標評價模型模擬效果，計算公式如下：

式中：Mt為模擬值；Ot為實測值。

4 模型應用

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

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

圖5 支流峰值時刻擾動在監(jiān)測斷面處的響應過程

圖6 理想案例支流反分析計算結果

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

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

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

表1 理想案例支流入流反分析結果誤差評價

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

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

由表1可以看出，入流邊界和觀測斷面誤差都對支流反演精度有一定負面影響，其中觀測斷面誤差較入流邊界誤差影響更大。盡管如此，在不同誤差擾動情形下，本文方法推算的支流入流過程各項評價指標都很接近，本文方法對入流邊界和觀測斷面流量誤差不敏感。

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

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

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

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

圖8 三實測場次洪水過程及潯江斷面Q-Z關系曲線

圖9 潯江實例中支流蒙江入流反分析計算結果

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

表2 潯江實例中支流蒙江入流反分析結果誤差評價

5 結論

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

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