CN109448046B - 一种基于多描述子的半自动道路中心线快速提取方法 - Google Patents

Abstract

本发明公开了一种基于多描述子的半自动道路中心线快速提取方法，步骤为：输入影像中道路起止点坐标；对原始影像进行L₀滤波；对原始影像进行线段提取；在滤波结果及线段提取结果的基础上，根据起止点坐标，建立多级线段方向直方图，以获得当前道路方向；建立扇形描述子；利用最小二乘法将所有得到的道路跟踪点进行拟合，进一步剔除错误提取的道路点，得到道路中心线。本发明有效减少了使用单一描述子进行道路提取过程中出现的错误提取的现象，只根据道路的起点及终点坐标信息进行跟踪，在跟踪过程中，根据具体的路面情况，自适应地确定道路方向及道路宽度信息，以达到提高算法的自动化程度及算法在复杂路况下的跟踪精度的目的。

Description

一种基于多描述子的半自动道路中心线快速提取方法

技术领域

本发明属于遥感影像信息提取的技术领域，尤其涉及一种基于多描述子的半自动道路中心线快速提取方法。

背景技术

利用遥感影像进行道路提取，在地理信息数据库的更新、军事、和导航方面发挥了重要作用。随着遥感影像的空间分辨率不断提高，这给道路提取带来了便利的同时也带来了巨大的挑战。相比于中分辨率的影像，在高分辨影像中，更多的细节清晰可见(斑马线、人行道及减速带等)，但道路和周围环境之间的差异尤其是水泥建筑物(如房屋、铁路和停车场)并不明显，道路噪声(如建筑物和树木的阴影、交通线路、行人、汽车、等)也不能忽视，这些因素使得基于高分辨率遥感影像进行道路提取更具复杂性和挑战性。

国内外学者对高分辨率遥感影像的道路提取做了很多研究，目前的道路提取方法可分为道路面提取和道路中心线提取。道路面提取主要基于分割和分类。具有代表性的影像分割方法主要有以下几种：主动轮廓模型、均值漂移算法和阈值分割等。然而，由于道路情况的复杂多变，很难找到一种适用于各种类型遥感影像的分割算法。道路中心线提取方法侧重于道路骨架的检测，通常采取细化和跟踪两种不同的方式。细化操作通常是在已提取道路区域上展开的，道路跟踪方法则通常需要提供初始信息，如种子点、起始跟踪方向、道路宽度或匹配规则等信息，通过迭代跟踪的方式得到道路中心线。其缺点在于跟踪过程中方向、位置难以控制，无法在确保道路提取精度的情况下提高算法的自动化程度。

发明内容

基于以上现有技术的不足，本发明所解决的技术问题在于提供一种基于多描述子的半自动道路中心线快速提取方法，达到提高算法的自动化程度及算法在复杂路况下的跟踪精度的目的。

为了解决上述技术问题，本发明通过以下技术方案来实现：

本发明提供一种基于多描述子的半自动道路中心线快速提取方法，包括以下步骤：

步骤1：输入影像中道路起止点坐标；

步骤2：对原始影像进行L₀滤波；

步骤3：对原始影像进行线段提取；

步骤4：在滤波结果及线段提取结果的基础上，根据起止点坐标，建立多级线段方向直方图，以获得当前道路方向：

步骤4.1：在起止点处建立圆形统计区域，根据区域内像素的梯度信息将起止点拟合至对应位置道路中心；

步骤4.2：以拟合后道路中心点为中心，在道路起止点处分别建立矩形搜索区；

步骤4.3：统计矩形搜索区内线段的方向信息，建立线段方向直方图，并根据直方图峰值信息确定当前位置道路跟踪方向；

步骤5：根据拟合后道路中心点及跟踪方向信息，建立扇形描述子；

步骤5.1：根据道路中心点信息设置候选点；

步骤5.2：以道路中心点及候选点为基础，建立基础三角形；

步骤5.3：设置其余6个与基础三角形等大的三角形，构建扇形描述子，设置其余6个候选点；

步骤5.4：统计扇形描述子内三角形内的梯度、灰度及角度信息，选择更符合条件的两个最优候选点，并自适应的确定当前位置道路半径；

步骤5.5：根据道路与非道路区域灰度对比度，确定最优道路点；

步骤5.6：根据较小的距离范围内，跟踪点夹角不会发生较大变化的现象，判断是否保留已确定跟踪点；

步骤5.7：执行步骤4.1-步骤5.6跟踪操作，若不符合跟踪结束条件，则继续跟踪，否则，跟踪结束；

步骤6：在道路提取过程中利用最小二乘法将所有得到的道路跟踪点进行拟合，进一步剔除错误提取的道路点，得到道路中心线。

优选的，所述步骤4.1包括：

步骤4.1.1：以起止点为中心，以一个像素为半径，建立圆形统计区；

