江福才 范慶波 馬全黨 張 帆 馬 勇

(武漢理工大學(xué)航運(yùn)學(xué)院1) 武漢 430063) (武漢理工大學(xué)內(nèi)河航運(yùn)技術(shù)湖北省重點(diǎn)實(shí)驗(yàn)室2) 武漢 430063)

0 引 言

隨著水上交通在智能化、信息化等方面的發(fā)展，實(shí)時(shí)、準(zhǔn)確的交通流量預(yù)測可為水上智能交通系統(tǒng)進(jìn)行交通誘導(dǎo)及控制提供有效的基礎(chǔ)數(shù)據(jù)，便于海事主管機(jī)關(guān)制定合理、有效的調(diào)度、交通組織方案.同時(shí)港口部門可根據(jù)船舶交通流量的變化趨勢，合理規(guī)劃碼頭、航道、錨地等水域的布局，提高港口資源的利用效率[1].

目前國內(nèi)外的學(xué)者在水上交通流預(yù)測方面取得了許多研究成果，但大多是借鑒陸地、航空等領(lǐng)域的預(yù)測方法，比較常用的方法有支持向量機(jī)(support vector machine，SVM)模型、灰色模型、時(shí)間序列模型、神經(jīng)網(wǎng)路模型、組合模型等[2-5]，這些模型有各自的優(yōu)勢，但很少考慮序列的非平穩(wěn)性，因此，預(yù)測精度相對較小.在進(jìn)行月、季、年度交通流量的預(yù)測時(shí)，往往受到趨勢性、季節(jié)性等因素的影響，季節(jié)性自回歸移動平均(seasonal auto regressive integrated moving average，SARIMA)模型可以很好的解決這些問題.在SARIMA模型研究方面，李曉磊等[6]運(yùn)用SARIMA模型對荊州港月度船舶交通流量進(jìn)行了預(yù)測，并將該方法與Holt-Winters加法模型、自回歸移動平均(auto regressive integrated moving average，ARIMA)模型進(jìn)行了對比，驗(yàn)證了SARIMA模型的可行性；郝軍章等[7]將SARIMA模型運(yùn)用到全國鐵路客運(yùn)量預(yù)測中，通過將預(yù)測值與實(shí)際值進(jìn)行對比，得到模型的預(yù)測能力為“優(yōu)良”；王瑩等[8]選擇SARIMA模型對北京地鐵進(jìn)站客流量進(jìn)行時(shí)間序列建模，從預(yù)測結(jié)果看，該模型可以較準(zhǔn)確的描述進(jìn)站客流量的變化趨勢.

馬爾科夫(Markov)模型具有無后效性，廣泛應(yīng)用到交通流量、瓦斯流量、交通事故等預(yù)測方面.Yang等[9]運(yùn)用Markov模型與小波神經(jīng)網(wǎng)絡(luò)模型結(jié)合，深入洞察樣本流量模式，提高預(yù)測的準(zhǔn)確性與魯棒性，對德克薩斯州休斯敦的高速公路流量進(jìn)行了預(yù)測；Lyu等[10]運(yùn)用BP神經(jīng)網(wǎng)絡(luò)與Markov模型組合，對防城港船舶交通流量進(jìn)行了預(yù)測；陳釗等[11]運(yùn)用GM(1,1)模型與Markov模型結(jié)合對礦井瓦斯流量進(jìn)行了預(yù)測等等.因傳統(tǒng)的Markov模型在求取灰色區(qū)間的白化系數(shù)時(shí)存在不足，本文選用具備較強(qiáng)全局搜索能力的粒子群算法求解最佳的白化系數(shù)，在粒子群算法優(yōu)化方面，樊冬雪[12]利用無偏灰色馬爾科夫模型對鐵路客流量進(jìn)行預(yù)測，并利用粒子群算法(particle swarm optimization，PSO)輔助求解馬爾科夫模型中最佳的白化系數(shù)；Eseye[13]運(yùn)用PSO進(jìn)行全局搜索，對支持向量機(jī)中的參數(shù)進(jìn)行優(yōu)化，以提高預(yù)測精度等.

鑒于以上分析，本文將SARIMA-Markov模型引入船舶交通流量預(yù)測領(lǐng)域，發(fā)揮了SARIMA模型與Markov模型所具有的優(yōu)勢，并且運(yùn)用PSO進(jìn)行優(yōu)化，求取最優(yōu)的白化系數(shù)，以期更好的表征船舶交通流量的發(fā)展趨勢.

1 SARIMA-Markov模型基本原理

1.1 SARIMA預(yù)測模型

ARIMA模型(B-J方法)，常用于解決隨機(jī)時(shí)序問題[14].SARIMA模型是由ARIMA模型衍生而來的，基本形式為SARIMA(p,d,q)(P,D,Q)s.

