随着人口数量的快速增长，土地资源变得日益重要，准确实时地掌握土地利用状况，可更为有效地保护和合理地开发土地资源.遥感技术由于可选用不同波段和遥感器，快速获取大范围的地球资源信息，且不受地面条件限制，具有人工实地测量和航空摄影测量所无法比拟的优点，成为了土地利用、变化检测的重要技术手段，在地物地层分析、灾情估计、道路交通分析、森林采伐等领域有着广泛的应用(范开国等，2012).变化检测实质上就是提取不同时相遥感图像的变化信息.目前多时相遥感图像的变化检测算法主要分为基于代数运算和基于图像变换两类.基于代数运算的变化检测算法包括图像差值法(范海生等，2001)、比值法(Durucan and Ebrahimi，2001)、图像回归分析(Jha and Unni，1994)等.差值法和比值法比较直观，变化检测速度快，但算法过于简单，不能定性地描述目标区的变化，很难考虑到光照、角度、大气条件等因素的影响；而图像回归分析法一定程度上减小了这些因素所带来的不利影响，但检测时需要选择准确的回归方程和合适的波段，所以实际应用精度不高.基于图像变换的变化检测算法大都利用主分量分析(PCA-Principal component analysis)(Qiu et al.，2003；张辉和王建国，2008)或独立分量分析(Hyvarinen et al.，2004；钟家强和王润生，2006)等实现.PCA是建立在二阶统计特性基础上的多维正交线性变换，变换得到的图像主分量包含了原始遥感图像的大部分背景信息，而第二主分量包含了大部分的变化信息.但只有当信号的统计分布符合高斯分布时，信号间的相关信息才能完全被消除.而遥感图像的地物光谱特性并不满足高斯分布，因此PCA不能完全消除各分量间的高阶相关信息，这会影响变化检测的精确度.ICA作为近年来提出的一种信源盲分离技术，不仅可以去除信源之间的二阶相关信息，而且能够消除高阶相关信息(谢德光等，2001；Berg et al.，2005).因此，在遥感图像变化检测中ICA能够比PCA获得更好的检测效果.

基于ICA的变化检测算法对图像背景比较敏感，在检测变化区域时，因辐射差异易将背景变化作为目标区域的变化，导致变化检测不够精确.因此，小波变换被应用于变化检测(黄世奇等，2010；黎萍，2011)，通过小波变换对图像数据进行分解，降低了图像中的噪声，同时保护了图像中的边缘和纹理细节结构，由此提高检测的精确度和运算效率.然而，小波变换不是图像的最优稀疏表示，缺乏方向性，仅仅限于水平、垂直和对角线三个方向，在一定程度上损失了图像的有效信息，无法准确地反映图像的边缘轮廓和纹理细节.而Contourlet变换作为小波变换的一种新扩展，具有多分辨率、局部性、方向性和各向异性等特性，它对图像中边缘轮廓和纹理细节的表达具有独特优势.因此，若将Contourlet变换用于变化检测可望改善检测精确度. 另一方面，现有的ICA算法存在计算量大、收敛速度慢等缺点，对其改进后，相比现有的ICA算法可以提高收敛的速度，满足实时处理要求，且无需调整学习速率因子等其他动态参数，进而增强可靠性.

基于上述分析，本文提出了一种基于Contourlet变换和ICA算法的多时相遥感图像变化检测算法.首先通过Contourlet变换对多时相遥感图像数据进行多尺度分解，利用它的特性更好地实现图像块的优化分割，以减少ICA算法的复杂程度；然后再利用改进的ICA算法对分解后的数据进行处理，提高变化检测的准确性和运算效率.文中给出了实验结果，并与现有的基于PCA、基于ICA、基于小波变换与ICA三种变化检测算法进行了比较，分别依据主观视觉和变化像元错判误差、非变化像元误判误差、总体精度和运行时间等客观定量指标，对上述算法进行了评价，显示了本文算法的优越性.

独立分量分析是近年来由盲源信号分离技术发展起来的一种多维信号处理算法(余先川等，2010；Khaparde，2012).其基本原理是依据等独立性度量准则建立目标函数，使分离出的独立分量最大程度地逼近各源信号.设一组观测信号X=[x₁，x₂，…，x_m]^T，相互独立的源信号S=[s₁，s₂，…，s_n]^T，则第i个观测信号x_i由n个独立源信号线性混合形成：

图 1

Fig. 1

图 1 ICA原理框图 Fig. 1 The principle diagram of ICA

