未知环境下无人机集群协同区域搜索算法

2019-03-05侯岳奇梁晓龙何吕龙刘流

北京航空航天大学学报 2019年2期

侯岳奇，梁晓龙,*，何吕龙，刘流

(1. 空军工程大学国家空管防相撞技术重点实验室, 西安 710051；2. 空军工程大学陕西省电子信息系统综合集成重点实验室, 西安 710051)

随着高精度影像设备与技术的快速发展，携带照相设备的无人机(Unmanned Aerial Vehicle，UAV)在民用和军事领域得到了广泛应用，例如环境监测、战场监视以及目标搜索等[1-2]。无人机集群通过无人机之间的协同合作，从而实现整体能力上的涌现，即系统涌现出的能力远超系统内单架无人机能力的总和[3-4]。相比于单架无人机，利用多架无人机协同执行区域搜索任务得到了越来越广泛的关注[5-6]。

针对协同搜索问题，诸多学者进行了深入的探索并取得丰硕的成果。离线规划方法[7-9]将任务区域进行分割，设计各个子区域的覆盖搜索航线，飞行航线固定。该方法的优势在于能够实现任务区域的全覆盖，而在无人机故障、火力威胁等突发情况下该方法受限。文献[5]使用笛卡儿栅格描述环境，赋予每个栅格一个值代表目标分布的不确定性，设计了搜索回报函数和禁飞区回避策略。在有先验信息的情况下，该方法可以实现重点侦察和覆盖搜索。文献[10]在文献[5]的基础上，考虑通信约束，分析了不同通信约束对协同搜索效率的影响。在通信距离和角度有约束的情况下，该方法能够较好地实现协同区域搜索。此类方法依赖先验信息来建立概率地图，然而在实际应用中，先验信息的获取和栅格量化是十分困难的。文献[11]设计了基于信息素的网格回访机制，引导无人机对目标存在概率较大的区域进行回访搜索，使得无人机能够尽早搜索到更多目标。上述方法[5,10-11]均假设无人机在相邻栅格之间运动，这种“粗粒度”的运动模型虽然简化了协同搜索决策的解空间，但却在一定程度上降低了决策结果的精细程度。文献[12-14]基于分布式模型预测控制框架，采用纳什最优和粒子群优化相结合的算法，有效地降低了协同搜索决策问题的求解规模和通信负担。虽然上述研究在一定程度上使多无人机具备了协同搜索的能力，但仍存在区域覆盖率较低的问题。存在这一问题的原因在于：缺少专门的引导机制，来引导无人机向未覆盖区域进行搜索。

在执行区域搜索任务时，无人机故障和环境中火力威胁等突发情况的影响是不容忽视的。考虑到传统的离线规划方法难以应对突发情况，本文提出了一种以覆盖率为实时搜索奖励的协同区域搜索算法，以改善一般在线规划搜索方法覆盖率不高的问题。

1 协同搜索问题建模

1.1 问题描述

多无人机协同搜索问题主要分为2类[15]：第1类是在环境信息已知的条件下，根据先验信息对区域进行搜索，以尽快发现目标，或者对重点区域进行监控；第2类是在环境信息未知的条件下，对区域进行覆盖。本文针对第2类问题展开研究。UAV集群携带通信设备和照相设备对特定任务区域展开搜索，区域内目标的位置分布未知，如图1所示。

UAV按照实时规划的航路对任务区域展开搜索，并通过通信组网模块进行实时通信，为实时决策提供信息支撑，通信内容包括UAV状态、环境信息和决策信息等。当某架UAV出现故障时，传统的离线规划方法难以适应性地完成搜索任务。相比而言，在线规划方法鲁棒性更强，在出现突发情况时，无人机集群能够继续保持协同，自组织地完成搜索任务。因此，本文主要研究如何建立一种高效的区域覆盖在线规划方法，确保UAV集群在尽可能短的时间内对区域进行覆盖。

图1 UAV集群协同搜索示意图Fig.1 Schematic diagram of UAV swarm cooperative search

1.2 环境模型

设任务区域Ω为Lx×Ly的矩形区域，将区域按照固定间隔Δd栅格化为M×N个栅格，如图2所示。

赋予每个栅格一个值μij(k)，用于描述截止到k时刻为止栅格(i,j)是否已被覆盖。为简化分析，做出如下假设：一旦栅格(i,j)处于UAV传感器探测范围内，则认为该栅格已被覆盖，且栅格内的目标存在情况完全已知。栅格的状态μij(k)表示为

(1)

式中：Ωc(k)为已覆盖的栅格集合；Ωnc(k)为未覆盖的栅格集合。任务区域Ω为Ωc(k)和Ωnc(k)的总和。构建覆盖分布地图(Coverage Distribution Map,CDM)来描述任务区域的覆盖分布情况，用矩阵形式来表示覆盖分布地图，定义环境矩阵为

