裴彬彬， 徐浩军， 薛源,*， 李哲， 刘东亮

1.空军工程大学 航空航天工程学院， 西安 710038 2.北京航空工程技术研究中心， 北京 100076

基于复杂动力学仿真的结冰情形下飞行安全窗构建方法

裴彬彬1， 徐浩军1， 薛源1,*， 李哲1， 刘东亮2

1.空军工程大学 航空航天工程学院， 西安 710038 2.北京航空工程技术研究中心， 北京 100076

目前针对结冰情形下增强驾驶员情景感知的研究比较有限，现有的手段一般为通过评估部分飞行安全关键参数是否超出其极限值来对风险事件是否发生进行预测。建立了驾驶员操纵-飞机本体-积冰影响的动力学模型，通过对单个飞行情形预测时间段内飞行参数风险度的叠加得到该情形下的飞行安全谱，并在此基础上得到该情形下的风险值。基于建立的并行飞行仿真平台，获取飞机在整个操纵范围内的风险拓扑图，即安全窗。分析了飞机在对称结冰情形和非对称结冰情形下飞行安全窗的变化，并对结冰的致灾机理进行了分析。仿真结果表明结冰导致安全飞行范围缩减，对于非对称结冰还会出现安全窗不对称的现象。安全谱的提出可以为事故的演化分析提供一种全面直观的分析方法，安全窗的构建可为飞机遭遇各种不利情形下的驾驶员操纵提供指示，也可为飞机设计人员优化飞机性能提供一定的参考。

飞机结冰； 情景感知； 计算飞行动力学； 非对称结冰； 并行飞行仿真； 安全窗

现代事故致因理论认为，飞行事故的发生往往是人-机-环复杂系统内部因素之间耦合作用的结果[1]。驾驶员在事故的演化过程中扮演了重要的角色，驾驶员在各种情形下的操纵正确与否直接关系到飞行安全。波音公司对全球商用喷气式飞机事故统计结果表明，超过50%的飞行事故主要由机组人员造成[2]，空客公司曾指出大约85%的飞行事件报告中均提到了驾驶员情景感知(Situational Awareness)能力的丧失[3]。提高驾驶员的情景感知能力对于飞行安全至关重要，这也是当今国际航空领域研究的热点问题之一[4-8]。

在众多的研究中，针对结冰情形下增强驾驶员情景感知的研究并不多见。Bragg[9]和Deters[10]等开创性地提出了智能结冰系统(Smart Icing System，SIS)的概念，该系统能够实时地感知积冰的存在及其影响，并计算出结冰状态下飞机的飞行安全边界提供给驾驶员及飞控系统。Gingras等[11-12]联合开发了积冰污染边界保护(Icing Contamination Envelope Protection，ICEPro)系统，该系统为驾驶员提供积冰信息及飞机的边界信息，具备提升驾驶员在结冰情形下情景感知的能力，其在地面模拟器试飞的结果得到了驾驶员的肯定[13]。

国内针对飞机结冰后飞行动力学特性改变对飞行安全影响的研究起步较晚。徐忠达等[14]就冰形对气动参数及操稳特性的影响进行了研究；Dong和Ai[15]研究了结冰后的参数辨识；应思斌和艾剑良[16]主要就结冰条件下的容冰控制进行研究。刘东亮等[17]提出了积冰条件下基于复杂系统动力学仿真的风险评估方法。王明丰等[18]就结冰对飞行动力学的影响进行了分析，张智勇[19]主要对智能防冰系统中的包线保护控制律设计问题进行了研究。由于研究手段等方面的限制，在致灾机理与防护方面，国内尚无系统性的研究。

当前结冰边界保护系统为驾驶员提供的都是实时的飞行状态及其边界信息，如迎角和舵面操纵量及其限制值等。NASA Langley研究中心近期(2013年)的一项调查结果表明[20]，驾驶员更倾向于得到可用滚转角和俯仰角等指示信息，特别是在安全边界改变的情形下，同时，提供这些信息并不会带来驾驶员负荷的显著增加。另一方面，可用滚转角和俯仰角等信息相对于实时的飞行状态参数及其极限值而言，更具预测性。而这正是当前飞机结冰研究所缺乏的。因此，研究结冰情形下飞机俯仰角和滚转角的可用安全范围对于提高驾驶员情景感知能力、保障飞行安全具有重要的作用。

本文提出了基于计算飞行动力学的结冰条件下飞行安全窗的构建方法，将飞机在机动空间中可能飞行轨迹的飞行风险以色彩化的形式呈现出来，揭示人机闭环系统在结冰情形下的安全边界以及不同结冰情形下的致灾机理，为驾驶员的安全操纵提供直观的指示告警。

1 模型的建立

对于结冰情形下飞行安全区域的确定要涉及临界飞行状态，往往伴随着强耦合、非线性等特点，因此采用6自由度全量非线性运动方程更能准确地描述飞机的运动规律。分别对飞机本体动力学模型、驾驶员操纵策略模型、积冰影响模型进行了研究，建立起不同结冰情形下的人-机-环闭环系统模型，为进行动力学仿真打下基础。

1.1 飞机动力学模型

1.1.1 飞机本体模型

飞机本体非线性动力学模型可用式(1)所示的向量形式表示[21-24]

(1)

式中：x为状态向量;u为控制向量；且

x=[Vαβq0q1q2q3pqr

xgygzg]T

(2)

其中:V、α和β分别为飞行速度、飞机迎角和侧滑角；xg、yg和zg为飞机在地面坐标系下的位置。

u=[δthδeδaδr]T

(3)

其中：δth为飞机油门偏度；δe、δa和δr分别升降舵、副翼以及方向舵偏转的大小。

为避免奇点的产生，采用四元数法来建立飞机动力学模型，得到向量f的表达式为

(4)

(5)

(6)

(7)

此外，

(8)

(9)

飞机的姿态角：滚转角φ、俯仰角θ和偏航角ψ可以通过四元数变换求得。

1.1.2 舵机模型

舵机模型可简化为由一阶惯性环节、速率限制器和位置限制器构成[25]，如图1所示。图中：T为一阶惯性环节的时间常数；1/(Ts+1)为舵机的传递函数。

图1 舵机动力学模型示意图
Fig.1 Schematic diagram of dynamics model for actuator

舵机速率限制器以及位置限制器的具体参数可参照文献[26]中的有关参数进行设置。

1.2 驾驶员模型

典型的驾驶员模型结构如图2所示[27]，输入信号是简单的误差信号，即基准的输入量与人机闭环系统实际响应量之间的偏差，虚线框中的4个模块构成了驾驶员主要行为特性的建模，各模块的说明如图2所示。

1) 含观测噪声的显示环节

驾驶员通过仪表盘等接收飞机当前飞行状态时引入的误差。

2) 延迟环节

驾驶员在读取信号、信号传输至大脑以及大脑做出决策的过程存在着一个不可控的最小反应延迟时间，常用e-τs表示，τ的范围一般为0.06～0.30 s。

3) 补偿操纵环节

补偿操纵环节相当于是反馈控制系统中的控制器，构成了驾驶员模型的核心。在指定飞机目标俯仰角、滚转角的情况下，根据文献[28-29]中采用的驾驶员模型，本文在各通道上建立了以下驾驶员补偿操纵策略模型。

俯仰角控制模型为

(10)

式中：Ixx、Iyy、Izz和Ixz为飞机的转动惯量；Kθ1和Kθ2为常数值；其他中间变量的计算参照文献[28-29]中所示。滚转角控制模型为

(11)

式中：φc为指令滚转角;Kpa、KDa和KIa为常数值，但在不同的φc区间段取值不同。

发动机油门控制受飞行速度影响，驾驶员操纵使得飞行速度尽量保持不变。

补偿操纵环节中驾驶员控制模型的参变量直接影响最终的安全窗范围，本文在进行参数的选取时，以尽量能够使飞机按照给定操纵指令飞行为准绳。

4) 神经肌肉动力学模型

当肌肉接受大脑发出的指令时，由于自身的惯性、黏性和收缩而呈现神经肌肉的时延特性，用一阶惯性环节1/(1+TNs)表示；此外驾驶员根据经验可采取提前的操纵量来补偿这种滞后而呈现导前特性，可用一阶微分环节1+TLs来表示[30]。TN和TL的值均约为0.1～0.2 s。另外，驾驶员在操纵过程中会出现自发性的手抖现象，同样会引入噪声信号。

图2 驾驶员模型结构示意图
Fig.2 Schematic diagram of structure of pilot model

观测噪声及操纵过程中的噪声信号对动力学仿真的结果影响并不大，为适当地简化计算，文中予以忽略。

1.3 结冰影响模型

飞行中根据气象条件来对冰形进行预测并估算对飞机空气动力特性的影响是对飞机飞行动力学特性改变进行计算的基础，也是最终实现飞行风险操纵区域划分的基础，国内外相关的研究也比较多。本文的研究重点在于在已经预测出冰形对飞机空气动力特性影响的前提下，计算飞机安全操纵范围，从而为驾驶员操纵提供有价值的信息。

1.3.1 对称结冰影响模型

利用6自由度(6-DOF)运动方程对飞机在积冰后的动力学特性进行研究，积冰前后，运动方程的形式并没有变，只是由于积冰影响了飞机的气动特性，使得方程中的力和力矩发生了变化。

文中采用Bragg等[31]提出的积冰影响模型对结冰后气动力进行建模。按照其理论，积冰前后气动参数的关系为

C(A)iced=(1+ηkC(A))C(A)

(12)

式中：C(A)与C(A)iced分别为任意的结冰前后飞机的性能、稳定性与控制参数或其导数；η为飞机结冰程度参数，取决于飞机的固有参数和结冰条件，其值越大，积冰后气动参数的变化也就越大，表明飞机积冰情况越严重，根据η随云层参数的变化曲线，积冰程度参数的变化范围大致在0～0.3之间[31]；kC(A)为结冰对飞机气动参数的影响参数，对于某架特定的飞机来说为一常数，其值可以通过数值仿真计算或飞行试验获得。

1.3.2 非对称结冰影响模型

在实际飞行过程中，考虑到冰形形成的随机性及除冰系统工作后冰形脱落的随机性，并不能保证左右机翼气动特性的一致。特别是当机翼除冰系统一侧发生故障，导致左右机翼出现升力、阻力差，并由此产生滚转及偏航力矩时，极易引发飞行事故。这种结冰情形不同于均匀结冰，因此，有必要单独构建出非对称结冰时的模型。

