王 爽，呂 靖，賴成壽

(大連海事大學(xué)交通運輸工程學(xué)院，遼寧大連116026)

0 引言

原油海運網(wǎng)絡(luò)是原油運輸?shù)闹匾d體，原油海運網(wǎng)絡(luò)中關(guān)鍵節(jié)點的連通狀況會影響原油的安全運輸.而原油海運網(wǎng)絡(luò)中關(guān)鍵節(jié)點的連通性受惡劣天氣、船舶交通事故、海盜及海上恐怖主義、節(jié)點隸屬國家政局不穩(wěn)定及軍事沖突等突發(fā)事件影響，且各突發(fā)事件之間存在相關(guān)關(guān)系，共同影響節(jié)點的連通性.如惡劣天氣很有可能增加發(fā)生船舶交通事故的可能性；在政局不穩(wěn)定及存在軍事沖突的地區(qū)船舶遭受海盜襲擊的可能性會增加.研究各突發(fā)事件相互關(guān)聯(lián)時的原油海運網(wǎng)絡(luò)中關(guān)鍵節(jié)點的連通可靠性，對優(yōu)化進(jìn)口原油海運路徑選擇，保障進(jìn)口原油運輸安全具有重要意義.

連通可靠性的概念由日本學(xué)者M(jìn)ine等于1982年提出，它反映的是交通網(wǎng)絡(luò)兩節(jié)點之間保持連通的概率，一般只研究路段的0和1兩種狀態(tài)，即連通或不連通[1].在Mine的基礎(chǔ)上，Iida等將兩點之間的連通性擴(kuò)展到k點及整個網(wǎng)絡(luò)的可靠性[2].分析連通可靠性的解析方法主要包括最小路/割算法[3]，其計算復(fù)雜度隨網(wǎng)絡(luò)節(jié)點數(shù)增加而指數(shù)增長，Wakayabashi[4]利用布爾代數(shù)，通過Esary-Proschan上下邊界法評估連通可靠性，簡化了計算.

Copula函數(shù)的概念由Sklar于1959年提出，它是將各單變量邊緣分布與聯(lián)合分布連接起來的函數(shù)，使得邊緣分布和聯(lián)合分布函數(shù)的選擇更加靈活[5].Copula函數(shù)被廣泛應(yīng)用于許多領(lǐng)域的相關(guān)關(guān)系建模中[6-7]，但當(dāng)維數(shù)增加時，copula的參數(shù)估計變得困難，因此 Joe[8]，Bedford[9]，Aas[10]等提出并發(fā)展了vine copula.Vine copula通過將高維copula分解成一系列成對二元copula的乘積來簡化高維數(shù)據(jù)聯(lián)合分布的構(gòu)建.Vine copula廣泛應(yīng)用于金融、水文等領(lǐng)域[11-12].

綜上可知，連通可靠性的研究主要集中于城市道路網(wǎng)絡(luò)等陸上交通網(wǎng)絡(luò)，關(guān)于海運網(wǎng)絡(luò)連通可靠性的研究較少，而vine copula函數(shù)能較好地描述變量間的相關(guān)關(guān)系，因此本文基于vine copula函數(shù)構(gòu)建各突發(fā)事件的聯(lián)合分布，從而分析各突發(fā)事件相互關(guān)聯(lián)時的原油海運網(wǎng)絡(luò)中關(guān)鍵節(jié)點的連通可靠性，為優(yōu)化原油海運路徑選擇提供依據(jù).

1 問題描述

根據(jù)國際海事組織全球綜合航運信息系統(tǒng)(GISIS)數(shù)據(jù)庫的統(tǒng)計，影響海運網(wǎng)絡(luò)中節(jié)點連通性的突發(fā)事件主要包括惡劣天氣、交通事故、海盜及海上恐怖主義、節(jié)點隸屬國家政局不穩(wěn)定及軍事沖突.而且這些突發(fā)事件間存在相關(guān)關(guān)系，一種突發(fā)事件的發(fā)生往往會對其他突發(fā)事件的發(fā)生產(chǎn)生影響.例如，惡劣天氣時船舶發(fā)生交通事故的概率很有可能提高；良好的天氣、海況是海盜活動成功的保障因素，因此惡劣天氣時發(fā)生海盜襲擊的可能性會降低；在政局不穩(wěn)定或存在軍事沖突的地區(qū)，遭受海盜襲擊的可能性會提高，如索馬里.這些相互關(guān)聯(lián)的突發(fā)事件共同影響節(jié)點的連通可靠性.

