吳睿奇，朱國榮，王 佩

(南京大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 210093)

雙種群進化粒子群算法求解地下水管理模型

吳睿奇，朱國榮，王 佩

(南京大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 210093)

為避免粒子群算法(PSO)早熟的缺點，設(shè)計了一種雙種群進化粒子群算法(DE-PSO)。DE-PSO是基于PSO，引入選擇、交叉及差分變異操作，并結(jié)合合理有效的粒子評價方法及越界處理方法之后形成的。將DE-PSO應(yīng)用于兩個地下水管理模型算例，第一個算例DE-PSO解的總抽水量分別比遺傳算法(GA)、模擬退火算法(SA)和PSO減少了64、256、207 m3/d，第二個算例DE-PSO解的總治理成本分別比GA、SA和PSO減少了57.74、151.93、76.59萬元。兩個算例中DE-PSO都表現(xiàn)出穩(wěn)定的進化趨勢，尋優(yōu)效率好于GA、SA和PSO，可以有效求解地下水管理模型問題。

地下水管理模型;粒子群算法;雙種群;差分變異

0 引言

地下水管理模型是實現(xiàn)地下水科學(xué)管理的有力工具(于福榮等，2010)，其實質(zhì)是尋求以最少的社會經(jīng)濟環(huán)境成本(包括安裝井、抽水和治理污染等方面)來滿足既定管理目標(biāo)(Wanger，1995)。目前，求解地下水管理模型的方法包括傳統(tǒng)優(yōu)化算法和智能算法。傳統(tǒng)優(yōu)化方法包括線性規(guī)劃(LP)和非線性規(guī)劃(NLP)，前者只能求解目標(biāo)函數(shù)和約束條件都是關(guān)于決策變量的線性函數(shù)的問題，后者要求目標(biāo)函數(shù)和約束條件連續(xù)、可導(dǎo)，而且容易陷入局部最優(yōu);智能算法包括遺傳算法(GA)、模擬退火算法(SA)、禁忌搜索(TS)等，這些智能算法具有啟發(fā)式尋優(yōu)特點，已被應(yīng)用于地下水優(yōu)化管理中，但由于所有智能算法都基于隨機搜索，無法確保獲得全局最優(yōu)解，因此如何提高全局尋優(yōu)能力仍是這些智能算法的改進方向。

粒子群算法(PSO)由 Kenney和 Eberhart于1995年提出，是基于群智能的隨機全局優(yōu)化技術(shù)，具有控制參數(shù)少、收斂快、原理簡單等優(yōu)點。然而，PSO容易早熟且在處理約束優(yōu)化時存在粒子越界飛行問題。目前，許多學(xué)者進行了傳統(tǒng)PSO的改進研究，提出了下列具體改進方法。①基于PSO公式中控制參數(shù)的改進。② 借鑒遺傳算法，引入選擇、交叉、變異(高斯變異、柯西變異、差分變異等)操作。③設(shè)立不同層次、功能的子種群。④ 將PSO與GA、SA、TS、模式搜索等算法嵌套成混合智能算法。⑤將PSO模型看成一種熱力學(xué)系統(tǒng)，引入熱力學(xué)思想。⑥基于粒子評價方法和越界處理方法的改進(針對約束優(yōu)化)。⑦ 以上方法的綜合。綜合考慮不同改進方法的復(fù)雜程度及有效性，筆者對PSO作了3個方面的改進。①設(shè)立2個平行的子種群進行優(yōu)化，2個子種群通過粒子更新公式及選擇、交叉操作交流信息。②對粒子更新公式作了調(diào)整，舍棄了慣性項，使粒子速度不再有繼承性，修正了粒子學(xué)習(xí)機制，使其不僅可以向自身個體歷史極值和整個群體的歷史極值學(xué)習(xí)，同時也向粒子所在子種群的歷史極值學(xué)習(xí)，并接受鄰群粒子的隨機差分變異擾動。③采用合理的粒子評價方法和越界處理方法，使其更適于求解約束優(yōu)化。經(jīng)過上述改進的算法借鑒了GA中生物進化(選擇、交叉、變異)思想，稱為雙種群進化粒子群算法(DE-PSO)，DE-PSO避免了PSO的早熟缺點，提高了全局尋優(yōu)能力。

