船舶是水面交通运输领域最为重要的交通工具，船舶航行安全是世界船舶领域长期的研究课题，其经济性与社会意义十分重大。船舶航行常规方式是操控人员凭借驾驶经验，开展航线规划和船舶操控，规划和操控效果往往过度依赖操控人员判断能力及经验，出现差错可能性高且潜在风险大。统计数据表明：80%的船舶碰撞事故是人为因素造成的^[1]。当前船舶航行要求越来越高，只采用人为决策与操控方式，难以满足船舶航行安全性要求，因此船舶航线自动规划具有重要的理论研究及实际应用意义。

船舶航线规划方法主要有传统控制方法和智能控制方法2类。文献[2]提出一种非线性遗传算法，用以弥补优化过程局部极小问题。Mostef等^[3]综合运用分支界限法、动态规划和遗传算法设计船舶安全航行路径，利用模糊集理论设计新型的船舶避碰系统。气象环境对船舶航线规划影响较大，一些学者^[4-5]采用A^*等算法设计安全航线以规避气象风险。文献[6]针对船舶避障路径规划问题，综合考虑本船负载、目标船运动状态、航行环境以及时间等因素。Kuczkowski等^[7]创新地提出单种群和多种群遗传算法，引入船舶航行相关约束条件，实现船舶避障路径规划。Lisowski等^[8]利用最优控制理论、博弈论、人工神经网络，模拟船舶驾驶人员的思维，设计了船舶航行规划及控制算法。一些学者^[9-10]采用蚁群算法等智能控制算法设计船舶航线，并利用仿真手段加以实现。文献[11]采用模糊控制设计船舶航行避碰算法，在模糊规则设计方面较具特色。Caldwell等^[12-13]采用MPC算法开展船舶避障与航迹规划，考虑推进系统动态、风浪流作用力以及运动障碍。

以上船舶航线规划方法存在如下问题：首先，航线规划结果并未考虑船舶操纵性约束，因此按照规划结果航行时易与目标船碰撞；其次，检测规划航线与陆地元素(岛屿等)相对位置时，往往将陆地元素简化为规则的几何形状，与现实海图不符。

本文建立船舶水动力模型，利用海图数据中陆地信息构建不规则多边形，并设计基于四叉树方法的航线与不规则多边形位置关系快速检测算法，精准计算船舶航行过程与目标船会遇的时间和位置，构建综合考虑航线长度、相邻航段转角、与静态障碍物距离和动态目标船距离的复合适应度遗传算法，通过2次遗传算法设计船舶航线，实现船舶避开静、动态障碍物的任务，解决以往航线规划方法的常见问题。

1 船舶运动学建模

本文涉及到的坐标系包括墨卡托坐标系和船舶本体系，具体如图 1所示。

墨卡托坐标系x_go_gy_g：原点位于经纬度(0, 0)处，横轴x_g指向东方，纵轴y_g指向北方；船舶本体系x_bo_by_b：原点o_b位于船舶几何重心，横轴x_b指向右舷方向，纵轴y_b指向船艏方向。本文采用致密的离散墨卡托地图来描述海图上任意形状的物标。

在不考虑风浪流干扰情况下，本文构建VLCC船舶的K-T方程^[14]，以此建立船舶舵角与艏向角速度传递函数，简化表示其动态运动特性：

$ f_{K T}(s)=\frac{T}{K s+1} $

(1)

式中：T是追随性指数; K是旋回性指数。

在船舶操纵性研究的文献中，通常将回转过程的速度设置为恒定值，这与船舶实际情况不符。本文构建船舶回转降速传递函数模型，根据车钟计算船舶的航速：

$ f_{D V}(s)=\frac{R_{d}}{R_{t} s+1} $

(2)

式中：R_d为速度下降百分比；R_t为速度下降过渡时间。

2 海图物标多边形检测算法

基于四叉树方法，设计线段与不规则多边形相对位置快速检测算法，实现规划航段与岛屿关系的高效快速检测任务，具体原理如图 2所示。

本文岛屿采用不规则多边形表示，图 2(a)为一个海图文件所示的矩形区域，利用四叉树基本原理不断切分所示矩形直至规定尺寸的最小矩形，根据目标航段位置，快速确定航段穿过树形结构的矩形，并且只处理含有岛屿元素的矩形，以图 2中D₁₃为例，放大效果如图 2(b)，开展如下预处理工作：

