CN110503631B - 一种遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种用低秩先验学习判别字典的遥感图像变化检测方法，主要提升无监督遥感图像变化检测的准确率，并通过学习变化判别字典和不变化判别字典代替自动阈值的选择。其步骤为：(1)分别输入待检测的两幅多时相遥感图像；(2)超像素分割差分图，得到超像素总样本；(3)对超像素总样本学习全局字典和稀疏系数；(4)根据稀疏系数矩阵的低秩表示挑选样本；(5)用变化样本和不变化样本分别学字典；(6)计算每个像素点邻域块在判别字典下的重构误差，比较重构误差得到变化区域二值图。本发明能无监督选择变化样本和不变化样本进行变化检测，可用于灾害评估，城市扩张检测，土地覆盖、利用检测等遥感领域。

Description

一种遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，更进一步涉及遥感图像中的变化检测方法，可用于灾害评估、城市扩张检测、土地覆盖等遥感领域。

背景技术

变化检测是对从不同时间获取的同一地理位置的遥感图像监测其变化区域的技术。随着空间分辨率的提高和卫星传感器监测范围的扩大，如何充分利用有限的光谱信息来突出变化区域，提升变化检测的准确性，对变化检测提出了巨大的挑战。

目前遥感图像变化检测的方法主要分为三类，监督、半监督和无监督方法，三类的主要区别在于检测时是否用到数据标签。监督方法需要用带标签的数据训练模型来实现变化检测，半监督方法在输入带标签数据训练的同时也输入无标签数据进行训练，无监督方法的训练数据都是无标签的。监督方法因为训练时使用大量的标签所以检测结果较好，但是数据标注工作对于遥感图像变化检测任务来说很困难，只能人工标记；在半监督的变化检测方法中，引入无标签信息来弥补标签信息的不足，但是训练集中难以分类的训练样本会削弱学习模型的辨识能力；无监督方法不需要标记好的数据来训练，但是变化检测结果上还有提升的空间，例如变化矢量分析法(CVA)、主成分分析法(PCA)都是无监督方法，原理简单，但是对噪声敏感，得到的结果往往不完整，而且存在细小的伪变化区域。

变化检测又分为深度神经网络获取图像特征检测变化以及采用稀疏编码的方式学习字典检测图像变化。其中基于深度学习的方法，由于是针对遥感图像的检测，数据集达不到训练量的要求，因此无法在实际应用中得到很好的结果。在有监督的变化检测方法中，因为加入了优先级信息，使得精度得到了很大的提升，但是标签信息标注困难；在半监督的变化检测方法中，引入无标签信息来弥补标签信息的不足，但是训练集中难以分类的训练样本会削弱学习模型的辨识能力。基于稀疏编码的遥感图像变化检测方法在检测结果和使用范围上都要优于基于深度学习的遥感图像变化检测方法。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种无监督的用低秩先验学习判别字典的遥感图像变化检测方法，提升变化检测的准确率，减少错检、漏检。

为实现上述目的，本发明采用的技术方案是：

一种遥感图像变化检测方法，该方法对两张同一地点不同时相的遥感图像差分图进行超像素分割，基于分割对象学习一个全局字典和稀疏表示；对稀疏系数矩阵进行低秩表示挑选变化的样本及不变化的样本；分别用变化的和不变化的样本学习变化字典和不变化字典；计算每个像素的重构误差，得到变化区域二值图，具体步骤包括如下：

(1)输入两幅不同时相的遥感图像：时相一X和时相二Y；

(2)计算差分图DI，DI＝X-Y,对差分图DI进行超像素分割，得到N个超像素块，从而得到超像素总样本M，具体方法为：

(2a)对DI、X、Y分别计算每一个超像素块各个波段上的平均像素值，公式为：

其中，i代表波段数，j代表超像素块序列，

分别代表在DI、X、Y上每一个超像素块各个波段上的平均像素值，R_j代表第j个超像素块，N_j表示每个超像素块中像素的总个数，x代表像素点，DI_i(x)、X_i(x)、Y_i(x)分别代表遥感图像DI、X、Y第i个波段上在x像素点的像素值；

(2b)对DI、X、Y，按照超像素分割的结果，对每一个超像素块各个波段上的像素值平均，按波段顺序排列成列向量，三个列向量合成一个列向量sp_j后，各个超像素列向量按超像素序号排列得到超像素总样本M，其中，sp_j表示如下：

该式中，spd_j代表从DI得到的每个超像素块列向量，spx_j代表从X得到的每个超像素块的列向量，spy_j代表从Y得到的每个超像素块的列向量；

(3)对超像素总样本M进行稀疏表示，学习全局字典D同时得到稀疏系数矩阵G；

(4)对稀疏系数矩阵G进行低秩表示，得到低秩表示矩阵Z，选出变化的超像素样本和不变化的超像素样本，具体方法为：

(4a)构建目标函数：

s.t.G＝GZ+E

其中Z为求得的低秩表示，E为稀疏噪声，min[·]为最小化，||·||_*表示核范数，即矩阵奇异值的和，λ为正则参数，||·||_2,1表示L2,1范数；

(4b)稀疏系数矩阵G的每一行对应着超像素总样本M的每一列，即每一个超像素块，对稀疏系数矩阵G进行低秩表示后得到的Z为分块矩阵，根据Z的分块结构，L2范数大的为不变化的超像素样本，L2范数小的为变化的超像素样本；

其中，挑选不变化样本的规则为：

其中，Z_i表示低秩表示矩阵的第i行，α为挑选不变化的超像素样本的参数；这个规则表示在总样本中挑选L2范数偏大的行序号i，作为不变化的样本对应的列序号，再从超像素总样本M中挑选列序号i对应的超像素样本列，组成不变化样本M_un；挑选的原则是低秩表示矩阵的行的L2范数若大于等于除掉该行之外，其他行的二范数的均值，则认为是不变化的样本对应的行；

挑选变化样本的规则为：

其中，β为挑选变化的超像素样本的参数，Zs表示Z去掉不变化样本对应的行得到的结果，Zs_i表示Zs的第i行，Ns表示去掉不变化超像素之后的超像素总个数；这个规则表示在去掉不变化超像素的样本中挑选L2范数偏小的行的序号i，作为变化样本的对应列序号，再从M中挑选对应的超像素样本列作为变化样本M_c；

(5)分别学习不变化样本M_un和变化样本M_c对应的字典D_un和D_c；

(6)分别计算每个像素点在字典D_un和D_c下的重构误差Re_un和Re_c，比较大小，得到变化区域，具体方法为：

(6a)每个像素点的判别加入了邻域信息，即取以像素点为中心的邻域窗口为w×w的像素块，将DI、X、Y的像素分别列向量化，各列向量按像素序号排列得到T；

(6b)计算重构误差，需要字典D_un和D_c对应的稀疏表示系数g_un和g_c，为求g_un和g_c，构建目标函数：

其中t为T的各列向量，即用来表示各个像素块的列向量，δ为正则化参数，g_un ^T和g_c ^T分别代表g_un和g_c的转置矩阵，e为全1列向量；

(6c)计算像素块在字典D_un和D_c及对应的稀疏表示系数g_un和g_c下的重构误差Re_un和Re_c，若Re_un>Re_c，则判定该像素点不变化，反之变化，将变化位置的像素值设为1，不变化位置的像素值设为0，得到变化图二值图CM。

特别地，所述步骤(6b)中，目标函数的解为：

其中，D_UN＝D_un-t，D_C＝D_c-t，I为单位矩阵，I的尺寸为d×d，其中d为字典原子数。

特别地，所述步骤(6c)中，求重构误差的公式为：Re_un＝||t-D_ung_un||₂，Re_c＝||t-D_cg_c||₂。

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

第一，本发明将基于对象和基于像素的检测相结合，基于对象学习全局字典，降低了噪声的影响，基于像素进行变化检测，提升了变化检测的准确率。

第二，通常无监督方法通过自动阈值的选择得到变化检测的效果，而阈值选择的不同对结果的影响非常大，本发明通过低秩表示，可以无监督的挑选变化样本和不变化的样本，学习到的判别字典包括变化字典和不变化字典，省去了自动阈值的选择过程，使得本发明比现有技术更加鲁棒。

第三，本发明在判断像素点是否变化时加入邻域信息，减少了噪声对结果的影响，使得检测结果更加准确。

本发明的技术过程和效果可结合以下附图详细说明。

附图说明

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

图2为用本发明实现变化检测和其他检测方法的结果对比图。

具体实施方式

下面参照图1，结合实施例对本发明的实现过程进行详细说明。

本发明实施例中，首先输入两幅在同一地点不同时间拍摄的遥感图像；本发明实施例中，采用landsat5卫星于2009年4月6日在河北省市拍摄的TM图像和landsat7卫星于2003年4月30日在对应位置拍摄的ETM图像，其中TM图像大小为400×400×6，ETM图像大小为400×400×6，空间分辨率为30m。具体实现步骤如下：

步骤1.两幅不同时相的遥感图像：时相一X、时相二Y；

步骤2.计算差分图DI，用ERS(Entropy Rate Superpixel)方法对差分图DI进行超像素分割，得到超像素N个，进而得到超像素总样本M；超像素分割同时也可以采用SILC(Simple Linear Iterative Cluster)算法本实施例中，设置N为1200；

(2a)对DI、X、Y分别计算每一个超像素块各个波段上的平均像素值，公式为：

其中，i代表波段数，j代表超像素块序列，

代表每一个超像素块各个波段上的平均像素值，R_j代表第j个超像素块，N_j表示每个超像素块中像素的总个数，x代表像素点，X_i(x)代表第i个波段上遥感图像X在x像素点的像素值；

(2b)对DI、X、Y，按照超像素分割的结果，对每一个超像素块各个波段上的像素值平均，按波段顺序排列成列向量，三个列向量合成一个列向量sp_j后，各个超像素列向量按超像素序号排列得到超像素总样本M；其中，sp_j表示如下：

该式中，spd_j代表从DI得到的每个超像素块列向量，spx_j代表从X得到的每个超像素块的列向量，spy_j代表从Y得到的每个超像素块的列向量；本实施例中，超像素总样本为M_18×N；

步骤3.利用K-SVD算法学习全局字典D和稀疏系数矩阵G；字典学习也可以采用OMP算法，BP算法，PCA算法等；

本发明实施例中，设置最大迭代次数为200次，字典原子数设为100，计算全局字典D和稀疏系数矩阵G；

步骤4.对稀疏稀疏矩阵G进行低秩表示，得到低秩表示矩阵Z，由于在总样本中，变化的样本少，不变化的样本多，根据Z的分块结构，可以通过每一行的L2范数进行挑选，其中L2范数大的为不变化的样本，L2范数小的为变化的样本，因此可以根据对地址表示矩阵Z的L2范数挑选变化超像素样本和不变化超像素样本；

(4a)构建目标函数：

s.t.G＝GZ+E

其中Z为求得的低秩表示，E为稀疏噪声，min[·]为最小化，||·||_*表示核范数，即矩阵奇异值的和，λ为正则参数，||·||_2,1表示L2,1范数；

(4b)稀疏系数矩阵G的每一行对应着超像素总样本M的每一列，即每一个超像素块，对稀疏系数矩阵进行低秩表示后得到的Z为分块矩阵，根据Z的分块结构，L2范数大的为不变化的超像素样本，L2范数小的为变化的超像素样本；

其中，挑选不变化样本的规则为：

其中，Z_i表示低秩表示矩阵的第i行，α为挑选不变化的超像素样本的参数；这个规则表示在总样本中挑选L2范数偏大的行序号i，作为不变化的样本对应的列序号，再从超像素总样本M中挑选列序号i对应的超像素样本列，组成不变化样本M_un；挑选的原则是低秩表示矩阵的行的L2范数若大于等于除掉该行之外，其他行的二范数的均值，则认为是不变化的样本对应的行；

挑选变化样本的规则为：

其中，β为挑选变化的超像素样本的参数，Zs表示Z去掉不变化样本对应的行得到的结果，Zs_i表示Zs的第i行，Ns表示去掉不变化超像素之后的超像素总个数；这个规则表示在去掉不变化超像素的样本中挑选L2范数偏小的行的序号i，作为变化样本的对应列序号，再从M中挑选对应的超像素样本列作为变化样本M_c；

本实施例中，设置α＝1.5，β＝0.5，迭代次数为1000，λ为10000；

步骤5.用K-SVD算法分别学习不变化样本M_un和变化样本M_c对应的字典D_un和D_c；字典学习也可以采用OMP算法，BP算法，PCA算法等；

本实施例中，设置最大迭代次数为30次，字典原子个数d为12；

