目前，无人系统已替代人类执行特定的任务，尤其无人机在军事领域取得的效果显著，将成为未来战争的核心武器。因此，它的发展前景无法估量。多无人机编队是无人系统的显著代表，在军用和民用领域发挥着至关重要的作用，也是当前研究的热点^[1-2]。在执行作战任务过程中，无人机编队队形保持能够很大程度提高作战效率，比如：协同侦察、追踪保持^[3-4]、感知识别以及协同打击等^[5-6]；在战场环境发生变化时，无人机编队则必须进行编队重构，才能保证编队的生存力^[7]；对于局部战争，携有无人机机群的大飞机在指定作战区域上方释放，队形控制极为重要，为任务分配和目标打击奠定基础。因此，无人机编队队形控制方法是实现无人机编队安全飞行的前提和必要的手段，其有效性、合理性以及科学性直接影响作战任务的成败。针对无人机编队队形控制的研究方法非常多，而且也非常丰富，比如：领航跟随法^[8-9]、虚拟结构法^[10]、基于行为法^[11]和图论法^[12]等，这些丰硕的理论成果虽发表在10年前，但解决了无人机近距编队控制问题，即队形保持控制。然而，对于舰载多无人机编队对海突击作战队形保持控制，却因其具有前沿性、复杂性等特点，导致这方面的研究较少。因此，多无人机在执行作战任务过程中，快速集结、形成期望的队形以及保持编队一致性等关键技术是科研工作者亟需解决的问题，尤其队形保持稳定。

无人机编队控制问题，已有大量相关领域的专家学者从理论到实验对其进行深入的研究。文献[13-16]，优化控制方法、图论法、导引方法和人工势场方法用于解决编队控制问题。文献[17]基于最优控制的航迹规划算法，设计了一种新的目标函数进行组合优化处理，使无人机搜索过程达到全局最优。文献[18]作者从队形保持和队形重构进行展开研究，一方面设计了一种分布式反馈控制器，采用虚拟仿真平台模拟3架无人机加速和转弯的场景，并验证所设计的控制器的有效性；另一方面，将队形重构问题转化为燃料最优控制问题，并采用最小安全距离和最大通信距离进行队形重构控制。文献[19]，作者基于“长机-僚机”结构的编队控制方法，结合军队急行军中同列依次替补和末排内向收拢原则，提出了机器人近距队形保持策略。文献[20]中，针对无人机编队感知能力有限和队形控制方法的不足，提出了一种基于规则队形控制方法，其核心思想：将队形控制问题转化为无人机追踪自身期望位置。然而，在实际飞行中，编队中的任意一架无人机出现故障或者战毁，队形会出现位置空缺，影响整个编队队形的紧凑性和稳定性。文献[21]中，作者针对无人机编队队形保持和队形重构问题，设计了一种分布式反馈优化控制器，对三角形和线型队形的加速度和转弯情景进行飞行控制仿真，并通过Lyapunov函数方法证明所设计控制器的有效性。但是，该文献未考虑有风场或者障碍物干扰的情形，因此，所设计的控制器具有局限性。针对无人机编队队形控制，已有的控制方法通常具有很高的计算成本而且很难实现局部子系统的控制。虽然在文献[18]中，作者也同样设计了一种分布式控制器，但由于实时数据计算量大，控制器的精度低，很难实现期望的队形。然而，Stipanovic等^[22]针对每架无人机设计了一种反馈控制器，并采用分布重叠控制技术，将分散的无人机队形变换为期望的队形。与此同时，茹常剑等^[23]采用模型预测方法解决编队队形控制问题，其基本思路都是采用滚动优化算法。

在文献[24]，针对“欧拉-拉格朗日(EL)”动力学描述的不确定非均匀的非线性智能体模型，提出了一种自适应分层编队队形控制方法。将编队队形控制问题转化为同步问题研究，采用分布式模型参考自适应控制对EL系统进行同步研究。其核心思想：每个智能体均收敛到由其分层次优越的邻居定义的模型。然而，该文献的作者未考虑EL系统动力学的约束条件和执行器饱和等因素。在同时期，文献[25]采用滚动时域控制方法研究复杂无人机编队队形控制。文献[26]基于分层机制和模型预测控制方法，提出了一种多无人机编队控制方法。这些方法，控制无人机编队队形的精度高，但是实时性差。Arcak等^[27]首次提出的Back-stepping控制方法，该方法用于解决非线性编队模型，其基本思想是首先设计子系统期望的虚拟输入信号，然后反馈给控制器得到真实子系统的控制输入；该方法常用于设计无人机编队控制器，通过构建Lyapunov函数证明其有效性。与线性反馈控制方法相比较，Back-stepping方法具有较强的灵活性，且不要求系统必须是线性闭环的。文献[28-29]将“长机-僚机”控制策略应用无人机编队控制中，将任意一架无人机作为长机，其余无人机作为僚机，简化为多机编队控制问题。然而，这种方法具有局限性。文献[30]，为解决模拟战场环境中战机编队队形在机动过程中避障后到达目标点问题，结合基于行为法思想，对跟随领航算法进行改进，并对战机集群实施控制。但是，战机集群仅仅有一个虚拟领航者，一旦战毁，势必会影响战场的态势。西弗吉尼亚大学研发团队^[31]提出虚拟“长机”方法，他们将无人机看作一个刚体，每架无人机保持既定的队形跟踪固定的点。Brad等^[32]采用基于行为控制方法，提出了编队队形几何中心的概念。在文献[33]，针对多无人机编队飞行过程队形保持问题，作者提出了一种二阶非线性切换拓扑结构的领航跟随一致性控制协议，该协议在速度一致性情况下可实现的编队队形协同飞行。然而，作者未考虑无人机之间信息交互过程中存在时延问题，且控制效果不明显，工程应用价值有限。因此，本文提出了一种Back-stepping控制方法，能够解决无人机编队队形控制精度低、时效性差的问题。

多无人机在集结期望队形过程中，任意一架无人机出现故障，使得整个编队的通讯中断，无人机出现不可控的局面，极容易发生事故。因此，任意一架无人机的姿态和位置的控制将直接影响整个编队的控制，进而影响集结效率。基于此，本文采用Back-stepping控制方法，解决无人机编队重构和达到稳态所需时间的问题。

