光谱数据预处理

一. 高光谱数据读取和特征

本文中数据用例使用
Pavia University HSI
以常见的公开高光谱数据mat格式为例,将mat数据转为字典类型，并读取其中的numpy

from scipy.io import loadmat
X = loadmat('data/PaviaU.mat')['paviaU']
Y = loadmat('data/PaviaU_gt.mat')['paviaU_gt'] #此为预测类型结果

该图像数据分辨率为610×340，同时有103个通道（波段）

二. 高光谱数据可视化

2.1 空间域可视化

import matplotlib.pyplot as plt
import earthpy.plot as ep
import numpy as np
data = np.moveaxis(X, 2, 0) #波段放到第一维度
ep.plot_bands(data[:10, :, :], title=[f'Band-{i}' for i in range(1, 11)], 
              cmap='gist_gray', cols=5, cbar=False, scale=True, figsize=(14, 10))
plt.show()

也可以通过遍历，随机预览十个波段的图像

fig = plt.figure(figsize = (20, 12))
ax = fig.subplots(2, 5)
for i in range(2):
  for j in range(5):
    c = np.random.randint(103)
    ax[i][j].imshow(X[:, :, c], cmap='gray')
    #plt.imshow(X[:,:,c],cmap='jet')  绘制热力图，查看各光谱具体反射情况
    ax[i][j].axis('off')
    ax[i][j].title.set_text(f"Band - {c}")
    c+=1
plt.tight_layout()
plt.show()

随机十个波段的灰度图片

随机十个波段的热力图片

通过对应三个波段合成彩色图像

ep.plot_rgb(X, rgb=(36, 17, 11), title='Composite Image of Pavia University', figsize=(10, 8))
plt.show()

全波段合成图像(直接线性合成，空间分辨率和对比度没有明显改善)

image_2d = np.sum(X,axis=2)
fig = plt.figure(figsize = (20, 12))
plt.imshow(image_2d, cmap='gray')

2.2 波段维度可视化

import matplotlib.pyplot as plt
fig = plt.figure(figsize = (10, 5))
for i in range(10):
        plt.plot(X[0][i])
plt.show()

三. 图像预处理

3.1 一般图像预处理

3.1.1 直方图均衡化

直方图均衡化通过拉伸图像的灰度范围，使得亮度级别更加均匀分布，从而增强图像的对比度
统计每个灰度值的频率，计算每个灰度值的累积概率，再将原始图像的每个灰度值映射到新的灰度值

import cv2
equalized_image = cv2.equalizeHist(image) #image可以为任意二维numpy，所以可以对高光谱任意一个波段图像进行处理

3.1.2 平滑

消除噪音的常用方式，通常有均值滤波、中值滤波、高斯滤波等
能够有效地去除椒盐噪声和其他类型的噪声，但也有可能损失细节
具体实现是使用一个滤波模板，通常是一个正方形的窗口，进行滑动并计算

blurred_image = cv2.blur(image, (5, 5))  # 使用5x5的卷积核进行均值滤波
blurred_image = cv2.GaussianBlur(image, (5, 5), 0)  # 使用5x5的高斯核进行高斯滤波
blurred_image = cv2.medianBlur(image, 5)  # 使用5x5的卷积核进行中值滤波

3.1.3 边缘检测

边缘检测可以用于图像增强，使图像中的边界更加清晰和突出
通过增强边缘，可以使图像更加鲜明，并减少噪声和不必要的细节
常用检测方法：Sobel算子、Canny、拉普拉斯算子、Roberts算子

# 使用Sobel算子计算X和Y方向的梯度
gradient_x = cv2.Sobel(image, cv2.CV_64F, 1, 0, ksize=3)
gradient_y = cv2.Sobel(image, cv2.CV_64F, 0, 1, ksize=3)
# 将X和Y方向的梯度进行合并
sobel_gradient = cv2.addWeighted(gradient_x, 0.5, gradient_y, 0.5, 0)

# 使用Canny边缘检测算法
canny_edges = cv2.Canny(image, 100, 200)

# 使用拉普拉斯算子进行边缘检测
laplacian_edges = cv2.Laplacian(image, cv2.CV_64F)

# 使用Roberts算子进行边缘检测
roberts_edges = cv2.filter2D(image, -1, np.array([[1, 0], [0, -1]]))

3.2 高光谱图像预处理（波段维度）

后面将图像展成一维，使数据变成二维数据，含义为每个像素对应的光谱
方便对样品（纵向）和光谱（横向）进行预处理
前者方便样品之间的比较，后者方便分析具体的光谱信息

data = np.reshape(X, (X.shape[0] * X.shape[1], X.shape[2])) #data为(207400, 103)的数据
#data = np.reshape(data, X.shape) 用于转回原来的形状

3.2.1 特征缩放（纵向样品）

光程变化和样品稀释等会对光谱产生一定的影响，所以需要进行特征缩放进行比较，方便样本之间的比较

from sklearn.preprocessing import MinMaxScaler, StandardScaler
#标准化
stand_data = StandardScaler().fit_transform(data)#平均值为0，标准差为1，从而使数据分布更加接近正态分布

