结合改进SWT和EMD的高分遥感影像道路提取

韦春桃，何蔚

(重庆交通大学 重庆智慧城市学院，重庆 400074)

0 引言

道路信息的提取与不断更新对于现代城市的发展具有重要意义。如今，基于遥感图像，尤其是高分辨率遥感图像的道路数据库已经广泛应用于车辆导航、城市规划等方面。

遥感影像道路提取方法主要可分为模板匹配、知识驱动、面向对象和深度学习四类[1]。连仁包等[2]提出了一种自适应圆形模板，结合形态学梯度图自动设定模板尺寸，利用道路显著图和几何形状信息搜索道路中心点。潘励等[3]发明了一种多层次知识驱动的全色遥感影像的道路变化信息提取方法，根据人的感知快速检测出道路的变化，将结果提供给用户，增强了地图修测的自动化程度。胡建青[4]基于易康软件，利用面向对象的方法，实现了对高分辨率遥感影像中的城市和乡村道路信息提取。戴激光等[5]针对神经网络训练存在分辨率降低和梯度消失，导致道路提取结果误提取率高和断裂问题，提出了一种基于多尺度卷积神经网络的遥感影像道路提取方法。

笔画宽度变换(stroke width transform，SWT)是Epshtein等[6]提出的一种用于自然图像文本检测的算法，它通过提取出具有一致宽度的带状目标来实现，因此也可以用于提取遥感图像中的带状地物。张霞等[7]考虑到SWT算法受周围地物影响较大的缺陷，提出了先将高分遥感图像通过均值漂移算法进行分割，再进行笔画宽度变换的框架，达到完整提取道路信息的目的。张祝鸿等[8]通过结合笔画宽度图与目标的几何特征，将高分二号近红外波段中的河流信息提取出来，较好地抑制了噪声。

陆地移动距离(earth mover’s distance，EMD)是Rubner等[9]提出的，从运输问题演化而来，是一种度量相似性的方法，可实现多种信息间的有效匹配[10]，对于遥感图像中的目标可以用作异质性分析。李建磊[11]利用EMD计算不同时相影像对象之间的颜色和边缘直线距离，为像斑的变化检测结果提供依据。Zhang等[12]提出了一个基于BOF模型的高分辨率遥感图像分类框架，通过EMD进行直方图匹配。

考虑到利用笔画宽度变换提取道路效果仍有待改进，本文提出了一种可自动求取最优阈值的笔画宽度变换方法，以优化道路的初步提取结果，并结合EMD进行多特征融合分析，对该方法进一步优化。

1 道路提取完整方法

首先，利用改进后的SWT初步提取有一定噪声的笔画宽度图，并确定道路参考区域；然后，利用连通区作为索引，通过计算光谱、纹理融合后的EMD值，分析道路参考区域与其他区域之间的特征相似度，设定合理的阈值过滤噪声，保留新的道路区域；最后，进行形态学后处理得到最终提取结果。总体技术流程如图1所示。

图1 道路提取流程

1.1 基于最优阈值的SWT

图2为整个SWT的主要过程。

图2 SWT过程

其步骤如下。

步骤1：利用Canny算子对灰度图像进行边缘检测，计算边缘图中每个边缘点的梯度值，得到梯度图。

步骤3：计算每个ray的长度，即ray上两个像素间的欧氏距离，并将该值分配给这个ray上的所有像素。该操作需另外建立一幅笔画宽度图像，并预先初始化(一般设为无穷大)所有像素值。

图3为高分影像中某段道路局部，在设定的不同最大笔画值作为阈值时获取的笔画宽度图。对于无阈值的笔画宽度变换，其结果依赖边缘检测，边缘点之间的连接长度没有限制，导致提取结果脱离了提取目标，出现了如图3(b)中红圈标记的大量长度过大的笔画，这对后续过滤噪声会有不利影响。但是，随意设定阈值同样无法得到最优粗提取结果。当阈值设置过小，会出现图3(c)中红圈标记的道路断裂情况，无法完整提取出来；当阈值设置过大，又会出现如图3(d)中红圈标记的多余笔画噪声；图3(e)为阈值合适时的SWT初步提取结果。

综上，有必要设计一个求取最佳阈值的方法，使得笔画宽度图为最优提取结果，为后续的精提取奠定基础。

图3 不同宽度阈值下的笔画宽度图

图4 道路参考区提取方法

自动求取最优阈值的具体步骤如下。

步骤1：在无阈值的情况下进行笔画宽度变换，得到具有最多噪声的笔画宽度图。以图3(b)为例，对该图进行连通区标记，图4(a)中每个红色矩形框内表示一个连通区。

步骤2：对笔画宽度图进行统计，并确定其峰值对应的笔画宽度，如图4(b)所示。通过图4(c)得到其峰值处笔画宽度占比最大的连通区，并将该区域作为道路参考区域，如图4(d)所示。

