一種平穩(wěn)時(shí)序數(shù)據(jù)的高效辨識(shí)改進(jìn)算法*

黃雄波

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

1 引 言

物理系統(tǒng)或工作部件在實(shí)際運(yùn)行過程中，各項(xiàng)性能狀態(tài)參數(shù)可通過安裝在某些部位的傳感器及專用測(cè)試設(shè)備加以測(cè)量，且系統(tǒng)的工作狀況在某種程度上可通過這些狀態(tài)參數(shù)來表征?；跁r(shí)序數(shù)據(jù)驅(qū)動(dòng)的系統(tǒng)辨識(shí)，就是以采集的過程數(shù)據(jù)為基礎(chǔ)，直接對(duì)物理系統(tǒng)的各類可用數(shù)據(jù)進(jìn)行分析和處理，以此判斷系統(tǒng)在未來可能發(fā)生故障的概率，或估計(jì)系統(tǒng)失效以至到達(dá)壽命閥值的時(shí)間節(jié)點(diǎn)[1-4]。通常，時(shí)序數(shù)據(jù)在經(jīng)過相關(guān)的函數(shù)變換后，總能分解出趨勢(shì)項(xiàng)、隨機(jī)項(xiàng)及周期項(xiàng)共三種數(shù)據(jù)成分，據(jù)此，基于時(shí)序數(shù)據(jù)驅(qū)動(dòng)的系統(tǒng)辨識(shí)就是對(duì)既有序列進(jìn)行數(shù)據(jù)成分的分離及建模[5]。

時(shí)序數(shù)據(jù)的隨機(jī)成分蘊(yùn)藏了豐富的信息，在系統(tǒng)失效和故障診斷方面有著重要的應(yīng)用價(jià)值。對(duì)于線性平穩(wěn)的隨機(jī)序列而言，美國(guó)學(xué)者博克斯與英國(guó)學(xué)者詹金斯根據(jù)序列的自相關(guān)函數(shù)、偏相關(guān)函數(shù)或兩者是否同時(shí)出現(xiàn)拖尾等情形，分別提出了3 種可用的辨識(shí)模型：自回歸（Auto-regressive, AR）模型、滑動(dòng)平均（Moving average, MA）模型和自回歸滑動(dòng)平均（Auto-regressive and moving average, ARMA）模型[6]。利用上述辨識(shí)模型進(jìn)行系統(tǒng)建模，就是需要確定辨識(shí)模型的階次和估計(jì)辨識(shí)模型的參數(shù)。由于自回歸模型具有模型求解簡(jiǎn)捷的優(yōu)點(diǎn)，且其它兩種模型亦可用一高階自回歸模型進(jìn)行逼近，故該模型在實(shí)際應(yīng)用中得到了更為廣泛的應(yīng)用。為了進(jìn)一步提升現(xiàn)有的自回歸模型辨識(shí)算法的相關(guān)效能，提出一種更為高效的改進(jìn)算法，并就算法的設(shè)計(jì)原理及實(shí)現(xiàn)過程做更深一層的分析和討論。

2 問題描述

設(shè)平穩(wěn)時(shí)序數(shù)據(jù)序列為xt(t=1,2,…,n)，其p 階線性自回歸模型的數(shù)學(xué)描述如下式所示：

式(1)中，p∈N+，為自回歸模型的階數(shù)；φ1，φ2，…，φp為自回歸模型的參數(shù)；ut是數(shù)學(xué)期望為0、方差為σ2的隨機(jī)白噪聲干擾，與xt-1，xt-2，…，xt-p不相關(guān)。

式(1)中的自回歸模型的具體求解過程為：首先基于某一準(zhǔn)則確定模型的階數(shù)p，然后對(duì)模型參數(shù)φ1，φ2，…，φp進(jìn)行估計(jì)[7]。對(duì)模型的定階方法進(jìn)行歸納，大致可以分為兩類[8]：信息量準(zhǔn)則法和線性代數(shù)法。著名的信息量準(zhǔn)則法有平均信息（AIC）準(zhǔn)則、最終預(yù)測(cè)誤差（FPE）準(zhǔn)則、最小描述長(zhǎng)度（MDL）準(zhǔn)則、自回歸傳遞（CAT）準(zhǔn)則；而典型的線性代數(shù)定階法有行列式檢驗(yàn)法、Gram－Schmidt 正交法和奇異值分解法。當(dāng)確定模型的階數(shù)后，便可用Levinson－Durbin 或Burg 算法對(duì)模型參數(shù)進(jìn)行遞推估計(jì)。

