CN106940782A - 基于变差函数的高分sar新增建设用地提取软件 - Google Patents

Abstract

本软件目的是对高分辨率SAR影像进行新增建设用地的自动提取，属于图像识别的领域。软件基于变差函数的方法，能够自动准确提取高分SAR影像中的新增建设用地。该软件的主要步骤如下：第一、对输入的图像进行预处理(滤波)；第二、计算变差函数获得两幅影像的纹理特征；第三、计算两幅纹理特征影像的比值，构造差异图像；第四、对比值影像作阈值分割处理，生成一幅二值影像，即初步新增建设用地提取结果；第五、对二值影像作后处理，生成最终提取结果，并进行精度评价。本软件的提取精度经内部测试可达到80％以上，可应用于城市规划、拆迁补偿界定等建筑变化区域的检测等方面。

Description

基于变差函数的高分SAR新增建设用地提取软件

技术领域

本软件属于信息技术图像识别的领域，可用于提取不同时相的高分辨率合成孔径雷达(Synthetic Aperture Radar，SAR)影像中新增加的建设用地。

背景技术

SAR以其全天时、全天候侦察和高分辨率成像的优势成为对地观测的重要手段之一。近年来，随着大量米级、亚米级高分辨率SAR图像的获取，基于SAR图像的城区环境研究成为当前SAR图像解译领域的重要课题之一。图像识别作为其中的关键环节，为城区土地利用情况调查、变化检测、地图绘制、灾害监测等专题解译提供了基础。目前，SAR图像识别以基于图像灰度分布模型的统计类方法和基于纹理分析的方法为主。城区存在大量建筑物等强散射目标，使得很多经典模型难以很好地拟合图像数据，导致统计类方法性能下降。此外，统计类方法大多采用逐像素分类方式，忽略图像的空间分布特性，分类结果存在明显的“椒盐”现象，进一步造成分类精度难以达到实用要求。纹理分析考虑了邻近像素间的空间信息而不仅是像素灰度信息，从而成为城区SAR图像识别的重要方法，变差函数法是近年来在遥感领域兴起的一种纹理分析的有效工具，广泛用于多光谱图像、DEM等遥感数据分类，在SAR图像的植被识别、建筑区提取等研究中亦得到应用。《信号处理》中《基于变差函数纹理特征的高分辨率SAR图像建筑区提取》首次将变差函数应用于高分辨率SAR影像的建筑区提取，取得了不错的效果。

变差函数用于纹理分析的常用计算方式是确定间距h、窗口大小w、计算方向，通过计算窗口w内所有间距为h的点对的半方差值，取平均作为窗口中心点的变差函数值，遍历全图即得到影像的变差函数特征图。一方面，窗口内中心元素的变差函数值是由窗口内所有距离是h的点对的半方差值取平均而得，这种取平均的计算方法很容易受噪声、孤立强反射点干扰，算法稳健性差；另一方面，由于以往的方法中参数确定只依靠经验，参数如果选取不当，对结果影响大，稳定性不高。

此外，SAR影像斑点噪声严重影响提取结果的精度。目前大多SAR图像斑点抑制滤波研究是针对中低分辨率的图像进行处理。传统的基于局部统计的自适应滤波器在处理中低分辨率的SAR图像时或者在不要求高分辨率的应用目的下，可以达到很好的去噪效果。但是对高分辨率SAR影像结构特征的保持能力不够，无法满足实际应用需求。

作为高分辨率SAR影像的显著特征，其结构信息是SAR图像解译和信息提取的重要依据。 本发明将结构检测技术引入到传统滤波中，针对单极化高分辨率SAR影像，研发了基于结构检测的相干斑抑制滤波(SDBSF)，并根据变差函数计算方法所存在的弊端，提出一种稳健的改进型变差函数方法，该算法继承了变差函数和标准中值滤波方法的优点，同时也提出一种确定最优参数的方法。然后利用该方法进行建筑区纹理特征提取，自动提取不同实相新增建设用地。

发明内容

本发明的目的在于提出一种自动提取高分辨率SAR影像新增建设用地的方法，提高传统方法的检测精度。

为实现上述目的，本发明提出的完整方法为：

第一步、影像滤波

