CN104700411A - 基于稀疏重建的两时相遥感影像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于稀疏重建的两时相遥感影像变化检测方法，主要解决现有技术去噪后图像信息丢失严重和变化检测精度低的问题。其实现过程为：输入两幅遥感影像并分别对其镜像延拓；将延拓后的两幅遥感影像X₁和X₂逐像素分别分割为重叠的图像块，构成两个训练样本集合Y₁和Y₂；通过稀疏字典训练得到两个训练字典D₁、D₂及其对应的稀疏表示系数S₁、S₂；由D₁S₁和D₂S₂重构出遥感影像和，进而得到差值图和比值图像，并分别进行三层小波分解；对分解后的高频带和低频带的小波系数采用不同的融合算子进行融合和逆变换得融合后的差异图X_D；对X_D进行聚类得到变化检测结果。本发明能有效地抑制噪声并保持影像边缘清晰度，提高变化检测的准确度。

Description

基于稀疏重建的两时相遥感影像变化检测方法

技术领域

本发明属于图像处理技术领域，特别涉及遥感影像变化检测方法，可用于对土地等地理资源的利用和植被覆盖情况的监测以及地震、洪水、泥石流、森林火灾等自然灾害灾情的评估等领域。

背景技术

遥感影像变化检测技术旨在检测同一地点不同时间的影像之间发生的变化。根据是否需要先验知识可将变化检测分为监督变化检测和非监督变化检测。由于实际中的先验知识比较难以获得，且需要人工参与，不能实现自动化；故非监督的可智能化的无监督变化检测方法被普遍采用。

无监督变化检测主要包括三个处理阶段：首先，对获得的影像进行预处理，包括辐射校正和几何配准等；然后，对校正过的影像进行比较，获得差异图像；最后，利用阈值或聚类的分割方法，通过对差异图像进行分析得到最终的变化检测结果图。其中遥感影像预处理的目的是将两幅遥感影像进行校准，预处理时首先需要对遥感影像进行去噪处理，对并尽可能地在抑制噪声的同时尽量地保留有用信号。

目前，对遥感影像的去噪方法主要分为基于空域和变换域两类。基于空域去噪主要是通过在图像上取一个滑动窗口对图像进行滤波，其滤波能力与滑动窗的大小成正比，随着窗口的增大其平滑效果变好，但是会导致图像边缘的模糊，使得细节信息丢失严重，在抑制噪声的同时很难保留边缘和纹理细节。基于变换域的去噪方法主要是利用多分辨分析的方法，常见的有基于小波变换和众多的后小波变换分析的方法。小波去噪效果的好坏取决于以下因素：小波基的选择、分解的层数、阈值函数的选取以及阈值估计方法的选取。其中最重要的环节是选择阈值函数和阈值选取。硬阈值函数的不连续性导致重构信号容易出现振铃、Pseudo-Gibbs现象；而软阈值函数虽整体连续性好，但估计值与实际值之间总存在恒定的偏差，具有一定的局限性，且可能会造成边缘模糊等失真现象。阈值的选取直接影响到最终的去噪效果，如何最大限度去除噪声的同时保留信号的原始特征是去噪过程中的一个难点，如果阈值选取过小，则会出现消噪不足，过多的保留了噪声，致使信号的弱特征成份被噪声淹没；如果阈值选取过大，则会出现过消噪，将信号中的弱特征成份误认为噪声消除。当信号所含噪声的水平不同，去噪时采用的阈值也应有所不同。

近年来，基于超完备稀疏分解的信号表示理论得到广泛关注，并取得巨大成果。其基本原理是利用超完备字典中的冗余基取代正交基，对字典的选择尽可能地包含分解信号的信息。对信号的稀疏分解就是从超完备字典当中选择出最佳线性组合的若干原子来表示信号，将该方法应用于遥感影像的去噪，能取得较好的去噪效果。

差异图像的生成也是变化检测过程中至关重要的一步，差异图像的质量直接决定后续分析处理的精确度，进而影响到整个变化检测系统的性能。差异图像的构造最为简单且常用的方法是差值法和比值法。差值法是对两时相影像对应像素点位置处的灰度值进行作差取绝对值，如果值为0则代表该像素点的信息未发生变化，否则判定该像素点的信息发生了变化。比值法是计算两时相影像对应像素点灰度值的比值，若一个像素点没有发生变化则该比值应接近于1，反之则远大于1或远小于1。

