CN113052130B - 基于深度残差网络和边缘保护滤波的高光谱影像分类方法 - Google Patents

Abstract

本发明公开了一种基于深度残差网络和边缘保护滤波的高光谱影像分类方法，涉及高光谱图像处理技术领域。该方法包括：利用主成分分析的方法对原始高光谱影像数据进行特征提取，将第一主成分特征作为引导图像，前三主成分特征作为输入图像；通过联合双边滤波提取出影像特征后进行图像叠加；构建一种深度残差网络模型并采用该模型对原始空谱特征进行处理提取出影像的深层空谱特征；将提取的空谱特征输入到Softmax分类器完成影像分类。该方法综合了联合双边滤波与深度学习模型的优势，缓解了卷积神经网络分类中的过拟合问题，在顾及边缘信息的基础上能够充分提取和利用影像的深层空谱特征，有效提高高光谱影像的地物分类精度。

Description

基于深度残差网络和边缘保护滤波的高光谱影像分类方法

技术领域

本发明涉及高光谱图像处理技术领域，涉及一种基于深度残差网络和边缘保护滤波的高光谱影像分类方法。

背景技术

高光谱遥感影像具有较高的光谱分辨率，包含了丰富的光谱信息和空间信息。高光谱影像分类作为高光谱遥感图像应用的重要方面，已广泛应用于环境检测、矿物填图、精细农业等领域。然而由于高光谱影像的相邻波段间相关度强、空间光谱结构复杂性高和训练样本数量有限等特点给分类任务带来了巨大的挑战。

在早期的遥感影像分类中，中值滤波、均值滤波、高斯滤波等一系列方法利用影像的光谱信息进行特征提取。大量研究发现，在高光谱影像分类任务中空间信息往往比光谱信息更为关键，仅使用光谱信息导致无法达到足够高的分类精度。因此，研究人员将空间信息融入到高光谱影像分类中，提出了空-谱特征融合的高光谱遥感影像分类方法，将空间特征作为光谱特征的一种互补。然而绝大多数特征提取方法虽然有效地利用了影像的空间特征，在一定程度上提高了影像的分类精度，但是不能很好地保护影像的边缘信息，而这些边缘信息很有可能是区别不同地物种类的关键特征，忽略边缘信息对地物进行分类往往引起地物边缘处的错分问题。

深度学习方法构建的网络模型包含复杂的多层非线性变换结构可以学习高度表达的特征，适用于处理复杂的地物分类问题，近年来广泛应用于高光谱遥感影像分类领域中。堆叠式自动编码器(SAE)最早用于提取高光谱影像的空-谱特征，由于SAE只能将输入数据转换为一维向量，导致部分空间信息受损。卷积神经网络(CNN)通过共享卷积核能够有效地处理高维数据，在减少模型参数的同时提高分类性能。由于高光谱影像存在一定的特殊性，其小样本、结构复杂等的特点导致难以为卷积神经网络提供足够多的训练样本，易产生过拟合现象，同时过多的网络层数会引起梯度消失、网络退化问题。

发明内容

本发明为了解决上述问题，提供一种基于深度残差网络和边缘保护滤波的高光谱影像分类方法，旨在顾及边缘信息的基础上能够充分提取和利用影像的深层空谱特征，有效提高高光谱影像的地物分类精度。

本发明的技术方案是：

一种基于深度残差网络和边缘保护滤波的高光谱影像分类方法，包括以下步骤：

步骤1：输入待分类的高光谱遥感影像原始数据集；

步骤2：对输入的高光谱遥感影像进行归一化处理，得到归一化后的高光谱影像数据集；

步骤3：利用主成分分析法对归一化后的影像进行降维，得到前三个主成分影像作为输入图像，第一主成分影像作为引导图像；

步骤4：通过步骤3得到的引导图像对输入图像进行联合双边滤波处理，提取出各个主成分的空间特征，得到每个主成分对应的输出图像；

步骤5：将步骤4滤波后的各个主成分输出图像进行叠加后，以每个目标像元为中心，提取大小为M×M×d的空间-光谱立方体及其类别l作为样本数据集；其中M×M为空间邻域大小，d为波段数目；

