CN101923711B - 基于邻域相似性及掩模增强的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于邻域相似性及掩模增强的SAR图像变化检测方法，它涉及遥感图像处理领域，主要解决SAR图像变化检测的差异图构造受SAR图像斑点噪声影响严重的问题，其实现步骤为：(1)对不同时相同一场景的两幅SAR图像I₁和I₂进行几何校正和配准；(2)用邻域相似度算子，构造图像I₁和I₂的差异影像图DI；(3)对差异影像图DI进行掩模增强处理，得到新的差异影像图NDI；(4)利用初始点确定的K均值聚类方法，对新的差异影像图NDI进行聚类，得到变化检测结果图CDI。本发明具有运行效率高，算法复杂度和时间复杂度低，实验效果好的特点，可用于两时相SAR图像的变化检测。

Description

基于邻域相似性及掩模增强的SAR图像变化检测方法

技术领域

本发明属于遥感图像处理领域，涉及遥感图像的变化检测，可用于对SAR图像变化的检测。 

背景技术

合成孔径雷达SAR是现代遥感领域的一项突破性技术，成为获取信息的一个重要手段。随着SAR技术的不断发展和成熟，及SAR系统的分辨率的不断提高，使得SAR系统具有了全天候、全天时及覆盖面积大等优点。不同时相同一场景的SAR图像的获得，使得通过SAR图像变化检测技术提供目标与场景实时动态信息成为可能。 

变化检测是遥感技术的主要应用之一，它通过对不同时相同一场景的各种图像进行对比分析，根据图像之间的差异来得到人们所需要的地物或目标的变化信息。变化检测技术可以检测出不同时期图像灰度值或局部纹理之间的变化，在此基础上获得感兴趣目标在形状、位置、数量、及其它属性的变化情况。这些变化可能是由图像场景的真实变化引起的，也可能是由入射角、大气条件、传感器精度、地面湿度等变化引起的。SAR图像的变化检测则特指研究不同时相同一场景的两幅或多幅SAR图像之间的变化信息，SAR图像变化检测技术在环境检测等已被广泛应用，如土地利用分析、森林采伐检测、灾情评估等。 

然而，SAR图像在具有上述优点的同时，受斑点噪声影响是其应用的一个主要瓶颈，如何从SAR图像自身特点出发，有效抑制斑点噪声是SAR图像实际应用中必须考虑的一个问题。 

随着计算机和成像技术的不断发展，SAR图像变化检测技术也在不断提高和改进中，近年来逐渐形成了以下几种较为经典的变化检测算法： 

图像差值法：其主要过程是将两时相SAR图像对应像素点的灰度值相减得到差异图，然后选取0～255的阈值进行分割，得到变化和非变化区域。图像差值法算法简单，方便易行，但这种方法的缺点也比较明显：容易受SAR成像质量、波谱特征不同等客观条件影响，易产生“伪变化”信息。图像比值法：其主要过程是计算两时相SAR图像对应像素点灰度值的比值得到差异图，如果一个像素没有发生变化，则比值 应该接近与1，反之，远大于或远小于1。该方法对SAR图像的乘性噪声是不敏感的，但得到的差异图精确度往往会受变化类型影响较大。分类后比较法：其主要过程是先用同一种分类方法对每一时相SAR图像进行分类，然后将分类结果图进行比较，以确定变化的类别和区域。如果对应像素的分类类别相同，则认为该像素没有发生变化，反之，发生了变化。该方法不仅可以检测变化区域，还可以确定变化类别信息，同时受配准误差影响较小，但它的缺点是受分类器的误差积累影响很大。 

以上三种方法是最为经典的变化检测方法，另外还有变化向量分析法，主分量分析法等等。但是这些方法均受SAR图像斑点噪声影响较大，必须预先对两时相SAR图像进行滤波去噪，才能得到较好的变化检测结果，而滤波在变化检测方法中需要占用较长的运行时间，增加了变化检测方法的复杂度，另外这些方法构造的差异图中变化类与非变化类区别都不够明显，增加了后续的分类的难度。 

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于邻域相似性及掩模增强的SAR图像变化检测方法，以在无需进行滤波去噪的情况下，抑制SAR图像的部分斑点噪声，提高变化检测结果的正确率。 

本发明的技术方案是：基于邻域相似度的原理，构建了一个由比值算子演化而来的邻域相似度NSR算子，根据NSR算子构造两时相SAR图像的差异影像图DI，再通过对DI的掩模增强处理后K均值聚类得到变化结果图。其具体实现步骤如下： 

