王 敬，韓 忠，梁 浩，劉樂軍，林艷竹，袁星芳

(1.山東省第六地質(zhì)礦產(chǎn)勘查院，山東 威海 264209；2.山東省地質(zhì)環(huán)境監(jiān)測總站，濟南 250014；3.中國地質(zhì)環(huán)境監(jiān)測院(地質(zhì)災(zāi)害技術(shù)指導(dǎo)中心)，北京 100081)

0 引 言

地下水是非常重要的飲用水源[1]。隨著經(jīng)濟的快速發(fā)展，我國地下水污染問題已經(jīng)非常嚴(yán)重[2]。在地下水污染發(fā)生后，如何準(zhǔn)確、快速通過監(jiān)測點污染物濃度識別出污染源排放強度及其動態(tài)變化就變得極為重要。此類問題屬于地下水污染反演識別問題，并已經(jīng)成為地下水科研領(lǐng)域一個非常重要的研究方向[3]。

面對這一問題，過去幾十年間國內(nèi)外水文地質(zhì)學(xué)家探討了各種方法來進行地下水污染源反演識別：Wagner開發(fā)了一種可以同時反演水文地質(zhì)參數(shù)和污染源特性的非線性極大似然優(yōu)化算法[4]；Alapati等設(shè)計了一種基于最小二乘法的一維地下水污染源反演方法[5]；Foddis等利用人工神經(jīng)網(wǎng)絡(luò)(ANNs)算法對二維含水層進行了污染源反演識別[6]；Jha等采用自適應(yīng)模擬退火算法進行了地下水污染源特征的反演[7]；Machiwal 結(jié)合多元統(tǒng)計分析方法和基于GIS的地質(zhì)模型，進行地下水污染的識別[8]。在國內(nèi)，江思珉等采用和聲搜索算法[9]、單純型模擬退火算法[10]和Hooke-Jeeves吸引擴散粒子群混合算法[11]等，進行地下水污染源相關(guān)特性的反演識別；顧文龍等利用改進的MCMC方法與貝葉斯推理進行地下水污染源釋放歷史反演[12]；肖傳寧等開發(fā)了一種基于徑向基函數(shù)模型的地下水污染反演優(yōu)化模型[13]。然而，由于地下水含水層的復(fù)雜性、地下水溶質(zhì)運移的非線性及人類活動影響等問題的存在，導(dǎo)致設(shè)計一種高效且可行的地下水污染源反演模型并非易事，如前面提到的各種優(yōu)化算法均或多或少存在一些弊端：局部尋優(yōu)算法(最小二乘法和線性規(guī)劃法等)往往由于容易陷入局部最優(yōu)解，而無法獲得全局最優(yōu)解；全局尋優(yōu)算法(遺傳算法、模擬退化算法等)能夠處理不連續(xù)和不可微的方程，并且能夠?qū)ふ业饺肿顑?yōu)解，但此類方法往往涉及到反復(fù)的模型調(diào)用，并且求解效率會隨著數(shù)值模型復(fù)雜程度的增大而降低[14]。

基于以上分析，亟待提出一種能夠高效獲得全局最優(yōu)解的地下水污染源反演優(yōu)化算法。由Duan等提出的SCE-UA算法結(jié)合了隨機搜索、生物競爭進化和單純型法等方法的優(yōu)點，是一種可以快速、有效、一致地搜索到全局最優(yōu)解的進化算法[15]。Kuczera等分別比較了SCE-UA算法、MSX算法(擬牛頓法和單純型法)和遺傳算法在進行地表水模型參數(shù)反演識別中的搜索效率，結(jié)果證明SCE-UA算法收斂效果更好、魯棒性更強[16]。然而，在利用SCE-UA算法進行水文地質(zhì)參數(shù)反演，特別是地下水污染強度及其動態(tài)變化反演方面的應(yīng)用還很少。為此，開發(fā)了一個耦合SCE-UA優(yōu)化算法、地下水水流和溶質(zhì)模型，適用于多污染源穩(wěn)定流和非穩(wěn)定流、定濃度和不定濃度排放等各種復(fù)雜條件下地下水污染源反演識別的優(yōu)化模型，并進行了不同案例情況下的反演驗證。

1 研究方法

1.1 SCE-UA優(yōu)化算法

最初由美國亞利桑那大學(xué)的Duan等在1992年進行降雨徑流模型參數(shù)優(yōu)化問題時提出的SCE-UA算法是一種全局優(yōu)化算法[17]。該算法能夠有效處理含水層各向異性及地下水溶質(zhì)運移強烈非線性所帶來的容易停留在局部最優(yōu)解、早熟收斂等弊端，可以高效、快速搜索到全局最優(yōu)解，且相較其他類似優(yōu)化算法效率更高、魯棒性更強。SCE-UA算法綜合了隨機搜索、生物競爭進化和單純型法等方法的優(yōu)點，引入了種群的概念，在可行域內(nèi)隨機生成復(fù)合型點，并依據(jù)生物進化規(guī)則不斷進行優(yōu)化。