步骤6：将步骤5得到的空谱特征划分为训练样本和测试样本，将训练样本输入到深度残差网络中并采用Softmax线性分类器进行分类。

与现有技术相比较，本发明具有如下有益效果：考虑到高光谱影像波段间相关度强、光谱和空间结构复杂性高和训练样本数量有限等特性，本发明方法在空间特征提取时采用联合双边滤波增强影像中地物的边缘结构，提取出高质量的空间特征并且与光谱特征融合得到影像的空谱特征，在有限样本条件下有效地缓解了CNN训练过程中的过拟合问题，降低了CNN对大量训练样本的需求。同时，为避免由网络层数增加导致的网络退化和梯度消失的问题，在卷积神经网络中加入残差网络实现通过跳层连接在保持参数量不变的同时提高网络性能。因此，本发明在顾及边缘信息的基础上能够充分提取和利用影像的深层空谱特征，有效提高了高光谱影像的地物分类精度。

附图说明

图1是本发明基于深度残差网络和边缘保护滤波的高光谱影像分类方法流程图；

图2(a)为Pavia University数据集拍摄区域的部分伪彩色图；(b)为(a)图所示区域的地面参考图；

图3(a)为采用SP-SVM方法对Pavia University数据集进行分类的随机一次仿真实验结果图；(b)为采用APs-2DCNN方法对Pavia University数据集进行分类的随机一次仿真实验结果图；(c)为采用Gabor-2DCNN方法对Pavia University数据集进行分类的随机一次仿真实验结果图；(d)采用本发明方法对Pavia University数据集进行分类的随机一次仿真实验结果图。

具体实施方式

为了便于理解本申请，下面将参照相关附图和实施例对本申请进行更全面的描述。附图中给出了本申请的较佳实施方式。但是，本申请可以以许多不同的形式来实现，并不限于本文所描述的实施方式。相反地，提供这些实施方式的目的是使对本申请的公开内容理解的更加透彻全面。

图1是本发明提供的基于深度残差网络和边缘保护滤波的高光谱影像分类方法的流程图，所述基于深度残差网络和边缘保护滤波的高光谱影像分类方法，包括如下步骤：

步骤1：输入待分类的高光谱遥感影像原始数据集；

在本实施例中，选取Pavia University数据集作为待分类的高光谱遥感影像原始数据集，Pavia University数据是通过ROSIS传感器在Pavia大学拍摄的校园建筑与环境情况，包含9类地物，42776个样本。影像大小为6105pixel×3405pixel，共103个波段(不含噪声波段)，波长范围为430-860nm，空间分辨率为1.3m。图2(a)为该区域的伪彩色图(合成波段：90、60和30)，图2(b)为该区域的地面参考图。

步骤2：对输入的高光谱遥感影像进行归一化处理，得到归一化后的高光谱影像数据集为：

其中，u表示样本总数，D表示波段个数；

步骤3：利用主成分分析(PCA)的方法对归一化后的数据进行降维，得到前三个主成分影像作为输入图像，第一主成分影像作为引导图像，具体方法为：

步骤3.1：计算归一化后的数据样本的协方差矩阵，具体方法为：

步骤3.1.1：计算第j个波段中所有样本的均值，如下公式所示：

其中，x_ij表示第j个波段的第i个样本，1≤i≤u，1≤j≤D；

步骤3.1.2：对样本X_norm进行中心化处理，得到数据集为：

步骤3.1.3：通过式(4)计算影像的协方差矩阵；

步骤3.2：计算协方差矩阵的特征值和特征向量，并将前三主成分并作为输入图像，将第一主成分作为引导图像，具体方法为：

步骤3.2.1：对协方差矩阵进行特征值分解，获得特征值λ₁≥λ₂≥…≥λ_d，对应的特征向量为W＝[ξ₁,ξ₂,…,ξ_d]；

步骤3.2.2：通过特征值λ₁,λ₂,λ₃对应的特征向量W₃＝[ξ₁,ξ₂,ξ₃]得到前三主成分分量，采用前三主成分对样本X_cen进行线性转换，如式(5)所示，得到前三主成分影像，并将其作为输入图像；

