多尺度直線擬合法在時(shí)間序列突變點(diǎn)檢測中的應(yīng)用

黃靜，李長春，延皓，趙旭昌，楊雪松

(北京交通大學(xué) 機(jī)械與電子控制工程學(xué)院，北京100044)

0 引言

電液伺服閥在戰(zhàn)機(jī)中存在較多應(yīng)用，對(duì)一些液壓動(dòng)作器件起控制作用，因此電液伺服閥在應(yīng)用之前必須要進(jìn)行性能測試和篩選。測試實(shí)驗(yàn)主要是依據(jù)控制電流-壓力曲線，獲取相關(guān)指標(biāo)。以往實(shí)驗(yàn)結(jié)果的讀取是通過人工的方法從曲線圖上判斷突變點(diǎn)或者轉(zhuǎn)折點(diǎn)，從而獲得結(jié)果，在有噪聲干擾的情況下存在較大的人工讀數(shù)誤差，效率也較低，不利于測試實(shí)驗(yàn)的自動(dòng)化。因此需要找到一種有效的突變點(diǎn)的計(jì)算和判斷方法，利用計(jì)算機(jī)技術(shù)來進(jìn)行自動(dòng)處理。

時(shí)間序列中當(dāng)數(shù)據(jù)從一種統(tǒng)計(jì)趨勢變化到另一種統(tǒng)計(jì)趨勢時(shí)存在的不連續(xù)點(diǎn)即突變點(diǎn)。目前在突變點(diǎn)的檢測上，國內(nèi)外很多學(xué)者都進(jìn)行了相關(guān)研究［1-3］，也提出了很多檢測的方法。突變點(diǎn)檢測在金融、氣候、水文、故障檢測等諸多領(lǐng)域中均有廣泛應(yīng)用［4-6］。

目前在突變點(diǎn)的檢測方法應(yīng)用較多的有Mann-Kendall 算法、累積和控制圖(CUSUM)算法［7］、最小均方差(MSE)法、小波變換法等。其中Mann-Kendall 算法及其改進(jìn)算法在氣候和水文領(lǐng)域應(yīng)用較多。Mann-Kendall 算法、CUSUM 算法和MSE 算法針對(duì)無趨勢變化的序列時(shí)檢測效果很好，但對(duì)于有趨勢變化的序列時(shí)檢測出的突變點(diǎn)相比于實(shí)際數(shù)據(jù)有較大出入。小波變換法中小波變換的系數(shù)選擇、所選用小波和噪聲干擾以及具體的一些參數(shù)確定上，會(huì)對(duì)檢測結(jié)果造成一定影響，并且采用小波變換法時(shí)，需要進(jìn)行多分辨率的逼近，因此算法計(jì)算量較大，耗時(shí)較多。

為了平衡計(jì)算效率和檢測精度，可用變換尺度的方法來逐步逼近時(shí)間序列中的突變點(diǎn)。基于此種思路，本文提出了一種多尺度直線擬合法。

1 多尺度直線擬合方法

1.1 初始尺度的計(jì)算方法

對(duì)于時(shí)間序列x1，x2，x3，…，xn，不論是有趨勢變換的，還是無趨勢變化，均可以按照一定的尺度，利用多段擬合直線來代替。只要尺度選擇適當(dāng)，就能很好地反映出原序列的信息。通過在一定尺度下利用多段擬合線段代替原始序列的方法，能提取出原時(shí)間序列的主要信息，舍棄一些不必要的信息，簡化問題，從而減少計(jì)算量。

初始尺度的選擇在擬合線段代替原時(shí)間序列的過程中比較重要。如果尺度選擇過大，則使得擬合線段過于粗糙，丟失信息過多，不能很好地反映原始信息，從而導(dǎo)致錯(cuò)誤的檢測結(jié)果。如果尺度過小，則會(huì)引入過多的計(jì)算，或者因?yàn)樵蛄兄心骋徊糠值木植刻卣鬟^于明顯而影響整個(gè)分析過程。初始尺度選擇過大或過小的效果如圖1(a)和圖1(b)所示。

圖1 尺度不當(dāng)擬合結(jié)果圖Fig.1 Improper scale diagram

對(duì)于初始尺度l1的選擇采用(1)式的計(jì)算方法:

式中:n 為整個(gè)時(shí)間序列的長度。

