張瑞城，周念清，江思珉，夏學(xué)敏

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

隨著現(xiàn)代工農(nóng)業(yè)生產(chǎn)的快速發(fā)展和污染物排放量的不斷增加，各種污染物通過不同的途徑入滲進(jìn)入到地下水體中，導(dǎo)致地下水污染問題日趨嚴(yán)重。在污染場(chǎng)地研究中，除了采用傳統(tǒng)的監(jiān)測(cè)、取樣、分析和評(píng)價(jià)外，通過建立數(shù)值模型對(duì)地下水污染物的遷移轉(zhuǎn)化過程進(jìn)行模擬預(yù)測(cè)已成為研究地下水污染問題的有效途徑和重要手段。準(zhǔn)確獲取模型參數(shù)是進(jìn)行數(shù)值模擬預(yù)測(cè)的關(guān)鍵，然而，由于參數(shù)存在空間變異性的特點(diǎn)，難以利用稀疏的觀測(cè)信息來準(zhǔn)確描述非均質(zhì)場(chǎng)地的模型參數(shù)［1］。近些年來，將地球物理方法（如高密度電阻率成像法、地質(zhì)雷達(dá)法等）和間接觀測(cè)手段引入污染水文地質(zhì)研究中，其應(yīng)用越來越廣泛［2-4］。

求解地下水逆問題通常需要進(jìn)行數(shù)據(jù)同化，利用濃度、水頭值以及采用地球物理方法獲取的間接觀測(cè)值來對(duì)模型參數(shù)進(jìn)行反演，從而對(duì)系統(tǒng)狀態(tài)提供更科學(xué)可靠的估計(jì)［5］。其中，基于集合的數(shù)據(jù)同化算法有集合卡爾曼濾波器（EnKF）和集合平滑器（ES）等，在地下水逆問題中已被廣泛使用［6-8］，二者的區(qū)別是EnKF需要同時(shí)更新模型參數(shù)和狀態(tài)，而ES僅進(jìn)行參數(shù)更新，就能有效避免EnKF中參數(shù)和狀態(tài)的不一致問題［9］。為了提高計(jì)算的運(yùn)行效率和準(zhǔn)確性，有學(xué)者提出了一種迭代局部更新集合平滑（iterative local updating ensemble smoother，ILUES）算法，通過采用局部更新策略，實(shí)現(xiàn)了對(duì)高維強(qiáng)非線性問題的求解［9-10］。

國內(nèi)外有關(guān)地下水污染參數(shù)反演研究，大多集中于污染源參數(shù)（源強(qiáng)和位置）或場(chǎng)地水文地質(zhì)參數(shù)（滲透系數(shù)場(chǎng)、孔隙度等）等單一類型參數(shù)的識(shí)別，而對(duì)污染源源強(qiáng)和滲透系數(shù)場(chǎng)的聯(lián)合反演少有研究［11-12］。與此同時(shí)，現(xiàn)有數(shù)據(jù)同化方法主要通過融合水頭與污染物濃度觀測(cè)值進(jìn)行參數(shù)反演，融合地球物理數(shù)據(jù)和水頭濃度數(shù)據(jù)的研究相對(duì)較少［1，13］。有鑒于此，本文綜合考慮污染源參數(shù)與場(chǎng)地水文地質(zhì)參數(shù)，將污染源源強(qiáng)和滲透系數(shù)場(chǎng)作為未知模型參數(shù)，構(gòu)建基于迭代局部更新平滑器算法（ILUES）的數(shù)據(jù)同化框架，并融合ERT（高密度電阻率法）數(shù)據(jù)進(jìn)行聯(lián)合參數(shù)反演，通過數(shù)值算例驗(yàn)證了該方法的有效性。

1 研究方法

1.1 地下水流與溶質(zhì)運(yùn)移模型

對(duì)于飽和多孔介質(zhì)中的三維地下水流運(yùn)動(dòng)，有如下基本微分方程：

式中：Kx，Ky，Kz分別為滲透系數(shù)在X，Y，Z方向上的分量，m·d-1；H為水頭，m；W為單位時(shí)間從單位體積含水層中流入或流出的水量，d-1；Ss為的貯水率，m-1；t為時(shí)間，d。

