王曉晨，常安定，王靜云，張 轉(zhuǎn)(長(zhǎng)安大學(xué) 理學(xué)院，陜西 西安 710064)

各向異性含水層參數(shù)估計(jì)的差分進(jìn)化算法

王曉晨，常安定，王靜云，張 轉(zhuǎn)(長(zhǎng)安大學(xué) 理學(xué)院，陜西 西安 710064)

【目的】 通過(guò)分析抽水試驗(yàn)數(shù)據(jù)進(jìn)行各向異性含水層參數(shù)估計(jì)，為各向異性含水層參數(shù)估計(jì)提供新的方法。【方法】 以各向同性情況下的泰斯公式為基礎(chǔ),將差分進(jìn)化算法應(yīng)用于各向異性條件下含水層參數(shù)的估計(jì)，并將其計(jì)算結(jié)果與其他方法的參數(shù)估計(jì)結(jié)果進(jìn)行比較，對(duì)不同初始參數(shù)范圍下的參數(shù)估計(jì)值進(jìn)行分析探討?！窘Y(jié)果】 差分進(jìn)化算法估計(jì)的參數(shù)值與其他方法估計(jì)的參數(shù)值之間有較小的差異，相對(duì)誤差不到4%，且目標(biāo)函數(shù)值相對(duì)更小，達(dá)到0.003 93；對(duì)于不同的初始參數(shù)范圍，利用該算法均能得到滿意的參數(shù)估計(jì)結(jié)果?！窘Y(jié)論】 基于抽水試驗(yàn)數(shù)據(jù)估計(jì)各向異性含水層參數(shù)的差分進(jìn)化算法估計(jì)結(jié)果有效且可靠，算法穩(wěn)定性高，尋優(yōu)能力強(qiáng)。

各向異性；含水層參數(shù)；差分進(jìn)化算法；參數(shù)估計(jì)；抽水試驗(yàn)

在研究與地下水運(yùn)動(dòng)有關(guān)的工程問(wèn)題中，含水層的導(dǎo)水性能具有各向異性的特征，且含水層導(dǎo)水系數(shù)和儲(chǔ)水系數(shù)作為最重要的參數(shù)，其準(zhǔn)確有效地求解極為重要。為此許多學(xué)者在利用抽水試驗(yàn)確定各向異性含水層參數(shù)問(wèn)題上做了很多的研究工作。Hantush[1]給出了有越流補(bǔ)給條件下，徑向各向異性含水層中非完整井流三維問(wèn)題的數(shù)學(xué)模型及相應(yīng)的計(jì)算公式。Papadopulos[2]推導(dǎo)出了徑向各向異性主值Kr與2個(gè)水平方向主滲透系數(shù)Kx、Ky大小和方向的關(guān)系式。之后，Way等[3]基于Hantush[1]和Papadopulos[2]的研究結(jié)果，提出了通過(guò)現(xiàn)場(chǎng)抽水試驗(yàn)確定徑向各向異性滲透系數(shù)主值(Kr、Kz和Kx、Ky)的方法。Neuman等[4]利用張量分析的方法推導(dǎo)出了平面各向異性滲透系數(shù)主值和主方向的解析公式。Chen等[5]討論了Hantush解中越流的影響及越流系數(shù)的計(jì)算方法。周志芳等[6]通過(guò)抽水試驗(yàn)，應(yīng)用圖解法和優(yōu)化算法結(jié)合的半解析方法確定平面各向異性滲透系數(shù)主值和主方向，之后又推導(dǎo)出了有越流情況下各向異性含水層中地下水三維井流問(wèn)題的計(jì)算公式[7]。劉燕等[8]以描述均質(zhì)各向同性與各向異性情況下的非穩(wěn)定流問(wèn)題的解析解為基礎(chǔ)，在滿足直線圖解法適用條件的情況下，分析了3組實(shí)際抽水試驗(yàn)數(shù)據(jù)，計(jì)算了各向同性與各向異性2種情況下的含水層參數(shù)。近年來(lái)，智能優(yōu)化算法在經(jīng)濟(jì)和工程等各個(gè)領(lǐng)域有著廣泛的應(yīng)用。康瑞龍[9]于2013年成功地應(yīng)用智能優(yōu)化算法——人工根系算法求解了各向異性含水層參數(shù)問(wèn)題。區(qū)別于傳統(tǒng)的估計(jì)方法，智能優(yōu)化算法求解速度快而且計(jì)算精度高。差分進(jìn)化算法(Differential evolution,DE)是一種隨機(jī)并行直接搜索的智能優(yōu)化算法[10]，該算法容易實(shí)現(xiàn)，對(duì)參數(shù)的依賴程度低。鑒于此，本研究根據(jù)抽水試驗(yàn)數(shù)據(jù)，嘗試將差分進(jìn)化算法應(yīng)用于確定各向異性條件下含水層參數(shù)的計(jì)算問(wèn)題，并對(duì)算法和計(jì)算結(jié)果進(jìn)行分析，以期為各向異性條件下含水層參數(shù)的準(zhǔn)確估計(jì)提供方法上的支持。