步骤6.分别计算每个像素点在字典D_un和D_c下的重构误差Re_un和Re_c，比较重构误差的大小，得到变化二值图CM；

(6a)每个像素点的判别加入了邻域信息，即取以像素点为中心的邻域窗口为w×w的像素块，将DI、X、Y的像素分别列向量化，各列向量按像素序号排列得到T；在本实施例中，设置邻域窗口为3×3大小，即w＝3；

(6b)计算重构误差，需要字典D_un和D_c对应的稀疏表示系数g_un和g_c，为求g_un和g_c，构建目标函数：

其中t为T的各列向量，即用来表示各个像素块的列向量，δ为正则化参数，在本实施例中，设置δ为500；g_un ^T和g_c ^T分别代表g_un和g_c的转置矩阵，e为全1列向量，求解此目标函数，可得g_un和g_c：

其中，D_UN＝D_un-t，D_C＝D_c-t，I为单位矩阵，I的尺寸为d×d，其中d为字典原子数，本实施例中，d为12；

(6c)计算像素块在字典D_un和D_c及对应的稀疏表示系数g_un和g_c下的重构误差Re_un和Re_c，其中，Re_un＝||t-D_ung_un||₂，Re_c＝||t-D_cg_c||₂，若Re_un>Re_c，则判定该像素点不变化，反之变化，将变化位置的像素值设为1，不变化位置的像素值设为0，得到变化图二值图CM。

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

1.仿真环境：

MATLAB R2017a,Intel(R)Core(TM)i7-8700CPU 3.20GHz,Window 10。

2.仿真内容：

仿真1，用本发明对landsat5卫星于2009年4月6日在河北省保定市拍摄的TM图像和landsat7卫星于2003年4月30日在对应位置拍摄的ETM图像进行变化检测，结果如图2所示，其中：

图2(a)是landsat7ETM图像，大小为400×400×6，

图2(b)是landsat5TM图像，大小为400×400×6，

图2(c)是一幅参考二值变化图，大小为400×400，

图2(d)为采用本发明对图2(a)和图2(b)进行变化检测后获得的变化检测二值图，大小为400×400。

由图2(d)和图2(c)的对比可见，本发明有检测变化检测区域的能力，而且准确率较高，错检漏检少。

仿真2，为了证明本发明的效果，分别用本发明的方法和现有技术差分(减法)阈值法、差分(除法)阈值法、IRMAD变换法对图2(a)和图2(b)的待融合图像进行融合，并对融合结果进行客观指标评价，评价指标如下：

1)假阳性FP(False Positive)，表示将不变化的像素点错误的检测成变化的像素点的总像素个数。FP越小，检测结果越好。

2)假阴性FN(False Negative)，表示将变化的像素点错误的检测成不变化的像素点的总像素个数。FN越小，检测结果越好。

3)准确率OA(OverallAccuracy)，表示检测的准确率，即正确检测的像素总个数/总像素个数，结果范围在区间[0,1]。OA越大，检测结果越好。

4)Kappa系数KC(Kappa Coefficient)，表示检测的精度，即像素点被正确分类的精度。Kappa系数越大，检测结果越好。

5)AUC(Area Under Curve)，表示在ROC(Receiver Operating Characteristic)曲线下与坐标轴围成的面积，结果范围在区间[0,1]。ROC曲线为受试者工作特征曲线，又称为感受性曲线，是以TP为纵坐标，FP为横坐标绘制的曲线，其中，TP表示将变化的像素正确分类到变化的总像素数。AUC的值越大，检测结果越好。

根据上述评价指标对本发明与现有技术的变化检测结果从客观评价指标上进行评价，结果如表1。

表1.各类方法变化检测结果的客观评价

Claims (3)

1.一种遥感图像变化检测方法，该方法对两张同一地点不同时相的遥感图像差分图进行超像素分割，基于分割对象学习一个全局字典和稀疏表示；对稀疏系数矩阵进行低秩表示挑选变化的样本及不变化的样本；分别用变化的和不变化的样本学习变化字典和不变化字典；计算每个像素的重构误差，得到变化区域二值图，具体步骤包括如下：

(1)输入两幅不同时相的遥感图像：时相一X和时相二Y；

(2)计算差分图DI，DI＝X-Y,对差分图DI进行超像素分割，得到N个超像素块，从而得到超像素总样本M，具体方法为：

