孫曉敏，孫秋紅

(西北農林科技大學 黃土高原與土壤侵蝕旱地農業(yè)國家重點實驗室， 陜西 楊凌 712100)

基于改進的基本粒子群算法的Van Genuchten方程參數優(yōu)化

孫曉敏，孫秋紅

(西北農林科技大學 黃土高原與土壤侵蝕旱地農業(yè)國家重點實驗室， 陜西 楊凌 712100)

針對目前土壤水分特征曲線Van Genuchten方程參數優(yōu)化的不足，引入動態(tài)調整慣性權重對基本粒子群算法進行改進，使慣性權重隨著迭代次數以及不同粒子與最優(yōu)粒子之間的距離大小而變化，并將其運用到Van Genuchten方程參數識別，最后進行了模型驗證和誤差分析。結果表明改進慣性權重的粒子群算法計算精度高，適用性強。從3種準則函數的優(yōu)化結果可以看出，加權耦合的絕對與相對標準差最小準則在參數優(yōu)化理論上值得進一步研究，在Van Genuchten方程參數擬合問題中值得推薦，從而為Van Genuchten方程參數的優(yōu)化求解提供了一條新途徑。

粒子群算法；慣性權重；Van Genuchten方程；準則函數；參數估計

在農田水利工程、節(jié)水灌溉和水土保持工程等設計中，土壤水分常數及水動力學參數極其重要，其中一個主要的參數獲取手段就是進行土壤持水特性的測定，即土壤水分特征曲線[1]。因此研究土壤持水曲線始終是土壤物理學界和農業(yè)節(jié)水領域的重中之重。土壤水分特征曲線反映了水分和化學物質在非飽和土壤中運移的特性，該曲線定義了土壤水分含量與壓力水頭之間的函數關系，該函數的可靠性可以直接影響土壤水分運動模型的預測結果。

目前存在許多經驗公式用于描述土壤持水曲線，較為常用的有Brooks-Corey模型[2]、Gardner模型[3-4]、Van Genuchten(簡稱VG)模型[5]和Gardner-Russo模型[6]。徐紹輝等[7]分析研究了這4個模型的適應性，他認為VG模型對較黏質地的土壤和粗質地土壤的擬合效果均較好。夏衛(wèi)生等[8]對國內外土壤水動力學參數研究結果做出分析，得出VG模型擬合效果較好，此外還得出其可以與土壤的機械組成、容重等聯(lián)系起來，從土壤本身特性上找到其含義。因此，在描述土壤持水曲線的全部模型中，VG模型由于線型與實測數據曲線擬合程度好而被廣泛應用。

通常VG模型參數的求解方法主要有：王金生等[9-10]采用最小二乘法和非線性單純形法相結合的方法，研究了土壤持水曲線的滯留特征；徐紹輝等[11]用最小二乘法結合Picard迭代法擬合出了砂質黏壤土持水曲線的VG方程參數；李春友等[12]借助單純形調優(yōu)方法研究了VG方程參數的擬合問題；馬英杰等[13]采用非線性阻尼最小二乘法擬合了VG方程的參數，以上方法對方程參數的求解均提供了較好的方法，但全局收斂性相對較弱。近年來，陳大春等[14]運用隨機粒子群算法求解VG方程的參數，該方法具有較強的全局搜索能力和較高的求解精度，但迭代次數較多，計算量較大。因此，尋求一種全局搜索能力較強，計算精度強的尋優(yōu)方法是VG模型參數求解亟待解決的科學問題，本文采用動態(tài)調整慣性權重的粒子群算法，并將優(yōu)化準則函數的選取和全局尋優(yōu)策略兩者相結合以期達到收斂性強、計算精度高特點。并通過實測數據驗證上述三種準則函數的優(yōu)化效果，為該問題的求解提供有益的參考。

1 材料和方法

1.1 Van Genuchten模型

Van Genuchten模型[4]方程的基本形式如下：

(1)

式中：θ為土壤體積含水率，cm3/cm3；h為土壤基質勢，cm；θs為土壤飽和含水率，cm3/cm3；θr為殘余含水率，cm3/cm3；α為進氣值的倒數，cm-1；n為孔徑分布指數；m=1-1/n。α參數取值為0～1 cm-1，通常為0.005 cm-1～0.05 cm-1；θs取值一般為：黏土0.5 cm3/cm3～0.6 cm3/cm3、壤土0.4 cm3/cm3～0.5 cm3/cm3、砂土0.3 cm3/cm3～0.4 cm3/cm3；θr取值為0～0.3 cm3/cm3；n取值為1～10，通常為1.1～3.5[5]。