步骤4.1.2：分别以起止点周围8邻域所有像素为中心，建立与步骤4.1.1中等大的圆形统计区，统计所有区域内像素梯度之和，取梯度之和最小的统计区域作为当前半径最合适的区域；

步骤4.1.3：将最合适统计区域内像素梯度之和与阈值进行比较，若小于阈值，则以1个像素为单位增大半径，以最合适统计区域对应的圆心为圆心，当前半径为半径，重复步骤4.1.2，直至最合适统计区域内像素梯度之和大于阈值，此时区域的圆心为道路中心，区域的半径为此位置的道路半径。

所述步骤4.2包括：

步骤4.2.1：根据跟踪点数目，确定矩形搜索区中心位置；

步骤4.2.2：根据矩形搜索区中心点位置，建立矩形搜索区。

进一步的，步骤4.2.1包括：

步骤4.2.1.1：当跟踪点数目为0时，直接以起点及终点位置的道路中心为中心，以道路两倍宽度为边长，建立矩形搜索区域。

步骤4.2.1.2：当跟踪点数目大于不为0时，根据跟踪点数确定矩形搜索区域中心位置，将道路的起点及终点分别记为P_S和P_E，起点和终点方向的跟踪点分别记为P_Si和P_Ei，θ_Si和θ_Ei分别为起始和结束方向上的第i跟踪点与第i-1跟踪点之间的夹角。当跟踪点为1时，θ_Si和θ_Ei分别为跟踪点与起点和终点的夹角；在起点方向，θ_i＝θ_Si，在终点方向，θ_i＝θ_Ei，θ＝θ_i。当跟踪点的个数为1-5时，θ是θ₀-θ_i的平均值，当跟踪点的个数大于或等于6时，不考虑θ₀，θ是θ_i-5-θ_i的平均值，通过这种方式，可以精确控制矩形搜索区域的位置，避免较远位置及单点误差对方向搜索位置的误导，使得矩形搜索区域可以始终停留在道路上。

可选的，所述步骤4.3包括：

步骤4.3.1：建立矩形搜索区域，统计该区域中每个方向上的线段的总长度，以15度为单位，将180度分为12等份，建立线段角度直方图；

步骤4.3.2：当出现多峰值情况时，根据各峰值所对应的角度与之前确定的道路方向之间的差值，确定道路方向；

步骤4.3.3：当出现未统计到线段方向的情况时，建立“线段金字塔”，确定道路方向。

可选的，所述步骤5.1的具体方法为：以道路中心点为底边中点(O)，沿着当前道路方向，在步长为三倍道路半径(S)处设置候选点(P)；

所述步骤5.2的具体方法为：以道路中心点为底边中心点，过中心点且垂直于道路方向的2倍道路半径的线段为底边，候选点为顶点，建立基础三角形。

所述步骤5.3的具体方法为：将基础三角形以道路中心点为中心，旋转±15°、±30°和±45°，得到6个与基础三角形等大的三角形组，与基础三角形共同构成扇形描述子，其中，每个三角形的顶点皆为候选点。

可选的，所述步骤5.4包括：

步骤5.4.1：根据角度信息，剔除最不符合条件的2个候选点；

步骤5.4.2：根据三角形内部像素灰度的均匀程度，再次剔除最不符合条件的2个候选点；

步骤5.4.3：根据三角形内部像素梯度信息，剔除最不符合条件的1个候选点。

进一步的，步骤5.4.2包括：

步骤5.4.2.1：以各候选点为圆心，以起点位置道路半径为直径，建立圆形统计区，计算每个圆形统计区内所有像素的灰度均值，记为统计灰度值；

步骤5.4.2.2：当跟踪点为0时，只考虑起止点。以起止点为圆心，以起止点位置道路半径为半径建立圆形参考区，计算参考区内像素灰度均值作为参考灰度值；当跟踪点不为0时，计算所有跟踪点及道路起止点参考区(半径为对应位置处的道路半径)内像素灰度均值，作为参考灰度值。

步骤5.4.2.3：计算每个候选点对应的统计灰度值与参考灰度的差值，将最大差值对应的两个候选点视为非道路点并将其剔除。

可选的，所述步骤5.5包括：

步骤5.5.1：以两最优候选点为圆心，区域半径为半径，建立圆形统计区，分别计算统计区内所有像素的灰度均值，分别记为M_S1，M_S2；

步骤5.5.2：分别计算两候选点与跟踪点之间的夹角，记为α_S1，α_S2，分别在以候选点为中心，垂直于α₁、α₂方向，候选点两侧2倍区域半径位置设置参考点；

步骤5.5.3：以参考点为圆心，区域半径为半径，建立参考区；分别计算统计参考区内所有像素的灰度均值，记为M_S1R1，M_S1R2，M_S2R1，M_S2R2；

