一种面向多重失效模式的民用飞机运行风险评估方法

本发明涉及航空器安全技术领域，具体为一种面向多重失效模式的民用飞机运行风险评估方法。

背景技术：

民用飞机因其高效性成为越来越多人出行的首选交通工具。对于民航运输业来说，安全性是第一位的，保证航空器在预期运行环境和使用限制下持续安全地飞行对所有在役运行的民用飞机十分重要。

现有的航空安全分析方法主要针对单个航空器，已形成包含定量与定性分析的较为完善的安全管理体制。运输类飞机风险评估方法(transportairplaneriskassessmentmethodology,taram)是美国联邦航空局(faa)指令8110.107a监控安全-分析数据(monitorsafety-analyzedata,masd)中的定量风险分析方法，该方法提供了针对机队和单机不同运行时期、不同风险值的量化方法，用于执行msad指令所必要的运输类机队持续运营安全的定量风险分析问题。针对机队运行的不同时期，方法分别提供了三种耗损失效风险量值：总体未纠正机队风险、90天机队风险和控制程序机队风险，其中总体未纠正机队风险(rt)用于对机队是否存在不安全状况的风险进行长期预测以及指导安全决策，可根据公式计算得出rt＝da×nd×cp×ir。其中，ir和nd通常由机队运营经验给出；针对cp，需要建立由初因事件导致不安全结果的因果链进行分析，并结合历史数据得到各失效事件发生概率，通过将不安全结果所有状况的条件概率累加得到最终概率值；对初因事件失效模式未知的情况，可利用威布尔分析工具或可靠性分析软件得到失效分布参数，根据公式计算da值，现有的民用飞机运行风险评估方式为通过运用taram方法对常失效率或单一耗损失效的情况进行研究。

但本申请发明人在实现本申请实施例中发明技术方案的过程中，发现上述技术至少存在如下技术问题：

1、民机风险分析方法仅考虑了单一失效对航空器运行安全的影响，然而随着民机结构复杂化和服役时间增加，产生的耗损失效通常呈现多重失效模式的特点，利用该方法无法准确评估机队运行风险；

2、当机队出现多重失效模式时，相应的失效分布函数相对复杂，利用解析法需要多次计算每架飞机当前时间与退役时间的失效累积分布函数值，计算过程复杂，效率较低；

3、在研究由初因事件到不安全结果的失效传播关系时，通过因果链分析逐条计算每条支链的风险值，并累加得到最终cp值，当因果链支链过多时计算过程异常繁琐。

基于此，本发明设计了一种面向多重失效模式的民用飞机运行风险评估方法，以解决上述问题。

技术实现要素：

本发明的目的在于提供一种面向多重失效模式的民用飞机运行风险评估方法，以解决上述背景技术中提出的现有的背景技术提及的技术问题。

为实现上述目的，本发明提供如下技术方案：一种面向多重失效模式的民用飞机运行风险评估方法，包括：

针对民用飞机结构件具有多重耗损失效的特点，利用机队运行中的失效数据，建立基于混合威布尔分布的飞机结构件失效模型，并基于pso优化的em算法进行模型参数估计；

执行蒙特卡洛仿真算法，根据失效模型模拟机队中出现多重耗损失效的缺陷飞机数量da；

构建贝叶斯网络对导致不安全结果的初因事件进行失效传播分析，计算初因事件发生条件下的不安全结果发生概率cp；

根据预定的损伤率ir和未检出率nd，计算得到机队风险水平rt＝da×nd×cp×ir。

进一步的，所述飞机结构件失效模型的构建包括：

利用失效数据t＝[ti]，构建完全数据的似然函数：

其中，n为样本数据总量，k为威布尔子分布数量，z＝[zik]，zik为隐变量，并且当第i个观测数据来自第k个威布尔子分布，则zik＝1，否则zik＝0；θ为分布的参数集合，表示为θ＝(π1，...，πk；θ1，...，θk)，πk为第k个子分布加权因子，θk＝(αk，βk)表示第k个子分布尺度参数αk和形状参数βk；f(·|θk)为第k个威布尔子分布的概率密度函数；

执行em算法，通过给定失效数据t和当前参数估计值θ^(j)求z的条件概率p(z|t，θ^(j))，获取该条件概率的似然函数期望值q：

将该q函数作为目标函数，建立混合威布尔分布参数估计的优化模型如下：

maxq(θ，θ^(j))

进一步的，执行em算法过程中还包括：

