基于频域特征的遥感图像城市道路绿地覆盖轮廓提取

赵好好，冯学智，肖鹏峰

(1.南京大学 地理信息科学系，南京 210023；2.南京信息工程大学 遥感学院，南京 210044)

1 引 言

道路绿地是城市绿地的重要组成部分，具有改善城市生态环境，调节小气候，提高城市居民生活舒适度等功能，城市道路绿地覆盖信息的提取对城市建设和管理具有重要意义。随着卫星遥感图像空间分辨率的提高，基于像元的遥感图像分类技术开始向基于对象的特征识别和提取技术过渡，这使得城市地物特征信息的提取成为可能[1]。

遥感图像城市绿地信息识别和提取研究中，光谱特性是主要的研究方法。植被指数是早期绿地特征识别的主要手段，其中NDVI应用最为广泛。也有学者则利用光谱混合分析进行城市绿地覆盖类型的识别[2-3]。随着图像分辨率的提高，地物纹理、轮廓、形状等特征得到重视，面向对象的方法在城市绿地信息提取中开始广泛应用，基于像元的分类方法逐渐向基于特征的提取技术过渡[4]。

频谱特征是进行图像识别的重要依据，通过建立空间纹理结构特征的频谱分析模型来计算不同大小地物类型的纹理频谱，可有效提高图像地物识别能力[5]。通过对不同地物进行频谱特征分析，建立地物频谱识别标志，可以有效提高地物信息提取的精度和自动化水平[6]。

在频率域遥感图像分析中，王珂等对城市河道线性信息进行了识别和提取[7]，吴桂平等对城市人工建筑进行了识别[8]。而城市道路绿地覆盖，与城市道路之间具有特殊的几何关系，本文研究在几何约束条件下的城市地物(道路绿地覆盖)频域信息提取方法。通过对高分辨率遥感图像进行傅里叶变换，建立对应道路绿地覆盖轮廓的截止频率，设计合适的滤波器，并借助道路中心线几何辅助信息，实现道路绿地覆盖的信息提取。

2 研究数据

本文所使用的高分辨率遥感数据为南京市主城区2007年7月获取的QuickBird图像，包含4个多光谱波段和1个全色波段，融合后的空间分辨率为0.61m。图像成像质量良好，无云层覆盖。

图1(a)为研究区近红外波段图像，位于中山北路与云南路交叉地带，图像大小为558×558像素。道路绿地与道路方向一致，在图像中沿着两个方向分布，且基本垂直。在几何关系上，道路绿地位于道路红线以内，如果道路中心线能够确定，做道路中心线的缓冲区，可提取道路绿地覆盖的位置信息。

图1(a)中，近红外波段道路绿地与道路的光谱差异显著，道路绿地的亮度值平均为650，而道路亮度均值为310。道路信息包含道路色调信息和道路边缘信息。本实验区中，由于道路被行道树覆盖比较严重，所以道路边缘信息不显著。而道路的色调信息在图像上亮度最低，与建筑、行道树差别显著，较易区别。

图1 研究区遥感图像及频谱特征

3 研究方法

3.1 频谱分析

假设一幅大小为M×N的遥感图像，其二维离散傅里叶变换如下式表示：

(1)

其中，x，y为空域坐标，u，v为频率。图像的频谱值可表示为|F(u,v)|。

频谱分析通常在极坐标系中进行，假设极坐标系频谱函数为S(r,θ)，r和θ是坐标变量。通过辐射扫描法可分析径向分布特征和角向分布特征，如下式所示：

(2)

(3)

