王當(dāng)利 王雪佳 呂 雪 楊馨穎

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

0 引 言

水上交通事故預(yù)測是水上交通安全領(lǐng)域的重點研究內(nèi)容之一，從技術(shù)上分為定性預(yù)測與定量預(yù)測兩種，而定量預(yù)測方法中又包括單一預(yù)測和組合預(yù)測.目前，國內(nèi)外關(guān)于交通事故預(yù)測使用的單一預(yù)測方法主要有灰色理論、時間序列、回歸分析、神經(jīng)網(wǎng)絡(luò)、馬爾科夫鏈等方法[1-2].Debnath等[3]基于時間序列預(yù)測方法對新加坡港航道內(nèi)碰撞事故進(jìn)行定量預(yù)測；Rahman[4]運用貝葉斯網(wǎng)絡(luò)分析水上交通事故影響因素與事故概率之間的關(guān)系，并對不同類別的事故數(shù)量進(jìn)行有效預(yù)測.這些方法均為水上交通事故預(yù)測提供了研究思路，但由于模型自身特點導(dǎo)致預(yù)測精度有限.組合預(yù)測具有能夠有效降低預(yù)測誤差和綜合利用單一預(yù)測方法的優(yōu)勢，但目前關(guān)于水上交通事故組合預(yù)測的研究并不多，現(xiàn)有的研究也僅限于灰色預(yù)測模型和其他預(yù)測模型結(jié)合的組合預(yù)測，李鈴鈴等[5]首次將灰色預(yù)測模型的預(yù)測值殘差運用BP神經(jīng)網(wǎng)絡(luò)加以修正，進(jìn)而得到水上交通事故量的組合預(yù)測值；趙佳妮等[6]運用馬爾可夫預(yù)測模型對灰色GM(1,1)預(yù)測模型的預(yù)測結(jié)果進(jìn)行優(yōu)化，對具有較強(qiáng)波動特性的數(shù)據(jù)能得到理想預(yù)測結(jié)果.與水上交通事故預(yù)測相比，道路交通對于組合預(yù)測研究更為深入，除了基于傳統(tǒng)灰色預(yù)測方法的組合預(yù)測外，有基于ARIMA模型和支持向量回歸機(jī)模型的時間序列組合預(yù)測、最優(yōu)加權(quán)組合預(yù)測、D-S證據(jù)理論組合預(yù)測等模型.可見，嘗試?yán)闷渌问降慕M合預(yù)測方法是水上交通事故預(yù)測研究的新思路.

水上交通系統(tǒng)具有典型的非線性、灰色和動態(tài)特性，水上交通事故的影響因素眾多且影響因素間相互影響、互相關(guān)聯(lián)；另外，受自然條件等外部環(huán)境的影響，水上交通事故在不同月份的發(fā)生量有所不同，具有明顯的季節(jié)變化趨勢.考慮到多元灰色預(yù)測模型MGM(1,n)能夠綜合考慮多個相關(guān)因素對事故發(fā)生的影響作用，季節(jié)差分自回歸滑動平均模型(seasonal autoregressive integrated moving average, SARIMA)能夠集合事故數(shù)據(jù)在周期間和周期內(nèi)的變化特征，更全面、真實的反映事故發(fā)生趨勢，本文將選用這兩種模型對事故數(shù)進(jìn)行組合預(yù)測[7-8].同時，為了突破原有組合預(yù)測的局限，改變各單預(yù)測模型權(quán)系數(shù)在待測時間點恒定不變的限制，引入誘導(dǎo)有序加權(quán)平均(induced ordered weighted averaging operate,IOWA)算子，以提高預(yù)測性能.以天津水域事故統(tǒng)計為原始數(shù)據(jù)，對比各預(yù)測模型的評價結(jié)果，驗證該模型的適用性及合理性，并對天津水域2018-2020年的事故數(shù)進(jìn)行預(yù)測，得到近3年水域內(nèi)事故發(fā)生的趨勢.

1 組合預(yù)測模型的構(gòu)建

1.1 基于背景值優(yōu)化的多元灰色預(yù)測模型MGM(1,n)

1) 確定水上交通事故原始數(shù)據(jù) 根據(jù)關(guān)聯(lián)度分析，確定水上交通事故原始數(shù)據(jù)序列.

(1)

式中：n為系統(tǒng)中變量個數(shù)；m為時間點個數(shù).

2) 獲取累加生成序列 為了弱化事故數(shù)據(jù)的隨機(jī)性，使其呈現(xiàn)出一定規(guī)律，對原始數(shù)據(jù)序列進(jìn)行一次累加處理，破壞歷史事故數(shù)據(jù)的離亂特性.一次累加生成序列為

(2)

k=1,2,…,m；i=1,2,…,n.

