CN102629378B - 基于多特征融合的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于多特征融合的遥感图像变化检测方法，主要解决现有变化检测方法漏检较多及变化检测整体精度不高的问题。其实现过程是：首先输入两时相遥感图像差值得到灰度差异图，计算灰度差异图的方差并设定阈值以决定是否对两时相图进行形态学预处理；接着计算邻域差分图和Gabor纹理特征差异图，使用Treelet算法融合三组差异图，对融合后的差异图像进行最大熵分割得到初始检测结果；最后通过面积阈值法去除部分孤立伪变化信息得到最终的检测结果图。实验表明，本发明能够较好的保持变化区域的边缘特征，在降低漏检率的同时有效抑制噪声的干扰，可用于环境保护、城市规划建设、自然灾害检测领域。

Description

基于多特征融合的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，具体来说是一种基于多特征融合的遥感图像变化检测方法，可用于环境保护、农业、水利及军事等诸多领域的遥感图像分析。

背景技术

遥感图像变化检测是遥感图像处理中的热点问题，它对同一地点在不同时期拍摄的两幅遥感图像进行比较和分析，从而检测出该地区的地物变化信息。该技术广泛应用国民经济和国防建设的诸多领域，如环境保护、城市规划建设、地理数据更新以及军事中道路、桥梁等战略目标的动态监视等。

地物变化在遥感图像上反映为图像灰度、形状和纹理等特征的变化，这些特征是进行遥感图像变化检测的主要依据。由于像素的灰度特征容易使用，与此相关的方法较多，其中最常见的方法就是对两时相图像的灰度特征进行比较，根据图像灰度信息的差别进行变化检测。Wei He于2010年提出一种基于欧式距离和最大熵分割的变化检测方法，首先采用欧式距离模型衡量两时相图像灰度值之间的相似性，构造差异影像，然后对该差异影像进行最大熵分割得到检测结果。图像的灰度信息提取较为方便，灰度差值法操作起来较为简单，然而不同时刻之间的遥感图像由于光照、辐射等因素造成图像灰度值存在偏差，仅仅利用像素灰度信息构造差异图像得到的检测结果往往精度不高。

相对于灰度特征，图像的边缘、纹理等特征更为稳定，不容易受到遥感图像时相变化的影响。随着对变化检测领域的深入研究，有关学者提出了许多利用其它特征或多特征融合的变化检测方法。S.Q.Huang等学者在文章“A New Change DetectionAlgorithm for SAR Images”中提出一种基于灰度共生矩阵和多数表决的变化检测技术，该方法首先构造五组基于灰度共生矩阵的纹理差异图像，分别对其进行分割，通过多数表决融合多组分割结果得到最终的检测结果。该方法容易累积分类误差，而且多数表决并不能保证得到最优的检测结果。Lingcao Huang等学者在文章“Anobject-based change detection approach by integrating intensity and texture differences”提出一种融合灰度特征和纹理特征的变化检测方法，其中纹理特征采用梯度矢量来描述，该方法得到的检测结果边缘保持较好，但是存在较多的虚警。

发明内容

本发明的目的在于针对上述已有遥感图像变化检测技术的不足，提出了一种基于多特征融合的遥感图像变化检测方法，以较好的保持变化区域的边缘信息，在抑制虚警的同时降低漏检率，提高变化检测的整体精度。

为实现上述目的，本发明的检测方法包括如下步骤：

(1)对输入的两幅已配准的大小均为m×n的两时相遥感图像I¹和I²进行差值运算，得到灰度差值图D1；

(2)计算灰度差值图D1的方差，设定方差阈值th1＝10，如果方差小于th1，则分别对两时相遥感图像I¹和I²进行形态学滤波，得到滤波后图像X¹和X²，，否则执行步骤(4)；

(3)对滤波后的两幅图像进行差分，用该差分结果替代步骤(1)中的灰度差值图D1；

(4)计算灰度差值图D1的邻域差分图D2；

(5)分别提取两时相遥感图像I¹和I²四个不同方向和四个不同频率下共16维的Gabor特征，再将两时相遥感图像I¹和I²相同方向、相同频率下的Gabor特征对应作差，得到16幅Gabor纹理特征差异图；

