基于U-D 分解卡爾曼濾波地下水污染源溯源辨識

賈順卿,盧文喜*,李久輝,白玉堃 (1.吉林大學地下水與資源環境教育部重點實驗室,吉林 長春 130012；2.吉林大學新能源與環境學院,吉林 長春 130012)

地下水作為人類生存與社會發展必要的自然資源,近年來遭受的污染問題愈發嚴重.地下水污染難以監測,不易治理.地下水污染溯源識別是地下水污染修復與治理過程中的重要一環,能提供污染源空間分布、污染物釋放歷史等信息,為其治理提供參考依據.

卡爾曼濾波是求解地下水污染反問題的一種有效方法[1].卡爾曼濾波作為一種最優狀態估計方法,可以應用于受隨機干擾的動態系統[2-3]. Herrera[4],Schmidt 等[5]應用卡爾曼濾波對地下水水質監測網進行時空優化;Dokou[6]嘗試使用卡爾曼濾波創建一個最優的搜索策略,使用最少的水質樣本來識別DNAPL 來源[6],Jiang 等[7],江思珉等[8-9],顧文龍等[10]應用基于卡爾曼濾波技術和模糊集理論識別污染羽的方法.常規卡爾曼濾波可能出現矩陣協方差非對稱或者負定的情況[11-12],崔尚進將卡爾曼濾波協方差矩陣進行U-D 分解,提高了卡爾曼濾波的穩定性[13].Chen 等[14], Xu 等[15-16]使用集合卡爾曼濾波對二維確定性含水層中的污染源進行了溯源辨識[14-16];場地信息存在不確定性[17-18],白玉堃等[19]采用靈敏度分析篩選出靈敏度最高的參數作為隨機變量進行污染源反演.

結合前人研究,本文運用U-D 分解卡爾曼濾波識別污染源的個數與位置,在此基礎上建立優化模型,采用克里格插值法建立替代模型并嵌入優化模型,運用遺傳算法求解優化模型得到污染源源強.

1 研究方法

1.1 技術路線

圖1 污染源反演技術路線Fig.1 Flow chart of contamination source identification

首先應用U-D 分解卡爾曼濾波反演污染源的個數與位置.在現場調查與動態監測的基礎上,結合專家意見初步估計污染源潛在的個數、位置以及初始權重,構建地下水溶質運移模擬模型.考慮場地信息存在不確定性,對場地參數進行靈敏度分析,篩選出靈敏度最高的參數作為隨機變量,其余參數作為確定型變量.給出隨機變量一個取值范圍,由拉丁超立方抽樣生成參數隨機場,采用蒙特卡羅法將參數隨機場輸入模擬模型生成溶質濃度場庫,結合各潛在污染源初始權重計算得到初始綜合濃度場與誤差協方差矩陣.對誤差協方差矩陣進行U-D 分解,結合采樣點數據運用U-D 分解卡爾曼濾波對污染濃度場和潛在污染源權重進行更新,權重穩定后判斷污染源的個數與位置.

在識別出污染源個數與位置的基礎上建立優化模型,進行污染源釋放強度反演.以污染物監測濃度與模擬計算濃度擬合誤差極小化為目標函數,污染源釋放強度為待求變量,替代模型作為優化模型等式約束條件,同時考慮源強上下限等約束條件,建立非線性規劃優化模型,運用遺傳算法求解優化模型反演出污染源釋放強度.為了減小求解優化模型時調用模擬模型產生的運算負荷,建立溶質運移模擬模型的克里格替代模型.

1.2 靈敏度分析與濃度場庫建立

靈敏度分析可反應模型輸出對參數變化的敏感程度,本文采用Morris 全局靈敏度分析方法評價參數變化對模型輸出的影響,選出影響程度最大的參數進行研究.Morris 設計采用“一次只改變一個參數”的抽樣取值方法,輪流計算各參數的目標函數值,從而得到各個參數全局靈敏度及參數間相關性和非線性的定性描述[20-21]:

