一种基于三角网的面向对象耕地信息自动提取方法与流程

本发明涉及一种耕地信息自动提取方法，特别是涉及一种面向对象的高效非监督耕地信息提取方法。

背景技术：

耕地信息的自动提取对于我国农村土地资源管理和土地规划有着至关重要的作用。当前,耕地信息制图的相关技术方法主要以遥感技术为主。在不同的空间尺度上，已经有很多关于耕地信息提取的研究，主要分为局部、地区和全球区域。局部区域的研究主要涉及单个地区或灌区盆地或特定区域；地区范围主要包括大河流域和大陆性的区域；而全球尺度则主要涉及全球范围的耕地信息制图。在不同空间尺度下，为了达到最佳的耕地信息提取效果，不同的传感器数据，需要采用不同的制图方式，根据不少学者的研究成果总结，局部区域的耕地信息提取采用Landsat TM/ETM+，SPOT，LISS，ASTER，CBERS数据能够达到较好的性价比，信息提取方法主要涉及相片解译、影像数字运算、影像分类、分割、影像融合。地区范围主要采用Landsat TM/ETM+，MODIS，MERIS，AVHRR，SPOT VGT数据，通过时间序列分析方法、监督或非监督分类方法、掩膜方法识别耕地信息。全球尺度主要采用MODIS，MERIS，AVHRR，SPOT VGT数据，信息提取方法主要涉及非监督聚类、针对时间序列数据的机器学习算法，以及采用其它辅助数据(包括统计数据，地面实时数据等)。

以上研究使用的遥感影像的空间分辨率较低，当识别小范围和重要地区的耕地往往比较困难，这就需要辅以更高分辨率的影像数据。特别地，中国始于“十一五”的高分辨率对地观测系统专项的实施，使得高空间分辨率影像数据(<10m)获取技术进一步得到发展。然而，由于高空间分辨率影像光谱信息更加丰富，异物同谱现象严重，很难用基于像素的分类方法提取耕地。另外，耕地不同物候的差异，例如，在无人机高分影像中，旱地与水田的色彩相差很大，旱地为黄色，水田为绿色，很难将两者同时提取。因此，高分影像中的耕地专题信息的提取难度更大。

针对高分影像的耕地专题信息提取，已有不少学者做了相关研究工作，Liu等2008年在“GIS与建筑环境联合会议”撰文“High quality prime farmland extraction pattern based on object-oriented image analysis”，利用航空影像数据(3m),通过边缘检测分割影像，随后计算对象特征值，最后利用C4.5构建耕地信息提取决策树，实现高质量标准农田的探测。Lu等2007年在《International Journal of Remote Sensing》撰文“Comparison between several feature extraction/classification methods for mapping complicated agricultural land use patches using airborne hyperspectral data”，利用航空高光谱数据(2m)提取复杂农业用地覆盖区域。Duro等2012年在《Remote Sensing of Environment》撰文“A comparison of pixel-based and object-based image analysis with selected machine learning algorithms for the classification of agricultural landscapes using SPOT-5HRG imagery”，使用SPOT影像，采用多尺度分割技术，分别利用随机森、支持向量机、决策树等监督分类方法，实现耕地信息的提取。诸如此类的基于分割后的对象，利用监督方法提取耕地信息的方法还很多，但是高分影像提供更加丰富的光谱和纹理特征信息，以上方法多为黑箱操作，并没有充分利用耕地地块信息自身独有的特点，导致提取效果不好。

目前，专门针对高分影像开展的非监督耕地信息提取方法很少，自动化程度并不高，对分割后的地理对象的语义信息利用并不充分。例如，Sun和Xu 2009年在“Transactions of the CSAE”撰文“Comer extraction algorithm for high-resolution imagery of agricultural land”，利用针对Quickbird全色波段影像，先对影像数据进行分割，结合角特征点所提供的形状信息(耕地外形规则，一般具有4个或4个以上的角特征点)提取分割后的耕地和坑塘等形状规则的农用地信息。该研究为高分辨率全色影像耕地信息提取开辟了一个新的思路，一定程度上利用高分辨率影像数据中耕地地块的几何和纹理特征，以及耕地地块的独特性质进行算法设计，从而快速高效的提取耕地信息。但是其仅利用了单一的形状特征，对复杂环境的提取效果并不好。考虑到耕地信息具有明显的集中连片特征，且由于高分影像分割后的耕地对象较大，而对于集中居民地区域，由于内部地物丰富，包括房屋、林地、道路、水泥地等，形成形态各异的分割对象，相对耕地区域，该分割区域看上去比较破碎。基于这一特征，我们考虑聚类的方法将破碎居民地、林地等地物剔除，从而自动提取耕地信息。

