0 引　言

载人潜水器作为在海洋勘察与开发中不可或缺的水下作业平台，在深海资源勘探、海洋科学研究等工作中发挥着重要作用。我国自主研制的“奋斗者”号全海深载人潜水器，具备海底自主避碰和贴海底自动匹配地形巡航等功能^[1]，提高了载人潜水器的智能程度。

然而，在载人潜水器航行的过程中，会遇到海山和深海热液喷口等障碍，可能还会受到外界各种不确定性因素的干扰，如海洋暗流、海洋生物群落等^[2]。在复杂的水下环境中，依据潜航员自身经验，在综合考虑路径长度和路径安全性等因素后，快速规划出一条安全无碰撞的路径有不小的难度。故而，在获取全局环境信息后，快速地自动完成载人潜水器的全局路径规划，能够辅助潜航员进行驾驶与决策，提高载人潜水器的智能程度。因此，复杂水下环境的智能路径规划技术是一个有意义的研究课题。

目前有多种全局路径规划算法，如粒子群算法、蚁群算法、A*算法、RRT算法等。冯炜等^[3]提出了基于量子行为的粒子群优化的路径规划方法，考虑了海流和障碍物的影响，但是其优化目标是航行时间代价，未考虑能量消耗代价。占银^[4]基于改进蚁群算法对水下机器人路径规划进行了研究，仿真结果表明该算法规划的路径更短、更平滑且收敛速度更快，但仅研究了二维平面内的路径规划。Garau等^[5]采用A*算法在海流环境下为AUV规划了运动路径。虽然 A*算法的执行效率高，但由于A*算法是启发式算法，无法保证所规划路径的全局最优性。RRT算法可以在短时间内规划出一条安全无碰撞的路径^[6]，Carreras等^[7]采用 RRT*算法在线进行二维 AUV 路径搜索，但其随机性强，往往难以得到全局最优路径。

为此，本文针对载人潜水器三维路径规划问题展开研究，提出一种基于改进人工蜂群算法的三维全局路径规划方法。在算法初始化阶段，采用混沌映射序列对蜜源进行初始化；为了跳出局部最优，采用高斯变异对部分适应度较差的个体进行扰动。最后，通过仿真证明了改进算法的有效性。

1 模型建立 1.1 三维海底模型

海洋环境模型的建立是进行路径规划的基础。基于全球海陆数据库(general bathymetric chart of the oceans，GEBCO)的开源海底地形数据对海洋环境建模。从GBECO地图上选定目标海域，下载后即可得目标海域范围内海底各部分的高度数据。参照二维栅格模型的建立方法，将其推广到三维环境中，三维海洋栅格模型示意图如图1所示。

图1中，x轴和y轴代表海域的长度和宽度，z轴代表海底地形高度。从原点开始，将oxyz空间按一定的长宽高分割为一个个小块，然后根据获取的海底地形高度数据值对每一个小立方体进行判断，当立方体栅格处的高度值小于或者等于该位置所对应的海底地形高度值时，则将其视为障碍栅格，载人潜水器无法通行（见图1中灰色立方体），否则视为可自由航行的栅格。

1.2 海流模型

海流是载人潜水器航行的重要影响因素。在地球自转运动的影响及海洋水平尺度远大于垂直尺度的限制下，海流在水平面内的运动规模远大于垂直面，因此对问题进行简化，将海流的运动看作二维水平面内的运动，并假设海流在目标海域中不同点上的幅度和速度方向时不变。

通常海流运动可以由二维Navier-Stokes方程进行判断和预测，然而实时计算不仅复杂而且成本高，而二维Navier-Stokes方程的近似解可以构造为粘性Lamb涡流，水平面内的海流运动可用多个粘性Lamb涡流叠加表示，其数学描述如下^[8]：

$ {V_x}({\boldsymbol{r}}) = - \varGamma \frac{{y - {y_0}}}{{2{\text{π}} {{({\boldsymbol{r}} - {{\boldsymbol{r}}_0})}^2}}}\left[1 - {e^{ - \left(\frac{{{{\left({\boldsymbol{r}} - {{\boldsymbol{r}}_0}\right)}^2}}}{{{\delta ^2}}}\right)}}\right]，$

(1)

$ {V_y}({\boldsymbol{r}}) = \varGamma \frac{{x - {x_0}}}{{2{\text{π}} {{({\boldsymbol{r}} - {{\boldsymbol{r}}_0})}^2}}}\left[1 - {e^{ - \left(\frac{{{{({\boldsymbol{r}} - {{\boldsymbol{r}}_0})}^2}}}{{{\delta ^2}}}\right)}}\right]，$