1 各向異性條件下的含水層參數(shù)問(wèn)題

1.1 抽水試驗(yàn)基本知識(shí)

在抽水試驗(yàn)中，假設(shè)在無(wú)限延伸含水層中，任一點(diǎn)處的水位降深隨時(shí)間t變化用各向異性表達(dá)式表示：

(1)

式中：s為水位降深，m；Q為抽水井中的抽水流量，m3/s；Te可以定義為含水層的等效導(dǎo)水系數(shù)，m2/s；W(uxy)為與泰斯公式[11]形式相同的井函數(shù)。Te的表達(dá)式為：

(2)

式中：Txx和Txy分別為導(dǎo)水系數(shù)在當(dāng)?shù)刈鴺?biāo)系x和y方向的張量分量，m2/s；函數(shù)W(uxy)中無(wú)量綱時(shí)間為：

(3)

式中：μ為含水層儲(chǔ)水系數(shù)；Tyy為導(dǎo)水系數(shù)在當(dāng)?shù)刈鴺?biāo)系下的張量分量，m2/s；x，y為當(dāng)?shù)刈鴺?biāo)系的坐標(biāo)分量，m。

在此，采用Srivastava[11]的泰斯井函數(shù)的近似表達(dá)式進(jìn)行計(jì)算，具體表達(dá)式如下。W(u)=-lnu+a0+a1u+a2u2+a3u3+a4u4+a5u5，u≤1。

(4)

(5)

其中，a0=-0.577 22，a1=0.999 99，a2=-0.249 91，a3=0.055 19，a4=-0.009 76，a5=0.001 08；b0=0.267 77，b1=8.634 76，b2=18.059 02，b3=8.573 33；c0=3.958 50，c1=21.099 65，c2=25.632 96，c3=9.573 32。Srivastava[11]指出，以上的近似表達(dá)式與泰斯公式精確解相比，相對(duì)誤差均小于 0.001%。

(6)

(7)

式中：X和Y為導(dǎo)水系數(shù)張量的全局坐標(biāo)。

1.2 目標(biāo)函數(shù)的建立

通常情況下，采用3個(gè)不在同一觀測(cè)線上的觀測(cè)孔中的降深-時(shí)間數(shù)據(jù)確定當(dāng)?shù)刈鴺?biāo)系下的導(dǎo)水系數(shù)張量分量Txx、Tyy和Txy的值，然后再利用式(6)和(7)計(jì)算全局坐標(biāo)系下的導(dǎo)水系數(shù)張量分量TX和TY值。因此，對(duì)于各向異性含水層參數(shù)的估計(jì)問(wèn)題，可以將含水層儲(chǔ)水系數(shù)μ及導(dǎo)水系數(shù)在當(dāng)?shù)刈鴺?biāo)系下的張量分量Txx、Tyy和Txy設(shè)定為待估參數(shù)。

