0 引　言

随着超大型油船及水下潜航器的出现和研制先进船舶运动控制器的需要，船舶运动数学模型研究工作受到越来越多的关注^[1-2]，现代控制理论、先进测量技术和系统辨识理论等进一步推动了船舶运动数学模型的发展^[3-4]。从模型结构的角度来看，船舶运动数学模型可分为水动力模型和响应型模型，而水动力模型又分为以Norrbin为代表的整体型模型和以MMG（Ship Manoeuvring Mathematical Model Group）为代表的分离型模型。Norrbin模型在分解流体动力的过程中更注重参数演绎的物理含义，给出了一种描述非线性流体动力的简洁形式^[5]，该成果在瑞典SSPA早期研制的船舶操纵模拟器中广泛应用。

近年来，随着船舶日益大型化、快速化、专业化和现代化，20世纪80年代建立的船舶运动整体型或分离型数学模型对于大型船舶（通常指8万载重吨以上或总长250 m及以上的船舶）来说其精度越来越低，不再适合航海仿真和船舶运动控制器设计的新要求。为了解决大型船舶非线性Norrbin整体型数学模型精度不高的难题，文献[6]针对压载排水量20万吨的矿砂船Vale Brasil旋回试验进行研究，在调节部分流体动力导数与Norrbin辨识的2个非线性力（矩）经验公式的基础上，对原Norrbin船舶运动数学模型进行改进，仿真结果优于改进前，与实船数据较为接近。本文在其基础上针对30万吨级船舶的非线性Norrbin模型精度展开后续研究。文献[7-8]采用MMG模型建模思想为满载30万吨级大型油船KVLCC2号建立了较精确的数学模型，但其建模过程中所需船舶参数繁多，且公式较为复杂。对于一般货轮来说，收集众多船舶参数并不容易，故精度高的MMG模型难以构造，而发展一套新的建模方法有很大难度。本文立足于已有的非线性Norrbin数学模型，挖掘其与船舶大型化密切相关的流体动力导数，对其进行适当放大，并在不改变Norrbin辨识的非线性力（矩）公式的基础上，给出30万吨级船舶的改进Norrbin数学模型。用改进的Norrbin模型对KVLCC2号油船进行仿真实验，对比旋回试验与Z形试验的结果可知改进的Norrbin模型具有较高的精度，且保留了原Norrbin整体型模型所需船舶参数少、物理意义明显的优点。继而将改进的Norrbin模型应用到30万吨级不同船型的船舶旋回预报上，仿真表明具有85.1%的精度，佐证了改进后的Norrbin模型具有较好的泛化性能。

1 Norrbin数学模型改进

考虑船舶的平面运动，沿用Abkowitz的研究方案^[9]，把流体动力X，Y，N展开成Taylor级数时只保留1阶小量，另考虑到船舶左右对称性，相关的导数如 ${X_{\dot v}}$ ， ${X_{\dot r}}$ ， ${X_v}$ ， ${X_r}$ ， ${X_\delta }$ ， ${Y_{\dot u}}$ ， ${Y_u}$ ， ${N_{\dot u}}$ ， ${N_u}$ 等为0，则可得到三自由度船舶平面运动数学模型为：

$ \begin{array}{l} m\Delta \dot u = {X_u}\Delta u + {X_{\dot u}}\Delta \dot u\text{，} \\ m\dot v + m{u_0}r + m{x_c}\dot r = {Y_v}v + {Y_r}r + {Y_{\dot v}}\dot v + {Y_{\dot r}}\dot r + {Y_\delta }\delta\text{，}\\ {I_{zz}}\dot r + m{x_c}\dot v + m{x_c}{u_0}r = {N_v}v + {N_r}r + {N_{\dot v}}\dot v + {N_{\dot r}}\dot r + {N_\delta }\delta \text{。}\\ \end{array} $

(1)

由式（1）可知，横漂与转首间存在着强耦合。从航速控制的角度，前进运动与其他2个自由度上的运动相互独立，可单独考虑。船舶运动用状态空间模型描述，便于处理控制作用下船舶的多变量运动问题，对风、浪、流干扰的引入较为直接和准确。将式（1）中的后2项化成矩阵形式，再叠加上Norrbin辨识出的非线性力（矩）公式及风浪干扰，采用“一撇”系统无量纲化，则考虑横漂速度 $v$ 和转首角速度 $r$ 的非线性Norrbin数学模型呈现如下^[10]：

