CN103093472A - 基于双字典交叉稀疏表示的光学遥感图像变化检测方法 - Google Patents

Abstract

本发明公开一种基于双字典交叉稀疏表示的光学遥感图像变化检测方法。主要解决现有方法检测结果不稳定和未能兼顾漏检和虚警的问题。实现步骤为：(1)读入两幅不同时相图像；(2)逐个像素构造特征向量；(3)逐个像素构造构造局部字典；(4)对两幅时相图像同一位置的像素，用时相1图像的局部字典稀疏表示时相2图像的特征向量，用时相2图像的局部字典稀疏表示时相1图像的特征向量；(5)利用所有特征向量的l₁范数构成两幅l₁范数图；(6)构造差异图；(7)对差异图进行最大熵阈值，得到初分类二值图；(8)对初分类二值图进行区域生长，得到最终变化检测结果。本发明的检测结果能够保持变化区域的边缘信息，减少伪变化区域，有效提高检测精度，可用于资源监测和灾害评估。

Description

基于双字典交叉稀疏表示的光学遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及光学遥感图像变化检测，具体地说是一种基于双字典交叉稀疏表示的光学遥感图像变化检测方法，适用于遥感图像处理和分析。

背景技术

遥感变化检测是检测同一地理位置在不同时期获取多幅遥感图像之间的变化信息，广泛应用于国民经济和国防建设的诸多领域，如国土资源和土地调查、森林资源监测、军事侦察、灾害预报与评估、国家重大生态工程监理等方面。

在已配准的多时相遥感图像变化检测方法中，常见的方法是先构造差异图，然后利用阈值或者分类方法确定变化与非变化类。该方法的优点是简单易行，不改变原始数据的信息，但它的缺点是对图像的几何校正、辐射校正等预处理要求较高，并且对精确阈值的自动选取和分类方法的要求更加苛刻，影响了检测的精度。

为了提高变化检测的精确度，至关重要的一步就是构建变化与非变化类别差异显著、能有效减少几何误差并降低辐射光照影响的差异图像。将两时相遥感图像信息交叉利用构建差异图像则是一个解决此问题的根本有效途径。目前此方面的研究才刚开始。西安电子科技大学在其专利申请“基于Treelets的遥感图像变化检测方法”（专利申请号:201110001584.O，公开号:CN102063720A）中提出了一种Treelets滤波和K-means聚类相结合的遥感图像变化检测方法。该方法能够减小因噪声、辐射和光照等因素对检测结果的影响，变化区域的一致性保持较好，减少了伪变化信息，提高了变化检测的精确性。不足之处是Treelets交叉滤波会导致检测结果中存在较多漏检信息，不能较好的保持变化区域边缘信息。2010年LamH.Nguyen和Trac D.Tran在文献“ASparsity-Driven Joint Image Registration And Change Detection Technique For SarImagery”(Proceedings 2010 IEEE International Conference on Acoustics,Speech andSignal Processing,ICASSP2010:2798-2801.)中引入字典和稀疏表示，提出了基于字典构造和稀疏逼近的变化检测方法。该方法利用一幅时相图中一个像素的若干局部图像块构造一个局部字典，用来稀疏逼近另一幅时相图的同一像素的图像块，其优点是利用原始图像数据自适应地构造局部字典、不需要训练字典，可以减小算法的执行时间和几何配准误差对变化检测的影响。不足之处是该方法将稀疏系数的l₁范数和逼近误差的l₂范数的乘积作为差异图的灰度幅值，而稀疏逼近误差的l₂范数对图像边缘敏感，对平滑区域不敏感，它在变化区域的边缘部分取得大值，在变化区域的非边缘部分取得小值，所以构造的差异图像并不合理，在变化区域的非边缘部分会产生很多的漏检；另外在非变化区域的图像边缘处也会产生虚警。

发明内容

本发明的目的在于针对上述遥感图像变化检测方法中的不足，提出了一种基于双字典交叉稀疏表示的光学遥感图像变化检测方法，以准确地检测出变化区域、减少伪变化信息，提高变化检测的精度。

本发明的实现方案，包括如下步骤：

（1）读入同一地区不同时刻获取的已配准的两幅遥感图像X₁和X₂，图像大小为I×J；

（2）对图像X₁和X₂分别进行N个像素的边界镜像拓展，得到边界拓展后的图像Y₁和Y₂，其中N为奇数，N∈{3,5,7}；

（3）对图像Y₁的像素(i+N,j+N)，以该像素为中心选取一个大小为N×N的正方形图像块，将该图像块中像素灰度值按照从左到右、再从上到下的顺序排列成一个列向量，该列向量作为像素(i+N,j+N)的特征向量

其中i和j为图像的行序号和列序号，i=1,2，...,I,j=1,2，...,J;

（4）对图像Y₂的像素(i+N,j+N)，以该像素为中心选取一个大小为N×N的正方形图像块，将该图像块中像素灰度值按照从左到右、再从上到下的顺序排列成一个列向量，该列向量作为像素(i+N,j+N)的特征向量

（5）对图像Y₁的像素(i+N,j+N)构造局部字典对图像Y₁的像素(i+N,j+N)，以该像素为中心选取一个大小为(2N+1)×(2N+1)的搜索窗

对搜索窗

中的像素，在保证不超出搜索窗

边界的情况下，逐个像素选取大小为N×N的正方形图像块，并将每一个正方形图像块中像素灰度值按照从左到右、再从上到下的顺序排列成一个列向量，每个列向量作为局部字典的一个原子，用总共(N+2)×(N+2)个原子构成一个维数大小为N²×(N+2)²的局部字典

（6）按照步骤（5），对图像Y₂的像素(i+N,j+N)构造局部字典

（7）利用局部字典

对特征向量

进行稀疏表示，使用正交匹配追踪算法计算稀疏表示系数

（8）利用局部字典

对特征向量

进行稀疏表示，使用正交匹配追踪算法计算稀疏表示系数

（9）计算稀疏表示系数

的l₁范数，并对l₁范数的结果取下整数，得到幅度值A₁(i,j)，由图像Y₁所有像素对应的A₁(i,j)得到l₁范数图A₁={A₁(i,j)|i=1,2，...,I，j=1,2,...,J}；

（10）计算稀疏表示系数的l₁范数，并对l₁范数的结果取下整数，得到幅度值A₂(i,j)，由图像Y₂所有像素对应的A₂(i,j)得到l₁范数图A₂={A₂(i,j)|i=1,2,...,I,j=,2，...,J}；

（11）构造差异图A：对两幅l₁范数图的对应像素点的幅度值进行减法运算，并对减法运算的结果取绝对值，得到一幅差异图像A=|A₁-A₂|；

（12）对差异图A进行最大熵阈值，得到初分类二值图C；

（13）对初分类二值图C进行区域生长，得到最终的变化检测结果图E。

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

1）本发明充分利用两个时相图像的局部信息给每一个像素构造两个不同的局部字典交互稀疏表示两个时相的图像块，可降低局部配准误差和弱噪声对变化检测结果的不良影响。

2）本发明不同于已有方法中利用一幅时相图像构造局部字典稀疏表示另外一幅时相图像的图像块，而是利用交互的稀疏表示方法，不会产生由第一时相构造局部字典稀疏表示第二时相的图像块与第二时相构造局部字典稀疏表示第一时相的图像块之间的变化检测结果不一致的情况，不会出现变化检测结果模棱两可、不稳定的现象，提高了检测的稳定性。

3）由于本发明中是对初分类二值图进行区域生长，初分类二值图中的像素灰度值只有0和1，因此，区域生长阈值可选初分类二值图中灰度值的0和1之间的任意一个值，最终的生长结果均相同且不受区域生长阈值选取的影响。因而，本发明方法在准确检测出变化信息的同时去掉了散杂的伪变化区域，弥补了现有技术中存在较多虚警的不足。

附图说明

图1是本发明的流程图；

图2是本发明使用的两组两时相遥感图像及其变化检测参考图；

图3是用本发明和现有方法对图2遥感图像的变化检测结果图。

具体实施方式

参照图1对本发明的步骤做进一步的详细描述。

