張雙圣，強(qiáng) 靜，劉漢湖，劉喜坤，孫韶華

(1.中國(guó)礦業(yè)大學(xué) 環(huán)境與測(cè)繪學(xué)院, 江蘇 徐州 221116； 2.中國(guó)礦業(yè)大學(xué) 數(shù)學(xué)學(xué)院, 江蘇 徐州 221116； 3.徐州市城區(qū)水資源管理處,江蘇 徐州 221018； 4.山東省城市供排水水質(zhì)監(jiān)測(cè)中心,山東 濟(jì)南 250100)

0 引言

目前,國(guó)內(nèi)外水污染溯源問(wèn)題的研究較多,其本質(zhì)是水質(zhì)模型參數(shù)的識(shí)別。 參數(shù)識(shí)別方法主要有馬爾科夫鏈蒙特卡羅算法(Markov chain Monte Carlo,MCMC)[1-7]、微分進(jìn)化算法[8-11]、遺傳算法[12-13]和最小二乘算法[14]等。陳海洋等[3]運(yùn)用經(jīng)典MCMC算法實(shí)現(xiàn)了水體污染源項(xiàng)的識(shí)別,得到了污染源強(qiáng)度、污染源位置和污染源排放時(shí)長(zhǎng)3個(gè)未知參數(shù)的估計(jì)值。閔濤等[11]利用遺傳算法求解對(duì)流-擴(kuò)散方程源項(xiàng)識(shí)別反問(wèn)題。這些算法可分為確定性方法和不確定性方法,其中確定性方法未考慮誤差因素,得到參數(shù)確定的估計(jì)值,但容易丟掉“真值”；不確定性方法具有較強(qiáng)的隨機(jī)性,收斂速度慢,計(jì)算量隨著參數(shù)的增多而呈指數(shù)增長(zhǎng),得到待估參數(shù)的一個(gè)小的可靠范圍。其中,MCMC方法應(yīng)用較為廣泛。Metropolis算法[13-14]是一種經(jīng)典的MCMC方法,應(yīng)用非常廣泛。大多數(shù)MCMC算法都是通過(guò)改進(jìn)Metropolis算法得到的,比如延遲拒絕算法(delayed rejection, DR)[15-16]、自適應(yīng)Metropolis算法(AM)[17]、延遲拒絕自適應(yīng)Metropolis算法(DRAM)[18]等。DRAM算法是一種高效自適應(yīng)MCMC算法,其本質(zhì)是把DR算法和AM算法二者組合起來(lái)。DR算法保證了馬爾科夫鏈的局部自適應(yīng),AM算法采用全局自適應(yīng)調(diào)整策略,因此DRAM算法能夠?qū)CMC鏈條進(jìn)行全局和局部的自適應(yīng),DRAM算法抽樣效率優(yōu)于DR算法和AM算法各自單獨(dú)使用的情況[19]。但是Metropolis算法和DRAM算法都只有一條鏈,容易陷入到局部最優(yōu)值[19-20],計(jì)算結(jié)果給決策者造成一定的負(fù)面影響。

運(yùn)用二維地下水水質(zhì)對(duì)流-擴(kuò)散方程,通過(guò)貝葉斯公式建立水體污染源識(shí)別模型,得到污染源強(qiáng)度、排放位置(x,y)和排放時(shí)長(zhǎng)4個(gè)未知參數(shù)的后驗(yàn)概率密度函數(shù),并運(yùn)用拉丁超立方抽樣方法[21]優(yōu)化DRAM算法,實(shí)現(xiàn)參數(shù)反演過(guò)程中準(zhǔn)確性與效率的雙提高,為突發(fā)性水污染反問(wèn)題研究提供借鑒。

1 研究方法

1.1 貝葉斯公式

貝葉斯公式[22-23]如下：

(1)

式中：α為模型的未知參數(shù)；d為實(shí)測(cè)數(shù)據(jù)；p(α|d)為參數(shù)的后驗(yàn)概率密度函數(shù)；p(α)為參數(shù)的先驗(yàn)概率密度函數(shù)；p(d|α)為條件概率密度函數(shù)；p(d)為歸一化的積分常數(shù)。

