CN103559493A - 一种线状地物的提取方法 - Google Patents

Abstract

本发明提供了一种线状地物的提取方法，对遥感影像进行预处理，得到预处理图像；从预处理图像中提取植被区域并将其剔除，得到目标区域；对目标区域的图像根据线状地物类型进行特征提取，并进行矢量化，得到的线状地物矢量化文件。该方法首先对遥感影像进行预处理和植被区域的提取，对线状地物周边的其它地物进行剔除，使得目标区域在影像中突出显示。再采用双阈值滤波对目标区域根据线状地物类型的不同进行特征提取，可以克服单阈值滤波提取时需要选择合适阈值的问题，具有较好的线状地物边缘检测能力，最大程度地避免椒盐现象，便于矢量化，弥补现有技术由于只能定性对线状地物进行检测和提取对耕地质量等级进行评价的不足。

Description

一种线状地物的提取方法

技术领域

本发明涉及遥感影像处理技术领域，尤其涉及一种线状地物的提取方法。

背景技术

随着高空获取遥感影响技术的迅速发展，目前在国土资源调查与监测、基础地理数据更新等方面广泛应用了高空遥感影像数据以及对遥感影像数据的分析和处理。其中，获取的底面遥感影像中灌排水渠、田间道路等耕地线状基础设施是影响农田识别的重要因素之一，因此对于灌排水渠、田间道路等耕地现状地物的正确识别与提取，对农田识别质量评价具有价值和意义。

传统的方式是根据遥感数据图像进行人工解译，但是误差很大，并且效率太低，需要耗费大量的人工劳动力。另外，在对遥感影像数据进行处理的过程中线状地物周边的植被等其它地物会对需要提取的线状地物产生干扰，如田间道路有行道树遮蔽覆盖，给目标区域线状地物的提取带来较大的困难。只能定性的进行地物的提取，还无法满足对土地整治区耕地质量等级进行快速评定的要求。

发明内容

（一）要解决的技术问题

针对上述缺陷，本发明要解决的技术问如何实现在对遥感影像中的线状地物进行提取的过程中尽量避免椒盐现象的发生，从而能更加有利于耕地质量等级的评价。

（二）技术方案

为解决上述问题，本发明提供了一种线状地物的提取方法，包括以下步骤：

S1、对遥感影像进行预处理，得到预处理图像；

S2、从所述预处理图像中提取植被区域并将其剔除，得到目标区域；

S3、对所述目标区域的图像根据线状地物类型进行特征提取，并进行矢量化，得到的线状地物矢量化文件。

进一步地，步骤S1中所述对遥感影像进行预处理具体包括：

S11、对所述遥感影像进行全色影像和多光谱影像的融合，得到融合影像；

S12、对所述融合影像进行灰度值计算，提取线状地物的最大灰度值和最小灰度值；

S13、根据所述最大灰度值和最小灰度值进行灰度二值化计算，得到所述预处理图像。

进一步地，所述融合影像的分辨率与全色波段一致，且所述融合影像的波段数与所述多光谱影像的波段数相同。

进一步地，步骤S2从所述预处理图像中提取植被区域具体包括：

S21、确定所述预处理图像中的植被区域并进行分类，包括：团状植被区域、簇状植被区域和片状植被区域；

S22、对团状植被区域和簇状植被区域采用面积阈值法进行剔除，对片状植被区域采用扫描方法进行剔除，得到所述目标区域。

进一步地，步骤S3对所述目标区域的图像根据线状地物类型进行特征提取具体包括：

S31、对所述目标区域根据不同线状地物的灰度值进行波段运算，得到二值化图；

S32、对所述二值化图进行去噪处理；

S33、进行边缘特征提取，得到边缘图像。

进一步地，步骤S32对所述二值化图进行去噪处理具体包括：

S321、将滤波模板在所述二值化图中进行移动，从所述二值化图中选取一个像素作为中心像素，并将所述滤波模板的中心与所述中心像素的位置重合；

