基于自适应性分段统计近似的时间序列相似性度量方法

摘要

本发明公开了一种基于自适应性分段统计近似的时间序列相似性度量方法，首先基于时间序列编码识别转折点，将时间序列切分为包含完整波动趋势的子序列；然后依次提取每条子序列的多种统计特征，构造局部模式特征向量；最后利用规范化距离计算局部模式特征向量之间的距离，实现局部模式匹配，并以此作为动态规划算法的子程序，实现全局模式匹配。本发明在度量精度和计算效率方面都以较大的程度优于其他度量方法，在人们的日常活动和工业生产中可发挥重要作用，如在金融交易、交通监管、空气质量和温度监测、工业流程监控、医疗诊断等应用中，对大规模采样数据或高速动态数据流进行相似性查询、分类、聚类、预测、异常检测、在线模式识别等处理。

说明书

技术领域

本发明涉及数据库、数据挖掘、机器学习、信息检索等领域，尤其涉及时间序列数据分 析和挖掘。

背景技术

时间序列广泛存在于人们的日常生活及工业生产中，如基金或股票的实时交易数据，零 售市场的日销量数据，流程工业的传感器监测数据，天文观测数据，航空航天雷达、卫星监 测数据，实时天气温度及空气质量指数等。工业界迄今提出了许多时间序列分析方法，包括 相似性查询方法、分类方法、聚类方法、预测方法、异常检测方法等。其中，许多方法都需 要对时间序列进行相似性判断，比如kNN分类器、k-means聚类方法等，因此，时间序列相 似性度量方法在工业界有着广泛的应用需求。

目前工业界最常用的时间序列相似性度量方法可分为锁步度量方法和弹性度量方法。前 者采用了一对一的度量方式，即时间序列T₁和T₂之间的距离是通过严格比较T₁和T₂在各自 第i个位置的点对，再累加所有点对的距离得到。该类方法最常见的有曼哈顿距离、欧氏距 离和切比雪夫距离，它们都是L_p-norms距离在p取不同值时的特例。该类方法具有易实现、 计算复杂度低、满足距离三角不等式、无参等优点；但是，其度量精度对噪声、异常点、幅 值伸缩和漂移、相位偏移等非常敏感，并且只能用于度量等长的时间序列。弹性度量方法采 用了一对多的度量方式，即时间序列T₁的一个点可以与T₂的多个连续点相对应，通过动态规 划方法遍历T₁和T₂的所有点对之间的距离。该类方法最常见的有动态时间弯曲距离(DTW) 和编辑距离的变种(如LCSS、EDR、ERP)等。与锁步度量相比，弹性度量能够实现两条时间 序列的最佳对齐匹配，可以有效处理时间弯曲、相位偏移、幅值伸缩和漂移等基本形态变化， 对噪声和异常点具有鲁棒性，因此，弹性度量具有较高的度量精度。但是，该类方法具有较 高的计算复杂度，当度量高维的时间序列时会导致高昂的时间开销，难以在工业生产中处理 大规模的时间序列或高速的动态数据流。

基于时间序列的特征计算弹性度量是改进其高计算复杂度的一种有效方法，即首先采用 数据表示方法将原始时间序列映射到低维的特征空间，然后进行弹性度量。目前工业界常用 的数据表示方法可分为非数据适应性方法和数据适应性方法。对于前者，变换参数不受单独 的时间序列影响，而始终保持不变；该类表示大多基于频谱分解实现，如离散傅里叶变换、 离散小波变换、离散余弦变换，它们主要通过对原始时间序列做相应的频域变换，提取主要 的频谱系数作为特征；该类方法各有缺陷，如离散傅里叶变换只能提取总体形态特征而忽略 了局部特征，离散小波变换只能处理长度为2的指数次的时间序列，离散余弦变换的信息丢 失较多，对原始数据的重构误差较大。数据适应性表示是指对变换参数的确定需要依赖数据 本身；通过增加数据敏感的选择处理过程，可以把大部分非数据适应性方法变为数据适应性 方法。该类方法有分段聚集近似、分段线性近似、符号化聚集近似、奇异值分解、主成分分 析等，前三种都需要先对原始时间序列进行分段，然后对每一子段单独处理：分段聚集近似 是对各段求平均值；分段线性近似是对各段做线段拟合；符号化聚集近似是在分段聚集近似 基础上将每段平均值离散化为符号；由于它们所提取的特征较为单一，使其对时间序列波动 模式的表达能力较弱。奇异值分解和主成分分析是通过对所有时间序列做统一的特征矩阵分 解实现的；这两类方法的典型缺陷是，它们具有很高的计算复杂度，而且分解过程只能在内 存完成，数据规模的可扩展性很低。