步骤5.5.4：计算M_S1分别于与M_S1R1，M_S1R2之间的差值并判断其是否大于阈值，如果至少有一个差值大于阈值，则将其视为待定点，同对M_S1进行相同操作，根据待定点信息选取跟踪点。

进一步的，步骤5.5.4包括：

步骤5.5.4.1：若待定点数目为0，则将此位置视为车辆等噪声干扰区域，将步长(step)设置为之前的两倍，重新执行步骤4.1-步骤5.5跟踪操作，最多可重复进行此操作两次，若待定点数目依然为0，则将此位置视为道路终点，终止跟踪。

步骤5.5.4.2：若待定点数目为1，则将此待定点视为跟踪点。

步骤5.5.4.3：若待定点数目为2，计算两个待定点与最新跟踪点之间的夹角与θ(步骤4.2.1.2)之间的差值，将差值较小者视为跟踪点。

可选的，所述步骤5.6的具体方法为：计算最新确定的跟踪点与前一个跟踪点之间的夹角与步骤4.2.1中的θ之间的差值，若差值小于阈值，则保留跟踪点；若大于阈值，则将此跟踪点视为错误提取的跟踪点或由于车辆等噪声干扰偏离道路中心，则将此跟踪点删除并将步长值(step)设置成之前的两倍，重新执行步骤4.1-步骤5.5跟踪操作，得到新的跟踪点。

可选的，所述步骤5.7的具体方法为：每循环一次，计算最新两个最优道路点之间的距离，若距离大于阈值，则继续跟踪，直至距离小于阈值或跟踪至道路断裂处，跟踪结束。

可选的，所述步骤6的具体方法为：

根据给定的M点，使用曲线y＝f(x)的近似曲线y＝φ(x)，具体过程下：

设置目标函数：y＝a₀+a₁x+L+a_kx^k；

计算每个点到该曲线的距离的总和，即平方差的和：

为了获得满足条件的值，确定方程右侧的a_i的偏导数，并以矩阵的形式表示它，可以得到以下矩阵：

对此矩阵进行化简，可得到以下矩阵：

即X*A＝Y，然后A＝(X′*X)-1*X′*Y，得到系数矩阵A，同时得到拟合曲线，即道路中心线。

由上，本发明的基于多描述子的半自动道路中心线快速提取方法至少具有如下有益效果：

(1)本发明提出了一种多级线段方向直方图(MLSOH)描述子，线性特征是道路的一种重要的结构特征，特别是在一个小的视觉区域，道路的中心线可以用线段表达，因此，本发明提出了一种在线段的基础进行道路跟踪的思想。

(2)在道路跟踪过程中，本发明将多种描述子结合在一起：利用MLSOH描述子确定道路方向，同时利用三角形描述子验证道路方向并确定最优道路点，有效减少了使用单一描述子进行道路提取过程中出现的错误提取的现象。

(3)本发明不设置初始道路的宽度及跟踪方向等信息，只根据道路的起点及终点坐标信息进行跟踪，在跟踪过程中，根据具体的路面情况，自适应地确定道路方向及道路宽度信息，以达到提高算法的自动化程度及算法在复杂路况下的跟踪精度的目的。

(4)本发明从道路两侧端点同时向中央进行跟踪，利用对向跟踪方向对跟踪方向进行动态约束，以避免发生错误跟踪问题。

上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，而可依照说明书的内容予以实施，并且为了让本发明的上述和其他目的、特征和优点能够更明显易懂，以下结合优选实施例，并配合附图，详细说明如下。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例的附图作简单地介绍。

图1为本发明的多描述子的半自动道路中心线快速提取方法的流程图；

图2为本发明中将种子点拟合至道路中心示意图；

图3为本发明中确定矩形搜索区中心点位置的角度约束示意图；

图4为本发明中理想条件下直方图的建立示意图，其中，(a)为线段提取结果与原图叠加显示图，(b)为矩形搜索区域内线段放大显示图，(c)为矩形搜索区对应的线段方向直方图；

图5为本发明中非理想条件下直方图的建立示意图，其中，(a)为线段提取结果与原图叠加图，(b)为矩形搜索区域内线段放大显示，(c)为矩形搜索区对应的线段方向直方图，Q1为噪声干扰造成的多峰值情况，Q2为阴影区域；

图6为本发明中金字塔的建立示意图；

图7为本发明中实验一的结果对比图，其中，(a)为原始影像图，(b)为真实道路路面结果，(c)为本发明道路提取结果及对应区域放大图，(d)为eCognition软件面向对象分类方法道路提取结果图，(e)为ERDAS软件道路提取结果图，(f)为T型模板道路提取结果图；

图8为本发明中实验二的结果对比图，其中，(a)为原始影像图，(b)为真实道路路面结果，(c)为本发明道路提取结果及对应区域放大图，(d)为eCognition软件面向对象分类方法道路提取结果图，(e)为ERDAS软件道路提取结果图，(f)为T型模板道路提取结果图；

