1 引言
变化检测是分析同一地区不同时刻获取的图像以识别某一目标或某种现象状态差异的过程[1]，并已成为开展环境监控[2]，土地利用/覆盖动态研究[3]等方面的关键技术。然而，相比于光学图像，以合成孔径雷达(synthetic aperture radar,SAR)图像为研究对象的变化检测技术较少。这主要是因为SAR图像数据的内在复杂性，既要求有效地去除斑点噪声的预处理，又要求能有效地处理乘性斑点噪声的分析技术。但SAR具有不受气候和光照条件影响的特点，这使得SAR图像具有潜在广泛的应用价值。
一般的，无监督的SAR图像变化检测技术由预处理，图像比较和阈值分割3步构成[4]。预处理的目的是去除斑点噪声，以加大差异图中两类的可分性。现有的滤波器，如Lee或Gamma Map滤波器均可用于降斑。图像比较是为生成可后续分析的差异图。差异图中两类的对比度越高，各类的类内一致性越强，越有利于后续分析。对SAR图像而言，由于其所含斑点噪声的乘性特点，比算子比差分算子更适合[5]，如比对数法(先取对数再作比)、对数比法(先作比再取对数)、区域均值比法等。差异图生成后，最简单易行的想法是阈值处理差异图生成变化检测结果，如手动的人工尝试法[6]、自动的分析技术[4, 7, 8, 9]等。
基于阈值的变化检测技术一般需要降斑预处理，以使变化检测效果更理想[7]，而预处理需人为设定滤波程度。为了避免该操作，文献[10, 11, 12]用多尺度分析工具，如平稳小波变换或对偶树复小波变换等替代滤波器，获取多个具有不同的噪声去除和细节保持平衡度的差异图表示。此外，基于阈值的变化检测技术一般都忽略邻域像素间的关系，而真实图像中邻域像素间具有很大的相关性。为此，文献[13, 14]通过马尔可夫随机场(Markov random field,MRF)模型或基于MRF的融合策略在决策阶段将空间关系引入。结果表明，空间关系的引入有助于提高变化检测性能，但文献[13]方法涉及分布模型的假设，文献[14]方法的变化检测性能依赖于阈值技术。而后，文献[15,16,17,18]利用马尔可夫链或结合模糊理论与马尔可夫链在分析过程中考虑空间关系。但变化检测结果仍是在像素层面上决策生成，这易导致变化检测结果中有孤立的像素点，连通区域内有孔洞及锯齿状边缘等缺陷[19]。
为克服滤波程度的设定及弥补在像素层面上决策生成的变化检测结果的缺陷，笔者提出一种基于抽取并处理感兴趣区域(region of interest,ROI)，在区域层面上决策生成变化检测结果的技术。抽取感兴趣区域的策略已成功地应用于提高阈值技术分割自然或医学图像的效果[20]。本文借用该思想并用于SAR图像变化检测，本文方法以使变化检测结果在区域层面上生成，定量性能有所提高的同时，视觉效果也有所改善。本文方法的关键在于无监督地求取合适的、用于抽取感兴趣区域的标签和在区域层面上处理感兴趣区域两方面。首先，本文在平稳小波变换(stationary wavelet transform,SWT)的帮助下，结合模糊C-均值(fuzzy C-means,FCM)聚类算法[21]和尺度间融合策略获取该标签。而后，为提高变化检测性能并在区域层面上生成变化检测结果，本文依据标签抽取感兴趣区域后，搜索感兴趣区域内所有的连通区域并把每个连通区域看成一个整体，由阈值技术处理生成最终的变化检测结果。对真实SAR图像数据的变化检测结果及比较，表明本文方法变化检测结果在定量性能和视觉效果都有所改善。
2 基于抽取感兴趣区域的变化检测对已配准、已校正的同一地区获取的两时相SAR图像I1和I2，构造出比对数法差异图Drl和区域均值比法差异图Dwr后，本文在基于抽取并处理感兴趣区域的策略上完成变化检测。该方法由两步构成：① 用SWT对Drl分解后，由FCM聚类算法对各层的逼近系数先分析后融合得到结果L，为消除因变换而产生的平移，用L′替代L，其中L′为FCM聚类算法处理依据结果L在Dwr内抽取的ROI得到的结果；② 依据标签L′，在Dwr内抽取感兴趣区域ROI′并由阈值技术(TH)在区域层面上处理生成最终的变化检测结果(change detection map,CDM)。具体的算法流程如图 1所示。
|
| 图 1 本文方法框架图 Fig. 1 The general scheme of the proposed method |
为得到抽取感兴趣区域所需的标签L′，先利用SWT对Drl进行多尺度分解并用FCM聚类算法对分解后每层的逼近系数分类，而后进行尺度间的融合。SWT分解Drl时，每层都产生1个低频子带和3个高频子带，且各层的低频和高频子带，大小与原图一致。鉴于分解后，斑点噪声一般处在高频系数内，实际工作中只利用低频子带系数。然而这样，细节信息的丢失会随着分解层数的增加而增多，因而本文采用尺度间“并”的融合策略，以使低分解层数下逼近系数内的细节信息能被利用。若用liDrl表示FCM聚类算法对第i层(条)低频子带的分类结果，则标签L可表示为
式中，S表示分解层数。
平稳小波变换本具有平移不变性，但标签L是分析变换域内的低频子带系数得到的。因而，相对于真实的变化区域，L指示的变化区域存有偏移。为消除此偏移，用FCM聚类算法对依据标签L抽取的感兴趣区域分类所得的结果L′替代L。为确保抽取的感兴趣区域包含几乎所有的真实变化的像素，先对L进行以3×3的全1矩阵为结构元素的膨胀后，再依据该标签在区域均值比值法差异图Dwr中抽取对应位置上的标签值为1的像素构成感兴趣区域，即
FCM聚类算法分类ROI所得结果即为L′。该结果相比于L，偏移被消除了，同时散落的杂点也有所减少。这是因为ROI，相对于差异图，部分属于非变化类的像素对结果的影响被排除。正如文献[19]所说，任何类型的先验信息加入图像处理的应用过程都有利于提高性能。然而，结果仍在像素层面上产生，其内存在如噪声般散落的杂点。
2.2 在区域层面上处理感兴趣区域生成变化检测结果为清除杂点，依据L′抽取感兴趣区域并搜索其内所有的连通区域后，再以连通区域为处理单元，使阈值技术在区域层面上生成变化检测结果。其中处理单元的特征设置为对应区域内所有像素灰度的平均值。但这样处理时存在一个问题，即L′中每个杂点都会被看做是一个连通区域，且具有较高的灰度特征，而这易导致杂点仍被判断为属于变化类。为此，引入一个参数：最小面积阈值(minimum area,minA)。对面积小于minA的区域，从非感兴趣区域内拉入与该区域相邻的像素以扩大其面积到最小面积阈值。由于被扩充的像素不在感兴趣区域内，通常具有不大于感兴趣区域内像素的灰度值，因而对应于小面积区域的数据点的灰度值会有所下降。这意味着数据点间的对比度被增强，因而更有利于清除小面积区域而提高变化检测性能。
依据标签L′抽取感兴趣区域
最终的变化检测结果由阈值技术在区域层面上处理ROI′得到。具体步骤如下：
(1) 根据标签L′搜寻ROI′内所有的连通区域，而后把每个连通区域看做为一个数据点并计算特征。若假设感兴趣区域ROI′内存在M片连通区域，并用{c1,c2,…，ci，…,cM}表示感兴趣区域ROI′内所含的连通区域，则对应的数据点的特征值Fi为
式中，c′i表示与ci对应的被扩展的连通区域；minA为预先设置的参数值。
(2) 将M个数据点分成两类。此处，采用从聚类分析的角度提出阈值标准的阈值技术[22]。该技术的核心思想为：每迭代一次，合并直方图上相邻且相似度最大的两类，直到只剩下想得到的类别数为止。初始时，若M个数据点的特征值均不等，则初始有M个类别模式；若有2个数据点的特征值相同，则有M－1个类别模式，依此类推。假设初始有K个类别模式，为将表征M片连通区域的数据点归类成2类(变化类和非变化类)，需重复下列两步K－2次：
① 按式(5)
计算在特征值上相邻的两个类别Cs和Ct间的相似性。该值越大，表示相似性越小。式中，σI表示类间差异，越大越好；σA表示类内方差，越小越好。
Zk是由类别Ck中数据点的特征值构成的真子集；PCk=∑z∈Zkhz为类Ck的概率;mCk=
为类Ck的均值,k=s,t；
类Cs和Ct合并后的均值。其中，直方图h按如下方式计算得到：若假设特征值μ的数据点有j个，则其概率
表示一个集合的势。
② 搜索使距离最小的两类，并将其合并，而后再重新分配相应的类别索引。
(3) 生成最终的变化检测结果CDM。迭代结束后，灰度值较低类中的最大灰度值为最优阈值t*，则变化检测结果为
即特征值大于该阈值t*的数据点对应的连通区域组成的区域即为本文方法检测出的变化区域。 3 参数分析及试验结果对比 3.1 试验数据描述
第1组真实的试验数据是关于Bern城市水灾的SAR图像，由欧空局ERS-2获得，大小为301×301。图 2(a)和图 2(b)分别为水灾前1999年4月的图像和水灾后1999年5月的图像，图 2(c)为其变化参考图，其中白色区域表示由1999年5月水灾引起的变化区域，面积为1155个像素单元。
|
| 图 2 Bern试验数据 Fig. 2 Bern dataset |
第2组真实的试验数据是有关Ottawa地区水灾的由加空局获取的RADARSAT-1 SAR图像(如图 3所示)，大小为290×350。图 3(a)为雨季中，1997年5月图像,图 3(b)为雨季过后，1997年8月图像，图 3(c)为变化参考图，其中白色区域表示洪水泛滥前后引起的变化，面积为16 049个像素单元。
|
| 图 3 Ottawa试验数据 Fig. 3 Ottawa dataset |
因本文方法涉及参数最小面积阈值minA和分解层数S，特组织两组试验分别说明本文技术对两参数的敏感度。为证明本文方法变化检测效果的有效性，将给出均为上下文敏感方法EM+MRF[13]，Ths+MRF[14]和FFL-ARS[10]的变化检测结果图，并作比较。本文方法的参数设置为S=2,minA=40。
为了更直观地分析各变化检测结果的效果，给出四个定量评价检测结果的分析指标：漏检数(missed alarms,MA，即变化的像素被判为未变化的像素个数)、虚警(false alarms,FA，即未变化的像素被判为变化的像素个数，也叫错误检测数)，总错误检测数(overall errors,OE，即漏检数＋虚警)和总错误检测率(总错误检测数与所处理的图像像素个数之比)和Kappa系数
式中，N为差异图中像素的个数；Fc为变化类像素被正确检测的个数；Tc为非变化类像素被正确检测的个数。若一个变化检测结果的漏检数越少，而总错误检测数也越少，Kappa系数越大，表明该技术的变化检测性能越好。 3.3 参数分析
第1组试验是为说明本文技术对参数minA的鲁棒性。此处固定参数S=2，分析参数minA对变化检测结果的影响。图 4给出两组试验数据变化检测结果的定量分析指标受参数minA影响的曲线。
|
| 图 4 参数minA对两组试验数据变化检测结果的影响 Fig. 4 Behaviors of quantitative analysis terms of two datasets changing with the parameter minA |
图 4(a)为Bern试验数据的变化检测结果的定量分析指标随参数minA的变化曲线。很明显，当minA＞12后，变化检测结果的定量分析指标基本保持不变。图 4(b)为Ottawa试验数据变化检测结果的定量分析指标随参数minA的变化曲线，当minA＞40后，变化检测效果基本不变。这说明当minA大于某一阈值后，本文技术受其影响很小。鉴于该参数是为去除小面积的杂点而引入的，一般将minA设置为小于100的数值即可。
第2组试验是为说明本文技术受参数分解层数S的影响。此处以Bern试验数据为例，固定minA=40，给出S取不同值时，两组试验数据的变化检测结果图及定量分析。表 1给出了分解层数分别从1到6下Bern试验数据变化检测结果的定量分析指标。
|
| S=1 | 83 | 236 | 319 | 0.868 2 |
| S=2 | 82 | 251 | 333 | 0.863 3 |
| S=3 | 78 | 260 | 338 | 0.861 9 |
| S=4 | 140 | 241 | 381 | 0.839 2 |
| S=5 | 135 | 256 | 391 | 0.836 3 |
| S=6 | 128 | 292 | 420 | 0.827 2 |
分析表 1给出的定量分析数据可知：随着S的增大，变化检测结果的总错误检测数呈增加趋势，Kappa系数呈下降趋势。当S较小时，由平稳小波变换引起的偏移较小，融合所得标签指示抽取的感兴趣区域ROI内含有较少的非变化类像素。经FCM聚类算法处理后，所得的感兴趣区域ROI′内亦含有较少的非变化类像素。这时，阈值技术能给出一个较好的变化检测结果。当S较大时，由平稳小波变换引起的偏移加大，融合后感兴趣区域ROI内存在较多的非变化类像素。这导致FCM聚类算法处理该感兴趣区域不能给出一个较理想的用于抽取感兴趣区域ROI′的标签，故最终的变化检测效果不甚理想。对本文测试的试验数据，S≤3都能给出较满意的结果。
3.4 结果对比及分析对Bern试验数据，图 5给出了不同变化检测技术的变化检测结果。图 5(a)为上下文敏感技术EM+MRF分析Drl的变化检测图；FFL-ARS法的变化检测结果图为图 5(b)；图 5(c)为Ths+MRF法的变化检测结果图；图 5(d)显示了本文方法的变化检测结果。表 2给出各变化检测结果的定量分析指标。
|
| 图 5 各技术对Bern试验数据的变化检测结果 Fig. 5 Results of Bern dataset obtained by different techniques |
| MA | FA | OE | Kappa | |
| EM+MRF | 242 | 105 | 347 | 0.837 6 |
| Ths+MRF | 2 | 13 619 | 13 621 | 0.123 6 |
| FFL-ARS | 399 | 87 | 486 | 0.752 9 |
| 本文方法 | 82 | 251 | 333 | 0.863 3 |
分析表 2的数据可看出：本文方法的变化检测效果最好。与EM+MRF法比较，本文方法的Kappa系数高0.025 6，这是因为本文方法避免了EM+MRF法中假设分布模型与真实情况不太符合的情况，且通过利用无固定形状和大小的邻域信息，弥补了MRF模型只能利用固定形状和大小的邻域内信息而带来的不足。与Ths+MRF法相比，效果有了较大的改善，这是因为本文方法克服了性能依赖于阈值技术变化检测结果的限制。相比于FFL-ARS法，本文方法的Kappa系数增加了0.110 4，这说明本文方法采取的“并”融合策略的有效性。正是因为此，本文方法的细节信息相比于FFL-ARS方法保持得要好，体现在漏检数比FFL-ARS方法的少317个像素单元。
上述从定量的角度说明了本文方法的有效性。事实上，本文方法的视觉效果也有较好的改善。比较图 5(d)与图 5(a)~图 5(c)中的任意一个，不难发现图 5(a)~图 5(c)显示的变化检测结果存在较多散落的杂点、连通区域内也有孔洞或边缘不光滑等现象。这是因为其他3种方法虽然都考虑了邻域信息，但仍在像素层面上决策生成变化检测结果，而本文方法借助于抽取感兴趣区域操作实现了区域层面上的决策。
分析Ottawa试验数据的变化检测结果，同样地能说明本文方法的有效性。图 6展示了不同的差异图分析技术对Ottawa试验数据的变化检测效果图。图 6(a)为EM+MRF方法处理Dwr的变化检测结果图；图 6(b)~图 6(c)分别为FFL-ARS法和Ths+MRF法的变化检测结果图。本文的变化检测结果在图 6(d)中显示。表 3给出了图 6所列各变化检测结果的定量分析指标。
|
| 图 6 各技术对Ottawa试验数据的变化检测结果 Fig. 6 Results of Ottawa dataset obtained by different techniques |
| MA | FA | OE | Kappa | |
| EM+MRF | 513 | 1293 | 1781 | 0.934 5 |
| Ths+MRF | 308 | 2418 | 2726 | 0.904 2 |
| FFL-ARS | 2260 | 941 | 3201 | 0.877 5 |
| 本文方法 | 775 | 630 | 1385 | 0.947 8 |
分析表 3的数据可得出与Bern试验数据一样的结论，即对Ottawa试验数据，本文方法的变化检测效果仍是最好的。一方面，本文方法的Kappa系数为0.947 8，比EM+MRF法的高0.013 3，比Ths+MRF法的高0.043 6，比FFL-ARS法的高0.070 3。这表明本文方法在提高定量性能方面的有效性。另一方面，比较各变化检测效果图的视觉效果，同样发现如图 6(a)~图 6(c)所示的、在像素层面上决策产生的变化检测结果，其内有噪声般的杂点或连通区域内出现孔洞或边缘呈锯齿状等现象出现，而本文方法的变化检测结果不存在类似的现象。这说明本文方法能给出具有较好的视觉效果的变化检测结果。总而言之，对两组真实SAR图像试验数据的变化检测结果表明，本文方法相比于其他相关的变化检测技术具有较好的定量性能和视觉效果。
3.5 时间复杂度比较本文所有的试验在主机为32字节的Windows XP系统，2GHz,Inter(R),Core(TM)2 CPU,2 GB的内存上安装的Matlab平台上实现。表 4给出了4种方法分别获取Bern和Ottawa试验数据变化检测结果所需的时间。
从表 4的数据显示，本文提出的方法，相比于其他的方法所需的运行时间最少。对Bern试验数据，本文方法需15.0 s，与EM+MRF法的相当，但比Ths+MRF法和FFL-ARS法的少很多。对Ottawa试验数据，本文方法需15.9 s，同样与EM+MRF法的相当，但比Ths+MRF法和FFL-ARS法的少很多。这表明本文方法具有一定的实用性。
4 结论本文给出了一种基于抽取并处理感兴趣区域的SAR图像变化检测方法。该方法考虑了来自无固定形状和大小的邻域内信息，并在区域层面上决策产生变化检测结果。试验表明本文方法相比于相关技术，变化检测结果在定量性能和视觉效果方面都有所提高。其中用抽取和处理感兴趣区域替代MRF模型的操作使本文方法能提高变化检测性能的同时，弥补了MRF模型的内在缺陷和在像素层面上生成的变化检测结果所含的缺陷。但在采用感兴趣区域策略时，必须确保所抽取的感兴趣区域包含绝大部分的变化类像素和较小部分的非变化类像素。否则，最终变化检测结果中的漏检数会增多，效果反而会有所降低。
本文方法较适合用于检测突变型变化，如水灾、火灾等引起的，只需两时相图像即可检测的变化。若将该方法用于渐变型变化检测，如绿地覆盖率、森林退化等，效果会不太理想。
| [1] | LU D,MAUSEL P,BRONADIZIO E,et al.Change Detection Techniques[J].International Journal of Remote Sensing,2004,25(12):2365-2407. |
| [2] | RIDD M,LIU J.A Comparison of Four Algorithms for Change Detection in an Urban Environment[J].Remote Sens Environ,1998,63:95-100. |
| [3] | MAS J F.Monitoring Land-cover Changes:a Comparison of Change Detection Techniques[J].International Journal of Remote Sensing,1999,20(1):139-152. |
| [4] | BAZI Y,BRUZZONE L,MELGANI F.An Unsupervised Approach Based on the Generalized Gaussian Model to Automatic Change Detection in Multitemporal SAR Images[J].IEEE Transaction Geosciences and Remote Sensing,2005,43(4):874-887. |
| [5] | RIGNOT E J M,VAN ZYL J J.Change Detection Techniques form ERS-1 SAR Data[J].IEEE Transaction on Geosciences and Remote Sensing,1993,31(4):896-906. |
| [6] | DEKKER R J.Speckle Filtering in Satellite SAR Change Detection Imagery[J].International Journal of Remote Sensing,1998,19(6):1133-1146. |
| [7] | MOSER G,SERPICO S B.Generalized Minimum-error Thresholding for Unsupervised Change Detection from SAR Amplitude Imagery[J].IEEE Transaction on Geosciences and Remote Sensing,2006,44(10):2972-2983. |
| [8] | COULON M,TOURNERET A Y.Bayesian Change Detection for Multi-temporal SAR Images[C]//Proceedings of IEEE International Conference on Acoustics,Speech,and Signal Processing.Orlando:IEEE,2002:1285-1288,2002. |
| [9] | BAZI Y,BRUZZONE L,MELGANI F.Change Detection in Multitemporal SAR Images Based on Generalized Gaussian Distribution and EM Algorithm[C]//Proceedings of the SPIE Conference on Image and Signal Processing for Remote Sensing X.Gran Canaria:SPIE,2004:364-375. |
| [10] | BOVOLO F,BRUZZONE L.A Detail-preserving Scale-driven Approach to Change Detection in Multitemporal SAR Images[J].IEEE Transaction on Geosciences and Remote Sensing,2005,43(12):2963-2972. |
| [11] | HUANG Shiqi,LIU Daizhi,HU Mingxing,et al,Multi-temporal SAR Images Change Detection Technique Based on Wavelet Transform[J].Acta Geodaetica et Cartographica Sinica,2010,39(2):180-186.(黄世奇,刘代志,胡明星,等.基于小波变换的多时相SAR图像变化检测技术[J].测绘学报,2010,39(2):180-186.) |
| [12] | CELIK T.Multiscale Change Detection in Multitemporal Satellite Images[J].IEEE Geoscience and Remote Sensing Letters,2009,6(4):820-824. |
| [13] | CELIK T,MA K K.Unsupervised Change Detection for Satellite Images Using Dual-tree Complex Wavelet Transform[J].IEEE Trans Geosci Remote Sens,2010,48(3):1199-1210. |
| [14] | BAZI Y,MELGANI F,BRUZZONE L,et al.A Genetic Expectation-maximization Method for Unsupervised Change Detection in Multitemporal SAR Imagery[J].International Journal of Remote Sensing,2009,30(24):6591-6610. |
| [15] | BRUZZONE L,PRIETO D F.Automatic Analysis of the Difference Image for Unsupervised Change Detection[J].IEEE Transaction on Geosciences and Remote Sensing,2000,38(3):1171-1182. |
| [16] | MELGANI F,BAZI Y.Markovian Fusion Approach to Robust Unsupervised Change Detection in Remotely Sensed Imagery[J].IEEE Transaction on Geosciences and Remote Sensing,2006,3(4):457-461. |
| [17] | BOUYAHIA Z,BENYOUSSEF L,DERRODE S.Change Detection in Synthetic Aperture Radar Images with a Sliding Hidden Markov Chain Model[J].Journal of Applied Remote Sensing,2008,2(1):1-13. |
| [18] | CARINCOTTE C,DERRODE S,BOURENNANE S.Unsupervised Change Detection on SAR Images Using Fuzzy Hidden Markov Chains[J].IEEE Transaction on Geosciences and Remote Sensing,2006,44(2):432-441. |
| [19] | RADKE R J.Image Change Detection Algorithms:A Systematic Survey[J].IEEE Trans on Image Process,2005,14(3):294-307. |
| [20] | HU Q,HOU Z J,NOWINSKI L.Supervised Range-constrained Thresholding[J].IEEE Trans on Image Process,2006,15(1):228-240. |
| [21] | BEZDEK J C.Pattern Recognition with Fuzzy Objective Function Algorithms[M].New York:Plenum Pub Corp,1981. |
| [22] | ARIFIN A Z,ASANO A.Image Segmentation by Histogram Thresholding Using Hierarchical Cluster Analysis[J].Pattern Recogn Lett,2006,27(13):1515-1521. |