Van Genuchten模型中共有4個相互獨立的參數，因此，土壤水分運動曲線VG模型的求解，實際上是一個4參數的非線性超越函數的參數優(yōu)化問題。

1.2 準則函數

在VG模型參數尋優(yōu)中，優(yōu)化準則的合理選取是關鍵之一。由于在實際尋優(yōu)過程中模擬值不可能與實測值完全吻合，當模擬土壤含水率點據分布于實際土壤含水率曲線的兩側，必然會存在一定的離差Δ=θ模擬-θ實測，由最小二乘原理，首先可按絕對標準差最小準則(準則1；公式(2))擬合：

(2)

式中：θi(hi)為實測土壤含水率，cm3/cm3；θp(hi，X)為模型模擬的土壤含水率，cm3/cm3；hi為實測土壤基質勢，cm；X為待優(yōu)選的參數(用行向量表示)；n為實測數據組數。

其次，可按相對標準差最小準則(準則2；公式(3))擬合：

(3)

在更一般的情況下，應該有：

(4)

顯然，λ∈[0，-2]當λ=0時，即為絕對標準差最小的擬合準則，其本質是取各經驗數據點據均為1的權重，當為λ=-2時，即為相對標準差最小的擬合準則，其本質是取各數據點據隨[θi(hi)]-2變化的權重，由此可以推斷，當取中間值，即λ=-1時的擬合準則為：

(5)

上述方程其本質是取各經驗數據點據隨1/[θi(hi)]變化的權重，這實際上是一種考慮了相對標準差最小與加權耦合的絕對標準差最小的擬合準則(準則3；公式(5))。

1.3 動態(tài)調整慣性權重的粒子群算法

1.3.1 基本粒子群算法

粒子群優(yōu)化(Particle Swarm Optimizer，PSO)算法作為一種新的基于群體智能的演化優(yōu)化方法，它第一次是由美國普渡大學的Kennedy和Eberhart博士于1995年提出的[15-18]。PSO的基本概念來源于對鳥群捕食行為的研究，因此也可稱作鳥群算法，是集群智能的代表性方法之一。該算法優(yōu)點在于可并行處理以及魯棒性好，既簡單易實現又具有深刻的智能背景。因此，PSO算法在提出后就引起了廣大學者的廣泛關注，作為一種新的并行優(yōu)化進化算法，PSO算法可以解決大量非線性、不可微、非連續(xù)和多峰的優(yōu)化問題，且已廣泛應用于科學和工程領域，如模糊系統(tǒng)控制、神經網絡訓練、函數優(yōu)化、模式分類等[17]。

假設在D維的目標搜索空間，存在m個粒子可以組成一個群落，其中第i個粒子表示為一個D維的向量xi=(xi1，xi2，xi2，…，xiD),i=1,2,3,…,m,即第i個粒子在D維搜索空間中的位置。將xi帶入目標函數即可計算出適應度值，進而以適應度值的大小判斷xi的優(yōu)劣。第i個粒子的“飛翔”速度也可作為一個D維的向量，記作vi=(vi1，vi2，vi3，…，viD)，將第i個粒子迄今為止搜索到的最優(yōu)位置記為pi=(pi1，pi2，pi3，…，piD)，在整個粒子群中，所有粒子所搜索到的最優(yōu)位置記為pgi=(pg1，pg2，pg3，…，pgD)。對于第k次迭代，每個粒子的進化采用公式(6)進行操作：

(6)

(7)

式中i=1，2，…，m是粒子數；d=1，2，…，D是向量維數；k=1，2，…，n為迭代次數；學習因子c1和c2為非負常數，分別表示調節(jié)向全局最好粒子、個體最好粒子方向所飛行的步長值；r1與r2是介于(0，1)之間的任意數，ω是慣性權重。根據具體問題確定迭代中止條件，一般選擇最大迭代次數或者粒子群所搜索到的最優(yōu)位置滿足預定值來作為終止條件?；綪SO算法中需要調整的參數不多，而且采用實數編碼可使其操作簡便, 因此使用起來較為方便，但其具有容易陷入局部極值點，搜索精度不高的缺點。

1.3.2 動態(tài)調整慣性權重的PSO算法

當公式(6)中的慣性權重ω比較大時，vid靠近全局最優(yōu)點pg和靠近歷史最優(yōu)點pi的機會較小，則vid會偏離當前最優(yōu)點，從而有可能探索更好的最優(yōu)點，這在早期中是非常有利的，但在晚期是不利的，這是由于晚期不利于全局最優(yōu)點的收斂。然而當ω較小時，正好與上面情況相反。所以，ω值的大小對于收斂與探索的影響是一對矛盾。故而在應用中先將ω的初始值作為最大值，并使其隨著迭代次數的增加而線性遞減至最小值，以達到所期望的優(yōu)化目標，權重函數ω的值用下式來確定：