技术实现要素：

本发明要解决的技术问题是：克服现有技术的上述不足，提供一种基于三角网的面向对象耕地信息自动提取方法，其结合了面向对象影像分类、三角网分析和AUTOCLUST聚类算法，基于高分遥感影像，自动、高效、准确地进行耕地信息的提取。

为了解决以上技术问题，本发明提供的一种基于三角网的面向对象耕地信息自动提取方法，包括以下步骤：

步骤1、利用多尺度分割法对高空间分辨率遥感影像进行分割；

步骤2、剔除长条状分割对象；

步骤3、提取剩余分割对象的中心点；

步骤4、利用中心点构建Delaunay三角网；

步骤5、对步骤4获得的Delaunay三角网进行剥皮操作，剥离Delaunay三角网的外三角；

步骤6、对步骤5中剥皮后的Delaunay三角网进行聚类，得到三角网节点分布密集和分布稀疏的两种类型的若干个聚类结果；

步骤7、利用Vornoio图约束对聚类结果进行修正，修正方法如下：利用Delaunay三角网节点构建Voronoi图，计算Voronoi图中Voronoi多边形面积A，Voronoi多边形面积小于下限阈值时，初步确定该Voronoi多边形对应的中心点属于破碎区域，Voronoi多边形面积大于上限阈值时，初步确定该多边形对应的中心点为耕地区域；将判定为破碎区域的中心点删去或归入节点分布密集的聚类中，将判定为耕地区域的中心点归入节点分布稀疏的聚类中；

步骤8、删除节点分布密集的聚类所对应的中心点；

步骤9、若所述步骤1分割后得到分割对象的最大光谱方差Max.diff大于预设阈值，则判定该分割对象为林地，剔除该分割对象所对应的中心点，最终保留下来的所有中心点所对应的分割对象为耕地，得到耕地提取结果。

本发明还具有以下进一步的特征：

1、步骤1中，使用影像分析软件对高空间分辨率影像进行反复的分割尝试，并观察分割对象的分布特征，选择对耕地地块分割较好的分割尺度进行影像分割。

2、矩形拟合度特征数值小于0.5的分割对象判定为长条状分割对象，并剔除。

3、步骤3中，计算剩余分割对象的多边形面积，并根据面积计算中心点的坐标，所述中心点为分割对象多边形的质心。

4、步骤5中，选择距离参数d，删除Delaunay三角网中三角形边长大于d的边，d的取值范围为[100m,150m]。

5、在步骤7中，下限阈值为上限阈值为其中，Voronio_mean为生成的泰森多边形的面积的平均值，和为系数，且

6、

7、在步骤7或步骤8执行之前，进行如下判断，如果第j个聚类满足则将该聚类归入节点分布稀疏的聚类，其中表示第j个聚类的三角形面积的均值，表示所有聚类的三角形集合的面积均值的平均值，表示第j个聚类的三角形面积集合的方差。

8、预设阈值为2倍所有分割对象的最大光谱方差均值与所有分割对象的最大光谱方差标准差之和。

9、步骤6中采用AUTOCLUST方法进行聚类。

本发明基于Delaunay三角网结合AUTOCLUST聚类算法和多种优化约束方法实现对耕地信息的自动提取。本发明能够在面向对象影像分类的基础上将三角网和聚类算法结合，自动准确地从中选取可靠的耕地信息，同时，简化了耕地信息提取过程中人机交互强度，保证了耕地提取的整体性。

本发明的步骤1中使用影像分割技术，利用高空分辨率影像，通过反复的分割尝试并观察分割对象的分布特征，发现耕地信息具有明显的集中连片特征，分割后的耕地块都比较大，且由于低空无人机影像的高分辨率特点，对于集中居民地区域，其内部地物丰富，包括房屋、林地、道路、水泥地等，分割后形成形态各异的分割对象，相对耕地区域，该分割区域看上去比较破碎。另一方面，对于中国广大农村，居民地区域呈棋盘状分布在耕地中，一般地，集中居住居民地和道路区域是耕地信息提取的主要干扰因素，而道路和沟渠等对象具有明显的长条状特征，与耕地信息不符，因此在步骤2中直接剔除此类对象。

