一种基于冬小麦面积指数的冬小麦种植面积计算方法

摘要

本发明涉及一种基于冬小麦面积指数的冬小麦种植面积计算方法，包括：对多时相遥感影像进行预处理，获取冬小麦生育期内的NDVI数据；将多时相遥感影像划分为规则网格，在规则网格上设置权值样本中心点，采用目视判读的方式确定网格中的样本类型为冬小麦和非冬小麦；获取多时相网格数据中权值样本中心点对应的NDVI值，分别计算出冬小麦和非冬小麦的NDVI平均值，比较两者大小并设置权值；根据权值设置结果对NDVI平均值进行加权操作，以获得冬小麦面积指数；利用自适应方式提取冬小麦面积指数最优阈值，据此求出冬小麦的种植面积。本发明提供的方法可以实现多时相高分辨率影像的区域尺度冬小麦面积高精度计算，并且方法的自动化程度高，识别结果一致性较强。

说明书

技术领域

本发明涉及农业遥感技术领域，特别涉及一种基于面积指数的多时相全国或省级尺度遥感影像冬小麦面积提取方法。

背景技术

大宗农作物的播种面积空间分布数据的准确获取是国家、省级政府制定粮食政策，保证粮食安全的重要依据，作为国家大宗农作物的冬小麦，占据国家粮食总量的较大比例。由于播种面积是作物估产的关键步骤，因此，开发目前区域尺度的冬小麦面积提取方法，对农情监测具有重要的意义。随着遥感技术的迅速发展，其为快速、准确获取作物种植空间分布提供了有效的手段，利用遥感技术进行冬小麦面积提取识别主要方法有单一时相法、多时相法、空间抽样法等。单一时相法通过选取关键物候期的单景遥感影像，利用监督分类、非监督分类、面向对象分类等多种技术方法，基于地物的光谱、文理特征等信息，对作物进行识别及面积提取；多时相法是根据待识别目标作物的物候特征，选取合适的多时相影像，通过构建植被指数或者光谱特征信息集，采用与单时相数据类似的方法进行作物分类及面积提取；空间抽样法是遥感与传统统计调查方法结合的农作物面积估算方法，针对区域尺度上全覆盖遥感影像获取困难、分类技术难以标准化而提出的，可以在统计抽样框架下，通过遥感技术获取抽样区域的作物面积推算区域农作物面积数量。以上这几种方法多基于中低分辨率遥感影像，且对专家知识较为依赖，缺少基于多时相中高分辨率遥感影像的区域尺度大宗农作物面积自动识别方法的研究。

归一化差分植被指数(Normalized Difference Vegetation Index，NDVI)是目前用于检测植被生长状态、植被覆盖度和消除部分辐射误差的主要指数，对于反映高分辨率遥感数据中的作物也具有一定的优势，然而根据各时相NDVI阈值进行作物类型的判别需要依赖足够的专家知识，人为主观性较强。

发明内容

本发明目的在于，针对如何利用多时相、中高分辨率遥感影像提取区域尺度下的冬小麦面积同时使得冬小麦面积遥感监测具备较高的自动化运行能力的问题提供解决方法。

为解决上述技术问题，本发明提供了一种基于冬小麦面积指数(Winter WheatArea Index，WWAI)的冬小麦种植面积提取方法，该方法通过对遥感影像中样本点的类型进行目视判读，统计样本点各类型的NDVI值来设置作物提取阈值，可实现更加高效快速的目标类型提取。因此，在归一化植被指数(NDVI)时序影像的基础上，设置类型识别的样本点，通过样本点目视判别结果设置作物提取阈值，以实现高分辨率遥感影像中区域尺度作物更加自动化的快速提取。

本发明提供的方法包括以下步骤：

S1：对研究区多时相遥感影像进行预处理，获取冬小麦生育期内的NDVI时序数据；

S2：按等间距原则将研究区多时相遥感影像划分为规则网格；

S3：在规则网格上设置权值样本和阈值样本中心点，采用目视判读的方式确定网格中的样本类型，并将类型分为冬小麦和其它两种地物类型；

S4：获取多时相网格数据中权值样本中心点对应的NDVI值，分别计算出冬小麦和其它地类两种类型的NDVI平均值，并比较两者大小进行权值的设置；

S5：根据权值设置结果对NDVI平均值进行加权操作，以构建冬小麦面积指数(WWAI)值；

