CN102360500A - 基于Treelet曲波域去噪的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于Treelet曲波域去噪的遥感图像变化检测方法，其实现步骤为：(1)读入数据；(2)中值滤波；(3)构造对数差异图像；(4)构造绝对值差异图像；(5)快速离散曲波分解；(6)曲波变换系数分类；(7)Fine尺度层置零；(8)Detail尺度层去噪；(9)曲波逆变换；(10)计算变化比例阈值；(11)分类；(12)获得变化检测结果图。本发明对噪声具有较好的鲁棒性，能够较好的保持变化区域的边缘信息，减少伪变化信息，具有较高的检测精度，可用于灾情监测、森林覆盖率评估、城市规划等领域。

Description

基于Treelet曲波域去噪的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，更进一步涉及一种基于Treelet曲波域去噪的遥感图像变化检测方法。该方法可应用于自然灾害中的灾情监测和评估、森林覆盖率的监测和评估、城市规划等领域。

背景技术

变化检测是通过观察不同时刻的一个物体或现象，鉴定其状态差异的过程。随着遥感技术和信息技术的发展，遥感图像变化检测已经成为当前遥感图像分析研究的一个重要方向。

遥感图像变化检测方法可以分为两类：先分类后比较法和先比较后分类法。先分类后比较法的优点在于可以克服由于大气、传感器、季节、地面环境、多时相遥感图像的分辨率等因素带来的不同时相图像间的差异。但该类方法存在误差累计问题，并且分类本身存在病态分割问题，因此会导致最终的变化检测精度不高。先比较后分类法的优势在于直观、简单、易行，没有先分类后比较法所存在的误差累计问题，是现在普遍应用的一类多时相遥感图像变化检测方法。

基于多尺度几何分析的遥感图像变化检测方法能够在不同尺度上进行检测，并通过最优尺度或者多尺度的融合确定最终的变化检测结果，在一定程度上避免了单一尺度下对边缘信息的保持和变化区域信息的检测难以兼顾的缺点。

西安电子科技大学在其专利申请“基于多尺度积和主成分分析的SAR图像变化检测方法”(专利申请号：200910023637.1，公开号：CN101634709)中提出了一种多尺度积去噪和主成分分析融合的SAR图像变化检测方法。该方法虽然能够减弱图像误配准的影响，但仍存在的不足是，该方法仅利用当前尺度层和其下一个尺度层进行去噪，没有考虑当前尺度层和其上一个尺度层的关系，使得该方法的去噪结果不够理想，对噪声信号较为敏感。此外，由于SAR图像并非完全服从广义高斯模型，因此基于广义高斯模型的最小错误率阈值法会使得检测结果中存在较多的伪变化信息，降低了遥感图像变化检测的精度。

2010年黄世奇等在文献“基于小波变换的多时相SAR图像变化检测技术”(测绘学报，2010，39(2)：180-186.)中提出了一种基于可靠尺度融合的遥感图像变化检测方法。该方法通过局部变化系数和全局变化系数的关系确定可靠尺度，并采用加入权重的可靠尺度特征级融合策略进行变化区域检测。该方法虽然能够对变化区域类进行了进一步的类型分析，得到变化区域减弱类和变化区域增强类，但仍存在的不足是，该方法在进行可靠尺度特征级融合时仅考虑了图像的低频信息，因此变化区域的边缘信息不能够得到较好的保持。此外，由于该方法使用EM算法进行分类，使得检测结果中存在较多的伪变化信息，影响了遥感图像变化检测的精度。

发明内容

本发明针对上述现有技术存在的不足，提出了一种基于Treelet曲波域去噪的遥感图像变化检测方法。本发明对噪声信号具有较好的鲁棒性，能够较好的保持变化区域的边缘信息，减少伪变化信息，具有较高的变化检测精度。

本发明实现上述目的的思路是：在对读入的数据进行中值滤波、构造差异图像、快速离散曲波分解和曲波变换系数分类处理之后，使用Treelet变换自适应的对Detail尺度层进行去噪，再利用曲波域低频变化图像及其边缘图计算变化比例阈值，最后进行图像分类和图像融合。

本发明的具体步骤包括如下：

(1)读入同一地区不同时刻获取的两幅遥感图像。

(2)中值滤波

2a)确定正方形窗口：选取步骤(1)中的一幅遥感图像，以该图像中的某一像素点为中心，选取一个大小为Nu×Nu的正方形窗口，其中，Nu为奇数；

2b)确定滤波值：将正方形窗口中全部像素点的灰度值按照由大到小的顺序排列，组成一个灰度序列，选取位于灰度序列中间位置的灰度值作为滤波值；

2c)滤波：用滤波值替代步骤2a)中像素点的灰度值；

2d)重复步骤2a)至步骤2c)，直至处理完图像中的全部像素点；

2e)按照步骤2a)至步骤2d)，对步骤(1)中的另外一幅遥感图像进行处理，得到滤波后的两幅图像。

(3)构造对数差异图像

3a)对步骤(2)中滤波后的两幅图像分别进行对数运算，得到两幅对数图像；

3b)对两幅对数图像对应像素点进行减法运算，并对减法运算的结果取绝对值，得到一幅对数差异图像。

(4)构造绝对值差异图像

