近几十年来，无人机在军事和民用方面发挥着越来越重要的作用。传统的空战中，主要通过派出单架或几架分散的无人机执行任务。然而，当执行复杂任务时，单架无人机有完成任务效率低、成功率低等缺点。因此，为了充分发挥单架无人机的性能，引入多架无人机以合作方式形成一定形状的编队，以扩大作战范围，提高执行任务的可靠性^[1]。编队飞行是指多架无人机在空间中能够保持一定队形飞行，并且当外界环境或任务发生变化时，编队队形能够根据要求进行变换，常见的编队队形有“人”字型、菱形、三角形等^[2]。不同的编队队形具有各自的功能，例如一般用于侦察的编队可采用“V”形、菱形、锲形或弓形编队；用于拦截的编队则采用行、列或反向圆弧编队形式。由于以上优点，军事上，无人机编队被用于侦察任务，在短时间内能够搜索广大区域^[3]。同样由于此优点，无人机编队被用于森林救火^[4]与农业播种^[5]，在一次飞行中能够覆盖大片区域。在以上任务中，为提高编队中无人机间的协作能力，需要对编队控制方法进行研究。

编队控制作为研究热点，目前已有大量学者在此方向进行了研究，形成了很多编队控制方法，下面对常见的编队控制方法进行讨论。第1种是基于“领导者-跟随者”的方法，编队中飞行状态能够被提前告知的无人机被指定为领导者，而其他无人机作为跟随者，其飞行状态通过跟踪领导者的飞行速度、偏航角、飞行高度以及队形要求进行调整^[6]。“领导者-跟随者”方法具有原理简单、易于操作等优点，但是其鲁棒性较差，跟踪误差会逐级向后传播并被放大。针对以上问题，学者们将“领导者-跟随者”方法与其他方法，例如鲁棒控制^[7]、极值搜索控制^[8]、自适应控制^[9]等相结合，能够减小跟踪误差，但当遇到紧急情况需要重新计算所有无人机的状态时，会增加计算负担。基于行为法是根据从传感器获取的信息来确定无人机的响应，是一种决策方法^[10]。通常无人机有以下4种响应方式，即机间碰撞避免、障碍物规避、目标靠近和队形保持，通过行为响应控制的平均权重来确定每架无人机的行为响应方式^[11]。文献[12]提出了一种基于分布式零空间的行为方法，在速度控制器中设置了多种响应方式的优先级，设计了一种考虑行为优先级的奇异-鲁棒反向运动学方法。基于行为法是通过模拟生物反应行为机制而形成的编队控制方法，具有较强的鲁棒性和灵活性，但是无法实现精确的队形保持，而且对基本行为融合的分析难度较大，这也限制了其发展与应用。虚拟结构法将编队视作一个虚拟的刚体结构，每架无人机都对应着虚拟结构上相对位置固定的一点^[13]。当每架无人机都与固定点保持一定距离时，编队便能形成固定队形。为精确保持编队队形，无人机间需要进行高质量的通信，并且具备较强的计算能力。该方法要求无人机整体保持刚性运动，这使得编队可靠性较差，实用性较低^[14]。

以上3种常用的编队控制方法各有其优缺点，它们都能够被统一到一致性理论的框架下^[15]。一致性理论的基本思想是编队中每架无人机都可以与其邻居无人机进行通信，通过设计一种分布式的算法使得无人机的状态变量达成一致^[16]，以此来实现无人机间的合作。一致性算法的优点是引入了邻居无人机的状态信息作为反馈，只需要邻居无人机间进行信息交换就能形成指定队形，无需复杂的稳定性分析。一致性算法已成为一个研究热点，由于其可靠性高、自愈性强、可扩展性好等优点，已被广泛应用于目标跟踪^[17]、传感器网络设计^[18]、编队控制^[19]等任务中。

当把一致性算法用于求解无人机编队控制问题中时，无人机运动模型、机动性能、飞行性能、飞行安全等因素会从不同方面对控制变量与飞行状态施加限制^[20]，而无人机在飞行过程中必须满足这些约束条件。一致性算法只能保证无人机编队在满足通信拓扑结构的条件下快速收敛到指定飞行状态与队形，但是对于收敛过程中无人机控制量与飞行状态违背约束条件的情况没有进行考虑。另外，虽然目前的研究针对某些飞行状态约束提出了处理策略，但仍不能完全保证无人机的飞行安全。为避免无人机在形成队形过程中发生相互碰撞，文献[21]设计了人工势场法。当2架飞机在水平面内或竖直方向上的距离小于设计阈值时，则在竖直方向上构造人工势场，并将计算出的人工势场力作为修正项添加到用一致性算法生成的控制指令中，在保证无人机飞行安全的前提下使一致性算法仍能收敛。但是，在计算人工势场力时，此方法未包含无人机飞行速度信息，也未对无人机间是否发生碰撞进行提前检测，存在算法启动时间晚、可靠性不高等问题。

针对以上问题，本文基于带自动驾驶仪的无人机三自由度运动学模型，对基本的一致性算法进行改进，使其满足无人机机动性与飞行性能约束条件，并能防止无人机之间的碰撞。本文的创新点如下：

1) 设计一致性算法对无人机编队队形进行控制。将无人机的运动统一到沿OX、OY与OZ轴3个方向进行描述，对编队状态控制律进行改写，并将无人机沿OX、OY与OZ轴3个方向的相对位置信息分别加入到相应方向的编队状态控制律中，形成编队队形控制律。

2) 提出“最小调整”策略对由一致性算法生成的控制指令进行调整，使调整后的控制指令与飞行状态满足无人机机动与飞行性能约束条件，且对控制指令的调整幅度尽量小。考虑到无人机的运动在XOY平面与OZ轴相对独立，调整策略分两步, 分别在XOY平面与沿OZ轴进行。

3) 为避免编队形成过程中无人机间发生碰撞，基于粒子群算法，对沿OZ轴的控制指令再次进行调整，即以改变飞行高度的方法防止无人机之间的碰撞。优化指标为对沿OZ轴的控制指令调整幅度最小。

1 无人机运动模型与一致性算法描述

本节首先建立坐标系描述无人机编队队形，在此基础上给出用于编队控制的无人机运动方程。最后对标准的一致性算法进行介绍，为后续设计编队控制算法做准备。

1.1 无人机编队队形描述

首先建立坐标系描述无人机的空中位置。本文将无人机视为质点，选取地面坐标系描述其运动状态。地面坐标系的原点O可根据无人机的初始位置在水平面上任意选取，OX轴方向也可根据无人机初始状态在水平面上进行选取。OZ轴竖直向上垂直于水平面，OY轴的方向可根据右手法则确定。在编队中，无人机之间的位置关系可通过相对位置关系矩阵R_x、R_y和R_z进行描述

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{R}}_x} = \left[ {\begin{array}{*{20}{c}} {{x_{11}}}&{{x_{12}}}& \cdots &{{x_{1n}}}\\ {{x_{21}}}&{{x_{22}}}& \cdots &{{x_{2n}}}\\ \vdots & \vdots &{}& \vdots \\ {{x_{n1}}}&{{x_{n2}}}& \cdots &{{x_{mn}}} \end{array}} \right]\\ {\mathit{\boldsymbol{R}}_y} = \left[ {\begin{array}{*{20}{c}} {{y_{11}}}&{{y_{12}}}& \cdots &{{y_{1n}}}\\ {{y_{21}}}&{{y_{22}}}& \cdots &{{y_{2n}}}\\ \vdots & \vdots &{}& \vdots \\ {{y_{n1}}}&{{y_{n2}}}& \cdots &{{y_m}} \end{array}} \right]\\ {\mathit{\boldsymbol{R}}_z} = \left[ {\begin{array}{*{20}{c}} {{z_{11}}}&{{z_{12}}}& \cdots &{{z_{1n}}}\\ {{z_{21}}}&{{z_{22}}}& \cdots &{{z_{2n}}}\\ \vdots & \vdots &{}& \vdots \\ {{z_{n1}}}&{{z_{n2}}}& \cdots &{{z_{mn}}} \end{array}} \right] \end{array} \right. $

（1）

