地下水污染強度反演識別的SCE-UA算法

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

(1.山東省第六地質礦產勘查院，山東 威海 264209；2.山東省地質環境監測總站，濟南 250014；3.中國地質環境監測院(地質災害技術指導中心)，北京 100081)

0 引 言

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

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

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

1 研究方法

1.1 SCE-UA優化算法

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

SCE-UA算法的具體實現步驟如下：

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

(2)產生樣本點。在可行域內隨機產生s個樣本點x1，…，xs，分別計算每一點xi處的目標函數值fi。

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

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

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

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

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

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

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

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

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

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

③如果r在可行解區域內，則計算該點的目標函數值fr并繼續進入第④步。否則在可行域范圍內隨機產生新個體z，使r=z(突變)。

④如果fr

⑤如果fc

⑥重復①～⑤步a次。

1.2 地下水控制方程

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

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

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

在本次優化模型的設計中，地下水流數值模型MODFLOW和溶質運移數值模型MT3DMS用來模擬污染物在地下水中運移的時空分布特征。由于MODFLOW和MT3DMS具有模塊化設計的特點，因此能夠將其非常方便的耦合到優化模型中，并且能夠針對不同的模塊進行單獨處理。

1.3 反演優化模型

1.3.1 數值模型與優化模型的耦合

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

圖1 數值模型與優化模型的耦合

1.3.2 目標函數

反演識別過程中的決策變量是地下水污染強度q，而狀態變量則是質量濃度C。目標函數E定義為：

(3)

圖2 優化模型求解流程

1.3.3 非穩定流反演問題

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

2 算例研究

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

表1 不同進化代數時反演結果

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結 論

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

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

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

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

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

□