高光谱图像数据？——What & How

文摘 2024-11-09 21:20 江苏

什么是高光谱图像

光谱分辨率在10-2λ数量级范围内的光谱图像称为高光谱图像(Hyperspectral Image)。

高光谱遥感的发展得益于成像光谱技术的发展与成熟。成像光谱技术最大特点是将成像技术与光谱探测技术结合，在对目标的空间特征成像的同时，对每个空间像元经过色散形成几十个乃至几百个窄波段以进行连续的光谱覆盖。

这样形成的数据可以用“三维数据块”来形象地描述。其中x和y表示二维平面像素信息坐标轴，第三维（λ轴）是波长信息坐标轴。

高光谱图像集样本的图像信息与光谱信息于一身。图像信息可以反映样本的大小、形状、缺陷等外部品质特征，由于不同成分对光谱吸收也不同，在某个特定波长下图像对某个缺陷会有较显著的反映，而光谱信息能充分反映样品内部的物理结构、化学成分的差异。

高光谱图像数据

遥感影像具有四种分辨率：时间分辨率、空间分辨率、辐射分辨率和光谱分辨率。

其中，对于光谱分辨率而言，按照成像传感器波谱通道划分数目的多少，分为多光谱、高光谱和超光谱。

通常认为，多光谱波段数目在100以下，高光谱波段数目在100~10000之间，超光谱波段数目在10000以上。

高光谱图像就是好多个灰度图像叠加到一起，每一个灰度图代表了一个光谱波段。

举个例子来说，一般的二维图像表示，这个图像像素为255*255，也就是说这幅图像有255*255个像素，可以是灰度图或者彩色图，相关信息包含在每一个像素中。

而同样大小的高光谱图像，若包含200个光谱波段的信息，那就该表示为255*255*200，可以理解为每一个像素上有200维的光谱域信息，就类似于，200幅255*255的二维图像叠加在一起，200幅图像中相同位置像素的灰度值画成曲线表示出来便是这一像素点的光谱域信息了。

也就是说，高光谱图像不仅包含丰富的光谱域信息，同时也跟一般的二维图像一样，包含相同的空间域信息。

Hyperspectral Remote Sensing Scenes - Grupo de Inteligencia Computacional (GIC) (ehu.eus)(点击下方阅读原文可跳转)

包含Indian Pines、Salinas、Pavia Centre and University、Cuprite、 Kennedy Space Center、Botswana、anomaly detection（7个地区），在每个地区都会有介绍类型数量、波段数，可直接下载。

这里，只下载了Botswana的数据，如下图所示：

高光谱图像数据的读取

1.Matlab读取.Mat格式的高光谱遥感数据：

load('Botswana.mat')InputMatImg=Botswana;b = size(InputMatImg);fprintf('输入图像宽度为 %d\n',b(1));fprintf('输入图像高度为 %d\n',b(2));fprintf('输入图像波段数为 %d\n',b(3));i=120;j=180;k=220;%自选三个波段InputImg_r= InputMatImg(:,:,i);%获取第i个波段InputImg_g= InputMatImg(:,:,j);%获取第j个波段InputImg_b= InputMatImg(:,:,k);%获取第k个波段InputImg_r= uint8(InputImg_r);%将i波段的灰度值转为0~255InputImg_g= uint8(InputImg_g);%将j波段的灰度值转为0~255InputImg_b= uint8(InputImg_b);%将k波段的灰度值转为0~255RGBImg=cat(3,InputImg_r,InputImg_g,InputImg_b);%将i、j、k三个波段进行合成figure;subplot(221);imshow(InputImg_r);title('红色波段');subplot(222);imshow(InputImg_g);title('绿色波段');subplot(223);imshow(InputImg_b);title('蓝色波段');subplot(224);imshow(RGBImg);title('合成波段');imwrite(InputImg_r,['MATBand',num2str(i),'.jpg']);imwrite(InputImg_g,['MATBand',num2str(j),'.jpg']);imwrite(InputImg_b,['MATBand',num2str(k),'.jpg']);imwrite(RGBImg,'compositeMATRGBimg.jpg');

2.python读取.Mat格式的高光谱遥感数据：

.mat是数据文件。mat数据文件是matlab的数据存储标准格式。mat数据文件是标准的二进制文件.

scipy是构建在numpy的基础之上的，它提供了许多的操作numpy的数组的函数。

https://scipy.io/(点击下方阅读原文可跳转)包提供了多种功能来解决不同格式的文件的输入和输出。

使用loadmat()函数加载.mat数据文件：

