李 琳，黃 銳，彭化偉

(1.密山市富密灌區(qū)水利水?？傉具B珠山分站，黑龍江密山158300;2.大慶地區(qū)防洪工程管理處，黑龍江大慶163000)

基于免疫遺傳算法估計(jì)Van-Genuchten方程參數(shù)

李 琳1，黃 銳2，彭化偉2

(1.密山市富密灌區(qū)水利水?？傉具B珠山分站，黑龍江密山158300;2.大慶地區(qū)防洪工程管理處，黑龍江大慶163000)

土壤水分特征曲線是模擬土壤水分運(yùn)動(dòng)和溶質(zhì)運(yùn)移的一個(gè)重要參數(shù)，Van Genuchten方程是最常用的土壤水分特征曲線方程。建立了2種Van Genuchten方程參數(shù)的優(yōu)化模型:基于非線性擬合函數(shù)lsqcurvefit與免疫遺傳算法共同改進(jìn)的混合免疫遺傳算法(LIGA)和非線性擬合函數(shù)lsqcurvefit。結(jié)果表明:分別針對(duì)粉壤土、細(xì)砂土、礫石的吸濕和脫濕曲線得到了Van Genuchten方程的相關(guān)參數(shù)，兩種方法計(jì)算的θ值與實(shí)測(cè)的θ值擬合較好，相對(duì)誤差較小，LIGA的計(jì)算精度高于lsqcurvefit函數(shù)的計(jì)算精度，改進(jìn)的LIGA算法收斂速度快，計(jì)算精度高，為該領(lǐng)域研究提供了新的思路。

遺傳算法;免疫算法;非線性擬合函數(shù);土壤水分特征曲線;土壤水吸力;土壤含水量;VG模型

土壤水分特征曲線(soil water characteristic curve，SWCC)是土壤的一個(gè)基本水力參數(shù)，表示土壤水吸力和土壤含水量之間的關(guān)系，反映了土壤水的能量與數(shù)量的關(guān)系，可以間接地反映土壤中孔隙大小的分布［1-2］。土壤水分特征曲線受多種因素的影響，其數(shù)學(xué)模型比較復(fù)雜，未知參數(shù)多由經(jīng)驗(yàn)得到，且參數(shù)較多，大部分用于描述土水特征曲線的數(shù)學(xué)模型都是根據(jù)經(jīng)驗(yàn)、土體結(jié)構(gòu)特征和曲線的形狀而建立起來(lái)的。準(zhǔn)確可靠的土壤水分特征曲線模型是分析土壤水勢(shì)的高低與植物根系吸水的速率的變化規(guī)律，研究土壤水分和相應(yīng)的吸力對(duì)作物生長(zhǎng)的有效性，對(duì)于非充分灌溉具有重要的實(shí)用價(jià)值。土壤水分特征曲線的確定方法總的來(lái)說(shuō)有直接法和參數(shù)法，直接法是指通過(guò)實(shí)驗(yàn)方法直接測(cè)定土壤水分特征曲線的方法，用數(shù)學(xué)表達(dá)(經(jīng)驗(yàn)公式)來(lái)描述水分特征曲線為參數(shù)估計(jì)法。目前，在已經(jīng)建立的眾多數(shù)學(xué)模型中，國(guó)內(nèi)外普遍使用的描述土壤水分特征曲線的方程是Van-Genuchten模型(簡(jiǎn)稱VG模型)［3-4］。根據(jù)不同土壤類型的的土壤含水率、土壤吸力的數(shù)據(jù)，運(yùn)用參數(shù)估計(jì)法求解VG模型的較多，如廖林仙等［5］運(yùn)用智能算法推求了VG方程參數(shù);馬英杰等［6］運(yùn)用阻尼最小二乘法，擬合了VG方程參數(shù);彭建平等［7］對(duì)VG模型中的參數(shù)利用Matlab軟件中的非線性擬合函數(shù)來(lái)確定，模型計(jì)算取得了較好的效果;郭向紅等［8］利用混合遺傳算法對(duì)VG方程進(jìn)行了求解;還有一些其他方法在推求VG方程的參數(shù)中取得了良好的效果。

