， ，，，

(1.南京大學地球科學與工程學院，江蘇 南京 210023；2.南京水利科學研究院， 江蘇 南京 210029)

隨著工業(yè)化的發(fā)展加速，地下水污染問題日益突出，極大地破壞著生態(tài)環(huán)境并威脅著民眾的健康[1-3]。根據污染物種類和含水層類型的不同，有多種不同的地下水污染治理與修復方法[4]，其中基于水力控制的污染物抽出-處理(pump and treat, PAT)技術最為典型。其主要優(yōu)點在于簡便易行，可對地下水污染事件做出快速反應，尤其適用于期望短期內迅速降低污染水平的地下水水源地含水層系統(tǒng)，缺點在于修復工程復雜、且運行成本高昂。因此，有效降低PAT系統(tǒng)的運行成本對于該技術的推廣應用具有重要意義。自20世紀80年代以來，地下水模擬-優(yōu)化模型就應用于確定地下水PAT系統(tǒng)的最優(yōu)策略[5]。然而，與國內外廣泛存在的模塊化地下水模擬軟件相比，地下水模擬-優(yōu)化模型方面的研究雖然很多[6-9]，但模塊化的模擬-優(yōu)化軟件相對很少，這種狀況極大地限制了地下水模擬-優(yōu)化模型在實際場地中地下水污染治理的優(yōu)化實踐。

MGO(Modular Groundwater Optimizer)[10]是一款模塊化的地下水模擬優(yōu)化軟件，它通過耦合地下水水流模型MODFLOW[11]和污染物運移模型MT3DMS[12]結合通用管理目標建立地下水PAT系統(tǒng)的模擬-優(yōu)化模型，并可根據需要選擇不同優(yōu)化算法對模型進行求解。本文基于MGO軟件，運用遺傳算法對美國猶他州的某地下水污染場地進行抽出-處理系統(tǒng)的優(yōu)化研究，確定抽水井群的單井最優(yōu)動態(tài)抽水量，實現以最低的系統(tǒng)運行成本達到修復含水層的最優(yōu)目標。

1 MGO軟件簡介

MGO(Modular Groundwater Optimizer)是美國阿拉巴馬大學教授Zheng和Wang基于地下水水流模型MODFLOW和污染物運移模型MT3DMS 開發(fā)出來的一款針對地下水PAT系統(tǒng)優(yōu)化的通用模塊化優(yōu)化管理軟件。利用MGO不僅能建立完全耦合地下水模擬模型(包括水流模型和運移模型)的嵌入式模擬-優(yōu)化模型，而且對于特別復雜的地下水優(yōu)化問題還能結合實際情況建立新的基于替代模型(或響應函數)的地下水模擬-優(yōu)化模型。無論是嵌入式優(yōu)化模型還是基于響應函數的優(yōu)化模型，MGO都可選擇遺傳算法(GA)、模擬退火(SA)和禁忌搜索(TS)等多種進化算法中的任一算法來求解，從而保證能找到復雜地下水非線性模擬-優(yōu)化問題的最優(yōu)解。

MGO既能應用于地下水系統(tǒng)的單一水量/水質優(yōu)化管理，也能應用于水量與水質的聯(lián)合優(yōu)化管理，尤其適用于基于水力控制的地下水污染治理的PAT系統(tǒng)中的抽注水方案優(yōu)化管理。同時，MGO采用模塊化設計(見圖1)，能直接利用地下水模擬模型MODFLOW和MT3DMS中的輸入文件，大大簡化了軟件的使用，便于推廣應用。

圖1 MGO代碼設計流程圖

2 基于MGO的地下水PAT系統(tǒng)模擬-優(yōu)化模型

為了能同時考慮各種不同的水量/水質優(yōu)化管理問題，基于MGO設計的地下水PAT系統(tǒng)優(yōu)化模型可表示為以下通用數學形式[10]：

目標函數：

(1)

約束條件：

(2)

Qmin≤Qi≤Qmax

(3)

hmin≤hm≤hmax

(4)

(5)

cmin≤cm≤cmax

(6)

(7)

