基于深度共享表示的SAR图像变化检测方法与流程

本发明属于图像处理领域，更进一步涉及遥感图像变化检测技术领域中的一种基于深度共享表示的合成孔径雷达sar(syntheticapertureradar)图像变化检测方法。本发明可用于农业调查、自然灾害检测、森林资源监测中检测两幅不同时相的sar图像变化区域。

背景技术：

变化检测是指通过观测某一物体或现象在不同时间的状态得出其变化信息的技术。sar图像变化检测是分析同一地区不同时段的sar图像，检测出该地区的变化信息。作为对地观测卫星的关键技术，合成孔径雷达sar图像变化检测已经被运用在包括农业调查、自然灾害检测、森林资源监测等诸多领域。变化检测指利用多时相获取的覆盖统一地表区域的遥感图像以及其他辅助数据来确定和分析地表变化。传统变化检测方法由以下经典三步流程范式：1)输入两张经过预处理后的合成孔径雷达sar图像；2)使用差异算子或其他方法获得差异图；3)分析差异图。早期人们使用差值算子求取差异图，对于传统的变化检测方法，性能好坏依赖于差异图的形成，而差异图的形成又有噪声的干扰，所以得不到很高的精度。由于采集图像有各种各样的噪声，后来人们相继提出对数比算子、均值比算子等。

广东工业大学在其申请的专利文献“一种sar图像变化检测的方法与装置”(专利申请号cn201710597812.2，公开号cn107451992a)中公开了一种基于一种sar图像变化检测的方法。该方法首先将两幅sar图像分为参考图像和测试图像；再将这两幅图像进行重叠分块处理，得到参考图像块和测试图像块；通过构建局部字典，并利用正交匹配追踪法，求解出测试图像块的稀疏估计值；然后依据测试图像块稀疏估计值和参考图像块稀疏值，构建比值图像，并对其进行中值滤波处理，得到sar图像变化检测结果。该方法存在的不足之处是，采用基于字典学习，正交匹配追踪求解稀疏值的线性分类方法，对于非线性分布的像素点无法很好地分类，影响分类结果。

西安电子科技大学在其申请的专利文献“基于最大边缘度量学习的sar图像变化检测方法”(专利申请号cn201810736484.4，公开号cn108846381a)中公开了一种基于最大边缘度量学习的sar图像变化检测方法。该方法先输入变化前后遥感图像采用链式结构构建全部样本，构建包含全部边界的训练样本，然后利用训练样本构建正负约束对作为输入，建立结构化的支撑矢量机模型获得映射矩阵，其次利用分解后的映射矩阵将全部样本映射到特征空间，在特征空间对全部样本进行sar图像变化检测分类。该方法存在的不足之处是，时相和差异图样本采取链式拼接，不能很好的提取时相图和差异图像素之间的信息，反而增加训练时间，鲁棒性较差，导致变化检测结果图中的噪声点较多。

技术实现要素：

本发明的目的在于针对上述现有技术的不足，提出了一种基于深度共享表示的sar图像变化检测方法，用于解决线性分类无法对非线性分布像素点进行准确分类，不能很好提取时相图和差异图像素之间的信息的问题。

实现本发明目的的思路是，构建三个深度神经网络网络来提取时相图和差异图的特征，利用深度神经网络对训练样本映射，提取后的特征能够对非线性分布的像素进行很好的分类；对深度神经网络映射结果进行广义相关分析，学习共享表示矩阵，利用共享矩阵能实现对两张时相图与差异图之间的信息融合；然后用训练好的神经网络对训练样本和测试样本进行映射，最后用支持向量机对映射结果进行分类，得到sar图像变化检测结果。

本发明的具体实施步骤如下：

(1)构建两个时相图像的差异图：

(1a)输入两个同一地区、不同时相的合成孔径雷达sar图像；

(1b)利用对数比值算子公式，计算第一个时相图像中每个像素与第二个时相图像对应位置像素的对数比值，将所有像素的对数比值按时相图对应位置组成差异图；

(2)生成训练样本和样本标签：

(2a)对差异图进行归一化处理；

(2b)分别以处理后的差异图、第一个时相图、第二个时相图、标签图像中的每个像素点为中心，按垂直方向选取每个像素点周围27×27大小的正方形的图像块；

(2c)将差异图中的每个图像块按列拉成列向量作为一个样本，将差异图中所有图像块形成的样本组成差异图测试样本集；