虽然这两种方法构造出的差异图像都可以用来有效地检测出变化信息，但仍存在各自的不足：差值法易受影像成像质量、波谱特性等客观条件影响，不能完全反应地面地物的辐射能量变化；变化检测结果漏检率较高；比值法在增强变化信息、抑制背景信息的同时可帮助减少大气条件造成的影响，但有时会过于夸大部分变化；变化检测结果误检率较高。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于稀疏重建的两时相遥感影像变化检测方法，以有效地抑制噪声，保持影像的边缘清晰度；同时尽可能降低漏检数及误检数，抑制背景信息，增强变化信息，提高变化检测的准确度。

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

(1)输入两幅大小均为I×J的遥感影像X_A和X_B，其中，I为遥感影像的行数，J为遥感影像的列数；

(2)对两幅遥感影像X_A和X_B分别进行N个像素的边界镜像延拓，得到边界延拓后的两幅遥感影像X₁和X₂，其中N为奇数，N∈{1,3,5,7,……}；

(3)以第一幅边界延拓后的遥感影像X₁的像素(i+N,j+N)为中心，选取一个大小为(2N+1)×(2N+1)的重叠图像块y_1k，并将其图像块拉成列向量作为像素(i+N,j+N)的特征向量y₁ ^ij，得到I×J个特征向量构成训练样本集合其中，i＝1,2,...,I,j＝1,2,...,J，k为正整数；

(4)以第二幅边界延拓后的遥感影像X₂的像素(i+N,j+N)为中心，选取一个大小为(2N+1)×(2N+1)的重叠图像块，并将其图像块拉成列向量作为像素(i+N,j+N)的特征向量y₂ ^ij，得到I×J个特征向量构成训练样本集合

(5)用训练样本集合Y₁进行稀疏字典训练，得到训练字典D₁及其对应的稀疏表示系数集合 $S_1={\left\{s_{1k}\right\}}_{k=1}^{I*J};$

(6)用训练样本集合Y₂进行稀疏字典训练，得到训练字典D₂及其对应的稀疏表示系数集合 $S_2={\left\{s_{2k}\right\}}_{k=1}^{I*J};$

(7)根据上述训练字典和系数表示集合，得到去噪后的两幅遥感影像的重叠图像块矩阵和其中，“*”表示矩阵相乘；

(8)将两幅遥感影像的重叠图像块矩阵和中的每个图像块放回与原始各自图像对应的像素位置处，并对这些图像块的重叠部分的灰度值进行平均，重构出两幅遥感影像和 ${\hat{X}}_2;$

(9)根据重构出的两幅遥感影像和得到两幅重构遥感影像的差值图像 $Y_D=|{\hat{X}}_1-{\hat{X}}_2|$ 和比值图像 $Y_R={\hat{X}}_1/{\hat{X}}_2;$

(10)对两幅重构遥感影像的差值图像Y_D和比值图像Y_R分别进行W层小波分解，得到各自图像在W个分解层上的高频带和低频带的小波系数，其中W＝3；

(11)对差值图像Y_D和比值图像Y_R的高频带和低频带的小波系数采用不同的融合算子进行融合，即对高频带的小波系数选取Y_D和Y_R这两幅图像区域标准差较小的小波系数作为高频带融合后的小波系数；对低频带的小波系数选取Y_D和Y_R这两幅图像边缘信息较大的小波系数作为低频带融合后的小波系数；

(12)对融合后的高频带和低频带的小波系数进行小波逆变换，得到融合后的差异图像X_D；

(13)运用加权的模糊C均值聚类算法对融合后差异图X_D的灰度值进行类别数为2的聚类，得到变化类和非变化类的二值变化检测结果图X_CD，分别用1、0表示，完成对两幅遥感影像变化信息的检测。

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

1)本发明由于对图像进行了重叠分块，利用重叠块的冗余信息，不仅保留了图像的结构信息，而且很好地修复了细节信息和边缘信息，有利于遥感影像的解译和后续的相关处理。

2)本发明由于采用了基于小波分解的图像融合策略，在时间和空间上分解频率，可以更容易地提取细节信息，保存图像细节信息，从而得到正确率更高的变化检测结果。

