CN102289807A - 基于Treelet变换和特征融合的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于Treelet变换和特征融合的遥感图像变化检测方法，其实现步骤为：(1)读入数据；(2)中值滤波；(3)构造差异图像；(4)分类；(5)判断差异图像的标准差是否小于先验阈值；(6)自适应空间信息填充；(7)Treelet模糊融合；(8)构造模糊差异图像；(9)K-means分类；(10)数学形态学后处理；(11)特征与运算。本发明既可以较好的保持变化区域的边缘信息，又可以较好的兼顾变化检测结果中的漏检信息和虚警信息，具有较好的实时性和较高的检测精度，可应用于环境变化中的湖泊水位动态监测、农作物生长状态的动态监测、军事侦察等领域。

Description

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

技术领域

本发明属于图像处理技术领域，更进一步涉及一种基于Treelet变换和特征融合的遥感图像变化检测方法。该方法可应用于环境变化中的湖泊水位动态监测、农作物生长状态的动态监测、军事侦察等领域，能够快速的检测出两时相遥感图像的变化信息。

背景技术

变化检测是通过分析同一地区不同时刻的多幅遥感图像，检测出该地区地物随时间发生变化的信息。随着遥感技术和信息技术的发展，多时相遥感图像变化检测已经成为当前遥感图像分析研究的一个重要方向。

在多时相遥感图像变化检测方法的研究中，常见的一种检测方法是先比较后分类法，即首先构造一幅差异图像，然后利用阈值或分类法确定变化类和非变化类。先比较后分类法的优点在于简单易行，没有先分类后比较法所存在的分类误差累计问题，但该方法存在的明显不足是，对图像的预处理如几何校正、辐射校正、图像滤波等要求较高，并且对阈值选择方法和分类方法要求更加严格，影响了遥感图像变化检测的精度。

西安电子科技大学在其专利申请“基于Treelets的遥感图像变化检测方法”(专利申请号：201110001584.0，公开号：CN102063720A)中提出了一种Treelets滤波和K-means聚类相结合的遥感图像变化检测方法。该方法虽然能够减小辐射校正和光照不均对检测结果的影响，但仍存在的不足是，Treelets交叉滤波会导致检测结果中存在较多漏检信息，不能较好的保持变化区域边缘信息，降低了遥感图像的变化检测精度。此外，由于该方法使用Treelet变换对两幅图像中的每个像素点都进行滤波，使得该方法的时间复杂度较高，不具有实时性。

为了更加准确、全面的获取多时相遥感图像的变化信息，许多学者将图像融合技术应用到多时相遥感图像变化检测中，使多时相遥感图像间的互补信息得到更好的结合，达到更好的检测结果和检测精度。

2010年Celik和Ma在文献“Unsupervised Change Detection for Satellite Images UsingDual-Tree Complex Wavelet Transform”(IEEE Transactions on Geoscience and RemoteSensing，2010，48(3)：1199-1210.)中提出了一种基于双树复小波变换和尺度融合的变化检测方法。该方法首先对传统的双线性插值方法进行了改进，提出了一种加权双线性插值方法，并对两时相遥感图像分别进行加权双线性插值和双树复小波分解，然后通过尺度内融合得到同一分解层的变化检测结果，最后通过加权双线性插值和不同尺度间融合得到最终的变化检测结果。该方法虽然能够较好的保持变化区域的边缘信息，存在较少漏检信息，但仍存在的不足是，该方法通过逻辑与运算进行融合，为了获得较为全面的变化信息，导致检测结果中存在较多虚警信息，降低了遥感图像变化检测精度，难以较好兼顾变化检测结果中的漏检信息和虚警信息。

发明内容

本发明针对上述现有技术存在的不足，提出了一种基于Treelet变换和特征融合的遥感图像变化检测方法。本发明既可以较好的保持变化区域的边缘信息，又可以较好的兼顾变化检测结果中的漏检信息和虚警信息，具有较好的实时性和较高的检测精度。

