简介针对高光谱数据预处理环节这套基于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), dtypebool) 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, window11, polyorder3): 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_lengthwindow, polyorderpolyorder) 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, transmittance0.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, axis0) refl atmospheric_correction_dark_object(rad, dark_mean)暗像元法在工程上的坑是np.percentile(flat, 2, axis0)选出的“暗像元”不一定是真实地面目标如果有云影或者传感器暗电流这个值会偏高或偏低。更稳的做法是人工标一块水体区域求均值。透过率 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(pathpreprocessing_config.yaml): with open(path, r, encodingutf-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(levellogging.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(fintermediate_{tag}.npy, cube) logger.info(fsaved 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_rows256): rows cube.shape[0] out np.empty_like(cube, dtypenp.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 降维把预处理后的立方体投影到前三个主成分上看看地物类别是否在特征空间里自然分离开如果分不开回来检查归一化方法是否选错了。很多人忽略了一个细节预处理对最终模型精度的提升往往比换一个更复杂的分类器更明显。我自己的习惯是永远保留预处理前和预处理后两份数据模型报告里同时给出两份数据的对比精度这样能直观看到预处理环节贡献了多少性能。这个习惯帮我在多次项目评审里免于被质疑。最后想提一个操作性很强的建议把预处理流水线的输入输出设计成标准化接口输入是原始影像输出是归一化后的光谱矩阵中间所有环节都封装成独立的可插拔函数。这样即使换一个传感器、换一套波段配置调整的只是配置文件和坏波段掩膜而不是重写整条流程。希望这篇笔记能帮你把高光谱预处理从玄学变成可复现的工程实践省下来的时间值得花在真正难啃的模型和地学解释上。本文还有配套的精品资源点击获取