C(k)=[μij(k)]M×N

(2)

在搜索开始前，由于整个区域的环境信息完全未知，且需要对整个区域展开覆盖搜索，因此每个栅格的值μij(0)=1；随着搜索的进行，μij(k)实时变化，覆盖分布地图实时更新，并在UAV集群内共享，为UAV的实时决策提供环境信息。

图2 任务区域栅格化Fig.2 Mission area rasterization

1.3 UAV模型

为简化分析，假设UAV在任务区域上空等高度飞行，并将UAV视为二维空间中运动的质点，其运动方程为

(3)

式中：[xi(k),yi(k)]为UAV的位置；φi(k)=φi(k-1)+Δφi(k)，φi为航向角，Δφi∈[-φmax,φmax]为航向偏转角，为决策输入，φmax为受机动性能限制下的最大转弯角；v0为平飞速度；Δt为时间步长。

将UAV集群视为一个控制系统，并将每架UAV视为一个子系统。k时刻，记子系统UAVi的状态为pi(k)=(xi(k),yi(k),φi(k))。状态方程为

pi(k+1)=f(pi(k),ui(k))

(4)

式中：ui(k)=Δφi(k)为子系统UAVi的输入；f(·)为状态转移函数，由式(3)确定。

根据式(4)，建立子系统UAVi的预测模型为

i=1,2,…,n;j=1,2,…,H

(5)

式中：n为UAV架数；H为预测周期；pi(k+j|k)为基于k+j-1时刻UAVi状态预测的k+j时刻的UAVi状态，其值取决于状态pi(k+j-1|k)和控制输入ui(k+j-1|k)。

从k时刻起，给定H步预测控制输入后，可以预测出未来H步以内UAV的航路，如图3所示。通过优化预测控制输入，来引导UAV尽可能向尚未被覆盖的区域进行搜索，以获得更高的区域覆盖率。

图3 H步预测航路图Fig.3 H step route prediction

2 覆盖分布地图更新方法

覆盖分布地图描述了任务区域的覆盖搜索情况，是UAV进行自主决策的重要依据。因此，如何快速更新覆盖分布地图，对在线实时决策具有重要意义。

覆盖分布地图的更新就是将传感器探测范围内覆盖部分的相应栅格置为0。通过逐一判断邻近栅格到UAV之间的距离是否小于传感器探测距离，并对覆盖范围内的栅格进行逐一赋值，可以实现覆盖分布地图的更新，但这种遍历的方法运算量大、算法复杂度高，不利于实时更新。本文利用Hadamard积进行覆盖分布地图更新，操作简单且易于实现，避免了遍历判断和逐一赋值，其流程如图4所示。

图4中探测矩阵和环境子矩阵的定义，以及对其进行Hadamard积的运算过程将在2.1节和 2.2节中详细介绍。

图4 覆盖分布地图更新流程Fig.4 Update process of coverage distribution map

2.1 探测矩阵

本文将传感器探测范围简化为：UAV所处位置为中心，半径为Rs的圆形区域。如图5所示，正方形区域记为Ψ，将区域Ψ栅格化为Q×Q个栅格。

(6)

式中：Rs为传感器探测半径；「⎤为向上取整函数。

赋予每个栅格一个值ηpq，用于描述UAV传感器能否探测到该栅格。若栅格(p,q)处于UAV传感器探测范围内，则令ηpq=0，反之，ηpq=1。ηpq表示为

(7)

式中：Ψc为传感器可探测区域；Ψnc为不可探测区域。区域Ψ为Ψc和Ψnc的总和。

图5 区域Ψ栅格化Fig.5 Region Ψ rasterization

以ηpq为元素建立探测矩阵D=[ηpq]Q×Q，表示UAV对邻近栅格的覆盖能力。

2.2 覆盖分布地图的更新

定义环境子矩阵C(mn)(k)：在环境矩阵C(k)中，以μmn(k)为中心元素，维度为Q×Q的子块矩阵，称为环境子矩阵。

(8)

(9)

当UAV处于栅格(m,n)内时，可近似认为UAV处于该栅格的中心。此时，环境子矩阵C(mn)(k)与探测矩阵D重合，且维度相等。对上述2个矩阵做Hadamard积

C(mn)(k+1)=C(mn)(k)∘D=

(10)

式中:“∘”为Hadamard积运算；p,q=1,2,…,Q。

假设传感器的探测周期为Ts，以Ts为时间间隔将UAV从k到k+1时刻的直线运动离散为运动点迹。对每个离散点做上述运算，即可实现一个步长内覆盖分布地图的更新。

2.3 覆盖分布地图的信息融合

