許循齊+陳悠超+茅晨昊+盧憶寧

摘 要：針對(duì)橋區(qū)水域船舶通航風(fēng)險(xiǎn)評(píng)估中，蒙特卡洛仿真收斂速度慢的缺點(diǎn)，在風(fēng)險(xiǎn)概率成對(duì)數(shù)正態(tài)分布的前提下，采用馬爾科夫鏈建立風(fēng)險(xiǎn)評(píng)估模型。將橋區(qū)水域船舶通航風(fēng)險(xiǎn)劃分為三種狀態(tài)，利用歷史數(shù)據(jù)以及蒙特卡洛仿真數(shù)據(jù)統(tǒng)計(jì)獲得狀態(tài)轉(zhuǎn)移矩陣，結(jié)合當(dāng)前船舶航行狀態(tài)，快速求解未來(lái)一段時(shí)間內(nèi)的狀態(tài)轉(zhuǎn)移概率，實(shí)現(xiàn)了風(fēng)險(xiǎn)的快速評(píng)估。算例結(jié)果表明：三種狀態(tài)的平穩(wěn)狀態(tài)概率為39.07%，59.43%，1.5%，轉(zhuǎn)移概率求解誤差率為1.24%，1.06%，16.84%，求解時(shí)間僅為0.053秒。船舶駕駛員可利用該模型進(jìn)行實(shí)時(shí)評(píng)估，以降低橋區(qū)水域通航安全風(fēng)險(xiǎn)。

關(guān)鍵詞：橋區(qū)水域；蒙特卡洛；馬爾科夫；快速評(píng)估

中圖分類(lèi)號(hào)：U692.3 文獻(xiàn)標(biāo)志碼：A 文章編號(hào)：2095-2945（2017）33-0011-02

引言

隨著航運(yùn)經(jīng)濟(jì)的不斷發(fā)展，船舶呈現(xiàn)大型化、高速化發(fā)展趨勢(shì)，與此同時(shí)，橋梁建造數(shù)量與日俱增，致使橋區(qū)水域船撞橋事故頻發(fā)。從上世紀(jì)60年代至今，國(guó)內(nèi)船撞橋事故近400起，其中重大事故近40起，每年平均1起重大事故[1]。由此可見(jiàn)，船撞橋事故不僅會(huì)造成船舶以及橋梁自身的損傷，人員傷亡，還會(huì)對(duì)橋區(qū)水域環(huán)境造成污染。

對(duì)于橋區(qū)水域船撞橋研究主要分為三類(lèi)。第一類(lèi)為船橋碰撞概率研究，有歐洲規(guī)范模型、昆茲模型等， 國(guó)內(nèi)主要為戴彤宇的簡(jiǎn)化模型[2]。第二類(lèi)為事故研究，Mastiglio統(tǒng)計(jì)了1960-2014年國(guó)際上由于船撞導(dǎo)致的橋梁毀壞事故[3]。第三類(lèi)為橋區(qū)水域主被動(dòng)防撞研究，陳明棟等[4]提出“攔-防”組合的非通航孔橋墩防撞方法。

顯然，對(duì)橋區(qū)水域風(fēng)險(xiǎn)評(píng)估方面研究?jī)H限于某一時(shí)間點(diǎn)的狀態(tài)研究，并不能滿(mǎn)足船舶穿越橋區(qū)通航孔整個(gè)過(guò)程中對(duì)時(shí)間的要求。

然而，掌握當(dāng)前船舶航行狀態(tài)風(fēng)險(xiǎn)是基礎(chǔ)，預(yù)測(cè)未來(lái)船舶航行狀態(tài)風(fēng)險(xiǎn)才是重中之重。

本文對(duì)橋區(qū)水域通航風(fēng)險(xiǎn)進(jìn)行量化，以歷史險(xiǎn)情與事故樣本為基礎(chǔ)，通過(guò)蒙特卡洛方法進(jìn)行仿真處理，統(tǒng)計(jì)得到通航狀態(tài)之間的相互轉(zhuǎn)移矩陣，利用橋區(qū)船舶當(dāng)前的狀態(tài)初始矩陣及馬爾科夫狀態(tài)轉(zhuǎn)移矩陣，快速判別出船舶通航風(fēng)險(xiǎn)概率。

1 馬爾科夫方法概述

設(shè){X（n），n=0，1，2，...}是一馬爾科夫鏈，狀態(tài)空間為E={0，1，2，...}，若其一步轉(zhuǎn)移概率與馬爾科夫鏈現(xiàn)在所在時(shí)刻無(wú)關(guān)，即滿(mǎn)足等式（1）的要求。

此時(shí)馬氏鏈從i狀態(tài)轉(zhuǎn)移到j(luò)狀態(tài)的概率與現(xiàn)在所在時(shí)刻k無(wú)關(guān)，只與現(xiàn)在所處狀態(tài)i有關(guān)，則稱(chēng)此為齊次馬爾科夫鏈，亦稱(chēng)之具有平穩(wěn)轉(zhuǎn)移概率的馬爾科夫鏈。