进一步地，对于非长条形的地物，本发明基于三角网分析对非耕地的干扰因素进行逐步剔除，在这个过程中将引入AUTOCLUST聚类、Voronio图约束和方差约束，多种聚类约束算法一步步提高耕地提取的精度。

可见，本发明以三角网为基础，使用AUTOCLUST聚类为主要提取算法，利用长短边聚类实现聚类信息的传递，具体地，采用局部和全局方差表征长短边信息，删除三角网中的短边集合和长边集合并进一步扩大范围。

相比现有技术，本发明充分结合了耕地遥感影像分割后的特征信息以及三角网分析，形成了一种自动、高效、准确的耕地信息提取方法。具体创新点及有益效果如下：

第一，本发明能够充分利用Delaunay三角网的能自动适应不同区域点的密度变化的特征，结合基于图论的AUTOCLUST算法，利用长短边聚类实现聚类信息的传递，具体地，采用局部和全局方差表征长短边信息，从而能够使用该方法实现三角网的聚类。

第二，本发明在AUTOCLUST算法的基础上，使用泰森多边形约束进一步处理被忽略三角网中三角形面积的特征，根据Delaunay三角网得到其对偶图Voronoi图，进一步通过计算泰森多边形面积的方法，表征点的分布密度，解决使用AUTOCLUST聚类算法时出现的欠聚类的现象；使用聚类三角网面积均值约束，解决AUTOCLUST聚类中较大三角网的误聚类问题；使用最大方差约束，进一步优化算法，剔除平原地区零星林地。

综上，本发明充分挖掘高分遥感影像的语义信息，提出了一种以三角网为基础，结合AUTOCLUST聚类的耕地提取方法，自动、高效、准确地进行耕地信息的提取。整个方法，自动化程度较高、精度较高，稳定性较高。不仅提高了高分遥感影像的耕地信息提取效率，同时，保证了耕地提取的整体性。

附图说明

下面结合附图对本发明作进一步的说明。

图1为本发明方法流程图。

图2a为无人机试验原始影像。

图2b为无人机试验原始影像分割尺度为50的分割结果。

图2c为无人机试验原始影像分割尺度为100的分割结果。

图2d为无人机试验原始影像分割尺度为180的分割结果。

图3a为Delaunay三角网示意图。

图3b为剥离Delaunay三角网的外三角示例图。

图4a为图2d分割结果构建的三角网。

图4b为AUTOCLUST聚类结果。

图4c为Voronio图下限阈值约束的聚类图。

图4d为Voronio图上限阈值约束的聚类图。

图5a为使用提出的方法的耕地提取结果。

图5b为使用图像分析软件eCognition进行耕地提取得到的结果。

具体实施方式

下面根据附图详细说明本发明，本发明的目的和效果将变得更加明显。

如图1所示为本发明实施例基于三角网的面向对象耕地信息自动提取方法的流程图，本实施方法包括以下步骤：

步骤1、利用多尺度分割法对高空间分辨率影像进行分割。本实施例中，使用eCognition对高空间分辨率影像进行反复的分割尝试，并观察分割对象的分布特征，选择对耕地地块分割较好的分割尺度进行影像分割。本发明建议分割尺度的取值范围为[150,200]，以便突出高分影像中的居民地与耕地对应的分割对象的差异。

图2为三景无人机试验影像。本步骤中，利用影像分析软件eCognition进行多尺度分割，分割参数Shape为0.4，Smoothness为0.5。分别实施三个尺度的分割(包括尺度50、100、180)，通过观察不同尺度的分割结果，选择对耕地地块分割较好的尺度进行后续分析。容易发现，随着分割尺度的增加，分割对象的面积逐渐增大。相对分割尺度50和100，当分割尺度在180时，能够更好将分辨率0.2米的无人机影像数据中具有相同像素特征的耕地地块分割出来。

图2a为原始影像，图2b是分割尺度为50的分割结果，图2c是分割尺度为100的分割结果，图2d是分割尺度为180的分割结果。

步骤2、剔除长条形状的分割对象(道路，沟渠等)。本实施例中，将矩形拟合度特征(Rectangular fit)数值小于0.5的分割对象剔除，矩形拟合度特征(Rectangular fit)可用于剔除道路、沟渠等与耕地地块特征明显差异的长条状分割对象。

