一种基于空间差异同化遥感数据与作物模型的区域农作物估产方法与流程

本发明涉及农业作物产量估算技术领域，更具体涉及一种基于空间差异同化遥感数据与作物模型的区域农作物估产方法。

背景技术：

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

现有的作物模型与遥感数据同化技术在实际应用上仍面临很大的局限性。这些同化技术往往要求高精度的遥感反演数据，例如准确的地表叶面积指数(lai)反演数据。然而，现有成熟的遥感反演产品，例如modislai,实际中对lai的反演存在明显误差，直接同化这些lai产品会导致模拟产量严重偏离实际情况。常用的解决办法包括利用地表观测数据修正低分辨率遥感lai数据产品，或使用地表观测数据协助高分辨遥感数据反演发展lai数据产品。这两种解决办法存在一个共同的不足之处就是要求高空间密度及时间连续性的地表观测数据，人力成本及时间成本非常高，这导致此类解决办法难以在大的空间范围及长的时间尺度上进行推广。

也即，现有的作物模型-遥感数据技术的主要缺点是要求大量的地表观测数据，数据收集的人力及时间成本高，而且点位数据在大面积范围的代表性存在较大的不确定性。使得现有作物模型-遥感数据同化技术难以在大的空间范围及长的时间尺度上进行推广应用。

因此，需要新的技术以至少部分解决现有技术中存在的局限。

技术实现要素：

本发明的发明人在实践以及研究中发现，通过对遥感产品所包含的信息进行提炼，剔除误差较大的信息，保留关键的、可以相对准确地反映地表特征的信息，并使用这些信息与作物模型进行同化。该同化技术将有效避免遥感数据误差导致的模拟结果不切实际地偏离实际观测的问题。同时，该方法可以基于成熟的遥感反演产品，不要求使用地表观测数据对遥感数据进行二次修正或单独反演，可以克服数据收集成本带来的同化技术应用能力不足的问题，并且模拟结果与实际产量之间能够实现良好的吻合。

根据本发明的一方面，提供一种基于空间差异同化遥感数据与作物模型的区域农作物估产方法，包括如下步骤：

s1：获取遥感数据，基于遥感数据提取各网格lai数据时序特征，按照作物生长期特征识别作物种植区并提取作物关键物候期数据，包括返青期、抽穗期和成熟期数据；

s2：基于获得的物候期数据以及种植区的样点产量记录数据分别对作物模型的作物生长模块参数以及生物物理过程模块参数进行校准，获得各网格各异的作物生长参数以及各网格相同的生物物理过程参数；

s3：基于校准的参数，运行作物模型进行逐日作物生长模拟；

s4：在具备遥感观测的时间点(称为同化时间点)提取种植区内作物模型模拟的各网格lai值及其对应的遥感反演lai值，分别构建模拟lai矩阵(也称为预报lai矩阵)及遥感lai矩阵(也称为观测lai矩阵)。

s5：对上述两个lai矩阵进行归一化处理以提取空间差异特征，然后对空间差异进行同化运算以获得目标空间差异矩阵，也即分析lai矩阵；

s6：在各网格对作物模型中的lai模拟参数进行分别修正，以使模拟lai空间差异与目标空间差异相符，完成修正后继续运行作物模型至下一同化时间点；以及

s7：重复步骤s4～s6直到设定的同化终点；

s8：将作物模型运行至成熟期，输出产量结果。

根据本发明的一个实施方案，步骤s1中，所述遥感数据选自网格分辨率不大于1km×1km的遥感产品。应该理解的是，本发明的原理也可适用于那些高网格分辨率的遥感产品。

根据本发明的一个实施方案，所述遥感产品为glasslai。

根据本发明的一个实施方案，步骤s2中，所述作物模型为mcwla作物模型，所述模型参数还包括叶面积生长模块参数。

根据本发明的一个实施方案，步骤s5中，所述归一化处理包括将lai矩阵除以其矩阵元素的最大值，使原矩阵归一化为取值范围为(0,1]的归一化矩阵，由下式(1)和(2)表示：

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

所述同化运算为固定增益卡尔曼滤波算法，由下式(3)表示：

