地下水是重要的飲用水源. 由于人類發(fā)展忽略了環(huán)境保護(hù)，導(dǎo)致地下水污染嚴(yán)重，其中硝酸鹽污染問題突出[1]. 過量的硝酸鹽會對人類健康造成嚴(yán)重危害，例如藍(lán)嬰病[2]. 因此，治理地下水硝酸鹽刻不容緩[3].

　　目前，國際上地下水硝酸鹽污染治理技術(shù)分為原位處理技術(shù)和抽出處理技術(shù). 抽出處理技術(shù)仍然是最常用的一種. 抽出處理技術(shù)的核心在于如何有效地把污染的地下水從含水層中抽出來，通常考慮的因素包括抽水井的位置、 抽水井?dāng)?shù)、 抽水速率、 安裝井和抽水處理成本，需要找到一個合適的方法得到最優(yōu)解，從而使治理成本最小化[4]. 因此，其關(guān)鍵是多目標(biāo)優(yōu)化問題. 目前，國外已有許多模擬優(yōu)化技術(shù)運(yùn)用到抽出處理技術(shù)設(shè)計[5]. 傳統(tǒng)的優(yōu)化方法包括線性規(guī)劃、 非線性規(guī)劃、 動態(tài)規(guī)劃、 二次規(guī)劃和整數(shù)規(guī)劃[6]. 較新的優(yōu)化方法包括遺傳算法(GA)[7, 8, 9, 10, 11]和模擬退火法(SA)[12, 13, 14, 15, 16]. 這些優(yōu)化方法可以耦合到地下水模擬模型中.

　　抽出處理在優(yōu)化過程中容易陷入局部優(yōu)化. GA和SA采用了多種全局優(yōu)化方法，能夠?qū)崿F(xiàn)全局優(yōu)化. Dougherty等[17]首次將模擬退火法(SA)運(yùn)用到地下水修復(fù)最優(yōu)設(shè)計，在其設(shè)計中包含井的安裝費用和抽水費用. McKinney等[18]運(yùn)用遺傳算法解決了最大抽水量和最低抽水成本等地下水優(yōu)化問題. Yoon等[19]運(yùn)用實數(shù)編碼遺傳算法來確定最優(yōu)抽水速率和最佳井的位置，最大限度減小了地下水修復(fù)系統(tǒng)成本. Erickson等[20]利用小生境技術(shù)遺傳算法解決了抽出處理地下水修復(fù)問題，實現(xiàn)了修復(fù)成本和剩余污染物濃度最小化. 吳劍鋒等[21]開發(fā)改進(jìn)了一種小生境Pareto遺傳算法(INPGA)，以美國麻省軍事保護(hù)區(qū)為實例，通過建立研究區(qū)復(fù)雜地下水污染治理的多目標(biāo)優(yōu)化管理模型，并添加消息傳遞接口(MPI)的并行計算和個體適應(yīng)度，提高了運(yùn)行庫計算速度.

　　本文以北京某場地生活垃圾填埋場硝酸鹽污染地下水為研究對象，采用Visual MODFLOW 4.6中Modular Groundwater Optimizer (MGO)[22]，將水流模擬軟件 MODFLOW、 溶質(zhì)運(yùn)移MT3DMS與遺傳算法和模擬退火法相結(jié)合，嘗試解決地下水硝酸鹽污染抽出處理優(yōu)化模擬，優(yōu)化井的位置、 井的抽水速率及最小化治理成本. 其執(zhí)行程序通過FORTRAN自動編碼器來匯編，屬于內(nèi)置模塊，在優(yōu)化過程中采用DOS界面運(yùn)行.

　　1 地下水硝酸鹽抽出處理優(yōu)化系統(tǒng)設(shè)計

　　1.1 遺傳算法