(6)使用Treelet算法对步骤(4)中得到的16幅Gabor纹理特征差异图进行降维，得到纹理差异图D3；

(7)使用Treelet算法融合三幅差异图D1、D2及D3，得到最终差异图D；

(8)对最终差异图D进行最大熵分割，得到初步变化检测结果图Z；再去除该检测结果Z中的孤立小噪声点，得到检测结果图ZF。

本发明与现有的技术相比具有以下优点：

(1)本发明采用基于Gabor变换的特征构造纹理差异图，Gabor特征对噪声有较强的免疫性，该纹理差异图能有效降低检测结果中的漏检率；

(2)本发明使用Treelet算法融合多组特征差异图，结合各个特征的不同优势，较好的保持了变化区域的边缘信息，并提高变化检测结果的整体精度；

(3)本发明通过设置面积阈值去除孤立的小噪声点，能抑制变化检测结果中存在的伪变化信息，降低虚警率。

附图说明

图1是本发明的实现流程图；

图2是本发明使用的两时相遥感图像及对应参考图；

图3是本发明通过Treelet融合得到的差异图；

图4是用本发明及对比方法对模拟遥感数据集实验得到的变化检测结果图；

图5是用本发明及对比方法对真实遥感数据集实验得到的变化检测结果图。

具体实施方式

参照图1，本发明的实施步骤如下：

步骤1，对输入的两个不同时相的已配准且大小均为m×n的遥感图像I¹和I²，如图2(a)和2(b)所示，进行差值得到灰度差值图D1；

D1＝|I¹-I²|。

步骤2，计算灰度差值图D1的方差，设定阈值th1＝10，如果方差小于th1，则分别对两时相图像I¹和I²进行形态学滤波，得到滤波后图像X¹和X²，否则执行步骤(4)。

对两时相图像I¹和I²分别进行形态学滤波的具体操作如下：

(2a)采用形态学开-闭滤波器对两时相遥感图像I¹进行滤波，即选取直线作为结构元素，定义长度分别为3和5两种尺寸、角度分别为0°、45°、90°、135°四种形状的八种结构元素，将相同角度不同尺寸的结构元素作为一类，共有四类结构元素，使用每一类结构元素对两时相遥感图像I¹进行串联滤波，即首先采用尺寸较小的结构元素进行滤波，接着采用尺寸较大的结构元素进行滤波，得到四幅滤波后图像，分别记为T¹¹，T¹²，T¹³，T¹⁴；

(2b)将T¹¹，T¹²，T¹³，T¹⁴均变换为m×n大小的列向量，组成初始样本X；

(2c)对初始样本X进行Treelet变换，具体步骤如下：

(2c1)定义Treelet变换的逐层聚类层数l＝0，1，...L-1，L为初始样本X中列向量的个数，L＝4，在第0层，每个变量采用初始样本X中的列向量表示，初始化基矩阵B⁰为L×L的单位矩阵及和变量下标集合Ω＝{1，2，...L}，计算初始样本X的初始协方差矩阵C⁰和初始相关系数矩阵M⁰：

C⁰＝[C_ij]

M⁰＝[M_ij]

其中C_ij＝E[(s_p-Es_p)(s_q-Es_q)]，表示初始协方差矩阵C⁰第i行第j列的计算值，表示初始相关系数矩阵M⁰第i行第j列的计算值，式中s_p，s_q表示初始样本X中两个不同的列向量，E表示求数学期望；

(2c2)当逐层聚类层数l≠0时，寻找相关系数矩阵M^l中最大的两个值，将最大值和次大值的对应位置序号分别记为α和β：

$(\mathrm{\α},\mathrm{\β})=\arg \underset{i,j\∈\mathrm{\Ω}}{\max }{M}_{\mathrm{ij}}^{l-1},$

这里i＜j，分别表示相关系数矩阵M^l中任意值的行和列，且只在和变量下标集合Ω内进行；

(2c3)计算Jacobi旋转矩阵J：

