时间序列如何分析周期性？

乔帮主

如何简单地分析一个时间序列的周期呢？

比如网络流量对时间的关系，很可能是有日周期的。那么怎么简单地根据一个时间序列来分析出它的周期性呢？

举个例子来说：

1, 11, 1, 10, 1, 12, 1, 11, 1, 10, 1, 12, 1, 11, 1, 10, 1, 12,1, 11, 1, 10, 1, 12
这个序列的周期确切为6， 但是看起来，周期似乎为2.

我需要如何得到2
原文地址：https://www.zhihu.com/question/23654854

继续前进

论文链接：

TimesNet: Temporal 2D-Variation Modeling for General Time Series...本文中了2023 ICLR，是清华软院龙明盛老师组的文章，一如既往的Solid。本文的作者就是Autoformer的作者，所以本文的很多思想都延续了Autoformer。Autoformer在知乎有作者团队官方的解析，如下：
游凯超：Autoformer:基于深度分解架构和自相关机制的长期序列预测模型不同于Autoformer只集中于时间序列预测，本文提出的TimesNet是一个通用的时间序列神经网络骨干，可处理各种不同的时间序列任务，如最常见的任务：预测、分类、异常检测等等。其实几乎所有的时间序列预测模型也可以当做是通用骨干，比如Autoformer，Informer，FEDformer，Preformer这些Transformer-based模型中只采用Encoder就相当于是一个时间序列的特征提取器，区别在于它们捕获时序依赖性的方式不同。比如Autoformer是用Auto-Correlation，Informer中的概率稀疏Attention，FEDformer的频域Attention，Preformer中的Multi-Scale Segment-Correlation。还有那些MLP-based模型比如DLinear也可以当做是通用骨干，它是直接采用线性层权重来表示时序依赖性。
Key Points

1D变2D

这是本文的核心。大部分现有方法都是作用于时间序列的时间维度，捕获时序依赖性。实际上，现实时间序列一般都有多种模式，比如不同的周期，各种趋势，这些模式混杂在一起。如果直接对原始序列的时间维度来建模，真正的时序关系很可能隐藏在这些混杂的模式中，无法被捕获。考虑到：现实世界的时间序列通常具有多周期性，比如每天周期、每周周期、每月周期；而且，每个周期内部的时间点是有依赖关系的（比如今天1点和2点），不同的相邻周期内的时间点也是有依赖关系的（比如今天1点和明天1点），作者提出将1D的时间维度reshape成2D的，示意图如下。下图左侧的时间序列具有三个比较显著的周期性（Period 1、Period 2、Period 3），将其reshape成三种不同的2D-variations，2D-variations的每一列包含一个时间段（周期）内的时间点，每一行包含不同时间段（周期）内同一阶段的时间点。变成2D-variations之后，就可以采用2D卷积等方式来同时捕获时间段内部依赖和相邻时间段依赖。


那么怎么确定时间序列中的周期性呢？采用傅里叶变换。给时间序列做傅里叶变换后，主要的周期会呈现对应的高幅值的频率分量。设定超参数k，然后只取top k个最大的幅值对应的频率分量，即可得到top k个主要的周期，这和Autoformer中的处理类似。具体操作如下图，左侧是确定top k个周期，在此只画了三个，然后将1D的时间序列reshape成3种不同的2D-variations（不能整除的可以用padding），对这三种2D-variations用2D卷积进行处理之后再聚合结果即可。


一般来说，对于一个多变量时间序列 ，其中 是变量维数， 是长度，虽然它是一个2D tensor，但作者将其称为是1D的，这是因为在时间维度上来看是1D的。可以通过上图中这种方式，先算出主要周期和频率，再根据主要周期和频率将时间维度上是1D的时间序列reshape成k个2D-variations。注意，对于 个变量，最终算得的主要周期是所有变量的主要周期的平均，这也说明输入的多变量时间序列中包含的不同单变量时间序列的周期模式需要相似。最后，第i个2D-variations即是 ，其中 和 分别表示第i个周期和频率，它们的关系如下式：