Lampton和Valasek[32-33]指出左右机翼的非对称结冰主要会导致左右机翼非对称升力与阻力的产生，从而产生附加的偏航与滚转力矩，基于此，可以建立非对称结冰情形下的积冰影响模型。文中假设右侧机翼的除冰系统出现故障不能正常除冰，则左右机翼的升力及阻力差值为

(13)

式中：CL和CLice分别为飞机结冰前后的升力系数；CD和CDice分别为结冰前后的阻力系数。此时，两侧机翼的升力差产生了正的滚转力矩，即使飞机向右滚转的力矩；阻力差产生了正的偏航力矩，即使飞机向右偏航的力矩。附加的滚转和偏航力矩可表示为

(14)

式中：dmgc为平均几何弦长位置到飞机中心线的距离；Q为动压；Sw为机翼面积。将式(13)和式(14)代入飞机6-DOF运动方程中的力和力矩项，通过仿真计算，即可得到左右机翼非对称结冰时的飞机状态变化。

1.3.3 结冰对其他参数的影响

结冰不仅会导致飞机气动参数的改变，同时还会改变一些飞行安全有关的参数的限制值。如结冰会引起失速迎角降低[34]、最小平飞速度增大[35]等。

这些与安全相关的参数改变可以事先通过计算以数据库的形式存储在计算机中，对驾驶员安全飞行边界进行计算时再进行调用。

需要指出的是上文所述的结冰影响模型是属于目前国内外比较常用的结冰后气动参数估算模型，由于其包含了一定的物理意义，既考虑了飞机特性又考虑了结冰条件的影响，具有一定的通用性，该方法适合作为结冰后气动参数尚未获取的情况下对飞行动力学特性进行初步的分析。若要较为准确地预测出结冰后的飞行动力学特性，则需要借助于高精度数值模型、风洞试验或真实试飞等手段来对特定的飞机预先获取结冰后准确的气动力方可进行。

2 单个飞行情形下飞行安全谱及风险量化

2.1 飞行安全谱计算

飞行风险的发生往往伴随着飞行参数的异常变化，通过对单个飞行情形下的飞行参数数据变化情况进行分析，可以预判出飞行风险事件的发生与否。本节将飞行参数数据在不同区间上的值以色彩化的形式呈现，通过对单个飞行情形下不同飞行参数的色谱信号叠加，来得到该飞行情形下的飞行安全谱。安全谱的获得可以很直观地看出单个飞行情形下飞行风险的演化过程。

当前对于飞行安全参数的限制的描述往往是确定性的，如某型飞机的手册中规定飞机的最大允许迎角α在马赫数Ma<0.55时，不得超过16°。其潜在的含义意味着迎角小于16° 即是安全的，大于16° 将会导致飞行事故的发生。Burdun[36]认为实际中人们对于这种限制的认识是带有模糊性质的，如迎角达到15° 时同样也是非常危险的状态。为此，他们提出了用色彩来对飞行风险进行表示的方法，但他们所提出的方法仅能用来表示飞行参数处于风险状态，而不能表示此时飞行参数是处于正的风险区间还是负的风险区间。为此，本文提出了一种考虑飞行参数风险区间正负性的表示方法，该方法可以更加合理地表示飞行参数变化及飞行风险的变化情况。

表1 与安全相关的飞行参数色彩化区间Table 1 Colored interval of safety related flight data

表格中舵面操纵量浅灰/深灰色对应的值为舵面的极限偏转角度。对于舵面操纵量δe、δa和δr而言，舵面偏转至极限位置飞行事故不一定会发生，其风险度设定为与浅红/深红一致，这里只是为了区分舵面操纵是否饱和，这一点与其他状态参数不同。

通过人-机-环系统仿真，得到每一个关键的飞行参数在预测时间段内的变化情况，根据表1中飞行参数的色彩化区间分布得到这些关键飞行参数的色谱图Ci(t)，图3所示为根据迎角在预测时间段内的变化情况得到的迎角安全谱。

飞行安全是典型的木桶理论，任何一个关键飞行参数的异常都可能会导致飞行风险事件的发生，故最终的飞行安全谱上每一时刻对应的风险色与所有飞行参数在该时刻处在最危险状态的飞行参数的风险色相同。以某型运输类飞机在平飞状态下以指定的航迹俯仰角(μc=8°)与滚转角(φc=35°)协调爬升转弯为例，计算出该飞行情形下各安全相关的飞行参数及整个飞行情形的飞行安全谱如图4所示。其中，每个飞行参数对应的安全谱为该飞行参数数据在预测时间段内的色谱，最下面一行表示该飞行情形整个时间段内的风险变化情况，定义为该飞行情形总的飞行安全谱，其在每一时刻的风险值与该时刻处在最危险状态的与安全相关的飞行参数风险值相同。

图3 迎角变化曲线及相应的风险色谱图
Fig.3 Time history and safety spectra for angle of
attack (AOA)

图4 单个飞行情形的飞行安全谱(μ=8°, φ=35°)
Fig.4 Safety spectra for single flight condition (μ=8°, φ=35°)

由于浅灰/深灰、浅红/深红、浅黄/深黄只是表示风险区间的正负，其代表的风险值是相同的，故在计算总的飞行安全谱时，浅灰/深灰同用黑色表示，浅红/深红同用红色表示，浅黄/深黄同用黄色表示。