1) 设目标航段与陆地元素的安全距离为D_safe，则四叉树最小矩形边长L_{r min}规定为：

$ L_{r \min } \geqslant D_{\text {safe }} $

(3)

2) 对原始海图中陆地元素进行插值处理，保证插值后相邻轮廓点的距离P₁P₂满足：

$ P_{1} P_{2} \leqslant 0.5 D_{\text {safe }} $

(4)

线段P_stP_et为目标航段，起点为(x_st, y_st)，终点为(x_et, y_et)，最小矩形平面内含有由若干线段收尾相连组成的岛屿轮廓，此处取其中一点岛屿轮廓点P_id=(x_id, y_id)，根据三角函数关系，可计算该点到目标航段距离D_Line：

$ \theta_{\mathrm{Line}}=\arccos \frac{\overline{P_{i d} P_{s t}} \cdot \overline{P_{i d} P_{s t}}+\overline{P_{s t} P_{e t}} \cdot \overline{P_{s t} P_{e t}}-\overline{P_{i d} P_{e t}} \cdot \overline{P_{i d} P_{e t}}}{2 \overline{P_{i d} P_{s t}} \cdot \overline{P_{s t} P_{e t}}} $

(5)

$ D_{\mathrm{is}}=\sqrt{\left(x_{i d}-x_{s t}\right)\left(x_{i d}-x_{s t}\right)+\left(y_{i d}-y_{s t}\right)\left(y_{i d}-y_{s t}\right)} $

(6)

$ D_{\text {Line }}=D_{\text {is }} \cdot \sin \left(\theta_{\text {Line }}\right) $

(7)

式中：θ_Line为$\overline {{P_{id}}{P_{et}}} $对应的三角形内角；D_is为$\overline {{P_{id}}{P_{st}}} $的长度。利用上式，可确定在目标航段穿过含有陆地元素最小矩形中陆地与目标航段最小距离。

3 航线规划目标函数

规划航线需综合考虑各影响因素，本文航线规划目标影响因素包括：航线里程、相邻航段夹角、与静态障碍物距离、与动态目标船距离。

1) 航线里程评价函数。

本文船舶规划航线里程评价函数L_traj为：

$ L_{\text {traj }}=\sum\limits_{i=1}^{n-1} \sqrt{\left(P_{x}(i)-P_{x}(i+1)\right)^{2}+\left(P_{y}(i)-P_{y}(i+1)\right)^{2}} $

(8)

式中：n为航迹点数量；P_x和P_y为各航迹点的横纵坐标。

2) 相邻航段夹角评价函数。

本文相邻规划航段的夹角评价函数θ_traj为：

$ \theta_{\text {traj }}=\sum\limits_{i=1}^{n-1} \theta_{i, i+1} $

(9)

式中：θ_{i, i+1}表示各航段夹角，该值可用平面几何关系计算。

3) 与静态障碍物距离评价函数。

本文将静态障碍物按照海图中真实数据将其表示为不规则多边形，按照第2节所描述的方法来规划航段与陆地的最小距离，评价函数D_sta为：

$ D_{\text {sta }}=\sum\limits_{j=1}^{N_{\text {st}a}}\left(\min \sum\limits_{k=1}^{N_{\text {tra }}-1} d_{j, k}\right) $

(10)

式中：N_sta为静态障碍物数量；N_tra-1为规划航段数；d_{j, k}为第j个障碍物与第k个航段距离。

4) 与动态目标船距离评价函数。

本文将目标船简化为质点处理，利用点与线段集合的几何关系来计算目标船与所有规划航段的最小距离，因此评价函数D_dyn表示为：

$ D_{\mathrm{dyn}}=\sum\limits_{j=1}^{N_{\mathrm{dyn}}}\left(\min \sum\limits_{k=1}^{N_{\mathrm{tra}}-1} d_{j, k}\right) $

(11)

式中：N_dyn为目标船数量；N_tra-1为规划航段数量；d_{j, k}表示第j个目标船与第k个航段距离。

5) 航线规划综合指标函数。

