CN103955860A - 一种基于集合卡尔曼滤波同化的区域作物估产方法 - Google Patents

Abstract

本发明提供一种基于集合卡尔曼滤波同化的区域作物估产方法，融合了遥感数据和作物模型的优势，把植被遥感中普遍使用的EVI作为观测变量，LAI作为同化变量，通过集合卡尔曼滤波算法进行了对模型LAI的最优化调整，并将PROSAIL模型作为观测算子，解决了观测变量和状态变量不一致的问题，实现了遥感信息和模型的同化，避免了用反射率反演LAI带来的误差。同化EVI后获得的作物产量，与未同化相比较，均方根误差RMSE减小而决定系数R²明显上升，同化后使作物模型产量估算的精度有显著提高，产量空间分布趋势与统计产量一致。

Description

一种基于集合卡尔曼滤波同化的区域作物估产方法

技术领域

本发明属于农业遥感领域，具体涉及一种基于集合卡尔曼滤波同化的区域作物估产方法。

背景技术

遥感技术是当今的前沿技术，能够帮助快速准确地收集农业资源和农业生产的信息，结合地理信息系统和全球定位系统等其他现代高新技术，可以实现信息收集和分析的定时、定量、定位，客观性强，不受人为干扰，方便决策。目前运用遥感技术开展农业监测工作，使农业决策科学化提高到了一个新的水平，同时也为农业生产提供了高质量的服务。现在业务上已经可以实现对我国冬小麦种植面积变化、长势状况、旱情、单产和总产进行常态化评价，在冬小麦单产评估方面，主要是利用长势和旱情遥感监测结果、农学模型和气象模型估测结果进行相关分析，以及地面样方实地调查结果综合评定后得出。与本发明最接近的技术有专利“基于模型的小麦产量遥感监测预报方法（CN102162850A）”和专利“一种大面积作物种植面积及其产量的遥感检测估算方法（CN100394212C）”。

现有的遥感估产技术主要分为两类：指数相关法和同化法。指数相关法是根据遥感植被指数、气象要素等数据和产量做相关分析来进行估测，因此机制简单粗糙，难以机理性地反映产量的形成过程。同化法能够耦合遥感观测和作物模型，能够实现两者的优势互补，但目前的同化法所要同化的变量都是叶面积指数LAI，高精度的LAI数据难以获得，如果使用高分辨率的影像反演则需要丰富的技术积累和巨大的工作量，同时还会引入误差。

发明内容

为解决现有技术中同化法存在误差的问题，本发明提供一种基于集合卡尔曼滤波同化的区域作物估产方法，基于同化法进行估产，但同化变量选择了基于国产高分辨率HJ-1A/B卫星多光谱CCD数据的增强型植被指数EVI，通过把PROSAIL冠层辐射传输模型作为观测算子带入集合卡尔曼滤波进行EVI的同化,避免了反演LAI带来的误差，提高了同化的精度。

本发明提供一种基于集合卡尔曼滤波同化的区域作物估产方法，具体步骤如下：

S1：对作物模型的参数进行标定，对遥感数据与地面参数进行空间匹配，使用地统计空间插值法实现大区域的参数确定，完成数据准备；

S2：对工作区整个生育期内HJ-1A/B卫星的EVI数据按时间序列合成，对每个网格单元生成时间序列曲线；

S3：逐网格单元运行作物模型，对遥感EVI添加高斯扰动，生成观测成员集合；

S4：逐天运行作物模型，对模型的状态变量LAI添加高斯扰动生成模型成员集合；

S5：当存在遥感观测时，将观测成员集合和模型成员集合输入到集合卡尔曼滤波中，计算同化后的LAI最优值；

S6：用S5中计算得到的最优LAI值替换作物模型中当天的状态变量LAI，重复进行步骤S4至S6，直到在该像素模型运行结束，输出产量结果；

