智能算法反演土壤水分特征曲线参数的应用

王 正，孙兆军，王 旭

(1.宁夏大学土木与水利工程学院，银川 750021；2.宁夏大学新华学院，银川 750021)

0 引 言

近年来，随着人们对环境质量重要性的认识的提高，研究水分和可溶性盐在非饱和土壤中的运移机理已成为众多领域的热点[1]。而可溶性盐则是通过水分在土壤中运移来进入地下水中的，因此，土壤中的水分运动对其有着至关重要的影响。为了更准确地从定量的方面来揭示水分和盐分的运移规律，就需要获取更为精确的输入参数，而土壤水分特征曲线(Water retention characteristics)是其中最为重要的参数之一[2]。

对于刻画土壤水分运动与含水率的土壤水分特征曲线问题，国内外学者先后提出了许多数学模型和计算方法，其中吻合度最好的就是Van Genuchten模型[3](简称VG模型)，而VG方程是一个高度非线性的数学模型，含有众多的未知参数，其表达式如下：

(1)

另外,求解方程(1)中的参数问题，本质就是利用实验实测数据来拟合这些参数的问题。一般情况下，实测数据的样本容量往往较小，故采用如下以回归估计值与实际观测值的偏差平方和最小为原则构造鲁棒性较好的优化正则函数[4]：

(2)

式中：θi为实测土壤含水率；N为实验次数。

从本质上讲，式(2)隶属于非线性问题的范畴，很多方法都可解决该问题。一般地，有直接法和间接法。通过实验数据拟定土壤水分特征曲线，也就是直接法，是一种最快捷的方式。而直接法中又有许多实验室方法和田间方法，例如测定土壤基水势就有负压力法、砂型漏斗法、压力仪法等。虽然这些实验方法可以使用，但会大量消耗人力物力，实验时间漫长，并且不能保证测定范围和实验精度[5]。间接法是指利用实验测得的有限的样本，在已知模型上结合智能算法和原先假定的参数、源汇项和边界条件等约束求得较优参数。最近，将待求解问题与智能算法相结合的方式被广泛应用，而且计算准确性较高。

1 参数反演算法现状

随着计算机和数学理论的不断发展，对于参数优化的问题，很多学者提出了不同的优化方法，例如仿生算法、层次贝叶斯方法、集合卡尔曼滤波算法等。本部分将着重介绍21世纪以来用于土壤水分特征曲线参数反演的不同算法，并研究它们在实验精度和算法性能方面的异同。

1.1 仿生算法

1.1.1 萤火虫算法

萤火虫算法(Firefly Algorithm，FA)是根据萤火虫独特的生理构造和生物习性总结出来的一种算法，其主要思想就是萤火虫之间会利用发光来相互吸引，它们聚集的地方就应是最优解。但是由于原始的算法会在最优值附近出现振荡的可能，随后付强[6]等人提出了一种改进的萤火虫算法。建立了一种约束判定条件：每两个萤火虫的距离lij小于某一特定值时，减小其步长λ。

设定λ的起始值为0.3，当lij大于0.3时，可用下式计算其位置：

(3)

其中，符号意义同参考文献[6]中表示。

其主要步骤如下：

(1)设定算法参数的初始值；

(2)随机生成种群的初始位置；

(3)决定萤火虫移动方向；

(4)用式(3)修正萤火虫位置，找寻最佳位置；

(5)在更改位置后，重新计算相对亮度；

(6)结果与判定条件进行对比，决定是否继续进行算法；

(7)算法结束，输出结果。

实验结果表明，用这种算法得到的计算结果与实测值有很好的拟合度，尤其是在较低的土壤水吸力条件下，误差相当小。但是萤火虫算法的理论支撑还不够，而且也没有一个完备的数学模型，所以仍需要完善现有的理论，与其他算法糅合也能够提高算法的性能。

1.1.2 粒子群算法

粒子群算法是一种基于进化的计算方法(Particle Swarm Optimization，PSO)，该算法是建立在已发现的飞鸟集群活动规律基础之上，然后利用人工智能得到的简化模型。该算法的基本思想是由个体与种群共享其获得的“食物”信息，使整个群体都向求解空间中的“食物”所处位置进行运动，从而产生从无序到有序的演化过程以获得最优解。