步骤1，读入同一地区在不同时刻获取的已配准的两幅遥感图像X₁和X₂，图像大小均为I×J。

步骤2，对图像X₁和X₂分别进行N个像素的边界镜像扩展，得到图像Y₁和Y₂。

2a）将图像X₁的前N行围绕X₁的上边界作镜像反射扩展，将图像X₁的最后N行围绕X₁的下边界作镜像反射扩展，将X₁的前N列围绕X₁的左边界作镜像反射扩展，将X₁的最后N列围绕X₁的右边界作镜像反射扩展；，

2b)将图像X₁的前N行与前N列形成的图像块围绕X₁的左上角点作镜像反射扩展，将图像X₁的前N行与最后N列形成的图像块围绕X₁的右上角点作镜像反射扩展，将X₁的前N列与最后N行形成的图像块围绕X₁的左下角点作镜像反射扩展，将X₁的最后N列与最后N行形成的图像块围绕X₁的右下角点作镜像反射扩展，得到大小为(2N+I)×(2N+J)的图像Y₁；

2c)将图像X₂做与X₁同样地镜像扩展，得到大小为(2N+I)×(2N+J)的图像Y₂,N为图像边界镜像扩展的像素点个数，其中N∈{3,5,7}，在本发明实施例中选取N＝3。

步骤3，对图像Y₁的像素点(i+N,j+N)构造特征向量

对图像Y₁的像素点(i+N,j+N)，以该像素点为中心，选取一个大小为N×N正方形图像块，将该图像块中像素灰度值按照从左到右、再从上到下的顺序排列成一个维数大小为(N×N)×1列向量，该列向量作为像素(i+N,j+N)的特征向量

其中i和j为图像的行序号和列序号，i=1,2，...,I,j=1,2，...,J。

步骤4，对图像Y₂的像素点(i+N,j+N)构造特征向量

对图像Y₂的像素点(i+N,j+N)，以该像素点为中心，选取一个大小为N×N正方形图像块，将该图像块中像素灰度值按照从左到右、再从上到下的顺序排列成一个维数大小为(N×N)×1列向量，该列向量作为该像素的特征向量

步骤5，对图像Y₁的像素点(i+N,j+N)构造局部字典

5a）在图像Y₁中，以像素点(i+N,j+N)为中心，选取一个大小为(2N+1)×(2N+1)像素的搜索窗

5b）对搜索窗中的像素，在保证不超出搜索窗

边界的情况下，逐个像素选取大小为N×N的正方形图像块，并将每一个正方形图像块中像素灰度值按照从左到右、再从上到下的顺序排列成一个列向量，每个列向量作为局部字典的一个原子，用总共(N+2)×(N+2)个原子构成一个维数大小为N²×(N+2)²的局部字典

步骤6，对图像Y₂的像素点(i+N,j+N)构造局部字典

6a）在图像Y₂中，以像素点(i+N,j+N)为中心，选取一个大小为(2N+1)×(2N+1)像素的搜索窗

6b）对搜索窗

中的像素，在保证不超出搜索窗边界的情况下，逐个像素选取大小为N×N的正方形图像块，并将每一个正方形图像块中像素灰度值按照从左到右、再从上到下的顺序排列成一个列向量，每个列向量作为局部字典的一个原子，用总共(N+2)×(N+2)个原子构成一个维数大小为N²×(N+2)²的局部字典

步骤7，利用局部字典

对特征向量

进行稀疏表示，使用正交匹配追踪算法计算稀疏表示系数

7a）设置初始特征向量残差r₀为

初始列空间矩阵Ψ₀为空矩阵，初始列向量指标集v₀为空集，初始化迭代次数t为1，初始化稀疏度m为1，初始化字典矩阵Φ为局部字典并设列向量

为字典矩阵Φ的第λ列，λ=1,2,...,d,其中d为字典矩阵Φ的列数；

7b）根据下式找出在字典矩阵Φ中与特征向量残差r_t-1相关性最大的列向量的指标λ_t：

其中，

为使目标函数f(λ)取最大值时的指标λ的值，符号|·|为取绝对值操作，< >为取内积操作；

7c）更新列向量指标集v_t=v_t-1∪{λ_t}和列空间矩阵

