基于復合型進化算法的地下水污染反演模型

黃林顯　劉治政　邢立立　楊麗芝　朱恒華　紀紋龍　張永勇

摘要：污染源位置和污染物排放濃度的快速確定直接決定著地下水污染的有效治理及修復，這屬于地下水反演問題。通過充分分析地下水污染反演問題，耦合地下水流模擬程序MODFLOW、溶質運移模擬程序MT3DMS和優化算法SCE-UA，設計了一種模擬一優化（ S/O）反演模型。通過實例驗證，反演結果表明：提出的網格遍歷CT算法可以自動驗證潛在污染區內所有可能的污染源位置組合方式：與傳統地下水污染反演模型相比，S/O模型不但能夠適用于穩定流條件，而且適用于非穩定流條件;所開發的S/O模型對于多污染源分別在穩定流和非穩定流下的反演均有非常高的精度，能夠準確反演污染源位置及污染物排放濃度。

關鍵詞：地下水污染;污染源位置：SCE-UA優化算法;S/O模型

中圖分類號：P641.2

文獻標志碼：A

doi：10.3969/j.issn. 1000- 1379.2019.04.013

地下水一旦遭受污染，其治理不但投入巨大，而且耗時長[1-3]。地下水污染防治的關鍵是及時確定污染源的位置，從而采取措施切斷污染途徑，并進行針對性修復，避免更大范圍的污染發生。如何快速、準確定位污染源位置和污染物排放濃度，已成為地下水科研領域一個非常重要的研究方向，這屬于地下水反演問題。國內外學者對地下水污染反演問題進行了大量研究：Prakash等[4]通過優化設計監測點網絡布局，利用污染物濃度的變化進行地下水污染源識別：Gorelick等[5]分別將最小二乘法、線性規劃法等與溶質運移模型結合，對比了兩種方法運用于地下水污染物運移特征識別的優缺點：Foddis等[6]設計了一種借助人工神經網絡（ ANNs）的優化反演模型，并將其應用于均質各向同性二維含水層反演中：江思珉等[7]采用Hooke -Jeeves吸引擴散粒子群混合算法、和聲搜索算法等對地下水污染物釋放強度進行了反演：顧文龍等[8]將污染源反演過程轉化為貝葉斯推理過程，并與克里格替代模型結合，提出了一種反演地下水污染源釋放濃度的新思路。此外，自適應模擬退火算法（ASA）[9]、自適應多尺度方法10]、正態轉換一集合卡爾曼濾波法（ NS-EnKF）[11]和蒙特卡洛[12]等方法也被應用于地下水污染反演問題的求解中。

上述方法雖然從不同角度對地下水污染反演問題進行了探討，并取得了很多成果，但也存在一定局限性。如監測點網絡優化布局法需要大量監測井的監測數據，采樣和測試過程耗費大量人力和財力，且只能大致確定污染源的方向，不能快速精確定位：傳統的最小二乘法、線性規劃法等，當研究區水文地質條件比較復雜時，反演過程容易陷入局部搜索，得不到全局最優解，影響反演精度：地質統計學方法的反演精度取決于對研究區的了解程度，當未知變量較多時，計算量大，且容易出現病態矩陣問題：近幾年興起的全局最優解啟發式搜索方法如人工神經網絡法等，需要大量數據進行樣本訓練，如果不能實現地下水流模型、溶質運移模型和優化算法的有效耦合，那么會引起計算過程的復雜化，增加計算負擔，或者只能模擬穩定流，不能模擬非穩定流，限制了其使用范圍。此外，對于一些傳統反演優化方法，需要已知污染源位置才能對排放濃度進行反演，而實際應用中污染源的位置是不確定的，或者只能確定在某一個范圍內，這限制了此類方法的應用。

