CN113063375A - 一种直线型耕作田埂的无人机遥感提取方法 - Google Patents

Abstract

本发明提供的一种直线型耕作田埂的无人机遥感提取方法，包括：1)先采集单时相耕地无人机序列图像，得到其正射影像和数字表面模型后，利用数字表面模型计算地表粗糙度，通过阈值分割将其二值化、形状指数滤波候选田埂线，再以正射影像进行植被分割的二值化图为掩膜剔除植被的影响，然后形态学运算获取耕作田埂线栅格二值图；2)对所提取的耕作田埂利用霍夫变换检测其上的点，再对这些点实施标记、排序、去冗余和修缮步骤，然后将点按照所述田埂连接为多段线(即耕地田埂线)。本发明步骤简洁明了，自动化程度较高，有利于快速提取耕地的耕作田埂。

Description

一种直线型耕作田埂的无人机遥感提取方法

技术领域

本发明属于土地信息技术、农业和遥感技术领域，特别涉及一种直线型耕作田埂的无人机遥感提取方法。

背景技术

耕地是人类赖以生存的基本资源，特别是立地条件良好的耕地所生产粮食占有极大比例。农业田间管理以耕作单元为依据，而当前全球以小斑块农田为基本单元的农业景观普遍存在，如在中国华北平原地区表现为耕地条带田块现象；精细的田间管理与决策需要具体田块信息的支持，耕作田埂是基本耕作单元的管理界线，当前的提取通常以实地测量与人工数字化为主，作业效率较低，而利用高分辨率卫星影像所获取的精度与预期使用精度存在差距。因此为精细监测与评价以耕地田块为基本单元，耕作田埂线的快速提取是极其必要的。

近几年来，轻小型无人机摄影测量技术以响应迅速、操作简便、重访周期短、成图精度高等诸多优势，得到了农业、林业、城市、矿业等诸多领域的广泛青睐。直线型耕作田埂在我国耕地景观较为常见，比如华北黄淮海平原耕作区、西南山地丘陵土地整治后的耕地区，是基本耕作管理单元的分隔界线。以华北平原为例，该地区耕地面积所占比例超80％，小麦和玉米产量分别占我国总产量的75％、32％，具有典型的耕地条带田块景观，是直线型田埂的典型地区之一。典型的耕地条带田块一般宽3～8米、长50～300米，由两个相邻的田埂线作为边界线；田块之间的田埂一般宽30～40cm、高10～20cm，是两个相邻田块的分界线，传统测量技术或高分辨率卫星影像难以实现高效、精确的测图。然而，轻小型无人机摄影测量技术获取厘米级的影像，能够清晰识别细长的田埂线，具有显著的图像特征。田埂与相邻田面存在一定的高差，在数字高程模型的图像上表现为显著的灰度值差异；且狭长的形态特征与道路的图像特征十分相似，因此借鉴道路检测方法的思路、挖掘厘米级无人机影像信息并实现直线型耕作田埂的快速提取具有挑战性。

发明内容

本发明的目的是为解决上述难题，提出一种直线型耕作田埂的无人机遥感提取方法，本发明利用单时相厘米级无人机影像的数字表面模型和正射影像，旨在将数字图像处理技术和地图制图学的方法结合，达到快速非人工地提取耕地基本耕作单元，实现较高自动化程度的直线型耕作田埂提取，有利于精细的田间管理与决策。

本发明提出一种直线型耕作田埂的无人机遥感提取方法，其特征在于，具体包括以下步骤：

1)在采集单时相耕地的无人机序列图像，并三维重建获取正射影像和数字表面模型基础上，利用数字表面模型计算地表粗糙度，进行阈值分割将其二值化为田埂、非田埂两类对象；再实施形状指数滤波，得到更纯净的田埂二值图；然后以正射影像进行植被分割的二值化图为掩膜剔除植被的影响，修正因植被覆盖影响的边缘突出的田埂二值图；接着通过形态学运算，获取图像边缘更加平滑的田埂二值图；

2)对步骤1)所提取的经平滑处理后的耕地田埂二值图利用霍夫变换检测其上的点，再对这些点标记、排序，并根据几何关系去冗余和修缮，然后将点按照所述田埂连接为多段线(即为耕地田埂线)。

