地效飞行器的纵向控制系统研究

罗 瑜，樊 赫

(1.陕西工业职业技术学院 电气工程学院，陕西 咸阳 712000; 2.西北机电工程研究所，陕西 咸阳 712000)

0 引言

地效飞行器(ground effect vehicle)，是一种利用地面效应来实现超低空高速飞行的特种飞行器[1]。国内外有关学者对其展开了一定的研究，国内的一些大学及科研机构在积极开发地效飞行器，如哈尔滨工程大学自动化学院测控技术与惯性导航教研室于2010年完成了原理样机的设计与制作[2]，并以地效飞行器为研究平台进行了相关实验，主要包括地效飞行器数学建模和仿真、导航系统软硬件设计、数据融合算法以及飞行控制系统设计[3]。总的来说，目前，我国已经形成了具有自主知识产权的地效飞行器设计技术，先后研制出若干种小型地效飞行器[4]。但是从地效飞行器飞行控制技术发展来看，地效飞行器的相关控制研究较普通固定翼飞机较为欠缺。如何判断地效飞行器的稳定性是地效飞行器设计的关键技术之一。它的稳定性跟常规的固定翼飞机有着不同之处，飞行高度会对地效飞行器的纵向气动系数产生影响，从而影响纵向稳定性[4]。本文就某型地效飞行器的纵向系统控制进行了研究和设计，分析了纵向系统的稳定性，并在低空飞行包线内的多状态点处，采用线性二次型调节器[2](LQR)进行纵向系统的增稳设计，并采用自抗扰控制方法进行了纵向高度的控制律设计。

1 地效飞行器的纵向运动方程

根据所研究的该型地效飞行器的发动机安装位置，对如图1所示的地效飞行器，进行巡航状态下的纵向运动受力分析。

图1 地效飞行器纵向受力指示图

图1中，L为升力，D为阻力,T为推力，为轴力矩。

可得到如下的方程：

(1)

在基准运动为定直平飞运动的条件下，根据小扰动原理，对式(1)无因次方程[5]进行线性化处理，用状态空间方程可描述为：

(2)

A=

(3)

从A阵可以看出，一般情况下，在地效区时，CLz，CDz，Cmz三者都不等于零，因此该阵的最后一列不全为零。可以看出，纵向的其他状态变量会受到高度变量的影响，这是因为高度变量会影响地效飞行器所受的升力，阻力和俯仰力矩[6]，可以看出与普通固定翼飞机有所不同。

2 纵向静动稳定性判定及分析

根据前述分析，地效飞行器纵向系统矩阵A的特征方程为：

|λIn-A|=0 (n=5)

(4)

将特征方程展开为：

A1λ5+B1λ4+C1λ3+D1λ2+E1λ+F1=0

(5)

忽略速度的变化，则特征方程变为：

Aλ4+Bλ3+Cλ2+Dλ+E=0

(6)

其中：

(7)

2.1 静稳定性判定

地效飞行器的静稳定性是指地效飞行器处于平衡状态时，受到外界小扰动而偏离平衡状态[7]，在干扰消失后，不加操纵能产生气动力和气动力矩，使其具有回到原来平衡状态的趋势。地效飞行器静稳定性的判定跟自身的气动导数有关。静态稳定性可用一定约束条件下的气动全导数来表示。根据数学知识可知，偏导数是约束条件下最全面最严格的全导数[8]。因此，静态稳定性就是求解不同约束下的气动偏导数的组合。

图2 高度对升力系数的影响

从图2和图3可知：在地效作用(高度小于7.5米)的区域内，DCL，DCm不为0，地效作用区外(自由空间)，DCL，DCm几乎等于0。地效区内的升力系数随高度的减小而增大；在地效区内升力系数是高度的非线性函数，且高度越小，升力上升速率越大；地效区内，俯仰力矩系数随着高度的增大而幅值减小，变化缓慢，直至离开地效区俯仰力矩不再随高度变化。

