基于treelet融合和水平集分割的遥感图像变化检测

摘要

本发明公开了一种基于treelet融合和水平集分割的遥感图像变化检测方法，主要解决现有变化检测方法存在较多伪变化信息的问题。其实现过程是：输入两时相遥感图像，对每幅图像分别进行均值漂移滤波，得到两时相滤波后图像并分别对其进行3次不同层数下的二维平稳小波分解，对相同分解层数对应方向子带的小波系数矩阵做差；采用sobel算子对得到的水平、垂直方向小波系数差矩阵进行增强并进行二维小波逆变换重构；采用treelet算法融合不同分解层数的重构图像得到最终的差异图，对该差异图进行水平集分割得到变化检测结果。本发明能够有效提高变化检测结果的精度，同时较好的保持变化区域的边缘特征，可用于对自然灾害的分析、土地资源监测等领域。

说明书

技术领域

本发明属于图像处理技术领域，具体来说是一种基于treelet融合和水平集分割遥 感图像变化检测方法，可用于土地资源监测、自然灾害分析、城市发展规划等诸多领 域的遥感图像分析。

背景技术

遥感图像变化检测技术是遥感图像研究的重要组成部分，它对同一地点不同时期 的图像进行比较分析，根据图像之间的差异来得到人们需要的地物变化信息。

在遥感图像的变化检测方法中，最简单常见的方法是直接对图像灰度值进行差分 得到差异影像，利用阈值区分变化类和非变化类。然而不同时刻之间的遥感图像由于 光照、辐射等因素造成图像灰度值存在偏差，简单对图像的灰度做差得到差异图的变 化检测方法结果中会存在较多的伪变化信息，同时，如何精确的估计阈值一直是一个 瓶颈问题。

为了提高变化检测结果的精度，有关学者提出了许多改进的方法：Yakoub Bazi 等学者于2009年在会议文章“A Variational Level-Set Method for Unsupervised Change  Detection in Remote Sensing Images”中提出基于CV模型的变化检测方法，2010年作 者在上述会议文章的基础上进一步完善理论部分并增加实验数据，又发表了期刊文章 “Unsupervised Change Detection in Multispectral Remotely Sensed Imagery With Level  Set Methods”。该方法不需要估计变化阈值，但是CV模型容易受到噪声的影响从而 影响结果的准确性。Turgay Celik等学者在文章“A Bayesian approach to unsupervised  multiscale change detection in synthetic aperture radar image”中提出一种基于双树复小 波变换的变化检测方法。由于双树复小波变换的下采样特性，该方法首先要对原始图 像进行插值，插值方法的选择对结果有一定影响，此外，分别对低频系数矩阵和高频 系数矩阵进行分割容易累积分类误差，同样影响变化检测结果的精度。

发明内容

本发明的目的在于克服上述已有的遥感图像变化检测技术的不足，提出了一种基 于treelet融合和水平集分割的遥感图像变化检测方法，以较好的保持变化区域的边缘， 减少伪变化信息，提高变化检测的精度。

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

(1)对输入的两幅已配准的大小均为m×n的多时相遥感图像I¹和I²分别进行均 值漂移滤波，得到滤波后图像X¹和X²；

(2)对滤波后图像X¹和X²分别进行3次二维平稳小波分解，且每一幅滤波后 图像3次分解的最高分解层数分别为s(s＝1，...3)，一次二维平稳小波分解得到四个 小波系数矩阵，即一个低频系数矩阵和三个分别表示水平、垂直、对角方向的高频小 波系数矩阵；

(3)对滤波后图像X¹和X²相同分解层数下对应方向子带的小波系数矩阵做差， 得到每一个分解层s的低频小波系数差矩阵和三个分别表示水平、垂直、对角方向的 高频小波系数差矩阵；

(4)对步骤(3)中的水平方向小波系数差矩阵和垂直方向小波系数差矩阵利用 sobel算子进行增强，保持低频小波系数差矩阵和对角方向小波系数差矩阵不变；