針對(duì)飽和多孔介質(zhì)中的溶質(zhì)運(yùn)移問題，其基本控制方程如下：

式中：C為濃度，g·m-3；t為時(shí)間，d；Dij為水動(dòng)力彌散系數(shù)，m2·d-1；vi為含水介質(zhì)中的實(shí)際水流速度，m·d-1；qs為單位體積含水層的源匯項(xiàng)流量，d-1；Cs為源匯項(xiàng)中溶質(zhì)濃度，g·m-3；∑Rn為化學(xué)反應(yīng)項(xiàng)總和，g·m3·d。

在建立地下水流三維運(yùn)動(dòng)方程和溶質(zhì)運(yùn)移模型的基礎(chǔ)上，分別采用MODFLOW程序［14］求解地下水流模型，用MT3DMS程序［15］求解溶質(zhì)運(yùn)移方程。

1.2 地層電阻率與溶質(zhì)濃度關(guān)系的構(gòu)建

溶質(zhì)具有導(dǎo)電性，不同的溶質(zhì)其導(dǎo)電性與電阻率差異比較大。Archie公式用于刻畫溶液電阻率與地層電阻率以及含水飽和度之間關(guān)系［16］，其有效性得到了許多學(xué)者的證實(shí)和采用，基本形式如下：

式中：φ為有效孔隙度；Sw為含水飽和度；Rw為溶液電阻率，Ω·m；Rt為地層總電阻率，Ω·m；m為黏結(jié)指數(shù)；n為飽和度指數(shù)。

溶液電阻率與溫度和溶液濃度等因素的有關(guān)，Sen［12］綜合考慮了這兩種因素對(duì)溶液電阻率的影響，將物理化學(xué)理論與室內(nèi)實(shí)驗(yàn)相結(jié)合，提出了如下模型公式。

式中：T為溶液的溫度，根據(jù)場(chǎng)地調(diào)查數(shù)據(jù)取為15℃，其余參數(shù)同前。

1.3 地球物理反演模型

高密度電阻率法以巖土體導(dǎo)電性差異為基礎(chǔ)，通過供電電極（A、B）向地下空間輸入穩(wěn)定的電流I，然后根據(jù)測(cè)量電極（M、N）的電位差分析電場(chǎng)的分布，進(jìn)而推算出視電阻率的空間分布。視電阻率正演模型可表示如下：

地球物理反演模型是利用反演方法通過對(duì)ERT數(shù)據(jù)（觀測(cè)值）的擬合實(shí)現(xiàn)正演模型參數(shù)（地層電阻率）的反演識(shí)別，可以利用開源程序Pygimli構(gòu)建地球物理模型并求解電阻率的反演問題［17］。

1.4 基于ILUES算法的數(shù)據(jù)同化框架

ILUES算法是一種迭代形式的ES算法，其協(xié)方差矩陣擾動(dòng)策略來源于Emerick和Reynolds［18］所提出的ES-MDA算法，通過添加擾動(dòng)的觀測(cè)誤差協(xié)方差矩陣對(duì)觀測(cè)數(shù)據(jù)進(jìn)行多次同化，從而實(shí)現(xiàn)對(duì)參數(shù)樣本的迭代更新。

為了更好地解決高維非線性參數(shù)反演問題，ILUES算法沒有采用傳統(tǒng)集合平滑方法的全局更新策略，而是使用了一種優(yōu)化的局部更新策略，即只對(duì)每個(gè)參數(shù)樣本的局部樣本集合進(jìn)行更新。在ILUES算法中，對(duì)于集合中的每一個(gè)樣本，其局部樣本集合通過以下距離評(píng)價(jià)因子進(jìn)行篩選：

本文構(gòu)建了融合ILUES算法、Archie公式、Sen模型和溶質(zhì)運(yùn)移模型的數(shù)據(jù)同化框架，其濃度和水頭及地球物理數(shù)據(jù)（ERT數(shù)據(jù)）作為數(shù)據(jù)同化的觀測(cè)值。算法流程如圖1所示。

