CN104134095A - 一种基于尺度转换和数据同化的农作物产量估测方法 - Google Patents

Abstract

本发明属于农业遥感领域，具体涉及一种基于尺度转换和数据同化的农作物产量估测方法及其指导农作物生产中的应用，具体步骤如下：S1收集参数完成WOFOST作物模型的空间化；S2获取作物类型分布图和待测农作物纯度百分比图；S3获得区域30米的TM LAI；S4构建二级尺度转换模型，生成时间序列尺度调整LAI；S5得待测农作物关键物候期的遥感观测误差和作物模型误差；S6构建四维变分代价函数，获得优化作物模型参数；S7输出县域待测农作物单产。本发明克服了遥感观测像元与作物模型模拟单元之间尺度不匹配的问题，提高了数据同化模型的精度，适合于县域尺度的农作物产量估测，尤其适合于县域尺度的冬小麦产量估测。

Description

一种基于尺度转换和数据同化的农作物产量估测方法

技术领域

本发明属于农业遥感领域，具体涉及一种基于尺度转换和数据同化的农作物产量估测方法。

背景技术

利用数据同化技术把遥感反演参数融合到作物机理过程模型，是当前改进区域作物生长模拟精度和提高作物估产精度的重要途径。在田间尺度上，基于作物光合、呼吸、蒸腾、营养等机理过程的作物生长模型依靠其内在的物理过程和动力学机制，可以准确模拟作物对象以“天”为时间步长的连续演进和单点尺度上作物的生长发育状况及产量。而当作物模型扩展应用到区域尺度时，由于地表、近地表环境非均匀性，导致了作物模型中的初始条件、土壤参数、作物参数、气象强迫因子的空间分布的不确定性和资料获取的困难。卫星遥感具有空间连续和时间动态变化的优势，能够有效解决区域作物参数获取困难这一瓶颈。然而遥感对地观测由于受卫星时空分辨率等因素的制约，还不能真正揭示作物生长发育和产量形成的内在过程机理、个体生长发育状况及其与环境气象条件的关系，而这正是作物模型的优势所在。数据同化技术通过耦合遥感观测和作物模型，能够实现两者的优势互补，对提高区域作物产量估测具有巨大的潜力。

由于高时间分辨率遥感数据能扑捉作物生长发育变化以及物候信息，因此对农作物长势监测与估产具有非常重要的意义，然而，目前在轨运行的高时间分辨率遥感卫星数据往往空间分辨率都很低，例如AVHRR，MODIS，MERIS和SPOT vegetation空间分辨率为250m-1km尺度。这种低空间分辨率遥感数据，进一步增大了像元异质性，造成了遥感观测尺度与作物模型模拟尺度之间的不匹配，极大限制了数据同化模型的精度。另外，标准的MODIS LAI产品是面向全球尺度所有的植被类型开发，并不是针对农作物设计开发的。研究表明，MODIS LAI在冬小麦区域存在严重的低估现象。因此，直接同化标准MODIS LAI会导致同化后的LAI和产量不切实际的偏低。基于地统计（变异函数）尺度纠正的方法和像元分解降尺度的方法是处理像元异质性和尺度效应的两种手段。然而，这两种方法的一个共同的不足之处是要求一系列的作物生育期的中高空间分辨率数据，并且对异质性像元建模是一个非常复杂的问题，要求严格的近似和一些先验知识，这在实际应用中，不易获取。一种有潜力的方法是基于地面采样LAI、中高分辨率遥感数据和低分辨率遥感数据，构建二级尺度转换方法，生成尺度修正后的LAI。同化尺度修正后的LAI，将会较大提高区域冬小麦产量的估测。以克服直接同化MODIS LAI产品带来的误差是应用数据同化模型估产实践中亟待解决的关键问题。

发明内容

为解决现有技术中存在的如下问题：“如何构建数据同化模型中的尺度转换方法，同化时间序列尺度修正后的LAI到作物生长模型，以减少遥感观测与作物模型之间尺度不匹配造成的数据同化模型误差”，本发明提供一种基于尺度转换和数据同化的农作物产量估测方法。

本发明提供一种基于尺度转换和数据同化的农作物产量估测方法，具体步骤如下：

