基于ARIMA-KF模型的船舶系統(tǒng)設(shè)備狀態(tài)參數(shù)預(yù)測(cè)

陳方圓， 鄒永久， 張 鵬， 張躍文， 孫培廷

(大連海事大學(xué)輪機(jī)工程學(xué)院， 大連 116026)

隨著大數(shù)據(jù)技術(shù)的快速發(fā)展和5G時(shí)代的到來，智能船舶已成為船舶行業(yè)未來發(fā)展的必然趨勢(shì)[1]。船舶系統(tǒng)設(shè)備狀態(tài)參數(shù)預(yù)測(cè)研究也將成為船舶智能化發(fā)展的重要基礎(chǔ)技術(shù)之一?；跀?shù)據(jù)驅(qū)動(dòng)的方法，如何有效地利用航參數(shù)據(jù)，建立預(yù)測(cè)模型，實(shí)現(xiàn)對(duì)船舶系統(tǒng)及設(shè)備工作狀態(tài)的智能分析和準(zhǔn)確預(yù)測(cè)，提高智能船舶的可靠性和船舶智能設(shè)備的可操作性，是目前亟需解決的問題[2]。此外，船舶系統(tǒng)設(shè)備狀態(tài)參數(shù)預(yù)測(cè)不僅可以實(shí)現(xiàn)船舶系統(tǒng)設(shè)備剩余壽命預(yù)測(cè)，還能為輪機(jī)員提供輔助決策支持[3]。

目前，在船舶領(lǐng)域中常用的預(yù)測(cè)方法主要有灰色模型、BP(back propagation)神經(jīng)網(wǎng)絡(luò)等。南曉雪[4]采用灰色模型實(shí)現(xiàn)了對(duì)有變化趨勢(shì)的港口船舶通航能力進(jìn)行了準(zhǔn)確預(yù)測(cè)。但灰色模型只能識(shí)別單一參數(shù)，限制了其使用范圍，且使用單一算法對(duì)參數(shù)預(yù)測(cè)的精度較低，無法滿足船舶系統(tǒng)設(shè)備狀態(tài)參數(shù)預(yù)測(cè)的實(shí)際需求。王子銘等[5]基于長(zhǎng)短期記憶神經(jīng)網(wǎng)絡(luò)中的時(shí)間序列預(yù)測(cè)模型構(gòu)建了主機(jī)油耗實(shí)時(shí)預(yù)報(bào)模型，實(shí)現(xiàn)了對(duì)船舶在航行過程中主機(jī)油耗參數(shù)的實(shí)時(shí)預(yù)測(cè)。呂鵬飛等[6]以船舶運(yùn)輸量為主要研究對(duì)象，結(jié)合馬爾科夫預(yù)測(cè)模型，對(duì)BP神經(jīng)網(wǎng)絡(luò)進(jìn)行優(yōu)化，實(shí)現(xiàn)了對(duì)運(yùn)輸量趨勢(shì)的準(zhǔn)確預(yù)測(cè)。王新全等[7]構(gòu)建了基于遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型，實(shí)現(xiàn)了對(duì)船舶主機(jī)排氣溫度變化趨勢(shì)的準(zhǔn)確預(yù)測(cè)。尹倩[8]基于長(zhǎng)短期記憶循環(huán)神經(jīng)網(wǎng)絡(luò)的溫度預(yù)測(cè)分析模型，實(shí)現(xiàn)了對(duì)壓縮機(jī)出風(fēng)口和渦輪機(jī)出口的未來溫度變化進(jìn)行了預(yù)測(cè)。上述基于神經(jīng)網(wǎng)絡(luò)進(jìn)行預(yù)測(cè)的方法雖然提高了預(yù)測(cè)的準(zhǔn)確度，但是基于神經(jīng)網(wǎng)絡(luò)方法需要從訓(xùn)練集中選取一些典型的樣本重新構(gòu)建訓(xùn)練集，容易陷入局部收斂，對(duì)樣本依賴性較大，只適用于個(gè)別條件下船舶系統(tǒng)設(shè)備狀態(tài)參數(shù)預(yù)測(cè)。近年來，也有不少專家和學(xué)者通過采用自回歸滑動(dòng)平均(auto-regressive moving-average model，ARMA)模型，實(shí)現(xiàn)了對(duì)船舶系統(tǒng)及設(shè)備狀態(tài)參數(shù)的預(yù)測(cè)。孫曉磊等[9]構(gòu)建了自回歸分布式滯后ARMA預(yù)測(cè)模型，實(shí)現(xiàn)了對(duì)船舶動(dòng)力裝置系統(tǒng)狀態(tài)參數(shù)的預(yù)測(cè)。優(yōu)化后的ARMA模型雖然提高了預(yù)測(cè)的準(zhǔn)確度，但由于船舶工作環(huán)境惡劣，受到的影響因素具有很大的不確定性，增大了優(yōu)化模型的構(gòu)建難度。因此，這在很大程度上限制了上述方法在船舶領(lǐng)域中的實(shí)際應(yīng)用，亟需研究新型的組合預(yù)測(cè)模型。

