随着航空技术的快速发展，固定翼无人机在航空领域扮演着越来越重要的角色，固定翼无人机的吸引力及其在实际中的应用在过去的十年间显著地增加。由于固定翼无人机强机动性的特点，使得固定翼无人机在军事、民用和商业领域的应用前景越来越广阔。同时，由于任务的多样性和复杂性，对固定翼无人机的控制性能的要求越来越高，控制性能的好坏极大影响着无人机飞行的安全性和稳定性。

近年来，针对具有强耦合、强非线性以及快速机动性特点的固定翼无人机，国内外的许多研究人员和学者采用多种控制方法开展了对控制问题的相关研究^[1]。2016年，上海交通大学的费爱玲等^[2]针对固定翼无人机的四元数模型，采用基于扩张观测器的反步法设计了无人机速度控制器和姿态控制器，实现了对无人机速度和姿态的高精度控制。2015年，哈尔滨工程大学的刘一鸣^[3]设计了基于反演法的固定翼无人机自适应控制律，并将该方法与智能自适应PID(Proportion Integral Derivative)方法进行了对比，仿真结果表明，文中所提出的自适应控制方法准确性更高。2016年，TOBB大学的Akyurek等^[4]采用H_∞方法设计了小型固定翼无人机的姿态控制器，提高了无人机控制系统的鲁棒性。2017年，挪威科技大学的Rotondo等^[5]采用线性参数变化模型参考方法针对复杂任务环境需求，设计了固定翼无人机控制器。2016年，美国肯塔基大学的Mullen等^[6]考虑无人机在姿态跟踪过程中的参数不确定以及风干扰的影响，采用FDI(Filtered Dynamic Inversion)方法针对小型固定翼无人机设计了姿态控制器，并在飞行试验中将该方法与传统的PI(Proportion Integral)方法进行了对比，试验数据表明，采用FDI方法设计的姿态控制器具有更为优越的控制性能。2016年，南洋理工大学的Hervas等^[7]针对固定翼无人机在舰船甲板的着陆问题，考虑随机风干扰的影响，采用基于扩张卡尔曼滤波的滑模控制方法的控制策略设计了控制器。2017年，墨西哥机器人和智能系统中心的Castañeda等^[8]采用扩张状态观测器对干扰以及未知状态进行了估计，并基于自适应超螺旋算法设计了姿态控制器和速度控制器，保证了系统的有限时间收敛并有效地减弱了滑模算法中的抖振效应。

当前，自适应滑模控制方法^[9-14]被广泛应用到航空航天飞行器的控制器设计中^[15-20]。本文考虑固定翼无人机实际飞行过程中的未知干扰、执行器饱和以及强耦合的影响，针对无人机姿态控制和速度控制的不同需求，设计一种新型的自适应多变量螺旋滑模姿态控制器和自适应多变量快速超螺旋滑模速度控制器，进而实现了对固定翼无人机姿态和速度的有限时间快速精确跟踪。主要创新点如下: ①文中设计的自适应多变量螺旋和自适应快速超螺旋算法无需预知系统所受干扰及其导数的上界便可实现对系统的精确稳定控制; ②针对典型的二阶姿态模型，采用自适应多变量螺旋滑模设计姿态控制器，自适应律的引入，有效地抑制了抖振; ③文中设计的新型自适应快速超螺旋速度控制器具有收敛速度快，收敛精度高的优点，适用于固定翼无人机对快速机动性的需求。

1 固定翼无人机模型

固定翼无人机姿态运动涉及惯性坐标系和机体坐标系等。相关坐标系的关系如图 1所示。

假设飞机为刚体模型、地面坐标系为惯性坐标系，固定翼无人机模型为^[8]

$ {\mathit{\boldsymbol{I\dot \omega }} = - \mathit{\boldsymbol{\omega }} \times \mathit{\boldsymbol{I\omega }} + \mathit{\boldsymbol{M}}} $	（1）

$ {\mathit{\boldsymbol{ \boldsymbol{\dot \varTheta} }} = {\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \cdot \mathit{\boldsymbol{\omega }}} $	（2）

$ {\mathit{\boldsymbol{\dot v}} = (\mathit{\boldsymbol{F}} + \mathit{\boldsymbol{T}})/m + {\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{I}}}\mathit{\boldsymbol{g}} - \mathit{\boldsymbol{\omega }} \times \mathit{\boldsymbol{v}}} $	（3）

