结冰条件下人-机-环系统的飞行风险概率

薛源， 徐浩军， 胡孟权

空军工程大学 航空航天工程学院， 西安 710038

结冰条件下人-机-环系统的飞行风险概率

薛源， 徐浩军*， 胡孟权

空军工程大学 航空航天工程学院， 西安 710038

以结冰条件下的飞行风险量化概率为研究对象，基于蒙特卡罗飞行仿真实验对结冰条件下人-机-环系统的耦合特性进行了分析, 并获取了飞行参数极值样本。构建了飞行风险发生的判定条件；对飞行参数极值样本进行了统计特性分析，验证了其厚尾分布特征。一维分布类型辨识结果表明广义极值分布对相对速度和迎角极值的描述精度最高。为描述二维变量对相关性的各自影响程度，提出了一种新的双参数变权重Copula模型；辨识结果表明该Copula模型能以较高的精度通过假设检验。相关性分析的结果表明相对速度和迎角同时出现极大值和极小值的概率较大。基于二维极值样本的Copula分布模型求出了不同结冰程度下的飞行风险概率值，探讨了飞行风险的非线性增长趋势。

蒙特卡罗仿真； 结冰条件； 人-机-环系统； 飞行风险； 极值分布； Copula模型

飞机结冰较复杂的随机性和不确定性容易引起驾驶员的耦合连锁反应，导致人-机-环系统失去稳定性，从而引发飞行风险或导致飞行事故。如何评估结冰条件下的飞行风险一直是一个难题。自从航空器出现以来，由于飞机内部部件失效、驾驶员操纵水平等内在原因导致的飞行风险及飞行事故均可以被有效地解决以及避免，但结冰导致的飞行风险一直伴随着航空器的发展史，每年均有类似的飞行事故发生。根据ICAO(International Civil Aviation Organization)和美国NTSB(National Transportation Safety Board)的飞行事故数据库统计，从2002—2014年，在所有气象因素带来的飞行事故中，12%是由于飞机结冰导致的，并且其中92%是在飞行中发生积冰。近年来，最严重的一起结冰事故发生在2009年6月 1日，法国航空447号班机的空客A330客机在大西洋上空坠毁，造成216名乘客以及12名机组人员死亡，事故的原因为结冰导致自动驾驶仪关闭，飞行员随后错误操作导致失速。

结冰条件下的飞行仿真与动力学计算一直是研究外部环境风险的重要组成部分[1]。如：Bragg等提出了物理意义明显、结构简单的描述积冰影响的模型[2]，可用于估算积冰后飞机的气动参数；Hui等通过对Convair 580飞机的建模研究指出，积冰气象条件的微观处理参数与积冰后气动导数的变化具有特定的相关性，可以使用无量纲的参数去体现积冰后配平以及稳定性和操纵性导数的变化[3]；Krzysztof等模拟了飞机在积冰环境的爬升过程，再现了该飞机失事的过程[4]；Frank等对三维机翼积冰后失速迎角与最大升力系数降低量的关系进行了分析[5]；Thomas等利用“双水獭”飞机积冰试飞数据，应用多种飞行策略对多项飞行导数、操纵导数进行了研究，并着重分析了机翼和平尾的失速问题[6]；王明丰等在突风扰动情况下，分别对无积冰、弱积冰和严重积冰情况下飞机飞行参数的时间响应进行了计算[7]；Lampton和Valasek利用简化的积冰影响模型，对左右机翼非对称积冰进行了仿真，同时还进行了纵向、横向和全机积冰条件的飞机仿真分析，通过比较积冰前后飞行参数的变化分析积冰对飞机稳定性和操纵性的影响[8-10]；美国伊利诺伊州大学的Bragg等提出并领导了飞机智能积冰系统研究[11]；Robert和Glen开发了结冰遭遇飞行仿真模拟器[12]；David等对结冰后的飞行包线保护系统进行了研究[13-15]。

