0 引　言

船舶运动数学模型是船舶运动仿真与控制问题的核心^[1]。运动数学模型可以提供船舶操纵性限制条件，辅助人或系统进行路径规划，可以嵌入控制系统中以提高控制的精度，更是船舶运动模拟器的基础。

传统的建立船舶运动数学模型的方法是根据MMG模型结构^[2]或者Abkowitz模型结构^[3]，求解其中的水动力导数。求解水动力导数的方法中，经验公式估算十分快速，但是准确性依赖于所积累的数据和船型；约束模试验数据准确，但是成本太高，也存在尺度效应问题。系统辨识可以从包含噪声的数据中建立数学模型，并且可以不依赖于经典船舶运动数学模型结构，因此成为目前主流的建模方法。传统的辨识方法如最小二乘法、极大似然法等由于精度过度依赖模型和初值，逐渐被以神经网络为代表的新型辨识方法取代。近年来，扩展卡尔曼滤波器（EKF）^[4]、撑向量回归（SVR）^[5]、萤火虫算法^[6]、人工蜂群算法^[7]、多目标群智能算法^[8]等都被应用在了船舶系统辨识领域。以神经网络为代表的不依赖经典船舶数学模型结构的辨识方法要想得到准确的辨识结果，需要输入输出数据中包含船舶完整的操纵性信息。以EKF为代表的对经典船舶数学模型中的水动力导数进行辨识的方法由于已经引入了数学模型结构这一先验信息，往往需要更少的数据进行辨识，但辨识的精度十分依赖经典数学模型结构，并且存在参数漂移和动力相消^[9]的问题。

扩展卡尔曼滤波器利用一阶泰勒级数将非线性项线性化，因此适用于非线性系统，也可以进行参数辨识，但是因为忽略了高阶项，在非线性较强时误差大。无迹卡尔曼滤波器^[10]（Unscented Kalman Filter，UKF）是无迹变换^[11]（Unscented Transformation，UT）与标准卡尔曼滤波器的结合，目的是为了提高非线性系统滤波的精度和效率。UKF利用UT来进行均值和协方差的非线性传递，其思想是近似非线性函数的概率分布而不是近似非线性函数。UKF相比EKF有着更高的精度，并且不用计算雅可比矩阵，因此被广泛应用于导航、目标跟踪等领域。

本文采用精度更高的UKF，结合成本相对较低的自航试验对目标三体船进行辨识。针对三体船自身特点对MMG模型结构进行化简，去掉了作用不大的水动力导数，并设计了不同的控制方式进行分步辨识，以减轻动力相消的问题。通过参数的试验测量值与辨识值的对比，验证了本文辨识出参数的准确性；通过模型预报数据与试验观测数据的对比，验证了本文建立的三体船数学模型的准确性，具有较高的工程意义和参考价值。

1 无迹卡尔曼滤波器介绍 1.1 无迹变换

随机向量 $ X $ 的均值和协方差分别为 $ \bar{X} $ 和 $ {P}_{xx} $ ，对于非线性变换 $ Y=f\left(X\right) $ ，无迹变换是利用采样法求取 $ Y $ 的均值 $ \bar{Y} $ 和协方差 $ {P}_{yy} $ 。首先构造 $ Sigma $ 点集 $ \left\{{\chi }_{i}\right\} $ ，如下式：

