顧爐華，賴錫軍

(1.江蘇省環(huán)境科學研究院，江蘇 南京 210036; 2.長江勘測規(guī)劃設計研究有限責任公司上海分公司，上海 200439; 3.中國科學院南京地理與湖泊研究所，江蘇 南京 210008)

為減少數(shù)值模型的不確定性，數(shù)據(jù)同化方法被引入到水文、水動力模型[1-2]。數(shù)據(jù)同化方法可分為確定性方法和不確定性方法，前者由果及因，根據(jù)已知結果反推未知參數(shù)，代表方法有變分法[3-4]；后者基于統(tǒng)計估值理論，通過實時融入觀測實現(xiàn)預報場的動態(tài)更新，代表方法有卡爾曼濾波和粒子濾波算法[5-6]?？柭鼮V波及其改進算法以其在高維度、非線性系統(tǒng)中的良好表現(xiàn)在水文預報領域得到應用與發(fā)展，成為實時校正的主流技術之一。但這些算法都假設噪聲在空間上為零均值高斯分布，在非高斯系統(tǒng)中應用會出現(xiàn)濾波發(fā)散問題。而粒子濾波算法能很好地解決這一問題，粒子濾波又叫順序蒙特卡羅濾波，是對貝葉斯估值理論的一種蒙特卡羅算法實現(xiàn)。與經(jīng)驗方法、優(yōu)化插值、連續(xù)方法等傳統(tǒng)實時校正方法及集合卡爾曼濾波方法相比，粒子濾波算法的優(yōu)勢主要體現(xiàn)為：不受高斯分布假設的約束、能更好表現(xiàn)非線性系統(tǒng)的變化信息、計算效率高、算法容易實現(xiàn)等[7-10]。

關于粒子濾波方法應用于水力學模型的研究目前主要集中在濾波參數(shù)選擇、算法結構完善、洪水預報等方面。Moradkhani等[11]運用粒子濾波對模型狀態(tài)變量及參數(shù)展開了不確定性分析；Nohs等[12]提出了基于粒子濾波的狀態(tài)變量和參數(shù)同步更新算法；畢海蕓等[13-16]提出了殘差重采樣、分層重采樣、聚合重采樣等方法，進一步完善了粒子濾波的性能；楊瑞祥等[17]和徐興亞等[18]分別將粒子濾波應用于水文模型和水動力模型進行洪水預報，對濾波的實用性進行了進一步驗證。這些成果豐富了粒子濾波方法的理論和應用。

本文將粒子濾波方法引入一維非恒定流模擬計算中，建立一維水量模擬的多變量校正方法，分析粒子數(shù)量、狀態(tài)變量擾動誤差標準差和邊界條件的不確定性對濾波性能的影響，比較粒子濾波與集合卡爾曼濾波的性能，為粒子濾波方法在河網(wǎng)水量分析、洪水預報等實際應用提供技術支撐。

1 粒子濾波基本原理

1.1 貝葉斯濾波原理

貝葉斯濾波利用觀測數(shù)據(jù)來構造系統(tǒng)狀態(tài)的后驗概率密度，利用狀態(tài)轉(zhuǎn)移模型估計系統(tǒng)狀態(tài)的先驗概率密度，再結合最新的觀測數(shù)據(jù)對先驗概率密度進行修正，繼而得到系統(tǒng)狀態(tài)的后驗概率分布。

首先，建立狀態(tài)方程和量測方程：

Xk=AXk-1+BUk+Wk

(1)

Zk=Hk+Vk

(2)

式中：Xk為k時刻系統(tǒng)的狀態(tài)量；Uk為k時刻系統(tǒng)的控制量；A、B為狀態(tài)系統(tǒng)參數(shù)矩陣；Zk為k時刻的觀測量；Hk為k時刻量測系統(tǒng)的觀測矩陣；Wk、Vk分別為k時刻狀態(tài)系統(tǒng)和量測系統(tǒng)的噪聲。

