基于河流示蹤實(shí)驗(yàn)的 Bayes 污染溯源:算法參數(shù)、影響因素及頻率法對(duì)比

姜繼平,董芙嘉,劉仁濤,袁一星

?

基于河流示蹤實(shí)驗(yàn)的 Bayes 污染溯源:算法參數(shù)、影響因素及頻率法對(duì)比

姜繼平1,2*,董芙嘉1,劉仁濤1,袁一星1

(1.哈爾濱工業(yè)大學(xué)環(huán)境學(xué)院,黑龍江哈爾濱150090；2.南方科技大學(xué)環(huán)境科學(xué)與工程學(xué)院,廣東深圳518055)

基于貝葉斯理論,結(jié)合濃度時(shí)間序列方差假定和Adaptive Metropolis MCMC后驗(yàn)采樣,建立了用于突發(fā)水污染應(yīng)急溯源的貝葉斯推理方法.該方法結(jié)合經(jīng)驗(yàn)知識(shí)和監(jiān)測(cè)事實(shí)對(duì)源項(xiàng)參數(shù)的分布進(jìn)行推理,直接對(duì)溯源結(jié)果的反向不確定性用概率分布形式進(jìn)行表征.依據(jù)河流實(shí)地示蹤劑實(shí)驗(yàn)案例,對(duì)Bayesian推理溯源的實(shí)際效果、后驗(yàn)參數(shù)相關(guān)性和影響因素等方面進(jìn)行了驗(yàn)證和測(cè)試,結(jié)果表明源項(xiàng)參數(shù)和方差的后驗(yàn)概率密度的偏度對(duì)方差假定敏感,且得到關(guān)鍵參數(shù)推薦值:使用RMSE為目標(biāo)函數(shù);異方差假定中穩(wěn)定化因子λ1為0, λ2為0.1~0.5;AM采樣建議比例因子sd選擇0.1~0.3.并對(duì)貝葉斯方法和傳統(tǒng)基于優(yōu)化的頻率法在求解思路、計(jì)算過(guò)程、溯源結(jié)果等角度進(jìn)行了深層次的辨析.本研究相關(guān)結(jié)果為貝葉斯推理技術(shù)在污染溯源的實(shí)際應(yīng)用中提供了較為重要的參考價(jià)值.

貝葉斯推理；污染源反演；突發(fā)水污染；AM采樣；河流示蹤劑試驗(yàn)

當(dāng)前,蓄意或非故意的化學(xué)品泄露造成的水污染事件仍不斷發(fā)生.作為典型發(fā)展中國(guó)家,我國(guó)地表水環(huán)境一直受到化學(xué)品泄露事故的威脅[1-2].其他發(fā)展中國(guó)家和發(fā)達(dá)國(guó)家,如南非[3]、印度[4]、美國(guó)[5]等,也面臨類似的困擾和挑戰(zhàn).許多河流污染發(fā)生時(shí),污染源未知,肉眼難以辨識(shí).同時(shí)不少企業(yè)依然存在夜間偷排行為,需要開(kāi)發(fā)基于監(jiān)測(cè)數(shù)據(jù)的污染源反演技術(shù)進(jìn)行污染溯源.確定污染排放信息是開(kāi)展應(yīng)急處置和風(fēng)險(xiǎn)防控的前提[6].

污染溯源技術(shù)的理論基礎(chǔ)是求解特定情景下的源項(xiàng)反演問(wèn)題,求解過(guò)程融入污染物輸移模型.理論上如果對(duì)污染輸移模型采用負(fù)時(shí)間步和空間步進(jìn)行數(shù)值計(jì)算,應(yīng)該能求得污染輸移歷史過(guò)程,然而在現(xiàn)實(shí)中由于監(jiān)測(cè)誤差和數(shù)值誤差的存在,其不可控的誤差反向傳播導(dǎo)致了反算結(jié)果極大地偏離實(shí)際值.該特征就是反問(wèn)題不適定性的典型表現(xiàn),也正因?yàn)榇?在各領(lǐng)域中出現(xiàn)了形形色色的反演技術(shù)來(lái)求解反問(wèn)題.

近幾年,河流污染源項(xiàng)反演問(wèn)題開(kāi)始得到關(guān)注,中國(guó)學(xué)者做了大量工作,提出各類求解方法.例如,解析法[7],基于遺傳算法、微分進(jìn)化算法等的優(yōu)化法[8-10],地學(xué)統(tǒng)計(jì)方法[11],反向位置概率密度函數(shù)法[12-13],伴隨方程法[14]等.都是基于相同的污染物輸移、稀釋和轉(zhuǎn)化機(jī)制,求得源項(xiàng)參數(shù)的唯一值,反演結(jié)果反向不確定性需要用攝動(dòng)法等多次計(jì)算進(jìn)行表征.還有從純數(shù)學(xué)角度進(jìn)行污染溯源反問(wèn)題研究[15-17].2015年后,有人開(kāi)始從實(shí)踐角度出發(fā),關(guān)注反演的影響因素分析,例如對(duì)寬淺河道[18]和二維河道概化下[19]瞬時(shí)排放源項(xiàng)反演精度的影響分析.

應(yīng)用Bayes方法解決科學(xué)和工程領(lǐng)域的參數(shù)估計(jì)研究發(fā)展迅速,例如水文學(xué)中的地表徑流模擬和參數(shù)估計(jì)[20-22]、土壤環(huán)境化學(xué)中物化機(jī)制的解釋[23]等等.近幾年,也有人開(kāi)始將Bayes方法用于水污染溯源,如朱崇[24-26]、毛星[27],陳海洋等[28],曹小群等[29],Wei等[30].有別于確定性的方法,Bayesian反演框架采用信息驅(qū)動(dòng),組合概率性的先驗(yàn)信息和觀測(cè)數(shù)據(jù)中所包含信息,通過(guò)貝葉斯推理更新先驗(yàn)信息,得到后驗(yàn)條件概率分布,是不確定性的反演方法,在應(yīng)急決策風(fēng)險(xiǎn)的不確定性量化方面更有優(yōu)勢(shì).

朱嵩等結(jié)合Metroplis-Hastings(MH)采樣算法,在假想的10m×4m區(qū)域內(nèi)進(jìn)行數(shù)值實(shí)驗(yàn),反演污染源位置和源強(qiáng)[24,26],這是Bayesian污染溯源較早的嘗試,但它更多是框架性、演示性.2012年陳海洋等[28]對(duì)Bayesian溯源方法進(jìn)行了較為細(xì)致的推導(dǎo)分析,對(duì)關(guān)鍵的誤差假設(shè)問(wèn)題進(jìn)行了初步交代.同樣采用MH采樣算法,基于更為接近真實(shí)水體參數(shù)的假想案例對(duì)污染源位置、強(qiáng)度河時(shí)刻同時(shí)反演,并同優(yōu)化算法進(jìn)行對(duì)比.曹小群等[29]首次采用Adative Metroplis(AM)算法基于假想案例,考慮多點(diǎn)源情形,對(duì)污染源位置和強(qiáng)度進(jìn)行反演.Wei等[30]結(jié)合了正向模型參數(shù)的不確定性分析,采用AM算法,基于實(shí)際河道假想排放情景和數(shù)據(jù),對(duì)源項(xiàng)三參數(shù)進(jìn)行了反演.

