随着遥感技术的发展，遥感图像的空间分辨率越来越高，为人们的生产生活带来了更多的应用，亚米级以上分辨率的遥感数据已经成为目前主流的消费指标^[1]，而遥感图像中的道路提取是其重要的一类应用。但高分辨率遥感图像同时也带来了更多的干扰信息，对机器自动识别地面目标提出了更高的挑战^[2]。对于道路识别而言，在中低分辨率的遥感图像中，道路目标呈现为线状或细条状的特征，之前以该特征为道路模型提出了一系列算法，但这些算法在高分辨率影像中几乎都不再适用，因此研究高分辨率遥感图像的道路提取算法成为目前的研究热点^[3]。

根据人工干预的程度，可以将遥感道路提取分为人工提取、半自动提取和自动提取^[4]。自动道路提取，可以分为识别、测量、校验3个步骤，由于遥感图像中存在大量异物同谱的现象^[5]，让计算机自动识别和校验道路目标存在很大的困难。目前比较现实的算法是半自动算法，该类算法充分利用了人工的识别、检验能力和计算机的精确定位能力，以确保在较快的速度下提取可靠的道路信息^[6]。目前，研究人员已经提出多种较为有效的半自动道路提取算法。常见的有以下几种：活动轮廓模型法^[7-12]、纹理与形状特征法^[13-14]、最短路径法^{[2, 15]}、模板匹配法^[16-20]、概率统计法^[21-22]等。文献[3]提出了一种平行角纹理跟踪算法，其本质是一种矩形模板算法。该算法充分利用道路基元的纹理特征来克服阴影和车辆遮挡对道路提取的不利影响，具有较高的稳健性，但该算法无法应对道路宽度和角度急剧变化的情况。文献[16]提出一种最小二乘相关匹配的矩形模版算法，对弯曲道路具有一定的适应性，但对人工初始种子比较敏感。文献[19]提出了一种曲折道路自动追踪的算法，将二维的像素灰度匹配转换为一维的区域灰度特征匹配，提高了模板匹配效率，同时较好地解决了道路方向突变的问题，但该算法过度依赖图像灰度信息，无法应对异物噪声干扰问题。文献[20]提出了一种圆形模板匹配算法。该算法结合图像灰度和形态学梯度信息，采用迭代内插的方法搜寻最优道路中心点，算法效率较高。但该算法也存在诸多问题，比如模板尺寸需要人工设定、搜索中心点过程中容易受到道路两边平滑区域的干扰，中心点的概率估计需要人为设定过多参数，且参数的比重难以估计等。本文提出了一种改进的圆形模板匹配算法。该算法基于改进的局部形态学梯度图自动计算模板尺寸，利用道路显著图和几何形状信息搜索道路中心点，克服了道路两侧的平滑区域对道路中心点搜索的干扰，且所需人工设定的参数更少，参数更具直观的物理意义，提高了算法的适应性和实用性。

1 道路提取算法 1.1 高分辨遥感道路特征

在高分辨率遥感图像中，道路具有明显特征：①道路的几何形状表现为一条带状目标；②道路表面的灰度具有较高的一致性，且与道路两旁的目标具有较强的反差，但不同传感器拍摄的图像可能呈现截然不同的光谱特性，此特性对某些道路识别算法(特别是基于机器学习的算法)会带来困扰^[23]，本文在试验部分对该问题给出了阐述；③道路的几何形状多变，城市道路段一般是直线或者以一定角度缓慢变化的平滑曲线，但乡村道路可能呈现较大的弯度突变；④城市道路更容易受到噪声干扰，即城市道路常因建筑物阴影、车辆和行道树覆盖等噪声而使道路灰度突变，而乡村道路干扰噪声相对较少。

1.2 改进的局部形态学梯度图

