CN103077525A - 基于Treelet图像融合的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于Treelet图像融合的遥感图像变化检测方法，主要解决现有技术中差值法抑制噪声能力差和对数比值法对边缘等细节信息保持不好的问题。其实现过程是：用不同方法提取变化前后图像的4幅差异图；然后用Treelet对4幅差异图构成的矩阵像进行降维，在降维时每一幅差异图数据作为待降维数据的一维；最后用FCM对降维后的数据进行聚类分割，得到变化监测的结果。本发明克服了差值图对噪声敏感和比值图边缘信息丢失的缺点，降低了错误率，更好的保存了细节信息，可用于森林资源的动态监测、土地覆盖和利用的变化监测、自然灾害评估。

Description

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

技术领域

本发明属于图像处理技术领域，主要是在遥感图像变化检测领域的应用，即根据不同时相的遥感图像来检测地表产生的变化，主要可应用于地物覆盖与利用、自然灾害监控，城区规划，地图更新以军事领域等。

背景技术

遥感图像的变化检测是指通过分析在不同时间来自同一地区的两幅或多幅遥感图像，检测出该地区的地物随时间发生的变化信息。遥感图像的变化检测已经广泛地应用于如森林资源的动态监测、土地覆盖和利用的变化监测、城市规划布局、环境监测分析、自然灾害评估、以及在军事中的人造目标监测和地面武装部署分析等许多领域。

在变化检测研究中，学者们将现有的多种变化检测方法归为不同的类别，其中最常见的是基于差异图像处理的变化检测方法。该方法通常包括3个关键的步骤：1)图像预处理；2)构造差异图像；3)分析差异图像提取变化检测结果。图像预处理包括图像几何校正，几何配准，辐射校正，图像去噪和图像增强等。其中对图像进行去噪处理则能很大程度上降低噪声，提高检测精度，但同时又会影响图像的细节信息，降低视觉效果。构造差异图像是基于差异图像分析方法的重要一步，尤其是构造视觉效果较好的差异图像难度比较大。其中，图像差值法和比值法是最常见的差异图提取方法。图像差值法对遥感图像的变化区域非常敏感，有利于保留变化信息的几何边缘，从而有利于精确地提取变化区域，但由于遥感图像噪声的影响，差值法并不适用于分析受噪声影响较大的遥感图像。而图像对数比能有效地抑制遥感图像的噪声，但却会对变化区域的边缘造成一定的模糊，特别是对于存在大量的微弱变化和小区域变化信息的图像，该方法不利于变化区域轮廓和细节信息的保持。

发明内容

本发明的目：本发明的目在于针对上述已有技术的不足，提出了一种基于Treelet图像融合的遥感图像变化检测方法，以减小信息丢失和累积误差，提高SAR图像变化检测的检测精度。

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

（1）输入变化前的图像P₁和变化后的图像P₂，分别进行中值滤波和均值滤波，得到中值滤波后的变化前图像P_1m，均值滤波后的变化前图像P_1a，中值滤波后的变化后图像P_2m，均值滤波后的变化后图像P_2a；

（2）用差值法得到中值滤波后的变化前图像P_1m和中值滤波后的变化后图像P_2m的差异图H₁，用对数比值法得到中值滤波后的变化前图像P_1m和中值滤波后的变化后图像P_2m的差异图H₂，用差值法得到均值滤波后的变化前图像P_1a和均值滤波后的变化后图像P_2a的差异图H₃，用对数比值法得到均值滤波后的变化前图像P_1a和均值滤波后的变化后图像P_2a的差异图H₄；

（3）分别将上述差异图H₁，H₂，H₃，H₄转化成一维的列向量，X₁，X₂，X₃，X₄，并将该4个列向量组成一个矩阵X=[X₁X₂X₃X₄]；

（4）对矩阵X用Treelet法进行降维，得到列向量P，并将该向量P转化成最终差异图D的灰度矩阵

（5）对所得的最终差异图D的灰度矩阵进行聚类分割，得到变化检测结果并输出。

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

1.本发明融合前的差值图主要采用差值法和对数比值法提取，通过差值法能很好的保持图像的边缘特征，通过比值法能有效地抑制了图像的噪声，并且这两种方法操作简单，易于理解，运算速度快；

2.本发明通过Treelet方法实现对4幅差异图的融合，综合了中值滤波平滑效果好和均值滤波细节保持好的优势，降低了错误率，更好的保存了细节信息。

附图说明

图1是本发明的流程图；

