CN103258203A - 遥感影像的道路中线自动提取方法 - Google Patents

Abstract

本发明涉及一种遥感影像的道路中线自动提取方法，分割影像，将光谱特征近似的道路区域与周围对比度明显的地物分割为不同区域；利用形状指数分析从分割斑块中分类识别道路区域；利用一种正交向排序法探测道路中心线基元；最后将不连续的中线道路基元迭代式共线链编组为道路网络。本发明具有易于实现、操作简单，准确率高等优点，可以为遥感影像中道路矢量中心线的提取提供一种快速便捷的自动获取方式。

Description

遥感影像的道路中线自动提取方法

技术领域

本发明涉及遥感测绘领域，尤其是一种遥感影像的道路中线自动提取方法。

背景技术

道路作为地理信息中重要的枢纽，与交通、水利、城市规划、电力、应急响应、能源等各行各业的建设密切相关。遥感影像将三维现实场景转换为二维场景，在计算机内存中以二维数组的形式存储，通过不直接接触现实物体的形式来研究现实世界。基于遥感影像的道路识别与提取研究一直是近二十年来的重要研究课题，有许多经典、有效的提取算法被用于实现道路特征的提取，如hough变换、主动轮廓模型、动态规划、模板匹配等。这些算法基本都是利用道路的光度和几何特性识别道路特征。Mena（2003）依据特征线索的不同将道路提取方法分为两类：①基于光谱和纹理分类的提取方法，先利用道路表面同一的光谱及纹理特性识别道路区域，再提取矢量，例如角度纹理法，利用矩形模板绕着固定点旋转，获取模板中纹理特征，从而判断是否为道路类别（张睿,2008;Haverkamp,2002;周绍光,2010）。道路特征聚类、图像分割与机器学习（Roggero,2002;Wenzhong Shi,2002;Sukhendu,2011）等同样属于这一类提取方法；②基于几何分析的方法，主要利用道路边缘、几何形状提取结果，例如边缘检测（唐亮,2005;Wenzhong Shi,2002）、模板匹配(Hu,2005;Taejung Kim,2004;XiangyunHu,2007)。但是综合现有的研究方法，提取完整正确的道路网仍处于研究探索阶段，这主要是因为道路特征的多样（不同道路类型，如形状及宽度）、场景的复杂（周围地物种类繁多：如房屋、停车场；路面地物多样：如汽车、标记线）以及严重遮挡与阴影导致。

发明内容

本发明的目的是提供一种遥感影像中道路中心线探测和提取的方法，该方法实现了基于遥感影像特征提取中道路中线的自动提取和矢量化。

本发明的技术方案是一种遥感影像的道路中线自动提取方法，包括以下步骤：

步骤1，输入原始遥感影像数据，分割为若干区域；

步骤2，通过对步骤1所得各区域进行形状指数分析，判断道路区域，并生成二值道路影像；

步骤3，从二值化道路影像中探测道路中线，并跟踪道路中线，获取道路基元；

步骤4，通过道路基元编组获得自动提取结果，包括以下子步骤，

步骤4.1，根据步骤3所得道路基元建立初始的道路基元连接矩阵，

设道路基元连接矩阵记为A，A为二维对称矩阵，矩阵元素为对应两个道路基元的连接概率，矩阵元素计算公式如下，

其中，道路基元个数为n+1，u的取值为0,1,…n，v的取值为0,1,…n，θu、θv分别为道路基元所在线段lu与lv的中点连线与lu、lv所形成的夹角，α是线段lu与lv之间的夹角，g是线段lu与lv之间的断裂长度；Max_|θu|+|θv|为|θu|和|θv|之和的可能最大值，T_g为断裂间隔的阈值，l_u是线段lu的长度，l_v是线段lv的长度；

步骤4.2，将矩阵元素大于给定阈值的道路基元编组连接，若没有能够连接的道路基元结束编组过程，否则进入步骤4.3；

步骤4.3，利用步骤4.2编组连接后新生成的道路基元重新建立道路基元连接矩阵，返回步骤4.2。

而且，步骤1中采用Meanshift算法进行分割。

而且，步骤2中所述形状指数包括区域的面积S、区域的最小面积外接矩形的长宽比R、区域的充满度F＝S/S_MAER，满足如下公式的区域被判断为道路区域，并生成二值道路影像，

S＜T_S∩(F＜T_F∪R＞T_R)

其中，T_S为预设的面积阈值，T_F为预设的长宽比阈值，T_R为预设的充满度阈值。

而且，步骤3实现方式为，