进一步地，所述步骤1)的具体过程如下：

11)采集耕地区域的无人机序列图像，地面同时布设和量测控制点，通过摄影测量软件处理得到数字表面模型和正射影像，再将目标耕地区域的影像按照典型区域的边界剪裁，然后将正射影像的色差空间从RGB模式变换为HSV(Hue-Saturation-Value)模式，对数字表面模型计算得到地表粗糙度f；

式中，地表粗糙度f的值在0～1之间，H为滑动分析窗口的中心像素值，

为滑动分析窗口内的像素均值，n为滑动分析窗口内的像素个数。

12)利用以下公式进行阈值分割，从步骤11)计算的地表粗糙度图像提取耕地田埂二值图f₁：

f(x,y)≥T₁ (2)

其中，地表粗糙度图像中满足公式(2)的像素为耕地田埂，不满足公式(2)的为非耕地田埂，f(x,y)是像素(x,y)处的二进制值，T₁是所采用的分割阈值。

13)利用以下公式进行4个形状指数滤波，对步骤12)所得的候选田埂线进行提纯：

其中，S₀代表地表粗糙度阈值分割后的二值图像，S₁代表第一次过滤后的二值图，S₂代表第二次过滤后的二值图，S₃代表第三次过滤后的二值图，S₄代表第四次过滤后的二值图。

其中，f₂是经过四个形状指数过滤后的候选田埂线图像。

14)将步骤11)所得的HSV模式下正射影像选用色调图像，再对经高斯曲线拟合的色调直方图进行阈值检测，并利用以下公式进行二值化分割，即可得到植被掩膜f₃：

f(x,y)≥T₂ (5)

其中，f(x,y)是像素(x,y)处的色调值，T₂为所检测到的阈值。

15)根据步骤14)所得的植被掩膜，利用以下公式对步骤13)滤波结果的候选田埂线进行图像点乘运算，得到剔除植被影响的候选田埂线f₄：

f₄＝f₂·f₃ (6)

其中，f₂是步骤13)形状指数滤波的结果，f₃是步骤14)植被掩膜的结果，f₄是本步骤剔除植被影响的候选田埂线结果。

16)依次利用基于多方向结构元素的图像开运算、小像素斑块移除、基于线性结构元素的图像闭运算和图像细化对步骤15)所得的候选田埂线进行平滑，得到候选田埂线f₅：

利用以下公式构造多方向性结构元素(g₁)：

其中，g₁(x_i,y_i)是像素(x_i,y_i)在g₁处的像素值，α_i是第i个方向的角度值，间隔长度i的取值范围从-90°至90°，L代表结构元素分析窗口的长度。

利用以下公式进行图像开运算：

其中，g₁代表多方向结构元素，f₄是步骤15)剔除植被影响的结果，f₅是本步骤图像开运算的候选田埂线。

利用以下公式对小像素斑块通过像素对象的面积实现过滤：

其中，T₃是像素斑块面积的过滤阈值，f₅是前一公式剔除植被影响的结果，f₆是本公式小斑块过滤的候选田埂线。

利用以下公式进行图像闭运算：

其中，g₂代表线性结构元素，f₆是前一公式小斑块过滤后的结果，f₇是本公式图像闭运算的候选田埂线。

然后利用图像细化算法，生成田埂线的骨架，即为经过形态学运算平滑的耕地田埂线图像(f₈)。

进一步地，所述步骤2)的具体过程如下：

21)利用霍夫变换算法对步骤16)所得的田埂线二值图进行检测，需要霍夫变换的参数空间定义、参数空间的峰值投票与识别和霍夫变换空间的交点参数确定三个步骤来提取直线段的两个端点，从而获取田埂线上的若干个点，得到图像坐标系下的点坐标，再根据图像的投影坐标信息，利用以下公式将其坐标由图像坐标系转为投影坐标系：

其中，n_x和n_y代表图像中的像素坐标，r₀代表图像分辨率，x和y代表对应的投影坐标，以及x₀和y₀代表图像左上角拐点的投影坐标。

