CN102651096A - 同化叶面积指数时序曲线特征的冬小麦估产方法 - Google Patents

Abstract

本发明公开了一种同化叶面积指数时序曲线特征的冬小麦估产方法，通过如下步骤实施：S1：对作物模型进行全局敏感性分析，并完成作物模型参数区域化；S2：合成MODIS LAI时间序列曲线；S3：对LAI时序曲线进行滤波处理；S4：进行曲线拟合并提取出关键特征点；S5：在区域范围运行作物模型，对模拟得到的LAI时间序列进行曲线拟合并提取出曲线上的关键特征点；S6：按S4和S5各自获得3类关键特征点的日期建立代价函数，并以遥感观测误差为权重求和得到总代价函数，对总代价函数进行最小化，使得代价函数快速收敛，最终在收敛条件被满足时，按行政边界汇总，输出产量结果。本发明提高了同化精度，克服了MODIS LAI产品系统偏低的影响，适合于区域尺度冬小麦产量估测。

Description

同化叶面积指数时序曲线特征的冬小麦估产方法

技术领域

本发明涉及农业遥感技术领域，特别涉及一种同化叶面积指数时序曲线特征的冬小麦估产方法。

背景技术

在作物估产应用过程中，基于作物光合、呼吸、蒸腾、营养等机理过程的作物模型依靠其内在的物理过程和动力学机制，可以准确模拟作物对象在时间和空间上的连续演进，能够准确地模拟单点作物的生长发育状况及产量。而在应用到区域尺度时，由于地表、近地表环境非均匀性，导致模型中一些宏观资料的获取和参数的区域化非常困难。卫星遥感方法具有空间连续和时间动态变化的优势，能够有效解决区域参数获取困难这一瓶颈。然而遥感对地观测由于受卫星时空分辨率等因素的制约，还不能真正揭示作物生长发育和产量形成的内在过程机理、个体生长发育状况及其与环境气象条件的关系，而这正是作物模型的优势所在。因此通过同化方法将二者结合，可以发挥各自的优势，从而快速、动态地监测大面积的作物长势与产量。

优化算法是同化方法的一类重要方法，其基本思想主要是采用优化算法通过多次迭代最小化遥感反演的状态变量与作物模型模拟的状态变量之间的差异，重新初始化作物模型参数达到对作物模型优化的目的。目前比较通用的方法主要有以下两种：

1、通过对比遥感反演和模型模拟获得的作物的叶面积指数LAI；

2、对比遥感反演和模型模拟获得的作物光谱信息或植被指数。

MODIS数据由于具有较高的时间分辨率，在农作物长势监测等方面具有优势。然而在我国华北地区的冬小麦主产区，由于MODIS数据空间分辨率较低，混合像元效应严重，导致MODIS LAI产品存在系统偏低的误差，直接同化MODIS LAI产品将导致更差的同化效果。而大量的野外观测实验发现，冬小麦生育期时间序列MODIS LAI产品的变化趋势与实测变化趋势十分吻合。因此，同化MODIS LAI时间序列曲线的关键特征到作物生长模型，以克服直接同化MODISLAI产品带来的误差是亟待解决的问题。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是：如何同化MODIS LAI时间序列曲线的关键特征到作物生长模型，以避免MODIS-LAI数据由于混合像元效应造成的LAI数值偏低的问题。

(二)技术方案

为解决上述技术问题，本发明提供了一种同化叶面积指数时序曲线特征的冬小麦估产方法，包括以下步骤：

S1：进行实验区作物模型的敏感性分析，得到需要标定和同化的模型参数，对遥感数据与地面参数进行空间匹配，使用插值法实现部分作物模型参数区域化标定，完成数据准备；

S2：对实验区整个生育期内的MODIS LAI数据按时间序列合成，对每个栅格单元生成时间序列曲线；

S3：对S2中获得的时间序列MODIS LAI曲线进行滤波，以消除数据缺失和云污染的影响；

S4：将S3中滤波后的时间序列曲线进行曲线拟合，得到每个栅格单元的MODIS LAI时间序列的曲线方程，并提取出该曲线上的关键特征点，所述关键特征点是曲率最大值点、斜率最大值点和叶面积最大值点；

S5：在S1进行作物模型标定的基础上，在实验区域内运行WOFOST作物模型，对模拟的LAI时间序列进行曲线拟合，并提取出曲线上的关键特征点，所述关键特征点是曲率最大值点、斜率最大值点和叶面积最大值点；