图2是本发明仿真使用的Feltwell遥感图像数据集；

图3是本发明仿真使用的意大利撒丁岛遥感图像数据集；

图4为对图2的变化检测仿真实验结果图；

图5为对图3的变化检测仿真实验结果图。

具体实施方式

参照图1，本发明基于Treelet图像融合的遥感图像变化检测方法，包括如下步骤：

步骤1:对输入图像进行去噪。

输入变化前的图像P₁，对变化前的图像P₁进行中值滤波得到中值滤波后的变化前图像P_1m，对变化前的图像P₁进行均值滤波得到均值滤波后的变化前图像P_1a；

输入变化后的图像P₂，对变化后的图像P₂进行中值滤波得到中值滤波后的变化后图像P_2m，对变化后的图像P₂进行均值滤波得到均滤波后的变化后图像P_2a。

步骤2:提取不同的差异图。

用差值法得到中值滤波后的变化前图像P_1m和中值滤波后的变化后图像P_2m的差异图H₁；

用对数比值法得到中值滤波后的变化前图像P_1m和中值滤波后的变化后图像P_2m的差异图H₂；

用差值法得到均值滤波后的变化前图像P_1a和均值滤波后的变化后图像P_2a的差异图H₃；

用对数比值法得到均值滤波后的变化前图像P_1a和均值滤波后的变化后图像P_2a的差异图H₄。

步骤3:产生待降维矩阵X。

将所述差异图H₁，H₂，H₃，H₄的灰度矩阵分别转化为化成列向量X₁，X₂，X₃，X₄，并将该4个列向量组成一个矩阵X：

其中，x_st为矩阵X第s行t列的值，t=1,2,3,4，s=1,2,3,…N，N为差异图的像素点数；

步骤4:对矩阵X用Treelet法进行降维，得到列向量P，并将该向量P转化成最终差异图D的灰度矩阵：

4.1）计算X的协方差矩阵C：

$C=\left[\begin{array}{cccc}{C}_{\mathrm{1,1}}& C_{\mathrm{1,2}}& ...& C_{\mathrm{1,4}}\\ C_{\mathrm{2,1}}& C_{\mathrm{2,2}}& ...& C_{\mathrm{2,4}}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ C_{\mathrm{4,1}}& C_{\mathrm{4,2}}& ...& C_{\mathrm{4,4}}\end{array}\right],$

其中协方差矩阵C中的每一个元素C_m，n按照如下公式计算：

$C_{m,n}=\frac{\underset{k=1}{\overset{N}{\mathrm{\Σ}}}(x_{k,m}-{\stackrel{\‾}{x}}_m)(x_{k,n}-{\stackrel{\‾}{x}}_n)}{\sqrt{\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{(x_{k,m}-{\stackrel{\‾}{x}}_m)}^2\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{(x_{k,n}-{\stackrel{\‾}{x}}_n)}^2}},$

式中，

和

分别表示矩阵X第m列和第n列元素的均值，x_k，m,为矩阵X在第k行第m列的值，x_k，n为矩阵X在第k行n列的值，m,n=1,2,3,4，k＝1,2,3,…N，N为差异图像素点数；

4.2）计算相关系数矩阵ρ，以此来度量两个变量之间的相似性

$\mathrm{\ρ}=\left[\begin{array}{cccc}{\mathrm{\ρ}}_{\mathrm{1,1}}& {\mathrm{\ρ}}_{\mathrm{1,2}}& ...& {\mathrm{\ρ}}_{\mathrm{1,4}}\\ {\mathrm{\ρ}}_{\mathrm{2,1}}& {\mathrm{\ρ}}_{\mathrm{2,2}}& ...& {\mathrm{\ρ}}_{\mathrm{2,4}}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\ρ}}_{\mathrm{4,1}}& {\mathrm{\ρ}}_{\mathrm{4,2}}& ...& {\mathrm{\ρ}}_{\mathrm{4,4}}\end{array}\right],$

其中

${\mathrm{\ρ}}_{a,b}=\frac{C_{a,b}}{\sqrt{C_{a,b}{C}_{a,b}}}$

Ρ_a，b表示矩阵ρ的第a行、第b列个元素，C_a，b为协方差矩阵C的第a行、第b列个元素,a，b=1,2,3,4；

4.3）初始化总的聚类层数为l=0，初始化正交Dirac基B⁽⁰⁾为4x4维的单位矩阵，C⁽⁰⁾=C；

