CN102663724A - 基于自适应差异图的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于自适应Treelet构造差异图的遥感图像变化检测方法，它属于图像处理技术领域，主要解决现有技术中变化检测精度不足的问题。其实现过程是：对输入的两幅不同时相遥感图像进行块处理，即计算第二幅图像搜索窗内的图像块与第一幅图像的中心图像块的差值，得到样本矩阵；利用Treelet算法对样本矩阵进行聚类，得到自适应差异图；计算差值差异图和自适应差异图的Otsu阈值，利用该阈值融合差值差异图和自适应差异图得到最终差异图，对最终差异图进行Otsu阈值分割，得到变化检测结果。本发明能够有效地提高变化检测精度，保持图像的边缘信息，可用于灾情监测和土地利用。

Description

基于自适应差异图的遥感图像变化检测方法

技术领域

本发明属于数字图像处理技术领域，具体的说是一种基于自适应Treelet构造差异图的遥感图像变化检测方法，适用于环境、城市规划中的数字图像处理。

背景技术

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

对配准后的两幅遥感图像的变化检测方法通常是通过直接运算获取差异图像，然后对差异图像进行分析处理，得到变化与非变化分类的检测结果图。这种方法具有相对简单、直接、易于实现的特点，缺点是获取的差异图像对遥感图像数据的预处理要求较高，差异图中的噪声对检测结果的影响很大。

为了改善直接差异图像的检测效果，提高变化检测的精度，Lingcao Huang等学者在文章“An Object-based Change Detection Approach by Integrating Intensity andTexture Differences”中提出了加权融合Sobel边缘纹理图和差值差异图得到概率差异图，利用阈值对融合后的概率差异图进行分割得到变化检测图像，在保持变化区域的同时引入了大量的噪声点，从而影响了变化检测的结果。Wei He在文章“Applicationof Euclidean norm in Multi-temporal Remote Sensing Image Change Detection”中提出的利用两幅遥感图像之间差值的欧式距离范数表示差异图，再利用最大熵阈值法提取变化部分，该方法抑制了噪声信息，但对边缘等细节保持不好，降低了变化区域的视觉效果。

发明内容

本发明的目的在于克服上述已有构造差异图像变化检测方法的不足，提出一种基于自适应差异图的遥感图像变化检测方法，以有效、准确的检测出变化区域边缘，减少伪变化信息，提高变化检测的准确度。

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

(1)输入两幅已配准的多时相遥感图像T1和T2，并将图像T1和T2空间位置对应像素的灰度值进行差值计算，得到差值差异图X；

(2)对第一幅图像T1以像素点i为中心取3×3大小的滑窗E_i，对第二幅图像T2同样以像素点i为中心取7×7像素大小搜索窗G_i，i为当前像素点，在搜索窗G_i内逐像素取3×3大小的滑窗，得到滑窗集合{F_j}，F_j为滑窗集合内以j为中心点的一个滑窗，定义滑窗集合{F_j}内的滑窗数目为Q个，此处Q＝25；

(3)分别将第二幅图像T2得到的Q个滑窗与第一幅图像T1的滑窗Ei进行差值计算，并将差值的结果用列向量进行表示，得到样本矩阵H，该样本矩阵H中包含Q个列向量且维数为9×Q，将第二幅图像T2的滑窗F_i与第一幅图像T1的滑窗E_i的差值列向量表示为h_ii；

(4)利用Treelet算法对样本矩阵H进行自适应聚类，得到与差值列向量h_ii相同类别标记的q个列向量，其中1≤q≤Q，将这q个列向量组成矩阵R，并计算该矩阵R的均值m_i，将均值m_i作为自适应差异图D中像素点i的灰度值；

(5)对两幅遥感图像T1和T2中的每个像素点，重复进行步骤2至步骤4，得到自适应差异图D；

(6)用Otsu阈值法分别计算差值差异图X和自适应差异图D的阈值K₁和K₂；

