基于MODIS数据的区域大尺度作物种植结构提取方法与流程

本发明涉及基于modis数据的区域大尺度作物种植结构提取方法，属于农业种植结构技术领域。

背景技术：

农作物的种植结构对高效作物规范管理具有重要的意义，农作物的遥感影像获取是种植结构调整的关键环节，在大尺度下准确提取作物种植结构关乎整个农业发展。

遥感技术的不断发展带动农业遥感的不断进步，目前农业遥感应用于作物估产、作物面积提取、灾害监测、土地覆盖分类、作物时空格局演变研究等诸多领域，是农业管理的一个重要基础，也是制定农业发展策略乃至整个国民经济发展方向的重要依据之一。但农用地遥感研究一般集中在两个尺度上：一是地块尺度的作物分类；二是土地利用中的耕地提取，部分研究划分了水田和旱地。地块尺度的作物分类对遥感数据的时间、空间分辨率以及地面样本数据要求高，且成本高、耗时耗力；耕地提取结果中，耕地作为最小信息单元，缺乏内部的作物种植结构信息，尽管获取相对容易，但难以用于农业生产管理、生态环境时空模型的定量模拟。

作物种植结构的研究一方面影响着作物长势遥感监测的准确性与发展方向，而且在分作物、分区的模型建立中也起到重要的作用。为了能够更好地、准确地获取农作物的种植结构、面积、产量等信息，并能及时、精确地取得农业灾情等关键信息，有关作物种植结构的研究也越来越受到了较多学者的关注。因此，如何获取较大尺度范围内农作物分布情况已经成为农业遥感、作物生长状况监测、农情预测的核心问题。传统作物类型及种植结构的提取方法通常是根据波谱特征的统计差异与训练区的样本相结合进行分类的。但对于混作区的定义却比较模糊，仅仅依据光谱特征进行分类使得农作物生长的季节性和物候特征的稳定性规律难以发挥优势。而且不同地区同一类型的作物生长节律不同，导致植被覆盖指数ndvi的时间序列曲线形状也不相同，会带来较大的分类误差。

技术实现要素：

本发明目的是为了解决现有农业遥感监测中农作物种植结构的提取方法分类误差大的问题，提供了一种基于modis数据的区域大尺度作物种植结构提取方法。

本发明所述基于modis数据的区域大尺度作物种植结构提取方法，它包括以下步骤：

步骤一：采集大尺度监测区域的全年影像数据，并进行预处理，将同一采样时刻的影像数据进行影像拼接，获得多张大尺度监测区域的完整影像；

步骤二：对所述的完整影像根据相应的耕地矢量图范围进行裁剪，建立提取年份的完整时序文档，并获取每个像元的ndvi时间序列曲线；

步骤三：对ndvi时间序列曲线趋势进行平滑重构，提取全年的11个物候数据，并对11个物候数据进行浮点化；

步骤四：对11个物候数据进行合成处理，获取综合信息数据，对综合信息数据进行标准化，再对标准化综合信息数据进行主成分分析，选取最优波段，提取每个物候数据的关键信息；

步骤五：根据物候数据的关键信息对所有完整影像进行分区获得多个物候分区；

步骤六：对物候分区进行多尺度分割获取种植结构单元；

步骤七：对种植结构单元采用最邻近分类方法进行基于样本的分类提取，获得大尺度监测区域的作物种植结构。

本发明的优点：本发明利用timesat进行物候数据提取，在分区基础上采用面向对象分割的种植结构提取方法获取种植结构信息。它采用modis遥感影像，空间分辨率为250m，时间分辨率为16天，在省级尺度，250m空间分辨率十分适合种植结构提取的研究。

本发明提出了作物种植结构单元(croppingpatternunit,crpu)这一概念，构建出一种快速、低成本、准确的区域尺度crpu提取方法。它应用遥感手段进行物候数据提取，对区域大尺度范围进行农业分区，避免由于同地区间生长季长短不一以及气候、地貌、植被、土壤等各种条件对同一区域内各种作物的配制结构制约等的影响，区域大尺度分区的方法更简单、直观、准确、快速；分区的界线更为合理，更能反映区域间的综合特征和微妙差异。同时该发明做到精准、及时、方便，很真实的反映作物的空间分布形式，为农业决策带来准确的分布结果。本发明的应用也是保证农业粮食安全、进行种植结构优化的依据，也对农业产业结构调整乃至制定进一步战略性目标具有一定理论意义。由此可见，开展基于遥感技术的作物种植结构监测具有深远的意义。