22)根据步骤21)所检测的点进行标记，具体过程为：利用以下公式对步骤21)所检测的点进行旋转变换为正南正北方向的田埂线，再根据这些点与每个田埂最小外接矩形的位置关系判断所述田埂，并按照自西向东方向标记点所在的田埂线序号属性；然后按照自南向北方向标记落在每条田埂线上点的序号：

其中，x和y为变换之前的坐标，x′和y′为变换之后的坐标，θ为旋转角度。

23)根据几何关系对步骤22)所标记和排序的点进行冗余点去除和缺陷点完善，具体过程为：对每个田埂线上的点计算相邻两点间距，若小于所设定的间距阈值，则这两点坐标被两者坐标均值替代，如此反复迭代直到所有相邻两点间距大于上述阈值为止；在坐标旋转变换为正南正北情况下进行三次检测与完善，包括：其一，如果田埂线的最小外接矩形质心与该田埂线一侧的端点距离更近，与该田埂线的最小外接矩形同侧次轴中心点更远，则该田埂线此侧端点坐标被最小外接矩形的次轴中心坐标替换；其二，拐角点坐标采用对应一侧所有田埂线端点的中值点坐标替换；其三，检测为异常值的端点坐标用该点对应的田埂线与同侧相邻两个田埂线端点连线的交点坐标替换；

24)根据地图制图学原理，对步骤23)所得田埂线上的点按照标记所在的田埂编号和属于每个田埂的点集序号依次连接为矢量田埂线，从而实现耕地耕作田埂的无人机遥感提取。

本发明的特点及有益效果：

本发明利用轻小型无人机摄影测量所得的数字地表模型计算地表粗糙度，实现对耕地区域田埂线栅格的提取，并采用数字图像处理技术和地图制图学原理的方法实现耕地田埂线矢量的提取。该发明操作简单，制图效率高，极大地实现了单次无人机摄影测量产品支持下的耕地耕作田埂自动化制图。

通过对数字表面模型计算的地表粗糙度阈值分割法对研究区田埂的提取，再通过形状指数滤波、植被影响的剔除和形态学运算平滑获取候选田埂二值图，然后通过霍夫变换检测田埂线上的点，并进行标记排序、去冗余和完善，最后按照地图制图学原理将点连接为矢量田埂线。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的其中五幅，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例得到的局部地表粗糙度阈值分割前图示；

图2为本发明实施例得到的局部地表粗糙度阈值分割后图示；

图3中为本发明实施例得到的耕地田埂线霍夫变换检测点结果、所检测点的标记的示意图；

图4为本发明实施例得到的耕地田埂线霍夫变换排序示意图；

图5为本发明实施例得到的耕地田埂线提取结果与参考值叠加图结果.

具体实施方式

本发明提出的一种直线型耕作田埂的无人机遥感提取方法，结合附图及实施例详细说明如下：

本发明提出一种直线型耕作田埂的无人机遥感提取方法，是基于数字图像处理和地图学原理，具体包括以下步骤：

1)选择作物适宜生长季实施单时相耕地的无人机序列图像采集，在三维重建获取正射影像和数字表面模型基础上，利用数字表面模型计算地表粗糙度，进行阈值分割将其二值化为田埂、非田埂两类对象；再实施形状指数滤波，得到更纯净的田埂二值图；然后以正射影像进行植被分割的二值化图为掩膜剔除植被的影响，修正因植被覆盖影响的边缘突出的田埂二值图；接着通过形态学运算，获取图像边缘更加平滑的田埂二值图，具体实现过程如下：

11)采集耕地区域的无人机序列图像，地面同时布设和量测控制点，通过摄影测量软件处理得到数字表面模型和正射影像，再将目标耕地区域的影像按照典型区域的边界剪裁，然后将正射影像的色差空间从RGB模式变换为HSV模式，对数字表面模型计算得到地表粗糙度f；

式中，地表粗糙度f的值在0～1之间，H为滑动分析窗口的中心像素值，

为滑动分析窗口内的像素均值，n为滑动分析窗口内的像素个数。

12)利用以下公式进行阈值分割，从步骤11)计算的地表粗糙度图像提取耕地田埂二值图f₁：

