无人机编队飞行可以提高完成任务的效率，拓宽无人机使用范围，在军用和民用领域均有着广阔的应用前景^[1]。如果编队中的无人机互相碰撞或与空域障碍物发生碰撞，不仅会造成经济损失，而且会严重影响整个战局，因此增强无人机的防碰撞能力是提高编队执行任务精确度和成功率的关键。

针对无人机编队的防碰撞控制问题，国内外学者进行了大量的研究，常见的防碰撞控制方法有人工势场^[2]、最优化理论^[3]、神经网络^[4]等。人工势场法因其实现简单、意义明了、计算量小、实时性好，引起了广泛关注^[5]。人工势场法最初是由Khatib提出的一种虚拟力法，使障碍物的分布情况及其形状等信息反映在环境每一点的势场值当中，根据势场值的大小，无人机决定行进方向和速度^[6]。人工势场实现防碰撞主要依靠斥力势，同时兼顾引力势，以确保人工势场有正确的平衡态^[7]。为解决无人机执行机构的饱和问题，分叉型人工势场被用来评估编队的防碰撞能力^[8]。对于空域障碍物，考虑障碍物闯入的突然性，使用非线性预测控制与人工势场相结合的方式，启动快速规避^[9]。障碍物较多时，使用规避时间、消耗能量等构建防碰撞代价函数，以高效、节能地规避障碍物^[10]。目前，人工势场用于编队防碰撞时存在严重缺陷：无人机编队不能充分利用机间通信，没有考虑通信结构对人工势场的影响，与实际情况不符。另外，大多数文献未清晰阐述如何将人工势场应用到无人机的飞行控制中，缺乏实际应用价值。

近年来，多智能体一致性理论逐渐应用于编队控制、群集运动、姿态同步、状态估计等领域，仅依靠局部信息交互的一致性算法日渐成熟，其智能化程度高、运算代价低，极大地促进了编队控制的发展^[11]。一致性理论深入研究通信机制，给出了通信结构对编队控制影响的普遍规律^[12]。MCLAIN将一致性算法应用于编队控制问题，指出对于任意有向信息流，只要存在有向生成树就能实现精确的编队控制^[13]。文献[14]通过有限时间一致性算法处理了有领导者和无领导者编队的防碰撞问题，以提高编队控制的效率。虽然这些文章只是初步地探索了利用一致性理论解决编队控制问题，未涉及防碰撞控制，也未深入探索通信权重对编队控制的影响，但是仍然催生了新的研究动机—使用一致性理论克服传统人工势场在处理编队防碰撞问题时的固有缺陷。将一致性理论应用于编队控制，引入通信拓扑和通信权重的概念，能够充分利用机间通信，而且依据不同的通信权重可以划分编队中无人机的避碰优先级别，优先保证处于重要地位无人机的安全。

1 机间人工势场

考虑由n架相同型号的无人机组成编队，第i架无人机表示为U_i。令J_i(ρ_i)为无人机U_i与它机U_j之间产生的势场J_ij(ρ_ij)的总和，对于空间中的每一个位置，J(ρ)和J_ij(ρ_ij)都是可微分的。ρ_i为U_i的三维位置，ρ_ij为从U_i到U_j的位置矢量。编队之间的人工势场如图 1所示。

1.1 引入一致性概念改进人工势场

传统人工势场认为无人机可以获得编队中所有它机的信息，这与现实中的通信拓扑往往不相符；而且没有设置通信权重，这就意味着不利于划分编队中无人机的防碰撞优先等级，难以保证重要成员的安全。

为了改进这两个缺陷，将多智能体一致性概念引入人工势场，形成新的人工势场。令所有向U_i通信的无人机构成的集合为N_i，U_i只能获得向其通信的其它无人机的信息，如位置、姿态等。令a_ij表示从U_j到U_i的通信权重，体现了U_j的信息对于U_i的重要程度。无人机U_i的总势与分势的关系定义为

${J_i}\left( {{\mathit{\boldsymbol{\rho }}_i}} \right) = \sum\limits_{j \in {N_i}} {{a_{ij}}{J_{ij}}\left( {\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|} \right)} $

(1)

式中：‖·‖为L₂范数。只有满足j∈N_i的无人机才向U_i通信，这就体现了真实的通信拓扑。同时，a_ij体现了U_j的信息对于U_i的重要程度，当U_i执行避碰动作时，它优先避开连接权重高的无人机。即使编队发生了撞机状况，也是先撞毁处于次要地位的无人机。