由图4的计算结果可知，在驾驶员按照目标飞行指令操纵时，整个飞行阶段操纵舵面δe、δa和δr均保持在一个较低的风险水平，飞行风险偏高的阶段，主要体现在操纵初期的迎角及法向过载值略微偏高；操纵后期爬升率以及滚转角度维持在一个比较高的水平，但总体来说飞行仍是安全的，但除非有必要驾驶员最好能够将滚转角与高度变化率维持在一个较低的值。

2.2 飞行风险的量化

图4中由各种颜色构成的飞行安全谱，代表了整个时间段内的风险度演变过程。然而，不同飞行情形下的飞行安全谱都是不同的，如何根据安全谱中所含的风险信息进行不同飞行情形风险度的对比，是本章的研究内容。提出了单个飞行情形风险量化的方法：采用对各种风险色所占百分比赋权相加的方法来得到该情形下飞行风险值。

假定风险色黑、红、黄、绿分别占整个预测时间段的百分比为Pk、Pr、Py和Pg，将每个风险色代表的风险值定义为Vk、Vr、Vy和Vg，那么在整个预测时间段内的风险值可表示为

R=PkVk+PrVr+PyVy+PgVg

(15)

考虑到安全谱中一旦有黑色区域出现表示至少有一个飞行参数的值超出了安全飞行的最大允许值，即认为飞行事故发生，这是飞行中要极力避免的，可将黑色代表的风险值Vk取一个较大的值，这样一旦出现飞行事故，得到的飞行风险明显要大于其他飞行情形。文中各风险色代表的风险值分别为：Vk=30、Vr=4、Vy=2、Vg=1。以图4所示的飞行情形为例，得到该飞行情形的风险值为1.982。

3 飞行安全窗的构建

安全谱的获得可以很直观地看出飞机以某个目标航向指令飞行的风险演变过程。而单个飞行情形下风险的量化则是构建飞机在整个操纵空间的飞行风险拓扑云图的基础。以指令滚转角与指令航向角构成的操纵空间为计算区域，将计算区域离散成许多计算单元，每个单元对应一个飞行情形，可以通过前述的方法计算出该飞行情形的风险值，进而可得到整个操纵空间风险场的拓扑云图。显然，每个计算单元之间是独立的，风险值的计算过程相互之间并不干扰，因此可以通过并行计算的方式来加快计算进程。

对整个飞行操纵范围的风险度进行计算需耗费一定的时长，采用并行仿真手段可以大大缩减计算时间，提高预测的时效性。本文采用的是基于MATLAB/Simulink平台的并行飞行仿真方法，来解决大量蒙特卡罗飞行仿真的耗时性问题，为计算飞行动力学的实施提供了便利。其实施过程如下：

1) 建立飞机本体运动学模型，并计算在设定初始飞行状态下的配平参数。

2) 在Simulink环境下搭建人-机-环复杂系统飞行动力学模型仿真平台，仿真初始条件设定为配平值。

3) 将Simulink的仿真模式设置为在Rapid Accelerator模式下运行，并指定可调变量，如φc和μc等。

4) 建立该模型的Rapid Accelerator目标对象，初始化可调变量构成的结构体。

5) 在并行计算命令parfor下运行仿真模型，得到在整个操纵空间上每个计算单元的安全谱。

6) 依据式(12)计算出飞机在整个操纵范围上每个计算单元对应的风险值，并在此基础上，做出飞机操纵范围的风险拓扑云图。

需要指出的是，飞行事故的风险值明显要大于一般飞行情形。例如在干净构型下，当指令滚转角φc=45°，指令航迹俯仰角μc=16° 时，通过上述方法得到的该飞行情形的风险值为18.447。过高的风险值会导致在建立风险拓扑云图时，风险值较低的飞行情形的风险色几乎为同一颜色，在拓扑图中区别不明显。为此，将所有飞行风险值大于某个限定值的飞行情形的飞行风险设定为该限定值，这样使得风险色的分离度更加明显。文中将该限定值取为4.5。根据2.2节中风险值的设定可知，只要总的风险值R>4，即意味着预测时间段内出现了参数超限的情形，将限定值定为4.5既保证了对事故状态点的区分，也确保了风险拓扑图中低风险状态点的风险色的分离度。

以飞机在H0=2 000 m，V0=120 m/s水平匀速直线飞行的初始飞行状态下来计算飞机的安全飞行范围。计算在曙光A840r-G服务器上实施，硬件为64核AMD Opteron 6276 2.3 GHz处理器，系统运行内存为128 GB，仿真运行环境为MATLAB 2014a，通过测试，得到并行和串行两种飞行仿真的计算信息对比如表2所示。

安全窗的计算结果如图5所示，图中不同的颜色代表了不同的风险度，图右边的色例为不同的风险值R对应的风险色。

表2并行与串行方法的计算条件与计算时间

Table2Computationconditionandtimeofparallelandserialmethods

ComputationmodeCalculationrangeNodenumberComputa⁃tiontime/sSerialcomputationμc∈[-6∶2∶18]ϕc∈[-55∶5∶55]299293．74Serialcomputationμc∈[-6∶0．5∶18]ϕc∈[-55∶2∶55]27442866．64Parallelcomputa⁃tion32workersμc∈[-6∶2∶18]ϕc∈[-55∶5∶55]29941．49Parallelcomputa⁃tion32workersμc∈[-6∶0．5∶18]ϕc∈[-55∶2∶55]2744232．79

