一种基于集合同化策略的区域农作物估产方法

摘要

本发明公开了一种基于集合同化策略的区域农作物估产方法，包括：获取研究区域内预测年Y以及之前M年的遥感LAI数据；针对第Y‑1年，进行差异同化，获得多个格网尺度模拟产量，以及相应的多组相关系数r以及均方根误差rmse；根据r以及rmse计算拟合指数F，共获得多个拟合指数；选择最大的拟合指数，作为最优结果，并记录对应的同化权重，作为所述第y‑1年的最优同化权重H

说明书

技术领域

本发明涉及农业遥感技术领域，具体而言涉及一种基于集合同化策略的区域农作物估产方法。

背景技术

作物模型在进行区域尺度的作物产量估算时常常面临输入数据不足的问题。地表特征、近地表环境以及作物管理措施往往存在明显的空间差异，而一些必要的模型输入数据，例如作物物候数据，气象数据，作物管理信息等均在站点尺度记录，这导致作物模型在应用到区域尺度时难以获取足够的数据以表现初始条件、作物参数、生长过程等关键因素的空间异质性。卫星遥感数据提供了大范围地表信息的连续监测数据，可以反映地表信息的空间连续性以及时序变化特征。这一优势有效地补充了作物模型在区域尺度上应用中的弱点，因此作物模型与遥感数据同化技术成为了当前改进作物模型区域模拟精度的重要途径，也是近年来的作物估产研究领域的热点方向之一。

在作物模型与遥感数据同化的过程中，一个对同化结果有关键影响的因子是遥感数据的不确定性。该不确定性的主要作用是量化模型模拟值与遥感数据的同化权重。同化权重的合理与否直接影响最终的同化结果精度，一旦同化权重设置出现偏差，则同化结果也会产生显著偏移，严重影响同化技术的应用效果。量化遥感数据不确定性的最佳方法是将遥感数据与地表观测结果进行对比分析，然而该方法要求时空密度的地表观测数据，人力成本及时间成本非常高，难以推广。现有的作物模型与遥感数据同化技术往往采用简化的方法，例如直接设置同化权重值，来降低对地表观测的需求。这种方法取得了一些积极的研究结果，但对于实际的同化估产应用帮助有限。其原因是在实际应用中，同化权重需要作为先验参数存在以驱动同化过程的进行，而现有研究并未涉及如何对同化权重进行准确先验。事实上，由于同化权重的影响因素多，年际变化复杂，因此通过历史经验无法直接推导出最优同化权重，这使得同化技术的实用性大打折扣。目前，最优同化权重无法先验获知是同化估产技术走向实际应用的一个主要限制，如何突破这一限制，在未知最优同化权重的前提下开展同化估产并获得足够精度的估产结果目前仍是技术空白，亟需突破。

因此，需要新技术来至少部分解决现有技术中上述局限性。

发明内容

本发明提出了基于集合同化策略的同化估产方案，该方案的估产精度与基于最优同化权重的估产精度相近，即在最优同化权重未知的情况下亦可获得足够准确的估产结果，极大地提升了同化估产技术的实用性。

根据本发明一方面，提供了一种基于集合同化策略的区域农作物估产方法，包括如下步骤：

S1、获取研究区域内预测年Y以及之前M年的遥感LAI数据，M为大于3的自然数；

S2、针对第Y-1年，利用作物模型以及所述遥感LAI数据进行空间差异同化运算，所述空间差异同化运算方程为：

其中是同化时间点k的归一化分析矩阵，和分别表示同化时间点k的归一化模拟LAI矩阵与归一化遥感反演LAI矩阵，H是同化权重，取值范围为0.01-1.00；

其中，在所述0.01-1.00的取值范围内，H以步长S分别取值，并在所得的1/S个不同的同化权重值情况下，分别进行同化估产运算，由此获得相应的1/S个格网尺度模拟产量；

S3，基于研究区域内相应的产量记录数据，对每一个所述格网尺度模拟产量进行精度评估，用相关系数r以及均方根误差rmse表征，共获得1/S组r以及rmse；