混合免疫遺傳算法是非線性擬合函數(shù)(lsqcurvefit)和免疫遺傳優(yōu)化算法(immune genetic algorithm，IGA)的耦合模型，具有免疫遺傳算法全局信息交換和局部深度搜索相結(jié)合的特點(diǎn)，利用lsqcurvefit函數(shù)的非線性擬合能力來(lái)捕獲變量間的非線性關(guān)系，進(jìn)一步提高模型的局部搜索能力，加快向最優(yōu)解的收斂速度。模型優(yōu)化結(jié)果證明，該模型具有較強(qiáng)的非線性逼近能力和較強(qiáng)的魯棒性、較快的收斂速度、較高的搜索精度。本文將混合免疫遺傳算法引用于VG方程中，為方程參數(shù)的確定提供參考依據(jù)。

1 Van Genuchten方程

Van Genuchten公式描述壓力水頭與含水量的關(guān)系，即土壤水分能量與數(shù)量的關(guān)系，方程形式為:

式中:θ為體積含水率，cm3/cm3;θr為殘留含水率，cm3/cm3;θS為飽和含水率，cm3/cm3;h 為負(fù)壓，cmH2O;α，n，m 分別表示土壤水分特征曲線形狀的參數(shù)。式(1)中含有4個(gè)參數(shù)α，n(其中n中含有m)，θS，θr，根據(jù)實(shí)測(cè)的土壤含水率θ和土壤水吸力h資料，采用最小二乘法擬合結(jié)果的殘差平方和最小為目標(biāo)，即式(2)，通過(guò)優(yōu)化得到這4個(gè)參數(shù):

式中:θi為實(shí)測(cè)含水率，cm3/cm3;hi為實(shí)測(cè)土壤含水率θi對(duì)應(yīng)的土壤水吸力，cm;θ( hi，X)為根據(jù)式(1)計(jì)算的土壤含水率，cm3/cm3;X 為待求參數(shù)向量(α，n，θS，θr);N 為實(shí)測(cè)數(shù)據(jù)個(gè)數(shù)。

2 算法的構(gòu)建

基于非線性擬合函數(shù)改進(jìn)的混合免疫遺傳算法(Lsqcurvefit Immune genetic algorithm，LIGA)將遺傳算法［9］(genetic algorithm，GA)、免疫算法 (Immune Algorithm，IA)、非線性擬合函數(shù)lsqcurvefit結(jié)合構(gòu)造的一種新的混合免疫遺傳算法。首先，在Matlab優(yōu)化工具箱中，lsqcurvefit函數(shù)采用高斯—牛頓(Gauss-Newton)方法或者麥夸特(Levenberg-Marquardt)方法，用于求解最小二乘非線性數(shù)據(jù)擬合問(wèn)題?；旌厦庖哌z傳算法以lsqcurvefit函數(shù)的非線性求解作為種群的初始染色體的最優(yōu)值，避免了算法開始階段搜索的隨機(jī)性與盲目性，使得混合免疫遺傳算法具有很強(qiáng)的非線性擬合能力與非線性逼近能力，群體具有的免疫機(jī)制對(duì)群體進(jìn)行控制和調(diào)節(jié)，把目標(biāo)函數(shù)和制約條件作為群體的抗原，保證生成的染色體直接與問(wèn)題相關(guān)聯(lián)，收斂方向得以控制，對(duì)與抗原親和力高的個(gè)體進(jìn)行記憶，并從中提取免疫疫苗對(duì)抗體群進(jìn)行注射，保證抗體群在更新過(guò)程中的多樣性，促進(jìn)求解的速度，提高算法的效率，并防止群體的退化，在很大程度上抑制混合免疫遺傳算法的未成熟收斂［10］。結(jié)合免疫算法、遺傳算法、lsqcurvefit函數(shù)的原理，下面給出具體實(shí)現(xiàn)步驟:

Step1:抗原輸入。輸入目標(biāo)函數(shù)式(2)和各種約束作為混合免疫遺傳算法的抗原。

Step2:產(chǎn)生初始群體。生成初始抗體群N、促進(jìn)記憶細(xì)胞庫(kù)N1、檢測(cè)記憶細(xì)胞庫(kù)N2。對(duì)初次應(yīng)答，初始抗體N、N1、N2全部隨機(jī)產(chǎn)生，而對(duì)再次應(yīng)答，則抗體群N借助免疫系統(tǒng)的記憶機(jī)制，部分初始抗體由免疫記憶單元獲取，即包含最優(yōu)抗體促進(jìn)記憶細(xì)胞庫(kù)，其余抗體隨機(jī)產(chǎn)生。

Step3:計(jì)算抗體適應(yīng)度。群體的抗原和抗體V之間的親和力:axv=optv，其中optv表示抗原和抗體之間的匹配程度，本文用抗原和抗體之間的適應(yīng)值函數(shù)fittness(v)來(lái)表示，即axv=fittness(v)。N個(gè)抗體的信息熵為:

