0 引言

1 研究方法

1.1 技术路线

首先对输入的原始遥感影像进行预处理，然后分析选定像元的角度纹理图，找出其灰度方差值最小的方向，并依据该方向设计对应的Gabor滤波器组，用滤波输出值作为分割对象结合k均值聚类分割策略进行分割，最终依据真实道路矢量数据对提取效果进行精度评价，研究步骤如图 1所示。

图 1 研究步骤

Fig. 1 The flow chart of study

1.2 影像预处理

1.3 角度纹理特征的获取

高分辨率遥感影像中，道路通常具有固定像素的宽度且其灰度常常沿道路方向分布均匀。利用这一特征，可以构造一个矩形的纹理滤波模板，对每个像元，让模板的几何中心与其重合并在[0°，360°]内旋转，每隔一定角度输出该模板区域内像素的灰度值均方差^[19-20]。

如图 2(a)中所示，用一个黄色的矩形模板使其中心点与某条南北向主干道的中心点重合，然后该模板在[0°，360°]内以15°为步长旋转(见图 2(b))，分别输出24个角度所对应的灰度值方差，由此得到一个角度纹理图，如图 3所示。其横坐标表示方向号，纵坐标表示灰度方差值。

图 2 矩形模板示意图

Fig. 2 Rectangle template chart

通过分析角度纹理图，可以判断其对应的像元点是否为可能的道路点。在角度纹理图中，每个灰度方差最小值所对应的方向(谷特征)通常都暗示着该方向可能存在道路。谷特征越明显，其所对应的方向存在道路的可能性越大。如果存在两个明显的谷特征，则表明该点位于连续道路段上。多数情况下，两个明显谷特征对应角度相差180°(如图 3中所示)。若存在3个或4个谷特征，则该点为道路交叉口。对每个像元点来说，角度纹理图提供两类信息：一是每个像元点对应有多少谷特征，二是谷特征所对应的方向。

图 3 角度纹理图

Fig. 3 Angle texture chart

运用角度纹理特征可以有效地将道路与建筑物、植被、水体等其他地物分离，从而获得道路网。Haverkamp等人^[19]设计了一个3×20的模板，用其找出所有具有2个或2个以上谷特征的像元点并将它们连通。如果存在只有一个谷特征的像元点且其谷特征方向与多谷特征点之间的连通方向一致，则该像元点被认为是道路点。该方法仅仅通过是否存在谷特征来判别一个像元点是否为道路点，这一直观的准则通常难以与其他基于影像特征的分割方案有机结合，但仍不失为为一种有意义的对纹理特征的描述。从中借鉴一个重要思路：通过比较不同方向上的灰度值方差，可以确定道路的方向并在该方向上构造用以分割道路像元点的特征量。

利用角度纹理特征的两个重要参数是矩形模板的大小及角度间隔。很多学者采用文献[19]中的模板。虽然道路在高分辨影像中具有部分面状地物的特征，但城市中小道路的宽度通常只有几个像素，如果采用宽度为3的模板，小道路往往难以分割出来，主道路的边缘易与周围的地物粘连。为此，模板的宽度通常应与拟提取的最窄道路宽度像元数相当，模板长度取最大道路宽度的1.5 2倍之间的某个值。

1.4 Gabor滤波器的设计

为了研究信号在局部范围内的频率特性，Gabor^[21]于1946年在高斯函数的基础上提出了著名的“窗口”傅里叶变换，即Gaobr变换。Gabor变换具有可以调制的中心频率和方向，由其衍生的Gabor滤波器在分析数字图像局部区域的频率特征和方向特征方面具有极其优异的性能。Daugma ^[22]将Gabor变换扩展到2维结构，其表达式为

$ \begin{array}{l} \begin{array}{*{20}{l}} {h\left( {x,y} \right) = \frac{1}{{2{\rm{ \mathsf{ π} }}{\sigma ^2}}}{\rm{exp}}\left[ { - \frac{{{x^2} + {y^2}}}{{2{\sigma ^2}}}} \right] \times } \end{array}\\ {\rm{exp}}[2{\rm{ \mathsf{ π} }}j(Ux + Vy)] \end{array} $

(1)