文献[16]提出了一种地图信息融合更新方法，该方法适用于通信理想的情况，在通信中断或数据丢包的情况下，部分地图信息无法融合更新，对搜索过程造成影响。因此，本文提出基于Hadamard积的地图信息融合方法(见图6)，能够在一定程度上减小通信中断或数据丢包对搜索过程的影响。

在协同搜索过程中，各UAV在本机进行覆盖分布地图更新，更新完毕后通过通信网络广播发送实现融合共享。k时刻，UAVi本地的覆盖分布地图对应的环境矩阵记为Ci(k)，并接收到其他UAV的广播消息记为Cj≠i(k)。基于Hadamard积可实现各UAV覆盖分布地图的信息融合：

(11)

与文献[16]方法相比，基于Hadamard积的地图信息融合方法共享的探测信息为环境矩阵Ci(k)，是地图的全局信息。假如k时刻数据丢包或通信中断，当前时刻的探测信息会丢失，地图信息无法更新。一旦k+x时刻通信恢复正常，环境矩阵Ci(k+x)中已包含了k到k+x时刻的历史探测信息，经过基于Hadamard积的地图信息融合方法更新后，丢失的历史探测信息就得以恢复。虽然k到k+x时刻通信中断会影响当前搜索决策，但是通信恢复后所有历史信息得到恢复，k+x时刻之后的搜索过程不会受到影响。

图6 基于Hadamard积的地图信息融合示意图Fig.6 Schematic diagram of map information fusion based on Hadamard product

3 奖励函数与搜索算法

3.1 搜索奖励函数

协同搜索的关键在于设计一个搜索奖励函数，对每条预测航路进行评估[5]。奖励的设定主要是一个步长内覆盖率增量的大小，并依据边界条件和转弯角度设计惩罚函数。在搜索过程中，UAV基于当前状态和覆盖分布地图，利用搜索奖励函数对预测航路进行评估，并自主地选择奖励值最大的航路作为决策输入。每架UAV使用搜索奖励函数J选择它的搜索航路：

J(pi(k),ui(k))=ω1γJC(k)+

ω2JT(k)+ω3JB(k)

(12)

式中：JC为覆盖率增量；JT和JB分别为转弯角度和边界距离的惩罚函数；ω1、ω2和ω3为相应权重；γ为重要性因子，覆盖率的重要性通过调整γ来体现，γ≥1。

k时刻的区域覆盖率O(k)是已搜索区域Ωc(k)占任务区域Ω的面积比，即已搜索栅格数量与总栅格数量的比值：

(13)

k到k+1时刻的覆盖率增量是指O(k+1)与O(k)的差值。实际意义为：k时刻，未搜索区域Ωnc(k)中，在k到k+1时刻内被搜索到的区域占任务区域Ω的面积比。覆盖率奖励函数为

JC(k)=O(k+1)-O(k)=

(14)

在执行任务过程中，若转弯角度过大，导致耗油增大，影响续航时间。因此，设计一个惩罚函数，尽可能减少UAV转弯角度过大引起的耗油代价。转弯角度的惩罚函数JT可以表示为

(15)

在搜索过程中，距离边界越近，传感器覆盖的有效区域越少，效率越低。借鉴虚拟势函数的思想，设计一个惩罚函数，靠近边界的UAV会受到边界的虚拟“斥力”，距离边界越近则“斥力”越大。因此，边界距离的惩罚函数JB可以表示为

(16)

(17)

3.2 基于DMPC和DE的协同搜索算法

模型预测控制(Model Predictive Control，MPC)是一种利用控制系统模型和优化技术设计预测周期内系统最优控制输入的方法，核心思想是滚动优化求解[17]。集中式MPC方法依赖中央节点进行决策，限制了系统规模的扩展和决策速度，在实际应用中有一定局限性[14,18]。考虑到UAV子系统之间不存在耦合性，即不同UAV的控制是相对独立的，它们在系统的动态特性上并没有关联。为提高整个系统的抗毁性和决策速度，其控制结构可以采用分布式模型预测控制(Distributed Model Predictive Control，DMPC)方式[19]，如图7所示。

MPCi决策流程分为3步。

步骤1预测

图7 DMPC框架示意图Fig.7 Schematic diagram of DMPC framework

图8 MPCi决策流程Fig.8 Decision-making process of MPCi

在预测阶段，每架UAV基于本地覆盖地图和自身状态进行优化求解，而不考虑其他UAV的运动，即UAV之间不进行协同。根据式(17)，求解UAVi的H步预测控制输入的优化模型可以描述为

(18)

步骤2通信

步骤3决策

(19)

搜索决策过程的算法伪代码如下：

1 初始化任务参数

2 fork= 1 tokmax

9 将Ci(k)更新为Ci(k+1)

10 ifO(k+1) ≥ 设定阈值

11 break

12 end if

13 end for

