CN110543885A - 一种交互式提取高分辨率遥感影像道路并生成路网的方法 - Google Patents

Abstract

为了解决高分辨率影像中道路半自动提取方法无法进行流程化快速提取并构建路网的问题，通过道路交互式提取、矢量化处理以及生成路网全流程处理，可以有效的提升道路提取效率。本发明公开了一种新的交互式提取高分辨率遥感影像道路并生成路网的方法，提出了利用缓冲区技术限制算法工作区域、利用形态学算法优化区域增长提取结果、利用矢量化跟踪算法、简化算法以及端点修饰算法对路段矢量进行优化、利用单条道路连接算法以及道路交叉口连接算法对阴影、路口情况进行合理连接和利用缓冲区技术生成路网的策略，实现了道路半自动交互式提取并生成路网全流程处理。

Description

一种交互式提取高分辨率遥感影像道路并生成路网的方法

技术领域

本方法属于图像处理与目标识别技术领域，特别是涉及一种新的流程化道路半自动提取并获取路网的方法。

背景技术

随着城市和乡村建设的加快，如何快速识别并提取道路、更新路网就成了至关重要的问题。道路作为城市交通的主要组成部分，在政治、经济、军事等领域发挥着重大作用。而目前，随着大规模高空间分辨率遥感器的投入使用，遥感影像的空间分辨率甚至可以达到亚米级。与道路在低分辨率影像中呈现的细线形状不同，高分辨率影像中道路是连续的同质区域，这意味着我们可以更加准确地提取道路。但由于“同谱异物，同物异谱”、影像分辨率不同、道路类型不同、道路遮挡、断裂、光照条件等情况的影响([文献1])，设计分类以及提取道路算法的难度也随之加大。

由于道路提取需要大量先验知识、环节复杂，全自动道路提取方法并不能取得很好的效果([文献2-4])。因而，在实际生产中，大多数企业仍旧通过传统人工勾画的方式提取路网。人工操作虽能保证拓扑关系建立准确，但工作量大、效率低，故结合计算机目视解译的半自动道路提取方法仍是当前研究的重点([文献5])。

根据提取算法的流程以及侧重点的不同，现有的半自动道路提取方法可分为以下几类：基于区域增长的方法([文献6-7])、基于动态规划的方法([文献8-9])、基于边缘特征分割的方法(如梯度计算边缘的方法([文献10])，边缘细化分割方法([文献11]))、基于图像分割的方法([文献2、4、12-13])、基于模板匹配的方法([文献14-15])、基于Snake等主动轮廓模型的方法([文献16-18])、基于机器学习以及神经网络的方法([文献2、10、19])。其中基于动态规划的方法需要手动选取多个点，耗时较长([文献8])。基于Snake等主动轮廓模型的方法近年来虽研究较多，且对传统Snake模型进行了较大改进，但其最大的缺点还是难以适应高分辨率遥感影像的复杂情况，需要人工设定参数，检测效率低，鲁棒性差([文献16])。基于区域增长、边缘以及图像分割的方法，容易因部分道路边缘模糊，将道路与周围地物混淆，尤其容易受到如遮挡、阴影、裸地等影响，出现错分割。而基于模板匹配以及神经网络的算法，受模板以及标记的样本影响大，且由于道路情况十分复杂，模型针对不同道路的适应能力不强([文献14])。

同时，上述方法大多都将重点放在影像栅格处理层面，对于中心线提取、矢量化以及路网生成讨论较少。而道路矢量对地理信息系统的构建以及信息综合分析有至关重要的作用，需要得到重视。并且，上述方法都无法攻克不同分辨率、不同宽度、遮挡、阴影等道路提取难点。针对以上问题，尤其是针对当前道路提取方法适应能力差、道路提取关键难题(如遮挡、阴影等)无法解决、矢量化与构建路网难的问题，本发明提出一种新的交互式道路提取以及自动生成路网的方法。首先由用户确定道路等级并通过交互给定相应等级道路在影像中的宽度，再给定种子点依次提取道路，最后通过路段以及交叉路口连接，生成正确路网。

该发明中，为解决不同道路宽度的问题，提出了分等级提取道路的思想，不同等级道路宽度不同，在提取道路前首先获得某一等级道路宽度，在这一等级道路提取结束后，再依次提取其他等级道路；为解决不同分辨率提取难的问题，本发明将影响提取效果的关键参数，如矢量化过程中删除多余线段的长度阈值、连接长度阈值、简化线参数等都与影像分辨率以及道路宽度联系起来，确保该方法在不同分辨率下都具有有效性；为解决阴影、遮挡等导致的道路断裂问题，采用了形态学方法以及连接算法填充不同尺度的断裂，保证单条道路连接正确；为解决地物光谱类型相似无法正确区分道路与其它地物的问题，利用种子点连线根据设置的道路宽度按一定准则生成缓冲区作为掩膜来确定道路提取算法的工作区域，防止周围地物因过于相似使提取效果变差；为解决路网交叉口连接问题，设计了道路交叉口连接算法，完成了对“T”型路口、“Y”型路口、“十”型路口的连接。通过实验分析，该系列策略对路网提取效果有很大的改善，可较好地用于高分辨遥感影像道路提取与路网更新。

参考文献：

[1]Zhao J Q,Yang J,Li P X,et al.Semi-automatic Road Extraction fromSAR images using EKF and PF[J].ISPRS-International Archives of thePhotogrammetry,Remote Sensing and Spatial Information Sciences,2015,XL-7/W4:227-230.

[2]Bakhtiari H R R,Abdollahi A,Rezaeian H.Semi automatic roadextraction from digital images[J].The Egyptian Journal of Remote Sensing andSpace Science,2017:S1110982317300820.

[3]Miao Z,Wang B,Shi W,et al.A Semi-Automatic Method for RoadCenterline Extraction From VHR Images[J].IEEE Geoscience and Remote SensingLetters,2014,11(11):1856-1860.