假設狀態(tài)系統(tǒng)Xk服從一階馬爾科夫過程且量測系統(tǒng)Zk獨立，那么根據(jù)條件概率理論，式(1)可表達為先驗概率p(Xk|Zk-1)，式(2)可表達為似然概率p(Zk|Xk)。由貝葉斯條件概率公式，可推出后驗概率p(Xk|Zk)的表達式：

(3)

在實際應用中，大多數(shù)的后驗概率密度都是很難直接進行解析的，為此引入多重積分的蒙特卡羅方法。

1.2 蒙特卡羅方法

蒙特卡羅方法的核心思想是把待求解的問題抽象為隨機事件，通過從已知的概率分布中隨機抽樣，實現(xiàn)對所求問題的解的近似估計。蒙特卡羅方法的理論基礎是大數(shù)問題，當樣本足夠大時，蒙特卡羅估計就等同于貝葉斯后驗概率密度。

從p(Xk|Zk)中隨機抽取獨立分布的N個樣本{Xk1,Xk2,…,XkN}，狀態(tài)系統(tǒng)的后驗概率密度分布可近似表示為

(4)

1.3 粒子濾波方法

粒子濾波方法是用蒙特卡羅方法求解貝葉斯后驗概率密度的方法，它將每一個樣本視作一個粒子，以似然值作為粒子的權重，通過對粒子加權平均得到后驗概率密度。粒子濾波可分為4步：模型預報、重要性采樣、重采樣和粒子擾動。

模型預報是在生成若干個水流狀態(tài)粒子后，基于數(shù)學模型驅(qū)動各個粒子平行計算，得到每一個時刻水流狀態(tài)的預報值。

粒子濾波基于序貫重要性采樣，即在蒙特卡羅思想中應用統(tǒng)計學理論中的序貫分析方法，以遞推得到后驗概率密度函數(shù)的最優(yōu)估計。其實現(xiàn)方式是計算每個粒子的似然值，并將其歸一化。本文采用如下似然函數(shù)：

(5)

式中：wki為k時刻第i個樣本的似然值；σ為相對標準差。

在粒子權重更新的過程中，不可避免地會出現(xiàn)“粒子退化”現(xiàn)象，造成大部分計算資源都用在權重很小的粒子上，因此需要用到重采樣技術，其本質(zhì)是復制大權重粒子，剔除小權重粒子，本文采用多項式重采樣，其原理見文獻[19]。

重采樣之后，需再次對粒子群施加擾動，保證粒子群在大權重粒子占多數(shù)的情況下仍有一定的多樣性。本文采用如下擾動方法:

X(i)=X(i)(1+N(i)σ)

(6)

式中：X(i)為第i個河道斷面的狀態(tài)量；N(i)為第i個河道斷面的偽光滑隨機擾動場，服從均值為0、方差為1的正態(tài)分布。

2 基于粒子濾波的一維非恒定流水量校正模型

2.1 一維非恒定流模型

本文建立的一維非恒定流模型基于圣維南方程組：

(7)

(8)

式中：Z為水位，m；Q為流量，m3/s；K為流量模數(shù)；m3/s；q為單位河長旁側入流，m2/s；A為過水斷面面積，m2；α為動量校正系數(shù)；g為重力加速度，m/s2；x為沿水流方向距離，m；t為時間，s；B為水面寬度，m。

采用Preissmann四點加權隱格式離散上述兩式得到：

A1jΔQj+B1jΔZj+C1jΔQj+1+D1jΔZj+1=E1j(9)

A2jΔQj+B2jΔZj+C2jΔQj+1+D2jΔZj+1=E2j(10)

式中符號的含義可見文獻[20]。對于緩流，在已知上下游邊界條件的情況下，可用追趕法求出各系數(shù)，進而求得河道各斷面每一時刻的水位值和流量值。

2.2 水量校正模型算法流程