研究区遥感图像的频谱图如图1(b)所示。从频谱图上可以看到，图像的中高频丰富，在与道路绿地垂直的方向上，存在两条显著的谱线。另外在水平与竖直方向上，也存在少量谱线，这与图像中的建筑等地物有关。图1(c)径向分布曲线表明，随着频率的增加，频谱逐渐降低，没有出现显著的谐波成分。在频率20周/像元，如图中箭头所示，能量下降的幅度出现了显著转折。图像的色调信息主要集中在直流中心，20周/像元反映了低频位置。图1(d)角向分布曲线中，出现了两个峰值能量方向，分别为45°和135°，两个方向均垂直道路绿地的方向。图像的径向分布和角向分布主要反映了在频率和方向上的能量累积值，而与道路绿地覆盖轮廓特征相关的频率细节信息则无法体现。

3.2 特征识别

为了详细分析对应道路绿地覆盖轮廓的截止频率，取图1(a)中小区域256×256尺寸的图像T1，如图2(a)所示。图2(b)为其频谱图，从角向分布曲线图2(c)可以发现，其峰值角向能量为135°，与道路绿地方向垂直。现做频谱图135°方向的剖面线，因频谱图的对称性质，这里仅取其一半显示。如图2(d)圆圈所示，在频率20、37、47周/图像出现了显著的波谷位置。频谱曲线中的波谷位置，反映了此频率处的频谱能量值较低，同时也是能量变化显著的区域。

图2 研究区频谱分析

进一步分析这3个频率所反映的图像特征，分别以20、37、47周/图像为中心频率，对原图像进行半径为5的频域理想滤波，并将滤波后的结果与原图像进行对比。取原图像T1的第89行以及滤波后图像的第89行比较，如图3所示，发现中心频率为37周/图像的滤波结果与道路绿地的轮廓位置符合最好。图2(d)中频谱剖面曲线的一阶导数曲线，如图4所示，频率37周/图像是曲线的最高峰值频率。

图3 不同中心频率的滤波结果

图4 图2(d)的一阶导数曲线

用同样的方法，截取原图像1(a)的另外3幅小区域256×256大小图像T2、T3、T4进行分析，如图5所示。发现3幅图像对应道路绿地覆盖轮廓的最佳截止频率分别为64、44、32周/图像。图5(b)、图5(d)、图5(f)分别为图像T2、T3、T4频谱最大能量方向一阶导数曲线，这几个频率均处于曲线的较高峰值位置。据此判断，做频谱图中最大能量方向一阶导数曲线，曲线的峰值位置可以作为判断图像轮廓截止频率的依据。37、64、44、32周/图像，位于频谱半径的中频0.32～0.64阶段，也即对应道路绿地覆盖轮廓的图像截止频率应主要分布在此频段范围。

对研究区原图像(图1(a))进行频谱特征识别，在45°、135°两个方向的频谱曲线如图6(a)、图6(b)所示，其一阶导数曲线如图6(c)、图6(d)所示。取中高频区域，也即频率半径的0.32～0.64频段的频谱能量峰值，那么45°方向的在频率122周/图像出现最高峰值；135°方向在频率为156周/图像出现最高峰值。因此，45°方向对应道路绿地覆盖轮廓的截止频率为122周/图像；135°方向截止频率为156周/图像。

图5 研究区截止频率建立

图6 研究区遥感图像频谱分析

3.3 滤波器设计

3.3.1 道路绿地覆盖轮廓特征滤波器

为了提取特定方向、特定频率的图像信息，需要设计适合图像特征的滤波器。Gabor滤波器与哺乳动物视网膜神经细胞的接收场模型相吻合，它能够同时在时域和频域获得最佳的局部化[9-10]，因此利用Gabor滤波器提取道路绿地的轮廓信息。一般情况下，二维Gabor函数h(x,y)可表示为：

h(x,y)=g(x′,y′)exp[2πj(Ux,Vy)]

(4)

其中，(U，V)表示特定的空间中心频率，g(x，y)为高斯函数，(x′，y′)是(x，y)旋转θ度角，即：

(5)

二维Gabor函数的傅里叶变换为：

(6)

其中，(U′，V′)和(u′，v′)表示类似于式(5)的角度旋转。

Gabor函数实际上是复数形式，将其分解可表示为：