TimesBlock

得到k个2D-variations之后该怎么处理呢？本文提出了TimesBlock，每层TimesBlock又分为两步。首先是要先对这些2D-variations分别用2D卷积（可以是ResNet、ConvNeXt等）或者其他的视觉骨干网络（比如Swin，Vit）处理；其次将k个处理后的结果再聚合起来。
对于第一步，本文采用了一种参数高效的Inception block。Inception block是GoogleNet中的模块，包含多个尺度的2D卷积核。如下图左侧蓝色区域，处理k个2D-variations的Inception block是参数共享的。因此，模型整体的参数量不会随着超参数k的增大而增大，因此本文将其称为参数高效的Inception block（Parameter-efficient Inception block）。


对于第二步，在处理完k个2D-variations之后，需要将其展平回1D-variations，并截断到原始长度 （这对应于前面不能整除时使用padding的情况，相当于把多余的padding给去掉）。总之，得到k个变换回去的1D-variations之后，该如何聚合这k个结果呢？如上图右侧所示，也是延续Autoformer的思路，根据傅里叶变换后频率周期对应的赋值大小来加权聚合，幅值大的证明该频率周期的分量越显著也越重要，给它较大的聚合权重，幅值小的则相反。直接用softmax归一化这些幅值 ，然后将归一化后幅值作为加权权重来聚合上面得到的k个1D-variations即可：


实验结果

作者在五种时间序列任务上做了实验，充分对比了一些其他的时间序列骨干。五边形战士：


作者也用了不同的视觉骨干来处理2D-variations：


在长时间序列预测上的效果：


Comments

文章真的写的很好，idea很清晰合理，实验很充分效果也很不错，在长时间序列预测上超越了很多很先进的Transformer-based模型和MLP-based模型。有些新中2023 ICLR的论文在长时间序列预测上的效果非常差，甚至是一些时序预测任务上中了oral的文章，写的花里胡哨，创新性也没有特别显著，常看这个领域的基本看一遍那些文章就知道大概啥水平，效果也不能打，根本不实用。

检验医师

从信号处理角度进行分析
简单的时间序列直接做各种谱分析（频谱，包络谱，平方包络谱，功率谱，倒谱等等）
比如一些简单的旋转机械振动时间序列信号


如果频谱不好分析，那可以分析如下图所示的时间序列的时频谱


给个简单的模拟信号的例子
t = 0:1/2000:1-1/2000;
dt = 1/2000;
x1 = sin(50*pi*t).*exp(-50*pi*(t-0.2).^2);
x2 = sin(50*pi*t).*exp(-100*pi*(t-0.5).^2);
x3 = 2*cos(140*pi*t).*exp(-50*pi*(t-0.2).^2);
x4 = 2*sin(140*pi*t).*exp(-80*pi*(t-0.8).^2);
x = x1+x2+x3+x4;
figure;
plot(t,x)
title('Superimposed Signal')

其连续小波变换时频谱如下


一个模拟的轴承内圈故障振动信号，带有明显的周期性


相应的频谱如下，红色虚线代表故障特征频率及相应的倍频


包络谱如下


看一下相应的CWT时频谱，很明显能看出冲击性


还可以试试小波相干与交叉小波分析
小波相干、交叉小波，可以很好地反映两个不同时间序列变化之间的“相关性”。小波相干分析，一般反映序列间周期性“变化趋势”的一致性，但不直接反映变化周期的强度关系；交叉小波分析，一般反映序列间“共有周期”的强度。






此外，如果时频谱线能量发散，时频脊线模糊，还可以试试同步压缩之类的算法
时间序列信号处理系列-基于Python的同步压缩变换 - 哥廷根数学学派的文章 - 知乎 https://zhuanlan.zhihu.com/p/554189692
当时间序列信号中噪声较大时，为了有利于周期性分析，不可避免地要进行降噪前处理
比如K-SVD降噪


样条框架降噪


Morlet小波降噪