值得注意的是，迄今為止，自回歸模型的定階問題仍然沒有得到根本性的解決，不同物理系統(tǒng)往往需要使用不同的定階準(zhǔn)則[9]。例如，段志善等根據(jù)辨識(shí)精度的需求，從頻域角度出發(fā)，分低階和高階兩種辨識(shí)情形提出了FFT－AR 譜對(duì)照的定階方法，并在實(shí)際的機(jī)械故障診斷中得到了較為滿意的結(jié)果[10]；楊帆等針對(duì)模型定階過程中的不確定性，利用冗余定理篩選基于多個(gè)回溯階預(yù)測(cè)方法的有效信息，設(shè)計(jì)實(shí)現(xiàn)了一種基于遺忘因子的變權(quán)組合定階方法[11]；黃雄波等通過對(duì)平穩(wěn)時(shí)序數(shù)據(jù)自相關(guān)函數(shù)矩陣的秩的下界值進(jìn)行估計(jì)，并以該估計(jì)值作為自回歸模型的起始階數(shù)對(duì)系統(tǒng)進(jìn)行依次的遞階辨識(shí)，得到一種自回歸模型的快速辨識(shí)算法[12]。

然而，以上方法并沒有把定階過程和估參過程有機(jī)地融合在一起，即整個(gè)模型的求解過程中仍然遵循著先定階后估參的次序，這嚴(yán)重制約了算法的計(jì)算效能的提升。為解決這一問題，張思宇等通過預(yù)留少量數(shù)據(jù)遞補(bǔ)進(jìn)入求解模型的方式，并結(jié)合矩陣遞推求逆的方法，提出了一種融合定階和估參同步進(jìn)行的自回歸模型求解算法[13]。該算法的主要思想與演算步驟如下：

基于最小二乘理論定義如下式所示的辨識(shí)殘差平方和：

將樣本序列排列為{x-v,x-v+1,…,x-1,x0,x1,…,xt}，其中子序列{x1,…,xt}稱為候補(bǔ)數(shù)據(jù)序列。對(duì)于AR(p)模型而言，選用的數(shù)據(jù)序列為{x-p+1,…,xp}，求式(2)的最小二乘解，有：

其中：

每次升階時(shí)，即AR(p)→AR(p+1)，則可由候補(bǔ)數(shù)據(jù)序列依次遞補(bǔ)一個(gè)數(shù)據(jù)，此時(shí)式(3)中的X、Y 矩陣應(yīng)調(diào)整為：

其中，D=H-GE-1F。

上述算法在自回歸模型的求解過程中，以升階的方式進(jìn)行定階，并利用矩陣分塊求逆的方式同步進(jìn)行估參，其優(yōu)點(diǎn)是獲得了較好的計(jì)算性能，然而此算法本身也存在兩點(diǎn)需要改進(jìn)的地方：1) 從式(1)可見，式(3)的自回歸模型僅考慮了數(shù)據(jù)樣本既有的觀測(cè)值，而忽視了白噪聲隨機(jī)干擾ut的作用，因此將引起一定的辨識(shí)誤差，事實(shí)上，根據(jù)Yule-Walker 方程可知，僅當(dāng)系數(shù)矩陣X 的元素為自相關(guān)函數(shù)時(shí)，式(3)才得以成立；2）以分塊矩陣的遞推求逆方式進(jìn)行自回歸模型求解，其矩陣求逆過程中所積累的誤差也不容忽略，應(yīng)予以消除。

針對(duì)上述問題，擬設(shè)計(jì)一種更為高效的平穩(wěn)時(shí)序數(shù)據(jù)的自回歸模型辨識(shí)改進(jìn)算法。

3 改進(jìn)算法

3.1 基于重抽樣的隨機(jī)白噪聲干擾估計(jì)

如前所述，式(3)在利用候補(bǔ)數(shù)據(jù)序列進(jìn)行升階辨識(shí)的過程中，由于忽略了隨機(jī)白噪聲ut的作用而導(dǎo)致了辨識(shí)結(jié)果的不可靠，對(duì)此，引入重抽樣方法，從辨識(shí)序列xt中估算并消除ut，進(jìn)而使得式(3)的辨識(shí)計(jì)算更為精確。

