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

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

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

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

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

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

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

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

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

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

1 研究方法

1.1 污染源反演識別流程

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

圖1 污染源反演流程

1.2 參數靈敏度分析

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

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

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

1.3 拉丁超立方方法

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

1.5 卡爾曼濾波更新

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

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

1.6 潛在污染源權重更新

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

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

對比各模糊集標準下綜合污染羽和單個污染羽的相似度,并計算全局相似度來更新潛在污染源的權重,全局相似度計算公式如下:

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

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

2 假想模型算例分析

2.1 模型建立

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

根據水文地質條件,建立地下水水流數學模型:

在地下水水流數學模型的基礎上,建立地下水溶質運移數學模型:

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

2.2 算例分析

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

表1 參數的概率分布及取值情況

圖3 污染場地平面圖

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

圖4 參數靈敏度分析

圖5 滲透系數場的某次抽樣

圖6 真實滲透系數場

圖7 真實污染羽

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

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

3 結論

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

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

[1] 王景瑞,胡立堂.地下水污染源識別的數學方法研究進展[J]. 水科學進展, 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]. 長春:吉林大學, 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]. 同濟大學學報(自然科學版), 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]. 水利學報, 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分解的卡爾曼濾波法在地下水污染源識別中的應用[J]. 中國環(huán)境科學, 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] 施小清,吳吉春,吳劍鋒,等.多個相關隨機參數的空間變異性對溶質運移的影響[J]. 水科學進展, 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] 李久輝,盧文喜,常振波,等.基于不確定性分析的地下水污染超標風險預警[J]. 中國環(huán)境科學, 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]. 水文地質工程地質, 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] 束龍倉,王茂枚,劉瑞國,等.地下水數值模擬中的參數靈敏度分析[J]. 河海大學學報(自然科學版), 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] 尹增謙,管景峰,張曉宏,等.蒙特卡羅方法及應用[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-),男,吉林大學碩士研究生,主要從事地下水污染源反演識別問題方面的研究

2018-12-27

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

* 責任作者, 教授, luwenxi@jlu.edu.cn