图9为本发明中实验三的结果对比图，其中，(a)为原始影像图，(b)为真实道路路面结果，(c)为本发明道路提取结果及对应区域放大图，(d)为eCognition软件面向对象分类方法道路提取结果图，(e)为ERDAS软件道路提取结果图，(f)为T型模板道路提取结果图；

图10为本发明中实验四的结果对比图，其中，(a)为原始影像图，(b)为真实道路路面结果，(c)为本发明道路提取结果及对应区域放大图，(d)为eCognition软件面向对象分类方法道路提取结果图，(e)为ERDAS软件道路提取结果图，(f)为T型模板道路提取结果图；

图11为本发明中的四组实验各种道路提取算法选取的种子点数量示意图。

图12为本发明中扇形描述子的构建示意图。

具体实施方式

下面结合附图详细说明本发明的具体实施方式，其作为本说明书的一部分，通过实施例来说明本发明的原理，本发明的其他方面、特征及其优点通过该详细说明将会变得一目了然。在所参照的附图中，不同的图中相同或相似的部件使用相同的附图标号来表示。

如图1至图12所示，本发明的一种基于多描述子的半自动道路中心线快速提取方法，包括如下步骤：

步骤1：输入影像中道路起止点坐标；

在道路起点及终点位置、无阴影及车辆遮挡、道路内部像素灰度均匀、道路边界清晰处，选取道路种子点。

步骤2：对原始影像进行L₀滤波；

L₀滤波算法通过去除小的非零梯度，去除不重要的细节信息，增强图像的显著性，从而在去除噪声的同时很好地保持了道路边界信息。

步骤3：对原始影像进行线段提取。

步骤4：在滤波结果及线段提取结果的基础上，根据起止点坐标，建立多级线段方向直方图(MLSOH)，以获得当前道路方向，具体步骤包括：

步骤4.1：如图2所示，在起止点处建立圆形统计区域，根据区域内像素的梯度信息将起止点拟合至对应位置道路中心；

步骤4.1.1：以起止点为中心，以一个像素为半径，建立圆形统计区；

步骤4.1.2：分别以起止点周围8邻域所有像素为中心，建立与步骤4.1.1中等大的圆形统计区，统计所有区域内像素梯度之和，取梯度之和最小的统计区域作为当前半径最合适的区域；

步骤4.1.3：将最合适统计区域内像素梯度之和与阈值进行比较，若小于阈值，则以1个像素为单位增大半径，以最合适统计区域对应的圆心为圆心，当前半径为半径，重复步骤4.1.2，直至最合适统计区域内像素梯度之和大于阈值，此时区域的圆心为道路中心，区域的半径为此位置的道路半径；

步骤4.2：以拟合后道路中心点为中心，在道路起止点处分别建立矩形搜索区；

步骤4.2.1：根据跟踪点数目，确定矩形搜索区中心位置；

步骤4.2.1.1：当跟踪点数目为0时，直接以起点及终点位置的道路中心为中心，以道路两倍宽度为边长，建立矩形搜索区域；

步骤4.2.1.2：当跟踪点数目大于不为0时，根据跟踪点数确定矩形搜索区域中心位置，如图3所示，将道路的起点及终点分别记为P_S和P_E，起点和终点方向的跟踪点分别记为P_Si和P_Ei，θ_Si和θ_Ei分别为起始和结束方向上的第i跟踪点与第i-1跟踪点之间的夹角。当跟踪点为1时，θ_Si和θ_Ei分别为跟踪点与起点和终点的夹角；在起点方向，θ_i＝θ_Si，在终点方向，θ_i＝θ_Ei，θ＝θ_i。当跟踪点的个数为1-5时，θ是θ₀-θ_i的平均值，当跟踪点的个数大于或等于6时，不考虑θ₀，θ是θ_i-5-θ_i的平均值，通过这种方式，可以精确控制矩形搜索区域的位置，避免较远位置及单点误差对方向搜索位置的误导，使得矩形搜索区域可以始终停留在道路上。根据以下公式计算搜索区域的中心点位置：

X_Rectangle＝x_reference+step×cosθ

Y_Rectangle＝y_reference+step×sinθ

其中，(X_Rectangle，Y_Rectangle)为矩形搜索区域中心坐标，(x_reference，y_reference)为已确定的最新跟踪点坐标，step为步长，初始步长长度为2倍道路半径。

步骤4.2.2：根据矩形搜索区中心点位置，建立矩形搜索区；

以矩形搜索区域中心为中心，以2倍道路宽度为边长，建立矩形搜索区域。

步骤4.3：统计矩形搜索区内线段的方向信息，建立线段方向直方图，并根据直方图峰值信息确定当前位置道路跟踪方向；