确定参数θ^(j)第j次迭代值对应的q函数，并对q(θ，θ^(j))进行极大化处理得到第j+1次迭代结果θ^(j+1)，重复以上操作，直至收敛。

进一步的，执行蒙特卡洛仿真算法，根据失效模型模拟机队中出现多重耗损失效的缺陷飞机数量da包括：

计算加权因子的累积概率其中q0＝0，k为威布尔子分布个数；

对每架飞机生成服从[0，1]均匀分布的随机数r，确定r属于的概率区间[qk-1，qk]，得到飞机服从第k个威布尔子分布失效模式；

执行逆变换模拟生成飞机故障时间其中，αk、βk分别为第k个子分布尺度参数、第k个子分布形状参数；

统计机队中所有飞机在失效预测期间内的tm数量，并将tm多次模拟后的平均值作为缺陷飞机数量da值。

进一步的，构建贝叶斯网络对导致不安全结果的初因事件进行失效传播分析包括：

分析飞机结构件由初因事件发展至不安全结果的失效传播链；

将每一级中间失效对应于bn中各节点事件；

选用特定符描述节点之间的因果关系，构建直观描述不安全结果发展的有向无环图。

进一步的，计算初因事件发生条件下的不安全结果发生概率cp包括：

通过变量消元法将bn联合分布分解为因子函数的乘积形式；

获取经消元顺序将每个因子函数中相应变量消除的只含查询变量的因子乘积式；

将初因事件设为证据变量，并将不安全结果设为查询变量，通过变量消元更新各节点分布并求出初因事件发生条件下的不安全结果cp值。

进一步的，通过pso优化em算法包括：

将pso中每一个粒子的d维位置信息对应于参数θ；

在算法初始化时，根据各参数的初始值确定粒子位置与速度的更新界限，并随机生成数量为popsize的粒子种群；

将q函数作为pso的适应度函数，通过比较每代粒子的适应度值更新粒子，并保留适应度值最优的粒子；

更新粒子至最大寻优次数，输出全局最优粒子位置作为本轮em迭代的极大化结果。

本发明实施例中提供的一个或多个技术方案，至少具有如下技术效果或优点：

1、针对多重失效模式开展运行风险评估，克服了现有风险评估方法只能应用于单一失效模式的缺陷；

2、通过利用bn分析初因事件失效传播关系，避免了使用因果链分析多种中间失效同时发生时的繁琐过程；

3、基于pso优化的em算法克服了常用的搜索和逼近算法难以解决混合分布参数估计问题的缺点。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明评估方法的流程框图；

图2为本发明基于pso优化em算法的流程框图；

图3为本发明实施机队初因事件失效概率密度图；

图4为本发明实施机队当前运行状况图；

图5为本发明基于蒙特卡洛模拟机队da流程框图；

图6为本发明实施机队da模拟结果图；

图7为本发明实施机队机翼裂纹化简前后bn图；

图8为本发明实施机队机翼翼肋寿命数据表图；

图9为本发明实施机队bn节点事件信息表图；

图10为本发明实施机队不同结果损伤率表图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

请参阅图1-10，本发明实施例提供一种技术方案：一种面向多重失效模式的民用飞机运行风险评估方法，包括：

针对民用飞机结构件具有多重耗损失效的特点，利用机队运行中的失效数据，建立基于混合威布尔分布的飞机结构件失效模型，并基于pso优化的em算法进行模型参数估计；

执行蒙特卡洛仿真算法，根据失效模型模拟机队中出现多重耗损失效的缺陷飞机数量da；

构建贝叶斯网络对导致不安全结果的初因事件进行失效传播分析，计算初因事件发生条件下的不安全结果发生概率cp；

根据预定的损伤率ir和未检出率nd，计算得到机队风险水平rt＝da×nd×cp×ir。

通过上述步骤，不难看出，在本申请的面向多重失效模式的民用飞机运行风险评估方法中，利用机队运行中的失效数据建立基于混合威布尔分布的飞机结构件多重失效模型，并基于pso优化em算法，可以实现对多重失效模型进行参数估计，而且基于pso优化em算法，可以有效克服常用的搜索和逼近算法难以解决混合分布参数估计问题的缺点。

威布尔分布是可靠性模型中应用最为广泛的一种寿命分布，主要优点为数据拟合性好、适用于小样本抽样以及对各类数据的适应性极强，利用混合威布尔分布拟合失效数据的寿命分布，可以实现对完全数据的似然函数的构建，所述飞机结构件失效模型的构建包括：