　　遺傳算法(Genetic Algorithm)是模擬達(dá)爾文生物進(jìn)化論的自然選擇和遺傳學(xué)機(jī)理的生物進(jìn)化過程的計算模型，是一種通過模擬自然進(jìn)化過程搜索最優(yōu)解的方法，它最初由美國Michigan大學(xué)Holland教授于1975年首先提出. 采用概率化的尋優(yōu)方法，能自動獲取和指導(dǎo)優(yōu)化的搜索空間，自適應(yīng)地調(diào)整搜索方向，不需要確定的規(guī)則. 遺傳算法的核心是選擇、 交叉、 變異. 更多遺傳算法的描述參見文獻(xiàn)[23, 24]. 從任一初始種群出發(fā)，通過選擇、 交叉、 變異操作，產(chǎn)生一群更適應(yīng)環(huán)境的個體，檢查其結(jié)果是否符合優(yōu)化的要求，如果不符合，重復(fù)其上述過程，直到優(yōu)化出滿意的結(jié)果.

　　1.2 模擬退火法

　　模擬退火法，源于統(tǒng)計熱力物理學(xué)，它模擬熔融狀態(tài)下物體逐漸冷卻達(dá)到結(jié)晶狀態(tài)的物理過程. 材料中的原子原來會停留在使內(nèi)能有局部最小值的位置，加熱使能量變大，原子會離開原來位置，而隨機(jī)在其他位置中移動. 退火冷卻時速度較慢，使得原子有較多可能可以找到內(nèi)能比原先更低的位置. 不依賴初始模型的選擇，能尋找全局最小點而不陷入局部極小[19]. 模擬退火算法的操作步驟如下.

　　步驟1：給定初始溫度T0解X0算f0=f(x0); 步驟2：隨機(jī)產(chǎn)生擾動Δx,計算f1=f(x0+Δx)和Δf=f1-f0; 步驟3：若Δf≤0則轉(zhuǎn)步驟5，否則產(chǎn)生隨機(jī)數(shù)R∈[0,1]; 步驟4：若p=exp(-Δf/T)≥R，則用x′=x0+Δx代替x0， f1代替f0; 步驟5：檢查M過程是否結(jié)束，若未結(jié)束則轉(zhuǎn)步驟2; 步驟6：降低溫度參數(shù)T; 步驟7：判斷終止準(zhǔn)則是否滿足，若滿足，則算法終止，否則，轉(zhuǎn)步驟2.

　　1.3 地下水多目標(biāo)優(yōu)化模型

　　地下水優(yōu)化模型通常包含兩組變量：決策變量和狀態(tài)變量. 在優(yōu)化過程中的目的是識別這些決策變量的最佳組合. 狀態(tài)變量依賴地下水流方程的水頭和地下水溶質(zhì)運(yùn)移方程的濃度. MODFLOW水流方程如下：

　　式中，Kx、 Ky、 Kz分別為滲透系數(shù)在x、 x、 y方向的分量; h表示水頭(m); qs表示單位時間從單位體積含水層流入或流出的水量(d-1); Ss表示貯水率(m-1); t表示時間(d). MT3DMS地下水溶質(zhì)運(yùn)移方程：

　　式中,θ表示含水層的孔隙度(無量綱); cK表示溶質(zhì)組分K的濃度(mg ·L-1); D ij是水動力彌散系數(shù)張量(m2 ·d-1); vi是孔隙中實際水流速度(m ·d-1); qs表示單位時間從單位體積含水層流入或流出的水量(d-1); cKs是源匯項溶質(zhì)組分K的濃度(mg ·L-1); �И禦n�П硎净瘜W(xué)反應(yīng)項總和[mg ·(L ·d)-1].

　　優(yōu)化設(shè)計過程中，目標(biāo)函數(shù)包含井的安裝費用、 抽水費用、 抽水地下水硝酸鹽治理費用等. 其數(shù)學(xué)表達(dá)式如下.

　　最小化：

　　約束條件：