可以看出,目前研究未有人基于實(shí)地示蹤劑實(shí)驗(yàn)和實(shí)際的試驗(yàn)監(jiān)測(cè)數(shù)據(jù)對(duì)Bayesian污染源反演方法進(jìn)行實(shí)踐測(cè)試.特別是需要針對(duì)不同特征河流、不同污染排放類型、不同濃度監(jiān)測(cè)精度等實(shí)際情形進(jìn)行對(duì)比分析.Bayeisan污染溯源推導(dǎo)過(guò)程中關(guān)鍵的誤差假定形式?jīng)]有深入探討分析,這將影響到似然函數(shù)的構(gòu)建,涉及方法的應(yīng)用范圍.在實(shí)際河流尺度和污染排放下,算法參數(shù)設(shè)置和影響因素也未有報(bào)道.Byesian方法同以優(yōu)化法為代表的確定性方法在計(jì)算原理、過(guò)程和結(jié)果等方面,也未有人進(jìn)行系統(tǒng)辨析,也有必要結(jié)合案例深入探討.

為此,本文針對(duì)突發(fā)水污染應(yīng)急中的污染溯源問(wèn)題,開(kāi)展基于河流示蹤劑實(shí)驗(yàn)的貝葉斯污染溯源研究.采用Adaptive Metropolis算法進(jìn)行后驗(yàn)概率密度采樣,驗(yàn)證與分析不同誤差假定下貝葉斯反演技術(shù)的有效性、反向不確定性和影響因素,給出操作參數(shù)推薦值.并同基于優(yōu)化的頻率法進(jìn)行對(duì)比,深入剖析兩者思路、過(guò)程和結(jié)果的差異.

1 基于Bayesian推理的不確定性應(yīng)急溯源方法構(gòu)建

1.1 正向模型

基于污染物對(duì)流擴(kuò)散模型,假設(shè)模型參數(shù)為常數(shù),污染物在水平和垂直方向上完全混合,通過(guò)特征線等方法可以得到不同條件下的水質(zhì)模型解析解.

在監(jiān)測(cè)點(diǎn)(,)處的污染物濃度可由下式計(jì)算:

式中:為河流斷面面積,m2;D為縱向平均擴(kuò)散系數(shù),m2/min;為河流平均流速,m;為衰減系數(shù),min-1;為監(jiān)測(cè)時(shí)刻,min;為監(jiān)測(cè)點(diǎn)位置,m.

考慮到污染溯源研究中源位置未知,未將坐標(biāo)原點(diǎn)建立在源排放處和排放時(shí)刻.慮到研究問(wèn)題的特征時(shí)間尺度,將時(shí)間單位設(shè)定為min或者h(yuǎn).

同時(shí)對(duì)于連續(xù)排放源的情境,也可推導(dǎo)出解析解,具體形式可參考文獻(xiàn)[31].為對(duì)比不同污染傳輸機(jī)制概化的影響,還采用了EPD-RIV1模型[32]和OTIS模型[33]兩個(gè)基于數(shù)值解的一維水質(zhì)模型.

1.2 應(yīng)急溯源Bayes理論框架的建立

1.2.1 污染源溯源問(wèn)題的廣義陳述 在河流某斷面檢出污染物嚴(yán)重超標(biāo),在該斷面上游開(kāi)展應(yīng)急監(jiān)測(cè),獲得個(gè)時(shí)空采樣點(diǎn)的濃度數(shù)據(jù),它們和污染物傳輸機(jī)理模型間通過(guò)正向模型算子建立聯(lián)系(式2).該算子即為將模型映射到數(shù)據(jù)空間的函數(shù).實(shí)際上,正向模型算子通常只能是真實(shí)物理過(guò)程的近似,當(dāng)使用函數(shù)來(lái)描述該物理過(guò)程時(shí)存在系統(tǒng)偏差,可用維向量model來(lái)表示.同樣也存在一個(gè)維向量meas,表示隨機(jī)測(cè)量誤差.

式中:為位置污染源參數(shù),為ISP問(wèn)題中的模型參數(shù).污染溯源的目標(biāo)是已知監(jiān)測(cè)數(shù)據(jù)向量估計(jì),或者是的函數(shù)(s)[34].需要注意的是,可能是模型邊界條件和初始條件而不是模型表達(dá)式中的參數(shù).

1.2.2 面向應(yīng)急溯源的Bayes理論框架 基于污染溯源問(wèn)題的Bayes公式如下:

式中:()是源項(xiàng)參數(shù)的后驗(yàn)分布,(|)是先驗(yàn)分布,()是似然函數(shù),()是證據(jù),為源項(xiàng),為濃度監(jiān)測(cè)數(shù)據(jù),為背景信息,即用于確定先驗(yàn)分布的信息.

1.2.3 Bayes推理應(yīng)急溯源方法的建立 基于模型的Bayes推理框架常分為四步:模型構(gòu)建、后驗(yàn)分布計(jì)算、分析后驗(yàn)分布和決策推理[35].

1.2.3.1 第一步:模型構(gòu)建

(1)溯源問(wèn)題的Bayesian形式化

假定這些誤差獨(dú)立(但并不一定同分布),有:

若進(jìn)一步假定這兩類誤差都是獨(dú)立同分布,即所有的meas,i都等于meas,所有的model,i等于model.可推導(dǎo)出似然函數(shù)計(jì)算公式(7),具體過(guò)程可參考[36-37].

對(duì)于異方差非相關(guān)誤差項(xiàng),可通過(guò)對(duì)監(jiān)測(cè)數(shù)據(jù)進(jìn)行如下轉(zhuǎn)換而穩(wěn)定化.Box and Cox[38]證實(shí)了它在水文數(shù)據(jù)中的適用性.

這里,1能從平均殘差中估計(jì)得到.如果區(qū)間平均殘差隨總平均殘差線性增加,就設(shè)置1為0.5.如果方差是二次方增加,就設(shè)置1為0[20].

基于異方差非相關(guān)誤差項(xiàng)假定的似然函數(shù)可以寫(xiě)為:

式中:是X的方差.

(2)設(shè)置先驗(yàn)概率

用Bayes方法處理統(tǒng)計(jì)推理問(wèn)題的前提是給出參數(shù)的先驗(yàn)分布.在文獻(xiàn)中提及的方法有客觀法、主觀概率法、同等無(wú)知原則法、共軛先驗(yàn)分布等.在水文學(xué)中,常用客觀法即經(jīng)驗(yàn)貝葉斯法,基于大量歷史觀測(cè)數(shù)據(jù),如徑流量等,得到先驗(yàn)分布,常用的有Gamma分布,正態(tài)分布,均一分布等.

