基于多变量灰色模型的航空安全预测方法

The aviation safety prediction method based on multivariate gray model

摘要

由于航空事故具有致因机理复杂且致因因素较强的灰色特性，传统的灰色预测模型只适用于单变量预测且具有预测精度低的缺陷。基于此，提出一种以遗传算法优化的多变量灰色模型开展航空安全预测的方法。首先，以航空事故致因理论为基础，从SHEL模型的角度运用鱼骨图确定航空安全影响因素，并以相关系数矩阵可视化图形进一步筛选关键致因因素；其次，构建以人为因素、环境因素、设备设施因素、外来影响因素等为强输入指标的多变量灰色模型，并利用遗传算法对模型的待定参数r全局搜索最优解；最后，以2007-2016年中国民用航空器事故征候万次率和航空不安全事件统计为对象进行仿真实验，并利用GM（1，1）和MGM（1，n）两种灰色预测模型进行预测对比。结果表明：提出的方法与传统的灰色模型相比，在短时间的航空安全预测中平均预测误差约为1.6%，验证了该方法的有效性和较高的预测精度。

关键词

航空安全 / 航空运输 / 短时预测 / 遗传算法 / 多变量灰色模型 / 参数寻优

Key words

aviation safety / air transportation / short-time prediction / genetic algorithm / multivariate gray model / parameter optimization

航空运输在建设交通强国的过程中具有先导性和标志性的作用，是交通运输行业不可或缺的重要组成部分。新型冠状病毒感染使全球航空运输的发展受到一定的影响，但航空运输正在逐步复苏。航空运输业在实现“量”的快速发展的同时，更要加强筑牢“安全质量”的界限，提升险情预警的能力。航空安全预测作为航空安全风险管理的重要环节，可通过深度挖掘致因因素与航空事故之间的内在联系得出航空安全状况的发展规律，实现航空安全的智能管理、提前决策及预警管理。航空运输系统运行的复杂性、非线性程度高等特点使得航空安全状态变化趋势复杂，加大了建立预测模型、开展长期精确安全状态预测的难度。同时，航空事故的发生是由多种不确定性因素相互交叉共同影响所致，而这些影响因素不仅具有随机性、时变性、高维性等特点，各因素之间的关联性也很强，这就增加了预测建模的难度。目前航空安全预测的方法主要分为时间序列预测、机器学习预测、计量模型预测3种。

时间序列预测法强调历史样本数据的时间序列特性，并从中挖掘变量的特征与发展规律，从而对变量的未来变化进行预测。在时间序列预测方法的研究中，整合移动平均自回归（auto-regressive integrated moving average，ARIMA）模型^［1-3］、自回归滑动平均（auto-regressive moving average，ARMA）^［4-5］模型、移动平均模型^［6］等几种方法应用较多。同时还有部分学者将时间序列预测方法和支持向量机^［7］、广义神经网络^［8］等方法相结合开展航空安全预测研究。显然该方法在应用时是以时间因子为单变量进行预测，然而航空事故致因机理复杂，航空安全状态的发展变化受很多因素的影响，因而存在着预测误差的缺陷。

在航空安全预测的研究中机器学习预测的应用同样较为广泛，主要方法有神经网络^［9-11］、贝叶斯网络^［12-13］、支持向量机^［14-15］、随机森林^［16］、机器学习^［17］等。这些方法是通过将大量数据输入到搭建的基于知识规则的映射模型中进行自我学习，并不断进行数据自我优化和改进。然而，航空运输的贫信息性限制了该类方法的使用。在航空安全预测的研究中，机器学习预测应用得同样较为广泛，主要方法有神经网络^［9-11］、贝叶斯网络^［12-13］、支持向量机^［14-15］、随机森林^［16］、机器学习^［17］等。该方法是通过将大量数据输入到搭建的基于知识规则的映射模型中进行自我学习中，并不断进行数据自我优化和改进。然而，航空运输的贫信息性同样限制了该类方法的使用。

计量模型预测主要以数理统计分析为基础，以数学方程的形式刻画预测对象与相关影响因素的作用机理。通过观察主要影响因素的变化趋势，估计模型参数以实现对预测对象未来变化趋势的推测。文献［18］采用灰色关联法对飞机火灾事故进行统计分析，但对影响因素的考量比较单一。文献［19］构建了一个灰色区间预测模型用于民航安全评估，可获得民航安全综合指数短时间的预测区间，却忽略了致因因素之间的相关性分析，不便于开展航空安全的智能管理。文献［20］将灰色理论和时间序列分析相结合开展航空装备事故组合预测模型，将两种模型计算得到的预测值进行相加得到组合预测结果，但所建立的组合模型无法根据系统实际状况对自身参数及时更新和调整，增加了预测的不确定性。文献［21］使用灰色马尔科夫模型状态预测方法对进近着陆不安全事件进行预测，结果表明所建立的模型预测精度高于经典GM（1，1）模型，而所建立的模型状态转移只取决于当前状态，不考虑历史状态，这种无关联性假设会增加预测结果的不确定性。