設(shè)原始時(shí)間序列是它的當(dāng)期和前期的隨機(jī)誤差項(xiàng)及前期值的線性函數(shù)，即可表示為

Xi=φ1xi-1+φ2xi-2+…+φpxi-p+

μi-θ1μi-1-θ2μi-2-…-θqμi-q

(1)

式中：時(shí)間序列Xi為自回歸移動平均序列；ui為相互獨(dú)立的白噪聲序列.上述為(p，q)階的自回歸移動平均模型，記為ARMA(p，q)；φ1,φ2,…,φp為自回歸系數(shù)；θ1,θ2,…,θq為移動平均系數(shù).

引入滯后算子，式(1)可簡記為

φ(B)xi=θ(B)μi

(2)

若序列在季節(jié)性差分前還進(jìn)行了d階的逐期差分才平穩(wěn)，則最終形成的SARIMA模型表達(dá)式為

(3)

軟件常被用于解決時(shí)間序列建模及數(shù)據(jù)處理問題[15].

1.2 Markov預(yù)測模型

Markov模型具有無后效性，對于一個(gè)具有無后效性的隨機(jī)序列{Y(t),t∈T}，T為離散的時(shí)間集合，其無后效性可以表示為

P{Yn+1=in+1|Y1=i1,Y2=i2,…，Yn=in}=

P{Yn+1=in+1|Yn=in}

(4)

Markov預(yù)測模型的構(gòu)建可以分為以下幾個(gè)步驟.

步驟1求取殘差序列e(i) 對原始數(shù)據(jù)Xi，運(yùn)用SARIMA模型進(jìn)行計(jì)算得到預(yù)測值fi，同時(shí)可以得到殘差序列e(i)

e(i)=xi-fi，i=0,1,2,…

(5)

步驟2狀態(tài)區(qū)間劃分 將殘差序列按照大小劃分為間隔寬度相等的r個(gè)狀態(tài)，殘差序列e(i)的第i步的第j個(gè)狀態(tài)的上下邊界分別用Qij和Hij表示，計(jì)算方法為

(6)

步驟3求取k步狀態(tài)轉(zhuǎn)移概率 設(shè)狀態(tài)空間為I={i1,i2,…，ir}，k步狀態(tài)轉(zhuǎn)移概率為

{ii→i1,ii→i2,…，ii→ir}

(7)

步驟4構(gòu)建狀態(tài)轉(zhuǎn)移概率矩陣 構(gòu)建能反映系統(tǒng)間各狀態(tài)概率轉(zhuǎn)移的一步轉(zhuǎn)移矩陣，其大小為r·r.該轉(zhuǎn)移矩陣描述為

(8)

步驟5求取預(yù)測值 Markov模型的預(yù)測值表示為

(9)

1.3 粒子群算法

選用粒子群算法進(jìn)行全局搜索來確定Markov算法中的白化系數(shù)λi(i=1,2,…，r)[16-17].

對于某個(gè)粒子i，其粒子群搜索公式為

vi(t+1)=wvi(t)+c1rand1(xpbesti-xi(t))+

c2rand2(xgbesti-xi(t))

(10)

xi(t+1)=xi(t)+vi(t+1)

式中：xi為在t時(shí)刻的位置；vi為粒子在t時(shí)刻的速度；xpbest為群體微粒搜索到的自身最佳位置；xgbest為群體微粒經(jīng)歷過的全局最好位置；rand1，rand2為服從U(0，1)的隨機(jī)數(shù)；w為慣性系數(shù)；c1，c2為加速度系數(shù).

用均方差fit來衡量每個(gè)粒子的適應(yīng)度函數(shù)

(11)

選擇軟件編寫相關(guān)程序求得最佳的白化系數(shù)λi，最終得到狀態(tài)區(qū)間中可能的預(yù)測值zi為

zi=(1-λi)Qij+λiHijλi∈[0,1]

(12)

2 仿真與數(shù)據(jù)分析

2.1 數(shù)據(jù)樣本的選取

為了更好的表征船舶交通狀況，實(shí)現(xiàn)對不同水域相同時(shí)段或同一水域不同時(shí)段的船舶交通流量的可比性，提出了標(biāo)準(zhǔn)船舶的概念[18].我國在評估海上交通安全時(shí)提出的標(biāo)準(zhǔn)船舶換算系數(shù)見表1.

表1 標(biāo)準(zhǔn)船舶換算系數(shù)表