本发明实现上述目的的思路是：在对读入的数据进行中值滤波后，先对构造的差异图像进行统计均值阈值分类和自适应空间信息填充，再对滤波后的图像利用Treelet变换、K-means聚类和数学形态学后处理生成模糊结果图，最后通过特征与运算进行融合。

本发明的步骤包括如下：

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

(2)中值滤波

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

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

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

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

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

(3)构造差异图像

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

3b)计算该差异图像的标准差。

(4)分类

4a)按照下式计算精细分类阈值；

$T=\arg \underset{G}{\min }\left(\right|P_G-\frac{\underset{G=0}{\overset{255}{\mathrm{\Σ}}}{P}_G}{256}\left|\right)$

其中，T为精细分类阈值，min为取最小值函数，P_G为差异图像中灰度值为G的像素点总个数，差异图像的大小为M×N，G∈{0，1，2…，255}，P_G∈{0，1，2…，M×N}；

4b)利用精细分类阈值按照下式对差异图像进行分类，得到精细分类图；

$D_A(m,n)=\left\{\begin{array}{c}1,A(m,n)\&GreaterEqual;T\\ 0,A(m,n)<T\end{array}\right.$

其中，D_A为精细分类图，A为差异图像(大小为M×N)，T为精细分类阈值，m和n为图像的行序号和列序号，m＝1，2，…，M，n＝1，2，…，N。

(5)判断差异图像的标准差是否小于先验阈值

若差异图像的标准差小于先验阈值T₁，则进行步骤(6)，否则步骤(4)中的精细分类图即为精细结果图，进行步骤(7)，其中，T₁∈{10，11，…，19，20}。

(6)自适应空间信息填充

6a)确定图像块：在精细分类图中选取一个待填充像素点，以该像素点为中心，以固定长度N₂个像素为半径，确定一个图像块，其中，N₂∈{1，2，3}；

6b)计算填充阈值：根据步骤6a)中的半径计算填充阈值，并统计图像块中灰度值为1的像素点总个数；

T_h＝N₂×(2×N₂+1)

其中，T_h为填充阈值，N₂为图像块半径；

6c)填充：当像素点总个数大于或等于填充阈值时，将待填充像素点的灰度值赋值为1，否则将待填充像素点的灰度值赋值为0；

6d)重复步骤6a)和步骤6c)，直至处理完精细分类图中的全部像素点，得到精细结果图。

(7)Treelet模糊融合

7a)构造低频序列：分别对步骤(2)中滤波后的两幅遥感图像进行Q层平稳小波分解，构造两个低频序列；

7b)构造Treelet基矩阵：对步骤7a)中的一个低频序列进行Treelet变换，得到Treelet基矩阵；

7c)获得模糊图像：分别将两个低频序列向Treelet基矩阵投影，得到两幅模糊图像。

(8)构造模糊差异图像

对步骤7c)中的两幅模糊图像对应像素点进行减法运算，并对减法运算的结果取绝对值，得到一幅模糊差异图像。

(9)采用K-means方法对模糊差异图像分类，得到模糊分类图。

(10)数学形态学后处理

10a)对模糊分类图采用一个半径为4个像素的圆形结构元素进行数学形态学闭合运算，得到闭合运算结果；

10b)对闭合运算结果采用步骤10a)中的结构元素进行数学形态学开启运算，得到模糊结果图。

(11)特征与运算

11a)确定空间位置：对精细结果图和模糊结果图进行逻辑与运算融合，确定待查找对象的空间位置；

11b)对象选择：在精细结果图中，查找步骤11a)中的空间位置所对应的W连通对象，即为最终的变化检测结果。

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

第一，本发明利用Treelet模糊融合的方法查找变化区域的模糊轮廓，克服了现有技术存在较多虚警信息的缺点，减少了本发明检测结果中的虚警信息。

第二，本发明通过统计均值阈值分类法得到精细分类图，克服了现有技术存在较多漏检信息的缺点，减少了本发明检测结果中的漏检信息，并使得本发明的变化区域边缘信息得到了较为完整的保持。

第三，本发明通过特征与运算进行图像融合，克服了现有技术难以同时较好的兼顾变化检测结果中的虚警信息和漏检信息的缺点，使得本发明的变化检测精度得到了提高。

