基于蒙特卡羅方法的廣州港航行風(fēng)險(xiǎn)評(píng)估

張靜 張志堅(jiān) 周芯玉

(1.廣州市氣象臺(tái)，廣東 廣州 511430； 2.廣州市突發(fā)事件預(yù)警信息發(fā)布中心，廣東 廣州 511430)

引言

廣州港地處珠江口中心地帶，是經(jīng)濟(jì)高度集中的區(qū)域，也是氣象災(zāi)害脆弱區(qū)和高暴露區(qū)，南海地區(qū)天氣復(fù)雜多變，災(zāi)害性天氣時(shí)常發(fā)生，夏季強(qiáng)對(duì)流頻發(fā)[1-2]，冬春季則容易發(fā)生海霧現(xiàn)象。廣州港幾乎每年都遭受強(qiáng)對(duì)流、海霧等惡劣天氣襲擊，給廣州港的運(yùn)行造成不利影響和經(jīng)濟(jì)損失[3]。海上的風(fēng)力等級(jí)、浪的高度、能見度距離都對(duì)海上航行有重大威脅。前人對(duì)海上航行氣象風(fēng)險(xiǎn)評(píng)估的主要做法是利用風(fēng)浪進(jìn)行的[4-7]。劉振[4]分析了黃渤海區(qū)因大風(fēng)浪導(dǎo)致的航行事故，構(gòu)造了黃渤海區(qū)船舶在大風(fēng)浪條件下的風(fēng)險(xiǎn)估算模型，該模型實(shí)現(xiàn)了在遇到大風(fēng)浪時(shí)對(duì)航線上重點(diǎn)船舶提前預(yù)警的功能，對(duì)減少海上事故和降低風(fēng)險(xiǎn)有顯著意義。李耿[5]統(tǒng)計(jì)分析了小型登陸艇受大風(fēng)浪影響的情況，搭建了大風(fēng)浪對(duì)小型登陸艇的風(fēng)險(xiǎn)評(píng)價(jià)系統(tǒng)，研究表明大風(fēng)浪下的風(fēng)險(xiǎn)預(yù)警方案對(duì)保障小型登陸艇航行有很大價(jià)值。前人大部分研究沒有充分考慮能見度對(duì)航行風(fēng)險(xiǎn)的影響。而大霧和強(qiáng)對(duì)流造成的低能見度事件對(duì)廣州港船舶航行和港口作業(yè)有重大影響，容易引發(fā)擱淺、碰撞或觸礁等類似險(xiǎn)情，給海上航行造成嚴(yán)重隱患[8-11]。為了廣州港安全作業(yè)，廣州海事局按照廣州港水域的狀況，制定了能見度低時(shí)的相關(guān)航行規(guī)則，在保障廣州港安全作業(yè)上發(fā)揮著積極的作用。但廣州港作為珠三角的航運(yùn)樞紐，船舶流量增長(zhǎng)迅速，而南方天氣復(fù)雜多變，天氣條件導(dǎo)致能見度變化頻繁，有關(guān)能見度低時(shí)的航行規(guī)則使得廣州港得不到最大化的利用[11]。為合理調(diào)配廣州港的資源，在充分考慮能見度對(duì)船舶航行和廣州港安全生產(chǎn)影響的同時(shí)，有必要對(duì)廣州港在不同等級(jí)能見度下的航行風(fēng)險(xiǎn)展開研究，將定性的廣州港航行風(fēng)險(xiǎn)進(jìn)行定量化，從而可以按照不同能見度條件采取更科學(xué)的風(fēng)險(xiǎn)控制措施。

