高光谱图像分类 - HybridSN

HybridSN: Exploring 3-D–2-D CNN Feature Hierarchy for Hyperspectral Image Classification

目录

前言

目前的研究方向还是以 HSIC (HyperSpectral Image Classificaition) 为主了。开始整理看过的论文,第一篇就从 HybridSN 开始。

HybridSN: Exploring 3-D–2-D CNN Feature Hierarchy for Hyperspectral Image Classification

S. K. Roy, G. Krishna, S. R. Dubey and B. B. Chaudhuri, "HybridSN: Exploring 3-D–2-D CNN Feature Hierarchy for Hyperspectral Image Classification," in IEEE Geoscience and Remote Sensing Letters, vol. 17, no. 2, pp. 277-281, Feb. 2020, doi: 10.1109/LGRS.2019.2918719.

说白了 HybridSN 其实就是混合了 2D-CNN 和 3D-CNN,而之所以要混合,就是因为要“取其精华去其糟粕”。正如 Paper 中引出 HybridSN 提出动机所讲:

" The 2D-CNN alone isn’t able to extract good discriminating feature maps from the spectral dimensions. Similarly, a deep 3D-CNN is more computationally complex and alone seems to perform worse for classes having similar textures over many spectral bands."

论文思路

高光谱图像基础

尽管这篇文章没有过多地去介绍高光谱图像(HSI)的相关知识,但是作为第一篇整理的论文,还是在这里简单说一下。如果要展开来说,那有太多的基础知识去补充了,因为当初去了解的时候就是“一层一层挖”,比如:高光谱图像什么意思?这个"高"高在哪里?多光谱图像又是什么?怎么区别高光谱图像和多光谱图像?高光谱图像怎么拍出来的?.....

image.png

我们先来看一下上面这张照片,这是我在 @condifood 这个公司官网找到的图,从左到右分别是单色图像(也就是我们常说的灰度图像,即只含亮度信息,不含色彩信息,该图像的每个像素的 R/G/B 数值都是相同的)、彩色图像、高光谱图像。单从颜色这个角度去看,好像第三张高光谱图像也有颜色(深板岩暗蓝灰色?),实际上这是 False Color,也就是伪彩色图像,实际上严格来说第二张图片有三个通道,而第三张图片则是由上百个通道,我们看到的只是选择了三个通道进行显示所呈现出来的图像,这是我想给大家带来高光谱图像的第一印象。

第二:高光谱图像是由成像光谱仪拍的。而高光谱遥感技术即:利用成像光谱仪,在电磁波的可见光,近红外,中红外和热红外波段范围内,以纳米级光谱仪在几十个或几百个波段同时对地物进行成像,并记录地物的连续光谱响应,实现了地物空间信息,辐射信息,光谱信息的同步获取。(PS:成像光谱仪的原理就不说了,有兴趣可以自行 search)

第三:因为物体反射、吸收、透射和辐射电磁波的能力不同,所以我用着成像光谱仪(里面有不同的电磁波)去照地物表面的时候,地物会对这些电磁波有不同的反射,我就可以根据这些反射的不同去区分地表有什么不同的物质。这些反射的信息就是高光谱图像中丰富的光谱信息,在高光谱图像中的第三个维度得以体现,最后就构成了高光谱立体数据。关于立体,看下面这张图就非常清晰了。

image.png

说白了,高光谱图像和我们平时看到的照片有什么不同呢?空间特征(宽、高)还是一样的,只是你所看到的颜色(在这里或者已经不能说是颜色了),图像的通道数已经不是普通照片的 3 个,而是几百个,我们所看到的很多高光谱图像数据集,虽然说是数据集,其实每一个数据集就是一张高光谱图像,但是里面包含的数据很丰富。

高光谱图像分类在分类什么?

经过上面的些许基础知识普及,可能又有一个疑问,数据集其实就是一张图片,那分类什么?

上面说过高光谱图像有着丰富的信息,而对应到每一个像素,也就是 $1×1×通道数$ ,我们所要做的就是搞清楚这个像素对应的是地物的什么类别。也就是说:高光谱分类其实就是对每一个像素进行分类。这个像素属于玉米?大豆?树林?最后的分类类别取决于你这个高光谱数据集在哪里拍的,那个地方有什么地物分类。

思路概括

像上面所说,HybridSN 的提出并不是凭空提出,是因为在 HybridSN 提出之前,2D-CNN 和 3D-CNN 的 Application 在高光谱图像分类中各有优点和弊端:

2D-CNN 简单但是只能 capture 空间特征,3D-CNN 虽然对空、光谱特征都能 capture 到,并且能够提高分类的准确率,但是计算量上来了。

所以 HybridSN 就有提出的 Motivation,我可以结合 2D-CNN 和 3D-CNN 的优点,在 HybridSN 中,先使用三维卷积再堆叠二维卷积,最后连接分类器。这样就能很好地提取空间特征和光谱特征,计算量又不会很大。

有了基本思路,作者就具体提出 HybridSN 的网络结构(怎么混合?),然后就是实验分析,最后便是总结。

网络结构

image.png

上图即为原论文中 HybridSN 的网络结构,原始高光谱立方体数据集输入,经过 PCA 降维得到空间特征相同,但是光谱通道数降为 $B$ 的立方体数据集,接着就是邻域像素提取,分割成一个个小立方体数据集。

接着就是 HybridSN 的核心部分:经过三层 3 维卷积,再经过一层 2 维卷积,2 层全连接层,1 层 SoftMax 分类层。

就拿 Indian Pines 数据集来说,原始输入为: $145×145×200$ ,首先进行 PCA 降维,降到 $145×145×30$ ,我们要对每个像素进行分类,所以要分割成每一个像素进行后续的卷积操作,可是分割成每一个像素进行卷积之后,卷着卷着之后大小就很小了,所以对原始图像每个像素点进行填充,将其封装成一个邻域像素块,接下来就是对这个邻域像素块进行卷积操作。

先不管这里面如何操作的,总之 Indian Pines 的输入,经过 PCA 降维后进行分割像素填充邻域像素块后,可以分割成:10249 个这样 $25×25×30$ 这样的小立方体,其中 25 是作者试验出来的参数,也就是说最后的小立方体宽高都是 25 个像素,30 也就是 PCA 降维之后保留的光谱维度。

初始准备工作完成(准备要卷积的小立方体块),接着就是先 3 层 三维卷积,接着就是 1 层 二维卷积,在二维卷积前,reshape 一下,堆叠其后两个维度,从而符合二维卷积的形状。接着经过二维卷积之后,进行一个 Flatten(拉平) 送入第一个全连接层,然后再送入第二个连接层,最后送入 Softmax 进行分类,因为最后 Indian Pines 最后地物分类有 16 种,所以最后一层有 16 个神经元。

实验结果

高光谱图像数据集介绍

在这篇 Paper 中,使用了高光谱图像的三个主要数据集:

趁这个机会,我也稍微介绍一下高光谱的数据集。高光谱数据集不多,因为采集难度高,标注也很麻烦,比较经典的其实就是这篇论文所用到的三个,当然还有其它数据集。在阅读 HSIC 的论文的时候,你会发现这几个数据集出镜率非常高,基本上以后我们自己发论文也是依靠这么几个数据集。

如你所见,这些数据集的命名都是地名,事实也确实如此,比如说:Indian Pines

其他数据集类似,另外我国这些年也有自己的高光谱数据集了,大家可以去搜索一下,也可以参考我粗略整理的高光谱数据集

下面是 HybridSN 论文中所使用的三个数据集的主要信息:

数据集 空间维度 光谱波段 地物类别数
Indian Pines(IP) $145 \times 145$ 200 16
University of Pavia(UP) $610 \times 340$ 103 9
Salinas Scene(SA) $512 \times 217$ 204 16

关于高光谱图像数据集,在这里想说的一点就是,因为拍摄难度大,标注麻烦,所以可能来来去去就那么十个左右,换句话来说就是可能你每次用来训练测试的数据集也就是这么“十张图片”,只是这 10 张图片内容非常丰富,数据特别大。此外每个数据集基本上都是用不同的成像仪去拍摄,拍摄的场景也不一样,所以地物类别也不一样,图像尺寸也不一样,光谱通道数也不一样。

评价指标

使用 OA、AA、KA 来判断高光谱分类算法的性能。

  1. Overal Accuracy(OA):正确分类的样本个数占总测试样本数的比例;
  2. Average Accuracy(AA):各类准确率的均值;
  3. Kappa:用来判定结果的一致性,即预测的结果是否与真实的结果一致,其值越高表明一致性越好;

分类结果分析

作者使用提出的 HybridSN 与其他比较经典的高光谱图像分类方法作比较(包括机器学习方法与深度学习方法),有:SVM、2D-CNN、3D-CNN、M3D-CNN、SSRN。实验结果如下(来自论文中的图片);

image.png

这是作者对结果的分析:

