徐廣波， 胡甚平

(1. 上海海事大學 商船學院，上海 201306； 2. 江陰海事局，江蘇 江陰 214431)

0 引 言

隨著經(jīng)濟建設發(fā)展和科技水平的不斷進步，航運業(yè)得到快速發(fā)展，船舶逐漸向大型化、高速化、綜合化、智能化方向邁進.大量頻繁活躍的船舶使港口、航道等水域交通流量急增，航行環(huán)境日益復雜.[1]因此，有必要對港口水域的交通風險進行分析和評估，為港航安全管理部門制定有針對性的預防措施提供依據(jù).目前，國內外已有很多學者采用定性、定量以及定性定量相結合的方法從不同角度研究水上交通風險問題，如綜合安全評估(FSA)法[2]、灰色系統(tǒng)理論法[3]、灰色馬爾可夫鏈方法[4]等.這些方法在一定程度上提供定性和定量的決策依據(jù)，然而，對動態(tài)、線性的數(shù)據(jù)進行實時評估，需要更多的數(shù)據(jù)推理及模擬實驗分析發(fā)展趨勢.

本文在港口交通風險定量化評估的基礎上給出交通事故率和事故后果的貝葉斯概率統(tǒng)計，接著建立港口交通系統(tǒng)風險的蒙特卡羅仿真模型，并通過馬爾可夫鏈蒙特卡羅(Markov Chain Monte Carlo，MCMC)方法對風險模型的參數(shù)進行推斷與優(yōu)化，最后進行仿真實驗.

1 港口交通系統(tǒng)風險建模

港口交通系統(tǒng)風險[5]是指在某一特定的客觀狀態(tài)下港口交通系統(tǒng)中人、船(貨)、環(huán)境受到傷害的可能性和這種傷害的嚴重程度，可表示為

R|S=f(P,C)|S

(1)

式中：P為某一時間事故發(fā)生的可能性；C為事件發(fā)生的后果；S代表某一特定的客觀狀態(tài)；R|S為在S狀態(tài)下分析對象的風險度；f為關于P和C的實數(shù)函數(shù).可能性指事故發(fā)生的機會，用于描述概率或頻率的性質.概率是理論值，由事件的本質決定，只能取唯一值，它能精確反映事件發(fā)生可能性的大小.稱隨機事件在一定時間內統(tǒng)計取得的發(fā)生次數(shù)為頻率. 頻率是試驗值或使用時的統(tǒng)計值，具有隨機性，可能取多個數(shù)值，因此只能近似反映事件發(fā)生可能性的大小.可能性

(2)

式中：i∈[1,n],i∈Z；P的物理意義是指在選取的樣本中發(fā)生事故的頻率；ai為年事故數(shù)；ai為每年船舶活動量.

后果是描述有害事件或非正常事件發(fā)生所造成損害的程度.在實際事故后果分析中，由于考慮不同的風險，常采用“事故等效后果”衡量.“事故等效后果”是經(jīng)歸一化以后事故的各種后果總和.

(3)

2 港口交通系統(tǒng)風險的MCMC仿真算法

2.1 MCMC方法及風險仿真模型

MCMC方法是最近發(fā)展起來的一種簡單且行之有效的貝葉斯計算方法.其基本思想是通過建立一個平穩(wěn)分布為π(z)的馬爾可夫鏈得到π(z)的樣本，基于這些樣本作各種統(tǒng)計推斷.其主要目的是借助貝葉斯概率估計，通過頻率數(shù)據(jù)獲取風險事件的概率.

將MCMC算法概括為3步：(1)在z上選一個“合適”的馬爾可夫鏈，使其轉移核為p(·,·)，這里“合適”的含義主要指π(z)應是其相應的平穩(wěn)分布；(2)由z中某一點Z(0)出發(fā)，用(1)中的馬爾可夫鏈產生點序列Z1,…,Zn；(3)對某個m和n，任一函數(shù)f(z)的期望估計

(4)

采用文獻[6]中所建立的模型作為風險仿真模型，即把要解決的目標問題抽象成一個概率模型R=f(P,C)(式中R表示風險度；P表示某一事件發(fā)生的概率；C表示事件發(fā)生的后果).分別對相關資料進行統(tǒng)計，獲取P和C的樣本信息；通過結合樣本信息和模擬風險值確定概率模型；最后選定適當?shù)哪M值個數(shù)N與次數(shù)M，對所獲得的M·N個樣本值進行統(tǒng)計分析，得到分布曲線和宏觀風險的特征.

2.2 基于MCMC的貝葉斯推斷及模型選擇

