一种细碎化农业区土地覆被与种植结构的遥感识别方法

1.本发明属于农业遥感领域，涉及一种细碎化农业区土地覆被与种植结构的遥感识别方法，采用机器学习算法提取土地利用，不依赖实测样本点，自动训练后得到高精度农业区土地利用栅格化数据。


背景技术：

2.准确及时地获取农业区土地覆被信息对监测农情，调整农业结构，估产和制定粮食政策等均具有重要意义。传统的田间调查方法费时费力，且难以获取大尺度区域土地利用信息。遥感技术能克服单点观测在复杂地表条件下空间代表性差的缺点，为区域土地利用识别提供多时相、多光谱及多角度的地表信息。
3.目前土地利用遥感反演主要分为基于多源数据物候和极化特征及基于实测样本点分类两种方式。基于样本点的识别算法只需输入样本训练集，利用卫星遥感图像，即可快速实现大尺度区域土地利用反演。其中支持向量机(svm)和随机森林等算法因其训练效率高，精度满足应用要求，被广泛应用于土地利用遥感识别，但该方法对样本精度要求高，且应用时需要大量实测样本，时间及金钱成本较高。同时对于光谱及纹理特征较为类似的作物来说，该类算法在识别时会出现作物类型混分的情况。基于多源物候及极化特征算法只需建立特定的分类规则，不需实测样本训练集，即可利用决策树算法实现土地利用反演。这种方法只针对于物候/极化特征差异十分明显的地物识别，且需对原始卫星图像光谱数据进行处理分析，计算得到物候/极化特征用于分类。总而言之，现有的分类算法在应用于卫星遥感图像时，单一分类算法运算时间长，对计算机配置要求高，需要大量实测样本，费时费力，因此难以实现大尺度土地利用识别。
4.在实际应用中，svm、随机森林、决策树等各种分类器都具有其优缺点，因此本发明结合各种算法的优势，针对细碎化农业区不同地物的特点应用不同算法，取长补短，同时提高分类的效率以及精度。


技术实现要素：

