结冰飞机非线性稳定域确定及安全操纵方法

周驰， 李颖晖， 郑无计， 武朋玮， 董泽洪

(空军工程大学航空工程学院, 西安 710038)

结冰不仅会恶化飞机的动力学特性，也会严重影响驾驶员正常操纵，因而结冰飞行事故持续发生[1]。驾驶员在结冰情况下仍继续保持常规操纵，飞行状态极易超出安全边界[2]。NASA对1976—1994年间的飞行事故进行了数据统计，其中由结冰引发的飞行事故就有16起，共造成139人死亡[3]。结冰不仅使飞机升力减小，阻力增大，失速迎角也会对应减小，这导致飞机更易失速。1994年，美国鹰航公司的一架ATR72-212飞机因机翼严重结冰在印第安纳州坠毁，当时飞行迎角仅为5°，远远低于飞控系统给出的限制值18.1°，但事故却仍因飞机失速而发生[4]。

对于刚体飞机而言，一旦控制结构确定下来，整个飞控系统也被固定。同时飞行指令也存在一个限制，只要机组下达的飞行指令超出这个限制，飞机将很难回到原有的平衡状态。当飞机的动力学特性发生变化，就可能出现抖动甚至失速。文献[5]指出通过分岔分析方法能够清晰地对环境影响下的飞机动力学状态变化进行观测，此时驾驶员可通过操纵将飞行状态调整到新的稳定平衡状态，这也是一种常用的边界恢复方法。分岔分析方法目前正被广泛应用于研究飞机的动力学特性变化，并用于指导飞行操纵[6-7]。Engelbrecht等[8]指出飞行系统的稳定性与分岔点有关。操纵过度，分岔就会出现，系统就可能失稳。Harry和Jean[9]利用分岔分析方法对飞机动力学特性进行了研究，但其只是研究了对飞行指令的限制，虽然提出了当飞机失稳后将飞行状态调整到某一稳定平衡点(Stable Equilibrium Point，SEP)的策略，但是该状态点附近稳定边界的范围却没有定义，如果找到的稳定平衡点附近的稳定边界较小，即使将飞行状态调整到该稳定状态，一个轻微的扰动都能使飞行状态重新失稳。为提升飞行安全，应该在此基础上找到可靠的边界确定方法。

近年来，国内外对安全边界的确定方法已经进行了大量研究。传统的边界确定方法是通过数值计算或飞行试验，给出迎角、侧滑角、俯仰角、飞行速度、过载等飞行参量的限制值。但由于传统方法无法将环境因素考虑在内，因此就不再适用于结冰条件安全边界确定。近年来，吸引域(Region of Attraction，ROA)方法[10]被用于安全边界的确定，但是ROA方法主要是基于Lyapunov能量函数[11]确定边界，精确性依赖于所选择的Lyapunov函数，因此保守性较强。Weekly和Tinka[12]提出用可达集方法确定飞行安全边界，并将最大可控不变集作为安全边界，但可达集方法基于最优控制，在不同控制下可能确定的边界也不同，保守性较强，同时该方法不是基于稳定平衡状态确定边界，因此在指导飞行操纵上存在一定难度。

本文提出了利用微分流形理论对结冰条件下的飞行系统进行非线性稳定域确定，该方法在稳定域的确定上主要是基于稳定平衡点。稳定域的确定是通过计算稳定平衡点稳定边界上的不稳定平衡点(Unstable Equilibrium Point，UEP)的稳定流形的并集得到，并将其作为飞行系统的安全边界[13-14]。只要将系统的飞行状态控制在稳定域内，飞行系统就不会有失稳风险，而在稳定域外的点则会振荡发散[15]。本文首先利用分岔分析方法对不同结冰程度条件下的飞机的动力学特性变化进行观测，给出了不同条件下驾驶员的可操纵安全范围，并得到了飞机的稳定平衡点；然后在选定的稳定平衡点处基于微分流形理论进行安全边界确定。通过将2种方法相结合，实现了结冰条件下，即使飞行状态超出原安全边界，仍然可以通过操纵指导使其到达新的安全边界。

1 结冰飞机非线性模型

1.1 飞机纵向通道动力学模型

结冰飞机非线性动力学方程可表示成如下微分方程形式：

(1)

式中：状态变量x中包含飞行迎角α、俯仰角θ以及俯仰角速度q；η为结冰因子，反映结冰程度。

建立的飞机三维非线性纵向通道动力学模型如下：

(2)

式中：m为飞机质量；L为飞机升力；V为飞行速度；M为飞机纵向通道力矩；Iy为飞机纵向通道转动惯量；g为重力加速度。

升力和力矩可以表示为

(3)