(5)使用(3)中低频小波系数差矩阵、对角方向小波系数差矩阵和(4)中增 强后的水平、垂直方向小波系数差矩阵，对每一个分解层s进行二维逆平稳小波变换， 得到每个分解层s的重构图像RIs；

(6)对每个分解层s的重构图像RIs使用treelet变换进行融合，得到融合后的 差异图D；

(7)对融合后的差异图D进行水平集分割，得到变化检测结果图Z。

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

(1)本发明由于通过对两幅滤波后图像分解得到的小波系数矩阵做差并对差矩 阵重构来构造差异图，因而能减轻不同时刻之间的遥感图像由于光照、辐射等因素造 成图像灰度值存在偏差对结果造成的影响，从而能减少结果中的伪变化信息，提高变 化检测结果的精度；

(2)本发明由于采用sobel算子的水平、垂直方向模板对水平、垂直方向的小波 系数差矩阵进行增强，使得图像中几何形状的边缘特征更加突出；

(3)本发明采用treelet变换对三张逆平稳小波变换重构得到的差异图进行融合， 该变换能得到一棵反映数据内部结构的层级聚类树和一组多尺度的正交基，本发明利 用这组正交基将三张不同分解层数下重构得到的差异图进行融合，相比于直接对融合 前的差异图进行分割得到的结果，对融合后的差异图进行分割得到的结果精度更高；

(4)本发明由于采用均值漂移法对原始图像进行滤波，不仅可以滤除图像中的 噪声信息，还具有良好的边缘保持特性；

(5)本发明由于采用二维平稳小波分解，平稳小波具有冗余性和平移不变性， 因此不必对原始图像进行插值，避免了插值算法带来的误差。

附图说明

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

图2是本发明使用的两时相遥感图像及对应参考图；

图3是本发明通过treelet融合得到的差异图；

图4是用本发明及对比方法对模拟遥感数据集实验得到的变化检测结果图；

图5是用本发明及对比方法对真实遥感数据集实验得到的变化检测结果图。

具体实施方式

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

步骤1，对输入的两个不同时相的已配准且大小均为m×n的遥感图像I¹和I²，如 2(a)和图2(b)所示，分别进行均值漂移滤波，得到滤波后图像X¹和X²；

(1a)对输入图像的每一个像素点，按照下式计算均值漂移矢量m_h(x)的值：

$m_h\left(x\right)=\frac{{\Sigma }_{i=1}^{n}K\left(\frac{x-x_i}{h}\right)x_i}{{\Sigma }_{i=1}^{n}K\left(\frac{x-x_i}{h}\right)}$

式中，K(x)为高斯核函数，h为高斯核函数K(x)的带宽矢量，本发明取h＝{5，5}， x为当前像素点的灰度值，x_i为当前像素点带宽h范围内像素点的灰度值，n为落入当 前像素点带宽h范围内像素点的总个数；

(1b)用计算得到的m_h(x)值更新当前像素点的灰度值；

(1c)如果满足条件||m_h(x)-x||≤ε_s，ε_s为给定的允许误差，则停止迭代，将m_h(x) 值作为当前像素点的灰度值输出，此时m_h(x)值即为经过均值漂移滤波后该像素点的 灰度值，否则返回步骤(1a)；

(1d)对输入图像中的所有像素点均进行步骤(1a)-(1c)的操作，输出滤波后图像。

步骤2，对滤波后图像X¹和X²分别进行3次二维平稳小波分解，且每一幅滤波 后图像3次分解的最高分解层数分别为s(s＝1，...3)，一次二维平稳小波分解得到4 个小波系数矩阵，即1个低频系数矩阵和3个分别表示水平、垂直、对角方向的高频 小波系数矩阵，则每一幅图像的3次分解均能得到12个小波系数矩阵，包括3个低 频系数矩阵和9个高频系数矩阵。