上述工作具有较高的学术价值，但在结冰条件下飞行风险发生概率的评估与预测方面尚缺乏相关的理论与方法，无法量化评估结冰条件下动态的飞行风险概率。而结冰条件下的飞行风险评估也正是当前飞机安全性和适航性工作重点关注的内容，是完善和改进现有的飞行安全分析理论亟待解决的问题。现有的SAE ARP 4761[16]、SAE ARP 4754A[17]、MIL-HDBK-516B[18]、MIL-STD-882E[19]和MIL-STD-1797B[20]等安全性规范与指南对于评估由于飞机内部硬件系统失效导致飞行事故概率有明确的思路及方法。国内外学术界对飞行安全与飞行风险评估的理论与方法也有较多的研究工作[21-25]。但其大多是定性的分析或者是静态的可靠性评估，此类方法在评估结冰条件下的飞行风险概率方面具有明显的局限性：一是由于此类方法基于确定性的模型，其不能反应结冰条件下多因素耦合复杂情形的随机性与不确定性；二是不能基于航空器的飞行过程评估动态的飞行风险大小，对于结冰导致的飞行风险较难进行量化概率指标的评估预测。

另一方面，Copula理论成为多元极值问题的研究热点[26]，在金融、灾害预警、信号处理等方面均有应用[27-36]。由于评估飞行风险需研究多元飞行参数极值(飞行参数极值)的尾部相关性，故较适合将多元极值Copula模型应用到其中。但尚未见到将这一有效理论应用到结冰飞行风险评估中的公开报道。鉴于此，本文首先对直观反映结冰后运动状态的飞行参数进行分析。由于结冰后飞行参数的变化可以直观反映人-机-环系统的运动状态，飞行参数极值的超限情形即是飞行风险发生的判定条件；因此，在对一维极值分布进行研究的基础上，对结冰条件下基于多元极值Copula模型的飞行风险概率量化评估方法进行探索；拟在有效解决一元极值模型局限性的基础上利用多元极值分布求得结冰空间中的飞行风险概率。

1 基于蒙特卡罗仿真实验提取结冰条件下的飞行参数极值

一方面，严重结冰条件下的试飞具有较大的飞行风险；另一方面，研究结冰条件下的随机性与不确定性需要的数据量较大。因此文中基于蒙特卡罗法考虑不同程度结冰条件下的物理特性，对结冰后的动态过程进行多次计算迭代，从而提取反映人-机-环系统在结冰后物理特性与随机不确定性的飞行参数极值样本点。为下文中基于多元极值Copula的飞行风险评估方法研究提供数据输入。飞行参数极值的提取系统基于某型飞机的地面试验系统改造而成，蒙特卡罗仿真实验的流程如图1所示。

首先，设定结冰程度，提取特定结冰程度下的结冰气动模型；然后，利用蒙特卡罗法将结冰条件下的内外部环境变量按照其出现频率进行随机抽样，从而对每次计算迭代过程中所使用的内外部环境参数产生影响；同时，飞行员对每一次蒙特卡罗仿真实验的迭代计算过程进行操纵，从而提供结冰后的舵面输入信号，对相关的气动参数及操纵信号产生量化影响。如此便可以通过多次的循环迭代反映人-机-环系统在结冰条件下的随机性与不确定性。飞机本体方程为基于四元数法的六自由度方程，微分算法为四阶龙格库塔(Runge-Kutta)算法，仿真实验的时间步长为20 ms。气动数据与基于某型机的全机测力风洞试验，所得到的各气动参数精度均满足合格或先进指标的要求。

结冰对飞机气动参数的影响模型是飞行力学建模的一个至关重要的环节。在飞行仿真计算中常用Bragg等提出的积冰影响模型：

C(A)ice=(1+η(t)KCA)C(A)

(1)

图1 基于蒙特卡罗仿真的飞行参数极值提取过程Fig.1 Extreme flight parameters extraction process based on Monte Carlo simulation

2 结冰条件下飞行参数分析与风险判定

蒙特卡罗仿真实验的初始条件为：高度5 000 m，速度120 m/s，某型飞机以配平状态进入中度结冰区域，η(t)的变化范围为0～0.35，开始结冰时全机质量为49 000 kg，蒙特卡罗仿真的次数为150次。选取两次较典型的飞行员分别在正确与错误操纵状态下飞机结冰后的响应为案例进行研究。一次为第53次，飞行员在结冰后介入延迟较大且操纵不当情形下飞机的响应；一次为第91次，飞行员及时介入且操纵恰当情形下飞机的响应。图2为第53次和第91次蒙特卡罗仿真实验中飞行参数的变化情况。从图2中可看到在第53次迭代计算中，当飞机进入结冰空域后，结冰情形逐渐加剧，飞行员在55.9 s时感受到过载和俯仰角的变化，开始操纵驾驶杆进行修正以期重新回到配平状态，但由于飞行员介入的延迟较大，加之修正时的操纵动作过猛，迎角与速度并未达到修正预期，同时引发了横向滚转效应，偏航角开始改变，飞行员在拉杆的同时还需横向压杆以平衡姿态。在111.2 s时，飞行员注意到速度有加速减小的趋势，开始操纵油门杆增加推力，138.7 s时，迎角超出极限值，飞机进入失速振荡，速度继续减小，飞行员此时将驾驶杆拉到极限并保持油门最大位置，但仍然无法改出失速，此刻即标志着飞行风险的发生。在第91次迭代计算中，飞机在50 s时结冰效应开始显现，飞行员在53.7 s时及时介入，平稳操纵驾驶杆稳定姿态，并提前推油门以保持速度和高度，迎角和速度在短暂变化后趋于平稳，飞机未进入风险状态。第53次和第91次迭代计算所提取的迎角和速度极值[αmax,vmax]分别为[21.332 3°,63.247 8 m/s]和[6.335 9°,116.788 3 m/s]，如图2所示。

