多尺度直线拟合法在时间序列突变点检测中的应用

黄静，李长春，延皓，赵旭昌，杨雪松

(北京交通大学 机械与电子控制工程学院，北京100044)

0 引言

电液伺服阀在战机中存在较多应用，对一些液压动作器件起控制作用，因此电液伺服阀在应用之前必须要进行性能测试和筛选。测试实验主要是依据控制电流-压力曲线，获取相关指标。以往实验结果的读取是通过人工的方法从曲线图上判断突变点或者转折点，从而获得结果，在有噪声干扰的情况下存在较大的人工读数误差，效率也较低，不利于测试实验的自动化。因此需要找到一种有效的突变点的计算和判断方法，利用计算机技术来进行自动处理。

时间序列中当数据从一种统计趋势变化到另一种统计趋势时存在的不连续点即突变点。目前在突变点的检测上，国内外很多学者都进行了相关研究［1-3］，也提出了很多检测的方法。突变点检测在金融、气候、水文、故障检测等诸多领域中均有广泛应用［4-6］。

目前在突变点的检测方法应用较多的有Mann-Kendall 算法、累积和控制图(CUSUM)算法［7］、最小均方差(MSE)法、小波变换法等。其中Mann-Kendall 算法及其改进算法在气候和水文领域应用较多。Mann-Kendall 算法、CUSUM 算法和MSE 算法针对无趋势变化的序列时检测效果很好，但对于有趋势变化的序列时检测出的突变点相比于实际数据有较大出入。小波变换法中小波变换的系数选择、所选用小波和噪声干扰以及具体的一些参数确定上，会对检测结果造成一定影响，并且采用小波变换法时，需要进行多分辨率的逼近，因此算法计算量较大，耗时较多。

为了平衡计算效率和检测精度，可用变换尺度的方法来逐步逼近时间序列中的突变点。基于此种思路，本文提出了一种多尺度直线拟合法。

1 多尺度直线拟合方法

1.1 初始尺度的计算方法

对于时间序列x1，x2，x3，…，xn，不论是有趋势变换的，还是无趋势变化，均可以按照一定的尺度，利用多段拟合直线来代替。只要尺度选择适当，就能很好地反映出原序列的信息。通过在一定尺度下利用多段拟合线段代替原始序列的方法，能提取出原时间序列的主要信息，舍弃一些不必要的信息，简化问题，从而减少计算量。

初始尺度的选择在拟合线段代替原时间序列的过程中比较重要。如果尺度选择过大，则使得拟合线段过于粗糙，丢失信息过多，不能很好地反映原始信息，从而导致错误的检测结果。如果尺度过小，则会引入过多的计算，或者因为原序列中某一部分的局部特征过于明显而影响整个分析过程。初始尺度选择过大或过小的效果如图1(a)和图1(b)所示。

图1 尺度不当拟合结果图Fig.1 Improper scale diagram

对于初始尺度l1的选择采用(1)式的计算方法:

式中:n 为整个时间序列的长度。

将图1中两组时间序列采用(1)式的计算方法后选择的初始尺度进行拟合后的效果图如图2所示。从图2中可以看出:采用(1)式后的尺度能较好地拟合原序列，并能有效排除局部特性过于强烈对整体的影响。对于无趋势变化的时间序列采用(1)式的计算方法后也具有较好的拟合替代效果，如图3所示。

图2 合适尺度拟合结果图Fig.2 Proper scale diagrams

图3 无趋势序列拟合效果图Fig.3 Fitting results without trend time series

1.2 直线拟合方法

采用(1)式计算得出初始尺度以后，可以将整个时间序列x1，x2，x3，…，xn分为m 段，每一段用sj表示，j=1，2，…，m，其长度为l1最后一段的长度如果小于l1，可将其纳入前一段数据中。

对于sj里面的数据序列，采用最小二乘法进行直线拟合。

设分段后的原始序列sj的拟合直线表达式为x'i=k·i+b，i=1，2，……，n，则序列sj和拟合直线的残差平方和为

(3)式的右边可以展开为

(5)式中除了k、b 为未知量以外，其他符号均为已知数据。根据最小二乘法的原理，要使得残差平方和RSS 为最小，必须且只需且其中分别为i 、xi的平均值。由此可以求出序列sj拟合直线的斜率k 和截距b，则求得序列s1，s2，s3，…，sm的每一段拟合直线的斜率依次为k1，k2，k3，…，km.

