CN102592181A - 一种作物种植面积统计数据空间分布的优化方法 - Google Patents

Abstract

本发明提供了一种作物种植面积统计数据空间分布的优化方法，包括以下步骤：A、收集作物遥感影像及行政单元内作物种植面积统计数据，并预处理作物遥感影像；B、在所述作物遥感影像上选择目标作物并建立一定区域内目标作物指数时间序列特征标准曲线；C、利用所述作物指数时间序列特征标准曲线确定作物物候信息，并根据该作物物候信息确定作物分类算法模型；D、利用全局优化算法并根据所述区域内作物实际种植面积农业统计数据，计算作物分类算法模型的特征参数的最优值；E、将所述特征参数的最优值代入所述作物分类算法模型求出作物种植面积空间分布信息。从而实现作物统计数据空间化信息精度的进一步提高。

Description

一种作物种植面积统计数据空间分布的优化方法

技术领域

本发明涉及一种作物种植面积统计数据空间分布的优化方法。

背景技术

作物种植面积统计数据是农业统计数据的核心指标之一。准确的作物种植面积统计数据不仅关系到国家农业政策的正确制定，而且对提高农业生产管理水平、增强农业资源监测能力、优化作物种植结构、准确预测作物产量、确保国家粮食安全具有重要意义。近年来，随着世界各国对全球气候变化、人口、资源、环境、生态和粮食安全等一系列问题定量研究的逐步深入，准确的农业统计数据在上述问题研究中发挥了重要作用。其中，作物种植面积统计数据成为不可或缺的基础数据之一。

然而，包括作物种植面积统计数据在内的传统农业统计大多以各级行政单元为统计单位，统计数据除仅仅反映行政单元内总量水平或平均水平外，数据存在统计单元不统一、定位不准确和空间分辨率低等缺点，尚未真正揭示统计数据在地理空间上的真实分布和空间变异情况，更不利于统计指标与其他自然、生态要素空间数据(如土地利用/土地覆盖数据、生态环境背景数据)实现交叉使用和叠置分析等功能。传统的区域农业统计指标空间化的方法是直接给区域行政单元赋以统计数值，然后生成图斑表现形式的专题地图，而图斑的边界就是行政单元的边界。通过这种方法获得的统计指标空间分布信息具有同一区域内指标均匀分布，不同区域间指标存在突变分布的特点，这与农业统计数据空间分布的真实状态不一致。此外，这种方法形成的空间分布结果难以满足跨行政边界的农业数据统计需要以及与自然、生态要素空间数据的交叉使用和叠置分析。可见，传统的农业统计数据空间化方法已经很难满足社会经济和科学技术发展的应用要求。

近年来，随着遥感、地理信息系统、全球定位系统和计算机技术的发展，以行政单元为空间表达方式的各类统计数据空间化已经成为空间信息技术、全球变化等社会交叉学科领域的研究热点，即：将基于社会行政单元(如省、市、县、乡镇等)的统计信息转化为基于一定尺寸网格单元的空间信息，实现在空间网格统计单元上与土地利用、生态环境背景数据等自然要素数据进行综合分析，真正实现自然、生态和人文要素的整合。目前，国内外社会经济统计数据空间化研究主要集中在人口统计数据空间化、国内生产总值(GDP)数据空间化以及空间化方法的改进等。统计数据空间化采用的方法主要包括面插值法和统计模型法，主要思想大都采用多源自然数据和社会经济数据融合的思路。近年来，社会经济统计数据空间化的研究内容和范围仍在不断拓展，如农业产值空间化、区域化肥投入空间化、作物产量空间化研究等。然而，与作物种植面积统计数据空间化研究相关的报道则很少，现有的这方面研究如历史耕地面积统计数据网格化、基于交叉信息熵的作物面积统计数据空间化方法等。上述作物或耕地面积空间化方法继续延续了基于自然、人文、社会经济等多源信息融合思想，对属性数据空间化研究具有开拓性意义，可为大尺度宏观研究提供很好的空间模拟信息，特别是对宏观研究中升尺度方法的成功利用发挥了重要作用。但上述作物种植面积统计数据空间化模型研究并未从社会经济统计作物自身生长变化规律出发进行研究。因此，导致空间化后的统计信息在统计数据自身尺度以下微观信息的准确性大大降低；此外，上述方法均采用集成遥感技术提取的土地利用、耕地分布或作物分布等作为辅助信息实现社会数据的空间化，但并未对作物自身遥感信息(特别是时序遥感信息)加以直接应用，在一定程度上影响了作物统计数据空间化信息精度的进一步提高。