其中c_n＝cos(θ_l)，s_n＝sin(θ_l)；旋转角θ_l由C^l＝J^T C^l-1J及计算得到，表示协方差矩阵C^l中行、列位置序号分别为α、β的元素值，表示协方差矩阵C^l中行、列位置序号分别为β、α的元素值，J^T表示旋转矩阵J的转置；

${\mathrm{\θ}}_l={\tan }^{-1}\left[\frac{C_{\mathrm{\α\α}}^l-C_{\mathrm{\β\β}}^l\±\sqrt{{(C_{\mathrm{\α\β}}^l-C_{\mathrm{\β\β}}^l)}^2+4{\left(C_{\mathrm{\α\β}}^l\right)}^2}}{2C_{\mathrm{\α\β}}^l}\right]$ 且 $\left|{\mathrm{\θ}}_l\right|\≤\frac{\mathrm{\π}}{4};$

式中，表示协方差矩阵C^l中行、列位置序号均为α的元素值，表示协方差矩阵C^l中行、列位置序号均为β的元素值，表示协方差矩阵C^l中行位置序号为α、列位置序号均为β的元素值；

(2c4)由矩阵J更新去相关后的正交基B^l＝B^l-1J，协方差矩阵C^l＝J^T C^l-1J，其中B^l-1为更新前的正交基，B^l为更新后的正交基，C^l-1为更新前的协方差矩阵，C^l为更新后的协方差矩阵；

(2c5)利用步骤(2c4)中计算得到的更新后协方差矩阵C^l更新相关系数矩阵M^l-1为M^l，M^l-1为更新前的相关系数矩阵，M^l为更新后的相关系数矩阵：

$M^l=\left[M_{\mathrm{ij}}^l\right]$

$M_{\mathrm{ij}}^l=\frac{C_{\mathrm{ij}}^l}{\sqrt{C_{\mathrm{ii}}^l{C}_{\mathrm{jj}}^l}}$

其中它表示更新后相关系数矩阵M^l中第i行j列的更新值，为更新后协方差矩阵C^l中第i行j列的值，为更新后协方差矩阵C^l中第i行i列的值，为更新后协方差矩阵C^l中第j行j列的值；

(2c6)定义更新后正交基矩阵B^l中位置序号为α和位置序号为β的两个列向量分别为尺度函数和细节函数ψ_l，定义当前l层的尺度向量集合是尺度函数和上一层的尺度向量集合的合集，将差变量的位置序号β从和变量的下标集合Ω中去除，即Ω＝Ω\{β}；

(2c7)重复步骤(2c2)至(2c6)直至聚类的最高层l＝L-1，得到最终的正交基矩阵其中，ψ₁为第一层Treelet变换得到的细节函数，ψ_L-1为最高层Treelet变换得到的细节函数；

(2d)将初始样本X沿正交基矩阵B转置方向进行投影，即P＝X·B^T，并将投影向量P变回m×n大小的图像，得到两时相遥感图像中I¹的滤波后图像X¹，B^T表示正交基矩阵B的转置；

(2e)对两时相遥感图像I²重复步骤(2a)到(2d)，得到两时相遥感图像中I²的滤波后图像X²。

步骤3，对步骤(2)中两时相遥感图像I¹和I²滤波后的两幅图像X¹和X²进行差值，用该差值结果代替步骤(1)中的灰度差值图D1；

D1＝|X¹-X²|。

步骤4，按照下式计算灰度差异图D1的邻域差分图D2：

D2＝[D2(i，j)]

D2(i，j)＝[D1(i，j-1)+D1(i-1，j)+D1(i，j+1)+D1(i+1，j)]/4式中，i、j分别表示灰度差异图D1中像素点的行和列；

步骤5，分别计算两时相滤波后图像X¹和X²四个不同方向及四个不同尺度下的Gabor纹理特征图，将两时相图像相同角度、相同频率下的Gabor纹理特征图对应差值得到16幅Gabor纹理特征差异图，具体步骤如下：

(5a)分别对两时相遥感图像I¹和I²按照下式进行四个不同方向θ＝0，π/4，π/2，3π/4和四个不同频率f＝1/2，1/4，1/8，1/16的Gabor滤波，则两时相遥感图像I¹和I²分别得到16组Gabor滤波响应，依次表示为和h＝1，2，...16：