基于以上分析，本文提出了一種基于全局搜索SCE-UA算法的模擬一優化S/O反演模型，并進行不同案例情況下的反演驗證。該模型采用SCE-UA優化算法，能夠更加有效、快速搜索全局最優解，避免了傳統優化算法容易早熟收斂、陷入局部最優解的弊端，魯棒性好。通過采用FORTRAN語言編寫接口程序，實現了地下水數值模擬程序MODFLOW（地下水流模擬程序）和MT3DMS（溶質運移模擬程序）與優化算法SCE-UA的對接，使得數據交換由傳統的文件讀取改進為內部變量傳遞，有效減小了計算負擔，且能夠適用于多污染源、穩定流和非穩定流等各種復雜情況。同時，提出了一種網格遍歷GT算法，可以在沒有或者只有很少關于污染源位置信息的情況下，通過驗證潛在污染區所有可能污染源的位置組合方式，對污染源進行精確定位。

1 方法原理

S/O反演模型利用MODFLOW和MT3DMS模擬污染物在地下水中的運移過程：SCE -UA優化算法根據觀測井處污染物濃度模擬值和觀測值的標準化差值NE，通過反射、變異和進化對污染源的污染物排放濃度進行反演：GT算法通過所設計的網格遍歷算法快速、有效地搜索研究區內所有可能的污染源位置組合，最終達到準確確定污染源位置、濃度以及污染物排放過程的目的。3個模塊通過FORTRAN程序實現相互之間的鏈接，并通過變量傳遞實現數據交換。S/O反演模型的主要結構及不同模塊之間的鏈接見圖1（其中BAS、RIV和LPF等分別為MODFLOW和MT3DMS的子程序包）。

1.1 地下水控制方程

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

研究中，地下水流模擬程序MODFLOW和溶質運移模擬程序MT3DMS被聯合使用，以模擬求解污染物運移的時空分布狀況。

1.2 SCE-UA算法

SCE-UA算法是一種全局搜索優化算法，最初由美國亞利桑那大學的Duan等在1992年提出。該算法可以有效解決地下水污染物運移強烈非線性特征所帶來的早熟收斂、容易陷入局部最優解等弊端，能夠快速搜索到全局最優解，且穩定可靠，較其他算法具有一定優越性。SCE-UA算法綜合了單純形法、隨機搜索和生物競爭進化等方法的優點，引人種群概念，復合形點在可行域內隨機生成，并根據生物進化規則不斷優化[15-16]。

1.3 遍歷尋優（GT）算法

進行地下水污染源反演時污染源的位置往往是未知的。解決此類問題通常的做法是通過實地調查和查閱有關資料預先劃定一個潛在子區域，確保所有可能的污染源組合落在這個子區域中，見圖2。

本次研究開發了一種可以自動搜索預定義子區域中多個污染源所有可能位置組合的網格遍歷GT算法。在GT算法中，開始單元格分別通過行、列、層號索引變量Rmin、Cmin、Lmin界定預定義子區域的上邊界，結束單元格分別通過行、列、層號索引變量Rmax、Cmax、Lmax界定預定義子區域的下邊界，所有可能的污染源組合通過以上6個變量遍歷。整個搜索過程通過FORTRAN語言編寫的計算機程序進行，并與S/O優化模型鏈接。對于非穩定流問題的反演，GT算法將遍歷每一個應力期內所有污染源位置的組合，并通過判斷目標函數RE值來獲得最優解。反演問題的復雜性隨著污染源個數的增加而增加，例如預定義子區域包含16個單元格，如果只有1個污染源，那么可能的污染源位置為16個;如果污染源的個數是2個，那么可能的污染源位置組合為120個：如果污染源的個數是3個，那么可能的污染源位置組合為560個。

2 模擬一優化模型

2.1 目標函數及標準化差公式

S/O模型在反演過程中通過目標函數RE值調整SCE-UA算法中種群的進化，通過標準化差值NE檢驗反演結果的精確性和魯棒性。

（1）目標函數。目標函數用來判斷監測井實測污染物濃度與S/O模型模擬污染物濃度的吻合程度，目標函數的選擇對反演結果的優劣至關重要，它是SCE-UA算法進化尋優的基礎。目標函數定義為