式中，$\left( U,V \right)$表示特定的空间频率，$g\left( x,y \right)$表示下式所示的高斯函数。由于滤波器的方向已经确定(即1.3节中最小灰度方差所在方向)，由式(1) 可知，如果要确定一个Gabor滤波器，只需确定其频率参数$\left( U,V \right)$和高斯参数$\sigma $。根据Gabor滤波器的带通特性，如果纹理能量集中的频率范围与Gabor滤波器的通频带吻合，则滤波器的输出值会很大，否则其输出将受到抑制。由此，对道路和非道路纹理进行傅里叶变换，找出区分它们的频率范围并设计出对应的Gabor滤波器。Gabor滤波器的半峰带宽$B$ (倍频程)为

$ B = {\rm{lo}}{{\rm{g}}_2}\left( {\frac{{{\rm{ \mathsf{ π} }}F\sigma + \alpha }}{{{\rm{ \mathsf{ π} }}F\sigma - \alpha }}} \right) $

(2)

式中，$F=\sqrt{{{U}^{2}}+{{V}^{2}}}\rm{ }$为Gabor滤波器的中心频率，$\alpha =\sqrt{\frac{(\rm{ln}2)\rm{ }}{2}}$。为了与人类视觉神经细胞的频率响应范围相吻合^[23]，通常取$B$的值为1，由此可以导出$\sigma =\frac{3\alpha }{\pi F}$。可见此时只需要确定$\left( U,V \right)$即可。假设输入的影像为$I\left( x,y \right)$，确定$\left( U,V \right)$具体步骤如下：

1) 用设计好的矩形模板随机选取大量的道路与非道路地物样本。注意，在选取道路地物样本时，矩形模板的方向尽量与道路中心线方向平行。如图 4中黄色所示。

图 4 选取不同地物的样本示意图

Fig. 4 Chart of selecting different samples

2) 假设道路类样本为${{I}_{r}}\left( x,y \right)$，非道路类样本为${{I}_{u}}\left( x,y \right)$。分别对它们进行傅里叶变换得到它们的功率谱密度${{P}_{r}}\left( u,v \right)$和${{P}_{u}}\left( u,v \right)$，其对应的极坐标形式为${{P}_{r}}\left( f,\theta \right)$和${{P}_{u}}\left( f,\theta \right)$。

3) 根据式(3) 找出两类纹理平均功率谱密度相差最大的3个频率点$\left( {{F}_{i,}}{{\theta }_{i}} \right)$，$i$ =1，2，3并将其转换为对应的频率参数$\left( {{u}_{i,}}{{v}_{i}} \right)$ =1，2，3。

$ {\left( {F,\theta } \right) = {\rm{arg}}\left\{ {{\rm{max}}\left[ {\frac{{{P_u}\left( {f,\theta } \right)}}{{{P_r}\left( {f,\theta } \right)}}} \right]} \right\}} $

(3)

4) 将得到的$\left( U,V \right)$和$\sigma $代入式(1)，则可以得到最佳的Gabor滤波器组。

1.5 特征向量的构造

采用设计出的Gabor滤波器组对原始遥感影像进行卷积滤波。如图 2所示，让Gabor滤波器模板的基准点对准当前像元点，然后以15°为步长绕当前点旋转模板，在每个角度位置上，计算模板内系数与下面像素灰度的乘积之和。3个模板顺次在24个方向上进行卷积滤波运算，每个像素得到72个滤波响应值。用$O_{s}^{ij}$表示第$i$个模板在$j$方向上对像元点$s$的滤波输出值。同时用${{O}_{\rm{max}}}$和$ {{O}_{\rm{min}~}}$分别表示所有滤波输出值中的最大值和最小值。为了保证图像的稳定分割，对滤波输出值进行归一化，即

$ \begin{align} &\left( O_{s}^{^{ij}} \right)\prime =\frac{O_{s}^{^{ij}}-{{O}_{\rm{min }\!\!~\!\!\rm{ }}}}{{{O}_{\rm{max}}}-{{O}_{\rm{min}}}} \\ &i\in \left\{ 1, 2, 3 \right\};j\in \{1, 2, 3, \ldots, 24\} \\ \end{align} $

(4)

接着，对同方向的3个归一化滤波值和归一化的影像灰度值 ${{g}_{s}}$ 进行加权，即

$ h_{s}^{j}=2\times {{g}_{s}}+\left( O_{s}^{1j} \right)\prime +\left( O_{s}^{2j} \right)\prime +\left( O_{s}^{3j} \right)\prime $