第四，本发明很少采用算法复杂度较大、逐点处理像素的方法，克服了现有技术时间复杂度较大的缺点，使得本发明的处理速度得到了提高，具有较好的实时性。

附图说明

图1为本发明的流程图；

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

具体实施方式

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

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

步骤2，中值滤波。

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

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

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

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

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

步骤3，构造差异图像。

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

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

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

3b)计算步骤3a)中的差异图像的标准差。

其中，σ为差异图像的标准差，

为求数学期望，A为差异图像。

步骤4，分类。

4a)按照下式计算精细分类阈值。

$T=\arg \underset{G}{\min }\left(\right|P_G-\frac{\underset{G=0}{\overset{255}{\mathrm{\&Sigma;}}}{P}_G}{256}\left|\right)$

其中，T为精细分类阈值，min为取最小值函数，P_G为差异图像中灰度值为G的像素点总个数，差异图像的大小为M×N，G∈{0，1，2…，255}，P_G∈{0，1，2…，M×N}。

4b)利用精细分类阈值按照下式对差异图像进行分类，得到精细分类图。

$D_A(m,n)=\left\{\begin{array}{c}1,A(m,n)\&GreaterEqual;T\\ 0,A(m,n)<T\end{array}\right.$

其中，D_A为精细分类图，A为差异图像(大小为M×N)，T为精细分类阈值，m和n为图像的行序号和列序号，m＝1，2，…，M，n＝1，2，…，N。

步骤5，判断差异图像的标准差是否小于先验阈值。

若差异图像的标准差小于先验阈值T₁，则进行步骤6，否则将步骤4中的精细分类图作为最终的精细结果图，进行步骤7，其中，T₁∈{10，11，…，19，20}，本发明实施例中选取先验阈值T₁＝15。

步骤6，自适应空间信息填充。

6a)确定图像块：在精细分类图中选取一个待填充像素点，以该像素点为中心，以固定长度N₂个像素为半径，确定一个图像块，其中，N₂∈{1，2，3}，本发明实施例中选取的固定长度为2个像素。

6b)计算填充阈值：根据步骤6a)中的半径计算填充阈值，并统计图像块中灰度值为1的像素点总个数。

T_h＝N₂×(2×N₂+1)

其中，T_h为填充阈值，N₂为图像块半径。

6c)填充：当像素点总个数大于或等于填充阈值时，将待填充像素点的灰度值赋值为1，否则将待填充像素点的灰度值赋值为0。

6d)重复步骤6a)和步骤6c)，直至处理完精细分类图中的全部像素点，得到精细结果图D_re。

步骤7，Treelet模糊融合。

7a)构造低频序列：分别对步骤(2)中滤波后的两幅遥感图像进行Q层平稳小波分解，构造两个低频序列，其中，Q∈{3，4，5}，本发明实施例中进行4层平稳小波分解。

a1)对步骤(2)中滤波后的一幅遥感图像进行4层平稳小波分解，得到4个低频子带L(level)和12个高频子带(水平方向子带H(level)、垂直方向子带V(level)、对角方向子带Z(level)，level＝1，2，3，4)。

a2)在第level(level∈1，2，3，4)分解层，将H(level)、V(level)和Z(level)的系数全部置为零，并与L(level)进行平稳小波逆变换，得到一幅低频图像。

a3)将该低频图像中的像素点按照从左到右、从上到下的顺序进行排列，组成一个列向量(大小为(M×N)×1)。

a4)重复步骤a2)和步骤a3)，直至处理完平稳小波分解的所有分解层，得到四个列向量。

a5)将第level分解层对应的列向量置于低频序列第level个列向量的位置，构造一个低频序列X₁(大小为(M×N)×4)。

a6)按照步骤a1)至步骤a5)，对步骤(2)中的另一幅遥感图像进行处理，得到一个低频序列X₂。

7b)构造Treelet基矩阵：对步骤7a)中的一个低频序列进行Treelet变换，得到Treelet基矩阵，其中，本发明实施例中选取低频序列X₂进行Treelet变换。