1 地下水管理模型

地下水管理模型由水流模型、溶質(zhì)運移模型及優(yōu)化模型耦合而成。前者用來不斷更新狀態(tài)變量，包括地下水流運動模型和溶質(zhì)運移模型;后者用來選擇最優(yōu)決策變量。狀態(tài)變量通常指水頭和溶質(zhì)濃度;決策變量包括抽水量(注水量)、井的數(shù)目、井的位置、治理周期末污染物的總量，還可能包括每口井的“開/關(guān)”狀態(tài)等(楊蘊等，2009)。

1.1 水流模型和溶質(zhì)運移模型

假設(shè)水的密度為常數(shù)，當(dāng)各向異性介質(zhì)的主方向與坐標(biāo)軸的方向一致時，三維地下水流的運動方程可表示為：

不考慮吸附作用及化學(xué)反應(yīng)，描述三維非穩(wěn)定滲流系統(tǒng)中的溶質(zhì)運移的方程可表示為：

式(1)、(2)、(3)中：Ki為滲透系數(shù)張量在 i方向的主值(LT-1);b為含水層飽和帶厚度(L);h為水頭(L);qs為源匯項(LT-1);μ為給水度(無量綱);t為時間(T);c為溶質(zhì)體積質(zhì)量分數(shù)(ML-3);Dij為水動力彌散系數(shù)(L2T-1);cs為源或匯組分中的溶質(zhì)體積質(zhì)量分數(shù)(ML-3);θ為含水介質(zhì)的孔隙度;νi為地下水沿不同方向i的達西流速(LT-1)。

式(1)與相應(yīng)邊界條件及初始條件結(jié)合構(gòu)成水流模型。式(2)為溶質(zhì)運移基本方程，其中的達西流速vi通過式(3)求得，而式(3)中的水頭h通過式(1)求得。因此，溶質(zhì)運移模型是由式(1)—式(3)與相應(yīng)邊界條件及初始條件構(gòu)成。

1.2 優(yōu)化模型

優(yōu)化模型由目標(biāo)函數(shù)及約束條件組成，目標(biāo)函數(shù)通常是總的價格、總的抽水量或總的抽出污染物質(zhì)量等。約束條件通常包含水位約束、水力梯度約束、濃度約束和流量約束等(楊蘊等，2009)。

2 雙種群進化粒子群算法(DE-PSO)求解地下水管理模型

2.1 粒子群算法(PSO)

粒子位置和速度根據(jù)以下公式進行更新：

式(4)、(5)中：vk為第k次迭代的粒子速度，xk為第k次迭代的粒子位置;Pk為第k次迭代粒子歷史最優(yōu)解，P為第k次迭代群體歷史最優(yōu)解;m為慣性權(quán)重，W1、W2為學(xué)習(xí)因子1和學(xué)習(xí)因子2，常取W1=W2=2.0;vmax表示粒子最大速度，K為系數(shù)，一般取0.1～1.0;S為粒子位置最大變化幅度，其值等于上下邊界值之差;r1和r2為0～1均勻分布的隨機數(shù)。

2.2 雙種群進化粒子群算法(DE-PSO)

基于子種群思想，DE-PSO設(shè)置兩個平行的子種群1和2進行優(yōu)化，子種群1和2的群體規(guī)模一致，兩個子種群的粒子同步進化，其粒子更新公式、粒子評價方法及越界處理方法一致。子種群間通過粒子更新公式和選擇、交叉操作交流信息。