$G_h^1=\frac{1}{2{\mathrm{\π\σ}}^2}\exp (-\frac{{U1}^2+{V1}^2}{2{\mathrm{\σ}}^2})\exp \left(2\mathrm{\π}{j}^{\′}\mathrm{fU}1\right)$

式中U1＝icosθ+jsinθ，V1＝-isinθ+jcosθ，式中i、j分别表示两时相遥感图像中I¹当前像素点的行和列；

$G_h^2=\frac{1}{2{\mathrm{\π\σ}}^2}\exp (-\frac{{U2}^2+{V2}^2}{2{\mathrm{\σ}}^2})\exp \left(2\mathrm{\π}{j}^{\′}\mathrm{fU}2\right)$

式中U2＝icosθ+jsinθ，V2＝-isinθ+jcosθ，式中i、j分别表示两时相遥感图像中I²当前像素点的行和列，j′表示虚部，exp表示指数，θ为角度，σ为标准差，f为空间频率；

(5b)根据步骤(5a)中得到的Gabor滤波响应和利用下式提取Gabor纹理特征，得到两时相遥感图像I¹和I²的Gabor纹理特征图，依次表示为和

$F_h^1=\left[F_h^1(i.j)\right],h=\mathrm{1,2},...16$

其中 $F_h^1(i,j)=\frac{1}{N^2}\underset{{\mathrm{\Ω}}_{i,j}^{\′}}{\∫\∫}\cot \left(\mathrm{\ρ}\frac{G_h^1(i,j)}{\mathrm{\σ}\left(G_h^1\right)}\right)\mathrm{didj},$ 表示Gabor纹理特征图中第i行第j列的特征值，表示Gabor滤波响应中第i行第j列的值；

$F_h^2=\left[F_h^2(i,j)\right],$

其中 $F_h^2(i,j)=\frac{1}{N^2}\underset{{\mathrm{\Ω}}_{i,j}^{\′}}{\∫\∫}\cot \left(\mathrm{\ρ}\frac{G_h^2(i,j)}{\mathrm{\σ}\left(G_h^2\right)}\right)\mathrm{didj},$ 表示Gabor纹理特征图中第i行第j列的特征值，表示Gabor滤波响应中第i行第j列的值，N为窗口大小，Ω′_i，j表示像素点(i，j)周围N×N大小的邻域，σ表示标准差，ρ为常参数，取N＝5，ρ＝0.25；

(5c)对两时相遥感图像I¹和I²相同频率、相同方向下的Gabor纹理特征图和对应作差，得到16幅Gabor纹理特征差异图D_h：

$D_h=F_h^1-F_h^2,h=\mathrm{1,2},...16$

式中为两时相遥感图像中I¹的Gabor纹理特征图，为两时相遥感图像中I²的Gabor纹理特征图。

步骤6，使用Treelet算法对步骤(5)中得到的16维Gabor纹理特征差异图进行降维，取前三维能量较大的Gabor纹理特征差异图组成纹理差异图D3。

(6a)将步骤(5)中得到的16维Gabor纹理特征差异图D_h均变换为m×n大小的列向量d_h(h＝1，2，...16)，组成初始样本X＝[d₁，d₂，d₃，d₄，d₅，d₆，d₇，d₈，d₉，d₁₀，d₁₁，d₁₂，d₁₃，d₁₄，d₁₅，d₁₆]；

(6b)设最大聚类层数L＝16，对步骤(6a)中的初始样本X进行Treelet变换，如步骤(2c)中所述；

(6c)将初始样本X沿正交基矩阵B前三行转置方向进行投影，即P＝X·B(1，：)^T+X·B(2，：)^T+X·B(3，：)^T，并将投影向量P变回m×n大小的图像，得到纹理差异图D3。

步骤7，对灰度差异图D1、邻域差分图D2及纹理差异图D3使用Treelet变换进行融合，得到本发明的差异影像D。

