基于Bayesian－MCMC方法的水體污染識(shí)別反問題＊

陳海洋，滕彥國，王金生，宋柳霆，周振瑤

（北京師范大學(xué)水科學(xué)研究院，北京 100875）

隨著我國經(jīng)濟(jì)的快速發(fā)展，各種自然災(zāi)害和人們的生產(chǎn)經(jīng)營活動(dòng)帶來的環(huán)境風(fēng)險(xiǎn)不斷加劇，突發(fā)性水體污染事件時(shí)有發(fā)生，給人民群眾的生命健康、財(cái)產(chǎn)安全和生態(tài)環(huán)境造成了巨大損失.迅速找出污染源位置、掌握其排放強(qiáng)度、了解事件發(fā)生的初始時(shí)間，對于開展污染應(yīng)急處置具有重要的現(xiàn)實(shí)意義.然而事件發(fā)生時(shí)，這些信息往往是不確知的，也是不適定的.為此，需根據(jù)水體環(huán)境系統(tǒng)的自身特點(diǎn)，構(gòu)建污染擴(kuò)散對流控制方程，利用已有水文水質(zhì)參數(shù)以及濃度場的部分信息，推求污染源位置、源強(qiáng)及事件發(fā)生初始時(shí)間，即為環(huán)境水力學(xué)反問題.

與已知水力水質(zhì)參數(shù)推算污染物濃度的正問題不同，反問題通常會(huì)由于測量數(shù)據(jù)有限且?guī)в性肼?，使得污染源識(shí)別成為不適定問題［1］.常見的反演方法有正則化方法［2］、優(yōu)化算法［3］、數(shù)值法［4］和統(tǒng)計(jì)反演法［5］等.其中建立在貝葉斯統(tǒng)計(jì)學(xué)基礎(chǔ)上的貝葉斯推理能以概率語言來描述和解決工程反問題，可以較好地用來實(shí)現(xiàn)污染源識(shí)別反問題的求解［6］.

當(dāng)前，環(huán)境水力學(xué)反問題研究尚處于起步階段，研究成果不多，且以參數(shù)反演居多，部分源項(xiàng)識(shí)別研究以數(shù)學(xué)方法證明為主，對模型邊界條件、初始條件假設(shè)過于嚴(yán)格［7］.本文基于貝葉斯推理和二維水質(zhì)擴(kuò)散對流模型，考慮有限寬度河流的瞬時(shí)岸邊污染泄漏，建立水體污染識(shí)別數(shù)學(xué)模型，以典型Metropolis算法構(gòu)建馬爾科夫鏈，取得污染泄漏源源強(qiáng)、泄漏位置和污染泄漏初始時(shí)間3個(gè)模型參數(shù)后驗(yàn)概率分布，并通過基于混合遺傳模式搜索算法的求解結(jié)果進(jìn)行比較，分析MCMC－Bayes方法的識(shí)別效果.

1 正問題

一個(gè)考慮二維水質(zhì)對流擴(kuò)散方程可以表述為：

式中：C為（x，y，t）處的污染物質(zhì)量濃度，mg／L；ux和uy分別為河流縱向和橫向的流速分量，m／s；Dx和Dy分別為河流縱向和橫向的混合系數(shù)，m2／s；K為污染物降解系數(shù)，d－1.

對于瞬時(shí)點(diǎn)源排放，方程解析解可表示為：

式中：h為河流深度；M為污染物排放量，g／s；B為點(diǎn)源到邊界的距離，m；其他同上.如在岸邊排放，則B＝0，有：

2 基于MCMC－Bayes方法的污染識(shí)別反問題

2.1 貝葉斯推理

貝葉斯推理基于如下貝葉斯公式：

式中：p（X｜y）為參數(shù)的后驗(yàn)概率密度函數(shù)；X為模型參數(shù)；y為觀測數(shù)據(jù)；p（X）為參數(shù)的先驗(yàn)概率密度函數(shù)；p（y｜X）為條件概率密度；p（y）為積分常數(shù).

通常，環(huán)境水力學(xué)的模型參數(shù)都分布在一個(gè)特定的范圍內(nèi)，滿足獨(dú)立的均勻分布，其先驗(yàn)概率密度函數(shù)可定義為：

則總的先驗(yàn)分布可表示為：