圖1 基于ILUES算法的數(shù)據(jù)同化框架Fig.1 The framework of data assimilation based on ILUES algorithm

2 算例研究

研究區(qū)為矩形區(qū)域（二維XZ剖面，100m×40m），概化為非均質(zhì)各向異性的承壓含水層，用邊長為2m的正方形網(wǎng)格將研究區(qū)剖分為20行50列有限差分網(wǎng)格（如圖2）。假設(shè)含水層水流運(yùn)動(dòng)為穩(wěn)定流，研究區(qū)下邊界為隔水邊界，左右邊界為定水頭邊界（左邊界水頭為50.0m，右邊界水頭為49.0m），上邊界無源匯項(xiàng)。含水介質(zhì)孔隙度為0.30，縱向彌散度為10m，垂直橫向彌散度為3m。初始時(shí)刻，研究區(qū)無污染物。含水層滲透系數(shù)K滿足對(duì)數(shù)正態(tài)分布，其均值（ˉˉˉˉˉlnK）與方差（σ2lnK）分別為2.0和1.0，x，z方向的相關(guān)長度分別為60m和24m，變差函數(shù)為指數(shù)型，真實(shí)滲透系數(shù)場(chǎng)見圖2。

圖2 污染場(chǎng)地概念模型Fig.2 The conceptual model of contaminated site

假定研究區(qū)存在3個(gè)潛在點(diǎn)源污染，分別標(biāo)記為S1、S2和S3（S2處實(shí)際上無污染物釋放），污染物進(jìn)入含水層的釋放強(qiáng)度如表1，場(chǎng)地內(nèi)布置有12個(gè)監(jiān)測(cè)點(diǎn)（圖2）。觀測(cè)點(diǎn)的濃度觀測(cè)值是利用地下水流和溶質(zhì)運(yùn)移模型正演計(jì)算得到。模擬總時(shí)長為200d，均勻離散為10個(gè)應(yīng)力期，且只在前6個(gè)應(yīng)力期內(nèi)發(fā)生污染物釋放，其污染源強(qiáng)度見表1。由正演模擬值和正態(tài)分布的測(cè)量噪聲計(jì)算得到觀測(cè)值。

表1 各應(yīng)力期的污染源源強(qiáng)真實(shí)值Tab.1 The reference true values of pollution sources strength in different stress periods

式中：Csimu表示模擬值；Cobs為觀測(cè)值；L為噪聲水平（取10%）；δ為滿足標(biāo)準(zhǔn)高斯分布的隨機(jī)偏差。

高密度電阻率成像法采用溫納裝置進(jìn)行探測(cè)，沿地表布置一排電極，參見圖2，共獲取了10組ERT觀測(cè)數(shù)據(jù)，ERT測(cè)量值相對(duì)誤差取10%。為了研究利用場(chǎng)地觀測(cè)數(shù)據(jù)進(jìn)行污染溯源和滲透系數(shù)場(chǎng)估計(jì)的可行性，設(shè)計(jì)了算例1和算例2，分別通過同化傳統(tǒng)觀測(cè)數(shù)據(jù)（水頭與濃度）、地球物理數(shù)據(jù)（ERT數(shù)據(jù)）來反演未知模型參數(shù)。算例1中選取12個(gè)觀測(cè)點(diǎn)（圖2中三角形）的水頭和t＝［20，30，40，50，60，70，80］day的濃度作為觀測(cè)值；算例2選取t＝［20，40，60，80］d的ERT監(jiān)測(cè)數(shù)據(jù)作為觀測(cè)值，算例的具體設(shè)置見表2。

表2 算例設(shè)置Tab.2 The setting of case studies

考慮對(duì)空間上連續(xù)變化的滲透系數(shù)場(chǎng)進(jìn)行反演時(shí)，滲透系數(shù)場(chǎng)的維度與模型的網(wǎng)格剖分程度有關(guān)（本算例為20×50＝1 000維）。高維地下水逆問題的求解效率較低，一般需要對(duì)滲透系數(shù)場(chǎng)進(jìn)行降維，采用Karhunen-Loève展開［19］的方法對(duì)其進(jìn)行降維。

