CN102169584B - 基于分水岭和treelet的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于分水岭和treelet的遥感图像变化检测方法，它属于图像处理技术领域。其实现过程是：输入两幅不同时相的遥感图像，并对每幅图像去噪，得到两时相的去噪后图像，并构造差值差异图像；应用分水岭算法得到差异图像的初始过分割标记图，并采用treelet算法对过分割标记区域的纹理特征进行合并，聚类得到变化类和非变化类的二值图像，利用变权马尔科夫随机场模型对二值图像进行空间约束得到后处理的边缘一致性图像，合并聚类二值图和边缘一致性图的连通区域，得到最终的变化检测结果。本发明能够有效地提高变化检测处理效率，保持图像的边缘信息，可用于灾情监测和土地利用对图像变化的检测。

Description

基于分水岭和treelet的遥感图像变化检测方法

技术领域

本发明属于数字图像处理技术领域，具体的说是一种基于分水岭和treelet的遥感图像变化检测方法，适用于在农业、环境、城市规划对图像的处理。

背景技术

遥感图像的变化检测是指对同一地理位置不同时期的遥感图像进行分析获得其中的变化信息，它是当前遥感数据处理技术的主要发展方向。

对配准后的两幅遥感图像的变化检测方法一般是先获取差异图，然后对差异图进行变化与非变化分类。这种方法由于对原始数据直接进行差异比较，不会改变数据本身，信息较为可靠。将差异图分为变化类与非变化类的一种常见方法是阈值法，人们尝试了各种精确估计阈值的方法，多阈值结果融合方法、统计分布建模阈值方法、基于优化算法等的阈值方法等，但变化阈值的自动、精确地获取一直是一个瓶颈问题。将差异图分为变化类与非变化类的另一种方法是分割方法。其中基于聚类的分割方法由于特征的相似度、类内和类间距离的度量等问题，效果不是十分理想。

为了提高变化检测的精确度，Bruzzone等学者在文章“Automatic Analysis of theDifference Image for Unsupervised Change Detection”中提出了自适应选取图像阈值的方法，但仍会出现阈值算法难以分辨的噪声区域，从而影响变化检测的结果。GuitingWang等学者在文章“A new change detection method based on non-parametric densityestimation and Markov random fields”中提出了基于马尔可夫场概率密度的变化检测分割方法，该方法是在随机统计样本的基础上获得总体数据的概率密度，需要大量的统计样本才能得到准确的结果，方法的正确性易受影响。

发明内容

本发明的目的在于克服上述已有遥感图像变化检测方法的不足，提出一种基于分水岭和treelet结合的遥感图像变化检测方法，以快速、准确的检测出连续的变化区域边缘，减少伪变化信息和运行时间，提高变化检测的准确度。

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

(1)对输入的两幅已配准的多时相遥感图像分别进行3×3像素的中值滤波，得到去噪后图像T1和T2；

(2)将去噪图像T1和T2空间位置对应的像素灰度值进行差值计算，得到一幅差值差异图X；

(3)对差值差异图X进行空间域到模糊域的变换，在模糊域修正隶属度后反变换回空间域，得到增强后的差异图像EX；

(4)采用标记分水岭算法将增强后的差异图像EX进行分割，得到L个区域的分水岭过分割图像Y，其中L为增强后的差异图像EX分水岭过分割后的区域数目，取值根据选取的图像确定，每一个标记区域记为lk(lk＝1，2，3，...，L)；

(5)在增强后的差异图像EX计算每一个区域lk内的均值、方差、平滑度、三阶矩、一致性和熵这六种纹理特征，组成特征列向量H_lk，并将所有L个区域的纹理特征向量H_lk合并，得到维数为6×L的特征矩阵H，再计算特征矩阵H的协方差矩阵C和相关系数矩阵A；

(6)采用treelet算法将特征矩阵H进行合并，得到纹理特征向量的合并标记S_lk；

(7)依据纹理特征向量合并标记S_lk，对步骤(4)所得的过分割图像Y中的L个区域中相同类别的区域进行合并，得到变化和非变化两分类图B_n；