ICA算法需要做如下假设(Marchesi and Bruzzone，2009)：源信号S为相互独立的平稳随机过程；观测信号的数目不少于源信号的数目；源信号中至多有一个服从高斯分布.ICA实际上是一个优化问题，因为没有唯一解，只能在某一衡量独立性的判据下寻求其近似解，使Y中各分量尽可能互相独立.

Contourlet变换(白蕊等，2009；Li et al.，2012)是在继承小波多尺度分析思想基础上提出的一种多方向多尺度分析算法.它通过LP(Laplacian Pyramid)变换对图像进行多尺度分解，以捕获奇异点；然后由方向滤波器组(DFB)将分布在同方向上的奇异点合成为一个系数.Contourlet变换具有多分辨率、局部定位、多方向性、近邻界采样和各向异性等性质，能更好地捕获二维图像中的线和面奇异，用更少的系数捕捉图像中的边缘轮廓和方向性纹理信息.

对于多时相遥感图像来说，每个像元由与之对应的多种地物信息混合而成.其中每个信号源(各种物质的光谱信息)彼此独立，这些信号源有些是静止的，有些是变化的.可以认为静止的信号源构成了背景图像，变化的信号源则构成变化图像.在ICA模型里，前后时相的遥感图像构成的混合信号由背景图像和变化图像两个独立信号源混合而成，分离的具体过程如下：

(1)首先通过按行扫描的方式将每幅遥感图像转化为一维信号x_i(i=1，2，…，n)，n为遥感图像的数目，由此构成观测信号X=[x₁，x₂，…，x_n]^T.

(2)解混过程分两步实现(Ma and Wang，2010)：球化过程(W)和正交变换(U)，如图 2所示.球化过程是对观测信号X进行白化预处理，使球化后向量z=[z₁，z₂，…，z_n]^T，其方差为1，期望值为0，以降低数据计算复杂度，减少噪声；而正交变换U可使得输出数据Y=[y₁，y₂，…，y_n]^T尽可能互相独立.

图 2

Fig. 2

图 2 ICA的解混框图 Fig. 2 The unmixing diagram of ICA

(3)解混后得到源信号的估计Y，再将分离后的信号分量y_i(i=1，2，…，n)转化为图像，从中选择出变化图像，经阈值分割转化为二值图像后，得到最终的变化结果.

整个变化检测过程中最主要的是步骤(2)，而在步骤(2)中如何快速寻求正交变换矩阵U=[u₁，u₂，…，u_n]^T最为关键.现有的ICA算法大多采用基于负熵的固定点算法，这是由于负熵是度量随机变量非高斯性的最优工具(王翠茹等，2006).利用它来寻找正交矩阵U的每个分量u_i，使其投影y_i=u_i^Tz具有最大的非高斯性.其基本思想如下：对于有近似高斯分布的函数p(y)可表示为

上述过程每次提取一个独立源之后，需在下次提取前再次正交化，以保证每次提取的都是尚未提取的信源.由于计算Jacobi矩阵过程中涉及到矩阵求逆，计算量较大，所以一般将所有的Jacobi矩阵均取值J_ac(u₀)，以减少运算量，但收敛速度较小，可能导致算法不能最终收敛.

针对上述问题，现给出改进的基于牛顿迭代的固定点ICA算法(曾生根等，2003;李芳芳等，2010).假定已经求得u_i(k)，则u_i(k+1)可通过下面过程得到：

改进后的算法计算一次Jacobi矩阵可以实现n次迭代，相比通常的ICA算法具有更快的收敛速度，因此可以满足实时处理要求，且无需调整学习速率因子等其他动态参数，可靠性增强.

由于图像数据矩阵的庞大和ICA算法的复杂性，往往在独立分量提取之前对图像数据进行预处理，优化数据块的组合，使后续数据处理算法更高效.本文通过Contourlet变换将图像矩阵分块，然后对分块后图像数据进行独立分量分析，得到互相独立的分量，最后利用Contourlet反变换重构图像数据矩阵，以获得变化图像.具体实现框图如图 3所示.

图 3

Fig. 3

图 3 基于Contourlet变换和ICA的变化检测流程图 Fig. 3 Flowchart of change detection based on contourlet transform and ICA

具体实现步骤如下：

步骤1：利用Contourlet变换对多时相遥感图像X₁进行多尺度分解，再将分解的结果转化为向量形式，得到以低频分量C_1j和高频分量D₁₁¹，D₁₁²，…，D_1j^2j组成的分块向量，同理可得X₂的分块向量，即

步骤2：将Contourlet变换后的数据向量作为ICA算法的输入信号，然后通过独立分量分析实现变化区域向量的提取，即