对步骤(2)中滤波后的两幅图像对应像素点进行减法运算，并对减法运算的结果取绝对值，得到一幅绝对值差异图像。

(5)对步骤(3)中的对数差异图像进行Sc层快速离散曲波分解。

(6)根据曲波(Curvelet)变换域频率的大小，对步骤(5)中的曲波变换系数进行分类。

(7)将Fine尺度层的曲波变换系数全部赋值为零。

(8)Detail尺度层去噪

8a)水平方向去噪：对Detail尺度层的任一方向子带，应用Treelet变换自适应的计算水平方向去噪阈值，并进行软阈值去噪，得到该方向子带的水平方向去噪结果；

8b)垂直方向去噪：对步骤8a)中选取的方向子带，应用Treelet变换自适应的计算垂直方向去噪阈值，并进行软阈值去噪，得到该方向子带的垂直方向去噪结果；

8c)平均加权求和：将水平方向去噪结果和垂直方向去噪结果进行平均加权求和，得到该方向子带的最终去噪结果；

8d)重复步骤8a)、8b)和8c)，直至处理完Detail尺度层的所有方向子带。

(9)曲波逆变换

对Coarse尺度层的Curvelet变换系数、经步骤(7)处理后的Fine尺度层的Curvelet变换系数和经步骤(8)处理后的Detail尺度层的Curvelet变换系数进行逆Curvelet变换，得到重构图像。

(10)计算变化比例阈值

10a)低频差异图像：将Fine尺度层和Detail尺度层的Curvelet变换系数全部置零，并与Coarse尺度层的Curvelet变换系数进行逆Curvelet变换，得到仅包含低频信息的低频差异图像；

10b)低频变化图像：对低频差异图像采用最大类间方差(Otsu)法分类，得到低频变化图像；

10c)最终边缘图：分别采用Canny、Sobel、Prewitt边缘算子提取低频变化图像的边缘，得到三幅边缘图，通过逻辑或运算将此三幅边缘图合成一幅最终边缘图；

10d)计算变化比例阈值：统计最终边缘图中的边缘像素点个数和低频变化图像中变化像素点的个数，计算变化比例阈值。

(11)分类

11a)通过变化比例阈值对步骤(4)获得的绝对值差异图像进行分类，得到一幅分类图像；

11b)通过变化比例阈值对步骤(9)获得的重构图像进行分类，得到另一幅分类图像。

(12)获得变化检测结果图

对步骤11a)获得的分类图像和步骤11b)获得的分类图像进行融合，得到变化检测结果图。

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

第一，本发明由于采用了Treelet变换自适应的对Detail尺度层进行去噪，克服了现有技术对噪声较为敏感的缺点，使得本发明的检测结果不容易受噪声信号的干扰，对噪声具有较好的鲁棒性。

第二，本发明根据Curvelet域低频变化图像及其边缘图计算变化比例阈值，克服了现有技术存在较多伪变换信息的缺点，使得本发明的检测精度得到了提高。

第三，本发明通过图像融合获得变化检测结果图像，克服了现有技术对变化区域边缘信息保持不理想的缺点，使得本发明的变化区域边缘信息得到了更加准确、全面的保持。

附图说明

图1为本发明的流程图；

图2为本发明的仿真效果图。

具体实施方式

下面结合附图1对本发明的步骤做进一步的详细描述。

步骤1，读入同一地区不同时刻获取的两幅遥感图像。

步骤2，中值滤波。

2a)确定正方形窗口：选取步骤1中的一幅遥感图像，以该图像中的某一像素点为中心，选取一个Nu×Nu的正方形窗口，其中，Nu为奇数，本发明实施例中选取一个3×3的正方形窗口。

2b)确定滤波值：将正方形窗口中全部像素点的灰度值按照由大到小的顺序排列，组成一个灰度序列，选取位于灰度序列中间位置的灰度值作为滤波值。

2c)滤波：将滤波值替代步骤2a)中像素点的灰度值。

2d)重复步骤2a)至步骤2c)，直至处理完图像中的全部像素点。

2e)按照步骤2a)至步骤2d)，对步骤(1)中的另一幅图像进行处理，得到滤波后的两幅图像。

步骤3，构造对数差异图像。

3a)对步骤2中滤波后的两幅图像分别进行对数运算，得到两幅对数图像。其对数运算公式为：

I₃(m，n)＝log(I₁(m，n)+1)

I₄(m，n)＝log(I₂(m，n)+1)

其中，I₃和I₄为两幅对数图像(大小为M×N)，I₁和I₂为滤波后的两幅图像，m和n为图像的行序号和列序号，m＝0，1，2，…，M-1，n＝0，1，2，…，N-1。

3b)对两幅对数图像对应像素点进行减法运算，并对减法运算的结果取绝对值，得到一幅对数差异图像L。

L(m，n)＝|I₃(m，n)-I₄(m，n)|

其中，L为对数差异图像，I₃和I₄为两幅对数图像(大小为M×N)，m和n为图像的行序号和列序号，m＝0，1，2，…，M-1，n＝0，1，2，…，N-1。

步骤4，构造绝对值差异图像。

对步骤2中滤波后的两幅图像对应像素点进行减法运算，并对减法运算的结果取绝对值，得到一幅绝对值差异图像A。