運(yùn)用以上換算系數(shù)對原始數(shù)據(jù)的船舶交通流量進(jìn)行加權(quán)，轉(zhuǎn)化為標(biāo)準(zhǔn)船舶，原始數(shù)據(jù)為從海事相關(guān)部門調(diào)研得到的赤壁長江公路大橋斷面2010年1月—2016年6月的船舶交通流量數(shù)據(jù)，經(jīng)過標(biāo)準(zhǔn)船舶轉(zhuǎn)化后的樣本數(shù)據(jù)見表2.

表2 赤壁長江公路大橋斷面標(biāo)準(zhǔn)船舶交通流量統(tǒng)計(jì)表

2.2 月度標(biāo)準(zhǔn)船舶交通流量預(yù)測

為驗(yàn)證SARIMA-Markov模型的適用性，運(yùn)用、軟件，選取赤壁長江公路大橋斷面2010年1月—2015年12月(編號1～36)的樣本數(shù)據(jù)用于建立模型，2016年1—6月(編號37～42)的樣本數(shù)據(jù)用于模型驗(yàn)證，并運(yùn)用相關(guān)評價(jià)指標(biāo)對模型進(jìn)行評價(jià).

2.2.1SARIMA模型

依據(jù)前面介紹的建模步驟，首先對序列進(jìn)行平穩(wěn)性處理，見圖1.各觀測值沒有完全圍繞其均值上下波動，曲線前后波動具有不一致性；表3為原始數(shù)據(jù)、一階差分后ACF-PACF分析圖，由圖3可知，原始數(shù)據(jù)的自相關(guān)系數(shù)(ACF)起初都落在顯著水平α=0.05的置信區(qū)間以外，綜上所述，得到序列xi為非平穩(wěn)時(shí)間序列.

對非平穩(wěn)序列通常采用逐期差分的方法，ACF很快的趨于0(滯后階數(shù)k>2或3時(shí))，落入置信區(qū)間內(nèi)；對序列xi進(jìn)行單位根檢驗(yàn)，結(jié)果見表4.顯著水平為1%的臨界值大于t統(tǒng)計(jì)值，得到序列xi轉(zhuǎn)化為平穩(wěn)序列，因此d=1.

圖1 時(shí)間序列圖(2010年1月—2015年12月)

由表3可知，ACF在12階，24階等以12為周期的階數(shù)位置存在較大的峰值，表明序列存在以12個(gè)月為周期的季節(jié)波動，為消除模型的季節(jié)性，本文應(yīng)用軟件中的seasonal adjustment-census X12功能.參照原始序列平穩(wěn)性的分析步驟，得到D=1.

ACF-PACF分析圖常被用來對SARIMA(p,d,q)(P,D,Q)s模型進(jìn)行識別和定階.根據(jù)以上分析，得到d=1，D=1，見表5，由表5可知，ACF在k=12,24處顯著超出置信區(qū)間，存在以12為周期的波動性，得到P=1、Q=1.參照相關(guān)文獻(xiàn)、書籍對p、q值的取值方法[19-20]，對比表3、表5可知，得到可供選擇的(p、q)組合有(3,1)、(4,0)、(2,1)、(3,0).

表3原始數(shù)據(jù)、一階差分后ACF-PACF分析圖

表4 一階差分后ADF檢驗(yàn)圖

根據(jù)不同的(p、q)組合建立SARIMA(3,1,1)(1,1,1)12，SARIMA(4,1,0)(1,1,1)12，SARIMA(2,1,1)(1,1,1)12，SARIMA(3,1,0)(1,1,1)12四組模型，運(yùn)用軟件中的Equation Estimation功能對模型進(jìn)行檢驗(yàn)，得到選取SARIMA(3,1,0)(1,1,1)12是比較合適的.

SARIMA(3,1,0)(1,1,1)12模型表達(dá)式為

(1+B12)(1+0.621B+0.376B2+0.138B3)

(1-B)(1-B12)lg(xi)=(1-0.009B12})μi

表5一階季節(jié)差分后ACF-PACF分析圖

運(yùn)用SARIMA(3,1,0)(1,1,1)12模型對2013年1月—2015年12月的月度標(biāo)準(zhǔn)船舶交通流量進(jìn)行預(yù)測，結(jié)果見表6.

2.2.2SARIMA-Markov模型

根據(jù)Markov模型的建模步驟，首先確定殘差序列的范圍e(i)∈[-1 455,1 442]，將殘差序列劃分為4個(gè)狀態(tài)區(qū)間，區(qū)間的上下邊界值為

表6 SARIMA模型預(yù)測值及殘差序列統(tǒng)計(jì)表 艘·次

同理，E2、E3、E4狀態(tài)的上下邊界區(qū)間分別為[-730.75,-6.5]，[-6.5,717.75]，[717.75,1 442].殘差序列狀態(tài)分布見表7.

表7 狀態(tài)分布表