步骤3：将变化区域的数据向量转化为矩阵，然后进行Contourlet逆变换得到变化图像，再经阈值分割后获得最终的变化检测结果.

Contourlet变换分解数据块，使得分解后的数据块保留了原始数据在各频率点及各方向上的特征，优化了ICA对源信号的估计，能够更真实地反映出图像块的变化.

为了验证本文算法的效果，分别利用本文算法及基于PCA、基于ICA、基于小波变换与ICA三种算法对大量多时相遥感图像进行了变化检测实验.实验环境为Intel(R)Core(TM)2，2.0 GHz，1GB内存，Matlab R2007b.在小波变换部分采用的是Daubechies正交规范小波基，分解到第三层.现选取其中两组实际的遥感图像数据集加以说明.第1组遥感图像真实数据集是分别于1996年8月和1998年8月拍摄的越南红河的图像，两幅图像为ERS-2遥感器获得的SAR图像，分别如图 4a和图 4b所示.两幅遥感图像的灰度级均为256，尺寸为512×512像素，变化区域主要为洪水淹没所致.

图 4

Fig. 4

图 4 第一组实际遥感图像 (a) 1996年8月越南红河洪水前图像； (b) 1996年8月越南红河洪水后图像. Fig. 4 The first group of real remote sensing images (a) Image before the flood in the red river, Vietnam, in August 1996; (b) Image after the flood in the red river, Vietnam, in August 1996.

第二组遥感图像数据集是Erlas8.5中自带的样品数据，由某地区的1987年和1992年的SPOT图像组成，分别如图 5a和图 5b所示.该地区的多个地方发生了变化，且变化类型复杂.两幅遥感图像的灰度级均为256，尺寸为512×512像素.

图 5

Fig. 5

图 5 第二组实际遥感图像 (a) 1987年SPOT影像数据； (b) 1992年SPOT影像数据. Fig. 5 The second group of real remote sensing images (a) SPOT imagery data in 1987; (b) SPOT imagery data in 1992.

针对上述两组实际遥感图像，分别利用基于PCA、基于ICA、基于小波变换与ICA以及本文算法分离出图像分量，如图 6和图 7所示.其中，a图是基于PCA算法分离出的图像分量，左图为图像的第一主分量，右图为图像的第二主分量；(b)是基于ICA算法分离出的图像分量；(c)是基于小波变换与ICA算法分离出的图像分量；(d)为本文算法分离出的图像分量.

图 6

Fig. 6

图 6 第一组遥感图像所分离出的分量 (a) 基于PCA分离出的第一主分量和第二主分量； (b) 基于ICA分离出的分量； (c) 基于小波变换与ICA分离出的分量； (d) 本文算法分离出的分量. Fig. 6 The separated components of the first group of remote sensing images (a) The first and second separated principal components based on PCA; (b) The separated components based on ICA; (c) The separated components based on wavelet transform and ICA; (d) The separated components based on the proposed algorithm.

图 7

Fig. 7

图 7 第二组遥感图像所分离出的分量 (a) 基于PCA分离出的第一主分量和第二主分量； (b) 基于ICA分离出的分量； (c) 基于小波变换与ICA分离出的分量； (d) 本文算法分离出的分量. Fig. 7 The separated components of the second group of remote sensing images (a) The first and second separated principal components based on PCA; (b) The separated components based on ICA; (c) The separated components based on wavelet transform and ICA; (d) The separated components based on the proposed algorithm.

由图 6a和图 7a可以看出，经过PCA分离出的第一主分量为背景图像，第二主分量为变化图像.从分离出的分量图像上可以看出图像中的地物信息部分混叠，背景图像中还包含部分未完全分离出的变化信息，说明PCA并不能完全消除遥感图像中各分量间的相关性.而基于ICA算法、基于小波变换与ICA算法和本文算法都能够较好地分离出变化图像，这说明ICA算法比PCA算法具有更好的去相关性和分离图像的能力.再比较图 6b、图 6c和图 6d，可以看出图 6b和图 6c的背景图像中还有部分变化区域的轮廓，而图 6d的背景图像中几乎没有混叠任何变化信息，变化量较少；图 7b和图 7c的背景图像中较暗区域为变化区域，而图 7d的背景图像几乎没有变化区域的信息，变化图像完全从背景图像中分离出来.因此，相比基于ICA算法、基于小波变换与ICA算法，本文算法能够更好地分离出变化信息.

