4. Savitzky-Golay滤波：多项式拟合原理、窗口与阶数选择、导数光谱计算

好，咱们进入第四章。这一章聊的是光谱预处理里一个非常经典的工具——Savitzky-Golay滤波，简称S-G滤波。说实话，我在做近红外光谱分析的头两年，几乎每个项目都会用到它。它既能平滑降噪，又能保留峰形细节，甚至还能直接算导数光谱，一举多得。

但很多人用S-G滤波就是调个参数，窗口大小和多项式阶数随便填，结果要么过平滑把峰磨没了，要么欠平滑噪声还在。今天我就把这里面的门道掰开揉碎了讲清楚。

4.1 多项式拟合原理：滑动窗口里的最小二乘

S-G滤波的核心思想其实很简单：在光谱上开一个滑动窗口，对窗口内的数据点做多项式最小二乘拟合，然后用拟合值替换窗口中心点的值。

你想想看，传统的移动平均滤波是把窗口内所有点取平均，这相当于用一个常数去拟合窗口内的数据。但光谱信号往往是有趋势的，用常数拟合会丢失细节。S-G滤波用多项式拟合，就能更好地保留峰的形状。

具体来说，假设窗口宽度为 \(2m+1\) 个点（m是半窗宽），窗口内的数据点为 \(x_{-m}, x_{-m+1}, ..., x_0, ..., x_{m-1}, x_m\)。我们用一个k阶多项式来拟合这些点：

\[ y_i = a_0 + a_1 i + a_2 i^2 + ... + a_k i^k \]

其中i是相对于窗口中心的偏移量（-m到m）。通过最小二乘法求出系数 \(a_0, a_1, ..., a_k\)，然后取 \(y_0 = a_0\) 作为滤波后的值。

我在项目中遇到过一个问题：有人把S-G滤波当黑盒用，结果发现光谱的峰位偏移了。后来一查，是因为窗口不对称导致的边界效应。嗯，这里要注意，S-G滤波对边界点的处理方式不同，有的库是补零，有的是镜像，有的是截断。我建议你使用前先确认一下边界处理方式。

4.2 窗口与阶数选择：参数调优的实战经验

窗口大小和多项式阶数是S-G滤波的两个核心参数。选不好，效果天差地别。我给大家总结一个经验法则：

参数	选小了会怎样	选大了会怎样	我的建议
窗口大小	降噪效果差，噪声残留多	过平滑，峰变矮变宽，细节丢失	窗口宽度应大于峰的半高宽对应的点数
多项式阶数	拟合能力不足，平滑过度	过拟合，噪声被保留甚至放大	阶数一般取2~4，不要超过5

为什么会这样？因为窗口大小决定了平滑的「力度」，阶数决定了拟合的「灵活度」。两者需要平衡。我个人习惯用3阶多项式配合窗口宽度为峰半高宽的1.5~2倍。比如一个峰的半高宽是20个数据点，那窗口宽度我选30~40。

还有一个实用技巧：如果你不确定参数，可以画一条「参数扫描曲线」。固定阶数，改变窗口大小，观察信噪比和峰高的变化。选一个拐点处的参数，往往就是最优解。

4.3 导数光谱计算：一阶导、二阶导的S-G实现

S-G滤波的一个杀手锏是：它不仅能平滑，还能直接计算导数光谱。而且这个导数是在平滑的同时计算的，比先平滑再差分要稳定得多。

原理很简单：多项式拟合得到系数后，导数就是多项式的导数。比如一阶导在窗口中心的值就是 \(a_1\)，二阶导是 \(2a_2\)。所以S-G滤波的卷积核，其实可以同时输出平滑值、一阶导、二阶导……只要你想。

导数光谱有什么用？

• 一阶导：消除基线漂移（常数基线被完全去除，线性基线被转化为常数偏移）
• 二阶导：消除线性基线，同时增强峰的分辨率（重叠峰可以被分开）

但注意，导数阶数越高，噪声放大越严重。我一般只用一阶导或二阶导，三阶以上很少用，除非数据信噪比极高。

4.4 实战案例：用Python实现S-G滤波与导数计算

好了，理论讲完了，咱们直接上代码。我用的是scipy.signal里的savgol_filter函数，它封装得很好，一行代码搞定。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter

# 生成模拟光谱：两个高斯峰 + 噪声 + 基线漂移
x = np.linspace(0, 100, 500)
peak1 = 5 * np.exp(-((x - 30) / 5)**2)
peak2 = 3 * np.exp(-((x - 60) / 4)**2)
baseline = 0.02 * x + 0.5  # 线性基线
noise = np.random.normal(0, 0.2, len(x))
y = peak1 + peak2 + baseline + noise

# S-G滤波：窗口31点，3阶多项式
y_sg = savgol_filter(y, window_length=31, polyorder=3)

# 一阶导：deriv=1
y_sg_1st = savgol_filter(y, window_length=31, polyorder=3, deriv=1, delta=x[1]-x[0])

# 二阶导：deriv=2
y_sg_2nd = savgol_filter(y, window_length=31, polyorder=3, deriv=2, delta=x[1]-x[0])

# 画图
plt.figure(figsize=(12, 8))
plt.subplot(3, 1, 1)
plt.plot(x, y, label='原始光谱（含噪声+基线）', alpha=0.5)
plt.plot(x, y_sg, label='S-G滤波后', linewidth=2)
plt.legend()
plt.title('S-G滤波效果')

plt.subplot(3, 1, 2)
plt.plot(x, y_sg_1st, label='一阶导', color='green')
plt.axhline(0, color='gray', linestyle='--')
plt.legend()
plt.title('一阶导数光谱')

plt.subplot(3, 1, 3)
plt.plot(x, y_sg_2nd, label='二阶导', color='red')
plt.axhline(0, color='gray', linestyle='--')
plt.legend()
plt.title('二阶导数光谱')

plt.tight_layout()
plt.show()

运行这段代码，你会看到：原始光谱有明显的基线漂移和噪声，S-G滤波后基线被保留但噪声被抑制。一阶导光谱中，基线变成了一个常数偏移（因为线性基线的导数是常数），两个峰变成了正负交替的波形。二阶导光谱中，基线几乎被完全消除，两个峰变成了向下的尖峰，峰位非常清晰。

我个人习惯在项目里先做一次S-G滤波（窗口小一点，比如15~21点），然后再做一次导数计算。这样既能看清峰位，又能保留定量信息。当然，如果你只需要定性分析，直接上二阶导也行。

最后，我用一张流程图来总结S-G滤波的核心逻辑，方便你理解整个流程：

好了，这一章的内容就到这里。S-G滤波是光谱预处理里最实用的工具之一，掌握了它，你就能在降噪和保留细节之间找到平衡点。下一章我们会聊另一种常用的降噪方法——小波变换，它和S-G滤波的思路不太一样，但同样强大。

---

公众号:蓝海资料掘金营，微信deep3321