CN102663730B - 基于Treelet和方向自适应滤波的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于Treelet和方向自适应滤波的遥感图像变化检测方法，其实现步骤为：(1)读入数据；(2)构造差异图像；(3)方向自适应滤波；(4)Treelet融合；(5)自适应阈值分类；(6)后处理。本发明既可以较好的保持变化区域的边缘信息，又能够减少变化检测结果中的漏检信息，具有较高的变化检测精度，可应用于湖泊水位的动态监测、农作物生长状态的动态监测、城区规划、军事侦察等领域。

Description

基于Treelet和方向自适应滤波的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，更进一步涉及一种基于Treelet和方向自适应滤波的遥感图像变化检测方法。该方法可应用于湖泊水位的动态监测、农作物生长状态的动态监测、城区规划、军事侦察等领域。

背景技术

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

在多时相遥感图像变化检测方法的研究中，常见的一种变化检测方法是先比较后分类法。它的优点在于简单易行，没有先分类后比较法存在的分类误差累计问题，但仍存在的不足是，该方法对分类阈值的选取十分敏感，如果分类阈值较大，会使变化检测结果存在较多漏检信息，如果分类阈值较小，会使变化检测结果存在较多虚警信息，降低变化检测精度。

光学遥感图像变化检测中图像滤波的方法已较为常见，如均值滤波、中值滤波、维纳滤波、形态学滤波等。这些滤波处理在一定程度上平滑了同质区域的噪声，但仍存在的不足是，该滤波处理难以在去除同质区噪声的同时保持图像的边缘信息。

西安电子科技大学在其专利申请“基于Treelets的遥感图像变化检测方法”(专利申请号：201110001584.0，公开号：CN102063720A)中提出了一种基于Treelets滤波的遥感图像变化检测方法。该方法虽然能够降低辐射校正和光照不均对变化检测结果的影响，但仍存在的不足是，Treelets滤波使得变化检测结果存在较多漏检信息，不能较好的保持变化区域边缘信息，降低了变化检测精度。

发明内容

本发明针对上述现有技术存在的不足，提出了一种基于Treelet和方向自适应滤波的遥感图像变化检测方法。本发明能够较好的保持变化区域的边缘信息，减少了检测结果中的漏检信息，具有较高的变化检测精度。

本发明实现上述目的的思路是：在构造完差异图像后，首先对差异图像进行方向自适应滤波，其次利用Treelet对滤波后的两幅图像和差异图像进行融合，然后对融合后的图像进行自适应阈值分类，最后对分类图进行基于面积阈值的后处理，得到变化检测结果图。

本发明的步骤包括如下：

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

(2)构造差异图像

将步骤(1)中的任意一幅遥感图像与另一幅遥感图像做减法运算，并将减法运算结果取绝对值，得到一幅差异图像。

(3)方向自适应滤波

3a)在差异图像中任意选取一个像素点，以该像素点为中心，以固定长度为边长，确定一个正方形邻域图像块；

3b)按照水平方向将正方形邻域图像块划分为2个方向模板，按照垂直方向将正方形邻域图像块划分为2个方向模板，按照对角方向将正方形邻域图像块划分为4个方向模板，按照原点位置将正方形邻域图像块划分为1个方向模板；

3c)按照标准差计算公式计算9个方向模板的标准差，对9个方向模板的标准差值按照从小到大的顺序进行排列；

3d)按照均值计算公式计算最小标准差值对应模板的灰度均值，将该值作为第一幅滤波图像在正方形邻域图像块中心处的灰度值；

3e)按照均值计算公式计算次小标准差值对应模板的灰度均值，将该值作为第二幅滤波图像在正方形邻域图像块中心处的灰度值；

3f)重复步骤3a)至步骤3e)，直至处理完差异图像中的全部像素点，得到一幅对应于最小标准差模板的滤波图像和一幅对应于次小标准差模板的滤波图像。

(4)Treelet融合

4a)将差异图像的像素按照从左到右、从上到下的顺序进行排列，组成一个列向量；

4b)将对应于最小标准差模板的滤波图像的像素按照从左到右、从上到下的顺序进行排列，组成一个列向量；

4c)将对应于次小标准差模板的滤波图像的像素按照从左到右、从上到下的顺序进行排列，组成一个列向量；