本文考虑无人机以不同航向角从不同位置起飞到集结期望队形过程中，因单机出现故障引起队形紊乱，致使整个无人机编队系统的通讯网络破坏，队形无法保持的问题，利用反步推演法对每架无人机模型进行虚拟控制律设计，以此类推，可获得整个无人机编队闭环系统的实际控制律(协同导引控制律)，用于解决无人机编队队形重构和快速达到期望的队形。同时，结合Lyapunov稳定性分析来保证闭环系统的收敛性和鲁棒性。

1 问题描述

假定有N(N>20)架固定翼无人机组成一个正多边形的多无人机编队，其中有一架虚拟长机位于正多边形的几何中心，其余N-1架无人机位于多边形的顶点处。每架无人机从不同起点起飞，追踪虚拟长机，并迅速集合编队，然后保持既定的队形编队飞行。在该过程中，建立恰当的无人机运动模型是控制无人机编队队形的关键因素之一，能够准确描述无人机之间的通信关系。多无人机编队从松散队形到预设队形直到队形保持过程中，历经外界扰动阶段，队形调整阶段以及队形保持阶段，如图 1所示。

1) 外界扰动阶段：编队中的虚拟长机携带无线通讯设备，充当预警机角色，与其余僚机进行信息交互，僚机跟随虚拟长机飞行。与此同时，受到三维空间风场扰动V_风，无人机队形处于杂乱无章的状态，此过程中，无人机编队处于外界(风场)扰动阶段。

2) 队形调整阶段：无人机编队中每架无人机在Back-stepping控制方法的作用下调整各自的位置和姿态，使得与邻近的无人机之间的相对距离误差在3个方向趋于零，保持稳定状态，此过程中，无人机编队处于队形调整阶段。

3) 队形保持阶段：无人机编队队形调整成正三角形编队后，任意两架无人机之间的相对距离a、b、c保持不变，按照此队形稳定飞行，此过程，无人机编队队形处于队形保持阶段。

2 建立无人机模型 2.1 无人机模型

在惯性坐标系下，假设无人机在惯性坐标系下的速度的v_g方向与机体轴重合，并且推力和阻力方向共线，如图 2所示。本文采用无人机三自由度的简化模型，且速度矢量、偏航角以及侧滑角都是一阶动力学模型^[34-37]，

$ \left\{ \begin{array}{l} \dot x = {v_{\rm{g}}}\cos \gamma \cos \chi \\ \dot y = {v_{\rm{g}}}\cos \gamma \sin \chi \\ \dot z = {v_g}\sin \chi \\ {{\dot v}_{\rm{g}}} = {c_1}\left( {v_{\rm{g}}^c - {v_{\rm{g}}}} \right)\\ \dot \gamma = {c_2}\left( {{\gamma ^c} - \gamma } \right)\\ \dot \chi = {c_3}\left( {{\chi ^c} - \chi } \right) \end{array} \right. $	（1）

式中：x、y、z表示无人机在惯性坐标系下的位置；v_g、γ、χ分别表示在惯性坐标系下的速度、俯仰角以及偏航角；c₁、c₂、c₃均表示时间常数；v_g^c、χ^c、γ^c分别表示速度指令、俯仰角指令以及偏航角指令。图 2中，V表示无人机在机体坐标系下的速度。

文献[38]采用虚拟长机的5个状态变量的一阶运动方程模型，研究无人机的路径规划。目前，无人机的一阶动力学模型仍然被学者广泛采用，因为一阶动力学模型是研究高阶的基础。基于此，虚拟长机的运动学模型为

$ \left\{ {\begin{array}{*{20}{l}} {x = {v_{\rm{g}}}\cos \gamma \cos \chi }\\ {\dot y = {v_{\rm{g}}}\cos \gamma \sin \chi }\\ {\dot z = {v_{\rm{g}}}\sin \chi }\\ {{{\dot v}_{\rm{g}}} = {c_1}\left( {v_{\rm{g}}^c - {v_{\rm{g}}}} \right)}\\ {\dot \gamma = \frac{g}{{{v_{\rm{g}}}}}\phi \eta \xi }\\ {\dot \chi = \frac{{g\cos \gamma }}{{{v_{\rm{g}}}}}\xi } \end{array}} \right. $	（2）

式中：ϕ、η、ξ分别表示虚拟长机滚转角、过载系数和侧风强度等级。根据几何关系，式(2)可简化为

$ \left\{ {\begin{array}{*{20}{l}} {\dot x = {v_{\rm{g}}}\cos \gamma \cos \chi }\\ {\dot y = {v_{\rm{g}}}\cos \gamma \sin \chi }\\ {\dot z = {v_{\rm{g}}}\sin \chi }\\ {{{\dot v}_{\rm{g}}} = {c_1}\left( {v_{\rm{g}}^c - {v_{\rm{g}}}} \right)}\\ {\dot \gamma = {a_{\rm{h}}}/{v_{\rm{g}}}}\\ {\dot \chi = {a_{\rm{v}}}/{v_{\rm{g}}} \cdot \cos \gamma } \end{array}} \right. $	（3）

式中：a_h=gϕηξ，a_v=gξ，a_v、a_h分别表示虚拟长机在横向和纵向的加速度。

2.2 无人机模型

为了用数学化语言描述该网络拓扑结构的模型，本文使用图论理论^[39]。本文中使用图G=(W, S)建立N架无人机之间信息交互模型，这里的W={w₁, w₂, …, w_N}表示有序点集，S∈W×W表示有序边集。图中的边(w_i, w_j)表示无人机i能够直接把信息传递给无人机j的有向路径，即无人机j能够直接获取并使用无人机i的信息。

假设虚拟长机嵌入到网络拓扑结构的子系统，虚拟长机与周围临近无人机之间进行信息共享，且至少有一个点(无人机)通讯是正常的，这样能够保证信息的共享和交互。网络拓扑结构的连接用矩阵A来表示，该矩阵的对角元素a_ii>0是追踪的增益。

3 控制器设计

本文提出的无人机编队模型，是由N-1(N>20)架无人机和一架虚拟长机构成。对于每架无人机期望的动力学模型，

$ \left\{ {\begin{array}{*{20}{l}} {{{\dot V}_{xi}} = a_{xi}^{\rm{d}}\cos \gamma _i^{\rm{d}}\cos \chi _i^{\rm{d}}}\\ {{{\dot V}_{yi}} = a_{yi}^{\rm{d}}\cos \gamma _i^{\rm{d}}\sin \chi _i^{\rm{d}}}\\ {{{\dot V}_{zi}} = a_{zi}^{\rm{d}}\sin \chi _i^{\rm{d}}} \end{array}} \right. $	（4）

式中：i=1, 2, …, N; a_xi^d、γ_i^d、χ_i^d分别表示期望的加速度指令、期望的俯仰角指令以及期望的偏航角指令；V_xi、V_yi、V_zi分别表示无人机速度在惯性坐标系下3个方向的速度大小。

对于任意一架无人机, 定义追踪误差为

$ \left\{ {\begin{array}{*{20}{l}} {{e_{xi}} = \sum\limits_{i,j = 1}^N {{a_{ij}}} \left( {{x_j} - {x_i}} \right) + k\left( {{x_{\rm{d}}} - {x_i}} \right)}\\ {{e_{yi}} = \sum\limits_{i,j = 1}^N {{a_{ij}}} \left( {{y_j} - {y_i}} \right) + k\left( {{y_{\rm{d}}} - {y_i}} \right)}\\ {{e_{zi}} = \sum\limits_{i,j = 1}^N {{a_{ij}}} \left( {{z_j} - {z_i}} \right) + k\left( {{z_{\rm{d}}} - {z_i}} \right)} \end{array}} \right. $

（5）

式中：i=1, 2, …, N；k为正系数；x_d、y_d、z_d分别表示无人机的所在的期望位置；e_xi、e_yi、e_zi分别表示无人机在空间3个方向的追踪误差。

基于矩阵论，式(5)可以简化为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{e}}_x} = - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)\left( {\mathit{\boldsymbol{x}} - {x_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)}\\ {{\mathit{\boldsymbol{e}}_y} = - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)\left( {\mathit{\boldsymbol{y}} - {y_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)}\\ {{\mathit{\boldsymbol{e}}_z} = - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)\left( {\mathit{\boldsymbol{z}} - {z_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)} \end{array}} \right. $

（6）

式中：L表示i×i阶拉普拉斯矩阵；K表示i×i阶正常数矩阵；I_i表示i×i阶的单位矩阵。

对式(6)时间求导，可得到进一步的简化后的误差动力学模型，

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{x}}} = - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)\left( {\mathit{\boldsymbol{\dot x}} - {{\dot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)}\\ {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{y}}} = - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)\left( {\mathit{\boldsymbol{\dot y}} - {{\dot y}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)}\\ {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{z}}} = - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)\left( {\mathit{\boldsymbol{\dot z}} - {{\dot z}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)} \end{array}} \right. $