(7a)将灰度差异图D1、邻域差分图D2及纹理差异图D3均变换为m×n大小的列向量d1′，d2′，d3′，组成初始样本X＝[d1′，d2′，d3′]；

(7b)设最大聚类层数L＝3，对步骤(7a)中的初始样本X进行Treelet变换，如步骤(2c)中所述；

(7c)将初始样本X沿正交基矩阵B转置的方向进行投影，即P＝X·B^T，并将投影向量P变回m×n大小的图像，得到本发明的差异影像D，如图3所示。

步骤8，对差异影像D进行最大熵分割，得到初始变化检测结果Z，再去除该检测结果Z中的孤立小噪声点，得到检测结果图ZF。

(8a)设分割阈值th可将差异影像D分成两个区域，目标区域和背景区域，分别按下式计算目标区域的熵H_O(th)和背景区域的熵H_B(th)：

$H_B\left(\mathrm{th}\right)=-\underset{s_B}{\mathrm{\Σ}}(p_{s_B}/p_k)\log (p_{s_B}/p_k)s_B=\mathrm{1,2},...\mathrm{th}$

$H_O\left(\mathrm{th}\right)=-\underset{s_O}{\mathrm{\Σ}}[p_{s_O}/(1-p_k)]\log [p_{s_O}/(1-p_k)]s_O=\mathrm{th}+1,\mathrm{th}+2,...Q-1$

$p_k=\underset{s_B=0}{\overset{\mathrm{th}}{\mathrm{\Σ}}}{p}_{s_B}$

$p_{s_B}=\frac{p_{s_B}}{m\×n}$

$p_{s_O}=\frac{n_{s_O}}{m\×n}$

其中Q为图像的最高灰度级，为图像中灰度级为s_B的像素点个数，为灰度级s_B出现的概率，为图像中灰度级为s_O的像素点个数，为灰度级s_O出现的概率，m、n均表示图像大小；

(8b)计算图像的总熵：H(th)＝H_O(th)+H_B(th)，将th在图像所有可能灰度级上进行迭代，使得总熵取最大值的th即为所求的最佳分割阈值；

(8c)将大于最佳分割阈值th的像素点标记为变化类，反之标记为非变化类，得到初步变化检测结果Z；

(8d)统计初步变化检测结果图Z中所有连续变化区域像素点的总个数，设置面积阈值th2＝20，将像素点总个数大于面积阈值的区域保留，反之则将其重新标记为非变化区域，生成变化检测结果图ZF，如图4和图5所示。

本发明的效果可通过以下实验结果与分析进一步说明：

1.实验数据

本发明的实验数据为一组模拟遥感数据集和一组真实遥感数据集，其中模拟遥感数据集的原始图像为470×335像素大小的位于英国Feltwell村庄农区的第3波段图像，配准误差为1.5个像素左右，模拟变化图像是通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工嵌入一些变化区域得到的，分别如图2(a)和图2(b)所示；真实遥感数据集是由2000年4月和2002年5月墨西哥郊外的两幅Landsat-7ETM+第4波段光谱图像组成，分别如图2(d)和图2(e)所示，两幅图像大小均为512×512，灰度级为256，它们之间发生的变化是由于大火破坏了大面积的当地植被引起的。

2.实验方法

方法1：Wei He在文章“Application of Euclidean norm in Multi-temporal RemoteSensing Images Change Detection”中提出的一种基于欧范数和最大熵分割的变化检测方法；

方法2：Lingcao Huang等学者在文章“An object-based change detection approachby integrating intensity and texture differences”提出的一种基于灰度特征和纹理特征融合的变化检测方法；

方法3：本发明方法。

3.实验内容与结果分析

实验1，用不同方法对如图2(a)和图2(b)所示的两幅不同时相的模拟数据集进行变化检测，变化检测参考图如图2(c)所示，对图2(a)和图2(b)的模拟遥感数据集进行变化检测得到的结果如图4所示，其中图4(a)为现有方法1得到的变化检测结果图，图4(b)为现有方法2得到的变化检测结果图，图4(c)为本发明方法得到的变化检测结果图，从图4可以看出，现有方法1和现有方法2结果图中存在较多的伪变化信息，本发明方法变化区域边缘保持较好，伪变化信息很少，比较接近参考图2(c)。

