从遥感影像上自动提取道路对影像匹配、地形图更新等具有重要意义，国内外学者对此进行了大量研究。从提取方法的理论基础来看，现有的道路提取方法主要分为基于边缘特征的方法和基于区域特征的方法两大类别。线特征检测[1, 2, 3]是最常用的基于边缘特征的方法。动态规划方法是基于边缘特征的另一大类方法，其基本思想是在最优化代价函数的过程中不断地令曲线趋近道路轮廓，典型方法包括Candy模型方法[4]、主动轮廓方法[5, 6]等。除此之外，很多学者提出了改进措施，包括：边缘连接与重组[7, 8]、张量投票[9]等。基于区域特征的方法主要考虑到道路表面体现出较强的一致性。该类方法的前期环节为通过分割得到道路的候选图斑，典型方法包括：阈值分割[10]、形态学分割[11]、均值漂移分割[12]等。在图斑基础上，采用监督[13]、非监督[14]等方法实现道路提取。
目前，针对道路提取尚缺乏公认的普适性方法，道路的宽度、方向等几何形态一致性差导致的匹配困难是主要原因之一。虽未有文献明确针对曲折道路遥感提取问题，但部分方法对曲折道路的提取也具有一定的适用性。典型方法包括双边界追踪[15]、道路剖面匹配[16]和旋转模板匹配[17]等。然而，这些方法往往计算量较大或对道路形态的一致性要求过高。本文提出一种基于圆投影匹配的道路追踪方法，通过将分布在圆周上的道路特征投影为一维特征向量，从而实现旋转不变特性，较大程度克服了弯曲道路几何形态多变的问题。同时，通过对模板半径的选择、圆投影向量的离散近似以及干扰校正等3个关键问题进行分析与改进，较大程度地降低了对模板一致性的要求，改善了模板匹配的计算效率，提高了方法的适用性。
2 圆投影匹配及其改进 2.1 标准圆投影匹配
模板匹配是一种在图像中寻找给定目标的方法。该方法对具有统一形态对象的识别较为有效，但难以实现旋转、变形目标的有效匹配。针对此问题，有学者提出了与方向无关的圆投影匹配方法，并将其应用于手写汉字识别[18]、影像匹配[19]等领域，文献[20]将该算法应用于遥感道路提取。圆投影匹配采用半径为R的圆形模板，通过逐层计算处于模板半径为r(0<r≤R)的同心圆环像素灰度平均值构造圆投影向量，后续相似度计算均在圆投影向量空间中进行，如图 1所示。
|
| 图 1 圆投影模板示意图 Fig. 1 Schematic of circular projection template |
上述圆投影向量的构造实现二维分布的影像灰度值向一维的圆投影向量的变换。当图像围绕圆心旋转时，其圆投影保持不变，故圆投影匹配具有方向无关的特性。
2.2 圆投影匹配算法分析与改进旋转无关性是圆投影匹配的优势。然而，圆投影匹配也具有明显的不足，主要表现在：
(1) 对模板半径的选择依赖过高，半径过小将无法反应道路的特征，过大将引入非道路对象的干扰。
(2) 圆投影的计算需遍历圆周上的全部像素，计算量极大。同时，离散造成的舍入误差将导致很多像素重复参与计算。
(3) 成像条件的不同造成模板与待匹配对象间存在较大差异，某些目标无法正确匹配。
针对上述问题，本文对算法进行了相应的改进：
(1) 通过研究模板半径与圆投影的关系，定量给出最优模板半径的选择方法。
(2) 提出圆投影向量近似度参数，利用该参数实现圆投影模板离散近似程度和计算量的控制。
(3)提出采用模板干扰校正和自适应匹配机制结合的方法，克服道路因亮度差异、宽度变化、杂物遮挡等因素导致的匹配困难。
文献[19]通过遍历整幅影像，匹配给定的小幅单一目标，作者针对第三点不足提出了针对亮度校正的改进措施，虽然在一定程度上克服了模板与目标因亮度差异造成的匹配困难，但难以适用道路这类线状连续目标。
2.2.1 最优模板半径考虑路面像素灰度值为n，宽度为D，背景像素灰度值为m的理想情况(如图 2(a)所示)。以道路中心作为模板中心，由定义可得半径r对应的圆投影P(r)为
从式(1)可知，选择模板半径r的取值应大于
，否则，模板将丢失道路重要的边界信息从而导致任何面积大于模板的单色对象都将被错误匹配。然而，模板半径过大将增大计算量，且增加覆盖不相关的对象的可能性，从而带来更大的匹配风险。
|
| 图 2 半径r与圆投影P(r)的关系示意图 Fig. 2 Schematic of relationship between r and P(r) |
图 2(b)、(c)为式(1)的函数示意图像，即表示在理想条件下半径r与圆投影P(r)的关系。从图中可以看出，(R′,P(R′))为函数的拐点，当r>R′，圆投影的变化减缓，即半径的增大对后续圆投影的贡献减小。故通过求解拐点横坐标，得到最优半径
。采用该方法选择圆投影模板半径，可在获得最大道路特征的条件下使计算相对较小。
将模板半径R等分为nR层，对于其中的第k(1≤k≤nR)层，将其按照角度等分为C(k)=mk个间隔(m为一正整数常量)，即从第k层圆环上取C(k)个像素点。k与C(k)的正比例关系使得在较大的圆环上取较多的像素，这与圆投影的原始定义是一致的。设P(k)为第k层所对应的圆投影，其离散计算公式为
式中，Δr=R/nR是模板半径的采样间隔；Δθ=2π/C(k)是第k层圆投影的角度间隔；(x,y)为模板中心在影像上的像素坐标。以π/2作为最大角度采样间隔、2像素作为最小半径，即该半径对应的圆投影包括4个像素，以此为基准，按照正比例关系可得：当给定R和nR，则m=2R/nR，故nR可作为圆投影向量近似度参数。离散化在基本保持圆投影计算精度的前提下降低了计算量。 2.2.3 模板干扰校正
本文考虑线性亮度差异和加性随机噪声两类干扰情况，受干扰圆投影向量的计算模型为
式中，P是未受干扰的圆投影向量；a为亮度放大因子；B为灰度偏移向量；aP+B为亮度差异影响；N为随机噪声向量。针对随机噪声干扰N，本文仅对圆投影离散点进行3×3的均值滤波，而非遍历整幅影像，从而降低了计算量。针对亮度条件干扰，采用如下方法重建圆投影，得到其修正结果P*(k)
式(4)中，P′(k)为一灰度平移不变量，消除了灰度偏移B的影响；式(5)将分布在正负区间的距平值平移至正区间，δ>0为防止对数计算错误添加的微小常量；式(6)的对数运算将放大因子a的乘性影响转换为加性影响，再进行距平运算从而将其影响消除，从而消除了亮度放大因子a的影响。重建后的圆投影P*(k)为一亮度差异不变量。
图 3对比了采用干扰校正前后的P(k)和P*(k)曲线。可以看出，亮度差异和随机噪声造成的干扰对P*(k)几乎无影响，本文采用的抗干扰机制对上述两种类型的干扰抑制有效。
3 基于追踪的道路提取考虑到道路往往占据影像的面积比例较小，整幅影像的遍历计算量极大，而且全局处理将形成很多道路片段，后期还涉及大量连接和重组工作，故本文采用双向追踪的道路提取思路。
|
| 图 3 不同干扰条件下的P(k)和P*(k)曲线(nR=35) Fig. 3 P(k) and P*(k) curves under different interferences (nR=35) |
道路追踪从一个给定的种子点开始，依次向着两个方向进行。其基本步骤为：从道路上选定一种子模板，将当前追踪位置设置在该模板中心位置。围绕当前模板，以一定的步长和角度间隔，计算每个位置的圆投影向量及与种子模板间的相似度，选择相似度最大的圆投影中心作为当前追踪位置，重复上一过程，直至追踪完毕，如图 4(a)所示。当一个方向道路追踪完毕后，交换第1、第2个已匹配模板，按照上述方式进行反向追踪。为避免追踪的回溯，除第一步追踪外，将追踪方向限制在一个对称前进方向、张角为β的扇形区域内，如图 4(b)所示。该方法在自适应道路转弯的同时，有效防止追踪的回溯。图 4(c)为急转弯道路的追踪示例。
|
| 图 4 基于圆投影匹配的道路追踪 Fig. 4 Road tracking based on circular projection matching |
为增加方法的实用性，本文采用栈作为基本数据结构进行分叉口处理，以提取全部道路分支。在进行每一步追踪时，将可匹配的模板及其方向进行压栈处理，追踪按照相似度最大的方向继续进行。在当前道路追踪完毕后，将栈顶的模板及其方向进行出栈操作，并以该种子为起点，沿对应方向进行道路追踪。持续上一步操作直到栈空为止。通过这种方法，能够在一次追踪中提取道路主干及与模板所在道路相连的全部道路分支。
以图 5所示情况为例，在追踪到A—B分叉口时，将可匹配的模板及其方向PB→压栈，追踪继续沿着A方向进行(假设A方向相似度最大，下同)；在A—C、A—D分叉口时，分别将PC→、PD→压栈并继续沿着A方向追踪。当道路A追踪完毕后，从栈中逐次弹出PD→、PC→、PB→，分别完成道路D、C、B的追踪。
|
| 图 5 分叉口处理示意 Fig. 5 Schematic of branch path processing |
传统模板匹配往往仅计算待匹配点与初始选择模板的相似度。当nR=2，任意圆投影向量都可以表示为P(1)-P(2)平面上的点，如图 6所示。以Euclidean距离作为度量，给定模板的圆投影向量PT1，给定距离阈值RT1>0作为匹配半径，则传统模板匹配中，可匹配区域为一圆形区域，如图 6(a)阴影区域所示。定义待匹配圆投影向量与模板的相似度S(S∈[0, 1])为
式中，dT1为待匹配圆投影向量与模板圆投影向量的距离；dT1越小表示其与模板相似度越高。S=0表示待匹配圆投影与模板完全不能匹配。图 6(a)中P′和P″分别为可匹配和不可匹配圆投影向量的示例。
|
| 图 6 二维圆投影匹配示意图 Fig. 6 Schematic of 2D circular projection matching |
本文提出自适应模板匹配的思路为：给后续匹配成功的每个圆投影向量也赋予可匹配半径，即构建独立的可匹配区域(图 6(b)中小圆)，最终的可匹配区域被定义为这些区域与初始可匹配区域(图 6(b)中大圆)的并集，如图 6(b)阴影区域所示。为防止模板变形过大，要求前者的可匹配区域小于模板的可匹配区域，图 6(a)中的P″得以成功匹配，则为其构建模板PT2，其可匹配半径可计算为RT2=RT1/(1+dT1)。可以看出，RT2的大小与PT2—PT1间的距离反相关，这就意味着PT2距PT1越远，其在最终的模板中占的权重越低，以此防止模板的误差积累。在图 6(b)中，待匹配圆投影向量P'''落入可接受区域，故可得成功匹配，而P′则无法匹配。P'''的相似度定义为S′=(d2S1+d1S2)/(d1+d2)，其中d1、d2分别为P'''到PT1、P2T的距离。S1、S2分别为按式(7)计算得到的相似度。
一般的，设PT1、PT2、…、PT 为N个成功匹配的圆投影向量，P为待匹配圆投影向量，RT>0为给定的可匹配半径，dk=||P－Pk||2,(k=1,2,…,N)为P与PT1、PT2、…、PNT的距离，则待匹配圆投影向量的相似度
4 试验与分析 4.1 试验
试验A选取大小为分辨率2.5 m的山区 SPOT-5真彩色融合影像，由两幅影像拼接得到，尺寸为1153像素×862像素，如图 7所示。任意选择圆投影模板(子窗口0)，从图 7(a)可以看出，模板与目标道路在亮度和几何形态存在较大差异：影像拼接线上存在较大的亮度差异和几何错位(子窗口1-5)；在道路转弯和分叉口处存在较大的宽度和形态变化(子窗口6、7)。影像中有两条相交道路，总长度约为3.2 km。设置提取参数为：模板半径R=4像素、近似度参数nR=5、追踪张角β=120°、追踪步长Δd=8像素。首次追踪(即选择一次模板，下同)结果如图 8(b)所示。
|
| 图 7 具有亮度和形变干扰SPOT-5影像及道路提取结果 Fig. 7 Road tracking from SPOT-5 image with brightness differences and geometric interferences |
|
| 图 8 具有局部干扰的山区GeoEye-1影像及道路提取结果 Fig. 8 Road tracking from GeoEye-1 mountain road image with local interferences |
高分辨率遥感影像能够体现更多的地物细节，干扰表现得更加明显，从而导致模板匹配的难度大大增加。故试验B选取分辨率0.41 m的山区GeoEye-1全色影像，大小为819像素×744像素，如图 8(a)所示。图 8(a)的局部窗口显示了选择模板和五种典型干扰：路面材质差异(子窗口1)、分叉及杂物覆盖(子窗口2)、建筑物阴影(子窗口3)、树枝遮挡(子窗口4)和山体阴影(子窗口5)。设置的提取参数为：模板半径R=18像素、近似度参数nR=7、追踪张角β=120°、追踪步长Δd=40像素。首次追踪结果如图 8(b)所示。
4.2 试验结果分析从图 7(b)所示的提取结果可以看出，对于以SPOT-5为代表的、分辨率相对较低的遥感影像，本文的干扰矫正机制可较为有效地应对亮度差异和几何形态差异；在影像的边界处和转弯处，道路追踪也可以流畅完成。从图 7(b)的局部窗口可以看出，按照上文提出的分叉口处理方法，道路A追踪完毕后，道路B得以追踪；按照反向追踪机制，道路C也可得到追踪。然而，由于长度小于所设置的追踪步长，图 7(b)局部窗口黄圈所示的小段环状道路未得到有效提取。针对这一问题可通在小段道路上再次选择模板，设置更小的步长实现其提取。从图 8(b)所示的提取结果可以看出，本文提出的自适应模板机制可较为有效地适应遮挡、阴影等干扰。然而，图 8所示影像由于亮度和道路宽度差异过大的原因，首次追踪有两个道路分叉没有得到有效提取，如图 8(b)子窗口1、2所示。需再次进行模板选择以完成其余道路追踪。而且，提取得到的道路中心线较真实道路中心线相比存在偏移和震荡现象，其原因在于：道路追踪通过围绕当前模板，在一定的搜索步长和角度间隔下进行最相似圆投影模板来实现。为方便计算，上例中本文采用的角度间隔同所在圆周的离散角度间隔。然而不同的角度间隔决定了不同的方向分辨能力，往往会导致所选择的最相似圆投影模板中心偏离道路中心，以至某些追踪结果产生系统性偏移。或在道路转弯处产生一定的震荡现象。这一问题可通过人工指定更密集的角度间隔和更小的追踪步长加以克服。
本文采用中心线提取的位置偏差作为精度衡量指标。在上述试验中，SPOT-5影像道路中心线提取的最大位置偏差为3像素，即实际位置误差约7.5 m左右，平均位置偏差为1.5像素，即实际平均位置误差约3 m左右。GeoEye-1影像道路中心线提取的最大位置偏差为5像素，即实际位置误差约2 m左右，平均偏差为3像素，即实际平均位置误差约1.2 m左右。
为了比较算法的准确率和效率，本文选择平行线特征方法作为基于边缘方法的代表(文献[1])和均值漂移方法(文献[12])作为基于区域方法的代表，与本文提出方法首次追踪结果进行试验对比，如图 9所示。对3种方法的提取结果定量分析，其结果如表 1所示。
|
| 图 9 平行线特征方法和均值漂移方法应用于试验A、B数据的道路提取结果 Fig. 9 Extract results of parallel line feature and mean-shift method for images used in experiment A and B |
| 方法 | 准确率/(%) | 效率/s | |||
| 试验A | 试验B | 试验A | 试验B | ||
| 平行线特征 | 52 | 68 | 4.5 | 5.2 | |
| 均值漂移 | 63 | 72 | 37.4 | 32.1 | |
| 本文方法 | 95 | 77 | 0.47 | 0.36 | |
| 注：准确率为提取道路总长度(像素)与测试影像中道路中心线总长度(像素)的比值 | |||||
由于曲折道路直线特征不够明显，山区道路往往宽度变化多变，平行线特征难以得到保证，故平行线特征方法针对道路转弯(图 9(a)子窗口2、图 9(c) 子窗口3所示)或道路单侧遮挡(图 9(c)子窗口2、6所示)等情况的提取结果往往不够理想，且当影像中含有非道路平行线边缘时(图 9(c) 子窗口4所示房屋边缘)，容易产生误提取现象。故平行线特征方法更适用于城市道路的提取。由于直线边缘检测效率较高，且第2步的边缘重组是在矢量层面上进行，故该方法具有较高的效率。均值漂移分割的基本思想为邻域聚类，对全局的亮度差异不敏感，且具有一定的滤波效果，故道路的微小干扰可得到较好克服，对较为细小的道路分支也具有一定的提取效果(图 9(b)子窗口4、6和图 9(d)子窗口4、6所示)。但是，均值漂移分割得到的道路板块往往过于琐碎，容易造成道路中断现象(图 9(b) 子窗口1所示)，不适合线状对象的提取。而且该方法的时间复杂度为O(n2)级，且后期需要对分割结果进行过滤和重组等后处理。本文提出的方法因直接针对道路对象进行追踪，避免了全局搜索，效率和准确率均有提高。本文所提方法是一种半自动化方法，需人工选择模板，暂时还不具备全自动提取完整道路网络的能力。然而，该方法允许用户再次进行模板选择，从而实现整幅影像的道路提取。
5 结论和展望本文提出了一种基于改进圆投影模板匹配和追踪算法的道路提取方法。该方法是一种具有旋转不变特性、计算量较低的模板匹配算法。在圆投影向量的离散近似、最优模板半径选择和干扰校正三方面进行了改进。采用追踪算法进行道路提取，以扇形范围限制模板匹配范围，以适应道路转弯并防止回溯追踪。在模板匹配时，改善了相似性度量的传统方法，允许模板具有微小变形，在某种程度上克服待匹配道路与选择模板之间的变化。SPOT-5和GeoEye-1遥感数据试验结果表明，该方法在曲折道路信息提取中较一些常规方法效率和准确率均有提高。
抗干扰机制中仅考虑了亮度和随机噪声两种干扰，而未考虑树木、房屋严重遮挡带来道路中断等问题。此外，本文方法需要人工完成模板的选择，每次仅能提取模板所在道路及其分岔路，尚未实现全自动提取整幅图像完整道路网络。因此，增加模板抗干扰能力，且在首次道路追踪完成之后，自动选择新的模板实现全自动操作是下一步需要解决的问题。
| [1] | ZHONG Jiaqiang,WANG Runsheng. A Road Network Change Detection Algorithm Based on Linear Feature[J]. Journal of Remote Sensing, 2007, 11(1): 27-32. (钟家强, 王润生. 一种基于线特征的道路网变化检测算法[J]. 遥感学报, 2007, 11(1): 27-32.) |
| [2] | XIAO Zhiqiang,BAO Guangshu,HUANG Jixian. Update of Road Network in GIS by Fusing SAR and TM Imagery[J]. Acta Geodaetica et Cartographica Sinica, 2006, 35(1): 46-51. (肖志强, 鲍光淑, 黄继先. 融合SAR和TM图像更新GIS道路网络数据[J]. 测绘学报, 2006, 35(1): 46-51.) |
| [3] | LIU Xiaodan,LIU Yan. Urban Road Extraction Based on Hough Transform and Path Morphology[J]. Computer Engineering, 2012, 38(6):265-268. (刘小丹, 刘岩. 基于Hough变换和路径形态学的城区道路提取[J]. 计算机工程,2012, 38(6): 265-268.) |
| [4] | LACOSTE C, DESCOMBES X, ZERUBIA J. Point Processes for Unsupervised Line Network Extraction in Remote Sensing[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2005, 27(10): 1568-1579. |
| [5] | SONG M J. CIVCO D. Road Extraction Using SVM and Image Segmentation[J]. Photogrammetric Engineering and Remote Sensing, 2004, 70(12):1365-1371. |
| [6] | WU Liang, HU Yunan. Balloon Snake with Adjustable Expansion Coefficient in Road Contour Extraction[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(1): 71-77. (吴亮, 胡云安. 膨胀系数可调的Balloon Snake方法在道路轮廓提取中的应用[J]. 测绘学报, 2011, 40(1): 71-77.) |
| [7] | TREASH K, AMARATUNGA K. Automatic Road Detection in Grayscale Aerial Images[J]. Journal of Computing in Civil Engineering, 2000, 14(1): 60-69. |
| [8] | SENGUPTA S K, LOPEZ A S, BRASE J M, et al. Phase-based Road Detection in Multi-source Images[C]//Proceedings of IGARSS ′04: 2004 IEEE International Geoscience and Remote Sensing Symposium: 6. New York: IEEE, 2004: 3833-3836. |
| [9] | SHEN D J, HU C, YANG B, et al. Road Detection in SAR Images Using a Tensor Voting Algorithm[C]//Proceedings of the SPIE: 5th International Symposium on Multispectral Image Processing and Pattern Recognition: 678703. Wuhan: SPIE, 2007: 1-6. |
| [10] | YAN D M, ZHAO Z M. Road Detection from QuickBird Fused Image Using IHS Transform and Morphology[C]//Proceedings of Geoscience and Remote Sensing Symposium, 2003. IGARSS ′03: IEEE International Geoscience and Remote Sensing Symposium: 6. New York: IEEE, 2003: 3967-3969. |
| [11] | ZHU Changqing, WANG Yaoge, MA Qiuhe, et al. Road Extraction from High-resolution Remotely Sensed Image Based on Morphological Segmentation[J]. Acta Geodaetica et Cartographica, 2004, 33(4): 347-351. (朱长青, 王耀革, 马秋禾, 等. 基于形态分割的高分辨率遥感影像道路提取[J]. 测绘学报, 2004, 33(4): 347-351.) |
| [12] | LONG H. ZHAO Z M. Urban Road Extraction from High-resolution Optical Satellite Images[J]. International Journal of Remote Sensing, 2005, 26(22): 4907-4921. |
| [13] | HUANG X, ZHANG L. Road Centreline Extraction from High-resolution Imagery Based on Multiscale Structural Features and Support Vector Machines[J]. International Journal of Remote Sensing, 2009, 30(8): 1977-1987. |
| [14] | YUAN J Y, WANG D L, WU B, et al. Automatic Road Extraction from Satellite Imagery Using LEGION Networks[C]//Proceedings of IJCNN 2009: International Joint Conference on Neural Networks. Atlanta: IEEE, 2009: 3471-3476. |
| [15] | MA H J, QIN Q M. Road Extraction from ETM Panchromatic Image Based on Dual-edge Following[C]//IEEE International Geoscience and Remote Sensing Symposium, 2007: 460-463. |
| [16] | BAUMGARTNER A, HINZ S, WIEDEMANN C. Efficient Methods and Interfaces for Road Tracking[J]. International Archives of Photogrammetry Remote Sensing and Spatial Information Sciences, 2002, 34(3/B):28-31. |
| [17] | KIM T J, PARK S R, KIM M G, et al. Tracking Road Centerlines from High Resolution Remote Sensing Images by Least Squares Correlation Matching[J]. Photogrammetric Engineering and Remote Sensing, 2004, 70(12): 1417-1422. |
| [18] | CHIU H P, TSEN D C, CHENG J C. Invariant Handwritten Chinese Character Recognition Using Weighted Ring-data Matrix[C]//Proceedings of 3rd International Conference on Document Analysis and Recognition:1. Montreal: IEEE, 1995: 116-119. |
| [19] | XU Yibin, WANG Jingdong, LI Peng. Research on Scene Matching Method Using Circular Projection[J]. Systems Engineering and Electronics, 2005. 10(27): 1725-1728. (徐亦斌, 王敬东, 李鹏. 基于圆投影向量的景象匹配方法研究[J]. 系统工程与电子技术, 2005, 10(27): 1725-1728.) |
| [20] | WANG Ximing. The Research and Application of Road Extraction from Remote Sensing Images[D]. Beijing: Tsinghua University, 2011. (汪夕明. 遥感影像道路提取方法研究与实现[D]. 北京: 清华大学, 2011.) |