近年來蒙特卡羅方法在小樣本事件分析中逐漸被廣泛運(yùn)用，趙濱等[12]在對(duì)T639和Grapes模式24降水預(yù)報(bào)效果評(píng)估中，采用蒙特卡羅方法重構(gòu)樣本，有效解決了評(píng)估差異是由于模式預(yù)報(bào)效果差異還是樣本數(shù)據(jù)的不確定性所造成的困擾。本文以提高廣州港水上交通氣象保障服務(wù)能力為核心，采用蒙特卡羅仿真對(duì)不同能見度下的航行氣象風(fēng)險(xiǎn)展開分析，為廣州港生產(chǎn)調(diào)度和船舶安全航行提供氣象科技支撐。本文統(tǒng)計(jì)分析了廣州港水域出現(xiàn)的航行事故信息，將船舶出現(xiàn)事故的頻率以及造成后果的嚴(yán)重性定量化，接著對(duì)二者建立概率分布模型，然后分別對(duì)不同等級(jí)能見度下的事故數(shù)據(jù)進(jìn)行蒙特卡羅仿真，綜合分析后將航行風(fēng)險(xiǎn)定量化。最終對(duì)能見度在0—12 km內(nèi)的風(fēng)險(xiǎn)值進(jìn)行擬合，得到能見度在0—12 km內(nèi)的風(fēng)險(xiǎn)分布特征，從而對(duì)廣州港水域航行風(fēng)險(xiǎn)作出合理的評(píng)估，為廣州港管理部門掌握不同能見度時(shí)的風(fēng)險(xiǎn)狀況提供參考，對(duì)水上交通，尤其港口(近海)的航行氣象保障技術(shù)是有力的完善和補(bǔ)充，是提升廣州港這一“海上絲綢之路”氣象保障能力的有效手段。

1 資料與方法

1.1 資料來源

使用的廣州港事故數(shù)據(jù)來源于廣州海事局，資料包括發(fā)生時(shí)間、遇險(xiǎn)狀況和救助結(jié)果；能見度數(shù)據(jù)來自于廣東省氣象局觀測(cè)站網(wǎng)，共有30個(gè)能見度自動(dòng)觀測(cè)站，儀器型號(hào)為美國(guó)Belfort公司生產(chǎn)的Model 6000能見度自動(dòng)觀測(cè)儀，數(shù)據(jù)時(shí)間分辨率為1 min，分析時(shí)能見度取10 min滑動(dòng)平均，資料起始時(shí)間為2014年6月5日至2019年12月31日。

1.2 航行風(fēng)險(xiǎn)建模

海上航行風(fēng)險(xiǎn)涉及兩個(gè)要素：一個(gè)是出現(xiàn)事故的可能性；另一個(gè)是事故導(dǎo)致?lián)p害的嚴(yán)重性。因此，風(fēng)險(xiǎn)可以看成是這二者的組合。由于事故樣本量比較少，引入蒙特卡羅方法可以將有限的事故數(shù)據(jù)在合理?xiàng)l件下進(jìn)行放大。蒙特卡羅(Monte Carlo)仿真是以概率模型為基本思想的統(tǒng)計(jì)實(shí)驗(yàn)方法[12-18]，在港口航行風(fēng)險(xiǎn)分析中，該方法被用作樣本數(shù)據(jù)的模擬仿真，簡(jiǎn)稱MC方法[13]。先構(gòu)造合適的概率模型來明確船舶航行風(fēng)險(xiǎn)事件，將船舶航行事故出現(xiàn)的可能性與造成后果的嚴(yán)重性分別轉(zhuǎn)化為具體的數(shù)學(xué)問題，使得航行風(fēng)險(xiǎn)值接近于風(fēng)險(xiǎn)事件出現(xiàn)的概率[12-14]。接著借助MC方法開展仿真實(shí)驗(yàn)，多次仿真得到更多的隨機(jī)樣本，然后用統(tǒng)計(jì)學(xué)方法分析仿真結(jié)果，將風(fēng)險(xiǎn)事件發(fā)生的概率近似為風(fēng)險(xiǎn)結(jié)果定量的數(shù)值解。