其中是同化时间点k的归一化分析lai矩阵，h是先验权重系数；通过归一化函数s^f的反函数进行反归一化转化以得到分析lai矩阵由下式(4)表示：

根据本发明的一个实施方案，步骤s6中，所述lai模拟参数包括返青期最大叶面积指数laimax,r1，抽穗期最大叶面积指数laimax以及“生殖生长期叶面积衰老系数laidg。

根据本发明的一个实施方案，步骤s6中，对lai模拟参数修正的修正包括：

将mcwla模型中对lai的模拟简化为下式(5):

式中是lai在同化时间点k网格i的lai模拟值，m代表mcwla模型中用来模拟lai的方程组，代表在同化时间点k-1到同化时间点k期间模型采用的lai生长参数(叶面积生长参数)，称为预报参数；代表同化时间点k-1网格i的分析lai；

在同化时间点k，当计算得到分析lai，也即后，将基于进行反向优化，由下式(6)和(7)表示：

其中，w是同化时间点k上的模型模拟lai值与分析lai值之间的误差绝对值；δp代表对的修正量，通过调整δp使得w最小化；如果模拟lai值与分析lai值的误差在0.05以下，则判断为达到优化目的，此时优化的调整量被标记为δpo；经过δpo修正之后得到同化时间点k网格i的分析参数

将该分析参数应用于同化时间点k-1到同化时间点k期间对网格i的lai模拟；此时，在同化后的模型中，式(5)转化为下式(8)：

同时，也作为下一个同化时段内，也即时间点k到时间点k+1期间的预报参数以模拟时间点k+1时的lai模拟值

根据本发明的一个实施方案，步骤s8包括将网格尺度的模拟产量按照行政边界汇总，输出行政单元尺度的产量模拟结果。

根据本发明的一个实施方案，所述作物为小麦。

本发明避免了同化技术对于地表观测数据的需求，仅使用中低分辨率遥感数据与作物模型同化，即可在大面积、长时间开展较高精度的作物估产。这一新的作物模型-遥感数据同化技术不仅能够提高大面积作物估产的精度，而且极大地降低了应用成本。

附图说明

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

图1是根据本发明一个实施方案的基于空间差异同化遥感数据与作物模型的区域农作物估产方法中的同化流程示意图。

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

图3为图2所示研究区中2001-2008年数据同化前后平均模拟产量与实际产量空间分布对比图。

具体实施方式

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

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

本发明提出一种新的基于lai时空格局特征的同化方案，具体实施步骤如下：

首先，基于遥感数据提取各网格lai数据时序特征，按照作物生长期特征识别作物种植区并提取作物关键物候期数据，包括返青期、抽穗期、和成熟期数据。所述遥感数据可以来源于成熟的遥感产品，特别是那些中低网格分辨率的产品，例如网格分辨率不大于1km×1km的那些。例如网格分辨率为1km×1km的glasslai遥感产品。当然，本发明亦可应用于那些高网格分辨率的产品。

其次，基于物候期数据对作物模型的作物生长模块参数进行校准，获得各网格各异的作物生长参数，并且基于种植区的样点产量记录数据对作物模型生物物理过程模块参数进行校准，获得各网格相同的生物物理过程参数。例如，所述作物模型可以为mcwla作物模型，应该理解是，本发明的原理也可应用于其他适当的作物模型。模型模拟在网格尺度运行，即对各网格分别模拟作物生长。模型参数可以分为三类，包括作物生长参数，用于控制物候模拟；叶面积生长参数，用于控制叶面积生长模拟；生物物理过程参数，用于控制水平衡、光合作用、呼吸作用、干物质累积等过程。这些具体参数本身为本领域所熟知，例如生物物理过程参数可以包括与土壤水分平衡、光合、呼吸作用以及产量形成等相关的参数。

接着，基于校准后的参数，运行该作物模型进行逐日作物生长模拟。所述模型运行的初始参数为基于研究区内(也即上述识别的作物种植区)典型种植区校准的参数，此套参数作为初始参数应用于所有研究区内网格的模拟初始化。

图1是根据本发明一个实施方案的基于空间差异同化遥感数据与作物模型的区域农作物估产方法中的同化流程示意图。下面结合附图对本发明的同化步骤进行详细的说明，图中以mcwla作物模型为例。