針對(duì)上述模型的局限性，現(xiàn)探究船舶系統(tǒng)設(shè)備狀態(tài)參數(shù)的變化規(guī)律，提出一種基于ARIMA-KF(auto-regressive integrated moving-average model,Kalman filter，KF)組合模型的船舶系統(tǒng)設(shè)備狀態(tài)參數(shù)的預(yù)測(cè)方法。該模型將兩種算法的優(yōu)勢(shì)結(jié)合在一起，能彌補(bǔ)單一模型的不足，提高模型的準(zhǔn)確度。此外，該模型將外界復(fù)雜的影響因素歸為單一的時(shí)間影響因素來考慮，可降低模型構(gòu)建的難度[10]，并依靠卡爾曼濾波算法強(qiáng)大的更新迭代能力，大大減少了工作量[3]。該模型對(duì)船舶系統(tǒng)狀態(tài)監(jiān)控和輔助決策具有一定的指導(dǎo)意義。

1 ARIMA模型原理以及建模過程

ARIMA(p,d,q)模型又稱整合移動(dòng)平均自回歸模型，其中，p、q為ARIMA模型的階數(shù)，d為模型差分次數(shù)，ARIMA(p,d,q)模型的實(shí)質(zhì)就是最原始的非平穩(wěn)序列做d次差分處理，然后再進(jìn)行ARMA(p,q)模型擬合[10]。ARIMA(p,d,q)模型形式為

(1)

當(dāng)d=0時(shí)，該模型為ARMA(p,q)模型，即

(2)

式中：yt為采集到的時(shí)序數(shù)據(jù)，t代表采樣時(shí)間；φi(i=0,1,2,…，p)、θj(j=1,2,…,q)為ARIMA模型的參數(shù)；{at}為t時(shí)刻白噪聲序列；εt為t時(shí)刻的零均值白噪聲；L為滯后算子。

建模具體過程如下。

(1)平穩(wěn)性檢驗(yàn)。采集到一組實(shí)船數(shù)據(jù)后，首先對(duì)該序列進(jìn)行平穩(wěn)性檢驗(yàn)，一般采用單位根檢驗(yàn)方法。若時(shí)序中存在單位根即為非平穩(wěn)性時(shí)間序列。反之，則為平穩(wěn)性時(shí)間序列，單位根檢驗(yàn)方法一共包含6種檢驗(yàn)方式，一般聯(lián)合使用ADF檢驗(yàn)(augmented dickey-fuller test)、KPSS檢驗(yàn) (kwiatkowski-phillips-schmidt-shin test)可以提高結(jié)果的可靠性[11]。

(2)數(shù)據(jù)預(yù)處理。若原始數(shù)據(jù)沒有通過平穩(wěn)性檢驗(yàn)，則需要對(duì)該序列進(jìn)行零均值處理，一般是通過將d階差分算子作用于原始序列yt，得到ARIMA(p,d,q)模型。一階差分示例為

?2yt=?(yt-yt-1)=yt-2yt-1+yt-2

(3)

式(3)中：?為差分算子；yt、yt-1、yt-2分別表示相鄰時(shí)刻采集的數(shù)據(jù)。

(3)模型判定。對(duì)ARIMA(p,d,q)模型的判定可以結(jié)合自相關(guān)函數(shù)和偏自相關(guān)函數(shù)性質(zhì)進(jìn)行初步估計(jì)，具體判定方法如表1所示。