J_i(ρ_i)包含引力势J_i^a(ρ_i)和斥力势J_i^r(ρ_i)：

$\begin{array}{*{20}{c}} {{J_i}\left( {{\mathit{\boldsymbol{\rho }}_i}} \right) = J_i^a\left( {{\mathit{\boldsymbol{\rho }}_i}} \right) + J_i^r\left( {{\mathit{\boldsymbol{\rho }}_i}} \right) = }\\ {\sum\limits_{j \in {N_i}} {{a_{ij}}J_{ij}^a\left( {\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|} \right)} + \sum\limits_{j \in {N_i}} {{a_{ij}}J_{ij}^r\left( {\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|} \right)} } \end{array}$

(2)

J_ij^a(‖ρ_ij‖)与J_ij^r(‖ρ_ij‖)为U_i与U_j之间的引力势和斥力势。存在唯一的平衡点‖ρ_ij^d‖使得两机之间的引力势与斥力势能量平衡，此时J_ij(ρ_ij)具有唯一的最小值。因此，通常将ρ_ij^d定义为预定的机间距离矢量。‖ρ_ij‖大于平衡点时引力势起主要作用，小于平衡点时斥力势起主要作用，它们的梯度为

$\left\{ \begin{array}{l} \nabla J_{ij}^a\left( {\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|} \right) = \frac{{\partial J_{ij}^a\left( {\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|} \right)}}{{\partial \left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}\nabla {\mathit{\boldsymbol{\rho }}_{ij}}\\ \nabla J_{ij}^r\left( {\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|} \right) = \frac{{\partial J_{ij}^r\left( {\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|} \right)}}{{\partial \left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}\nabla {\mathit{\boldsymbol{\rho }}_{ij}} \end{array} \right.$

(3)

其中，$\nabla \left( {{\mathit{\boldsymbol{\rho }}_{ij}}} \right) = \frac{{{\mathit{\boldsymbol{\rho }}_i} - {\mathit{\boldsymbol{\rho }}_j}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}$。

1.2 机间人工势场的设计

分别设计机间人工势场的斥力势和引力势，重点设计斥力势。首先，构建机间斥力势为广义Morse函数^[15]：

$J_i^r = \left\{ \begin{array}{l} \sum\limits_{j \in {N_i}} {{a_{ij}}\frac{b}{{{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}_{\min }}}}{c}}}}}} ,\;\;\;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\| \in D\\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\| \notin D \end{array} \right.$

(4)

式中：b、c均为常数，分别决定斥力势的幅值和变化速度，均为可调节参数。D=(‖ρ_ij‖_min, ‖ρ_ij‖_max]决定了机间人工势场的作用区域，‖ρ_ij‖_min＞0为机间最小安全距离，‖ρ_ij‖小于该值则发生机间碰撞；‖ρ_ij‖_max为机间人工势场的最大作用距离。

接下来，为确保人工势场有正确的平衡态‖ρ_ij^d‖，定义机间引力势：

$J_i^a = \left\{ \begin{array}{l} \sum\limits_{j \in {N_i}} {{a_{ij}}\frac{1}{2}{k_{ij}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}^2}} ,\;\;\;\;\;\;\;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\| \in D\\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\| \notin D \end{array} \right.$

(5)

式中：k_ij为正常数，改变k_ij的取值可以调节引力势的强度。

图 2给出了一组J_ij与‖ρ_ij‖的对应关系，取k_ij=0.006 1, b=1, c=1, a_ij=1, ‖ρ_ij‖_min=2, ‖ρ_ij^d‖=4，‖ρ_ij‖_max=10。

从图 2中可以发现，平衡点在‖ρ_ij‖=4处，此时机间斥力势与引力势平衡，机间人工势场为唯一的最小值。2＜‖ρ_ij‖＜4时，机间距离较小，斥力势起主要作用；4＜‖ρ_ij‖≤10时，机间距离较大，引力势起主要作用。‖ρ_ij‖＞10时，机间人工势场不再起作用。

给出机间人工势场后，需要定义一个与距离相关的机间速度场，来实现无人机编队的避碰。综合斥力势(4) 和引力势(5)，对于‖ρ_ij‖∈D，可得到机间速度场函数：