主要步驟[7]如下：(1)收集、分析主客觀先驗信息，確定合適的先驗分布形式以及先驗參數(shù).(2)結合試驗數(shù)據(jù)確定第i個模型，利用Gibbs抽樣對模型的后驗進行MCMC模擬.(3)判斷馬爾可夫鏈是否已收斂、MC誤差是否足夠小.如果馬爾可夫鏈已收斂、MC誤差足夠小則轉入下一步，否則需進一步調試模型.重新確定抽樣迭代次數(shù)及抽樣方法等，若效果依然不盡如人意，則返回(2)，重新考慮修改先驗參數(shù)和模型.(4)修改(3)中的模型、選擇更高一級的i+1個模型并返回(3)，比較模型、選擇相對更優(yōu)秀的后驗模型，進行模型的貝葉斯推斷，并根據(jù)有關準則得出正確的結論.

基于MCMC方法的模型運行流程見圖1.

2.3 參數(shù)優(yōu)化

構建風險的貝葉斯多層對數(shù)正態(tài)模型：

R[i]～d lnorm(μ,σ)//風險服從對數(shù)正態(tài)分布

R[i]<-P[i]*C[i]//風險由概率和后果確定

P[i]～d lnorm(μ,σ) //可能性服從對數(shù)正態(tài)分布

C[i]～d lnorm(μ,σ) //后果服從對數(shù)正態(tài)分布

μ～d norm(0,1.0E-6) //超參數(shù)服從正態(tài)分布

σ～d gamma(0.01,0.001) //超參數(shù)服從伽馬分布

sigma<-1/sqrt(σ) //符號函數(shù)

圖1 基于MCMC方法的模型運行流程

利用WinBUGS軟件對各參數(shù)進行迭代計算，從而獲得其后驗分布和數(shù)學特征值.

2.4 基于仿真數(shù)據(jù)的檢驗

仿真數(shù)據(jù)檢驗方法主要采用文獻[8]中的變異因數(shù)，即通過判斷估計值與待定值的收斂性檢驗仿真數(shù)據(jù)的可靠性，

V=σ/m

(5)

式中：σ為仿真樣本的風險標準差；m為仿真樣本的風險平均值.

3 應用算例

用我國某港口的交通事故數(shù)據(jù)[2]進行應用測試.由于事故數(shù)據(jù)樣本量較少，屬于小樣本事件，為獲取港口交通系統(tǒng)風險事件發(fā)生的可能性往往只能用頻率進行估計，通過蒙特卡羅方法獲取概率.后果采用因數(shù)進行換算.文獻[8]用基于對數(shù)正態(tài)的蒙特卡羅模型進行風險推理.本文則通過MCMC算法先修正概率，進而對風險事件數(shù)據(jù)進行仿真.

WinBUGS是英國劍橋公共衛(wèi)生研究所推出的用MCMC方法來分析復雜統(tǒng)計模型的貝葉斯推理軟件.其基本原理就是通過Gibbs sampling和Metropolis算法，從完全條件概率分布中抽樣，從而生成馬爾可夫鏈，通過迭代最終估計出優(yōu)化后的概率參數(shù)[9].

在對風險事件樣本數(shù)據(jù)整理之后，得到先驗分布數(shù)字特征，然后利用WinBUGS軟件進行MCMC參數(shù)計算[10]，獲得所建模型各參數(shù)的后驗分布數(shù)字特征，見表1.

表1 各參數(shù)的后驗分布數(shù)字特征

利用表1中的參數(shù)進行蒙特卡洛仿真，結果見圖2.圖中實心點為樣本數(shù)據(jù)，圓圈表示仿真獲得的隨機數(shù)據(jù).

3.1 港口交通系統(tǒng)風險MCMC仿真結果分析

從圖2(b)中可以看出，樣本量得到有效增加，新增加樣本對原樣本有較好的覆蓋.從圖2(c)和2(d)中可見高風險出現(xiàn)的概率明顯較小.通過擴充樣本信息量，在風險合成后可得到整體風險.

當模擬運行1 000次后(每次抽取樣本100個)，仿真結果相差極小.經(jīng)統(tǒng)計可得出該港口水域頻率均值為1.279 3×10-4，后果均值為0.764 2，后果均值與實際情況也較為符合，風險平均值為0.975 1×10-4，風險標準差小，為2.518 2×10-4，說明仿真的風險值穩(wěn)定性很好，仿真結果可信.

圖2經(jīng)優(yōu)化后的港口水域樣本及模擬樣本仿真風險分布