式中：(x_ij, y_ij, z_ij)(i, j=1, 2, …, n)为无人机i与j之间的相对位置，并且有x_ii=y_ii=z_ii=0。当编队达到稳定状态时，必须满足

$ \left\{ {\begin{array}{*{20}{l}} {|{x_i} - {x_j}| \to {x_{ij}}}\\ {|{y_i} - {y_j}| \to {y_{ij}}}\\ {|{z_i} - {z_j}| \to {z_{ij}}}\\ {|{V_i} - {V_j}| \to 0} \end{array}} \right. $	（2）

式中：x_i、x_j和z_i分别为编号i的无人机沿OX、OY和OZ轴坐标；V为无人机的飞行速度。

1.2 无人机运动学模型

在无人机编队控制问题中，通常采用带自动驾驶仪的三自由度运动学模型。无人机通过来自速度保持器、俯仰角保持器、航向角保持器的指令信号来保持飞行，即

$ \left\{ {\begin{array}{*{20}{l}} {{{\dot x}_i} = {V_i}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\varphi _i}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _i}}\\ {{{\dot y}_i} = {V_i}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\varphi _i}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _i}}\\ {{{\dot z}_i} = {V_i}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _i}}\\ {{{\dot V}_i} = - \frac{1}{{{\tau _V}}}{V_i} + \frac{1}{{{\tau _V}}}V_i^{\rm{c}}}\\ {{{\dot \varphi }_i} = - \frac{1}{{{\tau _\varphi }}}{\varphi _i} + \frac{1}{{{\tau _\varphi }}}\varphi _i^{\rm{c}}}\\ {{{\ddot \theta }_i} = - \frac{1}{{{\tau _{\dot \theta }}}}{{\dot \theta }_i} - \frac{1}{{{\tau _\theta }}}{\theta _i} + \frac{1}{{{\tau _\theta }}}\theta _i^{\rm{c}}} \end{array}} \right. $

（3）

式中：φ_i为航向角，θ_i为俯仰角；c_i=[V_i^c, φ_i^c, θ_i^c]为指令信号，分别对应于自动驾驶仪的速度指令、航向角指令和俯仰角指令；τ_V、τ_φ、${\tau _{\dot \theta }} $和τ_θ为对应飞行状态的常数。式(3)是无人机编队控制的基本运动学方程，但是其纵向和横航向运动是耦合的。文献[22]将横航向自动驾驶仪和纵向自动驾驶仪进行解耦，得到横纵向分离的运动学模型：

$ \left\{ {\begin{array}{*{20}{l}} {{{\dot x}_i} = {v_i}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\varphi _i}}\\ {{{\dot y}_i} = {v_i}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\varphi _i}}\\ {{{\dot \varphi }_i} = {\omega _i}}\\ {{{\dot v}_i} = \frac{1}{{{\tau _V}}}(v_i^{\rm{c}} - {v_i})}\\ {{{\dot \varphi }_i} = \frac{1}{{{\tau _\varphi }}}(\varphi _i^{\rm{c}} - {\varphi _i})}\\ {{{\ddot z}_i} = - \frac{1}{{{\tau _{\dot z}}}}{{\dot z}_i} + \frac{1}{{{\tau _z}}}(z_i^{\rm{c}} - {z_i})} \end{array}} \right. $

（4）

式中：v_i为无人机i在平面XOY内的速度；ẇ_i为航向角速度；ż_i和$ {\mathit{\boldsymbol{\ddot z}}_i}$分别为无人机的爬升率和爬升加速度; $ {\tau _{\dot z}}$和τ_z分别为爬升率和飞行高度的状态常数; z_i^c为高度指令。在飞行过程中，无人机的速度、加速度、爬升率、爬升加速度都必须在一定范围内变化，以满足机动性与飞行性能的要求。以上约束条件为

