基于卡爾曼濾波和模糊集的地下水污染羽識(shí)別

江思珉，張亞力，周念清，趙 姍

（同濟(jì)大學(xué) 土木工程學(xué)院，上海200092）

我國(guó)地下水污染狀況已十分嚴(yán)重.與地表水體相比，地下水一旦遭受污染其治理和修復(fù)通常需要更長(zhǎng)的時(shí)間（常持續(xù)30年以上），其治理和修復(fù)費(fèi)用往往極為昂貴.在地下水污染治理和修復(fù)過(guò)程中往往需要有污染羽的時(shí)空描述資料.對(duì)污染羽進(jìn)行時(shí)空描述通常需要通過(guò)建立相應(yīng)的地下水污染監(jiān)測(cè)網(wǎng)，以獲取地下水的物理特性、化學(xué)特性和生物特性在時(shí)間、空間上的變化特征數(shù)據(jù)［1］.

在地下水水質(zhì)的監(jiān)測(cè)過(guò)程中，過(guò)量取樣和取樣不足是比較常見(jiàn)的問(wèn)題.過(guò)度取樣會(huì)造成不必要的取樣和分析費(fèi)用，而取樣不足則不能獲得污染物隨時(shí)空變化的準(zhǔn)確信息，進(jìn)而不能確立地下水污染治理和修復(fù)的合理方案［2］.因此，文獻(xiàn)［3－7］通過(guò)對(duì)地下水污染監(jiān)測(cè)網(wǎng)進(jìn)行優(yōu)化設(shè)計(jì)，可以在不降低監(jiān)測(cè)精度和不影響評(píng)價(jià)風(fēng)險(xiǎn)的前提下減少不必要的取樣或增加必要的取樣，從而以盡可能少的費(fèi)用準(zhǔn)確地識(shí)別出地下水污染羽.

地下水系統(tǒng)是一個(gè)受多種因素影響的復(fù)雜系統(tǒng)，存在著各種不確定性.諸多不確定性的存在使得很難準(zhǔn)確地刻畫(huà)地下水中污染物的空間分布特征，主要的不確定性因素包括場(chǎng)地信息的不確定性和污染源（位置和強(qiáng)度）的不確定性等［8－10］.在污染源位置未知和場(chǎng)地信息不確定性的前提下，本文基于卡爾曼濾波技術(shù)和模糊集合理論提出一類地下水污染羽識(shí)別的新方法，能夠以盡可能少的采樣點(diǎn)達(dá)到污染羽準(zhǔn)確識(shí)別的目的.

1 污染羽識(shí)別算法

在進(jìn)行地下水污染治理和修復(fù)過(guò)程中，通常需要以盡可能少的水質(zhì)采樣點(diǎn)達(dá)到準(zhǔn)確識(shí)別地下水污染羽的目的.在場(chǎng)地信息不確定的前提下，可以選擇使污染濃度場(chǎng)的不確定度減少最大的點(diǎn)作為采樣點(diǎn)［11－12］.

1.1 基本假設(shè)

①污染源的強(qiáng)度已知；②地下水模型僅考慮滲透系數(shù)的不確定性.

1.2 算法流程

