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

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

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

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

針對上述模型的局限性，現(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ù)測方法。該模型將兩種算法的優(yōu)勢結(jié)合在一起，能彌補單一模型的不足，提高模型的準(zhǔn)確度。此外，該模型將外界復(fù)雜的影響因素歸為單一的時間影響因素來考慮，可降低模型構(gòu)建的難度[10]，并依靠卡爾曼濾波算法強大的更新迭代能力，大大減少了工作量[3]。該模型對船舶系統(tǒng)狀態(tài)監(jiān)控和輔助決策具有一定的指導(dǎo)意義。

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

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

(1)

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

(2)

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

建模具體過程如下。

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

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

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

(3)

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

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

表1 模型判定方法

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

(4)

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

2 卡爾曼濾波

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

狀態(tài)方程為

Xt=FtXt-1+ct+Rtwt

(5)

量測方程為

Yt+1=ZtXt+dt+vt

(6)

式中：Xt為t時刻的狀態(tài)估計值(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時刻k×1維的觀測向量；Zt為k×1維量程矩陣；dt是k×1維的干擾項；wt、vt代表的是均值為0、協(xié)方差矩陣分別為Gt、Qt的干擾項，即E(wt)=0，var(wt)=Qt。E(vt)=0，var(vt)=Gt，且Gt、Qt相互獨立。

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

狀態(tài)方程為

Xt=FtXt-1+Btut+wt

(7)

測量方程為

Yt=HtXt+vt

(8)

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

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

(9)

觀測量預(yù)測步為

(10)

式(9)中：

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

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

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

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

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

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

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

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ù)測方程來表示狀態(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和觀測初始值Yt=(y1t22t…ymt)T，根據(jù)式(8)、式(9)得到最優(yōu)估計狀態(tài)向量Xt|t。

Step 3采用MATLAB軟件，對ARIMA模型進行單步預(yù)測。

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

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

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

HFXt|t

(13)

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

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

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

4 實例分析

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

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

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

圖2 信號采集原始圖Fig.2 Original diagram of signal acquisition

Step 2采用MATLAB軟件進行數(shù)據(jù)處理和編程，利用KPSS函數(shù)和ADF函數(shù)檢驗數(shù)據(jù)的平穩(wěn)性，當(dāng)KSPP=0,ADF=1時說明被檢驗序列為平穩(wě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)計算不同組合階數(shù)對應(yīng)的AIC、BIC值，依照AIC準(zhǔn)則、BIC準(zhǔn)則，當(dāng)p=12、q=5時，AIC、BIC值相對最小，又d=2，故時間序列模型為ARIMA(12,2，5)。接下來對該模型的擬合殘差進行檢驗。結(jié)果如圖4所示，可知，擬合殘差均值為0，且服從正態(tài)分布，樣本正態(tài)分位數(shù)分布在紅線的周圍，滿足理想模型對擬合殘差的要求。

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

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

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)方程和觀測方程(令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，觀測初始值p(0|0)=10I，I為多維的單位矩陣，應(yīng)用MATLAB軟件根據(jù)式(8)～式(10)實現(xiàn)卡爾曼濾波優(yōu)化和遞推預(yù)測。

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

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

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

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

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

4.3 模型評價

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

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

5 結(jié)論

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