f(x,y)≥T₁ (2)

其中，地表粗糙度图像中满足公式(2)的像素为耕地田埂，不满足公式(2)的为非耕地田埂，f(x,y)是像素(x,y)处的二进制值，T₁是所采用的分割阈值(地表粗糙度的均值与标准差的一半之和)。

13)对步骤12)所得的候选田埂线计算斑块面积、最小外接矩形周长、最小外接矩形面积和主轴长度4个形状指数，再利用以下公式进行4个形状指数滤波，对步骤12)所得的候选田埂线进行提纯：

其中，S₀代表地表粗糙度阈值分割后的二值图像，S₁代表面积均值过滤后的二值图，S₂代表最小外接矩形的周长均值过滤后的二值图，S₃代表主轴长度均值过滤后的二值图，S₄代表最小外接矩形的面积均值过滤后的二值图，shape area代表像素斑块面积、MERperimeter代表最小外接矩形周长、major axis length代表主轴长度、area ofMERthreshold代表最小外接矩形面积。

其中，f₂是经过四个形状指数过滤后的候选田埂线图像。

14)将步骤11)所得的HSV模式下正射影像选用色调图像，再对经高斯曲线拟合的色调直方图进行阈值检测，并利用以下公式进行二值化分割，即可得到植被掩膜f₃：

f(x,y)≥T₂ (5)

其中，f(x,y)是像素(x,y)处的色调值，T₂为高斯曲线拟合色调直方图所检测到的阈值；

所述分割阈值方法根据色调统计直方图高斯拟合峰值而进行的阈值检测，该方法由Mohamed Hassanein等人2018年提出的(Hassanein,M.；Lari,Z.；El-Sheimy,N.A newvegetation segmentation approach for cropped fields based on thresholddetection from Hue histograms.Sensors 2018,18,1253)。

15)根据步骤14)所得的植被掩膜，利用以下公式对步骤13)滤波结果的候选田埂线进行图像点乘运算，得到剔除植被影响的候选田埂线f₄：

f₄＝f₂·f₃ (6)

其中，f₂是步骤13)形状指数滤波的结果，f₃是步骤14)植被掩膜的结果，f₄是本步骤剔除植被影响的候选田埂线结果。

16)依次利用基于多方向结构元素的图像开运算、小像素斑块移除、基于线性结构元素的图像闭运算和图像细化对步骤15)所得的候选田埂线进行平滑，得到候选田埂线f₅：

利用以下公式构造多方向性结构元素(g₁)：

其中，g₁(x_i,y_i)是像素(x_i,y_i)在g₁处的像素值，α_i是第i个方向的角度值(取45°)，间隔长度i的取值范围从-90°至90°，L代表结构元素分析窗口的长度(取田埂像素宽度)。

利用以下公式进行图像开运算：

其中，g₁代表多方向结构元素，f₄是步骤15)剔除植被影响的结果，f₅是本步骤图像开运算的候选田埂线。

利用以下公式对小像素斑块通过像素对象的面积实现过滤：

其中，T₃是像素斑块面积的过滤阈值(取1000个像素)，f₅是前一公式剔除植被影响的结果，f₆是本公式小斑块过滤的候选田埂线。

利用以下公式进行图像闭运算：

其中，g₂代表线性结构元素(线性结构元素的长度为研究区内田埂像素斑块的主轴长度的一半)，f₆是前一公式小斑块过滤后的结果，f₇是本公式图像闭运算的候选田埂线。

然后利用图像细化算法，生成田埂线的骨架，即为经过形态学运算平滑的耕地田埂线图像(f8)。

所述α_i为45°，所述L设置为田埂像素宽度，由研究区的田埂实际宽度和图像分辨率所决定；所述小斑块面积过滤的阈值选择1000个像素；所述线性结构元素的长度为研究区内田埂像素斑块的主轴长度的一半。

2)对步骤1)所提取的经平滑处理后的耕地田埂二值图利用霍夫变换检测其上的点，再对这些点标记、排序，并根据几何关系去冗余和修缮，然后将点按照所述田埂连接为多段线(即为耕地田埂线)，具体实现过程如下：