利用失效数据t＝[ti]，构建完全数据的似然函数：

其中，n为样本数据总量，k为威布尔子分布数量，z＝[zik]，zik为隐变量，并且当第i个观测数据来自第k个威布尔子分布，则zik＝1，否则zik＝0；θ为分布的参数集合，表示为θ＝(π1，...，πk；θ1，...，θk)，πk为第k个子分布加权因子，θk＝(αk，βk)表示第k个子分布尺度参数αk和形状参数βk；f(·|θk)为第k个威布尔子分布的概率密度函数；

执行em算法，通过给定失效数据t和当前参数估计值θ^(j)求z的条件概率p(z|t，θ^(j))，获取该条件概率的似然函数期望值q：

将该q函数作为目标函数，建立混合威布尔分布参数估计的优化模型如下：

maxq(θ，θ^(j))

由于混合威布尔分布参数的q函数通常为复杂的非线性方程，无法直接获得其解析值，利用基于pso优化的em算法实现混合威布尔分布参数估计的过程时，执行em算法过程中还包括：

确定参数θ^(j)第j次迭代值对应的q函数，并对q(θ，θ^(j))进行极大化处理得到第j+1次迭代结果θ^(j+1)，重复以上操作，直至收敛；

需要说明的是，em是一种对含有隐变量的概率分布参数的极大似然估计方法，主要包含：期望步(e-step)和极大化步(m-step)。在e步中，确定第j次迭代参数结果θ^(j)对应的q函数值，并在m步中对q(θ，θ^(j))进行极大化得到第j+1次迭代结果θ^(j+1)，重复以上两步直至达到收敛条件。

进一步的，执行蒙特卡洛仿真算法通过失效模型模拟机队中出现多重耗损失效的缺陷飞机数量da包括：

计算加权因子的累积概率其中q0＝0，k为威布尔子分布数；

对每架飞机生成服从[0，1]均匀分布的随机数r，确定r属于的概率区间[qk-1，qk]，得到飞机服从第k个威布尔子分布失效模式；

执行逆变换模拟生成飞机故障时间其中αk、βk分别为第k个子分布尺度参数、第k个子分布形状参数；

统计机队中所有飞机在失效预测期间内的tm数量，并将tm多次模拟后的平均值作为缺陷飞机数量da值；

在执行执行蒙特卡洛仿真算法通过失效模型模拟机队中出现多重耗损失效的缺陷飞机数量da时，1、输入机队飞机总量m以及每架飞机当前运行时间t0(mm)和退役时间tr，确定最大循环次数itermax，初始化da，开始循环；2、输入多重失效分布模型参数估计结果，针对权重不同的威布尔子分布，计算每个子分布加权因子的累积概率其中q0＝0，随机生成服从[0，1]均匀分布的随机数r，根据r属于的概率区间[qk-1，qk]，确定第m架飞机服从第k个威布尔子分布失效模式，模拟生成该架飞机故障时间为3、判断该架飞机tm是否在失效预测期间[t0(m)，tr]之内，若是，则da数量加1，否则返回到步骤2生成下一架飞机故障时间；4、遍历机队中每架飞机完成一次模拟，记录该次模拟得到的da值；重复2～4至最大循环次数itermax，输出全部da的平均值。

针对导致不安全结果的初因事件开展失效传播分析的bn实施方式为，构建贝叶斯网络对导致不安全结果的初因事件进行失效传播分析包括：

分析飞机结构件和导致不安全结果致因的失效传播链；

将每一级中间失效对应于bn中各节点事件；

选用特定符描述节点之间的因果关系，构建直观描述不安全结果发展的有向无环图；

需要说明的是，通过上述内容，首先分析飞机结构件由初因事件导致最终不安全结果的失效传播关系，每一级中间失效对应bn中各节点事件，并用例如“边”表示节点之间的因果关系，构建直观描述不安全结果发展的有向无环图。

运用变量消元法对bn推理时，计算初因事件发生条件下的不安全结果发生概率cp包括：

通过变量消元法将bn联合分布分解为因子函数的乘积形式；

获取经消元顺序将每个因子函数中相应变量消除的只含查询变量的因子乘积式；

将初因事件设为证据变量，并将不安全结果设为查询变量，通过变量消元更新各节点分布并求出初因事件发生条件下的不安全结果cp值；