3)本发明由于对高频带和低频带的小波系数采用不同的融合规则分别进行融合，可以更大程度上地抑制变化区域的背景信息，增强变化区域的变化信息，使差值法的漏检率和比值法的误检率互补，获得较高的检测精度。

附图说明

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

图2是本发明中的稀疏字典训练示意图；

图3是本发明中对差值图像Y_D和比值图像Y_R进行融合的子流程图；

图4是本发明仿真使用的Bern地区的两幅SAR遥感影像；

图5是用本发明对输入的两幅原始遥感影像进行稀疏重构后的遥感影像图；

图6是本发明对Bern地区遥感影像进行变化检测的结果与标准参考图的对比图。

具体实施方式

下面结合附图对本发明的技术方案和效果作进一步的详细描述。

参照图1，本发明的实现步骤如下：

步骤1，输入原始遥感影像。

这两幅原始遥感影像为经过辐射校正和几何配准的同一地区在不同时刻的遥感影像，记t1时刻的遥感影像为X_A，t2时刻的遥感影像为X_B，其影像大小均为I×J。

步骤2，对该两幅遥感影像X_A和X_B分别进行N个像素的边界镜像延拓。

分别以影像X_A和X_B的边界像素为中轴，向外对称复制N个像素，得到边界延拓后的两幅遥感影像X₁和X₂，其中N为奇数，N∈{1,3,5,7,……}。

步骤3，构造训练样本集合Y₁和Y₂。

3a)以t1时刻的遥感影像X₁在(i+N，j+N)位置处的像素为中心，选取一个大小为(2N+1)×(2N+1)的重叠图像块y_1k，并将其图像块拉成列向量作为像素(i+N,j+N)的特征向量y₁ ^ij，得到I×J个特征向量构成t1时刻的训练样本集合其中，I为遥感影像的行数，J为遥感影像的列数，i＝1,2,...,I,j＝1,2,...,J，k为正整数；

3b)以t2时刻的遥感影像X₂在(i+N，j+N)位置处的像素为中心，选取一个大小为(2N+1)×(2N+1)的重叠图像块y_2k，并将其图像块拉成列向量作为像素(i+N,j+N)的特征向量y₂ ^ij，得到I×J个特征向量构成t2时刻的训练样本集合

本发明仿真实验中取N＝1，即将遥感影像分割为大小为3×3的重叠块，由于实验中采用的影像大小为256×256，则分割后的影像块数目为256×256＝4096个，且将每个影像块拉成列向量作为一个训练样本，得到4096个训练样本。

步骤4，利用t1时刻的训练样本集合Y₁和t2时刻的训练样本集合Y₂，分别进行稀疏字典训练，得到最终的两个训练字典D₁和D₂及其对应的两个稀疏表示系数矩阵S₁和S₂：

参照图2，本步骤的具体实现步骤如下：

4a)用t1时刻的训练样本集合进行稀疏字典训练

4a1)初始化：初始化t1时刻的字典D₁为大小为((2N+1)*(2N+1))×L的随机高斯矩阵，并对其列向量进行归一化处理，L为所要构造的字典原子个数，且(2N+1)²≤L≤I*J；

4a2)稀疏编码：固定t1时刻的字典D₁，采用正交匹配追踪算法按照如下公式计算每个图像块y_1k的稀疏表示系数s_1k：

$\begin{array}{cc}{s}_{1k}=\underset{s_{1k}}{\arg \min }{\left|\right|y_{1k}-D_1\×s_{1k}\left|\right|}_2^2& s.t.{\left|\right|s_{1k}\left|\right|}_0\≤H,\end{array}$

其中，1<k≤I*J且为正整数，H为稀疏度，即为稀疏表示系数s_1k的最大非0值的个数，进而可得到稀疏表示系数矩阵S₁；

4a3)更新t1时刻的字典D₁的第m列原子d_1m及其对应的稀疏表示系数矩阵S₁的第m行元素，1≤m≤L：

4a3.1)定义集合为与d_1m相关的样本集合Y₁中的索引，即点的索引值，其中，为与d_1m对应的稀疏表示系数矩阵S₁的第m行元素；

4a3.2)将索引值的值置0，计算除去成分在所有I*J个样本中造成的表示误差矩阵：其中，d_1t为t1时刻字典D₁的第t个原子，为系数矩阵S₁中与d_1t相应的第t行，1≤t≤L；