S4、根据S3中所得r以及rmse计算拟合指数F，共获得1/S个拟合指数；

所述拟合指数F的计算公式为：

其中，Fi是第i次同化估产的拟合指数，r_i和rmse_i代表第i次同化估产的相关系数以及均方根误差，r_max,r_min,rmse_max以及rmse_min分别代表所述1/S组r以及rmse中r的最大值、最小值以及rmse的最大值、最小值；

S5、选择最大的拟合指数，作为最优结果，并记录对应的同化权重，作为所述第y-1年的最优同化权重H_{op_y-1}；

S6、对y-2，y-3......y-M年重复步骤S2～S5,获得相应各年的最优同化权重H_{op_y-2}，H_{op_y-3}......H_{op_y-M}；

S7、基于H_{op_y-1}，H_{op_y-2}，H_{op_y-3}......H_{op_y-M}中的最大值H_opmax以及最小值H_opmin划定扩张同化权重取值范围，扩张同化权重取值范围的下限称为H_low，上限称为H_up，H_low为所述最小值向外扩大10％所得，H_up为所述最大值向外扩大10％所得，且扩张同化权重取值范围不大于0.99，不小于0.01；

S8、基于S7中所得扩张同化权重取值区间H_low-H_up，针对所述预测年Y进行模拟估产，包括：利用扩张同化权重取值区间H_low-H_up替代取值范围0.01-1.00的情况下，重复步骤S2，展开空间差异同化运算，共进行(H_up-H_low)/S+1次同化运算，由此获得相应的(H_up-H_low)/S+1个格网尺度模拟产量，取平均值作为所述预测年Y的模拟产量。

根据本发明的一个实施方案，S1步骤中，M为4，当然也可以适当选择其他的数值。

根据本发明的一个实施方案，所述作物模型为MCWLA系列模型，或者也可以选择其他的作物模型。

根据本发明的一个实施方案，步长S为0.01，或者根据具体情况选择其他的数值，例如0.02等。

根据本发明的一个实施方案，遥感LAI数据可选自Copernicus LAI，GLASS LAI以及GLOBMAP LAI，或其他遥感LAI产品。

根据本发明的一个实施方案，所述作物为小麦，或者其他的作物例如玉米。

根据本发明的一个实施方案，步骤S2中，所述归一化分析矩阵、归一化模拟LAI矩阵和归一化遥感反演LAI矩阵通过归一化方法获得，使原矩阵归一化为取值范围为(0,1]的归一化矩阵，由下式(3)和(4)表示：

其中，表示同化时间点k上的模拟LAI矩阵，D_k代表同一时间点k上对应的遥感LAI矩阵；和分别表示由和D_k转化而来的归一化模拟LAI矩阵和归一化遥感LAI矩阵；S^f和S^D是用来对原矩阵进行归一化的归一化函数。

根据本发明的一个实施方案，步骤S2中，所述差异同化还包括：通过归一化函数S^f的反函数进行反归一化转化以得到分析LAI矩阵并由驱动模型向下一同化时间点k+1运行：

与现有技术相比，本发明在遥感数据与作物模型同化技术的基础上，融入了集合运算思想。为“无法先验获知最优同化权重的情况下如何开展同化估产并获得足够精度的估产结果”这一实际问题提供了解决策略。运用本发明提出的基于集合同化策略的同化估产方案，估产精度与基于最优同化权重的估产精度相近，即在最优同化权重未知的情况下亦可获得足够准确的估产结果。本发明弥补了同化估产技术面向业务化应用时的不足，极大地提升了同化估产技术的实际应用能力。

附图说明

附图中相同的附图标记标示了相同或类似的部件或部分。本发明的目标及特征考虑到如下结合附图的描述将更加明显，附图中：

图1是根据本发明一个实施方案的基于集合同化策略的区域农作物估产方法的流程示意图。

图2是根据本发明一个实施方案的实施本发明方法的研究区的示意图；

图3为根据本发明方法在应用不同遥感LAI产品时的实施效果图。

具体实施方式