刘利斌[4]等、马亮[7]完善了最初的算法，虽二人提出的算法结构有些不同，但其大致步骤是相类似的，其主要步骤总结如下：

(1)参数初始化。计算各粒子对应的目标向量、确定种群规模、随机产生θs、θr、α、m等参数的初始值。

(2)从初始值开始进行迭代，得到位置更优的初始值；

(3)计算粒子适应度；

(4)持续更新适应度和位置；

(5)对比适应度值，得到全局最优位置；

(6)比对结果和判定条件，决定是否继续进行。

计算结果表明，改进后的粒子群算法有了更快的搜索速度，并且其精度也有很大的提高，因此，该算法具有较强的全局优化能力和计算结果准确的特点；另外，迭代次数少、鲁棒性强、效率高也是该算法最主要的优点。

1.2 层次贝叶斯方法

层次贝叶斯方法(Hierarchical Bayesian Method，HBM)于1996年由Berliner[8]提出，将复杂问题利用条件概率理论，逐级分解为不同层次的简单后验概率求解问题[9]。而且各层次之间并不是孤立的，它们是由条件概率联系起来的。层次贝叶斯方法在数据分析方面具有很大的潜力，现已经成功应用于其他领域[10-13]。

一般地，层次贝叶斯方法有数据S、过程Q和参数C等三层，数据S包括实地实测数据、反演同化数据和模型计算数据等；过程Q表示土壤湿度等参数的时空分布。其主要原理为：定义数据、过程和参数三层为随机变量，并构建条件概率模型，在已有数据条件下，将数据同化问题转化为推理过程和参数的后验概率分布：

(4)

式中：p(×)为某事件的概率；p(a|b)表示b条件下a发生概率。

根据式(4)，层次贝叶斯方法数据模型为第一层，过程模型为第二层，参数模型为第三层，因此，定义了这三层模型之后，就能推导后验概率。另外，这三个模型还能够进行拆分，这样做的目的是用其他的观测数据或过程就能完全确定这些模型，这就实现了用简单条件概率确定复杂的模型，层次贝叶斯方法也可以用于确定复杂模型的参数[14]。

(1)数据模型。数据模型表示X和θD在真实过程条件下与数据Y的概率关系，其模型为：

(5)

式中，数据Y=(Y1,Y2,…,Yn)表示不同数据源数据，X=(X1,X2,…,Xn)为每个数据源对应尺度的过程，θD为参数，p(Yi|Xi,θD)表示为第i个数据源的数据模型。当由不同分辨率的多源数据构成输入数据时，层次贝叶斯方法可分别建立相对应的数据模型，并依赖于真实过程和参数，这就是说层次贝叶斯方法在解决多源数据问题方面具有一定优势。

(2)过程模型。层次贝叶斯方法中最为关键一步就是构建过程模型，而过程模型需要对物理过程有大量的经验信息和深刻的理解，最后建立过程X的时空分布模型p(X|θp)。通常来说，可将过程模型拆分为若干层次的子模型，其形式为：

(6)

式中：X=(X1,X2,…,Xn)为多尺度不同源数据所对应的过程；θp为参数。

(3)参数模型。假定参数之间均相互独立，则模型可定义为：

p(θD,θp)=p(θD)p(θp)

(7)

将式(5)～(7)代入式(4)中，数据同化要推理的概率为：

(8)

用层次贝叶斯方法来解决数据同化问题是一种非常有效的方法[15]，层次贝叶斯方法有以下几个明显区别于其他数据同化算法的特点：

(1)以贝叶斯理论为基础，将估计状态的后验概率作为目标，同时完全考虑给定的观测数据和模型数据；

(2)同化的数据中难免存在不确定性，而且对结果有重要影响，该方法能克服到这些缺点，也能脱离线性和高斯假设的束缚；

(3)分层考虑各个模型，并整合为一体，这样处理的方式便于计算，而且在理论和实际应用中更加合理方便。

1.3 集合卡尔曼滤波算法

Evensen[16]基于卡尔曼滤波提出了集合卡尔曼滤波(Ensemble Kalman Filter，En-KF)算法，用集合卡尔曼滤波研究土壤水分时，主要关注的是状态变量的估计和预测[17-21]。

En-KF是基于蒙特卡洛的方法形成每一个元素的数据向前积分来预报误差的统计量集合。一般用均值来近似最佳估计，用集合的取值范围来表征参数估计的不确定性。

