考虑船位预测不确定性的船舶碰撞危险度计算方法

王少博，张英俊，胡鑫

大连海事大学 航海学院，辽宁 大连 116026

0 引 言

智能船舶的核心功能之一是自主航行能力。船舶在海上航行过程中，首先需要收集当前所处的航行场景信息；然后判断本船与他船之间是否存在碰撞危险，如果存在，则分析两船间的会遇态势并制定避让决策[1-2]。由此可见，船舶碰撞危险度的计算占有重要地位，它作为中间环节，直接决定了两船间由自由航行局面转变为会遇局面的时机。传统航海中，一般采用自动雷达标绘仪（automatic radar plotting aid, ARPA)实现对碰撞危险的预警，即通过设定最近会遇时间（time to the closest point of approach, TCPA)与最近会遇距离(distance of the closest point of approach, DCPA)的安全阈值来实现对潜在碰撞危险的探测。但是，阈值的设定并没有统一标准，如果设定值过大，则会造成频繁的碰撞预警；如果设定值过小，则会影响船舶航行安全[3]。

多年来，国内外学者对船舶碰撞危险度计算模型进行了许多研究，根据研究思路的不同，大致分为2 类：一是指标因素法，二是船舶领域法。指标因素法就是采用多种典型的航海指标对船舶潜在的碰撞危险进行评判。Hilgert 等[4]基于两船间的最近会遇距离将两船之间的碰撞危险划分为4 个等级；吴兆麟等[5]和郑中义等[6]提出了时间与空间船舶碰撞危险度计算模型；严庆新[7]应用模糊数学综合评判理论，把目标船的距离、相对方位、最近会遇距离、至最近点的时间、船速比、碰角等6 个基本参数作为碰撞危险度的判断参数；刘茹茹等[8]在DCPA 和TCPA 的基础上加入了船舶安全距离圈和最晚施舵时间等主观因素；Bukhari 等[9]为了减轻船舶交通管理(VTS)值班员的工作负担，提出了一种综合考虑罗经方位变化、DCPA 和TCPA 的船舶碰撞危险度计算模型；文元桥等[10]选取船舶会遇方式、风速和能见度作为控制变量，运用序次Probit 模型拟合出航道船舶碰撞危险度的状态方程；汤国瑞等[11]针对复杂水域船舶碰撞风险较高的问题，分别选取航速、航向、距离、能见度、风、浪、流等要素，设计了碰撞危险度的计算方法。而在船舶领域法方面，Fujii 等[12]首先提出了船舶领域的定义，它指的是绝大多数后继船舶驾驶员避免进入前一艘在航船舶周围的领域；Goodwin[13]和Coldwell[14]分别根据来船方位和会遇局面构造了不同的船舶领域模型；Szlapczynski 等[15]在船舶领域的基础上提出了计算碰撞危险的2 个参数，即领域侵入程度(degree of domain violation, DDV)和领域侵入前剩余时间(time to domain violation, TDV)；Im 等[16]提出了基于领域的潜在风险等级划分模型；周丹等[17]使用主成分回归算法和粗糙集算法分析了影响船舶领域因素的重要性。

通过对上述船舶碰撞危险度研究现状的分析发现，指标因素法使用的多种指标会共同影响船舶碰撞危险计算的结果，各指标权重的确定受人为主观因素影响较大，计算过程不如船舶领域法直观；而在船舶领域，大部分研究仍集中在如何进一步完善船舶领域模型，对船舶领域在计算碰撞危险度中的具体应用描述较少。考虑到智能船舶在海上航行过程中需要更加可靠且稳定的碰撞危险度计算模型，因此，本文将基于Szlapczynski等[15]提出的DDV 和TDV 的概念，分析基于船舶领域的碰撞危险参数计算方法，对船位预测不确定性进行建模，提出考虑船位预测不确定性的碰撞危险度计算模型。

1 基于船舶领域的碰撞危险参数计算模型

1.1 参数定义

学者们构造了多种形状的船舶领域，其中大部分为多边形和椭圆形的船舶领域[18]。本文采用藤井模型[12]，即以本船位置为中心的椭圆形领域，其中椭圆的长半轴a等 于4 倍船长，短半轴b等于1.6 倍船长，如图1 所示。

图1 椭圆船舶领域Fig. 1 Elliptical ship domain

根据Szlapczynski 等[15]提出的DDV 和TDV概念，TDV 被定义为本船侵入他船领域之前剩余的时间，如果TDV 有合理值，则表示DDV 存在。本文用K表示DDV，用T表示TDV。当本船侵入他船领域之后，原椭圆领域长短轴乘以缩放因子r，使得缩小后的领域边界刚好经过本船，在本船侵入到离开他船领域期间，取最小缩放因子rmin，则K可以表示为