1.3 突变点的求取

求出每一段序列的斜率k1，k2，k3，…，km以后，依次求出其差值的绝对值ej= |kj+1- kj|，1 ≤j≤m-1，假设emax=ej，说明原始序列s1，s2，s3，…，sm在sj与sj+1之间的变化趋势最大，则相应的突变点一定落在sj、sj+1所包含的区间(j×l1，j×l1+l1)内。因此在区间(j×l1，j×l1+l1)内变换尺度，将尺度由l1变为l2，l2=floor(l1×0.5). 在尺度l2下针对新的区间再次采用1.2 节和1.3 节的方法逐渐逼近突变点，直至最终尺度ln=2，此时搜寻区间长度为3，区间内的中点便是时间序列x1，x2，x3，…，xn的突变点。

2 同其他方法的对比

在寻找突变点的计算方法上有很多应用广泛的成熟算法，如滑动t-检验算法、Cramer’s 算法、Yamamoto 算法、Pettitt 算法、Lepage 算法、Mann-Kendall 算法及其改进算法［8］、CUSUM 算法、小波变换法等，其中应用较多的有Mann-Kendall 算法、CUSUM 算法和小波变换法。

2.1 Mann-Kendall 算法

Mann-Kendall 算法是一种非参数统计检验方法，即无分布检验法，其优点是待处理序列不需要遵从一定的分布规律，不受少数异常值的干扰，计算也较为简便。具体计算方法如下［4，8］。

对于时间序列x1，x2，x3，…，xn，可以构造一秩序列Si，i=1，2，3，…，n.

式中:

可见Si是第i 时刻值大于前各个时刻数值个数的累加数。在时间序列独立随机的的假设下，定义一个统计量UFi:

式中:UF1=0;E(Si)、Var(Si)分别为Si的均值和方差，可由(8)式和(9)式进行计算:

UFi为标准正态分布，给定显著性水平α，查对应的正态分布表，若|UFi| ＞Uα，则说明该序列存在明显的趋势变化。UFi，i=1，2，3，…，n，可组成一条曲线c1.

针对时间序列的逆序xn，xn-1，xn-2，…，x1，再重复上述过程，求得UFp，令UBp= -UFp，p =n，n -1，n-2，…，1，同样UBp，p =n，n -1，n -2，…，1，可以组成一条曲线c2. 如果c1，c2曲线有交叉，且交点位于显著性水平α 所在的置信区间内，说明该交叉点就是突变点。

利用该方法，检测1900年～1990年上海年平均气温突变点的结果如图4所示，显著性水平α 取值0.05. 从图4中可以看出Mann-Kendall 算法可以很好地将1900年～1990年上海年平均气温突变点检测出来。

图4 Mann-Kendall 算法处理结果图Fig.4 The results from Mann-Kendall method

2.2 CUSUM 算法

CUSUM 算法是利用累积时间序列中每个点与序列平均值的差值的和的方法来进行突变点的判断。其计算方法如下［7］:

1)首先计算出整个序列x1，x2，x3，…，xn的平均值

2)令累积和Ti为

3)判断突变点的位置。进行判断处理时有两种方法:第1 种方法是从Ti中找出最大值Tmax，其对应的点xmax就是突变点发生开始的位置。即|Tmax| =max(|Ti|);第2 种方法是采用MSE 的方式来进行判断。

定义

式中:

若从MSEg中取得最小值MSEmin，突变点就在xmin点处。

2.3 小波变换法

小波变换法检测突变点的基本思想是:对待分析信号进行多级小波分解，并只对小波的第一层高频系数进行信号重构，重构后的信号可以明显地表现出信号的局部特征，在局部范围内变换后的信号极大值或者过零值就可以看作是信号的突变点［9］。如果检测t0点为奇异点，那么在该点处小波变换取得极大值，可以通过对模量极大值点的检测来确定突变点发生的时间点［10］。

其具体做法是:第1 步采用小波去噪的方法，去除信号中的噪声信号;第2 步是利用小波变换来进行突变点的检测。

值得注意的是，在利用小波变换来检测信号的突变点过程中，小波变换系数的选择、所选用小波和噪声干扰，以及具体的一些参数确定上，会对检测结果有一定的影响。检测结果的读取需要人工参与判断，其去噪方法和如何更好地确定小波模极大值或过零值都需要再深入研究。

2.4 算法比较