[4]Nunes D M,Medeiros,Nilcilene das Santos A D P D.Semi-Automatic Road Network Extraction from Digital Images Using Object-BasedClassification and Morphological Operators[J].Boletim de Ciências Geodésicas,2018,24(4):485-502.

[5]周绍光.高分辨率遥感影像中道路提取方法的研究[M].北京：科学出版社，2012：3-10.

[6]余洁,余峰,张晶,等.结合区域生长与道路基元的高分辨率遥感影像道路提取[J].武汉大学学报(信息科学版),2013,38(7):761-764.

[7]李建飞,文志强,胡永祥,et al.基于改进区域生长的遥感影像道路提取[J].计算机工程与应用,2016,52(2):209-213.

[8]曹帆之.利用动态规划半自动提取高分辨率遥感影像道路中心线[J].测绘科学技术学报.2015:1673-6338(2015)06-0615-04

[9]Gruen A,Li H.Road extraction from aerial and satellite images bydynamic programming[J].Isprs Journal of Photogrammetry&Remote Sensing,1995,50(4):11-20.

[10]Hormese J,Saravanan C.Automated Road Extraction From HighResolution Satellite Images[J].Procedia Technology,2016,24:1460-1467.

[11]Xiao Y,Tan T S,Tay S C.Utilizing edge to extract roads in HIGH-RESOLUTION SATELLITE IMAGERY[C].IEEE International Conference on ImageProcessing.IEEE,2005.

[12]Chaudhuri D,Kushwaha N K,Samal A.Semi-Automated Road DetectionFrom High Resolution Satellite Images by Directional MorphologicalEnhancement and Segmentation Techniques[J].IEEE Journal of Selected Topics inApplied Earth Observations and Remote Sensing,2012,5(5):1538-1544.

[13]Jianhua Wang,Qiming Qin,Zhongling Gao.A New Approach to UrbanRoad Extraction Using High-Resolution Aerial Image[J].International Journalof Geo-Information.Geo-Inf.2016,5,114；doi:10.3390/ijgi5070114

[14]陈光,眭海刚,涂继辉,et al.基于P-N学习的高分遥感影像道路半自动提取方法[J].武汉大学学报(信息科学版),2017(6).

[15]阳春花,李芹,李春燕.面向地理国情监测的道路信息半自动提取[J].地理空间信息,2018,v.16；No.105(05):9+82-84.

[16]王峰萍,王卫星,薛柏玉,et al.GVF Snake与显著特征相结合的高分辨率遥感图像道路提取[J].测绘学报,2017,46(12).

[17]刘少创,林宗坚.航空遥感影像中道路的半自动提取[J].武汉大学学报·信息科学版,1996,21(3):258-264.

[18]Abdelfattah R,Chokmani K.A semi automatic off-roads and trailsextraction method from Sentinel-1data[C].IEEE 2017IEEE InternationalGeoscience and Remote Sensing Symposium(IGARSS)-Fort Worth,TX(2017.7.23-2017.7.28)

[19]Yanan Wei,Zulin Wang,and Mai Xu.Road Structure Refined CNN forRoad Extraction in Aerial Image[J].IEEE GEOSCIENCE AND REMOTE SENSINGLETTERS,VOL.14,NO.5,MAY 2017

发明内容

针对目前高分辨率遥感影像道路提取无法满足流程化提取路网的要求，且无法满足已有的半自动交互式提取方法的问题，本发明提出一种新的交互式半自动高分辨率遥感影像道路提取并生成路网的方法，从而实现快速提取道路并更新路网。

本发明所采用的技术方案是：一种交互式提取高分辨率遥感影像道路并生成路网的方法，其特征在于，根据本发明的技术流程和所涉及的主要技术环节，如图1所示，包括用户交互、栅格影像处理、单条道路矢量化处理以及路网生成四大步骤，具体为：

步骤1，通过用户交互确定当前等级道路宽度以及种子点，同时确定算法工作区域，即缓冲区；

步骤2，栅格影像处理，包括如下子步骤；

步骤2.1，在步骤1中的缓冲区内进行区域增长，获得道路初始提取结果；

步骤2.2，利用形态学闭运算填充空洞；

步骤2.3，利用形态学细化操作获取道路中心线；

步骤3，单条道路矢量化处理，包括如下子步骤；

步骤3.1，利用矢量化处理获得单条道路路段点集；

步骤3.2，首先对单条道路路段点集进行预处理，然后利用道格拉斯-普克算法对点集进行简化；

步骤3.3，利用端点删除算法对路段进行矢量修饰；

步骤3.4，利用单条道路连接算法对单条道路进行连接；

步骤4，通过交叉路口连接算法将路口处连接，并利用缓冲区技术结合步骤1获取的道路宽度生成路网。

进一步的，步骤2.2中采用了Zhang-Suen并行算法进行细化处理，具体实现方式如下，

(2.2.a)如若当前中心点p₁满足以下四个要求，p₁的像素值为255，则删除，即将其灰度值改为0：

2≤N(p₁)≤6 (2)

T(p₁)＝1 (3)

p₂·p₄·p₆＝0 (4)

p₄·p₆·p₈＝0 (5)

其中，N(p₁)表示p₁八邻域的非零点数量，p_i(i＝2，3，...，9)表示为八邻域中第i个像素，T(p₁)表示八邻域内由0到255的跃迁次数；

(2.2.b)保留(2)(3)条件不变，添加以下条件，若满足，则删除：

p₂·p₄·p₈＝0 (6)

p₂·p₆·p₈＝0 (7)

重复上述两步，直至没有符合上述要求的像素为止。

进一步的，步骤3.1的具体实现方式如下，

(3.1.a)新建一个储存点集的双向队列，从上到下、从左到右，全局搜索灰度值为255的初始点，记为FirstPoint，令当前点CurrentPoint＝FirstPoint，将此点坐标储存进队列中，并将灰度值置为0，进入(3.1.b)，若全局未搜索到灰度值为255的像素，则全部道路已被跟踪，跳出循环并结束算法；