2.2.1 引入差分變異的粒子更新公式 PSO中慣性權(quán)重的確定沒有普適的方法，對于復(fù)雜的非線性優(yōu)化問題，即使采用自適應(yīng)方法調(diào)整慣性權(quán)重，搜索精度提高效果也不明顯(羅德相，2009)。許多研究者發(fā)現(xiàn)PSO算法中速度的概念并不是必須的，徐星等(2011)取消了粒子速度的概念，保留了PSO的學(xué)習(xí)機制。筆者設(shè)計的DE-PSO雖保留了粒子速度概念，但取消了慣性權(quán)重;粒子速度無需初始化，且不再具有繼承性。差分變異算子對粒子群算法性能的改進有很大的促進作用。因此，DE-PSO將其引入到粒子更新公式中，改進后的DE-PSO粒子更新公式為：

式(6)、(7)中：α為分配系數(shù)，取值在0～1之間，經(jīng)實驗α取0.5較合適;X、X為從鄰群中隨機選擇的兩個不同粒子;P表示第k次迭代整個群體(包含兩個子種群)的歷史最優(yōu)解，P表示第k次迭代粒子所在子種群的歷史最優(yōu)解;r1、r2、r3、r4為0～1均勻分布的隨機數(shù);其余符號含義同式(4)、式(5)。

2.2.2 選擇與交叉操作 選擇操作可以不斷淘汰較差的粒子，使算法對于當(dāng)前較好區(qū)域有更好的開發(fā)能力;而交叉操作與子種群結(jié)合，可以使粒子受益于父母雙方，增強搜索能力，易于跳出局優(yōu)(汪定偉等，2007)。DE-PSO在每次粒子評價結(jié)束后，分別對兩個子種群的粒子進行排序，選擇兩個子種群中較好的一半粒子作為父代粒子。設(shè)子種群的群體規(guī)模均為N，按公式(8)進行N次交叉操作，所得2N個子代粒子最后放回子種群1和2，取代原來粒子構(gòu)成新種群。

式中：parent1和parent2表示父代粒子，分別隨機選自子種群1和2，且是所在子種群排名前N/2的粒子;child1和child2表示由父代粒子 parent1和 parent2交叉產(chǎn)生的兩個子代粒子，分別被放回子種群1和2;xi表示粒子第i維的值;pi為0～1均勻分布的隨機數(shù)。

2.2.3 粒子評價與越界處理 當(dāng)粒子違反約束條件時，一般可以采用罰函數(shù)法或拒絕策略，罰函數(shù)法會限制算法求解精度，拒絕策略則很難執(zhí)行，因為對于一個復(fù)雜問題，即使重復(fù)多次初始化也很難獲得較多的可行解，完全摒棄或不評價不可行解將導(dǎo)致優(yōu)化難以進行(汪定偉等，2007)。劉衍民等(2011)根據(jù)粒子違背約束條件的程度，提出一種新的粒子間比較準(zhǔn)則，DE-PSO借鑒其思想進行粒子評價。設(shè)有兩個粒子X、Y，則在以下兩種情況下DE-PSO判定X比Y更優(yōu)：①約束條件的滿足個數(shù)相同，但X的目標(biāo)函數(shù)值比Y小(最小值問題);②X約束條件的滿足個數(shù)比Y多。

PSO在處理約束問題時，決策變量越界是經(jīng)常發(fā)生的，目前，處理粒子越界主要有吸收墻、反射墻、不可見墻3種方法，吸收墻和反射墻能使粒子達到邊界值，但同時會加大算法在邊界附近陷入局優(yōu)的風(fēng)險，不可見墻降低了算法在邊界附近的早熟風(fēng)險，但會導(dǎo)致邊界值不可到達，也會導(dǎo)致額外升級和進化不均。DE-PSO對越界粒子作如下處理：若粒子在某維越界前的取值與邊界的差值很小，則判定粒子原本就處于邊界附近，采用吸收墻辦法，將其值設(shè)為邊界值，反之判定粒子在該維的更新無效，粒子在該維與當(dāng)前群體歷史最優(yōu)解進行交叉操作，具體越界處理公式為：

