利用边缘密度特征提取高分辨率遥感影像中的居民区

陈 洪， 陶 超， 邹峥嵘， 邵 磊

中南大学地球科学与信息物理学院，长沙410083

随着遥感影像分辨率的逐步提高，居民区提取已成为一个热门课题.准确而快速地提取居民区可以为土地管理、城市规划等市政部门在做土地利用现状调查和宏观规划等方面的工作时提供重要的决策支持.在遥感影像中，居民区不仅具有覆盖面广、地物信息丰富的特点，而且它是一个随着时间动态变化的区域.人工提取居民区虽然可以保持较高的准确度，但是仅仅依靠人工的监测和分割提取居民区不但耗时而且成本高.因此，自动居民区提取方法越来越受到关注.

目前，许多研究学者提出了关于居民区提取的方法.文献[1]首先通过数学形态学方法从卫星影像上提取结构信息，然后利用神经网络分类居民区，该方法需要人工训练样本集的支持.文献[2]融合傅里叶变换和神经网络分类算法提出了一种居民区提取方法，该方法同样需要构建训练样本.文献[3]采用灰度共现矩阵方法从高分辨率SAR影像上提取居民区.文献[4]利用多级分类器在Ikonos灰度影像上提取居民区.文献[5]结合结构特征与图形学理论从卫星影像提取居民区，该方法使用了统计分类器.文献[6]利用建筑物密度分类居民区，该方法是基于影像纹理信息而提出的.文献[7]以纹理单元为模型检测高分辨率遥感影像上的居民区，该方法利用支持向量机进行分类.文献[8]融合了多种不同的分类器从全色卫星影像上提取居民区.文献[9]利用局部边缘分布信息检测高分辨率遥感影像上的居民区.分析现有居民区提取方法，发现存在以下两方面的局限性：1)大多数方法基于监督分类机制，需要大量训练样本保证分类精度，干扰因素较多，自动化程度有限；2)现有的居民区提取方法强调影像整体的纹理、光谱等特征，却忽略了居民区的局部特征，如边缘特征.在影像中，居民区内部包含大量的人工地物(如建筑物、道路等)，相对于其他非居民区(如耕地、林地等)具有更为明显和丰富的边缘特征，因此可以根据边缘密度差异提取居民区.

本文提出一种无监督的高分辨率遥感影像居民区自动提取方法，首先将居民区内的边缘特征直线化表达，然后以高斯函数为模型将边缘密度特征量化，最后通过阈值分割提取居民区.

1 居民区边缘密度特征提取

居民区与非居民区不同，其内部包含了密集的建筑物和道路等地物，这些地物都具有很明显的边缘特征；非居民区(如耕地、林地等)内的纹理较单一，没有明显的边缘特征.因此，居民区具有更高的边缘密度，根据边缘密度特征差异就可以区分识别居民区与非居民区.本文按照以下4个步骤提取影像上的边缘密度特征：

步骤1 采用Mean Shift算法平滑原始影像，并检测平滑影像上的边缘；

步骤2 将所有边缘拟合成直线段；

步骤3 计算每个像素点到所有边缘线段的空间距离，并用集合存储；

步骤4 根据距离集合，利用高斯函数构建空间投票矩阵，最终获取影像上的边缘密度特征.

在复杂的高分辨率遥感影像中具有大量的无关噪声，因此在提取边缘之前需要做平滑去噪处理.Mean Shift滤波与常用的去噪算法(如中值滤波、高斯滤波等)相比具有良好的保边缘性[10-11]，于是本文采用Mean Shift滤波平滑原始影像，平滑结果如图1中的(b)所示.由图1可知，Mean Shift具有良好的保边缘性，因此在滤波后的影像中，大量纹理噪声被有效地抑制，同时居民区内的边缘细节也保留得较好.在有效地滤除掉影像中的噪声后，采用Canny算子提取影像中的边缘.图1中的(c)和(d)分别显示了在原始影像和Mean Shift滤波影像中的边缘提取结果，从中可以看出：直接对原始影像做边缘提取，无论是居民区还是非居民区，其边缘都非常密集，边缘密度没有区分度；而在Mean Shift滤波后影像中，非居民区内的边缘比较稀疏，但居民区仍然拥有较高的边缘密度，因此可以利用影像中边缘密度分布推断影像中居民区的位置.