假設(shè)模型中未知參數(shù)共有m個(gè),則α=(α1,α2,…,αm)。環(huán)境水力學(xué)的模型參數(shù)都分布在一個(gè)特定的范圍內(nèi),可以認(rèn)為每個(gè)參數(shù)都服從均勻分布,且α1,α2,…,αm相互獨(dú)立。模型參數(shù)αi的先驗(yàn)概率密度函數(shù)可定義為：

(2)

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

(3)

假設(shè)模型中實(shí)測(cè)值共有n個(gè),即d=(d1,d2,…,dn)。di表示第i個(gè)實(shí)測(cè)數(shù)據(jù),Ci(x,y,t|α)表示相應(yīng)的第i個(gè)預(yù)測(cè)值,則εi=di-Ci(x,y,t|α)為測(cè)量誤差,i=1,2,…,n。假設(shè)測(cè)量誤差服從均值為0、標(biāo)準(zhǔn)偏差為σ的正態(tài)分布,且相互獨(dú)立,則條件概率密度函數(shù)p(d|α)可表示如下：

(4)

聯(lián)合式(1)～(4),可得后驗(yàn)概率密度函數(shù)p(α|d)為：

(5)

(6)

在實(shí)測(cè)數(shù)據(jù)d固定的條件下,式(6)是關(guān)于參數(shù)α的函數(shù)。式(6)通過(guò)積分求解參數(shù)α的概率分布,計(jì)算較為復(fù)雜,難以得出明確的解析表達(dá)式。并且隨未知參量維數(shù)的增加,數(shù)值積分算法的計(jì)算量將呈指數(shù)增長(zhǎng),實(shí)現(xiàn)復(fù)雜且難度較大。

1.2 延遲拒絕自適應(yīng)Metropolis算法

式(6)一般可采用馬爾科夫鏈蒙特卡羅方法(MCMC)近似求解。MCMC方法的核心是Monte Carlo模擬方法和Markov Chain抽樣方法。當(dāng)“概率事件”樣本點(diǎn)足夠多時(shí),頻率可以近似代表概率,這是Monte Carlo模擬方法的本質(zhì),而Markov Chain抽樣方法能夠保證Markov鏈花更多的時(shí)間在概率大的區(qū)域,節(jié)省Monte Carlo模擬方法的工作量。

Metropolis算法是一種經(jīng)典的MCMC算法,應(yīng)用非常廣泛,具體算法參見(jiàn)文獻(xiàn)[13-14],其提議分布一旦設(shè)定便保持不變,當(dāng)提議分布遠(yuǎn)離目標(biāo)分布時(shí),馬爾科夫鏈?zhǔn)諗克俣葧?huì)變慢。為了提高馬爾科夫鏈?zhǔn)諗啃?Haario等[18]將DR算法和AM算法結(jié)合提出了高效自適應(yīng)DRAM算法,其具體步驟如下。

1.2.1 非自適應(yīng)階段

(1)在模型未知參數(shù)α先驗(yàn)范圍內(nèi)隨機(jī)產(chǎn)生初始的參數(shù)樣本Xt(t=1),非自適應(yīng)階段迭代次數(shù)為N0,提議分布為高斯分布qt=N(Xt,C0),其中C0為協(xié)方差矩陣,要求C0滿足對(duì)稱正定條件,并且高斯分布滿足對(duì)稱隨機(jī)游走。

(3)重復(fù)過(guò)程(2),直至達(dá)到迭代次數(shù)N0。

1.2.2 自適應(yīng)階段

(3)重復(fù)過(guò)程(2),直至達(dá)到迭代次數(shù)N。

1.3 拉丁超立方抽樣

采用Metropolis算法和DRAM算法求解多維參數(shù)反演問(wèn)題,Markov Chain容易受到樣本初始點(diǎn)的影響,可能在局部最優(yōu)處達(dá)到穩(wěn)定或者產(chǎn)生難以收斂的問(wèn)題。