常記一步轉(zhuǎn)移概率為Pij=Pij（1）。

齊次馬爾科夫鏈的一步轉(zhuǎn)移概率Pij具有如下性質(zhì)：

若將全部的一步轉(zhuǎn)移概率表示為矩陣的形式，則有

對(duì)于相互獨(dú)立且同分布隨機(jī)變量列，其一步轉(zhuǎn)移概率與時(shí)刻無(wú)關(guān)，且與現(xiàn)在狀態(tài)i也無(wú)關(guān)，即

它的一步轉(zhuǎn)移概率矩陣為

（5）

顯然，當(dāng)歷史樣本數(shù)據(jù)足夠多時(shí)，可根據(jù)系統(tǒng)狀態(tài)的變化，統(tǒng)計(jì)獲得誤差足夠小的馬爾科夫狀態(tài)轉(zhuǎn)移概率矩陣P。通常采用頻率表示概率p，即樣本點(diǎn)在單位時(shí)間△t內(nèi)事件發(fā)生率來(lái)表示事故概率。當(dāng)歷史樣本數(shù)量趨近于無(wú)窮時(shí)，轉(zhuǎn)移概率矩陣P的統(tǒng)計(jì)值將無(wú)限趨近于其真實(shí)值[5]。

2 橋區(qū)水域風(fēng)險(xiǎn)劃分

基于船橋碰撞概率，文獻(xiàn)調(diào)研以及專(zhuān)家評(píng)測(cè)，本文將橋區(qū)水域風(fēng)險(xiǎn)分為三種狀態(tài)。狀態(tài)一為安全狀態(tài)，狀態(tài)二為臨界狀態(tài)，狀態(tài)三為危險(xiǎn)狀態(tài)。安全狀態(tài)是指船舶根據(jù)既定航向航速可順利穿越通航孔。臨界狀態(tài)是指船橋之間已構(gòu)成碰撞危險(xiǎn)，船舶只有通過(guò)改速改向才可導(dǎo)致船橋在安全距離上駛過(guò)。危險(xiǎn)狀態(tài)是指單憑船舶自身的避讓行動(dòng)已無(wú)法避免碰撞的局面。

3 橋區(qū)水域風(fēng)險(xiǎn)模型構(gòu)建

橋區(qū)水域風(fēng)險(xiǎn)的三種狀態(tài)空間E={0，1，2}，其狀態(tài)轉(zhuǎn)移矩陣如式（6）所示。

假設(shè)橋區(qū)水域初始狀態(tài)為？仔（0），經(jīng)過(guò)△t時(shí)間后轉(zhuǎn)移到下一個(gè)狀態(tài)？仔（1），則？仔（1）=？仔（0）P，再經(jīng)過(guò)△t時(shí)間后，轉(zhuǎn)移到下一個(gè)狀態(tài)？仔（2），則？仔（2）=？仔（1）P=？仔（0）P2，依次類(lèi)推，即可得到橋區(qū)水域風(fēng)險(xiǎn)狀態(tài)轉(zhuǎn)移關(guān)系方程？仔（n）=？仔（n-1）P=...=？仔（1）Pn-1=？仔（0）Pn，由此，在已知初始狀態(tài)的情況下即可快速分析出間隔任意時(shí)刻的風(fēng)險(xiǎn)狀態(tài)水平。

對(duì)于橋區(qū)風(fēng)險(xiǎn)馬爾科夫鏈模型，當(dāng) 時(shí)， 最終趨近于一個(gè)穩(wěn)定值，滿(mǎn)足式（7）

本文將狀態(tài)空間E={0，1，2}分為兩類(lèi)，一類(lèi)為可接受狀態(tài)A={1，2}，另一類(lèi)為不可接受狀態(tài)U={3}。風(fēng)險(xiǎn)狀態(tài)轉(zhuǎn)移矩陣也相應(yīng)地劃分為四類(lèi)。第一類(lèi)二階轉(zhuǎn)移矩陣為可接受狀態(tài)之間的轉(zhuǎn)移概率，第二類(lèi)1×2的轉(zhuǎn)移矩陣為不可接受狀態(tài)向可接受狀態(tài)轉(zhuǎn)移的概率，第三類(lèi)2×1的轉(zhuǎn)移矩陣為可接受狀態(tài)向不可接受狀態(tài)轉(zhuǎn)移的概率，第四類(lèi)一階矩陣為不可接受狀態(tài)之間的轉(zhuǎn)移概率。

（8）

4 算例分析

