当前位置：首页 > news >正文

保姆级教程：用Python的NumPy和Matplotlib一步步拆解时间序列（含SSA算法完整代码）

news 2026/6/2 8:14:25

从零实现时间序列奇异谱分析：用Python完整复现SSA算法

时间序列分析是数据科学中的核心技能之一，而奇异谱分析（Singular Spectrum Analysis, SSA）作为一种无需预设模型的方法，近年来在金融预测、气候研究、信号处理等领域展现出独特优势。本文将带您从零开始，用Python完整实现SSA算法，并通过可视化手段直观理解其工作原理。

1. 环境准备与数据生成

在开始SSA算法实现前，我们需要搭建Python环境并生成一组包含趋势、周期和噪声的合成数据。这不仅能验证算法效果，也为后续分析提供明确参照。

import numpy as np import matplotlib.pyplot as plt plt.style.use('seaborn') # 使用更美观的绘图样式 # 参数设置 days = 180 # 总天数 period = 30 # 周期天数 np.random.seed(42) # 固定随机种子保证可复现性 # 生成趋势项（线性下降） trend = np.linspace(2, -2, num=days) # 生成周期项（正弦波） time_points = np.linspace(0, 2*np.pi*days/period, num=days) seasonality = np.sin(time_points) * 1.5 # 生成随机噪声 noise = np.random.normal(0, 0.5, days) # 组合信号 signal = trend + seasonality noisy_signal = signal + noise # 可视化各成分 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8)) ax1.plot(trend, label='Trend', color='orange') ax1.plot(seasonality, label='Seasonality', color='green') ax1.plot(noise, label='Noise', alpha=0.5, color='red') ax1.set_title('Signal Components') ax1.legend() ax2.plot(noisy_signal, label='Noisy Signal', color='blue') ax2.plot(signal, label='True Signal', linestyle='--', color='black') ax2.set_title('Composite Signal') ax2.legend() plt.tight_layout() plt.show()

这段代码生成了三个关键组成部分：

趋势项：表现为线性下降的黄色线条
周期项：呈现规则波动的绿色正弦曲线
噪声项：随机分布的红色点状波动

可视化输出将清晰展示原始信号与加噪信号的对比，为后续去噪效果评估提供基准。

2. SSA算法核心四步实现

SSA算法的核心流程可分为四个关键步骤，我们将逐一实现并解释每个步骤的数学意义和实现细节。

2.1 构建轨迹矩阵（嵌入阶段）

轨迹矩阵构建是SSA的第一步，其本质是将一维时间序列升维到矩阵空间。关键参数是窗口长度L的选择，它直接影响分析结果。

def build_trajectory_matrix(series, L): """构建轨迹矩阵""" N = len(series) K = N - L + 1 X = np.zeros((L, K)) for i in range(L): X[i, :] = series[i:i+K] return X # 窗口长度选择（通常取序列长度的1/3到1/2） L = 60 traj_matrix = build_trajectory_matrix(noisy_signal, L) # 可视化轨迹矩阵 plt.figure(figsize=(10, 6)) plt.imshow(traj_matrix, cmap='viridis', aspect='auto') plt.colorbar(label='Value') plt.title(f'Trajectory Matrix (L={L})') plt.xlabel('Window Position') plt.ylabel('Window Length') plt.show()

窗口长度L的选择至关重要：

L值过小：难以捕捉长周期模式
L值过大：可能导致过度拟合和计算效率下降
经验法则：通常取N/3到N/2之间，其中N为序列长度

2.2 SVD分解与特征分析

对轨迹矩阵进行奇异值分解（SVD），这是SSA算法的核心数学操作，能够揭示时间序列的内在结构。

def perform_svd(X): """执行SVD分解并计算贡献率""" U, sigma, VT = np.linalg.svd(X, full_matrices=False) # 计算各分量的能量贡献 energy = sigma**2 / np.sum(sigma**2) return U, sigma, VT, energy U, sigma, VT, energy = perform_svd(traj_matrix) # 绘制奇异值谱 plt.figure(figsize=(12, 5)) plt.subplot(121) plt.plot(sigma, 'o-') plt.title('Singular Values (Log Scale)') plt.yscale('log') plt.xlabel('Component Index') plt.ylabel('Singular Value') plt.subplot(122) plt.plot(np.cumsum(energy), 'o-') plt.title('Cumulative Energy') plt.xlabel('Component Index') plt.ylabel('Explained Variance Ratio') plt.axhline(0.9, color='red', linestyle='--', alpha=0.5) plt.tight_layout() plt.show()

SVD分解后我们得到三个矩阵：

U矩阵：代表时间域的特征模式
Σ矩阵：对角矩阵，包含奇异值（反映各成分重要性）
V矩阵：代表序列域的特征模式

奇异值谱图可帮助我们确定重要成分的数量，通常前几个成分对应趋势和主要周期。

2.3 成分分组与重构

根据奇异值谱的分析结果，我们需要对成分进行合理分组，这是SSA最具技巧性的环节。

def reconstruct_component(U, sigma, VT, component_idx, L): """重构指定成分""" s = np.zeros(len(sigma)) s[component_idx] = sigma[component_idx] S = np.diag(s) X_comp = U @ S @ VT return diagonal_averaging(X_comp, L) def diagonal_averaging(X, L): """对角平均将矩阵转换回时间序列""" N = X.shape[1] + L - 1 reconstructed = np.zeros(N) for k in range(N): if k < L - 1: reconstructed[k] = np.mean([X[i, k-i] for i in range(k+1)]) elif k < X.shape[1]: reconstructed[k] = np.mean([X[i, k-i] for i in range(L)]) else: reconstructed[k] = np.mean([X[i, k-i] for i in range(k-X.shape[1]+1, L)]) return reconstructed # 重构前4个主要成分 components = [] for i in range(4): components.append(reconstruct_component(U, sigma, VT, i, L)) # 可视化各成分 plt.figure(figsize=(12, 8)) for i, comp in enumerate(components): plt.subplot(4, 1, i+1) plt.plot(comp) plt.title(f'Component {i+1} (Energy: {energy[i]:.2%})') plt.tight_layout() plt.show()

