CN102163333A - 谱聚类的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种谱聚类的SAR图像变化检测方法，它属于图像处理技术领域。其实现步骤为：(1)选择同一场景中大小相同，时段不同的两幅SAR图像T₁和T₂作为测试图像；(2)使用差值方法构造两幅测试图像T₁和T₂的差分子图D₁；(3)使用领域比值法构造两幅测试图像T₁和T₂的差分子图D₂；(4)使用乘积变换融合法对差分子图D₁和差分子图D₂进行融合，生成差分图D；(5)使用谱聚类算法，将差分图D的像素聚类，得到变化检测结果图。本发明很好的抑制差分图的斑点噪声，有效去除变化检测结果图中的杂点，得到的变化检测结果图比较准确，可用于遥感图像变化的检测处理。

Description

谱聚类的SAR图像变化检测方法

技术领域

本发明属于图像处理领域，涉及一种谱聚类的变化检测方法。该方法可用于从不同时间获取的同一地理区域的多时相遥感影像中，定性或定量地分析和确定地表变化特征和过程，得到我们所需要的地物变化信息的问题。

背景技术

现代遥感技术的飞速发展为变化检测提供了一种便捷的途径，遥感数据成为变化检测的主要数据源，这就对现有的遥感图像变化检测技术提出了新的挑战。合成孔径雷达SAR是遥感检测的一种雷达，他具有全天时检测地貌的优点，而不受天气影响。近年逐出现了许多种多时相SAR变化检测方法，最常用的SAR图像变化检测技术由两步组成：首先对多时相SAR图像构造差分图，然后在差分图基础上提取变化区域。其中：

对多时相SAR图像构造差分图是变化检测的一个关键问题，目前构造差分图的传统方法包括影像差值法、影像比值法和变化向量分析法等，这些方法的不足是：没有考虑多时相图像之间的相关性，再加上SAR图像受斑点噪声影响，因此不能很好的抑制斑点噪声，从而降低了在差分图的基础上提取变化区域的准确率。

在差分图的基础上提取变化区域是变化检测的另一个关键问题，传统的在差分图的基础上提取变化区域的算法是在差分图上确定变化阈值来提取差分图的变化区域。然而变化阈值如何确定一直没有好的方法解决，现有在差分图基础上确定阈值来提取变化区域的算法，包括Nakagami-Ratio GKIT，Log-Normal GKIT和Weibull-Ratio GKIT；这三种方法分别利用Nakagami-Gamma，Weibull和log-normal三种模型来选取最优变化阈值，理论上使得变化检测的总体错误率达到最低；这三种方法的缺陷在于都是假定差分图近似服从高斯分布，但是由于实际地物光谱特征的复杂性，用单一的高斯分布模型，会产生较大的偏差，得到错误的变化检测结果。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种谱聚类的变化检测方法，以抑制斑点噪声，提高变化检测的准确率。

本发明的技术方案是先选择一组常用于SAR图像变化检测精度评价的SAR图像数据作为实验数据，具体步骤包括如下：

(1)选择同一场景中大小相同，时段不同的两幅SAR图像T₁和T₂作为测试图像；

(2)使用差值方法构造两幅测试图像T₁和T₂的差分子图D₁：

D₁＝255-|T₁-T₂|，

(3)使用领域比值法构造两幅测试图像T₁和T₂的差分子图D₂：

3a)分别取得两幅测试图像在同一位置x上的像素值邻域集合S₁(x)和S₂(x)，其大小为M×M，M∈{3，5}；

3b)比较两个邻域集合S₁(x)和S₂(x)的相似性，得到差分子图D₂在位置x上的像素灰度值D₂(x)：

$D_2\left(x\right)=255\×\frac{\underset{i=1}{\overset{M\×M}{\mathrm{\Σ}}}\min \{S_1\left(x_i\right),S_2\left(x_i\right)\}}{\underset{i=1}{\overset{M\×M}{\mathrm{\Σ}}}\max \{S_1\left(x_i\right),S_2\left(x_i\right)\}}$

其中，S₁(x_i)和S₂(x_i)分别表示测试图像T₁和T₂在位置x上邻域集合S₁(x)和S₂(x)的第i个元素，其中0≤i≤M×M；