$\begin{array}{*{20}{c}} {{V_i}\left( {{\mathit{\boldsymbol{\rho }}_i}} \right) = - \nabla J_i^a\left( {{\mathit{\boldsymbol{\rho }}_i}} \right) - \nabla J_i^r\left( {{\mathit{\boldsymbol{\rho }}_i}} \right) = }\\ {\sum\limits_{j \in {N_i}} {{a_{ij}}\left[ { - {k_{ij}}\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\| + \frac{b}{c}\frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}_{\min }}}}{c}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}}} \right]\nabla \left( {{\mathit{\boldsymbol{\rho }}_{ij}}} \right)} = }\\ {\sum\limits_{j \in {N_i}} {{a_{ij}}\left[ { - {k_{ij}}\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\| + \frac{b}{c}\frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}_{\min }}}}{c}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}}} \right]\frac{{{\mathit{\boldsymbol{\rho }}_i} - {\mathit{\boldsymbol{\rho }}_j}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}} } \end{array}$

(6)

式(6) 是改进的机间人工势场方法，充分考虑了通信拓扑、通信权重对构建人工势场的影响。只要通信连通，编队就可以避免机间碰撞，实现稳定的群集运动。机间人工势场在平衡点的速度场为0，所以参数k_ij、b、c的取值应满足：

$ - {k_{ij}}\left\| {\mathit{\boldsymbol{\rho }}_{ij}^d} \right\| + \frac{b}{c}\frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {\mathit{\boldsymbol{\rho }}_{ij}^d} \right\|}}{c}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}_{\min }}}}{c}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {\mathit{\boldsymbol{\rho }}_{ij}^d} \right\|}}{c}}} = 0$

(7)

当‖ρ_ij‖≠‖ρ_ij^d‖时，‖V_i(ρ_i)‖≠0。V_i(ρ_i)的大小与方向除了与‖ρ_ij‖有关外，根据式(6)，还取决于无人机U_i与U_j的位置关系ρ_ij。

2 空域障碍物的人工势场

基于人工势场的防碰撞控制还可以扩展到规避空域障碍物，但是与机间避碰机制有所差别。在机间防碰撞过程中，相互靠近的无人机均可做机动调整；而在障碍物规避过程中，只有无人机能主动躲避障碍物，而障碍物不能主动躲避无人机。所以，规避障碍物时，只能单方面调整无人机的轨迹和姿态。另外，不用设计引力势，因为无人机无需靠近障碍物。

假设障碍物的位置为ρ_o=(x_o, y_o, z_o)^T，定义无人机与障碍物间的斥力势为

$\begin{array}{*{20}{c}} {J_i^r\left( {\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|} \right) = }\\ {\left\{ \begin{array}{l} \left[ {1 + k\left( {{V_o}} \right)} \right] \cdot \frac{{{b_o}}}{{{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_o}}}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}_{\min }}}}{{{c_o}}}}}}},\;\;\;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\| \in E\\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\| \notin E \end{array} \right.} \end{array}$

(8)

其中，

$k\left( {{V_o}} \right) = \left\{ \begin{array}{l} {{\rm{e}}^{ - \frac{1}{{{V_o}}}}},{V_o} ＞ 0,\;\;\;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\| \in E\\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right.$

(9)

式中：b_o、c_o均为常数，分别决定斥力势的幅值和变化速度。E=(‖ρ_io‖_min, ‖ρ_io‖_max]决定了与障碍物斥力势的作用区域，‖ρ_io‖_min为无人机与障碍物之间的最小安全距离，‖ρ_io‖_max是避障边界。V_o代表障碍物相对于无人机的运动速度，当障碍物向无人机靠近时，V_o＞0；否则，V_o≤0。无人机周边没有障碍物或者不与障碍物靠近时，k(V_o)=0；当障碍物靠近无人机时，k(V_o)＞0，且靠近速度越高时，k(V_o)值越大。当V_o→0时，k(V_o)→0；当V_o→+∞时，k(V_o)→1。

总之，引入相对速度可以使无人机更具效率地规避障碍物，以补偿障碍物不能主动规避无人机的缺点。

令$\nabla \left( {{\mathit{\boldsymbol{\rho }}_{io}}} \right) = \frac{{{\mathit{\boldsymbol{\rho }}_i} - {\mathit{\boldsymbol{\rho }}_o}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}$，设计‖ρ_io‖∈E时的无人机避障速度场：

$\begin{array}{l} {V_i}\left( {{\mathit{\boldsymbol{\rho }}_{io}}} \right) = - \nabla J_i^r\left( {\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|} \right) = \left( {1 + {{\rm{e}}^{ - \frac{1}{{{V_o}}}}}} \right) \cdot \\ \frac{{{b_o}}}{{{c_o}}} \cdot \frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_o}}}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}_{\min }}}}{{{c_o}}}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_o}}}}}\nabla \left( {{\mathit{\boldsymbol{\rho }}_{io}}} \right) = \\ \left( {1 + {{\rm{e}}^{ - \frac{1}{{{V_o}}}}}} \right) \cdot \frac{{{b_o}}}{{{c_o}}} \cdot \frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_o}}}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}_{\min }}}}{{{c_o}}}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_o}}}}} \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{{\mathit{\boldsymbol{\rho }}_i} - {\mathit{\boldsymbol{\rho }}_o}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}} \end{array}$