综合以上分析，灰色预测模型所呈现的不确定性特性能够较好地反映航空事故/事故征候数据小样本、贫信息的特点及航空事故致因因素的灰色特性。但目前的航空安全灰色预测方法以单因素度量居多，缺乏对多变量信息的挖掘，使得预测的准确性和适用性受到限制。鉴于此，提出基于遗传算法优化的多变量MGM（1，n，GA（r））模型开展航空安全预测。其中，MGM（1，n）模型不仅能够从多个角度综合考虑各个因素的变化，反映系统整体状态，而且可以建立多变量之间的动态相关性，弥补单变量预测的局限性。为了提高预测模型的适应性，用遗传算法对多变量灰色预测模型中的参数进行优化求解，并以实例去验证其可行性与有效性。

1 基于MGM（1，n，r）的预测建模

1.1 MGM（1，n）模型

（1）累加生成数列

假设通过相关性分析得到有n个以m为周期的灰数据序列 \{ x_i^{(\mathrm{0})} ( k ) \}（k＝1，2，…，m，i＝1， 2，…，n）。定义 x_{\mathrm{1}}^{(\mathrm{0})} ( k )为航空事故序列； x_i^{(\mathrm{0})} ( k )（i=2，3，…，n）为致因因素序列，其相应的一次累加生成数列为 \{ x_i^{(\mathrm{1})} ( k ) \}，其中

x_i^{(\mathrm{1})} ( k ) = {\displaystyle \sum _{j=\mathrm{1}}^k} x_i^{\left(\mathrm{0}\right)}(j) ，i=1，2，…，n，k=1，2，…，m （1）

记