上述優(yōu)化模型中，式(1)為通用目標函數，J是管理目標，可以是總成本、總抽水量或含水層中污染物的總去除量等。用戶可以根據實際情況選取目標函數右端項中的一項或若干項(亦可表示為最大化)。Qi是某個管理期內的第i個井在單位時間內的抽/注水量，正值代表注水、負值代表抽水。對于單個管理期優(yōu)化問題，任何井的抽/注水量均為是常數，只需用一個參數表示；而對于多管理期的優(yōu)化問題，由于任何一個井在不同管理期的抽/注水量可以不同，則需要多個參數來表示井在不同管理期的抽/注水量。Mi第i井的溶質移除量。F(q, h, c)用戶自定義成本函數，可以依賴于流量q,水頭h,濃度c。N是待優(yōu)化參數的總數目，下同。yi是個二進制數，當考慮參數i活躍時，yi=1；否則yi=0。di是第i參數相關的井徑深度。Δti是參數i的抽/注水持續(xù)時長期。α1每個井的固定資產成本。α2是井筒單位深度的安裝、打鉆成本。α3是單位體積的抽水和處理成本。α4是單位質量污染物的處理成本。α5是外部用戶指定的成本系數。

式(2)-(7)為通用約束條件。式(2)要求任何時間段內的決策井數目不得超過一個給定的約束值(NW)。 式(3)是對任一管理期內單位時間的抽水量約束，Qmin、Qmax分別是對應井的最小、最大單位抽水量。式(4)是水頭約束，hm是監(jiān)測點位置的水頭，hmin、hmax分別是其最小、最大約束水頭。式(5)是水力約束，Δhmin是外部監(jiān)測點與內部監(jiān)測點之間所允許的最小水頭差。式(6)是濃度約束，cm是任一監(jiān)測點的濃度，cmin、cmax是該監(jiān)測點的最小、最大濃度允許值。式(7)可作為均衡約束，Qm是與從井I1到井I2的抽/注水量之和成比例的抽水量，其中A、B為比例系數。同樣地，用戶可以根據優(yōu)化模型中目標函數的組成來確定和選擇必要的約束條件。

優(yōu)化模型中涉及到兩種變量：決策變量和狀態(tài)變量。決策變量包括抽/注水井的數量、位置和流量；狀態(tài)變量就是指地下水的水位和污染物的濃度。水位和濃度狀態(tài)變量對應于決策變量的響應分別由水流模型MODFLOW和溶質運移模型MT3DMS來獲得更新。因此，對于嵌入式模擬-優(yōu)化模型來說，采用進化算法求解時，都要針對決策變量(井流量、井位等)不斷重復調用水流模型MODFLOW與溶質運移模型MT3DMS來計算相應的狀態(tài)變量(水位、濃度)，由此進一步計算目標函數值并判斷約束條件是否滿足，從而搜索得到模型的最優(yōu)策略。

圖2 研究區(qū)場地位置及模型第1層的初始TCE污染羽分布

3 實例應用

3.1 污染場地概述與模擬-優(yōu)化模型

3.1.1 場地概述及模擬模型

研究場址位于美國猶他州大鹽湖以南約10 km的圖埃爾(Tooele)山谷。含水層主要由沖積層構成，但在其中部有一向上抬升的基巖塊，使地下水從沖積沉積物中流向斷裂風化的基巖，然后又回到沖積沉積物中。區(qū)域地下水流大致向東南往西北方向流動。場地主要污染來自于其東南角的一個工業(yè)區(qū)，主要污染物為TCE。監(jiān)測結果表明，污染物主要分布在含水層的淺部，但由于場地地下水流的復雜性，淺部污染羽和深部污染羽的范圍并不一致，甚至東北角的污染羽已延伸到特征邊界之外(見圖2)。

現已確定采用基于水力技術的抽出-處理技術來治理和修復該含水層，該PAT系統(tǒng)由美國某陸軍工程兵團設計，包含16口抽水井、13口注水井和1口附加抽水井(E-fixed)。相應的模擬模型已實現了抽水/不抽水條件下的多次校正，該模型是一個三維、穩(wěn)態(tài)的MODFLOW模型，有4層，165行，和99列。平面上剖分為正方形單元網格，其邊長為60.96 m，垂向上自上而下4個模型層的平均厚度分別為45.72 m、30.48 m、45.72和91.44 m?；贛T3DMS構建與之對應的污染物運移模型。模型時間離散化為有7個3年的應力期來匹配優(yōu)化方案中定義的管理期。

現以校正后的模擬模型計算出的2003年1月的水頭分布與污染羽分布分別作為初始水頭、初始濃度分布。由于現存PAT系統(tǒng)模擬結果顯示，沒有可行解能滿足污染羽截斷邊界的約束條件。經過前期的專家論證，需新增4口注水井(IN1～IN4)，原有PAT系統(tǒng)中的16口抽水井、13口注水井和1口附加抽水井(E-fixed)，僅需保留抽水井E11和E15、注水井I4及E-fixed，但E11和E15均以最大抽水能力運行，各井位置見圖2。需要說明的是，為了維持PAT系統(tǒng)的穩(wěn)定性，所有抽水井與注水井的水量總量相等。