　　式中,J表示目標(biāo)函數(shù)，a1表示安裝一口井的總費用，W表示研究區(qū)備選井的數(shù)量，N是總的優(yōu)選井?dāng)?shù)量，即從備選井中優(yōu)化出來井的數(shù)量，a2是抽出和處理單位體積的水的費用[4]，Qi單位時間抽水的體積(m3 ·d-1)，Qmin和Qmax分別是最小抽水和最大抽水量，Δti是井i持續(xù)抽水時間(d)，hi是井i的水頭(m)，hmin和hmax分別是最小水頭和最大水頭(m); ci是井i的溶質(zhì)的濃度(mg ·L-1)，c*表示目標(biāo)治理濃度(mg ·L-1). 上述定義多目標(biāo)優(yōu)化方案也可以通過試錯法解決. 雖然試錯法簡單、 易用，但測試和檢驗煩瑣，且不能保證組合的最優(yōu)方案，容易產(chǎn)生局部優(yōu)化，往往很難找到最優(yōu)的抽出處理方案. 所以在滿足所有約束條件下，應(yīng)用遺傳算法和模擬退火法比試錯法更有效解決井的位置及組合問題[22].

　　圖 2 研究區(qū)三維地質(zhì)剖面

　　2 材料與方法

　　2.1 研究區(qū)水文地質(zhì)條件

　　研究區(qū)位于北京某場地，所處區(qū)域地勢平坦. 場地東側(cè)為養(yǎng)殖場，西南側(cè)為當(dāng)?shù)剞r(nóng)民耕作農(nóng)田，農(nóng)村垃圾和畜禽養(yǎng)殖廢水均排放到池塘(圖 1).

　　圖 1 研究區(qū)初始流場示意

　　依據(jù)場地調(diào)查資料，利用CTECH對研究區(qū)構(gòu)建三維地質(zhì)剖面圖，第一層主要含水層為非承壓含水層，地下水標(biāo)高16~23 m，厚度約6.92 m，以粉細(xì)砂為主，局部夾有砂質(zhì)粉土、 黏質(zhì)粉土和粉質(zhì)黏土(圖 2). 研究區(qū)第一含水層下面主要為砂質(zhì)粉土、 黏質(zhì)粉土，滲透性很差，能夠起到較好的隔水作用.

　　在研究區(qū)進(jìn)行單孔抽水試驗. 通過Aquifer Test軟件分析，結(jié)果見表 1. 研究區(qū)含水層的平均滲透系數(shù)為3.37 m ·d-1.

　　表 1 不同求參方法中滲透系數(shù)和導(dǎo)水系數(shù)結(jié)果表

　　2.2 研究區(qū)地下水污染狀況

　　利用Geoprobe鉆機(jī)對生活垃圾填埋場地下水及池塘污水進(jìn)行取樣調(diào)查. 現(xiàn)場測試數(shù)據(jù)：池塘地表水的NH+4-N高達(dá)46.8 mg ·L-1，同時地表 50 cm滲坑NO-3-N高達(dá)65 mg ·L-1. 研究區(qū)地下水硝酸鹽仍處于不斷變化的情況，同時場地調(diào)查尚不夠充分，不足以刻畫出實際地下水硝酸鹽的污染羽狀況. 因此，本研究以農(nóng)村垃圾和畜禽養(yǎng)殖廢水為污染源，硝酸鹽入滲濃度為65 mg ·L-1，模擬在天然地下水流動情況下釋放5 a形成的硝酸鹽污染羽作為抽出處理的初始污染羽. 研究區(qū)硝酸鹽初始污染羽總量為36.99 kg，見圖 3.

　　2.3 抽出處理備選井方案

　　為了設(shè)計合理的抽出處理治理優(yōu)化方案，根據(jù)模擬生成硝酸鹽污染羽形成的范圍及通過試錯法確定的較為有效的抽水井位置，在研究區(qū)布置了18口備選井方案(圖 3). 目的是在滿足治理成本最小情況下，抽出處理后達(dá)到目標(biāo)治理濃度，同時確定最佳井?dāng)?shù)及其井的位置和抽水速率.