当待分析的时间序列过于复杂时，那可能要引入模态分解（多分辨分析），比如小波分解，经验模态分解及其变体，变分模态分解，经验小波变换，局部均值分解，辛几何模态分解，各种各样的自适应分解算法
基于小波脊线的时间序列分解










好多同学都对各种模态分解方法的时间序列处理感兴趣，那就随便说一下
实际上，时间序列通常由多个具有物理意义的分量组成，在很多时候，为了更容易的研究信号，我们希望在与原始数据相同的时间尺度上单独研究这些分量中的一个或多个，理想情况下，我们希望这些经MRA分解到的多个分量在物理上是有意义的，可容易解释的。多分辨率分析MRA通常与小波或小波包相关联，但诸如经验模态分解EMD，变分模态分解VMD等模态分解方法也可以构成MRA。
先给一个简单的合成信号，信号以1000Hz的频率采样1秒钟。
Fs = 1e3;
t = 0:1/Fs:1-1/Fs;
comp1 = cos(2*pi*200*t).*(t>0.7);
comp2 = cos(2*pi*60*t).*(t>=0.1 & t<0.3);
trend = sin(2*pi*1/2*t);
rng default
wgnNoise = 0.4*randn(size(t));
x = comp1+comp2+trend+wgnNoise;
plot(t,x)
xlabel('Seconds')
ylabel('Amplitude')
title('Synthetic Signal')

该信号由3个主要分量组成：频率为 60Hz的时间局部振荡分量、频率为 200 Hz的时间局部振荡分量和趋势项分量。趋势项分量为正弦曲线，频率为0.5Hz。60Hz的振荡分量发生在 0.1到 0.3 秒之间，而 200Hz的振荡分量发生在 0.7 到 1 秒之间。
但这些分量从时域波形中无法分辨，因此进行频域变换。
xdft = fft(x);
N = numel(x);
xdft = xdft(1:numel(xdft)/2+1);
freq = 0:Fs/N:Fs/2;
plot(freq,20*log10(abs(xdft)))
xlabel('Cycles/second')
ylabel('dB')
grid on

从频率中可以更容易地辨别振荡分量的频率，但时间局部性信号却丢失。为了同时定位时间和频率信息，使用连续小波变换进行分析。


从CWT时频谱图中可以看出60Hz和200Hz分量的时间范围，但没有发现趋势项分量。为了分离出信号的分量并单独进行分析，接下来使用多分辨分析，直接在时域中进行相关操作。
多分辨分析通过将信号分成不同分辨率的分量进而缩小分析范围，而提取不同分辨率的信号分量相当于分解数据在不同时间尺度上的变化，或等效地在不同频带上进行分析。首先，采用离散小波变换的变体最大重叠离散小波变换对信号进行多分辨分析，分解层数为8。关于最大重叠离散小波变换的相关内容，请查看如下文献。


最大重叠离散小波变换的8层多分辨分析分解如下：


如果从上向下看，会看到所分解的分量变得越来越平滑，即分量频率越来越低。回想一下，原始信号包含3个主要分量，一个 200 Hz 的高频振荡成分、一个 60 Hz 的低频振荡成分和一个趋势成分，它们都被加性噪声破坏了。
从D2 图中可以看出时间局部化的高频分量被分解出来，而下面的两个图包含较低频率的振荡分量，这是多分辨率分析的一个重要方面，最后S8子图中包含了趋势项分量。
除了小波多分辨分析，经验模态分解 (EMD) 是一种所谓的数据自适应多分辨技术。 EMD 在不使用固定基函数的情况下递归地从数据中提取不同的分辨率成分，关于EMD相关文献浩如烟海，不做赘述了。EMD的多分辨分析分解如下所示：


