[机器学习] 奇异谱分析(SSA)原理及Python实现 |
您所在的位置:网站首页 › ssa中文翻译 › [机器学习] 奇异谱分析(SSA)原理及Python实现 |
最近做时间序列分析的时候需要用到奇异谱分析,发现网上可以查到的资料很有限,看paper的时候发现大部分也说得有些简略,所以这里看完之后总结一下。 奇异谱分析(Singular Spectrum Analysis, SSA)是一种处理非线性时间序列数据的方法,通过对所要研究的时间序列的轨迹矩阵进行分解、重构等操作,提取出时间序列中的不同成分序列(长期趋势,季节趋势,噪声等),从而进行对时间序列进行分析或去噪并用于其他一些任务。 奇异谱分析主要包括四个步骤:嵌入——分解——分组——重构。 1. 嵌入SSA的分析对象是有限长一维时间序列 [ x 1 , x 2 , . . . , x N ] [x_1, x_2,...,x_N] [x1,x2,...,xN], N N N 为序列长度。首先需要选择合适的窗口长度 L L L 将原始时间序列进行滞后排列得到轨迹矩阵: X = [ x 1 x 2 ⋯ x N − L + 1 x 2 x 3 ⋯ x N − L + 2 ⋮ ⋮ ⋮ x L x L + 1 ⋯ x N ] \boldsymbol{X}=\left[\begin{array}{cccc}{x_{1}} & {x_{2}} & {\cdots}& {x_{N- L+1}} \\ {x_{2}} & {x_{3}} & {\cdots} & {x_{N-L+2}} \\ {\vdots} & {\vdots} & {} & {\vdots} \\ {x_{L}} & {x_{L+1}} & {\cdots} & {x_{N}}\end{array}\right] X=⎣⎢⎢⎢⎡x1x2⋮xLx2x3⋮xL+1⋯⋯⋯xN−L+1xN−L+2⋮xN⎦⎥⎥⎥⎤通常情况下取 L < N / 2 L ⋯ > λ L ⩾ 0 \lambda_{1}>\lambda_{2}>\cdots>\lambda_{L} \geqslant 0 λ1>λ2>⋯>λL⩾0 和对应的特征向量 U 1 , U 2 , ⋯ , U L U_{1}, U_{2}, \cdots, U_{L} U1,U2,⋯,UL。此时 U = [ U 1 , U 2 , ⋯ , U L ] \boldsymbol{U} =[U_{1}, U_{2}, \cdots, U_{L}] U=[U1,U2,⋯,UL], λ 1 > λ 2 > ⋯ > λ L ⩾ 0 \sqrt{\lambda_{1}}>\sqrt{\lambda_{2}}>\cdots>\sqrt{\lambda_{L}} \geqslant 0 λ1 >λ2 >⋯>λL ⩾0为原序列的奇异谱 。并且有 X = ∑ m = 1 L λ m U m V m T , V m = X T U m / λ m , m = 1 , 2 , . . . , L \boldsymbol{X}=\sum_{m=1}^{L} \sqrt{\lambda_{m}} U_{m} V_{m}^{T}, \quad V_{m}=\boldsymbol{X}^{\mathrm{T}} U_{m} / \sqrt{\lambda_{m}}, \quad m=1,2,...,L X=m=1∑Lλm UmVmT,Vm=XTUm/λm ,m=1,2,...,L 这里 λ i \lambda_{i} λi 对应的特征向量 U i U_{i} Ui 反映了时间序列的演变型,称为时间经验正交函数(T-EOF)。 实际上python已经提供了奇异值分解的函数np.linalg.svd()可以很方便的计算。关于奇异值分解更详细的介绍可以看这篇博客。 3. 分组关于分组,文献中很常见的叙述是下面这样: 简单来说将所有的 L L L 个成分分为 c c c 个不相交的组,代表着不同的趋势成分。这样接下来选择主要的成分进行重构得到重构序列。Emmm。。。。这样介绍可真是太简洁明了导致动手实现的时候真是一脸懵。 因此在实现的时候参考了另一个版本,这里将分组和重构放到一块吧。。。。。这个版本有助于实现但是ran半天ran不清哪里是分组,被自己菜哭。。。。。。。。。。 ![]() 所以这里接分解步。首先计算迟滞序列 X i X_i Xi 在 U m U_m Um 上的投影: a i m = X i U m = ∑ j = 1 L x i + j U m , j , 0 ≤ i ≤ N − L a_{i}^{m}=\boldsymbol{X}_{i} U_m=\sum_{j=1}^{L} x_{i+j} U_{m,j}, \quad 0\leq{i}\leq{N-L} aim=XiUm=j=1∑Lxi+jUm,j,0≤i≤N−L X i X_i Xi 表示轨迹矩阵 X \boldsymbol{X} X 的第 i i i 列, a i m a_{i}^{m} aim 是 X i \boldsymbol{X}_{i} Xi 所反映的时间演变型在原序列的 x i + 1 , x i + 2 , … , x i + L x_{i +1} , x_{i +2} ,…, x_{i +L} xi+1,xi+2,…,xi+L时段的权重, 称为时间主成分(TPC)。看到这里应当发现了,由 a i m a_{i}^{m} aim 构成的矩阵实际上就是没有归一化的右矩阵, 即 λ m V m \sqrt{\lambda_{m}}V_{m} λm Vm ! 接下来就可以通过时间经验正交函数和时间主成分来进行重建,具体重构过程如下: x i k = { 1 i ∑ j = 1 i a i − j k U k , j , 1 ⩽ i ⩽ L − 1 1 L ∑ j = 1 L a i − j k U k , j , L ⩽ i ⩽ N − L + 1 1 N − i + 1 ∑ j = i − N + L L a i − j k E k , j , N − L + 2 ⩽ i ⩽ N x_{i}^{k}=\left\{\begin{array}{l}{\frac{1}{i} \sum_{j=1}^{i} a_{i-j}^{k} U_{k, j}, \quad 1 \leqslant i \leqslant L-1} \\ \\{\frac{1}{L} \sum_{j=1}^{L} a_{i-j}^{k} U_{k, j}, \quad L \leqslant i \leqslant N-L+1} \\ \\ {\frac{1}{N-i+1} \sum_{j=i-N+L}^{L} a_{i-j}^{k} E_{k, j}, \quad N-L+2 \leqslant i \leqslant N}\end{array}\right. xik=⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧i1∑j=1iai−jkUk,j,1⩽i⩽L−1L1∑j=1Lai−jkUk,j,L⩽i⩽N−L+1N−i+11∑j=i−N+LLai−jkEk,j,N−L+2⩽i⩽N 这样,所有重构序列的和应当等于原序列,即 x i = ∑ k = 1 L x i k i = 1 , 2 ⋯ , N x_{i}=\sum_{k=1}^{L} x_{i}^{k} \quad i=1,2 \cdots, N xi=k=1∑Lxiki=1,2⋯,N 通常情况下我们使用SSA只是为了提取原序列的主要成分,以去噪为例,我们只需要根据奇异值的大小选择前 k ( k ≤ L ) k(k \leq L) k(k≤L) 个贡献大的成分重构原序列即可。 python程序 #!/usr/bin/python3 # -*- coding: utf-8 -*- ''' @Date : 2019/11/11 @Author : Rezero ''' import numpy as np import matplotlib.pyplot as plt path = "xxxx" # 数据集路径 series = np.loadtxt(path) series = series - np.mean(series) # 中心化(非必须) # 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 # 重组 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) # 选择重构的部分,这里选了全部 plt.figure() for i in range(10): ax = plt.subplot(5,2,i+1) ax.plot(rec[i, :]) plt.figure(2) plt.plot(series) plt.show()运行程序结果如下,左边是原始序列,右边是按奇异值排序的前十个成分序列,可以看到除了前几个剩余的基本都可以视为噪声序列。 https://www.cnblogs.com/endlesscoding/p/10033527.html 基于SSA的GPS坐标序列去噪及季节信号提取 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |