黃雄波

(佛山職業(yè)技術(shù)學(xué)院 電子信息系, 佛山 528000)

0 引言

所謂時(shí)序數(shù)據(jù)，就是將某一物理量在不同時(shí)間上的實(shí)際觀測(cè)值，按照時(shí)間的先后順序排列而成的數(shù)列。 “平穩(wěn)”是時(shí)序數(shù)據(jù)的重要特性，它描述了時(shí)序數(shù)據(jù)的數(shù)字統(tǒng)計(jì)特性是否隨著時(shí)間的推移而變化，通常，若時(shí)序數(shù)據(jù)的平均功率存在，其均值為常數(shù)，且自相關(guān)函數(shù)與起始時(shí)間無(wú)關(guān)，則稱該時(shí)序數(shù)據(jù)為廣義平穩(wěn)時(shí)序數(shù)據(jù)或弱平穩(wěn)時(shí)序數(shù)據(jù)，簡(jiǎn)稱為平穩(wěn)時(shí)序數(shù)據(jù)[1-2]。經(jīng)過(guò)眾多專家學(xué)者的長(zhǎng)期研究，平穩(wěn)時(shí)序數(shù)據(jù)現(xiàn)已形成了一系列成熟可靠的辨識(shí)模型及建模方法，以自回歸模型(Auto-regressive model，AR)為例，該模型因結(jié)構(gòu)簡(jiǎn)單且具有高效的參數(shù)迭代估計(jì)算法，故在實(shí)際應(yīng)用中得到了廣泛的應(yīng)用[3-4]。然而，在現(xiàn)實(shí)生活中所獲取的時(shí)序數(shù)據(jù)普遍具有非平穩(wěn)特性，又由于非平穩(wěn)時(shí)序數(shù)據(jù)尚未有完整和統(tǒng)一的描述方法，據(jù)此，人們往往需要對(duì)其進(jìn)行相關(guān)的平穩(wěn)預(yù)處理[5-7]后，再采用平穩(wěn)時(shí)序數(shù)據(jù)建模方法進(jìn)行后續(xù)的辨識(shí)。例如，博克思(Box)和詹金斯(Jenkins)于上世紀(jì)70年代提出的自回歸積分滑動(dòng)平均模型(Auto-regressive Integrated Moving Average Model，ARIMA)，其核心就是通過(guò)有限次的差分處理，把原來(lái)的非平穩(wěn)時(shí)序數(shù)據(jù)轉(zhuǎn)化為平穩(wěn)時(shí)序數(shù)據(jù)[8]；黃雄波基于自相關(guān)函數(shù)理論，對(duì)非平穩(wěn)時(shí)序數(shù)據(jù)中的趨勢(shì)和周期成分的分離次序和方法作了系統(tǒng)而深入的研究，并得到了一種有效的平穩(wěn)化轉(zhuǎn)換算法[9]；王文華等基于貝葉斯框架對(duì)非平穩(wěn)時(shí)序數(shù)據(jù)的分段平穩(wěn)問(wèn)題作了更為深刻的研究，并推導(dǎo)出具有遞歸關(guān)系的高效建模算法[10]；張海勇等基于經(jīng)驗(yàn)?zāi)B(tài)分解法(Empirical mode decomposition, EMD)把待處理的非平穩(wěn)時(shí)序數(shù)據(jù)分解成有限個(gè)基本的模式分量，然后分別為這些模式分量建立時(shí)變的AR模型，進(jìn)而得到一種新的非平穩(wěn)時(shí)序數(shù)據(jù)的自回歸模型分析方法[11-12]。

為了改善自回歸模型的辨識(shí)精度并給出相應(yīng)的置信空間，近年來(lái)，國(guó)內(nèi)外一些專家學(xué)者開(kāi)始將Bootstrap方法應(yīng)用到平穩(wěn)時(shí)序數(shù)據(jù)的自回歸辨識(shí)過(guò)程中。例如，軒建平等基于小樣本統(tǒng)計(jì)的 Bootstrap 方法，對(duì)車削振顫自回歸模型參數(shù)的方差進(jìn)行了估計(jì)，從而較好地解決了機(jī)床故障診斷中所需的大量樣本和重復(fù)多次試驗(yàn)的難題[13]；楊曉蓉等在單位根存在的條件下，證明了自回歸模型所構(gòu)造的單位根檢驗(yàn)統(tǒng)計(jì)量的極限分布，可以通過(guò)對(duì)最小二乘殘差進(jìn)行Bootstrap重抽樣方法來(lái)逼近[14]。目前，基于Bootstrap方法的自回歸辨識(shí)算法在實(shí)際應(yīng)用中還存在著一些問(wèn)題，據(jù)此，本文擬設(shè)計(jì)相應(yīng)的改進(jìn)算法予以解決，并以相關(guān)的實(shí)驗(yàn)來(lái)證明改進(jìn)算法的有效性和先進(jìn)性。