(8)

式中：ωmax、ωmin分別表示ω的最大值和最小值；niter、nitermax分別表示當前迭代次數和最大迭代次數。但以上權重函數的設置也存在不利的一面，假設在迭代初期就已找到全局最優(yōu)點，可能會因其權重過大而跳出這個最優(yōu)點，因而不在其周圍探索，導致算法收斂速度降低。若想提高收斂速度就必須賦予靠近最優(yōu)點的部分粒子較小的權重，使粒子在靠近pg領域范圍進行微小探索，而不是承擔更大范圍的探索。賦予遠離最優(yōu)點的其它粒子以較大權重，令其承擔更大范圍的探索任務，從而進一步去探索更優(yōu)的點。因此，慣性權重ω的值不但隨著迭代次數而變化，而且對不同粒子是根據其與最優(yōu)粒子之間的距離大小而變化。假設Iig表示第i粒子與全局最優(yōu)粒子Pg的距離，則慣性權重可采用公式(9)調整：

(9)

式中：Imax、Imin分別為粒子距全局最優(yōu)值的最大距離和最小距離。

1.3.3 改進PSO算法在VG方程參數估計中的實現

對VG方程參數可理解為在土壤水力參數范圍內隨機生成的一組水深序列，以此作為初始粒子種群，再根據不同的優(yōu)化準則設定目標函數，計算每個粒子的適應度，然后按照式(6)、式(9)進行迭代計算，直至達到最大迭代次數，具體的迭代步驟為：

步驟1 確定需要優(yōu)化變量的取值范圍，θs∈[0.4,0.6]，θr∈[0,3]，α∈[0,1]，n∈[0,10]。

步驟2 設定PSO的計算參數：粒子群體大小、學習因子c1和c2、最大迭代次數nitermax、慣性權重的最大值ωmax和最小值ωmin，粒子最大距離Imax和最小距離Imin。

步驟3 在優(yōu)化變量θs，θr，α，n的取值范圍內，采用蒙特卡洛隨機方法產生初始群體粒子個體。

步驟4 根據不同的優(yōu)化準則，計算出群體中每個粒子的適應度值。

步驟5 依據每個粒子計算出來的適應值大小，判斷該粒子迄今為止所經歷的最好位置(即準則函數最小的土壤水力學參數)以及整個粒子群迄今為止搜索到的最好位置(全局最優(yōu)解)。

步驟6 計算出每個粒子相對全局最優(yōu)解的距離并根據式(9)計算慣性權重。

步驟7 根據式(6)，更新所有粒子的速度和位置，若新速度超出最大速度限制[vmin,vmax]，則設新速度為最大速度。

步驟8 結束條件的判斷。若達到最大迭代次數，則停止計算，輸出結果；否則，轉至步驟4繼續(xù)計算，直至滿足結束條件為止。

為驗證本文提出算法的全局收斂性及不同優(yōu)化準則的模擬效果，特以文獻[3]中所給的實測資料作為本文方法的計算原始數據，對以上過程編程計算。取計算參數為：粒子群體大小20，最大迭代次數1000，初始慣性權重0.8，最終慣性權重0.2,學習因子c1=c2=2，最大粒子距離Imax=1最小粒子距離Imax=0.2，經過MATLAB 7.1對以上算法進行編程計算[19]，采用改進PSO算法連續(xù)進行10次獨立計算，從中選取最優(yōu)結果作為最終優(yōu)化結果。

由以上分析可知，VG模型參數尋優(yōu)中要想達到計算速度快、精度高的雙重要求就必須將準則函數與算法兩者相結合。為凸顯改進PSO算法的優(yōu)點及準則函數的適應性，特選取目前比較典型的幾種方法，包括文獻[10]的非線性單純形法、文獻[13]的非線性阻尼最小二乘法和文獻[14]的隨機粒子群算法(上述求解方法的準則函數均采用準則1)與本文改進的PSO算法和準則1(改進的PSO算法1)、改進的PSO算法和準則2(改進的PSO算法2)和改進的PSO算法和準則3(改進的PSO算法3)共6種方法進行比較研究，以期檢驗改進PSO算法的優(yōu)越性及準則函數的適用性。

1.4 試驗數據