4d)将步骤4a)、4b)和4c)中的列向量依次按照从左到右的顺序进行排列，组成一个图像序列，对该图像序列进行Treelet变换，得到一个基矩阵；

4e)将图像序列向基矩阵投影，得到一幅融合图像。

(5)自适应阈值分类

5a)按照标准差计算公式计算融合图像的标准差；

5b)按照均值计算公式计算融合图像的均值；

5c)采用均值和标准差计算分类阈值；

5d)采用分类阈值对融合图像进行分类，得到一幅分类图像，将像素灰度值为1的像素点作为变化信息，将像素灰度值为0的像素点作为非变化信息。

(6)后处理

6a)将分类图像中像素灰度值为1的像素点按照八连通方式进行连接，得到分类图像的八连通区域；

6b)统计分类图像中八连通区域中像素灰度值为1的像素点个数；

6c)判断八连通区域像素点个数是否大于面积阈值，若满足，则将八连通区域视为变化信息区域，否则视为非变化信息区域，将该区域内的像素灰度值赋值为0。

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

第一，本发明利用方向自适应滤波进行图像滤波，克服了现有滤波技术难以在消除同质区噪声的同时保持图像细节信息的缺点，使得通过本发明获得的变化区域边缘信息保持较为完整。

第二，本发明采用自适应阈值分类和后处理方法检测变化区域，克服了现有图像分类技术存在较多漏检信息的缺点，提高了变化检测精度。

附图说明

图1为本发明的流程图；

图2为本发明中方向自适应滤波的9个方向模板示意图；

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

具体实施方式

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

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

步骤2，构造差异图像。

将步骤1中的任意一幅遥感图像与另一幅遥感图像做减法运算，并将减法运算结果取绝对值，得到一幅差异图像。

步骤3，方向自适应滤波。

选取正方形邻域图像块：在差异图像中任意选取一个像素点，其空间位置记为(m，n)，以该像素点为中心，以固定长度Nu为边长，确定一个正方形邻域图像块，其中，本发明实施例中固定长度Nu为5个像素点。

参照附图2，对正方形邻域图像块进行方向模板划分。附图2(a)和2(b)为正方形邻域图像块按照水平方向划分的2个方向模板。附图2(c)和2(d)为正方形邻域图像块按照垂直方向划分的2个方向模板。附图2(e)、2(f)、2(g)和2(h)为正方形邻域图像块对角水平方向划分的4个方向模板。附图2(i)为正方形邻域图像块按照原点位置划分的1个方向模板。在正方形邻域图像块中，将对应于空间位置(m，n)、(m-1，n-2)、(m-1，n-1)、(m，n-2)、(m，n-1)、(m+1，n-2)和(m+1，n-1)的7个像素点划分为第1个方向模板，如附图2(a)。将对应于空间位置(m，n)、(m-1，n+1)、(m-1，n+2)、(m，n+1)、(m，n+2)、(m+1，n+1)和(m+1，n+2)的7个像素点划分为第2个方向模板，如附图2(b)。将对应于空间位置(m，n)、(m-1，n-1)、(m-1，n)、(m-1，n+1)、(m-2，n-1)、(m-2，n)和(m-2，n+1)的7个像素点划分为第3个方向模板，如附图2(c)。将对应于空间位置(m，n)、(m+1，n-1)、(m+1，n)、(m+1，n+1)、(m+2，n-1)、(m+2，n)和(m+2，n+1)的7个像素点划分为第4个方向模板，如附图2(d)。将对应于空间位置(m，n)、(m，n-1)、(m-1，n)、(m-1，n-1)、(m-1，n-2)、(m-2，n-2)和(m-2，n-1)的7个像素点划分为第5个方向模板，如附图2(e)；将对应于空间位置(m，n)、(m，n+1)、(m-1，n)、(i-1，n+1)、(i-1，n+2)、(m-2，n+1)和(m-2，n+2)的7个像素点划分为第6个方向模板，如附图2(f)。将对应于空间位置(m，n)、(m，n-1)、(m+1，n-2)、(m+1，n-1)、(m+1，n)、(m+2，n-2)和(m+2，n-1)的7个像素点划分为第7个方向模板，如附图2(g)。将对应于空间位置(m，n)、(m，n+1)、(m+1，n)、(m+1，n+1)、(m+1，n+2)、(m+2，n+1)和(m+2，n+2)的7个像素点划分为第8个方向模板，如附图2(h)。将对应于空间位置(m，n)、(m-1，n-1)、(m-1，n)、(m-1，n+1)、(m+1，n-1)、(m+1，n)、(m+1，n+1)、(m，n-1)和(m，n+1)的9个像素点划分为第9个方向模板，如附图2(i)。