對(duì)于第i(i=1，2，3)個(gè)觀測(cè)孔，其tj時(shí)刻的水位降深sj可由式(1)求得。為了使所求得的sj與真實(shí)水位降深相接近，所以構(gòu)造目標(biāo)函數(shù)：

(8)

對(duì)于各向異性條件下的含水層，u、Txx、Tyy、Txy為待估參數(shù)，其向量為θ=(u,Txx,Tyy,Txy)。為此，使得目標(biāo)函數(shù)最小的參數(shù)值θ=(μ,Txx,Tyy,Txy)即為所求。

2 差分進(jìn)化算法

2.1 DE算法的思想

DE算法是由Storn等[12]于1995年提出的,它是一種隨機(jī)并行直接搜索算法，該算法的基本思想是應(yīng)用當(dāng)前種群個(gè)體的差異重組得到中間種群，然后應(yīng)用子代個(gè)體與父代個(gè)體競(jìng)爭(zhēng)獲得新一代種群[13]。算法中采用由變異到交叉再到選擇的操作順序更新下一代種群。其中，變異操作是差異演化的關(guān)鍵步驟，交叉操作可以增加種群的多樣性，選擇決定著試驗(yàn)向量是否會(huì)成為下一代中的成員，確保下一代中的所有個(gè)體都比當(dāng)前種群的對(duì)應(yīng)個(gè)體更佳或者至少某一方面要好[14]。

2.2 DE算法的步驟及流程

步驟1：隨機(jī)初始化種群，設(shè)置相關(guān)參數(shù)。

步驟2：計(jì)算各個(gè)體最優(yōu)適應(yīng)度f(wàn)i，及其相應(yīng)的位置gi，從中找出當(dāng)代全局最優(yōu)適應(yīng)度f(wàn)及其相應(yīng)的位置g。

步驟4：若滿足終止條件，則迭代停止；否則，k=k+1，轉(zhuǎn)步驟3。

步驟5：輸出全局最優(yōu)適應(yīng)度f(wàn)、相應(yīng)的位置g和迭代次數(shù)t。

差分進(jìn)化算法流程圖如圖1所示。

3 數(shù)值試驗(yàn)

3.1 試驗(yàn)數(shù)據(jù)與試驗(yàn)條件

3.1.1 試驗(yàn)數(shù)據(jù) 試驗(yàn)數(shù)據(jù)引自文獻(xiàn)[15]。假定在全局某正交各向異性含水層中，有一完整井以定流量Q=0.012 57 m3/s進(jìn)行非穩(wěn)定流抽水試驗(yàn)，有3個(gè)觀測(cè)孔OW1、OW2和OW3，其在坐標(biāo)系中的位置分別為(28.3，0)、(9.0，33.5)、(-19.3，-5.2)。當(dāng)抽水井以定流量非穩(wěn)定流抽水后，在3個(gè)觀測(cè)孔中觀測(cè)到的地下水降深隨時(shí)間變化的試驗(yàn)數(shù)據(jù)如表1所示。

圖1 差分進(jìn)化算法流程圖Fig.1 Flow diagram of DE

3.1.2 試驗(yàn)條件 初始種群規(guī)模為50；最大迭代次數(shù)2 000；達(dá)到精度eps=0.004；參數(shù)限定范圍參考文獻(xiàn)[9]：0≤μ≤1，0≤Txx≤600，0≤Tyy≤600，-600≤Txy≤0。對(duì)于差分進(jìn)化算法中的參數(shù)CR的選取與文獻(xiàn)[16]不同，其選取原則是使差分進(jìn)化算法在解決各向異性條件下含水層參數(shù)問(wèn)題具有較好的全局搜索能力，在此選取CR=0.8。