(7)利用Otsu阈值K₁和K₂融合差值差异图X和自适应差异图D，得到最终差异图I；

(8)对最终差异图I利用Otsu阈值法进行分割，得到变化检测的结果图像Z。

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

(1)本发明通过对相似窗口进行Treelet聚类构造自适应差异图，能有效的实现变化信息的分散，突出对比度，提高了检测准确性。

(2)本发明通过对差值差异图和自适应差异图进行图像融合，保留了差值差异图的细节信息和自适应差异图中的变化区域，在保证准确性的同时降低了漏检影响。

附图说明

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

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

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

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

具体实施方式

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

步骤1，输入两幅已配准的多时相遥感图像T1和T2，将图像T1和图像T2对应位置(x，y)处的像素点灰度值

和

进行差值计算，得到差值差异图1≤x≤M，1≤y≤N，其中M为遥感图像T1的长，N为遥感图像T1的宽。

步骤2，对第一幅图像T1以像素点i为中心取3×3大小的滑窗E_i，对第二幅图像T2同样以像素点i为中心取7×7像素大小搜索窗G_i，i为当前像素点，在搜索窗G_i内逐像素取3×3大小的滑窗，得到滑窗集合{F_j}，F_j为滑窗集合内以j为中心点的一个滑窗，j为搜索窗G_i内任意像素点，设滑窗集合{F_j}内的滑窗数目为Q个，此处Q＝25。

步骤3，分别将第二幅图像T2得到的Q个滑窗与第一幅图像T1的滑窗E_i进行差值计算，并将差值的结果用列向量进行表示，得到列向量组成的样本矩阵H，该样本矩阵H中包含Q个列向量且维数为9×Q：

<1>

其中x_p，t为滑窗集合{F_j}内第p个滑窗与滑窗Ei进行差值计算后的列向量表示，

1≤p≤Q，t为下标，1≤t≤9，h_ii表示将第二幅图像T2的滑窗F_i与第一幅图像T1的滑窗E_i的差值列向量，F_i为第二幅图像T2内以i为中心点的一个滑窗。

步骤4，利用Treelet算法对样本矩阵H进行自适应聚类，得到自适应差异图D：

(4a)计算样本矩阵H的协方差系数矩阵C：

$C=\left[\begin{array}{cccc}{C}_{\mathrm{1,1}}& C_{\mathrm{1,2}}& \&CenterDot;\&CenterDot;\&CenterDot;& C_{1,Q}\\ C_{\mathrm{2,1}}& C_{\mathrm{2,2}}& \&CenterDot;\&CenterDot;\&CenterDot;& C_{2,Q}\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& C_{u,v}& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ C_{Q,1}& C_{Q,2}& \&CenterDot;\&CenterDot;\&CenterDot;& C_{Q,Q}\end{array}\right]$ <2>

其中，协方差系数矩阵C中的任意一个元素C_uv的计算公式如下

$C_{u,v}=\frac{\underset{k=1}{\overset{9}{\mathrm{\&Sigma;}}}(x_{k,u}-{\stackrel{\&OverBar;}{x}}_u)(x_{k,v}-{\stackrel{\&OverBar;}{x}}_v)}{\sqrt{\underset{k=1}{\overset{9}{\mathrm{\&Sigma;}}}{(x_{k,u}-{\stackrel{\&OverBar;}{x}}_u)}^2\underset{k=1}{\overset{9}{\mathrm{\&Sigma;}}}{(x_{k,v}-{\stackrel{\&OverBar;}{x}}_v)}^2}}$ <3>

其中C_u，v表示协方差矩阵C中第u行第v列的计算值，

和

分别表示第u个和第v个差值列向量中像素灰度差值的均值，这里1≤u≤Q，，1≤v≤Q；

(4b)计算样本矩阵H的的相关系数矩阵A：

