首页 > 编程学习 > 时间序列转二维图像方法及其应用研究综述

目录

1 前言

2 方法综述

2.1 时频分析法

2.1.1 短时傅里叶变换

2.1.2 小波变换

2.1.3 希尔伯特-黄变换

2.2 图像编码方法

2.2.1 格兰姆角场

2.2.2 马尔可夫转移场

2.2.3 递归图

2.2.4 图形差分场

2.2.5 相对位置矩阵

3 应用研究现状

4 方法仿真与分析

5 总结与分析

6 参考文献


1 前言

时间序列是指按时间的顺序关于事件变化发展的过程记录,它保存了观测数据的时间结构性。因此时间序列常常被当成一个整体进行研究分析,而不是一个个独立的数值。许多工程领域的数据都以时间序列的形式存在,如医学上的心率记录、心电图,公司的财务报表,股票市场价格,交通拥堵指数,还有流行病的传播情况等,这些都是典型的时间序列数据。时间序列有数据量大、维度高、实时更新等特点,针对时间序列数据分析的目的主要有预测、聚类、分类(包括异常检测)、关联分析等。在这些任务中,时间序列分类问题是现有相关研究中备受关注的[2]。

现如今,随着时间序列数据的不断增加,分析目的的不断深入和拓展,使用传统的方法提高时间序列分类问题的精度越来越困难。然而,近年来由于数据量和计算量的不断增加,深度学习快速发展并迅速成为数据研究的重要法宝。其中,卷积神经网络(CNN)是发展最快、技术相对最成熟、且应用最成功的深度学习模型。从经典的LeNet[3]、AlexNet[4]、ResNet[5]到EfficentNet[6],CNN网络结构推陈出新,不断发展。CNN通过共享权重,一方面大幅度减少模型大小,使得网络优化快速迭代,另一方面降低了过拟合风险,有更好的泛化能力。这些使得CNN在多个领域中成功应用,并取得了巨大的效果,尤其在图像领域的分类任务。

考虑到时间序列分类任务面对的困难和CNN在图像分类任务中取得的巨大成功,相关领域的研究者提出了这样的观点:将一维时间序列数据转化为二维图像,把时间序列分类任务转化为图像分类任务,再用CNN模型进行训练学习。从宏观的角度分析,时间序列的分类任务与图像领域的分类任务具有一定的相似性,都是从数据中提取关键特征,再根据特征进行分类。正是这种相似性,促使一些新的研究出现,比如将时间序列数据转化为二维图像的方法研究、将转为方法实际应用到相关领域的研究等等。本文将系统综述时间序列数据转化为二维图像的方法以及其相关应用,最后采用案例研究分析各个方法的优劣。

2 方法综述

在计算机中图像可以用矩阵的形式表示。因此,将时间序列转化为图像等价于将时间序列转化为矩阵。根据这个原则,如表1所示,本论文整理了两大类方法:第一大类为时频分析方法,将时间序列作为信号分析对象,采用时频分析方法分析求解其时频图,主要有短时傅里叶变换[7],小波变换[8],希尔伯特-黄变换[9];第二大类为图像编码方法,这一大类方法指的是通过其他方法对时间序列进行图像编码,主要有格兰姆角场[10],马尔可夫转移场[10],递归图[11],图形差分场[12],相对位置矩阵[13]。下面详细介绍各种方法。

表1:时间序列转二维图像方法归类
归类方法简单描述
时频分析方法短时傅里叶变换选择适当的窗函数;计算窗口内信号的相频信息;推动窗口,得到信号的时频图。
小波变换构造一个快速衰减的母小波;通过缩放和平移得到子小波;进行参数匹配和叠加拟合信号。
希尔伯特-黄变换通过经验模态分解信号得到有限数量的模态函数;对IMF分量进行希尔伯特变换,得到瞬时频率。
图像编码方法格兰姆角场通过空间变换来消除时间序列噪声;采用向量内积来保存时间信息。
马尔可夫转移场将时间顺序看成是马尔可夫过程;分成Q个分位箱,构造马尔可夫转移矩阵,生成马尔可夫转移场。
递归图将时域空间变换位相空间,即每个点变换成相空间中的对应状态;计算状态之间的距离,从而得到相应的图像特征。
图形差分法通过截取图形不同长度,进行差分等处理,构造图形差分场;保留序列熵(复杂性测量、动态系统特征)
相对位置矩阵先进行z-分值标准化,采用PAA进行降维,计算每两个时间戳之间的相对大小关系,最后转换为灰度值矩阵。

2.1 时频分析法

将时间序列看成是一个频率和相位随时间变化的信号,根据时频分析方法,可将时间序列数据转换成时频图,其主要方法有短时傅里叶变换、小波变换、希尔伯特-黄变换。

2.1.1 短时傅里叶变换

短时傅里叶变换(Short Time Fourier Transform, STFT)是傅里叶变换(Fourier Transform, FT)的一种变形,以确认时变信号的频率和相位在时间轴上的变化情况。其基本思想是:选择一个合适的窗函数(常见有方形、三角形、高斯函数等),假设时变信号在这个短时间间隔的窗内是一个平稳的信号,通过傅里叶变换计算出信号的频率和相位信息,然后在时间轴上移动窗函数,计算出不同短时间段内信号的频率和相位,最终实现信号的时频分析,得到其对应的时频图。