將差分進(jìn)化算法程序反復(fù)執(zhí)行100次，如果計(jì)算的目標(biāo)函數(shù)值小于精度0.004，認(rèn)為算法成功，即利用差分進(jìn)化算法成功地估計(jì)了各向異性條件下的含水層參數(shù)；否則，認(rèn)為算法不成功。記錄算法每次的運(yùn)行時(shí)間和算法成功的次數(shù)。

對(duì)于上述試驗(yàn)數(shù)據(jù)和試驗(yàn)條件，采用Matlab語(yǔ)言，根據(jù)算法流程實(shí)現(xiàn)DE算法。

表1 觀測(cè)孔時(shí)間與降深數(shù)據(jù)Table 1 The time-drawdown data of the observation

3.2 試驗(yàn)結(jié)果與分析

3.2.1 試驗(yàn)結(jié)果 應(yīng)用差分進(jìn)化算法估計(jì)各向異性條件下的含水層參數(shù)，將其計(jì)算結(jié)果(迭代第80次)與相關(guān)文獻(xiàn)中的結(jié)果進(jìn)行對(duì)比，結(jié)果如表2所示。

表2 基于差分進(jìn)化算法的各向異性含水層參數(shù)估計(jì)結(jié)果(迭代第80次)與相關(guān)文獻(xiàn)結(jié)果的比較Table 2 The results (Running the 80th time) of the anisotropy aquifer parameters based on DE compare with the results of the related literature

3.2.2 算法有效性 由表2可以看出，在將差分進(jìn)化算法應(yīng)用于估計(jì)各向異性條件下含水層參數(shù)時(shí)，程序迭代到80次時(shí)達(dá)到精度要求，其目標(biāo)函數(shù)值為0.003 93，計(jì)算出的相應(yīng)參數(shù)與文獻(xiàn)[15]中的參數(shù)相差較小，且計(jì)算出的相應(yīng)參數(shù)間的相對(duì)誤差均小于4%。因此，差分進(jìn)化算法不論對(duì)待估參數(shù)μ、Txx、Tyy、Txy，還是通過(guò)式(6)、式(7)計(jì)算的TX和TY，均能得到可靠的計(jì)算結(jié)果。

圖2是差分進(jìn)化算法程序執(zhí)行100次過(guò)程中，目標(biāo)函數(shù)值的分布情況。由圖2可以看出，在差分進(jìn)化算法程序執(zhí)行100次中，有93次成功，其尋優(yōu)率達(dá)到93%，平均運(yùn)行時(shí)間為1.928 35 s。由此可知，差分進(jìn)化算法的尋優(yōu)效率較高且運(yùn)行時(shí)間較短，具有高效性。因此，運(yùn)用差分進(jìn)化算法能有效地估計(jì)各向異性條件下的含水層參數(shù)。

圖2 差分進(jìn)化算法執(zhí)行100次過(guò)程中目標(biāo)函數(shù)值的分布Fig.2 Distribution of the objective function value when DE is executed 100 times

3.2.3 算法穩(wěn)定性 表3的試驗(yàn)結(jié)果是在算法參數(shù)設(shè)定不變，參數(shù)限定范圍為文獻(xiàn)[15]的不同倍數(shù)的情況下，應(yīng)用差分進(jìn)化算法得到的參數(shù)估計(jì)值。從表3中的時(shí)間數(shù)據(jù)可以看出，隨著初始狀態(tài)參數(shù)范圍的不斷增大，差分進(jìn)化算法均能較好地收斂，得到較為滿意的計(jì)算結(jié)果，即差分進(jìn)化算法對(duì)初值選取的敏感性較低，尋優(yōu)能力較強(qiáng)，具有較好的穩(wěn)定性。在初始狀態(tài)參數(shù)范圍倍數(shù)增大的情況下，相應(yīng)的所需運(yùn)行時(shí)間有所增加，但是時(shí)間波動(dòng)不大，其增加量不超過(guò)0.3 s，且運(yùn)行時(shí)間較短，這說(shuō)明差分進(jìn)化算法能以較短的時(shí)間收斂于全局值。因此，差分進(jìn)化算法解決各向異性條件下含水層參數(shù)問(wèn)題是穩(wěn)定的。