圖1為污染羽識(shí)別算法的流程圖.具體計(jì)算步驟如下：①確定污染源可能存在的位置，并對(duì)每個(gè)可能位置賦以概率值作為初始權(quán)重.權(quán)重值介于0～1之間，0表示該位置不存在污染源，1表示該位置存在污染源；②利用地質(zhì)統(tǒng)計(jì)軟件SGeMS生成場(chǎng)地的滲透系數(shù)隨機(jī)場(chǎng)［13］；③構(gòu)建地下水流模型與溶質(zhì)運(yùn)移模型，通過(guò)蒙特卡羅（Monte Carlo）隨機(jī)模擬得到污染濃度場(chǎng)庫(kù)；④根據(jù)污染源權(quán)重隨機(jī)選取污染濃度場(chǎng)，提取濾波區(qū)的濃度場(chǎng)，經(jīng)統(tǒng)計(jì)分析得到初始的綜合污染羽和濃度誤差協(xié)方差矩陣；⑤結(jié)合已有采樣點(diǎn)資料，利用卡爾曼濾波技術(shù)對(duì)綜合污染羽及誤差協(xié)方差矩陣進(jìn)行更新；⑥通過(guò)誤差協(xié)方差矩陣尋找新采樣點(diǎn)；⑦結(jié)合新采樣點(diǎn)利用卡爾曼濾波技術(shù)更新綜合污染羽；⑧利用模糊集合理論將綜合污染羽與單個(gè)的污染源所對(duì)應(yīng)的污染羽進(jìn)行對(duì)比，得到污染羽的相似度；⑨根據(jù)相似度更新污染源權(quán)重.○10判斷權(quán)重值和不確定度是否趨于穩(wěn)定，滿足則停止計(jì)算，否則回到步驟④.

圖1 污染羽識(shí)別算法流程Fig.1 Flow chart of the algorithm for identification of contaminant plume

1.3 綜合污染羽生成與誤差協(xié)方差估計(jì)

借鑒文獻(xiàn)［14］的個(gè)體適應(yīng)值庫(kù)的概念，在利用蒙特卡羅隨機(jī)模擬方法生成污染濃度場(chǎng)時(shí)將模擬結(jié)果存入污染濃度場(chǎng)庫(kù).如果污染源可能存在于m個(gè)位置，且滲透系數(shù)隨機(jī)場(chǎng)有q次實(shí)現(xiàn)，則將可能的污染源與滲透系數(shù)場(chǎng)進(jìn)行組合，可生成m×q個(gè)污染濃度場(chǎng).

根據(jù)污染源權(quán)重大小，從各個(gè)污染源所對(duì)應(yīng)的q個(gè)濃度場(chǎng)（對(duì)應(yīng)于滲透系數(shù)場(chǎng)的q次實(shí)現(xiàn)）分別隨機(jī)選取濃度場(chǎng).在所有被選取出的污染濃度場(chǎng)中將相同滲透系數(shù)場(chǎng)所對(duì)應(yīng)的濃度場(chǎng)疊加，得到當(dāng)前污染源權(quán)重下的濃度隨機(jī)場(chǎng)的q次實(shí)現(xiàn)，即統(tǒng)計(jì)意義上的q個(gè)樣本.

對(duì)濃度隨機(jī)場(chǎng)數(shù)據(jù)（q次實(shí)現(xiàn)）進(jìn)行統(tǒng)計(jì)分析，可以估計(jì)出綜合污染羽的濃度值以及濃度誤差協(xié)方差矩陣.濃度隨機(jī)場(chǎng)（q次實(shí)現(xiàn)）的均值即為綜合污染羽的濃度，節(jié)點(diǎn)j處的濃度值記為.濃度誤差協(xié)方差矩陣中的元素為R（ci，cj），其表達(dá)式為

式中：ci，k與cj，k分別表示在第k個(gè)污染濃度場(chǎng)中位置i和位置j處的濃度值分別表示在位置i和位置j處的平均濃度值（即綜合污染羽的濃度值）.需要說(shuō)明的是，上述對(duì)濃度隨機(jī)場(chǎng)（q次實(shí)現(xiàn)）數(shù)據(jù)的統(tǒng)計(jì)分析（綜合污染羽生成和誤差協(xié)方差矩陣估計(jì)）并未在整個(gè)研究區(qū)（含水層）進(jìn)行，而只在濾波區(qū)（即應(yīng)用卡爾曼濾波技術(shù)的區(qū)域）進(jìn)行.

2 污染羽的更新和新采樣點(diǎn)的選取

2.1 綜合污染羽的更新