S322、读取所述滤波模板中的所有像素的灰度值，按照所述灰度值的等级进行排序；

S323、从中排序的灰度值中选择平均值作为当前灰度值进行中值滤波，去出孤立噪声。

进一步地，步骤S33所述对所述目标区域的图像根据线状地物类型进行特征提取时采用自适应双阈值滤波，具体包括：

S331、对所述目标区域经过二值化和去噪处理得到的图像进行高斯滤波，通过利用两个一维高斯核两次加权实现或者利用一个二维高斯核和一次卷积实现，其中一维高斯核为：

$K_1=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}{e}^{-\frac{x*x}{2\mathrm{\σ}*\mathrm{\σ}}},$

二维高斯核为： $K_2=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^{*}\mathrm{\σ}}{e}^{-\frac{x*x+y*y}{2\mathrm{\σ}*\mathrm{\σ}}},$

卷积算子为： $s_x=\left[\begin{array}{cc}-1& 1\\ -1& 1\end{array}\right],s_y=\left[\begin{array}{cc}1& 1\\ -1& -1\end{array}\right],$ s_x为卷积算子沿x向的一阶偏导数矩阵，s_y为卷积算子沿y向的一阶偏导数矩阵；

S332、计算梯度幅值P[i,j]、Q[i,j]、M[i,j]以及梯度方向θ[i,j]的计算公式为：

P[i,j]=(f[i,j+1]-f[i,j]+f[i+1,j+1]-f[i+1,j])/2；

Q[i,j]=(f[i,j]-f[i+1,j]+f[i,j+1]-f[i+1,j+1])/2；

$M[i,j]=\sqrt{{P[i,j]}^2+{Q[i,j]}^2};$

θ[i,j]=arctan(Q[i,j]/P[i,j])；

S333、计算梯度幅值对所述梯度幅值经非极大值抑制，剔除所述灰度图中的非边缘点；

S334、再利用双阈值法检测边缘，并进行边缘连接，得到闭合的边缘图像。

进一步地，S3中进行矢量化时通过进行栅格图像矢量化操作得到的，具体包括：

S34、对所述边缘图像进行细化，得到栅格数据；

S35、跟踪所述栅格数据，并整理成为线条，得到弧段和多边形；

S36、去除所述弧段和所述多边形上的多余数据点；

S37、判断所述弧段和所述多边形之间的空间关系，形成线状地物的矢量化文件。

进一步地，得到所述线状地物的矢量化文件之后进一步还包括：根据不同线状地物的矢量化文件对耕地质量等级进行评价，对于所述矢量化文件表示的线状地物的基本属性进行提取，并根据耕地等级检测评定指标进行耕地质量等级的评价。

进一步地，进行双阈值滤波时，根据需要提取线状地物的光谱特征选取最大阈值和最小阈值，用所述最小阈值的图像去匹配所述最大阈值的图像，对所述最高阈值的图像中缺失的线状地物进行补充，使得提取的线状地物具有完整的闭合曲线。

（三）有益效果

本发明提供了一种线状地物的提取方法，对遥感影像进行预处理，得到预处理图像；从预处理图像中提取植被区域并将其剔除，得到目标区域；对目标区域的图像根据线状地物类型进行特征提取，并进行矢量化，得到的线状地物矢量化文件。该方法首先对遥感影像进行预处理和植被区域的提取，对线状地物周边的其它地物进行剔除，使得目标区域在影像中突出显示。再采用双阈值滤波对目标区域根据线状地物类型的不同进行特征提取，可以克服单阈值滤波提取时需要选择合适阈值的问题，具有较好的线状地物边缘检测能力，最大程度地避免椒盐现象，便于矢量化，弥补现有技术由于只能定性对线状地物进行检测和提取对耕地质量等级进行评价的不足。

附图说明

图1为本发明实施例提供的一种线状地物的提取方法的步骤流程图；

图2为本发明实施例提供的步骤S1中对遥感影像进行预处理的步骤流程图；

图3为本发明实施例提供的步骤S2进行植被区域的提取的步骤流程图；