$ {{ I}'_{(2)}}{\dot X_{(2)}} = {{ P}'_{(2)}}{X_{(2)}} + {{ Q}'_{(2)}}U + {F'_{{\rm{NON}}}} + {F'_{{\rm{WIND}}}} + {F'_{{\rm{WAVE}}}}\text{，} $

(2)

其中： ${{ I}'_{(2)}},{{ P}'_{(2)}},{{ Q}'_{(2)}}$ 分别为惯性力导数矩阵，黏性力导数矩阵及舵力导数矩阵，其数值计算涉及10个流体动力导数 ${Y_{\dot v}},{Y_v},{N_{\dot v}},{N_v},{Y_{\dot r}},{Y_r},{N_{\dot r}},{N_r},{Y_\delta },{N_\delta }$ 的求解，Clarke通过整理大量的船模试验数据，给出了上述10个流体动力导数无量纲化后的回归公式^[10-11]； ${\dot X_{(2)}} = {\left[ {\begin{array}{*{20}{c}} v&r \end{array}} \right]^{\rm{T}}}$ 为状态变量， $U = \delta $ 为控制输入，即舵角； ${F'_{{\rm{WIND}}}} = $ $\left[ {{Y_{{\rm{WIND}}}}/\left( {0.5\rho {L^3}} \right)} \right.$ ${\left. {{N_{{\rm{WIND}}}}/\left( {0.5\rho {L^4}} \right)} \right]^{\rm{T}}}$ ， ${F'_{{\rm{WAVE}}}}$ = ${\left[ {\begin{array}{*{20}{c}} {{Y_{{\rm{WAVE}}}}/\left( {0.5\rho {L^3}} \right)}& {{N_{{\rm{WAVE}}}}/\left( {0.5\rho {L^4}} \right)} \end{array}} \right]^{\rm{T}}}$ ，具体表达式详见文献[10]； ${F'_{{\rm{NON}}}}$ 的求解采用Norrbin研究船舶参数辨识提出的简洁表达式，具体如下：