A(m，n)＝|I₁(m，n)-I₂(m，n)|

其中，A为绝对值差异图像，I₁和I₂为滤波后的两幅图像(大小为M×N)，m和n为图像的行序号和列序号，m＝0，1，2，…，M-1，n＝0，1，2，…，N-1。

步骤5，快速离散曲波分解。

对步骤3中的对数差异图像L进行Sc层快速离散曲波分解。

其中，C′(sc，d)为曲波变换系数，sc为分解层参数，d为方向参数，L为对数差异图像，

为曲波函数，q为方向子带内的位置参数。

本发明实施例中，对数差异图像L进行5层快速离散曲波分解。对于第sc＝1分解层，方向参数d＝1，即无方向信息；对于第sc＝2分解层，方向参数d＝{1，2，…，16}，即包含16个方向子带；对于第sc＝3分解层，方向参数d＝{1，2，…，16}，即包含16个方向子带；对于第sc＝4分解层，方向参数d＝{1，2，…，32}，即包含32个方向子带；对于第sc＝5分解层，方向参数d＝1，即无方向信息。

步骤6，曲波变换系数分类。

根据曲波变换域频率的大小，对步骤5中的曲波(Curvelet)变换系数进行分类。将第Sc分解层的低频系数归类为Coarse尺度层，第Sc-1、Sc-2、…、2分解层的高频系数归类为Detail尺度层，第1分解层的高频系数归类为Fine尺度层。

在本发明实施例中，将第5分解层的低频系数归类为Coarse尺度层；将第4、3、2分解层的高频系数归类为Detail尺度层，共包含64个不同的方向子带；将第1分解层的高频系数归类为Fine尺度层。

步骤7，Fine尺度层置零。

将Fine尺度层的Curvelet变换系数全部赋值为零。

步骤8，Detail尺度层去噪。

8a)水平方向去噪。对Detail尺度层的任一方向子带，应用Treelet变换自适应的计算其水平方向去噪阈值，并进行软阈值去噪，得到该方向子带的水平方向去噪结果。

(a1)第l＝0分解层数据初始化。

将Detail尺度层任一方向子带(大小为K×P)的每一行视为一个P维行向量，则共有K个行向量。记

为水平方向初始化结果，和向量下标集δ_H＝{1，2，…，K}，正交基矩阵为B₀＝[Φ_0，1，Φ_0，2，…，Φ_0，K]。其中，B₀是一个K×K的单位矩阵，上标T表示转置。

计算

的协方差矩阵

其中，

为求数学期望，h＝{1，2，…，K}和v＝{1，2，…，K}分别为协方差矩阵的位置索引。

计算的相似度矩阵

相似度矩阵

中的每个元素均为相关系数

${\mathrm{\Θ}}_{H(h,v)}^{\left(0\right)}={\mathrm{\ρ}}_{\mathrm{hv}}^{\left(0\right)}=\left|\frac{{\mathrm{\Σ}}_{\mathrm{hv}}^{\left(0\right)}}{\sqrt{{\mathrm{\Σ}}_{\mathrm{hh}}^{\left(0\right)}{\mathrm{\Σ}}_{\mathrm{vv}}^{\left(0\right)}}}\right|$

(a2)第l＝1，2，…，K-1分解层Treelet变换。

由相似度矩阵

找出最相似的两个向量：

$(\mathrm{\α},\mathrm{\β})=\arg \underset{h,v\∈{\mathrm{\δ}}_H}{\max }{{\mathrm{\Θ}}_H^{(l-1)}}_{(h,v)}$

其中，h＜v。

对上述最相似的两个向量进行局部PCA变换：

其中，c＝cos(θ_l)，s＝sin(θ_l)。旋转角θ_l由以下三个式子计算得到：

|θ_l|≤π/4

${\mathrm{\Σ}}_H^{\left(l\right)}=J^T{\mathrm{\Σ}}_H^{(l-i)}J$

${\mathrm{\Σ}}_{\mathrm{\α\β}}^{\left(l\right)}={\mathrm{\Σ}}_{\mathrm{\β\α}}^{\left(l\right)}=0$

利用Jacobi旋转矩阵J更新第l分解层的基矩阵B_l＝B_l-1J＝[Φ_l，1，Φ_l，2，…，Φ_l，K]和第l分解层的 $C_H^{\left(l\right)}=J^T{C}_H^{(l-1)}.$

经过Jacobi旋转后，使和

满足关系式定义第l分解层的和向量与差向量分别为

和并定义第l分解层的尺度向量Ω_l和细节向量Ψ_l分别为基矩阵B_l的第α列和第β列。去除和向量下标集中的差向量下标，即δ_H＝δ_H\{β}。则第l分解层的

的表达式为：

$C_H^{\left(l\right)}=\underset{z_1=1}{\overset{K-l}{\mathrm{\Σ}}}{s}_{l,z_1}{\mathrm{\Ω}}_{l,z_1}+\underset{z_1=1}{\overset{l}{\mathrm{\Σ}}}{d}_{z_1}{\mathrm{\Ψ}}_{z_1}$

(a3)计算水平方向去噪阈值。

计算第l分解层第i₁个尺度向量的归一化能量

