卡爾曼濾波方法在地下水污染源反演中的應(yīng)用

白玉堃,盧文喜*,李久輝

卡爾曼濾波方法在地下水污染源反演中的應(yīng)用

白玉堃1,2,盧文喜1,2*,李久輝1,2

(1.吉林大學(xué)地下水與資源環(huán)境教育部重點實驗室,吉林 長春 130012；2.吉林大學(xué)新能源與環(huán)境學(xué)院,吉林 長春 130012)

采用卡爾曼濾波方法反演識別地下水污染源的個數(shù)和大概位置.借助一個假想算例,建立地下水系統(tǒng)水流和溶質(zhì)運移模擬模型,利用靈敏度分析方法篩選出對模擬結(jié)果影響最大的參數(shù)作為隨機變量,對該參數(shù)進行抽樣,運用蒙特卡羅方法將抽樣結(jié)果輸入模擬模型,生成污染質(zhì)濃度場.采用卡爾曼濾波方法構(gòu)造迭代過程,逐個利用采樣點處濃度的實測值不斷更新綜合濃度場.引入模糊集理論表示污染羽,對比綜合污染羽和單個污染羽的模糊集來更新各潛在污染源的權(quán)重,根據(jù)潛在污染源權(quán)重大小和綜合污染羽收斂形狀判斷真實污染源的個數(shù)和大概位置.算例結(jié)果表明:采用卡爾曼濾波方法可以成功反演識別出地下水污染中真實污染源的準(zhǔn)確個數(shù)和大概位置;引入模糊集理論表示污染羽,通過對比綜合污染羽和單個污染羽的模糊集,可以確定各潛在污染源的權(quán)重.

污染源識別；卡爾曼濾波；模糊集；蒙特卡羅；靈敏度分析

與地表水不同,地下水污染一般具有存在的隱蔽性和發(fā)現(xiàn)的滯后性等特點,使人們對地下水污染源的狀況缺乏了解,包括污染源的個數(shù)、空間位置以及釋放歷史等.展開對地下水污染源反演識別的研究工作,可以了解地下水污染源的特征,為設(shè)計污染修復(fù)方案提供決策依據(jù).

地下水污染源反演識別是,根據(jù)污染場地的水文地質(zhì)條件和動態(tài)監(jiān)測數(shù)據(jù)以及一些專家經(jīng)驗等輔助信息,對描述地下水污染的模擬模型進行反演求解,從而辨識確定含水層中地下水污染源的特征.目前已有的應(yīng)用于地下水污染源反演識別問題的方法有:地球物理探測法、污染因子分析法以及數(shù)理方程反演法[1-3]00,其中又以數(shù)理方程反演法應(yīng)用最為廣泛.

卡爾曼濾波法是基于概率統(tǒng)計的一種數(shù)理方程反演方法.早先Dokou等[4]以蒙特卡洛隨機地下水水流模型和水質(zhì)模型通過卡爾曼濾波方法反演識別了DNAPL源的最佳位置;之后又在污染源強度未知的三維模型中驗證了該方法的可行性[5].江思珉等[6]基于卡爾曼濾波方法和模糊集合理論成功反演識別了污染源的位置;后又結(jié)合單純形法對地下水污染源強度進行了識別[7].為適應(yīng)不同問題的需要,卡爾曼濾波方法逐漸發(fā)展出擴展、集合、無跡卡爾曼濾波等形式,Xu等[8]利用集合卡爾曼濾波方法在一個二維合成確定性含水層中有效的反演識別了污染源的位置、釋放時間和初始釋放濃度.為了解決反演過程中濾波不穩(wěn)定或發(fā)散問題,崔尚進等[9]通過對常規(guī)卡爾曼濾波方程中的協(xié)方差矩陣進行U-D分解,增加了該方法的數(shù)值穩(wěn)定性,并用一個二維模型成功反演識別了污染源位置.