（7）

由式(7)可知，动力学误差模型分别由3个位置通道复合而成。令$\mathit{\boldsymbol{\dot x}} = {\mathit{\boldsymbol{u}}_x}, \mathit{\boldsymbol{\dot y}} = {\mathit{\boldsymbol{u}}_y}, \mathit{\boldsymbol{\dot z}} = {\mathit{\boldsymbol{u}}_z}$，则式(7)可简化为

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{x}}} = - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)\left( {{\mathit{\boldsymbol{u}}_x} - {{\dot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)}\\ {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{y}}} = - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)\left( {{\mathit{\boldsymbol{u}}_y} - {{\dot y}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)}\\ {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{z}}} = - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)\left( {{\mathit{\boldsymbol{u}}_z} - {{\dot z}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)} \end{array}} \right. $

（8）

解得:

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}_x} = {{\left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}{{\mathit{\boldsymbol{\dot e}}}_x} + {{\dot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}}\\ {{\mathit{\boldsymbol{u}}_y} = {{\left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}{{\mathit{\boldsymbol{\dot e}}}_y} + {{\dot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}}\\ {{\mathit{\boldsymbol{u}}_z} = {{\left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}{{\mathit{\boldsymbol{\dot e}}}_z} + {{\dot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \end{array}} \right. $

（9）

令P=(L+K)^-1，式(9)可简化为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}_x} = \mathit{\boldsymbol{P}}{{\mathit{\boldsymbol{\dot e}}}_x} + {{\dot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}}\\ {{\mathit{\boldsymbol{u}}_y} = \mathit{\boldsymbol{P}}{{\mathit{\boldsymbol{\dot e}}}_y} + {{\dot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}}\\ {{\mathit{\boldsymbol{u}}_z} = \mathit{\boldsymbol{P}}{{\mathit{\boldsymbol{\dot e}}}_z} + {{\dot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \end{array}} \right. $

（10）

假如图G有向的，可以推断出矩阵(L+K)是正定的，易得矩阵P是正定的。

此处，定义误差期望的值为：

$ \left\{ {\begin{array}{*{20}{l}} {v_{\rm{g}}^{\rm{d}} = \sqrt {{{\left( {{u_{xi}}} \right)}^2} + {{\left( {{u_{yi}}} \right)}^2} + {{\left( {{u_{zi}}} \right)}^2}} }\\ {\gamma _i^{\rm{d}} = \arctan \frac{{{u_{zi}}}}{{\sqrt {{{\left( {{u_{xi}}} \right)}^2} + {{\left( {{u_{yi}}} \right)}^2}} }}}\\ {\chi _i^{\rm{d}} = \arctan \frac{{{u_{yi}}}}{{\sqrt {{{\left( {{u_{xi}}} \right)}^2} + {{\left( {{u_{yi}}} \right)}^2}} }}} \end{array}} \right. $

（11）

同时，定义e_vg，e_γ，e_χ的误差为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}} = {\mathit{\boldsymbol{v}}_{\rm{g}}} - \mathit{\boldsymbol{v}}_{\rm{g}}^{\rm{d}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}} = \mathit{\boldsymbol{\gamma }} - {\mathit{\boldsymbol{\gamma }}^{\rm{d}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}} = \mathit{\boldsymbol{\chi }} - {\mathit{\boldsymbol{\chi }}^{\rm{d}}}} \end{array}} \right. $

（12）

式中：v_g=[v_g1, v_g2, …, v_gN]^T，γ=[γ₁, γ₂, …, γ_N]^T，χ=[χ₁, χ₂, …, χ_N]^T。

欲保持多无人机按照期望的队形飞行，则e_vg、e_γ、e_χ均大于零。为保持所设计的控制器收敛，误差的导数需小于零，有以下关系：

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot e}}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}} = {\lambda _{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}}\\ {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{\gamma }}} = {\lambda _\mathit{\boldsymbol{\gamma }}}{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}}}\\ {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{\chi }}} = {\lambda _\mathit{\boldsymbol{\chi }}}{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}} \end{array}} \right. $

（13）

式中：λ_vg, λ_γ, λ_χ均小于零。

对式(11)中v_g^d、γ_i^d、χ_i^d分别求导，

$ \left\{ \begin{array}{l} \dot v_{\rm{g}}^{\rm{d}} = \frac{{{{\dot u}_{xi}}{u_{xi}} + {{\dot u}_{yi}}{u_{yi}} + {{\dot u}_{zi}}{u_{zi}}}}{{\sqrt {{{\left( {{u_{xi}}} \right)}^2} + {{\left( {{u_{yi}}} \right)}^2} + {{\left( {{u_{zi}}} \right)}^2}} }}\\ \dot \gamma _i^{\rm{d}} = \frac{{{{\dot u}_{zi}}\left[ {{{\left( {{u_{xi}}} \right)}^2} + {{\left( {{u_{yi}}} \right)}^2}} \right] - {u_{zi}}\left( {{{\dot u}_{xi}}{u_{xi}} + {{\dot u}_{yi}}{u_{yi}}} \right)}}{{{{\left[ {{{\left( {{u_{xi}}} \right)}^2} + {{\left( {{u_{yi}}} \right)}^2} + {{\left( {{u_{zi}}} \right)}^2}} \right]}^2} \cdot \sqrt {{{\left( {{u_{xi}}} \right)}^2} + {{\left( {{u_{yi}}} \right)}^2}} }}\\ \dot \chi _i^{\rm{d}} = \frac{{{{\dot u}_{yi}}{{\left( {{u_{xi}}} \right)}^2} - {{\dot u}_{xi}}{u_{xi}}{u_{yi}}}}{{{{\left[ {{{\left( {{u_{xi}}} \right)}^2} + {{\left( {{u_{yi}}} \right)}^2} + {{\left( {{u_{zi}}} \right)}^2}} \right]}^2} \cdot \sqrt {{{\left( {{u_{xi}}} \right)}^2} + {{\left( {{u_{yi}}} \right)}^2}} }} \end{array} \right. $

（14）

由式(14)可知，v_g^d、γ_i^d、χ_i^d的导数大小取决于闭环系统中位置信号的数值，即u_x、u_y、u_z。

对式(9)时间求导，可得

$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\dot u}}}_\mathit{\boldsymbol{x}}} = \mathit{\boldsymbol{P}}{{\mathit{\boldsymbol{\ddot e}}}_\mathit{\boldsymbol{x}}} + {{\ddot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}\\ {{\mathit{\boldsymbol{\dot u}}}_\mathit{\boldsymbol{y}}} = \mathit{\boldsymbol{P}}{{\mathit{\boldsymbol{\ddot e}}}_\mathit{\boldsymbol{y}}} + {{\ddot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}\\ {{\mathit{\boldsymbol{\dot u}}}_\mathit{\boldsymbol{z}}} = \mathit{\boldsymbol{P}}{{\mathit{\boldsymbol{\ddot e}}}_\mathit{\boldsymbol{z}}} + {{\ddot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i} \end{array} \right. $

（15）

此处，虚拟长机位置追踪误差的控制需得到u_x、u_y、u_z。为实现无人机快速追踪虚拟长机，需要得到$\left({{{\mathit{\boldsymbol{\dot u}}}_x}, {{\mathit{\boldsymbol{\dot u}}}_y}, {{\mathit{\boldsymbol{\dot u}}}_z}} \right)$。根据虚拟长机的位置，可以得到无人机追踪的参考轨迹。假如无人机的期望的速度v，俯仰角γ以及偏航角χ是正常数，可得到式(13)另一种形式，

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot e}}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}} = {\lambda _{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}} + \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{v}}_{\rm{g}}}{\mathit{\boldsymbol{I}}_i}}\\ {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{\gamma }}} = {\lambda _\gamma }{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}} + \mathit{\boldsymbol{B\gamma }}{\mathit{\boldsymbol{I}}_i}}\\ {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{\chi }}} = {\lambda _\chi }{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}} + \mathit{\boldsymbol{B\chi }}{\mathit{\boldsymbol{I}}_i}} \end{array}} \right. $

（16）

其中：式(18)中的矩阵B=P^－1；由于矩阵P为正定，易知矩阵B也是正定的。

联立式(1)和式(12)，整理可得

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{v}}_{\rm{g}}^c = {\mathit{\boldsymbol{v}}_{\rm{g}}} + \frac{1}{{{c_1}}}\left( {{{\mathit{\boldsymbol{\dot e}}}_{{v_{\rm{g}}}}} + \mathit{\boldsymbol{\dot v}}_{\rm{g}}^{\rm{d}}} \right)}\\ {{\mathit{\boldsymbol{\gamma }}^c} = \gamma + \frac{1}{{{c_2}}}\left( {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{\gamma }}} + {{\mathit{\boldsymbol{\dot \gamma }}}^{\mathop {\rm{d}}\limits^. }}} \right)}\\ {{\mathit{\boldsymbol{\chi }}^c} = \chi + \frac{1}{{{c_3}}}\left( {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{\chi }}} + {{\mathit{\boldsymbol{\dot \chi }}}^{\rm{d}}}} \right)} \end{array}} \right. $

（17）