步骤3：以道路参考区作为对象，不断初始其最大笔画宽度值并计算其连通区域数量，以连通区域数量N发生改变(8邻域)作为道路断裂的判断依据，视其为迭代终止条件，输出此时的最大笔画宽度M，即最优阈值。

求取阈值流程图如图5所示。

图5 最优阈值确定方法

1.2 EMD

在进行SWT初步提取后，通过1.1节中的方法对该结果再次标记连通区并确定道路参考区域。通过计算特征距离，反映道路参考区与其他区域之间的特征相似度(图6)，从而确定新的道路区域以及噪声(非道路)区域。

图6 特征相似度示意图

EMD作为用来计算特征相似度而存在的特征距离，相比于范式距离(如曼哈顿距离和欧氏距离)，可以更准确地量化直方图元组之间的相似性。设现在需要求直方图P={(p1，wp1)，…，(pm，wpm)}和Q={(q1，wq1)，…，(qn，wqn)}之间的相似度，其中pi和qj分别表示两个直方图不同的数值，wpi和wqj则分别表示它们所占的权重，则可以先求得一个距离矩阵，如式(1)所示。

D=[dij]

(1)

式中：dij表示P、Q中任意两项pi和qj之间的代价(距离)。

此时，从P转化为Q的工作量可以用WEMD表示(式(2))。

(2)

F=[fij]

(3)

式中：fij表示从pi到qj的权重流动量。EMD的本质就是找到一个最佳F，使得WEMD最小，它的实质是一个线性规划问题，如图7所示。对于以上公式有以下四个需要注意的限制条件。

1)fij不能小于0，即fij≥0。

图7 EMD示意图

1.3 基于直方图提取光谱、纹理特征

在进行特征相似度比较前，需先提取出适用于后续基于EMD识别道路目标的特征，本文选择通过光谱、纹理直方图特征进行表达。

光谱特征的表达对于遥感影像不同地物的识别有较大影响，颜色直方图一定程度上能反映地物光谱信息的分布情况，因此可作为光谱特征进行提取[13]。本文选择在HSV颜色空间中提取直方图，该空间颜色分为H、S、V三个分量(色调、饱和度、明度)，相对于常用的RGB颜色空间，它的色相和饱和度接近于人类的视觉系统。这里将遥感影像的H、S、V三个分量提取出来进行量化，通过构造HSV累计颜色直方图表达影像光谱特征。

纹理特征主要用来描述地表的空间分布状态和粗糙程度，同样也是判别遥感图像中地物的重要标志[14]。局部二值模式(local binary pattern，LBP)[15]是一种用来描述图像局部纹理特征的算子，通过比较中心像元与邻域像元的灰度值表达纹理，计算简单，能够清晰表达纹理。通过对获取的LBP图谱进行直方图统计，可得到影像的LBP纹理特征。图8为3×3窗口的LBP值计算示例。其中，pi和pc分别表示邻域像素值和中心像素值，T的取值由pi与pc的差值确定。其主要原理是确定一个中心像素，比较一定窗口内相邻像素与中心像素的大小，灰度值大于中心像素的设为1，否则为0。

图8 3×3窗口LBP计算示例

Ojala等[16]又提出了一种旋转不变均匀的LBP算子，它相对于传统的LBP算子包含了更多的局部纹理特征，并且减轻了由图像旋转引起的LBP值改变，此处采用改进的LBP算子提高去噪的精度。

1.4 基于多特征融合EMD去噪

(4)

(5)

(6)

(7)

式中：φ1和φ0分别表示判定为保留区域(道路)和删除区域(噪声)。

1.5 形态学后处理

根据EMD 分析去噪后可去除掉大部分非道路区域，完整保留道路区域。将此时的道路提取结果转为二值图像，在不损坏提取出的道路信息情况下，通过形态学后处理可去除剩余的前文方法无法识别的其他噪声，针对道路区域也可以有效去除其毛刺并且填充孔洞部分，达到优化其最终成果的目的。具体的形态学操作如图9 所示，采用边长为3 的方形结构元。

图9 形态学后处理

2 实验结果与分析

2.1 实验数据及参数设置

实验使用分辨率为0.8 m的高分二号影像。获得多幅不同场景影像，道路的分布、形态也有所不同，如图10(a)～图10(e)所示。其中，道路1(359像素×361像素)为山区道路，其光谱属性明显区别于背景，亮度较高，弯曲程度较大，但周围有较多与其光谱特征相似的细小岩石；道路2(444像素×392像素)、道路3(788像素×387像素)为耕地间道路，相对其他场景道路较窄，且与耕地的纹理分布有一定相似性；道路4(251像素×210像素)为厂区道路，其周围地物相对较多，有些许建筑物稀疏分布，纹理相对道路1、2、3复杂，但道路与其他地物的特征差别较大；道路5(244像素×252像素)为住宅区道路，房屋分布密集，其他地物大多与道路相连，且与道路间有相似度较高的纹理及光谱特征，因此其提取道路难度较高。