建立三个窗口Ls、Le、Lb，其中Ls表示一维的探测窗口大小，Le表示统计窗口大小，Lb表示待判断像素pi的邻近范围，Lb<Le<Ls；执行如下步骤，

步骤3.1，输入二值道路影像，平滑处理；

步骤3.2，探测窗口沿水平与垂直方向依次遍历探测道路区域的中心线像素，并取两个结果的并集；探测方式如下，

设像素pi为探测窗口遍历到的中心像素，将像素pi与邻近范围内的Lb个像素都存入数组T_w(k)，其中k=0,1..Lb；对探测窗口内Ls个像素灰度值排序，将排在前面Le+1位的像素存入数组S_w(t)，t=0,1…Le；比较T_w(k)和S_w(t)，如果满足T_w(k)∈S_w(t)，k=0,1..Lb，则将像素pi判断为道路区域的中心线像素；

步骤3.3，对步骤3.2所得道路区域的中心线像素进行细化，获取单像素的道路中线；

步骤3.4，跟踪步骤3.3所得道路中线，获取道路基元的矢量。

本发明具有易于实现、操作简单，准确率高等优点，可以为遥感影像中道路矢量中心线的提取提供一种快速便捷的自动获取方式。

附图说明

图1为本发明实施例的形状指数示意图之一。

图2为本发明实施例的形状指数示意图之二。

图3为本发明实施例的共线链编组中连接参数示意图。

图4为本发明实施例的共线链编组中共线参数示意图。

具体实施方式

具体实施时，本发明可采用计算机软件技术实现自动运行流程。下面结合附图和实施例对本发明提供方法作进一步的说明。

本发明实施例的流程包括以下步骤：

步骤1，影像输入与分割；

读入原始遥感影像数据，利用meanshift算法影像分割。基于区域的影像分割用于道路提取的目的就是依据光谱特性将道路区域与周围地物分割，用于识别和探测道路特征。

meanshift算法为现有技术，影像的二维坐标和光谱通道组成多维特征空间。对于彩色影像，像素的三通道值对应3维色彩域r，像素行列号则对应二维空间域s，分割时核函数带宽分别为hs，hr。如果对于处理的当前点x，用x_i、z_i分别表示图像平滑前后的像素，则提供meanshit算法实现分割过程如下，以便实施参考：

(1)令变量y_i,j=x_i，其中i为像素号，j为迭代次数，从1开始；

(2)计算y_i,j+1=y_i,j+m_h(x)，直到算法收敛，记收敛后计算结果值为y_i,c，其中m_h(x)表示在空间域hs范围中的所有样本点相对于当前点x的偏移向量的平均值，它会指向高维球内样本点概率密度最大的方向，即梯度方向；

(3)赋值z_i=y_i,c；

(4)将空间域hs范围内，色彩域hr范围内的所有平滑后像素z_i编组为类{C_p}_p=1...m，也就是将收敛于一点的所有像素归为一类C_p，m为类的数目；

(5)标记分割区域，并将面积小于给定阈值的分割块合并，得到若干区域。

步骤2，道路区域形状指数分析；

道路的材料一致，表面平滑，在影像中具有近似的光谱特征。步骤1进行影像分割时将光谱相似的邻近像素标记为同一区域，所以道路与其他不同光谱的地物被分割为若干互不连通的区域，如植被、部分建筑物等。从这些分割后的区域中识别道路，需要利用道路的其他特性。道路的几何形状具有自身特殊性质，表现为狭长的、具有一定宽度的带状特征，利用形状指数辅助识别道路，可以将基于光谱信息和形状指数识别的道路区域赋予灰度值255，非道路区域灰度置零，生成道路的二值化影像。

实施例的实现方式为：

(1)面积S

城市中的道路供各类车辆行驶，连通不同的目的地，所以道路面积不会很小，利用面积阈值可以去除较小的干扰区域，如树木、汽车、房屋等。具体实施时，面积阈值可根据影像分辨率预先设定。

(2)长宽比R

道路狭长，包围道路区域的最小面积外接矩形的长宽比应该比较大。

R＝L_MAER/W_MAER

其中，L_MAER表示最小面积外接矩形的长，W_MAER表示最小面积外接矩形的宽。多边形最小外接矩形（Minimum Bounding Rectangle,MBR）是指以多边形所有顶点中最大及最小坐标确定的矩形，如图1中虚线矩形框。最小面积外接矩形（MAER）如图1、图2中点划线矩形框。对于具有一定倾斜角度的道路区域，MABR的长和宽比MBR更符合实际的区域形状。

(3)充满度F

F＝S/S_MAER