为清楚地说明本发明中的方案，下面给出优选的实施例并结合附图详细说明。以下的说明本质上仅仅是示例性的而并不是为了限制本公开的应用或用途

应该理解的是，本发明所引用的作物模型和遥感数据本身是已知的，例如模型的各个子模块、各种参数、运行机制等等，因此本发明重点阐述作物模型和遥感数据之间集合同化过程。

图1是根据本发明的一个实施方案的示意图。参考图1，根据本发明的基于集合同化策略的区域农作物估产方法，包括如下步骤：

S1、获取研究区域内预测年Y以及之前M年的遥感LAI数据，M为大于3的自然数；图中示出了M为4，也即获取第y-1,y-2,y-3,y-4的遥感LAI数据，其中遥感LAI数据可选自Copernicus LAI，GLASS LAI以及GLOBMAP LAI，或者其他适当的数据。也即，本方法针可应用于多种遥感LAI数据，同化估产运算的空间分辨率及时间频率与遥感LAI数据的时空分辨率保持一致即可。

S2、针对第Y-1年，利用作物模型(图中示出了作物模型为MCWLA-Wheat模型)以及所述遥感LAI数据进行差异同化，包括进行空间差异同化运算，所述空间差异同化运算方程为：

其中是同化时间点k的归一化分析矩阵，和分别表示同化时间点k的归一化模拟LAI矩阵与归一化遥感反演LAI矩阵，H是同化权重，取值范围为0.01-1.00；应该理解，式中的上标是为了区分矩阵，a是analysis缩写，表示分析矩阵；f是forecast缩写，表示预报矩阵即模型模拟矩阵；n是normalization缩写，表示归一化。

其中，在所述0.01-1.00的取值范围内，H以步长S(例如0.01)分别取值，并在所得的1/S(例如100)个不同的同化权重值情况下，分别进行同化估产运算，由此获得相应的1/S(100)个格网尺度模拟产量。

更具体地，差异同化的一次同化过程包括：

1)在存在遥感数据的时间点(称为同化时间点)对模拟的LAI矩阵以及对应的遥感LAI矩阵进行归一化运算使原矩阵归一化为取值范围为(0,1]的归一化矩阵(下式(3)和(4))，归一化运算方法较多，例如可以通过将LAI矩阵除以其矩阵元素的最大值来获得：

其中，表示同化时间点k上的模拟LAI矩阵，D_k代表同一时间点k上对应的遥感LAI矩阵；和分别表示由和D_k转化而来的归一化模拟LAI矩阵和归一化遥感LAI矩阵；S^f和S^D是用来对原矩阵进行归一化的归一化函数。

2)使用同化方程(1)对归一化模拟LAI矩阵及归一化遥感LAI矩阵进行同化：

其中是同化时间点k的归一化分析矩阵，H是同化权重，取值范围为0.01-1.00，在此H以步长0.01分别取值，得到100个不同的同化权重值。

3)通过归一化函数S^f的反函数进行反归一化转化以得到分析LAI矩阵并由驱动模型向下一同化时间点k+1运行：

其中，该步骤中一次完整同化估产运算过程包括：模型从播种期开始逐日模拟，至同化起点开始与遥感数据进行同化，同化过程持续至同化终点，模型继续运行至成熟期输出格网尺度模拟产量。

S3，基于研究区域内相应的产量记录数据，对每一个所述格网尺度模拟产量进行精度评估，用相关系数r以及均方根误差rmse表征，共获得1/S组r以及rmse；所述估产精度评估基于实际统计产量的空间尺度，同化运算所得估产结果为格网尺度，需聚合至统计产量的尺度，例如县级。

S4、根据S3中所得r以及rmse计算拟合指数F，共获得1/S个拟合指数；

所述拟合指数F的计算公式为：

其中，F_i是第i次同化估产的拟合指数，r_i和rmse_i代表第i次同化估产的相关系数以及均方根误差，r_max,r_min,rmse_max以及rmse_min分别代表所述1/S组r以及rmse中r的最大值、最小值以及rmse的最大值、最小值；F值越大代表估产结果精度越高。