1-1)输入不同时相的高分辨率SAR图像

1-2)检测输入图像的尺寸是否相同

1-3)进行滤波处理(SDBSF)

第二步、改进变差函数提取新增建设用地

2-1)对滤波后的图像计算变差函数生成纹理特征图

2-2)不同时相的纹理特征图作比值

2-3)确定阈值提取新增建设用地

2-4)对提取的结果进行后处理，得到最终的结果

2-5)精度评价

附图说明

图1，本发明主流程示意图

图2，SDBSF滤波流程图

图3，传统的比值检测模板

图4，改进的比值检测模板

图5，改进后的比值线检测模板

图6，快速选择最优参数流程图

图7，窗口w＝7，步长h＝1，四方向0°，45°，90°，135°的模板

图8，不同时相的两幅RADARSAT2图像

图9，Google Earth对应的Quickbird卫星的光学影像

图10，人工解译出的新增建设用地的真值图像

图11，未改进算法提取结果

图12，改进算法提取结果

具体实施方式

一、图像滤波

具体步骤如图2所示：

1.强点目标标记与保留

SAR图像中的强点目标是SAR图像中一类常见有重要的目标，它们往往对应于一些类似于角反射器的人造地物，这些点对某些特定目标的精确定位和图像同名点的选取有重要意义，因此应予以保留。如果不考虑这些点目标的特殊性，滤波器很容易将其视为噪声而平滑掉。此外，这些强反射点还会严重影响其周边像素的滤波过程，因此，应标记强点目标，使其不参与周边像素的滤波计算。

强点目标的识别方法如下：

采用5*5矩形窗口，若中心像素值与窗口内其他像素值的均值比小于阈值T，则标记为强点目标，保留其原始灰度值不变。遍历每个像素，得到一幅标记有点目标的掩膜图像，标记的所有的点目标即不参与周围像素的结构检测过程，也不参与它们的滤波过程。

2.基于局域统计特性的纹理区与均质区的划分

SAR图像中均质区内与非均质区可以依据局域方差系数进行划分。局域方差系数(标准差/均值)C_ij是有效衡量图像局部的均匀性的指标，常用来反映窗口内的局部灰度特征。当局域方差系数C_ij＞阈值Cu时，视为非均质区。反之，视为均质区。

3.纹理区内的线、边缘及微纹理区的检测及相应滤波模版的选择

比值边缘检测方法是一种具有恒虚警概率(constant false-alarm rate，CFAR)的检测算法，因此适用于SAR图像的应用。传统的比值检测模板如图3。

通过计算4个方向上，边缘两侧区域灰度平均值的比值(r1，r2，r3，r4)与阈值的关系，进行判别。设阈值记为Tt，令rmin＝min(r1，r2，r3，r4)。

当rmin＜Tt时，该中心像素属于边缘区，否则，该中心像素属于均质区。

传统比值检测法的缺点是它只检测了边缘的方向，却没考虑中心像素更接近哪一侧的边缘。因此会造成边缘的定位不准。因此本文在传统方法上进行了改进，改进后的比值边缘检测模板见图4。(其中灰色部分D1像素参与计算过程)。同理计算16个模板的D1与D2两个 区域平均灰度值的比值得到一个边缘检测比值向量(r1，r2，...r16)，求其中的最小值记r_edge＝min(r1，r2，...r16)。

改进后的算法优点有两个：首先原先的四个方向扩展到了8个，更加细化了边缘检测的方向。其次是考虑到中心像素更接近那一侧边缘问题，针对每个方向设计了两个模板，因此能更精确定位边缘。

类似的，将改进后的边缘检测方法扩展到线检测中，得到线检测比值向量(r1，r2，...r8)和线检测最小比值r_line＝min(r1，r2，...r8)。改进后的比值线检测模板见图5。

然后，对纹理区内的像素同时采用边缘检测模板及线检测模板，计算比值r_edge和r_line。设阈值为Tt。

若r_edge＜r_line且r_edge＜Tt，中心像素被视为边缘，此时该中心像素的滤波模板选取比值向量中最小值对应的那个边缘检测模版。

若r_edge＞＝r_line且r_line＜Tt，中心像素被视为线，此时该中心像素的滤波模板选取比值向量中最小值对应的那个线检测模版。

