内河水域沿岸经济繁荣[1]，航运需求旺盛，交通密度持续增加，增大了船舶碰撞等事故概率[2]，事故易造成航道堵塞、经济损失甚至人员伤亡[3]．内河水域通常通航条件复杂，航道狭窄，船舶自动化、信息化条件有限．但长江干线江苏段等水域，航道宽阔且水流平缓、水深条件较好、航行规则完善，研究这类水域(称为“内河宽水域”)航行辅助或自主航行决策方法，具备一定可行性．航行辅助及自主航行技术是保障安全、提高船舶航行自动化、减少人为事故的有效手段[4]．船舶航行决策研究领域成果丰富．文献[5]提出了一种能够动态适应系统中残差和目标船舶随机操纵的智能航行决策方法．文献[6]提出一种基于船舶运动数学模型的开阔水域船舶自主避碰方法，将模型预测控制算法和模糊数学理论结合，构建避碰算法的评价函数，实现船舶自主避碰[7]．基于深度强化学习(DRL)的船舶智能避碰方法可实现内河避碰规则下对目标船和其他碍航物的避让[8]．文献[9]研究了一种基于双延迟深度确定性策略梯度(TD3)的船舶自主避碰方法．文献[10-12]基于改进的人工势场法，提出了能适应复杂动态环境下的自主船舶避碰方法；在分析不同会遇情境下船舶自主避碰决策航向和航速变化的耦合机制的基础上，补充和优化了船舶自主避碰决策理论．文献[13]提出了一种基于船舶操纵过程演绎的自主船舶动态自适应避碰决策支持方法．文献[14]基于Coldwell船舶领域改进了多船避碰危险度评估模型．融合《1972年国际海上避碰规则》和良好船艺的直航船智能避碰决策方法[15]，可实现复杂水域下的自动航行．以浏河口至江阴大桥段分道通航制水域为例，尚有以下关键问题亟待解决：满足决策需要的特殊水域数字化交通环境构建；《中华人民共和国内河避碰规则》《长江江苏段船舶定线制管理规定》和良好船艺(合称“航行规则”)的量化解析；特殊环境下耦合本船变速变向运动，速度矢量与避让效果之间的关系规律(避碰机理)；融入量化解析航行规则的决策框架、航向航速控制和碰撞危险量化判断等航行决策方法．围绕这些问题，基于非线性船舶操纵运动模型及航向航速控制方法，在构建数字化交通环境的基础上，提出针对该类水域的航行决策方法．1 动态自适应航行决策框架船舶自主航行中环境感知是航行决策的关键．环境感知是监测动静态物标运动状态，为决策提供基础信息，但感知误差、船舶运动模型的不确定性及控制参数等因素会导致航行决策存在误差．通过选用合适的数学模型、提升传感器精度和优化控制器参数可缩小误差，但误差仍存在．这种误差定义为剩余误差，可设计自适应反馈校正框架补偿．自适应航行决策流程见图1．决策框架根据t时刻的航行信息计算t+Δt时刻的航行决策，确保Δt时间内航行安全，每个Δt周期内的决策会根据实时航行信息进行更新，补偿Δt内决策误差．Δt维持在一个适宜的范围内，即使存在误差，通过连续的决策更新，误差可保持在可控范围．根据以上思路，构建时序滚动自适应航行决策框架．以本船和交通环境信息为输入，本船时序操纵方案为输出．以步长Δt获取动静态目标信息，∆t取1 s[5]，实时判断其运动状态，分析航行态势，反馈调节、执行操纵方案，主要步骤如下．a．构建数字化交通环境．b．输入目标船、其他物标、航道及本船运动状态等信息．c．基于船位推算方法推算本船、目标船船位，判断是否存在碰撞危险，若存在危险则依据“航行规则”约束下避碰机理确定操纵方案；不存在危险则根据航线跟踪方法计算目标航向和目标车速，并输入航向航速控制推演船舶运动过程，确定航线跟踪是否安全，是则执行航线跟踪方案，否则执行当前方案．d．更新信息，重复步骤b～d．10.13245/j.hust.241281.F001图1航行决策框架图上述步骤根据时序滚动计算方法，快速更新航行信息，补偿剩余误差对操纵方案的影响．在实践中对不按“航行规则”行动的船舶，本船也会立即根据其动向更新决策方案，直至本船安全．2 数字化交通环境2.1　数字化交通环境模型通过解构、提取研究水域的动、静态环境要素，建立满足航行决策需要的数字化交通环境模型，数字化交通环境具体分类见表1．10.13245/j.hust.241281.T001表1物标分类环境要素类型要素模型要素状态浮标、明(暗)礁点状静态浅滩、岸线多边形静态常规机动船普通机动船动态拖带船、拖网渔船条形动态操限船、失控船圆形动态非机动船点状动态研究水域由M个分道通航制，s个沿岸通航带，G个航标及W个不规则区域(如航道附近的锚地、码头和汊河口)构成．每个分道通航制由n个通航分道和D条边界线组成，每条边界线含有e个点．分道通航制、通航分道、边界线、边界线包含的点、通航分道划分的航段按由西向东、由北向南依次递增．某个分道通航制的d条边可由向量Z表示，Z=(Z1d,Z2d,…,Zed)，（1）式中Zed为某个分道通航制的d条边上第e个点．分道通航制的n个通航分道可由向量C表示，C=[Bdn,Bd+1n]T，(2)式中n与d的对应关系为d=2n-1．计划航线由通航分道两条边界线的中点依次连接组成．Lij为第i个分道通航制的第j个通航分道上的计划航线，Lij=line|P1j,P2j,…,Pgj，(3)式中：line为点的连线；pgj为第i个分道通航制的第j个通航分道的两条边界线上第e个点连线的中点，其中g为点的个数．根据通航分道中船舶位置和流向的不同，将某个通航分道划分为多个航段．Qnw为第j个通航分道上的第w个航段，Qnw=Fed      Fe+1dFed+1   Fe+1d+1．(4)将航标处理成一定半径的圆形区域，MI=[area (x+rcos ∂),(y+rsin ∂)]，(5)式中：(x,y)和r为区域圆心和半径；∂为参数且∂∈[0,2π]．汊河口、沿岸通航带、隔离带、附近锚地、码头、危险标志等形状不一且不规则，用多个点连线表示．第a个不规则区域为Ua，Ua=[p1,p2,…,po]，(6)式中：po为组成该区域的点；o为点的数量．2.2　数字化交通环境仿真浏河口至江阴大桥段数字化交通环境如图2所示．将数字化交通环境模型应用到研究水域中，10.13245/j.hust.241281.F002图2浏河口至江阴大桥段数字化交通环境M=3，G=43，W=4，D=6．Z，C，Lij，Qnw，MI，Ua均可用通用模型表示，其参数取值如下：当m=1，2，3时，n=4，灯浮半径r=15 m．第a个不规则区域Ua，当a=1，2，3，4时，j=4；当a=5时，j=9；当a=6时，j=6．3 航向航速控制系统与航线跟踪方法3.1　船舶操纵运动模型考虑外部环境因素如风、浪和流的干扰，建立船舶三自由度运动模型如下    (m+my)v˙+(m+mx)ur=YH+YP+YR+Ywind +Ywave ;    (m+mx)u˙-(m+my)vr=XH+XP+XR+Xwind +Xwave ;    (Izz+Jzz)r˙=NH+NR+NP+Nwind +Nwave ,(7)式中：m，mx和my分别为船舶质量以及在x轴和y轴方向产生的附加质量；Izz和Jzz分别为船舶惯性矩和附加惯性矩；u，v和r及其导数分别为船舶在x轴、y轴方向上的速度和加速度及其角速度和角加速度；X，Y和N分别为船体、桨和舵在x轴和y轴方向上所受的力与力矩；下标H，P，R，wind和wave分别表示船体、桨、舵、风和流．当船舶匀速直线运动时，有效推力表达式为Xp=(1-tp0)L，(8)式中：L为螺旋桨推力；tp0为推力减额系数，与船型、螺旋桨尺度、负荷有关．tp0受船舶运动状态及螺旋桨负荷影响，tp0=0.04+(tp0'-0.04)u/u0    (u≤u0);tp0'    (uu0),(9)式中：u0和u分别为额定、瞬时速度；tp0'为额定速度时的减额系数，受旋转和横向运动影响．又有f=tp0'-tp0，(10)式中f为减额系数的影响函数．风流求解公式见文献[16]．3.2　自适应最优航向控制方法传统比例-积分-微分(PID)控制的参数根据经验设置，无法随外界环境实时调整，因此引入最优控制理论，自适应调整PID参数．为求取船舶操纵性指数K和T，将式(7)经拉普拉斯变化得到二阶响应模型：    (m+mx)su(s)=[Xr-(m+my)v0]r(s)+Xuv(s)+Xδ(s);    (Izz+Jzz)sr(s)=Nu(s)+Nrr(s)+Nδδ(s).(11)经过拉普拉斯逆换算，频域内转首、横向加速度和舵角的关系可转化为时域方程：T1T2r¨+(T1+T2)r˙+r=K(δ+T3δ˙);T1T2u¨+(T1+T2)u˙+u=Ku(δ+T3uδ˙),(12)式中K，T1，T2，T3均为船舶操纵性指数，T=T1+T2-T3．统筹考虑航向保持和节省燃油，选取性能优化函数J=(φe2+λδ2)/2,(13)式中：φe为航向偏差；λ为加权系数，其取值与风速大小相关．根据式(13)，通过水动力导数可求出K和T，具体计算方法详见文献[17]．为使式(13)二次型性能优化函数值最小，根据线性二次型非定点输出理论，可解得舵角的控制律δ=1λ1K(1-1+2KTλ)EEC,则KP，KI，KD计算公式为：KP=1/λ;KD=(1-1+2KT/λ)/K；KI=K/(Tλ)/10．根据上述所求的PID参数实现航向自适应最优控制．根据实验船舶在浏河口至江阴大桥段航行的航速数据，可计算出随航速V实时变化的K和T值，变化曲线如图3所示．10.13245/j.hust.241281.F003图3K和T值随航速变化曲线3.3　分档车钟航速控制方法船舶通过改变螺旋桨转速来实现变速，在实践中螺旋桨转速微小变化对避碰效果影响不大．选用分档车钟来发出转速指令．考虑“航行规则”中有航速限制，故只考虑船舶前进转速．以“华洋理想”轮为例，车钟转速表见表2．10.13245/j.hust.241281.T002表2车钟转速表车令转速/(r∙min-1)航速/kn车令转速/(r∙min-1)航速/kn海上速度9314.0后退三8512.0前进三8512.0后退二7010.5前进二7010.5后退一485.0前进一485.0微速后退353.0微速前进353.0停车00.0若t1时刻发出车令，完成变速时间为t2，变速前后螺旋桨转速分别为P1和P2，则t时刻螺旋桨转速表达式为Pt=P1    (tt1);P2    (tt2);P1+(t-t1)k    (t1tt2,P1P2);P1-(t-t1)k    (t1tt2,P1P2),式中k为螺旋桨转速变化率，根据该轮船长经验取k=0.25 r/s．3.4　航线跟踪方法为解决航线跟踪问题，引入视线(line-of-sight，LOS)制导算法，算法原理如图4所示．通过推导船舶在跟踪期望航向点时的目标航向，结合自适应航向控制方法，控制船舶沿计划航线跟踪．为便于计算船舶航向，以坐标系y轴为正北，x轴为正东，船舶航向即本船船首向与y轴的夹角．10.13245/j.hust.241281.F004图4LOS制导示意图由图4可知期望航向点PLOS=(xLOS,yLOS)与船舶在期望航迹上的投影点相距nLpp,Lpp为船长，n=2～5．Δ为船舶的可视距离，船舶期望航向角为αϕ=αk+arctan(-e/Δ)，式中：αk为大地坐标系正北方向与期望航迹线pkpk+1的夹角；e为船舶的横向跟随误差．若船舶位置坐标为(x0,y0)，此时跟踪的期望航向点为Pk=(xk,yk)，由上述制导算法得到期望航向角αϕ，船舶航向角偏差ψe=αϕ-ψ，ψe将传输至自适应最优航向控制系统，计算出对应的舵角δ，传递至船舶运动模型．同时更新此时船舶航向和位置，并反馈给制导算法，依靠闭环的数据迭代，使船舶稳定跟踪航线，跟踪方法原理如图5所示．10.13245/j.hust.241281.F005图5航线跟踪方法原理图当ψe趋近0°时，船舶即向目标航向点航行；当船舶距下一航向点小于设定的距离R时，船舶将自动切换航向点跟踪．4 航行决策方法4.1　避碰机理4.1.1　基于船位推算的危险度模型为了实时推算本船与目标船的位置，考虑到沿通航分道航行船舶航迹向与推荐交通流方向一致时，一般会继续沿着航道航行．若航迹向与航道走向有明显差别，或者是航道外船舶，则按照当前航速矢量推算船位．本船若因避让或其他操作偏离计划航线，则当条件允许时，应尽快调整航向，恢复与推荐交通流同向航行．本船按照船舶运动模型推算，目标船具体推算方法如下：x=x0+vtsin(TC)    (t∈[0,tj)),y=y0+vtcos(TC)    (t∈[0,tj));x=x1+v(t-tj)sin αi+1    (t≥tj),y=y1+v(t-tj)cos αi+1    (t≥tj);x=x0+vtsin(TC),y=y0+vtcos(TC).将航道边界离散化为若干个点，航标视为静态碍航物．在本船运动过程中，若目标船或其他碍航物某一时刻进入本船领域，则视为存在潜在碰撞危险(potential collision risk，PCR)．碰撞危险度指数u(collision risk index，CRI)是反映碰撞危险是否存在和采取避让措施紧迫程度的物理量，可分为空间碰撞危险度指数us(space collision risk index，SCRI)和时间碰撞危险度指数ut (time collision risk index，TCRI)[18]，其表达式为u=usut．4.1.2　会遇局面辨识模型研究水域的“航行规则”除规定存在碰撞危险时两船在避碰行动中的责任和义务外，还对汊河口有特殊规定．船舶是否处于汊河口，可由下式判断o=(Xk2-Xk3)(Yt0-Yk3)-(Yk2-Yk3)(Xt0-Xk3);p=(Xk1-Xk4)(Yt0-Y4)-(Yk1-Yk4)(Xt0-Xk4);q=(Xk2-Xk1)(Yt0-Yk1)-(Yk2-Yk1)(Xt0-Xk1);r=(Xk3-Xk2)(Yt0-Yk2)-(Yk3-Yk2)(Xt0-Xk2),式中：(Xki,Yki)为研究水域汊河口第k个航段坐标；(Xt0,Yt0)为t时刻本船船位，当且仅当o，p，q，r同时大于或小于0时，船舶位于汊河口航段内．为明确船舶之间的会遇局面，通过对“航行规则”和良好船艺的解析，提出了基于交通流向与船舶航向之间差值判断船舶会遇局面的方法，并构建了该研究水域的会遇局面辨识模型，方法如下．当存在危险时，本船和目标船航向分别为φt0和φti，对应的船舶所在航段主流向为φm0和φni，其中m和n为上下行航路编号．a．若φt0-φm0≤5°，φti-φni≤5°且m=n，则此时两船为追越局面．b．若φt0-φm0≤5°，φti-φni≤5°且m≠n，则此时两船为对驶局面．c．若φt0-φm0≤5°，85°≤φti-φni≤95°，则此时两船为横越局面．d．若目标船位于汊河口且φti-φni≤5°，则此时两船为汊河口危险局面．e．若不构成追越、对驶、横越、汊河口局面，则此时两船为交叉相遇局面．4.1.3　规则融入与解析本船为研究水域顺交通流方向行驶的机动船，规则要求如下：a．船舶应在规定的通航分道或航路内行驶；b．进出汊河口、支流及专用航道的船舶，应主动避让在规定的通航分道、推荐航路内正常行驶的船舶；c．通航分道内航速较低的大型船舶，应尽可能沿通航分道右侧边缘行驶，在保障安全的情况下，可进入推荐航路航行．为确定会遇局面，引入上述局面辨识模型，按照会遇态势确定两机动船的避让关系如图6所示．10.13245/j.hust.241281.F006图6避让原则流程图4.1.4　可行操纵决策方案求取船舶可采用改向、变速、变速和改向相结合等操纵方式避让．若存在碰撞危险，则计算改向角度和变速档位．避让行动优先顺序为转向、变速、转向和变速相结合．可行操纵决策求取方案如下．方案1 单凭转向方案．当根据“航行规则”单凭转向时，将Caim 输入航向航速控制系统，依次验证是否目标船及其他碍航物(航道外边界、航标、浅滩等)会进入本船领域．此时Caim 即为单凭转向的目标航向，α为最小改向角．若所有改向角都不可行，则采取方案2；若船舶领域到达航道边界线，航向更新为当前航段的航道边界线方向．有Caim =C0±α    (1°≤α≤30°)，式中：Caim 和C0为目标航向角和初始航向角；α为改向角，α=1°～30°．方案2 单凭变速方案．将船舶车令档位Naim 和航线跟踪方法中目标航向Caim 输入航向航速控制系统，逐步验证各车令档位是否可以避让，档位幅度变化最小的方案即为单凭变速的最优避让方案．若所有档位都不可行，则采用方案3．有Naim =N0±ΔN    (Naim ∈[Ns,NF])，式中：Naim 为目标转速；Ns和NF为本船停车、全速前进时的螺旋桨转速；N0和ΔN为本船初始转速和转速改变量．方案3 转向变速方案．根据方案1和方案2中所列本船的Caim 和螺旋桨转速Naim ，从较小改向角改变车令档位，将Caim 和Naim 输入航向航速控制系统，依次验证是否目标船及其他碍航物进入船舶领域．改向幅度和档位变化幅度最小的组合即为最优避让方案．若未搜寻到有效方案，则采取人工干预．4.2　航行决策方法在时变的航行环境中，以Δt为周期实时获取目标船和环境信息．通过构建动态自适应航行决策框架，搜寻最优决策方案，并循环执行，如图7所示．10.13245/j.hust.241281.F007图7航行决策方法5 实验与分析5.1　仿真环境设置风、浪、流等外界气象条件会直接影响航线跟踪精度与避碰决策方案的求取．取2023年9月17日14:00时刻浏河口至江阴大桥段真实风、流数据．风数据为风向313°，风速2.1 m/s；流数据为流向150°，流速0.2 m/s．本船吃水为10.5 m，在上述环境中进行旋回实验，如图8所示．10.13245/j.hust.241281.F008图8风、流环境下旋回实验在上述环境下，设置目标船舶参数见表3．实验中除目标船9和目标船10为锚泊船外，其余均为机动船．机动船中除目标船18外，其他均保向保速．10.13245/j.hust.241281.T003表3仿真实验船舶初始数据船舶类别位置坐标/(°)航向/(°)航速/kn船长/m本船(121.330，31.551)31812.5225目标船1(121.320，31.562)3162.080目标船2(121.314，31.573)3263.080目标船3(121.259，31.646)2267.085目标船4(121.253，31.588)706.0130目标船5(121.335，31.562)2295.8135目标船6(121.325，31.582)2289.775目标船7(121.319，31.590)2375.895目标船8(121.267，31.649)15613.0110目标船9(121.296，31.620)2300.0120目标船10(121.293，31.628)2320.0145目标船11(121.259，31.646)2267.0165目标船12(121.243，31.656)2266.0160目标船13(121.225，31.664)2255.3125目标船14(121.241，31.637)1236.080目标船15(121.141，31.674)339.4105目标船16(121.140，31.701)1225.0110目标船17(121.121，31.702)354.0100目标船18(121.144，31.718)2047.01205.2　仿真实验根据表3的初始数据在Python平台进行模拟仿真．图9(a)、(b)展示了追越避让过程，图9(c)、(d)展示了交叉避让过程，图9(e)、(f)分别显示了目标船机动下危险的识别与避让过程．10.13245/j.hust.241281.F009图9仿真实验结果图10为仿真实验本船的航速与航向变化曲线，该场景下操纵决策方案如下．10.13245/j.hust.241281.F010图10仿真实验过程航速航向变化图当t1=0 s时，本船航向由318.4°左转5°(车钟前进三)追越目标船1，t2=286 s时刻开始跟踪航线．当t3=1 595 s时，本船减速避让(降速至前进二)目标船11；当t4=1 766 s时，本船(车钟前进二)按照推荐交通流方向左转9°；当t5=2 005 s时，开始跟踪航线(车钟前进三)．当t6=3 188 s时，变速至前进二避让目标船15；当t7=3 685 s时，开始跟踪航线(车钟前进三)．当t8=4 063 s时，目标船18开始运动；当t9=4 126 s时目标船18开始机动转向；当t10=4 194 s时，本船识别危险并判定为追越局面，本船航向由312.7°左转6°(车钟前进三)避让目标船18；当t11=4 393 s时开始跟踪航线．通过仿真实验验证了提出的航行决策方法的有效性，能够安全避让所有目标并及时跟踪航线．操纵决策方案经过了穷举验证和重复计算，但在考虑环境、船舶操纵性、航行规则和船艺约束的基础上，事先排除了不符合航海实践和海员通常做法的方案，以降低计算时间复杂度．利用Python编译器可实现方案的实时计算，计算时间约为0.1 s左右．6 结语针对内河宽水域分道通航制下的航行决策问题，通过构建静态数字化交通环境，建立航向航速控制系统；结合水域环境、“航行规则”与良好船艺等约束下的避碰机理，依靠时序滚动的闭环思路，提出一种能动态自适应系统剩余误差和目标船随机运动的内河宽水域变速航行决策方法；以浏河口至江阴大桥的深水航道为例进行了仿真实验，结果表明该方法满足在此类水域安全航行的需求，可为该类水域自主航行提供理论和技术支持．但研究中未考虑极端情形下的“倒车”和“停车”及变速与变向不同步时的耦合和浅水效应等问题，将在后续研究中进一步完善．