步骤4.3.1：建立矩形搜索区域，统计该区域中每个方向上的线段的总长度，如图4所示，以15度为单位，将180度分为12等份，建立线段角度直方图。理想情况下，直方图峰值即为此时的道路方向。

步骤4.3.2：当出现多峰值情况时，根据各峰值所对应的角度与之前确定的道路方向之间的差值，确定道路方向。

当直方图出现多峰值情况时，由于即使在大弯道的道路中，道路方向在较小步长范围内也不会发生剧烈变化，因此计算峰值对应的角度与之前确定的道路方向之间的差。如果差值小于阈值，则将该角度用作新的道路方向；否则，将该峰值视为干扰值，并且类似地判断次级峰值，直到确定道路方向。

步骤4.3.3：当出现未统计到线段方向的情况时，建立“线段金字塔”，确定道路方向。

当遇到阴影等情况时，直方图会出现没有统计到道路方向的情况，此时，通过构建“线段金字塔”解决此问题。与图像金字塔类似，线段金字塔将由0层图像检测到的线段信息投影到第一层和第二层。在0层，由于阴影的存在，矩形搜索区域无法检测线段，于是利用相同大小的矩形搜索区域在1级线段金字塔中进行线段检测。如果检测到线段，则建立线段角度直方图，确定道路方向。否则，继续进行二级线段金字塔的检测。如果仍然没有检测到线段，它们此位置视为道路的终点，停止跟踪。

步骤5：根据拟合后道路中心点及跟踪方向信息，建立扇形描述子，以验证道路方向及确定最佳道路点；

步骤5.1：根据道路中心点信息设置候选点；

如图12所示，以道路中心点为底边中点(O)，沿着当前道路方向，在步长为三倍道路半径(S)处设置候选点(P)；

步骤5.2：以道路中心点及候选点为基础，建立基础三角形；

以道路中心点为底边中心点，过中心点且垂直于道路方向的2倍道路半径的线段为底边，候选点为顶点，建立基础三角形；

步骤5.3：设置其余6个与基础三角形等大的三角形，构建扇形描述子，设置其余6个候选点；

将基础三角形以道路中心点为中心，旋转±15°、±30°和±45°，得到6个与基础三角形等大的三角形组，与基础三角形共同构成扇形描述子，其中，每个三角形的顶点皆为候选点；

步骤5.4：统计扇形描述子内三角形内的梯度、灰度及角度信息，选择更符合条件的两个最优候选点，并自适应的确定当前位置道路半径；

步骤5.4.1：根据角度信息，剔除最不符合条件的2个候选点；

以起点方向为例，将终点方向的最后一个跟踪点标记为P_end，计算三角形底边中心点O和P_end之间的夹角，记为θ_reference，同时，计算每个候选点与P_end之间的夹角，记为θ_1-7。分别计算θ_reference与θ₁-θ₇之间的角度差，剔除与θ_reference差值最大两个角度对应的候选点。

步骤5.4.2：根据三角形内部像素灰度的均匀程度，再次剔除最不符合条件的2个候选点；

一条道路内部像素的灰度值趋于稳定，因此，根据候选点位置道路像素灰度均值与已确定的所有跟踪点位置的灰度值进行比较，即可进一步剔除不符合道路条件的候选点。具体过程如下：

步骤5.4.2.1：以各候选点为圆心，以起点位置道路半径为直径，建立圆形统计区，计算每个圆形统计区内所有像素的灰度均值，记为统计灰度值；

步骤5.4.2.2：当跟踪点为0时，只考虑起止点。以起止点为圆心，以起止点位置道路半径为半径建立圆形参考区，计算参考区内像素灰度均值作为参考灰度值；当跟踪点不为0时，计算所有跟踪点及道路起止点参考区(半径为对应位置处的道路半径)内像素灰度均值，作为参考灰度值。

步骤5.4.2.3：计算每个候选点对应的统计灰度值与参考灰度的差值，将最大差值对应的两个候选点视为非道路点并将其剔除。

步骤5.4.3：根据三角形内部像素梯度信息，剔除最不符合条件的1个候选点；

道路内部像素灰度均匀，因此，在三角形统计区内，所有像素梯度之和越小，对应的候选点属于道路的可能性越大。因此，将剔除梯度之和最大的三角形区域所对应的候选点，余下两点即为可能性最大的道路点。

步骤5.4.4：将两个点利用步骤4.1提到的方法拟合至对应区域中心，同时，记录统计区域半径信息。

步骤5.5：根据道路与非道路区域灰度对比度，确定最优道路点。

通常情况下，由于道路与其两侧非道路位置具有明显的灰度差异，根据此信息，确定最优道路点。

步骤5.5.1：以两最优候选点为圆心，区域半径为半径，建立圆形统计区，分别计算统计区内所有像素的灰度均值，分别记为M_S1，M_S2。