(8)统计两分类图B_n的邻域像素信息，利用最大后验概率的变权高斯马尔科夫模型得到变化区域边缘细化的检测结果图像B_E；

(9)以变化检测结果图像B_E参考图，在图B_E中利用区域连通性保留B_E和B_n两幅图像中的相同区域，去除仅在B_E图中存在的噪声杂点，得到最终的变化检测结果图Z。

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

(1)本发明由于通过两个主要步骤来实现遥感图像变化检测的，即先对初始差异图像进行标记分水岭分割，再利用treelet算法获得标记区域特征有效聚类后对应的变化检测二值图，因而可在保证检测效果的同时有效减少运行时间，并通过聚类的方法降低分类错误数。

(2)本发明由于对获得的初始差异图像进行了模糊增强，实现变化信息的分散，且增强后的差异图能大幅减少分水岭变换的过分割图像标记数，继而减少后续treelet算法的运算量。

(3)本发明使用treelet算法对初始过分割标记图各区域的纹理特征进行聚类合并，该方法根据特征数据的相似性最大化实现，具有更强的抗噪性能。

(4)本方法对变化检测的结果使用马尔可夫模型进行后处理，得到边缘细化的后处理检测结果图，并通过初始检测结果图和后处理检测结果图的连通区域合并，删除多余噪声区域，得到更加有效的检测结果。

附图说明

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

图2是本发明使用的两时相遥感图像及其变化检测参考图；

图3是用本发明进行标记分水岭分割后的过分割图像；

图4是用本发明进行treelet合并区域后的二值图像；

图5是用本发明对模拟遥感图像实验的变化检测结果图；

图6是用本发明对真实遥感图像实验的变化检测结果图。

具体实施方式

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

步骤1，对输入的如图2(a1)和图2(a2)所示的两幅不同时相的已配准遥感图像，分别进行窗口大小为3×3像素的中值滤波，得到两幅滤波去噪后的图像T1和T2。步骤2，将去噪后的图像T1和T2在空间位置为(m，n)的像素点灰度值

和做差：

其中m＝1，2，…，M，n＝1，2，…，N，得到一幅差值差异图X。

步骤3，对差异图X进行模糊增强，得到增强后的差异图EX。

(3a)采用G变换将差异图X从空间域变换到模糊域：

${\mathrm{\υ}}_{\mathrm{mn}}=G\left(f_{\mathrm{mn}}\right)={[1+(f_{\max }-f_{\mathrm{mn}})/P_d]}^{-P_c}---\left(1\right)$

其中υ_mn表示差异图像X中在像素点(m，n)的像素f_mn的隶属度，f_max为图像像素的最大灰度值，P_c和P_d分别是指数型和倒数型模糊因子，在此取P_c＝1、P_d＝f_max-f_min，f_min为图像像素的最小灰度值；

(3b)采用下式将隶属度υ_mn修正为υ′_mn：