图4为本发明实施例中经过植被区域剔除得到的图像；

图5为本发明实施例提供的步骤S3对目标区域的图像根据线状地物类型进行特征提取的步骤流程；

图6为本发明实施例中经过灰度二值化得到的二值化图；

图7为本发明实施例提供的步骤S32中进行去噪处理的步骤流程图；

图8为本发明实施例提供的对经过二值化和去噪处理后的目标线性地物进行边缘特征提取的步骤流程图；

图9为本发明实施例中非极大值抑制的原理示意图；

图10为本发明实施例中经过边缘特征提取之后得到的边缘图像；

图11为本发明实施例提供的边缘提取之后的矢量化操作步骤流程图；

图12为本发明实施例中删除多余点的算法处理示意图；

图13为本发明实施例中最后得到的矢量化图像。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

面对现在各级监管部门实时监管的要求，传统的人工解译难以满足大面积线状地物的识别提取以及耕地质量等级评价应用的需求。利用航空遥感和卫星航天遥感可以在短期内获得大范围的高空间分辨率遥感影像，可以清晰辨识田间道路、灌排沟渠等线状地物，这为耕地质量部分指标的计算提供了基础。

本发明实施例提供了一种线状地物的提取方法，步骤流程如图1所示，具体包括以下步骤：

步骤S1、对遥感影像进行预处理，得到预处理图像。

步骤S2、从预处理图像中提取植被区域并将其剔除，得到目标区域。

步骤S3、对目标区域的图像根据线状地物类型进行特征提取，并进行矢量化，得到的线状地物矢量化文件。

针对目前耕地、灌排沟渠等线状地物在耕地质量等级评价中存在的问题，本发明实施例提供了一种线状地物的提取方法，通过对遥感影像进行预处理和植被区域的剔除，使得目标区域在影像中突出显示。再采用双阈值滤波对目标区域根据线状地物类型的不同进行特征提取，可以克服单阈值滤波提取时需要选择合适阈值的问题，具有较好的线状地物边缘检测能力，最大程度地避免椒盐现象，便于矢量化。

优选地，上述步骤S1中对遥感影像进行预处理的步骤流程如图2所示，具体包括以下步骤：

步骤S11、对遥感影像进行全色影像和多光谱影像的融合，得到融合影像。选择需要研究和进行地物提取的高分辨率遥感影像数据，该遥感影像数据一般是航空或航天多光谱影像数据，能够覆盖农田线状地物及周边地物，例如田间道路和灌溉设施等。同时根据耕地线状地物以及周边还有行道树等其它地物生长。除了多光谱影像还需要全色影像（即全色波段数据）进行融合，得到的融合影像的分辨率与全色波段一致，且融合影像的波段数与多光谱影像的波段数相同。

步骤S12、对融合影像进行灰度值计算，提取线状地物的最大灰度值和最小灰度值。

步骤S13、根据最大灰度值和最小灰度值进行灰度二值化计算，得到预处理图像。根据步骤S12得到的地物最大灰度值和最小灰度值，编写波段运算语句，进行遥感影像数据的灰度二值化，使得能够在影响上突出想要研究的目标地物。

对遥感影像数据的预处理还包括：精确集合配准、辐射校正以及影像降噪处理，以便于遥感影像中的线状地物的提取。

之后步骤S2中对步骤S1得到的预处理图像采用NDVI（Normalized Difference Vegetation Index，归一化植被指数，也称为标准差异植被指数）进行植被区域的提取，以得到目标区域，步骤流程如图3所示，具体包括以下步骤：

步骤S21、确定预处理图像中的植被区域并进行分类，包括：团状植被区域、簇状植被区域和片状植被区域。根据植被区域与其他地物的不同特性进行确定。

步骤S22、对团状植被区域和簇状植被区域采用面积阈值法进行剔除，对片状植被区域采用扫描方法进行剔除，得到目标区域，以一幅覆盖有耕地、田间道路、行道树等其他周边地物的遥感影像经过植被区域剔除得到的图像如图4所示。

