CN109948708A - 基于迭代隐式正则的多时多光谱图像特征级信息融合方法 - Google Patents

Abstract

本发明公开了一种基于迭代隐式正则的多时多光谱图像特征级信息融合方法，包括：(1)输入两幅多时多光谱图像；(2)超像素提取；(3)对超像素进行无监督聚类，产生分组信息；(4)特征提取；(5)对提取的特征进行操作生成伪标签；(6)根据样本数据、伪标签数据和组信息构成训练样本集；(7)选择并初始化分类器；(8)判断是否满足终止条件，如果满足则执行(10)，否则执行(9)；(9)用基于迭代隐式正则的方法训练分类器，分类器训练结束后返回步骤(8)；(10)用分类器对所有特征进行分类得到最终的融合结果图像。本发明不仅可以通过融合辨别出两幅高光谱图像的变化信息，还可以区分出不同的变化类型。

Description

基于迭代隐式正则的多时多光谱图像特征级信息融合方法

技术领域

本发明属于图像处理技术领域，涉及多时多光谱图像特征融合方法，具体是一种基于迭代隐式正则的多时多光谱图像特征级信息融合方法，该方法可以得到更好的融合结果图像用于消除其中的互斥的、冗余的信息，减少其中的模糊的不确定信息，最终得到这个场景中的有用的、可靠的图像信息，可以通过融合结果图像进一步得到同一地区多时多光谱图像中的变化部分，可以应用于洪水灾害前后的受灾情况评估、医疗诊断、冰川变化检测、地震灾害评估、土地资源检测、城市规划、环境检测、视频监督等。

背景技术

多时多光谱图像特征级信息融合是一门针对同一场景下多个时刻的多个波段信息或者将不同类型的图像信息，消除其中的互斥的、冗余的信息，减少其中的模糊的不确定信息，最终得到这个场景中的有用的、可靠的图像信息的技术。通过多时图像特征级信息融合技术进一步分析和理解之后可以为后续的图像中的目标识别、检测、跟踪、变化检测等相关研究提供更加可靠的数据来源。随着遥感技术的不断发展，多光谱图像已经成为图像融合技术的重要来源，多光谱图像具有多个波段信息且具有较高的分辨率，因此多时多光谱图像特征集信息融合技术的应用需求不断提高。

随着国内外多时多光谱图像特征级信息融合技术的不断发展，相关的研究人员相继提出了许多有效的解决方法，其中应用最广泛的是利用变化矢量分析的方法先对图像进行特征提取，然后利用无监督方法对特征进行聚类最终产生结果图像。无监督方法中最常用的是阈值法和聚类法，其中最广泛应用的是模糊C均值聚类方法。然而这种方法的融合结果效果并不理想，分析原因是因为一方面图像中存在噪声，另一方面，模糊C均值聚类方法本身处理能力有限，因为没有利用图像中的空间信息，例如邻域像素的信息等。

随着机器学习和深度学习的不断发展，一些研究人员提出了利用监督学习的方法来解决此类问题，其主要思路是：首先对图像进行特征提取，然后利用无监督方法进行预分类得到一个伪标签集，通过选取特征和这个伪标签集构成训练样本集，然后将训练样本送入分类器进行训练，最后利用训练的分类器对所有特征进行分类，得到最后的融合结果。这种方法可以明显改善图像融合的结果，并且可以一定程度上克服图像中的噪声影响的问题。然而对于一些高分辨率的图像，该方法得到的结果依然受到噪声的严重影响，产生了较差的融合结果。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于迭代隐式正则的多时多光谱图像特征级信息融合方法，以提高多时多光谱图像特征集信息融合的结果。本发明既可以较好地融合得到两幅图像的变化区域，又可以进一步辨别出变化区域的变化类型，具有较广泛的使用范围和较高的检测精度。

本发明的技术方案是：将多时多光谱图像融合问题看作一个分类问题，分类的标签来自无监督预融合结果，分类的样本数据来自两幅图像的特征提取。分类器的训练是通过在原始分类器的目标函数中添加一个隐式正则项，其实质是通过该正则项来不断地为分类器选取可靠样本，从而得到一个性能良好的分类器模型。然后通过分类器对特征进行分类，得到最后的融合图像，即表明两幅图像变化信息的图像。

其实现步骤包括如下：

(1)输入多时多光谱图像X₁和X₂：输入在同一地区，不同时刻获取的两幅多光谱图像，分别标记为X₁和X₂。

(2)两幅原始图像进行超像素提取操作：采用SLIC方法对原始图像进行超像素提取操作；