S6：S4和S5各自获得3类关键特征点，按类别对点的日期建立代价函数，并以遥感观测误差为权重求和得到总代价函数，对总代价函数进行最小化，使其不断重新初始化S1获得的需进行同化的作物敏感性的模型参数，在初始参数空间中搜索全局最优解，使得代价函数快速收敛，最终在收敛条件被满足时，同化过程结束，按行政边界汇总，输出产量结果。

其中，所述步骤S1中所采用的作物模型为WOFOST模型，采集研究区内的土壤、气象和作物参数，对影像和采集的参数进行空间位置的匹配，对于不敏感模型参数直接使用默认值，对于只对产量敏感，对LAI时序曲线特征不敏感的模型参数使用反距离权重插值法对遥感影像的每个像素赋予参数值，完成参数区域化标定；采用扩展的傅里叶振幅敏感度检验法EFAST对WOFOST模型的作物参数和土壤与管理参数进行全局敏感性分析，按如下公式(1)求解一阶敏感性指数和总敏感性指数，

$S_i=\frac{V_i}{V}$ $S_{\mathrm{ij}}=\frac{V_{\mathrm{ij}}}{V}$ $S_{\mathrm{ijm}}=\frac{V_{\mathrm{ijm}}}{V}...S_{12...i...k}=\frac{V_{12...i...k}}{V}$

S_T.i＝S_i+S_ij+S_ijm+…+S_12…i…k    (1)

$V\left(X\right)={\mathrm{\Σ}}_i{V}_i+{\mathrm{\Σ}}_{i\≠j}{V}_{\mathrm{ij}}+{\mathrm{\Σ}}_{i\≠j\≠m}{V}_{\mathrm{ijm}}+V_{12...i...k}$

式中，i，j，m为作物模型参数的序列号，k为模型参数的总个数，V(X)为模型结果X的总方差；V_i为参数x_i的方差；V_ij～V_12...k为各参数相互作用的方差，S_i、S_ij、S_ijm为参数x_i的一阶、二阶、三阶敏感性指数，S_T.i为参数x_i的总敏感性指数。

其中，所述步骤S3中采用上包络线滤波，滤波过程按如下公式(2)：

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

式(2)中Y_j+i表示原始LAI曲线上的一块窗口内的值，m为窗口的半径，N为卷积数目，窗口的宽度为2m+1，

表示滤波后窗口中心的LAI值，C_i表示第i个LAI值的滤波系数。

其中，所述步骤S4或S5中所采用的曲线拟合为Logistic曲线拟合，曲线方程按如下公式(3)：

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

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

其中，所述步骤S4或S5中所要提取的关键特征点按如下公式(4)计算曲率：

$K=\frac{\mathrm{d\α}}{\mathrm{ds}}=\frac{b^2\mathrm{cz}(1-z){(1+z)}^3}{{[{(1+z)}^4+{\left(\mathrm{bcz}\right)}^2]}^{\frac{3}{2}}}---\left(4\right)$

其中，z＝e^a+bt，按如下公式(5)计算斜率：

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

其中，所述步骤S6中遥感观测误差表示的权重按如下公式(6)求得：

$w_i=1-\frac{v_i}{\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{v}_i}---\left(6\right)$

其中，

$v_i=\frac{|\left(\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{LAI}}_{\mathrm{kj}}^i\right)-{\mathrm{LAI}}_{\mathrm{RS}}^i|}{\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{LAI}}_{\mathrm{kj}}^i}$

其中m是在研究区内土地单元个数，n是在每个样地内布设的采样块个数，LAIⁱ _kj是在第i个曲线特征点处第k个单元中第j个采样块内的单点测量LAI值，LAIⁱ _RS是土地单元的MODIS LAI值，w_i是三个曲线特征点处的遥感观测误差权重。

其中，所述步骤S6中采用SCE_UA算法对代价函数进行最小化，并结合所述权重，代价函数按如下公式(7)：

$J=w_1{(T_1^R-T_1^M)}^2+w_2{(T_2^R-T_2^M)}^2+w_3{(T_3^R-T_3^M)}^2---\left(7\right)$

式(7)中，T₁，T₂，T₃分别代表曲率最大值点、斜率最大值点、叶面积最大值点处的儒略日时间，加上标R代表遥感观测值、加上标M代表模型模拟值，w_i代表该时间点i处以遥感观测误差表示的权重，J为代价函数值。