1 問(wèn)題描述

時(shí)序數(shù)據(jù)辨識(shí)的過(guò)程就是對(duì)序列的數(shù)字統(tǒng)計(jì)特征及其分布的不確定性進(jìn)行評(píng)價(jià)，在實(shí)際的辨識(shí)過(guò)程中，我們可以采用先驗(yàn)知識(shí)、概率模型、可能性及置信區(qū)間等來(lái)求解不確定關(guān)系的表達(dá)式或者對(duì)其進(jìn)行估計(jì)。然而，在大多數(shù)情況下，由于我們并不掌握待分析序列的總體分布，故很可能會(huì)得到錯(cuò)誤的辨識(shí)結(jié)果。為克服因樣本數(shù)量不足而引起的分布估計(jì)誤差，一種有效的方法就是重抽樣，即通過(guò)某種模型對(duì)原序列進(jìn)行反復(fù)采樣，使得原序列容量得到擴(kuò)充，從而可以較為準(zhǔn)確地估計(jì)某一統(tǒng)計(jì)量的標(biāo)準(zhǔn)差、均值及概率分布等。

(1)

式(1)中，n為樣本長(zhǎng)度，N(y)指時(shí)間y發(fā)生的總次數(shù)。

如圖1所示。

圖1 基于Bootstrap方法的自回歸辨識(shí)的算法流程圖

基于Bootstrap方法的自回歸辨識(shí)其核心思想就是對(duì)殘差序列進(jìn)行m次重抽樣，然后構(gòu)造出一系列的Bootstrap序列，通過(guò)對(duì)這些殘差序列進(jìn)行參數(shù)估計(jì)，并用所得的參數(shù)均值來(lái)修正原有的自回歸模型。

2 改進(jìn)算法的設(shè)計(jì)與實(shí)現(xiàn)

針對(duì)上述問(wèn)題，本文擬對(duì)原有算法作如下改進(jìn)：(1) 以自回歸模型的階數(shù)為Bootstrap的作用范圍區(qū)間，對(duì)殘差序列進(jìn)行重抽樣處理；(2) 基于矩陣奇異值的迭代分解理論，對(duì)Bootstrap序列的參數(shù)進(jìn)行求解。

2.1 改進(jìn)的Bootstrap序列生成算法

在理想化的情形下，利用自回歸模型對(duì)平穩(wěn)時(shí)序數(shù)據(jù)進(jìn)行辨識(shí)，其剩余的殘差序列Et應(yīng)為一白噪聲，即Et～N(0,σ2)成立。然而，由于數(shù)據(jù)污染和計(jì)算誤差等原因，殘差序列Et往往不能被白化，其各時(shí)刻之間的取值仍然存在著某種的統(tǒng)計(jì)關(guān)聯(lián)性。據(jù)此，在Bootstrap重抽樣的過(guò)程中，其作用范圍區(qū)間不應(yīng)在整個(gè)序列內(nèi)進(jìn)行，而應(yīng)以自回歸模型階數(shù)p作為滑動(dòng)窗口的寬度，對(duì)殘差序列Et進(jìn)行重抽樣，以便保持既有的相關(guān)性。根據(jù)上述的分析，可設(shè)計(jì)一種改進(jìn)的Bootstrap序列生成算法。

算法1：保持既有相關(guān)性的Bootstrap序列生成算法

輸入：殘差序列Et，自回歸模型階數(shù)p；

步驟1：用Random()函數(shù)在1～n(n為序列的長(zhǎng)度)中產(chǎn)生一個(gè)隨機(jī)整數(shù)R；