用EnKF来估计参数的主要步骤如下：

(1)形成初始变量和初始猜想集合。

(2)预报。集合中的每个元素向前积分，得到土壤持水率以矩阵形式存在的预报结果，并存储在搜索空间中。

(3)产生观测值的集合。

(4)分析。利用观测算子分析其结果，将分析的结果存储到搜索空间中。

(5)判定结果是否满足要求。如果不满足，返回第2步，否则停止并输出。

Evensen[22,23]等给出更详细的程序。

EnKF方法优化的结果不受初始条件的影响，且具有很强的适应性和更快的收敛性，还能够同时考虑土壤水分函数模型以及数据观测不确定性，大大降低了获得伪优参数值的风险。

2 存在的问题

进入21世纪以来，我国在土壤水分渗透和水盐运移方面取得了许多成绩，但是在国际科技不断发展的大环境下，我国在土壤水分数值模拟方面仍不能跟上世界的步伐。急需一大批科技工作者投入到这项事业当中来，同时，也需要从更高的角度冷静地思考仍未解决的问题。存在的主要问题归纳如下：

(1)严重依赖仿真技术，缺乏地质考察经验和理论。由于实验不仅可以发展理论还可以验证理论，因此，本文认为实验是提升理论最有效的方式之一。建议国家有关部门要加大投入，要让更多的人才参与到这项事业中的，同时让他们戒除浮躁，安心工作。大量的经验表明，只有通过大量的试验，获取更多的数据样本，才有利于发掘新的认识。

(2)轻视具体地质条件，随意主观地删改数据，为了达到较好的拟合效果，随意忽视一些现象或边界条件。这种潮流，在现在的学术界尤为盛行。没有深入考察当地的地质情况下，就生搬硬套，必然只会产生片面的、带主观臆想的甚至是错误的结论来。无数实践证明，拟合效果不佳往往都是模型有缺失，而这又是不认真考察地质条件的结果。所以，结合当地真实情况建立符合实际的模型，对客观反映地区具体条件大有裨益。

(3)盲目追求参数拟合软件。现在有很多商用软件都可以进行参数反演的研究，但是不能过分依赖这些软件，不是说有了软件就万无一失了，而是要在软件的基础上加入符合当地特色的因素，只有这样拟合出来的结果才能令人信服。

(4)模型虽多，但是没有证明模型的可信度。早期一些学者提出来的模型和现如今的实际情况已相去甚远，原因就在于当时的地质条件已发生了变化，而且在主观认识上也有了不同，因此，以早期的模型来拟合现在的参数，恐怕也只能是刻舟求剑、朝花夕拾了。所以，要敢于挑战权威，在已有的模型基础上加以改进优化，这样才是科学发展的规律。

3 发展方向与展望

近年来，已将众多智能算法和技术应用到土壤水分曲线参数反演的研究当中来，但是随着更多更好的、更新的方法的提出，在这方面的探索仍然任重道远。下面提出几点发展的方向与展望：

(1)微波遥感技术可实时测量观测，而且对地表环境的要求不高，能避免其他方式在土壤水分监测应用中的缺点和短板，这就为研究土壤水分提供了一个全新的方法和途径。但是，由于目前的研究尚未成熟，所以有许多问题仍未完全解决，例如，一旦有植被附着在土壤的表面，微波会发生时空变异，会直接影响水分模型的参数；再如，雷达监测虽然具有一定的穿透能力，但是也只是针对浅层土壤而言，对于深层土壤水分的探测仍然无能为力，无法探测到作物的根系，这样会得到片面的结果。

(2)影响土壤含水率的因素很复杂，包括温度、湿度、有机质含量、矿物类型、滞后、阳离子交量和阳离子类型等，还需要进一步考虑这些因素对水分曲线的影响机理。

(3)土壤结构影响土壤储水性能。土壤入渗能力、持水能力及土壤水分特征参数直接受土壤孔隙数量与质量的影响，而团聚体含量又影响土壤孔隙数量与质量，所以要进一步来研究土壤结构性与土壤水分特征曲线的关系，也就要建立基于土壤传递函数的可以反映土壤入渗、蓄水能力的计算模型，为后续的研究提供简便的方法。

4 结 语