卡爾曼濾波是由卡爾曼在1960年推導(dǎo)出的一種線性、無(wú)偏且誤差方差最小的優(yōu)化估計(jì)方法，它在已知系統(tǒng)模型和量測(cè)數(shù)學(xué)模型、量測(cè)噪聲統(tǒng)計(jì)特征及系統(tǒng)狀態(tài)的初始情況下利用量測(cè)數(shù)據(jù)獲取系統(tǒng)狀態(tài)變量的最優(yōu)估計(jì)值［15］.

利用卡爾曼濾波技術(shù)通過(guò)連續(xù)作用采樣點(diǎn)更新綜合污染羽與誤差協(xié)方差矩陣.卡爾曼濾波的更新方程如下.卡爾曼增益矩陣為

濃度更新方程為

誤差協(xié)方差更新方程為

式中：K為卡爾曼增益矩陣；P為n階濃度誤差協(xié)方差矩陣，下標(biāo)－，＋分別表示先驗(yàn)估計(jì)和后驗(yàn)估計(jì)，對(duì)應(yīng)利用測(cè)量數(shù)據(jù)更新之前和更新之后的狀態(tài)；H為量測(cè)矩陣，通常為l×n階的矩陣，l為取樣點(diǎn)的數(shù)目，n為用于計(jì)算的網(wǎng)格點(diǎn)個(gè)數(shù)（濾波區(qū)內(nèi)）；r為采樣誤差的方差；＾c是n維狀態(tài)向量，表示濃度估計(jì)值；z為采樣值；I為n階單位矩陣.通過(guò)連續(xù)添加單個(gè)采樣點(diǎn)對(duì)綜合污染羽及誤差協(xié)方差矩陣進(jìn)行更新，此時(shí)H為一個(gè)1×n階的矩陣，其形式如下：

其中，向量元素僅在采樣點(diǎn)處取1，其他位置取0.若已知采樣點(diǎn)的位置和采樣值，可以利用式（2）～（5）更新綜合污染羽.需要說(shuō)明的是，初始的＾c－和P－的計(jì)算方法見(jiàn)1.3節(jié).

2.2 新采樣點(diǎn)的選取

隨著采樣點(diǎn)的增加，污染濃度場(chǎng)的不確定度（即誤差協(xié)方差矩陣的對(duì)角線之和）會(huì)逐漸降低.新采樣點(diǎn)的選取方法如下：首先將濾波區(qū)已有采樣點(diǎn)之外的其他網(wǎng)格點(diǎn)視為備選采樣點(diǎn)，分別利用卡爾曼濾波的更新方程（式（4）和式（5））得到更新后的誤差協(xié)方差矩陣，然后計(jì)算各備選采樣點(diǎn)對(duì)應(yīng)的不確定度，最后選擇不確定度最小的點(diǎn)作為新采樣點(diǎn).

本文采用一種更為簡(jiǎn)便的方法，就是尋找更新后的不確定度下降最大的點(diǎn)，將其作為新采樣點(diǎn).更新后的不確定度可以表示為更新前的不確定度減去新增采樣點(diǎn)使不確定度下降的量［11］，即

3 污染羽的模糊對(duì)比

單個(gè)污染源的污染羽定義為該污染源在滲透系數(shù)隨機(jī)場(chǎng)的所有實(shí)現(xiàn)下的濃度均值.如果有5個(gè)污染源，則可以生成5個(gè)單個(gè)污染源的污染羽.通過(guò)將單個(gè)污染源的污染羽與利用卡爾曼濾波更新后的綜合污染羽進(jìn)行對(duì)比，對(duì)污染源可能所在位置的權(quán)重進(jìn)行更新，從而確定污染源位置.