已有的地下水污染源反演研究多是借助假想算例反演識別污染源的位置和強度,且算例中真實污染源的個數(shù)往往設(shè)定為1個, 本文針對同時存在有多個真實污染源的算例問題,采用卡爾曼濾波方法對真實污染源的個數(shù)及其大概位置進行識別.另外,在實際問題中,由于地下水系統(tǒng)的復(fù)雜性,地下水污染源反演問題往往存在諸多場地信息不確定性,如水文地質(zhì)參數(shù)的不確定性以及由調(diào)查數(shù)據(jù)的誤差或缺失引起的不確定性[9-12]00;為表示場地信息的不確定性,已有的研究通常事先確定模型中的某種參數(shù)作為隨機參數(shù),本文通過對所建立模擬模型中的參數(shù)進行靈敏度分析,篩選出靈敏度最高的參數(shù)作為隨機參數(shù), 其他參數(shù)作為確定性參數(shù),使隨機參數(shù)的選取有據(jù)可依.

1 研究方法

1.1 污染源反演識別流程

在現(xiàn)場調(diào)查、動態(tài)監(jiān)測和定性分析的基礎(chǔ)上,根據(jù)研究區(qū)的水文地質(zhì)條件和專家經(jīng)驗,確定污染源可能出現(xiàn)的范圍,初步估計污染源存在的位置、個數(shù)及每個潛在污染源的初始權(quán)重(權(quán)重表示該污染源為真實污染源的可能性大小,取值0~1之間).寫出地下水污染的水流和溶質(zhì)運移數(shù)學(xué)模型的一般表達式,對模型中的各項輸入變量均給出經(jīng)驗估計值(包括污染源釋放歷史),初步建立地下水水流數(shù)值模擬模型和溶質(zhì)運移數(shù)值模擬模型.考慮場地信息的不確定性,對輸入模擬模型中的參數(shù)進行靈敏度分析,篩選出靈敏度最高的參數(shù)作為隨機變量,其它參數(shù)作為確定性變量,根據(jù)現(xiàn)場調(diào)查和專業(yè)經(jīng)驗,對于隨機變量給出其取值區(qū)間,對于確定性變量給出其具體經(jīng)驗估計值;對于篩選出的隨機變量,在其取值區(qū)間內(nèi)采用拉丁超立方方法進行抽樣生成參數(shù)隨機場,運用蒙特卡羅方法將參數(shù)隨機場輸入模擬模型,生成溶質(zhì)濃度場,結(jié)合各潛在污染源的初始權(quán)重計算得到初始綜合濃度場和初始誤差協(xié)方差矩陣.根據(jù)采樣點實測數(shù)據(jù),運用卡爾曼濾波更新方程構(gòu)造合適的迭代過程,對綜合濃度場和協(xié)方差矩陣進行修正更新,對比用模糊集表示的綜合污染羽和各潛在污染源的單個污染羽,運用全局相似度公式,計算各潛在污染源的全局相似度來更新潛在污染源的權(quán)重.判斷權(quán)重是否穩(wěn)定,若穩(wěn)定則停止更新,否則選取新的采樣點繼續(xù)更新.污染源反演識別流程見圖1.

圖1 污染源反演流程

1.2 參數(shù)靈敏度分析

靈敏度分析方法有兩種,一是局部靈敏度分析,它的特點是只改變某一個參數(shù)的值而保持其他參數(shù)值不變,在此前提下分析某個參數(shù)發(fā)生變化時對模擬模型輸出結(jié)果的影響;二是全局靈敏度分析,它的特點是可以考慮不同參數(shù)共同作用對模擬模型輸出結(jié)果的影響,能夠計算出所有參數(shù)的總靈敏度[13].

本文運用局部靈敏度分析方法,篩選出對模擬模型影響較大的參數(shù):首先給待分析參數(shù)賦值,輸入模擬模型并得到輸出結(jié)果,然后將參數(shù)增加和減少一定的幅度并輸入模擬模型,可以獲得與輸入對應(yīng)的輸出結(jié)果.運用下式計算出該參數(shù)的靈敏度:

在實際問題中可以用下式來近似求得某一特定參數(shù)的靈敏度系數(shù):

1.3 拉丁超立方方法

1.4 計算初始綜合濃度場和初始誤差協(xié)方差矩陣

1.5 卡爾曼濾波更新

運用卡爾曼濾波狀態(tài)更新方程,結(jié)合采樣點處的濃度值,對綜合濃度場和協(xié)方差矩陣進行迭代更新,使綜合污染羽逐漸向真實污染羽收斂,達到反演識別真實污染源個數(shù)及其大概位置的目的.卡爾曼濾波的狀態(tài)更新方程如下:

在對初始綜合濃度場進行更新時,第一步迭代更新所需要的-和-既是1.4節(jié)計算得到的初始誤差協(xié)方差矩陣和初始綜合濃度場.

1.6 潛在污染源權(quán)重更新

為了判斷潛在污染源是否為真實污染源,每運用卡爾曼濾波方程更新一次綜合濃度場,都要計算新的綜合污染羽中各潛在污染源的權(quán)重大小.引入模糊集表示污染羽,即將綜合濃度場和各個潛在污染源的均值濃度場進行標(biāo)準(zhǔn)化(場內(nèi)各濃度值除以其中濃度最大值),設(shè)定4個模糊集標(biāo)準(zhǔn):

1=0.2,2=0.4,3=0.6,4=0.8

對比各模糊集標(biāo)準(zhǔn)下綜合污染羽和單個污染羽的相似度,并計算全局相似度來更新潛在污染源的權(quán)重,全局相似度計算公式如下:

式中:as表示某一給定的標(biāo)準(zhǔn)值;Ss表示在某一標(biāo)準(zhǔn)as下的綜合污染羽與單個污染羽比較所得兩者的重合面積(圖2灰色區(qū)域);g即是單個污染源的全局相似度.

計算所有潛在污染源的全局相似度后對它們進行標(biāo)準(zhǔn)化(各潛在污染源的全局相似度除以其中最大的全局相似度),標(biāo)準(zhǔn)化后的全局相似度即為新的各潛在污染源的權(quán)重.

2 假想模型算例分析

2.1 模型建立

假定一個1000m′1000m大小、含水層厚30m的研究區(qū),含水層非均質(zhì)各向同性,巖性主要為粗砂,水流為潛水非穩(wěn)定流,以底板為零基準(zhǔn)面,初始水位23m,建立水流模型時將研究區(qū)北部邊界和南部邊界概化為已知水頭邊界,水頭分別為27m(北)、22m(南),東、西兩側(cè)概化為隔水邊界;忽略蒸發(fā)騰散作用,研究區(qū)垂向上均勻接受降水入滲補給,年平均降水量為 550mm.研究區(qū)內(nèi)估計存在潛在污染源個數(shù)有3個(圖3),設(shè)定其中1號和2號為真實污染源.將污染質(zhì)看作為不發(fā)生化學(xué)轉(zhuǎn)化與生物遷移的保守污染物,污染泄漏量500m3/d,污染泄漏濃度40mg/L;模擬時長為800d;在建立溶質(zhì)運移模型時將南北兩側(cè)的已知水頭邊界概化為已知濃度邊界,東西兩側(cè)的隔水邊界概化為零通量邊界;根據(jù)專家經(jīng)驗,給3個潛在污染源賦予初始權(quán)重,分別為0.6,0.5, 0.7.

根據(jù)水文地質(zhì)條件,建立地下水水流數(shù)學(xué)模型:

在地下水水流數(shù)學(xué)模型的基礎(chǔ)上,建立地下水溶質(zhì)運移數(shù)學(xué)模型:

式中:為介質(zhì)孔隙度,無量綱;為污染質(zhì)濃度, mg/L;為水動力彌散系數(shù),m2/d;為滲流速度, m/d;表示單位時間單位液相體積內(nèi)污染質(zhì)質(zhì)量的增減量,1,3為已知濃度邊界;2,4為零通量邊界.

2.2 算例分析

將研究區(qū)剖分為10′10的單元網(wǎng)格.運用 GMS中的MODFLOW和MT3DMS工具箱完成研究區(qū)的水流和污染質(zhì)運移模擬.考慮到場地信息的不確定性,運用參數(shù)靈敏度分析方法篩選出靈敏度最高的參數(shù)作為隨機變量,其他參數(shù)作為確定性變量輸入模擬模型.篩選時先將各參數(shù)取均值(表1)輸入模擬模型,輸出對應(yīng)的模擬結(jié)果,然后將其中的某一個參數(shù)加上或減去其均值的10%,20%,其他參數(shù)不變輸入模擬模型,可以分別得到對應(yīng)參數(shù)取值下的模擬結(jié)果.

表1 參數(shù)的概率分布及取值情況

圖3 污染場地平面圖

利用(2)式計算該參數(shù)的靈敏度.同理,計算出所有參數(shù)的靈敏度后,對比篩選出靈敏度最高的參數(shù).由圖4可以看出,對模型結(jié)果影響最大的參數(shù)為滲透系數(shù).

圖4 參數(shù)靈敏度分析

圖5 滲透系數(shù)場的某次抽樣

圖6 真實滲透系數(shù)場

圖7 真實污染羽

利用拉丁超立方方法對篩選出的滲透系數(shù)進行抽樣,本次抽樣120次,即得到120個滲透系數(shù)隨機場,滲透系數(shù)場的某次抽樣結(jié)果見圖5.考慮到場地相關(guān)性,本文對抽出的滲透系數(shù)隨機場進行相關(guān)性排序[20],將排序后的滲透系數(shù)隨機場輸入模擬模型.圖6所示為真實滲透系數(shù)場,結(jié)合給定的真實污染源信息,輸入模型可獲得真實污染羽(圖7).算例中有3個潛在污染源,結(jié)合120個滲透系數(shù)隨機場可生成3′120個溶質(zhì)濃度場,在此基礎(chǔ)上按照1.4 節(jié)表述過程計算得到初始綜合濃度場和初始誤差協(xié)方差矩陣,然后運用卡爾曼濾波狀態(tài)更新方程更新綜合濃度場,圖8(a)~8(d)給出綜合污染羽的更新結(jié)果.

圖8(a)所示為未進行更新前的初始綜合污染羽,圖8(b),(c),(d)分別為利用1個、3個、5個采樣點數(shù)據(jù)更新綜合濃度場后的綜合污染羽,可以看出,在5次采樣更新完成之后綜合污染羽不斷收斂,最終相似于圖7所示的真實污染羽,且3號污染源已被完全排出污染羽范圍.圖8(e),(f)為更新前后潛在污染源權(quán)重分布圖,發(fā)現(xiàn)1號、2號、3號污染源權(quán)重分別由初始的0.6、0.5、0.7經(jīng)過5次采樣更新后變化為0.9、1、0.09.由5次采樣更新后的綜合污染羽形狀以及污染源的權(quán)重分布可以判斷:3個潛在污染源中1號和2號為真實存在的污染源,研究區(qū)內(nèi)真實污染源個數(shù)為2個,且污染羽內(nèi)濃度最高處即為污染源最可能存在的位置.

3 結(jié)論

3.2 對于含有多個真實污染源的反演識別問題,卡爾曼濾波方法仍然適用;本文經(jīng)過5次采樣更新后,準(zhǔn)確識別出研究區(qū)內(nèi)真實污染源的個數(shù)為2個,且真實污染源的位置最可能在污染羽內(nèi)濃度最高處.