3c)对测试图像T₁和T₂的每一个像素点从左到右，从上到下重复步骤3a)和3b)，得到差分子图D₂

(4)使用乘积变换融合法对差分子图D₁和差分子图D₂进行融合，生成差分图D；

(5)使用谱聚类算法，将差分图D的像素聚类，得到变化检测结果图。

本发明与现有的技术相比具有如下优点：

1、本发明由于采用邻域比值法构造了测试图像T₁和T₂的差分子图D₂，可以最大程度地抑制背景，增强变化区域，有效减少了处于变化类和非变化类之间的区域对变化检测结果的影响。

2、本发明由于采用乘积变换融合法，融合了差分子图D₁和D₂，得到了差分图像D，同时由于在融合时充分考虑了差分图的像素之间的相似性和邻域信息，从而很好的抑制差分图的斑点噪声。

3、仿真结果表明，本发明采用的乘积变换融合法的差分像构造方法和谱聚类方法与比值算子，对数比值算子构造差分图方法和Nakagami-Ratio GKIT，Log-Normal GKIT和Weibull-Ratio GKIT方法相比，正确检测率高，错检和漏检个数低。

附图说明

图1是本发明的主要流程框图；

图2是本发明仿真使用的第一组SAR图像；

图3是用本发明与现有方法对图2构造的差分图像；

图4是对图3的差分图分别用谱聚类算法进行聚类的仿真结果图；

图5是本发明仿真使用的第二组SAR图像；

图6是用谱聚类方法与Nakagami-Ratio GKIT，Log-Normal GKIT和Weibull-Ratio GKIT三种方法对图5提取变化区域得到的变化检测结果图。

具体实施方式

参照图1，本发明的SAR图像变化检测方法包括如下步骤：

步骤1.选择在同一雷达拍摄下的同一地区场景中大小相同，时段不同的两幅SAR图像T₁和T₂作为测试图像。

步骤2.使用差值方法构造两幅测试图像T₁和T₂的差分子图D₁：

D₁＝255-|T₁-T₂|，

如果测试图像T₁和T₂在某一个位置上的像素值没有发生变化，则差分子图D₁在该位置上的像素值为255，如果测试图像T₁和T₂在某一个位置上的像素值发生的变化程度越大，则差分子图D₁在该位置像素值越接近0。

步骤3.使用领域比值法构造两幅测试图像T₁和T₂的差分子图D₂：

3a)分别取得两幅测试图像T₁和T₂在同一位置x上的像素点邻域集合S₁(x)和S₂(x)，其大小均为M×M，M∈{3，5}，其中邻域集合S₁(x)以测试图像T₁的位置x为中心取一个M×M大小的窗口，邻域集合S₂(x)以测试图像T₂的位置x为中心取一个M×M大小的窗口；

3b)比较两个邻域集合S₁(x)和S₂(x)的相似性，得到差分子图D₂在位置x上的像素值D₂(x)：

$D_2\left(x\right)=255\×\frac{\underset{i=1}{\overset{M\×M}{\mathrm{\Σ}}}\min \{S_1\left(x_i\right),S_2\left(x_i\right)\}}{\underset{i=1}{\overset{M\×M}{\mathrm{\Σ}}}\max \{S_1\left(x_i\right),S_2\left(x_i\right)\}}$

其中，S₁(x_i)和S₂(x_i)分别表示测试图像T₁和T₂在位置x上邻域集合S₁(x)和S₂(x)的第i个元素，其中0≤i≤M×M；D₂(x)值越小，说明测试图像T₁和T₂在位置x上的差别越大，则位置x上的像素点属于变化区域的可能性也越大，反之，该位置像素点属于非变化区域的可能性越大；

3c)对测试图像T₁和T₂的每一个像素点从左到右，从上到下重复步骤3a)和3b)，得到差异子图D₂。

步骤4.使用乘积变换融合法对差分子图D₁和差分子图D₂进行融合，生成差分图D：

4a)对大小为H×W的差分子图D₁分别取第i个像素点邻域集合N_i和第j个像素点邻域集合N_j，1≤i≤H×W，1≤j≤H×W；其中邻域集合N_i以位置i为中心取一个P×P大小的窗口，邻域集合N_j以位置j为中心取一个P×P大小的窗口，P∈{3，5}；