首先將污染羽用模糊集來(lái)表示，隸屬函數(shù)為標(biāo)準(zhǔn)化的濃度值c′（定義為濃度值除以最大濃度值）.在考慮數(shù)個(gè)模糊集時(shí)，每個(gè)模糊集中的值均小于一個(gè)給定的標(biāo)準(zhǔn)αs.若αs分別取值為0.2，0.4，0.6，0.8，則共可得到4個(gè)模糊集.

在某一αs下，將單一污染源的污染羽的模糊集分別與綜合污染羽的模糊集進(jìn)行對(duì)比，計(jì)算交叉面積為Ss，如圖2所示的交叉面積表征單一污染源的污染羽與綜合污染羽的相似度.全局相似度g為各個(gè)標(biāo)準(zhǔn)αs下的相似度的加權(quán)和［12］，即

圖2為2個(gè)模糊集的比較，它一方面考慮了污染羽截面的形狀，另一方面也考慮對(duì)污染羽中濃度較高的區(qū)域賦以較大的權(quán)重.污染源新的權(quán)重通過(guò)污染源全局相似度進(jìn)行標(biāo)準(zhǔn)化后得到，具體就是將各相似度除以它們中的最大值，得到的就是標(biāo)準(zhǔn)化的相似度.

圖2 綜合污染羽與單個(gè)污染源污染的模糊比較Fig.2 Comparison of composite plume and individual source plume

4 算例

4.1 問(wèn)題概述

圖3 含水層結(jié)構(gòu)平面Fig.3 Plan view of confined aquifer

如圖3所示，研究區(qū)存在二維非均質(zhì)各向同性的承壓含水層，尺寸為400m×500m，用邊長(zhǎng)為5m的正方形網(wǎng)格將研究區(qū)域剖分為100行、80列的有限差分網(wǎng)格.假設(shè)含水層水流運(yùn)動(dòng)為穩(wěn)定流，含水層厚度為10m.上邊界為二類邊界，單寬流量為1.0 m2·d－1，下邊界為定水頭邊界，水頭為0 m，左、右邊界為隔水邊界.含水介質(zhì)孔隙度為0.3，縱向彌散度為40.0m，水平橫向彌散度為1.2m.已知污染源體積注入率和污染物濃度，分別為1m3·d－1和500 mg·L－1，但污染源位置未知.

含水層的滲透系數(shù)為一隨機(jī)變量，且服從對(duì)數(shù)正態(tài)分布，其均值和標(biāo)準(zhǔn)差σlnK分別為2.62和0.53，在x方向和y方向的相關(guān)長(zhǎng)度均為140m，變差函數(shù)為指數(shù)型.圖4為滲透系數(shù)lnK場(chǎng)的一次實(shí)現(xiàn).根據(jù)已有的場(chǎng)地信息，污染源可能存在的位置如圖5所示，其中位置2為污染源的真實(shí)位置（事先未知）.本算例是在含水層參數(shù)不確定、污染源位置未知條件下通過(guò)擇優(yōu)選擇采樣點(diǎn)識(shí)別地下水污染羽.

圖4 滲透系數(shù)隨機(jī)場(chǎng)的任意一次實(shí)現(xiàn)Fig.4 One randomly selected stochastic hydraulic conductivity field

圖5 真實(shí)的污染源位置和對(duì)應(yīng)的污染羽分布Fig.5 True source location and the corresponding plume distribution

本算例中，滲透系數(shù)隨機(jī)場(chǎng)共有100次實(shí)現(xiàn)，污染源可能存在的位置為5處，將滲透系數(shù)隨機(jī)場(chǎng)的每次實(shí)現(xiàn)分別與各個(gè)可能的污染源結(jié)合，進(jìn)行蒙特卡羅隨機(jī)模擬共得到500個(gè)污染濃度場(chǎng).為了減少調(diào)用模擬計(jì)算的次數(shù)，將這500個(gè)污染濃度場(chǎng)存入污染濃度場(chǎng)庫(kù).之后的計(jì)算中，只需要根據(jù)每次更新后的污染源權(quán)重隨機(jī)選取濃度場(chǎng)即可.卡爾曼濾波區(qū)為圖3虛線包圍的網(wǎng)格.