4a3.3)从E_m中只选择与集合ω_m对应的列构成残差矩阵

4a3.4)对残差矩阵进行奇异值分解，即其中，U为左奇异矩阵，Δ为奇异值矩阵，V^T为右奇异矩阵的转置；

4a3.5)用左奇异矩阵U的第一列替换t1时刻字典D₁的原子得到更新后的原子用右奇异矩阵V的第一列乘以奇异值矩阵的第一行第一列的元素值Δ(1,1)得到系数矢量其中，为的去掉零值元素后的收缩结果；

4a3.6)根据系数矢量和集合ω_m对应的索引值更新稀疏表示系数矩阵S₁的第m行元素

4a4)重复步骤4a3)，对t1时刻字典D₁进行逐列更新得到t1时刻的新字典

4a5)返回步骤4a2)，重复执行步骤4a2)-4a4)共200次，得到t1时刻最终的字典和其对应的稀疏表示系数矩阵S₁；

4b)用t2时刻的训练样本集合进行稀疏字典训练

4b1)初始化：初始化t2时刻的字典D₂为大小为((2N+1)*(2N+1))×L的随机高斯矩阵，并对其列向量进行归一化处理，L为所要构造的字典原子个数，且(2N+1)²≤L≤I*J；

4b2)稀疏编码：固定t2时刻的字典D₂，采用正交匹配追踪算法按照如下公式计算每个图像块y_2k的稀疏表示系数s_2k：

$\begin{array}{cc}{s}_{2k}=\underset{s_{2k}}{\arg \min }{\left|\right|y_{2k}-D_2\&times;s_{2k}\left|\right|}_2^2& s.t.{\left|\right|s_{2k}\left|\right|}_0\&le;H,\end{array}$

其中，1<k≤I*J且为正整数，H为稀疏度，即为稀疏表示系数s_2k的最大非0值的个数，进而得到稀疏表示系数矩阵S₂；

4b3)更新t2时刻的字典D₂的第m列原子d_2m及其对应的稀疏表示系数矩阵S₂的第m行，1≤m≤L：

4b3.1)定义集合为与d_2m相关的样本集合Y₂中的索引，即点的索引值，其中，为与d_2m对应的稀疏表示系数矩阵S₂的第m行元素；

4b3.2)将索引值的值置0，计算除去成分在所有I*J个样本中造成的表示误差矩阵：其中，d_2t为t2时刻字典D₂的第t个原子，为系数矩阵S₂中与d_2t相应的第t行，1≤t≤L；

4b3.3)从E_m中只选择与集合ω_m对应的列构成残差矩阵

4b3.4)对残差矩阵进行奇异值分解，即其中，U为左奇异矩阵，Δ为奇异值矩阵，V^T为右奇异矩阵的转置；

4b3.5)用左奇异矩阵U的第一列替换t2时刻字典D₂的原子得到更新后的原子用右奇异矩阵V的第一列乘以奇异值矩阵的第一行第一列的元素值Δ(1,1)得到系数矢量其中，为的去掉零值元素后的收缩结果；

4b3.6)根据系数矢量和集合ω_m对应的索引值更新稀疏表示系数矩阵S₂的第m行元素

4b4)重复步骤4b3)，对t2时刻字典D₂进行逐列更新得到t2时刻的新字典

4b5)返回步骤4b2)，重复执行步骤4b2)-4b4)共200次，得到t2时刻最终的字典和其对应的稀疏表示系数矩阵S₂。

步骤5，稀疏重构两幅遥感影像和

5a)根据t1时刻最终的字典和其对应的稀疏表示系数矩阵S₁得到t1时刻去噪后的重叠图像块矩阵其中，“*”表示矩阵相乘；

5b)根据t2时刻最终的字典和其对应的稀疏表示系数矩阵S₂得到t2时刻去噪后的重叠图像块矩阵 ${\hat{Y}}_2=D_2*S_2;$

5c)将t1时刻去噪后的重叠图像块矩阵的每列代表的图像块放回与其原始遥感影像X₁对应的像素位置处，并对这些图像块的重叠部分的灰度值进行平均，重构出t1时刻的遥感影像

5d)将t2时刻去噪后的重叠图像块矩阵的每列代表的图像块放回与其原始遥感影像X₂对应的像素位置处，并对这些图像块的重叠部分的灰度值进行平均，重构出t2时刻的遥感影像