图5 干净外形下的飞行安全窗(初始状态H=2 000 m，V=120 m/s)
Fig.5 Flight safety windows for clean configuration
(trimmed at H=2 000 m，V=120 m/s)

4 案例分析

结冰会导致飞机稳定性和操纵性变差，在飞机设计初期，借助于风洞试验或者流场仿真等手段，获得在不同结冰情形下飞机空气动力学特性改变的前提下，通过前文描述的方法便可获取飞行安全操纵范围，进而为结冰条件下的边界保护系统、飞行控制律、驾驶员情景感知增强系统等的设计提供一定的指导意义。本节以飞机遭遇机翼对称结冰和非对称结冰两种飞行情形为例进行研究。飞行初始条件均设定为H0=2 000 m，V0=120 m/s，飞机保持平飞。

4.1 对称结冰情形下的安全窗

4.1.1 安全窗的计算

假定结冰严重程度参数η=0.1，计算范围φc∈[-55∶2∶55]，μc∈[-6∶0.5∶18]，可得计算节点数为46×49=2 254个，机翼对称结冰时飞机的安全窗如图6所示。

图6 对称结冰情形下的飞行安全窗(初始状态H=2 000 m，V=120 m/s)
Fig.6 Flight safety windows for symmetric icing condition (trimmed at H=2 000 m，V=120 m/s)

从仿真结果可以很明显地看出，结冰导致了飞机安全操纵范围缩减。当飞机处于干净构型时，纵向上，飞机最大正向航迹俯仰角大于18°，最大滚转角约为53.5°。当飞机表面积聚冰形时，飞机在纵向和横向上的安全操纵范围分别缩减至小于12°和小于42°。

4.1.2 致灾机理分析

为研究对称结冰情形下飞行事故发生的致灾机理，在图6中选取4个典型状态点((a)～(d))来进行研究。这4个状态点的参数及其对应飞行

安全谱如图7所示。

从图7所示的仿真结果来看，发生对称结冰时，飞机不管是仅目标航迹俯仰角设定过高(如图7(a) 所示)，还是仅目标滚转角设定过高(如图7(c)所示)，或者是目标航迹俯仰角与目标滚转角同时过高(如图7(b)和图7(d)所示)，往往伴随着：迎角超过其右边界，即大于失速迎角；飞行速度超过其左边界，即小于最小飞行速度；高度变化率超过其左边界，即飞机下降速度过快。即使目标航迹俯仰角设定的是正值，飞机最终也会陷入高度下降率变化过快的情形。

从事故演化过程来看，往往是由于飞机迎角超过结冰后的失速迎角，同时由于结冰带来的气动性能的恶化，飞机的推力不足以维持安全飞行的需要，导致飞行速度减小到其极限值，随后飞机开始急速下降，此时很难避免事故的发生。因此，为确保对称结冰情形下的飞行安全，首先要确保飞机处在安全的飞行范围，驾驶员或者飞控系统还应当重点关注飞行过程中迎角、速度和高度变化率这3个安全关键参数的变化情况。

图7 对称结冰情形下典型事故状态点的飞行安全谱
Fig.7 Flight safety spectrum of typical accident state points under symmetric icing conditions

4.2 非对称结冰情形下的安全窗

4.2.1 安全窗的计算及分析

假定飞机在上述初始飞行条件及结冰气象条件下，右侧机翼除冰系统出现故障时，根据1.3.2节中所建立的非对称结冰模型，通过文中所提出的飞行安全窗的计算理论，得到飞机在该情形下飞行风险的安全窗如图8所示。

图8 非对称结冰情形下的飞行安全窗——右侧除冰系统故障(初始状态H=2 000 m，V=120 m/s)
Fig.8 Flight safety windows for asymmetry icing conditions—right wing half fatigue (trimmed at H=2 000 m，V=120 m/s)

从计算的结果来看，非对称结冰情形下，结冰不但导致安全范围的缩减，同时当结冰的不对称性超过一定程度时，飞机的安全飞行范围还出现了明显的不对称现象：对于绿色安全飞行区域，当飞机爬升时，飞机偏向除冰系统故障那一侧的安全飞行范围相对大一些，而当飞机下降时，飞机偏向除冰系统正常工作的一侧进行操纵安全性相对大一些；对于安全窗整体而言，飞机向除冰系统正常一侧滚转时的安全飞行范围比飞机向除冰系统故障一侧滚转时的安全飞行范围要大些。下面分析造成结冰安全操纵范围不对称现象的原因。

1) 绿色飞行区域不对称原因分析。

为分析非对称结冰情形下，绿色飞行区域不对称现象的原因，在图8中选择4个状态点进行分析，它们的位置如图8中蓝色三角点所示，其对应的安全谱如图9所示。

从图9中对比分析的结果来看，造成机翼非对称结冰时的非对称安全飞行范围的主要因素是由于飞机向左或向右滚转时产生的高度变化率不同导致的：当飞机滚转爬升时，向左滚转爬升时的高度变化率(正值)要大于向右滚转爬升时的高度变化率，导致向左滚转爬升时的高度变化率更容易超限，因此飞机向右滚转爬升时的安全飞行范围要大些；而当飞机滚转下滑时，向左滚转下滑时的高度变化率(负值)要慢于向右滚转下滑时的高度变化率，导致向右滚转下滑时的高度变化率更容易超限，因此飞机向左滚转下滑时的安全性相对而言更高些。

2) 安全窗整体不对称原因分析。