(5)

式中，加权系数2、1、1、1来自经验数据。在构造特征向量时，$j$通常取24方向中对应的$h_{s}^{j}$方差最小的那一个。即使是非道路点，该方向亦有较好的一致性。最终用以分割的特征向量为

$ \left[ {{\sigma _s}_{{\rm{min}}},\frac{{O_s^{1j}}}{{{g_s}}},{\rm{ }}\frac{{O_s^{2j}}}{{{g_s}}},{\rm{ }}\frac{{O_s^{3j}}}{{{g_s}}}} \right] $

(6)

式中，${{\sigma }_{s}}_{\rm{min}}$表示的$h_{s}^{j}$最小方差，$h_{s}^{1j},h_{s}^{2j},h_{s}^{3j}$除以自身灰度值的目的是为了消除灰度值大小对分割结果的影响。

文献[19]计算最小灰度方差，本文方法计算加权数据最小方差，对这两个最小方差可以通过一个简单实验做个对比。以图 4作为原始影像，图 5给出了文献[19]方法和本文方法计算得到的最小方差。显然，图 5(a)中道路与周边地物存在着一定混淆，而图 5(b)中区分较为明显。

图 5 不同方法方差对比图

Fig. 5 Chart of Variance contrast((a)reference[19]; (b)ours)

1.6 k均值聚类分割

1.7 精度评价

2 实验结果与讨论

2.1 研究数据及区域

本文研究数据是高分一号全色单波段遥感影像，成像时间为2015年4月20日。高分一号卫星是中国高分辨率对地观测系统的第一颗卫星，它的特点是增加了高分辨率多光谱相机，全色分辨率是2 m，多光谱分辨率为8 m，该相机的性能在国内投入运行的对地观测卫星中最强。“高分一号”在具有类似空间分辨率的同时，可以在更短的时间内对一个地区重复拍照，其重复周期只有4 d，而世界上同类卫星的重复周期大多为10余天。可以说，高分一号实现了高空间分辨率和高时间分辨率的完美结合。使用高分一号数据并辅以适当的道路提取算法，就不难实现对所观测区域路网信息的实时更新。

本文选择的研究区域位于郑州市中心区(图 6)。该区域路网密集，主干道路轮廓清晰、干扰较少，是较为理想的研究区域。

图 6 原始影像

Fig. 6 The original image

2.2 实验结果

图 7中白色部分为提取出来的主干道路面状矢量，为了方便比较与原始影像进行了叠加。

图 7 城市主干道路提取结果

Fig. 7 The result of urban trunk road extraction

2.3 精度评价

图 8(a)为河南省测绘地理信息局提供的郑州市道路矢量与原始影像的叠加，图 8(b)是借助ARCGIS的中心线提取功能将图 7中的面状道路矢量转换为道路中心线并与原始影像叠加后的结果。图 8(c)—(f)分别表示通过表 1中方法2—方法6获取道路矢量后提取中心线并与原始影像的叠加。

图 8 不同方法获取的道路中心线矢量

Fig. 8 Road centerline vector obtained by different methods

((a)ground truth; (b) ours; (c)method two; (d)method three; (e)method four; (f)method five)

表 1 城市主干道路提取精度评价
Table 1 Quantitative evaluation of urban trunk road extraction

下载CSV

方法	方法名称	完全率	正确率	提取质量
1	本文方法	0.72	0.78	0.65
2	分割对象为灰度值	0.45	0.51	0.37
3	分割对象为灰度共生矩阵	0.62	0.70	0.54
4	分割对象为多通道2维Gabor滤波器	0.58	0.66	0.52
5	Hough变换法	0.41	0.56	0.34
6	多尺度分割法	0.66	0.73	0.58

图 8中所有道路中心线均以红色线状矢量表示。

