基于贝叶斯理论的单一机型事故预测方法研究

崔利杰 丁野 任博 李泽

贝叶斯统计方法是一种以贝叶斯公式为核心,以先验信息和后验信息为综合依据,以“辩证”推断为主要特征的统计方法[1].运用贝叶斯方法进行统计推断时,在获取数据信息前,人们将主观经验与样本信息进行综合作为先验信息,实现样本信息的更新,再根据贝叶斯定理推导出参数的后验概率分布,最后以后验分布作为基础进行参数推测.由于利用了先验知识,贝叶斯推断对小样本数据具有较好的统计推断效果[2−4].

贝叶斯理论经常被用于各类事故预测.如Shahram Heydari等人改进了传统的道路交通事故预测方法[5−6],给出了小子样条件下交通事故贝叶斯预测模型,极大地缩小了预测误差;LouiseK.Lloyd和JonathanJ.Forster通过分析交通流量数据[7],以注册汽车型号、路况和汽车行驶里程等数据作为分布信息,建立了一套预测交通事故发生率的贝叶斯数学方法,取得了良好的效果;贾丰源等人基于贝叶斯理论对上海快速路采集的事故数据和相关检测器数据进行数学建模[8],事故预测准确率达到了82.78%.

随着现代航空装备的技术含量不断增加,装备的制造维护成本也水涨船高.对于各国空军而言,由于装备事故所造成的损失也愈加难以承担.因此,能否有效对航空装备的事故进行预测,并及时采取措施加以防范,是摆在装备研究人员面前的重要课题.Ralph R Fullwood和Robert E Hall等人结合航空服务困难报告中的事故数据[9],运用线性回归的方法对航空安全趋势进行了预测;中国民航大学的王永刚等人针对民航事故征候的预测[10],建立了基于灰色神经网络的预测模型;端木京顺和甘旭升等人将灰色模型和均生函数结合[11],建立了一套航空装备的事故预测理论.以上模型虽然能够取得不错的事故预测效果,但均存在对事故数据变化不敏感的问题,当事故数据更新时,无法对模型进行快速修正,容易出现较大偏差.另外,这些预测模型对随机波动较大的数据序列预测效果不佳.而贝叶斯方法由于结合样本信息不断对预测模型进行更新,因而能够及时反映出数据的变化趋势,迅速对模型参数进行修正,提高预测准确率.本文将以贝叶斯统计决策理论作为基础,针对军用航空装备的事故预测建立数学模型,并以美国某型战斗机为例,对其事故发生情况进行预测.

1 军用飞机事故特点分析

要对军机事故进行预测,首先要对军机的事故数据特点有定性认识.一般来说,投入使用后,其事故规律一般遵循类似抛物线的“先升后降”的趋势.这是因为装备在刚交付部队时,由于装备本身的技术不成熟以及使用维护人员对装备不了解、不熟悉,导致在这个期间装备事故呈逐渐上升的趋势.随着设计方对装备的不断改进和使用方对装备性能的掌握,装备的事故数据会逐渐下降.通过查阅文献[12],我们对美国的F-15与欧洲的“狂风”战机的事故数据(不计作战损失)进行统计,得出其事故数据曲线,在一定程度上印证这种军机事故的发生发展趋势.(见图1、图2).

另外,与民航客机常年飞固定航线不同,由于事故发生具有极强的偶然性,且执行任务多样、飞行环境复杂,导致军用飞机的飞行事故数据在呈现上述总体趋势的同时还具有较大的波动性,从图中也可以看出,事故数据浮动的随机性较强.

基于以上思想,我们认为军机事故的发生发展是遵循一定数学规律的,通过对事故规律进行数学描述,就能够以此为基础建立起关于事故预测的数学模型.

2 贝叶斯事故预测理论数学模型

2.1 贝叶斯统计决策理论基础

假定θ是参数空间Θ上的随机变量,它的概率分布为参数 Θ 的先验分布,记为 {π(θ):θ∈Θ}.又因为θ是随机变量,因此在参数统计模型中,样本X=(X1,X2,···,X n)T的分布族应理解为条件分布族{F(x|θ):θ∈ Θ},不论连续型还是离散型变量,样本分布族均可表为条件密度函数{f(x|θ):θ∈Θ},其中f(x|θ)为样本密度函数.先验分布 {π(θ):θ∈ Θ}与样本分布族{f(x|θ):θ∈Θ}构成贝叶斯参数统计模型.