发明内容

本发明要解决的问题是如何高效及高精度地度量时间序列之间的相似性。为了解决该问 题，本发明提出了一种基于自适应性分段统计近似的时间序列相似性度量方法。

本发明的目的是通过以下技术方案来实现的：一种基于自适应性分段统计近似的时间序 列相似性度量方法，包括以下步骤：

(1)自适应性分段，具体包括以下子步骤：

(1.1)读取原始时间序列T和Q；

(1.2)对T和Q做Z-规范化处理，得到规范化的时间序列T'和Q'；

(1.3)对规范化的时间序列T'和Q'做移动平滑处理，得到平滑时间序列T"和Q"；

(1.4)基于滑动窗口依次截取T"和Q"的相邻3点，并计算平均值，通过判断各点与相 应平均值的大小关系对其编码，得到T和Q的编码序列C_T和C_Q，并定义转折模式表TP_table；

(1.5)顺序扫描C_T和C_Q，对每对相邻编码组合查询TP_table中的转折模式，如果模式 匹配，则将该编码组合所在位置作为分段点；

(1.6)扫描完毕，分别将T和Q分为M和N段子序列，得到子序列集合S_T＝{S₁,...,S_M} 和S_Q＝{S'₁,...,S'_N}；

(2)特征提取，具体包括以下子步骤：

(2.1)依次扫描S_T'和S_Q，依次读取T和Q的每条子序列S_i和S'_i；

(2.2)依次计算S_i和S'_i的平均值μ、标准差σ、离散系数CV、偏态SK、峰态K，构造 局部模式特征向量LPV＝[μ,σ,CV,SK,K]；

(2.3)扫描完毕，得到T和Q的自适应性分段统计近似表示APSA(T)和APSA(Q)；

(3)动态模式匹配，具体包括以下子步骤：

(3.1)初始化动态规划表Table＝cell(M,N)；

(3.2)依次计算APSA(T)的第1个局部模式特征向量LPV₁与APSA(Q)的N个局部模式 特征向量LPV'₁～LPV'_N之间的加权欧氏距离{dist(LPV₁,LPV'₁),...,dist(LPV₁,LPV'_N)}，并存入 Table的第1行Table(1,1:N)；