式中：ξi是獨(dú)立標(biāo)準(zhǔn)高斯分布的隨機(jī)數(shù)；τi和si(x，y)是特征值和特征函數(shù)；NKL是KL展開的項(xiàng)數(shù)。

為了定量評(píng)估ILUES算法同化2類不同觀測(cè)數(shù)據(jù)的參數(shù)反演效果，引入均方根誤差（root-meansquare error，RMSE）作為評(píng)價(jià)標(biāo)準(zhǔn)，該值反映了由ILUES算法反演得到的模型參數(shù)估計(jì)值與其真實(shí)值之間的差距，RMSE的值越?。ㄚ吔?），說明參數(shù)反演的精確度越高。RMSE的計(jì)算如下：

式中：y?n為參數(shù)反演得到的模型參數(shù)估計(jì)值；yn為模型參數(shù)的真實(shí)值；Ne為ILUES算法中集合的大小。

3 結(jié)果與討論

在2個(gè)算例中，滲透系數(shù)場(chǎng)的KL展開項(xiàng)數(shù)NKL＝60，可以保留lnK場(chǎng)約95%的變異性，從而對(duì)滲透系數(shù)場(chǎng)（20×50＝1 000維）的識(shí)別降維為對(duì)KL展開的60個(gè)高斯隨機(jī)數(shù)的反演。對(duì)于潛在污染點(diǎn)源S1、S2和S3，污染物釋放發(fā)生在前6個(gè)應(yīng)力期，共有18個(gè)污染源強(qiáng)度參數(shù)需要識(shí)別，其先驗(yàn)分布為［0，20］的均勻分布。因此，對(duì)于2個(gè)算例總共需要反演識(shí)別的未知模型參數(shù)數(shù)量均為78個(gè)。為了獲得準(zhǔn)確的參數(shù)反演結(jié)果，根據(jù)以往數(shù)值算例的經(jīng)驗(yàn)，將ILUES算法的參數(shù)設(shè)置為：集合數(shù)Ne＝1 000，迭代次數(shù)NIter＝10，α＝0.1。

3.1 污染源源強(qiáng)反演結(jié)果

圖3為算例1和算例2的污染源源強(qiáng)的反演結(jié)果。圖3a～3c是利用水頭和［20，30，40，50，60，70，80］d的濃度觀測(cè)值，得到的源強(qiáng)反演結(jié)果，經(jīng)10次迭代后基本收斂，但與各應(yīng)力期的源強(qiáng)真值仍有較大的偏差；相較而言，圖3d～3f融合［20，40，60，80］d的ERT時(shí)變數(shù)據(jù)反演得到的污染源源強(qiáng)，經(jīng)7次迭代后已經(jīng)收斂于真實(shí)值，其中S2各應(yīng)力期的源強(qiáng)值均接近于0（即無污染發(fā)生）。

圖3 污染源源強(qiáng)的反演結(jié)果Fig.3 The inversion results of contaminant source strength

因此，融合ERT數(shù)據(jù)的ILUES算法可以更好地估計(jì)污染源源強(qiáng)。值得注意的是，算例2中ERT數(shù)據(jù)測(cè)量時(shí)間為［20，40，60，80］d，而算例1中污染物濃度觀測(cè)時(shí)間為［20，30，40，50，60，70，80］d，在較少時(shí)間點(diǎn)采集的ERT數(shù)據(jù)的反演結(jié)果反而更好，其原因在于ERT等地球物理方法雖然測(cè)量精度不如水頭、濃度數(shù)據(jù)，但可以便捷地獲取大量具有空間連續(xù)性的觀測(cè)數(shù)據(jù)，其采集的觀測(cè)數(shù)據(jù)量（620個(gè)）遠(yuǎn)多于算例1中的水頭與濃度數(shù)據(jù)（96個(gè)）。