步骤5.5.2：分别计算两候选点与跟踪点之间的夹角，记为α_S1，α_S2，分别在以候选点为中心，垂直于α₁、α₂方向，候选点两侧2倍区域半径位置设置参考点。

步骤5.5.3：以参考点为圆心，区域半径为半径，建立参考区。分别计算统计参考区内所有像素的灰度均值，记为M_S1R1，M_S1R2，M_S2R1，M_S2R2。

步骤5.5.4：计算M_S1分别于与M_S1R1，M_S1R2之间的差值并判断其是否大于阈值，若至少有一个差值大于阈值，则将其视为待定点，同对M_S1进行相同操作，根据待定点信息选取跟踪点。

步骤5.5.4.1：若待定点数目为0，则将此位置视为车辆等噪声干扰区域，将步长(step)设置为之前的两倍，重新执行步骤4.1-步骤5.5跟踪操作，最多可重复进行此操作两次，若待定点数目依然为0，则将此位置视为道路终点，终止跟踪。

步骤5.5.4.2：若待定点数目为1，则将此待定点视为跟踪点。

步骤5.5.4.3：若待定点数目为2，计算两个待定点与最新跟踪点之间的夹角与θ(步骤4.2.1.2)之间的差值，将差值较小者视为跟踪点。

步骤5.6：根据较小的距离范围内，跟踪点夹角不会发生较大变化的现象，判断是否保留已确定跟踪点。

计算最新确定的跟踪点与前一个跟踪点之间的夹角与θ(步骤4.2.1.2)之间的差值，若差值小于阈值，则保留跟踪点；若大于阈值，则将此跟踪点视为错误提取的跟踪点或由于车辆等噪声干扰偏离道路中心，则将此跟踪点删除并将步长值(step)设置为之前的两倍，重新执行步骤4.1-步骤5.5跟踪操作，得到新的跟踪点。

步骤5.7：执行步骤4.1-步骤5.5跟踪操作，若不符合跟踪结束条件，则继续跟踪，否则，跟踪结束。

每循环一次，计算最新两个最优道路点之间的距离，若距离大于阈值，则继续跟踪，直至距离小于阈值或跟踪至道路断裂处，跟踪结束。

步骤6：由于道路的复杂性和成像过程中的噪声，在道路提取过程中不可避免地产生错误提取。因此，利用最小二乘法将所有得到的道路跟踪点进行拟合，进一步剔除错误提取的道路点，得到道路中心线。

根据给定的M点，该方法不要求曲线f(x)精确地传递这些点。而是使用曲线y＝f(x)的近似曲线y＝φ(x)。具体过程下：

设置目标函数：y＝a₀+a₁x+L+a_kx^k；

计算每个点到该曲线的距离的总和，即平方差的和：

为了获得满足条件的值，确定方程右侧的a_i的偏导数，并以矩阵的形式表示它。可以得到以下矩阵：

对此矩阵进行化简，可得到以下矩阵：

即X*A＝Y，然后A＝(X′*X)-1*X′*Y，得到系数矩阵A，同时得到拟合曲线，即道路中心线。

根据公式，多项式次数越多，结果就越准确。然而，精度高，计算量也随之增大。在实验中，在考虑精度的同时，还应考虑实验效率。经过实验验证，当多项式数为k＝9时，得到了较好的实验结果。

以下通过实验验证本发明所提方法的性能，并与现有经典软件及算法进行比较，分别选取高分2号影像和法国PLEIADES卫星影像，共四幅不同侧重点及分辨率的影像数据作为实验数据，以验证本发明算法结果的可靠性。

图7为高分2号影像数据，覆盖城市区域，影像大小为2000*2000像素，空间分辨率为0.8m。此影像道路情况复杂，路面与背景的差异不明显，许多道路路面被建筑物和树木的阴影阻挡，部分路段甚至被建筑物阴影完全遮挡，单纯依靠像素信息提取道路是非常困难的。从图7c可以看出，本发明提出的方法对被遮挡道路有很好的提取效果，可以对半遮挡道路和大部分完全遮挡道路的中心线进行正确提取。但是，对于图7c的第一张放大图，阴影面积太大；道路的准确方向不能由MLSOH描述符确定；同时，阴影的宽度过大，缺乏参考纹理信息，因此，扇区描述符无法获得准确的道路点。在通过此阴影区域后，本发明算法通过角度控制，及时修正跟踪点，使整体的道路提取结果误差不大。图7d为利用eCognition软件中的面向对象方法进行道路提取的结果图。在此图中，很多道路没有完全提取出来，且阴影区域均未检测到路面，整张影像的道路提取效果较差。