3) 建立n元一階常微分方程 根據(jù)生成數(shù)列建立預(yù)測模型的一階常微分方程組，并轉(zhuǎn)換成矩陣形式，為

(3)

4) 求解模型中的參數(shù)A和B將式(2)進(jìn)行離散化處理，利用最小二乘法得：

(4)

k=1,2,…,m；i=1,2,…,n.

(5)

H=(BTB)-1BTY

(6)

式中：

(7)

采用非齊次指數(shù)函數(shù)對MGM(1,n)模型一次累加生成序列進(jìn)行擬合，得到優(yōu)化后的背景值為

(8)

從式(7)的矩陣H中可以得到A和B的辨識值和

(9)

5) 建立時間響應(yīng)函數(shù) 將A的辨識值帶入所構(gòu)造的微分方程組內(nèi)，同時假定邊界條件為然后按一般微分方程組求解的過程進(jìn)行求解，由此得到時間響應(yīng)函數(shù)為

X(1)(t)=eAtX(1)(0)+A-1B-A-1B

(10)

6) 時間響應(yīng)函數(shù)離散化 將X(0)(1)作為第一個數(shù)據(jù)的前提，經(jīng)過k-1個間隔到達(dá)X(0)(k)，得到：

k=2,3,…,n.

(11)

1.2 SARIMA時間序列預(yù)測模型

SARIMA模型即SARIMA(p,d,q)(P,D,Q)s模型，一般形式為

(12)

式中：S為季節(jié)性周期長度；φp(Bs)和ΗQ(Bs)ut分別為季節(jié)P階自回歸算子和Q階移動平均算子；P為季節(jié)性自回歸階數(shù)；Q為季節(jié)性移動平均階數(shù)；D為季節(jié)性差分階數(shù).當(dāng)P,D,Q的值均為0時，此時的模型即為ARIMA模型.模型的具體建模步驟為：

步驟1平穩(wěn)性檢驗 通過分析時序圖和序列相關(guān)圖，判斷序列的平穩(wěn)性.當(dāng)序列出現(xiàn)周期性的變化趨勢時，需要對其平穩(wěn)性進(jìn)行檢驗.根據(jù)單位根(augmented dickry-fuller，ADF)檢驗的結(jié)果判定待測事故數(shù)據(jù)序列是否平穩(wěn)，若序列出現(xiàn)明顯的波動性，則需要通過差分和季節(jié)差分處理，使之變?yōu)槠椒€(wěn)序列.

步驟2模型的識別和參數(shù)估計 模型識別過程即通過分析自相關(guān)和偏相關(guān)函數(shù)，確定SARIMA模型階數(shù)p，d，q及P，D，Q進(jìn)而選取預(yù)測效果最好的預(yù)測模型的過程.當(dāng)原始事故序列經(jīng)過 階差分和步長w的n階季節(jié)差分后平穩(wěn)，那么d=m，D=n，s=w.通過分析自相關(guān)函數(shù)和偏相關(guān)函數(shù)，確定模型的階數(shù)，通過設(shè)定不同的預(yù)測模型，并根據(jù)AIC準(zhǔn)則，確定最佳的SARIMA模型.

步驟3白噪聲檢驗 對SARIMA預(yù)測模型進(jìn)行檢驗時，需要判斷殘差序列的獨立性.當(dāng)殘差序列是白噪聲序列時，通過序列t檢驗或殘差的Q統(tǒng)計量判斷模型的合理性.白噪聲序列是不存在相關(guān)關(guān)系的序列.

步驟4模型預(yù)測 利用Eviews8.0軟件對事故數(shù)據(jù)進(jìn)行預(yù)測，根據(jù)上述步驟建立預(yù)測模型，實現(xiàn)對樣本內(nèi)事故數(shù)據(jù)序列的靜態(tài)預(yù)測和樣本外事故數(shù)據(jù)序列的動態(tài)預(yù)測.

1.3 基于IOWA算子的組合預(yù)測模型

1.3.1IOWA算子

設(shè)某實際情況的指標(biāo)序列值為xt(t=1,2,…,N)，預(yù)測時利用m個單預(yù)測模型，xit為第i個單預(yù)測模型在t時刻根據(jù)預(yù)測模型計算得出的預(yù)測值，i=1,2,…,m,t=1,2,…,N.ait為第i個單項預(yù)測模型在第t時刻的預(yù)測精度，同時a的值在0與1的范圍之間.可將ait視為xit的誘導(dǎo)值，第t時刻m個單預(yù)測模型的預(yù)測精度ait及其預(yù)測值xit即組成了m個二維數(shù)組(a1t,x1t)，(a2t,x2t),…,(amt,xmt).

(13)

1.3.2組合預(yù)測模型計算方法

