CN110321917A - 无参数的农田数据同化方法 - Google Patents

Abstract

本发明公开了一种无参数的农田数据同化方法，先预测，利用高斯过程替代物理模型，对已有的观测数据进行训练学习，预测下一时刻的状态向量，观测数据包括土壤含水量；再更新，基于集合卡尔曼滤波思想同化当前时刻的观测，对当前时刻的观测值和高斯过程预测值进行加权且权重根据二者的误差确定，得到当前时刻状态最优估计值并将其作为下一时刻的先验状态向量；最后重复预测和更新，直到所有可用的观测全部被同化。该方法利用高斯过程对时序土壤水观测数据进行顺序性训练并预测，摆脱了物理模型，可以对多种观测信息进行训练学习，可以根据需要预测得到任意分辨率的信息，应用简单，具有普适性。

Description

无参数的农田数据同化方法

技术领域

本发明属于农田水土环境领域，具体涉及一种无参数的农田数据同化方法。

背景技术

非饱和带(非饱和带又称包气带或通气带，指地下水位以上、土壤水分含量小于饱和含水量的土层)对人类活动有着非常重要的意义，涉及农业、林业、水文和工程等多个领域，因此利用非饱和土壤水模型对土壤水运移过程进行模拟十分重要，然而，这种物理模型往往受限于多种不确定性(其来源可能包括输出和胁迫条件数据的误差)以及对于内在机理模型的不完美描述等等。

为了解决上述问题，数据同化方法被提出，它将可用的实际观测融入到所构建的物理模型中，基于两者的误差来确定权重，加权得到状态变量的最优估计值。该方法已在多个领域成功应用，如海洋和大气科学，水文学等。但是，标准的数据同化也存在一些限制，一方面是应用到实际的高维系统时对计算成本要求过高，另一方面是当物理模型不可用或难以构建时，这种传统的基于物理模型的参数型滤波受到极大的挑战。

近些年，为避免构建复杂的物理机理模型，大量的数据驱动(无参数)方法被提出，并被融合到数据同化框架中。其中最为常用的是基于Takens理论的数据同化方法，Hamilton et al.2016提出一种混合的Kalman-Takens方法，利用迟豫坐标向量来重构非线性系统并进行预测，然后再采用Kalman更新同化观，在该方法的基础上，Hamilton etal.2017考虑了动态系统噪声的存在，同时利用自适应滤波过程将动态噪声与观测噪声进行分离，Khaki et al.2018则首次将Kalman-Takens方法应用到实际的陆地水文模型的数据同化中。

发明内容

本发明的目的是提供一种无参数的农田数据同化方法，该方法利用高斯过程对时序土壤水观测数据进行顺序性训练并预测，摆脱了物理模型，可以对多种观测信息进行训练学习，可以根据需要预测得到任意分辨率的信息，应用简单，具有普适性。

本发明采用的技术方案是：

一种无参数的农田数据同化方法，包括步骤：

S1、预测—利用高斯过程(Gaussian Process，简称GP)替代物理模型，对已有的观测数据(如土壤含水量等)进行训练学习，预测下一时刻的状态向量；

S2、更新—基于集合卡尔曼滤波(简称EnKF)思想同化当前时刻的观测，对当前时刻的观测值和高斯过程预测值进行加权且权重根据二者的误差确定，得到当前时刻状态最优估计值并将其作为下一时刻的先验状态向量；

S3、重复S1和S2，直到所有可用的观测全部被同化。

在S1中，首先，给实际观测一定噪声，将实际观测离散为N个符合高斯分布的观测样本；其次，通过以下三种高斯过程建模方法得到t时刻的状态向量预测值：1)利用[1:t-1]时刻的所有观测，构建N个高斯过程模型；2)利用[1:t]时刻的所有观测，构建N个高斯过程模型；3)利用整个模拟时间T内的所有观测，构建N个高斯过程模型；三种高斯过程模型的输入信息均为相关观测信息(如时空信息、胁迫条件等，时空信息为时刻、深度等，胁迫条件为降雨、土壤温度等)，输出为目标变量(如土壤含水量等)；最后，利用训练后的高斯过程模型，预测t时刻的目标状态向量。

S1的具体方法为，

假设当前时刻为t，总模拟时长为T，当前时刻的观测为先验状态向量为三种构建高斯过程的方法为分别利用 和进行训练学习，构建高斯过程模型，得到当前时刻状态向量的预测值其中上标f表示预测。

利用时，

高斯过程模型的训练数据输入项X包括时空信息(时刻、深度等)和胁迫条件(降雨，土壤温度等)，构成维度为n的训练集n＝单次观测数目×观测次数，动力学过程(土壤水运移)高斯过程模型在n维训练集内建立先验分布，然后在n′维测试集下转变为后验分布，n′＝剖面总节点数，则训练样本观测值和测试数据的输出向量之间构成联合高斯分布

