基于改进的基本粒子群算法的Van Genuchten方程参数优化

孙晓敏，孙秋红

(西北农林科技大学 黄土高原与土壤侵蚀旱地农业国家重点实验室， 陕西 杨凌 712100)

基于改进的基本粒子群算法的Van Genuchten方程参数优化

孙晓敏，孙秋红

(西北农林科技大学 黄土高原与土壤侵蚀旱地农业国家重点实验室， 陕西 杨凌 712100)

针对目前土壤水分特征曲线Van Genuchten方程参数优化的不足，引入动态调整惯性权重对基本粒子群算法进行改进，使惯性权重随着迭代次数以及不同粒子与最优粒子之间的距离大小而变化，并将其运用到Van Genuchten方程参数识别，最后进行了模型验证和误差分析。结果表明改进惯性权重的粒子群算法计算精度高，适用性强。从3种准则函数的优化结果可以看出，加权耦合的绝对与相对标准差最小准则在参数优化理论上值得进一步研究，在Van Genuchten方程参数拟合问题中值得推荐，从而为Van Genuchten方程参数的优化求解提供了一条新途径。

粒子群算法；惯性权重；Van Genuchten方程；准则函数；参数估计

在农田水利工程、节水灌溉和水土保持工程等设计中，土壤水分常数及水动力学参数极其重要，其中一个主要的参数获取手段就是进行土壤持水特性的测定，即土壤水分特征曲线[1]。因此研究土壤持水曲线始终是土壤物理学界和农业节水领域的重中之重。土壤水分特征曲线反映了水分和化学物质在非饱和土壤中运移的特性，该曲线定义了土壤水分含量与压力水头之间的函数关系，该函数的可靠性可以直接影响土壤水分运动模型的预测结果。

目前存在许多经验公式用于描述土壤持水曲线，较为常用的有Brooks-Corey模型[2]、Gardner模型[3-4]、Van Genuchten(简称VG)模型[5]和Gardner-Russo模型[6]。徐绍辉等[7]分析研究了这4个模型的适应性，他认为VG模型对较黏质地的土壤和粗质地土壤的拟合效果均较好。夏卫生等[8]对国内外土壤水动力学参数研究结果做出分析，得出VG模型拟合效果较好，此外还得出其可以与土壤的机械组成、容重等联系起来，从土壤本身特性上找到其含义。因此，在描述土壤持水曲线的全部模型中，VG模型由于线型与实测数据曲线拟合程度好而被广泛应用。

通常VG模型参数的求解方法主要有：王金生等[9-10]采用最小二乘法和非线性单纯形法相结合的方法，研究了土壤持水曲线的滞留特征；徐绍辉等[11]用最小二乘法结合Picard迭代法拟合出了砂质黏壤土持水曲线的VG方程参数；李春友等[12]借助单纯形调优方法研究了VG方程参数的拟合问题；马英杰等[13]采用非线性阻尼最小二乘法拟合了VG方程的参数，以上方法对方程参数的求解均提供了较好的方法，但全局收敛性相对较弱。近年来，陈大春等[14]运用随机粒子群算法求解VG方程的参数，该方法具有较强的全局搜索能力和较高的求解精度，但迭代次数较多，计算量较大。因此，寻求一种全局搜索能力较强，计算精度强的寻优方法是VG模型参数求解亟待解决的科学问题，本文采用动态调整惯性权重的粒子群算法，并将优化准则函数的选取和全局寻优策略两者相结合以期达到收敛性强、计算精度高特点。并通过实测数据验证上述三种准则函数的优化效果，为该问题的求解提供有益的参考。

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参数的非线性超越函数的参数优化问题。

1.2 准则函数

在VG模型参数寻优中，优化准则的合理选取是关键之一。由于在实际寻优过程中模拟值不可能与实测值完全吻合，当模拟土壤含水率点据分布于实际土壤含水率曲线的两侧，必然会存在一定的离差Δ=θ模拟-θ实测，由最小二乘原理，首先可按绝对标准差最小准则(准则1；公式(2))拟合：

(2)