$\mathrm{\ϵ}\left({\mathrm{\Φ}}_{i,i_1}\right)=\frac{{\mathrm{\Σ}}_{j_1=1}^K{|{\mathrm{\Φ}}_{l,i_1}\·x_{j_1}^{\left(0\right)}|}^2}{{\mathrm{\Σ}}_{j_1=1}^K{\left|\right|x_{j_1}^{\left(0\right)}\left|\right|}^2}$

其中，

为第l分解层第i₁个尺度向量的归一化能量，为Treelet变换第l分解层基矩阵的第i₁个尺度向量，K为该方向子带的行向量个数，为该方向子带的第j₁个行向量，i₁＝1，2，…，K，l＝1，2，…，K-1。

计算第l分解层的总能量

$E_H^{\left(l\right)}={\mathrm{\Σ}}_{i_1=1}^{K-l}\mathrm{\ϵ}\left({\mathrm{\Φ}}_{l,i_1}\right)$

其中，

为第l分解层的总能量，K为该方向子带的行向量个数，为第l分解层第i₁个尺度向量的归一化能量。

选择最优分解层lb。

$\mathrm{lb}=\arg \underset{l=\mathrm{1,2},...,K-3}{\max }\frac{E_H^{(l+1)}-E_H^{(l+2)}}{E_H^{\left(l\right)}-E_H^{(l+1)}}$

其中，lb为最优分解层，l为分解层数，K为该方向子带的行向量个数，

为第l分解层的总能量，

为第l+1分解层的总能量，

为第l+2分解层的总能量，max为取最大值函数。

构造水平方向最优投影基矩阵PB。

PB＝[Φ_lb，1，…，Φ_lb，f，…，Φ_lb，(K-lb)]

其中，PB为水平方向最优投影基矩阵，lb为最优分解层，Φ_lb，1为第lb分解层基矩阵B_lb的第1个尺度向量，Φ_lb，f为B_lb的第f个尺度向量，Φ_lb，(K-lb)为B_lb的第K-lb个尺度向量。

水平方向主要成分的计算公式为：

$C_1=C_H^{\left(0\right)}\·\mathrm{PB}$

其中，C₁为水平方向主要成分，

为该子带的水平方向初始化结果，PB为水平方向最优投影基矩阵。

水平方向去噪阈值的计算公式为：

T₁＝mean(C₁)

其中，T₁为水平方向去噪阈值，C₁为水平方向主要成分，mean为均值运算。

(a4)软阈值去噪。

对

进行软阈值去噪，得到水平方向去噪结果C_re1。

$C_{\mathrm{rel}}=\left\{\begin{array}{cc}\mathrm{sign}\left(C_H^{\left(0\right)}\right)\·\left(\right|C_H^{\left(0\right)}|-T_1),& \left|C_H^{\left(0\right)}\right|\≥T_1\\ 0,& \mathrm{else}\end{array}\right.$

其中，C_re1为水平方向去噪结果，sign(·)为符号函数，

为该子带的水平方向初始化结果，T₁为水平方向去噪阈值。

8b)垂直方向去噪。对步骤8a)中选取的方向子带，应用Treelet变换自适应的计算其垂直方向去噪阈值，并进行软阈值去噪，得到该方向子带的垂直方向去噪结果。

(b1)第g＝0分解层数据初始化。

将步骤8a)中选取的方向子带的每一列视为一个K维列向量，则共有P个列向量。记

为垂直方向初始化结果，和向量下标集δ_V＝{1，2，…，P}，正交基矩阵为B′₀＝[Φ′_0，1，Φ′_0，2，…，Φ′_0，P]。其中，B′₀是一个P×P的单位矩阵。

计算

的协方差矩阵

其中，

为求数学期望，w＝{1，2，…，P}和z＝{1，2，…，P}分别为协方差矩阵

的位置索引。

计算

的相似度矩阵

相似度矩阵中的每个元素均为相关系数

${\mathrm{\Θ}}_{V(w,z)}^{\left(0\right)}={{\mathrm{\ρ}}^{\′}}_{\mathrm{wz}}^{\left(0\right)}=\left|\frac{{{\mathrm{\Σ}}^{\′}}_{\mathrm{wz}}^{\left(0\right)}}{\sqrt{{{\mathrm{\Σ}}^{\′}}_{\mathrm{ww}}^{\left(0\right)}{{\mathrm{\Σ}}^{\′}}_{\mathrm{zz}}^{\left(0\right)}}}\right|$

(b2)第g＝1，2，…，P-1分解层Treelet变换。

由相似度矩阵

找出最相似的两个向量：

$(\mathrm{\η},\mathrm{\λ})=\arg \underset{w,z\∈{\mathrm{\δ}}_V}{\max }{{\mathrm{\Θ}}_V^{(g-1)}}_{(w,z)}$

其中，w＜z。

对上述最相似的两个向量进行局部PCA变换：

其中，x＝cos(θ_g)，y＝sin(θ_g)。旋转角θ_g由以下三个式子计算得到：

|θ_g|≤π/4

${\mathrm{\Σ}}_V^{\left(g\right)}=Y^T{\mathrm{\Σ}}_V^{(g-1)}Y$

