多尺度直線擬合法在時間序列突變點檢測中的應用

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

(北京交通大學 機械與電子控制工程學院，北京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)