注意到隨機(jī)白噪聲ut為一具有高斯正態(tài)分布的完全隨機(jī)型序列，據(jù)此，需要對(duì)其數(shù)學(xué)期望和方差進(jìn)行估計(jì)。對(duì)于一個(gè)容量為n 且來自未知分布F的數(shù)據(jù)樣本而言，近年來出現(xiàn)了一種名為重抽樣（Bootstrap）的實(shí)用統(tǒng)計(jì)推斷方法[14-15]，其推斷過程如圖1 所示。

圖1 Bootstrap重抽樣的統(tǒng)計(jì)推斷過程

由圖中可見，此算法是先從既有數(shù)據(jù)樣本中按放回抽樣的方法抽取一個(gè)容量為n 的樣本，并稱這種抽樣樣本為Bootstrap 樣本；之后相繼地、獨(dú)立地自原始樣本中抽取一定數(shù)量的Bootstrap 樣本，并利用這些Bootstrap 樣本對(duì)總體F 進(jìn)行統(tǒng)計(jì)推斷。

根據(jù)以上的Bootstrap 重抽樣的統(tǒng)計(jì)推斷分析，估計(jì)隨機(jī)白噪聲ut的數(shù)學(xué)期望和方差的流程應(yīng)為：

1）用隨機(jī)數(shù)發(fā)生器對(duì)既有樣本進(jìn)行N 次重抽樣，得到N 組新樣本；

3) 分別計(jì)算總體樣本的Bootstrap 數(shù)學(xué)期望估值和方差估值。

得到了隨機(jī)白噪聲ut的Bootstrap 數(shù)學(xué)期望估值和Bootstrap 方差估值后，便可依下式所示的概率密度函數(shù)生成完全隨機(jī)型序列：

在此基礎(chǔ)上，運(yùn)用下式即可從辨識(shí)序列xt中消除隨機(jī)白噪聲ut的作用，得到相依隨機(jī)型序列：

3.2 融合迭代和遞推機(jī)制的

自回歸模型求解方法

文獻(xiàn)[13]用矩陣分塊求逆的遞推方法對(duì)升階后的自回歸模型進(jìn)行求解，但矩陣求逆過程中所積累的誤差不能忽視，故此引入融合迭代和遞推機(jī)制的求解方法，通過迭代法對(duì)模型進(jìn)行求解，在迭代計(jì)算收斂的條件下便可以將計(jì)算誤差控制在某一特定范圍內(nèi)。

則有[16-17]：

而剩余向量為：

滿足：

從候補(bǔ)數(shù)據(jù)序列遞補(bǔ)一個(gè)新的數(shù)據(jù)xp+1，得到第p+1 個(gè)方程：

于是，有：

對(duì)式(14)的系數(shù)矩陣進(jìn)行如下的分塊：

其中，X11∈Rp×p，X12∈Rp×1，X21∈R1×p，X22∈R 。令,Y1, f1∈Rp，,Y2, f2∈R，則從式(11)～(12)有：

將式(15)的方程組擴(kuò)展為增廣系統(tǒng)，有：

設(shè)：

則式（16）的超松弛迭代（SOR）求解表達(dá)式為[18]：

式(17)中，ω∈(0,2)，稱為松弛因子；而迭代矩陣則為：

給定某一迭代求解誤差值，便可運(yùn)用式(17)迭代求解升階后的模型參數(shù)，然后從候補(bǔ)數(shù)據(jù)序列遞補(bǔ)一個(gè)新的數(shù)據(jù)繼續(xù)進(jìn)行升階辨識(shí)。如果出現(xiàn)rank(Xp+1) ＜p+1，則確定當(dāng)前的自回歸模型的階數(shù)為p，并同時(shí)結(jié)束整個(gè)辨識(shí)過程。

3.3 算法設(shè)計(jì)

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

算法名稱：平穩(wěn)時(shí)序數(shù)據(jù)的高效辨識(shí)改進(jìn)算法

輸入：平穩(wěn)時(shí)序數(shù)據(jù)序列xt(t=1,2,…,n)、Bootstrap 重抽樣次數(shù)N、迭代求解誤差e；

步驟1：對(duì)xt進(jìn)行N 次Bootstrap 重抽樣，并估算總體樣本的Bootstrap 數(shù)學(xué)期望估值和Bootstrap方差估值；

步驟2：根據(jù)步驟1 的計(jì)算結(jié)果，按照式(7)生成完全隨機(jī)型序列，并利用式(8)析出相依隨機(jī)型序列，同時(shí)置p=0；

