有限差分卡爾曼濾波算法估計(jì)污染物在地下水中遷移轉(zhuǎn)化
——在遼河流域的應(yīng)用

王新民，陳慧穎，馬旭東

(長春工業(yè)大學(xué) 基礎(chǔ)科學(xué)學(xué)院，吉林 長春 130012)

?

有限差分卡爾曼濾波算法估計(jì)污染物在地下水中遷移轉(zhuǎn)化
——在遼河流域的應(yīng)用

王新民，陳慧穎*，馬旭東

(長春工業(yè)大學(xué) 基礎(chǔ)科學(xué)學(xué)院，吉林 長春 130012)

采用有限差分方法,用差商代替微商導(dǎo)出卡爾曼濾波的轉(zhuǎn)移矩陣，從而估計(jì)其主要污染物氨氮的遷移轉(zhuǎn)化，為遼河流域治理提供技術(shù)支持。

地下水污染；污染物遷移轉(zhuǎn)化；有限差分；卡爾曼濾波

0 引 言

在過去的幾十年里，地表水和地下水污染問題已成為地質(zhì)科學(xué)研究的主要領(lǐng)域,因?yàn)樗鼘θ祟惿鐣?huì)經(jīng)濟(jì)生活有著重要影響。包括污染物的跟蹤、地下水水質(zhì)的預(yù)測。目前已經(jīng)建立了很多針對不同情況的地下水污染數(shù)值模型，在大多數(shù)情況下，這些模型許多都是帶有參數(shù)的偏微分方程，模型是否可靠很多情況下取決于水文地質(zhì)參數(shù)以及邊界條件的概化是否準(zhǔn)確。為此，許多學(xué)者們提出了大量關(guān)于水文地質(zhì)參數(shù)估計(jì)的方法，用以提高模型的預(yù)測能力，如蒙特卡羅法、正則化方法等。這些方法在某種程度上可以處理模型參數(shù)的不確定性，但忽略了構(gòu)建模型時(shí)模型本身存在的誤差以及獲得數(shù)據(jù)時(shí)的測量誤差。參數(shù)估計(jì)的另一種選擇是連續(xù)的數(shù)據(jù)同化，卡爾曼濾波器(KF)是最佳的順序同化方案之一，它可以用來濾除觀測數(shù)據(jù)中的噪聲和干擾，以達(dá)到最佳的狀態(tài)估計(jì)[1-2]。當(dāng)模型和測量噪聲在非線性模型情況下，可使用擴(kuò)展卡爾曼濾波器[3](EKF)將模型線性化。并且可將狀態(tài)變量(即水體中待求的污染物濃度)和相關(guān)參數(shù)合并為一個(gè)向量進(jìn)行計(jì)算，這樣可將模型參數(shù)和污染物濃度同時(shí)確定。但由于實(shí)際應(yīng)用中擴(kuò)展卡爾曼濾波線性化過程的不穩(wěn)定且Jacobian矩陣的推導(dǎo)十分復(fù)雜，文中采用有限差分的方法代替擴(kuò)展卡爾曼濾波中的Taylor展開，導(dǎo)出卡爾曼濾波的轉(zhuǎn)移矩陣。

1 水污染模型的差分格式

有限差分法[4]的基本思想是用差商近似代替方程中的微商，將連續(xù)的問題離散化，然后耦合初始條件及邊界條件，求解封閉的線性代數(shù)方程組。這里只需完成第一部分。

考慮如下地下水污染模型:

式中：C——污染物濃度;

Dx，Dy——彌散系數(shù);

u——流速。

2 卡爾曼濾波技術(shù)在地下水污染模型中的應(yīng)用

卡爾曼濾波有兩個(gè)基本方程[5]：狀態(tài)方程和測量方程。對于地下水污染系統(tǒng)而言，狀態(tài)方程是描述地下水污染物運(yùn)移過程的方程式，傳統(tǒng)的卡爾曼濾波是針對線性系統(tǒng)的一種算法，而實(shí)際的觀測模型往往是非線性的，包括地下水污染模型。文中通過上述有限差分法將模型線性化，加上誤差隨機(jī)項(xiàng)，即可得到卡爾曼濾波的系統(tǒng)狀態(tài)方程：

系統(tǒng)的測量方程為：

式中：Z(i,j,n)——n時(shí)刻節(jié)點(diǎn)觀測的矢量值；

H(n)——測量矩陣，由監(jiān)測孔相對于模型分節(jié)點(diǎn)的位置確定；

W(n)，V(n)——誤差項(xiàng)，一般假設(shè)相互獨(dú)立，服從高斯分布，并具有下列協(xié)方差矩陣