SCE-UA算法的具體實現(xiàn)步驟如下：

(1)算法初始化。設(shè)定維數(shù)n，參與進化的復(fù)合型p個(p≥1)和每個復(fù)合型包含的頂點數(shù)目為m(m≥n+1)，計算樣本點數(shù)目s=pm。

(2)產(chǎn)生樣本點。在可行域內(nèi)隨機產(chǎn)生s個樣本點x1，…，xs，分別計算每一點xi處的目標(biāo)函數(shù)值fi。

(3)樣本點排序。把生成的s個樣本點按目標(biāo)函數(shù)值f升序排列。

(4)劃分復(fù)合型群體。將s個樣本點劃分為p個復(fù)合型，A1…AK，k=1,2…p。復(fù)合型按照如下標(biāo)準(zhǔn)劃分，按照排序第一個復(fù)合型包含p(k-1)+1，k=1,2…m位置處的樣本點，第二個復(fù)合型包含p(k-1)+2，k=1,2,…,m位置處的樣本點。以此類推。

(5)復(fù)合型進化。按復(fù)合型進化算法(CCE)分別進化各個復(fù)合型，在后面將詳細(xì)敘述。

(6)復(fù)合型混合。把進化后的每個復(fù)合型的所有頂點組合成新的點集，再次按目標(biāo)函數(shù)值f升序排列。

(7)收斂性判斷。如果滿足收斂條件或者累計進化代數(shù)達到總進化代數(shù)，終止循環(huán)；否則返回步驟(4)。若累計進化代數(shù)達到總進化代數(shù)，但未得到優(yōu)化結(jié)果，終止循環(huán)后調(diào)整參數(shù)重新進行計算。

SCE-UA方法一個關(guān)鍵的組成部分是步驟(5)中提到的復(fù)合型競爭進化法則。其基本步驟如下：

(1)初始化，選取參數(shù)q，a，b，其中2≤q≤m，a≥1，b≥1。

(3)根據(jù)權(quán)重在復(fù)合型Ak中選擇q個父個體(因此q即是每個子復(fù)合型中樣本點的個數(shù))，u1-uq記作集合B，并將父個體在Ak中的位置記作集合L。

(4)通過如下幾步產(chǎn)生子個體：

②計算新個體r=2g-uq(反射)。

③如果r在可行解區(qū)域內(nèi)，則計算該點的目標(biāo)函數(shù)值fr并繼續(xù)進入第④步。否則在可行域范圍內(nèi)隨機產(chǎn)生新個體z，使r=z(突變)。

④如果fr

⑤如果fc

⑥重復(fù)①～⑤步a次。

1.2 地下水控制方程

(1)地下水流動方程。根據(jù)Darcy定律以及質(zhì)量守恒，不考慮水的密度變化條件下，三維空間孔隙介質(zhì)中地下水流動的偏微分方程[18]：

式中：Kx、Ky、Kz分別為滲透系數(shù)在x、y、z方向上的分量；W為源匯項，用以表示流進匯或來自源的水量；h為含水層水位標(biāo)高；Ss為含水介質(zhì)的貯水率；t為時間。

(2)污染物遷移的數(shù)學(xué)表達式[19]：

在本次優(yōu)化模型的設(shè)計中，地下水流數(shù)值模型MODFLOW和溶質(zhì)運移數(shù)值模型MT3DMS用來模擬污染物在地下水中運移的時空分布特征。由于MODFLOW和MT3DMS具有模塊化設(shè)計的特點，因此能夠?qū)⑵浞浅７奖愕鸟詈系絻?yōu)化模型中，并且能夠針對不同的模塊進行單獨處理。

1.3 反演優(yōu)化模型

1.3.1 數(shù)值模型與優(yōu)化模型的耦合

