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

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

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

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

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

1 污染羽識別算法

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

1.1 基本假設(shè)

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

1.2 算法流程

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

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

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

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

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

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

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

2 污染羽的更新和新采樣點的選取

2.1 綜合污染羽的更新

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

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

濃度更新方程為

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

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

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

2.2 新采樣點的選取

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

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

3 污染羽的模糊對比

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

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

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

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

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

4 算例

4.1 問題概述

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

如圖3所示，研究區(qū)存在二維非均質(zhì)各向同性的承壓含水層，尺寸為400m×500m，用邊長為5m的正方形網(wǎng)格將研究區(qū)域剖分為100行、80列的有限差分網(wǎng)格.假設(shè)含水層水流運動為穩(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ī)變量，且服從對數(shù)正態(tài)分布，其均值和標(biāo)準(zhǔn)差σlnK分別為2.62和0.53，在x方向和y方向的相關(guān)長度均為140m，變差函數(shù)為指數(shù)型.圖4為滲透系數(shù)lnK場的一次實現(xiàn).根據(jù)已有的場地信息，污染源可能存在的位置如圖5所示，其中位置2為污染源的真實位置（事先未知）.本算例是在含水層參數(shù)不確定、污染源位置未知條件下通過擇優(yōu)選擇采樣點識別地下水污染羽.

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

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

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

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

4.2 算例分析

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

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

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

5 結(jié)論

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

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

圖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ù)的不確定性，并且污染源的強(qiáng)度假設(shè)是已知的.在下一步的研究中要考慮其他水文地質(zhì)參數(shù)的不確定性影響下污染源的位置和強(qiáng)度的同步識別.

［1］ 吳劍鋒，鄭春苗.地下水污染監(jiān)測網(wǎng)的設(shè)計研究進(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］ 王開麗，黃冠華.二維強(qiáng)非均質(zhì)含水層中滲透系數(shù)空間變異對污染物遷移的影響［J］.水利學(xué)報，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ù)的空間變異性對污染物運移的影響研究［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］ 施小清，吳吉春，吳劍鋒，等.多個相關(guān)隨機(jī)參數(shù)的空間變異性對溶質(zhì)運移的影響［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］ 吳劍鋒，彭偉，錢家忠，等.基于INPGA的地下水污染治理多目標(biāo)優(yōu)化管理模型：I——理論方法與算例驗證［J］.地質(zhì)論評，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)估計理論［M］.北京：科學(xué)出版社，2011.LIU Sheng，ZHANG Hongmei.Theory of optimal estimation［M］.Beijing：Science Press，2011.