飞机发生非对称结冰时，除冰系统工作正常一侧的可用飞行范围明显大于除冰系统故障的一侧，为分析原因，在图8中选择4个状态点进行对比，它们的位置如图8中紫色三角点所示。这4个状态点对应的飞行安全谱如图10所示。

图9 非对称结冰情形下典型状态点对应的飞行安全谱(φ=±25°)
Fig.9 Flight safety spectrum of typical state point under asymmetry icing conditions (φ=±25°)

图10(a)和图10(c)中，仿真时长均为60 s，而图10(b)和图10(d)中仿真均停止于滚转角大于150°(滚转角达到不可逆转的值)。对比分析这4种飞行状态的飞行安全谱可知，对于非对称结冰情形而言，飞机向着除冰系统故障的一侧滚转机动时，当目标滚转角设定较大时飞机很难达到并保持该目标滚转角，极易引发滚转角迅速发散的情形，进而导致飞行事故的发生。当飞机向着除冰系统工作正常的一侧滚转机动时，对于同样大小的目标滚转角，飞机仍可达到并保持该飞行状态进行飞行，各安全关键参数变化稳定。因此，对于非对称结冰情形而言，尤其要关注飞机向除冰系统故障一侧滚转角的变化情况。

图10 非对称结冰情形下典型状态点对应的飞行安全谱(φ=±42°)
Fig.10 Flight safety spectrum of typical state point under asymmetry icing conditions (φ=±42°)

需要指出的是，虽然飞机向左滚转的飞行范围与干净构型下的情况相差不大，但相比于干净构型的飞机，驾驶员模型中增益值相对较小，也就是说驾驶员在实际操作中要采用比较柔和的操纵方式才能保证飞行安全。

4.2.2 致灾机理分析

为研究非对称结冰情形下飞行事故发生的致灾机理，选取4个典型的飞行事故状态点来进行研究，它们的位置如图8中白点所示。这4个状态点对应的飞行安全谱如图11所示。

图11 非对称结冰情形下典型事故状态点的飞行安全谱
Fig.11 Flight safety spectrum of typical accident state points under asymmetry icing conditions

图11(a)和图11(b)中，指令滚转角设定为向除冰系统正常一侧以较大角度偏转，仿真结果表明滚转角出现了大幅发散振荡，直至超出边界；副翼正偏(即右副翼下偏)的操纵效能不足，使得驾驶员很难使滚转角维持在一个较大的角度；同时还伴随着迎角的超限。图11(c)表明当目标航迹俯仰角设定过大时，飞机的横向滚转会向着除冰系统故障的一侧而无法纠正，同时迎角也会一直保持在一个较危险的区间直至飞行事故的发生。对于指令滚转角设定为向除冰系统故障一侧以较大角度偏转时，如图10(b)、图10(d)和图11(d)中所示(对应不同的航迹俯仰角)，滚转角很快便发散至不可挽回的地步。

仿真结果表明，飞机出现非对称结冰时，引起飞行事故的原因主要是滚转角超限，而且向着除冰系统故障的一侧滚转角很容易超限，即使是目标航迹滚转角设定为向除冰系统工作正常的一侧偏转，飞机也会由于副翼正向偏转操纵效能的不足而最终演化为向着除冰系统故障的一侧异常偏转。此外，迎角也很容易在事故演化过程中达到其极限值。而此时飞机飞行速度、高度变化率基本处在安全的范围之内，只有在仿真的末期，才会出现超限的可能，这与对称结冰情形下的致灾机理有所差异。

因此，在非对称结冰情形下，驾驶员或者飞行控制系统还应当首要关注飞行过程中的滚转角和迎角这两个安全关键参数的变化情况。

5 结 论

基于人-机-环复杂系统动力学仿真，研究了飞行安全窗的构建方法。建立了蒙特卡罗并行飞行仿真平台，计算了不同结冰情形下的飞行安全窗口，揭示了不同结冰情形下，由于驾驶员可能的操纵失误带来的飞行事故发生机理，并给出驾驶员应对策略建议。

结冰对于不同飞机的气动性能影响必然是不同的，文中只是以某典型构型的运输类飞机为例对安全窗的构建方法及结冰致灾机理进行研究分析。对于不同的飞机其在结冰情形下安全窗的变化情况必然会受到结冰影响模型准确性以及飞机固有特性的影响。本文侧重于对风险谱的建立及风险的量化以及安全窗的建立方法进行研究，所提出的理论方法及研究结论对于结冰条件下致灾与防护方面的研究具有一定的借鉴意义和价值。

计算飞机在不同结冰情形下的风险安全窗，可以以数据库的形式存储在机载计算机中。当飞机处于不同的飞行状态下，遭遇不同程度的结冰时，可以通过插值计算方式得到飞机在该状态下的安全窗，从而为驾驶员提供指示，提升结冰情形下驾驶员的情景感知能力。

驾驶员的操纵具有随机性，在安全的飞行范围内也存在着驾驶员的粗暴操纵导致飞行事故发生的可能。文中在计算时并没有考虑这一点。在对驾驶员模型进行建模时，考虑的只是尽可能地使飞机按照指定的飞行指令来飞行。因此得到的安全飞行范围只是一个参考范围，驾驶员在进行操纵时根据提供的安全窗可以对自己的操纵裕度产生直观的感受，从而保证飞行安全。

