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

结冰条件下的飞行仿真与动力学计算一直是研究外部环境风险的重要组成部分^[1]。目前，国内外在结冰后飞机动力学特性的分析方面有较多的研究工作。如：Bragg等^{[2, 3]}提出了可用于估算结冰后飞机的气动参数的积冰影响的模型；文献^[4]研究了Convair 580飞机在结冰后稳定性和操纵性导数的计算方法，分析了气象条件和结冰气动导数的相关性；Krzysztof等在文献^[5]中模拟了飞机在积冰环境的爬升过程，再现了该飞机失事的过程；Frank等^[6]对三维机翼积冰后失速迎角与最大升力系数的降低量的关系进行了分析；文献^{[7, 8]}基于试飞数据研究了多项飞行导数、操纵导数，分析了结冰条件下的机翼和平尾失速问题；文献^[9]依据气象条件、飞机特征以及飞行条件建立了结冰飞机的飞行动力学模型，在此基础上，对飞机结冰后的动力学特性进行了计算；文献^[10-11]从建立结冰翼型气动系数增量函数出发，通过对全机在不同冰形下的气动性能的计算，研究了不同冰形对飞机飞行性能影响的程度，并得到了结冰后的纵向操稳特性；Lampton^{[12, 13, 14]}采用飞机结冰参量模型，以小扰动理论为基础研究了结冰对轻型飞机纵向及横航向稳定性和操纵性的影响，有助于驾驶员实施合理的操纵策略；美国伊利诺伊州大学的Bragg提出并领导了飞机智能积冰系统（Smart Icing System-SIS）研究^[15]。Robert和Glen等人开发了结冰遭遇飞行仿真模拟器^[16]。David和Richard等人对结冰后的飞行包线保护系统进行了研究^{[17, 18, 19]}。

但是上述研究大多是对特定结冰情形下的动力学特性进行单次的仿真计算分析，因此其结论并不能完全反映结冰条件下人-机-环系统中的随机性和不确定性。鉴于此，论文基于蒙特卡罗法构建结冰条件下的飞行动力学仿真计算系统，通过多次的迭代计算获得飞行参数样本，以达到研究飞行安全和飞行风险所需统计学上的可靠度。从而为结冰条件下的飞行风险评估提供相对准确的数据输入。

结冰条件下内外部影响因素具有复杂随机的特性，主要体现在：一是复杂的物理特性，这里所说的物理特性是指人员操纵、航空器运动能用解析方法描述部分的特性，一般具有确定性的数学模型，其复杂性主要是指各部分中的非线性；二是复杂的随机特性，是指人-机-环系统中不能用解析方法描述部分的特性，包括外部环境因素、驾驶员个体操纵的差异、航空器故障的随机发生、以及结冰状态的多样性等等，其复杂性主要指各部分的不确定性或随机性。研究结冰条件下内外部影响因素的复杂随机性需要的数据量较大，试飞数据与人在回路地面实验数据无法满足数据量的要求；同时，空中结冰飞行风险较大，而且需模拟条件众多，由此导致所需实验条件苛刻。以上原因使得采用蒙特卡罗法仿真的手段成为较为有效的途径。本章首先基于分布式框架构建了人-机-环仿真系统，而后研究了飞行仿真中典型模型的构建技术。基于蒙特卡罗法考虑不同程度结冰条件下的物理特性，对结冰后的动态过程进行多次计算迭代，从而提取反映人-机-环系统在结冰后物理特性与随机不确定性的飞行参数极值样本点。为下文中基于多元极值Copula的飞行风险评估方法研究提供数据输入。飞行参数极值的提取系统基于某型飞机的地面试验系统改造而成，实物如图 1所示。

图 1 提取极值参数的仿真系统 Fig. 1 The simulation system of extracting extreme parameters

图选项

蒙特卡罗仿真实验的流程如图 2图蒙特卡洛流程所示。

图 2 结冰条件下基于蒙特卡罗仿真实验的飞行参数极值提取过程 Fig. 2 Flight parameter extraction process under icing conditions based on Monte Carlo simulation

图选项

首先，设定结冰程度，提取特定结冰程度下的结冰气动模型；而后，利用蒙特卡罗方法将结冰条件下的内外部环境变量按照其出现频率进行随机抽样，从而对每次计算迭代过程中所使用的内外部环境参数产生影响；同时，飞行员对每一次蒙特卡罗仿真实验的迭代计算过程进行操纵，从而提供结冰后的舵面输入信号，对相关的气动参数及操纵信号产生量化影响。如此便可以通过多次的循环迭代反映人-机-环系统在结冰条件下的随机性与不确定性。在本文中，为保证数据准确性，驾驶员输入采用真实人在回路的补偿操纵模型。补偿操纵模型基于飞机的操稳特性与驾驶员的操纵行为反应，可以对驾驶员遇到突发情况下的反映进行精确建模，是人-机-环系统研究中最常用的模型，如图 3所示，其适用于需精确控制飞机飞行姿态的结冰情形。