成功对无人机影像数据进行分割后，将使用三角网聚类剔除破碎居民地区域，以提取耕地信息。首先利用分割对象提取中心点(步骤3)，提取的中心点个数和该区域分割对象数目相同。图2d耕地虽然分布较集中，但是难免夹杂有道路、沟渠、田埂等地物，在高分辨率的影像中，这些地物都能被分割出来，且它们都属于长条形状地物，与耕地的矩形有较大差别。在步骤3提取的中心点中，多了这些地物的中心点，它们会影响方法精度，为提高方法可靠性，在构建三角网之前，需要将这些条状地物剔除。所以，首先使用空间特征矩形拟合度(Rectangular fit)阈值剔除条状地物(Rectangular fit<0.5)对应的中心点。

步骤3、提取剩余分割对象的中心点。本例中，首先计算剩余分割对象的多边形面积，并根据面积计算其中心点的坐标。

为计算分割对象多边形的中心点(重心或质心)，首先计算多边形的面积，如果某个分割对象的分割线有N个节点(xi,yi)组成，其中0≤i≤N-1。假设最后一个节点与第一个节点相同，那么多边形封闭。那么，多边形的面积计算公式表示为：

其中xN＝x0，则表示多边形封闭。在计算得到面积以后，多边形中心点(这里实际上表示重心或质心)的坐标(Cx，Cy)能够通过下式计算：

步骤4、利用步骤3提取的中心点构建Delaunay三角网,Delaunay三角网如图3a所示。

Delaunay三角网是相互不重叠且相互邻接的三角形集合，一般需满足两个准则：Delaunay三角网中的各三角形的外接圆中不存在其它点。在所有可能的三角网中，Delaunay三角网中的各三角形的最小角一般最大。根据以上特征，许多研究提出不同的三角网构建方法，流行的方法包括：分治法、扫描线算法、增量算法、快速增量构造算法、凸包算法。本发明重点不在于三角网构建算法，本实施例利用成熟的分治算法(Divide-and-Conquer Algorithm)构建Delaunay三角网，首先，尽量划分点集到足够小的子集，对不同子集构建子三角网，随后整合各子三角网，最后，利用局部优化法保证构建的三角网为Delaunay三角网，这样构建三角网的优点是时间效率较高，缺点是需要的递归运算量大，对内存空间的要求较高。

步骤5、对步骤4获得的Delaunay三角网进行剥皮操作，剥离Delaunay三角网的外三角。本例中，利用步骤4中获得的Delaunay三角网，选择距离参数d，距离参数d应当不小于三角网内部最长边，本发明建议d的取值范围为[100m,150m]，删除三角形边长大于d的边，并将仅与一个三角形关联的边作为点群的边界。如图3b所示为剥离Delaunay三角网的外三角的结果。

点群分布范围属于不确定性的问题。在GIS的研究中，常用作替代点群分布范围。但是，当点群的分布凸凹有致，那么可能在凸壳中包含很多不存在点的凹部区域，这样形成的凸壳并不能真实的代表点群的分布范围。一般地，基于Delaunay三角网的聚类步骤为：①利用点群构建三角网，如图3a；②选择一个距离参数d，通过参数d实现对三角网的“剥皮”，删除三角形边长大于d的边；③将与一个三角形关联的边提取作为点群的边界。图3b显示，通过“剥皮”操作，删除了边长很大的三角形，使得“剥皮”后的三角网更接近真实的点群分布。

步骤6、对步骤5中剥皮后的Delaunay三角网进行AUTOCLUST聚类，得到三角网节点分布密集和分布稀疏的两种类型的若干个聚类结果。AUTOCLUST聚类算法为现有成熟算法，可参阅Estivill-Castro和Lee于2002年发表在《Computers,Environment and urban systems》上的文献“Argument free clustering for large spatial point data sets via boundary extraction from Delaunay Diagram”。如图4a中A、B、C、D、E、F四块区域为高密度聚类所在位置。

除了使用AUTOCLUST聚类方法以外，也可以使用K-means、基于密度的空间聚类等方法对中心点进行聚类，获得分布密集和分布稀疏的两种类型的若干个聚类。

步骤7、利用Vornoio图约束对AUTOCLUST聚类结果进行修正，避免过度聚类和欠聚类，如图4b所示，图中黑色四边形和三角形内代表未被聚类的居民点区域，即欠聚类。本实施例中，利用Delaunay三角网节点(也就是第4步中获得的中心点)构建Voronoi图，计算Voronoi图中Voronoi多边形面积，Voronoi多边形面积大于上限阈值时，初步确定该多边形对应的中心点为耕地，Voronoi多边形面积小于下限阈值时，初步确定该Voronoi多边形对应的中心点属于破碎区域(非耕地)，将破碎区域的中心点删去或归入节点分布密集的聚类中。完成对聚类结果的修正，避免过度聚类和欠聚类现象。