图10 道路原始影像

由于实现本文方法及对比方法的实验需要，确定了 1.1节中改进SWT自动求取的阈值参数M，并人工设定了1.4节中基于多特征EMD去噪的阈值参数λ，具体数值如表1所示。其中，通过改进SWT操作后的各道路笔画宽度如图11所示。

表1 道路提取参数

图11 基于改进SWT的笔画宽度图

2.2 本文方法实验结果与分析

在默认都使用形态学后处理的情况下，使用本文方法与仅使用SWT提取道路的结果对比如图12所示。由图12(a1)～图12(e1)可以看出，SWT算法可以较为完整地提取道路信息，但同时也提取出了与道路一样宽度变化小的其他地物。对于如道路2和道路3的田间道路，有明显与道路几何特征类似的区域被同时提取；对于道路4和道路5，容易提取出排列、形状较为规律的房屋。由图12(a2)～图12(e2)可以看出，本文方法在完整提取道路信息的同时，也结合纹理和光谱特征，滤除了其他通过SWT提取出的多余区域，说明本文方法不仅成功将SWT运用在提取道路方面，还弥补了SWT无法高准确性提取道路的不足。

图12 SWT与本文方法提取结果

实验采用Bowyer等[17]提出的道路提取精度评价方法，对本文与其对比的其他方法进行定量评估。评估标准包括提取出道路的完整率Cp、准确率Cr、提取质量Ql三项。

表2 SWT与本文方法的提取精度对比 %

对于SWT与本文方法精度定量评估结果如表2所示。由表2可以看出，本文方法提取道路的完整率都要稍低于SWT，但都控制在3%以内，对于提取难度较小的道路1，提取完整率可达到97.507%，对于提取难度较大的道路4、5，提取完整率也高于85%。但SWT方法提取道路的准确率却远小于本文方法，尤其在住宅区道路5的准确率只有44.068%，而使用本文方法道路1、2、3、4则保持准确率在90%以上，道路5的提取准确率也能达到72.228%，因此本文方法的提取质量均要高于SWT方法。

2.3 对比方法实验结果与分析

实验选择使用文献[7]中提到的将中值滤波、均值漂移算法与SWT相结合的方法，以及利用巴特沃斯高通滤波[18]图像增强的方法，与本文方法提取道路的结果进行比较，提取结果如图13所示。

对比图13与图12可以看出，文献[7]和文献[18]相对于SWT都能够更少提取出非道路区域，两种方法的提取结果都存在少许的多余地物，但文献[18]由于可以突出图像中模糊的边缘轮廓，使线条变得清晰，因此明显比文献[7]保留了更多的道路信息，文献[7]针对SWT受Canny边缘检测结果影响较大的局限性，加入中值滤波及均值漂移分割算法，减少了大量非道路的细碎区域。相比之下，本文方法与文献[18]一样也保留了大量的道路信息，但其不仅减少了多余的细碎区域，也删除了大面积的其他地物，尤其体现在道路3及道路4的提取结果。

图13 文献[7]和文献[18]提取结果

文献[7]和文献[18]方法提取结果精度定量评价结果如表3所示。对比表2和表3可以看出，对于轮廓清晰，与背景光谱及纹理差异较大的道路1，本文方法与对比方法的提取完整率都可以达到96%以上；对于其他道路，文献[18]与本文方法的提取完整率相差不多，但文献[7]相对较低。除道路3以外，文献[7]的提取准确率均高于文献[18]，但两种对比方法对道路2、3、4、5的提取准确率均小于85%，远低于本文方法。使用文献[7]方法提取道路4、5，提取质量仅为58.639%、54.864%；使用文献[18]方法提取道路5，提取质量仅为55.533%；而对比方法提取道路2、3的提取质量均低于75%。综上，使用本文方法提取道路，可以在保证道路信息完整的同时，大量减小错提概率，得到质量较好的道路区域。

表3 文献[7]和文献[18]高通滤波提取精度对比 %

3 结束语

本文提出了一种结合改进后SWT和多特征EMD的道路提取方法。通过实验结果的定性与定量分析，本文方法表现出了较好的效果。与深度学习方法相比，此方法的普适性更高，可用于各种不同传感器的遥感影像及不同场景中道路识别，并且无需大量训练及测试样本准备，人工干涉量小。但由于SWT的效果本身依赖边缘检测结果，因此在应用于城区、住宅区等有阴影、遮挡物且边缘模糊场景中的道路提取时，精度相对较低，在实验结果中也有所体现。另外，考虑到要计算道路参考区与其他区域之间特征的EMD值，选用影像场景中的道路材质需要一致或者特征相似，因此目前本文方法一定程度上限制了场景大小。实验中选取的EMD阈值是建立在经验总结的情况下，后续研究会找到科学、严谨的确定阈值方法。