(3.1.b)在八邻域内搜索下一个点，其中八邻域搜索顺序为：以p₁为中心点，设1-8为搜索顺序，其中1位于p₁的左上方，2位于p₁的正上方，3位于p₁的正左方，4位于p₁的右上方，5位于p₁的左下方，6位于p₁的正右方，7位于p₁的正下方，8位于p₁的右下方，通过规定先搜索上邻域(1-4)，再搜索下邻域(5-8)，使搜索具有方向性，达到减少搜索次数的目的；若搜索到灰度为255的点，则将此点记为CurrentPoint，将此点坐标储存进双向队列尾部，并将灰度值置为0，重复(3.1.b)；直到当前像素八邻域未搜索到灰度为255的点，则进入(3.1.c)。

(3.1.c)由于搜索的初始点并非是端点，此时令CurrentPoint＝FirstPoint，进入步骤(3.1.d)，开始跟踪FirstPoint另一个方向的道路；

(3.1.d)依照八邻域搜索顺序搜索当前点的八邻域，若搜索到灰度为255的点，则将此记为当前点，将此点加入双向队列头部，并将灰度值置为0；重复(3.1.d)，直到当前像素八邻域未搜索到灰度为255的点，则完成一条道路的跟踪，即获得单条道路点集，然后再进入(3.1.a)。

进一步的，步骤3.2的具体实现方式如下，

首先对单条道路路段点集进行预处理，具体为：对跟踪得到的多个队列，通过长度阈值筛选将部分队列删除，长度阈值T_delete定义为：

T_delete＝(50.0*resolution)*L/3.0 (8)

其中resolution指影像分辨率，若影像无地理编码，则将resolution赋值为1即可，L是用户交互中提取的道路宽度；

然后利用道格拉斯-普克算法对点集进行简化，具体实现如下：

(3.2.a)连接两端点A和B并记为保留点，计算两端点间各点到两端点连线的欧式距离，并求出距离最远的点P，若两端点间无点，则所有点已简化完毕，退出循环；

(3.2.b)比较P至直线的距离与阈值的大小，若大于阈值，则保留P，否则删去当前点；

(3.2.c)以A、P为端点重复(3.2.a)-(3.2.c)；同时，以P、B为端点重复(3.2.a)-(3.2.c)。

进一步的，步骤(3.2.b)中阈值T定义为：

T＝3.0*resolution (9)

其中resolution是指影像分辨率，3.0是在影像无地理编码或者分辨率为1米时得到的最优阈值。

进一步的，步骤3.3的具体实现方式如下，

对每条线段，取与其端点C相邻的两点A与B，求向量与的夹角θ，当θ＜10°时，保留C点，否则删除C点。

进一步的，步骤3.4中单条道路连接的连接规则由矢量端点间距离规则以及矢量间夹角规则共同组成，若两矢量路段满足连接规则，则将两条矢量路段连接，并将属性ID修改为同一个值；同时，建立ID查找表，记录每条矢量线原始ID与修改后ID，方便追踪每条路段的修改，保证了多条路段搜索连接正确；

其中连接规则中，矢量间端点距离规则由两部分构成；首先，矢量间端点距离必须为全局最小距离；其次，规定阈值d_t，若两端点间距离小于d_t，则可以连接，如公式10，阈值d_t的定义为，

d_t＝(50.0*resolution)*L/3.0 (10)

其中resolution指影像分辨率，若影像无地理编码，则将resolution赋值为1即可，L是用户交互中提取的道路宽度；

矢量间夹角规则定义如下，若已找到满足距离规则的矢量路段端点A与B，则寻找与端点距离最近的节点A’与B’，求出A’A与B’B的方位角θ₁与θ₂(θ₁，θ₂∈(0°，180°])，计算θ＝|θ₁-θ₂|，若θ＞90°，即满足夹角规则；

单条道路连接的具体步骤如下，

(3.4.a)建立端点容器Ver，将提取到的道路端点坐标依次压入Ver中；假设端点个数为n，端点编号为i(0，1，2，3，...，m)，此时，当端点编号i为奇数时，此点与编号为(i-1)的点共处于同一条路段上，对每一个端点寻找离其最近的节点，并计算出两点构成矢量的方位角；路段属性ID从1开始编号，值为int(i/2)+1，同时，初始化ID查找表；

(3.4.b)对端点Ver[i]，当i为偶数时，通过对比Ver[i]与Ver[i+2]，Ver[i+3]，...，Ver[n]的距离，找出最短距离d_min对应的端点Ver[j]；当i为奇数时，则对比Ver[i+1]，Ver[i2]，...，Ver[n]的距离，找出最短距离点Ver[j]；

(3.4.c)当Ver[i]与Ver[j]间距离d_min小于阈值d_t时，说明已满足矢量间端点距离规则，则进行矢量间夹角规则判断，此时判断Ver[i]与Ver[j]所在路段的角度关系，根据公式θ＝|θ₁-θ₂|求出θ，若θ＞90°，则进入(3.4.d)，开始连接两条路段，若两端点间距离d_min大于阈值d_t，则不连接路段；

(3.4.d)根据ID＝int(i/2)+1定位当前处理路段a，根据ID＝int(j/2)+1定位当前处理路段b，若无法直接找到该路段，说明路段ID已在之前连接时被修改，需进入ID查找表找到新的路段编号；当查找到路段a、b后，连接路段，并将路段属性ID值改为二者中较小值，同时，在ID查找表内添加此次修改记录；

(3.4.e)此时令(i+1)为i，返回(3.4.b)，直至所有端点遍历完毕。

进一步的，步骤4中通过交叉路口连接算法将路口处连接的具体实现方式如下，

4.1，对于端点Ver[i]，寻找其所在线段最近节点B，沿B与Ver[i]的方向按阈值d_t为长度做延长线；

4.2，判断是否有道路与当前延长线有交点，若有交点，则记录此点E，并连接Verpi]与E点，否则，直接进入4.3；

4.3，令(i+1)为i，重复搜索道路端点，直至无点可以连接，则循环结束。

相对于现有技术，本发明的有益效果是：