X′_i＝W₃ ^TX_cen (5)

步骤3.2.3：通过特征值λ₁对应的特征向量W₁＝[ξ₁]得到第一主成分分量，采用第一主成分对样本X_cen进行线性转换，如式(6)所示，得到第一主成分影像，并将其作为引导图像；

X′_g＝W₁ ^TX_cen (6)

步骤4：通过步骤3得到的引导图像X′_g对输入图像X′_i进行联合双边滤波处理，提取出各个主成分的空间特征，得到每个主成分对应的输出图像，具体方法为；

步骤4.1：对于输入图像X′_i中第i个像素点r_i，计算其空域距离权重和值域距离权重/>具体方法为：

步骤4.1.1：在保证输出结果中存在引导图像X′_g的边缘信息基础上，设定空域距离和值域距离高斯函数标准差σ_c和σ_s；

步骤4.1.2：按照式(7)计算r_i的空域距离权重

其中，r_j∈η(r_i)表示在邻域η中，r_j是以r_i为中心的空间邻域点；

步骤4.1.3：按照式(8)计算r_i的值域距离权重

其中，和/>是引导图像中与r_i和r_j对应的像素点；

步骤4.2：由步骤4.1得到的空域距离权重和值域距离权重/>根据公式(9)计算归一化因子K_i；

步骤4.3：按照式(10)计算像素点r_i经滤波后对应的输出值；

其中，I_j为像素点r_i的邻域像素加权和；

步骤4.4：重复步骤4.1-4.3，直到输入图像中所有像素点都得到对应的输出值。由所有的输出J_i得出最终的滤波图像J；

步骤5：将步骤4滤波后的各个主成分输出图像进行叠加后，以每个目标像元为中心，提取大小为M×M×d的空间-光谱立方体及其类别l作为样本数据集；其中M×M为空间邻域大小，在本实施例中M的值为9；d为波段数目，在本实施例中该值为3。

步骤6：将步骤5得到的空谱特征划分为训练样本和测试样本，将训练样本输入到深度残差网络中并采用Softmax线性分类器进行分类，具体方法为：

步骤6.1：将步骤5得到的空谱特征随机划分为训练样本和测试样本；

在本实施例中，从Pavia University数据集的每类地物样本中选取50个样本(地物总数量不足50时，选取一半)作为训练样本，剩余样本作为测试样本。

步骤6.2：搭建卷积神经网络，提取影像的空谱特征，按照式(11)计算网络的第p个卷积层的第q个特征图在(a,b)处的特征值

其中，t为第p-1层中特征图的个数；M_p和N_q分别为卷积核的长度和宽度；为与(m,n)处相连的第t个特征图的权值；b_pq为第p层中第q个特征图的偏置；f(·)为激活函数；

步骤6.3：在每个卷积层、池化层和全连接层后引入ReLU作为网络的激活函数，如下公式所示：

ReLU(z)＝max(0,z) (12)

其中，z表示输入数据；

步骤6.4：在经过一个卷积层后连接一个自定义的残差网络模块，该残差网络模块为六个卷积层的密集连接，具体方法为：

步骤6.4.1：按照式(13)计算基本残差单元；

H(x_i)＝F(x_i)+x_i (13)

其中，x_i表示第i层的输入，H(x_i)为期望输出；

步骤6.4.2：在步骤6.3.1的基础上，引入残差模块中的权重矩阵{W_i}，在忽略偏差的前提下，按照式(14)计算残差映射F(x)；

F(x,{W_i})＝w₂μ(W₁x) (14)

其中，μ为ReLU激活函数；

步骤6.4.3：按照式(15)计算输出结果z；

z＝F(x,{W_i})+x (15)

其中，x和F(x,{W_i})需要保持相同的维数；

步骤6.4.4：在维数不相同时，引入线性投影W_j保持维数的一致，如下公式所示：

z＝F(x,{W_i})+W_jx (16)

在本实施例中，经过一个步长为1，大小为3×3卷积层后连接一个残差网络模块，该残差网络模块为六个卷积层的密集连接，其中卷积核的步长为1，大小为3×3；