发明内容

有鉴于此，本发明的主要目的在于提供一种作物种植面积统计数据空间分布优化的方法，以实现作物统计数据空间化信息精度的进一步提高。

本发明提供了一种作物种植面积统计数据空间分布的优化方法，包括以下步骤：

A、收集作物遥感影像及行政单元内作物实际种植面积农业统计数据，并预处理作物遥感影像；

B、在所述作物遥感影像上选择目标作物并建立一定区域内目标作物指数时间序列特征标准曲线；

C、利用所述作物指数时间序列特征曲线确定作物物候信息，并根据该作物物候信息确定作物分类算法模型；

D、利用全局优化算法并根据所述区域内作物实际种植面积农业统计数据，计算作物分类算法模型的特征参数的最优值；

E、将所述特征参数的最优值代入所述作物分类算法模型求出作物种植面积空间分布信息。

由上可以看出，本发明方法从作物自身规律出发，利用作物物候信息和作物生长特征曲线实现了作物统计数据空间化，该方法大大提高了统计数据单元内微观尺度(像元尺度)上作物分布信息的准确性，从而增强了作物分布空间化信息在模拟单元像元尺度上的可用性；其次本发明方法直接并充分地利用与作物直接相关的时间序列动态遥感信息进行作物空间化信息的提取，比传统方法利用单一时间点的辅助数据或历史数据具有较大改进，同时，该方法侧重于大范围覆盖的遥感信息的直接应用，利用该方法可在更大范围内应用和获取相关数据。

在上述方法中，步骤A中所述遥感影像的收集包括：

收集高时间分辨率的最大值合成的一定时期时间序列数据。

由上可以看出，上述遥感数据利于与作物生长规律相结合，更利于作物空间化信息的提取。

在上述方法中，步骤A中所述作物遥感影像的预处理包括：

对所述作物遥感影像进行影像拼接、投影转换和裁切，并将所有所述作物遥感影像按照时间先后顺序排序后合成时间序列影像；

采用S-G滤波平滑法对时间序列影像进行去云平滑处理。

由上可以看出，采用S-G滤波平滑法可有效去除多时相遥感数据受云、气溶胶影响造成的噪音残留，从而进一步提高基础数据质量。

在上述方法中，所述步骤B包括：

作物生长期内地面样方调查；

利用差分GPS对调查样区中心点位置进行定位；

根据地面调查样方中心点位置结果，在时间序列的合成影像上确定作物种植地块中心位置所对应的像元，并获取地面调查样方的作物指数时间序列特征曲线；

采用植被指数的上包络线来拟合每个作物地面调查样方的作物指数时间序列特征曲线；

对上述作物地面调查样方所对应的作物指数特征曲线进行平均，获得一定区域内的作物指数时间序列特征标准曲线。

在上述方法中，步骤C所述作物物候信息的确定为：

根据所述时间序列影像中的几个关键时期的作物指数值的变化趋势，提取出作物物候信息。

由上可以看出，只根据所述时间序列影像中的几个关键时期的作物指数值的变化趋势提取作物物候信息，可缩短作物物候信息提取时间并减少非关键时期冗余数据对作物物候信息提取的干扰。

在上述方法中，其特征在于，步骤C所述作物分类算法模型的确定为：

根据所述作物物候信息反映出的作物指数时间序列规律，利用指数阈值法确定作物分类算法模型。

