简介这是一套面向高光谱数据分析与建模的Python预处理方法集合尤其适合毕业设计、课程设计与相关课题研究。资源以pretreatment.py为核心集中实现了标准正态变换MSC、多元散射校正SNV、Savitzky-Golay平滑滤波SG、滑动平均滤波、一阶与二阶差分、小波变换、均值中心化、标准化、最大最小归一化、矢量归一化等十余种经典预处理算法并附带demo.py演示调用过程与peach_spectra_brix.csv真实光谱样本数据。整个压缩包共17个文件包含2个Python脚本、1个测试数据集、1个Markdown说明文档以及12张算法效果示意图整体仅2.48MB结构轻量却覆盖全面。目前已有487人学习使用代码经过严格测试配套文档与注释能帮助读者快速理解各方法的原理与适用场景并可在其基础上直接扩展或集成至自己的光谱分析流程。1. 高光谱数据预处理没那么玄先把 DN 数字值变成可分析的光谱把 .hdr 和 .dat 读出来满屏的数字直接拿去当反射率用训练出来的模型十有八九只对当天的光照有效。高光谱相机给的是 DN 数字值里面混着暗电流、白板起伏、坏像元、波段噪声和颗粒散射的多种影响。所谓基于 Python 的高光谱数据预处理就是把原始 cube 一步步清洗成能进分类器或回归模型的光谱矩阵这也是这个方向高分优秀项目里源码和代码解析通常围绕的核心。下面按我跑过多次的流程拆开新手能跟着步骤走熟手可以对照参数边界和几个非常容易翻车的细节。2. 数据组织与元数据解析动手前先弄懂 Cube 的排列与 .hdr 里的信息2.1 先弄清 interleave 的方向别等画图才发现波段错位高光谱数据在磁盘上不是我们想象的“立方体”而是按一定顺序拍平的一维文件ENVI 格式最典型其它 .raw 也多沿用这套约定。BSQ 是一个波段一整块BIL 是“一行内所有波段”BIP 是“一个像元内所有波段”。同样是物理上的 (rows, cols, bands) 立方体三种摆法读出来的临时数组维度分别是 (bands, rows, cols)、(rows, bands, cols)、(rows, cols, bands)。这个差别不是小细节。我接手过别人导出的 .raw没看 interleave 直接 fromfilereshape结果特征波段错位后续导数光谱长得像噪声。用 spectral 库读 ENVI 能自动按 hdr 里的 interleave 转轴但如果想搞清楚源码里每一步在做什么还是建议手动读一次。三种 interleave 的转换关系见下表。interleave文件内大致顺序直接 reshape 后的形状转成 (rows, cols, bands)BSQ波段 1 整块波段 2 整块(bands, rows, cols)transpose(1, 2, 0)BIL第 1 行的全部波段再第 2 行(rows, bands, cols)transpose(0, 2, 1)BIP每个像元的全部波段连续(rows, cols, bands)不需要转换在环境里装齐依赖后面代码直接跑pip install numpy scipy scikit-learn spectral matplotlib pandas安装完成后spectral 库可以帮你把 hdr 解析和读取都做掉但理解上面表格仍然值得因为很多“高分优秀项目”的源码其实是裸读写 .raw 的没有依赖这个库。2.2 .hdr 是自带的结构化文档解析 wavelength 比解析 shape 更重要ENVI 的 .hdr 看着像文本说明其实它同时决定 dtype、samples、lines、bands 和 interleave也带着 wavelength 与 fwhm。很多预处理流程只取了 shape把 wavelength 落掉等到需要画导数光谱、按波长裁剪波段时才发现还得回来补。我一般会把它当结构化文档解析优先取出波长列表后续 savgol_filter 的 delta 参数、按特征波段裁剪都靠它。高光谱 .hdr 里的 wavelength 经常写在一对大括号里而且可能跨行直接按行读会漏内容。下面是解析波长的一个小函数重点是把回车先替换掉再配合正则取大括号内容import re import numpy as np def parse_wavelength(hdr_text: str) - np.ndarray: # 关键点wavelength 可能被换行拆开先把回车换成空格再匹配大括号 text hdr_text.replace(\n, ) m re.search(rwavelength\s*\s*\{(.*?)\}, text) if m is None: return np.array([]) # 大括号内是逗号分隔的浮点也可能混着空项逐一过滤 values [float(v.strip()) for v in m.group(1).split(,) if v.strip()] return np.asarray(values, dtypenp.float32) with open(scene.hdr, r, encodingutf-8, errorsignore) as f: hdr_text f.read() wavelength parse_wavelength(hdr_text) print(len(wavelength), wavelength[:5], wavelength[-5:])这段代码里正则里的.*?做了非贪婪匹配避免把 hdr 后面其它大括号内容也吞进来。解析出来之后顺手检查首尾波长很多国产高光谱相机给出的单位是微米不是纳米数值在 0.4~2.5 区间的话记得乘以 1000否则后面按 nm 裁剪波段会全部偏移。2.3 裸读 .dat 的最小函数自动按 interleave 转成统一 shape这里给一个由 numpy 完成读取的函数也是后面全部预处理步骤的数据入口。手写一遍虽然 spectral.load 一行能解决但这是最容易排查的一个黑匣子数据读出来不对先怀疑 interleave而不是先怀疑滤波参数。from pathlib import Path import numpy as np def load_envi_cube(data_path, hdr_path): # 把 hdr 按行解析成简单 dict只取关键字段 meta {} for line in Path(hdr_path).read_text(encodingutf-8, errorsignore).splitlines(): if in line: key, val line.split(, 1) meta[key.strip()] val.strip().strip({}).strip() rows int(meta[lines]) cols int(meta[samples]) bands int(meta[bands]) # ENVI 常见 data type 映射1uint8 2int16 12uint16 4float32 5float64 dtype_map {1: u1, 2: i2, 12: u2, 4: f4, 5: f8} dtype np.dtype(dtype_map[int(meta[data type])]) interleave meta.get(interleave, BSQ).upper() raw np.fromfile(data_path, dtypedtype) if interleave BSQ: cube raw.reshape(bands, rows, cols).transpose(1, 2, 0) elif interleave BIL: cube raw.reshape(rows, bands, cols).transpose(0, 2, 1) elif interleave BIP: cube raw.reshape(rows, cols, bands) else: raise ValueError(funknown interleave: {interleave}) # 源数据哪怕是 uint16也立刻提升到 float32避免后续减法发生无符号溢出 return cube.astype(np.float32), meta cube, meta load_envi_cube(scene.dat, scene.hdr) print(cube.shape, cube.dtype)逻辑说明.dat 文件里没有维度信息维度全在 .hdr。BSQ 先按“波段→行→列”切再做一次 transpose 得到“行→列→波段”BIL 先按“行→波段→列”切再转成“行→列→波段”BIP 本身就是行、列、波段顺序不需要转。最后统一转 float32是为了让后面白板校正、坏像元检测里的减法都在有符号浮点上进行。实际使用中还有两个边界值得留意一是 .dat 前面可能带一段字节偏移常见是 1024 字节的文件头numpy.fromfile 要加上 offset 参数二是少数 hdr 会声明 byte order大端序数据要在 fromfile 前用 dtype.newbyteorder() 处理。这两种情况在公共数据集里不常见但自采数据时经常遇到。读进来之后我还会顺手做三个初检shape 是否符合 (rows, cols, bands)、全数据 min/max 是否在预期范围、NaN 占比是否异常。高光谱相机掉线或标定失败时NaN 或全零波段会直接污染后面的均值谱早发现比晚发现好处理得多。3. 预处理主流程黑帧、白板、坏像元、去噪照着这套代码跑3.1 黑帧白板校正把 DN 转成反射率的第一步DN 与地物反射率之间不是简单的线性系数关系。传感器有暗电流即使镜头全黑也有底噪光源和相机响应还会随时间漂移。只有先把暗帧减掉、再用白板归一才能产出可跨场景复用的反射率。常见做法是R (DN_sample - DN_dark) / (DN_white - DN_dark)严格一点还要乘白板标称反射率聚四氟乙烯白板一般标称 0.98~0.99。采集暗帧和白帧时通常连续拍多帧取平均单帧噪声会被摊薄。下面这个校正函数带防除零和防溢出适合直接抄进项目。import numpy as np def calibrate_to_reflectance(sample, dark_frames, white_frames, white_ref0.99): # dark_frames / white_frames 形状都是 (n_frames, rows, cols, bands) sample sample.astype(np.float32, copyFalse) dark np.mean(dark_frames, axis0).astype(np.float32) white np.mean(white_frames, axis0).astype(np.float32) denom white - dark # where 条件里的绝对值阈值避免除零denom 为 0 的像元直接赋值 0 reflectance np.divide( sample - dark, denom, outnp.zeros_like(sample, dtypenp.float32), wherenp.abs(denom) 1e-6 ) reflectance reflectance * white_ref # 白板欠曝或过曝时反射率会越界clip 到 [0,1] 当最后一道保险 return np.clip(reflectance, 0.0, 1.0).astype(np.float32)参数说明white_ref 常见取 0.99如果你的白板是别的反射率查厂商标称值传进来就行。暗帧和白帧我一般要求不低于 16 帧取平均之后随机噪声能压下去一截。np.divide里刻意写了 where 条件否则 w-d 为 0 的像元会出现 inf 或 nan后面的 SG 滤波会把 nan 扩散成整片黑洞。有些源码会省略这一步直接用 DN 做归一化。如果只是单张影像的快速分类模型也能学到东西但一旦要跨影像、跨时间复用模型不校正的模型几乎注定泛化失败。3.2 坏像元检测与修复宁可掩膜不要硬填高光谱探测器上总有少数像元响应异常表现是亮点、暗点或恒值点。它们不只是占几个像素而是会让每波段的均值、标准差和后续归一化全部偏掉。检测上可以对每个波段做“相对中位数的偏离度”判断修复时用空间邻域中值替代。为什么用中值不用均值因为均值遇到邻域里另一个坏点会被带偏中值天然抗拒离群值。from scipy.ndimage import median_filter def repair_bad_pixels(cube, dead_maskNone, n_sigma8, k3): # cube: (rows, cols, bands)按波段逐个处理 if dead_mask is None: # 每个波段的中心和离散度要分别算不能拿整张 cube 的全局方差凑数 med np.median(cube, axis(0, 1), keepdimsTrue) std np.std(cube, axis(0, 1), keepdimsTrue) dead_mask (cube med n_sigma * std) | (cube med - n_sigma * std) # size(k, k, 1)只在空间邻居里取中值不做光谱维平滑保护吸收谷 repaired np.where(dead_mask, median_filter(cube, size(k, k, 1)), cube) return repaired, dead_mask参数说明n_sigma8 是经验值常规地物反射率动态范围不大8 倍标准差能覆盖大多数孤立坏点又不容易把真实亮目标比如镜面反射误杀。k3 表示 3×3 邻域空间分辨率高的数据可以放到 5但邻域太大修复结果会发糊。mid_filter 的 size(k, k, 1) 是关键第三个维度是 1表示不做光谱维平滑否则会把相邻波段的吸收峰差异也抹掉。一个值得记住的边界坏点连片或呈条带状时不要用邻域中值硬修最好的做法是生成 mask 交给后续建模时直接忽略。硬填会制造出“看起来很整齐、实际上全是假数据”的光谱模型精度越高越要警惕。3.3 SG 滤波光谱去噪里最值得调的两个参数高光谱去噪首推 Savitzky-Golay 滤波它用一个小窗口内的多项式拟合中心点比移动平均更保峰形。需要调的参数只有两个但影响最大window_length窗口点数和 polyorder多项式阶数。窗口必须是奇数一般从 7 开始试阶数常见 2 或 3。窗口太长会把真正的吸收谷当噪声抹掉窗口太短噪声依旧。from scipy.signal import savgol_filter def spectral_sg(cube, win9, poly2): # 输入 (rows, cols, bands)对最后一个轴做事 return savgol_filter(cube, window_lengthwin, polyorderpoly, axis-1) smoothed spectral_sg(reflectance, win9, poly2)为什么 win9 是常用起点假设采样间隔 5 nm9 个点覆盖 40 nm和多数矿物 20~50 nm 的吸收带比较匹配。如果目标吸收峰很窄窗口要降到 5如果只做平滑不做导数11 也能接受但要先看一眼曲线。poly 选 3 更贴近缓变背景选 2 更保守两者差异在导数光谱里会被放大。验证方式不复杂把处理前后的同一条曲线叠画吸收谷位置发生偏移说明窗口选大了曲线还是毛刺说明窗口或阶数需要加。SG 滤波在整个预处理流程里的位置应当是黑帧白板校正之后、散射校正之前。因为校正后的反射率才具备物理意义这时候平滑才有意义反过来先平滑再校白板会把白板帧里的噪声也平滑进去。4. 散射校正与光谱增强SNV、MSC、导数光谱的用法和边界4.1 颗粒度与光程变化带来的“伪差异”SNV 和 MSC 在解决什么反射率数据里常见的一类非目标噪声是散射效应同一种物料颗粒大小、表面粗糙度、装填紧密程度不同会让光谱出现基线抬高或整体斜率变化。这不是化学吸收造成的却最影响建模泛化。SNV标准正态变量变换的思路是对每条光谱做一次“逐样本标准化”减掉该样本所有波段的均值再除以标准差把乘性强度变化和加性基线偏移都压到同一尺度。MSC多元散射校正的思路是用全体样本的平均光谱做参考对每个样本做线性回归 x a b * ref再用 (x - a) / b 作为校正结果等于把每个样本都搬到参考光谱的“姿势”上。def snv(x): # x: (n_pixels, n_bands)每个像素一条光谱不要使用全局统计量 x np.asarray(x, dtypenp.float32) return (x - x.mean(axis1, keepdimsTrue)) / x.std(axis1, keepdimsTrue) def msc(x, refNone): x np.asarray(x, dtypenp.float32) if ref is None: # 行业惯例参考谱 全体样本的均值谱 ref x.mean(axis0) # 构造常数项与 ref 的设计矩阵一次性对所有样本做最小二乘 design np.column_stack([np.ones(ref.size), ref]) # (n_bands, 2) coef, *_ np.linalg.lstsq(design, x.T, rcondNone) # 结果是 (2, n_samples) coef coef.T # (n_samples, 2) correct (x - coef[:, 0][:, None]) / coef[:, 1][:, None] return correct.astype(np.float32)SNV 不依赖参考谱MSC 依赖全体样本均值谱样本数量少时 MSC 会因参考谱不稳定而抖动。SNV 的副作用是会把真实强度信息去掉如果后面要做定量反演就尽量别用 SNV。我常做的顺序是快速验证用 SNV模型精度有提升但物理解释变差时换成 MSC 对比一次。上面 MSC 的批量写法相当于一次性对所有样本求解比 for 循环快很多design 的形状是 (n_bands, 2)x.T 是 (n_bands, n_samples)lstsq 解出来直接就是所有样本的截距和斜率。4.2 导数光谱基线漂移的杀手同时也是噪声放大器导数光谱是对光谱沿波长方向求微分。一阶导数消去常数基线偏移二阶导数能分辨重叠吸收峰。预处理中不建议用 np.diff 直接求导它会把噪声放大得没法看。常见做法是仍然用 SG 滤波求导把 deriv 参数设成 1 或 2并用 delta 指定波长间隔让导数带有“每纳米变化量”的物理单位。from scipy.signal import savgol_filter def spectral_derivative(x, deriv1, win7, poly2, delta5): # delta 相邻波段中心波长间隔nm return savgol_filter( x, window_lengthwin, polyorderpoly, derivderiv, deltadelta, axis-1 )窗口选择比普通平滑更敏感win7、poly2 在 5 nm 采样间隔下是比较稳的组合。如果波形仍有高频抖动把 win 加到 9~11但这时要回头检查吸收峰位置有没有位移。另一个注意点导数光谱没有真实的 0-1 尺度后续接分类模型几乎必须再做一次标准化否则数值范围差异会盖过真正有用的光谱特征。4.3 组合策略不是所有方法都该同时上有些源码把 SG、SNV、MSC、导数依次全部调用一遍这种“全家桶”方案容易让模型过拟合到处理假象。可复现的组合方式通常是这样SG 平滑永远放在最前面SNV 与 MSC 二选一不要同时用导数作为最后一步特征增强最终归一化建议用分位数裁剪后的 min-max而不是裸 min-max因为裸 min-max 在坏像元残余存在时会被单个异常值带偏。组合适用场景不建议使用的场景SG → SNV → 标准化土壤、粉末、颗粒物分类需要保留反射率强度的定量任务SG → MSC → 导数 → 标准化叶片探测、重叠吸收峰解析导数导致信噪比不足的小样本场景只做 SG → 分位数归一化数据本身质量好快速验证跨光照、跨仪器复用的模型端到端的流程通常不超过十行读 cube、反射率校正、坏像元修复、SG、SNV/MSC、导数、分位数归一化每一步都有对应数组状态。保持每步输出为 float32 的二维或三维数组不要在中途变回整型坑会少很多。5. 高光谱预处理的五个踩坑现场从 uint16 溢出到过平滑5.1 坑一无符号整型直接相减负反射率卷成 65535现象黑帧校正后图像出现大量亮白色像素反射率直方图尾部顶到 65535 附近。原因数据是 uint16Python 里 uint16 减去 uint16 会发生无符号溢出-1 会变成 65535。高光谱相机输出多是无符号整型这是最隐蔽的数据坑之一。解决在任何减法之前先 astype(np.float32)。import numpy as np dark np.array([100], dtypenp.uint16) dn np.array([50], dtypenp.uint16) print(dn - dark) # 失控无符号溢出 # 65546 print(dn.astype(np.float32) - dark.astype(np.float32)) # -50.0重点在于统一入口自定义 load 函数里读进来就转 float32后面所有中间数据都不要掉回整型。黑帧校正、坏像元检测、均值计算这些步骤一旦出现溢出后续所有统计量都会错而且很难肉眼发现。5.2 坑二interleave 判断失误光谱在空间维和波段维之间串位现象读出来的 cube shape 是对的但首波段显示出来像横纹多数波段边缘出现错位。原因手动裸读 .dat 时没查 hdr 的 interleave默认按 BSQ 处理实际采集软件存的是 BIL 或 BIP。解决用前面 2.3 的函数读或者至少打印 hdr 里的 interleave。拿到数据后找一个强吸收特征验证排列是否正确比如 1450 nm 附近的水汽吸收带看该波段二值化后是否在预期区域成片。如果特征形状和空间分布对不上多半是 interleave 错了。5.3 坑三SG 窗口比吸收峰还宽吸收谷被抹成平肩现象平滑后曲线很干净但分类精度反而下降特征波段的贡献度消失。原因SG 窗口过大把有意义的吸收谷当噪声拟合掉了。高光谱吸收峰半高宽通常 10~50 nm窗口覆盖宽度超过半峰宽两倍时就危险。解决用波长间隔乘以窗口点数估算覆盖宽度从覆盖半峰宽 1~1.5 倍的窗口开始调。比如采样间隔 5 nm、吸收峰半宽 20 nmwin9 覆盖 40 nm 还算安全win21 覆盖 100 nm基本不可接受。正确做法是画几张不同 win 的导数谱对比吸收峰位置稳定时再定参。5.4 坑四坏像元修得太干净等于给模型伪造特征现象修复后的训练集精度大幅提升换一个批次数据直接崩。原因邻域中值修复在坏像元占比高时会把异常值替换成周边均值相当于人为制造低方差区域模型学到的是“这里有修复痕迹”而不是地物差异。解决修复前统计 dead_mask 比例超过 5% 的波段落不要修直接用 mask 排除。修复后保留 mask 文件模型验证时把测试数据对应区域一并排除避免信息泄漏。宁可让模型少几个像素也不要让它学修复痕迹。5.5 坑五白板和暗帧采集条件不一致反射率出现负数或大面积超 1现象校准后反射率大量为负值或者整体大于 1。原因白帧和样本不是同一照明条件下采集的或者白板过曝导致 DN_white 饱和。白板饱和时分母偏小反射率会被放大到大于 1灯光变化使暗帧不匹配也会出现负值。解决采集时养成“每拍一组样本立刻拍白板和暗帧”的习惯。代码里 clip 到 [0,1] 只是后悔药真正的修复是从采集流程上保证条件一致。如果白板最大 DN 已经接近饱和值比如 uint16 超过 60000这组白帧要作废重采而不是继续用。6. 预处理完了别急着训练三项验证与参数复现习惯6.1 先验证曲线、统计量和特征空间第一项是曲线叠加把同一条光谱处理前和处理后叠在一起看吸收谷位置不能位移平滑后的曲线不应该出现新拐点如果吸收峰变了回第 3 章调 SG 窗口。第二项是统计量检查对每个波段计算均值、标准差和信噪比预期处理后 SNR 有抬升、坏波段比例下降如果标准差反而增大多半是白板校正时引入了放大噪声。第三项是特征空间检查对某个地物类别的像元做主成分分析看同类是否更聚拢、异类是否拉开。import numpy as np from sklearn.decomposition import PCA def variance_ratio(x, k5): # 前 k 个主成分的解释方差比预处理后一般应高于原始数据 return PCA(n_componentsk).fit(x).explained_variance_ratio_.sum() raw_flat cube.reshape(-1, cube.shape[2]) clean_flat pipeline_output.reshape(-1, cube.shape[2]) print(raw explain ratio:, variance_ratio(raw_flat)) print(clean explain ratio:, variance_ratio(clean_flat))这个比值不是越高越好但如果处理后的解释方差比显著下降说明预处理过程把信号也当噪声去掉了需要回头检查是哪一步过平滑。6.2 参数配置写进 JSON别靠记忆重跑我一般会在项目里放一个 preprocess_config.json记录暗帧数、white_ref、SG 窗口、poly、是否用 SNV/MSC、导数阶数、波段裁剪范围。换数据或论文返修要重新生成图时直接加载这套配置不靠回忆。每跑一次把预处理前后平均光谱 CSV 留存一份方便随时回看。这个习惯救过我几次特别是项目搁置几个月后再回来重跑时参数表和结果文件能让所有步骤快速对齐。希望帮到你。本文还有配套的精品资源点击获取