其中，符号∪为取并集操作；

7d）利用最小二乘法计算稀疏系数b_t：

$b_t=\underset{b}{\arg \min }{\left|\right|y_1^{\mathrm{ij}}-{\mathrm{\&Psi;}}_{t}b\left|\right|}_2,$

其中，

为使目标函数f(b)取最小值时的稀疏矩阵b的值，符号||·||₂为取l₂范数操作；

7e）按照下式更新残差r_t：

$r_t=y_1^{\mathrm{ij}}-{\mathrm{\&Psi;}}_t{b}_t;$

7f）将t的值加1然后再赋给t，若t<m则返回步骤7b），否则将稀疏系数b_t的值作为最终的稀疏系数

的值。

步骤8，利用局部字典

对特征向量

进行稀疏表示，使用正交匹配追踪算法计算稀疏表示系数

8a）设置初始特征向量残差r₀为

初始列空间矩阵Ψ₀为空矩阵，初始列向量指标集v₀为空集，初始化迭代次数t为1，初始化稀疏度m为1，初始化字典矩阵Φ为局部字典

并设列向量

为字典矩阵Φ的第λ列，λ=1,2，...,d，其中d为字典矩阵Φ的总列数；

8b）根据下式找出在字典矩阵Φ中与特征向量残差r_t-1相关性最大的列向量的指标λ_t：

其中，

为使目标函数f(λ)取最大值时的指标λ的值，符号|·|为取绝对值操作，< >为取内积操作；

8c）更新列向量指标集v_t=v_t-1∪{λ_t}和列空间矩阵

其中，符号∪为取并集操作；

8d）利用最小二乘法计算稀疏系数b_t；

$b_t=\underset{b}{\arg \min }{\left|\right|y_1^{\mathrm{ij}}-{\mathrm{\&Psi;}}_{t}b\left|\right|}_2,$

其中，

为使目标函数f(b)取最小值时的稀疏矩阵b的值，符号||·||₂为取l₂范数操作；

8e）按照下式更新残差r_t：

$r_t=y_1^{\mathrm{ij}}-{\mathrm{\&Psi;}}_t{b}_t,$

8f）将t的值加1然后再赋给t，若t<m则返回步骤8b），否则将稀疏系数b_t的值作为最终的稀疏表示系数

的值。

步骤9，计算稀疏表示系数

的l₁范数，并对l₁范数的结果取下整数，得到幅度值A₁(i,j)，由图像Y₁所有像素对应的A₁(i,j)构成l₁范数图A₁={A₁(i,j)|i=1,2，...,I,j=1,2,...,J}。

按照下式计算图像Y₁的像素(i+N,j+N)对应的幅度值A₁(i,j)：

其中，符号

为取下整操作，符号||·||₁为取l₁范数操作。

步骤10，计算稀疏表示系数

的l₁范数，并对l₁范数的结果取下整数，得到幅度值A₂(i,j)，由图像Y₂所有像素对应的A₂(i,j)构成l₁范数图A₂={A₂(i,j)|i=1,2,...,I,j=1,2,...,J}。

按照下式计算图像Y₂的像素(i+N,j+N)对应的幅度值A₂(i,j)：

其中，符号

为取下整操作，符号||·||₁为取l₁范数操作。

步骤11，构造差异图A。

对两幅l₁范数图的对应像素点的幅度值进行相减并取绝对值，得到一幅差异图像A=|A₁-A₂|。

步骤12，对差异图像A计算最大熵阈值T,根据阈值T对差异图像A进行分割，得到初分类二值图C。

步骤13，对初分类二值图C进行区域生长，得到最终的变化检测结果图E。

13a)选取变化区域的生长种子：选取一个种子阈值T₁=λ×T，利用阈值T₁按照下式对差异图A进行分类，得到生长种子结果图D；