由上可以看出，该方法适合所有依靠阈值算法控制作物分类结果的作物面积提取方法，如决策树等。此外，该方法克服了传统作物面积提取中分类阈值自动优化提取的难题，进而增强了该方法的自动化程度以及在更大范围内应用的可能性和可操作性。

在上述方法中，所述作物指数时间序列为作物植被指数时间序列。

在上述方法中，所述步骤D包括以下子步骤：

D1、基于时间序列影像和所述作物分类算法模型进行作物分类，模拟出作物种植面积；

D2、由上述作物分类算法模型统计出所述区域内作物的模拟种植面积；

D3、建立所述区域内作物的模拟种植面积与该区域内作物实际种植面积农业统计数据之间的代价函数；

D4、当临近特定数量个代价函数值之差的绝对值小于特定值时，确定与最小代价函数值相对应的作物分类算法模型特征参数的初始值为特征参数的最优值。

在上述方法中，步骤D4所述特定数量为5，所述特定值为0.001。

在上述方法中，所述代价函数如下：

$y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{sim}}-X_{\mathrm{obs}})}^2$

其中，n是所述区域内作物实际种植面积统计数据的个数，X_sim是所述作物分类算法模型模拟出的所述区域内作物种植面积，X_obs是所述区域内作物实际种植面积的农业统计数据。

附图说明

图1为本发明提供的作物种植面积统计数据空间分布优化方法框架图；

图2为本发明提供的作物种植面积统计数据空间分布优化方法的研究区位置和调查样区点位示意图；

图3为本发明提供的作物种植面积统计数据空间分布优化方法的具体实施流程图；

图4为研究区域内冬小麦从播种到成熟期间MODIS-NDVI时间序列数据的平滑效果图；

图5为研究区域内冬小麦NDVI时间序列特征曲线图；

图6a为通过本发明方法获得的整个研究区域内的冬小麦种植面积统计数据空间化结果示意图；

图6b为通过高分辨率遥感影像获取的冬小麦分布数据示意图。

具体实施方式

本发明是一种在高时间分辨率遥感植被指数构成的时序数据支持下展开作物面积统计数据空间分布优化方法。由于植被生长以年为周期，在这个生长周期内不同植被类型有着各自生繁衰枯的物候节律，表现出不同的生长规律，而且规律性极强。这种规律性可以作为植被分类的出发点。作物生长周期内具有其独特的生长规律(如播种期、出苗期、起身期、拔节、抽穗、乳熟、收获期等)，但不同作物的叶子内部构造存在一定的差异。因此，不同作物的物候期会存在一定差异，导致不同作物的波谱形态时间曲线会有所不同，但同一种作物在同一地区具有相对稳定的作物生长发育规律和NDVI-时间变化曲线(NDVI，Normalized Difference Vegetation Index，归一化植被指数，是一个能够通过遥感获取且能够直接、有效地反映作物绿度长势、盖度、生物量和作物单产的最佳参考量之一)。因此，本发明采用多时相植被指数遥感信息，根据绿色植物独特的波谱反射特征，将植被(作物)与其他地物区分开，以便进行作物空间分布信息提取。

所述优化方法的主要框架如图1所示，包括：

首先，利用平滑过滤技术提高时序植被指数自身质量，以便消除云、气溶胶及传感器不稳定等因素产生的噪音对数据质量的影响；然后，在地面调查数据支持下，通过遥感影像获取作物时序植被指数构成的标准时序特征曲线；再利用该时序特征曲线反映的作物物候信息，通过NDVI阈值法确定提取作物分类算法模型，并指定阈值初始值和变化范围；在优化算法支持下，以作物区域统计面积数据作为优化比较对象，当模拟的区域作物面积与区域实际作物面积农业统计数据间差异达到最小时，得到作物分类模型最优参数；最后，在最优参数驱动下，利用作物分类算法模型提取作物空间分布信息，从而实现作物面积统计数据的空间分布优化。

