简介：高光谱数据预处理的Python代码与文档说明压缩包，面向课程设计、期末大作业和毕业设计等场景，也适合刚接触高光谱分析的新手，核心解决数据清洗、波段选择与降维等预处理环节的实现问题。整个包共16个文件，包含预处理和演示两个Python脚本、一份Markdown说明文档、一个CSV光谱样本数据，以及12张结果示意图，压缩后仅2.48MB，结构清晰、方便部署。项目经过严格调试并配有逐段注释，曾作为高分作业获得导师认可，目前已有533人浏览学习。使用时可借助示例数据与截图对照查看每一步输出，快速理解代码逻辑；Markdown文档还对数据格式、运行流程与常见参数作了说明，便于在此基础上替换数据集或扩展新的预处理算法，非常适合课程设计与毕业设计的基础参考。
1. 高光谱数据预处理：一个让期末大作业从“能跑”到“能用”的 Python 代码包
高光谱数据预处理在遥感、农业和食品检测方向的课程设计里，永远是那个看起来最简单、实际翻车率最高的环节——数据明明能读出来、能画图，进了模型精度就崩。这份资源是一个完整的高光谱预处理 Python 代码包，配套文档说明，覆盖从 ENVI/Mat 原始数据读取到特征矩阵构建的整条链路：坏波段剔除、Savitzky-Golay 平滑、三种归一化、PCA 降维，每个步骤都有现成函数和参数注释。适合正在做期末大作业的学生，也适合想快速复现一套预处理流程的工程师。核心思路只有一句话：把原始 DN 值变成模型真正能用的光谱特征，而不是给人看图用的漂亮曲线。
2. 环境搭建与数据读取：从三维数据立方体到二维光谱矩阵
2.1 为什么预处理要先解决数据读取顺序
高光谱数据不是普通彩色图片，它是一个三维数据立方体，三个维度分别是高度、宽度和波段数。以公开数据集为例，Indian Pines 是 145×145×200，每个像素点上有 200 个波段的值；Salinas 是 512×217×204。这种三维结构人看着直观，但绝大多数机器学习模型只吃二维输入——样本数×特征数。预处理的第一个动作，就是把三维立方体压成二维矩阵，让每一行代表一个像元的光谱曲线，每一列代表一个波段。
这个 reshape 操作看似简单，却是后面所有步骤的地基。如果行和列的顺序没理清，后面做坏波段剔除时索引就会错位，直接导致特征和标签对不上。我一般会先打印出原始数组的 shape，再确认波段维度在最后一个轴，然后再做变形。这个习惯帮我避免了至少三次数据错位的低级事故。
2.2 Python 环境安装与关键库版本选择
这个代码包基于 Python 3.8 及以上版本，核心依赖是 numpy、scipy、matplotlib 和 spectral。spectral 库用于读 ENVI 格式，如果你拿到的数据是 Mat 格式，则用 scipy.io 读取，不需要 spectral。先确认自己机器上的 Python 环境，没装的话先去 Python 官网下个 3.8 以上的版本，安装时记得勾选“Add Python to PATH”。
# 建议直接用 pip 安装，豆瓣源速度更快 pip install numpy scipy matplotlib -i https://pypi.douban.com/simple pip install spectral逻辑说明：numpy 负责数组变形和数学运算，scipy 提供 savgol_filter 平滑函数，matplotlib 用在预处理前后对比图，spectral 专读 ENVI 格式的 .hdr/.dat 文件。参数说明：spectral 这个库比较老，Python 3.10 以上偶尔会编译报错，装不上就用 gdal 或者手动解析 .hdr 头文件读取二进制数据，这个坑我踩过，后面会细说。
2.3 把 ENVI 格式读成光谱矩阵的完整代码
ENVI 格式在高光谱数据里太常见了，几乎每个公开数据集都有这个格式的影子。它的文件结构是 .hdr 头文件加上 .dat（或 .raw、.bsq）二进制主体，头文件里记录了波段数、行数、列数、数据类型、波长等信息。下面的代码演示了如何完整读取它并转换成二维矩阵。
# 读取ENVI格式高光谱影像并转换成二维光谱矩阵 import numpy as np from spectral.io import envi img = envi.open('envi_image.hdr', 'envi_image.dat') array = img.asarray() print('原始三维形状:', array.shape) # 输出类似 (145, 145, 200) rows, cols, bands = array.shape spectra_2d = array.reshape(rows * cols, bands) print('光谱矩阵形状:', spectra_2d.shape) # 输出类似 (21025, 200)逻辑说明：envi.open 需要把 .hdr 和 .dat 文件放在同一目录下，它会自动读取元数据中记录的行列波段信息。asarray 方法把数据加载成 numpy 数组，reshape 将前两个维度展平，保留波段维度。参数说明：如果你的数据是 BIP 或 BIL 格式（波段按像素或按行存储），asarray 返回的数组维度顺序可能不同，先打印 shape 再 reshape，别盲拆。
如果你拿到的是 .mat 文件，读取方式完全换一套：
# 读取Matlab格式的高光谱数据 from scipy.io import loadmat mat = loadmat('hyperspectral_data.mat') keys = [k for k in mat.keys() if not k.startswith('__')] print('文件包含的变量:', keys) # 假设变量名叫 data，形状可能是 (rows, cols, bands) 或 (bands, rows, cols) data = mat['data'] # 如果波段在第一维，需要转置成 (rows, cols, bands) if data.shape[0] < data.shape[2]: data = np.transpose(data, (1, 2, 0)) # 再展平成二维矩阵 rows, cols, bands = data.shape spectra_2d = data.reshape(rows * cols, bands)逻辑说明：loadmat 把 .mat 文件读成字典结构，变量名需要自己确认，很多数据集变量名不是 data 而是 raw_data 或 image。转置那一步是关键，Matlab 里习惯把波段放第一维，Python 里则要放在最后一维。参数说明：判断条件用维数大小去猜顺序是一个土办法，更保险的是查看数据集的说明文档，确认变量名和维度顺序。
3. 核心预处理操作：坏波段剔除、SG平滑与三种归一化方式的选型
3.1 坏波段剔除：水汽吸收带与噪声波段的区间设计
高光谱数据的前后波段通常信噪比很低，接近传感器响应极限的位置噪声非常大。此外，大气中水汽的吸收会形成几个典型吸收带，例如 1350~1420nm 和 1800~1950nm 附近的波段反射率接近零，几乎没有信息量。这些坏波段如果混进特征矩阵，会像黑匣子一样干扰后续分类器的权重计算。
剔除的方式有两种：直接按索引删列，或者将这些波段的值设为 NaN。直接删列更干净，但要注意保留波段连续性，否则后面 SG 平滑时窗口跨过断点会产生假峰；设 NaN 的做法在需要保持波段数量固定的场景下更稳妥。下面的函数按索引做删除。
# 按波段索引剔除坏波段 import numpy as np def remove_bad_bands(spectra_2d, bad_bands_indices): """ spectra_2d: 二维光谱矩阵，形状 (样本数, 波段数) bad_bands_indices: 要剔除的波段索引列表 """ mask = np.ones(spectra_2d.shape[1], dtype=bool) mask[bad_bands_indices] = False return spectra_2d[:, mask]逻辑说明：先构建一个全 True 的布尔掩膜，把要剔除的位置置为 False，然后用布尔索引过滤列。这样比直接 np.delete 更直观，而且后续要调整保留区间时，只需要改 bad_bands_indices 列表。参数说明：建议先画出每个波段的均值曲线，目视确认噪声区间，再确定索引范围。比如 200 个波段的数据，索引 0~4 和 180~199 通常直接删，具体范围看传感器型号。
3.2 Savitzky-Golay 平滑：窗口宽度与多项式阶次的搭配经验
光谱曲线里的高频噪声会掩盖真实的吸收特征，Savitzky-Golay 滤波是目前预处理里最常用的平滑方法。它本质上是一个移动窗口内的最小二乘多项式拟合，在去噪的同时能保留信号的峰谷形状——这比普通移动平均要强得多，不会把吸收峰抹平。
from scipy.signal import savgol_filter def sg_smooth(spectra_2d, window_length=15, polyorder=3): """ 对每个像元的光谱曲线做SG平滑 window_length: 窗口宽度，必须是奇数 polyorder: 多项式阶次，必须小于 window_length """ smooth_data = np.apply_along_axis( lambda x: savgol_filter(x, window_length, polyorder), axis=1, arr=spectra_2d ) return smooth_data逻辑说明：apply_along_axis 沿着行方向逐条处理光谱曲线，等价于写一个 for 循环但更简洁。savgol_filter 的 window_length 控制着平滑强度，值越大曲线越平滑，但窗口过大时原本的窄吸收峰会失真。参数选取经验：窗口宽度建议从 9 起步，按波段总数的 5%~10% 来试；polyorder 用 2 或 3 就够，超过 4 阶容易产生振荡伪影。窗口 15 阶次 3 是光谱预处理论文里的标配，但拿到自己的数据后先画两张原始和滤波后的对比图，确认没有把吸收峰压扁再批量跑。
3.3 三种归一化方式对比：最大最小与 Z-Score 还是矢量归一化
不同传感器、不同批次数据的反射率基线差异很大，同一地物在不同光照条件下的光谱值会整体平移或缩放。归一化的目标就是消除这些非化学成分的干扰，让模型学到的差异来自物质本身。
# 三种归一化方式的实现与对比 import numpy as np from sklearn.preprocessing import MinMaxScaler, StandardScaler def normalize_spectra(spectra_2d, method='minmax'): """ method: minmax -> [0,1]区间缩放 zscore -> 均值0标准差1 vector -> 矢量归一化，每个样本模长为1 """ if method == 'minmax': # 按行做最大最小缩放 min_vals = spectra_2d.min(axis=1, keepdims=True) max_vals = spectra_2d.max(axis=1, keepdims=True) return (spectra_2d - min_vals) / (max_vals - min_vals + 1e-8) elif method == 'zscore': # 标准分数归一化 mean_vals = spectra_2d.mean(axis=1, keepdims=True) std_vals = spectra_2d.std(axis=1, keepdims=True) return (spectra_2d - mean_vals) / (std_vals + 1e-8) elif method == 'vector': # 矢量归一化：每个样本除以自身的L2范数 norms = np.linalg.norm(spectra_2d, axis=1, keepdims=True) return spectra_2d / (norms + 1e-8)逻辑说明：代码里加 1e-8 是为了防止除零，这在光谱数据里很常见——某个坏波段如果被剔除后残留常数区域，std 或 norm 可能为 0。minmax 会放大噪声对动态范围的影响；zscore 在数据接近正态分布时最稳定；vector 归一化在光谱库匹配场景比较常用。选型建议：做分类实验时我最常用 zscore，因为大多数分类器对特征的尺度不敏感；minmax 偶尔会导致模型在个别噪声波段上过拟合。如果后续接 PCA，则先 zscore 再 PCA 是标准流程。
4. 光谱变换与降维：把冗余波段压缩成有效特征
4.1 导数光谱：一阶导数与二阶导数的差异增强方法
光谱的导数变换可以放大谱形细节，一阶导数反映反射率随波长的变化率，能减弱基线漂移的影响；二阶导数对窄带吸收特征特别敏感，但也会放大高频噪声。做导数变换前必须先平滑，否则噪声会被导数额外放大，翻车概率极高。
# 对平滑后的光谱计算一阶和二阶导数 import numpy as np def derivative_spectra(smooth_data, order=1, wavelen=None): """ smooth_data: SG平滑后的光谱矩阵 order: 1表示一阶导数，2表示二阶导数 wavelen: 可选的波长列表，默认为均匀间隔 """ if wavelen is None: # 默认均匀间隔，直接用np.gradient deriv = np.gradient(smooth_data, axis=1) else: # 波长不均匀时传入波长数组，结果按物理单位计算 deriv = np.gradient(smooth_data, wavelen, axis=1) if order == 2: # 二阶导数：对一阶结果再求一次梯度 if wavelen is None: deriv = np.gradient(deriv, axis=1) else: deriv = np.gradient(deriv, wavelen, axis=1) return deriv逻辑说明：np.gradient 用中心差分算法计算梯度，边界点自动退化到单边差分。传入波长数组时得到的是 dR/dλ，单位是反射率每纳米；不传则按波段间隔 1 计算。参数说明：如果做完导数后曲线像锯齿一样毛刺密布，说明先前的 SG 平滑窗口设小了，回到 3.2 把窗口加大到 21 或 25 再跑。很多论文里先做一阶导数再归一化，这个顺序会明显改变特征分布，值得尝试。
4.2 PCA 降维：主成分数量的判定标准
高光谱的波段相关性极高，相邻波段之间相关系数经常超过 0.9。200 个波段直接进模型，不仅计算量大，而且多重共线性会让分类器不稳定。PCA 通过正交变换把原始特征压缩成互不相关的几个主成分，前 10~20 个主成分往往能解释 99% 以上的方差。
# PCA降维并查看方差解释率 from sklearn.decomposition import PCA import numpy as np def pca_reduce(spectra_2d, variance_threshold=0.99, n_components=None): """ variance_threshold: 累计方差解释率阈值，达到即截断 n_components: 指定主成分数量，不传则按阈值自动选 """ pca = PCA(n_components=n_components) pca_result = pca.fit_transform(spectra_2d) if n_components is None: # 按累计方差解释率自动筛选成分数 cumsum_ratio = np.cumsum(pca.explained_variance_ratio_) n_selected = int(np.argmax(cumsum_ratio >= variance_threshold) + 1) pca = PCA(n_components=n_selected) pca_result = pca.fit_transform(spectra_2d) print('自动选择的主成分数量:', n_selected) print('累计方差解释率:', cumsum_ratio[n_selected - 1]) return pca_result, pca逻辑说明：PCA 先拟合再变换，explained_variance_ratio_ 给出每个主成分的方差解释比例。阈值设 0.99 是常见选择，但实际实验中主成分取 20~30 个时分类精度往往达到峰值，取太多反而过拟合。参数说明：PCA 对输入特征做了中心化，所以前面归一化推荐 zscore，效果比 minmax 好。如果跑完 PCA 后发现前几个主成分的可视化分布混在一起，需要结合地物标签着色查看是否真的可分，别只盯着累计贡献率。
4.3 训练集与测试集划分：按光谱相似性避免数据泄漏
预处理做完后，数据要以样本为单位划分训练集和测试集。这里一个常见误区是：如果原始数据是影像，直接从二维矩阵里随机抽样划分，会导致同一个地块的相邻像元同时出现在训练集和测试集里，精度虚高。更严谨的做法是按空间位置划分——先按行或按块切分，再展平成样本。
# 按空间区域划分训练集和测试集 from sklearn.model_selection import train_test_split def spatial_split(spectra_3d, labels_3d, test_size=0.3, random_state=42): """ spectra_3d: 原始三维影像，shape=(rows, cols, bands) labels_3d: 标签三维图，shape=(rows, cols) """ rows, cols, bands = spectra_3d.shape # 把三维数据展平成二维 X = spectra_3d.reshape(rows * cols, bands) y = labels_3d.reshape(rows * cols) # 过滤掉标签为0的区域（通常是背景） valid_mask = y != 0 X_valid, y_valid = X[valid_mask], y[valid_mask] # 常规随机划分作为对比 X_train, X_test, y_train, y_test = train_test_split( X_valid, y_valid, test_size=test_size, stratify=y_valid, random_state=random_state ) return X_train, X_test, y_train, y_test逻辑说明：stratify 参数让划分后各类别比例与全量数据保持一致，这在地物类别不平衡时非常重要。如果做空间划分，则应该引入空间坐标信息，将影像按矩形块切分而不是逐像元随机抽。参数说明：random_state 固定为 42 是保证实验可复现的底线——没有固定随机种子的话，每次跑出来的精度都不同，老师验收时你没法解释波动，这是血泪教训。
5. 避坑：高光谱预处理内存崩溃与波段匹配的常见问题排查
5.1 内存溢出：reshape 之后程序直接卡死
现象：代码跑到一个 512×512×200 的影像做 reshape 时，内存占用飙到十几个 GB，程序被系统杀死。
原因：高光谱影像展开后的样本数太大。512×512 的影像有 26 万个像元，每个像元 200 维 float64，数据量超过 400MB，这还只是单张影像；如果堆叠了多个文件，内存翻倍。reshape 本身不复制数据，但后续归一化、SG 平滑都会产生中间数组，内存峰值可能达到原始数据的三倍。
解决：先转成 float32 再计算，float64 到 float32 精度损失对光谱分类影响极小，但内存立刻减半。此外，用小批量处理替代一次性全量计算——SG 平滑可以按行循环处理，不把整个矩阵读进内存，或者用 np.memmap 做内存映射。
# 转float32降低内存占用 import numpy as np array = array.astype(np.float32)5.2 归一化之后分类精度反而下降
现象：原始数据进 SVM 精度有 82%，做了 minmax 归一化后掉到 75%。
原因：高光谱数据里的噪声波段经过 minmax 归一会把噪声的动态范围放大到和有用信号相同量级，分类器被迫学习了噪声模式。另外，如果原始数据已经做过传感器定标和大气校正，不同像元间的值本身具有物理可比性，强行归一化反而抹掉了差异。
解决：先尝试不做归一化，直接做标准化（zscore）。zscore 对离群值容忍度更高，不容易被单个噪声波段带偏。如果 zscore 也没提升，干脆跳过归一化，只做坏波段剔除和 SG 平滑后直接进模型。预处理不是每步都要上，按结果检验。
5.3 波段匹配错位：训练和测试的波段对不上
现象：同一数据集，第一次跑训练准确率 98%，第二次用不同的预处理脚本跑，准确率骤降到 60%。
原因：两次处理时坏波段剔除的索引列表不一致，或者 ENVI 文件读取时波段顺序不同。比如数据在转置后波段维在第一维还是最后一维没统一，造成训练和测试特征列错位。
解决：每次预处理结束后保存一个波段索引文件（npy 格式），带着波形图确认关键吸收带位置。我在预处理流程最后强制加一段检查：把处理后的光谱矩阵的 shape 和保留波段数写入日志，跟原始 head 文件里记录的波段数做减法验证。
# 保存和校验保留波段索引 np.save('bad_bands_removed.npy', bad_bands_indices) np.save('selected_bands.npy', np.arange(bands)[mask])5.4 读取 ENVI 文件报错 segment fault
现象：执行 envi.open 时程序直接崩溃，连报错信息都没有。
原因：spectral 库依赖底层 C 扩展，在 Python 3.10 以上版本或者 Windows 上没有预编译包时，会触发内存错误。常见场景是 in-place pip install 装到了不匹配的 wheel。
解决：降级 Python 到 3.8 解决得最彻底。或者绕开 spectral，手工解析 .hdr 里记录的字节顺序和数据类型，配合 np.fromfile 读 .dat。这个场景下我会优先用后者，避免为了读一个文件重装环境。手工读取的代码我一般放在工具脚本里，读法不复杂但很容易在字节顺序上翻车——提前确认 hdr 里 byte order 字段是 0 还是 1。
6. 进阶验证：把预处理后的数据接上模型，量化每一步的提升幅度
6.1 用可视化对比预处理前后光谱曲线的实际差异
很多同学做完预处理直接拿模型跑分，从不看光谱曲线长什么样。这不对。我习惯把处理前和处理后的光谱曲线叠加画在一起，观察三件事：曲线是否平滑、吸收特征是否保留、不同类别间是否有可分性趋势。一张清晰的可视化图，比一百组精度数字更能说明预处理做没做对。
import matplotlib.pyplot as plt def plot_spectra_before_after(raw_spectra, processed_spectra, sample_index=0, title=''): plt.figure(figsize=(12, 6)) plt.plot(raw_spectra[sample_index], label='raw', alpha=0.6, linewidth=0.8) plt.plot(processed_spectra[sample_index], label='processed', alpha=0.9, linewidth=1.2) plt.xlabel('Band Index') plt.ylabel('Reflectance (normalized)') plt.title(title) plt.legend() plt.grid(True, linestyle='--', alpha=0.3) plt.show()逻辑说明：选一个典型的样本，把处理前后的曲线放一起，能直观看到平滑程度和基线偏移是否被修正。如果两条曲线几乎重合，说明预处理没生效，大概率是归一化方式没有把尺度差异消掉。参数说明：多抽几个不同类别的像元对比，别只看一个样本——不同地物的吸收特征是评价预处理效果的钥匙。
6.2 五折交叉验证：用实际分类精度收口验证
比单次实验更可靠的验证方式是五折交叉验证。它把数据分成五份，轮流取四份训练一份测试，五次结果取平均值。这样做能避开因为某次划分恰好好坏带来的性能虚高或虚低，是课程设计报告里让老师信服的扎实做法。
from sklearn.model_selection import cross_val_score from sklearn.svm import SVC from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler # 五折交叉验证评估预处理后的光谱分类效果 X_train, X_test, y_train, y_test = spatial_split(spectra_3d, labels_3d) clf = make_pipeline( StandardScaler(), SVC(kernel='rbf', C=100, gamma=0.01) ) scores = cross_val_score(clf, X_train, y_train, cv=5) print('五折交叉验证精度: %.4f ± %.4f' % (scores.mean(), scores.std()))逻辑说明：pipeline 把标准化和分类器串起来，避免在交叉验证过程中由于缩放器拟合整个训练集造成数据泄漏。SVC 参数里 C 控制正则强度，gamma 控制 RBF 核的宽度，这两个参数对高光谱分类影响很大。参数说明：C 和 gamma 的取值不是固定的，我先用默认参数跑一遍，再看准确率决定是否用网格搜索调参。调参时用 GridSearchCV 配合五折验证，是课程设计里最稳妥的收尾方案。
从那以后，我每次接手一份高光谱数据，都会强制自己先画原始光谱曲线，再决定处理顺序；每做完一步预处理，就把中间结果存档并输出 shape 校验。这成了我个人的一个习惯——眼见的曲线图永远比数字诚实的多。希望这次把流程拆开讲一遍，能让你在课设里少走几步弯路。
本文还有配套的精品资源，点击获取