“Table II shows the results in terms of the OA, AA and Kappa coefficient for different methods. It can be observed from Table II that the HybridSN outperforms all the compared methods over each dataset while maintaining the minimum standard deviation. The proposed HybridSN is based on the hierarchical representation of spectral-spatial 3D CNN followed by a spatial 2D CNN, which are complementary to each other. It is also observed from these results that the performance of 3D-CNN is poor than 2D CNN over Salinas Scene dataset. To the best of the our knowledge, this is probably due to the presence of two classes in the Salinas dataset (namely Grapes-untrained and Vinyard-untrained) which have very similar textures over most spectral bands. Hence, due to the increased redundancy among the spectral bands, the 2D-CNN outperforms the 3D-CNN over Salinas Scene dataset. Moreover, the performance of SSRN and HybridSN is always far better than M3D-CNN. It is evident that 3D or 2D convolution alone is not able to represent the highly discriminative feature compared to hybrid 3D and 2D convolutions

整体来说有这么几点:

  1. HybridSN 在三个数据集上的表现优于其他分类方法(其次就是 SSRN 了),并且保持最小的标准差,另外可以看到 HybridSN 在 Salinas Scene 数据集中上三个评价指标达到了 100%(离谱!);
  2. 并不是说 3D-CNN 就比 2D-CNN 好,这一点在 Salina Scene 数据集上有所表现:可以看到 2D-CNN 方法在 Salina Scene 数据集中表现得更好。造成这样的原因是虽然高光谱数据集如果考虑光谱特征整体分类效果会更加,但是 Salina Scene 数据集中有两个类的样本(Grapes-untrained and Vinyard-untrained)在大多数波段上有非常相似的纹理特征,而这两个类的样本数在总体可用样本数中又占比较大,所以 3D-CNN 在提取光谱特征的时候可能导致了信息冗余;
  1. 此外,作者拿 SSRN 和 HybridSN 结果与 M3D-CNN 作对比,更加证明了作者提出 HybridSN 的动机,单独使用二维卷积或者三维卷积,还不如混合两者:Mixing is better than using alone.

其他实验

👉 而在可视化分类效果这一部分,作者对比了 HybridSN 和其他方法的分类 Map

image.png

可以看出最后分类效果最佳的还是 SSRN 和 HybridSN,并且作者说在某些小部分 HybridSN 分类显示效果比 SSRN 更好。我们拿真实标签值与 SSRN 和 HybridSN 的 classification map 来看看有什么区别:

如作者所说,SSRN在某些像素分类上没有 HybridSN 强(如上图箭头所指区域像素所示等),但是总体情况上确实是不相上下。


👉 在分析 HybridSN 在 Indian Pines 数据集上训练过程中 100 个 epoch 的 Accuracy 和 Loss ,给出了如下变化曲线:

image.png

收敛大概是在 50 个 epoch 内实现的,从而表现出了 HybridSN 的快速收敛能力。


👉 在探讨 HybridSN 的计算效率上,作者对比 HybridSN 与其他分类方法的训练和测试时间:

image.png

作者想表达的是混合了 2D-CNN 和 3D-CNN 的 HybridSN 计算效率是比单独使用 3D-CNN 要快的,当然肯定是说比单独使用 2D-CNN 是要慢的(毕竟多了一个维度),这也证明了作者提出的 HybridSN 的效率可行性。


👉 在分割大立方体得到小立方体的时候,作者对像素进行了邻域填充分割,而这个邻域尺寸并不是说凭空得出是:patchSize=25,作者是提出了几个尺寸然后做实验得出最佳尺寸是 25,如下图所示:

image.png

虽然有些尺寸不是 25 的 Window Size 在某些数据集上的表现的比 Window Size 为 25 的要好,但是总体来说,Window Size 为 25 是最适合 HybridSN 的。


👉 在考虑使用更少的训练数据下(仅使用了总样本的 10% ):

image.png

其实仔细看可以发现,在使用了更少的训练数据下,SSRN 在 Indian Pines 的 OA 竟然比 HybridSN 高,在其他数据集上也有一些指标会比 HybridSN 要高,这就有点意思了。这其实归功于 SSRN 的 residual block,我们在整理 SSRN 这篇论文的时候详细说说。

整体来看,在样本数更少的情况下,HybridSN 和其他网络的性能都会下降。但是HybridSN 方法整体依然比其他方法好。

代码复现

作者提供的代码使用 Tensorflow 实现的,这里使用 Pytorch 进行复现。

导入库