(1)对不同时相同一场景的两幅SAR图像I₁和I₂进行几何校正和配准； 

(2)根据邻域相似度NSR算子，按照如下步骤构造图像I₁和I₂的差异影像图DI； 

2a)分别取得两时相SAR图像I₁和I₂在同一位置x上的像素灰度值邻域集合N₁(x)和N₂(x)，其大小为N×N，N∈{3，5，7，9}； 

2b)比较两个邻域集合N₁(x)和N₂(x)的相似性，得到差异图DI在位置x上的像素灰度值DI(x)： 

$\mathrm{DI}\left(x\right)=\frac{\underset{i=1}{\overset{N\×N}{\mathrm{\Σ}}}\min (N_1\left(x_i\right),N_2\left(x_i\right))}{\underset{i=1}{\overset{N\×N}{\mathrm{\Σ}}}\max (N_1\left(x_i\right),N_2\left(x_i\right))},$

其中，N₁(x_i)和N₂(x_i)分别表示图像I₁和I₂在位置x上邻域集合N₁(x)和N₂(x)的第i个元素，集合N₁(x)与N₂(x)越相似，则DI(x)值越大，图像I₁和I₂在位置x上的像 素点属于非变化区域的可能性越大，反之，属于变化区域的可能性越大； 

2c)对图像I₁和I₂的每一位置x从左到右，从上到下重复步骤2a)和2b)，得到差异影像图DI； 

(3)对差异影像图DI进行掩模增强处理，得到新的差异影像图NDI， 

所述的掩模增强处理包括如下步骤： 

3a)根据差异影像图DI的图像大小选用掩模，若DI小于1000×1000，则选用中心系数为-8的掩模，若DI大于1000×1000，则选用中心系数为-9的掩模； 

3b)利用所选的掩模，增大DI在位置x上像素点与其周围相邻像素点间的差异，得到新的差异影像图在位置x上的值NDI(x)： 

$\mathrm{NDI}\left(x\right)=\underset{i=1}{\overset{\mathrm{LN}}{\mathrm{\Σ}}}(\mathrm{DI}\left(x_i\right)-L\left(x_i\right)),1\≤i\≤\mathrm{LN},$

其中LN表示掩模的长度，取值为9，L(x_i)为掩模中的第i个数的值，DI(x_i)表示DI中大小为3×3的x邻域集合中的第i个像素的灰度值； 

3c)对DI的从左到右，从上到下重复步骤3b)，得到新的差异影像图NDI； 

(4)利用初始点确定的K均值方法，对新的差异影像图NDI进行聚类，得到变化检测结果图CDI。 

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

1、本发明采用邻域相似度NSR算子构造的差异影像图，抑制了SAR图像的部分斑点噪声，不需要对原始SAR图像滤波去噪，降低了变化检测方法的复杂度； 

2、本发明针对邻域相似度算子构造的差异影像图，创建的掩模增强处理方法，进一步减弱了差异影像图中的噪声影响，同时有效地增大了该差异影像图中的变化类中心和非变化类中心的距离，有利于对差异影像图的分类，另外，该掩模增强处理方法运行时间少，此优点在检测两时相大图时尤为重要； 

3、仿真结果表明，本发明采用的NSR算子构造差异影像图方法和掩模增强处理方法较比值R算子，对数比值LR算子构造差异影像图方法和双边滤波处理方法的正确检测率高，错检和漏检个数低，时间复杂度小。 

附图说明

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

图2是本发明设计的掩模； 

图3是第一组实验仿真图； 

图4是不同算子构造的差异影像图； 

图5是对应图4的K均值聚类结果图； 

图6是掩模增强处理前后的差异影像图； 

图7是对应图6的K均值聚类结果图； 

图8是对应图6的数据分布直方图； 

图9是第二组实验仿真图； 

图10是双边滤波处理DI后的K均值聚类结果图； 

图11是掩模增强处理DI后的K均值聚类结果图； 

图12是双边滤波处理和掩模增强处理的运行时间对比图； 

图13是分别将图9中所标记的A，B和C区域放大后的原图； 

图14是对应图13的双边滤波处理DI后的K均值聚类结果图； 

图15是对应图13的掩模增强处理DI后的K均值聚类结果图。 

具体实施方式

参照附图1，本发明包括如下步骤： 

步骤一：通过遥感图像处理软件ERDAS对不同时相同一场景的两幅SAR图像I₁和I₂进行几何校正和配准。 

步骤二：根据邻域相似度NSR算子，按照如下步骤构造图像I₁和I₂的差异影像图DI。 

