使用Python 实现奇异谱分析SSA（Singular Spectrum Analysis）算法

简介

SSA是一种近几年新兴的时序分析算法，通过对时间序列进行升维后做SVD，进而分组重构来实现对时间序列的成分分析，目前用于预测、插值补全以及去噪
具体步骤包括四个部分：
1.构建轨迹矩阵(Embedding)
2.轨迹矩阵SVD分解
3.分组
4.对角平均

1.准备工作

这里我们生成一组间隔带有噪声的序列,包含了周期项、趋势项和随机噪声

import numpy as np
import matplotlib.pyplot as plt
number = 30
days = 180
tend_sequence = np.linspace(2,-2,num=days)
time = np.linspace(0,2*np.pi*days/number,num=days)
sin_sequence = np.sin(time)
signal = sin_sequence+tend_sequence
noise = np.random.randn(days)
sequence = noise+sin_sequence+tend_sequence
plt.plot(tend_sequence,color='yellow',label='Tendency')
plt.plot(sin_sequence,color='blue',label='Periodicity')
plt.plot(noise,color='red',label='noise')
plt.legend()
plt.show()
plt.plot(sequence,color='red',label='with noise')
plt.plot(signal,color='blue',label='original')
plt.legend()
plt.show()

2.SSA具体实现

在实际操作中，我们将从以下步骤进行计算：
1. 对于长度为 $N$ 的时间序列，以 $L$ 为滑动窗口，依次进行采样，获得大小为 $L\times K$ 的轨迹矩阵 $X$，其中 $K=N-L+1$
2. 对 $X$ 做SVD分解，文献中经常给出完整的计算方法（先计算$X{X}^T$等），但是经过检验直接使用numpy.linalg.svd()是可行的。
分解后获得 $U、\Sigma 、V^T$
3. 构建$X$的基本矩阵（elementary matrices），这里$X_i={\sigma }_i{U}_i{V}_i^T$，这里认为轨迹矩阵$X$是由$L$个这样的基本矩阵加和而成的。
4. 将基本矩阵$X_i$通过对角平均的操作，还原到长度为$N$的时间序列，即为原序列的一个构成成分

这一部分代码参考了

[机器学习] 奇异谱分析(SSA)原理及Python实现

windowlen = 45
serieslen = days
#1.嵌入
K = serieslen-windowlen+1
X = np.zeros((windowlen,K))
for i in range(windowlen):
    X[i,:] = sequence[i:i+K]
#2.SVD分解
U,sigma,VT = np.linalg.svd(X,full_matrices=False)
ZT = np.zeros((VT.shape))
for n in range(len(sigma)):
    ZT[n,:] = sigma[n]*VT[n,:]

#3.对每个Ui*ZTi所得到的Xi进行对角平均获得重构子序列，一般我们假设L<K
RC = np.zeros((windowlen,serieslen))
L = windowlen
for num in range(U.shape[1]):
    Xi = np.outer(U[:,num],ZT[num,:])
    for k in range(serieslen):
        sum=0
        if(k<=L-1):
            for p in range(k+1):
                sum+=Xi[p,k-p]
            RC[num,k] = sum/(k+1)
        elif(k<=K-1):
            for p in range(L):
                sum+=Xi[p,k-p]
            RC[num,k] = sum/L
        elif(k<=serieslen-1):
            for p in range(k-K+1,serieslen-K+1):
                sum+=Xi[p,k-p]
            RC[num,k] = sum/(serieslen-k)

3.结果分析

现在我们可以看到，所获得的RC是由$L$个成分序列构成的矩阵，每一行对应原序列的一个成分

showN = 8
for m in range(showN):
    plt.subplot(showN,1,m+1)
    plt.plot(RC[m,:])
    plt.ylabel(('RC {}'.format(m)))
plt.show()

查看周期项分离效果

plt.plot(RC[1,:]+RC[2,:],color='blue',label = 'predicted')
plt.plot(sin_sequence,color='red',label = 'periodicity')
plt.legend()
plt.show()

查看趋势项分离效果

plt.plot(RC[0,:],color='blue',label = 'predicted')
plt.plot(tend_sequence,color='red',label = 'trend')
plt.legend()
plt.show()

查看整体去噪效果

选用RC0\1\2\3四项进行重组

plt.plot(RC[0,:]+RC[1,:]+RC[2,:]+RC[3,:],color='blue',label = 'predicted')
plt.plot(signal,color='red',label = 'original')
plt.legend()
plt.show()

目前所掌握的知识还不足以将SSA发展为一种通用的方法，原因是窗口长度 $L$ 的长度确定具有较强的主观性，重现原序列时对各个成分的取舍也没有具体的标准，还需要进一步的学习。

博主不是信号处理及数学方面的专业人士，科研项目中用到了这一块，简单记录下学习过程，如果有不正确的地方麻烦指出，欢迎交流~

参考文献：
[1] GOLYANDINA N, ZHIGLJAVSKY A. Singular Spectrum Analysis for Time Series[M/OL]. 2013. DOI:10.1007/978-3-642-34913-3.