$ \left\{ {\begin{array}{*{20}{l}} {{v_i} \in ({v_{{\rm{min}}}},{v_{{\rm{max}}}})}\\ {{{\dot v}_i} \in ({a_{{\rm{min}}}},{a_{{\rm{max}}}})}\\ {{{\dot z}_i} \in ({{\dot z}_{{\rm{min}}}},{{\dot z}_{{\rm{max}}}})}\\ {{{\ddot z}_i} \in ({{\ddot z}_{{\rm{min}}}},{{\ddot z}_{{\rm{max}}}})}\\ {{{\dot \varphi }_i} \in ({\omega _{{\rm{min}}}},{\omega _{{\rm{max}}}})} \end{array}} \right. $	（5）

式中：下标min和max分别表示相应变量的最小值与最大值。

1.3 一致性算法基本原理

对于任意一个航行体，其运动状态可用微分方程描述为

$ {\dot x_i} = v,{\dot v_i} = u\quad i = 1,2, \cdots ,n $	（6）

式中：n为航行体的数量；x_i∈R^m为航行体的状态；v_i∈R^m为航行体状态的微分；u_i∈R^m为控制输入。由式(6)给出基本的一致性算法为

$ {u_i} = - \sum\limits_{j = 1}^n {{a_{ij}}} [({x_i} - {x_j}) + \alpha ({v_i} - {v_j})] $	（7）

式中：i与j表示2个不同的航行体；α为一个表示增益的正数；a_ij为邻接矩阵A_n的元素，表示航行体之间的通信拓扑关系。如果航行体i能接收到航行体j的信息，则a_ij=1，否则a_ij=0。本文中无人机间的通信拓扑关系不是主要的研究内容，假设a_ij=1(i, j=1, 2, …, n, i≠j)，即任意2架无人机间均可以互相向对方发送和接收对方的信息，通信拓扑为一个有向完全图，有关有向完全图的定义和性质可参考文献[23]。

对于n架无人机，定义其通信拓扑结构为G_n，拉普拉斯矩阵L_n中的元素l_ij可定义为

$ {l_{ij}} = \left\{ {\begin{array}{*{20}{l}} { - {a_{ij}}}&{i \ne j}\\ {\sum\limits_{i \ne j,j = 1}^n {{a_{ij}}} }&{i = j} \end{array}} \right. $	（8）

假设无人机在飞行过程中x_i和v_i的通信拓扑关系保持一致，给出一致性算法的收敛条件：

定理  如果G_n中含有一簇有向生成树，并且α>α，则无人机编队的状态能够渐进达到一致。其中，如果-L_n的所有n-1个非零特征值均为负数，则=0，否则

$ \bar \alpha = \mathop {{\rm{max}}}\limits_{\forall {\rm{Re}} ({\mu _i}) < 0, {\rm{Im}} ({\mu _i}) > 0} \sqrt {\frac{2}{{|{\mu _i}|{\rm{cos}}\left( {{\rm{arctan}}\frac{{ {\rm{Im}} ({\mu _i})}}{{ - {\rm{Re}} ({\mu _i})}}} \right)}}} $	（9）

式中:μ_i为-L_n的第i个特征值，Re(μ_i)和Im(μ_i)为μ_i的实数和虚数部分。以上定理的证明过程可参见文献[15]。

以上定理表明，当无人机编队的状态根据一致性算法能够收敛时，即对于任意的初始状态x_i(0)和v_i(0)，当t→∞时，有|x_i(t)-x_j(t)|→0和v_i(t)-v_j(t)→0。要使无人机编队的状态收敛到一致，则要设置合理的通信拓扑结构与α值。

2 改进的一致性算法

鉴于一致性算法式(7)中未考虑对无人机机动与飞行性能的约束(即式(5)中的约束条件)与无人机之间发生碰撞的可能，需要对一致性算法进行改进，以满足无人机编队实际飞行的要求。

首先基于一致性算法对无人机的飞行状态(航向角φ_i、XOY平面内速度v_i、飞行高度z_i)进行控制。由于在飞行状态控制律中未包括编队队形信息(R_x、R_y和R_z)，为了使控制指令与不同方向编队队形信息相对应，设计各运动方向上的一致性算法，使得无人机同时满足状态一致与指定队形要求。

随后对由所设计的一致性算法生成的控制指令进行修正，使其自身和所对应的状态输出均满足无人机机动与飞行性能约束，且无人机之间不发生碰撞。修正分为2个层次进行：先使得无人机控制指令与飞行状态满足机动性与飞行性能约束，再进一步修正OZ轴方向的控制指令防止无人机之间的碰撞。以上修正步骤不能交换的原因在于将满足防碰撞约束的控制指令进一步修正时可能会使得修正后的控制指令再次不满足防碰撞约束，影响无人机飞行安全。

2.1 编队状态控制

状态控制是指让编队中无人机的某些飞行状态分别收敛到相同的值。无人机飞行状态收敛是编队队形收敛的必要条件，故本节基于一致性算法式(7)，首先讨论编队状态控制律的设计问题。

由式(4)可以看出，航向角φ_i、XOY平面内速度v_i、飞行高度z_i是3个独立的运动自由度，当φ_i^c、v_i^c与z_i^c的值(即状态的指令信号值)给定后，无人机的飞行状态便可完全确定。由一致性算法式(7)，可推导出用于求解φ_i^c、v_i^c与z_i^c的控制律为

$ \left\{ \begin{array}{l} \varphi _i^{\rm{c}} = {\varphi _i} + {\tau _\varphi }{u_{\varphi i}}\\ {u_{\varphi i}} = - {b_i}({\varphi _i} - {\varphi ^*}) - \sum\limits_{j = 1}^n {{a_{ij}}} ({\varphi _i} - {\varphi _j}) \end{array} \right. $	（10）

$ \left\{ \begin{array}{l} v_i^{\rm{c}} = {v_i} + {\tau _v}{u_{vi}}\\ {u_{vi}} = - {c_i}({v_i} - {v^*}) - \sum\limits_{j = 1}^n {{a_{ij}}} ({v_i} - {v_j}) \end{array} \right. $	（11）

$ \left\{ {\begin{array}{*{20}{l}} {z_i^{\rm{c}} = {z_i} + \frac{{{\tau _z}}}{{{\tau _{\dot z}}}}{{\dot z}_i} + {\tau _z}{u_{zi}}}\\ {{u_{zi}} = - {d_i}({z_i} - {z^*}) - k{{\dot z}_i} - \sum\limits_{j = 1}^n {{a_{ij}}} [({z_i} - }\\ {{\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} {z_j}) + \gamma ({{\dot z}_i} - {{\dot z}_j})]} \end{array}} \right. $

（12）

式中：φ^*、v^*、z^*分别为无人机编队的理想航向角、XOY平面内理想速度、理想飞行高度，即要求无人机的飞行状态收敛到以上值。τ_v为XOY平面内速度常数; b_i、c_i、d_i是控制系数且要求满足b_i, c_i, d_i>0。可以看出，u_φi、u_vi、u_zi控制律的表达形式与式(7)一致性算法相同，只是多了第1项用于表示对应状态的跟踪目标，而第2项均是保证所有无人机的该状态都能收敛到同一个值。u_φi、u_vi、u_zi可以被视作无人机对当前编队飞行状态的反馈信号，可为生成飞行状态的指令信号(φ_i^c、v_i^c与z_i^c)提供重要信息，见式(10)~式(12)中的第1个等式。

2.2 编队队形控制

式(10)~式(12)的控制律可以使编队中各无人机以相同的飞行速度在指定高度上朝相同方向飞行。但由于上述控制律中未考虑无人机之间的相对位置关系，故无人机无法以指定队形飞行。应当加入无人机位置关系矩阵R_x、R_y和R_z信息，进一步设计编队队形控制律。考虑到R_x和R_y是沿OX与OY轴方向对无人机间的相对位置进行描述的，与式(10)、式(11)中涉及的飞行状态描述方式不同，但相互间又存在转换关系，本节首先将无人机的运动统一在沿OX、OY与OZ轴3个方向进行描述，对式(10)与式(11)的控制律在形式上进行改写。最后将无人机沿OX、OY与OZ轴3个方向的相对位置信息分别加入到相应方向的编队状态控制律中，形成编队队形控制律。

状态控制律中φ_i和v_i与矩阵R_x和R_y中元素物理意义不同，不能简单地直接进行运算，需要根据φ_i和v_i与无人机沿OX和OY轴速度之间的关系进行转换

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{tan}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\varphi _i} = \frac{{{v_{yi}}}}{{{v_{xi}}}}}\\ {{v_i} = \sqrt {v_{xi}^2 + v_{yi}^2} } \end{array}} \right. $	（13）

式中：v_xi与v_yi分别为无人机i沿OX轴与OY轴的速度。通过式(13)，可将φ_i与v_i用v_xi与v_yi来表示，带自动驾驶仪的无人机运动学方程式(4)也可转换为

$ \left\{ {\begin{array}{*{20}{l}} {{{\dot x}_i} = {v_{xi}}}\\ {{{\dot y}_i} = {v_{yi}}}\\ {{{\dot z}_i} = {v_{zi}}}\\ {{{\dot v}_{xi}} = \frac{1}{{{\tau _v}}}(v_{xi}^{\rm{c}} - {v_{xi}})}\\ {{{\dot v}_{yi}} = \frac{1}{{{\tau _v}}}(v_{yi}^{\rm{c}} - {v_{yi}})}\\ {{{\ddot z}_i} = - \frac{1}{{{\tau _{\dot z}}}}{{\dot z}_i} + \frac{1}{{{\tau _z}}}(z_i^{\rm{c}} - {z_i})} \end{array}} \right. $	（14）

式中：v_zi为沿OZ轴的速度，即无人机i的爬升率。式(14)中，无人机控制指令变为c_i=[v_xi^c, v_yi^c, z_i^c]，状态量为x_i=[x_i, y_i, z_i, v_xi, v_yi, v_zi]，而约束条件仍然保持式(5)的形式。

当无人机编队形成稳定队形时，对于每架无人机都有x_i-x_j=x_ij、y_i-y_j=y_ij与z_i-z_j=z_ij，其中x_ij、y_ij与z_ij分别为无人机i与j之间在OX轴、OY轴、OZ轴方向上的理想距离，x_ij、y_ij与z_ij也分别对应于R_x、R_y和R_z中第i行第j列的元素。由基本的一致性算法式(7)，可得到无人机编队控制算法

$ \left\{ \begin{array}{l} v_{xi}^{\rm{c}} = {v_{xi}} + {\tau _v}{u_{xi}}\\ {u_{xi}} = - {b_i}({v_{xi}} - v_x^*) - \mathop \sum \limits_{j = 1}^n [{a_{ij}}({v_{xi}} - {v_{xj}}) - \\ {\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} k_i^x({x_i} - {x_j} - {{\bar x}_{ij}})] \end{array} \right. $

（15）

$ \left\{ \begin{array}{l} v_{yi}^{\rm{c}} = {v_{yi}} + {\tau _v}{u_{yi}}\\ {u_{yi}} = - {c_i}({v_{yi}} - v_y^*) - \mathop \sum \limits_{j = 1}^n [{a_{ij}}({v_{yi}} - {v_{yj}}) - \\ {\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} k_i^y({y_i} - {y_j} - {{\bar y}_{ij}})] \end{array} \right. $

（16）

$ \left\{ \begin{array}{l} z_i^{\rm{c}} = {z_i} + \frac{{{\tau _z}}}{{{\tau _{\dot z}}}}{{\dot z}_i} + {\tau _z}{u_{zi}}\\ {u_{zi}} = - {d_i}(zi - z_{}^*) - k{{\dot z}_i}\mathop \sum \limits_{j = 1}^n {a_{ij}}[({{\dot z}_i} - \\ {\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} {z_j} - {{\bar z}_{ij}}) + \gamma ({{\dot z}_i} - {{\dot z}_j})] \end{array} \right. $

（17）

式中：k_i^x、k_i^y与γ为可调的参数。与状态控制律式(10)~式(12)相比，编队队形控制算法中加入了无人机间的相对位置信息。经式(13)的变换后，无人机在XOY平面内的飞行速度v_i与航向角可分别收敛到$\sqrt {v_x^{*2} + v_y^{*2}} $与$ {\text{arctan}}\frac{{v_y^*}}{{v_x^*}}$，并且能保持指定队形。

需要指出的是，式(15)~式(17)的编队队形控制算法可用于无人机编队集结、队形变换等情形。编队集结是指无人机从不同的初始状态经编队控制算法形成稳定队形的过程。队形变换是根据任务需求，将编队队形由一个稳定状态转换到另一个稳定状态的过程，即以某一稳定队形为初始状态，并改变矩阵R_x、R_y和R_z，利用式(15)~式(17)的编队控制算法便可实现。

2.3 基于“最小调整”策略的约束条件处理方法

2.2节在设计编队队形控制算法时并没有考虑式(5)中的约束条件，使得生成的控制指令及相应的飞行状态可能不满足无人机机动性与飞行性能要求。针对以上问题，本节提出“最小调整”策略对由2.2节编队队形控制算法生成的控制指令进行调整。由于无人机在XOY平面与OZ轴的运动相对独立，故调整分为两步进行。首先对XOY平面内的控制指令u_xi和u_yi进行调整，使其满足XOY平面内加速度${{\dot v}_i} $、速度v_i、航向角速度${{\dot \varphi }_i} $的约束。然后将u_xi和u_yi的值固定，对OZ轴方向的控制指令u_zi进行调整，使其满足OZ轴方向对爬升速度${{\dot z}_i} $与爬升加速度${\ddot z_i}$的约束。

首先对式(15)与式(16)中的u_xi与u_yi进行调整，与这两个变量有关的约束条件为加速度${{\dot v}_i} $、速度v_i以及航向角速度${{\dot \varphi }_i} $。假设当前时刻为t，计算的时间步长为Δt，通过当前无人机飞行状态可计算得到以下信息：

$ \left\{ {\begin{array}{*{20}{l}} {{a_i}(t) = \sqrt {u_{xi}^2(t) + u_{yi}^2(t)} }\\ {{v_i}(t + \Delta t) = {v_i}(t) + {a_i}(t)\Delta t} \end{array}} \right. $	（18）

式中:a_i(t)为t时刻无人机i在XOY平面内的加速度(与${{\dot v}_i} $(t)的含义相同)，v_i(t+Δt)为下一采样时刻无人机i在XOY平面内的速度。如果v_i(t+Δt)不满足约束条件v_i(t+Δt)∈(v_min, v_max)，则定义以下变量：

$ \left\{ {\begin{array}{*{20}{l}} {a_{{\rm{min}},i}^\prime (t) = \frac{{{v_{{\rm{min}}}} - {v_i}(t)}}{{\Delta t}}}\\ {a_{{\rm{max}},i}^\prime (t) = \frac{{{v_{{\rm{max}}}} - {v_i}(t)}}{{\Delta t}}} \end{array}} \right. $	（19）

式中：a′_{min, i}(t)与a′_{max, i}(t)对应于t+Δt时刻飞行速度分别为v_min与v_max时无人机i在t时刻的加速度，即为满足XOY平面内飞行速度v_i(t+Δt)约束，加速度a_i(t)允许的取值范围为[a′_{min, i}(t), a′_{max, i}(t)]。将a′_{min, i}(t)与a_min、a′_{max, i}(t)与a_max分别进行比较，得到更新后的关于加速度a_i(t)的约束条件为

$ \left\{ {\begin{array}{*{20}{l}} {a_{{\rm{min}},i}^{{\rm{ new }}}(t) = {\rm{max}}({a_{{\rm{min}}}},a_{{\rm{min}},i}^\prime (t))}\\ {a_{{\rm{max}},i}^{{\rm{ new }}}(t) = {\rm{min}}(a_{{\rm{max}},i}^\prime (t),{a_{{\rm{max}}}})} \end{array}} \right. $	（20）

式中：a_{min, i}^new(t)与a_{max, i}^new(t)分别为加速度a_i(t)允许的下限与上限。式(20)中实际上也包括了对速度v_i(t+Δt)的约束，因此，只要无人机i的加速度a_i(t)满足式(20)，便可同时满足对XOY平面内加速度与速度的两个约束条件。如果v_i(t+Δt)的值满足约束条件v_i(t+Δt)∈(v_min, v_max)，则a_min与a_max的值无需用式(21)进行更新。下一步将利用约束a_i(t)∈[a_{min, i}^new(t), a_{max, i}^new(t)]对u_xi(t)与u_yi(t)的值进行调整，使它们满足约束条件要求。

从式(15)与式(16)可以看出，u_xi(t)与u_yi(t)即为无人机i在t时刻分别沿OX轴与OY轴的加速度值。由于对加速度a_i(t)的约束条件是相对于XOY平面的，所以对u_xi(t)与u_yi(t)的调整需同步进行，且对原加速度值的影响要尽量小。提出调整u_xi(t)与u_yi(t)的策略为

$ \left\{ \begin{array}{l} u_{xi}^\prime (t) = a_{{\rm{max}},i}^{{\rm{ new }}}(t)\frac{{{u_{xi}}(t)}}{{{a_i}(t)}},u_{yi}^\prime (t) = a_{{\rm{max}},i}^{{\rm{ new }}}(t)\frac{{{u_{yi}}(t)}}{{{a_i}(t)}}\\ {\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} {a_i}(t) > a_{{\rm{max}},i}^{{\rm{ new }}}(t)\\ u_{xi}^\prime (t) = a_{{\rm{min}},i}^{{\rm{ new }}}(t)\frac{{{u_{xi}}(t)}}{{{a_i}(t)}},u_{yi}^\prime (t) = a_{{\rm{min}},i}^{{\rm{ new }}}(t)\frac{{{u_{yi}}(t)}}{{{a_i}(t)}}\\ {\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} {a_i}(t) < a_{{\rm{min}},i}^{{\rm{ new }}}(t) \end{array} \right. $

（21）

式中:u′_xi(t)与u′_yi(t)分别为对u_xi(t)与u_yi(t)进行调整后的值。利用式(21)，可将加速度a_i(t)调整至其约束条件阈值a_{min, i}^new(t)或a_{max, i}^new(t)，并且对u_xi(t)与u_yi(t)进行调整的比例相同(同时放大或缩小相同倍数)，保证了加速度a_i(t)经调整后其方向仍保持不变。经过以上调整，u′_xi(t)与u′_yi(t)的值满足XOY平面内的加速度与速度约束条件，且调整幅度最小(调整至其约束条件阈值)。

下一步将对XOY平面内的航向角速度${{\dot \varphi }_i} $(t)约束进行处理，需要对u′_xi(t)与u′_yi(t)的值进行进一步调整。由于u′_xi(t)与u′_yi(t)的值已满足XOY平面内的加速度与速度约束，因此后续调整必须在继续满足以上约束条件的前提下进行。根据${{\dot \varphi }_i} $(t)∈(ω_min, ω_max)这一约束条件，可得出在下一采样时刻，航向角φ_i(t+Δt)允许的取值范围为

$ \left\{ {\begin{array}{*{20}{l}} {{\varphi _{{\rm{max}},i}}(t + \Delta t) = {\varphi _i}(t) + {\omega _{{\rm{max}}}}\Delta t}\\ {{\varphi _{{\rm{min}},i}}(t + \Delta t) = {\varphi _i}(t) + {\omega _{{\rm{min}}}}\Delta t} \end{array}} \right. $	（22）

同时，下一采样时刻航向角φ_i(t+Δt)的值可由u′_xi(t)与u′_yi(t)的值确定，即

$ {\varphi _i}(t + \Delta t) = {\rm{arctan}}\frac{{{v_{yi}}(t) + u_{yi}^\prime (t)\Delta t}}{{{v_{xi}}(t) + u_{xi}^\prime (t)\Delta t}} $	（23）

式中：v_xi(t)与v_yi(t)分别为t时刻无人机i沿OX轴与OY轴的速度。若φ_i(t+Δt)∉[φ_{min, i}(t+Δt), φ_{max, i}(t+Δt)]，则u′_xi(t)与u′_yi(t)的值可通过求解式(24)或式(25)中的方程组进行调整。

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{v_{yi}}(t) + u_{yi}^{\prime \prime }(t)\Delta t}}{{{v_{xi}}(t) + u_{xi}^{\prime \prime }(t)\Delta t}} = {\rm{tan}}({\varphi _{{\rm{max}},i}}(t + \Delta t))}\\ {u_{xi}^{\prime \prime 2}(t) + u_{yi}^{\prime \prime 2}(t) = a_i^{\prime 2}(t)}\\ {{\varphi _i}(t + \Delta t) > {\varphi _{{\rm{max}},i}}(t + \Delta t)} \end{array}} \right. $	（24）

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{v_{yi}}(t) + u_{yi}^{\prime \prime }(t)\Delta t}}{{{v_{xi}}(t) + u_{xi}^{\prime \prime }(t)\Delta t}} = {\rm{tan}}({\varphi _{{\rm{min}},i}}(t + \Delta t))}\\ {u_{xi}^{\prime \prime 2}(t) + u_{yi}^{\prime \prime 2} = a_i^{\prime 2}(t)}\\ {{\varphi _i}(t + \Delta t) < {\varphi _{{\rm{min}},i}}(t + \Delta t)} \end{array}} \right. $	（25）

式中：u″_xi(t)与u″_yi(t)为经u′_xi(t)与u′_yi(t)调整后的值。式(24)和式(25)为二元二次方程组，通常情况下有两组不同的解，分别记为u″_xi1(t)、u″_yi1(t)与u″_xi2(t)、u″_yi2(t)。由于两组解均满足航向角速度_i的约束条件，需要根据“最小调整”策略进一步确定最终选取哪一组作为调整后的结果。

令$ \eta {\prime _{ai}}\left( t \right) = {\text{arctan}}\frac{{u{\prime _{yi}}\left( t \right)}}{{u{\prime _{xi}}\left( t \right)}}$表征加速度a′_i(t)的方向，调整后的u″_xi(t)与u″_yi(t)值应当不仅保持a′_i(t)的值不变，a′_i(t)的方向也应当变化最小。定义$ {{\eta ''}_{ai1}}\left( t \right) = {\text{arctan}}\frac{{{{u''}_{yi1}}\left( t \right)}}{{{{u''}_{xi1}}\left( t \right)}}$与${{\eta ''}_{ai2}}\left( t \right) = {\text{arctan}}\;\frac{{{{u''}_{yi2}}\left( t \right)}}{{{{u''}_{xi2}}\left( t \right)}} $表征经式(24)或式(25)调整后无人机i在XOY平面内加速度的方向。在两组解中，选取min(|η″_ai1(t)-η′_ai(t)|, |η″_ai2(t)-η′_ai(t)|)所对应的那组解作为最终的u″_xi(t)与u″_yi(t)值。

以上对u_xi(t)与u_yi(t)进行调整的过程中，先将对XOY平面内速度的约束转换到对加速度的约束中，在保证加速度不超出阈值的前提下等比例调整u_xi(t)与u_yi(t)，调整后的u_xi(t)和u_yi(t)值分别为u″_xi(t)和u″_yi(t)。再根据航向角速度约束再次对u″_xi(t)与u″_yi(t)进行调整，最后根据对加速度方向调整最小的原则确定u″_xi(t)与u″_yi(t)的值。以上过程使得对u_xi(t)与u_yi(t)的调整最小且满足XOY平面内的约束条件。

对于沿OZ轴方向的约束条件，即对无人机的爬升率$ {{\dot z}_i}$(t)与爬升加速度${{\ddot z}_i} $(t)进行限制，也按照上述方法对式(17)中的u_zi(t)进行调整。首先通过当前的无人机飞行状态计算下一采样时刻的爬升率$ {{\dot z}_i}$(t+Δt)，即

$ {\dot z_i}(t + \Delta t) = {\dot z_i}(t) + {u_{zi}}(t)\Delta t $	（26）

若$ {{\dot z}_i}$(t+Δt)∉[$ {\dot z_{\min }},{\dot z_{\max }}$]，则更新关于爬升加速度$ {{\ddot z}_i}$(t)的约束条件为

$ \left\{ {\begin{array}{*{20}{l}} {\ddot z_{{\rm{min}},i}^\prime (t) = \frac{{{{\dot z}_{{\rm{min}}}} - {{\dot z}_i}(t)}}{{\Delta t}}}\\ {\ddot z_{{\rm{max}},i}^\prime (t) = \frac{{{{\dot z}_{{\rm{max}}}} - {{\dot z}_i}(t)}}{{\Delta t}}} \end{array}} \right. $	（27）

式中: $ {{\ddot z'}_{{\text{min, }}i}}\left( t \right)$与${{\ddot z'}_{{\text{max, }}i}}\left( t \right) $分别为将t时刻爬升率约束转换为t时刻爬升加速度约束后的约束下限和上限。若${{\dot z}_{\min }} \leqslant {{\dot z}_i}\left( {t + \Delta t} \right) \leqslant {{\dot z}_{\max }} $，则$ {{\ddot z}_{\min }}$与$ {{\ddot z}_{\max }}$的值无需用式(27)进行更新。同样，通过式(28)，也把t时刻对爬升率的约束转换为对爬升加速度的约束。随后比较$ {{\ddot z'}_{{\text{min, }}i}}\left( t \right)$与${{\ddot z}_{\min }} $、$ {{\ddot z'}_{{\text{max, }}i}}\left( t \right)$与$ {{\ddot z}_{\max }}$，确定更新后的爬升加速度约束条件为

$ \left\{ {\begin{array}{*{20}{l}} {\ddot z_{{\rm{min}},i}^{{\rm{new}}}(t) = {\rm{max}}({{\ddot z}_{{\rm{min}}}},\ddot z_{{\rm{min}},i}^\prime (t))}\\ {\ddot z_{{\rm{max}},i}^{{\rm{new}}}(t) = {\rm{min}}(\ddot z_{{\rm{max}}}^\prime ,{{\ddot z}_{{\rm{max}},i}}(t))} \end{array}} \right. $	（28）

式中: $ \ddot z_{{\text{min, }}i}^{{\text{new}}}\left( t \right)$与$ \ddot z_{{\text{max, }}i}^{{\text{new}}}\left( t \right)$为同时包含爬升率与爬升加速度约束条件后爬升加速度最终的下限与上限值。最后将当前的爬升加速度u_zi(t)限制在允许范围内，即

$ u_{zi}^\prime (t) = {\rm{max}}(\ddot z_{\min ,i}^{{\rm{new}}}(t),\min ({u_{zi}}(t),\ddot z_{\max ,i}^{new}(t))) $	（29）

式中:u″_zi(t)为调整后的爬升加速度值。从式(29)可以看出，当u_zi(t)∈ $ \left[ {\ddot z_{{\text{min, }}i}^{{\text{new}}}\left( t \right), \ddot z_{{\text{max, }}i}^{{\text{new}}}\left( t \right)} \right]$时，u″_zi(t)=u_zi(t)；当u_zi(t) < $ {\ddot z_{{\text{min, }}i}^{{\text{new}}}\left( t \right)}$时，u′_zi(t)=$ {\ddot z_{{\text{min, }}i}^{{\text{new}}}\left( t \right)}$；当u_zi(t)> $ {\ddot z_{{\text{max, }}i}^{{\text{new}}}\left( t \right)}$时，u″_zi(t)=$ {\ddot z_{{\text{max, }}i}^{{\text{new}}}\left( t \right)}$。与处理沿OZ轴约束条件相比，处理XOY平面内的约束条件的步骤更复杂。这是因为在XOY平面内约束条件与状态量之间没有直接的微分或积分关系，需要对状态量进行转换才能根据约束条件对其进行调整。对于OZ轴，对爬升率与爬升加速度的约束与需要调整的控制量u_zi(t)直接相关，故调整过程更简单。

2.4 基于粒子群优化的无人机防撞策略

2.3节中设计的对控制指令的“最小调整”策略虽然能保证每架无人机的运动满足机动性和飞行性能的约束条件，但是却未考虑无人机之间发生碰撞的可能。因此，需要对u″_xi(t)、u″_yi(t)与u″_zi(t)进行进一步调整，以满足无人机飞行时的防碰撞要求，使任意2架无人机之间的距离均大于安全距离d_safe。从2.3节中处理XOY平面内约束条件的过程可以看出，OX轴与OY轴运动具有关联性，如果再对u″_xi(t)与u″_yi(t)的值进行调整，则会使得在XOY平面内本来已经满足的约束条件可能被破坏。而无人机沿OZ轴方向的运动相对独立，因此，只需进一步调整u″_zi(t)即可使无人机避免相互碰撞，即通过改变无人机飞行高度的方法来实现防碰撞。设经调整后u″_zi(t)的值变为u″_zi(t)，本节采用离子群算法对u″_zi(t)的值进行优化，优化的目标是使调整最小，即u″_zi(t)-u″_zi(t)最小。约束条件为由无人机机动性与飞行性能确定u″_zi(t)的取值范围，以及无人机之间的距离必须大于安全距离d_safe。

首先根据当前u″_zi(t)值，利用u″_xi(t)与u″_yi(t)值与式(14)计算出t+Δt时刻无人机的位置，判断无人机之间的距离是否小于d_safe，若任意2架无人机之间的距离均大于d_safe，则无需执行防撞策略，u″_zi(t)=u″_zi(t)，否则必须对u″_zi(t)的值进行调整，用粒子群算法优化u″_zi(t)的值。

在此优化问题中，待优化变量为u″_zi(t)，记待优化变量组成的向量为u″_z=[u″_z1(t), u″_z2(t), …u″_zi(t), …, u″_zn(t)](i=1, 2, …, n)，u″_zi(t)的取值在区间$ \left[ {\ddot z_{{\text{min, }}i}^{{\text{new}}}\left( t \right), \ddot z_{{\text{max, }}i}^{{\text{new}}}\left( t \right)} \right]$内。在粒子群算法中的u″_z更新为

$ \begin{array}{l} \Delta {\mathit{\boldsymbol{u}}^{\prime \prime }}_{ze}(o) = \bar \omega \Delta {\mathit{\boldsymbol{u}}^{\prime \prime }}_{ze}(o - 1) + {c_1}{\rm{ rand 1}} \cdot \\ {\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} ({\mathit{\boldsymbol{p}}_{ze}}(o) - {\mathit{\boldsymbol{u}}^{\prime \prime }}_{ze}(o)) + {c_2}{\rm{ rand 2}} \cdot \\ {\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} ({\mathit{\boldsymbol{p}}_{zw}}(o) - {\mathit{\boldsymbol{u}}^{\prime \prime }}_{ze}(o)) \end{array} $

（30）

$ {\mathit{\boldsymbol{u}}^{\prime \prime }}_{ze}(o + 1) = {\mathit{\boldsymbol{u}}^{\prime \prime }}_{ze}(o) + \Delta {\mathit{\boldsymbol{u}}^{\prime \prime }}_{ze}(o) $	（31）

式中:u″_ze(o)(e=1, 2, …, N_p; o=1, 2, …, N_iter)表示在第o次迭代过程中粒子e的优化变量；N_p为粒子的数量，N_iter为最大迭代次数；Δu″_ze(o)为优化变量的增量。另外，在粒子群优化算法中，w为惯性权重; c₁和c₂为加速常数; rand1和rand2为区间[0, 1]内服从均匀分布的随机数; p_ze(o)为粒子e的个体极值对应的优化变量；p_zw(o)为全体极值对应的优化变量。u″_z的初始值在各架无人机允许值的范围内，即区间$ \left[ {\ddot z_{{\text{min, }}i}^{{\text{new}}}\left( t \right), \ddot z_{{\text{max, }}i}^{{\text{new}}}\left( t \right)} \right]$内随机选取。而在标准的粒子群算法中，Δu″_ze的初值在区间[0, 1]内进行随机选取。本问题的优化指标是对u_zi(t)的调整尽量小，可表示为

$ J = \sum\limits_{i = 1}^n | {u^{\prime \prime }}_{zi}(t) - u_{zi}^\prime (t)| $	（32）

此优化问题中，需要满足的约束条件是任意2架无人机间的距离均大于d_safe。如果在迭代过程中某个粒子不满足以上约束条件，则将式(32)中的J值设置成一个很大的数，以表示其对应的解是不可行的，应当被舍弃。

根据以上描述，无人机防撞策略的实质就是用粒子群算法对无人机爬升加速度在其允许变化范围内进行优化，找出使无人机避免碰撞且对爬升加速度调整最小的解。

2.5 改进一致性算法的收敛性证明

2.3与2.4节对由一致性算法生成的控制指令u_xi(t)、u_yi(t)与u_zi(t)进行了调整，以满足飞机机动性、飞行性能与防碰撞的要求。本节通过设计李雅普诺夫函数对改进一致性算法的收敛性进行证明。

本文对一致性算法进行改进的实质是对控制指令在某些时刻进行“限幅”和改变值操作，使调整后控制指令和状态量满足所有约束条件。以上改进过程将由式(6)与式(7)构成的线性系统变成了非线性系统。对于非线性系统的稳定性，需要通过构造李雅普诺夫函数进行证明。在构造李雅普诺夫函数时，通常可通过引入状态量的能量函数，设计李雅普诺夫函数V如式(33)所示，对其求导可得式(34)。

$ V = \frac{1}{2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n n } {({x_i} - {x_j})^2} + \frac{1}{2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{{({v_i} - {v_j})}^2}} } $	（33）

$ \begin{array}{*{20}{l}} {\dot V = \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n n } ({x_i} - {x_j})({v_i} - {v_j}) + }\\ {{\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} \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {({v_i} - {v_j})} } ({u_i} - {u_j})} \end{array} $

（34）

式(34)中包含项u_i-u_j，可利用一致性算法式(7)将其表示成状态变量的形式，推导过程为

$ \begin{array}{l} {u_i} - {u_j} = - \sum\limits_{b = 1}^n {{a_{ik}}} [({x_i} - {x_k}) + \alpha ({v_i} - {v_k})] + \\ {\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} \begin{array}{*{20}{l}} {\sum\limits_{k = 1}^n {{a_{jk}}} [({x_j} - {x_k}) + \alpha ({v_j} - {v_k})] = }\\ {\sum\limits_{k = 1}^n {{a_{jk}}} ({x_j} - {x_k}) - {a_{ik}}({x_i} - {x_k}) + }\\ {\alpha \sum\limits_{k = 1}^n {{a_{jk}}} ({v_j} - {v_k}) - {a_{ik}}({v_i} - {v_k}) = } \end{array}\\ {\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} \begin{array}{*{20}{l}} {{x_j}\sum\limits_{k = 1}^n {{a_{jk}}} - \sum\limits_{k = 1}^n {{a_{jk}}} {x_k} - {x_i}\sum\limits_{k = 1}^n {{a_{ik}}} + \sum\limits_{k = 1}^n {{a_{ik}}} {x_k} + }\\ {\alpha \left( {{v_j}\sum\limits_{k = 1}^n {{a_{jk}}} - \sum\limits_{k = 1}^n {{a_{jk}}} {v_k} - {x_i}\sum\limits_{k = 1}^n {{a_{ik}}} + } \right.} \end{array}\\ {\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} \begin{array}{*{20}{l}} {\left. {\sum\limits_{k = 1}^n {{a_{ik}}} {v_k}} \right) = {x_j}\sum\limits_{k = 1}^n {{a_{jk}}} - {x_i}\sum\limits_{k = 1}^n {{a_{ik}}} - }\\ {\sum\limits_{k = 1}^n {({a_{jk}} - {a_{ik}})} {x_k} + \alpha \left[ {{v_j}\sum\limits_{k = 1}^n {{a_{jk}}} - } \right.}\\ {\left. {{v_i}\sum\limits_{k = 1}^n {{a_{ik}}} - \sum\limits_{k = 1}^n {({a_{jk}} - {a_{ik}})} {v_k}} \right]} \end{array} \end{array} $

（35）

考虑编队通信拓扑结构为有向完全图的情况，则有：

$ {u_i} - {u_j} = n({x_j} - {x_i}) + n\alpha ({v_j} - {v_i}) $	（36）

在改进的一致性算法中，对由式(7)计算出的u_i值进行了修正，其数学形式可表示为u_i-Δu_i。则有：

$ ({u_i} - \Delta {u_i}) - ({u_j} - \Delta {u_j}) = n({x_j} - {x_i}) + n\alpha ({v_j} - {v_i}) $	（37）

将式(37)代入式(36)，有

$ \begin{array}{l} \dot V = \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n n } ({x_i} - {x_j})({v_i} - {v_j}) + \\ {\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} \begin{array}{*{20}{l}} {\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {({v_i} - {v_j})} } [ - n({x_i} - {x_j}) - }\\ {n\alpha ({v_i} - {v_j}) + \Delta {u_i} - \Delta {u_j}] = }\\ {\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {({v_i} - {v_j})} } [ - n\alpha ({v_i} - {v_j}) + } \end{array}\\ {\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} \begin{array}{*{20}{l}} {\Delta {u_i} - \Delta {u_j}] = - n\alpha \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{{({v_i} - {v_j})}^2}} } + }\\ {\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {({v_i} - {v_j})} } (\Delta {u_i} - \Delta {u_j}) \le }\\ { - (n\alpha - 0.5)\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{{({v_i} - {v_j})}^2}} } + }\\ {0.5\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{{(\Delta {u_i} - \Delta {u_j})}^2}} } } \end{array} \end{array} $