步骤3，对滤波后图像X¹和X²相同分解层数下对应方向子带的小波系数矩阵做 差，得到每一个分解层s的低频小波系数差矩阵和三个分别表示水平、垂直、对角方 向的高频小波系数差矩阵，对滤波后图像X¹和X²相同分解层数下对应方向子带的小 波系数矩阵做差，包括对低频小波系数矩阵做差和对高频小波系数矩阵做差，即

$D_s^L=|L_s^1-L_s^2|$

$D_{\theta ,s}^H={|H}_{\theta ,s}^1-H_{\theta ,s}^2|$

式中和分别表示滤波后图像X¹和X²第s分解层得到的低频小波系数矩阵， 为和的低频小波系数差矩阵；和分别表示滤波后图像X¹和X²第s 分解层得到的高频小波系数矩阵，θ取值为0°，45°和90°，分别表示水平、对角和垂 直方向，为和的高频小波系数差矩阵。

步骤4，对水平方向小波系数差矩阵和垂直方向小波系数差矩阵利用sobel算子 进行增强，保持低频小波系数差矩阵和对角方向小波系数差矩阵不变。

(4a)采用sobel算子的水平方向模板与水平方向的高频小波系数差矩阵进 行卷积，得到增强后的水平方向小波系数差矩阵

(4b)采用sobel算子的垂直方向模板与垂直方向的高频小波系数差矩阵进 行卷积，得到增强后的垂直方向小波系数差矩阵

步骤5，对每个分解层s，采用步骤(3)中低频小波系数差矩阵对角方向 小波系数差矩阵以及步骤(4)中增强后的水平、垂直方向小波系数差矩阵和这四个小波系数矩阵进行二维逆平稳小波变换，得到三张重构后的图像RIs；

步骤6，对每个分解层s的重构图像RIs使用treelet变换进行融合，得到融合后 的差异图D，具体执行如下；

(6a)将重构后的图像RIs均变换为m×n大小的列向量RI1′，RI2′，RI3′，组 成初始样本X＝[RI1′，RI2′，RI3′]；

(6b)定义treelet变换的逐层聚类层数l＝0，1，...L-1，L为初始样本X中列向量的 个数，L＝3，在第0层，每个变量采用初始样本X中的列向量表示，初始化基矩阵 B⁰为L×L的单位矩阵及和变量下标集合Ω＝{1，2，...L}，计算样本X的初始协方差矩阵 C⁰和初始相关系数矩阵M⁰，计算公式如下：

C_ij＝E[(s_p-Es_p)(s_q-Es_q)]

$M_{\mathrm{ij}}=\frac{C_{\mathrm{ij}}}{\sqrt{C_{\mathrm{ii}}{C}_{\mathrm{jj}}}}$

式中C_ij表示初始协方差矩阵C⁰第i行第j列的计算值，s_p，s_q表示样本X中两 个不同的列向量，E表示求数学期望，M_ij表示初始相关系数矩阵M⁰第i行第j列的 计算值；

(6c)当逐层聚类层数l≠0时，寻找相关系数矩阵M^l中最大的两个值，将最大 值和次大值的对应位置序号分别记为α和β：

$(\alpha ,\beta )=\arg \underset{i,j\in \Omega }{\max }{M}_{\mathrm{ij}}^{l-1}$

这里i＜j，分别表示相关系数矩阵M^l中任意值的行和列，且只在和变量下标集合Ω内 进行；

(6d)计算Jacobi旋转矩阵J：

$J=\left(\begin{array}{ccccccc}1& L& 0& L& 0& L& 0\\ M& O& M& & M& & M\\ 0& L& c_n& L& s_n& L& 0\\ M& & M& O& M& & M\\ 0& L& -s_n& L& c_n& L& 0\\ M& & M& & M& O& M\\ 0& L& 0& L& 0& L& 1\end{array}\right)$