150次蒙特卡罗仿真实验后所提取的150组极值样本散点分布如图3所示，其中相对速度为飞机配平速度与结冰速度的比值，将这个比值的极大值记为vstable/vmin，文中vstable=120 m/s。如此可以将对速度极小值的研究转换为对相对速度极大值的研究。图3中标注为上文中第53次和第91次蒙特卡罗仿真实验的极值样本。

图2 第53次和第91次蒙特卡罗仿真实验中飞行参数变化情况Fig.2 Flight parameter changes in the 53rd and 91st Monte Carlo simulation experiments

图3 150组极值样本散点分布Fig.3 Scatter distributions of 150 extreme samples

根据150次蒙特卡罗迭代计算的结果可确定迎角或者速度超限为飞行风险发生的判定条件。对文中所涉及到的飞行风险进行定义如下：以超过95%的概率极易引起MIL-STD-882E[19]中所定义的风险范畴中评估值为1-5的灾难性飞行事故。即不能安全飞行和着陆的失效情况，引起飞机结构损伤并导致至少一人的伤亡。给出判断文中定义的尾流飞行风险是否发生的判定条件为

(2)

3 一维飞行参数极值样本的统计特性

从图3中可以看出滚转角极值、下降高度极值和迎角极值的分布都存在明显的厚尾特性，这种分布形式在低频高危事件(如地震、海啸、金融风险、飞行事故等)中较为常见，文中涉及遭遇尾流情形下的飞行风险亦属于低频高危事件的范畴。针对此种分布形式，目前较有效的描述方法为采用极值理论。下文将针对飞行参数极值样本分布模型的不确定性，对尾流遭遇情形下一维飞行参数极值的尾部分布规律进行研究，采用最优化辨识算法确定分布模型中的未知参数，并利用辨识后的分布模型对一维飞行参数极值进行拟合优度检验，从而找出最合适的分布类型。在此之前，首先对飞行参数极值变量的统计特性进行研究。对飞行参数极值变量进行统计分析，结果如表1所示。

表1 极值样本统计量Table 1 Statistics of extreme samples

观察表1可以发现，相对速度极值与迎角极值的最大值均比最小值偏离均值与中位数的程度要大，说明其分布形式并不是左右对称的。继续分析表1可以看到三组极值样本滚转角、下降高度、迎角的峰度系数均大于3，说明三组极值样本均有比正态分布更长的尾部；其偏斜度均大于0°，表明分布类型在右侧具有较长尾部。

图4为三组极值样本的盒形图，盒形中加粗的线为样本的中位数，上下边界分别为样本正态分布范围的25%与75%界限。可以看出两组极值参数在上尾均有多个样本点超出正态分布的界限范围。显然，亦可得到样本点的分布具有上厚尾特征的结论。

图4 极值样本盒形图Fig.4 Box graph of extreme samples

4 一维飞行参数极值样本的分布类型辨识

4.1 分布类型假设

针对图3中的极值样本，初步推断广义极值分布(式(3))较适合描述其厚尾分布特性。为验证广义极值分布的有效性，采用广义极值分布模型和另外5种可以描述厚尾特性的分布模型来辨识一维极值参数，从而寻找对一维极值参数样本分布的描述最准确的分布模型。所采用的分布模型主要有：广义极值分布(GEV)、正态分布模型(Normal)、对数正态分布(Lognormal)、威布分布(Weibull)、指数分布(EXP)、极值分布(EV)。其中广义极值分布模型为

(3)