2.2.4 DE-PSO算法步驟 ①初始化粒子位置、個體極值、子種群極值及群體極值。②分別對兩個子種群的粒子進行排序，按式(8)執(zhí)行選擇、交叉操作(完成一次迭代)。③按式(6)—(7)更新粒子位置(進入新的迭代)，如有越界按式(9)—(10)處理。④評價粒子，更新個體極值、子種群極值及群體極值。⑤判斷是否滿足迭代終止條件(或達到最大迭代次數(shù))，如滿足，算法終止，輸出群體極值，否則轉(zhuǎn)②繼續(xù)優(yōu)化。

3 算例研究

3.1 水頭約束問題

3.1.1 問題概述 算例改編自文獻(于福榮等，2010)，如圖1所示，某均質(zhì)、各項同性承壓含水層中擬設(shè)9口抽水井用于疏干建設(shè)區(qū)內(nèi)的地下水，使得建設(shè)區(qū)范圍內(nèi)的地下水位降至50 m以下。含水層底板水平;南北部為隔水邊界;東西部為一類邊界，給定水頭值分別為60 m和80 m;地下水流為二維穩(wěn)定流，含水層導(dǎo)水系數(shù)為 50 m2/d，采用MODFLOW模擬水流。每口井的最大抽水量為3 000 m3/d。用邊長為100 m的正方形網(wǎng)格將研究區(qū)剖分為20行、30列的網(wǎng)格。

圖1 水頭約束問題研究區(qū)示意圖

約束條件：

式中，J為目標(biāo)函數(shù)，為總抽水量(m3/d);Qi為井i的抽水量(m3/d);Hj表示水頭約束點j的水頭;hk+1表示研究區(qū)地下水在k+1時段的水頭分布;Transh代表水頭由k時段到k+1時段的狀態(tài)轉(zhuǎn)移函數(shù);Q為決策變量，Q=［Q1，Q2，Q3，…，Q9］。

分別用GA、SA、PSO及筆者提出的DE-PSO求解算例，其中PSO對約束條件及越界的處理與DE-PSO一致。GA參數(shù)為：編碼長度45，交叉概率0.5，變異概率0.01，種群規(guī)模100，最大迭代次數(shù)100。SA參數(shù)為：每個變量的離散化區(qū)間數(shù)為32，鄰域搜索決策變量單步長中增加量為1，初始溫度e10，衰減因子0.75，每個T值的迭代次數(shù)為100，外循環(huán)最大迭代次數(shù)為100。PSO參數(shù)為：W1=W2=2.0，m=0.5，K=0.1，β =0.01，種群規(guī)模 100，最大迭代次數(shù)為100。DE-PSO參數(shù)為：W1=W2=2.0，α=0.5，K=0.1，β=0.01，兩個子種群規(guī)模均為50，最大迭代次數(shù)為100。

3.1.2 優(yōu)化結(jié)果對比 表1為優(yōu)化結(jié)果對比表，其中，與GA、SA及PSO解相比，DE-PSO解得到的總抽水量最少，分別比GA、SA和PSO減少了64、256、207 m3/d。將幾種智能算法解代入水流模型，發(fā)現(xiàn)所有解均滿足約束條件。圖2為相應(yīng)的地下水位等值線優(yōu)化結(jié)果圖，其中，DE-PSO地下水位等值線圖中50 m等水位線所包圍的區(qū)域集中在建設(shè)區(qū)附近。圖3對比了各智能算法運行時，10～100代最優(yōu)解總抽水量的變化趨勢。由圖3可知，PSO進化緩慢，近于停滯;GA、SA在相同迭代次數(shù)下的總抽水量總體高于DE-PSO;DE-PSO曲線呈較穩(wěn)定的下降趨勢，尋優(yōu)效率最好，分析其原因在于DE-PSO引入的選擇、交叉及差分變異操作避免了PSO的早熟現(xiàn)象。

表1 水頭約束問題優(yōu)化結(jié)果對比

3.2 質(zhì)量分數(shù)約束問題

