基于集合卡尔曼滤波法的二维土壤水流状态变量和参数联合估计

刘 琨，黄冠华

（1.中国农业大学 水利与土木工程学院，北京 100083；2.中国-以色列国际农业研究培训中心，北京 100083）

1 研究背景

土壤水流模型的预测精度很大程度上依赖于土壤水力特性参数的准确估计。常用的模型参数反求方法有Levenberg-Marquardt算法［1］，Shuffled complex evolution算法［2］等，它们均以观测值与模拟值差异最小为目标，得到一组最优的参数为估计值。然而这些反求优化算法将预测结果的不确定性归因于模型参数，可能导致参数的次优解，进而降低了模型的预测精度。

集合卡尔曼滤波（EnKF）方法［3］可以显式地考虑模型输入、输出以及模型结构等因素的不确定性，与传统的卡尔曼滤波相比，EnKF方法引入了集合预报技术，通过蒙特卡洛方法计算预报误差的协方差，因此可以处理高维和非线性的问题。此外EnKF方法对观测数据集进行时间序贯处理，降低了数据的存储量和计算量。因此，EnKF方法被广泛应用于水文模型参数估计的研究中［4-6］。

非饱和带水力参数与状态变量之间的高度非线性关系增加了其参数估计的难度，近年来基于En-KF方法已开展了一些非饱和带水力模型的参数估计研究。Brandhorst等［7］研究了在模型参数高度不确定且不可识别的条件下非饱和区模型的集合卡尔曼滤波系统的数据同化的潜力。Li和Ren［8］采用En-KF方法研究了不同土壤条件下状态变量同化以及非均匀土壤条件下的多参数估计问题，并分析了样本数量和观测误差等因素对同化效果的影响。Song等［9］评估了三种迭代集合卡尔曼滤波算法的同化效果，分析了不同观测数据类型，观测误差方差，集合大小和阻尼因子等因素对同化效果的影响。

目前关于集合卡尔曼滤波法的非饱和带参数估计的研究大多考虑一维问题［10-12］，然而在实际问题中由于复杂的边界条件（例如沟灌、滴灌等）以及土壤剖面非均质性，土壤水流运动多为二维问题。此外与一维问题相比，二维土壤水流模型模拟区域离散后通常产生较多的计算节点，而在实际中考虑成本问题，布置观测点的数量有限，因此如何合理的布置观测网以实现观测点信息的高效利用具有重要意义。总之，EnKF方法在二维土壤水流运动条件下的土壤水力特性参数估计还需进一步研究。本研究基于EnKF方法，开展二维土壤水流模型状态和参数联合估计研究，探讨在线源入渗条件下EnKF方法对粉壤土、壤土和砂壤土的饱和导水率和进气值参数估计精度以及对土壤压力水头的同化效果，分析了观测点布置方式和观测点数量对同化效果的影响，以期指导田间观测网点的合理布置。

2 理论与方法

2.1 数据同化系统数据同化系统由模型算子、观测算子和同化算法组成。首先模型算子根据上一时刻的状态向量向前预测得到状态向量的预测值；观测算子将状态向量转化为观测值所对应的模型预测值；最后同化系统利用观测信息，对预测值进行同化运算，得到状态向量的分析值。分析值作为下一时刻的模型算子输入变量，重复以上过程直至模型模拟结束。

2.1.1 模型算子 本研究中的模型算子为二维土壤水流运动模型，即二维Richards方程:

式中:θ为体积含水率，cm3/cm3；h为压力水头，cm；K为非饱和导水率，cm/min；KijA为无量纲各向异性张量KA的分量；KizA为垂向方向无量纲张量KA的分量；xi和xj是空间坐标，cm；S为源汇项。非饱和导水率函数采用van Genuchten-Mualem模型表示［13-14］。本研究中使用HYDRUS-2D［15］的水流运动模块模拟二维土壤水流运动。由于本研究观测数据和模型预测值均为压力水头，因此观测算子为对应观测点位置处为1，其余元素为0的矩阵。

2.1.2 同化算法 采用扩展状态变量法将模型状态变量和土壤水力特性参数包含在一个矩阵中形成状态向量，使用EnKF方法同时更新状态变量和待估计参数。EnKF方法中通过向模型状态变量和参数加入高斯白噪声生成状态向量集合，利用集合计算出模型预报误差协方差矩阵，再基于蒙特卡洛抽样方法来估计高维非线性动力学模型中的状态变量和参数。在具体的实施过程中，EnKF方法包括预测和分析两个步骤，首先根据前一时刻的模型状态向量生成当前时刻状态向量的预报值:式中:Xi,tf+1为t+1时刻第i个集合成员模型向量的预测值；Xi,at为t时刻第i个集合成员模型向量的分析值；M为模型算子，即二维土壤水流运动模型；ei为模型误差向量，其满足均值为0、方差为Q的高斯分布。

当有观测数据时，利用观测数据计算得到当前时刻状态向量的分析值:

式中:H为观测算子；Pef为预测值观测误差协方差；Re为观测误差协方差；di，t+1为t+1时刻各集合成员的观测向量。详细的EnKF计算流程请参见文献［3，16］。

基于EnKF土壤水力参数估计流程如下:（1）模型初始设置，包括模型参数，初始条件和边界条件设置；（2）根据集合数量通过随机数函数抽样得到待估计参数集合；（3）每一时刻，HYDRUS模型向前计算，得到模型向量的预测值；（4）当该时刻有观测数据时，对观测值进行扰动得到观测向量；（5）利用公式（3）计算得到状态向量的分析值，即更新压力水头和水力参数，作为下一时刻模型预测的初始条件；重复步骤（3）—（5）直至模拟时间结束。

2.2 数值试验采用数值实验的方法探究EnKF方法的同化效果。设置人工假想算例，实验方案如下:地表设置长度为20cm的线形入渗源，模拟区域为50 cm×50 cm的均质土壤剖面。考虑粉壤土、壤土和砂壤土三种土壤类型，三种土壤的土壤水力特性参数根据Carsel和Parrish取值［17］（见表1）。粉壤土和壤土的模拟时间设为400 min；由于砂壤土的导水率较大，入渗时间较短，模拟时间设为200 min。模型的初始条件为均匀分布的土壤压力水头剖面，水头值为-100cm；上边界条件为定水头边界，压力水头为2 cm；下边界条件为自由排水边界。

表1 土壤类型及其土壤水力特性参数

待估计参数选取土壤水力特性参数中变异性最大的饱和导水率和进气值参数［17］。为了探究EnKF方法的参数估计效果，参数的初值设置均大于参数的真值，假定待估计参数服从正态分布，初始统计特征为lnKS～N（1.5a，0.01）、α～N（1.2a，0.0001）。其中，a为参数的真值。由于参数KS的方差较大，为保证KS均为正值，因此对KS进行了取对数处理。观测数据通过以参数真值模拟水分入渗过程，每隔20 min记录观测点的压力水头，添加一定的扰动作为具有观测误差的观测值，假设扰动服从正态分布N（0，0.001）。研究中样本数目设置为150以保证同化效果。

2.3 观测点布置方案为了研究观测点位置对同化效果的影响，设置了不同观测点布置方案，观测点布置见图1。其中方案1—4中观测点垂向布置，分别距离模拟区域中心轴0、5、10和15 cm；方案5—8中观测点水平向布置，分别距离地表10、20、30和40 cm。

为研究观测点数量对同化效果的影响，选择每种土壤最优的三种观测点布置方案进行组合确定最终的布置方案。采用平均相对误差（MRE）评价参数估计效果，MRE值越小，表示参数越快收敛于参数真值。采用绝对误差（AE）评价土壤压力水头分布剖面的同化效果。

式中:Ya为同化后的参数样本均值；Y为参数真值；N为同化次数。

图1 不同观测点位置布置方案（实心点为观测点，单位:cm）

3 结果分析与讨论

3.1 观测点位置影响图2为不同观测点布置方案饱和导水率的演化。由图2可知粉壤土条件下，垂向布置观测点的方案比水平向布置的方案对参数KS估计效果更好，最优的3种布置方案为1、2和3（见表2）。壤土条件下，除了方案8不能较好地估计参数KS，其他方案参数KS都可以收敛于真值，最优的3种布置方案为1、5和6。砂壤土条件下，8种布置方案均可以较好地估计KS，最优的3种布置方案为2、3和4。这可能是由于粉壤土的导水率较小，水平布置观测点时当水分运动还未到达观测点深度，同化系统未能利用观测信息更新参数，参数估计偏差较大，尤其随着观测点埋深的增加，估计效果较差。而对于壤土和砂壤土由于其导水率较大，入渗水分可以很快到达观测点深度，因此水平向布置方案也可以得到较好的参数估计结果。

参数KS标准差随着同化次数增加逐渐减小（见图2）。垂向布置方案的KS标准差变化趋势类似，经过3次同化标准差值已经很小；水平向布置条件下，不同方案的标准差变化差异较大，随着观测点埋深的增加，参数估计效果变差，标准差降低的速度变小。这是由于当入渗水分没有运动到观测点埋深时，观测点的压力水头不变，同化系统尚未利用观测信息更新状态变量和参数，因此参数估计结果较差。图3为不同观测点布置方案进气值参数的演化，进气值参数α的演化与KS的演化趋势一致，KS同化效果较好的方案，对α的估计效果也较好，这是由于同化计算中参数KS和α是同时更新的。与参数KS不同，8种方案参数α的标准差均下降较快，这是由于参数α的初始方差设置较小。

图2 不同观测点布置方案饱和导水率演化

图3 不同观测点布置方案进气值参数演化

表2 不同观测点位置KS估计MRE

表3 不同观测点位置α估计MRE