本发明交互简单，只需要确定道路宽度以及给定种子点即可流程化提取道路并获得路网。提取时算法工作区域的设定，能够很大提高提取道路的效率。道路形态学处理、连接算法、交叉路口连接算法能够将道路阴影、交叉路口等中断处进行合理连接。且简化算法的使用保证了提取的道路矢量可以根据道路弯曲程度自动排布点的密度。综上所述，本发明完成了道路半自动提取以及生成路网流程化处理，策略有效且效率较高。

附图说明

图1为发明的技术方案示意图；

图2为区域增长算法说明图，其中(a)(b)为一次交互的完整图，(c)(d)(e)(f)为局部放大图。

图3为闭运算处理前后对比图；

图4为Zhang-Suen细化算法中八邻域像素排布图；

图5为细化处理效果对比图；

图6为矢量化跟踪算法八邻域搜索顺序；

图7为矢量化跟踪算法处理效果图；

图8为剔除过短道路结果对比图，其中下图是上图圆圈内的局部放大图；

图9为道格拉斯-普克算法矢量简化原理图；

图10为不同简化阈值效果对比图；

图11为端点删除算法角度计算图；

图12为端点偏移原因示例图以及端点处理前后对比图；

图13为单条道路连接算法矢量间夹角规则原理图；

图14为单条道路连接效果图；

图15为交叉路口连接原理图，其中(a)为“T”型与“Y”型连接原理解释，(b)(c)为根据交互方式所得到的两种不同“十”字型连接策略。

图16为交叉路口连接前后对比图；

图17为不同类型单条道路提取及连接结果示例图，其中(a)(c)(e)为用户交互结果，(b)(d)(f)为单条道路提取结果。

图18为不同类型路网连接结果示例图，其中(a)(b)为交叉路口连接前的情况，(c)(d)为连接后结果。

图19为不同道路路网成图局部放大图，其中(a)(b)(c)为单条道路局部成图放大图，(d)(e)(f)为住宅区周边道路成图放大图，(g)(h)(i)为遮挡区域道路成图放大图，(j)(k)(l)为交叉路口道路成图放大图；

图20为用分辨率为0.2m的遥感影像进行全图道路提取，影像大小为5000像素*5000像素的提取结果，其中(c)(d)为(b)中提取的细节图；

图21为用分辨率为0.5m的遥感影像进行全图道路提取，影像大小为1500像素*2000像素的提取结果，其中(c)和(d)为(b)中局部细节展示。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图实施示例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本发明提供的一种交互式半自动高分辨率遥感影像道路以及路网提取由用户交互提取道路宽度及道路种子点、栅格影像处理、矢量处理、路网生成四个步骤构成。

步骤1：通过用户交互确定道路宽度以及种子点，同时确定算法工作范围。

本发明要求用户交互画出道路宽度，通过在道路上点击两点确定道路宽度。要求用户尽可能放大影像使道路宽度准确。此处将获得的道路宽度记为L。

本发明要求用户给定道路种子点，并根据种子点以及道路宽度计算出算法工作区。对种子点的要求：第一，种子点必须保证准确，即必须标记在正确的道路上。第二，随种子点标记，会同时生成一个缓冲区域，缓冲区域确定了算法的工作范围，故必须保证缓冲区域完全覆盖所要提取的道路。

制定缓冲区掩膜主要有两个原因：一是限制区域增长范围，避免因过度增长导致算法效率降低；二是减少种子点数量，仅需缓冲区覆盖道路即可，使提取方法相较于传统的目视解译更为简单。

缓冲区建立的基本思想是对给定的空间对象集确定一邻域，邻域大小由邻域半径R决定。实体对象O_i生成的缓冲区B_i可定义为：B_i＝{x：d(x，O_i)≤R}，即对象O_i缓冲距为R的缓冲区是距O_i全部距离d小于等于R的点集，d通常指最小欧氏距离。对于多个对象的集合O＝{O_i：i＝1，2，3，...，n}，其半径为R的缓冲区是单个对象的缓冲区的并集，即为线缓冲区的生成算法，常用的方法是双线圆弧法。通常以线为中心轴线，沿其法线方向，分别向两侧平移一定距离得到的两条线，并与在端点处形成的光滑曲线(多为半圆弧)接合形成的封闭区域。

根据种子点以及缓冲区技术得到的道路提取工作区如图2(a)所示，记缓冲区宽度为Buffer_Width，令Buffer_Width＝L。此时以种子点的连线为中线，以Buffer_Width为半径生成掩膜。道路周围灰色区域即为掩膜，其内部灰色线条即种子点连线。

步骤2：栅格影像处理并获得道路中心线，具体包括如下子步骤。

(2.1)区域增长

本发明采用区域增长技术获取初始道路提取结果。区域增长通过以某一像素为增长点，比较八邻域灰度是否与当前像素满足合并阈值要求，若满足即合并，并将合并点作为新的种子点继续搜索，若不满足即停止增长。其具有规则简单、运算速度快、可交互等特点，很适合交互式提取道路要求。本发明中，考虑到城区建筑物在R、G、B三波段与道路灰度差异不大，故将区域增长合并阈值设置为可调节阈值。在建筑物比较少的郊区，建议设定为40，在建筑物较多的城市，可以相应的降低阈值，建议为20。为满足较长道路提取要求，本发明采用多个种子点同时增长。同时考虑到区域增长全图搜索效率低下，易出现过增长问题，本发明在步骤一生成的缓冲区内进行区域增长的策略。

图2(b)为区域增长后结果。图2(c)、图2(e)为局部放大影像，图2(d)、2(f)为相应增长后影像。从图中可以看出，区域增长可以得到大致道路范围，但是边缘存在毛刺、不平滑，且道路中央因为阴影、遮挡，或者由于未增长到而存在很多空洞。故区域增长结果需要进一步优化。对于内部空洞，我们利用形态学闭运算进行优化。对于道路断裂问题，我们在得到道路矢量后，通过道路矢量连接算法来进一步优化。

(2.2)形态学闭运算

为解决区域增长存在的问题，得到区域增长结果后，对结果进行形态学闭运算处理，从而对初始提取结果进行进一步优化，填补部分空洞以及细碎图斑，从而得到更加平滑的效果。