步驟3：從候補(bǔ)數(shù)據(jù)序列遞補(bǔ)一個(gè)新的數(shù)據(jù)xp+1，在迭代求解誤差e 的約束條件下，依照式(17)的迭代和遞推機(jī)制計(jì)算升階后的參數(shù)向量

步驟4：若rank(Xp+1)＜p+1 成立，則跳轉(zhuǎn)至步驟6，否則跳轉(zhuǎn)至步驟5；

步驟5：p=p+1，跳轉(zhuǎn)至步驟3；

步驟6：打印輸出計(jì)算結(jié)果并結(jié)束算法。

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

為驗(yàn)證上述改進(jìn)算法的有效性及先進(jìn)性，選取不同階數(shù)，并疊加不同高斯正態(tài)分布的完全隨機(jī)型序列的平穩(wěn)時(shí)序數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)在PC 機(jī)上進(jìn)行，其硬件配置為：

Intel 酷睿i5 4570 四核CPU，

Kingmax DDR3 16GB RAM，

Samsung 850PRO SSD 固態(tài)硬盤；

操作系統(tǒng)與開發(fā)環(huán)境為：

Microsoft Windows 10，

Microsoft Visual Studio 2010 集成開發(fā)環(huán)境中的C++。

在實(shí)驗(yàn)過程中，著重從辨識(shí)精度和計(jì)算成本等技術(shù)指標(biāo)方面與現(xiàn)有算法進(jìn)行對(duì)比，并就相關(guān)結(jié)果加以詳細(xì)分析和討論。

4.1 實(shí)驗(yàn)過程與方法

為不失一般性，使用階數(shù)分別為2、4 及6 共三個(gè)平穩(wěn)時(shí)序數(shù)據(jù)模型進(jìn)行實(shí)驗(yàn)。具體的模型參數(shù)如表1 所示。此外，在實(shí)驗(yàn)的過程中，分別在上述模型疊加了數(shù)學(xué)期望和方差為{(0,1),(1,2)}的完全隨機(jī)型序列。

表1 實(shí)驗(yàn)選用數(shù)據(jù)模型具體參數(shù)

對(duì)上述不同階數(shù)的自回歸模型、不同數(shù)學(xué)期望和方差的完全隨機(jī)型序列進(jìn)行組合，便可得到不同的實(shí)驗(yàn)組合，實(shí)驗(yàn)過程中將對(duì)文獻(xiàn)[13]算法和本改進(jìn)算法的辨識(shí)性能指標(biāo)進(jìn)行相關(guān)的比對(duì)。

4.2 實(shí)驗(yàn)的結(jié)果與分析

為了全面考察各算法的辨識(shí)精度，引入平均絕對(duì)百分誤差（Mean Absolute Percentage Error, MAPE）的性能評(píng)價(jià)標(biāo)準(zhǔn)，其具體定義如式(19)所示[19]：

各種實(shí)驗(yàn)組合的辨識(shí)結(jié)果如表2 所示，它們的辨識(shí)誤差曲線則如圖2～圖4 所示。

表2 各種實(shí)驗(yàn)組合的辨識(shí)結(jié)果

圖2 AR（2）模型實(shí)驗(yàn)組合的辨識(shí)誤差曲線

圖3 AR（4）模型實(shí)驗(yàn)組合的辨識(shí)誤差曲線

圖4 AR（6）模型實(shí)驗(yàn)組合的辨識(shí)誤差曲線