首先對水位和流量施加隨機擾動，生成N個水流狀態(tài)粒子；沒有觀測數(shù)據(jù)融入的時刻，驅(qū)動粒子進行并行計算，形成粒子群搜索空間；在有觀測數(shù)據(jù)輸入的時刻，啟動粒子濾波模塊，對觀測數(shù)據(jù)的類型及誤差進行分析；計算每個粒子與觀測之間的似然值，得到粒子權重；接著對粒子群進行多項式重采樣，計算粒子的平均值并輸出。圖1為第k時段內(nèi)粒子濾波水量校正模型算法流程。

圖1 基于粒子濾波的水量校正模型算法流程(第k時段內(nèi))

3 合成數(shù)值試驗

3.1 案例選擇

為分析粒子濾波在非恒定流計算中的應用性能，以單一河道洪水演進過程為例開展合成數(shù)值試驗。該河道全長60 km，上游給定流量邊界，下游給定水位邊界。共設置4個觀測站點，分別布置在距離上游11 km、23 km、35 km和47 km處，觀測時間間隔為0.5 h。

3.2 試驗設計

選擇沿程斷面水位和流量組成的矩陣作為水流狀態(tài)粒子，將上游來流過程作為主要的誤差源，將其按比例縮放作為上游預報洪水過程，研究粒子濾波在不同預報誤差和粒子擾動誤差標準差下的校正性能。將模型在實際流量過程下的計算結果作為真實值，從中抽取水位、流量作為觀測數(shù)據(jù)。此外，本文計算了上游預報洪水過程下純水動力模型的計算結果作為河道預報水位，與數(shù)據(jù)同化結果進行對比，以檢驗粒子濾波對水動力模型精度的改善效果。從粒子數(shù)量、邊界條件縮放系數(shù)、水位和流量擾動誤差標準差4個方面設置計算條件，開展情景計算，各試驗相關參數(shù)見表1。

表1 各試驗參數(shù)值

3.3 試驗結果與分析

3.3.1粒子數(shù)量對粒子濾波的影響

試驗A1～A4的結果顯示，粒子數(shù)量越多，流量殘差和水位殘差越小。從圖2可以看出，當粒子數(shù)量小于100時，流量和水位殘差隨粒子數(shù)量呈快速下降趨勢；當粒子數(shù)量大于100后，流量和水位殘差減小的速度減緩。粒子數(shù)量為100時，流量和水位殘差分別為17.54 m3/s和0.083 m。更大的粒子數(shù)量雖然可以帶來更好的效果，但是也會意味著更大的計算量。在本算例中，粒子數(shù)量取100可同時滿足計算精度與效率的要求。

圖2 不同粒子數(shù)量下水位和流量殘差

3.3.2擾動誤差標準差對粒子濾波的影響

狀態(tài)變量擾動誤差決定了粒子群搜索空間大小和粒子群分布密度。搜索空間太小將限制粒子活動范圍，使其無法到達真實值附近；搜索空間太大又會造成粒子群分布密度降低，真實值周邊缺乏有效粒子。采用純水動力模型的計算結果(預報值)獲得的水位殘差為0.160 m，流量殘差為35.45 m3/s，水位精度為0.950，流量精度為0.942。測試結果顯示，試驗B3中參數(shù)設置對發(fā)揮粒子濾波效果最有利，水位殘差和流量殘差分別為0.078 m和17.37 m3/s，較預報值情況下降低了約50%，水位和流量模擬精度達到99%(表2)。

表2 不同擾動誤差標準差下水位和流量殘差統(tǒng)計

3.3.3邊界條件對粒子濾波的影響

由表3可知，當邊界條件縮放系數(shù)取0.9和1.1時，水位殘差分別只有0.08 m和0.07 m，對應的流量殘差均小于20 m3/s；當邊界條件縮放系數(shù)取0.7和1.3時，水位殘差均達到0.26 m，流量殘差分別高達54.54 m3/s和64.57 m3/s。平均納什效率系數(shù)(NSE)也顯示，當邊界條件相對誤差在20%以內(nèi)，即邊界條件縮放系數(shù)在0.8～1.2之間時，粒子濾波效果較好，模擬精度可達94%以上；當邊界條件相對誤差大于20%時，濾波后精度小于90%(表3)。

