顧爐華，賴(lài)錫軍

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

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

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

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

1 粒子濾波基本原理

1.1 貝葉斯濾波原理

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

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

Xk=AXk-1+BUk+Wk

(1)

Zk=Hk+Vk

(2)

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

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

(3)

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

1.2 蒙特卡羅方法

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

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

(4)

1.3 粒子濾波方法

粒子濾波方法是用蒙特卡羅方法求解貝葉斯后驗(yàn)概率密度的方法，它將每一個(gè)樣本視作一個(gè)粒子，以似然值作為粒子的權(quán)重，通過(guò)對(duì)粒子加權(quán)平均得到后驗(yàn)概率密度。粒子濾波可分為4步：模型預(yù)報(bào)、重要性采樣、重采樣和粒子擾動(dòng)。

模型預(yù)報(bào)是在生成若干個(gè)水流狀態(tài)粒子后，基于數(shù)學(xué)模型驅(qū)動(dòng)各個(gè)粒子平行計(jì)算，得到每一個(gè)時(shí)刻水流狀態(tài)的預(yù)報(bào)值。

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

(5)

式中：wki為k時(shí)刻第i個(gè)樣本的似然值；σ為相對(duì)標(biāo)準(zhǔn)差。

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

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

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

(6)

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

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

2.1 一維非恒定流模型

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

(7)

(8)

式中：Z為水位，m；Q為流量，m3/s；K為流量模數(shù)；m3/s；q為單位河長(zhǎng)旁側(cè)入流，m2/s；A為過(guò)水?dāng)嗝婷娣e，m2；α為動(dòng)量校正系數(shù)；g為重力加速度，m/s2；x為沿水流方向距離，m；t為時(shí)間，s；B為水面寬度，m。

采用Preissmann四點(diǎn)加權(quán)隱格式離散上述兩式得到：

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

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

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

2.2 水量校正模型算法流程

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

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

3 合成數(shù)值試驗(yàn)

3.1 案例選擇

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

3.2 試驗(yàn)設(shè)計(jì)

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

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

3.3 試驗(yàn)結(jié)果與分析

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

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

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

3.3.2擾動(dòng)誤差標(biāo)準(zhǔn)差對(duì)粒子濾波的影響

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

表2 不同擾動(dòng)誤差標(biāo)準(zhǔn)差下水位和流量殘差統(tǒng)計(jì)

3.3.3邊界條件對(duì)粒子濾波的影響

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

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

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

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

圖3 不同河道斷面流量過(guò)程比較

4 實(shí)例驗(yàn)證

4.1 案例選擇

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

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

4.2 計(jì)算結(jié)果與分析

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

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

圖5 主要水文站水位過(guò)程對(duì)比

5 結(jié) 語(yǔ)

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