為了更好地量化和對(duì)比使用ERT數(shù)據(jù)和水頭濃度數(shù)據(jù)的污染源源強(qiáng)反演效果，將算例1和算例2中源強(qiáng)反演結(jié)果的RMSE值列于表3。由表3可知，使用ERT數(shù)據(jù)作為觀測(cè)值的算例2，任一污染源的任一應(yīng)力期對(duì)應(yīng)的均方根誤差均小于使用水頭與濃度數(shù)據(jù)的算例1，這進(jìn)一步說明了利用ILUES算法融合ERT數(shù)據(jù)進(jìn)行參數(shù)反演可以獲得更加精確的結(jié)果。

表3 污染源參數(shù)反演結(jié)果的均方根誤差Tab.3 Root-mean-square error of Case 1 and Case 2

3.2 滲透系數(shù)場(chǎng)的反演結(jié)果

圖4和圖5是經(jīng)過10次迭代后lnK的反演結(jié)果。對(duì)于算例1和算例2，作為參照?qǐng)龅膌nK真實(shí)場(chǎng)如圖4a和圖5a；初始時(shí)刻隨機(jī)生成的lnK場(chǎng)如圖4b和圖5b，為了有效對(duì)比算例1與算例2的反演精度，兩算例采用相同的初始隨機(jī)lnK場(chǎng)。分別隨機(jī)選取2個(gè)后驗(yàn)估計(jì)，見圖4c～4d和圖5c～5d，lnK的估計(jì)均值場(chǎng)與方差場(chǎng)（反映集合的離散型和估計(jì)的不確定性）分別見圖4e、4f和圖5e、5f。對(duì)比污染物分布形態(tài)及高、低值區(qū)域，發(fā)現(xiàn)算例2擬合程度更好，且對(duì)滲透系數(shù)場(chǎng)的刻畫更為精細(xì)。

在使用水頭與濃度數(shù)據(jù)的算例1中（見圖4），其2組后驗(yàn)估計(jì)場(chǎng)在輪廓上基本反映了lnK場(chǎng)的形態(tài)及高、低值區(qū)域，但是在細(xì)節(jié)刻畫上遠(yuǎn)不如使用ERT數(shù)據(jù)的算例2，其反演得到的lnK場(chǎng)與參照?qǐng)鋈杂胁恍〉牟罹啵瑫r(shí)如圖4f所示算例1的估計(jì)方差已然趨近于0，說明在現(xiàn)有基礎(chǔ)上繼續(xù)進(jìn)行迭代對(duì)反演結(jié)果的提升非常有限。在使用ERT數(shù)據(jù)的算例2中（見圖5），兩組后驗(yàn)估計(jì)場(chǎng)與參照?qǐng)龈鼮榻咏?，方差?chǎng)（圖5f）在研究區(qū)左右兩側(cè)較高，其原因是利用溫納裝置獲得的呈倒梯形分布的ERT數(shù)據(jù)無法覆蓋到相關(guān)區(qū)域。由此可知，融合ERT數(shù)據(jù)的ILUES方法可以更加精確地估計(jì)污染場(chǎng)地的滲透系數(shù)場(chǎng)。

圖4 滲透系數(shù)場(chǎng)的反演結(jié)果（算例1）Fig.4 The inversion results of hydraulic conductivity in Case 1

圖5 滲透系數(shù)場(chǎng)的反演結(jié)果（算例2）Fig.5 The inversion results of hydraulic conductivity in Case 2

在此基礎(chǔ)上，將算例1與算例2反演所得ln K場(chǎng)代入地下水流與溶質(zhì)運(yùn)移模型，獲取其末時(shí)刻的污染物濃度分布，并與作為參照的真實(shí)濃度分布進(jìn)行對(duì)比，如圖6。圖6a反映了研究區(qū)末時(shí)刻真實(shí)的污染物濃度分布；圖6b是使用初始隨機(jī)lnK場(chǎng)正演得到的末時(shí)刻濃度分布；圖6c、6d分別為算例1和算例2反演所得lnK場(chǎng)代入地下水流與溶質(zhì)運(yùn)移模型中所獲取的污染物濃度分布情況。經(jīng)過對(duì)比可以得到以下結(jié)論：雖然初始隨機(jī)生成的lnK場(chǎng)與真實(shí)場(chǎng)有較大的差距，但是經(jīng)過算例1與算例2的參數(shù)反演后，其分布基本接近于真實(shí)場(chǎng)，證明了ILUES算法在求解參數(shù)反演問題上的有效性；算例2末時(shí)刻的濃度分布情況與真實(shí)濃度分布更為接近，在細(xì)節(jié)的刻畫上更加精確，進(jìn)一步證明了融合ERT數(shù)據(jù)的ILUES方法的優(yōu)越性。

