基于高分辨率影像的道路中心线提取技术研究

周绍光，向 晶，邱 伟，孙金彦，凡 莉

(河海大学地球科学与工程学院测绘科学与工程系，南京 210098)

0 引言

道路的检测与提取为地理信息系统(GIS)数据库的更新以及 GIS应用［1－2］提供了基础数据，具有十分重要的现实意义。针对遥感影像道路提取问题，国内外学者提出了很多算法和模型，朱长青等［3］提出了基于形态分割技术的道路网络提取方法;李卉［4］结合LiDAR数据和遥感影像信息提取城市道路并进行建模;Zhang等［5］使用拉东变换(radon transform)探测道路中心线和估计道路宽度;周绍光等［6］结合图割(graph cuts)理论和形状先验的思想，从高分辨率遥感影像中检测道路段。此外，曾发明等［7］、陈卓等［8］、蔡红玥等［9］等也给出了很好的道路提取策略。

部分研究着力于从获取的道路条带中提取道路中心线［10－11］，其中细化算法［12］是最普遍的，但该算法会产生很多“毛刺”，降低了道路中心线的光滑性和准确性。近期，一种叫做子空间约束的均值飘移算法［13］被应用到由道路条带获取精确中心线的研究中，该方法可以去除噪声点，获取的道路中心线比较光滑、没有“毛刺”;但该方法是有偏估计［14］，严重影响了道路中心线提取结果的精确性。

本文提出了一种基于道路条带自动检测道路中心线的方法。首先基于概率增进树算法获得道路候选点，并结合形态学运算获得道路条带;然后对获取的道路条带进行细化、获得道路骨架数组，通过迭代提取最长测地线，获得初始道路中心线;再根据Dijkstra最短路径算法优化每条道路中心线;最后断开道路中心线中角度变化大的部位，并根据方向一致性和道路连续性重新连接道路中心线段，获得最终结果。

1 道路中心线提取算法

图1示出道路中心线提取的技术流程。

图1 道路中心线提取技术流程Fig.1 Technical flowchart of road centerline extraction

1.1 道路条带的获取

利用概率增进树算法［15］从遥感影像中提取感兴趣目标的轮廓线是一种独具特色的方法，相关的研究成果已证明了这一策略的效果。本文借鉴该思想，探索了从遥感影像中提取道路的方法。该方法分为学习和测试2个阶段:①在学习阶段，选择并组合大量不同的特征，将概率增进树与决策树相结合，通过训练得到分类模型;②在测试阶段，利用训练好的模型计算出影像中各点属于道路的概率，即P(S/IN(c))，其中:IN(c)为实验影像I中以c为圆心、具有一定半径的图斑(本文中半径取25个像素的长度);S为判断道路点函数，若S=1则属于道路点，S=0则为非道路点。通过以上所建模型对测试影像进行实验，获得了道路候选点。

概率增进树算法几乎能检测出所有的道路点。但鉴于概率点不能作为分割结果，并且不可避免地会有些灰度、纹理与道路相似的房屋点被误判为道路点，而且路面上存在“孔洞”;故本文结合路径形态学［16－17］以及数学形态学中的一些其他运算对道路候选点进行处理，获得了令人满意的道路条带。

1.2 初始道路中心线的获取

细化算法是一种从道路条带中提取中心线的有效方法，但提取结果会伴随着许多“毛刺”，并且局部方向不一致。本文在保证道路网拓扑关系正确的前提下，采取求各对端点间测地距离的策略来消除“毛刺”。初始道路中心线获取的流程如图2所示。

图2 初始道路中心线获取流程Fig.2 Flow chart of initial road centerline acquisition

图2中提到的阈值是根据骨架图中最短道路中心线的长度确定的，如果阈值过大，大于该阈值的“毛刺”将被错分为道路中心线。经过该流程，起始道路中心线被存入一个新的数组，而骨架图中只剩下“毛刺”。

1.3 道路中心线的优化

上文去除了骨架图中的“毛刺”并获得了起始道路中心线，但局部道路中心线不够平滑。本文期望获得的道路中心线走向与参考方向一致，位置必须尽量接近道路条带中心。为此，引入能量最小［18］框架来优化道路中心线。

1.3.1 直线段获取