数学形态学是修整边缘，填充孔隙最常用的方法。形态学基础的算子包括腐蚀、膨胀、开运算和闭运算。形态学闭运算定义如下：

其中，运算符·表示闭运算，A表示二值影像，B表示结构元素。和分别表示膨胀和腐蚀算子。

由于本发明区域增长后，存在的主要问题是同质区域空洞过多，故而在此发明中，我们将结构元素尺寸设计较大，选取了大小为15*15像素的方形结构元素。图3对比了使用闭运算前后的结果。

(2.3)形态学细化算法

道路提取的最终目的是为了获得矢量化的道路网。为了避免道路边缘不平滑对矢量化的影响，本发明提出通过细化获得道路中心线进而构建路网的策略。

图像细化是将二值影像中的道路提取结果进行骨架化处理，通过对道路像素的层层剥离，从而只留下一组保持其连通性的单像素宽度道路。在本发明中，道路中心线提取直接影响到路网的生成。而好的细化结果往往满足以下几个条件：

1、细化中心线的连通性不变；

2、线状道路只用单像素宽度点表示；

3、道路两端点不能缺失；

4、速度快。

针对以上几个方面，本发明采用了Zhang-Suen并行算法进行细化处理。Zhang-Suen算法是一种快速并行迭代算法，由于其迭代少，速度快，在交叉口具有优越的细化能力，广泛应用于图像细化领域。对于形态学闭运算得到的二值影像(其中道路灰度值为255，背景灰度值为0)，通过逐点搜索道路点，即像素值为255的点，判断当前点八邻域内是否符合规则，符合规则即删除当前点。具体规则如下：

1、不是内部点；

2、不是孤立点；

3、不是直线端点；

4、如果当前点是边界点，去掉后，连通分量不会增加，则可以删除。

具体实现主要通过以下两个过程来完成：

(2.2.a)如若当前点p₁(像素值为255)满足以下四个要求，则删除，即将其灰度值改为0：

2≤N(p₁)≤6 (2)

T(p₁)＝1 (3)

p₂·p₄·p₆＝0 (4)

p₄·p₆·p₈＝0 (5)

其中，N(p₁)表示p₁八邻域的非零点数量，p_i(i＝2，3，...，9)表示为八邻域中第i个像素，如图4所示，T(p₁)表示八邻域内由0到255的跃迁次数。

(2.2.b)保留(2)(3)条件不变，添加以下条件，若满足，则删除：

p₂·p₄·p₈＝0 (6)

p₂·p₆·p₈＝0 (7)

重复上述两步，直至没有符合上述要求的像素为止。图5展示了细化效果。

步骤3：矢量化处理并将矢量化结果进行修饰，连接路段获得单条道路，具体包括如下子步骤。

(3.1)矢量化处理

经过上述栅格影像层面处理，得到单像素宽的道路骨架。此时需要对道路骨架进行矢量化处理。矢量图相较于栅格图更易保存与修改，也更易于数据更新与分析。矢量化第一步需要对道路骨架进行跟踪，得到道路线的链码串，跟踪的完整性与准确性会很大程度影响后续矢量化结果。本发明道路跟踪算法主要由以下四个步骤构成：

3.1.a、新建一个储存点集的双向队列，从上到下、从左到右，全局搜索灰度值为255的初始点，记为FirstPoint，令当前点CurrentPoint＝FirstPoint，将此点坐标储存进队列中，并将灰度值置为0，进入3.1.b。若全局未搜索到灰度值为255的像素，则全部道路已被跟踪，跳出循环并结束算法。

3.1.b、在八邻域内搜索下一个点。其中八邻域搜索顺序如图6，其中p₁为当前点，1-8为搜索顺序，其中1位于p₁的左上方，2位于p₁的正上方，3位于p₁的正左方，4位于p₁的右上方，5位于p₁的左下方，6位于p₁的正右方，7位于p₁的正下方，8位于p₁的右下方。通过规定先搜索上邻域(1-4)，再搜索下邻域(5-8)，使搜索具有方向性，达到减少搜索次数的目的。若搜索到灰度为255的点，则将此点记为CurrentPoint，将此点坐标储存进双向队列尾部，并将灰度值置为0，重复步骤二。直到当前像素八邻域未搜索到灰度为255的点，则进入3.1.c。

3.1.c、由于搜索的初始点并非是端点，此时令CurrentPoint＝FirstPoint，进入步骤四，开始跟踪FirstPoint另一个方向的道路。

3.1.d、依照图6的规则搜索当前点的八邻域，若搜索到灰度为255的点，则将此记为当前点，将此点加入双向队列头部，并将灰度值置为0。重复3.1.d，直到当前像素八邻域未搜索到灰度为255的点，则完成一条道路的跟踪，即获得单条道路点集，然后再进入3.1.a。

图7为矢量化处理前后效果图，分别为原始影像，细化后的二值图以及矢量化后的结果。矢量化后，点用×符号来表示。从结果上明显可以看出矢量化处理后的点很密集，不利于道路后处理。

(3.2)简化算法

在上述道路跟踪完毕后，会存在一些过短道路，如图7(b)以及图8(a)中，细化会在道路两端产生短枝丫，且在道路周边，由于受到同谱异物、道路交叉口等影响，也会存在短小分叉，这些过短道路对构成路网无效且易对后续连接算法造成干扰。因而，本发明对跟踪得到的多个队列，通过长度阈值筛选将部分队列删除，方便后续矢量处理。筛选阈值T_delete定义为：

T_delete＝(50.0*resolution)*L/3.0 (8)

其中resolution指影像分辨率，若影像无地理编码，则将resolution赋值为1即可。L是用户交互中提取的道路宽度。对于无地理编码或者分辨率为1米的影像，通过经验判定当道路长度大约为3个像素时，50为最佳阈值，而针对有地理编码的影像，则根据上述公式，按比例对阈值进行调整，从而为不同分辨率影像寻找到其对应阈值。图8展示了剔除过短道路前后的结果，其中下图为上图圆圈内的放大图。