(2)

$ w({\boldsymbol{r}}) = \varGamma \frac{1}{{{\text{π}} {\delta ^2}}}{e^{ - \left(\frac{{{{({\boldsymbol{r}} - {{\boldsymbol{r}}_0})}^2}}}{{{\delta ^2}}}\right)}} 。$

(3)

式中： ${{\boldsymbol{r}}_0}$ 为涡流中心的位置向量； $\delta 和\varGamma$ 代表涡流半径和涡流强度； ${V}_{x}(r)和{V}_{y}(r)$ 分别代表横向、纵向的涡流速度分量； $ w({\boldsymbol{r}}) $ 代表涡量。式(1)和式(2)中的参数设置如下： $\delta = 1.0$ ， $\varGamma {\text{ = 10}}$ ，将海流平面设置为长宽均为100 km的二维平面，将其按照1 km×1 km的大小进行分割，涡流中心矢量 ${{\boldsymbol{r}}_0}$ 随机取海流平面内的10个点。由上述涡流方程可得海流模型如图2所示。

1.3 路径规划策略

给定起点和终点，再加入若干个控制点，使用直线进行连接便可以得到一条折线路径，这条不与障碍相碰撞的折线就是路径规划问题中的可行解。然而这样一条折线路径在实际航行过程中不够光滑，不具有可行性，因此本文使用三次均匀B样条曲线来对路径进行平滑处理。

k阶B样条函数表达式为：

$ P(t) = \mathop \sum \limits_{i = 0}^n {B_{i,k}}(u){P_i}。$

(4)

式中： ${P_i}(i = 0,1,2,3\cdots )$ 表示控制点； $ {B_{i,k}}(u) $ 表示 $k$ 次B样条基函数。目前广泛使用的B样条基函数是由DeBoor^[9]和 Cox^[10]分别提出的递推定义：