\left\{\begin{array}{l}{x}_{\mathrm{1}}^{(\mathrm{0})}=\left\{x_{\mathrm{1}}^{(\mathrm{0})}(\mathrm{1}),x_{\mathrm{1}}^{(\mathrm{0})}(\mathrm{2}),\cdots ,x_{\mathrm{1}}^{(\mathrm{0})}(k)\right\}\\ x_{\mathrm{2}}^{(\mathrm{0})}=\left\{x_{\mathrm{2}}^{(\mathrm{0})}(\mathrm{1}),x_{\mathrm{2}}^{(\mathrm{0})}(\mathrm{2}),\cdots ,x_{\mathrm{2}}^{(\mathrm{0})}(k)\right\}\\ \begin{array}{cccc}& & & \begin{array}{cc}& ⋮\end{array}\end{array}\\ x_n^{(\mathrm{0})}=\left\{x_n^{(\mathrm{0})}(\mathrm{1}),x_n^{(\mathrm{0})}(\mathrm{2}),\cdots ,x_n^{(\mathrm{0})}(k)\right\}\\ x_{\mathrm{1}}^{(\mathrm{1})}=\left\{x_{\mathrm{1}}^{(\mathrm{1})}(\mathrm{1}),x_{\mathrm{1}}^{(\mathrm{1})}(\mathrm{2}),\cdots ,x_{\mathrm{1}}^{(\mathrm{1})}(k)\right\}\\ x_{\mathrm{2}}^{(\mathrm{1})}=\left\{x_{\mathrm{2}}^{(\mathrm{1})}(\mathrm{1}),x_{\mathrm{2}}^{(\mathrm{1})}(\mathrm{2}),\cdots ,x_{\mathrm{2}}^{(\mathrm{1})}(k)\right\}\\ \begin{array}{cccc}\begin{array}{cc}& \end{array}& & & ⋮\end{array}\\ x_n^{(\mathrm{1})}=\left\{x_n^{(\mathrm{1})}(\mathrm{1}),x_n^{(\mathrm{1})}(\mathrm{2}),\cdots ,x_n^{(\mathrm{1})}(k)\right\}\end{array}\right.

（2）构建微分方程组

\left\{\begin{array}{l}\frac{\mathrm{d}{x}_{\mathrm{1}}^{(\mathrm{1})}}{\mathrm{d}t}=a_{\mathrm{11}}{x}_{\mathrm{1}}^{(\mathrm{1})}+a_{\mathrm{12}}{x}_{\mathrm{2}}^{(\mathrm{1})}+\cdots a_{\mathrm{1}n}{x}_n^{(\mathrm{1})}+b_{\mathrm{1}}\\ \frac{\mathrm{d}{x}_{\mathrm{2}}^{(\mathrm{1})}}{\mathrm{d}t}=a_{\mathrm{21}}{x}_{\mathrm{1}}^{(\mathrm{1})}+a_{\mathrm{22}}{x}_{\mathrm{2}}^{(\mathrm{1})}+\cdots a_{\mathrm{2}n}{x}_n^{(\mathrm{1})}+b_{\mathrm{2}}\\ \begin{array}{cccc}& & & \begin{array}{cc}& ⋮\end{array}\end{array}\\ \frac{\mathrm{d}{x}_n^{(\mathrm{1})}}{\mathrm{d}t}=a_{n\mathrm{1}}{x}_{\mathrm{1}}^{(\mathrm{1})}+a_{n\mathrm{2}}{x}_{\mathrm{2}}^{(\mathrm{1})}+\cdots a_{nn}{x}_n^{(\mathrm{1})}+b_n\end{array}\right. （2）

令

{\mathit{X}}^{(\mathrm{0})} ( k ) = ( x_{\mathrm{1}}^{(\mathrm{0})} ( k ) , x_{\mathrm{2}}^{(\mathrm{0})} ( k ) , \cdots , x_n^{(\mathrm{0})} {(k))}^{\mathrm{T}}

{\mathit{X}}^{(\mathrm{1})} ( k ) = ( x_{\mathrm{1}}^{(\mathrm{1})} ( k ) , x_{\mathrm{2}}^{(\mathrm{1})} ( k ) , \cdots , x_n^{(\mathrm{1})} {(k))}^{\mathrm{T}}

\mathit{A} = \left[\begin{array}{cccc}{a}_{\mathrm{11}}& a_{\mathrm{12}}& ...& a_{\mathrm{1}n}\\ a_{\mathrm{21}}& a_{\mathrm{22}}& ...& a_{\mathrm{2}n}\\ & & & \\ ⋮& ⋮& & ⋮\\ a_{n\mathrm{1}}& a_{n\mathrm{2}}& ...& a_{nn}\end{array}\right] , \mathit{B} = {\left[\begin{array}{cccc}{b}_{\mathrm{1}}& b_{\mathrm{2}}& \cdots & b_n\end{array}\right]}^{\mathrm{T}}

则动态微分方程组可表示为

\frac{\mathrm{d}\mathit{X}}{\mathrm{d}t} = \mathit{A} {\mathit{X}}^{(\mathrm{1})} + \mathit{B} （3）

（3）建立时间响应函数

由积分生成变换原理得到连续的时间响应函数为

{\mathit{X}}^{(\mathrm{1})} ( t ) = e^{At} {\mathit{X}}^{(\mathrm{1})} ( \mathrm{0} ) + {\mathit{A}}^{-\mathrm{1}} ( e^{\mathit{A}\mathit{t}} - I ) \mathit{B} （4）

式中： e^{\mathit{A}k} = \mathit{I} + \mathit{A} t + \frac{{\mathit{A}}^{\mathrm{2}}}{\mathrm{2}!} t^{\mathrm{2}} + . . . = \mathit{I} + {\displaystyle \sum _{k=\mathrm{1}}^{\mathrm{\infty }}} \frac{{\mathit{A}}^n}{k!} t^n；

式中： I 为单位矩阵。

（4）辨识参数向量

为求得参数 A 和 B，将式（3）离散化后得到

x_i^{(\mathrm{1})} ( k + \mathrm{1} ) = {\displaystyle \sum _{j=\mathrm{1}}^n} \frac{a_{ij}}{\mathrm{2}} ( x_j^{(\mathrm{1})} ( k + \mathrm{1} ) + x_j^{(\mathrm{1})} ( k ) ) + b_i （5）

记 {\mathit{a}}_i = {\left[a_{i\mathrm{1}},a_{i\mathrm{2}},\cdots ,a_{in},b_i\right]}^{\mathrm{T}}，则由最小二乘法得到 a _i 的辨识值 {\widehat{\mathit{a}}}_i，即

{\widehat{\mathit{a}}}_i = {\left[{\widehat{a}}_{i\mathrm{1}},{\widehat{a}}_{i\mathrm{2}},\cdots ,{\widehat{a}}_{in},{\widehat{b}}_i\right]}^{\mathrm{T}} = ( {\mathit{L}}^{\mathrm{T}} {\mathit{L})}^{-\mathrm{1}} {\mathit{L}}^{\mathrm{T}} {\mathit{Y}}_i，i=1，2，…，n （6）

式中： {\mathit{Y}}_i = {\left[x_i^{(\mathrm{0})}(\mathrm{2}),x_i^{(\mathrm{0})}(\mathrm{3}),\cdots ,x_i^{(\mathrm{0})}(m)\right]}^{\mathrm{T}}

\mathit{L} = \left[\begin{array}{ccccc}\frac{\mathrm{1}}{\mathrm{2}}(x_{\mathrm{1}}^{(\mathrm{1})}(\mathrm{2})+x_{\mathrm{1}}^{(\mathrm{1})}(\mathrm{1}))& \frac{\mathrm{1}}{\mathrm{2}}(x_{\mathrm{2}}^{(\mathrm{1})}(\mathrm{2})+x_{\mathrm{2}}^{(\mathrm{1})}(\mathrm{1}))& \dots & \frac{\mathrm{1}}{\mathrm{2}}(x_n^{(\mathrm{1})}(\mathrm{2})+x_n^{(\mathrm{1})}(\mathrm{1}))& \mathrm{1}\\ \frac{\mathrm{1}}{\mathrm{2}}(x_{\mathrm{1}}^{(\mathrm{1})}(\mathrm{3})+x_{\mathrm{1}}^{(\mathrm{1})}(\mathrm{2}))& \frac{\mathrm{1}}{\mathrm{2}}(x_{\mathrm{2}}^{(\mathrm{1})}(\mathrm{3})+x_{\mathrm{2}}^{(\mathrm{1})}(\mathrm{2}))& \dots & \frac{\mathrm{1}}{\mathrm{2}}(x_n^{(\mathrm{1})}(\mathrm{3})+x_n^{(\mathrm{1})}(\mathrm{2}))& \mathrm{1}\\ ⋮& ⋮& & ⋮& \\ \frac{\mathrm{1}}{\mathrm{2}}(x_{\mathrm{1}}^{(\mathrm{1})}(m)+x_{\mathrm{1}}^{(\mathrm{1})}(m-\mathrm{1}))& \frac{\mathrm{1}}{\mathrm{2}}(x_{\mathrm{2}}^{(\mathrm{1})}(m)+x_{\mathrm{2}}^{(\mathrm{1})}(m-\mathrm{1}))& \dots & \frac{\mathrm{1}}{\mathrm{2}}(x_n^{(\mathrm{1})}(m)+x_n^{(\mathrm{1})}(m-\mathrm{1}))& \mathrm{1}\end{array}\right]

（5）建立预测模型

由参数估计离散化时间响应函数为

{\widehat{\mathit{X}}}^{(\mathrm{1})} = e^{\widehat{\mathit{A}}k} {\mathit{X}}^{(\mathrm{1})} ( \mathrm{1} ) + {\widehat{\mathit{A}}}^{-\mathrm{1}} ( e^{{}^{\widehat{\mathit{A}}k}} - I ) \widehat{\mathit{B}}，k=1， 2，…，n （7）

{\mathit{X}}^{(\mathrm{0})} ( \mathrm{1} ) = {\widehat{\mathit{X}}}^{(\mathrm{1})} ( \mathrm{1} ) = {\mathit{X}}^{(\mathrm{0})} ( \mathrm{1} ) （8）

{\widehat{\mathit{X}}}^{(\mathrm{0})} ( k ) = {\mathit{X}}^{(\mathrm{1})} ( k ) - {\mathit{X}}^{(\mathrm{1})} ( k - \mathrm{1} )， k=2，3，…，n （9）

其中

\widehat{\mathit{A}} = \left[\begin{array}{cccc}{\widehat{a}}_{\mathrm{11}}& {\widehat{a}}_{\mathrm{12}}& ...& {\widehat{a}}_{\mathrm{1}n}\\ {\widehat{a}}_{\mathrm{21}}& {\widehat{a}}_{\mathrm{22}}& ...& {\widehat{a}}_{\mathrm{2}n}\\ & & & \\ ⋮& ⋮& & ⋮\\ {\widehat{a}}_{n\mathrm{1}}& {\widehat{a}}_{n\mathrm{2}}& ...& {\widehat{a}}_{nn}\end{array}\right]

\widehat{\mathit{B}} = {\left[\begin{array}{cccc}{\widehat{b}}_{\mathrm{1}}& {\widehat{b}}_{\mathrm{2}}& \dots & {\widehat{b}}_n\end{array}\right]}^{\mathrm{T}}

1.2 MGM（1，n，r）模型的构建

（1）对动态微分方程组进行差分，将公式 （3）进行向前差分得到

\frac{{\mathit{X}}_{t\mathrm{1}}^{(\mathrm{1})}-{\mathit{X}}_{t\mathrm{2}}^{(\mathrm{1})}}{t_{\mathrm{1}}-t_{\mathrm{2}}} - \mathit{A} {\mathit{X}}_{t\mathrm{2}}^{(\mathrm{1})} = \mathit{B} （10）

由于时间序列中t ₁-t ₂=1，故有

{\mathit{X}}_{t\mathrm{1}}^{(\mathrm{1})} - {\mathit{X}}_{t\mathrm{2}}^{(\mathrm{1})} - \mathit{A} {\mathit{X}}_{t\mathrm{2}}^{(\mathrm{1})} = \mathit{B} （11）

向后差分得

{\mathit{X}}_t^{(\mathrm{1})} - {\mathit{X}}_{t-\mathrm{1}}^{(\mathrm{1})} - \mathit{A} {\mathit{X}}_t^{(\mathrm{1})} = \mathit{B}或{\mathit{X}}_{t+\mathrm{1}}^{(\mathrm{1})} - {\mathit{X}}_t^{(\mathrm{1})} - \mathit{A} {\mathit{X}}_{t+\mathrm{1}}^{(\mathrm{1})} = \mathit{B} （12）

（2）建立组合差分格式

序列的差异性导致差分格式的不同，基于此，建立两种差分方式的组合格式，即

( {\mathit{X}}_{t\mathrm{1}}^{(\mathrm{1})} - {\mathit{X}}_{t\mathrm{2}}^{(\mathrm{1})} - \mathit{A} {\mathit{X}}_{t\mathrm{2}}^{(\mathrm{1})} ) \times r + \left[({\mathit{X}}_{t+\mathrm{1}}^{(\mathrm{1})}-{\mathit{X}}_t^{(\mathrm{1})}-\mathit{A}{\mathit{X}}_{t+\mathrm{1}}^{(\mathrm{1})})\times (\mathrm{1}-r)\right] = \mathit{B} （13）

（3）更新辨识值 \widehat{\dot{\mathit{A}}}和 \widehat{\dot{\mathit{B}}}

很显然，当r=0.5时，即为GM（1，1）模型。对于给定的参数r对应的矩阵 L ₁和 \widehat{\dot{\mathit{a}}} _i 分别为

{\mathit{L}}_{\mathrm{1}} = \left[\begin{array}{ccccc}{r}_{\mathrm{0}}{x}_{\mathrm{1}}^{(\mathrm{1})}(\mathrm{2})+(\mathrm{1}-r_{\mathrm{0}})x_{\mathrm{1}}^{(\mathrm{1})}(\mathrm{1})& r_{\mathrm{0}}{x}_{\mathrm{2}}^{(\mathrm{1})}(\mathrm{2})+(\mathrm{1}-r_{\mathrm{0}})x_{\mathrm{2}}^{(\mathrm{1})}(\mathrm{1}))& \dots & r_{\mathrm{0}}{x}_n^{(\mathrm{1})}(\mathrm{2})+(\mathrm{1}-r_{\mathrm{0}})x_n^{(\mathrm{1})}(\mathrm{1}))& \begin{array}{l}\mathrm{1}\\ \end{array}\\ r_{\mathrm{0}}{x}_{\mathrm{1}}^{(\mathrm{1})}(\mathrm{3})+(\mathrm{1}-r_{\mathrm{0}})x_{\mathrm{1}}^{(\mathrm{1})}(\mathrm{2})& r_{\mathrm{0}}{x}_{\mathrm{2}}^{(\mathrm{1})}(\mathrm{3})+(\mathrm{1}-r_{\mathrm{0}})x_{\mathrm{2}}^{(\mathrm{1})}(\mathrm{2}))& \dots & r_{\mathrm{0}}{x}_n^{(\mathrm{1})}(\mathrm{3})+(\mathrm{1}-r_{\mathrm{0}})x_n^{(\mathrm{1})}(\mathrm{2}))& \begin{array}{l}\mathrm{1}\\ \end{array}\\ ⋮& ⋮& & ⋮& \\ r_{\mathrm{0}}{x}_{\mathrm{1}}^{(\mathrm{1})}(m)+(\mathrm{1}-r_{\mathrm{0}})x_{\mathrm{1}}^{(\mathrm{1})}(m-\mathrm{1})& r_{\mathrm{0}}{x}_{\mathrm{2}}^{(\mathrm{1})}(m)+(\mathrm{1}-r_{\mathrm{0}})x_{\mathrm{2}}^{(\mathrm{1})}(m-\mathrm{1})& \dots & r_{\mathrm{0}}{x}_n^{(\mathrm{1})}(m)+(\mathrm{1}-r_{\mathrm{0}})x_n^{(\mathrm{1})}(m-\mathrm{1})& \mathrm{1}\end{array}\right] （14）