${\mathrm{\υ}}_{\mathrm{mn}}^{\′}=T\left({\mathrm{\υ}}_{\mathrm{mn}}\right)=\left\{\begin{array}{cc}2*{\left({\mathrm{\υ}}_{\mathrm{mn}}\right)}^2,& 0\≤{\mathrm{\υ}}_{\mathrm{mn}}\≤0.5\\ 1-2*{(1-{\mathrm{\υ}}_{\mathrm{mn}})}^2,& 0.5\≤{\mathrm{\υ}}_{\mathrm{mn}}\≤1\end{array}\right.---\left(2\right)$

其中T(υ_mn)是修正隶属度函数的表达式，隶属度υ_mn＝0.5所对应的像素点为渡越点，若隶属度大于渡越点则经过修正后会变得更大，若隶属度小于渡越点则经过修正后会变得更小，使得变化与非变化的两类更为分散；

(3c)将修正后隶属度υ′_mn进行G^-1反变换回到图像空间域，生成新的像素灰度值f′_mn：

$f_{\mathrm{mn}}^{\′}=G^{-1}\left({\mathrm{\υ}}_{\mathrm{mn}}^{\′}\right)=f_{\max }-P_d({\left({\mathrm{\υ}}_{\mathrm{mn}}^{\′}\right)}^{-1/P_c}-1)---\left(3\right)$

得到增强后的差异图像EX。

步骤4，采用标记符分水岭算法对增强后的差异图EX进行分割，得到过分割图像Y。

(4a)寻找增强后的差异图像EX中所有相同的局部极小值像素，将该局部极小值像素组成的连通区域作为内部标记；以内部标记的灰度值极小值点为中心进行分水岭变换，将得到的分水岭标记线作为外部标记；

(4b)计算增强后的差异图像EX的形态学梯度图像，删除梯度图像中内部和外部标记以外的局部最小区域，对删除后的梯度图像进行分水岭变换得到分水岭边界线，对边界线上的像素按梯度最大值进行八邻域连接得到轮廓线，即分水岭分割的图像，如图3所示；

(4c)以分水岭分割图像中的轮廓线上每个像素点为中心，统计窗口大小为3×3像素的邻域内其他各像素的类别标记，将该点像素邻域窗口内数目最多的标记赋给轮廓线上的当前像素点，当轮廓线上的所有像素点均赋予类别标记时，就得到了去除轮廓线的分水岭过分割标记图像Y，图像Y中过分割区域的总数为L，该L取值根据选取的图像确定，每一个标记区域记为lk(lk＝1，2，3，...，L)。

步骤5，计算增强后的差异图像EX中每一个标记区域lk内的纹理特征向量H_lk，将所有L个区域的纹理特征向量H_lk合并，得到特征矩阵H，再计算特征矩阵H的协方差矩阵C和相关系数矩阵A。

(5a)在增强后的差异图像EX中分别计算每一个标记区域lk内所有Q个像素的均值μ、方差σ、平滑度R、三阶矩w、一致性I和熵E这六种纹理特征，具体计算公式如下：

$\mathrm{\μ}=\underset{q=0}{\overset{Q-1}{\mathrm{\Σ}}}{z}_{q}P\left(z_q\right)$

$\mathrm{\σ}=\sqrt{\underset{q=0}{\overset{Q-1}{\mathrm{\Σ}}}{(z_q-\mathrm{\μ})}^{2}P\left(z_q\right)}$

R＝1-1/(1+σ²)

$w=\underset{q=0}{\overset{Q-1}{\mathrm{\Σ}}}{(z_q-\mathrm{\μ})}^{3}P\left(z_q\right)$

$I=\underset{q=0}{\overset{Q-1}{\mathrm{\Σ}}}{P}^2\left(z_q\right)$

$E=-\underset{q=0}{\overset{Q-1}{\mathrm{\Σ}}}P\left(z_q\right){\log }_{2}P\left(z_q\right)$

其中，z_q是表示像素第q个灰度级的灰度值，P(z_q)是该标记区域lk内所以像素的灰度级直方图；

(5b)将每一个标记区域lk的六个纹理特征用一个6×1维的列向量H_lk＝[μ_lk，σ_lk，R_lk，w_lk，U_lk，E_lk]^T表示，则整个增强后差异图像EX的L个区域的纹理特征特征向量H_lk构成一个6×L维的特征矩阵H，表示为：

$H=\left[\begin{array}{ccccc}{x}_{\mathrm{1,1}}& L& x_{1,h}& L& x_{1,L}\\ M& M& M& M& M\\ x_{k,1}& L& x_{k,h}& L& x_{k,L}\\ M& M& M& M& M\\ x_{6,1}& L& x_{6,h}& L& x_{6,L}\end{array}\right]---\left(4\right)$

其中x_k，h表示第h个纹理特征向量H_h的第k个特征，这里1≤h≤L，1≤k≤6；(5c)计算特征矩阵H的协方差系数矩阵C：

$C_{i,j}=\frac{\underset{k=1}{\overset{6}{\mathrm{\Σ}}}(x_{k,i}-{\stackrel{\‾}{x}}_i)(x_{k,j}-{\stackrel{\‾}{x}}_j)}{\sqrt{\underset{k=1}{\overset{6}{\mathrm{\Σ}}}{(x_{k,i}-{\stackrel{\‾}{x}}_i)}^2\underset{k=1}{\overset{6}{\mathrm{\Σ}}}{(x_{k,j}-{\stackrel{\‾}{x}}_j)}^2}}---\left(5\right)$

其中C_i，j表示方差矩阵C第i行第j列的计算值，

和

分别表示第i个和第j个纹理特征向量的均值，这里1≤i，j≤L；

(5d)计算特征矩阵H的相关系数矩阵A：

$A=\left[\begin{array}{cccc}{A}_{\mathrm{1,1}}& A_{\mathrm{1,2}}& L& A_{1,L}\\ A_{\mathrm{2,1}}& A_{\mathrm{2,2}}& L& A_{2,L}\\ M& M& M& M\\ A_{L,1}& A_{L,2}& L& A_{L,L}\end{array}\right]---\left(6\right)$

其中，相关系数矩阵A中的元素A_i，j的计算公式如下：

$A_{i,j}=\frac{C_{i,j}}{\sqrt{C_{i,i}{C}_{j,j}}}---\left(7\right)$

其中A_i，j为相关系数矩阵A中第i行第j列的计算值，C_i，j为协方差矩阵C中第i行第j列的值，C_i，i为协方差矩阵C中第i行第i列的值，C_j，j为协方差矩阵C中第j行第j列的值。

步骤6，使用treelet算法将区域特征矩阵H由L类逐层聚类合并至两类。

(6a)定义聚类的层数l＝0，1，…，L-2，当l＝0时，初始化和变量，即初始特征矩阵为S⁽⁰⁾＝H，差变量为D⁽⁰⁾为空集，并将和变量的集合Ω可用下标集表示为Ω＝{1，2，...，L}，将差变量集合Φ作为空集，将正交Dirac基B⁽⁰⁾＝[φ_0，1，φ_0，2，φ_0，3，…φ_0，L]作为L×L维单位矩阵；

(6b)当l≠0时，寻找相关系数矩阵A中最大的两个值，将最大值和次大值的对应位置序号分别记为α和β：

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

(6c)计算Jacobi旋转矩阵J：

$J=\left[\begin{array}{ccccccc}1& L& 0& L& 0& L& 0\\ M& O& M& & M& & M\\ 0& L& c_n& L& s_n& L& 0\\ M& & M& O& M& & M\\ 0& L& {-s}_n& L& c_n& L& 0\\ M& & M& & M& O& M\\ 0& L& 0& L& 0& L& 1\end{array}\right]---\left(9\right)$

其中c_n＝cos(θ_l)，s_n＝sin(θ_l)；旋转角θ_l由C^(l)＝J^TC^(l-1)J及

计算得：

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

由矩阵J更新去相关后的正交基B^(l)＝B^(l-1)J，协方差矩阵C^(l)＝J^TC^(l-1)J，其中B^(l)为更新后的正交基，C^(l)为更新后的协方差矩阵，并利用协方差矩阵C^(l)根据(7)式更新相关系数矩阵A^(l)；

(6d)定义正交基矩阵B^(l)中第α列和第β列的列向量分别为尺度函数φ_l和细节函数ψ_l，定义当前l层的尺度向量集合{φ_l}是尺度函数φ_l和上一层的尺度向量集合{φ_l-1}的合集，将β差变量从和变量集合Ω中去除，即Ω＝Ω\{β}，加入到差变量集合Φ中，即Φ＝{β}；

(6e)重复步骤(6b)至(6d)直至聚类的最高层l＝L-2，可得最终的正交基矩阵B＝[φ_L-2，1，φ_L-2，2，ψ₁，...ψ_L-2]，其中φ_L-2，1，φ_L-2，2∈{φ_L-2}；

(6f)将特征矩阵H沿正交基矩阵B的方向进行投影，计算投影后矩阵的能量矩阵En：

En＝(H×B′)·(H×B′)^T    (11)

提取能量阵En中的前两行得到尺度向量{φ_L-2}对应的能量矩阵Em(2×L)即为聚类得到的结果，按照各自对应尺度的能量大小进行排序，得到特征的合并标记S_lk。

步骤7，根据过分割图像Y中的任意标记区域lk对应的类别标记S_lk，将Y中具有相同标记值的区域进行合并，得到一幅两分类图像；再分别统计该两分类图像中各自类别的像素数目，将数目较多的一类标记为0，即确定为非变化类，将数目较少的另一类标记为1，即确定为变化类，最终得到变化和非变化两分类图B_n，如图4所示。

步骤8，对增强后的差异图像EX利用变权马尔科夫随机场模型进行后处理，得到边缘细化的变换检测结果图像B_E。

(8a)定义增强后差异图像EX的像素集合ρ＝{c＝(m，n)|1≤m≤M，1≤n≤N}中的像素点c的灰度值为f_c，用随机场F作为增强后差异图像EX，将检测结果图像记为类标随机场ω＝{ω_c|1≤c≤M×N}，并由马尔科夫随机场模型求解满足最小化能量函数U(ω|F)的类标随机场ω，其中能量函数U(ω|F)由似然能量U_data和先验能量U_context两部分构成；

(8b)初始化ω为检测结果图像B_n，似然能量U_data由二维高斯概率密度计算：

$U_{\mathrm{data}}=\ln p\left(F\right|\mathrm{\ω})=\underset{c\∈\mathrm{\ρ}}{\mathrm{\Σ}}(\frac{1}{2}\ln \left(2\mathrm{\π}{\mathrm{\σ}}_{{\mathrm{\ω}}_c}^2\right)+\frac{{(f_c-{\mathrm{\μ}}_{{\mathrm{\ω}}_c})}^2}{{2\mathrm{\σ}}_{{\mathrm{\ω}}_c}^2})---\left(\text{12}\right)$

其中ω_c表示增强后差异图像EX的像素点c处的类别标记，

和

分别表示集合ρ内所有具有相同标记的像素灰度均值和方差；

(8c)计算先验能量U_context：

$U_{\mathrm{context}}=\ln p\left(\mathrm{\ω}\right)=\underset{d\∈O}{\mathrm{\Σ}}-V\left({\mathrm{\ω}}_d\right)---\left(13\right)$

其中O表示ρ中所有为3×3像素大小的邻域窗口的集合，d为O中任意邻域窗口的中心像素点，V(ω_d)是与邻域信息相关的势函数，计算公式为：

$V\left({\mathrm{\ω}}_d\right)=\left\{\begin{array}{cc}-W,& {\mathrm{\ω}}_l={\mathrm{\ω}}_d\\ 0,& {\mathrm{\ω}}_l\≠{\mathrm{\ω}}_d\end{array}\right.---\left(14\right)$

其中ω_l为邻域窗口内中心像素d点以外像素的类别标记，W为惩罚因子，当W取1时即为原始马尔科夫随机场模型，先验能量U_context的取值范围为-8≤U_context≤0。本发明采用的是变权马尔科夫随机场模型，先验能量根据邻域变化信息获得，该邻域区域的变化信息计算公式为：

其中t_d表示当前像素点d的O邻域内各像素的变化信息统计值，f_d为邻域O内中心像素点d的灰度值，μ_O为邻域窗口内像素点的均值。统计t_d的最大值最小值，将邻域变化信息统计值t_d映射至[-8，0]上即得先验能量U_context；

(8d)利用最大后验概率算法求解满足最小化能量函数U(ω|F)的类别标记ω：

U(ω|F)＝U_data+U_context              (15)

当求解达到迭代的指定次数K，或迭代前后类别标记ω的变化满足收敛阈值T时，将计算得到的标记随机场ω输出作为检测结果图像B_E。

步骤9，以步骤(8)的检测结果图像B_E作为参考图，在检测结果图像B_E中利用区域连通性保留检测结果图像B_E和检测结果图像B_n两幅图像中的相同区域，去除仅在检测结果图像B_E中存在的噪声杂点即得最终的检测结果图Z。

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

1.实验数据

本发明的实验数据为一组模拟变化遥感图像和一组真实多时相遥感图像共两组图像进行实验，每一组遥感图像都有变化检测结果参考图。模拟原始遥感图像为470×335像素大小、位于英国Feltwell村庄农区的ATM(Airborne Thematic Mapper)第3波段的图像，模拟变化遥感图像是通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工的嵌入一些变化区域得到的。真实遥感图像为243×306像素大小的Landsat5TM第4波段光谱图像，由1992年8月和1994年8月意大利Elba岛西部地区的两幅图像组成，变化区域是由于大火破坏了当地植被所致。

2.对比实验及实验评价指标

本发明使用的对比实验方法如下所述：

对比方法1，是Lorenzo Bruzzone等学者在文章“Automatic Analysis of theDifference Image for Unsupervised Change Detection”中提出了的自适应选取图像阈值的方法，对差异图像考虑像素空间邻域信息的马尔科夫随机场模型的差异图像阈值分类。

对比方法2，是Guiting Wang等学者在文章“Anew change detection method basedon non-parametric density estimation and Markov random fields”中提出的基于马尔可夫场的非参数概率密度的变化检测方法。首先将两时相遥感图像进行差值运算得到差异图像，对差异图像进行模糊C均值分类得到变化检测结果，在此基础上采用变权马尔可夫模型分类得到最终的检测结果。

实验的最后是对变化检测结果进行评价和分析。将变化检测结果图与参考图进行主观视觉比较及客观比较，评价指标包括虚警像素数、漏检像素数和总错误像素数，及各方法的实验运行时间。

3.实验的内容与分析

实验1，用不同方法对两幅不同时相的模拟变化遥感图像进行变化检测，其中原始的两幅遥感图像分别如图2(a1)和图2(a2)所示，变化检测参考图如图2(a3)所示，对图2(a1)和图2(a2)的模拟变化遥感图像进行变化检测得到的结果如图5所示，其中图5(a)为用本发明方法得到的变化检测结果图，图5(b)为对比方法1得到的变化检测结果图，图5(c)为对比方法2得到的变化检测结果图。从图5可以看出三种方法在细节上的保持都比较好，接近于变化检测参考图的图2(a3)。

实验2，用不同方法对两幅不同时相的真实遥感图像进行变化检测，其中原始的两幅遥感图像分别如图2(b1)和图2(b2)所示，其变化检测参考图如图2(b3)所示，对图2(b1)和图2(b2)的真实遥感图像进行变化检测得到的结果如图6所示，其中图6(a)为用本发明方法得到的变化检测结果图，图6(b)为对比方法1得到的变化检测结果图，图6(c)为对比方法2得到的变化检测结果图。从图6中可看出，两种对比方法的检测结果都存在较多的噪声区域，而本发明的检测结果图与变化检测参考图的图2(b3)最为接近。

表1列出了两组遥感图像变化检测结果的评价指标，表2为两组遥感图像变化检测实验中本发明和对比方法的运行时间比较。

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

从表1的评价指标中可以更为直观的看出，本发明对遥感图像变化检测的有效性。

表2.遥感图像变化检测运行时间对比

(单位：s)	模拟图像	真实图像
			对比方法1	31.12	15.51
对比方法2	90.36	56.44
			本发明方法	30.21	12.56

从表2中可以看出，本发明对遥感图像的变化检测能有效地减少运行时间。

Claims (1)

1.一种基于分水岭和treelet的遥感图像变化检测方法，包括如下步骤：

（1）对输入的两幅已配准的多时相遥感图像分别进行3×3像素的中值滤波，得到去噪后图像T1和T2；

（2）将去噪图像T1和T2空间位置对应的像素灰度值进行差值计算，得到一幅差值差异图X；

（3）对差值差异图X进行空间域到模糊域的变换，在模糊域修正隶属度后反变换回空间域，得到增强后的差异图像EX；

（4）采用标记分水岭算法将增强后的差异图像EX进行分割，得到L个区域的分水岭过分割图像Y，其中L为增强后的差异图像EX分水岭过分割后的区域数目，取值根据选取的图像确定，每一个标记区域记为lk，lk＝1,2,3，…,L；

（5）在增强后的差异图像EX计算每一个区域lk内的均值、方差、平滑度、三阶矩、一致性和熵这六种纹理特征，组成特征列向量H_lk，并将所有L个区域的特征列向量H_lk合并，得到维数为6×L的特征矩阵H，再计算特征矩阵H的协方差矩阵C和相关系数矩阵A；

（6）采用treelet算法将特征矩阵H进行合并，得到纹理特征向量的合并标记S_lk：

（6a）定义聚类的层数l＝0,1,…,L-2，当l＝0时，初始化和变量为S⁽⁰⁾=H，差变量D⁽⁰⁾为空集，并将和变量的集合Ω用下标集表示为Ω={1,2,…,L}，将差变量集合Φ作为空集，将正交基B⁽⁰⁾取为L×L维单位矩阵；

（6b）当l≠0时，寻找相关系数矩阵A中最大的两个值，将最大值和次大值的对应位置序号分别记为α和β：

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

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

（6c）计算Jacobi旋转矩阵J：

其中c_n=cos(θ_l)，s_n=sin(θ_l)；旋转角θ_l由C^(l)＝J^TC^(l-1)J及

计算得到：

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

由矩阵J更新去相关后的正交基B^(l)＝B^(l-1)J，协方差矩阵C^(l)＝J^TC^(l-1)J，其中B^(l)为更新后的正交基，C^(l)为更新后的协方差矩阵，并利用C^(l)更新相关系数矩阵A^(l-1)为A^(l)：

$A_{i,j}^{\left(l\right)}=\frac{C_{i,j}^{\left(l\right)}}{\sqrt{C_{i,i}^{\left(l\right)}{C}_{j,j}^{\left(l\right)}}},$

其中为相关系数矩阵A^(l)中第i行j列的更新值，

为协方差矩阵C^(l)中第i行j列的值，

为协方差矩阵C^(l)中第i行i列的值，

为协方差矩阵C^(l)第中j行j列的值；

（6d)定义正交基矩阵B^(l)中第α列和第β列的列向量分别为尺度函数φ_l和细节函数ψ_l，定义当前l层的尺度向量集合{φ_l}是尺度函数φ_l和上一层的尺度向量集合{φ_l-1}的合集，将β差变量从和变量集合Ω中去除，即Ω=Ω\{β}，加入到差变量集合Φ中，即Φ={β}；

（6e）重复步骤（6b）至（6d）直至聚类的最高层l＝L-2，可得最终的正交基矩阵B=[φ_L-2，1,φ_L-2，2,ψ₁，…ψ_L-2]，其中φ_L-2，1,φ_L-2，2∈{φ_L-2}；

（6f）将特征矩阵H沿正交基矩阵B的方向进行投影，计算投影后矩阵的能量矩阵En＝(H×B′)·(H×B′)^T，提取能量阵En中的前两行得到尺度向量{φ_L-2}对应的能量矩阵Em(2×L)即为聚类得到的结果，按照各自对应尺度的能量大小进行排序，得到特征的合并标记S_lk；

（7）依据纹理特征向量合并标记S_lk，对步骤（4）所得的过分割图像Y中的L个区域中相同类别的区域进行合并，即先根据过分割图像Y中的任意区域lk对应的类别标记S_lk，将Y中具有相同标记值的区域进行合并，得到一幅两分类图像；再分别统计该两分类图像中各自类别的像素数目，将数目较多的一类标记为0，确定为非变化类，将数目较少的另一类标记为1，确定为变化类，最终得到变化和非变化两分类图B_n；

（8）统计两分类图B_n的邻域像素信息，利用最大后验概率的变权高斯马尔科夫模型得到变化区域边缘细化的检测结果图像B_E；

（9）以变化检测结果图像B_E参考图，在图B_E中利用区域连通性保留B_E和B_n两幅图像中的相同区域，去除仅在B_E图中存在的噪声杂点，得到最终的变化检测结果图Z。

Images