将式(16)代入式(17)，整理可得

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{v}}_{\rm{g}}^c = {\mathit{\boldsymbol{v}}_{\rm{g}}} + \frac{1}{{{c_1}}}\left( {{\lambda _{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}{{\mathit{\boldsymbol{\dot e}}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}} + \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{v}}_{\rm{g}}}{\mathit{\boldsymbol{I}}_i} + \mathit{\boldsymbol{\dot v}}_{\rm{g}}^{\rm{d}}} \right)}\\ {{\mathit{\boldsymbol{\gamma }}^c} = \gamma + \frac{1}{{{c_2}}}\left( {{\lambda _\mathit{\boldsymbol{\gamma }}}{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{\gamma }}} + \mathit{\boldsymbol{B\gamma }}{\mathit{\boldsymbol{I}}_i} + \mathit{\boldsymbol{\dot v}}_{\rm{g}}^{\rm{d}}} \right)}\\ {{\mathit{\boldsymbol{\chi }}^c} = \mathit{\boldsymbol{\chi }} + \frac{1}{{{c_3}}}\left( {{\lambda _\mathit{\boldsymbol{\chi }}}{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}} + \mathit{\boldsymbol{B\chi }}{\mathit{\boldsymbol{I}}_i} + {{\mathit{\boldsymbol{\dot \chi }}}^{\rm{d}}}} \right)} \end{array}} \right. $

（18）

对式(18)进行等价变化，可得协方差控制律，

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{v}}_{\rm{g}}^{ci} = {\mathit{\boldsymbol{v}}_{gi}} + \frac{1}{{{c_{1i}}}}\left( {{\lambda _{{\mathit{\boldsymbol{v}}_{{\rm{g}}i}}}}{{\mathit{\boldsymbol{\dot e}}}_{{\mathit{\boldsymbol{v}}_{{\rm{g}}i}}}} + \mathit{\boldsymbol{B}}{v_{{\rm{g}}i}}{\mathit{\boldsymbol{I}}_i} + \dot v_{{\rm{g}}i}^{\rm{d}}} \right)}\\ {{\mathit{\boldsymbol{\gamma }}^{ci}} = {\mathit{\boldsymbol{\gamma }}_i} + \frac{1}{{{c_{2i}}}}\left( {{\lambda _{\gamma i}}{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{\gamma }}} + \mathit{\boldsymbol{B}}{\gamma _i}{\mathit{\boldsymbol{I}}_i} + \dot \gamma _i^{\rm{d}}} \right)}\\ {{\mathit{\boldsymbol{\chi }}^{ci}} = {\mathit{\boldsymbol{\chi }}_i} + \frac{1}{{{c_{3i}}}}\left( {{\lambda _{\mathit{\boldsymbol{\chi }}i}}{\mathit{\boldsymbol{e}}_{\mathit{\boldsymbol{\chi }}i}} + \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{\chi }}_i}{\mathit{\boldsymbol{I}}_i} + \dot \chi _i^{\rm{d}}} \right)} \end{array}} \right. $

（19）

式中：i=1, 2, …, N。

4 稳定性分析

为验证所设计无人机编队队形控制器的稳定性，构建Lyapunov函数，

$ \begin{array}{l} V = \frac{1}{2}\left[ {{{\left\| {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}} \right\|}^2} + {{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}}} \right\|}^2} + {{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}} \right\|}^2} + {{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}} \right\|}^2} + } \right.\\ \;\;\;\;\;\;\left. {{{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \right\|}^2} + {{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}} \right\|}^2}} \right] \end{array} $

（20）

由式(20)，可知V≥0。对式(20)中V的时间求导，可得

$ \begin{array}{l} \dot V = \frac{1}{2}\left[ {{{\left\| {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}} \right\|}^2} + {{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}}} \right\|}^2} + {{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}} \right\|}^2} + {{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}} \right\|}^2} + } \right.\\ \;\;\;\;\;\;\left. {{{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \right\|}^2} + {{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}} \right\|}^2}} \right] \end{array} $

（21）

将式(6)、式(7)、式(12)以及(16)代入式(21)，可得