(3)对提取的超像素进行无监督聚类产生组信息：用模糊C均值聚类(FCM)方法，对超像素进行聚类；

(4)根据原始两幅图像进行特征提取；

(5)对提取的特征进行无监督聚类(FCM)操作生成伪标签；

(6)根据样本数据、伪标签数据和组信息构成训练样本集；对提取的特征进行采样生成训练样本数据：以图像中的像素及它的ρ×ρ邻域特征组成训练样本和测试样本数据；

(7)选择一个合适的分类器并初始化分类器的模型参数；

(8)判断是否满足终止条件，如果满足则执行(10)，否则执行(9)；

(9)用基于迭代隐式正则的方法训练分类器，训练完后执行步骤(8)；

(10)用分类器对所有特征进行分类得到最终的融合结果图像。

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

第一，本发明通过使用一个隐式正则项，可以自主地为分类器选取可靠样本，而不需要人工去选择，克服了现有技术中样本难以获取可靠的训练样本的问题。

第二，由于本发明采用了有监督分类方法对差异特征进行分类，克服了现有技术在有噪声存在的情况下的精度低的问题，可以有效地解决图像中存在噪声时的图像特征集融合问题，提高了噪声存在情况下的利用融合结果的效果。

第三，本发明通过使用多分类器对不同地物变化前后的图像进行融合可以得到两幅图像的变化信息，克服了现有技术中无法得到变不同的化类型的缺点，使得本发明可以获得具有不同地物信息的融合结果图像，可以显示出不同地物的不同的变化类型。

附图说明

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

图2为本发明仿真采用的多光谱遥感图像；

图3为本发明与现有的四种方法对数据集进行特征级信息融和的结果。

具体实施方式

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

步骤1，输入多时图像X₁和X₂：输入在同一地区，不同时刻获取的融合前后的两幅图像，分别标记为X₁和X₂。

步骤2，对两幅原始多光谱图像进行超像素提取操作：采用SLIC方法对原始图像进行超像素提取操作；

2a)随机初始化种子点，即聚类中心：按照设定的超像素个数，在两幅原始的多光谱图像内均匀地分配种子点，通常的方法是随机选取图像内的一定数量的像素点作为初始化种子点。

2b)计算每一个种子的n×n邻域(一般n＝3)内所有像素点的梯度值，并用该种子点移动到该邻域内梯度最小点替代该种子点的像素值作为新的种子点。

2c)在每个种子点周围的邻域内为每个像素点分配类标签(即属于哪个聚类中心)。

2d)对每个像素点进行颜色距离和空间距离度量；

2e)判断是否达到预设的最大迭代次数(一般设定为10，因为实践发现一般在10次迭代内绝大部分图片都可以你得到比较理想的分割效果)，如果达到最大迭代次数，则进行2f)，否则返回2b)；

2f)进行区域连通操作。新建一张标记表，表内元素均为-1，按照“Z”型走向(从左到右，从上到下顺序)将不连续的超像素、尺寸过小超像素重新分配给邻近的超像素，遍历过的像素点分配给相应的标签，直到所有的像素点都遍历完毕为止。

步骤3，对提取的超像素进行无监督聚类，产生组信息：用模糊C均值聚类(FCM)方法，对超像素进行聚类；

3a)随机初始化聚类中心C。通常的做法是从所有的数据点钟任取C个点作为初始聚类中心。

3b)根据聚类中心C计算当前所有样本数据对每个聚类中心的属度u_ij并组成隶属度矩阵U，其中隶属度u_ij的计算公式为：

其中c为聚类类别数目，m为模糊因子，一般取值为m＝2,d_ij表示第i个样本数据到第j个聚类中心的欧式距离，u_ij表示第i个样本数据相对于第j类聚类中心的隶属度。

3c)根据聚类中心C和隶属度矩阵U计算当前迭代次数的代价函数值J，并得到相邻两次迭代的目标函数值J的改变量σ，其中J的计算公式如下：

其中N表示需要聚类的样本数据的总个数。

3d)根据改变量σ和当前迭代次数判断是否满足终止条件(总代价函数值小于预设阈值或迭代次数达到预设最大迭代次数)，满足则输出隶属度矩阵U，否则执行3e)；

3e)采用梯度下降法计算新的聚类中心C，并将迭代次数加1；返回步骤3b)。其中第j个聚类中心c_j的计算方式如下：

步骤4，根据原始两幅图像进行特征提取；

4a)对X₁和X₂图像做相对差值运算得到特征矢量。记图像中的每个像素灰度矢量分别表示为X₁＝(g₁,g₂,...,g_n)^T，X₂＝(h₁,h₂,...,h_n)^T，其中n为波段数。变化矢量用公式可以表示为：