下面，以位于中国北方粮食生产基地黄淮海平原内河北省石家庄地区13个地处平原的县(市)作为研究区域(如图2所示，37.60°-38.77°N，113.81°-115.35°E，覆盖面积8230km²)，详细介绍本发明提供的作物种植面积统计数据空间分布优化方法在该研究区域内的具体实施过程。由于山区县市地形复杂，作物种植面积破碎，如采用中低分辨率遥感影像，则很难满足山区或半山区作物种植面积提取要求。因此，在本实施例中仅考虑了研究区域内的部分平原县市。该区域属温带半湿润季风气候，≥0℃年积温4200℃～5500℃，年累计辐射量5.0*10⁶～5.2*10⁶KJ·m^-2，无霜期170～220d，年均降水量500～900mm，主要粮食作物为冬小麦、夏玉米，一年两熟轮作制度。研究区域内冬小麦集中种植时间为上年9月下旬至10月上旬，11月下旬至12月上旬达到冬前分蘖峰值，12月中旬开始停止生长，第二年返青时间为下年2月下旬至3月上旬，起身期为3月中下旬，拔节期为4月上旬至4月中旬，孕穗期为4月下旬，抽穗期为5月上旬，开花期为5月中旬，乳熟期为5月下旬，成熟期为6月上旬。其中，NDVI通常在孕穗-抽穗期达到最大值。

以图1示出的所述作物种植面积统计数据空间分布优化方法框架为基础，结合图3具体介绍本发明提供的一种作物种植面积统计数据空间分布优化方法。其中，所述方法包括以下步骤：

第一步骤：收集作物遥感影像及行政单元内作物种植面积统计数据，并预处理高时间分辨率作物遥感影像。

在本实施例中，高时间分辨率作物遥感影像为16天最大值合成的250米MODIS-NDVI时间序列数据(MODIS数据由terra和aqua卫星的一个重要的传感器实时观测得到的，由terra和aqua卫星通过x波段向全世界免费实时广播发送)。MODIS数据具有波段多(36个)、中低空间分辨率(250～1000米)和时间分辨率高(1～2天)等特点，非常利于与作物的生长规律相结合，因此被广泛应用于农作物动态监测业务之中，如作物面积提取、作物估产等。而NDVI的计算公式为

其中R_n为近红外波段的反射率，R_r为红光波段的反射率。所述遥感影像的收集过程为：从EOS Data Gateway中下载2008年第273天(9月29日)到2009年第177天(6月26日)期间每16天合成的MODIS-NDVI数据。所述遥感影像的预处理包括：对所述遥感影像进行影像拼接、投影转换和裁切，并将所有遥感影像按照时间先后顺序排序后合成时间序列影像；接着，采用S-G滤波平滑法(Savitzky-Golay)对时间序列影像进行去云平滑处理。

在上述预处理过程中，可采用公知的MRT(MODIS ReprojectionTool，MODIS投影变换工具)工具将所述遥感影像从SIN投影转换成Albers等面积投影，中央经线105°E，第一条纬线25°N，第二条纬线47°N。同时，为了方便使用数据，通过将下载的上述MODIS-NDVI数据缩小10000倍，将其缩小至-1～1正常范围。关于MRT工具的使用方法在此不再赘述。

而S-G滤波平滑法是利用最小二乘积拟合方法来平滑和计算一组相邻值或光谱导数，可简单理解为是一种权重滑动平均滤波，其权重取决于滤波窗口范围内做最小二乘拟合的多项式次数。该滤波平滑法可以应用于任何相同时间间隔、连续且具有一定平滑特征的数据。由于植被NDVI时间序列具有上述特征。因此，可利用S-G滤波平滑法对上述NDVI时间序列影像进行平滑处理。平滑时，采用NDVI的上包络线来拟合NDVI时间序列的变化趋势，通过迭代使S-G平滑达到最好的效果。图4示出了上述研究区域内冬小麦从播种到成熟期间MODIS-NDVI时间序列数据的平滑效果，其中，横坐标格式为年份+每16天数据的起始天。

第二步骤：在地面调查数据(作物野外地面样方调查信息)的支持下，在所述作物遥感影像上选择目标作物并建立目标作物植被指数时间序列特征曲线。

