曹懷火，歐陽(yáng)艾嘉，艾海男

1.池州學(xué)院 數(shù)學(xué)與計(jì)算機(jī)科學(xué)系，安徽 池州 247000

2.湖南大學(xué) 計(jì)算機(jī)與通信學(xué)院，長(zhǎng)沙 410082

3.重慶大學(xué) 三峽庫(kù)區(qū)生態(tài)環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，重慶 400045

1 引言

土壤水分特征曲線(xiàn)是表征土壤含水率與土壤水吸力的關(guān)系曲線(xiàn)，是定量研究土壤水分運(yùn)動(dòng)的重要參考。影響土壤水分特征曲線(xiàn)的因素非常復(fù)雜，通常情況下難以從理論上得出確定的關(guān)系式。因此，人們尋求一些經(jīng)驗(yàn)公式來(lái)數(shù)值模擬土壤水分特征曲線(xiàn)，目前國(guó)內(nèi)外最為普遍使用的描述土壤水分特征曲線(xiàn)的Van Genuchten方程[1]（簡(jiǎn)稱(chēng)VG方程）是美國(guó)學(xué)者Van Genuchten在1980年提出的，其具體表達(dá)式為：

其中θ為土壤含水率（m3/cm3）；h為土壤水吸力（cm）；θs為土壤飽和含水率（cm3/cm3）；θr為土壤殘余含水率（cm3/cm3）；α,m,n為土壤水分特征曲線(xiàn)形狀參數(shù)，m=1-1/n，n>1。

通常對(duì)于求解VG方程式（1）的參數(shù)估計(jì)問(wèn)題采用間接法，即轉(zhuǎn)化為如下最優(yōu)化模型（2）的求解問(wèn)題：

其中θi為實(shí)測(cè)土壤含水率；hi為實(shí)測(cè)土壤基質(zhì)勢(shì)；N為實(shí)測(cè)數(shù)據(jù)組數(shù)。傳統(tǒng)的優(yōu)化方法[2-5]大多基于梯度計(jì)算，計(jì)算效率比較高，但由于其固有的局部?jī)?yōu)化性，并不適合于全局優(yōu)化問(wèn)題的求解。于是很多學(xué)者對(duì)求解VG方程的參數(shù)估計(jì)問(wèn)題引入了智能進(jìn)化算法，如遺傳算法、粒子群算法等[6-8]，同時(shí)也取得了一定的效果。

由于遺傳算法和粒子群算法受初值影響，易于產(chǎn)生局部最優(yōu)解等缺陷，目前很多學(xué)者給出了一些改進(jìn)的遺傳算法[9]和改進(jìn)的粒子群算法[10]。為了提高粒子群算法的收斂概率和減少迭代次數(shù)，綜合文獻(xiàn)[11-12]方法的思想，將lsqcurvefit擬合方法與粒子群優(yōu)化算法相結(jié)合，構(gòu)造一種新的混合型粒子群優(yōu)化算法，并用于求解VG方程的參數(shù)估計(jì)問(wèn)題。該算法不僅包含粒子群算法的全局搜索能力，還具有l(wèi)sqcurvefit擬合函數(shù)較強(qiáng)的局部?jī)?yōu)化能力。

2 算法構(gòu)造

2.1 粒子群優(yōu)化（PSO）算法簡(jiǎn)介

粒子群優(yōu)化算法是一種進(jìn)化計(jì)算技術(shù)，它是模仿鳥(niǎo)群或魚(yú)群的社會(huì)行為而提出的算法。首先初始化一群隨機(jī)粒子，然后粒子們就追隨當(dāng)前的最優(yōu)粒子在解空間中搜索，即通過(guò)迭代找到最優(yōu)解。假設(shè)d維搜索空間中的第i個(gè)粒子的位置和速度分別為Xi=(xi1,xi2,…,xid)和Vi=(vi1,vi2,…,vid)，在每一次迭代中，粒子通過(guò)跟蹤兩個(gè)最優(yōu)解來(lái)更新自己，第一個(gè)就是粒子本身所找到的最優(yōu)解，即個(gè)體極值pbest，Pi=(pi1,pi2,…,pid)；另一個(gè)是整個(gè)種群目前找到的最優(yōu)解，即全局最優(yōu)解gbest，在找到這兩個(gè)最優(yōu)解時(shí)，粒子根據(jù)如下公式來(lái)更新自己的速度和新的位置：