表3 不同縮放系數(shù)下水位和流量殘差統(tǒng)計

3.3.4粒子濾波與集合卡爾曼濾波的比較

將粒子濾波與文獻[6]中基于集合卡爾曼濾波的多變量交替校正方法進行對比，為保證結果的合理性，模型邊界條件、狀態(tài)變量擾動誤差標準差、集合與粒子個數(shù)等主要參數(shù)均保持一致。選取距上游邊界10 km(斷面1)、30 km(斷面2)、50 km(斷面3)3個斷面流量過程作為比較對象，總體而言，粒子濾波與集合卡爾曼濾波都具有很好的濾波效果，3個斷面的平均流量預報精度分別達到了98.4%和98.8%，集合卡爾曼濾波算法在結果的準確性和穩(wěn)定性上更具優(yōu)勢，洪峰時刻采用集合卡爾曼濾波分析出的洪峰流量比粒子濾波更接近真實值；洪水后期，采用粒子濾波分析出的流量圍繞真實值上下波動，而采用集合卡爾曼濾波獲得的流量與真實值緊密貼合(圖3)。粒子濾波由于無需進行復雜的矩陣計算，從而具有更高的計算效率。在本文使用的案例中，使用集合卡爾曼濾波算法完成一次洪水過程的水量校正需要的CPU時間為883 s，而粒子濾波只要379 s，計算效率約是前者的2.3倍。

圖3 不同河道斷面流量過程比較

4 實例驗證

4.1 案例選擇

選取太湖流域河網(wǎng)水量模型對粒子濾波的實際應用性能進行驗證，驗證時段為2012年全年。以河網(wǎng)節(jié)點水位作為狀態(tài)變量，以太湖(西山)、常州、金壇、溧陽等13個站點的逐日實測水位作為更新模型狀態(tài)的觀測資料，水位擾動誤差標準差取為實際值的5%，考慮到模型計算量過大，為減少計算時間，本案例中將粒子數(shù)取為50。這種選擇雖會在一定程度上減弱濾波效果，但是不妨礙驗證粒子濾波算法在太湖河網(wǎng)水量同化中應用的有效性。太湖流域河網(wǎng)及水文站分布見圖4。

圖4 太湖流域河網(wǎng)及水文站分布

4.2 計算結果與分析

平原河網(wǎng)地區(qū)河道多為往復流，流量站少于水位站，且大型水利樞紐多以水位作為調(diào)度控制指標，水位對流域水資源的管理與調(diào)度更具意義，因此本次分析僅關注流域水位過程。

河網(wǎng)概化、邊界條件、閘站調(diào)度等方面的不確定性使得太湖流域河網(wǎng)水量模型模擬出的太湖(西山)、溧陽、蘇州、新市4站的水位與實測值存在較大差距。濾波后的模擬效果有了很大改善，但是不同站點的改善效果存在差異(圖5)。太湖(西山)站、溧陽站和蘇州站濾波效果較好，濾波后的水位與實測水位十分吻合，水位殘差從濾波前的0.17 m、0.40 m、0.10 m分別降為0.05 m、0.10 m和0.05 m，納什效率系數(shù)也分別達到了91%、82%和84%；新市站的濾波效果稍差，濾波后的水位仍低于實測水位，但相比于河網(wǎng)水量模型預報結果，水位殘差從0.29 m降至0.18 m，降幅為38%?？傮w而言，粒子濾波應用于實際河網(wǎng)水量校正可有效改善模擬精度。

圖5 主要水文站水位過程對比

5 結 語

本文建立了基于粒子濾波的一維非恒定流水量校正模型，通過對水流狀態(tài)粒子進行統(tǒng)計估值，實現(xiàn)模型狀態(tài)變量的最優(yōu)估計。探討了多變量校正方法在一維非恒定流計算中的性能，并分析了粒子數(shù)量、狀態(tài)變量擾動誤差和邊界條件的不確定性對濾波效果的影響。結果證實了粒子濾波可顯著改善模擬精度，且具有較高的計算效率。