其中，所述步骤S6中的收敛条件为以下条件之一：

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

(2)目标函数值在5次循环后无法提高0.0001％；

(3)计算代价函数的次数超过10000次，

WOFOST模型将按照同化获得的模型参数运行，输出产量结果。

(三)有益效果

本发明所提供的同化叶面积指数时序曲线特征的冬小麦估产方法，通过建立生育期内冬小麦的LAI时间序列曲线，并对其进行S_G滤波处理，采用与以往同化方法不同的策略，提取其关键曲线特征点的日期来构建代价函数，从而有效地避免了MODIS-LAI数据由于混合像元效应造成的LAI数值偏低的问题，为大面积区域尺度冬小麦或其它农作物产量估测提供了办法。

附图说明

图1是本发明实施例的一种同化叶面积指数时序曲线特征的冬小麦估产方法流程图；

图2是时序曲线特征点提取效果示意图；

图3为河北省南部地区区域验证效果示意图。

具体实施方式

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

本实施例中，以估测河北南部研究区的冬小麦的产量为例进一步阐述本发明的技术方案。本实施例的同化叶面积指数时序曲线特征的冬小麦估产方法的流程如图1所示，包括：

步骤S1，进行实验区作物模型的敏感性分析，对遥感数据与地面参数进行空间匹配，并完成作物模型参数区域化。

在分析敏感性之前需要采集数据，选择河北南部地区为研究区域，该区为洪冲积平原，地处太行山脉以东，北京市以南。获取以下数据：2000年1∶10万的研究区土地利用数据；根据研究区外包络范围，选取12个气象站点的最高/最低气温、总辐射、水汽压、风速、降水模型所需的6个气象要素；研究区内农业气象试验站采集的土壤参数和作物参数；经纬度、高程等控制参数；农业气象数据和08年河北省分县冬小麦产量数据；MODIS 15A2数据产品，将含有地理位置信息的数据统一坐标，完成空间匹配。

利用研究区同时相的美国专题制图仪(Thematic Mapper，TM)影像解译得到冬小麦种植区域，生成1公里的数据掩膜，然后利用MODIS-LAI数据中包含的MOD12土地利用信息并设定阈值剔除不符合冬小麦物候规律的像素生成新的掩膜，将两种掩膜进行与运算，生成最终的运行区域。对模型的产量输出和气象数据统一度量，将气象数据和作物参数反距离插值生成1公里每像素的栅格数据。

根据荷兰瓦赫宁根大学开发的世界粮食研究模型(WOFOST模型)在河北地区冬小麦的初步标定结果，借助于敏感性与不确定性专业分析软件Simlab的EFAST模块实现WOFOST模型的作物参数全局敏感性分析：首先，根据设定WOFOST模型各输入参数在取值范围内均匀分布，利用蒙特卡洛方法分别随机采样3000次和1000次；然后，将生成的参数组合输入到WOFOST模型中获得冬小麦模拟产量；最后，利用扩展的傅里叶振幅敏感度检验法EFAST计算各个参数的一阶敏感型指数和总敏感性指数，得到对于冬小麦模拟产量敏感度较高的WOFOST作物参数(敏感性分析就是为了找到这些参数进行本地化调整，这些参数调整确定后就可以进行同化，同化过程中模型就会生成LAI曲线)，对于不敏感模型参数直接使用默认值，对于只对产量敏感，对LAI时序曲线特征不敏感的模型参数，利用地面样点实测数据采用反距离插值法获得整个实验区每个栅格单元的值进行标定，完成模型的本地化调整。对于产量和LAI时序曲线特征都敏感的参数进行后续步骤的同化过程。

按如下公式求解一阶敏感性指数和总敏感性指数，

$S_i=\frac{V_i}{V}$ $S_{\mathrm{ij}}=\frac{V_{\mathrm{ij}}}{V}$ $S_{\mathrm{ijm}}=\frac{V_{\mathrm{ijm}}}{V}...S_{12...i...k}=\frac{V_{12...i...k}}{V}$

S_T.i＝S_i+S_ij+S_ijm+…+S_12…i…k

V(X)＝∑_iV_i+∑_i≠jV_ij+∑_i≠j≠mV_ijm+V_12...k

