徐广波, 胡甚平
(1. 上海海事大学 商船学院,上海 201306; 2. 江阴海事局,江苏 江阴 214431)
随着经济建设发展和科技水平的不断进步,航运业得到快速发展,船舶逐渐向大型化、高速化、综合化、智能化方向迈进.大量频繁活跃的船舶使港口、航道等水域交通流量急增,航行环境日益复杂.[1]因此,有必要对港口水域的交通风险进行分析和评估,为港航安全管理部门制定有针对性的预防措施提供依据.目前,国内外已有很多学者采用定性、定量以及定性定量相结合的方法从不同角度研究水上交通风险问题,如综合安全评估(FSA)法[2]、灰色系统理论法[3]、灰色马尔可夫链方法[4]等.这些方法在一定程度上提供定性和定量的决策依据,然而,对动态、线性的数据进行实时评估,需要更多的数据推理及模拟实验分析发展趋势.
本文在港口交通风险定量化评估的基础上给出交通事故率和事故后果的贝叶斯概率统计,接着建立港口交通系统风险的蒙特卡罗仿真模型,并通过马尔可夫链蒙特卡罗(Markov Chain Monte Carlo,MCMC)方法对风险模型的参数进行推断与优化,最后进行仿真实验.
港口交通系统风险[5]是指在某一特定的客观状态下港口交通系统中人、船(货)、环境受到伤害的可能性和这种伤害的严重程度,可表示为
R|S=f(P,C)|S
(1)
式中:P为某一时间事故发生的可能性;C为事件发生的后果;S代表某一特定的客观状态;R|S为在S状态下分析对象的风险度;f为关于P和C的实数函数.可能性指事故发生的机会,用于描述概率或频率的性质.概率是理论值,由事件的本质决定,只能取唯一值,它能精确反映事件发生可能性的大小.称随机事件在一定时间内统计取得的发生次数为频率. 频率是试验值或使用时的统计值,具有随机性,可能取多个数值,因此只能近似反映事件发生可能性的大小.可能性
(2)
式中:i∈[1,n],i∈Z;P的物理意义是指在选取的样本中发生事故的频率;ai为年事故数;ai为每年船舶活动量.
后果是描述有害事件或非正常事件发生所造成损害的程度.在实际事故后果分析中,由于考虑不同的风险,常采用“事故等效后果”衡量.“事故等效后果”是经归一化以后事故的各种后果总和.
(3)
MCMC方法是最近发展起来的一种简单且行之有效的贝叶斯计算方法.其基本思想是通过建立一个平稳分布为π(z)的马尔可夫链得到π(z)的样本,基于这些样本作各种统计推断.其主要目的是借助贝叶斯概率估计,通过频率数据获取风险事件的概率.
将MCMC算法概括为3步:(1)在z上选一个“合适”的马尔可夫链,使其转移核为p(·,·),这里“合适”的含义主要指π(z)应是其相应的平稳分布;(2)由z中某一点Z(0)出发,用(1)中的马尔可夫链产生点序列Z1,…,Zn;(3)对某个m和n,任一函数f(z)的期望估计
(4)
采用文献[6]中所建立的模型作为风险仿真模型,即把要解决的目标问题抽象成一个概率模型R=f(P,C)(式中R表示风险度;P表示某一事件发生的概率;C表示事件发生的后果).分别对相关资料进行统计,获取P和C的样本信息;通过结合样本信息和模拟风险值确定概率模型;最后选定适当的模拟值个数N与次数M,对所获得的M·N个样本值进行统计分析,得到分布曲线和宏观风险的特征.
主要步骤[7]如下:(1)收集、分析主客观先验信息,确定合适的先验分布形式以及先验参数.(2)结合试验数据确定第i个模型,利用Gibbs抽样对模型的后验进行MCMC模拟.(3)判断马尔可夫链是否已收敛、MC误差是否足够小.如果马尔可夫链已收敛、MC误差足够小则转入下一步,否则需进一步调试模型.重新确定抽样迭代次数及抽样方法等,若效果依然不尽如人意,则返回(2),重新考虑修改先验参数和模型.(4)修改(3)中的模型、选择更高一级的i+1个模型并返回(3),比较模型、选择相对更优秀的后验模型,进行模型的贝叶斯推断,并根据有关准则得出正确的结论.
基于MCMC方法的模型运行流程见图1.
构建风险的贝叶斯多层对数正态模型:
R[i]~d lnorm(μ,σ)//风险服从对数正态分布
R[i]<-P[i]*C[i]//风险由概率和后果确定
P[i]~d lnorm(μ,σ) //可能性服从对数正态分布
C[i]~d lnorm(μ,σ) //后果服从对数正态分布
μ~d norm(0,1.0E-6) //超参数服从正态分布
σ~d gamma(0.01,0.001) //超参数服从伽马分布
sigma<-1/sqrt(σ) //符号函数
图1 基于MCMC方法的模型运行流程
利用WinBUGS软件对各参数进行迭代计算,从而获得其后验分布和数学特征值.
仿真数据检验方法主要采用文献[8]中的变异因数,即通过判断估计值与待定值的收敛性检验仿真数据的可靠性,
V=σ/m
(5)
式中:σ为仿真样本的风险标准差;m为仿真样本的风险平均值.
用我国某港口的交通事故数据[2]进行应用测试.由于事故数据样本量较少,属于小样本事件,为获取港口交通系统风险事件发生的可能性往往只能用频率进行估计,通过蒙特卡罗方法获取概率.后果采用因数进行换算.文献[8]用基于对数正态的蒙特卡罗模型进行风险推理.本文则通过MCMC算法先修正概率,进而对风险事件数据进行仿真.
WinBUGS是英国剑桥公共卫生研究所推出的用MCMC方法来分析复杂统计模型的贝叶斯推理软件.其基本原理就是通过Gibbs sampling和Metropolis算法,从完全条件概率分布中抽样,从而生成马尔可夫链,通过迭代最终估计出优化后的概率参数[9].
在对风险事件样本数据整理之后,得到先验分布数字特征,然后利用WinBUGS软件进行MCMC参数计算[10],获得所建模型各参数的后验分布数字特征,见表1.
表1 各参数的后验分布数字特征
利用表1中的参数进行蒙特卡洛仿真,结果见图2.图中实心点为样本数据,圆圈表示仿真获得的随机数据.
从图2(b)中可以看出,样本量得到有效增加,新增加样本对原样本有较好的覆盖.从图2(c)和2(d)中可见高风险出现的概率明显较小.通过扩充样本信息量,在风险合成后可得到整体风险.
当模拟运行1 000次后(每次抽取样本100个),仿真结果相差极小.经统计可得出该港口水域频率均值为1.279 3×10-4,后果均值为0.764 2,后果均值与实际情况也较为符合,风险平均值为0.975 1×10-4,风险标准差小,为2.518 2×10-4,说明仿真的风险值稳定性很好,仿真结果可信.
图2经优化后的港口水域样本及模拟样本仿真风险分布
利用上面得到的模拟风险值(1 000个),可以画出一个统计直方图,见图3.
(a)直方图 (b)三维离散图
图3多次模拟样本的仿真风险
按照文献[8]的模型,对数正态下的概率参数与MCMC优化参数后的概率参数的风险统计特征值见表2.
表2 风险统计特征值
从图3可得到以下结论:(1)该港口不同水域不同事故类型的平均风险为0.975 1×10-4,风险值标准差较小,稳定性较好.(2)风险值主要集中在0.4×10-4~1.5×10-4之间,峰值主要集中在0.95×10-4处,与风险均值接近重合.该水域的风险值是较为典型的对数正态分布的随机变量.(3)利用统计数据,对风险进行分级.根据二八定律,可将低风险、合理可行风险、高风险等3类风险度的阈值进行重新划分,进一步保证相对风险分级的可行性.(4)与参数优化前MC仿真的结果相比较,参数优化对仿真结果起到一定的优化作用,主要表现在所获得的参数本身具有更强的稳定性,仿真所得的频率均值、风险均值及后果均值与实际结果更为接近,同时风险值区间上下限值有所扩大,较高风险可能性增大,更能提醒相关职能部门和船舶驾驶人员防范交通风险.
从第3.1节风险信息(N=100,M=1 000)中随机抽取100次仿真数据进行变异因数计算,得到图4(其中稳定数据线为样本的变异因数,波动数据线为仿真样本变异因数).
图4 基于WinBUGS软件的仿真风险的变异因数
图4中,虽然仿真样本的数据曲线振荡幅度较为显著,但是样本的统计结果基本落在仿真样本的变异因数的变化范围之内,两种样本的关联性较强.
应用MC方法建立的仿真模型能够有效扩展小样本问题分析数据,为确立风险分级标准提供依据.利用WinBUGS软件通过MCMC方法能够实现对风险模型中概率参数的优化,更加接近风险事件发生的可能性.实例表明,参数优化后的风险模型具有更强的稳定性,与实际结果具有较好的一致性.如何利用MCMC方法实现风险的动态转移是后续的研究工作.
参考文献:
[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] 胡甚平. 海上交通系统风险蒙特卡洛仿真[J].上海海事大学学报, 2011, 32(4): 7-11.
[9] 余芳. 基于 MCMC 方法的 WinBUGS 软件应用[J]. 经营管理者, 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.