其中c_n＝cos(θ_l)，s_n＝sin(θ_l)；旋转角θ_l由C^l＝J^TC^l-1J及计算得到， J^T表示旋转矩阵J的转置：

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

式中，表示协方差矩阵C^l中行、列位置序号均为α的元素值，表示协方 差矩阵C^l中行、列位置序号均为β的元素值，表示协方差矩阵C^l中行位置序号为 α、列位置序号均为β的元素值；

(6e)由矩阵J更新去相关后的正交基B^l＝B^l-1J，协方差矩阵C^l＝J^TC^l-1J，其 中B^l-1为更新前的正交基，B^l为更新后的正交基，C^l-1为更新前的协方差矩阵，C^l为 更新后的协方差矩阵，并利用C^l更新相关系数矩阵M^l-1为M^l，M^l-1为更新前的相关 系数矩阵，M^l为更新后的相关系数矩阵：

$M_{\mathrm{ij}}^l=\frac{C_{\mathrm{ij}}^l}{\sqrt{C_{\mathrm{ii}}^l{C}_{\mathrm{jj}}^l}}$

其中为更新后相关系数矩阵M^l中第i行j列的更新值，为更新后协方差矩 阵C^l中第i行j列的值，为更新后协方差矩阵C^l中第i行i列的值，为更新后协 方差矩阵C^l中第j行j列的值；

(6f)定义更新后正交基矩阵B^l中位置序号为α和位置序号为β的两个列向量分 别为尺度函数φ_l和细节函数ψ_l，定义当前l层的尺度向量集合{φ_l}是尺度函数φ_l和上一 层的尺度向量集合{φ_l-1}的合集，将差变量的位置序号β从和变量的下标集合Ω中去 除，即Ω＝Ω{β}；

(6g)重复步骤(6c)至(6f)直至聚类的最高层l＝L-1，可得最终的正交基矩 阵B＝[φ_L-1，ψ₁，...ψ_L-1]，其中φ_L-1∈{φ_L-1}，ψ₁为第一层treelet变换得到的细节函数，ψ _L-1为最高层treelet变换得到的细节函数；

(6h)将初始样本X沿正交基矩阵B转置的方向进行投影，即P＝X×B^T，并将投 影向量P变回m×n大小的图像，得到融合后的差异图D，如图3所示。

步骤7，对融合后差异图D进行水平集分割，得到变化检测结果Z，具体执行如 下：

(7a)初始化水平集函数为ξ₀：

${\xi }_0(x,y)=\left(\begin{array}{cc}-\rho & (x_a,y_a)\in {\Lambda }_0-\partial {\Lambda }_0\\ 0& (x_a,y_a)\in \partial {\Lambda }_0\\ \rho & \Lambda -{\Lambda }_0\end{array}\right)$

式中Λ表示待分割差异图D的整个图像域，Λ₀为Λ中任意一个子区域，为该 子区域所有边界上的点，ρ为常数，本发明中取ρ＝6，(x_a，y_a)表示图像中任意像素点的 坐标；

(7b)对初始水平集函数ξ₀进行演化，使得总能量函数ε(ξ)不断减小直到达到规 定的迭代次数15时停止，此时水平集函数ξ对应的零水平集就是变化区域的边界，将 变化区域的像素灰度值赋值1，非变化区域的像素灰度值赋值0，得到的二值图就是 最终的变化检测结果图Z，如图4(c)所示，演化公式如下：

$\frac{\partial \xi }{\partial t}=\mu [\mathrm{\Delta \xi }-\mathrm{div}\left(\frac{▿\xi }{|▿\xi |}\right)]+\mathrm{\lambda \delta }\left(\xi \right)\mathrm{div}\left(g\frac{▿\xi }{|▿\xi |}\right)+\mathrm{vg\delta }\left(\xi \right)$