(2d)将第一个时相图中每个图像块按列拉成列向量作为一个样本，将第一时相图中所有图像块形成的样本组成第一时相图的测试样本集；

(2e)将第二个时相图中每个图像块按列拉成列向量作为一个样本，将第二时相图所有图像块形成的样本组成第二时相图的测试样本集；

(2f)将标签图像每个图像块按列拉成列向量作为一个标签，将第二时相图所有图像块形成的样本组成标签图像的测试样本标签；

(2g)从第二个时相图的测试样本集中选取30％的样本组成第一时相训练样本集；

(2h)从第一个时相图的测试样本集中选取30％的样本组成第二时相训练样本集；

(2i)从差异图的测试样本集中选取30％的样本组成差异图训练样本集；

(2j)从测试样本标签中选取30％的标签组成训练样本标签；

(3)构建深度神经网络：

(3a)搭建三个深度为五层的神经网络，每个神经网络的结构依次为，输入层→第一层隐藏层→第二层隐藏层→第三层隐藏层→输出层；

(3b)设置神经网络的参数；

(4)生成映射矩阵：

将第一时相训练样本集、第二时相训练样本集、差异图训练样本集分别输入到第一神经网络、第二神经网络、第三神经网络中，分别输出第一时相、第二时相、差异图的映射矩阵，每个矩阵的行均为30，列均与训练样本的总数相等；

(5)更新网络：

(5a)利用相关映射矩阵方法，生成合相关映射矩阵和经验协方差矩阵；

(5b)对合相关映射矩阵进行奇异值分解，得到合相关映射矩阵的特征值和对应特征值的特征向量，将其特征值由大到小排列，提取前30个特征值对应的特征向量组成共享矩阵；

(5c)对经验协方差矩阵、映射矩阵和共享矩阵进行链式乘处理，得到三个投影矩阵；

(5d)计算第一时相，第二时相、差异图的损失函数梯度值；

(5e)利用反向传播算法分别对三个神经网络进行迭代更新，直到三个神经网络的损失函数之和收敛，得到训练好的三个神经网络；

(6)生成标签向量：

(6a)将第一时相训练样本集、第二时相训练样本集、差异图训练样本集分别输入到第一神经网络、第二神经网络、第三神经网络中，分别输出三个训练样本映射矩阵，将三个训练样本映射矩阵按行拼接，得到一个合训练样本矩阵；

(6b)将第一时相测试样本集、第二时相测试样本集、差异图测试样本集分别输入到第一神经网络、第二神经网络、第三神经网络中，分别输出三个测试样本映射矩阵，将三个测试样本映射矩阵按行拼接，得到一个合测试样本矩阵；

(6c)输入合训练样本和训练样本标签，训练支持向量机，利用训练好的支持向量机对合测试样本矩阵进行分类，生成标签向量，该标签向量的每个元素为总测试样本矩阵每一行的标签；

(7)获得sar图像变化检测结果：

(7a)构建一个行的总数等于第一时相图的宽，列的总数等于第一时相图的高的标签矩阵；

(7b)从标签向量中依次选取与标签矩阵列的总数相等的元素，填入到标签矩阵的列中；

(7c)将标签矩阵中的每个元素乘以255后按图像形式输出，得到输入时相图所在地区是否发生变化的sar图像变化的检测结果。

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

第一，由于本发明构建了三个深度神经网络用于提取图像特征，克服了现有技术中采用基于字典学习，正交匹配追踪求解稀疏值的线性分类方法，对于非线性分布的像素点无法很好地分类，影响分类结果的缺点，使得本发明具有能够将非线性分布的像素点进行很好的分类，提高了检测变化精度的优点。

第二，由于本发明学习第一时相、第二时相、差异图训练样本的共享矩阵，克服了现有技术中采用链式拼接时相和差异图样本的方法，不能很好的提取时相图和差异图像素之间的信息，反而增加训练时间，鲁棒性较差，导致变化检测结果图中的噪声点较多的缺点，使得本发明能够很好的提取第一时相图，第二时相图和差异图之间的相关信息，提高训练效率，减少样本图像中噪声点的影响，鲁棒性强。

附图说明：

图1为本发明的流程图；

图2为本发明仿真图。

具体实施方式：

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

参照图1，对本发明的具体步骤做进一步的描述。

步骤1，构建两个时相图像的差异图。

输入两个同一地区、不同时相的合成孔径雷达sar图像。

利用对数比值算子公式，计算第一个时相图像中每个像素与第二个时相图像对应位置像素的对数比值，将所有像素的对数比值按时相图对应位置组成差异图。

所述的对数比值算子公式如下：

其中，idl表示第一个时相图像中第l个像素与第二个时相图像对应位置像素的对数比值，|·|表示取绝对值操作，lg表示以10为底的对数操作，i2m表示第二个时相图中第m个像素，i1n表示第一个时相图中第n个像素，l、m、n的位置对应相同。

步骤2，生成训练样本和样本标签。

对差异图进行归一化处理。

分别以处理后的差异图、第一个时相图、第二个时相图、标签图像中的每个像素点为中心，按垂直方向选取每个像素点周围27×27大小的正方形的图像块。

所述标签图像是指，根据输入的第一时相图对应地区的光学图像和第二时相图对应地区的光学图像人工标注的该地区已经变化和没有非变化部分的参考图，或者去实地考察标注的该地区已经变化和没有非变化部分的参考图。

将差异图中的每个图像块按列拉成列向量作为一个样本，将差异图中所有图像块形成的样本组成差异图测试样本集。

将第一个时相图中每个图像块按列拉成列向量作为一个样本，将第一时相图中所有图像块形成的样本组成第一时相图的测试样本集。

将第二个时相图中每个图像块按列拉成列向量作为一个样本，将第二时相图所有图像块形成的样本组成第二时相图的测试样本集。

将标签图像每个图像块按列拉成列向量作为一个标签，将第二时相图所有图像块形成的样本组成标签图像的测试样本标签。

从第二个时相图的测试样本集中选取30％的样本组成第一时相训练样本集。

从第一个时相图的测试样本集中选取30％的样本组成第二时相训练样本集。

从差异图的测试样本集中选取30％的样本组成差异图训练样本集。

从测试样本标签中选取30％的标签组成训练样本标签。

步骤3，构建深度神经网络。

搭建三个深度为五层的神经网络，每个神经网络的结构依次为，输入层→第一层隐藏层→第二层隐藏层→第三层隐藏层→输出层。

设置神经网络的参数。

所述设置神经网络的参数如下：

将多层神经网络的输入层的神经元个数设置为729；

将第一个、第二个、第三个隐藏层均设置25个神经元，输出层设置30个神经元；

将每层的激活函数设置为relu函数；

按照下式，设置每个神经网络的损失函数如下：

其中，g表示共享矩阵，u表示投影矩阵，x表示该神经网络的输出矩阵。

步骤4，生成映射矩阵。

将第一时相训练样本集、第二时相训练样本集、差异图训练样本集分别输入到第一神经网络、第二神经网络、第三神经网络中，分别输出第一时相、第二时相、差异图的映射矩阵，每个矩阵的行均为30，列均与训练样本的总数相等。

步骤5，更新网络。

利用相关映射矩阵方法，生成合相关映射矩阵和经验协方差矩阵。

所述相关映射矩阵方法的步骤如下：

第一步，用第一时相的映射矩阵x1、第二时相的映射矩阵x2、差异图的映射矩阵x3分别与其转置后的映射矩阵x1^t、x2^t、x3^t相乘，得到第一时相的经验协方差矩阵c11、第二时相的经验协方差矩阵c22、差异图经验协方差矩阵c33；

第二步，用第一时相的映射矩阵的转置矩阵x1^t与其对应的经验协方差矩阵的逆矩阵相乘，在乘上第一时相的映射矩阵x1，得到第一时相的相关映射矩阵p1；用第二时相的映射矩阵的转置矩阵x2^t与其对应的经验协方差矩阵的逆矩阵相乘，在乘上第二时相的映射矩阵x2，得到第二时相的相关映射矩阵p2；用差异图映射矩阵的转置矩阵x3^t与其对应的经验协方差矩阵的逆矩阵相乘，在乘上差异图映射矩阵x3，得到差异图的相关映射矩阵p3；

第三步，将第一时相的相关映射矩阵p1，第二时相的相关映射矩阵p2与差异图的相关映射矩阵p3相加，计算得到合相关映射矩阵。

对合相关映射矩阵进行奇异值分解，得到合相关映射矩阵的特征值和对应特征值的特征向量，将其特征值由大到小排列，提取前30个特征值对应的特征向量组成共享矩阵。

对经验协方差矩阵、映射矩阵和共享矩阵进行链式乘处理，得到三个投影矩阵。