表 1利用1.7节的评价指标对本文方法及引言提到的其他几种方法进行了比较，列表总结如下：方法1—方法4都以k均值聚类法作为分割策略。方法2的分割对象为像元灰度值^[13-14]；方法3以共生矩阵生成的对比度、熵和相关性为分割对象，采用5×5窗口，步长取1，4个方向分别为0°、45°、90°、135°^[15-16]；方法4以多通道的2维Gabor滤波结果作为分割对象，使用3个频率、6个方向^[18]。Hough变换法是一种典型的基于线状特征提取方法，该方法首先对影像二值化，然后对非道路信息去噪，最后用Hough变换的方法检测提取道路^{[1, 2]}。由于难以提取弯曲的边缘信息，该方法的效果较差。多尺度分割法是基于面状特征提取方法中极具代表性的一种方法，首先使用邻域均值作为阈值二值化，影像以改进一般LBP算子；再根据局部方差与尺度间的关系，找到适合道路LBP纹理值的最佳尺度；最后利用区域生长和形状指数提取最佳尺度下的影像道路信息，并通过尺度收缩关联到原始影像^[7-8]。从表 1中不难看出，在使用相同分割策略的前提下，使用纹理值作为分割对象其分类精度比使用灰度值要高。如果按照本文方法，沿着像元点最小灰度方差的方向来构造滤波器并以其输出值作为分类对象，则可以获得更好的分类精度。由于Gabor滤波器优越的带通特性，较于Hough变换法和多尺度分割法，本文方法体现出更好的分类效果。

3 结论

参考文献

