基于河流示蹤實驗的 Bayes 污染溯源:算法參數、影響因素及頻率法對比

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

?

基于河流示蹤實驗的 Bayes 污染溯源:算法參數、影響因素及頻率法對比

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

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

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

貝葉斯推理；污染源反演；突發水污染；AM采樣；河流示蹤劑試驗

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

污染溯源技術的理論基礎是求解特定情景下的源項反演問題,求解過程融入污染物輸移模型.理論上如果對污染輸移模型采用負時間步和空間步進行數值計算,應該能求得污染輸移歷史過程,然而在現實中由于監測誤差和數值誤差的存在,其不可控的誤差反向傳播導致了反算結果極大地偏離實際值.該特征就是反問題不適定性的典型表現,也正因為此,在各領域中出現了形形色色的反演技術來求解反問題.

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

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

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

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

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

1 基于Bayesian推理的不確定性應急溯源方法構建

1.1 正向模型

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

在監測點(,)處的污染物濃度可由下式計算:

式中:為河流斷面面積,m2;D為縱向平均擴散系數,m2/min;為河流平均流速,m;為衰減系數,min-1;為監測時刻,min;為監測點位置,m.

考慮到污染溯源研究中源位置未知,未將坐標原點建立在源排放處和排放時刻.慮到研究問題的特征時間尺度,將時間單位設定為min或者h.

同時對于連續排放源的情境,也可推導出解析解,具體形式可參考文獻[31].為對比不同污染傳輸機制概化的影響,還采用了EPD-RIV1模型[32]和OTIS模型[33]兩個基于數值解的一維水質模型.

1.2 應急溯源Bayes理論框架的建立

1.2.1 污染源溯源問題的廣義陳述 在河流某斷面檢出污染物嚴重超標,在該斷面上游開展應急監測,獲得個時空采樣點的濃度數據,它們和污染物傳輸機理模型間通過正向模型算子建立聯系(式2).該算子即為將模型映射到數據空間的函數.實際上,正向模型算子通常只能是真實物理過程的近似,當使用函數來描述該物理過程時存在系統偏差,可用維向量model來表示.同樣也存在一個維向量meas,表示隨機測量誤差.

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

1.2.2 面向應急溯源的Bayes理論框架 基于污染溯源問題的Bayes公式如下:

式中:()是源項參數的后驗分布,(|)是先驗分布,()是似然函數,()是證據,為源項,為濃度監測數據,為背景信息,即用于確定先驗分布的信息.

1.2.3 Bayes推理應急溯源方法的建立 基于模型的Bayes推理框架常分為四步:模型構建、后驗分布計算、分析后驗分布和決策推理[35].

1.2.3.1 第一步:模型構建

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

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

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

對于異方差非相關誤差項,可通過對監測數據進行如下轉換而穩定化.Box and Cox[38]證實了它在水文數據中的適用性.

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

基于異方差非相關誤差項假定的似然函數可以寫為:

式中:是X的方差.

(2)設置先驗概率

用Bayes方法處理統計推理問題的前提是給出參數的先驗分布.在文獻中提及的方法有客觀法、主觀概率法、同等無知原則法、共軛先驗分布等.在水文學中,常用客觀法即經驗貝葉斯法,基于大量歷史觀測數據,如徑流量等,得到先驗分布,常用的有Gamma分布,正態分布,均一分布等.

河流污染應急溯源中的先驗信息常常有限,因此選擇均一分布作為優先考慮的先驗概率密度函數.當然,如果河流某區段存在大量風險源,如化工區、養殖場等,可優先考慮該污染源出現在該區段的可能性,從而設置綜合的概率密度函數.

均一分布的先驗概率密度可如下描述:

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

(3)后驗概率密度函數

基于同方差非相關誤差項假設的后驗概率密度函數可由如下公式描述:

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

1.2.3.2 第二步:計算后驗分布-MCMC采樣

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

Marshall等[20]經過一個全面的MCMC算法對比實驗后,建議在降雨-徑流機理模型中使用Adaptive Metropolis(AM)算法.本文也采用AM算法作為采樣方法.

1.2.3.3 第三步:分析后驗概率密度

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

1.2.3.4 第四步:對結果進行推斷

依據上述統計分析結果,對污染源信息做最后推斷.本文采用中值和貝葉斯區間對源參數結果進行推斷.此外,在實際應急響應過程中,還需要結合物理方法進行實地驗證,才能最后確定污染源信息.

1.3 后驗概率密度采樣的AM算法

1.3.1 Markov鏈及其收斂 Markov鏈是一個未來狀態只和當前狀態有關而與過去狀態相獨立的隨機過程.數學上的定義及Markov鏈收斂到穩態分布的嚴格數學分析和證明可以參見文獻[42].

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

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

建議協方差的計算需要定義一個任意的初始協方差0.為讓該過程自動進行,設置這個初始協方差為先驗分布條件下的參數的初始協方差,例如最初的5%迭代步長的參數的協方差.

綜上,實現AM算法的步驟如下:

(1)初始化=0

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

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

c. 計算接受率:

式中:(|) 是似然函數,()是的先驗分布.

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

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

(3)重復步驟(2).

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

1.3.3 Markov鏈的收斂診斷 診斷Markov鏈收斂方法已經很多[45],其中最為常用是Gelman 和Rubin 1992年開發以及Raftery和Lewis于同年開發的統計方法.本文采用Gelman-Rubin潛在規模減縮因子(potential scale reduction factor, PSRF)來診斷模型的收斂.R的計算如下所示:

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

2 基于實地示蹤實驗的污染溯源

2.1 示蹤劑實驗

結合論文研究目標,尋找到USGS開展的三個典型河流示蹤劑實驗用于污染溯源研究,重新構建溯源情景并整理其監測數據.其中以Truckee River示蹤劑實驗為主,用于算法驗證; Lagan River 和West Hobolochitto Creek示蹤劑實驗為輔,用于對比和影響因素分析.這三個案例有著不同的河道形態與水動力特征及污染物排放特征,案例分別標識為Case-T1,T2和T3.

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

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

2.1.2 Lagan River示蹤劑實驗 2004年11月與2005年4月間,在北愛爾蘭的Lagan River進行了8次示蹤劑投放實驗,用來測量河流的復氧系數[49],試驗河段長2.6km.2004年12月10號進行的B號測試監測數據用于溯源反演.惰性氣體氙、氪在6個大氣壓下通入水中耗時30s,31.8min后,在下游1200m 和2600m 處的站點1和2監測穿透曲線,采樣每隔數min進行1次.更多細節見文獻[49].

2.1.3 West Hobolochitto Creek示蹤劑實驗 Rathbun于1975年7月在美國密西西比州West Hobolochitto Creek進行了多次氣體示蹤劑實驗來測量大氣復氧速率系數[46].West Hobolochitto Creek為一條水文特征復雜既有死水區又有激流的小溪.研究河段長4190m,實驗時河流處在低流量條件.

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

2.2 正向模型參數化

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

表1 Truckee River案例源項參數和正向模型設置

2.3 采樣算法參數設置

表2 Truckee River案Bayesian溯源中AM采樣參數設置

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

2.4 溯源反演結果

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

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

表3 Case-T1(Truckee River)反演結果的概要統計量

圖1 監測濃度和正向預測濃度值對比

PDF曲線通過高斯核對Markov鏈中后80,000個樣本進行估計得到.PDF曲線形態也定性佐證了上述分析結論.

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

至此Bayesian推理應急溯源技術完成了反演計算,給出了源項后驗信息概率分布形式的結果.在應急響應中,響應人員可以在均值正負一個標準差的范圍內,或者貝葉斯區間的范圍內進行污染源位置物理搜索和化學確認.對污染物排放總量和時刻也可進行同樣的推理.

圖2 源項參數和方差的后驗概率密度曲線

圖3 源項參數和方差σ2的迭代軌跡

圖4 源項參數和方差σ2在MCMC采樣過程的均值變化

3 Bayesian推理應急溯源的影響因素分析

3.1 參數相關性分析

源項參數間的相關性必然會對反演結果的分布形式造成影響.對Case-T1案例后驗樣本中源項參數相關性分析發現s和s的相關系數高達0.999,說明兩者存在依賴關系.高相關性容易導致收斂到后驗分布變慢,前面研究可以看到Markov鏈經過上萬步迭代才收斂.

實際上,分析無量綱對流擴散方程(式17)即可發現兩者確實存在著依賴關系.該線性偏微分方程中系數只有Peclet準數,,該線性系統結構的差異只依賴于參數.小Peclet準數流動是擴散控制而大Peclet準數流動則是對流控制.

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

3.2 濃度時間序列方差假定的影響

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

可初步推斷,s和s更適合于用同方差假設而s適合異方差假設,可能也是由混合區影響所致.對于s,不同監測斷面距混合區距離不同,造成實測數據同正向模型計算間偏差不同,從而導致異方差的分布形式.這一結論同樣被小尺度Lagan River案例所證實.需要注意的是,在具體應用中異方差假設需依據歷史觀測數據去調節穩定化因子λ1和λ2,這在數據稀缺的環境下存在一定困難.

3.3 后驗概率密度采樣方法的影響

考察了傳統的Metropolis采樣和AM采樣方法對Case-T1和Case-T2反演結果的影響.其后驗參數的差別不大,但是由于傳統的Metroplis采樣需要對建議概率密度進行測試,比較麻煩.而AM采樣方法調節的參數少,不需要考慮尋找合適的建議概率密度函數.

3.4 污染物傳輸機制概化的影響

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

OTIS模型和解析解模型溯源計算結果差異不大,主要是因為Truckee River河段不存在死水區.但是由于OTIS模型多了三個參數,河流斷面面積、死水區斷面面積和死水區同主流區間的傳質交換系數,需要更多的事前率定工作,在實際操作過程中比較麻煩.

4 Bayesian推理溯源方法和確定性優化法的對比辨析

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