图8为高分2号影像数据，影像大小为2000*2000像素，空间分辨率为0.8m，主要覆盖高速公路区域。这幅影像由许多变分道路组成，在这些路段，道路由双向行驶的道路变为单向行驶道路，道路宽度会相应减小。同时，较大的道路弯度也增加了道路提取的难度。本发明利用扇区描述符和自适应道路宽度方法来提取道路，有效地避免了变分道路错误提取的出现。从图8c中可以看出，在处理变分道路时，本发明方法优先提取小曲率的道路，然后重新选择种子点，实现变分道路中另一条路的提取。

图9为法国PLEIADES卫星影像数据，覆盖农村区域，影像大小为2000*2000像素，空间分辨率0.5m。这幅影像中的主要的道路线是由一条的弯路和两条直路组成，弯路由一个90°的弯路，一个135°的弯路和两个小弯路组成(图9c放大图)，这些弯路的提取是对本发明算法的一个相当大的挑战，路边树木的阴影也对本算法造成了很大的干扰。本发明使用线性约束的算法确定道路方向，利用起点和终点同时向中间进行跟踪的方式，动态约束跟踪点，有效解决了此问题。通过观察影像中间曲率较小的垂直道路的特征，发现这条道路的宽度很小，不属于主干道，不进行提取。图9c显示了本发明算法的道路提取效果。在另外两种半自动算法中，需要在弯曲处选择更多的点来控制方向。

图10为法国PLEIADES卫星影像数据，空间分辨率为1.0m，影像大小为2000*2000像素。在这幅影像中，本发明主要讨论了算法对环形交叉口的提取效果。在小范围内，环形交叉口的曲率有很大的变化。环形交叉口连接了所有的道路，进入环形交叉口后跟踪方向很难控制。例如，经过分析，一条道路更符合跟踪条件的要求，而在现实中，另一条路才是需要跟踪道路。此外，环形交叉口的边界不清晰，沿MLSOH预测的方向进行跟踪可能会导致跟踪点在不进入转盘的情况下直接向另一个方向进行跟踪，从而导致偏移预设跟踪方向的现象的发生；同时，道路斑马线和减速带对道路提取也有一定的影响。在本发明中，通过增加搜索步长的方式跨越减速带。通过观察可以发现，影像下方的水平道路在环形交叉口两侧呈现不同的特征：环形交叉口左侧道路为双向四车道，中间有栅栏，因此，这段路被认为是两条路；环形交叉口右侧为双向双车道，中间没有围栏，因此，这段路被认为是一条路。本发明是基于路的右侧进行提取，因此在环形交叉口下方存在缺失现象。

表1为四种不同方法对四组实验影像进行试验的统计结果，“COM”是指“完整性”、“CORR”是指“正确性”，“RMS”体现的是检测结果到道路中心线的均方差。在四组实验中，四幅影像中道路的宽度大约是14至15个像素，部分道路宽度达到25个像素，因此，2-3像素的偏移量对实验结果影响不大。因此，本发明提出的算法是完全可靠的。在四种不同侧重点的不同类型影像中，与其他三种算法相比，本发明的算法在自动化和提取精度方面具有很大优势。

表1：各种算法的道路提取效果评价表

以上所述是本发明的优选实施方式而已，当然不能以此来限定本发明之权利范围，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和变动，这些改进和变动也视为本发明的保护范围。

Claims (9)

1.一种基于多描述子的半自动道路中心线快速提取方法，其特征在于，包括以下步骤：

步骤1：输入影像中道路起止点坐标；

步骤2：对原始影像进行L₀滤波；

步骤3：对原始影像进行线段提取；

步骤4：在滤波结果及线段提取结果的基础上，根据起止点坐标，建立多级线段方向直方图，以获得当前道路方向：

步骤4.1：在起止点处建立圆形统计区域，根据区域内像素的梯度信息将起止点拟合至对应位置道路中心；

步骤4.2：以拟合后道路中心点为中心，在道路起止点处分别建立矩形搜索区；

步骤4.3：统计矩形搜索区内线段的方向信息，建立线段方向直方图，并根据直方图峰值信息确定当前位置道路跟踪方向；

步骤4.3.1：建立矩形搜索区域，统计该区域中每个方向上的线段的总长度，以15度为单位，将180度分为12等份，建立线段角度直方图；

步骤4.3.2：当出现多峰值情况时，根据各峰值所对应的角度与之前确定的道路方向之间的差值，确定道路方向；

步骤4.3.3：当出现未统计到线段方向的情况时，建立“线段金字塔”，确定道路方向；

步骤5：根据拟合后道路中心点及跟踪方向信息，建立扇形描述子；

步骤5.1：根据道路中心点信息设置候选点(P)；

步骤5.2：以道路中心点及候选点为基础，建立基础三角形；

步骤5.3：设置其余6个与基础三角形等大的三角形，构建扇形描述子，设置其余6个候选点；

步骤5.4：统计扇形描述子内三角形内的梯度、灰度及角度信息，选择更符合条件的两个最优候选点，并自适应的确定当前位置道路半径；