提供試驗土壤為陜西楊凌地區(qū)的黏壤土，自然風干后經過1 mm孔徑的土篩選后利用沉降法進行顆粒分析(見表1)。土壤水分特征曲線采用高速離心機測定(見表2)。

表1 試驗土壤顆粒級配

表2 土壤水分特性曲線實測數據

2 結果與討論

2.1 不同優(yōu)化準則的比較

由表3可知，不管采用哪種優(yōu)化準則，即本文述及的3種優(yōu)化準則，改進的PSO算法明顯優(yōu)于非線性單純性法和阻尼最小二乘法，同時也優(yōu)于隨機粒子群法的優(yōu)化結果。由此可知，本文提出的改進慣性權重模型的PSO算法能更好的尋求全局最優(yōu)解，是一種很好的全局搜索策略。此外，還可以看出，3種優(yōu)化準則其擬合精度均較高，這說明在現行VG方程參數優(yōu)化中采用絕對標準誤差最小準則基本上是科學、合理、可行的。但仔細分析表3數據可以看出，采用準則3雖然計算誤差不能保證絕對標準誤差和相對標準誤差同時達到最小，但其絕對標準誤差小于準則1的計算結果，且相對標準誤差小于準則2擬合的結果。

從表4中最大相對誤差絕對值和平均相對誤差絕對值的比較可知，從模型求參方法來看，改進慣性權重粒子群算法模擬精度均高于文獻[10,13]的方法。由平均相對誤差來看，改進PSO與文獻[14]相差不大，但以相對誤差的絕對值來看，改進PSO方法優(yōu)于隨機粒子群算法。從優(yōu)化準則來看，準則3擬合結果中平均誤差為5.675%，小于準則1優(yōu)化結果(平均相對誤差為7.899)，略大于準則2的優(yōu)化結果(平均相對誤差為5.338)，而最大相對誤差的絕對值明顯小于準則1和準則2。

通過以上分析可知，加權耦合的絕對與相對標準差最小準則，解決了僅以絕對標準差最小準則或僅以相對標準差最小準則進行參數擬合所存在的某些不足，故加權耦合的絕對與相對標準差最小準則是一種在土壤VG方程參數擬合問題中值得推薦的優(yōu)化準則。

表3 不同方法確定的粉壤土Van Genuehten方程參數

表4 不同方法模擬θ-h數據與實測數據對比

2.2 VG模型參數的優(yōu)化

以上分析可知，加權耦合優(yōu)化準則優(yōu)于絕對標準差最小準則和相對標準差最小準則，由此可將加權耦合優(yōu)化準則運用于本文試驗數據的優(yōu)化擬合中，并采用改進慣性權重的粒子群優(yōu)化算法進行參數的尋優(yōu)。

由表5和表6可知，最大相對誤差的絕對值僅為3.135，且分布在土壤水吸力1000 cm以上，該土壤水吸力區(qū)間對應的含水率變化區(qū)間很小僅為0.2228～0.1636之間，在實際中應用相對也較少，而在土壤水吸力小于1000 cm以下時，擬合精度相對較高。總體來說，采用加權耦合優(yōu)化準則的改進PSO算法比非線性單純形法和阻尼最小二乘法及隨機粒子群算法要優(yōu)越，且不受初值的影響。

表5 根據加權耦合優(yōu)化準則確定的黏壤土Van Genuehten方程參數

表6 模擬值與實測數據對比分析

3 結 論

(1) 本文提出了一種全局優(yōu)化策略—改進慣性權重粒子群優(yōu)化算法，通過實例計算表明，該算法全局搜索能力強，優(yōu)于現有幾種優(yōu)化方法。

(2) 本文對Van Genuchten方程參數尋優(yōu)中準則函數進行比較研究表明，加權耦合絕對標準差和相對標準差最小準則雖然計算誤差不能滿足絕對標準誤差和相對標準誤差同時達到最小，但可以解決僅以絕對標準差最小準則或僅以相對標準差最小準則進行參數擬合所存在的某些不足，可作為Van Genuchten方程參數擬合中較優(yōu)的優(yōu)化準則。

(3) 加權耦合的絕對與相對標準差最小準則是一種在回歸理論上值得進一步研究的問題，同時可將其運用到其它非線性問題中，驗證其適宜性，可望成為非線性參數擬合問題中值得推薦的準則。

[1] 雷志棟，楊詩秀，謝森傳．土壤水動力學[M].北京：清華大學出版社，1988：18-24.

[2] Miay P C D. Estimation of the brooks corey parameters from water retention data[J]. Water Resource Research，1987，23(6)：1085-1089.

