HybridSN 高光譜分類

S. K. Roy, G. Krishna, S. R. Dubey, B. B. Chaudhuri HybridSN: Exploring 3-D–2-D CNN Feature Hierarchy for Hyperspectral Image Classification, IEEE GRSL 2020

這篇論文構建了一個 混合網絡 解決高光譜圖像分類問題，首先用 3D卷積，然后使用 2D卷積，代碼相對簡單，下面是代碼的解析。

取得數據

! 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 類

三維卷積部分：

• 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 個節點，是最終的分類類別數。

下面是 HybridSN 類的代碼

class_num = 16

class HybridSN(nn.Module):
  def __init__(self, num_classes=16):
    super(HybridSN, self).__init__()
    # conv1：（1, 30, 25, 25）， 8個 7x3x3 的卷積核 ==>（8, 24, 23, 23）
    self.conv1 = nn.Conv3d(1, 8, (7, 3, 3))
    # conv2：（8, 24, 23, 23）， 16個 5x3x3 的卷積核 ==>（16, 20, 21, 21）
    self.conv2 = nn.Conv3d(8, 16, (5, 3, 3))
    # conv3：（16, 20, 21, 21），32個 3x3x3 的卷積核 ==>（32, 18, 19, 19）
    self.conv3 = nn.Conv3d(16, 32, (3, 3, 3))
    # conv3_2d （576, 19, 19），64個 3x3 的卷積核 ==>（（64, 17, 17）
    self.conv3_2d = nn.Conv2d(576, 64, (3,3))
    # 全連接層（256個節點）
    self.dense1 =  nn.Linear(18496,256)
    # 全連接層（128個節點）
    self.dense2 =  nn.Linear(256,128)
    # 最終輸出層(16個節點)
    self.out = nn.Linear(128, num_classes)
    #  Dropout（0.4)
    self.drop = nn.Dropout(p=0.4)
    # 這是個坑， 下面細說
    # self.soft = nn.LogSoftmax(dim=1)
    # 激活函數ReLU
    self.relu = nn.ReLU()
  def forward(self, x):
    out = self.relu(self.conv1(x))
    out = self.relu(self.conv2(out))
    out = self.relu(self.conv3(out))
    # 進行二維卷積，因此把前面的 32*18 reshape 一下，得到 （576, 19, 19）
    out = out.view(-1, out.shape[1] * out.shape[2], out.shape[3], out.shape[4])
    out = self.relu(self.conv3_2d(out))
    # flatten 操作，變為 18496 維的向量，
    out = out.view(out.size(0), -1)
    out = self.dense1(out)
    out = self.drop(out)
    out = self.dense2(out)
    out = self.drop(out)
    out = self.out(out)
    out = self.soft(out)
    return out

當我第一次寫完代碼的時候，跑完發現准確率只有84%，並且avg loss一直不怎么下降。我仔細檢查了 HybridSN類，通過增添bn層，發現效果也並不好。最終我把目標放在了唯一可能優化的地方，self.soft()。

在第一次編寫時，我是定義的softmax函數，即self.soft()=nn.softmax()，我在百度查詢softmax時，發現了還有其他的例如logsoftmax，在使用logsoftmax代替softmax之后，效果顯著提升，准確率穩定在96%以上。

那么為什么使用logsoftmax之后效果會好這么多呢？？？我也不知道，先挖個坑以后填。

經過反復測試和查閱資料找到原因了，其實不寫self.soft()就能得到正確的准確率，多寫了以后反而多此一舉。

原因是 CrossEntropyLoss() = softmax + 負對數損失（已經包含了softmax) 。如果多寫一次softmax，則結果會發生錯誤。

至於BN層的添加，我對conv3d，2d，混合都進行了添加，發現效果並不明顯。。

# 隨機輸入，測試網絡結構是否通
x = torch.randn(1, 1, 30, 25, 25)
net = HybridSN()
y = net(x)
print(y.shape)
'''
torch.Size([1, 16])
'''

創建數據集

首先對高光譜數據實施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)

模型訓練

# 使用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.0005)

# 開始訓練
total_loss = 0
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')

模型測試

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)