附图说明

图1是步骤三中采用savitsky-golay算法平滑前后各作物ndvi时间序列图；图中点数据为平滑前数据，线数据为平滑后数据；

图2timesat中的物候参数起始时间图；

图中：a表示作物生长期起始点；

e表示作物生长期终止点；

b表示上升幅度为整体增幅的20％时的某一时间；

d表示下降幅度为整体降幅的20％时的某一时间；

c表示最大值；

g表示生长过程中的振幅；

f表示生长季的长度；

h和h+i表示两种积分。

图3是耕地物候分区图；图中轮廓线为分区边界；

图4是黑龙江省耕地物候分区最优分割尺度评价图；

图5是采用本发明方法进行多尺度分割，获取种植结构单元后进行评价的曲线图；

图6是本发明提取种植结构分类体系的特征空间构建图；

图7是以黑龙江省为例的种植结构提取结果图。

具体实施方式

具体实施方式一：下面结合图1至图7说明本实施方式，本实施方式所述基于modis数据的区域大尺度作物种植结构提取方法，它包括以下步骤：

步骤一：采集大尺度监测区域的全年影像数据，并进行预处理，将同一采样时刻的影像数据进行影像拼接，获得多张大尺度监测区域的完整影像；

步骤二：对所述的完整影像根据相应的耕地矢量图范围进行裁剪，建立提取年份的完整时序文档，并获取每个像元的ndvi时间序列曲线；

步骤三：对ndvi时间序列曲线趋势进行平滑重构，提取全年的11个物候数据，并对11个物候数据进行浮点化；

步骤四：对11个物候数据进行合成处理，获取综合信息数据，对综合信息数据进行标准化，再对标准化综合信息数据进行主成分分析，选取最优波段，提取每个物候数据的关键信息；

步骤五：根据物候数据的关键信息对所有完整影像进行分区获得多个物候分区；

步骤六：对物候分区进行多尺度分割获取种植结构单元；

步骤七：对种植结构单元采用最邻近分类方法进行基于样本的分类提取，获得大尺度监测区域的作物种植结构。

步骤一中的预处理可以包括以下至少一种方式：格式与投影转换及利用双线性内插的采样方法提取ndvi。

步骤五中获得多个物候分区的具体方法可以为：

对完整影像利用多尺度间隔执行分割，在多次分区的基础上，利用各区域间同质性及异质性计算评价指数，对每次的分区结果进行有效评价，获取最优分区尺度，将按照最优分区尺度进行分区获得的结果作为物候分区结果。

步骤一中，采集的全年影像数据，可以是下载的modis16天合成的陆地产品(modis/terravegetationindices16-day/mod13q1/16d)，包括生长季在内的全年影像，对遥感影像利用modis转换工具——modisreprojectiontool(mrt)进行格式与投影转换、采用双线性内插的采样方法提取植被指数(normaldifferentialvegetationindex/ndvi)，利用erdas软件拼接工具(mosaictool)进行影像拼接等预处理。

步骤二中，将步骤一中处理后的影像根据区域大尺度耕地矢量图范围进行裁剪，在envi软件中存为标准格式，建立一个完整提取年份的时序的txt文档，从而获取每个像元的ndvi时间序列曲线。

步骤三中，根据时间序列曲线趋势，选择savitzky-golay(s-g)滤波法进行平滑重构，如图1所示，选取动态阈值法，将作物生长期开始时间定义为左侧部分ndvi曲线急剧上升点，即上升幅度为整体增幅的20％时对应的时间。同时，定义作物生长期结束时间为急剧下降点，即下降幅度为整体降幅的20％时对应的时间，见附图2，利用timesat软件进行全年11个物候数据的提取，进而将利用envi软件中编辑功能(editheader)进行浮点化，将11个数据分别设置为浮点格式(floatingpoint)后存为标签图像文件格式——tiff(tagimagefileformat)文件。

所述11个物候数据为作物生长起始期、作物生长结束期、振幅、ndvi的平均值、生长期长度、ndvi的积分、ndvi最大值、左侧上升曲线间斜率、右侧下降曲线间斜率、整个时期的中间点和整个时期ndvi的积分。