{\widehat{\dot{\mathit{a}}}}_i = \left[\begin{array}{c}{\widehat{\dot{a}}}_{i\mathrm{1}}\\ {\widehat{\dot{a}}}_{i\mathrm{2}}\\ ⋮\\ {\widehat{\dot{a}}}_{in}\\ {\widehat{\dot{b}}}_i\end{array}\right] = ( {{\mathit{L}}_{\mathrm{1}}}^{\mathrm{T}} {\mathit{L}}_{\mathrm{1}} {)}^{-\mathrm{1}} {{\mathit{L}}_{\mathrm{1}}}^{\mathrm{T}} {\mathit{Y}}_i （15）

因此，从上述分析可知，在参数r确定的基础上，即可通过式（14）和（15）确定辨析值 \widehat{\dot{\mathit{A}}}和 \widehat{\dot{\mathit{B}}}。

（4）建立预测模型

通过以上推论分析，得到MGM（1，n，r）模型中的辨识值 \widehat{\dot{\mathit{A}}}、 \widehat{\dot{\mathit{B}}}及参数r，然后由时间响应函数得出模型的预测值。

1.3 基于遗传算法的MGM（1，n，r）模型参数优化

（1）遗传算法

遗传算法是通过模拟自然选择、复制、交叉、变异等过程，得到待解问题优质解的随机全局搜索优化方法 ^［22］。详细算法流程如图1所示。