式中：ξ∈R，μ∈R，σ>0，1+ξ(x-μ)/σ>0。μ为位置参数；σ为尺度参数；ξ为形状参数。

4.2 辨识结果

利用自适应粒子群优化(Adaptive Range Particle Swarm Optimization，ARPSO)算法辨识出一维极值分布模型的未知参数，如表2所示。在表2中，GEV的待辨识依次为ξ、μ和σ，Normal的待辨识参数为μ和σ，Lognormal的待辨识参数为μ和σ，Weibull的待辨识参数为a和b，EXP的待辨识参数为μ，EV的待辨识参数为μ和σ。

4.3 GEV描述一维极值分布的准确性验证

表2 极值参数的辨识结果Table 2 Identification results of extreme parameters

表3 相对速度极值的拟合优度检验Table 3 Fitting test for extreme relative velocity

表4 迎角极值的拟合优度检验Table 4 Fitting test for extreme angle of attack

从图5和图6中亦可看出GEV对极值样本的概率密度和累积概率的拟合最为准确，其他分布形式不能充分反映极值样本的分布信息。

图5 相对速度极值的概率密度图与累积概率图Fig.5 Probability density graph and cumulative probability graph of extreme relative velocity

图6 迎角极值的概率密度图与累积概率图Fig.6 Probability density graph and cumulative probability graph of extreme angle of attack

(4)

(5)

根据式(2)的飞行风险发生条件，将表2中的GEV分布参数ξ、μ和σ代入到式(6)的F2(αmax)中，再根据式(7)求得以迎角αmax为风险评估参数时发生飞行风险的概率为0.118 0。

(6)

Pr=1-F2(αmax/αc(δf,Ma,η)>1)

(7)

由上文可知，式(5)和式(7)求出的飞行风险概率值是不相同的，这说明了由于飞行风险的判定需考虑多个具有决定性作用的极值参数，利用一维极值参数与一元极值GEV分布去求飞行风险概率不能全面地考虑其他飞行参数极值对飞行风险的影响，具有一定的局限性。而多元极值理论可以考虑到所有对飞行风险发生起决定作用的飞行参数极值，这也正是下面所探讨的内容。

5 结冰条件下的Copula模型构建

上文对一维极值分布进行了辨识，并分别给出了相对速度和迎角作为单变量临界极值时的飞行风险概率。但根据式(2)，文中结冰后飞行风险的判定牵扯到二维极值参数的相关性结构。有关单变量的结论并不能平行推广到多元情形，其主要原因是多维随机变量的联合分布除了与各分量的边缘分布有关之外，更重要的是与变量之间相关性的关系。当随机变量的个数比较大时，单个分量的极值行为未必含有整个向量的联合极值行为。因此，下文基于Copula理论构建二维极值的相关结构，由于Copula理论是针对极值分布的评估而提出的，故在对极值相关性的描述上，Copula极值分布模型能较好地反映极值参数之间的联系和发展趋势。

5.1 二维分布模型假设

(8)

(9)

式中：Copula模型的未知参数ω≥1。当ω=1时，表示极值参数u与v相互独立；ω=+∞表示u与v完全相关。

(10)

式中：ω≠0;如果ω>0，表示u与v正相关;ω→0，表示u与v相互独立;ω<0，表示u与v负相关。

C(u,v)=(u-ω+v-ω-1)-1/ω

(11)

式中：ω∈(0,∞)，当ω→0时，u与v趋于独立；ω=+∞表示u与v完全相关。

(12)

式中：ω≥1，当ω=1时，极值参数u与v相互独立；当ω=+∞时，极值参数u与v完全相关。

C(u,v)=1-[(1-u)ω+(1-v)ω-

(1-u)ω(1-v)ω]1/ω

(13)

式中：ω≥1，当ω=1时，极值参数u与v相互独立；当ω=+∞时，极值参数u与v完全相关。

5.2 提出二维极值变量的DAWP Copula模型

根据上文可以初步判定对上尾变化敏感的Gumbel Copula模型及Joe Copula模型最能反映文中极值的分布情况，但其单参数的特征使得在描述二维变量对相关性的各自影响程度时具有一定的局限性，故本文在Gumbel模型的基础上提出一种新的双参数变权重Copula模型(Double Adaptive Weight Parameter Copula,DAWP Copula)，如式(14)所示。

(14)

5.3 二维Copula模型的未知参数辨识

(15)

(16)

2) 求出Copula的密度函数为