虽然MRA分解分量的数目不同，但 EMD MRA和小波 MRA会产生相似的信号波形，在 EMD MRA分解中，高频振荡成分位于第1个本征模态函数中 (IMF1)，低频振荡成分主要位于IMF2和IMF3中，IMF6 中的趋势项分量与小波技术提取的趋势分量非常相似。
自适应多分辨分析的另一种技术是变分模态分解 (VMD)，VMD 从信号中提取固有模式函数或振荡模式，并不使用固定基函数进行分析。EMD在时域上递归，以逐步提取低频IMF分量，而VMD 首先识别频域中的信号峰值并同时提取所有模式，相关文献如下：
Dragomiretskiy, Konstantin, and Dominique Zosso. “Variational Mode Decomposition.” IEEE Transactions on Signal Processing 62, no. 3 (February 2014): 531–44. https://doi.org/10.1109/TSP.2013.2288675.
VMD的多分辨分析分解如下所示：


由上图可知，与小波和EMD类似，VMD将3个分量基本分离了出来。
还有一种数据自适应多分辨分析技术：经验小波变换 (EWT) ，EWT根据分析信号的频率构造 Meye小波进而进行自适应小波，之前写过EWT相关的内容：
经验小波变换在信号处理及轴承故障诊断中的应用 - 哥廷根数学学派的文章 - 知乎https://zhuanlan.zhihu.com/p/53
EWT的多分辨分析分解如下所示：


与之前的EMD和小波MRA类似，EWT分解出了相关的振荡分量，用于执行分析的滤波器及其通带信息如下：


下面考虑一段神户地震信号，源于1995 年 1 月 16 日在澳大利亚霍巴特的塔斯马尼亚大学记录，从 20:56:51 (GMT) 开始，以 1 秒的间隔持续 51分钟。
figure
plot(T,kobe)
title('Kobe Earthquake Seismograph')
ylabel('Vertical Acceleration (nm/s^2)')
xlabel('Time')
axis tight
grid on

以最大重叠离散小波变换为例，其8层MRA分解如下：


从D4和D5子图中可以看出初级与延迟次级波分量，地震波中的分量以不同的速度传播，初级波比次级（剪切）波传播的更快。
将信号分解为若干分量的目的通常是去除某些分量以减轻对信号分析的影响，MRA技术的关键是重建原始信号的能力，如下：


每种方法的最大重建误差约为10^(-12) 或更小，表明它们可以完美地对信号进行重建。在很多研究中我们对趋势项不感兴趣，由于趋势项一般位于最后一个 MRA 分解分量中，因此只需将该分量去除，然后进行重建。


此外，再删除第1个MRA分解分量（看起来主要是噪声）


在前面我们将趋势项删除，然而在许多应用中，趋势项可能是我们的主要研究部分，因此可视化几种MRA方法所提取的趋势项分量。


根据以上的分析，小波MRA技术可以更平滑且最准确地提取趋势项，EMD提取了一个平滑的趋势项，但它相对于真实趋势幅度发生了偏移，而 VMD似乎比小波和EMD更偏向于提取振荡分量。
在前面的示例中，强调了多分辨分析在检测数据中的振荡分量和总体趋势中的作用，然而MRA还可以定位和检测信号中的瞬态成分。为了说明这一点，以1947年第一季度至 2011 年第四季度美国实际国内生产总值 (GDP) 数据，垂直的黑线标志着“大缓和”的开始，标志着从 1980 年代中期开始，美国宏观经济波动性减弱的时期，很难从原始数据中辨别出来。

长长的路

时间序列预测的主流模型结构一直以来被RNN、CNN、Transformer三大模型主体主导。而Nbeats的出现，让纯全连接的模型结构在时间序列预测问题上也能取得非常好的效果。这篇文章梳理了Nbeats系列工作，从最基础的Nbeats版本，到可以引入外部变量的Nbeats版本，再到能够处理时空预测的Nbeats版本。


1
基础版本Nbeats
论文题目：N-BEATS: NEURAL BASIS EXPANSION ANALYSIS FOR INTERPRETABLE TIME SERIES FORECASTING
下载地址：https://arxiv.org/pdf/1905.10437.pdf