为了验证本文算法的有效性，采用Otsu阈值法(陈琪等，2010)分别对图 6和图 7中的变化图像进行分割，所得到的变化检测结果如图 8和图 9所示.其中，(a)(b)(c)(d)分别是基于PCA、基于ICA、基于小波变换与ICA三种算法和本文算法的最终变化检测结果.可以看出，基于PCA算法的变化检测结果存在很多由不同时相中光照等因素导致的伪变化信息，而基于ICA算法、基于小波变换与ICA算法及本文算法都获得了较好的变化检测结果.此外，对于图 8中标出的区域，本文算法、基于小波变换与ICA算法两种算法都比基于传统ICA算法在凹面边缘轮廓处的检测更为细腻，边缘连续性较好，说明了通过Contourlet变换和小波变换分解能够优化ICA方法对边缘细节的处理，并且本文算法检测结果中的孤立点比基于小波变换与ICA算法少，变化检测效果更好.对于图 9中标出的区域，由上方的标出区域可以看出，基于传统ICA算法在此边缘处存在多个断点，而本文算法在边缘轮廓处比较完整，说明了Contourlet变换更好地保留了图像的边缘细节特征；此外，由中间的标出区域可以看出，本文算法检测出了一条完整的直线状的未变化区域，这与实际遥感图像的变化相吻合，而其他算法却未完整检测出来.这是由于本文算法利用Contourlet变换的多尺度、多方向和各向异性等特性，使得分解后的数据块保留了原始数据在各频率点及各方向上的特征，优化了ICA对源信号的估计，能更准确地反映变化信息.

图 8

Fig. 8

图 8 第一组遥感图像的变化检测结果 (a) 基于PCA的变化检测结果； (b) 基于ICA的变化检测结果； (c) 基于小波变换与ICA的变化检测结果； (d) 本文算法的变化检测结果. Fig. 8 Change detection results of the first group of remote sensing images (a) Change detection result based on PCA; (b) Change detection result based on ICA; (c) Change detection result based on wavelet transform and ICA; (d) Change detection result based on the proposed algorithm.

图 9

Fig. 9

图 9 第二组遥感图像的变化检测结果 (a) 基于PCA的变化检测结果； (b) 基于ICA的变化检测结果； (c) 基于小波变换与ICA的变化检测结果； (d) 本文算法的变化检测结果. Fig. 9 Change detection results of the second group of remote sensing images (a) Change detection result based on PCA; (b) Change detection result based on ICA; (c) Change detection result based on wavelet transform and ICA; (d) Change detection result based on the proposed algorithm.

依据这两组实际遥感图像的变化像元错判误差、非变化像元误判误差、总体精度和运行时间等客观定量指标对四种算法进行了比较，第2组遥感图像变化检测结果如表 1所示.本文算法与基于PCA、基于ICA、基于小波变换与ICA算法相比，精确度分别提高了1.8%、0.5%、0.3%，由于Contourlet变换能更好地描述图像中的轮廓和方向性纹理信息，很好地弥补了小波变换的不足，因此本文算法在边缘细节上处理得更好，错判误差较少，总体精度更高.在运算时间上，本文算法比其他三种算法分别减少了12 s、5 s、6 s.基于小波变换与ICA算法比基于ICA算法耗时，是因为小波变换只是对数据进行了优化排列，再加上小波变换自身所消耗的时间，整体时间上就没有优势.而本文算法通过Contourlet变换对数据进行多方向多尺度分解和优化组合，并采用改进的基于牛顿迭代的固定点ICA算法，大大缩短了整体的运算时间.综上所述，本文算法与基于PCA、基于ICA、基于小波变换与ICA三种算法相比，具有更高的精确度和更快的运算速度.

表 1(Table 1)

表 1 四种变化检测算法的精确度比较 Table 1 Comparison of four algorithms in accuracy 不同算法 变化像元错判误差 非变化像元误判误差 总体精度 运行时间/s 基于PCA算法 3.15% 0.43% 88.58 % 24 基于ICA算法 2.94% 0.36% 89.79% 17 基于小波变换与ICA算法 2.83% 0.31% 89.95% 18 本文算法 2.69% 0.26% 90.21% 12

表 1 四种变化检测算法的精确度比较 Table 1 Comparison of four algorithms in accuracy

本文提出了利用Contourlet变换和独立分量分析的多时相遥感图像的变化检测算法.通过Contourlet变换的多分辨率分解、较好的方向性和各向异性，提高了变化检测的精确度，降低了计算的复杂性，并通过改进的ICA算法提高了运算效率.实验结果表明：本文算法与现有的基于PCA、基于ICA、基于小波变换与ICA三种算法相比，在边缘细节方面处理得更好，获得的变化图像具有更高的精确度和较强的稳健性，且运算时间较少，可以扩展到更多时相的遥感图像序列，它是一种准确高效的多时相遥感图像变化检测算法.