其中w為慣性權(quán)因子，c1和c2為正的學(xué)習(xí)因子，r1和r2為0到1之間均勻分布的隨機(jī)數(shù)。

2.2 基于Lsqcurvefit擬合的粒子群優(yōu)化算法

PSO算法的搜優(yōu)速度在一定程度上快于其他進(jìn)化算法，但卻存在全局搜索能力極強(qiáng)而局部尋優(yōu)能力比較差，也就是說(shuō)該算法可以用于極快的速度搜索到最優(yōu)解附近，但要進(jìn)一步達(dá)到最優(yōu)解則速度很慢。為了加快局部尋優(yōu)能力，將lsqcurvefit擬合方法與PSO算法相結(jié)合，構(gòu)造一種新的混合型粒子群優(yōu)化算法（LPSO算法），其基本思路是在進(jìn)化過(guò)程當(dāng)中利用非線(xiàn)性擬合函數(shù)lsqcurvefit優(yōu)化初始種群，通過(guò)lsqcurvefit函數(shù)較強(qiáng)的局部?jī)?yōu)化能力來(lái)提高收斂速度，進(jìn)而同時(shí)具有PSO算法的全局優(yōu)化能力和lsqcurvefit擬合方法的高效局部?jī)?yōu)化能力。主要步驟如下：

步驟1隨機(jī)初始化種群中各微粒的位置和速度，并確定粒子種群規(guī)模N、學(xué)習(xí)因子c1，c2以及慣性權(quán)因子w。

步驟2用Matlab中的非線(xiàn)性擬合函數(shù)lsqcurvefit優(yōu)化初始種群位置，其調(diào)用格式為：

其中X為優(yōu)化后的初始種群，resnorm為要返回的殘差平方和，fun為擬合函數(shù)，x0為初始種群，xdata、ydata分別為實(shí)測(cè)數(shù)據(jù)。

步驟3計(jì)算每個(gè)微粒的適應(yīng)值，并將當(dāng)前各微粒的位置和適應(yīng)值存儲(chǔ)在各微粒的pbest中，將所有pbest中適應(yīng)值最優(yōu)個(gè)體的位置和適應(yīng)值存儲(chǔ)于gbest中。

步驟4用式（3）～（4）更新各微粒的速度和位移。

步驟5對(duì)于每個(gè)微粒，將其適應(yīng)值與其經(jīng)歷過(guò)的最好位置作比較，如果較好，則將其作為當(dāng)前的最好位置。

步驟6比較當(dāng)前所有pbest和gbest的值，更新gbest。

步驟7若滿(mǎn)足停止條件，搜索停止，輸出結(jié)果，否則返回步驟4繼續(xù)搜索。

3 數(shù)值實(shí)驗(yàn)與結(jié)果分析

3.1 數(shù)值實(shí)驗(yàn)

為了驗(yàn)證所提出的LPSO算法的優(yōu)化性能和優(yōu)化精度，以石家莊農(nóng)業(yè)現(xiàn)代化研究所欒城實(shí)驗(yàn)站提供的土壤數(shù)據(jù)[3]為例。算法的參數(shù)設(shè)置為：學(xué)習(xí)因子c1=c2=1.2，粒子數(shù)N=200，慣性權(quán)因子w=0.45，最大迭代次數(shù)為1 000次。在計(jì)算過(guò)程當(dāng)中，參數(shù)α，n，θs，θr的取值范圍為[0，5]。

驗(yàn)證分成兩部分。第一部分是分別采用LPSO算法與PSO算法進(jìn)行求解VG方程參數(shù)的40次數(shù)值實(shí)驗(yàn)，以分析其能夠收斂的次數(shù)、收斂效率以及收斂迭代次數(shù)來(lái)比較兩者的優(yōu)缺點(diǎn)。第二部分從最優(yōu)化準(zhǔn)則函數(shù)值的角度，分別采用lsqcurvefit擬合方法、PSO算法和LPSO算法對(duì)VG方程的參數(shù)進(jìn)行估計(jì)，以比較LPSO算法與PSO算法以及l(fā)sqcurvefit擬合函數(shù)的最優(yōu)值。

3.2 結(jié)果分析