${{\mathrm{\Σ}}^{\′}}_{\mathrm{n\λ}}^{\left(g\right)}={{\mathrm{\Σ}}^{\′}}_{\mathrm{\λ\η}}^{\left(g\right)}=0$

更新第g分解层的基矩阵B′_g＝B′_g-1Y＝[Φ′_g，1，Φ′_g，2，…，Φ′_g，P]和第g分解层的 $C_V^{\left(g\right)}=Y^T{C}_V^{(g-1)}.$

经过Jacobi旋转后，使

和

满足关系式

定义第g分解层的和向量与差向量分别为

和

并定义第g分解层的尺度向量Ω′_g和细节向量Ψ′_g分别为基矩阵B′_g的第η列和第λ列。去除和向量下标集中的差向量下标，即δ_V＝δ_V\{λ}。则第g分解层的的表达式为：

$C_V^{\left(g\right)}=\underset{z_2=1}{\overset{P-g}{\mathrm{\Σ}}}{s^{\′}}_{g,z_2}{{\mathrm{\Ω}}^{\′}}_{g,z_2}+\underset{z_2=1}{\overset{g}{\mathrm{\Σ}}}{d^{\′}}_{z_2}{{\mathrm{\Ψ}}^{\′}}_{z_2}$

(b3)计算垂直方向去噪阈值。

计算第g分解层第i₂个尺度向量的归一化能量

${\mathrm{\ϵ}}_V\left({{\mathrm{\Φ}}^{\′}}_{g,i_2}\right)=\frac{{\mathrm{\Σ}}_{j_2=1}^P{|{{\mathrm{\Φ}}^{\′}}_{g,i_2}\·{x^{\′}}_{j_2}^{\left(0\right)}|}^2}{{\mathrm{\Σ}}_{j_2=1}^P{\left|\right|{x^{\′}}_{j_2}^{\left(0\right)}\left|\right|}^2}$

其中，

为第g分解层第i₂个尺度向量的归一化能量，

为Treelet变换第g分解层基矩阵的第i₂个尺度向量，P为该方向子带的行向量个数，为该方向子带的第j₂个行向量，i₂＝1，2，…，P，g＝1，2，…，P-1。

计算第g分解层的总能量

$E_V^{\left(g\right)}={\mathrm{\Σ}}_{i_2=1}^{P-g}{\mathrm{\ϵ}}_V\left({{\mathrm{\Φ}}^{\′}}_{g,i_2}\right)$

其中，

为第g分解层的总能量，P为该方向子带的列向量个数，

为第g分解层第i₂个尺度向量的归一化能量。

选择最优分解层lb′

${\mathrm{lb}}^{\′}=\arg \underset{g=\mathrm{1,2},...,P-3}{\max }\frac{E_V^{(g+1)}-E_V^{(g+2)}}{E_V^{\left(g\right)}-E_V^{(g+1)}}$

其中，lb′为最优分解层，g为分解层数，P为该方向子带的列向量个数，

为第g分解层的总能量，

为第g+1分解层的总能量，

为第g+2分解层的总能量，max为取最大值函数。

构造水平方向最优投影基矩阵PB′。

PB′＝[Φ′_lb′，1，…，Φ′_lb′，k，…，Φ′_{lb′，(P-lb′)}]

其中，PB′为垂直方向最优投影基矩阵，lb′为最优分解层，Φ′_lb′，1为第lb分解层基矩阵B′_lb′的第1个尺度向量，Φ′_lb′，k为B′_lb′的第k个尺度向量，Φ′_lb′，_(P-lb′)为B′_lb′的第P-lb′个尺度向量。

垂直方向主要成分的计算公式为：

$C_2=C_V^{\left(0\right)}\·{\mathrm{PB}}^{\′}$

其中，C₂为垂直方向主要成分，

为该子带的垂直方向初始化结果，PB′为垂直方向最优投影基矩阵。

垂直方向去噪阈值的计算公式为：

T₂＝mean(C₂)

其中，T₂为垂直方向去噪阈值，C₂为垂直方向主要成分，mean为均值运算。

(b4)软阈值去噪。

对

进行软阈值去噪，得到垂直方向去噪结果C_re2。

$C_{\mathrm{re}2}=\left\{\begin{array}{cc}\mathrm{sign}\left(C_V^{\left(0\right)}\right)\·\left(\right|C_V^{\left(0\right)}|-T_2),& \left|C_V^{\left(0\right)}\right|\≥T_2\\ 0,& \mathrm{else}\end{array}\right.$

其中，C_re2为垂直方向去噪结果，sign(·)为符号函数，为该子带的垂直方向初始化结果，T₂为垂直方向去噪阈值。

8c)平均加权求和。将水平方向去噪结果C_re1和垂直方向去噪结果C_re2进行平均加权求和，得到该方向子带的最终去噪结果C_dn。

$C_{\mathrm{dn}}=\frac{1}{2}(C_{\mathrm{re}1}+C_{\mathrm{re}2})$

其中，C_dn为Detail尺度层任一方向子带的最终去噪结果，C_re1为该方向子带的水平方向去噪结果，C_re2为该方向子带垂直方向的去噪结果。

8d)重复步骤8a)、8b)和8c)，直至处理完Detail尺度层的所有方向子带。

步骤9，曲波逆变换。

