Python高光谱数据预处理实战：ENVI读取、反射率转换与光谱平滑

简介这是一份基于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 里只存紧凑的像素流。这种拆法的好处很直接高光谱波段数常到几百单文件轻松上 GBTIFF 会被 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, errorsreplace) 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, dtypedtype, counttotal, offsetoffset) # 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 处理时直接 axis2不用来回 transpose。提示整景数据非常大时把 np.fromfile 换成 np.memmap(dat_path, dtypedtype, moder) 再 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_ref0.98, solar_zenith_degNone): 把原始 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(axis2) if dark.ndim 3 else dark white_m white.mean(axis2) 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, requiredTrue, help目标影像 .hdr) ap.add_argument(--dat, requiredTrue, help目标影像 .dat) ap.add_argument(--dark, requiredTrue, help暗电流影像) ap.add_argument(--white, requiredTrue, 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_deg30.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_length7, polyorder3, axis2)两个参数决定一切。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_threshold30): 返回坏波段布尔掩码True 表示该波段不可靠。 wavelengths np.asarray(wavelengths) bad np.zeros_like(wavelengths, dtypebool) # 大气水汽吸收强波段物理上不可用 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, step40): plt.figure(figsize(10, 4)) plt.plot(wavelengths, spectrum, lw0.8) # 每隔 step 个波段显示一个标签否则全部标签会叠成黑块 idx list(range(0, len(wavelengths), step)) plt.xticks(wavelengths[idx], [f{wavelengths[i]:.0f} for i in idx], rotation45) plt.xlabel(wavelength (nm)) plt.ylabel(reflectance) plt.grid(alpha0.3) plt.tight_layout()step 取多少要看波段总数。400 个波段我取 40正好显示 10 个标签如果波段间隔宽也可以把 step 换成固定的波长间隔比如每 100nm 标一个点。加 rotation45 和 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 设为 7polyorder 设为 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(fNDVI 范围: {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也避免每次跑代码都用不同版本的预处理参数。最后再说一句经验预处理这东西参数没绝对标准关键是把每个中间结果存下来、画出来再进入下一步。把验证这一步养成肌肉记忆很多玄学精度的坑都能提前绕开。希望帮到你。本文还有配套的精品资源点击获取