图1 平滑及边缘影像Figur e 1 Smoothing and edge images

假设经过Canny边缘提取得到一组边缘特征集合E0={e1,e2,···,en}，影像中边缘密度特征提取分为以下3个步骤：

步骤1 如图2所示，将检测到的边缘特征拟合成一条或多条直线段，依次按照先左至右，后从上至下的顺序，记录下每条边缘上的所有像素点，并令P0,P1,···,Pn为某一边缘上的像素点.根据以下步骤进行直线拟合：

步骤1-1 令边缘的起始像素点P0为固定点，P2为浮动点，两点构成直线段P0P2，计算顶点P1到P0P2的垂直距离a；

步骤1-2 设定最大容许距离值为ε，若a≤ε，则P3定义为新的浮动点，计算P1、P2到P0P3垂直距离的最大值，假设为b，并将其与ε再次进行比较；若a＞ε，则用直线段P0P2替代边缘折线P0P1P2，然后将P2定义为新的固定点；

步骤1-3 重复执行上一步直至遍历所有边缘像素点，最终获取边缘的直线表示形式并记录下每条直线段两个端点的像素坐标.直线拟合后的边缘集合表示为，其中表示第i条边缘上的第j条直线段，分别将它的两个端点坐标记为和

图2 直线拟合算法Figure 2 Line straight algorithm

如图3所示，将所有边缘拟合成直线段后，每条边缘均用随机的颜色显示.边缘特征经过直线拟合后，边缘像素点总数减少了，边缘直线表达也得到了简化，因此出现了局部边缘信息缺失的现象，但总体而言，居民区内的边缘密度仍然保持较高的水平，满足居民区提取的要求.

图3 边缘直线拟合影像Figure 3 Edge straight line image

步骤2 计算图像上每个像素点(假设共有R个像素点)到所有边缘直线段(假设共有K条边缘线段)的距离，用集合D={d11,d12,···,d1K,d21,d22,···,d2K,···,dR1,dR2,···,dRK}表示，其中dij表示第i个像素点到第j条边缘的距离，计算公式如下：

式中，(xi,yi)为第i个像素点的坐标，和分别为第j条边缘直线段的两个端点的像素坐标.

步骤3 利用步骤2计算得到的距离集合衡量边缘密度特征.假设原始影像大小为M×N，定义投票矩阵V(M,N).如果像素点p(xi,yi)是居民区内的像素点，由于居民区内的边缘更为密集，因此从统计学的角度来说，该像素点到所有边缘直线段的空间距离的总和会相对较小，然后通过式(2)将空间距离转化为投票数V(xi,yi)，则V(xi,yi)的值就更大；而非居民区内的像素点到所有边缘线段的空间距离的总和相对较大，故投票数就更少

式中，K为边缘总条数，dij为根据式(1)计算得到的像素点到边缘直线段的距离，σ为邻近投票参数，3.2节对该参数进行了具体分析.

图4为根据上述方法计算得到的边缘密度特征图，图中红色区域为拥有投票数较多的区域，即为边缘密度较高的地方；蓝色区域获得的投票数最少，即为边缘密度较低的地方.由图4可知，居民区内具有更高的投票数，离居民区中心越远，投票数就越少.

图4 边缘密度特征图Figure 4 Edge density features image

2 基于边缘密度特征居民区提取

如图4所示，边缘密度特征图中投票值较大的区域就可能为居民区，而非居民内投票值普遍偏低，所以影像上所有像素点的投票值服从一个双峰分布，则在两个峰值之间必然存在一个最佳的分割阈值，本文参考Ostu方法[12]自动获取该阈值：

1)对所有投票值进行排序，用集合Vs={v1,v2,···,vL}表示，令投票值为vi的像素点的数目为ni，则总的像素数位N=n1+n2+···+nL；