上述步骤对预处理图像中的制备区域进行初步确定，对制备区域的提取主要是对行道树和防护林等植被区域进行剔除。根据耕地线状地物和周边其它地物的实际情况，田间道路两边会有行道树的遮蔽和覆盖，灌溉沟渠周边也会有植被覆盖。对于植被区域根据不同的分类采用不同的方法进行处理能够更加具有针对性地对植被区域进行剔除，例如灌溉沟渠中生长的植被或者田间孤立的植被一般是团状或者簇状呈小块区域分布的植被，采用面积与执法进行剔除，而对于成片的行道树等呈片状分布的植被需要采用扫描法进行剔除。通过采用分类的方式进行植被区域的剔除，能够有效地解决行道树对田间道路等线状地物的提取产生较大干扰的区域，从而在得到的影像中更加突出显示线状地物，即目标地物。

进一步地，步骤S3对目标区域的图像根据线状地物类型进行特征提取的步骤流程如图5所示，具体包括以下步骤：

步骤S31、对目标区域根据不同线状地物的灰度值进行波段运算，得到二值化图，如图6所示。对目标区域的图像进行二值化，突出所要提取的田间道路、灌溉沟渠等线状地物，并根据线状地物的不同类型，选取不同的灰度值进行波段运算，得到包含田间道路、灌溉沟渠等线状地物的二值化图。

步骤S32、对二值化图进行去噪处理。采用中值滤波对上述步骤得到的二值化图进行去噪处理，能够有效去除二值化图像中的椒盐现象，为下一步边缘提取提供良好的精确度。

中值滤波是一种非线性的信号处理方法，在一定条件下可以克服线性滤波器处理图像细节模糊的问题，并且对于滤除脉冲干扰和图像扫面噪声非常有效，中值滤波的基本思想就是把局部区域的像素按灰度等级进行排序，取该领域中灰度中值作为当前像素的灰度值，通过中值滤波去除孤立的高亮噪点，能够在保证二值化图像地物基础特性不丢失的情况下很好地去除椒盐现象。另外，采用中值滤波将二值化图像中的局部极值消除掉，不会产生新的区域极值，处理后的重建图像中保存下来的区域轮廓就不会发生位置偏移。

上述步骤S32对所述二值化图进行去噪处理的步骤流程如图7所示，具体包括以下步骤：

步骤S321、将滤波模板在二值化图中进行移动，从二值化图中选取一个像素作为中心像素，并将滤波模板的中心与中心像素的位置重合。其中的滤波模板是含有若干个点的滑动窗口。

步骤S322、读取滤波模板中的所有像素的灰度值，按照灰度值的等级进行排序。排序的时候可以按照从小到大的顺序进行排列，还可以是按照从大到小的顺序进行排列。

步骤S323、从中排序的灰度值中选择平均值作为当前灰度值进行中值滤波，去出孤立噪声。在选取平均值（即中值）时，根据上述步骤排序结果选取序列的中间数据将其赋给滤波模板中心位置的像素，如果窗口中有奇数个元素，则按灰度值的大小排序后的中间元素灰度值作为中值。如果窗口中有偶数个元素，则按灰度值大小排序后选取中间两个元素灰度的平均值作为中值。

通过使用中值滤波的方法进行去噪处理，对于保护边缘的特性很有帮助，尽量减少出现边缘模数的现象。二值化图像为二维信号，中值滤波的窗口形状和尺寸对滤波器效果影响很大，不同图像内容和不同应用要求往往选用不同的窗口形状和尺寸。本发明实施例中优先选取5*5的窗口进行滤波。

步骤S33、进行边缘特征提取，得到边缘图像。此步骤中采用只适应Canny算子进行滤波，能够对二值化和去噪后的目标线性地物的边缘进行提取，步骤流程如图8所示，具体包括以下步骤：

步骤S331、对所述目标区域经过二值化和去噪处理得到的图像采用Canny算法进行高斯滤波。