式中：Tx和Tz分别为推力在机体轴上的分量，T为飞机推力；Fx和Fz分别为气动力在机体轴上的分量；ρ为大气密度；S为飞机参考面积；c为平均空气动力学弦长；Cx、Cz、Cm为力及力矩系数，具体见参考文献[9]。

1.2 结冰飞机多项式系数模型

对飞行数据进行多项式参数拟合能够精确地得到气动数据随飞行迎角α、俯仰角速度q以及升降舵偏角δe变化的多项式形式的解析表达式[16]，如下：

(4)

式中：多项式系数xi(i=1, 2, 3, 4) 和zi(i=1, 2, 3, 4)为结冰因子η的函数。

可通过如下表达式体现结冰程度对飞机气动特性的影响[2]：

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

(5)

式中：KC(A)表示飞行参数的变化；C(A)为结冰前力或力矩系数。结冰因子η是一个时变量，随结冰程度变化，变化范围从0到ηmax。η=0表示飞机不受结冰影响，η=ηmax表示飞机完全结冰。另外，式(5)所用的结冰程度影响模型是由Bragg提出的初步模型，虽然模型较为简单，但在一定程度上能够反映结冰的趋势，根据公开发表的文献[4,17-18]，可知该方法得到的结果对于本文使用的GTM飞机具有较强的通用性，且比较精确可靠。

1.3 增稳控制器设计

为了提升飞机的稳定性，对飞行系统进行反馈增稳控制。选取飞行迎角α、俯仰角θ以及俯仰角速度q作为控制变量，对整个飞行系统进行增稳设计，具体结构如图1所示。图中：θc为驾驶员操纵指令输入；Kα、Kq为相应的反馈控制律参数；θe为指令输入误差；KP为指令放大系数。

图1 飞机纵向通道控制器结构Fig.1 Structure for longitudinal channel controller of aircraft

2 安全性分析方法

本文的主要研究内容是将微分流形理论与分岔分析方法相结合，并共同用于指导结冰条件下的驾驶员操纵。下面对这2种方法的基础理论进行介绍。

2.1 分岔分析方法

分岔分析方法可用于如下的非线性微分系统：

(6)

式中：f(·)为非线性微分系统的函数；x为状态变量；u为控制输入。分岔分析方法可以用于展示全维飞行状态，但对于本文研究的飞机纵向通道而言，主要展示(V，α，θ，q)这4维。在分岔分析中，分岔点的形成主要是由于系统平衡点的稳定性质发生变化所导致。当出现结冰情况，飞机的动力学系统将被破坏，这也会对平衡点的稳定性质造成影响。伴随着结冰程度的不断增加，原本稳定的平衡点也可能变成不稳定平衡点，这样就形成了分岔点。同样的，当驾驶员操纵过猛也会造成动力学系统失稳，特别是对于结冰飞机，驾驶员的操纵会变得更加敏感。因此，本文分别将操纵指令输入δe及结冰因子η作为分岔参数进行研究。

2.2 微分流形理论

稳定性是衡量飞行安全的一项主要指标，只要能将飞行状态保持在稳定域内，飞机就不会有失稳的风险[19-20]。本文提出的基于微分流形理论方法求稳定域是通过以稳定平衡点稳定边界上的不稳定平衡点的稳定流形作为安全边界，该方法能精确地刻画出飞行系统的稳定域。

2.2.1 基本理论

非线性动力学系统可表示为以下形式：

(7)

对于系统式(7)，若存在点x0，使f(x0)=0，则x0称为系统的平衡点。若x0处的Jacobian矩阵的特征值实部全不为零，则x0为双曲平衡点。若特征值实部全为负，则x0为稳定平衡点；反之，若其中存在实部大于零的特征值，则x0是不稳定平衡点。对于不稳定平衡点，实部小于零的部分的特征向量{v1,v2,…,vl}张成稳定的特征空间Es，大于零的部分对应的特征向量{vl+1,vl+2,…,vn}张成不稳定的特征空间Eu，全空间表示为Rn=Es+Eu。

双曲平衡点x0的稳定流形Ws(x0)和不稳定流形Wu(x0)可表示如下：

(8)

式中：φ(t,x)为系统式(7)始于x的解。

2.2.2 非线性系稳定域具体构造方法

本文的稳定域构造主要是基于轨道弧长法[21]，流程如图2所示。具体步骤如下：

步骤1确定初始圆。利用不稳定平衡点xi所对应的稳定特征向量组成特征子空间(即切平面)，然后在切平面上以xi为中心，r为半径作初始圆，并在圆上均匀取N个点{p1,1,p1,2,…,p1,N}。

