CN106296655A - 基于自适应权值和高频阈值的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于自适应权值和高频阈值的SAR图像变化检测方法，主要解决现有技术检测效果不理想，单一类型差异图检测精度低的问题。其实现步骤为：1.读入不同时刻从同一地区获取的两幅SAR图像I₁和I₂；2.分别计算两幅图像I₁和I₂对应像素的均值比差异图D_m和邻域对数比差异图D_l；3.对均值比差异图D_m和邻域对数比差异图D_l进行小波融合，得到融合后的差异图像X_d；4.将融合差异图X_d聚成不同的两类，得到变化检测结果。本发明操作简单、抗噪性好、检测精度高，对不同类型的SAR图像均能获得较好的效果，可应用于环境监测、海洋观测、灾害评估、资源勘察、城区规划和地理测绘。

Description

基于自适应权值和高频阈值的SAR图像变化检测方法

技术领域

本发明属于雷达图像处理技术领域，特别涉及一种SAR图像的变化检测方法，可用于环境监测、海洋观测、灾害评估、资源勘察、城区规划和地理测绘。

背景技术

合成孔径雷达SAR出自20世纪50年代，作为一种主动式相干微波成像雷达，相比可见光和红外遥感，微波遥感具有得天独厚的优点：一、合成孔径雷达采用侧视方式成像，覆盖面积大，分辨率高；二、受气候等外界自然环境影响小，能全天候24小时工作；三、具有很强的对地穿透能力，能有效地识别伪装和穿透掩盖物。由于SAR图像具备的这些无法比拟的优点，加之近年来星载、机载和无人机载合成孔径雷达成像技术的飞速发展，使得SAR图像日益成为变化检测研究的重要数据源。SAR图像的变化检测是通过对同一地域不同时间所获取的两幅或者多幅不同图像的对比分析，根据图像之间的差异信息来得到所需地物或目标的变化情况的过程。该技术可以检测出不同时期图像灰度值或局部纹理之间的变化，在此基础上获得所需目标在形状、位置、数量及其它属性方面的变化情况，在民事领域和军事领域都得到了广泛应用。

一个完整的SAR变化检测过程包括：(1)图像预处理(2)变化区域检测(3)精度评估，其中变化区域检测是核心的步骤。目前，最常见的SAR图像变化检测方法是基于差异图像分析的变化检测方法，其变化区域检测研究通常包括2个关键的步骤：即差异图的构造和变化信息的提取,其中差异图构造的质量高低对变化检测的结果有着决定性的影响。近年来，以差值法、比值法等基于代数运算的差异图构造方法颇受广大研究者的喜爱，通常情况下，由于差值法构造差异图时易受SAR图像本身质量和相干斑乘性噪声的影响，产生“伪变化”信息，因此多采用比值法构造差异图，其中以对数比值法和均值比值法应用最为广泛。Dekker于1997年在“Speckle filtering in satellite SAR change detection imagery”一文中提出使用对数比的方法构造差异图，对数比法就是对比值法取对数，它除了可以将乘性噪声转化为加性噪声外，还对比值图像的变化范围进行了压缩，因此削弱了高像素值区域的变化，部分真实变化区域轮廓的像素灰度值与非变化类像素的灰度值比较接近,有效抑制了背景信息，即非变化区域，使得对数比值图的背景区域比较平坦。但是该方法对比值图像变化范围的压缩，不利于保留变化区域的轮廓，不能在最大程度上反应变化区域的真实情况。

图像的均值比值法是假设图像中的变化会引起局部均值的变化，在对应像素灰度值相除之前加入了邻域信息，起到平滑噪声的作用，使得差异影像图的变化区域与实际情况相接近。但是由于加入了邻域信息，其背景区域比较粗糙，存在过多的伪变化信息，不利于差异图的结果分类，而且对于噪声较大的图像，均值比法的抑噪能力有限。因此，使用单一类型的差异图进行变化检测会导致检测精度低以及适用范围窄等问题。

近年来，将单一类型的差异图进行融合的策略在变化检测研究中起到至关重要的作用。对SAR图像变化检测来说，图像融合是从两幅或者更多的图像中结合相关信息而得到一幅比任何一幅输入图像信息都多的图像的过程，目的是为了构造比单一差异图质量更高的融合差异图，从而提高变化检测的精度。其中利用小波变换进行图像融合的方法得到了广泛应用，小波变换应用于图像融合的优势在于：小波变换的紧支性、对称性和正交性赋与它良好的图像融合性能，克服了局部相邻图像信息弱联系的缺点，而且还反映了原始图像的局部变化，可以将图像分解到不同的频域，在不同的频域运用不同的融合规则，得到合成图像的多分辨率分析，从而在合成图像中保留原图像在不同频域的特征，另外小波变换还有较低的计算复杂度。一幅图像经3层平稳小波变换之后，被分解为1个低频子带系数和3个高频子带系数，其中低频子带系数包含了图像的轮廓信息，能够充分体现出变化区域的信息；高频子带系数包含了图像的细节部分，对SAR图像来说，由于乘性噪声是随机分布的，因此高频子带系数也包含了大量噪声。

发明内容

本发明目的在于提出一种基于自适应权值和高频阈值的SAR图像变化检测方法，以解决现有技术检测效果不理想，单一类型差异图检测精度低、适用范围窄的问题，更好地检测出SAR图像的变化区域。