廣州港航行風(fēng)險(xiǎn)是指在人、船、環(huán)境、管理組成的船舶交通系統(tǒng)中，在某種狀態(tài)下，導(dǎo)致系統(tǒng)中的各部分遭到損害的可能性和這種損害嚴(yán)重程度的組合[19-22]?？赡苄灾甘鹿食霈F(xiàn)的機(jī)會(huì)，用頻率來表示，損害嚴(yán)重程度用后果來表示。頻率p用式(1)表示，是指每年出現(xiàn)的事故次數(shù)n占總的船舶活動(dòng)量N的比例。后果涉及的內(nèi)容一般分為：人身傷亡、直接經(jīng)濟(jì)虧損、環(huán)境污染損害和社會(huì)效應(yīng)等[20-22]，可以簡(jiǎn)單地概括為人命損失、財(cái)產(chǎn)損失。根據(jù)廣州港航行事故信息，按照遇險(xiǎn)情況和救助結(jié)果，考慮采取“事故等效后果”來表述[19-21]，例如一箱貨物落水與直接經(jīng)濟(jì)損失10萬元的后果等效[20]，按照損失程度可用“1—10”來衡量，從而將后果定量化?！笆鹿实刃Ш蠊眂用式(2)表示，是指單位時(shí)間內(nèi)事故造成的等效后果總和K與出現(xiàn)次數(shù)n的比值[19-21]。

p=n/N

(1)

式(1)中，p為事故出現(xiàn)的頻率；n為每年出現(xiàn)的事故數(shù)；N為每年的船舶活動(dòng)量。

c=K/n

(2)

式(2)中，c為事故造成的后果；n為每年出現(xiàn)的事故數(shù)；K為每年事故造成的各種后果總和。

國(guó)際海事組織(IMO)在舉辦的海上安全委員會(huì)(MSC)第74屆會(huì)議上通過的《FSA應(yīng)用正式指南》[21]將風(fēng)險(xiǎn)定義為具有頻率和后果程度的雙重特性，其函數(shù)關(guān)系式見式(3)：

R=f(p,c)

(3)

從式(3)可以看出，風(fēng)險(xiǎn)可以看成系統(tǒng)內(nèi)事故出現(xiàn)的頻率；風(fēng)險(xiǎn)也可以看成一次事故導(dǎo)致?lián)p害的后果程度[21]，所以式(3)可表示為：

(4)

1.3 概率模型構(gòu)建

在MC中確定合適的概率分布模型是很重要的工作。在現(xiàn)有的樣本信息中，確定各種隨機(jī)樣本的性質(zhì)，選用合適的概率分布模型對(duì)廣州港風(fēng)險(xiǎn)事件進(jìn)行仿真，充分增加廣州港船舶風(fēng)險(xiǎn)的數(shù)據(jù)樣本，進(jìn)而對(duì)風(fēng)險(xiǎn)定量化展開分析。

1.3.1 船舶發(fā)生事故的可能性

當(dāng)不知道船舶發(fā)生事故的可能性是多少時(shí)，可以認(rèn)為事故率是[0，1]區(qū)間里的任意數(shù)字。而貝塔(Beta)分布常用來表示一個(gè)事件出現(xiàn)的概率的概率分布。Beta分布是一組定義在[0，1]區(qū)間里的連續(xù)概率分布，具有兩個(gè)參數(shù)α、β，且均為正值。Beta分布的概率密度函數(shù)為：

(5)

1.3.2 事故的后果

各事故造成的后果之間相互獨(dú)立，且都為正值，有可加性，所以可認(rèn)為事故后果服從伽馬(Gamma)分布。Gamma分布的概率密度函數(shù)為：

(6)

式(6)中，α為形狀參數(shù)；β為尺度參數(shù)，均為正值。