圖6 ln K正演所得末時(shí)刻濃度分布Fig.6 The final concentration distribution from forward modeling of ln K

4 結(jié)論

基于ILUES方法構(gòu)建融合溶質(zhì)運(yùn)移模型和地球物理模型的地下水污染溯源的數(shù)據(jù)同化框架，并通過數(shù)值算例分析，得到以下結(jié)論：

（1）使用ILUES算法同化水頭與濃度數(shù)據(jù)和ERT數(shù)據(jù)，可以實(shí)現(xiàn)對(duì)各應(yīng)力期污染源源強(qiáng)和滲透系數(shù)場(chǎng)的聯(lián)合反演，并且隨著迭代次數(shù)的增加，反演結(jié)果越來越趨近于真實(shí)值，證明基于ILUES算法的數(shù)據(jù)同化框架對(duì)于求解地下水參數(shù)反演問題的有效性。

（2）在ILUES數(shù)據(jù)同化框架下，設(shè)計(jì)了兩個(gè)數(shù)值算例對(duì)比同化濃度與水頭數(shù)據(jù)（算例1）和ERT數(shù)據(jù)（算例2）在相同條件下對(duì)于反演結(jié)果精度的影響。在對(duì)污染源源強(qiáng)的反演上，算例2結(jié)果的RMSE值在3個(gè)潛在污染源的全部應(yīng)力期均小于算例1。并且在對(duì)滲透系數(shù)場(chǎng)的反演上，算例2對(duì)于lnK場(chǎng)細(xì)節(jié)的刻畫上明顯優(yōu)于算例1。充分說明了融合ERT數(shù)據(jù)的ILUES方法在求解地下水參數(shù)反演問題上的優(yōu)越性。

（3）算例2通過在較少的觀測(cè)時(shí)間點(diǎn)采集信息，得到的反演結(jié)果精度反而優(yōu)于觀測(cè)時(shí)間點(diǎn)較多的算例1，其原因在于使用高密度電法一次可以獲取大量的具有空間連續(xù)性的觀測(cè)信息，而傳統(tǒng)的濃度與水頭數(shù)據(jù)，只能在不連續(xù)的少量觀測(cè)井處獲取。在獲取觀測(cè)信息的數(shù)量與規(guī)模上，結(jié)合地球物理方法的算例2要明顯優(yōu)于使用傳統(tǒng)觀測(cè)方法的算例1。并且由于利用高密度電法等地球物理方法進(jìn)行觀測(cè)的便利性和觀測(cè)數(shù)據(jù)的連續(xù)性，其在實(shí)際場(chǎng)地應(yīng)用中將更有優(yōu)勢(shì)。

（4）本文僅研究了對(duì)于非均質(zhì)滲透系數(shù)場(chǎng)的反演，而在實(shí)際應(yīng)用中，孔隙度、彌散度等水文地質(zhì)參數(shù)也存在著不同程度的非均質(zhì)性，且各參數(shù)之間存在著一定的因果關(guān)系和相關(guān)性。后續(xù)研究計(jì)劃結(jié)合機(jī)器學(xué)習(xí)和深度學(xué)習(xí)等方法同時(shí)反演多種水文地質(zhì)參數(shù)。

作者貢獻(xiàn)聲明：

張瑞城：論文的主要撰寫，數(shù)據(jù)處理與分析。

周念清：指導(dǎo)論文的總體框架，協(xié)助修改與定稿。

江思珉：協(xié)助完善數(shù)值方法的理論部分。

夏學(xué)敏：協(xié)助完善圖表的繪制與處理。