奇异谱分析


奇异谱分析简介


奇异谱分析(SSA)是一种谱估计方法，可以作为一种自适应信号分解工具，将信号分解为若干个分量，在不剔除任何信息的情况下，所有SSA分量之和可以完美重构回之前的时间序列。

但奇异谱分析的亮点在于，它可以将时间序列分解成趋势项、季节项、噪音等等，将其中的噪音进行提取并剔除是最常用的一种应用，因此奇异谱分析在时间序列预测中多用于降噪处理。


奇异谱分析算法过程


总述

奇异谱分析的原理大致可以分为四步：构成轨迹矩阵——奇异值分解——分组——通过对角线平均法重构信号。接下来对每一步进行通俗的介绍，以便理解。

嵌入

将一维原始时间序列数据构建轨迹矩阵

输入：1.原始时间序列  2.窗口长度

输出：轨迹矩阵

将原始序列(一维数据)类似样本熵分组那样，依次往后平移一个进行区间划分，将划分出来的数据作为一列，组合拼成轨迹矩阵。

注意：分组的个数一定要保证窗口长度大于1小于等于N/2。

# step1 嵌入windowLen = 20  # 嵌入窗口长度    一定确保窗口长度小于seriesLen = len(series)  # 序列长度K = seriesLen - windowLen + 1  #轨迹矩阵一共有多少列数据X = np.zeros((windowLen, K))for i in range(K):    X[:, i] = series[i:i + windowLen]
# step1 嵌入
windowLen = 20  # 嵌入窗口长度    一定确保窗口长度小于
seriesLen = len(series)  # 序列长度
K = seriesLen - windowLen + 1  #轨迹矩阵一共有多少列数据
X = np.zeros((windowLen, K))

for i in range(K):
    X[:, i] = series[i:i + windowLen]




分解

进行奇异值及其对应特征向量的分解

输入：轨迹矩阵  

输出：多个可以加和生成轨迹矩阵的分量

通过奇异值分解，求奇异值以及对应的特征向量，其中的个特征值由大到小对角线排列。

通过把特征向量矩阵单拎出来，用单列特征项向量和轨迹矩阵进行处理，生成轨迹矩阵的分量(我的理解：每个分量分别代表每个特征值对应的在轨迹矩阵中的信息)

# step2: SVD奇异值分解， U和sigma已经按升序排序U, sigma, VT = np.linalg.svd(X, full_matrices=False)for i in range(VT.shape[0]):    VT[i, :] *= sigma[i]A = VT
# step2: SVD奇异值分解， U和sigma已经按升序排序
U, sigma, VT = np.linalg.svd(X, full_matrices=False)

for i in range(VT.shape[0]):
    VT[i, :] *= sigma[i]
A = VT




分组

将分解得到的矩阵单拎出来，进行分组。去噪的环节。

输入：1.分组数量 2.如何分组 3.其中的哪几组用于重构信号  

输出：保留了哪几组轨迹矩阵的分量的加和

可以不等距的将轨迹矩阵的分量进行分组，然后决定哪几组是用来重构的。

最常见的方法是，将轨迹矩阵的分量划分为两组，第一组有多少个是根据提前设置的阈值，当特征值的依次加和/总特征值的和 达到阈值时，以此为界限，进行两组的划分。

# 分组的代码部分和重构混合在一起，见重构部分代码。

重构

对角平均法重构信号，将矩阵转换成多个一维数据，然后进行相加成一条一维数据。

输入：保留了哪几组轨迹矩阵的分量的加和  

输出：重构之后的一维数据(已去噪)

使用对角平均法，将进行重构的每组轨迹矩阵分量的加和分别进行对角平均，每个对应一个一维的数据，将多个一维的数据进行加和，完成重构。

# 分组+重组rec = np.zeros((windowLen, seriesLen))for i in range(windowLen):    for j in range(windowLen - 1):        for m in range(j + 1):            rec[i, j] += A[i, j - m] * U[m, i]        rec[i, j] /= (j + 1)    for j in range(windowLen - 1, seriesLen - windowLen + 1):        for m in range(windowLen):            rec[i, j] += A[i, j - m] * U[m, i]        rec[i, j] /= windowLen    for j in range(seriesLen - windowLen + 1, seriesLen):        for m in range(j - seriesLen + windowLen, windowLen):            rec[i, j] += A[i, j - m] * U[m, i]        rec[i, j] /= (seriesLen - j)rrr = np.sum(rec, axis=0)  # 选择重构的部分，这里选了全部# 根据rec的索引选择需要的部分就可以了，比如rrr = np.sum(rec[1:5, :], axis=0)
# 分组+重组
rec = np.zeros((windowLen, seriesLen))
for i in range(windowLen):
    for j in range(windowLen - 1):
        for m in range(j + 1):
            rec[i, j] += A[i, j - m] * U[m, i]
        rec[i, j] /= (j + 1)
    for j in range(windowLen - 1, seriesLen - windowLen + 1):
        for m in range(windowLen):
            rec[i, j] += A[i, j - m] * U[m, i]
        rec[i, j] /= windowLen
    for j in range(seriesLen - windowLen + 1, seriesLen):
        for m in range(j - seriesLen + windowLen, windowLen):
            rec[i, j] += A[i, j - m] * U[m, i]
        rec[i, j] /= (seriesLen - j)

rrr = np.sum(rec, axis=0)  # 选择重构的部分，这里选了全部
# 根据rec的索引选择需要的部分就可以了，比如rrr = np.sum(rec[1:5, :], axis=0)




至此，奇异谱的分析的降噪工作已完成。

注：分组环节如果不剔除噪音的话，可以完美重构回之前的原始序列。