基于自相關(guān)函數(shù)的非平穩(wěn)時序數(shù)據(jù)的辨識改進(jìn)*

黃雄波

(佛山職業(yè)技術(shù)學(xué)院 電子信息系，廣東 佛山 528137)

?

基于自相關(guān)函數(shù)的非平穩(wěn)時序數(shù)據(jù)的辨識改進(jìn)*

黃雄波

(佛山職業(yè)技術(shù)學(xué)院 電子信息系，廣東 佛山 528137)

摘要：由于自相關(guān)函數(shù)刻畫了時序數(shù)據(jù)在不同時刻取值的線性相關(guān)程度，故其在時序數(shù)據(jù)的統(tǒng)計分析中得到了廣泛的應(yīng)用。討論了基于FFT變換的自相關(guān)函數(shù)計算原理，結(jié)合非平穩(wěn)時序數(shù)據(jù)的辨識需求，基于自相關(guān)函數(shù)理論對趨勢和周期成份的分離次序以及殘留序列的隨機(jī)類型識別等問題進(jìn)行了深入分析，進(jìn)一步提出了一種改進(jìn)的非平穩(wěn)時序數(shù)據(jù)的辨識算法。實(shí)驗(yàn)驗(yàn)證了改進(jìn)算法的合理性和有效性。

關(guān)鍵詞：自相關(guān)函數(shù)； FFT變換；非平穩(wěn)時序數(shù)據(jù)；系統(tǒng)辨識

引用格式：黃雄波. 基于自相關(guān)函數(shù)的非平穩(wěn)時序數(shù)據(jù)的辨識改進(jìn)[J].微型機(jī)與應(yīng)用，2016,35(13)：10-14，18.

0引言

在工程、經(jīng)濟(jì)、自然和社會科學(xué)等很多領(lǐng)域中，被考查對象在其歷史的演變過程中，其相關(guān)的物理量常常以一系列隨時間而變化的數(shù)據(jù)序列而被人們記載，這種序列通稱為時間序列或時序數(shù)據(jù)[1-3]。事實(shí)上，由于受到諸多偶然因素的影響，時序數(shù)據(jù)很難用一個精確的數(shù)學(xué)表達(dá)式來描述，但由于它們大都具有統(tǒng)計規(guī)律的特性，因此，可以通過對時序數(shù)據(jù)的辨識和建模，進(jìn)而達(dá)到認(rèn)識事物、掌握其內(nèi)在變化規(guī)律的目的[4-6]。

一般地說，事物在演變過程中往往具有某種趨勢或周期規(guī)律的特性，故在現(xiàn)實(shí)中所獲得的時序數(shù)據(jù)也就具有非平穩(wěn)的特點(diǎn)，即序列的均值(一階矩)為非常數(shù)且自相關(guān)函數(shù)(二階矩)與起始時間有關(guān)[7-8]。由于非平穩(wěn)時序數(shù)據(jù)具有時變的統(tǒng)計結(jié)構(gòu)，故其辨識和建模的過程比平穩(wěn)時序數(shù)據(jù)要復(fù)雜得多，目前比較有效的做法是[9-13]：首先從原始序列中分離出趨勢和周期成份并分別對它們進(jìn)行辨識，然后對殘留的平穩(wěn)隨機(jī)序列建立相應(yīng)的AR、MA或ARMA模型。本文基于自相關(guān)函數(shù)的理論，對諸如趨勢和周期成份的分離次序以及殘留序列的隨機(jī)類型識別等問題進(jìn)行了分析，在此基礎(chǔ)上，提出一種改進(jìn)的非平穩(wěn)時序數(shù)據(jù)的辨識算法，實(shí)驗(yàn)證明改進(jìn)算法是有效的。

1問題描述

1.1時序數(shù)據(jù)的自相關(guān)函數(shù)

(1)

自相關(guān)函數(shù)是時序數(shù)據(jù)的二階統(tǒng)計特性，它表征了序列數(shù)據(jù)項(xiàng)之間的依賴程度及趨勢的變化情況，據(jù)此，可以利用自相關(guān)函數(shù)的上述性質(zhì)來識別序列的趨勢和周期成份的顯著性及它們之間的關(guān)聯(lián)信息等。

1.2基于FFT變換的自相關(guān)函數(shù)計算

(2)

為了能用離散傅里葉變換(Discrete Fourier Transform，DFT)來計算式(2)的線性卷積運(yùn)算，故需要對時序數(shù)據(jù)y(t)補(bǔ)充n位零數(shù)值，得y′(t)，即：

(3)

設(shè)y′(t)的Fourier變換為y′(ejω)，則由于

(4)

故有

(5)

2非平穩(wěn)時序數(shù)據(jù)的辨識改進(jìn)

2.1自相關(guān)函數(shù)在非平穩(wěn)時序數(shù)據(jù)的辨識應(yīng)用

非平穩(wěn)時序數(shù)據(jù)y(t)的辨識模型如式(6)所示:

(6)

2.1.1趨勢和周期成份的分離

對于既有趨勢又有周期影響的非平穩(wěn)時序數(shù)據(jù)而言，當(dāng)趨勢成份很強(qiáng)時，則其趨勢特點(diǎn)表現(xiàn)突出，甚至?xí)阎芷谔攸c(diǎn)淹沒掉；而當(dāng)周期成份很強(qiáng)，其趨勢特點(diǎn)也有可能顯現(xiàn)不出來。 據(jù)此，有必要在現(xiàn)有的非平穩(wěn)時序數(shù)據(jù)辨識算法的基礎(chǔ)上，根據(jù)趨勢和周期成份在序列中的輕重地位進(jìn)行依次分離和辨識。

以趨勢和周期成份同時顯著為例，其基于最小二乘法的分離過程如下：

(7)

求解式(7)的線性方程組[15]，便可得到趨勢和周期成份同時顯著時的辨識參數(shù)向量[α0,α1,…,αu,λ1,λ2,…,λv]T。

2.1.2殘留序列的平穩(wěn)性判別及辨識

通常，非平穩(wěn)時序數(shù)據(jù)y(t)在去掉趨勢和周期成份的影響后，其殘留序列x(t)可能是一個分段局部平穩(wěn)、整體平穩(wěn)或完全隨機(jī)型的序列，故在殘留序列x(t)的辨識之前還應(yīng)進(jìn)行隨機(jī)類型的識別。完全隨機(jī)型序列其各階次的自相關(guān)函數(shù)具有零值的特點(diǎn)；而分段局部平穩(wěn)或整體平穩(wěn)的隨機(jī)型序列在平穩(wěn)范圍內(nèi)其數(shù)據(jù)項(xiàng)之間是允許存在自相關(guān)的，但數(shù)據(jù)項(xiàng)xt和xt+p的相關(guān)性會隨著間隔項(xiàng)p的增大而減少。在工程應(yīng)用中，若某序列的各階次的自相關(guān)函數(shù)滿足式(8)，則可把該序列視為平穩(wěn)時序數(shù)據(jù)[16]。

(8)

由于完全隨機(jī)型序列不包含任何形式的模型，故無需對它進(jìn)行任何的處理；而平穩(wěn)隨機(jī)序列則可以采用式(9)的AR模型進(jìn)行建模：

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)的Durbin-Levinson遞推公式求得，而遞推的具體計算流程則如圖1所示。

(10)

圖1　自回歸參數(shù)的Durbin-Levinson遞推計算流程

對于分段局部平穩(wěn)的殘留序列而言，在辨識之前還應(yīng)完成各段平穩(wěn)子序列的劃分，即需要確定平穩(wěn)子序列的段數(shù)及各段之間的分割點(diǎn)。根據(jù)平穩(wěn)時序數(shù)據(jù)的自相關(guān)函數(shù)特征，可得到如下平穩(wěn)子序列的分割方法：以1個數(shù)據(jù)項(xiàng)為步長自左至右地搜索分割點(diǎn)，若新增數(shù)據(jù)項(xiàng)后的子序列其自相關(guān)函數(shù)仍滿足式(8)中的邏輯表達(dá)式，則合并該新增數(shù)據(jù)項(xiàng)并繼續(xù)向右搜索；否則以該新增數(shù)據(jù)項(xiàng)為分割點(diǎn)，新建另一平穩(wěn)子序列并繼續(xù)向右搜索；重復(fù)上述過程直至遍歷整個殘留序列為止。完成平穩(wěn)子序列的劃分后，便可根據(jù)式(10)中的遞推公式求得它們各自的辨識參數(shù)向量。

2.2非平穩(wěn)時序數(shù)據(jù)的辨識改進(jìn)算法

綜上所述，可設(shè)計如下的非平穩(wěn)時序數(shù)據(jù)的辨識改進(jìn)算法。

算法：基于自相關(guān)函數(shù)的非平穩(wěn)時序數(shù)據(jù)的辨識算法。

輸入：樣本長度為n的非平穩(wěn)時間數(shù)據(jù)y(t)。

輸出：趨勢成份的辨識參數(shù)向量[α0,α1,…,αu]T，周期成份的辨識參數(shù)向量[λ1,λ2,…,λv]T，隨機(jī)成份的辨識參數(shù)向量[φ1,φ2,…,φk]T。

步驟1設(shè)計相關(guān)的公用函數(shù)模塊。

模塊1依照上述的式(2)～式(5)編寫序列的FFT自相關(guān)函數(shù)計算模塊。

模塊2編寫趨勢成份顯著的識別模塊，具體的識別流程為：

{

標(biāo)識趨勢成份顯著；

}

模塊3編寫周期成份顯著的標(biāo)識模塊，具體的識別流程為：

for(L=2;L<=int(2/n);L++)

{//遍歷搜索顯著的周期成份

{

標(biāo)識周期成份顯著；

}

}

模塊4線性方程組(系數(shù)矩陣A=[aij]∈Rn×n)的Gauss-Seidel迭代求解模塊，其迭代求解公式為：

步驟2從非平穩(wěn)時序數(shù)據(jù)中分離趨勢和周期成份。

(2)調(diào)用公用模塊2和模塊3進(jìn)行趨勢和周期成份的識別；

(3)對(2)的識別結(jié)果進(jìn)行處理，分“僅趨勢成份顯著”、“僅周期成份顯著”、“趨勢和周期成份同時顯著”、“趨勢和周期均不顯著”4種情況；

case1://限于篇幅，這里只給出“僅趨勢成份顯著”的處理流程

{

(1)基于最小二乘法計算趨勢成份顯著時y(t)所對應(yīng)的線性方程組的系數(shù)矩陣A；

(2)調(diào)用公用模塊4求得趨勢成份的辨識參數(shù)向量[α0,α1,…,αu]T；

(5)調(diào)用公用模塊3對殘留序列x(t)進(jìn)行周期成份的識別；

(6)對殘留序列x(t)的周期成份識別結(jié)果進(jìn)行處理：

if(殘留序列x(t)的周期成份顯著)

{

①基于最小二乘法計算周期成份顯著時x(t)所對應(yīng)的線性方程組的系數(shù)矩陣A；

②調(diào)用公用模塊4求得周期成份的辨識參數(shù)向量[λ1,λ2,…,λv]T；

}

④轉(zhuǎn)步驟3；

}

步驟3對殘留序列x(t)進(jìn)行辨識。

(1)遍歷殘留序列x(t)，完成各局部平穩(wěn)子序列的劃分；