步驟3：應(yīng)用Bootstrap方法獨(dú)立地對(duì)步驟2中的各個(gè)子空間中進(jìn)行重抽樣處理；

2.2 改進(jìn)的Bootstrap序列估參算法

Bφ=b.

(2)

式(2)中，

(3)

奇異值分解方法是最小二乘問(wèn)題的有效求解方法，該方法的主要優(yōu)點(diǎn)是處理病態(tài)和不相容線性方程組的能力強(qiáng)、運(yùn)算過(guò)程中不放大誤差且具有良好的穩(wěn)定性；而缺點(diǎn)就是運(yùn)算量大。據(jù)此，這里擬引入矩陣奇異值分解方法對(duì)Bootstrap序列的自回歸模型參數(shù)進(jìn)行求解，同時(shí)，為了克服運(yùn)算量大的問(wèn)題，還需要對(duì)奇異值分解方法進(jìn)行如下的迭代計(jì)算改進(jìn)。

設(shè)B(B∈R(n-1-p)×P,n>2p+1)的奇異值分解為式(4)。

(4)

其中，U=[u1,…,un-1-p]和V=[v1,…,vp]是正交矩陣,∑r=diag(σ1,…,σr),σ1≥…≥σr>0。根據(jù)矩陣Moore-Penrose廣義逆的定義，式(2)的最小二乘解為式(5)。

(5)

為了能使用迭代法對(duì)系數(shù)矩陣B進(jìn)行奇異值分解，可以運(yùn)用Householder變換求得正交矩陣U,V,從而實(shí)現(xiàn)系數(shù)矩陣B的二對(duì)角化，即式(6)。

(6)

其中，

1) 初始化：

①k=1;

④c=α1-μ;

⑤d=β1;

2) 求解如下的矩陣方程，得到奇異值σk：

3)k=k+1，求解如下的矩陣方程，更新θk：

4) 若k

綜上所述，可設(shè)計(jì)如下的Bootstrap序列的估參改進(jìn)算法。

算法2：基于奇異值迭代分解的Bootstrap序列估參算法

步驟2：對(duì)步驟(1)中的線性方程組的系數(shù)矩陣施行Householder變換，求得對(duì)應(yīng)的正交矩陣U,V，把系數(shù)矩陣B化為式(6)所示的二對(duì)角化矩陣B′；

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

為了驗(yàn)證上述改進(jìn)算法的有效性及先進(jìn)性，這里選取了一個(gè)自回歸仿真模型來(lái)進(jìn)行相關(guān)的Bootstrap重抽樣辨識(shí)。實(shí)驗(yàn)在PC機(jī)上進(jìn)行，其硬件配置為，Intel 酷睿i5 4570四核CPU、Kingmax DDR3 16GB RAM、Western Digital 500G Hard Disk；操作系統(tǒng)與開(kāi)發(fā)環(huán)境為，Microsoft Windows 10、Microsoft Visual Studio 2010集成開(kāi)發(fā)環(huán)境中的C++。在實(shí)驗(yàn)過(guò)程中，著重關(guān)注改進(jìn)算法的辨識(shí)誤差和計(jì)算開(kāi)銷等技術(shù)指標(biāo)的改善情況，并對(duì)相關(guān)結(jié)果加以詳細(xì)的分析和討論。

3.1 實(shí)驗(yàn)過(guò)程與方法

實(shí)驗(yàn)選用了一個(gè)五階的自回歸仿真模型，具體數(shù)學(xué)模型如式(7)。

yt=-0.81yt-1+0.21yt-2+0.19yt-3-0.17yt-4+0.33yt-5+zt.

(7)

式(7)中，zt為服從標(biāo)準(zhǔn)高斯分布的白噪聲。

基于MATLAB軟件中，利用randn()隨機(jī)數(shù)發(fā)生函數(shù)和filter()數(shù)字濾波器函數(shù)為式(7)所示的模型生成長(zhǎng)度為600的樣本序列。首先基于Burg算法對(duì)式(7)的模型進(jìn)行辨識(shí)，然后分別應(yīng)用現(xiàn)有的自回歸模型Bootstrap辨識(shí)算法、本文算法1、本文算法2及本文算法1+算法2對(duì)殘差序列進(jìn)行辨識(shí)，并在{200,250,300,…,600}等不同抽樣次數(shù)的情形下，比對(duì)各種算法的辨識(shí)精度和辨識(shí)耗時(shí)的性能表現(xiàn)。