(10)

综合机间速度场(6) 和避障速度场(10)，当‖ρ_ij‖∈D且‖ρ_io‖∈E时，U_i的总速度场为

$\begin{array}{*{20}{c}} {v_d^i\left( {{\mathit{\boldsymbol{\rho }}_i}} \right) = {v_i}\left( {{\mathit{\boldsymbol{\rho }}_i}} \right) + {V_i}\left( {{\mathit{\boldsymbol{\rho }}_{io}}} \right) = }\\ {\sum\limits_{j \in {N_i}} {{a_{ij}}\left[ { - {k_{ij}}\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\| + \frac{b}{c}\frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}_{\min }}}}{c}}}} \right)}^2}}}} \right.} }.\\ {\left. {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}}} \right]\frac{{{\mathit{\boldsymbol{\rho }}_i} - {\mathit{\boldsymbol{\rho }}_j}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}} + \left( {1 + {{\rm{e}}^{ - \frac{1}{{{V_o}}}}}} \right) \cdot \frac{{{b_o}}}{{{c_o}}} \cdot }\\ {\frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_0}}}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}_{\min }}}}{{{c_0}}}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_0}}}}}\frac{{{\mathit{\boldsymbol{\rho }}_i} - {\mathit{\boldsymbol{\rho }}_o}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}} \end{array}$

(11)

3 无人机编队防碰撞控制系统设计

得到无人机的速度场后，需根据速度场设计编队防碰撞控制系统。编队防碰撞控制系统，包含防碰撞控制算法和相应的飞行控制器。防碰撞控制算法是外环，产生防碰撞指令；飞行控制器是内环，跟踪防碰撞指令来调整无人机姿态。无人机防碰撞控制系统结构如图 3所示。

3.1 防碰撞控制算法

防碰撞控制算法综合本机、它机和障碍物的信息，经过人工势场的计算，产生速度指令、俯仰角指令、偏航角指令，并将这些防碰撞指令发送给飞行控制器。

防碰撞控制是通过调整各无人机的速度矢量实现的。因此，可将期望的速度矢量定义为

$V_i^d = V_i^d\left( {{\mathit{\boldsymbol{\rho }}_i}} \right) + {V_i}$

(12)

其中，V_i为无人机U_i的速度。根据式(11)，将V_i^d拆分到三个通道(地面坐标系下)，有

$\begin{array}{*{20}{c}} {V_{xi}^d = \sum\limits_{j \in {N_i}} {{a_{ij}}\left[ { - {k_{ij}}\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\| + \frac{b}{c}\frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}_{\min }}}}{c}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}}} \right]} \frac{{{x_i} - {x_j}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}} + \left( {1 + {{\rm{e}}^{ - \frac{1}{{{V_o}}}}}} \right)\frac{{{b_o}}}{{{c_o}}} \cdot }\\ {\frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_o}}}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}_{\min }}}}{{{c_o}}}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_o}}}}}\frac{{{x_i} - {x_o}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}} + {V_{xi}}} \end{array}$

(13)

$\begin{array}{*{20}{c}} {V_{yi}^d = \sum\limits_{j \in {N_i}} {{a_{ij}}\left[ { - {k_{ij}}\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\| + \frac{b}{c}\frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}_{\min }}}}{c}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}}} \right]} \frac{{{y_i} - {y_j}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}} + \left( {1 + {{\rm{e}}^{ - \frac{1}{{{V_o}}}}}} \right)\frac{{{b_o}}}{{{c_o}}} \cdot }\\ {\frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_o}}}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}_{\min }}}}{{{c_o}}}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_o}}}}}\frac{{{y_i} - {y_o}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}} + {V_{yi}}} \end{array}$

(14)