步骤6，构造差值图像Y_D和比值图像Y_R。

6a)将重构出的两幅遥感影像和对应位置处像素的灰度值作差并取绝对值，得到两幅遥感影像的差值图像

6b)将重构出的两幅遥感影像和对应位置处像素的灰度值作比值处理，得到两幅遥感影像的比值图像

步骤7，对构造出的差值图像Y_D和比值图像Y_R分别进行W层小波分解，得到各自图像在W个分解层上的高频带和低频带的小波系数。

设W＝3，本步骤的具体步骤如下：

7a)在第一层，即W＝1的分解：

7a1)将差值图像Y_D分解为一个低频带子图像和三个高频带子图像其对应的小波系数记作和其中，ξ＝{LH,HL,HH}，LH表示垂直方向的高频带，HL表示水平方向的高频带，HH表示对角方向上的高频带；

7a2)将比值图像Y_R分解为一个低频带子图像和三个高频带子图像其对应的小波系数记作和

7b)在第二层，即W＝2的分解：

7b1)将差值图像Y_D经步骤7a1)得到的低频带子图像继续分解为一个低频带子图像和三个高频带子图像其对应的小波系数记作和其中，ξ＝{LH,HL,HH}；

7b2)将比值图像Y_R经步骤7a2)得到的低频带子图像继续分解为一个低频带子图像和三个高频带子图像其对应的小波系数记作和

7c)在第三层，即W＝3的分解：

7c1)将差值图像Y_D经步骤7b1)得到的低频带子图像继续分解为一个低频带子图像和三个高频带子图像其对应的小波系数记作和其中，ξ＝{LH,HL,HH}；

7c2)将比值图像Y_R经步骤7b2)得到的低频带子图像继续分解为一个低频带子图像和三个高频带子图像其对应的小波系数记作和

步骤8，对高频带和低频带的小波系数采用不同的融合算子进行处理。

8a)对低频带小波系数采用取边缘信息最大的方法得到低频带融合后的小波系数

$W_{\mathrm{LL}}^F(i,j)=\left\{\begin{array}{cc}{W}_{\mathrm{LL},3}^D(i,j),& I_{\mathrm{LL\&xi;},3}^D(i,j)\&GreaterEqual;I_{\mathrm{LL\&xi;},3}^R(i,j)\\ W_{\mathrm{LL},3}^R(i,j),& I_{\mathrm{LL},3}^D(i,j)\&le;{\mathrm{I\&sigma;}}_{\mathrm{LL},3}^R(i,j)\end{array}\right.,$

其中，表示融合后的小波系数在像素点(i,j)处即第i行第j列位置处的小波系数值；

指差值图像Y_D在第三分解层的低频带小波系数在像素点(i,j)处的值；

指比值图像Y_R在第三分解层的低频带小波系数在像素点(i,j)处的值；

表示差值图像第三分解层上的小波系数在像素点(i,j)处的总边缘信息，其计算公式为 $I_{\mathrm{LL},3}^D(i,j)={(F_1*W_{\mathrm{LL},3}^D)}^2(i,j)+{(F_2*W_{\mathrm{LL},3}^D)}^2(i,j)+{(F_3*W_{\mathrm{LL},3}^D)}^2(i,j);$

表示比值图像第三分解层上的小波系数在像素点(i,j)处的总边缘信息，其计算公式为： $I_{\mathrm{LL},3}^R(i,j)={(F_1*W_{\mathrm{LL},3}^R)}^2(i,j)+{(F_2*W_{\mathrm{LL},3}^R)}^2(i,j)+{(F_3*W_{\mathrm{LL},3}^R)}^2(i,j);$

式中“*”代表卷积运算符，F₁、F₂、F₃分别为垂直、水平和对角方向上的梯度算子，其值分别为：F₁＝{{-1,-1,-1},{2,2,2},{-1,-1,-1}}，F₂＝{{-1,2,-1},{-1,2,-1},{-1,2,-1}}，F₃＝{{-1,0,-1},{0,4,0},{-1,0,-1}}；

8b)对高频带小波系数采用取区域标准差最小的方法，得到高频带融合后的小波系数