由于3种土壤条件下压力水头误差结果类似，我们仅给出了壤土条件下不同观测点布置方案土壤压力水头误差分布（见图4）。垂向布置条件下，方案1—4随着观测点到中轴距离的增加，误差分布也随之偏移，观测点两侧部分的AE较大。水平布置条件下，土壤剖面0～30 cm区域的AE较大，其中方案7和8的同化效果较差。王文等［10］和Montzka等［18］也得到类似的研究结果，其研究表明使用土壤表层观测信息，同化效果仅限于一定深度的土壤。随着观测点埋深增大，参数估计偏差增大，导致土壤压力水头的同化结果较差，并且这种参数估计偏差引起的误差不会在同化过程中被消除［19］。以上结果表明有观测点的区域误差较小，说明同化系统可以有效地利用观测信息改善土壤压力水头的预测结果。

图4 壤土不同观测点布置方案土壤压力水头误差分布（单位:cm）

3.2 观测点数量影响在观测点位置影响结果的基础上选择每种土壤最优的三种布置方案组合得到新的布置方案，观测点数量分别为4、8和12，观测点布置见图5。

图6给出了不同观测点数量条件下参数KS和α的演化。随着观测点数量增加，参数更快趋近于其真值。观测点数量为12时参数的MRE值最小值，粉壤土、壤土和砂壤土KS的MRE分别为0.097、0.065和0.021，α的MRE分别为0.075、0.073和0.059。这是由于同化系统利用更多的观测信息更新参数，提高了参数的估计精度。图7给出了土壤压力水头误差分布剖面。土壤压力水头误差随着观测点数量的增加而减小，这是由于同化系统利用更多的观测信息实时校正土壤压力水头剖面；此外增加观测点数量可以得到更好的参数估计值，进而提高了土壤剖面压力水头的预测精度。这与兰天等［20］结果一致，其研究表明增加观测点数量可以有效地改善同化效果，参数收敛速度也更快。

本研究在数值试验中仅考虑了3种土壤类型。这是由于质地更粗的土壤如砂土，土壤导水率较大导致入渗时间较短，当整个土壤剖面饱和时样本之间无差异，导致在同化过程中无法进行矩阵求逆运算，导致计算过程终止；而对于细质土壤如黏土，在长时间的积水入渗条件下，如果参数的初始均值和方差以及观测误差方差等变量设置不合理会导致模拟计算失败，这是由于同化过程中部分集合的土壤参数偏差较大，导致模型求解不收敛。对于细质土壤宜采用较小的参数方差以及较大的观测值方差以保证同化过程中模型的顺利运行。而本研究重点关注观测点布置方式和数量对同化效果的影响，因此为了保证参数统计变量和观测误差方差一致以及模型的正常运行，选取了粉壤土、壤土和砂壤土3种土壤类型进行试验。在之后的研究中我们将考虑更多的同化影响因素，如集合数量和同化时间间隔等，并通过设置合理的参数统计变量和观测误差方差等方法以涵盖更多土壤类型开展进一步研究。

图5 不同观测点数量布置方案（实心点为观测点 单位:cm）

图6 不同观测点数量条件下饱和导水率和进气值参数演化

图7 不同观测点数目条件下土壤压力水头误差分布（单位:cm）

与传统的参数估计方法相比，EnKF方法在处理不确定性方面有着明显的优势。例如在Levenberg-Marquardt算法中，观测数据被认定为是衡量预测精度的唯一标准，模型的不确定性完全归因于参数，并没有考虑观测数据误差。而在模型构建过程中对物理过程的认识不足和简化以及观测手段的局限，模型结构和观测数据都具有不确定性，因此仅将不确定归因于参数是不恰当的。EnKF方法综合观测数据和模型预测的不确定性，可以得到参数的优化值。本研究中我们已通过数值实验从理论上证明了状态变量与参数联合估计效果，在未来的研究中还需结合田间实验进一步检验该方法的精度。

4 结论

本文采用EnKF方法开展二维土壤水流运动模型状态变量和参数联合估计研究，采用数值实验方法研究了粉壤土、壤土和砂壤土在线源入渗条件下EnKF方法对饱和导水率和进气值参数估计以及压力水头的同化效果，分析了观测点布置方式和观测点数量对同化效果的影响。主要结论如下:

（1）观测点的布置方式应捕捉二维土壤水流运动过程，才能得到较好的参数估计值。粉壤土条件下，宜采用垂向的观测点布置方式；壤土和砂壤土条件下，在0～30 cm深土壤中水平向布置观测点可以得到较好的参数估计值。

（2）入渗条件下观测点布置应尽量靠近地表，使得同化系统可以尽早利用观测信息更新状态向量，参数可以更快地收敛于真值，但压力水头的同化效果仅限于一定深度的土壤。

（3）增加观测点数量可以使参数更快地收敛于其真值，进而提高土壤剖面压力水头的预测精度。