$\begin{array}{*{20}{c}} {V_{zi}^d = \sum\limits_{j \in {N_i}} {{a_{ij}}\left[ { - {k_{ij}}\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\| + \frac{b}{c}\frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}_{\min }}}}{c}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}}{c}}}} \right]} \frac{{{z_i} - {z_j}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{ij}}} \right\|}} + \left( {1 + {{\rm{e}}^{ - \frac{1}{{{V_o}}}}}} \right) \cdot \frac{{{b_o}}}{{{c_o}}}}\\ {\frac{1}{{{{\left( {{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_o}}}}} - {{\rm{e}}^{\frac{{{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}_{\min }}}}{{{c_o}}}}}} \right)}^2}}}{{\rm{e}}^{\frac{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}}{{{c_o}}}}}\frac{{{z_i} - {z_o}}}{{\left\| {{\mathit{\boldsymbol{\rho }}_{io}}} \right\|}} + {V_{zi}}} \end{array}$

(15)

式中：V_xi、V_yi、V_zi为V_i的三轴分量，V_xi^d、V_yi^d、V_zi^d为V_i^d的三轴分量。V_xi^d、V_yi^d、V_zi^d继而转化为一系列的无人机防碰撞指令，包括速度指令V_i^d、俯仰角指令θ_i^d、偏航角指令ψ_i^d，具体形式如下：

