王燕红，程博，尤淑撑，武盟盟

（1.中国科学院遥感与数字地球研究所，北京 100094；2.中国科学院大学，北京 100049；3.中国土地勘测规划院，北京 100035）

1 引 言

近些年，城市扩张和城市建设的迅猛发展，运用遥感技术进行城市环境研究具有重要应用价值。星载合成孔径雷达（Synthetic Aperture Radar，SAR）具有全天时、全天候、穿云透雾的优势，对我国幅员辽阔、气象气候和下垫面条件复杂、且急需快速获取监测数据的前提下，具有明显的优势。2007年以来，加拿大 RADARSAT－2、德国 TerraSAR－X及意大利COSMO－SkyMed等分辨率达1m雷达卫星的陆续成功发射，高分辨率城区SAR图像的获取能力大大提升，一方面为准确提取建筑区提供了前提，但是另一方面却导致高分辨率条件下城区环境更为复杂，给高分辨率SAR影像应用带来了新的挑战。

建筑区提取是城区遥感图像解译的重要研究课题之一。建筑区是城市环境最基本的一部分，一般由住宅区、公共建筑、道路与绿地和场地组成。在高分辨率SAR图像上，灰度分布比较复杂，用传统的基于图像灰度的提取方法很容易受到居民地内部要素以及噪声的影响而得不到满意的效果，针对高分辨率SAR影像中含有丰富的纹理信息特点，纹理分析方法是对像素间空间分布关系进行研究，因而可利用纹理分析的方法来解决。国内外学者也做了大量的研究，如Dell’Acqua等［1］利用灰度共生纹理度量来刻画城区SAR图像上不同区域的建筑覆盖密度，从而区分城市中心、居住区和郊区。Tison等［2］利用马尔科夫随机场纹理分类方法对城区进行精细的分类。强永刚等［3］利用小波变换和数学形态学提取人工建筑区。赵凌君等［4］利用变差函数的方法，实现了建筑区的快速提取。变差函数方法是一种较新的纹理提取方法，在建筑区和非建筑区上的特征表现差异大，算法简单、快速，具有重要应用价值，然而目前常用的变差函数模型易受噪声干扰，稳健性差［5］。本文提出一种基于中值滤波的变差函数纹理分析方法提取建筑区，通过改进变差函数的计算方法，使窗口的大小不再受步长的约束，不仅保留了变差函数区分建筑区与非建筑区的优势，而且改善了高分辨率SAR提取建筑区时受强反射点、噪声干扰大的缺点，提取效果更完整。

2 变差函数纹理分析

2.1 变差函数理论

变差函数理论由数学专家G.Maberon教授1962年创立，作为地质统计学的重要工具，变差函数被应用于空间随机场的统计特性研究。随着变差函数理论的不断发展，其应用领域从早期的空间数据分析、地理数据位置分析等逐渐扩展到遥感数据分析，尤其是在20世纪90年代，Miranda等人成功将变差函数纹理分析方法应用于SAR图像植被识别之后［6-7］，变差函数成为了一种重要的纹理特征分析方法，被广泛应用于遥感图像信息提取领域。

变差函数又称半变差函数（Semivariogram Function），定义为区域化变量Z（x）和Z（x＋h）（同时包含两点距离和方向信息）两点之差的方差之半，即

对于离散的栅格数据，实验中的变差函数定义为：

其中，N（h）表示观测数据中间距为h的点对数目，估计值γ＊（h）通常称为实验变差函数。变差函数用于度量区域化变量的空间相关性，能充分反映图像数据的随机性和结构性。图1是一个理想的变差函数模型，当h＝0时的变差函数值C0为块金值（Nugget Value），代表不确定的变化，随着h的增大，空间样本之间的相关性越来越小（非相关性越来越大），当样本变为相互独立时的距离a称为变程，相应的变差函数值逐渐增大并在间距为a时趋于平稳值C0＋C1，这一平稳值称为基台值（Sill）。

图1 理想的变差函数模型

2.2 改进的变差函数计算方法

变差函数用于纹理分析的常用计算方式是首先确定步长h、窗口w、计算方向，然后计算窗口w内所有间距为h的点对的半方差值，然后取平均（如公式（2））作为窗口中心点的变差函数值，遍历全图即得到影像的变差函数特征图。对于高分辨率的SAR图像而言，窗口w的选取取决于h，要保证窗口w内间距h的点对数目足够，窗口w至少应为3h～5h［4］。然而，窗口w过大会造成图像整体模糊、边缘虚警率高，且窗口内取平均的方法易受噪声、孤立强反射点干扰，算法稳健性差。

