基于小波和动态时间弯曲的时间序列相似匹配

D01:10.13374j.isml00103x.206.04.21 第28卷第4期 北京科技大学学报 Vol.28 Na 4 2006年4月 Journal of University of Science and Technology Beijing Apr.2006 基于小波和动态时间弯曲的时间序列相似匹配 曲文龙张德政杨炳儒 北京科技大学信息工程学院，北京100083 摘要提出了一种基于小波和动态时间弯曲(DTW距离的时间序列索引和相似匹配方法.该 方法采用小波变换进行数据降维，利用R一re建立多维索引结构.给出了查询序列的DTW距 离边界和其在小波空间的查询超矩形的计算方法.从而将原始空间的基于DTW距离的相似匹配 转换为小波空间基于欧氏距离的相似匹配.证明了此匹配方法不会产生漏报，给出了基于DTW 距离的范围查询算法和近邻查询算法.实验结果表明该方法具有较高匹配精度和其较低的计算代 价. 关键词时间序列：相似匹配：动态时间弯曲：小波变换 分类号TP311.13 时间序列分析与挖掘广泛应用于包括金融、 出了范围查询和近邻查询算法，最后通过对比试 商业、气象、医学、电力、水文、工业等众多领域，具 验验证了该方法的优越性. 有重要的研究价值。相似模式匹配是时间序列分 1动态时间弯曲距离 类、聚类、规则获取和预测等挖掘方法的基础，建 立索引是实现时间序列相似匹配的有效方法. 11动态时间弯曲概念 典型的相似性测度多采用欧几里德距离或其 设有两个时间序列Q=(q1,；q,9n) 改进，但欧氏距离测度存在局限性，对数据在时间 和C=(c1,,G,,cm),长度分别为n和m(如 轴上的形变缺乏辨识能力和对噪声的鲁棒性.动 图1)，为利用DTW将两个时间序列对准，首先定 态时间弯曲(Dy namic Time Warping,DTW)川可 义DTW对准矩阵M. 以获得很高的识别、匹配精度.由于DTW距离 不满足三角不等式，在低维特征空间中无法保证 (a) 检索的完整性，因此基于DTW的索引和相似性 搜索有待研究.Yi给出了一个基于FASTMAP 映射的DTW索引方法？，但只是近似索引，无法 保证检索完整性.Pak采用线性分段表示建立 DTW索引3，但无法保证无漏报且查询精度较 (c) 低.Km采用序列的四个特征建立索引，保证 图1(a)时间序列O和C:(b)DTW对准方式：(cDTW对 无漏报，但未能有效的利用多维索引结构缩减搜 准矩阵和最佳弯曲路径 索空间，误报率很高.Keogh等给出了局部DTW Fig.I (a)Time serise O and C;(b)DTW aligned:(c)DTW 的包围边界，并采用分段近似表示匀，保证无漏 aligned matrix and optimum wraping path 报并使搜索空间得到一定程度缩减.由于小波变 定义1n行m列矩阵M,矩阵中的元素 换可以有效约简特征空间、具有多尺度特性，因此 (i,j)为两时间序列数据中对准点9：和G之间的 本文采用小波方法进行时间序列的DTW距离索 引和相似匹配，给出了DTW的小波变换边界，进 距离d(q,c(d(q,G)=(q-c92),定义M为 时间序列Q和C的DTW对准矩阵. 一步缩减搜索空间，并证明了该方法无漏报且给 定义2对于两个时间序列的DTW对准异 收稿日期：2005-0221修回日期：200504-24 矩阵定义矩阵中一组连续的矩阵元素的集合W 基金项目：北京市自然科学基金资助项目(No.4022008):国家科 ={w1,wk,;wK,(wk=d(9,G),称满 技攻关项目(No.2004BA6164一11) 作者简介：曲文龙(1970一），男.博士研究生 足如以下条件的W为时间序列Q和C的弯曲路

基于小波和动态时间弯曲的时间序列相似匹配 曲文龙 张德政 杨炳儒 北京科技大学信息工程学院, 北京 100083 摘 要 提出了一种基于小波和动态时间弯曲( DTW) 距离的时间序列索引和相似匹配方法.该 方法采用小波变换进行数据降维, 利用 R ＊-tree 建立多维索引结构.给出了查询序列的 DTW 距 离边界和其在小波空间的查询超矩形的计算方法, 从而将原始空间的基于 DTW 距离的相似匹配 转换为小波空间基于欧氏距离的相似匹配.证明了此匹配方法不会产生漏报, 给出了基于 DTW 距离的范围查询算法和近邻查询算法.实验结果表明该方法具有较高匹配精度和其较低的计算代 价. 关键词 时间序列;相似匹配;动态时间弯曲;小波变换 分类号 TP311.13 收稿日期:2005 02 21 修回日期:2005 04 24 基金项目:北京市自然科学基金资助项目( No .4022008) ;国家科 技攻关项目( No .2004BA616A-11) 作者简介:曲文龙( 1970—) , 男, 博士研究生 时间序列分析与挖掘广泛应用于包括金融、 商业 、气象 、医学 、电力、水文、工业等众多领域, 具 有重要的研究价值.相似模式匹配是时间序列分 类、聚类 、规则获取和预测等挖掘方法的基础, 建 立索引是实现时间序列相似匹配的有效方法 . 典型的相似性测度多采用欧几里德距离或其 改进, 但欧氏距离测度存在局限性, 对数据在时间 轴上的形变缺乏辨识能力和对噪声的鲁棒性.动 态时间弯曲( Dy namic Time Warping , DTW) [ 1] 可 以获得很高的识别 、匹配精度 .由于 DTW 距离 不满足三角不等式, 在低维特征空间中无法保证 检索的完整性, 因此基于 DTW 的索引和相似性 搜索有待研究.Yi 给出了一个基于 FASTMAP 映射的 DTW 索引方法[ 2] , 但只是近似索引, 无法 保证检索完整性.Park 采用线性分段表示建立 DTW 索引 [ 3] , 但无法保证无漏报且查询精度较 低.Kim 采用序列的四个特征建立索引[ 4] , 保证 无漏报, 但未能有效的利用多维索引结构缩减搜 索空间, 误报率很高 .Keog h 等给出了局部 DTW 的包围边界, 并采用分段近似表示 [ 5] , 保证无漏 报并使搜索空间得到一定程度缩减.由于小波变 换可以有效约简特征空间 、具有多尺度特性, 因此 本文采用小波方法进行时间序列的 DTW 距离索 引和相似匹配, 给出了 DTW 的小波变换边界, 进 一步缩减搜索空间, 并证明了该方法无漏报, 且给 出了范围查询和近邻查询算法, 最后通过对比试 验验证了该方法的优越性 . 1 动态时间弯曲距离 1.1 动态时间弯曲概念 设有两个时间序列 Q =( q1, …, qi , …, qn ) 和 C =( c1, …, cj , …, cm ), 长度分别为 n 和m(如 图 1), 为利用 DTW 将两个时间序列对准, 首先定 义 DTW 对准矩阵 M . 图1 ( a) 时间序列 Q 和C ;( b) DTW 对准方式;( c) DTW对 准矩阵和最佳弯曲路径 Fig.1 ( a) Time serise Q and C ;( b) DTW aligned;( c) DTW aligned matrix and optimum wraping path 定义 1 n 行 m 列矩阵 M, 矩阵中的元素 ( i, j) 为两时间序列数据中对准点 qi 和cj 之间的 距离d ( qi , cj)( d ( qi , cj) =( qi -cj) 2 ), 定义 M 为 时间序列Q 和C 的 DTW 对准矩阵 . 定义 2 对于两个时间序列的 DTW 对准异 矩阵, 定义矩阵中一组连续的矩阵元素的集合 W ={w 1, …, wk , …, wK}, ( wk =d ( qi , cj )), 称满 足如以下条件的 W 为时间序列Q 和C 的弯曲路 第 28 卷 第 4 期 2006 年 4 月 北 京 科 技 大 学 学 报 Journal of University of Science and Technology Beijing Vol .28 No.4 Apr.2006 DOI :10.13374/j .issn1001 -053x.2006.04.021

Vol.28 No.4 曲文龙等：基于小波和动态时间弯曲的时间序列相似匹配 ·397。 径 定义5对任意两个时间序列Q和C,U和 (1)有界性：max(m,n)≤K≤m十n十1： L为Q的DTW上边界序列和DTW下边界序 (2边界性：w1=d(q1,c1,wk=d(9a, 列，定义Q和C的低限距离为： Cm); (c-u2, G≥ui (3)连续性：给定wk=d(qa,c6),w-1=d diB(O,C) (h,-c)2, G≤1， (qa°，cb),必有a-a*≤1且b-b*≤1： 0, li<ci<ui (4)单调性：给定wk=d(qa,Cb),w-1=d (2) (ga°，cb,必有a-a≥0且b-b≥0. 定理1对任意两个时间序列Q和C,r为 除上述限制条件外，在实际应用中还需要限 允许弯曲范围，则其dB距离是两序列DTW距离 制弯曲路径的宽度，防止病态弯曲和提高DTW 算法的速度，通常对弯曲窗口加以限制，即对于路 的低限距离，即dLB(Q,C)≤D6Tw(Q,C)(定理 1证明参见文献[可) 径中任一点(wk=d(q,G),要求|i-j川≤r. 定理2任给时间序列Q和C,U和L为Q 弯曲路径存在多解，DTW距离取弯曲路径总长 的DTW上边界序列和DTW下边界序列，必存在 度的最小值，即dorw(2,C)=min L≤T≤U(即l≤t:≤山：)使得dB(Q,C)= d(T,C)(d(*,*)表示两时间序列的欧距 最佳路径可以由时间起始点(1,1)到终点(m,n) 离，下同). 之间的局部最优解通过递归获得，公式如下： 证明：设U(ui=max(q-n,…,q什r)和L S(1,1)=d(q1,c1) (l=min(q-r,q+r)为Q的DTW上、下边 S(i,j)=d(gi,g)+min(S(i-1,j)(1) 界序列，构造序列T: S(i,-1)S(i-1,j-1)) ui,ciui 式中，S(i,)为累积距离，由当前对准点的距离 和相邻点的累积DTW距离计算得到，则 LCi, lC斯 dmw(Q,C)=JS(n,m).欧几里德距离可视为 显然lh≤t≤i成立，容易验证dB(Q,C)= DTW距离的特例，即第k路径为wk=d(qk,G, d(T,C). 欧氏距离要求两个序列长度相等. 时间序列DTW相似匹配的基本思想是使 12动态时间弯曲距离低限边界 用低限距离函数剪除那些不可能匹配的序列以缩 在基于DTW的相似匹配中，为防止路径病 减搜索空间，并保证不产生漏报，然后对候选序列 态弯曲和提高计算速度，通常需对弯曲路径进行 再采用DTW距离进一步精确匹配.采用的低限 限制. 距离要比使用原距离计算复杂性低。其次要使用 定义3设wk=d(i,j)为DTW的第k个 紧低限距离（定义的低限距离尽可能接近原距 路径，设弯曲路径限制为i一r≤≤i十r,r称为 离)，以减少使用低限距离检索结果中的候选数据 弯曲范围(r为常数或为i和j的函数). 的数量 定义4=(q1,…；9,；9n)为待匹配的 小波变换具有保矩性，采用小波变换将时间 时间序列，r为序列中每个点允许的最大弯曲范 序列变换到时频域，利用低频段的小波系数建索 围，构造如下两个新序列U(4=mr(q-r,, 引结构，可以提高检索效率. 9+）和L(l=min(qi-r,…,9+r)(图2)，称 U和L分别为序列Q的DTW上边界序列和 2离散小波变换 DTW下边界序列. 21基本概念 小波变换是一种非平稳信号分析方法9，它 通过一个满足条件J。9(x)dx=0的基本小波函 数P(x)的平移和伸缩构成一族小波函数系去表 图2Q的，弯曲上边界序列和下边界序列 示或逼近一个函数.二进制小波是由伸缩因子和 Fig.2 warping upper and lower bound sequences of O 平移因子满足一定条件的一组函数：

径: ( 1) 有界性:max( m, n) ≤K ≤m +n +1 ; ( 2) 边界性 :w1 =d ( q1, c1 ), wK =d ( qn, cm) ; ( 3) 连续性:给定 wk =d ( qa, cb) , wk-1 =d ( qa·, cb·), 必有 a -a ＊ ≤1 且 b -b ＊≤1 ; ( 4) 单调性:给定 wk =d ( qa, cb) , wk-1 =d ( qa·, cb·), 必有 a -a ＊ ≥0 且 b -b ＊≥0 . 除上述限制条件外, 在实际应用中还需要限 制弯曲路径的宽度, 防止病态弯曲和提高 DTW 算法的速度, 通常对弯曲窗口加以限制, 即对于路 径中任一点( wk =d ( qi , cj)) , 要求 i -j ≤r . 弯曲路径存在多解, DTW 距离取弯曲路径总长 度的最小值, 即 d DTW( Q, C) =min ∑ K k =1 wk , 最佳路径可以由时间起始点( 1, 1) 到终点( m, n) 之间的局部最优解通过递归获得, 公式如下: S ( 1, 1) =d( q1, c1) S ( i, j) =d( qi , cj) +min( S ( i -1, j) S ( i, j -1), S ( i -1, j -1) ) ( 1) 式中, S ( i, j) 为累积距离, 由当前对准点的距离 和相 邻 点 的 累 积 DTW 距 离 计 算 得 到, 则 dDTW( Q, C) = S ( n, m) .欧几里德距离可视为 DTW 距离的特例, 即第 k 路径为wk =d( qk , ck ), 欧氏距离要求两个序列长度相等. 1.2 动态时间弯曲距离低限边界 在基于 DTW 的相似匹配中, 为防止路径病 态弯曲和提高计算速度, 通常需对弯曲路径进行 限制 . 定义 3 设 wk =d ( i, j) 为 DTW 的第 k 个 路径, 设弯曲路径限制为 i -r ≤j ≤i +r, r 称为 弯曲范围( r 为常数或为i 和j 的函数) . 图 2 Q 的r 弯曲上边界序列和下边界序列 Fig.2 r-warping upper and lower bound sequences of Q 定义 4 Q =( q1, …, qi , …, qn ) 为待匹配的 时间序列, r 为序列中每个点允许的最大弯曲范 围, 构造如下两个新序列 U ( ui =max ( qi-r, …, qi+r))和 L( li =min( qi -r, …, qi+r)) ( 图 2), 称 U 和 L 分别为序列 Q 的 DTW 上边界序列和 DTW 下边界序列 . 定义 5 对任意两个时间序列 Q 和C, U 和 L 为 Q 的 DTW 上边界序列和 DTW 下边界序 列, 定义 Q 和C 的低限距离为: d LB( Q, C) = ∑ n i =1 ( ci -ui) 2 , ci ≥ui ( li -ci) 2 , ci ≤li 0, li <ci < ui ( 2) 定理 1 对任意两个时间序列 Q 和C, r 为 允许弯曲范围, 则其 d LB距离是两序列DTW 距离 的低限距离, 即 d LB( Q, C) ≤D r DTW ( Q, C) (定理 1 证明参见文献[ 5] ) . 定理 2 任给时间序列 Q 和C, U 和L 为Q 的DTW 上边界序列和 DTW 下边界序列, 必存在 L ≤T ≤U ( 即 li ≤ti ≤ui ) 使得 dLB ( Q, C) = d( T, C)( d ( ＊, ＊) 表示两时间序列的欧氏距 离, 下同) . 证明:设 U ( ui =max ( qi -r, …, qi+r)) 和 L ( li =min( qi-r, …, qi+r))为 Q 的 DTW 上、下边 界序列, 构造序列 T : ti = ui , ci >ui li , ci <li ci , li ≤ci ≤ui 显然 li ≤ti ≤ui 成立, 容易验证 d LB( Q, C) = d( T, C) . 时间序列 DTW 相似匹配的基本思想, 是使 用低限距离函数剪除那些不可能匹配的序列以缩 减搜索空间, 并保证不产生漏报, 然后对候选序列 再采用 DTW 距离进一步精确匹配 .采用的低限 距离要比使用原距离计算复杂性低, 其次要使用 紧低限距离( 定义的低限距离尽可能接近原距 离), 以减少使用低限距离检索结果中的候选数据 的数量. 小波变换具有保矩性, 采用小波变换将时间 序列变换到时频域, 利用低频段的小波系数建索 引结构, 可以提高检索效率. 2 离散小波变换 2.1 基本概念 小波变换是一种非平稳信号分析方法[ 4] , 它 通过一个满足条件∫R φ( x ) d x =0 的基本小波函 数 φ( x )的平移和伸缩构成一族小波函数系去表 示或逼近一个函数.二进制小波是由伸缩因子和 平移因子满足一定条件的一组函数 : Vol.28 No.4 曲文龙等:基于小波和动态时间弯曲的时间序列相似匹配 · 397 ·

。398· 北京科技大学学报 2006年第4期 9.k(x)=229(2x-k,j,k∈Z(3) 索引并进行相似匹配，需给出两时间序列Q和C 对任意平方可积函数∫(x)来说，其离散小波变 的DTW距离在小波域的低限距离定义d(O, 换为： C),并证明d"s(Q,C)为其DTW距离的低限距 十 离，从而保证相似匹配无漏报.Haar小波二进制 Wgfj,k材= =f x)9(x)dx=(f.9),f(x) 平移系为： (4) 2?2≤K(2k+1)21 其中，j,k∈Zj为尺度因子，k为平移因子.信 hj.k(t)= -2"?,(2k+)2-≤K(k+1)2 号可以用各尺度的小波系数重构： 0. K2k或⊙≥(k+1)2 fx)=∑胸，k)9(x) (5) (8) i.k 给定时间序列X={x(t)}(i=L,,n)(设 其中，j,k∈z,j为尺度因子，k为平移因子.对 n=2,不足位补0)，利用Har小波变换对X进 于离散时间序列Q=(q1,,qi,；qm人，其Haar 行J尺度分解，分解系数为山，D,;Di,;D1 的小波系数为： (2+)2广1 (J为最大分解尺度).A山为原始序列在尺度J下 Wo(j,k)= 的逼近信号，D为原始序列在尺度1≤≤J)下 =2k 的细节信号.j尺度信号的长度为2一（即k=0, (9) -1,小波分解系数的总长度为 =H2+)21 +召-%与原序列张度相等. 定义7设时间序列Q的DTW上下包围边 界序列为U和L,按如下方法构造两个序列U" 小波分解系数的排列次序为：J尺度近似信 和L": 号、按尺度j从大到小的各尺度细节信号（对于同 . 一尺度的小波系数按平移因子k从小到大排 =1+24 列).由于在时间序列相似匹配中一般首先对序 (10 列进行归一化，J尺度近似值A均为0，因此只取 各尺度细节信号即可. L(,k)= =+2k =1H2+1)2 定义6对长为n(n=2)的序列X进行J (11) 层小波分解，小波系数按尺度值从大到小排列 式中，=1，；lbg2n,k=0,；j广-1. (同一尺度的小波系数按平移因子k从小到大排 将U(j,k)按j从大到小（相同时按k从小 列)，取前m(1≤m≤n一1)个系数可得到序列 到大)排列，得到序列U"={};将L(j,k)按j X"={x}(i=l,,m),定义为X的小波变换 从大到小(j相同时按k从小到大)排列，取前m 序列. 个系数得到序列L"={}(=1,…,m).则称 利用小波变换序列近似表示原序列，构建索 U"和L"为序列Q在小波变换域的DTW上下 引结构，使索引维数得到约简，对得到的匹配候选 边界序列. 序列.可通过后处理（计算DTW距离）来滤除匹 定理3对任意两个时间序列Q和X,U和 配误报.研究中采用简单易用的Har小波，也可 L为Q的DTW上边界序列和下边界序列，U" 以采用其他正交小波，如Daubechies紧支集小波 和L"为序列Q在小波变换域的DTW上、下边 等，Har小波是正交小波，其小波函数P(x)为： 界序列，X"为X的小波变换序列.如果L≤X 1,0x≤0.5 ≤U(≤x≤)，必有L"≤x"≤U". 9(x)= -1, 0.5<x≤1 (6) 证明：设X的小波系数为Wx(,k),由于 0. x≤0或x>1 l≤x≤u,所以有 尺度函数为： (2+1)2 wX(j.)= x)= L,0x≤1 7) =2 =H21)2 0,x≤0或x心1 2.2小波变换域的动态弯曲低限距离 利用时间序列的前m个小波变换系数建立 =1+H2+1)Y

φj , k ( x ) =2 -j/2 φ( 2 -j x -k ), j, k ∈ Z ( 3) 对任意平方可积函数 f ( x ) 来说, 其离散小波变 换为 : Wφf( j, k) =∫ +∞ -∞ =f( x) φj , k ( x) d x =〈f , φj , k〉, f( x) ( 4) 其中, j, k ∈ Z, j 为尺度因子, k 为平移因子 .信 号可以用各尺度的小波系数重构: f ( x ) = ∑ j, k Wφf ( j, k ) φj, k ( x ) ( 5) 给定时间序列 X ={x ( ti)}( i =1, …, n) (设 n =2 J , 不足位补 0) , 利用 Haar 小波变换对 X 进 行J 尺度分解, 分解系数为 AJ , DJ , …, Dj , …, D1 ( J 为最大分解尺度) .AJ 为原始序列在尺度J 下 的逼近信号, Dj 为原始序列在尺度j( 1 ≤j ≤J )下 的细节信号.j 尺度信号的长度为2 J-j (即 k =0, …, 2 J -j -1 ) , 小 波 分 解 系 数 的 总 长 度 为 1 + ∑ J j =1 2 J -j =2 J =n, 与原序列长度相等. 小波分解系数的排列次序为 :J 尺度近似信 号、按尺度 j 从大到小的各尺度细节信号(对于同 一尺度的小波系数按平移因子 k 从小到大排 列) .由于在时间序列相似匹配中一般首先对序 列进行归一化, J 尺度近似值 AJ 均为0, 因此只取 各尺度细节信号即可 . 定义 6 对长为 n ( n =2 J ) 的序列 X 进行J 层小波分解, 小波系数按尺度 j 值从大到小排列 (同一尺度的小波系数按平移因子 k 从小到大排 列), 取前 m ( 1 ≤m ≤n -1) 个系数可得到序列 X w ={x w i }( i =1, …, m), 定义为 X 的小波变换 序列 . 利用小波变换序列近似表示原序列, 构建索 引结构, 使索引维数得到约简, 对得到的匹配候选 序列.可通过后处理(计算 DTW 距离) 来滤除匹 配误报 .研究中采用简单易用的 Haar 小波, 也可 以采用其他正交小波, 如 Daubechies 紧支集小波 等, Haar 小波是正交小波, 其小波函数 φ( x )为: φ( x ) = 1, 0 <x ≤0.5 -1, 0.5 <x ≤1 0, x ≤0 或 x >1 ( 6) 尺度函数为: ψ( x ) = 1, 0 <x ≤1 0, x ≤0 或 x >1 ( 7) 2.2 小波变换域的动态弯曲低限距离 利用时间序列的前 m 个小波变换系数建立 索引并进行相似匹配, 需给出两时间序列 Q 和C 的 DTW 距离在小波域的低限距离定义 d w LB( Q, C), 并证明 d w LB( Q, C) 为其 DTW 距离的低限距 离, 从而保证相似匹配无漏报 .Haar 小波二进制 平移系为 : hj, k ( t) = 2 -j/ 2 , 2 j k ≤t <( 2k +1) 2 j-1 -2 -j/ 2 , ( 2k +1) 2 j -1 ≤t <( k +1)2 j 0, t <2 j k 或 t ≥( k +1)2 j ( 8) 其中, j, k ∈ z, j 为尺度因子, k 为平移因子.对 于离散时间序列 Q =( q1, …, qi , …, qn ), 其 Haar 的小波系数为 : WQ( j, k ) = ∑ ( 2 k+1) 2 j-1 i =1+2 j k q( i) 2 -j 2 - ∑ ( k+1) 2 j i =1+( 2 k+1) 2 j-1 q( i) 2 -j 2 ( 9) 定义 7 设时间序列 Q 的 DTW 上下包围边 界序列为 U 和L, 按如下方法构造两个序列 U w 和 L w : U( j , k ) = ∑ ( 2k +1) 2 j -1 i =1 +2 j k u( i) 2 - j 2 - ∑ ( k+1) 2 j i =1+( 2 k+1) 2 j-1 l( i) 2 - j 2 ( 10) L ( j , k ) = ∑ ( 2 k+1) 2 j-1 i =1+2 j k l( i) 2 - j 2 - ∑ ( k+1) 2 j i =1 +( 2k +1) 2 j-1 u( i) 2 - j 2 ( 11) 式中, j =1, …, log2 n , k =0, …, j -1 . 将 U( j, k )按 j 从大到小( j 相同时按k 从小 到大)排列, 得到序列 U w ={u w i };将 L( j, k )按 j 从大到小( j 相同时按 k 从小到大) 排列, 取前 m 个系数得到序列 L w ={l w i }( i =1, …, m ) .则称 U w 和 L w 为序列 Q 在小波变换域的 DTW 上下 边界序列 . 定理 3 对任意两个时间序列 Q 和 X , U 和 L 为 Q 的 DTW 上边界序列和下边界序列, U w 和 L w 为序列 Q 在小波变换域的 DTW 上 、下边 界序列, X w 为 X 的小波变换序列.如果 L ≤X ≤U( li ≤x i ≤ui) , 必有 L w ≤X W ≤U w . 证明 :设 X 的小波系数为 WX ( j, k ), 由于 li ≤xi ≤ui , 所以有 WX ( j, k ) = ∑ ( 2 k+1) 2 j-1 i =1+2 j k x( i) 2 -j 2 - ∑ ( k+1) 2 j i =1+( 2 k+1) 2 j-1 x( i ) 2 -j 2 ≤ ∑ ( 2 k+1) 2 j-1 i =1+2 j k u( i) 2 - j 2 - ∑ ( k +1) 2 j i =1 +( 2k +1) 2 j -1 l( i) 2 - j 2 =U( j , k) · 398 · 北 京 科 技 大 学 学 报 2006 年第 4 期

Vol.28 No.4 曲文龙等：基于小波和动态时间弯曲的时间序列相似匹配 ·399。 和 由于Haar为正交小波满足parseval定理，可 得 d(T",c")≤d(T,C). =什2 H2k+1)2 所以， - d"sQ,C)≤diB(Q,C). =H2+1)2广1 所以L(j,k)≤Wx(j,k)≤U(j,k), 又由定理L,可得 即le≤xae≤ue, dsgC)≤Dbw(Q,C), 可证L“≤X"≤U" 所以， 定义8对任意时间序列Q和C,C"为C diB(g,C)≤DDTw(Q,C). 的小波变换序列，U"和L"为序列Q在小波域 定理6采用小波域的dB距离进行DTW 的DTW上、下边界序列.如下定义Q和C的 距离相似匹配不会产生漏报. d距离： 证明：设相似阈值为飞，Q为待查询序列，任 给时间序列C,其弯曲范围为r的DTW距离为 (c-u)2, c≥u D'orw(O,C). diB(2,C)= (1-c)2, c"≤I 如果Dbw(Q,C)<e,由定理5，必有diB 0, 俏<c"<u (Q,C)<e,即使用ds距离在小波域进行相似 (12) 匹配得到的候选序列必包含所有Dw距离的匹 定理4任给时间序列Q,C和X,U和L 配序列，不会产生漏报（即在小波域没有检索到却 为Q的DTW上边界序列和下边界序列，C"为C 满足查询条件的序列) 的小波变换序列，X“为X的小波变换序列.如 定理6保证对于采用DTW距离的相似匹 果L≤X≤U(lh≤x≤)，则： 配可以使用小波变换域定义的d"距离检索，不 disg,C)≤d(X",C". 会产生漏报，但可能有误报（即在小波检索到却不 证明： 满足查询条件的序列)，误报可以通后处理来过滤 以得到到正确的检索结果。该方法定义的低限距 d(C",X")= (ci-x2= = 离比先前方法相比更紧，可最大限度的减小搜索 x)只， 空间.采用小波系数方法建立索引，利用其降维 (c"- c≥u 能力和多尺度特性提高搜索效率. (c-x)2,c≤片 (ei-xi),li<ei<ui 3基于小波和DTW距离的相似匹 因为L≤X≤U,由定理3可得L"≤X"≤ 配算法 U",即≤x"≤，所以当c≥u时，有(c”- 3.1建立索引 x)2≥c曾-u)2.当c≤u时，有(c-x)2 首先对长度为n(设n=2,不足位补0)时间 ≥(c-u)2.当<c<:时，显然有(c- 序列X实施归一化，以校正序列在Y轴的平移 x)2≥0.对照式(12)，可得d"(C,2)≤d 和幅度伸缩.即，=二x,其中x= (x",C". xmax一xmin 定理5任给时间序列Q和C,必有d" 另，xm和xn为序列的最大值和最小值。 (2,C)≤Dbw(2,C). 然后对标准化后的序列实施离散小波变换，取前 证明：设U和L为Q的DTW上边界序列和 m个小波系数组成小波系数序列，视为m维空 下边界序列，由定理2，必存在L≤T≤U(即l≤ 间的一个点，并将降维之后的小波系数序列组织 t≤u),使得dsQ,C)=d(T,C). 为多维索引结构R'-tree索引. 设C"为C的小波变换序列，T“为T的小 R一tree是一种多级平衡树9，树中的每个 波变换序列，由定理4可得： 非叶节点对应一个多维超矩形，该超矩形为其子 diBO,C)≤d(T",C"). 节点代表超矩形的最小包围超矩形(MBR).非叶

和 WX ( j, k ) = ∑ ( 2 k+1) 2 j-1 i =1+2 j k x( i) 2 -j 2 - ∑ ( k+1 ) 2 j i =1+( 2 k+1) 2 j-1 x( i) 2 -j 2 ≥ ∑ ( 2 k+1) 2 j-1 i =1+2 j k l( i) 2 - j 2 - ∑ ( k+1) 2 j i =1+( 2 k+1) 2 j-1 u( i) 2 - j 2 =L( j , k) , 所以 L ( j, k ) ≤WX( j, k) ≤U( j, k) , 即 l wave l ≤x wave i ≤u wave i , 可证 L w ≤X w ≤U w . 定义 8 对任意时间序列 Q 和 C, C w 为 C 的小波变换序列, U w 和 L w 为序列 Q 在小波域 的DTW 上、下边界序列 .如下定义 Q 和 C 的 d w LB距离: d w LB( Q, C) = ∑ m i =1 ( c w i -u w i ) 2 , c w i ≥ u w i ( l w i -c w i ) 2 , c w i ≤ l w i 0, l w i <c w i < u w i ( 12) 定理 4 任给时间序列 Q, C 和 X , U 和 L 为Q 的DTW 上边界序列和下边界序列, C w 为 C 的小波变换序列, X w 为 X 的小波变换序列 .如 果 L ≤X ≤U( li ≤x i ≤ui) , 则: d w LB( Q, C) ≤d( X w , C w ) . 证明 : d ( C w , X w ) = ∑ m i =1 ( c w i -x w i ) 2 = ∑ m i =1 ( c w i -x w i ) 2 , c w i ≥u w i ( c w i -x w i ) 2 , c w i ≤l w i ( c w i -x w i ) 2 , l w i <c w i <u w i 因为 L ≤X ≤U, 由定理 3 可得 L w ≤X w ≤ U w , 即 l w i ≤x w i ≤u w i , 所以当 c w i ≥u w i 时, 有( c w i - x w i ) 2 ≥( c w i -u w i ) 2 .当 c w i ≤u w i 时, 有( c w i -x w i ) 2 ≥( c w i -u w i ) 2 .当 l w i <c w i <u w i 时, 显然有( c w i - x w i ) 2 ≥0 .对照式 ( 12), 可得 d w LB ( C, Q ) ≤d ( X W , C W ) . 定理 5 任给时间序列 Q 和 C, 必有 d w LB ( Q, C) ≤D r DTW( Q, C) . 证明 :设 U 和L 为Q 的 DTW 上边界序列和 下边界序列, 由定理 2, 必存在 L ≤T ≤U (即 li ≤ ti ≤ui), 使得 d LB( Q, C) =d( T , C) . 设 C w 为 C 的小波变换序列, T w 为 T 的小 波变换序列, 由定理 4 可得: d w LB( Q, C) ≤d( T w , C w ) . 由于 Haar 为正交小波满足 parseval 定理, 可 得 d( T w , C w ) ≤d( T, C) . 所以, d w LB( Q, C) ≤d w LB( Q, C) . 又由定理 1, 可得 d r LB( Q, C) ≤D r DTW( Q, C), 所以, d w LB( Q, C) ≤D r DTW( Q, C) . 定理 6 采用小波域的 d w LB距离进行 DTW 距离相似匹配不会产生漏报. 证明 :设相似阈值为 ε, Q 为待查询序列, 任 给时间序列 C, 其弯曲范围为 r 的 DTW 距离为 D r DTW( Q, C) . 如果 D r DTW ( Q, C) <ε, 由定理 5, 必有 d w LB ( Q, C) <ε, 即使用 d w LB距离在小波域进行相似 匹配得到的候选序列必包含所有 D r DTW距离的匹 配序列, 不会产生漏报(即在小波域没有检索到却 满足查询条件的序列) . 定理 6 保证对于采用 DTW 距离的相似匹 配, 可以使用小波变换域定义的 d w LB距离检索, 不 会产生漏报, 但可能有误报(即在小波检索到却不 满足查询条件的序列), 误报可以通后处理来过滤 以得到到正确的检索结果 .该方法定义的低限距 离比先前方法相比更紧, 可最大限度的减小搜索 空间 .采用小波系数方法建立索引, 利用其降维 能力和多尺度特性提高搜索效率. 3 基于小波和 DTW 距离的相似匹 配算法 3.1 建立索引 首先对长度为 n( 设 n =2 J , 不足位补0)时间 序列 X 实施归一化, 以校正序列在 Y 轴的平移 和幅 度 伸缩 .即 x i = x i -x x max -x min , 其 中 x = 1 n ∑ n -1 i =0 xi , x max和 x min为序列的最大值和最小值. 然后对标准化后的序列实施离散小波变换, 取前 m 个小波系数组成小波系数序列, 视为 m 维空 间的一个点, 并将降维之后的小波系数序列组织 为多维索引结构 R ＊-tree 索引. R ＊-tree 是一种多级平衡树[ 6] , 树中的每个 非叶节点对应一个多维超矩形, 该超矩形为其子 节点代表超矩形的最小包围超矩形( MBR) .非叶 Vol.28 No.4 曲文龙等:基于小波和动态时间弯曲的时间序列相似匹配 · 399 ·

。400 北京科技大学学报 2006年第4期 节点由多个(SREC T,P)结构组成，其中P为子 离并按从小到大的顺序排列为(C,C2,,C). 节点指针，SRECT为与子节点P相关的MBR Stp4取(Ci,C2,；Ck)作为k最近邻查 叶节点含有多个(SRECT,O)结构，其中O为空 询结果 间对象的标识号.R"一tre中每个节点最多包含 一个点在小波域中基于欧氏距离的k最近 的结构个数称为Rte的度.索引建立过程如 邻点未必是基于DTW距离的k最近邻点（两序 下： 列的欧氏距离即弯曲度=0的DTW距离，不小 (1)对每一序列实施归一化. 于两序列弯曲度>0的DTW距离).因此要使 (2)每一序列实施离散Haar小波变换，取前 用范围查询进行修正.Step3中得到小波域的最 m个小波系数，得到降维之后的小波系数序列， 大欧氏距离e=d(Q,C),为进一步的DTW范 视为多维空间中的一个点 围查询提供了足够小的上界，可以避免过多的剪 (3)利用小波系数序列构建R*tee多维索 枝并保证不出现漏报. 引结构. 3.2范围查询算法 4实验结果 范围查询定义为对一个查询点Q=(91,92, 实验采用上交所股票数据，选取2003年上证 ；qm),检索Q与的DTW距离不超过e的点 180指数股和其它随机选取的70种股票2003-01 集，算法如下. 02到200412-31共241d的交易数据.选取每 Step1对长为n的查询序列(n维点)Q, 日收盘价得到250个维数为256的时间序列，进 按定义7的公式计算其在小波域的DTW上下边 行归一化、小波降维并构建R一tee索引.由于 界序列U"和L". 早期的基于DTW的索引和相似匹配算法效率很 Step2给定范围查询阈值e,则对应查询超 低不列入比较范围，只将本文方法和PAA方 矩形为： 法9进行比较.采用Haar小波，尺度取4，得到小 SRect(O)= 波系数为16维.使用同样的R*一tee维数比较 ([li-e,i+g,…,[m-e,um+e]）(12) 两方法的低限距离函数、范围查询精度、所需 Step3在R*-一tree索引结构内检索与查询 DTW距离计算比率. 超矩形SRect(Q)相交的m维小波域点的集合， 41低限距离比较 设共有L个点，并且它们对应的原始空间的(n 两序列的低限距离值小于其DTW距离，低 维)点集{C1,C2,,C}作为查询候选集. 限距离越接近DTW距离，则该低限距离的剪枝 Step4对查询候选集每一个点C(l≤i≤ 能力越强.随机选取50对序列，计算每两序列的 L),计算其与Q的DTW距离判断是否D5rw Haar低限距离、PAA低限距离和DTW距离，取 (C,Q~e,如果为假则认为C;是误报，为真则 50次计算的平均值.图3是弯曲路径P取不同 将C:加入到查询结果集. 值时的结果，可以看出低限距离的紧度随弯曲路 3.3最近邻查询算法 径增大而减小，但Haar低限距离值高于PAA低 k最近邻查询定义为对一个查询点Q=(91, 限距离，表明Haar低限距离可以剪除更多得误 92,“,qm),检索与Q的DTW距离最近的k个 报. 点，算法如下 3.f Step1对长度为n的查询序列Q进行小波 3.0 变换，取前m个小波系数组成Q的小波变换序 2.5 列Q. DTW 1.5 -e-Haar Step2在R一tree索引中调用快速最近邻 ★PAA 100 2468101214 算法了，得到小波域与Q的欧氏距离最近的k个 I)Tw弯h长度 点，对应到原始空间，按照它们与Q的距离从小 图3H圈r和PAA低限距离 到大的顺序排列为(C,C2,Ck). Fig.3 Haar and PAA lower bound distances Step3取e=d(Q,Ck),调用范围查询算 法进行点Q的DTW距离的e范围查询，得到原 42 范围查询精度比较 始空间中的L个点，分别计算其与Q的DTW距 范围查询精度定义为基于DTW的e查询的

节点由多个( SREC T, P ) 结构组成, 其中 P 为子 节点指针, SRECT 为与子节点 P 相关的 M BR. 叶节点含有多个( SRECT, O)结构, 其中 O 为空 间对象的标识号 .R ＊-tree 中每个节点最多包含 的结构个数称为 R ＊-tree 的度.索引建立过程如 下: ( 1) 对每一序列实施归一化 . ( 2) 每一序列实施离散 Haar 小波变换, 取前 m 个小波系数, 得到降维之后的小波系数序列, 视为多维空间中的一个点 . ( 3) 利用小波系数序列构建 R ＊-tree 多维索 引结构. 3.2 范围查询算法 范围查询定义为对一个查询点 Q =( q1, q2, …, qn ) , 检索 Q 与的 DTW 距离不超过 ε的点 集, 算法如下. Step 1 对长为 n 的查询序列( n 维点) Q, 按定义 7 的公式计算其在小波域的 DTW 上下边 界序列 U w 和 L w . Step 2 给定范围查询阈值 ε, 则对应查询超 矩形为: SRect( Q) = ([ l w 1 -ε, u w 1 +ε] , …,[ l w m -ε, u w m +ε] ) ( 12) Step 3 在 R ＊-tree 索引结构内检索与查询 超矩形 SRect( Q) 相交的 m 维小波域点的集合, 设共有 L 个点, 并且它们对应的原始空间的( n 维)点集{C1, C2, …, CL}作为查询候选集. Step 4 对查询候选集每一个点 Ci ( 1 ≤i ≤ L), 计算其与 Q 的 DTW 距离判断是否 D r DT W ( Ci , Q) <ε, 如果为假则认为 Ci 是误报, 为真则 将 Ci 加入到查询结果集. 3.3 最近邻查询算法 k 最近邻查询定义为对一个查询点 Q =( q1, q2, …, qn) , 检索与 Q 的 DTW 距离最近的 k 个 点, 算法如下. Step 1 对长度为 n 的查询序列Q 进行小波 变换, 取前 m 个小波系数组成 Q 的小波变换序 列Q w . Step 2 在 R ＊-tree 索引中调用快速最近邻 算法[ 7] , 得到小波域与 Q 的欧氏距离最近的 k 个 点, 对应到原始空间, 按照它们与 Q 的距离从小 到大的顺序排列为( C ＊ 1 , C ＊ 2 , …, C ＊ k ) . Step 3 取 ε=d ( Q, C ＊ k ), 调用范围查询算 法进行点 Q 的 DTW 距离的 ε范围查询, 得到原 始空间中的 L 个点, 分别计算其与 Q 的 DTW 距 离并按从小到大的顺序排列为( C1, C2, …, CL ) . S tep 4 取( C1, C2, …, Ck )作为 k 最近邻查 询结果. 一个点在小波域中基于欧氏距离的 k 最近 邻点未必是基于 DTW 距离的 k 最近邻点(两序 列的欧氏距离即弯曲度 r =0 的 DTW 距离, 不小 于两序列弯曲度 r >0 的 DTW 距离) .因此要使 用范围查询进行修正 .Step 3 中得到小波域的最 大欧氏距离 ε=d ( Q, C ＊ k ), 为进一步的 DTW 范 围查询提供了足够小的上界, 可以避免过多的剪 枝并保证不出现漏报 . 4 实验结果 实验采用上交所股票数据, 选取 2003 年上证 180 指数股和其它随机选取的 70 种股票 2003-01 -02 到 2004-12-31 共241 d 的交易数据 .选取每 日收盘价得到 250 个维数为 256 的时间序列, 进 行归一化、小波降维并构建 R ＊-tree 索引 .由于 早期的基于 DTW 的索引和相似匹配算法效率很 低, 不列入比较范围, 只将本文方法和 PAA 方 法[ 5] 进行比较.采用Haar 小波, 尺度取 4, 得到小 波系数为 16 维 .使用同样的 R ＊-tree 维数比较 两方法的低限距离函数 、范围查询精度、所需 DTW 距离计算比率. 4.1 低限距离比较 两序列的低限距离值小于其 DTW 距离, 低 限距离越接近 DTW 距离, 则该低限距离的剪枝 能力越强.随机选取 50 对序列, 计算每两序列的 Haar 低限距离、PAA 低限距离和 DTW 距离, 取 50 次计算的平均值 .图 3 是弯曲路径 ρ取不同 值时的结果, 可以看出低限距离的紧度随弯曲路 径增大而减小, 但 Haar 低限距离值高于 PAA 低 限距离, 表明 Haar 低限距离可以剪除更多得误 报. 图 3 Haar 和 PAA 低限距离 Fig.3 Haar and PAA lower bound distances 4.2 范围查询精度比较 范围查询精度定义为基于 DTW 的 ε查询的 · 400 · 北 京 科 技 大 学 学 报 2006 年第 4 期