$D(i,j)=\left\{\begin{array}{cc}1,& A(i,j)\&GreaterEqual;T_1\\ 0,& A(i,j)<T_1\end{array}\right.,$

其中，T为差异图A的最大熵阈值，λ为阈值放大因子，1.5<λ<1.8，本发明实施例中选取λ=1.6；A(i,j)为生长种子结果图A像素(i,j)的灰度值，D(i,j)为生长种子结果图D像素(i,j)的灰度值，生长种子结果图D={D(i,j)|i=1,2，...,I,j=1,2，...,J}；

13b)将图D中所有灰度值为1的像素作为变化区域的生长种子S；

S={(i,j)|D(i,j)=1，i=1,2，...,I,j=1,2,...,J}，

13c)利用步骤13b)的生长种子S对初分类二值图C进行区域生长，得到最终的变化检测结果E；

其中区域生长算法中的区域生长阈值为α，0<α<1，由于初分类二值图C中的像素灰度值只有0和1，α选取小于1的任一正数，其获得的最终变化检测结果均相同，在本发明实施例中选取α=0.5。

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

1.实验数据及评价指标

本发明仿真实验所用数据为两组真实遥感数据集。第一组真实遥感数据集是意大利撒丁岛的两幅Landsat-5卫星TM第四波段光谱图像，两幅图像大小均为300×412像素，它们之间发生的变化是由湖中水位上升所致，包括7626个变化像素和115974个非变化像素；其两时相原始图像和变化检测参考图分别如图2(a)、图2(b)、图2(c)所示。第二组真实遥感数据集是墨西哥郊外的两幅Landsat-7EM+第四波段光谱图像，两幅图像的大小均为512×512像素，它们之间发生的变化是由火灾破坏了大面积的当地植被所致，包括25589个变化像素和236555个非变化像素；其两时相原始图像和变化检测参考图分别如图2(d)、图2(e)、图2(f)所示。

本发明中衡量变化检测算法的客观评价指标采用虚警数、漏检数、总错误数和正确率。

2.本发明使用的对比实验方法如下所述：

对比方法1，是Lam H.Nguyen等学者在文章“A Sparsity-Driven Joint ImageRegistration And Change Detection Technique For Sar Imagery”中提出的基于字典构造和稀疏逼近的变化检测方法，记为SDJIR法。由于Lam H.Nguyen等学者在文章中仅提出了构造差异图的方法而没有给出对差异图分类的方法，为了对比验证的公平性和客观性，其差异图的分类方法采用本发明方法中的最大熵阈值，并采用区域生长法做后处理。

对比方法2，是西安电子科技大学的专利申请“基于Tree1ets的遥感图像变化检测方法”（专利申请号:201110001584.O，公开号:CN102063720A）中提出的一种Tree1ets滤波和K-means聚类相结合的遥感图像变化检测方法，记为TC法。

2.实验内容和分析

为了验证本发明方法使用局部字典交互稀疏表示图像块的有效性，将本发明方法与SDJIR法进行对比；为了验证本发明方法交互稀疏表示和区域生长策略的有效性，将本发明方法与TC法进行对比。

仿真1，对附图2中两组真实遥感数据集，用SDJIR法进行变化检测，变化检测结果的评价指标如表1的第一行和第四行所示，其变化检测结果图如图3(b)和图3(e)所示。

仿真2，对附图2中两组真实遥感数据集，用TC法进行变化检测，变化检测结果的评价指标如表1的第二行和第五行所示，其变化检测结果图如图3(a)和图3(d)所示。

仿真3，对附图2中两组真实遥感数据集，用本发明方法进行变化检测，变化检测结果的评价指标如表1的第三行和第六行所示，其变化检测结果图如图3(c)和图3(f)所示。

表1两组真实遥感图像数据集采用不同方法变化检测结果的性能评价