import numpy as np
import scipy.io as sio
import spectral
from sklearn.decomposition import PCA
from sklearn.metrics import classification_report, accuracy_score, confusion_matrix, cohen_kappa_score
from sklearn.model_selection import train_test_split
import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
from operator import truediv

HybridSN 模型定义

首先就是对整个模型的定义,这部分代码非常清晰,在这边直接给出,不再过多赘述。Pytorch 比较麻烦的是你要知道每层卷积之后的具体计算,这里也不再过多计算,顺便贴上论文中作者给出的每一层的 Output Shape,方便查看

# 3 层 三维卷积,一层二维卷积,三层全连接层
class HybridSN(nn.Module):
def __init__(self, num_classes=16):
super(HybridSN, self).__init__()
# 第一层三维卷积:Conv3d_1, Input:(1, 30, 25, 25), 8 个 7 × 3 × 3 的卷积核
# Output: (8, 24, 23, 23)
self.conv1 = nn.Sequential(
nn.Conv3d(1, 8, (7, 3, 3)),
nn.ReLU()
)
# 第二层三维卷积:Conv3d_2,Input:(8, 24, 23, 23),16 个 5 × 3 × 3 的卷积核
# Output: (16, 20, 21, 21)
self.conv2 = nn.Sequential(
nn.Conv3d(8, 16, (5, 3, 3)),
nn.ReLU(inplace=True)
)
# 第三层三维卷积:Conv3d_3, Input(16, 20, 21, 21),32 个 3 × 3 × 3 的卷积核
# Output: (32, 18, 19, 19)
self.conv3 = nn.Sequential(
nn.Conv3d(16, 32, (3, 3, 3), stride=1, padding=0),
nn.ReLU(inplace=True)
)
# 第四层是二维卷积,但是要二维卷积之前就要把第三层的输出 reshape 成可以二维卷积的形状
# 其实就是样本数和维度数堆叠:32 × 18 = 576,所以这一层二维卷积的输入为 576
# Output: (64, 17, 17)
self.conv4_2d = nn.Sequential(
nn.Conv2d(574, 64, (3, 3)),
nn.ReLU(inplace=True)
)
# 二维卷积层之后就是一个全连接层,而在送入这个全连接层之前,要把二维卷积之后的输出降为一维才能送入全连接层
# 所以在 forward 里直接使用 flatten(),那么输入全连接层参数即为:64 × 17 × 17 = 18496
# 这一全连接层有 256 个节点
self.dense1 = nn.Linear(18496, 256)
# 接着又是一个全连接层,有 128 个节点
self.dense2 = nn.Linear(256, 128)
# 最终输出层,有 num_classes 这么多节点,因为我们最后就是要分 num_classes 这么多地物类别(Indian Pines)数据集
self.dense_out = nn.Linear(128, num_classes)
# 论文中使用了 Dropout ,参数为 0.4
self.drop = nn.Dropout(0.4)
# 定义前向传播函数
def forward(self, x):
out = self.conv1(x) # 第一层三维卷积
out = self.conv2(out) # 第二层三维卷积
out = self.conv3(out) # 第三层三维卷积
out = self.conv4_2d(out.reshape(out.shape[0], -1, 19, 19)) # 第四层二维卷积
out = out.reshape(out.shape[0], -1)
out = F.relu(self.dropout(self.dense1(out)))
out = F.relu(self.dropout(self.dense2(out)))
out - self.dense_out(out)
return out
image.png

数据预处理

上面说过高光谱原始数据集是一个大立方体,我们要对其中每个像素值进行分类,所以要切割成一个个小的立方体。当然,如果按照 $1×1$ 的像素去切割的话,经过卷积之后就“所剩无几”了,所以我们选择邻域像素提取成一个个小立方体。

具体步骤是:

上面的预处理已经说的非常详细了,下面的代码我也写上了非常详细的注释:

# 对高光谱数据集 X 应用 PCA 变换
def applyPCA(X, num_components):
newX = np.reshape(X, (-1, X.shape[2]))
# print(type(newX))
# print(newX.shape) # (145*145=21025, 200)
pca = PCA(n_components=num_components, whiten=True)
newX = pca.fit_transform(newX)
# print(newX.shape) # (21025, 30)
newX = np.reshape(newX, (X.shape[0], X.shape[1], num_components))
# print(newX.shape)
# print("------")
return newX
# 对单个像素周围提取 patch 时,边缘像素就无法取了,因此,给这部分像素进行 padding 操作
# 经过 padWithZeros 之后,宽、高会变成 169,169, 也就是左右各 padding 了 12
def padWithZeros(X, margin=12):
# print(X)
newX = np.zeros((X.shape[0] + 2 * margin, X.shape[1] + 2 * margin, X.shape[2]))
#print(newX.shape) # (169, 169, 30) 全是 0
x_offset = margin
y_offset = margin
newX[x_offset:X.shape[0] + x_offset, y_offset:X.shape[1] + y_offset, :] = X
return newX
# 对单个像素周围提取 patch 时,边缘像素就无法取了,因此,给这部分像素进行 padding 操作
def createImageCubes(X, y, windowSize=5, removeZeroLabels=True):
# 给 X 做 padding
margin = int((windowSize - 1) / 2) # 12
zeroPaddedX = padWithZeros(X, margin=margin) # (169, 30, 30)
# split patches
patchesData = np.zeros((X.shape[0] * X.shape[1], windowSize, windowSize, X.shape[2]))
# print(patchesData.shape) # (21025, 25, 25, 200)
patchesLabels = np.zeros((X.shape[0] * X.shape[1]))
# print(patchesLabels.shape) # (21025, )
patchIndex = 0
for r in range(margin, zeroPaddedX.shape[0] - margin): # for r in range(12, 169-12)
for c in range(margin, zeroPaddedX.shape[1] - margin): # for c in range(12, 169-12)
patch = zeroPaddedX[r - margin:r + margin + 1, c - margin:c + margin + 1]
# print(patch.shape)
# zeroPaddedX[r - 12 : r + 12 + 1, c -12 : c + 12 + 1]
patchesData[patchIndex, :, :, :] = patch
patchesLabels[patchIndex] = y[r - margin, c - margin]
patchIndex = patchIndex + 1
# 如果没有 removeZeroLabels
# patchesData 的 shape 将会是 (21025, 25, 25, 30)
# patchesLabels 的 shape 将会是 (21025, )
# removeZeroLabels 的作用就是去掉 gt 标签集 groundtruth 中为 0 的数据,因为里面的数据值有 0~16,而刚好 1~16 刚好对应地物分类的 16 类
if removeZeroLabels:
patchesData = patchesData[patchesLabels > 0, :, :, :]
patchesLabels = patchesLabels[patchesLabels > 0]
patchesLabels -= 1
return patchesData, patchesLabels
# 随机划分训练集和测试集
# 这里最后取 testradio=0.9 也就是随机抽 90% 为测试集,10% 为训练集
def splitTrainTestSet(X, y, testRatio, randomState=345):
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=testRatio, random_state=randomState, stratify=y)
return X_train, X_test, y_train, y_test

构建训练集和测试集加载类

上面已经成功划分了训练集和测试集,但是不能直接拿来训练,因为他是不可迭代的,所以我们在下面创建可以迭代的训练集和测试集类

# 构建训练数据集
""" Training dataset"""
class TrainDS(torch.utils.data.Dataset):
def __init__(self):
self.len = Xtrain.shape[0]
self.x_data = torch.FloatTensor(Xtrain)
self.y_data = torch.LongTensor(ytrain)
def __getitem__(self, index):
# 根据索引返回数据和对应的标签
return self.x_data[index], self.y_data[index]
def __len__(self):
# 返回文件数据的数目
return self.len
# 构建测试数据集
""" Testing dataset"""
class TestDS(torch.utils.data.Dataset):
def __init__(self):
self.len = Xtest.shape[0]
self.x_data = torch.FloatTensor(Xtest)
self.y_data = torch.LongTensor(ytest)
def __getitem__(self, index):
# 根据索引返回数据和对应的标签
return self.x_data[index], self.y_data[index]
def __len__(self):
# 返回文件数据的数目
return self.len
# 创建 trainloader 和 testloader
trainset = TrainDS()
testset = TestDS()
train_loader = torch.utils.data.DataLoader(dataset=trainset, batch_size=128, shuffle=True, num_workers=2)
test_loader = torch.utils.data.DataLoader(dataset=testset, batch_size=128, shuffle=False, num_workers=2)

训练过程

# 选择 GPU 进行训练
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
# 网络放到GPU上
net = HybridSN().to(device)
criterion = nn.CrossEntropyLoss()
optimizer = optim.Adam(net.parameters(), lr=0.001)
# 开始训练
total_loss = 0
# net.train()
for epoch in range(100):
for i, (inputs, labels) in enumerate(train_loader):
inputs = inputs.to(device)
labels = labels.to(device)
# 优化器梯度归零
optimizer.zero_grad()
# 正向传播 + 反向传播 + 优化
outputs = net(inputs)
loss = criterion(outputs, labels)
loss.backward()
optimizer.step()
total_loss += loss.item()
print('[Epoch: %d] [loss avg: %.4f] [current loss: %.4f]' % (epoch + 1, total_loss / (epoch + 1), loss.item()))
print('Finished Training')