如图4c所示，为使用下限阈值进行约束的结果，图中线条较粗的Voronoi多边形为面积小于下限值的多边形，初步确定为非耕地。如图4d所示，为实用上限阈值进行约束的结果，图中线条较粗的Voronoi多边形为面积大于下限值的多边形，初步确定为耕地。

基于Delaunay三角网构建Voronoi图，计算多边形面积，表示为A。所以，Voronoi多边形的点密度可以通过1/A表征。一般地，该值越小，点群密度越小，该值越大，表示点群密度越大。泰森多边形面积约束策略：一般地，AUTOCLUST进行聚类时，有不少区域出现欠聚类的现象，那么在采用泰森多边形进行约束时，使用上限阈值与下限阈值综合约束，一般地，当泰森多边形面积大于上限阈值时，该多边形对应的中心点一定是耕地，当泰森多边形面积小于下限阈值时，该多边形对应的中心点一定是属于破碎区域。这里的上下限阈值尽量设置松散。用以下公式过滤高密度点区域和低密度点区域：

这里Voronio_mean为生成的泰森多边形的面积的平均值，和为系数，且明显地，这里的系数决定了点的聚类程度，通过实验，为了避免过度聚类和欠聚类，设置作为一般值，对AUTOCLUST聚类结果进行修正。

在AUTOCLUST聚类的过程中，一些聚类主要还包含较大的三角网，且其面积较周围三角网面积还大，这可能是因为极少数的两三个耕地分割对象对应的中心点距离很近，而距离相邻的其它点较远，这种点往往不是居民区域分割对象的中心点。这种聚类中的三角网个数往往少，且三角形面积的平均值比需要识别的破碎区域的三角网面积的平均值大很多。

所以进一步引入聚类三角网面积均值的概念来去除AUTOCLUST聚类过程中这类错误聚类，方程如下：

其中表示第j个聚类的三角形面积的均值，Ti为第i个三角形的面积，nj为第j个聚类的三角形个数，表示所有聚类的三角形面积均值的平均值，表示第j个聚类的三角形面积集合的方差，m表示聚类的个数，i表示第j个聚类中的三角形序号。

如果第j个聚类满足则将该聚类归入分布稀疏的聚类中去。

步骤8、删除节点分布密集的聚类所对应的中心点，至此完成了非耕地区域的中心点的剔除。

步骤9、利用最大方差约束，剔除剩余的零星林地，若所述步骤1分割后得到分割对象的最大光谱方差Max.diff大于预设阈值，则判定该分割对象为林地，剔除该分割对象所对应的中心点，最终保留下来的所有中心点所对应的分割对象为耕地，得到耕地提取结果并进行精度评估。本实施例中，所述预设阈值为2倍最大光谱方差均值与最大光谱方差标准差之和。

本实施例提出的方法虽然能够聚类居民地聚集区域，并将其剔除，但是对于平原地区，还零星存在不少林地，由于林地的纹理特征相差不大，对应的分割对象往往较大，很难用三角网聚类方法将其剔除。而林地的纹理粗糙，方差较大，耕地较为光滑，方差较小，所以这里提出方差约束的方法用以对耕地信息提取结果做修剪。本实施例采用2倍最大光谱方差均值与最大光谱方差标准差之和作为阈值进行识别。如下式：

Max.diff>2×mean_Max.diff+Standard Deviation_Max.diff (9)

其中mean_Max.diff表示分割图层中所有对象的最大光谱方差特征的均值；特征Max.diff的标准差表示为Standard Deviation_Max.diff。

图5a为采用基于三角网的面向对象耕地信息自动提取方法提取的耕地信息，图5a中，白色表示提取的耕地，黑色表示识别的非耕地，图5b表示利用eCognition软件进行监督分类提取的耕地信息，图中白色为耕地，黑色为非耕地。比较发现：本发明提取方法在保证精度要求的前提下，不仅能够在整体上保持耕地的集中连片特性，还能在一定程度上避免由于播种类型差异(物候差异)，导致的零星错分现象。同时，避免监督分类选择训练样本的过程，提高了自动化处理程度。

除上述实施例外，本发明还可以有其他实施方式。凡采用分割对象中心点聚类等形式的耕地信息提取技术方案，均落在本发明要求的保护范围。