根据变差函数计算方法所存在的弊端，本文提出一种稳健的变差函数计算方法，该算法继承了变差函数和标准中值滤波方法的优点，窗口w取值不再受h约束，主要思想是：在计算像素点的变差函数值时，不再取窗口w内间距为h的点对来求取，而是计算像素点（x，y）与固定方向（如0°方向）上间距h的点（x＋h，y）的半方差值，将此值作为点（x，y）的半方差值；然后取窗口内所有像素的半方差值的中值代替均值作为变差函数值。

计算公式如下：

其中，w为二维模板。

下面以窗口w＝3，h＝1，全方向为例，如图2所示。

图2 变差函数纹理特征计算示意图

3 基于改进的变差函数纹理分析的建筑区提取方法

3.1 高分辨率SAR影像的建筑区纹理特征

遥感影像上可以分为建筑区与非建筑区两大类。建筑区可分为城镇建筑区和农村建筑区。在高分辨率SAR中，建筑物由于角反射器及屋顶强反射而产生的亮斑，通常灰度（亮度）要高于周边地区，然而城镇建筑区和农村建筑区由于建筑结构、材质、分布的差异，纹理特征也不尽相同。城镇建筑物分布整齐，楼房之间间隔大，大多数为平顶整齐的高层楼房，使用材料大多具有良好反射率，当屋顶在垂直于雷达观测角的情况下产生的镜面发射是强反射，因此多表现为亮斑。而高层楼房之间道路、草坪等散射回波较弱，在影像上表现为亮斑之间的暗斑块，加之建筑物通常排列较为整齐，形成有一定周期规律的明暗相间的纹理特征。农村建筑区房屋间隔小，分布相对密集，通常为房屋聚集的村庄，在影像上表现为棋盘式分布的亮斑，纹理特征不如城镇建筑区有规律。非建筑区主要包括植被、水体等，整体亮度较低。

3.2 建筑区的变差函数曲线

高分辨率SAR图像中建筑区、植被区和水域的实验变差函数曲线如图3所示（A为建筑区，B为植被区，C为水域）。从图中可以看出，首先，建筑区的变差函数曲线整体上都远远高于植被区和水域，且随着间距h的增加周期性地出现峰值和谷点，但是在数值上，峰值与峰值、谷点与谷点不一定严格相等，随着h的增大会表现出衰减的现象［8-9］，这是由于建筑区是富含大量强散射点，且建筑区与街道等间隔分布；植被区的变差函数曲线随着间距h的增大，变差函数曲线先缓慢上升，随后趋于平稳。水体由于镜面反射回波很弱，像素灰度值很低且非常均匀，变差函数曲线近似为一条直线。建筑区的变差函数曲线具有特殊性，与非建筑区特征有明显差异，因此可以作为有效区分建筑区和非建筑区的依据。

图3 实验区各地物的变差函数曲线

3.3 变差函数参数选取

（1）步长h

如图3（b）所示，变差函数曲线在h1＝6的位置达到第一个峰值，h2＝10的位置达到波谷。在0～h1区间内，曲线上升非常快，不具有稳定性，因此步长h应大于h1；另外考虑到随着步长增大，曲线逐渐衰减，相关性减小，因此本文选取稍大于h1，小于h2的位置作为步长。这样既具有稳健性，又能保证相关性比较强。

（2）窗口w

一般来说，窗口选择太大可能包含多余的无用信息；窗口太小又不能有效、全面地描述地物的特征，而且当窗口选择的较小时，噪声对纹理的计算影响较大。所以，合适的窗口大小对于纹理的计算十分重要。在实际应用中，一般根据影像地物的实际特征来确定窗口的大小。

（3）方向

要考虑的因素是变差函数的各向同性和各向异性，即变差函数计算方向的选择。一般采用全方向。

3.4 提取流程

建筑区提取流程如图4所示，首先进行变差函数实验，确定变差函数参数h、w，然后计算目标影像的变差函数纹理图，利用EM算法对得到的纹理图进行阈值分割，得到初步的提取结果，利用形态学方法进行后处理，包括填洞，去除小面积区域等，最后得到建筑区提取结果。

图4 算法提取流程

4 实 验

4.1 实验数据

为验证算法的有效性，选取了高分辨率图像中农村建筑区（实验区1）和城镇建筑区（实验区2）进行了实验。实验区1和实验区2均采用2007年8月19日获取北京市升轨右视的TerraSAR－X数据，极化方式为HH，空间分辨率为1.25m，大小均为1600×1200像素。

4.2 变差函数参数的确定

根据图3的变差曲线特征（本文所有实验影像均来自同一景影像），本文设定步长h＝8。一般来说，窗口w的选定根据实验区特征选取，目前多数依靠经验选取。为了定量的分析窗口大小对提取结果的影响，本文通过固定步长h，改变窗口w的大小，计算实验区的检测率、虚警率、漏警率，以此选择最优的窗口w。