为了获取准确的作物植被指数时序特征曲线，在2009年3月开展了作物返青-拔节期(生长期)地面样方调查。其中，在石家庄地区13个县选择了45个冬小麦地面样方进行观测。考虑到冬小麦长势、种植结构的代表性以及调查点分布的均匀性，样区间隔不小于5km，每个样区面积不小于500m*500m，利用差分GPS对调查样区中心点位置进行精确定位。

然后，根据野外地面调查样方中心点位置结果，在MODIS-NDVI时间序列的合成影像上确定冬小麦种植地块中心位置所对应的像元，并获取每个野外地面调查样方冬小麦NDVI时间序列特征曲线。最后，求多个冬小麦野外地面样方所对应的NDVI特征曲线值的平均值，从而获得研究区域内具有一定代表性的冬小麦NDVI时间序列特征曲线。如图5所示，其中横轴为时间序列图像的代号，其格式为年份+每16天数据的起始天，从2008年第273天至2009年第177天，每隔16天一景合成NDVI数据，共18张NDVI影像。纵轴为拉伸后的-1～1范围内的NDVI值。

第三步骤：根据所述作物植被指数时间序列特征曲线反映出的作物物候信息，通过NDVI阈值法确定作物分类算法模型(作物面积提取模型)，并指定待优化参数(阈值)的初始值和变化范围。

根据图5示出的研究区域内冬小麦NDVI时间序列特征曲线以及研究区域内多年冬小麦物候可知，黄淮海地区9月下旬至10月上旬开始播种冬小麦，10月上中旬开始出苗。随着冬小麦的出苗生长，NDVI值开始升高，11月下旬至12月上旬冬前分蘖达到峰值，然后开始逐步停止生长。第二年2月下旬冬小麦开始返青，NDVI值逐步升高，4月上中旬开始拔节，4月下旬开始孕穗，5月上旬抽穗，其中4月下旬至5月上旬NDVI值达到峰值。然后，随着生育期的延续，冬小麦从抽穗、灌浆到乳熟，NDVI值逐步下降，一直到6月上旬冬小麦成熟，NDVI值出现大幅度的降低。因此，利用MODIS影像的几个关键时期的NDVI值的变化趋势，可以提取出冬小麦物候信息。假设用NDVI(r)表示不同时相的NDVI值，r为图5中18张时间序列图像中的某一张，取值范围为1～18。针对每个像元找出作物光谱曲线中两个波峰所对应的时相号r1、r2，并利用NDVI阈值法(公知技术，在此不再赘述)确定下列作物分类算法模型：

a、由冬小麦出苗期的NDVI时间序列规律可知，满足判别函数NDVI(r1)-NDVI(2)＞T1，NDVI(r1)＞NDVI(1)的像元为冬小麦；

b、由冬小麦刚进入返青期的NDVI时间序列规律可知，满足判别函数NDVI(12)-NDVI(10)＞T2，NDVI(11)＞NDVI(10)的像元为冬小麦；

c、由冬小麦收获期的NDVI时间序列规律可知，满足判别函数NDVI(r2)-NDVI(r2+2)＞T3的像元为冬小麦；

d、由冬小麦NDVI时间序列双峰波形规律可知，满足判别函数NDVI(r2)-NDVI(r1)＞T4的像元为冬小麦；

其中，T1、T2、T3和T4为作物分类算法模型的待优化的特征参数，其初始值为0～1内的任意值，其变化范围可以从MODIS时间序列的影像上获取，即变化范围为0～1。

第四步骤：在全局优化算法支持下，以作物种植面积统计数据作为外部优化数据比较对象，当模拟的区域作物面积与区域实际作物面积农业统计数据间差异达到最小时，得到作物分类模型的最优特征参数。

在本步骤中，利用优化效率较高、运算速度快且对参数初始值设置不敏感的SCE-UA优化算法(Shuffled Complex Evolution-Universityof Arizona，复合形混合演化算法)对作物分类算法模型的待优化的特征参数T1、T2、T3和T4进行优化。具体优化过程如下：