實(shí)驗(yàn)過(guò)程中，為了客觀準(zhǔn)確地評(píng)價(jià)各種算法的辨識(shí)精度，這里以式 (8)所示的平均絕對(duì)百分誤差(Mean absolute percent error, MAPE)作為評(píng)價(jià)指標(biāo)，對(duì)各種實(shí)驗(yàn)組合的辨識(shí)精度進(jìn)行評(píng)價(jià)。如式(8)。

(8)

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

用Burg算法對(duì)式(7)所示的模型進(jìn)行辨識(shí)，其結(jié)果如式(9)所示，于是，便可得到式(10)所示的殘差序列Et如式(9)、(10)。

(9)

(10)

分別用現(xiàn)有算法、本文算法1、本文算法2及本文算法1+算法2對(duì)殘差序列Et進(jìn)行自回歸辨識(shí)，取每次重抽樣其辨識(shí)參數(shù)的均值與式(9)中對(duì)應(yīng)參數(shù)進(jìn)行疊加求和，從而實(shí)現(xiàn)對(duì)原Burg算法辨識(shí)模型的修正。實(shí)驗(yàn)所得的最終辨識(shí)參數(shù)如表1所示。

表1 各種算法所得的最終辨識(shí)參數(shù)

而各種算法的MAPE與重抽樣次數(shù)的關(guān)系則如圖2所示，限于篇幅，這里僅列出抽樣次數(shù)為200,400,600共3種情況。

圖2 各種算法的MAPE與重抽樣次數(shù)的關(guān)系

從表1可以發(fā)現(xiàn)，由于本文算法1在重抽樣的過(guò)程中保持了殘差序列之間既有的相關(guān)性，故其辨識(shí)精度較現(xiàn)有算法有了3%左右的提升；同樣地，基于奇異值迭代分解法的本文算法2因其數(shù)值計(jì)算更為精確，相應(yīng)地，辨識(shí)精度較現(xiàn)有算法也有了一定的提升。而本文算法1+算法2則是融合了兩種改進(jìn)算法的優(yōu)點(diǎn)，在保持殘差序列之間既有相關(guān)性的同時(shí)，其估參過(guò)程中所產(chǎn)生的計(jì)算誤差也得到了較好的控制，故其對(duì)應(yīng)的辨識(shí)精度也是最高的。如圖2所示。

MAPE與重抽樣次數(shù)之間的關(guān)系則表明，各種算法的重抽樣次數(shù)越多，其辨識(shí)精度也就越高，但當(dāng)重抽樣次數(shù)達(dá)到序列長(zhǎng)度的2/3 (400次)之后，辨識(shí)精度的提升效果便不再顯著，究其原因是因?yàn)榇藭r(shí)的重抽樣操作已基本覆蓋了殘差序列的未知分布。

各種算法的計(jì)算耗時(shí)，如圖3所示。

圖3 各種算法的計(jì)算耗時(shí)

從圖3可以得知，本文算法1花費(fèi)的計(jì)算耗時(shí)最多，現(xiàn)有算法次之，本文算法1+算法2排第三，而本文算法2則為最小。

以實(shí)驗(yàn)中的模型為例，現(xiàn)有算法的每一次Bootstrap過(guò)程，均需要花費(fèi)一定的時(shí)間進(jìn)行重抽樣并生成Bootstrap序列，而且還需基于最小二乘法求解一個(gè)行數(shù)和列數(shù)分別為600和5的線性方程組，以便估算出對(duì)應(yīng)的自回歸模型參數(shù)φ。實(shí)驗(yàn)表明，現(xiàn)有算法的計(jì)算耗時(shí)主要是消耗在求解高度不相容線性方程組的過(guò)程中，據(jù)此，其計(jì)算耗時(shí)與重抽樣次數(shù)也就存在著明顯的線性關(guān)系。由于本文算法1的Bootstrap重抽樣過(guò)程較現(xiàn)有算法復(fù)雜，且自回歸模型參數(shù)φ的求解方法又與現(xiàn)有算法相同，故其計(jì)算耗時(shí)在略高于現(xiàn)有算法的同時(shí)又與其有著相同的特征。相對(duì)地，本文算法2基于奇異值迭代分解法改進(jìn)了自回歸模型參數(shù)φ的求解過(guò)程，故其計(jì)算耗時(shí)有了大幅度的減少，并表現(xiàn)出對(duì)重抽樣次數(shù)具有良好的負(fù)載能力，即計(jì)算耗時(shí)不隨重抽樣次數(shù)的增加而顯著增加。而本文算法1+算法2的計(jì)算耗時(shí)在略高于本文算法2的同時(shí)，卻能保持著本文算法2的良好特性。