污染羽識(shí)別算法的主要參數(shù)設(shè)置如下：污染源初始權(quán)重分別設(shè)為0.80，0.85，0.90，0.95，1.00.污染羽的模糊對(duì)比采用3 個(gè)模糊集，對(duì)應(yīng)的標(biāo)準(zhǔn)為0.2，0.4和0.6.

4.2 算例分析

圖6a為根據(jù)污染源初始權(quán)重計(jì)算得到的綜合污染羽，圖6c—6h為利用1個(gè)、6個(gè)和7個(gè)采樣點(diǎn)更新后的污染羽及其對(duì)應(yīng)的污染源權(quán)重.

由圖6c—6h可知，經(jīng)過(guò)7 次采樣后，位置2 權(quán)重穩(wěn)定于1，其他位置權(quán)重接近于零，同時(shí)將圖6g與真實(shí)的污染羽形態(tài)進(jìn)行比較，兩者形態(tài)十分接近，表明只需要7個(gè)采樣數(shù)據(jù)就可以確定污染源的位置，并能夠達(dá)到較為準(zhǔn)確地識(shí)別污染羽的目的.此外，由圖6g可以看出，4個(gè)采樣點(diǎn)（點(diǎn)2，4，5，6）分布在臨近污染源的下游方向，其余3個(gè)采樣點(diǎn)（點(diǎn)1，3，6）起到控制污染羽形態(tài)的作用.上述結(jié)果表明，采樣點(diǎn)的分布主要考慮2個(gè)因素，一是距離污染源的位置，二是污染濃度場(chǎng)的不確定度.

圖7為濃度場(chǎng)的不確定度隨采樣點(diǎn)數(shù)目的變化情況.從圖中可以看出，隨著采樣點(diǎn)數(shù)目的增加，濃度場(chǎng)的不確定度隨之降低，并且當(dāng)采樣點(diǎn)的數(shù)目增加至7時(shí)，濃度場(chǎng)的不確定度降到了較低的水平，為237.45（mg·L－1）2，即平均標(biāo)準(zhǔn)差為0.24 mg·L－1，已達(dá)到污染羽識(shí)別的精度要求，結(jié)束計(jì)算.

5 結(jié)論

（1）基于卡爾曼濾波技術(shù)和模糊集合理論提出一類地下水污染羽識(shí)別的方法，算例研究表明該方法是一種有效的污染羽識(shí)別算法，能夠以盡可能少的采樣點(diǎn)確定污染源的位置，并較為準(zhǔn)確地達(dá)到識(shí)別污染羽的目的.

（2）該方法應(yīng)用了污染濃度場(chǎng)庫(kù)的概念，只需要在計(jì)算開(kāi)始時(shí)利用數(shù)值模型生成多個(gè)濃度場(chǎng)，此后的計(jì)算中不需要重新計(jì)算濃度場(chǎng)，節(jié)省了計(jì)算成本.

圖6 綜合污染羽與污染源權(quán)重隨采樣的變化Fig.6 The change of composite plume and its source location weights with sampling

圖7 不確定度隨采樣次數(shù)的變化Fig.7 Convergence of uncertainty with the procedure of sampling

（3）所建立的污染羽識(shí)別方法目前只是考慮了滲透系數(shù)的不確定性，并且污染源的強(qiáng)度假設(shè)是已知的.在下一步的研究中要考慮其他水文地質(zhì)參數(shù)的不確定性影響下污染源的位置和強(qiáng)度的同步識(shí)別.

［1］ 吳劍鋒，鄭春苗.地下水污染監(jiān)測(cè)網(wǎng)的設(shè)計(jì)研究進(jìn)展［J］.地球科學(xué)進(jìn)展，2004，19（3）：429.WU Jianfeng，ZHENG Chunmiao.Contaminant monitoring network design：recent advances and future directions ［J］.Advance in Earth Sciences，2004，19（3）：429.