表1 模型判定方法

根據(jù)估計(jì)得到的不同p、q值，采用AIC準(zhǔn)則(akaike information criterion)和BIC準(zhǔn)則(bayesian information criterion)給模型定階。準(zhǔn)則公式為

(4)

(4)參數(shù)估計(jì)和模型檢驗(yàn)。參數(shù)估計(jì)一般采用yule-walker逆矩陣法。模型檢驗(yàn)一般采用白噪聲檢驗(yàn)方式。

2 卡爾曼濾波

時(shí)間序列模型可以轉(zhuǎn)換成狀態(tài)空間模型的一種形式，并通過卡爾曼濾波進(jìn)行求解[12]。狀態(tài)空間模型由狀態(tài)方程和觀測(cè)方程兩部分組成，運(yùn)用狀態(tài)空間分析方法將時(shí)間作為動(dòng)態(tài)系統(tǒng)的隱藏變量的一種動(dòng)態(tài)時(shí)域模型[13]。

狀態(tài)方程為

Xt=FtXt-1+ct+Rtwt

(5)

量測(cè)方程為

Yt+1=ZtXt+dt+vt

(6)

式中：Xt為t時(shí)刻的狀態(tài)估計(jì)值(t=1,2,3，…，T)；Ft代表m×m維的狀態(tài)轉(zhuǎn)移矩陣；ct為m×1維的常數(shù)向量；Rt是m×g維的噪聲矩陣；Yt+1是系統(tǒng)在t+1時(shí)刻k×1維的觀測(cè)向量；Zt為k×1維量程矩陣；dt是k×1維的干擾項(xiàng)；wt、vt代表的是均值為0、協(xié)方差矩陣分別為Gt、Qt的干擾項(xiàng)，即E(wt)=0，var(wt)=Qt。E(vt)=0，var(vt)=Gt，且Gt、Qt相互獨(dú)立。

卡爾曼濾波(Kalman filtering，KF)由狀態(tài)方程和測(cè)量方程兩部分構(gòu)成，在狀態(tài)方程和測(cè)量方程已知的情況下，從一系列存在觀測(cè)噪聲數(shù)據(jù)中可以迭代出系統(tǒng)的最優(yōu)狀態(tài)估計(jì)[14]。

狀態(tài)方程為

Xt=FtXt-1+Btut+wt

(7)

測(cè)量方程為

Yt=HtXt+vt

(8)

式中：Xt、Ft、Yt分別代表狀態(tài)向量、狀態(tài)轉(zhuǎn)移矩陣和量測(cè)向量；Bt、ut、Ht分別為控制矩陣、狀態(tài)控制向量、轉(zhuǎn)換矩陣；wt、vt分別為系統(tǒng)噪聲和測(cè)量噪聲且都服從高斯分布。在已知t-1時(shí)刻Xt-1的情況下，通過狀態(tài)方程和測(cè)量方程不斷更新和預(yù)測(cè)，就可以得到t時(shí)刻的最優(yōu)狀態(tài)值Xt|t-1。

系統(tǒng)狀態(tài)更新步為

(9)

觀測(cè)量預(yù)測(cè)步為

(10)

式(9)中：

分別為t時(shí)刻的后驗(yàn)狀態(tài)估計(jì)、卡爾曼增益、后驗(yàn)估計(jì)協(xié)方差、先驗(yàn)估計(jì)協(xié)方差；Yt|t-1、St、et分別為單步預(yù)測(cè)結(jié)果、預(yù)測(cè)誤差、預(yù)測(cè)誤差的協(xié)方差。

若已知t-1時(shí)刻系統(tǒng)的最優(yōu)狀態(tài)服從高斯分布的情況下，即Xt-1～(xt-1,Pt-1)，根據(jù)預(yù)測(cè)方程(9)可以直接推導(dǎo)出t時(shí)刻的狀態(tài)預(yù)測(cè)值Xt|t-1，也為t+1時(shí)刻的先驗(yàn)估計(jì)，然后再由測(cè)量方程[式(10)]得到t時(shí)刻的量測(cè)狀態(tài)Yt。將t時(shí)刻的預(yù)測(cè)值和量測(cè)方程服從的高斯分布相乘，得到t時(shí)刻狀態(tài)的最優(yōu)估計(jì)服從的高斯分布，即Xt|t～(xt,Pt)，再返回卡爾曼濾波更新步即式(10)便得到了卡爾曼濾波的結(jié)果最優(yōu)狀態(tài)估計(jì)Xt和最優(yōu)狀態(tài)協(xié)方差估計(jì)Pt。