S7：重复步骤S3至S6，最终在完成所有像素上的运行后，按行政边界汇总，输出平均产量，指导作物生产。

其中，S1所述作物模型为SWAP作物模型。

其中，S1具体为：采用SWAP作物模型，采集工作区内的土壤参数、气象参数和作物参数，对遥感影像和采集的参数进行空间位置的匹配，对于不敏感的模型参数直接使用经验值，对于六种气象参数和两种积温参数使用反距离权重插值法对每个网格单元赋予参数值，完成参数区域化标定。

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

其中，所述两种积温参数为播种到出苗的积温TSUMEM和出苗到开花的积温TSUM1。

其中，S1所述作物参数为SWAP模型中描述作物生长和发育特性的参数，如作物初始干物质量TDWI、叶片在35℃时的生命期SPAN等。

其中，S2所述EVI数据是由环境系列中的A、B两颗卫星的传感器载荷（HJ-1A/B CCD）所采集数据经过FLAASH大气校正后计算得到的，大气校正过程中输入FLAASH的标准输入参数以及卫星的观测天顶角和方位角，EVI计算如公式（1）：

$\mathrm{EVI}=2.5\×\frac{{\mathrm{\ρ}}_{\mathrm{NIR}}-{\mathrm{\ρ}}_{\mathrm{RED}}}{{\mathrm{\ρ}}_{\mathrm{NIR}}+0.6{\mathrm{\ρ}}_{\mathrm{RED}}-7.5{\mathrm{\ρ}}_{\mathrm{BLUE}}+1}---\left(1\right)$

ρ_NIR为传感器载荷所测得的近红外波段反射率，ρ_RED为红光波段反射率，ρ_BLUE为蓝光波段反射率。

其中，S3、S4所述添加高斯扰动，计算如公式（2）：

S=s+aε                           （2）

其中S表示增加高斯扰动后的观测或模型成员组成的向量；s表示初始观测值或模型状态变量；a为常数；ε则表示与S同维的随机数向量，其成员符合标准正态分布，即ε_i～N(0,1)。

其中，S5所述将观测成员集合和模型成员集合输入到集合卡尔曼滤波中所采用的同化方法集合卡尔曼滤波以公式（3）（4）（5）进行计算：

B_t=HA_t+v_t                                        （3）

A^f _t=MA^a _t-1+w_t                                      （4）

A^a _t=A^f _t+K_t(B-HA^f _t)                               （5）

A_t表示t时刻下作物模型中的LAI状态变量集合；A^f _t表示LAI的预报集合；B_t为t时刻的一个观测数据集合；该时刻状态的最优估计集合为A^a _t，A^a _t的均值即为该时刻状态的最优估计值；H为观测算子；M是状态变换方程，即SWAP模型；v_t是测量噪声；w_t是过程误差；K_t是卡尔曼增益，表示观测数据的权重。

其中，S5中所要采用的集合卡尔曼滤波中所使用的观测算子是叶片辐射传输模型PROSAIL，作用是将同化变量LAI转换为观测变量EVI。PROSAIL由叶片反射率模型PROSPECT和辐射传输模型SAIL组合而成，模型假定植物冠层是由方位随机分布的水平均一及无限扩展的各向同性叶片组成的混合体，叶片均具有漫散射的发射和透射特性，其透射率和吸收系数是叶片结构参数和生物化学参数的函数。在给定太阳高度和地表植被理化、几何参数后，即可获得特定观测方向上的全波段冠层反射率，从而根据蓝光、红光和近红外波段反射率计算得到EVI。

其中，所述K_t按照如下公式（6）-（11）进行计算：

$\stackrel{\‾}{a_t^f}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{a}_{i,t}^f---\left(6\right)$

$\stackrel{\‾}{b_t}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{b}_{i,t}---\left(7\right)$

$P_t^f=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(a_{i,t}^f-\stackrel{\‾}{a_t^f})\·{(a_{i,t}^f-\stackrel{\‾}{a_t^f})}^T---\left(8\right)$