对Coarse尺度层的Curvelet变换系数、经步骤7处理后的Fine尺度层的Curvelet变换系数和经步骤8处理后的Detail尺度层的Curvelet变换系数进行逆Curvelet变换，得到重构图像。

步骤10，计算变化比例阈值。

10a)低频差异图像。将Fine尺度层和Detail尺度层的Curvelet变换系数全部置零，并与Coarse尺度层的Curvelet变换系数进行逆Curvelet变换，得到仅包含低频信息的低频差异图像。

10b)低频变化图像。对低频差异图像采用最大类间方差(Otsu)法分类，得到低频变化图像。

10c)最终边缘图。分别采用Canny、Sobel、Prewitt边缘算子提取低频变化图像的边缘，得到三幅边缘图，通过逻辑或运算将此三幅边缘图合成一幅最终边缘图。

10d)计算变化比例阈值。

统计最终边缘图中的边缘像素点总个数和低频变化图像中的变化像素点总个数，计算变化比例阈值P。

其计算变化比例阈值P的公式为：

P＝P₁+P₂

其中，P为变化比例阈值，P₁为最终边缘图中的边缘像素点总个数，P₂为低频变化检测图像中的变化像素点总个数。

步骤11，分类。

11a)将步骤4获得的绝对值差异图像中所有像素点的灰度值按照由大到小的顺序排列，将前P个最大灰度值对应的像素点视为变化类，赋值为1；其他像素点视为非变化类，赋值为0，得到分类图像X₁，其中，P为变化比例阈值。

11b)将步骤9获得的重构图像中所有像素点的灰度值按照由大到小的顺序排列，将前P个最大灰度值对应的像素点视为变化类，赋值为1；其他像素点视为非变化类，赋值为0，得到分类图像X₂，其中，P为变化比例阈值。

步骤12，获得变化检测结果图。

对步骤11a)获得的分类图像和步骤11b)获得的分类图像进行融合，得到变化检测结果图M_P。

其融合方式为：

M_P(m，n)＝X₁(m，n)∩X₂(m，n)

其中，M_P(m，n)为变化检测结果图，X₁和X₂为两幅分类图像(大小为M×N)，m和n为图像的行序号和列序号，m＝0，1，2，…，M-1，n＝0，1，2，…，N-1，∩表示逻辑与运算，即当且仅当X₁(m，n)和X₂(m，n)的灰度值都为1时，将M_P(m，n)的灰度值赋值为1，否则将M_P(m，n)的灰度值赋值为0。

下面结合附图2对本发明的仿真效果做进一步的描述。

1.仿真条件

本发明的仿真是在主频2.5GHZ的Pentium Dual_Core CPU E5200、内存1.98GB的硬件环境和MATLAB R2008a的软件环境下进行的。

本发明实施例中，采用3×3的中值滤波和5层快速离散曲波分解。对于快速离散曲波分解的第1分解层，其方向参数d＝1，即无方向信息；对于第2分解层，其方向参数d＝{1，2，…，16}，即包含16个方向子带；对于第3分解层，其方向参数d＝{1，2，…，16}，即包含16个方向子带；对于第4分解层，其方向参数d＝{1，2，…，32}，即包含32个方向子带；对于第5分解层，其方向参数d＝1，即无方向信息。

2.仿真内容

本发明仿真实验所用数据为两组真实遥感数据集。第一组真实遥感数据集是1995年9月和1996年7月意大利撒丁岛Mulargia湖泊区域的两幅Landsat5TM+5波段图像，两幅图像的大小均为256×384像素，它们之间发生的变化是由于湖水水位上升引起的，包括7613个变化像素和90691个非变化像素。第二组真实遥感数据集是1994年8月和1994年9月意大利Elba岛西部地区的两幅Landsat-5TM第4波段光谱图像，两幅图像的大小均为384×320像素，它们之间发生的变化是由于森林火灾破坏大量植被引起的，包含2415个变化像素和120465个非变化像素。

本发明采用虚警数、漏检数和总错误数三个指标来评价变化检测方法的好坏。

3.仿真效果分析

本发明使用Treelet变换自适应的对Detail尺度层进行去噪，使得去噪结果更为理想。为了验证使用Treelet变换的有效性，本发明进行了对比实验。将采用Treelet变换获得Detail尺度层各个方向子带去噪阈值的方法与直接采用Detail尺度层各个方向子带均值作为去噪阈值的方法(S_Denoise法)进行了对比。

本发明根据Curvelet域低频变化检测图像及其边缘图计算变化比例阈值，可以减少检测结果中的伪变化信息。为了验证变化比例阈值选取策略的正确性和有效性，本发明进行了对比实验。将采用Otsu法对绝对值差异图像和重构图像进行分类的方法简记为Otsu法。

图2为本发明的仿真效果图。其中，图2(a)为第一组真实遥感数据集Otsu法的效果图，图2(b)为第一组真实遥感数据集S_Denoise法的效果图，图2(c)为第一组真实遥感数据集本发明的效果图，图2(d)为第二组真实遥感数据集Otsu法的效果图，图2(e)为第二组真实遥感数据集S_Denoise法的效果图，图2(f)为第二组真实遥感数据集本发明的效果图。

表1.变化检测结果性能评价