其他情况，则视为微纹理区，此时该中心像素的滤波模板选取M*M的矩形窗口，窗口的大小应小于初始检测窗口，可根据具体图像适当调整。

4.均质区内寻找最大同质区及相应滤波模版的选择

对于均质区域的处理，窗口大小的选取对平滑噪声的效果有至关重要的作用，窗口越大，噪声的抑制效果越好。因此可以采用自适应改变窗口大小和区域增长结合的方法寻找最大的均质区域。

首先采用自适应增加窗口大小的方法，假设通过步骤2计算得出中心像素在均质区内，此时的初始窗口大小为M*M，则增加窗口大小为(M+1)*(M+1)，再次计算局域方差系数C，若C＜阈值T时，继续扩大窗口，直到不满足阈值条件停止。记此时的窗口大小为M’*M’，该窗口为最大均质矩形。

由于实际中的均质区域往往是不规则的形状，因此随后使用区域增长的方法寻找更准确的均质区域。传统的区域增长是基于梯度进行的，同样不具有恒虚警率，不适用于符合乘性斑噪模型的SAR影像，因此用比值的方法代替。以M’*M’窗口内的像素为种子点，根据需要可选用4邻接和8邻接，以8邻接为例，若某像素在种子点的8邻接范围内，且该像素值与种子像素值的比值大于阈值T时，则将该像素合并为种子点，以该种子像素继续进行区域增长，直到不满足阈值条件停止。依次遍历每个种子点，得到区域增长后的一块最大同质区，则该最大同质区内的全部像素参与滤波。

5.使用选择的模板滤波

最后，将选择后的模板对当前像素进行滤波，对于强点目标，保留原值；对于非均质区内的像素，选择MMSE方法进行滤波；对于均质区内的像素选择均值滤波器。遍历整幅影像，最后得到一幅滤波后的影像。

二、计算滤波后影像的变差函数

变差函数理论由数学专家G.Maberon教授1962年创立，作为地质统计学的重要工具，变差函数被应用于空间随机场的统计特性研究。变差函数又称半变差函数(Semivariogram Function)，定义为区域化变量Z(x)和Z(x+h)(同时包含两点距离和方向信息)两点之差的方差之半：

(式1)

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

(式2)

其中，N(h)表示观测数据中间距为h的点对数目，估计值γ*(h)通常称为实验变差函数。变差函数用于度量区域化变量的空间相关性，能充分反映图像数据的随机性和结构性。

变差函数用于纹理分析的常用计算方式是首先确定步长h、窗口w、计算方向，然后计算窗口w内所有间距为h的点对的半方差值，然后取平均(如式(2))作为窗口中心点的变差函数值，遍历全图即得到影像的变差函数特征图。对于高分辨率的SAR图像而言，窗口w的选取取决于h，要保证窗口w内间距h的点对数目足够，窗口w至少应为3h～5h。然而，窗口w过大会造成图像整体模糊、边缘虚警率高，且窗口内取平均的方法易受噪声、孤立强反射点干扰，算法稳健性差。发明根据变差函数计算方法所存在的弊端，提出一种稳健的变差函数计算方法，该算法继承了变差函数和标准中值滤波方法的优点，窗口w取值不再受h约束，主要思想是：在计算像素点的变差函数值时，不再取窗口w内间距为h的点对来求取，而是计算像素点(x，y)与固定方向(如0°方向)上间距h的点(x+h，y)的半方差值，将此值作为点(x，y)的半方差值；然后取窗口内所有像素的半方差值的中值代替均值作为变差函数值。计算公式如下：

(式3)

其中，ω为二维模板，以窗口w＝7，步长h＝1，四方向0°，45°，90°，135°为例，如图所示。

根据图5的流程快速选择最优参数。

分别计算两幅图像的变差函数，生成各自的纹理图像。之后以新时相数据的变差函数图与旧时相数据的变差函数图作比值，得到一副比值图，并根据阈值分割算法得到一个划分阈值a，从而根据a对比值图像进行阈值分割，最终得到新增建设用地的初步提取结果。