把研究區用網格剖分,將篩選出的參數進行拉丁超立方抽樣得到n 組隨機變量參數場,只考慮一個潛在污染源時,輸入模擬模型得到對應的n 組溶質濃度場,m 個潛在污染源共m ×n 組.計算每個網格上參數場取不同值時濃度值的均值,就得到了該污染源的均值濃度場.在同一隨機變量參數場下,m個污染源的溶質濃度場按照初始權重進行加權疊加得到疊加濃度場.將m 個污染源的均值濃度場按初始權重進行加權疊加得到初始綜合濃度場.

1.3 U-D 分解卡爾曼濾波

常規卡爾曼濾波存在數值穩定性較差,計算復雜等問題,針對這一點可引入U-D 分解改進.

常規卡爾曼濾波更新方程[8-9]:

式中:K 為卡爾曼增益矩陣,P 為n 階誤差協方差矩陣,C 為n 維狀態向量,表示濃度估計值;-、+號表示先驗估計與后驗估計;Z 為采樣值,r 是采樣誤差的方差,H 是1×n 階量測矩陣,采樣點處矩陣元素為1,其余位置為0,H = [0…0,1,0 … 0];I 是n 階單位向量.

P 矩陣更新時矩陣相減易導致矩陣穩定性降低,甚至引發濾波發散,將P 矩陣分解為P=UDUT的形式,方程表示如下:

其中

式中:U 為單位上三角矩陣;D 為正定對角矩陣.

1.4 權重更新

每個潛在污染源均值濃度場對應的污染羽為單個污染羽,綜合濃度場對應的污染羽為綜合污染羽.確定污染源的個數與位置需要通過污染源權重判斷,將綜合污染羽與單個污染羽標準化,用模糊集表示(所有濃度值除以最大濃度值),每個模糊集合的元素隸屬度都大于等于給定的αi,本文取值為0.2、0.4、0.6、0.8.

將以模糊集形式表示的綜合污染羽與各單個污染羽進行比較,記錄二者的公共面積Si,計算出全局相似度并用各污染源全局相似度除以各污染源全局相似度中的最大者,即可算得各污染源的權重.全局相似度為[7-9]:

式中:g 為全局相似度;αi為模糊集標準值,i=1,2,…;Si為綜合污染羽與單個污染羽在同一αi下的交叉面積.

圖2 污染羽對比示意Fig.2 Comparison of pollution plumes

1.5 優化模型建立

以污染物監測濃度與模擬計算濃度擬合誤差極小化為目標函數,污染源釋放強度為待求變量,替代模型作為優化模型等式約束條件,同時考慮源強上下限等約束條件,建立非線性規劃優化模型

目標函數:

約束條件:

式中, q 為污染源釋放強度,n 為采樣點個數,ck為替代模型輸出的模擬值,c?k為實測值,f 為替代模型.

1.6 替代模型

應用遺傳算法求解優化模型的過程中會大量調用溶質運移模擬模型,使用克里格方法建立溶質運移模擬模型的替代模型能有效減少運算負荷.替代模型回歸方程為

2 假想算例

2.1 問題概述

圖3 污染場地Fig.3 The plan of contaminated site

假定研究區大小為1000m×1000m,潛水含水層厚度為20m,介質為粗砂,含水層非均質各向同性.忽略非飽和帶的影響.南北邊界為已知水頭邊界,以含水層底板為基準面,北邊界水位為16m,南邊界水位為14m;東西邊界為隔水邊界.模型模擬時間長度為500d,分為5 個時段,每個時段為100d.研究區接受降雨入滲補給,5 個時段降水量分別為100,300, 100,50,100mm.忽略研究區蒸發蒸騰作用的影響.結合場地信息與專家經驗,估計潛在污染源個數為4 個(圖3),并給出潛在污染源的初始權重(權重表示潛在污染源為真實污染源的可能性大小,取值0～1 之間),Ⅰ、Ⅱ、Ⅲ、Ⅳ號污染源分別為0.6、0.7、0.7、0.6.設定Ⅲ號為真實污染源,污染源泄漏流量為2000mg/d.