控制输入Ui(k)中包含H个未知变量，此优化问题是一个非线性优化问题。考虑到DE算法在求解优化问题，尤其是非线性优化问题中的优势[20]，采用DE算法进行子系统本地优化求解，算法细节不再赘述。

4 仿真分析

为验证本文方法的有效性，本节对其进行仿真验证。仿真环境为I7-4960，主频2.60 GHz，16 GB内存，基于MATLAB 2014a为平台进行仿真实验。

4.1 实验 1

设任务区域为8 km×8 km的矩形区域，每个栅格大小为20 m×20 m。执行搜索任务的4架UAV从不同位置进入搜索区域，其进入点坐标分别为(200,0) m，(2 200,0) m，(4 200,0) m，(6 200,0) m。UAV之间的通信均为理想条件。UAV平飞速度v0=30 m/s，传感器探测半径Rs=200 m，最大转弯角φmax=60°，仿真步长Δt=10 s，预测步长H=3，ω1=0.9，ω2=0.05，ω3=0.05，γ=200。设定仿真环境条件：任务环境为无遮挡的平地区域，任务区域中存在分布未知的火力威胁，有一定几率造成UAV设备故障。在仿真中加入突发情况来模拟这一环境条件：运行至20步长时，假定UAV1设备故障，停止执行任务；运行至100步长时，假定UAV3设备故障，停止执行任务。为体现本文算法的优势，分别运用平行搜索、随机搜索和本文算法进行对比仿真。仿真结果如图9所示(图中黑色区域为UAV传感器未覆盖的区域)。

图9 3种搜索方法的仿真结果Fig.9 Simulation results of three search methods

如图9(a)所示，在2架友机相继发生故障的情况下，UAV2和UAV4只完成了各自预先分配的搜索任务，故障UAV未完成的搜索任务得不到继续执行。如图9(b)所示，随机搜索方法作为一种无引导机制的在线规划方法，在任务过程中多处存在UAV航迹交叉重叠的情况，搜索效率较低。如图9(c)所示，在搜索初始阶段，各UAV之间保持传感器探测范围尽可能不重叠，实现较高的覆盖率增长。在2架友机相继故障的情况下，UAV2和UAV4通过实时协同，保持原有的搜索策略，充分发挥各自的搜索能力，继续完成搜索任务。总的来看，本文算法各UAV探测区域之间重叠部分较少，在出现突发情况时，能够继续完成搜索任务，体现了无人机集群在线协同的优势。

为消除随机因素的影响，在相同仿真条件下，针对本文算法和随机搜索方法使用蒙特卡罗方法进行500次仿真，得到平均覆盖率随时间变化曲线如图10所示。

图10 3种搜索方法的覆盖率变化曲线Fig.10 Coverage rate changing curve of three search methods

如图10所示，平行搜索方法和本文算法覆盖率高于随机搜索方法。平行搜索方法的覆盖率变化曲线呈折线状：当某架UAV发生故障时，覆盖率曲线的斜率随之降低；当UAV到达边界时，执行转弯程序，覆盖率斜率为零。在搜索初期，本文算法与平行搜索方法覆盖率曲线斜率基本一致，体现了较高的搜索效率。当UAV2和UAV4完成各自的搜索任务后，平行搜索方法的覆盖率保持不变。相比之下，本文算法在UAV1和UAV3发生故障后，仍能保持覆盖率的稳定增长，最终任务结束时的覆盖率远高于平行搜索方法。

4.2 实验 2

针对本文算法在集群规模较大时的有效性进行验证，采用10架无人机组成的无人机集群进行仿真。设任务区域为15 km×15 km的矩形区域，执行搜索任务的10架UAV的初始位置和航向由程序随机产生，UAV故障随机指定：设定UAV2、UAV3、UAV8、UAV9和UAV10分别于仿真步数为185、120、175、170和165时发生故障，其他仿真条件同实验1。用本文算法进行仿真，仿真步数step为100和230时的仿真结果如图11所示。

由图11(a)可以看出，在搜索初期，10架UAV保持传感器探测范围尽可能不重叠，以获得较高的覆盖率增长。由图11(b)可以看出，在搜索后期，由于未搜索区域被分割成多个不规则的形状，UAV航路出现了部分交叉，但本文算法还是能够引导UAV尽可能向未搜索区域移动。搜索过程中覆盖率变化曲线如图12所示。

如图12所示，在搜索前期，覆盖率增长速度较快，覆盖率随时间变化曲线的斜率保持稳定。到搜索后期，由于出现重复搜索的情况，覆盖率增长放缓，斜率逐渐降低。当搜索时间到达38.33 min时，覆盖率达到了90.13%，有效地对区域进行了覆盖。仿真结果表明，在集群规模达到10 架的情况下，本文算法能够较好地完成搜索任务。