从上表可以看出，在三种方法的变化检测结果中，本发明方法对第一组和第二组真实遥感数据集的整体评价是最好的。本发明方法对第一组真实遥感数据集的总错误数比TC法的总错误数少194个像素点，比SDJIR法的总错误数少5118个像素点；对第二组真实遥感数据集的总错误数比TC法的总错误数少925个像素点，比SDJIR法的总错误数少24476个像素点；SDJIR法的总错误数过多，是其它方法的约八倍，本发明方法的总错误数是已有方法中最少的。对第一组遥感数据集的漏检数比TC法的漏检数分别少154像素点，比SDJIR法的漏检数少6495个像素点；对第二组遥感数据集的漏检数比TC法的漏检数分别少2074像素点，比SDJIR法的漏检数少21774个像素点；SDJIR法的漏检数过多，是其它方法的近十倍，本发明方法的漏检数是已有方法中最少的。由此可以看出，本发明方法能够较为全面、准确地检测出变化信息，减少伪变化信息，较好的兼顾检测结果中的漏检和虚警，具有较高的检测精度。从两组实验数据集的效果图中可以看出，与TC法和SDJIR法相比，本发明方法能够较好的保持变化区域的边缘信息，误检的孤立像素点也是较少的。

Claims (4)

1.一种基于双字典交叉稀疏表示的光学遥感图像变化检测方法，包括如下步骤：

（1）读入同一地区不同时刻获取的已配准的两幅遥感图像X₁和X₂，图像大小为I×J；

（2）对图像X₁和X₂分别进行N个像素的边界镜像拓展，得到边界拓展后的图像Y₁和Y₂，其中N为奇数，N∈{3,5,7}；

（3）对图像Y₁的像素(i+N,j+N)，以该像素为中心选取一个大小为N×N的正方形图像块，将该图像块中像素灰度值按照从左到右、再从上到下的顺序排列成一个列向量，该列向量作为像素(i+N,j+N)的特征向量

其中i和j为图像的行序号和列序号，i=1,2,...,I，j＝1,2,...,J；

（4）对图像Y₂的像素(i+N,j+N)，以该像素为中心选取一个大小为N×N的正方形图像块，将该图像块中像素灰度值按照从左到右、再从上到下的顺序排列成一个列向量，该列向量作为像素(i+N,j+N)的特征向量

（5）对图像Y₁的像素(i+N,j+N)构造局部字典对图像Y₁的像素(i+N,j+N)，以该像素为中心选取一个大小为(2N+1)×(2N+1)的搜索窗

对搜索窗

中的像素，在保证不超出搜索窗

边界的情况下，逐个像素选取大小为N×N的正方形图像块，并将每一个正方形图像块中像素灰度值按照从左到右、再从上到下的顺序排列成一个列向量，每个列向量作为局部字典的一个原子，用总共(N+2)×(N+2)个原子构成一个维数大小为N²×(N+2)²的局部字典

（6）按照步骤（5），对图像Y₂的像素(i+N，j+N)构造局部字典

（7）利用局部字典

对特征向量

进行稀疏表示，使用正交匹配追踪算法计算稀疏表示系数

（8）利用局部字典

对特征向量

进行稀疏表示，使用正交匹配追踪算法计算稀疏表示系数

（9）计算稀疏表示系数

的l₁范数，并对l₁范数的结果取下整数，得到幅度值A₁(i,j)，由图像Y₁所有像素对应的A₁(i,j)得到l₁范数图A₁={A₁(i,j)i=1,2,...,I，j＝1,2,...,J}；

（10）计算稀疏表示系数

的l₁范数，并对l₁范数的结果取下整数，得到幅度值A₂(i,j)，由图像Y₂所有像素对应的A₂(i,j)得到l₁范数图A₂={A₂(i,j)|i=1,2,...,I，j=1,2,...,J}；

（11）构造差异图A：对两幅l₁范数图的对应像素点的幅度值进行减法运算，并对减法运算的结果取绝对值，得到一幅差异图像A=|A₁-A₂|；

（12）对差异图A进行最大熵阈值，得到初分类二值图C；

（13）对初分类二值图C进行区域生长，得到最终的变化检测结果图E。

2.根据权利要求1所述的基于双字典交叉稀疏表示的光学遥感图像变化检测方法，其中步骤(7)所述的使用正交匹配追踪算法计算稀疏表示系数

按如下步骤进行：

7a）设置初始特征向量残差r₀为

初始列空间矩阵Ψ₀为空矩阵，初始列向量指标集v₀为空集，初始化迭代次数t为1，初始化稀疏度m为1，初始化字典矩阵Φ为局部字典

并设列向量

为字典矩阵Φ的第λ列，λ=1,2,...,d，其中d为字典矩阵Φ的总列数；