本文的研究目的是要確定各突發(fā)事件相互關(guān)聯(lián)時的原油海運網(wǎng)絡(luò)中關(guān)鍵節(jié)點的連通可靠性，為優(yōu)化原油運輸路徑提供支持.不同突發(fā)事件相互關(guān)聯(lián)時的原油海運網(wǎng)絡(luò)中關(guān)鍵節(jié)點連通可靠性屬于多維聯(lián)合概率問題，令X=(X1,X2,…,Xd)表示各突發(fā)事件，即通過求各突發(fā)事件的聯(lián)合分布F=F(x1,…,xd)=P{X1≤x1,X2≤x2,…,Xd≤xd} ，進(jìn)而計算各關(guān)鍵節(jié)點的連通可靠性.

2 原油海運網(wǎng)絡(luò)中關(guān)鍵節(jié)點連通可靠性度量的vine copula模型構(gòu)建

Vine copula函數(shù)能較好地求解上述存在相關(guān)關(guān)系的多變量聯(lián)合概率問題，通過vine copula函數(shù)，多維聯(lián)合概率密度函數(shù)可以分解為一系列二元copula概率密度函數(shù)的乘積.其中二元copula函數(shù)有多種類型可選擇，可采用擬合優(yōu)度檢驗方法擇優(yōu).基于最優(yōu)的vine copula函數(shù)構(gòu)建的聯(lián)合分布，可以計算各關(guān)鍵節(jié)點的連通可靠性.

2.1 基于高斯混合模型的突發(fā)事件邊緣分布估計

建立vine copula函數(shù)的前提是估計各突發(fā)事件的邊緣分布.高斯混合模型是用有限個高斯密度函數(shù)的線性組合來估計隨機(jī)變量的概率密度分布，如果高斯分布的個數(shù)足夠多，通過調(diào)整線性組合的權(quán)重系數(shù)及每個高斯分布的數(shù)字特征，可以精確描述任何連續(xù)隨機(jī)變量的概率密度分布.因此本文使用高斯混合模型來估計各突發(fā)事件的邊緣分布.

各突發(fā)事件的邊緣分布可以表示為

式中：x是各突發(fā)事件的一維觀測值；M為高斯混合模型的階數(shù)；αl(l=1,2,…,M)是高斯混合模型中子高斯分布的權(quán)重，滿足

Nl(x;μl;σl)(l=1,2,…,M)是高斯分布，可以表示為

式中：μl,σl分別為第l個高斯分布的均值和方差.αl,μl,σl為待估參數(shù)，本文采用EM算法估計高斯混合模型中的參數(shù).

2.2 關(guān)鍵節(jié)點連通可靠性度量的vine copula模型構(gòu)建

Copula是指“連接”多維隨機(jī)變量的聯(lián)合分布和邊緣分布的函數(shù)，它可以將一個聯(lián)合分布分解為它的k個邊際分布和一個copula函數(shù)，其中的copula函數(shù)形式描述了變量間的相關(guān)結(jié)構(gòu)[5].

根據(jù)Sklar定理：令惡劣天氣、交通事故、海盜及海上恐怖主義、節(jié)點隸屬國家政局不穩(wěn)定及軍事沖突等突發(fā)事件的聯(lián)合分布為F=F(x1,…,xd)，邊緣分布分別為F1=F1(x1),…,Fd(xd)，則存在1個Copula函數(shù)C，使得

對于擁有連續(xù)且嚴(yán)格遞增邊緣分布的聯(lián)合分布函數(shù)F，聯(lián)合概率密度函數(shù)f可表示為

Vine copula可以將式(5)分解為d(d-1)/2個二元copula密度函數(shù)[10]，vine copula有兩種比較常用的形式，分別為C-vines和D-vines.C-vines適用于多變量間存在明顯的引導(dǎo)變量的情況，不存在明顯的引導(dǎo)變量時D-vines更適用.本文不存在一個突發(fā)事件引導(dǎo)另一個或兩個突發(fā)事件的情況，因此選擇D-vines copula分解式(5)[9].