$A=\left[\begin{array}{cccc}{A}_{\mathrm{1,1}}& A_{\mathrm{1,2}}& \&CenterDot;\&CenterDot;\&CenterDot;& A_{1,Q}\\ A_{\mathrm{2,1}}& A_{\mathrm{2,2}}& \&CenterDot;\&CenterDot;\&CenterDot;& A_{2,Q}\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& A_{u,v}& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ A_{Q,1}& A_{Q,2}& \&CenterDot;\&CenterDot;\&CenterDot;& A_{Q,Q}\end{array}\right]$ <4>

其中，相关系数矩阵A中的任意一个元素A_u，v的计算公式如下：

$A_{u,v}=\frac{C_{u,v}}{\sqrt{C_{u,u}{C}_{v,v}}}$ <5>

其中A_u，v为相关系数矩阵A中第u行第v列的计算值，C_u，v为协方差矩阵C中第u行第v列的值，C_u，u为协方差矩阵C中第u行第v列的值，C_v，v为协方差矩阵C中第v行第v列的值。

(4c)定义分解的层数l＝0，1，…，Q-1。当l＝0时，初始化和变量为样本矩阵H，即S⁽⁰⁾＝H，差变量D⁽⁰⁾为空集，并将和变量的下标由集合Ω表示，Ω＝{1，2，…，Q}，将差变量下标集合Φ作为空集，将正交Dirac基B⁽⁰⁾＝[φ_0，1，φ_0，2，φ_0，3，…φ_0，L]作为Q×Q维的单位矩阵；

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

$(\mathrm{\&alpha;},\mathrm{\&beta;})=\arg \underset{u,v\&Element;\mathrm{\&Omega;}}{\max }{A}_{u,v}^{(l-1)}$ <6>

此时仅考虑相关系数矩阵A的上三角部分，u和v分别表示相关系数矩阵A中任意值的行和列，且只在和变量集合Ω内进行；

(4e)计算Jacobi旋转矩阵J：

<7>

其中c_n＝cos(θ_l)，s_n＝sin(θ_l)；

旋转角θ_l通过C^(l)＝J^TC^(l-1)J和

这两个公式计算得到：

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

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

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

(4g)计算样本矩阵H沿当前正交基矩阵B^(l)的方向上的归一化能量矩阵Fn^(l)：

${\mathrm{Fn}}^{\left(l\right)}=[\underset{s\&Element;\mathrm{\&Omega;}}{{\mathrm{Fn}}_s^{\left(l\right)}},\underset{d\&Element;\mathrm{\&Phi;}}{{\mathrm{Fn}}_d^{\left(l\right)}}]=\frac{E{\left|\right|H\&CenterDot;B^{\left(l\right)}\left|\right|}^2}{E{\left|\right|H\left|\right|}^2}$ <9>

其中E||·||²为计算数学期望，

为归一化能量矩阵Fn^(l)中和变量下标集合Ω所对应向量组成的矩阵，

为归一化能量矩阵Fn^(l)中差变量下标集合Φ所对应向量组成的矩阵；

计算和变量下标集合Ω所对应向量组成的矩阵

的能量总和

以及差变量下标集合Φ所对应向量组成的矩阵的能量总和

当满足

将当前分解层数l作为最终分解层数L；

(4h)利用正交基矩阵B^(L)对样本矩阵H进行投影，得到投影矩阵

En＝(H×B^(L)′)·(H×B^(L)′)^T，

提取投影矩阵En中的前Q-L行，得到尺度对应矩阵Em，将尺度对应矩阵Em中各列向量中最大值对应的行标号作为该列向量的类别标记；

(4i)寻找与差值列向量h_ii相同类别标记的q个列向量，其中1≤q≤Q，将这q个列向量组成矩阵R并计算该矩阵R的均值m_i，将均值m_i作为自适应差异图D中像素点i的灰度值。

步骤5，对两幅遥感图像T1和T2中的每个像素点，重复进行步骤2至步骤4，得到自适应差异图D。

步骤6，用Otsu阈值法分别计算差值差异图X和自适应差异图D的Otsu阈值K₁和K₂：