利用多尺度直线拟合方法寻找突变点的过程其实就是通过变换尺度，从大到小、从粗到细不断逼近突变点的过程。在求解过程中因为会不断缩小计算范围，因此计算量较小，其花费时间比较少。假设时间序列的长度为n，初始尺度为则计算次数极端情况下为

而Mann-Kendall 算法，针对序列中的每一个点都要循环同当前点之前和之后的点进行数值大小的比较并计算累积和，正序计算完成后还要对序列的逆序进行计算，因此运算量较大，花费时间较长。针对长度为n 的时间序列，其计算次数为

CUSUM 算法在使用最小MSE 检测突变点也同样存在类似的运算量较大问题，因为该算法也会在当前点处循环计算之前和之后点的累积和。CUSUM 算法使用最大累积和判断法时计算次数为

使用最小MSE 判断法时其计算次数为

因此从效率方面来讲，多尺度直线拟合法的计算效率要高于Mann-Kendall 算法和CUSUM 算法。针对同样一段长度为400 的时间序列，各种算法的耗时如表1所示。从表1中可以看出，多尺度直线拟合法的耗时要多于CUSUM 算法采用最大累积和法的时间，这是因为多尺度直线拟合法在每次计算的时候还需要利用最小二乘法进行直线拟合，因此耗时要多一些。但即使这样，其耗时也远小于Mann-Kendall 算法和CUSUM 算法采用最小方差法。

表1 算法耗时对比表Tab.1 Comparison of time consuming

除了效率的对比，还需要从精度上来对以上4 种算法进行对比，为了使对比效果更加明显，采用了4 组不同的数据来进行验证。第1 组为无趋势变化无噪声干扰的序列，第2 组为无趋势变化有噪声干扰的序列，第3 组为有递增趋势的无噪声干扰信号，第4 组信号为有递增趋势有噪声干扰的信号。为了检验噪声对4 种算法的干扰程度，对比实验中加入较强的噪声干扰，其信噪比SNR 为29.68. 4 种方法求得的突变点位置如表2所示。第1 组和第2 组数据的效果如图5所示。第3 组和第4 组数据的效果如图6所示。从图5、图6和表2中可以看出，多尺度直线拟合法对于有趋势变化序列和无趋势变化序列都能有很好地检测结果，在有较强噪声干扰的情况下也能有很好地检测效果。通过实验结果的对比，说明了多尺度直线拟合法的检测精度和效果都要优于其他3 种算法。

表2 检测结果对比表Tab.2 Comparison of detection results

图5 无趋势变化及噪声干扰效果图Fig.5 Time series without trend change and noise

图6 有趋势变化及噪声干扰效果图Fig.6 Time series with trend change and noise

3 应用实例

在针对伺服阀性能的实际测量中，需要测得一个完整的伺服阀的“电流-压力”曲线，即控制电流从零值开始到最大，再从最大恢复到零值的过程，通过测得的曲线可以看出控制电流与压力之间的关系，并可以从曲线图上测得伺服阀的死区区间和分辨率区间，从而进一步评价被测伺服阀的相关性能。而死区和分辨率的判断都是从曲线图上寻找相关突变点作为区间边界来进行判定。

在实验测试系统中采用本文介绍的多尺度直线拟合法进行突变点的检测，其检测结果如图7所示。在实际实验情况下，连续做了5 组实验，针对5 组实验曲线，多尺度直线拟合法的检测突变点的结果如表3所示。从图7和表3可以看出，在有噪声存在的情况下，多尺度直线拟合法具有很好的检测效果和检测精度。

图7 实际数据检测效果图Fig.7 Real data detection results

表3 实验数据检测结果Tab.3 Detection results of experimental data

为了进一步验证该方法的有效性，对行星齿轮箱的振动信号进行检测处理。该行星齿轮箱的时域波形图，如图8所示，间隔为0.1 s，即频率为10 Hz.从时域振动信号可以看出，时域中因为行星架转动引起的周期性的冲击。从图8中可以看出，在时间区间［0.3 s 0.4 s］和［0.7 s 0.8 s］之间有较强的冲击产生，属于故障区间。为了找出其中的故障信号以及定位故障区间，使用本文提出的多尺度直线拟合法进行查找。查找出的信号突变点在图8中标出。从查找结果可以看出，该方法能较好地定位到故障区间。但是由于噪声信号干扰较大及信号较复杂，导致某些区间故障点的定位不是特别精准。