河流污染應(yīng)急溯源中的先驗(yàn)信息常常有限,因此選擇均一分布作為優(yōu)先考慮的先驗(yàn)概率密度函數(shù).當(dāng)然,如果河流某區(qū)段存在大量風(fēng)險(xiǎn)源,如化工區(qū)、養(yǎng)殖場(chǎng)等,可優(yōu)先考慮該污染源出現(xiàn)在該區(qū)段的可能性,從而設(shè)置綜合的概率密度函數(shù).

均一分布的先驗(yàn)概率密度可如下描述:

(|)=常數(shù),∈Real (10)

(3)后驗(yàn)概率密度函數(shù)

基于同方差非相關(guān)誤差項(xiàng)假設(shè)的后驗(yàn)概率密度函數(shù)可由如下公式描述:

式中:(.)為指示函數(shù).

1.2.3.2 第二步:計(jì)算后驗(yàn)分布-MCMC采樣

后驗(yàn)分布常常同先驗(yàn)分布是非共軛的,難以近似求解或者是全條件分布,不像已知的任何一個(gè)分布.Markov Chain Monte Carlo(MCMC)常被用來(lái)快速估計(jì)后驗(yàn)分布.最常用的兩種MCMC方法是:Metroplis-Hastings(MH)算法[39-40]和Gibbs采樣[41].

Marshall等[20]經(jīng)過(guò)一個(gè)全面的MCMC算法對(duì)比實(shí)驗(yàn)后,建議在降雨-徑流機(jī)理模型中使用Adaptive Metropolis(AM)算法.本文也采用AM算法作為采樣方法.

1.2.3.3 第三步:分析后驗(yàn)概率密度

MCMC采樣最后收斂到后驗(yàn)概率密度(或者其對(duì)數(shù)形式),對(duì)樣本進(jìn)行統(tǒng)計(jì),如均值、方差、中值、分位數(shù)、偏度等可以得到對(duì)后驗(yàn)概率密度的描述性統(tǒng)計(jì)分析結(jié)果.另外,也可以采用數(shù)值積分的方法,對(duì)各參數(shù)的邊緣概率密度進(jìn)行統(tǒng)計(jì).例如,對(duì)于s有:

1.2.3.4 第四步:對(duì)結(jié)果進(jìn)行推斷

依據(jù)上述統(tǒng)計(jì)分析結(jié)果,對(duì)污染源信息做最后推斷.本文采用中值和貝葉斯區(qū)間對(duì)源參數(shù)結(jié)果進(jìn)行推斷.此外,在實(shí)際應(yīng)急響應(yīng)過(guò)程中,還需要結(jié)合物理方法進(jìn)行實(shí)地驗(yàn)證,才能最后確定污染源信息.

1.3 后驗(yàn)概率密度采樣的AM算法

1.3.1 Markov鏈及其收斂 Markov鏈?zhǔn)且粋€(gè)未來(lái)狀態(tài)只和當(dāng)前狀態(tài)有關(guān)而與過(guò)去狀態(tài)相獨(dú)立的隨機(jī)過(guò)程.數(shù)學(xué)上的定義及Markov鏈?zhǔn)諗康椒€(wěn)態(tài)分布的嚴(yán)格數(shù)學(xué)分析和證明可以參見(jiàn)文獻(xiàn)[42].

1.3.2 Adaptive Metropolis算法 Adaptive Metropolis (AM) 算法是一種較好的改進(jìn)后的MH算法[26,43].其特征是MH算法中的建議PDF基于參數(shù)后驗(yàn)協(xié)方差得到.該后驗(yàn)協(xié)方差矩陣是由過(guò)去的迭代結(jié)果計(jì)算得到,而且每一步都會(huì)計(jì)算該協(xié)方差矩陣.這樣,建議分布通過(guò)剛剛獲取的后驗(yàn)分布信息進(jìn)行更新.在第步[43],考慮用一個(gè)多變量正態(tài)分布來(lái)表示建議PDF.其均值用當(dāng)前樣本的平均值,而協(xié)方差采用矩陣.該協(xié)方差矩陣在前面0步中已經(jīng)固定了一個(gè)值0并且通過(guò)下面的方式進(jìn)行后續(xù)的更新:

式中:是一個(gè)數(shù)量很小的參數(shù),被用來(lái)確保B不會(huì)變得奇異,sd是比例參數(shù),依賴于參數(shù)向量的維數(shù),確保建議狀態(tài)有合理的接受率(例如25%~75%).第+1個(gè)迭代步的協(xié)方差計(jì)算成本較低,它符合下面的公式:

建議協(xié)方差的計(jì)算需要定義一個(gè)任意的初始協(xié)方差0.為讓該過(guò)程自動(dòng)進(jìn)行,設(shè)置這個(gè)初始協(xié)方差為先驗(yàn)分布條件下的參數(shù)的初始協(xié)方差,例如最初的5%迭代步長(zhǎng)的參數(shù)的協(xié)方差.

綜上,實(shí)現(xiàn)AM算法的步驟如下:

(1)初始化=0

(2)a. 為當(dāng)前迭代步選擇B,采用公式(13)計(jì)算.

b. 為生成建議參數(shù)值*,其中*~N(θ,B).

c. 計(jì)算接受率:

式中:(|) 是似然函數(shù),()是的先驗(yàn)分布.

d. 生成~U[0,1].

e. 如果<,接受θ+1=*,否則設(shè)置θ+1=θ.

(3)重復(fù)步驟(2).

對(duì)于參數(shù)值在約束區(qū)間外的情形,設(shè)置接受率為零.AM算法的一個(gè)最大優(yōu)點(diǎn)是整個(gè)參數(shù)集同時(shí)被更新,減少了計(jì)算時(shí)間和復(fù)雜度.這對(duì)于參數(shù)高度相關(guān)的情形尤為適用(后面可以看到,溯源計(jì)算會(huì)出現(xiàn)此情形).它與大多數(shù)經(jīng)典算法不一樣,后者對(duì)每一采樣得到的參數(shù)需要有不同的建議分布,致使在參數(shù)分塊采樣(block)時(shí)極大增加了復(fù)雜度.AM采樣的有效性和遍歷性特征,可參考文獻(xiàn)[44].

1.3.3 Markov鏈的收斂診斷 診斷Markov鏈?zhǔn)諗糠椒ㄒ呀?jīng)很多[45],其中最為常用是Gelman 和Rubin 1992年開(kāi)發(fā)以及Raftery和Lewis于同年開(kāi)發(fā)的統(tǒng)計(jì)方法.本文采用Gelman-Rubin潛在規(guī)模減縮因子(potential scale reduction factor, PSRF)來(lái)診斷模型的收斂.R的計(jì)算如下所示:

式中:是次平行采樣鏈條的平均值的方差,是個(gè)鏈內(nèi)方差平均值, df是漸進(jìn)Student t分布的自由度.更多信息可參見(jiàn)文獻(xiàn)[45].需要注意的是,PSRF統(tǒng)計(jì)是基于多條鏈平行采樣的結(jié)果,假如是單一長(zhǎng)鏈采樣,可以把鏈分為三段進(jìn)行統(tǒng)計(jì).另外通過(guò)樣本軌跡圖也可以對(duì)收斂性做出定性分析.

2 基于實(shí)地示蹤實(shí)驗(yàn)的污染溯源

2.1 示蹤劑實(shí)驗(yàn)

結(jié)合論文研究目標(biāo),尋找到USGS開(kāi)展的三個(gè)典型河流示蹤劑實(shí)驗(yàn)用于污染溯源研究,重新構(gòu)建溯源情景并整理其監(jiān)測(cè)數(shù)據(jù).其中以Truckee River示蹤劑實(shí)驗(yàn)為主,用于算法驗(yàn)證; Lagan River 和West Hobolochitto Creek示蹤劑實(shí)驗(yàn)為輔,用于對(duì)比和影響因素分析.這三個(gè)案例有著不同的河道形態(tài)與水動(dòng)力特征及污染物排放特征,案例分別標(biāo)識(shí)為Case-T1,T2和T3.

2.1.1 Truckee River示蹤劑實(shí)驗(yàn) Truckee River位于美國(guó)加利福尼亞州Tahoe City和內(nèi)華達(dá)州Vista之間.Rivord等在2006~2007年間進(jìn)行了三次示蹤劑實(shí)驗(yàn),為管理河流發(fā)生的突發(fā)污染事件和污染傳輸預(yù)測(cè)提供數(shù)據(jù)支持.當(dāng)時(shí)河流流量介于143ft3/s和2660ft3/s之間.熒光羅丹明B(RWT),即示蹤的‘污染物’,沿此河流三個(gè)位置注入,在下游監(jiān)測(cè)其流經(jīng)過(guò)程,具體操作和說(shuō)明在文獻(xiàn)[47]和[48]中有詳細(xì)介紹.基于上述示蹤劑實(shí)驗(yàn),Rivord等人率定并建立了Truckee River的OTIS模型.

本文采用Truckee River中段進(jìn)行的示蹤劑實(shí)驗(yàn).中段長(zhǎng)約44km,設(shè)置4個(gè)采樣點(diǎn)(6~9號(hào)站點(diǎn)).2006年7月29日6點(diǎn)5分,1.3kg RWT在5號(hào)點(diǎn)釋放,接近加利福尼亞特拉基的Glenshire Driver,示蹤劑濃度數(shù)據(jù)在6~9號(hào)點(diǎn)檢測(cè)得到.

2.1.2 Lagan River示蹤劑實(shí)驗(yàn) 2004年11月與2005年4月間,在北愛(ài)爾蘭的Lagan River進(jìn)行了8次示蹤劑投放實(shí)驗(yàn),用來(lái)測(cè)量河流的復(fù)氧系數(shù)[49],試驗(yàn)河段長(zhǎng)2.6km.2004年12月10號(hào)進(jìn)行的B號(hào)測(cè)試監(jiān)測(cè)數(shù)據(jù)用于溯源反演.惰性氣體氙、氪在6個(gè)大氣壓下通入水中耗時(shí)30s,31.8min后,在下游1200m 和2600m 處的站點(diǎn)1和2監(jiān)測(cè)穿透曲線,采樣每隔數(shù)min進(jìn)行1次.更多細(xì)節(jié)見(jiàn)文獻(xiàn)[49].

2.1.3 West Hobolochitto Creek示蹤劑實(shí)驗(yàn) Rathbun于1975年7月在美國(guó)密西西比州West Hobolochitto Creek進(jìn)行了多次氣體示蹤劑實(shí)驗(yàn)來(lái)測(cè)量大氣復(fù)氧速率系數(shù)[46].West Hobolochitto Creek為一條水文特征復(fù)雜既有死水區(qū)又有激流的小溪.研究河段長(zhǎng)4190m,實(shí)驗(yàn)時(shí)河流處在低流量條件.

本研究采用了Rathbun的“2號(hào)實(shí)驗(yàn)”. Rathbun通過(guò)多孔管擴(kuò)散器在100分鐘左右的時(shí)間內(nèi)將1.9kg乙烯示蹤劑氣體注入到河流中.同時(shí),水溶性RWT在同一地點(diǎn)注入.兩示蹤劑的濃度在注入點(diǎn)下游五個(gè)站點(diǎn)進(jìn)行監(jiān)測(cè).有關(guān)示蹤實(shí)驗(yàn)更多細(xì)節(jié)可參考[46].本研究將用于溯源反演的應(yīng)急監(jiān)測(cè)情景設(shè)置為:在示蹤劑釋放點(diǎn)下游2680m、3370m、4190m處設(shè)置3個(gè)監(jiān)測(cè)站點(diǎn)在示蹤劑排放后的389分鐘時(shí)開(kāi)始監(jiān)測(cè),持續(xù)了3.5h,每十分鐘采樣1次.

2.2 正向模型參數(shù)化

采用Truckee River中段示蹤劑實(shí)驗(yàn)案例(Case-T1)進(jìn)行溯源反演.正向模型首先采用點(diǎn)源污染瞬時(shí)排放一維模型(式1),模型參數(shù)設(shè)置見(jiàn)表1.

表1 Truckee River案例源項(xiàng)參數(shù)和正向模型設(shè)置

2.3 采樣算法參數(shù)設(shè)置

表2 Truckee River案Bayesian溯源中AM采樣參數(shù)設(shè)置

污染源參數(shù)分布設(shè)定s~U(100,5000),s= U(-30000, -1000),s=U(-300,-10),其中U為均一分布.先驗(yàn)概率密度就是它們的聯(lián)合分布.同時(shí)先假定誤差均方差非相關(guān),似然函數(shù)用(7)式表示.采樣次數(shù)設(shè)定為100,000,起初20,000次不參與最后樣本統(tǒng)計(jì),只取后面80,000個(gè)Markov鏈樣本進(jìn)行分析,考察后驗(yàn)概率密度及不確定性.如果在初次運(yùn)行中接受率未處于25%~75%,需要人工交互通過(guò)調(diào)節(jié)建立比例因子(proposal scaling factor, sd)來(lái)調(diào)整建議密度.另外,算法實(shí)現(xiàn)中對(duì)似然函數(shù)進(jìn)行了對(duì)數(shù)處理,其他參數(shù)設(shè)置在表2中.

2.4 溯源反演結(jié)果

Case-T1溯源計(jì)算結(jié)果顯示總體PBIAS(平均值百分比偏差)為95%,接受率為35%,在可接受范圍內(nèi).污染源反演后驗(yàn)分布的概要統(tǒng)計(jì)量列在表3中.用樣本均值表示對(duì)污染源的參數(shù)估計(jì).

可見(jiàn)溯源結(jié)果與真實(shí)污染源情況非常接近而且標(biāo)準(zhǔn)差很小,Bayesian算法很好的完成了應(yīng)急溯源.從偏度(skewness)、均值(mean)和中值(0.5)的一致性可以看出,s,s,s的邊緣PDF基本對(duì)稱,但方差2呈正偏態(tài).為0.05的最高概率密度區(qū)間,即貝葉斯區(qū)間(Bayesian interval),表明s不確定性最大,而s貝葉斯區(qū)間最小.