给定了贝叶斯参数统计模型,即给定了先验分布 {π(θ):θ∈ Θ}与样本分布族 {f(x|θ):θ∈ Θ}.其中,π(θ)为密度函数.通过给定的模型,可以确定(θ,X)的联合分布.它们的联合密度函数为

X的边缘密度函数为

在X=x时,可以求得θ的条件密度函数为

式(3)即贝叶斯公式.h(θ|x)反映了得到样本观测x后θ的后验分布.h(θ|x)为后验密度函数.Bayes学派认为:后验分布族{h(θ|x),θ∈Θ}是统计推断出发点,对参数θ作的任何统计推断必须且只能基于θ的后验分布[2].

得到了参数θ的后验分布,还需要对未知参数θ给出一个合适的估计量,作为问题的解.根据贝叶斯统计决策理论[13−14],用样本对总体进行推断时,误差总是不可避免的,因此还应确定损失函数定量表达统计过程中的优劣程度.

2.2 贝叶斯事故预测数学模型的建立

2.2.1 计算对象的确定

本文中,建立数学模型的基本思想是“拟合+分布预测”的方法.想要找到事故数据的数学规律,就必须对事故数据进行统计学规律的探索.最常用的方法是对已有的事故数据进行拟合,得到一条能够由数学方程描述的曲线.前文中,分析了军机事故数据大体符合“先升后降”的趋势,这种趋势基本符合二次曲线的形状特点,因此,首先对已有的事故数据进行二次拟合.

假设掌握了某型飞机n年的事故数据

对事故数据进行拟合后,得到该机型的“事故－年份”拟合曲线以及拟合方程,这条曲线反映了该型飞机的事故趋势.同时,根据曲线方程,可以得到每年事故的拟合值

据此能够计算每年事故的实际发生数与拟合值的残差

在对具体数据进行分析时,可以对数据拟合状况进行拟合优度检验,以确认选择的拟合方式是合适的.这里s i表示的是实际值与拟合值之间的离散程度.由此得到了表达数据拟合效果的统计量S.在对第(n+m)年的事故进行预测时(m=1,2,···),该年发生事故的拟合值µn+m可以由拟合方程得出,只要确定了该年事故数据的残差s n+m,就能够给出该年事故数据的确切范围,从而对事故做出推断.

要想对残差统计量S作出推断,有必要了解其统计规律.以“狂风”战斗机的事故数据为例,对图2中的事故数据进行拟合,得到拟合曲线,并据此计算每一年的数据残差s i.为得到军机残差分布的统计规律,将所统计的残差值在频率分布直方图上表示出来.见图3.

通过“狂风”战机的残差分布图可以看出,事故数据拟合值与实际值的残差s i分布的随机性较强,分布规律基本符合均值为0的正态分布.除此之外,我们还对F-15、米格-29、F-14等多型军机的事故数据进行了拟合统计,发现在样本容量足够大的条件下,数据的分布都更加接近这种分布规律.因此在这里作出假设:在对军机进行统计规律分析时,认为统计量S关于拟合数据点是满足正态分布的,且以每一个拟合点µi(i=1,2,···,n)作为数据统计点,s i关于µi都满足S～N(0,σ2)的正态分布律.这里,标准差σ为未知量.由此,我们对于未知量S的求解就聚焦到σ的确定上来.

2.2.2 先验分布的确定

通过上文的分析,明确了研究的对象σ.这里就将σ作为先验信息,通过样本观测,对σ进行分析,进而对统计量S进行确定.

在很多情况下,军方关于军机事故的数据都只有具体事故数,对于数据的其他信息我们所知甚少.根据已有数据提供的信息,无法确定σ的具体分布形式,根据贝叶斯假设,先验分布是无信息先验分布时,在参数的取值区域内“均匀分布”[2],即

其中,Σ为σ的参数空间,b为σ的取值上限,a为σ的取值下限,即σ在区间[a,b]上服从均匀分布.接下来要确定σ的取值区域.由于我们掌握了前n年的数据残差s i(i=1,2,···,n),统计量S满足分布S～N(0,σ2).根据正态分布的“3σ”法则,变量s i落在[−3σ,3σ]区间的概率为99.7%,因此我们近似认为的取值全部落入[−3σ,3σ]区间中,又因为σ＞0,故可以得到σ的区间为(0,|s i|max/3].因此,就确定了σ的先验分布形式