图8 行星齿轮箱振动信号及检测结果Fig.8 Vibration signal of planetary gear and detection result

针对傅里叶变换后的频谱信号，也可以用该检测方法进行频谱中最大分量的检测。对于两个正弦信号的叠加信号x，x=0.6 ×sin(2π×50t)+0.35 ×sin(2π ×120t)，再叠加一个比较大的噪声信号，其信噪比为-9.269 9 dB. 对添加噪声以后的信号进行快速傅里叶变换(FFT)，得到其频谱图，如图9所示。从图中可以看出，经过FFT 后，能很明显看到在50 Hz 和120 Hz 有较大的频率分量。使用本文的“多尺度直线拟合法”进行突变点的检测，从检测结果来看，该方法可以很准确地检测到频谱中最大的频率分量。

图9 频谱图及检测结果Fig.9 Frequency spectrogram and detection result

4 结论

多尺度直线拟合法检测突变点其本质就是通过尺度变换，并不断缩小检测范围，利用拟合线段斜率的变化从而最终逼近时间序列中的突变点。而初始尺度的选择对检测结果有较大影响，本文给出了一种有效的初始尺度计算方法。将多尺度直线拟合法同目前应用较多的其他3 种方法进行比较，证明了该方法的计算效率、准确度以及抗干扰方面都具有一定的优越性。将多尺度直线拟合检测方法应用于实际的实验系统中，也取得了很好的效果，提供了检测效率和检测精度。对于初始尺度的计算方法还可以再进一步进行研究。

References)

［1］Burke M D，Bewa G. Change-point detection for general nonparametric regression models［J］. Open Journal of Statistics，2013，3(4):261 -267.

［2］Habibi R. A note on change point detection using weighted least square［J］.Applied Mathematics，2011，2 (10):1309 -1312.

［3］Inoue S，Yamada S. A bivariate software reliability model with change-point and its applications［J］. American Journal of Operations Research，2011(1):1 -7.

［4］王金花，张荣刚，张攀，等.两种水沙系列突变点算法的对比分析［J］. 中国水土保持，2009(12):43 -44.WANG Jin-hua，ZHANG Rong-gang，ZHANG Pan，et al. Comparison on catastrophe point calculation of two water and sediment series［J］. Soil and Water Conservation in China，2009(12):43 -44.(in Chinese)

［5］刘嘉琦，龚政，张长宽.重大水利工程运用对长江入海径流量的影响［J］.水道港口，2013，34(6):461 -466.LIU Jia-qi，GONG Zheng，ZHANG Chang-kuan. Impact of largescale water projects on the Yangtze river［J］.Journal of Waterway and Harbor，2013，34(6):461 -466.(in Chinese)

［6］陈远中，陆宝宏，张育德，等.改进的有序聚类分析法提取时间序列转折点［J］.水文，2011，31(1):41 -44.CHEN Yuan-zhong，LU Bao-hong，ZHANG Yu-de，et al. Improvement of sequential cluster analysis method for extracting turning point of time series［J］.Journal of China Hydrology，2011，31(1):41 -44.(in Chinese)

［7］Khan A，Chatterjee S，Bisai D，et a1.Analysis of change point in surface temperature time series using cumulative sum chart and bootstrapping for Asansol weather observation station，West Bengal，India［J］.American Journal of Climate Change，2014(3):83-94.

［8］魏凤英.现代气候统计诊断与预测技术［M］. 北京:气象出版社，2009:62 -73.WEI Feng-ying. Modern climatological statistical diagnosis and prediction technology［M］. Beijing:Weteorological Press，2009:62 -73.(in Chinese)

［9］吴凡，熊高君，叶志婵.小波变换在信号突变点检测中的应用［J］.计算机与现代化，2008(8):133 -135.WU Fan，XIONG Gao-jun，YE Zhi-chan. Applications of wavelet transform on signals catastrophe-points detection［J］.Computer and Modernization，2008(8):133 -135.(in Chinese)

［10］赖富文，张志杰，刘景江，等. 基于炮口脉冲噪声信号的射速测试方法研究［J］. 兵工学报，2013，34(9):1180 -1186.LAI Fu-wen，ZHANG Zhi-jie，LIU Jing-jiang，et al. Research on testing firing rate based on muzzle impulse noise［J］.Acta Armamentarii，2013，34(9):1180 -1186.(in Chinese)