单条道路细、长且平滑，若逐像素去描述，不仅占用空间，且大量冗余也不利于地理空间分析。道路这一特性决定其更适合用关键点去描述。本发明针对矢量化后的道路，利用道格拉斯-普克算法(Douglas–Peucker Algorithm)，对点集进行简化。道格拉斯-普克算法具有平移和旋转不变性，通过设定一个距离阈值，判断每点是否满足阈值要求来删除点，从而用少量点描述曲线，以达到数据压缩的目的。该算法主要由以下三个步骤组成：

3.2.a、连接两端点A和B并记为保留点，计算两端点间各点到两端点连线的欧式距离，并求出距离最远的点P。若两端点间无点，则所有点已简化完毕，退出循环。

3.2.b、比较P至直线的距离与阈值的大小，若大于阈值，则保留P。否则删去当前点。

3.3.c、以A、P为端点重复3.2.a-3.3.c；同时，以P、B为端点重复3.2.a-3.3.c。

图9中展示了道格拉斯-普克算法原理，由左上到右下，根据箭头方向，依次将端点进行了判断并删除，最终用五个点的多段线来描述初始七个点的多段线。

简化阈值对算法效果影响很大。如图10，影像空间分辨率为1米，比较了不同简化阈值对结果的影响。其中图10(a)是矢量化跟踪结果，(b)是简化阈值T＝1时的简化结果，圈出的部分为本应该省略的点集，(c)是简化阈值T＝3时的简化结果，可以看出简化结果很好，既删除了不必要的节点，又保留了关键点，很好的贴合了路面。(d)是简化阈值T＝5的情况，其中部分关键点缺失，圈内中心线已经偏出路面。综上所述，当分辨率为1米时，简化阈值T＝3为最佳阈值。

在此发明中，考虑到不同分辨率的影响，我们将阈值T定义为：

T＝3.0*resolution (9)

其中resolution是指影像分辨率，3.0是在影像无地理编码或者分辨率为1米时得到的最优阈值。

(3.3)端点删除算法

如图12所示，(a)为区域增长及形态学闭运算处理后的结果，(b)为细化结果，(c)为简化结果，(e)、(f)与(g)为(c)中圆圈内的局部放大图，(i)、(j)与(k)为其方框内的局部放大图。从图(e)、(f)、(g)中可以明显看出，道路在细化过程中，首尾容易偏移。且如图13(i)、(j)、(k)所示，本发明采用的矢量线跟踪算法导致了较短路线剔除法无法删除全部路口分支。因此，本发明针对简化后的矢量进行端点分支、偏移处理。本发明通过判断道路的首尾端点方向的方式来剔除偏移过大的点，从而达到修剪分支及删除偏移的目的，也为单条道路连接、路网构成奠定了基础。

具体思路为，对每条线段，如图11所示，取与其端点C相邻的两点A与B，求向量与的夹角θ，当θ＜10°时，保留C点，否则删除C点。

(3.4)单条道路连接算法

一次画线，默认提取的是一条连续的道路。为解决如图14(a)、(b)、(c)所示的道路断裂问题，需要对端点处理后的矢量进行连接操作。连接规则由矢量端点间距离规则以及矢量间夹角规则共同组成。若两矢量路段满足连接规则，则将两条矢量路段连接，并将属性ID修改为同一个值。同时，建立ID查找表，记录每条矢量线原始ID与修改后ID，方便追踪每条路段的修改，保证了多条路段搜索连接正确。

连接规则中，矢量间端点距离规则由两部分构成。首先，矢量间端点距离必须为全局最小距离；其次，规定长度阈值d_t，若两端点间距离小于d_t，则可以连接。如公式10，长度阈值d_t与筛选阈值T_delete定义相同，这样可以使筛选掉的短小路段空缺得以正确连接。

d_t＝(50.0*resolution)*L/3.0 (10)

同时，根据同一条道路既长又平滑的特点，矢量间夹角规则定义如下，如图13，若已找到满足距离规则的矢量路段端点A与B，则寻找与端点距离最近的节点A’与B’。求出A’A与B’B的方位角θ₁与θ₂(θ₁，θ₂∈(0°，180°])，计算θ＝|θ₁-θ₂|，若θ＞90°，即满足夹角规则。

单条道路连接步骤如下：

3.4.a、建立端点容器Ver。将提取到的道路端点坐标依次压入Ver中。假设端点个数为n，端点编号为i(0，1，2，3，...，n)。此时，当端点编号i为奇数时，此点与编号为(i-1)的点共处于同一条路段上。对每一个端点寻找离其最近的节点，并计算出两点构成矢量的方位角，如图13中的θ₁与θ₂。路段属性ID从1开始编号，值为int(i/2)+1。同时，初始化ID查找表。

3.4.b、对端点Ver[i]，当i为偶数时，通过对比Ver[i]与Ver[i+2]，Ver[i+3]，...，Ver[n]的距离，找出最短距离d_min对应的端点Ver[j]；当i为奇数时，则对比Ver[i+1]，Ver[i+2]，...，Ver[n]的距离，找出最短距离点Ver[j]。

3.4.c、当Ver[i]与Ver[j]间距离d_min小于阈值d_t时，说明已满足矢量间端点距离规则，则根据图13进行矢量间夹角判断。此时判断Ver[i]与Ver[j]所在路段的角度关系，根据公式θ＝|θ₁-θ₂|求出θ，若θ＞90°，则进入3.4.d，开始连接两条路段。若两端点间距离d_min大于阈值d_t，则不连接路段。

3.4.d、根据ID＝int(i/2)+1定位当前处理路段a，根据ID＝int(j/2)+1定位当前处理路段b，若无法直接找到该路段，说明路段ID已在之前连接时被修改，需进入ID查找表找到新的路段编号。当查找到路段a、b后，连接路段，并将路段属性ID值改为二者中较小值。同时，在ID查找表内添加此次修改记录。

3.4.e、此时令(i+1)为i，返回3.4.b，直至所有端点遍历完毕。