表3 參數(shù)初始范圍為文獻(xiàn)[15]參數(shù)值的不同倍數(shù)時(shí)DE算法對(duì)各向異性含水層參數(shù)的計(jì)算結(jié)果Table 3 The calculation results of DE under the initial range of parameters are setted to multiple the parameters of the literature [15]

3.2.4 算法收斂性 在不限制算法精度的條件下，讓算法以迭代2 000次為終止條件，可以得到算法在執(zhí)行189次時(shí)會(huì)收斂到0.003 402 81(圖3)，由此可知差分進(jìn)化算法有較快的收斂速度。

圖3 差分進(jìn)化算法執(zhí)行2 000次時(shí)的收斂曲線Fig.3 Convergent curve of DE when DE is executed 2 000 times

4 結(jié)束語(yǔ)

鑒于各向異性含水層參數(shù)估計(jì)的復(fù)雜性，本文嘗試應(yīng)用差分進(jìn)化算法對(duì)各參數(shù)進(jìn)行估計(jì)。根據(jù)試驗(yàn)結(jié)果可以知道，用差分進(jìn)化算法能夠成功確定各向異性條件下的含水層參數(shù)，并且可以得到可靠的優(yōu)化計(jì)算結(jié)果。本試驗(yàn)結(jié)果表明：(1)差分進(jìn)化算法估計(jì)各向異性含水層參數(shù)時(shí)，收斂速度快，運(yùn)行時(shí)間短，精度高而且具有較好的尋優(yōu)率；(2)算法對(duì)不同的初值范圍均能得到滿意的參數(shù)計(jì)算結(jié)果，具有較高的穩(wěn)定性和較強(qiáng)的尋優(yōu)能力，且對(duì)初值的敏感程度較低。因此，差分進(jìn)化算法可以成功、有效地用于求解各向異性含水層的參數(shù)問(wèn)題。

[1] Hantush M S.A method for analyzing a drawdown test in anisotropic aquifers [J].Water Resources,1966,2(2):281-285.

[2] Papadopulos I S.Nonsteady flow to a well in an infinite anisotropic aquifer [C]//United Nations Educational Scientific and Cultural Organization.Symposium on Hydrology of Fractured Rocks,Paris:U.S.Geological Survey,1965:21-31.

[3] Way S C,Mckee C R.In-situ determination of three-dimensional aquifer permeabilities [J].Ground Water,1982,20(5):594-603.

[4] Neuman S P,Walter G R,Bentley W,et al.Determination of horizontal anisotropy with three wells [J].Ground Water,1984,22(1):66-72.

[5] Chen X,Ayers J F.Utilization of the Hantush solution for the simultaneous determination of aquifer paramenter [J].Ground Water,1997,35(5):751-756.

[6] 周志芳,朱學(xué)愚,李 艷.巖體滲透系數(shù)張量的半解析計(jì)算 [J].水利學(xué)報(bào),1997(9):06-11. Zhou Z F,Zhu X Y,Li Y.Semi-analytic method for conductivity tensor in rock mass [J].Journal of Hydraulic Engineering,1997(9):6-11.(in Chinese)

[7] 周志芳.任意各向異性巖體滲透系數(shù)張量的半解析計(jì)算 [J].水利學(xué)報(bào),1999(3):65-70. Zhou Z F.Semi-analytic method for conductivity tensor in random anisotropic rock mass [J].Journal of Hydraulic Engineering,1999(3):65-70.(in Chinese)

[8] 劉 燕,辛璐軍,郭建青，等.抽水試驗(yàn)確定各向異性含水層參數(shù)的實(shí)例討論 [J].勘察科技技術(shù),2012(6):5-9. Liu Y,Xin L J,Guo J Q,et al.Discussion of determining anisotropic aquifer parameters by pumping tests data [J].Site Investigation Science and Technology,2012(6):5-9.(in Chinese)