(6a)设差值差异图X的灰度等级为W，差值差异图X中像素取值为[0，…，W]，灰度值为r的像素个数是n_r，则灰度值r的概率p_r＝n_r/Np，Np为差值差异图X的总像素数，灰度值K为阈值将差值差异图X分为两类，即灰度为[0，…，K]的像素构成一类，记为D₀；灰度值为(K，…，W]的像素构成另一类，记为D₁，且0≤K≤W；

(6b)计算整幅差值差异图X图像的灰度均值μ：

μ＝P₀(K)μ₀(K)+P₁(K)μ₁(K)

其中P₀(K)为D₀类的概率，P₁(K)为D₁类的概率，μ₀(K)为D₀类所有像素的灰度均值，μ₁(K)为D₁类所有像素的灰度均值，具体计算公式如下：

$P_0\left(K\right)=\underset{r=1}{\overset{K}{\mathrm{\&Sigma;}}}{p}_r$

P₁(K)＝1-P₀(K)

${\mathrm{\&mu;}}_0\left(K\right)=\frac{1}{P_0\left(K\right)}\underset{r=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{rp}}_r$

${\mathrm{\&mu;}}_1\left(K\right)=\frac{1}{P_1\left(K\right)}\underset{r=K+1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{rp}}_r$

根据差值差异图X的灰度均值μ，计算两类的类间距离平方

${\mathrm{\&sigma;}}_b^2\left(K\right)=P_0\left(K\right){({\mathrm{\&mu;}}_0\left(K\right)-\mathrm{\&mu;})}^2+P_1\left(K\right){({\mathrm{\&mu;}}_1\left(K\right)-\mathrm{\&mu;})}^2$ <10>

(6c)当满足时，K₁即为差值差异图X的Otsu阈值；

(6d)根据步骤(6a)至步骤(6c)，计算自适应差异图D的Otsu阈值K₂。

步骤7，利用Otsu阈值K₁和K₂，按如下融合规则对差值差异图X和自适应差异图D进行融合：

<11>

其中X(f)和D(f)分别为差值差异图X和自适应差异图D空间对应位置f处的像素点灰度值，K₁为差值差异图X的Otsu阈值，K₂为自适应差异图D的Otsu阈值，I(f)为融合后最终差异图I空间对应位置f处的像素点灰度值。

步骤8，根据步骤(6a)至步骤(6c)，计算最终差异图I的Otsu分割阈值K₃，将最终差异图I中像素灰度值大于K₃的像素点标记为1，其余像素点标记为0，该标记结果即为变化检测结果图Z。

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

1.实验数据

本发明的实验数据为一组模拟变化遥感图像和一组真实多时相遥感图像共两组图像，每一组遥感图像都有变化检测结果参考图。模拟原始遥感图像为470×335像素大小、位于英国Feltwell村庄农区的Airborne Thematic Mapper第3波段的图像，模拟变化遥感图像是通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工的嵌入一些变化区域得到的。真实遥感图像为大小为512×512像素的Landsat7Enhanced Thematic Mapper Plus第4波段光谱图像，由2000年4月和2002年5月的墨西哥郊外的两幅图像组成，变化部分是由于火灾破坏了大面积的当地植被所致。

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

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

对比方法1，是Lingcao Huang学者在文章“An Object-based Change DetectionApproach by Integrating Intensity and Texture Differences”中提出的加权融合Sobel边缘纹理图和差值差异图得到概率差异图，利用阈值对融合后的概率差异图进行分割得到变化检测图像。

对比方法2，是Wei He在文章“Application of Euclidean norm in Multi-temporalRemote Sensing Image Change Detection”中提出的利用两幅遥感图像之间差值的欧式距离范数表示差异图，再利用最大熵阈值法提取变化区域。

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

3.实验的内容与分析

