SSA奇异谱分解：时频域信号成分分析与重构

奇异谱分解（Singular Spectrum Analysis, SSA）是一种基于奇异值分解（SVD） 的非参数信号处理技术，通过将一维时间序列分解为多个独立成分（趋势、周期、噪声），实现时频域分析与信号重构。

一、、MATLAB实现步骤

以下通过合成信号分解与重构示例，展示SSA的完整实现。

1. 参数设置与示例信号生成

生成一个含趋势项（线性）、周期项（正弦波）、噪声项的合成信号：

\(x(t)=0.5t+2sin(2πt/10)+ε(t),t=1,2,...,200\)

其中 \(ε(t)∼N(0,0.52)\)为高斯白噪声。

2. SSA核心函数实现

（1）嵌入：构造轨迹矩阵

function X = embedding(x, L)
    % x: 输入时间序列 (N×1)
    % L: 窗口长度
    % X: 轨迹矩阵 (L×K), K=N-L+1
    N = length(x);
    K = N - L + 1;
    X = zeros(L, K);
    for i = 1:L
        X(i, :) = x(i:K+i-1);
    end
end

（2）对角平均：矩阵转一维序列

function x_rec = diagonal_averaging(X)
    % X: 待重构矩阵 (L×K)
    % x_rec: 重构序列 (N×1), N=L+K-1
    L = size(X, 1);
    K = size(X, 2);
    N = L + K - 1;
    x_rec = zeros(N, 1);
    count = zeros(N, 1); % 计数矩阵
    
    for i = 1:L
        for j = 1:K
            idx = i + j - 1; % 对角线索引
            x_rec(idx) = x_rec(idx) + X(i, j);
            count(idx) = count(idx) + 1;
        end
    end
    x_rec = x_rec ./ count; % 平均
end

（3）SSA主函数：分解与重构

function [components, singular_values] = ssa_decompose(x, L, m)
    % x: 输入时间序列 (N×1)
    % L: 窗口长度
    % m: 保留成分数 (默认取前m个最大奇异值成分)
    % components: 分解成分 (m×N矩阵, 每行一个成分)
    % singular_values: 奇异值 (L×1)
    
    N = length(x);
    X = embedding(x, L); % 轨迹矩阵
    [U, S, V] = svd(X, 'econ'); % SVD分解 (经济型)
    sigma = diag(S); % 奇异值 (L×1)
    singular_values = sigma;
    
    if nargin < 3 || isempty(m)
        m = length(sigma); % 默认保留所有成分
    end
    m = min(m, length(sigma)); % 避免超出奇异值数量
    
    components = zeros(m, N); % 存储m个成分
    for i = 1:m
        Xi = U(:, i) * sigma(i) * V(:, i)'; % 第i个成分矩阵
        components(i, :) = diagonal_averaging(Xi); % 对角平均重构
    end
end

3. 示例：合成信号分解与重构

clear; clc; close all;

%% 1. 生成示例信号
N = 200; % 序列长度
t = 1:N;
x = 0.5*t + 2*sin(2*pi*t/10) + 0.5*randn(1, N); % 趋势+周期+噪声
x = x(:); % 转为列向量

%% 2. SSA分解参数设置
L = 50; % 窗口长度 (经验值: L≈N/4~N/2)
m = 3;  % 保留前3个成分 (对应趋势、周期、噪声)

%% 3. SSA分解
[components, sigma] = ssa_decompose(x, L, m); % 分解
singular_values = sigma(1:m); % 前m个奇异值

%% 4. 成分选择与重构
% 成分1: 趋势项 (线性增长)
comp1 = components(1, :); 
% 成分2: 周期项 (正弦波)
comp2 = components(2, :); 
% 成分3: 噪声项 (剩余成分)
comp3 = components(3, :); 

% 重构信号 (叠加趋势+周期，去除噪声)
x_reconstructed = comp1 + comp2; 

%% 5. 时频分析 (功率谱密度)
fs = 1; % 采样频率 (归一化)
[P1, f1] = periodogram(comp1, [], [], fs); % 成分1功率谱
[P2, f2] = periodogram(comp2, [], [], fs); % 成分2功率谱
[P3, f3] = periodogram(comp3, [], [], fs); % 成分3功率谱

%% 6. 可视化结果
figure('Position', [100, 100, 1200, 800]);