而污染物濃度場的測量噪聲一般可以認(rèn)為是白噪聲，其誤差服從均值為零、標(biāo)準(zhǔn)偏差為σ的正態(tài)分布，ε～N（0，σ2），條件概率密度定義為：

式中：n為測量數(shù)據(jù)個(gè)數(shù).

把式（6）～（7）代入式（4），并把分母表示的常數(shù)項(xiàng)正規(guī)化為1可得：

理論上，利用式（8）可以獲得任何環(huán)境模型參量的統(tǒng)計(jì)矩，但事實(shí)上，后驗(yàn)概率密度的求解算法較為復(fù)雜，難以得出明確的解析表達(dá)式.另一方面，隨未知參量維數(shù)的增加，數(shù)值積分算法的計(jì)算量將呈指數(shù)增長，實(shí)現(xiàn)復(fù)雜而難度較大.為此，需要采取其他改進(jìn)的方法以實(shí)現(xiàn)后驗(yàn)概率密度的求解.

2.2 馬爾科夫鏈蒙特卡羅法

馬爾科夫鏈蒙特卡羅法（Markov Chain Monte Carlo，簡稱MCMC）是近期發(fā)展起來、可很好地用來求解貝葉斯后驗(yàn)概率密度的統(tǒng)計(jì)方法.其基本思想是對定義在高維隨機(jī)向量空間的無明確數(shù)學(xué)表達(dá)式的概率分布p進(jìn)行抽樣，產(chǎn)生大量服從分布p的隨機(jī)向量序列，通過建立合適的抽樣算法使得該序列滿足馬爾科夫鏈性質(zhì)，以保證向量Vi＋1的產(chǎn)生僅依賴于向量Vi，而與Vi－1，Vi－2，…，V1等向量無關(guān).這樣，抽樣算法完全可以由轉(zhuǎn)移概率矩陣P來描述，矩陣元素pij表示當(dāng)前訪問向量Vj的條件下接著將要訪問Vi的條件概率［8］.

可見，MCMC方法利用Markov鏈機(jī)制探索狀態(tài)空間以生成樣本的方法能夠保證Markov鏈花更多的時(shí)間在最重要的區(qū)域，產(chǎn)生的樣本更能夠模仿目標(biāo)分布樣本.常見的構(gòu)造Markov鏈轉(zhuǎn)移概率矩陣的方法有Gibbs抽樣算法、Metropolis－Hastings算法、Metropolis算法和自適應(yīng)Metropolis算法等.本研究采用Metropolis算法.

2.3 Metropolis算法

Metropolis是Metropolis－Hastings算法的特例，其Proposal分布函數(shù)滿足對稱隨機(jī)游走，即q（X＊｜X（i））＝q（X（i）｜X＊）.從當(dāng)前位置X（i）獲得候選抽樣值X＊，Markov鏈從X（i）位置移動(dòng)到X＊的接受概率為［9］：

式中：p（X＊），p（X（i））為目標(biāo)概率密度函數(shù)，在貝葉斯推理中為條件概率密度.

2.4 基于Metropolis算法的反演思路

1）在模型參數(shù)先驗(yàn)范圍內(nèi)隨機(jī)產(chǎn)生模型參數(shù)初始點(diǎn)X（i）（i＝1）；

2）利用水質(zhì)模型計(jì)算出模型參數(shù)初始點(diǎn)對應(yīng)的污染物濃度值，并計(jì)算出該模型參數(shù)對應(yīng)的條件概率密度L＿M(jìn)ODEL；

3）根據(jù)Proposal分布在當(dāng)前參數(shù)X（i）狀態(tài)下提取新的測試參數(shù)X（＊），利用水質(zhì)模型計(jì)算出測試參數(shù)X（＊）對應(yīng)的污染物濃度及條件概率密度L＿TEST；

4）產(chǎn)生一個(gè)0～1之間的隨機(jī)數(shù)u，如果u＜min（1，L＿TEST／L＿M(jìn)ODEL），則接受該測試參數(shù)并設(shè)定為當(dāng)前模型參數(shù)，即X（i＋1）＝X（＊），否則不接受該測試參數(shù)，X（i＋1）＝X（i）；

5）重復(fù)步驟3），4）直至達(dá)到迭代次數(shù).

3 實(shí)例與討論

模擬案例假定如下：某河段斷面A上游發(fā)現(xiàn)污染帶前鋒，經(jīng)檢測污染物為B危險(xiǎn)化學(xué)品，現(xiàn)需利用斷面A帶噪聲的濃度場部分分布數(shù)據(jù)反演上游污染泄漏源位置、初始泄漏時(shí)間及泄漏強(qiáng)度.假定該河流河段均勻，河流水文條件穩(wěn)定，河段斷面積、平均流速穩(wěn)定，滿足二維水質(zhì)模型應(yīng)用條件，污染排放為岸邊瞬時(shí)泄漏，相關(guān)參數(shù)及先驗(yàn)信息見表1；A斷面污染物濃度監(jiān)測數(shù)據(jù)見圖1；初始監(jiān)測時(shí)間為發(fā)現(xiàn)污染帶后的第一時(shí)間，某日下午15時(shí)，后每隔30 min監(jiān)測一次，直至當(dāng)日下午19時(shí).

表1 某河流斷面水文水質(zhì)參數(shù)Tab.1 Model parameters of study area

圖1 斷面A污染物觀測濃度序列Fig.1 Simulation observe concentration serial of site A

3.1 反演分析

假定污染物從岸邊瞬間排入河道，泄漏量為M，泄漏地點(diǎn)距離A斷面S0，泄漏時(shí)間距離首次監(jiān)測時(shí)間為T0，參數(shù)“真值”見表1.利用表1中的3個(gè)模型參數(shù)真值及其他參數(shù)計(jì)算出斷面A的污染物濃度分布，初始監(jiān)測時(shí)間為發(fā)現(xiàn)污染帶后的第一時(shí)間，某日下午15時(shí)，后每隔30min監(jiān)測一次，直至當(dāng)日下午19時(shí)，見圖1.

從圖1可以看出，9組A斷面污染物濃度序列近似正態(tài)分布，最高值發(fā)生于該日16時(shí)30分，污染物質(zhì)量濃度達(dá)到1.62mg／L，到19時(shí)已經(jīng)降低為0.19 mg／L.據(jù)先驗(yàn)信息，該突發(fā)污染事故可能為交通事故導(dǎo)致的污染物泄漏或化學(xué)品企業(yè)儲(chǔ)罐爆炸泄漏瞬間排入河流造成的污染.

設(shè)坐標(biāo)中心為污染源泄漏點(diǎn)，測點(diǎn)A坐標(biāo)為（S0，0），待反演參數(shù)定義為X（M，S0，T0）.根據(jù)先驗(yàn)信息，3個(gè)待反演模型參數(shù)先驗(yàn)分布為均勻分布，其對應(yīng)的先驗(yàn)概率密度函數(shù)分別為：

假定觀測誤差為白噪聲N（0，σ2），σ＝0.01，條件概率密度表示為：

根據(jù)貝葉斯定理，模型參數(shù)的后驗(yàn)概率密度函數(shù)為：

式中：λ為一比例常數(shù).

3.2 反演結(jié)果

取Proposal分布為均勻分布q（X＊｜X（i））＝U（X（i）－ΔX，X（i）＋ΔX），步長為先驗(yàn)范圍的5%，取ΔX（M，S0，T0）＝（60 000，450，400），按照Metropolis算法反演思路迭代10 000次.3個(gè)參數(shù)的迭代曲線見圖2，參數(shù)后驗(yàn)概率直方圖見圖3.

圖2 參數(shù)反演迭代曲線Fig.2 Iteration curve of model parameters

從圖2可看出，迭代約2 000次后，M，S0，T0均達(dá)到穩(wěn)定，馬爾科夫鏈進(jìn)入穩(wěn)定收斂區(qū)域.從圖3可看出，3個(gè)模型參數(shù)后驗(yàn)概率密度呈現(xiàn)不太明顯的正態(tài)分布；污染源源強(qiáng)M在1 000 000g附近的頻度最大，完全符合預(yù)定設(shè)定初值；污染源距離A斷面距離S0在5 000m附近的頻度最大，與預(yù)先設(shè)定值完全符合；而污染泄漏時(shí)間距A斷面首次監(jiān)測時(shí)間T0在7 000～8 000s區(qū)間頻度最大，與預(yù)先設(shè)定值較為接近.