c(ui,vi;ω1,ω2，…,ωd)=

(17)

3) 对于式(17)中未知参数ω1,ω2,…,ωd的估计，结合似然函数与χ2分布检验构建目标函数：

H(u1,u2,…,un;v1,v2,…,vn;ω1,ω2,…,ωd)=

(18)

式中：Ai为极值样本点在区间(u≤ui,v≤vi)内的个数。

4) 利用非线性辨识算法求出未知参数的值，如表5所示。

表5Copula模型中未知参数辨识结果

Table5IdentificationresultsofunknownparametersinCopulamodels

CopulamodelUnknownparameterGumbelω=2.3353Frankω=6.6261Claytonω=1.3951GSω=1.5058Joeω=3.0506DAWPω1=-0.9159,ω2=-0.6137

5.4 二维Copula模型的拟合优度检验

对于涉及到的Copula模型，分别应用AIC(Akaike Information Criteria)准则、BIC(Bayesian Information Criteria)准则、χ2检验、K-S检验评价其拟合优度，如表6所示。

表6Copula模型对极值样本的拟合优度检验

Table6FittingtestofextremesamplefordifferentCopulamodels

CopulamodelAICBICχ2K-SPvalueGumbel-149.6779-146.66730.22330.09690.2653Frank-118.0216-115.01100.47240.21720.0764Clayton-86.7373-83.72660.36930.36470.0081GS-122.7043-119.69360.13010.20130.0235Joe-143.6251-140.61441.43800.11560.2693DAWP-149.2077-146.19710.18680.08120.3139

从表6中可以看出，Gumbel Copula、Joe Copula、DAWP Copula的P值均大于显著性水平0.01，0.02，0.05，即这3种Copula在99%、98%、95%的置信水平下均能通过检验；而Clayton Copula甚至在99%的置信水平下亦未能通过检验。继续比较发现Gumbel Copula和DAWP Copula的AIC与BIC检验值比其他Copula的值均要小，说明二者的辨识精度较高。而DAWP Copula的K-S值为0.081 2，比Gumbel Copula的K-S值还要小，同时其P值亦是所有 Copula 模型中最大，说明其对极值参数相关结构的描述更为准确。

5.5 相关性分析

5.6 基于DAWP Copula模型求飞行风险概率

根据风险判定式(2)及构建的DAWP Copula模型求出风险概率为

F2(αmax/αc(δf,Ma,η)>1))

(19)

CDAWP代表基于式(14)的DAWP Copula模型。根据表5中辨识出的DAWP Copula未知参数求得在第1节提到的仿真特定点上的飞行风险概率Pr=0.130 4。对比通过一维飞行参数极值GEV模型(式(5)和式(7))求出的飞行风险概率，可发现基于多元极值Copula模型求出的飞行风险概率值比式(5)和式(7)求出的风险概率值要大，说明其充分考虑了不同飞行参数极值超限时的飞行风险状况，相对于一维飞行参数极值判定飞行风险更全面、更精确。

6 不同结冰程度下的飞行风险

不同冰型、不同结冰位置、不同结冰厚度等会对气动参数产生不同的量化影响，但其均会使得气动参数恶化，使飞机的飞行品质降低，操控性变差。针对文中所研究的背景飞机，在相同的高度速度以及外部环境参数下，按照不同结冰情形对气动参数的影响程度，将无结冰状态到极严重结冰状态均分为10个等级。无结冰等级为0，此时对气动参数无影响；极严重结冰等级为9，此时结冰导致气动参数极具恶化，对于式(1)来说，等级9意味着参数η将大于0.3，阻力系数增加300%以上，升力系数减小30%以上。采用上文方法基于蒙特卡罗仿真实验提取每个状态下的迎角与速度极值样本，利用DAWP Copula模型对样本进行辨识后计算出每个状态发生飞行风险的概率如表7和图7所示。表7中的等级7对应于5.6节案例中所求的飞行风险0.130 4。

表7不同结冰等级下的飞行风险概率

Table7Flightriskprobabilitiesunderdifferenticinglevels

IcinglevelFlightriskprobability00.00001910.00020620.00091530.00138140.00617250.01559360.05682970.13041180.25337590.360572

图7 不同结冰等级下的飞行风险概率曲线Fig.7 Flight risk probability curve under different icing levels

