邓鸿儒,崔宸洋,单文龙,徐梦竹
(1. 河海大学 地球科学与工程学院,江苏 南京 211100;2. 南京林业大学 土木工程学院,江苏 南京 210037)
建筑区是人类进行土地规划、城区扩张、城市监测、灾后评估等研究分析的必要信息,为提高研究效率并降低成本,利用遥感影像自动或半自动提取建筑区已经成为城市研究的重要技术手段。合成孔径雷达(Synthetic Aperture Radar,SAR)由于其全天时、全天候、高分辨率等特点,能够有效监测提取建筑区。近年来,随着SAR影像分辨率不断提高,数据资源不断丰富,使其成为城市研究的重要数据源之一。目前,国内外学者利用SAR影像开展了有关城市区域提取的很多研究。如Ban等人在Radarsat影像上对城市目标进行了有效分类[1],吴樊等从分辨率为3 m的机载SAR影像中利用灰度特征提取出近80%的居民区[2],赵凌君等利用分辨率为1 m的机载SAR影像采用变差函数的方法得到精度为78%的建筑区[3],张舞燕等通过相干系数利用ENVISAT ASAR的两景影像,提取出城市建筑区边界与光学影像吻合度达到817%[4],孙萍[5]等通过NASA/JPL AIRSAR的L波段数据提取出82%的建筑物。
但对于我国最新发射升空的高分三号卫星,由于刚发射不久,国内利用其进行相关应用研究还很少,在建筑区提取方面研究更少,仅有少量关于海浪定量遥感[6]、土壤水分反演的研究。盐城主要以平原为主,地区平均海拔不足5 m,地表起伏比较平缓,其城区比较适合进行建筑区提取的研究。因此,本文以盐城地区为例,开展基于高分三号SAR影像的城市建筑区提取研究。首先对建筑区的散射机制和其SAR影像特征进行分析,然后基于灰度共生矩阵提取纹理特征,针对高分三号SAR影像,确定灰度共生矩阵参数,并筛选最优纹理特征组合,最后综合利用灰度和纹理特征,从高分三号SAR影像中有效提取出建筑区范围,为高分三号卫星数据在城市研究中的应用提供参考。
表1 不同成像模式参数[7]Tab.1 Parameters of different imaging modes
高分三号卫星是我国在2016年发射的C波段多极化合成孔径雷达成像卫星,其分辨率可达1 m,能够全天候不间断获取全球海洋和陆表信息,能长期提供高质量高精度观测数据用于监测国土资源、维护海洋权益以及灾害天气预警等,对国家经济建设和国际区域合作都有促进作用。
高分三号卫星是目前世界上成像模式最多的合成孔径雷达卫星[7],包含了聚束、条带、扫描等12种成像模式,空间分辨率在1 m到500 m之间,幅宽可从10 km调整到650 km。因此,其应用范围较广,既可用于海面船只和陆地建筑物等精细目标的分辨,也可用于资源环境及生态状况的大范围普查。
高分三号对地观测所有成像模式见表1。
与可见光或红外遥感相比,微波遥感具有一定的穿透性,能够获取更多的目标信息。通常来说,可见光和近红外遥感通过地物反射率获取地物信息,热红外遥感通过地物热辐射获取地物信息,而合成孔径雷达中的主动微波遥感是通过后向散射回波获取地物信息[8]。由于自然地表几何特征复杂且分布随机,因此雷达入射微波与地表之间的散射作用比较复杂。建筑区包括建筑物、道路、公共设施、绿地等,SAR影像上建筑区的表现是各种地物综合反映的结果[9]。建筑区不同地物由于其表面粗糙程度不同,在微波的照射下会发生不同的散射,产生的回波强度也会有差异。根据建筑物表面粗糙程度以及立体结构,建筑区可能发生的散射情况有以下5种[10]:①由建筑物屋顶引起的镜面反射,如图1中a所示;②由光滑地面与建筑物墙面引起的二次散射,如图1中b所示;③在相邻建筑物墙面及建筑物之间的地面发生的多次散射,如图1中c所示;④由粗糙地表引起的表面散射,如图1中d所示;⑤由粗糙地表与建筑物墙面引起的二次散射,如图1中e所示。
图1 建筑区散射机制示意图Fig.1 The scattering mechanism diagram in building area
一般来说,屋顶表面可以被认为是光滑的,会发生很强的镜面反射,其回波受到建筑屋顶和墙面朝向的影响。如图1所示,对于斜顶建筑物,垂直于雷达观测角的倾斜屋顶产生的镜面反射本身是强反射;而无论斜顶还是平顶建筑物,当其墙面朝向雷达入射波方向时,墙面和地面形成的各向同性二次散射会产生较强的回波;此外,在密集建筑区,建筑墙面还会与地面之间发生多次散射[11]。因此建筑区在SAR影像上呈现出的整体亮度一般要高于周边地物,如图2a所示。在建筑区内,建筑物会引起亮斑或亮线,粗糙植被会引起浅灰色斑块,再加上黑色的道路、阴影等,整体就会形成较明显的区域边界,同时建筑区内建筑往往排列较为整齐,因此建筑区在SAR影像中会呈现出明暗相间的丰富纹理,如图2b所示。
图2 建筑区SAR影像特征Fig.2 SAR image features in building areas
根据上一节的分析,若仅仅利用SAR影像的灰度特征来提取建筑区,难以取得较好的效果,可以综合利用灰度和纹理特征来提高信息提取精度。通过纹理特征进行影像分析的方法有很多,而本文采用通过研究灰度的空间相关特性来描述纹理的灰度共生矩阵(Gray Level Co-occurrence Matrix,GLCM)的方法,是因为该方法被认为是描述影像纹理特征最有效的方法之一。
灰度共生矩阵是用来统计影像灰度变化的数学度量,它从影像(x,y)灰度为i的像元出发,统计与其灰度为j、距离为d的像元 (x+Δx,y+Δy)出现的概率P(i,j,d,θ)。用数学表达则为:
式中,i、j分别表示灰度共生矩阵的行列号;x、y表示影像中的像元坐标;Nx、Ny表示影像的行列数;θ表示两像素连线向量的角度,一般有0o,45o,90o和135o4个角度;d表示两像元之间的距离,为较好地分析结果,d取为1。
在实际应用中,通常需要根据灰度共生矩阵计算出的一系列特征量来表达纹理特征。在本文中,主要采用均匀性、对比度、方差、均值、差异性、熵、能量、相关度8种特征,其计算方法可参见文献[2],本文不逐一列出。
2.1.1 灰度共生矩阵参数确定
灰度共生矩阵主要受影像的量化级、方向、步长以及窗口大小等因素的影响,需要首先根据影像纹理特点选择合适的参数[12]。以往的研究表明,过多的影像灰度级数会导致计算量大大增加,然而灰度级数对纹理特征计算影响不大,因此本文选择把影像灰度级数压缩为16级;为了考虑更加全面,我们取0°,45°,90°和135°共4个角度方向计算,最后计算4个方向的平均值作为特征值;步长取值一般较小,根据经验本文取为1。通常来说,窗口大小对灰度共生矩阵有较大的影响,窗口过大,纹理特征会变模糊,窗口太小,又不能很好地反映影像的纹理特征,而且对于不同影像,适用的窗口大小也不同。因此,本文重点探讨窗口大小的确定。
由于本文的目标是从SAR影像中提取出建筑区,对于城市区域而言,非建筑区主要有植被、水体等典型地物,因此截取多个包含典型地物的非建筑区样本以及建筑区样本开展实验。以4为间隔依次采用从5×5到49×49之间不同的窗口大小进行纹理特征提取实验,针对建筑区、非建筑区两类地物,分别计算不同窗口大小下的8种纹理特征统计量,并画出每一种纹理特征统计量随窗口大小变化的曲线,据此分析出这种统计量用于提取建筑区的最佳窗口尺寸。下面以能量和熵为例进行说明,如图3所示,带方框的曲线代表非建筑区,带三角形的曲线代表建筑区,在窗口尺寸接近41时,如图3黄色椭圆内,非建筑区和建筑区特征值相差较明显,大于41时,特征值又再次接近甚至交叉。
图3 不同窗口建筑区与非建筑区纹理特征比较Fig.3 Comparison of texture features of different windows in building area and non - building area
由于本实验是为了区分大块的建筑区和非建筑区,无需内部细节,所以由实验的特征随窗口大小变化的曲线直接获得,窗口取为41×41较合适。
2.1.2 纹理特征选择
不同的特征值反映了影像不同的纹理特性,参与分类的特征过多,不利于得到好的分类效果,所以需要进行特征选择,根据前人研究,利用巴氏距离可以较好地对纹理特征进行选择[13]。巴氏距离定义为:
式中,μ1,σ1和μ2,σ2分别表示两个不同类别在某一特征影像上的均值和标准差。BD值越大表明该特征区分类别的能力越强。将每种特征的BD值除以最大值得到归一化后的BD值见表2。根据得到的BD值大小,可知均匀性、相关性、惯性矩区分类别能力最强。
表2 归一化BD值Tab.2 Normalized the BD value
根据上一节的分析,本文利用惯性矩、能量、相关性、差异性来描述纹理特征,并与灰度特征一起组合成一幅新影像,从中提取建筑区。该方法主要实现过程如下:
1)先用窗口大小为3×3的LEE滤波方法对影像进行斑点噪声抑制处理;
2)对滤波后的影像进行基于灰度共生矩阵的纹理分析,采用上一节确定的参数计算得到8个常用的纹理统计特征;
3)按巴氏距离最大原则挑选出3个最优的纹理特征,采用主成分分析方法去除相关性,并将信息量最大的前两个主成分与灰度影像进行波段组合;
4)利用K均值聚类分析的方法对波段组合后的影像分类;
5)最后对分类后的影像进行形态学滤波,去除单个孤立的面积小的点,得到建筑区。
研究区位于江苏省盐城市,地形主要以平原为主,本次实验数据是选取2017年4月20日精细条带1模式下的高分三号影像,极化方式为HH,标称分辨率为5 m,选取一幅大小500×500和一幅大小为700×700的子影像作为测试影像,如图4所示,同时选择同一地区对应的Landsat8影像作为参考影像,如图5所示。
图4 研究区高分三号影像Fig.4 GF-3 Satellite SAR image of study area
图5 研究区Landsat8影像Fig.5 Landsat8 image of study area
对滤波后的整幅影像进行纹理分析,得到均匀性、相关性、惯性矩三幅特征影像,再使用主成分分析法去除它们之间的相关性,选取前两幅信息量最大的波段,再与原灰度影像进行波段组合。利用K均值聚类分析方法,对组合后的影像进行分类,得到分类影像,通过一定的形态学处理,将孤立的点去掉,填充面积较小的孔洞,最后得到分类结果[14]。如图6是结合灰度纹理的非监督分类及后处理结果,a和b分别是实验区Ⅰ分类结果和与光学影像叠加后的结果,c、d是实验区Ⅱ分类结果和与光学影像叠加后的结果;图7则是对应的无纹理的分类结果。
从图5和图6可以看出,分类得到的影像与对应的光学影像比较,无论是区域Ⅰ,还是区域Ⅱ,对于建筑比较密集或者SAR影像纹理比较明显的区域,利用纹理和灰度特征能够更加容易地提取出建筑区。从图7看出没有纹理信息的提取结果相对较差。
图6 结合灰度纹理的非监督分类及后处理结果Fig.6 Unsupervised classif i cation and post-processing results combined with grayscale texture
图7 无纹理的非监督分类及后处理结果Fig.7 Untextured unsupervised classif i cation and post-processing results
根据研究区的地理范围,为了更精确的统计,均匀分布15行,15列,共225个检验样本点,检验样本点分布如图8、图9所示,以研究区的Landsat8 OLI影像为参考数据,对这些检验样本点逐一进行人工目视解译,确认检验点分类结果以及参考数据的类别属性,本文采用检测率,漏检率和错检率来进行结果评价,精度评价见表3。
图8 基于灰度纹理的精度评定Fig.8 Accuracy evaluation based on gray texture
图9 无纹理的精度评定Fig.9 Non-textured accuracy assessment
表3 总体分类精度统计Tab.3 Overall classif i cation accuracy statistics
本文针对高分三号影像,基于灰度共生矩阵提取纹理特征,再进行特征选择,通过主成分分析消除特征间相关性,与原灰度影像叠加后,采用非监督分类的方法,实现了高分三号研究区建筑区的自动提取。实验结果表明,对于如区域Ⅰ和区域Ⅱ小范围建筑区特征较为明显的影像,综合纹理灰度特征的分类方法精度较高,综合纹理灰度特征分类的检测率接近72%。而无纹理的分类方法检测率相对较低、漏检率和错检率较高。
从精度评定的结果看,这种方法仍然存在一些错分的情况,由于特征区分不够明显,可利用特征数偏少,将一些平原错分为建筑区的情况较为频繁,但总体提取效果较好,基本可以将大片的建筑区块提取出来。实验证明,这种方法除了可以对高分辨率SAR影像适用外,对高分三号卫星的影像同样适用。考虑到SAR影像中建筑区的干涉相干特性与其他地物相比有较明显的区别,后续研究将把相干性特征引入到SAR建筑区提取中,以期进一步提高SAR影像建筑区提取的精度。