3.3 本文采用一個含有多個真實污染源的假想算例,旨在展示卡爾曼濾波方法在該情境下的適用性,為減少工作量,研究區(qū)網(wǎng)格剖分數(shù)目較少,建議在實際應(yīng)用中適當(dāng)增加網(wǎng)格剖分數(shù)目,以提高采樣點數(shù)據(jù)對網(wǎng)格所表示區(qū)域的代表性.

[1] 王景瑞,胡立堂.地下水污染源識別的數(shù)學(xué)方法研究進展[J]. 水科學(xué)進展, 2017,28(6):943-952. Wang J R, Hu L T.Advances in mathematical methods of groundwater pollution source identification [J]. Advances in Water Science, 2017,28(6):943-952.

[2] 董海彪.基于克里格替代模型和改進的Bayesian-MCMC方法的地下水污染源反演識別研究[D]. 長春:吉林大學(xué), 2016. Identification of groundwater pollution sources based on kriging alternative model and improved Bayesian-MCMC methods [D]. Changchun: Jilin University, 2016.

[3] Milnes E, Perrochet P. Simultaneous identification of a single pollution point-source location and contamination time under known flow field conditions [J]. Advances in Water Resources, 2007,30(12):2439-2446.

[4] Dokou Z, Pinder G F. Optimal search strategy for the definition of a DNAPL source [J]. Journal of Hydrology, 2009,376(3):542-556.

[5] Dokou Z, Pinder G F. Extension and field application of an integrated DNAPL source identification algorithm that utilizes stochastic modeling and a Kalman filter [J]. Journal of Hydrology, 2011,398(3):277-291.

[6] 江思珉,張亞力,周念清,等.基于卡爾曼濾波和模糊集的地下水污染羽識別[J]. 同濟大學(xué)學(xué)報(自然科學(xué)版), 2014(3):435-440.Jiang S M, Zhang Y L, Zhou N Q, et al.Groundwater contaminant plume identification by using Kalman filter technique and fuzzy set theory [J].Journal of Tongji University (Natural Science), 2014(3):435-440.

[7] 江思珉,張亞力,周念清,等.基于污染羽形態(tài)對比的地下水污染源識別研究[J]. 水利學(xué)報, 2014,45(6):735-741.Jiang S M, Zhang Y L, Zhou N Q, et al.Groundwater contaminant source identification based on the morphological comparison of pollution plume [J].Journal of Hydraulic Engineering, 2014,45(6):735-741.

[8] Xu T, G MEZ-HERN NDEZ J J. Joint identification of contaminant source location, initial release time and initial solute concentration in an aquifer via ensemble Kalman filtering: Identification of a contaminant source [J]. Water Resources Research, 2016,52(8):6587-6595.

[9] 崔尚進,盧文喜,顧文龍,等.基于U-D分解的卡爾曼濾波法在地下水污染源識別中的應(yīng)用[J]. 中國環(huán)境科學(xué), 2016,36(9):2633-2637. Cui S J, Lu W X, Gu W L, et al.Application of U-D factorization- based Kalman filter to identify the groundwater pollution source [J]. China Environmental Science, 2016,36(9):2633-2637.

[10] 施小清,吳吉春,吳劍鋒,等.多個相關(guān)隨機參數(shù)的空間變異性對溶質(zhì)運移的影響[J]. 水科學(xué)進展, 2012,23(4):509-515. Shi X Q, Wu J C, Wu J F, et al.Effects of the heterogeneity of multiple correlated random parameters on solute transport [J]. Advances in Water Science, 2012,23(4):509-515.

[11] 李久輝,盧文喜,常振波,等.基于不確定性分析的地下水污染超標(biāo)風(fēng)險預(yù)警[J]. 中國環(huán)境科學(xué), 2017,37(6):2270-2277. Li J H, Lu W X, Chang Z B, et al.Risk prediction of groundwater pollution based on uncertainty analysis [J]. China Environmental Science, 2017,37(6):2270-2277.