3 ARIMA-KF組合模型預(yù)測(cè)方法

根據(jù)以往研究表明，在不需要以往存儲(chǔ)的歷史數(shù)據(jù)和已知狀態(tài)向量初始值、狀態(tài)協(xié)方差矩陣初始值的情況下，卡爾曼濾波可通過計(jì)算機(jī)程序?qū)顟B(tài)空間模型參數(shù)進(jìn)行估計(jì)，并且調(diào)整卡爾曼濾波的參數(shù)可以尋求和真實(shí)曲線吻合度更高的仿真曲線，達(dá)到提高預(yù)測(cè)精度的目的[11]。因時(shí)間序列ARIMA可以實(shí)現(xiàn)單步預(yù)測(cè)和多步預(yù)測(cè)，故下面將從單步預(yù)測(cè)和多步預(yù)測(cè)兩方面構(gòu)建線性狀態(tài)參數(shù)的組合預(yù)測(cè)模型，對(duì)船舶系統(tǒng)與設(shè)備狀態(tài)參數(shù)進(jìn)行預(yù)測(cè)研究。

3.1 ARIMA-KF單步預(yù)測(cè)方法

利用ARIMA模型對(duì)船舶系統(tǒng)及設(shè)備狀態(tài)參數(shù)建立預(yù)測(cè)模型，得到初始參數(shù)估計(jì)。具體步驟如下。

Step 1建立時(shí)間序列ARIMA模型。對(duì)狀態(tài)參數(shù)時(shí)序進(jìn)行平穩(wěn)性檢驗(yàn)和季節(jié)性調(diào)整，確定差分階數(shù)，依照AIC準(zhǔn)則、BIC準(zhǔn)則給模型定階，用yule-walker逆矩陣法對(duì)模型系數(shù){φ1,φ2,…,φp}，{θ1，θ2，…，θq}進(jìn)行求解。

Step 2構(gòu)建狀態(tài)空間模型的狀態(tài)轉(zhuǎn)移矩陣。將ARIMA模型系數(shù)值引入狀態(tài)空間模型的m×m維的狀態(tài)轉(zhuǎn)移矩陣Ft中，其中m=max(p,q)，p、q為ARIMA模型階數(shù)。用模型預(yù)測(cè)方程來表示狀態(tài)空間模型的狀態(tài)轉(zhuǎn)移方程。綜上所述，ARIMA模型系數(shù)來表示狀態(tài)空間模型的狀態(tài)轉(zhuǎn)移矩陣可以寫為

(12)

H=[1 0 0 … 0m]，m=max(p,q)。其中，X1t的初始值為{y1y2…ym}。X(t)其他元素通過遞推公式Xmt=φmYt-1-θm-1at+1，最后得到狀態(tài)初始值X1|0=(s1ts2t…smt)T和觀測(cè)初始值Yt=(y1t22t…ymt)T，根據(jù)式(8)、式(9)得到最優(yōu)估計(jì)狀態(tài)向量Xt|t。

Step 3采用MATLAB軟件，對(duì)ARIMA模型進(jìn)行單步預(yù)測(cè)。

Step 4構(gòu)建狀態(tài)空間模型的觀測(cè)方程。

ARIMA模型和組合模型單步預(yù)測(cè)值作為卡爾曼濾波的觀測(cè)量，在已知Xt|t的情況下，狀態(tài)空間模型參數(shù)估計(jì)后，就能夠得到t+1時(shí)刻的預(yù)測(cè)。

Yt+1|t=E(Yt+1|Yt)=HXt+1|t+E(vt+1)=

HFXt|t

(13)

3.2 ARIMA-KF外推預(yù)測(cè)方法