图14为单条道路连接效果展示，其中(a)、(b)、(c)分别为原始影像、闭运算结果、细化结果，可见道路中间因树木遮挡而缺失，且无法通过形态学闭运算修复。(d)为实施单条道路连接算法后效果，(e)为根据道路宽度生成缓冲区，结合(f)属性表显示，可见该条道路已完全合并，属性ID正确。

步骤4：交叉路口连接，并生成路网。

路网构成即在多次提取多条道路，且单条道路都已准确连接的基础上，对所有道路端点进行道路交叉口判断并连接。如图16(a)、(b)、(c)所示，道路交叉口主要分“T”型、“Y”型与“十”字型路口。

其中，“T”型与“Y”型交叉路口，皆为一条道路提取准确，另一条道路末端留有一端点未连接。如图15(a)，此时的连接策略为反向延长留有端点的道路，在一定范围内寻找是否有其他道路与之相交，若相交，则连接。

对于“十”字型道路，本发明根据提取方式的不同分为两类，第一种如图15(b)，即在提取单条道路时，首先将横向道路全部提取，之后，再将纵向道路全部提取，二者直接通过单条道路处理方式完成连接，这种情况已是连接正确，在交叉路口连接算法中不予考虑。但在道路不规则地区，如图15(c)所示，由于两条道路无法直接连接，“十”字型道路呈现为一条主干道提取准确，另两条道路末端留有未连接端点。此时，如图15(c)，分别反向延长两条留有端点的道路，在一定范围内判断是否有道路与其相交，若有，则连接。

根据上述分析，道路交叉路口连接步骤如下：

4.1、对于端点Ver[i]，寻找其所在线段最近节点B，沿B与Ver[i]的方向按阈值d_t为长度做延长线。

4.2、判断是否有道路与当前延长线有交点，若有交点，则记录此点E，并连接Ver[i]与E点。否则，直接进入3。

4.3、令(i+1)为i，重复搜索道路端点，直至无点可以连接，则循环结束。

图16展示了道路交叉口连接结果，以典型的“Y”型路口、“T”型路口，“十”字型路口为例，可以明显看到，通过单条道路连接以及单条道路反向延长连接两个策略后，路口有很好的连接效果。

最后，利用缓冲区技术，按步骤1提取的道路宽度以及步骤2提取的道路中心线，生成最终路网。

分别以不同分辨率的高分辨率航空影像为数据源进行不同道路提取实验，原始影像以及提取结果如图17至21所示。图17展示了单条道路提取以及连接后的结果对比；图18展示了不同类型道路路网连接前后对比；图19展示了不同环境下道路路网成图后的局部放大图；图20展示了0.2m分辨率影像全图道路提取结果以及细节；图21展示了0.5m分辨率影像全图道路提取结果以及细节。实验结果表明本发明提出的交互式提取高分辨率遥感影像道路并生成路网的方法能够快速准确地提取出高分辨率道路，能够合理地连接道路断裂处和交叉口并生成路网。

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

Claims (8)

1.一种交互式提取高分辨率遥感影像道路并生成路网的方法，其特征在于，包括以下步骤：

步骤1，通过用户交互确定当前等级道路宽度以及种子点，同时确定算法工作区域，即缓冲区；

步骤2，栅格影像处理，包括如下子步骤；

步骤2.1，在步骤1中的缓冲区内进行区域增长，获得道路初始提取结果；

步骤2.2，利用形态学闭运算填充空洞；

步骤2.3，利用形态学细化操作获取道路中心线；

步骤3，单条道路矢量化处理，包括如下子步骤；

步骤3.1，利用矢量化处理获得单条道路路段点集；

步骤3.2，首先对单条道路路段点集进行预处理，然后利用道路拉普拉斯-普克算法对点集进行简化；

步骤3.3，利用端点删除算法对路段进行矢量修饰；

步骤3.4，利用单条道路连接算法对单条道路进行连接；

步骤4，通过交叉路口连接算法将路口处连接，并利用缓冲区技术结合步骤1获取的道路宽度生成路网。

2.根据权利要求1所述的一种交互式提取高分辨率遥感影像道路并生成路网的方法，其特征在于：步骤2.2中采用了Zhang-Suen并行算法进行细化处理，具体实现方式如下，

(2.2.a)如若当前中心点p₁满足以下四个要求，p₁的像素值为255，则删除，即将其灰度值改为0：

2≤N(p₁)≤6 (2)

T(p₁)＝1 (3)

p₂·p₄·p₆＝0 (4)

p₄·p₆·p₈＝0 (5)

其中，N(p₁)表示p₁八邻域的非零点数量，p_i(i＝2，3，...，9)表示为八邻域中第i个像素，T(p₁)表示八邻域内由0到255的跃迁次数；

(2.2.b)保留(2)(3)条件不变，添加以下条件，若满足，则删除：

p₂·p₄·p₈＝0 (6)

p₂·p₆·p₈＝0 (7)

重复上述两步，直至没有符合上述要求的像素为止。

3.根据权利要求1所述的一种交互式提取高分辨率遥感影像道路并生成路网的方法，其特征在于：步骤3.1的具体实现方式如下，

(3.1.a)新建一个储存点集的双向队列，从上到下、从左到右，全局搜索灰度值为255的初始点，记为FirstPoint，令当前点CurrentPoint＝FirstPoint，将此点坐标储存进队列中，并将灰度值置为0，进入(3.1.b)，若全局未搜索到灰度值为255的像素，则全部道路已被跟踪，跳出循环并结束算法；

(3.1.b)在八邻域内搜索下一个点，其中八邻域搜索顺序为：以p₁为中心点，设1-8为搜索顺序，其中1位于p₁的左上方，2位于p₁的正上方，3位于p₁的正左方，4位于p₁的右上方，5位于p₁的左下方，6位于p₁的正右方，7位于p₁的正下方，8位于p₁的右下方，通过规定先搜索上邻域(1-4)，再搜索下邻域(5-8)，使搜索具有方向性，达到减少搜索次数的目的；若搜索到灰度为255的点，则将此点记为CurrentPoint，将此点坐标储存进双向队列尾部，并将灰度值置为0，重复(3.1.b)；直到当前像素八邻域未搜索到灰度为255的点，则进入(3.1.c)。