2)将集合Vs拆分两个类V0={v1,v2,···,vk}和V1={vk+1,vk+2,···,vL}，然后分别计算两个类的类内方差、类间方差和总方差；

3)进行循环运算，直至式(3)取得最大值时判定取得最佳阈值k∗

通过上述方法获取边缘密度特征图的最佳分割阈值，然后进行图像分割，如图5中的(a)为阈值分割后的二值图，白色区域即为居民区，5中的(b)用红色曲线标记居民区外围边界，可见提取结果与影像上的居民区具有很高的吻合度，说明本文方法具有可行性.另外，对于不同的实验影像，利用Ostu方法获取的分割阈值各不相同，因此本文方法是一个自适应的居民区分割提取方法.

图5 居民区二值及边界线图Figur e 5 Built-up area binary and boundary line images

3 实验结果与分析

3.1 居民区提取结果分析及评价

整个算法流程已由前面的例子完整演示.为验证本文算法的准确性及鲁棒性，下面采用多幅高分辨率遥感影像进行实验.另外通过人工方式准确地标记出居民区.与此同时，采用Paul Rosin的基于边缘特征的显著区域检测方法[13]进行对比试验，实验结果如图6所示，黄色区域人工标记居民区结果，红色区域为本文方法实验结果，蓝色区域为Paul Rosin方法实验结果.用于实验的10幅影像分辨率在0.5～4.0 m之间，场景均较复杂，包含各种具有密集边缘特征的地形地貌情况，如规则田地、山地起伏地区，以及林地、耕地、河流交织区域，居民区提取结果良好，说明本文方法具有良好的实际性能.

为了度量居民区提取精度，参考文献[14]的评价方法，首先分别统计以上10幅影像经本文方法及Paul Rosin方法所提取的居民区像素总数Sauto，人工提取的理想居民区像素总数Smanual，两者相同的像素总数Scommon，然后分别以Scommon与Smanual之比为准确率Pd，以两种方法错检像素总数(Sauto-Scommon)与Sauto之比为虚警率因子Pf，评价结果如图7所示：

图7中的(a)为本文方法和Paul Rosin方法提取以上10幅影像的居民区的精度对比图，本文方法提取精度稳定保持在一个较高的水平，而Paul Rosin方法则波动性较大；图7中的(b)为对应的虚警率因子图，Paul Rosin方法虚警率普遍高于本文方法，说明本文方法具有较高的准确度性.

3.2 实验参数分析

图7 提取结果精度及虚警率图Figur e 7 Images of accuracy and false alarm rate

本文涉及的可调参数可用于边缘特征量化的邻近投票参数σ、Mean Shift平滑参数(包括空间域带宽Bs和频率域带宽Br)以及Canny边缘检测参数.其中Mean Shift平滑主要应用在影像预处理阶段，因此相关的平滑参数取值选择在合理的取值范围内，提取结果精度能稳定保持在较高的水平.另外，本文方法不要求获取准确的边缘特征，且要求居民区内边缘密度明显比非居民区高，因此Canny边缘检测的参数选择对本文提取结果影响不大.综上所述，对本文算法提取结果影响最大的参数为用于量化边缘特征密度的邻近投票参数σ.

以实验的10幅影像为例，分辨率为0.5～4.0 m，逐一分析上述各参数对居民区检测成功率的影响：1)图8(a)给出了邻近投票参数σ取值为10～50时，10幅实验影像的平均居民区检测成功率的变化情况，σ分别取值34时，平均居民区检测成功率取得最大值，σ取值为26～34时，居民区检测成功率变化不明显，稳定保持在88%以上.σ取值不能过大或者过小，因为不同分辨率的影像均有其相适应的最佳σ取值范围；2)图8(b)为空间域带宽Bs和频率域带宽Br与10幅实验影像平均检测民区成功率的关系，其中红色折线表示Br=20,Bs取值10～30之间时，平均检测成功率的变化情况，蓝色折线表示Bs=18,Br取值10～30之间时平均检测成功率的变化情况，分析可知，两者取值在16～22之间时，居民区检测成功率变化不大；3)图8(c)为Canny边缘检测的相关参数与10幅影像的平均居民区检测成功率之间的关系，边缘检测的阈值分割参数在MATLAB编程环境下自动获取，另外一个可变参数为高斯滤波标准方差α，默认取值为，图8(c)中α取值在0.5～3.0之间时，平均检测率的变化情况，整体波动较小，均能保持在88%以上，说明参数α对实验结果影响不大.综合上述分析，当以上参数选取在合理的取值区间内时，本文方法提取居民区的精度能够稳定保持在较高的水平，分析结果表明本文算法对于参数变化具有良好的鲁棒性.