测试过程

# 模型测试 (添加 net.eval())
net.eval()
count = 0
for inputs, _ in test_loader:
inputs = inputs.to(device)
outputs = net(inputs)
outputs = np.argmax(outputs.detach().cpu().numpy(), axis=1)
if count == 0:
y_pred_test = outputs
count = 1
else:
y_pred_test = np.concatenate( (y_pred_test, outputs))
# 生成分类报告
classification = classification_report(ytest, y_pred_test, digits=4)
print(classification)

之前每次测试的时候结果都不一样,其实应该要在测试之前添加 net.eval() 进入测试模式,为什么会有这个影响呢?因为在全连接层中我们使用了 Dropout 随机丢弃神经元,如果测试模式中不添加 net.eval() 那么这个随机丢弃也会继续随机,相反加入了 net.eval() 之后可以固定住训练完的 Dropout 那么分类后的结果也是一样了。

可视化分类结果

def AA_andEachClassAccuracy(confusion_matrix):
counter = confusion_matrix.shape[0]
list_diag = np.diag(confusion_matrix)
list_raw_sum = np.sum(confusion_matrix, axis=1)
each_acc = np.nan_to_num(truediv(list_diag, list_raw_sum))
average_acc = np.mean(each_acc)
return each_acc, average_acc
def reports(test_loader, y_test, name):
count = 0
# 模型测试
for inputs, _ in test_loader:
inputs = inputs.to(device)
outputs = net(inputs)
outputs = np.argmax(outputs.detach().cpu().numpy(), axis=1)
if count == 0:
y_pred = outputs
count = 1
else:
y_pred = np.concatenate((y_pred, outputs))
if name == 'IP':
target_names = ['Alfalfa', 'Corn-notill', 'Corn-mintill', 'Corn'
, 'Grass-pasture', 'Grass-trees', 'Grass-pasture-mowed',
'Hay-windrowed', 'Oats', 'Soybean-notill', 'Soybean-mintill',
'Soybean-clean', 'Wheat', 'Woods', 'Buildings-Grass-Trees-Drives',
'Stone-Steel-Towers']
elif name == 'SA':
target_names = ['Brocoli_green_weeds_1', 'Brocoli_green_weeds_2', 'Fallow', 'Fallow_rough_plow', 'Fallow_smooth',
'Stubble', 'Celery', 'Grapes_untrained', 'Soil_vinyard_develop', 'Corn_senesced_green_weeds',
'Lettuce_romaine_4wk', 'Lettuce_romaine_5wk', 'Lettuce_romaine_6wk', 'Lettuce_romaine_7wk',
'Vinyard_untrained', 'Vinyard_vertical_trellis']
elif name == 'PU':
target_names = ['Asphalt', 'Meadows', 'Gravel', 'Trees', 'Painted metal sheets', 'Bare Soil', 'Bitumen',
'Self-Blocking Bricks', 'Shadows']
classification = classification_report(y_test, y_pred, target_names=target_names)
oa = accuracy_score(y_test, y_pred)
confusion = confusion_matrix(y_test, y_pred)
each_acc, aa = AA_andEachClassAccuracy(confusion)
kappa = cohen_kappa_score(y_test, y_pred)
return classification, confusion, oa * 100, each_acc * 100, aa * 100, kappa * 100
classification, confusion, oa, each_acc, aa, kappa = reports(test_loader, ytest, 'IP')
classification = str(classification)
confusion = str(confusion)
file_name = "classification_report.txt"
with open(file_name, 'w') as x_file:
x_file.write('\n')
x_file.write('{} Kappa accuracy (%)'.format(kappa))
x_file.write('\n')
x_file.write('{} Overall accuracy (%)'.format(oa))
x_file.write('\n')
x_file.write('{} Average accuracy (%)'.format(aa))
x_file.write('\n')
x_file.write('\n')
x_file.write('{}'.format(classification))
x_file.write('\n')
x_file.write('{}'.format(confusion))
X = sio.loadmat('Indian_pines_corrected.mat')['indian_pines_corrected']
y = sio.loadmat('Indian_pines_gt.mat')['indian_pines_gt']
height = y.shape[0]
width = y.shape[1]
X = applyPCA(X, num_components = num_components)
X = padWithZeros(X, patch_size//2)
# 逐像素预测类别
outputs = np.zeros((height,width))
for i in range(height):
for j in range(width):
if int(y[i,j]) == 0:
continue
else :
image_patch = X[i:i+patch_size, j:j+patch_size, :]
image_patch = image_patch.reshape(1,image_patch.shape[0],image_patch.shape[1], image_patch.shape[2], 1)
X_test_image = torch.FloatTensor(image_patch.transpose(0, 4, 3, 1, 2)).to(device)
prediction = net(X_test_image)
prediction = np.argmax(prediction.detach().cpu().numpy(), axis=1)
outputs[i][j] = prediction+1
if i % 20 == 0:
print('... ... row ', i, ' handling ... ...')
predict_image = spectral.imshow(classes = outputs.astype(int),figsize =(5,5))
image.png

模型改进

单独加入 BatchNorm 层

在原论文以及作者提供的代码中,在每层卷积之后(无论是三维卷积还是二维卷积)都没有 BatchNorm 层,而在深度学习中,BN 的确是奏效的,BN 层的引入可以加快收敛速度,防止“梯度弥散”,在这里就涉及到 BatchNorm 的很多基础知识,包括 BatchNorm 的原理等等,不在这里一一赘述。

在 CNN 中,BN 层引入应该要在非线性激活层之前,也就是在本代码中,应该在每层卷积中的 ReLU() 函数之前。

也就是对于 HybridSN 模型定义更改为:

# 3 层 三维卷积,一层二维卷积,三层全连接层
class HybridSN(nn.Module):
def __init__(self, num_classes=16):
super(HybridSN, self).__init__()
# 第一层三维卷积:Conv3d_1, Input:(1, 30, 25, 25), 8 个 7 × 3 × 3 的卷积核
# Output: (8, 24, 23, 23)
self.conv1 = nn.Sequential(
nn.Conv3d(1, 8, (7, 3, 3)),
nn.BatchNorm3d(8),
nn.ReLU()
)
# 第二层三维卷积:Conv3d_2,Input:(8, 24, 23, 23),16 个 5 × 3 × 3 的卷积核
# Output: (16, 20, 21, 21)
self.conv2 = nn.Sequential(
nn.Conv3d(8, 16, (5, 3, 3)),
nn.BatchNorm3d(16),
nn.ReLU(inplace=True)
)
# 第三层三维卷积:Conv3d_3, Input(16, 20, 21, 21),32 个 3 × 3 × 3 的卷积核
# Output: (32, 18, 19, 19)
self.conv3 = nn.Sequential(
nn.Conv3d(16, 32, (3, 3, 3), stride=1, padding=0),
nn.BatchNorm3d(32),
nn.ReLU(inplace=True)
)
# 第四层是二维卷积,但是要二维卷积之前就要把第三层的输出 reshape 成可以二维卷积的形状
# 其实就是样本数和维度数堆叠:32 × 18 = 576,所以这一层二维卷积的输入为 576
# Output: (64, 17, 17)
self.conv4_2d = nn.Sequential(
nn.Conv2d(574, 64, (3, 3)),
nn.BatchNorm2d(64),
nn.ReLU(inplace=True)
)
# 二维卷积层之后就是一个全连接层,而在送入这个全连接层之前,要把二维卷积之后的输出降为一维才能送入全连接层
# 所以在 forward 里直接使用 flatten(),那么输入全连接层参数即为:64 × 17 × 17 = 18496
# 这一全连接层有 256 个节点
self.dense1 = nn.Linear(18496, 256)
# 接着又是一个全连接层,有 128 个节点
self.dense2 = nn.Linear(256, 128)
# 最终输出层,有 num_classes 这么多节点,因为我们最后就是要分 num_classes 这么多地物类别(Indian Pines)数据集
self.dense_out = nn.Linear(128, num_classes)
# 论文中使用了 Dropout ,参数为 0.4
self.drop = nn.Dropout(0.4)
# 定义前向传播函数
def forward(self, x):
out = self.conv1(x) # 第一层三维卷积
out = self.conv2(out) # 第二层三维卷积
out = self.conv3(out) # 第三层三维卷积
out = self.conv4_2d(out.reshape(out.shape[0], -1, 19, 19)) # 第四层二维卷积
out = out.reshape(out.shape[0], -1)
out = F.relu(self.dropout(self.dense1(out)))
out = F.relu(self.dropout(self.dense2(out)))
out - self.dense_out(out)
return out

加入 BatchNorm 之后更快收敛了:

单独加入注意力机制

添加通道注意力机制和空间注意力机制(待补充)...

加入 BatchNorm 和 注意力机制

(待补充)

加入 SE-Net 块

(待补充)

使用其他数据集的小问题

这是我在翻看网上关于 HybridSN 资料发现的一个问题,原链接在此:HybridSN代码修改 , 大概意思就是说无论是作者所提供的代码还是本文章中上面所复现的代码,都是在拿 Indian Pines 数据集进行训练和测试的,可是上面的文章指出在 Pavid University 数据集上进行训练和测试的时候,内存会爆满,究其原因其实是在 patchesData 这部分出错了。我们重新拿分割大立方体的代码来看看:

# 对单个像素周围提取 patch 时,边缘像素就无法取了,因此,给这部分像素进行 padding 操作
def createImageCubes(X, y, windowSize=5, removeZeroLabels=True):
# 给 X 做 padding
margin = int((windowSize - 1) / 2) # 12
zeroPaddedX = padWithZeros(X, margin=margin) # (169, 30, 30)
# split patches
patchesData = np.zeros((X.shape[0] * X.shape[1], windowSize, windowSize, X.shape[2]))
# print(patchesData.shape) # (21025, 25, 25, 200)
patchesLabels = np.zeros((X.shape[0] * X.shape[1]))
# print(patchesLabels.shape) # (21025, )
patchIndex = 0
for r in range(margin, zeroPaddedX.shape[0] - margin): # for r in range(12, 169-12)
for c in range(margin, zeroPaddedX.shape[1] - margin): # for c in range(12, 169-12)
patch = zeroPaddedX[r - margin:r + margin + 1, c - margin:c + margin + 1]
# print(patch.shape)
# zeroPaddedX[r - 12 : r + 12 + 1, c -12 : c + 12 + 1]
patchesData[patchIndex, :, :, :] = patch
patchesLabels[patchIndex] = y[r - margin, c - margin]
patchIndex = patchIndex + 1
# 如果没有 removeZeroLabels
# patchesData 的 shape 将会是 (21025, 25, 25, 30)
# patchesLabels 的 shape 将会是 (21025, )
# removeZeroLabels 的作用就是去掉 gt 标签集 groundtruth 中为 0 的数据,因为里面的数据值有 0~16,而刚好 1~16 刚好对应地物分类的 16 类
if removeZeroLabels:
patchesData = patchesData[patchesLabels > 0, :, :, :]
patchesLabels = patchesLabels[patchesLabels > 0]

对于 patchesData 我们最先做的是初始化这么一个四维张量,并且 sample 数为 X.shape[0] * X.shape[1],Indian Pines 的 spatial demension 是 145 × 145,没有问题;而 Pavia University 的 spatial demension 是 512×217 ,已经达到了 111104,而 Indian Pines 才 21025,这无疑为内存带来了巨大的压力。

解决的方法如上面文章里面所说,上面我们讲解代码的时候也有所提及,最后真正的 patchesData 其实并非有 21025,因为最后的 Label 不是每一个 split 出来的小立方体都是在 1~16 这个类别中,也就是 label = 0,这部分我们需要进行丢弃,即代码中的 removeZeroLabels。所以,从这个角度去思考,我们是不是可以先去看看真正有用的 label 有多少,也不是先让内存去占 height * weight 这么多空间呢?

这就是解决的方法,先对这么多个小立方体进行一个循环,找到真正有类别的小立方体,再去开辟空间,改善后的代码如下:

# 在每个像素周围提取 patch
def createImageCubes(X, y, windowSize=5, removeZeroLabels = True):
# 给 X 做 padding
margin = int((windowSize - 1) / 2)
zeroPaddedX = padWithZeros(X, margin=margin)
# 获得 y 中的标记样本数
count = 0
for r in range(0, y.shape[0]):
for c in range(0, y.shape[1]):
if y[r, c] != 0:
count = count+1
# split patches
patchesData = np.zeros([count, windowSize, windowSize, X.shape[2]])
patchesLabels = np.zeros(count)
count = 0
for r in range(margin, zeroPaddedX.shape[0] - margin):
for c in range(margin, zeroPaddedX.shape[1] - margin):
if y[r-margin, c-margin] != 0:
patch = zeroPaddedX[r - margin:r + margin + 1, c - margin:c + margin + 1]
patchesData[count, :, :, :] = patch
patchesLabels[count] = y[r-margin, c-margin]
count = count + 1
return patchesData, patchesLabels

参考