式中，i，j，m为作物模型参数的序列号，k为模型参数的总个数。V(X)为模型结果X的总方差；V_i为参数x_i的方差；V_ij～V_12...k为各参数相互作用的方差。S_i、S_ij、S_ijm为参数x_i的一阶、二阶、三阶敏感性指数，S_T.i为参数x_i的总敏感性指数。

步骤S2，对实验区整个生育期内的MODIS LAI数据按时间序列合成，对每个栅格单元生成时间序列曲线。具体通过NASA提供的MRT投影转换工具，将河北地区的数据进行镶嵌、投影转换和格式转换，中选取1至177天(儒略日)，共23景影像进行叠加，生成每个栅格单元的时间序列曲线。

步骤S3，对步骤S2中获得的时间序列MODIS LAI曲线进行滤波，以消除云污染造成的数据缺失。具体利用上包络线(Savitzky-Golay，SG)滤波解决叶面积指数变化过程不连续的问题(见附图2)，其原理表达如下：

$Y_j^{*}=\frac{\underset{i=-m}{\overset{i=m}{\mathrm{\Σ}}}{C}_i{Y}_{j+i}}{N}$

式中：Y_j+i表示原始LAI曲线上的一块窗口内的值，m为窗口的半径，N为卷积数目，窗口的宽度为2m+1；

表示滤波后的LAI值；C表示第i个LAI值的滤波系数。这个多项式的设计是为了保留高的数值而减少异常值。具体操作如下：

(1)对初始的LAI时间序列进行SG滤波，得到的平滑后的结果，并分别保存平滑前和平滑后的序列。

(2)对比上一步的两个序列，用如下公式生成新的序列，作为初始序列，

$N_i^t=\left\{\begin{array}{cc}{O}_i& \mathrm{if}{O}_i^{t-1}\&GreaterEqual;N_i^{t-1}\\ N_i^{t-1}& \mathrm{if}{O}_i^{t-1}<N_i^{t-1}\end{array}\right.$

其中O和N分别为初始和滤波后的LAI值，上标t代表迭代次数，下标i是LAI时间序列索引。

(3)重复(1)(2)两步，直到整个序列∑N_i ^t小于指定的阈值0.1。

步骤S4，将S3中滤波后的时间序列曲线进行曲线拟合，得到每个栅格单元的MODIS LAI时间序列的曲线方程，并提取出该曲线上的关键特征点。关键特征点是曲率最大值点、斜率最大值点和叶面积最大值点。

使用Logistic曲线拟合法对滤波后的时间序列曲线进行拟合(见附图2)，其公式表达如下：

$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时间序列中的第一个值。拟合完成后，再对拟合了的Logistic方程求其曲率，得到如下公式：

$K=\frac{\mathrm{d\&alpha;}}{\mathrm{ds}}=\frac{b^2\mathrm{cz}(1-z){(1+z)}^3}{{[{(1+z)}^4+{\left(\mathrm{bcz}\right)}^2]}^{\frac{3}{2}}}$

其中，z＝e^a+bt。同时还要对方程求导以获得曲线斜率：

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

其中，a，b，c与公式(3)中含义相同。输出曲率最大值点的儒略日为T₁，斜率最大值点的儒略日为T₂，另外拟合曲线中的最大值点的儒略日为T₃(见附图2)。

步骤S5，在S1进行作物模型标定的基础上，在实验区域内运行WOFOST作物模型，对模拟的LAI时间序列进行曲线拟合，并提取出曲线上的关键特征点。关键特征点是曲率最大值点T₁、斜率最大值点T₂和叶面积最大值点T₃。

在区域范围运行作物模型。针对生成的最终运行区域中的所有像素运行作物模型，和S4类似，对结果中的LAI时间序列利用Logistic曲线拟合来提取T₁、T₂、T₃。

步骤S6，S4和S5各自获得3类关键特征点，按类别对点的日期建立代价函数，并以遥感观测误差为权重求和得到总代价函数，对总代价函数进行最小化，使其不断重新初始化S1获得的需进行同化的敏感性的作物模型参数，在初始参数空间中搜索全局最优解，使得代价函数快速收敛，最终在收敛条件被满足时，同化过程结束，按行政边界汇总，输出产量结果。

在研究区内布设的3至4个有代表性的1km×1km土地单元，在每个样地内布设9～12个30m×30m采样块，每个采样块内进行一定数量的单点测量，结合中、高分辨率遥感影像，建立基本采样块与样地之间的尺度转换关系，按公式(6)获得遥感观测误差在三个点处的误差权重：

$w_i=1-\frac{v_i}{\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}{v}_i}---\left(6\right)$