2.2 問題分析

2.2.1 位置個數反演 將研究區剖分成20×20 的網格,應用GMS 軟件的MODFLOW 和MT3DMS 模塊建立溶質運移模擬模型,考慮場地信息的不確定性,對模型中的參數進行Morris 全局靈敏度分析,篩選出滲透系數K、縱向彌散系數Ld、給水度Sy、降雨入滲補給系數α靈敏度最高的參數作為隨機變量.共有15 種組合方式,具體方法參照文獻[21].

表1 參數組合方式Table 1 Combination-patterns of parameters

圖4 參數靈敏度分析Fig.4 Sensitivity analysis of parameters

從圖4 可以看出,滲透系數靈敏度最高,本研究將滲透系數作為隨機變量,其余參數作為確定性變量,應用拉丁超立方抽樣生成110 組滲透系數隨機場并進行相關性排序,場地存在4 個潛在污染源,輸入模型得到4×110 組溶質濃度場,計算得出初始綜合濃度場和誤差協方差矩陣,Ⅲ號污染源的均值濃度場作為真實濃度場.

圖5 真實污染羽Fig.5 True pollution plume

將采樣點數據帶入U-D 分解卡爾曼濾波更新方程對綜合濃度場和污染源權重進行更新,圖為濃度場和權重更新結果.

圖5 為污染源真實污染羽,圖6 為初始綜合污染羽.從圖7 可以看出,經過4 次采樣更新后,綜合污染羽逐漸收斂于Ⅲ號污染源,與真實污染羽相似.4 個污染源權重從0.6,0.7,0.7,0.6 變為0,0.2,1,0.05,可以判斷出Ⅲ號污染源為真實污染源,更新完成后污染羽中濃度最高處就是污染源位置.

圖6 初始綜合污染羽Fig.6 Initial composite pollution plume

圖7 污染羽與權重更新Fig.7 Update composite pollution plume and source location weight

2.2.2 源強反演 在運用U-D 分解卡爾曼濾波反演出污染源的個數和位置的基礎上,建立優化模型并應用遺傳算法求解.采用克里格插值法建立模擬模型的替代模型,應用拉丁超立方抽樣獲得80 組源強,輸入模擬模型得到4 口采樣井的輸出,生成80 組輸入-輸出樣本,前60 組作為訓練樣本,后20 組作為檢驗樣本,將平均相對誤差(MRE)作為克里格替代模型近似精度的評價指標,結果見表2.

表2 替代模型精度評價Table 2 Accuracy assessment of surrogate model

由表2 可知,誤差精度不超過1%,替代模型符合精度要求.

以污染物監測濃度與模擬計算濃度擬合誤差極小化為目標函數,污染源釋放強度為待求變量,替代模型作為優化模型等式約束條件,同時考慮源強上下限等約束條件,建立非線性規劃優化模型.采用MATLAB 軟件中的遺傳算法工具箱求解優化模型,經過 46 次迭代求得源強 1972mg/d,與真實值2000mg/d 接近,相對誤差1.4%,源強反演精度較高.

3 結論

3.1 將協方差矩陣進行U-D 分解,能夠避免在迭代計算過程中卡爾曼濾波因為協方差矩陣可能出現不對稱或者負定,導致迭代過程發散的情況,保證了卡爾曼濾波的數值穩定性,使迭代過程更容易收斂.算例表明基于U-D 分解的卡爾曼濾波,能夠準確地辨識含水層中污染源的個數與位置.

3.2 采用局部靈敏度分析方法篩選出滲透系數K對模型輸出影響最大,作為反映場地信息不確定性的隨機變量輸入模型建立濃度場庫.

3.3 針對污染源釋放強度辨識問題,建立非線性規劃優化模型.應用遺傳算法求解非線性規劃優化模型,能快速計算出污染源釋放強度.應用克里格方法建立地下水污染質數值模擬模型的替代模型,將克里格替代模型代替模擬模型嵌入優化模型,能夠在保證精度的同時減少大量的運算負荷和運算時間.