Canny算法通常处理的图像为灰度图，因此如果摄像机获取的是彩色图像，就需要首先进行灰度化。对一幅彩色图进行灰度化就是根据图像各个通道的采样值进行加权平均。由于本发明在步骤S1中已经进行过图像的灰度化处理，因此这一步已经完成。直接进行高斯滤波，利用两个一维高斯核两次加权实现或者利用一个二维高斯核和一次卷积实现，其中一维高斯核为：

$K_1=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}{e}^{-\frac{x*x}{2\mathrm{\σ}*\mathrm{\σ}}}$                     公式（1）

上式为离散化的一维高斯函数，确定参数就可以得到一维高斯核向量。

二维高斯核为：

$K_2=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^{*}\mathrm{\σ}}{e}^{-\frac{x*x+y*y}{2\mathrm{\σ}*\mathrm{\σ}}}$                        公式（2）

上式为离散化的二维高斯函数，确定参数就可以得到二维高斯核向量。

卷积算子为： $s_x=\left[\begin{array}{cc}-1& 1\\ -1& 1\end{array}\right],s_y=\left[\begin{array}{cc}1& 1\\ -1& -1\end{array}\right],$ s_x为卷积算子沿x向的一阶偏导数矩阵，s_y为卷积算子沿y向的一阶偏导数矩阵。

步骤S332、计算梯度幅值P[i,j]、Q[i,j]、M[i,j]以及梯度方向θ[i,j]的计算公式为：

P[i,j]=(f[i,j+1]-f[i,j]+f[i+1,j+1]-f[i+1,j])/2        公式（3）

Q[i,j]=(f[i,j]-f[i+1,j]+f[i,j+1]-f[i+1,j+1])/2         公式（4）

$M[i,j]=\sqrt{{P[i,j]}^2+{Q[i,j]}^2}$                         公式（5）

θ[i,j]=arctan(Q[i,j]/P[i,j])                        公式（6）

其中，i、j分别表示所在梯度的像素点坐标值。

步骤S333、计算梯度幅值对梯度幅值经非极大值抑制，剔除灰度图中的非边缘点。在Canny算法中，非极大值抑制是进行边缘检测的重要步骤，一般是指寻找像素点局部最大值，将非极大值点所对应的灰度值置为0，这样可以剔除掉一大部分非边缘的点。

具体的非极大值抑制的工作原理为：

以图9为例，要进行非极大值抑制，就首先要确定像素点C的灰度值在其8值邻域内是否为最大。图9中线条L1的方向为C点的梯度方向，θ表示上述步骤得到的梯度方向，g1、g2、g3、g4分别代表像素点八邻域的其中四个方向点。这样就可以确定其局部的最大值肯定分布在这条线上，也即出了C点外，梯度方向的交点dTmp1和dTmp2这两个点的值也可能会是局部最大值。因此，判断C点灰度与这两个点灰度大小即可判断C点是否为其邻域内的局部最大灰度点。如果经过判断，C点灰度值小于这两个点中的任一个，那就说明C点不是局部极大值，那么则可以排除C点为边缘。

步骤S334、再利用双阈值法检测边缘，并进行边缘连接，得到闭合的边缘图像。

进行双阈值滤波时，根据需要提取线状地物的光谱特征选取最大阈值和最小阈值，用最小阈值的图像去匹配最大阈值的图像，对最高阈值的图像中缺失的线状地物进行补充，使得提取的线状地物具有完整的闭合曲线。

Canny算法中减少假边缘数量的方法是采用双阈值法进行边缘检测。选择两个阈值，根据高阈值（即值较大的阈值）得到一个边缘图像，这样一个图像含有很少的假边缘，但是由于阈值较大，产生的图像边缘可能不闭合。为了解决这样一个问题采用了另外一个低阈值（即值较小的阈值）。在高阈值图像中把边缘链接成轮廓，当到达轮廓的端点时，该算法会在断点的8邻域点中寻找满足低阈值的点，再根据此点收集新的边缘，直到整个图像边缘闭合，至此得到闭合的边缘图像。进行边缘提取之后得到的图如图10所示。

