结合路径形态学的高分遥感影像道路提取方法

戴激光，苗志鹏，葛连茂，王晓桐，朱婷婷

(1.东华理工大学 江西省数字国土重点实验室，南昌 330013；2.辽宁工程技术大学 测绘与地理科学学院，辽宁 阜新 123000)

0 引言

利用遥感影像进行道路信息采集在交通管理、城市规划、自动车辆导航、地理信息系统数据库更新、制作电子地图等方面有重要应用价值，针对该领域的研究一直是国内外学者研究的热点。随着空间分辨率的提高，高分辨率影像可对道路几何光谱特征进行更加精细地刻画；但另一方面，冗余的细节信息诸如房屋、隔离带、阴影、车辆等因素均会对道路提取构成巨大的干扰，进而大幅增加了道路提取的难度。

Vosselman等[1]对道路几何与纹理信息进行了描述，道路路面存在一定的纹理相似性，并且道路具有长度较长的几何特性。作为一种形态学方法，路径形态学可提取图像中长而窄的结构，这一特点亦符合道路的特征，故很多学者将路径形态学方法应用于道路提取领域[2]。例如Talboth等人[3]提出不完整的路径开运算，能越过一些“缺失的像素”搜索路径，可应用于阴影遮挡、纹理变化较大或破损道路的提取，但该方法的缺点是效率低并且误提取率较高。刘小丹等人[4]提出一种将Hough变换与路径形态学结合的道路提取方法，通过Hough变换构造合适的结构元，再利用路径形态学提取道路。但是Hough变换产生的结构元方向单一，难以提取曲率较大的道路。Appleton等[5]在Heijmans等[2]方法基础上设计了一种有效的分解方法，将空间复杂度降低到对数级。同时，设定了4种三邻接结构元，可提取曲率较大的道路。Schubert等[6]改进了Appleton等[5]提出的方法结构，增加4种对角线方向的结构元，可提取更多方向的道路。王双等[7]设计了2种四邻接结构元，可提取曲率更大的道路，但同时也会连接较多的非道路区域。

路径开运算应用于道路提取工作中，虽然取得了较好的实验成果，但由于该方法仅考虑道路纹理和长度特征，并未顾及道路宽度一致性和连通性，因而提取结果中存在道路与非道路区域粘连和断裂问题。通过深入研究后发现引发这些问题的主要因素是：(1)道路两侧界限不清晰。当道路与其周围地物间存在“异物同谱”现象时，路径开运算便将其视为一个整体，导致道路与非道路区域粘连问题。(2)路面纹理信息不均一。路面上的车辆、人行横道、建筑物或树木阴影会造成局部纹理信息发生突变，这使得路面纹理一致性的特点不再突出，直接影响影像分割结果，从而易引发道路断裂。有鉴于此，针对上述问题，在传统路径形态学的基础上，本文提出了一种改进道路提取方法。

1 原理

1.1 数学形态学

数学形态学(mathematical morphology)是建立在集合论、拓扑论、随机集等众多学科的理论基础之上描述形状结构的科学。其基本思想是将携带形状、大小、灰度等信息的结构元作为“探针”遍历图像，从图像中“探测”出相应信息。根据输入图像的不同，数学形态学可以分为二值形态学和灰度形态学，两者都包含膨胀、腐蚀、开运算、闭运算、顶帽运算、黑帽运算、击中击不中变换等运算。

由于路面上存在路面油渍、油返砂等问题，这使得道路内部光谱相似性降低，并将直接影响道路提取结果，所以本文使用灰度形态学对影像进行预处理，以减弱或消除这些干扰。根据上述对路面问题的分析，对于深色路面上的浅色干扰信息，本文采用腐蚀或开运算处理；对于浅色路面上的深色干扰信息，本文利用膨胀或闭运算处理。

1.2 Otsu双阈值分割

当深色路面存在阴影等深色干扰时，路面灰度值的范围介于浅色地物和深色干扰(阴影等)的灰度值之间，因此需用区间阈值才能获得道路区域，故本文采用Otsu双阈值分割方法对影像进行分割。利用Otsu方法得到2个阈值，将灰度值处于两阈值之间的像素视为道路(灰度值为255)，其余视为非道路(灰度值为0)，从而得到二值图像。

Otsu双阈值分割是基于灰度直方图得到各分割特性值发生的概率，用两个阈值变量T1、T2将影像分为三类(C0、C1和C2)，其灰度值范围分别是[0，T1]、[T1+1，T2]、[T2+1，m]，然后求出每一类的类内方差及类间方差，选取类间方差最大或类内方差最小的T1、T2作为最佳阈值[8]。其公式如下：

(1)

1.3 路径开运算

路径开运算(path opening)是路径形态学的一种基本算法，属于形态学的范畴。形态学中的开运算利用结构元遍历图像，可以剔除小于结构元大小的像素块。而路径开运算则利用结构元设定连接规则，将灰度值相近且邻接的像素点连接成长度不同的路径，然后剔除长度小于阈值L路径点。定义如下：在二值图像E中定义一种连接关系：ab，表明由a到b存在一条路径。其中，a是b的后继节点，b是a的前驱节点，符号“”表明路径方向由a到b。将这种连接关系和图像E中所有像素点Ω结合在一起(Ω，)，即为邻接图(directed acyclic graphs，DAG)。邻接图可以清晰地表示像素点间的连接关系。图1(a)中，黑色部分为一种结构元定义的连接关系(方向为S-N)，灰色部分为该结构元构成的邻接图。其中，任一像素点F0最多拥有S1、S2、S3 3个后继节点，任一点S0最多拥有F1、F2、F3 3个前驱节点。

定义1将每个点x关于连接关系“”的后继节点的集合，记为δ({x})，即δ({x})={y∈Ω|yx}；x关于连接关系的前驱节点的集合，记为即y}。设X为前景图像，对于∀X⊆Ω，有：

δ(X)={y∈Ω|xyfor somex∈X}

(2)

(3)

定义2设Ta=(a1，a2，…，ak)为a的k元组，那么当且仅当aiai+1，且i∈[1，k- 1]时称Ta为a长度为k的路径。

定义3设Пk为所有长度为k的路径集合，设路径Ta里的所有像素点的集合为σ(Ta)。路径开运算的定义为：

Ak(X)=∪{σ(Ta)|Ta∈Пkandσ(Ta)⊆X}

(4)

Ak给出前景图像X中路径长度大于或等于k的像素点集合，即为路径开运算结果。图1为路径开运算示意图，图1(b)是一些前景黑色像素点在图1(a)给出连接关系下的所有路径。其中，每个点旁边的数字是该点所在路径的长度值(在一种连接关系下，若一个点同时存在于多条路径上，规定其属于最长的路径)。图1(c)为长度阈值L为6时路径开运算结果，对比图1(b)可以看出，路径长度大于或等于6的像素点被保留下来，其余像素点被去除。当L为7时，仅保留了路径长度大于或等于7的像素点(图1(d))。

图1 路径开运算原理

基于以上分析，结合影像中道路区域与非道路区域的特点，发现道路区域的路径长度一般长于非道路区域。因此将路径开运算应用于道路提取中，通过人工设定合适的路径长度阈值L，就可以有效地剔除非道路区域，保留道路区域。

本文的路径开运算方法来自Schubert等方法，该方法包含8种三邻接结构元(图2)。其中每个点有3个连接方向，表示在搜索路径时有角度方向上的容差，可提取具有一定曲率的道路。每种结构元可以探测一定方向的道路，不同的探测结果进行组合即可获得完整道路网(图3(c))。但由于前文所述的原因，有些道路提取结果不尽人意。如图3(d)所示，影像中部存在的树木阴影导致了路径开运算结果中的道路断裂问题；影像中的房屋与道路灰度近似并且邻接，导致结果中存在道路与非道路区域的粘连问题(图3(f))。

图2 邻接图

图3 路径开运算提取道路

1.4 非道路粘连削减方法

针对上述路径开运算中存在的问题，本文提出一种非道路粘连削减方法。道路与非道路粘连处(以下简称粘连处)可分为3个区域：道路、非道路、连接区域。连接区域即道路与非道路衔接区域。如图4(a)所示，可以看到粘连处存在以下2个特征：(1) 连接区域宽度较小。(2) 非道路区域的路径长度小于道路区域的路径长度。

基于以上特征，本文首先设定大小为3×3的矩形结构元，通过形态学开运算将连接区域断开，得到孤立的非道路区域；其次设定路径长度阈值L进行路径开运算。该阈值可依据实际影像区域道路目视分析，由人工进行设定；最后使用闭运算将道路网复原。

该方法有一个重要的前提：连接区域宽度应当很小，这样形态学开运算才能在不影响道路区域的前提下断开连接区域。若连接区域很宽，那么此处粘连无法去除。图4展示了一条道路线上的两处粘连问题的处理结果。图4(a)中矩形框内的粘连，其连接区域宽度较小，本文方法可以将其去除。而圆圈内的粘连由于其连接区域宽度较大，因而无法进行删除。

图4 削减粘连示意图

1.5 断裂端点判读

由于路径开运算只能探测相邻区域灰度值近似的像素，无法跨越灰度值相差很大的像素，故仅依靠形态学方法难以解决道路断裂问题。由于模板跟踪在道路提取过程中对道路光谱变化有一定的鲁棒性[10]。因此本文在路径开运算的结果中，取其断裂处的端点作为种子点，并由此提出改进的圆形模板跟踪方法，以此完成断裂处道路中心线的跟踪，但在开展该工作前首先需要确定断裂处端点。

考虑到路径开运算提取的道路结果中存在很多端点，因此在道路跟踪前首先需要判断哪2个端点属于同一断裂区域。仔细观察道路断裂处存在下述2个特征：(1) 道路断裂是由阴影、人行横道等造成的，因而道路断裂处长度较短。因此可依据这一特征，比较所有端点的间距，找出间距较小的端点对。(2) 断裂处曲率变化较小。即断裂处两侧线段的角度差异应当较小。依据上述特征，可采用最小二乘法将两端点所在区域进行拟合，提取2条最佳拟合线段，同时将2个端点连接起来形成连接线段，分析3条线段的角度差异是否低于π/4。若2个端点同时满足上述2个特征，则认为这2个端点应属于同一断裂区域。

1.6 改进的圆形模板匹配跟踪方法

本文提出的跟踪方法是对文献[10]方法的改进。文献[10]方法的思想是：人工选取若干种子点，在两两种子点间以一定规则内插新种子点，不断重复内插过程直至达到阈值。该方法用直线段逼近曲线，对曲率较小的道路提取效果较好。其优点在于，若前几次内插的种子点精度较高，随着种子点的不断增多，其形状会不断逼近道路形状，其计算过程会越来越快，精度会越来越高。其缺点在于：(1) 若初始种子点间道路曲率过大，第一次内插过程的难度就会比较高。(2) 已有内插点的精度会影响后续内插点精度，进而影响道路提取质量。

基于上述问题，本文提出改进的圆形模板提取方法。与原方法主要不同之处是：(1) 改进了内插种子点的方式，将其对精度影响均匀分布到每个点上。(2) 加入步长增加措施，可越过“障碍”继续搜索道路。具体流程如下：

(1)采用L0滤波方法[11]处理输入影像，利用3×3矩形结构元计算每个点的形态学梯度。

(2)将上文得到的断裂处两端点作为初始种子点。以一个初始种子点及其八邻域内的点为圆心，分别建立半径为1像素的圆形模板，计算模板内点的梯度之和Gradsum。将同一半径下Gradsum最小的模板视为最优模板，并判断其Gradsum是否超过阈值δ(经验阈值，本文设为200)，若不超过δ，则以1个像素的步长逐步扩大模板半径。继续该搜索过程直至最优模板Gradsum超过δ，将前一半径的最优模板设为标准模板，其圆心p1为道路中心点，其半径r为道路宽度的一半(图5(a)过程①)。按照上述流程，在第二个初始种子点处寻找半径为r的标准模板，获得道路中心点p2(如图5(a)过程④)。将p1到p2的方向设为初始跟踪方向。

(3)若p1到p2的距离小于r，连接上述两点并停止跟踪。否则从初始端点p1开始，在跟踪方向上长度为r的点处，垂直于跟踪方向上建立一系列标准模板(图5(a)中过程②的Ti)。定义模板角度θmod为当前模板圆心P与前一个模板圆心M连线的角度(图5(a)中的θ)。

定义灰度均值差为模板P中每个点q的灰度值与前一个模板M灰度均值之差的绝对值之和，公式如下所示：

(5)

定义梯度方差为模板P中每个点q的梯度Gq与模板P梯度均值之差的平方之和的平方根，公式如下所示：

(6)

定义模板角度差为模板P的θmod与其前6个模板θmod均值(不足6个按实际数量计算)之差的绝对值，公式如下所示：

(7)

分别计算Ti的灰度均值差、梯度方差、以及模板角度差，将3个量归一化后加权求和(公式(8))，求得权重函数取得最大值时的圆心坐标，将其作为待定道路中心点。若2个模板的权重函数值相同，取模板角度差较小的一个；

(4)为验证待定道路中心点是否是道路中心，本文采用边缘判断法对其进行判断。由于路面与道路两侧存在灰度差异，所以在待定道路中心点模板A两侧，垂直于跟踪方向上建立标准模板B、C(图5(b))，分别计算B与A、C与A的灰度均值之差的绝对值。若任一差值大于阈值Dgary(本文设为30)，说明该点为道路中心点，并将其作为新的初始端点，其模板角度为跟踪方向，重复步骤(3)；若差值均小于阈值Dgary，说明该步长下求得的待定点存在异常(图5(a)过程③中障碍物)，需要进一步判断，进入步骤(5)；

(5)将步长变为二倍，重复步骤(3)、步骤(4)。若二倍步长无法满足要求，则变为三倍步长(图5(a)过程③)。若三倍步长也无法满足要求，停止跟踪。

图5 道路边缘判断、断裂处跟踪示意图

(8)

1.7 方法步骤

通过上述对不同方法原理的分析，本文给出了具体的方法处理流程(图6)。

图6 技术流程图

(1) 采用灰度形态学对输入高分辨率遥感影像进行处理，以消除路面噪声；

(2) 利用Otsu方法分割影像，得到包含道路区域和非道路区域的二值图；

(3) 采用路径开运算消除大部分的非道路区域；

(4) 设定道路长度最低阈值L，利用形态学和路径开运算削减粘连的非道路区域；

(5) 对断裂端点进行判断，若存在则转入步骤(6)，否则转入步骤(7)；

(6) 在L0滤波影像基础上，以断裂端点为种子点，运用改进圆形模板匹配方法跟踪道路区域，进入步骤(7)；

(7) 输出道路结果。

2 实验结果分析和算法评价

2.1 实验影像及实验方法介绍

本文基于VS2013平台，采用C++编程语言对3幅遥感影像进行实验。其中，Pleiades影像大小为632像素×679像素，空间分辨率为0.5 m；2幅高分二号影像大小均为1 000像素×1 000像素，空间分辨率为0.8 m。本文首先分别利用TALBOT 等方法中“不完整路径开运算”、Appleton等方法中“有效的路径开运算和闭运算”以及Schubert等方法中“有效的灰度路径开运算”进行路径开运算实验，阐述断裂问题和粘连问题。然后在Schubert等方法的结果中，加入本文的后处理方法，以验证本文方法的有效性。其中，TALBOT 等方法提供了该方法的C++代码，Schubert等方法提供了Schubert等方法和Appleton等方法的C++代码，因而对比实验结果是可信的。

2.2 实验一

实验一选取Pleiades拍摄的农村遥感影像。该影像中道路与居民区有一定的粘连，路面有树木阴影遮挡；一些居民区呈带状分布，灰度值与道路相近，形成了类似道路的结构，对于道路识别构成很大的干扰。

基于如图7(b)所示二值图，3种方法开展的道路提取实验，结果分别如图7(c)、图7(d)、图7(e)所示。可以看到，由于当前影像覆盖区域比较复杂，TALBOT 等方法虽然将道路转弯处提取出来，也解决了树木阴影遮挡对道路遮挡的问题，但仍然存在很多粘连，右侧道路误提取率较大；而Appleton等方法在道路转弯区域提取效果不佳，并且将成排的房屋视为道路；而Schubert等方法提取道路结果精度较高，不仅能够提取出弯曲的道路并且误提取率低，但依然存在非道路与道路粘连问题和道路断裂问题。因此本文以Schubert等方法为基础，实验证明是可行的。如图7(f)所示，依据本文提出的削减粘连方法，经过一系列形态学与路径开运算相结合的操作(图8)，削减了大部分的粘连区域，仍有一些粘连由于连接处较宽而无法削减，并且产生了一个新的道路断裂(图8(c)左下)。如图7(g)断裂处放大显示，改进圆形模板跟踪的道路中心线与实际道路中心线基本吻合。最后用道路中心线将断裂处连接起来即可得到该影像的道路网(图7(h))。

图7 Pleiades影像实验

图8 削减粘连

2.3 实验二