21)利用霍夫变换算法对步骤16)所得的田埂线二值图进行检测，需要霍夫变换的参数空间定义、参数空间的峰值投票与识别和霍夫变换空间的交点参数确定三个步骤来提取直线段的两个端点，从而获取田埂线上的若干个点，得到图像坐标系下的点坐标，再根据图像的投影坐标信息，利用以下公式将其坐标由图像坐标系转为投影坐标系：

其中，n_x和n_y代表图像中的像素坐标，r₀代表图像分辨率，x和y代表对应的投影坐标，以及x₀和y₀代表图像左上角拐点的投影坐标。

所述的霍夫变换三个步骤参数设置如下：定义空间参数两项参数：其一检测到的线段角度为研究区所有候选田埂角度中值的±3°范围，其二霍夫空间沿着角度轴线的距离为0.25；投票与识别参数空间的峰值设置为300个，据此确定霍夫变换空间的交点参数，从而将该峰值转为坐标空间的直线段。

22)根据步骤21)所检测的点进行标记，具体过程为：利用以下公式对步骤21)所检测的点进行旋转变换为正南正北方向的田埂线，再根据这些点与每个田埂最小外接矩形的位置关系判断所述田埂，并按照自西向东方向标记点所在的田埂线序号属性；然后按照自南向北方向标记落在每条田埂线上点的序号：

其中，x和y为变换之前的坐标，x′和y′为变换之后的坐标，θ为旋转角度(取研究区所有田埂线的方向角中值)。

所述的田埂线旋转角度为研究区所有候选田埂线的方向角中值。

23)根据几何关系对步骤22)所标记和排序的点进行冗余点去除和缺陷点完善，具体过程为：对每个田埂线上的点计算相邻两点间距，若小于所设定的间距阈值，则这两点坐标被两者坐标均值替代，如此反复迭代直到所有相邻两点间距大于上述阈值为止；在坐标旋转变换为正南正北情况下进行三次检测与完善，包括：其一，如果田埂线的最小外接矩形质心与该田埂线一侧的端点距离更近，与该田埂线的最小外接矩形同侧次轴中心点更远，则该田埂线此侧端点坐标被最小外接矩形的次轴中心坐标替换；其二，拐角点坐标采用对应一侧所有田埂线端点的中值点坐标替换；其三，检测为异常值的端点坐标用该点对应的田埂线与同侧相邻两个田埂线端点连线的交点坐标替换；

24)根据地图制图学原理，对步骤33)所得田埂线上的点按照标记所在的田埂编号和属于每个田埂的点集序号依次连接为矢量田埂线，从而实现直线型耕作田埂的无人机遥感提取。

实施例：

本实施例为华北平原中部某耕地条带田块景观，地势平坦，地表高程约30m。该地区为一年两熟耕作制度，包括冬小麦和夏玉米，粮食产量分别占全国的75％和32％；以水浇地为主，形成了狭长的基本耕作单元景观，田埂线宽度约为0.35m、高度约为0.1m，条带田块长度约为180m、宽度约为3.8～7.5m。

本实施例的一种直线型耕作田埂的无人机遥感提取方法，包括：

1)选择冬小麦生长早期采集单时相耕地的无人机序列图像，在三维重建获取正射影像和数字表面模型基础上，利用数字表面模型计算地表粗糙度，进行阈值分割将其二值化为田埂、非田埂两类对象；再实施形状指数滤波，得到更纯净的田埂二值图；然后以正射影像进行植被分割的二值化图为掩膜剔除植被的影响，修正因植被覆盖影响的边缘突出的田埂二值图；接着通过形态学运算，获取图像边缘更加平滑的田埂二值图，具体实现过程如下：

11)选择2016年11月的冬小麦采集研究区耕地的无人机图像，在Pix4D mapper软件中加入图像和地面控制点，实现三维重建获取正射影像DOM和数字表面模型DSM，两者空间分辨率为2.5cm；再将目标区域耕地影像按研究区边界剪裁，然后将正射影像的色差空间从RGB模式变换为HSV模式，对数字表面模型以本例田埂宽度(0.35m)对应的像素窗口17*17为滑动分析尺寸，利用上述公式(1)计算得到地表粗糙度f；