S5、选择最大的拟合指数，作为最优结果，并记录对应的同化权重，作为所述第y-1年的最优同化权重H_{op_y-1}；

S6、对y-2，y-3......y-M年重复步骤S2～S5,获得相应各年的最优同化权重H_{op_y-2}，H_{op_y-3}......H_{op_y-M}，图中示出了M为4，也即H_{op_y-4}；

S7、基于H_{op_y-1}，H_{op_y-2}，H_{op_y-3}......H_{op_y-M}中的最大值H_opmax以及最小值H_opmin划定扩张同化权重取值范围，扩张同化权重取值范围的下限称为H_low，上限称为H_up，H_low为所述最小值向外扩大10％所得，H_up为所述最大值向外扩大10％所得，且扩张同化权重取值范围不大于0.99，不小于0.01；也即，

H_low＝max(0.9*H_opmin,0.01)>

H_up＝min(1.1*H_opmax,0.99)>

S8、基于S7中所得扩张同化权重取值区间H_low-H_up，针对所述预测年Y进行模拟估产，包括：利用扩张同化权重取值区间H_low-H_up替代取值范围0.01-1.00的情况下，重复步骤S2，展开空间差异同化运算，共进行(H_up-H_low)/S+1(101)次同化运算，由此获得相应的(H_up-H_low)/S+1个格网尺度模拟产量(也即附图中n组结果)，取平均值作为所述预测年Y的模拟产量。

下面结合具体实施例，来进一步阐述本发明方法：

下面以华北平原地区冬小麦估产为例，示例性说明本发明的方法的具体应用。

选择华北平原中部作为研究区域，研究时间为2008-2015年。采用MCWLA-Wheat模型，使用的遥感LAI数据包括Copernicus LAI(空间分辨率1km×1km，时间频率为每10天一期),GLASS LAI(空间分辨率1km×1km，时间频率为每8天一期),以及GLOBMAP LAI(空间分辨率0.08°×0.08°，时间频率为每8天一期)，每种遥感LAI产品分别与MCWLA-Wheat模型进行同化，对冬小麦产量进行估算。研究区如图2所示。

基于本发明提出的集合同化策略，同化估产实施过程如下(以Copernicus LAI为例，GLASS LAI以及GLOBMAP LAI的实施过程与Copernicus LAI一致)：

在本例中，将2012-2015年作为估产时期。针对每一年的估产，以其前4年作为历史时期进行同化估产以提供关于最优同化权重的先验知识。即针对年y的同化估产运算，需要对其前4年(即y-1,y-2,y-3,y-4)先进行同化估产运算以对年y提供关于最优同化权重的先验知识。例如对于2012年的同化估产，需要对2008～2011年先行采用同化估产运算以获取逐年的最优同化权重。同化采用空间差异同化算法进行。

下面结合附图1对本发明的操作步骤进行详细说明。

如图所示，在年y-1，使用MCWLA-Wheat模型以及遥感LAI数据进行空间差异同化运算，同化权重H取值0.01～1.00，步长为0.01，即H有100个取值，由此可以驱动100次同化估产运算，得到100组估产结果。每次同化估产运算的完整过程如下：

从冬小麦播种期开始使用MCWLA-Wheat对冬小麦生长进行逐日模拟。同化起点为冬小麦返青期，同化终点为冬小麦成熟期。在返青期至成熟期之间，在每个具备遥感观测的时间点(称为同化时间点)使用空间差异同化算法对模型模拟LAI以及遥感LAI数据进行同化。在成熟期结束同化及模型模拟过程，逐格网输出冬小麦模拟产量。

一次同化估产运算完成后，将格网尺度的估产结果聚合至县级尺度并与统计产量进行对比，以相关系数r以及均方根误差rmse表征估产精度。使用不同的同化权重H重复100次同化运算，可得到100组r以及rmse，每个同化权重H均对应一组r及rmse值。在此基础上，计算每次同化估产结果的拟合指数F，其计算公式为上述公式(2)：其中，r_max,r_min,rmse_max以及rmse_min分别代表100组r以及rmse中r的最大值、最小值以及rmse的最大值、最小值。