基于Canny算发得到的线状地物进行边缘提取后，图像中地物仍然仅停留在定性分析阶段，因此，还需要对进行边缘提取的地物轮廓进行矢量提取和输出，因此步骤S3还包括矢量化操作阶段，进行矢量化时通过进行栅格图像矢量化操作得到的，步骤流程如图11所示，具体包括：

步骤S34、对边缘图像进行细化，得到栅格数据。

传统的矢量化工具主要有ArcMap中的矢量化工具，而如果对于上述步骤中提取的轮廓再进行人工矢量化，就有可能造成二次误差，影响后期耕地质量等级的评定，因此，本发明采用如下方法进行栅格图像矢量化操作：

1）二值化：一般情况下，栅格数据是按0～255的不同灰度值表达的。为了简化追踪算法，需把256个灰阶压缩为2个灰阶，即0和1两级。本发明实施例中的图像已经在之前的步骤中进行了二值化操作，因此可以直接进行下一步。

2）细化：细化是消除线划横断面栅格数的差异，使得每一条线只保留代表其轴线或周围轮廓线位置的单个栅格的宽度。对于栅格线划的细化方法，本发明实施例中采用"剥皮法"。剥皮法的实质是剥掉等于一个栅格宽的一层，直到最后留下彼此连通的由单个栅格点组成的图形。

步骤S35、跟踪栅格数据，并整理成为线条，得到弧段和多边形。

跟踪的目的是把细化后的栅格数据整理为从结点出发的线段或闭合的线条，并以矢量形式加以存储。跟踪时，从起始点开始，在保证趋势的情况下对八个邻域进行搜索，依次得到相邻点，最终得到完整的弧段或多边形。

步骤S36、去除弧段和多边形上的多余数据点。

由于搜索是逐个栅格进行的，所以，弧段或多边形的数据列十分密集。为了减少存储量，在保证线段精度的情况下可以删除部分多余的数据点。

本发明实施例采用如下算法删除多余点：计算当前点A与相邻点B、C组成的线段BA、AC的夹角，如果夹角大于某一固定值（如150度），就删除B点。当然，如果删除B点后，AC距离太长，那么，删除B点是不合适的，需要保留B点，上述算法的处理示意图如图12所示。

另外，去除多余数据点之后还包括曲线光滑处理，采用曲线张力样条函数，张力样条函数是Schweikert为消除三次样条插值函数有时会出现多余的拐点而引入的。其基本构思是分段插值函数为直线插值和两个双曲函数shξx和chξx的线性组合：

f(x)=C1+C2X+C3shξx+C4chξx                公式（7）

其中公式（7）中的ξ是张力系数，作用是控制拐点的位置和曲线的形态。C1和C2是直线插值的相关系数，C3和C4是曲线插值的相关系数，通过有关复杂的数学推导，可以得到张力样条函数的函数表达式：

$\begin{array}{c}f\left(x\right)=\frac{1}{{\mathrm{\σ}}^2\mathrm{sh}\left({\mathrm{\σh}}_i\right)}[f^{\′\′}\left(x_i\right)\mathrm{sh}\left(\mathrm{\σ}(x_{i+1}-x)\right)+f^{\′\′}\left(x_{i+1}\right)\mathrm{sh}\left(\mathrm{\σ}(x-x_i)\right)]+\\ [y_i-\frac{f^{\′\′}\left(x_i\right)}{{\mathrm{\σ}}^2}]\frac{x_{i+1}-x}{h_i}+[y_{i+1}-\frac{f^{\′\′}\left(x_{i+1}\right)}{{\mathrm{\σ}}^2}]\frac{x-x_i}{h_i}\end{array}$     公式（8）

其中x_i≤x≤x_i+1，其中i=0,1,2,...n-1，h_i=x_i+1-x_i。

从上述公式（8）中张力样条函数优良的数学特性可以看出，只要选择合适的张力系数，就可以尽量避免合乎要求的光滑曲线相交。