#最大值最小值归一化
minmax_data = MinMaxScaler().fit_transform(data) #通过线性变化缩放到指定范围

#均值中心化
centered_data = data - np.mean(data, axis=0)

#正态变换
normalized_data = (data - np.mean(data, axis=0)) / np.std(data, axis=0)

3.2.2 导数（横向光谱）

注意前面的算子和卷积操作都是在空间域求导
样品不同组分的相互干扰会导致吸收光谱谱线重叠的现象
常用一阶求导和二阶求导减少光谱干扰，提高光谱的分辨率和灵敏度
同时获得如光谱峰值、波段间的距离和光谱曲线的形状的主要特征

#一阶导数，光谱特征减少一个（变成102）
first_derivative = np.diff(data, axis=1)

#二阶导数，光谱特征再减少一个（变成101）
second_derivative = np.diff(first_derivative, axis=1)

3.2.3 基线校正（横向光谱）

基线通常是指光谱中不包含样品信号的基础背景，由于光谱仪器漂移、光源变化、采样条件发生改变等原因
可能会导致光谱数据中存在趋势效应，使得数据在不同的波长或位置上呈现出逐渐增加或递减的趋势
基线校正的目标是将光谱的基线调整到一个平坦的水平，使得光谱信号在各个波长或位置上的峰值和谷值更清晰可见

#使用多项式进行拟合校正
x = np.linspace(0, data.shape[1]-1, data.shape[1])
for i in range(data.shape[0]):
    coefficients = np.polyfit(x, data[i], deg=1)  #这里是一阶，等价于线性拟合
    trend = np.polyval(coefficients, x)
    data[i] = data[i] - trend

3.2.4 多元散射校正（纵向样品和横向光谱）

由于样品的不均衡性，常常导致待测的样品光谱具有很大的差异性
多元散射矫正用来校正每个光谱的散射，以期得到较为理想的光谱，有效提高光谱的信噪比

#线性拟合，校正再除以散射因子
x = np.mean(data,axis=0) #每一列的均值
for i in range(data.shape[0]):
    coefficients = np.polyfit(x, data[i], deg=1)
    data[i] = (data[i] - coefficients[1])/coefficients[0]

3.2.5 小波变换(横向光谱)

相比傅里叶变换，小波变换是时域和频域的局部变换，能够提供更精细的频率分析
通过小波变换，我们可以将光谱信号分解成多个不同频率的小波系数，这些小波系数可以表示光谱中的特征

import pywt
wavelet = pywt.Wavelet('db8') ## 选用Daubechies8小波,更多的保留信号细节
maxlev = pywt.dwt_max_level(data.shape[1], wavelet.dec_len) #最大变换级别
for i in range(data.shape[0]):
    coefficients = pywt.wavedec(data[i], wavelet,level=maxlev) #小波变换
    threshold = 0.04
    for j in range(1, len(coefficients)): #对高频部分进行软阈值处理，消除噪音，更加平滑
        coefficients[j] = pywt.threshold(coefficients[j], threshold * max(coefficients[j]))
    data[i] = pywt.waverec(coefficients, wavelet)[:data.shape[1]] #重构信号，这里重构多了一个，截掉

3.2.6 其它

一维信号平滑处理（针对光谱信号横向处理）

from scipy import signal
#移动平均滤波
smoothed_signal = np.convolve(signal, np.ones(window_size)/window_size, mode='same')
#中值滤波
smoothed_signal = signal.medfilt(signal, kernel_size=window_size)
#Savitzky-Golay滤波
smoothed_signal = signal.savgol_filter(signal, window_size, order)

3.3 针对性预处理（以遥感图像为例）

3.3.1 传感器矫正

不同原理的采集装置在使用的时候会存在一定的误差

3.3.2 大气辐射矫正

消除大气散射、吸收、反射引起的误差
将辐射亮度值或表现反射率转化为地表实际反射

3.3.3 几何矫正

遥感图像成像由于受到卫星位置和运动状态变化的影响
以及地球自转、地球表面曲率、地形起伏、大气折射等因素
这样会造成采集到的图像发生几何畸变

3.4 聚类

根据处理的数据进行聚类（需要对数据进行较好的预处理结果才会好）

from sklearn.cluster import KMeans
cluster = KMeans(n_clusters=10, random_state=0, max_iter=30)
cluster.fit(np.concatenate((first_derivative, second_derivative), axis=1))
label = cluster.labels_  # 对原数据表进行类别标记
label = np.reshape(label,Y.shape)

绘制对比图

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
cmap = ListedColormap(['red', 'blue', 'green', 'yellow', 'purple', 'orange', 'cyan', 'magenta', 'gray', 'black'])
axs[0].imshow(Y, cmap=cmap)
axs[0].axis('off')
axs[0].set_title('real classification')

axs[1].imshow(label, cmap=cmap)
axs[1].axis('off')
axs[1].set_title('predict classification')

plt.tight_layout()
plt.show()