式中：θi(hi)为实测土壤含水率，cm3/cm3；θp(hi，X)为模型模拟的土壤含水率，cm3/cm3；hi为实测土壤基质势，cm；X为待优选的参数(用行向量表示)；n为实测数据组数。

其次，可按相对标准差最小准则(准则2；公式(3))拟合：

(3)

在更一般的情况下，应该有：

(4)

显然，λ∈[0，-2]当λ=0时，即为绝对标准差最小的拟合准则，其本质是取各经验数据点据均为1的权重，当为λ=-2时，即为相对标准差最小的拟合准则，其本质是取各数据点据随[θi(hi)]-2变化的权重，由此可以推断，当取中间值，即λ=-1时的拟合准则为：

(5)

上述方程其本质是取各经验数据点据随1/[θi(hi)]变化的权重，这实际上是一种考虑了相对标准差最小与加权耦合的绝对标准差最小的拟合准则(准则3；公式(5))。

1.3 动态调整惯性权重的粒子群算法

1.3.1 基本粒子群算法

粒子群优化(Particle Swarm Optimizer，PSO)算法作为一种新的基于群体智能的演化优化方法，它第一次是由美国普渡大学的Kennedy和Eberhart博士于1995年提出的[15-18]。PSO的基本概念来源于对鸟群捕食行为的研究，因此也可称作鸟群算法，是集群智能的代表性方法之一。该算法优点在于可并行处理以及鲁棒性好，既简单易实现又具有深刻的智能背景。因此，PSO算法在提出后就引起了广大学者的广泛关注，作为一种新的并行优化进化算法，PSO算法可以解决大量非线性、不可微、非连续和多峰的优化问题，且已广泛应用于科学和工程领域，如模糊系统控制、神经网络训练、函数优化、模式分类等[17]。

假设在D维的目标搜索空间，存在m个粒子可以组成一个群落，其中第i个粒子表示为一个D维的向量xi=(xi1，xi2，xi2，…，xiD),i=1,2,3,…,m,即第i个粒子在D维搜索空间中的位置。将xi带入目标函数即可计算出适应度值，进而以适应度值的大小判断xi的优劣。第i个粒子的“飞翔”速度也可作为一个D维的向量，记作vi=(vi1，vi2，vi3，…，viD)，将第i个粒子迄今为止搜索到的最优位置记为pi=(pi1，pi2，pi3，…，piD)，在整个粒子群中，所有粒子所搜索到的最优位置记为pgi=(pg1，pg2，pg3，…，pgD)。对于第k次迭代，每个粒子的进化采用公式(6)进行操作：

(6)

(7)

式中i=1，2，…，m是粒子数；d=1，2，…，D是向量维数；k=1，2，…，n为迭代次数；学习因子c1和c2为非负常数，分别表示调节向全局最好粒子、个体最好粒子方向所飞行的步长值；r1与r2是介于(0，1)之间的任意数，ω是惯性权重。根据具体问题确定迭代中止条件，一般选择最大迭代次数或者粒子群所搜索到的最优位置满足预定值来作为终止条件。基本PSO算法中需要调整的参数不多，而且采用实数编码可使其操作简便, 因此使用起来较为方便，但其具有容易陷入局部极值点，搜索精度不高的缺点。

1.3.2 动态调整惯性权重的PSO算法

当公式(6)中的惯性权重ω比较大时，vid靠近全局最优点pg和靠近历史最优点pi的机会较小，则vid会偏离当前最优点，从而有可能探索更好的最优点，这在早期中是非常有利的，但在晚期是不利的，这是由于晚期不利于全局最优点的收敛。然而当ω较小时，正好与上面情况相反。所以，ω值的大小对于收敛与探索的影响是一对矛盾。故而在应用中先将ω的初始值作为最大值，并使其随着迭代次数的增加而线性递减至最小值，以达到所期望的优化目标，权重函数ω的值用下式来确定：

(8)