从表1中可以看出，本发明对第一组真实遥感数据集的总错误数比Otsu法的总错误数少57个像素点，比S_Denoise法的总错误数少11个像素点；对第一组真实遥感数据集的虚警数比Otsu法的虚警数多483个像素点，比S_Denoise法的虚警数少34个像素点；对第一组真实遥感数据集的漏检数比Otsu法的漏检数少534个像素点，比S_Denoise法的漏检数多23个像素点。本发明对第二组真实遥感数据集的总错误数比Otsu法的总错误数少131个像素点，比S_Denoise法的总错误数少45个像素点；对第二组真实遥感数据集的虚警数比Otsu法的虚警数多135个像素点，比S_Denoise法的虚警数少20个像素点；对第二组真实遥感数据集的漏检数比Otsu法的漏检数多4个像素点，比S_Denoise法的漏检数少25个像素点。由此可以看出，本发明能够较为全面、准确的检测出变化信息，减少伪变化信息，具有较高的检测精度。从两组实验数据集的效果图中可以看出，与S_Denoise法和Otsu法相比，本发明能够较好的保持变化区域的边缘信息，误检的孤立像素点也是较少的。

Claims (9)

1.基于Treelet曲波域去噪的遥感图像变化检测方法，包括如下步骤：

(1)读入同一地区不同时刻获取的两幅遥感图像；

(2)中值滤波

2a)确定正方形窗口：选取步骤(1)中的一幅遥感图像，以该图像中的某一像素点为中心，选取一个大小为Nu×Nu的正方形窗口，其中，Nu为奇数；

2b)确定滤波值：将正方形窗口中全部像素点的灰度值按照由大到小的顺序排列，组成一个灰度序列，选取位于灰度序列中间位置的灰度值作为滤波值；

2c)滤波：用滤波值替代步骤2a)中像素点的灰度值；

2d)重复步骤2a)至步骤2c)，直至处理完图像中的全部像素点；

2e)按照步骤2a)至步骤2d)，对步骤(1)中的另外一幅遥感图像进行处理，得到滤波后的两幅图像；

(3)构造对数差异图像

3a)对步骤(2)中滤波后的两幅图像分别进行对数运算，得到两幅对数图像；

3b)对两幅对数图像对应像素点进行减法运算，并对减法运算的结果取绝对值，得到一幅对数差异图像；

(4)构造绝对值差异图像

对步骤(2)中滤波后的两幅图像对应像素点进行减法运算，并对减法运算的结果取绝对值，得到一幅绝对值差异图像；

(5)对步骤(3)中的对数差异图像进行Sc层快速离散曲波分解；

(6)根据曲波变换域频率的大小，对步骤(5)中的Curvelet变换系数进行分类；

(7)将Fine尺度层的Curvelet变换系数全部赋值为零；

(8)Detail尺度层去噪

8a)水平方向去噪：对Detail尺度层的任一方向子带，应用Treelet变换自适应的计算水平方向去噪阈值，并进行软阈值去噪，得到该方向子带的水平方向去噪结果；

8b)垂直方向去噪：对步骤8a)中选取的方向子带，应用Treelet变换自适应的计算垂直方向去噪阈值，并进行软阈值去噪，得到该方向子带的垂直方向去噪结果；

8c)平均加权求和：将水平方向去噪结果和垂直方向去噪结果进行平均加权求和，得到该方向子带的最终去噪结果；

8d)重复步骤8a)、8b)和8c)，直至处理完Detail尺度层的所有方向子带；

(9)曲波逆变换

对Coarse尺度层的Curvelet变换系数、经步骤(7)处理后的Fine尺度层的Curvelet变换系数和经步骤(8)处理后的Detail尺度层的Curvelet变换系数进行逆Curvelet变换，得到重构图像；

(10)计算变化比例阈值

10a)低频差异图像：将Fine尺度层和Detail尺度层的Curvelet变换系数全部置零，并与Coarse尺度层的Curvelet变换系数进行逆Curvelet变换，得到仅包含低频信息的低频差异图像；

10b)低频变化图像：对低频差异图像采用最大类间方差(Otsu)法分类，得到低频变化图像；

10c)最终边缘图：分别采用Canny、Sobel、Prewitt边缘算子提取低频变化图像的边缘，得到三幅边缘图，通过逻辑或运算将此三幅边缘图合成一幅最终边缘图；

10d)计算变化比例阈值：统计最终边缘图中的边缘像素点个数和低频变化图像中变化像素点的个数，计算变化比例阈值；

(11)分类

11a)通过变化比例阈值对步骤(4)获得的绝对值差异图像进行分类，得到一幅分类图像；

11b)通过变化比例阈值对步骤(9)获得的重构图像进行分类，得到另一幅分类图像；

(12)获得变化检测结果图

对步骤11a)获得的分类图像和步骤11b)获得的分类图像进行融合，得到变化检测结果图。

2.根据权利要求1所述的基于Treelet曲波域去噪的遥感图像变化检测方法，其特征在于：步骤3a)所述的对数运算公式为：

I₃(m，n)＝log(I₁(m，n)+1)

I₄(m，n)＝log(I₂(m，n)+1)