4b)比较差分子图D₁两个邻域集合N_i和N_j的相似性，得到差分子图D₁在第i个像素点和第j个像素点的相似性u(i，j)：

$u(i,j)=\frac{1}{Z\left(i\right)}{e}^{-\frac{{\left|\right|d_1\left(N_i\right)-d_1\left(N_j\right)\left|\right|}_{2,a}^2}{h^2}}$

其中

d₁(N_i)表示在差分子图D₁(i)中，以i为中心的窗口，N_i内所有像素点的灰度值依次构成的向量，所有窗口取定同一方向；

表示带权的欧拉距离，a＞0，它是高斯核的标准方差，h是一个控制衰减程度的参数，h≠0；u(i，j)满足：0≤u(i，j)≤1，∑_ju(i，j)＝1，两个像素点i和j之间的相似性由像素值向量d(N_i)与d(N_j)之间的相似性决定，通过高斯加权的欧氏距离来衡量，相似性越高，u(i，j)的值越大。

4c)利用如下公式得到差分图D在位置i上的像素值D(i)：

D(i)＝∑u(i，j)d₂(j)

其中d₂(j)表示差分子图D₂的第j个元素的值；

4d)对两个差分子图D₁和D₂的每一个像素点从左到右，从上到下重复步骤4a)-4c)，得到差分图D。

步骤5.使用谱聚类算法，将差分图D的像素点聚类，得到变化检测结果图：

5a)将大小为H×W差分图D的每一个像素点看成是一个数据集合X上一个点，X＝(x₁，x₂，...，x_n)∈R^n×1，n＝H×W，X上每一个点与差分图D的一个像素点对应；

5b)构造数据集合X的相似度矩阵：G＝(g_ij)∈R^n×n，当i≠j时，g_ij＝exp(-|x_i-x_j|²/σ²)，当i＝j时，g_ij＝0，1≤i≤n，1≤j≤n，其中x_i是数据集合X的第i个点，x_j是数据集合X第j个点，σ是一个控制衰减程度的参数，σ≠0；

5c)根据相似度矩阵G，构造拉布拉斯矩阵：L＝A^-1/2GA^-1/2，其中A是一个对角矩阵：A＝(a_ij)∈R^n×n，当i≠j时，a_ij＝0，当i＝j时，

1≤i≤n，1≤j≤n；

5d)根据拉布拉斯矩阵L，构造特征矩阵V＝[v₁，v₂]∈R^N×2，其中v₁，v₂是拉布拉斯矩阵L的两个最大特征值所对应的特征向量，v₁，v₂都为列向量；

5e)根据特征矩阵V，构造标准矩阵Y：

①利用如下公式得到标准矩阵Y的第i行第j列上的数据Y_ij：

$Y_{\mathrm{ij}}=V_{\mathrm{ij}}/{\left(\underset{j}{\mathrm{\Σ}}{V}_{\mathrm{ij}}^2\right)}^{1/2};$

其中V_ij表示特征矩阵V的第i行第j列上的数据；

②对特征矩阵V的每一个数据点从左到右，从上到下重复步骤①，得到标准矩阵Y。

5f)将标准矩阵Y的每一行看成是R²空间上一点，使用k均值将标准矩阵Y聚为2类；如果标准矩阵Y的第i行数据为第b类，则将原数据点x_i也划为第b类，其中b∈{0，1}，得到变化检测结果图。

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

1、仿真条件：

在CPU为core2 2.4HZ、内存2G、WINDOWS XP***上使用MATLAB7.0进行了仿真。

本发明仿真使用的第一组SAR测试图像，是图2中的Ottaw组图像，其中图2(a)和图2(b)的拍摄时间分别为1997.05和1997.08，大小均为290×350，图2(c)为参考图。

本发明仿真使用的第二组SAR测试图像，是图5中的Berma组图像，Berma组图是ERS-2卫星拍摄的两幅大小均为301×301的SAR图像，图5错误！未找到引用源。(a)和图5(b)的拍摄时间分别为1999.4和1999.5，图5(c)为参考图。

2、仿真参数

对于具有参考图的实验仿真图，可进行定量的变化检测结果分析，评价变化检测结果的指标有4个：