S6：采用自适应方式提取冬小麦面积指数最优阈值，利用所述冬小麦面积指数最优阈值求出所述区域内冬小麦的种植面积。

具体而言，所述方法包括以下步骤：

S1：在待计算的冬小麦种植区域内，获取涵盖多个冬小麦生育期的多时相遥感影像；其中，每个生育期对应一幅影像，每幅影像均来源于同一卫星的同一传感器；

对所述影像进行预处理，获取与影像中每个像素点一一对应的NDVI；

S2：将每个生育期对应的影像分别按等间距原则进行划分，得到多个规则排列的形状、大小相同的基准网格；

S3：在每个所述基准网格内相同的位置上选择一个像素点，作为权值样本中心点，采用目视判读方式判定各个权值样本中心点为冬小麦样本或为非冬小麦样本；

S4：统计每个生育期对应的影像中权值样本中心点对应的冬小麦样本总数和非冬小麦样本总数，依据每个权值样本中心点对应的NDVI，分别求出该生育期内冬小麦样本的NDVI平均值和非冬小麦样本的NDVI平均值；

若所述生育期内冬小麦样本的NDVI平均值大于所述非冬小麦样本的NDVI平均值，则设定该生育期内冬小麦权值为1，否则设定为-1；

S5：将每个生育期内冬小麦样本的NDVI平均值乘以该生育期内冬小麦权值，得到每个生育期内的冬小麦加权值；求出全部所述多个生育期内的冬小麦加权值的平均值，即得冬小麦面积指数；

S6：采用自适应方式提取某一生育期影像中冬小麦面积指数的最优阈值，逐一比较该生育期影像中每个像素点对应的归一化差分植被指数与所述生育期冬小麦面积指数最优阈值的大小，若像素点对应的归一化差分植被指数不小于所述最优阈值，则判断该像素点为冬小麦，否则为非冬小麦；统计所述影像中冬小麦像素点的总个数，用所述冬小麦像素点的总个数×每个像素代表的面积，即得所述生育期内冬小麦的种植面积。

上述方法中：

步骤S1所述多个生育期包括从播种至成熟的多个时期；优选为冬小麦播种、出苗、分蘖、越冬、返青、拔节、抽穗、乳熟和成熟的全部生育期。

其中，所述生育期内每月至少一幅影像需为云量小于10％的晴空影像，拍摄于各个生育期内的任意一个时间点，优选均拍摄于各个生育期的时间中点。

为了提高步骤S1所采集影像的效果，从而确保面积提取结果的准确性，影像的获取根据冬小麦的生长时间确定，且保证每月1景无云影像。拍摄步骤S1所述影像的遥感卫星相机所覆盖的光谱波段包括红光波段和近红外光波段，优选为波长0.63μm～0.69μm的红光波段和波长为0.77μm～0.89μm的近红外光波段。拍摄步骤S1所述影像的为中高空间分辨率遥感卫星影像，分辨率优选为1～50m，进一步优选为1～20m；相应的，所述每个像素点对应的种植区域面积为1m²～2500m²，优选为1m²～500m²。

其中，所述获取NDVI的计算方法为：NDVI＝(NIR-R)/(NIR+R)；式中，NIR为近红外波段处的反射率，R为红光波段处的反射率。所述NDVI的绝对值≤1，即-1<＝NDVI<＝1，若NDVI为负值，则表示地面覆盖为云、水、雪等，对可见光高反射；若NDVI为0，则表示有岩石或裸土等，NIR和R近似相等；若NDVI为正值，表示有植被覆盖，且随覆盖度增大而增大。

步骤S1所述预处理包括辐射定标、大气校正、几何校正中的一种或多种。为了提高预处理效果，所述辐射定标方法优选为：L(λ)＝Gain·DN+Bias；式中，L(λ)为传感器入瞳处辐射亮度值(可采用如下公式进行计算：W/(m²·sr·μm))，Gain为增益系数，DN为观测灰度值，Bias为偏置系数。所述大气校正优选使用6S模型进行。所述几何校正优选使用区域网平差模型进行。

所述步骤S2中，影像网格划分是将按照一定等间距原则进行的。为了确保计算的准确性，每个基准网格对应的种植面积优选为1km²～50km²，更优选为1km²～20km²。所述基准网格优选为矩阵，即影像中基准网格的行、列数相等。