表1給出了LPSO算法與PSO算法的優(yōu)化性能結(jié)果比較。由表1可以看出，LPSO算法的收斂次數(shù)為40次，收斂概率為100%，而PSO算法的收斂次數(shù)僅有7次，收斂概率為17.5%。由此可見(jiàn)，在收斂性能上，由于LPSO算法引入lsqcurvefit擬合函數(shù)，使其收斂效率比PSO算法大大提高。從收斂時(shí)迭代次數(shù)的分布來(lái)看，LPSO算法的收斂迭代次數(shù)在區(qū)間[101，200]的次數(shù)最大達(dá)到22次，而PSO算法在區(qū)間[101，200]的收斂迭代次數(shù)僅有3次。因此，LPSO算法與PSO算法相比，不但收斂效率提高，而且收斂迭代次數(shù)也大大增加。

表2給出了lsqcurvefit擬合方法、PSO算法和LPSO算法估計(jì)VG方程參數(shù)的函數(shù)最優(yōu)值結(jié)果比較。由表2可以看出，LPSO算法的結(jié)果要比lsqcurvefit擬合[3]方法、PSO算法的結(jié)果都小一些。另外，lsqcurvefit擬合方法對(duì)參數(shù)初值的依賴(lài)很強(qiáng)，而LPSO算法沒(méi)必要給出各參數(shù)的初值。綜合可見(jiàn)，LPSO算法對(duì)參數(shù)估計(jì)的結(jié)果要優(yōu)于lsqcurvefit擬合方法和PSO算法計(jì)算的結(jié)果。

表1 VG方程的粒子群算法優(yōu)化性能比較

表2 VG方程的參數(shù)估計(jì)結(jié)果比較

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

本文將lsqcurvefit擬合方法與粒子群（PSO）算法相結(jié)合，實(shí)現(xiàn)了不同土樣的土壤水分特征曲線(xiàn)VG方程的參數(shù)估計(jì)，并與PSO算法和lsqcurvefit擬合方法相比較，結(jié)果表明了LPSO算法求解VG方程參數(shù)效果更好，避免了初值選取和陷入局部最優(yōu)解的問(wèn)題。另外該算法也可適用于其他類(lèi)似的非線(xiàn)性參數(shù)估計(jì)問(wèn)題。

[1]Van Genuchten M.A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J].Soil Sci Soc Am J，1980，44：892-898.

[2]馬英杰，虎膽·吐馬爾拜，沈冰.利用阻尼最小二乘法求解Van Genuchten方程參數(shù)[J].農(nóng)業(yè)工程學(xué)報(bào)，2005，21（8）：l79-180.

[3]彭建平，邵愛(ài)軍.用MATLAB確定土壤水分特征曲線(xiàn)參數(shù)[J].土壤，2007，39（3）：433-438.

[4]楊改強(qiáng)，霍麗娟，楊國(guó)義，等.利用MATLAB擬合van Genuchten方程參數(shù)的研究[J].土壤，2010，42（2）：268-274.

[5]王小華，賈克力，劉景輝，等.Van Genuchten模型在土壤水分特征曲線(xiàn)擬合分析中的應(yīng)用[J].干旱地區(qū)農(nóng)業(yè)研究，2009，27（2）：179-183.

[6]許小健，黃小平，張金輪.用遺傳算法優(yōu)化估計(jì)Van Genuchten方程參數(shù)[J].巖土工程技術(shù)，2008，22（2）：75-78.

[7]郭向紅，孫西歡，馬娟娟.基于混合遺傳算法估計(jì)Van Genuchten方程參數(shù)[J].水科學(xué)進(jìn)展，2009，20（5）：677-682.

[8]陳大春，馬英杰.基于隨機(jī)粒子群算法的Van Genuchten方程參數(shù)優(yōu)化求解[J].農(nóng)業(yè)工程學(xué)報(bào)，2006，22（12）：82-84.

[9]譚躍，譚冠政，葉勇，等.具有混沌局部搜索策略的雙種群遺傳算法[J].計(jì)算機(jī)應(yīng)用研究，2011，28（2）：469-471.

[10]張艷瓊.改進(jìn)的云自適應(yīng)粒子群優(yōu)化算法[J].計(jì)算機(jī)應(yīng)用研究，2010，27（9）：3250-3252.

[11]艾海男，劉利斌，張永，等.基于Lingo和單純形算法的綜合暴雨強(qiáng)度公式參數(shù)解析[J].中國(guó)給水排水，2011，27（17）：71-74.

[12]劉利斌，歐陽(yáng)艾嘉，樂(lè)光學(xué)，等.基于混合粒子群的土壤水分特征曲線(xiàn)參數(shù)優(yōu)化[J].計(jì)算機(jī)工程與應(yīng)用，2011，47（35）：218-221.