3.2.1 問題概述 算例選自文獻(徐星等，2011)如圖4所示，一承壓含水層存在硝酸鹽(以氮計)污染，擬用抽出處理方法治理，治理周期為5 a，4口預(yù)選治理抽水井的位置已在圖中標(biāo)出。水流模型為二維非穩(wěn)定流模型。研究區(qū)南北部為隔水邊界，東西部為定水頭邊界，水頭值分別為25 m和35 m;研究區(qū)西部、北部和南部為溶質(zhì)零通量邊界，東部為溶質(zhì)對流給定通量邊界。用邊長為150 m的正方形網(wǎng)格將研究區(qū)剖分為17行23列有限差分網(wǎng)格，采用MODFLOW模擬流場，MT3DMS模擬質(zhì)量分數(shù)場。模型設(shè)1個應(yīng)力期，5個時間步長，每個時間步長進一步被分解為若干個運移步長。模型用到的水文地質(zhì)參數(shù)：孔隙度為0.2，導(dǎo)水系數(shù)為500 m2/d，貯水系數(shù)為0.000 1，縱向彌散度為50 m，橫向彌散度為10 m。管理模型為：

約束條件：

圖2 水頭約束問題地下水位等值線優(yōu)化結(jié)果圖

圖3 水頭約束問題各智能算法的迭代運行

圖4 含水層PAT系統(tǒng)設(shè)計的結(jié)構(gòu)示意圖

式中，J為目標(biāo)函數(shù)，為總治理成本;NA為非零流量井的數(shù)目;N為總的優(yōu)選井的數(shù)量;a1為安裝井的費用系數(shù)，算例中單口井安裝費用10萬元;a2為處理費用系數(shù)，算例中取0.8元/m3;Qi為井i的抽水量(m3/d);cm是治理周期末體積質(zhì)量分數(shù)約束區(qū)域任意一點的溶質(zhì)體積質(zhì)量分數(shù)(mg/L);c*為治理周期結(jié)束后體積質(zhì)量分數(shù)約束區(qū)域的體積質(zhì)量分數(shù)約束條件(mg/L)，算例中設(shè)置為20 mg/L，即達到地下水III類水標(biāo)準(zhǔn);hk+1、ck+1分別表示研究區(qū)地下水在k+1時段的水頭和溶質(zhì)體積質(zhì)量分數(shù)分布;Transh，Transc分別代表水頭和溶質(zhì)體積質(zhì)量分數(shù)由k時段到k+1時段的狀態(tài)轉(zhuǎn)移函數(shù);Q為決策變量，Q=［Q1，Q2，Q3，Q4］。

分別用GA、SA、PSO及筆者提出的DE-PSO求解算例，其中PSO對約束條件及越界的處理與DE-PSO一致。GA參數(shù)為：編碼長度20，交叉概率0.5，變異概率0.025，種群規(guī)模40，最大迭代次數(shù)50;SA參數(shù)為：每個變量的離散化區(qū)間數(shù)為32，每個T值的迭代次數(shù)40，外循環(huán)最大迭代次數(shù)50，其他參數(shù)值同算例1;PSO參數(shù)為：種群規(guī)模40，最大迭代次數(shù)50，其他參數(shù)值同算例1;DE-PSO參數(shù)為：兩個子種群規(guī)模均為20，最大迭代次數(shù)50，其他參數(shù)值同算例1。

3.2.2 優(yōu)化結(jié)果對比 表2為優(yōu)化結(jié)果對比。其中，DE-PSO解得到的總治理成本最少，分別比GA、SA和PSO減少了57.74、151.93、76.59萬元。

表2 濃度約束問題優(yōu)化結(jié)果對比

圖5對比了各智能算法運行時，10～50代最優(yōu)解總治理成本的變化趨勢，其中GA、SA、PSO在相同迭代次數(shù)下的總治理成本總體高于DE-PSO，DE-PSO曲線呈較穩(wěn)定的下降趨勢，尋優(yōu)效率最好，避免了PSO的早熟現(xiàn)象。

圖5 質(zhì)量分數(shù)約束問題各智能算法的迭代運行

4 結(jié)論