h(x,y)=R(x,y)+jI(x,y)R(x,y)=g(x,y)cos[ϖ(xcosθ+ysinθ)]I(x,y)=g(x,y)sin[ϖ(xcosθ+ysinθ)]

(7)

如式(7)所示，Gabor函数可以分成奇函数和偶函数两个部分，实部是个偶函数，用实部对图像滤波仅能起到图像平滑，为了减少运算量，实现边缘的快速提取，仅使用Gabor函数的奇部进行道路绿地覆盖轮廓的图像滤波。一个θ方向的二维Gabor奇函数可具体表示为：

(8)

其中，ω为角频率矢量，ω=2πf，f为频率，σ为空间常数。

根据频谱分析结果，45°方向边缘截止频率为122周/图像，为了计算方便，将其除以频率半径，取其归一化频率，即122/279=0.44周/像元，那么奇Gabor滤波器角频率取值为ω=2πf=2π×0.44=2.76。135°方向的边缘截止频率为156/279=0.56周/像元，则奇Gabor滤波器角频率取值为ω=2πf=2π×0.56=3.52。

除了方向和角频率，空间常数也是奇Gabor滤波器设置的一个重要参数。

假设理想的阶跃边缘模型为：

U(x,y)=1(y≤0)U(x,y)=0(y>0)

(9)

做奇Gabor函数与理想边缘的卷积，并固定角频率ω和方向θ，则卷积函数随着空间常数σ变化曲线如图7所示。45°方向卷积函数随着空间常数的变化曲线见图7(a)，当σ=0.70，函数取最大值，因此45°方向滤波器参数设置为σ=0.70，ω=2.76，θ=π/4；135°方向卷积函数随着空间常数的变化曲线见图7(b)，当σ=0.55时函数取得最大值，因此135°方向滤波器参数为：σ=0.55，ω=3.52，θ=3π/4。

图7 两个方向的空间常数取值

3.3.2 道路特征滤波器

根据研究区近红外波段图像道路的基本特征，道路信息提取将主要提取色调信息，然后提取其脊线作为道路的中心线。据参考文献[7]，色调信息集中在频谱图的低频直流中心，因此道路色调信息通过对低频信息进行滤波得出。结合图1(a)，道路信息的直线方向特征显著，但色调最暗。因此，本文对原图像进行反变换，突出道路色调信息，然后进行低频滤波提取道路色调信息。

用式(6)Gabor滤波器进行频域滤波，滤波器中心设置为(0，0)，方向设置分别为45°和135°。在图1(c)中看到，随着频率增加，频谱能量迅速下降，当频率为20周/图像时，下降速度显著减慢，因此沿着长轴方向的σu取值为20周/图像。为了突出线性信息，令垂直长轴方向的σv为6周/图像。

4 研究结果与评价

对图像进行基于3.3.1参数设置的频域滤波，结果如图8所示。从滤波结果可以看出，道路绿地覆盖方向上的轮廓信息提取较为完整，闭合性好。同时部分非道路绿地、人工建筑等其他地物频率与确立的截止频率接近，也有部分提取出来。

图8 Gabor滤波结果

为了限制道路绿地的位置信息，利用3.3.2关于道路特征的滤波器，进行道路中心线辅助信息的提取。滤波结果如图9所示，图9(c)中道路的线性特征显著，通过二值分割和脊线提取，道路中心线如图9(d)所示。文中方法提取的道路中心线连续，没有出现断点，是直线型道路中心线提取较为实用的方法。

根据实地调查，研究区道路45°方向宽20m，约32像素，135°方向宽27m，约45像素。分别在道路中心线两侧设置缓冲区，对边缘提取图像进行缓冲区限制之后的结果如图10所示。图10(a)和图10(b)为两个方向的道路绿地覆盖轮廓，图10(c)为两个方向之和，图10(d)为二值化后的结果。从图中看出，建筑、道路等其他地物完全去除，道路绿地很好地检测出来，45°方向轮廓完整清晰，整体闭合性较好；135°方向轮廓清晰，较为完整，受图像道路绿地覆盖特征影响，右下角有部分边缘没有闭合，如图中圆圈所示。

