水上运输具有运费低、运量大的优势，在全球经济发展中作用十分重要[1]．然而，海事不断发生给人命、财产及环境带来了极大风险[2]，其中超过80%由人为因素造成[3]．自主航行有效减少人为失误，日益受到重视．通航密集区、狭水道和分道通航制等复杂水域是海事多发区域，在这些水域实现自主航行意义重大．近年来，自主航行方法方面成果丰硕．文献[4]基于改进速度障碍和人工势场法提出路径规划模型．文献[5]提出了广义速度障碍算法用于船舶智能避碰．文献[6]设计了开阔水域多船会遇局面下的动态智能航行决策方法．文献[7]基于《1972年国际海上避碰规则》(以下简称“避碰规则”)和良好船艺提出一种直航船智能避碰决策方法．文献[8]设计一种基于操纵过程推演的自主避碰决策支持模型．文献[9]基于速度障碍法提出一种最优路径和船舶运动控制的自动避碰决策模型．考虑到目标船不确定性，文献[10]提出了一种实时风险分析的时变碰撞风险模型．针对船舶发展趋势，文献[11]提出了一套人机协作的船舶避碰系统．文献[12]建立基于优化遗传算法和非线性规划理论的避碰路径规划模型．针对碰撞危险度计算精度不足和危险目标船辨识性差的问题，文献[13]基于Coldwell船舶领域改进了多船避碰危险度评估模型．定义船舶慢速域为全速与能维持航向航速间的区间．复杂环境是指航行水域内交通流密集，船舶航行受多种因素限制的通航环境．船舶在复杂水域需随时变速避让．在慢速域船舶航行方面，已有一些小型无人船变速避让决策研究，但针对大型船舶自主航行的研究大多考虑密度较低的开阔水域，倾向定速、变向避让．考虑操纵特性的慢速域自主航行成果缺乏[14-16]，很难满足大型船舶复杂水域自主航行决策要求．为实现复杂水域多船会遇自主航行决策：a．考虑到大型船舶操纵延迟性，引入MMG (mathematical model group)时增加推力随转速和航速变化的影响因子，更精确地表征了船舶慢速域变速变向操纵运动过程，当转速和航速不匹配时有很好的适用性；b．提出了符合航海实践的按车钟档位分档控制船舶速度的思想和方法，该方法容易被理解和接受，方便传统船舶升级改造；c．引入了适应变速过程中船舶操纵性能的最优PID(proportion integration differentiation)航向控制方法；d．在避碰规则和良好船艺的约束下，结合导航方法，提出了复杂水域自主航行决策方法．1 基础理论1.1　船舶会遇四阶段理论存在碰撞危险两船会遇过程可分为无碰撞危险、碰撞危险、紧迫局面和紧迫危险4个阶段，如图1所示．船舶避让责任随着阶段发展而变化[17]．为确定避让时机，根据四阶段理论建立碰撞危险度(uT)模型，具体为uT=usut，其中，uT∈[0,1]；us和ut为空间和时间碰撞危险度，分别表示两船空间碰撞的可能性和采取避让行动的紧迫性[18]．10.13245/j.hust.240406.F001图1船舶会遇四阶段理论1.2　船舶领域模型船舶领域是避免他船进入确保航行安全的区域，是避碰规则中安全距离概念的延伸[19]．根据良好船艺，不同会遇局面的安全距离不同，交叉局面下采用左偏心圆形船舶领域(图2(b))，对遇和追越局面下采用左偏心椭圆领域(图2(a)和(c))．10.13245/j.hust.240406.F002图2不同会遇局面的船舶领域a=2(LO+LT)和b=2(LO+LT)表示椭圆长和短半轴，c=0.5(LO+LT)表示虚拟船舶和实际船舶距离，其中LO和LT为本船和目标船的船长．偏心角θ为19°或161°，R=2(LO+LT)表示圆形船舶领域半径．2 改进速度障碍方法2.1　变速MMG模型避碰时须推演船舶的操纵运动过程．实践中的避让一般不考虑船舶纵摇、垂荡、横摇，选用三自由度MMG模型推演船舶运动过程，具体为：(m+my)v˙+(m+mx)ur=YH+YP+YR;(m+mx)u˙-(m+my)vr=XH+XP+XR;(Iz+Jz)r˙=NH+NR+NP, (1)式中符号含义及主要参数见文献[20]．船舶直航匀速运动时，有效推力为XP=(1-tP0)T，式中：T为螺旋桨推力；tP0为推力减额系数，与船型、螺旋桨尺度、负荷有关．对特定船舶，tP0受船舶运动状态及螺旋桨负荷影响．螺旋桨负荷变化引起前后吸入排出流变化，进而影响tP0，具体为tP0=0.04+(tP'-0.04)u/u0   (u≤u0);tP'                 (uu0)  ,  式中：u0和u分别为额定、瞬时速度；tP'为额定速度时的减额系数，受旋转和横向运动影响，即(1-tP0)=(1-tP')+f，其中f=ktβR为减额系数的影响函数，kt=2.3×10-4(γAL/Dρ)-0.028，γA=(B/d){1.3(1-Cb)-3.1lc},lc=100xc/L,βR=β-lRr/V,βR和β分别为舵和船中心处漂角，B，d，L，Dρ和xc分别为船宽、吃水、船长、螺旋桨直径和浮心纵向坐标，kt，γA，Cb，lc和lR分别为修正系数、船体整流系数、方形系数、浮心系数和漂角系数；r和V为旋转角速度和纵向速度．且T=ρNP2Dρ2KT，其中：KT为螺旋桨推力系数；ρ为水密度；NP为船舶实时转速．船舶变速通过改变螺旋桨转速实现．转速微小变化对避碰效果影响不大，大型船舶选用分档(一般10档)车钟来发出车种指令．表1为散货船“华洋理想”车钟分档与匹配航速．10.13245/j.hust.240406.T001表1车钟转速表车令转速/(r∙min-1)航速/kn海上速度9314前进三8512前进二7010.5前进一485微速前进353后退三8512后退二7010.5后退一485微速后退353停车00如车钟命令在t1发出，船舶变速前、后转速分别为N1和N2，t2为完成螺旋桨变速时间，t时刻转速为Nt=N2               (  t≥t2);N1               (t≤t1);N1-(t-t1)k              (t1tt2,N1N2);N2+(t-t1)k              (t1tt2,N1N2).根据该轮船长经验，转速变化率k取0.25 r/s．2.2　自适应最优PID控制采用自适应PID控制方法，通过工况实时调节PID参数．由式(1)拉普拉斯变化得到二阶响应模型，具体为     (m+mx)su(s)=(Xr-(m+my)v0)r(s)+  Xuv(s)+Xδδ(s);(Iz+Jz)sr(s)=Nuu(s)+Nrr(s)+Xδδ(s),式中：u(s)=L(u(t))；r(s)=L(r(t))；δ(s)=L(δ(t))，舵角δ到转首角速度r的传递函数以及到横移速度u的传递函数为H(s)=r(s)δ(s)=K(1+T3S)(1+T1S)(1+T2S)，G(s)=u(s)δ(s)=Ku(1+TuS)(1+T1S)(1+T2S)．经过拉普拉斯逆换算，频域内转首、横向加速度和舵角的关系，可转换为时域方程为：T1T2r¨+(T1+T2)r˙+r=K(δ+T3δ˙);T1T2u¨+(T1+T2)u˙+r=Ku(δ+Tuδ˙).为要大舵角避让，引入非线性一阶KT方程，具体为Tr˙+r+ar3=Kδ，其中T=T1+T2-T3．根据上式通过水动力导数求出K和T．考虑航向保持和燃油消耗性能的优化函数为J=0.5(e2+λδ2)，其中：e为航向偏差；λ为参数．根据二次型非零定点输出理论求得不同航速、螺旋桨转速下的自适应最优PID参数为：Kp=1/λ;Kd=1/K-λ+2KT/(Kλ);Ki=0.1K/Tλ.2.3　改进速度障碍方法经典速度障碍的原理如图3所示，图中：Lt和Lt'为本船到目标船的射线；Vt为相对速度．以目标船T为圆心的圆形表示船舶领域，VO和VT为本船和目标船的速度．Vr=VO-VT表示相对速度．10.13245/j.hust.240406.F003图3经典速度障碍原理当VOT落入相对碰撞区(relative collision cone，RC)时，两船会发生碰撞，即RC={Vr|Lr⋂T≠0}．绝对碰撞区(absolute collision cone，AC)为RC和目标船速度矢量合成的集合，即AC=RC⊗VT(t)，其中：⊗为闵可夫斯基矢量和运算；VO(t)和VT(t)为本船和目标船的瞬时速度，当VO(t)∈AC时两船碰撞．传统速度障碍模型不考虑船舶变速变向过程，如图4(a)红色网格区域表征传统速度障碍区间．考虑到航向和速度变化时船舶延迟性，现提出结合船舶操纵运动特性的改进速度障碍模型．复杂水域中备车航行时，为保证保向能力，驾驶员常用前进三、前进二、前进一和微速前进等4挡车钟(分别用F,H,S和D表示)控制速度．10.13245/j.hust.240406.F004图4传统速度障碍和改进速度障碍算法如图4(b)所示，圆弧表示车钟相对应的改向范围，速度障碍区间为Oi=AC⊕Ti，其中：Ti表示车钟；i∈{F,H,S,D}；Oi表示车钟i下的速度障碍区间，用不同颜色表示．3 避碰决策3.1　动态可行操纵区间动态可行操纵区间是当前时刻船舶操纵不同车钟时的可行变向区间．图5中本船和目标船存在碰撞危险，本船作为让路船，应采取碰撞行动．10.13245/j.hust.240406.F005图5特定车钟下的船舶可行航向区间根据改进速度障碍可得出，向左和向右最小可行改向角分别为C0和C1，L0和L1为对应运动航迹．如果改向角度在(不在)[-C0,C1]区间内，目标将进入(不进入)本船船舶领域．图6和图7表示不同局面船舶避碰策略示意图，红色虚线箭头表示避碰规则允许的避让方向．因此，多因素约束下包含变速和变向的动态可行操纵区间为Dt=(C1t,TFt),(C2t,TFt),...,(Cmt,TFt)(C1t,THt),(C2t,THt),...,(Cmt,THt)(C1t,TSt),(C2t,TSt),...,(Cmt,TSt)(C1t,TDt),(C2t,TDt),...,(Cmt,TDt)，式中：C表示改向角度，追越局面下，C∈(-90，10.13245/j.hust.240406.F006图6让路船不同会遇情况下的避碰策略10.13245/j.hust.240406.F007图7直航船不同会遇情况下的避碰策略90)，对遇和交叉局面下，C∈(0,90)；Dt表示t时刻的动态可行操纵区间；m，n，p和q分别表示F,H,S和D档车钟对应避碰决策方案的数量．3.2　航线跟踪船舶在偏离航线后，须回到计划航线，该过程称为航线跟踪．航线跟踪分为确定跟踪航段和跟踪决策两步．计划航线有S个转向点，依次编号为P1，P2，…，PS，计算本船和各转向点的距离．记录和距离本船最近的转向点为Pm．通过对比本船与Pm相邻的计划航段的距离来确定跟踪航段．本船与航段垂直距离为d=|Xo-Xm|     (δ=0.5π⋃δ=1.5π);|Xotanδ-Xmtanδ+Ym-Yo|tan2δ+1    (其他),式中：d表示与航段的垂直距离；δ表示计划航向；(xm,ym)为Pm的位置．两个航段和本船的垂直距离中短者为跟踪航段．跟踪决策关键在于确定目标航向(CI)，如图8所示，图中NTRUE表示真北方向．可由跟踪航段计划航向δ、终点方位(TI)等确定，即CI=λ(TI-δ)+δ，其中λ为跟踪系数，决定船舶改向幅度．系数越大(小)，转向幅度越大(小)，回到计划航线时间短(长)，但航程损耗多(少)，由航行环境决定，暂取λ=1.2．10.13245/j.hust.240406.F008图8计划航线跟踪决策3.3　优化评价函数动态可行操纵区间包含多个操纵决策，最终决策由优化评价函数决定．通过询问有经验的船长，结合实际需要，在保证安全(本研究安全的含义为保证所有目标船舶领域外通过)的前提下，选择时效性和经济性因素作为优化评价函数指标．船舶的时效性通过避让开始到返回计划航线的时间衡量．时间越短经济性越好，时效性函数为f1=ςi/maxς，式中：ςi为第i个决策时间；max ς为所有决策中的最长用时．经济性通过本船避让和航线跟踪的转向幅度表示，具体为f2=0.5Ci/maxCi+0.5CI/maxCI，其中：Ci和CI为i方案避让和航线跟踪决策改向角；max Ci和max CI为避让和航线跟踪最大改向角．优化评价函数为F=af1+bf2，其中a和b表示时效性和经济性指标权重，满足a+b=1．F越小，方案越优．3.4　自主航行决策结合以上模型，融入避碰规则与良好船艺要求的自主航行决策方法流程如图9所示．10.13245/j.hust.240406.F009图9自主航行决策流程图具体步骤如下所示．步骤1 获取本船、目标船位置、航向、航速等初始信息；步骤2 判断所有目标船碰撞危险度是否大于0，若不大于0则跳转到步骤3，否则跳转到步骤5；步骤3 判断本船是否偏离航线，若不偏离则保向保速，反之转到步骤4；步骤4 基于航线跟踪模型的决策，并根据变速MMG模型和自适应模糊PID模型推演船舶运动过程，结合目标船预测轨迹，确定回到计划航线后碰撞危险度是否大于0，否则执行航线跟踪决策，反之继续等待时机；步骤5 结合改进速度障碍求出动态船舶操纵可行区间和区间内每个元素对应的航线跟踪方案，利用优化评价函数得出自主航行方案；步骤6 并判断目标船状态是否发生变化，若变化则重复步骤1~6，反之执行原方案．4 仿真实验与结果分析4.1　运动模型的验证选择“华洋理想”号散货船作为模拟船舶．平静海面上，初始速度12.5 kn，仿真船舶右舵35°旋回轨迹如图10(a)．旋转360°需要927 s，定常旋回速度为5.75 kn，进距和横距分别为3.48Lo和1.91Lo．旋回初始速度下降迅速，定常旋回时速度恒定，如图10(b)所示，图中t表示时间．图11为不同螺旋桨转速R下变速MMG模型航速与实船航速，差值较小精度满足变速操纵推演要求．10.13245/j.hust.240406.F010图10旋回实验模拟10.13245/j.hust.240406.F011图11实船和模拟船舶的速度对比4.2　仿真实验设置两组实验验证船舶自主航行决策方法的可行性，参数如表3和表4所示．其中，实验1中目标船保向保速；实验2中目标船1和2在运动过程分别为变速变向和变向．10.13245/j.hust.240406.T002表3仿真实验1的初始数据船舶类别位置/nm航向/(°)航速/knDCPA/mTCPA/s本船(0，0)04512.5——目标船1(5.0，-1.0)32515.020.91 032.4目标船2(6.0，6.1)22510.0765.51 414.0目标船3(2.0，1.8)0456.0658.81 467.3目标船4(3.0，6.0)18111.0452.61 107.2目标船5(1.6，0.2)0308.51 194.91 007.110.13245/j.hust.240406.T003表4仿真实验2的初始数据船舶类别位置/nm航向/(°)航速/knDCPA/mTCPA/s本船(0，0)00012.5——目标船1(1.0，5.0)21014.5794.7701.2目标船2(0.0，7.5)17514.0640.21 018.8目标船3(-1.0，3.5)1559.0647.6621.0目标船4(1.7，3.2)2705.5495.0952.64.2.1　仿真实验1根据表3的初始信息，决策是120 s时右转向44°，车钟降到前进二，图中12(a)表示本船和目标船的运动轨迹．t2时刻执行航线跟踪决策，左转53°车钟回到前进3．从图12(b)和12(c)可看出目标船在一直处于船舶领域外,其中D表示两船之间的距离，C表示两船之间的最小会遇距离．图12(d)表示速度变化，在t1时刻开始避让，速度减小．t2时刻开始左转跟踪航线，阻力增大速度继续下降，t3时刻速度最低．随后船舶加速并回到计划航线．10.13245/j.hust.240406.F012图12实验1仿真结果4.2.2　仿真实验2基于碰撞危险度模型，t1=200 s存在碰撞危险，决策结果为此时右转42°避让，图中13(a)显示本船和目标船的运动轨迹．当t3=500 s时，目标船1和2分别变速变向和变向，参数见表5．本船决策方案为左转28°．t4=760 s完成转向，然后直行285 s后让清目标船，本船左转25°开始航线跟踪．图13(b)和(c)表征距离和DCPA的变化．图13(d)显示：t1时刻船速大幅度下降并在t2时刻到达最低点，t3时执行新方案，转向幅度小，降速小．t3到t4船舶直行，t5开始航线跟踪，速度逐渐增加，t6恢复正常速度．算法充分考虑了变速下的船舶操纵性，可适应目标船运动不确定性，符合现有船舶的操纵做法，容易被船员和研究人员理解，也方便在传统船舶上安装自主航行系统，推动自主航行实现．多船会遇局面仿真实验验证了自主航行决策方法的可行性和有效性．然而，变速避让范围未包括极端情形下的“倒车”和“停车”，也未考虑风浪流等影响．未来将进一步完善全车钟变速、风浪流等环境约束下的自主航行决策方法．10.13245/j.hust.240406.F013图13试验2仿真结果10.13245/j.hust.240406.T004表5目标船舶变速变向参数船舶类别位置/m航向/(°)航速/kn目标船1(-12.8，6 029.9)1908.5目标船2(313.9，12 302.0)18514.0目标船4(1 733.7，5 926.4)2708.5