SAR图像NSCT域显著图去噪变化检测

慕彩红， 吴生财， 刘　逸， 彭　鹏， 刘若辰

(1. 西安电子科技大学 智能感知与图像理解教育部重点实验室，陕西 西安 710071; 2. 西安电子科技大学 电子工程学院，陕西 西安 710071; 3. 西安电子科技大学 物理与光电工程学院，陕西 西安 710071)

合成孔径雷达(Synthetic Aperture Radar， SAR)因其全天候、全天时探测能力和强大的对地穿透的能力，已成为变化检测的热门信息源．SAR图像变化检测技术[1-4]在气候变化、自然灾害评估、战场目标打击效果评估和军事侦察等领域有着广阔的应用前景．随着研究人员对SAR图像变化检测技术的深入研究，陆续提出了很多方法．通常情况下，根据是否利用先验知识，可以将SAR图像变化检测方法分为有监督变化检测和无监督变化检测方法两类，其中以无监督变化检测方法较为常见．文献[5]提出一种基于小波融合和主成分分析(Principal Component Analysis，PCA)-核模糊聚类的遥感图像变化检测方法，对光谱遥感图像和SAR图像均可取得良好的变化检测效果; 文献[6]将邻域信息粒子群聚类用于SAR图像变化检测，充分利用了像素的邻域信息，增强了抗噪性能，提高了变化检测精度; 文献[7]提出了一种基于显著图导向的SAR图像变化检测方法，将视觉感知技术应用到变化检测领域，取得了不错的效果．近年来，利用图像显著性进行SAR图像变化检测的方法引起了一些研究人员的兴趣．所谓图像的显著性，就是模拟生物视觉关注机制的一种选择性注意模型，当人观察一幅图片时总是第一时间无意识地将目光聚焦在一些具备显著特征的区域，这一生物观察机制可以很好地应用于图像处理的许多领域．文献[8]将显著性提取的方法应用到探测高分辨卫星图像中的感兴趣区域，文献[7]将显著性信息应用到了SAR图像变化检测领域．但是完全以显著信息为导向的变化检测技术还是有其固有的缺点．图像显著性获取的原理与人眼在关注图像时把注意力集中在特征明显区域的原理相类似，当图像的鲜明特征区域相对集中时，就能提取到完整的显著图，当图像的鲜明特征区域比较分散时，往往只会提取到鲜明特征最集中和最显著的区域，而忽略偏远和边缘的一些区域，从而造成显著信息的缺失．如果以这样的显著图引导变化检测，就会从源头上漏检部分像素，影响检测精度．为解决上述问题，笔者采用图像融合的策略，在非下采样轮廓变换(NonSubsampled Contourlet Transform， NSCT)[9]域运用显著图信息对邻域对数比图的低频子带进行范围限定，提高均值比图低频系数所占比例，以突出融合差异图的变化区域; 对这两幅差异图的方向子带则采用显著图进行选择性的去噪，再采用局部能量最小原则进行融合，以抑制融合差异图的背景区域; 最后对融合差异图进行k均值聚类，得到检测结果图．该方法综合利用了均值比图的变化区域接近实际的特点和邻域对数比图的背景区域平滑的特点，且对NSCT分解后的低频子带和方向子带运用显著信息设计了有利于提高融合差异图质量的融合策略，确保了变化区域的完整，为获取精确的变化检测结果图打下了良好的基础．图1是文中方法的流程图，图形⊙表示引入显著图信息，从图1中可以看出，在低频子带和方向子带的融合过程中均引入了显著图信息．

图1　文中方法的流程图

1　显著图的提取

用已配准和校正的两时刻SAR图像I1和I2分别构造均值比差异图[5]D1，对数比差异图[5]D和邻域对数比差异图D2，其中邻域对数比差异图D2的构造方法如下: 首先对两幅输入图像中以对应像素(i，j)为中心的 3×3 局部邻域进行大小为3、标准偏差为5的高斯滤波，得到邻域矩阵l1和l2；其次，计算局部邻域内对应像素灰度值的对数比，得到局部邻域对数比矩阵Zlog在像素(m，n)的值，其中，m，n是局部邻域内的像素的横纵坐标；接着，遍历l1和l2邻域内所有对应像素，得到局部邻对数比矩阵Zlog，以Zlog的均值d2(i，j)作为邻域对数比差异图D2在像素(i，j)处的值．按上述步骤遍历两幅图像所有对应像素，得到邻域对数比差异图D2= {d2(i，j)}．