式中：fk(xk)是突發(fā)事件的邊緣概率密度函數(shù)；ci,i+j|i+1,…,i+j-1是二元copula的概率密度函數(shù)；βi,i+j|i+1,…,i+j-1為參數(shù)；θ為所有參數(shù)的集合.

在上述vine copula的構(gòu)建過程中包含了邊緣條件分布F(x|v)，Joe[8]指出

式中：vj是向量v的任意一個分量；v-j是向量v除去vj后的分量；Cx,vj|v-j表示二元copula條件分布函數(shù).特別的當(dāng)v是單變量時有

當(dāng)x和v是均勻分布，即f(x)=f(v)=1，F(xiàn)(x)=x，F(xiàn)(v)=v時，h(x,v,Θ)可以用來表示這一條件分布函數(shù)，即有

式中：Θ為參數(shù)集合.

在進(jìn)行D-vine copula結(jié)構(gòu)選擇時，本文通過計算兩兩變量間的秩相關(guān)系數(shù)Kendall’sτ，選擇令每一層的τ值絕對值總和最大的結(jié)構(gòu)[10].式(6)中的二元copula不需要是同一類型，需要選擇能最優(yōu)刻畫每對突發(fā)事件相關(guān)關(guān)系的copula.本文分別選擇Gaussian,Student,Frank,Clayton及Gumbel copula這5種類型的二元copula來刻畫突發(fā)事件間的相關(guān)關(guān)系.并采用逐步半?yún)?shù)法(Stepwise Semiparametric Estimators)估計D-vine copula的參數(shù)[10].本文采用基于經(jīng)驗copula的Cramér-von Mises統(tǒng)計量進(jìn)行擬合優(yōu)度檢驗[13].

3 實證分析

3.1 數(shù)據(jù)選取及處理