如图所示，在具备遥感观测的时间点(称为同化时间点)提取研究区内作物模型模拟的各网格lai值及其对应的遥感反演lai值，分别构建模拟lai矩阵(也可称为预报lai矩阵，也即模拟所得的各lai值的集合)及遥感lai矩阵(也可称为观测lai矩阵，也即遥感反演所得的各lai值的集合)。图中示出了k+1个同化时间点，其中k可以为大1的自然数，可以根据具体情况在作物的生长期内进行适当的选择，图中示出了在播种期至成熟期之间。

对两个lai矩阵进行归一化处理以提取空间差异特征，对空间差异进行同化运算以获得目标空间差异矩阵(称为分析lai矩阵)。

所述归一化方法为将lai矩阵除以其矩阵元素的最大值，使原矩阵归一化为取值范围为(0,1]的归一化矩阵：

其中，表示时间节点k上的lai预报矩阵，dk代表由同一时间点上对应的观测矩阵。和分别表示由和dk转化而来的归一化的预报lai矩阵和观测lai矩阵。s^f和s^d是用来对原始矩阵进行归一化的归一化函数。

其中，s5所述同化方法为固定增益卡尔曼滤波算法：

其中是同化时间点k的归一化分析矩阵，和分别表示同时间的归一化模拟lai矩阵与归一化遥感反演lai矩阵。h是一个固定的先验权重系数,取值范围可以为(0,1)。通过归一化函数s^f的反函数(记做)进行反归一化转化以得到分析lai矩阵如下式(4)所示。

之后，在各网格对模型中的lai生长参数进行分别修正以使模拟lai空间差异与目标空间差异相符。

更具体地，所述lai生长参数包括返青期最大叶面积指数laimax,r1，抽穗期最大叶面积指数laimax以及“生殖生长期叶面积衰老系数laidg等。

所述lai参数的修正并非更改模型输入参数的初始状态，而是对当前时间点上的参数值进行修正，该参数值仅应用于当前同化时间点到下一同化时间点之间的模型模拟。对参数的修正基于以下步骤：

将mcwla模型中对lai的模拟简化为如下式(5):

式中是lai在同化时间点k网格i的lai模拟值，m代表mcwla模型中用来模拟lai的方程组，代表在同化时间点k-1到同化时间点k期间模型采用的lai生长参数，称为预报参数。代表同化时间点k-1网格i的分析lai。在同化时间点k，当计算得到分析lai之后，将基于进行反向优化，如下式(6)和(7)：

其中，w是同化时间点k模型模拟lai值与分析lai值之间的误差绝对值。δp代表对的修正量。对的反向修正的目的是通过调整δp使得w最小化。如果模拟lai值与分析lai值的误差在0.05以下，则认为达到优化目的，满足精度的要求，此时优化的调整量被标记为δpo。经过δpo修正之后可以得到被称作同化时间点k网格i的分析参数。该参数将应用于同化时间点k-1到同化时间点k期间对网格i的lai模拟。此时，在同化后的模型中，式(5)可以被转化为下式(8)：

同时，也作为下一个同化时段内(时间点k到时间点k+1期间)的预报参数模拟时间点k+1时的lai值

也即，在完成修正后继续运行模型至下一同化时间点。如此循环运行，直到设定的同化终点。

最后，将模型运行至成熟期，输出模拟产量结果。例如可以将网格尺度的模拟产量按照行政边界汇总，输出行政单元尺度的产量模拟结果。

本发明的方法可以适用于多种作物例如小麦、大豆等等。

实施例

本案例以估测华北平原地区冬小麦产量为例进一步阐述本发明的技术方案。包括以下步骤：

选择华北平原中部作为研究区域，研究时间为2001-2008年。基于1km×1kmglasslai遥感数据，提取各网格lai数据时序特征，按照作物生长期特征识别作物种植区并提取作物关键物候期数据，研究区如图2所示。

基于物候期数据对mcwla作物模型的作物生长模块参数进行校准，获得各网格各异的作物生长参数。获取样点县的记录产量数据，并依据产量数据对生物物理过程参数进行校准。校准所得生物物理过程参数应用到整个研究区的模拟中。对于叶面积生长参数(“返青期最大叶面积指数laimax,r1”，“抽穗期最大叶面积指数laimax”以及“生殖生长期叶面积衰老系数laidg”)采用其模型默认值。