下面介绍显著图的提取．利用上下文感知的显著性检测方法[10]从对数比值图D中提取显著图．将对数比值图D切分成多个块，将以i像素为中心的块xi与所有的其他块xj在Lab颜色空间作对比，如果xi和xj之间的差距大，则说明是显著性特征．dv(xi，xj)表示xi和xj之间的颜色距离，取值范围为[0，1]，即为两图像块在Lab空间的欧几里得颜色距离;dp(xi，xj)表示向量xi和xj之间的空间距离，取值范围为[0，1]，即为两图像块的欧几里得位置距离．一对图像块xi和xj之间的颜色距离越大，位置距离越小，则它们的差异值越大，得到颜色距离和空间距离以后，则两个块向量之间的非相似性FS可表示为

FS(xi,xj)=dv(xi,xj)/(1+cdp(xi,xj)),

(1)

其中，c=3．为了判断精准，以任意像素i为中心的图像块xi需要与多个尺度的A个图像块去计算显著值，背景区域在多个尺度上都是和其他块相似的，所以在多个尺度上计算可以降低背景区域的显著性．以任意像素i为中心的图像块xi在单尺度h上的显著值为

(2)

按下式整合所有尺度h上的显著值，并加入上下文修正，便得到该像素在对数比图D中的显著值，即

(3)

其中，H包括{100%， 80%， 50%， 30%}4个等级，M是等级数 (M=4);[·]表示把当前规模的图像尺寸恢复到原图像的尺寸大小;dfoci，h(i)表示在h尺度上像素i与离它最近的焦点像素之间的欧式距离，其值的取值范围为[0，1]; 文中图像块数量A=64，图像块尺寸为 7×7，选取时有50%的像素重叠，遍历所有像素便得到了对数比值图D的显著区域图S= {Si}，更多详细内容见文献[8]．最后采用最大类间方差法[11]从显著区域图S中提取显著图Y，并将Y中显著区的像素值转化为1，非显著区像素值转化为0．当用显著图Y与任意图像矩阵相乘后，显著区内的像素保留原图像像素的灰度值，显著区外的像素值则全部为零，有判别性信息的显著区被完好保存，而非显著区中的噪声则被过滤．

图2　NSCT的结构示意图

2　融合差异图的构造及变化检测

2.1　非下采样轮廓变换

NSCT是一种利用迭代非下采样滤波器组来实现一系列多分辨率、多方向和平移不变的频域子图像的算法，具有良好的多方向性、各向异性和平移不变性，能够很好地保持图像的奇异信息，细节保持能力强．其分解过程可表述为: 先利用非下采样塔式滤波器组对图像进行多尺度的分解，得到各种不同频率的高频子带图像和一个低频子带图像，然后再利用非下采样方向滤波器组对各高频子带图像进行多方向分解，从而得到不同尺度、不同方向的子带系数，如图2所示．SAR图像经NSCT分解后除低频子带中包含的少量噪声外，绝大部分噪声主要包含在方向子带中，且分解级数越高，其方向子带包含的细节信息和噪声越多．

2.2　融合差异图的构造

对均值比图D1和邻域对数比图D2进行3级NSCT变换，分别产生1个低频子带和14个方向子带．其中低频子带包含了图像的轮廓信息，能够充分体现出变化区域，邻域对数比图的变化区域与实际的变化不相符，而均值比图的变化区域接近于实际情况，因此在低频融合时就增加均值比图的低频子带系数的比例，减少邻域对数比图的低频子带系数的比例，以确保融合差异图的变化区域更接近于实际情况; 方向子带包含的是图像的细节信息和噪声，由于不同的SAR图像受噪声污染的程度不同，引入滤波参数K来决定对哪几层对应的方向子带进行显著图去噪，设定当K=1 时，对第3层的8个方向子带去噪;K=2 时，对第2层和第3层的12个方向子带进行去噪;K=3 时，对所有层的14个方向子带进行去噪，每幅图像经过3次试验就可以得到K的最佳取值．

2.2.1低频子带融合

首先对均值比图的低频子带系数进行高斯滤波以过滤低频子带中的噪声邻域对数比图在构造过程中已进行了3次滤波，且所占比例小，因此不进行滤波．其次对邻域对数比图的低频子带按下式用显著图进行范围限定:

D2, LY=D2,L⊙Y,

(4)

其中，D2, L表示邻域对数比图的低频子带系数，Y为显著图，⊙为矩阵对应元素相乘运算，D2, LY为经显著图限定后的邻域对数比图的低频子带系数，该系数在显著区内为邻域对数比图的低频子带的系数值，非显著区的系数值均变成0．然后按下式对低频子带系数进行融合:

DL=0.5D1, L+0.5D2, LY,

(5)

其中，D1, L为均值比图的低频子带系数，DL为融合后的低频子带系数，该系数在显著区外为均值比图的低频子带系数，在显著区域内为融合系数．均值比图的低频子带系数占有比例多，利于突出变化区域的信息．

2.2.2方向子带融合