（2）参数优化

遗传算法具有内在的隐并行性、自适应性及强鲁棒性等优点，已广泛应用于航迹优化、信号处理和人工智能等领域。基于此，本文将遗传算法取代灰色系统理论中的最小二乘法进行求解，以全局并行搜索的方式获得MGM （1，n，r）模型中参数r的最优解。当用遗传算法求出MGM（1，n，r）模型中参数r的最优值时，代入到式（14）和（15）中，获得辨析值 \widehat{\dot{\mathit{A}}}和 \widehat{\dot{\mathit{B}}}的数值，然后在建立的MGM （1，n，r）航空安全预测模型中开展分析计算，如式（16）~（18）所示。

{\widehat{\dot{\mathit{X}}}}^{(\mathrm{1})} = e^{\widehat{\dot{A}}k} {\mathit{X}}^{(\mathrm{1})} ( \mathrm{1} ) + {\widehat{\dot{\mathit{A}}}}^{-\mathrm{1}} ( e^{\widehat{\dot{A}}k} - I ) \widehat{\dot{\mathit{B}}} , k=1，2，…，n （16）

{\widehat{\mathit{X}}}^{(\mathrm{0})} ( k ) = {\widehat{\mathit{X}}}^{(\mathrm{1})} ( k ) - {\widehat{\mathit{X}}}^{(\mathrm{1})} ( k - \mathrm{1} )，k=2，3，…，n （17）