b1)Treelet变换第l＝0层：

初始化X₂为

和变量集为δ＝{1，2，3，4}，正交基矩阵为B₀＝[Φ_0，1，Φ_0，2，Φ_0，3，Φ_0，4]，其中，B₀是一个4×4的单位矩阵。

计算的协方差矩阵∑⁽⁰⁾。

其中，

为求数学期望，h＝{1，2，3，4}和v＝{1，2，3，4}分别为协方差矩阵∑⁽⁰⁾的位置索引。

计算的相似度矩阵Θ⁽⁰⁾。相似度矩阵Θ⁽⁰⁾中的每个元素均为相关系数：

${\mathrm{\&Theta;}}_{(h,v)}^{\left(0\right)}=\left|\frac{{\mathrm{\&Sigma;}}_{\mathrm{hv}}^{\left(0\right)}}{\sqrt{{\mathrm{\&Sigma;}}_{\mathrm{hh}}^{\left(0\right)}{\mathrm{\&Sigma;}}_{\mathrm{vv}}^{\left(0\right)}}}\right|$

b2)Treelet变换第l＝1，2，3层：

根据相似度矩阵Θ^(l-1)找到两个最相似的和变量，即：

$(\mathrm{\&alpha;},\mathrm{\&beta;})=\arg \underset{h,v\&Element;\mathrm{\&delta;}}{\max }{\mathrm{\&Theta;}}_{(h,v)}^{(l-1)}$

其中，h＜v。

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

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

|θ_l|≤π/4

∑^(l)＝J^T∑^(l-1)J

${\mathrm{\&Sigma;}}_{\mathrm{\&alpha;\&beta;}}^{\left(l\right)}={\mathrm{\&Sigma;}}_{\mathrm{\&beta;\&alpha;}}^{\left(l\right)}=0$

利用Jacobi旋转矩阵J更新第l分解层的基矩阵B_l＝B_l-1J＝[Φ_l，1，Φ_l，2，Φ_l，3，Φ_l，4]和第l分解层的

其中，上标T表示转置。

经过Jacobi旋转后，使

和

满足关系式

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

和

并定义第l分解层的尺度向量Ω_l和细节向量Ψ_l分别为基矩阵B_l的第α列和第β列。

b3)提取第l＝3层基矩阵B₃的尺度向量，该尺度向量即为Treelet基矩阵PB。

PB＝[Φ_3，1]

7c)获得模糊图像：分别将两个低频序列向Treelet基矩阵投影，得到两幅模糊图像。

F₁＝X₁·PB

F₂＝X₂·PB

其中，F₁和F₂为两幅模糊图像，PB为Treelet基矩阵，X₁和X₂为两个低频序列。

步骤8，构造模糊差异图像。

对步骤7c)中的两幅模糊图像对应像素点进行减法运算，并对减法运算的结果取绝对值，得到一幅模糊差异图像。

F_d(m，n)＝|F₁(m，n)-F₂(m，n)|

其中，F_d为模糊差异图像，F₁和F₂为两幅模糊图像(大小为M×N)，m和n为图像的行序号和列序号，m＝1，2，…，M，n＝1，2，…，N。

步骤9，采用K-means方法对模糊差异图像进行分类，得到模糊分类图。

步骤10，数学形态学后处理。

10a)对模糊分类图采用一个半径为4个像素的圆形结构元素进行数学形态学闭合运算，得到闭合运算结果。

10b)对闭合运算结果采用步骤10a)中的结构元素进行数学形态学开启运算，得到模糊结果图F_re。

步骤11，特征与运算。

11a)确定空间位置：对精细结果图D_re和模糊结果图F_re进行逻辑与运算融合，确定待查找对象的空间位置。

Γ＝{(x，y)|(D_re(x，y)＝1)∩(F_re(x，y)＝1)}

其中，Γ为待查找对象空间位置的集合，D_re为精细结果图(大小为M×N)，F_re为模糊结果图，x和y为图像的行序号和列序号，x∈{1，2，…，M}，y∈{1，2，…，N}，∩表示逻辑与运算，即当且仅当同时满足D_re(x，y)＝1和F_re(x，y)＝1两个条件时，空间位置坐标(x，y)才属于集合Γ。