（38）

为表达方便，令$ \Delta U = \frac{1}{2}\sum\limits_{i = 1}^n {\sum\limits_{i = 1}^n {{{(\Delta {u_i} - \Delta {u_j})}^2}} } $，则$\dot V $的上界可表示为

$ \dot V \le - \left( {n\alpha - \frac{1}{2}} \right)\sum\limits_{l = 1}^n {\sum\limits_{j = 1}^n {{{({v_i} - {v_j})}^2}} } + \Delta U $	（39）

式(39)中，若$n\alpha - \frac{1}{2} > 0{\text{ }} $，为简化表达，令$ {y_i} = \sum\limits_{j = 1}^n {({v_i} - {v_j})} $，向量y=[y₁, y₂, …, y_n]^T，则式(39)可改写为

$ \dot V \le - \left( {n\alpha - \frac{1}{2}} \right){\mathit{\boldsymbol{y}}^{\rm{T}}}\mathit{\boldsymbol{y}} + \Delta U $	（40）

则当$ \left\| \mathit{\boldsymbol{y}} \right\| \geqslant \sqrt {\frac{{\Delta U}}{{n\alpha - \frac{1}{2}}}} $时，有$\dot V \leqslant 0 $，此时非线性系统稳定且收敛。因此，v_i-v_j最终收敛到的区间为$\left\{ {{v_i} - {v_j} \in {\bf{R}}, {\text{ }}\left\| \mathit{\boldsymbol{y}} \right\| < \sqrt {\frac{{\Delta U}}{{n\alpha - \frac{1}{2}}}} } \right\} $。

