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

　　目前，國際上地下水硝酸鹽污染治理技術(shù)分為原位處理技術(shù)和抽出處理技術(shù). 抽出處理技術(shù)仍然是最常用的一種. 抽出處理技術(shù)的核心在于如何有效地把污染的地下水從含水層中抽出來，通常考慮的因素包括抽水井的位置、 抽水井?dāng)?shù)、 抽水速率、 安裝井和抽水處理成本，需要找到一個(gè)合適的方法得到最優(yōu)解，從而使治理成本最小化[4]. 因此，其關(guān)鍵是多目標(biāo)優(yōu)化問題. 目前，國外已有許多模擬優(yōu)化技術(shù)運(yùn)用到抽出處理技術(shù)設(shè)計(jì)[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è)計(jì)，在其設(shè)計(jì)中包含井的安裝費(fèi)用和抽水費(fèi)用. McKinney等[18]運(yùn)用遺傳算法解決了最大抽水量和最低抽水成本等地下水優(yōu)化問題. Yoon等[19]運(yùn)用實(shí)數(shù)編碼遺傳算法來確定最優(yōu)抽水速率和最佳井的位置，最大限度減小了地下水修復(fù)系統(tǒng)成本. Erickson等[20]利用小生境技術(shù)遺傳算法解決了抽出處理地下水修復(fù)問題，實(shí)現(xiàn)了修復(fù)成本和剩余污染物濃度最小化. 吳劍鋒等[21]開發(fā)改進(jìn)了一種小生境Pareto遺傳算法(INPGA)，以美國麻省軍事保護(hù)區(qū)為實(shí)例，通過建立研究區(qū)復(fù)雜地下水污染治理的多目標(biāo)優(yōu)化管理模型，并添加消息傳遞接口(MPI)的并行計(jì)算和個(gè)體適應(yīng)度，提高了運(yùn)行庫計(jì)算速度.

　　本文以北京某場地生活垃圾填埋場硝酸鹽污染地下水為研究對象，采用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è)計(jì)

　　1.1 遺傳算法

　　遺傳算法(Genetic Algorithm)是模擬達(dá)爾文生物進(jìn)化論的自然選擇和遺傳學(xué)機(jī)理的生物進(jìn)化過程的計(jì)算模型，是一種通過模擬自然進(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)境的個(gè)體，檢查其結(jié)果是否符合優(yōu)化的要求，如果不符合，重復(fù)其上述過程，直到優(yōu)化出滿意的結(jié)果.

　　1.2 模擬退火法

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

　　步驟1：給定初始溫度T0解X0算f0=f(x0); 步驟2：隨機(jī)產(chǎn)生擾動Δx,計(jì)算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表示單位時(shí)間從單位體積含水層流入或流出的水量(d-1); Ss表示貯水率(m-1); t表示時(shí)間(d). MT3DMS地下水溶質(zhì)運(yùn)移方程：

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

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

　　最小化：

　　約束條件：

　　式中,J表示目標(biāo)函數(shù)，a1表示安裝一口井的總費(fèi)用，W表示研究區(qū)備選井的數(shù)量，N是總的優(yōu)選井?dāng)?shù)量，即從備選井中優(yōu)化出來井的數(shù)量，a2是抽出和處理單位體積的水的費(fèi)用[4]，Qi單位時(shí)間抽水的體積(m3 ·d-1)，Qmin和Qmax分別是最小抽水和最大抽水量，Δti是井i持續(xù)抽水時(shí)間(d)，hi是井i的水頭(m)，hmin和hmax分別是最小水頭和最大水頭(m); ci是井i的溶質(zhì)的濃度(mg ·L-1)，c*表示目標(biāo)治理濃度(mg ·L-1). 上述定義多目標(biāo)優(yōu)化方案也可以通過試錯(cuò)法解決. 雖然試錯(cuò)法簡單、 易用，但測試和檢驗(yàn)煩瑣，且不能保證組合的最優(yōu)方案，容易產(chǎn)生局部優(yōu)化，往往很難找到最優(yōu)的抽出處理方案. 所以在滿足所有約束條件下，應(yīng)用遺傳算法和模擬退火法比試錯(cuò)法更有效解決井的位置及組合問題[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)行單孔抽水試驗(yà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，同時(shí)地表 50 cm滲坑NO-3-N高達(dá)65 mg ·L-1. 研究區(qū)地下水硝酸鹽仍處于不斷變化的情況，同時(shí)場地調(diào)查尚不夠充分，不足以刻畫出實(shí)際地下水硝酸鹽的污染羽狀況. 因此，本研究以農(nóng)村垃圾和畜禽養(yǎng)殖廢水為污染源，硝酸鹽入滲濃度為65 mg ·L-1，模擬在天然地下水流動情況下釋放5 a形成的硝酸鹽污染羽作為抽出處理的初始污染羽. 研究區(qū)硝酸鹽初始污染羽總量為36.99 kg，見圖 3.

　　2.3 抽出處理備選井方案

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

　　研究區(qū)周圍地下水硝酸鹽背景濃度約為1 mg ·L-1. 將研究區(qū)紅色框作為濃度約束區(qū)，在100 d抽出處理修復(fù)時(shí)間內(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. 西北方向和東南方向?yàn)楦羲吔? 地下水硝酸鹽的阻滯因子設(shè)為1. 地下水流模擬設(shè)置為穩(wěn)定流模擬(兩種優(yōu)化方法只適合穩(wěn)定流[4])，模型輸入?yún)?shù)見表 2.

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

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

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

　　表 3 費(fèi)用系數(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后的地下水硝酸鹽(以氮計(jì))達(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è)コ糠謩e為28.45 kg和31.42 kg，其總量去除率分別76.89%和84.92%，見圖 4. 優(yōu)化過程中發(fā)現(xiàn)，總抽水速率和硝酸鹽總?cè)コ看嬖诮凭€性增加的趨勢.

　　GA、 SA同時(shí)優(yōu)化了井的位置和抽水量. 圖 4可以看出研究區(qū)中18口備選井方案，兩種優(yōu)化方法最優(yōu)的兩口井都是8號井和14號井. 說明在該場地抽出處理中，只需要2口井就可以滿足約束條件，到達(dá)治理目標(biāo)，比初始方案減少16口井，節(jié)省了井的安裝成本. 同時(shí)說明最優(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具體參見 污水處理技術(shù)資料或污水技術(shù)資料更多相關(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，同時(shí)其硝酸鹽總量去除率分別為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è)計(jì)合理的抽出處理方案和符合管理目標(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)境工程北京市重點(diǎn)實(shí)驗(yàn)室 姜烈、何江濤、姜永海、劉菲)