为了实现上述目的，本发明的技术方案包括如下：

(1)读入不同时刻从同一地区获取的两幅SAR图像I₁和I₂,并计算这两幅图像对应像素的均值比差异图D₀；

(2)对均值比差异图D₀进行归一化，得到归一化后的均值比差异图D₁；

(3)将归一化后的均比值差异图D₁进行3*3窗口的中值滤波，抑制均值比图中的孤立像素点，得到中值滤波后的均值比差异图D_m；

(4)计算两幅输入图像I₁和I₂对应像素的邻域对数比差异图D₂；并对该邻域对数比差异图D₂进行归一化，得到归一化后的邻域对数比差异图D_l；

(5)根据均值比差异图D_m和邻域对数比差异图D_l，获得融合后的差异图像X_d：

5a)对均值比差异图D_m和邻域对数比差异图D_l分别进行3层平稳小波变换，得到每一幅差异图的1个低频子带系数和3个高频子带系数；

5b)对2幅差异图的低频子带系数按照最大-最小自适应权值策略进行融合，生成1个融合后的低频子带系数；

5c)对2幅差异图的高频子带系数按照局部能量最小原则进行融合，生成3个融合后的高频子带系数；

5d)对融合后的3个高频子带系数进行自适应阈值滤波操作；

5e)对融合后的1个低频子带系数和3个经自适应阈值滤波后的高频子带系数进行平稳小波逆变换，得到融合后的差异图像X_d；

(6)用k-means聚类方法将融合后的差异图X_d聚类成不同的两类，分别计算这两个不同类别的均值，定义均值较大的那一类为变化类，均值较小的那一类为非变化类，得到最终的变化检测结果。

本发明与现有技术相比具有以下优点:

1、本发明采用基于小波变换的图像融合方法构造差异图，融合后的差异图综合利用了均值比差异图的变化区域接近实际情况的特点和对数比差异图的背景区域相对平滑的特点，达到了抑制背景区域信息的同时增加变化区域信息的目的，解决了单一类型差异图检测精度低、适用范围窄的问题,用此方法构造的融合差异图的质量要好于单幅差异图。

2、本发明采用最大-最小自适应权值策略对两幅差异图的低频分量进行融合，该融合权值能根据不同的像素点自适应调整，无需手动调整或者人工设定，操作简单。

3、本发明采用自适应阈值滤波方法对融合后的高频子带系数进行去噪处理，能够按照高频系数所处的位置自主调节阈值，从而实现去噪的同时又在最大程度上保留细节信息的目的，有利于提高融合差异图的抗噪能力。

实验表明，本发明对多幅SAR图像都能得到满意的变化检测结果，具有较好的鲁棒性。

附图说明

图1是本发明的实现流程图；

图2是Bern地区SAR图像数据集和该数据集的标准变化检测图；

图3是本发明和对比方法对Bern地区SAR图像数据集的变化检测结果图；

图4是Ottawa地区SAR图像数据集和该数据集的标准变化检测图；

图5是本发明和对比方法对Ottawa地区SAR图像数据集的变化检测结果图；

图6是黄河入海口地区SAR图像数据集和该数据集的标准变化检测图；

图7是本发明和对比方法对黄河入海口地区SAR图像数据集的变化检测结果图。

具体实施方式

参照图1，本发明的具体步骤如下：

读入不同时刻从同一地区获取的两幅SAR图像I₁和I₂，并计算这两幅图像对应像素的均值比差异图D₀。

同一地区不同时刻的两幅SAR图像I₁和I₂通常是由机载或者星载体SAR对某一固定地区按照不同的时间间隔采集所得到，由于均值比值法在对应像素灰度值相除之前加入了邻域信息，起到平滑噪声的作用，所得到的差异影像图的变化区域与实际情况相接近，因此本发明采用均值比值法构造第一幅差异图，其构造方法按如下步骤进行：

1a)对第一幅输入图像I₁在像素(i,j)处的灰度值I₁(i,j)与第二幅输入图像I₂在像素(i,j)处的灰度值I₂(i,j)进行均值比运算，得到均值比矩阵D₀在像素(i,j)处的值d₀(i,j)：

$d_0(i,j)=1-min(\frac{{\mathrm{\μ}}_1(i,j)}{{\mathrm{\μ}}_2(i,j)},\frac{{\mathrm{\μ}}_2(i,j)}{{\mathrm{\μ}}_1(i,j)})$

上式中μ₁(i,j)和μ₂(i,j)分别表示第一幅输入图像I₁和第二幅输入图像I₂中以像素(i,j)为中心的3*3邻域窗口内所有像素灰度值的平均值；

1b)遍历两幅输入图像所有对应的像素，得到均值比差异图D₀＝{d₀(i,j)}；

1c)对均值比差异图D₀中每一像素按照如下公式进行归一化,得到均值比差异图D₀在像素(i,j)处归一化后的值d₁(i,j)；

$d_1(i,j)=\frac{d_0(i,j)-min\left(D_0\right)}{max\left(D_0\right)-min\left(D_0\right)},$

其中，d₁(i,j)表示均值比差异图D₀中的像素d₀(i,j)被归一化后的值；

1d)遍历均值比差异图D₀中的所有像素，得到归一化后的均值比差异图D₁＝{d₁(i,j)}。

步骤2：对均值比差异图D₁进行3*3窗口的中值滤波，得到中值滤波后的均值比差异图D_m。