$ \begin{array}{l} \dot V = {\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}{{\mathit{\boldsymbol{\dot e}}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}} + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}}{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{\gamma }}} + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{\chi }}} + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{x}}} + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{y}}} + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{z}}} = \\ \;\;\;\;\;{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}\left( {{\lambda _{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}} + \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{v}}_{\rm{g}}}{\mathit{\boldsymbol{I}}_i}} \right) + {\mathit{\boldsymbol{e}}_\gamma }\left( {{\lambda _\mathit{\boldsymbol{\gamma }}}{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}} + \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{\gamma }}_r}{\mathit{\boldsymbol{I}}_i}} \right) + \\ \;\;\;\;\;{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}\left( {{\lambda _\mathit{\boldsymbol{\chi }}}{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}} + \mathit{\boldsymbol{B}}{\chi _r}{\mathit{\boldsymbol{I}}_i}} \right) + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}\left[ { - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right) \cdot } \right.\\ \;\;\;\;\;\left. {\left( {{\mathit{\boldsymbol{u}}_x} - {{\dot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)} \right] + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}\left[ { - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right) \cdot \left( {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{y}}} - {{\dot y}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)} \right] + \\ \;\;\;\;\;{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}\left[ { - \left( {\mathit{\boldsymbol{L}} + \mathit{\boldsymbol{K}}} \right)\left( {{\mathit{\boldsymbol{u}}_z} - {{\dot z}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right)} \right] \end{array} $

（22）

根据图论，如果无人机编队之间的通讯链路是单向的，则矩阵B是正定的。假设矩阵B是奇异矩阵，其行列式的最大值和最小值分别为ζ和ζ；(L+K)也是奇异矩阵，其行列式的最大值和最小值分别为ω和ω；矩阵[λ_vg, λ_γ, λ_χ]的下界为[λ_vg, λ_γ, λ_χ]，可得

$ \begin{array}{l} \dot V \le {\underline \lambda _{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}{\left\| {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}} \right\|^2} + \bar \zeta \left\| {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}} \right\|\left\| {{\mathit{\boldsymbol{v}}_{\rm{g}}}{\mathit{\boldsymbol{I}}_i}} \right\| + {\underline \lambda _\mathit{\boldsymbol{\gamma }}}{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}}} \right\|^2} + \\ \;\;\;\;\;\;\bar \zeta \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}}} \right\|\left\| {\mathit{\boldsymbol{\gamma }}{\mathit{\boldsymbol{I}}_i}} \right\| + {\underline \lambda _\mathit{\boldsymbol{\chi }}}{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}} \right\|^2} + \bar \zeta \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}} \right\|\left\| {\mathit{\boldsymbol{\chi }}{\mathit{\boldsymbol{I}}_i}} \right\| - \\ \;\;\;\;\;\;\underline {\zeta \omega } \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}} \right\|\left\| {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{x}}}} \right\| + \underline \zeta \bar \omega \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}} \right\|\left\| {{{\dot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right\| - \underline {\zeta \omega } \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \right\|\left\| {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{y}}}} \right\| + \\ \;\;\;\;\;\;\underline \zeta \bar \omega \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \right\|\left\| {{{\dot y}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right\| - \underline {\zeta \omega } \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}} \right\|\left\| {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{z}}}} \right\| + \underline \zeta \bar \omega \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}} \right\|\left\| {{{\dot z}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right\| \end{array} $

（23）

由于无人机编队飞行轨迹是有界的，易知${\left\| {{\mathit{\boldsymbol{v}}_{\rm{g}}}{\mathit{\boldsymbol{I}}_i}} \right\| \le \max {\mathit{\boldsymbol{v}}_{\rm{g}}}, \left\| {\mathit{\boldsymbol{\gamma }}{\mathit{\boldsymbol{I}}_i}} \right\| \le \max \mathit{\boldsymbol{\gamma }}, \left\| {\mathit{\boldsymbol{\chi }}{\mathit{\boldsymbol{I}}_i}} \right\| \le \max \mathit{\boldsymbol{\chi, }}}$${\left\| {{{\dot x}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right\| \le \max {{\dot x}_{\rm{d}}}, \left\| {{{\dot y}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right\| \le \max {{\dot y}_{\rm{d}}}, \left\| {{{\mathit{\dot z}}_{\rm{d}}}{\mathit{\boldsymbol{I}}_i}} \right\| \le \max \; {{\mathit{\dot z}}_{\rm{d}}}}$。

于是，式(23)可进一步简化，

$ \begin{array}{l} \dot V \le {\underline \lambda _{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}{\left\| {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}} \right\|^2} + \bar \zeta \left\| {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}} \right\| + {\underline \lambda _\gamma }{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}}} \right\|^2} + \bar \zeta \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}}} \right\| + \\ {\underline \lambda _\mathit{\boldsymbol{\chi }}}{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}} \right\|^2} + \bar \zeta \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}} \right\| - \underline {\zeta \omega } \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}} \right\|\left\| {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{x}}}} \right\| + \underline \zeta \bar \omega \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}} \right\| - \\ \underline {\zeta \omega } \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \right\|\left\| {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{y}}}} \right\| + \underline \zeta \bar \omega \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \right\| - \underline {\zeta \omega } \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}} \right\|\left\| {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{z}}}} \right\| + \underline \zeta \bar \omega \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}} \right\| \end{array} $

（24）

根据柯西不等式定理，由式(10)可得以下不等式：

$ \left\{ \begin{array}{l} \left\| {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{x}}}} \right\| \le \left\| {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{x}}}} \right\| \le \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}} \right\|\\ \left\| {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{y}}}} \right\| \le \left\| {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{y}}}} \right\| \le \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \right\|\\ \left\| {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{z}}}} \right\| \le \left\| {{{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{z}}}} \right\| \le \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}} \right\| \end{array} \right. $

（25）

因此，式(24)进一步化简为：

$ \begin{array}{l} \dot V \le {\underline \lambda _{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}{\left\| {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}} \right\|^2} + \bar \zeta \left\| {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}} \right\| + {\underline \lambda _\gamma }{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}}} \right\|^2} + \bar \zeta \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}}} \right\| + \\ \;\;\;\;\;{\underline \lambda _\mathit{\boldsymbol{\chi }}}{\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}} \right\|^2} + \bar \zeta \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}} \right\| + \underline {\zeta \omega } {\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}} \right\|^2} + \underline \zeta \bar \omega \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}} \right\| + \\ \;\;\;\;\;\underline {\zeta \omega } {\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \right\|^2} + \underline \zeta \bar \omega \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \right\| + \underline {\zeta \omega } {\left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}} \right\|^2} + \underline \zeta \bar \omega \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}} \right\| \le \\ \;\;\;\;\;a{\left\| \mathit{\boldsymbol{s}} \right\|^2} + b\left\| \mathit{\boldsymbol{s}} \right\| \end{array} $

（26）

定义P=$ \operatorname{diag}\left(\underline{\lambda}_{v_{\mathrm{g}}}, \underline{\lambda}_{\gamma}, \underline{\lambda}_{\chi}, \underline{\zeta \omega}, \underline{\zeta \omega}, \underline{\zeta \omega}\right) $，a是P的最小特征值。${\mathit{\boldsymbol{s}} = {{\left[{\left\| {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{v}}_{\rm{g}}}}}} \right\|, \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\gamma }}}} \right\|, \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{\chi }}}} \right\|, \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}} \right\|, \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \right\|, \left\| {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{z}}}} \right\|} \right]}^{\rm{T}}}}$，b=max{max(ζv_g), max(ζγ), max(ζχ), ωe_x, ωe_y, ωe_z}。

于是式(26)可进一步简写为

$ \dot V \le 2aV + b\sqrt {2V} $	（27）

令方程$\dot{V}-2 a V-b \sqrt{2 V}=0$，其正数解为$\sqrt{V}=\sqrt{V(0)} \mathrm{e}^{(a t)}-\frac{b}{a}\left(1-\mathrm{e}^{(a t)}\right)$，则$\sqrt V \le \sqrt {V\left(0 \right)} - \frac{b}{{\sqrt 2 a}}$。

于是，$2 a V+b \sqrt{2 V} \leqslant \frac{b^{2}(1-2 a)}{2 a^{2}}<0$，可得$\dot{V} \leqslant 0$，证毕。

在设计无人机编队队形控制器过程中，有以下4点假设：