定义R为离开目标船领域的时间，则领域侵入时间 δ为

以上参数的定义如图2 所示。图中，OS 为本船，TS 为他船，两船处于对遇状态，OS 从t0时刻起沿相对运动速度矢量方向运动，在t1时刻到达TS 椭圆领域的边界，并在t2时刻达到椭圆的最小缩放因子rmin，在t3时 刻离开椭圆领域，其中t0到t1时 间为T，t0到t3时 间为R，t1到t3为 领域侵入时间δ。

图2 参数示意图Fig. 2 Schematic diagram of parameters

根据以上参数的定义，可以得到参数变化对碰撞危险的影响（表1）。

表1 参数变化影响Table 1 Influence of parameters variation

1.2 场景分割模型

以上3 个参数的计算都是面向单船领域，考虑船舶在海上航行的真实场景：在开阔水域，船舶分布比较分散，一般只考虑本船与他船两船之间的碰撞危险；而在狭窄水域或繁忙水域，经常会出现部分船舶分布较为密集的情况。由于这些船舶距离较近，互相之间存在一定的干扰。因此，本文尝试将分布密集的船舶作为一个船舶群组，从整体角度考虑本船相对于此船舶群组的碰撞危险。

根据当前海面所有他船的位置关系，首先采用层次聚类法[19]对船舶进行聚类。根据层次聚类法的基本原理，通过设定距离阈值，可以得到当前场景下的船舶层次聚类结果如图3(a)所示，聚类过程如图3(b)所示，其中蓝色表示本船，黑色表示他船。

得到虚拟船中心坐标之后，获取Uts中所有船舶领域组合构成的最大横向长度Llat和最大纵向长度Llong，并以此构造外围矩形，则此矩形的内切椭圆存在且唯一，之后构造与此内切椭圆具有相同离心率的矩形外接椭圆E，则此外接椭圆的方程存在且唯一，表达式为

图3 层次聚类结果图Fig. 3 Hierarchical clustering results

该椭圆即为以虚拟船中心O(Ox,Oy)构建的新椭圆领域。新椭圆领域的构建过程如图4 所示。图4 中3 艘船舶构成一个船舶群组，每一艘船的领域范围用黄色椭圆表示，虚拟船位置用红色十字星表示，航向用红色箭头表示，根据式(4)可得组合后的新椭圆领域为图中右侧的黑色实线椭圆。

至此，将本船当前航行场景分成了本船与他船之间的单船会遇局面以及本船与他船群组的会遇局面。无论哪一种情况，每一艘船舶或者每一船舶群组都会有围绕自己的领域。假定所有船舶均保向保速航行，如果以未来是否侵入对方领域作为判断碰撞危险出现的指标，则可分为4 种场景[18]，如图5 所示，其中蓝色表示本船，黑色代表目标船或目标船群组。

图4 新组合领域构建过程Fig. 4 Construction process of new ship domain

图5 领域侵入的4 种场景Fig. 5 Four scenarios of ship domain invasion

4 种场景分别描述如下：

1） 本船与他船均未侵入对方领域；

2） 本船侵入他船领域，他船未侵入本船领域；

3） 他船侵入本船领域，本船未侵入他船领域；

4） 本船与他船互相侵入对方领域。

设定当本船侵入他船领域时的各参数用下标α 表示，当他船侵入本船领域时的各参数用下标β表示，可以分别得到4 种场景下他船相对于本船的碰撞风险参数表达式。对于场景1），由于双方均未侵入对方领域，因此K1= 0，T1没有合理值，δ1= 0；对于场景2），可以按照参数定义得到K2=Kα，T2=Tα， δ2= δα；对于场景3），尽管本船并未侵入他船领域，但是对于他船而言，认为与本船之间存在碰撞危险，设定本船椭圆领域的长轴为lOS，他船或他船群组的椭圆领域长轴为lTS，各参数定义如下：

对于情形4），本船与他船互相侵入对方领域，各参数表达式为：

图6 船位预测偏差Fig. 6 Deviation of ship position prediction

2 船位预测不确定性模型

上文中给出的描述船舶碰撞危险的3 个参数K,T和 δ都是基于本船和他船当前的运动参数计算得到，它们均假定当前海上场景下的所有目标在未来一段时间内的运动状态保持不变。但在实际航行中，尽管本船与他船保向、保速航行，但由于受到海上若干不确定性因素的影响，比如海面气象水文条件、船舶的操纵性能等，实际船位并非完全等于预测船位，总是存在一定的偏差。因此，定义初始船位为Xo(xo,yo)，理论预测船位为Xt(xt,yt) ，实际船位为Xr(xr,yr)，则位置偏差为