实验1，用不同方法对两幅不同时相的模拟变化遥感图像进行变化检测，如图2所示，其中原始的两幅遥感图像分别如图2(a)和图2(b)所示，变化检测参考图如图2(c)所示，对图2(a)和图2(b)的模拟变化遥感图像进行变化检测得到的结果如图3所示，其中图3(a)为对比方法1得到的变化检测结果图，图3(b)为对比方法2得到的变化检测结果图，图3(c)为用本发明方法得到的变化检测结果图。从图3可以看出对比方法1得到的结果图含有大量的噪声点，虚警数非常多，对比方法2和本发明方法得到的结果图在细节上的保持都比较好，但本发明方法得到的结果图最接近于变化检测参考图的图2(c)。

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

表1列出了两组遥感图像变化检测结果的评价指标。

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

从表1的评价指标中可以更为直观的看出，本发明对遥感图像变化检测的各项指标均优于现有检测方法。

Claims (4)

1.一种基于自适应差异图的遥感图像变化检测方法，包括如下步骤：

(1)输入两幅已配准的多时相遥感图像T1和T2，并将图像T1和T2空间位置对应像素的灰度值进行差值计算，得到差值差异图X；

(2)对第一幅图像T1以像素点i为中心取3×3大小的滑窗E_i，对第二幅图像T2同样以像素点i为中心取7×7像素大小搜索窗G_i，i为当前像素点，在搜索窗G_i内逐像素取3×3大小的滑窗，得到滑窗集合{F_j}，F_j为滑窗集合内以j为中心点的一个滑窗，定义滑窗集合{F_j}内的滑窗数目为Q个，此处Q＝25；

(3)分别将第二幅图像T2得到的Q个滑窗与第一幅图像T1的滑窗Ei进行差值计算，并将差值的结果用列向量进行表示，得到样本矩阵H，该样本矩阵H中包含Q个列向量且维数为9×Q，将第二幅图像T2的滑窗F_i与第一幅图像T1的滑窗E_i的差值列向量表示为h_ii；

(4)利用Treelet算法对样本矩阵H进行自适应聚类，得到与差值列向量h_ii相同类别标记的q个列向量，其中1≤q≤Q，将这q个列向量组成矩阵R，并计算该矩阵R的均值m_i，将均值m_i作为自适应差异图D中像素点i的灰度值；

(5)对两幅遥感图像T1和T2中的每个像素点，重复进行步骤2至步骤4，得到自适应差异图D；

(6)用Otsu阈值法分别计算差值差异图X和自适应差异图D的阈值K₁和K₂；

(7)利用Otsu阈值K₁和K₂融合差值差异图X和自适应差异图D，得到最终差异图I；

(8)对最终差异图I利用Otsu阈值法进行分割，得到变化检测的结果图像Z。

2.根据权利要求1所述的遥感图像变化检测方法，其中步骤(4)所述的利用Treelet算法对样本矩阵H进行自适应聚类，按如下步骤进行：

(4a)计算样本矩阵H的协方差系数矩阵C：

$C=\left[\begin{array}{cccc}{C}_{\mathrm{1,1}}& C_{\mathrm{1,2}}& \&CenterDot;\&CenterDot;\&CenterDot;& C_{1,Q}\\ C_{\mathrm{2,1}}& C_{\mathrm{2,2}}& \&CenterDot;\&CenterDot;\&CenterDot;& C_{2,Q}\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& C_{u,v}& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ C_{Q,1}& C_{Q,2}& \&CenterDot;\&CenterDot;\&CenterDot;& C_{Q,Q}\end{array}\right]$

其中，协方差系数矩阵C中的任意一个元素C_uv的计算公式如下

$C_{u,v}=\frac{\underset{k=1}{\overset{9}{\mathrm{\&Sigma;}}}(x_{k,u}-{\stackrel{\&OverBar;}{x}}_u)(x_{k,v}-{\stackrel{\&OverBar;}{x}}_v)}{\sqrt{\underset{k=1}{\overset{9}{\mathrm{\&Sigma;}}}{(x_{k,u}-{\stackrel{\&OverBar;}{x}}_u)}^2\underset{k=1}{\overset{9}{\mathrm{\&Sigma;}}}{(x_{k,v}-{\stackrel{\&OverBar;}{x}}_v)}^2}}$