图5为固定h＝8，检测率、虚警率、漏警率随窗口变化的曲线。两实验区的曲线变化大致相同，检测率曲线先上升后下降，漏警率呈下降趋势，虚警率则是先下降后上升。当窗口在0～20时，检测率一直上升，大约达到23～28的时候，检测率最高，此时虚警率和漏警率较小。因此本文选取w＝23作为最优窗口。

4.3 实验结果及评价

图5 检测率、虚警率、漏警率随窗口变化曲线图

根据上节提出的方法对原始影像进行建筑区提取（h＝8，窗口w＝23，全方向），本文将结合光学影像人工提取的建筑区作为标准结果。实验区的初始提取结果如图6（b）、图7（b）所示，然后将初始结果进行了填洞、去除小区域（小于300像素），与标准影像进行叠加如图6（c）、图7（c）所示，并且与未改进变差函数的方法的后处理结果进行对比（白色部分为正确区域，蓝色部分为漏警区，黄色部分为虚警区）。

从对比图可以明显看出，本文提出的方法更能准确、完整地提取出建筑区的边界，边缘虚警区小，建筑区内部的空洞很少，一定程度改善了受孤立强反射点、噪声干扰影响的缺点，提取效果优于未改进的方法。为了使得对比结果清晰，本文在实验区1截取建筑区边缘部分进行了放大，如图8所示，首先观察建筑区边界，本文的方法边界的虚警很小，而未改进方法明显地扩大的边界范围，其次对于影像中的道路旁边的强散射点，未改进方法错误地提取为建筑区，造成了虚警，而本文方法一定程度改善了这一缺陷。

为了定量地评价提取结果，本文采用检测率（Detection Rate），虚警率（False Alarm Rate）［10］和漏警率（Missing Alarm Rate）3个度量值进行结果评价，评价结果如表1所示。从结果可以看出，改进的方法检测率较高，虚警率和漏警率较为改进方法均有明显降低，总体明显优于未改进方法提取的结果。

图6 实验区1提取结果

图7 实验区2提取结果

图8 提取结果细节图

表1 实验区精度评价表

5 结束语

本文提出了一种改进变差函数的建筑区提取方法，该方法以建筑区及非建筑区的变差函数特征为依据，通过建筑区变差函数实验提取计算间距h，进而确定窗口w，计算此参数下的变差函数值作为分类特征，采用EM算法实现建筑区和非建筑区的两分类问题。相比于传统的变差函数计算方法，本文在计算方法上提出改进：计算像素点（x，y）与固定方向（如0°方向）上间距h的点（x＋h，y）的半方差值，将此值作为点（x，y）的半方差值，然后计算窗口w内的所有像素点的半方差值的中值作为变差函数值。结果表明，改进的方法能够更完整、更高效地提取建筑区。但是对于纹理特征不明显，灰度值较低的建筑区，提取的效果不是很理想，下一步将尝试引入其他特征对纹理特征不明显的建筑区进行提取算法研究。

［1]DELL’ACQUA F，GAMBA P.Texture－based characterization of urban environments on satellite SAR images［J］.IEEE Transactions on Geoscience and Remote Sensing，2003，41（1）：153－159.

［2]TISON C，NICOLAS J M，TUPIN F，et al.A new statistical model for Markovian classification of urban areas in highresolution SAR images［J］.IEEE Transactions on Geosicence and Remote Sensing，2004，42（10）：2046－2057.

［3]强永刚，殷建平，祝恩，等.基于小波变换和数学形态学的遥感图像人工建筑区提取［J］.中国图象图形学报，2008，13（8）：1459－1464.

［4]赵凌君，高贵，匡纲要.基于变差函数纹理特征的高分辨率SAR图像建筑区提取［J］.信号处理，2009，（25）9：1433－1442.

［5]金毅，潘懋，姚凌青，等.一种稳健变差函数计算方法［J］.北京大学学报（自然科学版），2009（6）：1033－1038.

［6]MIRANDA F P，FONSECA L E N，CARR J R，et al.Analysis of JERS－1（Fuyo－1）SAR data for vegetation discrimination in northwestern Brazil using the semivariogram textural classifier［J］.International Journal of Remote Sensing，1996，17（17）：3523－3529.

［7]MIRANDA F P，FONSECA L E N，CARR J R.Semivariogram textural classification ofJERS－1（Fuyo－1）SAR data obtained over a flooded area of the Amazonrainforest［J］.International Journal of Remote Sensing，1998，19（3）：549－556.

［8]CURRAN P J.The semivariogram in remote sensing：An introduction［J］.Remote Sensing of Environment，1988（24）：493－507.

［9]OLIVER M，WEBSTER R，GERRARD J.Geostatistics in physical geography.Part I：Theory［J］.Transactions of the Institute of British Geographer，1989，41（3）：259－269.

［10]SHUFELT A.Performance evaluation and analysis of monocular building extraction from aerial imagery ［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1999，21（4）：311－326.