S1收集研究区内的气象参数、作物参数、土壤参数和管理参数并以之进行WOFOST作物模型的单点尺度标定；之后基于气象参数进行区域化，完成WOFOST作物模型的空间化；

S2基于多时相的遥感影像，结合地面调查获得的作物类型样本点，结合待测农作物的物候特征，构建分类规则，获取作物类型分布图和待测农作物纯度百分比图；

S3对待测农作物生育期时间序列MODIS LAI曲线进行滤波，以消除数据缺失和云污染的影响；基于地面实测样本点的实测LAI与多时相植被指数（TM VI）构建回归统计模型，获得区域30米的TM LAI；

S4融合地面实测样本点LAI、多时相TM LAI和滤波后的MODISLAI信息，构建二级尺度转换模型，生成时间序列尺度调整LAI；

S5基于地面样本点关键物候期的LAI，计算遥感观测和作物模型模拟的标准差，获得待测农作物关键物候期的遥感观测误差和作物模型误差；

S6以S4获得的尺度调整LAI和作物模型模拟LAI，并引入遥感观测误差和作物模型误差，构建四维变分代价函数，使用最优化算法最小化代价函数，通过多次迭代，获得优化作物模型参数；

S7将S6获得的优化作物模型参数代入作物模型，选择待测农作物纯度大于60%的像元，逐像元单元运行作物模型模拟产量，然后汇总到县域行政单元，输出县域待测农作物单产，指导粮食生产。

其中，S1所述气象参数为最高气温、最低气温、总辐射量、水汽压、风速、降水等参数。

其中，S1所述作物参数为农作物分布参数、农作物的物候特征等。

其中，S1所述基于气象参数进行区域化，完成作物模型的空间化是指结合遥感影像对收集的参数进行空间位置的匹配，对气象参数和作物模型所需的积温参数，采用Kriging插值算法，完成参数区域化标定。

其中，S2所述多时相的遥感影像为多时相的Landsat TM。

其中，S2所述获取作物类型分布图和待测农作物纯度百分比图为采用C5.0的决策树分类算法，结合农作物的物候特征和光谱特征，构建分类规则，获取作物类型分布图，采用1km格网套合，生成1km的待测农作物纯度百分比空间分布图，为了减少非待测农作物区域对数据同化结果的影响，选择待测农作物纯度大于60%的像元，作为后续县域尺度产量估测的像元。

其中，S3所述进行滤波，采用上包络线（Savitzky-Golay，S-G）滤波算法，其表达式见公式（1）：

$y_J^{*}=\frac{\underset{i=-m}{\overset{i=m}{\mathrm{\Σ}}}{C}_i{Y}_{j+i}}{N}---\left(1\right)$

式中：Y_j+i表示原始LAI曲线上的一块窗口内的值，m为窗口的半径，N为卷积数目，窗口的宽度为2m+1；Y_j ^*表示滤波后的LAI值；C表示第i个LAI值的滤波系数。

其中，S3所述进行滤波，在对S-G滤波改进的基础上，形成上包络线滤波方法，通过LAI时序变化趋势将LAI值分为真值和假值，用局部循环迭代的方式将S-G滤波值代替假值，与真值合成新的LAI平滑曲线，逐步拟合形成LAI时序数据的上包络线。

其中，S4所述构建二级尺度转换模型由2个关键步骤构成：

1）从S3生成区域的30米待测农作物TM LAI中选择3个关键生育期的TM LAI；

2）利用已选的关键生育期TM LAI与对应时期的S-G MODIS LAI之间的比值系数，调整其他关键生育期S-G MODIS LAI，在待测农作物生育期的上升阶段和下降阶段，分别利用Logistic曲线进行拟合，

Logistic曲线方程见公式（2）：

$y\left(t\right)=\frac{c}{1+e^{a+\mathrm{bt}}}+d---\left(2\right)$

其中，t是MODIS LAI时间序列的索引，y(t)是t时间对应的LAI值，a和b是拟合参数，c+d是最大LAI值，d为LAI初始值，即LAI时间序列中的第一个值。

其中，S5所述待测农作物关键物候期的遥感观测误差是基于实测LAI与尺度修正LAI之间的标准方差计算的；S5所述作物模型误差来源于优化的2个模型参数，TDWI和SPAN。