由公式(2)计算可得到对应100个同化权重H的100个拟合指数F。选择最大的拟合指数作为最优拟合指数，并将其对应的同化权重H作为年y-1的最优同化权重，命名为H_{op_y-1}。

针对y-1年计算Hop_y-1的过程需在y-2，y-3,以及y-4三年分别重复，获得逐年的最优同化权重H_{op_y-2}，H_{op_y-3}，H_{op_y-4}。

_{op_y-1}，H_{op_y-2}，H_{op_y-3}，H_{op_y-4}中的最大值H_opmax以及最小值H_opmin划定同化权重取值区间，取值区间设定为最大、最小值区间向外扩大10％形成，且取值范围不大于0.99，不小于0.01。区间下限称为H_low，上限称为H_up，具体参见上述公式(6)和(7)。

完成上述计算后，对年y进行集合同化估产运算，运算所采用的同化权重H为基于得到的H_low及H_up所围成的取值范围[H_low,H_up]，取值以0.01步长递增。共运行(H_up-H_low)/0.01+1次同化运算。集合运算所得估产结果的均值作为最终估产结果，模拟结果参见附图3。

针对不同的遥感LAI数据，估产方式与上述过程一致，区别在于使用不同的遥感LAI数据时，模型运行的格网尺度的空间分辨率需与遥感LAI数据保持一致，同化时间点间隔需与遥感LAI数据的时间频率保持一致。另外，使用不同遥感LAI数据时，每一年所得的最优同化权重会发生变化。

为了验证本发明的技术效果，采用如下方法来评价本发明的实施效果：

首先，将应用于年y-1～y-4的最优同化权重计算方法同样应用于年y，获得年y的如下参数：包括r_max，r_min，rmse_max，rmse_min以及最大拟合指数F，称为F_op。接着，将年y的集合同化估产结果与实际统计产量进行对比，计算得到其r及rmse，并结合前述r_max，r_min，rmse_max，rmse_min计算集合同化估产结果的F值，称为F_s。根据F_op及F_s，计算准确度指标pc：

准确度指标pc用于表征集合同化的估产准确程度相比应用最优同化权重的估产准确度的百分比，pc越大表示准确度越高。

此外，为了对比估产准确度，可以计算一些典型的估产方法的准确度指标pc，这些方法包括：

a.使用年y-1的最优同化权重H驱动年y的同化估产运算；

b.使用年y-1～年y-2的最优同化权重H的均值驱动年y的同化估产运算；

c.使用年y-1～年y-3的最优同化权重H的均值驱动年y的同化估产运算；

d.使用年y-1～年y-4的最优同化权重H的均值驱动年y的同化估产运算；

e.使用年y-1～年y-4的最优同化权重H驱动年y的同化估产运算并以结果均值作为估产结果；

f.使用年y-1～年y-4的最优同化权重H的最大最小值区间,即[H_opmin,H_opmax]，步长0.01，驱动年y的同化估产运算并以结果均值作为估产结果；

g.使用本发明所提出方法对年y进行估产。

使用上述方法a～g在2012-2015年分别进行同化估产运算，四年估产结果的准确度指标pc的均值被用来衡量该方法的综合实施效果。结果如图3所示，显示在应用不同遥感LAI产品时，本发明所提出的集合同化策略均表现出优秀的实施效果。例如对于CopernicusLAI，估产准确度可以达到使用最优同化权重时估产准确度的98.1％，使用GLASS LAI以及GLOBMAP LAI时，估产准确度分别达到了使用最优同化权重估产准确度的90.8％以及85.0％。使用分辨率较高的遥感产品可以获得相对更好的准确度，更接近使用最优同化权重时的估产效果。

此外，本发明提出的集合同化策略的估产准确度在应用不同遥感LAI数据进行的同化估产中稳定优于其他估产方法，具有普适性。上述结果表明，本发明所提出的集合同化策略可以在未知最优同化权重的情况下开展同化估产并获得足够精度的估产结果。本发明弥补了同化估产技术面向业务化应用时的不足，极大地提升了同化估产技术的实际应用能力。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的装置及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。