[3] Gardner W R, Hill D, Benyamini Y. Post irrigation movement of soil water redistribution[J]. Water Resource Research，1970，6(3)：851-861.

[4] Gardner W R, Hillel D, Benyamini Y. Post irrigation movement of soil waterⅡ Simultaneous redistribution and evaporation[J]. Water Resolve Research，1970，6(4)：1148-1153.

[5] Genuchten M T V. A closed form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Sci．See．Am．J, 1980，44：892-898.

[6] Russo D. Determining soil hydraulic properties by parameter：On the selection of model for the hydraulic properties estimation[J]. Water Resource Research，1988，24(3)：453-459.

[7] 徐紹輝，張佳寶，劉建立，等．表征土壤水分持留曲線的幾種模型的適應性研究[J]．土壤學報，2000，37(3)：271-274．

[8] 夏衛(wèi)生，雷廷武，潘英華，等.土壤水動力學參數研究與評價[J].灌溉排水，2002，21(1):72-75.

[9] Shao M A, Horton Robert. Integral method for estimating soil hydraulic properties[J] .Soil Sci Soc Am J，1998，62(3)：585-592.

[10] 王金生，楊志峰，陳家軍，等.包氣帶土壤水分滯留特征研究[J].水利學報,2000,31(2):1-6.

[11] 徐紹輝，張佳寶．求土壤水力特征的一種迭代方法[J].土壤學報，2000，37(3) :271-274.

[12] 李春友，任 理，李保國.利用優(yōu)化方法求解Van Genuchen方程參數[J].水科學進展，2001，12(4)：473-478.

[13] 馬英杰，虎膽.吐馬爾拜，沈 冰.利用阻尼最小二乘法求解Van Genuchten方程參數[J].農業(yè)工程學報，2005，21(8)：179-180.

[14] 陳大春，馬英杰.基于隨機粒子群算法的Van Genuchten方程參數優(yōu)化求解[J].農業(yè)工程學報，2006，22(12)：82-84.

[15] Parker J C, Kool J B, Genuchten M T V. Determining soil hydraulic properties from one2step outflow experiments parameter estimation: I. Theory and numerical studies[J]. Soil Sci Soc Am J, 1985,49(6)：1348-1354.

[16] Eberhart R, Kennedy J. A new optimizer using particle swarm theory[C]//Proc 6th Int Symposium on Micro Machine and Human Science, Nagoya, 1995:39-43.

[17] Kennedy J, Eberhart R. Particle swarm optimization[C]//Proc of IEEE Int Conf on Neural Networks Perth, 1995:1942-1948.

[18] Shi Y, Eberhart R C. Empirical study of particle swarm optimization[C]//Proceedings of the IEEE World Congress on Evolutionary Computation, IEEE Press, 1999:1945-1950.

[19] 龔 純，王正林.精通MATLAB最優(yōu)化計算[M].北京：電子工業(yè)出版社，2009：258-268.

ParameterCalibrationoftheVanGenuchtenEquationbasedonModifiedParticleSwarmOptimization

SUN Xiaomin, SUN Qiuhong

(InstituteofSoilandWaterConservation,NorthwestA&FUniversity,Yangling,Shaanxi712100,China)

Focusing on the issue of parameter calibration of the Van Genuchten equation, an innovated dynamic regulation of inertia weight was introduced into the normal particle swarm optimization (PSO) to ensure inertia weight change both with the number of iterations and with the distance between various particles to optimal particles. Furthermore the new method was applied in the comparisons of criterion functions in parameters identification of the Van Genuchten equation. Finally the new method was verified in model validation and error analysis. The results showed that the new PSO method with improved inertia weights displayed higher accuracy and wider applicability. Besides, optimization results based on the three criterion functions revealed that the criterion of minimum absolute and relative standard deviations in weighted coupling proved to be an optimization criterion deserving further study and prior selection in parameter estimation for the Van Genuchten equation in soil science. The new method proposed an innovated approach for the calibration of Van Genuchten equation parameters.

particleswarmoptimization;inertiaweight;VanGenuchtenequation;criterionfunction;parameterestimation

10.3969/j.issn.1672-1144.2017.05.015

2017-05-05

2017-06-25

國家自然科學基金資助項目(51579214；41001159)；國家973計劃課題(2015CB040441)；黃土高原土壤侵蝕與旱地農業(yè)國家重點實驗室主任基金項目(A314021402-1619)

孫曉敏(1986—)，女，陜西寶雞人，主要從事水土保持與生態(tài)修復方面的研究工作。 E-mail：sunxiaomin428@126.com

S152.7

A

1672—1144(2017)05—0088—06