预测土壤水分特征曲线的间接方法在今后的研究中有很广泛的前景。而且，尤其在较大范围的模型研究中，间接方法的优势就突显出来。就目前来说，土壤水分特征曲线参数反演主要依靠的是智能算法，例如仿生算法、数据同化算法等，这些算法在不同程度上都具有很好的反演精度，但是在模型的使用和算法的结构上仍然需要进一步提高。另外，其他新技术的摸索和使用，例如遥感技术，在很大程度上也能提供更广阔的平台。

[1] 李 峰，缴锡云，李盼盼，等.田间土壤水分特征曲线参数反演[J].河海大学学报(自然科学版)，2009，(4)：373-376.

[2] 刘建立，徐绍辉，刘 慧.估计土壤水分特征曲线的间接方法研究进展[J].水利学报，2004，(2)：68-73.

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

[4] 刘利斌，欧阳艾嘉，乐光学，等.基于混合粒子群的土壤水分特征曲线参数优化[J].计算机工程与应用，2011,47(35)：218-221.

[5] 唐亚莉，董文明，赵晶晶，等.优化算法确定土壤水分特征曲线的分析[J].新疆大学学报(自然科学版)，2006，23(2)：240-243.

[6] 付 强，蒋睿奇，王子龙，等.基于改进萤火虫算法的土壤水分特征曲线参数优化[J].农业工程学报，2015,31(11)：117-122.

[7] 马 亮.基于改进粒子群的土壤水分特征曲线参数优化研究[J].节水灌溉，2014，(11):17-20.

[8] Berliner L M，Hanson K，Silver R. Hierarchical bayesian time series models[M]. Maximum Entropy and Bayesian Method.Dordrecht：Kluwer Academic Public，1996:15-23.

[9] Cressie N，Wikle C K. Statistics for Spatio-temporal Data[M]. New Jersey：Wiley，John & Sons，2010.

[10] Gelman A，Carlin J B，Stern H S，et al.Bayesian Data Analysis[M].New York：Chapman & Hall/CRC，2003.

[11] Gelfand A E, Diggle P, Guttorp P, et al. Handbook of Spatial Statistics[M].New York：Chapman and Hall /CRC，2010.

[12] Gelfand A E, Sahu S K. Combining monitoring data and computer model output in assessing environmental exposure[M].The Oxford Handbook of Applied Bayesian Analysis. Oxford：Oxford University Press，2009.

[13] Sahu S K，Yip S，Holland D M. Improved space-time forecasting of next day ozone concentrations in the eastern US[M]. Atmospheric Environment，2009，43( 3)：494-501.

[14] Liu Y，Gupa H V. Uncertainty in hydrologic modeling toward an integrated data assimilation framework[J].Water Resources Research，2007，43：10.1029 /2006WR005756.

[15] Wikle C K，Berliner L M.A Bayesian tutorial for data assimilation[J]. Physica D，2006，230( 1 /2) ： 1-16.

[16] Evensen G. Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics [J]. Journal of Geophyscial Research，1994，99：10 143-10 162.

[17] Reichle R，Walker J，Koster R，et al.Extended versus ensemble Kalman filtering for land data assimilation [J]. Journal of Hydrometeorology，2002，(3)：728-740.

[18] Crow W T. Correcting land surface model predictions for the impact of temporally sparse rainfall rate measurements using an ensemble Kalman filter and surface brightness temperature observation [J].Journal of Hydrometeorology, 2003，(4)： 960-973.

[19] Zhang S W，Li H R，Zhang W D，et al. Estimating the soil moisture profile by assimilating near-surface observations with the Ensemble Kalman Filter (EnKF) [J]. Advances in Atmospheric Sciences，2006，22(6)： 936-945.

[20] 黄春林，李 新.土壤水分同化系统的敏感性试验研究[J].水科学进展， 2006，17(4)： 457-465.

[21] 苟浩锋，刘彦华，张述文，等.评估集合卡曼滤波反演土壤湿度廓线的性能[J].地球科学进展，2010，25(4)：400-407.

[22] Evensen G.The ensemble Kalman filter：Theoretical formulation and practical implementation [J].Ocean Dynamic，2003，53：343-367.

[23] Evensen G. Sampling strategies and square root analysis schemes for the EnKF[J].Ocean Dynamic，2004，54：539-560.