式中:Hj(N)為N個(gè)抗體第j位的信息熵，Pij為N個(gè)抗體中的第j位為數(shù)值ki的概率?？贵wv和抗體w的親和力為:

抗體v的濃度cv為:

式中:Tac1為一個(gè)預(yù)先確定的親和力閾值;

Step4:免疫記憶細(xì)胞庫(kù)更新。根據(jù)一定比例，從抗體群N中選出親和力高的抗體，用它們替換促進(jìn)記憶細(xì)胞庫(kù)N1中親和力低的抗體;將抗體群N中親和力低的抗體選入檢測(cè)記憶細(xì)胞庫(kù)N2，用它們替換檢測(cè)記憶細(xì)胞庫(kù)N2中親和力高的抗體。

Step5:抗體群的促進(jìn)與抑制。當(dāng)記憶群體N1中抗體濃度的最大值cmax低于抗體濃度閾值th時(shí)，記憶群體中的抗體處于多樣化，為未飽和狀態(tài);否則抗體趣于一致化，為飽和狀態(tài)，濃度大的抗體將被淘汰，隨機(jī)產(chǎn)生的新個(gè)體代替被淘汰的個(gè)體。用檢測(cè)記憶群體N2去檢測(cè)N中是否含有已搜索過(guò)的抗體(即N2中記錄的抗體)，如果有就用隨機(jī)抗體取代它。

Step6:將促進(jìn)記憶群體N1與檢測(cè)后的抗體群N相結(jié)合生成新的抗體種群，即父代群體y。

Step7:利用Matlab優(yōu)化工具箱中的非線性擬合函數(shù)lsqcurvefit對(duì)目標(biāo)函數(shù)進(jìn)行優(yōu)化，將優(yōu)化的非線性結(jié)果作為初始個(gè)體最優(yōu)值置于父代群體y中，提高各染色體對(duì)初始最優(yōu)解的搜索能力，非線性擬合的數(shù)學(xué)模型［11］［12］:

根據(jù)輸入數(shù)據(jù)xdata和得到的輸出數(shù)據(jù)ydata找到與目標(biāo)函數(shù)F(x，xdata)最佳的擬合參數(shù)向量x，以x作為群體中每條染色體的初始最優(yōu)個(gè)體值。

Step8:對(duì)于父代群體y進(jìn)行選擇操作、雜交操作和變異操作［9-13］。

Step9:終止條件是否滿足?如果滿足，結(jié)束迭代，否則轉(zhuǎn)向Step3。

3 實(shí)例分析

以文獻(xiàn)［14］中的粉壤土、細(xì)沙土、礫石3種介質(zhì)的吸濕和脫濕曲線的實(shí)測(cè)資料為例，說(shuō)明混合免疫遺傳算法在Van Genuchten方程參數(shù)擬合中的應(yīng)用效果。表1為分別采用基于非線性擬合函數(shù)lsqcurvefit改進(jìn)的LIGA和非線性擬合函數(shù)lsqcurvefit兩種方法得到的不同土壤類型的Van Genuchten方程參數(shù)結(jié)果。采用MATLAB7.0編程處理，選取選定父代初始群體規(guī)模為n=400，抗體濃度閾值th=0.85，交叉概率pc=0.8，變異概率pm=0.8，優(yōu)秀個(gè)體數(shù)目選定20個(gè)，Van-Genuchten方程的估計(jì)參數(shù)取值范圍，α取值為0～1、n取值為1～10、θS取值為 0.2～0.7、θr取值為 0～0.2;通過(guò)MATLAB程序設(shè)計(jì)與其工具箱命令的條用實(shí)現(xiàn)非線性曲線擬合函數(shù)lsqcurvefit的計(jì)算。

表1 土壤水分特征曲線Van Genuchten方程參數(shù)計(jì)算結(jié)果

從表1中結(jié)果可以看出，LIGA與lsqcurvefit函數(shù)優(yōu)化的參數(shù)估計(jì)結(jié)果比較接近，進(jìn)而比較2種方法的求解擬合值與實(shí)測(cè)值之間的殘差平方和范數(shù)可知，對(duì)于目標(biāo)函數(shù)值優(yōu)化計(jì)算，LIGA優(yōu)化結(jié)果計(jì)算結(jié)果θ值與實(shí)測(cè)的θ的殘差平方和范數(shù)小于lsqcurvefit函數(shù)的計(jì)算，可見(jiàn)改進(jìn)的LIGA優(yōu)化的參數(shù)結(jié)果優(yōu)于lsqcurvefit函數(shù)優(yōu)化的參數(shù)結(jié)果，并且這兩種方法的計(jì)算均達(dá)到了很高的精度。