其中S为区域面积，S_MAER为最小面积外接矩形的面积。

充满度主要针对曲率较大的道路区域，或是纵横交错的道路连通区域。这类区域的形状较为复杂，外接矩形内有较大面积的空洞，如果只用长宽比和面积形状指标难以判断。如图1，道路连通区域MABR的长宽近似，但由于轮廓弯曲的形状导致区域实际面积远小于MABR的包围面积，针对此类型区域，利用充满度代替长宽比作为形状判断指标，充满度越小，道路区域可能性越大。实施例利用形状指数识别道路区域，满足如下公式的区域被判断为道路区域，并生成二值道路影像。

S＜T_S∩(F＜T_F∪R＞T_R)

其中，T_S为预设的面积阈值，T_F为预设的长宽比阈值，T_R为预设的充满度阈值。

步骤3，道路区域中线探测；

从二值化道路影像中探测道路中线。线特征探测法可分为两个层次：局部窗口分析和考虑光谱和上下文信息约束的模型法。局部窗口分析法常使用小影像区域，对局部线特征建模，依此模型查找线特征中心点；而第二种方法在第一种方法的基础上引入光谱和几何信息的约束，方法提高了探测的准确率，但增加了模型复杂度，降低了算法的效率。本发明从分类后的二值道路影像中探测道路，无需光谱等信息，只需从带状区域中获取中心点像素。实施例利用一种正交向排序法探测道路区域中心像素，并跟踪提取道路矢量基元。

现有技术中的基于正交方向排序的道路中线探测法假设道路中心点像素灰度值高于邻近像素，而算法的目的就是探测带状道路的中心区域像素。其核心思想是建立一维探测窗口，分别沿着水平和垂直方向在影像中移动探测，最终的道路中心像素为两个结果的并集。

建立三个窗口Ls、Le、Lb，其中Ls表示一维的探测窗口大小，Le表示统计窗口大小，Lb表示待判断像素pi的邻近范围，Lb<Le<Ls。算法步骤如下：

(1)输入二值道路影像，平滑处理；

(2)探测道路区域中心线像素：①将像素pi作为探测窗口的中心像素，将其与邻近范围内的Lb个像素都存入数组T_w(k)，其中k=0,1..Lb。对探测窗口内Ls个像素灰度值排序，将排在前面Le+1位的像素存入数组S_w(t)，t=0,1…Le；②比较T_w(k)和S_w(t)，如果满足T_w(k)∈S_w(t)，k=0,1..Lb，也就是说，像素pi与周围邻近像素是窗口内灰度最亮几个像素，则将像素pi判断为道路区域的中心线像素。探测窗口沿水平与垂直方向依次遍历，并取两个结果的并集；

(3)经过(2)中算法探测的道路中心线并非单像素，与中心像素邻近的像素同样被判断为中线点，针对这一问题，算法需要细化处理（可采用现有技术中Hilditch法）以获取单像素的道路中线；

(4)跟踪道路中线，获取矢量道路基元。

步骤4，道路基元编组。

矢量跟踪从道路区域中探测的道路中心线受多种因素的干扰，难以形成连续完整的矢量线。例如路面的阴影、遮蔽、路面的汽车等，导致道路线中断，甚至破碎。这些不连续的道路线是组成完整道路的基元，为获取完整道路线，则需要将基元编组。采用现有技术中的共线链编组法，可依据基元之间的共线参数及连接参数确定连接概率，将断开基元连接为连续的道路线。

(1)连接参数

道路基元之间的连接参数包括两线段之间的夹角α及断裂长度g，如图3，线段l1与l2之间的断裂长度为g，夹角α=α₁+α₂。其中，α₁为线段l1与l2之间的连线与线段l1的夹角，α₂为线段l1与l2之间的连线与线段l2的夹角。

(2)共线参数

共线参数指线段之间共线的指标参数θ1+θ2。如图4，θ1、θ2分别为线段l1与l2的中点（即基元中点）连线与l1、l2所形成的夹角。θ1+θ2越小，基元共线的可能性越大。此外，候选的道路基元越长，其为道路的可能性越大，越有可能与其他基元编组为完整道路线。所以实施例道路基元l1和l2之间的连接概率表示如下：

其中c₀、c₁分别表示连接参数与共线参数的权系数，c₀越大，c₁越小越考虑断裂项，c₀越小，c₁越大基元连接是对共线参数要求越高，两个系数可根据经验设定。Max_|θ1|+|θ2|为|θ1|和|θ2|之和的可能最大值，T_g为断裂间隔的阈值。l₁是线段l1的长度，l₂是线段l2的长度。