i、基于MODIS-NDVI时间序列数据和上述作物分类算法模型进行冬小麦分类(模拟出冬小麦种植面积)；

ii、由上述作物分类算法模型统计出研究区域内各县冬小麦的模拟种植面积；

iii、建立研究区域内各区县冬小麦的模拟种植面积与研究区域内各区县冬小麦的实际种植面积统计数据之间的代价函数如下：

$y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{sim}}-X_{\mathrm{obs}})}^2$

其中，n是研究区域内各县实际作物面积统计数据的个数(在本实施例中n＝1)，X_sim是上述作物分类算法模型模拟出的研究区域内各县冬小麦种植面积，X_obs是研究区域内各县冬小麦实际种植面积农业统计数据；

iv、当模拟的冬小麦种植面积与各区县冬小麦的实际种植面积农业统计数据间差异达到最小时，即临近5个最优目标函数值之差的绝对值小于0.001时，认为与最小目标函数值相对应的T1、T2、T3、T4的初始值为T1、T2、T3和T4的最优值，并退出优化过程；否则，调整特征参数T1、T2、T3和T4，并返回步骤i；当计算目标函数(上述代价函数)的次数超过10000次或待优化特征参数的值收缩到预定的较小值域内时，退出优化过程，优化失败。

第五步骤：在最优特征参数驱动下，利用作物分类模型提取作物空间分布信息，从而实现作物面积统计数据空间分布的优化。

在本步骤中，将特征参数T1、T2、T3和T4的最优值代入上述作物分类算法模型，由此得到最优的冬小麦面积统计数据空间分布。

最终，通过上述方法步骤分别提取了13个县市的冬小麦空间分布结果，进而将这13个冬小麦空间分布结果进行拼接，从而获得了整个研究区域内的冬小麦种植面积统计数据空间化结果(如图6a所示)。通过与高分辨率遥感影像获取的冬小麦分布数据(如图6b所示)对比可知，利用本发明提供的作物种植面积统计数据空间分布优化方法获得的冬小麦面积统计数据空间分布优化结果与高分辨率冬小麦分布数据整体分布规律一致。

为了进一步验证利用本发明方法获得的冬小麦种植面积统计数据空间化结果的精度，可利用冬小麦种植面积统计数据和高分辨率冬小麦分布图对冬小麦面积统计数据空间化分布结果进行精度分析，其中包括作物面积提取区域精度验证和作物面积提取空间位置精度验证。所述区域精度即一定大小的行政单元内冬小麦种植面积的总量精度；所述空间位置精度即研究区域内像元识别正确的百分比，也称像元精度。

在所述区域精度验证中，将各行政单元作物面积统计数据作为冬小麦种植面积总量A₀基准值，将基于时间序列数据NDVI获得的作物种植面积空间分布优化结果统计得到的冬小麦种植面积总量A与A₀进行对比，以获得各个县市行政单元内区域总量精度k_r如下：

$k_r=1-\frac{|A-A_0|}{A_0}\×100\%$

此外，为了表明在本实施例中选取的研究区域内总的区域精度，在所述区域精度验证中还计算了研究区域内所有县市的平均精度如下：

$K_r=1-\frac{|\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{A}_i-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{A}_{0i}|}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{A}_{0i}}\×100\%$

其中，k_r为研究区域范围内作物面积统计数据空间化结果的总量精度，A_i为研究区域内各县冬小麦面积分布优化总量结果，A_0i为研究区域内各县作物面积统计数据总量基准值，n为研究区域内县市的总个数。