Nbeats是Element AI发表于ICLR 2020的一篇工作，目前引用量200多，在时间序列预测这个领域还是比较有影响力的。Nbeats开创了一个全新的时间序列预测backbone，仅通过全连接实现时间序列预测。Nbeats的核心思路是，通过多层全连接进行时间序列分解，每层拟合时间序列部分信息（之前层拟合的残差），有点类似于GBDT的啥思路。
Nbeats的具体模型结构如下图所示。整个模型包括多个stack，每个stack包括多个block，每个block是Nbeats的最基础结构模块，由多个全连接层组成。每个block包含两个主要部分，第一个部分将输入的时间序列映射成expansion coefficients，第二部分将expansion coefficients映射回时间序列。




什么是expansion coefficients呢？expansion coefficients可以理解为存储了时间序列内在的信息形成的一个低维向量。在模型实现上，其实就是一个向量映射过程：将输入的时间序列（维度为length）映射成低维向量（维度为dim），第二部分再将其映射回时间序列（长度为length）。这个步骤也类似于AutoEncoder，将时间序列映射成一个低维向量保存核心信息，再还原回来。假设每个block模块的输入为
，将其映射为expansion coefficients的过程可以表示为




每个模块会生成两组expansion coefficients，一组用来预测未来（forecast），另一组用来预测过去（backcast）。这个过程可以表示为如下公式：




最终，每个block对输入的序列进行处理后，输出一个预测未来的序列，以及一个预测过去的序列。每个block的输入，是上一层block的输入减去上一层block的输出。通过这种方式，模型每层需要处理的是之前层无法正确拟合的残差，也起到了一个将时间序列进行逐层分解，每层预测时间序列一部分的作用。最终的预测结果，是各个block预测结果的加和。

为了能让模型的分解具有可解释性，文中也提出了在各个层引入一些先验知识，强制让某些层学习某种类型的时间序列特性，实现可解释的时间序列分解。实现的方法是通过约束expansion coefficients到输出序列的函数形式来实现。例如想让某层block主要预测时间序列的季节性，就可以用下面的公式强制输出是季节性的：




下图是作者用这种思路约束不同层学习不同信息的可视化结果，有的层学习了趋势性，有的层学习了周期性。




2
引入外部变量的Nbeats
论文题目：Neural basis expansion analysis with exogenous variables:Forecasting electricity prices with NBEATSx
下载地址：https://arxiv.org/pdf/1905.10437.pdf

第一版本的Nbeats，输入只能是单一的时间序列，无法输入额外的特征。而在时间序列预测问题中，诸如日期信息、节日信息、属性信息等外部特征也是非常重要的。因此基于初版Nbeats，该团队又提出了可以引入外部特征的Nbeatsx，和初版Nbeats的主要区别是引入了外部特征X。




模型的主体结构和Nbeats基本一致，每个block除了输出序列外，还会输入外部特征，二者一起通过全连接层得到隐状态，再基于隐状态生成expansion coefficients。




此外，文中还提出了另一种引入time-dependent特征的方法：采用一个TCN（时间序列卷积模块）作为encodeer，对外部特征进行编码，将该编码作为生成预测结果的因素：




3
用于时空预测的Nbeats
论文题目：GAGA: Fully Connected Gated Graph Architecture for Spatio-Temporal Traffic Forecasting
下载地址：https://arxiv.org/pdf/2007.15531.pdf

该团队提出的第三个版本的Nbeats将Nbeats扩展到了时空预测领域，能够处理存在空间关系的多个时间序列的建模。模型总共分为Graph Edge Weight、Time Gate、Graph Gate，以及和Nbeats类似的多block、stack嵌套的主体全连接结构，模型整体结构如下图所示。




Graph Edge Weight：将每个节点都表示成一个embeddiing，对于两个节点的关系，使用这两个节点对应的embedding内积表示，最终可以用一个矩阵W表示整个空间上两两节点之间的关系：




Time Gate Block（紫色部分）：Time Gate模块主要用来对时间特征进行编码。这里会将时间特征和每个节点的embedding进行拼接后，通过全连接生成时间特征相关的表示。这个时间特征信息会通过乘法、除法的方式，实现将时间信息从原始序列中剥离（对应除法）再融合（对应乘法）的目的。Time Gate的输出对应两个Linear映射结果，一个用于从历史序列中剥离时间因素，一个用于最终的还原，因为历史序列和未来序列处于不同时间窗口，因此这里采用两套参数分别建模时间信息。Time Gate这种将时间因素从序列中剥离的方法，让模型能够更专注于非时间因素的纯序列建模，降低模型拟合难度。