2a)分别取得两时相SAR图像I₁和I₂在同一位置x上的像素灰度值邻域集合N₁(x)和N₂(x)，其大小为N×N，N∈{3，5，7，9}； 

2b)比较两个邻域集合N₁(x)和N₂(x)的相似性，得到差异图DI在位置x上的像素灰度值DI(x)： 

$\mathrm{DI}\left(x\right)=\frac{\underset{i=1}{\overset{N\×N}{\mathrm{\Σ}}}\min (N_1\left(x_i\right),N_2\left(x_i\right))}{\underset{i=1}{\overset{N\×N}{\mathrm{\Σ}}}\max (N_1\left(x_i\right),N_2\left(x_i\right))},$

其中，N₁(x_i)和N₂(x_i)分别表示图像I₁和I₂在位置x上邻域集合N₁(x)和N₂(x)的第i个元素，集合N₁(x)与N₂(x)越相似，则DI(x)值越大，图像I₁和I₂在位置x上的像素点属于非变化区域的可能性越大，反之，属于变化区域的可能性越大； 

2c)对图像I₁和I₂的每一位置x从左到右，从上到下重复步骤2a)和2b)，得到 差异影像图DI。 

步骤三：对差异影像图DI按如下步骤进行掩模增强处理，得到新的差异影像图NDI。 

3a)根据差异影像图DI的图像大小选用掩模，如附图2所示，其中图2(a)是中心系数为-8的掩模，适用于DI小于1000×1000的情况，图2(b)是中心系数为-9的掩模，适用于DI大于1000×1000的情况； 

3b)利用所选的掩模，增大DI在位置x上像素点与其周围相邻像素点间的差异，得到新的差异影像图在位置x上的值NDI(x)： 

$\mathrm{NDI}\left(x\right)=\underset{i=1}{\overset{\mathrm{LN}}{\mathrm{\Σ}}}(\mathrm{DI}\left(x_i\right)-L\left(x_i\right)),1\≤i\≤\mathrm{LN},$

其中LN表示掩模的长度，取值为9，L(x_i)为掩模中的第i个数的值，DI(x_i)表示DI中大小为3×3的x邻域集合中的第i个像素的灰度值； 

3c)对DI的从左到右，从上到下重复步骤3b)，得到新的差异影像图NDI。 

步骤四：利用初始点确定的K均值方法，对新的差异影像图NDI进行聚类，得到变化检测结果图CDI。 

4a)选用两个初始聚类中心为： 

c₁＝gray_min+(gray_max-gray_min)/4 

c₂＝gray_min-(gray_max-gray_min)/4， 

其中，gray_min和gray_max分别表示NDI中所有像素点灰度的最小值和最大值； 

4b)根据初始聚类中心，对NDI中所有的像素点进行首次分类，计算每个像素点到两个聚类中心c₁和c₂之间的欧氏距离，并将该像素点归到与它最近的那一类； 

4c)根据分类结果，计算新的聚类中心c_i， 

$c_i=\frac{1}{\left|N_i\right|}\underset{j=1}{\overset{N_i}{\mathrm{\Σ}}}{x}_j,i\∈\left\{\mathrm{1,2}\right\},1\≤j\≤N_i,$

其中c_i表示第i类的新的聚类中心，N_i表示第i类中像素点的总个数，x_j表示第i类中的第j个像素点； 

4d)根据新的聚类中心c_i，对NDI中所有的像素点再次分类，计算每个像素点到新的聚类中心c_i，i∈{1，2}之间的欧氏距离，并将该像素点归到与自己最近那一类； 

4e)重复步骤4c)和4d)直到聚类结束，聚类结束的停机准则如下： 

满足 $J=\underset{i=1}{\overset{K=2}{\mathrm{\&Sigma;}}}\left(\underset{j=1}{\overset{N_i}{\mathrm{\&Sigma;}}}{\left|\right|x_j-c_i\left|\right|}^2\right)$ 最小，且 $\underset{i=1}{\overset{K=2}{\mathrm{\&Sigma;}}}|c_i^{\&prime;}-c_i|<\mathrm{\&epsiv;},(\mathrm{\&epsiv;}=0.01),$

其中K表示聚类的类别数，取值为2，||x_j-c_i||表示x_j和c_i的欧式距离，c′_i表示下一次迭代后的聚类中心； 

4f)标记最终的聚类结果，第一类标记为0，第二类标记为1，并将标记为1的像素转换为255，得到变化检测结果图。 

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

1仿真参数 

对于具有参考图的第一组实验仿真图，进行定量的变化检测结果分析： 

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

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

总检测错误数OE等于漏检个数和错检个数的和； 

