简介面向人工智能与机器学习场景下的高光谱数据预处理需求这份Python代码包提供了可直接调用的实践方案适用于遥感、农业光谱分析等任务的建模前准备。压缩包共15个文件整体仅2.48MB包括两个核心Python脚本预处理函数与演示程序、一个高光谱csv样例数据以及12张过程与结果对比图便于拆解每一步处理效果。代码覆盖光谱校正、去噪、平滑、光谱指数计算、特征选择、异常检测、标准化与降维等常见环节适合有Python基础、从事遥感或光谱数据分析的算法工程师与研究者参考学习。目前已有1375人学习下载作者为admin_maxin。通过阅读两个脚本的调用逻辑可快速掌握PCA、z-score等在高光谱数据上的落地写法并配合图示检查效果直接迁移到自身项目中为后续分类或回归任务提供更干净的输入特征。1. 高光谱数据预处理为什么说拿到手的第一件事不是跑模型高光谱数据预处理的Python代码本质是把遥感影像从“机器读数”翻译成“物理量”的一套标准化流程。业内常说要“先定标、再去噪、后校正、最后降维”这四步决定了后续分类或定量反演的天花板。我见过太多人拿到高光谱数据直接丢进分类器结果精度虚高换一块区域或换一景影像就全线崩盘问题几乎都出在预处理环节辐射定标系数没查、坏波段没剔除、大气校正参数填错。这篇笔记就按一个可复现的Python预处理流程展开从ENVI元数据读取到PCA降维把参数、坑点和验证方法都讲清楚适合遥感专业学生、刚转行的算法工程师以及被高光谱数据折磨过的科研人员。2. 辐射定标与坏波段剔除预处理的第一步是把单位搞对2.1 为什么必须做辐射定标DN值不是反射率高光谱传感器记录的是原始DN值Digital Number它受传感器增益、偏移、暗电流和光照条件共同影响。以推扫式成像光谱仪为例同一地物在不同扫描行、不同探测元件上的DN值都可能不一致这是传感器硬件特性决定的。如果不做辐射定标后续的反射率反演、植被指数计算、光谱角匹配全部建立在错误的数值基础上等于把误差放大了好几倍。辐射定标的数学表达式是L gain × DN offset其中gain增益和offset偏移通常可以从影像自带的元数据文件或头文件里读取。高光谱数据常见格式为ENVI Standard.dat .hdr头文件里会记录calibrate gain和calibrate offset这两个关键字段但不同传感器命名可能不同有的叫scale_factor有的叫reflectance_coefficients需要具体看hdr里的Data_Ignore_Value和Band_Names。2.2 用Python实现辐射定标最小可运行代码下面是辐射定标的核心代码按标准流程读取头文件参数并转换为辐射亮度import numpy as np from osgeo import gdal def read_envi_hdr(hdr_path): 解析ENVI头文件返回参数字典 params {} with open(hdr_path, r, encodingutf-8, errorsignore) as f: for line in f: if in line: key, value line.strip().split(, 1) params[key.strip()] value.strip() return params def radiometric_calibration(dn_img_path, hdr_path, out_path): DN值转辐射亮度 gain和offset从hdr中读取如果没有则使用默认值1和0 hdr read_envi_hdr(hdr_path) ds gdal.Open(dn_img_path) img ds.ReadAsArray().astype(np.float32) # 读取增益和偏移没有就默认1和0相当于不做定标 try: gain float(hdr.get(calibrate_gain, 1.0)) offset float(hdr.get(calibrate_offset, 0.0)) except ValueError: gain, offset 1.0, 0.0 radiance img * gain offset driver gdal.GetDriverByName(ENVI) out_ds driver.Create(out_path, ds.RasterXSize, ds.RasterYSize, ds.RasterCount, gdal.GDT_Float32) out_ds.WriteArray(radiance) # 复制地理参考信息 out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_ds.FlushCache() del ds, out_ds return radiance这段代码先把所有波段读入内存转换为float32类型防止溢出然后逐波段执行L gain × DN offset。注意ReadAsArray不带波段参数时返回三维数组形状是波段数, 行数, 列数如果影像太大内存不够需要分块处理用ReadAsArray(col_off, row_off, col_size, row_size)循环读写。写完数据后务必检查输出值域辐射亮度的典型数值范围在0~100之间如果出现几千上万的值说明gain和offset的单位理解错了很可能是把整数增益当成浮点增益用了需要回看hdr里的Units字段。2.3 坏波段识别信噪比低到一定程度就该放弃高光谱数据不是每个波段都能用。以水汽吸收波段为例1350~1450nm和1800~1950nm区间的大气透过率几乎为零传感器记录的基本是噪声另外传感器的首尾波段常因探测器响应不稳而质量很差。把这些波段保留在数据里会显著干扰归一化处理和后续的降维结果。识别坏波段最稳妥的方法是计算每个波段的信噪比SNR筛选出SNR低于阈值的波段。计算方法是用影像中均匀地物的均值除以标准差def bad_band_detection(img_array, threshold10.0): 通过SNR识别坏波段 img_array形状为(波段数, 行数, 列数) 返回坏波段索引列表和每个波段的SNR值 bands img_array.shape[0] snr_values [] bad_bands [] for i in range(bands): band_data img_array[i] # 均匀区域的标准差代表噪声水平 mean_val np.mean(band_data) std_val np.std(band_data) if std_val 0: snr float(inf) else: snr mean_val / std_val snr_values.append(snr) if snr threshold: bad_bands.append(i) return bad_bands, snr_values这段代码有个隐含假设整个影像的均值和标准差能代表地物信号与噪声的比值。实际情况中影像里如果有大面积水体或阴影标准差会被拉高SNR被低估。更稳妥的做法是从影像中心裁剪一个纯地物区域来计算SNR或者在确定好波段后做一次目视检查。我一般会把每个波段的均值和标准差画成折线图坏波段的曲线会出现明显的尖峰或断崖比阈值判断更直观。坏波段剔除不是简单地删掉索引要考虑波段顺序是否影响后续处理。如果做光谱重采样或融合删除波段后需要同步更新波长数组如果直接做分类保持剩余波段的连续性和波长对应关系即可。3. 去噪处理的三种主流方法MNF变换与Savitzky-Golay滤波3.1 去噪在高光谱预处理中的位置不是越干净越好高光谱影像的噪声来源复杂包括探测器暗电流噪声、散粒噪声、条纹噪声坏线以及大气校正残留的噪声。去噪做得好能显著提升分类精度做得过度会把细微的光谱吸收特征抹掉反而降低地物识别能力。我在实际项目中见过有人把光谱曲线磨得像镜面一样光滑结果矿物填图时方解石和白云石的光谱差异完全分不开了这就是过度平滑的代价。主流的去噪方法分两类一类是变换域方法把影像变换到特征空间再截断噪声分量代表是MNF最小噪声分离变换对高光谱数据效果最好另一类是平滑滤波方法在光谱维度上做卷积平滑代表是Savitzky-GolaySG滤波计算快、参数直观适合批量处理。3.2 MNF变换去噪用噪声协方差矩阵分离信号与噪声MNF变换的原理是先估计影像的噪声协方差矩阵然后做两次PCA第一次对噪声白化第二次在去除噪声相关的空间中做标准PCA。经过MNF变换后信号集中在前面几个分量里噪声分散在后面分量里通过截取前若干个MNF分量再逆变换即可实现去噪。Python实现MNF变换需要借助scikit-learn的PCA类from sklearn.decomposition import PCA def mnf_denoise(img_array, n_components30): MNF变换去噪 img_array: 形状为(bands, rows, cols)的影像数组 n_components: 保留的MNF分量数一般取总波段的20%~30% bands, rows, cols img_array.shape # 把影像展平为(像素数, 波段数)矩阵 pixels img_array.reshape(bands, -1).T # shape: (rows*cols, bands) # 计算噪声协方差矩阵用相邻像元差分的均值代表噪声 # 对每个像元计算与右/下相邻像元的差 diff_right np.diff(img_array, axis2) # 沿列方向差分 diff_down np.diff(img_array, axis1) # 沿行方向差分 noise np.concatenate([ diff_right.reshape(bands, -1).T, diff_down.reshape(bands, -1).T ], axis0) # 对噪声做PCA得到白化矩阵 noise_pca PCA() noise_pca.fit(noise) white_mat noise_pca.components_ # 噪声主成分方向 # 白化数据投影到噪声主成分空间并归一化 centered pixels - np.mean(pixels, axis0) whitened centered white_mat.T / np.sqrt(noise_pca.explained_variance_ 1e-10) # 在白化空间中做PCA取前n_components个分量 signal_pca PCA(n_componentsn_components) transformed signal_pca.fit_transform(whitened) # 逆变换回原空间 restored signal_pca.inverse_transform(transformed) restored restored * np.sqrt(noise_pca.explained_variance_ 1e-10) restored restored white_mat np.mean(pixels, axis0) # 还原为影像形状 denoised restored.T.reshape(bands, rows, cols) return denoised这个实现有几个关键细节噪声估计用了相邻像元的差分这是高光谱去噪的常见做法假设真实地物在局部区域变化平缓差分结果主要代表噪声白化后数据的PCA本质上是标准化的MNF变换。n_components的取值需要凭经验调试取太少会丢失弱信号取太多则去噪不彻底。MNF去噪的缺点是计算量大band×pixel的矩阵做特征分解在数据量大时非常吃内存。我处理200波段、512×512影像时大概要用8GB内存如果影像更大建议分块处理或直接改用下一种方案。3.3 Savitzky-Golay滤波光谱维平滑的轻量方案SG滤波是一种局部多项式拟合的平滑方法它不像滑动平均那样简单取窗口均值而是在滑动窗口内做多项式最小二乘拟合用拟合值替代中心点。这样做的好处是能保持光谱曲线的峰谷形状不会把吸收特征直接抹平。scipy.signal库直接提供了savgol_filter函数参数只有两个窗口长度和多项式阶数from scipy.signal import savgol_filter def sg_denoise(img_array, window_length11, polyorder3): SG滤波去噪 逐像元在光谱维度上平滑 window_length必须是奇数一般取值5~15 polyorder必须小于window_length一般取2~4 bands, rows, cols img_array.shape # 展平为(像素数, 波段数) pixels img_array.reshape(bands, -1).T # 对每个像元的光谱曲线做SG平滑 denoised_pixels np.apply_along_axis( lambda x: savgol_filter(x, window_length, polyorder), axis1, arrpixels ) # 还原影像形状 denoised denoised_pixels.T.reshape(bands, rows, cols) return denoised这段代码写成np.apply_along_axis行数太大时效率较低。更高效的做法是直接对三维数组沿波段轴做卷积但scipy的sg接口不支持多维直接处理所以我通常用循环逐像元处理在代码可读性和性能之间取平衡。SG滤波的参数选择有规律可循窗口越大平滑程度越高但光谱细节丢失越多多项式阶数越高越能拟合尖锐特征但阶数过高会导致过拟合噪声。对于一般的高光谱影像我习惯先用window_length11, polyorder3跑一遍观察典型地物光谱曲线的变化再根据保真度调整。3.4 两种方法怎么选数据量大用SG精度要求高用MNF实际项目中MNF和SG的选择取决于数据规模和应用场景。MNF去噪更彻底能处理条纹噪声和随机噪声的混合体适合矿物填图这种对光谱细节要求高的任务SG滤波速度快、参数少适合批量预处理的中间环节尤其适合用于后续需要计算光谱指数的场景。我个人的处理习惯是先用SG滤波做轻量平滑再做大气校正最后用MNF去除大气校正引入的残留噪声。这样能把两种方法的优势叠加又不至于过度处理。4. 大气校正的工程化落地从辐射亮度到地表反射率4.1 高光谱数据为什么必须做大气校正辐射定标之后的数据仍然是“表观反射率”它包含了大气分子散射、气溶胶散射、水汽吸收对辐射信号的衰减和增强。如果不做大气校正直接拿表观反射率去建谱库到另一个季节或另一个大气条件下同样地物的光谱曲线可能差出20%以上分类器或光谱匹配算法的结果就不具普适性。高光谱大气校正的常用方法是辐射传输模型法代表是6SSecond Simulation of a Satellite Signal in the Solar Spectrum它在ENVI里的实现是FLAASH模块。Python生态里没有现成的6S封装但可以通过几个替代方案实现一是调用Py6S库封装6S模型二是用ENVI的FLAASH模型独立运行三是针对特定传感器使用简化校正公式。4.2 用Py6S调用辐射传输模型基础参数怎么填Py6S是6S模型的Python封装安装后可以构建大气参数并逐波段计算校正系数from Py6S import SixS, Parameters, Wavelength, AtmosProfile, AeroProfile def atmospheric_correction_py6s(wavelengths, atmosphere_typemidlatitude): 计算每个波段的辐射传输参数 返回每个波段的透射率、大气路径辐射等 s SixS() # 设置大气模式 if atmosphere_type midlatitude: s.atmos_profile AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer) elif atmosphere_type tropical: s.atmos_profile AtmosProfile.PredefinedType(AtmosProfile.Tropical) # 设置气溶胶模式大陆型、乡村型、海洋型等 s.aero_profile AeroProfile.PredefinedType(AeroProfile.Continental) # 设置气溶胶光学厚度AOT没有实测值就用0.2 s.aot550 0.2 # 设置传感器高度和观测几何 s.sensor_z 705 # 以千米为单位的星下点高度 s.satellite_altitude 705 s.geometry Geometry.User() s.geometry.solar_z 30 # 太阳天顶角 s.geometry.solar_a 190 # 太阳方位角 s.geometry.view_z 0 # 观测天顶角 s.geometry.view_a 0 # 观测方位角 correction_params [] for wl in wavelengths: s.wavelength Wavelength(wl) # 单位纳米 s.run() # 6S输出transmittance、大气路径辐射等 params { transmittance: s.outputs.transmittance_total, path_radiance: s.outputs.atmospheric_intrinsic_radiance, downward_flux: s.outputs.downward_flux } correction_params.append(params) return correction_paramsPy6S的输出参数里最关键的有三个大气总透射率transmittance_total、大气路径辐射atmospheric_intrinsic_radiance和向下辐射通量。地表反射率的计算公式是ρ (L_sat - L_path) / (T_up × F_down / π T_down × ρ_prev × S)这一套公式里没有完全标定好的传感器参数时直接把Py6S输出的透射率、路径辐射用一个线性回归模型拟合到地面实测反射率上是工程上更省事的做法。如果区域内有实测光谱点用这个思路做精度远高于纯模型参数。4.3 FLAASH大气校正的调用参数与输入文件准备ENVI的FLAASH是工程上最常用的高光谱大气校正工具它需要的输入参数包括影像中心经纬度、传感器高度、飞行时间、大气模式、气溶胶模式、气溶胶反演方法、初始能见度或气溶胶光学厚度。这些参数设置得当与否直接决定校正效果。关键参数的典型取值如下表参数典型取值说明大气模式Mid-Latitude Summer / Tropical根据纬度和季节选择可查FLAASH手册的大气模式表水汽反演开启高光谱数据必须开启水汽反演用940nm和1130nm吸收带气溶胶模式Urban / Rural城市区域选Urban郊野选Rural初始能见度40km如果怀疑大气浑浊可调小到20km光谱平滑关闭只做气溶胶校正时关闭避免过度平滑调用FLAASH有两种方式一是在ENVI图形界面里操作二是准备好输入文件和参数后用ENVI的批处理功能自动跑。批处理模式下输入文件必须是BIL或BSQ格式的浮点辐射亮度数据单位W/m²/μm/sr先转成ENVI标准格式再调用FLAASH。# 使用ENVI的FLAASH批处理接口通过命令行或ENVI服务沙箱调用 envi_flaash_engine \ -input radiance.dat \ -input_header radiance.hdr \ -output reflectance.dat \ -lat 39.9 -lon 116.4 \ -sensor_height 705 \ -flight_date 2024-06-15 \ -flight_time_gmt 03:30:00 \ -atmos_model Mid-Latitude Summer \ -aerosol_model Rural \ -aot 0.2 \ -water_retrieval true温度、湿度对高光谱大气校正的影响不可忽视。夏季潮湿大气的水汽含量高940nm附近的吸收很强不开启水汽反演校正结果会明显偏暗。我遇到过几次因为忘记设置飞行时间导致水汽反演失败的状况排查了两天才发现是UTC和本地时间换算出了问题。4.4 大气校正后的验证用光谱曲线形状做快速判断大气校正是否成功可以在校正后的影像上找几个已知地物做光谱形状验证。植被的典型光谱曲线应该是“绿峰-红边-近红外高原-红边下降”的形态水体在近红外波段应该是低反射率裸土在短波红外波段有缓慢上升趋势。如果植被在红波段680nm的反射率高于绿波段550nm说明大气校正的路径辐射扣除得不够彻底如果植被的近红外反射率异常低可能是水汽反演失败或气溶胶光学厚度设置过大。用Python可视化验证光谱曲线是一个好习惯import matplotlib.pyplot as plt def plot_spectra(img_array, wavelengths, pixel_coords): 可视化典型地物的光谱曲线 pixel_coords: 列表每个元素是(行, 列)坐标 bands img_array.shape[0] plt.figure(figsize(10, 6)) for idx, (row, col) in enumerate(pixel_coords): spectra img_array[:, row, col] plt.plot(wavelengths, spectra, labelfPixel {idx1}: ({row},{col})) plt.xlabel(Wavelength (nm)) plt.ylabel(Reflectance) plt.legend() plt.title(Spectral Curves After Atmospheric Correction) plt.grid(True) plt.show()这段代码特别适合在大气校正前后对比同一坐标的像元光谱。校正前和校正后曲线形状的差异如果不明显就要怀疑参数设置是否真的生效了如果曲线整体被平移说明路径辐射校正力度过大或不足。5. 高光谱预处理避坑指南四个反复翻车的场景5.1 数据格式陷阱BSQ/BIL/BIP混用导致波段顺序错乱现象读完数据后做坏波段剔除发现剔除的波段对不上号或者某个波段看起来特别模糊。原因高光谱数据常见的三种存储格式——BSQ波段顺序存储、BIL波段交叉存储、BIP像元交叉存储——在重新组织数据时会直接影响内存中的波段顺序。用gdal打开ENVI格式文件时gdal会自动处理格式转换但用numpy的fromfile直接读二进制文件时不会一旦数据换了一种格式就得手动重排波段轴。解决处理ENVI格式的二进制数据时优先使用gdal或spectral库来打开不要直接读.dat裸二进制。spectral库的open_image接口能正确读取ENVI头文件并返回band-interleaved格式的数组省去格式转换的麻烦。5.2 内存爆炸问题大影像预处理时带宽分配不当现象在8GB内存的电脑上跑MNF变换程序直接闪退或报MemoryError。原因高光谱影像动辄上百波段如果影像尺寸是1000×1000一个float32数组就占400MBMNF变换要把影像展平成100万×200的矩阵做特征分解内存直接翻了好几倍。解决分块处理是唯一靠谱的方案。用gdal的ReadAsArray按行块读取每处理完一个块就释放内存同时对中间变量用del和gc.collect()及时清理。阈值要留到实际数据的1.5倍以上否则Python的垃圾回收跟不上数组创建和销毁的速度。5.3 波长单位不统一nm和μm混用导致索引偏移现象坏波段剔除后光谱曲线的波峰位置比标准谱库偏移了20nm左右。原因ENVI头文件里波长单位有时是nm有时是μm如果代码里默认波长单位是nm遇到μm的数据会把所有波长值放大1000倍导致波段索引错位。另一种情况是FLAASH输入时需要μm而辐射定标输出是nm中间少了换算。解决在读取头文件后立即统一波长单位写成一个小函数做强制转换def normalize_wavelength_unit(wavelengths, unit_str): 统一波长单位为nanometer unit_str可能是Nanometers、Micrometers、nm、μm等 if unit_str.lower() in [micrometers, μm, um, microns]: return np.array(wavelengths) * 1000.0 elif unit_str.lower() in [nanometers, nm]: return np.array(wavelengths) else: raise ValueError(fUnknown wavelength unit: {unit_str})这个函数虽然简单但能避免90%的波长错位问题。5.4 FLAASH输入数据的单位坑辐射亮度单位换算失误现象FLAASH校正结果整体偏亮或偏暗而且所有地物的反射率变化规律相同。原因FLAASH要求输入辐射亮度单位是W/m²/μm/sr但有些传感器输出的是W/m²/nm/sr两者之间差1000倍。如果用了错误的单位大气校正算法会认为图像异常亮或暗从而自动调整模型参数导致整体偏向。解决在做辐射定标时明确单位换算因子。如果原始数据单位是W/m²/nm/sr在写入FLAASH输入文件前乘以0.001转换成μm单位。这个换算因子建议写在代码注释里方便同事复用代码时不踩坑。6. 数据降维与特征提取PCA的波段选择实战6.1 PCA降维后保留多少个分量用累计方差贡献率做判断标准高光谱数据在预处理完成后仍然保持着上百个波段的维度如果直接作为分类器的输入不仅计算量巨大还会受到维数灾难的困扰。PCA降维是高光谱数据预处理链条的最后一环它的核心价值是把高度相关的波段压缩成少数几个主成分保留数据的主要方差信息。主成分数量的选择不能拍脑袋。我习惯用累计方差贡献率达到99%这个标准来确定分量数同时观察特征值的陡坎位置取“肘部”之前的分量。下图是判断逻辑的代码from sklearn.decomposition import PCA import numpy as np def pca_dimension_reduction(img_array, var_threshold0.99): PCA降维根据累计方差贡献率自动选择分量数 img_array: 形状(bands, rows, cols) 返回降维后的影像数组和分量数 bands, rows, cols img_array.shape pixels img_array.reshape(bands, -1).T # 先做一次完整PCA计算累计方差 pca_full PCA() pca_full.fit(pixels) cum_ratio np.cumsum(pca_full.explained_variance_ratio_) n_components np.argmax(cum_ratio var_threshold) 1 # 检查拐点如果拐点出现在前面就取拐点 diffs np.diff(pca_full.explained_variance_ratio_) elbow_candidates np.where(diffs[:-1] / (diffs[1:] 1e-10) 5)[0] if len(elbow_candidates) 0: elbow elbow_candidates[0] 1 n_components min(n_components, elbow) # 用选定的分量数重新做PCA pca_final PCA(n_componentsn_components) reduced pca_final.fit_transform(pixels) # 还原为三维影像形状 rows_cols (rows, cols) reduced_img reduced.T.reshape(-1, rows_cols[0], rows_cols[1]) return reduced_img, n_components, pca_final这段代码做了两件事先通过累计方差贡献率找到满足阈值的最小分量数再通过特征值陡坎判断是否存在更激进的降维空间。高光谱地物类别多时可以多保留几个分量作为冗余分类时再靠正则化来抑制噪声。6.2 加载PCA后的主成分影像验证降维效果降维后的影像在ENVI或Python里查看时前三个主分量通常可以用RGB合成的方式显示视觉效果比原始波段自然得多。如果前三个分量的合成图有明显的地物边界说明降维保留了空间结构信息如果图像看起来像随机噪声说明预处理环节出了问题多半是大气校正失败或坏波段剔除不彻底。我习惯在PCA降维后做一个小实验用降维数据训练一个随机森林分类器比较全波段数据与降维数据的分类精度差异。通常降维后的精度不会显著下降甚至略微提升同时训练时间缩短一个数量级这就是降维策略有效的证据。6.3 进阶用特征重要性反向验证预处理质量Python的随机森林分类器可以输出特征重要性把重要性最高的前30个波段与已知的矿物或植被吸收特征波段做对比。如果重要波段集中在已知的地物诊断特征附近说明预处理保持了光谱的真实性如果重要性急剧偏高或偏低没有跟诊断特征对应多半是坏波段去除过度或平滑过度。这一步其实是在用结果验证过程的合理性。我的血泪经验是永远不要跳过验证步骤直接进入建模等模型效果不好再回头查预处理排查成本会高出好几倍。每次处理完数据都顺手做一次光谱曲线可视化、一次特征波段合理性检查养成了习惯后面翻车的概率能减少八成。最后给一个我在项目中坚持的细节把所有预处理的中间产物辐射亮度、去噪数据、反射率数据、PCA主成分数据的保存格式和参数记录到一个文本文件里包括hdr文件中追加的处理历史字段。这样每次换传感器或换数据源时只需要对照历史记录调整参数即可不用从零开始摸索。希望帮到你。本文还有配套的精品资源点击获取