图3 高度对俯仰力矩系数的影响

由此可见，地效飞行器在低空飞行时，气动系数随高度变化明显，是高度的非线性函数，且气动力和力矩是气动系数的函数，因此高度会作用于升力、阻力、俯仰力矩，从而影响地效飞行器纵向系统的特征根的分布，进而对系统的稳定性有影响[13]。

广义上系统稳定性的数学表达式为：

(8)

式中，Y为广义恢复力，X为广义位移量。

根据上述定义，在定直平飞的基准运动条件下，对于地效飞行器有：

(9)

分别表示其高度稳定性，迎角稳定性，速度稳定性。

这三项的无因次[14]表达式为：

(10)

(11)

由于地效飞行器进行低空飞行时，忽略速度变化，它的纵向运动方程中的气动系数本质上是迎角α和无因次高度z的非线性函数，则有：

dCL(α,z)=CLαdα+CLzdz

(12)

又因为：

Cm(α,z)=0⟺Cmαdα+Cmzdz=0

(13)

所以：

CLzCmα-CLαCmz>0

(14)

令：

E′=CLzCmα-CLαCmz

(15)

则有：

(16)

同理可以得到：

(17)

因此可以看出：E′=CLzCmα-CLαCmz>0表示地效飞行器具有定速稳定性。同时这说明，在不考虑速度的变化情况下，从单一的偏导数Cmα<0和CLz<0都不能表示地效飞行器的纵向静稳定性，需要保持E′>0。若平衡位置x0(α0,z0)受某种扰动后到达x1(α1,z1)，假设有α1>α0，z1>z0，由气动系数公式可知，尽管Cmα<0和CLz<0都成立，但当Cα1>Cα0，CL1>CL0时，使得在x1(α1,z1)不具备恢复到原平衡位置x0(α0,z0)的趋势，也就是说系统在x0处不具有静稳定性。同时区别于普通固定翼飞机，普通的固定翼飞机可通过Cmα<0即能判断静稳定性，从而仅靠迎角焦点与重心的相对位置就能判断其静稳定性[16]。根据上述分析，相比普通固定翼飞机，地效飞行器的静稳定性判定就较为复杂。

在忽略速度变化静稳定性分析的基础上，假设推力为定值时，对式(10)中的第一式进行分析可知：

(18)

类似的，可以得到：

(19)

(20)

由地效飞行器的特性可知，CLz<0，CDz>0，所以式(19)第一项小于零，因此在E′=CLzCmα-CLαCmz>0成立的条件下，高度稳定性，迎角稳定性，速度稳定性判定可以归结为：

(21)

因此根据对地效飞行器的静稳定性定义分析可知，可将地效飞行器具有高度稳定性，迎角稳定性，速度稳定性的判定条件表达为：当表达式E′>0，F′>0同时成立时，可以判定地效飞行器具有高度稳定性，迎角稳定性，速度稳定性。另外从状态方程来看，时不变线性系统静稳定的充要条件是特征方程的常数项大于零。也就是方程式(5)中的常数项系数F1>0，忽略速度变化时为式(6)中的常数项E>0。所以从状态方程判定静稳定性的条件为E>0，F1>0，其中，E为定速静稳定性，F1为变速静稳定性。

综上分析，可以看出E与E′，F1与F′仅差一个系数，即从定义判断与从状态方程中的常数项为正判断静稳定性是一致的。因为E>0与E′>0等价，F1>0于F′>0等价。可以求得E，F1是随飞行高度变化的，因此得到结论：即静稳定性是随着高度的变化而变化，不同飞行高度下，地效飞行器纵向系统的静稳定性是不同的。

2.2 动稳定性判定

(22)

在忽略速度变化时，式(22)判据为：

(23)

因此可以判定静稳定判据E>0,F1>0是动稳定判据的必要条件。由于直接根据稳定性定义计算系数过于复杂，这里采取根据配平状态，直接求系统的特征根在复平面内的分布情况直接判定系统的动稳定性。当系统的全部特征根都分布在左半s平面时，可以判定系统具有动稳定性。