根據(jù)表7中狀態(tài)劃分情況，對從Ei經(jīng)一步轉(zhuǎn)移到Ej的頻數(shù)統(tǒng)計(jì)見表8.

表8 一步狀態(tài)轉(zhuǎn)移頻數(shù)統(tǒng)計(jì)表

狀態(tài)轉(zhuǎn)移概率：P11=P(E1→E1)=P(E1|E1)=0.4

同理可得到：P12=0.3，P13=0.2，P14=0.1，P21=P22=0.364，P23=0.091，P24=0.182，P31=0.2,P32=0.2，P33=0.3，P34=0.1，P43=0.75，P44=0.25.

狀態(tài)轉(zhuǎn)移概率矩為

運(yùn)用編制粒子群算法程序，進(jìn)行全局搜索，以找到最佳的白化系數(shù)，計(jì)算結(jié)果為

λ1=0.015 4，λ2=0.142 0，

λ3=0.327 8，λ4=0.988 4

根據(jù)式(12)，得到：z1=-741.90、z2=-109.34、z3=480.34、z4=726.15.根據(jù)式(9)計(jì)算得到SARIMA-Markov模型的預(yù)測值，見表9和圖2.

表9 SARIMA-Markov模型預(yù)測值及殘差序列統(tǒng)計(jì)表 艘·次

圖2 SARIMA-Markov模型擬合示意圖

由圖2可知，觀察兩種模型的擬合曲線，SARIMA-Markov模型的擬合效果要優(yōu)于SARIMA模型.

2.2.3模型驗(yàn)證

為對所建模型的預(yù)測精度進(jìn)行驗(yàn)證，擬運(yùn)用構(gòu)建好的模型對赤壁長江公路大橋斷面2016年1—6月的月度標(biāo)準(zhǔn)船舶交通流量進(jìn)行預(yù)測.

以第37、38組數(shù)據(jù)為例，運(yùn)用軟件，計(jì)算SARIMA(3,1,0)(1,1,1)12模型2016年1月、2月的預(yù)測值為6 015,4 213/(艘·次)，然后根據(jù)式(8)～(9)計(jì)算得到第36組數(shù)據(jù)轉(zhuǎn)移到37、38組數(shù)據(jù)的1步、2步狀態(tài)轉(zhuǎn)移概率矩陣，結(jié)果為

第36組數(shù)據(jù)處于E4，得到第37組數(shù)據(jù)的預(yù)測值為：5 800/(艘次)，38組數(shù)據(jù)的預(yù)測值為：6 009/艘次.

同理可得到后面幾個(gè)月的預(yù)測值，結(jié)果見表10.

表10 預(yù)測結(jié)果統(tǒng)計(jì)表 艘·次

本文擬采用平均絕對誤差、平均絕對百分誤差對模型擬合效果進(jìn)行評價(jià)，見表11.優(yōu)化后模型的兩項(xiàng)指標(biāo)均小于優(yōu)化前，SARIMA-Markov模型的擬合精度達(dá)到了92.0849%，優(yōu)于SARIMA模型84.755 2%的擬合精度.SARIMA-Markov模型預(yù)測精度為95.786 1%，優(yōu)于單純利用SARIMA模型87.445 4%的預(yù)測精度，SARIMA-Markov模型預(yù)測精度較高，滿足要求.

表11 模型檢驗(yàn)表

3 結(jié) 論

1) 將SARIMA-Markov模型應(yīng)用到月度船舶交通流量預(yù)測，樣本數(shù)據(jù)為經(jīng)過加權(quán)后的標(biāo)準(zhǔn)船舶交通流量.該模型綜合考慮了SARIMA模型對月度船舶交通流量季節(jié)相關(guān)性的表達(dá)，Markov模型運(yùn)用殘差序列、狀態(tài)劃分、狀態(tài)轉(zhuǎn)移矩陣等方法處理樣本近期狀況對預(yù)測值的影響，并運(yùn)用具有全局搜索能力的粒子群算法求取模型中的最佳白化系數(shù)，為船舶交通流量預(yù)測提供了新方法.

2) 本文模型與單純利用SARIMA模型相比，能確保更高精度的預(yù)測結(jié)果，驗(yàn)證了該模型的科學(xué)性、可行性.

3) 對于該模型的可推廣性及最佳狀態(tài)劃分方案等，在后續(xù)的研究工作仍有待進(jìn)一步完善.

4) 預(yù)測結(jié)果可為港口水域的合理布局、水上智能交通系統(tǒng)、海事主管機(jī)關(guān)等提供基礎(chǔ)性數(shù)據(jù)，該方法對受季節(jié)性影響的港口、航道等水域船舶交通流量預(yù)測也有很好的借鑒意義.