对均值比差异图D₁中的每一个像素点，选取一个大小为3*3的矩形窗口，将窗口中全部像素点的灰度值按从大到小的顺序排列，组成一个灰度序列，选取序列中间位置的灰度值作为滤波后的输出值，从而抑制均值比图中的孤立像素点，得到中值滤波后的均值比差异图D_m。

步骤3：计算两幅输入图像I₁和I₂对应像素的邻域对数比差异图D₂，并对该邻域对数比差异图D₂进行归一化，得到归一化后的邻域对数比差异图D_l。

对数比值法对比值图像的变化范围进行了压缩，因此削弱了高像素值区域的变化，部分真实变化区域轮廓的像素灰度值与非变化类像素的灰度值比较接近,有效抑制了背景信息，即非变化区域，使得对数比值图的背景区域比较平坦,鉴于以上优点，本发明采用基于邻域的对数比值法构造第二幅差异图。基于邻域的对数比值法首先将两幅输入图像对应象素点的3*3局部邻域内所有像素进行高斯滤波，滤除该局部邻域内中的乘性噪声,由于SAR图像噪声类型复杂，分布随机，高斯滤波操作并不能完全过滤噪声，因此，再通过对数比值法，把乘性噪声转换为加性噪声，进一步滤除部分加性噪声,最后对3*3局部邻域内的元素值取均值，消除图像中孤立的噪声点。经过以上操作步骤，使该方法具备很高的抗噪能力，同时能够有效抑制背景信息即非变化区域，所构造的差异图的背景区域比较平坦。本步骤的具体实现如下：

3a)对第一幅输入图像I₁中以像素(i,j)为中心的3*3局部邻域l₁(m,n)进行大小为3，标准偏差为5的高斯滤波操作,得到滤波后的局部邻域l₁'(m,n)；对第二幅输入图像I₂中以像素(i,j)为中心的3*3局部邻域l₂(m,n)进行大小为3，标准偏差为5的高斯滤波操作,得到滤波后的局部邻域l'₂(m,n)，其中m是局部邻域内的像素的横坐标，n是局部邻域内的像素的纵坐标；由于高斯滤波器是根据高斯函数的形状来选择权值的线性滤波器，对去除服从正态分布的噪声是很有效果的，而SAR图像中本身就存在许多服从正态分布的乘性噪声，因此，采用高斯滤波器进行滤波后，能进一步降低噪声的负面影响，高斯滤波的具体操作是：用一个模板扫描图像中的每一个像素，然后用模板确定的邻域内像素的加权平均灰度值去替代模板中心像素点的值；

3b)按照下式，计算局部邻域l₁'(m,n)和局部邻域l'₂(m,n)内对应像素灰度值的商的对数l′_log(m,n)：

$l_{\log }^{\′}(m,n)=\log \frac{\left(l_2^{\′}\right(m,n)+1)}{\left(l_1^{\′}\right(m,n)+1)}=\log \left(l_2^{\′}\right(m,n)+1)-\log \left(l_1^{\′}\right(m,n)+1),$

3c)遍历局部邻域l₁'(m,n)和局部邻域l'₂(m,n)中所有对应像素，得到局部邻域对数比矩阵Z_log(m,n)＝{l′_log(m,n)}；

3d)按照下式，计算局部邻域对数比值矩阵Z_log(m,n)的均值d₂(i,j)，将其作为邻域对数比差异图D₂在像素(i,j)处的值：

$d_2(i,j)=\left|\frac{1}{N}\underset{(m,n\∈B)}{\Σ}{l}_{\log }^{\′}(m,n)\right|$

其中B是以像素(i,j)为中心的3*3局部邻域窗口的大小，N是局部邻域窗口内的元素数量；

3e)遍历两幅输入图像I₁和I₂中所有对应的像素，得到邻域对数比差异图D₂＝{d₂(i,j)}。

3f)对邻域对数比差异图D₂中每一像素按照如下公式进行归一化，得到邻域对数比差异图D₂在像素(i,j)处归一化后的值d_l(i,j)；

$d_l(i,j)=\frac{d_2(i,j)-min\left(D_2\right)}{max\left(D_2\right)-min\left(D_2\right)}$

其中，d_l(i,j)表示邻域对数比差异图D₂中的像素d₂(i,j)被归一化后的值；min(D₂)表示取邻域对数比差异图D₂中灰度值最小的像素值；max(D₂)表示取邻域对数比差异图D₂中灰度值最大的像素值；

3g)遍历邻域对数比差异图D₂中所有像素，得到归一化后的邻域对数比差异图D_l＝{d_l(i,j)}。

步骤4：对均值比差异图D_m和邻域对数比差异图D_l进行小波融合。

由于均值比差异图所反映的变化区域与实际情况相接近，但是背景区域却很粗糙，而对数比差异图所反映的背景区域比较平坦，但是变化区域往往与真实情况不相符。因此，采取小波融合的方法将两幅差异图进行融合，综合利用均值比差异图的变化区域接近实际情况的特点和对数比差异图的背景区域相对平滑的特点，达到抑制背景区域信息的同时增加变化区域信息的目的，构造比单幅差异图质量更高的差异图。其具体操作过程如下：

4a)对均值比差异图D_m进行3层平稳小波分解：