Graph Gate Block（绿色部分）：Graph Gate的目的是让每个节点都能融合整个空间中其他节点的信息。核心思路是，对于每一个节点，利用它和其他节点的Graph Edge Weight融合其他节点的时间序列信息。Graph Gate模块输出一个矩阵G，矩阵中每个元素的计算过程如下，相当于用节点i和节点j之间的距离对节点j的时间序列的每个时刻k加权：




由于使用了节点i的最大值进行减法和除法的处理后再过Relu激活函数，对于不相关的节点对，Relu激活函数可以起到过滤的作用，让这两个节点的输出结果为0。

上面生成的矩阵G、每个节点的embedding以及每个节点的序列，最终作为后续Nbeats模型的输入。Time Gate部分对输入进行信息分离，再在最终的产出结果上进行时间信息融合。

4
总结
本文介绍了Nbeats的序列模型，包括最基础的Nbeats模型，以及在此基础上衍生出来的引入外部特征的Nbeats（Nbeatsx）和用于时空预测的Nbeats（GAGA）。Nbeats相比其他时间序列预测模型，独创了一种全部为全连接的backbone，核心思路是通过序列信息分解、分而治之的方法，实现准确的时间序列预测。

长长的路

一种方法是傅立叶变换，该方法快速且计算效率高(O(nlogn))，但在存在多季节性模式时效果不佳。 另外，它对任何丢失的数据都非常敏感。 与傅立叶变换相比，Serial correlation （也称信号自相关，自相关图或ACF）既更准确，又对丢失数据不那么敏感（不过在计算上也更昂贵）。
附1：
from statsmodels import api as sm
import numpy as np
import scipy.io as sio
import numpy as np
import matplotlib.pyplot as plt
from scipy import fft, arange, signal
tempNorm=np.array([1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4])
acf = sm.tsa.acf(tempNorm, nlags=len(tempNorm))
plt.figure(figsize = (10, 8))
lag = arange(len(tempNorm))
plt.plot(lag, acf)
plt.xlim((0, 10))
plt.xlabel('Lags (days)')
plt.ylabel('Autocorrelation')

除0之外的尖峰都是周期，最尖峰是最小正周期
附2：
时间序列数据集网站https://www.cs.ucr.edu/~eamonn/time_series_data_2018/
时间序列数据挖掘百问百答https://www.cs.ucr.edu/~eamonn/100_Time_Series_Data_Mining_Questions__with_Answers.pdf
附3：
Anodot公司一家做时间序列异常检测的公司则使用自己的专有算法“ Vivaldi”检测周期性。Vivaldi的基础是ACF，但使用智能子采样以减少计算量。 通过对同一时间序列的多个过滤版本（multiple filtered versions）应用此方法，可以检测到多个季节模式。
https://www.anodot.com/blog/closer-look-time-series-anomaly-detection/

卡卡

楼上有人提到频谱图，我弄一下
假设抽样频率为1000Hz，频谱图如下


看出主要频率成分是500Hz，即以每两个数据为一个周期
matlab代码如下

Fs=1000;
tmp=[1, 11, 1, 10, 1, 12, 1, 11, 1, 10, 1, 12, 1, 11, 1, 10, 1, 12,1, 11, 1, 10, 1, 12];
tmp=2*(tmp-min(tmp))/(max(tmp)-min(tmp))-1;%归一化到[-1,1]
plot(tmp)
L=length(tmp);
NFFT=2^nextpow2(L);
Y=fft(tmp,NFFT)/L;
f=Fs/2*linspace(0,1,NFFT/2+1);
figure
plot(f,2*abs(Y(1:NFFT/2+1)));
title('幅度谱')
xlabel('频率 Hz');
ylabel('|Y(f)|')