4.4)设l＝l+1，寻找相关系数矩阵ρ中最大的值，将该值的座标序号分别标记为α和β，α>β；

4.5)利用得到的位置序号α和β，计算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{\α}}^{(l-1)}-C_{\mathrm{\β},\mathrm{\β}}^{(l-1)}\±\sqrt{{[C_{\mathrm{\α},\mathrm{\α}}^{(l-1)}-C_{\mathrm{\β},\mathrm{\β}}^{(l-1)}]}^2+4{\left[C_{\mathrm{\α},\mathrm{\β}}^{(l-1)}\right]}^2}}{2C_{\mathrm{\α},\mathrm{\β}}^{(l-1)}}\right]$ 且 $\left|{\mathrm{\θ}}_l\right|\≤\frac{\mathrm{\π}}{4},$

为第l-1层协方差矩阵C在第α行第α列的值，

为第l-1层协方差矩阵C在第α行第β列的值，

为第l-1层协方差矩阵C在第β行第β列的值；

4.6)利用旋转矩阵J得到去相关后的正交基B^(l)=B^(l-1)J，更新协方差矩阵C^(l)=J^TC^(l-1)J，其中B^(l)为更新后的正交基，C^(l)为更新后的协方差矩阵，根据协方差矩阵C^(l)更新相关系数矩阵ρ^(l)；

4.7）判断l的值，如果l＝3，则执行步骤4.8），否则返回步骤4.4）；

4.8）将矩阵X沿正交基矩阵B⁽³⁾的转置方向进行投影，得到投影后的矩阵P_n：

P_n=X×B^(3)T；

4.9）将矩阵P_n的第一列P转化为与输入图像同等大小的矩阵，该矩阵即为融合后最终差异图D的灰度矩阵。

步骤5：利用模糊C均值聚类算法FCM对差异图像D的灰度矩阵进行聚类分割，得到变化检测结果。

5.1)将最终差异图D的灰度矩阵转化为一个列向量T=[t_i]，t_j为该向量第i个像素的灰度值，i=1,2,3,…N，N为总像素点数，将第i个像素属于第k类的隶属度记为u_ki，u_ki∈[0,1]且

k＝1,2，随机初始化隶属度矩阵U^b=[u_ki]的每一个元素，设定循环计数器b=0；

5.2)按如下公式计算第k类的聚类中心c_k：

$c_k=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{u}_{\mathrm{ki}}^2{t}_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{u}_{\mathrm{ki}}^2};$

5.3)按如下公式计算隶属度矩阵U^（b+1）:

U^（b+1)={u_ki}，

其中，k＝1,2, $u_{1i}=\frac{1}{1+{\left(\frac{t_i-c_1}{t_i-c_2}\right)}^2},$ $u_{2i}=\frac{1}{{\left(\frac{t_i-c_2}{t_i-c_1}\right)}^2+1};$

5.4)如果满足max{U^(b)-U^(b+1)}<0.001则转至步骤5.5），否则设置b=b+1,转至步骤5.2)；

5.5)通过比较c₁与c₂的大小确定变化类与非变化类，如果c₁≥c₂，则以c₁为聚类中心的那一类为变化类，以c₂为聚类中心的那一类为非变化类；如果c₁<c₂，以c₂为聚类中心的那一类为变化类，以c₁为聚类中心的一类为非变化类。

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

1.实验条件

实验环境为：windows XP，SPI，CPU Pentium（R）4，基本频率为2.4GHZ，软件平台为MatlabR2010a。

第一组数据集为英格兰Feltwell村庄遥感图像数据集,如图2所示，其中2（a）图像是ATM（Airborne Thematic Mapper）第3波段的位于英国Feltwell村庄农的田区图像，，图2(b)是通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工的嵌入一些变化区域得到的。图像大小均为470×335像素，灰度级为256。检测的标准结果图采用如图2（c）所示的对Feltwell SAR图像数据集变化检测的结果图，包括153214个非变化像素和4236个变化像素。

第二个数据集为意大利撒丁岛遥感图像数据集，如图3所示，其中图3(a)为意大利撒丁岛变化前的图像，图3(b)为意大利撒丁岛变化后的图像。该组真实遥感图像由1995年9月和1996年7月的Landsat-5卫星TM（Thematic Mapper）第5波段的两幅意大利撒丁岛Mulargia湖泊区域的光谱组成。图像大小均为300×412像素，灰度级为256。检测的标准结果图采用如图3（c）所示的对撒丁岛SAR图像数据集变化检测的结果图，包括115974个非变化像素和7626个变化像素。