以上稳定性证明过程表明，只要ΔU有界，改进的一致性算法就能收敛到一个0附近的区间，即以牺牲一定的稳定误差为代价来保证算法的收敛性。

3 仿真验证与结果分析

为验证所提出改进一致性算法的有效性，本文以8架无人机组成的编队为例，分别设计了3组仿真算例。首先验证状态控制算法的有效性，同时也表明状态控制算法不能收敛到指定队形。第2组算例对基于改进一致性算法的编队控制进行仿真，对比了是否考虑无人机机动性与飞行性能约束条件2种情况，表明了所提出的约束条件处理策略及防碰撞策略的有效性。最后，考虑无人机编队形成后需要变换队形的情况，表明所提出的编队控制算法可用于多种任务。以上算例中，假设编队通信拓扑结构为完全图，即任意2架无人机间均可以向对方发送和接收对方的信息。在第1、2组仿真算例中，无人机队形信息与相对位置关系矩阵R_x、R_y和R_z如图 1与式(41)所示。

$ \begin{array}{l} {\mathit{\boldsymbol{R}}_x} = \left[ {\begin{array}{*{20}{c}} 0&{ - 40}&{ - 20}&{ - 80}&{ - 40}&{ - 100}&{ - 80}&{ - 120}\\ {40}&0&{20}&{ - 40}&0&{ - 60}&{ - 40}&{ - 80}\\ {20}&{ - 20}&0&{ - 60}&{ - 20}&{ - 80}&{ - 60}&{ - 100}\\ {80}&{40}&{60}&0&{40}&{ - 20}&0&{ - 40}\\ {40}&0&{20}&{ - 40}&0&{ - 60}&{ - 40}&{ - 80}\\ {100}&{60}&{80}&{20}&{60}&0&{20}&{ - 20}\\ {80}&{40}&{60}&0&{40}&{ - 20}&0&{ - 40}\\ {120}&{80}&{100}&{40}&{80}&{20}&{40}&0 \end{array}} \right]\\ {\mathit{\boldsymbol{R}}_y} = \left[ {\begin{array}{*{20}{c}} 0&{20}&{40}&{40}&{80}&{80}&{100}&{120}\\ { - 20}&0&{20}&{20}&{60}&{60}&{80}&{100}\\ { - 40}&{ - 20}&0&0&{40}&{40}&{60}&{80}\\ { - 40}&{ - 20}&0&0&{40}&{40}&{60}&{80}\\ { - 80}&{ - 60}&{ - 40}&{ - 40}&0&0&{20}&{40}\\ { - 80}&{ - 60}&{ - 40}&{ - 40}&0&0&{20}&{40}\\ { - 100}&{ - 80}&{ - 60}&{ - 60}&{ - 20}&{ - 20}&0&{20}\\ { - 120}&{ - 100}&{ - 80}&{ - 80}&{ - 40}&{ - 40}&{ - 20}&0 \end{array}} \right]\\ {\mathit{\boldsymbol{R}}_z} = \left[ {\begin{array}{*{20}{l}} 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0 \end{array}} \right] \end{array} $