（2）標準化差公式。為了檢驗S/O反演模型的精度和魯棒性，引入了標準化差公式。假設污染源污染物排放濃度真實值Cact已知，可以通過MODFLOW和MT3DMS計算監測井處的濃度Cobs：假設Cact未知，可以通過監測濃度Cobs和S/O模型反演污染源污染物排放濃度Cide。S/O模型的反演精度可以通過以下標準化差公式進行衡量：

2.2 S/O模型結構

S/O模型可以劃分為模擬模型和優化模型兩部分。模擬模型包含地下水流模擬程序MODFLOW和溶質運移模擬程序MT3DMS.主要用來模擬地下水流和污染物的有關運移特征：優化模型主要通過SCE -UA優化算法產生種群樣本點（污染源污染物濃度值），并根據RE值進行變異、反射和進化樣本點值，在全局范圍內搜索可能的污染源位置和污染物排放濃度，最終達到快速、準確鎖定污染源位置和污染物濃度的目的。需要注意的是當RE值為0時，說明反演的污染源位置和污染物濃度值與真實情況完全一致，但由于實際操作過程中模擬誤差和觀測誤差的存在，因此這種情況很難發生。GT遍歷算法通過自動搜索所有可能的污染源位置.提高了S/O優化模型的反演效率和準確度。地下水模擬模型、優化模型和GT算法通過FORTRAN語言編寫的接口程序內部鏈接，使得數據交換由傳統的文件讀取轉化為內部變量傳遞，極大提高了計算效率。

地下水模擬模型和優化模型的鏈接：①SCE -UA算法和接口程序重設MODFLOW和MT3DMS的相關輸入文件，如WEL、SSM和BTN文件等。②MODFLOW和MT3DMS模擬計算監測井污染物濃度，并與監測濃度相比較，計算RE值，如果RE值小于收斂指標，認為反演值與實際情況一致，則進行下一個應力期的反演：如果RE值大于收斂指標，SCE -UA算法根據RE值，通過變異、反射和進化對污染源污染物濃度值優化，繼續重設MODFLOW和MT3DMS的相關輸入文件。在一個應力期內當所有可能的污染源位置都通過GT算法被驗證后，反演過程將會移向下一個應力期。S/O模型鏈接過程見圖3。

3 案例研究

S/O反演模型的反演效果分別通過以下兩個案例進行驗證：案例1為穩定流條件，有兩個污染源但位置未知，只知道可能存在的范圍，通過4個監測井的監測濃度值反演污染源的位置和污染物濃度排放值：案例2更接近于實際情況，為包含3個應力期的非穩定流條件，有兩個位置未知的污染源，且污染源污染物排放濃度在每個應力期均不相同，通過6個監測井的監測濃度值反演不同應力期污染源的位置和污染物排放濃度。案例研究參考RAJ等[17]進行地下水污染源反演驗證時所設計的模擬河間地塊地下水流動情況的實例模型（見圖4）。本次案例研究的主要目的是評價S/O模型的反演效果，因此在地下水數值建模時對研究區水文地質條件進行了一定程度的概化。數值模型在橫向上剖分成規格為100 mxl00 m的正方形網格，縱向上剖分為1層。模型區東、西邊界分別為給定水頭邊界（東邊界為88 m，西邊界為100 m），南、北邊界分別為隔水邊界。S/O模型相關參數見表1。

SCE-UA算法進化代數為10時的反演結果見表2。需要注意是，由于276個位置組合需要被搜索，因此在實際反演過程中產生了很多結果，表2僅僅根據NE值選取了幾組具有代表性的結果進行展示。從表2可以看出，真實污染源的位置可以從276個可能的位置組合中被準確搜索出來，如結果排序1-3，且反演結果最好時兩個污染源污染物排放濃度分別為47.950mg/L和35.805 mg/L，與真實值48 mg/L和36 mg/L較為接近，NE值為0.323%，能夠滿足精度要求。