[9] 康瑞龍.一種新型的智能優(yōu)化算法：人工根系算法 [D].西安:長(zhǎng)安大學(xué),2013. Kang R L.A new intelligent optimization method:Artificial root optmization [D].Xi’an:Chang’an University,2013.(in Chinese)

[10] Stroll R,Price K.Differentifial evolution:A survey of the state of the art [J].IEEE Transaction on Evolutionary Computation,2011,15(1):4-31.

[11] Srivastava R.Implications of using approximation expressions for well function [J].Journal of Irrigation and Drainage Engineering,1995,121(6):459-462.

[12] Storn R,Price K.Differential evolution:A simple and efficient adaptive scheme for global optimization over continuous spaces [M].Berkeley:ICSI,1995:22-25.

[13] 萬(wàn) 東.差分進(jìn)化算法研究及其應(yīng)用 [J].科學(xué)技術(shù)與工程,2009,9(22):6672-6676. Wan D.Research and application based on differential evolution algorithm [J].Science Technology and Engineering,2009,9(22):6672-6676.(in Chinese)

[14] Storn R,Price K.Differential evolution:A simple and efficient heuristic for global optimization over continuous spaces [J].Journal of Global Optimization,1997,11(4):341-359.

[15] Kruseman G P,Ridder N A,Verweij J M.Analysis and evaluation of pumping test data [M].2ndEdition.Wageningen Netherlands:International Institute for Land Reclamation and Improvement,1994:133-145.

[16] 周艷平,顧幸生.差分進(jìn)化算法研究進(jìn)展 [J].化工自動(dòng)化及儀表,2007,34(3):1-5. Zhou Y P,Gu X S.Development of differential evolution algorithm [J].Chemical Industry Automation and Instrumentation,2007,34(3):1-5.(in Chinese)

Differential evolution method for estimating the anisotropy aquifer parameters

WANG Xiao-chen,CHANG An-ding,WANG Jing-yun,ZHANG Zhuan

(CollegeofScience,ChanganUniversity,Xi’an,Shaanxi710064,China)

【Objective】 This study estimated the anisotropy aquifer parameters by the analysis of pumping test data. In order to provide a new method for estimating the parameter under the condition of anisotropy aquifer.【Method】 On the basis of the theis formula in the isotropic condition,differential evolution algorithm was applied to estimate the parameter under the condition of anisotropy aquifer.On the one hand,the calculation compare with other methods of parameter estimation and analyzes the results.On the other hand,discussing the estimated value of parameter that calculate under different initial scope.【Result】 The relative error that the values of the parameters are estimated between differential evolution algorithm and other method and the objective function value are relatively smaller.The relative error less than 4% and the objective function value is 0.003 93;For different initial parameters,differential evolution algorithm can get a satisfactory parameter estimation result.【Conclusion】 Based on pumping test data,the results of differential evolution method are more effective and reliable.It has higher stability and the stronger ability of optimization than other methods.

random anisotropic;aquifer parameters;differential evolution;parameter estimate;pumping test

2013-12-09

中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金(310829130225)

王曉晨(1988-),女,遼寧鞍山人,在讀碩士,主要從事最優(yōu)化理論與方法研究。E-mail:wangxiaochenbest@126.com

常安定(1964-),男,陜西大荔人,教授,碩士生導(dǎo)師,主要從事水文地質(zhì)的數(shù)學(xué)分析方法研究。E-mail:chdanding@126.com

時(shí)間:2015-03-12 14:17

10.13207/j.cnki.jnwafu.2015.04.023網(wǎng)絡(luò)出版地址:http://www.cnki.net/kcms/detail/61.1390.S.20150312.1417.023.html

1671-9387(2015)04-0223-06

TV211.1+2

A