$ \left\{\begin{array}{l}{B}_{i,0}(x)=\left\{\begin{array}{l}1,{u}_{i}\leqslant x\leqslant {u}_{i+1}，\\ 0，{\rm{otherwise}}。\end{array}\right.\\ {B}_{i,k}(x)=\dfrac{x-{u}_{i}}{{u}_{i+k}-{u}_{i}}{B}_{i,k-1}(x)+\\ \qquad\quad\;\;\; \dfrac{{u}_{i+k+1}-x}{{u}_{i+k+1}-{u}_{i+1}}{B}_{i+1,k-1}(x),k > 0 ，\\ 约定\text{0/0=0}。\end{array}\right.$

(5)

式中：u为节点；下标i为B样条的序号；k为B样条的幂次。

对于均匀B样条曲线，因节点等间距分布，对于任意i恒有：u_i+1−u_i=c，c为常数。故经过推导计算可得三次均匀B样条基函数如下：

$ \left\{ \begin{gathered} {B_{0,3}}(u) = \frac{1}{6}{(1 - u)^3} ，\\ {B_{1,3}}(u) = \frac{1}{6}{(3{u^3} - 6{u^2}{{ + 4}})^3}，\\ {B_{2,3}}(u) = \frac{1}{6}{({{ - 3}}{u^3} + 3{u^2} + 3u + 1)^3} ，\\ {B_{3,3}}(u) = \frac{1}{6}{u^3} ，\\ \end{gathered} \right.\quad u \in [0,1]，$

(6)

因此，可得三次B样条曲线段如下：

$ \begin{split}{P_{i,3}}(u) = &{P_i}*{B_{0,3}}(u) + {P_{i + 1}}*{B_{1,3}}(u) + \\ &{P_{i + 2}}*{B_{2,3}}(u) + {P_{i + 3}}*{B_{3,3}}(u),u \in [0,1]。\end{split}$

(7)

在给定控制点P_i后，使用式(7)便可求得一条由若干满足三次B样条曲线点所组成的光滑曲线路径。该方法把三维空间路径规划问题转化为控制点搜索问题，降低了问题求解难度。

1.4 路径的代价函数

从给定起点R₀开始，经过中间路径点R₁,···,R_n-1，最终到达终点R_n，该路径可表示如下：

$ \begin{split} & {\rm{Path}} = \{ {{\boldsymbol{R}}_0},{{\boldsymbol{R}}_1},\cdots,{{\boldsymbol{R}}_n}\} ，\\ & {{\boldsymbol{R}}_i} = ({x_i},{y_i},{z_i})。\\ \end{split} $

(8)

对于所有可行路径，需要建立一个路径评价模型对路径质量进行评估，在路径长度的基础上，加入地形代价和路径能耗这2个指标，得到路径规划的综合代价函数，其表达式如下：

$ W = {P_{\text{c}}} + {S_{\text{c}}} + {E_{\text{c}}}。$

(9)

式中：P_c为路径长度；S_c为地形代价；E_c为能耗代价。

1.4.1 路径长度Pc

路径长度P_c表示载人潜水器从起点到终点的行驶距离，使用欧几里得距离度量，其计算过程可表示如下：

$ {P_c} = \mathop \sum \limits_{i = 1}^{n - 1} {D_i} = \mathop \sum \limits_{i = 1}^{n - 1} ||{{\boldsymbol{N}}_{i + 1}} - {{\boldsymbol{N}}_i}||。$

(10)

式中：

$ ||{{\boldsymbol{N}}_{i + 1}} - {{\boldsymbol{N}}_i}|| = \sqrt {{{({x_{i + 1}}{\text{ - }}{x_i})}^2} + {{({y_{i + 1}} - {y_i})}^2} + {{({z_{i + 1}}{\text{ - }}{z_i})}^2}}。$

(11)

式中：D_i为当前路径点N_i和下一路径点N_i+1之间的距离；n为路径点总个数。

1.4.2 地形代价Sc

地形代价S_c表示载人潜水器在航行过程中来自海山等静态障碍物的威胁，通过它，可以使载人潜水器避开航行路径中的障碍，地形代价S_c可表示如下：

$ \left\{ \begin{gathered} {S_c} = \mathop \sum \limits_{i = 1}^n {S_{{c_i}}}，\\ {S_{{c_i}}} = \left\{ \begin{gathered} {\text{M}},{z_i} < {z_2}({x_i},{y_i}) ，\\ 0,{\rm{otherwise}}。\\ \end{gathered} \right. \\ \end{gathered} \right. $

(12)

1.4.3 能耗代价E_c

能耗代价E_c表示载人潜水器沿规划路径的航行过程中为克服海流影响消耗的总能量，计算表达式如下：

$ {E_c} = \mathop \sum \limits_{i = 1}^{n - 1} {W_i}，$

(13)

式中，W_i为第i段路径的能耗。为了计算这一项，首先对潜器在海流环境下的航行速度进行合成，其航行速度合成如图3所示^[11]。

图中，V为载人潜水器在静态海洋中的航行速度；V_oc为海流速度；V_h代表合成后载人潜水器的速度； ${\theta _{{oc}}}$ 和 ${\psi _{{{oc}}}}$ 分别为海流速度在水平和垂直方向上的夹角； $\theta $ 和 $\psi $ 分别为V_h在水平和垂直方向上的夹角；V_h与V_oc的夹角为Ф。由图可得，V_oc在3个方向上的速度分量 $({u_{{{oc}}}},{v_{{{oc}}}},{w_{{{oc}}}})$ 可表示如下：

$ \left\{ \begin{gathered} {u_{{{oc}}}} = |{{\boldsymbol{V}}_{{{oc}}}}|\cos {\theta _{{{oc}}}}\cos {\psi _{{{oc}}}} ，\\ {v_{{{oc}}}} = |{{\boldsymbol{V}}_{{{oc}}}}|\cos {\theta _{{{oc}}}}\sin {\psi _{{{oc}}}}，\\ {w_{{{oc}}}} = |{{\boldsymbol{V}}_{{{oc}}}}|\sin {\theta _{{{oc}}}}。\\ \end{gathered} \right. $

(14)

载人潜水器航行速度在3个方向上的速度分量 $(u,v,w)$ 和夹角Ф可表示如下：

$ \left\{\begin{array}{l}u=\left|V\right|\mathrm{cos}\theta \mathrm{cos}\psi +{u}_{{oc}}，\\ v=\left|V\right|\mathrm{cos}\theta \mathrm{sin}\psi +{v}_{{oc}}，\\ w=\left|V\right|\mathrm{sin}\theta +{w}_{{oc}}，\\ \varPhi=\mathrm{arccos}\left(\dfrac{({u}_{{oc}},{v}_{{oc}},{w}_{{oc}})\cdot {({R}_{i}{R}_{i+1})}^{{\rm{T}}}}{\left|{V}_{{oc}}\right|\left|{R}_{i}{R}_{i+1}\right|}\right)。\end{array}\right. $

(15)

结合式(14)和式(15)，推导出合成后载人潜水器的速度V_h如下：

$ \left|{V}_{{h}}\right|=\left|\sqrt{{V}^{2}-{({V}_{{oc}}\mathrm{sin}\Phi)}^{2}}+{V}_{{oc}}\mathrm{cos}\right| ，$

(16)

参照Alvarez等^[12]的研究，路径中第i段的能量消耗W_i可表示如下：

$ \left\{ \begin{gathered} {J_i} = \iiint\limits_{{R_{i - 1}}{R_i}} {||{{\boldsymbol{V}}_h}(x,y,z)|{|^3}{\rm{d}}x{\rm{d}}y{\rm{d}}z} ，\\ {W_i} = \frac{{\rho {J_i}}}{{{V_{{{oc}}}}}}。\\ \end{gathered} \right. $

(17)

式中：ρ为依赖于潜器尺寸和水流特性的常数；V_oc为潜器的静水速度值。结合式(13)可计算出路径的总能耗代价E_c。

2 人工蜂群算法 2.1 基本人工蜂群算法

人工蜂群算法是一种模拟蜜蜂采蜜的群智能优化算法，由Karaboga于2005年提出^[13]，其主要思路是通过蜜蜂不同个体的分工及个体间的信息交流协作寻找高质量的蜜源。人工蜂群算法包括蜜源、采蜜蜂、观察蜂和侦察蜂4个要素，以及种类转化、放弃蜜源和搜索新蜜源3种基本行为。人工蜂群算法示意如图4所示。

一个蜜源代表所求问题的一个解。初始时，蜜源 ${{\boldsymbol{g}}_i} = ({g_i}_1,{g_i}_2,\cdots ,{g_i}_d)$ 可由下式随机生成：

$ {g_i}_j = {L_{ij}} + {r_1} \times ({U_{ij}} - {L_{ij}})。$

(18)

式中： ${g_i}_j$ 表示第 $i$ 个蜜源的第 $j$ 维变量， $j \in \{ 1,2,\cdots,d\} $ ， $d$ 为问题维数； $i \in \{ 1,2,\cdots,N\} $ ，N为蜜源个数； ${U_{ij}}$ ， $ {L_{ij}} $ 为 ${g_{ij}}$ 的上界和下界；r₁为(0,1)内的随机数。

初始化N个蜜源后，采蜜蜂在第i个蜜源周围依照下式搜索产生一个新的蜜源：

$ {{\boldsymbol{g}}_{new}} = {{\boldsymbol{g}}_i} + {r_2} \times ({{\boldsymbol{g}}_i} - {{\boldsymbol{g}}_k})。$

(19)

式中： ${{\boldsymbol{g}}_n}_{ew}$ 为搜索得到的新蜜源；r₂为随机数； $k \in \{ 1,2,\cdots,N\} $ ，且 $k \ne i$ 。通过比较新旧蜜源的适应度大小，基于贪婪原则保留更优质的蜜源。在人工蜂群算法中，蜜源的适应度值计算方法如下：

$ fi{t}_{i}=\left\{\begin{array}{ll}\dfrac{1}{1+re{s}_{i}}\text{，}&re{s}_{i}\geqslant 0，\\ 1+abs(re{s}_{i})\text{，}&re{s}_{i} < 0。\end{array}\right. $

(20)

式中， $re{s_i}$ 为所求解问题的函数值。

在采蜜蜂搜索完成后，飞回蜂巢跳摆尾舞与其他蜜蜂共享蜜源信息。观察蜂根据蜜源适应度，采用轮盘赌法选择接下来要搜索的蜜源，按照下式计算对蜜源i的选择概率：

$ {p_i} = \frac{{fi{t_i}}}{{\displaystyle\sum\limits_{i = 1}^N {fi{t_i}} }}，$

(21)

在搜索过程中，如果经过t次迭代搜索到达陷入局部最优后的次数判断阈值Jgelim后仍未找到更优质的蜜源，则会放弃该蜜源。同时与之对应的雇佣蜂会转换为侦察蜂，侦察蜂会在搜索空间中随机搜索产生一个新的蜜源，该过程可表示为如下：

$ {{g}_{i}}^{t+1}=\left\{\begin{array}{ll}{L}_{i}+{r}_{1}\times ({U}_{i}-{L}_{i})\text{，}&t\geqslant J{\rm{gelim}}，\\ {{g}_{i}}^{t}\text{，}&t < J{\rm{gelim}}。\end{array}\right. $

(22)

2.2 改进人工蜂群算法 2.2.1 蜜源初始化

蜜源的初始化对求解的速度和质量有很大影响。在基本人工蜂群算法中，蜜源的位置通过随机初始化产生，这样会导致蜜源在整个目标空间中的位置分布不均匀。而基于混沌理论的混沌序列具有全空间遍历性和动态性，因此使用混沌序列对人工蜂群算法的蜜源初始化进行改进以克服蜜源随机初始化时存在的分布不均匀的不足。

采用Logistic映射式产生混沌序列，公式如下：

$ {x_k} = \mu {x_{k - 1}}(1 - {x_{k - 1}}) 。$

(23)

式中： $\ \mu $ 为增长率，当 $\ \mu = 4$ 时，系统处于全混沌状态； ${x_0}$ 为初值， ${x_k}$ 指算法在第k次迭代时的值，根据混沌理论，当初值 ${x_0}$ 不等于0，0.25，0.5，0.75，1.0时，序列处于混沌态，可以扩大算法的搜索范围。

在引入混沌序列后，可对式(18)进行改进得到下式：

$ {{\boldsymbol{g}}_i} = {{\boldsymbol{L}}_i} + {{\boldsymbol{x}}_i} \times ({{\boldsymbol{U}}_i} - {{\boldsymbol{L}}_i})。$

(24)

式中： $i \in \{ 1,2,...,N\} $ ，N为蜜源个数； ${{\boldsymbol{U}}_i}$ ， $ {{\boldsymbol{L}}_i} $ 为 ${{\boldsymbol{g}}_i}$ 的上界和下界； $ {{\boldsymbol{x}}_i} $ 为迭代过一定次数后的d维混度序列。

2.2.2 蜜源搜索机制

为了提高人工蜂群算法的局部搜索能力，引入高斯变异机制，在每次迭代后对蜜源的适应度函数值进行排序，选择最差的 $n \times \eta $ 个蜜源，进行高斯变异，其中 $n$ 是蜜源个数， $\eta $ 为变异的比例，依据经验， $\eta $ 取为1/12； $\mu $ 为均值， $\sigma $ 为标准差。

最差的 $n \times \eta $ 个蜜源进行高斯变异后，新蜜源产生公式如下：

$ {{\boldsymbol{g}}_{new}}^t = {{\boldsymbol{g}}_i}^t + Gauss(\mu ,{\sigma ^2})。$

(25)

式中， ${{\boldsymbol{g}}_i}^t$ 和 ${{\boldsymbol{g}}_{new}}^t$ 分别为变异前和变异后的蜂群位置。

改进人工蜂群算法流程如图5所示。

3 仿真及结果分析 3.1 试验环境及算法参数设置

将式(9)作为优化目标为潜器搜索全局最优路径。采用的仿真平台为Matlab 2021a。基于GEBCO数据库，选取14.7601°N～15.8°N，158.4009°E～159.4428°E作为试验的目标海域建立三维海底地形模型，目标海域大小为100 km×100 km。取海域中最深为Z=0的点，其余部分的高度值根据该值确定。起点为(1,1,50)，终点为(100,100,80)。

选取基本人工蜂群算法、遗传算法和粒子群算法与改进人工蜂群算法进行对比。这4种算法均采用相同的最大迭代次数IterMax=100, 这些算法的其余参数如表1所示。

3.2 试验结果

使用改进人工蜂群算法、基本人工蜂群算法、遗传算法和粒子群算法进行三维路径规划，最终规划得到的路径如图6～图8所示。

可以看出，采用这4种算法为潜器所规划的路径均能避开环境中的障碍物，并成功到达终点。使用上述4种算法均分别迭代100次后，所规划的路径长度见表2，所得到的适应度函数值随着迭代次数变化曲线如图9所示。

从表2可以看出，改进人工蜂群算法相比于其他3种算法能够规划出一条更短的全局最优路径。从图9可以看出，改进人工蜂群算法、人工蜂群算法和粒子群算法均能较快速地收敛，而遗传算法的收敛速度则相对较慢。随着迭代次数的增加，相比于其他3种算法，改进人工蜂群算法能够不断跳出局部最优解，最终收敛在一个更小的适应度函数值。相比之下，改进人工蜂群算法的路径规划效果最好。

4 结　语

本文提出一种基于改进人工蜂群算法的载人潜水器路径规划方法，使用三维栅格模型为海洋环境建模，综合考虑了路径长度、地形代价和能耗代价三方面因素建立了路径评价模型。基于Logistic混沌映射在算法的初始化阶段对蜜源进行初始化，使用高斯变异对每次路径搜索中适应度较差的一部分个体进行扰动，提高了其局部搜索能力，并采用三次 B 样条曲线对路径进行平滑处理。仿真试验表明，相比于其他群智能算法，本文提出的改进算法收敛较快，跳出局部最优的能力更强，可以高效地为载人潜水器规划出一条满足性能要求的最优路径。