G_{f}(\omega, \tau)=\int_{-\infty}^{+\infty} g(t-\tau) f(\tau) e^{-j w t} d \tau

其中, $g(t)$ 为时间窗函数。 根据短时傅里叶变化的定义,给定时变信号 $f(t)$, 其编码步骤如下:

  1. 确定参数:信号长度 f_{len}, 采样频率sample_{rate}, 窗口长度 window_{len }, 步长time_{step}, 窗函数,重叠点数over_{lap}
  2. 计算窗口滑定初始滑动次数为n_{\max }=f i x\left(\frac{f_{len}-over_{lap }}{(window_{len}- over_{lap}) \times time_{step }}\right), 构造一个大小为window_{len} \times n、元素值都为0的初始时频矩阵S
  3. 根据当前窗口位蝠,截取信号片段,进行快速傅里叶变换,得到当前窗口中信号片段所对应的频率分布, 记为列向量P;
  4. 更新当前滑动次数 : n_{now }=n_{now }+1, 更新时频矩阵 S\left(n_{now}\right)=P;
  5. 判断当前滑动次数 n_{now} 是否等于总滑动次数 n_{\max }, 若是 : 编码完成, 输出时频矩阵S; 否则,窗口滑动一个步长, 然后跳回步骤 2 ; 

短时傅里叶变换的实现流程简单, 处理时间短且应用广泛, 但其时频分析后的时间-频 率窗口大小不变, 很难捕捉到一些细小的局部信息。

2.1.2 小波变换

小波变换(Wavelet Transform, WT)是1984年由Morlet和Grossman提出的概念,该方法承袭了短时间窗变换的局部化思想,且克服了时间-窗口大小不变的缺陷,提供了一个窗口宽度可随频率变化而变宽变窄的时频窗,从而充分突出信号的某些特征。其基本思想是:先构造一个有限长或快速衰减的母小波,然后通过缩放和平移生成多个子小波,再叠加以匹配输入信号。将其缩放尺度和平移参数对应频率和时间参数,最终得到信号的时频图。

W_{f}(a, b)=<f, \psi_{a, b}>=\frac{1}{\sqrt{|a|}} \int_{-\infty}^{+\infty} f(t) \psi\left(\frac{t-b}{a}\right) d t

其中, \psi_{a, b} 为母小波函数(Morlet、Ricker 等), a为尺度給数,b为平移給数。 根据小波变换的定义, 给定时变信号f(t), 其编码步骤如下:

  1. 确定参数:信号长度f_{len}, 采样频率sample_{rate}, 母小波函数, 中心频率滑动步长cenfre_{step};
  2. 计算最大中心频率cenfre_{\max}=\frac{sample_{rate}}{2},设置当前中心频率cenfre_{now}=1,初始化时频矩阵S;
  3. 根据中心频率和小波函数,构造小波曲线,再与原信号卷积,得到当前频率的时间分布向量,更新时频矩阵;
  4. 判断当前中心频率是否大于最大中心频率,若是,输出时频矩阵S;否则,更新当前频率cenfre_{now}=cenfre_{now}+cenfre_{step},然后跳回步骤2。

小波变换相较于短时傅里叶变换,具有较好的时频分辨率自适应能力,更能突出实际信号的局部特征,即高频处采用低频率分辨率和高时间分辨率,低频处采用高频率分辨率和低时间分辨率。因而,小波变换在信号处理、语音处理、图像处理等领域得到广泛应用。

2.1.3 希尔伯特-黄变换

前面提到的信号处理方法基本都受到傅里叶理论的影响,不能很好的处理不规则的信号,因此,1998年Norden E. Huang 等人[9]提出经验模态分解方法,并引入Hilbert谱的概念和Hilbert谱分析方法,称为希尔伯特-黄变换(Hilbert-Huang Transform, HHT)。

希尔伯特-黄变换主要包括两个阶段,分别是经验模态分解(EMD)和Hilbert变换(HT)。经验模态分解流程为:

  1. 找到信号的极大值和极小值找到信号f(t)的极大值和极小值,通过三次样条拟合得到上、下包络线,计算其均值得m_1 (t);
  2. 得到第一个分量h_1 (t)=f(t)-m_1 (t), 检查其是否满足模态分量的条件:
    • h_1 (t)得极大值点与过0点数量相差不超过1个;
    • h_1 (t)的上、下包络线均值恒为0。如不满足,重复操作1、2直至得到满足模态函数(IMF)条件的模态分量c_1 (t)
  3. 原始信号减去第一个模态分量,得到信号r_1 (t)=f(t)-c_1 (t),将r_1 (t)当成新的“原始信号”,重复以上操作,直至筛选条件SD=\frac{\sum_{t=0}^{T}\left|h_{k-1}(t)-h_{k}(t)\right|^{2}}{\sum_{t=0}^{T}h_{k-1}^{2}}小于预设值时,经验模态分解结束。这样原始信号便分成若干经验模态分量和一个残余信号:f(t)=\sum_{i=0}^{n} c_{i}+r_{n}(t)

相较于短时傅里叶变换和小波变换,HHT得到时频图分辨率比较低,具有较好的自适应性。

2.2 图像编码方法

2.2.1 格兰姆角场

格拉姆角场(Gramian Angular Field, GAF)是结合坐标变换和格拉姆矩阵的相关知识,实现将时间序列变换成图像的一种编码方法。