基于本发明提出的技术方案的叶面积同化实施过程如下：

从播种时期开始使用mcwla模型进行模拟，在具备遥感观测的时间点(称为同化时间点)提取研究区内作物模型模拟的各网格lai值及其对应的遥感反演lai值，分别构建模拟lai矩阵(也可称为预报lai矩阵)及遥感lai矩阵(也可称为观测lai矩阵)。

对两个lai矩阵进行归一化处理以提取空间差异特征，对空间差异进行同化运算以获得目标空间差异矩阵(称为分析lai矩阵)。

所述归一化方法为将lai矩阵除以其矩阵元素的最大值，使原矩阵归一化为取值范围为(0,1]的归一化矩阵：

其中，表示时间节点k上的lai预报矩阵，dk代表由同一时间点上对应的观测矩阵。和分别表示由和dk转化而来的归一化的预报lai矩阵和观测lai矩阵。s^f和s^d是用来对原始矩阵进行归一化的归一化函数。

其中，s5所述同化方法为固定增益卡尔曼滤波算法：

其中是同化时间点k的归一化分析矩阵，和分别表示同时间的归一化模拟lai矩阵与归一化遥感反演lai矩阵。h是一个固定的先验权重系数，本案例中取值0.53。通过归一化函数s^f的反函数(记做)进行反归一化转化以得到分析lai矩阵如下式(4)所示。

之后，在各网格对上述提及的模型中的lai生长参数进行分别修正以使模拟lai空间差异与目标空间差异相符。对参数的修正基于以下步骤：

将mcwla模型中对lai的模拟简化为如下式(5):

式中是lai在同化时间点k网格i的lai模拟值，m代表mcwla模型中用来模拟lai的方程组，代表在同化时间点k-1到同化时间点k期间模型采用的lai生长参数，称为预报参数。代表同化时间点k-1网格i的分析lai。在同化时间点k，当计算得到分析lai之后，将基于进行反向优化，如下式(6)和(7)：

其中，w是同化时间点k模型模拟lai值与分析lai值之间的误差绝对值。δp代表对的修正量。对的反向修正的目的是通过调整δp使得w最小化。如果模拟lai值与分析lai值的误差在0.05以下，则认为达到优化目的，满足精度的要求，此时优化的调整量被标记为δpo。经过δpo修正之后可以得到被称作同化时间点k网格i的分析参数。该参数将应用于同化时间点k-1到同化时间点k期间对网格i的lai模拟。此时，在同化后的模型中，式(5)可以被转化为下式(8)：

同时，也作为下一个同化时段内(时间点k到时间点k+1期间)的预报参数模拟时间点k+1时的lai值

也即，在完成修正后继续运行模型至下一同化时间点。如此循环运行，直到设定的同化终点。在本案例中，同化步长同glasslai数据步长，为8天。同化起点设置为返青期，终点设置为灌浆期。同化运行至灌浆期后结束，模型继续运行至成熟期，输出作物产量并聚合至县级尺度与产量记录数据进行比较。

本发明实施例所述的估产方法，融合了遥感数据和作物模型的优势，将lai的空间差异作为同化变量，避免了遥感数据中lai偏低的不利影响，实现了对区域产量估算精度的提升。与未同化相比，2001-2008年均方根误差(rmse)平均减小了33.1％,最大减小50.2％。r²平均提升0.13，最高提升0.27。估产精度明显提升，区域空间分布趋势也更为符合现实。图3为2001-2008年研究区内同化前后平均模拟产量与实际产量空间分布对比图。

本发明融合了遥感数据和作物模型的优势，利用成熟的lai遥感反演产品以及mcwla作物模型进行同化估产。使用空间差异同化的方法对模型模拟的lai进行了调整，实现了遥感数据与作物模型的同化。同化后的模拟产量与未同化相比较，均方根误差(rmse)减小而决定系数(r²)明显上升，作物模型产量估算的精度有显著提高，产量空间分布趋势与统计产量一致。同时，本发明不要求地表观测数据，不要求对遥感产品进行二次修正，显著降低了同化技术的成本，提升了同化技术的实际应用能力。

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