def testLoadmat(is_plot):    """
    作用：        loadmat()函数用于加载.mat数据文件
    参数:        file_name: str        .mat文件的文件名（当appendmat = True时，不用加.mat后缀），也可以打开类似文件的对象（file-like object）
        mdict: dict, 可选，插入文件变量的字典
        appendmat: bool , 可选，若为真则在文件名后添加后缀
        byte_order: str / None, 可选        默认为none。表示从mat文件中猜测的字节顺序，可以是“native”，“=”，“little”，“<”，“>”，“BIG”之一。
        mat_dtype: bool，可选        若为真，则返回与加载到MATLAB中相同的dtype的数组，而不是保存数组的dtype。
        squeeze_me: bool，可选，判断是否压缩单位矩阵的维数
        chars_as_string: bool，可选，将char数组转string数组
        matlab_compatible: bool，可选        返回被MATLAB读取的数组（当squeeze_me = false, chars_as_string = false, mat_dtype = true, struct_as_record = true）
        struct_as_record: bool，可选        设置flag来判断加载MATLAB以numpy记录数组还是以原形式numpy数组（dtype为对象）。
        verify_compressed_data_integrity: bool，可选，MATLAB文件的长度是否已确认。
        variable_names: None / sequence        若为None（默认），则读取文件中的所有变量。否则variable_names将是一个string序列，表示需要从mat文件中读取的变量名。此时读取器将跳过不是这个变量名的序列，一定程度上减少读取时间。
    返回:        mat_dict: dict        以变量名为key，数组为values的字典
# 加载数据X = sio.loadmat('data/Botswana/Botswana.mat')  print(X)此时输出的X是一个字典：{'__header__': b'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Thu Feb 20 15:12:42 2014', '__version__': '1.0', '__globals__': [], 'Botswana': array([[[3996, 3952, 3698, ...,   60,   53,   54], [8471, 8255, 8288, ...,    0,    0,    0]]], dtype=uint16)}若要获取其中的数据部分Botswana，则可直接输出X[‘Botswana’]，

Python代码如下：

import scipy.io as sio
# 加载高光谱数据和其对应的标签X = sio.loadmat('data/Botswana/Botswana.mat')['Botswana']  # 取字典中该项的值y = sio.loadmat('data/Botswana/Botswana_gt.mat')['Botswana_gt']
# 打印高光谱数据及其标签的形状print('Hyperspectral data shape: ', X.shape)  # (1476, 256, 145)print('Label shape: ', y.shape)  # (1476, 256)

高光谱图像数据的预处理及使用（python）

PCA（Principal Component Analysis）是一种常见的数据分析方式，即主成分分析技术，又称主分量分析，常用于高维数据的降维，把多指标转化为少数几个综合指标,可用于提取数据的主要特征分量。

在统计学中，主成分分析PCA是一种简化数据集的技术。它是一个线性变换。这个变换把数据变换到一个新的坐标系统中，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上，依次类推。主成分分析经常用于减少数据集的维数，同时保持数据集的对方差贡献最大的特征。这是通过保留低阶主成分，忽略高阶主成分做到的。这样低阶成分往往能够保留住数据的最重要方面。

关于 PCA降维的数学原理可参考这篇文章，非常详细。【机器学习】降维——PCA（非常详细） - 知乎 (zhihu.com)(点击下方阅读原文可跳转)

代码如下：

def applyPCA(X, numComponents):    newX = np.reshape(X, (-1, X.shape[2]))  # 行数未知，列数为145 将高光谱数据展开成一维数据序列(377856, 145)    # print(newX.shape)    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
# 对高光谱数据 X 应用 PCA 变换X_pca = applyPCA(X, numComponents=pca_components)  print('Data shape after PCA: ', X_pca.shape)

sklearn.decomposition.PCA(n_components=None, copy=True, whiten=False)

参数：

n_components:

意义：PCA算法中所要保留的主成分个数n，也即保留下来的特征个数n

类型：int 或者 string，缺省时默认为None，所有成分被保留。

copy:

类型：bool，True或者False，缺省时默认为True。

意义：表示是否在运行算法时，将原始训练数据复制一份。

若为True，则运行PCA算法后，原始训练数据的值不会有任何改变，因为是在原始数据的副本上进行运算；

若为False，则运行PCA算法后，原始训练数据的值会改，因为是在原始数据上进行降维计算。

whiten:

类型：bool，缺省时默认为False

意义：白化，使得每个特征具有相同的方差。

PCA方法：fit_transform(X)
对部分数据先拟合fit，找到该part的整体指标，如均值、方差、最大值最小值等等，然后对该X进行转换transform，从而实现数据的标准化、归一化等等。
用X来训练PCA模型，同时返回降维后的数据。
newX=pca.fit_transform(X)，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
# 在每个像素周围提取 patch ，然后创建成符合 keras 处理的格式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

对单个像素周围提取 patch 时，边缘像素就无法取了，因此，给边缘像素进行 padding 操作。在原有数据的基础上给周围增加了7列0元素。

同时，对原有数据进行取样，取样大小为15*15.在一个特征点上，共取出了1476*256=377856个样本，最后将其中标签为0的去掉，共得到3248个样本，数据块大小为(3248, 15, 15, 30)。

创建训练数据和测试数据集:

def splitTrainTestSet(X, y, testRatio, randomState=345):    # 能够将数据集按照用户的需要指定划分为训练集和测试集    # test_size    若在0~1之间，为测试集样本数目与原始样本数目之比；若为整数，则是测试集样本的数目    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
Xtrain, Xtest, ytrain, ytest = splitTrainTestSet(X_pca, y, test_ratio)print('Xtrain shape: ', Xtrain.shape)print('Xtest  shape: ', Xtest.shape)

train_test_split随机划分训练集和测试集：

X_train,X_test, y_train, y_test =train_test_split(

train_data,train_target,test_size=0.25, random_state=0,stratify=y)

# train_data：所要划分的样本特征集

# train_target：所要划分的样本结果

# test_size：样本占比，如果是整数的话就是样本的数量

# random_state：是随机数的种子,随机数种子：其实就是该组随机数的编号，在需要重复试验的时候，保证得到一组一样的随机数。比如你每次都填1，其他参数一样的情况下你得到的随机数组是一样的。但填0或不填，每次都会不一样。

# stratify: 是为了保持split前类的分布，通常在这种类分布不平衡的情况下会用到stratify。比如有100个数据，80个属于A类，20个属于B类。如果train_test_split(... test_size=0.25, stratify = y), 那么split之后数据如下：
training: 75个数据，其中60个属于A类，15个属于B类。
testing: 25个数据，其中20个属于A类，5个属于B类。

用了stratify参数，training集和testing集的类的比例是 A：B= 4：1，等同于split前的比例（80：20）。

将stratify=X就是按照X中的比例分配，将stratify=y就是按照y中的比例分配。

改变 Xtrain, Ytrain 的形状，以符合 keras 的要求：

Xtrain = Xtrain.reshape(-1, patch_size, patch_size, pca_components, 1)# (974, 15, 15, 30, 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 结构，数据要做转置：

# 一个数的下标为[z, y, x]Xtrain = Xtrain.transpose(0, 4, 3, 1, 2)  # 对数据做转置Xtest = Xtest.transpose(0, 4, 3, 1, 2)print('after transpose: Xtrain shape: ', Xtrain.shape)  # (162, 1, 30, 15, 15)print('after transpose: Xtest  shape: ', Xtest.shape)

因为输入数据后，将首先通过三维卷积层，所以要将数据调整为三维卷积层可以接受的形式，

即输入层现在的维度是（Channel×Depth1×Height1×Width1）。

对Xtrain、Xtest数据做调整：Channel-1、Depth1-30、Height1-15、Width1-15

制作数据集：

""" 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 和 testloadertrainset = TrainDS()testset = TestDS()