12)对步骤11)计算的地表粗糙度图像利用上述公式(2)进行阈值分割，阈值为地表粗糙度的均值与标准差的一半之和，从而提取耕地田埂二值图f₁；

13)对步骤12)所得的候选田埂线计算斑块面积、最小外接矩形周长、最小外接矩形面积和主轴长度4个形状指数，再利用上述公式(3)进行4个形状指数滤波，对步骤12)所得的候选田埂线进行提纯，得到二值图f₂；

14)将步骤11)所得的HSV模式下正射影像选用色调图像，再对经高斯曲线拟合的色调直方图进行阈值检测，并利用上述公式(5)以所检测的阈值进行二值化分割，即可得到植被掩膜f₃；

15)根据步骤14)所得的植被掩膜，利用上述公式(6)对步骤13)滤波结果的候选田埂线进行图像点乘运算，得到剔除植被影响的候选田埂线f₄；

16)依次利用基于多方向结构元素的图像开运算公式(7)和公式(8)(方向角为45°和分析窗口的长度取田埂像素宽度17)、小像素斑块移除公式(9)(阈值取1000个像素)、基于线性结构元素的图像闭运算公式(10)(所述线性结构元素的长度为研究区内田埂像素斑块的主轴长度的一半)和图像细化对步骤15)所得的候选田埂线进行平滑，依次得到候选田埂线f₅、f₆、f₇、f₈，并以最后候选田埂线图像f₈为准。

2)对步骤1)所提取的经平滑处理后的耕地田埂二值图利用霍夫变换检测其上的点，再对这些点标记、排序，并根据几何关系去冗余和修缮，然后将点按照所述田埂连接为多段线(即为耕地田埂线)，具体实现过程如下：

21)利用霍夫变换算法对步骤16)所得的田埂线二值图进行检测，经过三个步骤：定义霍夫变换的参数空间(检测到的线段角度θ为研究区所有候选田埂角度中值的±3°范围，霍夫空间沿着角度轴线的距离ρ为0.25)，投票与识别参数空间的峰值(设置为300个)，据此确定霍夫变换空间的交点参数，从而将该峰值转为坐标空间的直线段，可以提取直线段的两个端点，获取田埂线上的若干个点，得到图像坐标系下的点坐标，再根据图像的投影坐标信息，利用上述公式(11)将其点坐标由图像坐标系转为投影坐标系；

22)利用上述公式(12)将步骤21)所检测的点按耕地田埂方位角旋转(取研究区所有田埂线的方向角中值)为正南正北方向，根据这些点与每个田埂最小外接矩形的位置关系判断所述田埂，并按照自西向东方向标记点所在的田埂线序号属性，然后按照自南向北方向标记落在每条田埂线上点的序号；

23)根据几何关系对步骤22)所标记和排序的点进行冗余点去除和缺陷点完善，具体过程为：对每个田埂线上的点计算相邻两点间距，若小于所设定的间距阈值(设定为研究区耕地田埂主轴长度的3％)，则这两点坐标被两者坐标均值替代，如此反复迭代直到所有相邻两点间距大于上述阈值为止；在坐标旋转变换为正南正北情况下进行三次检测与完善，包括：其一，如果田埂线的最小外接矩形质心与该田埂线一侧的端点距离更近，与该田埂线的最小外接矩形同侧次轴中心点更远，则该田埂线此侧端点坐标被最小外接矩形的次轴中心坐标替换；其二，拐角点坐标采用对应一侧所有田埂线端点的中值点坐标替换；其三，检测为异常值的端点坐标用该点对应的田埂线与同侧相邻两个田埂线端点连线的交点坐标替换；

34)根据地图制图学原理，对步骤23)所得田埂线上的点按照标记所在的田埂编号和属于每个田埂的点集序号依次连接为矢量田埂线，从而实现直线型耕作田埂的无人机遥感提取。

本发明有效性验证：

以空间分辨率2.5cm的正射影像为底图，利用目视解译方法提取田埂线和田块面，以此作为验证数据。结果表明，耕地田埂线提取的长度相对误差为-0.16％、召回率98.6％、准确率98.8％，证明了本发明的有效性。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (8)