2.2.3 后验分布的确定

随着时间的推移,军机的事故数据也在不断积累.意味着我们掌握的先验信息已经逐渐“过时”了,原有的预测模型无法保证对事故预测的准确性.这时需要结合更新的事故数据,综合先验信息得到后验分布,调整样本信息的参数,保证事故预测严格遵循事故数据的统计规律,从而提高预测精度.

假定现在更新了m年的事故数据(x n+1,x n+2,···,x n+m),新样本的加入可能会改变 “事故－年份”曲线的拟合状况,因此需要对数据样本(x1,x2,···,x n,···,x n+m)进行重新拟合,得到了一组新的拟合值相应地,也得到了一组关于样本观测值的残差2,···,n+m).这里,将这m年的残差看作独立同分布样本.由2.2.1的分析可知应当同样满足分布律S～N(0,σ2),则这m年的残差数据的联合分布概率密度为:

那么由贝叶斯理论,根据式(3)有

其中,Θ为θ的参数空间.

分母积分后为常量,故有

由此给出了σ的后验分布形式,根据式(12),给出了对σ进行信息更新后的数学表达,只要据此求出σ的取值,就能够对预测数据的残差进行计算,进而给出预测数据的取值范围.

2.2.4 损失函数的确定

要对σ进行预测,需要能够得到σ的具体取值.因此在这里对σ采取点估计的估计方法.二次损失函数是参数点估计中常用的一种损失函数[15],其形式为

其中,θ为参数的实际值,a为估计值.在二次损失函数的条件下,参数的贝叶斯估计为

根据式(14),可以确定σ的点估计值,进而确定预测年份s i范围,对事故发生区间进行预测.

2.2.5 预测模型建立流程

根据前文分析,建立起预测模型的建立流程,见图4.

3 典型军机事故预测算例

某型战斗机是外军的主力机种之一,该机型积累了丰富的事故数据资料,研究它的事故预测具有很强的代表意义,对于我国军机事故预测研究也具有较好的指导作用.

3.1 历史事故数据拟合统计

对事故进行预测,是建立在已经掌握了一定历史数据的基础之上的.假定现在掌握了该型飞机装备部队后第1年到第10年的历史数据,由前文分析,首先应当对数据进行拟合.见图5.

得到拟合方程

得到拟合曲线后,对拟合效果进行了拟合优度检验,结果证明在显著性水平α=0.01的条件下,判定成立.因此我们认为选取二次拟合是合适的.对这10年的残差进行统计,如图6.

统计这10年中残差绝对值的最大值

3.2 先验信息确定

前文中,做出了先验信息满足均匀分布的假设.但在实际运算过程中,确定样本的分布信息仍需要积累一定的样本数据作为佐证.为此通过对F-15、“狂风”战机等大量军机事故数据的拟合统计,结合该机型的先验样本信息,认为对其先验分布的客观规律是满足2.2.2中的选取原则的.由式(8)知:

3.3 后验分布确定及预测

掌握了事故数据的先验信息,还需要一定的样本数据对信息进行更新.贝叶斯方法的优势在于对信息的样本容量要求不高,在样本数量较少的条件下依然能够进行预测.假定我们现在仅掌握了第11年、第12年两年的事故数据,要对第13年的事故进行预测.由于新数据的添加,首先需要对这12年的事故数据进行重新拟合,并统计残差,如图7.

得到新的拟合方程

与式(15)相比,拟合方程的参数确实发生了变化.统计第11、12年的残差s11、s12:

依据式(9),有

又因为:

因此,σ的后验分布为:

由此就确定了σ的后验分布密度函数形式.由式(14),平方损失下σ的点估计为:

计算得到σ的期望ˆσ为0.92作为σ的预测值,根据拟合方程,将年份13代入式(18)中,计算第13年事故拟合值为µ13=18.8.则根据“3σ”法则:

则第13年事故数的预测值

即我们有99.7%的把握认为该机型在第13年发生的事故数范围在[16,22]区间之间.查阅文献,该机型在第13年实际发生事故为20起,确实在预测区间之内.