將圖1中兩組時(shí)間序列采用(1)式的計(jì)算方法后選擇的初始尺度進(jìn)行擬合后的效果圖如圖2所示。從圖2中可以看出:采用(1)式后的尺度能較好地?cái)M合原序列，并能有效排除局部特性過于強(qiáng)烈對(duì)整體的影響。對(duì)于無趨勢變化的時(shí)間序列采用(1)式的計(jì)算方法后也具有較好的擬合替代效果，如圖3所示。

圖2 合適尺度擬合結(jié)果圖Fig.2 Proper scale diagrams

圖3 無趨勢序列擬合效果圖Fig.3 Fitting results without trend time series

1.2 直線擬合方法

采用(1)式計(jì)算得出初始尺度以后，可以將整個(gè)時(shí)間序列x1，x2，x3，…，xn分為m 段，每一段用sj表示，j=1，2，…，m，其長度為l1最后一段的長度如果小于l1，可將其納入前一段數(shù)據(jù)中。

對(duì)于sj里面的數(shù)據(jù)序列，采用最小二乘法進(jìn)行直線擬合。

設(shè)分段后的原始序列sj的擬合直線表達(dá)式為x'i=k·i+b，i=1，2，……，n，則序列sj和擬合直線的殘差平方和為

(3)式的右邊可以展開為

(5)式中除了k、b 為未知量以外，其他符號(hào)均為已知數(shù)據(jù)。根據(jù)最小二乘法的原理，要使得殘差平方和RSS 為最小，必須且只需且其中分別為i 、xi的平均值。由此可以求出序列sj擬合直線的斜率k 和截距b，則求得序列s1，s2，s3，…，sm的每一段擬合直線的斜率依次為k1，k2，k3，…，km.

1.3 突變點(diǎn)的求取

求出每一段序列的斜率k1，k2，k3，…，km以后，依次求出其差值的絕對(duì)值ej= |kj+1- kj|，1 ≤j≤m-1，假設(shè)emax=ej，說明原始序列s1，s2，s3，…，sm在sj與sj+1之間的變化趨勢最大，則相應(yīng)的突變點(diǎn)一定落在sj、sj+1所包含的區(qū)間(j×l1，j×l1+l1)內(nèi)。因此在區(qū)間(j×l1，j×l1+l1)內(nèi)變換尺度，將尺度由l1變?yōu)閘2，l2=floor(l1×0.5). 在尺度l2下針對(duì)新的區(qū)間再次采用1.2 節(jié)和1.3 節(jié)的方法逐漸逼近突變點(diǎn)，直至最終尺度ln=2，此時(shí)搜尋區(qū)間長度為3，區(qū)間內(nèi)的中點(diǎn)便是時(shí)間序列x1，x2，x3，…，xn的突變點(diǎn)。

2 同其他方法的對(duì)比

在尋找突變點(diǎn)的計(jì)算方法上有很多應(yīng)用廣泛的成熟算法，如滑動(dòng)t-檢驗(yàn)算法、Cramer’s 算法、Yamamoto 算法、Pettitt 算法、Lepage 算法、Mann-Kendall 算法及其改進(jìn)算法［8］、CUSUM 算法、小波變換法等，其中應(yīng)用較多的有Mann-Kendall 算法、CUSUM 算法和小波變換法。

2.1 Mann-Kendall 算法

Mann-Kendall 算法是一種非參數(shù)統(tǒng)計(jì)檢驗(yàn)方法，即無分布檢驗(yàn)法，其優(yōu)點(diǎn)是待處理序列不需要遵從一定的分布規(guī)律，不受少數(shù)異常值的干擾，計(jì)算也較為簡便。具體計(jì)算方法如下［4，8］。

對(duì)于時(shí)間序列x1，x2，x3，…，xn，可以構(gòu)造一秩序列Si，i=1，2，3，…，n.

式中:

可見Si是第i 時(shí)刻值大于前各個(gè)時(shí)刻數(shù)值個(gè)數(shù)的累加數(shù)。在時(shí)間序列獨(dú)立隨機(jī)的的假設(shè)下，定義一個(gè)統(tǒng)計(jì)量UFi:

式中:UF1=0;E(Si)、Var(Si)分別為Si的均值和方差，可由(8)式和(9)式進(jìn)行計(jì)算:

UFi為標(biāo)準(zhǔn)正態(tài)分布，給定顯著性水平α，查對(duì)應(yīng)的正態(tài)分布表，若|UFi| ＞Uα，則說明該序列存在明顯的趨勢變化。UFi，i=1，2，3，…，n，可組成一條曲線c1.

針對(duì)時(shí)間序列的逆序xn，xn-1，xn-2，…，x1，再重復(fù)上述過程，求得UFp，令UBp= -UFp，p =n，n -1，n-2，…，1，同樣UBp，p =n，n -1，n -2，…，1，可以組成一條曲線c2. 如果c1，c2曲線有交叉，且交點(diǎn)位于顯著性水平α 所在的置信區(qū)間內(nèi)，說明該交叉點(diǎn)就是突變點(diǎn)。

利用該方法，檢測1900年～1990年上海年平均氣溫突變點(diǎn)的結(jié)果如圖4所示，顯著性水平α 取值0.05. 從圖4中可以看出Mann-Kendall 算法可以很好地將1900年～1990年上海年平均氣溫突變點(diǎn)檢測出來。

圖4 Mann-Kendall 算法處理結(jié)果圖Fig.4 The results from Mann-Kendall method

2.2 CUSUM 算法

CUSUM 算法是利用累積時(shí)間序列中每個(gè)點(diǎn)與序列平均值的差值的和的方法來進(jìn)行突變點(diǎn)的判斷。其計(jì)算方法如下［7］:

1)首先計(jì)算出整個(gè)序列x1，x2，x3，…，xn的平均值

2)令累積和Ti為

3)判斷突變點(diǎn)的位置。進(jìn)行判斷處理時(shí)有兩種方法:第1 種方法是從Ti中找出最大值Tmax，其對(duì)應(yīng)的點(diǎn)xmax就是突變點(diǎn)發(fā)生開始的位置。即|Tmax| =max(|Ti|);第2 種方法是采用MSE 的方式來進(jìn)行判斷。

定義

式中:

若從MSEg中取得最小值MSEmin，突變點(diǎn)就在xmin點(diǎn)處。

2.3 小波變換法

小波變換法檢測突變點(diǎn)的基本思想是:對(duì)待分析信號(hào)進(jìn)行多級(jí)小波分解，并只對(duì)小波的第一層高頻系數(shù)進(jìn)行信號(hào)重構(gòu)，重構(gòu)后的信號(hào)可以明顯地表現(xiàn)出信號(hào)的局部特征，在局部范圍內(nèi)變換后的信號(hào)極大值或者過零值就可以看作是信號(hào)的突變點(diǎn)［9］。如果檢測t0點(diǎn)為奇異點(diǎn)，那么在該點(diǎn)處小波變換取得極大值，可以通過對(duì)模量極大值點(diǎn)的檢測來確定突變點(diǎn)發(fā)生的時(shí)間點(diǎn)［10］。

其具體做法是:第1 步采用小波去噪的方法，去除信號(hào)中的噪聲信號(hào);第2 步是利用小波變換來進(jìn)行突變點(diǎn)的檢測。

值得注意的是，在利用小波變換來檢測信號(hào)的突變點(diǎn)過程中，小波變換系數(shù)的選擇、所選用小波和噪聲干擾，以及具體的一些參數(shù)確定上，會(huì)對(duì)檢測結(jié)果有一定的影響。檢測結(jié)果的讀取需要人工參與判斷，其去噪方法和如何更好地確定小波模極大值或過零值都需要再深入研究。

2.4 算法比較

利用多尺度直線擬合方法尋找突變點(diǎn)的過程其實(shí)就是通過變換尺度，從大到小、從粗到細(xì)不斷逼近突變點(diǎn)的過程。在求解過程中因?yàn)闀?huì)不斷縮小計(jì)算范圍，因此計(jì)算量較小，其花費(fèi)時(shí)間比較少。假設(shè)時(shí)間序列的長度為n，初始尺度為則計(jì)算次數(shù)極端情況下為

而Mann-Kendall 算法，針對(duì)序列中的每一個(gè)點(diǎn)都要循環(huán)同當(dāng)前點(diǎn)之前和之后的點(diǎn)進(jìn)行數(shù)值大小的比較并計(jì)算累積和，正序計(jì)算完成后還要對(duì)序列的逆序進(jìn)行計(jì)算，因此運(yùn)算量較大，花費(fèi)時(shí)間較長。針對(duì)長度為n 的時(shí)間序列，其計(jì)算次數(shù)為