$\left\{ \begin{array}{l} V_i^d = \sqrt {V_{xi}^2 + V_{zi}^2 + V_{yi}^2} \\ \theta _i^d = \arctan \left( {\frac{{{V_{zi}}}}{{{V_{xi}}}}} \right)\\ \psi _i^d = \arctan \left( {\frac{{{V_{yi}}}}{{{V_{xi}}}}} \right) \end{array} \right.$

(16)

式(16) 即为基于改进人工势场的防碰撞控制算法。为了方便起见，令v_i=[V_i^d, θ_i^d, ψ_i^d]^T，无人机的飞行控制器跟踪v_i，就可最终实现编队机间防碰撞和对障碍物的规避。

3.2 飞行控制

无人机模型使用六自由度线性模型，无耦合的纵向和横侧向方程分别为

$\begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {{{\dot V}_i}}\\ {{{\dot \alpha }_i}}\\ {{{\dot \theta }_i}}\\ {{{\dot q}_i}}\\ {{{\dot \delta }_{ti}}}\\ {{{\dot \delta }_{ei}}} \end{array}} \right] =\\ \left[ {\begin{array}{*{20}{c}} {{X_V} + {X_{tV}}\cos {\alpha _e}}&{{X_a}}&{ - g\cos {\mu _e}}&0&0&0\\ {\frac{{{Z_V} - {X_{tV}}\sin {\alpha _e}}}{{V - {Z_{\dot \alpha }}}}}&{\frac{{{Z_\alpha }}}{{V - {Z_{\dot \alpha }}}}}&{\frac{{ - g\sin {\mu _e}}}{{V - {Z_{\dot \alpha }}}}}&{\frac{{V + {Z_q}}}{{V - {Z_{\dot \alpha }}}}}&0&0\\ 0&0&0&1&0&0\\ {\frac{{{M_{\dot \alpha }}\left( {{Z_V} - {X_{tV}}\sin {\alpha _e}} \right)}}{{V - {Z_{\dot \alpha }}}} + {M_V} + {M_{tV}}}&{\frac{{{M_{\dot \alpha }}{Z_\alpha }}}{{V - {Z_{\dot \alpha }}}} + {M_\alpha }}&{\frac{{ - g{M_{\dot \alpha }}\sin {\mu _e}}}{{V - {Z_{\dot \alpha }}}}}&{\frac{{{M_{\dot \alpha }}\left( {V + {Z_q}} \right)}}{{V - {Z_{\dot \alpha }}}} + {M_q}}&0&0\\ 0&0&0&0&{ - \frac{1}{{{T_t}}}}&0\\ 0&0&0&0&0&{ - \frac{1}{{{T_e}}}} \end{array}} \right] }\\ {\cdot \left[ {\begin{array}{*{20}{c}} {{V_i}}\\ {{\alpha _i}}\\ {{\theta _i}}\\ {{q_i}}\\ {{\delta _{ti}}}\\ {{\delta _{ei}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{X_{\delta t}}\cos {\alpha _e}}&{{X_{\delta e}}}\\ {\frac{{ - {X_{\delta t}}\sin {\alpha _e}}}{{V - {Z_{\dot \alpha }}}}}&{\frac{{{X_{\delta e}}}}{{V - {Z_{\dot \alpha }}}}}\\ 0&0\\ { - \frac{{{M_{\dot \alpha }}{X_{\delta t}}\sin {\alpha _e}}}{{V - {Z_{\dot \alpha }}}} + {M_{\delta t}}}&{\frac{{{M_{\dot \alpha }}{Z_{\delta e}}}}{{V - {Z_{\dot \alpha }}}} + {M_{\delta e}}}\\ { - \frac{1}{{{T_t}}}}&0\\ 0&{ - \frac{1}{{{T_e}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\delta _{ti}^{\rm{d}}}\\ {\delta _{te}^{\rm{d}}} \end{array}} \right]} \end{array}$

(17)

$\left[ {\begin{array}{*{20}{c}} {{{\dot \beta }_i}}\\ {{{\dot \varphi }_i}}\\ {{{\dot \psi }_i}}\\ {{{\dot p}_i}}\\ {{{\dot r}_i}}\\ {{{\dot \delta }_{ai}}}\\ {{{\dot \delta }_{ri}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\frac{{{Y_B}}}{V}}&{\frac{{g\cos {\mu _e}}}{V}}&0&{\frac{{{Y_p}}}{V}}&{\frac{{{Y_r} - V}}{V}}&0&0\\ 0&0&0&{\frac{{\cos {\mu _e}}}{{\cos {\theta _e}}}}&{\frac{{\sin {\mu _e}}}{{\sin {\theta _e}}}}&0&0\\ 0&0&0&0&1&0&0\\ {L_\beta ^ * }&0&0&{L_p^ * }&{L_r^ * }&0&0\\ {N_\beta ^ * }&0&0&{N_p^ * }&{N_r^ * }&0&0\\ 0&0&0&0&0&{ - \frac{1}{{{T_a}}}}&0\\ 0&0&0&0&0&0&{ - \frac{1}{{{T_r}}}} \end{array}} \right]\\\left[ {\begin{array}{*{20}{c}} {{\beta _i}}\\ {{\varphi _i}}\\ {{\psi _i}}\\ {{p_i}}\\ {{r_i}}\\ {{\delta _{ai}}}\\ {{\delta _{ri}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\frac{{{Y_{\delta a}}}}{V}}&{\frac{{{Y_{\delta r}}}}{V}}\\ 0&0\\ 0&0\\ {L_{\delta a}^ * }&{L_{\delta r}^ * }\\ {N_{\delta a}^ * }&{N_{\delta r}^ * }\\ { - \frac{1}{{{T_a}}}}&0\\ 0&{ - \frac{1}{{{T_r}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\delta _{ai,d}}}\\ {{\delta _{ri,d}}} \end{array}} \right]$

(18)

式中：V_i、α_i、β_i分别为空速、迎角、侧滑角，θ_i、ψ_i、φ_i分别为俯仰角、偏航角、滚转角，q_i、r_i、p_i分别为俯仰角速度、偏航角速度、滚转角速度，δ_ti、δ_ei、δ_ai、δ_ri和δ_{ti, d}、δ_{ei, d}、δ_{ai, d}、δ_{ri, d}分别表示实际和期望的油门舵、升降舵、副翼、方向舵偏转量，T_t、T_e、T_a、T_r为相应舵机的时间常数。g为重力系数，一般取g=9.8 N/kg。气动参数X_V、X_TV、α_e、X_a、Z_a、${Z_{\dot a}}$、M_α、${M_{\dot \alpha }}$、M_V、M_TV、μ_e、Z_q、θ_e、Y_β、Y_p、Y_r、L_β^*、L_p^*、L_r^*、N_β^*、N_p^*、N_r^*的定义详见文献[16]。

将纵向和横侧向方程写为统一形式：

$\left\{ \begin{align} & {{{\mathit{\boldsymbol{\dot{\bar{x}}}}}}_{i}}=\mathit{\boldsymbol{A}}{{{\mathit{\boldsymbol{\bar{x}}}}}_{i}}+\mathit{\boldsymbol{B}}{{\mathit{\boldsymbol{u}}}_{i}}=\mathit{\boldsymbol{A}}\left[ \begin{matrix} \mathit{\boldsymbol{x}}_{i}^{\rm{long}} \\ \mathit{\boldsymbol{x}}_{i}^{\rm{lat}} \\ \end{matrix} \right]+\mathit{\boldsymbol{B}}\left[ \begin{matrix} \mathit{\boldsymbol{u}}_{i}^{\rm{long}} \\ \mathit{\boldsymbol{u}}_{i}^{\rm{lat}} \\ \end{matrix} \right] \\ & {{{\mathit{\boldsymbol{\bar{y}}}}}_{i}}=\mathit{\boldsymbol{C}}{{{\mathit{\boldsymbol{\bar{x}}}}}_{i}}=\mathit{\boldsymbol{C}}\left[ \begin{matrix} \mathit{\boldsymbol{x}}_{i}^{\rm{long}} \\ \mathit{\boldsymbol{x}}_{i}^{\rm{lat}} \\ \end{matrix} \right] \\ \end{align} \right.$

(19)

式中：$\mathit{\boldsymbol{x}}_{i}^{\log }={{\left[ {{V}_{i}},{{\alpha }_{i}},{{\theta }_{i}},{{q}_{i}},{{\delta }_{ei}},{{\delta }_{ti}} \right]}^{\rm{T}}}={{\left[ {{\beta }_{i}},{{\varphi }_{i}},{{\psi }_{i}},{{p}_{i}},{{q}_{i}},{{\delta }_{ai}},{{\delta }_{ri}} \right]}^{\rm{T}}},\mathit{\boldsymbol{u}}_{i}^{\log }={{\left[ {{\delta }_{ei,d}},{{\delta }_{ti,d}} \right]}^{\rm{T}}}$和$\mathit{\boldsymbol{u}}_{i}^{\rm{lat}}={{\left[ {{\delta }_{ai,d}},{{\delta }_{ri,d}} \right]}^{\rm{T}}}$是系统的输入，y_i=[V_i, θ_i, ψ_i]^T是系统的输出。

然后，给出一种PID控制器来实现飞行控制，即图 3中的飞行控制器，跟踪防碰撞指令。PID控制器的结构如图 4所示。

根据图 4设计的结构，令控制误差y_i为

${\mathit{\boldsymbol{e}}_{\bar yi}} = {{\mathit{\boldsymbol{\bar y}}}_i} - {\mathit{\boldsymbol{v}}_i}$

(20)

稳定状态下${{{\mathit{\boldsymbol{\dot{v}}}}}_{i}}=0$，计算式(19)、(20) 的微分：

$\left\{ \begin{align} & \frac{\rm{d}}{\rm{d}t}{{{\mathit{\boldsymbol{\dot{\bar{x}}}}}}_{i}}=\mathit{\boldsymbol{A}}{{{\mathit{\boldsymbol{\dot{\bar{x}}}}}}_{i}}+\mathit{\boldsymbol{B}}{{{\mathit{\boldsymbol{\dot{u}}}}}_{i}} \\ & \frac{\rm{d}}{\rm{d}t}{{\mathit{\boldsymbol{e}}}_{\bar{y}i}}=\mathit{\boldsymbol{C}}{{{\mathit{\boldsymbol{\dot{\bar{x}}}}}}_{i}} \\ \end{align} \right.$

(21)

进而得到

$\frac{\rm{d}}{\rm{d}t}\left[ \begin{matrix} {{{\mathit{\boldsymbol{\dot{\bar{x}}}}}}_{i}} \\ {{\mathit{\boldsymbol{e}}}_{yi}} \\ \end{matrix} \right]=\left[ \begin{matrix} \mathit{\boldsymbol{A}} & 0 \\ \mathit{\boldsymbol{C}} & 0 \\ \end{matrix} \right]\left[ \begin{matrix} {{{\mathit{\boldsymbol{\dot{\bar{x}}}}}}_{i}} \\ {{\mathit{\boldsymbol{e}}}_{\bar{y}i}} \\ \end{matrix} \right]+\left[ \begin{matrix} B \\ 0 \\ \end{matrix} \right]{{{\mathit{\boldsymbol{\dot{u}}}}}_{i}}$

(22)

如果系统(22) 是可控的，需要满足

$\rm{rank}\left[ \begin{matrix} \mathit{\boldsymbol{A}} & 0 \\ \mathit{\boldsymbol{C}} & 0 \\ \end{matrix} \right]={{n}_{A}}+{{n}_{C}}$

(23)

式中：n_A是矩阵A的秩，n_C是矩阵C的秩。

最后，设计PID控制器：

${{\mathit{\boldsymbol{u}}}_{i}}={{\mathit{\boldsymbol{K}}}_{\rm{P}}}{{\mathit{\boldsymbol{e}}}_{\bar{y}i}}+{{\mathit{\boldsymbol{K}}}_{\rm{I}}}\int_{0}^{t}{{{\mathit{\boldsymbol{e}}}_{\bar{y}i}}\rm{d}t}+{{\mathit{\boldsymbol{K}}}_{\rm{D}}}\frac{\rm{d}}{\rm{d}t}{{\mathit{\boldsymbol{e}}}_{\bar{y}i}}$

(24)

式中：K_P、K_I、K_D分别为比例、积分、微分系数矩阵。式(24) 即所设计的飞行控制器，用来跟踪防碰撞指令，与式(16) 共同构成了防碰撞控制系统。

4 仿真实例及分析

仿真采用3架相同的小型固定翼无人机来验证基于改进人工势场的无人机编队防碰撞方法。每架无人机质量85 kg，翼展2.4 m。图 5给出了编队的通信拓扑，从U_j到U_i的边代表U_i可以获得U_j的位置、姿态等信息。

无人机编队为紧密队形，要求编队在水平面内的投影为等边三角形，每条边的长度为4 m；同时，无人机U₂和U₃在无人机U₁下方4 m。给定一条航线的起始和终止坐标位置分别为A(0, 0, 500) m，B(30 000, 30 000, 500) m。采用式(16) 的基于改进人工势场的防碰撞控制算法和飞行控制器(24)，进行仿真，仿真时间为500 s。几个重要参数分别为：‖ρ_ij‖_min=2 m，‖ρ_ij^d‖=4 m；人工势场参数取k_ij=0.006 1, b=1, c=1, b_o=0.1, c_o=0.1。‖ρ_ij‖_min=2 m, ‖ρ_ij‖_max=10 m, ‖ρ_io‖_min=2 m, ‖ρ_io‖_max=1 000 m。

这是一组机间距离较近的初始仿真条件，其中，无人机U₂与U₃之间的距离只有2 m，是机间距离的最小值，需要实施紧急避碰；U₁与U₂之间、U₁与U₃之间也只有2.8 m的间距，同样需要避碰。编队的其他初始仿真参数如表 1所示。

给定不同的通信权重，令a₁₂=2, a₁₃=2, a₂₁=1, a₂₃=1, a₃₁=1, a₃₂=1。U₁在编队中处于最重要位置，防碰撞时优先保证其安全。使用防碰撞控制算法(16) 和飞行控制器(24) 进行仿真。

图 6为速度、俯仰角、偏航角的响应曲线，以及机间平面距离和高度差曲线。为了清楚体现机间防碰撞控制效果，只选取了前5 s。机间距离从不足到回到预定队形，均用了不到2 s的时间，都能满足防碰撞控制的要求。但是，从图 6(a)、(b)的曲线中可以发现，U₂、U₃的速度、俯仰角的变化尺度比U₁更剧烈，优先躲避重要度更高的U₁。图 6(d)、(e)反映了机间距离的变化情况，U₁与U₂、U₁与U₃之间分离地快，U₂与U₃之间分离地慢，也体现了U₁重要度更高。所以，引入不同的通信权重可以有效优先保证高重要度无人机的安全。

图 7为防碰撞控制算法(16) 和飞行控制器(24) 长时间作用下的仿真结果，仿真初始条件与图 6相同，仿真时间为500 s。障碍物坐落在编队的预定航线上，在图上以圆球表示。编队同时对速度、俯仰角、偏航角进行调整，合理地调整航线，实现了对空间障碍物的快速规避。所以，防碰撞控制算法(16) 和飞行控制器(24) 不仅能够实现机间防碰撞，还能够实现对空域障碍物的规避。

5 结论

1) 本文针对无人机编队飞行过程中的机间防碰撞问题，提出了基于改进人工势场的防碰撞控制方法。利用一致性理论中通信拓扑和通信权重的概念，给出了改进的机间人工势场函数。三维飞行仿真结果表明，改进的机间人工势场有效利用了机间通信，不同的通信权重优先保障了编队中处于重要地位无人机的安全。

2) 针对障碍物规避该问题，构建了受无人机与空域障碍物相对运动速度影响的辅助斥力势，提高了无人机编队规避障碍物的效率。仿真结果表明，由于引入了障碍物与无人机的相对运动速度，弥补了障碍物不能主动规避无人机的缺点。

3) 综合本机、它机和障碍物的信息，经过改进人工势场的计算，产生防碰撞指令。并设计了一种PID飞行控制器来跟踪防碰撞指令，进而构建了含防碰撞控制算法、飞行控制器的完整编队系统。该编队系统简单有效，具有良好的工程应用前景。