无迹加权集合卡尔曼滤波及其在土壤湿度同化中的应用

付晓雷，余钟波，丁永建，蒋晓蕾，杨传国，鞠 琴

(1.福州大学土木工程学院，福州 350116;2.中国科学院西北生态环境资源研究院冰冻圈科学国家重点实验室，兰州 730000；3.水文水资源与水利工程科学国家重点实验室，南京 210098;4.河海大学水文水资源学院，南京 210098)

0 引 言

土壤湿度是气象、水文、农业等领域重要的物理状态变量，它对地表能量[1]以及降雨下渗[2,3]的再分配和水资源管理、农业灌溉等[4,5]方面具有重要影响。尽管土壤湿度在不同学科领域的重要性越来越受到关注，但是由于缺乏长期的全球高密度的土壤湿度观测数据，使其土壤湿度的研究受到了影响[6]。随着卫星遥感技术的发展，极大地提高了全球土壤湿度数据库的空间分辨率。然而，遥感反演产品仅限于距地表几厘米[7,8]，且易受土壤类型等影响，精度不高。此外，陆面水文模型可以利用气象、土壤、植被等信息模拟土壤水分的动态过程，且能够得到高时空分辨率的区域土壤湿度分布数据[9]。但是，由于模型物理过程概化的不精确性(包括大气模型和陆面模式及参数化方案)、土壤参数和初始状态的不准确性以及气象信息、陆面模式本身的不确定性，使得土壤湿度的模拟精度受到限制。

为了满足对高精度高分辨率土壤湿度数据的需求，综合观测、反演、模拟数据的优点，数据同化方法逐渐在土壤湿度研究领域得到了应用发展[10-15]。尤其是Kalman[16]提出针对线性系统的卡尔曼滤波方法之后，数据同化方法得到了快速有效的发展以及在土壤湿度研究中的应用。Jazwinski[16]提出了扩展卡尔曼滤波EKF，有效解决了非线性系统问题。Yeh等[10]基于扩展卡尔曼滤波评估气候和水文变化对土壤湿度的影响；Lü等[5]利用扩展卡尔曼滤波同化表层土壤含水量来预测根域土壤含水量。但是，扩展卡尔曼滤波中需要建立的动力模型切线性算子，在强非线性的动力系统中可能并不存在。为了更好地解决非线性问题，Evensen[18]提出了集合卡尔曼滤波(EnKF)，并取得了巨大成功。Brandhorst等[19]通过处理土壤水力参数的不确定性，结合EnKF来预测土壤湿度；Liu等[20]利用EnKF和支持向量机分析土壤湿度模拟不确定性的因素。

然而，在利用集合卡尔曼滤波对土壤湿度进行同化时，也存在一些不足：由于集合成员取样的随机性，使得每次取样的随机误差不同，导致每次的同化结果存在一定差别；另外，每个集合成员赋予相同的权重，降低了影响性较大的集合成员的作用。本文提出了无迹加权集合卡尔曼滤波UWEnKF，增加了较大影响性集合成员的权重，并结合一维Richards方程进行了土壤湿度同化实验，验证其有效性和适用性。

1 方 法

1.1 Richards方程

一维Richards方程描述了非饱和区域土壤水分的运动过程[5,21]，该方程可以表达为：

(1)

式中：θ是土壤水分含量，m3/m3；t是时间，s；z是土壤深度，mm；ψ是土水势，mm；K是土壤水力传导度，mm/s；e是源汇项，即根系吸水后的蒸散发损失率，mm/s。

源汇项的求解参见LÜ等[5]，对方程离散求解过程中，上边界条件取决于下渗率、雨强(当雨强小于下渗率，按雨强下渗，反之按下渗率下渗)；下边界条件采用最底层的土壤水力传导度。方程的时间和空间离散分别采用Crank-Nicolson方法和向后的欧拉方法。

1.2 无迹加权集合卡尔曼滤波UWEnKF

集合卡尔曼滤波EnKF[8,11,18,21-25]于1994年由Evensen[18]提出，之后便被广泛用于提高土壤湿度的模拟精度。该方法包含预测和分析两个阶段，在分析阶段，通过引入噪声扰动，生成状态变量的初始集合，利用模型算子得到每个集合成员的预测值，以及状态变量的误差协方差矩阵；在分析阶段，利用观测值对状态变量的误差协方差矩阵和预测值进行更新，得到状态变量的分析值。许多学者对其理论过程进行了详细的阐述[8,11,18,21-25]。因此，本文不在对EnKF进行描述。

(2)

本文在EnKF的基础上，结合无迹变换[26,27]理论，提出了无迹加权集合卡尔曼滤波UWEnKF，该方法增加了影响性较大的集合成员的权重，且集合成员关于均值对称。

基于公式(2)和随机误差粒子，可以得到UWEnKF中对称性的 个Sigma集合成员：

(3)

UWEnKF的预测和更新过程如下：

预测:

每个集合成员的状态预测过程为：

(4)

更新:

k时刻状态变量X的误差协方差矩阵的计算如下：