本文对以上各个目标函数进行综合加权，建立航线规划综合指标函数F_fitness：

$ F_{\text {fitness }}=\omega_{1} L_{\text {traj }}+\omega_{2} \theta_{\text {traj }}+\omega_{3} D_{\text {sta }}+\omega_{4} D_{\text {dyn }} $

(12)

式中：ω₁、ω₂、ω₃、ω₄为各指标函数的权值系数。本文制定权值系数制定原则：

1) 由于本船需避免出现与静态障碍物和动态目标船撞击情况，本文设置当D_sta或D_dyn为0时，ω₃或ω₄为极大值；

2) 本船航线里程与相邻航段夹角相比，更能够体现油耗情况，因此设置ω₁比ω₂略大。

4 考虑静态障碍物的一次航线规划

本文设计二次遗传算法寻优方式，开展航线优化工作，本节将仅考虑静态障碍物对航线规划的影响，开展船舶航线一次优化，即式(12)中ω₄=0。

4.1 遗传算法种群元素设置

定义船舶航行区域是边长为L_area的正方形，设航迹点种群总数为M_pop，种群每个个体染色体节点数为N_ind，每个染色体节点为二维变量[P_{chorm_x}, P_{chorm_y}]，表示航路点平面坐标。种群进化迭代次数为N_iter，种群进化交叉和变异概率分别为p_cross和p_mut。考虑到遗传算法寻优结果可能为局部最优情况，本文强制算法完成多次进化次数，从而避免局部最优问题，本文中用N_iter变量表示。

4.2 基于变概率策略的交叉、变异和选择操作

定义如下形式的交叉和变异概率p_cross和p_mut：

$ \left\{\begin{array}{l} p_{\text {cross }}=\frac{p_{c 0}}{N_{\text {iter }}\left(F_{\text {opt }}-F_{\text {ave }}\right)} \\ p_{\text {mut }}=\frac{p_{m 0}}{N_{\text {iter }}\left(F_{\text {opt }}-F_{\text {ave }}\right)} \end{array}\right. $

(13)

式中：p_c0和p_m0分别为交叉和变异的初始概率；F_opt为当前代最优个体适应度；F_ave为当前代的平均适应度。

1) 交叉操作。

针对当前代每个个体I_i，根据p_cross判断是否进行交叉操作，如需要，则在[1, I_i)和(I_i, M_pop]中得到均匀分布的随机数R，将个体I_i和I_R开展交叉。

2) 变异操作。

针对当前代每个个体I_i，根据p_mut判断是否进行变异操作，变异操作采用航迹点位置增加或减少方式，分别在[1, N_ind]、[0, 1]、[-1, 1]中得到均匀分布的随机数Q、增加或减少权值随机数N_rand、增加或减少概率p_mx和p_my，确定个体I_i的染色体节点Q所代表的航迹点需要变异，利用p_mx和p_my决定航迹点增加或减少，航迹点位置变异操作后坐标P_{Q_x}和P_{Q_y}计算式为：

$ \left\{\begin{array}{l} P_{Q\_{x}}=P_{Q\_{x0}}+\operatorname{sgn}\left(p_{m x}\right)\left(P_{Q\_x}-L_{\text {area }}\right)\left(1-\frac{N_{\text {cur }}}{N_{\text {iter }}}\right) \\ P_{Q\_y}=P_{Q\_{y0}}+\operatorname{sgn}\left(p_{{my }}\right)\left(P_{Q\_{y}}-L_{\text {area }}\right)\left(1-\frac{N_{\text {cur }}}{N_{\text {iter }}}\right) \end{array}\right. $

(14)

式中N_cur为在N_iter次迭代过程中的代数; P_{Q_x0}、P_{Q_y0}分别为变异前的染色体节点变量。

3) 选择操作。

计算上一代种群中最优个体染色体，建立初始可行航迹点个体I_fsb，其染色体节点为[P_{fsb_x1}, P_{fsb_y1}]，[P_{fsb_x2}, P_{fsb_y2}]，…，[P_{fsb_xN}, P_{fsb_yN}]，逐一判断当前代每个个体I_j是否存在穿越静态障碍物的情况，如存在(即D_sta=0)，则使用最优染色体I_fsb来替换I_j，否则保留个体I_j，将替换筛选后的个体组成临时种群M_{pop_tmp}，设定适应度阈值F_{fit_td}为：