标准差排序：按照标准差计算公式计算9个方向模板的标准差，对9个方向模板的标准差值按照从小到大的顺序进行排列，将排序后的标准差对应的方向模板记为BM₁、BM₂、BM₃、BM₄、BM₅、BM₆、BM₇、BM₈、BM₉。

最小标准差对应方向模板滤波：按照下列均值计算公式计算最小标准差对应模板的灰度均值，将该值作为第一幅滤波图像在正方形邻域图像块中心处的灰度值：

ξ₁＝mean(BM₁)

其中，ξ₁为最小标准差对应方向模板的均值，mean为均值运算符，BM₁为最小标准差对应的方向模板。

次小标准差对应方向模板滤波：按照下列均值计算公式计算次小标准差对应模板的灰度均值，并将该值作为第二幅滤波图像在正方形邻域图像块中心处的灰度值：

ξ₂＝mean(BM₂)

其中，ξ₂为次小标准差对应方向模板的均值，mean为均值运算符，BM₂为次小标准差对应的方向模板。

重复步骤选取正方形邻域图像块、对正方形邻域图像块进行方向模板划分、标准差排序、最小标准差对应模板滤波和次小标准差对应模板滤波，直至处理完差异图像中的全部像素点，得到一幅对应于最小标准差模板的滤波图像和一幅对应于次小标准差模板的滤波图像。

步骤4，Treelet融合。

对差异图像的像素点按照从左到右、从上到下的顺序进行排列，组成一个列向量S₁。

对第一幅滤波图像中的像素点按照从左到右、从上到下的顺序进行排列，组成一个列向量S₂。

对第二幅滤波图像中的像素点按照从左到右、从上到下的顺序进行排列，组成一个列向量S₃。

对列向量S₁、S₂和S₃依次按照从左到右的顺序进行排列，组成一个图像序列X，对该图像序列进行Treelet变换，得到一个基矩阵。

在Treelet变换第l＝0分解层，将图像序列X初始化为和变量集初始化为δ＝{1，2，3}，正交基矩阵初始化为B₀＝[Φ_0，1，Φ_0，2，Φ_0，3]，其中，B₀是一个3×3的单位矩阵。

按照下式计算X⁽⁰⁾的协方差矩阵∑⁽⁰⁾：

${\mathrm{\Σ}}^{\left(0\right)}=\left[\begin{array}{ccc}{\mathrm{\Σ}}_{11}^{\left(0\right)}& {\mathrm{\Σ}}_{12}^{\left(0\right)}& {\mathrm{\Σ}}_{13}^{\left(0\right)}\\ {\mathrm{\Σ}}_{21}^{\left(0\right)}& {\mathrm{\Σ}}_{22}^{\left(0\right)}& {\mathrm{\Σ}}_{23}^{\left(0\right)}\\ {\mathrm{\Σ}}_{31}^{\left(0\right)}& {\mathrm{\Σ}}_{32}^{\left(0\right)}& {\mathrm{\Σ}}_{33}^{\left(0\right)}\end{array}\right]$

其中， 为求数学期望，η＝{1，2，3}和λ＝{1，2，3}为协方差矩阵∑⁽⁰⁾的位置索引。

按照下式计算相似度矩阵Θ⁽⁰⁾的每个元素：

${\mathrm{\Θ}}_{(\mathrm{\η},\mathrm{\λ})}^{\left(0\right)}=\left|\frac{{\mathrm{\Σ}}_{\mathrm{\η\λ}}^{\left(0\right)}}{\sqrt{{\mathrm{\Σ}}_{\mathrm{\η\η}}^{\left(0\right)}{\mathrm{\Σ}}_{\mathrm{\λ\λ}}^{\left(0\right)}}}\right|$

其中，为第0分解层相似度矩阵Θ⁽⁰⁾的第η行第λ列矩阵元素，为协方差矩阵∑⁽⁰⁾的第η行第λ列矩阵元素，为协方差矩阵∑⁽⁰⁾的第η行第η列矩阵元素，为协方差矩阵∑⁽⁰⁾的第λ行第λ列矩阵元素。