图10 研究区道路绿地覆盖轮廓

将二值化道路绿地轮廓图像经过数学形态学开运算，去除斑点噪声，结合图像NDVI，最后结果如图11(a)所示。道路绿地覆盖轮廓曲线闭合，完整，内部基本无噪声影响。非道路绿地部分已被较好地去除。

为了进行比较分析，做图像的Canny边缘提取(图11(b))。由于Canny算法不是针对图像特征提取的边缘检测算法，没有顾及道路绿地覆盖轮廓的线性特征，因此检测结果中，部分受建筑物影响的道路绿地覆盖轮廓不连续，没有较好地提取出来，如图11(b)中圆圈所示部分。而且Canny算子检测的是图像的边缘点，这不仅包含道路绿地覆盖轮廓，还包含许多具有边缘特征的内部噪声。

研究区NDVI结果的二值化图像如图11(c)所示。NDVI主要依据绿地的光谱特征，能将大多数与绿地光谱差异大的地物去除。但是没有道路中心线辅助信息，利用NDVI提取的绿地结果中，存在大片的非道路绿地。尽管NDVI在中低分辨率绿地覆盖信息提取中具有优势，但对于高分辨率遥感图像说来，NDVI结果中的绿地轮廓过于光滑，不能充分显示与图像结构相关的细节信息。

图11 3种结果的比较

5 结束语

本文通过对高分辨率遥感图像城市道路绿地覆盖进行频谱分析，找到对应道路绿地覆盖轮廓的截止频率，设计奇Gabor滤波器提取两个方向上的线性轮廓信息。奇Gabor滤波器提取出的道路绿地轮廓线性特征显著，但在提取过程中受到了噪声影响，携带了部分其他地物方向特征信息。通过道路中心线进行道路绿地覆盖范围的限制，最终得到了道路绿地覆盖轮廓图。本文的研究方法还可以用到其他受线性条件约束的地物特征提取中。

本文的研究区位于城市的主、次干道，道路的类型为直线型。其他类型的道路如弯曲形，在进行滤波时需要调整滤波参数，因此关于其他类型的城市道路绿地覆盖还有待进一步研究。

参考文献：

[1] 宫鹏，黎夏，徐冰.高分辨率影像解译理论与应用方法中的一些研究问题[J].遥感学报，2006，10(1)：1-5.

[2] SMALL C.High spatial resolution spectral mixture analysis of urban reflectance [J].Remote Sensing of Environment，2003，88(1)：170-186.

[3] TOOKE T R，COOPS N C，GOODWIN N R.Extracting urban vegetation characteristics using spectral mixture analysis and decision tree classifications[J].Remote Sensing of Environment，2009，113(2)：398-407.

[4] 申广荣，钱振华，徐敬敬，等.基于eCognition的城镇绿地信息动态监测研究[J].上海交通大学学报，2009，27(1)：1-6.

[5] 童庆禧，薛永祺.空间遥感技术发展与应用的深化[C].中国空间科学学会第七次学术年会，2009.

[6] 叶泽田.频谱段图像及其应用的探讨[J].环境遥感，1993，8(2)：139-146.

[7] 王珂，肖鹏峰，冯学智，等.基于频域滤波的高分辨率遥感图像城市河道信息提取.遥感学报，2013，17(2)：269-285.

[8] 吴桂平，肖鹏峰，冯学智，等.基于光谱空间变换的遥感图像目标探测方法研究.光谱学与光谱分析，2013，33(3)：741-745.

[9] MEHROTRA R，NAMUDURI K R，RANGANATHAN N.Gabor filter-based edge detection [J].Pattern Recognition，1992，25(12)：1479-1494.

[10] 傅一平，李志能，袁丁.基于优化设计Gabor滤波器的边缘提取方法[J].计算机辅助设计与图形学学报，2004，16(4)：481-486.