基于EnKF的無實測資料區(qū)間支流反分析

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

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

1 研究背景

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

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

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

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關(guān)系曲線作為分段河道的邊界條件，具體步驟如下：

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

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

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

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

（3）以監(jiān)測斷面G處實測流量過程QG(t)為上邊界條件，步驟（2）計算的斷面j處水位過程Zj(t)（反向水流計算對邊界條件敏感，若將步驟（1）得到的斷面j處Q-Z關(guān)系曲線作為邊界條件，反向結(jié)果不穩(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處觀測到，存在滯時響應(yīng)的問題[27]。由數(shù)據(jù)同化方法對支流入流進行校正時，采用的觀測斷面G處實測值應(yīng)盡可能地包含當前時刻支流誤差擾動信息。用于校正的監(jiān)測值QG(tG)與支流入流值QL(tj)存在如下時間滯后關(guān)系：

式中，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ù)對模型預(yù)測值進行動態(tài)校正。根據(jù)3.1節(jié)支流初步計算結(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為模型預(yù)測誤差協(xié)方差；負上標代表校正前狀態(tài)變量。

（4）支流入流校正。

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

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

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

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

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

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

4 模型應(yīng)用

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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