為了保證樣本初始點(diǎn)的隨機(jī)性和均勻性,本研究基于拉丁超立方抽樣方法[21]對(duì)DRAM算法進(jìn)行改進(jìn)。拉丁超立方抽樣方法是一種多維分層隨機(jī)抽樣方法,具有良好的散布均勻性和代表性。假設(shè)要在m維模型參數(shù)α先驗(yàn)范圍[Ai,Bi](i=1,2,…,m)內(nèi)抽取q組樣本,具體步驟如下：

(1)把m維模型參數(shù)的m個(gè)先驗(yàn)范圍[Ai,Bi](i=1,2,…,m)都均分成q個(gè)小區(qū)間,小區(qū)間可記為[Aij,Bij](i=1,2,…,m;j=1,2,…,q),總共產(chǎn)生m×q個(gè)小區(qū)間。

(2)在任一小區(qū)間[Aij,Bij]中隨機(jī)抽取一個(gè)數(shù),記為αij,總共產(chǎn)生m×q個(gè)數(shù),組成矩陣為

(3)對(duì)矩陣Ψmq中行向量[αi1,αi2,…,αiq](i=1,2,…,m)進(jìn)行隨機(jī)排序，所得向量記為[βi1,βi2,…,βiq](i=1,2,…,m);由此得到矩陣

(4) 矩陣Φmq中每個(gè)列向量是一組樣本,共抽取得到q組樣本。

1.4 改進(jìn)型多鏈DRAM算法

改進(jìn)型多鏈DRAM算法具體步驟：

(1)運(yùn)用拉丁超立方抽樣方法,在模型參數(shù)先驗(yàn)范圍內(nèi)隨機(jī)抽取q組初始樣本。

(2)分別以q組初始樣本作為初始點(diǎn)采用DRAM算法進(jìn)行迭代運(yùn)算,得到q條Markov Chains。

(3)對(duì)上面q條Markov Chains求平均值作為最終結(jié)果。

2 算例應(yīng)用

2.1 算例概述

均質(zhì)各向同性地下含水層中,二維地下水水質(zhì)對(duì)流-擴(kuò)散方程[24-25]可以表述為：

(7)

式中:C為測(cè)點(diǎn)(x,y)在t時(shí)長(zhǎng)的污染物的質(zhì)量濃度,mg/L；t為污染物排放后開(kāi)始預(yù)測(cè)的時(shí)長(zhǎng),d；Dx、Dy分別為縱向、橫向的彌散系數(shù), m2/d；u為地下含水層水流平均流速, m/d。

對(duì)于污染物瞬時(shí)點(diǎn)源排放模式,方程的解析解可表示為：

(8)

式中:M為污染物的排放量,g；h為含水層的厚度,m；x、y分別為預(yù)測(cè)點(diǎn)距排放點(diǎn)的縱向、橫向距離,m；k為含水層的有效孔隙率。

假設(shè)研究區(qū)域?yàn)榫|(zhì)各向同性地下含水層(建立坐標(biāo)系，如圖1)。含水層厚度為h=1 m, 孔隙率k=0.3，x軸代表含水層中水流的方向, 水流速度為u=5 m/d；縱向彌散系數(shù)為1.5 m2/d, 橫向彌散系數(shù)為0.3 m2/d。待求參數(shù)真值:污染物泄漏點(diǎn)坐標(biāo)為(X0=200 m,Y0=200 m), 污染源強(qiáng)度M=1 000 g, 污染源排放時(shí)長(zhǎng)T0=110 d。觀測(cè)點(diǎn)D坐標(biāo)為(X1=800 m,Y1=215 m)。

圖1 算例模型示意圖Figure 1 Sketch of example model

2.2 監(jiān)測(cè)數(shù)據(jù)