式中λ、μ、v均为常数，本发明中取λ＝5，μ＝0.2，v＝3，δ为狄拉克函数，Δ表示 拉普拉斯运算符，表示求梯度，div表示求散度，g为边界检测函数，

$g=\frac{1}{1+{|▿G_{\sigma }*D|}^2}$

式中G_σ表示标准差为σ高斯核，D为融合后差异图；

总能量函数ε(ξ)表示为内部能量函数与外部能量函数之和，如下式所示：

ε(ξ)＝μP(ξ)+ε_g，λ，v(ξ)

式中P(ξ)为内部能量函数，ε_g，λ，v(ξ)为外部能量函数，

$P\left(\xi \right)={\int }_{\Lambda }\frac{1}{2}{\left(\right|▿\xi |-1)}^{2}d{x}_{a}d{y}_a$

ε_g，λ，v(ξ)＝λL_g(ξ)+vA_g(ξ)

这里L_g(ξ)为水平集函数ξ对应的零水平集曲线长度计算函数，A_g(ξ)为曲线演化加 速函数：

$L_g\left(\xi \right)={\int }_{\Lambda }\mathrm{g\delta }\left(\xi \right)|▿\xi |d{x}_{a}d{y}_a$

A_g(ξ)＝∫_ΛgH(-ξ)dx_ady_a

式中H表示阶梯函数。

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

1.实验数据

本发明的实验数据为一组模拟数据集和一组真实多时相遥感数据集，其中模拟数 据集的原始图像为470×335像素大小的位于英国Feltwell村庄农区的ATM(Airborne  Thematic Mapper)第3波段图像，配准误差为1.5个像素左右，模拟变化图像是通过模 拟地球的天气变化和电磁波的辐射特性等因素影响并人工嵌入一些变化区域得到的， 真实遥感数据集为意大利撒丁岛的两时相Landsat5 TM+5波段图像，大小为 412×300，256灰度级。

2.实验方法

方法1：Turgay Celik等学者在文章“A Bayesian approach to unsupervised multiscale  change detection in synthetic aperture radar image”中提出的基于双树复小波变换的变化 检测方法；

方法2：Yakoub Bazi等学者2010年在文章“A Variational Level-Set Method for  Unsupervised Change Detection in Remote Sensing Images”中提出的基于CV模型的变 化检测方法；

方法3：本发明方法。

3.实验内容与结果分析

实验1，用不同方法对如图2(a)和图2(b)所示的两幅不同时相的模拟数据集进行变 化检测，变化检测参考图如图2(c)所示，对图2(a)和图2(b)的模拟遥感数据集进行变化 检测得到的结果如图4所示，其中图4(a)为现有方法1得到的变化检测结果图，图4(b) 为现有方法2得到的变化检测结果图，图4(c)为本发明方法得到的变化检测结果图，从 图4可以看出，现有方法1结果图中存在较多的伪变化信息，现有方法2结果图中存在 明显的漏检，本发明方法则比较接近参考图2(c)。

实验2，用不同方法对如图2(d)和图2(e)所示的两幅不同时相的真实遥感数据集进 行变化检测，变化检测参考图如图2(f)所示，对图2(d)和图2(e)的模拟遥感数据集进行 变化检测得到的结果如图5所示，图5(a)为现有方法1得到的变化检测结果图，图5(b) 为现有方法2得到的变化检测结果图，图5(c)为本发明方法得到的变化检测结果图，从 图5可以看出，现有方法1和现有方法2的结果图中都存在许多明显的噪声区域，本发 明方法的结果图和参考图2(f)最为接近。

对实验变化检测结果进行定量分析，即比较变化检测结果图与标准参考图，计算 总错误数、虚景数、漏检数等性能指标，结果如表1所示。

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

从表1上可以更为直观的看出，本发明方法能有效稳定地的获得遥感图像的变化 区域，与其它两种现有方法相比较变化检测结果的精度最高。