其中C_u，v表示协方差矩阵C中第u行第v列的计算值，

和

分别表示第u个和第v个差值列向量中像素灰度差值的均值，这里1≤u≤Q，，1≤v≤Q；

(4b)计算样本矩阵H的相关系数矩阵A：

$A=\left[\begin{array}{cccc}{A}_{\mathrm{1,1}}& A_{\mathrm{1,2}}& \&CenterDot;\&CenterDot;\&CenterDot;& A_{1,Q}\\ A_{\mathrm{2,1}}& A_{\mathrm{2,2}}& \&CenterDot;\&CenterDot;\&CenterDot;& A_{2,Q}\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& A_{u,v}& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ A_{Q,1}& A_{Q,2}& \&CenterDot;\&CenterDot;\&CenterDot;& A_{Q,Q}\end{array}\right]$

其中，相关系数矩阵A中的任意一个元素A_u，v的计算公式如下：

$A_{u,v}=\frac{C_{u,v}}{\sqrt{C_{u,u}{C}_{v,v}}}$

其中A_u，v为相关系数矩阵A中第u行第v列的计算值，C_u，v为协方差矩阵C中第u行第v列的值，C_u，u为协方差矩阵C中第u行第v列的值，C_v，v为协方差矩阵C中第v行第v列的值；

(4c)定义分解层数l＝0，1，…，Q-1，当l＝0时，初始化和变量为样本矩阵H，即S⁽⁰⁾＝H，差变量D⁽⁰⁾为空集，并将和变量的下标由集合Ω表示，Ω＝{1，2，…，Q}，将差变量下标集合Φ作为空集，将正交Dirac基B⁽⁰⁾＝[φ_0，1，φ_0，2，φ_0，3，…φ_0，L]作为Q×Q维的单位矩阵；

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

$(\mathrm{\&alpha;},\mathrm{\&beta;})=\arg \underset{u,v\&Element;\mathrm{\&Omega;}}{\max }{A}_{u,v}^{(l-1)}$

此时仅考虑相关系数矩阵A的上三角部分，u和v分别表示相关系数矩阵A中任意值的行和列，且只在和变量集合Ω内进行；

(4e)计算Jacobi旋转矩阵J：

其中cn＝cos(θ_l)，sn＝sin(θ_l)；

旋转角θl通过C(l)＝JTC(l-1)J和

这两个公式计算得到：

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

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

$A_{u,v}^{\left(l\right)}=\frac{C_{u,v}^{\left(l\right)}}{\sqrt{C_{u,u}^{\left(l\right)}{C}_{v,v}^{\left(l\right)}}};$

(4f)定义正交基矩阵B^(l)中第α列和第β列的列向量分别为尺度函数θ_l和细节函数ψ_l，定义当前l层的尺度向量集合{φ_l}是尺度函数φ_l和上一层的尺度向量集合{φ_l-1}的合集，将β差变量从和变量下标集合Ω中去除，即Ω＝Ω\{β}，加入到差变量下标集合Φ中，即Φ＝{β}，此时B(l)＝[φ_l，1，...φ_l，Q-l，ψ₁，...Ψ_l]，其中φ_l，1，...φ_l，Q-l∈{φ_l}；

(4g)计算样本矩阵H沿当前正交基矩阵B^(l)方向上的归一化能量矩阵Fn^(l)：

${\mathrm{Fn}}^{\left(l\right)}=[\underset{s\&Element;\mathrm{\&Omega;}}{{\mathrm{Fn}}_s^{\left(l\right)}},\underset{d\&Element;\mathrm{\&Phi;}}{{\mathrm{Fn}}_d^{\left(l\right)}}]=\frac{E{\left|\right|H\&CenterDot;B^{\left(l\right)}\left|\right|}^2}{E{\left|\right|H\left|\right|}^2}$

其中E||·||²为计算数学期望，为归一化能量矩阵Fn^(l)中和变量下标集合Ω所对应向量组成的矩阵，