表3 Case-T1(Truckee River)反演結(jié)果的概要統(tǒng)計(jì)量

圖1 監(jiān)測(cè)濃度和正向預(yù)測(cè)濃度值對(duì)比

PDF曲線通過(guò)高斯核對(duì)Markov鏈中后80,000個(gè)樣本進(jìn)行估計(jì)得到.PDF曲線形態(tài)也定性佐證了上述分析結(jié)論.

圖1表示在后驗(yàn)概率密度極值處計(jì)算得到的污染源參數(shù)濃度值與實(shí)測(cè)濃度值,數(shù)據(jù)對(duì)稱集中地分布在=直線兩側(cè).表3也對(duì)后驗(yàn)分布百分位數(shù)0.025、0.5和0.975進(jìn)行了統(tǒng)計(jì).全部迭代過(guò)程后驗(yàn)參數(shù)及其均值迭代軌跡繪于圖3和圖4中.可以看出,在大約40000步后,所有參數(shù)Markov鏈?zhǔn)諗坑跇O限分布.Gelman-Rubin PSRF診斷得到4個(gè)參數(shù)的值分別為1.0470, 1.0361, 1.0362, 1.0017,接近于1,也說(shuō)明Markov鏈?zhǔn)諗?自相關(guān)函數(shù)(Autocorrelation function, ACF)分析表明當(dāng)遲滯大于30個(gè)迭代步后,ACF較小,殘差均值相關(guān)性弱.高ACF意味著鏈內(nèi)低混合,常會(huì)導(dǎo)致收斂到后驗(yàn)分布較慢.

至此Bayesian推理應(yīng)急溯源技術(shù)完成了反演計(jì)算,給出了源項(xiàng)后驗(yàn)信息概率分布形式的結(jié)果.在應(yīng)急響應(yīng)中,響應(yīng)人員可以在均值正負(fù)一個(gè)標(biāo)準(zhǔn)差的范圍內(nèi),或者貝葉斯區(qū)間的范圍內(nèi)進(jìn)行污染源位置物理搜索和化學(xué)確認(rèn).對(duì)污染物排放總量和時(shí)刻也可進(jìn)行同樣的推理.

圖2 源項(xiàng)參數(shù)和方差的后驗(yàn)概率密度曲線

圖3 源項(xiàng)參數(shù)和方差σ2的迭代軌跡

圖4 源項(xiàng)參數(shù)和方差σ2在MCMC采樣過(guò)程的均值變化

3 Bayesian推理應(yīng)急溯源的影響因素分析

3.1 參數(shù)相關(guān)性分析

源項(xiàng)參數(shù)間的相關(guān)性必然會(huì)對(duì)反演結(jié)果的分布形式造成影響.對(duì)Case-T1案例后驗(yàn)樣本中源項(xiàng)參數(shù)相關(guān)性分析發(fā)現(xiàn)s和s的相關(guān)系數(shù)高達(dá)0.999,說(shuō)明兩者存在依賴關(guān)系.高相關(guān)性容易導(dǎo)致收斂到后驗(yàn)分布變慢,前面研究可以看到Markov鏈經(jīng)過(guò)上萬(wàn)步迭代才收斂.

實(shí)際上,分析無(wú)量綱對(duì)流擴(kuò)散方程(式17)即可發(fā)現(xiàn)兩者確實(shí)存在著依賴關(guān)系.該線性偏微分方程中系數(shù)只有Peclet準(zhǔn)數(shù),,該線性系統(tǒng)結(jié)構(gòu)的差異只依賴于參數(shù).小Peclet準(zhǔn)數(shù)流動(dòng)是擴(kuò)散控制而大Peclet準(zhǔn)數(shù)流動(dòng)則是對(duì)流控制.

只需確定幾何特征尺度,其無(wú)量綱時(shí)間就可確定,定義應(yīng)急監(jiān)測(cè)斷面間的平均距離為. Truckee River 案例中,以10km計(jì)可以計(jì)算出為329,為對(duì)流控制.對(duì)于Case-T2,特征尺度取600m求得為11.7,反演參數(shù)相關(guān)性將要比Case-T1小.事實(shí)上,通過(guò)計(jì)算發(fā)現(xiàn)Case-T2溯源反演源項(xiàng)參數(shù)的后驗(yàn)樣本中s和s相關(guān)系數(shù)為0.8,較Case-T1小很多.

3.2 濃度時(shí)間序列方差假定的影響

首先考察不同誤差項(xiàng)假設(shè)的影響.對(duì)基于異方差假設(shè)似然函數(shù)的方法進(jìn)行計(jì)算,結(jié)果也列于表3中.可以看出除s的反演結(jié)果比均方差假設(shè)的結(jié)果大外,兩者的相差不太大.但是后者的偏度都比較大,正偏斜,說(shuō)明參數(shù)后驗(yàn)PDF的偏度對(duì)誤差假定是敏感的.Bayes區(qū)間也要較同方差假設(shè)要小一些.

可初步推斷,s和s更適合于用同方差假設(shè)而s適合異方差假設(shè),可能也是由混合區(qū)影響所致.對(duì)于s,不同監(jiān)測(cè)斷面距混合區(qū)距離不同,造成實(shí)測(cè)數(shù)據(jù)同正向模型計(jì)算間偏差不同,從而導(dǎo)致異方差的分布形式.這一結(jié)論同樣被小尺度Lagan River案例所證實(shí).需要注意的是,在具體應(yīng)用中異方差假設(shè)需依據(jù)歷史觀測(cè)數(shù)據(jù)去調(diào)節(jié)穩(wěn)定化因子λ1和λ2,這在數(shù)據(jù)稀缺的環(huán)境下存在一定困難.

3.3 后驗(yàn)概率密度采樣方法的影響

考察了傳統(tǒng)的Metropolis采樣和AM采樣方法對(duì)Case-T1和Case-T2反演結(jié)果的影響.其后驗(yàn)參數(shù)的差別不大,但是由于傳統(tǒng)的Metroplis采樣需要對(duì)建議概率密度進(jìn)行測(cè)試,比較麻煩.而AM采樣方法調(diào)節(jié)的參數(shù)少,不需要考慮尋找合適的建議概率密度函數(shù).

3.4 污染物傳輸機(jī)制概化的影響

對(duì)比EPD-RIV1數(shù)值解模型和解析解模型的溯源計(jì)算可知,前者基于有限差分求解線性代數(shù)方程組而后者基于一個(gè)代數(shù)方程,計(jì)算復(fù)雜度差一個(gè)數(shù)量級(jí),所以前者溯源計(jì)算要慢些.EPD- RIV1模型可考慮不同曲線形式的污染源輸入,更適用于排放歷史重構(gòu).