1.一种直线型耕作田埂的无人机遥感提取方法，其特征在于，包括以下步骤：

1)在采集单时相耕地的无人机序列图像，并三维重建获取正射影像和数字表面模型基础上，利用数字表面模型计算地表粗糙度，进行阈值分割将其二值化为田埂、非田埂两类对象；再实施形状指数滤波，得到更纯净的田埂二值图；然后以正射影像进行植被分割的二值化图为掩膜剔除植被的影响，修正因植被覆盖影响的边缘突出的田埂二值图；接着通过形态学运算，获取图像边缘更加平滑的田埂二值图；

2)对步骤1)所提取的经平滑处理后的耕地田埂二值图利用霍夫变换检测其上的点，对这些点标记、排序，并根据几何关系去冗余和修缮，将点按照所述田埂连接为多段线，该多段线即为耕作田埂线。

2.根据权利要求1所述的一种直线型耕作田埂的无人机遥感提取方法，其特征在于，所述步骤1)的具体过程如下：

11)采集耕地区域的无人机序列图像，地面同时布设和量测控制点，通过摄影测量软件处理得到数字表面模型和正射影像，再将目标耕地区域的影像按照典型区域的边界剪裁，然后将正射影像的色差空间从RGB模式变换为HSV模式，对数字表面模型计算得到地表粗糙度f；

式中，地表粗糙度f的值在0～1之间，H为滑动分析窗口的中心像素值，

为滑动分析窗口内的像素均值，n为滑动分析窗口内的像素个数；

12)利用以下公式进行阈值分割，从步骤11)计算的地表粗糙度图像提取耕地田埂二值图f₁：

f(x,y)≥T₁ (2)

其中，地表粗糙度图像中满足公式(2)的像素为耕地田埂，不满足公式(2)的为非耕地田埂，f(x,y)是像素(x,y)处的二进制值，T₁是所采用的分割阈值；

13)利用以下公式进行4个形状指数滤波，对步骤12)所得的候选田埂线进行提纯：

其中，S₀代表地表粗糙度阈值分割后的二值图像，S₁代表第一次过滤后的二值图，S₂代表第二次过滤后的二值图，S₃代表第三次过滤后的二值图，S₄代表第四次过滤后的二值图，shape area代表像素斑块面积、MER perimeter代表最小外接矩形周长、major axislength代表主轴长度、area ofMER threshold代表最小外接矩形面积；

其中，f₂是经过四个形状指数过滤后的候选田埂线图像；

14)将步骤11)所得的HSV模式下正射影像选用色调图像，再对经高斯曲线拟合的色调直方图进行阈值检测，并利用以下公式进行二值化分割，即可得到植被掩膜f₃：

f(x,y)≥T₂ (5)

其中，f(x,y)是像素(x,y)处的色调值，T₂为所检测到的阈值；

15)根据步骤14)所得的植被掩膜，利用以下公式对步骤13)滤波结果的候选田埂线进行图像点乘运算，得到剔除植被影响的候选田埂线f₄：

f₄＝f₂·f₃ (6)

其中，f₂是步骤13)形状指数滤波的结果，f₃是步骤14)植被掩膜的结果，f₄是本步骤剔除植被影响的候选田埂线结果；

16)依次利用基于多方向结构元素的图像开运算、小像素斑块移除、基于线性结构元素的图像闭运算和图像细化对步骤15)所得的候选田埂线进行平滑，得到候选田埂线f₅：

利用以下公式构造多方向性结构元素(g₁)：

其中，g₁(x_i,y_i)是像素(x_i,y_i)在g₁处的像素值，α_i是第i个方向的角度值，间隔长度i的取值范围从-90°至90°，L代表结构元素分析窗口的长度；

利用以下公式进行图像开运算：

其中，g₁代表多方向结构元素，f₄是步骤15)剔除植被影响的结果，f₅是本步骤图像开运算的候选田埂线；

利用以下公式对小像素斑块通过像素对象的面积实现过滤：

其中，T₃是像素斑块面积的过滤阈值，f₅是前一公式剔除植被影响的结果，f₆是本公式小斑块过滤的候选田埂线；