具体的为：将bn联合分布分解为因子函数的乘积形式，得到由因子函数组成的概率分布集合ζ；设定初因事件为证据变量、不安全结果为查询变量，为证据变量设置观测值并更新含有证据的因子式；对除证据和查询变量外的剩余变量排序，在更新后的ζ中按顺序依次对所含变量进行边缘化消去；将消元后的ζ中的因子式相乘得到只含查询变量的概率分布，再对其进行归一化得到所求变量的后验概率；在风险评估中，将初因事件设为证据变量，不安全结果设为查询变量，运用变量消元法更新各节点分布即可求出最终不安全结果发生的cp值。

由于混合威布尔分布参数的q函数通常为复杂的非线性方程，无法直接获得其解析值，因此，通过pso优化em算法包括：

将pso中每一个粒子的d维位置信息对应于参数θ；

在算法初始化时，根据各参数的初始值确定粒子位置与速度的更新界限，并随机生成数量为popsize的粒子种群；

将q函数作为pso的适应度函数，通过比较每代粒子的适应度值更新粒子；

判断更新粒子是否达到最大寻优次数，并在达到最大寻优次数，输出全局最优粒子位置作为极大化结果。

本申请的具体实施过程为：

如图1所示，s1、结合机队运行中得到的飞机部件失效数据，构建基于优化em算法的混合威布尔分布参数估计模型。以某个共有500架飞机的客运机队为例，在一份msad报告中指出，该机队在定期检修中发现多架飞机机翼结构存在裂纹，经初步分析，该机队裂纹是由两种原因导致的疲劳损伤的可能性较大。随着每架飞机上服役时间的累积，翼肋裂纹增长至断裂会导致多余的载荷将传递给由周围的结构部件，最终会产生飞机失控等严重后果，因此需要对受影响机队进行风险评估，由历史经验获知该疲劳裂纹失效的形状参数为β＝2.78，则构建的两重混合威布尔分布参数em模型可简化为：

maxq(θ，θ^(j))

目前可获得的裂纹增长到引发失效的寿命样本如图8中的表1所示。将寿命数据输入参数估计模型中，利用图2给出的基于pso优化的em算法流程进行模型求解，具体步骤如下：

(1)参数及粒子初始化：pso中每一个粒子的d维位置信息对应所需估计参数θ，通过线性化方法确定各参数的近似值，根据每个参数的初始值设置粒子位置界限和粒子种群数popsize以及粒子寻优最大循环次数maxiter。本案例将参数向量θ＝(α1，α2，π)作为粒子的三维位置信息，第d维的粒子位置和速度初始化公式为

xd＝xd+rand·(xd^max-xd^min)

vd＝0.1·(xd^max-xd^min)；

将得到的粒子位置同时设定为粒子个体最优位置pbest与全局最优位置gbest的初值；

(2)e-step：在期望步中需要由当前参数确定隐变量z以及q函数。首先对各粒子位置和速度进行更新，第s代粒子的速度和位置更新公式为

式中，c1、c2为调整系数，其目的是调整粒子自身经验和全局经验在运动中的作用，本案例取加速度系数均为0.5；r1、r2是(0，1)之间的随机数；w为惯性权重，常见的有线性调整权重和随机调整权重，本模型采用随机调整法为

w＝0.5+rand/2.0

依据粒子生成的模型参数，计算隐变量估计值为：

其中，ti为输入的案例失效样本，数据总数为20。记θ^(j)为第j次em迭代后得到的参数的估计值，在下一次的迭代中代入θ^(j)与计算q(θ，θ^(j))值。

(3)m-step：极大化步的核心问题是求解q(θ，θ^(j))极大化的参数θ，得到新一轮迭代的模型参数估计值。pso优化的目标函数为：

通过比较每一代粒子q函数值的大小来更新其个体与全局位置，粒子寻优至最大循环次数，得到的粒子全局最优位置即为本轮em迭代的估计结果θ^(j+1)。为了提高加权因子的估计精度，可根据当前隐变量求得：

比较更新前后q函数值的大小，取较大值作为下一轮迭代的输入加权因子。

(4)判断迭代结束条件：重复第(2)、(3)步直到到达收敛条件

||q(θ^j+1，θ^(j))-q(θ^j，θ^(j))||＜ε

式中，ε为较小的正数。满足收敛条件后输出模型最优估计参数值。

通过上述算法得到分布尺度参数估计值分别为加权因子估计值为机队裂纹失效事件的概率密度函数如图3所示，图中散点表示真实失效样本值的概率密度值。需要说明地是，当α小于飞机退役时间表明损伤发生在寿命早期阶段，与之相反，当α大于退役时间则表明损伤发展较为缓慢；