CUSUM 算法在使用最小MSE 檢測突變點(diǎn)也同樣存在類似的運(yùn)算量較大問題，因?yàn)樵撍惴ㄒ矔?huì)在當(dāng)前點(diǎn)處循環(huán)計(jì)算之前和之后點(diǎn)的累積和。CUSUM 算法使用最大累積和判斷法時(shí)計(jì)算次數(shù)為

使用最小MSE 判斷法時(shí)其計(jì)算次數(shù)為

因此從效率方面來講，多尺度直線擬合法的計(jì)算效率要高于Mann-Kendall 算法和CUSUM 算法。針對(duì)同樣一段長度為400 的時(shí)間序列，各種算法的耗時(shí)如表1所示。從表1中可以看出，多尺度直線擬合法的耗時(shí)要多于CUSUM 算法采用最大累積和法的時(shí)間，這是因?yàn)槎喑叨戎本€擬合法在每次計(jì)算的時(shí)候還需要利用最小二乘法進(jìn)行直線擬合，因此耗時(shí)要多一些。但即使這樣，其耗時(shí)也遠(yuǎn)小于Mann-Kendall 算法和CUSUM 算法采用最小方差法。

表1 算法耗時(shí)對(duì)比表Tab.1 Comparison of time consuming

除了效率的對(duì)比，還需要從精度上來對(duì)以上4 種算法進(jìn)行對(duì)比，為了使對(duì)比效果更加明顯，采用了4 組不同的數(shù)據(jù)來進(jìn)行驗(yàn)證。第1 組為無趨勢變化無噪聲干擾的序列，第2 組為無趨勢變化有噪聲干擾的序列，第3 組為有遞增趨勢的無噪聲干擾信號(hào)，第4 組信號(hào)為有遞增趨勢有噪聲干擾的信號(hào)。為了檢驗(yàn)噪聲對(duì)4 種算法的干擾程度，對(duì)比實(shí)驗(yàn)中加入較強(qiáng)的噪聲干擾，其信噪比SNR 為29.68. 4 種方法求得的突變點(diǎn)位置如表2所示。第1 組和第2 組數(shù)據(jù)的效果如圖5所示。第3 組和第4 組數(shù)據(jù)的效果如圖6所示。從圖5、圖6和表2中可以看出，多尺度直線擬合法對(duì)于有趨勢變化序列和無趨勢變化序列都能有很好地檢測結(jié)果，在有較強(qiáng)噪聲干擾的情況下也能有很好地檢測效果。通過實(shí)驗(yàn)結(jié)果的對(duì)比，說明了多尺度直線擬合法的檢測精度和效果都要優(yōu)于其他3 種算法。

表2 檢測結(jié)果對(duì)比表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 應(yīng)用實(shí)例

在針對(duì)伺服閥性能的實(shí)際測量中，需要測得一個(gè)完整的伺服閥的“電流-壓力”曲線，即控制電流從零值開始到最大，再從最大恢復(fù)到零值的過程，通過測得的曲線可以看出控制電流與壓力之間的關(guān)系，并可以從曲線圖上測得伺服閥的死區(qū)區(qū)間和分辨率區(qū)間，從而進(jìn)一步評(píng)價(jià)被測伺服閥的相關(guān)性能。而死區(qū)和分辨率的判斷都是從曲線圖上尋找相關(guān)突變點(diǎn)作為區(qū)間邊界來進(jìn)行判定。

在實(shí)驗(yàn)測試系統(tǒng)中采用本文介紹的多尺度直線擬合法進(jìn)行突變點(diǎn)的檢測，其檢測結(jié)果如圖7所示。在實(shí)際實(shí)驗(yàn)情況下，連續(xù)做了5 組實(shí)驗(yàn)，針對(duì)5 組實(shí)驗(yàn)曲線，多尺度直線擬合法的檢測突變點(diǎn)的結(jié)果如表3所示。從圖7和表3可以看出，在有噪聲存在的情況下，多尺度直線擬合法具有很好的檢測效果和檢測精度。

圖7 實(shí)際數(shù)據(jù)檢測效果圖Fig.7 Real data detection results

表3 實(shí)驗(yàn)數(shù)據(jù)檢測結(jié)果Tab.3 Detection results of experimental data