步骤2求解轨线。以点集{p1，1,p1，2,…,p1，N}为初始点，分别对系统进行反时间积分，当轨线长度达到设定值L0时停止，得到第一代轨线{T1,1,T1,2,…,T1,N}，各条轨线的终点记为{p2,1,p2,2,…,p2,N}。

步骤3轨线疏密程度判断。首先检查各条轨线间的距离，若2条轨线间的距离大于Dmax，则在这2条轨线的初始点间重新插入一个新初始点。若2条轨线间的距离小于Dmin，则删除其中一条轨线的初始点。最后重新调整初始点集，返回步骤2，直到满足要求，进入步骤4。

图2 轨道弧长法构造稳定边界流程Fig.2 Flowchart of stability boundary with track arc length method

步骤4以步骤3中确定的轨线终点为二代初始点，重复步骤2和步骤3，迭代达到Zmax次时停止。

步骤5连接相邻代的轨线。将始于相邻初始点的轨线连接，形成边界面。

步骤6将稳定边界上所有的不稳定平衡点所确定的稳定流形组合在一起形成系统的稳定域。

3 结冰飞机案例分析

本文以NASA的GTM为研究对象，考虑增稳控制结构的影响，利用1.1节给出的结冰飞机纵向通道模型，对飞行系统的动力学特性及非线性稳定域进行分析。飞机的相关动力学参数如下：飞机质量m=90 703 kg，大气密度ρ=0.538 4 kg/m3，重力加速度g=9.806 3 m/s2，飞行高度H=6 000 m，俯仰轴的转动惯量Iy=14 679 000 kg/m2，机翼参考面积S=185.25 m2，平均空气动力学弦长c=4.87 m，飞机具体的气动力及力矩参数见参考文献[9]。

3.1 分岔分析

本文中分岔分析方法主要是用于研究结冰后飞机的动力学特性变化，通过展示飞行状态平衡点的形式揭示。为了便于分析，研究过程中各个控制参数保持不变，而驾驶员操纵指令作为变化量，其中平衡点的计算是通过计算飞机的动力学方程得到。仿真分别针对不同结冰条件给出(V，α，θ，q)的分岔分析图，如图3～图5所示。

图3 无结冰(η=0)分岔图Fig.3 Bifurcation diagram without icing (η=0)

图4 结冰因子η=0.1分岔图Fig.4 Bifurcation diagram for icing factor η=0.1

图5 结冰因子η=0.3分岔图Fig.5 Bifurcation diagram for icing factor η=0.3

图3为飞机不结冰(η=0)时的分岔图，其中红色和蓝色部分分别表示稳定平衡点和不稳定平衡点，指令输入δe的范围为-10°～10°之间变化。可以看出，当飞机不结冰时指令输入δe在-2.2°～10°之间都存在稳定的平衡点，当操纵指令值低于2.2°，飞机就会存在失速的风险。文献[9]研究了一旦飞机失稳，通常是将驾驶杆控制到对应的指令为δe=0°时的稳定平衡状态附近进行边界恢复，这也是一种最直接有效的方法。

图4为飞机轻度结冰(η=0.1)时的分岔图。通过与图3对比可知，当飞机遭受结冰影响，飞机动力学模型结构也发生变化，从而平衡点分布也被改变。相对于不结冰情形，稳定平衡点的范围变化不大，驾驶员指令δe在-1.8°～10°之间都存在稳定的平衡点，因此轻度结冰对飞机的动力学特性的影响较小。

图5为飞机重度结冰(η=0.3)时的分岔图。可以看出，随着结冰程度增加，在同一操纵指令下不仅平衡点的分布发生巨大变化，稳定平衡点的稳定性质也发生较大变化，原本稳定的平衡点也因结冰的影响变得不稳定。同时结冰程度的加剧也使驾驶员的可操纵范围相应减小，驾驶员指令δe只有在0°～6.8°范围内存在稳定平衡点，从而造成驾驶员的操纵变得更加敏感，轻易推拉杆都容易导致飞机失稳。

3.2 微分流形理论确定安全边界

3.2.1 精确性验证

本节中蒙特卡罗方法被用于对微分流形方法进行精确性验证。蒙特卡罗方法作为一种传统方法，已经广泛运用于飞行仿真验证。其原理是：首先在工作点(任意取的某一稳定平衡点)所确定的状态空间内大量取初始点，然后在初始点处进行正向积分，最后将能够收敛于工作点的初始点保留下来构成稳定域。该方法的特点是初始点取的越多，精度越高。但是初始点取得过多，计算所消耗的时间也越长。下面分别利用微分流形方法及蒙特卡罗方法计算稳定域，仿真环境在Intel (R) Core (TM) i7-4790 处理器下进行，其中微分流形方法选择的圆弧向外增长长度(步长)为0.01，向外扩展圈数为10圈。蒙特卡罗方法在满足精度要求的前提下，将其初始点数控制在50×50×50。分别用蒙特卡罗方法和微分流形方法所计算的稳定域如图6所示。