建立道路基元连接矩阵A，A为二维对称矩阵，矩阵元素为对应两个道路基元的连接概率，矩阵元素计算公式如下，

其中，

u的取值为0,1,…n，v的取值为0,1,…n，θu、θv分别为基元所在线段lu与lv的中点连线与lu、lv所形成的夹角，α是线段lu与lv之间的夹角，g是线段lu与lv之间的断裂长度。Max_|θu|+|θv|为|θu|和|θv|之和的可能最大值，T_g为断裂间隔的阈值。l_u是线段lu的长度，l_v是线段lv的长度。

其中基元个数为n+1。利用所有原始道路基元建立初始的道路基元连接矩阵，将连接概率，即矩阵元素大于给定阈值的道路基元编组连接，利用新生成的道路基元重新建立连接矩阵，计算矩阵元素，如此迭代直到没有可以连接的道路基元时结束编组过程。具体实施时，可设计为以下流程：

步骤4.1，根据步骤3所得道路基元建立初始的道路基元连接矩阵；

步骤4.2，将矩阵元素大于给定阈值的道路基元编组连接，若没有能够连接的道路基元结束编组过程，否则进入步骤4.3；

步骤4.3，利用步骤4.2编组连接后新生成的道路基元重新建立道路基元连接矩阵，返回步骤4.2。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (4)

1.一种遥感影像的道路中线自动提取方法，其特征在于，包括以下步骤：

步骤1，输入原始遥感影像数据，分割为若干区域；

步骤2，通过对步骤1所得各区域进行形状指数分析，判断道路区域，并生成二值道路影像；

步骤3，从二值化道路影像中探测道路中线，并跟踪道路中线，获取道路基元；

步骤4，通过道路基元编组获得自动提取结果，包括以下子步骤，

步骤4.1，根据步骤3所得道路基元建立初始的道路基元连接矩阵，

设道路基元连接矩阵记为A，A为二维对称矩阵，矩阵元素为对应两个道路基元的连接概率，矩阵元素计算公式如下，

其中，道路基元个数为n+1，u的取值为0,1,…n，v的取值为0,1,…n，θu、θv分别为道路基元所在线段lu与lv的中点连线与lu、lv所形成的夹角，α是线段lu与lv之间的夹角，g是线段lu与lv之间的断裂长度；Max_|θu|+|θv|为|θ_u|和|θ_v|之和的可能最大值，T_g为断裂间隔的阈值，l_u是线段lu的长度，l_v是线段lv的长度；

步骤4.2，将矩阵元素大于给定阈值的道路基元编组连接，若没有能够连接的道路基元结束编组过程，否则进入步骤4.3；

步骤4.3，利用步骤4.2编组连接后新生成的道路基元重新建立道路基元连接矩阵，返回步骤4.2。

2.根据权利要求1所述遥感影像的道路中线自动提取方法，其特征在于：步骤1中采用Meanshift算法进行分割。

3.根据权利要求1或2所述遥感影像的道路中线自动提取方法，其特征在于：步骤2中所述形状指数包括区域的面积S、区域的最小面积外接矩形的长宽比R、区域的充满度F=S/S_MAER，满足如下公式的区域被判断为道路区域，并生成二值道路影像，

S＜T_S∩(F＜T_F∪R＞T_R)

其中，T_S为预设的面积阈值，T_F为预设的长宽比阈值，T_R为预设的充满度阈值。

4.根据权利要求1或2所述遥感影像的道路中线自动提取方法，其特征在于：步骤3实现方式为，

建立三个窗口Ls、Le、Lb，其中Ls表示一维的探测窗口大小，Le表示统计窗口大小，Lb表示待判断像素pi的邻近范围，Lb<Le<Ls；执行如下步骤，

步骤3.1，输入二值道路影像，平滑处理；

步骤3.2，探测窗口沿水平与垂直方向依次遍历探测道路区域的中心线像素，并取两个结果的并集；探测方式如下，

设像素pi为探测窗口遍历到的中心像素，将像素pi与邻近范围内的Lb个像素都存入数组T_w(k)，其中k=0,1..Lb；对探测窗口内Ls个像素灰度值排序，将排在前面Le+1位的像素存入数组S_w(t)，t=0,1…Le；比较T_w(k)和S_w(t)，如果满足T_w(k)∈S_w(t)，k=0,1..Lb，则将像素pi判断为道路区域的中心线像素；

步骤3.3，对步骤3.2所得道路区域的中心线像素进行细化，获取单像素的道路中线；

步骤3.4，跟踪步骤3.3所得道路中线，获取道路基元的矢量。

Images