在小波分解的第一层用均值比差异图D_m与小波基做内积，再经过水平和垂直方向的两次2倍间隔采样后，均值比差异图D_m就被分解为4个具有不同分辨率的图像，包括1个低频子带图像和水平、垂直、对角线方向的3个高频子带图像；

在小波分解的第二层对第一层分解所得到的1个低频子带图像与小波基做内积，再经过水平和垂直方向的两次2倍间隔采样后，第一层的低频子带图像又被分解为4个具有不同分辨率的图像，包括1个低频子带图像和水平、垂直、对角线方向的3个高频子带图像；

在小波分解的第三层对第二层分解所得到的1个低频子带图像与小波基做内积，再经过水平和垂直方向的两次2倍间隔采样后，第二层的低频子带图像再被分解为4个具有不同分辨率的图像，包括1个低频子带图像，水平、垂直和对角线方向各一的个高频子带图像；

经过以上三层小波分解后，就得到均值比差异图D_m在第3分解层上的低频子带系数在第1个分解层上的3个方向的高频子带系数在第2个分解层上的3个方向的高频子带系数在第3个分解层上的3个方向的高频子带系数其中d＝(HH,HV,HD)表示高频子带的方向，HH表示水平方向，HV表示垂直方向，HD表示对角方向；m表示均值比差异图，L表示低频子带系数，mH表示均值比差异图的高频子带系数；

4b)对邻域对数比差异图D_l进行3层平稳小波分解：

在小波分解的第一层用邻域对数比差异图D_l与小波基做内积，再经过水平和垂直方向的两次2倍间隔采样后，邻域对数比差异图D_l就被分解为4个具有不同分辨率的图像，包括1个低频子带图像和水平、垂直、对角线方向的3个高频子带图像；

在小波分解的第二层对第一层分解所得到的1个低频子带图像与小波基做内积，再经过水平和垂直方向的两次2倍间隔采样后，第一层的低频子带图像又被分解为4个具有不同分辨率的图像，包括1个低频子带图像和水平、垂直、对角线方向的3个高频子带图像；

在小波分解的第三层对第二层分解所得到的1个低频子带图像与小波基做内积，再经过水平和垂直方向的两次2倍间隔采样后，第二层的低频子带图像再被分解为4个具有不同分辨率的图像，包括1个低频子带图像，水平、垂直和对角线方向各一的个高频子带图像；

经过以上三层小波分解后，就得到邻域对数比差异图D_l在第3分解层上的低频子带系数在第1个分解层上的3个方向的高频子带系数在第2个分解层上的3个方向的高频子带系数在第3个分解层上的3个方向的高频子带系数其中d＝(HH,HV,HD)表示高频子带的方向，HH表示水平方向，HV表示垂直方向，HD表示对角方向；l表示邻域对数比差异图，L表示低频子带系数，lH表示邻域对数比差异图的高频子带系数；

4c)按下式对均值比差异图D_m的低频子带系数和邻域对数比差异图D_l的低频子带系数进行融合，得到融合后的低频子带系数

$D_L^{ml}=k\×D_L^m+(1-k)\×D_L^l$

上式中，和分别表示均值比差异图D_m和邻域对数比差异图D_l经3层平稳小波分解后的低频子带系数，表示融合后的低频系数。k是最大-最小自适应权值，其构造原理如下所述：

均值比差异图D_m和邻域对数比差异图D_l经3层平稳小波分解后的低频子带系数均反映了图像的轮廓信息，能够体现出变化区域的信息，而均值比值法差异图所反映的变化区域相比对数比值法差异图更接近于实际情况，因此，在低频融合时就给均值比差异图的低频分配最大权重，而给对数比差异图的低频分配最小权重，这样的融合策略可以确保最后的融合差异图所反映的变化区域更接近于实际情况。

为满足上述权重分配的理论条件并确保权重值能够自适应，本发明按如下步骤构造最大-最小自适应权值k：

首先，从均值比差异图D_m和邻域对数比差异图D_l的对应象素中取出灰度值最小的象素值和灰度值最大的象素值；

其次，用最小象素值与最大象素值做商得到商值；

然后，用1与商值做差得到差值，这样做的目的是使两幅差异图所分配的权值之和永远为1，从而保证了合理的权值分配比；

最后，从商值和差值中选择最大值k作为均值比差异图低频子带系数的融合权值，选择最小值(1-k)作为邻域对数比差异图低频子带系数的融合权值，

该最大-最小自适应权值k的数学表达式如下：

$k=max(\frac{\min (D_m,D_l)}{\max (D_m,D_l)},1-\frac{\min (D_m,D_l)}{\max (D_m,D_l)}),$

上式中min(D_m,D_l)表示从均值比差异图D_m和邻域对数比差异图D_l的对应象素值中取最小值，max(D_m,D_l)表示从均值比差异图D_m和邻域对数比差异图D_l的对应象素值中取最大值。

通过以上操作，就实现了给均值比差异图的低频分配最大权重，而给对数比差异图的低频分配最小权重的目的，使得融合差异图所反映的变化区域更接近于实际情况。

4d)按局部能量系数最小原则的方法对均值比差异图D_m的3个高频子带系数和邻域对数比差异图D_l的3个高频子带系数进行融合：

考虑到邻域对数比值法差异图的背景信息相对均值比差异图要平滑，采用局部能量最小原则的方法可以合并来自两幅差异图高频部分的同质区域，当把邻域对数比差异图的高频部分融入到均值比差异图中后，就会在一定程度上抑制背景信息，使得融合差异图的背景区域保持平滑。局部能量系数最小原则方法的数学表达式如下：