從表2可以看出，當反演位置正確的時候，反演的濃度值均較接近真實值，NE值為0.323% -1 .945%.基本滿足精度要求：但當反演位置不準確時，反演的污染物濃度值與真實值差別很大，NE值為8.236% -90.765%.不能滿足精度要求，因此可以看出污染源位置的準確確定是污染物濃度反演成功的前提條件。此外，對進化代數對反演結果的影響進行了測試，從表3可以看出，當進化代數分別為10、20和50時，反演結果與真實值均非常接近，NE值為0.008% - 0.323%，完全能夠滿足精度要求。同時，隨著進化代數的增加，NE值逐漸減小，說明反演精度提高。但需要注意的是，進化代數的增加會帶來很大的計算負擔，反演時間會成倍增加。

（2）案例2。現實中污染物的運移往往是在非穩定流條件下進行的，并且不同應力期的排放濃度是變化的。在案例2中，對具有3個應力期（應力期為la）的非穩定流情形進行驗證。污染源的真實位置分別位于A1（行=3、列=2）和A2（行=7、列=3），且假設未知，可能的位置位于圖6陰影區域。6個監測井分別位于01（行=2、列=4）、02（行=3、列=4）、O3（行=4、列=5）、O4（行=5、列=5）、05（行=6、列=4）、06（行=7、列=4）。每個應力期的真實濃度和監測濃度見表4。

由于案例2有3個應力期，因此每個應力期都需要搜索276個可能的位置組合。S/O反演模型在處理非穩定流問題時，當前應力期反演的進行要基于前面應力期的反演結果：反演第一個應力期并獲得污染源位置和污染物排放濃度的最優解：第二個應力期的反演，首先運行第一個應力期反演的最優解，然后再反演第二個應力期并獲得最優解;第三個應力期以此類推。所以，本應力期的反演結果取決于前面應力期反演結果的準確性，NE值也會因前面應力期反演誤差的累積而不斷增大。

當進化代數為10時，每個應力期的最優反演值見表5。從表5可以看出，第一個應力期NE值是2.91%，反演結果尚在可以接受的精度范圍內，但隨著應力期的增加，NE值在第三個應力期增大為58.00%（第二個應力期由于A2點污染物排放濃度為0，因此NE值無法計算），且反演的位置與真實位置不一致。為了獲得更準確的反演結果，把進化代數分別增加為20和50，反演結果見表6。

從表6可以看出，當進化代數為20時，只有前兩個應力期能夠獲得可以接受的結果，到第三個應力期，反演結果與真實值相差較大;當進化代數為50時，全部3個應力期污染源位置和污染物排放濃度反演值都與真實值較為接近，能夠獲得令人滿意的結果。由此可以得出：通過增加進化代數.S/O反演模型的反演精度可以得到提高;隨著非穩定流應力期的增加，需要相應增加進化代數以獲得更加準確的反演結果。

4 結語

在已知監測井監測濃度的情況下推求污染源位置和污染物排放濃度是典型的地下水數值模擬反演問題。反演過程中可轉化成決策變量為污染源位置和污染物排放濃度的最優化問題進行求解。將優化算法SCE - UA和地下水數值模擬程序MODFLOW和MT3DMS結合起來，建立了S/O優化搜索模型。

（1）通過接口程序，實現了SCE -UA優化算法和MODFLOW、MT3DMS的鏈接，使得數據交換由文件讀取改進為內部變量傳遞，不但大大提高了S/O反演模型的計算效率，而且使其同時具備在穩定流和非穩定流條件下反演的能力。

（2）穩定流案例研究表明，S/O優化模型對于雙污染源、污染源位置未知情況下的反演，能一致、高效地收斂到全局最優解，收斂速度較快、穩定性好，在進化代數較少的情況下就可以實現對污染源位置和污染物排放濃度的準確反演。

（3）對于非穩定流案例研究表明，S/O優化模型的反演誤差會隨著應力期的增加不斷累積，但通過增加進化代數，同樣可以獲得令人滿意的反演結果，能夠對污染源位置和污染物排放濃度進行精確反演。

（4）隨著應力期和進化代數的增加，計算負擔不斷加大，計算時間甚至是以天或月計，下一步要重點研究并行算法在S/O反演模型中的應用，以進一步提高計算效率，并采用更加接近實際情況的實例模型進行驗證，使S/O反演模型更加具有實用性。

參考文獻：