2.3 典型状态点处系统的动稳定性分析

在典型状态点处对系统进行动稳定性分析。在巡航状态下，纵向线性小扰动方程为：

(24)

式中，x=[△v△α△q△θ△h]T，u=[δeδT]T，y=[△v△α△q△θ△h]T分别为纵向小扰动方程的状态变量，控制量和系统输出。△v为速度(m/s)；△α为迎角(rad)；△q为俯仰角速度(rad/s)；△θ为俯仰角(rad)；△h为高度(m)；△δe为升降舵偏角(rad/s)；δT为油门杆偏角(rad)。

从时域角度出发分析系统的动稳定性。根据前述的内容，考虑到重点研究地效飞行器在低空巡航状态下的稳定性，因此选取飞行高度为3 m,5 m,7 m,8 m,10 m,30 m，速度为78 m/s处分别进行配平，分析地效飞行器纵向运动的特征根如表1所示。

表1 典型状态下纵向线性系统的特征根

从表1中可以得到以下结论：在这6个典型状态点处，飞行高度从h6到h1依次减小，可以看出，当地效飞行器飞离地效区(飞向非地效区，非地效区指的是飞行高度大于7.5 m)的过程中，短周期运动由发散变成稳定，长周期运动由原来的阻尼震荡变得不稳定。在非地效区，气动导数不受高度影响，运动方程中△v,△α,△q,△θ四个状态变量不受飞行高度的影响，不考虑它在高度上的变化，这时也不具备高度稳定性，对应特征根中λ5=0。随着接近水面，地效作用使得气动导数受高度影响作用明显，运动方程△v,△α,△q,△θ,△h五个状态变量受到高度的影响，这符合前面分析的小扰动模型中的A阵最后一列不为零。另外从长短周期根的变化可以看出，地效影响随着高度的减小而增强。

3 纵向LQR增稳控制设计及仿真

经上述分析可知，地效飞行器在巡航状态时不具备动稳定性，为了解决地效飞行器在受到外界干扰时，能通过控制的涉入[21]，使其恢复原定的航行状态，使得地效飞行器具有良好的稳定性[22]，因此需要设计增稳来提高飞机的动稳定性。通常飞机的纵向增稳都是采用迎角和俯仰角反馈来增大系统阻尼，从而增强系统的动稳定性。该方法设计的增稳虽然简单直观，但是鲁棒性不强[23]。线性二次型技术由于具有良好的鲁棒性，在常规布局的固定翼飞机上已有广泛应用，但在地效飞行器上目前应用较少。这里采用常规LQ方法[24]设计巡航时典型状态下的地效飞行器增稳，来改善地效飞行器的静稳定性和动稳定性。

3.1 纵向LQR增稳控制设计

如何设计控制律实现受控系统性能指标最小的最优控制问题即为LQR问题。LQR方法是目前算法中最为成熟应用最为广泛的一种控制方法，设线性时变系统的状态方程和输出方程为：

(25)

式中，状态变量x(t)∈Rn×1，输出变量y(t)∈Rl×1，控制变量u(t)∈Rm×1，时变系统矩阵、增益矩阵和输出矩阵分别为A(t)∈Rn×n、B(t)∈Rn×m和C(t)∈Rl×1。令系统的输出期望向量为yr(t)∈Rl×1。并定义系统的输出误差向量为：

e(t)=yr(t)-y(t)

(26)

设受控系统的性能指标函数为：

(27)

式中，F∈Rl×l为半正定常数加权矩阵，Q(t)∈Rl×l为半正定时变加权矩阵，R(t)∈Rm×m为正定时变加权矩阵，t0和tf分别为系统响应初始时刻和终止时刻。