为了自动计算圆形模板的尺寸，本文引入了一种改进的局部形态学梯度图。形态学梯度算子是一种非线性差分边缘检测算子，利用结构算子E(s, t)对灰度图像G(x, y)进行形态学膨胀、腐蚀、开运算、闭运算等操作(式(1)—(4))。传统的形态学梯度=膨胀-腐蚀，如文献[20]的做法。本文采用改进的形态学梯度算法^[24]，如式(5)所示，具有更好的性能，使得平滑区域响应更低，同时增强了边缘信息，如图 1所示。图 1(a)为原始图像，比较图 1(b)与图 1(c)的效果，可以发现改进的梯度图假边缘信息更少，同时提升了边缘的强度，使本文后续的模板尺寸自动估算更精确。

形态学梯度图的计算复杂度较大，一张N×N的图像，使用M×M的结构算子，式(5)的计算量为7×M²×N²，而工业生产中，道路矢量化的图像尺寸一般是5000×5000像素(按0.5 m分辨率的航拍图像6.25 km²的图幅计)，计算全图的形态学梯度图效率极低。在本文算法中，形态学梯度图仅用于计算每次提取的初始模板尺寸。为此本文提出局部形态学梯度图的概念，即以用户初始种子点为中心，截取一小片方形子图，按式(5)计算子图的形态学梯度图。为方便自动计算圆形模板半径，本文对形态学梯度图进行了归一化处理，如式(6)所示，后文提到的形态学梯度值均为归一化的值。

1.3 自适应模板生成

圆形模板匹配主要依靠模板内的局部灰度信息来识别道路位置，模板尺寸太小或太大均不利于道路中心位置的精确定位，合理的模板直径应约等于道路的宽度。本文利用改进的局部形态学梯度图，自动估算模板半径，同时自动将人工种子点调整到道路中心位置。首先人工粗略地在道路区域选取一个种子点，算法生成局部形态学梯度图，然后以种子点附近一定范围内(范围大小可以根据图像分辨率自动设定)的每个像素为圆心，生成半径为一个像素的等大圆形模板，记录所有模板投影内的形态学梯度值之和(简称梯度和)的最小值，然后半径增加一个像素，重复以上步骤。经过对大量图片的试验发现，当梯度和最小值刚超过半径值的模板为最佳模板。其原理为, 路面灰度具有较高的一致性，在梯度图中呈现为黑色区域，而道路两旁的其他目标与路面灰度具有较大的反差，在梯度图中呈现为强烈的边缘，当不断扩大的模板碰到边缘信息后，梯度和的值将会出现急剧的增大，很容易超过半径值。图 2显示了某次模板自动生成过程，其中图 2(a)为人工选择的道路种子点，图 2(b)为自动生成的圆形模板，图 2(c)显示了半径与梯度和的关联关系。可以发现当半径为8时，梯度和发生了剧烈的变化，因此取8为最佳半径。

自适应模板生成算法的步骤如下：

输入：

归一化局部形态学梯度图，初始道路种子点。

处理：

步骤1——初始化模板半径为1。

步骤2——以种子点附近邻域的每个像素为圆心，计算每个模板内的形态学梯度的累计值。

步骤3——获取所有梯度累计值的最小值，如果梯度值累计值小于半径，则半径长度增加1个像素，转到步骤2，否则转到步骤4。

步骤4——记录当前的模板半径和模板中心坐标。

输出：

模板半径和模板中心坐标。

1.4 改进的显著图

为了搜索道路中心点，本文定义如式(7)所示的显著图。该显著图改进自文献[2]，其中，avg(template)为最佳模板内的图像灰度平均值，因此显著图信息不再依赖于道路上的某个点的灰度，使得显著图信息更加稳定可靠。原图中与道路目标灰度相似的区域在显著图中具有更低的灰度值，如图 3(b)所示。显著图中灰度越低的区域，属于道路目标的概率越大，为算法后续的迭代内插道路中心点搜索提供了更加精确的概率估计。

1.5 中间点搜索算法