作为一种优选方案，所述步骤S2中，将每个生育期对应的影像按等间距原则进行3划分的结果作为基准网格，再将每个基准网格进行四等分，形成二级网格；所述步骤S3中，将每个基准网格中位于左上方的二级网格的中心点作为权值样本中心点。

所述步骤S4的具体定义与构建公式如下：

$>P^i=\left(\begin{array}{cc}1& ({\mathrm{NDVI}}_w^i>{\mathrm{NDVI}}_o^i)\\ -1& ({\mathrm{NDVI}}_w^i<{\mathrm{NDVI}}_o^i)\end{array}\right);$>

$>WWAI=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}({\mathrm{NDVI}}^i\times P^i);$>

式中，P为权值，i表示第i期影像，下标w表示冬小麦类型，下标o表示非冬小麦类型。

为了提高本所述方法提取冬小麦种植面积的准确性，本发明优选将设置阈值样本。

具体而言，所述步骤S3还包括：在每个所述基准网格内的另一相同位置设置一个阈值样本中心点，采用所述目视判读方式判定各个阈值样本中心点为冬小麦样本或为非冬小麦样本。

作为一种优选方案，所述步骤S2中，将每个生育期对应的影像按等间距原则进行划分的结果作为基准网格，再将每个基准网格进行四等分，形成二级网格；将每个基准网格中位于右下方的二级网格的中心点作为阈值样本中心点。

步骤S6所述自适应方式具体为：将所述冬小麦面积指数值进行等分，获得多个面积指数阈值。为了提高准确性，优选将所述冬小麦面积指数自0％到100％等分为101个面积指数阈值，各相邻面积指数阈值之间的差均为1％；举例而言，当冬小麦面积指数为50时，由所述面积指数的0％至所述面积指数的100％，等分为101个面积指数阈值，具体为0、0.5、1、1.5、2、2.5……49、49.5、50，各相邻面积指数阈值之间的差均为所述面积指数的1％(即0.5)。

对于每一个面积指数阈值：将该面积指数阈值与步骤S6所述生育期对应影像中各个阈值样本中心点的归一化差分植被指数进行比较，若阈值样本中心点的归一化差分植被指数不小于所述面积指数阈值时，判断该阈值样本中心点为冬小麦，否则为非冬小麦；以阈值样本中心点的目视判读方式判定结果为标准检验上述判断结果，计算该面积指数阈值的判断正确率；统计所述多个面积指数阈值的判断正确率，以正确率最高面积指数阈值作为所述生育期冬小麦面积指数的最优阈值。

本发明提供的方法可以实现多时相高分辨率影像的区域尺度冬小麦面积高精度提取，并且方法的自动化程度高，识别结果一致性较强。

附图说明

图1为本发明提供的冬小麦面积指数法流程图；

图2为2014年5月17日安平县及周边地面网格划分及样本点选取示意图；

图3为2014年5月17日安平县及周边冬小麦面积指数影像；

图4为2014年5月17日安平县及周边冬小麦空间分布图。

具体实施方式

以下实施例用于说明本发明，但不用来限制本发明的范围。

实施例1

以提取安平县及周边冬小麦面积为例进一步阐述本发明的技术方案。本实施例的流程如图1所示，包括：

步骤S1，试验区多时相遥感数据获取及预处理，并计算出各时相的NDVI。

试验区主要位于安平县及周边地区，安平县隶属河北省衡水市，地理坐标为E115°19′至115°40′，N38°14′至38°21′，海拔高度26m。属半干旱半湿润大陆性季风气候区，春季少雨多风，夏季高温多雨，秋季气温适中，冬季寒冷少雪，四季分明，冷暖干湿差异较大，年日照时数2400～3100h，年均降水量300～800mm，1月平均气温在3℃以下，7月平均气温18～27℃。

安平县及周边试验区内的冬小麦一般在每年秋季10月10日前后陆续开始播种，下一年6月20日前完成收获，生育期近8个月。按照生长季内冬小麦发育特点，选取覆盖试验区的2013年10月26日、11月27日、12月29日、2014年1月3日、2月9日、3月29日、4月14日、5月17日、6月27日的高分一号卫星9景WFV影像，涵盖了冬小麦从种植到收获的整个生长发育期。高分一号卫星共有4台16m分辨率多光谱相机(WFV1～WFV4)，每台相机包含蓝(0.45～0.52μm)、绿(0.52～0.59μm)、红(0.63～0.69μm)和近红外(0.77～0.89μm)4个波段，4台相机组合宽度可达800km，重访周期2d。在大气校正和几何校正的基础上，计算各期影像的NDVI值，计算公式为：