圖3 參數(shù)后驗(yàn)概率直方圖Fig.3 Posterior histogram of model parameters

3.3 反演效果分析

為進(jìn)一步驗(yàn)證反演結(jié)果，剔除前面2 000次結(jié)果，對其他8 000次模型參數(shù)進(jìn)行后驗(yàn)統(tǒng)計(jì)分析，分析結(jié)果見表2.可以看出，3個(gè)參數(shù)的數(shù)學(xué)期望誤差、中值誤差均較小，而以污染源強(qiáng)度M的誤差最小.

表2 污染源參數(shù)后驗(yàn)統(tǒng)計(jì)結(jié)果Tab.2 Posterior statistical results of model parameters

為了檢查基于MCMC－Bayes方法的反演穩(wěn)定性，采取多次間隔性反演提取模型參數(shù)反演誤差序列，見圖4.可以看出，反演結(jié)果較穩(wěn)定，3個(gè)參數(shù)的相對誤差均在9%以內(nèi)，平均誤差分別為2.94%，2.50%和2.21%.比較而言，污染源源強(qiáng)M的誤差起伏較大，位于［0.48%，8.48%］；污染源位置S0的誤差起伏次之，位于［2.50%，5.90%］；污染源泄漏時(shí)間T0的誤差起伏最小，位于［2.21%，4.65%］.

圖4 基于MCMC－Bayes問題反演誤差分析Fig.4 Error serial with method of MCMC－Bayes

作為比較，本文采用混合遺傳模式搜索算法對該問題進(jìn)行求解，反演誤差見圖5.可以看出，采用混合遺傳模式搜索算法求解的問題解誤差較大，最高誤差達(dá)到24.7%，3個(gè)參數(shù)的平均誤差分別達(dá)到3.20%，4.08%，7.45%.其中污染源泄漏時(shí)間T0的誤差起伏最大，位于［1.21%，24.7%］；污染源源強(qiáng)M的誤差起伏次之，位于［0.54%，10.69%］；污染源位置S0的誤差起伏最小，位于［0.68%，13.33%］.

圖5基于混合遺傳模式搜索算法的問題反演誤差分析Fig.5 Error serial with method of hybrid GA－Patter search

兩種方法反演的模型參數(shù)誤差比較見表3.很明顯，基于MCMC－Bayes的源識(shí)別誤差遠(yuǎn)小于優(yōu)化算法，3個(gè)參數(shù)的最小誤差、最大誤差和平均誤差都呈現(xiàn)明顯的特征，可見，相對優(yōu)化算法而言，基于MCMC－Bayes的源識(shí)別方法具有良好的可靠性和穩(wěn)定性.

表3 基于MCMC－Bayes 和混合遺傳模式搜索算法（GA－PS）的反演誤差對照Tab.3 Error comparison with method of MCMC－Bayes and GA－PS%

4 結(jié) 論

污染源識(shí)別對于開展突發(fā)污染事故應(yīng)急處置、加強(qiáng)水質(zhì)管理和控制具有重要的現(xiàn)實(shí)意義，而由于測量數(shù)據(jù)有限且?guī)в性肼暤仍蚴沟梦廴驹醋R(shí)別成為一個(gè)不適定問題.本研究基于貝葉斯推理和二維水質(zhì)模型建立水體污染識(shí)別數(shù)學(xué)模型，運(yùn)用馬爾科夫鏈蒙特卡羅法抽樣獲得了污染源源強(qiáng)、污染源位置和污染泄漏時(shí)間3個(gè)模型參數(shù)的后驗(yàn)概率分布和統(tǒng)計(jì)結(jié)果.實(shí)例研究結(jié)果表明，基于馬爾科夫鏈蒙特卡羅抽樣算法的貝葉斯推理可以較好地用來實(shí)現(xiàn)水體污染識(shí)別，具有識(shí)別精度高、誤差小的特點(diǎn)，其可靠性和穩(wěn)定性高于混合遺傳模式搜索優(yōu)化算法.

