ACTA Scientiarum Naturalium Universitatis Pekinensis
我国石油化工行业事故风险分析
赵岩 徐建华
北京大学环境科学与工程学院, 北京 100871; † 通信作者, E-mail: jianhua.xu@pku.edu.cn
摘要 为了针对石油化工行业事故风险制定管理政策, 需要了解事故的发生状况和地区差异。通过多个来源, 收集我国石油化工行业 2006—2015 年发生的 1479 起事故数据, 用分层贝叶斯模型估算事故的发生频率和事故后果的统计分布函数的参数。结果表明, 我国石油化工行业事故频率服从均值参数为 156.00 的泊松分布, 事故的死亡人数服从尺度参数为 1.40、形状参数为 0.36 的广义帕累托分布, 我国不同地区间石油化工行业的事故频率和事故后果分布有较大差异。关键词 石油化工事故; 事故风险; 风险分析中图分类号 X915
石油化工行业是以石油和天然气为起始原料的有机合成工业[1], 在国民经济中占有重要的地位,为多个行业提供能源和基础原材料。依据《中国工业统计年鉴》, 我国内地 31 个省份 2015 年石油化工行业中规模以上企业共计 27393 家, 工业销售产值占全国规模以上工业企业销售产值总额的11.4%。同时, 石油化工行业也是事故风险相对较高的行业, 主要事故后果包括人员伤亡、环境污染、经济损失等[2]。例如, 2005 年 11 月 13 日, 中国石油天然气股份有限公司吉林石化分公司双苯厂
发生爆炸, 造成 8 人死亡, 60 人受伤, 直接经济损失近 7000 万元, 并引发松花江水污染事件[3]。不论是通过政府出台规制措施还是通过保险等经济手段管理石油化工行业的事故风险, 都需要对石油化工行业事故风险的基本状况有准确的了解,以便估算管理措施的成本和收益。在风险研究中,风险常被定义为事件发生的概率和事件带来的后果的乘积, 所以在研究石油化工行业的事故风险时,我们也主要关注: 事故发生的频率有怎样的分布规律? 事故带来的后果有怎样的分布规律? 事故发生
的频率和带来的后果在地区间有怎样的差异? 但是, 现有研究没有就这些问题给出答案。
本研究从多个来源收集 2006—2015 年我国石油化工行业的事故数据, 利用分层贝叶斯模型估算事故频率和事故后果的统计分布函数的参数, 以期为我国石油化工行业的事故风险管理工作提供基础信息。
1 研究现状
现有的关于石油化工行业的事故频率的研究,多数采用简单的描述性统计分析方法, 给出一定地理范围内的事故频率随时间的变化状况[4–7]。例如,
[5] Vilchez 等 用柱状图展示 20 世纪初到 1992 年间95 个国家每十年的化工厂事故和危化品运输过程中发生事故的数目, 发现事故频率随时间呈明显增高趋势, 其中 95%的事故集中发生在 1963—1992
[7]年。叶永锋等 绘制了 1974—2010 年我国化工行业事故发生的趋势图, 发现 20 世纪末至 21 世纪初是我国化工事故的高发时期。在一些试图刻画事故频率分布规律的研究中, 大多假设事故发生的频率符合泊松分布, 并根据历史事故信息对统计分布函数的参数进行估计[8–9]。比如, Meel 等[8]假设化工事故的发生符合泊松分布, 并根据美国国家响应中心数据库的历史事故记录中两家企业的化工事故数据, 分别估算这两家企业事故发生频率分布的均值参数。然而, 目前尚无研究对我国石油化工行业的事故频率所呈现的统计分布规律进行量化分析, 给出其统计分布函数并估算参数值。石油化工行业的事故后果主要包括人员伤亡、
[5,7]经济损失和环境破坏三类 。现有关于我国石油化工行业事故后果的研究对于人员伤亡的讨论相对较多, 主要是对比不同地区平均每起事故死亡人数或者同一地区平均每起事故死亡人数随年份的变化情况。比如, He 等[10]以事故的平均死亡人数作为测量指标, 研究我国不同省份化学品事故后果的严重程度, 发现重庆、吉林、湖南等省份是化学品事
[11]故后果较为严重的地区。杨宗政等 给出 1972— 2010 年我国化工行业突发性大气污染事故平均死亡人数的年际变化状况。也有个别研究比较我国与其他国家石油化工行业事故的平均死亡人数。比如, 杜红岩等[12]的研究表明, 2012 年我国石油化工行业事故的平均死亡人数远高于美国。由于事故后 果的严重程度呈“金字塔”型分布[13], 即后果严重的事故发生概率非常低, 后果轻微的事故发生概率则相对较高, 因此仅以平均每起事故死亡人数衡量事故后果, 不能全面地反映石油化工行业事故后果的分布情况。有研究采用 FN (fatality-frequency)曲线描述事故死亡人数与事故频率之间的关系[14–15]。比如, Fabiano 等[15]绘制了 20 世纪 30 年代到 2010年间荷兰石油工业事故的 FN 曲线。到目前为止,对我国石油化工行业的事故后果统计规律的研究甚少。
已有研究表明, 不同地区间石油化工行业事故风险有很大的差异。比如, 针对 20 世纪 95 个国家的化工厂事故和危化品运输过程中发生事故的研究表明, 大部分事故发生在人口稠密地区[5]。杜红岩等[12]基于 2012 年国内外石油化工行业事故信息的分析, 发现发展中国家和最不发达国家化工事故占总事故数目的 73%, 平均每起事故的死亡人数与国家的经济发展水平之间有一定的关联。针对我国1970—2009 年化学品事故的研究表明, 山东、浙江
[10]和江苏是我国化学品事故多发的省份 。虽然上述研究在一定程度上揭示了石油化工行业事故风险的地区差异, 但仍然没有给出各个地区事故分布规律的差异。
综上所述, 对石油化工行业的事故风险管理需要了解该行业的事故风险状况。但是, 已有的关于我国石油化工行业事故风险的研究都流于简单的描述性统计分析, 没有将事故的规律用统计分布函数描述出来。然而, 不论以保险等经济手段还是通过政府规制来管理风险, 这种统计分布函数都是非常重要的基础信息。本研究拟弥补现有研究中这方面的不足。
2 研究方法
本研究从 4 个来源收集我国石油化工行业2006—2015 年的事故数据: 国家安全生产监督管理总局、化学品安全协会、化学品登记中心以及该时段媒体报道的该行业事故信息, 并选择分层贝叶斯模型来估算我国石油化工行业事故风险的统计分布函数的参数。
2.1 数据来源
我国负责事故调查和事故统计分析工作的主管部门是国家安全生产监督管理总局, 其事故查询系
统记录了我国所有行业中较大及较大以上的事故。根据《国民经济行业分类(GB/T 4754—2011)》, 石油化工行业包括石油和天然气开采业、石油及炼焦加工业、化学原料及化学制品制造业。依据此范畴, 我们对安全生产监督管理总局事故查询系统中2006—2015 年的事故进行筛选和整理, 共得到关于225 起事故的记录。安全生产监督管理总局事故查询系统中记录的事故大部分为较大及较大以上的事故(引起 3 人以上死亡), 如果仅针对上述 225 起事故进行分析, 不仅会低估石油化工行业事故的发生频率, 而且会高估石油化工行业事故的平均损失。Kirchsteiger[16]指出, 在企业安全管理中, 一般事故同样可以揭露安全缺陷。为补充石油化工行业的一般事故数据,我们整理了发布在化学品安全协会网站(http://www. chemicalsafety.org.cn/)和化学品事故信息网(http:// accident.nrcc.com.cn:9090/safeweb/)上的事故信息,分别获得 849 和 418 起事故数据。鉴于上述两个网站数据库的建立没有明确的标准, 且数据来源杂乱, 记录形式和内容差异较大,我们在整理石油化工行业相关事故的基础上, 利用爬虫软件在新浪搜索(http://search.sina.com.cn/)中检索相关事故的新闻, 加以验证和补充。选取两组关键词, 分别为“化工、石化”与“事故、爆炸、火灾、泄漏、中毒”, 在实际检索时, 每组关键词任意选取一个进行搭配, 限定关键词位于全文中, 限定检索的时间区间为 2006 年 1 月 1 日至 2015 年 12 月 31日。共检索到约 3.78 万条网页信息, 通过人工筛查整理, 获得关于 598 起事故的记录。
剔除 4种途径收集到的事故信息中的重复事故后, 整理得到 1479 条石油化工事故记录(表 1)。
2.2 统计模型的选择
事故风险可以通过事故频率和事故后果的严重程度进行定量表征。本研究对事故频率和事故后果分布函数的参数估计均采用分层贝叶斯模型。
在贝叶斯方法中, 待估参数被视为随机变量且具有先验分布。贝叶斯参数估计使用概率分布刻画数据模型, 并利用参数的先验分布和样本似然函数,通过贝叶斯定理给出参数的后验分布。分层贝叶斯模型的特点体现在对先验分布的选取采用分层先验, 即对先验分布中的超参数再给定一个超先验分布。在参数估计过程中, 样本数据的结构会使模型纳入额外的变异性(如不同时段或不同来源的数据引起的变异性)。分层贝叶斯模型可以通过分层先验来刻画这些额外的变异性, 增强估计模型的稳健程度[17]。石油化工行业中后果严重的事故发生概率非常低, 因此传统的统计方法并不适用, 而贝叶斯理论即使在样本数据很少的情况下, 仍能对事故风险进行合理的量化, 同时增加概率模型的置信水平[18]。此外, 事故风险存在不确定性, 贝叶斯参数估计可以使用概率分布表达不确定性。由于事故信息来源的多样性以及事故数据较长的时间跨度使得模型存在额外的变异性, 因此本研究在贝叶斯方法的基础上, 选用分层贝叶斯模型。
我们整理了多个来源的事故数据, 所得信息较充分, 比任何合理的先验分布包含的信息都占据优势, 并且我们希望先验分布对后验分布的影响尽可能地小, 因此模型中参数的先验分布均采用无信息
[17]先验分布 。
马尔科夫链的蒙特卡罗方法(Markov Chain Monte Carlo Modeling, MCMC)允许对复杂模型的后验分布进行抽样分析, 因此, 本研究选取采用MCMC 算法的 OPENBUGS 软件对统计模型进行参数估计[19]。
2.2.1 关于事故频率的统计模型
在统计学中, 常用泊松分布描述单位时间内随机事件的发生次数。实际上,可以认为事故是罕见的随机事件, 其频率服从泊松分布:
发生次数[8,20]。泊松分布的参数 λ 为单位时间内事故平均发生次数, 可以据此进行时间趋势和空间差异分析。本研究采用 Eckle 等[21]对石油消费链中事故频率分析的模型来模拟石油化工事故的频率分布和时间趋势。事故频率服从泊松分布, 其均值参数λt 受时间变量 t影响。模型假定事故频率的均值参数 λt以指数形式增加或减少:
log(t ) 0 1t， (2)其中, λ1为时间趋势参数。
上述模型中, 假定参数 λ0, λ1服从正态分布 λ0 ~ N(μ0, σ0), λ1 ~ N(μ1, σ1), 其中参数 μ0 和 μ1的先验分布均服从正态分布 N(0, 0.01), σ0 和 σ1 的先验分布均服从伽马分布 Ga(0.1, 0.1)。
2.2.2 关于事故后果的统计模型
对石油化工行业事故后果的研究中, 我们更关心的是重大事故。根据事故金字塔理论[13], 后果严重的事故发生概率非常低, 使得事故的后果分布呈现厚尾特征, 是非正态的。在对分布的尾部特征的研究中, 极值理论是一种有效且相当稳定的工具。依据极值理论, 观察值中所有达到或超过某给定阈值的数据近似地服从广义帕累托分布。广义帕累托分布广泛用于金融市场、保险索赔、自然灾害
[22–23]等 领域极端事件的模拟研究。本研究采用广义帕累托分布模拟石油化工行业事故的损失分布。
在整理的事故信息中, 仅有少数事故记录给出直接经济损失和环境影响。另外, 不同数据来源中的事故受伤人数存在误差。因此, 在分布模拟中,事故后果仅考虑死亡人数。利用广义帕累托分布模拟死亡人数分布, 其累积概率函数为
1 F ( y ) 1 1 , y , (3) y
其中, ξ 为形状参数, φ 为尺度参数, θ为阈值参数。
为了模拟形状参数和尺度参数的先验分布, 应使 3 个参数彼此独立, 进行参数变换[21]:
, (1 )。(4)上述模型中, 假定 ξ和修正的尺度参数 χ 服从正态分布 ξ ~ N(μξ, σξ)和 log(χ) ~ N(μχ, σχ), 其中参数 μξ 和 μχ的先验分布均服从正态分布N(0, 0.01), σξ和 σχ的先验分布均服从伽马分布Ga(0.1, 0.1)。
3 结果分析3.1 事故的整体状况
图 1 给出我国石油化工行业 2006—2015 年的年事故数目(用于表征事故频率)和年死亡人数总数(用于表征事故后果), 其中 2006 年事故发生次数较低, 可能是由于部分一般事故未记录或未经新闻报道。每年事故死亡人数虽有所波动, 但整体呈下降趋势。
3.2 事故频率分布及地区差异
由于事故信息来源的局限性导致整理的 2006年事故数目明显低于实际水平, 为保证分析结果的可靠性, 剔除 2006 年的数据信息。以每年的事故数目作为我国石油化工行业的事故频率, 利用分层贝叶斯模型对我国石油化工行业事故频率的泊松分布进行参数估计, 其概率密度函数为156x 156 f ( x ) e , x 0, 1, x!其中, 频率的均值参数为 156.00。模型中对时间趋势参数的估计值为−0.01, 表明 2007—2015 年我国石油化工行业事故频率呈下降趋势, 但降幅较小(图 2)。
如何测量风险决定了可以从风险分析中获得哪些信息, 以及得到的结论是否合理。Wilson 等[24]对美国煤矿事故的风险分析证明, 不同的风险度量指标会使风险分析产生不同的结论。为全面地比较我国石油化工行业事故频率的地区差异, 我们选取事故频率、每千家石油化工企业事故频率和每千亿元行业销售产值事故频率 3 个指标, 对不同省份的 3个指标进行参数估计。由于 2011 年规模以上企业标准从年主营业务收入 500 万元提高到 2000 万元,故本研究基于 2011—2015 年数据分析我国不同地区石油化工行业的事故频率。
按照事故频率递增顺序, 图 3~5 给出我国石油化工行业不同省份的事故频率参数、每千家企业事故频率参数和每千亿元销售产值事故频率参数的均值、5%分位点和 95%分位点。若仅考虑事故发生情况, 我国石油化工行业事故主要集中于江苏、广东、浙江、山东等东部沿海省份。从每千家企业事
故发生频率来看, 石油化工行业事故频率较高的省份主要集中在西北地区, 如陕西、内蒙古、甘肃和新疆等。从每千亿元销售产值事故发生频率来看,石油化工行业事故频率的地区分布没有明显的规律, 较高的省份有宁夏、浙江、内蒙古和安徽。
3.3 事故后果分布及地区差异
如前所述, 在石油化工行业事故后果的分布中,我们更关心重大事故的后果, 即全部事故后果分布的尾部特征。在本研究整理的全部事故中, 有 30%的石油化工事故引起人员死亡。在事故后果仅考虑死亡人数的情况下, 一人死亡已是较严重的事故后果, 因此事故死亡人数的广义帕累托分布阈值参数选为 1。通过贝叶斯参数估计, 我国石油化工行业2006—2015 年死亡事故中每次事故死亡人数的广义帕累托分布为1 0.36 y 1 0.36 , y 1, F ( y ) 1 1 1.40
其中, 尺度参数 φ 为 1.40, 形状参数 ξ 为 0.36。图6 以 FN 曲线的形式展示我国石油化工行业 2006— 2015年事故死亡人数的分布情况。已有研究假定, 石油化工行业事故后果可能与
[12]当地的经济发展水平有一定的关联 。在此基础上, 本研究着重关注不同经济发展水平地区的事故后果具有怎样的特征。基于 2006—2015 年各省人均国内生产总值的年均值和年均增长率, 通过聚类分析, 将我国 29 个省份(不含西藏自治区、海南省和港澳台地区)分为 4 组(表 2)。第一组和第二组省份经济发展水平较高, 人均国内生产总值的年均值在 3.4 万以上。第三组省份经济发展速度较快, 人均国内生产总值的年均增速在 14.7%以上。对各组数据采用贝叶斯方法估计其参数。
由于第一组省份死亡事故数量较少, 模型未给出有效的参数估计结果, 其他组省份石油化工行业事故死亡人数分布的尺度参数 φ、形状参数 ξ 和阈值参数 θ 见表 3。图 7 以 FN 曲线的形式分别展示上述地区事故死亡人数的分布情况。为了比较不同地区事故后果的严重程度, 根据死亡人数的概率分 布, 计算单起事故平均死亡人数和 1% 概率的极端事故死亡人数两个指标。从单起死亡事故平均死亡人数看, 在我国经济发展水平较高的东部地区, 石油化工行业事故后果更严重; 在经济发展水平一般的其他省份中, 后果更严重的极端事故可能与经济的快速发展有关。
4 结论与讨论
本研究综合国家安全生产监督管理总局事故查询系统、化学品安全协会网站、化学品事故信息网以及网络新闻检索 4 种信息来源, 整理得到我国石油化工行业 2006—2015 年发生的 1479 条事故信息, 利用贝叶斯理论模拟了我国石油化工行业事故频率的泊松分布和后果损失的广义帕累托分布, 给
出相应分布的后验参数估计。
根据参数估计结果, 我国石油化工行业的事故频率服从均值参数为156.00 的泊松分布, 每起死亡事故的死亡人数服从尺度参数为 1.40, 形状参数为 0.36 的广义帕累托分布。另外, 时间趋势参数的估计结果表明, 我国石油化工行业 2007—2015 年事故频率略有下降趋势。
受石油化工行业地域分布影响, 我国石油化工行业事故主要集中于江苏、浙江、山东等东部沿海省份, 每年事故的频率均值均超过 15。从每千家企业事故发生频率来看, 石油化工行业事故频率较高的省份主要集中在西北地区, 如陕西、内蒙古、甘肃等, 这在一定程度上也反映出我国西北地区石油化工行业的企业安全管理水平较差。从每千亿元销售产值事故发生频率来看, 石油化工行业事故频率的地区分布没有明显的规律, 频率较高的省份有宁夏、浙江、内蒙古和安徽。
在死亡事故后果方面, 我国经济发展水平较高的东部地区石油化工行业事故后果最严重, 单起死亡事故平均死亡人数和极端事故死亡人数均为最多。杜红岩等[12]对国内外石油化工行业事故的研究认为, 较高的经济发展水平会带来安全管理水平的改善, 而本研究表明, 我国经济发展水平较高的 地区事故后果更严重, 可能是由于我国经济发展处于转型期间, 这段时期安全管理水平的发展相对较慢, 而经济的快速发展导致行业规模不断扩大, 事故不断增多, 极端事故时有发生。在其他经济发展水平一般的省份中, 经济发展较快的地区极端事故的死亡人数更多, 也在一定程度上说明, 经济的快速发展导致石油化工行业安全管理水平的发展相对滞后, 致使极端事故后果更严重。为控制并降低我国石油化工行业的事故风险,可以从事故频率和事故后果两方面采取相应的治理措施: 一方面, 通过加强企业的安全管理工作, 降低事故发生频率, 例如工人的安全培训教育、生产设备的定期检修和维护、安全设备的配备及有效使用; 另一方面, 通过完善相应的应急措施, 减轻或修复事故的负面后果, 如厂内和厂外应急计划的完备及有效执行、事故区域的修复及重建、医疗应急响应机制的建立等。
参考文献
[1] 吴建卿. 石油化工的主要作用[EB/OL]. (2014–10– 29) [2017–04–12]. http://www.cnpc.com.cn/syzs/lyhg/ 201410/3fb1385b9ccd42faa0a37b172315e87f.shtml [2] 黄岳元, 保宇. 化工环境保护与安全技术概论. 北
京: 高等教育出版社, 2006 [3] 国家安全生产监督管理总局. 国务院对中石油吉林石化分公司双苯厂“11·13”爆炸事故及松花江水污染事件作出处理[EB/OL]. (2006–11–24) [2017–04– 12]. http://www.chinasafety.gov.cn/zuixinyaowen/2006-11 /24/content_205732.htm [4] Konstandinidou M, Nivolianitou Z, Markatos N, et al. Statistical analysis of incidents reported in the Greek petrochemical industry for the period 1997– 2003. Journal of Hazardous Materials, 2006, 135 (1): 1–9 [5] Vilchez J A, Sevilla S, Montiel H, et al. Historical analysis of accidents in chemical plants and in the transportation of hazardous materials. Journal of Loss Prevention in the process Industries, 1995, 8(2): 87– 96 [6] Nivolianitou Z, Konstandinidou M, Michalis C. Statistical analysis of major accidents in petrochemical industry notified to the major accident reporting system (MARS). Journal of Hazardous Materials, 2006, 137(1): 1–7 [7] 叶永峰, 夏昕, 李竹霞. 化工行业典型安全事故统计分析. 工业安全与环保, 2012, 38(9): 49–51 [8] Meel A, O’neill L M, Levin J H, et al. Operational risk assessment of chemical industries by exploiting accident databases. Journal of Loss Prevention in the Process Industries, 2007, 20(2): 113–127 [9] Meel A, Seider W D. Plant-specific dynamic failure assessment using Bayesian theory. Chemical Engineering Science, 2006, 61(21): 7036–7056 [10] He G, Zhang L, Lu Y, et al. Managing major chemical accidents in China: towards effective risk information. Journal of Hazardous Materials, 2011, 187(1): 171–181 [11] 杨宗政, 熊发, 袁雪竹, 等. 化工行业突发大气污染事故统计与趋势分析. 天津科技大学学报, 2015, 30(4): 8–12 [12] 杜红岩, 王延平, 卢均臣. 2012 年国内外石油化工行业事故统计分析. 中国安全生产科学技术, 2013, 9(6): 184–188 [13] Rausand M. Risk assessment: theory, methods, and applications. Hoboken: John Wiley & Sons, 2013 [14] Khan F I, Abbasi S A. Major accidents in process industries and an analysis of causes and consequences. Journal of Loss Prevention in the Process Industries, 1999, 12(5): 361–378 [15] Fabiano B, Currò F. From a survey on accidents in the downstream oil industry to the development of a detailed near-miss reporting system. Process Safety and Environmental Protection, 2012, 90(5): 357–367 [16] Kirchsteiger C. Summary of JRC/ESREDA seminar on safety investigation of accidents, 12–13 may 2003, European Commission, DG JRC-IE, Petten, the Netherlands. Journal of Hazardous Materials, 2004, 111(1/2/3): 167–170 [17] Kelly D, Smith C. Bayesian inference for probabilistic risk assessment: a practitioner's guidebook. London: Springer, 2011 [18] Robert C. The Bayesian choice: from decisiontheoretic foundations to computational implementation. New York: Springer, 2007 [19] Lunn D, Spiegelhalter D, Thomas A, et al. The BUGS project: evolution, critique and future directions. Statistics in Medicine, 2009, 28(25): 3049–3067 [20] Carson J, Mannering F. The effect of ice warning signs on ice-accident frequencies and severities. Accident Analysis & Prevention, 2001, 33(1): 99–109 [21] Eckle P, Burgherr P. Bayesian data analysis of severe fatal accident risk in the oil chain. Risk Analysis, 2013, 33(1): 146–160 [22] Embrechts P, Schmidli H. Modelling of extremal events in insurance and finance. Zeitschrift für Operations Research, 1994, 39(1): 1–34 [23] Coles S, Casson E. Extreme value modelling of hurricane wind speeds. Structural Safety, 1998, 20(3): 283–296 [24] Wilson R, Crouch E A C. Risk-benefit analysis. Cambridge: Harvard University Press, 2001