Python SciPy时间序列分析：滤波去趋势与周期检测实战

脚本专家

时间序列分析处理的是按时间排列的数据：判断趋势、周期性、噪声强弱，或者两个序列之间的延迟关系。SciPy 在信号处理和统计检验方面提供了不少现成工具，配合 NumPy 与 Matplotlib，可以覆盖日常分析的大部分需求。下面直接用一组模拟的全年每日温度数据，依次演示移动平均、低通滤波、去趋势、频谱分析、自相关和互相关的代码实现与参数含义。

1. 构造带季节性和趋势的温度数据

原文用 NumPy 模拟 365 天温度，由三部分组成：周期为 365 天的正弦波表示季节性，0.05 * days 表示缓慢线性上升趋势，再叠加正态随机噪声。代码先固定随机种子，保证结果可复现。

2. import numpy as np
3. import matplotlib.pyplot as plt
5. plt.rcParams['font.sans-serif'] = ['SimHei']
6. plt.rcParams['axes.unicode_minus'] = False
7. np.random.seed(0)
9. days = np.arange(365)
10. seasonal = 10 * np.sin(2 * np.pi * days / 365)
11. trend = 0.05 * days
12. noise = np.random.normal(0, 2, 365)
13. temperature_data = 20 + seasonal + trend + noise
15. plt.plot(days, temperature_data)
16. plt.xlabel('天数')
17. plt.ylabel('温度')
18. plt.title('模拟的每日温度数据')
19. plt.show()

复制代码

绘制后是一条上下波动且整体略微上升的曲线，噪声让线条带有毛刺。这组 temperature_data 就是后续所有分析的输入。

2. 用移动平均压制高频噪声

原始数据中的高频噪声会干扰趋势判断。移动平均的思路是取一个窗口，例如 7 天，计算窗口内均值，再向后滑动。NumPy 的 convolve 可以简洁实现。

2. def moving_average(data, window_size):
3.     return np.convolve(data, np.ones(window_size) / window_size, mode='valid')
5. smoothed = moving_average(temperature_data, 7)
6. plt.plot(days, temperature_data, label='原始数据', alpha=0.5)
7. plt.plot(days[6:], smoothed, label='7天移动平均', linewidth=2)
8. plt.legend()
9. plt.show()

复制代码

窗口取 7 天，是因为温度数据通常以周为周期存在波动，7 天平均能抹掉大部分随机噪声，同时保留季节性的轮廓。平滑后的曲线更干净，上升趋势和正弦波动都更清晰。需要留意，mode='valid' 会使输出长度变为 len(data)-window_size+1，因此画图时横轴使用了 days[6:]。

3. Butterworth 低通滤波：更灵活地控制截止频率

移动平均本质上是一个简单低通滤波器，但如果想更精确地控制截止频率，可以改用 SciPy 的 signal.butter 和 filtfilt。Butterworth 滤波器会衰减高于指定频率的成分，filtfilt 做双向滤波，避免相位偏移。

2. from scipy.signal import butter, filtfilt
4. def lowpass_filter(data, cutoff, fs, order=4):
5.     nyquist = 0.5 * fs
6.     normal_cutoff = cutoff / nyquist
7.     b, a = butter(order, normal_cutoff, btype='low', analog=False)
8.     return filtfilt(b, a, data)
10. filtered = lowpass_filter(temperature_data, cutoff=0.1, fs=1)
11. plt.plot(days, temperature_data, alpha=0.5, label='原始数据')
12. plt.plot(days, filtered, label='低通滤波后', linewidth=2)
13. plt.legend()
14. plt.show()

复制代码

这里采样频率 fs=1，表示每天一个采样点；截止频率 cutoff=0.1，表示保留周期大于 10 天的成分，把更快的波动滤掉。结果和移动平均类似，但边缘处理更平滑，没有移动平均两端缺数据的问题。

4. 去趋势：分离长期上升与季节性波动

温度数据里包含线性上升趋势。如果只想观察季节性的正弦波动，可以用 scipy.signal.detrend 把线性趋势减掉。

2. from scipy.signal import detrend
4. detrended = detrend(temperature_data)
5. plt.plot(days, detrended)
6. plt.xlabel('天数')
7. plt.ylabel('去趋势后的温度')
8. plt.title('去趋势结果')
9. plt.show()