利用以下公式进行图像闭运算：

其中，g₂代表线性结构元素，f₆是前一公式小斑块过滤后的结果，f₇是本公式图像闭运算的候选田埂线；

然后利用图像细化算法，生成田埂线的骨架，即为经过形态学运算平滑的耕地田埂线图像(f₈)。

3.根据权利要求2所述的一种直线型耕作田埂的无人机遥感提取方法，其特征在于，所述步骤12)，所述分割阈值为研究区地表粗糙度的均值与标准差的一半之和。

4.根据权利要求2所述的一种直线型耕作田埂的无人机遥感提取方法，其特征在于，所述步骤13)，所述像素斑块面积、最小外接矩形周长、主轴长度和最小外接矩形面积，且过滤阈值都取均值。

5.根据权利要求2所述的一种直线型耕作田埂的无人机遥感提取方法，其特征在于，所述步骤16)，所述α_i为45°，所述L设置为田埂像素宽度，由研究区的田埂实际宽度和图像分辨率所决定；所述小斑块面积过滤的阈值选择1000个像素；所述线性结构元素的长度为研究区内田埂像素斑块的主轴长度的一半。

6.根据权利要求1所述的一种直线型耕作田埂的无人机遥感提取方法，其特征在于，所述步骤2)的具体过程如下：

21)利用霍夫变换算法对步骤16)所得的田埂线二值图进行检测，需要霍夫变换的参数空间定义、参数空间的峰值投票与识别和霍夫变换空间的交点参数确定三个步骤来提取直线段的两个端点，从而获取田埂线上的若干个点，得到图像坐标系下的点坐标，再根据图像的投影坐标信息，利用以下公式将其坐标由图像坐标系转为投影坐标系：

其中，n_x和n_y代表图像中的像素坐标，r₀代表图像分辨率，x和y代表对应的投影坐标，以及x₀和y₀代表图像左上角拐点的投影坐标；

22)根据步骤21)所检测的点进行标记，具体过程为：利用以下公式对步骤21)所检测的点进行旋转变换为正南正北方向的田埂线，再根据这些点与每个田埂最小外接矩形的位置关系判断所述田埂，并按照自西向东方向标记点所在的田埂线序号属性；然后按照自南向北方向标记落在每条田埂线上点的序号：

其中，x和y为变换之前的坐标，x′和y′为变换之后的坐标，θ为旋转角度；

23)根据几何关系对步骤22)所标记和排序的点进行冗余点去除和缺陷点完善，具体过程为：对每个田埂线上的点计算相邻两点间距，若小于所设定的间距阈值，则这两点坐标被两者坐标均值替代，如此反复迭代直到所有相邻两点间距大于上述阈值为止；在坐标旋转变换为正南正北情况下进行三次检测与完善，包括：其一，如果田埂线的最小外接矩形质心与该田埂线一侧的端点距离更近，与该田埂线的最小外接矩形同侧次轴中心点更远，则该田埂线此侧端点坐标被最小外接矩形的次轴中心坐标替换；其二，拐角点坐标采用对应一侧所有田埂线端点的中值点坐标替换；其三，检测为异常值的端点坐标用该点对应的田埂线与同侧相邻两个田埂线端点连线的交点坐标替换；

24)根据地图制图学原理，对步骤23)所得田埂线上的点按照标记所在的田埂编号和属于每个田埂的点集序号依次连接为矢量田埂线，从而实现直线型耕作田埂的无人机遥感提取。

7.根据权利要求7所述的一种直线型耕作田埂的无人机遥感提取方法，其特征在于，所述步骤21)，所述的霍夫变换三个步骤参数设置如下：定义空间参数两项参数：其一检测到的线段角度为研究区所有候选田埂角度中值的±3°范围，其二霍夫空间沿着角度轴线的距离为0.25；投票与识别参数空间的峰值设置为300个，据此确定霍夫变换空间的交点参数，从而将该峰值转为坐标空间的直线段。

8.根据权利要求7所述的一种直线型耕作田埂的无人机遥感提取方法，其特征在于，所述步骤22)，所述的田埂线旋转角度为研究区所有候选田埂线的方向角中值。

Images