本文選擇馬六甲海峽進(jìn)行案例分析.馬六甲海峽是原油海運網(wǎng)絡(luò)中最重要的戰(zhàn)略節(jié)點，其連通可靠性的評估對原油海上運輸安全至關(guān)重要.由于惡劣天氣時海面風(fēng)速較大，因此本文中的惡劣天氣用海面風(fēng)速表示，海面風(fēng)速數(shù)據(jù)來源于Remote Sensing Systems(http://www.remss.com).節(jié)點隸屬國家政局不穩(wěn)定及軍事沖突情況用國家風(fēng)險指數(shù)(International Country Risk Guide)表示，該風(fēng)險指數(shù)由PRS集團(tuán)發(fā)布(https://www.prsgroup.com/).對于馬六甲海峽，其隸屬國家包括新加坡、馬來西亞和印度尼西亞，因此馬六甲海峽的隸屬國家政局不穩(wěn)定及軍事沖突情況用上述3個國家風(fēng)險指數(shù)的平均值表示.海盜及海上恐怖主義用海盜及持械搶劫數(shù)量表示，它和船舶事故數(shù)據(jù)均來源于國際海事組織的全球綜合航運信息系統(tǒng)(GISIS)數(shù)據(jù)庫.

由于海盜和船舶事故并不總是同時發(fā)生，因此本文將突發(fā)事件分成兩種情形.情形1包括3種突發(fā)事件，分別是：海盜及持械搶劫數(shù)量，海盜及持械搶劫發(fā)生時的海面風(fēng)速，以及國家風(fēng)險指數(shù)；情形2也包括3種突發(fā)事件，分別是：船舶事故數(shù)量、船舶事故發(fā)生時的海面風(fēng)速及國家風(fēng)險指數(shù).本文收集了1993—2013年馬六甲海峽發(fā)生的海盜及持械搶劫數(shù)、船舶事故數(shù)，新加坡、馬來西亞和印度尼西亞國家風(fēng)險指數(shù)及海面風(fēng)速.其中海面風(fēng)速指的是海盜及持械搶劫發(fā)生時或船舶事故發(fā)生時海面風(fēng)速的年度平均值.本文分別評價了情形1和情形2中突發(fā)事件相互關(guān)聯(lián)時的馬六甲海峽的連通可靠性.

3.2 結(jié)果及結(jié)果分析

3.2.1 突發(fā)事件邊緣分布

影響馬六甲海峽的各突發(fā)事件邊緣分布的高斯混合模型參數(shù)估計值如表1所示.

表1 突發(fā)事件邊緣分布的參數(shù)估計Table 1 Parameter estimates of marginal distribution functions for extreme events

3.2.2 Vine copula結(jié)構(gòu)選擇及參數(shù)估計

估計各突發(fā)事件邊緣分布函數(shù)之后，計算秩相關(guān)系數(shù)Kendall’sτ的值如表2和表3所示，Kendall’sτ表示各突發(fā)事件間的相關(guān)結(jié)構(gòu).從表2和表3可以看出，國家風(fēng)險指數(shù)和海盜及持械搶劫數(shù)之間有較強(qiáng)的負(fù)相關(guān)關(guān)系，這是因為較高的國家風(fēng)險指數(shù)表示國家的潛在風(fēng)險較低，因此海盜及持械搶劫數(shù)量較少.平均海面風(fēng)速與海盜及持械搶劫數(shù)及船舶事故數(shù)之間也存在相關(guān)關(guān)系.

表2 各突發(fā)事件秩相關(guān)系數(shù)矩陣(情形1)Table 2 Kendall’s Tau dependence matrix of extreme events in Scenario 1

表3 各突發(fā)事件秩相關(guān)系數(shù)矩陣(情形2)Table 3 Kendall’s Tau dependence matrix of extreme events in Scenario 2

根據(jù)表2和表3及D-vine copula的結(jié)構(gòu)選擇準(zhǔn)則，情形1中三維D-vine copula變量的順序為：國家風(fēng)險指數(shù)(x1)，海盜及持械搶劫數(shù)(x2)及海盜及持械搶劫發(fā)生時的海面風(fēng)速(x3).情形2中變量的順序為：國家風(fēng)險指數(shù)(y1)，船舶事故發(fā)生時的海面風(fēng)速(y2)及船舶事故數(shù)(y3).采用逐步半?yún)?shù)法估計D-vine copula的參數(shù)，并根據(jù)AIC值選擇最優(yōu)的二元copula.情形1和2中三維D-vine copula包含的二元copula的類型及其參數(shù)如表4和表5所示.計算基于經(jīng)驗copula的Cramér-von Mises統(tǒng)計量Sn，進(jìn)行擬合優(yōu)度檢驗，如表6所示.

表4 三維D-vine copula的參數(shù)估計(情形1)Table 4 Parameter estimates of the 3-dimensional D-vine copula in Scenario 1

表5 三維D-vine copula的參數(shù)估計(情形2)Table 5 Parameter estimates of the 3-dimensional D-vine copula in Scenario 2

表6 擬合優(yōu)度檢驗Table 6 Goodness-of-fit tests

Sn越小，擬合效果越好，因此由表6可以發(fā)現(xiàn)Sn的檢驗結(jié)果與根據(jù)AIC值選擇的最優(yōu)copula結(jié)果一致.

表4～表6表明，Gaussian和 Clayton copula在刻畫兩兩突發(fā)事件相關(guān)關(guān)系上具有優(yōu)勢.這些二元copula的參數(shù)表明國家風(fēng)險指數(shù)和海盜及持械搶劫數(shù)，以及海盜及持械搶劫數(shù)和海盜及持械搶劫發(fā)生時的海面風(fēng)速間有較強(qiáng)的相關(guān)性.而船舶事故數(shù)與船舶事故發(fā)生時的海面風(fēng)速，以及在海盜及持械搶劫數(shù)條件下，國家風(fēng)險指數(shù)與海盜及持械搶劫發(fā)生時的海面風(fēng)速間的局部相關(guān)性相對較弱.