对于地效飞行器，为了保证系统在受到外界干扰时，使其恢复干扰前的状态，必须进行增稳控制。那么通过适当的反馈保证地效飞行器具有良好的稳定性问题可以看作是状态调节器的设计问题。对于式(2)描述的地效飞行器的纵向线性定常系统，其中系统的状态方程和输出方程为：

(28)

式中，A∈Rn×n为系统矩阵，B∈Rn×m为控制矩阵，C∈Rl×n为输出矩阵，D为对应的l×m零阵。x(t)为系统的状态变量，u(t)为系统的控制变量。

对地效飞行器设计增稳，选用式(27)中的tf→∞的状态调节器模型，终端时刻tf→∞，是为了得到常值反馈增益矩阵，这时F=0，式(27)中第一项终端性能指标失去意义，又因为系统是时不变的，性能指标的权矩阵Q(t)和R(t)就为常矩阵Q和R，同时系统完全可控，那么选定二次型能指标函数为：

(29)

为了使地效飞行器的纵向系统满足性能指标J最小，首先构造Hamilton函数：

(30)

其次通过求微分的方法求出最优控制信号u(t)：

u(t)=-R-1BTP(t)x(t)

(31)

可组成对应矩阵Riccati微分方程：

(32)

因为上述Riccati微分方程求解比较困难，这里假设P(t)是一个常数矩阵，所以P(t)的一阶微分趋于零，所以上式变成：

0=PA+ATP+Q-PBR-1BTP

(33)

那么，在已知A、B、Q、R，就较容易求解上式Riccati微分方程的对称正定解P。再依据u(t)=-R-1BTPx(t)=Kx(t)，求得u(t)。其中，K为最优常值的反馈增益阵。这里采用Matlab中lqr函数来求解Riccati微分方程。求解K的过程中，Q的选取是根据对自然模态特性的分析，依据各状态变量在各模态中所起的重要程度的不同来进行加权，R的选取是根据操纵面的限制条件来定。

针对式(2)所示的地效飞行器纵向线性系统，搭建如图4的模型，进行多状态反馈，采用最优二次型的设计方法来求取反馈矩阵K。

图4 地效飞行器纵向线性模型

3.2 纵向增稳控制仿真结果分析

为了验证LQR状态调节器的增稳效果，这里选取地效飞行器飞行高度为3 m，飞行速度为78 m/s的巡航状态为例，采用了LQR状态调节器方法设计地效飞行器的纵向增稳控制律。可以求得该状态下的线性方程系数矩阵为：

易知，该状态下是动不稳定的，所以用状态调节器增稳。根据系统的速度和俯仰角速率在模态中起主要作用，以及保证产生合适的舵偏角。这里分别选取加权阵Q和R为：Q=diag([20,5,50,10,0.5])，R=100。

可求得反馈矩阵K为：K=[0.41242.7765-1.4301-4.79740.0458]。增稳后闭环系统特征值分别为λ1=-5.3196，λ2=-0.1025，λ3,4=-2.1954±2.4214i，λ5=-0.8185。则闭环系统在受迎角为5°的扰动时，纵向各个状态的响应曲线如图5所示。

图5 反馈后受迎角扰动纵向各参数的响应曲线

从图5和闭环系统的特征根可以看出，采用LQR增稳后，地效飞行器的纵向稳定性增强，对于外界干扰能够很快恢复，控制舵偏角也符合实际要求，不超过±25°，增稳后的地效飞行器纵向静、动稳定性得到了明显的改善。经过调参，可以得到在不同高度下的最优反馈阵K如表2示。

表2 不同高度下的最优反馈阵

可以看到，尽管高度从h1到h6不断增加，地效飞行器飞行在h1到h3之间时，处于地效区，最优反馈阵可选定为一个值；h4及其以上高度属于非地效区，最优反馈阵K可选定为另外一个值，因此可以考虑在地效区和非地效区分别选定一个最优反馈阵K，来满足不需要在不同高度下不停地切换最优反馈控制律。在飞行高度为3 m，飞行速度78 m/s所设计LQR增稳的Q，R阵用于高度为5 m，7 m，速度均为78 m/s的状态下施加迎角为5度的扰动时，地效飞行器纵向闭环系统各状态的响应曲线如图6所示。

