光谱数据预处理
一. 高光谱数据读取和特征
本文中数据用例使用
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()
![图片名称](https://img2023.cnblogs.com/blog/2868879/202309/2868879-20230920170708644-1920134786.png)
全波段合成图像(直接线性合成,空间分辨率和对比度没有明显改善)
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()
![图片名称](https://img2023.cnblogs.com/blog/2868879/202309/2868879-20230921170644172-1880314556.png)
三. 图像预处理
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()
![效果一般](https://img2023.cnblogs.com/blog/2868879/202309/2868879-20230921173520784-723467612.png)