(1)基于PSO算法，設(shè)立了兩個平行子種群進行優(yōu)化，增加選擇、交叉操作，在粒子更新公式中引入差分變異操作，并采用有效的粒子評價、越界處理方法，提出DE-PSO算法。

(2)算例模擬結(jié)果表明：DE-PSO對PSO作改進后，增強了全局尋優(yōu)能力且進化穩(wěn)定，克服了PSO的早熟缺點，優(yōu)化結(jié)果及尋優(yōu)效率優(yōu)于GA、SA及PSO，適于求解地下水管理模型。

李相勇，田澎，孔民.2007.解約束優(yōu)化問題的新粒子群算法［J］.系統(tǒng)管理學(xué)報，16(2)：120-129.

劉衍民，隋常玲，牛奔.2011.解決約束優(yōu)化問題的改進粒子群算法［J］.計算機工程與應(yīng)用，47(12)：23-26.

羅德相.2009.粒子群算法改進方法研究［D］.廣西：廣西民族大學(xué).

汪定偉，王俊偉，王洪峰，等.2007.智能優(yōu)化方法［M］.北京：高等教育出版社.

吳蓉，周志芳.2002.改進的前饋型神經(jīng)網(wǎng)絡(luò)模型在承壓水漏斗水位預(yù)報中的應(yīng)用［J］.地質(zhì)學(xué)刊(原《江蘇地質(zhì)》)，26(1)：19-21.

徐星，李元香，吳昱.2010.基于擴散機制的雙種群粒子群優(yōu)化算法［J］.計算機應(yīng)用研究，27(8)：2882-2898.

徐星，吳昱，李元香.2011.基于布朗運動的改進粒子群算法［J］.計算機應(yīng)用研究，28(7)：2439-2442.

楊蘊，吳劍鋒，吳吉春.2009.兩種智能算法在求解地下水管理模型中的對比［J］.吉林大學(xué)學(xué)報：地球科學(xué)版，39(3)：474-481.

于福榮，盧文嘉，李平，等.2010.地下水管理模型求解方法綜述［J］.水文地質(zhì)工程地質(zhì)，37(2)：27-31.

AHLFELD D P，MULLIGAN A E.2000.Optimal management of flow in groundwater systems［M］.San Diego：Academic Press.

MCKINNEY D C，LIN M D.1995.Mixed-integer nonlinear programming methods for optimal aquifer remediation design［J］.Water Resources Research，31(3)：731-740.

WANGER B J.1995.Recent advances in simulation optimization groundwater management modeling［J］.Review of Geophysics，(Suppl)：1021-1028.

Double-population evolution-particle swarm optimization algorithm for solving groundwater management model

WU Rui-qi，ZHU Guo-rong，WANG Pei

(School of Earth Sciences and Engineering，Nanjing University，Nanjing 210093，China)

To solve the premature convergence problem of Particle Swarm Optimization(PSO)，a new algorithm named DE-PSO(Double-population Evolution-Particle Swarm Optimization)was designed.DE-PSO introduced selection，crossover and differential mutation into PSO，and adopted a new evaluation method to evaluate swarms and a new control method to ensure all swarms can fly inside search space.DE-PSO was applied to solve two groundwater management cases.In the first case，DE-PSO produced a design with respectively 64，256 and 207 m3/d less pumping rate than those of GA 、SA and PSO;in the second case，DE-PSO produced a design with respectively ￥577，400，￥1，519，300 and ￥765，900 less remedial design cost than those of GA，SA and PSO.Two case studies indicated that DE-PSO could evolve steadily，the searching efficiency was better than GA，SA and PSO.DE-PSO could solve groundwater management model problems effectively.

Groundwater management model;Particle swarm optimization;Double population;Differential variation

TP18;P641.8

A

1674-3636(2012)01-0037-07

10.3969/j.issn.1674-3636.2012.01.37

2011-10-26;編輯：陸李萍

國家自然科學(xué)基金項目(J0830522)資助

吳睿奇(1988— )，女，碩士研究生，水文與水資源專業(yè)，E-mail：gwflow@nju.edu.cn