(2a)对DI、X、Y分别计算每一个超像素块各个波段上的平均像素值，公式为：

其中，i代表波段数，j代表超像素块序列，

分别代表在DI、X、Y上每一个超像素块各个波段上的平均像素值，R_j代表第j个超像素块，N_j表示每个超像素块中像素的总个数，x代表像素点，DI_i(x)、X_i(x)、Y_i(x)分别代表遥感图像DI、X、Y第i个波段上在x像素点的像素值；

(2b)对DI、X、Y，按照超像素分割的结果，对每一个超像素块各个波段上的像素值平均，按波段顺序排列成列向量，三个列向量合成一个列向量sp_j后，各个超像素列向量按超像素序号排列得到超像素总样本M，其中，sp_j表示如下：

该式中，spd_j代表从DI得到的每个超像素块列向量，spx_j代表从X得到的每个超像素块的列向量，spy_j代表从Y得到的每个超像素块的列向量；

(3)对超像素总样本M进行稀疏表示，学习全局字典D同时得到稀疏系数矩阵G；

(4)对稀疏系数矩阵G进行低秩表示，得到低秩表示矩阵Z，选出变化的超像素样本和不变化的超像素样本，具体方法为：

(4a)构建目标函数：

s.t.G＝GZ+E

其中Z为求得的低秩表示，E为稀疏噪声，min为最小化，||·||_*表示核范数，即矩阵奇异值的和，λ为正则参数，||·||_2,1表示L2,1范数；

(4b)稀疏系数矩阵G的每一行对应着超像素总样本M的每一列，即每一个超像素块，对稀疏系数矩阵G进行低秩表示后得到的Z为分块矩阵，根据Z的分块结构，L2范数大的为不变化的超像素样本，L2范数小的为变化的超像素样本；

其中，挑选不变化样本的规则为：

其中，Z_i表示低秩表示矩阵的第i行，α为挑选不变化的超像素样本的参数；这个规则表示在总样本中挑选L2范数偏大的行序号i，作为不变化的样本对应的列序号，再从超像素总样本M中挑选列序号i对应的超像素样本列，组成不变化样本M_un；挑选的原则是低秩表示矩阵的行的L2范数若大于等于除掉该行之外，其他行的二范数的均值，则认为是不变化的样本对应的行；

挑选变化样本的规则为：

其中，β为挑选变化的超像素样本的参数，Zs表示Z去掉不变化样本对应的行得到的结果，Zs_i表示Zs的第i行，Ns表示去掉不变化超像素之后的超像素总个数；这个规则表示在去掉不变化超像素的样本中挑选L2范数偏小的行的序号i，作为变化样本的对应列序号，再从M中挑选对应的超像素样本列作为变化样本M_c；

(5)分别学习不变化样本M_un和变化样本M_c对应的字典D_un和D_c；

(6)分别计算每个像素点在字典D_un和D_c下的重构误差Re_un和Re_c，比较大小，得到变化区域，具体方法为：

(6a)每个像素点的判别加入了邻域信息，即取以像素点为中心的邻域窗口为w×w的像素块，将DI、X、Y的像素分别列向量化，各列向量按像素序号排列得到T；

(6b)计算重构误差，需要字典D_un和D_c对应的稀疏表示系数g_un和g_c，为求g_un和g_c，构建目标函数：

其中t为T的各列向量，即用来表示各个像素块的列向量，δ为正则化参数，g_un ^T和g_c ^T分别代表g_un和g_c的转置矩阵，e为全1列向量；

(6c)计算像素块在字典D_un和D_c及对应的稀疏表示系数g_un和g_c下的重构误差Re_un和Re_c，若Re_un>Re_c，则判定该像素点不变化，反之变化，将变化位置的像素值设为1，不变化位置的像素值设为0，得到变化图二值图CM。

2.根据权利要求1所述的遥感图像变化检测方法，其特征在于：所述步骤(6b)中，目标函数的解为：

其中，D_UN＝D_un-t，D_C＝D_c-t，I为单位矩阵，I的尺寸为d×d，其中d为字典原子数。

3.根据权利要求1所述的遥感图像变化检测的方法，其特征在于：所述步骤(6c)中，求重构误差的公式为：Re_un＝||t-D_ung_un||₂，Re_c＝||t-D_cg_c||₂。

Images