所述对经验协方差矩阵、映射矩阵和共享矩阵进行链式乘处理是指，将第一时相的经验协方差矩阵的逆矩阵与其映射矩阵x1相乘，在乘上共享矩阵g，得到第一时相的投影矩阵u1；将第二时相的经验协方差矩阵的逆矩阵与其映射矩阵x2相乘，在乘上共享矩阵g，得到第二时相的投影矩阵u2；将差异图的经验协方差矩阵的逆矩阵与其映射矩阵x3相乘，在乘上共享矩阵g，得到差异图的投影矩阵u3。

计算第一时相，第二时相、差异图的损失函数梯度值。

所述第一时相，第二时相、差异图的损失函数梯度值是由下述公式计算得到的：

第一步，按照下式，计算第一时相的损失函数梯度值：

其中，p1表示第一时相的损失函数梯度值，u1表示第一时相的投影矩阵，g表示共享矩阵，t表示转置操作，x1表示第一时相的映射矩阵；

第二步，按照下式，计算第二时相的损失函数梯度值：

其中，p2表示第二时相的损失函数梯度值，u2表示第二时相的投影矩阵，x2表示第二时相的映射矩阵；

第三步，按照下式，计算差异图的损失函数梯度值：

其中，p3表示差异图的损失函数梯度值，u3表示差异图的投影矩阵，x3表示差异图的映射矩阵。

利用反向传播算法分别对三个神经网络进行迭代更新，直到三个神经网络的损失函数之和收敛，得到训练好的三个神经网络。

步骤6，生成标签向量。

将第一时相训练样本集、第二时相训练样本集、差异图训练样本集分别输入到第一神经网络、第二神经网络、第三神经网络中，分别输出三个训练样本映射矩阵，将三个训练样本映射矩阵按行拼接，得到一个合训练样本矩阵。

将第一时相测试样本集、第二时相测试样本集、差异图测试样本集分别输入到第一神经网络、第二神经网络、第三神经网络中，分别输出三个测试样本映射矩阵，将三个测试样本映射矩阵按行拼接，得到一个合测试样本矩阵。

输入合训练样本和训练样本标签，训练支持向量机，利用训练好的支持向量机对合测试样本矩阵进行分类，生成标签向量，该标签向量的每个元素为总测试样本矩阵每一行的标签。

步骤7，获得sar图像变化检测结果。

构建一个行的总数等于第一时相图的宽，列的总数等于第一时相图的高的标签矩阵。

从标签向量中依次选取与标签矩阵列的总数相等的元素，填入到标签矩阵的列中。

将标签矩阵中的每个元素乘以255后按图像形式输出，得到输入时相图所在地区是否发生变化的sar图像变化的检测结果。

下面结合仿真实验对本发明的效果做进一步的说明：

1.仿真实验条件：

本发明的仿真实验的硬件平台为：处理器为inteli9-9700kcpu，主频为3.6ghz，内存16g。

本发明的仿真实验的软件平台为：ubuntu16.04操作系统和python3.6。

2.仿真内容及其结果分析：

本发明仿真实验是采用本发明和现有技术支持向量机svm分类对输入的两组分别对应sendai和velcea两个地区的图像进行变化检测，获得变化检测结果图。

在仿真实验中，采用的现有技术支持向量机svm分类方法是指，melgani等人在“classificationofhyperspectralremotesensingimageswithsupportvectormachines,ieeetrans.geosci.remotesens.,vol.42,no.8,pp.1778–1790,aug.2004”中提出的高光谱图像分类方法，简称支持向量机svm分类方法。

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

图2(a)为2011年日本大地震发生之前sendai地区sar图，成像于2010年10月20号，图2(b)为日本大地震发生之后sendai地区sar图，图2(c)为标记日本大地震发生之前后sendai地区标签图，图2(d)为由图2(a)和图2(b)得到的差异图。四个图像大小都为590×687，图像格式为jpg格式。图2(e)为本发明和现有技术支持向量机对sendai地区变化检测的结果图。

从图2(d)中可以看出，在右下角有很大一片白色干扰噪声，从图2(e)图中可以看出，变化检测的结果图很好的抑制了噪声，只有在变化区域周围有一部分错检，同时从结果图也可以看出对于变化和非变化边界区域检测很好，在对噪声鲁棒的同时很好的保留了边界的细节，达到了很好的性能。