以構(gòu)建的算例模型的計(jì)算值作為觀測(cè)值, 設(shè)定2017年5月1日上午10∶00點(diǎn)作為初始監(jiān)測(cè)時(shí)間,其后每隔2 d測(cè)量一次,直至2017年5月21日上午10∶00點(diǎn)結(jié)束,將研究區(qū)水文地質(zhì)參數(shù)代入式(8),可計(jì)算得出觀測(cè)點(diǎn)D的污染物濃度di,以此作為觀測(cè)點(diǎn)D的污染物的觀測(cè)值,如表1所示。

表1 觀測(cè)點(diǎn)D處污染物觀測(cè)濃度序列Table 1 Sequence value of pollutant concentration at observation point D

2.3 結(jié)果及比較

根據(jù)先驗(yàn)信息,污染源強(qiáng)度M、污染源位置坐標(biāo)(X0,Y0)和污染源排放時(shí)長(zhǎng)T04個(gè)待求參數(shù)的取值范圍分別為：400 g≤M≤1 600 g, -200 m≤X0≤600 m, 185 m≤Y0≤215 m, 60 d≤T0≤160 d。運(yùn)用MATLAB軟件編程,對(duì)比分析不同算法的求解結(jié)果。

2.3.1 Metropolis算法求解結(jié)果

采用經(jīng)典Metropolis算法構(gòu)建馬爾科夫鏈對(duì)模型進(jìn)行求解,馬爾科夫鏈選取初始點(diǎn)不同,4個(gè)待估參數(shù)M、X0、Y0、T0的迭代曲線差別很大。以待估參數(shù)M為例,隨機(jī)選取兩組不同初始樣本,如樣本點(diǎn)1:(M,X0,Y0,T0)=(827.44,592.05,190.59,80.77)；樣本點(diǎn)2:(M,X0,Y0,T0)=(687.50,141.15,213.89,146.97),待估參數(shù)M迭代曲線見(jiàn)圖2。

圖2 基于Metropolis算法的不同初始點(diǎn)模型參數(shù)M迭代曲線圖Figure 2 Iterative curves of model parameter M based on Metropolis algorithm with different initial values

由圖2(a)可知,在初始樣本點(diǎn)為M1=827.44 g Metropolis算法迭代到5 000次時(shí),參數(shù)M值逐漸趨于穩(wěn)定,并最終在M=1 450 g處收斂,但是該值與真值(M=1 000 g)相差較大,誤差達(dá)45%。另由圖2(b)可知,在樣本初始點(diǎn)為M2=687.50 g,Metropolis算法迭代到50 000次時(shí),參數(shù)M值仍未收斂。因此運(yùn)用Metropolis算法求解污染源強(qiáng)度時(shí),求解結(jié)果受樣本初始點(diǎn)的影響較大,容易出現(xiàn)局部最優(yōu)或者難以收斂的問(wèn)題。同理,在求解另外3個(gè)待求參數(shù)X0、Y0、T0時(shí)也會(huì)出現(xiàn)類(lèi)似問(wèn)題。

2.3.2 多鏈Metropolis算法求解結(jié)果

為防止反演結(jié)果的局部最優(yōu),構(gòu)建基于拉丁超立方抽樣的多鏈Metropolis算法求解算例(設(shè)定40條鏈, 單鏈長(zhǎng)度為50 000)。污染源強(qiáng)度M、污染源位置坐標(biāo)(X0,Y0)和污染源排放時(shí)長(zhǎng)T0的反演迭代曲線見(jiàn)圖3。

圖3 基于多鏈Metropolis算法的模型參數(shù)迭代曲線Figure 3 Iterative curves of model parameters based on multi-chain Metropolis algorithm

由圖3可知, 當(dāng)基于拉丁超立方抽樣的多鏈Metropolis算法迭代到25 000次時(shí),4個(gè)參數(shù)反演結(jié)果均基本在真值處達(dá)到穩(wěn)定且收斂,而且通過(guò)程序驗(yàn)證迭代結(jié)果不受樣本初始點(diǎn)選取的影響,有效解決了經(jīng)典Metropolis算法受初始點(diǎn)的影響導(dǎo)致不收斂或局部最優(yōu)的問(wèn)題,表明該算法具有較強(qiáng)的穩(wěn)定性。剔除前25 000次不穩(wěn)定結(jié)果,對(duì)剩余的25 000次計(jì)算結(jié)果進(jìn)行后驗(yàn)統(tǒng)計(jì)分析,見(jiàn)表2。