PXk|k-1=A(I·×WT)AT

(5)

(6)

(7)

(8)

(9)

(10)

(11)

PYk=S(I·×WT)ST

(12)

(13)

PXYk=A(I·×WT)ST

(14)

Kk=PXYk(PYk+R)-1

(15)

(16)

(17)

(18)

PXk=E(I·×WT)ET

(19)

1.3 评价准则

本研究采用均方根误差RMSE、绝对误差AE和平均绝对误差MAE三个标准对UWEnKF的有效性进行评估，具体公式如下：

(20)

AE=|Xsim,i-Xobs,i|

(21)

(22)

式中：Xsim,i和Xobs,i分别是第i时刻各层土壤湿度的模拟值/同化值和观测值；m为总模拟时刻数。

1.4 研究区和实验数据

流域面积约70 hm2的梅林小流域(31°20′ N,119°51′ E)属于封闭型小流域，地处江苏省宜兴市东南，位于太湖以西约9 km处，如图1[29,30]。流域属湿润气候区域，多年平均降雨量约1 150 mm，年平均气温15.5 ℃，土壤类型主要为红壤和黄棕壤，植被主要为林地、茶树、草地[30]。本文收集了该流域2007年8月1日(第244 d)至9月30日(第273 d)的实时(步长为1 h)气象数据，期间最大日降雨量65 mm，日平均气温23.6 ℃，平均相对湿度86.5%，平均风速8.1 m/s。收集了流域内杨家山水文实验站A站点2007年8月10日至22日、以及8月26日至9月30日的深度分别为5、30、50、100 cm的土壤湿度实时观测数据等，观测步长为1 h，以及土壤孔隙率、残差土壤含水量等参数值。根据Lü等[5]关于梅林流域不同深度的土壤属性数值，本文采用的各层土壤参数见表1。

图1 梅林流域及土壤湿度观测站点[28,29]Fig.1 Meilin study area and soil moisture observation sites

表1 梅林流域各层土壤参数的设定[5]Tab.1 The values of soil parameters at different depth in Meilin area

2 实验结果与分析

在实验设置中，首先利用2007年8月1日至8月31日的实时气象数据对一维Richards方程进行校正，然后，利用9月1日至30日的数据进行土壤湿度数据同化实验。在同化过程中，仅利用EnKF和UWEnKF同化土壤表层(5 cm)的站点观测土壤湿度，进而更新不同深度的土壤湿度模拟值。在同化实验中，随机产生的误差粒子数N设为80；方程模拟了4层不同土壤层的土壤湿度，所以，状态变量的维数n为4，则Sigma集合成员数为2nN+1=641。根据气象数据和土壤湿度数据的观测步长，同化时间步长和模型模拟步长均设为1 h。图2显示了不同深度的土壤湿度同化及模拟值，其中模拟值是基于表1的土壤参数和气象数据等，单独求解Richards方程得到。

由图2可知，不同深度(5、30、50、100 cm)的土壤湿度变化趋势基本一致，都在降雨量较大的第246 d(9月3日)和第262 d(9月19日)有明显增高，其中，30 cm深的土壤湿度变化最为剧烈，100 cm深的土壤湿度变化最为缓慢。校正后的一维Richards方程在5、30、50 cm土壤层的模拟结果与实测变化趋势一致，在强降雨发生时(第246 d和第262 d)，一维Richards方程可以反映土壤湿度的突增变化，但对增长幅度的模拟精度不够。尽管如此，经过一段时间(第248～261 d，第263～273 d)的稳定期(无强降雨)后，土壤湿度的模拟结果与观测趋于一致。在底层(100 cm)，一维Richards方程的模拟结果无法反映土壤湿度的变化过程，尤其是在降雨量较大的第262 d(9月19日)，原因可能是下边界采用最底层的土壤水力传导度，而非实际下边界的土壤水力传导度，导致土壤湿度模拟结果与实际情况存在偏差。

图2还表明，经过同化后(EnKF和UWEnKF)的土壤湿度比模拟值更接近观测值，而UWEnKF方法的同化效果比EnKF方法更好。在上述4种不同深度的土壤层中，UWEnKF同化值比EnKF同化值更趋近实测土壤湿度。当有强降雨输入时(第246 d和第262 d)，虽然Richards方程能够较好地反映土壤湿度的突然升高变化，但是，土壤湿度模拟值与观测值之间相差较大，导致同化值与观测值之间的误差也较大。在第264～271 d的第3层(50 cm)，同化值远离观测值，表明在此期间，同化方法在该层没有改善土壤湿度模拟结果。主要原因可能在于经过9月19日的强降雨后，土壤湿度增量较大，模拟结果虽然反映出土壤湿度变化过程，但是增量幅度较小(第3层)，而EnKF和UWEnKF使得同化结果大于模拟结果，以至于在264～273 dEnKF和UWEnKF同化后的结果在此期间大于模拟结果，但是，UWEnKF同化后的结果与模拟结果越来越接近。在底层(100 cm)，由于Richards方程的模拟结果较稳定，导致同化结果也趋于稳定，但是，整体上UWEnKF同化值比EnKF同化值更接近观测值。综上所述：在土壤湿度同化中，UWEnKF方法比EnKF方法更有效。