利用上述公式在获得冬小麦种植面积统计数据空间化结果的基础上，结合ArcGIS软件的使用按照县级行政单元对冬小麦空间化分布结果进行统计，得到各县内提取冬小麦面积对应的像元数量。然后将提取冬小麦面积对应的像元数量与冬小麦种植面积统计数据折合的像元数量进行比较，即得到各县研究区域的总量精度。同理，通过比较研究区域内各县的冬小麦统计面积像元数量和作物面积提取像元数量，最终获得研究区域内作物面积统计数据空间化总量精度。通过计算可知，通过上述本发明方法获得的研究区域冬小麦种植面积统计数据的总量精度为99.70％，如表1所示，研究区域内除石家庄市冬小麦提取面积总体精度为92.65％外，其他地区的精度都在98％以上或者接近100％。石家庄冬小麦提取面积总体精度略低的主要原因是由于城市周边地区耕地较少，冬小麦种植地块破碎，250米分辨率的MODIS数据不能完全满足地形、种植结构复杂或地块破碎地区高精度提取冬小麦种植面积的要求。

表1.冬小麦种植面积统计数据空间化结果总量精度评价

由上可见，上述本发明方法可获得较高精度的作物面积统计数据空间化结果，且各县作物面积统计数据空间化结果与区域空间化结果的精度均能够满足统计数据空间化单元内容量不变性的基本原则和精度要求。

在所述空间位置精度验证中，将10米ALOS(Advanced LandObserving Satellite，陆地观测卫星)高分辨率影像获得的冬小麦分布数据作为空间化结果验证的“位置真值”。因为不同分辨率数据间由于尺度问题导致数据结果间的对比性不强，故而将10米空间分辨率的冬小麦分布数据转化为与MODIS遥感数据一致的空间分辨率(250米)，从而便于开展逐像元比较法验证作物面积统计数据空间化结果的位置精度。在本实施例中，空间位置精度采用Kappa分析法获得作物分布空间化结果与作物分布影像间吻合度或精度指标K_hat，其公式为：

$K_{\mathrm{hat}}=\frac{N\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{x}_{\mathrm{ii}}-\underset{i=1}{\overset{r}{\mathrm{\Σ}}}\left(x_{i+}{x}_{+i}\right)}{N^2-\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{x}_{i+}{x}_{+i}}$

其中，r是错误矩阵中的总列数(即总的类别数)，x_ii是错误矩阵中第i行、第i列上像元数量(即正确分类的数目)，x_i+和x_+i分别是第i行和第i列的总像元数量，N是总的用于精度验证的像元数量。

如表2所示，利用上述公式将MODIS-NDVI时间序列影像获得的冬小麦统计数据空间化结果与高分辨率影像获得的冬小麦分布数据进行对比可知，研究区域内表示位置总体精度的kappa系数为0.898，其内分县位置精度评价结果中，除鹿泉市、元氏县、灵寿县、栾城县和石家庄市外，其他8个显示的kappa系数均在0.90以上。可见，在各县区域内也取得了较好的统计数据空间化效果。其中，鹿泉市、元氏县、灵寿县kappa系数稍低的原因是上述县市为山前平原地区，境内地形稍有起伏，地块相对破碎，因此，中低分辨率MODIS影像提取冬小麦分布结果位置精度稍差。而石家庄市和栾城县由于处于城市及城市边缘地带，土地利用变化剧烈且类型多样，种植结构相对复杂，因此，MODIS提取冬小麦空间分布信息位置精度也略低。总体看，利用本发明作物种植面积统计数据空间分布的优化方法获取的冬小麦统计数据空间化结果与高分影像获得的作物分布数据间吻合度较好，具有较高的位置精度。

表2.冬小麦种植面积统计数据空间化结果位置精度评价

综上，本发明提出的在高时间分辨率遥感植被指数构成的时间序列数据支持下展开作物面积统计数据空间分布优化方法是一种通用体系框架(如图1所示)。在本文中仅以黄淮海地区石家庄平原区县冬小麦统计数据空间化为例利用SCE-UA优化算法对本发明方法进行阐述和验证。但是，从适合的作物来看，本发明方法除适合简单的一种作物的统计数据空间优化之外，还适合复杂种植结构下提取同季多种作物的空间分布，如玉米、水稻、大豆、棉花等；从应用的遥感数据来看，除MODIS植被指数时间序列数据外，还可以应用其他遥感传感器的高时间分辨率植被指数(如NOAA-AVHRR、SPOT-VGT、风云系列、环境小卫星等)；从应用的区域来看，本发明方法除适合平原区外，在适当提高使用遥感数据的空间分辨率的情况下，同样适合于其他复杂地形或破碎地块地区的作物种植面积统计数据空间化工作；从应用的优化算法来看，本发明方法也同样适合SCE-UA优化算法以外的其他全局优化算法(如模拟退火、遗传算法、四维变分和Kalman滤波算法等)；本发明方法适合所有依靠阈值算法控制作物分类结果的作物面积提取方法(如决策树法等)。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种作物种植面积统计数据空间分布的优化方法，其特征在于，包括以下步骤：

