专利名称：基于邻域相似性及掩模增强的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图像I1和I2进行几何校正和配准；(2)根据邻域相似度NSR算子，按照如下步骤构造图像I1和I2的差异影像图DI ；2a)分别取得两时相SAR图像IJPI2在同一位置χ上的像素灰度值邻域集合N1(X) 和队⑴，其大小为NXN，N e {3,5,7,9}；2b)比较两个邻域集合N1(X)和队00的相似性，得到差异图DI在位置χ上的像 素灰度值DI(X) 其中，N1(Xi) ^PN2(Xi)分别表示图像I1和I2在位置χ上邻域集合N1(X) ^P N2(χ) 的第i个元素，集合N1(X)与N2(X)越相似，则DI(X)值越大，图像I1和I2在位置χ上的像 素点属于非变化区域的可能性越大，反之，属于变化区域的可能性越大；2c)对图像I1和I2的每一位置χ从左到右，从上到下重复步骤2a)和2b)，得到差 异影像图DI ；(3)对差异影像图DI进行掩模增强处理，得到新的差异影像图NDI，所述的掩模增强处理包括如下步骤3a)根据差异影像图DI的图像大小选用掩模，若DI小于1000X 1000，则选用中心 系数为_8的掩模，若DI大于1000 X 1000，则选用中心系数为-9的掩模；3b)利用所选的掩模，增大DI在位置χ上像素点与其周围相邻像素点间的差异，得 到新的差异影像图在位置X上的值NDI(X) 其中LN表示掩模的长度，取值为9，L(Xi)为掩模中的第i个数的值，DI(Xi)表示 DI中大小为3X3的χ邻域集合中的第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图像I1和 I2进行几何校正和配准。步骤二 根据邻域相似度NSR算子，按照如下步骤构造图像I1和I2的差异影像图 DI。2a)分别取得两时相SAR图像I1和I2在同一位置χ上的像素灰度值邻域集合N1 (χ) 和队⑴，其大小为NXN，N e {3,5,7,9}；2b)比较两个邻域集合N1(X)和队00的相似性，得到差异图DI在位置χ上的像 素灰度值DI(X) 其中，N1(Xi)和N2 (Xi)分别表示图像I1和I2在位置χ上邻域集合N1(X)和N2(X) 的第i个元素，集合N1(X)与N2(X)越相似，则DI(X)值越大，图像I1和I2在位置χ上的像 素点属于非变化区域的可能性越大，反之，属于变化区域的可能性越大；2c)对图像I1和I2的每一位置χ从左到右，从上到下重复步骤2a)和2b)，得到差 异影像图DI。步骤三对差异影像图DI按如下步骤进行掩模增强处理，得到新的差异影像图 NDI。3a)根据差异影像图DI的图像大小选用掩模，如附图2所示，其中图2 (a)是中心 系数为_8的掩模，适用于DI小于1000X 1000的情况，图2 (b)是中心系数为_9的掩模，适 用于DI大于1000X1000的情况；3b)利用所选的掩模，增大DI在位置χ上像素点与其周围相邻像素点间的差异，得 到新的差异影像图在位置X上的值NDI(X)
LNNDI(X) = X(Z)/(x,) - L(X1)), \<i<LN
'=I’其中LN表示掩模的长度，取值为9，L(Xi)为掩模中的第i个数的值，DI(Xi)表示 DI中大小为3X3的χ邻域集合中的第i个像素的灰度值；3c)对DI的从左到右，从上到下重复步骤3b)，得到新的差异影像图NDI。步骤四利用初始点确定的K均值方法，对新的差异影像图NDI进行聚类，得到变 化检测结果图⑶I。4a)选用两个初始聚类中心为C1 = graymin+ (graymax-graymin) /4c2 = graymin- (graymax-graymin) /4,其中，graymin和graymax分别表示NDI中所有像素点灰度的最小值和最大值；4b)根据初始聚类中心，对NDI中所有的像素点进行首次分类，计算每个像素点到 两个聚类中心C1和C2之间的欧氏距离，并将该像素点归到与它最近的那一类；4c)根据分类结果，计算新的聚类中心Ci, 其中Ci表示第i类的新的聚类中心，Ni表示第i类中像素点的总个数，Xj表示第 i类中的第j个像素点；4d)根据新的聚类中心Ci,对NDI中所有的像素点再次分类，计算每个像素点到新 的聚类中心Ci，i e {1,2}之间的欧氏距离，并将该像素点归到与自己最近那一类；4e)重复步骤4c)和4d)直到聚类结束，聚类结束的停机准则如下 其中K表示聚类的类别数，取值为2，Mxj-CiII表示XjPci的欧式距离，C' 1表 示下一次迭代后的聚类中心；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，大小均为290X 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所示。
8
表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，大小均为7692X7666。双边滤波 处理和掩模增强处理的时间对比如图12所示，实线和虚线分别代表了双边滤波处理与掩 模增强处理随着图像不断增大时的运行时间。从图12可以看出，在图像小于800X800时，掩模增强处理的时间与双边滤波处 理时间相差不大，但是随着图像的不断增大，双边滤波处理时间几乎成指数级增长，而掩 模增强处理时间却以非常缓慢的速率在增加，如对于第二组仿真图的双边滤波处理时间为
91596. 844s，而掩模增强处理时间仅为18. 562s，且取得了与双边滤波一样的，甚至更好的检 测效果，如图10和图11所示。放大它们的代表性区域A，B和C后的原图如图13所示，它 们的变化检测结果图，如图14和15所示，其中A，B和C区域分别反应的是2008年到2009 年间黄河入海口地域的陆地，港口和海岸的变化。
权利要求
一种基于邻域相似性及掩模增强的SAR图像变化检测方法，包括如下步骤(1)对不同时相同一场景的两幅SAR图像I1和I2进行几何校正和配准；(2)根据邻域相似度NSR算子，按照如下步骤构造图像I1和I2的差异影像图DI；2a)分别取得两时相SAR图像I1和I2在同一位置x上的像素灰度值邻域集合N1(x)和N2(x)，其大小为N×N，N∈{3，5，7，9}；2b)比较两个邻域集合N1(x)和N2(x)的相似性，得到差异图DI在位置x上的像素灰度值DI(x) <mrow><mi>DI</mi><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow><mo>=</mo><mfrac> <mrow><munderover> <mi>Σ</mi> <mrow><mi>i</mi><mo>=</mo><mn>1</mn> </mrow> <mrow><mi>N</mi><mo>×</mo><mi>N</mi> </mrow></munderover><mi>min</mi><mrow> <mo>(</mo> <msub><mi>N</mi><mn>1</mn> </msub> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>)</mo> </mrow> <mo>,</mo> <msub><mi>N</mi><mn>2</mn> </msub> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>)</mo> </mrow> <mo>)</mo></mrow> </mrow> <mrow><munderover> <mi>Σ</mi> <mrow><mi>i</mi><mo>=</mo><mn>1</mn> </mrow> <mrow><mi>N</mi><mo>×</mo><mi>N</mi> </mrow></munderover><mi>max</mi><mrow> <mo>(</mo> <msub><mi>N</mi><mn>1</mn> </msub> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>)</mo> </mrow> <mo>,</mo> <msub><mi>N</mi><mn>2</mn> </msub> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>)</mo> </mrow> <mo>)</mo></mrow> </mrow></mfrac><mo>,</mo> </mrow>其中，N1(xi)和N2(xi)分别表示图像I1和I2在位置x上邻域集合N1(x)和N2(x)的第i个元素，集合N1(x)与N2(x)越相似，则DI(x)值越大，图像I1和I2在位置x上的像素点属于非变化区域的可能性越大，反之，属于变化区域的可能性越大；2c)对图像I1和I2的每一位置x从左到右，从上到下重复步骤2a)和2b)，得到差异影像图DI；(3)对差异影像图DI进行掩模增强处理，得到新的差异影像图NDI，所述的掩模增强处理包括如下步骤3a)根据差异影像图DI的图像大小选用掩模，若DI小于1000×1000，则选用中心系数为 8的掩模，若DI大于1000×1000，则选用中心系数为 9的掩模；3b)利用所选的掩模，增大DI在位置x上像素点与其周围相邻像素点间的差异，得到新的差异影像图在位置x上的值NDI(x) <mrow><mi>NDI</mi><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow><mo>=</mo><munderover> <mi>Σ</mi> <mrow><mi>i</mi><mo>=</mo><mn>1</mn> </mrow> <mi>LN</mi></munderover><mrow> <mo>(</mo> <mi>DI</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>)</mo> </mrow> <mo>-</mo> <mi>L</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>)</mo> </mrow> <mo>)</mo></mrow><mo>,</mo><mn>1</mn><mo>≤</mo><mi>i</mi><mo>≤</mo><mi>LN</mi><mo>,</mo> </mrow>其中LN表示掩模的长度，取值为9，L(xi)为掩模中的第i个数的值，DI(xi)表示DI中大小为3×3的x邻域集合中的第i个像素的灰度值；3c)对DI的从左到右，从上到下重复步骤3b)，得到新的差异影像图NDI；(4)利用初始点确定的K均值聚类方法，对新的差异影像图NDI进行聚类，得到变化检测结果图CDI。
2.根据权利要求1所述的SAR图像变化检测方法，其中步骤(1)所述的对不同时相同 一场景的两幅SAR图像I1和I2进行几何校正和配准，是通过遥感图像处理软件ERDAS完 成。
3.根据权利要求1所述的SAR图像变化检测方法，其中步骤(4)所述的利用初始点确 定的K均值方法，对新的差异影像图NDI进行聚类，按如下步骤进行；4a)选用两个初始聚类中心为 其中，grayfflin和graymax分别表示NDI中所有像素点灰度的最小值和最大值；4b)根据初始聚类中心，对NDI中所有的像素点进行首次分类，计算每个像素点到两个聚类中心C1和C2之间的欧氏距离，并将该像素点归到与它最近的那一类； 4c)根据分类结果，计算新的聚类中心Ci，Kl M,其中Ci表示第i类的新的聚类中心，Ni表示第i类中像素点的总个数，Xj表示第i类 中的第j个像素点；4d)根据新的聚类中心Ci，对NDI中所有的像素点再次分类，计算每个像素点到新的聚 类中心Ci，i e {1,2}之间的欧氏距离，并将该像素点归到与自己最近那一类； 4e)重复步骤4c)和4d)直到聚类结束，聚类结束的停机准则如下满足 ” 最小，且Σ K I (ε = 0.01) ’ 其中K表示聚类的类别数，取值为2，Mxj-CiII表示\和4的欧式距离，C' i表示下 一次迭代后的聚类中心；4f)标记最终的聚类结果，第一类标记为0，第二类标记为1，并将标记为1的像素转换 为255，得到变化检测结果图。
全文摘要
本发明公开了一种基于邻域相似性及掩模增强的SAR图像变化检测方法，它涉及遥感图像处理领域，主要解决SAR图像变化检测的差异图构造受SAR图像斑点噪声影响严重的问题，其实现步骤为(1)对不同时相同一场景的两幅SAR图像I1和I2进行几何校正和配准；(2)用邻域相似度算子，构造图像I1和I2的差异影像图DI；(3)对差异影像图DI进行掩模增强处理，得到新的差异影像图NDI；(4)利用初始点确定的K均值聚类方法，对新的差异影像图NDI进行聚类，得到变化检测结果图CDI。本发明具有运行效率高，算法复杂度和时间复杂度低，实验效果好的特点，可用于两时相SAR图像的变化检测。
文档编号G06T5/00GK101923711SQ20101023065
公开日2010年12月22日 申请日期2010年7月16日 优先权日2010年7月16日
发明者公茂果, 刘芳, 吴巧娣, 焦李成, 王桂婷, 王爽, 钟桦, 马文萍 申请人:西安电子科技大学