3.1.2 優(yōu)化模型及其求解

本文在確定了井的位置基礎上，令每一個管理期的所有注/抽水井(E-fixed與I4除外)的流量為一個決策變量，以制定出一種最優(yōu)的動態(tài)解。選擇使用GA進行最優(yōu)解的搜索，相關的參數設置如下：POPSIZE=200，PCROSS=0.5， PMUTATE=0.05，NPOSIBL=32。

根據目標函數的一般形式，由于井位、井數已定，總成本項主要來自抽水用電成本項和化學處理項，其中本文把抽水用電成本項當作安裝成本，調整后的通用目標函數可簡化為如下PAT系統(tǒng)的最小運行成本：

(8)

相應的約束條件如下：

|Qsumin|≤|Q*|≤|Qsumax|

(9)

(10)

QI4=-(QEi+QNIi)-QE-fixed

(11)

Cm≤Cmax

(12)

0≤QEi≤│Li│

(13)

0≤QNIi≤Li

(14)

式(8)中，α1取值為34.5K$(1K$表示1 000美元)，α3取值為1.01e-05K$/(m3)；M表示管理期，M=7；N=6；yji是二進制數，當第j管理期新建或優(yōu)化第i個井時，yji=1，否則yji=0。式(9)、(10)都作為全局約束條件：其中式(9)是對總抽水量的約束；Q*是參與模擬的總抽水量(包括E-fixed的抽水量)，Qsumax為最大總抽水量，取值為-41 427.42 m3/d，Qsumin為最小總抽水量，取值為0；式(10)是對抽/注水量和的約束，Qi為參與模擬的第i個井的抽/注水量，i=1,2,…,8；Qc為最大均衡誤差值，取值為5.45 m3/d。NW為各管理期參與模擬的井數，本次優(yōu)化各管理期NW均取值為8。式(11)作為均衡約束條件：QI4是現用注水井I4的注水量；QEi是第i個現用抽水井的抽水量，i=1,2；QNIi是第i個新建注水井的注水量，i=1,2,3,4；QE-fixed是E-fixed的抽水量，值為-7 767.68 m3/d。式(12)為濃度約束條件：Cm是第m個濃度約束邊界單元上TCE的濃度值，Cmax為不同管理期末的約束濃度上限，取值為5 μg/L。式(13)、(14)分別為現用井的抽/注水量、新建注水井的注水量的約束條件：QEi為第i個現存井的抽/注水量，i=1, 2, 3，Li為第i個井的最大抽/注水能力。PAT系統(tǒng)中參與優(yōu)化的井的最大抽/注水能力如表1所示。

圖3 執(zhí)行最優(yōu)策略最終管理期末得到模型第1層中的污染羽分布

圖4 執(zhí)行最優(yōu)策略最終管理期末得到模型第2層中的污染羽分布表1 優(yōu)化前PAT系統(tǒng)中待優(yōu)化抽/注水井的抽/注水能力及其模擬層

井Li/m3/d模擬層E11-3 365.992E15-3 314.212, 3I44 163.472IN13 107.071IN23 107.072IN33 107.071IN43 107.071

如表1所示，抽水井E15抽采2、3兩層，因此只要以其第2層的抽水量作為決策變量，而將第3層的抽水量處理為依賴變量，抽水量的分配依據兩層的導水系數可由優(yōu)化程序自動計算得到。

3.2 優(yōu)化結果與分析

最優(yōu)抽水策略滿足所有約束條件，包括在4個模擬層的污染約束邊界單元處，污染羽濃度在每個管理期末均接近或小于5 μg/L。第1、2含水層是污染物主要分布層，這兩層在第7個管理期的污染羽分布見圖3和圖4。動態(tài)策略下各管理期的抽/注水量見圖5、表2。原穩(wěn)態(tài)策略下的抽/注水量同見圖5。

圖5 最優(yōu)策略下各管理期的抽水(負值)/注水(正值)量分布表2 動態(tài)策略下各管理期內優(yōu)化得到PAT系統(tǒng)的抽水(負值)和注水(正值)方案

井號優(yōu)化前優(yōu)化后管理期抽/注水量/m3/d1234567E11-3 365.99-3 365.99-3 245.85-3 365.99-3 365.99-3 365.99-3 257.14-3 149.81E15-3 314.210-2 876.33-3 314.21-3 314.21-3 314.21-3 314.21-2 886.52NI13 107.071 202.871 109.663 107.071 904.3072 004.552 606.012 706.20NI23 107.072 906.852 407.983 107.073 107.073 107.073 107.072 906.85NI33 107.073 005.743 107.073 107.073 107.073 107.073 107.072 405.42NI43 107.073 005.743 107.073 107.073 005.743 005.742 606.013 107.07I42 019.601 010.784 158.302 019.543 322.113 223.232 912.522 677.75