[12] 束龍倉,朱元生.地下水資源評價中的不確定性因素分析[J]. 水文地質(zhì)工程地質(zhì), 2000,27(6):6-8.Shu L C, Zhu Y S.Analysis of uncertainties in groundwater resource evaluation [J].Hydrogeology and Engineering Geology, 2000,27(6):6-8.

[13] 束龍倉,王茂枚,劉瑞國,等.地下水?dāng)?shù)值模擬中的參數(shù)靈敏度分析[J]. 河海大學(xué)學(xué)報(自然科學(xué)版), 2007,35(5):491-495. Shu L C, Wang M M, Liu R G, et al.Sensitivity analysis of parameters in numerical simulation of groundwater [J]. Journal of Hohai University (Natural Sciences), 2007,35(5):491-495.

[14] Mckay M D, Beckman R J, Conover W J. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code [J]. Technometrics, 2000,42(1):55-61.

[15] Rajabi M M, Ataie-Ashtiani B, Janssen H. Efficiency enhancement of optimized Latin hypercube sampling strategies: Application to Monte Carlo uncertainty analysis and meta-modeling [J]. Advances in Water Resources, 2015,76:127-139.

[16] Helton J C, Davis F J. Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems [J]. Reliability Engineering & System Safety, 2003,81(1):23-69.

[17] 尹增謙,管景峰,張曉宏,等.蒙特卡羅方法及應(yīng)用[J]. 物理與工程, 2002,12(3):45-49. Yin Z Q, Guan J F, Zhang X H, et al. The Monte Carlo method and its application [J].Physics and Engineering, 2002,12(3):45-49.

[18] Densmore J D, Urbatsch T J, Evans T M, et al. A hybrid transport- diffusion method for Monte Carlo radiative-transfer simulations [J]. Journal of Computational Physics, 2007,222(2):485-503.

[19] Herrera G S, Pinder G F. Space-time optimization of groundwater quality sampling networks [J]. Water Resources Research, 2005, 41(12):410-411.

[20] Zhang Y Q, Pinder G. Latin hypercube lattice sample selection strategy for correlated random hydraulic conductivity fields [J]. Water Resources Research, 2003,39(8):472-477.

Application of Kalman filter to identify the groundwater contaminant sources.

BAI Yu-kun1,2, LU Wen-xi1,2*, LI Jiu-hui1,2

(1.Key Laboratory of Groundwater Resources and Environment Ministry of Education, Jilin University, Changchun 130012, China；2.College of Environment and Resources, Jilin University, Changchun 130012, China)., 2019,39(8)：3450~3456

Kalman filter was used to identify the number and approximate location of groundwater contaminant sources. Based on a hypothetical example, a groundwater flow and transport simulation model was established. The parameter that had the greatest impact on the simulation results was selected as a random variable by sensitivity analysis method. Then it was sampled and the sampling results were input into the simulation model by Monte Carlo method to generate the contaminant concentration field. Kalman filter method was used to update the composite concentration field one by one by using the measured concentration values at the sampling point. The fuzzy set theory was introduced to represent the pollution plume, and the weight of each potential contaminant source was updated by comparing the fuzzy sets of composite plume and individual plume. The number and approximate location of the real contaminant sources were judged according to the weight of potential contaminant sources and the convergence shape of composite plume. The results showed that the Kalman filter method can successfully identify the exact number and approximate location of the real contaminant sources in groundwater pollution; the fuzzy set theory was introduced to represent the pollution plume, and the weight of each potential contaminant source can be determined by comparing the fuzzy sets of the composite plume and the individual plume.

contaminant sources identification；Kalman filter；fuzzy set；Monte Carlo；sensitivity analysis

X523

A

1000-6923(2019)08-3450-07

白玉堃(1995-),男,吉林大學(xué)碩士研究生,主要從事地下水污染源反演識別問題方面的研究

2018-12-27

國家自然科學(xué)基金資助項目(41672232)

* 責(zé)任作者, 教授, luwenxi@jlu.edu.cn