步骤四中，利用arcgis中的波段合成(compositebands)模块，将步骤三获取的11个物候数据进行合成处理获取综合信息数据，在envi中利用波段运算(bandmath)模块，利用标准分数(standardscore)对合成数据进行标准化，从而降低由各参数的数量级不同而造成的异常影响。利用envi5.3对标准化数据进行主成分分析选取最优波段，提取出每个物候信息中的关键信息，降低波段冗余现象。

z＝(x-μ)/σ，

式中z为某物候的具体分数，x为某一具体物候参数，μ为平均数，σ为标准差。

步骤五中，为了避免由于尺度过大造成的区域间物候差异，利用econgnitiondeveloper8.7软件，将步骤四处理后的物候数据再进行物候分区，获得多个物候区，如附图3所示，利用多尺度间隔执行分割，在多次分区的基础上，利用各区域间同质性及异质性计算评价指数，进而对分区结果进行有效评价，获取最优分区尺度，如附图4所示，最优点为曲线最高点；具体方法如下：

1)将合成后待分区的物候数据加入econgnitiondeveloper8.7中；

2)利用多尺度分割功能(multiresolutionsegmentation)进行分区；

3)定义不同尺度(scale)下的形状及紧致度，进行差异比较。

4)输出同质性(standarddeviation)和异质性(meandiffto)进行评价计算。

步骤六中，在步骤五的基础上对获取的各物候区进行多尺度分割，获取种植结构单元，充分利用对象(单元)特征，基于某种尺度分割出多个对象，将对象标准差(standarddeviation)作为同质性σp的标准，与邻域的平均差分的绝对值(meandiff.to)作为异质性δcw的标准，在二者之间找到一个最优点，即分割评价指数sei，保证同类作物种植结构的聚类，不同类间差异显著。分割后，对种植结构单元的提取结果进行评价，如附图5所示，本发明的最优尺度为急剧下降后趋于平缓的点。

式中，σp为物候指标l的标准差；cpi为像元i的灰均值，cp为灰度均值；n为crpu内像元的个数；n为与当前crpu邻接crpu的个数；s为周长，si为与第i相邻crpu公共边长；δcw为与相邻对象的均差分绝对值，cw为当前物候指标的灰度均值，cwi为第i个相邻crpu的灰度均值；sei为分割评价指数；m为波段个数；ql为不同物候波段l的权重。

通过单元面积变量的引入，减少由于个别区域面积较小使得面积较大的单元对权重造成的影响而带来的评价结果的偏离。最终利用平均分割评价指数进行分割评价。平均分割评价指数(averagesegmentationevaluationindex,asei)公式为：

式中，r表示所有种植结构单元——crpu(croppingpatternunit)的总面积；ri表示i对象的crpu的面积；z表示crpu数量；seii表示i对象的crpu的分割评价指数。

步骤七中，在步骤六获取种植结构对象后，需将种植结构各对象赋予的种植信息属性进行分类提取，本发明建立提取种植结构的分类体系，是利用各物候主成分的均值、标准差、亮度等光谱特征；形状指数、长宽比等形状特征构建特征空间，形成分类规则，如附图6所示，选取最邻近分类方法进行分类提取。具体方式如下：

1)在econgnitiondeveloper8.7软件的分类结构(classhierarchy)栏中建立详细的分类体系；

2)选择分类样本，在菜单中选择分类功能(classification)，进行样本选择(selectsamples)的执行；

3)然后是构建特征空间，如上所述，利用特征空间模块功能(standardnnfeaturespace)选择一些特征信息进行构建，为避免分类特征冗余，该发明最大特点就是可以进行特征优选，在特征空间优化(featurespaceoptimazation)里面进行特征维数的计算，当维数为5时，种植结构对象与对象之间的可分性最大，最后是执行分类。

本发明以黑龙江省为研究区，将通过物候数据分割获取的种植结构单元视为单独对象，采用最邻近分类方法对各分割后的对象进行基于样本的分类提取，最终获取区域大尺度全省的种植结构，如附图7所示，为种植结构调整提供科学依据。

本发明利用全年250m-modis时间序列遥感影像，采用非对称高斯函数拟合方法，对植被指数时间序列曲线进行平滑重构处理，进而根据曲线变化特征提取了区域大尺度作物物候参数；在提取物候信息的基础上进行耕地物候分区；对各物候区进行多尺度分割，最终提取区域大尺度作物种植结构分布信息。

本发明中所述的区域大尺度一般指市、省级尺度或更大尺度。