11b)对象选择：在精细结果图中，查找步骤11a)中的空间位置所对应的W连通对象，即为最终的变化检测结果，其中，W∈{4，8}，本发明实施例中采用8连通对象。

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

1.仿真条件

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

本发明实施例中采用了3×3的中值滤波和4层平稳小波分解，自适应空间信息填充步骤中选取5×5的图像块，并且先验阈值T₁＝15，计算Treelet基矩阵步骤中选取低频序列X₂，数学形态学后处理步骤中采用半径为4个像素的圆形结构元素，特征与运算步骤中选取的是8连通对象。

2.仿真内容

本发明仿真实验所用数据为两组真实遥感数据集。第一组真实遥感数据集是墨西哥郊外的两幅Landsat7 ETM+第4波段光谱图像，两幅图像的大小均为512×512像素，它们之间发生的变化是由火灾破坏了大面积的当地植被所致，包括25589个变化像素和236555个非变化像素。第二组真实遥感数据集是1994年8月和1994年9月意大利Elba岛西部地区的两幅Landsat-5 TM第4波段光谱图像，两幅图像的大小均为384×320像素，它们之间发生的变化是由森林火灾破坏大量植被引起的，包含2415个变化像素和120465个非变化像素。

本发明采用虚警数、漏检数、总错误数和仿真运行时间四个指标来评价变化检测方法的好坏。

3.仿真效果分析

本发明中，精细结果图具有漏检信息少、边缘信息保持较好的优点，模糊结果图具有虚警信息少、变化区域定位准确的优点。因此，采用特征与运算对精细结构图和模糊结果图进行融合，能够较好的兼顾变化检测结果中的漏检信息和虚警信息，并且较好的保持了变化区域边缘信息，具有较好的实时性和较高的检测准确度。

为了说明本发明的有效性和优越性，将背景技术中提到的两种遥感图像变化检测方法作为对比方法进行验证。其中，将专利“基于Treelets的遥感图像变化检测方法”中的方法简记为TC法，将文献“Unsupervised Change Detection for Satellite Images UsingDual-Tree Complex Wavelet Transform”中的方法简记为DTCWT法。

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

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

从表1中可以看出，在三种方法的变化检测结果中，本发明对第一组和第二组真实遥感数据集的整体评价是最好的。相对于TC法的8个评价指标，本发明有6个评价指标是最好的，即：第一组和第二组真实遥感数据集的总错误数，亦即变化检测精度，比TC法的总错误数分别少1937像素点和60个像素点；第一组和第二组真实遥感数据集的漏检数，比TC法的漏检数分别少2456像素点和192个像素点；第一组和第二组真实遥感数据集的运行时间，亦即实时性，比TC法的运行时间分别少46.6s和14.1s。相对于DTCWT法的8个评价指标，本发明有7个评价指标是最好的，即：第一组和第二组真实遥感数据集的总错误数，亦即变化检测精度，比DTCWT法的总错误数分别少1814像素点和168个像素点；第一组和第二组真实遥感数据集的虚警数，比DTCWT法的虚警数分别少2326像素点和43个像素点；第二组真实遥感数据集的漏检数，比DTCWT法的漏检数少125个像素点；第一组和第二组真实遥感数据集的运行时间，亦即实时性，比DTCWT法的运行时间分别少79s和37.2s。由此可以看出，本发明能够较好的兼顾变化检测结果中的漏检信息和虚警信息，并且具有较高的检测准确度和较好的实时性。从两组实验数据集的效果图中可以看出，与TC法和DTCWT法相比，本发明能够较好的保持变化区域的边缘信息。

Claims (5)

1.基于Treelet变换和特征融合的遥感图像变化检测方法，包括如下步骤：

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

(2)中值滤波

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

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

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

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

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

(3)构造差异图像

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

3b)计算该差异图像的标准差；

(4)分类

4a)按照下式计算精细分类阈值；