其中，μ为Y的均值；

∑＝k(x，x)表示Y的协方差矩阵，其元素∑_i，j＝k(x_i，x_j)；

∑′＝∑^′T为Y与Y′之间的协方差矩阵，∑′_i，j＝k(x_i，x′_j)；

μ′为Y′的先验均值；

∑″为Y′的先验自协方差，∑″_i，j＝k(x′_i，x′_j)；

为观测Y的误差方差；

l为n×n的单位矩阵；

因此得出预测值Y′(即)的后验分布为

后验均值向量为高斯过程所构建的无参数模型的预测均值，Cov(Y′)为相应的方差，那么基于此得到满足高斯分布的后验观测向量

利用时，将换为其余公式完全相同；利用时，将换为其余公式完全相同。

在S2中，假设通过S1得到t时刻的状态向量预测值为当前时刻的观测向量为那么得到当前时刻状态向量的更新值

式中，上标f和a分别表示预测过程和同化过程，i表示样本标号；是基于和观测噪声扰动生成的第i个观测样本；H为观测算子；K_t是t时刻的卡尔曼增益，表示为

式中，R_t为t时刻的观测误差矩阵；为t时刻的状态向量协方差矩阵，表示为

式中，N为样本总数。

在S3中，将t时刻的后验状态向量作为t+1时刻的先验状态向量，然后再次运行S1和S2，重复运行，直到所有观测数据被同化。

本发明的有益效果是：

1.考虑物理模型构建较为复杂，本方法从饱和土壤水问题出发，不再构建土壤水运动方程，利用机器学习方法替代物理模型进行状态预测，本方法不同于Takens理论利用迟豫坐标重构动态过程，而是利用高斯过程对时序土壤水观测数据进行顺序性训练并预测，将其预测值与观测进行加权分析，同时将该方法耦合于EnKF(集合卡曼滤波)框架中构建无参数的农田数据同化方法，本方法基于EnKF(集合卡曼滤波)思想，但摆脱了物理模型，完全依赖于机器学习对未来状态向量进行预测，一方面可以有效地避免由于先验知识不足而造成的难以建立物理模型的困境，另一方面则通过摆脱物理方程的求解而大大提高计算效率。

2.本方法基于Kalman-GP框架，提出三种将高斯过程与集合卡尔曼滤波进行融合的方法，并首次将该框架应用到实际的强烈非线性土壤水问题中。

3.Takens方法要求预测变量和观测的类型需要是一致的，高斯过程可以对多种观测信息(甚至非直接观测)进行训练学习，为非直接观测融入到数据同化框架提供了平台。

4.Takens方法仅能得到与观测在空间上相应的预测值，而高斯过程可以根据需要预测得到任意分辨率的信息，如可以根据某几个深度的含水量预测整个剖面的含水量信息。

5.Takens方法构建迟豫坐标存在较多的随机性，如延迟数、近邻点数和局部替代模型的形式等等很难较为客观地确定，而高斯过程的构建和应用较为简单，更具有普适性。

6.Kalman-Takens方法所提出的动力替代模型是由延迟数决定的局部替代模型，而本发明的Kalman-GP方法则是基于所有先验数据的构建的全局替代模型。

附图说明

图1为试验Falkenberg站点的具体地理位置。

图2为试验站点的日降雨量、不同深度土壤含水量、空气温度和不同深度土壤温度随时间变化图。

图3为采用不同方法预测的土壤含水量与实测值的RMSE(均方根误差)。

具体实施方式

下面结合附图和实施例对本发明作进一步的说明。

采用本方法的试验过程如下。

研究区及数据

本发明实施例的数据来源于全球土壤水网络(International Soil MoistureNetwork)，选取其中的Falkenberg站点作为验证数据，该站点具体位置见图1。如图2所示，数据主要包括不同深度(0.08m，0.15m，0.30m，0.45m，0.60m和0.90m)土壤含水量数据、日降雨数据和空气即土壤温度数据，数据收集于2005年4月25日开始，2015年8月2日结束，共历史100天。

算例设计

本实施例将试验站点简化为高为1m的一维土柱，并将其离散为52个单元(共53个节点)。我们利用不同深度(0.08m，0.15m，0.30m，0.45m，0.60m和0.90m)的含水量观测进行高斯过程(GP)训练建模，得到N＝50个GP模型，对当前时刻含水量剖面进行预测，并采用集合卡尔曼滤波(EnKF)的更新公式对观测进行同化。其中基于三种不同的GP建模方式，将三种EnKF与GP耦合的方法分别表示为：

(1)EnKF-GPA:利用[1:t-1]时刻的所有观测；