卡爾曼濾波算法主要包括兩步[6]：預(yù)測和更新。在預(yù)測階段，用運(yùn)移模型根據(jù)t-1時(shí)刻狀態(tài)向量和誤差協(xié)方差矩陣預(yù)測t時(shí)刻的狀態(tài)和協(xié)方差矩陣。

式中：“∧”——估計(jì)值；

“—”——前面式子算出來的估計(jì)值。

更新階段是對模型狀態(tài)和參數(shù)進(jìn)行時(shí)間更新。更新過程可以表示為：

為卡爾曼增益，類似一個(gè)信心指數(shù)，用來說明對測量值還是預(yù)測值的信心程度,I是單位矩陣。經(jīng)過足夠多步的計(jì)算，模型狀態(tài)、參數(shù)的估計(jì)值就會(huì)逐漸逼近真實(shí)值，從而達(dá)到狀態(tài)參數(shù)估計(jì)的目的[7]。

3 實(shí)例分析

研究區(qū)位于沈陽市區(qū)西南部，渾河北岸漫灘地區(qū)，地理坐標(biāo)為北緯41°30′～42°00′，東經(jīng)123°00′～123°40′，面積約37.88 km2。研究區(qū)內(nèi)地表水體主要有渾河及其支流細(xì)河，屬于遼河水系。地下水水質(zhì)以“三氮”(氨氮、硝酸根離子、亞硝酸根離子)污染為主。其中氨氮主要來源于人和動(dòng)物的排泄物、降水的垂直入滲、工業(yè)廢水以及河流的側(cè)向補(bǔ)給。主要排泄方式為人工開采。研究區(qū)地下徑流方向?yàn)橛蓶|北向西南，東南部受渾河影響，西北部臨界細(xì)河，可概化為一類定濃度邊界；北部和西部是地下水側(cè)向徑流補(bǔ)給和排泄段，可概化為二類定通量邊界[8-9]。依據(jù)上述特征建立數(shù)學(xué)模型：

式中：R——阻滯系數(shù);

c——溶質(zhì)氨氮的濃度;

D——彌散系數(shù);

V——流速;

λ——衰變常數(shù)，氨氮在運(yùn)移過程中會(huì)不斷衰變，發(fā)生硝化反應(yīng);

Ω——整個(gè)模型區(qū)域;

Γ1——定濃度邊界;

Γ2——通量邊界;

ck(x,y,t)——已知的初始濃度分布;

g(x,y,t)——已知函數(shù)，表示對流彌散通量。

研究區(qū)采用100×100等間距有限差分法離散，共剖分矩形單元格17 402個(gè)。研究區(qū)網(wǎng)格剖分圖如圖1所示。

圖1 研究區(qū)網(wǎng)格剖分圖

運(yùn)用前述的卡爾曼濾波最優(yōu)估計(jì)算法，選取已知監(jiān)測井45個(gè)，統(tǒng)一測定時(shí)間為2011年3月至2012年9月，測量誤差的方差取0.002 5。

校正后的彌散系數(shù)和噪聲標(biāo)準(zhǔn)差見表1。

表1 校正前后的模型參數(shù)表

研究區(qū)2011年3月實(shí)測初始氨根濃度等值線圖和估計(jì)協(xié)方差等值線圖分別如圖2和圖3所示。

圖2 研究區(qū)初始氨根濃度等值線圖

圖3 研究區(qū)初始濃度估計(jì)協(xié)方差等值線圖

研究區(qū)2012年9月實(shí)測氨根濃度等值線圖和估計(jì)協(xié)方差等值線圖分別如圖4和圖5所示。

圖4 實(shí)測氨根濃度等值線圖

圖5 實(shí)測氨根濃度估計(jì)協(xié)方差等值線圖

研究區(qū)氨根濃度最優(yōu)估計(jì)等值線圖及估計(jì)誤差協(xié)方差等值線圖分別如圖6和圖7所示。

圖6 研究區(qū)氨根濃度最優(yōu)估計(jì)等值線圖

圖7 研究區(qū)氨根濃度估計(jì)協(xié)方差等值線圖

比較圖4和圖6，以及5和圖7可以看出，通過有限差分卡爾曼濾波算法計(jì)算出的濃度場與實(shí)測濃度場基本吻合，可以反映研究區(qū)氨根離子的空間分布。

研究區(qū)部分采樣點(diǎn)氨根濃度實(shí)測值與最優(yōu)估計(jì)值比較如圖8所示。

圖8 研究區(qū)部分采樣點(diǎn)氨根濃度實(shí)測值與最優(yōu)估計(jì)值比較