總體看來，MCMC－Bayes通過構(gòu)造Markov鏈進(jìn)行隨機(jī)動(dòng)態(tài)模擬，實(shí)現(xiàn)對高維空間無明確數(shù)學(xué)表達(dá)式概率分布密度函數(shù)的數(shù)值計(jì)算，且具備計(jì)算速度快、復(fù)雜度不依賴于計(jì)算空間維數(shù)等特點(diǎn).同時(shí)，它可以方便地將各種先驗(yàn)信息和誤差信息高效地融合到問題求解過程中，減小問題的不確定性，獲得全局最可能解，可以很好地用來開展同屬問題不適定的水體污染反問題研究.

［1］ 戴會(huì)超，王玲玲.工程水力學(xué)反問題研究及應(yīng)用［J］.四川大學(xué)學(xué)報(bào)：工程科學(xué)版，2006，38（1）：15－19.DAI Hui－chao，WANG Ling－ling.Study and application of the inverse problem on hydraulics［J］.Journal of Sichuan University：Engineering Science Edition，2006，38（1）：15－19.（In Chinese）

［2］ 劉繼軍.不適定問題的正則化方法及應(yīng)用［M］.北京：科學(xué)出版社，2005：42－46.LIU Ji－jun.Regularization methods of ill－posed problems and its application［M］.Beijing：Science Press，2005：42－46.（In Chinese）

［3］ 閔濤，周孝德，張世梅，等.對流擴(kuò)散方程源項(xiàng)識(shí)別反問題的遺傳算法［J］.水動(dòng)力學(xué)研究與進(jìn)展，A輯，2004，19（4）：520－524.MIN Tao，ZHOU Xiao－de，ZHANG Shi－mei，et al.Genetic algorithm to an inverse problem of source term identification for convection－diffusion equation［J］.Journal of Hydrodynamic Ser A，2004，19（4）：520－524.（In Chinese）

［4］ 王澤文，邱淑芳.一類流域點(diǎn)污染源識(shí)別的穩(wěn)定性與數(shù)值模擬［J］.水動(dòng)力學(xué)研究與進(jìn)展，A輯，2008，23（4）：364－371.WANG Ze－wen，QIU Shu－fang.Stability and numerical simulation of pollution point source identification in a watershed［J］.Journal of Hydrodynamic Ser A，2008，23（4）：364－371.（In Chinese）

［5］ 嚴(yán)齊斌.河流水質(zhì)參數(shù)估計(jì)的蒙特卡羅方法［J］.水利水電技術(shù)，2006，37（10）：14－16.YAN Qi－bin.Monte Carlo method for estimation on parameters of river water quality［J］.Water Resources and Hydropower Engineering，2006，37（10）：14－16.（In Chinese）

［6］ 朱嵩，劉國華，毛根海，等.利用貝葉斯推理估計(jì)二維含源對流擴(kuò)散方程參數(shù)［J］.四川大學(xué)學(xué)報(bào)：工程科學(xué)版，2008，40（2）：38－43.ZHU Song，LIU Guo－h(huán)ua，MAO Geng－h(huán)ai，et al.Application of Bayesian inference to estimate the parameters in 2Dconvection－diffusion equation with source［J］.Journal of Sichuan University：Engineering Science Edition，2008，40（2）：38－43.（In Chinese）

［7］ 劉曉東，姚琪，薛紅琴，等.環(huán)境水力學(xué)反問題研究進(jìn)展［J］.水科學(xué)進(jìn)展，2009，20（6）：885－890.LIU Xiao－dong，YAO Qi，XUE Hong－qin，et al.Advance in inverse problems of environmental hydraulics［J］.Advance in Water Science，2009，20（6）：885－890.（In Chinese）

［8］ 曹小群，宋君強(qiáng)，張衛(wèi)民，等.對流擴(kuò)散方程源項(xiàng)識(shí)別反問題的MCMC方法［J］.水動(dòng)力學(xué)研究與進(jìn)展，A輯，2010，25（2）：129－135.CAO Xiao－qun，SONG Jun－qiang，ZHANG Wei－min，et al.MCMC method on an inverse problem of source term identification for convection－diffusion equation［J］.Journal of Hydrodynamic Ser A，2010，25（2）：129－135.（In Chinese）

［9］ ANDRIEU C，F(xiàn)REITAS N D，DOUCET A，et al.An introduction to MCMC for machine learning［J］.Machine Learning，2003，50：5－43.