复制代码

去趋势后的曲线围绕 0 上下波动，正弦形状更明显。要注意，detrend 默认只去线性趋势；如果原始数据中的趋势不是线性的，效果可能不够。这里的数据恰好是线性趋势，所以结果比较干净。

从分析流程看，做频谱分析、自相关或寻找周期性之前，通常先去趋势。否则线性趋势会在低频段产生很强的能量，掩盖真正的周期信号；去趋势也能让数据均值趋于稳定，满足平稳性要求。换句话说，detrend 没有改变数据的形状，只是把数据整体平移并拉平，让它围绕 0 波动，相当于减掉了 20 + 0.05 * days 这条背景直线。因此去趋势结果图的 Y 轴围绕 0 左右波动，而初始模拟数据图围绕 30 左右波动。

5. 频谱分析：用 periodogram 找主要周期

想知道数据里有哪些周期成分，频谱分析很直接。scipy.signal.periodogram 可以计算功率谱密度，横轴是频率，纵轴是功率，频率的倒数就是周期。

2. from scipy.signal import periodogram
4. frequencies, power = periodogram(temperature_data, fs=1)
5. plt.semilogy(frequencies, power)
6. plt.xlabel('频率（周期/天）')
7. plt.ylabel('功率')
8. plt.title('周期图')
9. plt.show()

复制代码

图上会在频率 1/365 ≈ 0.00274 附近出现明显峰值，对应 365 天的年周期。使用对数纵轴，是因为不同频率处的功率差异可能很大。如果还看到其他小峰值，可能是噪声或谐波，但主峰位置直接告诉你数据里最强的周期。

6. 自相关：观察当前值与过去值的关系

自相关衡量时间序列和它自己延迟若干天后的相关性。如果数据有季节性，自相关会在滞后 365 天附近出现高点。下面用 scipy.stats.pearsonr 计算滞后 1 到 30 天的相关系数。

2. from scipy.stats import pearsonr
4. lags = range(1, 31)
5. autocorrs = [pearsonr(temperature_data[:-lag], temperature_data[lag:])[0] for lag in lags]
6. plt.stem(lags, autocorrs)
7. plt.xlabel('滞后（天）')
8. plt.ylabel('自相关')
9. plt.title('30天内的自相关')
10. plt.show()

复制代码

由于数据包含正弦季节性，自相关会呈现波浪形，滞后越接近 365 的整数倍，相关性越高。这里只画了 30 天，能看到相关性逐渐下降然后可能回升，这反映了短期的周期结构。

7. 互相关：分析两个序列之间的延迟关系

有时关心两个序列是否相关，以及一个序列是否领先另一个序列。互相关就是处理这类问题的工具。原文生成了一组模拟湿度数据，让它与温度存在 5 天延迟关系，然后计算不同偏移下的相关系数。

2. humidity_data = 70 + 0.5 * temperature_data[:-5] + np.random.normal(0, 1, 360)
3. shifts = range(-10, 11)
4. cross_corrs = [pearsonr(temperature_data[:360], np.roll(humidity_data, shift))[0] for shift in shifts]
5. plt.stem(shifts, cross_corrs)
6. plt.xlabel('偏移量')
7. plt.ylabel('互相关')
8. plt.title('温度与湿度的互相关')
9. plt.show()

复制代码

结果会在偏移 +5 或 -5 附近出现峰值，具体方向取决于 np.roll 的定义。这告诉我们湿度对温度有大约 5 天的滞后响应。在实际分析中，互相关能帮助判断因果方向，或者找到两个序列的最佳对齐方式。

8. 小结与参数选择建议

整个流程可以概括为：构造数据 → 移动平均 → 低通滤波 → 去趋势 → 频谱分析 → 自相关 → 互相关。移动平均和低通滤波适合平滑与去噪；去趋势用于分离长期趋势和周期成分；频谱和自相关用于揭示周期；互相关用于两个序列的延迟分析。这些技术也是更高级建模，例如 ARIMA、傅里叶拟合的基础。实际使用时，可以在自己的数据上替换 temperature_data，并调整窗口大小、截止频率、滤波器阶数、滞后范围和偏移范围，观察结果如何变化。