实验2，用不同方法对如图2(d)和图2(e)所示的两幅不同时相的真实遥感数据集进行变化检测，变化检测参考图如图2(f)所示，对图2(d)和图2(e)的真实遥感数据集进行变化检测得到的结果如图5所示，图5(a)为现有方法1得到的变化检测结果图，图5(b)为现有方法2得到的变化检测结果图，图5(c)为本发明方法得到的变化检测结果图，从图5可以看出，现有方法1结果图中都存在许多明显的噪声区域，现有方法2中同样存在大量的伪变化信息，本发明方法的结果图和参考图2(f)最为接近，变化检测结果的正确率最高。

对实验变化检测结果进行定量分析，即比较变化检测结果图与标准参考图，计算总错误数、虚景数、漏检数等性能指标，结果如表1所示。

表1.两组遥感图像数据变化检测结果评价指标

从表1上可以更为直观的看出，本发明方法能有效稳定地的获得遥感图像的变化区域，与其它两种现有方法相比较变化检测结果的精度最高。

Claims (4)

1.一种基于多特征差异图融合的遥感图像变化检测方法，包括如下步骤：

（1）对输入的两幅已配准的大小均为m×n的两时相遥感图像I¹和I²进行差值运算，得到灰度差值图D1；

（2）计算灰度差值图D1的方差，设定方差阈值th1=10，如果方差小于th1，则分别对两时相遥感图像I¹和I²进行形态学滤波，得到滤波后图像X¹和X²，否则执行步骤（4）；

（3）对滤波后的两幅图像进行差分，用该差分结果替代步骤（1）中的灰度差值图D1；

（4）计算灰度差值图D1的邻域差分图D2；

（5）分别提取两时相遥感图像I¹和I²四个不同方向和四个不同频率下共16维的Gabor特征，再将两时相遥感图像I¹和I²相同方向、相同频率下的Gabor特征对应作差，得到16幅Gabor纹理特征差异图；

（6）使用Treelet算法对步骤(4)中得到的16幅Gabor纹理特征差异图进行降维，得到纹理差异图D3；

（7）使用Treelet算法融合三幅差异图D1、D2及D3，得到最终差异图D；

（8）对最终差异图D进行最大熵分割，得到初步变化检测结果图Z；再去除该检测结果图Z中的孤立小噪声点，得到检测结果图ZF；

步骤（2）所述的分别对两时相遥感图像I¹和I²进行形态学滤波，得到滤波后图像X¹和X²，按如下步骤执行：

（2a）采用形态学开-闭滤波器对两时相遥感图像中I¹进行滤波，即选取直线作为结构元素，定义长度分别为3和5两种尺寸、角度分别为0°、45°、90°、135°四种形状的八种结构元素，将相同角度不同尺寸的结构元素作为一类，共有四类结构元素，使用每一类结构元素对两时相遥感图像中I¹进行串联滤波，即首先采用尺寸较小的结构元素进行滤波，接着采用尺寸较大的结构元素进行滤波，得到四幅滤波后图像，分别记为T¹¹，T¹²，T¹³，T¹⁴；

（2b）将T¹¹，T¹²，T¹³，T¹⁴均变换为m×n大小的列向量，组成初始样本X；

（2c）对初始样本X进行Treelet变换：

（2c1）定义Treelet变换的逐层聚类层数l=0,1,…L-1，L为初始样本X中列向量的个数，L=4，在第0层，每个变量采用初始样本X中的列向量表示，初始化基矩阵B⁰为L×L的单位矩阵及和变量下标集合Ω={1,2，…L}，计算初始样本X的初始协方差矩阵C⁰和初始相关系数矩阵M⁰：

C⁰＝[C_ij]

M⁰＝[M_ij]

其中C_ij＝E[(s_p-Es_p)(s_q-Es_q)]，表示初始协方差矩阵C⁰第i行第j列的计算值，表示初始相关系数矩阵M⁰第i行第j列的计算值，式中s_p，s_q表示初始样本X中两个不同的列向量，E表示求数学期望；