$R_t=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(b_{i,t}-\stackrel{\‾}{b_t})\·{(b_{i,t}-\stackrel{\‾}{b_t})}^T---\left(9\right)$

$K_t=P_t^f{H}^T{(H{P}_t^f{H}^T+R_t)}^{-1}---\left(10\right)$

${\mathrm{HP}}_t^f{H}^T=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}[H\left(a_{i,t}^f\right)-H\left(\stackrel{\‾}{a_t^f}\right)]\·{[H\left(a_{i,t}^f\right)-H\left(\stackrel{\‾}{a_t^f}\right)]}^T---\left(11\right)$

N表示集合大小，a^f _i,t表示t时刻模型第i个集合成员的一个预报状态，b_i,t表示t时刻观测集合的第i个成员，P_t ^f是表示t时刻的预报集合的方差，R_t表示t时刻观测集合的方差，H为观测算子PROSAIL模型，其他符号表示计算过程中的中间变量。

其中，所述作物优选冬小麦。

本发明还提供所述基于集合卡尔曼滤波同化的区域作物估产方法在指导作物生产中的应用。

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

本发明融合了遥感数据和作物模型的优势，把植被遥感中普遍使用的EVI作为观测变量，LAI作为同化变量，通过集合卡尔曼滤波算法进行了对模型LAI的最优化调整，将PROSAIL模型作为观测算子，解决了观测变量和状态变量不一致的问题，实现了遥感和模型的同化，避免了遥感反射率反演LAI误差。同化EVI后获得的作物产量，与未同化相比较，均方根误差RMSE减小而决定系数R²明显上升，同化后使作物模型产量估算的精度有显著提高，产量空间分布趋势与统计产量一致。

附图说明

图1为本发明一种基于集合卡尔曼滤波同化的区域作物估产方法的流程图。

图2为潜在水平下同化后的产量结果图。

图3为水分胁迫水平下同化后的产量结果图。

具体实施方式

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

实施例1

本发明的流程图见图1。

步骤S1，选择河北保定地区冬小麦主产区为研究区域，该区域地处115°10′E–116°20′E，38°15′N–39°40′N，地形以平原为主，耕地占总面积的60%以上，大部分地区适宜小麦生长，也是河北省冬小麦主产区。气候属温带季风性气候，年日照时数2400～3100h，年均降水量300～800mm。获取以下数据：根据研究区外包络范围，选取21个国家级气象站的日最高/最低气温、日照总辐射、水汽压、风速、降水模型所需的6个气象要素；研究区内农业气象试验站采集的土壤参数和作物参数；经纬度、高程等控制参数；农业气象数据和09年河北省分县冬小麦产量数据；冬小麦关键生育期3月至6月的时间序列HJ-1A/B CCD数据，将含有地理位置信息的数据统一坐标，完成数据的空间匹配。

利用研究区的遥感影像和野外调查数据解译得到冬小麦种植区域，生成1公里的网格，计算网格单元内冬小麦种植百分比，设定阈值剔除冬小麦种植比例低于20%的单元。对模型的产量输出和气象数据统一度量，根据荷兰瓦赫宁根大学开发的SWAP作物模型在河北地区冬小麦的初步标定结果，将气象数据和作物参数反距离插值生成1公里每像素的栅格数据。

步骤S2，将收集的作物生育期内的HJ-1A/B CCD数据进行FLAASH大气校正，然后利用红光和近红外波段的反射率计算EVI，大气校正过程中输入FLAASH的标准输入参数以及卫星的观测天顶角和方位角，EVI计算如公式（1）：

$\mathrm{EVI}=2.5\×\frac{{\mathrm{\ρ}}_{\mathrm{NIR}}-{\mathrm{\ρ}}_{\mathrm{RED}}}{{\mathrm{\ρ}}_{\mathrm{NIR}}+0.6{\mathrm{\ρ}}_{\mathrm{RED}}-7.5{\mathrm{\ρ}}_{\mathrm{BLUE}}+1}---\left(1\right)$