$ F_{\text {fit_}td }=\left(1-\frac{1}{N_{\text {iter }}}\right) F_{j\_\min }+\frac{1}{N_{\text {iter }}} F_{j\_\max } $

(15)

式中：F_{j_min}和F_{j_max}分别表示当前j代种群M_{pop_tmp}中最小适应度和最大适应度。

5 考虑动态目标船的二次航线规划

本节在第1次航线规划基础上开展第2次航线规划任务，处理本船一次规划结果与目标船遭遇问题。

5.1 船舶会遇时间预测与安全性判定

引入船舶运动方程用以预测本船与目标船的会遇状况，预测原理如下：设规划航线为图 3中A-B-C段，本船起点为航迹点A，目的地为航迹点C，其间经过航迹点B来避开静态障碍物，第1次遗传算法寻优得到图中ADEC所示航线，船舶直航和转弯过程的航速变化情况，可简化表示为图 3中的下部所示曲线。

假设船舶初始为直航状态，速度稳定在最大航速，然后左转90°后直航，再右转90°后直航，航行过程位置和速度曲线如图 4、5所示的7个阶段：阶段1：维持最大航速；阶段2：左转过程速度下降；阶段3：由左转到直航过程，速度上升；阶段4：速度上升至最大并保持；阶段5：右转过程速度下降；阶段6：由右转到直航过程，速度上升；阶段7：速度上升至最大并保持。

根据图 4和图 5可知，在第1次规划航线的基础上，代入船舶操纵性约束后，可处理直航航段和回转航段，其中回转航段圆弧半径为本船的转弯半径决定，本文将根据船舶运动学方程(1)和(2)，分别针对直航阶段和回转阶段，计算本船由航迹点A运动到会遇点的总时间。

在图 3中，第1次规划航线AB段和BC段夹角为2θ，设本船完成该转弯过程的转弯半径为R，圆弧与AB段和BC段的切点为D和E，设航迹点A、B、C的坐标分别为(P_ax, P_ay)、(P_bx, P_by)、(P_cx, P_cy)，根据几何关系，可计算求得D和E的坐标(P_dx, P_dy)、(P_ex, P_ey) 为：

$ \left\{\begin{array}{l} P_{d x}=P_{b x}-\frac{R}{\tan \theta}\left(\frac{P_{b x}-P_{a x}}{\sqrt{\left(P_{b x}-P_{a x}\right)^{2}+\left(P_{b y}-P_{a y}\right)^{2}}}\right) \\ P_{d y}=P_{b y}-\frac{R}{\tan \theta}\left(\frac{P_{b y}-P_{a y}}{\sqrt{\left(P_{b x}-P_{a x}\right)^{2}+\left(P_{b y}-P_{a y}\right)^{2}}}\right) \end{array}\right. $

(16)

$ \left\{\begin{array}{l} P_{e x}=P_{b x}-\frac{R}{\tan \theta}\left(\frac{P_{b x}-P_{c x}}{\sqrt{\left(P_{c x}-P_{b x}\right)^{2}+\left(P_{c y}-P_{b y}\right)^{2}}}\right) \\ P_{e y}=P_{b y}-\frac{R}{\tan \theta}\left(\frac{P_{b y}-P_{c y}}{\sqrt{\left(P_{c x}-P_{b x}\right)^{2}+\left(P_{c y}-P_{b y}\right)^{2}}}\right. \end{array}\right) $

(17)

经过第1次遗传算法寻优后的实际航段由直航段AD、转弯段$\overset\frown{DE}$和直航段EC，在AD段，船舶将维持当前航速V₀，$\overset\frown{DE}$段一个降速直到稳定速度过程，EC是一个升速直到稳定过程。

船舶航行过程往往保持发动机功率，航速变化往往是改变舵角所致，本文基于船舶运动学方程开展仿真试验，获得船舶在不同舵角情况下的航速、位置和时间变化数据，仿真工况如下：船舶最大航速16 kn；舵角分别为5°、10°、15°、20°、25°。

转弯减速仿真曲线如图 6所示。按照图 6所示方式，可建立航速、舵角、时间以及航程的关系表，由于篇幅限制，此处不列出。为实现对船舶转弯降速及直航升速过程的量化计算，本文采用如下步骤加以实现：

1) 船舶转弯时，根据船舶回转跟踪的相邻航段夹角，确定船舶舵角角度；船舶即将直航时，确定船舶直航前舵角角度。

2) 根据船舶转弯或直航航程，确定船舶转弯或直航过程航速变化及用时情况，在图 6中查找船舶航速和时间。