4 結(jié)論

Van Genuchten方程是最常用的土壤水分特征曲線方程，方程中的參數(shù)擬合具有非線性特點(diǎn)，本文將免疫算法、遺傳算法和lsqcurvefit函數(shù)三者相結(jié)合構(gòu)建新的混合免疫遺傳算法(LIGA)，并且嘗試性地將改進(jìn)后的LIGA用于Van Genuchten方程參數(shù)擬合中，并進(jìn)行了數(shù)值試驗(yàn)，結(jié)果表明:改進(jìn)的LIGA的計(jì)算精度比lsqcurvefit模型高，LIGA模型用于Van Genuchten方程參數(shù)求解問(wèn)題中是可行的，既拓展了LIGA的應(yīng)用領(lǐng)域，又為解決Van Genuchten方程參數(shù)擬合問(wèn)題提供了新的思路和方法。

［1］趙愛(ài)輝，黃明斌，史竹葉.兩種土壤水分特征曲線間接推求方法對(duì)黃土的適應(yīng)性評(píng)價(jià)［J］，農(nóng)業(yè)工程學(xué)報(bào)，2008，24(9):11.

［2］盧小慧，靳孟貴，汪丙國(guó).欒城農(nóng)業(yè)生態(tài)系統(tǒng)試驗(yàn)站土壤水分特征曲線分析［J］，中國(guó)農(nóng)村水利水電，2006(12):30.

［3］李春友，任理，李保國(guó).利用優(yōu)化方法求解VanGenuchten方程參數(shù)［J］，水科學(xué)進(jìn)展，2001，12(4):473-475.

［4］呂殿青，邵明安.非飽和土壤水力參數(shù)的模型及確定方法，應(yīng)用生態(tài)學(xué)報(bào)［J］，應(yīng)用生態(tài)學(xué)報(bào)，2004，15(1):63-165.

［5］廖林仙，邵孝侯，徐俊增.基于智能算法推求VanGenuchten方程的參數(shù)［J］，水利學(xué)報(bào)，2007(1):696-697.

［6］馬英杰，虎膽·吐馬爾拜，沈 冰.利用阻尼最小二乘法求解VanGenuchten方程參數(shù)［J］，農(nóng)業(yè)工程學(xué)報(bào)，2005，21(8):179-180.

［7］王小華，賈克力，劉景輝，李立軍.VanGenuchten模型在土壤水分特征曲線擬合分析中的應(yīng)用［J］，干旱地區(qū)農(nóng)業(yè)研究，2009，27(2):179-180.

［8］郭向紅，孫西歡，馬娟娟.基于混合遺傳算法估計(jì)van Genuchten方程參數(shù)［J］，水科學(xué)進(jìn)展，2009，20(5):677-680.

［9］付強(qiáng)，金菊良，梁川.基于實(shí)碼加速遺傳算法的投影尋蹤分類模型在水稻灌溉制度優(yōu)化中的應(yīng)用［J］，水利學(xué)報(bào)，2002(10):39-42.

［10］葛紅，毛宗源.免疫算法的改進(jìn)［J］，計(jì)算機(jī)工程與應(yīng)用，2002，14(47):47-48.

［11］彭建平，邵愛(ài)軍.用MatLab確定土壤水分特征曲線參數(shù)［J］，土壤，2007，39(3):433-436.

［12］唐家德.基于MATLAB的非線性曲線擬合［J］，計(jì)算機(jī)與現(xiàn)代化，2008(6):15-18.

［13］金菊良，楊曉華，丁晶.標(biāo)準(zhǔn)遺傳算法的改進(jìn)方案—加速遺傳算法.系統(tǒng)工程理論與實(shí)踐［J］.2001(4):9-10.

［14］王金生，楊志峰，陳家軍，等.包氣帶土壤水分滯留特性研究［J］. 水利學(xué)報(bào)，2000(2):1-6.

S15

A

1007-7596(2012)02-0046-03

2011-12-19

李琳(1977-)，女，黑龍江密山人，工程師;董銳(1978-)女，黑龍江密山人，工程師;彭化偉(1976-)，男，黑龍江呼瑪人，工程師。。