如何實現(xiàn)地下水?dāng)?shù)值模型(MODFLOW、MT3DMS)和優(yōu)化模型的有效耦合是決定地下水污染強度優(yōu)化識別能否高效執(zhí)行的關(guān)鍵因素。傳統(tǒng)的優(yōu)化模型往往是將MODFLOW和MT3DMS作為單獨的外部可執(zhí)行文件進行調(diào)用，并且數(shù)據(jù)傳遞是通過文件讀寫進行。由于符合精度要求的反演結(jié)果往往需要成百上千次調(diào)用數(shù)值模型，通過調(diào)用外部可執(zhí)行文件的方式會造成優(yōu)化模型的反演效率非常低下。特別是當(dāng)數(shù)值模型模擬范圍大、網(wǎng)格剖分精細(xì)和反演問題復(fù)雜時，反演識別的時間通常需要幾天甚至十幾天，這種優(yōu)化識別時間是無法忍受的。為了解決這個問題，本次研究將SCE-UA程序改寫為主程序，將MODFLOW和MT3DMS分別作為兩個子程序，在此基礎(chǔ)上編寫了進行子程序調(diào)用的接口程序。優(yōu)化過程中，當(dāng)主程序需要獲得地下水流和污染物的有關(guān)運移特征時，則通過接口程序?qū)蓚€子程序進行調(diào)用(如圖1所示)，以此模擬求解污染物運移的時空分布狀況，并將模擬結(jié)果通過變量返回到主程序中。由于實現(xiàn)了數(shù)值模型和優(yōu)化模型的內(nèi)部耦合及數(shù)據(jù)交換傳遞由文件讀寫改進為內(nèi)部變量傳遞，可以顯著提高優(yōu)化模型的執(zhí)行效率。

圖1 數(shù)值模型與優(yōu)化模型的耦合

1.3.2 目標(biāo)函數(shù)

反演識別過程中的決策變量是地下水污染強度q，而狀態(tài)變量則是質(zhì)量濃度C。目標(biāo)函數(shù)E定義為：

(3)

圖2 優(yōu)化模型求解流程

1.3.3 非穩(wěn)定流反演問題

本次優(yōu)化模型對非穩(wěn)定流數(shù)值模型的反演采用如下策略，即先反演識別第一個應(yīng)力期污染源強度信息，然后對所有的反演結(jié)果進行排序，用第一個應(yīng)力期污染源排放濃度的最優(yōu)反演結(jié)果來反演第二個應(yīng)力期污染源排放濃度；然后在第二個應(yīng)力期最優(yōu)解的基礎(chǔ)上反演第三個應(yīng)力期，依次類推，直到所有的應(yīng)力期反演識別結(jié)束。此種設(shè)計能夠確保每一個應(yīng)力期均使用優(yōu)化模型反演的最優(yōu)解，可以充分保證反演的精度。但是需要指出的是由于當(dāng)前應(yīng)力期的反演是基于上一應(yīng)力期的反演結(jié)果，因此上一應(yīng)力期的反演誤差也會累積到當(dāng)前應(yīng)力期，這樣造成反演誤差會隨著應(yīng)力期的增多不斷被放大。為了減少誤差的累積，對于非穩(wěn)定流問題的反演需要通過相應(yīng)增多反演代數(shù)或者種群個數(shù)來提高每一應(yīng)力期的反演精度。

2 算例研究

如圖3所示，研究區(qū)為一二維均質(zhì)各向同性的矩形承壓含水層(100 m×100 m)，用邊長為10 m的正方形網(wǎng)格將整個研究區(qū)域剖分為10行10列的有限差分網(wǎng)格(案例模型參考Singh等的污染物反演模型)[20]。含水層厚度為10 m，上下邊界為隔水邊界，左右邊界為水頭邊界(左邊界水頭由上到下線性變化，上側(cè)水頭為9.7m，下側(cè)水頭為0.4m；右邊界水頭為7 m)。含水介質(zhì)滲透系數(shù)為1 m/d，孔隙度為0.3，彌散度為20 m；反演時污染物排放濃度變化范圍設(shè)定介于0～100 mg/L之間。由于本次案例研究的主要目的是評價優(yōu)化模型的反演效果，因此在進行地下水?dāng)?shù)值建模時對研究區(qū)水文地質(zhì)條件進行了一定程度的概化。

表1 不同進化代數(shù)時反演結(jié)果

通過表1可以看出，所開發(fā)的優(yōu)化模型可以迅速反演識別出污染物排放濃度：在進化代數(shù)為10代時，反演誤差(|反演值-真實值|/真實×100)為3.13，基本可以滿足精度的需要；隨著進化代數(shù)的增加，反演誤差值逐漸變小，說明反演精度提高；在進化代數(shù)為50代時，可以精確的反演出排放濃度。但需要指出的是，隨著進化代數(shù)的增加會帶來很大的計算負(fù)擔(dān)，造成反演時間增加，因此反演過程中要很好的在反演精度和反演效率之間做出平衡。

在進化代數(shù)為300代時，反演濃度隨代數(shù)增加的變化趨勢如圖4所示。通過圖4可以看出優(yōu)化模型濃度反演最優(yōu)值基本在15代后已達到穩(wěn)定；濃度反演平均值一開始偏離真實值，但在50代左右基本穩(wěn)定；濃度反演最差值一開始也偏離真實值，但在100代左右也可以達到穩(wěn)定狀態(tài)。說明優(yōu)化模型已經(jīng)求解到最優(yōu)解，并且與其他相似優(yōu)化模型相比，該SCE-UA優(yōu)化模型求解的收斂速度也較快。

