奇异谱分析（SSA）MATLAB实现：遥感时序数据分解与重构
算法核心流程
奇异谱分析（SSA）是一种非参数化的时间序列分解方法，尤其适用于遥感影像的时空变化模式挖掘。其核心步骤包括：将一维序列嵌入到轨迹矩阵、进行奇异值分解（SVD）、基于成分分组重构序列。
%% 主函数：SSA分解与重构
function [reconstructed, components, singular_vals] = ssa_decompose(series, window_len, group_def)
% 输入参数说明：
% series : 输入的一维时间序列
% window_len : 窗口长度（嵌入维度），需为正整数且小于序列长度的一半
% group_def : 可选，分组定义，例如 {[1:2], [3:5], [6:end]} 表示三组
% 输出参数说明：
% reconstructed : 按照分组重构后的时间序列
% components : 每个奇异成分对应的时间序列（列向量形式）
% singular_vals : 奇异值向量
%% 参数校验
if nargin < 3
group_def = [];
end
N = length(series);
if window_len > N/2
error('窗口长度应小于序列长度的一半');
end
if floor(window_len) ~= window_len || window_len <= 0
error('窗口长度必须为正整数');
end
%% 第一步：构建轨迹矩阵（嵌入）
K = N - window_len + 1;
trajectory = zeros(window_len, K);
for idx = 1:K
trajectory(:, idx) = series(idx:idx+window_len-1);
end
%% 第二步：奇异值分解
% 使用经济型SVD分解
[U, S, V] = svd(trajectory, 'econ');
singular_vals = diag(S);
%% 第三步：对角平均获取各成分序列
num_comp = min(window_len, K);
components = zeros(N, num_comp);
for comp_idx = 1:num_comp
% 重构第comp_idx个成分矩阵
comp_matrix = S(comp_idx, comp_idx) * U(:, comp_idx) * V(:, comp_idx)';
% 对角平均操作
components(:, comp_idx) = diagonal_reconstruct(comp_matrix, N);
end
%% 第四步：分组重构
if isempty(group_def)
% 默认重构所有成分
reconstructed = sum(components, 2);
else
reconstructed = zeros(N, 1);
for g = 1:length(group_def)
group_indices = group_def{g};
group_series = sum(components(:, group_indices), 2);
reconstructed = reconstructed + group_series;
% 可视化分组结果
figure;
plot(series, 'k', 'LineWidth', 1.2); hold on;
plot(group_series, 'LineWidth', 1.5);
title(sprintf('第 %d 组重构 (索引: %s)', g, mat2str(group_indices)));
legend('原始序列', '本组重构');
xlabel('时间索引'); ylabel('值');
grid on;
end
end
%% 第五步：结果展示
plot_ssa_results(series, reconstructed, components, singular_vals, window_len);
end
%% 对角平均辅助函数
function ts = diagonal_reconstruct(matrix, len)
% 将轨迹矩阵沿对角线平均还原为时间序列
[L, K] = size(matrix);
ts = zeros(len, 1);
count = zeros(len, 1);
for i = 1:L
for j = 1:K
idx = i + j - 1;
if idx <= len
ts(idx) = ts(idx) + matrix(i, j);
count(idx) = count(idx) + 1;
end
end
end
count(count == 0) = 1;
ts = ts ./ count;
end
%% 可视化函数
function plot_ssa_results(original, recon, comps, svals, wlen)
figure('Name', 'SSA分析结果', 'Position', [100, 100, 1200, 800]);
% 子图1：原始与重构对比
subplot(3,2,1);
plot(original, 'b-', 'LineWidth', 1); hold on;
plot(recon, 'r--', 'LineWidth', 1.5);
title('原始时间序列 vs 重构结果');
xlabel('时间'); ylabel('值');
legend('原始', '重构');
grid on;
% 子图2：奇异值分布
subplot(3,2,2);
stem(svals, 'filled', 'LineWidth', 1.5);
title('奇异值分布');
xlabel('成分编号'); ylabel('奇异值');
grid on;
% 子图3：累积能量
subplot(3,2,3);
cum_energy = cumsum(svals.^2) / sum(svals.^2);
plot(cum_energy, 'g-o', 'LineWidth', 1.5);
title('累积能量百分比');
xlabel('成分数量'); ylabel('累积能量(%)');
grid on;
% 子图4：前几个主要成分
subplot(3,2,4);
plot(comps(:, 1:min(5, size(comps,2))), 'LineWidth', 1.2);
title('主要成分时间序列');
xlabel('时间'); ylabel('值');
legend(cellstr(num2str((1:min(5, size(comps,2)))', '成分%d')));
grid on;
% 子图5-6：功率谱密度
subplot(3,2,[5,6]);
periodogram(recon, [], [], 1024);
title('重构序列功率谱');
grid on;
saveas(gcf, 'SSA_results_visualization.png');
end
%% 批量处理遥感影像的函数
function batch_process_remote_sensing(image_cube, wlen, group_definition)
% image_cube: 三维数据 (行数 × 列数 × 时间步)
% wlen: 窗口长度
% group_definition: 分组定义
[rows, cols, steps] = size(image_cube);
output_cube = zeros(rows, cols, steps);
h = waitbar(0, '遥感影像SSA处理中...');
for r = 1:rows
for c = 1:cols
pixel_series = squeeze(image_cube(r, c, :))';
[recon, ~, ~] = ssa_decompose(pixel_series, wlen, group_definition);
output_cube(r, c, :) = recon;
end
waitbar(r/rows, h, sprintf('行进度: %d/%d', r, rows));
end
close(h);
figure;
montage(reshape(output_cube, [rows, cols, 1, steps]), 'DisplayRange', []);
title('SSA重构后的影像序列');
save('ssa_processed_data.mat', 'output_cube');
end
%% 示例演示函数
function demo_ssa_analysis()
% 构造模拟数据：趋势 + 季节 + 噪声 + 异常
t = 1:120;
trend = 0.03 * t;
seasonal = 8 * sin(2*pi*t/12);
noise = 1.5 * randn(1, 120);
signal = trend + seasonal + noise;
signal(40) = signal(40) + 15; % 模拟异常
signal(80) = signal(80) - 12;
% 选择窗口长度和分组
W = 30;
groups = {1:2, 3:6, 7:12, 13:end};
[rec, comps, svals] = ssa_decompose(signal, W, groups);
% 输出统计
rmse_val = sqrt(mean((signal - rec).^2));
fprintf('重构均方根误差 (RMSE): %.4f\n', rmse_val);
fprintf('最大奇异值: %.4f\n', max(svals));
fprintf('前三成分能量占比: %.2f%%\n', sum(svals(1:3).^2)/sum(svals.^2)*100);
end
关键功能模块
- 嵌入操作：将一维序列映射为L×K的轨迹矩阵，K = N-L+1。
- SVD分解：利用奇异值分解分离出主要成分。
- 对角平均：将每个重构矩阵还原为时间序列。
- 分组重构：允许用户自定义分组，例如分离趋势、周期和噪声。
- 批量影像处理：逐像元应用SSA，适合多波段遥感数据。
分组配置策略参考
| 成分类型 | 典型索引范围 | 物理意义 |
|---|---|---|
| 趋势项 | 1-2 | 长期缓慢变化 |
| 季节项 | 3-4 | 年/季周期波动 |
| 短期事件 | 5-6 | 异常或突变信号 |
| 噪声 | 7:end | 残余随机噪声 |
遥感应用案例
1. NDVI数据去噪与趋势提取
% 假设已有NDVI 3D数据 ndvi_cube
W = 36; % 3年窗口
group_def = {1:2, 3:8, 9:12, 13:end};
[recon, comps, ~] = ssa_decompose(ndvi_cube(:), W, group_def);
seasonal_part = comps(:, group_def{2});
deseasoned_ndvi = ndvi_cube - reshape(seasonal_part, size(ndvi_cube));
% 趋势可视化
trend_part = comps(:, group_def{1});
figure; plot(trend_part); title('NDVI长期趋势');
2. 地表温度异常检测
% 温度数据 temp_cube (h×w×t)
W = 60; % 2个月
group_def = {1, 2:4, 5:10, 11:end};
[~, comps, ~] = ssa_decompose(temp_cube(:), W, group_def);
residual = comps(:, group_def{end});
anomaly_map = reshape(residual, size(temp_cube));
% 二值化异常（2倍标准差）
threshold = 2 * std(residual);
binary_anomaly = anomaly_map > threshold;
3. 降水周期模式分析
% 降水序列
W = 365; % 一年
group_def = {1, 2:3, 4:12, 13:end};
[~, comps, ~] = ssa_decompose(precip_series, W, group_def);
cyclical_comp = sum(comps(:, group_def{2}), 2);
% 重构周期图
figure; contourf(reshape(cyclical_comp, [h,w])); colorbar;
性能优化建议
- 分块处理：对于超大影像，按块处理减少内存压力：
block_h = 256; block_w = 256;
for r = 1:block_h:rows
for c = 1:block_w:cols
block = image_cube(r:min(r+block_h-1,rows), c:min(c+block_w-1,cols), :);
% 处理此块...
end
end
- GPU加速：将数据移至GPU运算：
X_gpu = gpuArray(X);
[U_g, S_g, V_g] = svd(X_gpu, 'econ');
U = gather(U_g); S = gather(S_g); V = gather(V_g);
- 并行循环：使用parfor处理每个像元：
parfor r = 1:rows
for c = 1:cols
% 并行处理
end
end
常见问题与解决
窗口长度选择
推荐取序列长度的1/3左右，或基于主周期估算：
[pxx, f] = periodogram(series);
[~, max_idx] = max(pxx);
main_period = 1/f(max_idx);
L = round(main_period * 1.5);
成分选取
使用累积能量准则：选取前k个成分使得累积能量超过90%：
cum_energy = cumsum(svals.^2) / sum(svals.^2);
k = find(cum_energy >= 0.9, 1, 'first');
端点效应处理
通过镜像延拓减少边界失真：
extended = [flipud(series(1:L)); series; flipud(series(end-L+1:end))];
[recon, ~, ~] = ssa_decompose(extended, L);
recon = recon(L+1:end-L);
扩展功能
多变量SSA
function multi_ssa_analysis(data_matrix, L)
% data_matrix: 变量×时间
[U, S, V] = svd(data_matrix', 'econ');
% 后续处理与单变量类似
end
SSA-ARIMA混合模型
denoised_series = ssa_decompose(ts, L);
model = arima(2,1,2);
fit = estimate(model, denoised_series);
pred = forecast(fit, 12);
SSA-LSTM预测
% 使用前几个成分为特征输入LSTM
features = comps(:, 1:5);
% 构建LSTM数据集并训练...
总结
本SSA实现提供了从嵌入、分解到重构的完整流程，特别针对遥感影像时间序列优化，支持灵活分组、批量处理、可视化输出。通过合理选择窗口和分组，可以有效分离趋势、周期与噪声，适用于NDVI、地表温度、降水等多类遥感数据。