实验二选取高分二号城区遥感影像。高大建筑物产生的阴影将影像上的路面区域截为两段(图9(a)下部)，因而道路内部光谱信息呈现出不一致的现象。同时由于道路上存在大量的车辆，以及较多的人行横道信息，这均会导致二值图中出现孔洞。并且右下角横向道路两侧存在与道路灰度相近的植物以及部分低矮建筑的阴影，这将会造成道路与非道路的粘连。

图9 高分二号影像实验

采用灰度形态学对图像预处理，采用Otsu方法对影像进行二值化处理，采用3种方法进行路径开运算，其结果分别如图9(c)、图9(d)、图9(e)所示。由于该影像道路结构比较简单，3种路径开运算方法提取道路的完整性很高，但是提取质量各不相同。其中，TALBOT 等方法在能够较好地提取部分道路，但误提取率很高。Appleton等方法由于其结构元比较适合提取直线，其误提取率最低，但完整性略低(图9(d)下部道路缺失)。Schubert等方法误提取率介于其他2种方法之间，完整性较好。对Schubert等方法的结果进行削减粘连后，发现路面仍存在空洞，用二值形态学闭运算对空洞进行填补，并通过圆形模板匹配跟踪方法对道路断裂处进行连接，得到道路提取结果(图9(f))。

2.4 实验三

实验三选取高分二号城市郊区遥感影像。对影像的目视分析可以发现，场景中存在的建筑物阴影会造成道路线断裂，一些道路两侧界线不清晰可能会造成道路与非道路的粘连。

采用Otsu方法对影像进行二值化处理后，如图10(b)可以发现道路边缘比较规整，横向部分的道路虽然曲率较小，但断裂情况依然非常明显。3种路径开运算提取结果分别如图10(c)、图10(d)、图10(e)所示。其中，TALBOT 等方法虽对断裂问题有所改善，但误提取较高并且粘连问题较严重；Appleton等方法和Schubert等方法完整性相似，但Schubert等方法的粘连较少。如图10(f)所示，针对道路三岔口处的断裂问题，利用本文的方法，不仅能够道路粘连处进行削减，同时能够解决由阴影导致的道路断裂问题，并如图10(f)所示由此得到完整的道路提取结果。

图10 高分二号影像实验

2.5 算法的精度和效率评价

本文利用文献[12]提出的完整性(complete percent，CP)、正确率(correct percent，CR)、提取质量(quality percentage，QL)3个指标对实验结果进行量化分析评判，3个指标的计算方法如下所示。

(9)

(10)

(11)

式中：TP表示匹配的道路长度；FP表示未匹配的提取道路的长度；TN表示匹配的参考实际道路长度；FN表示未匹配的实际道路长度。本文通过统计结果中各部分对应的像素的个数，根据公式(9)、公式(10)、公式(11)，得到各类方法的精度统计如表1所示。

表1 不同方法精度统计表

由表1可以看出，Appleton等方法由于其采用的结构元数量较少(4个)，因而运行速度较快，但提取的完整性不高。Schubert等方法和TALBOT 等方法包含较多方向的结构元(8个)，因而提取道路的完整性较高，效率较低。另外，TALBOT 等方法还包含对道路断裂问题的改善，因此运行速度最低，完整度最高但出现很多的误提取问题。本文所用的Schubert等方法相较于Appleton等方法，优势在于可提取弯曲的道路并且道路提取完整性较好；相较于TALBOT 等方法，优势在于误提取率较低、运行速度快。而本文方法在Schubert等方法基础上进行改进，后处理结果的完整性、正确率和提取质量均表现良好。本文方法主要涉及3个步骤，分别是路径开运算、削减粘连和连接断裂。其中路径开运算方法效率与二值图像的像素数和结构元的数量有关。削减粘连过程效率取决于形态学运算和路径开运算的方法效率。连接断裂方法效率与断裂区域的灰度稳定程度有关。从整体来看，本文提出的改进步骤速度均控制在1 s之内，效率较高。

3 结束语

本文首先对基于路径开运算的高分辨率道路提取结果存在的问题进行分析探讨；然后提出削减道路粘连的处理策略以及解决道路断裂问题的改进圆形模板跟踪方法；最后选用3幅有典型意义的影像进行实验。实验结果表明本文方法对于解决道路粘连问题和道路断裂问题有较好的效果。但本文方法也有不足之处，即需要人工设定路径开运算阈值L，这也是路径开运算难以避免的问题。因此如何改进路径开运算方法，自适应设定阈值L，这将是未来研究的重点。