5.2 附加航迹点的遗传算法二次规划

针对与目标船会遇导致一次规划航线不合理问题，本节提出遗传算法二次规划策略，原理如图 7。

如图 7所示，如目标船在第1次规划航线的航迹点B和C之间与本船会遇，在遗传算法中增加染色体F，以第1次规划航迹点为初始种群，开展第2次寻优以获得F坐标。为了增加第2次寻优速度，并且不对已有航迹点做较大改变，采用如下策略：

策略1：仅对由于动态目标船产生的附加航迹点及其前面1个航迹点开展动态寻优，即图中航迹点F和航迹点B；

策略2：对于航迹点B，在第1次规划航线时确定的位置上正方形范围，正方形的边长为航段AB和航段BC长度的较小长度的一半；对于航迹点F，在目标船航向的反向延长线上寻优。

6 仿真分析与验证

本文在已构建的船舶航行半物理仿真平台上开展算法验证工作，仿真平台如下图所示：由船舶动力学仿真系统、实船船桥软硬件系统、三维视景系统以及评估系统组成，以上各系统按照实船通信方式，利用网络和串口完成通信，仿真步长为0.1 s。

仿真初始工况如下：船舶航行海域为东经122.6°~123.5°，北纬28.2°~29.75°，本船起点为29.911 7°N，122.071°E，终点为29.705 6°N，122.440 7°E，本船航速16 kn。航线规划区域存在2个目标船，均直线航行，航速为12 kn，出发地分别为29.705 6°N，122.440 7°E和29.705 6°N，122.440 7°E，目的地均为29.705 6°N，122.440 7°E。遗传算法中航迹点种群变量：M_pop=100，N_ind=5，N_iter=30，N_inip=5，开展一次寻优算法航线规划。

在海图上绘制规划航线如下图所示：VLCC所在虚线为本船航线，其余虚线为目标船航线。

如图 9所示，按照本文一次规划算法可得到由起点到终点的航线，算法仅考虑规划航程、相邻航段航向转角以及与静态障碍物距离。但是由于规划航区存在2个目标船，航线规划并未考虑船舶操纵性约束，导致本船按照规划航线行驶时，与目标船“Bulk carrier”和“Merchant ship”出现会遇情况，如图 10所示，本船与目标船“Bulk carrier”存在碰撞可能性。虽然本船可通过会遇前调整航向实施避碰，但此时的航线并不是最优的，仍需改进。

按照本文第5节算法，开展第2次航线规划任务，在与“Bulk carrier”和“Merchant ship”会遇的航段中增加航迹点。在海图上绘制规划航线如下。

如图 11所示，本文在与目标船会遇的第2航段和第3航段中增加2个新航迹点，本船按照二次规划后航线行驶，利用新增航迹点可避开目标船“Bulk carrier”，具体如图 12所示。由于篇幅原因，此处不赘述避开目标船“Merchant ship”内容。

7 结论

1) 构建船舶运动学模型，可精准计算船舶在转弯和直航期间速度变化情况。引入船舶操纵性约束，与采用传统遗传算法共同构建航线规划算法，在航线跟踪过程中可避免与动态目标船出现会遇情况，提高了航行的安全性；

2) 基于四叉树思想，建立的线段与不规则多边形的位置关系检测算法，可快速准确判定规划航段与海图陆地元素的位置关系；

3) 本文利用的2次遗传算法的思想，可综合考虑船舶航程、相邻航段转向、与静态障碍物距离以及动态躲避目标船，并通过优化上述指标的适应度函数，得到可行的优化航线。

本文并未考虑船舶航线规划对海图中等深线信息，简化为与陆地元素的距离，将来将就此问题开展更进一步的研究工作。