由于卡爾曼濾波無法單獨(dú)進(jìn)行向后預(yù)測(cè)，因此需要借助模型擬合的預(yù)測(cè)數(shù)據(jù)進(jìn)行卡爾曼濾波分析，在模型預(yù)測(cè)的基礎(chǔ)上提供新的觀測(cè)值，利用卡爾曼濾波強(qiáng)大的迭代功能對(duì)狀態(tài)向量估計(jì)不斷修正進(jìn)而得到最優(yōu)估計(jì)，提高模型的預(yù)測(cè)精度。建模流程如圖1所示。

圖1 建模流程圖Fig.1 Modeling flow chart

4 實(shí)例分析

4.1 基于ARIMA-KF的設(shè)備主要狀態(tài)參數(shù)預(yù)測(cè)仿真

以小時(shí)為單位，將采集到312 h的實(shí)船海水冷卻系統(tǒng)淡水出口溫度值作為原始數(shù)據(jù)樣本。利用MATLAB軟件，將前300個(gè)數(shù)據(jù)作為訓(xùn)練樣本輸入預(yù)測(cè)模型中，預(yù)測(cè)未來12 h的溫度值。將預(yù)測(cè)值和采集到的實(shí)船測(cè)量值進(jìn)行誤差對(duì)比分析，具體步驟內(nèi)容如下。

Step 1確保在船舶正常運(yùn)行的情況下，對(duì)實(shí)船冷卻水系統(tǒng)淡水出口的溫度進(jìn)行采集，采集到的原始樣本數(shù)據(jù)如圖2所示。

圖2 信號(hào)采集原始圖Fig.2 Original diagram of signal acquisition

Step 2采用MATLAB軟件進(jìn)行數(shù)據(jù)處理和編程，利用KPSS函數(shù)和ADF函數(shù)檢驗(yàn)數(shù)據(jù)的平穩(wěn)性，當(dāng)KSPP=0,ADF=1時(shí)說明被檢驗(yàn)序列為平穩(wěn)序列。運(yùn)行后的結(jié)果顯示原始序列需經(jīng)過2階差分才轉(zhuǎn)換為平穩(wěn)序列，故該模型初定為ARIMA(p,2,q)模型。該模型的自相關(guān)函數(shù)圖形與偏自相關(guān)函數(shù)圖形如圖3所示。

圖3 自相關(guān)函數(shù)、偏自相關(guān)函數(shù)圖像Fig.3 Autocorrelation function、Partial autocorrelation function image

Step 3確定階數(shù)p、q。根據(jù)表1和圖3可知模型的自相關(guān)函數(shù)與偏自相關(guān)函數(shù)均為拖尾性。初步將階數(shù)的范圍定為5≤p≤12，5≤q≤12。根據(jù)式(6)計(jì)算不同組合階數(shù)對(duì)應(yīng)的AIC、BIC值，依照AIC準(zhǔn)則、BIC準(zhǔn)則，當(dāng)p=12、q=5時(shí)，AIC、BIC值相對(duì)最小，又d=2，故時(shí)間序列模型為ARIMA(12,2，5)。接下來對(duì)該模型的擬合殘差進(jìn)行檢驗(yàn)。結(jié)果如圖4所示，可知，擬合殘差均值為0，且服從正態(tài)分布，樣本正態(tài)分位數(shù)分布在紅線的周圍，滿足理想模型對(duì)擬合殘差的要求。

圖4 模型檢測(cè)結(jié)果圖Fig.4 Model test results

Step 4在已知p=12，d=2，q=5。由式(4)得到ARMA模型的預(yù)測(cè)方程為

x(t+1)=1.983x(t)+…-0.153x(t-13)-α(t+1)-…-0.189α(t-4)

(14)

根據(jù)式(14)轉(zhuǎn)化為狀態(tài)空間模型的狀態(tài)方程和觀測(cè)方程(令H=[1,0,0,0,0,0,0,0,0,0,0,0])分別為

(15)

ω(t+1)

(16)

式中：α(t)、ω(t)均為白噪聲序列。一般取狀態(tài)初始值X1|0=(s1ts2t…smt)T=(y1ty2t…ymt)T，觀測(cè)初始值p(0|0)=10I，I為多維的單位矩陣，應(yīng)用MATLAB軟件根據(jù)式(8)～式(10)實(shí)現(xiàn)卡爾曼濾波優(yōu)化和遞推預(yù)測(cè)。

4.2 仿真結(jié)果與分析