OTIS模型和解析解模型溯源計(jì)算結(jié)果差異不大,主要是因?yàn)門(mén)ruckee River河段不存在死水區(qū).但是由于OTIS模型多了三個(gè)參數(shù),河流斷面面積、死水區(qū)斷面面積和死水區(qū)同主流區(qū)間的傳質(zhì)交換系數(shù),需要更多的事前率定工作,在實(shí)際操作過(guò)程中比較麻煩.

4 Bayesian推理溯源方法和確定性優(yōu)化法的對(duì)比辨析

傳統(tǒng)優(yōu)化算法,如遺傳算法,也能夠快速求解簡(jiǎn)單情景下的污染反演問(wèn)題,其策略如圖5所示.然而優(yōu)化算法和貝葉斯算法進(jìn)行溯源分別對(duì)應(yīng)了頻率主義(Frequentism)和貝葉斯主義(Bayesianism)兩個(gè)方法論,目前尚未有人對(duì)兩者在污染溯源問(wèn)題上結(jié)合案例進(jìn)行對(duì)比辨析,有必要在思路、過(guò)程和結(jié)果上進(jìn)行剖析,有益于應(yīng)用實(shí)踐.

頻率主義和貝葉斯主義的本質(zhì)區(qū)別是對(duì)概率定義的差異.前者認(rèn)為概率表示有限次的重復(fù)測(cè)量,后者認(rèn)為概率是對(duì)事件認(rèn)識(shí)程度的度量.這兩學(xué)派圍繞此問(wèn)題在上世紀(jì)60年代展開(kāi)了長(zhǎng)久的爭(zhēng)論[50-51].關(guān)于頻率學(xué)派和Bayes 學(xué)派方法論上區(qū)別的更多探討和解釋可參考相關(guān)文獻(xiàn)和網(wǎng)絡(luò)資料[50-52].

本文進(jìn)一步采用遺傳優(yōu)化算法對(duì)Case-T1進(jìn)行溯源計(jì)算,重復(fù)運(yùn)行20次統(tǒng)計(jì)其結(jié)果,見(jiàn)表4. 需要注意的是,也可以采用bootstrap或jackknife重采樣進(jìn)行統(tǒng)計(jì).

圖5 基于優(yōu)化法的污染溯源求解策略

表4 基于遺傳算法優(yōu)化求解的Case-T1反演結(jié)果

4.1 求解思路的差異分析

在優(yōu)化方法中,源項(xiàng)參數(shù)是未知固定量,通過(guò)優(yōu)化可以找到這個(gè)值或者說(shuō)可能值.一次優(yōu)化模擬求得一個(gè)最優(yōu)解,就可以看成一次試驗(yàn).運(yùn)行20次對(duì)最優(yōu)解進(jìn)行統(tǒng)計(jì),用中值和標(biāo)準(zhǔn)差來(lái)描述反演結(jié)果就是采用了頻率的觀點(diǎn).在此,觀測(cè)數(shù)據(jù)和正向模型捆綁在一起,參與重復(fù)隨機(jī)試驗(yàn).

對(duì)于本文的Bayesian框架,把源項(xiàng)參數(shù)看成未知隨機(jī)變量而觀測(cè)數(shù)據(jù)是確定的知識(shí)和信息.溯源反演中估計(jì)的是隨機(jī)變量在特定場(chǎng)合下所取的特定值.例如,可進(jìn)一步解釋:“根據(jù)以前對(duì)的了解(先驗(yàn)分布)以及現(xiàn)在觀察結(jié)果(濃度樣本C)推斷:未知s有60%的可能性在-20km到-40km之間,有30%的可能性-60km到-40km之間.”當(dāng)然,如果需要一個(gè)更確定的解釋,可采用后驗(yàn)分布均值(本論文采用該方法)和廣義極大似然估計(jì)(似然函數(shù)最大的位置)等.

此外,Bayesian框架中MCMC只不過(guò)是對(duì)后驗(yàn)分布的采樣,如果可以直接求解似然函數(shù)積分方程就可不用MCMC采樣.這和圖5中直接Monte Carlo法作用完全不一樣.兩個(gè)方法中樣本(濃度數(shù)據(jù))所起的作用不一樣.在Bayes學(xué)派中,樣本的唯一作用是在于它使對(duì)所要估計(jì)參數(shù)()的認(rèn)識(shí)起了變化.

4.2 計(jì)算過(guò)程的差異分析

可以說(shuō),優(yōu)化應(yīng)急溯源是為解決直接蒙特卡洛暴力求解低效耗時(shí)的缺點(diǎn)而建立,而B(niǎo)ayesian推理應(yīng)急溯源是為解決反向不確定性量化困難而建立,兩者的計(jì)算過(guò)程也顯然存在差異.從計(jì)算復(fù)雜度上看,優(yōu)化算法控制步驟通常為正向計(jì)算過(guò)程,而B(niǎo)ayes方法控制步驟為后驗(yàn)采樣.兩者計(jì)算收斂都需要調(diào)節(jié)算法參數(shù)來(lái)獲得.另外,Bayes方法要求應(yīng)用者對(duì)算法參數(shù)有一定了解.

在計(jì)算過(guò)程中,Bayesian推理更多從信息論的角度分析和求解問(wèn)題,而優(yōu)化方法更多是數(shù)學(xué)角度.例如Bayeisan框架中似然函數(shù)設(shè)置時(shí)需要對(duì)誤差項(xiàng)做出假設(shè),不同假設(shè)就附帶了對(duì)給定問(wèn)題的不同理解和分析.對(duì)先驗(yàn)分布所做的假設(shè)也是如此.在Bayesian框架求解過(guò)程中注入了不同信息.而對(duì)后驗(yàn)分布的分析和推理也是提取信息的過(guò)程,優(yōu)化法中只考慮選擇什么樣的正向模型,一旦確定后就不考慮其他可用信息.

4.3 溯源結(jié)果的差異分析

對(duì)Truckee River示蹤實(shí)驗(yàn)案例的溯源結(jié)果進(jìn)行對(duì)比,兩者在均值上表現(xiàn)基本一致,盡管Bayesian方法得到的反演結(jié)果要略接近真實(shí)值.優(yōu)化法沒(méi)有給出置信區(qū)間(也可通過(guò)前述的bootstrap法近似給出),Bayesian推理通過(guò)后驗(yàn)樣本得到Bayes信任區(qū)間為決策者提供決策信心.

在兩個(gè)方法對(duì)標(biāo)準(zhǔn)差的解釋也有差異.優(yōu)化法是20次運(yùn)行中最優(yōu)解統(tǒng)計(jì)得到的標(biāo)準(zhǔn)差,不能認(rèn)為是源項(xiàng)參數(shù)概率分布的描述.Bayesian方法基于MCMC采樣通過(guò)參數(shù)后驗(yàn)分布而得到標(biāo)準(zhǔn)差,源項(xiàng)參數(shù)反向不確定性定量的客觀描述,因此更有意義.從而在確定溯源反演方法后,在真實(shí)污染事件的溯源反演過(guò)程中提高污染物傳輸建模的精度是減少結(jié)構(gòu)不確定性提高反演準(zhǔn)確度的關(guān)鍵之一.