其中，S6所述使用最优化算法最小化代价函数采用SCE_UA算法对四维变分代价函数进行最小化，四维变分代价函数按公式（3）计算：

$J\left(x\right)={\mathrm{\Σ}}_{k=1}^2{(x_k-x_{k0})}^T{B}^{-1}(x_k-x_{k0})+{\mathrm{\Σ}}_{i=1}^N{(y_i-H_i\left(x\right))}^T{Q}_0^{-1}(y_i-H_i\left(x\right))---\left(3\right)$

其中，k表示代价函数中优化模型参数的个数，x_k代表优化的WOFOST参数在区间范围的数值，x_k0代表优化的WOFOST参数个数的经验值，B表示模型误差，以2个优化模型参数的误差表示，T表示矩阵的转置，N代表时间序列遥感观测数据的次数，y_i表示遥感观测LAI，H_i(x)表示WOFOST模型模拟LAI，Q₀代表遥感LAI的误差。

其中，S7具体步骤为：选择待测农作物纯度大于60%的像元通过不断迭代，重新初始化2个作物模型参数，使得作物模型输出的LAI和产量不断发生变化，然后在代价函数中对比MODIS LAI和作物模型输出LAI的最小差异下，当以下三个收敛条件满足其一即可结束同化。获得优化参数的数值：

①连续5次循环后待优化参数值已收缩到指定的值域范围；

②目标函数值在5次循环后无法提高0.0001%；

③计算代价函数的次数超过10000次；

将优化后的模型参数代入作物模型运行，获得产量，结合像元的待测农作物分布百分比，按县域行政边界汇总，输出区域上的单产结果。

其中，所述待测农作物优选为冬小麦。

本发明还提供所述基于尺度转换和数据同化的农作物产量估测方法在指导农作物生产中的应用。

本发明与现有技术相比，有益效果为：

本发明克服了遥感观测像元与作物模型模拟单元之间尺度不匹配的问题，提高了数据同化模型的精度，适合于县域尺度的冬小麦产量估测。

具体实施方式

下面结合实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

实施例1

以河北保定和衡水地区的冬小麦产量估测为例进一步阐述本发明的技术方案。本实施例的基于尺度转换和数据同化的区域冬小麦产量估测方法的流程包括：

S1收集研究区内的气象参数、作物参数、土壤参数和管理参数并以之进行WOFOST作物模型的单点尺度标定；之后基于气象参数进行区域化，完成作物模型的空间化。具体如下：

选择河北保定和衡水地区为研究区域，该区为河北省重要的冬小麦产区。获取以下数据：选取6个气象站点的最高/最低气温、总辐射、水汽压、风速、降水模型所需的6个气象要素；研究区内农业气象试验站采集的土壤参数和作物参数等；2009年河北省分县冬小麦产量统计数据；MODIS15A2数据产品，统一所有数据的坐标投影为AlbersConical Equal Area，空间参考为WGS1984。

对于WOFOST模型中不敏感的作物参数，或敏感程度较高但是变化范围较小的参数，通过查阅文献或者直接采用WOFOST模型的默认值来确定（如发育最适光长DLO，出苗最低温度下限TBASE，叶的维持呼吸最用RML等）；对于敏感程度高、变化范围大的作物参数，能够通过查阅本研究区的相关文献或者利用实测数据计算获得的（如积温TSUM1，TSUM2，叶同化物转化率CVL，比叶面积SLATB等）。也可以研究区的先验知识等确定参数可能的取值范围和数值（如叶片在35℃时的生命期SPAN，最大CO2同化速率AMAXTB，单叶片同化CO2的光能利用率EFFTB等）。

利用研究区6个气象站点的观测数据，结合克里金空间插值方法对2008～2009年逐日最高气温、最低气温、总辐射量、水汽压、风速、降水量六个要素进行插值，并将其转换为WOFOST模型的标准格式，得到研究区每个格点相对应的气象数据。对于空间范围内容易获取的参数，例如播种期到出苗期的积温、出苗期到开花期的积温、开花期到成熟期的积温等，利用气象观测台站获得的日平均气温与冬小麦的生物学零度相减得到每日的积温，通过累加不同物候阶段的每日积温得到以上三个积温变量，同样采用与气象数据插值相同的克里金插值法，获得整个研究区内空间变化的积温参数。

