3. 时间序列重构：Savitzky-Golay滤波、Whittaker平滑、时间序列插值方法
好，咱们进入第三章。说实话，时间序列重构这块，是我在遥感数据处理中花时间最多的地方之一。你拿到的NDVI、EVI这些时序数据，看着挺规整，其实里面全是坑——云遮挡、大气扰动、传感器角度变化……说白了，原始数据就像一张被揉皱了的纸，咱们得想办法把它熨平，同时还得保留有用的细节。
这一章，我带你过三种最常用的重构方法：Savitzky-Golay滤波、Whittaker平滑，以及时间序列插值。每种方法我都踩过坑，咱们边聊边避雷。
3.1 Savitzky-Golay滤波：局部拟合的智慧
先说说Savitzky-Golay滤波，简称S-G滤波。这玩意儿最早是化学领域搞光谱分析的，后来被遥感界“借”过来处理时序数据。它的核心思想很简单：用一个滑动窗口，在每个窗口内做多项式拟合，然后用拟合值代替原始值。
嗯，这里要注意：S-G滤波有两个关键参数——窗口大小和多项式阶数。窗口越大，平滑效果越强，但细节丢失也越多；多项式阶数越高，对局部变化的拟合能力越强，但过拟合风险也越大。
我个人习惯：处理MODIS NDVI时序数据时，窗口大小通常设为5-7个时间点，多项式阶数选2或3。阶数超过3之后，效果提升不明显，反而容易引入振荡。
来看一段Python代码，我用的是scipy里的savgol_filter函数：
import numpy as np
from scipy.signal import savgol_filter
import matplotlib.pyplot as plt
# 模拟一段NDVI时序数据（含噪声）
np.random.seed(42)
t = np.arange(1, 366, 8) # 约46个时间点
ndvi_true = 0.3 + 0.4 * np.sin(2 * np.pi * t / 365) # 理想曲线
ndvi_noisy = ndvi_true + np.random.normal(0, 0.05, len(t)) # 加噪声
# S-G滤波
window_size = 7 # 窗口大小，必须为奇数
polyorder = 2 # 多项式阶数
ndvi_sg = savgol_filter(ndvi_noisy, window_size, polyorder)
# 可视化
plt.figure(figsize=(10, 5))
plt.plot(t, ndvi_true, 'g-', label='True', linewidth=2)
plt.plot(t, ndvi_noisy, 'r.', label='Noisy', alpha=0.5)
plt.plot(t, ndvi_sg, 'b-', label='S-G filtered', linewidth=2)
plt.legend()
plt.xlabel('Day of Year')
plt.ylabel('NDVI')
plt.title('Savitzky-Golay Filter Example')
plt.show()这段代码跑完之后，你会看到蓝色曲线基本跟绿色曲线重合，噪声被有效抑制了。但注意，如果窗口开太大（比如15以上），蓝色曲线会变得过于平滑，连植被生长的季节性拐点都被抹平了——这就是过平滑。
我曾经踩过的坑：有一次处理高时间分辨率数据（每天一个点），我习惯性地用了窗口大小5，结果发现滤波后的曲线跟原始数据几乎没区别。后来才意识到，窗口太小根本起不到平滑作用。对于日尺度数据，窗口至少得11-15才行。记住：窗口大小要跟数据的时间分辨率匹配。
3.2 Whittaker平滑：兼顾保真度与平滑度
Whittaker平滑，也叫“惩罚最小二乘平滑”。它跟S-G滤波的思路不太一样——S-G是局部拟合，Whittaker是全局优化。它通过一个目标函数来平衡两个目标：
- 保真度：平滑后的曲线要尽量接近原始数据
- 平滑度：平滑后的曲线要尽量光滑（用二阶差分来度量）
目标函数长这样：
Q = Σ(y_i - z_i)² + λ Σ(Δ²z_i)²
其中y_i是原始值，z_i是平滑后的值，λ是平滑参数。λ越大，曲线越平滑；λ越小，越忠实于原始数据。
说白了，Whittaker平滑就是让你在“跟着数据走”和“把曲线拉直”之间做个权衡。这个λ参数，我建议你从0.5开始试，然后根据效果调整。
我的经验：处理植被物候数据时，λ取1-5之间效果比较好。如果数据噪声特别大（比如热带雨林区云遮挡严重），λ可以调到10以上。但别超过100，否则曲线会变成一条直线，那就失去意义了。
Python实现Whittaker平滑，我一般自己写个函数，或者用现成的库：
def whittaker_smooth(y, lambda_=2, d=2):
"""
Whittaker平滑
y: 原始序列
lambda_: 平滑参数
d: 差分阶数（通常用2阶）
"""
m = len(y)
E = np.eye(m)
D = np.diff(E, n=d, axis=0)
w = np.ones(m)
W = np.diag(w)
z = np.linalg.solve(W + lambda_ * D.T @ D, W @ y)
return z
# 使用示例
ndvi_whit = whittaker_smooth(ndvi_noisy, lambda_=2)
plt.plot(t, ndvi_noisy, 'r.', label='Noisy')
plt.plot(t, ndvi_whit, 'b-', label='Whittaker', linewidth=2)
plt.legend()
plt.show()这段代码里，我用了2阶差分（d=2），这是最常用的设置。1阶差分会让曲线过于“直”，3阶以上计算量增大但效果提升有限。
3.3 时间序列插值：填补缺失值
说完平滑，咱们聊聊插值。遥感时序数据里，缺失值太常见了——云遮挡、传感器故障、数据下载失败……你想想看，如果直接拿有空洞的数据去做趋势分析，结果肯定不靠谱。
插值方法有很多，我挑三种最常用的：
| 方法 | 原理 | 适用场景 | 注意事项 |
|---|---|---|---|
| 线性插值 | 用相邻两个有效值的连线填补中间 | 缺失段较短（1-2个点） | 会破坏曲线的平滑性 |
| 三次样条插值 | 用分段三次多项式拟合 | 缺失段中等长度（3-5个点） | 边界处可能振荡 |
| 时间序列分解插值 | 先分解出趋势+季节成分，再填补 | 缺失段较长（>5个点） | 计算量较大 |
我个人最常用的是三次样条插值，因为它能在保持平滑的同时，较好地拟合植被生长的非线性变化。来看代码：
from scipy.interpolate import CubicSpline
# 模拟有缺失的数据
t_full = np.arange(1, 366, 8)
ndvi_full = 0.3 + 0.4 * np.sin(2 * np.pi * t_full / 365)
# 随机删除20%的数据
np.random.seed(42)
missing_idx = np.random.choice(len(t_full), size=int(0.2*len(t_full)), replace=False)
t_obs = np.delete(t_full, missing_idx)
ndvi_obs = np.delete(ndvi_full, missing_idx)
# 三次样条插值
cs = CubicSpline(t_obs, ndvi_obs)
ndvi_interp = cs(t_full)
# 对比
plt.figure(figsize=(10, 5))
plt.plot(t_full, ndvi_full, 'g-', label='True', linewidth=2)
plt.plot(t_obs, ndvi_obs, 'ro', label='Observed')
plt.plot(t_full, ndvi_interp, 'b--', label='Interpolated', linewidth=2)
plt.legend()
plt.show()我曾经犯过的错：有一次做长时间序列分析，我直接用线性插值填补了连续10个缺失点。结果在植被生长旺季，插值出来的曲线完全偏离了实际物候规律——因为线性插值根本不知道植被在夏天应该快速生长。后来我改用时间序列分解插值，先提取出季节分量，再填补，效果就好多了。
3.4 三种方法的对比与选择
好了，三种方法都讲完了。你可能会问：到底该用哪个？
我的建议是：
- 数据噪声大、但缺失少 → 用S-G滤波。窗口调小一点，保留细节。
- 数据噪声大、缺失也多 → 先用插值填补缺失，再用Whittaker平滑。顺序别搞反了。
- 需要保留物候拐点 → 优先用Whittaker平滑，λ调小一点（0.5-2）。
- 数据量超大（比如几十万个像元） → 用S-G滤波，计算速度快。Whittaker需要解线性方程组，像元多了会慢。
核心原则：没有万能的方法。我每次拿到新数据，都会先画个原始曲线看看，再试两三种方法对比效果。你想想看，连数据长什么样都不知道，就直接套方法，那不是盲人摸象吗？
3.5 本章知识体系
下面这张图，是我自己画的知识结构图，帮你理清三种方法的关系：
这张图里，我把三种方法并列展示，下面列出了各自的关键参数。你可以把它当成一个“决策树”——根据你的数据特点，沿着分支找到最合适的方法。
好了，这一章就到这里。三种方法各有千秋，关键是多动手试。我当年也是反复调参、对比、踩坑，才慢慢摸出门道的。你只要把代码跑一遍，再换几组参数看看效果，很快就能找到感觉。