步骤S37、判断弧段和多边形之间的空间关系，形成线状地物的矢量化文件。判断弧段与多边形间的空间关系，以形成完整的拓扑结构并建立与属性数据的关系，得到的矢量化图像如图13所示。

优选地，得到线状地物的矢量化文件之后进一步还包括：根据不同线状地物的矢量化文件对耕地质量等级进行评价，对于矢量化文件表示的线状地物的基本属性进行提取，并根据耕地等级检测评定指标进行耕地质量等级的评价。对于矢量化文件通过测定线状地物的长度、宽度、周长和面积等基本属性，引入耕地等级检测评定相关指标，进行耕地等级评价，如道路通达度和灌溉保证率。

通过以上步骤，获得了线状地物的基本矢量文件，下一步将进行耕地质量等级评价。由于耕地的好坏评价指标除了通过耕地本身的亩产外，还需要其周边线状地物如灌溉沟渠、田间道路建立相关评价指标进行分等定级的，如果仅以边缘提取后的图像观察线状地物分布进行定性评价，将会造成很大的误差；另一方面，灌溉沟渠和田间道路等地物在图像上提取时其光谱特征并不相同，如果一次提取多个类别的地物，容易造成错分和混分，对于后期的质量等级评定影像较大。基于以上两点考虑，本发明实施例中以不同地物建立不同特征提取方法为目标，通过逐个地物的提取以及矢量化，最后统一进行耕地质量评价。

针对不同的线状地物，如何选取合适的特征阈值进行提取是进行耕地质量评定的关键。本发明实施例提供的方法，根据不同田间线状地物的不同特征，以及植被覆盖存在干扰性的特点，前期采用面积阈值法和扫描法进行植被区域的剔除，然后采用灰度选择后进行二值化的方法，在图上充分突出所要提取的地物，辅以Canny算子进行边缘轮廓的提取，最后进行矢量输出后进行耕地质量的相关评价。基于土地整治区高分辨遥感影像，充分利用线状地物的形状特征、光谱特征以及最大灰度值、最小灰度值，对后期研究土地整治区新增田间道路、灌溉沟渠、排水渠和农田林网等耕地基础设施信息提供了一种自动获取线状地物进行评价的方法。

综上所述，通过本发明实施例提供的现状地物的提取方法，首先对线状地物以外的其他地物进行先行剔除的办法，最大程度上减少提取时的干扰，同时，易于快速识别和提取，提高算法实现效率，通过选取不同线状地物的灰度值范围，利用边缘提取算子（Canny算子）进行提取，通过对提取后的线状地物矢量化，得到具有长度、宽度、周长和面积的线状地物基本属性，通过建立相应缓冲区引入耕地质量评价指标，如道路通达度和灌溉保证率进行耕地质量评定，实现了在土地整治区土地质量等级评价中线状地物部分的自动识别提取与应用，有助于对整治区耕地质量等级进行快速评定。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (10)

1.一种线状地物的提取方法，其特征在于，所述方法包括：

S1、对遥感影像进行预处理，得到预处理图像；

S2、从所述预处理图像中提取植被区域并将其剔除，得到目标区域；

S3、采用双阈值滤波对所述目标区域的图像根据线状地物类型进行特征提取，并进行矢量化，得到的线状地物矢量化文件。

2.如权利要求1所述的提取方法，其特征在于，步骤S1中所述对遥感影像进行预处理具体包括：

S11、对所述遥感影像进行全色影像和多光谱影像的融合，得到融合影像；

S12、对所述融合影像进行灰度值计算，提取线状地物的最大灰度值和最小灰度值；

S13、根据所述最大灰度值和最小灰度值进行灰度二值化计算，得到所述预处理图像。

3.如权利要求2所述的提取方法，其特征在于，所述融合影像的分辨率与全色波段一致，且所述融合影像的波段数与所述多光谱影像的波段数相同。

4.如权利要求1所述的方法，其特征在于，步骤S2从所述预处理图像中提取植被区域具体包括：

S21、确定所述预处理图像中的植被区域并进行分类，包括：团状植被区域、簇状植被区域和片状植被区域；