格拉姆矩阵是两两向量的内积组成,可以保存时间序列的时间依赖性,却不能有效的区分价值信息和高斯噪声。因此,在进行格拉姆矩阵变换之前,时间序列需要进行空间变换,普遍的方法是将笛卡尔坐标系转换成极坐标系(半径、角度)。

所以对于一个时间序列X=(x_t,t=1,2,...,N),可以通过以下步骤得到GAF图:

  1. 使用最小-最大定标器(Min-Max scaler),将原始时间序列数据缩放到[-1,1];

    \tilde{x}_{-1}^{i}=\frac{\left(x_{i}-\max (X)+\left(x_{i}-\min (X)\right)\right.}{\max (X)-\min (X)}

    \text { or } \quad \tilde{x}_{0}^{i}=\frac{x_{i}-\min (X)}{\max (X)-\min (X)}

  2. 将第一步得到的数据进行极坐标系变换,得到每一个数据点对应的半径和角度:

    \left\{\begin{array}{c} \phi=\arccos \left(\tilde{x}_{i}\right),-1 \leq \tilde{x}_{i} \leq 1, \tilde{x}_{i} \in \tilde{X} \\ r=\frac{t_{i}}{N}, t_{i} \in \mathbb{N} \end{array}\right.

  3. 利用和角关系和差角关系,得到对应的 GASF 图和 GADF 图:

    \begin{aligned} G A S F &=\left[\cos \left(\phi_{i}+\phi_{j}\right)\right] \\ &=\tilde{X}^{\prime} \cdot \tilde{X}-\sqrt{I-\tilde{X}^{2}}^{\prime} \cdot \sqrt{I-\tilde{X}^{2}} \\ G A D F &=\left[\sin \left(\phi_{i}-\phi_{j}\right)\right] \\ &=\sqrt{I-\tilde{X}^{2}}^{\prime} \cdot \tilde{X}-\tilde{X}^{\prime} \cdot \sqrt{I-\tilde{X}^{2}} \end{aligned}

2.2.2 马尔可夫转移场

 马尔可夫转移场(Markov Transition Field, MTF)是基于马尔可夫转移矩阵的一种时间序列图像编码方法。该方法将时间序列的时间推移看成是一个马尔可夫过程,即:在已知目前状态的条件下,它未来的演变不依赖于它以往的演变,由此构造马尔可夫转移矩阵,进而拓展为马尔可夫转移场,实现图像编码。

对于时间序列X=(x_t,t=1,2,...,T),其图像编码步骤如下:

  1. 将时间序列X(t)分成Q个分位箱(标记为1,2,...,Q,每个分位箱内的数据量相同);
  2. 将时间序列中每一个数据更改为其对应的分位箱的序号;
  3. 构造转移矩阵W\omega_{i j}表示分位箱i转移到分位箱j的频率):

    W=\left[\begin{array}{ccc} \omega_{11} & \ldots & \omega_{1 Q} \\ \omega_{21} & \cdots & \omega_{2 Q} \\ \vdots & \ddots & \vdots \\ \omega_{Q 1} & \cdots & \omega_{Q Q} \end{array}\right] \quad \text { s.t. } \quad \sum_{j} \omega_{i j}=1

  4. 构造马尔可夫转移场M

    M=\left[\begin{array}{ccc} \omega_{i j} \mid x_{1} \in q_{i}, x_{1} \epsilon q_{j} & \ldots & \omega_{i j} \mid x_{1} \in q_{i}, x_{N} \in q_{j} \\ \omega_{i j} \mid x_{2} \in q_{i}, x_{1} \epsilon q_{j} & \cdots & \omega_{i j} \mid x_{2} \in q_{i}, x_{N} \in q_{j} \\ \vdots & \ddots & \vdots \\ \omega_{i j} \mid x_{N} \epsilon q_{i}, x_{1} \epsilon q_{j} & \cdots & \omega_{i j} \mid x_{N} \in q_{i}, x_{N} \in q_{j} \end{array}\right]

2.2.3 递归图

递归图(Recurrence Plots, RP)是由Eckmann等人[14]在1995年提出的,用来使动态系统的递归特性可视化。

将递归图应用在时间序列上,首先将时间序列的时域空间变换到相空间,从而将时域中的每个点x_i变换成相空间的对应状态\overrightarrow{S_{i}};接着计算每两个状态(向量)之间的距离(向量范数);然后进行阈值二值化,得到递归图中对应两个状态之间的特征。

递归图可用一系列递归矩阵来表示,如下式所示:

R_{i, j}(\epsilon)=\Theta\left(\epsilon-\left\|\overrightarrow{s_{l}}-\overrightarrow{s_{j}}\right\|\right), i, j=1, \ldots, N

其中R_{i,j}是一个N\times N的方阵,\|\cdot\|示向量范,\epsilon为距离阈值使得R_{i,j} \in\{0,1\}\Theta(\cdot)表示Heaviside函数。

其算法流程如下:

  1. 由时间序列得到相空间状态集;
  2. 计算每两个状态之间的距离(向量范数);
  3. 进行阈值二值化,得到递归图矩阵。

2.2.4 图形差分场

时间序列图形表达了时间序列的时间结构信息,基于时间序列图形的序列熵则可以用于时间序列的复杂性测量、动态系统表征等。因此,作者提出了基于不同图形长度的图形差分场,实现时间序列到图像的转换。

图形差分场(Motif Difference Field, MDF)的基本思想是:给定一个时间序列X=(x_t,t=1,2,...,N),设定需要得到的图像数量n(1<n<N)​​​​​​​,再设定不同步长d(1\leq d\leq d_{\max},d_{\max}=[(N-1)/(n-1)]),不同长度s(s=1,2,3,...,N-(n-1)d)的时间窗口,以此多次提取时间序列原始图形的某一段,然后经过组合变换获得n个图像。具体的流程如下:

  1. 根据时间序列X,得到图形集:
        

    M_{d}^{n}=\left\{M_{d, s}^{n}, s=1,2,3, \ldots, T-(n-1) d\right\}

    其中M_{d, s}^{n}=\left(x_{t}, t=s, s+d, s+2 d, \ldots, s+(n-1) d\right)s相当于一个时间窗口,d相当于步长(时移),t相当于移动窗口(有n个时刻点),t相当于移动窗口(有n个时刻点),共有n\times d_{\max}个图形;
  2. 得到图形差分集:

    d M_{d}^{n}=\left\{d M_{d, s}^{n}, s=1,2,3, \ldots, T-(n-1) d\right\}

    其中d M_{d, s}^{n}=\left(x_{s+d}-x_{s}, x_{s+2 d}-x_{s+d}, \ldots, x_{s+(n-1) d}-x_{s+(n-2) d}\right),共有(n-1) \times d_{\max }个图形;
  3. 定义新序列I_{d,s}^n,通过补零操作使得序列长度一致:

    I_{d, s}^{n}=\left\{\begin{array}{cc} d M_{d, s}^{n}, & 1 \leq s \leq N-(n-1) d \\ 0, & N-(n-1) d<s \leq N-(n-1) \end{array}\right.

  4. 定义图形差分场:

    MDF^{n}=\left\{I_{1}^{n}, I_{2}^{n}, \ldots, I_{d_{m a x}}^{n}\right\}

    I_d^n代表步长d对应的(n-1)个序列集,这样图形差分场可以生成对应的n-1个通道图像。

  5. 针对第i个通道,定义图像数组为:

    G_{i}^{n}=\left[\vec{I}_{1}^{n}(i), \vec{I}_{2}^{n}(i), \ldots, \vec{I}_{d}^{n}(i), \ldots \vec{I}_{d_{m a x}}^{n}(i)\right]^{N}

    其中\vec{I}_{d}^{n}(i)=\left[\vec{I}_{d, 1}^{n}(i), \vec{I}_{d, 2}^{n}(i), \ldots, \vec{I}_{d, T-n+1}^{n}(i)\right]^{N}, 1 \leq i \leq n-1
  6. 针对第i个通道,定义图像数组为:

    G_{i}^{n}=\left[\vec{I}_{1}^{n}(i), \vec{I}_{2}^{n}(i), \ldots, \vec{I}_{d}^{n}(i), \ldots \vec{I}_{d_{m a x}}^{n}(i)\right]^{N}

    其中\vec{I}_{d}^{n}(i)=\left[\vec{I}_{d, 1}^{n}(i), \vec{I}_{d, 2}^{n}(i), \ldots, \vec{I}_{d, T-n+1}^{n}(i)\right]^{N}, 1 \leq i \leq n-1
  7. 填补G_i^n中的零元素,定义MDF图像的每一个通道:

    IMG_{i}^{n}=G_{i}^{n}+K^{n} \odot G_{i}^{\prime n}

    其中K_{d, s}^{n}=\left\{\begin{array}{lr}0, & 1 \leq s \leq N-(n-1) d \\ 1, & N-(n-1) d<s \leq N-(n-1)\end{array}\right.G_{i}^{\prime n}是由G_{i}^{n}旋转而来,\odot为哈达玛积[15]。

 最终实现时间序列到图像的转换。

2.2.5 相对位置矩阵

相对位置矩阵(Relative Position Matrix, RPM)包含了原始时间序列的冗余特征,使转换后的图像中,类间和类内的相似度信息更容易被捕捉。 对于一个时间序列X=(x_t,t=1,2,\ldots,N),可以通过以下步骤得到RPM图:

  1. 针对原始时间序列,通过以下z-分值标准化的方法得到一个标准正态分布Z

    z_{t}=\frac{x_{t}-\mu}{\sigma}, t=1,2, \ldots, N

    其中\mu表示X的平均值,\sigma表示X的标准差。
  2. 采用分段聚合近似(PAA)方法,选择一个合适的缩减因子k,生成一个新的平滑时间序列\tilde X,将维度N减少到m

    \begin{aligned} &\tilde x_{i}=\left\{\begin{array}{l} \frac{1}{k} \sum_{j=k *(i-1)+1}^{k * i} z_{j}, i=1,2, \ldots, m,\left\lceil\frac{N}{k}\right\rceil-\left\lfloor\frac{N}{k}\right\rfloor=0 \\ {\left\{\begin{array}{l} \frac{1}{k} \sum_{j=k *(i-1)+1}^{k * i} z_{j}, i=1,2, \ldots, m-1 \\ \frac{1}{N-k *(m-1)} \sum_{j=k *(m-1)+1}^{N} z_{j}, i=m \end{array},\left\lceil\frac{N}{k}\right\rceil-\left\lfloor\frac{N}{k}\right\rfloor>0\right.} \end{array}\right.\\ &m=\left\lceil\frac{N}{k}\right\rceil \end{aligned}

    通过计算分段常数的平均值进行降维,可以保持原始时间序列的近似趋势,最终新的平滑时间序列\tilde X的长度为m
  3. 计算两个时间戳之间的相对位置,将预处理后的时间序列X转换为二维矩阵M

    M=\left[\begin{array}{cccc} \tilde x_{1}-\tilde x_{1} & \tilde x_{2}-\tilde x_{1} & \cdots & \tilde x_{m}-\tilde x_{1} \\ \tilde x_{1}-\tilde x_{2} & \tilde x_{2}-\tilde x_{2} & \cdots & \tilde x_{m}-\tilde x_{2} \\ \vdots & \vdots & \ddots & \vdots \\ \tilde x_{1}-\tilde x_{m} & \tilde x_{2}-\tilde x_{m} & \cdots & \tilde x_{m}-\tilde x_{m} \end{array}\right]

    如上所示,该矩阵表征了时间序列中每两个时间戳之间的相对位置关系。其每一行和每一列都以某一个时间戳为参考,进一步表征整个序列的信息;
  4. 最后利用最小-最大归一化将M​​​​​​​转换为灰度值矩阵,最终得到相对位移矩阵F

    F=\frac{M-\min (M)}{\max (M)-\min (M)} \times 255

3 应用研究现状

表2:时间序列转二维图像的应用研究
方法文献应用场景主要工作
格拉姆角场[16]交通事件检测采用分段聚合近似(PAA)将车辆速度时间序列数据不同长度转化为同一长度,通过格拉姆角场编码为图像,放入CNN模型中进行训练和检测。
[17]输电线路故障检测和分类采用离散小波变换对故障电压和电流信号进行去噪处理,得到时间序列数据,通过格拉姆角场编码为图像,采用提出的SAT-CNN模型进行训练检测。
[18]基于可穿戴传感器的人类活动识别将传感器数据的时间序列通过格拉姆角场编码为双通道图像,利用CNN模型进行训练识别;采用融合残差网络融合训练异构数据。
[19]输送机电机的预测性维护采用主成分分析方法将多变量时间序列减少到最多两个通道,通过格拉姆角场编码为图像,采用带有PReLU的CNN模型进行预测维护。
[20]动作识别采用格拉姆角场将可佩戴传感器的一维时间序列信号转换为二维图像,然后提出语义感知自适应知识蒸馏网络,提升视觉传感器模态中的动作识别。
[21]表面肌电信号分类过分段聚合近似(PAA)和格拉姆角场将表面肌电信号编码为图像,引入胶囊网络和卷积神经网络进行特征提取和融合,实现不同条件下的手势识别。
[22]金融预测通过格拉姆角场将与标准普尔500指数期货相关的时间序列编码为二维图像,采用简化的VGG网络作为识别网络,输出交易策略(做多、做空、持有)。
马尔可夫转移场[23]住宅建筑能耗异常识别首先将能耗数据进行分割,然后计算基本统计、执行频谱分析或应用信号处理技术进行特征提取,再通过马尔可夫转移场将特征编码为二维图像,采用一类支持向量机完成分类任务。
递归图[24]时间序列预测首先将时间序列转换成递归图,然后使用计算机视觉算法从递归图中提取局部特征,用于预测模型平均。
[25]旋转机械状态监测将三相电流信号的时间序列数据编码为递归图,然后分类器从纹理图像中自动学习特征以分类电机的故障状态。
图形差分场[26]机械故障检测在故障工业机器调查和检查的声音数据集(MIMII数据集)中,采用图形差分场将时间序列编码为图像,再用全连接网络进行识别分类。
相对位置矩阵[13]时间序列分类采用相对位置矩阵将时间序列数据表示为二维图像,基于VGG的改进模型进行分类。

时间序列数据转化为二维图像的目的,是为了更好的提取特征,从而取得更佳的表现效果。如表2所示,是近年来相关研究者在不同领域中应用该类方法的研究现状,主要介绍各项研究的研究场景、使用的方法以及主要工作。在各项研究中,该类方法进行实际应用主要有以下两个问题:

  1. 分析任务多样,如分类、异常检测、预测等,其需要分析的数据长度可能长度不一,则转化后的图像大小可能也不一样。
  2. 以时间序列数据形式存在的场景多样,比如机械状态监控、金融、交通流量等,因此其属性可能有一种,也可能有多种。

针对以上问题,主要有以下解决方案:

  1. 直接应用VS数据预处理:有的数据本身长度同一,可以直接应用时间序列转二维图像的方法;而有的数据长度不一样,则可以用分段聚合近似等方法同一化数据长度,再进行图像转化。此外,有的数据本身存在很大的噪声,或者需要提取出其他有代表性的时序数据特征,这些情况都可以在应用该类方法前进行数据预处理。
    类别文献
    直接应用[18]、[19]、[20]、[22]、[24]、[25]、[26]
    数据预处理[16]、[17]、[21]、[23]、[13]
  2. 单通道VS多通道:时间序列数据的每一种属性都可以通过图像化转化为一张图,作为整体转化后的图当中的一个通道。因而如果只有一种属性,其转化后为一张单通道的图;如果有多种属性,则转化为一张多通道的图。此外,如果属性种类过多且属性相关性较大,也可以通过主成分分析方法进行降维。
    类别文献
    单通道[16]、[17]、[20]、[21]、[22]、[23]、[24]、[25]、[26]、[13]
    多通道[18]、[19]

4 方法仿真与分析

为了进一步分析不同时间序列转图像方法的性能,这里对所有方法进行逐一实现。如图1所示有两条时间序列数据,该数据为广州机场高速公路上某一路段某一段时间内的汽车速度数据,左图为正常情况下的速度数据,右图为发生异常事件情况下的速度数据:

图1:左图为非事件下的数据,右图为事件下的数据。

采用以上详细介绍的各种时间序列转二维图像的方法,其时频图转换结果如图2所示,图编码转换结果如图3所示:

图2:时频图转换方法的结果图,第一栏为异常场景下的图,第二栏为正常场景下的图。
图3:图编码转换方法的结果图,第一栏和第三栏为异常场景下的图,第二栏和第四栏为正常场景下的图

可以看出时频方法与其他方法的实验结果差别很大,时频分析方法的实验结果从图像上看没有明显的差别,这是因为时频分析方法分析的是数据的波动频率,而波动频率不是数据的异常特征;相比之下,其他的方法有明显的图像特征,会出现成块的异常特征,其中GAF、MTF、RP、RPM的实验结果相似。

部分代码为:

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import math

def GAF(X):
    '''
    the code of GAF
    input: the time series in raw
    output: the heatmap of GASF and GADF
    '''
    # normalization to [-1,1]
    X = ((X - np.max(X)) + (X - np.min(X))) / (np.max(X) + np.min(X))

    # generate the GASF
    gasf = X.T * X - np.sqrt(1 - np.square(X)).T * np.sqrt(1 - np.square(X))
    sns.heatmap(gasf, cbar=False, square=True, cmap='GnBu',xticklabels=False, yticklabels=False)
    # plt.show()
    plt.savefig('picture/GASF_1.jpg')

    # generate the GADF
    gadf = np.sqrt(1 - np.square(X)).T * X +  X.T * np.sqrt(1 - np.square(X))
    
    # plot the heatmap
    sns.heatmap(gadf, cbar=False, square=True, cmap='GnBu', xticklabels=False, yticklabels=False)
    # plt.show()
    plt.savefig('picture/GADF_1.jpg')

    return 0

def MTF(X):
    '''
    the code of MTF
    input: the time series in raw
    output: the heatmap of MTF
    '''
    # normalization to [0,1]
    X = (X - np.min(X)) / (np.max(X) - np.min(X))
    # the length of X
    Xlen = X.shape[1]
    # design the number of bins
    Q = 4
    # compute the temp matrix X_Q
    X_Q = np.ones([1, Xlen]) * 4
    # print(X_Q)
    temp = 0
    threshold = np.zeros([1, Q])
    for i in range(Q):
        # print((Xlen * i / Q))
        # print(np.sum(X < temp))
        while np.sum(X < temp) < (Xlen * i / Q):
            temp += 0.01
        # print(threshold.shape)
        threshold[0][i] = temp
        X_Q[np.where(X<temp)] -= 1
        # print(X_Q)
    # print(threshold)

    # generate the Markov matrix
    sum_MM = np.zeros([4,4])

    # compute the probability of Markov
    for i in range(Xlen-1):
        if X_Q[0][i] - X_Q[0][i+1] == -3:
            sum_MM[0][3] = sum_MM[0][3] + 1
        elif X_Q[0][i] - X_Q[0][i+1] == -2:
            if X_Q[0][i] == 1:
                sum_MM[0][2] = sum_MM[0][2] + 1
            elif X_Q[0][i] == 2:
                sum_MM[1][3] = sum_MM[1][3] +1
        elif X_Q[0][i] - X_Q[0][i+1] == -1:
            if X_Q[0][i] == 1:
                sum_MM[0][1] = sum_MM[0][1] + 1
            elif X_Q[0][i] == 2:
                sum_MM[1][2] = sum_MM[1][2] + 1
            elif X_Q[0][i] == 3:
                sum_MM[2][3] = sum_MM[2][3] + 1
        elif X_Q[0][i] - X_Q[0][i+1] == 0:
            if X_Q[0][i] == 1:
                sum_MM[0][0] = sum_MM[0][0] + 1
            elif X_Q[0][i] == 2:
                sum_MM[1][1] = sum_MM[1][1] + 1
            elif X_Q[0][i] == 3:
                sum_MM[2][2] = sum_MM[2][2] + 1
            elif X_Q[0][i] == 4:
                sum_MM[3][3] = sum_MM[3][3] + 1
        elif X_Q[0][i] - X_Q[0][i+1] == 1:
            if X_Q[0][i] == 2:
                sum_MM[1][0] = sum_MM[1][0] + 1
            elif X_Q[0][i] == 3:
                sum_MM[2][1] = sum_MM[2][1] + 1
            elif X_Q[0][i] == 4:
                sum_MM[3][2] = sum_MM[3][2] + 1
        elif X_Q[0][i] - X_Q[0][i+1] == 2:
            if X_Q[0][i] == 3:
                sum_MM[2][0] = sum_MM[2][0] + 1
            elif X_Q[0][i] == 4:
                sum_MM[3][1] = sum_MM[3][1] + 1
        elif X_Q[0][i] - X_Q[0][i+1] == 3:
            sum_MM[3][0] = sum_MM[3][0] + 1
    W = sum_MM
    W = W / np.sum(W,axis=1)
    # print(W)

    # generate the Markov Transform Field
    mtf = np.zeros([Xlen, Xlen])
    for i in range(Xlen):
        for j in range(Xlen):
            # print(X_Q[0][i])
            mtf[i][j] = W[int(X_Q[0][i])-1][int(X_Q[0][j])-1]
    mtf = (mtf - mtf.min()) / (mtf.max()- mtf.min()) * 4
    # generate the heatmap
    sns.heatmap(mtf,cbar=False,square=True,cmap='GnBu',xticklabels=False,yticklabels=False)
    # plt.show()
    plt.savefig('picture/MTF_1.jpg')
    return 0

def RP(X):
    # normalization to [0,1]
    X = (X - np.max(X))  / (np.max(X) + np.min(X))
    Xlen = X.shape[1]
    # convert to the phase space(第一元素是此时高度,第二个给元素为下一时刻的高度)
    S = np.zeros([Xlen-1,2])

    S[:,0] = X[0,:-1]
    S[:,1] = X[0,1:]

    # compute RRP matrix
    R = np.zeros([Xlen-1,Xlen-1])
    for i in range(Xlen-1):
        for j in range(Xlen-1):
            R[i,j] = sum(pow(S[i,:]-S[j,:],2))
    # normalization to [0,4] of RP
    R = (R - R.min()) / (R.max() - R.min()) * 4

    # show the heatmap(bwr,coolwarm,GnBu)
    sns.heatmap(R, cbar=False, square=True,cmap='GnBu',xticklabels=False,yticklabels=False)
    # plt.show()
    plt.savefig('picture/RP_1.jpg')


    return 0

def MDF(X,n):
    # normalization to [0,1]
    X = (X - np.max(X))  / (np.max(X) + np.min(X))

    # compute the length of time series and the range of d
    T = X.shape[1]
    dMax = math.floor((T-1)/(n-1))

    # initial the M,IMG
    M = np.zeros([n,T-n+1,dMax])
    
    
    # initial the dM,K
    dM = np.zeros([n-1,T-n+1,dMax])
    K = np.ones([n-1,T-n+1,dMax])   
    for d in range(dMax):
        d = d + 1
        # initial the s
        s = np.zeros([T-(n-1)*d])
        for i in range(T-(n-1)*d):
            s[i] = i
        for ImageIndex in range(n):
            # print(s+ImageIndex*d)
            s_index = (s+ImageIndex*d).astype(np.int16)
            s = s.astype(np.int16)
            # print(X[0,s_index])
            M[ImageIndex,s,d-1] = X[0,s_index]

            
            if ImageIndex >= 1:
                # motif difference
                dM[ImageIndex-1,s,d-1] = M[ImageIndex,s,d-1] - M[ImageIndex-1,s,d-1]
                # K
                K[ImageIndex-1,s,d-1] = np.zeros([T-(n-1)*d])
    
    IMG = np.zeros([n-1,T-n+1,dMax])
    for ImageIndex in range(n-1):
        # G
        G = dM[ImageIndex]
        # the rot180 of G
        G_1 = G.reshape(1,(T-n+1)*dMax)
        G_2 = G_1[0][::-1]
        G_rot = G_2.reshape(T-n+1,dMax)
        IMG[ImageIndex,:,:] = G + K[ImageIndex] * G_rot
        sns.heatmap(IMG[ImageIndex,:,:].T, cbar=False, square=True, cmap='GnBu', xticklabels=False, yticklabels=False)
        # plt.show()
        # print(IMG)
        plt.savefig('picture/MDF'+'%d'%ImageIndex+'_1.jpg')
    return 0

if __name__ == "__main__":

    # get the data
    data = pd.read_csv(filename)

    # dataframe to numpy
    data_numpy = data.values

    # show the motif of the data
    fig1 = plt.figure(1)
    plt.plot(range(len(data_numpy)), data_numpy)
    plt.xlim([0, len(data_numpy)-1])
    plt.ylim([20, 120])
    plt.ylabel('speed')
    plt.xlabel('time index')
    plt.title('time series of normal')
    # plt.show()
    plt.savefig('picture/normal.jpg')

    # GAF(data_numpy.T)

    # MTF(data_numpy.T)

    # RP(data_numpy.T)

    n = 4
    MDF(data_numpy.T,n)

    # x=list()
    # temp = list()
    # for t in np.arange(0,10,0.01):
    #     temp.append(t)
    #     x.append(math.sin(t))
    # print(x)
    # plt.plot(temp,x)
    # plt.show()

5 总结与分析

本文详细介绍了目前时间序列转二维图像的方法及其应用现状;同时以交通异常检测场景为例,对各个方法做了案例分析。各项研究结果表明,将时间序列数据转化为二维图像,利用成熟的计算机视觉技术进行特征提取和识别,可以有效提高效果,这将有助于时间序列数据的研究。总的来说,此类方法之所以更优于传统的时间序列分析方法,主要是因为有以下几点:

  1. 转化的过程本质上是一种升维的过程,这意味着信息量增加,比如增加了不同时刻数据值之间的关系表征;
  2. 转化的过程同时也是一种数据处理过程,凸显了数据的主要特征;
  3. 可以借助计算机视觉领域中的特征提取机制进行高效地特征提取;

如今,随着数据采集技术的不断提升和数据采集设备的广泛铺设,众多场景下的时间序列数据呈大体量涌现,急需高效的数据分析方法。而将时间序列转化为二维图像,借助现有的图像领域的研究成果进行分析,是一个新的突破途径。该领域还有一些待挖掘的研究点,比如转化过程中的数据不丢失问题。目前,时间序列分析领域相对于图像领域,在研究进展上还未有较大的突破,希望本文可以提供给相关领域的研究者一个新的视角。

6 参考文献


[1] Fu T. A review on time series data mining[J]. Engineering Applications of Artificial Intelligence, 2011, 24(1): 164-181.

[2] Fawaz H I, Forestier G, Weber J, et al. Deep learning for time series classification: a review[J]. Data mining and knowledge discovery, 2019, 33(4): 917-963.

[3] LeCun Y, Bottou L, Bengio Y, et al. Gradient-based learning applied to document recognition[J]. Proceedings of the IEEE, 1998, 86(11): 2278-2324.

[4] Krizhevsky A, Sutskever I, Hinton G E. Imagenet classification with deep convolutional neural networks[J]. Advances in neural information processing systems, 2012, 25: 1097-1105.

[5] He K, Zhang X, Ren S, et al. Deep residual learning for image recognition[C]//Proceedings of the IEEE conference on computer vision and pattern recognition. 2016: 770-778.

[6] Tan M, Le Q. Efficientnet: Rethinking model scaling for convolutional neural networks[C]//International Conference on Machine Learning. PMLR, 2019: 6105-6114.

[7] Gabor D. Theory of communication. Part 1: The analysis of information[J]. Journal of the Institution of Electrical Engineers-Part III: Radio and Communication Engineering, 1946, 93(26): 429-441.

[8] Morlet J, Arens G, Fourgeau E, et al. Wave propagation and sampling theory—Part II: Sampling theory and complex waves[J]. Geophysics, 1982, 47(2): 222-236.

[9] Huang N E, Shen Z, Long S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the Royal Society of London. Series A: mathematical, physical and engineering sciences, 1998, 454(1971): 903-995.

[10] Wang Z, Oates T. Imaging time-series to improve classification and imputation[C]//Twenty-Fourth International Joint Conference on Artificial Intelligence. 2015.

[11] Hatami N, Gavet Y, Debayle J. Classification of time-series images using deep convolutional neural networks[C]//Tenth international conference on machine vision (ICMV 2017). International Society for Optics and Photonics, 2018, 10696: 106960Y.

[12] Zhang Y, Chen X. Motif Difference Field: A Simple and Effective Image Representation of Time Series for Classification[J]. arXiv preprint arXiv:2001.07582, 2020.

[13] Chen W, Shi K. A deep learning framework for time series classification using Relative Position Matrix and Convolutional Neural Network[J]. Neurocomputing, 2019, 359: 384-394.

[14] Eckmann J P, Kamphorst S O, Ruelle D. Recurrence plots of dynamical systems[J]. World Scientific Series on Nonlinear Science Series A, 1995, 16: 441-446.

[15] Davis C. The norm of the Schur product operation[J]. Numerische Mathematik, 1962, 4(1): 343-344.

[16] Liu X, Cai H, Zhong R, et al. Learning traffic as images for incident detection using convolutional neural networks[J]. IEEE Access, 2020, 8: 7916-7924.

[17] Fahim S R, Sarker Y, Sarker S K, et al. Self attention convolutional neural network with time series imaging based feature extraction for transmission line fault detection and classification[J]. Electric Power Systems Research, 2020, 187: 106437.

[18] Qin Z, Zhang Y, Meng S, et al. Imaging and fusing time series for wearable sensor-based human activity recognition[J]. Information Fusion, 2020, 53: 80-87.

[19] Fahim S R, Sarker Y, Sarker S K, et al. Self attention convolutional neural network with time series imaging based feature extraction for transmission line fault detection and classification[J]. Electric Power Systems Research, 2020, 187: 106437.

[20] Kiangala K S, Wang Z. An effective predictive maintenance framework for conveyor motors using dual time-series imaging and convolutional neural network in an industry 4.0 environment[J]. Ieee Access, 2020, 8: 121033-121049.

[21] 骆俊锦, 王万良, 王铮, 等. 基于时序二维化和卷积特征融合的表面肌电信号分类方法[J]. 模式识别与人工智能, 2020, 33(7): 588-599.

[22] Barra S, Carta S M, Corriga A, et al. Deep learning and time series-to-image encoding for financial forecasting[J]. IEEE/CAA Journal of Automatica Sinica, 2020, 7(3): 683-692.

[23] Fahim M, Fraz K, Sillitti A. TSI: Time series to imaging based model for detecting anomalous energy consumption in smart buildings[J]. Information Sciences, 2020, 523: 1-13.

[24] Li X, Kang Y, Li F. Forecasting with time series imaging[J]. Expert Systems with Applications, 2020, 160: 113680.

[25] Hsueh Y, Ittangihala V R, Wu W B, et al. Condition monitor system for rotation machine by CNN with recurrence plot[J]. Energies, 2019, 12(17): 3221.

[26] Zhang Y, Gan F, Chen X. Motif Difference Field: An Effective Image-based Time Series Classification and Applications in Machine Malfunction Detection[C]//2020 IEEE 4th Conference on Energy Internet and Energy System Integration (EI2). IEEE, 3079-3083.


Copyright © 2010-2022 dgrt.cn 版权所有 |关于我们| 联系方式