①计算漏检个数：统计实验结果图中发生变化区域的像素个数，与参考图中变化区域的像素个数进行对比，把参考图中发生变化但实验结果图中检测为未变化的像素个数，称为漏检个数FN；

②计算错检个数：统计实验结果图中未发生变化区域的像素个数，与参考图中非变化区域的像素个数进行对比，把参考图中未发生变化但实验结果图中检测为变化的像素个数，称为错检个数FP；

③总错误数：漏检个数FN与错检个数FP之和为总错误数TE；

④计算实验检测结果图的正确检测率：

正确检测率PCC的定义为：正确测得的目标数与实验数据图的总像素个数之比，如公式1)所示；其中正确测得的目标数等于参考图中发生变化且实验结果图中检测为变化的像素个数TP与参考图中未发生变化且实验结果图中检测为未变化的像素个数TN的和；

PCC＝((TP+TN)/(TP+FP+TN+FN))              1)

这四个指标集中反映了对变化区域的检测性能，具有很强的针对性。

3、仿真内容

3.1)用本发明方法与现有比值算子、对数比值算子方法，对图2进行差分图像仿真构造，仿真结果如图3所示，其中图3(a)为用比值算子构造的差分图像，图3(b)为用对数比值算子构造的差分图像，图3(c)为用本发明方法构造的差分图像。

3.2)对图3的差分图分别用谱聚类算法进行聚类，其聚类结果如图4所示，其中图4(a)是对图3(a)聚类后的仿真结果图，图4(b)是对图3(b)聚类后的仿真结果图，图4(c)是对图3(c)聚类后的仿真结果图。

3.3)用谱聚类算法、Nakagami-Ratio GKIT算法、Log-Normal GKIT算法和Weibull-Ratio GKIT算法分别对图5所述的第二组SAR图像进行变化检测，变化检测结果图如图6所示，其中图6(a)是用Log-Normal GKIT方法的仿真结果图，图6(b)是用Nakagami-Ratio GKIT方法的仿真结果图，图6(c)是用Weibull-Ratio GKIT方法的仿真结果图，图6(d)是用谱聚类方法的仿真结果图。

4、仿真实验结果及分析

从图3可以看出，本发明方法的仿真结果图的斑点噪声较少，对斑点噪声有很好的抑制；

从图4可以看出，本发明方法的仿真结果图杂点较少，变化检测比较好。对图4变化检测结果定量分析如表1所示。

表1对图4变化检测结果定量分析

从表1中可以看出，本发明方法在变化检测总错误数上，比其他对比实验方法都要少，而且在误检数与漏检数上相对都比较均衡，同时与图(2c)参考图像对比可看出，本发明对斑点噪声的抑制效果比其他两个方法都要好。

从图6可以看到，本发明方法的仿真结果图杂点较少，变化检测比较好，并且对细小边缘部分的检测结果非常好，能有效的检测出细小边缘。对图6的变化检测结果定量分析如表2所示。

表2对图6的变化检测结果定量分析

由表2可以看出，本发明方法，漏检个数比较少，变化检测总错误数比其他方法都要少，本发明方法比其他三种方法要更优越。

Claims (3)

1.一种谱聚类的SAR图像变化检测方法，包括如下步骤：

(1)选择同一场景中大小相同，时段不同的两幅SAR图像T₁和T₂作为测试图像；

(2)使用差值方法构造两幅测试图像T₁和T₂的差分子图D₁：

D₁＝255-|T₁-T₂|，

(3)使用领域比值法构造两幅测试图像T₁和T₂的差分子图D₂：

3a)分别取得两幅测试图像在同一位置x上的像素值邻域集合S₁(x)和S₂(x)，其大小为M×M，M∈{3，5}；

3b)比较两个邻域集合S₁(x)和S₂(x)的相似性，得到差分子图D₂在位置x上的像素灰度值D₂(x)：

$D_2\left(x\right)=255\×\frac{\underset{i=1}{\overset{M\×M}{\mathrm{\Σ}}}\min \{S_1\left(x_i\right),S_2\left(x_i\right)\}}{\underset{i=1}{\overset{M\×M}{\mathrm{\Σ}}}\max \{S_1\left(x_i\right),S_2\left(x_i\right)\}}$