道路大多数呈直线状，或者稍微弯曲，对于后者也可将其视为由若干连续的直线段组成。这里要用到直线匹配法［19］获取一条初始道路中心线上所有像素的局部方向。直线匹配法的核心思想是将矩形道路视为由一系列的匹配直线构成，提取道路段就是检测出具有一定方向和长度的匹配直线。该方法最终获得每个像素的距离和方向矩阵，此处只用到了方向数据。直线段获取的具体步骤为:

1)用一个5像素×5像素的移动窗口，对所有位置的值为1的结构元素膨胀起始道路中心线(直线匹配法不适用于单像素宽曲线);

2)利用步骤1)得到的数组获取方向数据，取出起始道路中心线上每个点的局部方向数据;

3)将具有相同局部方向(差异在5°范围内)的连续点归为同一条直线段;

4)一条起始道路中心线被分为若干连续的直线段，设邻近直线段的交点为关键点(见图3中紫色的点)。值得指出的是，一条道路中心线的2个端点必须被设为关键点;

图3 沿道路中心线的关键点和道路截面Fig.3 Crucial points and road sections determ ined along road centerlines

5)对中心线上相邻关键点间的像素进行直线拟合，并将拟合直线的方向作为相应中心线段的参考方向。

1.3.2 道路截面确定

沿着关键点处曲线的法线方向取截面(如图3中的蓝色线段)，截面的长度需包含最优中心点。为便于理解，本文用Vi={v1，…，vj，…，vm}代表第i个截面的点集，一个截面中的所有点都分别与其邻近截面中的任意点相连，这些连接的短线段组成了所有的候选线。优化的过程实际上是从每个截面点集中获取最优点，并顺次连接成线。该思想充分利用了道路的形状特征，降低了初始中心线上局部方向波动的影响。

1.3.3 路径代价设计

考虑到中心线的走向应与参考方向一致、位置应尽量接近道路条带的中心这2条优化规则，路径代价将受角度约束和影像约束这2个因素制约，故将能量函数定义为

式中:Eangle为测定候选线方向与相应参考方向的差异;Eimage为影像约束能量;α为用于权衡2个约束能量相对重要性的参数。如果α值很大，相邻两截面间最优候选线的方向和相应的参考方向可能存在很大的差异;如果α值过小，则优化的道路中心线可能会偏离条带中心很多。

图4示出角度约束的含义。

图4 角度能量的含义Fig.4 M eaning of angle energy

图4中vj和vl分别为截面点集Vi和Vi+1中的任意一点，红色实线代表第i个截面和第i+1个截面之间候选线的参考方向 φi，i+1，φvjvl代表候选线段vj vl的方向。则第一个约束条件可定义为

角度差异越大，代价越大。

除了角度约束，影像约束也至关重要，它决定道路中心线能否尽量位于道路条带的中心。本文首先计算影像中任一像素最邻近非零点到该点的欧式距离，并将该距离值作为该像素的灰度值;则道路条带部分的点像素值大于0，并且越靠近条带中心、点像素的值越大;非道路条带部分点的像素值都为0，因为其到自身的距离最近。为计算方便，对距离转换结果取反，作为最终的参考影像(图5)。

图5 影像约束示意图Fig.5 Sketch map of image constraints

假设图4中候选线段vj vl由m个点组成，则局部影像代价定义为

式中I为参考影像。影像灰度值越大，影像代价消耗越大。

1.3.4 基于Dijkstra的最短路径搜索

Dijkstra 算法是典型的最短路径算法［20－21］，用于计算一个节点到其他所有节点的最小代价。定义C为带权有向图的点集，设定源点并初始化S=0，S为用于存放已求出最短路径的顶点集合，则C－S代表尚未求出最短路径的顶点集合。

本文中点集C由所有的截面点集{V1，…，Vi，…，VN}组成。如图6所示，假设在中心线的一端有一点设为源点，该点与V1点集中的所有点相连;此外，第二个截面中的一点vl分别与V1中的所有点相连。事实上，截面点集V2中的任意一点(如点vl)都与第一个截面中所有点相连，而其它任意2个截面之间点的连接方式也与此类似。

图6 连接图Fig.6 Connection image

候选线的传播代价按式(2)计算。非邻近的截面点不相连，权重设为∞。本文的目的在于搜索出源点到各个节点的最短路径，当到达截面Vi时，便得到了从源点到当前截面各点的最小代价，将该代价存于{D1，…，Di，…，DN}，其中 Di={di1，…，dij，…，dim}。

初始化中，将连接源点与V1中各点的候选线的权重设为0;添加V1至S。当S=C，则执行以下步骤:

1)对于Vi+1中的任一点vl，计算该点与Vi中各点相连所得候选线的传播代价，记为El={e1l，…，ejl，…，em l};

2)计算Di+El，如果所得向量的最小值为dij+ejl，则代表从源点到Vi+1中点vl的最短路径经过Vi中的vj点，更新di+1，l=dij+ejl;

3)当顺序遍历完Vi+1中的所有点，则添加Vi+1至S，并更新i=i+1。

由以上步骤获取了从源点至所有截面点的最短路径，则DN中的最小值即为优化该条中心线的最小代价。找出该代价所对应的各截面点并顺次连接，即为优化结果。

1.3.5 道路中心线断开与重新连接

1.3.1 节中获得的起始中心线确定了道路骨架的连续性和延展性，后续优化过程是顺序针对每条中心线的。以上操作并没有考虑道路的方向一致性，因此在交叉口处存在不合理连接(图7(a)中红色椭圆内)。本文的策略是，先找到各个交叉口的近似中心点，然后以一定的半径去除错误的交叉口连接(本文半径取5个像素，图7(b)给出了去除交叉口后的道路点方向图)，最后根据方向的一致性和道路的连续性将这些断开的道路段重新连接，得到最终的道路中心线(如图7(c)所示，不同的颜色表示不同方向的道路)。

图7 道路中心线断开与重新连接Fig.7 Disconnection and reconnection of road central lines

2 实验结果与分析

选取1景500像素×500像素的高分辨率航空影像作为实验图像(图8)，该图像覆盖的区域为英国伦敦某一城区。

图8 原始图像Fig.8 Original image

首先采用概率增进树算法获得初始道路点。由于得到的初始道路点是属于道路的概率，所以非道路点也会有值，但都比较小(图9(a)中绿色椭圆内);然后结合路径形态学和其他形态学运算进行处理，得到光滑完整的道路条带(图9(b))。

图9 道路条带获取Fig.9 Road stripe acquision

图10(a)给出对图9(b)细化运算得到的骨架图;图10(b)为去除“毛刺”后的结果，但道路局部方向不一致，线条不够平滑。根据方向信息构建截面，并拟合得到参考方向;基于Dijkstra算法解算最小代价路径，优化后的道路中心线与原始图像叠置(图10(c))的效果很好。

图10 道路中心线获取Fig.10 Road centralline acquision

利用下列公式进行道路条带和中心线提取精度的定量评价，即

式(4)—(6)中:Rm为人工解译的道路;Rc，Rl和Re分别为正确提取的道路、遗漏的道路和错误提取的道路。

根据上述评价指标，对本文方法的道路提取精度做出统计(表1)。

表1 本文方法道路提取结果统计Tab.1 Statistic results of road extracted with method proposed in this paper (%)

从表1可看出，提取出的道路条带和中心线的准确率均大于90%，提取结果较为理想。

3 结论

本文给出的从高分辨率遥感影像中提取道路中心线的方法，其创新点在于借鉴概率增进树算法，并将其应用于道路提取领域，获得了完整的道路段;基于迭代求取最长测地距离的思想，获取初始道路中心线，构建能量函数;利用Dijkstra算法优化道路中心线，最终得到了精确的道路中心线。试验结果证明了该策略的有效性和可行性。

本文尚未验证该方法在大幅影像中的性能;另外，将本文策略用于道路边缘提取领域，还没有考虑立交路、环形路等形状特别复杂的道路。这些将是今后进一步研究的内容。

［1］ 史文中，朱长青，王 昱.从遥感影像提取道路特征的方法综述与展望［J］.测绘学报，2001，30(3):257－262.ShiW Z，Zhu CQ，Wang Y.Road feature extraction from remotely sensed image:Review and prospects［J］.Acta Geodaetica et Cartographica Sinica，2001，30(3):257－262.

［2］ 赵晓峰.高分辨率遥感影像城区道路提取方法研究［D］.武汉:华中科技大学，2010:1－2.Zhao X F.A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of Master of Engineering［D］.Wuhan:Huazhong University of Science and Technology，2010:1－2.

［3］ 朱长青，王耀革，马秋禾，等.基于形态分割的高分辨率遥感影像道路提取［J］.测绘学报，2004，33(4):347－351.Zhu C Q，Wang Y G，Ma Q H，et al.Road extraction from high－resolution remotely sensed image based on morphological segmentation［J］.Acta Geodaetica et Cartographic a Sinica，2004，33(4):347－351.