(3.3)依次计算APSA(Q)的第1个局部模式特征向量LPV'₁与APSA(T)的M个局部模式 特征向量LPV₁～LPV_M之间的加权欧氏距离{dist(LPV₁,LPV'₁),...,dist(LPV_M,LPV'₁)}，并存入Table 的第1列Table(1:M,1)；

(3.4)利用动态规划方法，依次扫描APSA(T)的第2到第M个局部模式特征向量 LPV₂～LPV_M和APSA(Q)的第2到第N个局部模式特征向量LPV'₂～LPV'_N，基于加权欧氏距离 计算Table(2:M,2:N)的每个单元值；

(3.5)返回Table(M,N)的值作为最终的度量结果。

进一步地，所述步骤3.4包括以下子步骤：

(3.4.1)顺序扫描LPV₂～LPV_M，对于第i个局部模式特征向量LPV_i，依次计算它与 LPV'₂～LPV'_N之间的加权欧氏距离{dist(LPV_i,LPV'₂),...,dist(LPV_i,LPV'_N)}；

(3.4.2)根据先行后列的顺序扫描Table(2:M,2:N)，在每个单元Ta ble(i,j)中，首先比较 Table(i-1,j)、Table(i,j-1)和Table(i-1,j-1)的大小，选择最小值记为min，然后计算dist(LPV_i, LPV'_j)+min的值赋予Table(i,j)。

本发明的有益效果是：

1、在自适应性分段阶段，采用了简单有效的编码方法和转折模式识别方法，可高效识别 转折点，保证了切分出的子序列具有完整的波动趋势。

2、在特征提取阶段，对每条子序列提取多种统计特征，从多方面反映了时间序列的波动 特性，可全面捕捉时间序列的局部波动模式，实现了较高的时间序列局部模式匹配精度。

3、在动态模式匹配阶段，基于局部模式层次的动态规划计算，克服了时间弯曲造成的局 部模式之间的相位偏移问题，实现了较高的时间序列全局模式匹配精度。

附图说明

图1为基于自适应性分段统计近似的时间序列相似性度量方法流程图；

图2为自适应性分段时间序列的流程图；

图3为采用自适应性分段统计近似表示时间序列的流程图；

图4为时间序列相似性计算的动态模式匹配过程。

具体实施方式

下面结合附图对本发明作进一步详细说明。

如图1所示，本发明一种基于自适应性分段统计近似的时间序列相似性度量方法，包括 以下步骤：

(1)自适应性分段，如图2所示，具体包括以下子步骤：

(1.1)读取原始时间序列T＝{t₁,t₂,…,t_i,…,t_n}和Q＝{q₁,q₂,…,q_i,…,q_n}；

(1.2)对时间序列T和Q，分别计算T的采样点的平均值m'和标准差σ'，Q的采样点的 平均值m'和标准差σ'，根据公式(1)对T和Q做Z-规范化处理，得到规范化的时间序列T'＝{t'₁, t'₂,…,t'_i,…,t'_n}和Q'＝{q'₁,q'₂,…,q'_i,…,q'_n}；

      ${t^{\prime }}_i=\frac{t_i-m}{\sigma }---\left(1\right)$      

(1.3)依次计算T'和Q'相邻3点的平均值，对其做移动平滑处理，得到平滑时间序列 T"＝{t"₁,t"₂,…,t"_i,…,t"_n}和Q"＝{q"₁,q"₂,…,q"_i,…,q"_n}；

(1.4)基于滑动窗口依次截取T"和Q"的相邻3点，并计算平均值，通过判断各点与相 应平均值的大小关系对其编码，得到T和Q的编码序列C_T和C_Q，并定义转折模式表TP_table， 该过程包括以下子步骤：

(1.4.1)采用滑动窗口W，依次截取T"和Q"的相邻3点<t"_i-1,t"_i,t"_i+1>和<q"_i-1,q"_i,q"_i+1>， 并计算平均值m^t_i和m^q_i；

(1.4.2)判断<t"_i-1,t"_i,t"_i+1>和<q"_i-1,q"_i,q"_i+1>的各点与相应平均值m^t_i和m^q_i的关系，若 t"_i>m^t_i，则code(t"_i)＝1；否则code(t"_i)＝0，由此将<t"_i-1,t"_i,t"_i+1>和<q"_i-1,q"_i,q"_i+1>编码为d^t_i＝<c^t_i-1,c^t_i,c^t_i+1>和d^q_i＝<c^q_i-1,c^q_i,c^q_i+1>；由此得到T和Q的编码序列C_T＝{d^t₁,d^t₂,...,d^t_n}和C_Q＝{d^q₁,d^q₂,...,d^q_n}；

(1.4.3)根据编码定义所有转折模式TP，得到转折模式表TP_table＝{上升-下降：001-100, 001-110,011-100,011-110,001/011-010-100/110；下降-上升：100-001,100-011,110-001,110-011, 100/110-101-001/011}；

(1.5)顺序扫描C_T和C_Q，对每对相邻编码组合<d^t_i,d^t_i+1>和<d^q_i,d^q_i+1>查询TP_table，如 果模式匹配，则将该将i作为分段点，得到T和Q的第i条子序列S_i和S'_i；

(1.6)扫描完毕，对T和Q完成分段，得到子序列集合S_T＝{S₁,S₂,...,S_M}和S_Q＝{S'₁,S'₂,..., S'_N}；

(2)特征提取，如图3所示，具体包括以下子步骤：

(2.1)依次扫描S_T'和S_Q，依次读取T和Q的每条子序列S_i和S'_i；

(2.2)依次对S_i和S'_i计算多种统计特征，构造局部模式特征向量LPV_i和LPV'_i，该过程 包括以下子步骤：

(2.2.1)初始化T和Q的分段统计近似表示APSA(T)和APSA(Q)为空集；

(2.2.2)根据公式(2)，计算长度为l的子序列S_i和S'_i的平均值μ_i和μ'_i；

      $\mu =\left(\underset{i}{\overset{i+l-1}{\Sigma }}{s}_i\right)/l---\left(2\right)$      

(2.2.3)依次根据公式(3)～(7)，计算S_i和S'_i的方差D、标准差σ、离散系数CV、偏态 SK、峰态K，分别构造局部模式特征向量LPV_i＝[μ_i,D_i,σ_i,CV_i,SK_i,K_i]和LPV'_i＝[μ'_i,D'_i,σ'_i,CV'_i, SK'_i,K'_i]，并分别插入APSA(T)和APSA(Q)；

      $D=\underset{i}{\overset{i+l-1}{\Sigma }}{(s_i-\mu )}^2/(l-1)---\left(3\right)$      

      $\sigma =\sqrt{\underset{i}{\overset{i+l-1}{\Sigma }}{(s_i-\mu )}^2/(l-1)}---\left(4\right)$      

      $\mathrm{CV}=\frac{\sigma }{\mu }---\left(5\right)$      

      $\mathrm{SK}=\frac{l\underset{i=1}{\overset{i+l-1}{\Sigma }}{(s_i-\mu )}^3}{(l-1)(l-2){\sigma }^3}---\left(6\right)$      

      $K=\frac{l(l+1)\underset{i=1}{\overset{i+l-1}{\Sigma }}{(s_i-\mu )}^4-3(l-1){\left[\underset{i=1}{\overset{i+l-1}{\Sigma }}{(s_i-\mu )}^2\right]}^2}{(l-1)(l-2)(l-3){\sigma }^4}---\left(7\right)$      

(2.3)扫描完毕，得到T和Q的自适应性分段统计近似表示APSA(T)和APSA(Q)；

(3)动态模式匹配，如图4所示，具体包括以下子步骤：

(3.1)初始化动态规划表Table＝cell(M,N)；

(3.2)根据公式(8)，依次计算APSA(T)的第1个局部模式特征向量LPV₁与APSA(Q)的N 个局部模式特征向量LPV'₁～LPV'_N之间的加权欧氏距离{dist(LPV₁,LPV'₁),...,dist(LPV₁,LPV'_N)}， 并依次存入Table的第1行Table(1,1:N)；

      $\mathrm{dist}\{L{\mathrm{PV}}_i,{\mathrm{LPV}}_j^{\prime }\}=\sqrt{\underset{k=1}{\overset{6}{\Sigma }}{a}_k·{(v_k-{v^{\prime }}_k)}^2}---\left(8\right)$      

其中，a_k表示局部模式特征向量第k个特征的权重系数，v_k和v'_k分别表示LPV和LPV' 的第k个元素。

(3.3)根据公式(8)，依次计算APSA(Q)的第1个局部模式特征向量LPV'₁与APSA(T)的 M个局部模式特征向量LPV₁～LPV_M之间的加权欧氏距离 {dist(LPV₁,LPV'₁),...,dist(LPV_M,LPV'₁)}，并依次存入Table的第1列Table(1:M,1)；

(3.4)利用动态规划方法，基于公式(8)计算Table(2:M,2:N)的每个单元值，该过程包括 以下子步骤：

(3.4.1)顺序扫描LPV₂～LPV_M，对于APSA(T)的第i个局部模式特征向量LPV_i，依次计 算它与LPV'₂～LPV'_N之间的加权欧氏距离{dist(LPV_i,LPV'₂),...,dist(LPV_i,LPV'_N)}；

(3.4.2)当扫描LPV_i与LPV'_j时，首先比较Table(i-1,j)、Table(i,j-1)和Table(i-1,j-1)的大 小，选择最小值记为min，然后计算dist(LPV_i,LPV'_j)+min的值赋予Table(i,j)。

(3.5)返回Table(M,N)的值作为最终的度量结果。

时间序列相似性度量，在人们的日常活动及工业生产中可发挥重要作用，有着广泛的应 用需求。本发明针对工业界当前提出的众多时间序列分析方法，提出了一种基于自适应性分 段统计近似表示的时间序列相似性度量方法，可以对时间序列进行数据适应性地分段，并实 现高效及高精度地相似性度量，由此实现对大规模采样数据或高速动态数据流进行相似性查 询、分类、聚类、预测、异常检测、在线模式识别等处理，以满足工业生产的应用需求。