ρ_NIR为传感器载荷所测得的近红外波段反射率，ρ_RED为红光波段反射率，ρ_BLUE为蓝光波段反射率。

根据解译得到冬小麦种植区域制作掩膜，对处于掩膜内的EVI数值按照1km网格单元计算平均值，得到升尺度后的1km分辨率EVI影像，对整个生育期内的EVI数据进行同样的操作最后按时间序列叠加合成，对每个网格单元生成该点的EVI时间序列曲线。

步骤S3，在模型运行过程中，如果当天有遥感观测则进行同化，此时对该天的HJ-1A/B EVI添加高斯扰动，计算如公式（2）：

S=s+aε                            （2）

其中S表示增加高斯扰动后的观测或模型成员组成的向量；s表示初始观测值或模型状态变量；a为常数；ε则表示与S同维的随机数向量，其成员符合标准正态分布，即ε_i～N(0,1)。

步骤S4，将该网格单元的土壤、气候、作物数据和驱动数据输入SWAP模型，从冬小麦播种期为起点开始运行。SWAP模型运行到有观测数据的日期时，对模型模拟的LAI进行扰动，生成1个预报集合MLAI{mLAI₁,mLAI₂...mLAI_N}。

步骤S5，当存在遥感观测时，将观测成员集合和模型成员集合输入到集合卡尔曼滤波中，计算同化后的LAI最优值。集合卡尔曼滤波应用以公式（3）（4）（5）进行计算：

B_t=HA_t+v_t                                        （3）

A^f _t=MA^a _t-1+w_t                                     （4）

A^a _t=A^f _t+K_t(B-HA^f _t)                               （5）

A_t表示t时刻下作物模型中的LAI状态变量集合；A^f _t表示LAI的预报集合；B_t为t时刻的一个观测数据集合；该时刻状态的最优估计集合为A^a _t，A^a _t的均值即为该时刻状态的最优估计值；H为观测算子；M是状态变换方程，即SWAP模型；v_t是测量噪声；w_t是过程误差；K_t是卡尔曼增益，表示观测数据的权重。

其中，K_t是按照公式（6）-（11）进行计算：

$\stackrel{\‾}{a_t^f}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{a}_{i,t}^f---\left(6\right)$

$\stackrel{\‾}{b_t}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{b}_{i,t}---\left(7\right)$

$P_t^f=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(a_{i,t}^f-\stackrel{\‾}{a_t^f})\·{(a_{i,t}^f-\stackrel{\‾}{a_t^f})}^T---\left(8\right)$

$R_t=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(b_{i,t}-\stackrel{\‾}{b_t})\·{(b_{i,t}-\stackrel{\‾}{b_t})}^T---\left(9\right)$

$K_t=P_t^f{H}^T{(H{P}_t^f{H}^T+R_t)}^{-1}---\left(10\right)$

${\mathrm{HP}}_t^f{H}^T=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}[H\left(a_{i,t}^f\right)-H\left(\stackrel{\‾}{a_t^f}\right)]\·{[H\left(a_{i,t}^f\right)-H\left(\stackrel{\‾}{a_t^f}\right)]}^T---\left(11\right)$

N表示集合大小，a^f _i,t表示t时刻模型第i个集合成员的一个预报状态，b_i,t表示t时刻观测集合的第i个成员，P_t ^f是表示t时刻的预报集合的方差，R_t表示t时刻观测集合的方差，H为观测算子PROSAIL模型，其他符号表示计算过程中的中间变量。