进一步,要想在预测区间的基础上再缩小预测范围,依然可以依照“3σ”法则来缩小区间.

从图8可以知道,将s13的预测区间缩小到[−2σ,2σ]时,仍有95.5%的把握认为预测是有效的.在此基础上,若要将预测区间缩小至[−σ,σ],预测的准确度就会下降.本算例中,按照[−σ,σ]区间对s13进行预测时,得出x13区间为[18,20],尽管可靠性只有68.3%,但第13年的事故数仍在预测区间内.这还是建立在样本信息只有两年的基础之上的.贝叶斯预测模型的准确性可见一斑.

该算例证实了本文基于贝叶斯统计决策理论建立的预测方法切实可行.在对本模型进行实际应用时,甚至可以以最近一年的事故数据作为样本信息对历史数据进行更新,从而对下一年事故进行预测.

4 结论

本文基于贝叶斯统计决策理论,建立了一套针对单一机型的事故预测模型,并以某型战斗机事故数据为例,对模型的可行性进行了验证,取得了良好的效果.但是,我们也应看到,军机的事故数据整理是一项严谨的工作,囿于事故资料来源的单一,我们无法保证能够完全掌握外军事故资料的完整数据.本文重点探讨运用贝叶斯统计决策理论建立起事故预测模型,即便在有限的数据条件下,也能够佐证预测模型的可行性.

随着我军装备技术的发展,装备成本也在不断提高.高发的事故率会导致军方难以承受高昂的代价,极大地影响部队对新装备的使用.因此,根据以往成熟机型摸索军机事故规律,指导新机型的事故防范工作,具有广泛的现实意义.

如今,我军新一代主战战机正在大规模装备航空兵部队,通过摸索典型军机从交付部队开始到装备发展成熟全过程的事故发生规律,建立起一套行之有效的事故预测方法,对于我国的军机研制单位和使用部队具有充分的指导作用,推动军机事故研究工作不断发展,从而为我军的军机安全保障工作提供数据支撑和理论支持.

1 KRUSCHKE J K.Doing Bayesian data analysis(影印版)[M].北京:机械工业出版社,2015.

2 王哲.现代贝叶斯统计学理论与方法研究[J].价值工程,2012,31(1):267−277.

3 SCHMELTERA M L,ERWINB S O,WILCOCKC P R.Accounting for uncertainty in cumulative sediment transport using Bayesian statistics[J].Geomorphology,2012,175(15):1−13.

4 韦程东.贝叶斯统计分析及其应用[M].北京:科学出版社,2015.

5 HEYDARI S,MIRANDA-MORENO L F,LORD D,et al.Bayesian methodology to estimate and update safety performance functions under limited data conditions:A sensitivity analysis[J].Accident Analysis and Prevention,2014,64(3):41−51.

6 LLOYDLK,FORSTERJJ.Modellingtrendsinroadaccidentfrequency-Bayesian inference for rates with uncertain exposure[J].Computational Statistics and Data Analysis,2014,73(2):189−204.

7 LIN L,WANG Q,ADEL W S.A novel variable selection method based on frequent pattern tree for real-time traffic accident risk prediction[J].Transportation Rearch Part C:Emerging Technologies,2015,55(6):444−459

8 贾丰源,孙杰,孙剑.快速路交通流运行安全关键参数识别与评估[J].同济大学学报(自然科学版),2015,43(2):221−225.

9 FULLWOOD R R,HALL R E,MARTINEZ-GURIDI G,et al.Relating aviation service difficulty reports to accident data for safety trend prediction[J].Reliability Engineering＆System Safety,1998,60(10):83−87.

10王永刚,李辉.基于灰色神经网络的民航事故征候预测模型研究[J].中国安全科学学报,2012,22(3):10−15.

11甘旭升,端木京顺,卢永祥.灰色均生函数模型及其在航空装备事故预测中的应用[J].中国安全科学学报,2010,20(6):40−45.

12刘选民,李凡.国外现代军机飞行事故[M].北京:航空工业出版社,2011.

13师义民,许勇,周丙常.近代统计方法[M].北京:高等教育出版社,2011.

14 PREVIC F H.Spatial disorientation in aviation[M].Reston:American Instritute of Aeronautics and Astronautics,Inc.,2004.

15赵宇.可靠性数学分析[M].北京:国防工业出版社,2011.