式中：ωmax、ωmin分别表示ω的最大值和最小值；niter、nitermax分别表示当前迭代次数和最大迭代次数。但以上权重函数的设置也存在不利的一面，假设在迭代初期就已找到全局最优点，可能会因其权重过大而跳出这个最优点，因而不在其周围探索，导致算法收敛速度降低。若想提高收敛速度就必须赋予靠近最优点的部分粒子较小的权重，使粒子在靠近pg领域范围进行微小探索，而不是承担更大范围的探索。赋予远离最优点的其它粒子以较大权重，令其承担更大范围的探索任务，从而进一步去探索更优的点。因此，惯性权重ω的值不但随着迭代次数而变化，而且对不同粒子是根据其与最优粒子之间的距离大小而变化。假设Iig表示第i粒子与全局最优粒子Pg的距离，则惯性权重可采用公式(9)调整：

(9)

式中：Imax、Imin分别为粒子距全局最优值的最大距离和最小距离。

1.3.3 改进PSO算法在VG方程参数估计中的实现

对VG方程参数可理解为在土壤水力参数范围内随机生成的一组水深序列，以此作为初始粒子种群，再根据不同的优化准则设定目标函数，计算每个粒子的适应度，然后按照式(6)、式(9)进行迭代计算，直至达到最大迭代次数，具体的迭代步骤为：

步骤1 确定需要优化变量的取值范围，θs∈[0.4,0.6]，θr∈[0,3]，α∈[0,1]，n∈[0,10]。

步骤2 设定PSO的计算参数：粒子群体大小、学习因子c1和c2、最大迭代次数nitermax、惯性权重的最大值ωmax和最小值ωmin，粒子最大距离Imax和最小距离Imin。

步骤3 在优化变量θs，θr，α，n的取值范围内，采用蒙特卡洛随机方法产生初始群体粒子个体。

步骤4 根据不同的优化准则，计算出群体中每个粒子的适应度值。

步骤5 依据每个粒子计算出来的适应值大小，判断该粒子迄今为止所经历的最好位置(即准则函数最小的土壤水力学参数)以及整个粒子群迄今为止搜索到的最好位置(全局最优解)。

步骤6 计算出每个粒子相对全局最优解的距离并根据式(9)计算惯性权重。

步骤7 根据式(6)，更新所有粒子的速度和位置，若新速度超出最大速度限制[vmin,vmax]，则设新速度为最大速度。

步骤8 结束条件的判断。若达到最大迭代次数，则停止计算，输出结果；否则，转至步骤4继续计算，直至满足结束条件为止。

为验证本文提出算法的全局收敛性及不同优化准则的模拟效果，特以文献[3]中所给的实测资料作为本文方法的计算原始数据，对以上过程编程计算。取计算参数为：粒子群体大小20，最大迭代次数1000，初始惯性权重0.8，最终惯性权重0.2,学习因子c1=c2=2，最大粒子距离Imax=1最小粒子距离Imax=0.2，经过MATLAB 7.1对以上算法进行编程计算[19]，采用改进PSO算法连续进行10次独立计算，从中选取最优结果作为最终优化结果。

由以上分析可知，VG模型参数寻优中要想达到计算速度快、精度高的双重要求就必须将准则函数与算法两者相结合。为凸显改进PSO算法的优点及准则函数的适应性，特选取目前比较典型的几种方法，包括文献[10]的非线性单纯形法、文献[13]的非线性阻尼最小二乘法和文献[14]的随机粒子群算法(上述求解方法的准则函数均采用准则1)与本文改进的PSO算法和准则1(改进的PSO算法1)、改进的PSO算法和准则2(改进的PSO算法2)和改进的PSO算法和准则3(改进的PSO算法3)共6种方法进行比较研究，以期检验改进PSO算法的优越性及准则函数的适用性。

1.4 试验数据

提供试验土壤为陕西杨凌地区的黏壤土，自然风干后经过1 mm孔径的土筛选后利用沉降法进行颗粒分析(见表1)。土壤水分特征曲线采用高速离心机测定(见表2)。

表1 试验土壤颗粒级配

表2 土壤水分特性曲线实测数据

2 结果与讨论

2.1 不同优化准则的比较