$>NDVI=\frac{{\mathrm{Ref}}_4-{\mathrm{Ref}}_3}{{\mathrm{Ref}}_4+{\mathrm{Ref}}_3}$>

式中，Ref₄和Ref₃分别为WFV第4和第3波段的反射率。

步骤S2，按照等间距原则将覆盖试验区的影像划分成规则格网。考虑到试验区的大小和作物种植分布情况，根据等间距的规则获取10×10个网格，作为一个基准网格，试验区每个网格的大小为4.1km×4.1km；在每个基准网格的基础上，再将其划分为2×2个格网，分别获取左上、右下2个子格网的中心，称为左上中心、右下中心；在试验区范围内，总共生成100个左上网格和100个右下网格(见图2)。

步骤S3，在规则格网上设置权值样本和阈值样本，将样本类型分为冬小麦和非冬小麦两种类型，通过目视识别的方式判断样本的类型。以每个基准网格中的100个左上中心作为权值样本，根据目视判断出47个冬小麦点，53个非冬小麦点，冬小麦比例47％；以基准网格中的100个右下中心作为阈值样本，计有44个冬小麦，56个非冬小麦点，冬小麦比例44％。

步骤S4，根据权值样本中心点的NDVI值，分别计算出冬小麦和非冬小麦两种类型的NDVI平均值，并比较两者大小进行权值的设置。分别计算每期NDVI影像两种地物类型的权值样本平均值，具体到安平县及周边地区，就是计算9期影像中每期47个冬小麦类型以及53个非冬小麦权值样本的平均值，分别获得两类各9个平均值；比较每期两种类型NDVI平均值的大小，若当前冬小麦NDVI平均值大于非冬小麦NDVI平均值，则权值设为1，若当期冬小麦NDVI平均值小于非冬小麦NDVI平均值，则权值设为-1。

步骤S5，将设定的权值赋予到各期对应的冬小麦NDVI平均值上，构建冬小麦面积指数(WWAI)。将每期所设定的权值分别应用到对应时期的影像，即将各期冬小麦NDVI平均值与权值相乘，并取9期乘积加和的平均值，得到的平均值作为冬小麦面积指数(WWAI)。

步骤S6，冬小麦面积指数(WWAI)最优阈值的自适应确定。利用自适应WWAI提取方式，实现不同区域最优WWAI阈值的自动提取。将冬小麦面积指数以1％的比例等分，从0％比例开始至100％结束，可获取101个对应的WWAI阈值，将该阈值用于抽穗期影像中阈值样本点的冬小麦识别，并将识别结果与目视判别结果对比，统计每个阈值对应的准确率，选择101个结果中正确率最高的阈值作为抽穗期冬小麦面积指数的最优阈值；

取2014年5月17日拍摄的冬小麦生育期(即抽穗期)影像作为样本，逐一比较该影像中每个像素点对应的NDVI与所述最优阈值的大小，若像素点对应的NDVI不小于所述最优阈值，则判断该像素点为冬小麦，否则为非冬小麦；统计所述影像中冬小麦像素点的总个数，用所述冬小麦像素点的总个数2491766个×每个像素代表的面积(16m×16m)，计算得到抽穗期内冬小麦的种植面积约为637.89km²。

图3为2014年5月17日安平县及周边冬小麦面积指数影像，从图中可以看出，冬小麦类型的冬小麦面积指数(WWAI)值要明显高于非冬小麦类型。在9期影像中，与原始影像NDVI平均值相比，冬小麦面积指数将冬小麦与非冬小麦的差异倍数从1.43倍增加到2.08倍，说明冬小麦面积指数(WWAI)具有特征性的放大冬小麦NDVI值以及提高冬小麦与非冬小麦类型差异的能力。

图4为利用冬小麦面积指数提取冬小麦面积的结果，在最优阈值下，阈值样本具有最高的总体识别精度，达到了94.0％，冬小麦的识别精度为95.5％。从图中可以看出安平县及周边冬小麦种植分布广泛且均匀，与实际种植情况相吻合。

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