$ {{{\mathit{\boldsymbol{\dot p}}}^n} = {\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{I}}}\mathit{\boldsymbol{v}}} $	（4）

式中：m为无人机质量；ω=[p q r]^T为无人机的角速度，p为滚转角速度，q为俯仰角速度，r为偏航角速度；Θ=[ϕ θ ψ]^T为姿态角，ϕ为滚转角，θ为俯仰角，ψ为偏航角；v=[u v w]^T为飞行速度，u为无人机机体坐标系x轴方向的速度，v为无人机机体坐标系y轴方向的速度，w为无人机机体坐标系z轴方向的速度；pⁿ=[x_g y_g h]^T为无人机的位置信息；M=[M₁ M₂ M₃]^T为调整无人机姿态的控制力矩, 包含滚转M₁、俯仰M₂和偏航M₃3个力矩；F∈R^3×1和T=[T_x 0 0]^T分别为无人机所受的气动力和发动机提供的推力, T_x为沿机体坐标系x方向推力；R_Θ∈R^3×3，为转换矩阵，I∈R^3×3, 为惯性矩阵，R_I∈R^3×3表示为

$ {\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} = \left[ {\begin{array}{*{20}{c}} 0&{{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi }&{ - {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi }\\ 0&{\frac{{{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi }}{{{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta }}}&{\frac{{{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi }}{{{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta }}}\\ 1&{{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi {\rm{tan}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta }&{{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi {\kern 1pt} {\rm{tan}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta } \end{array}} \right],\mathit{\boldsymbol{I}} = \left[ {\begin{array}{*{20}{c}} {{I_{xx}}}&0&{{I_{xz}}}\\ 0&{{I_{yy}}}&0\\ {{I_{zx}}}&0&{{I_{zz}}} \end{array}} \right] $

$ {\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{I}}} = \left[ {\begin{array}{*{20}{c}} {{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \psi }&{{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \psi {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi - {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \psi {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi }&{{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \psi {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi + {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta }\\ {{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \psi }&{{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \psi {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi + {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \varphi {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi }&{{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \psi {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi - {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \psi {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi }\\ { - {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta }&{{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi }&{{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi } \end{array}} \right] $

2 基于自适应滑模算法的控制器设计

由于固定翼无人机的推进系统仅能提供沿机体坐标系x轴方向的推力T_x。所以速度控制器基于空速方程式(5)设计。

$ \begin{array}{*{20}{l}} {m\dot V = {T_x}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \alpha {\rm{cos}}{\kern 1pt} {\kern 1pt} \beta - D - mg[ - {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \alpha {\rm{cos}}{\kern 1pt} {\kern 1pt} \beta {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{sin}}{\kern 1pt} {\kern 1pt} \beta {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta + {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \alpha {\rm{cos}}{\kern 1pt} {\kern 1pt} \beta {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \phi {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta ]} \end{array} $

（5）

式中：V为无人机的空速；α=arctan(w/u)和β=arcsin(u/v)分别为迎角和侧滑角; D为阻力。

在固定翼无人机的控制中，设计飞行控制器的目的是精确跟踪航路点指令p_dⁿ，而在实际的控制过程中，通常将p_dⁿ解析为空速参考指令V_d和Θ_d，进而，控制设计的目的可转化为精确跟踪参考指令V_d和Θ_d。本节针对固定翼无人机控制中的强干扰以及输入饱和问题，基于新型自适应多变量螺旋算法设计二阶姿态控制器；基于与速度子系统相关的空速动态方程，采用自适应快速超螺旋算法设计速度控制器，实现对参考指令的有限时间精确跟踪。控制框图如图 2所示。

2.1 姿态控制器设计

本节针对姿态子系统式(1)和式(2)，考虑执行器饱和以及外界干扰Δd_ω=[Δd_ω1 Δd_ω2 Δd_ω3]^T的影响，设计姿态控制器，定义姿态角跟踪误差为

$ {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} = \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} - {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_{\rm{d}}} $	（6）

式中:Θ_d为姿态参考指令，基于式(2)，对式(6)求导得：

$ {\mathit{\boldsymbol{\dot e}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} = \mathit{\boldsymbol{ \boldsymbol{\dot \varTheta} }} - {\mathit{\boldsymbol{ \boldsymbol{\dot \varTheta} }}_{\rm{d}}} = {\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \cdot \mathit{\boldsymbol{\omega }} - {\mathit{\boldsymbol{ \boldsymbol{\dot \varTheta} }}_{\rm{d}}} $	（7）

设$\dot{\boldsymbol{e}}_{\boldsymbol{\varTheta}}=\boldsymbol{z}_{\boldsymbol{\varTheta}} $，则

$ {\mathit{\boldsymbol{\dot z}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} = {\mathit{\boldsymbol{\dot R}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \cdot \mathit{\boldsymbol{\omega }} + {\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}{\mathit{\boldsymbol{I}}^{ - 1}}( - \mathit{\boldsymbol{\omega }} \times \mathit{\boldsymbol{I\omega }} + \mathit{\boldsymbol{M}} + \Delta {\mathit{\boldsymbol{d}}_\omega }) - {\mathit{\boldsymbol{ \boldsymbol{\ddot \varTheta} }}_{\rm{d}}} $

（8）

考虑固定翼无人机输入饱和的影响，式(8)可转化为式(9)的形式

$ {\mathit{\boldsymbol{\dot z}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} = {\mathit{\boldsymbol{\dot R}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \cdot \mathit{\boldsymbol{\omega }} + {\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}{\mathit{\boldsymbol{I}}^{ - 1}}( - \mathit{\boldsymbol{\omega }} \times \mathit{\boldsymbol{I\omega }} + {\bf{sa}}{{\bf{t}}_M} + \Delta {\mathit{\boldsymbol{d}}_\omega }) - {\mathit{\boldsymbol{ \boldsymbol{\ddot \varTheta} }}_{\rm{d}}} $

（9）

式中:sat_M∈R^3×1为非线性饱和特性函数。令sat_M=f+M=[sat_M1 sat_M2 sat_M3]^T，其中，f=[f₁ f₂ f₃]^T，表达式为

$ \begin{array}{*{20}{l}} {{f_i}(i = 1,2,3) = }\\ {\;\:\left\{ {\begin{array}{*{20}{c}} { {\rm{sign}} ({M_i}){M_{{\rm{max}}}} - {M_i}}\\ 0 \end{array}\;\:\begin{array}{*{20}{l}} {|{M_i}| \ge {M_{{\rm{max}}}}}\\ {|{M_i}| < {M_{{\rm{max}}}}} \end{array}} \right.} \end{array} $	（10）

其中：M_max为执行机构所允许的最大力矩值，通过式(10)可得：当控制器产生的力矩|M_i|大于M_max时，sat_Mi=M_maxsign(M_i), i=1, 2, 3；当控制器产生的力矩|M_i|小于M_max时，sat_Mi=M_i，通过函数sat_M∈R^3×1将输出的控制力矩限制在了执行器允许的范围内。

通过上述变换，得到姿态子系统的面向控制模型式(11)

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} = {\mathit{\boldsymbol{z}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}}\\ {{{\mathit{\boldsymbol{\dot z}}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} = {\mathit{\boldsymbol{G}}_0} + {\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}{\mathit{\boldsymbol{I}}^{ - 1}}\mathit{\boldsymbol{M}} + \mathit{\boldsymbol{\rho }}} \end{array}} \right. $

（11）

式中:G₀和ρ分别为

$ {{\mathit{\boldsymbol{G}}_0} = {{\mathit{\boldsymbol{\dot R}}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \cdot \mathit{\boldsymbol{\omega }} + {\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}{\mathit{\boldsymbol{I}}^{ - 1}}( - \mathit{\boldsymbol{\omega }} \times \mathit{\boldsymbol{I\omega }} + \mathit{\boldsymbol{f}}) - {{\mathit{\boldsymbol{ \boldsymbol{\ddot \varTheta} }}}_{\rm{d}}}} $

$ {{\eta _1} = {\rm{min}}\left( {{\eta _0},\frac{{{\alpha _1}}}{{\sqrt {2{\beta _1}} }}} \right)} $

||ρ||≤D，D>0且未知。

基于式(11)，采用自适应多变量螺旋算法设计姿态控制器如式(12)所示

$ \mathit{\boldsymbol{M}} = \mathit{\boldsymbol{R}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}^{ - 1}\mathit{\boldsymbol{I}}\left[ { - {\mathit{\boldsymbol{G}}_0} - {k_1}\left( {\frac{{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}}}{{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}} \right\|}} + 0.5 \times \frac{{{\mathit{\boldsymbol{z}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}}}{{\left\| {{\mathit{\boldsymbol{z}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}} \right\|}}} \right)} \right] $

（12）

式中：k₁>0为自适应增益，满足式(13)中所示的自适应律：

$ {\dot k_1} = \left\{ {\begin{array}{*{20}{c}} {{\alpha _1}\sqrt {\frac{{{\beta _1}}}{2}} }&{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}} \right\| \ge {\varepsilon _1}}\\ 0&{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}} \right\| < {\varepsilon _1}} \end{array}} \right. $	（13）

其中:α₁，β₁，ε₁为大于零的常数。

将姿态控制律式(12)代入式(11)中，得式(14)

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} = {\mathit{\boldsymbol{z}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}}\\ {{{\mathit{\boldsymbol{\dot z}}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} = - {k_1}\left( {\frac{{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}}}{{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}} \right\|}} + 0.5 \times \frac{{{\mathit{\boldsymbol{z}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}}}{{\left\| {{\mathit{\boldsymbol{z}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}} \right\|}}} \right) + \mathit{\boldsymbol{\rho }}} \end{array}} \right. $

（14）

定理1  对于姿态子系统式(1)和式(2)，在自适应多变量螺旋滑模控制器式(12)和自适应律式(13)的作用下，选取适当的α₁、β₁和ε₁，姿态子系统可实现有限时间稳定，姿态角跟踪误差可以在有限时间内收敛到原点附近的邻域内。

证明  为证明系统在控制器式(12)~式(13)的作用下可以实现有限时间稳定，选择如下的Lyapunov候选函数

$ {V_1} = \underbrace {\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}} \right\|{\mathit{\boldsymbol{z}}^{\rm{T}}}\mathit{\boldsymbol{Az}} + \frac{1}{4}{{(\mathit{\boldsymbol{z}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}^{\rm{T}}{\mathit{\boldsymbol{z}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}})}^2}}_{{V_{10}}} + \underbrace {\frac{1}{{{\beta _1}}}{{({k_1} - k_1^*)}^4}}_{{V_{{\rm{1ad}}}}} $

（15）

式中:V₁₀和V_1ad为组成V₁的2个部分，方便后续进行稳定性分析。$z=\left[\frac{\boldsymbol{e}_{\boldsymbol{\varTheta}}^{\mathrm{T}}}{\left\|\boldsymbol{e}_{\boldsymbol{\varTheta}}\right\|^{1 / 2}} \quad \boldsymbol{z}_{\boldsymbol{\varTheta}}^{\mathrm{T}}\right]^{\mathrm{T}} $。其中$\gamma<2 k_{1}^{\frac{3}{2}}, k_{1}^{*}>0, \beta_{1}>0, \boldsymbol{A}=\left[\begin{array}{cc} k_{1}^{2} & \frac{\gamma}{2} \\ \frac{\gamma}{2} & k_{1} \end{array}\right] $为正定对称矩阵。对于系统式(14)的稳定性证明可以分为2步：①对V₁₀进行分析; ②对V₁进行分析。

对V₁₀进行稳定性分析。由不等式λ_min{A}||z||²≤z^TAz≤λ_max{A}||z||²(λ_min{A}为矩阵的A最小特征值，λ_max{A}为矩阵A的最大特征值)可得

$ \begin{array}{*{20}{l}} {{V_{10}} \le \frac{3}{2}{\lambda _{{\rm{max}}}}(\mathit{\boldsymbol{A}}){{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}} \right\|}^2} + \frac{1}{2}({\lambda _{{\rm{max}}}}(\mathit{\boldsymbol{A}}) + \frac{1}{2}){{\left\| {{\mathit{\boldsymbol{z}}_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}} \right\|}^4}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \le {\mathit{\boldsymbol{\chi }}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\chi }} \le {\lambda _{{\rm{max}}}}({\mathit{\boldsymbol{P}}_1})({{\left\| {{\mathit{\boldsymbol{e}}_\theta }} \right\|}^2} + {{\left\| {{\mathit{\boldsymbol{z}}_\theta }} \right\|}^4})} \end{array} $

（16）

式中: $ \boldsymbol{P}_{1}=\left[\begin{array}{cc} \frac{3}{2} \lambda_{\max }(\boldsymbol{A}) & 0 \\ 0 & \frac{1}{2} \lambda_{\max }(\boldsymbol{A})+\frac{1}{4} \end{array}\right]>0$; $ \boldsymbol{\chi}=\left[\begin{array}{ll} \mid e_{\boldsymbol{\varTheta}} \| & \left\|\boldsymbol{z}_{\boldsymbol{\varTheta}}\right\|^{2} \end{array}\right]^{\mathrm{T}}$。

基于文献[9]和[17]，可以得到:

$ {\dot V_{10}} \le - \frac{\vartheta }{{{2^{2/3}}\lambda _{{\rm{max}}}^{3/4}\{ {\mathit{\boldsymbol{P}}_1}\} }}V_{10}^{3/4} = - {\eta _0}V_{10}^{3/4} $	（17）

式中：$\eta_{0}=\frac{\vartheta}{2^{2 / 3} \lambda_{\max }^{3 / 4}\left\{\boldsymbol{P}_{1}\right\}}, \lambda_{\max }\left\{\boldsymbol{P}_{1}\right\} $为矩阵P₁的最大特征值，$ \vartheta=\min \left\{\lambda_{\min }(\boldsymbol{P}), \gamma\left(\frac{1}{2} k_{1}-D\right)\right\}$。

由式(17)可得

$ {\dot V_1} \le - {\eta _1}V_1^{3/4} + \tilde \omega $	（18）

式中：

$ {\tilde \omega = \left( { - \frac{1}{{4{\beta _1}}}{{\dot k}_1} + \frac{{{\alpha _1}}}{{4\sqrt {2{\beta _1}} }}} \right)|{k_1} - k_1^*{|^3}} $

当$\dot{k}_{1}=\alpha_{1} \sqrt{\beta_{1} / 2} $时，$ \tilde{\omega}$为零，则

$ {\dot V_1} \le - {\eta _1}V_1^{3/4} $	（19）

当自适应增益k₁按照式(13)中的自适应律增长到满足k₁>2D时，e_Θ可在有限时间内收敛到限定区间ε₁内。

2.2 速度控制器设计

基于式(5)，考虑未知外界干扰Δd_v的影响，得到空速误差方程:

$ {\dot e_v} = \dot V - {\dot V_d} = \gamma u + {\varPhi _1} + {\varPhi _2} $	（20）

式中: $ \gamma=\frac{\cos \alpha \cos \beta}{m} ; u=T_{x} ; \varPhi_{1}=-\frac{D}{m}-$g[-cosαcosβsinθ+sinβsinϕcosθ+sinαcosβcosϕcosθ]; Φ₂=Δd_v为速度环所受的外界干扰，且满足$ \left\|\dot{\varPhi}_{2}\right\| \leqslant \bar{\delta}$。

根据式(20)，设计速度控制器为

$ \left\{ {\begin{array}{*{20}{l}} {u = - \frac{1}{\gamma }{\varPhi _1} - \frac{1}{\gamma }({k_{v1}}{\xi _1}({\mathit{\boldsymbol{e}}_v}) + {z_1})}\\ {{{\dot z}_1} = - {k_{v2}}{\xi _2}({\mathit{\boldsymbol{e}}_v})} \end{array}} \right. $	（21）

ξ₁(e_v)和ξ₂(e_v)表示如下：

$ \left\{ \begin{array}{l} {\xi _1}({\mathit{\boldsymbol{e}}_v}) = \frac{{{\mathit{\boldsymbol{e}}_v}}}{{{{\left\| {{\mathit{\boldsymbol{e}}_v}} \right\|}^{1/2}}}} + {k_{v3}}{\mathit{\boldsymbol{e}}_v}\;\:{k_{v3}} > 0\\ {\xi _2}({\mathit{\boldsymbol{e}}_v}) = \frac{{{\rm{d}}{\xi _1}({\mathit{\boldsymbol{e}}_v})}}{{{\rm{d}}{\mathit{\boldsymbol{e}}_v}}}{\xi _1}({\mathit{\boldsymbol{e}}_v}) = \frac{1}{2} \cdot \frac{{{\mathit{\boldsymbol{e}}_v}}}{{\left\| {{\mathit{\boldsymbol{e}}_v}} \right\|}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{3}{2}{k_{v3}}\frac{{{\mathit{\boldsymbol{e}}_v}}}{{{{\left\| {{\mathit{\boldsymbol{e}}_v}} \right\|}^{1/2}}}} + k_{v3}^2{\mathit{\boldsymbol{e}}_v} \end{array} \right. $

（22）

k_v3>0为固定增益；k_v1, k_v2为大于0的自适应增益，

$ \left\{ {\begin{array}{*{20}{l}} {{{\dot k}_{v1}} = \left\{ {\begin{array}{*{20}{c}} {{{\bar k}_{v1}} \cdot \left\| {{\mathit{\boldsymbol{e}}_v}} \right\| \cdot {\rm{sign}} (\left\| {{\mathit{\boldsymbol{e}}_v}} \right\| - {\varepsilon _2})}&{{k_{v1}} > \mu }\\ \mu &{{k_{v1}} \le \mu } \end{array}} \right.}\\ {{k_{v2}} = {\eta _v}{k_{v1}}} \end{array}} \right. $	（23）

其中:k_v1、ε₂、μ和η_v均为大于零的常数。

将式(21)代入式(20)可得：

$ \left\{ {\begin{array}{*{20}{l}} {{{\dot e}_v} = - {k_{v1}}{\xi _1}({\mathit{\boldsymbol{e}}_v}) + {z_v}}\\ {{{\dot z}_v} = - {k_{v2}}{\xi _2}({\mathit{\boldsymbol{e}}_v}) + \Delta {{\dot d}_v}} \end{array}} \right. $	（24）

定理2  对于速度子系统，在速度控制器式(21)和自适应增益式(23)的作用下，选取合适的参数ε₂、k_v3、μ、η_v和k_v1，速度子系统的跟踪误差e_v可以在有限时间内收敛到零。

证明

采用Lyapunov方法对系统式(24)进行稳定性分析，选择Lyapunov候选函数为

$ {V_2} = \underbrace {{\zeta ^{\rm{T}}}\mathit{\boldsymbol{P \boldsymbol{\varsigma} }}}_{{V_{20}}} + \underbrace {\frac{1}{{2{\gamma _1}}}{{({k_{v1}} - k_{v1}^*)}^2} + \frac{1}{{2{\gamma _2}}}{{({k_{v2}} - k_{v2}^*)}^2}}_{{V_{{\rm{2ad}}}}} $	（25）

式中:ς=[ξ₁ z_v]^T，k_v1^*和k_v2^*为大于零的常数，$ \boldsymbol{P}=\left[\begin{array}{cc} \lambda^{2}+4 \varepsilon & -\lambda \\ -\lambda & 1 \end{array}\right], \lambda>0, \varepsilon>0$为对称正定矩阵; V₂₀和V_2ad为V₂的2个部分。

对式(25)求导得：

$ \begin{array}{*{20}{l}} {{{\dot V}_2} = {{\mathit{\boldsymbol{ \boldsymbol{\dot \varsigma} }}}^{\rm{T}}}\mathit{\boldsymbol{P \boldsymbol{\varsigma} }} + {{\mathit{\boldsymbol{ \boldsymbol{\dot \varsigma} }}}^{\rm{T}}}\mathit{\boldsymbol{P \boldsymbol{\varsigma} }} + \frac{1}{{{\gamma _1}}}({k_{v1}} - k_{v1}^*){{\dot k}_{v1}} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{{{\gamma _2}}}({k_{v2}} - k_{v2}^*){{\dot k}_{v2}}} \end{array} $

（26）

设$ k_{2}=k_{v 2}-\frac{\Delta \dot{d}_{v}}{\xi_{2}}$，则式(26)可简化为

$ {\dot V_{20}} = - 2\xi _1^\prime {\mathit{\boldsymbol{ \boldsymbol{\varsigma} }}^{\rm{T}}}\mathit{\boldsymbol{Q\zeta }} $	（27）

式中: $ \boldsymbol{Q}=\left[\begin{array}{ll} q_{1} & q_{2} \\ q_{2} & q_{3} \end{array}\right]$为对称矩阵，

$ {{q_1} = {k_{v1}}{\lambda ^2} + 4{k_{v1}}\varepsilon - \lambda {k_2} = {a_1}{k_{v1}} - \lambda {k_2}} $

$ \begin{array}{l} {q_2} = \frac{1}{2}( - {\lambda ^2} - 4\varepsilon - \lambda {k_{v1}} + {k_2}) = \\ \;\:{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2} \cdot ( - {a_1} - \lambda {k_{v1}} + {k_2}) \end{array} $

$ {{q_3} = \lambda } $

其中:a₁=λ²+4ε₀。

若$\dot{V}_{20} $为负定，则矩阵Q为正定，即矩阵的行列式需要大于零。为保证det(Q)>0，求根判别式λk₂(k_v1-λ)要大于零。由于$ \left|\Delta \dot{d}_{v}\right| \leqslant \bar{\delta}$，$ \left\|\frac{1}{\xi_{2}}\right\| \leqslant 2$，则可得到：

$ {k_2} > 0 \Rightarrow {k_{v2}} > 2\bar \delta \;\:{k_{v1}} > \lambda $	（28）

k₂的取值范围为

$ {{k}_{2}}\in [{{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{k}}_{2}},{{\bar{k}}_{2}}]=[{{k}_{v2}}-2\bar{\delta },{{k}_{v2}}+2\bar{\delta }] $	（29）

设

$ {{k}_{v1}}=\lambda +\tau \ \;\tau >0 $	（30）

方程det[Q]=0的解为

$ \left\{ \begin{array}{*{35}{l}} p_{1}^{+}=\lambda {{k}_{v1}}+{{k}_{2}}+2\sqrt{\lambda {{k}_{2}}({{k}_{v1}}-\lambda )} \\ p_{1}^{-}=\lambda {{k}_{v1}}+{{k}_{2}}-2\sqrt{\lambda {{k}_{2}}({{k}_{v1}}-\lambda )} \\ \end{array} \right. $	（31）

当p₁∈(p_1max^-, p_1min⁺)时，det[Q]>0。为保证根的存在性，需要满足p_1max^- < p_1min⁺。设k₂=κ²，得到

$ {{\bar{\delta }}^{2}}<\lambda \tau {{\kappa }^{2}}\Rightarrow {{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{k}}_{2}}>\frac{{{{\bar{\delta }}}^{2}}}{\lambda \tau }\Rightarrow {{k}_{2}}>\frac{{{{\bar{\delta }}}^{2}}}{\lambda \tau }+2\bar{\delta } $	（32）

设矩阵P的最小特征值和最大特征值分别为λ_min{P}和λ_max{P}，则式(27)可以转换为式(33)的形式

$ {\dot V_{20}} = - 2\xi _1^\prime {\left\| \mathit{\boldsymbol{ \boldsymbol{\varsigma} }} \right\|^2}{\lambda _{{\rm{min}}}}\{ \mathit{\boldsymbol{Q}}\} \le - {\gamma _1}V_{20}^{1/2} - {\gamma _2}{V_20} $	（33）

式中: $ \gamma_{1}=\frac{\lambda_{\min }\{Q\} \lambda \lambda_{\min }^{1 / 2}\{\boldsymbol{P}\}}{\lambda_{\max }\{\boldsymbol{P}\}}, \gamma_{2}=2 k_{v 3} \frac{\lambda_{\min }\{\boldsymbol{Q}\} \lambda}{\lambda_{\max }\{\boldsymbol{P}\}}$。进一步得到

$ {\dot V_{20}} \le - {\gamma _1}V_{20}^{1/2} $	（34）

式(26)可转化为

$ \begin{array}{l} \begin{array}{*{20}{l}} {{{\dot V}_2} \le - {\beta _\eta }V_2^{1/2} + }\\ {\underbrace {( - \frac{1}{{{\gamma _1}}} \cdot k_{v1}^ - \cdot |{\mathit{\boldsymbol{e}}_v}| \cdot {\rm{sign}} (|{\mathit{\boldsymbol{e}}_v}| - \varepsilon ) + {\beta _1}) \cdot |{k_{v1}} - k_{v1}^*|}_{{\zeta _1}} + } \end{array}\\ \underbrace {( - \frac{1}{{{\gamma _2}}} \cdot \eta \cdot k_{v1}^ - \cdot |{\mathit{\boldsymbol{e}}_v}| \cdot {\rm{ sign}} (|{\mathit{\boldsymbol{e}}_v}| - \varepsilon ) + {\beta _2}) \cdot |{k_{v2}} - k_{v2}^*|}_{{\zeta _2}} \end{array} $

（35）

式中:β_η=min{η₁, β₁, β₂}。

当|e_v|>ε₂时，k_v1和k_v2分别以k_v1·|e_v|和η_vk_v1·|e_v|的斜率增长，若满足式(36)中所示的条件时，可得ζ₁>0，ζ₂>0。

$ {\gamma _1} < \frac{{{{\bar k}_{v1}} \cdot \varepsilon }}{{{\beta _1}}},{\gamma _2} < \frac{{\eta \cdot {{\bar k}_{v1}} \cdot \varepsilon }}{{{\beta _2}}} $	（36）

则$ \dot{V}_{2} \leqslant-\beta_{\eta} V_{2}^{1 / 2}+\zeta_{1}+\zeta_{2} \leqslant-\beta_{\eta} V_{2}^{1 / 2}$，即|e_v|可在有限时间内收敛到|e_v|>ε₂区间内；当|e_v| < ε₂时，ζ₁、ζ₂小于零，$\dot{V}_{2} $的正负不确定，增益变化率会变为-k_v1·|e_v|和-η_vk_v1·|e_v|，即增益会逐渐减小直至|e_v|>ε₂后，再以k_v1·|e_v|和η_vk_v1·|e_v|的斜率增长。

3 仿真分析 3.1 仿真参数设定

仿真中，采用文献[8]中提供的固定翼无人机的六自由度模型以及气动数据，转动惯量I为

$ \mathit{\boldsymbol{I}} = \left[ {\begin{array}{*{20}{c}} {0.552{\kern 1pt} {\kern 1pt} {\kern 1pt} 8}&0&{0.001{\kern 1pt} {\kern 1pt} {\kern 1pt} 5}\\ 0&{0.633{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 5}&0\\ {0.001{\kern 1pt} {\kern 1pt} {\kern 1pt} 5}&0&{1.078{\kern 1pt} {\kern 1pt} {\kern 1pt} 3} \end{array}} \right] $

姿态控制器式(12)和式(13)和速度控制器式(21)~式(23)的控制参数如表 1所示，仿真周期为1 ms。

设定姿态干扰d_ω=[40 40 40]^T·sin(3πt/19)，速度环的干扰为Δd_v=30sin(4πt/19)，姿态角的初值为Θ₀=[0.6° 1.1°1.7°]，姿态角参考值为Θ_d=[40°51°69°]，角速率的初值为ω₀=[0.7 0.750.8](°)/s，空速初值为V₀=4 m/s，空速参考值为V_d=12 m/s。为了验证本文所提算法的有效性和优越性，将本文所提的算法分别与单变量螺旋算法和文献[16]中的自适应超螺旋算法进行对比。控制参数如下所示：姿态环的单变量螺旋滑模：k_1ϕ=k_1θ=k_1ψ=45；速度环的自适应超螺旋滑模：σ_T=0.001，ω₁=5，r₁=5，ε₁^[16]=2。固定翼无人机控制力矩上界为M_max=200 N·m。

3.2 仿真结果

图 3和图 4分别为自适应多变量螺旋控制器和传统的固定增益螺旋控制器(在图中下标为fixed)的姿态角误差曲线以及姿态角的跟踪曲线，由图 3可以看出，在姿态控制器式(12)和式(13)的作用下，姿态角误差可以在有限时间(< 1.5 s)内收敛到0附近的有限邻域(-0.002°~0.002°)内，与传统螺旋控制器相比，本文所提出的自适应多变量超螺旋控制器精度更高，收敛速度更快，超调更小。

图 5为力矩变化的曲线，由图中曲线知，采用自适应多变量螺旋滑模产生的力矩M相比于M_fixed抖振更小，且力矩最大值≤M_max。

图 6和图 7分别为在自适应快速超螺旋和自适应超螺旋^[16]控制下的速度误差曲线和速度跟踪曲线，由图 6可以看出，在速度控制器式(21)~式(23)的作用下，空速V在1 s内收敛，并且超调很小，空速稳定时的收敛误差在10^-3级，相比于自适应超螺旋，本文所提的算法具有更好的控制性能。

图 8为姿态控制器控制增益和速度控制器控制增益的变化曲线，当姿态角跟踪误差小于ε₁=0.002时，增益保持在30不变；当速度误差e_v小于ε₂时，控制增益以较小的斜率减小。

图 9为推力变化曲线，由图可以看出，推力在开始的1 s内变化较为剧烈，待无人机状态稳定后，推力变化平缓，根据外界扰动，实时变化。而在自适应超螺旋的作用下，后期稳定时，推力的变化更为剧烈。

4 结论

本文考虑固定翼无人机飞行过程中的强干扰、不确定以及执行器饱和等问题，基于自适应多变量螺旋滑模和自适应快速超螺旋滑模设计了姿态控制器和速度控制器，有效地解决了干扰以及不确定上界未知时的控制问题，有效地抑制了抖振，简化了控制结构，对参考指令实现快速稳定精确跟踪。在理论上，能较好地适用于固定翼无人机的飞行控制需求。