（41）

由式(41)可以看出，当编队稳定时，所有无人机都处于同一高度飞行。无人机的初始状态如表 1所示，约束条件及算法中的参数设置如表 2所示。

在无人机编队飞行过程中，假设理想飞行状态分别为$ {\varphi ^*} = \frac{{\text{π }}}{4}$、v^*=150 m/s、z^*=100 m，无人机之间的安全距离d_safe=15 m。一致性算法中参数设置为k_i^x=k_i^y=γ=1与b_i=c_i=d_i=1(i=1, 2, …, 8)。粒子群算法中其他参数的设置参考了文献[24]。需要注意的是，在使用离子群算法时，由于本问题中需要优化的变量个数为8(8架无人机沿OZ轴的控制量)，且优化变量的取值范围较小，使得搜索范围有限。经过几次尝试，算法均在迭代15次左右收敛，因此设置最大迭代次数为20。

3.1 无人机编队状态控制

利用2.1节提出的编队状态控制算法，并且采用2.3节的约束条件处理策略，得到仿真结果如图 2与图 3所示。

从图 2中可以看出，无人机编队的速度与航向角可以快速分别趋于一致，并最终等于期望的速度与航向角值。相比于航向角速度，飞行高度的响应速度较慢，但最终还是可以收敛到期望高度。这是因为与速度和航向角的控制律相比，高度控制律是基于二阶微分方程的，即对速度和航向角的控制通过改变加速度与航向角速度来实现，而对高度的控制则通过改变爬升加速度来实现，需要经过两次积分过程才能得到高度值，故响应过程较慢。另外，从图 2(d)XOY平面内飞行轨迹可以看出无人机没有形成任何队形，说明只运用状态控制算法还不能使无人形成队形，需要加入无人机间的相互位置关系，利用编队队形控制算法进行求解。