優(yōu)化后的動態(tài)策略抽水井E15在第1管理期內抽水量為零，并且在第2、第7管理期內低于最大負荷抽水率運行，在整個項目期內抽水井E15抽水量相比優(yōu)化前減少了18.02%，抽水井E11在第2、6、7管理期也是低于最大抽水量運行，項目期內抽水井E11的抽水量減少了1.89%，管理期的總抽水量減小了9.89%，此外，抽注水量之和小于給定最大值5.45 m3/d，說明來自三口抽水井的抽水量都從5口注水井注入含水層中。因此動態(tài)策略下整個項目期的總抽水量減小，使得抽水用電成本和化學處理成本相應減少。

在穩(wěn)定策略下，總成本為1 258.85K$，其 中抽水用電成本、化學處理成本各占73.79%、26.21%。而優(yōu)化后動態(tài)策略下抽水用電成本相比優(yōu)化前減少10.62%；化學處理成本相比優(yōu)化前減少13.08%?？偝杀緸? 117.02K$，相比優(yōu)化前總成本減少了11.27%。

4 結語

本文通過MGO軟件對美國猶他州某處TCE污染場地的PAT修復系統(tǒng)進行模擬優(yōu)化管理，實現了對修復系統(tǒng)的動態(tài)優(yōu)化管理。優(yōu)化結果表明優(yōu)化后的修復策略減少了11.27%的修復成本，同時保證了修復效率。實例研究的結果表明MGO可以有效地為地下水修復提高經濟、高效的管理措施，為今后的地下水污染修復提供一種可靠的工具，具有較好的應用前景。

[1]呂書君. 我國地下水污染分析[J]. 地下水.2009.31(1): 1-5.

[2]Soderstrom M, Boldemann C, Sahlin U, et al. The quality of the outdoor environment influences childrens health-across-sectional study of preschools. Acta Paediatrica, 2013.102: 83-91.

[3]Ozdemir S, Johnson FR. Estimating willingness to pay: do health and environmental researchers have different methodological standards Applied Economics, 2013.45: 2215-2229.

[4]李瑋,王明玉,韓占濤,等. 棕地地下水污染修復技術篩選方法研究—以某廢棄化工廠污染場地為例[J]. 水文地質工程地質.2016.43(3): 131-141.

[5]Minsker BS editor. Long-Term Groundwater Monitoring: The State of the Art[M]. ASCE/EWRI, 0-7844-0678-2, Reston, VA; 2003.

[6]肖傳寧,盧文喜,安永凱,等. 基于兩種耦合方法的模擬-優(yōu)化模型在地下水污染源識別中的對比[J]. 中國環(huán)境科學.2015.35(8): 2393-2399.

[7]吳劍鋒,彭偉, 錢家忠, 等. 基于INPGA的地下水污染治理多目標優(yōu)化管理模型:I—理論方法與算例驗證[J]. 地質評論.2011.57(2): 277-283.

[8]吳劍鋒,彭偉, 錢家忠, 等. 基于INPGA的地下水污染治理多目標優(yōu)化管理模型: II.實例應用[J]. 地質評論.2011.57(3): 1-7.

[9]Zhang JL, Li YP, Huang GH. A robust simulation-optimization modeling system for effluent trading—a case study of nonpoint source pollution control[J]. Environmental Science and Pollution Research.2014.21: 5036-5053.

[10]Zheng C, Wang PP. MGO-A Modular Groundwater Optimizer Incorporating MODFLOW/MT3DMS, Documentation and User’s Guide[R]. The University of Alabama and Groundwater Research Ltd., Tuscaloosa, AL, 2003.

[11]McDonald MG, Harbaugh AW. A modular three-dimensional finite-difference groundwater flow model. Techniques of Water Resources Investigations of the U.S. Geological Survey, Book 6[R], U.S. Government Printing Office, Washington, D.C., 1988.

[12]Zheng C, Wang PP. MT3DMS: A Modular Three-Dimensional Multispecies Transport Model for Simulation of Advection, Dispersion and Chemical Reactions of Contaminants in Groundwater Systems; Documentation and User’s Guide, Contract Report SERDP-99-1[R]. U.S. Army Engineer Research and Development Center, Vicksburg, MS. 1999.