每個單預(yù)測模型在組合預(yù)測中的有序加權(quán)平均向量為W=(w1,w2,…,wm)T，按由大至小的排序依次列出m個單預(yù)測模型在第t時刻的預(yù)測精度序列a1,a2,…,am，根據(jù)IOWA算子，得到組合預(yù)測模型預(yù)測值的計算公式為

fw{(a1t,x1t),(a2t,x2t),…,(amt,xmt)}=

(14)

由于各單預(yù)測模型的預(yù)測精度隨著預(yù)測時間點的改變而發(fā)生變化，且組合預(yù)測方法的權(quán)系數(shù)僅與單預(yù)測方法的預(yù)測精度相關(guān).因此，組合預(yù)測權(quán)重系數(shù)確定的關(guān)鍵在于單預(yù)測模型在不同時間點上預(yù)測精度的高低[9-10].

設(shè)ea-index(it)=xt-xa-index(it)，N期組合預(yù)測誤差的總平方和S是：

(15)

綜上所述，基于IOWA算子的組合預(yù)測模型的計算公式為

W=minS(w1,w2,…,wm)=

(16)

1.3.3基于IOWA-MC的組合預(yù)測模型

鑒于IOWA算子組合預(yù)測模型中權(quán)系數(shù)與不同預(yù)測時間點的預(yù)測精度存在一定相關(guān)性，且不能預(yù)先獲知未來預(yù)測時間點的預(yù)測精度，因此，對未知序列進(jìn)行預(yù)測時，需要依據(jù)各單預(yù)測模型預(yù)測精度的排序確定單預(yù)測模型的權(quán)系數(shù).馬爾可夫鏈(MC)是一種能夠按照事物原始狀態(tài)和可能狀態(tài)之間的轉(zhuǎn)移概率預(yù)測下一步發(fā)展方向的馬爾可夫過程.因此，將MC與IOWA算子組合預(yù)測模型相結(jié)合，在對將來時間進(jìn)行預(yù)測時，利用MC定性確定每個單預(yù)測模型在預(yù)測時間點上的預(yù)測精度狀態(tài).通過確定預(yù)測精度排序，計算組合預(yù)測模型中的權(quán)系數(shù)，建立組合預(yù)測模型[11-12].設(shè)狀態(tài)轉(zhuǎn)移概率為

(17)

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

(18)

式中：Mi,j(k)為k步轉(zhuǎn)移后狀態(tài)i變?yōu)闋顟B(tài)j的樣本總數(shù)；Mi為處于狀態(tài)i的樣本數(shù)；Pi,j(k)為k步轉(zhuǎn)移后狀態(tài)i變?yōu)闋顟B(tài)j的可能性.構(gòu)建IOWA-MC模型的具體步驟如下.

步驟1計算各單預(yù)測模型的預(yù)測值，并計算各模型的預(yù)測精度.

步驟2按照聚類分析方法，將預(yù)測精度相差不大的預(yù)測時間點視為精度狀態(tài)相同，并對各單預(yù)測模型在預(yù)測期間的精度進(jìn)行狀態(tài)區(qū)間的劃分.根據(jù)各單預(yù)測模型預(yù)測精度范圍，劃分為若干個精度狀態(tài)區(qū)間.

步驟3計算各預(yù)測模型預(yù)測精度的狀態(tài)轉(zhuǎn)移概率矩陣.

步驟4根據(jù)最近時間點所處的精度狀態(tài)和轉(zhuǎn)移概率矩陣可以預(yù)測出待測時間點預(yù)測精度可能的狀態(tài).采用滑動轉(zhuǎn)移概率矩陣方法，得到各預(yù)測模型全部待測時間點的預(yù)測精度狀態(tài).

步驟5根據(jù)IOWA組合預(yù)測模型的建模步驟，確定組合預(yù)測模型的權(quán)系數(shù)，計算組合預(yù)測模型的預(yù)測值.

2 模型應(yīng)用

為了驗證組合預(yù)測模型的合理性，以天津水域2003-2017年水上交通事故的統(tǒng)計數(shù)據(jù)進(jìn)行實例研究.首先，依次利用優(yōu)化的MGM(1,4)預(yù)測模型、SARIMA預(yù)測模型和IOWA算子組合預(yù)測模型對事故數(shù)進(jìn)行預(yù)測；其次，通過對比分析三種預(yù)測模型的預(yù)測誤差，證明本文建立的組合預(yù)測模型相比于單預(yù)測模型來說，能夠作為提高預(yù)測準(zhǔn)確性的有效途徑；最后，對天津水域的事故數(shù)進(jìn)行預(yù)測，得到近三年水域內(nèi)事故發(fā)生趨勢.

2.1 背景值優(yōu)化的多元灰色預(yù)測模型MGM(1,4)