　　研究區(qū)周圍地下水硝酸鹽背景濃度約為1 mg ·L-1. 將研究區(qū)紅色框作為濃度約束區(qū)，在100 d抽出處理修復(fù)時間內(nèi)，其硝酸鹽目標(biāo)修復(fù)濃度低于10 mg ·L-1.

　　2.4 模型概化及參數(shù)確定

　　研究區(qū)場地范圍較小，可以將研究區(qū)假設(shè)為均質(zhì)各向異性含水層. 地下水硝酸鹽污染濃度主要分布在第一含水層，該層下部主要為砂質(zhì)粉土、 黏質(zhì)粉土，滲透系數(shù)較差，能夠起到較好的隔水作用，故 而將含水層概化成一層. 研究場地為225 m×200m，將其剖分成2.25 m×2 m. 地下水從北東方向向南西方向流(圖 1). 研究區(qū)的初始水力坡度約為0.01. 西北方向和東南方向為隔水邊界. 地下水硝酸鹽的阻滯因子設(shè)為1. 地下水流模擬設(shè)置為穩(wěn)定流模擬(兩種優(yōu)化方法只適合穩(wěn)定流[4])，模型輸入?yún)��?shù)見表 2.

　　圖 3 地下水硝酸鹽初始污染羽分布示意和備選井位置

　　表 2 模型輸入?yún)��?shù)

　　表 3列出了安裝一口井的總費用、 抽出和處理單位體積水的費用[4]. 最小與最大水頭值、 最小與最大抽水量及硝酸鹽目標(biāo)治理濃度等約束條件.

　　表 3 費用系數(shù)和約束變量

　　3 優(yōu)化結(jié)果分析與討論

　　3.1 抽出處理優(yōu)化

　　本研究中GA設(shè)定的相關(guān)參數(shù)如下：種群數(shù)目為100，最大遺傳代數(shù)為100，每口井的抽水速率離散化區(qū)間數(shù)為32，交叉概率為0.90，變異概率為0.05. SA設(shè)定的相關(guān)參數(shù)如下：初始溫度為1，抽樣次數(shù)10000，降溫率為0.9，最大迭代次數(shù)為100.

　　兩種優(yōu)化方法都要求濃度約束區(qū)治理100 d后的地下水硝酸鹽(以氮計)達(dá)到10 mg ·L-1. GA優(yōu)化結(jié)果：8號井優(yōu)化抽水速率為155 m3 ·d-1，14號井優(yōu)化抽水速率為10 m3 ·d-1; SA優(yōu)化結(jié)果：8號井優(yōu)化抽水速率為82 m3 ·d-1，14號井優(yōu)化抽水速率為39 m3 ·d-1. 其他井都為零(表 4). GA在8號井優(yōu)化的抽水速率比SA在8號井的要大，在14號井優(yōu)化的抽水速率比SA在14號優(yōu)化的速率要小. 但SA總抽水速率(121 m3 ·d-1)比GA(165 m3 ·d-1)小.

　　表 4 優(yōu)化井位置與抽水速率

　　基于GA和SA優(yōu)化的抽水速率，100 d后硝酸鹽總?cè)コ�量分別為28.45 kg和31.42 kg，其總量去除率分別76.89%和84.92%，見圖 4. 優(yōu)化過程中發(fā)現(xiàn)，總抽水速率和硝酸鹽總?cè)コ�量存在近似線性增加的趨勢.

　　GA、 SA同時優(yōu)化了井的位置和抽水量. 圖 4可以看出研究區(qū)中18口備選井方案，兩種優(yōu)化方法最優(yōu)的兩口井都是8號井和14號井. 說明在該場地抽出處理中，只需要2口井就可以滿足約束條件，到達(dá)治理目標(biāo)，比初始方案減少16口井，節(jié)省了井的安裝成本. 同時說明最優(yōu)井的位置在污染羽的中軸線上中游和中下游抽水布置最為合理.