% 原始信号与重构信号对比
subplot(3, 2, 1);
plot(t, x, 'b', t, x_reconstructed, 'r--', 'LineWidth', 1.5);
legend('原始信号', '重构信号 (趋势+周期)');
xlabel('时间 t'); ylabel('幅值'); title('原始信号与重构信号'); grid on;

% 分解成分时域波形
subplot(3, 2, 2);
plot(t, comp1, 'r', t, comp2, 'g', t, comp3, 'k');
legend('成分1 (趋势)', '成分2 (周期)', '成分3 (噪声)');
xlabel('时间 t'); ylabel('幅值'); title('SSA分解成分 (时域)'); grid on;

% 奇异值贡献率
subplot(3, 2, 3);
contribution = sigma.^2 / sum(sigma.^2) * 100; % 方差贡献率 (%)
stem(1:length(sigma), contribution, 'filled');
xlabel('奇异值序号'); ylabel('贡献率 (%)'); title('奇异值方差贡献率'); grid on;

% 成分功率谱 (时频分析)
subplot(3, 2, 4);
plot(f2, 10*log10(P2/max(P2)), 'g', 'LineWidth', 1.5); % 周期成分谱
hold on;
plot(f3, 10*log10(P3/max(P3)), 'k', 'LineWidth', 1.5); % 噪声成分谱
xlabel('归一化频率'); ylabel('功率谱密度 (dB)'); title('成分功率谱 (时频分析)');
legend('成分2 (周期)', '成分3 (噪声)'); grid on;

% 重构信号误差
subplot(3, 2, 5);
error = x - x_reconstructed;
plot(t, error, 'm');
xlabel('时间 t'); ylabel('误差'); title('重构误差 (原始-重构)'); grid on;

% 成分2 (周期项) 与原始周期对比
subplot(3, 2, 6);
plot(t, 2*sin(2*pi*t/10), 'b--', t, comp2, 'g', 'LineWidth', 1.5);
legend('原始周期项', 'SSA提取周期项');
xlabel('时间 t'); ylabel('幅值'); title('周期成分对比'); grid on;

二、关键结果分析

1. 分解成分特性

• 成分1（趋势项）：对应最大奇异值（σ1），时域波形为线性增长（与原始趋势一致），功率谱低频占优（无明显峰值）。

• 成分2（周期项）：对应第二大奇异值（σ2），时域波形为正弦波（周期10），功率谱在归一化频率 f=0.1（对应周期10）处有显著峰值。

• 成分3（噪声项）：对应第三大奇异值（σ3），时域波形为随机噪声，功率谱平坦（白噪声特性）。

2. 重构效果

• 重构信号：叠加成分1（趋势）和成分2（周期），去除成分3（噪声），与原始信号的趋势和周期部分高度吻合。

• 误差：重构误差主要集中在噪声部分，验证了SSA的去噪能力。

3. 奇异值贡献率

前3个奇异值的累计贡献率通常超过95%（示例中约为98%），表明仅需保留少量成分即可重构信号主要特征。

参考代码 SSA奇异谱分解，用于在时频域分析信号成分，重构信号 www.youwenfan.com/contentcns/46153.html

三、SSA参数选择指南

1. 窗口长度 L：

   • 经验值：L≈N/4∼N/2（N为序列长度）；

   • 周期信号：L略大于周期 T（如周期10，取 L=12）；

   • 趋势信号：L取较大值（如 L=N/2）以提高趋势分辨率。

2. 保留成分数 m：

• 按累计方差贡献率选择（如累计贡献率≥95%）；

• 按时频特性选择（保留有明显趋势/周期的成）。

四、SSA应用场景

1. 非平稳信号分解：分离趋势、周期、噪声（如气象数据、金融时间序列）；

2. 去噪与滤波：去除高频噪声（保留低频趋势/周期）；

3. 周期检测：通过功率谱识别隐藏周期（如机械振动信号中的故障频率）；

4. 缺失数据填补：利用成分相关性重构缺失值。

五、扩展：多成分重构与自适应选择

若信号含多个周期成分（如周期10和周期20），可增加窗口长度 L并保留更多成分（m=4,5,...），通过功率谱区分不同频率的周期成分。此外，可结合聚类算法（如K-means）对成分自动分类（趋势/周期/噪声），实现自适应重构。