简介：这是一份基于Python实现的高光谱数据预处理项目资料，面向毕业设计、课程设计及项目开发场景，适合需要完成光谱数据清洗、校正、降噪与归一化处理的开发者参考。资料以预处理算法集合为核心，覆盖标准正态变换(MSC)、多元散射校正(SNV)、SG平滑滤波、滑动平均、一阶差分、二阶差分、小波变换、均值中心化、标准化、最大最小归一化、矢量归一化等十一种常用方法，每种算法均附带可运行源码、说明文档和逐行代码解析，读者可结合示例清晰掌握各方法的数学原理、适用场景与调用方式。压缩包共含17个文件，其中2个Python脚本是核心实现，1个CSV文件提供光谱示例数据，12个PNG图片直观展示各算法处理前后的对比效果与整体流程，并配有License与Markdown说明文档，整体大小仅2.48MB，目录清晰、检索方便。目前该项目已有486人浏览学习，代码经过严格测试，可直接运行验证，既能帮助高光谱领域初学者快速建立预处理流程认知，也能作为相关课题或工程方案的基础框架进行二次扩展与集成。
1. 高光谱数据预处理用 Python 做，到底在解决什么问题
如果你拿一台推扫式高光谱相机去拍农田，采集卡写进磁盘的是 16 位整数 DN 值。这组数据看起来像“偏暗的灰度视频”，里面却混着暗电流噪声、传感器响应不均、镜头渐晕、天空光和大气水汽吸收。直接用原始 DN 建分类模型，精度往往不如一道简单归一化。用 Python 开发高光谱数据预处理方法并发成源码，价值就在于把从文件读取到反射率输出的流程固化成一套可复用代码，再配合文档和代码解析，让接手的人不用重新猜。这篇按我平时落地的顺序讲：读数据、定标转反射率、平滑去坏波段、验证。适合刚拿到无人机高光谱数据找不到入手点的人，也适合做课程项目时参考整体结构。
2. 高光谱数据在 Python 里怎么读：从 ENVI/BSQ 到内存数组
原始高光谱数据十有八九以 ENVI 格式存储。整套 Python 高光谱工程的第一步，是把二进制文件解析成形状确定的三维数组；这一步错了，后面所有反射率计算和建模都是在脏数据上做文章，再好的算法也救不回来。
2.1 为什么高光谱数据默认长在 ENVI 格式上
ENVI 格式不是一个单一图像文件，而是 .hdr 文本头加 .dat 裸二进制组合。.hdr 里记录 samples、lines、bands、data type、interleave、byte order 和 wavelength 这些元信息，.dat 里只存紧凑的像素流。这种拆法的好处很直接：高光谱波段数常到几百，单文件轻松上 GB，TIFF 会被 4GB 限制卡住，或者波段组织方式混乱，而 .dat 这种裸字节流在磁盘上更快、更可控。
无人机高光谱农业数据集现在很多，公开数据常见一景 400×400×250，原始 DN 文件约 400MB；如果是推扫式条带，行数可能上万，文件直接冲到几个 GB。遇到这种数据，第一件事永远是打开 .hdr 看 interleave，而不是急着写训练代码。interleave 决定文件里三个维度的排列顺序，读错会让你拿到一组“看起来像图，实际完全错位”的数组。
2.2 两条读取路线：现成 spectral 库和手写解析
常见做法是直接调 spectral 库，读 ENVI 格式只需两行：
from spectral.io import envi img = envi.open('scene.hdr', 'scene.dat') arr = img.asarray() # shape 是 (samples, lines, bands)spectral 库封装了绝大多数细节，适合快速看一眼公开数据。但实际工程里我更多是手写解析，原因有两个：一是现场环境 pip 装 spectral 时，scipy 二进制轮子经常出问题，明明 python 安装 numpy 库的方法很简单，一装 scipy 反而翻车；二是 spectral 返回的数组带着自己的 EnviImage 对象和坐标系统，要做分块读取、暗电流对齐、波长重采样时，多包一层对象反而碍事。手写解析虽然代码长，但每个字段都由你控制，出问题排查起来也直接。
2.3 一套解析代码，把 BSQ/BIL/BIP 统一成三维数组
下面这个模块我命名为 read_envi.py，可以直接抄进项目。它做三件事：读 .hdr 元信息、解析波长列表、把 .dat 读成 (samples, lines, bands) 的三维数组。
import numpy as np # ENVI header 里的 data type 编号，对应 numpy dtype ENVI_DTYPE_MAP = { 1: 'uint8', 2: 'int16', 3: 'int32', 4: 'float32', 5: 'float64', 12: 'uint16', # 大多数高光谱相机原始 DN 是 uint16 13: 'uint32', } def parse_envi_header(hdr_path): """解析 ENVI .hdr 文本，返回键值字典。""" info = {} with open(hdr_path, 'r', errors='replace') as f: for raw_line in f: line = raw_line.strip() if not line or '=' not in line: continue key, val = line.split('=', 1) info[key.strip()] = val.strip() return info def extract_wavelengths(info): """把 hdr 里的 wavelength 字符串变成 float 数组。""" if 'wavelength' not in info: return None try: return np.array([float(x) for x in info['wavelength'].split(',')]) except ValueError: return None def read_envi_binary(hdr_path, dat_path): """读 ENVI 二进制，统一返回 (samples, lines, bands) 和元信息。""" info = parse_envi_header(hdr_path) samples = int(info['samples']) lines = int(info['lines']) bands = int(info['bands']) dtype = np.dtype(ENVI_DTYPE_MAP[int(info['data type'])]) interleave = info.get('interleave', 'bsq').lower() offset = int(info.get('data offset', 0)) # 字节序：0 是小端，1 是大端，不处理会读出雪花噪声 if info.get('byte order', '0') == '1': dtype = dtype.newbyteorder('>') else: dtype = dtype.newbyteorder('<') total = samples * lines * bands raw = np.fromfile(dat_path, dtype=dtype, count=total, offset=offset) # reshape 顺序由 interleave 决定，这是全函数最核心的分支 if interleave == 'bsq': # 波段优先 cube = raw.reshape(bands, lines, samples).transpose(1, 2, 0) elif interleave == 'bil': # 行内波段连续 cube = raw.reshape(lines, bands, samples).transpose(0, 2, 1) elif interleave == 'bip': # 像元内波段连续 cube = raw.reshape(lines, samples, bands) else: raise ValueError(f'不认识的 interleave: {interleave}') wavelengths = extract_wavelengths(info) return cube, info, wavelengths这段代码有三个关键点。第一是 byte order，我在这上面花过一整个下午：某国产相机导出的 .hdr 里 byte order 写 1，也就是大端序，np.fromfile 没做字节序翻转时，图像像元呈雪花状跳动，当时还以为是镜头问题。第二是 data offset，有些数据文件头部嵌着一段产品信息，不跳过会让第一个波段整体错位。第三是统一返回 (samples, lines, bands)，让最后一个轴永远对应波段，后面用 scipy 的 savgol_filter 处理时直接 axis=2，不用来回 transpose。
提示：整景数据非常大时，把 np.fromfile 换成 np.memmap(dat_path, dtype=dtype, mode='r') 再 reshape，按需切片读取，不会一次性把几个 GB 塞进内存。对连续条带数据这招几乎必备。
2.4 建议的源码结构：方法、源码、文档怎么摆
既然标题写的是“方法+源码+文档+代码解析”，我一般把工程组织成下面这种结构，交接和自用都清楚:
hsi_preprocess/ ├── read_envi.py # ENVI 读取器 ├── calibrate.py # 辐射定标、暗电流校正、反射率转换 ├── smooth.py # Savitzky-Golay 平滑与坏波段掩码 ├── visualize.py # 画光谱曲线、假彩色图、NDVI 验证 ├── main.py # 命令行入口，串起整条流水线 ├── docs/ │ └── preprocessing.md # 每步使用的参数和现场采集要求 └── data/ # 原始数据不入库，只放样例和结果这个结构不复杂，但符合“源码+文档”的交付形态：read_envi.py 负责数据读取，calibrate.py 负责预处理核心，smooth.py 负责质量和坏波段，visualize.py 负责验证。docs/preprocessing.md 里写清每一步为什么这么做、参数怎么调，这样别人拿到源码能跑通，也能看懂。做课程项目时，建议在 docs 里补一张“处理前后光谱曲线对比”图，这是高分项目常见的评分亮点。
3. 高光谱如何转反射率：暗电流、白板和太阳高度角校正一条龙
读完数据只是把二进制变成了数组，真实的地物反射率还藏在 DN 值里。高光谱如何转反射率，核心是一条线性公式，但公式背后的现场数据决定成败。
3.1 一条公式把 DN 到反射率讲清楚
R = (DN_target - DN_dark) / (DN_white - DN_dark) × ρ_white
其中 DN_target 是目标像元原始值，DN_dark 是暗电流响应，DN_white 是标准白板在相同光照和积分时间下的响应，ρ_white 是白板的已知反射率，聚四氟乙烯白板一般在 0.98 左右。如果做地面光谱测量，还要在最后除以 cos(太阳天顶角)，把斜射光改正到垂直入射等效值。
这段公式看着简单，三个细节却常被写反。第一，先减暗电流再除白板差，顺序不能反，因为白板 DN 里同样带着暗电流；第二，暗电流要取多帧均值而不是用单帧，否则影阴影区域会大幅抖动；第三，最后乘以白板反射率而不是除以它，很多新手把 0.98 写在分母上，结果整幅反射率整体偏高。有了这条光路，后面代码正确与否就有参照系了。
3.2 现场采集习惯：暗电流、白板和太阳高度角
预处理能做好，一半功夫在采集现场。我一般到测区第一件事是盖镜头盖，用当前积分时间保存 100 帧暗电流；然后找一块平整的聚四氟乙烯白板放在目标附近，保证白板朝向天空方向，再保存 50 帧亮场。白板不能被影子挡到，也不能被太阳直射得过曝。很多公开的无人机高光谱农业数据集会自带白板帧，但格式五花八门，导入代码后仍然要手动选白板区域 ROI 去取均值。
太阳高度角这个参数，地面光谱仪用得多，无人机影像很少逐像元算。严格说机载推扫影像的每个像元对应不同的观测几何，需要做逐行逐帧的 BRDF 校正；项目里少有人做得那么细，常见做法是取整个架次的平均太阳高度角作为一个标量简单修正。如果目标是获取相对反射率而不是绝对反射率，这一项可以直接跳过，我在代码里把它做成可选项。
3.3 核心代码：dn_to_reflectance 的完整实现
这部分是整个预处理源码里最关键的一段，直接给了再往下走的信心：
import numpy as np def dn_to_reflectance(raw, dark, white, white_ref=0.98, solar_zenith_deg=None): """把原始 DN 立方体转成反射率。 参数 raw : (samples, lines, bands) 原始 DN 数组 dark : 暗电流，可以是 (samples, lines, bands) 或二维均值 white : 白板响应，与 dark 同尺寸 white_ref : 白板出厂反射率，常见 0.98 solar_zenith_deg : 太阳天顶角，None 表示不做余弦校正 """ raw = raw.astype(np.float32) # 转 float32，防止整型相减溢出 # 三维暗电流/白板先压成单通道均值，能明显压随机噪声 dark_m = dark.mean(axis=2) if dark.ndim == 3 else dark white_m = white.mean(axis=2) if white.ndim == 3 else white denom = white_m - dark_m # 白板在某些波段响应可能接近 0，先置 NaN 而不是直接除零 denom = np.where(denom <= 0, np.nan, denom) refl = (raw - dark_m) / denom * white_ref if solar_zenith_deg is not None: refl = refl / np.cos(np.deg2rad(solar_zenith_deg)) return np.clip(refl, 0.0, 1.0)这个函数的参数说明值得细看。raw 必须在函数内 astype(np.float32)，否则 uint16 相减时出现负值会向上溢出成很大的正数，反射率曲线会莫名其妙在暗区爆高。dark 和 white 支持三维原始图像，也支持预先算好的二维均值，这样在内存紧张的机器上可以先在外面压维再调用。white_ref 默认 0.98，但不同厂家白板可能是 0.99 或 0.95，拿到设备先看一眼出厂证书，别默认。
回到命令行入口，我在 main.py 里把它串成完整流程：
import argparse import numpy as np from read_envi import read_envi_binary from calibrate import dn_to_reflectance def main(): ap = argparse.ArgumentParser() ap.add_argument('--hdr', required=True, help='目标影像 .hdr') ap.add_argument('--dat', required=True, help='目标影像 .dat') ap.add_argument('--dark', required=True, help='暗电流影像') ap.add_argument('--white', required=True, help='白板影像') args = ap.parse_args() raw, _, _ = read_envi_binary(args.hdr, args.dat) dark, _, _ = read_envi_binary(args.dark, args.dat) white, _, _ = read_envi_binary(args.white, args.dat) refl = dn_to_reflectance(raw, dark, white, solar_zenith_deg=30.0) np.save('reflectance.npy', refl) print('reflectance shape:', refl.shape) if __name__ == '__main__': main()注意这里一次性载入 raw、dark、white 三个立方体，内存占用是单景的三倍。处理连续条带数据时我会把函数改成边读边算，或者在 read_envi_binary 里用 memmap 返回，然后逐行做定标。跑 OOM 并不总是因为数据大，很多时候就是进程里堆着三个没释放的数组。
3.4 全图定标和分帧定标的差别
如果一景影像是在很短时间内拍完的，光照变化可以忽略，全图用同一组暗电流和白板就行。但如果拍摄跨了半小时，太阳位置移动明显，最好把条带按时间切段，每段分别用邻近时刻的白板定标。我常用的做法是：飞行前采一次暗电流，飞行中每隔 10 分钟采一次白板，然后在代码里按文件时间戳或 GPS 时间对齐。否则正午前后的两段数据，反射率会出现肉眼可见的接缝，地面上没有变化，处理结果却像打了补丁。
4. 光谱平滑与坏波段剔除：让曲线跳变别再骗过模型
定标之后，反射率曲线仍然会有高频毛刺，还会在一些波段出现离谱的峰谷。如果不处理，分类模型可能把噪声当成特征，回归模型的曲线拟合也会被带偏。这一步属于“数据质量工程”，不算建模，但直接影响建模效果。
4.1 定标后要不要平滑，先看任务
平滑不是必须步骤。如果做的是矿物识别这类看重窄吸收特征的分类，过度平滑反而会抹掉关键诊断光谱；如果做植被参数反演或回归，适度平滑能明显提升稳定性。我判断标准是：看单像元光谱曲线的相邻波段差分，毛刺幅度超过信号幅度的 5% 就平滑，否则不动。千万别整景图无脑套滤波器，那是拿噪声换信息。
4.2 Savitzky-Golay 平滑的两个参数怎么定
光谱平滑领域我用得最多的是 Savitzky-Golay 滤波，它的特点是在局部做多项式拟合，比移动平均更能保留吸收峰形状。调用非常短：
from scipy.signal import savgol_filter # refl 是 (samples, lines, bands) 反射率数组 smooth_img = savgol_filter(refl, window_length=7, polyorder=3, axis=2)两个参数决定一切。window_length 必须是奇数，它代表参与拟合的相邻波段数；对 400-1000nm、波段间隔 2-5nm 的数据，一般取 7 到 17 之间。聚酯类矿物识别取小一点 5-7，植被光谱取 9 左右合适。polyorder 是多项式阶数，通常取 3；如果取到 5-7，拟合容易把噪声的局部起伏也当成趋势放大，曲线会出现新的假峰。经验法则：window_length 要比吸收特征宽度窄，这样才能保住 700nm 的红边和各类吸收谷。
4.3 坏波段掩码：水汽波段和低信噪比波段
高光谱数据里有几段是物理上注定不能用的：1350-1450nm 和 1800-2000nm 是大气水汽强吸收段，传感器在那里的信噪比几乎为 0。另外，每台相机的首尾波段因为探测器响应弱，噪声也大。我用下面这段代码自动生成坏波段掩码：
import numpy as np def find_bad_bands(wavelengths, reflectance_img, snr_threshold=30): """返回坏波段布尔掩码，True 表示该波段不可靠。""" wavelengths = np.asarray(wavelengths) bad = np.zeros_like(wavelengths, dtype=bool) # 大气水汽吸收强波段，物理上不可用 bad |= (wavelengths > 1350) & (wavelengths < 1450) bad |= (wavelengths > 1800) & (wavelengths < 2000) # 全图统计每个波段的信噪比：均值 / 标准差 band_mean = reflectance_img.mean(axis=(0, 1)) band_std = reflectance_img.std(axis=(0, 1)) snr = band_mean / (band_std + 1e-10) bad |= snr < snr_threshold return bad, snr注意这里生成的是掩码而不是直接删波段。删除波段会破坏数组连续性，后续波长重采样时要重新建索引；掩码则可以一直带在元数据里，建模时只在特征选择阶段把坏波段排除。信噪比阈值 30 是我常用的起点，如果整景数据整体偏暗，可以降到 20；如果是高配相机，调到 50 更合适。
4.4 画光谱曲线时横坐标太密集的解决办法
验证平滑效果时，很多人画光谱曲线发现横轴波长标签挤成一团，这是高光谱画图的经典问题，搜索“python 画图横坐标太密集”基本都是这个场景。解决思路是只挑一部分波长显示标签：
import matplotlib.pyplot as plt def plot_spectrum(wavelengths, spectrum, step=40): plt.figure(figsize=(10, 4)) plt.plot(wavelengths, spectrum, lw=0.8) # 每隔 step 个波段显示一个标签，否则全部标签会叠成黑块 idx = list(range(0, len(wavelengths), step)) plt.xticks(wavelengths[idx], [f'{wavelengths[i]:.0f}' for i in idx], rotation=45) plt.xlabel('wavelength (nm)') plt.ylabel('reflectance') plt.grid(alpha=0.3) plt.tight_layout()step 取多少要看波段总数。400 个波段我取 40，正好显示 10 个标签；如果波段间隔宽，也可以把 step 换成固定的波长间隔，比如每 100nm 标一个点。加 rotation=45 和 tight_layout 是为了避免标签互相压住，这是最常被忽略的一步。
5. Python 高光谱预处理避坑清单：5 个让数据作废的典型案例
这几条是我在真实数据上踩过的坑，每条都是“现象→原因→解决”的结构，希望能让你少走一段夜路。
5.1 内存崩了：大影像在 32 位 Python 里根本没得跑
现象：读入 .dat 后执行 reshape，直接 MemoryError，程序当场退出。
原因：32 位 Python 单进程内存上限约 2GB，而一景 400×400×250 的 uint16 数据就要 400MB，再加上 dark、white 两个数组，轻松超过 1.2GB，还没算反射率结果。
解决：换 64 位 Python，这是最根本的一步。如果机器确实只能跑 32 位，就把 read_envi_binary 里的 np.fromfile 换成 np.memmap，让数组按页从磁盘读，不要一次性进内存。处理连续条带数据时，哪怕 64 位环境我也用 memmap，防止并行任务把内存占满。
5.2 反射率越界：负值和大于 1 的现象从哪来
现象：白板区域反射率算出来 1.04，阴影区域是 -0.2，整体值域完全超出物理范围。
原因：暗电流没有做多帧均值，白板 DN 又没减暗电流就做分母；另一个常见原因是把白板 DN 直接当成 0-65535 去归一化，而不是除以真实辐射响应。
解决：回到 dn_to_reflectance 的顺序：先对暗电流求均值，再算 white-dark 差，最后乘以白板反射率。如果某些波段白板响应为 0，除以这些波段会造成异常，所以代码里先 np.where 把分母置为 NaN，再统一 clip 到 0-1。
5.3 数组形状对了但波段顺序错：模型会默默买单
现象：读数组后打印 shape 是 (400, 400, 250)，看着正常，但画单像元光谱曲线时，曲线像被撕碎了一样，毫无相邻波段连续性。
原因：ENVI 文件的 interleave 写的是 bip，代码里却按 bil 做了 reshape；数组维度一样，内部排列完全不同。人眼很难看出，分类模型照样能收敛，但特征含义已经乱了。
解决：不要只看 shape，要画一条已知地物的光谱曲线做检验。可以用植被像元验证 700nm 附近是否出现红边，或者直接用 find_bad_bands 的信噪比曲线看是否连续平滑。我把“画一段随机像元光谱”写进了 read_envi 的单元测试里，每次改解析代码都会触发。
5.4 暗电流采集时镜头盖没盖严
现象：反射率全场偏高，连水体的近红外波段都有 0.3，明显不合常理。
原因：暗电流帧采集时，镜头盖没盖严或镜头缝隙漏光，暗电流均值被抬到很高，随后 (raw - dark_m) 的结果被压缩。
解决：现场采集时用遮光罩再盖一层镜头盖，多保存 100 帧备选。代码里也可以加一道防线：如果 dark_m 的整体均值超过 raw 最小值的 20%，直接报错提示，不让流程静默往下走。
5.5 Savitzky-Golay 窗口调到 31，吸收谷全被擦掉
现象：平滑后的植被光谱曲线非常顺滑，但 700nm 红边坡度变缓，970nm 水汽吸收小谷消失。
原因：window_length 设为 31，相当于把 31 个波段做局部多项式拟合，而典型吸收特征的宽度才 10-20 个波段，特征被当成噪声抹平了。
解决：先把 window_length 设为 7，polyorder 设为 3，然后对比平滑前后的差分曲线。如果差分在吸收峰位置明显抬高，说明窗口太大；如果曲线依旧毛刺明显，再逐步加到 9 或 11。务必每次只调一个参数并保存中间结果，调参不后悔。
6. 用最近波段索引和 NDVI 对预处理结果做最后验证
预处理做完不是直接丢给模型，我习惯先做一个 30 秒的快速验证：找两三个典型波段算 NDVI，再画一条植被和一条土壤的光谱曲线。这一步能拦住多数定标错误，比如白板过期、暗电流漏光、波段顺序错位。
NDVI 验证用最近波段索引来选波段，比教条式地写死“第 67 波段”可靠得多，因为不同相机的波段中心波长不一样。代码是这样：
import numpy as np def nearest_idx(wavelengths, target): """在波长列表里找离目标最近的波段索引。""" return int(np.argmin(np.abs(np.asarray(wavelengths) - target))) nir_idx = nearest_idx(wavelengths, 800) red_idx = nearest_idx(wavelengths, 670) green_idx = nearest_idx(wavelengths, 550) nir = refl[:, :, nir_idx] red = refl[:, :, red_idx] ndvi = (nir - red) / (nir + red + 1e-6) print(f'NDVI 范围: {ndvi.min():.2f} ~ {ndvi.max():.2f}') print(f'中心波长为 {wavelengths[nir_idx]:.1f} / {wavelengths[red_idx]:.1f}')如果影像里包含农田或林地，NDVI 最大值应该在 0.5 到 0.9 之间；如果全场最高只有 0.2，反射率定标大概率偏了。我还会顺手用 visualize.py 里的 plot_spectrum 画出两类地物平均光谱，在 670-730nm 看有没有明显红边，在 970nm 附近看有没有水汽吸收小谷。这两个特征同时出现，说明暗电流、白板和太阳高度角校正都落在合理范围。
我的习惯是每个项目都输出一份 npz 格式的结果文件，里面存 wavelengths.npy、reflectance.npy、bad_bands_mask.npy 三样东西。这样建模阶段不需要重新读原始 ENVI，也避免每次跑代码都用不同版本的预处理参数。最后再说一句经验：预处理这东西，参数没绝对标准，关键是把每个中间结果存下来、画出来，再进入下一步。把验证这一步养成肌肉记忆，很多玄学精度的坑都能提前绕开。希望帮到你。
本文还有配套的精品资源，点击获取
实现精准ROI提取)