2.实验内容和实验结果

为了验证本发明所提出方法的可靠性，将本发明方法与中值融合法和均值融合法进行比较。

中值融合法是指图像融合部分用Treelet法进融合差异图H₁和H₂。

均值融合法是指图像融合部分用Treelet法进融合差异图H₃和H₄。

实验一：用本发明方法与中值融合法和均值融合法分别对第一组数据进行变化检测，其实验结果如图4所示，其中4(a)为中值融合法的检测结果图，4(b)为均值融合法的检测结果图，4(c)为本发明方法的检测结果图。

实验二：用本发明方法与中值融合法和均值融合法分别对第二组数据进行变化检测，其实验结果如图5所示，其中5(a)为中值融合法的检测结果图，5(b)为均值融合法的检测结果图，5(c)为本发明方法的检测结果图。

中值融合法，均值融合法和本发明方法对图2和图3进行变化检测的数据结果，如表1所示：

表1试验数据结果

表1中列出了三种评价指标：分别为漏检数，误检数和总错误数，其中，漏检数为没有检测出来的实际发生了变化的像素，误检数为实际没有发生变化但被当作变化的像素检测出来，总错误数=漏检数+误检数。

从表1可以看出，本发明与所述对比方法相比，可以获得较少的漏检数和误检数，以及最少的总错误数，提高了变化检测的检测精度。

Claims (3)

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

（1）输入变化前的图像P₁和变化后的图像P₂，分别进行中值滤波和均值滤波，得到中值滤波后的变化前图像P_1m，均值滤波后的变化前图像P_1a，中值滤波后的变化后图像P_2m，均值滤波后的变化后图像P_2a；

（2）用差值法得到中值滤波后的变化前图像P_1m和中值滤波后的变化后图像P_2m的差异图H₁，用对数比值法得到中值滤波后的变化前图像P_1m和中值滤波后的变化后图像P_2m的差异图H₂，用差值法得到均值滤波后的变化前图像P_1a和均值滤波后的变化后图像P_2a的差异图H₃，用对数比值法得到均值滤波后的变化前图像P_1a和均值滤波后的变化后图像P_2a的差异图H₄；

（3）分别将上述差异图H₁，H₂，H₃，H₄转化成一维的列向量，X₁，X₂，X₃，X₄，并将该4个列向量组成一个矩阵X=[X₁X₂X₃X₄]；

（4）对矩阵X用treelet法进行降维，得到列向量P，并将该向量P转化成最终差异图D的灰度矩阵；

（5）对所得的最终差异图D的灰度矩阵进行聚类分割，得到变化检测结果并输出。

2.根据权利要求书1中所述的方法，其中步骤(4)所述的对矩阵X用treelet法进行降维处理，降为一维并将该向量转化成最终差异图D的灰度矩阵，按如下步骤进行：

4.1）计算X的协方差矩阵C：

$C=\left[\begin{array}{cccc}{C}_{\mathrm{1,1}}& C_{\mathrm{1,2}}& ...& C_{\mathrm{1,4}}\\ C_{\mathrm{2,1}}& C_{\mathrm{2,2}}& ...& C_{\mathrm{2,4}}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ C_{\mathrm{4,1}}& C_{\mathrm{4,2}}& ...& C_{\mathrm{4,4}}\end{array}\right]$

其中协方差矩阵C中的每一个元素C_m，n按照如下公式计算：

$C_{m,n}=\frac{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_{k,m}-{\stackrel{\&OverBar;}{x}}_m)(x_{k,n}-{\stackrel{\&OverBar;}{x}}_n)}{\sqrt{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{(x_{k,m}-{\stackrel{\&OverBar;}{x}}_m)}^2\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{(x_{k,n}-{\stackrel{\&OverBar;}{x}}_n)}^2}}$

和

分别表示矩阵X第m列和第n列元素的均值,x_k，m,为矩阵X第k行第m列的值，x_k，n为矩阵X第k行第n列的值，m,n=1,2,3,4，k＝1,2,3,…N，N为差异图像素点数；

4.2）计算相关系数矩阵ρ，以此来度量两个变量之间的相似性