在Treelet变换第l＝1，2分解层，按照下式找到两个最相似的和变量：

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

其中，α和β为两个最相似的和变量，arg为取参数运算符，max为取最大值运算符，η和λ为第l-1分解层相似度矩阵Θ^(l-1)的位置索引，δ为和变量集，为第l-1分解层相似度矩阵Θ^(l-1)的第η行第λ列矩阵元素。

按照下式对最相似的两个和变量进行局部主成分分析：

其中，J为旋转矩阵，α和β为两个最相似的和变量，θ_l为旋转角，l为分解层数，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$

其中，θ_l为旋转角，l为分解层数，∑^(l)为第l分解层的协方差矩阵，J为旋转矩阵，T为转置符号，∑^(l-1)为第l-1分解层的协方差矩阵，为第l分解层协方差矩阵∑^(l)的第α行第β列矩阵元素，为第l分解层协方差矩阵∑^(l)的第β行第α列矩阵元素，α和β为两个最相似的和变量。

利用Jacobi旋转矩阵J更新第l分解层基矩阵B_l＝B_l-1J＝[Φ_l，1，Φ_l，2，Φ_l，3]和第l分解层X^(l)＝J^TX^(l-1)，其中，上标T表示转置。

在Treelet变换的第2分解层，提取基矩阵PB。

PB＝[Φ_2，1]

其中，PB为基矩阵，Φ_2，1为第2分解层基矩阵B₂的第1个列向量。

按照下式将图像序列X向基矩阵PB进行投影，得到一幅融合图像Fus。

Fus＝X·PB

其中，Fus为融合图像，X为图像序列，PB为基矩阵。

步骤5，自适应阈值分类。

按照常用的标准差计算公式计算融合图像的标准差sig。

按照常用的均值计算公式计算融合图像的均值ε。

按照下式，采用融合图像的标准差和均值计算融合图像的分类阈值：

$T=\left\{\begin{array}{cc}\mathrm{\&epsiv;}+\mathrm{sig},& \mathrm{sig}>\mathrm{th}\\ \mathrm{\&epsiv;}+2\&times;\mathrm{sig},& \mathrm{sig}\&le;\mathrm{th}\end{array}\right.$

其中，T为分类阈值，ε为融合图像的均值，sig为融合图像的标准差，th为一个先验阈值，本发明实施例中th＝15。

采用分类阈值对融合图像进行分类，得到一幅分类图像，将像素灰度值为1的像素点作为变化信息，将像素灰度值为0的像素点作为非变化信息。

步骤6，后处理。

按照八连通方式对分类图像中像素灰度值为1的像素点进行连接，得到分类图像的八连通区域。

统计分类图像中的八连通区域中像素灰度值为1的像素点个数。

判断八连通区域像素点个数是否大于面积阈值，若满足，则将八连通区域视为变化信息区域，否则视为非变化信息区域，将该区域内的像素灰度值赋值为0，其中，本发明实施例中面积阈值为75个像素点。

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

1.仿真条件

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

2.仿真内容

本发明仿真实验所用数据为两组真实遥感数据集。第一组真实遥感数据集是墨西哥郊外的两幅Landsat7ETM+第4波段光谱图像，两幅图像的大小均为512×512像素，两幅图像之间发尘的变化是由火灾破坏了大面积的当地植被所致，包括25589个变化像素和236555个非变化像素。第二组真实遥感数据集是1995年9月和1996年7月意大利撒丁岛Mulargia湖泊区域的两幅Landsat5TM+5波段图像，两幅图像的大小均为256×384像素，两幅图像之间发生的变化是由于湖水水位上升引起的，包括7613个变化像素和90691个非变化像素。

3.仿真效果分析

本发明中，将方向自适应滤波、自适应阈值分类法和基于面积阈值的后处理法相结合检测变化区域，不仅可以保持变化区域的边缘信息，而且降低了变化检测结果中的漏检信息，提高了变化检测精度。为了验证本发明的有效性和优越性，将本发明和上述盖超提出的变化检测方法进行比较。

附图3(a)为第一组真实遥感数据集采用盖超提出的变化检测方法的效果图，附图3(b)为第一组真实遥感数据集采用本发明的效果图，附图3(c)为第二组真实遥感数据集采用盖超提出的变化检测方法的效果图，附图3(d)为第二组真实遥感数据集采用本发明的效果图，四幅效果图中的白色区域均为变化区域，黑色区域均为非变化区域。从附图3可以看出，盖超提出的变化检测方法的检测结果图中存在较多漏检信息，平滑了变化区域的边缘信息，通过本发明得到的检测结果图存在较少漏检信息，较好的保持了变化区域的边缘信息。

本发明通过下表中的总错误数、虚警数和漏检数三个指标评价变化检测方法的好坏。

从上述表格可以看出，相对于盖超提出的变化检测方法的评价指标，本发明对第一组真实遥感数据集和第二组真实遥感数据集的变化检测结果均具有最少的总错误数，即具有最高的变化检测精度，并且具有最少的漏检数。

Claims (3)

1.基于Treelet和方向自适应滤波的遥感图像变化检测方法，包括如下步骤：

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

(2)构造差异图像

将步骤(1)中的任意一幅遥感图像与另一幅遥感图像做减法运算，并将减法运算结果取绝对值，得到一幅差异图像；

(3)方向自适应滤波

3a)在差异图像中任意选取一个像素点，以该像素点为中心，以固定长度为边长，确定一个正方形邻域图像块；