(3.1.c)由于搜索的初始点并非是端点，此时令CurrentPoint＝FirstPoint，进入步骤(3.1.d)，开始跟踪FirstPoint另一个方向的道路；

(3.1.d)依照八邻域搜索顺序搜索当前点的八邻域，若搜索到灰度为255的点，则将此记为当前点，将此点加入双向队列头部，并将灰度值置为0；重复(3.1.d)，直到当前像素八邻域未搜索到灰度为255的点，则完成一条道路的跟踪，即获得单条道路点集，然后再进入(3.1.a)。

4.根据权利要求3所述的一种交互式提取高分辨率遥感影像道路并生成路网的方法，其特征在于：步骤3.2的具体实现方式如下，

首先对单条道路路段点集进行预处理，具体为：对跟踪得到的多个队列，通过长度阈值筛选将部分队列删除，长度阈值T_delete定义为：

T-delete＝(50.0*resolution)*L/3.0 (8)

其中resolution指影像分辨率，若影像无地理编码，则将resolution赋值为1即可，L是用户交互中提取的道路宽度；

然后利用道路拉普拉斯-普克算法对点集进行简化，具体实现如下：

(3.2.a)连接两端点A和B并记为保留点，计算两端点间各点到两端点连线的欧式距离，并求出距离最远的点P，若两端点间无点，则所有点已简化完毕，退出循环；

(3.2.b)比较P至直线的距离与阈值的大小，若大于阈值，则保留P，否则删去当前点；

(3.2.c)以A、P为端点重复(3.2.a)-(3.2.c)；同时，以P、B为端点重复(3.2.a)-(3.2.c)。

5.根据权利要求4所述的一种交互式提取高分辨率遥感影像道路并生成路网的方法，其特征在于：步骤(3.2.b)中阈值T定义为：

T＝3.0*resolution (9)

其中resolution是指影像分辨率，3.0是在影像无地理编码或者分辨率为1米时得到的最优阈值。

6.根据权利要求4所述的一种交互式提取高分辨率遥感影像道路并生成路网的方法，其特征在于：步骤3.3的具体实现方式如下，

对每条线段，取与其端点C相邻的两点A与B，求向量与的夹角θ，当θ＜10°时，保留C点，否则删除C点。

7.根据权利要求1所述的一种交互式提取高分辨率遥感影像道路并生成路网的方法，其特征在于：步骤3.4中单条道路连接的连接规则由矢量端点间距离规则以及矢量间夹角规则共同组成，若两矢量路段满足连接规则，则将两条矢量路段连接，并将属性ID修改为同一个值；同时，建立ID查找表，记录每条矢量线原始ID与修改后ID，方便追踪每条路段的修改，保证了多条路段搜索连接正确；

其中连接规则中，矢量间端点距离规则由两部分构成；首先，矢量间端点距离必须为全局最小距离；其次，规定阈值d_t，若两端点间距离小于d_t，则可以连接，如公式10，阈值d_t的定义为，

d_t＝(50.0*resolution)*L/3.0 (10)

其中resolution指影像分辨率，若影像无地理编码，则将resolution赋值为1即可，L是用户交互中提取的道路宽度；

矢量间夹角规则定义如下，若已找到满足距离规则的矢量路段端点A与B，则寻找与端点距离最近的节点A’与B’，求出A’A与B’B的方位角θ₁与θ₂(θ₁，θ₂∈(0°，180°])，计算θ＝|θ₁-θ₂|，若θ＞90°，即满足夹角规则；

单条道路连接的具体步骤如下，

(3.4.a)建立端点容器Ver，将提取到的道路端点坐标依次压入Ver中；假设端点个数为n，端点编号为i(0，1，2，3，...，n)，此时，当端点编号i为奇数时，此点与编号为(i-1)的点共处于同一条路段上，对每一个端点寻找离其最近的节点，并计算出两点构成矢量的方位角；路段属性ID从1开始编号，值为int(i/2)+1，同时，初始化ID查找表；

(3.4.b)对端点Ver[i]，当i为偶数时，通过对比Ver[i]与Ver[i+2]，Ver[i+3]，...，Ver[n]的距离，找出最短距离d_min对应的端点Ver[j]；当i为奇数时，则对比Ver|i+1]，Ver[i+2]，...，Ver[n]的距离，找出最短距离点Ver[j]；

(3.4.c)当Ver[i]与Ver[j]间距离d_mmin小于阈值d_t时，说明已满足矢量间端点距离规则，则进行矢量间夹角规则判断，此时判断Ver[i]与Ver[j]所在路段的角度关系，根据公式θ＝|θ₁-θ₂|求出θ，若θ＞90°，则进入(3.4.d)，开始连接两条路段，若两端点间距离d_mmin大于阈值d_t，则不连接路段；

(3.4.d)根据ID＝int(i/2)+1定位当前处理路段a，根据ID＝int(j/2)+1定位当前处理路段b，若无法直接找到该路段，说明路段ID已在之前连接时被修改，需进入ID查找表找到新的路段编号；当查找到路段a、b后，连接路段，并将路段属性ID值改为二者中较小值，同时，在ID查找表内添加此次修改记录；

(3.4.e)此时令(i+1)为i，返回(3.4.b)，直至所有端点遍历完毕。

8.根据权利要求7所述的一种交互式提取高分辨率遥感影像道路并生成路网的方法，其特征在于：步骤4中通过交叉路口连接算法将路口处连接的具体实现方式如下，

4.1，对于端点Ver[i]，寻找其所在线段最近节点B，沿B与Ver[i]的方向按阈值d_t为长度做延长线；

4.2，判断是否有道路与当前延长线有交点，若有交点，则记录此点E，并连接Ver[i]与E点，否则，直接进入4.3；

4.3，令(i+1)为i，重复搜索道路端点，直至无点可以连接，则循环结束。