［2］ Chien C C，M A Medina，G F Pinder，et al.Environmental modeling and management：theory，practice and future directions［M］.［S.l.］：Today Media Inc，2002.

［3］ Wu J，Zheng C，Chen C C.Cost－effective sampling network design for contaminant plume monitoring under general hydrogeological conditions ［J］.Journal of Contaminant Hydrology，2005，77（1／2）：41.

［4］ Dhar A，Datta B.Global optimal design of ground water monitoring network using embedded Kriging ［J］.Ground Water，2009，47（6）：806.

［5］ Owlia R R，Abrishamchi A，M Tajrishy.Spatial－temporal assessment and redesign of groundwater quality monitoring network：a case study ［J］.Environmental Monitoring and Assessment，2011，172（1／4）：263.

［6］ Datta B，Chakrabarty D，Dhar A.Optimal dynamic monitoring network design and identification of unknown groundwater pollution sources［J］.Water Resources Management，2009，23（10）：2031.

［7］ Kim K，Lee K. Optimization of groundwater－monitoring networks for identification of the distribution of a contaminant plume ［J］.Stochastic Environmental Research and Risk Assessment，2007，21（6）：785.

［8］ 王開(kāi)麗，黃冠華.二維強(qiáng)非均質(zhì)含水層中滲透系數(shù)空間變異對(duì)污染物遷移的影響［J］.水利學(xué)報(bào)，2010，41（4）：437.WANG Kaili， HUANG Guanhua. Effect of hydraulic conductivity on contaminant transport in highly heterogeneous two－dimensional aquifers［J］.Journal of Hydraulic Engineering，2010，41（4）：437.

［9］ 閻婷婷，吳劍鋒.滲透系數(shù)的空間變異性對(duì)污染物運(yùn)移的影響研究［J］.水科學(xué)進(jìn)展，2006，17（1）：29.YAN Tingting，WU Jianfeng.Impacts of the spatial variation of hydraulic conductivity on the transport fate of contaminant plume［J］.Advances in Water Science，2006，17（1）：29.

［10］ 施小清，吳吉春，吳劍鋒，等.多個(gè)相關(guān)隨機(jī)參數(shù)的空間變異性對(duì)溶質(zhì)運(yùn)移的影響［J］.水科學(xué)進(jìn)展，2012，23（4）：509.SHI Xiaoqing，WU Jichun，WU Jianfeng，et al.Effects of the heterogeneity of multiple correlated random parameters on solute transport［J］.Advances in Water Science，2012，23（4）：509.

［11］ Herrera G S，Pinder F G. Space－time optimization of groundwater quality sampling networks［J］.Water Resources Research，2005，41（12）：407.

［12］ Dokou Z，Pinder F G.Optimal search strategy for the definition of a DNAPL source［J］.Journal of Hydrology，2009，376：542.

［13］ Remy N，Boucher A，Wu B J.Applied geostatistics with SGeMS：a user’s guide［M］.New York：Cambridge University Press，2009.

［14］ 吳劍鋒，彭偉，錢(qián)家忠，等.基于INPGA的地下水污染治理多目標(biāo)優(yōu)化管理模型：I——理論方法與算例驗(yàn)證［J］.地質(zhì)論評(píng)，2011，57（2）：277.WU Jianfeng，PENG Wei，QIAN Jiazhong，et al.INPGAbased multi－objective management model for optimal design of groundwater remediation system：Ⅰ—methodology and its experimental validation ［J］.Geological Review，2011，2：277.

［15］ 劉勝，張紅梅.最優(yōu)估計(jì)理論［M］.北京：科學(xué)出版社，2011.LIU Sheng，ZHANG Hongmei.Theory of optimal estimation［M］.Beijing：Science Press，2011.