根據(jù)2014—2019年的樣本信息可以得到能見度不良(Vis≤5 km)情況下平均發(fā)生的事故率為2.287×10-5，后果為8.653，能見度受限(5 km

表1 Beta分布和Gamma分布的參數(shù)Table 1 Parameters of Bela distribution and Gamma distribution

各參數(shù)值見表1。實(shí)現(xiàn)MC仿真的基本工具就是根據(jù)已知樣本產(chǎn)生隨機(jī)數(shù)。風(fēng)險(xiǎn)是船舶事故可能性和后果程度的組合。二者的任意變量之間相互獨(dú)立，且對(duì)總體風(fēng)險(xiǎn)影響都較小，因此,可分別對(duì)事故可能性和后果進(jìn)行隨機(jī)抽樣。按照三種能見度條件下樣本數(shù)據(jù)的均值μ，標(biāo)準(zhǔn)差σ，分別生成服從Beta分布模型的事故率的隨機(jī)數(shù)和服從Gamma分布模型的后果隨機(jī)數(shù)。

1.3.3 模型的檢驗(yàn)

模型仿真數(shù)據(jù)的檢驗(yàn)方法主要采用文獻(xiàn)[13-15]中的變異系數(shù)。根據(jù)大數(shù)法則和中心極限定理，在抽取一定的樣本后，需要判斷估計(jì)值與待求真值的收斂性檢驗(yàn)仿真數(shù)據(jù)的可靠性。一般地，引入變異系數(shù)作為計(jì)算收斂判斷依據(jù)。

(7)

式(7)中，σ為風(fēng)險(xiǎn)標(biāo)準(zhǔn)差；μ為風(fēng)險(xiǎn)平均值。

由表2可以看出，原始數(shù)據(jù)和MC仿真模擬的數(shù)據(jù)均值、標(biāo)準(zhǔn)差及變異系數(shù)都相差不大，可見仿真數(shù)據(jù)有效，模型參數(shù)合理。

表2 風(fēng)險(xiǎn)統(tǒng)計(jì)特征值Table 2 Statistical characteristic value of risk

1.4 航行風(fēng)險(xiǎn)MC仿真流程

(1)確定隨機(jī)變量。確定對(duì)航行風(fēng)險(xiǎn)有影響的因素，把航行風(fēng)險(xiǎn)構(gòu)造成一個(gè)概率模型，模型主要涉及兩個(gè)隨機(jī)變量，即事故出現(xiàn)的頻率及其造成的后果。

(2)數(shù)據(jù)收集及統(tǒng)計(jì)。收集并分析與航行風(fēng)險(xiǎn)相關(guān)的樣本信息，得到事故率和事故后果的樣本數(shù)據(jù)。

(3)分析上述兩個(gè)變量的性質(zhì)，確定二者的概率分布模型。

(4)為這兩個(gè)變量分別抽取隨機(jī)數(shù)。

(5)根據(jù)隨機(jī)抽樣的數(shù)據(jù)計(jì)算出航行風(fēng)險(xiǎn)值。

(6)重復(fù)(4)和(5)兩個(gè)步驟，直到滿足需要的模擬次數(shù)。

(7)對(duì)仿真實(shí)驗(yàn)數(shù)據(jù)作統(tǒng)計(jì)處理，分析得到分布曲線，并檢驗(yàn)其概率分布，計(jì)算平均值和標(biāo)準(zhǔn)差及其他統(tǒng)計(jì)特征；最終以平均值、標(biāo)準(zhǔn)差等其他統(tǒng)計(jì)特征作為航行風(fēng)險(xiǎn)評(píng)估的結(jié)果。

具體仿真流程如圖1。

圖1 廣州港航行風(fēng)險(xiǎn)仿真流程Fig.1 Simulation process of navigation risk of Guangzhou port

2 結(jié)果分析

根據(jù)廣州港2014—2019年的航行事故信息統(tǒng)計(jì)，按照三種等級(jí)能見度距離(能見度不良，能見度受限，能見度良好)對(duì)其分類。由于廣州港事故數(shù)據(jù)樣本量較少，屬于小樣本事件，出現(xiàn)風(fēng)險(xiǎn)事件的可能性和后果都屬于有限個(gè)離散數(shù)據(jù)，風(fēng)險(xiǎn)在不同特征條件下呈現(xiàn)較大的波動(dòng)性。為獲取航行風(fēng)險(xiǎn)事件發(fā)生的可能性只能用頻率進(jìn)行估計(jì)，后果采取“事故等效后果”進(jìn)行換算。分別對(duì)事故可能性建立Beta分布模型，對(duì)事故后果建立Gamma分布模型，進(jìn)行蒙特卡羅仿真，每次仿真抽取500個(gè)隨機(jī)數(shù)，求均值得出一組數(shù)據(jù)，連續(xù)仿真得出2000組數(shù)據(jù)。最終對(duì)不同等級(jí)能見度下的風(fēng)險(xiǎn)值進(jìn)行擬合，得出廣州港在不同能見度下的總體風(fēng)險(xiǎn)分布。

2.1 不同能見度下的風(fēng)險(xiǎn)仿真

廣州港能見度低絕大多數(shù)是天氣狀況造成的，通常是由于海霧、連陰雨、強(qiáng)對(duì)流、霾等導(dǎo)致能見度受到限制。能見度不佳時(shí)，視線受到阻礙，對(duì)周圍物標(biāo)、航標(biāo)和船舶動(dòng)態(tài)辨別不清，給船舶的避讓、導(dǎo)航及定位帶來了較大的不便，航行風(fēng)險(xiǎn)很大。根據(jù)廣州港水域的情況，將廣州港能見度分成3個(gè)等級(jí)，即能見度不良(Vis≤5 km)、能見度受限(5 km

由圖2可知，能見度不良(Vis≤5 km)時(shí)，廣州港航行平均事故發(fā)生率為2.753×10-5，平均后果嚴(yán)重程度為7.046，平均風(fēng)險(xiǎn)值為1.92×10-4，風(fēng)險(xiǎn)值標(biāo)準(zhǔn)差7.197×10-5。能見度不良時(shí)廣州港出現(xiàn)險(xiǎn)情的幾率最小，而造成的后果最為嚴(yán)重，整體的風(fēng)險(xiǎn)值最低。樣本數(shù)據(jù)波動(dòng)性較小。

圖2 能見度不良時(shí)事故后果值頻率(a)、風(fēng)險(xiǎn)值占比(b)、仿真次數(shù)平均風(fēng)險(xiǎn)值(c)和平均風(fēng)險(xiǎn)值占比(d)仿真圖Fig.2 Simulation result of frequency of accident consequence (a),percentage of risk value (b),average risk value of the number of simulation times (c),and percentage of average risk value (d) for poor visibility

由圖3可知，能見度受限(5 km

圖3 能見度受限時(shí)時(shí)事故后果值頻率(a)、風(fēng)險(xiǎn)值占比(b)、仿真次數(shù)平均風(fēng)險(xiǎn)值(c)和平均風(fēng)險(xiǎn)值占比(d)仿真圖Fig.3 Simulation result of frequency of accident consequence (a),percentage of risk value (b),average risk value of the number of simulation times (c),and percentage of average risk value (d) for limited visibility

由圖4可知，能見度良好(Vis≥10 km)的情況下，廣州港航行平均事故發(fā)生率為5.830×10-5，平均后果嚴(yán)重程度為3.615，平均風(fēng)險(xiǎn)值為2.647×10-4，風(fēng)險(xiǎn)值標(biāo)準(zhǔn)差為1.278×10-4。結(jié)果表明，當(dāng)能見度良好時(shí)，廣州港出現(xiàn)險(xiǎn)情的幾率相對(duì)較小，且造成后果的嚴(yán)重性最小，因而整體的風(fēng)險(xiǎn)值不大。說明廣州港船舶航行受能見度的影響很大。樣本數(shù)據(jù)總體波動(dòng)不大。

圖4 能見度良好時(shí)事故后果值頻率(a)、風(fēng)險(xiǎn)值占比(b)、仿真次數(shù)平均風(fēng)險(xiǎn)值(c)和平均風(fēng)險(xiǎn)值占比(d)仿真圖Fig.4 Simulation result of frequency of accident consequence (a),percentage of risk value (b),average risk value of the number of simulation times (c),and percentage of average risk value (d) for good visibility

2.2 風(fēng)險(xiǎn)曲線擬合

采用蒙特卡羅方法對(duì)廣州港航行事故樣本的仿真，得出在三種等級(jí)能見度下船舶航行的風(fēng)險(xiǎn)值。在實(shí)際航行中，12 km以下的能見度可對(duì)船舶航行造成一定的影響，因此，為了得到能見度在0—12 km連續(xù)狀態(tài)下的航行風(fēng)險(xiǎn)分布曲線，在0—12 km范圍內(nèi)逐1 km仿真得出風(fēng)險(xiǎn)值和能見度的隨機(jī)數(shù)，采用數(shù)值分析方法對(duì)0—12 km的能見度風(fēng)險(xiǎn)進(jìn)行擬合，利用高斯分布獲得風(fēng)險(xiǎn)分布曲線(圖5)。

圖5 0—12 km的風(fēng)險(xiǎn)擬合曲線Fig.5 Risk fitting curve between 0 and 12 km

從圖5可以得到0—12 km能見度下的廣州港航行風(fēng)險(xiǎn)的估算方程:

y=-0.793+5.909e-0.034(x-6.829)2

(7)

由風(fēng)險(xiǎn)曲線可以得出以下結(jié)果：

(1)在一定范圍內(nèi)，即能見度在0—7 km，船舶航行風(fēng)險(xiǎn)與能見度呈正比，即能見度越高，風(fēng)險(xiǎn)越大；而當(dāng)能見度在7—12 km時(shí)，二者呈反比，即能見度越高，風(fēng)險(xiǎn)越小。

(2)不同的能見度條件航行風(fēng)險(xiǎn)也不同。當(dāng)能見度很低時(shí)(如0 km

(3)當(dāng)能見度處于受限水平時(shí)(如6 km

3 結(jié)論與討論

(1)廣州港航行事故出現(xiàn)的頻率表現(xiàn)為能見度不良時(shí)最小，良好次之，而能見度受限時(shí)出現(xiàn)的頻率幾乎為前兩者的1.3—2.5倍。說明在能見度特別差時(shí)，廣州港船舶航行事故出現(xiàn)的幾率并不高。廣州港航行事故造成的后果嚴(yán)重性表現(xiàn)為能見度不良時(shí)后果最嚴(yán)重，能見度受限時(shí)次之，能見度良好時(shí)后果最小。廣州港航行的平均風(fēng)險(xiǎn)值表現(xiàn)為能見度不良時(shí)風(fēng)險(xiǎn)最小，良好次之；能見度受限時(shí)風(fēng)險(xiǎn)最大，近似為前兩者的1.7—2.4倍左右。

(2)在0—7 km內(nèi)，能見度處于較低水平時(shí)，船舶航行風(fēng)險(xiǎn)與能見度呈正比，即能見度越低，風(fēng)險(xiǎn)越??；能見度在7—12 km內(nèi)，二者呈反比，即能見度越高，風(fēng)險(xiǎn)越小。

(3)不同能見度下船舶航行的風(fēng)險(xiǎn)狀況不一樣，能見度在受限水平時(shí)，風(fēng)險(xiǎn)相對(duì)最高；能見度處于不良與良好水平時(shí)，風(fēng)險(xiǎn)較小。由此可知，能見度受限時(shí)，廣州港相關(guān)部門在日常引航、調(diào)度及風(fēng)險(xiǎn)控制中應(yīng)該更為重視。

(4)廣州港船舶事故發(fā)生的可能性采用的是Beta分布模型，后果采用的是Gamma分布模型進(jìn)行MC仿真，需要進(jìn)一步研究其他概率分布模型從而獲取較小的偏差和最佳的擬合風(fēng)險(xiǎn)；同時(shí)本文僅考慮能見度單因子對(duì)船舶航行風(fēng)險(xiǎn)的影響，實(shí)際航行中容易受到強(qiáng)降水、雷電、大風(fēng)及其伴隨的大浪等多種惡劣天氣的影響，這些需今后進(jìn)一步探討。