綜上所述，本文算法1+算法2在結(jié)合了兩種改進(jìn)算法的基礎(chǔ)上，其計(jì)算精度和計(jì)算耗時(shí)均較現(xiàn)有算法有了顯著的提升。

4 總結(jié)

對(duì)現(xiàn)有的基于Bootstrap方法的自回歸模型辨識(shí)算法，

進(jìn)行了兩點(diǎn)有意義的改進(jìn)，改進(jìn)后的算法具有更優(yōu)異的辨識(shí)性能指標(biāo)。下一步的主要工作有，研究更為合理的Bootstrap重抽樣約束機(jī)制，同時(shí)，也需要研究并行的奇異值迭代分解法，以便進(jìn)一步提升算法的適用范圍和計(jì)算性能。

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

[2] 冀振元.時(shí)間序列分析與現(xiàn)代譜估計(jì)[M].哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2016.

[3] 黃雄波，胡永健.利用自回歸模型的平穩(wěn)時(shí)序數(shù)據(jù)快速辨識(shí)算法[J/OL，網(wǎng)絡(luò)優(yōu)先出版].計(jì)算機(jī)應(yīng)用研究，2018,35(9).

[4] 蘇志銘，陳靚影.基于自回歸模型的動(dòng)態(tài)表情識(shí)別[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2017, 29(6):1085-1092.

[5] 王宏禹，邱天爽，陳喆.非平穩(wěn)隨機(jī)信號(hào)分析與處理(第2版)[M]. 北京：國(guó)防工業(yè)出版社，2008.

[6] 王宏禹，邱天爽.非平穩(wěn)確定性信號(hào)與非平穩(wěn)隨機(jī)信號(hào)統(tǒng)一分類法的探討[J]. 通信學(xué)報(bào)，2015, 36(2):2801-2810.

[7] 王宏禹，邱天爽.確定性信號(hào)分解與平穩(wěn)隨機(jī)信號(hào)分解的統(tǒng)一研究[J]. 通信學(xué)報(bào)，2016, 37(10):1891-1898.

[8] 博克思，詹金斯，萊因澤爾.時(shí)間序列分析：預(yù)測(cè)與控制(第4版)[M].王成璋譯.北京：機(jī)械工業(yè)出版社，2011.

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

[10] 王文華，王宏禹.分段平穩(wěn)隨機(jī)過(guò)程的參數(shù)估計(jì)方法[J]. 電子科學(xué)學(xué)刊，1997, 19(3):311-317.

[11] 張海勇，馬孝江，蓋強(qiáng).一種新的時(shí)變參數(shù)AR模型分析方法[J]. 大連理工大學(xué)學(xué)報(bào)，2002, 42(2):238-241.

[12] 張海勇，李勘.非平穩(wěn)隨機(jī)信號(hào)的參數(shù)模型分析方法[J]. 系統(tǒng)工程與電子技術(shù)，2003, 25(3):386-390.

[13] 軒建平，史鐵林，楊叔子. AR模型參數(shù)的Bootstrap方差估計(jì)[J]. 華中科技大學(xué)學(xué)報(bào)，2001, 29(9):81-83.

[14] 楊曉蓉.自回歸時(shí)間序列的極限理論及其應(yīng)用[D].浙江大學(xué)博士學(xué)位論文，2008.

[15] 徐禮文.復(fù)雜數(shù)據(jù)的bootstrap統(tǒng)計(jì)推斷及其應(yīng)用[M].北京：科學(xué)出版社，2016.

[16] James W.Demmel.應(yīng)用數(shù)值線性代數(shù)[M] 王國(guó)榮譯.北京：人民郵電出版社，2007.

[17] 徐樹(shù)方.數(shù)值線性代數(shù)(第二版)[M].北京：北京大學(xué)出版社，2013.