A、收集作物遥感影像及行政单元内作物实际种植面积农业统计数据，并预处理作物遥感影像；

B、在所述作物遥感影像上选择目标作物并建立一定区域内目标作物指数时间序列特征标准曲线；

C、利用所述作物植被指数时间序列特征标准曲线确定作物物候信息，并根据该作物物候信息确定作物分类算法模型；

D、利用全局优化算法并根据所述区域内作物实际种植面积农业统计数据，计算作物分类算法模型的特征参数的最优值；

E、将所述特征参数的最优值代入所述作物分类算法模型求出作物种植面积空间分布信息。

2.根据权利要求1所述的方法，其特征在于，步骤A中所述遥感影像的收集包括：

收集高时间分辨率最大值合成的一定时期作物指数时间序列数据。

3.根据权利要求1所述的方法，其特征在于，步骤A中所述作物遥感影像的预处理包括：

对所述作物遥感影像进行影像拼接、投影转换和裁切，并将所有所述作物遥感影像按照时间先后顺序排序后合成时间序列影像；

采用S-G滤波平滑法对时间序列影像进行去云平滑处理。

4.根据权利要求3所述的方法，其特征在于，所述步骤B包括：

作物生长期内地面样方调查；

利用差分GPS对地面调查样方中心点位置进行定位；

根据地面调查样方中心点位置结果，在时间序列的合成影像上确定作物种植地块中心位置所对应的像元并获取地面调查样方的作物指数时间序列特征曲线；

采用植被指数的上包络线来拟合每个作物地面调查样方的作物指数时间序列特征曲线；

对上述作物地面调查样方所对应的作物指数特征曲线进行平均，获得一定区域内的作物指数时间序列特征标准曲线。

5.根据权利要求3所述的方法，其特征在于，步骤C所述作物物候信息的确定为：

根据所述时间序列影像中的几个关键时期的作物指数值的变化趋势，提取出作物物候信息。

6.根据权利要求3所述的方法，其特征在于，步骤C所述作物分类算法模型的确定为：

根据所述作物物候信息反映出的作物指数时间序列规律，利用指数阈值法确定作物分类算法模型。

7.根据权利要求1、2、4、5或6所述的方法，其特征在于，所述作物指数时间序列为作物植被指数时间序列。

8.根据权利要求1所述的方法，其特征在于，所述步骤D包括以下子步骤：

D1、基于时间序列影像和所述作物分类算法模型进行作物分类，模拟出作物种植面积；

D2、由上述作物分类算法模型统计出所述区域内作物的模拟种植面积；

D3、建立由该区域内作物的模拟种植面积与该区域内作物实际种植面积农业统计数据之间的代价函数；

D4、当临近特定数量个代价函数值之差的绝对值小于特定值时，确定与最小代价函数值相对应的作物分类算法模型特征参数的初始值为特征参数的最优值。

9.根据权利要求8所述的方法，其特征在于，步骤D4所述特定数量为5，所述特定值为0.001。

10.根据权利要求9所述的方法，其特征在于，所述代价函数如下：

$y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{sim}}-X_{\mathrm{obs}})}^2$

其中，n是所述区域内作物实际种植面积统计数据的个数，X_sim是所述作物分类算法模型模拟出的所述区域内作物种植面积，X_obs是所述区域内作物实际种植面积的农业统计数据。

Images