為進(jìn)一步驗(yàn)證組合預(yù)測(cè)模型的有效性，選取圖2前300個(gè)實(shí)測(cè)數(shù)據(jù)作為組合預(yù)測(cè)模型的輸入，采用MATLAB軟件編程預(yù)測(cè)未來12 h的數(shù)據(jù)，和實(shí)測(cè)數(shù)據(jù)以及單一的ARIMA模型的預(yù)測(cè)值進(jìn)行誤差對(duì)比分析。如圖5～圖7所示，從整體上看，無論是單一的ARIMA模型還是ARIMA-KF組合模型，預(yù)測(cè)結(jié)果都是呈現(xiàn)線性下降趨勢(shì)，符合真實(shí)值的變化趨勢(shì)，但ARIMA-KF組合模型較單一的ARIMA模型明顯更接近實(shí)測(cè)值。圖5為ARIMA模型和ARIMA-KF組合模型的多步預(yù)測(cè)結(jié)果表明，組合模型的單步預(yù)測(cè)結(jié)果較單一的ARIMA模型預(yù)測(cè)結(jié)果的準(zhǔn)確性高。圖6、圖7是兩種模型對(duì)實(shí)測(cè)數(shù)據(jù)的單步預(yù)測(cè)整體圖和分時(shí)圖，從中可看出：組合模型的單步預(yù)測(cè)基本與是實(shí)測(cè)數(shù)據(jù)保持一致，能夠較好地反映冷卻水淡水出口溫度的變化趨勢(shì)。從圖5和圖6、圖7對(duì)比來看，無論是ARIMA-KF組合模型還是單一ARIMA預(yù)測(cè)模型，單步預(yù)測(cè)的準(zhǔn)確度更高。

圖5 多步預(yù)測(cè)結(jié)果對(duì)比Fig.5 Multi step predictionresult comparison

圖6 單步預(yù)測(cè)結(jié)果整體對(duì)比Fig.6 Overallcomparison of single step prediction result

圖7 單步預(yù)測(cè)結(jié)果分時(shí)對(duì)比Fig.7 Time-sharing comparison of single-step prediction results

4.3 模型評(píng)價(jià)

為進(jìn)一步分析不同預(yù)測(cè)方法的預(yù)測(cè)效果，本文對(duì)每種模型算法進(jìn)行了均方誤差(MSE)、平均絕對(duì)誤差(MAE)、均方根誤差(RMSE)計(jì)算分析。最終結(jié)果如表2所示。

從圖7和表2中可以看出，單一的ARIMA模型多步預(yù)測(cè)誤差是最大。ARIMA-KF組合模型的多步預(yù)測(cè)誤差和ARIMA模型的單步預(yù)測(cè)誤差結(jié)果差別很小，但均比單一的ARIMA模型預(yù)測(cè)效果好。ARIMA-KF組合模型單步預(yù)測(cè)誤差最低，預(yù)測(cè)效果最好。

5 結(jié)論

介紹了ARIMA模型單步預(yù)測(cè)和多步預(yù)測(cè)方法，分析了現(xiàn)有船舶設(shè)備狀態(tài)參數(shù)預(yù)測(cè)模型的不足之處，利用卡爾曼濾波算法強(qiáng)大迭代能力，對(duì)ARIMA模型的參數(shù)進(jìn)行尋優(yōu)，彌補(bǔ)了單一算法的不足的同時(shí)降低了建模難度。該算法對(duì)監(jiān)測(cè)信號(hào)樣本少、呈線性的非平穩(wěn)船舶系統(tǒng)及設(shè)備狀態(tài)參數(shù)進(jìn)行預(yù)測(cè)研究，將外界的復(fù)雜影響因素歸為單一的時(shí)間因素來考慮。通過對(duì)實(shí)船海水冷卻系統(tǒng)設(shè)備的淡水出口溫度進(jìn)行單步和多步預(yù)測(cè)，驗(yàn)證了模型有效性。將組合模型與單一的ARIMA模型的預(yù)測(cè)值進(jìn)行對(duì)比，結(jié)果表明，組合模型的預(yù)測(cè)精度更高，預(yù)測(cè)效果更為理想。為發(fā)展船舶系統(tǒng)設(shè)備的智能維修決策提供了理論支持。