$W_{\mathrm{\&xi;},W}^F(i,j)=\left\{\begin{array}{cc}{W}_{\mathrm{\&xi;},W}^R(i,j),& {\mathrm{\&sigma;}}_{\mathrm{\&xi;},W}^D(i,j)\&GreaterEqual;{\mathrm{\&sigma;}}_{\mathrm{\&xi;},W}^R(i,j)\\ W_{\mathrm{\&xi;},W}^D(i,j),& {\mathrm{\&sigma;}}_{\mathrm{\&xi;},W}^D(i,j)\&le;{\mathrm{\&sigma;}}_{\mathrm{\&xi;},W}^R(i,j)\end{array}\right.,$

其中，ξ＝{LH,HL,HH}，W＝{1,2,3}；表示第W分解层上的像素点(i,j)位置处融合后的小波系数值；

${\mathrm{\&sigma;}}_{\mathrm{\&xi;},W}^D(i,j)=\sqrt{\frac{1}{3*3}\underset{m=-1}{\overset{1}{\mathrm{\&Sigma;}}}\underset{n=-1}{\overset{1}{\mathrm{\&Sigma;}}}{[W_{\mathrm{\&xi;},W}^D(i+m,j+n)-\mathrm{\&mu;}(i,j)]}^2}$ 表示差值图像在第W分解层上的在像素点(i,j)处3*3区域窗口内的标准差，式中，为3*3区域窗口内小波系数的均值，为3*3区域窗口内像素点(i+m,j+n)处的小波系数值；

${\mathrm{\&sigma;}}_{\mathrm{\&xi;},W}^R(i,j)=\sqrt{\frac{1}{3*3}\underset{m=-1}{\overset{1}{\mathrm{\&Sigma;}}}\underset{n=-1}{\overset{1}{\mathrm{\&Sigma;}}}{[W_{\mathrm{\&xi;},W}^R(i+m,j+n)-\mathrm{\&mu;}(i,j)]}^2}$ 表示比值图像在第W分解层上位于像素点(i,j)处的3*3区域窗口内的标准差。

步骤9，对融合后高频带和低频带的小波系数进行小波逆变换，得到融合后的差异图像X_D。

上述步骤7-步骤9的内容如图3所示。

步骤10，运用加权的模糊C均值聚类算法WFCM对融合后的差异图X_D的灰度值进行类别数为2的聚类，得到分别用1、0表示的变化类和非变化类的二值变化检测结果图X_CD，完成对两时相遥感影像的变化信息的检测。

本发明的效果可以通过以下仿真实验作进一步的说明。

1.仿真条件

本发明是在中央处理器为Intel(R)Core i3-5302.93GHZ、内存4G、WINDOWS 7操作系统上，运用MATLAB软件进行的仿真。

本发明的仿真数据如图4所示的两幅Bern城市的大小均为301×301的ERS-2SAR遥感影像，图4(a)和图4(b)的拍摄时间分别为1999年4月和1999年5月，图6(a)为标准的变化检测参考图。

2.仿真参数

对具有标准的变化检测参考图的实验数据可进行定量的变化检测结果分析，主要的评价指标有：

漏检数：统计实验结果图中发生变化区域的像素个数，与参考图中变化区域的像素个数进行对比，把参考图中发生变化但实验结果图中检测为未变化的像素个数，称为漏检数FN；

误检数：统计实验结果图中未发生变化区域的像素个数，与参考图中未发生变化区域的像素个数进行对比，把参考图中未发生变化而实验结果图中检测为变化的像素个数，称为误检个数FP；

正确分类的概率PCC：其中TP、TN分别为正确检测到非变化及变化的像素个数。

衡量检测结果图与参考图一致性的Kappa系数：其中，这里，N表示像素总个数，Nc和Nu分别表示实际的变化像素数和未变化像素数。

3.仿真内容

用本发明对图4所示的两幅SAR遥感影像进行变化检测：即先对这两幅遥感影像分别进行稀疏字典训练并重构，得到去噪后的遥感影像和其结果分别如图5(a)和5(b)所示；再根据去噪后的遥感影像和构造差异图像，并用WFCM聚类方法对差异图像进行聚类生成最终的变化检测结果图如图6(b)所示。

4.仿真实验结果及分析

将图6(b)所示的变化检测结果与图6(a)所示的标准变化检测参考图进行比较，可以看出，用本发明方法得到的仿真实验结果具有较好的主观视觉效果，变化检测结果图噪声点较少。

进一步以图6(a)的标准变化检测参考图为基准对图6(b)所示的变化检测结果进行定量的指标评价，其指标如下表1所示。

表1 本发明变化检测结果的定量指标评价

从表1中可以看出，本发明方法的正确分类概率PCC和Kappa系数值都比较理想，误检数FP和漏检数FN比较低，检测效果优良。通过这些指标可以定量地看出本发明应用于遥感影像变化检测产生了较好的效果。

Claims (4)

1.一种基于稀疏重建的两时相遥感影像变化检测方法，其特征在于：包括如下步骤：

(1)输入两幅大小均为I×J的遥感影像X_A和X_B，其中，I为遥感影像的行数，J为遥感影像的列数；

(2)对两幅遥感影像X_A和X_B分别进行N个像素的边界镜像延拓，得到边界延拓后的两幅遥感影像X₁和X₂，其中N为奇数，N∈{1,3,5,7,……}；

(3)以第一幅边界延拓后的遥感影像X₁的像素(i+N,j+N)为中心，选取一个大小为(2N+1)×(2N+1)的重叠图像块y_1k，并将其图像块拉成列向量作为像素(i+N,j+N)的特征向量y₁ ^ij，得到I×J个特征向量构成训练样本集合其中，i＝1,2,...,I,j＝1,2,...,J，k为正整数；

(4)以第二幅边界延拓后的遥感影像X₂的像素(i+N,j+N)为中心，选取一个大小为(2N+1)×(2N+1)的重叠图像块y_2k，并将其图像块拉成列向量作为像素(i+N,j+N)的特征向量y₂ ^ij，得到I×J个特征向量构成训练样本集合

(5)用训练样本集合Y₁进行稀疏字典训练，得到训练字典D₁及其对应的稀疏表示系数集合 $S_1={\left\{s_{1k}\right\}}_{k=1}^{I*J};$

(6)用训练样本集合Y₂进行稀疏字典训练，得到训练字典D₂及其对应的稀疏表示系数集合 $S_2={\left\{s_{2k}\right\}}_{k=1}^{I*J};$

(7)根据上述训练字典和稀疏表示系数集合，得到去噪后的两幅遥感影像的重叠图像块矩阵和其中，“*”表示矩阵相乘；

(8)将两幅遥感影像的重叠图像块矩阵和中的每个图像块放回与原始各自图像对应的像素位置处，并对这些图像块的重叠部分的灰度值进行平均，重构出两幅遥感影像和

(9)根据重构出的两幅遥感影像和得到两幅重构遥感影像的差值图像 $Y_D=|{\hat{X}}_1-{\hat{X}}_2|$ 和比值图像 $Y_R={\hat{X}}_1/{\hat{X}}_2;$

(10)对两幅重构遥感影像的差值图像Y_D和比值图像Y_R分别进行W层小波分解，得到各自图像在W个分解层上的高频带和低频带的小波系数，其中W＝3；

(11)对差值图像Y_D和比值图像Y_R的高频带和低频带的小波系数采用不同的融合算子进行融合，即对高频带的小波系数选取Y_D和Y_R这两幅图像区域标准差较小的小波系数作为高频带融合后的小波系数；对低频带的小波系数选取Y_D和Y_R这两幅图像边缘信息较大的小波系数作为低频带融合后的小波系数；

(12)对融合后的高频带和低频带的小波系数进行小波逆变换，得到融合后的差异图像X_D；

(13)运用加权的模糊C均值聚类算法对融合后差异图X_D的灰度值进行类别数为2的聚类，得到变化类和非变化类的二值变化检测结果图X_CD，分别用1、0表示，完成对两幅遥感影像变化信息的检测。

2.根据权利要求1所述的基于稀疏重建的两时相遥感影像变化检测方法，其特征在于：步骤(2)所述的对两幅遥感影像X_A和X_B分别进行N个像素的边界镜像延拓，是分别以第一幅遥感影像X_A和第二幅遥感影像X_B的边界像素为中轴，向外对称复制N个像素，得到边界延拓后的两幅遥感影像X₁和X₂。

3.根据权利要求1所述的基于稀疏重建的两时相遥感影像变化检测方法，其特征在于：步骤(5)所述的用训练样本集合Y₁进行稀疏字典训练，按如下步骤进行：

(5a)初始化字典D₁为大小为((2N+1)*(2N+1))×L的随机高斯矩阵，并对其列向量进行归一化处理，其中L为所要构造的字典原子个数，且(2N+1)²≤L≤I*J；

(5b)固定字典D₁，采用正交匹配追踪算法计算每个图像块y_1k的稀疏表示系数s_1k：

$\begin{array}{cc}{s}_{1k}=\underset{s_{1k}}{\mathrm{atg}\min }{\left|\right|y_{1k}-D_1\&times;s_{1k}\left|\right|}_2^2& s.t.{\left|\right|s_{1k}\left|\right|}_0\&le;H,\end{array}$

其中，1<k≤I*J且为正整数，H为稀疏度，即稀疏表示系数s_1k的最大非0值个数；

(5c)由稀疏表示系数s_1k，得稀疏表示系数矩阵S₁；

(5d)更新字典的第m列d_1m及其对应的稀疏表示系数矩阵S₁的第m行，1≤m≤L：

(5d1)定义集合为与d_1m相关的样本集合Y₁中的索引，即点的索引值，其中，为与d_1m对应的稀疏表示系数矩阵S₁的第m行元素；

(5d2)将的值置0，计算除去成分在所有I*J个样本中造成的表示误差矩阵：其中，d_1t为字典D₁的第t个原子，为系数矩阵S₁中与d_1t相应的第t行；

(5d3)从E_m中只选择与ω_m对应的列构成残差矩阵

(5d4)对进行奇异值分解，即其中，U为左奇异矩阵，Δ为奇异值矩阵，V^T为右奇异矩阵的转置；

(5d5)用左奇异矩阵U的第一列更新字典的原子用右奇异矩阵V的第一列乘以Δ(1,1)更新系数矢量其中为中去掉零输入后的收缩结果；

(5e)重复步骤(5d)，对字典D₁进行逐列更新得到新的字典

(5f)返回步骤(5b)，重复执行步骤(5b)-(5e)共200次，得到优化的字典D₁和与训练样本集Y₁对应的稀疏表示系数矩阵S₁。

4.根据权利要求1所述的基于稀疏重建的两时相遥感影像变化检测方法，其特征在于：步骤(6)所述的用训练样本集合Y₂进行稀疏字典训练，按如下步骤进行：

(6a)初始化字典D₂为大小为((2N+1)*(2N+1))×L的随机高斯矩阵，并对其列向量进行归一化处理，其中L为所要构造的字典原子个数，且(2N+1)²≤L≤I*J；

(6b)固定字典D₂，采用正交匹配追踪算法计算每个图像块y_2k的稀疏表示系数s_2k：

$\begin{array}{cc}{s}_{2k}=\underset{s_{2k}}{\mathrm{atg}\min }{\left|\right|y_{2k}-D_2\&times;s_{2k}\left|\right|}_2^2& s.t.{\left|\right|s_{2k}\left|\right|}_0\&le;H,\end{array}$

其中，1<k≤I*J且为正整数，H为稀疏度，即稀疏表示系数s_2k的最大非0值个数；

(6c)由稀疏表示系数s_2k，得稀疏表示系数矩阵S₂；

(6d)更新字典的第m列d_2m及其对应的稀疏表示系数矩阵S₂的第m行，1≤m≤L：

(6d1)定义集合为与d_2m相关的样本集合Y₂中的索引，即点的索引值，其中，为与d_2m对应的稀疏表示系数矩阵S₂的第m行元素；

(6d2)将的值置0，计算除去成分在所有I*J个样本中造成的表示误差矩阵：其中，d_2t为字典D₂的第t个原子，为系数矩阵S₂中与d_2t相应的第t行；

(6d3)从E_m中只选择与ω_m对应的列构成残差矩阵

(6d4)对进行奇异值分解，即其中，U为左奇异矩阵，Δ为奇异值矩阵，V^T为右奇异矩阵的转置；

(6d5)用左奇异矩阵U的第一列更新字典的原子用右奇异矩阵V的第一列乘以Δ(1,1)更新系数矢量其中为中去掉零输入后的收缩结果；

(6e)重复步骤(6d)，对字典D₂进行逐列更新得到新的字典

(6f)返回步骤(6b)，重复执行步骤(6b)-(6e)共200次，得到优化的字典D₂和与训练样本集Y₂对应的稀疏表示系数矩阵S₂。