集合卡尔曼滤波中所使用的观测算子是冠层辐射传输模型PROSAIL，作用是将同化变量LAI转换为观测变量EVI。PROSAIL由叶片反射率模型PROSPECT和辐射传输模型SAIL组合而成，模型假定植物冠层是由方位随机分布的水平均一及无限扩展的各向同性叶片组成的混合体，叶片均具有漫散射的发射和透射特性，其透射率和吸收系数是叶片结构参数和生物化学参数的函数。在给定太阳高度和地表植被理化、几何参数后，即可获得特定观测方向上的全波段冠层反射率，从而根据蓝光、红光和近红外波段反射率计算得到EVI。

步骤S6，t时刻的LAI集合的均值，即为该时刻的LAI最优估计值，用同化得到的最优LAI替换SWAP模型中的LAI，进行下一时刻的运行。如果没有观测数据，则SWAP模型直接运行到下一时刻。重复上述过程直到冬小麦成熟期。在同化过程中，每个集合成员都是独立运行。

步骤S7，重复步骤S4至S6，最终运行完成所有冬小麦像素网格单元后，按行政边界汇总，输出县域冬小麦单产。

按行政区划计算该县域冬小麦单产，输出区域上的产量结果见图2、图3，图2为潜在水平下同化后的产量结果图，图3为水分胁迫水平下同化后的产量结果图。

本发明实施例所述的一种基于集合卡尔曼滤波同化的区域冬小麦估产方法，融合了遥感数据和作物模型的优势，把植被遥感中普遍使用的EVI作为观测变量，LAI作为同化变量，通过集合卡尔曼滤波算法进行了对模型LAI的最优化调整，将PROSAIL模型作为观测算子，解决了观测变量和状态变量不一致的问题，实现了遥感和模型的同化，避免了遥感反射率反演LAI误差。同化EVI后获得的冬小麦产量，与未同化相比较，均方根误差RMSE从2168kg ha^-1减小到775kg ha^-1，决定系数R²从0.25增加到0.51，而在水分胁迫水平下的RMSE和R²分别为3005kg ha^-1和0.53。结果表明，同化后使作物模型产量估算的精度有显著提高，产量空间分布趋势与统计产量一致。

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

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

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

Claims (9)

1.一种基于集合卡尔曼滤波同化的区域作物估产方法，其特征在于，具体步骤如下：

S1：对作物模型的参数进行标定，对遥感数据与地面参数进行空间匹配，使用地统计空间插值法实现大区域的参数确定，完成数据准备；

S2：对工作区整个生育期内HJ-1A/B卫星的EVI数据按时间序列合成，对每个网格单元生成时间序列曲线；

S3：逐网格单元运行作物模型，对遥感EVI添加高斯扰动，生成观测成员集合；

S4：逐天运行作物模型，对模型的状态变量LAI添加高斯扰动生成模型成员集合；

S5：当存在遥感观测时，将观测成员集合和模型成员集合输入到集合卡尔曼滤波中，计算同化后的LAI最优值；

S6：用S5中计算得到的最优LAI值替换作物模型中当天的状态变量LAI，重复进行步骤S4至S6，直到在该像素模型运行结束，输出产量结果；

S7：重复步骤S3至S6，最终在完成所有像素上的运行后，按行政边界汇总，输出平均产量，指导作物生产。

2.权利要求1所述基于集合卡尔曼滤波同化的区域作物估产方法，其特征在于，S1所述作物模型为SWAP作物模型。

3.权利要求2所述基于集合卡尔曼滤波同化的区域作物估产方法，其特征在于，S1具体为：采用SWAP作物模型，采集工作区内的土壤参数、气象参数和作物参数，对遥感影像和采集的参数进行空间位置的匹配，对于不敏感的模型参数直接使用经验值，对于六种气象参数和两种积温参数使用反距离权重插值法对每个网格单元赋予参数值，完成参数区域化标定。