其中：

$v_i=\frac{|\left(\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{LAI}}_{\mathrm{kj}}^i\right)-{\mathrm{LAI}}_{\mathrm{RS}}^i|}{\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{LAI}}_{\mathrm{kj}}^i}$

其中，m是在研究区内土地单元个数，n是在每个样地内布设的采样块个数，LAIⁱ _kj是在第i个曲线特征点处第k个单元中第j个采样块内的单点测量LAI值，LAIⁱ _RS是该时刻土地单元的MODIS LAI值，w_i是三个曲线特征点处的遥感观测权重。建立代价函数开始优化过程，代价函数表达如下：

$J=w_1{(T_1^R-T_1^M)}^2+w_2{(T_2^R-T_2^M)}^2+w_3{(T_3^R-T_3^M)}^2$

式中，T₁、T₂、T₃分别代表曲率最大值点、斜率最大值点、叶面积最大值点处的儒略日时间；加上标R代表遥感观测值、加上标M代表模型模拟值；w_i代表该时间点处以遥感观测误差表示的误差权重；J为代价函数值。

对每个像素通过不断重新初始化步骤S1中确定的WOFOST模型敏感参数，使得模型输出的物候信息和产量不断发生变化。然后在代价函数中对比MODIS提取和模型输出的两组物候日期，SCE_UA算法用来在初始参数空间中搜索全局最优解，使得代价函数快速收敛。当以下三个收敛条件满足其一即可结束同化，WOFOST模型将按照同化后的模型参数运行，获得产量结果。按行政边界汇总，输出区域上的产量结果(见附图3，左边为同化前的产量结果图示，右边为同化后的产量结果图示)。

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

(2)目标函数值在5次循环后无法提高0.0001％；

(3)计算代价函数的次数超过10000次。

本发明实施例所述的同化叶面积指数时序曲线特征的冬小麦估产方法，融合了遥感数据和作物模型的优势，把遥感观测的时间序列LAI变化特征提取出来作为同化变量，通过优化算法实现了遥感和模型的同化，解决了MODIS LAI偏低不能直接应用于数据同化中的问题。通过LAI时间序列变化特征的同化得到的冬小麦产量，与未同化相比较，均方根误差RMSE在2007减小了5.3％，在2008年减小了15.8％，同化后使作物模型产量估算的精度有显著提高，区域空间分布趋势也符合规律。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (8)

1.一种同化叶面积指数时序曲线特征的冬小麦估产方法，其特征在于，包括以下步骤：

S1：进行实验区作物模型的敏感性分析，得到需要标定和同化的模型参数，对遥感数据与地面参数进行空间匹配，使用插值法实现部分作物模型参数区域化标定，完成数据准备；

S2：对实验区整个生育期内的MODIS LAI数据按时间序列合成，对每个栅格单元生成时间序列曲线；

S3：对S2中获得的时间序列MODIS LAI曲线进行滤波，以消除数据缺失和云污染的影响；

S4：将S3中滤波后的时间序列曲线进行曲线拟合，得到每个栅格单元的MODIS LAI时间序列的曲线方程，并提取出该曲线上的关键特征点，所述关键特征点是曲率最大值点、斜率最大值点和叶面积最大值点；

S5：在S1进行作物模型标定的基础上，在实验区域内运行WOFOST作物模型，对模拟的LAI时间序列进行曲线拟合，并提取出曲线上的关键特征点，所述关键特征点是曲率最大值点、斜率最大值点和叶面积最大值点；

S6：S4和S5各自获得3类关键特征点，按类别对点的日期建立代价函数，并以遥感观测误差为权重求和得到总代价函数，对总代价函数进行最小化，使其不断重新初始化S1获得的需进行同化的敏感性的作物模型参数，在初始参数空间中搜索全局最优解，使得代价函数快速收敛，最终在收敛条件被满足时，同化过程结束，按行政边界汇总，输出产量结果。

2.如权利要求1所述的同化叶面积指数时序曲线特征的冬小麦估产方法，其特征在于，所述步骤S1中所采用的作物模型为WOFOST模型，采集研究区内的土壤、气象和作物参数，对影像和采集的参数进行空间位置的匹配，对于不敏感模型参数直接使用默认值，对于只对产量敏感，对LAI时序曲线特征不敏感的模型参数使用反距离权重插值法对遥感影像的每个像素赋予参数值，完成参数区域化标定；采用扩展的傅里叶振幅敏感度检验法EFAST对WOFOST模型的作物参数和土壤与管理参数进行全局敏感性分析，按如下公式(1)求解一阶敏感性指数和总敏感性指数，