图 3 驾驶员补偿操纵模型 Fig. 3 The pilot compensation control model

图选项

图 2中蒙特卡罗仿真实验的迭代计算步骤可概括如下：

（1）设置飞机的机体参数、飞机初始飞行状态；

（2）设置i=1；

（3）利用蒙特卡罗法根据各个条件变量的统计频率特征抽样结冰条件下人-机-环境系统中的随机变量值；

（4）进行第i次飞行仿真计算；

（5）记录第i次计算结果中的飞行参数，提取极值点数据，存入数据库；

（6）i=i+1,回到步骤4，循环迭代，直至i=n。

其中n为设置的仿真特定点上的蒙特卡罗计算次数。

为避免奇点的产生，飞机本体模型采用基于四元数法的六自由度方程，微分算法为四阶龙格库塔（Runge-Kutta）算法，仿真实验的时间步长为20 ms。四元数是基于飞行器在空间的姿态是通过机体轴系相对地面轴系的三个欧拉角［φ,θ,ψ］来构建的。文中四元数与欧拉角的关系如公式(1)。将四元数代入到运动方程中，对各项求导，并考虑欧拉角速率与旋转角速度的关系，构建最终的六自由度全量方程如公式(2)-公式(7)；在求解方程时，四元数的修正公式如公式(8)，四元数的导数形式最终可写成公式(9)的形式。

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

速度分量u,v,w与角速度p,q,r是基于体轴系给出的。$\bar{X}$,$\bar{Y}$,$\bar{Z}$为气动力在体轴系内的三个分量。$\bar{L}$,$\bar{M}$,$\bar{N}$为分别为滚转、俯仰、偏航气动力矩。发动机推力F_T=f(H,V,T,Thro)。

式中I_xx，I_yy和I_zz为飞行器对O_x轴、O_y轴和O_z轴的惯性矩，分别为

点击浏览公式

而I_xy，I_yz和I_xz则为对O_x与O_y轴、O_y轴与O_z轴和O_z轴与O_x轴的惯性积，分别为

点击浏览公式

方程中需说明的参数表达式如下

气动数据与基于某型机的全机测力风洞试验，所得到的各气动参数精度均满足合格或先进指标的要求。图 2中所涉及的结冰模型如公式(12)所示，该公式基于Bragg等人提出的积冰影响模型构建。

点击浏览公式

其中C_(A)与C_(A)ice分别为结冰前后的气动参数，η为飞机的结冰严重程度参量，K_CA为飞机的结冰因子常数，可以认为是飞机的固有参数。在本文中，为反应结冰的动态过程，将参数η转变为随时间变化的函数，表示为η(t)。η在不同气象条件下的变化规律可采取实验或数值仿真的手段获得。文中所使用的η(t)基于某型飞机冰风洞实验数据拟合而来，具有较高的准确性。

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

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

图选项

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

图 5 150组极值样本散点分布 Fig. 5 The scatter distributions of 150 extreme samples

图选项

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

点击浏览公式

其中α_max为迎角极值，$\bar{v}$_max为相对速度极值；α_c为结冰时的失速迎角，它是襟翼角度δ_f（本文中δ_f（为0）、马赫数Ma与结冰程度因子η的函数；$\bar{v}$_c代表相对失速速度。

从图 5中可以看出相对速度和迎角极值的分布都存在明显的厚尾特性，这种分布形式在低频高危事件（如：地震、海啸、金融风险、飞行事故等）中较为常见。为进一步验证极值样本的分布特性，对飞参极值变量的统计特性进行研究。飞参极值变量的统计分析结果如表 1所示。

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

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

图 6 极值样本盒形图 Fig. 6 The box graph of sample extremes

图选项

（1）基于复杂多因素耦合系统仿真理论与方法，构建了结冰条件下的人-机-环实时仿真系统；构建了基于四元数法的某型飞机六自由度运动方程；考虑结冰条件下的不确定性与随机性因素，提出了提取飞参极值样本的蒙特卡罗仿真法。

（2）对中度结冰条件下的人-机-环系统耦合特性进行了分析。结果表明在中度结冰条件下某型飞机的动力学特性逐步恶化，迎角和法向过载增量较大，易于超出临界值以致飞机失控；同时由于阻力系数的增加和升力系数的减小，飞行速度和飞行高度逐渐降低。如驾驶员在此时介入延迟较大且操纵不当，极易引起飞行参数极值超限从而引发飞行风险；驾驶员需在中度结冰条件下迅速介入并平稳操纵以脱离飞行风险。

（3）对飞参极值样本进行了统计学分析，确定了其具有显著的厚尾特性。本文思路及方法不仅局限于结冰条件下人-机-环仿真计算，也可以用来研究其他具有多因素耦合特性的情形，比如：危险科目下的试飞、复杂外部气流下的飞行、飞机软件或硬件故障下的飞行等等。