S22、对团状植被区域和簇状植被区域采用面积阈值法进行剔除，对片状植被区域采用扫描方法进行剔除，得到所述目标区域。

5.如权利要求1所述的提取方法，其特征在于，步骤S3对所述目标区域的图像根据线状地物类型进行特征提取具体包括：

S31、对所述目标区域根据不同线状地物的灰度值进行波段运算，得到二值化图；

S32、对所述二值化图进行去噪处理；

S33、进行边缘特征提取，得到边缘图像。

6.如权利要求5所述的提取方法，其特征在于，步骤S32对所述二值化图进行去噪处理具体包括：

S321、将滤波模板在所述二值化图中进行移动，从所述二值化图中选取一个像素作为中心像素，并将所述滤波模板的中心与所述中心像素的位置重合；

S322、读取所述滤波模板中的所有像素的灰度值，按照所述灰度值的等级进行排序；

S323、从中排序的灰度值中选择平均值作为当前灰度值进行中值滤波，去出孤立噪声。

7.如权利要求5所述的提取方法，其特征在于，步骤S33所述对所述目标区域的图像根据线状地物类型进行特征提取时采用自适应双阈值滤波，具体包括：

S331、对所述目标区域经过二值化和去噪处理得到的图像进行高斯滤波，通过利用两个一维高斯核两次加权实现或者利用一个二维高斯核和一次卷积实现，其中一维高斯核为：

$K_1=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}{e}^{-\frac{x*x}{2\mathrm{\σ}*\mathrm{\σ}}},$

二维高斯核为： $K_2=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^{*}\mathrm{\σ}}{e}^{-\frac{x*x+y*y}{2\mathrm{\σ}*\mathrm{\σ}}},$

卷积算子为： $s_x=\left[\begin{array}{cc}-1& 1\\ -1& 1\end{array}\right],s_y=\left[\begin{array}{cc}1& 1\\ -1& -1\end{array}\right],$ s_x为卷积算子沿x向的一阶偏导数矩阵，s_y为卷积算子沿y向的一阶偏导数矩阵；

S332、计算梯度幅值P[i,j]、Q[i,j]、M[i,j]以及梯度方向θ[i,j]的计算公式为：

P[i,j]=(f[i,j+1]-f[i,j]+f[i+1,j+1]-f[i+1,j])/2；

Q[i,j]=(f[i,j]-f[i+1,j]+f[i,j+1]-f[i+1,j+1])/2；

$M[i,j]=\sqrt{{P[i,j]}^2+{Q[i,j]}^2};$

θ[i,j]=arctan(Q[i,j]/P[i,j])；

S333、计算梯度幅值对所述梯度幅值经非极大值抑制，剔除所述灰度图中的非边缘点；

S334、再利用双阈值滤波检测边缘，并进行边缘连接，得到闭合的边缘图像。

8.如权利要求1所述的提取方法，其特征在于，S3中进行矢量化时通过进行栅格图像矢量化操作得到的，具体包括：

S34、对所述边缘图像进行细化，得到栅格数据；

S35、跟踪所述栅格数据，并整理成为线条，得到弧段和多边形；

S36、去除所述弧段和所述多边形上的多余数据点；

S37、判断所述弧段和所述多边形之间的空间关系，形成线状地物的矢量化文件。

9.如权利要求1所述的提取方法，其特征在于，得到所述线状地物的矢量化文件之后进一步还包括：根据不同线状地物的矢量化文件对耕地质量等级进行评价，对于所述矢量化文件表示的线状地物的基本属性进行提取，并根据耕地等级检测评定指标进行耕地质量等级的评价。

10.如权利要求1或7所述的提取方法，其特征在于，进行双阈值滤波时，根据需要提取线状地物的光谱特征选取最大阈值和最小阈值，用所述最小阈值的图像去匹配所述最大阈值的图像，对所述最高阈值的图像中缺失的线状地物进行补充，使得提取的线状地物具有完整的闭合曲线。

Images