(2)EnKF-GPB:利用[1:t]时刻的所有观测；

(3)EnKF-GPC利用整个模拟时间100天内的所有观测。

为证明新提出的EnKF-GP方法的有效性，增加EnKF-Physics(即传统的基于土壤水运动物理模型的EnKF)和GPA(即纯粹的数据驱动方法，不进行数据同化)两个算例。

为量化含水量的预测性能，均方根误差(RMSE)这一指标被引入：

其中，T为总模拟时间，这里T＝100d，E(θ)为同化后的含水量值。

实施结果

如图3所示为采用不同方法时含水量预测的RMSE结果。很明显可以出结论，采用数据驱动构建的替代模型进行数据同化可以达到利用物理方程同化的预测性能，甚至优于其性能，这充分证明了本发明所提出的EnKF-GP方法(包括EnKF-GPA、EnKF-GPB和EnKF-GPC方法)在农田土壤水系统中的有效性。同时，与EnKF-GP方法相比，单纯的GP方法性能不佳，这表明本发明中数据驱动方法与数据同化结合是有必要的。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (5)

1.一种无参数的农田数据同化方法，其特征在于：包括步骤，

S1、预测—利用高斯过程替代物理模型，对已有的观测数据进行训练学习，预测下一时刻的状态向量，观测数据包括土壤含水量；

S2、更新—基于集合卡尔曼滤波思想同化当前时刻的观测，对当前时刻的观测值和高斯过程预测值进行加权且权重根据二者的误差确定，得到当前时刻状态最优估计值并将其作为下一时刻的先验状态向量；

S3、重复S1和S2，直到所有可用的观测全部被同化。

2.如权利要求1所述的无参数的农田数据同化方法，其特征在于：在S1中，首先，给实际观测一定噪声，将实际观测离散为N个符合高斯分布的观测样本；其次，通过以下三种高斯过程建模方法得到t时刻的状态向量预测值，1)利用[1:t-1]时刻的所有观测，构建N个高斯过程模型；2)利用[1:t]时刻的所有观测，构建N个高斯过程模型；3)利用整个模拟时间T内的所有观测，构建N个高斯过程模型；三种高斯过程模型的输入信息均为相关观测信息，输出为目标变量，相关观测信息包括时空信息和胁迫条件，时空信息包括时刻和深度，胁迫条件包括降雨和土壤温度，目标变量包括土壤含水量；最后，利用训练后的高斯过程模型，预测t时刻的目标状态向量。

3.如权利要求2所述的无参数的农田数据同化方法，其特征在于：S1的具体方法为，

假设当前时刻为t，总模拟时长为T，当前时刻的观测为先验状态向量为三种构建高斯过程的方法为分别利用 和进行训练学习，构建高斯过程模型，得到当前时刻状态向量的预测值其中上标f表示预测；

利用时，

高斯过程模型的训练数据输入项X包括时空信息和胁迫条件，构成维度为n的训练集n＝单次观测数目×观测次数，动力学过程高斯过程模型在n维训练集内建立先验分布，然后在n′维测试集下转变为后验分布，n′＝剖面总节点数，则训练样本观测值和测试数据的输出向量之间构成联合高斯分布

其中，μ为Y的均值；

∑＝k(x，x)表示Y的协方差矩阵，其元素∑_i，j＝k(x_i，x_j)；

∑′＝∑′^T为Y与Y′之间的协方差矩阵，∑′_i，j＝k(x_i，x′_j)；

μ′为Y′的先验均值；

∑″为Y′的先验自协方差，∑″_i，j＝k(x′_i，x′_j)；

为观测Y的误差方差；

l为n×n的单位矩阵；

因此得出预测值Y′(即)的后验分布为

后验均值向量为高斯过程所构建的无参数模型的预测均值，Cov(Y′)为相应的方差，那么基于此得到满足高斯分布的后验观测向量

利用时，将换为其余公式完全相同；利用时，将换为其余公式完全相同。

4.如权利要求1所述的无参数的农田数据同化方法，其特征在于：在S2中，假设通过S1得到t时刻的状态向量预测值为当前时刻的观测向量为那么得到当前时刻状态向量的更新值

式中，上标f和a分别表示预测过程和同化过程，i表示样本标号；是基于和观测噪声扰动生成的第i个观测样本；H为观测算子；K_t是t时刻的卡尔曼增益，表示为

式中，R_t为t时刻的观测误差矩阵；为t时刻的状态向量协方差矩阵，表示为

式中，N为样本总数。

5.如权利要求1所述的无参数的农田数据同化方法，其特征在于：在S3中，将t时刻的后验状态向量作为t+1时刻的先验状态向量，然后再次运行S1和S2，重复运行，直到所有观测数据被同化。