图 3中检验了所提出的“最小调整”约束条件处理策略的有效性，可以看到，在XOY平面内，加速度、航向角速度均满足约束条件，而且在个别点处加速度与航向角速度接近于临界值，说明了该策略能够起到“限幅”作用。在OZ轴方向，爬升率与爬升加速度也满足约束条件要求。

3.2 无人机编队队形控制

本节采用的参数设置与无人机初始条件均与3.1节相同，在加入无人机间的相对位置信息后，设置无人机间的安全距离d_safe=15 m，采用所提出的编队队形控制算法，得到无人机编队飞行轨迹如图 4所示。

图 4(b)中，通过放大的细节可以看出，采用改进的一致性算法，处于不同初始状态的无人机能够收敛到指定的队形，并保持该队形飞行。图 5与图 6给出了与图 4结果相对应的无人机编队飞行状态和约束条件值。

同样，无人机编队的速度、航向角、飞行高度均可收敛到指定的理想值，且各项约束条件均能满足。图 5中“无人机与其他飞机间最小距离”指的是在某一时刻某1架无人机与其他7架无人机距离中的最小值，最小距离值(最小值为17.77 m)均大于机间最小安全距离，能够防止无人机间相互碰撞。

为进一步说明所提出“最小调整”约束条件处理策略的有效性，给出了在不采用此策略下无人机编队的飞行轨迹与状态，如图 7~图 9所示。