［4］ 李 卉.集成LiDAR和遥感影像城市道路提取与三维建模［J］.测绘学报，2011，40(1):133－133.Li H.Road extraction and modeling with LiDAR and RS image in urban area［J］.Acta Geodaetica et Cartographic a Sinica，2011，40(1):133－133.

［5］ Zhang Q P，Couloigner I.Accurate centerline detection and line width estimation of thick lines using the radon transform［J］.IEEE Transactions on Image Processing，2007，16(2):310－316.

［6］ 周绍光，陈 超，岳建平.形状先验和图割的高分辨率遥感影像道路段提取［J］.测绘学报，2014，43(1):60－65.Zhou SG，Chen C，Yue JP.Extracting roads from high－resolution RS images based on shape priors and graph cuts［J］.Acta Geodaetica et Cartographic a Sinica，2014，43(1):60－65.

［7］ 曾发明，杨 波，吴德文，等.基于Canny边缘检测算子的矿区道路提取［J］.国土资源遥感，2013，25(4):72－78.doi:10.6046/gtzyyg.2013.04.12.Zeng FM，Yang B，Wu DW，et al.Extraction of roads in mining area based on Canny edge detection operator［J］.Remote Sensing for Land and Resources，2013，25(4):72－78.doi:10.6046/gtzyyg.2013.04.12.

［8］ 陈 卓，马洪超，李云帆.结合角度纹理信息和Snake方法从Li-DAR点云数据中提取道路交叉口［J］.国土资源遥感，2013，25(4):79－84.doi:10.6046/gtzyyg.2013.04.13.Chen Z，Ma H C，Li Y F.Extraction of road intersection from Li-DAR point cloud data based on ATSand Snake［J］.Remote Sensing for Land and Resources，2013，25(4):79－84.doi:10.6046/gtzyyg.2013.04.13.

［9］ 蔡红玥，姚国清.基于分水岭算法的高分遥感图像道路提取优化方法［J］.国土资源遥感，2013，25(3):25－29.doi:10.6046/gtzyyg.2013.03.05.CaiH Y，Yao GQ.Optimized method for road extraction from high resolution remote sensing image based on watershed algorithm［J］.Remote Sensing for Land and Resources，2013，25(3):25－29.doi:10.6046/gtzyyg.2013.03.05.

［10］ Zhang Q P，Couloigner I.Accurate centerline detection and line width estimation of thick lines using the Radon transform［J］.IEEE Transactions on Image Processing，2007，16(2):310－316.

［11］ Miao Z L，Shi W Z.Road centreline extraction from classified images by using the geodesic method［J］.Remote Sensing Letters，2014，5(4):367－376.

［12］ Soille P.Morphological Image Analysis:Principles and Applications［M］.New York:Springer Verlag，1998.

［13］ Ozertem U，Erdogmus D.Locally defined principal curves and surfaces［J］.Journal of Machine Learning Research，2011，12:1249－1286.

［14］ Genovese C R，Perone－Pacifico M，Verdinelli I，et al.Nonparametric ridge estimation［J］.Institute of Mathematical Statistics，2014，42(4):1511－1545.

［15］ Tu ZW.Probabilistic boosting－tree:Learning discriminative models for classification，recognition，and clustering［C］//Proceedings of the Tenth IEEE International Conference on Compute Vision.Beijing:IEEE，2005:1589－1596.

［16］ Randen T，Husoy JH.Filtering for texture classification:A comparative study［J］.IEEE Transactions on Pattern Analysis ＆ Machine Intelligence，1999，21(4):291－310.

［17］ Meli L，Santiago JM，Lodge TP.Path－dependent morphology and relaxation kinetics of highly amphiphilic diblock copolymer micelles in lonic liquids［J］.Macromolecules，2010，43(4):2018－2027.

［18］ Deo N，Pang C Y.Shortest－path algorithms:Taxonomy and annotation［J］.Networks，1984，14(2):275－323.

［19］ Sperling G.Binocular vision:A physical and a neural theory［J］.The American Journal of Psychology，1970，83(4):461－534.

［20］ ShiW Z，Zhu CQ.The line segment match method for extracting road network from high－resolution satellite images［J］.IEEE Transaction on Geoscience and Remote Sensing，2002，40(2):511－514.

［21］ Ahuja R K，Magnanti T L，Orlin JB.Network Flows:Theory，Algorithms and Applications［M］.Englewood Cliffs，NJ:Prentice－Hall，1993.