3b)按照水平方向将正方形邻域图像块划分为2个方向模板，按照垂直方向将正方形邻域图像块划分为2个方向模板，按照对角方向将正方形邻域图像块划分为4个方向模板，按照原点位置将正方形邻域图像块划分为1个方向模板；

3c)按照标准差计算公式计算9个方向模板的标准差，对9个方向模板的标准差值按照从小到大的顺序进行排列；

3d)按照均值计算公式计算最小标准差值对应模板的灰度均值，将该值作为第一幅滤波图像在正方形邻域图像块中心处的灰度值；

3e)按照均值计算公式计算次小标准差值对应模板的灰度均值，将该值作为第二幅滤波图像在正方形邻域图像块中心处的灰度值；

3f)重复步骤3a)至步骤3e)，直至处理完差异图像中的全部像素点，得到一幅对应于最小标准差模板的滤波图像和一幅对应于次小标准差模板的滤波图像；

(4)Treelet融合

4a)将差异图像的像素按照从左到右、从上到下的顺序进行排列，组成一个列向量；

4b)将对应于最小标准差模板的滤波图像的像素按照从左到右、从上到下的顺序进行排列，组成一个列向量；

4c)将对应于次小标准差模板的滤波图像的像素按照从左到右、从上到下的顺序进行排列，组成一个列向量；

4d)将步骤4a)、4b)和4c)中的列向量依次按照从左到右的顺序进行排列，组成一个图像序列，对该图像序列进行Treelet变换，得到一个基矩阵；

4e)将图像序列向基矩阵投影，得到一幅融合图像；

(5)自适应阈值分类

5a)按照标准差计算公式计算融合图像的标准差；

5b)按照均值计算公式计算融合图像的均值；

5c)采用均值和标准差计算分类阈值，所述的分类阈值的计算公式如下：

$T=\left\{\begin{array}{cc}\mathrm{\&epsiv;}+\mathrm{sig},& \mathrm{sig}>\mathrm{th}\\ \mathrm{\&epsiv;}+2\&times;\mathrm{sig},& \mathrm{sig}\&le;\mathrm{th}\end{array}\right.$

其中，T为分类阈值，ε为融合图像的均值，sig为融合图像的标准差，th为一个先验阈值，th=15；

5d)采用分类阈值对融合图像进行分类，得到一幅分类图像，将像素灰度值为1的像素点作为变化信息，将像素灰度值为0的像素点作为非变化信息；

(6)后处理

6a)将分类图像中像素灰度值为1的像素点按照八连通方式进行连接，得到分类图像的八连通区域；

6b)统计分类图像中八连通区域中像素灰度值为1的像素点个数；

6c)判断八连通区域像素点个数是否大于面积阈值，若满足，则将八连通区域视为变化信息区域，否则视为非变化信息区域，将该区域内的像素灰度值赋值为0。

2.根据权利要求1所述的基于Treelet和方向自适应滤波的遥感图像变化检测方法，其特征在于：步骤3a)所述的固定长度为5个像素点。

3.根据权利要求1所述的基于Treelet和方向自适应滤波的遥感图像变化检测方法，其特征在于：步骤6c)所述的面积阈值为75个像素点。