步骤6.5：将步骤6.4得到的特征图连接一个预设的最大池化层；

在本实施例中，是将步骤6.4得到的特征图连接一个步长为2、池化窗口大小为3×3的最大池化层；

步骤6.6：将步骤6.5得到的特征图连接一个全连接层并采用Dropout正则化方法随机丢弃一定比例的节点避免过拟合的风险，最终输出一个256维的向量D；

步骤6.7：网络的输出层采用Softmax线性分类器进行逻辑回归分类，具体过程为：

步骤6.7.1：假设输入的数据集为X＝{(a₁,b₁),(a₂,b₂),...,(a_n,b_n)}，n为样本总数，对于第i个样本a_i所对应的类别标签表示为b_i(b_i∈{1,2,...,N}，其中，N表示类别总数)，按照式(17)计算样本a_i属于类别l的概率；

其中，为模型的参数矩阵，用于归一化图像输出特征，使输出概率之和为1；

步骤6.7.2：按照式(18)计算损失函数J(θ)；

步骤6.7.3：选择最大概率值所对应的类别作为最终结果，输出类别标签；

步骤6.8：在输出类别标签后使用分类交叉熵函数(categorical_crossentropy)构造目标函数，通过自适应矩估计(Adam)优化器对网络进行优化，初始学习率设置为0.001。

本实施例中，卷积神经网络中卷积核的数量、dropout概率值和网络迭代次数epoch，联合双边滤波中的参数有邻域大小η、空域距离高斯函数的标准差σ_c和值域距离高斯函数的标准差σ_s会影响分类效果，当第一个卷积层卷积核数量为32、dropout的值为0.6、当epoch的值为150、当η＝4、σ_c＝3和σ_s＝0.5时分类精度最优为97.83％。分类精度随着窗口半径η的增加而不断提高，但是随着η的增加空间信息达到饱和状态后出现了信息冗余，导致分类精度不再提高甚至有下降的趋势。

图3(a)、(b)、(c)、(d)依次为采用SP-SVM方法、APs-2DCNN方法、Gabor-2DCNN方法和本发明方法对Pavia University数据集进行分类的随机一次仿真实验结果图，由图3可以看出，SP-SVM由于仅利用了光谱特征，分类结果图中出现严重的“椒盐现象”；APs-2DCNN、Gabor-2DCNN和本发明方法引入了影像的空间特征，分类结果有明显的提升，但在地物的边缘处有较多的错分结果；本发明方法中噪声点更少、分类边界更平滑，较好的反映出地物的真实分布情况。

各方法对应的各类地物分类精度、总体分类精度(OA)、平均分类精度(AA)和Kappa系数如表1所示。本发明方法的分类结果OA，AA和Kappa系数分别约为99.35％，98.84％，99.13％，相比SP-SVM、APs-2DCNN和Gabor-2DCNN方法的分类结果均有所提高，OA分别提高了约10.77％、3.50％和2.57％，AA分别提高了约8.03％、4.08％和3.42％，Kappa系数分别提高了约14.60％、4.65％和3.40％，同时获取了8个最佳类别的地物分类精度。

表1四种方法的分类结果(％)

类别	SP-SVM	APs-2DCNN	Gabor-2DCNN	本发明方法
					沥青路面	81.44	95.09	96.08	99.38
草地	90.26	97.83	99.23	99.92
					砾石	82.63	93.27	94.96	96.32
树	95.65	99.24	97.74	99.45
					金属板	99.15	99.23	100	98.02
裸土	94.11	97.74	97.41	100
					柏油路面	94.77	83.96	85.22	98.84
砖块	79.26	86.52	88.15	97.67
					阴影	100	100	100	100
OA	88.58	95.85	96.78	99.35
					AA	90.81	94.76	95.42	98.84
Kappa	84.53	94.48	95.73	99.13

应当理解的是，本领域技术人员在本发明技术构思的启发下，在不脱离本发明内容的基础上，可以根据上述说明做出各种改进或变换，这仍落在本发明的保护范围之内。

Claims (4)

1.一种基于深度残差网络和边缘保护滤波的高光谱影像分类方法，其特征在于，包括以下步骤：