通過對(duì)比表2 中各種實(shí)驗(yàn)組合的辨識(shí)結(jié)果，結(jié)合圖2～圖4 的辨識(shí)誤差曲線，便可發(fā)現(xiàn)：模型階數(shù)和完全隨機(jī)型序列均在不同程度上影響著算法的辨識(shí)精度。其中，文獻(xiàn)[13]算法的辨識(shí)精度與模型階數(shù)有著高度的非線性增長(zhǎng)關(guān)系，在階數(shù)從2 升至4、6 的過程中，其辨識(shí)精度也從高精度辨識(shí)(MAPE≤10 )向良好辨識(shí)（10＜MAPE≤20%）、可行辨識(shí)（20＜MAPE≤30%）以至粗略辨識(shí)（30＜MAPE≤50%）切換；經(jīng)過對(duì)比明顯可見，改進(jìn)算法的辨識(shí)精度具有良好的穩(wěn)定性，在升階過程中，其辨識(shí)精度始終維持在高精度辨識(shí)的范圍內(nèi)。事實(shí)上，本改進(jìn)算法在通過Bootstrap 重抽樣對(duì)完全隨機(jī)型序列進(jìn)行估計(jì)和消除后，其自回歸模型方程有了更為精確的數(shù)學(xué)描述，從而保證了估參過程中的誤差沒有隨著階數(shù)的增加而增加。另一方面，具有較大的數(shù)學(xué)期望和方差的完全隨機(jī)型序列對(duì)算法的辨識(shí)精度也有一定的制約。對(duì)于文獻(xiàn)[13]算法而言，由于忽略了的作用，其自回歸模型方程的數(shù)學(xué)描述出現(xiàn)了更大的偏差；對(duì)于本改進(jìn)算法而言，由于完全隨機(jī)型序列在原序列中的成分比重加大，在某一程度上也影響了Bootstrap 重抽樣估計(jì)的精度，但根據(jù)表2 的MAPE數(shù)值可知，該影響對(duì)改進(jìn)算法而言還是較為輕微的。

為了評(píng)價(jià)融合迭代和遞推機(jī)制的模型求解方法的計(jì)算性能，在實(shí)驗(yàn)過程中，利用Bootstrap 重抽樣從待辨識(shí)序列中析出完全隨機(jī)型序列后，分別用文獻(xiàn)[13]的分塊矩陣遞推求逆方法、改進(jìn)后的融合迭代和遞推機(jī)制的求解法對(duì)同一辨識(shí)序列進(jìn)行模型求解，結(jié)果如表3 所示。融合迭代和遞推機(jī)制的求解法較分塊矩陣遞推求逆方法更為精確。此外，對(duì)比表2 和表3 的MAPE 數(shù)值，也可發(fā)現(xiàn)：文獻(xiàn)[13]算法的辨識(shí)誤差主要是由忽略了完全隨機(jī)型序列而引起的。

表3 不同模型估參法的辨識(shí)精度對(duì)比

在實(shí)驗(yàn)序列的樣本長(zhǎng)度n 為300，Bootstrap 重抽樣次數(shù)N 為1000 的條件下，各種實(shí)驗(yàn)組合的計(jì)算耗時(shí)情況如圖5 所示。由圖易知，不同類型的完全隨機(jī)型序列并不影響文獻(xiàn)[13]算法及本改進(jìn)算法的計(jì)算耗時(shí)。

從圖5 可見，在不同模型階數(shù)的辨識(shí)過程中，本改進(jìn)算法的計(jì)算耗時(shí)較文獻(xiàn)[13]算法增加了12%～15%左右。值得注意的是，改進(jìn)算法由于設(shè)計(jì)了融合迭代和遞推機(jī)制的模型參數(shù)估計(jì)方法，其計(jì)算耗時(shí)并沒有隨著模型階數(shù)的增加而大幅增長(zhǎng)，事實(shí)上，本算法所增加的計(jì)算耗時(shí)是由于引入了Bootstrap 重抽樣處理而引起的，該耗時(shí)增量為定值，且僅與抽樣次數(shù)N 有關(guān)。

圖5 各種實(shí)驗(yàn)?zāi)Ｐ徒M合的計(jì)算耗時(shí)

綜上所述，改進(jìn)算法僅花費(fèi)了少量的額外計(jì)算耗時(shí)，其辨識(shí)精度穩(wěn)定性便有了顯著提升。據(jù)此，本改進(jìn)算法是正確和有效的。

5 結(jié)束語

基于時(shí)序數(shù)據(jù)驅(qū)動(dòng)的系統(tǒng)辨識(shí)，在實(shí)際應(yīng)用中有著廣泛的適用場(chǎng)合，在綜合考慮已有算法的優(yōu)劣的基礎(chǔ)上，提出一種平穩(wěn)時(shí)序數(shù)據(jù)的自回歸模型高效辨識(shí)改進(jìn)算法。經(jīng)過實(shí)驗(yàn)得知，改進(jìn)算法的辨識(shí)精度具有較強(qiáng)的魯棒性。此改進(jìn)算法還有繼續(xù)完善的空間，下一步的主要工作應(yīng)包括：改進(jìn)現(xiàn)有的Bootstrap 重抽樣的估計(jì)效率，研究更為高效的迭代和遞推的求解機(jī)制，以便進(jìn)一步提升算法的預(yù)測(cè)精度和計(jì)算效能。