（2c2）当逐层聚类层数l≠0时，寻找相关系数矩阵M^l中最大的两个值，将最大值和次大值的对应位置序号分别记为α和β：

$(\mathrm{\α},\mathrm{\β})=\arg \underset{i,j\∈\mathrm{\Ω}}{\max }{M}_{\mathrm{ij}}^{l-1},$

这里i<j，分别表示相关系数矩阵M^l中任意值的行和列，且只在和变量下标集合Ω内进行；

（2c3）计算Jacobi旋转矩阵J：

其中c_n=cos(θl)，s_n=sin(θ_l)；旋转角θ_l由C^l＝J^TC^l-1J及计算得到，表示协方差矩阵C^l中行、列位置序号分别为α、β的元素值，表示协方差矩阵C^l中行、列位置序号分别为β、α的元素值，J^T表示旋转矩阵J的转置；

${\mathrm{\&theta;}}_l={\tan }^{-1}\left[\frac{C_{\mathrm{\&alpha;\&alpha;}}^l-C_{\mathrm{\&beta;\&beta;}}^l\&PlusMinus;\sqrt{{(C_{\mathrm{\&alpha;\&beta;}}^l-C_{\mathrm{\&beta;\&beta;}}^l)}^2+4{\left(C_{\mathrm{\&alpha;\&beta;}}^l\right)}^2}}{2C_{\mathrm{\&alpha;\&beta;}}^l}\right]$ 且 $\left|{\mathrm{\&theta;}}_l\right|\&le;\frac{\mathrm{\&pi;}}{4};$

式中，表示协方差矩阵C^l中行、列位置序号均为α的元素值，表示协方差矩阵C^l中行、列位置序号均为β的元素值，表示协方差矩阵C^l中行位置序号为α、列位置序号均为β的元素值；

（2c4）由矩阵J更新去相关后的正交基B^l＝B^l-1J，协方差矩阵C^l＝J^TC^l-1J，其中B^l-1为更新前的正交基，B^l为更新后的正交基，C^l-1为更新前的协方差矩阵，C^l为更新后的协方差矩阵；

（2c5）利用步骤（2c4）中计算得到的更新后协方差矩阵C^l更新相关系数矩阵M^l-1为M^l，M^l-1为更新前的相关系数矩阵，M^l为更新后的相关系数矩阵：

$M^l=\left[M_{\mathrm{ij}}^l\right]$

$M_{\mathrm{ij}}^l=\frac{C_{\mathrm{ij}}^l}{\sqrt{C_{\mathrm{ii}}^l{C}_{\mathrm{jj}}^l}}$

其中它表示更新后相关系数矩阵M^l中第i行j列的更新值，为更新后协方差矩阵C^l中第i行j列的值，为更新后协方差矩阵C^l中第i行i列的值，为更新后协方差矩阵C^l中第j行j列的值；

（2c6）定义更新后正交基矩阵B^l中位置序号为α和位置序号为β的两个列向量分别为尺度函数φ_l和细节函数ψ_l，定义当前l层的尺度向量集合{φ_l}是尺度函数φ_l和上一层的尺度向量集合{φ_l-1}的合集，将差变量的位置序号β从和变量的下标集合Ω中去除，即Ω=Ω\{β}；

（2c7）重复步骤（2c2）至（2c6）直至聚类的最高层l=L-1，得到最终的正交基矩阵B=[φ_L-1,ψ₁,...ψ_L-1]，其中φ_L-1∈{φ_L-1}，ψ₁为第一层Treelet变换得到的细节函数，ψ_L-1为最高层Treelet变换得到的细节函数；

（2d）将初始样本X沿正交基矩阵B转置方向进行投影，即P=X·B^T，并将投影向量P变回m×n大小的图像，得到两时相遥感图像中I¹的滤波后图像X¹，B^T表示正交基矩阵B的转置；

（2e）对两时相遥感图像中I²重复步骤（2a）到（2d），得到两时相遥感图像中I²的滤波后图像X²。

2.根据权利要求1所述的遥感图像变化检测方法，所述步骤（5），按照以下步骤执行：