train_loader = torch.utils.data.DataLoader(dataset=trainset,                                           batch_size=batch_size,  # 256                                           shuffle=True,                                           num_workers=0)
test_loader = torch.utils.data.DataLoader(dataset=testset,                                          batch_size=batch_size,                                          shuffle=False,                                          num_workers=0)

torch.FloatTensor()的作用就是把给定的list或者numpy转换成浮点数类型的tensor。

torch.LongTensor（）的作用就是把给定的list或者numpy转换成64位整型的tensor。

# torch.utils.data.DataLoader数据加载器
# 该接口主要用来将自定义的数据读取接口的输出或者PyTorch已有的数据读取接口的输入按照batch size封装成Tensor，后续只需要再包装成Variable即可作为模型的输入，因此该接口有点承上启下的作用，
# dataset (Dataset) – 加载数据的数据集。
# batch_size (int, optional) – 每个batch加载多少个样本(默认: 1)。
# shuffle (bool, optional) – 设置为True时会在每个epoch重新打乱数据(默认: False).
# sampler (Sampler, optional) – 定义从数据集中提取样本的策略，即生成index的方式，可以顺序也可以乱序
# num_workers (int, optional) – 用多少个子进程加载数据。0表示数据将在主进程中加载(默认: 0)
# collate_fn (callable, optional) –将一个batch的数据和标签进行合并操作。
# pin_memory (bool, optional) –设置pin_memory=True，则意味着生成的Tensor数据最开始是属于内存中的锁页内存，这样将内存的Tensor转义到GPU的显存就会更快一些。
# drop_last (bool, optional) – 如果数据集大小不能被batch size整除，则设置为True后可删除最后一个不完整的batch。如果设为False并且数据集的大小不能被batch size整除，则最后一个batch将更小。(默认: False)
# timeout，是用来设置数据读取的超时时间的，但超过这个时间还没读取到数据的话就会报错。

这样对高光谱图像的数据就做好处理了，就可以输入搭建好的卷积网络模型中进行训练了。

参考资料：《HybridSN: Exploring 3-D–2-DCNN Feature Hierarchy for Hyperspectral Image Classification》

本文仅做学术分享，如有侵权，请联系删文。

http://mp.weixin.qq.com/s?__biz=MzAxNTc4MTc1Ng==&mid=2649472007&idx=2&sn=3f2edb7f510dd7ac0ff4cf0a4b99833d

新机器视觉

一个值得关注的AI视觉技术公众号，主要涉及人工智能领域机器视觉、计算机视觉、机器学习、深度学习等前沿知识干货和资源！致力于为您提供切实可行的AI学习线路。