for(i=0;i

{

sl=i-d[m];

//計算當(dāng)前子序列的序列長度，d[m]為第m段平穩(wěn)子序列的起點(diǎn)位置

if(sl>=min)

{

{

標(biāo)識當(dāng)前子序列為完全隨機(jī)型序列；

}

else

{

{

標(biāo)識當(dāng)前的第m段子序列為局部平穩(wěn)子序列；

}

else//第i項(xiàng)數(shù)據(jù)為分割點(diǎn)

{

m++;//局部平穩(wěn)子序列的段數(shù)增1；

d[m]=i;//保存第m段局部平穩(wěn)子序列的分割點(diǎn)；

表1　各次實(shí)驗(yàn)所得的辨識函數(shù)表達(dá)式

}

}

}

}

(2)對平穩(wěn)序列的辨識參數(shù)向量進(jìn)行估計：

for(i=0;i

{

while(abs(φi_kk)>(2/sqrt(n)))

{

利用式(10)對第i段局部平穩(wěn)子序列的辨識參數(shù)向量[φi_1，φi_2，…，φi_k]T進(jìn)行遞推計算；

}

步驟4輸出趨勢成份、周期成份及隨機(jī)成份的辨識參數(shù)向量，算法結(jié)束。

3實(shí)驗(yàn)及結(jié)果分析

為了驗(yàn)證上述算法的合理性及有效性，這里將分別對趨勢成份顯著、周期成份顯著、殘留序列為分段局部平穩(wěn)的非平穩(wěn)時序數(shù)據(jù)進(jìn)行相關(guān)的辨識實(shí)驗(yàn)。實(shí)驗(yàn)的硬件環(huán)境為惠普ProDesk 490 G2 MT商用臺式機(jī)(CPU:i5-45704×3.2 GHz;內(nèi)存:4 GB DDR3 1600)，軟件環(huán)境及開發(fā)工具為Windows 8.1+Microsoft Visual C++2010。實(shí)驗(yàn)的主要目的是考察改進(jìn)算法與現(xiàn)有算法之間的辨識精度及計算效能。

3.1實(shí)驗(yàn)設(shè)計

不失一般性，實(shí)驗(yàn)所用的非平穩(wěn)時序數(shù)據(jù)約定滿足如下假設(shè)：(1)趨勢成份為二次多項(xiàng)式α0+α1t+α2t2,且α1?α0和α1?α2成立，即以直線趨勢成份為主；(2)周期成份為一次諧波λ1sinωt；(3)序列的樣本長度n=200。

分別利用改進(jìn)算法和現(xiàn)有算法進(jìn)行如下實(shí)驗(yàn)：

實(shí)驗(yàn)1：趨勢成份顯著的非平穩(wěn)時序數(shù)據(jù)的辨識，所選的數(shù)據(jù)模型如式(11):

y(t)=3t+0.011t2+0.3sin0.785t-0.28y(t-1)

(11)

實(shí)驗(yàn)2：周期成份顯著的非平穩(wěn)時序數(shù)據(jù)的辨識，所選的數(shù)據(jù)模型如式(12)：

y(t)=0.4t+0.001t2+200sin0.785t-0.28y(t-1)

(12)

實(shí)驗(yàn)3：殘留序列為分段局部平穩(wěn)的非平穩(wěn)時序數(shù)據(jù)的辨識，所選的數(shù)據(jù)模型如式(13)：

(13)

3.2實(shí)驗(yàn)結(jié)果與討論

為了能全面考察改進(jìn)算法與現(xiàn)有算法的辨識精度，這里引入了均方標(biāo)準(zhǔn)誤差(Root Mean Squared Error，RMSE)和平均絕對百分誤差(Mean Absolute Percentage Error，MAPE)共兩個性能評價標(biāo)準(zhǔn)，具體定義如式(14)和式(15)所示。

(14)

(15)

實(shí)驗(yàn)1～實(shí)驗(yàn)3的辨識結(jié)果如表1、表2及圖2～圖4所示。其中，表1為辨識序列的函數(shù)表達(dá)式，表2為辨識性能的具體數(shù)值，而圖2～圖4則是各次實(shí)驗(yàn)的辨識擬合曲線。

表2　各次實(shí)驗(yàn)的辨識性能的具體數(shù)值

圖2　實(shí)驗(yàn)1對應(yīng)的辨識擬合曲線

圖3　實(shí)驗(yàn)2對應(yīng)的辨識擬合曲線

圖4　實(shí)驗(yàn)3對應(yīng)的辨識擬合曲線

從表1的辨識函數(shù)表達(dá)式易知，對于趨勢或周期成份顯著的非平穩(wěn)時序數(shù)據(jù)而言，由于現(xiàn)有算法沒有按照一定的次序進(jìn)行分離，故不顯著的數(shù)據(jù)成份未能有效地被識別，進(jìn)而導(dǎo)致了殘留隨機(jī)序列存在一定的偏差。在表2中，根據(jù)辨識性能的具體數(shù)值可發(fā)現(xiàn)，改進(jìn)算法在增加了一定的計算耗時后其辨識精度有了顯著的提高。從圖4的辨識擬合曲線則不難發(fā)現(xiàn)，現(xiàn)有算法對分段局部平穩(wěn)序列的辨識效果較差，究其原因是因?yàn)橛昧似椒€(wěn)隨機(jī)模型來對整體不平穩(wěn)的序列進(jìn)行辨識，故不僅存在較大的辨識誤差而且辨識序列的模型階數(shù)也出現(xiàn)了增加；而改進(jìn)算法雖然花費(fèi)了相當(dāng)?shù)挠嬎愫臅r，但能從根本上保證辨識精度。綜上所述，本文提出的改進(jìn)算法對非平穩(wěn)時序數(shù)據(jù)具有良好的辨識性能。

4結(jié)論

本文提出了一種非平穩(wěn)時序數(shù)據(jù)的辨識改進(jìn)算法，由于對趨勢和周期成份的分離次序以及殘留序列的隨機(jī)類型識別等細(xì)節(jié)問題均作了相應(yīng)的處理，故改進(jìn)算法的辨識性能有了明顯的提升。下一步的主要工作是引入自相關(guān)函數(shù)的并行快速變換運(yùn)算，同時研究更為有效的平穩(wěn)子序列劃分方法，以便進(jìn)一步提升改進(jìn)算法的計算效能。

參考文獻(xiàn)

[1] (德)蓋哈德·克西蓋斯納，(德)約根·沃特斯，(德)烏沃·哈斯勒.現(xiàn)代時間序列分析導(dǎo)論[M].張延群，劉曉飛，譯.北京: 中國人民大學(xué)出版社,2015．

[2] (美)漢密爾頓.時間序列分析[M]．夏曉華,譯.北京:中國人民大學(xué)出版社，2015．

[3] 張鵬.飛行數(shù)據(jù)的時間序列分析方法及其應(yīng)用[M].北京:國防工業(yè)出版社,2013．

[4] 趙雪巖，李衛(wèi)華.系統(tǒng)建模與仿真[M].北京:國防工業(yè)出版社,2015．

[5] 蕭德云.系統(tǒng)辨識理論及應(yīng)用[M].北京:清華大學(xué)出版社,2014．

[6] 黃介武.線性與廣義線性模型中參數(shù)估計的一些研究[D].重慶: 重慶大學(xué),2014．

[7] 楊繼生.非平穩(wěn)綜列數(shù)據(jù)分析：理論與應(yīng)用[M].北京:中國社會科學(xué)出版社,2010．

[8] 王宏禹，邱天爽，陳哲.非平穩(wěn)隨機(jī)信號分析與處理(第2版)[M]．北京:國防工業(yè)出版社,2008．

[9] 黃雄波.時序數(shù)據(jù)趨勢項(xiàng)的分段擬合[J]．計算機(jī)系統(tǒng)應(yīng)用，2015，24(2):174-179．

[10] 黃雄波.多周期時序數(shù)據(jù)的傅氏級數(shù)擬合算法[J]．計算機(jī)系統(tǒng)應(yīng)用，2015，24(7):142-148．

[11] 黃雄波.時序數(shù)據(jù)的周期模式發(fā)現(xiàn)算法的遞推改進(jìn)[J]．計算機(jī)技術(shù)與發(fā)展，2016，26(2):47-51．

[12] MAHESWARAN R, KHOSA R.Wavelet volterra coupled models for forecasting of nonlinear and non-stationary time series[J]．Neurocomputing, 2015, 149(Part B, 3):1074-1084.

[13] SPIRIDONAKOS M D, FASSOIS S D. Non-stationary random vibration modelling and analysis via functional series time-dependent ARMA (FS-TARMA) models-A critical survey[J]．Mechanical Systems and Signal Processing, 2014, 47(1-2):175-224.

[14] (美) 奧本海姆.離散時間信號處理(第3版)[M]．黃建國，劉樹棠，張國海，譯.北京:電子工業(yè)出版社，2015．

[15] 李慶揚(yáng)，王能超，易大義.數(shù)值分析[M].北京:清華大學(xué)出版社,2010．

[16] 侯文超.經(jīng)濟(jì)預(yù)測——理論、方法及應(yīng)用[M].北京:商務(wù)印書館出版,1993．

*基金項(xiàng)目：廣東省科技計劃工業(yè)攻關(guān)項(xiàng)目(2011B010200031)；佛山職業(yè)技術(shù)學(xué)院校級重點(diǎn)科研項(xiàng)目(2015KY006)

中圖分類號：TP311

文獻(xiàn)標(biāo)識碼：A

DOI：10.19358/j.issn.1674- 7720.2016.13.004

(收稿日期：2016-03-14)

作者簡介：

黃雄波(1975-)，通信作者，男，碩士，副教授，CCF會員(E200011115M)，主要研究方向：時間序列分析及數(shù)字圖像處理。E-mail：xb-huang@hotmail.com。

Improved identification of non-stationary time series data based on autocorrelation function

Huang Xiongbo

(Department of Electronic and Information Engineering ,Foshan Professional Technical College, Foshan 528137，China)

Abstract:Because the autocorrelation function is used to describe the linear dependence of the time series data at different time, it has been widely used in the statistical analysis of the time series data. This paper discussed the calculation principle of FFT transform autocorrelation function, combined with identification requirements of non stationary time series data, based on autocorrelation function theory, analysed the trend and cycle component data separation order and residual sequence of random type recognition problem indepth. Furthermore, an improved non stationary time series data identification algorithm is proposed. The rationality and validity of the improved algorithm are verified by experiments.

Key words:autocorrelation function; Fast Fourier Transform (FFT); non-stationary time series data; system identification