s2、根据所得分布参数，利用蒙特卡洛仿真算法模拟机队中受多重耗损失效模式影响的da数量。机队风险中对da的具体定义是如果不采取措施纠正已确定的潜在不安全状况，则在规定的时间内(通常是飞机当前役龄到退役役龄的时间)预计发生的权重事件或死亡人数。本案例中机队的飞机数量m为500架，每架飞机当前役龄状况如图4所示，役龄单位为飞行循环(flightcycle，fc)，机队的平均退役年龄tr＝50000fc。将已知参数输入蒙特卡洛模拟算法中，如图5所示，设定itermax为5000次，初始化da为0并开始循环。本案例中各子分布加权因子的累积概率分别为q1＝0.26，q1＝1，根据r属于的概率区间[q0，q1]或[q1，q2]，确定第m架飞机服从的威布尔子分布失效模式，利用逆变换法生成飞机结构件失效时间tm，判断其是否在失效预测期间[t0(m)，tr]之内，最终得到模拟结果图如图6所示，模拟得到机队da值约为153架。

s3、基于bn综合分析机队产生不安全结果条件概率。在建立bn前，必须将失效事件的基本情况明确定义为所研究的状况，通常为msad指令中的潜在不安全状况，如本案例发现的机翼裂纹事件。随着飞机服役时间的累积，翼肋裂纹增长至断裂，从而使其自身承担的载荷传递给与之相连的桁条和翼梁腹板，累积的多余载荷可能导致蒙皮以及机翼壳体失效，最终产生飞机失控等严重后果。

构建bn时，每个节点与失效传播相关事件一一对应，并用有向“边”表示从初因事件到损伤率已知的最终不安全结果的因果关系。本案例的bn由于中间节点数量较多，现实中无法获知每个节点的条件概率表，可通过增加节点状态合并相同层级的中间事件将其化简，化简前后的bn为如图7所示，图中的节点信息如图9中的表2所示，在简化后的bn中，节点a～b均为二值，“1”表示失效状态，“0”表示工作状态；d、e为多值状态。运用ve法对初始裂纹失效进行因果推理求解后验概率即cp值，具体过程为：

(1)分解bn联合分布ζ＝{p(a),p(ba),p(cb),p(dc),p(ed)}；

(2)设定证据变量a＝1，e为查询变量，得到ζ＝{p(a＝1),p(b|a＝1),p(c|b),p(d|c),p(e|d)}；

(3)对a、e外的剩余变量排序，设消元顺序为c、d、b，消去c得ζ＝{p(a＝1),p(b|a＝1),p(e|d),λ1(b,d)}，其中下一个消去的变量为d，得到ζ＝{p(a＝1),p(ba＝1),λ2(b,e)}，其中更新后的因子最后一个要消去变量是b，同理可得到ζ＝{p(a＝1),λ3(e)}，其中

(4)将消元后的ζ中相乘并归一化得到所求变量的后验概率为代入基于经验获取得相应条件概率，得到最终不安全结果全部发生的累加概率为cp＝0.00364。

s4、根据所述事件与结果因子，计算总体未纠正机队(rt)，判断其风险等级，输出定量风险评估结果。taram定义的总体未纠正机队风险是指已知潜在不安全状况时，若未对机队采取纠正措施，权重事件在受影响机队的剩余寿命中预计发生的统计次数。rt对应的商业运营的客机风险阈值为0.02，通过比较机队风险值与阈值以判断机队是否存在危险状况，其结果可以指导安全决策。

nd(notdetected)是在导致不安全的结果或情况之前，基于历史经验得到的未检测到缺陷发生的平均概率，本案例中nd取0.2；ir(injuryratio)表示结果的损伤率即严重度的量化形式，由死亡人数除以飞机上的总人数求得，通常基于运输飞机事故的历史数据，可根据美国联邦航空管理局运输飞机理事会给出的如图10中的表3确定；由于各不安全结果的损伤率不同，可按照步骤s3分别计算每个后果的条件概率值，代入rt＝da×nd×cp×ir得到最终风险结果为rt＝0.035。

根据对应的风险等级，可以判断本案例中的运输机队总体未纠正风险高于可接受风险值0.02，说明机翼裂纹初因事件影响机队运行安全，因此需要对机队制定相应的风险管理方案使其恢复到可接受水平。

在本说明书的描述中，参考术语“一个实施例”、“示例”、“具体示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节，也不限制该发明仅为所述的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。