$ {F'_{{\rm{NON}}}} = \left[ {\begin{array}{*{20}{c}} {{{Y'}_{{\rm{NON}}}}} \\ {{{N'}_{{\rm{NON}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {C{f_Y}(v,r)} \\ {C{f_N}(v,r)} \end{array}} \right]\text{，} $

(3)

$ {f_Y}(v,r) = \left\{ {\begin{array}{*{20}{c}} {T \cdot r\left| r \right|[ - \displaystyle\frac{1}{{12}} - \frac{1}{{{L^2}}}{{(\dfrac{v}{r})}^2}]\text{，}\quad - \infty < - \frac{1}{L}\frac{v}{r} < - \frac{1}{2}}\text{，} \\ {T \cdot r\left| r \right|[ - \displaystyle\frac{1}{2}\frac{1}{L}\dfrac{v}{r} - \frac{2}{3}\frac{1}{{{L^3}}}{{(\dfrac{v}{r})}^3}]\text{，} - \frac{1}{2} \leqslant - \frac{1}{L}\frac{v}{r} \leqslant \frac{1}{2}} \text{，}\\ {T \cdot r\left| r \right|[\displaystyle\frac{1}{{12}} + \frac{1}{{{L^2}}}{{(\frac{v}{r})}^2}]\text{，}\quad \quad \quad \frac{1}{2} < - \frac{1}{L}\frac{v}{r} < \infty } \text{，} \end{array}} \right. $

(4)

$ {f_N}(v,r) = \left\{ {\begin{array}{*{20}{l}} {T \cdot r\left| r \right|[ - \displaystyle\frac{1}{6}\frac{1}{{{L^{}}}}(\frac{v}{r})]}\text{，} \\ {T \cdot r\left| r \right|[ - \displaystyle\frac{1}{{32}} - \frac{1}{4}\frac{1}{{{L^2}}}{{(\frac{v}{r})}^2} + \frac{1}{6}\frac{1}{{{L^4}}}{{(\frac{v}{r})}^4}]} \text{，}\\ {T \cdot r\left| r \right|[\displaystyle\frac{1}{6}\frac{1}{{{L^{}}}}(\frac{v}{r})]} \text{，} \\ - \infty < - \frac{1}{L}\displaystyle\frac{v}{r} < - \frac{1}{2}\text{，} \\ -\displaystyle \frac{1}{2} \leqslant - \frac{1}{L}\frac{v}{r} \leqslant \frac{1}{2}\text{，} \\ \displaystyle\frac{1}{2} < - \frac{1}{L}\frac{v}{r} < \infty \text{。} \end{array}} \right. $

(5)

其中：C为无量纲横流系数，取0.3～0.8； $L$ 为船长； $T$ 为船舶吃水。

假定海流干扰为定常流或具有慢时变特性，只影响船舶的平面位置和速度，不影响船舶姿态角，则有速度平衡方程：

$ \begin{array}{l} {{\dot x}_0} = u\cos \psi - v\sin \psi + {U_c}\cos {\psi _c}\text{，} \\ {{\dot y}_0} = u\sin \psi + v\cos \psi + {U_c}\sin {\psi _c}\text{。} \\ \end{array} $

(6)

其中： $\psi $ 为船舶航向； ${U_c}$ ， ${\psi _c}$ 分别为绝对流速和绝对流向，具体参见图1。运用4阶龙格-库塔算法求解微分方程组式（2）和式（6），即可求得船舶任意操舵时刻的三自由度运动数学模型。

从式（4）和式（5）可知，非线性力（矩）的求解只需要知道 $L$ 和 $T$ ，而不需其他船舶结构参数，使用方便。考虑到式（2）的求解，其等式两边的单位均为加速度的单位，而式（4）和式（5）所代表的非线性力和非线性力矩也已经化为加速度的单位，此处不需要改动，求解过程更加便捷。

根据航海实践，结合理论分析可知船舶大型化后将增大与转首角速度 $r$ 和转首角加速度 $\dot r$ 有关的流体动力导数。在基于VB与Matlab混合编程的船舶运动模型实验平台上反复测试，发现将与 $r$ ， $\dot r$ 有关的流体动力导数 ${Y_{\dot r}},{Y_r},{N_{\dot r}},{N_r}$ 均放大1.4倍，无量纲横流系数C取0.5，仿真结果与实船试验吻合度较好。

2 模型仿真实验与讨论

旋回试验与Z形试验能判定船舶操纵性优劣，被称为“标准操纵试验”，在船舶操纵试验中受到广泛的重视。本节运用原Norrbin模型、改进的Norrbin模型分别进行满载35°右旋回仿真和Z形试验，并与文献[7-8]中的实船试验、MMG模型仿真实验数据进行对比，具体描述见表1。

2.1 旋回速降

1979年278 000 t油船ESSO OSAKA号操纵性海试表明，在不改变主机转速、功率的情况下，船舶35°舵角右旋回时，其前进航速由初始动舵时的12 kn下降到定常旋回的3 kn^[12]，这一显著降速所代表的非线性已不能用 ${X_u}$ 中关于 $\Delta u$ 的三次项所描述，必须寻求替代方法。本文采用国际标准船模之一的30万吨级油船KVLCC2号^[13]作为仿真对象，由文献[7]可知其满载时右满舵35°旋回时船速的速降系数 $\mu $ 随时间 $t$ 变化，如表2所示。

使用指数下降法拟合速降公式，令响应曲线表达式为：

$ k(t) = {{\rm{e}}^{ - t/{T_0}}},\begin{array}{*{20}{c}} {}&{t \geqslant 0} \text{，} \end{array} $

(7)

由控制理论知识可知，当 $t$ 分别取0， ${T_0}$ ，2 ${T_0}$ ，3 ${T_0}$ ，4 ${T_0}$ 时，响应曲线各处斜率为：

$ \begin{array}{l} {\left. {\displaystyle\frac{{{\rm{d}}k(t)}}{{{\rm{d}}t}}} \right|_{t = 0}} = - \displaystyle\frac{1}{{{T_0}}},{\left. {\displaystyle\frac{{{\rm{d}}k(t)}}{{{\rm{d}}t}}} \right|_{t = {T_0}}} = - 0.368\displaystyle\frac{1}{{{T_0}}},\\ {\left. {{{\left. {\displaystyle\frac{{{\rm{d}}k(t)}}{{{\rm{d}}t}}} \right|}_{t = 2{T_0}}} = - 0.135\displaystyle\frac{1}{{{T_0}}},\displaystyle\frac{{{\rm{d}}k(t)}}{{{\rm{d}}t}}} \right|_{t = 3{T_0}}} = - 0.05\displaystyle\frac{1}{{{T_0}}},\\ {\left. {\displaystyle\frac{{{\rm{d}}k(t)}}{{{\rm{d}}t}}} \right|_{t = 4{T_0}}} = - 0.018\displaystyle\frac{1}{{{T_0}}},{\left. {\displaystyle\frac{{{\rm{d}}k(t)}}{{{\rm{d}}t}}} \right|_{t = \infty }} = 0\text{。} \end{array} $

(8)

若定义该响应曲线衰减到初始值2%所需的时间为调节时间，则有 ${t_{{s}}}$ =4 ${T_0}$ 。对于KVLCC2号油船旋回速降，取 ${t_{{s}}}$ =1 000 s，则有 ${T_0}$ =250 s，速降系数 $\mu $ 随时间变化的关系为：

$\mu = 0.25 + 0.75{{\rm{e}}^{ - t/250}}\text{。} $

(9)

图2所示为速降系数的拟合情况，易知式（9）拟合效果良好。

2.2 35°右旋回试验

首先运用原Norrbin模型进行无风流影响下的船舶旋回仿真实验。KVLCC2号油船满载时的主要船型参数^[7]如表3所示。考虑到舵机伺服系统的特性，将其等效为一个一阶惯性环节，其时间常数为 ${T_r}$ ，通常取7～14 s，本文取 ${T_r}$ =8 s，则舵机伺服模型为：

$ \dot \delta = - \frac{1}{{{T_r}}}\delta + \frac{1}{{{T_r}}}{\delta _r}\text{。} $

(10)

式中， ${\delta _r}$ 为指令舵角。

式（3）中无量纲横流系数C取0.5，图3给出了原Norrbin模型仿真，及按照第1节中Norrbin模型改进策略所做的仿真实验与实船试验的旋回圈对比，易知原Norrbin模型仿真的旋回圈图形符合度较低，对大型船舶来说精度有所欠缺；将与转首角速度 $r$ 和转首角加速度 $\dot r$ 有关的流体动力导数分别放大1.3，1.4，1.5倍后，其旋回圈横纵尺度均变大，更贴近于实船试验。由于放大1.4倍后的旋回圈曲线介于放大1.3倍与1.5倍之间，与实船旋回圈图形符合度较好，故选取将与 $r$ ， $\dot r$ 有关的流体动力导数 ${Y_{\dot r}},{Y_r},{N_{\dot r}},{N_r}$ 均放大1.4倍为最终改进方案。

针对改进的Norrbin模型进行仿真实验，具体数据如图4所示。虚线代表KVLCC2号实船旋回数据，点线代表MMG模型旋回仿真数据，实线代表改进的Norrbin模型旋回仿真数据，实船旋回数据与MMG模型仿真数据来源于文献[7-8]，图4（b）～图4（d）中横坐标时间 $t$ 除以了 $L/{V_0}$ ，无因次时间小刻度显示，便于观察因变量的历时变化。

图4（c）和图4（d）说明在旋回的初始阶段Ⅰ，船舶横漂速度 $v$ 及转首角速度 $r$ 急剧增大；然后在向定常回转转变的过渡阶段Ⅱ， $v$ 和 $r$ 增速变缓；最后在定常回转阶段Ⅲ时， $v$ 和 $r$ 趋于常值。表4给出了KVLCC2号各模型旋回试验的特征参数统计，其中旋回直径采用旋回圈横向与纵向的平均直径来估算，图5用柱状图的形式对比了各模型旋回特征参数与实船的差异，由柱状图可以较清晰地看出改进后的Norrbin模型旋回时在各项旋回指标上与实船较为接近。参照文献[14]中关于符合度的定义，定义平均符合度

$ \begin{array}{l} {{\bar C}_M} = \displaystyle\frac{1}{4}\left( {\frac{{\min \left( {{A_{d{\text{实船}}}},{A_{d{\text{模型}}}}} \right)}}{{\max \left( {{A_{d{\text{实船}}}},{A_{d{\text{模型}}}}} \right)}} + \frac{{\min \left( {{T_{r{\text{实船}}}},{T_{r{\text{模型}}}}} \right)}}{{\max \left( {{T_{r{\text{实船}}}},{T_{r{\text{模型}}}}} \right)}}} \right.{\rm{ + }}\\ \;\;\;\;\;\;\;\left. {\displaystyle\frac{{\min \left( {{D_{T{\text{实船}}}},{D_{T{\text{模型}}}}} \right)}}{{\max \left( {{D_{T{\text{实船}}}},{D_{T{\text{模型}}}}} \right)}} + \frac{{\min \left( {{D_{{\text{实船}}}},{D_{{\text{模型}}}}} \right)}}{{\max \left( {{D_{{\text{实船}}}},{D_{{\text{模型}}}}} \right)}}} \right) \times 100\% \text{。} \end{array} $

(11)

其中： $A_{d{\text{实船}}}$ 为实船旋回进距； $A_{d{\text{模型}}}$ 为模型旋回仿真的进距； $T_{r{\text{实船}}}$ 为实船旋回横距； $T_{r{\text{模型}}}$ 为模型旋回横距； $D_{T{\text{实船}}}$ 为实船旋回初径； $D_{T{\text{模型}}}$ 为模型旋回初径； $D_{{\text{实船}}}$ 为实船旋回直径； $D_{{\text{模型}}}$ 为模型旋回直径。于是可计算出原Norrbin模型、MMG模型、改进的Norrbin模型平均符合度分别为76.0%，92.3%，93.6%，改进的Norrbin模型旋回试验精度高于原Norrbin模型和MMG模型。

2.3 Z形试验

采用第1节中改进的Norrbin模型做10°/10°与20°/20°Z形试验仿真，因Z形试验时船舶降速较少，故此处不考虑船舶速降问题。实船数据与MMG模型仿真数据来源于文献[7-8]。图6和图7为Z形试验的仿真图，同时给出了船舶平面位置、横漂速度和转艏角速度的历时曲线，揭示了Z形试验中一些船舶参量的变化机理。在船舶由一侧舵角最大值转向另一侧的瞬间， $v$ 和 $r$ 同时达到该阶段的极值，且船舶偏离原始位置的最大横向距离时刻相应地晚于第一、第二超越角出现的时刻。表5为2次Z形试验所对应的超越角数值。仿照式（11）定义Z形试验的平均符合度 ${\bar C'_M}$ ，将第一、第二超越角的角度数值与超越角出现的时刻纳入到符合度指标，则可计算出MMG模型、改进的Norrbin模型仿真时相应的 ${\bar C'_M}$ 分别为85.2%，84.8%，两者精度相当，但改进的Norrbin模型仿真时超越角出现的时刻与实船试验更接近，即图形的相位符合度更好。

根据MSC.137（76）（《船舶操纵性标准》）规定，因本次试验 $L/V$ =39 s，则10°/10°Z行操纵试验的第一超越角应不超过20°，第二超越角应不超过40°；20°/20°Z行操纵试验的第一超越角应不超过25°。表5中各超越角数值均满足IMO《船舶操纵性标准》^[15]。

3 泛化性验证

由前文讨论可知，改进的Norrbin模型对于30万吨级船舶KVLCC2号具有较好的仿真精度。本节验证改进的Norrbin模型是否适用于排水量等级相差不大的其他船型船舶。

表6给出了大型矿砂船Stellar Topaz号的满载船型参数，图8为改进的Norrbin模型用于矿砂船Stellar topaz号旋回试验的结果对比，对应的海洋环境为风速7.6 m/s，风向040°，3级海况。考虑旋回速降时采用2.1节中的速降模型。图8中虚线为Stellar topaz号实船旋回数据，点线为原Norrbin模型旋回仿真数据，实线为改进的Norrbin模型旋回仿真数据，具体的旋回试验特征参数见表6。经计算可得改进的Norrbin模型用于Stellar topaz号旋回试验仿真平均符合度为85.1%，高于原Norrbin模型73.4%的精度。

因泛化性验证选取的船舶为不同船型船舶，充分考虑了由于船舶结构不同带来的操纵特性差异化，且仿真实验表明改进的Norrbin模型具有较高的精度，故有较好的泛化性能。

4 结　语

本研究立足于原有的非线性Norrbin模型，借鉴已有的研究成果，进一步拓展改进后的非线性Norrbin模型的适用范围。根据大型船舶的特点，挖掘出与其大型化密切相关的流体动力导数。通过在VB与Matlab混合编程的船舶运动模型实验平台反复测试，在不改变Norrbin辨识的非线性力（矩）公式的基础上，将与 $r$ 和 $\dot r$ 有关的4个流体动力导数放大1.4倍，改进的Norrbin模型与实船平均符合度高，在旋回试验上精度优于MMG模型，在Z形试验上与MMG模型精度相当，但图形的相位符合度更好。通过不同船型的泛化性验证实验证明改进的Norrbin模型具有较好的泛化性能，可应用到相近吨位级别的其他船型船舶的操纵性预报。后续可继续完成其他吨位的大型船舶的改进研究，确定参数调整的规律及其与船舶排水体积和船长的内在联系，给出通用的模型修正方法。