[1]KANEL S R，MALLA G B，CHOI H.Modeling and Study ofthe Mechanism of Mobilization of Arsenic Contamination inthe Croundwater of Nepal in South Asia[J].Clean Technol-ogies and Environmental Policy， 2013， 15（6）：1077-1082.

[2]DUCCI D， SELLERINO M.Vulnerability Mapping ofCroundwater Contamination Based on 3D LithostratigraphicalModels of Porous Aquifers[J].Science of the Total Environ-ment， 2013， 447： 315-322.

[3]江思珉，王佩，施小清，等，地下水污染源反演的Hooke -Jeeves吸引擴散粒子群混合算法[J].吉林大學學報：地球科學版，2012，42（6）：1866-1872.

[4] PRAKASH 0， DATTA B.Optimal Monitoring NetworkDesign for Efficient Identification of Unknown CroundwaterPollution Sources[ J]. Intemational Journal of Ceomate，2013， 23（ 10）： 2031-2049.

[5] CORELICK S M， EVANS B，REMSON I.IdentifyingSources of Croundwater Pollution： an Optimization Approach[J]. Water Resources Research， 1983， 19（3）： 779-790.

[6] FODDIS M L，ACKERER P，MONTISCI A， et al.ANN-based Approach for the Estimation of Aquifer PollutantSource Behavior[J].Water Science and Technology： WaterSupply， 2015， 15（6）：1285-1294.

[7] 江思珉，蔡奕，王敏，等，基于和聲搜索算法的地下水污染源與未知含水層參數的同步反演研究[J].水利學報，2012，43（ 12）：1470-1477.

[8] 顧文龍，盧文喜，張宇，等，基于貝葉斯推理與改進的MCMC方法反演地下水污染源釋放歷史[J].水利學報，2016，47（6）：772-779.

[9] JHA M， DATTA B.Three-Dimensional Croundwater Con-tamination Source Identification Using Adaptive SimulatedAnnealing[J].Journal of Hydrologic Engineering， 2012， 18（3）：307-317.

[10]MAJDALAM S，ACKERER P.Identification of CroundwaterParameters Using an Adaptive Multiscale Method[J].Ground-water， 2011， 49（4）： 548-559.

[11]LL，ZHOU H， HENDRICKS F H J，et al.CroundwaterFlow Inverse Modeling in Non-Multigaussian Media： Per-formance Assessment of the Normal - Score EnsembleKalman Filter[J].Hydrology and Earth System Sciences，2012， 16（2）：573-590.

[12]HOSSEINI A H， DEUTSCH C V， MENDOZA C A， et al.Inverse Modeling for Characterization of Uncertainty inTransport Parameters Under Uncertainty of Source Ceometryin Heterogeneous Aquifers[J].Joumal of Hydrology， 2011（3）：402-416.

[13]HARBAUCH A W. Modflow-2005， the US Ceological SurveyModular Croundwater Model： the Croundwater Flow Process[R]. Reston， VA： US Ceological Survey， 2005： 1-21.

[14]ZHENC C M， WANC P P.MT3DMS：A Modular Three-Diruensional Multi-Species Transport Model for Simulationof Advection， Dispersion and Chemical Reactions of Con-taminants in Croundwater Systems[ R]. Alabama：AlabamaUniv University， 1999：3-15.

[15]DUAN Q Y， GUPTA V K， SOROOSHIAN S.Suffled Com-plex Evolution Approach for Effective and Efficient ClobalMinimization[J].Journal of Optimization Theory and Ap-plications， 1993， 76（3）： 501-521.

[16] SOROOSHIAN S， DUAN Q， CUPTA V K. Calibration ofRainfall - Runoff Models： Application of ClobalOptimization to the Sacramento Soil Moisture AccountingModeI[J]. Water Resources Research， 1993， 29 （ 4） ：1185-1194.

[17] RAJ Mohan Singh， BITHIN Datta. Croundwater PoUutionSource Identification and Simultaneous Parameter EstimationUsing Pattem Matching by Artificial Neural Network [ J] . En-vironmental Forensics， 2004， 5（ 3） ： 143-153.