道路一般是直线或者缓慢变化的光滑曲线，而光滑曲线可以用直线段无限逼近，因此本文沿用并改进了文献[20]的迭代内插道路中间点搜索方案。首先种子点按用户点击顺序排列，连接每对相邻种子，在连接线的中垂线上搜寻最优的道路中心点，然后将该道路中心点视为新的人工种子点，重新搜索每对种子点间的最优道路中心点，直到所有相邻种子点的欧氏距离小于指定阀值为止，最后用直线段连接所有种子点，如图 4所示。

最优道路中心点的判定如式(8)—(12)所示，其中W为显著图，PL为中垂线，a、b为比例参数。

参数a、b在公式中具有明显的物理意义，即增大a表示目标之间的灰度相似更重要，增大b表示最佳道路中心点应更靠近首尾种子点的连线上。参数b有两个作用：第1，避免两条距离较近的道路在提取中间点的时候“串道”；第2，增加参数b的值，可以有效解决车辆、行道树、建筑物阴影等产生的干扰问题。算法的具体步骤如下：

输入：

种子点序列pts，显著图W，参数a、b，模板半径r，中间点间距阀值d。

处理：

步骤1——根据式(8)统计W_STD，设p₁、p₂为序列pts第1、第2个种子点。

步骤2——计算p₁p₂连线的中点p₀及p₁p₂的欧氏距离L，并更新种子点间距的最大值。

步骤3——在p₁p₂的中垂线上[-L/2~L/2]的范围内移动p₀的坐标，按式(9)计算L个测试模板内显著图平均值W_i，并计算对应的L个中心夹角θ_i=∠p₁p₀p₂。

步骤4——按式(10)计算L个测试模板的显著图差异最大值W_max。

步骤5——按式(11)计算L个测试模板的似然概率f_i。

步骤6——取L个测试模板的似然概率最大值，如式(12)，以该模板的圆心作为新的种子插入到pts，放在p₁、p₂之间。

步骤7——如果p₂为最后一个种子点，则转到步骤8，否则令p₁=p₂，p₂为下一个种子点，转到步骤2。

步骤8——如果种子点的间距最大值小于d，则转向输出，否则转到步骤2。

输出：

所有道路中心点序列。

本文的道路中心点搜索算法相比文献[18, 20]的算法具有如下3个优点：第1，搜索过程不再依赖形态学梯度图，而改为显著图，因道路两旁的平滑区域与路面的平滑区域在形态学梯度图中具有相似的特征，容易产生识别错误；第2，本文算法减少了一个算法参数，且参数的物理意义更加明显；第3，在一个路段提取过程中，如果存在大量中间点的W_i与W_STD差距较大，则可判定本次的中间点搜索失败，提醒用户介入，提高了算法实用性。

2 试验分析与对比

为了检验本文算法的性能和改进效果，分别采用文献[2]、文献[20]和本文算法对多张高分辨卫星遥感图像和航拍图像进行了试验和对比分析。

2.1 种子点敏感度试验

为了测试本文算法对用户种子点的敏感程度，对一张空间分辨率为0.3 m的航拍图像做了多次试验，结果如图 5所示。图 5(b)中首尾的白色圆形模板是自动生成的，无须人工设定阀值，黑色圆圈为迭代内插搜索得到的最佳中间点模板。从试验结果图可以看出，不论用户选点的初始位置如何偏离道路中心线，本文算法均能将模板调整到道路中心位置，说明本文算法对用户选点不敏感，用户只需将种子点落子道路上即可。

2.2 不同光谱道路的提取试验

为检验算法在不同光谱特性情况下的提取效果，分别对QuickBird分辨率为0.6 m的城市卫星图像和分辨率为0.5 m的乡村航空图像进行试验，效果如图 6所示。虽然不同类型传感器捕获的道路光谱特性不同，只要满足路面灰度平滑、道路目标与背景有明显的反差，就不影响形态学梯度图和显著图的生成，进而保证了算法的有效性。

2.3 曲折道路的提取试验