1) 如果无人机的期望位置(x_d, y_d, z_d)是常数且b=0，则图G是有向生成树，同时编队中的所有无人机均能实现一致追踪的目的。

2) 如果虚拟长机期望的参考位置是恒定的且追踪误差方程的二阶导数存在，则该方程是有界的。

3) 如果无人机编队系统包含具有生成树的有向图，则该系统是渐进稳定的。

4) 如果该系统的追踪误差方程组中至少有一个方程的二阶导数是有界的，则该系统的所有追踪误差也是一致有界的。基于以上的假设，本文所提出的控制方法是有效的，且能够实现无人机编队快速集结和队形保持。

5 仿真实验

为了验证无人机协同编队飞行队形保持控制器的有效性，基于MATLAB R2014a搭建Simulink仿真模块，并在同一台电脑进行虚拟仿真实验。本文采用了N(N>20)架无人机进行仿真实验，其中以无人机编队的最小单位3架无人机作为被控对象，然后扩展到N架。在5.1节，针对20架以上无飞机验证算法的有效性，即所提出的算法适用于无人机编队的规模；在5.2节和5.3节，20架无人机作为编队进行研究，由于他们具有相同的特性，如速度、姿态以及相对位置保持不变。在文中，最小被控对象由3架无人机和一架虚拟长机组成；虚拟长机携带雷达侦察设备，其余无人机携带作战武器。在集结过程中，所设计的控制器能够使得任意2架无人机之间相对距离误差收敛于零，即期望位置和实际位置重合，始终保持相对距离恒定，这样可以保证任意2架无人机避免发生碰撞。4架无人机编队的拓扑结构，如图 3所示。

网络拓扑结构图能够有效描述无人机之间的通讯关系，并也是多系统建模的工具。在4架无人机拓扑图中，虚拟长机嵌入该拓扑中，作为正三角形的几何中心。邻接矩阵A、度矩阵D、拉普拉斯矩阵L以及关联矩阵B，分别为

$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{l}} 0&1&1&1\\ 1&0&1&1\\ 1&1&0&1\\ 1&1&1&0 \end{array}} \right],\;\;\;\;\mathit{\boldsymbol{D}} = \left[ {\begin{array}{*{20}{l}} 3&0&0&0\\ 0&3&0&0\\ 0&0&3&0\\ 0&0&0&3 \end{array}} \right] $

$ \mathit{\boldsymbol{L}} = \left[ {\begin{array}{*{20}{c}} 3&{ - 1}&{ - 1}&{ - 1}\\ { - 1}&3&{ - 1}&{ - 1}\\ { - 1}&{ - 1}&3&{ - 1}\\ { - 1}&{ - 1}&{ - 1}&3 \end{array}} \right] $

$ \mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{l}} 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right] $

在仿真过程中，每架无人机的初始条件如表 1所示。其他初始条件：v_g^c(0)=12 m/s，a_h^c=a_v^c=0.25 m/s²，c₁=c₂=c₃=5，λ_vg=λ_γ=λ_χ=1。同时，无人机编队的协调转弯半径R取值范围(0≤R≤10 m)。从松散队形到集结队形直至期望队形，无人机编队的速度范围(0≤V≤47.5 m/s)，俯仰角θ变化范围(－11°≤θ≤11°)。

5.1 动力学特性分析

基于以上的假设和初始条件，无人机以不同的速度和偏航角从不同起点位置出发到集结编队，直至期望编队形稳定飞行。其动力特性如图 4~图 7所示。

图 4给出了无人机编队以不同速度从不同初始位置起飞到集结编队再到队形保持过程中速度变化曲线图。随着时间的变化，该曲线呈现一种先快速上升，之后缓慢上升直至稳态趋势。在10 s之前，曲线呈现快速上升。这是由于无人机以不同偏航角且快速追踪虚拟长机，由于虚拟长机按照预设的轨迹飞行，无人机不断修正姿态并保持速度方向始终朝向虚拟长机；在10~25 s，曲线缓慢上升。由于此阶段无人机处于集结状态，他们进行微调速度和姿态完成期望的编队；在25 s之后，无人机编队处于队形保持阶段，此阶段无人机编队保持速度和相对距离不变，实现稳态飞行。

图 5表示了无人机从初始位置到集结位置过程中转弯半径曲线图。随着时间的变化，曲线呈现一种先上升后缓慢上升，之后瞬间接近零的稳态趋势。根据无人机不发生侧翻或者横向滑移，无人机速度的平方与转弯半径成正比，则该曲线的变化趋势与图 4基本一致。无人机在10 s之前以不同的航向角加速追赶虚拟长机，由于速度增大，所需的转弯半径也相应增大，实现协同转弯；在10~25 s，曲线处于缓慢增大，这是由于无人机即将进入集结阶段，速度减小，则相应的转弯半径也减小；在25 s之后，曲线瞬间接近于零，这是由于无人机队形保持定高平飞状态，此时无人机不需要大幅度转弯。

图 6给出了无人机编队从松散到集结过程俯仰角速率变化曲线图。随着时间的变化，该曲线呈现一种先上升后缓慢下降，直至稳态的趋势。在4.8 s(集结)之前，无人机以较短时间且较大俯仰角追踪虚拟长机；在4.8~10 s(集结)，无人机微调姿态进行修正，实现期望的队形；在10 s之后，无人机保持定高平飞状态。在无人机集结整个过程中，爬升(下降)能力处于先增大后缓慢增大，直至最大。

图 7给出了无人机编队规模对稳态飞行影响曲线图。随着时间变化，无人机编队规模N=5, 10, 15, 20的曲线呈现一种先小幅振荡，之后达到稳态的趋势，无人机编队数量大于20的曲线呈现等幅增大后发散的趋势。由于Back-stepping控制方法是一种逐步递推方法，当无人机的数量超过20架，在25 s之后稳态误差逐渐增大，直至无限大，该方法可能会失效。这是由于无人机数量增大后，他们之间的信息交互数据量庞大，易出现延迟导致单个无人机出现指令执行误差，误差的累积使得整个无人机编队系统队形出现故障，甚至发生碰撞引起坠机。因此，本文针对无人机编队队形控制提出的算法使用无人机编队的规模范围为(3＜N≤20)。