（5a）分别对两时相遥感图像I¹和I²按照下式进行四个不同方向θ=0，π/4，π/2，3π/4和四个不同频率f=1/2，1/4，1/8，1/16的Gabor滤波，则两时相遥感图像I¹和I²分别得到16组Gabor滤波响应，依次表示为和h=1，2，…16：

$G_h^1=\frac{1}{2\mathrm{\&pi;}{\mathrm{\&sigma;}}^2}\exp (-\frac{U{1}^2+V{1}^2}{2{\mathrm{\&sigma;}}^2})\exp \left(2\mathrm{\&pi;}{j}^{\&prime;}\mathrm{fU}1\right)$

式中U1=icosθ+jsinθ，V1=-isinθ+jcosθ，式中i、j分别表示两时相遥感图像中I¹当前像素点的行和列；

$G_h^2=\frac{1}{2\mathrm{\&pi;}{\mathrm{\&sigma;}}^2}\exp (-\frac{U{2}^2+V{2}^2}{2{\mathrm{\&sigma;}}^2})\exp \left(2\mathrm{\&pi;}{j}^{\&prime;}\mathrm{fU}2\right)$

式中U2=icosθ+jsinθ，V2=-isinθ+jcosθ，式中i、j分别表示两时相遥感图像中I²当前像素点的行和列，j'表示虚部，exp表示指数，θ为角度，σ为标准差，f为空间频率；

（5b）根据步骤（5a）中得到的Gabor滤波响应和利用下式提取Gabor纹理特征，得到两时相遥感图像I¹和I²的Gabor纹理特征图，依次表示为F_h ¹和F_h ²，h=1,2,…16：

${F_h}^1=\left[{F_h}^1(i,j)\right]$

其中表示Gabor纹理特征图F_h ¹中第i行第j列的特征值，表示Gabor滤波响应中第i行第j列的值；

${F_h}^2=\left[{F_h}^2(i,j)\right]$

其中表示Gabor纹理特征图F_h ²中第i行第j列的特征值，表示Gabor滤波响应中第i行第j列的值，N为窗口大小，Ω′_i，j表示像素点（i，j）周围N×N大小的邻域，σ表示标准差，ρ为常参数，取N=5，ρ=0.25；

（5c）对两时相遥感图像I¹和I²相同频率、相同方向下的Gabor纹理特征图F_h ¹和F_h ²对应作差，得到16幅Gabor纹理特征差异图D_h，h=1,2,…16：

$D_h={F_h}^1-{F_h}^2,$

式中F_h ¹为两时相遥感图像中I¹的Gabor纹理特征图，F_h ¹为两时相遥感图像中I²的Gabor纹理特征图。

3.根据权利要求1所述的遥感图像变化检测方法，所述步骤（6），具体执行如下：

（6a）将步骤（5）中得到的16幅Gabor纹理特征差异图D_h均变换为m×n大小的列向量d_h(h=1,2,…16)组成初始样本X=[d₁,d₂,d₃,d₄,d₅,d₆,d₇,d₈,d₉,d₁₀,d₁₁,d₁₂,d₁₃,d₁₄,d₁₅,d₁₆]；

（6b）设最大聚类层数L=16，对步骤（6a）中的初始样本X进行Treelet变换；

（6c）将初始样本X沿正交基矩阵B前三行转置方向进行投影，即P=X·B(1,:)^T+X·B(2,:)^T+X·B(3,:)^T，并将投影向量P变回m×n大小的图像，得到纹理差异图D3。

4.根据权利要求1所述的遥感图像变化检测方法，其中步骤（7）所述的使用Treelet算法融合三幅差异图D1、D2及D3，得出最终差异图D，按如下步骤进行：

（7a）将灰度差异图D1、邻域差分图D2及纹理差异图D3均变换为m×n大小的列向量d1'，d2'，d3'，组成初始样本X=[d1′，d2′，d3′]；

（7b）设最大聚类层数L=3，对步骤（7a）中的初始样本X进行Treelet变换；

（7c）将初始样本X沿正交基矩阵B转置的方向进行投影，即P=X·B^T，并将投影向量P变回m×n大小的图像，得出最终差异影像D。