图像经NSCT分解后的方向子带包含的是图像的细节信息和噪声，其中第3层的8个方向子带包含的细节和噪声信息最多，前2层相对较少．首先，按滤波参数K对两幅差异图的方向子带进行显著图去噪以滤除方向子带的非显著区域中的噪声．其表达式为

(6)

其中，D(1,2), H, r, n表示均值比图D1或邻域对数比图D2中第r层的第n个方向子带系数．当滤波参数K的取值为{1， 2， 3}时，就对K所对应分解层数的方向子带系数进行显著图去噪．

其次，按局部能量最小原则对均值比图D1和邻域对数比图D2的方向子带系数D1, H, r, n和D2, H, r, n进行融合．采用局部能量最小原则可以合并来自两幅差异图方向子带的同质区域，起到抑制背景信息和去噪的作用．使融合差异图的背景区域保持平滑局部能量最小原则的表达式为

(7)

其中，r={1，2，3}，表示分解的层数;n={1， 2， 3，…，13， 14};DH，r，n(i，j)表示融合后的方向子带系数在像素(ij)处的值;D1，H，r，n(i，j)和D2，H，r，n(i，j)分别表示均值比图D1和邻域对数比图D2在第r层的第n个方向子带系数在像素(ij)处的值;E1，r，n(i，j)和E2，r，n(i，j)分别表示两幅差异图在第r层的第n个方向子带系数中以像素(i,j)为中心的局部能量系数值，按下式计算:

(8)

其中，Nij表示以像素(i,j)为中心的3×3局部邻域窗口的大小，D(1,2), H, r, n(B)表示均值比图D1和邻域对数比图D2中第r层第n个方向子带系数在邻域窗口Nij内的第B个元素的值．

2.2.3NSCT反变换和k均值聚类

对融合后的1个低频子带DL和3个分解层的14个方向子带DH, r, n进行NSCT反变换后，得到融合差异图Xd，该图兼顾了均值比图的变化区域接近实际情况的特点和邻域对数比图的背景区域(非变化区域)相对平滑的特点，且抗噪能力好．最后采用k均值聚类方法对融合差异图Xd进行两分类，输出变化检测结果图．

3　对比实验结果分析与有效性验证

3.1　实验数据集介绍

第1个数据集是分别于1999年4月和1999年5月获得的伯尔尼地区两时相SAR图像，尺寸均为 301× 301，像素灰度级为256，图3(a)是洪水侵袭前的图像，图3(b)是洪水侵袭后的图像，图3(h)是标准检测结果图，包括 1 269 个变化像素点．第2个数据集是星载SAR传感器分别于1997年5月和1997年8月拍摄的渥太华地区的两时相SAR图像，尺寸均为 290× 250，灰度级为256，图4(a)是洪水侵袭前的图像，图4(b)是洪水侵袭后的图像，图4(h)为标准检测结果图，包括 16 049 个变化像素点．笔者用漏检数误检数总错误数和Kappa系数共4个定量分析指标[12]来评价检测结果．

图3　伯尔尼地区SAR图像数据集变化检测结果

图4　渥太华地区SAR图像数据集变化检测结果

3.2　对比实验结果和分析

用文中方法NSCT域显著图去噪(Saliency map Denoising in Non subsampled contourlet transform Domain，SDND)与文献[7]提出的“Unsupervised saliency-guided SAR image change detection”(USG)和“Using Combined Difference Image and K-means Clustering for SAR Image Change Detection[13]”(CDI-K)，文献[14]提出的“Wavelet Fusion on Ratio Images for Change Detection in SAR Images”(WFD-FL)以及文献[15]提出的“混合的SAR图像变化检测算法[15]”(NSCTF)在两个数据集上进行对比实验，其实验结果如表1、图3和图4所示．

表1　文中方法与现有方法的变化检测结果比较

从图3、图4和表1可以得出如下结论，文中方法在2个SAR图像数据集上的所检测结果图在直观效果上与标准检测图最接近，USG法和CDI-K法的漏检像素数多变化区域不完整; WFD-FL法和NSCT法的检测结果图在直观效果上都存在边缘区域和细节信息不清晰，特别是NSCTF方法的检测结果图基本看不出变化区域的边缘和细节，只能看清变化区域的轮廓，与标准检测图差距很大．文中方法的漏检数与误检数相对均衡，总错误数和Kappa系数均优于4种现有的对比方法，主要原因: 一是文中方法采用融合策略，综合利用了均值比图的变化区域接近实际的特点和邻域对数比图的背景区域相对来说平滑的特点; 二是采用NSCT变换构造融合差异图，利于增强去噪效果和细节结构的保持能力; 三是构造融合差异图的过程中对低频子带系数所占有的比例用显著图进行分配调整以突出变化区域，对方向子带中的噪声用显著图进行过滤并采用局部能量最小原则进行融合以抑制背景信息，从而使融合后的差异图既接近于实际变化情况，又具备很强的抗噪能力，有利于提高变化检测精度．