5.2 收敛性分析

采用相同初始条件和参数，对无人机编队快速集结和队形重构进行仿真实验，如图 8~图 11所示。

图 8为无人机编队队形保持运动轨迹图。随着时间的变化，无人机编队运动轨迹呈现一种缓慢螺旋上升趋势。4架无人机在不同的初始位置以不同的速度、俯仰角以及偏航角起飞，与此同时，无人机在虚拟长机的引导下以恒定的速度和追踪误差编队飞行。在15 s时，无人机收敛于虚拟长机的飞行轨迹，之后保持相对距离和速度一致稳定飞行。

图 9为无人机编队追踪误差曲线图。由图可知，该曲线呈现一种先缓慢减小，之后缓慢增大直至趋于稳态的趋势。在3 s之前，无人机从不同的起点集结编队飞行，他们与虚拟长机之间的距离逐渐缩小；在3 s之后，无人机之间的相对距离小于安全距离(任意2架无人机之间的相对距离大于翼展长度的2倍)，此时无人机之间的追踪误差缓慢增大，直至大于安全距离；在10 s之后，无人机编队在三维空间的误差均趋于零，按照既定的队形稳定飞行。

图 10为无人机编队相对运动曲线图。随着时间的变化，该曲线呈现一种先迅速增大，之后缓慢减小直至稳态趋势。在集结编队之前，4架无人机的初始速度不同，他们在刚开始2 s相对距离增大；由于虚拟长机的速度最大，其他3架无人机在横向提前达到稳定状态，而横侧向稍有延迟，但最终实现期望的队形。

图 11为无人机编队稳态误差。由图可知，无人机从不同初始位置集结编队，到微调队形，然后保持期望队形飞行。在此过程中，无人机编队稳态误差曲线呈现一种先缓慢增大，然后稍微振荡最后趋于稳态。在1 s之前，曲线呈现缓慢增大，这是由于无人机以不同的速度追踪虚拟长机，无人机之间的相对距离增大；在1~3 s，该曲线呈现波动，这是由于无人机处于从即将达到既定队形到期望队形之间的调整阶段；在3 s之后，无人机编队实现期望的队形，并稳定飞行。

5.3 鲁棒分析

为了进一步验证无人机编队系统在偏航角、俯仰角以及相对距离误差的动态响应和追踪效果，本文以相同的初始条件采用模型预测控制(MPC)方法、拉普拉斯方法以及所提出的方法进行模拟仿真。其仿真结果，如图 12~图 15所示。

图 12为无人机编队偏航角曲线图。随着时间的变化，该曲线呈现一种先快速增大而后缓慢增大，最后趋于稳态的趋势。在无人机编队集结过程中，每架无人机均增大偏航角朝向期望的队形方向追踪虚拟长机；在即将集结既定队形时，无人机调整偏航角实现期望的队形。由图 12可知，采用3种方法均使无人机编队收敛于虚拟长机的飞行轨迹，即达到稳定状态，然而所需要的时间不同。模型预测控制方法所需要时间是最短的，但是无人机偏航角变化较大，不利于地面站人员操控；拉普拉斯方法使得无人机偏航角在刚开始编队较大且所需要的时间最长，会延迟集结所需要的时间；所提方法需要的时间处于其他两种方法之间，且航向角变化较为缓慢，有利于地面站人员操控并缩短集结时间。

图 13为无人机编队俯仰角曲线图。在3种控制方法的作用下，该曲线呈现一种先迅速增大后趋于稳态的趋势。由图易知，3种方法均能使无人机编队达到稳定状态，但所需时间不同，而且拉普拉斯方法使得无人机出现波动现象。模型预测控制方法所需时间最短，但俯仰角编队变化较大；拉普拉斯方法所需时间较短，但是无人机俯仰角变化大；相比较前两种方法，所提方法需要的时间较短且俯仰角变化最小，有利于地面人员操控。

图 14为无人机编队横侧向距离误差曲线图。随着时间的变化，该曲线呈现一种小幅度振荡后趋于稳态的趋势。3种控制方法均能使无人机收敛期望的队形，但需要时间不同。模型预测方法所需的时间较短，无人机之间的横侧向距离误差最大且收敛时间也是最小的，但是集结期望队形所需的时间很长；拉普拉斯方法横侧向距离误差最小，但是所需时间是最长的；所提方法所需时间和相对距离误差均较小，是比较理想的队形控制方法。

图 15为无人机编队俯仰角速率曲线图。由图 15可知，模型预测控制方法、拉普拉斯方法以及所提方法的俯仰角速率变化量分别是10、7和3(°)/s，稳态所需时间分别为1、10和3 s。模型预测控制方法使无人机收敛稳态所需时间最短，但俯仰角速率变化最大；拉普拉斯方法使无人机俯仰角速率变化最小，但收敛所需时间最长；所提方法使无人机俯仰角速率变化较小且收敛时间较短，有利于地面人员操控。

从以上分析可知，通过模型预测控制方法、拉普拉斯方法以及所提方法均能使无人机收敛于期望队形，实现稳定飞行状态。然后，所提方法在动态响应和追踪误差效果均优于模型预测方法和拉普拉斯方法。

6 结论

基于反步推演法，本文提出一种无人机编队队形协同导引控制方法，该方法采用李雅普诺夫稳定性理论进行证明。主要贡献和存在不足为：

1) 通过比较模型预测控制方法、拉普拉斯方法以及所提方法，所提方法不仅使得无人机编队快速集结并形成期望的队形，而且在追踪虚拟长机时，动态响应快和追踪误差小。

2) 本文研究仅仅研究从无人机起飞到构成队形过程的控制，将每架无人机当作质点处理，未考虑每架飞机执行的任务(侦察，打击，评估)、形状以及质量。

3) 本文不是以任务为导向进行建立模型，并对队形进行控制。

4) 无人机编队以期望的队形飞行过程中，未考虑近距编队涡旋效应的影响。

5) 在下一步工作中，将把2)、3)以及4)的因素考虑进去，使多无人机编队以近距高速的队形稳定飞行，并进一步在样机上验证所提出控制方法的有效性。