此外，文中所提出的方法还可以用在飞行控制系统设计(控制律参数优化)、验证控制律、飞机性能计算、飞机遭遇各种不利情形时(如单发失效、多不利因素情形)的操纵范围。

[1] BURDUN I Y. A method for accident reconstruction and neighborhood analysis using an autonomous situational model of flight and flight recorder data[C]//Advances in Aviation Safety Conference and Exposition. Warrendale, PA: SAE International, 1999: 1-13.

[2] Boeing. Airplanes statistical summary of commercial jet airplane accidents: Worldwide operations 1959-2005[R]. Chicago, IL: The Boeing Company, 2006: 1-25.

[3] Airbus Customer Services. Human performance: Enhancing situational awareness: FLT_OPS-HUM-PERF-SEQ 06-REV 01[R]. Blagnac Cedex: Airbus, 2007.

[4] BURDUN I Y. The intelligent situation awareness and forecasting environment (The S.A.F.E. Concept) a case study[C]//Advances in Aviation Safety Conference and Exposition. Warrendale, PA: SAE International, 1998: 1-15.

[5] ZOLGHADRI A. Early warning and prediction of flight parameter abnormalities for improved system safety assessment[J]. Reliability Engineering and System Safety, 2002, 76(1): 19-27.

[6] ROUWHORST W F J A, MARSMAN A P L A. A piloted investigation of an integrated situation awareness system (ISAS): AIAA-2006-6723[R]. Reston: AIAA, 2006.

[7] BORST C, SJER F A, MULDER M, et al. Ecological approach to support pilot terrain awareness after total engine failure[J]. Journal of Aircraft, 2008, 45(1): 159-171.

[8] BORST C, GROOTENDORST F H, BROUWER D I K, et al. Design and evaluation of a safety augmentation system for aircraft[J]. Journal of Aircraft, 2014, 51(1): 12-22.

[9] BRAGG M B, BASAR T, PERKINS W R, et al. Smart icing systems for aircraft icing safety[C]//40th AIAA Aerospace Sciences Meeting & Exhibit. Reston: AIAA, 2002.

[10] DETERS R, DIMOCK G A, SELIG M S. Icing encounter flight simulator with an integrated smart icing system: AIAA-2002-4599[R]. Reston: AIAA, 2002.

[11] GINGRAS D R, BARNHART B, RANAUDO R, et al. Envelope protection for in-flight ice contamination[C]//47th AIAA Aerospace Sciences Meeting including The New Horizons Forum and Aerospace Exposition. Reston: AIAA, 2009.

[12] GINGRAS D R, BARNHART B, RANAUDO R, et al. Development and implementation of a model-driven envelope protection system for in-flight ice contamination[C]//AIAA Guidance, Navigation, and Control Conference. Reston: AIAA, 2010.

[13] RANAUDO R, MARTOS B, NORTON B, et al. Piloted simulation to evaluate the utility of a real time envelope protection system for mitigating in-flight icing hazards[C]//AIAA Atmospheric and Space Environments Conference. Reston: AIAA, 2010.

[14] 徐忠达, 苏媛, 曹义华.平尾积冰对飞机纵向气动参数的影响[J]. 航空学报, 2013, 34(7): 1563-1571.