$D_{d,n}^H(i,j)=\left\{\begin{array}{c}{D}_{d,n}^{mH}(i,j),E_{d,n}^m(i,j)<E_{d,n}^l(i,j)\\ D_{d,n}^{lH}(i,j),E_{d,n}^m(i,j)\&GreaterEqual;E_{d,n}^l(i,j)\end{array}\right.,$

上式中，n＝(1,2,3)表示分解的层数；d＝(HH,HV,HD)分别表示差异图经3层平稳小波分解后高频子带系数的方向，其中HH表示水平方向、HV表示垂直方向、HD表示对角方向；表示融合后的高频子带系数在像素(i,j)处的值；和分别表示均值比差异图和邻域对数比差异图经3层平稳小波分解后的高频子带系数在像素(i,j)处的值；表示均值比差异图的高频子带系数中以像素(i,j)为中心的局部能量系数值，按如下公式计算：

$E_{d,n}^m(i,j)=\underset{s\&Element;N_{i,j}}{\&Sigma;}{\&lsqb;D_{d,n}^{mH}\left(s\right)\&rsqb;}^2,$

上式中N_i,j表示以像素(i,j)为中心的3*3局部邻域窗口的大小，表示均值比差异图的高频子带系数在邻域窗口N_i,j内第S个系数的值；

表示在邻域对数比差异图的高频子带系数中以像素(i,j)为中心的局部能量系数值，按如下公式计算：

$E_{d,n}^l(i,j)=\underset{s\&Element;N_{i,j}}{\&Sigma;}{\&lsqb;D_{d,n}^{lH}\left(s\right)\&rsqb;}^2,$

上式中N_i,j表示以像素(i,j)为中心的3*3局部邻域窗口的大小，表示邻域对数比差异图的高频子带系数在邻域窗口N_i,j内第S个系数的值；

4e)对融合后的高频子带系数进行自适应阈值滤波：

图像经3层平稳小波分解后的高频子带系数包含的是图像的细节部分，但是由于SAR图像的相干斑乘性噪声随机分布，致使高频系数中也包含大量的噪声，为了尽可能消除这些噪声对融合差异图的影响，对融合后高频子带系数进行自适应阈值滤波操作。小波变换后的高频系数在图像细节和边缘区域的系数能量集中，幅值较大，所以在这些区域内系数的绝对值之和也就大，而由于噪声的分布是随机的，噪声能量分布相对分散，幅值较小，则在噪声区域内系数绝对值之和就小。因此，要实现去噪和保留细节信息的双重目的，只需要对高频子带的每个系数在细节和边缘区域通过减小阈值来保留更多系数，对噪声区域通过增加阈值来去掉更多噪声。为此，本发明引用Kang Li等在文章“AdaptiveShrinkage for Image Denoising Based on Contourlet Transform”中的提出的自适应收缩阈值法进行去噪，其表达式如下：

$D_{d,n}^H(i,j)=\left\{\begin{array}{c}{D}_{d,n}^H(i,j);D_{d,n}^H(i,j)\&GreaterEqual;T\\ 0;D_{d,n}^H(i,j)<T\end{array}\right.,$

其中，T表示自适应阈值参数，其值按下式计算：

$T=(\mathrm{\&lambda;}-e^{\left(\frac{mean\left(A\right)-\min \left(A\right)}{A(i,j)-\min \left(A\right)}\right)\&times;{\log }_{10}(\mathrm{\&lambda;}-1)})\&times;(c\&times;\mathrm{\&sigma;})$

其中λ是介于1与2之间的常数值；C是噪声能量分布标准差的调和系数，取值为3；σ是融合后高频子带系数的噪声估计值，其值按下式计算：

$\mathrm{\&sigma;}=\frac{median\left(D_{d,n}^H\right(i,j))}{0.6745},$

A(i,j)表示不同方向高频子带系数在局部邻域内的绝对值之和的均值，其值按下式计算：

$A(i,j)=\frac{1}{N}\underset{(i,j\&Element;B)}{\&Sigma;}\left|D_{d,n}^H\left(s\right)\right|,$

其中B是以像素(i,j)为中心的3*3局部邻域窗口的大小，N是局部邻域窗口内的元素数量，表示融合后的高频子带系数在邻域窗口B内第S个系数的值。

从自适应阈值参数T的计算公式中可以看出其值是随着融合后高频系数所处的位置而变化的：当融合后高频系数处于细节和边缘区域时，其系数能量集中，幅值较大，所以在这些区域内系数的绝对值之和也就大，A(i,j)的值也增大；同理，当融合后高频系数处于噪声区域时，其系数能量分散，幅值较小，A(i,j)的值也相应减小。由于λ介于1到2之间时，由对数函数特性可知，参数T的计算公式的分子为一负常数，分母是随着A(i,j)变化的负指数函数，按照指数函数的特性可知，其取值介于(0，1)之间，当融合后高频子带系数处于细节和边缘区域时，A(i,j)增大，指数函数值趋向1，T的值减小，即减小阈值，刚好能满足边缘细节区域减小阈值来保留更多系数的条件；同理，当高频系数处于噪声区域时，A(i,j)减小，T值增大，即增大阈值，刚好能满足噪声区域增加阈值来去除更多噪声的条件，从而实现即去除噪声又在最大程度上保留边缘和细节信息的目的。

步骤5：对融合后的低频子带系数和滤波后的高频子带系数进行平稳小波反变换，得到融合后的图像X_d。

融合生成的1个低频子带系数所包含的变化区域的信息接近于实际情况，融合后的3个高频子带系数经过自适应阈值滤波处理后去除了包含在两幅差异图像中的噪声，因此，将1个低频子带和3个高频子带通过平稳小波反变换所得到的融合差异图能更好的反应变化区域的真实情况。

步骤6：用k-means聚类方法对融合差异图X_d进行聚类，得到最终变化检测结果。

6a)从融合差异图X_d中随机选择2个像素点作为原始聚类中心；