頻率主義和貝葉斯主義的本質區別是對概率定義的差異.前者認為概率表示有限次的重復測量,后者認為概率是對事件認識程度的度量.這兩學派圍繞此問題在上世紀60年代展開了長久的爭論[50-51].關于頻率學派和Bayes 學派方法論上區別的更多探討和解釋可參考相關文獻和網絡資料[50-52].

本文進一步采用遺傳優化算法對Case-T1進行溯源計算,重復運行20次統計其結果,見表4. 需要注意的是,也可以采用bootstrap或jackknife重采樣進行統計.

圖5 基于優化法的污染溯源求解策略

表4 基于遺傳算法優化求解的Case-T1反演結果

4.1 求解思路的差異分析

在優化方法中,源項參數是未知固定量,通過優化可以找到這個值或者說可能值.一次優化模擬求得一個最優解,就可以看成一次試驗.運行20次對最優解進行統計,用中值和標準差來描述反演結果就是采用了頻率的觀點.在此,觀測數據和正向模型捆綁在一起,參與重復隨機試驗.

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

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

4.2 計算過程的差異分析

可以說,優化應急溯源是為解決直接蒙特卡洛暴力求解低效耗時的缺點而建立,而Bayesian推理應急溯源是為解決反向不確定性量化困難而建立,兩者的計算過程也顯然存在差異.從計算復雜度上看,優化算法控制步驟通常為正向計算過程,而Bayes方法控制步驟為后驗采樣.兩者計算收斂都需要調節算法參數來獲得.另外,Bayes方法要求應用者對算法參數有一定了解.

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

4.3 溯源結果的差異分析

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

在兩個方法對標準差的解釋也有差異.優化法是20次運行中最優解統計得到的標準差,不能認為是源項參數概率分布的描述.Bayesian方法基于MCMC采樣通過參數后驗分布而得到標準差,源項參數反向不確定性定量的客觀描述,因此更有意義.從而在確定溯源反演方法后,在真實污染事件的溯源反演過程中提高污染物傳輸建模的精度是減少結構不確定性提高反演準確度的關鍵之一.

5 結論

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

5.2 污染傳輸的物理機制中源項參數排放時間s和排放時刻s存在固有的相關性,對反演準確度造成影響.源項參數和方差的后驗概率密度的偏度對方差假定敏感.對算法使用的目標函數,異方差假定中穩定化因子λ1、λ2,AM采樣建議比例因子sd等參數都給出了推薦值.這為貝葉斯推理技術在污染溯源的實際應用提供了較為重要的參考價值.

[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] 曲建華,孟憲林,尤 宏.兩階段評估體系篩選水源突發污染應急最優技術方案 [J]. 中國環境科學, 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] 王圣瑞,張 蕊,過龍根,等.洞庭湖水生態風險防控技術體系研究 [J]. 中國環境科學, 2017,37(5):1896-905.

[7] 陳媛華,王 鵬,姜繼平.基于相關系數優化法的河流突發污染源項識別 [J]. 中國環境科學, 2011,31(11):1802-1807.

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

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

[10] 牟行洋.基于微分進化算法的污染物源項識別反問題研究 [J]. 水動力學研究與進展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] 王家彪,雷曉輝,廖衛紅.基于耦合概率密度方法的河渠突發水污染溯源 [J]. 水利學報, 2015,46(11):1280-1289.

[14] 吳自庫,范海梅,陳秀榮.對流－擴散過程逆過程反問題的伴隨同化研究 [J]. 水動力學研究與進展, 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] 吳一亞,金文龍,吳云波.寬淺河道瞬時源源項反問題及反演精度主要影響因子分析 [J]. 水資源保護, 2015,31(5):58-61.

[19] 高 琦,韓龍喜,陳麗娜.平面二維河道瞬時源反演及反演精度影響分析 [J]. 四川環境, 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] 朱 嵩.基于貝葉斯推理的環境水力學反問題研究 [D]. 浙江大學, 2008.

[25] 朱 嵩,劉國華,毛根海.利用貝葉斯推理估計二維含源對流擴散方程參數 [J]. 四川大學學報(工程科學版), 2008,40(2): 38-43.

[26] 朱 嵩,劉國華,王立忠.水動力-水質耦合模型污染源識別的貝葉斯方法 [J]. 四川大學學報(工程科學版), 2009,41(5):30-35.

[27] 毛 星.基于貝葉斯理論的事故場景重建技術 [D]. 天津:南開大學, 2009.

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

[29] 曹小群,宋君強,張衛民.對流-擴散方程源項識別反問題的MCMC方法 [J]. 水動力學研究與進展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] 謝更新.水環境中的不確定性理論與方法研究—以三峽水庫為例 [D]. 長沙:湖南大學, 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.

致謝：感謝新南威爾士大學Ashish Sharma教授對論文的指導, 感謝莫納什大學馬來西亞分校邱順添教授對論文英文部分的審閱和修訂.

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-),男,浙江金華人,講師,博士,主要從事環境數學模型與決策支持系統方向研究.發表論文30余篇.

2017-03-07

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

* 責任作者, 講師, jjp_lab@sina.com