曲折道路提取是道路提取算法的难点，本文算法也不例外。图 6的试验结果表明，在弯曲度较小时，本文算法依然有效，但弯曲度一旦过大，将出现模板匹配错误，如图 7(a)所示。本文算法为了提高执行效率，在中间点搜索的过程中，限制了中垂线的总长度为相邻种子点的间距，如果弯道过于弯曲，真实道路将会出现在中垂线范围之外。如图 7(a)所示，其中白色曲线为错误提取的道路信息，虚线为起始点与终止点之间连线，点虚线为第一次迭代搜索中间点的中垂线范围，可以看出真实的道路中心点(白色圆圈)超出了中垂线搜索范围，造成模板匹配错误。增长中间点的搜索范围即可解决此问题，但必然会降低算法在平缓道路情况下的搜索效率，用户可在弯道处多设一个种子点来解决该问题，如图 7(b)所示。

2.4 噪声干扰试验

行道树遮盖、建筑物的阴影阻隔以及道路上的车辆，是构成遥感图像中道路提取的主要干扰信息，任何以局部灰度信息为主要依据的算法，均会受到不同程度的影响^[23]。本文算法综合利用了局部灰度信息和道路几何信息，外加人工介入辅助，可以较好地克服这些噪声的影响。图 8(a)—图 8(c)为3组噪声干扰试验，每组试验左侧为原始图，右侧为模板生成情况，白色圆圈代表算法根据人工种子点自适应生成的模板，黑色圆圈代表算法迭代内插搜索得到的最佳中间点模板。图 8(a)是一幅分辨率为0.6 m的QuickBird遥感图像，道路上分布大量的阴影噪声，算法参数为a=1、b=25。从图中可以看出，在加大参数b的比重后，迫使中间点夹角加大，将中间点“逼向”真实道路区域，但受阴影噪声影响，提取的道路信息稍微偏离了道路中心线。图 8(b)是一幅分辨率为0.6 m的pleiades-1遥感图像，道路上布满了车辆噪声，算法参数为a=1、b=1，可以看出路面车辆噪声对本文算法的影响较低。图 8(c)为一幅分辨率为0.6 m的Google地图，道路被大量的行道树覆盖。为了能够正确找到真实道路，在道路拐弯处多设了一个人工种子点，算法参数为a=1、b=10，说明绿化覆盖对算法影响较大，需要更多的人工介入才能解决。

2.5 改进效果试验

本文算法从多个方面对文献[20]的算法进行改进，为了检验本文算法的改进效果，用相同的人工种子点和模板尺寸，分别用文献[20]算法和本文算法对同一幅航空影像进行道路中心线提取，试验结果如图 9所示。图 9(a)为一幅分辨率为0.6 m的乡村道路航拍图像；图 9(b)是文献[20]算法的提取结果，算法参数α=1、β=1、γ=1；图 9(c)为本文算法的提取结果，算法参数a=1、b=1。文献[20]的算法出现了大量的模板匹配错误，因该算法在搜索道路中心点的过程中依赖于形态学梯度图，但本图像的道路两侧出现较多的平滑区域，且形态学梯度图对平滑区域的响应是相同的，从而无法区分不同类型的平滑区域。本文算法依赖显著图搜索道路中间点，图中道路和背景反差较大，算法很容易根据显著图识别真实的道路位置。试验结果表明本文算法具有更好的识别效果。

3 路网提取试验 3.1 图像去雾预处理

遥感图像通常受到空气中的细小颗粒、水滴等因素的影响，使得图像呈现“朦胧”的效果，从而影响算法的识别效率。图 10(a)是路网提取试验原始图像的一个局部，从中可以看出图像的雾化严重。为了得到更好的道路显著图，本文利用文献[25]的算法，对原始航拍图像进行去雾处理，效果如图 10(b)所示。视觉对比可以看出去雾后的图像使道路更加清晰，且道路与周围的其他物体反差更强。

3.2 路网提取结果