6b)分别计算融合差异图X_d中每一像素点(i,j)与两个原始聚类中心的距离，并根据最小距离标准重新对相应像素点进行划分，把(i,j)归到与两个原始聚类中心距离最近的那一类中去；

6c)计算归类后将每一类的均值作为新的聚类中心；

6d)重复步骤6b)-6c)，直到每个聚类中心不再发生变化为止，得到两个类别的均值；

6e)将两个类别中均值小的那一类像素转化为0，即为未变化类，均值大的那一类像素转化为255，即为变化类，得到最终的变化检测结果。

本发明的效果可以通过以下实验结果进一步说明：

1.实验条件

实验环境为：windows 7，Service Pack 1，CPU：Intel(R)Core(TM)i3M380，基本频率2.53GHZ，软件平台为MatlabR2012a。

第一个数据集为瑞士Bern地区SAR图像数据，该数据集反映的是瑞士首都伯尔尼市郊区附近的水灾情况，分别于1999年4月和1999年5月通过ERS-2获得的Bern地区的前后两时相SAR图像，两幅图像的尺寸均为301×301像素，灰度级为256，如图2所示，其中，图2(a)是洪水侵袭之前的图像，图2(b)是洪水侵袭之后的图像，图2(c)为该数据集标准变化检测结果图，包括1269个变化像素点，89332个未变化像素点。

第二个数据集为渥太华地区SAR图像数据，该数据集是用RADARSAT型号卫星的星载SAR传感器分别于1997年5月和1997年8月拍摄的渥太华地区的前后两时相SAR图像，两幅图像的尺寸均为290×250，灰度级为256，如图4所示，其中，图4(a)是洪水侵袭之前的图像，图4(b)是洪水侵袭之后的图像，图4(c)为该数据集的标准变化检测结果图，包括85451个未变化像素点和16049个变化像素点。

第三个数据集是中国黄河入海口地区的SAR图像中截取的部分图像数据，该数据集是是由Radarsat-2型号卫星的星载SAR传感器于2008年6月和2009年6月拍摄的黄河入海口地区的前后两时相SAR图像，两幅图像的尺寸均为291×444，灰度级为256，如图6所示，其中，图6(a)为2008年6月的图像，图6(b)为2009年6月的图像，图6(c)为标准变化检测结果图，包括7295个变化像素点和121909个未变化像素点。

2.实验评价指标

实验使用的评价指标是漏检数FN、误检数FP、总错误数OE、正确检测率PCC和Kappa系数。其中，漏检数为实际发生了变化而没有检测出来的像素的总和；误检数为实际没有发生变化但被当作变化所检测出来的像素的总和；总错误数是漏检数和误检数之和。正确检测率PCC按下式计算：

$PCC=1-\frac{OE}{N}$

上式中N为总像素数，OE为总错误数。

Kappa系数按下式计算:

$Kappa=\frac{PCC-PR}{1-PR}$

其中Nc为总变化数，Nu为总未变化数，TP为正确检测出的变化数，TN为正确检测出的未变化数，TP为总变化数与漏检数的差，TN为总未变化数与误检数的差。

3.实验内容

用本发明方法与两种基于单一类型的差异图进行变化检测的方法和现有的三种基于图像融合的变化检测方法在3个图像数据集上进行对比。其中，第一种基于单一类型差异图的变化检测对比方法记为均值比值法，其主要步骤为：

首先，采用均值比值运算构造均值比差异图；

其次，对均值比差异图进行高斯滤波，滤除差异图中的噪声；

最后，用K-means聚类方法将滤波后的均值比差异图聚成变化类和非变化类，得到变化检测结果图；

第二种基于单一类型差异图的变化检测对比方法记为对数比值法，其主要步骤为：

首先，采用对数比值运算构造对数比差异图；

其次，对对数比差异图进行高斯滤波，滤除差异图中的噪声；

最后，用K-means聚类方法将滤波后的对数比差异图聚成变化类和非变化类，得到变化检测结果图；

另外3种基于图像融合的变化检测的对比方法分别为:1)Yaoguozheng在文章“Using Combined Difference Image and K-means Clustering for SAR Image ChangeDetection”中提出的方法，记做CDI-K法；2)Jingjing Ma在文章“Wavelet Fusion onRatio Images for Change Detection in SAR Images”中提出的方法，记做WFD-FL法；3)Mao guo Gong在文章“Change Detection in Synthetic Aperture Radar Images basedon Image Fusion and Fuzzy Clustering”中提出的方法，记做WFD-RFL法。

实验1

用本发明方法和所述的现有5种变化检测对比方法对第一个SAR图像数据集中的图2(a)和图2(b)进行变化检测，实验结果如图3所示，其中

图3(a)是用均值比值法的变化检测结果图，

图3(b)是用对数比值法的变化检测结果图，

图3(c)是用CDI-K法的变化检测结果图，

图3(d)是用WFD-FL法的变化检测结果图，

图3(e)是用WFD-RFL法的变化检测结果图，

图3(f)是本发明的变化检测结果图。

实验2

用本发明方法和所述的现有5种变化检测对比方法对第二个SAR图像数据集中的图4(a)和图4(b)进行变化检测，实验结果如图5所示，其中

图5(a)是用均值比值法的变化检测结果图，

图5(b)是用对数比值法的变化检测结果图，

图5(c)是CDI-K法的变化检测结果图，

图5(d)是WFD-FL法的变化检测结果图，

图5(e)是WFD-RFL法的变化检测结果图，

图5(f)是本发明的变化检测结果图。

实验3

用本发明方法和所述的现有5种变化检测对比方法对第三个SAR图像数据集中的图6(a)和图6(b)进行变化检测，实验结果如图7所示，其中

图7(a)是用均值比值法的变化检测结果图，

图7(b)是用对数比值法的变化检测结果图，

图7(c)是用CDI-K法的变化检测结果图，

图7(d)是用WFD-FL法的变化检测结果图，

图7(e)是用WFD-RFL法的变化检测结果图，

图7(f)是用本发明的变化检测结果图。

对本发明与现有5种对比方法在3个SAR图像数据集上的变化检测结果的评价指标数值进行统计，结果如表1。

表1本发明与现有5种对比方法的变化检测结果评价指标数值统计

4.实验结果分析

从图3、图5、图7和表1可以看出：本发明对3个SAR图像数据集进行变化检测所得到检测结果图在直观效果上与3个数据集的标准检测图2(c)、图4(c)和图6(c)最接近，而基于均值比差异图、对数比差异图等传统单一类型差异图的变化检测结果图含噪多，细节和边缘模糊不清晰，特别是图3(a)和图7(a)基本看不出变化区域，与标准检测图差距很大；CDI-K方法在3个数据集上的检测结果图在直观效果上均漏检了部分像素，变化区域不够完整；WWFD-FL方法和WFD-RFL方法在3个数据集上的检测结果图在直观效果上都存在边缘区域和细节信息不清晰，噪声多和漏检像素多的现象，特别是WWFD-FL方法的检测结果图基本看不出变化区域的边缘和细节，只能看清变化区域的轮廓，与标准检测图差距很大；

表1也进一步证明本发明对前3个SAR图像数据集进行变化检测所得到的变化检测结果在总错误数和正确检测率两个评价指标上均优于现有的5种对比方法，Kappa系数评价指标在前两个数据集上均优于现有的5种对比方法，在第三个数据集上与WFD-RFL方法基本持平。

综上，本发明采用小波融合的方法构造差异图，综合利用了两种不同差异图的优点，并且对融合后高频子带系数进行了自适应阈值滤波，使得融合后的差异图既接近实际变化情况又具备很强的抗噪能力，检测精度优于现有的5种对比方法。

Claims (9)

1.基于自适应权值和高频阈值的SAR图像变化检测方法，包括如下步骤：

(1)读入不同时刻从同一地区获取的两幅SAR图像I₁和I₂,并计算这两幅图像对应像素的均值比差异图D₀；

(2)对均值比差异图D₀进行归一化，得到归一化后的均值比差异图D₁；

(3)将归一化后的均比值差异图D₁进行3*3窗口的中值滤波，抑制均值比图中的孤立像素点，得到中值滤波后的均值比差异图D_m；

(4)计算两幅输入图像I₁和I₂对应像素的邻域对数比差异图D₂；并对该邻域对数比差异图D₂进行归一化，得到归一化后的邻域对数比差异图D_l；

(5)根据均值比差异图D_m和邻域对数比差异图D_l，获得融合后的差异图像X_d：

5a)对均值比差异图D_m和邻域对数比差异图D_l分别进行3层平稳小波变换，得到每一幅差异图的1个低频子带系数和3个高频子带系数；

5b)对2幅差异图的低频子带系数按照最大-最小自适应权值策略进行融合，生成1个融合后的低频子带系数；

5c)对2幅差异图的高频子带系数按照局部能量最小原则进行融合，生成3个融合后的高频子带系数；

5d)对融合后的3个高频子带系数进行自适应阈值滤波操作；

5e)对融合后的1个低频子带系数和3个经自适应阈值滤波后的高频子带系数进行平稳小波逆变换，得到融合后的差异图像X_d；

(6)用k-means聚类方法将融合后的差异图X_d聚类成不同的两类，分别计算这两个不同类别的均值，定义均值较大的那一类为变化类，均值较小的那一类为非变化类，得到最终的变化检测结果。

2.根据权利要求书1中所述的方法，其中步骤(1)中两幅SAR图像I₁和I₂对应像素的均值比差异图D₀，按如下步骤进行：

1a)对第一幅输入图像I₁在像素(i,j)处的灰度值I₁(i,j)与第二幅输入图像I₂在像素(i,j)处的灰度值I₂(i,j)进行均值比运算，得到均值比矩阵D₀在像素(i,j)处的值d₀(i,j)：

上式中μ₁(i,j)和μ₂(i,j)分别表示第一幅输入图像I₁和第二幅输入图像I₂中以像素(i,j)为中心的3*3邻域窗口内所有像素灰度值的平均值；