$T=\arg \underset{G}{\min }\left(\right|P_G-\frac{\underset{G=0}{\overset{255}{\mathrm{\&Sigma;}}}{P}_G}{256}\left|\right)$

其中，T为精细分类阈值，min为取最小值函数，P_G为差异图像中灰度值为G的像素点总个数，差异图像的大小为M×N，G∈{0，1，2…，255}，P_G∈{0，1，2…，M×N}；

4b)利用精细分类阈值按照下式对差异图像进行分类，得到精细分类图；

$D_A(m,n)=\left\{\begin{array}{c}1,A(m,n)\&GreaterEqual;T\\ 0,A(m,n)<T\end{array}\right.$

其中，D_A为精细分类图，A为差异图像(大小为M×N)，T为精细分类阈值，m和n为图像的行序号和列序号，m＝1，2，…，M，n＝1，2，…，N；

(5)判断差异图像的标准差是否小于先验阈值

若差异图像的标准差小于先验阈值T₁，则进行步骤(6)，否则步骤(4)中的精细分类图即为精细结果图，进行步骤(7)，其中，T₁∈{10，11，…，19，20}；

(6)自适应空间信息填充

6a)确定图像块：在精细分类图中选取一个待填充像素点，以该像素点为中心，以固定长度N₂个像素为半径，确定一个图像块，其中，N₂∈{1，2，3}；

6b)计算填充阈值：根据步骤6a)中的半径计算填充阈值，并统计图像块中灰度值为1的像素点总个数；

T_h＝N₂×(2×N₂+1)

其中，T_h为填充阈值，N₂为图像块半径；

6c)填充：当像素点总个数大于或等于填充阈值时，将待填充像素点的灰度值赋值为1，否则将待填充像素点的灰度值赋值为0；

6d)重复步骤6a)和步骤6c)，直至处理完精细分类图中的全部像素点，得到精细结果图；

(7)Treelet模糊融合

7a)构造低频序列：分别对步骤(2)中滤波后的两幅遥感图像进行Q层平稳小波分解，构造两个低频序列；

7b)构造Treelet基矩阵：对步骤7a)中的一个低频序列进行Treelet变换，得到Treelet基矩阵；

7c)获得模糊图像：分别将两个低频序列向Treelet基矩阵投影，得到两幅模糊图像；

(8)构造模糊差异图像

对步骤7c)中的两幅模糊图像对应像素点进行减法运算，并对减法运算的结果取绝对值，得到一幅模糊差异图像；

(9)采用K-means方法对模糊差异图像分类，得到模糊分类图；

(10)数学形态学后处理

10a)对模糊分类图采用一个半径为4个像素的圆形结构元素进行数学形态学闭合运算，得到闭合运算结果；

10b)对闭合运算结果采用步骤10a)中的结构元素进行数学形态学开启运算，得到模糊结果图；

(11)特征与运算

11a)确定空间位置：对精细结果图和模糊结果图进行逻辑与运算融合，确定待查找对象的空间位置；

11b)对象选择：在精细结果图中，查找步骤11a)中的空间位置所对应的W连通对象，即为最终的变化检测结果。

2.根据权利要求1所述的基于Treelet变换和特征融合的遥感图像变化检测方法，其特征在于：步骤7a)所述的平稳小波分解的层数Q∈{3，4，5}。

3.根据权利要求1所述的基于Treelet变换和特征融合的遥感图像变化检测方法，其特征在于：步骤7b)所述的低频序列为X₂。

4.根据权利要求1所述的基于Treelet变换和特征融合的遥感图像变化检测方法，其特征在于：步骤11a)所述的确定空间位置的公式为：

Γ＝{(x，y)|(D_re(x，y)＝1)∩(E_re(x，y)＝1)}

其中，Γ为待查找对象空间位置的集合，D_re为精细结果图(大小为M×N)，F_re为模糊结果图，x和y为图像的行序号和列序号，x∈{1，2，…，M}，y∈{1，2，…，N}，∩表示逻辑与运算。

5.根据权利要求1所述的基于Treelet变换和特征融合的遥感图像变化检测方法，其特征在于：步骤11b)所述对象的连通方式为W∈{4，8}。

Images