圖4 反演濃度隨代數(shù)增加的變化趨勢

圖5 案例2含水層示意圖(單位：m)

(2)案例2。為一非穩(wěn)定流模型(共有12個應(yīng)力期)，有兩口井污染井1和污染井2分別位于模型(3, 2)和(5, 5)處，并在12個應(yīng)力期分別以30和80 mg/L的定濃度向地下水中排放污染物；在其下游(3, 8)，(4, 6)，(6, 9)和(7, 8)處有四口監(jiān)測井，以四口監(jiān)測井的污染物監(jiān)測濃度，利用優(yōu)化模型反演識別污染源處12個應(yīng)力期的污染物排放強度。

反演識別結(jié)果如圖6、圖7所示。從圖6可以看出，當(dāng)進化代數(shù)為20代時，兩口污染井的反演結(jié)果都出現(xiàn)比較劇烈的震蕩[特別是(3, 2)處井]，這主要是因為對于非穩(wěn)定流的模擬，下一應(yīng)力期的模擬需要用到上一應(yīng)力期的反演結(jié)果，因此會造成反演誤差的累積，使得模擬結(jié)果變差；隨著進化代數(shù)增加到50代，雖然反演結(jié)果在某些應(yīng)力期仍然有輕微的震蕩，但已基本能滿足反演精度的需要；當(dāng)進化代數(shù)增大到100代，濃度反演值與真實值已非常接近，能夠取得滿意的反演結(jié)果。由此也可以看出隨著非穩(wěn)定流應(yīng)力期的增加，需要相應(yīng)增加進化代數(shù)以獲得更加準(zhǔn)確的反演結(jié)果。

圖6 不同進化代數(shù)條件下污染井1(3, 2)處反演識別結(jié)果

圖7 不同進化代數(shù)條件下污染井2(5, 5)處反演識別結(jié)果

圖8 案例3含水層示意圖(單位：m)

(3)案例3。案例3更接近實際情況，兩口污染井污染物的排放濃度不是定濃度排放，而是隨著時間的變化而變化，其中旱季排放的少，雨季排放的多(如圖9所示)；同時與案例2相比，在(4, 6)和(9, 9)處分別增加了兩口監(jiān)測井；通過六口監(jiān)測井的污染濃度監(jiān)測值，反演兩口污染井在不同應(yīng)力期污染物排放強度的變化。

圖9 進化100代時(3, 2)和(5, 5)處污染井真實值與反演值

通過圖9可以看出，進化代數(shù)為100代時，優(yōu)化模型可以非常好的反演識別不同應(yīng)力期污染源的排放強度，反演值與真實值都較為接近，能夠滿足反演精度的要求。說明通過增加進化代數(shù)可以有效減少反演誤差在不同應(yīng)力期之間的累積；優(yōu)化模型對于多污染源不定濃度排放時污染源強度的反演同樣可以取得令人滿意的結(jié)果。

通過以上3個不同污染物排放情景實例的反演，進一步說明在足夠進化代數(shù)前提下，所提出的基于SCE-UA算法的地下水污染強度反演優(yōu)化模型均能夠獲得令人滿意的結(jié)果。

3 結(jié) 論

通過污染物濃度監(jiān)測數(shù)據(jù)推求污染源排放強度是典型的求逆問題。本次研究將優(yōu)化算法和地下水?dāng)?shù)值模擬程序結(jié)合起來，建立了一種優(yōu)化搜索模型，并進行了不同案例情況下的反演驗證。

(1)研究將優(yōu)化算法SCE-UA和地下水?dāng)?shù)值模擬程序MODFLOW和MT3DMS結(jié)合起來，大大豐富了搜索行為，提高了反演識別的能力和效率。

(2)通過案例1的測試，在進化代數(shù)較少的情況下優(yōu)化模型即可高效、準(zhǔn)確的搜索污染源污染物排放強度；且反演結(jié)果的準(zhǔn)確度隨著進化代數(shù)的增大而增加。但需要注意的是，進化代數(shù)的增加同樣會帶來很大的計算負(fù)擔(dān)，反演時間亦會相應(yīng)增加。

(3)對于非穩(wěn)定流的測試，由于不同應(yīng)力期之間反演誤差的累積，反演精度會隨著應(yīng)力期的增多而減小，需要相應(yīng)增加進化代數(shù)以獲得更加準(zhǔn)確的反演結(jié)果。

(4)在足夠進化代數(shù)的前提下，優(yōu)化模型可以非常好的反演識別不同應(yīng)力期多污染源排放強度的變化。

□