简介：针对高光谱数据预处理环节，这套基于Python开发的完整项目源码与配套文档，适合进行毕业设计、课程设计或相关算法研究的学生和开发者使用。压缩包共17个文件，包含2个Python脚本、1个CSV样例光谱数据、12张说明图片以及License和Markdown文档，整体仅2.48MB，轻量易部署。资源内置标准正态变换(MSC)、多元散射校正(SNV)、Savitzky-Golay平滑、滑动平均滤波、一阶/二阶差分、小波变换、均值中心化、标准化、最大最小归一化、矢量归一化等十余种预处理算法，每种算法均配有可运行的代码实现与图示解析。代码已经过严格测试，可直接参考或扩展，方便读者对比不同预处理方法对光谱数据的影响，快速完成算法验证与论文实验。该资源已有485人学习下载，项目结构清晰，是入门高光谱数据分析与预处理的实用参考。
1. 高光谱数据预处理：为什么说它比建模更决定项目成败
拿到一个 .mat 或者 ENVI 格式的高光谱影像，里面有 100 多个连续波段，第一件事不是上模型，而是把数据“收拾干净”。高光谱数据预处理这件事，很多做遥感反演或者地物分类的 Python 开发者都栽过跟头：坏波段没剔、波长顺序对不上、反射率算出大于 1 的怪值，后边分类精度再高也是自欺欺人。我见过太多人在这上面反复返工，最后才意识到，预处理代码写得好不好，直接决定整个项目能不能落地。这篇笔记把一个基于 Python 的高光谱预处理方案拆开来讲，从数据读取、坏波段剔除、辐射定标到平滑归一化，每一步都给出可复现的代码和参数说明，适合正在做高光谱分类、矿物填图或者植被指数反演的从业者照着改。
2. 预处理链路拆解：从 DN 值到反射率，每一步在解决什么问题
2.1 高光谱数据的存储形态：ENVI、Mat 和 TIFF 三种载体的差异
高光谱数据最常见的三种载体是 ENVI 标准格式（.dat + .hdr）、MATLAB 的 .mat 文件，以及 GeoTIFF 多波段影像。ENVI 格式在遥感圈子里最通用，它的 .hdr 头文件里记录了波段数、行数、列数、数据位深和波长列表，但实际的数据排列还有 BSQ、BIL、BIP 三种区别，分别对应波段顺序优先、行顺序优先和像元顺序优先。用 Python 读取的时候，如果没按头文件里的 interleave 字段去解析，后面所有波段索引全部错位。
.mat 文件在学术数据集里很常见，比如 ICVL 高光谱数据集就是一堆 .mat 文件，每个文件里存了一个三维数组和对应的波长向量。用 h5py 读这类文件要注意，MATLAB 存储的数组默认是列优先（Fortran order），直接 NumPy 读出来维度是对的，但内存布局和 Python 的行优先不同，做逐波段操作时性能差异明显。TIFF 相对省心，但 16 位无符号整型存反射率时，经常有一个整体的缩放系数（比如 10000 倍），不做尺度还原，后面算植被指数全是错。
2.2 预处理四步链路：坏波段剔除、辐射定标、大气校正、平滑归一化
高光谱预处理的顺序不是随便定的，每步解决一个具体的信号问题。坏波段剔除放在最前面，因为水汽吸收波段（比如 1350~1450 nm、1800~1950 nm 附近）的信号基本是噪声，留着它们会让后面归一化的均值和方差被带偏。常见的做法是手动设定波长黑名单，或者用信噪比阈值自动识别。
辐射定标把传感器记录的 DN 值转成辐射亮度，这在数据头文件里通常有增益和偏置参数。然后是大气校正，把辐射亮度转成地表反射率，这一步是物理链路里最重的。Python 生态里没有像 ENVI 那样一键完成的 FLAASH 工具，常见做法是用 6S 模型的 Python 封装——比如 Py6S，或者自己实现简化暗像元法。如果研究区域和成像条件比较单一，暗像元法反而比完整大气模型更稳定，因为它依赖经验常数少，参数容易控制。
平滑和归一化放在链路最后。灰度共生矩阵或者光谱角匹配这类算法对噪声敏感，Savitzky-Golay 平滑可以在保留吸收峰的前提下滤掉高频噪声。归一化让每个像元的光谱向量变成单位长度或者零均值单位方差，保证后续分类器不因为波段的绝对辐亮度差异而偏向某个波段。链路顺序一旦颠倒，比如先归一化再剔除坏波段，坏波段的异常高方差会摊到每个波段上，后面排查的时候根本看不出是哪个环节出的问题。
3. 用 Python 从零跑通高光谱预处理：核心代码与参数调优
3.1 读取高光谱数据：spectral 和 h5py 的最小读取入口
Python 里读高光谱数据有两条主流路径。ENVI 格式用 spectral 库读最省事，它能直接解析 .hdr 头文件里的 interleave、波长信息和数据类型；.mat 格式则用 h5py 读，关键是把数组转成行优先并同步取出波长向量。下面是最小读取代码：
import numpy as np import spectral.io.envi as envi import h5py # 读取 ENVI 格式：envi.open 接收 .hdr 路径和数据文件路径 img = envi.open('scene.hdr', 'scene.dat') # 转成 (rows, cols, bands) 的 numpy 数组，数据类型由头文件自动决定 cube = np.array(img.load()) wavelengths = np.array([float(w) for w in img.metadata.get('wavelength', [])]) print('ENVI cube shape:', cube.shape, 'dtype:', cube.dtype) print('wavelength range:', wavelengths.min(), wavelengths.max()) # 读取 MAT 格式：ICVL 数据集常见结构是 rad 和 wavelength 两个 key with h5py.File('icvl_1.mat', 'r') as f: rad = np.array(f['rad']) # (cols, rows, bands) 的列优先存储 cube = np.transpose(rad, (1, 0, 2)) # 转成 (rows, cols, bands) wavelengths = np.array(f['wavelength']) print('MAT cube shape:', cube.shape, 'wavelengths:', wavelengths[:5], '...')这段代码的关键在np.transpose(rad, (1, 0, 2))。MATLAB 存数组是列优先，h5py 读出来后第一个维度实际对应的是列数，第二个维度才是行数，不做转置的话影像会旋转 90 度，逐波段操作时索引全乱。另外envi.open的 metadata 里 wavelength 是字符串列表，直接转 float 之前要检查有没有空字符串，否则整行报错。读大文件时不要直接.load()，先用img.read_band(b)逐波段读，把数据量控制在内存能承受的范围。
3.2 坏波段剔除与 Savitzky-Golay 平滑：参数怎么调
坏波段剔除的常见策略是维护一个波长区间黑名单，同时配合信噪比阈值自动检测。我一般先看头文件里的波长列表，再对照水汽吸收带的位置，把明显异常的波段索引记下来。水汽吸收带的判断可以参考公开的 HITRAN 数据库，但工程上更快的办法是直接看每个波段全影像的方差，方差异常低的波段基本就是坏波段。
from scipy.signal import savgol_filter # 按波长区间构造坏波段掩膜 def build_bad_band_mask(wavelengths, bad_ranges): mask = np.ones(len(wavelengths), dtype=bool) for lo, hi in bad_ranges: mask &= ~((wavelengths >= lo) & (wavelengths <= hi)) return mask bad_ranges = [(1350, 1450), (1800, 1950)] # 水汽吸收带，具体范围需按传感器调整 good_mask = build_bad_band_mask(wavelengths, bad_ranges) # 每个像元的光谱做 Savitzky-Golay 平滑 def smooth_cube(cube, window=11, polyorder=3): rows, cols, bands = cube.shape out = np.empty_like(cube) for r in range(rows): for c in range(cols): out[r, c] = savgol_filter(cube[r, c], window_length=window, polyorder=polyorder) return outSG 平滑的两个参数是最容易翻车的点。窗口长度必须小于光谱段有效长度，而且必须是奇数；多项式阶数一般取 2 或 3，阶数太高会把噪声当信号拟合。窗口取 11、阶数取 3 对多数 5~10 nm 分辨率的星载高光谱数据是安全的，但如果你的数据波段间隔很密，比如 1 nm 分辨率，窗口可以放大到 21。一个快捷的判断标准是：平滑后某条典型地物光谱的吸收峰深度如果明显变浅，说明窗口太大，已经磨掉了有效信息。
3.3 辐射定标与大气校正的工程替身：增益偏置与暗像元法
辐射定标本质是线性变换。大多数高光谱传感器的头文件里带了 gains 和 offsets，直接对坏波段剔除后的立方体做乘加即可。大气校正我倾向于用简化暗像元法：影像里找一块反射率近似零的暗像元，通常是清洁水体或者山体阴影，用它的每个波段均值作为大气路径辐射的近似，从所有像元里减掉，再除以每个波段的大气透过率估计。
# 辐射定标：DN -> 辐射亮度 def radiance_calibrate(cube, gains, offsets): cube = cube.astype(np.float64) return cube * gains[None, None, :] + offsets[None, None, :] # 简化暗像元法大气校正 def atmospheric_correction_dark_object(radiance, dark_mean, transmittance=0.9): # dark_mean 是暗像元区域每个波段的平均辐射亮度 corrected = (radiance - dark_mean) / transmittance return np.clip(corrected, 0, None) # 提取暗像元：选整幅影像第 2 百分位作为暗像元光谱 flat = rad.reshape(-1, rad.shape[2]) dark_mean = np.percentile(flat, 2, axis=0) refl = atmospheric_correction_dark_object(rad, dark_mean)暗像元法在工程上的坑是np.percentile(flat, 2, axis=0)选出的“暗像元”不一定是真实地面目标，如果有云影或者传感器暗电流，这个值会偏高或偏低。更稳的做法是人工标一块水体区域求均值。透过率 0.9 是经验值，对平原地区、晴空条件大致可用，山区或者有薄雾时，透过率会明显下降，建议按成像时间和大气能见度查 6S 模型的标准大气表来修正。
4. 把预处理代码工程化：模块划分、配置管理与文档落地
4.1 配置驱动：用 YAML 管理参数而不是改代码
预处理脚本写多了以后会发现，真正让人崩溃的不是算法，是参数散落在代码里。坏波段区间、SG 窗口、暗像元百分位、归一化方式，这些参数一旦写死在函数里，换一个数据集就要改代码、找函数、重新跑，出了结果对比不上。把参数抽到 YAML 配置文件里，用同一个脚本处理不同的数据源，是工程化的第一步。
# preprocessing_config.yaml data: format: envi # envi / mat / tiff scene_path: scene.dat header_path: scene.hdr good_bands: null # null 表示自动检测，也可以给 [0, 1, 2, ...] preprocess: radiation: apply: true gains: [0.001, 0.001, 0.001] offsets: [0.0, 0.0, 0.0] bad_band_mask: auto: true ranges: [[1350, 1450], [1800, 1950]] smooth: apply: true window: 11 polyorder: 3 normalization: method: minmax # minmax / zscore / none output: save_path: preprocessed.npy save_intermediate: truePython 侧用 PyYAML 加载配置，再把配置传给预处理函数：
import yaml def load_config(path='preprocessing_config.yaml'): with open(path, 'r', encoding='utf-8') as f: cfg = yaml.safe_load(f) return cfg cfg = load_config() cube = read_data(cfg['data']) processed = run_pipeline(cube, cfg['preprocess'])配置驱动的好处是换数据集不用动代码，只改 YAML 里的文件路径、坏波段区间和归一化方式。这个习惯越早养成越好，因为你做完一个项目三个月后回来看代码，唯一能帮你找回记忆的就是配置文件和日志，而不是那些写满魔法数字的函数。
4.2 日志记录与中间产物保存：给调试留后悔药
预处理链路长，中间状态多，调试的时候最怕的是不知道哪一步出了问题。我建议每个环节都把中间结果写盘，文件名带环节名和参数摘要。这样一旦最终结果不对，可以做二分排查——比如反射率出现负值，就去查辐射定标前的数据是否正常，而不是从头开始重跑。
import logging import numpy as np logging.basicConfig(level=logging.INFO, format='%(asctime)s %(levelname)s %(message)s') logger = logging.getLogger(__name__) def save_intermediate(cube, tag, cfg): if cfg['output']['save_intermediate']: np.save(f'intermediate_{tag}.npy', cube) logger.info(f'saved intermediate_{tag}.npy, shape={cube.shape}') # 在流水线各环节之间插入 save_intermediate cube_raw = read_data(cfg['data']) save_intermediate(cube_raw, '00_raw', cfg) cube_cal = radiance_calibrate(cube_raw, ...) save_intermediate(cube_cal, '01_radiance', cfg) cube_refl = atmospheric_correction_dark_object(cube_cal, ...) save_intermediate(cube_refl, '02_reflectance', cfg)中间产物都用 numpy 的.npy格式存，读写快而且不丢元信息。日志里除了记录文件路径，还应该记录每个环节的数值范围，比如处理前的 DN 值范围和处理后的反射率范围。如果哪天发现反射率最大值是 1.7，看日志就能快速定位是哪一步引入了异常放大的系数。
4.3 模块划分：把单脚本改造成可复用的小库
单个脚本写到底，最后一定会膨胀成一坨几百行的“面条代码”，后面的人根本不敢动。常见做法是拆成五个小模块：reader（负责不同格式的读取）、calibration（辐射定标）、correction（大气校正）、smoothing（平滑与归一化）、pipeline（组装整条链路）。模块之间只通过 NumPy 数组传递数据，不共享全局状态。
# pipeline.py from reader import read_envi, read_mat from calibration import radiance_calibrate from correction import dark_object_correction from smoothing import smooth_and_normalize def run_pipeline(cfg): if cfg['data']['format'] == 'envi': cube, wavelengths = read_envi(cfg['data']) else: cube, wavelengths = read_mat(cfg['data']) cube = radiance_calibrate(cube, cfg['preprocess']['radiation']) cube = dark_object_correction(cube, cfg['preprocess']['bad_band_mask']) cube = smooth_and_normalize(cube, cfg['preprocess']['smooth'], cfg['preprocess']['normalization']) return cube, wavelengths这样拆分的好处是单元测试好写。比如 calibration 模块可以用一组构造好的模拟 DN 值验证输出是否正确，correction 模块可以单独用只有暗像元的影像测试减法逻辑。代码解析这份文档的价值也在这里——每个模块的输入输出清晰，读者不需要在整个脚本里跳来跳去就能理解单个环节做了什么。
5. 高光谱预处理避坑指南：5 条典型的翻车记录与排查路径
5.1 反射率大于 1 或者出现负值：大气校正后的数值范围校验
现象：大气校正完的影像，反射率最大值跑到 1.5 以上，有些波段全是负值。 原因：暗像元法的暗像元光谱选得不纯，水体区域混了悬浮泥沙或耀斑，减掉的路径辐射偏小，校正后高亮像元就溢出 1；而某些波段暗像元本身有传感器暗电流，减多了就出负值。 解决：先做波段级别的暗像元置信区间筛选——对每个波段单独取第 2 百分位，而不是用整幅影像所有像元拉平；随后对校正结果做硬裁剪，反射率一律限制在 [0, 1] 之间。裁剪前务必看一眼直方图，如果大量像元挤在 0 或 1 的边界，说明暗像元选取或透过率参数还是不对。
5.2 波段顺序被悄悄打乱：ENVI 和 Mat 读出来的数据对不上
现象：从同一景数据分别用 ENVI 和 Mat 格式读取，两条光谱曲线在可见光段对不上，波峰位置错位几十纳米。 原因：ENVI 的 .dat 是按波段号从小到大存的，但 .hdr 里的波长列表有时是反序的；而 .mat 文件里数组的前两个维度在 MATLAB 中是列优先，读出来之后没有转置，影像变成了转置后的形状，逐波段操作时索引全错。 解决：不管哪种格式，读取之后立刻打印 shape 和波长数组的前三个值、后三个值，和已知参考数据做一次交叉核对。写一个assert wavelengths[0] < wavelengths[-1]，不满足就反转。这行断言救过我很多次，看起来简单但特别容易忘。
5.3 SG 平滑把吸收峰磨平了：窗口长度和阶数的搭配误区
现象：平滑后光谱变得很圆润，但叶绿素吸收峰深度从 0.35 降到 0.2，植物光谱的特征完全失真。 原因：窗口长度选得太大，比如把 200 个波段的整条光谱用一个 51 窗口去平滑，相当于做了一个大尺度的低通滤波；或者多项式阶数太高，把噪声的局部起伏也保留了下来，导致平滑效果名存实亡。 解决：先按传感器波段间隔估算窗口。经验值是窗口长度不超过有效波段数的 1/10，阶数取 2~3。每次平滑后选一条典型光谱（植被、水体、裸土各一条），把平滑前后的吸收峰深度打印出来对比，峰深变化超过 5% 就减小窗口。
5.4 整幅影像一次性加载把内存打爆：大影像的贪心加载陷阱
现象：一景 2000x2000x128 的 uint16 影像，load 进来直接内存占用超过 1 GB，做平滑时内存溢出不谈，程序直接卡死。 原因：全部波段的影像如果一次性进内存，再叠加 float64 转换和多个中间结果的拷贝，内存轻松翻 3~4 倍。Python 的内存管理在这种情况下不会自动释放前一个变量的数据。 解决：用分块处理思路，把影像沿行方向切成长度为 256 的块，逐块做预处理再拼回去。SG 平滑本身是逐像元的操作，不存在跨块依赖，所以分块几乎不影响结果。
def process_in_blocks(cube, block_rows=256): rows = cube.shape[0] out = np.empty_like(cube, dtype=np.float64) for start in range(0, rows, block_rows): end = min(start + block_rows, rows) block = cube[start:end].astype(np.float64) out[start:end] = smooth_and_normalize_block(block) return out5.5 坏波段掩膜与波段索引错位：掩膜长度对不上波段数
现象：事先准备了一个坏波段索引列表，应用到数据时却提示索引越界，或者剔除后波段数和波长列表不一致。 原因：波段索引是从 0 开始还是从 1 开始，不同数据集的元数据描述不一样。有些 .hdr 里的波段号是从 1 开始的，直接拿来当 numpy 索引，最后一个波段永远访问不到，反而把第一个波段重复访问了一次。 解决：统一在读取阶段把波段号强制转换成 numpy 从 0 开始的索引，并且生成的坏波段掩膜长度必须等于 cube.shape[2]。在应用掩膜之后加一条断言：assert cube.shape[2] == len(wavelengths[mask])。如果不等，优先回去检查数据读取阶段有没有做转置或者索引偏移。
6. 预处理效果验证：光谱曲线对比与信噪比评估
预处理做完不能直接扔给模型，得先验证结果在物理上合理。我最常用的一套验证流程分三步。第一步，随机抽 50 个像元，把处理前后的光谱曲线画在同一张图上，重点观察三条特征谱线——植被的绿峰、红边和叶绿素吸收峰是否完整保留，水体在近红外波段是否断崖式下跌。第二步，计算每个波段的信噪比，公式是波段均值除以波段标准差，如果某个波段信噪比低于 10，那它在后续分析里基本就是噪声源，宁可继续剔除。第三步，做一次 PCA 降维，把预处理后的立方体投影到前三个主成分上，看看地物类别是否在特征空间里自然分离开，如果分不开，回来检查归一化方法是否选错了。
很多人忽略了一个细节：预处理对最终模型精度的提升，往往比换一个更复杂的分类器更明显。我自己的习惯是永远保留预处理前和预处理后两份数据，模型报告里同时给出两份数据的对比精度，这样能直观看到预处理环节贡献了多少性能。这个习惯帮我在多次项目评审里免于被质疑。
最后想提一个操作性很强的建议：把预处理流水线的输入输出设计成标准化接口，输入是原始影像，输出是归一化后的光谱矩阵，中间所有环节都封装成独立的可插拔函数。这样即使换一个传感器、换一套波段配置，调整的只是配置文件和坏波段掩膜，而不是重写整条流程。希望这篇笔记能帮你把高光谱预处理从玄学变成可复现的工程实践，省下来的时间值得花在真正难啃的模型和地学解释上。
本文还有配套的精品资源，点击获取