$ {\chi _i} = \left\{ {\begin{array}{l} {\bar X,i = 0}\text{，}\\ {\bar X + \left( {\sqrt {(n + \lambda ){P_{xx}}} } \right)_i^{\rm{T}},i = 1, \cdots ,n}\text{，}\\ {\bar X - \left( {\sqrt {\left( {n + \lambda } \right){P_{xx}}} } \right)_i^{\rm{T}},i = n + 1, \cdots ,2n}\text{。} \end{array}} \right. $

(1)

式中： $ n $ 为系统的维数； $ \lambda ={\alpha }^{2}\left(n+\kappa \right)-n $ 为尺度参数； $ \alpha $ 决定了 $ Sigma $ 点围绕 $ \bar{X} $ 的分布，通常取一个小正值， $ 1>\alpha >{e}^{-4} $ ； $ \kappa $ 为第2个尺度参数，通常取为0； $ {\left(\sqrt{(n+\lambda ){P}_{xx}}\right)}_{i}^{\rm{T}} $ 表示矩阵开方运算后选择其中第 $ i $ 行并转置后得到的 $ n $ 维向量。

对 $ \left\{{\chi }_{i}\right\} $ 中的每一个点进行非线性变换 $ {Y}_{i}=f\left({\chi }_{i}\right) $ ，变换后的点集 $ \left\{{Y}_{i}\right\} $ 可以近似表示 $ Y=f\left(X\right) $ 的分布，对 $ {\{Y}_{i}\} $ 进行加权处理，得到 $ \bar{Y} $ 和 $ {P}_{yy} $ ，如下式：

$ \left\{ {\begin{array}{l} {{Y_i} = f\left( {{\chi _i}} \right),i = 0, \cdots ,2n}\text{，}\\ {\bar Y = \displaystyle\sum\nolimits_{i = 0}^{2n} {W_i^{\left( m \right)}} {Y_i}}\text{，}\\ {{P_{yy}} = \displaystyle\sum\nolimits_{i = 0}^{2n} {W_i^{\left( c \right)}} ({Y_i} - \bar Y){{({Y_i} - \bar Y)}^{\rm{T}}}}\text{，}\\ {W_0^{\left( m \right)} = \lambda /(n + \lambda )}\text{，}\\ {W_0^{\left( c \right)} = \lambda /(n + \lambda ) + 1 - {\alpha ^2} + \beta }\text{，}\\ {W_i^{\left( m \right)} = W_i^{\left( c \right)} = \lambda /\left[ {2\left( {n + \lambda } \right)} \right],i = 1, \cdots ,2n}\text{。} \end{array}} \right. $

(2)

式中： $ {W}_{i}^{\left(m\right)} $ 和 $ {W}_{i}^{\left(c\right)} $ 分别为计算均值和方差时的权重； $ \alpha $ ， $ \;\beta $ 为系数。

为叙述方便，将上述无迹变换过程用函数 $ UT $ 表示。对于随机向量 $ X $ 和一个非线性变换 $ Y=f\left(X\right) $ ，以及另一个非线性变换 $ Z=h\left(X\right) $ ，互协方差 $ {P}_{zy} $ 的计算用函数 $ {UT}_{c} $ 表示，如下式：

$ \left\{ {\begin{array}{l} {\bar Y,{P_{yy}} = UT(\bar X,{P_{xx}},f)}\text{，}\\ {\bar Z,{P_{zz}} = UT(\bar X,{P_{xx}},h)}\text{，}\\ {{P_{zy}} = U{T_c}\left( {\bar X,{P_{xx}},f,h} \right) = \displaystyle\sum\nolimits_{i = 0}^{2n} {W_i^{\left( c \right)}} ({Z_i} - \bar Z){{({Y_i} - \bar Y)}^{\rm{T}}}}\text{。} \end{array}} \right. $

(3)

1.2 无迹卡尔曼滤波器

对n维离散的非线性系统，状态空间模型如下式：

$ \left\{ {\begin{array}{*{20}{c}} {{X_k} = f\left( {{X_{k - 1}},{U_k}} \right) + {w_k}}\text{，}\\ {{Z_k} = h\left( {{X_k}} \right) + {v_k}}\text{。} \end{array}} \right. $

(4)

式中： $ {X}_{k} $ 和 $ {X}_{k-1} $ 分别为k时刻和k-1时刻的系统状态向量； $ {Z}_{k} $ 为系统观测向量； $ f $ 为状态转移函数； $ h $ 为系统观测函数； $ {U}_{k} $ 为输入控制； $ {w}_{k} $ 为服从高斯分布的过程噪声，其均值为0，协方差矩阵为 $ {{Q}}_{k} $ ； $ {v}_{k} $ 为服从高斯分布的观测噪声，其均值为0，协方差矩阵为 $ {{R}}_{k} $ 。

UKF的基本过程与标准卡尔曼滤波器一致，在误差协方差矩阵预报和更新时都用无迹变换来计算，递推公式如式（5）~式（8）所示。

状态估计预报、误差协方差矩阵预报：

$ \left\{ {\begin{array}{*{20}{c}}{X}_{k|k-1},{P}_{k|k-1}^{{'}}=UT({X}_{k-1},{P}_{k-1},f)\text{，}\\ {P}_{k|k-1}={P}_{k|k-1}^{{'}}+{Q}_{k}\text{，}\end{array}} \right. $

(5)

滤波增益：

$ \left\{ {\begin{array}{*{20}{c}}{Z}_{k|k-1},{P}_{zz}^{{'}}=UT\left({X}_{k|k-1},{P}_{k|k-1},h\right)\text{，}\\ {P}_{zz}={P}_{zz}^{{'}}+{R}_{k}\text{，}\\ {P}_{xz}={UT}_{c}\left({X}_{k|k-1},{P}_{k|k-1},f,h\right)\text{，}\\ {G}_{k}={P}_{xz}{{P}_{zz}}^{-1}\text{，}\end{array}} \right. $

(6)

状态估计更新：

$ {X}_{k}={X}_{k|k-1}+{G}_{k}({Z}_{k}-{Z}_{k|k-1})\text{，} $

(7)

误差协方差矩阵更新：

$ {P}_{k}={P}_{k|k-1}-{G}_{k}{P}_{zz}{{G}_{k}}^{T}\text{。} $

(8)

2 目标船运动数学模型构建与参数辨识方法 2.1 目标船概况

安装完成的目标三体船如图1所示，3个船体都采用了S标准型，其主要参数如表1所示。三体船的推进装置为船体两侧的1对螺旋桨，采用差分推力来控制前进和转向。

2.2 目标船运动数学模型结构

基于MMG模型建立目标三体船的运动数学模型结构，三自由度操纵性方程如下式：

$ \left\{ {\begin{array}{*{20}{c}} {\left( {m + {m_x}} \right)\dot u - \left( {m + {m_y}} \right)rv = {X_H} + {X_P}}\text{，}\\ {\left( {m + {m_y}} \right)\dot v - \left( {m + {m_x}} \right)ru = {Y_H} + {Y_P}}\text{，}\\ {\left( {{I_{zz}} + {J_{zz}}} \right)\dot r = {N_H} + {N_P}}\text{。} \end{array}} \right. $

(9)

式中： $ m $ ， $ {I}_{zz} $ 分别为船舶的质量和绕z轴的转动惯性矩； $ {m}_{x} $ ， $ {m}_{y} $ ， $ {J}_{zz} $ 分别为船舶纵向附加质量，横向附加质量，附加转动惯性矩； $ u $ ， $ v $ ， $ r $ 分别为纵向速度，横向速度，转首角速度； $ {X}_{H} $ ， $ {Y}_{H} $ ， $ {N}_{H} $ 分别为船体受到的纵向力，横向力和转首力矩； $ {X}_{P} $ ， $ {Y}_{P} $ ， $ {N}_{P} $ 分别为桨受到的纵向力，横向力和转艏力矩。

$ {X}_{H} $ ， $ {Y}_{H} $ ， $ {N}_{H} $ 是船体受到的粘性阻力和力矩，由于三体船没有横向主动力， $ v $ 的产生靠的是坐标变换的科里奥利力，即 $ r $ 和 $ u $ 是产生 $ v $ 的必要条件，所以与 $ v $ 和 $ r $ 相关的水动力导数耦合在一起，产生了动力相消的问题，无法分开辨识。另外，由于 $ v $ 很小，横向信噪比高， $ {X}_{H} $ ， $ {Y}_{H} $ ， $ {N}_{H} $ 中与 $ v $ 耦合的高阶项作用很小。因此可以在井上模型^[1]的基础上对与 $ v $ 相关的高阶项和耦合项进行化简，如下式：

$ \left\{ {\begin{array}{*{20}{c}} {{X_H} = {X_{uu}}\left| u \right|u + {X_{rr}}{r^2}}\text{，}\\ {{Y_H} = {Y_v}v + {Y_r}r}\text{，}\\ {{N_H} = {N_v}v + {N_r}r + {N_{rr}}\left| r \right|r}\text{。} \end{array}} \right. $

(10)

式中： $ {X}_{uu} $ ， $ {X}_{rr} $ ， $ {Y}_{v} $ ， $ {Y}_{r} $ ， $ {N}_{v} $ ， $ {N}_{r} $ ， $ {N}_{rr} $ 为水动力导数。

假设推力系数为进速系数的一次函数关系，如下式：

$ \left\{ {\begin{array}{*{20}{c}} {{K_T} = {a_0} + {a_1}{J_p}}\text{，}\\ {{J_p} = \dfrac{{(1 - {w_p})u}}{{nD}}}\text{。} \end{array}} \right. $

(11)

式中： $ {J}_{p} $ 为进速系数； $ {w}_{p} $ 为伴流分数；u为船速；n为螺旋桨转速；D为螺旋桨直径； $ {a}_{0} $ ， $ {a}_{1} $ 常系数。

单个螺旋桨的推力 $ T\left(n,u\right) $ 如下式：

$ \left\{ {\begin{array}{*{20}{c}} {T\left( {n,u} \right) = \left( {1 - t} \right)\rho {n^2}{D^4}{K_T} = {T_{nn}}{n^2} + {T_{un}}un}\text{，}\\ {{T_{nn}} = {a_0}\rho {D^4}(1 - t)}\text{，}\\ {{T_{un}} = {a_1}\rho {D^3}(1 - t)(1 - {w_p})}\text{。} \end{array}} \right. $

(12)

式中： $ n $ 表示螺旋桨转速； $ t $ 为推力减额； $ \ \rho $ 为液体密度； $ {T}_{nn} $ ， $ {T}_{un} $ 为系数； $ {T}_{nn} $ 也是螺旋桨系柱推力系数。

对三体船，2个螺旋桨X方向位置与船重心X位置一致，距离船中的距离为 $ l=0.45m $ ，则螺旋桨的推力和力矩 $ {X}_{P} $ ， $ {Y}_{P} $ ， $ {N}_{P} $ 的计算如下式：

$ \left\{ {\begin{array}{*{20}{c}} {{X_P} = T\left( {{n_l},u + rl} \right) + T({n_r},u - rl)}\text{，}\\ {{Y_P} = 0}\text{，}\\ {{N_P} = [T\left( {{n_l},u + rl} \right) - T({n_r},u - rl)]l}\text{。} \end{array}} \right. $

(13)

式中： $ {n}_{l} $ ， $ {n}_{r} $ 分别表示左右螺旋桨的转速。

上述模型结构中，需要确定的参数为 $ {m}_{x} $ ， $ {m}_{y} $ ， $ {J}_{zz} $ ， $ {X}_{uu} $ ， $ {X}_{rr} $ ， $ {Y}_{v} $ ， $ {Y}_{r} $ ， $ {N}_{v} $ ， $ {N}_{r} $ ， $ {N}_{rr} $ ， $ {T}_{nn} $ ， $ {T}_{un} $ 。

2.3 目标船参数辨识方法

构建系统状态向量 $ X $ 是纵向速度 $ u $ ，横向速度 $ v $ ，转艏角速度 $ r $ ，大地坐标系下正北坐标 $ x $ ，大地坐标系下正东坐标 $ y $ ，艏向角 $ \psi $ 与待辨识的参数组合而成的向量；输入控制 $ U $ 是左右螺旋桨转速组成的向量；系统观测向量 $ Z $ 是 $ X $ 前六维变量的测量值，如下式：

$ { \left\{ {\begin{array}{l} {X\! \!=\!\! {{\left[ {u,\!v,\!r,\!x,\!y,\psi ,\!{X_{uu}},\!{X_{rr}},\!{Y_v},\!{Y_r},\!{N_v},\!{N_r},\!{N_{rr}},\!{T_{nn}},\!{T_{un}}} \right]}^{\rm{T}}}}\!\text{，}\\ {U = {{[{n_l},{n_r}]}^{\rm{T}}}}\text{，}\\ {Z = {{\left[ {u,v,r,x,y,\psi } \right]}^{\rm{T}}}}\text{。} \end{array}} \right.} $

(14)

状态转移函数和系统观测函数如下式：

$ { \left\{ {\begin{array}{*{20}{c}} {f\left( {X,U} \right) = \left[ {\begin{array}{*{20}{c}} {u + \dot u \cdot dt}\\ {v + \dot v \cdot dt}\\ {r + \dot r \cdot dt}\\ {x + (u \cdot {\rm{cos}}\psi - v \cdot {\rm{sin}}\psi ){\rm{d}}t}\\ {y + (u \cdot {\rm{sin}}\psi + v \cdot {\rm{cos}}\psi ){\rm{d}}t}\\ {\psi + r \cdot {\rm{d}}t}\\ {{X_{uu}}}\\ \cdots \\ {{T_{un}}} \end{array}} \right]}\\ {h\left( X \right) = \left[ {{I_{6 \times 6}},{{\bf{0}}_{6 \times 9}}} \right]X} \end{array}} \right.\text{。}} $

(15)

其中 $ dt=0.2s $ ， $ \dot{u} $ ， $ \dot{v} $ ， $ \dot{r} $ 由式（9）计算得到。无迹变换参数取值为 $ \alpha =0.3 $ ， $\ \beta =2 $ ， $\kappa =-12$ 。 $ {Q}_{k} $ 为对角矩阵， $ u,v,r $ 分量分别取 $ {0.01}^{2} $ ， $ {0.01}^{2} $ ， $ {0.05}^{2} $ ，其余分量取0； $ {R}_{k} $ 为对角矩阵，各分量按顺序依次取 $ {0.055\;4}^{2} $ ， $ {0.055\;4}^{2} $ ， $ {0.004\;5}^{2} $ ， $ {0.2}^{2} $ ， $ {0.2}^{2} $ ， $ {0.1}^{2} $ 。

系统辨识过程中，可能会出现参数漂移或动力相消^[9]的情况，例如在式（9）中如果 $ m+{m}_{x} $ ， $ m+{m}_{y} $ ， $ {X}_{uu} $ ， $ {X}_{rr} $ ， $ {T}_{nn} $ ， $ {T}_{un} $ 同时增大同样的倍数，纵向的运动方程并不会发生变化。为避免参数漂移， $ {m}_{x} $ ， $ {m}_{y} $ ， $ {J}_{zz} $ 不参与辨识，而是由由周昭明回归公式^[12]计算得到，如表2所示。分步辨识法也可以有效地缓解动力相消问题，其思想是分批辨识参数，将前步辨识出的参数作为已知值再进行后面其他参数的辨识。对于目标船，由于对称性，在进行直航运动时 $ v $ ， $ r $ 为0， $ {X}_{rr} $ ， $ {Y}_{v} $ ， $ {Y}_{r} $ ， $ {N}_{v} $ ， $ {,N}_{r} $ ， $ {N}_{rr} $ 不产生作用，因此可以优先辨识出 $ {X}_{uu} $ ， $ {T}_{nn} $ ， $ {T}_{un} $ 这3个参数，然后将其设为已知值，再进行剩下参数的辨识。

为获取目标船参数辨识所需的观测数据序列，进行目标船静水中的自航试验。目标船上搭载了AHRS，GPS传感器，可以实时地获取位置、速度、姿态信息。目标船由岸上的遥控器进行操控，控制信号通过WIFI发送给船上的电机驱动器，以控制2个螺旋桨的转速。首先进行多组目标船直航试验，再进行多组目标船S型回转试验，将传感器观测数据记录下进行离线辨识。

3 目标船运动数学模型辨识结果

直航试验和S型回转试验螺旋桨转速信号如图2所示。将多组直航试验观测数据辨识的得到的参数 $ {X}_{uu} $ ， $ {T}_{nn} $ ， $ {T}_{un} $ 取均值，然后将其作为已知值，再利用S型回转试验观测数据对 $ {X}_{rr} $ ， $ {Y}_{v} $ ， $ {Y}_{r} $ ， $ {N}_{v} $ ， ${N}_{r}$ ， $ {N}_{rr} $ 进行辨识，结果如表3所示。

将辨识出的参数代入运动数学模型，并将模型预报的结果与实际观测数据进行比较，直航试验对比结果如图3所示，S型回转试验对比结果如图4所示。由对比结果可以得出，直航试验观测的 $ u $ 和模型预报值比较吻合；S型回转试验观测的 $ u $ 和 $ r $ 与模型预报值十分吻合， $ v $ 存在一定的误差，原因是三体船没有横向主动力，横向速度比较小，因此信噪比较大。

辨识的参数中， $ {T}_{nn} $ 是可以由螺旋桨系柱推力试验获得精确值的，为进一步验证参数辨识的准确性，在上海交通大学循环水槽进行了目标船螺旋桨的系柱推力标定试验。将左右螺旋桨转速设为一样，转速依次为5，10，15，20，25，30 r/min，测量对应的系柱推力，将系柱推力对转速平方进行线性拟合，拟合直线的斜率即为2 $ {T}_{nn} $ ，如图5所示。 $ {T}_{nn} $ 由系统辨识得到的结果与系柱推力测量的结果对比如表4所示，结果验证了辨识参数的准确性。

4 结　语

本文针对目标三体船，基于MMG模型建立了其运动数学模型结构，利用无迹卡尔曼滤波器结合自航试验数据对模型中的参数进行辨识。为减轻动力相消带来的影响，结合三体船特点对模型结构进行了化简，并设计了不同的航行方式和分步辨识方法。模型预报的速度、转首角速度与实际观测值十分吻合，对其中可以由试验测量的系柱推力系数，试验测量值与辨识值误差很小，验证了本文辨识方法的有效性和优越性。