以上仿真结果中，在不考虑约束条件时，无人机编队仍然可以收敛到指定飞行状态，并保持指定队形。但是在达到收敛状态前，在XOY平面内，无人机的加速度与航向角速度在某些时刻均超出了允许范围。对比图 6与图 9中XOY平面内加速度，发现图 6中，在采用本文提出的约束条件调整策略后，a_xy的值并不都处于最小值或最大值处，与图 9中不考虑约束条件时的结果不能相互佐证。原因在于对XOY平面内约束中，不仅包括对加速度和速度的约束，也包括对航向角的约束，而且这些约束间是相互影响的。对于加速度而言，一方面需要进行调整满足加速度与速度约束条件，根据本文提出的调整策略，调整后的加速度值应处于最大值或最小值处；另一方面，为了满足航向角的约束条件，OX轴方向与OY轴方向的加速度需要进行再次调整，调整后加速度就可能出现没有维持在最大值或最小值处的现象。在沿OZ轴方向，无人机爬升率与爬升加速度在某些时刻也违背了约束条件，这使得生成的指令信号对于无人机来说在实际飞行过程中是不可行的。另外，无人机间的最小距离在少数时刻也略小于最小安全距离(最小值为14.05 m)，使得无人机间存在相互碰撞的风险。

下面将用文献[21]中提出的基于人工势场的避障方法与本文提出的防碰撞策略进行比较。在式(17)的基础上，在OZ轴方向加入了人工势场力以防止无人机之间的碰撞，控制律的表达式为

$ \begin{array}{*{20}{l}} {{u_{zi}} = - {d_i}({z_i} - {z^*}) - k{{\dot z}_i} - }\\ {{\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} \sum\limits_{i = 1}^n {{a_{ij}}} [{z_i} - {z_j} - {{\bar z}_{ij}} + \gamma ({{\dot z}_i} - {{\dot z}_j})] + {f_{c{a_i}}}} \end{array} $

（42）

$ {\mathit{\boldsymbol{f}}_{ca}} = - \nabla {U_c} = {\left[ {\frac{{\partial {U_c}}}{{\partial {h_1}}},\frac{{\partial {U_c}}}{{\partial {h_2}}}, \cdots ,\frac{{\partial {U_c}}}{{\partial {h_n}}}} \right]^{\rm{T}}} $	（43）

式中：f_cai为定义的每架无人机的势场力，其具体的推导过程可详见文献[21]中式(7)~式(11)。设置无人机间在XOY平面内的安全距离ΔR=10 m，无人机沿OZ轴方向的安全距离ΔH=10 m。需要注意的是，人工势场力仅当2架无人机间在XOY平面内的距离小于2ΔR或沿OZ轴方向的距离小于2ΔH时才产生，否则仍用式(17)来计算u_zi的值。计算人工势场力时参数K_h的值设为1(见文献[21]中式(8))。在以上设置下，仿真得到无人机与其他飞机间最小距离”随时间变化曲线如图 10所示。

从图 10可以看出，无人机之间的距离最小值为14.88 m，比未考虑机间碰撞的情形(图 8)相比无人机间的距离有所增加，但仍小于本文所规定的安全距离(d_safe=15 m)。以上结果说明人工势场法在某些情况下并不能完全使无人机避免碰撞，原因在于未考虑无人机飞行速度的影响，可靠性不高，有待进一步改进。

通过以上仿真结果的对比可以看出，本文提出的“最小调整”约束条件处理策略与基于粒子群算法的防碰撞策略有效，能够保证生成的指令信号在无人机实际飞行过程中有效，并能防止机间碰撞。

3.3 无人机编队队形变换飞行控制仿真

为说明本文设计的一致性算法可用于无人机编队队形变换任务，本节用一个综合算例进行仿真。仿真中无人机编队的飞行过程设置如下：

无人机的初始状态与3.2节中设置相同，理想编队队形描述如图 11与式(44)所示。当t=30 s时，要求编队变换到图 1与式(41)描述的队形，队形的理想速度、航向角、飞行高度仍保持不变。以上情形如图 11所示。

在以上设置下，仿真结果如图 12~图 14所示。

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{R}}_x} = \left[ {\begin{array}{*{20}{c}} 0&{ - 20}&0&{ - 40}&{ - 20}&{ - 40}&{ - 20}&{ - 40}\\ {20}&0&{20}&{ - 20}&0&{ - 20}&0&{ - 20}\\ 0&{ - 20}&0&{ - 40}&{ - 20}&{ - 40}&{ - 20}&{ - 40}\\ {40}&{20}&{40}&0&{20}&0&{20}&0\\ {20}&0&{20}&{ - 20}&0&{ - 20}&0&{ - 20}\\ {40}&{20}&{40}&0&{20}&0&{20}&0\\ {20}&0&{20}&{ - 20}&0&{ - 20}&0&{ - 20}\\ {40}&{20}&{40}&0&{20}&0&{20}&0 \end{array}} \right]\\ {\mathit{\boldsymbol{R}}_y} = \left[ {\begin{array}{*{20}{c}} 0&0&{20}&0&{20}&{20}&{40}&{40}\\ 0&0&{20}&0&{20}&{20}&{40}&{40}\\ { - 20}&{ - 20}&0&{ - 20}&0&0&{20}&{20}\\ 0&0&{20}&0&{20}&{20}&{40}&{40}\\ { - 20}&{ - 20}&0&{ - 20}&0&0&{20}&{20}\\ { - 20}&{ - 20}&0&{ - 20}&0&0&{20}&{20}\\ { - 40}&{ - 40}&{ - 20}&{ - 40}&{ - 20}&{ - 20}&0&0\\ { - 40}&{ - 40}&{ - 20}&{ - 40}&{ - 20}&{ - 20}&0&0 \end{array}} \right],{\mathit{\boldsymbol{R}}_z} = \left[ {\begin{array}{*{20}{l}} 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0 \end{array}} \right] \end{array} \right. $

（44）

图 12中，当t=28 s时，无人机编队队形为五边形，到t=40 s时，无人机编队队形已变换为指定的菱形形状。在图 13中，当t=30 s发出编队队形变换指令时，无人机速度与航向角会立即发生变化，随后会很快收敛到理想状态。队形变换过程中，无人机飞行高度不发生变化。由图 14中的结果可知，在整个飞行过程中，无人机的控制信号、飞行状态与机间距离均满足约束条件，保证了所生成的指令信号在实际过程中可行。以上仿真结果表明，本文提出的基于改进一致性算法的编队控制算法可以实现队形形成与变换功能。

4 结论

1) 利用无人机各飞行状态间的几何关系，基于带自动驾驶仪的三自由度运动学方程，对基本一致性算法进行改进，在状态控制算法的基础上设计了编队队形控制算法。

2) 提出“最小调整”策略对由编队控制算法生成的指令信号进行调整，使得指令信号满足无人机机动性与飞行性能约束条件，在实际飞行过程中可用。对指令信号的调整分别在XOY平面与OZ轴两个方向进行。

3) 设计基于粒子群算法的无人机防碰撞策略。无人机之间的防碰撞措施在OZ轴方向进行，在允许范围内优化各无人机爬升加速度值，避免碰撞，优化指标为对各无人机原爬升加速度调整之和最小。

4) 仿真结果表明，所提出的改进一致性无人机编队控制算法可用于无人机编队队形成形、变换等任务，具有良好的通用性。