$\mathrm{\&rho;}=\left[\begin{array}{cccc}{\mathrm{\&rho;}}_{\mathrm{1,1}}& {\mathrm{\&rho;}}_{\mathrm{1,2}}& ...& {\mathrm{\&rho;}}_{\mathrm{1,4}}\\ {\mathrm{\&rho;}}_{\mathrm{2,1}}& {\mathrm{\&rho;}}_{\mathrm{2,2}}& ...& {\mathrm{\&rho;}}_{\mathrm{2,4}}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\&rho;}}_{\mathrm{4,1}}& {\mathrm{\&rho;}}_{\mathrm{4,2}}& ...& {\mathrm{\&rho;}}_{\mathrm{4,4}}\end{array}\right]$

其中

${\mathrm{\&rho;}}_{a,b}=\frac{C_{a,b}}{\sqrt{C_{a,b}{C}_{a,b}}}$

Ρ_a，b表示矩阵ρ的第a行、第b列个元素，C_a，b为协方差矩阵C的第a行、第b列个元素,a，b=1,2,3,4；

4.3）初始化总的聚类层数为l=0，初始化正交Dirac基B⁽⁰⁾为4x4维的单位矩阵，C⁽⁰⁾=C；

4.4)设l＝l+1，寻找相关系数矩阵ρ中最大的值，将该值的座标序号分别标记为α和β，α>β；

4.5)利用得到的位置序号α和β，计算Jacobi旋转矩阵J：

其中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;}}^{(l-1)}-C_{\mathrm{\&beta;},\mathrm{\&beta;}}^{(l-1)}\&PlusMinus;\sqrt{{[C_{\mathrm{\&alpha;},\mathrm{\&alpha;}}^{(l-1)}-C_{\mathrm{\&beta;},\mathrm{\&beta;}}^{(l-1)}]}^2+4{\left[C_{\mathrm{\&alpha;},\mathrm{\&beta;}}^{(l-1)}\right]}^2}}{2C_{\mathrm{\&alpha;},\mathrm{\&beta;}}^{(l-1)}}\right]$ 且 $\left|{\mathrm{\&theta;}}_l\right|\&le;\frac{\mathrm{\&pi;}}{4},$

为第l-1层协方差矩阵C在α行α列的值，

为第l-1层协方差矩阵C在α行β列的值，

为第l-1层协方差矩阵C在β行β列的值；

4.6)利用旋转矩阵J得到去相关后的正交基B^(l)=B^(l-1)J，并更新协方差矩阵C^(l)=J^TC^(l-1)J，其中B^(l)为更新后的正交基，C^(l)为更新后的协方差矩阵，并根据协方差矩阵C^(l)更新相关系数矩阵ρ^(l);

4.7）判断l的值，如果l＝3执行步骤4.8），否则返回步骤4.4）；

4.8）将矩阵X沿正交基矩阵B⁽³⁾的转置方向进行投影，得到投影后的矩阵P_n：

P_n=X×B^(3)T；

4.9）将投影后的矩阵P_n的第一列P转化为与输入图像同等大小的矩阵，该矩阵即为融合后最终差异图D的灰度矩阵。

3.根据权利要求书1中所述的方法，其中步骤(5)所述的对所得的最终差异图D的灰度矩阵进行聚类分割，得到变化检测结果，按如下步骤进行：

5.1)将最终差异图D的灰度矩阵转化为一个列向量T=[t_i]，t_i为该向量第i个像素的灰度值，i=1,2,3,…N，N为总像素点数，将第i个像素属于第k类的隶属度记为u_ki，u_ki∈[0,1]且

k＝1,2，随机初始化隶属度矩阵U^b=[u_ki]的每一个元素，设定循环计数器b=0；

5.2)按如下公式计算第k类的聚类中心c_k：

$c_k=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ki}}^2{t}_i}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ki}}^2};$

5.3)按如下公式计算隶属度矩阵U^（b+1）:

U^（b+1)={u_ki}，

其中，k＝1,2, $u_{1i}=\frac{1}{1+{\left(\frac{t_i-c_1}{t_i-c_2}\right)}^2},$ $u_{2i}=\frac{1}{{\left(\frac{t_i-c_2}{t_i-c_1}\right)}^2+1};$

5.4)如果满足max{U^(b)-U^(b+1)}<0.001则转至步骤5.5），否则设置b=b+1,转至步骤5.2)；

5.5)通过比较c₁与c₂的大小确定变化类与非变化类，如果c₁≥c₂，则以c₁为聚类中心的那一类为变化类，以c₂为聚类中心的那一类为非变化类；如果c₁<c₂，以c₂为聚类中心的那一类为变化类，以c₁为聚类中心的一类为非变化类。

Images