由表2可知,基于拉丁超立方抽樣的多鏈Metropolis算法得出的4個(gè)待估參數(shù)M、X0、Y0、T0的均值與真值的誤差分別為4.36%、0.18%、0.35%、0.08%,中值誤差分別為4.30%、0.15%、0.35%、0.06%,表明基于拉丁超立方抽樣方法的Metropolis算法可實(shí)現(xiàn)計(jì)算結(jié)果的全局最優(yōu),有效提高反演結(jié)果的準(zhǔn)確性。

表2 基于多鏈Metropolis算法的模型參數(shù)后驗(yàn)統(tǒng)計(jì)結(jié)果Table 2 Posterior statistical results of model parameters based on multi-chain Metropolis algorithm

2.3.3 改進(jìn)型多鏈DRAM算法求解結(jié)果

多鏈Metropolis算法雖然有效提高了反演結(jié)果的準(zhǔn)確性,但是迭代次數(shù)需達(dá)到25 000次時(shí)才能收斂,反演效率相對(duì)偏低,為此筆者構(gòu)建基于拉丁超立方抽樣的多鏈DRAM算法(設(shè)定40條鏈，單鏈長(zhǎng)度為12 000)進(jìn)行參數(shù)求解。 4個(gè)參數(shù)的反演遍歷均值圖見(jiàn)圖4。

由圖4可知,改進(jìn)型多鏈DRAM算法迭代到8 000次時(shí),4個(gè)參數(shù)取值均基本達(dá)到穩(wěn)定,與基于拉丁超立方抽樣的多鏈Metropolis算法(25 000次時(shí)迭代穩(wěn)定且收斂)相比,速率提高68%,并且通過(guò)程序驗(yàn)證迭代結(jié)果不受樣本初始點(diǎn)的影響,表明基于拉丁超立方抽樣的DRAM方法具有較高的反演效率及穩(wěn)定性。

圖4 基于改進(jìn)型多鏈DRAM算法的模型參數(shù)遍歷均值圖Figure 4 Ergodic mean plots of model parameters based on improved multi-chain DRAM algorithm

對(duì)于改進(jìn)型多鏈DRAM算法的反演結(jié)果,剔除前10 000次不穩(wěn)定結(jié)果,對(duì)剩余的2 000次模型參數(shù)進(jìn)行后驗(yàn)統(tǒng)計(jì)分析,見(jiàn)表3。由表3可知,該算法在提高反演效率的前提下,同樣可保證反演結(jié)果的準(zhǔn)確性。

表3 基于改進(jìn)型多鏈DRAM算法的模型參數(shù)后驗(yàn)統(tǒng)計(jì)結(jié)果Table 3 Posterior statistical results of model parameters based on improved multi-chain DRAM algorithm

3 結(jié)論

(1)利用Bayesian公式以概率語(yǔ)言解決突發(fā)性地下水污染反問(wèn)題,可有效獲取污染源信息(污染源強(qiáng)度、污染源位置和污染源排放時(shí)長(zhǎng))。

(2)運(yùn)用經(jīng)典Metropolis算法求解污染源信息時(shí),求解結(jié)果受樣本初始點(diǎn)的影響較大,容易出現(xiàn)局部最優(yōu)或者難以收斂的問(wèn)題,反演得到的未知參數(shù)估計(jì)值誤差較大。

(3)基于拉丁超立方抽樣的Metropolis算法可實(shí)現(xiàn)反演結(jié)果的全局最優(yōu),有效提高反演結(jié)果的準(zhǔn)確性,但是效率相對(duì)低下。 基于拉丁超立方抽樣的DRAM算法在保證反演結(jié)果準(zhǔn)確性的前提下,可顯著提高反演效率,其可靠性和穩(wěn)定性均優(yōu)于經(jīng)典Metropolis算法。