$S_i=\frac{V_i}{V}$ $S_{\mathrm{ij}}=\frac{V_{\mathrm{ij}}}{V}$ $S_{\mathrm{ijm}}=\frac{V_{\mathrm{ijm}}}{V}...S_{12...i...k}=\frac{V_{12...i...k}}{V}$

S_T.i＝S_i+S_ij+S_ijm+…+S_12…i…k    (1)

V(X)＝∑_iV_i+∑_i≠jV_ij+∑_i≠j≠mV_ijm+V_12...i...k

式中，i，j，m为作物模型参数的序列号，k为模型参数的总个数，V(X)为模型结果X的总方差；V_i为参数x_i的方差；V_ij～V_12...k为各参数相互作用的方差，S_i、S_ij、S_ijm为参数x_i的一阶、二阶、三阶敏感性指数，S_T.i为参数x_i的总敏感性指数。

3.如权利要求1所述的同化叶面积指数时序曲线特征的冬小麦估产方法，其特征在于，所述步骤S3中采用上包络线滤波，滤波过程按如下公式(2)：

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

式(2)中Y_j+i表示原始LAI曲线上的一块窗口内的值，m为窗口的半径，N为卷积数目，窗口的宽度为2m+1，

表示滤波后窗口中心的LAI值，C_i表示第i个LAI值的滤波系数。

4.如权利要求1所述的同化叶面积指数时序曲线特征的冬小麦估产方法，其特征在于，所述步骤S4或S5中所采用的曲线拟合为Logistic曲线拟合，曲线方程按如下公式(3)：

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

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

5.如权利要求4所述的同化叶面积指数时序曲线特征的冬小麦估产方法，其特征在于，所述步骤S4或S5中所要提取的关键特征点按如下公式(4)计算曲率：

$K=\frac{\mathrm{d\&alpha;}}{\mathrm{ds}}=\frac{b^2\mathrm{cz}(1-z){(1+z)}^3}{{[{(1+z)}^4+{\left(\mathrm{bcz}\right)}^2]}^{\frac{3}{2}}}---\left(4\right)$

其中，z＝e^a+bt，按如下公式(5)计算斜率：

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

6.如权利要求1所述的同化叶面积指数时序曲线特征的冬小麦估产方法，其特征在于，所述步骤S6中遥感观测误差表示的权重按如下公式(6)求得：

$w_i=1-\frac{v_i}{\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}{v}_i}---\left(6\right)$

其中，

$v_i=\frac{|\left(\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{LAI}}_{\mathrm{kj}}^i\right)-{\mathrm{LAI}}_{\mathrm{RS}}^i|}{\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{LAI}}_{\mathrm{kj}}^i}$

其中m是在研究区内土地单元个数，n是在每个样地内布设的采样块个数，LAIⁱ _kj是在第i个曲线特征点处第k个单元中第j个采样块内的单点测量LAI值，LAIⁱ _RS是土地单元的MODIS LAI值，w_i是三个曲线特征点处的遥感观测误差权重。

7.如权利要求6所述的同化叶面积指数时序曲线特征的冬小麦估产方法，其特征在于，所述步骤S6中采用SCE_UA算法对代价函数进行最小化，并结合所述权重，代价函数按如下公式(7)：

$J=w_1{(T_1^R-T_1^M)}^2+w_2{(T_2^R-T_2^M)}^2+w_3{(T_3^R-T_3^M)}^2---\left(7\right)$

式(7)中，T₁，T₂，T₃分别代表曲率最大值点、斜率最大值点、叶面积最大值点处的儒略日时间，加上标R代表遥感观测值、加上标M代表模型模拟值，w_i代表该时间点i处以遥感观测误差表示的权重，J为代价函数值。

8.如权利要求1所述的同化叶面积指数时序曲线特征的冬小麦估产方法，其特征在于，所述步骤S6中的收敛条件为以下条件之一：

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

(2)目标函数值在5次循环后无法提高0.0001％；

(3)计算代价函数的次数超过10000次，

WOFOST模型将按照同化获得的模型参数运行，输出产量结果。

Images