　　圖 4 GA和SA抽水治理100 d后地下水硝酸鹽分布示意

　　3.2 GA與SA對比分析

　　GA、 SA不僅可以優(yōu)化井的位置、 抽水速率及井?dāng)?shù)，而且還可以優(yōu)化系統(tǒng)治理成本(井安裝成本與抽出處理成本之和)，見表 5. GA和SA優(yōu)化井?dāng)?shù)一樣，其安裝成本也就一樣，但SA優(yōu)化總的抽水量比GA要低，所以SA優(yōu)化的系統(tǒng)成本要比GA要低，說明SA在多目標(biāo)地下水硝酸鹽修復(fù)優(yōu)化管理要比GA體現(xiàn)出一定的優(yōu)越性，體現(xiàn)在節(jié)省系統(tǒng)治理成本.

　　圖 5為GA與SA遺傳代數(shù)與目標(biāo)函數(shù)優(yōu)化收斂對比. GA在40代后優(yōu)化目標(biāo)函數(shù)趨于收斂，而SA在20代后趨于穩(wěn)定. GA在前20代目標(biāo)函數(shù)優(yōu)化波動性比SA要大. 兩者最終都趨于穩(wěn)定. 兩種優(yōu)化方法識別決策變量的最佳組合. 在最小化治理成本下，全局優(yōu)化井?dāng)?shù)、 井的抽水量、 安裝成本、 抽水處理成本，使得優(yōu)化結(jié)果可信. SA的系統(tǒng)治理成本要比GA少3338.45元，節(jié)省成本約6.8%，見表 5具體參見污水寶商城資料或http://www.northcarolinalenders.com更多相關(guān)技術(shù)文檔。

　　圖 5 GA、 SA目標(biāo)函數(shù)優(yōu)化收斂對比

　　表 5 GA、 SA優(yōu)化系統(tǒng)治理成本/元

　　4 結(jié)論

　　(1)本研究依據(jù)北京某場地調(diào)查資料，利用遺傳算法(GA)和模擬退火法(SA)對地下水硝酸鹽污染羽治理區(qū)中布置18口備選井方案進(jìn)行優(yōu)化. 為達(dá)到100 d將研究區(qū)硝酸鹽濃度降低到10 mg ·L-1的目標(biāo)，GA和SA優(yōu)化的井都為8號井和14號井，優(yōu)化總抽水速率分別為165 m3 ·d-1和121 m3 ·d-1，同時其硝酸鹽總量去除率分別為76.89%和84.92%.

　　(2)優(yōu)化結(jié)果說明，最佳井的位置應(yīng)位于硝酸鹽污染羽中游及中下游的中軸線上，且中游抽水速率比下游要大. 兩種優(yōu)化算法對比表明模擬退火法(SA)優(yōu)化系統(tǒng)治理成本比遺傳算法(GA)節(jié)省6.8%，且波動性小，更收斂.

　　(3)模擬退火法(SA)和遺傳算法(GA)耦合到地下水流模擬程序MODFLOW和溶質(zhì)運(yùn)移程序MT3DMS中，可以設(shè)計合理的抽出處理方案和符合管理目標(biāo)修復(fù)系統(tǒng)，為地下水污染現(xiàn)場修復(fù)提供科學(xué)的指導(dǎo).

　　致謝： 本研究場地基礎(chǔ)信息調(diào)查由中國環(huán)境科學(xué)研究院、 北京建工環(huán)境修復(fù)有限責(zé)任公司、 輕工環(huán)保研究所等共同完成. 此外，南京大學(xué)吳劍鋒老師對遺傳算法優(yōu)化方法工作給予指導(dǎo)!在此一并表示感謝.(來源及作者：中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院, 水資源與環(huán)境工程北京市重點實驗室 姜烈、何江濤、姜永海、劉菲)