4b)：计算变化强度：

4c)：计算变化方向(cosθ₁,cosθ₂,...,cosθ_n)：

4d)：由特征强度和特征方向构成特征(|△V|，cosθ₁,cosθ₂,...,cosθ_n)。

步骤5，对提取的特征进行无监督聚类(FCM)操作生成预融合结果，并且将这个结果记为伪标签；

步骤6，根据样本数据、伪标签数据和组信息构成训练样本集；对提取的进行采样生成训练样本数据：以图像中的像素及它的ρ×ρ邻域特征组成训练样本数据；

步骤7，选择一个合适的分类器并初始化分类器的模型参数；

步骤8，判断是否满足终止条件，如果满足则执行(10)，否则执行(9)；

步骤9，用基于迭代隐式正则的方法训练分类器；

9a)随机初始化神经网络的参数W，b，正则项参数；

9b)根据神经网络分类器的模型参数计算每个训练样本x_i的代价值l_i，并根据隐式正则项和l_i、λ、μ计算每个训练样本x_i的样本权重v_i；

9b1)根据模型参数计算每个训练样本的损失值计算如下：

其中为第j组内位于第i个训练样本的标签信息的one-hot编码形式，为第j组内位于第i个训练样本经过神经网络的输出层的输出结果。

9b2)将所有样本损失值进行归一化处理并在各自的组内按照升序方式进行排列；

9b3)根据隐式正则项推导出样本权重的计算方式。带隐式正则的神经网络的目标函数为：

其中表示第j组内第i个样本的代价值，第二项为隐式正则，计算方式为：

求新的目标函数E(W,v)并求对v的偏导数令该偏导数为零，即得到样本的权重计算公式为：

9b4)根据该权重计算方法计算所有训练样本的权重。

9c)根据带权重的样本训练神经网络，神经网络的训练过程如下：

9c1)初始化神经网络的模型参数W,b；

9c2)神经网络进行前向传播，根据每个训练样本计算网络的每层结点输出并且计算当前迭代次数的总的代价值l，其中l的计算方式为其中N为训练样本总数；

9c3)神经网络进行反向传播计算模型参数的梯度；

其中δ_L＝h_L-y_i＝O_i-y_i,δ_l＝W_lδ_l+1⊙f'(W_lh_l-1+b_l),l＝1,2,...,L-1,⊙表示点积运算。

9c4)根据梯度和样本权重V更新模型参数，对于每个训练样本x_i神经网络模型的参数的更新公式如下：

9c5)判断是否当前迭代次数是否达到最大迭代次数，若达到，则输出模型参数W,

b,否则，返回5c2)。

步骤10，用分类器对所有特征进行分类得到最终的融合结果图像。

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

1.仿真条件

本实例在Intel(R)Core(TM)i5-4210U CPU@1.70GHz 2.10GHZ Windows 10系统下，Matlab 2016a运行平台上，完成本发明与现有方法CVA，CVA+NN，CVA+SVM和CVA+Softmax的仿真实验。

2.仿真实验内容

用本发明和现有的C2VA方法对多光谱图像数据集进行变化检测，其结果如图2和图3所示，其中：

图2(a)为待处理的第一幅图像X1；

图2(b)为待处理的第一幅图像X2；

图2(c)为特征级融合结果变化检测的真实图像X0；

图3(a)为CVA方法的变化检测结果图像；

图3(b)为CVA+NN方法的变化检测结果图像；

图3(c)为CVA+SVM方法的变化检测结果图像；

图3(d)为CVA+Softmax方法变化检测结果图像；

图3(e)为本发明的融合变化检测结果图像；

从图3可以看出，与其他的对比方法相比，本发明的方法变化检测结果图像更加接近真实的结果图像：错分点明显减少，各类变化检测的检测准确率增加，各类变化的精度明显提高，证明了本方法可以得到较好的融合结果。

将本发明方法和其他四种对比方法的结果图，即图3中(a)(b)(c)(d)(e)与真实的变化检测结果图进行对比，统计每种方法的OE，PCC，KC和NMI值，其中：

OE表示所用方法分类错误的像素总数；

PCC表示所用方法的分类正确率，计算方式如下：

KC表示为综合评价指标，可以很好的评价融合结果对变化检测结果的好坏，其值越大表明结果越好，计算如下：

NMI表示标准化的互信息，可以衡量融合结果与参考图像的一致性程度，其值越大表明融合结果与参考图越相似，结果越好。NMI的计算方法如下：