此外，由于Gaussian copula用來刻畫變量間對稱的尾部相關(guān)性，由此可知國家風(fēng)險指數(shù)與海盜及持械搶劫數(shù)之間具有對稱的尾部相關(guān)性，即無論國家風(fēng)險指數(shù)增加還是減少，只要變化的絕對值相等，他們之間的相關(guān)程度就相同.而Clayton copula對變量在分布下尾部的變化很敏感，因此可以發(fā)現(xiàn)海盜及持械搶劫數(shù)與海盜及持械搶劫發(fā)生時的海面風(fēng)速，船舶事故數(shù)與船舶事故發(fā)生時的海面風(fēng)速間具有下尾相關(guān)性.即當(dāng)海面風(fēng)速降低時，海盜及持械搶劫數(shù)和船舶事故數(shù)減少的可能性會明顯增加，此時他們間的相依關(guān)系很可能會增強(qiáng).通過對各突發(fā)事件相關(guān)關(guān)系的分析，進(jìn)而可以計算各突發(fā)事件相互關(guān)聯(lián)時的節(jié)點連通可靠性.

3.2.3 連通可靠性分析

根據(jù)式(6)可知，由估計出的各二維copula的參數(shù)、邊緣條件分布及邊緣分布求解各突發(fā)事件的聯(lián)合概率密度，進(jìn)而求解聯(lián)合分布、計算聯(lián)合概率是十分困難的.由于海盜及持械搶劫和船舶事故是影響節(jié)點連通可靠性最直接的突發(fā)事件，而且由上面的分析可知海面風(fēng)速及國家風(fēng)險狀況會對海盜及持械搶劫和船舶事故的發(fā)生產(chǎn)生影響，因此本文將給定海面風(fēng)速及國家風(fēng)險指數(shù)時，海盜及持械搶劫數(shù)和船舶事故數(shù)的條件概率作為馬六甲海峽的連通可靠性.

基于上面建立的vine copula，可以獲得條件分布函數(shù)，進(jìn)而可以分析馬六甲海峽的連通可靠性.給定國家風(fēng)險指數(shù)和海面風(fēng)速的值，海盜及持械搶劫數(shù)的條件分布可以表示為F(x2|x1,x3)，根據(jù)式(7)～式(9)，條件分布F(x2|x1,x3)可以表示為

也可以表示為

類似的，給定國家風(fēng)險指數(shù)和海面風(fēng)速值，船舶事故數(shù)的條件分布F(y3|y1,y2)可以表示為

基于vine copula的參數(shù)估計值，可以求出F(x2|x1,x3)和F(y3|y1,y2).本文選擇各突發(fā)事件2013年的取值，計算條件概率，即馬六甲海峽的連通可靠性.給定x1=6.875 0，x3=3.234 4，計算得P(x2≤20|x1=6.875 0,x3=3.234 4)=0.239 5.同樣的，給定y1=6.875 0，y2=10.820 6，計算得P(y3≤ 9|y1=6.875 0,y2=10.820 6)=0.906 1.可以看出在節(jié)點隸屬國家政局不穩(wěn)定及軍事沖突、惡劣天氣和海盜及海上恐怖主義的作用下，馬六甲海峽的連通可靠性較低，僅為0.239 5；而在節(jié)點隸屬國家政局不穩(wěn)定及軍事沖突、惡劣天氣和船舶交通事故的作用下，馬六甲海峽的連通可靠性相對較高，為0.906 1.

4 結(jié)論

本文引入vine copula刻畫各突發(fā)事件間的相關(guān)關(guān)系，建立突發(fā)事件的聯(lián)合分布，并基于各突發(fā)事件的聯(lián)合分布和條件分布進(jìn)行節(jié)點連通可靠性分析，本文以馬六甲海峽為例進(jìn)行案例分析.結(jié)果表明，節(jié)點隸屬國家政局不穩(wěn)定及軍事沖突與海盜及海上恐怖主義間具有對稱的尾部相關(guān)性，而海盜及海上恐怖主義和船舶交通事故分別在惡劣天氣轉(zhuǎn)好時與其具有更強(qiáng)的相關(guān)性，即具有下尾相關(guān)性.此外，在節(jié)點隸屬國家政局不穩(wěn)定及軍事沖突、惡劣天氣、海盜及海上恐怖主義的作用下馬六甲海峽的連通可靠性較低，僅為0.239 5；而在節(jié)點隸屬國家政局不穩(wěn)定及軍事沖突、惡劣天氣、船舶交通事故的作用下馬六甲海峽的連通可靠性相對較高，研究結(jié)果可以為制定原油海上運輸策略提供依據(jù).