0 引言
图像变化检测是研究不同时段同一场景图像之间发生的变化[1-2],在灾害评估、土地利用、城市规划及军事侦察等方面具有广阔的应用前景。合成孔径雷达(Synthetic Aperture Radar,SAR)是一种主动式微波成像传感器[3-6],与传统光学遥感和高光谱遥感相比,SAR具有全天时、全天候、大区域对地观测的特点,因而在变化检测领域的应用越来越广。
在SAR图像变化检测技术中,无监督变化检测由于不需要先验变化信息的支持,成为国内外研究的重点[7]。目前,SAR图像变化检测主要有两类思想:1)分类后比较法,即先将两幅SAR图像分别进行分类,然后比较分类结果,从而得到变化信息。2)先获取差异图,然后对差异图进行分析,从而获取最终的变化信息。在这两类思想中,第一类思想的变化检测精度受分类结果的影响很大,且容易导致分类误差积累问题,而第二类思想的思路简单明了,且具有较大的研究空间,被广泛采用。但第二类思想仍存在以下3个问题:1)如何抑制斑点噪声。SAR是相干系统,斑点噪声是其固有特性,SAR图像中斑点噪声的存在无疑会影响变化检测的结果。SAR图像的这一特性使SAR图像变化检测成为变化检测领域极具挑战的难点和热点问题。2)如何构造性能良好的差异图。差异图可分离性的好坏直接影响变化检测的最终结果,因此构造性能良好的差异图成为这类方法的研究重点之一。3)如何设计一种有效的分类方法,用以区分变化信息和非变化信息。这是SAR图像变化检测中的重要步骤之一,会直接影响变化检测的结果。
传统的SAR图像变化检测方法仅使用单一特征,未能充分利用图像的细节信息,因此检测效果并不理想。近年来,多通道特征如ICF(Integral Channel Features),ACF(Aggregate Channel Features)等被成功应用于行人检测[8]及人脸检测[9-10]领域。多通道特征从不同角度集合了多种通道特征,图像信息更加丰富。基于上述内容,本文提出一种基于多通道特征的无监督SAR图像变化检测方法。首先,引入引导滤波来抑制相干斑噪声并尽可能地保留SAR图像的边缘及局部信息;其次,提取8个通道特征来获得性能良好的差异图;最后,利用主成分分析(PCA)和K-means聚类得到最终的变化信息。
1 本文算法
本文算法可以分为3个部分:引导滤波、基于多通道特征的差异图获取和差异图分析。算法整体流程如图1所示。
1.1 引导图像滤波
为了在抑制相干斑噪声的同时尽可能多地保留SAR图像的边缘及局部信息,本文引入引导图像滤波[11]对SAR图像进行预处理。
引导图像滤波的本质是将空域滤波模型推广到基于图像信息的滤波,它包括引导图像I、输入图像p、输出图像q,其中I和p是根据具体应用事先就给定的。
假设输出图像q与引导图像I在以k为窗后中心的窗口wk中满足线性关系,如下所示:
图1 本文算法流程图
式中,(ak,bk)是当窗口为wk时的线性系数。假定窗口大小为r,这一个线性的关系就保证了只有当I有边缘时,输出q才有边缘,这是因为∇q=a∇i。这一模型在图像抠图[12]、超分辨率图像[13]、图像去雾[14]等应用中得到有效的证实。
为了确定线性系数的大小,使得输出q和输入p之间的差异最小,即将下式最小化:
式中,ε是用来防止ak过大的正则化因子,则上式得解为
式中,uk和
分别为I在窗口wk下的均值和方 差,|w|为窗口中像素的个数,
接下来计算整幅图像每个窗口的线性系数,但是由于一个像素会被多个窗口包含,使得式(1)中qi的值并不一样,因此在计算每一个输出时可以采用计算包含该点的所有线性系数的平均值。计算公式如下:
式中,
。改进之后,由于线性系数
分布于不同的空间,使得∇q不再与∇i成比例,但是由于是均值滤波的输出,使得在I中明显边缘处的梯度变小了许多,在这种情况下可以认为∇q≈a∇i,这也证明了为什么引导滤波可以较好地保持边缘特性。
实际上,可以将式(3)重新写成p的加权和:ak=∑jAkj(I)pj,其中Aij是仅依赖于I的权重,同理:bk=∑jBkj(I)pj,qi=∑jWij(I)pj。当把引导滤波用作边缘保持滤波器时,往往有I=p,如果ε=0,显然ak=1,bk=0是E(ak,bk)为最小值的解,这时的滤波器没有任何作用,输出原始图像。如果ε＞0,在像素强度变化小的区域(或单色区域),有a近似于(或等于)0,而b近似于(或等于)
,相当于对图像作了一个加权均值滤波;而在像素强度变化大的区域,a近似于1,b近似于0,对图像的滤波效果很弱,有助于保持边缘。ε的作用是决定滤波效果,在窗口大小不变的情况下,随着ε的增大,滤波效果越明显。
综上所述,引导滤波具有与双边滤波那样保持边缘平滑的特性,不受梯度反转伪影影响,所以它在细节处理上比双边滤波具有更好的效果,被广泛应用于光学图像的去噪、细节平滑、图像增强和图像去雾等领域。引导滤波最大的优势在于,可以写出时间复杂度与窗口大小无关的算法,因此在使用大窗口处理图片时,其效率更高。由于其良好的边缘保持特性和较快的计算效率(时间复杂度为O(N)),本文将其应用到SAR图像去噪,并取得了良好的效果。
1.2 基于多通道特征的差异图获取
在变化检测中,差异图获取是一个关键环节,差异图的好坏影响最终的变化检测结果。为获取性能良好的差异图,本文提出了基于多通道特征的差异图获取。首先对预处理后的SAR图像提取8个通道特征,然后利用对数比值法获取8个通道特征的差异图,最后对各通道差异图进行融合,得到最终的差异图。采用多个通道特征获取差异图,可以从不同角度集成多种特征,图像信息更丰富,更有利于区分变化信息和非变化信息。
1.2.1 通道特征
图像的通道广义上可以看作是原始图像的某种映射,它的像素值是由原始图像的对应区域通过某种计算得到的[8]。给定一个输入图像I,其通道图像C的计算公式如下:
式中,Ω代表图像的某种通道计算函数。
通道类型有多种,如颜色通道、灰度、线性滤波、非线性变换、逐点变换、积分直方图、梯度直方图等[8]。本文选取8个通道特征用于获取差异图,包括1个灰度通道、1个归一化的梯度幅值通道和6个梯度方向直方图通道。在选取的8个通道特征中:
1)灰度是最简单的通道。对于SAR图像而言,由于其本身就是一张灰度图,所以,SAR图像的灰度通道C=I,即SAR图像的原图,如图2(a)所示。
2)梯度幅值通道属于非线性变换通道的一种,如图2(b)所示,梯度幅值通道可以捕获边缘强度信息。
3)梯度方向直方图是一个加权直方图,如图2(c)～(h)所示。梯度方向取0°～360°,然后将梯度分布平均分成6个方向角度,即每60°取一个方向角度。梯度方向直方图计算公式如下:
式中,G(x,y)代表图像I(x,y)处的梯度幅值,Θ(x,y)代表图像I(x,y)处量化的梯度方向。通过对原始SAR图像I进行不同尺度的模糊,可以计算出不同尺度的梯度信息。此外,借助于梯度幅值信息,可以对计算出来的梯度直方图进行归一化。
本文所使用的多通道特征提取流程如图2所示。
从图2可以看出,灰度通道和梯度通道的大小和原图相同,而6个梯度方向直方图的大小仅为原图的1/3,这是由于在计算梯度方向直方图时使用了3×3的单元格。为了得到和原图大小相同的梯度方向直方图,方便后期进行差异图融合,本文对其进行插值处理。从而最终提取到所使用的8个通道特征。
之所以选取这8个通道特征,是因为它们不仅提取方法操作简单,可以快速高效地计算出来,而且包含丰富的图像细节信息,可以从不同角度收集输入图像的信息,有利于抑制斑点噪声,获得性能良好的差异图,从而提高检测精度。
图2 多通道特征提取流程图
1.2.2 差异图融合
对数比值法在获取差异图方面有独特的优势,它不仅能够很大程度上消除乘性噪声的影响,减少定标处理引起的额外误差,而且能够突显出SAR图像的相对变化区域[15]。为充分利用不同的通道特征信息,生成基于多通道特征的融合差异图,本文利用对数比值法获取每一个通道特征的差异图,计算公式如下:
式中,Ri表示各通道特征的差异图,I0i表示t0时刻第i个通道图像的幅度,I1i表示t1时刻第i个通道图像的幅度。
在获取各通道特征的差异图后,取各通道差异图的平均值作为最终的差异图,计算公式如下:
1.3 差异图分析
在获取性能良好的差异图后,为了有效地区分变化信息和非变化信息,得到精确、可靠的变化检测结果,本文采用PCA和K-means聚类对差异图进行分析[16]。具体步骤如下:
1)将差异图R分成n×n的互不重叠的块(假设分成了m个小块),将每一小块转化为行向量Xi(i=1,2,…,m),用这些行向量构建矩阵F,矩阵大小为m×n2。分块的思想使得算法具有一定的抗噪性能。
2)用PCA将步骤1)中矩阵F的行向量维度从n2维降成S维,得到对应的变换矩阵Z。
3)取差异图R中每个像素所在的n×n邻域小块,将每个小块转化成行向量,共有H×W个行向量(H×W为差异图R的大小)。用得到的变换矩阵Z将所有的行向量变换到特征向量空间Q上(每个像素就由S维的向量表示),Q的维度为S×HW。
4)用K-means聚类算法将特征向量空间Q聚为两类:变化类和非变化类,从而得到最终的变化检测结果。
2 实验结果
2.1 实验数据集
本文选取两组配准后的数据集进行算法验证及实验分析。
数据集I是由RADARSAT-1卫星拍摄的加拿大渥太华地区的两景SAR图像(分别如图3(a)、图3(b)所示),拍摄时间分别为1997年7月和1997年8月,图像大小为290×350,图3(c)为变化信息参考图。
数据集II是由C波段ERS-2卫星拍摄的伯尔尼地区的两景SAR图像(分别如图4(a)、图4(b)所示),拍摄时间分别为1999年4月和1999年5月,图像大小为301×301,图4(c)为变化信息参考图。
图3 数据集Ⅰ:渥太华地区SAR图像
图4 数据集Ⅱ:伯尔尼地区SAR图像
之所以选取这两组数据集,首先是因为这两组数据集都存在较多的斑点噪声,其次是因为这两组数据集具有对比性,在数据集I的图像中有较多的变化像素,而数据集II的图像中变化像素较少。
2.2 评价指标
在评价指标方面,本文选取错检数(FP)、漏检数(FN)、总体错误数(OE)、百分比修正(PCC)、Kappa系数这5个检测指标来评价算法性能,其中:
1)FP:未变化像素被错检为变化像素的个数,即错检数(或称虚警数)。
2)FN:变化像素被检测为未变化像素的个数,即漏检数。
3)OE:总体错误数,其值为错检数与漏检数之和。
4)PCC:百分比修正,计算公式如下:
式中,TP为正确检测到的变化像素的数目,TN为正确检测到的未变化像素的数目。PCC的值越高,说明算法检测效果越好。
5)Kappa系数:计算公式如下:
式中,
式中,Mc为实际的变化像素的数目,Mu为实际的非变化像素的数目。
Kappa系数的值越高,说明算法的检测精度越高。
2.3 实验结果及分析
2.3.1 检测结果分析
为验证算法的性能,本文算法GF_MCF将与以下3种算法进行对比:文献[16]中的算法命名为PCAKM;文献[17]中的算法命名为MRFFCM;文献[18]中的算法命名为MOHO。
两组数据集的变化检测结果如图5、图6所示。其中,图5(a)和图6(a)显示的是不同算法提取的变化信息。图5(b)和图6(b)展示了更具体的检测结果。具体评价指标如表1所示。
图5 渥太华地区不同算法检测结果
图6 伯尔尼地区不同算法的检测结果
表1 算法性能指标
1)数据集I(渥太华地区)实验结果分析
从图5可以看出,本文算法GF_MCF在对数据集I的检测效果上要明显优于另外3种算法,不仅有较少的漏检和错检像素,且保持了更多的图像细节。这是由于GF_MCF算法结合了多个通道特征的优点,图像细节信息更加丰富,因此检测效果更好。从表1可以看出,GF_MCF算法在OE,PCC和Kappa系数上要明显优于另外3种算法,进一步说明了GF_MCF算法的检测优势。
2)数据集II(伯尔尼地区)实验结果分析
数据集II与数据集I的特点不同,数据集II的变化区域较小且比较集中,图像细节更容易被错检和漏检。从图6可以看出,由于斑点噪声的影响,PCAKM,MRFFCM和MOHO这3种算法的检测结果中都出现了孤立点,而GF_MCF算法的检测结果无孤立点产生,对斑点噪声的抑制明显强于另外3种算法。这不仅是因为GF_MCF算法中使用了引导滤波,还因为多通道特征的思想对斑点噪声也有一定的抑制作用,从而提高了算法的抗噪性能。同时,从表1可以看出,GF_MCF的FN和OE最低,PCC和Kappa系数最高,检测效果最好。
3)通道特征分析
本文采用了3种类型的通道特征:灰度通道、梯度幅值通道和梯度直方图通道(6个)。图7展示了不同通道特征之间相互组合后对检测结果的影响,其中n代表差异图分析部分中的窗口大小。从图7可以看出,结合了3种类型通道特征的检测效果最好。此外,从图中还可以看出,当窗口n=3时,PCC的值最高,检测效果最好。此组实验再次证明了通道特征在SAR图像变化检测领域的有效性。
图7 通道特征组合分析
2.3.2 引导图像滤波性能分析
为了验证引导图像滤波在SAR图像变化检测领域的适用性,本文在PCAKM,MRFFCM和MOHO算法中加入了引导滤波,分别命名为GF_PCAKM,GF_MRFFCM和GF_MOHO,检测性能指标如表1所示。从表1可以看出,加入引导滤波后,PCAKM,MRFFCM和MOHO算法的检测性能都有明显的提升。
同时,为进一步说明引导滤波的性能及有效性,本文将引导图像滤波与其他滤波算法在检测性能和滤波时间等指标方面进行对比,具体结果如表2所示。
表2 各滤波算法性能指标
从表2可以看出,引导滤波和其他滤波算法相比,OE最低,PCC和Kappa系数最高,并且滤波所耗时间最短,说明引导图像滤波的检测性能明显优于另外几种滤波算法。
在引导滤波中,可调参数有2个:窗口大小r和正则化因子ε。图8展示了引导滤波的2个参数对检测效果的影响。从图8可以看出,当窗口大小r=13,正则化因子ε=0.007时,PCC最高,检测效果最好。
图8 引导图像滤波参数对检测效果的影响
综上所述,引导图像滤波对斑点噪声有一定的抑制作用,适用于SAR图像变化检测,且检测效果较好。
2.3.3 算法抗噪性能分析
以上两部分的实验分析表明,本文算法具有良好的抗噪性能。为了进一步验证本文算法GF_MCF对斑点噪声的抑制能力,本文给出了不同算法的抗噪性能曲线,如图9所示。
给定一个输入SAR图像X,大小为H×W,加入斑点噪声后的图像计为
。加入的噪声水平用峰值信噪比(PSNR)来衡量,PSNR计算公式如下(假设输入SAR图像的最大像素值为1):
给定两个SAR图像X1,X2,用变化检测算法得到检测结果图CM1,然后往X1中加入斑点噪声,X2保持不变,用同一种变化检测算法得到检测结果图CM2。CM1和CM2之间的差异可用参数τ表示,计算公式如下:
τ的大小代表算法对斑点噪声的抑制能力,其值越高则算法的抗噪性能越好。
图9为不同算法的抗噪性能曲线。从图9可以看出,本文算法GF_MCF对斑点噪声具有较好的抑制能力,这也是本文算法的优势之一。
图9 算法抗噪性能曲线
3 结束语
本文提出了一种基于多通道特征的SAR图像变化检测方法。实验结果表明,该方法对SAR图像的变化具有良好的检测效果,检测精度进一步提高且具有较好的抗噪性能。在未来的工作中,拟进一步对通道特征进行探索,并考虑用两级聚类分析差异图,以得到更为精确的检测结果并进一步提高算法的鲁棒性。
参考文献:
[1]GONG Maoguo,LI Yu,JIAO Licheng,et al.SAR Change Detection Based on Intensity and Texture Changes[J].ISPRS Journal of Photogrammetry and Remote Sensing,2014,93(1):123-135.
[2]LI H,CELIK T,LONGBOTHAM N,et al.Gabor Feature Based Unsupervised Change Detection of Multitemporal SAR Images Based on Two-Level Clustering[J].IEEE Geoscience and Remote Sensing Letters,2015,12(12):2458-2462.
[3]邓云凯,赵凤军,王宇.星载SAR技术的发展趋势及应用浅析[J].雷达学报,2012,1(1):1-10.
[4]NAIDOO L,MATHIEU R,MAIN R,et al.The Assessment of Data Mining Algorithms for Modelling Savannah Woody Cover Using Multi-Frequency(X-,C-and L-Band)Synthetic Aperture Radar(SAR)Datasets[C]∥2014 IEEE International Geoscience and Remote Sensing Symposium,Quebec City,QC,Canada:IEEE,2014:1049-1052.
[5]BHATTACHARYA C,ROY A,NAVNEET S,et al.MicroSAR:Calibration of X-Band High Resolution FMCW Synthetic Aperture Radar(SAR)[C]∥2014 IEEE International Microwave and RF Conference,Bangalore,India:IEEE,2014:377-380.
[6]KRIEGER G,MOREIRA A,FIEDLER H,et al.TanDEM-X:A Satellite Formation for High-Resolution SAR Interferometry[J].IEEE Trans on Geoscience and Remote Sensing,2007,45(11):3317-3341.
[7]WANG Yan,DU Lan,DAI Hui.Unsupervised SAR Image Change Detection Based on SIFT Keypoints and Region Information[J].IEEE Geoscience and Remote Sensing Letters,2016,13(7):931-935.
[8]DOLLÁR P,TU Z,PERONA P,et al.Integral Channel Features[C]∥British Machine Vision Conference,London:BMVA Press,2009:1-11.
[9]YANG Bin,YAN Junjie,LEI Zhen,et al.Aggregate Channel Features for Multi-View Face Detection[C]∥2014 IEEE International Joint Conference on Biometrics,Clearwater,FL,USA:IEEE,2014:1-8.
[10]WANG Shuo,YANG Bin,LEI Zhen,et al.A Convolutional Neural Network Combined with Aggregate Channel Feature for Face Detection[C]∥6th International Conference on Wireless,Mobile and Multi-Media,Beijing:IET,2015:304-308.
[11]HE Kaiming,SUN Jian,TANG Xiaoou.Guided Image Filtering[J].IEEE Trans on Pattern Analysis and Machine Intelligence,2013,35(6):1397-1409.
[12]LEVIN A,LISCHINSKI D,WEISS Y.A Closed-Form Solution to Natural Image Matting[J].IEEE Trans on Pattern Analysis and Machine Intelligence,2008,30(2):228-242.
[13]ZOMET A,PELEG S.Multi-Sensor Super-Resolution[C]∥Sixth IEEE Workshop on Applications of Computer Vision,Orlando,FL,USA:IEEE,2002:27-31.
[14]HE Kaiming,SUN Jian,TANG Xiaoou.Single Image Haze Removal Using Dark Channel Prior[J].IEEE Trans on Pattern Analysis and Machine Intelligence,2011,33(12):2341-2353.
[15]吴涛,陈曦,牛蕾,等.非监督SAR图像变化检测研究最新进展[J].遥感信息,2013,28(1):110-118.
[16]CELIK T.Unsupervised Change Detection in Satellite Images Using Principal Component Analysis and k-Means Clustering[J].IEEE Geoscience and Remote Sensing Letters,2009,6(4):772-776.
[17]GONG Maoguo,SU Linzhi,JIA Meng,et al.Fuzzy Clustering with a Modified MRF Energy Function for Change Detection in Synthetic Aperture Radar Images[J].IEEE Trans on Fuzzy Systems,2014,22(1):98-109.
[18]万红林,汪洋,江凯,等.基于模糊贴近度和非紧凑邻域的变化检测[J].雷达科学与技术,2014,12(3):229-234.WAN Honglin,WANG Yang,JIANG Kai,et al.Fuzzy Neartude and Information in a Non-Compact Neighborhood Based Change Detection[J].Radar Science and Technology,2014,12(3):229-234.(in Chinese)
