高光谱数据预处理实战：从DN值到反射率的Python全流程

简介：这是一套面向高光谱数据分析与建模的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 直接 fromfile+reshape，结果特征波段错位，后续导数光谱长得像噪声。用 spectral 库读 ENVI 能自动按 hdr 里的 interleave 转轴，但如果想搞清楚源码里每一步在做什么，还是建议手动读一次。三种 interleave 的转换关系见下表。

在环境里装齐依赖，后面代码直接跑：

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(r'wavelength\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, dtype=np.float32) with open('scene.hdr', 'r', encoding='utf-8', errors='ignore') 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(encoding='utf-8', errors='ignore').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 映射：1=uint8 2=int16 12=uint16 4=float32 5=float64 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, dtype=dtype) 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(f'unknown 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_ref=0.99): # dark_frames / white_frames 形状都是 (n_frames, rows, cols, bands) sample = sample.astype(np.float32, copy=False) dark = np.mean(dark_frames, axis=0).astype(np.float32) white = np.mean(white_frames, axis=0).astype(np.float32) denom = white - dark # where 条件里的绝对值阈值避免除零；denom 为 0 的像元直接赋值 0 reflectance = np.divide( sample - dark, denom, out=np.zeros_like(sample, dtype=np.float32), where=np.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_mask=None, n_sigma=8, k=3): # cube: (rows, cols, bands)，按波段逐个处理 if dead_mask is None: # 每个波段的中心和离散度要分别算，不能拿整张 cube 的全局方差凑数 med = np.median(cube, axis=(0, 1), keepdims=True) std = np.std(cube, axis=(0, 1), keepdims=True) 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_sigma=8 是经验值，常规地物反射率动态范围不大，8 倍标准差能覆盖大多数孤立坏点，又不容易把真实亮目标（比如镜面反射）误杀。k=3 表示 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, win=9, poly=2): # 输入 (rows, cols, bands)，对最后一个轴做事 return savgol_filter(cube, window_length=win, polyorder=poly, axis=-1) smoothed = spectral_sg(reflectance, win=9, poly=2)

为什么 win=9 是常用起点？假设采样间隔 5 nm，9 个点覆盖 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, dtype=np.float32) return (x - x.mean(axis=1, keepdims=True)) / x.std(axis=1, keepdims=True) def msc(x, ref=None): x = np.asarray(x, dtype=np.float32) if ref is None: # 行业惯例：参考谱 = 全体样本的均值谱 ref = x.mean(axis=0) # 构造常数项与 ref 的设计矩阵，一次性对所有样本做最小二乘 design = np.column_stack([np.ones(ref.size), ref]) # (n_bands, 2) coef, *_ = np.linalg.lstsq(design, x.T, rcond=None) # 结果是 (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, deriv=1, win=7, poly=2, delta=5): # delta = 相邻波段中心波长间隔（nm） return savgol_filter( x, window_length=win, polyorder=poly, deriv=deriv, delta=delta, axis=-1 )

窗口选择比普通平滑更敏感，win=7、poly=2 在 5 nm 采样间隔下是比较稳的组合。如果波形仍有高频抖动，把 win 加到 9~11，但这时要回头检查吸收峰位置有没有位移。另一个注意点：导数光谱没有真实的 0-1 尺度，后续接分类模型几乎必须再做一次标准化，否则数值范围差异会盖过真正有用的光谱特征。

4.3 组合策略：不是所有方法都该同时上

有些源码把 SG、SNV、MSC、导数依次全部调用一遍，这种“全家桶”方案容易让模型过拟合到处理假象。可复现的组合方式通常是这样：SG 平滑永远放在最前面；SNV 与 MSC 二选一，不要同时用；导数作为最后一步特征增强；最终归一化建议用分位数裁剪后的 min-max，而不是裸 min-max，因为裸 min-max 在坏像元残余存在时会被单个异常值带偏。

端到端的流程通常不超过十行：读 cube、反射率校正、坏像元修复、SG、SNV/MSC、导数、分位数归一化，每一步都有对应数组状态。保持每步输出为 float32 的二维或三维数组，不要在中途变回整型，坑会少很多。

5. 高光谱预处理的五个踩坑现场：从 uint16 溢出到过平滑

5.1 坑一：无符号整型直接相减，负反射率卷成 65535

现象：黑帧校正后图像出现大量亮白色像素，反射率直方图尾部顶到 65535 附近。

原因：数据是 uint16，Python 里 uint16 减去 uint16 会发生无符号溢出，-1 会变成 65535。高光谱相机输出多是无符号整型，这是最隐蔽的数据坑之一。

解决：在任何减法之前先 astype(np.float32)。

import numpy as np dark = np.array([100], dtype=np.uint16) dn = np.array([50], dtype=np.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 nm，win=9 覆盖 40 nm 还算安全；win=21 覆盖 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, k=5): # 前 k 个主成分的解释方差比；预处理后一般应高于原始数据 return PCA(n_components=k).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 留存一份，方便随时回看。这个习惯救过我几次，特别是项目搁置几个月后再回来重跑时，参数表和结果文件能让所有步骤快速对齐。希望帮到你。

本文还有配套的精品资源，点击获取