1 引 言
合成孔径雷达(synthetic aperture radar,SAR)图像的变化检测，由于数据的内在复杂性，既要求一个可有效去除斑点噪声的预处理，又要求一个能有效处理乘性斑点噪声的分析技术，这使得相对光学图像以SAR图像为研究对象的变化检测技术较少。然而，与光学遥感相比，SAR具有的不受气候和太阳光照条件影响的特点，使之具有潜在的实际应用价值，比如，SAR能有规律地监测常出现云雨天气的地区，而光学遥感则易被中断。目前，SAR图像的变化检测技术已成功用于地震灾区的标识[1]、水灾监控[2]和水稻栽种面积的评估[3]等方面。
一般地，若检测由水灾、火灾或地震造成的突变性变化，只需事件发生前后两时相的图像，而检测如绿地覆盖率、湿地变化等逐渐发生的变化，则需一系列不同时刻获取的图像[4]。本文检测由水灾引起的突变性变化，两时相图像即可。对已配准、已校正的两时相SAR图像，无监督的变化检测技术一般基于以下3步：预处理、图像比较和阈值分割[5]。预处理的目的是去除斑点噪声，现有的自适应滤波器，如Lee或Gamma MAP滤波器等，均可用于去除斑点噪声，但都需事先设定滤波程度。图像比较是以对应位置有变化发生的可能性大小作为像素的灰度值生成差异图。由于SAR图像所含斑点噪声的乘性特点，最广泛使用的方法是比值法[6]，但该法只利用了单个像素的信息。现有基于空间信息的差异图构造方法，如：交叉互相关测度[7]、区域均值比等[8]，这种比较策略能较好地抑制因辐射等原因引起的灰度值突变的点。本文即采用区域均值比构造差异图。
差异图生成后，比较经典的方法是通过阈值分割差异图确定变化类和非变化类。一般的，基于阈值的差异图分析技术在基于像素间相互独立的假设下，利用单个像素的信息人工[9]或自动[10, 11, 12, 13]地给出决策阈值。然而，真实图像中邻域像素间具有很大相关性。为利用空间关系，文献[12]通过马尔可夫随机场(Markov random field,MRF)模型在决策阶段将其引入，而后，文献[13]利用邻域像素间的空间关系纠正广义高斯模型下EM算法的分类结果。结果表明变化检测性能有所提高，但与文献[10, 11]中的方法一样，变化检测结果仍受所选统计模型与真实分布函数间拟合度的影响。为避免分布模型的假设，文献[14]提出基于DCT特征抽取，文献[15]提出基于纹理分类器的差异图分析技术。然而，上述方法的变化检测性能均依赖于预处理，这意味着滤波程度的设定问题仍存在。为此，文献[5, 16, 17, 18]利用多尺度分析工具代替滤波器，免去了预处理操作，但文献[16, 17, 18]中的技术仍涉及分布模型的假设，而文献[5]涉及多个参数的组合优化。
为了免去预处理及克服分布模型假设的限制，本文借助静态小波变换，结合差异图的特点及交互式的分割技术，给出了一种新的SAR图像变化检测方法。利用静态小波变换将每个像素的特征设置为对噪声具有一定抗差性的由灰度值构成的矢量后，再由一种不涉及分布模型的交互式分割技术[19]结合由差异图的特点给出的“种子点”，并通过像素间的竞争(特征间的比较)将差异图分成两类，得到不同“种子点”下的变化检测结果；而后，再对不同“种子点”的变化检测结果进行决策级融合，完成变化检测。该方法无需对已配准、已校正的两时相SAR图像作降斑预处理，能产生优于文献[5]中最好方法的变化检测结果，并克服文献[10, 11, 12, 13, 14, 15, 16, 17, 18]需为类分布选择合适的统计模型的限制。对真实SAR图像试验数据的变化检测结果及比较，说明本文方法优于其他相关技术。
2 差异图的分析对已配准、已校正的两时刻SAR图像It1和It2，构造差异图DA
式中，N(p)表示位于(i,j)的像素p的邻域像素构成的集合，p′为其中的一个元素。考虑到窗口越大，噪声去除得越干净，但变化区域扩大得也越厉害，本文依据实际情况选取大小为3×3的窗口。给定差异图后，若把变化检测看成是一个对差异图二分类的问题，则利用交互式的分割技术生成变化检测结果是可行的，但存在以下两个问题：如何给出交互式分割技术所需的“种子点”及克服差异图中噪声对其分割性能的影响。对交互式分割技术所需的“种子点”，本文根据差异图的灰度分布特点，给出分别属于变化类和非变化类的子集，并将其作为分割技术所需的种子点；对差异图中残留的部分噪声，通过设置对噪声具有一定抗差性的特征矢量来减少残留的部分噪声对分割性能的影响。最后，由于不同“种子点”对应不同的变化检测结果，对不同“种子点”下得到的变化检测结果，采用投票竞争策略对其进行决策级的融合生成最终的变化检测结果。决策级融合的策略使最终的变化检测结果对“种子点”的选取不敏感的同时，也加强了对噪声的抗差性。流程图如图 1所示。
|
| 图 1 本文方法的框架图 Fig. 1 The general scheme of our proposed method |
将交互式分割技术用于分析差异图DA时，其实质为对差异图中每个像素赋予一个属性向量后，由指定的“种子点”开始，根据竞争规则δ，通过邻域N逐渐向外扩张，直到所有像素的状态子集s不再变化。其中每个像素p的属性向量由标签lp、强度θp、特征Vp、邻域系统N及竞争规则δ构成，状态子集由标签lp、强度θp、特征矢量Vp构成，即sp={lp,θp,Vp}。对给定的“种子点”，迭代结束后，所有像素的标签即构成该“种子点”下的变化检测结果。
如前所述，迭代前需解决种子点的选取和像素特征的设置问题。对分割方法所需的种子点，根据差异图的分布特点[13]给出。利用差异图的直方图hD计算出两个阈值Tn和Tc后，由该阈值确定分别属于变化类和非变化类的“种子”。若两阈值为Tn=MD(1－α)，Tc=MD(1+α)，则变化类的种子有Sc={DA(p)|DA(p)＞Tc}，非变化类的种子有Sn={DA(p)|DA(p)＜Tn}，其中MD为直方图hD的中间值，即MD=[max{DA}－min{DA}]/2；α∈(0,1)是人为设定的一个参数，表征灰度值介于以MD为中心，MD×α为邻域半径内的像素，其类别是有待确定的。
关于设置像素p的特征Vp问题，其关键在于部分残留在差异图中的斑点噪声对分割技术的性能有影响，因此需设置对噪声具有一定的抗差性的特征。受文献[5]启发，根据静态小波变换对差异图分解后再重构所得的对差异图的不同层表示，其内噪声被去除程度不同的特点设置特征Vp。利用静态小波变换获得对差异图的不同层表示后，将像素p在差异图及其各层表示内对应位置上的灰度值构成的矢量作为特征Vp，以消除斑点噪声对分割技术性能的影响。静态小波变换分解差异图时，每层都产生4个大小与原图一致的图像，其中1个表示低频信息，3个表示高频信息。鉴于斑点噪声都处在高频系数内，将分解后表示高频信息的变换系数置0，只利用低频系数进行逆静态小波变换，得到对差异图的某个分解层数下的表示。随着分解层数的增加，该层表示内噪声去除的越多，但同时丢失的细节信息也越多。考虑到噪声去除和细节信息保持间的两难，分别求得对差异图的1层和2层表示，将特征设置为由三维灰度值构成的矢量。若Dk表示对差异图的k层表示，则对位置(i,j)上的像素p，其特征矢量Vp为Vp=[DA(i,j) D1(i,j) D2(i,j)]。
初始时，给定一个α，确定出变化类wc的种子
，非变化类wn的种子
及类别不确定的像素集Sα后，将差异图中每个像素的状态子集设置为
不同的α代表选取不同的种子点，而不同的种子点产生不同的变化检测结果，这表明变化检测结果的效果与α的取值有关。为了使最终的变化检测结果较稳定，同时增强其对斑点噪声的抗差性，对不同种子点下的变化检测结果进行决策级融合。将每个α对应的变化检测结果都作为一个专家的判断，根据少数服从多数的原则，通过投票策略[20]决定最终的变化检测结果。对位置(i,j)上的像素p，其类别为获取投票数最多的一类。若以
,k=c,n表示专家认为该像素属于变化类(或非变化类)的票数，则差异图DA中每个像素的类别可按式(4)确定
第一组真实的试验数据是关于Bern城市水灾的SAR图像，由欧洲遥感2号卫星SAR传感器获得，大小为301×301像素，如图 2所示。图 2(a)、(b)分别为水灾前1999年4月的图像和水灾后1999年5月的图像，图 2(c)为变化参考图，其中真实变化目标数为1155个像素单元。要检测的是1999年5月水灾引起的变化。
|
| 图 2 第1组真实的有关Bern城市的试验数据 Fig. 2 The first real dataset to be experimented on Bern city |
第2组真实的试验数据是Ottawa地区水灾的Radarsat SAR图像，大小为290×350像素，如图 3所示。图 3(a)为正处雨季，1997年5月的图像，图 3(b)为雨季过后，1997年8月的图像，图 3(c)为变化参考图，其中变化目标数为16 049个像素单元。两幅图像已经过配准，要检测的是因为雨季来临，洪水泛滥造成的变化。
|
| 图 3 第2组真实的有关Ottawa地区的试验数据 Fig. 3 The second real dataset on Ottawa district |
虽然已知α的值域为(0,1)，且已知决策级融合可以消除不同的α对变化检测性能的影响，但融合哪些α下的变化检测结果得到最终的变化检测结果，即存在α该如何取值的问题。为简单起见，本文采取均分[0.05,0.95]的取值方式。这样，最终变化检测结果转化为受取值间隔Δα的影响。为验证本文方法的有效性，先分析参数Δα对变化检测结果性能的影响，而后，从本文方法具有的3个特点组织了以下的比较试验证明本文方法的有效性：① 与无监督方式的变化检测技术GGKIT[10]和LN-GKIT(或WR-GKIT)[11]比较；② 与上下文敏感分析技术(EM+MRF)[13]比较；③ 与无分布假设及利用上下文信息的变化检测技术(FFL-ARS)[5]比较。其中，本文方法的变化检测结果为投票策略融合α={0.05:Δα:0.95}，Δα=0.05下的变化检测结果所得。
为了更客观地说明变化检测结果的效果，给出4个定量评价检测结果的分析指标：正确检测数(right detections(RD)，变化的像素被判断为变化的像素个数)、漏检数(missed alarms(MA)，变化的像素被判为未变化的像素个数)、虚警(false alarms(FA)，未变化的像素被判为变化的像素个数，也叫错误检测数)和总错误检测数(overall error(OE)，即漏检数＋虚警)。一个好的变化检测技术在于使正确检测数尽可能地多，而总错误检测数尽可能地少。
3.3 参数分析为分析参数α的取值间隔Δα对变化检测结果的影响，以Bern试验数据为例，给出以不同取值间隔Δα={0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1}均分区间[0.05,0.95]下各α对应的变化检测结果，并在图 4给出其定量分析指标随参数Δα的变化曲线。
|
| 图 4 本文方法的变化检测结果受参数Δα的影响曲线 Fig. 4 Behaviors of quantitative analysis terms changing with the parameter Δα |
给出不同Δα下变化检测结果的定量分析指标后，得出正确检测数在90.4±1.562 2个像素单元内波动，而虚警数在289.1±5.557 1个像素单元范围内波动，这说明本文方法具有较好的变化检测精度及较稳定的性能。图 4展示的是取值间隔Δα对变化检测结果性能的影响曲线，也表明本方法的变化检测性能受取值间隔Δα的影响较小。
3.4 试验结果及分析对Bern试验数据，各差异图分析技术的变化检测结果在图 5展示，其中图 5(a)、(b)分别为文献[10, 11]给出的对Bern 试验数据最好的变化检测结果图，图 5(e)为本文方法的变化检测结果图。表 1给出了图 5中各变化检测结果的定量分析指标。图 5(a)为对数比值法差异图，假设其类分布模型符合广义高斯分布(GGKIT)产生的变化检测结果图，图 5(b)为对比值法差异图，假设其分布函数符合对数正态分布模型(LN-GKIT)得到的变化检测结果；比较GGKIT和LN-GKIT法，GGKIT法的漏检数比本文的多33个像素单元，LN-GKIT法的漏检数比本文的多175个像素单元，这一方面是因为GGKIT和LN-GKIT法均需降斑预处理，而滤波程度设置得过小漏检数较多，反之，则错误检测数较多；另一方面是因为GGKIT和LN-GKIT法均涉及统计模型的选择，对相同的试验数据，相同的差异图构造，不同的分布模型对应不同的变化检测结果，即变化检测结果的好坏取决于所选与真实的分布模型间的拟合度，若选取的分布模型不合适的话，变化检测效果会不如本文方法的。
|
| 图 5 差异图分析技术对Bern试验数据均值比差异图的变化检测结果 Fig. 5 The results on Bern dataset of different techniques |
| RD | MA | FA | OE | |
| GGKIT | 1026 | 129 | 274 | 403 |
| LN-GKIT | 884 | 271 | 65 | 333 |
| EM+MRF | 1107 | 48 | 3386 | 3434 |
| FFL-ARS | 756 | 399 | 87 | 486 |
| Ours | 1059 | 96 | 283 | 379 |
| MTEP | 929 | 226 | 88 | 314 |
图 5(c)和图 5(d)分别为上下文敏感技术EM+MRF和FFL-ARS法的变化检测结果。比较EM+MRF法与本文方法，虽然EM+MRF法利用了邻域信息，但EM+MRF法多3085个像素单元。这主要是因为：① 差异图中的斑点噪声影响了EM+MRF法的变化检测性能；② 暗含地认为差异图的类分布符合高斯模型，而这与Bern试验数据的差异图中两类的真实分布不太符合。比较FFL-ARS法与本文方法，FFL-ARS法的总错误检测数比本文方法多107个像素单元。虽然FFL-ARS法不涉及分布模型的选择且利用了邻域信息，但该方法丢弃了差异图自身所含的信息，只利用了由静态小波变换对差异图进行分解后再丢弃高频系数重构的图像数据(即对差异图的k(k≥1)层表示)信息，而无论k取值多少，细节信息难免会有所丢失。本文方法利用了差异图自身的信息，在一定程度上抑制了细节信息的丢失。比较两方法的漏检数，本文方法的漏检数比FFL-ARS方法的少303个像素单元，这也定量地说明了FFL-ARS法细节信息丢失的较多。
通过上述比较，结果显示本文方法优于上述变化检测技术。本文方法的变化检测结果拥有较多的正确检测数的同时，含有的总错误检测数也较少，这说明在总错误检测数和正确检测数这个两难的问题上达到了一个较好的平衡。究其原因，主要是因为：① 本文方法不涉及降斑预处理，避免了因滤波程度设置不当引起效果降低的可能性，且设置特征时，利用了差异图本身，在一定程度上抑制了信息的丢失；② 不涉及分布模型的选择，变化检测结果只与初始选取的“种子点”有关，而通过决策级融合后，变化检测结果已相对稳定，且进一步克服了斑点噪声对变化检测结果的影响；③ 在变化检测过程利用邻域内像素的竞争完成变化检测，将邻域信息引入提高了变化检测精度。图 5(f)为给定差异图下能获取的最优的变化检测结果，即MTEP技术的变化检测结果，与之相比进一步说明本文方法效果的有效性。
对Ottawa试验数据，文献[10]中GGKIT法的变化检测结果较好，即对数比值法构造的差异图，假设其类条件分布符合广义高斯分布更合适，如图 6(a)所示；对比值法构造的差异图，假设其类分布符合“Weibull-Ratio”模型更合适，即WR-GKIT法的变化检测结果较好，如图 6(b)所示。图 6(c)为上下文敏感技术EM+MRF法的变化检测结果图，而图 6(d)为上下文敏感和无分布假设的FFL-ARS法的变化检测图；本文方法及MTEP技术的变化检测结果分别在图 6(e)和图 6(f)展示。表 2给出了各变化检测结果的定量分析结果。
|
| 图 6 差异图分析技术对Ottawa试验数据均值比差异图的变化检测结果 Fig. 6 The results on Ottawa dataset of different techniques |
| RD | MA | FA | OE | |
| GGKIT | 13 630 | 2 419 | 288 | 2 707 |
| WR-GKIT | 15 530 | 519 | 2 074 | 2 593 |
| EM+MRF | 15 535 | 514 | 1 250 | 1 764 |
| FFL-ARS | 13 789 | 2 260 | 641 | 3 201 |
| Ours | 15 686 | 363 | 836 | 1 199 |
| MTEP | 15 201 | 848 | 826 | 1 674 |
分析表 2的数据可知：与Bern试验数据相同的是，对Ottawa试验数据，本文方法相比于其他技术，变化检测效果有所改善。对Ottawa试验数据，本文方法变化检测结果的总错误检测数为1199个像素单元，比GGKIT法少1508个像素单元，比WR-GKIT法的少1394个像素单元；比上下文敏感的EM+MRF法少547个像素单元；比同时具有上下文敏感和无分布假设的FFL-ARS法少2002个像素单元。与MTEP技术的变化检测结果相比，可进一步说明本文方法的有效性。而与Bern试验数据不同的是，对Ottawa试验数据，文献[11]中WR-GKIT法的变化检测效果较好，而对Bern试验数据，相同方法构造的差异图，LN-GKIT法的变化检测效果较好，这进一步说明本文提出的变化检测技术，不涉及分布模型假设的变化检测技术的有效性。
4 结 论本文利用差异图的特点给出种子点后，结合种子点和一种基于邻域像素竞争的交互式分割技术及决策级融合策略，给出了一种无监督、上下文敏感及无分布假设SAR图像变化检测技术。该方法无需要求对SAR图像作降斑预处理，且性能优于其他相关的变化检测技术，但当处理的图像所含噪声严重时，可以利用静态小波变换对差异图作多次分解后再重构，将特征设置为由灰度值构成的高维矢量，以增强对斑点噪声的抗差性；也可以先作预处理去除部分斑点噪声再由本文技术产生变化检测结果，其中，滤波程度不够强、仍有残留的斑点噪声对变化检测结果有影响的情况可以不考虑，因为本文方法可以较好地处理含噪声不严重的图像。
分割技术中邻域像素间竞争实为特征间的相似性比较，本文以灰度为特征、欧式范数为测度。实际上，使用者可以考虑任何具有辨别力的特征，诸如纹理、方差、局部直方图等以及相应的相似性度量测度，如自适应的距离测度[22]等。
| [1] | RIDD M, LIU J. A Comparison of Four Algorithms for Change Detection in an Urban Environment [J]. Remote Sensing of Environment,1998, 63(2): 95-100. |
| [2] | MAS J F. Monitoring Land-cover Changes: a Comparison of Change Detection Techniques [J]. International Journal of Remote Sensing, 1999, 20(1): 139-152. |
| [3] | WEYDAHL D J. Detecting Seasonal Changes at Nordic Latitudes Using ERS-1 SAR Images [C]//Proceedings of IEEE International Geoscience and Remote Sensing Symposium, IGARSS’ 93. Tokyo: IEEE Computer Society, 1993: 1462-1464. |
| [4] | COPPIN P, JONCKHEERE I, NACKAERTS K, et al. Digital Change Detection Methods in Ecosystem Monitoring: a Review [J]. International Journal of Remote Sensing, 2004, 25(9): 1565-1596. |
| [5] | BOVOLO F, BRUZZONE L. A Detail-preserving Scale-driven Approach to Change Detection in Multitemporal SAR Images [J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(12): 2963-2972. |
| [6] | RIGNOT E J M, VAN ZYL J J. Change Detection Techniques form ERS-1 SAR Data [J]. IEEE Transactions on Geoscience and Remote Sensing,1993, 31(4): 896-906. |
| [7] | ONANA V P, TROUVE E, MAURIS G, et al. Change Detection in Urban Context with Multitemporal ERS-SAR Images by Using Data Fusion Approach[C]//Proceedings of IEEE International Geoscience and Remote Sensing Symposium, IGARSS’ 03. Toulouse: IEEE Computer Society, 2003: 3650-3652. |
| [8] | WU F, WANG C, Zhang H, et al. Change Detection and Analysis with Radarsat-1 SAR Image[C]//Proceedings of IEEE International Geoscience and Remote Sensing Symposium, IGARSS’07. Barcelona: IEEE Computer Society, 2007: 2601-2604. |
| [9] | DEKKER R J. Speckle Filtering in Satellite SAR Change Detection Imagery [J]. International Journal of Remote Sensing, 1998, 19(6): 1133-1146. |
| [10] | 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 Transactions on Geoscience and Remote Sensing, 2005, 43(4): 874-887. |
| [11] | MOSER G, SERPICO S B. Generalized Minimum-error Thresholding for Unsupervised Change Detection from SAR Amplitude Imagery [J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(10):2972-2983. |
| [12] | BRUZZONE L, PRIETO D F. Automatic Analysis of the Difference Image for Unsupervised Change Detection[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(3): 1171-1182. |
| [13] | 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. |
| [14] | YU H T, LIN A P, CHAI G L, et al.A DCT-based Change Detection Method for Multi-temporal SAR Images [C]//Proceedings of IEEE International Geoscience and Remote Sensing Symposimu, IGARSS’08, Boston: IEEE Computer Society, 2008: 166-169. |
| [15] | ZHANG S Q, LU H Q. Learning Texture Classifier for Flooded Region Detection in SAR Images [C]//Proceedings of the International Conference on Computer Graphics, Imaging and Visualization, CGIV’05. Washington DC: IEEE Computer Society, 2005: 93-98. |
| [16] | CELIK T, MA K K. Unsupervised Change Detection for Satellite Images Using Dual-tree Complex Wavelet Transform[J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(3): 1199-1210. |
| [17] | CELIK T. A Bayesian Approach to Unsupervised Multiscale Change Detection in Synthetic Aperture Radar Images[J]. Signal Processing, 2010, 90(5):1147-1471. |
| [18] | 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.) |
| [19] | VEZHNEVETS V, KONOUCHINE V. “GrowCut”—Interactive Multi-label N-D Image Segmentation by Cellular Automata[C]//Proceedings of Graphicon. ]Novosibirsk:[s.n., 2005. |
| [20] | KITTLER J, HATEF M, DUIN R P W, et al. On Combining Classifier [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1998, 20(3): 226-239. |
| [21] | KITTLER J, ILLINGWORTH J. Minimum Error Thresholding[J]. Journal of Pattern Recognition, 1986, 19(1): 41-47. |
| [22] | DAVIS J V, KULIS B, JAIN P, et al. Information—Theoretic Metric Learning[C]//Proceedings of the 24th International Conference on Machine learning,ICML’07. New York: ACM, 2007:209-216. |