為了進(jìn)一步驗(yàn)證該方法的有效性，對(duì)行星齒輪箱的振動(dòng)信號(hào)進(jìn)行檢測處理。該行星齒輪箱的時(shí)域波形圖，如圖8所示，間隔為0.1 s，即頻率為10 Hz.從時(shí)域振動(dòng)信號(hào)可以看出，時(shí)域中因?yàn)樾行羌苻D(zhuǎn)動(dòng)引起的周期性的沖擊。從圖8中可以看出，在時(shí)間區(qū)間［0.3 s 0.4 s］和［0.7 s 0.8 s］之間有較強(qiáng)的沖擊產(chǎn)生，屬于故障區(qū)間。為了找出其中的故障信號(hào)以及定位故障區(qū)間，使用本文提出的多尺度直線擬合法進(jìn)行查找。查找出的信號(hào)突變點(diǎn)在圖8中標(biāo)出。從查找結(jié)果可以看出，該方法能較好地定位到故障區(qū)間。但是由于噪聲信號(hào)干擾較大及信號(hào)較復(fù)雜，導(dǎo)致某些區(qū)間故障點(diǎn)的定位不是特別精準(zhǔn)。

圖8 行星齒輪箱振動(dòng)信號(hào)及檢測結(jié)果Fig.8 Vibration signal of planetary gear and detection result

針對(duì)傅里葉變換后的頻譜信號(hào)，也可以用該檢測方法進(jìn)行頻譜中最大分量的檢測。對(duì)于兩個(gè)正弦信號(hào)的疊加信號(hào)x，x=0.6 ×sin(2π×50t)+0.35 ×sin(2π ×120t)，再疊加一個(gè)比較大的噪聲信號(hào)，其信噪比為-9.269 9 dB. 對(duì)添加噪聲以后的信號(hào)進(jìn)行快速傅里葉變換(FFT)，得到其頻譜圖，如圖9所示。從圖中可以看出，經(jīng)過FFT 后，能很明顯看到在50 Hz 和120 Hz 有較大的頻率分量。使用本文的“多尺度直線擬合法”進(jìn)行突變點(diǎn)的檢測，從檢測結(jié)果來看，該方法可以很準(zhǔn)確地檢測到頻譜中最大的頻率分量。

圖9 頻譜圖及檢測結(jié)果Fig.9 Frequency spectrogram and detection result

4 結(jié)論

多尺度直線擬合法檢測突變點(diǎn)其本質(zhì)就是通過尺度變換，并不斷縮小檢測范圍，利用擬合線段斜率的變化從而最終逼近時(shí)間序列中的突變點(diǎn)。而初始尺度的選擇對(duì)檢測結(jié)果有較大影響，本文給出了一種有效的初始尺度計(jì)算方法。將多尺度直線擬合法同目前應(yīng)用較多的其他3 種方法進(jìn)行比較，證明了該方法的計(jì)算效率、準(zhǔn)確度以及抗干擾方面都具有一定的優(yōu)越性。將多尺度直線擬合檢測方法應(yīng)用于實(shí)際的實(shí)驗(yàn)系統(tǒng)中，也取得了很好的效果，提供了檢測效率和檢測精度。對(duì)于初始尺度的計(jì)算方法還可以再進(jìn)一步進(jìn)行研究。

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］王金花，張榮剛，張攀，等.兩種水沙系列突變點(diǎn)算法的對(duì)比分析［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］劉嘉琦，龔政，張長寬.重大水利工程運(yùn)用對(duì)長江入海徑流量的影響［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］陳遠(yuǎn)中，陸寶宏，張育德，等.改進(jìn)的有序聚類分析法提取時(shí)間序列轉(zhuǎn)折點(diǎn)［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］魏鳳英.現(xiàn)代氣候統(tǒng)計(jì)診斷與預(yù)測技術(shù)［M］. 北京:氣象出版社，2009:62 -73.WEI Feng-ying. Modern climatological statistical diagnosis and prediction technology［M］. Beijing:Weteorological Press，2009:62 -73.(in Chinese)

［9］吳凡，熊高君，葉志嬋.小波變換在信號(hào)突變點(diǎn)檢測中的應(yīng)用［J］.計(jì)算機(jī)與現(xiàn)代化，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］賴富文，張志杰，劉景江，等. 基于炮口脈沖噪聲信號(hào)的射速測試方法研究［J］. 兵工學(xué)報(bào)，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)