其中t＝[1,1,...,1]_1×M，M为总的类别数，N为每幅图像中的像素总数，n_r表示在参考图X0中属于类别r的像素总数，n_p在结果图中属于类别p的像素总，n_rp表明在相同位置上属于参考图中r类，属于结果图中p类的像素总数。D为混淆矩阵：

D的统计方法为：

表1：混淆矩阵D的统计方法：

将本发明与CVA，CVA+NN，CVA+SVM和CVA+Softmax方法在灞河数据集上融合的结果的OE，PCC，KC和NMI值进行统计，如表2所示。其中OE的值越小越好，PCC，KC，NMI的值越大越好。

表2 不同方法对灞河数据集的融合变化检测结果

方法	OE	PCC	KC	NMI
					CVA	25172	0.79451	0.41246	0.34278
CVA+NN	21417	0.82517	0.55789	0.45907
					CVA+SVM	24576	0.79938	0.43312	0.35945
CVA+Softmax	79447	0.35145	0.21505	0.32272
					本发明	20317	0.83415	0.61241	0.51081

从表2可以看出，对于同一组高广谱图像数据集，CVA和CVA+Softmax的效果最差，其他的结果要比这两种方法的好很多，但是其误检数OE的值依然很大。而本发明方法具有最低的误检数OE和最高的正确率PCC，且具有最高的KC和NMI，比其他的4中方法的KC都要高出至少5个百分点，表明采用本发明方法对高广谱图像进行特征级融合可以取得比较理想的结果。

Claims (10)

1.一种基于迭代隐式正则的多时多光谱图像特征级信息融合方法，包括下列步骤：

(1)输入两幅原始图像X₁和X₂：输入在同一地区，不同时刻获取的两幅多光谱图像，标记为X₁和X₂；

(2)对这两幅原始图像进行超像素提取操作：采用简单线性迭代聚类法(SLIC方法)，对两幅原始图像分别进行超像素提取操作，得到这两幅原始图像的超像素；

(3)对提取的超像素进行无监督聚类产生组信息：用模糊C均值聚类(FCM)方法，对超像素进行聚类，产生组信息；

(4)对两幅原始图像进行特征提取，得到提取的特征；

(5)对提取的特征进行无监督聚类：利用FCM对步骤(4)提取的特征进行聚类产生伪标签；

(6)由步骤(4)提取的特征组成的样本数据、步骤(5)生成的伪标签和步骤(3)产生的组信息构成训练样本集；对提取的特征进行采样生成训练样本数据：以这两幅原始图像中的像素及所述像素的ρ×ρ邻域内的像素所对应的步骤(4)中提取到的特征数据组成训练样本数据和测试样本数据，其中ρ取奇数值，可取值为1,3,5；

(7)选择分类器并初始化分类器的模型参数；

(8)判断是否满足终止条件，如果满足则执行(10)，否则执行(9)；

(9)根据步骤(6)得到的训练样本数据及对应的伪标签和组信息，用基于迭代隐式正则的方法训练分类器，训练完后执行步骤(8)；

(10)用分类器对步骤(6)得到的测试样本数据进行分类得到最终的融合结果图像。

2.根据权利要求1所述的方法，其中步骤(2)所述的采用SLIC进行超像素提取的步骤如下：

2a)随机初始化种子点，所述种子点称为聚类中心：按照设定的超像素个数，在这两幅原始图像内均匀地分配种子点；

2b)计算每一个种子的n×n邻域内所有像素点的梯度值，并用该种子点移动到该邻域内梯度最小点替代该种子点作为新的种子点，该该邻域内梯度最小点的像素值作为新的种子点的像素值；

2c)给每个种子点周围的邻域内的每个像素点分配类标签，即确定该邻域内每个像素点属于哪个聚类中心；

2d)对每个种子点周围的邻域内的每个像素点进行颜色距离和空间距离度量，分别计算每个种子点周围的邻域内的每个像素点到对应种子点的颜色距离和空间距离；

2e)判断是否达到预设的最大迭代次数，如果达到最大迭代次数，则执行步骤2f)，否则执行步骤2b)；

2f)进行区域连通操作：新建一张标记表，表内元素设置为-1，按照“Z”型走向将不连续的超像素、尺寸过小超像素重新分配给邻近的超像素，直到所有的像素点都遍历完毕并获得相应的标签。

3.根据权利要求1所述的方法，其中步骤(3)所述的采用模糊C均值聚类(FCM)方法对超像素进行聚类的步骤如下：

3a)随机初始化C个聚类中心；