[1] An L.Extraction of city road from high resolution sensing images based on high transform[D]. Nanchang:East China Institute of Technology, 2015. [安丽. 基于Hough变换的高分辨率遥感影像道路提取[D]. 南昌: 东华理工大学, 2015.]
[2] Yu F, Yu J, Jian X. Study on the attraction methods of high precision remote sensing road images based on Hough transformation[J]. Railway Investigation and Surveying, 2011, 37(5): 42–45. [余峰, 余洁, 简骁. 基于Hough变换的高分辨率遥感影像道路提取方法研究[J]. 铁道勘察, 2011, 37(5): 42–45. ] [DOI:10.3969/j.issn.1672-7479.2011.05.013]
[3] Fu X Y, Zhang F L, Wang G J, et al. Road extraction from SAR images using tensor voting and Snakes model[J]. Journal of Image and Graphics, 2015, 20(10): 1403–1411. [符喜优, 张风丽, 王国军, 等. 结合张量投票和Snakes模型的SAR图像道路提取[J]. 中国图象图形学报, 2015, 20(10): 1403–1411. ] [DOI:10.11834/jig.20151014]
[4] Tang W, Zhao S H. Road extraction in quaternion space from high spatial resolution remotely sensed images basing on GVF Snake model[J]. Journal of Remote Sensing, 2011, 15(5): 1040–1052. [唐伟, 赵书河. 基于GVF和Snake模型的高分辨率遥感图像四元数空间道路提取[J]. 遥感学报, 2011, 15(5): 1040–1052. ] [DOI:10.11834/jrs.20110239]
[5] Cao X G.Research on unstructured road recognition and fitting algorithm[D]. Nanchang:East China Institute of Technology, 2015. [曹旭光. 非结构化道路识别及拟合算法研究[D]. 南昌: 东华理工大学, 2015.]
[6] Xu G, Zhang D W, Liu X Y.Road extraction in high resolution images from Google earth[C]//Proceedings of IEEE International Conference on Information, Communications and Signal Processing.Macau:IEEE, 2009:1-5.[DOI:10.1109/ICICS.2009.5397470]
[7] Huang Z J.Research on multiscale methods in object-based image analysis[D]. Changsha:National University of Defense Technology, 2014. [黄志坚. 面向对象影像分析中的多尺度方法研究[D]. 长沙: 国防科学技术大学, 2014.]
[8] Chen J, Yang M H, Zang Z, et al. Research of improved automatic high resolution image extraction method based on LBP operator and multi-scale analysis[J]. Remote Sensing Information, 2010(2): 3–7, 11. [陈杰, 杨敏华, 臧卓, 等. 基于LBP算子与多尺度分析的高分辨率遥感影像道路自动提取方法改进研究[J]. 遥感信息, 2010(2): 3–7, 11. ] [DOI:10.3969/j.issn.1000-3177.2010.02.001]
[9] Yin W T.Research on high resolution remote sensing image building extraction based on the multi-kernel support vector machine[D]. Ganzhou:Jiangxi University of Science and Technology, 2015. [尹文亭. 基于多核支持向量机的高分辨率遥感影像建筑物提取研究[D]. 赣州: 江西理工大学, 2015.]
[10] Wang M, Luo C J, Zhou C H, et al. Extraction of road network from high resolution remote sensed imagery with the combination of Gaussian Markov random field texture model and support vector machine[J]. Journal of Remote Sensing, 2005, 9(3): 271–276. [汪闽, 骆剑承, 周成虎, 等. 结合高斯马尔可夫随机场纹理模型与支撑向量机在高分辨率遥感图像上提取道路网[J]. 遥感学报, 2005, 9(3): 271–276. ] [DOI:10.3321/j.issn:1007-4619.2005.03.008]
[11] Lin P, Ruan R Z, Wang Y Q, et al. Research on extraction of road based on object oriented in an urban context[J]. Geography and Geo-Information Science, 2016, 32(1): 49–55. [林鹏, 阮仁宗, 王玉强, 等. 一种基于面向对象的城镇道路自动提取方法研究[J]. 地理与地理信息科学, 2016, 32(1): 49–55. ] [DOI:10.3969/j.issn.1672-0504.2016.01.010]
[12] Chen C, Qin Q M, Chen L, et al. Extraction of bridge over water from high-resolution remote sensing images based on spectral characteristics of ground objects[J]. Spectroscopy and Spectral Analysis, 2013, 33(3): 718–722. [陈超, 秦其明, 陈理, 等. 基于地物光谱特征分析的高分辨率遥感图像水上桥梁提取[J]. 光谱学与光谱分析, 2013, 33(3): 718–722. ] [DOI:10.3964/j.issn.1000-0593(2013)03-0718-05]
[13] Shi W Z, Zhu C Q. The line segment match method for extracting road network from high-resolution satellite images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(2): 511–514. [DOI:10.1109/36.992826]
[14] Zhang Q P, Couloigner I. Benefit of the angular texture signature for the separation of parking lots and roads on high resolution multi-spectral imagery[J]. Pattern Recognition Letters, 2006, 27(9): 937–946. [DOI:10.1016/j.patrec.2005.12.003]
[15] Song M J, Daniel Civco. Road extraction using SVM and image segmentation[J]. Photogrammetric Engineering & Remote Sensing, 2004, 70(12): 1365–1371. [DOI:10.14358/PERS.70.12.1365]
[16] Wu J, Zhao Z M. Scale co-occurrence matrix for texture analysis using wavelet transform[J]. Journal of Remote Sensing, 2001, 5(2): 100–103. [吴均, 赵忠明. 利用基于小波的尺度共生矩阵进行纹理分析[J]. 遥感学报, 2001, 5(2): 100–103. ] [DOI:10.11834/jrs.20010205]
[17] Weldon T P, Higgins W E.An algorithm for designing multiple Gabor filters for segmenting multi-textured images[C]//Proceedings of International Conference on Image Processing, 1998.ICIP 98.Chicago, IL:IEEE, 1998:333-337.[DOI:10.1109/ICIP.1998.999023]
[18] Jain A K, Farrokhnia F. Unsupervised texture segmentation using Gabor filters[J]. Pattern Recognition, 1991, 24(12): 1167–1186. [DOI:10.1016/0031-3203(91)90143-S]
[19] Haverkamp D S. Extracting straight road structure in urban environments using IKONOS satellite imagery[J]. Optical Engineering, 2002, 41(9): 2107–2110. [DOI:10.1117/1.1496785]
[20] Haverkamp D S, Poulsen R.Complementary methods for extracting road centerlines from IKONOS imagery[C]//Proceedings of the SPIE Volume 4885, Image and Signal Processing for Remote Sensing Ⅷ.Crete, Greece:SPIE, 2003, 4485:501-511.[DOI:10.1117/12.468220]
[21] Gabor D. Theory of communication[J]. Proceedings of the Institute of Electrical Engineers, 1946, 93(26): 429–457.
[22] Daugman J G. Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by two-dimensional visual cortical filters[J]. Journal of the Optical Society of America A-Optics Image Science and Vision, 1985, 2(7): 1160–1169. [DOI:10.1364/JOSAA.2.001160]
[23] Pollen D A, Ronner S E. Visual cortical neurons as localized spatial frequency filters[J]. IEEE Transactions on Systems, Man, and Cybernetics, 1983. [DOI:10.1109/TSMC.1983.6313086]
[24] Wiedemann C, Heipke C, Mayer H.Empirical evaluation of automatically extracted road axes[C]//Proceedings of 9th Australasian Remote Sensing and Photogrammetry Conference.Sydney:IEEE, 1988:1-16.