5.针对上述技术问题，本发明的目的是提供一种细碎化农业区土地覆被与种植结构的遥感识别方法，通过卫星遥感图像光谱及纹理特征，以及识别农业区和非农业区效果较好的lswi指数，利用机器学习算法实现农业区和非农业区的识别。在此基础上，针对农业区不同作物所具有的物候特征建立分类规则，利用ndvi时序数据计算得到的物候指标作为输入数据，基于决策树算法实现农业区不同作物的进一步划分，最终实现细碎化农业区高精度土地利用遥感识别。
6.为了实现上述目的，本发明提供了如下技术方案：
7.一种细碎化农业区土地覆被与种植结构的遥感识别方法，包括如下步骤：
8.s1、利用支持向量机svm监督学习算法对原始卫星图像进行农业区及非农业区识别：
9.s1.1、采用双线性插值法将多张不同时期不同气象条件的原始卫星图像的每个波段重采样为统一的空间分辨率；选择一张八月初作物生长旺盛的原始卫星图像，并选取该图像中能较好分辨地物的可见光波段以及能判断植物生长状况的红边波段、近红外波段或短波红外波段作为分类识别波段；利用二阶概率统计方法用一个灰色调空间相关性矩阵，基于3
×
3的窗口和64灰度量化级别计算分类识别波段的纹理值；
10.s1.2、利用生育期内多期晴空无云条件下的原始卫星图像，通过公式一得到陆地水体指数lswi时序数据集，采用线性插值的方法将陆地水体指数lswi时序数据集重采样为统一的时间分辨率，并以一预设阈值对数据集中的陆地水体指数lswi时序数据进行重新赋值，大于预设阈值的像元赋值为0，小于预设阈值的像元赋值为1，将重新赋值的陆地水体指数lswi时序数据依次相加求和，获得一个陆地水体指数lswi重分类数据；
[0011][0012]
式中：lswi为陆地水体指数；ρ
nir
表示近红外波段的反射率；ρ
swir
表示短波红外波段的反射率；
[0013]
s1.3、从google earth中选取土地利用样本集作为训练数据，利用支持向量机svm监督分类算法对步骤s1.1得到的分类识别波段及纹理值和步骤s1.2获得的陆地水体指数lswi重分类数据进行训练学习，计算得到农业区及非农业区的土地利用分类结果；
[0014]
s2、应用决策树算法对农业区进行进一步种植结构划分，得到最终农业区高精度土地利用数据：
[0015]
s2.1、将步骤s1得到的土地利用分类结果转变为栅格数据并对其进行重新赋值，将农业区赋值为1，非农业区赋值为nodata；将重新赋值的栅格数据转变为矢量数据，获得只包含农业区的shp矢量数据；
[0016]
s2.2、基于步骤s2.1得到的只包含农业区的shp矢量数据对生育期内多期晴空无云条件下的原始卫星图像进行裁剪，得到只包含农业区的原始卫星图像，并通过公式二计算得到归一化植被指数ndvi时序数据集，采用线性插值的方法将归一化植被指数ndvi时序数据集重采样为统一的时间分辨率；
[0017][0018]
式中：ndvi为归一化植被指数；ρ
nir
表示近红外波段的反射率；ρ
red
表示红光波段的反射率；
[0019]
s2.3、利用savitzky-golay滤波器对步骤s2.2重采样后的归一化植被指数ndvi时序数据集进行平滑处理；基于平滑后的ndvi生育期过程曲线建立物候指标判断标准，对农业区每一个像元进行判断，最终得到物候指标栅格数据；所述物候指标栅格数据中包含生育期起始时间、生育期终止时间和生育期长度三个物候指标；
[0020]
所述物候指标判断标准的建立过程如下：
[0021]
确定在平滑后的ndvi生育期过程曲线的上升阶段、下降阶段以及ndvi设定值，当ndvi＝设定值时，上升阶段所对应的时刻为生育期起始时间，下降阶段所对应的时刻为生育期终止时间，生育期起始时间与生育期终止时间内的时间长度为生育期长度；
[0022]
s2.4、根据不同作物的三个物候指标的值域范围，建立作物分类规则；基于作物分
类规则，利用决策树算法对步骤s2.3获得的物候指标栅格数据中的每一个栅格进行判断，获得农业种植结构分类结果，最终使得每个栅格都被划分为某一类具体的作物类型，从而实现对农业区的作物识别；
[0023]
s3、将步骤s1得到的土地利用分类结果中非农业区与步骤s2获得的农业种植结构分类结果镶嵌到一起，得到一张带有分类标签色彩的农业种植区土地利用及种植结构分类图。
[0024]
所述步骤s1.2中，预设阈值为0.2。
[0025]
所述步骤s1.3中，土地利用分类结果中，农业区为农田，非农业区包括水体、沙丘、居民区和天然地。
[0026]
所述步骤s2.3中，设定值为0.3。
[0027]
所述步骤s2.4得到的农业种植结构分类结果包括小麦、玉米、葵花和果蔬。
[0028]
所述步骤s2.4中，不同作物的三个物候指标的值域范围基于当地实际经验和实地调研分别确定。
[0029]
所述步骤s3中，农业种植区土地利用及种植结构分类图中的每一个像元都有一个分类，不同的值和颜色代表着不同的分类。
[0030]
与现有技术相比，本发明的有益效果在于：
[0031]
1、本发明解决了实地样本获取数量有限，且实测样本质量要求较高，获取困难的问题，在无需实测样本的情况下，保证分类精度及稳定性。
[0032]
2、本发明与传统算法相比，降低了时间成本及金钱成本，且可操作性较强，在应用于大区域大面积高精度遥感图像时，也可用普通硬件设备实现高精度土地利用识别，使得大区域的农业区土地利用自动获取成为可能。
附图说明
[0033]
图1为本发明的细碎化农业区土地覆被与种植结构的遥感识别方法的流程图；
[0034]
图2为本发明实施例的原始卫星遥感图像(在本描述中举例使用sentinel-2历史图像，其他卫星原始图像均可)；
[0035]
图3为本发明实施例的原始卫星遥感图像纹理值计算过程示意图；
[0036]
图4为本发明实施例监督算法支持向量机svm分类的计算过程及经过分类后得到的土地利用图像；
[0037]
图5为本发明选用的物候指标示意图；
[0038]
图6为本发明实施例用于作物分类的决策树分类规则；
[0039]
图7为本发明实施例对农业区进行具体种植结构分类后的图像；
[0040]
图8为本发明实施例最终得到的土地利用及种植结构分类图；
[0041]
图9为分类结果中各个类别的图例。
具体实施方式
[0042]
下面结合附图和实施例对本发明进行进一步说明。
[0043]
如图1所示，一种细碎化农业区土地覆被与种植结构的遥感识别方法，包括如下步骤：
[0044]
s1、利用支持向量机svm监督学习算法对原始卫星图像进行农业区及非农业区识别：
[0045]
s1.1、采用双线性插值法将多张不同时期不同气象条件的原始卫星图像的每个波段重采样为统一的空间分辨率；选择一张八月初作物生长旺盛的原始卫星图像，并选取该图像中能较好分辨地物的可见光波段以及能判断植物生长状况的红边波段、近红外波段或短波红外波段作为分类识别波段；利用二阶概率统计方法用一个灰色调空间相关性矩阵，基于3
×
3的窗口和64灰度量化级别计算分类识别波段的纹理值；
[0046]
s1.2、利用生育期内多期晴空无云条件下的原始卫星图像，通过公式一得到陆地水体指数lswi时序数据集，采用线性插值的方法将陆地水体指数lswi时序数据集重采样为统一的时间分辨率，并以一预设阈值对数据集中的陆地水体指数lswi时序数据进行重新赋值，大于预设阈值的像元赋值为0，小于预设阈值的像元赋值为1，将重新赋值的陆地水体指数lswi时序数据依次相加求和，获得一个陆地水体指数lswi重分类数据；
[0047][0048]
式中：lswi为陆地水体指数；ρ
nir
表示近红外波段的反射率；ρ
swir
表示短波红外波段的反射率。
[0049]
所述预设阈值为0.2。
[0050]
s1.3、从google earth中选取土地利用样本集作为训练数据，利用支持向量机svm监督分类算法对步骤s1.1得到的分类识别波段及纹理值和步骤s1.2获得的陆地水体指数lswi重分类数据进行训练学习，计算得到农业区及非农业区的土地利用分类结果；
[0051]
所述土地利用分类结果中，农业区为农田，非农业区包括水体、沙丘、居民区和天然地。
[0052]
s2、应用决策树算法对农业区进行进一步种植结构划分，得到最终农业区高精度土地利用数据：
[0053]
s2.1、将步骤s1得到的土地利用分类结果转变为栅格数据并对其进行重新赋值，将农业区赋值为1，非农业区赋值为nodata；将重新赋值的栅格数据转变为矢量数据，获得只包含农业区的shp矢量数据；
[0054]
s2.2、基于步骤s2.1得到的只包含农业区的shp矢量数据对生育期内多期晴空无云条件下的原始卫星图像进行裁剪，得到只包含农业区的原始卫星图像，并通过公式二计算得到归一化植被指数ndvi时序数据集，采用线性插值的方法将归一化植被指数ndvi时序数据集重采样为统一的时间分辨率；
[0055][0056]
式中：ndvi为归一化植被指数；ρ
nir
表示近红外波段的反射率；ρ
red
表示红光波段的反射率。
[0057]
s2.3、利用savitzky-golay(s-g)滤波器对步骤s2.2重采样后的归一化植被指数ndvi时序数据集进行平滑处理；基于平滑后的ndvi生育期过程曲线建立物候指标判断标准，对农业区每一个像元进行判断，最终得到物候指标栅格数据；所述物候指标栅格数据中包含生育期起始时间、生育期终止时间和生育期长度三个物候指标；
[0058]
所述物候指标判断标准的建立过程如下：
[0059]
确定在平滑后的ndvi生育期过程曲线的上升阶段、下降阶段以及ndvi设定值，当ndvi＝设定值时，上升阶段所对应的时刻为生育期起始时间，下降阶段所对应的时刻为生育期终止时间，生育期起始时间与生育期终止时间内的时间长度为生育期长度。
[0060]
所述设定值为0.3。
[0061]
s2.4、根据不同作物的三个物候指标的值域范围，建立作物分类规则；基于作物分类规则，利用决策树算法对步骤s2.3获得的物候指标栅格数据中的每一个栅格进行判断，获得农业种植结构分类结果，最终使得每个栅格都被划分为某一类具体的作物类型，从而实现对农业区的作物识别。
[0062]
所述步骤s2得到的农业种植结构分类结果包括小麦、玉米、葵花和果蔬。
[0063]
不同作物的三个物候指标的值域范围基于当地实际经验和实地调研分别确定。
[0064]
s3、将步骤s1得到的土地利用分类结果中非农业区与步骤s2获得的农业种植结构分类结果镶嵌到一起，得到一张带有分类标签色彩的农业种植区土地利用及种植结构分类图。
[0065]
所述农业种植区土地利用及种植结构分类图中的每一个像元都有一个分类，不同的值和颜色代表着不同的分类。
[0066]
实施例
[0067]
s1、利用支持向量机(svm)监督学习算法对原始卫星图像进行农业区及非农业区识别；
[0068]
将sentinel-2原始卫星图像每个波段采用双线性插值法重采样为统一的10m空间分辨率，如图2所示；选择一张八月初作物生长旺盛的原始卫星图像，选取该图像中能较好分辨地物的可见光(红绿蓝)波段，以及能判断植物生长状况的红边波段、近红外波段以及短波红外波段作为分类识别波段。利用二阶概率统计方法用一个灰色调空间相关性矩阵，基于3
×
3的窗口和64灰度量化级别对最优波段的纹理值进行计算，最终得到基于八个纹理滤波(均值、方差、协同性、对比度、相异性、信息熵、二阶矩和相关性)依次计算得到的纹理值，如图3所示；
[0069]
利用生育期内多期晴空无云条件下的sentinel-2原始卫星图像，经过公式一计算得到陆地水体指数lswi时序数据集，采用线性插值的方法将陆地水体指数lswi数据集重采样为统一的5天时间分辨率，并以0.2为阈值对其进行重新赋值，大于0.2的像元赋值为0，小于0.2的像元赋值为1，将重新赋值的陆地水体指数lswi时序数据依次相加求和，获得一个陆地水体指数lswi重分类数据；
[0070]
从google earth中人工选取土地利用样本集(居民地、沙丘、荒地、水体和农田)作为训练数据，利用支持向量机svm监督分类算法对最优波段及纹理值以及lswi重分类数据进行训练学习，计算得到农业区及非农业区土地利用分类结果，如图4所示。
[0071]
s2、应用决策树算法对农业区进行进一步种植结构划分，得到最终农业区高精度土地利用数据；
[0072]
将土地利用分类结果转变为栅格数据并对其进行重新赋值，将农业区赋值为1，非农业区统一赋值为nodata，将重新赋值的栅格数据转变为矢量数据，得到农业区的shp矢量数据；
[0073]
基于农业区矢量数据对生育期内多期晴空无云条件下的原始卫星图像进行裁剪，得到只包含农业区的原始卫星图像，通过公式二计算得到归一化植被指数ndvi时序数据集，将时序数据集采用线性插值的方法重采样为统一的5天时间分辨率；
[0074]
利用savitzky-golay(s-g)滤波器对生育期ndvi时序数据集进行平滑处理。如图5所示，基于平滑后的ndvi生育期过程曲线，确定在曲线的上升阶段以及下降阶段，当ndvi＝0.3时，上升阶段以及下降阶段所对应的时刻分别为生育期起始时间及生育期终止时间，生育期起始时间与生育期终止时间内的时间长度为生育期长度，基于上述物候指标判断标准对农业区每一个像元进行判断，最终得到包含生育期起始时间、生育期终止时间以及生育期长度三个物候指标的物候指标栅格数据，用于进一步的作物分类。
[0075]
基于当地实际经验和实地调研，分别确定不同作物的三个物候指标的值域范围，建立作物分类规则，如图6所示。基于分类规则，利用决策树算法对物候指标栅格数据中的每一个栅格进行判断，最终使得每个栅格都被划分为某一类具体的作物类型，从而实现对农业区的作物识别，得到如图7所示的农业区种植结构分类结果。
[0076]
将土地利用结果中的非农业区与种植结构识别结果镶嵌到一起，如图8所示，进而得到农业区高精度土地利用分类结果，分类的结果是得到带有分类标签色彩的一张图像，包括水体、沙丘、居民区、天然地非农业区土地利用类型以及小麦、玉米、葵花和果蔬种植结构分类结果，各类地物及作物的图例如图9所示。