其中，I₃和I₄为两幅对数图像(大小为M×N)，I₁和I₂为滤波后的两幅图像，m和n为图像的行序号和列序号，m＝0，1，2，…，M-1，n＝0，1，2，…，N-1。

3.根据权利要求1所述的基于Treelet曲波域去噪的遥感图像变化检测方法，其特征在于：步骤(5)所述的快速离散曲波分解的分解层数为5。

4.根据权利要求1所述的基于Treelet曲波域去噪的遥感图像变化检测方法，其特征在于：步骤(6)所述的曲波变换系数分类的步骤如下：将第Sc分解层的低频系数归类为Coarse尺度层，第Sc-1、Sc-2、…、2分解层的高频系数归类为Detail尺度层，第1分解层的高频系数归类为Fine尺度层，其中，Sc为快速离散曲波分解的分解层数。

5.根据权利要求1所述的基于Treelet曲波域去噪的遥感图像变化检测方法，其特征在于：步骤8a)所述的计算水平方向去噪阈值的步骤如下：

计算第l分解层第i₁个尺度向量的归一化能量

$\mathrm{\ϵ}\left({\mathrm{\Φ}}_{i,i_1}\right)=\frac{{\mathrm{\Σ}}_{j_1=1}^K{|{\mathrm{\Φ}}_{l,i_1}\·x_{j_1}^{\left(0\right)}|}^2}{{\mathrm{\Σ}}_{j_1=1}^K{\left|\right|x_{j_1}^{\left(0\right)}\left|\right|}^2}$

其中，为第l分解层第i₁个尺度向量的归一化能量，为Treelet变换第l分解层基矩阵的第i₁个尺度向量，K为Detail尺度层任一方向子带的行向量个数，为该方向子带的第j₁个行向量，i₁＝1，2，…，K，l＝1，2，…，K-1；

计算第l分解层的总能量

$E_H^{\left(l\right)}={\mathrm{\Σ}}_{i_1=1}^{K-l}\mathrm{\ϵ}\left({\mathrm{\Φ}}_{l,i_1}\right)$

其中，

为第l分解层的总能量，K为该方向子带的行向量个数，为第l分解层第i₁个尺度向量的归一化能量；

选择最优分解层lb：

$\mathrm{lb}=\arg \underset{l=\mathrm{1,2},...,K-3}{\max }\frac{E_H^{(l+1)}-E_H^{(l+2)}}{E_H^{\left(l\right)}-E_H^{(l+1)}}$

其中，lb为最优分解层，l为分解层数，K为该方向子带的行向量个数，

为第l分解层的总能量，

为第l+1分解层的总能量，

为第l+2分解层的总能量，max为取最大值函数；

构造水平方向最优投影基矩阵PB：

PB＝[Φ_lb，1，…，Φ_lb，f，…，Φ_lb，(K-lb)]

其中，PB为水平方向最优投影基矩阵，lb为最优分解层，Φ_lb，1为第lb分解层基矩阵B_lb的第1个尺度向量，Φ_lb，f为B_lb的第f个尺度向量，Φ_lb，(K-lb)为B_lb的第K-lb个尺度向量；

水平方向主要成分的计算公式为：

$C_1=C_H^{\left(0\right)}\·\mathrm{PB}$

其中，C₁为水平方向主要成分，

为Detail尺度层任一方向子带的水平方向初始化结果，PB为水平方向最优投影基矩阵；

水平方向去噪阈值的计算公式为：

T₁＝mean(C₁)

其中，T₁为水平方向去噪阈值，C₁为水平方向主要成分，mean为均值运算。

6.根据权利要求1所述的基于Treelet曲波域去噪的遥感图像变化检测方法，其特征在于：步骤8c)所述的平均加权求和的计算公式为：

$C_{\mathrm{dn}}=\frac{1}{2}(C_{\mathrm{re}1}+C_{\mathrm{re}2})$

其中，C_dn为Detail尺度层任一方向子带的最终去噪结果，C_re1为该方向子带的水平方向去噪结果，C_re2为该方向子带垂直方向的去噪结果。

7.根据权利要求1所述的基于Treelet曲波域去噪的遥感图像变化检测方法，其特征在于：步骤(10)所述的变化比例阈值的计算公式如下：

P＝P₁+P₂

其中，P为变化比例阈值，P₁为最终边缘图中的边缘像素点总个数，P₂为低频变化检测图像中的变化像素点总个数。

8.根据权利要求1所述的基于Treelet曲波域去噪的遥感图像变化检测方法，其特征在于：步骤11a)所述的分类的步骤如下：将绝对值差异图像中所有像素点的灰度值按照由大到小的顺序排列，将前P个最大灰度值对应的像素点视为变化类，赋值为1；其他像素点视为非变化类，赋值为0，得到分类图像，其中，P为变化比例阈值。

9.根据权利要求1所述的基于Treelet曲波域去噪的遥感图像变化检测方法，其特征在于：步骤(12)所述的融合方式如下：

M_P(m，n)＝X₁(m，n)∩X₂(m，n)

其中，M_P(m，n)为变化检测结果图，X₁和X₂为两幅分类图像(大小为M×N)，m和n为图像的行序号和列序号，m＝0，1，2，…，M-1，n＝0，1，2，…，N-1，∩表示逻辑与运算。

Images