步骤1：输入待分类的高光谱遥感影像原始数据集；

步骤2：对输入的高光谱遥感影像进行归一化处理，得到归一化后的高光谱影像数据集；

步骤3：利用主成分分析法对归一化后的影像进行降维，得到前三个主成分影像作为输入图像，第一主成分影像作为引导图像；

步骤4：通过步骤3得到的引导图像对输入图像进行联合双边滤波处理，提取出各个主成分的空间特征，得到每个主成分对应的输出图像；

步骤5：将步骤4滤波后的各个主成分输出图像进行叠加后，以每个目标像元为中心，提取大小为M×M×d的空间-光谱立方体及其类别l作为样本数据集；其中M×M为空间邻域大小，d为波段数目；

步骤6：将步骤5得到的空谱特征划分为训练样本和测试样本，将训练样本输入到深度残差网络中并采用Softmax线性分类器进行分类；

所述步骤4包括如下步骤：

步骤4.1：对于输入图像X_i′中第i个像素点r_i，计算其空域距离权重和值域距离权重

步骤4.2：由步骤4.1得到的空域距离权重和值域距离权重/>计算归一化因子K_i；

步骤4.3：计算像素点r_i经滤波后对应的输出值；

步骤4.4：重复步骤4.1-4.3，直到输入图像中所有像素点都得到对应的输出值，从而得到最终的滤波图像J；

所述步骤4.1包括：

步骤4.1.1：在保证输出结果中存在引导图像X_g′的边缘信息基础上，设定空域距离和值域距离高斯函数标准差σ_c和σ_s；

步骤4.1.2：按照式(7)计算r_i的空域距离权重

其中，r_j∈η(r_i)表示在邻域η中，r_j是以r_i为中心的空间邻域点；

步骤4.1.3：按照式(8)计算r_i的值域距离权重

其中，和/>是引导图像中与r_i和r_j对应的像素点；

所述M的值为9，d为3；

所述步骤6包括如下步骤：

步骤6.1：将步骤5得到的空谱特征随机划分为训练样本和测试样本；

步骤6.2：搭建卷积神经网络，提取影像的空谱特征，计算网络的第p个卷积层的第q个特征图在(a,b)处的特征值

步骤6.3：在每个卷积层、池化层和全连接层后引入ReLU作为网络的激活函数；

步骤6.4：在经过一个卷积层后连接一个自定义的残差网络模块，该残差网络模块为六个卷积层的密集连接；

步骤6.5：将步骤6.4得到的特征图连接一个预设的最大池化层；

步骤6.6：将步骤6.5得到的特征图连接一个全连接层并采用Dropout正则化方法随机丢弃一定比例的节点避免过拟合的风险，最终输出一个256维的向量D；

步骤6.7：网络的输出层采用Softmax线性分类器进行逻辑回归分类。

2.根据权利要求1所述的基于深度残差网络和边缘保护滤波的高光谱影像分类方法，其特征在于，在步骤6.4中，在经过一个步长为1，大小为3×3卷积层后连接一个自定义的残差网络模块，所述六个卷积层的卷积核的步长为1，大小为3×3。

3.根据权利要求1所述的基于深度残差网络和边缘保护滤波的高光谱影像分类方法，其特征在于，在步骤6.5中所述的最大池化层的步长为2、池化窗口大小为3×3。

4.根据权利要求1所述的基于深度残差网络和边缘保护滤波的高光谱影像分类方法，其特征在于，所述步骤6.7包括如下步骤：

步骤6.7.1：假设输入的数据集为X＝{(a₁,b₁),(a₂,b₂),...,(a_n,b_n)}，n为样本总数，对于第i个样本a_i所对应的类别标签表示为b_i，b_i∈{1,2,...,N}，N表示类别总数，按照式(17)计算样本a_i属于类别l的概率；

其中，为模型的参数矩阵，用于归一化图像输出特征，使输出概率之和为1；

步骤6.7.2：按照式(18)计算损失函数J(θ)；

步骤6.7.3：选择最大概率值所对应的类别作为最终结果，输出类别标签。