黃雄波
?。ǚ鹕铰殬I(yè)技術學院 電子信息系，廣東 佛山 528137）
摘要：由于自相關函數(shù)刻畫了時序數(shù)據在不同時刻取值的線性相關程度，故其在時序數(shù)據的統(tǒng)計分析中得到了廣泛的應用。討論了基于FFT變換的自相關函數(shù)計算原理，結合非平穩(wěn)時序數(shù)據的辨識需求，基于自相關函數(shù)理論對趨勢和周期成份的分離次序以及殘留序列的隨機類型識別等問題進行了深入分析，進一步提出了一種改進的非平穩(wěn)時序數(shù)據的辨識算法。實驗驗證了改進算法的合理性和有效性。
關鍵詞： 自相關函數(shù)； FFT變換；非平穩(wěn)時序數(shù)據；系統(tǒng)辨識
0引言
在工程、經濟、自然和社會科學等很多領域中，被考查對象在其歷史的演變過程中，其相關的物理量常常以一系列隨時間而變化的數(shù)據序列而被人們記載，這種序列通稱為時間序列或時序數(shù)據 ［13］。事實上，由于受到諸多偶然因素的影響，時序數(shù)據很難用一個精確的數(shù)學表達式來描述，但由于它們大都具有統(tǒng)計規(guī)律的特性，因此，可以通過對時序數(shù)據的辨識和建模，進而達到認識事物、掌握其內在變化規(guī)律的目的［46］。
一般地說，事物在演變過程中往往具有某種趨勢或周期規(guī)律的特性，故在現(xiàn)實中所獲得的時序數(shù)據也就具有非平穩(wěn)的特點，即序列的均值(一階矩)為非常數(shù)且自相關函數(shù)(二階矩)與起始時間有關［78］。由于非平穩(wěn)時序數(shù)據具有時變的統(tǒng)計結構，故其辨識和建模的過程比平穩(wěn)時序數(shù)據要復雜得多，目前比較有效的做法是［913］：首先從原始序列中分離出趨勢和周期成份并分別對它們進行辨識，然后對殘留的平穩(wěn)隨機序列建立相應的AR、MA或ARMA模型。本文基于自相關函數(shù)的理論，對諸如趨勢和周期成份的分離次序以及殘留序列的隨機類型識別等問題進行了分析，在此基礎上，提出一種改進的非平穩(wěn)時序數(shù)據的辨識算法，實驗證明改進算法是有效的。
1問題描述
1.1時序數(shù)據的自相關函數(shù)
設樣本長度為n的非平穩(wěn)時序數(shù)據{y(t),t=0,1,…,n－1}在相鄰p時刻的數(shù)據取值分別為yt和yt+p，則yt的p階樣本自相關函數(shù)的估算值r^(p)可用式(1)進行計算：
自相關函數(shù)是時序數(shù)據的二階統(tǒng)計特性，它表征了序列數(shù)據項之間的依賴程度及趨勢的變化情況，據此，可以利用自相關函數(shù)的上述性質來識別序列的趨勢和周期成份的顯著性及它們之間的關聯(lián)信息等。
1.2基于FFT變換的自相關函數(shù)計算
利用式(1)計算r^(p)時，若序列的樣本長度n比較大，則所需要的運算量也隨之增加，為了提高運算速度，可以利用快速傅里葉變換(Fast Fourier Transformation，F(xiàn)FT)算法來實現(xiàn)對r^(p)的快速運算［14］。
對r^(p)施行Fourier變換，由式(1)有：
為了能用離散傅里葉變換(Discrete Fourier Transform，DFT)來計算式(2)的線性卷積運算，故需要對時序數(shù)據y(t)補充n位零數(shù)值，得y′(t)，即：
綜合式(2)~式(5)可知，基于FFT變換的自相關函數(shù)的計算流程為：(1)對y(t)補充n位零數(shù)值，得y′(t)；(2) 對y′(t)做DFT變換，得到y(tǒng)′(m)；(3) 對1n|y′(m)|2做逆DFT變換，得到r^′(p)；(4)對r^′(p)中0~n的點進行整體右移n位，得到r^(p)。
2非平穩(wěn)時序數(shù)據的辨識改進
2.1自相關函數(shù)在非平穩(wěn)時序數(shù)據的辨識應用
非平穩(wěn)時序數(shù)據y(t)的辨識模型如式(6)所示:
其中，∑ui=0αiti是趨勢成份，由一個u階的多項式進行辨識，αi為多項式的各階次的辨識參數(shù)；∑vi=1λisin(iωt)是周期成份，ω為基波角頻率，λi為各諧波的幅值；x(t)為殘留的平穩(wěn)隨機序列。
2.1.1趨勢和周期成份的分離
對于既有趨勢又有周期影響的非平穩(wěn)時序數(shù)據而言，當趨勢成份很強時，則其趨勢特點表現(xiàn)突出，甚至會把周期特點淹沒掉；而當周期成份很強，其趨勢特點也有可能顯現(xiàn)不出來。 據此，有必要在現(xiàn)有的非平穩(wěn)時序數(shù)據辨識算法的基礎上，根據趨勢和周期成份在序列中的輕重地位進行依次分離和辨識。
對于趨勢成份明顯的時序數(shù)據，由于緊鄰的數(shù)據項具有相同的取值方向，故其自相關函數(shù)的前3~4階通常會出現(xiàn)正數(shù)且具有大于2n的特點；而當序列存在周期長度為L的數(shù)據成份時，則由于各個周期內相關位置的數(shù)據具有同時為大或同時為小的特點，故其L,2L,…,int(nL)L(int()為取整函數(shù))等各階次的自相關函數(shù)均為正值。依照上述規(guī)則，可以通過計算各階次的自相關函數(shù)來判別待辨識的非平穩(wěn)時序數(shù)據是否存在著趨勢和周期成份，當僅發(fā)現(xiàn)趨勢時，則首先分離趨勢成份，然后再對殘留序列進行周期識別；反之，當僅發(fā)現(xiàn)周期時，則首先分離周期成份，然后再對殘留序列進行趨勢識別；而當趨勢和周期同時顯著時，則應從原始序列中同時分離這兩種數(shù)據成份。
以趨勢和周期成份同時顯著為例，其基于最小二乘法的分離過程如下：
為滿足arg mine=∑nt=0|y(t)－f(t)|2，依最小二乘法有：
求解式(7)的線性方程組［15］，便可得到趨勢和周期成份同時顯著時的辨識參數(shù)向量［α0,α1,…,αu,λ1,λ2,…,λv］T。
2.1.2殘留序列的平穩(wěn)性判別及辨識
通常，非平穩(wěn)時序數(shù)據y(t)在去掉趨勢和周期成份的影響后，其殘留序列x(t)可能是一個分段局部平穩(wěn)、整體平穩(wěn)或完全隨機型的序列，故在殘留序列x(t)的辨識之前還應進行隨機類型的識別。完全隨機型序列其各階次的自相關函數(shù)具有零值的特點；而分段局部平穩(wěn)或整體平穩(wěn)的隨機型序列在平穩(wěn)范圍內其數(shù)據項之間是允許存在自相關的，但數(shù)據項xt和xt+p的相關性會隨著間隔項p的增大而減少。在工程應用中，若某序列的各階次的自相關函數(shù)滿足式(8)，則可把該序列視為平穩(wěn)時序數(shù)據［16］。
式(8)中，符號“&&”為邏輯與運算，“int( )”為取整函數(shù)，該邏輯表達式描述了序列的自相關函數(shù)的絕對值在1~2階次時均大于2n，但其后便衰減至－2n,2n區(qū)間內。
由于完全隨機型序列不包含任何形式的模型，故無需對它進行任何的處理；而平穩(wěn)隨機序列則可以采用式(9)的AR模型進行建模：
x(t)=φ1xt－1+φ2xt－2+…+φkxt－k+un(9)
式(9)中，［φ1,φ2,…,φk］T為自回歸參數(shù)向量，un是均值為零、方差為σ2的正態(tài)分布白噪聲。自回歸參數(shù)向量［φ1,φ2,…,φk］T可用式(10)的DurbinLevinson遞推公式求得，而遞推的具體計算流程則如圖1所示。
對于分段局部平穩(wěn)的殘留序列而言，在辨識之前還應完成各段平穩(wěn)子序列的劃分，即需要確定平穩(wěn)子序列的段數(shù)及各段之間的分割點。根據平穩(wěn)時序數(shù)據的自相關函數(shù)特征，可得到如下平穩(wěn)子序列的分割方法：以1個數(shù)據項為步長自左至右地搜索分割點，若新增數(shù)據項后的子序列其自相關函數(shù)仍滿足式(8)中的邏輯表達式，則合并該新增數(shù)據項并繼續(xù)向右搜索；否則以該新增數(shù)據項為分割點，新建另一平穩(wěn)子序列并繼續(xù)向右搜索；重復上述過程直至遍歷整個殘留序列為止。完成平穩(wěn)子序列的劃分后，便可根據式(10)中的遞推公式求得它們各自的辨識參數(shù)向量。
2.2非平穩(wěn)時序數(shù)據的辨識改進算法
綜上所述，可設計如下的非平穩(wěn)時序數(shù)據的辨識改進算法。
算法：基于自相關函數(shù)的非平穩(wěn)時序數(shù)據的辨識算法。
輸入：樣本長度為n的非平穩(wěn)時間數(shù)據y(t)。
輸出：趨勢成份的辨識參數(shù)向量［α0,α1,…,αu］T，周期成份的辨識參數(shù)向量［λ1,λ2,…,λv］T，隨機成份的辨識參數(shù)向量［φ1,φ2,…,φk］T。
步驟1設計相關的公用函數(shù)模塊。
模塊1依照上述的式(2)～式(5)編寫序列的FFT自相關函數(shù)計算模塊。
模塊2編寫趨勢成份顯著的識別模塊，具體的識別流程為：
if((r^(1)>(2/sqrt(n)))&&(r^(2)>(2/sqrt(n)))&&(r^(3)>(2/sqrt(n))))
{
標識趨勢成份顯著；
}
模塊3編寫周期成份顯著的標識模塊，具體的識別流程為：
for(L=2;L<=int(2/n);L++)
{//遍歷搜索顯著的周期成份
if((r^(L)>0)&&(r^(2*L)>0)&&…&&(r^(int(n/L)*L)>0))
{
標識周期成份顯著；
}
}
模塊4線性方程組(系數(shù)矩陣A=［aij］∈Rn×n)的GaussSeidel迭代求解模塊，其迭代求解公式為：
步驟2從非平穩(wěn)時序數(shù)據中分離趨勢和周期成份。
(1)調用公用模塊1計算y(t)的r^(1)、r^(2)、…、r^(2n)等各階次的自相關函數(shù)；
(2)調用公用模塊2和模塊3進行趨勢和周期成份的識別；
(3)對(2)的識別結果進行處理，分“僅趨勢成份顯著”、“僅周期成份顯著”、“趨勢和周期成份同時顯著”、“趨勢和周期均不顯著”4種情況；
case1://限于篇幅，這里只給出“僅趨勢成份顯著”的處理流程
{
(1)基于最小二乘法計算趨勢成份顯著時y(t)所對應的線性方程組的系數(shù)矩陣A；
(2)調用公用模塊4求得趨勢成份的辨識參數(shù)向量［α0,α1,…,αu］T；
(3)計算殘留序列，x(t)=y(t)－∑ui=0αiti－1；
(4)調用公用模塊1計算殘留序列x(t)的r^(1)、r^(2)、…、r^(2n)等各階次的自相關函數(shù)；
(5)調用公用模塊3對殘留序列x(t)進行周期成份的識別；
(6)對殘留序列x(t)的周期成份識別結果進行處理：
if(殘留序列x(t)的周期成份顯著)
{
①基于最小二乘法計算周期成份顯著時x(t)所對應的線性方程組的系數(shù)矩陣A；
?、谡{用公用模塊4求得周期成份的辨識參數(shù)向量［λ1,λ2,…,λv］T；
?、鄹職埩粜蛄?，x(t)=x(t)－∑vi=1λi(siniωt)；
}
?、苻D步驟3；
}
步驟3對殘留序列x(t)進行辨識。
(1)遍歷殘留序列x(t)，完成各局部平穩(wěn)子序列的劃分；
for(i=0;i<n;i++)
{
sl=i－d［m］;
//計算當前子序列的序列長度，d［m］為第m段平穩(wěn)子序列的起點位置
if(sl>=min)
{
if((r^(1)==0)&&(r^(2)==0)&&…&&(r^(int(2/sl))==0))
{
標識當前子序列為完全隨機型序列；
}
else
{
if((abs(r^(1))>(2/sqrt(sl)))&&(abs(r^(2))>(2/sqrt(sl)))&&…&&(abs(r^(3))<(2/sqrt(sl)))&&…&&(abs(r^(int(2/n)))<(2/sqrt(sl))))
{
標識當前的第m段子序列為局部平穩(wěn)子序列；
}
else//第i項數(shù)據為分割點
{
m++;//局部平穩(wěn)子序列的段數(shù)增1；
d［m］=i;//保存第m段局部平穩(wěn)子序列的分割點；
(2)對平穩(wěn)序列的辨識參數(shù)向量進行估計：
for(i=0;i<m;i++)//m為(1)確定的局部平穩(wěn)子序列的段數(shù)
{
while(abs(φi_kk)>(2/sqrt(n)))
//取顯著性水平α=0.05，若|φi_kk|<2n成立，則認為k階自相關函數(shù)等于0；
{
利用式(10)對第i段局部平穩(wěn)子序列的辨識參數(shù)向量［φi_1，φi_2，…，φi_k］T進行遞推計算；
}
步驟4輸出趨勢成份、周期成份及隨機成份的辨識參數(shù)向量，算法結束。
3實驗及結果分析
為了驗證上述算法的合理性及有效性，這里將分別對趨勢成份顯著、周期成份顯著、殘留序列為分段局部平穩(wěn)的非平穩(wěn)時序數(shù)據進行相關的辨識實驗。實驗的硬件環(huán)境為惠普ProDesk 490 G2 MT商用臺式機(CPU:i5-45704×3.2 GHz;內存:4 GB DDR3 1600)，軟件環(huán)境及開發(fā)工具為Windows 8.1+Microsoft Visual C++2010。實驗的主要目的是考察改進算法與現(xiàn)有算法之間的辨識精度及計算效能。
3.1實驗設計
不失一般性，實驗所用的非平穩(wěn)時序數(shù)據約定滿足如下假設：(1)趨勢成份為二次多項式α0+α1t+α2t2,且α1α0和α1α2成立，即以直線趨勢成份為主；(2)周期成份為一次諧波λ1sinωt；(3)序列的樣本長度n=200。
分別利用改進算法和現(xiàn)有算法進行如下實驗：
實驗1：趨勢成份顯著的非平穩(wěn)時序數(shù)據的辨識，所選的數(shù)據模型如式(11):
y(t)=3t+0.011t2+0.3sin0.785t－0.28y(t－1)(11)
實驗2：周期成份顯著的非平穩(wěn)時序數(shù)據的辨識，所選的數(shù)據模型如式（12）：
y(t)=0.4t+0.001t2+200sin0.785t－0.28y(t－1)(12)
實驗3：殘留序列為分段局部平穩(wěn)的非平穩(wěn)時序數(shù)據的辨識，所選的數(shù)據模型如式（13）：
3.2實驗結果與討論
為了能全面考察改進算法與現(xiàn)有算法的辨識精度，這里引入了均方標準誤差(Root Mean Squared Error，RMSE)和平均絕對百分誤差(Mean Absolute Percentage Error，MAPE)共兩個性能評價標準，具體定義如式(14)和式(15)所示。
式(14)和式(15)中，y(t)為原始序列，y^(t)為辨識序列，n為序列的樣本長度。
實驗1~實驗3的辨識結果如表1、表2及圖2~圖4所示。其中，表1為辨識序列的函數(shù)表達式，表2為辨識性能的具體數(shù)值，而圖2~圖4則是各次實驗的辨識擬合曲線。
從表1的辨識函數(shù)表達式易知，對于趨勢或周期成份顯著的非平穩(wěn)時序數(shù)據而言，由于現(xiàn)有算法沒有按照一定的次序進行分離，故不顯著的數(shù)據成份未能有效地被識別，進而導致了殘留隨機序列存在一定的偏差。在表2中，根據辨識性能的具體數(shù)值可發(fā)現(xiàn)，改進算法在增加了一定的計算耗時后其辨識精度有了顯著的提高。從圖4的辨識擬合曲線則不難發(fā)現(xiàn)，現(xiàn)有算法對分段局部平穩(wěn)序列的辨識效果較差，究其原因是因為用了平穩(wěn)隨機模型來對整體不平穩(wěn)的序列進行辨識，故不僅存在較大的辨識誤差而且辨識序列的模型階數(shù)也出現(xiàn)了增加；而改進算法雖然花費了相當?shù)挠嬎愫臅r，但能從根本上保證辨識精度。綜上所述，本文提出的改進算法對非平穩(wěn)時序數(shù)據具有良好的辨識性能。
4結論
本文提出了一種非平穩(wěn)時序數(shù)據的辨識改進算法，由于對趨勢和周期成份的分離次序以及殘留序列的隨機類型識別等細節(jié)問題均作了相應的處理，故改進算法的辨識性能有了明顯的提升。下一步的主要工作是引入自相關函數(shù)的并行快速變換運算，同時研究更為有效的平穩(wěn)子序列劃分方法，以便進一步提升改進算法的計算效能。
參考文獻
?。?］ (德)蓋哈德·克西蓋斯納，(德)約根·沃特斯，(德)烏沃·哈斯勒.現(xiàn)代時間序列分析導論［M］.張延群，劉曉飛，譯.北京: 中國人民大學出版社,2015．
?。?］ (美)漢密爾頓.時間序列分析［M］．夏曉華,譯.北京:中國人民大學出版社，2015．
?。?］ 張鵬.飛行數(shù)據的時間序列分析方法及其應用［M］.北京:國防工業(yè)出版社,2013．
?。?］ 趙雪巖，李衛(wèi)華.系統(tǒng)建模與仿真［M］.北京:國防工業(yè)出版社,2015．
?。?］ 蕭德云.系統(tǒng)辨識理論及應用［M］.北京:清華大學出版社,2014．
?。?］ 黃介武.線性與廣義線性模型中參數(shù)估計的一些研究［D］.重慶: 重慶大學,2014．
?。?］ 楊繼生.非平穩(wěn)綜列數(shù)據分析：理論與應用［M］.北京:中國社會科學出版社,2010．
?。?］ 王宏禹，邱天爽，陳哲.非平穩(wěn)隨機信號分析與處理(第2版)［M］．北京:國防工業(yè)出版社,2008．
［9］ 黃雄波.時序數(shù)據趨勢項的分段擬合［J］．計算機系統(tǒng)應用，2015，24(2):174179．
?。?0］ 黃雄波.多周期時序數(shù)據的傅氏級數(shù)擬合算法［J］．計算機系統(tǒng)應用，2015，24(7):142148．
［11］ 黃雄波.時序數(shù)據的周期模式發(fā)現(xiàn)算法的遞推改進［J］．計算機技術與發(fā)展，2016，26(2):4751．
?。?2］ MAHESWARAN R, KHOSA R.Wavelet volterra coupled models for forecasting of nonlinear and nonstationary time series［J］．Neurocomputing, 2015, 149(Part B, 3):10741084.
?。?3］ SPIRIDONAKOS M D, FASSOIS S D. Nonstationary random vibration modelling and analysis via functional series timedependent ARMA (FSTARMA) modelsA critical survey［J］．Mechanical Systems and Signal Processing, 2014, 47(12):175224.