本发明涉及遥感影像特征提取领域，具体涉及一种多维度多尺度遥感与地形特征融合的耕地非农化提取方法。
背景技术：
1、我国人多地少，耕地是最为宝贵的资源之一。近年来，严厉查处违法违规占用耕地，严格管控耕地转为其他农用地，强化耕地用途管制。因此利用遥感技术精准解译耕地，实现耕地非农化智能监测，对于耕地保护和粮食安全具有重要作用。
2、目前耕地非农化遥感监测方法主要分为目视判读和自动化算法这两种。耕地非农化监测遥感目视判读是通过前后期影像，依据人工目视全域查找耕地非农化图斑，需要一定的先验知识，人力和时间成本巨大，在遥感大数据时代面对高频次、大范围的耕地保护监测，目视判读显然无法满足实际需求。自动化耕地非农化监测算法则是利用遥感智能识别技术从影像上自动识别耕地及耕地变化区域，主要包括耕地非农化区域智能变化检测、耕地先智能解译再找非农化区域这两种方式。耕地非农化区域智能变化检测需要先制作海量优质耕地非农化的变化样本，再构建变化检测模型，最后利用训练后的模型实现自动查找耕地非农化图斑，该方法依赖于采集制作海量优质耕地非农化的变化样本，整个周期较长，难以满足实际业务的急切需求。耕地先智能解译再找非农化区域是通过对前后期影像进行解译提取耕地，再进行耕地变化及非农化判断，或者仅仅对后期影像解译耕地，再利用历史成果耕地矢量，得到非农化图斑。
3、特别是在西南山区，地形起伏大、地表垂直分异严重，西南山区耕地在遥感影像上特征多变、块小零散，耕地与林草地、种植园地、水草覆盖坑塘等的特征较为相似，“同物异谱”和“异物同谱”问题严重，造成目前主流的单一运用高分辨率影像的耕地自动化识别方法精度不高。从目前的方法技术来看，仅仅利用高分辨率影像的空间特征，难以精准提取复杂山区环境的耕地。
技术实现思路
1、针对现有技术的不足，本发明提供一种多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，解决复杂山地环境耕地精准识别和耕地非农化智能监测的难题，该方法监测结果更准确，具有实操性、便捷性和可行性。
2、为达到上述目的，本发明采用的技术方案如下：
3、一种多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，包括如下步骤：
4、步骤1、划分训练区和测试区，获取训练区的高分辨率影像、高光谱影像、数字高程模型dem数据、以及耕地真值矢量数据，获取测试区的高分辨率影像、高光谱影像、数字高程模型dem数据，以及已调查监测的耕地矢量数据；
5、步骤2、先对所述高分辨率影像利用嵌入边缘检测canny的简单非迭代聚类snic算法进行超像素分割，再利用高斯混合算法gmm，结合高分辨率影像对超像素对象进行聚类，获得超像素地物对象；
6、步骤3、以所述步骤2获得的超像素地物对象为单位，提取各超像素地物对象范围内的特征，所述特征至少包括高分辨率影像的空间域特征、高分辨率影像的频率域高频特征、高光谱影像的光谱域特征及dem数据的地形特征；
7、步骤4、将所述训练区的高分辨率影像、高分辨率影像、dem数据和耕地真值矢量数据根据所述步骤2、所述步骤3中提取的各特征组合成每个超像素地物对象的特征向量，根据超像素地物对象内类别标签像素面积占比最大的确定超像素地物对象标签，构建超像素耕地对象训练样本数据；
8、步骤5、将所述步骤4中构建的超像素耕地对象训练样本数据输入随机森林分类器，设置随机森林算法相关超参数，训练耕地对象识别模型；
9、步骤6、将所述测试区的高分辨率影像、高光谱影像和dem数据根据所述步骤2、所述步骤3中提取的各特征组合成每个超像素地物对象的特征向量，再将所述特征向量输入所述步骤5中耕地对象识别模型进行预测，得到所述测试区的每个超像素对象的类别，提取测试区耕地。
10、步骤7、基于所述步骤6提取的测试区耕地和所述已调查统计的耕地矢量数据，利用空间叠置分析自动提取耕地变为非耕地的图斑，先对所述图斑通过图斑面积、长宽比、紧致度阈值删除掉细碎狭长图斑，再利用图斑边界简化圆滑、孔洞填充、边界指数、形状指数对所述图斑进行消除孔洞、锯齿和尖角过滤、筛选出候选耕地非农化图斑矢量，最后依据实际业务设定面积阈值，获得耕地非农化图斑。
11、进一步地，所述步骤1中划分训练区和测试区，获取训练区的高分辨率影像、高光谱影像和数字高程模型dem数据、以及耕地真值矢量数据，获取测试区的高分辨率影像、高光谱影像、数字高程模型dem数据和已调查统计的耕地矢量数据，包括：
12、选取具有典型山区耕地特色的区域划分出训练区和测试区，在训练区和测试区采集获取相同区域时相相近的高分辨率影像和高光谱遥感影像，对所述高分辨率影像和高光谱遥感影像进行配准和几何精校正，并采集对应区域dem数据，所述耕地真值矢量数据是将与所述高分辨率影像对应时相的调查监测成果数据提取的耕地范围矢量，所述耕地范围矢量为耕地真值矢量；所述已调查监测的耕地成果矢量数据是以当前高分辨率影像时相之前已调查监测的耕地成果矢量数据。
13、进一步地，所述步骤2中先对所述高分辨率影像利用嵌入边缘检测canny的snic算法进行超像素分割，具体过程如下：
14、步骤2.1、对高分辨率的红、绿、蓝三个波段，记为red、green、blue，分别利用canny算子提取边缘，得到三个波段的边缘检测结果，记为redcanny、greencanny、bluecanny；
15、步骤2.2、将高分辨率的红、绿、蓝三个原始波段和三个波段的边缘检测结果进行乘性和加性运算，得到带边缘信息的红、绿、蓝三个波段，记为rededge、greenedge、blueedge，计算公式如下：
16、
17、步骤2.3、将rededge、greenedge、blueedge组成的三波段图像，转为cielab颜色空间，l为cielab颜色空间下的亮度通道，a、b为cielab颜色空间下的两个颜色通道，a为从洋红色至绿色的范围，b为从黄色至蓝色的范围；
18、步骤2.4、将cielab颜色空间下的l、a、b三通道作为snic算法的处理对象，依据高分辨影像的分辨率、耕地最小地块面积和待分割区域的面积设置snic算法分割输出的超像素数量为n，平衡颜色优先性和空间优先性的权重为c，根据影像分辨率和待解译地物的面积尺寸调整n和c的值生成紧凑、近似均匀的超像素，使得超像素的边界与地物的边界尽可能吻合，按snic算法分割输出的超像素数量n采用如下公式：
19、
20、其中s为待分割高分辨率影像区域的面积，单位为平方米；pnmin为依据影像分辨率和耕地最小地块面积设定的像素个数。
21、进一步地，所述步骤3中以所述步骤2获得的超像素地物对象为单位，提取各超像素地物对象范围内的特征，具体包括：以超像素地物对象为单位，
22、a、利用灰度共生矩阵glcm计算每个超像素地物对象范围内高分辨率影像的空间域特征，所述空间域特征至少包括角二阶矩asm、对比度con、相关性cor、熵ent和同质度hom五种空间纹理特征，记为meanasm,meancon,meancor,meanent,meanhom；
23、b、利用傅里叶变换获得频谱图，经过布特沃斯高通滤波器过滤低频信息，再经过傅里叶逆变换得到每个超像素地物对象范围内高分辨率影像的频率域高频特征，记为meanhighfrequency；
24、c、利用主成分分析获得的前三个主成分、改良土壤调整植被指数msavi、耕地敏感响应波段特征得到每个超像素地物对象范围内高光谱影像的光谱域特征，记为meanpca1,meanpca2,meanpca3,meanmsavi,meanbandtop1,meanbandtop2,meanbandtop3；
25、d、利用arcgis软件的栅格表面分析，计算每个超像素地物对象范围内dem数据的地形特征，所述地形特征至少包括坡度slope、坡向aspect、山体阴影hillshade、地形起伏度undulation和地表粗糙度roughness特征，记为meanslope,meanaspect,meanhillshade、meanundulation、meanroughness。
26、进一步地，所述利用灰度共生矩阵glcm计算每个超像素地物对象范围内高分辨率影像的空间域特征，具体过程如下：
27、首先、将高分辨率影像多波段取平均值转为单波段灰度图像，并压缩灰度级至l(l一般为16、32、64))，以每个像素的w*w邻域为窗口，w根据影像分辨率和待解译地物的面积确定，计算横向、纵向和两个斜向这四个不同角度方向、一个像素间隔的glcm，每个邻域窗口得到四个glcm，glcm为l行l列的矩阵；
28、其次、每种纹理特征，对依据四个glcm计算出的四个值取平均作为该邻域窗口中心像素的纹理特征值，滑动窗口遍历所有像素，得到该种纹理特征图像，最终获得五种纹理特征图像；
29、所述基于glcm计算五种纹理特征的公式分别为：
30、
31、其中，p(i,j)为glcm中第i行、第j列的值，μ1和μ2为均值，σ1和σ2为方差，计算公式如下：
32、
33、最后、对每个超像素地物对象，计算其内部每种纹理特征图像的像素平均值，记为meanasm,meancon,meancor,meanent,meanhom，获得该超像素地物对象的五种空间纹理特征值。
34、进一步地，所述利用傅里叶变换获得频谱图，经过布特沃斯高通滤波器过滤低频信息，再经过傅里叶逆变换得到每个超像素地物对象范围内高分辨率影像的频率域高频特征，具体过程如下：
35、首先、将高分辨率影像多波段取平均值转为单波段灰度图像f(x,y)，利用傅里叶变换得到频谱图f(u,v)：
36、
37、其中x和y是单波段灰度图像横向纵向像素位置，变量u和v是频谱图中频率变量，j为虚数单位，m和n为单波段灰度图像f(x,y)的横向与纵向像素尺寸；
38、其次、对频谱图f(u,v)，经过布特沃斯高通滤波器过滤低频信息得到g(u,v)：
39、
40、其中d0为频谱图中心为圆点的半径阈值，d(u,v)为频谱图中任意位置距离中心的距离；
41、然后、对布特沃斯高通滤波后的频谱图g(u,v)，进行傅里叶逆变换得到高频信息图f(x,y)h：
42、
43、其中变量u和v是频谱图中频率变量，j为虚数单位，x和y是高频信息图f(x,y)h横向纵向像素位置；
44、最后、对每个超像素地物对象，计算其内部高频信息图f(x,y)h的平均值，记为meanhighfrequency，获得该超像素地物对象的频率域高频特征值。
45、进一步地，所述利用主成分分析获得的前三个主成分、改良土壤调整植被指数msavi、耕地敏感响应波段特征得到每个超像素地物对象范围内高光谱影像的光谱域特征，具体包括：
46、首先、对高光谱影像进行主成分分析pca，获取前三个主成分pca1、pca2、pca3；
47、其次、利用高光谱影像的红波段red、近红外波段nir计算msavi；
48、
49、然后、在高光谱影像上选取耕地和背景的纯像元区域，计算每个波段的耕地与背景之间的像素差异性pd，所有波段按pd值降序排列，最高的前三个波段为耕地敏感响应波段；
50、pd＝1-cos(clhistogram,bghistogram)
51、其中clhistogram和bghistogram为某波段的耕地与背景的像素灰度直方图，cos为余弦距离；
52、最后、对每个超像素地物对象，计算其内部pca1、pca2、pca3、msavi、和三个耕地敏感响应波段的平均值，记为meanpca1、meanpca2、meanpca3、meanmsavi、meanbandtop1、meanbandtop2、meanbandtop3，获得该超像素地物对象的光谱域光谱特征值。
53、进一步地，所述步骤4中根据所述步骤2、所述步骤3中提取的各特征组合成每个超像素地物对象的特征向量features＝[meanasm,meancon,meancor,meanent,meanhom,meanhighfrequency,meanpca1,meanpca2,meanpca3,mean msavi,meanbandtop1,meanbandtop2,meanbandtop3,meanslope,meanaspect,meanhillshade、meanundulation、meanroughness]。
54、进一步地，所述步骤4中根据超像素地物对象内类别标签像素面积占比最大的确定超像素地物对象标签，具体包括：
55、所述超像素对象的标签label根据该对象内类别标签面积占比最大确定，1表示为耕地类，0表示为背景类，公式如下：
56、
57、其中，s耕地和s背景分别为对象内耕地和背景的像素面积。
58、进一步地，所述步骤6中还包括将得到所述测试区的每个超像素对象的类别进行平滑滤波，提取测试区耕地。
59、进一步地，所述步骤7，具体过程如下：
60、步骤7.1、以已调查监测的耕地矢量数据成果为基准，将利用随机森林模型在测试区解译的耕地矢量与其进行空间擦除分析，获得在已调查监测的耕地矢量数据成果上是耕地、而在解译耕地矢量中不是耕地的图斑，得到初步疑似耕地非农化图斑矢量shapedata1；
61、步骤7.2、对于初步疑似耕地非农化图斑矢量shapedata1，计算每个图斑的面积s、长宽比b、紧致度c，设定阈值s0、b0、c0，若s<s0或b>b0或c<c0，则删除掉细碎狭长图斑，再填充图斑内的细小孔洞，对每个图斑内孔洞面积sk<sk0的进行填充，得到疑似耕地非农化图斑矢量shapedata2；
62、所述图斑的长宽比b计算式为：
63、
64、其中，len和wid分别为图斑最小外接矩形的长度和宽度；
65、所述图斑的紧致度c的计算式为：
66、
67、其中，s和p分别为该图斑的面积和周长；
68、步骤7.3、对于疑似耕地非农化图斑矢量shapedata2，利用道格拉斯-普克算法douglas–peucker进行图斑简化，消除图斑边界的锯齿，再利用贝塞尔插值方法bezierinterpolation进行图斑边界圆滑，消除图斑的尖角，然后计算每个图斑的边界指数borderindex、形状指数shapeindex，设定阈值bi0、si0，若borderindex<bi0且shapeindex<si0，则保留满足条件的图斑，得到候选耕地非农化图斑矢量shapedata3；
69、所述边界指数borderindex和所述形状指数shapeindex的计算公司分别为：
70、
71、其中，s和p分别为该图斑的面积和周长，extentp为为图斑最小外接矩形的周长；
72、步骤7.4、根据实际中满足一定面积的耕地减少图斑才是耕地非农化图斑这一要求，对于候选耕地非农化图斑矢量shapedata3，依据业务需求选取阈值st0，按照s>st0得到耕地非农化图斑成果矢量。
73、本发明的显著效果是：
74、本发明鉴于对象级特征提取和融合的相关方法，先通过利用嵌入边缘检测canny的简单非迭代聚类snic分割高分辨率影像得到超像素，再运用高斯混合算法gmm聚类超像素获得地物对象，对高分辨率影像进行分割和聚类，获得准确的耕地对象，且相较于已有对象分割方法，分割出对象的边界与实际地物更吻合。
75、然后以超像素耕地对象为处理单元，综合顾及空间域特征、频率域特征、光谱域特征及地形特征，从高分辨率影像提取超像素地物对象的空间域特征、频率域高频特征，从高光谱影像提取超像素地物对象的光谱特征，从数字高程dem提取超像素地物对象的地形特征，并将获得的多维度多尺度特征组合为特征向量；以随机森林为分类器，利用随机森林分类算法构建超像素耕地对象的提取模型，进而精准解译耕地，实现了多维度多尺度遥感与地形特征融合的耕地提取及耕地非农化图斑提取，提高了复杂山区耕地识别精度。本方法充分考虑到耕地和其他地类的光谱特性存在差异，将高光谱影像和高分辨率影像协同，充分利用高光谱影像的光谱特征和高分辨率影像的空间结构特征，以提高耕地的识别精度。另外西南山区的耕地具有较强的地形特征，如坡耕地、槽谷耕地具有显著的坡度、坡向、山体阴影特征，因此地形特征对于西南山区耕地提取和耕地非农化监测将有重要促进作用。本方法利用在超像素对象层面进行特征提取与融合，有效解决了复杂山地环境耕地精准识别和耕地非农化智能监测的难题。
76、最后以调查监测统计的耕地矢量数据成果为基准，如：国土三调耕地历史成果矢量数据；对随机森林模型解译耕地的结果矢量进行擦除叠加分析，自动提取耕地变为非耕地的图斑，实现西南山区耕地非农化监测，并且对该图斑通过优化、筛选、面积阈值过滤出耕地非农化图斑，获得的数据更为准确。本发明具有显著的实操性、便捷性、可行性等优势。
1.一种多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，其特征在于，包括如下步骤：
2.根据权利要求1所述的多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，其特征在于：所述步骤1：划分训练区和测试区，获取训练区的高分辨率影像、高光谱影像和数字高程模型dem数据、以及耕地真值矢量数据，获取测试区的高分辨率影像、高光谱影像、数字高程模型dem数据以及已调查监测的耕地成果矢量数据，包括：
3.根据权利要求1所述的多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，其特征在于：所述步骤2中先对所述高分辨率影像利用嵌入边缘检测canny的snic算法进行超像素分割，具体过程如下：
4.根据权利要求1所述的多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，其特征在于：所述步骤3中以所述步骤2获得的超像素地物对象为单位，提取各超像素地物对象范围内的特征，具体包括：以超像素地物对象为单位，
5.根据权利要求4所述的多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，其特征在于：所述利用灰度共生矩阵glcm计算每个超像素地物对象范围内高分辨率影像的空间域特征，具体过程如下：
6.根据权利要求4所述的多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，其特征在于：所述利用傅里叶变换获得频谱图，经过布特沃斯高通滤波器过滤低频信息，再经过傅里叶逆变换得到每个超像素地物对象范围内高分辨率影像的频率域高频特征，具体过程如下：
7.根据权利要求4所述的多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，其特征在于：所述利用主成分分析获得的前三个主成分、改良土壤调整植被指数msavi、耕地敏感响应波段特征得到每个超像素地物对象范围内高光谱影像的光谱域特征，具体过程：
8.根据权利要求4所述的多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，其特征在于：所述步骤4中根据所述步骤2、所述步骤3中提取的各特征组合成每个超像素地物对象的特征向量features＝[meanasm,meancon,meancor,meanent,meanhom,meanhighfrequency,meanpca1,meanpca2,meanpca3,meanmsavi,meanbandtop1,meanbandtop2,meanbandtop3,meanslope,meanaspect,meanhillshade、meanundulation、meanroughness]。
9.根据权利要求1所述的多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，其特征在于：所述步骤6中还包括将得到所述测试区的每个超像素对象的类别进行平滑滤波，提取测试区耕地。
10.根据权利要求1所述的多维度多尺度遥感与地形特征融合的耕地非农化智能监测方法，其特征在于：所述步骤7，具体包括：