S2利用研究区多时相的美国专题制图仪（Thematic Mapper，TM）影像，结合地面调查获得的作物类型样本点，采用C5.0的决策树分类算法，结合农作物的物候特征，构建分类规则，获取作物类型分布图，与1km格网套合，生成1km的冬小麦纯度百分比空间分布图，为了减少非冬小麦区域对数据同化结果的影响，选择冬小麦纯度大于60%的像元，作为后续县域尺度产量估测的基础。

S3对实验区整个生育期内的MODIS LAI数据按时间序列合成，对每个栅格单元生成时间序列曲线。具体通过NASA提供的MRT投影转换工具，将河北地区的数据进行镶嵌、投影转换和格式转换，选取1到177天（儒略日），共23个时相影像进行叠加，生成基于栅格单元的时间序列曲线。时间序列MODIS LAI曲线进行滤波，以消除云污染造成的数据缺失。具体使用上包络线（Savitzky-Golay，S-G）滤波解决叶面积指数变化过程不连续的问题，其原理表达如下：

$y_J^{*}=\frac{\underset{i=-m}{\overset{i=m}{\mathrm{\Σ}}}{C}_i{Y}_{j+i}}{N}---\left(1\right)$

式中：Y_j+i表示原始LAI曲线上的一块窗口内的值，m为窗口的半径，N为卷积数目，窗口的宽度为2m+1；Y_j ^*表示滤波后的LAI值；C表示第i个LAI值的滤波系数。

在对S-G滤波改进的基础上，形成上包络线滤波方法，通过LAI时序变化趋势将LAI值分为真值和假值，用局部循环迭代的方式将S-G滤波值代替假值，与真值合成新的LAI平滑曲线，逐步拟合形成LAI时序数据的上包络线。具体步骤如下：

（1）对有云污染的LAI值点进行线性插值。由于云覆盖或是大气效应的影响，会使一些像元值带有较大的噪声，对这些像元采用其周边可信像元值进行空间线性插值，生成初始的LAI数据。

（2）用S-G滤波进行LAI时序变化趋势的提取。根据假设，LAI时序数据应遵循植被物候动态变化的渐进特征，在LAI时序中突变被认为是由于云或大气产生的噪音造成的。因此，对研究区内每个冬小麦像元应用S-G滤波进行平滑，得到LAI时间序列的变化趋势数据LAItr。

（3）确定LAI时间序列中每个值点的权重Wi，为后续的重构结果提供判断依据。Wi的计算公式为：