\mathrm{相对} \mathrm{误差} = \frac{\left|{\widehat{x}}^{(\mathrm{0})}-x_i^{(\mathrm{0})}\right|}{x_i^{(\mathrm{0})}} （18）

1.4 基于遗传MGM（1，n，GA（r））模型的航空安全预测流程

（1）致因因素确定

航空运输是一个复杂且庞大的系统，保持航空运输的持续安全性需要机场、飞行机务等多方面的共同协作。从航空安全特征的角度分析，可理解为航空运输的安全性受众多因素的影响，且影响作用机理复杂。由航空事故致因理论可知，影响因素又是诱发各类航空事故/事故征候的风险源。因此，可得出各影响因素与航空事故/事故征候之间满足一定的因果映射关系，用模型的形式可表示为

x ₁=G（x_i ），i=2，3，…n （19）

式中： x ₁表示实际发生的航空事故或事故征候； x _i 表示导致各类航空事故的致因因素。

在本文中，以航空事故致因理论为基础，从SHEL模型的角度将航空事故致因因素分为人为因素、环境因素、管理因素、设备设施因素、外来影响因素5个方面。通过查阅文献、调研分析及咨询专家意见，以鱼骨图的方法呈现对航空安全风险因素的分析结果，如图2所示。

（2）数据处理

为了消除各致因因素量纲不一致对数据分析带来的干扰，对各指标进行归一化处理。将各因素的数值按比例缩放在区间（0，1）内，如式（20）所示。