XU Z D, SU Y, CAO Y H. Effects of tailplane icing on aircraft longitudinal aerodynamic parameters[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(7): 1563-1571 (in Chinese).

[15] DONG Y, AI J. Research on inflight parameter identification and icing location[J]. Aerospace Science and Technology, 2013, 29(1): 305-312.

[16] 应思斌， 艾剑良. 飞机结冰包线保护对开环飞行性能影响与仿真[J]. 系统仿真学报, 2010, 22(10): 2273-2301.

YING S B, AI J L. Simulation of aircraft flight envelope protect in icing encounters effects on open loop dynamic [J]. Journal of System Simulation, 2010, 22(10): 2273-2301 (in Chinese).

[17] 刘东亮， 徐浩军， 李嘉林， 等. 飞行结冰后复杂系统动力学仿真与风险评估[J]. 系统仿真学报, 2011, 23(4): 643-647.

LIU D L, XU H J, LI J L, et al. Dynamic simulation study of stalling in wing icing conditions and risk evaluation[J]. Journal of System Simulation, 2011, 23(4): 643-647 (in Chinese).

[18] 王明丰, 王立新, 黄成涛. 积冰对飞机纵向操稳特性的量化影响[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 Beijing University of Aeronautics and Astronautics, 2008, 34(5): 592-595 (in Chines).

[19] 张智勇. 结冰飞行动力学特性与包线保护控制律研究[D]. 南京： 南京航空航天大学, 2006.

ZHANG Z Y. Research on iced aircraft flight dynamics characteristics and envelope protection control law[D]. Nanjing: Nanjing University of Aeronautics and Astronautics， 2006 (in Chinese).

[20] TRUJILLO A, GREGORY I. Pilot preferences on displayed aircraft control variables[J]. Lecture Notes in Computer Science, 2013, 8020(1): 203-211.

[21] SIBILSKI K, LASEK M, LADYZYNSKA-KOZDRAS E, et al. Aircraft climbing flight dynamics with simulated ice accretion: AIAA-2004-4948[R]. Reston: AIAA, 2004.

[22] SONNEVELDT L. Nonlinear F-16 model description[R]. The Netherlands: Delft University of Technology, 2010.

[23] COOK M V. Flight dynamics principles[M]. Burlington, MA: Butterworth-Heinemann, 2007: 66-95.

[24] 吴森堂. 飞行控制系统[M]. 北京: 北京航空航天大学出版社, 2013: 54-64.

WU S T. Flight control system[M]. Beijing: Beihang University Press, 2013: 54-64 (in Chinese).

[25] HUESCHEN R M. Development of the transport class model (TCM) aircraft simulation from a sub-scale generic tra nsport model (GTM) simulation: NASA/TM-2011-217169[R]. Washington, D.C.: NASA, 2011.

[26] HANKE C R, NORDWALL D R. The simulation of a jumbo jet transport aircraft. Volume 2: Modeling data/detail: NASA-CR-114494[R]. Washington, D.C.: NASA, 1970.

[27] DOGAN A, KAEWCHAY K. Probabilistic human pilot approach: Application to microburst escape maneuver[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(2): 357-369.

[28] GAO Z, GU H. Simulation of microburst escape with probabilistic pilot model[J]. Lecture Notes in Computer Science, 2011, 7030(1): 663-670.

[29] DOGAN A, KABAMBA P T. Escaping microburst with turbulence altitude, dive, and pitch guidance strategies[J]. Journal of Aircraft, 2000, 37(3): 417-426.

[30] 高金源， 李陆豫， 冯亚昌. 飞机飞行品质[M]. 北京: 国防工业出版社, 2003: 135-140.

GAO J Y, LI L Y, FENG Y C. Aircraft handling quality[M]. Beijing: National Defence Industry Press, 2003: 135-140 (in Chinese).

[31] BRAGG M, HUTCHISON T, MERRET J . Effect of ice accretion on aircraft flight dynamics[C]//38th AIAA Aerospace Sciences Meeting & Exhibit. Reston: AIAA, 2000.

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

[33] LAMPTON A, VALASEK J. Prediction of icing effects on the coupled dynamic response of light airplanes[C]//AIAA Atmospheric Flight Mechanics Conference and Exhibit. Reston: AIAA: 2007.

[34] SHARMA V, VOULGARIS P G, FRAZZOLI E. Aircraft autopilot analysis and envelope protection for operation under icing conditions[J]. Journal of Guidance, Control, and Dynamics, 2004, 27(3): 454-465.

[35] Federal Aviation Administration. Pilot guide flight in icing conditions: AC 91-74A[R]. Washington, D.C.: FAA, 2007.

[36] BURDUN I Y. Automated planning, exploration and mapping of complex operational domains of flight using multifactor situational trees[J]. SAE International Journal of Aerospace, 2011, 4(2): 1149-1175.

(责任编辑： 李明敏)

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

Developmentofflightsafetywindowinicingconditionsbasedon

complexdynamicssimulationPEIBinbin1,XUHaojun1,XUEYuan1,*,LIZhe1,LIUDongliang2

1.AeronauticsandAstronauticsEngineeringCollege,AirForceEngineeringUniversity,Xi’an710038,China2.BeijingAeronauticalEngineeringTechnologyResearchCenter,Beijing100076,China

Studiesonenhancingsituationalawarenessofpilotsinicingconditionsarerelativelyrare.Currentmethodsnormallypredictoccurrenceofaccidentsbyestimatingwhetherthesafetyrelatedparametersexceedtheirlimitations.Complexdynamicsofthepilot-vehicle-icingeffectmodelisestablished.Safetyspectrumofasingleflightconditionisobtainedbycomparingtheriskdegreeoftheflightdatathroughthepredicttimeinterval,andthecoloredriskvaluefortheflightconditionisfurtheracquired.Safetywindowforflightsafetyinthewholemanipulationcanbecalculatedusingparallelflightsimulationplatform.Thesafetywindowsofsymmetricandasymmetricicingconditionsareresearched,andthedisaster-causingmechanismisanalyzed.Simulationresultsshowthaticewillleadtoshrinkofthesafetyflightscope,andtheasymmetricalicewillleadtotheunevensafetywindow.Theproposedmethodcanprovidetheoreticalsupportforenhancementofsituationalawarenessindifferentadverseevents,andengineeringtoolforoptimizingtheaircraftperformanceforaircraftdesigners.

aircrafticing;situationalawareness;computationalflightdynamics;asymmetricicing;parallelflightsimulation;safetywindow

2016-08-23；Revised2016-09-20；Accepted2016-10-25；Publishedonline2016-11-031633

s：NationalBasicResearchProgramofChina(2015CB755800);NationalNaturalScienceFoundationofChina(61374145,U1333131)

.E-mailszxy1986@163.com

2016-08-23；退修日期2016-09-20；录用日期2016-10-25； < class="emphasis_bold">网络出版时间

时间：2016-11-031633

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

国家“973”计划 (2015CB755800)； 国家自然科学基金 (61374145，U1333131)

.E-mailszxy1986@163.com

裴彬彬， 徐浩军， 薛源， 等． 基于复杂动力学仿真的结冰情形下飞行安全窗构建方法J． 航空学报,2017,38(2):520695.PEIBB,XUHJ,XUEY,etal.DevelopmentofflightsafetywindowinicingconditionsbasedoncomplexdynamicssimulationJ.ActaAeronauticaetAstronauticaSinica,2017,38(2):520695.

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

10.7527/S1000-6893.2016.0274

V212

A

1000-6893(2017)02-520695-14