图2 2007年9月1日至30日(第273 d)梅林实验区域观测站点A的土壤湿度模拟/同化结果Fig.2 The soil moisture simulations and assimilations from Sep.1 (day 244)to Sep.30 (day 273),2007 at Point A in Meilin watershed

本文计算了一维Richards方程、EnKF方法和UWEnKF方法在不同土壤层的均方根误差RMSE、绝对误差AE和平均绝对误差MAE来评价方法对土壤湿度的模拟/同化效果，图3显示了模拟/同化结果与观测值之间的均方根误差RMSE，图4和图5分别显示了土壤湿度模拟/同化值与观测值之间的绝对误差AE变化过程和平均绝对误差MAE。

由图3可知，EnKF和UWEnKF同化的结果与观测值之间的RMSE在前3层要明显小于模拟结果与观测值之间的RMSE，说明EnKF和UWEnKF可以改善土壤湿度的模拟结果，提高土壤湿度的模拟精度。然而，在底层，EnKF同化后的结果与观测值之间的RMSE与模拟结果与观测值之间的RMSE相差不大，表明EnKF仅同化表层土壤含水量时，在底层改善土壤湿度模拟结果的效果不明显；UWEnKF在该层仍可显著提高土壤湿度的模拟精度。另外，UWEnKF方法的RMSE值在各土壤层均明显小于EnKF方法，表明UWEnKF方法比EnKF方法更好地提高了土壤湿度的模拟精度。这是因为EnKF方法中集合成员等权重，而在UWEnKF方法中，影响性大的集合成员被赋予了更高的权重，这提高了滤波方法的有效性。

由图4可知，在上述4个土壤层，同化值(EnKF和UWEnKF)的AE整体上而言均小于模拟值，且模拟结果与观测值之间的AE越大，基本上同化结果与模拟结果之间的AE也越大，即模型误差越大，同化结果的误差越大。但在第3层(50 cm)，EnKF和UWEnKF同化值的AE在264～273 d要大于模拟值的AE；在第四层(100 cm)，EnKF同化值的AE在第268～273 d高于模拟值的AE(见前文分析)。另外，UWEnKF同化结果的AE整体上小于EnKF同化结果的AE。其中在第2～4层，除了降雨量较大的前后几天，相比EnKF，UWEnKF同化结果的AE明显小于0.02 m3/m3，甚至小于0.01 m3/m3。因此，UWEnKF方法比EnKF方法更好地改善了土壤湿度的模拟结果。

图3 2007年9月1日(第244 d)至30日梅林实验区域观测站点A不同深度的土壤湿度同化/模拟结果与观测值之间的均方根误差RMSEFig.3 The RMSE between the soil moisture assimilations/simulations and observations Sep.1 to Sep.30,2007 at Point A in Meilin watershed

图4 2007年9月1日至30日梅林实验区域观测站点A的土壤湿度同化/模拟结果的绝对误差变化过程Fig.4 The AE changes of soil moisture assimilations/simulations from Sep.1 to Sep.30,2007 at Point A in Meilin watershed

由图5可知，同化和模拟的MAE类似于图3中RMSE的结论，即EnKF和UWEnKF同化的MAE在各层要明显小于模拟结果的MAE，且UWEnKF方法的MAE值在各土壤层均明显小于EnKF方法的MAE。

图5 2007年9月1日至30日梅林实验区域观测站点A不同深度的土壤湿度同化/模拟结果与观测值之间的平均绝对误差MAEFig.5 The MAE between the soil moisture assimilations/simulations and observations Sep.1 to Sep.30,2007 at Point A in Meilin watershed

3 结 论

本文基于无迹转移方程和集合卡尔曼滤波，提出了无迹加权集合卡尔曼滤波UWEnKF。该滤波方法中集合成员的影响性越大，其相应的权重越大，且集合成员的分布关于均值对称。结合一维Richards方程在梅林流域开展了同化表层土壤湿度的实时同化实验。具体的结论如下：

(1)一维Richards方程可以较好地模拟出土壤湿度的动态变化过程，但在底层，该方程的模拟效果不佳，无法较好地反应土壤湿度的变化过程。

(2)在土壤湿度模拟过程中，采用数据同化方法(UWEnKF和EnKF)对表层(5 cm)土壤湿度进行同化，可以提高不同土壤层的土壤湿度模拟精度。同时，同化效果与模型联系紧密，模型误差越大，同化结果的误差越大。

(3)通过土壤湿度模拟的同化实验可知，相较于EnKF方法，无迹加权集合卡尔曼滤波UWEnKF能够更显著地改善土壤湿度的模拟结果，提高土壤湿度的模拟精度。

总之，无迹加权集合卡尔曼滤波UWEnKF比EnKF能够更显著提高土壤水分模拟精度，是一种有效的、实用的数据同化方法。