图6 流形方法精确性验证Fig.6 Accuracy verification for manifold method

图6中，蓝色点为蒙特卡罗方法构成的稳定域，彩色曲面为基于微分流形理论确定的安全边界。可以看出，用微分流形方法确定的稳定域能够与蒙特卡罗方法高度吻合，因此证明了基于微分流形理论构造稳定域具有较高的精度。另外在仿真时间上微分流形方法耗时较短仅需2 min，而蒙特卡罗方法虽然同样能够得到精确的稳定域，但耗时却超过了1 h。

3.2.2 结冰飞机安全边界确定及操纵保护

本节主要研究基于微分流形理论确定结冰飞机非线性稳定域，并将其作为飞机的安全边界。对于式(2)的模型，从3.1节的分岔分析图中选驾驶杆操纵值为δe=5°时的稳定平衡点(0.01，-0.35,0)作为工作点在(α，θ，q)3维进行稳定域确定。针对不同的结冰程度，所确定的稳定域如图7所示。

图7中红色、绿色、黑色区域分别代表无结冰(η=0)、轻度结冰(η=0.1)以及重度结冰(η=0.3)时飞机的稳定域。可以看出，当飞机只是轻度结冰时，稳定域变化不大，飞行状态仍然容易维持在稳定域内。但当飞机发生重度结冰，稳定域已经严重收缩，原先在稳定域内的平衡状态可能会超出稳定域，这给飞机带来失稳的风险，下面给出重度结冰时的飞机操纵保护策略。

当飞机遭遇重度结冰影响，可能会导致飞行状态超出原有安全边界，由分岔分析可知，可通过改变驾驶杆，使飞行状态进入新的安全边界，如图8所示。图中黑色区域为重度结冰(η=0.3)，驾驶杆操纵量δe=5°时飞机的稳定域。A(0.17，0.08，0.06) rad为原先稳定域内的一点，但当飞机重度结冰，飞行状态已经超出原稳定域，如图中黑色区域所示。在点A处进行时域仿真，解轨线为图中蓝色曲线1所示，从解轨线可看出，最终飞行状态振荡发散。通过拉杆将操纵杆拉回至δe=0°处，此时状态点A进入了新的稳定域内，如图中红色区域所示。在A点处进行时域仿真，解轨线为图中绿色曲线2，从解轨线可看出，轨线最后将会收敛于新的稳定平衡点，最终飞行状态达到稳定。

图7 不同结冰程度稳定域Fig.7 Stability region for different degrees of icing

图8 不同操纵下稳定域Fig.8 Stability region for different manipulations

4 结 论

本文以GTM为案例飞机，将飞机气动参数进行多项式拟合并结合结冰因子模型，建立了结冰条件下的飞机纵向通道动力学模型。提出了基于微分流形理论确定飞机的安全边界，然后结合分岔分析方法对结冰后飞机的稳定性质变化进行了研究，最后对于重度结冰情形给出了驾驶员操纵策略。主要结论如下：

1) 分岔分析可以得到飞机在不同结冰条件下的平衡点，并给出安全操纵范围，而微分流形理论正是基于稳定平衡点确定安全边界，相比在状态空间内进行大量搜索寻点，简化了工作量的同时能够更加清晰直观地对飞行状态变化进行观测。

2) 不同结冰程度对于飞机的操稳性能影响不同。对于轻度结冰，稳定域向内收缩程度低于10%，此时在同一稳定平衡状态下飞机仍具有较强抗干扰能力，同时对于驾驶员的可操纵范围影响较小，操纵范围相较不结冰情形只减少5%左右。但对于严重结冰，稳定域向内收缩程度已经超过50%，飞机的抗扰动能力也严重降低，甚至原先的稳定平衡状态将变得不稳定，同时驾驶员的可操纵范围也减少了45%左右，这给飞行安全带来巨大风险。

3) 微分流形理论在不同程度结冰条件下都能够精确地确定飞机高维稳定域，将其与分岔分析方法相结合能够实现即使驾驶员的飞行状态在稳定域以外，仍能够通过指导操纵使其到达一个新的稳定域，最终使飞行状态达到稳定。

另外针对结冰因子模型的局限性，将进一步完善飞机结冰后的动力学非线性模型建立，并针对不同冰型、不同部位结冰以及不同飞机结冰进行研究，同时对飞机结冰保护控制方法进行更加深入研究。