文章目录
论文基本信息
- 标题：HybridSN: Exploring 3-D–2-DCNN Feature Hierarchy for Hyperspectral Image Classification
- 作者：Swalpa Kumar Roy, Student Member , IEEE, Gopal Krishna, Shiv Ram Dubey, Member , IEEE, and Bidyut B. Chaudhuri, Life Fellow, IEEE
- 时间：2019
- 论文地址：arXiv:1902.06701v3 [cs.CV] 3 Jul 2019
- code：https://github.com/OUCTheoryGroup/colab_demo/blob/master/202003_models/HybridSN_GRSL2020.ipynb
研究背景
高光谱图像分类广泛应用于遥感图像的分析。高光谱图像包括不同波段的图像。卷积神经网络(CNN)是一种常用的基于深度学习的视觉数据处理方法。CNN用于HSI分类的使用在最近的文章中也可以看到。这些方法大多基于2D CNN，然而，HSI分类性能高度依赖于空间和光谱信息。由于计算复杂度的增加，很少有方法利用3D CNN。本文提出了一种用于HSI分类的混合光谱卷积神经网络(HybridSN)。基本上，HybridSN是一个光谱-空间 3D-CNN，然后是空间 2D-CNN。3D-CNN促进了从光谱波段堆栈的联合空间-光谱特征表示。在3D-CNN之上的2D-CNN进一步学习了更抽象的层次空间表示。此外，使用混合CNN比单独使用3D-CNN降低了模型的复杂性。
读完摘要后的疑问
- 高光谱是啥？
- HybridSN的模型结构
读完文章后对上述问题的回答
高光谱基础知识
- 单色光：单一波长（或频率）的光，不能产生色散。
- 复色光：几种单色光合成的光。
- 色散系统：复色光分解为单色光而形成光谱的现象。
- 光谱（光学频谱，spectrum）：复色光经过色散系统（如棱镜、光栅）分光后，被色散开的单色光按波长（或频率）大小而依次排列的图案。
- 光栅：由大量等宽等间距的平行狭缝构成的光学器件。
- 高光谱图像（Hyperspectral Image）：在光谱的维度进行了细致的分割，不仅仅是传统的黑，白或者RGB的区别，而是在光谱维度上也有N个通道。例如：我们可以把400nm-1000nm分为300个通道，一次，通过高光谱设备获取的是一个数据立方，不仅有图像的信息，并且在光谱维度上进行展开，结果不仅可以获得图像上每个点的光谱数据，还可以获得任意一个谱段的影像信息。
- 高光谱图像成像原理：空间中的一维信息通过镜头和狭缝后，不同波长的光按照不同程度的弯散传播，一维图像上的每个点，再通过光栅进行衍射分光，形成一个谱带，照射到探测器上，探测器上的每个像素位置和强度表征光谱和强度。
- HSI：数字图像模型，反映了人的视觉系统感知彩色的方式，以色调H（Hue）、饱和度S（Saturation）、亮度I（Intensity）三种基本特征量来感知颜色。
- 色调H（Hue）：与光波的频率有关，它表示人的感官对不同颜色的感受，也可以表示一定范围的颜色；
- 饱和度S（Saturation）：表示颜色的纯度，纯光谱色是完全饱和的，加入白光会稀释饱和度。饱和度越大，颜色看起来就越鲜艳。
- 亮度I（Intensity）：对应成像亮度和图像灰度，是颜色的明亮程度。
HybridSN
根据上边的描述，2D的CNN很难充分地提取高光谱的空间信息，而3D的CNN则是由于计算量庞大而很少被使用。
本文将2D和3D的CNN进行了结合，使用3D的CNN高光谱的空间特征的提取，使用2D的CNN进行更抽象的层次空间表示。
网络结构如上图所示（3个三维卷积，1个二维卷积，3个全连接层）：
- 三维卷积中，卷积核的尺寸为8×3×3×7×1、16×3×3×5×8、32×3×3×3×16（16个三维核，3×3×5维）
- 二维卷积中，卷积核的尺寸为64×3×3×576（576为二维输入特征图的数量）
实验验证
数据集描述
三种开源高光谱图像数据集，Indian Pines（IP）, University of Pavia（UP） and Salinas Scene（SA）
IP图像空间维度为：145×145，波长范围为400-2500nm，共有224个光谱波段
UP图像空间维度为：610×340，波长范围430-860nm，共有103个光谱波段
SA图像空间维度为：512×217，波长范围360-2500nm，共有224个光谱波段
分类结果
使用了Overall Accuracy (OA)、Average Accuracy (AA)和Kappa Coefficient (Kappa)评价指标来判断HSI分类性能。其中，OA表示总体测试样本中正确分类的样本数;AA表示分类精度的平均值;Kappa是一种统计测量的度量，它提供了关于地面真相图和分类图之间的强一致性的相互信息。将所提出的HybridSN模型的结果与最广泛使用的监督方法如SVM[33]、2D-CNN[34]、3D-CNN[35]、M3D-CNN[36]和SSRN[27]进行比较。30%和70%的数据分别随机分为训练组和测试组2。我们使用了比较方法的公开代码3来计算结果。
从表二可以看出，混合算法在保持最小标准偏差的同时，在每个数据集上都优于所有的比较方法。
图3显示了混合网络在IP、UP和SA数据集上的HSI分类性能的混淆矩阵。
训练集和验证集100个epoch的准确性和损失收敛如图4所示。可以看出，收敛时间约为50个epoch，这表明作者的方法收敛速度较快。
混合动力网络模型的计算效率体现在训练和测试次数方面，如表3所示。
代码分析
三维卷积部分：
conv1：（1, 30, 25, 25）， 8个 7x3x3 的卷积核 ==>（8, 24, 23, 23）
conv2：（8, 24, 23, 23）， 16个 5x3x3 的卷积核 ==>（16, 20, 21, 21）
conv3：（16, 20, 21, 21），32个 3x3x3 的卷积核 ==>（32, 18, 19, 19）
接下来要进行二维卷积，因此把前面的 32*18 reshape 一下，得到 （576, 19, 19）
二维卷积部分：
（576, 19, 19） 64个 3x3 的卷积核，得到 （64, 17, 17）
接下来是一个 flatten 操作，变为 18496 维的向量，
接下来依次为256，128节点的全连接层，都使用比例为0.4的 Dropout，
最后输出为 16 个节点，是最终的分类类别数。
! wget http://www.ehu.eus/ccwintco/uploads/6/67/Indian_pines_corrected.mat
! wget http://www.ehu.eus/ccwintco/uploads/c/c4/Indian_pines_gt.mat
! pip install spectral
#获取数据，并导入基本库
import numpy as np
import matplotlib.pyplot as plt
import scipy.io as sio
from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix, accuracy_score, classification_report, cohen_kappa_score
import spectral
import torch
import torchvision
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
定义HybridSN 类：
class_num = 16
class HybridSN(nn.Module):
''' your code here '''
def __init__(self):
super(HybridSN, self).__init__()
self.conv3d_1 = nn.Sequential(
nn.Conv3d(1, 8, kernel_size=(7, 3, 3), stride=1, padding=0),
nn.BatchNorm3d(8),
nn.ReLU(inplace = True),
)
self.conv3d_2 = nn.Sequential(
nn.Conv3d(8, 16, kernel_size=(5, 3, 3), stride=1, padding=0),
nn.BatchNorm3d(16),
nn.ReLU(inplace = True),
)
self.conv3d_3 = nn.Sequential(
nn.Conv3d(16, 32, kernel_size=(3, 3, 3), stride=1, padding=0),
nn.BatchNorm3d(32),
nn.ReLU(inplace = True)
)
self.conv2d_4 = nn.Sequential(
nn.Conv2d(576, 64, kernel_size=(3, 3), stride=1, padding=0),
nn.BatchNorm2d(64),
nn.ReLU(inplace = True),
)
self.fc1 = nn.Linear(18496,256)
self.fc2 = nn.Linear(256,128)
self.fc3 = nn.Linear(128,16)
self.dropout = nn.Dropout(p = 0.4)
def forward(self,x):
out = self.conv3d_1(x)
out = self.conv3d_2(out)
out = self.conv3d_3(out)
out = self.conv2d_4(out.reshape(out.shape[0],-1,19,19))
out = out.reshape(out.shape[0],-1)
out = F.relu(self.dropout(self.fc1(out)))
out = F.relu(self.dropout(self.fc2(out)))
out = self.fc3(out)
return out
# 随机输入，测试网络结构是否通
# x = torch.randn(1, 1, 30, 25, 25)
# net = HybridSN()
# y = net(x)
# print(y.shape)
创建数据集
首先对高光谱数据实施PCA降维；然后创建 keras 方便处理的数据格式；然后随机抽取 10% 数据做为训练集，剩余的做为测试集。
首先定义基本函数：
# 对高光谱数据 X 应用 PCA 变换
def applyPCA(X, numComponents):
newX = np.reshape(X, (-1, X.shape[2]))
pca = PCA(n_components=numComponents, whiten=True)
newX = pca.fit_transform(newX)
newX = np.reshape(newX, (X.shape[0], X.shape[1], numComponents))
return newX
# 对单个像素周围提取 patch 时，边缘像素就无法取了，因此，给这部分像素进行 padding 操作
def padWithZeros(X, margin=2):
newX = np.zeros((X.shape[0] + 2 * margin, X.shape[1] + 2* margin, X.shape[2]))
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 ，然后创建成符合 keras 处理的格式
def createImageCubes(X, y, windowSize=5, removeZeroLabels = True):
# 给 X 做 padding
margin = int((windowSize - 1) / 2)
zeroPaddedX = padWithZeros(X, margin=margin)
# split patches
patchesData = np.zeros((X.shape[0] * X.shape[1], windowSize, windowSize, X.shape[2]))
patchesLabels = np.zeros((X.shape[0] * X.shape[1]))
patchIndex = 0
for r in range(margin, zeroPaddedX.shape[0] - margin):
for c in range(margin, zeroPaddedX.shape[1] - margin):
patch = zeroPaddedX[r - margin:r + margin + 1, c - margin:c + margin + 1]
patchesData[patchIndex, :, :, :] = patch
patchesLabels[patchIndex] = y[r-margin, c-margin]
patchIndex = patchIndex + 1
if removeZeroLabels:
patchesData = patchesData[patchesLabels>0,:,:,:]
patchesLabels = patchesLabels[patchesLabels>0]
patchesLabels -= 1
return patchesData, patchesLabels
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
下面读取并创建数据集：
# 地物类别
class_num = 16
X = sio.loadmat('Indian_pines_corrected.mat')['indian_pines_corrected']
y = sio.loadmat('Indian_pines_gt.mat')['indian_pines_gt']
# 用于测试样本的比例
test_ratio = 0.90
# 每个像素周围提取 patch 的尺寸
patch_size = 25
# 使用 PCA 降维，得到主成分的数量
pca_components = 30
print('Hyperspectral data shape: ', X.shape)
print('Label shape: ', y.shape)
print('\n... ... PCA tranformation ... ...')
X_pca = applyPCA(X, numComponents=pca_components)
print('Data shape after PCA: ', X_pca.shape)
print('\n... ... create data cubes ... ...')
X_pca, y = createImageCubes(X_pca, y, windowSize=patch_size)
print('Data cube X shape: ', X_pca.shape)
print('Data cube y shape: ', y.shape)
print('\n... ... create train & test data ... ...')
Xtrain, Xtest, ytrain, ytest = splitTrainTestSet(X_pca, y, test_ratio)
print('Xtrain shape: ', Xtrain.shape)
print('Xtest shape: ', Xtest.shape)
# 改变 Xtrain, Ytrain 的形状，以符合 keras 的要求
Xtrain = Xtrain.reshape(-1, patch_size, patch_size, pca_components, 1)
Xtest = Xtest.reshape(-1, patch_size, patch_size, pca_components, 1)
print('before transpose: Xtrain shape: ', Xtrain.shape)
print('before transpose: Xtest shape: ', Xtest.shape)
# 为了适应 pytorch 结构，数据要做 transpose
Xtrain = Xtrain.transpose(0, 4, 3, 1, 2)
Xtest = Xtest.transpose(0, 4, 3, 1, 2)
print('after transpose: Xtrain shape: ', Xtrain.shape)
print('after transpose: Xtest shape: ', Xtest.shape)
""" 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)
开始训练
这里使用cpu训练，只跑了50个epoch
# 使用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
for epoch in range(50):
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')
模型测试
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)
备用函数¶
下面是用于计算各个类准确率，显示结果的备用函数，以供参考
from operator import truediv
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))
下面代码用于显示分类结果：
# load the original image
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, numComponents= pca_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 ... ...')