4　结 束 语

文中提出了一种在NSCT域结合显著图信息构造融合差异图的SAR图像变化检测方法．首先分别构造均值比、对数比和邻域对数比差异图，并用对数比图获取显著图; 然后对均值比图和邻域对数比图进行3级NSCT分解，低频融合时对邻域对数比图的低频子带用显著图进行范围限定，以突出融合差异图的变化区域; 高频融合时对两幅差异图的方向子带进行选择性的显著图限定去噪，再采用局部能量最小的原则进行融合，以抑制融合差异图的背景区域; 最后经过NSCT反变换，得到融合差异图并进行k均值聚类，输出检测结果图．实验证明，文中方法抗噪能力强，边缘细节保持好，检测精度高．

参考文献：

[1] VU V T, PETTERSSON M I, MACHADO R, et al. False Alarm Reduction in Wavelength-resolution SAR Change Detection Using Adaptive Noise Canceler[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(1): 591-599.

[2]LI H, GONG M, WANG Q, et al. A Multiobjective Fuzzy Clustering Method for Change Detection in SAR Images[J]. Applied Soft Computing Journal, 2016, 46: 767-777.

[3]GONG M G, ZHAO J, LIU J, et al. Change Detection in Synthetic Aperture Radar Images Based on Deep Neural Networks[J]. IEEE Transactions on Neural Networks and Learning Systems, 2016, 27(1): 125-138.

[4]MU C H, LI C, LIU Y, et al. Change Detection in SAR Images Based on the Salient Map Guidance and an Accelerated Genetic Algorithm[C]//Proceedings of the 2017 IEEE Congress on Evolutionary Computation. Piscataway: IEEE, 2017: 1150-1157.

[5]慕彩红, 霍利利, 刘逸, 等. 基于小波融合和PCA-核模糊聚类的遥感图像变化检测[J]. 电子学报, 2015, 43(7): 1375-1381.

MU Caihong, HUO Lili, LIU Yi, et al. Change Detection for Remote Sensing Images Based on Wavelet Fusion and PCA-kernel Fuzzy Clustering[J]. Acta Electronica Sinica, 2015, 43(7): 1375-1381.

[6]刘逸, 慕彩红, 刘敬. 结合邻域信息粒子群聚类用于SAR图像变化检测[J]. 西安电子科技大学学报, 2015, 42(1): 187-193.

LIU Yi, MU Caihong, LIU Jing. Change Detection for SAR Images Based on the Particle Swarm Clustering Algorithm Using Neighborhood in Formation[J]. Journal of Xidian University, 2015, 42(1): 187-193.

[7]ZHENG Y G, JIAO L C, LIU H Y, et al. Unsupervised Saliency-guided SAR Image Change Detection[J]. Pattern Recognition, 2017, 61: 309-326.

[8]ZHANG F, DU B, ZHANG L. Saliency-guided Unsupervised Feature Learning for Scene Classification[J]. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53(4): 2175-2184.

[9]da CUNHA A L, ZHOU J, DO M N. The Nonsubsampled Contourlet Transform: Theory, Design, and Applications[J]. IEEE Transactions on Image Processing, 2006, 15(10): 3089-3101.

[10]GOFERMAN S, ZELNIK-MANOR L, TAL A. Context-aware Saliency Detection[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2012, 34(10): 1915-1926.

[11]OTSU N. A Threshold Selection Method from Gray-level Histograms[J]. IEEE Transactions on System, Man and Cybemetics, 1979, 9(1): 62-66.

[12]付明柏. 基于UDWT与Snake模型的多时相SAR图像变化检测方法[J]. 计算机应用与软件, 2014, 31(4): 254-257.

FU Mingbai. Multitemporal SAR Images Change Detection Method Based on UDWT and Snake Models[J]. Computer Applications and Software, 2014, 31(4): 254-257.

[13]ZHENG Y G, ZHANG X, HOU B, et al. Using Combined Difference Image and,k-means Clustering for SAR Image Change Detection[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 11(3): 691-695.

[14]MA J J, GONG M, ZHOU Z. Wavelet Fusion on Ratio Images for Change Detection in SAR Images[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(6): 1122-1126.

[15]李玲玲, 贾振红, 覃锡忠, 等. 混合的SAR图像变化检测算法[J]. 计算机工程与设计, 2015, 36(5): 1256-1259.

LI Lingling, JIA Zhenhong, QIN Xizhong, et al. Hybrid SAR Image Change Detection Algorithm[J]. Computer Engineering and Design, 2015, 36( 5): 1256-1259.