根據(jù)陸夢[13]基于DEMATEL與ISM集成方法的水上交通系統(tǒng)脆性影響因素研究，以及指標(biāo)參數(shù)最好不超過三個的選取原則，最終確定選取日交通流量、船員任職資歷,以及能見度不良天數(shù)三個影響因素與水上交通事故數(shù)建立背景值優(yōu)化的MGM(1,4)預(yù)測模型進(jìn)行實例預(yù)測.原始數(shù)據(jù)序列見表1.

表1 原始事故數(shù)據(jù)序列

根據(jù)式(2)，基于背景值優(yōu)化的MGM(1,4)模型的一次累加生成序列為

根據(jù)式(10)～(11)，利用Matlab編程進(jìn)行計算，得到天津水域2008-2015年水上交通事故的預(yù)測值，并與實際值進(jìn)行曲線擬合，見圖1.

圖1 背景值優(yōu)化的MGM(1,4)預(yù)測模型預(yù)測值與實際值擬合圖

2.2 SARIMA 時間序列預(yù)測模型

根據(jù)時間序列預(yù)測方法的適用條件，選取天津水域2003-2015年實際水上交通事故月度統(tǒng)計量作為實例研究對象.按照1.1.2的建模步驟重新建立預(yù)測模型，通過對模型進(jìn)行識別和檢驗，最終確定SARIMA(1,0,1)(1,1,1)12預(yù)測模型為最佳預(yù)測模型.最后將預(yù)測值和實際值進(jìn)行擬合，見圖2.

圖2 SARIMA(2,0,1)(1,1,1)12模型預(yù)測值與真實值擬合圖

利用該模型對2016-2017年的事故進(jìn)行動態(tài)預(yù)測，將得到的月度預(yù)測值相加，得到年事故量預(yù)測值.

2.3 組合預(yù)測模型及預(yù)測結(jié)果分析

2.3.1組合預(yù)測模型

根據(jù)式(16)，分別計算各單預(yù)測模型得到的2008-2015年預(yù)測事故量及預(yù)測精度，將各單預(yù)測模型的預(yù)測值和預(yù)測精度代入式(16)中，得到基于IOWA算子的最優(yōu)化組合預(yù)測模型為

(19)

利用MATLAB對得到的預(yù)測模型計算式進(jìn)行求解，得到最優(yōu)權(quán)系數(shù)分別為：w1=0.993 6；w2=0.006 4.根據(jù)IOWA算子預(yù)測模型的建模步驟，將求得的w1，w2代入式(19)，由此，即可得到基于IOWA算子的水上交通事故組合預(yù)測模型的預(yù)測值，見表2.將三種預(yù)測模型的預(yù)測值和實際值進(jìn)行擬合，見圖3.

表2 各預(yù)測模型的預(yù)測值

圖3 各預(yù)測模型預(yù)測值與實際值擬合圖

2.3.2組合預(yù)測模型評價

為了進(jìn)一步驗證預(yù)測模型的合理性，分別對所選用的三種預(yù)測模型的預(yù)測誤差進(jìn)行比較，檢驗?zāi)Ｐ偷念A(yù)測效果.根據(jù)組合預(yù)測效果評價原則，分別對各預(yù)測模型的五項誤差指標(biāo)進(jìn)行統(tǒng)計，具體結(jié)果見表3.由表3可知，組合預(yù)測模型的各項誤差指標(biāo)值最小，在整個事故預(yù)測過程中預(yù)測性能最佳.

表3 預(yù)測效果評價

2.3.3未來預(yù)測

根據(jù)背景值優(yōu)化的MGM(1,4)模型，得到2018-2020年水上交通事故的預(yù)測值分別為23.4,23.6,23.9起.根據(jù)SARIMA(1,0,1) (1,1,1)12預(yù)測模型，得到2018-2020年的水上交通事故的預(yù)測值分別為20.6,23.1,18.5起.

將各單一預(yù)測模型2008-2020年的預(yù)測值和預(yù)測精度代入式(19)中，對數(shù)據(jù)進(jìn)行相應(yīng)整理后，得到如下的最優(yōu)化組合預(yù)測模型.

min(w1，w2)=0.942w1+

(20)

利用MATLAB求解上式，得到最優(yōu)權(quán)系數(shù)分別為：w1=0.999 5；w2=0.000 5.將w1，w2代入殘差檢驗公式中，，得到2018，2019及2020年的水上交通事故預(yù)測值分別為21，23，19起.

3 結(jié) 束 語

本文將誘導(dǎo)有序加權(quán)平均算子引入到水上交通事故組合預(yù)測中，所建模型能較好利用不同單一模型信息，突破傳統(tǒng)組合預(yù)測模型的現(xiàn)狀，有效降低預(yù)測誤差，可作為水上交通事故預(yù)測研究的新方法，短期預(yù)測效果顯著，中長期預(yù)測精度仍需深入研究.