4 结语

本文利用居民区边缘密度特征并结合空间投票机制，提出一种新的高分辨率遥感影像居民区自动提取方法.实验表明，该算法具有较好的准确性和鲁棒性，可以从场景复杂的高分辨率遥感影像上提取居民区，主要创新点是以边缘特征为基础，发展了一种新的直线拟合方法，并以高斯函数为模型构建了一种新的空间投票机制，从而准确有效地度量了直线拟合后的边缘密度特征.

对于地物信息过于复杂的影像，例如当居民区周围存在过多零散分布的灌木丛或独立树木时，或者当实验影像上居民区与非居民区比例相差较大时，本文算法的提取结果仍然存在一些错判问题.因此，下一步工作考虑在居民区提取过程中加入除边缘特征外更多的特征信息以提高检测的准确度，同时进一步提高算法的运行效率，以实现既快速又准确的高分辨率遥感影像居民区自动提取方法.

[1]BENEDIKTSSON J A,PESARESI M,ARNASON K.Classif ication and feature extraction for remote sensing images from built-up areas based on morphological transformations[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(9):1940-1949.

[2]WEIGL K,GIRAUDON G,BERTHOD M.Application of projection learning to the detection of builtup areas in spot satellite images[R].Rapports de recherche-INRIA,1993.

[3]FANG W,CHAO W,HONG Z.Residential area information extraction by combining China airborne SAR and optical images[C]//IEEE International Geoscience and Remote Sensing Symposium,2004.

[4]ZHONG P,WANG R.A multiple conditional random f ields ensemble model for built-up area detection in remote sensing optical images[J].IEEE Transactions on Geoscience and Remote Sensing,2007,45(12):3978-3988.

[5]UNSALAN C,BOYER K L.A theoretical and experimental investigation of graph theoretical measures for land development in satellite imagery[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2005,27(4):575-589.

[6]KARATHANASSI V,IOSSIFIDIS C,ROKOS D.A texture-based classif ication method for classifying built areas according to their density[J].International Journal of Remote Sensing,2000,21(9):1807-1823.

[7]BRUZZONE L,CARLINL.A multilevel context-based system for classif ication of very high spatial resolution images[J].IEEE Transactions on Geoscience and Remote Sensing,2006,44(9):2587-2600.

[8]FAUVEL M,CHANUSSOT J,BENEDIKTSSONJ A.Decision fusion for the classif ication of urban remote sensing images[J].IEEE Transactions on Geoscience and Remote Sensing,2006,44(10):2828-2838.

[9]HU Xiangyun,SHENJiajie,SHAN Jie,PAN Li.Local edge distributions for detection of salient structure textures and objects[J].IEEE Geoscience and Remote Sensing Letters.

[10]COMANICIU D.Mean Shift:a robust approach toward feature space analysis[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2002,24(5):1-18.

[11]COMANICIU D,RAMESH V,MEER P.Real-time tracking of non-rigid objects using Mean Shift[C]//IEEE Conference on Computer Vision and Pattern Recognition,2000.

[12]OTSU N.A threshold selection method from graylevel histograms[J].IEEE Transactions on Systems,Man and Cybernetics,1979,SMC-9(1):62-66.

[13]ROSIN P L.A simple method for detecting salient regions[J].Pattern recognition,2009,42(11):263-2371.

[14]SIRMACEK B,UNSALAN C.Built-up area detection using local feature pointsand spatial voting[J].IEEE Geoscience and Remote Sensing Letters,2010,7(1):146-150.