图6 反馈后迎角扰动下纵向参数的响应曲线

仿真结果表明：飞行高度为3 m，速度为78 m/s处所设计LQR增稳控制具有较强的鲁棒性，对于高度为7 m时的控制性能稍弱，是因为该状态点位于临界地效区。同理，可以对高度为30 m，速度为78 m/s时所设计的Q，R阵进行鲁棒性验证。

4 高度控制设计及仿真分析

对地效飞行器的高度控制是在俯仰角稳定回路的基础上加上外回路实现的。这里外回路采用一阶自抗扰控制(ADRC)[25]，将实际高度信息与指令高度信号作差，将高度差通过高度控制输入到俯仰角控制系统，改变航迹倾斜角u来控制地效飞行器的飞行高度，直至稳定到指令高度，实现高度的稳定与控制。

4.1 高度控制设计

地效飞行器的纵向外环高度运动可以用以高度Δh为输出，俯仰角Δθ为输入的一阶微分方程来描述，将迎角α看做扰动控制量，引入一阶自抗扰控制器[26]，实现高度控制。由于：

ΔH=-ΔVΔα+ΔVΔθ

(34)

对此一阶系统引入二阶线性扩张状态观测器观测纵向外环可能存在的扰动总和，引入线的一阶ADRC由3个部分构成[27]。

其中：扩张状态观测器：

(35)

线性误差反馈控制律：

(36)

控制量：

(37)

4.2 高度控制仿真结果分析

在内环俯仰姿态控制的基础上，采用S-Function函数建立高度环一阶扩张状态观测器，则纵向高度控制仿真模型如图7所示。

图7 地效飞行器高度通道自抗扰控制仿真模型

这里选定飞行速度为78 m/s，飞行高度为3 m的状态进行配平得到该状态点下系统的线性模型，给定指令爬升高度5 m，分别采用PID和自抗扰控制进行高度控制。采用PID控制时，由于高度通道是在俯仰通道的基础上进行设计的，所以必须把高度输入值转化成相应的俯仰角度。当得到俯仰角需要改变的角度后，进行俯仰通道控制，则飞行高度也会随之而发生改变。地效飞行器爬升下滑阶段的被控制量为爬升高度Δh，则控制量Δθ则为PID控制的控制量可示为：

(38)

其中：Δhc为俯仰角的指令信号，kp为比例系数，ki为积分系数。

采用自抗扰对高度进行控制是在俯仰通道的基础上对高度进行一节ADRC控制，对高度环的一阶ADRC的三个模块进行调参，进行控制器效果仿真验证。其中设置一阶高度环ADRC参数：h=0.002，boh=65，β01h=9 000，KhP=10.5，KhD=1.73。综合PID和自抗扰仿真，可得高度的阶跃响应如图8所示。

图8 PID控制和自抗扰控制的爬升高度响应

对比图8两条响应曲线，可以发现，相比于PID控制，自抗扰控制下的高度响应迅速准确，且超调较小，在15 s以内能很快地进入稳态。且于阶跃响应，具有良好的控制效果。

5 结束语

本文针对地效飞行器纵向系统不稳定的问题，分析了地效飞行器纵向的静稳定性和动稳定性判定条件，并在典型状态下状态点进行了动稳定性的分析。依据LQR控制进行了LQR纵向增稳系统的设计，针对选取的状态点纵向运动的动不稳定，设计了LQR增稳控制，并在飞行包线内选取了其他的状态点进行仿真验证。针对纵向高度稳定控制的问题，采用自抗扰进行了高度控制律设计并进行仿真验证。仿真结果表明，所设计的LQR增稳和自抗扰高度控制具有良好的鲁棒性、准确性和快速性。