步骤5.5：根据道路与非道路区域灰度对比度，确定最优道路点；

步骤5.6：根据距离范围内，跟踪点夹角不会发生变化的现象，判断是否保留已确定跟踪点；

步骤5.7：执行步骤4.1-步骤5.6跟踪操作，若不符合跟踪结束条件，则继续跟踪，否则，跟踪结束；

步骤6：在道路提取过程中利用最小二乘法将所有得到的道路跟踪点进行拟合，进一步剔除错误提取的道路点，得到道路中心线。

2.如权利要求1所述的基于多描述子的半自动道路中心线快速提取方法，其特征在于，所述步骤4.1包括：

步骤4.1.1：以起止点为中心，以一个像素为半径，建立圆形统计区；

步骤4.1.2：分别以起止点周围8邻域所有像素为中心，建立与步骤4.1.1中等大的圆形统计区，统计所有区域内像素梯度之和，取梯度之和最小的统计区域作为当前半径最合适的区域；

步骤4.1.3：将最合适统计区域内像素梯度之和与阈值进行比较，若小于阈值，则以1个像素为单位增大半径，以最合适统计区域对应的圆心为圆心，当前半径为半径，重复步骤4.1.2，直至最合适统计区域内像素梯度之和大于阈值，此时区域的圆心为道路中心，区域的半径为此位置的道路半径。

3.如权利要求1所述的基于多描述子的半自动道路中心线快速提取方法，其特征在于，所述步骤4.2包括：

步骤4.2.1：根据跟踪点数目，确定矩形搜索区中心位置；

步骤4.2.2：根据矩形搜索区中心点位置，建立矩形搜索区。

4.如权利要求1所述的基于多描述子的半自动道路中心线快速提取方法，其特征在于，所述步骤5.1的具体方法为：以道路中心点为底边中点(O)，沿着当前道路方向，在步长为三倍道路半径(S)处设置候选点(P)；

所述步骤5.2的具体方法为：以道路中心点为底边中心点，过中心点且垂直于道路方向的2倍道路半径的线段为底边，候选点为顶点，建立基础三角形。

5.如权利要求1所述的基于多描述子的半自动道路中心线快速提取方法，其特征在于，所述步骤5.3的具体方法为：将基础三角形以道路中心点为中心，旋转±15°、±30°和±45°，得到6个与基础三角形等大的三角形组，与基础三角形共同构成扇形描述子，其中，每个三角形的顶点皆为候选点。

6.如权利要求1所述的基于多描述子的半自动道路中心线快速提取方法，其特征在于，所述步骤5.4包括：

步骤5.4.1：根据角度信息，剔除最不符合条件的2个候选点；

步骤5.4.2：根据三角形内部像素灰度的均匀程度，再次剔除最不符合条件的2个候选点；

步骤5.4.3：根据三角形内部像素梯度信息，剔除最不符合条件的1个候选点。

7.如权利要求1所述的基于多描述子的半自动道路中心线快速提取方法，其特征在于，所述步骤5.5包括：

步骤5.5.1：以两最优候选点为圆心，区域半径为半径，建立圆形统计区，分别计算统计区内所有像素的灰度均值，分别记为M_S1，M_S2；

步骤5.5.2：分别计算两候选点与跟踪点之间的夹角，记为α_S1，α_S2，分别在以候选点为中心，垂直于α₁、α₂方向，候选点两侧2倍区域半径位置设置参考点；

步骤5.5.3：以参考点为圆心，区域半径为半径，建立参考区；分别计算统计参考区内所有像素的灰度均值，记为M_S1R1，M_S1R2，M_S2R1，M_S2R2；

步骤5.5.4：计算M_S1分别于与M_S1R1，M_S1R2之间的差值并判断其是否大于阈值，如果至少有一个差值大于阈值，则将其视为待定点，同对M_S1进行相同操作，根据待定点信息选取跟踪点。

8.如权利要求1所述的基于多描述子的半自动道路中心线快速提取方法，其特征在于，所述步骤5.6的具体方法为：计算最新确定的跟踪点与前一个跟踪点之间的夹角与步骤4.2.1中的θ之间的差值，若差值小于阈值，则保留跟踪点；若大于阈值，则将此跟踪点视为错误提取的跟踪点或由于车辆噪声干扰偏离道路中心，则将此跟踪点删除并将步长值(step)设置成之前的两倍，重新执行步骤4.1-步骤5.5跟踪操作，得到新的跟踪点。

9.如权利要求1所述的基于多描述子的半自动道路中心线快速提取方法，其特征在于，所述步骤5.7的具体方法为：每循环一次，计算最新两个最优道路点之间的距离，若距离大于阈值，则继续跟踪，直至距离小于阈值或跟踪至道路断裂处，跟踪结束。

Images