7b）根据下式找出在字典矩阵Φ中与特征向量残差r_t-1相关性最大的列向量的指标λ_t：

其中，

为使目标函数f(λ)取最大值时的指标λ的值，符号|·|为取绝对值操作，< >为取内积操作；

7c）更新列向量指标集v_t=v_t-1∪{λ_t}和列空间矩阵

其中，符号∪为取并集操作；

7d）利用最小二乘法计算稀疏系数b_t；

$b_t=\underset{b}{\arg \min }{\left|\right|y_1^{\mathrm{ij}}-{\mathrm{\&Psi;}}_{t}b\left|\right|}_2,$

其中，

为使目标函数f(b)取最小值时的稀疏矩阵b的值，符号||·||₂为取l₂范数操作；

7e）按照下式更新残差r_t：

$r_t=y_1^{\mathrm{ij}}-{\mathrm{\&Psi;}}_t{b}_t,$

7f）将t的值加1然后再赋给t，若t<m则返回步骤7b），否则将稀疏系数b_t的值作为最终的稀疏系数

的值。

3.根据权利要求1所述的基于双字典交叉稀疏表示的光学遥感图像变化检测方法，其中步骤(8)所述的使用正交匹配追踪算法计算稀疏表示系数

按如下步骤进行：

8a）设置初始特征向量残差r₀为

初始列空间矩阵Ψ₀为空矩阵，初始列向量指标集v₀为空集，初始化迭代次数t为1，初始化稀疏度m为1，初始化字典矩阵Φ为局部字典

并设列向量

为字典矩阵Φ的第λ列，λ=1,2,...,d，其中d为字典矩阵Φ的总列数；

8b）根据下式找出在字典矩阵Φ中与特征向量残差r_t-1相关性最大的列向量的指标λ_t：

其中，

为使目标函数f(λ)取最大值时的指标λ的值，符号|·|为取绝对值操作，< >为取内积操作；

8c）更新列向量指标集v_t=v_t-1∪{λ_t}和列空间矩阵

其中，符号∪为取并集操作；

8d）利用最小二乘法计算稀疏系数b_t；

$b_t=\underset{b}{\arg \min }{\left|\right|y_1^{\mathrm{ij}}-{\mathrm{\&Psi;}}_{t}b\left|\right|}_2,$

其中，

为使目标函数f(b)取最小值时的稀疏矩阵b的值，符号||·||₂为取l₂范数操作；

8e）按照下式更新残差r_t：

$r_t=y_1^{\mathrm{ij}}-{\mathrm{\&Psi;}}_t{b}_t,$

8f）将t的值加1然后再赋给t，若t<m则返回步骤8b），否则将稀疏系数b_t的值作为最终的稀疏表示系数

的值。

4.根据权利要求1所述的基于双字典交叉稀疏表示的光学遥感图像变化检测方法，其中步骤(13)所述的对初分类二值图C进行区域生长，得到最终的变化检测结果图E，按如下步骤进行：

13a)选取变化区域的生长种子：选取一个种子阈值T₁=λ×T，利用阈值T₁按照下式对差异图A进行分类，得到生长种子结果图D；

$D(i,j)=\left\{\begin{array}{cc}1,& A(i,j)\&GreaterEqual;T_1\\ 0,& A(i,j)<T_1\end{array}\right.,$

其中，T为差异图A的最大熵阈值，λ为阈值放大因子，1.5<λ<1.8，A(i,j)为差异图A中像素(i,j)的灰度值，D(i,j)为生长种子结果图D中像素(i,j)的灰度值，生长种子结果图D={D(i,j)|i=1,2,...,I，j=1,2,...,J}；

13b)将结果图D中所有灰度值为1的像素作为变化区域的生长种子S；

S={(i,j)|D(i,j)=1，i=1,2,...,I，j=1,2,...,J}；

13c)利用步骤13b)的生长种子S对图C进行区域生长，得到最终的变化检测结果E，

其中区域生长算法中的区域生长阈值为α，0<α<1，由于初分类二值图C中的像素灰度值只有0和1，α选取小于1的任一正数，其获得的最终变化检测结果均相同。

Images