從圖8中可以看出，45個(gè)采樣點(diǎn)中除2、4、11、12、40五個(gè)采樣點(diǎn)擬合較差外,其余40個(gè)采樣點(diǎn)的最優(yōu)估計(jì)值與實(shí)際測量值擬合的比較好，可能由于觀測精度或誤差協(xié)方差過大引起的，但能表現(xiàn)出算法的可靠性。通過文中提出的有限差分卡爾曼濾波算法校正后的模型，基本能反映實(shí)際氨根離子在研究區(qū)的遷移轉(zhuǎn)化特征，最優(yōu)估計(jì)值與實(shí)際測量值擬合良好。

由圖中可知，通過文中提出的有限差分卡爾曼濾波算法校正后的模型，基本能反映實(shí)際氨根離子在研究區(qū)的遷移轉(zhuǎn)化特征，最優(yōu)估計(jì)值與實(shí)際測量值擬合較好。

4 結(jié) 語

綜合分析卡爾曼濾波算法的實(shí)際功能和地下水污染物運(yùn)移的實(shí)際情況，文中提出一種基于有限差分改進(jìn)的卡爾曼濾波算法聯(lián)合估計(jì)地下水污染物運(yùn)移模型的狀態(tài)及相關(guān)參數(shù)的新方法，通過對遼河流域地下水污染物遷移轉(zhuǎn)換問題的分析，驗(yàn)證了該方法在地下水污染物運(yùn)移模型上的可行性，使監(jiān)測數(shù)據(jù)得到充分利用，并且充分考慮到各種誤差干擾，有效提高模型準(zhǔn)確性。

[1] Hendricks Franssen H J,Kinzelbach W,Hendricks Franssen，et al. Real-time groundwater flow modeling with the ensemble Kalman Filter:Joint estimation of states and parameters and the filter inbreeding problem [J]. Water Resources Research,2008,44(9):108-128.

[2] 王昕.卡爾曼濾波在模型參數(shù)估計(jì)中的應(yīng)用[D].哈爾濱：哈爾濱工業(yè)大學(xué),2008.

[3] Assumaning G A，Chang S Y. Use of simulation filters in three-dimensional groundwater contaminant transport modeling [J]. Journal of Environmental Engineering，2012,138(11):1122-1129.

[4] 王新民,董小剛.計(jì)算方法簡明教程[M].北京：科學(xué)出版社，2010.

[5] 孫永泉,劉梅俠,董學(xué)良,等.運(yùn)用Kalman濾波技術(shù)校正地下水流系統(tǒng)確定-隨機(jī)性數(shù)值模型[J].黑龍江水專學(xué)報(bào),2010,37(1):110-114.

[6] 王廣德,余蕓生,賀宜達(dá)日.地下水可溶性污染物遷移的最優(yōu)估算[J].環(huán)境科學(xué)學(xué)報(bào)，1988,8(4):385-392.

[7] Gharamti M E,Hoteit I. Complex step-based low-rank extended Kalman filtering for state-parameter estimation in subsurface transport models [J]. Journal of Hydrology,2014,509(4):588-600.

[8] 李波，王建帥，任媛，等.傍河型水源地水質(zhì)預(yù)測解析解模型[J].長春工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2012,33(3):245-249.

[9] 許可,陳鴻漢.地下水中三氮污染物遷移轉(zhuǎn)化規(guī)律研究進(jìn)展[J].中國人口：資源與環(huán)境,2011(S2):421-424.

Study on the estimation of pollutant transfer and transformation of pollutants in groundwater with finite difference Kalman filter——Application of the Liaohe River basin

WANG Xinmin,CHEN Huiying*,MA Xudong

(School of Basic Science,Changchun University of Technology,Changchun 130012，China)

Derivative is replaced with difference quotient in Finite difference method to deduce the transfer matrix of the Kalman filter for estimating the main pollutants,ammonia nitrogen migration and transformation. It provides some technical support for the Liaohe River basin management.

groundwater pollution; pollutant transfer; finite difference; Kalman filter.

2016-04-20

國家自然科學(xué)基金資助項(xiàng)目(51278065)

王新民(1957-)，男，漢族，山東牟平人，長春工業(yè)大學(xué)教授，博士，主要從事水環(huán)境系統(tǒng)正反問題數(shù)學(xué)模型與數(shù)值方法研究,E-mail:wxm@jlu.edu.cn. *通訊作者：陳慧穎(1991-)，女，漢族，黑龍江雙鴨山人，長春工業(yè)大學(xué)碩士研究生，主要從事水環(huán)境系統(tǒng)正反問題數(shù)值模擬方向研究,E-mail:8497871436@qq.com.

O 242

A

1674-1374(2016)05-0417-05