③计算正确检测率： 

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

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

2仿真内容 

①将邻域相似度NSR算子构造的差异影像图，与比值R算子，对数比值LR算子构造的差异影像图的效果对比，及对应的变化检测结果对比和定量分析，通过第一组实验仿真图完成； 

②对差异影像图的掩模增强处理MEP前后的直方图数据分布对比，及对应的变化检测结果对比和定量分析，通过第一组实验仿真图完成； 

③对差异影像图的双边滤波处理BF和掩模增强处理MEP的运行时间对比分析，及对应的变化检测效果对比，通过第二组实验仿真图完成。 

3仿真实验结果及分析 

①反应Ottawa地区水灾的RADARSAT SAR图像如图3所示，图3(a)和图3(b)的拍摄时间分别为1997.05和1997.08，大小均为290×350，图3(c)为参考图，其中真实变化目标数为16049。 

用不同算子构造的差异影像图如图4所示，其中图4(a)，图4(b)和图4(c)分别表示由R算子，LR算子和NSR算子构造的差异影像图。从图4可以看出本发明的NSR算子构造的差异影像图最为平滑，抑制了图像的部分噪声，差异影像图的变化类与非变化类之间的区别最为明显。用不同算子构造的差异影像图的K均值聚类结果如图5所示，其中图5(a)，图5(b)和图5(c)分别表示由R算子，LR算子和NSR算子构造的差异影像图的K均值聚类结果图，从图5可以看出，NSR算子构造的差异影像图的K均值聚类结果图杂点最少，变化检测结果最好。用不同算子构造的差异影像图的类内方差如表1所示，用不同算子构造差异影像图的变化检测结果分析如表2所示。 

表1不同算子构造的差异影像图的类内方差 

从表1可以看出NSR算子构造的差异影像图的类内方差最小，类内的数据分布最为均匀，有利于降低FP。 

表2不同算子构造差异影像图的变化检测结果分析 

从表2可以看出NSR算子构造差异影像图的K均值聚类结果的FN，FP和OE最小，PCC最高。 

②用NSR算子构造第一组仿真数据的差异影像图如图6所示，其中图6(a)是掩模增强处理前的差异影像图，图6(b)是掩模增强处理后的新差异影像图，它们的数据分布直方图如图8所示，实线和虚线分别代表的是掩模处理前后的差异影像图数据分布。从图8可以看出，虚线的两个波峰距离较远，表明掩模增强处理后差异影像图的类内中心距离增大了，有利于后续的分类。掩模增强处理差异影像图前后的K均 值聚类结果如图7所示，从图7可以看出，差异影像图经过掩模增强处理后的检测结果杂点有效减少，边缘更加平滑。掩模增强处理差异影像图前后的变化检测结果分析如表3所示。 

表3掩模处理差异影像图前后的变化检测结果分析 

从表3可以看出差异影像图经过掩模增强处理后，FN，FP和OE明显减少，PCC得到提高。 

③关于黄河入海口地域变化的图像如图9所示，由radarSAR-2卫星获得，图9(a)和图9(b)的拍摄时间分别为2008.06.18和2009.06.19，大小均为7692×7666。双边滤波处理和掩模增强处理的时间对比如图12所示，实线和虚线分别代表了双边滤波处理与掩模增强处理随着图像不断增大时的运行时间。 

从图12可以看出，在图像小于800×800时，掩模增强处理的时间与双边滤波处理时间相差不大，但是随着图像的不断增大，双边滤波处理时间几乎成指数级增长，而掩模增强处理时间却以非常缓慢的速率在增加，如对于第二组仿真图的双边滤波处理时间为1596.844s，而掩模增强处理时间仅为18.562s，且取得了与双边滤波一样的，甚至更好的检测效果，如图10和图11所示。放大它们的代表性区域A，B和C后的原图如图13所示，它们的变化检测结果图，如图14和15所示，其中A，B和C区域分别反应的是2008年到2009年间黄河入海口地域的陆地，港口和海岸的变化。 

Claims (3)

1.一种基于邻域相似性及掩模增强的SAR图像变化检测方法，包括如下步骤：

(1)对不同时相同一场景的两幅SAR图像I₁和I₂进行几何校正和配准；

(2)根据邻域相似度NSR算子，按照如下步骤构造图像I₁和I₂的差异影像图DI；

2a)分别取得两时相SAR图像I₁和I₂在同一位置x上的像素灰度值邻域集合N₁(x)和N₂(x)，其大小为N×N，N∈{3，5，7，9}；

2b)比较两个邻域集合N₁(x)和N₂(x)的相似性，得到差异图DI在位置x上的像素灰度值DI(x)：