3b)计算当前所有样本数据对C个聚类中心的每个聚类中心的隶属度，组成隶属度矩阵U，其中隶属度u_ij的计算公式为：

其中，C为聚类数目；m为模糊因子，取值可为2；d_ik表示第i个样本数据到第k个聚类中心的欧式距离；

3c)根据聚类中心C和隶属度矩阵U计算当前迭代次数的代价函数值J，并得到相邻两次迭代的代价函数值J的改变量σ，其中J的计算公式如下：

其中，N表示需要聚类的样本数据的总个数，C表示聚类数目，u_ij表示第i个数据属于第j类的隶属度，d_ik表示第i个样本数据到第j个聚类中心的欧式距离；

3d)根据改变量σ和当前迭代次数判断是否满足代价函数值J的改变量σ小于预设阈值/迭代次数达到预设最大迭代次数，满足则输出隶属度矩阵U，否则执行步骤3e)；

3e)采用梯度下降法计算新的C个聚类中心，并将迭代次数加1，返回步骤3b)，其中第j个聚类中心c_j的计算方式如下：

其中，x_i表示第i个待聚类数据；u_ij表示第i个数据属于第j类的隶属度；m表示模糊因子，取值可为2；N表示待聚类数据总数目。

4.根据权利要求1所述的方法，其中步骤(4)所述的根据原始两幅图像进行特征提取的步骤如下：

4a)对X₁和X₂图像做相对差值运算得到特征矢量，X₁和X₂图像中的每个像素的像素值矢量分别表示为X₁＝(g₁,g₂,...,g_n)^T，X₂＝(h₁,h₂,...,h_n)^T，其中n为波段数，变化矢量用公式可以表示为：

4b)：计算变化强度|△V|，其计算方式为：

4c)：计算变化方向(cosθ₁,cosθ₂,...,cosθ_n)：

4d)：由步骤4b)提取的强度特征和步骤4c)提取的方向特征构成最终的总特征(|△V|，cosθ₁,cosθ₂,...,cosθ_n)，

其中，|△V|为步骤4b)提取的强度特征，cosθ_i表示第i个波段处特征的变化方向的余弦值。

5.根据权利要求1所述的方法，其中步骤(9)所述的用基于迭代隐式正则的方法训练神经网络分类器的步骤如下：

9a)随机初始化神经网络的参数W，b，初始化正则项参数λ、μ；

9b)根据神经网络分类器的模型参数计算每个训练样本x_i的代价值l_i，并根据隐式正则项和l_i、λ、μ计算每个训练样本x_i的样本权重v_i；

9b1)根据模型参数计算每个训练样本的损失值计算如下：

其中为第j组内位于第i个训练样本的标签信息的one-hot编码形式，为第j组内位于第i个训练样本经过神经网络的输出层的输出结果；

9b2)将所有样本损失值进行归一化处理并在各自的组内按照升序方式进行排列；

9b3)根据隐式正则项推导出样本权重的计算方式，带隐式正则的神经网络的目标函数为：

其中表示第j组内第i个样本的代价值，b为总的组数，n_j表示第j组样本总数，

第二项为隐式正则项，计算方式为：

将新的目标函数E(W,v)并求对v的偏导数令该偏导数为零，即得到样本的权重计算公式为：

其中，表示第j组内第i个样本的权重值；λ和γ分别表示两个步长参数；表示第j组的第i个样本的损失值；C为一个常数值，用来控制组内样本的权重分布，一般取值为1；π为圆周率值；

9b4)根据该权重计算方法计算所有训练样本的权重；

9c)使用带权重的样本训练神经网络，神经网络的训练过程如下：

9c1)初始化神经网络的模型参数W,b；

9c2)神经网络进行前向传播，根据训练样本计算神经网络的每层结点输出，计算当前迭代次数的总的代价值l，其中N为训练样本总数，l的计算方式为

9c3)神经网络进行反向传播计算模型参数的梯度；

9c4)根据梯度和样本权重V更新模型参数W；

9c5)判断是否当前迭代次数是否达到最大迭代次数，若达到，则输出模型参数W,b,否则，返回执行5c2)。

6.根据权利要求1所述的方法，在两幅原始的多光谱图像内均匀地分配种子点采用随机选取两幅原始的多光谱图像内的预定设置数量的像素点作为种子点。

7.根据权利要求1所述的方法，n＝3。

8.根据权利要求1所述的方法，预设的最大迭代次数为10。

9.根据权利要求1所述的方法，其中，“Z”型走向指从左到右，从上到下顺序。

10.根据权利要求1所述的方法，随机初始化C个聚类中心包括从所有的数据点中任取C个点作为初始聚类中心。