4、高光谱数据预处理（二）：光谱平滑与去噪（Savitzky-Golay滤波）、坏波段剔除、数据降维（MNF变换）
好，咱们接着聊高光谱数据预处理。上一章我们把辐射定标和大气校正搞定了，数据总算有了物理意义。但别高兴太早，你拿到的光谱曲线，十有八九还是像得了帕金森一样——抖得厉害。
为什么会这样？说白了，传感器本身就有噪声，再加上大气残留、仪器暗电流这些因素，光谱曲线上的毛刺和尖峰是家常便饭。我刚开始做高光谱分析时，拿到一条玉米叶片的光谱，那曲线锯齿得跟心电图似的，差点以为作物生病了。后来才发现，嗯，是噪声在捣鬼。
4.1 光谱平滑与去噪：Savitzky-Golay滤波
平滑去噪的方法很多，移动平均、中值滤波、小波去噪……但我个人最常用、也最推荐的是 Savitzky-Golay（SG）滤波。为什么？因为它能在去噪的同时，尽量保留光谱的细节特征，比如吸收谷、反射峰的位置和宽度。
SG滤波的核心思想很简单：用一个滑动窗口，对窗口内的数据点进行多项式拟合，然后用拟合值代替原始值。你想想看，移动平均只是简单求平均，会把峰谷都抹平；而SG滤波是用多项式去逼近局部曲线，相当于给每个小段配了一条「定制曲线」，保形能力好得多。
- 窗口宽度（window_length）：必须是奇数。窗口越大，平滑效果越强，但细节丢失也越多。我一般从5开始试，最多不超过21。
- 多项式阶数（polyorder）：一般取2或3。阶数太高容易过拟合，反而把噪声也拟合进去了。
我在项目中遇到过一件事：有一次处理小麦病害的高光谱数据，窗口设了15，阶数取了4，结果把病害特征吸收谷给平滑没了。后来改成窗口9、阶数2，效果就对了。所以参数调优，一定要结合你的实际数据来试。
# Python示例：Savitzky-Golay滤波
from scipy.signal import savgol_filter
import numpy as np
# 假设 spectra 是形状为 (n_samples, n_bands) 的数组
# 窗口宽度=11，多项式阶数=2
spectra_smooth = savgol_filter(spectra, window_length=11, polyorder=2, axis=1)
# 可视化对比
import matplotlib.pyplot as plt
plt.plot(wavelengths, spectra[0, :], label='原始光谱', alpha=0.6)
plt.plot(wavelengths, spectra_smooth[0, :], label='SG平滑后', linewidth=2)
plt.legend()
plt.xlabel('波长 (nm)')
plt.ylabel('反射率')
plt.show()4.2 坏波段剔除
平滑完了，接下来要干一件事——剔除坏波段。什么叫坏波段？就是那些信噪比极低、或者受水汽吸收严重影响、或者传感器本身有问题的波段。
我记得有一次处理国产某高光谱卫星的数据，发现1450nm和1950nm附近反射率直接掉到负值。一开始我还以为是大气校正没做好，折腾了半天，后来查资料才发现，这两个波段正好是水汽强吸收带，传感器在这两个位置基本是「瞎的」。嗯，这就是典型的坏波段。
坏波段剔除通常分两步：
- 基于物理知识剔除：比如1350-1450nm、1800-1950nm的水汽吸收带，以及传感器边缘波段（如首尾各10个波段）。
- 基于统计指标剔除：计算每个波段的信噪比（SNR）或方差，把明显异常的波段去掉。
| 波段范围 (nm) | 常见问题 | 处理建议 |
|---|---|---|
| 350-400 | 信噪比极低 | 通常剔除 |
| 1350-1450 | 水汽强吸收 | 必须剔除 |
| 1800-1950 | 水汽强吸收 | 必须剔除 |
| 2400-2500 | 传感器边缘，信噪比低 | 建议剔除 |
# 坏波段剔除示例
bad_bands = [0,1,2,3,4, # 前5个波段信噪比低
135,136,137,138,139,140,141,142,143,144, # 1350-1450nm
180,181,182,183,184,185,186,187,188,189,190,191,192,193,194] # 1800-1950nm
# 保留好波段
good_bands = [i for i in range(spectra.shape[1]) if i not in bad_bands]
spectra_clean = spectra[:, good_bands]
wavelengths_clean = wavelengths[good_bands]4.3 数据降维：MNF变换
好，现在数据平滑了，坏波段也剔除了。但还有一个问题——高光谱数据的维度太高了。几百个波段，你直接扔进分类模型，不仅计算慢，还容易过拟合。这就是所谓的「维度灾难」。
降维的方法很多，PCA（主成分分析）是最常见的。但做高光谱的人都知道，PCA有个毛病：它不分噪声和信号，一股脑儿全压到前几个主成分里。结果就是，前几个主成分里既有信号也有噪声，后面的主成分里全是噪声。
MNF（Minimum Noise Fraction，最小噪声分数）变换就是来解决这个问题的。它先把数据分成「信号空间」和「噪声空间」，然后只保留信号占主导的成分。说白了，MNF是PCA的升级版，专门为高光谱数据设计的。
- 估计噪声协方差矩阵：通常用空间差分法（比如计算相邻像素的差值）来估计噪声。
- 白化噪声：对数据进行变换，使噪声变成单位方差的白噪声。
- 标准PCA：在白化后的数据上做PCA，得到MNF成分。
- 选择有效成分：保留特征值大于1的成分（因为噪声的特征值理论上为1）。
我个人习惯用MNF而不是PCA，尤其是在做作物精细分类时。有一次做水稻品种识别，用PCA降维后分类精度只有78%，换成MNF后直接跳到89%。为什么？因为MNF把噪声分离出去了，模型学到的都是「干净」的信号。
# MNF变换示例（使用spectral库）
from spectral import open_image, MNF
# 读取高光谱数据
img = open_image('hyperspectral_image.hdr')
data = img.load()
# 执行MNF变换
mnf = MNF(data)
mnf_data = mnf.reduce() # 自动选择有效成分
# 查看特征值
print('特征值:', mnf.eigenvalues)
# 手动选择前N个成分
n_components = 10 # 根据特征值>1的原则选择
mnf_selected = mnf_data[:, :, :n_components]4.4 本章知识体系
说了这么多，咱们用一张图把本章的核心逻辑串起来。你想想看，整个预处理流程就像一条流水线：原始数据进来，先平滑去噪，再剔除坏波段，最后降维压缩。每一步都有它的道理，缺一不可。
好了，这一章的内容就到这里。SG滤波怎么调参、坏波段怎么判断、MNF为什么比PCA更适合高光谱，这些我都结合自己的经验讲清楚了。你回去拿自己的数据跑一遍，遇到问题随时回来翻翻这一章。