高光谱数据预处理Python代码包：从辐射定标到光谱平滑的完整实现

简介本资源是一套面向人工智能与机器学习初学者及遥感图像处理实践者的高光谱数据预处理Python代码集聚焦解决高维光谱数据噪声大、维度高、特征隐含深等典型建模前瓶颈问题。压缩包共15个文件2.48MB含2个核心Python脚本pretreatment.py实现光谱校正、去噪、平滑、标准化与PCA降维等9类主流预处理逻辑demo.py提供端到端调用示例、1个真实桃子光谱CSV数据集peach_spectra_brix.csv及12张可视化结果图涵盖原始/去噪/平滑/降维前后对比、光谱曲线与特征分布热力图等直观呈现各步骤效果。已有1370人学习下载代码结构清晰、注释完整可直接运行复现全部流程配套图像覆盖从数据加载、异常检测到特征选择的全链路操作为后续分类、回归或目标识别任务提供高质量输入基础。 做高光谱的朋友一定清楚数据预处理这步在整个流程里有多磨人。我手里的这份“高光谱数据预处理方法python代码.zip”就是把那些绕来绕去的标准和经验整理成了一套可以直接跑的Python脚本。这不是什么花哨的框架就是一个按照遥感数据常规流程组织的代码包从最原始的DN值读取到辐射定标、反射率转换、坏波段剔除、去噪平滑再到数据导出每一步都单独做了封装。对于刚入门的遥感、地质或农业方向的同学以及被一堆.hdr和.dat文件搞到头大的从业者来说直接对着这份代码改改路径和参数比从零开始查文档要省太多时间。这篇内容我会把这套代码包的设计思路、核心函数逐个拆开讲包括各个模块的“为什么这么做”以及实际运行中一定会遇到的坑。强烈建议你自己把代码跑一遍再对照文章去理解每个参数的意义这样后面换到自己的数据上才不会只是“能跑”而不知道“怎么调”。1. 高光谱数据预处理的“为什么”与整体流程拆解1.1 原始影像为什么不能直接拿来用绝大多数卫星或无人机高光谱传感器直接下传的是量化后的DN值即传感器记录的数字量化值。DN值本身没有物理意义它受传感器增益、暗电流、光照条件、大气散射吸收等一系列因素影响。你拿DN值去做分类、反演、矿物填图出来的结果不仅不同时间、不同传感器之间没法对比就算是同一景影像内部的暗目标和亮目标也没法用同一套阈值去处理。传感器定标文件里通常会给一组辐射定标系数把DN值转换成传感器入瞳处的辐射亮度单位一般是W·m⁻²·sr⁻¹·μm⁻¹。这只是第一步。辐射亮度还包含了大气路径辐射、大气吸收散射的影响真正到地表像元反射率还得做大气校正或至少一个简化的经验线性校正。这套代码里把这两个步骤都做了“可运行版本”的封装辐射定标放在第一步反射率转换放在第二步这样你至少能在不依赖ENVI或其它专业软件的情况下得到一个物理意义相对合理的结果。另外还有一点经常被忽略就是传感器自身的坏波段。像Hyperion这类早期星载高光谱传感器数据里有很多波段是坏的要么全零要么噪声极大尤其在1350-1450nm和1800-1950nm附近的水汽吸收带。这些波段如果不提前剔除后续无论是主成分分析还是一阶微分变换都会被它们带偏所以代码包里有一个专门的坏波段列表文件你自己拿到数据后一定要根据传感器的波段定义去修改它。1.2 一个标准的预处理流程应该是什么样从拿到原始高光谱影像到把数据喂给分类模型或光谱分析工具我建议按下面这条主线来走原始数据读取同时读取头文件里的波段信息、波长单位、数据存储格式如果是多景拼接或镶嵌前先做波段子集裁剪和坏波段剔除辐射定标把DN值转成辐射亮度大气校正或经验线性法把辐射亮度转成地表反射率光谱平滑去噪例如Savitzky-Golay滤波或移动均值平滑可选的光谱增强变换例如标准正态变量变换SNV、一阶导数、连续统去除按研究区裁剪、掩膜、重采样最后输出为ENVI标准格式或GeoTIFF这套代码包基本就是照着这条主线写的但为了保持脚本独立性模块之间都用中间文件或numpy数组衔接你可以只跑其中的某一步不必每次都从头来。这对调试非常友好。1.3 这份代码包的整体模块划分打开压缩包后你会看到脚本是按数字前缀排好的这样运行顺序很清晰00_check_data.py检查头文件和原始数据是否完整01_read_envi.py读取ENVI标准格式同时可视化几个关键波段的灰度图02_radiometric_calibration.py辐射定标根据增益和偏置计算辐射亮度03_atmospheric_correction.py简化大气校正输出地表反射率04_band_removal.py剔除坏波段和水汽吸收带05_savitzky_golay.py光谱平滑06_snv_derivative.py标准化和导数变换07_export_tiff.py结果输出这种模块化拆法是从实际项目里沉淀下来的。早期我自己写处理流程时图省事把代码全塞在一个大脚本里结果每换一个传感器就得重写一遍波段参数、定标文件、波长范围全揉在一起改一个地方就担心碰坏另一处。后来干脆按流程节点拆开每个脚本只干一件事参数统一放一个config.py里这才算真正解脱。这套代码包也沿用了这个习惯建议你拿到后不要破坏这个结构。2. 环境准备与核心工具库选型解析2.1 别急着pip install先把依赖理清楚这套代码包依赖的主要库是numpy、spectral、rasterio、scipy、matplotlib再加一个可选但强烈推荐的scikit-learn。下面我把它们分别用来干什么、是否必需写得再明白一点库名用途是否必需numpy底层数组运算所有处理的基础必需spectral读取ENVI格式的.hdr和.dat文件这是遥感数据最常用格式之一必需rasterio输出GeoTIFF以及读取带地理坐标的TIFF数据必需scipySavitzky-Golay滤波以及一些信号处理函数必需matplotlib可视化单波段图像和光谱曲线调试时非常好用建议安装scikit-learn如果后续做PCA、MNF降维或简单分类可复用可选安装时特别注意spectral这个库在Python 3.10以上版本有时会有些小兼容问题如果你pip install spectral之后import报错多半是版本冲突。我遇到比较稳的组合是Python 3.9或3.10配numpy 1.23-1.26。如果你用的是较新的Python 3.11或3.12建议优先用conda创建一个独立环境把python版本钉在3.10能省掉很多莫名其妙的报错。2.2 关于PRISMA、Hyperion等不同传感器数据的读取差异代码包默认按ENVI标准格式读取因为它最通用。但高光谱数据实际来源很杂常见几种HyperionL1T产品是HDF或GeoTIFF需要用USGS提供的工具转成ENVI格式再读取PRISMA官方下载的产品通常是HDF5格式里面有多个数据集需要先提取出VNIR和SWIR的辐射亮度数据再做波段合并EnMAP也常见HDF5格式但波段定义和填充值跟PRISMA不一样机载或地面高光谱很多厂商软件输出直接就是ENVI格式或者带一个自定义头文件如果你拿到的是PRISMA或EnMAP的HDF5产品代码包里的01_read_envi.py就用不了需要额外写一个HDF5读取脚本。我个人的建议是不管原始格式是什么第一步先统一转成ENVI标准格式也就是生成一个.hdr文本头文件和一个.dat二进制文件之后所有处理都走同一套流程。这样做的好处是后续脚本完全不用改坏波段列表、定标参数也更好统一管理。HDF5转ENVI格式时最容易出问题的是波长单位。高光谱头文件里wavelength单位通常要求是纳米而有些新传感器产品给的是微米或厘米⁻¹转换时一定要先看清楚元数据说明。我这里踩过坑当年处理一批数据时没注意单位反射率曲线整个平移了几百纳米起初还以为是大气校正参数设错了查了好久才定位到是波长单位的问题。2.3 目录组织与配置文件说明拿到压缩包解压后建议保持以下目录结构preprocess/ ├── config.py ├── 00_check_data.py ├── 01_read_envi.py ├── 02_radiometric_calibration.py ├── 03_atmospheric_correction.py ├── 04_band_removal.py ├── 05_savitzky_golay.py ├── 06_snv_derivative.py ├── 07_export_tiff.py ├── bad_bands.txt ├── calibration.txt └── data/ ├── input.hdr ├── input.dat └── output/config.py就是所有参数的集中地包括输入文件路径、输出路径、传感器名称、波长范围、定标文件路径、是否做坏波段剔除等。这样做的好处是其他脚本基本不需要改通用逻辑写在代码里个性化差异全部收敛到配置文件里。这其实是我后来不管做什么项目都坚持的习惯参数分散在多个脚本里是对自己耐心的最大消耗。bad_bands.txt的格式很简单每行一个波段序号或一个区间例如“1-4”表示第1到第4波段或直接写“150”表示第150波段。calibration.txt里则保存每个波段的增益和偏置按“波段号 增益 偏置”的格式一行行写用于辐射定标。如果你手上暂时没有精确定标系数也可以先用一个粗略的全局增益偏置跑通流程后续再替换。3. 核心模块代码解析与实操要点3.1 读取ENVI标准格式同时别忽略头文件里的隐藏信息ENVI格式的特点是二进制数据文件和一个文本头文件两者同名扩展名不同。头文件里记录了行数、列数、波段数、数据类型、文件存储格式BSQ/BIL/BIP、波长信息、以及地图投影信息。使用spectral库读取非常简单import spectral.io.envi as envi img envi.open(data/input.hdr, data/input.dat) data img.asarray() print(f数据维度: {data.shape}) print(f波段数: {img.nbands}) print(f波长列表: {img.bands.centers})这段代码看起来简单但有几个隐藏点需要说明。asarray()返回的是一个三维数组维度顺序是(行, 列, 波段)。但如果你处理的是BSQ格式数据实际存储顺序可能是(波段, 行, 列)spectral库在读取时会自动处理这个差异返回给用户的都是(行, 列, 波段)。如果你图省事直接自己用numpy的fromfile读二进制那就必须搞清楚存储顺序否则数据维度对不上图像会完全错乱。再一个容易忽略的是img.bands.centers返回的是波长中心位置列表。有些头文件里波段信息可能缺失或者不是按波长排列spectral会直接抛异常。所以00_check_data.py里我写了个预检查先打印前几个和后几个波长值帮你确认波段顺序是否正常。这个环节看似多余实际能救你命尤其当你从网上下载的数据头文件本身被改过时。3.2 辐射定标与反射率转换的实现细节辐射定标是让高光谱数据从“记录值”变成“物理量”的关键一步。核心公式很简单[ L gain \times DN offset ]在代码里对应的是def radiometric_calibration(image, gains, offsets): # image: (rows, cols, bands) int16数组 # gains和offsets: 每个波段的浮点数组长度等于波段数 calibrated image.astype(np.float32) * gains[np.newaxis, np.newaxis, :] offsets[np.newaxis, np.newaxis, :] return calibrated这里有个维度技巧gains是长度为波段数的一维数组通过np.newaxis把它扩展成(1, 1, bands)实现对整个三维数组的广播计算。这样做比写for循环遍历每个波段快非常多尤其是数据量大的时候。大气校正比较复杂要做到精度很高一般得用FLAASH、6S、MODTRAN等专业模型需要输入大气水汽含量、气溶胶光学厚度、成像日期、经纬度、传感器高度等参数。但在这套代码包里我提供一个简化的经验线性法实现适合在没有完整大气参数时的快速处理。核心思路是用两个或多个已知反射率的地面目标建立辐射亮度和反射率之间的线性关系[ \rho a \times L b ]a和b通过地面实测反射率与影像对应像元的辐射亮度来拟合。如果你手头没有实测数据也可以用暗像元和亮像元的经验值近似精度差一点但至少能让光谱曲线形态变得合理。代码实现如下def empirical_line_correction(radiance, dark_ref, bright_ref, dark_dn, bright_dn): # dark_ref/bright_ref: 地面实测反射率 # dark_dn/bright_dn: 影像中暗亮目标对应的辐射亮度 slope (bright_ref - dark_ref) / (bright_dn - dark_dn) intercept dark_ref - slope * dark_dn reflectance slope * radiance intercept return np.clip(reflectance, 0, 1)实操中特别注意clip到0和1之间看似粗暴但确实必要。因为简化校正后有些像元反射率会出现负值或大于1的情况如果不处理后面做对数变换或比值的算法会直接崩溃。3.3 坏波段剔除与光谱平滑哪些波段该留哪些该扔坏波段剔除不能一视同仁它取决于你的应用场景。如果是做矿物识别VNIR和SWIR区域的硅酸盐、碳酸盐、铁氧化物特征吸收带必须尽量保留即使有些噪声比较大的波段也只能平滑不能整段删除。如果是做植被监测红光和近红外区域尤其重要而水汽吸收带则基本是废的。代码里的bad_bands.txt建议你按下面这张表的规律来组织波段范围处理建议400-450nm信噪比一般偏低视情况保留用于大气校正验证1350-1450nm水汽强吸收通常直接剔除1800-1950nm水汽和二氧化碳吸收噪声大通常剔除2400-2500nm部分传感器信噪比急剧下降视应用决定传感器特定坏波段全零或固定值必须剔除剔除操作在代码里其实是靠一个布尔索引数组来完成的def remove_bad_bands(data, bad_bands): valid_mask np.ones(data.shape[2], dtypebool) for band in bad_bands: if isinstance(band, range): valid_mask[band] False else: valid_mask[band] False return data[:, :, valid_mask], np.where(valid_mask)[0]平滑方面我默认给的是Savitzky-Golay滤波。这个滤波器比简单移动平均好的地方在于它通过局部多项式拟合来保留光谱的峰谷细节不会把细微的吸收特征抹得太平。scipy里直接有现成函数from scipy.signal import savgol_filter smoothed savgol_filter(spectrum, window_length9, polyorder2)window_length和polyorder怎么选经验是window_length取奇数一般7-15之间polyorder取2-3不要超过4。窗口太窄去噪效果差窗口太宽会把有用的窄吸收带也磨掉。我实测下来对大多数地物光谱9点窗口配合2阶多项式是个比较稳妥的起点。如果你处理的是矿物光谱吸收带比较尖锐建议把窗口缩小到5-7如果是植被光谱曲线相对平缓窗口可以放宽到11-15。3.4 标准正态变量变换、一阶导数与数据归一化代码包里的第6个脚本主要做两件可选操作SNV和导数变换。SNV最初是近红外光谱分析里常用的一种预处理方法目的是消除固体颗粒大小、表面散射以及光程变化对光谱的影响。它的算法很简单就是每个光谱减去该光谱均值再除以该光谱的标准差def snv_transform(spectrum): mean_val np.mean(spectrum) std_val np.std(spectrum) if std_val 0: return spectrum return (spectrum - mean_val) / std_val对高光谱影像来说SNV对每个像元逐波段执行处理完后每条光谱曲线都被标准化到均值0、方差1附近不同亮度的像元之间变得更可比。但它也有缺点就是把原始反射率的绝对幅值信息抹掉了如果你后期要做生物量或理化参数定量反演SNV要慎用因为幅值本身可能是关键信息。一阶导数变换能有效去除基线漂移和线性背景干扰同时突出光谱曲线的斜率变化特征。导数计算其实就是一个有限差分但直接相邻波段差分会放大噪声所以通常搭配平滑一起用或者在差分前先对光谱做插值derivative np.gradient(smoothed_spectrum, wavelength_spacing)这里wavelength_spacing是波段间隔单位是纳米你在config.py里要提前配好。如果数据是重采样后的等间隔波段这个参数就是常数如果波段间隔不均匀那最好改成逐点计算的差分否则导数曲线会有高频噪声。归一化在这里我默认按波段最大值归一化也就是每个波段全图所有像元的反射率除以该波段的全局最大值。这样做的好处是不同波段之间动态范围统一后续做聚类或分类时各波段的权重更均衡。当然如果你要做的是光谱角匹配那么逐像元向量归一化更合适这两个别搞混了。3.5 数据导出保存成GeoTIFF保留坐标系和分辨率处理完的数据如果只存在numpy数组里换台电脑就全没了。所以07_export_tiff.py里我用rasterio把结果写回带地理参考的GeoTIFF文件。这里需要从原始ENVI头文件里读取地理变换参数和投影信息import rasterio from rasterio.transform import from_origin transform from_origin(origin_x, origin_y, pixel_width, pixel_height) with rasterio.open( data/output/reflectance.tif, w, driverGTiff, heightrows, widthcols, countbands, dtypefloat32, transformtransform, crsEPSG:32650 ) as dst: for i in range(bands): dst.write(reflectance[:, :, i], i 1)如果你输入的数据本身没有地理坐标信息那就没法直接套用这步需要你在config.py里手动指定一个原点坐标和像元尺寸或者干脆输出成不带坐标的TIFF后续在GIS软件里再手动配准。我一般建议在处理前就确认好坐标系不然后面做空间分析时还得重新投影。4. 常见问题与排查技巧实录4.1 压缩包解压失败问题可能出在下载环节这个代码包是以zip形式分发的但有不少人反馈解压时报“file is not a zip file”或“invalid zip archive: could not find eocd”。这两个报错几乎都不是代码本身的问题而是文件下载不完整或传输过程损坏。zip文件末尾有一个称为EOCDEnd Of Central Directory的固定结构如果下载过程被中断或浏览器把它当成文本文件保存了解压工具就找不到这个结束标记。解决办法很简单重新下载整个zip文件不要用下载工具的多线程加速有时候反而把文件拆坏。下载完成后先看一下文件大小跟服务端标注的是否一致。如果下载的是压缩包的一部分哪怕只差几个字节也可能解压失败。还有一种情况是文件名里有中文或特殊符号导致某些老旧解压工具识别异常建议把zip文件放到纯英文路径下再解压例如D:\hs_preprocess\。如果你发现压缩包确实完整但双击打开还是报错那就换一个解压工具。Windows自带的“资源管理器”和WinRAR的处理方式不同我用7-Zip遇到这种问题最少。至于网上常说的zip密码问题这个代码包没有加密不需要考虑。4.2 spectral库读取ENVI文件时报错通常是头文件格式不规范spectral库对ENVI头文件的解析比较严格像波长字段必须是“wavelength {400.0, 401.0, ...}”这种花括号形式有些软件生成的头文件里可能没有这个字段或者用了“wavelength units Nanometers”但数据实际是微米。如果你打开时遇到类似“unknown header parameter”的警告不要忽略有些只是警告不影响读取但有些会导致波段信息丢失。另外如果你读取的是BIP或BIL格式spectral库本身支持得很好但前提是头文件里的“interleave”字段必须写正确。一个常见的坑是原始数据明明是BSQ格式头文件里却写成了BIL读取后图像看起来全是条纹状噪声。遇到这种问题可以先打开.hdr文件看一眼interleave字段再结合数据生产者提供的说明去核对。我建议遇到任何读取异常第一步都先用文本编辑器打开.hdr文件对照spectral官方文档检查必填字段。这一步看似原始但定位问题比在Python里反复试错快得多。4.3 处理大影像时内存占用过高如何分块处理高光谱影像动辄几千乘几千像素波段数上百个如果用float32把所有数据一次性读进内存一个5GB的文件转成numpy数组直接就是爆炸性的内存占用。我跑机载高光谱数据时曾把32GB内存吃满然后被系统强制杀掉进程。好在这套代码包里的大部分函数都是逐波段或逐像元操作的可以很方便地改成按块处理。最简单的写法是这样block_rows 256 for i in range(0, rows, block_rows): block data[i:iblock_rows, :, :] # 对block做定标、平滑、导数等操作 # 将结果写入输出文件对应位置block_rows取多少取决于你机器的内存大小和数据的列数。256行对大多数情况是个安全值块太大内存容易爆块太小I/O开销增加。另一个实用技巧是如果你只需要做光谱维的处理可以先按波段读取数据也就是“band-interleaved”的方式逐波段处理这样内存占用更可控。4.4 参数设置与结果验证的几点心得预处理做完不代表万事大吉我建议至少做三件事来验证结果第一随机抽取几个像元绘制处理前后的光谱曲线用肉眼确认没有明显的突变、负值、或全零波段。正常的植被光谱应该表现出红边效应即红光波段吸收明显近红外波段反射率上升正常的矿物光谱应该在2.2-2.3微米附近有特征吸收谷。第二计算全图的均值光谱和标准差如果均值光谱在某个波段突然跳变到异常值多半是坏波段没剔除干净或定标系数写错了。第三和你已知的地面实测点做对比哪怕只有几个点也能精确评估预处理流程是否需要调整。精度上不要追求一次到位光谱预处理本来就是一个多步骤迭代调优的过程前几步做粗了后面再调也就只能在小范围内修补了。说到底这套代码包提供的是流程框架和常用算法但每个传感器的数据特性差很远拿到数据后磨参数的过程是逃不掉的。所谓“经验”其实就是把这些坑踩一遍之后下次能更快定位问题在哪。希望这包里已经帮你填掉一部分坑剩下的就靠你手上的数据和实际需求去迭代了。本文还有配套的精品资源点击获取