x^{\text{'}} = \frac{x-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}{x_{\mathrm{m}\mathrm{a}\mathrm{x}}-x_{\mathrm{m}\mathrm{i}\mathrm{n}}} （20）

式中： x_{\mathrm{m}\mathrm{a}\mathrm{x}}、 x_{\mathrm{m}\mathrm{i}\mathrm{n}}分别为样本数据中的最大值和最小值； x^{\text{'}}为归一化后的样本数据。

（3）相关性分析

以定量分析的方式确定各致因变量对航空安全水平的影响程度，从而在减少噪声干扰的同时排除弱相关变量，降低运算成本。同时，在利用多变量灰色模型开展航空安全预测时，也需要考虑各变量之间的关联度是否满足进行联合预测的要求。在文中采用皮尔逊、斯皮尔曼、肯德尔分别以定距、定序、定秩的尺度衡量各致因变量的相关性，从而实现线性相关与单调相关的兼容性，如式（21）~（23）所示。

{\gamma }_{xy} = \frac{{\displaystyle \sum _{i=\mathrm{1}}^n}(x_i-\overline{x})(y_i-\overline{y})}{\sqrt[]{{\displaystyle \sum _{i=\mathrm{1}}^n}(x_i{-\overline{x})}^{\mathrm{2}}(y_i{-\overline{y})}^{\mathrm{2}}}} （21）

\rho = \frac{\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}(R(x_i)-\overline{R(x)})(R(y_i)-\overline{R(y)})}{\sqrt[]{(\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}(R(x_i{)-\overline{R(x)})}^{\mathrm{2}}{(\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}(R(y_i)-\overline{R(y)})}^{\mathrm{2}})}} （22）

\tau = \frac{N_c-N_d}{n(n-\mathrm{1})/\mathrm{2}} （23）

式中： R ( x_i )、 R ( y_i )分别为两属性变量x_i 、y_i 的秩； \overline{R(x)}、 \overline{R(y)}为平均秩；N_c 为两属性度量下观测点排序一致的对数；N_d 为排序不一致的对数。

（4）航空安全预测建模

构建的MGM（1，n，GA（r））航空安全预测模型流程如图3所示。

（5）模型精度评价

为了验证所建立预测模型的精度，利用GM（1，1）和MGM（1，n）两种灰色预测模型进行拟合预测，对比三者之间的预测值及误差绝对值。

2 案例分析

（1）数据分析与处理

文中数据选自2007-2016年中国民用航空器事故征候统计，包括航空不安全事件统计和航空事故征候万次率统计^［23］。其中，航空不安全事件主要指诱发各类航空器事故征候的风险源，以航空安全风险因素分析为基础，主要从5个方面建立数据清单，并对数据进行归一化处理，如表1所示。

（2）相关性分析

利用式（21）～（23）定量分析航空不安全事件与航空事故征候之间的相关性。其中，x ₁为航空事故征候万次率，x_i （i=2，3，4，5，6）表示外来影响因素、人为因素、环境因素、设备设施因素、管理因素5类航空不安全事件，具体计算结果如表2所示。通过R语言的corrplot程序包画出相关系数二维矩阵图直观地显示计算结果。在二维矩阵图中，不同指标之间的相关性可以通过矩阵单元格灰度的深浅来表示。颜色越深表示关联性越强，反之，相关性越弱，如图4~6所示。

通过3种相关性分析可知，在斯皮尔曼、皮尔逊、肯德尔3种指数上，人为因素对航空事故征候的关联程度相对于其他因素是最大的，分别为0.854 8、0.900 5、0.750 0。而在斯皮尔曼和皮尔逊指数上，外来影响因素和设备设施因素对事故征候的关联程度均超过了0.75。据此初步确定人为因素、外来影响因素、设备设施因素对事故征候的影响程度较大，存在较强的正相关性，将作为航空安全预测模型的强输入指标。管理因素和环境因素在3种相关性的分析中，关联度均较小于0.5。尤其是管理因素，在3种指数的分析上，关联度出现了最小值0.271，因此，可以判断管理因素与航空事故征候之间是弱相关。此外，考虑到环境因素在斯皮尔曼、皮尔逊两种指数上，相关度均超过了0.4，说明尽管Kendall指数较低，其与事故征候的线性独立性、非线性单调相关性较强。所以经过以上综合分析，判定人为因素、设备设施因素、外来影响因素、环境因素与航空事故征候之间存在强相关性，将作为航空安全预测模型的强输入指标。

（3）模型参数优化

作为解决结构优化设计问题的一种高级算法，遗传算法各操作参数的选择（如种群大小N、适应度函数、交叉概率P_c 、变异概率P_m ）直接影响优化结果。本文针对各操作参数的选择不再进行详细描述，根据问题的固有知识及相关经验，在利用遗传算法寻优时采用二进制编码，群体数目12，编码长度22、交叉概率0.9、变异概率0.05、最大代数200。观察进化过程，当迭代次数为60时，参数r已达到收敛，最终寻优参数为r=0.643 32.

（4）结果分析与比较

将MGM（1，5）和GM（1，1）两种预测模型进行对比分析，得到3种模型的预测结果如表3所示。

由表3和图7可知，灰色算法GM（1，1）在航空安全预测时误差最大，平均约为7%。而在MGM（1，5）算法中除了航空事故征候自身发展规律外，充分考虑了人为因素、外来影响因素等对航空事故征候的影响，使得预测精度明显高于单变量度量的GM（1，1）预测模型。针对本文提出的基于遗传算法优化的多变量灰色MGM（1，5，GA（r））预测模型预测精度最高，平均预测误差约为1.6%。可以看出，灰色理论中的模糊特性能够较好地体现航空安全预测的不确定性因素。同时，若提高预测精度，需要对预测模型中所涉及不确定性因素进行耦合分析处理。

3 结论

为了实现航空事故发生前的精确预测，以航空事故致因理论为基础考虑人为因素、外来影响因素、设备设施因素、环境因素等对航空事故的影响机理。提出多变量灰色模型来反映各致因因素之间的相互影响关系，同时利用遗传算法的全局并行搜索能力来优化模型参数，建立基于遗传算法优化的多变量灰色MGM（1，n，GA（r））航空安全预测模型。以中国民航2007-2016年中国民航运行相关数据为例进行验证并与MGM（1，n）、GM（1，1）两种算法进行比较。

通过验证发现，灰色模型能够较好地解决航空安全影响因素的随机性以及致因因素之间的关联性问题，从而能够反映航空安全状况的整体趋势。同时，经过对比分析得到，文中所建立的多变量灰色遗传MGM（1，n，GA（r））航空安全预测模型相较于MGM（1，n）、GM（1，1）两种预测模型，有效降低了相对误差，平均预测误差约为1.6%。因此，MGM（1，5，GA（r））预测模型适用性更强、预测精度较高，更适合作为未来航空安全预警管理的新方法。

参考文献

[1] 程明，梁文娟.MLR和ARIMA模型在民航安全业绩预测中的应用［J］.中国安全科学学报，2016，26（2）：25-30.

[2] Li C， Chiang T W.Complex neurofuzzy ARIMA forecasting：a new approach using complex fuzzy sets［J］.IEEE Transactions on Fuzzy Systems，2013，21（3）：567-584.

[3] 甘旭升，端木京顺，高建国，基于ARIMA模型的航空装备事故时序预测［J］.中国安全科学学报，2012，22（3）：97-102.

[4] 杜毅.民用航空事故征候万时率的回归分析［J］.中国民航飞行学院学报，2010，21（1）：42-44.

[5] Wei P F， Lu Z Z， SONG J W. Variable importance analysis： a comprehensive review［J］. Relia-bility Engineering and System Safety，2015，142（3）：399-432.

[6] Shyur H J.A quantitative model for aviation safety risk assessment［J］.Computers and Industrial Engineering，2008，54（1）：34-44.

[7] 戴蓉，黄成.飞机飞行事故率预测建模与仿真研究［J］.计算机仿真，2011，28（7）：120-123.

[8] 郭娜.民航不安全事件预测研究［D］.天津：中国民航大学，2019.

[9] 曾航，张红梅，任博，基于改进LSTM模型的航空安全预测方法研究［J］.系统工程与电子技术，2022，44（2）：569-576.

[10] 于洪霞，李兴.一种基于MLP神经网络的大额损失飞行事故预测模型［J］.上海电力学院学报，2016，32（5）：504-506.

[11] Santhosh T V， Gopika V， Ghosh A K，et al.An approach for reliability prediction of instrumentation & control cables by artificial neural networks and Weibull theory for probabilistic safety assessment of NPPs［J］.Reliability Engineering & System Sa-fety，2018，170（4）：31-44.

[12] Arnaldo V R M， Comendador V F， Perez S L，et al.Prediction of aircraft safety incidents using Bayesian inference and hierarchical structures［J］.Safety Science，2018，104（6）：216-230.

[13] Bauranov A， Rakas J.Bayesian network model of aviation safety：impact of new communication technologies on mid-air collisions［J］.Reliability Engineering & System Safety，2023，243（3）：1-15.

[14] 李大伟，徐浩军，胡良谋，基于支持向量机的飞行事故率预测模型［J］.数学的实践与认识，2009，39（8）：124-128.

[15] 梁文娟，李雪艳.基于ARIMA、LS-SVM和BP神经网络组合模型的航空运输飞行事故征候预测［J］.安全与环境工程，2018，25（1）：130-136.

[16] 任博，岳珠峰，司勇，基于随机森林的航空安全因果预测新方法［J］.系统工程与电子技术，2023，45（3）：762-768.

[17] Zhang X G， Mahadevan S.Ensemble machine learning models for aviation incident risk prediction［J］.Decision Support Systems，2019，116（11）：48-63.

[18] 凤四海，李枣，贺元骅.基于灰色关联法的飞机火灾事故统计分析与启示［J］.安全与环境工程，2017，24（3）：138-143，149.

[19] 王衍洋，曹义华.民航安全指数结果分析与预测［J］.北京航空航天大学学报，2011，37（10）：1223-1227.

[20] 甘旭升，端木京顺，王青.航空装备事故的灰色时序组合预测模型［J］.中国安全科学学报，2012，22（4）：32-37.

[21] 陈农田，张宏钰，曹玉宽.基于灰色马尔科夫理论的进近着陆不安全事件预测［J］.数学的实践与认识，2020，50（4）：306-314.

[22] 李延梅.一种改进的遗传算法及应用［D］.广州：华南理工大学，2012.

[23] 中国民用航空局航空安全办公室.中国民航航空安全报告［R］.北京：中国民用航空局，2016.

基金资助

山东省自然科学基金(ZR2022QE284)

山东省重点研发计划项目(2023RKY02016)

滨州学院科研基金项目(BZXYLG202204)