由于船橋碰撞事故為小樣本事件，樣本特征的波動(dòng)性強(qiáng)，因此，對(duì)樣本數(shù)據(jù)進(jìn)行概率分析后，引入蒙特卡洛方法對(duì)船橋碰撞風(fēng)險(xiǎn)進(jìn)行仿真，以獲得大樣本信息，對(duì)整體態(tài)勢(shì)進(jìn)行評(píng)價(jià)。根據(jù)海事局事故及險(xiǎn)情報(bào)告[6]，通過(guò)上述方法進(jìn)行事故狀態(tài)可能性分析[7]?！鱰時(shí)間間隔內(nèi)，發(fā)生事故概率0.058%，采取措施可避免事故的概率為1.521%。本文運(yùn)用蒙特卡洛方法進(jìn)行仿真10000次，每次仿真200個(gè)時(shí)間間隔[8]。

通過(guò)表1中的數(shù)據(jù)，可以獲得轉(zhuǎn)移矩陣P。

（9）

由轉(zhuǎn)移矩陣可知，經(jīng)過(guò)時(shí)間△t，船橋碰撞風(fēng)險(xiǎn)保持安全狀態(tài)的概率為0.986699，從安全狀態(tài)轉(zhuǎn)變?yōu)榕R界狀態(tài)的概率為0.013267，從安全狀態(tài)轉(zhuǎn)變?yōu)槲ｋU(xiǎn)狀態(tài)的概率為0.000034。假設(shè)，初始狀態(tài)為安全狀態(tài)，表示為？仔（0）=[1 0 0]，經(jīng)過(guò)200次的矩陣計(jì)算之后，得到？仔（0）～？仔（200）的值，運(yùn)行時(shí)間為0.053秒，幾乎可以忽略不計(jì)。統(tǒng)計(jì)每次仿真，三種狀態(tài)的次數(shù)，即？仔'（0）～？仔'（200）。將三種狀態(tài)下？仔（0）～？仔（200）與？仔'（0）～？仔'（200）相比較，如圖1所示。

根據(jù)式（10）可計(jì)算出200個(gè)時(shí)間間隔內(nèi)演算值與統(tǒng)計(jì)值之間的誤差。安全狀態(tài)誤差為1.24%，臨界狀態(tài)誤差為1.06%，危險(xiǎn)狀態(tài)誤差為16.84%。由于危險(xiǎn)狀態(tài)的樣本量相對(duì)較小，即？仔（j）較小，從而導(dǎo)致危險(xiǎn)狀態(tài)的誤差較大。

利用式（7）可以求出橋區(qū)水域風(fēng)險(xiǎn)平穩(wěn)狀態(tài)概率？仔（∞）=[0.3900 0.5943 0.015]，200個(gè)時(shí)間間隔后，狀態(tài)概率矩陣開(kāi)始收斂于？仔（∞）。

5 結(jié)束語(yǔ)

借助MCMC算法，對(duì)船橋碰撞風(fēng)險(xiǎn)進(jìn)行了劃分，相比于傳統(tǒng)的評(píng)估方法，更加簡(jiǎn)便快捷。用蒙特卡洛方法產(chǎn)生的大樣本情況下，統(tǒng)計(jì)獲得狀態(tài)轉(zhuǎn)移矩陣P，將其與船舶航行初始狀態(tài)相結(jié)合，可以快速求得未來(lái)一段時(shí)間內(nèi)船橋碰撞的概率變化，并可求解出橋區(qū)水域船舶通航風(fēng)險(xiǎn)的平穩(wěn)狀態(tài)概率。因此，本文研究具有理論意義，而且對(duì)降低橋區(qū)水域船舶通航安全風(fēng)險(xiǎn)有著重要的顯示價(jià)值。

參考文獻(xiàn)：

[1]黃常海，胡甚平，高德毅，等.大橋船撞動(dòng)態(tài)風(fēng)險(xiǎn)評(píng)估系統(tǒng)的設(shè)計(jì)與實(shí)現(xiàn)[J].中國(guó)安全科學(xué)學(xué)報(bào)，2013，23（4）：120-126.

[2]Kunz CU. Ship bridge collision in river traffic， analysis and design practice [C]. Ship Collision Analysis，1998：13-21.

[3]Mastaglio L. Bridge Bashing[J].Civil Engineering-ASCE，2014，67：38-40.

[4]陳明棟，向 .非通航孔橋墩“攔-防”組合防撞方案應(yīng)用研究[J].水利水運(yùn)工程學(xué)報(bào)，2014，8（04）：98-104.

[5]李裕奇.隨機(jī)過(guò)程[M].北京：國(guó)防工業(yè)出版社，2012：219-230.

[6]上海市航海學(xué)會(huì).東海大橋船撞風(fēng)險(xiǎn)評(píng)估與橋區(qū)安全管理對(duì)策專(zhuān)項(xiàng)論證報(bào)告[R].上海：洋山港海事局，2016.

[7]胡甚平.海上交通系統(tǒng)風(fēng)險(xiǎn)蒙特卡洛仿真[J].上海海事大學(xué)學(xué)報(bào)，2011，32（4）：7-16.

[8]徐鐘濟(jì).蒙特卡羅方法[M].上海：上?？茖W(xué)技術(shù)出版社，1985：8.endprint