专利名称：一种作物种植面积统计数据空间分布的优化方法
技术领域：
本发明涉及一种作物种植面积统计数据空间分布的优化方法。
背景技术：
作物种植面积统计数据是农业统计数据的核心指标之一。准确的作物种植面积统计数据不仅关系到国家农业政策的正确制定，而且对提高农业生产管理水平、增强农业资源监测能力、优化作物种植结构、准确预测作物产量、确保国家粮食安全具有重要意义。近年来，随着世界各国对全球气候变化、人口、资源、环境、生态和粮食安全等一系列问题定量研究的逐步深入，准确的农业统计数据在上述问题研究中发挥了重要作用。其中，作物种植面积统计数据成为不可或缺的基础数据之一。然而，包括作物种植面积统计数据在内的传统农业统计大多以各级行政单元为统计单位，统计数据除仅仅反映行政单元内总量水平或平均水平外，数据存在统计单元不统一、定位不准确和空间分辨率低等缺点，尚未真正揭示统计数据在地理空间上的真实分布和空间变异情况，更不利于统计指标与其他自然、生态要素空间数据(如土地利用/ 土地覆盖数据、生态环境背景数据)实现交叉使用和叠置分析等功能。传统的区域农业统计指标空间化的方法是直接给区域行政单元赋以统计数值，然后生成图斑表现形式的专题地图， 而图斑的边界就是行政单元的边界。通过这种方法获得的统计指标空间分布信息具有同一区域内指标均勻分布，不同区域间指标存在突变分布的特点，这与农业统计数据空间分布的真实状态不一致。此外，这种方法形成的空间分布结果难以满足跨行政边界的农业数据统计需要以及与自然、生态要素空间数据的交叉使用和叠置分析。可见，传统的农业统计数据空间化方法已经很难满足社会经济和科学技术发展的应用要求。近年来，随着遥感、地理信息系统、全球定位系统和计算机技术的发展，以行政单元为空间表达方式的各类统计数据空间化已经成为空间信息技术、全球变化等社会交叉学科领域的研究热点，即将基于社会行政单元(如省、市、县、乡镇等)的统计信息转化为基于一定尺寸网格单元的空间信息，实现在空间网格统计单元上与土地利用、生态环境背景数据等自然要素数据进行综合分析，真正实现自然、生态和人文要素的整合。目前，国内外社会经济统计数据空间化研究主要集中在人口统计数据空间化、国内生产总值(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 = Y^iXsim-Xobsf
/=1其中，η是所述区域内作物实际种植面积统计数据的个数，Xsiffl是所述作物分类算法模型模拟出的所述区域内作物种植面积，Xobs是所述区域内作物实际种植面积的农业统计数据。
图1为本发明提供的作物种植面积统计数据空间分布优化方法框架图；图2为本发明提供的作物种植面积统计数据空间分布优化方法的研究区位置和调查样区点位示意图3为本发明提供的作物种植面积统计数据空间分布优化方法的具体实施流程图；图4为研究区域内冬小麦从播种到成熟期间M0DIS-NDVI时间序列数据的平滑效果图；图5为研究区域内冬小麦NDVI时间序列特征曲线图；图6a为通过本发明方法获得的整个研究区域内的冬小麦种植面积统计数据空间化结果示意图；图6b为通过高分辨率遥感影像获取的冬小麦分布数据示意图。
具体实施例方式本发明是一种在高时间分辨率遥感植被指数构成的时序数据支持下展开作物面积统计数据空间分布优化方法。由于植被生长以年为周期，在这个生长周期内不同植被类型有着各自生繁衰枯的物候节律，表现出不同的生长规律，而且规律性极强。这种规律性可以作为植被分类的出发点。作物生长周期内具有其独特的生长规律(如播种期、出苗期、起身期、拔节、抽穗、乳熟、收获期等)，但不同作物的叶子内部构造存在一定的差异。因此，不同作物的物候期会存在一定差异，导致不同作物的波谱形态时间曲线会有所不同，但同一种作物在同一地区具有相对稳定的作物生长发育规律和NDVI-时间变化曲线(NDVI， Normalized Difference Vegetation hdex，归一化植被指数，是一个能够通过遥感获取且能够直接、有效地反映作物绿度长势、盖度、生物量和作物单产的最佳参考量之一)。因此， 本发明采用多时相植被指数遥感信息，根据绿色植物独特的波谱反射特征，将植被(作物) 与其他地物区分开，以便进行作物空间分布信息提取。所述优化方法的主要框架如图1所示，包括首先，利用平滑过滤技术提高时序植被指数自身质量，以便消除云、气溶胶及传感器不稳定等因素产生的噪音对数据质量的影响；然后，在地面调查数据支持下，通过遥感影像获取作物时序植被指数构成的标准时序特征曲线；再利用该时序特征曲线反映的作物物候信息，通过NDVI阈值法确定提取作物分类算法模型，并指定阈值初始值和变化范围；在优化算法支持下，以作物区域统计面积数据作为优化比较对象，当模拟的区域作物面积与区域实际作物面积农业统计数据间差异达到最小时，得到作物分类模型最优参数；最后，在最优参数驱动下，利用作物分类算法模型提取作物空间分布信息，从而实现作物面积统计数据的空间分布优化。下面，以位于中国北方粮食生产基地黄淮海平原内河北省石家庄地区13个地处平原的县(市)作为研究区域(如图2所示，37.60° -38.77° N，113. 81° -115.35° E， 覆盖面积8230km2)，详细介绍本发明提供的作物种植面积统计数据空间分布优化方法在该研究区域内的具体实施过程。由于山区县市地形复杂，作物种植面积破碎，如采用中低分辨率遥感影像，则很难满足山区或半山区作物种植面积提取要求。因此，在本实施例中仅考虑了研究区域内的部分平原县市。该区域属温带半湿润季风气候，彡0°C年积温4200°C 5500°C，年累计辐射量5. 0*106 5. 2*106KJ · πΓ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米 M0DIS-NDVI时间序列数据(M0DIS数据由terra和aqua卫星的一个重要的传感器实时观测得到的，由terra和aqua卫星通过χ波段向全世界免费实时广播发送)。MODIS数据具有波段多(36个)、中低空间分辨率Q50 1000米)和时间分辨率高(1 2天)等特
点，非常利于与作物的生长规律相结合，因此被广泛应用于农作物动态监测业务之中，如作
R -R
物面积提取、作物估产等。而NDVI的计算公式为-W =’其中I^n为近红外波段的
ηr
反射率，艮为红光波段的反射率。所述遥感影像的收集过程为从EOS Data Gateway中下载2008年第273天(9月四日)到2009年第177天(6月沈日)期间每16天合成的 M0DIS-NDVI数据。所述遥感影像的预处理包括对所述遥感影像进行影像拼接、投影转换和裁切，并将所有遥感影像按照时间先后顺序排序后合成时间序列影像；接着，采用S-G滤波平滑法(Mvitzky-Golay)对时间序列影像进行去云平滑处理。在上述预处理过程中，可采用公知的MRT(MODIS Reprojection Tool，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示出了上述研究区域内冬小麦从播种到成熟期间M0DIS-NDVI时间序列数据的平滑效果，其中，横坐标格式为年份+每16天数据的起始天。第二步骤在地面调查数据(作物野外地面样方调查信息)的支持下，在所述作物遥感影像上选择目标作物并建立目标作物植被指数时间序列特征曲线。为了获取准确的作物植被指数时序特征曲线，在2009年3月开展了作物返青-拔节期(生长期)地面样方调查。其中，在石家庄地区13个县选择了 45个冬小麦地面样方进行观测。考虑到冬小麦长势、种植结构的代表性以及调查点分布的均勻性，样区间隔不小于5km，每个样区面积不小于500m*500m，利用差分GPS对调查样区中心点位置进行精确定位。然后，根据野外地面调查样方中心点位置结果，在M0DIS-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。针对每个像元找出作物光谱曲线中两个波峰所对应的时相号rl、r2，并利用NDVI阈值法(公知技术，在此不再赘述)确定下列作物分类算法模型a、由冬小麦出苗期的NDVI时间序列规律可知，满足判别函数NDVI (rl)-NDVI⑵
>Tl，NDVI (rl) > 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 (rl)
>T4的像元为冬小麦；其中，T1、T2、T3和Τ4为作物分类算法模型的待优化的特征参数，其初始值为0 1内的任意值，其变化范围可以从MODIS时间序列的影像上获取，即变化范围为0 1。第四步骤在全局优化算法支持下，以作物种植面积统计数据作为外部优化数据比较对象，当模拟的区域作物面积与区域实际作物面积农业统计数据间差异达到最小时， 得到作物分类模型的最优特征参数。在本步骤中，利用优化效率较高、运算速度快且对参数初始值设置不敏感的 SCE-UA 优化算法(Shuffled Complex Evolution-University of Arizona，复合形混合演化算法)对作物分类算法模型的待优化的特征参数T1、T2、T3和T4进行优化。具体优化过程如下
i、基于M0DIS-NDVI时间序列数据和上述作物分类算法模型进行冬小麦分类(模拟出冬小麦种植面积)；ii、由上述作物分类算法模型统计出研究区域内各县冬小麦的模拟种植面积；iii、建立研究区域内各区县冬小麦的模拟种植面积与研究区域内各区县冬小麦的实际种植面积统计数据之间的代价函数如下
权利要求
1.一种作物种植面积统计数据空间分布的优化方法，其特征在于，包括以下步骤A、收集作物遥感影像及行政单元内作物实际种植面积农业统计数据，并预处理作物遥感影像；B、在所述作物遥感影像上选择目标作物并建立一定区域内目标作物指数时间序列特征标准曲线；C、利用所述作物植被指数时间序列特征标准曲线确定作物物候信息，并根据该作物物候信息确定作物分类算法模型；D、利用全局优化算法并根据所述区域内作物实际种植面积农业统计数据，计算作物分类算法模型的特征参数的最优值；E、将所述特征参数的最优值代入所述作物分类算法模型求出作物种植面积空间分布 fn息ο
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 = YjixSim ~xobs)2 /=1其中，η是所述区域内作物实际种植面积统计数据的个数，Xsiffl是所述作物分类算法模型模拟出的所述区域内作物种植面积，Xobs是所述区域内作物实际种植面积的农业统计数据。
全文摘要
本发明提供了一种作物种植面积统计数据空间分布的优化方法，包括以下步骤A、收集作物遥感影像及行政单元内作物种植面积统计数据，并预处理作物遥感影像；B、在所述作物遥感影像上选择目标作物并建立一定区域内目标作物指数时间序列特征标准曲线；C、利用所述作物指数时间序列特征标准曲线确定作物物候信息，并根据该作物物候信息确定作物分类算法模型；D、利用全局优化算法并根据所述区域内作物实际种植面积农业统计数据，计算作物分类算法模型的特征参数的最优值；E、将所述特征参数的最优值代入所述作物分类算法模型求出作物种植面积空间分布信息。从而实现作物统计数据空间化信息精度的进一步提高。
文档编号G06Q10/04GK102592181SQ20111045624
公开日2012年7月18日 申请日期2011年12月30日 优先权日2011年12月30日
发明者任建强, 刘佳, 周清波, 唐华俊, 王晋臣, 王道龙, 陈仲新 申请人:中国农业科学院农业资源与农业区划研究所