5 結(jié)論

5.1 針對(duì)貝葉斯推理技術(shù)在河流突發(fā)污染溯源實(shí)際應(yīng)用中存在的問(wèn)題,結(jié)合實(shí)地示蹤劑實(shí)驗(yàn)數(shù)據(jù),開(kāi)展算法參數(shù)、影響因素、同頻率法對(duì)比辨析的研究. 針對(duì)傳統(tǒng)MH算法建議密度選擇的問(wèn)題,采用AM算法進(jìn)行后驗(yàn)概論密度采樣;針對(duì)監(jiān)測(cè)污染物濃度時(shí)間序列不同方差情形,建立不同的似然函數(shù)進(jìn)行溯源.三個(gè)示蹤劑實(shí)驗(yàn)案例,基于AM采樣的貝葉斯反演都取得了較好的結(jié)果.驗(yàn)證了貝葉斯方法在實(shí)際污染溯源中的適用性.

5.2 污染傳輸?shù)奈锢頇C(jī)制中源項(xiàng)參數(shù)排放時(shí)間s和排放時(shí)刻s存在固有的相關(guān)性,對(duì)反演準(zhǔn)確度造成影響.源項(xiàng)參數(shù)和方差的后驗(yàn)概率密度的偏度對(duì)方差假定敏感.對(duì)算法使用的目標(biāo)函數(shù),異方差假定中穩(wěn)定化因子λ1、λ2,AM采樣建議比例因子sd等參數(shù)都給出了推薦值.這為貝葉斯推理技術(shù)在污染溯源的實(shí)際應(yīng)用提供了較為重要的參考價(jià)值.

[1] Xue P, Zeng W. Trends of environmental accidents and impact factors in China [J]. Frontiers of Environmental Science & Engineering in China, 2011,5(2):266-276.

[2] 曲建華,孟憲林,尤 宏.兩階段評(píng)估體系篩選水源突發(fā)污染應(yīng)急最優(yōu)技術(shù)方案 [J]. 中國(guó)環(huán)境科學(xué), 2015,35(10):3193-200.

[3] Staff Reporter. Crocodile River water affected by toxic spill [M/OL]. 2005[2017-03-05].http://mg.co.za/article/2005-12-23- crocodile-river-water-affected-by-toxic-spill.

[4] India Environment Portal [M/OL]. 2017[2017-03-05]http://www. indiaenvironmentportal.org.in/category/thesaurus/water-pollution.

[5] NRC. National Response Center (NRC) data download [M/OL]. [2017-03-05] http://www.nrc.uscg.mil/download.html.

[6] 王圣瑞,張 蕊,過(guò)龍根,等.洞庭湖水生態(tài)風(fēng)險(xiǎn)防控技術(shù)體系研究 [J]. 中國(guó)環(huán)境科學(xué), 2017,37(5):1896-905.

[7] 陳媛華,王 鵬,姜繼平.基于相關(guān)系數(shù)優(yōu)化法的河流突發(fā)污染源項(xiàng)識(shí)別 [J]. 中國(guó)環(huán)境科學(xué), 2011,31(11):1802-1807.

[8] 彭亞綿,劉春鳳,楊愛(ài)民.二維對(duì)流一擴(kuò)散方程反問(wèn)題的遺傳算法求解 [J]. 河北理工大學(xué)學(xué)報(bào)(自然科學(xué)版), 2008,30(2):84- 87.

[9] 辛小康,韓小波,李 建.基于遺傳算法的水污染事故污染源識(shí)別模型 [J]. 水電能源科學(xué), 2014,32(7):52-55.

[10] 牟行洋.基于微分進(jìn)化算法的污染物源項(xiàng)識(shí)別反問(wèn)題研究 [J]. 水動(dòng)力學(xué)研究與進(jìn)展A輯, 2011,26(1):24-30.

[11] Boano F, Revelli R, Ridolfi L. Source identification in river pollution problems: A geostatistical approach [J]. Water Resources Research, 2005,41(7):1-13

[12] Cheng W P, Jia Y. Identification of contaminant point source in surface waters based on backward location probability density function method [J]. Advances in Water Resources, 2010,33(4): 397-410.

[13] 王家彪,雷曉輝,廖衛(wèi)紅.基于耦合概率密度方法的河渠突發(fā)水污染溯源 [J]. 水利學(xué)報(bào), 2015,46(11):1280-1289.

[14] 吳自庫(kù),范海梅,陳秀榮.對(duì)流－擴(kuò)散過(guò)程逆過(guò)程反問(wèn)題的伴隨同化研究 [J]. 水動(dòng)力學(xué)研究與進(jìn)展, 2008,23(2):111-115.

[15] Hamdi A. Inverse source problem in a 2D linear evolution transport equation: detection of pollution source [J]. Inverse Problems in Science and Engineering, 2012,20(3):401-421.

[16] Hamdi A. Identification of point sources in two-dimensional advection-diffusion-reaction equation: application to pollution sources in a river. Stationary case [J]. Inverse Problems in Science & Engineering, 2007,15(8):855-870.

[17] Hamdi A. The recovery of a time-dependent point source in a linear transport equation: application to surface water pollution [J]. Inverse Problems, 2009,25(7):6-23.

[18] 吳一亞,金文龍,吳云波.寬淺河道瞬時(shí)源源項(xiàng)反問(wèn)題及反演精度主要影響因子分析 [J]. 水資源保護(hù), 2015,31(5):58-61.

[19] 高 琦,韓龍喜,陳麗娜.平面二維河道瞬時(shí)源反演及反演精度影響分析 [J]. 四川環(huán)境, 2016,35(3):67-72.

[20] Marshall L, Nott D, Sharma A. A comparative study of Markov chain Monte Carlo methods for conceptual rainfall-runoff modeling [J]. Water Resources Research, 2004,40(2):1-11.

[21] Campbell E P, Fox D R, Bates B C. A Bayesian Approach to parameter estimation and pooling in nonlinear flood event models [J]. Water Resources Research, 1999,35(1):211-220.

[22] Bates B C, Campbell E P. A Markov Chain Monte Carlo Scheme for parameter estimation and inference in conceptual rainfall- runoff modeling [J]. Water Resources Research, 2001,37(4):937- 947.

[23] Loos M, Krauss M, Fenner K. Pesticide Nonextractable Residue Formation in Soil: Insights from Inverse Modeling of Degradation Time Series [J]. Environmental Science & Technology, 2012,46(18):9830-9837.

[24] 朱 嵩.基于貝葉斯推理的環(huán)境水力學(xué)反問(wèn)題研究 [D]. 浙江大學(xué), 2008.

[25] 朱 嵩,劉國(guó)華,毛根海.利用貝葉斯推理估計(jì)二維含源對(duì)流擴(kuò)散方程參數(shù) [J]. 四川大學(xué)學(xué)報(bào)(工程科學(xué)版), 2008,40(2): 38-43.