由表3可知，不管采用哪种优化准则，即本文述及的3种优化准则，改进的PSO算法明显优于非线性单纯性法和阻尼最小二乘法，同时也优于随机粒子群法的优化结果。由此可知，本文提出的改进惯性权重模型的PSO算法能更好的寻求全局最优解，是一种很好的全局搜索策略。此外，还可以看出，3种优化准则其拟合精度均较高，这说明在现行VG方程参数优化中采用绝对标准误差最小准则基本上是科学、合理、可行的。但仔细分析表3数据可以看出，采用准则3虽然计算误差不能保证绝对标准误差和相对标准误差同时达到最小，但其绝对标准误差小于准则1的计算结果，且相对标准误差小于准则2拟合的结果。

从表4中最大相对误差绝对值和平均相对误差绝对值的比较可知，从模型求参方法来看，改进惯性权重粒子群算法模拟精度均高于文献[10,13]的方法。由平均相对误差来看，改进PSO与文献[14]相差不大，但以相对误差的绝对值来看，改进PSO方法优于随机粒子群算法。从优化准则来看，准则3拟合结果中平均误差为5.675%，小于准则1优化结果(平均相对误差为7.899)，略大于准则2的优化结果(平均相对误差为5.338)，而最大相对误差的绝对值明显小于准则1和准则2。

通过以上分析可知，加权耦合的绝对与相对标准差最小准则，解决了仅以绝对标准差最小准则或仅以相对标准差最小准则进行参数拟合所存在的某些不足，故加权耦合的绝对与相对标准差最小准则是一种在土壤VG方程参数拟合问题中值得推荐的优化准则。

表3 不同方法确定的粉壤土Van Genuehten方程参数

表4 不同方法模拟θ-h数据与实测数据对比

2.2 VG模型参数的优化

以上分析可知，加权耦合优化准则优于绝对标准差最小准则和相对标准差最小准则，由此可将加权耦合优化准则运用于本文试验数据的优化拟合中，并采用改进惯性权重的粒子群优化算法进行参数的寻优。

由表5和表6可知，最大相对误差的绝对值仅为3.135，且分布在土壤水吸力1000 cm以上，该土壤水吸力区间对应的含水率变化区间很小仅为0.2228～0.1636之间，在实际中应用相对也较少，而在土壤水吸力小于1000 cm以下时，拟合精度相对较高。总体来说，采用加权耦合优化准则的改进PSO算法比非线性单纯形法和阻尼最小二乘法及随机粒子群算法要优越，且不受初值的影响。

表5 根据加权耦合优化准则确定的黏壤土Van Genuehten方程参数

表6 模拟值与实测数据对比分析

3 结 论

(1) 本文提出了一种全局优化策略—改进惯性权重粒子群优化算法，通过实例计算表明，该算法全局搜索能力强，优于现有几种优化方法。

(2) 本文对Van Genuchten方程参数寻优中准则函数进行比较研究表明，加权耦合绝对标准差和相对标准差最小准则虽然计算误差不能满足绝对标准误差和相对标准误差同时达到最小，但可以解决仅以绝对标准差最小准则或仅以相对标准差最小准则进行参数拟合所存在的某些不足，可作为Van Genuchten方程参数拟合中较优的优化准则。

(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] 夏卫生，雷廷武，潘英华，等.土壤水动力学参数研究与评价[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] 李春友，任 理，李保国.利用优化方法求解Van Genuchen方程参数[J].水科学进展，2001，12(4)：473-478.

[13] 马英杰，虎胆.吐马尔拜，沈 冰.利用阻尼最小二乘法求解Van Genuchten方程参数[J].农业工程学报，2005，21(8)：179-180.

[14] 陈大春，马英杰.基于随机粒子群算法的Van Genuchten方程参数优化求解[J].农业工程学报，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最优化计算[M].北京：电子工业出版社，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)；黄土高原土壤侵蚀与旱地农业国家重点实验室主任基金项目(A314021402-1619)

孙晓敏(1986—)，女，陕西宝鸡人，主要从事水土保持与生态修复方面的研究工作。 E-mail：sunxiaomin428@126.com

S152.7

A

1672—1144(2017)05—0088—06