为归一化能量矩阵Fn^(l)中差变量下标集合Φ所对应向量组成的矩阵；

计算和变量下标集合Ω所对应向量组成的矩阵的能量总和

以及差变量下标集合Φ所对应向量组成的矩阵

的能量总和

当满足 $\underset{s\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{\mathrm{Fn}}_s^{\left(l\right)}-\underset{s\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{\mathrm{Fn}}_s^{(l-1)}<\underset{d\&Element;\mathrm{\&Phi;}}{\mathrm{\&Sigma;}}{\mathrm{Fn}}_d^{\left(l\right)}$ 时，将当前分解层数l作为最终分解层数L；

(4h)利用正交基矩阵B^(L)对样本矩阵H进行投影，得到投影矩阵：

En＝(H×B^(L)′)·(H×B^(L)′)T，

提取投影矩阵En中的前Q-L行，得到尺度对应矩阵Em，将尺度对应矩阵Em中各列向量中最大值对应的行标号作为该列向量的类别标记；

(4i)寻找与差值列向量h_ii相同类别标记的q个列向量，其中1≤q≤Q，将这q个列向量组成矩阵R，并计算该矩阵R的均值m_i，将均值m_i作为自适应差异图D中像素点i的灰度值。

3.根据权利要求1所述的遥感图像变化检测方法，其中步骤(7)所述的利用Otsu阈值K₁和K₂融合差值差异图X和自适应差异图D，按如下融合规则进行：

其中X(f)和D(f)分别为差值差异图X和自适应差异图D空间对应位置f处的像素点灰度值，K₁为差值差异图X的Otsu阈值，K₂为自适应差异图D的Otsu阈值，I(f)为融合后最终差异图I空间对应位置f处的像素点灰度值。

4.根据权利要求1所述的遥感图像变化检测方法，其中步骤(8)所述的对最终差异图I利用Otsu阈值法进行分割，按如下步骤进行：

(8a)设最终差异图I的灰度等级为W，最终差异图I中像素取值为[0，…，W]，灰度值为r的像素个数是n_r，则灰度值r的概率p_r＝n_r/Np，Np为最终差异图I的总像素数，灰度值K为阈值将最终差异图I分为两类，即灰度为[0，…，K]的像素构成一类，记为D₀；灰度值为(K，…，W]的像素构成另一类，记为D₁，且0≤K≤W；

(8b)计算整幅最终差异图I的灰度均值μ：

μ＝P₀(K)μ₀(K)+P₁(K)μ₁(K)

其中P₀(K)为D₀类的概率，P₁(K)为D₁类的概率，μ₀(K)为D₀类所有像素的灰度均值，μ₁(K)为D₁类所有像素的灰度均值，具体公式如下：

$P_0\left(K\right)=\underset{r=1}{\overset{K}{\mathrm{\&Sigma;}}}{p}_r$

P₁(K)＝1-P₀(K)

${\mathrm{\&mu;}}_0\left(K\right)=\frac{1}{P_0\left(K\right)}\underset{r=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{rp}}_r$

${\mathrm{\&mu;}}_1\left(K\right)=\frac{1}{P_1\left(K\right)}\underset{r=K+1}{\overset{W}{\mathrm{\&Sigma;}}}{\mathrm{rp}}_r$

根据最终差异图I的灰度均值μ，计算两类的类间距离平方

${\mathrm{\&sigma;}}_b^2\left(K\right)=P_0\left(K\right){({\mathrm{\&mu;}}_0\left(K\right)-\mathrm{\&mu;})}^2+P_1\left(K\right){({\mathrm{\&mu;}}_1\left(K\right)-\mathrm{\&mu;})}^2$

(8c)当满足

时，K3即为最终差异图I的Otsu算法分割阈值；

(8d)将最终差异图I中像素灰度值大于K₃的所有像素点标记为1，其余像素点标记为0，该标记结果即为变化检测结果图Z。

Images