[26] 朱 嵩,劉國(guó)華,王立忠.水動(dòng)力-水質(zhì)耦合模型污染源識(shí)別的貝葉斯方法 [J]. 四川大學(xué)學(xué)報(bào)(工程科學(xué)版), 2009,41(5):30-35.

[27] 毛 星.基于貝葉斯理論的事故場(chǎng)景重建技術(shù) [D]. 天津:南開(kāi)大學(xué), 2009.

[28] 陳海洋,滕彥國(guó),王金生.基于Bayesian-MCMC方法的水體污染識(shí)別反問(wèn)題 [J]. 湖南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2012,39(6):74- 78.

[29] 曹小群,宋君強(qiáng),張衛(wèi)民.對(duì)流-擴(kuò)散方程源項(xiàng)識(shí)別反問(wèn)題的MCMC方法 [J]. 水動(dòng)力學(xué)研究與進(jìn)展A輯, 2010,25(2):127- 136.

[30] Wei G, Chi Z, Yu L, et al. Source identification of sudden contamination based on the parameter uncertainty analysis [J]. Journal of Hydroinformatics, 2016,18(6):919-927.

[31] Thomann R V, Mueller J A. Principal of surface water quality modelling and control [M]. Prentice Hall, 1987.

[32] 謝更新.水環(huán)境中的不確定性理論與方法研究—以三峽水庫(kù)為例 [D]. 長(zhǎng)沙:湖南大學(xué), 2005.

[33] Runkel R L. One-dimensional transport with inflow and storage (OTIS): a solute transport model for streams and rivers [M/OL]. Water-Resource Investigations Report, 1998.

[34] Scales J A, Tenorio L. Prior information and uncertainty in inverse problems [J]. Geophysics, 2001,66(2):389-397.

[35] Ntzoufras I. Bayesian modeling using WinBUGS [M]. Hoboken, New Jersey: John Wiley&Sons, 2009.

[36] Keats A, Yee E, Lien F-S. Bayesian inference for source determination with applications to complex urban environment [J]. Atmospheric Environment, 2007,41(3):465-479.

[37] Keats A, Yee E, Lien F S. Information-driven receptor placement for contaminant source determination [J]. Environmental Modelling & Software, 2010,25(9):1000-1013.

[38] Box G E P, Cox D R. An Analysis of Transformations [J]. Journal of the Royal Statistical Society, 1964,26(2):211-252.

[39] Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. Equation of state calculations by fast computing machines [J]. Journal of Chemical Physics, 1953,21(10):87-92.

[40] Hasting W K. Monte Carlo sampling methods using Markov chains and their applications [J]. Biometrika, 1970,57(1):97-109.

[41] Geman S, Geman D. Stochastic relaxtion, Gibbs distirubtions and the Bayesian restoration of images [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1984,6:721-741.

[42] Liu J S. Monte Carlo Strategies in Scientific Computing [M]. New York: Springer-Verlag, 2001.

[43] Haario H, Saksman E, Tamminen J. An adaptive Metropolis algorithm [J]. Bernoulli, 2001,7(2):223-242.

[44] Haario H, Saksman E, Tamminen J. Componentwise adaptation for high dimensional MCMC [J]. Computational Statistics, 2005,20(2):265-273.

[45] Cowles M K, Carlin B P. Markov chain Monte Carlo convergence diagnostics: a comparative review [J]. Journal of the American Statistical Association, 1996,91(434):883-904.

[46] Rathbun R E, Shultz D J, Stephens D W. Preliminary experiments with a modified tracer technique for measuring stream reaeration coefficients [M/OL]. USGS Open-File Report, 1975. http: //pubs.er.usgs.gov/publication/ofr75256.

[47] Crompton J. Traveltime Data for the Truckee River Between Tahoe City, California, and Vista, Nevada, 2006 and 2007. USGS OFR 2008-1084 [M]. 2008.

[48] Rivord J, Saito L, Miller G, et al. Modeling Contaminant Spills in a Regulated River in the Western United States [J]. Journal of Water Resources Planning and Management, 2014,140(3): 343-354.

[49] Reid S E, Mackinnon P A, Elliot T. Direct measurements of reaeration rates using noble gas tracers in the River Lagan, Northern Ireland [J]. Water and Environment Journal, 2007, 21(3):182-191.

[50] Efron B. Why Isn't Everyone a Bayesian? [J]. The American Statistician, 1986,40(1):1-5.

[51] Lindley D V. The Future of Statistics: A Bayesian 21st Century [J]. Advances in Applied Probability, 1975:7.

[52] Efron B. Bayesian inference and the parametric bootstrap [J]. The Annals of Applied Statistics, 2012,6(4):1971-1997.

致謝：感謝新南威爾士大學(xué)Ashish Sharma教授對(duì)論文的指導(dǎo), 感謝莫納什大學(xué)馬來(lái)西亞分校邱順添教授對(duì)論文英文部分的審閱和修訂.

Applicability of Bayesian inference approach for pollution source identification of river chemical spills: A tracer experiment based analysis of algorithmic parameters, impacts and comparison with Frequentist approaches.

JIANG Ji-ping1,2*, DONG Fu-jia1, LIU Ren-tao1, YUAN Yi-xing1

(1.School of Environment, Harbin Institute of Technology, Harbin 150090, China；2.School of Environmental Science and Engineering, Southern University of Science and Technology, Shenzhen 518055, China)., 2017,37(10)：3813~3825

Based on Bayes theorem and combined variance assumptions on pollutant concentration time series with Adaptive-Metropolis sampling, a modular Bayesian approach was established targeting at pollutant source identification during spills. This probability approach updated the prior knowledge on source information by combining experiments and monitoring and was able to directly characterize uncertainty due to the inversion process by probability distribution. Source inversion test results from field tracer experiments were investigated to determine the validity of this Bayesian inference approach, correlation of posterior parameter and impact factors. Results indicate that Bayesian approach was successful in identifying the source parameters and could effectively reduce the emergency decision risks. It is shown that the skewness of posterior distribution of source parameters and variation were sensitive to assumed variance. Using RMSE as objective function, test results also suggested that the default parameters for the established Bayesian source inversion method, were as follows: heteroscedasticity setting stabilization factors λ1= 0, and λ2= 0.1-0.5, and AM sampling proposal scale factor sd=0.1-0.3. Comparisons between the Bayesian approach and optimization approach on aspects of solution methodology, computing process and inverse results were made and differentiation were highlighted. This work provides valuable references for the practical usage of Bayesian approach in surface water pollution source identification.

Bayesian inference；source inversion；river chemical spill；Adaptive-Metropolis sampling；river tracer experiments

X52

A

1000-6923(2017)10-3813-13

姜繼平(1986-),男,浙江金華人,講師,博士,主要從事環(huán)境數(shù)學(xué)模型與決策支持系統(tǒng)方向研究.發(fā)表論文30余篇.

2017-03-07

國(guó)家水體污染控制與治理科技重大專項(xiàng)基金(2012ZX07205-005);中國(guó)博士后科學(xué)基金(2014M551249)

* 責(zé)任作者, 講師, jjp_lab@sina.com