1b)遍历两幅图像所有对应的像素，得到均值比差异图D₀＝{d₀(i,j)}。

3.根据权利要求1所述的方法，其中步骤(2)中对均值比差异图D₀进行归一化,按如下步骤进行：

2a)对均值比差异图D₀中每一像素按照如下公式进行归一化：

其中，d₁(i,j)表示均值比差异图D₀中的像素d₀(i,j)被归一化后的值；

2b)遍历均值比差异图D₀中的所有像素，得到归一化后的均值比差异图D₁＝{d₁(i,j)}。

4.根据权利要求1所述的方法，其中步骤(4)中两幅SAR图像I₁和I₂对应像素的邻域对数比差异图D₂，按如下步骤进行：

4a)对第一幅输入图像I₁中以像素(i,j)为中心的3*3局部邻域l₁(m,n)和第二幅输入图像I₂中以像素(i,j)为中心的3*3局部邻域l₂(m,n)分别进行大小为3，标准偏差为5的高斯滤波操作,得到滤波后的局部邻域l′₁(m,n)和l'₂(m,n)；

4b)按下式计算局部邻域l₁'(m,n)和l'₂(m,n)内对应像素灰度值的商的对数l′_log(m,n)：

4c)遍历两个局部邻域内的所有对应像素，得到局部邻域对数比矩阵Z_log(m,n)＝{l′_log(m,n)}；

4d)按下式取局部邻域对数比值矩阵Z_log(m,n)的均值d₂(i,j)作为邻域对数比差异图D₂在像素(i,j)处的值：

其中B是以像素(i,j)为中心的3*3局部邻域窗口的大小，N是局部邻域窗口内的元素数量；

4e)遍历两幅输入图像所有对应的像素，得到邻域对数比差异图D₂＝{d₂(i,j)}。

5.根据权利要求1所述的方法，其中步骤(4)中对邻域对数比差异图D₂进行归一化，按如下步骤进行：

4f)对邻域对数比差异图D₂中每一像素按照如下公式进行归一化：

其中，d_l(i,j)表示邻域对数比差异图D₂中的像素d₂(i,j)被归一化后的值；

4g)遍历邻域对数比差异图D₂中所有像素，得到归一化后的邻域对数比差异图D_l＝{d_l(i,j)}。

6.根据权利要求1所述的方法，其中步骤(5b)中对2幅差异图的低频子带系数按照最大-最小自适应权值策略进行融合,按如下公式进行：

上式中，和分别表示均值比差异图D_m和邻域对数比差异图D_l经3层平稳小波分解后的低频子带系数，表示融合后的低频系数。k表示低频子带系数的最大-最小自适应融合权值，其计算按下式进行：

7.根据权利要求1所述的方法，其中步骤(5c)中对2幅差异图的高频子带系数按照局部能量最小原则进行融合，按如下公式进行：

上式中，n＝(1,2,3)表示分解的层数；d＝(HH,HV,HD)分别表示差异图经3层平稳小波分解后高频子带系数的方向，其中HH表示水平方向、HV表示垂直方向、HD表示对角方向；表示融合后的高频子带系数在像素(i,j)处的值；和分别表示均值比差异图和邻域对数比差异图经3层平稳小波分解后的高频子带系数在像素(i,j)处的值；表示均值比差异图的高频子带系数中以像素(i,j)为中心的局部能量系数值，按如下公式计算：

上式中N_i,j表示以像素(i,j)为中心的3*3局部邻域窗口的大小，表示均值比差异图的高频子带系数在邻域窗口N_i,j内第S个系数的值；

表示在邻域对数比差异图的高频子带系数中以像素(i,j)为中心的局部能量系数值，按如下公式计算：

上式中N_i,j表示以像素(i,j)为中心的3*3局部邻域窗口的大小，表示邻域对数比差异图的高频子带系数在邻域窗口N_i,j内第S个系数的值。

8.根据权利要求1中所述的方法，其中步骤(5d)中对融合后的3个高频子带系数进行自适应阈值滤波操作，按下式进行：

其中，T表示自适应阈值参数，其值按下式计算：

其中λ是介于1与2之间的常数值；C是噪声能量分布标准差的调和系数，取值为3；σ是融合后高频子带系数的噪声估计值，其值按下式计算：

A(i,j)表示不同方向高频子带系数在局部邻域内的绝对值之和的均值，其值按下式计算：

其中B是以像素(i,j)为中心的3*3局部邻域窗口的大小，N是局部邻域窗口内的元素数量， 表示融合后的高频子带系数在邻域窗口B内第S个系数的值。

9.根据权利要求1中所述的方法，其中步骤(6)中用k-means聚类方法将融合差异图X_d聚成不同的两类，按如下步骤进行：

6a)从融合差异图X_d中随机选择2个像素点作为原始聚类中心；

6b)分别计算融合差异图X_d中每一像素点(i,j)与两个原始聚类中心的距离，并根据最小距离标准重新对相应像素点进行划分，把(i,j)归到与两个原始聚类中心距离最近的那一类中去；

6c)计算归类后将每一类的均值作为新的聚类中心；

6d)重复步骤6b)-6c)，直到每个聚类中心不再发生变化为止，得到两个类别的均值；

6e)将两个类别中均值小的那一类像素转化为0，即为未变化类，均值大的那一类像素转化为255，即为变化类，得到最终的变化检测结果。