图 11是利用本文算法提取路网的试验结果。试验选取一幅地面分辨率为0.6 m(航摄比例尺1:15 000)的航拍图像，其中白色虚线为算法提取的道路，灰色的条状为人工提取的参考路网(含道路位置和缓冲宽度)，黑色小点为算法提取过程中人工设置的种子点。图中道路宽度不一、噪声类型多样、干扰程度迥异，可充分检验算法的优缺点。试验共分21个路段分别进行道路提取。对于被树木严重覆盖的路段，适当提高算法参数b的值，算法将提取到接近直线的道路信息，使得提取的道路信息可以穿越遮挡的异物；对于低噪声的路段，提高算法参数a的值，使之主要依靠显著图信息定位道路位置。算法最后对提取的曲线施加Bézier曲线平滑，使提取的结果更加贴近真实的道路形状。

3.3 算法性能比较

为了检验算法性能，通过人工提取参考路网，并按文献[23]定义如式(13)的评估参数，其中TP(true positive)、FP(false positive)、FN(false negative)分别表示正确提取、错误提取和未提取的道路。道路缓冲宽度不固定，以实际宽度为准(图 11灰色条形)，提取的道路坐标落在道路缓冲范围内即为正确。表 1给出了本文和文献[2, 20]的结果对比，其中总耗时包含人工选点、参数设置、算法识别以及部分路段返工操作的时间。从表中可以看出，文献[2]的算法虽然种子点数、参数数目最少，且完成度最高，但该算法在每个路段提取时都需要两次计算最短测地路径，耗时较长，且算法过度依赖图像的灰度信息，受限于试验素材的强干扰信息，正确率较低。文献[20]的算法需要设置更多的参数，且参数的配比难以人为判定，因此增加了返工次数，也增加了模板误匹配情况，降低了路网提取的正确率。本文算法在模板匹配过程中主要依赖道路显著图和道路几何信息，更容易根据路段的实际情况设置合理的参数配比。具体而言，对于具有大量干扰因素的平缓路段适当提高参数b的比重，对于平滑的弯曲路段适当提高参数a的比重，对于强噪声的弯曲道路需适当增加人工种子点，并提高参数b的比重。另外，由于模板具有一定的尺寸，造成路段首尾存在小段道路未能有效提取，因此本文和文献[20]的算法在完成度上均不如文献[2]，该问题可以通过自然延伸路段到图像边界或其他路段来解决。

3.4 路网准确性分析

从图 11的结果可以看出，大部分路段的提取结果均较好地贴近真实的道路中心线，但仍然存在部分路段偏离严重。如图 11下部的一个路段(标①处), 由于该路段宽度过大，且初始种子点处的路面灰度不均匀，造成自动生成的圆形模板无法完全覆盖整条道路宽度，从而使提取的结果与真实的道路中心线有较大的偏离。另外，由于本试验的素材存在严重的绿化噪声干扰，导致提取的部分路段出现一定偏差和震荡现象(标②处)。路网准确性问题的根本原因在于本算法主要根据图像的显著图来寻找道路中心点，在图像条件较好情况下，可以提取到较为精确的道路中心线，如图 6和图 9所示，而当道路灰度不均匀或异物干扰严重情况下将无法保证路网的准确性，通过适当增设人工种子点并提高参数b的比重，可在一定程度上解决中心线偏离或震荡问题，这也体现了半自动提取方式的优势所在。相较于文献[2]和文献[20]的道路提取结果(如图 12所示，分别用虚线和实线表示)，本文算法在路网准确性方面有一定程度的提高。

4 总结

模板匹配是高分辨率遥感图像中道路信息提取的常用手段，矩形模板及其变种因模板更贴近道路条状模型而更为常见，但矩形模板存在更多的人工介入、难以适应弯曲道路及抗干扰性差等缺点，而圆形模板在这些方面表现出更好的性能。本文提出了一种改良的圆形模板算法。与前人的圆形模板及矩形模板相比，该算法可根据形态学梯度图自动计算模板尺寸，同时在内插道路中间点搜索过程中引入显著图，使得算法更加有效，且算法参数更少，参数物理意义更直观，提高了算法的实用性。本文算法在道路弯度过大或干扰噪声过于复杂的情况下需要增加少量的人工干预，如何克服这方面不足，有待进一步的研究。