由于初步提取的结果比较破碎，并且有许多误提的小斑块，由实际经验可知这种小斑块并非建设用地，可能是传感器成像角度不同引起的灰度变化而被算法误提出的小斑块。因此进行后处理，膨胀与腐蚀连接建设用地区域的外部轮廓，去除小面积区域删除误提取的碎小斑块。

精度验证

测试数据采用时相分别为2009年8月6日、2011年6月30日的浙江省杭州市北部城区的RADARSAT2影像，影像具体参数信息如表1。两幅影像首先进行了几何配准和斑点噪声抑制、亮度归一化等预处理工作，处理后的图像均为1000*1000像素大小，256级灰度。利用Google Earth上对应的Quickbird卫星的光学影像作为精度评价数据。图8，图9即为经过图像配准后的两时相SAR影像及Quickbird影像，图10为人工解译出的新增建设用地的真值图像。

表1影像参数信息表

从结果可以看出变化区域大致被检测出来，整体精度较高，漏警率较小。两种方法作对比，改进方法在整体精度及虚警率都优于未改进方法。但是整体的检测结果要比真实变化区域范围大，而且包含很多虚警区域，即虚警率比较大。这主要是由于两幅影像本身差异以及入射角相差较大，在一些未变化的区域两幅影像在灰度上表现出差异，因此造成结果的虚警率大。

方法	总体精度/％	虚警率/％	漏警率/％
				未改进方法	75.8	35.5	24.2
改进方法	81.1	26.9	18.9

表2变化检测精度表。

Claims (1)

1.一种基于高分SAR新增建设用地检测方法，其步骤为：

第一步、数据预处理

1-1)输入不同时相的高分辨率SAR图像

1-2)检测输入图像的尺寸是否相同

若相同进行滤波处理，若不同则终止程序

1-3)进行滤波处理(SDBSF)。首先进行强点目标检测，对于强点目标滤波算法保持其原值不作处理，对于非强点目标初始窗口大小为N，计算局域方差系数C_ij(标准差/均值)，设置阈值Cu，当C_ij≥Cu时，认定窗口区域为同质区域，进一步寻找最大同质区域大小为N×N矩形，然后进行区域增长，选择最终参与滤波的像素，进行均值滤波；当C_ij＞Cu时，认定窗口区域为非同质区域，进行边缘与线检测，对于微纹理区选择N×N矩形窗口计算滤波参数b进行MMSE滤波，对于边缘选择基于边缘的滤波窗口计算滤波参数b进行MMSE滤波，同样对于线选择基于线的滤波窗口计算滤波参数b后进行MMSE滤波，最终输出滤波后的影像。

第二步、改进变差函数提取新增建设用地

2-1)对滤波后的图像计算变差函数生成纹理特征图

变差函数 $\mathrm{\γ}*\left(h\right)=med\{\frac{1}{2}\[Z\left(X\right)-Z(X+h),X\∈\mathrm{\ω}\]\}$

其中，ω为二维模板，本发明采用0°，45°，90°，135°四个方向，在计算像素点的变差函数值时，计算像素点(x，y)与模板四个方向上间距h的点(x+h，y)的半方差值，将此值作为点(x，y)的半方差值；然后取窗口内所有像素的半方差值的中值代替均值作为变差函数值。

2-2)不同时相的纹理特征图作比值

以新时相的图像除以老时相的图像，老时相图像可能存在为0的数据点，在相除之前将所有为0的数据点赋值为0.000001，相除后的结果进行归一化至0-255之间。

2-3)确定阈值提取新增建设用地

根据经验取得阈值a＝200，当归一化后的数据点值大于a时，认为此数据点为新增建设用地的数据点。

2-4)对提取的结果进行后处理，得到最终的结果

搜索初步提取结果中所有的联通区，通过膨胀腐蚀来连接破碎的新增建设用地边界，并通过填洞获取完整的建筑区，最后利用删除去掉非新增建设用地的小面积区域。

2-5)精度评价

可以通过输入真值进行提取结果的精度评价。

检测率：DR＝TP/(TP+FN)*100

虚警率：Far＝FP/(TP+FP)*100

漏警率：MAR＝FN/(TP+FN)*100

其中TP为正确提取的像元数量，FN为真值中是新增建设用地但程序未提取的像元数量，FP为真值中是非建设用地但程序错误提取的像元数量。