利用上面得到的模擬風險值(1 000個)，可以畫出一個統(tǒng)計直方圖，見圖3.

(a)直方圖 (b)三維離散圖

圖3多次模擬樣本的仿真風險

按照文獻[8]的模型，對數(shù)正態(tài)下的概率參數(shù)與MCMC優(yōu)化參數(shù)后的概率參數(shù)的風險統(tǒng)計特征值見表2.

表2 風險統(tǒng)計特征值

從圖3可得到以下結論：(1)該港口不同水域不同事故類型的平均風險為0.975 1×10-4，風險值標準差較小，穩(wěn)定性較好.(2)風險值主要集中在0.4×10-4～1.5×10-4之間，峰值主要集中在0.95×10-4處，與風險均值接近重合.該水域的風險值是較為典型的對數(shù)正態(tài)分布的隨機變量.(3)利用統(tǒng)計數(shù)據(jù)，對風險進行分級.根據(jù)二八定律，可將低風險、合理可行風險、高風險等3類風險度的閾值進行重新劃分，進一步保證相對風險分級的可行性.(4)與參數(shù)優(yōu)化前MC仿真的結果相比較，參數(shù)優(yōu)化對仿真結果起到一定的優(yōu)化作用，主要表現(xiàn)在所獲得的參數(shù)本身具有更強的穩(wěn)定性，仿真所得的頻率均值、風險均值及后果均值與實際結果更為接近，同時風險值區(qū)間上下限值有所擴大，較高風險可能性增大，更能提醒相關職能部門和船舶駕駛人員防范交通風險.

3.2 仿真樣本差異性分析

從第3.1節(jié)風險信息(N=100，M=1 000)中隨機抽取100次仿真數(shù)據(jù)進行變異因數(shù)計算，得到圖4(其中穩(wěn)定數(shù)據(jù)線為樣本的變異因數(shù)，波動數(shù)據(jù)線為仿真樣本變異因數(shù)).

圖4 基于WinBUGS軟件的仿真風險的變異因數(shù)

圖4中，雖然仿真樣本的數(shù)據(jù)曲線振蕩幅度較為顯著，但是樣本的統(tǒng)計結果基本落在仿真樣本的變異因數(shù)的變化范圍之內，兩種樣本的關聯(lián)性較強.

4 結束語

應用MC方法建立的仿真模型能夠有效擴展小樣本問題分析數(shù)據(jù)，為確立風險分級標準提供依據(jù).利用WinBUGS軟件通過MCMC方法能夠實現(xiàn)對風險模型中概率參數(shù)的優(yōu)化，更加接近風險事件發(fā)生的可能性.實例表明，參數(shù)優(yōu)化后的風險模型具有更強的穩(wěn)定性，與實際結果具有較好的一致性.如何利用MCMC方法實現(xiàn)風險的動態(tài)轉移是后續(xù)的研究工作.

參考文獻：

[1] 徐廣波, 軒少永, 尤慶華. 基于熵權的模糊集對模型在港口水域通航風險評價中的應用[J].上海海事大學學報, 2012, 33(1): 7-11.

[2] 方泉根, 胡甚平. FSA 在船舶引航風險評估中的應用[J]. 哈爾濱工程大學學報, 2006, 27(3): 329-334.

[3] 周麗麗, 胡甚平. 船舶引航風險成因灰色綜合評價模型[J]. 上海海事大學學報, 2008, 29(3): 21-25.

[4] 趙佳妮, 吳兆麟. 基于灰色馬爾可夫模型的水上交通事故預測[J]. 大連海事大學學報, 2005, 31(4): 15-18.

[5] 胡甚平. 不確定性信息下的海上交通風險評估方法及應用研究[D]. 上海: 上海海事大學, 2010.

[6] 高帥. 基于蒙特卡羅方法的沿海水上交通風險仿真[D]. 上海: 上海海事大學, 2011.

[7] 林靜. 基于MCMC的貝葉斯生存分析理論及其在可靠性評估中的應用[D]. 南京: 南京理工大學, 2008.

[8] 胡甚平. 海上交通系統(tǒng)風險蒙特卡洛仿真[J].上海海事大學學報, 2011, 32(4): 7-11.

[9] 余芳. 基于 MCMC 方法的 WinBUGS 軟件應用[J]. 經(jīng)營管理者, 2010(15): 15.

[10] LUNN D J, THOMAS A, BEST N,etal. WinBUGS-a Bayesian modelling framework: concepts, structure, and extensibility[J]. Stat & Computing, 2000, 10(4): 325-337.