由表7和图7可以看出，无结冰状态发生飞行风险的概率为0.000 019，基于MIL-STD-882E中飞行事故率的划分可认为此时没有飞行风险。从状态0到状态5，飞行风险概率平稳增大，状态5 的飞行风险为0.015 593，已经高于SAE ARP-4761所规定系统安全性水平中的B级危险性失效状态概率，在MIL-STD-882E从1～20的风险评级中定为2级。从状态5到状态9，也就是中度结冰到极严重结冰这一区间，飞行风险概率呈指数形态激增。状态8严重结冰时的风险概率已经高达25.337 5%，此概率值远远高于MIL-STD-882E中所规定A级(经常发生)事件10%的出现概率，在MIL-STD-882E从1～20的风险评级中定为1级。这一阶段需要驾驶员提前做出反应并正确平稳操纵飞机才能避免飞行风险的发生。从图7中也可以看出，从无结冰到极严重结冰，飞行风险概率的增长不是一个线性的过程，其在轻度结冰时的数值较小，在中度结冰状态之后开始激增，为一种新的飞行风险变化趋势。

7 结 论

1) 考虑结冰条件下的不确定性与随机性因素，利用蒙特卡罗仿真实验研究了人-机-环系统在结冰条件下的动力学特性。结果表明在中度结冰条件下某型飞机的动力学特性逐步恶化，迎角和法向过载增量较大，易于超出临界值以致飞机失控；同时由于阻力系数的增加和升力系数的减小，飞行速度和飞行高度逐渐降低。如驾驶员在此时介入延迟较大且操纵不当，则极易引起飞行参数极值超限从而引发飞行风险；驾驶员需在中度结冰条件下迅速介入并平稳操纵以脱离飞行风险。

2) 验证了所提取极值样本具有厚尾分布的特性。拟合优度检验的结果表明GEV模型能以较高的精度对飞行参数极值样本的分布特性进行描述，从而验证了迎角极值和相对速度极值均符合GEV分布。结合多元极值Copula模型提出了一种量化评估飞行风险概率的新方法，并将其用到结冰条件下的风险概率量化评估中。为了解决单参数Copula模型在描述多维变量对相关性的各自影响程度时具有的局限性，提出了二维的DAWP Copula模型。基于DAWP Copula模型对某型飞机不同结冰等级下的飞行风险进行了计算。结果表明从无结冰到极严重结冰状态，飞行风险概率的增长不是一个线性的过程，为一种新的变化趋势。

本文思路可对其他内外部环境因素影响下的飞行风险概率评估提供参考。其中的飞行风险概率评估方法是对现有各类飞行安全规范中风险评估理论的有效补充，对于飞行安全与适航性管理具有积极的作用。需注意的是，由于飞行事故的发生是一个多因素影响的不确定过程，不可能将所有内外部随机因素考虑完全。因此，本文得到的飞行风险量化概率在多数状况下是一种类似于飞行安全指南(如SAE ARP-4761与MIL-STD-882E)中事故率的参考值，与真实值必然有一定的误差。但其在不同状况下飞行风险的横向对比分析、风险程度的归类划分中具有积极的意义。同时，由于本文侧重于风险评估方法，主要考虑了结冰对气动参数的影响，没有全面考虑如发动机结冰、传感器结冰以及舵面结冰卡滞等其他结冰风险因素。因此，下一步工作将利用文中方法综合考虑飞机结冰后多种不利因素相互耦合交联的情形。

[1] CONNOR P, DEA A, KENNEDY Q, et al. Measuring safety climate in aviation: A review and recommendations for the future[J]. Safety Science, 2011, 49(2): 128-138.

[2] BRAGG M B, PERKINS W R, SARTER N B, et al. An interdisciplinary approach to inflight aircraft icing safety: AIAA-1998-0095[R]. Reston: AIAA, 1998.

[3] HUI K, WOLDE M, BROWN A. Flight dynamics model of turboprop transport aircraft icing effects based on preliminary flight data: AIAA-2005-1068[R]. Reston: AIAA, 2005.

[4] KRZYSZTOF S, MACIEJ L, EDYTA L, et al. Aircraft flight dynamics with simulated ice accretion: AIAA-2004-4948[R]. Reston: AIAA, 2004.

[5] FRANK T L, ADBOLLAH K. Effects of ice accretions on aircraft aerodynamics[J]. Progress in Aerospace Sciences, 2001, 37(8): 669-767.

[6] THOMAS P R, BILLY P B, SAM L. Current methods for modeling and simulating icing effects on aircraft performance stability and control: AIAA-2008-6204[R]. Reston: AIAA, 2008.