$W_i=\left\{\begin{array}{cc}1& {\mathrm{LAI}}_i\&GreaterEqual;{\mathrm{LAI}}_i^{\mathrm{tr}}\\ 1-d_i/d_{\max }& {\mathrm{LAI}}_i<{\mathrm{LAI}}_i^{\mathrm{tr}}\end{array}\right.$

其中， $d_i=|{\mathrm{LAI}}_i-{\mathrm{LAI}}_i^{\mathrm{tr}|},$ d_max为di的最大值。

（4）生成新的LAI时间序列。由以上理论可知，含有噪声的LAI值要比时间变化趋势曲线上的LAItr值小，对这些点进行简单的替换生成新的LAI1时序曲线。

${\mathrm{LAI}}_i^1=\left\{\begin{array}{cc}{\mathrm{LAI}}_i& {\mathrm{LAI}}_i\&GreaterEqual;{\mathrm{LAI}}_i^{\mathrm{tr}}\\ {\mathrm{LAI}}_i^{\mathrm{tr}}& {\mathrm{LAI}}_i<{\mathrm{LAI}}_i^{\mathrm{tr}}\end{array}\right.$

（5）对LAI1进行S-G滤波，再将滤波后的数据与空间插值后的LAI值进行比较。同样，如果空间插值后的LAI值小于新生成的LAI j，用LAI j代替，生成新序列LAI j+1，迭代进行滤波重构，滤波效果由F_k的大小来评定：

$F_k=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\right|{\mathrm{LAI}}_i^{k+1}-{\mathrm{LAI}}_i|\&times;W_i)$

LAIik为第i个像元在第k次滤波重构后的LAI值，Wi为步骤（3）计算得到的权重，Fk为第k次滤波后的迭代终止判断参数。退出迭代的条件为F_k-1≥F_k≤F_k+1。

获取了实验区3景无云的TM数据，时相分别为：2009-3-14、2009-5-17、2009-6-10，与对应的实测LAI分别构建回归统计模型。

2009年3月14日：由于3月份中旬冬小麦处于返青阶段，裸土对植被指数的影响较大，因此选择SAVI植被指数来建立LAI与TM VI之间的统计模型为：

$\mathrm{LAI}=\frac{\ln [(1-\mathrm{SAVI}/1.2581)/0.9130]}{-0.8377}$

决定系数R²=0.849。

2009年5月17日，建立实测LAI与TM NDVI之间的统计模型为：

$\mathrm{LAI}=\frac{\ln [(1-\mathrm{NDVI}/1.0866)/3.3790]}{-0.3994},$

决定系数R²=0.742。

2009年6月10日，建立实测LAI与TM NDVI之间的统计模型为:

$\mathrm{LAI}=\frac{\ln [(1-\mathrm{NDVI}/9.7639)/1.0081]}{-0.0155}$

决定系数R²=0.874。

基于上述的3个回归统计模型，获得研究区30米冬小麦的空间分布图。

S4融合地面实测样本点LAI、多时相TM LAI和滤波后的MODISLAI信息，构建二级尺度转换模型，生成时间序列尺度调整LAI。

二级尺度转换模型包括2个步骤，首先，按S3步骤中生成的3个物候期的实测LAI和TM VI回归统计模型，生成区域30米冬小麦TMLAI；第二步，利用已有的关键生育期TM LAI与对应时期的S-GMODIS LAI之间的比值系数，调整其他关键生育期S-G MODIS LAI，在冬小麦生育期的上升阶段（返青到拔节期），获得4个物候期的调整LAI（孕穗期，拔节期，起身期，返青期）；在冬小麦的下降阶段（拔节到成熟期），获得4个物候期的调整LAI（孕穗期，抽穗期，开花期，成熟期），分别利用Logistic曲线进行拟合。

Logistic曲线方程如下：

$y\left(t\right)=\frac{c}{1+e^{a+\mathrm{bt}}}+d$

其中，t是MODIS LAI时间序列的索引，y(t)是t时间对应的LAI值，a和b是拟合参数，c+d是最大LAI值，d为LAI初始值，即LAI时间序列中的第一个值。

S5基于地面样本点关键物候期的LAI，计算遥感观测和作物模型模拟的标准差，获得冬小麦关键物候期的遥感观测误差和作物模型误差。

模型误差认为是来源于优化的2个参数（TDWI和SPAN）,经过分析计算获得尺度修正LAI在7个物候期的观测误差为0.34；模型误差设置为0.52。

S6以S4获得的尺度调整LAI和作物模型模拟LAI，并引入遥感观测误差和作物模型误差，构建四维变分代价函数，使用最优化算法最小化代价函数，通过多次迭代，获得优化作物模型参数。

在S1进行作物模型标定的基础上，在实验区域内，运行WOFOST作物模型模拟LAI，结合S4中生成的尺度修正LAI，引入S5中的遥感观测与模型模拟误差,构建四维变分代价函数，如下：

$J\left(x\right)={\mathrm{\&Sigma;}}_{k=1}^2{(x_k-x_{k0})}^T{B}^{-1}(x_k-x_{k0})+{\mathrm{\&Sigma;}}_{i=1}^N{(y_i-H_i\left(x\right))}^T{Q}_0^{-1}(y_i-H_i\left(x\right))$

其中，k表示代价函数中优化模型参数的个数，x_k代表优化的WOFOST参数在区间范围的数值，x_k0代表优化的WOFOST参数个数的经验值，B表示模型误差，以2个优化模型参数的误差表示，T表示矩阵的转置，N代表时间序列遥感观测数据的次数，y_i表示遥感观测LAI，H_i(x)表示WOFOST模型模拟LAI，Q₀代表遥感LAI的误差。

选择冬小麦纯度大于60%的像元，运用优化算法，不断迭代重新初始化2个WOFOST模型参数，使得模型输出的LAI和产量发生变化，然后在代价函数中对比MODIS LAI和WOFOST模型输出LAI的最小差异下，获得最优化参数。SCE_UA算法用来在初始参数空间中搜索全局最优解，使得代价函数快速收敛。当以下三个收敛条件满足其一即可结束同化，获得优化参数的数值。

（1）连续5次循环后待优化参数值已收缩到指定的值域范围；

（2）目标函数值在5次循环后无法提高0.0001%；

（3）计算代价函数的次数超过10000次。

S7将S6获得的优化作物模型参数代入作物模型，选择冬小麦纯度大于60%的像元，逐像元单元运行作物模型模拟产量，然后汇总到县域行政单元，输出县域冬小麦单产。将优化后的TDWI和SPAN参数代入WOFOST模型运行，获得产量。结合像元的冬小麦百分比，按县域行政边界汇总，输出区域上的冬小麦单产结果。

本发明通过利用S-G滤波算法对冬小麦生育期时间序列MODISLAI曲线进行滤波，消除了数据缺失和云污染的影响，生成了平滑的LAI时间序列曲线，为二级尺度转换模型提供了高精度的数据源。然后通过二级尺度转换模型生成了时间序列尺度调整LAI，有效融合了MODIS LAI时间序列的趋势信息及TM LAI信息，提高了LAI的精度，为同化过程提供了高精度的数据源。最后，通过同化调整后的LAI与作物模型模拟的LAI，提高了区域上冬小麦产量估测的精度。

粮食的生产者与消费者都需要及时准确地了解粮食产量信息，根据本方法，可以在冬小麦成熟期大面积地获得产量数据，为国家有关部门进行粮情判断、粮食调控等科学决策等提供重要的科学依据，并且可以作为粮食贸易的重要依据。

本发明的方法也可以用于其他农作物产量的估测。

虽然，上文中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。

Claims (10)

1.一种基于尺度转换和数据同化的农作物产量估测方法，其特征在于，具体步骤如下：

S1收集研究区内的气象参数、作物参数、土壤参数和管理参数并以之进行WOFOST作物模型的单点尺度标定；之后基于气象参数进行区域化，完成WOFOST作物模型的空间化；

S2基于多时相的遥感影像，结合地面调查获得的作物类型样本点，结合待测农作物的物候特征，构建分类规则，获取作物类型分布图和待测农作物纯度百分比图；

S3对待测农作物生育期时间序列MODIS LAI曲线进行滤波，以消除数据缺失和云污染的影响；基于地面实测样本点的实测LAI与多时相植被指数TM VI构建回归统计模型，获得区域30米的TM LAI；

S4融合地面实测样本点LAI、多时相TM LAI和滤波后的MODISLAI信息，构建二级尺度转换模型，生成时间序列尺度调整LAI；

S5基于地面样本点关键物候期的LAI，计算遥感观测和作物模型模拟的标准差，获得待测农作物关键物候期的遥感观测误差和作物模型误差；

S6以S4获得的尺度调整LAI和作物模型模拟LAI，并引入遥感观测误差和作物模型误差，构建四维变分代价函数，使用最优化算法最小化代价函数，通过多次迭代，获得优化作物模型参数；

S7将S6获得的优化作物模型参数代入作物模型，选择待测农作物纯度大于60%的像元，逐像元单元运行作物模型模拟产量，然后汇总到县域行政单元，输出县域待测农作物单产，指导粮食生产。

2.权利要求1所述基于尺度转换和数据同化的农作物产量估测方法，其特征在于，S1所述基于气象参数进行区域化，完成作物模型的空间化是指结合遥感影像对收集的参数进行空间位置的匹配，对气象参数和作物模型所需的积温参数，采用Kriging插值算法，完成参数区域化标定。

3.权利要求1所述基于尺度转换和数据同化的农作物产量估测方法，其特征在于，S2所述多时相的遥感影像为多时相的LandsatTM。

4.权利要求1所述基于尺度转换和数据同化的农作物产量估测方法，其特征在于，S2所述获取作物类型分布图和待测农作物纯度百分比图为采用C5.0的决策树分类算法，结合农作物的物候特征和光谱特征，构建分类规则，获取作物类型分布图，采用1km格网套合，生成1km的待测农作物纯度百分比空间分布图，为了减少非待测农作物区域对数据同化结果的影响，选择待测农作物纯度大于60%的像元，作为后续县域尺度产量估测的像元。

5.权利要求1所述基于尺度转换和数据同化的农作物产量估测方法，其特征在于，S3所述进行滤波，采用上包络线（Savitzky-Golay，S-G）滤波算法，其表达式见公式（1）：

$y_J^{*}=\frac{\underset{i=-m}{\overset{i=m}{\mathrm{\&Sigma;}}}{C}_i{Y}_{j+i}}{N}---\left(1\right)$

式中：Y_j+i表示原始LAI曲线上的一块窗口内的值，m为窗口的半径，N为卷积数目，窗口的宽度为2m+1；Y_j ^*表示滤波后的LAI值；C表示第i个LAI值的滤波系数。

6.权利要求1所述基于尺度转换和数据同化的农作物产量估测方0法，其特征在于，S4所述构建二级尺度转换模型由2个关键步骤构成：

1）从S3生成区域的30米待测农作物TM LAI中选择3个关键生育期的TM LAI；

2）利用已选的关键生育期TM LAI与对应时期的S-G MODIS LAI之间的比值系数，调整其他关键生育期S-G MODIS LAI，在待测农作物生育期的上升阶段和下降阶段，分别利用Logistic曲线进行拟合，Logistic曲线方程见公式（2）：

$y\left(t\right)=\frac{c}{1+e^{a+\mathrm{bt}}}+d---\left(2\right)$

其中，t是MODIS LAI时间序列的索引，y(t)是t时间对应的LAI值，a和b是拟合参数，c+d是最大LAI值，d为LAI初始值，即LAI时间序列中的第一个值。

7.权利要求1所述基于尺度转换和数据同化的农作物产量估测方法，其特征在于，S5所述待测农作物关键物候期的遥感观测误差是基于实测LAI与尺度修正LAI之间的标准方差计算的；S5所述作物模型误差来源于优化的2个模型参数，TDWI和SPAN。

8.权利要求1所述基于尺度转换和数据同化的农作物产量估测方法，其特征在于，S6所述使用最优化算法最小化代价函数采用SCE_UA算法对四维变分代价函数进行最小化，四维变分代价函数按公式（3）计算：

$J\left(x\right)={\mathrm{\&Sigma;}}_{k=1}^2{(x_k-x_{k0})}^T{B}^{-1}(x_k-x_{k0})+{\mathrm{\&Sigma;}}_{i=1}^N{(y_i-H_i\left(x\right))}^T{Q}_0^{-1}(y_i-H_i\left(x\right))---\left(3\right)$

其中，k表示代价函数中优化模型参数的个数，x_k代表优化的WOFOST参数在区间范围的数值，x_k0代表优化的WOFOST参数个数的经验值，B表示模型误差，以2个优化模型参数的误差表示，T表示矩阵的转置，N代表时间序列遥感观测数据的次数，y_i表示遥感观测LAI，H_i(x)表示WOFOST模型模拟LAI，Q₀代表遥感LAI的误差。

9.权利要求1所述基于尺度转换和数据同化的农作物产量估测方法，其特征在于，S7具体步骤为：选择待测农作物纯度大于60%的像元通过不断迭代，重新初始化2个作物模型参数，使得作物模型输出的LAI和产量不断发生变化，然后在代价函数中对比MODIS LAI和作物模型输出LAI的最小差异下，当以下三个收敛条件满足其一即可结束同化。获得优化参数的数值：

①连续5次循环后待优化参数值已收缩到指定的值域范围；

②目标函数值在5次循环后无法提高0.0001%；

③计算代价函数的次数超过10000次；

将优化后的模型参数代入作物模型运行，获得产量，结合像元的待测农作物分布百分比，按县域行政边界汇总，输出区域上的单产结果。

10.权利要求1-9任一项所述基于尺度转换和数据同化的农作物产量估测方法在指导农作物生产中的应用。