分组策略通常包括：

趋势成分：通常对应第一个奇异值
周期成分：对应后续几个显著奇异值
噪声成分：剩余的小奇异值

2.4 结果分析与可视化

最后一步是将重构的成分与原始信号进行对比，评估SSA的分解效果。

# 组合趋势和周期成分 reconstructed_trend = components[0] reconstructed_seasonality = components[1] + components[2] # 可视化对比 plt.figure(figsize=(12, 8)) plt.subplot(311) plt.plot(trend, label='True Trend', color='black', linestyle='--') plt.plot(reconstructed_trend, label='Reconstructed Trend') plt.title('Trend Component Comparison') plt.legend() plt.subplot(312) plt.plot(seasonality, label='True Seasonality', color='black', linestyle='--') plt.plot(reconstructed_seasonality, label='Reconstructed Seasonality') plt.title('Seasonal Component Comparison') plt.legend() plt.subplot(313) plt.plot(signal, label='True Signal', color='black', linestyle='--') plt.plot(reconstructed_trend + reconstructed_seasonality, label='Reconstructed Signal') plt.title('Full Signal Comparison') plt.legend() plt.tight_layout() plt.show()

通过对比图可以直观评估：

趋势成分的捕捉是否准确
周期信号的振幅和相位是否匹配
噪声抑制的效果如何

3. 参数优化与实用技巧

SSA算法的效果很大程度上依赖于参数选择和操作技巧，本节将分享实际应用中的经验。

3.1 窗口长度L的优化选择

窗口长度L是SSA最重要的参数，没有统一选择标准，但有一些实用准则：

选择方法	说明	适用场景
1/3规则	取N/3，N为序列长度	通用场景
周期倍数	取主要周期的整数倍	已知周期特征
试错法	尝试多个L值比较结果	复杂序列

# 测试不同L值的效果 L_values = [30, 60, 90, 120] fig, axes = plt.subplots(len(L_values), 1, figsize=(12, 10)) for ax, L in zip(axes, L_values): traj_matrix = build_trajectory_matrix(noisy_signal, L) U, sigma, VT, _ = perform_svd(traj_matrix) comp1 = reconstruct_component(U, sigma, VT, 0, L) ax.plot(trend, label='True Trend', linestyle='--') ax.plot(comp1, label=f'L={L}') ax.legend() ax.set_title(f'Window Length L = {L}') plt.tight_layout() plt.show()

实验表明，L=60（序列长度的1/3）在本案例中表现最佳，既能捕捉趋势又不会引入过多噪声。

3.2 成分选择与噪声处理

如何选择有效成分、舍弃噪声成分是SSA应用的关键决策点。常用方法包括：

能量阈值法：保留累计能量达到90%的成分
奇异值拐点法：选择奇异值下降曲线的拐点
频率分析法：通过FFT检查各成分的频率特征

# 成分能量分析示例 cum_energy = np.cumsum(energy) threshold = 0.9 # 90%能量阈值 n_components = np.argmax(cum_energy >= threshold) + 1 print(f'保留前{n_components}个成分可解释{threshold:.0%}的方差')

在实际项目中，建议结合多种方法综合判断，并通过业务知识验证成分的合理性。

4. 完整SSA算法封装与扩展

为了便于实际应用，我们将前述步骤封装为完整的Python类，并添加实用功能。

class SingularSpectrumAnalysis: def __init__(self, series, window_length=None): self.series = series self.N = len(series) self.L = window_length if window_length else self.N // 3 self.K = self.N - self.L + 1 def decompose(self, n_components=None): """执行完整SSA分解""" # 构建轨迹矩阵 self.X = build_trajectory_matrix(self.series, self.L) # SVD分解 self.U, self.sigma, self.VT, self.energy = perform_svd(self.X) # 确定成分数量 if n_components is None: cum_energy = np.cumsum(self.energy) n_components = np.argmax(cum_energy >= 0.9) + 1 # 重构成分 self.components = [] for i in range(n_components): comp = reconstruct_component(self.U, self.sigma, self.VT, i, self.L) self.components.append(comp) return self.components def reconstruct(self, component_indices): """重构指定成分组合""" return np.sum([self.components[i] for i in component_indices], axis=0) def plot_components(self, n_components=6): """可视化主要成分""" plt.figure(figsize=(12, 8)) for i in range(min(n_components, len(self.components))): plt.subplot(n_components, 1, i+1) plt.plot(self.components[i]) plt.title(f'Component {i+1} (Energy: {self.energy[i]:.2%})') plt.tight_layout() plt.show() # 使用示例 ssa = SingularSpectrumAnalysis(noisy_signal, L=60) components = ssa.decompose() ssa.plot_components() # 重构信号 reconstructed = ssa.reconstruct([0, 1, 2]) # 组合趋势和前两个周期成分 # 性能评估 from sklearn.metrics import mean_squared_error mse = mean_squared_error(signal, reconstructed) print(f'Reconstruction MSE: {mse:.4f}')

这个封装类提供了以下扩展功能：