如图6 所示，将实际的船位偏差量沿船舶椭圆领域的横纵向进行分解，偏差量的最终结果，可以看作船舶在航行过程中，沿椭圆领域横纵向航行时，速度矢量随时间累积的误差值。

由此，分析得到关于实际船位偏移量 ∆D的2 点性质：

1） 船位偏移量的大小与船舶沿椭圆领域长短轴向的速度分量直接相关；

2） 船位偏移量的大小与累积时间成正比，偏移增量仅与某一时间段的区间长度有关，而与初始时刻无关。

考虑到上述船位偏移量的2 点性质，用随机过程中经典的维纳过程[20]来对船位预测的不确定性进行近似建模。维纳过程 {W(j),j≥0}属于典型的独立增量过程，对于 ∀j＞s≥0，有增量：

式中：j和s代表2 个不同时间点；M为船位误差的协方差矩阵。维纳过程是一个马尔科夫过程，对未来船位的预测所需的信息就是该过程当前的数据值。维纳过程在任何有限时间上的变化服从正态分布。其中N(0,M)是由椭圆长短轴2 个方向组成的二维正态分布，其均值和协方差矩阵分别为：

式中， (µ1,σ1)和 (µ2,σ2)分别为沿椭圆长、短轴方向的船位偏移量的均值和标准差。k[x,y] 表示沿椭圆长、短轴方向船位偏移量组成的二维随机变量。由此，可以得到维纳过程遵循的二维正态分布概率密度函数:

该函数图像示意图如图7(a)所示，将图像投影到XOY 平面上可以得到一个概率椭圆，如图7(b)所示。其中，椭圆的中心坐标为 (µ1,µ2)，椭圆方程可表示为：

式中， ε为概率椭圆的规模， ε的选择代表所选取的置信水平。由于船位偏移量在椭圆领域长、短轴方向上的分布均为零均值的正态分布，且互相关系数为零，因此式(11)等号左侧实际上代表了2 个独立正态分布的随机变量的平方和，其分布规律服从卡方分布(Chi-Square)[21]。根据二自由度卡方分布数据表，选取6 个置信等级，可以得到置信度与椭圆长短轴长度的对应关系（表2）。

将椭圆船舶领域与表示船位预测不确定性的置信椭圆进行叠加，即将原椭圆领域的长短半轴分别与置信椭圆的长短半轴相加，可以得到融合船位不确定性置信度的新椭圆领域。如图8 所示，其中l表示船长。

3 考虑船位预测不确定性的碰撞危险参数计算方法

船位预测不确定性椭圆的大小与预测时长成正比，令 σ1(t)与 σ2(t)分别为沿椭圆领域长、短轴方向船位偏移量的标准差关于预测时间t的函数，可以得到融合船位不确定性置信度的新椭圆领域：

图7 概率密度函数图Fig. 7 Image of the probability density function

表2 置信度与椭圆长短轴的对应关系表Table 2 Table of correspondence between confidence level and ellipse's major-minor axis

图8 叠加领域模型Fig. 8 Superposition ship domain model

式中：a和b分别为椭圆领域的长、短半轴长度；ω为对应置信系数下不确定椭圆长、短半轴的乘积因子，如在表2 中，当 γ为10%时， ω=0.46。

图9 参数计算示意图Fig. 9 Schematic diagram of parameter calculation

根据旋转平移变换可得到目标船在初始位置P时的椭圆领域方程为：

该椭圆中心沿相对运动线移动，则目标船的位置方程为：

当达到T时，本船刚好过目标船椭圆领域边界，由于本船位于坐标原点 (0,0)点，结合式(12)可得T的解为式(15)中t的计算结果：

式(15)的计算结果如果存在正实数根，则说明T有合理值，即本船侵入了他船领域，如果存在2 个不同正实数根t1和t2， 且t1＞t2， 则t1为R，t2为T，由式(2)可得领域侵入时间 δ为：

椭圆领域的缩放因子r是关于时间t的方程r(t) ，其中t∈(t2,t1)，满足式(17)：

整理后可得r(t)表达式为：

4 结 语

本文首先将海上会遇局面分为单船会遇以及本船与船舶群组会遇的局面，并根据条件“是否侵入对方领域”将会遇局面分为4 种场景，提出了基于船舶领域侵入程度、剩余侵入时间以及领域侵入时间的船舶碰撞危险参数计算模型；同时，根据维纳过程对船位预测不确定性进行建模；最后给出了考虑船位预测不确定性的碰撞危险参数的计算方法。本文考虑了船位预测不确定性给船舶碰撞危险带来的影响，提出的模型能够为智能船舶碰撞危险分析等相关领域的研究提供参考。