[7] 王明丰, 王立新, 黄成涛. 积冰对飞机纵向操稳特性的量化影响[J]. 北京航空航天大学学报, 2008, 34(5): 592-595.

WANG M F, WANG L X, HUANG C T. Computational effects of ice accretion on aircraft longitudinal stability and control[J]. Journal of Beijng University of Aeronautics and Astronautics, 2008, 34(5): 592-595 (in Chinese).

[8] LAMPTON A, VALASEK J. Prediction of icing effects on the dynamic response of light airplanes[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(3): 722-732.

[9] LAMPTON A, VALASEK J. Prediction of icing effects on the coupled dynamic response of light airplanes[J]. Journal of Guidance, Control, and Dynamics, 2008, 31(3): 656-673.

[10] LAMPTON A, VALASEK J. Prediction of icing effects on the lateral directional stability and control of light[J]. Aerospace Science and Technology, 2012, 23: 305-311.

[11] BRAGG M B, BASAR T, PERKINS W R, et al. Smart icing systems for aircraft icing safety: AIAA-2002-0813[R]. Reston: AIAA, 2002.

[12] ROBERT W D, GLEN A D. Icing encounter flight simulator[J]. Journal of Aircraft, 2006, 43(5): 1528-1537.

[13] DAVID R G, BILLY B, RICHARD R, et al. Development and implementation of a model-driven envelope protection system for in-flight ice contamination: AIAA-2010-8141[R]. Reston: AIAA, 2010.

[14] RICHARD R, BORJA M, BILLY N, et al. Piloted simulation to evaluate the utility of a real time envelope protection system for mitigating in-flight icing hazards: AIAA-2010-7987[R]. Reston: AIAA, 2010.

[15] DAVID R G. Requirements and modeling of in-flight icing effects for flight training: AIAA-2013-5075[R]. Reston: AIAA, 2013.

[16] Society of Automotive Engineers. Guidelines and methods for conducting the safety assessment process on civil airborne systems and equipment： ARP 4761[S]. Washington, D.C.: SAE,1996.

[17] Society of Automotive Engineers. Certification considerations for high-integrated or complex aircraft systems： ARP 4754[S]. Washington, D.C.: SAE, 2010.

[18] USA Department of Defense. Airworthiness certification criteria： MIL-HDBK-516B[S]. Washington, D.C.: DOD, 2005.

[19] USA Department of Defense. Standard practice for system safety： MIL-STD-882E[S]. Washington, D.C.: DOD, 2012.

[20] USA Department of Defense. Flying qualities of piloted aircraft： MIL-STD-1797B[S]. Washington, D.C.: DOD, 2012:673-695.

[21] BROOKER P. Experts, Bayesian belief networks, rare events and aviation risk estimates[J]. Safety Science, 2011, 49(8): 1142-1155.

[22] WANG W H, JIANG X B, XIA S C. Incident tree model and incident tree analysis method for quantified risk assessment: An in-depth accident study in traffic operation[J]. Safety Science, 2010, 48(10): 1248-1262.

[23] MATTHEWS B, DAS S, BHADURI K, et al. Discovering anomalous aviation safety events using scalable data mining algorithms[J]. Journal of Aerospace Information Systems, 2013, 10(10): 467-475.

[24] OCAMPO J, MILLWATER H, SINGH G, et al. Development of a probabilistic linear damage methodology for small aircraft[J]. Journal of Aircraft, 2011， 48(6): 2090-2106.

[25] BALACHANDRAN S, ATKINS E M. A constrained Markova decision process framework for flight safety assessment and management: AIAA-2015-0115[R]. Reston: AIAA, 2015.

[26] NELSEN R B. An introduction to copulas[M]. 2rd ed. New York: Springer, 2006: 51-108.

[27] DIKS C, PANCHENKO V, SOKOLINSKIY O, et al. Comparing the accuracy of multivariate density forecasts in selected regions of the copula support[J]. Journal of Economic Dynamics & Control, 2014, 48: 79-94.

[28] SUKCHAROEN K, ZOHRABYAN T, LEATHAM D, et al. Interdependence of oil prices and stock market indices: A copula approach[J]. Energy Economics, 2014, 44: 331-339.

[30] YASMIN S, ELURU N, ABDUL R, et al. Examining driver injury severity in two vehicle crashes—A copula based approach[J]. Accident Analysis and Prevention, 2014, 66(3): 120-135.

[31] MASIN M, LAMBERTI A, ARCHETTI R. Coastal flooding: A copula based approach for estimating the joint probability of water levels and waves[J]. Coastal Engineering, 2015, 97: 37-52.

[32] BESSA R J, MIRANDA V, BOTTERUD A, et al. Time-adaptive quantile-copula for wind power probabilistic forecasting[J]. Renewable Energy, 2012, 40(1): 29-39.

[33] ERYILMAZ S. Estimation in coherent reliability systems through Copulas[J]. Reliability Engineering and System Safety, 2011, 96(5): 564-568.

[34] BERGER T. Forecasting value-at-risk using time varying copulas and EVT return distributions[J]. International Economics, 2013, 133: 93-106.

[35] MOAZAMI S, GOLIAN S, KAVIANPOUR M R, et al. Uncertainty analysis of bias from satellite rainfall estimates using copula method[J]. Atmospheric Research, 2014, 137(2): 145-166.

[36] XUE Y, XU H J, WANG X L. Build probability distribution maps of flight risk during wake encountering[J]. Journal of Aircraft, 2015, 52(3): 805-818.

薛源男， 博士， 讲师。主要研究方向： 飞行仿真与飞行安全。

Tel.: 029-84787637

E-mail: xuhaojun@xjtu.edu.cn

徐浩军男， 硕士， 教授， 博士生导师。主要研究方向： 飞行安全与作战效能。

Tel.: 029-84787637

E-mail: xuhaojun@xjtu.edu.cn

*Correspondingauthor.Tel.：029-84787637E-mail:xuhaojun@xjtu.edu.cn

Flightriskprobabilityofpilot-aircraft-environmentsystemundericingconditions

XUEYuan,XUHaojun*,HUMengquan

AeronauticsandAstronauticsEngineeringCollege,AirForceEngineeringUniversity,Xi’an710038,China

Thequantitativeflightriskprobabilityundericingconditionsissetastheresearchobject.BasedonMonteCarloflightsimulationexperiment,couplingcharacteristicsofpilot-aircraft-environmentsystemaftericingarestudied,andtheextremeflightparametersareextracted.Thejudgementconditionsforflightriskaregiven.Andstatisticalpropertiesoftheextremesamplesareanalyzed.Thefactthatextremeflightparametershaveheavytaildistributioncharacteristicsisverified.Identificationresultsindicatethatgeneralizedextremevaluemodelcanbestdescribethedistributioncharacteristicsofrelativevelocityandangleofattack.Inordertodescribetheinfluencelevelofeachtwo-dimensionalparameteronthecorrelations,anewCopulamodelthathastwochangeableweightsisproposed;identificationresultshowsthatthisnewCopulamodelcanpassthehypothesistestingwithhighaccuracy.Correlationanalysisresultsrevealthatrelativevelocityandangleofattackwouldbothshowupextremevalueswithgreaterprobability.FlightriskprobabilitiesunderdifferenticinglevelsarecalculatedbasedonCopuladistributionmodelwithtwo-dimensionalextremevalues.Thenonlineargrowthtrendofflightriskisalsodiscussed.

MonteCarlosimulation;icingconditions;pilot-aircraft-environmentsystem;flightrisk;extremevaluedistribution;Copulamodel

2016-01-03；Revised2016-02-22；Accepted2016-03-16；Publishedonline2016-05-041416

URL：www.cnki.net/kcms/detail/11.1929.V.20160504.1416.008.html

s：NationalBasicResearchProgramofChina(2015CB755802)；NationalNaturalScienceFoundationofChina(61503406)

2016-01-03；退修日期2016-02-22；录用日期2016-03-16； < class="emphasis_bold">网络出版时间

时间：2016-05-041416

www.cnki.net/kcms/detail/11.1929.V.20160504.1416.008.html

国家“973”计划 (2015CB755802)； 国家自然科学基金 (61503406)

*

.Tel.：029-84787637E-mailxuhaojun@xjtu.edu.cn

薛源， 徐浩军， 胡孟权． 结冰条件下人-机-环系统的飞行风险概率J． 航空学报,2016,37(11):3328-3339.XUEY,XUHJ,HUMQ.Flightriskprobabilityofpilot-aircraft-environmentsystemundericingconditionsJ.ActaAeronauticaetAstronauticaSinica,2016,37(11):3328-3339.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0086

V212; X949

A

1000-6893(2016)11-3328-12