$\mathrm{DI}\left(x\right)=\frac{\underset{i=1}{\overset{N\&times;N}{\mathrm{\&Sigma;}}}\min (N_1\left(x_i\right),N_2\left(x_i\right))}{\underset{i=1}{\overset{N\&times;N}{\mathrm{\&Sigma;}}}\max (N_1\left(x_i\right),N_2\left(x_i\right))},$

其中，N₁(x_i)和N₂(x_i)分别表示图像I₁和I₂在位置x上邻域集合N₁(x)和N₂(x)的第i个元素，集合N₁(x)与N₂(x)越相似，则DI(x)值越大，图像I₁和I₂在位置x上的像素点属于非变化区域的可能性越大，反之，属于变化区域的可能性越大；

2c)对图像I₁和I₂的每一位置x从左到右，从上到下重复步骤2a)和2b)，得到差异影像图DI；

(3)对差异影像图DI进行掩模增强处理，得到新的差异影像图NDI，

所述的掩模增强处理包括如下步骤：

3a)根据差异影像图DI的图像大小选用掩模，若DI小于1000×1000，则选用中心系数为-8的掩模，若DI大于1000×1000，则选用中心系数为-9的掩模；

3b)利用所选的掩模，增大DI在位置x上像素点与其周围相邻像素点间的差异，得到新的差异影像图在位置x上的值NDI(x)：

$\mathrm{NDI}\left(x\right)=\underset{i=1}{\overset{\mathrm{LN}}{\mathrm{\&Sigma;}}}(\mathrm{DI}\left(x_i\right)-L\left(x_i\right)),1\&le;i\&le;\mathrm{LN},$

其中LN表示掩模的长度，取值为9，L(x_i)为掩模中的第i个数的值，DI(x_i)表示DI中大小为3×3的x邻域集合中的第i个像素的灰度值；

3c)对DI的从左到右，从上到下重复步骤3b)，得到新的差异影像图NDI；

(4)利用初始点确定的K均值聚类方法，对新的差异影像图NDI进行聚类，得到变化检测结果图CDI。

2.根据权利要求1所述的SAR图像变化检测方法，其中步骤(1)所述的对不同时相同一场景的两幅SAR图像I₁和I₂进行几何校正和配准，是通过遥感图像处理软件ERDAS完成。

3.根据权利要求1所述的SAR图像变化检测方法，其中步骤(4)所述的利用初始点确定的K均值方法，对新的差异影像图NDI进行聚类，按如下步骤进行；

4a)选用两个初始聚类中心为：

c₁＝gray_min+(gray_max-gray_min)/4

c₂＝gray_min-(gray_max-gray_min)/4

其中，gray_min和gray_max分别表示NDI中所有像素点灰度的最小值和最大值；

4b)根据初始聚类中心，对NDI中所有的像素点进行首次分类，计算每个像素点到两个聚类中心c₁和c₂之间的欧氏距离，并将该像素点归到与它最近的那一类；

4c)根据分类结果，计算新的聚类中心c_i，

$c_i=\frac{1}{\left|N_i\right|}\underset{j=1}{\overset{N_i}{\mathrm{\&Sigma;}}}{x}_j,i\&Element;\left\{\mathrm{1,2}\right\},1\&le;j\&le;N_i,$

其中c_i表示第i类的新的聚类中心，N_i表示第i类中像素点的总个数，x_j表示第i类中的第j个像素点；

4d)根据新的聚类中心c_i，对NDI中所有的像素点再次分类，计算每个像素点到新的聚类中心c_i，i∈{1，2}之间的欧氏距离，并将该像素点归到与自己最近那一类；

4e)重复步骤4c)和4d)直到聚类结束，聚类结束的停机准则如下：

满足 $J=\underset{i=1}{\overset{K=2}{\mathrm{\&Sigma;}}}\left(\underset{j=1}{\overset{N_i}{\mathrm{\&Sigma;}}}{\left|\right|x_j-c_i\left|\right|}^2\right)$ 最小，且 $\underset{i=1}{\overset{K=2}{\mathrm{\&Sigma;}}}|c_i^{\&prime;}-c_i|<\mathrm{\&epsiv;},(\mathrm{\&epsiv;}=0.01),$

其中K表示聚类的类别数，取值为2，||x_j-c_i||表示x_j和c_i的欧式距离，c′_i表示下一次迭代后的聚类中心；

4f)标记最终的聚类结果，第一类标记为0，第二类标记为1，并将标记为1的像素转换为255，得到变化检测结果图。

Images