4.权利要求1所述基于集合卡尔曼滤波同化的区域作物估产方法，其特征在于，S2所述EVI数据是由环境系列中的A、B两颗卫星的传感器载荷（HJ-1A/B CCD）所采集数据经过FLAASH大气校正后计算得到的，大气校正过程中输入FLAASH的标准输入参数以及卫星的观测天顶角和方位角，EVI计算如公式（1）：

$\mathrm{EVI}=2.5\×\frac{{\mathrm{\ρ}}_{\mathrm{NIR}}-{\mathrm{\ρ}}_{\mathrm{RED}}}{{\mathrm{\ρ}}_{\mathrm{NIR}}+0.6{\mathrm{\ρ}}_{\mathrm{RED}}-7.5{\mathrm{\ρ}}_{\mathrm{BLUE}}+1}---\left(1\right)$

ρ_NIR为传感器载荷所测得的近红外波段反射率，ρ_RED为红光波段反射率，_ρBLUE为蓝光波段反射率。

5.权利要求1所述基于集合卡尔曼滤波同化的区域作物估产方法，其特征在于，S3、S4所述添加高斯扰动，计算如公式（2）：

S=s+aε                             （2）

其中S表示增加高斯扰动后的观测或模型成员组成的向量；s表示初始观测值或模型状态变量；a为常数；ε则表示与S同维的随机数向量，其成员符合标准正态分布，即ε_i～N(0,1)。

6.权利要求1所述基于集合卡尔曼滤波同化的区域作物估产方法，其特征在于，S5所述将观测成员集合和模型成员集合输入到集合卡尔曼滤波中所采用的同化方法集合卡尔曼滤波以公式（3）（4）（5）进行计算：

B_t=HA_t+v_t                                        （3）

A^f _t=MA^a _t-1+w_t                                     （4）

A^a _t=A^f _t+K_t(B-HA^f _t)                               （5）

A_t表示t时刻下作物模型中的LAI状态变量集合；A^f _t表示LAI的预报集合；B_t为t时刻的一个观测数据集合；该时刻状态的最优估计集合为A^a _t，A^a _t的均值即为该时刻状态的最优估计值；H为观测算子；M是状态变换方程，即SWAP模型；v_t是测量噪声；w_t是过程误差；K_t是卡尔曼增益，表示观测数据的权重。

7.权利要求6所述基于集合卡尔曼滤波同化的区域作物估产方法，其特征在于，所述K_t按照如下公式（6）-（11）进行计算：

$\stackrel{\‾}{a_t^f}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{a}_{i,t}^f---\left(6\right)$

$\stackrel{\‾}{b_t}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{b}_{i,t}---\left(7\right)$

$P_t^f=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(a_{i,t}^f-\stackrel{\‾}{a_t^f})\·{(a_{i,t}^f-\stackrel{\‾}{a_t^f})}^T---\left(8\right)$

$R_t=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(b_{i,t}-\stackrel{\‾}{b_t})\·{(b_{i,t}-\stackrel{\‾}{b_t})}^T---\left(9\right)$

$K_t=P_t^f{H}^T{(H{P}_t^f{H}^T+R_t)}^{-1}---\left(10\right)$

${\mathrm{HP}}_t^f{H}^T=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}[H\left(a_{i,t}^f\right)-H\left(\stackrel{\‾}{a_t^f}\right)]\·{[H\left(a_{i,t}^f\right)-H\left(\stackrel{\‾}{a_t^f}\right)]}^T---\left(11\right)$

N表示集合大小，a^f _i,t表示t时刻模型第i个集合成员的一个预报状态，b_i,t表示t时刻观测集合的第i个成员，P_t ^f是表示t时刻的预报集合的方差，R_t表示t时刻观测集合的方差，H为观测算子PROSAIL模型，其他符号表示计算过程中的中间变量。

8.权利要求1所述基于集合卡尔曼滤波同化的区域作物估产方法，其特征在于，所述作物为冬小麦。

9.权利要求1-7任一项所述基于集合卡尔曼滤波同化的区域作物估产方法在指导作物生产中的应用。