其中，S₁(x_i)和S₂(x_i)分别表示测试图像T₁和T₂在位置x上邻域集合S₁(x)和S₂(x)的第i个元素，其中0≤i≤M×M；

3c)对测试图像T₁和T₂的每一个像素点从左到右，从上到下重复步骤3a)和3b)，得到差分子图D₂

(4)使用乘积变换融合法对差分子图D₁和差分子图D₂进行融合，生成差分图D；

(5)使用谱聚类算法，将差分图D的像素聚类，得到变化检测结果图。

2.根据权利要求1所述的SAR图像变化检测方法，其中步骤(4)所述的使用乘积变换融合法对差分子图D₁和差分子图D₂进行融合，生成差分图D，按如下步骤进行：

4a)对大小为H×W的差分子图D₁分别取第i个像素值邻域集合N_i和第j个像素值邻域集合N_j，其大小为P×P，P∈{3，5}，1≤i≤H×W，1≤j≤H×W；

4b)比较差分子图D₁两个邻域集合N_i和N_j的相似性，得到差分子图D₁在第i个像素和第j个像素的相似性u(i，j)：

$u(i,j)=\frac{1}{Z\left(i\right)}{e}^{-\frac{{\left|\right|d_1\left(N_i\right)-d_1\left(N_j\right)\left|\right|}_{2,a}^2}{h^2}}$

其中d₁(N_i)表示在差分子图D₁(i)中，以i为中心的窗口，N_i内所有像素点的灰度值依次构成的向量，所有窗口取定同一方向，

表示带权的欧拉距离，a＞0，它是高斯核的标准方差，h是一个控制衰减程度的参数，h≠0；

4c)利用如下公式得到差分图D在位置i上的像素灰度值D(i)：

D(i)＝∑u(i，j)d₂(j)

其中d₂(j)表示差分子图D₂的第j个元素的值；

4d)对两个差分子图D₁和D₂的每一个像素点从左到右，从上到下重复步骤4a)-4c)，得到差分图D。

3.根据并权利要求1所述的SAR图像变化检测方法，其中步骤5所述的使用谱聚类算法，将差分图D的像素聚类，按如下步骤进行：

5a)将大小为H×W差分图D的每一个像素看成是一个数据集合X上一个点，X＝(x₁，x₂，...，x_n)∈R^n×1，n＝H×W，其中X上每一个点与差分图D的一个像素点对应；

5b)对X构造相似度矩阵G＝(g_ij)∈R^n×n，当i≠j时，g_ij＝exp(-|x_i～x_j|²/σ²)，当i＝j时，g_ij＝0，1≤i≤n×n，1≤j≤n×n，其中x_i是数据集合X的第i个点，x_j是数据集合X第j个点，σ是一个控制衰减程度的参数，σ≠0；

5c)根据相似度矩阵G，构造拉布拉斯矩阵：L＝A^-1/2GA^-1/2，其中A是一个对角矩阵：A＝(a_ij)∈R^n×n，当i≠j  时，a_ij＝0，当i＝j时，

1≤i≤n×n，1≤j≤n×n；

5d)根据拉布拉斯矩阵L，构造特征矩阵V＝[v₁，v₂]∈R^n×2，其中v₁，v₂是拉布拉斯L的两个最大特征值所对应的特征向量，v₁，v₂都为列向量；

5e)根据特征矩阵V，构造标准矩阵Y：

①利用如下公式得到标准矩阵Y的第i行第j列上的数据Y_ij：

$Y_{\mathrm{ij}}=V_{\mathrm{ij}}/{\left(\underset{j}{\mathrm{\Σ}}{V}_{\mathrm{ij}}^2\right)}^{1/2};$

其中V_ij表示特征矩阵V的第i行第j列上的数据；

②对特征矩阵V的每一个数据点从左到右，从上到下重复步骤①，得到标准矩阵Y。

5f)将标准矩阵Y的每一行看成是R²空间上一点，对标准矩阵Y使用k均值将它聚为2类；如果标准矩阵Y的第i行数据为第b类，则将原数据点x_i也划为到第b类，其中b∈{0，1}，得到变化检测结果图。

Images