CN108846207B - 一种基于非线性模型辨识的船舶运动预报装置及其预报方法 - Google Patents

Abstract

本发明公开一种基于非线性模型辨识的船舶运动预报装置及其预报方法，属于船舶工程、运动预报领域。本发明涉及的船舶运动预报装置包括：数据采集模块、输入处理模块、基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块、预报模块、预报误差分析模块和输出模块；其特征在于，数据采集模块、输入处理模块、基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块、预报模块、预报误差分析模块、输出模块依次串行连接。本发明涉及的船舶运动预报方法包括：数据采集；输入处理；基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正；船舶运动预报；预报效果分析；预报结果评价；预报结果输出。

Description

一种基于非线性模型辨识的船舶运动预报装置及其预报方法

技术领域

本发明属于船舶工程、运动预报领域，具体涉及一种基于非线性模型辨识的船舶运动预报装置及其预报方法。

背景技术

船舶在海浪中的运动姿态预报通常指极短期预报，或称实时预报，它是指利用某种理论和技术预测出船舶在未来较短时间内(3秒至10秒)的运动姿态。舰船运动极短期预报研究具有很高的军事和民用价值，用于舰载机起降作业指导系统，指导舰载机在飞行甲板上安全起降；用于舰船海上补给、舰船在风浪中转向时机的优选，并提供了在较高海情下实施上述作业的可能性；提高舰载武器系统、尤其是导弹的发射精度；用于船舶动力定位；有利于舰船在大风浪中航行的控制及用于舰载减摇设备减摇。另一方面通过对船舶运动姿态预报也可为航行过程中的决策提供依据。舰船平台作为武器系统的载体，是武器发挥作战能力的基础，无论是论证阶段进行的作战能力评估，还是舰载武器系统的论证、设计、仿真研究或实验研究，都需要较准确的描述舰艇的运动特性。特别是对舰船进行作战性能评估时，必须考虑海洋环境条件造成的武器系统衰减，准确预报舰船在不同海浪环境下摇荡运动幅度，对于评价舰船平台运动是否超过舰载机起降等空中活动限值的时间期望值即平台适用度是必需的。预报未来极短期船舶运动姿态，不仅有利于舰载机的起降，更能提前加以相反的缓冲力或是操纵，进而降低船舶运动中的横摇、纵摇以及艏摇的摇摆幅度，从而提高船舶的耐波性。

发明内容

本发明的目的是这样实现的：

一种基于非线性模型辨识的船舶运动预报装置，包括数据采集模块、输入处理模块、基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块、预报模块、预报误差分析模块和输出模块，其特征在于，数据采集模块对船舶运动的数据进行采集；输入处理模块对输入数据进行滤波处理；基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块是根据环境的变化自适应修正辨识模型的参数；预报模块自适应修正模型参数对船舶运动的预报；预报误差分析模块检验预报精度的各种误差指标；输出模块显示和输出船舶运动的报结果和预报误差结果。

数据采集模块、输入处理模块、基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块、预报模块、预报误差分析模块、输出模块依次串行连接；采集模块获取船舶运动数据后，经过输入处理模块进行数据处理；处理后的数据进入基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块，得到自适应修正辨识模型的参数；自适应修正辨识模型的参数作为预报模块的输入；预报误差分析模块对预报模块的预报结果进行精度要求检验；最后是输出预报结果和误差结果。

一种基于非线性模型辨识的船舶运动预报方法，其特征在于，包括以下步骤：

步骤一数据采集；

需要的数据包括：航速、横荡速度、横摇角速度、艏摇角速度、舵角以及螺旋桨转速的量测数据，各量测数据的采样周期相等且为1秒；

步骤二输入处理；

数据的采样周期为1秒，仿真时间为500秒，对辨识数据进行卡尔曼滤波平滑处理；

适应的系统方程为：

X(k+1)＝ΦX(k)+w(k)

Z(k)＝HX(k)+n(k) (1)

式中，X(k)分别取横荡速度v，横摇角速度p，艏摇角速度r进行分批滤波平滑，或增加状态变量的维度，进行整批的滤波平滑；Z(k)为观测矩阵；w(k)为零均值的白噪声序列，方差为Q；n(k)是零均值的白噪声序列，方差为T为采样时间间隔，T＝1秒；

卡尔曼滤波递推公式：

X(k+1|k)＝ΦX(k)

P(k+1|k)＝ΦP(k)Φ^T+Q(k)

P(k+1)＝[I-K(k+1)H]P(k+1|k)

X(k+1)＝X(k+1|k)+K(k+1)·[Z(k+1)-HX(k+1|k)] (2)

式中：卡尔曼滤波递推计算的初值选取

步骤三基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正；包括以下步骤：

步骤1建立横向运动船舶模型化的动量及动量矩方程；

其中m为质量，u为纵荡速度，为纵荡加速度，v为横荡速度，/>为横荡加速度，r为首摇角速度，/>为首摇角加速度，p为横摇角速度，/>为横摇角加速度，x_G为船舶重心到x轴距离，z_G为船舶重心到z轴距离，I_zz为船体对O_x轴的惯性矩，I_xx为船体对O_z轴的惯性矩，I_xz为关于zx平面的惯性积；将船体坐标系建立在船舶重心处，即x_G和z_G为0；而X，Y，K和N分别代表作用于船体的水动力及力矩，它们是船舶运动变量和控制量的非线性函数；

步骤2由分离型模型的思想建立方程右端作用于船舶上的流体惯性力或力矩、流体黏性力或力矩、螺旋桨推力或力矩、舵力或力矩方程；

船舶在流体中作变速或旋转运动时，将迫使周围的流体介质加速或是减速，流体必然对船体施以反作用力，即流体惯性力；

流体惯性力(力矩)表达式：

式中X_I，Y_I，K_I和N_I分别是纵荡惯性力、横荡惯性矩、横摇惯性力矩和艏摇惯性力矩；J_xx和J_zz为横摇惯性力矩和艏摇惯性力矩系数；由式(3)和(4)得，流体惯性力或力矩作用的结果相当于物体的质量和惯性矩均增加某一值，也称流体惯性力或力矩为附加质量或附加惯性矩；

船舶流体粘性力或力矩；

分别考虑涉及的粘性力或力矩，X_H，Y_H，K_H和N_H分别表示船舶所受的纵荡粘性力、横荡粘性力、横摇粘性力矩和艏摇粘性力矩；

X_H采用如下的估算公式来表示：

X_H＝X_uuu²+X_vvv²+X_vrvr+X_rrr² (5)

式中流体粘性力系数的无量纲估算公式为：

其中S为船体浸湿面积，船体浸湿面积根据船体的型线图，通过微元迭加求得；采用单桨商船的船型进行估算，估算公式如下：

S＝▽^2/3(3.432+0.305L/B+0.443B/d-0.643C_b) (7)

C_t为船舶总阻力系数，船舶的总阻力求取十分困难，需要对其进行分化求取，最后再相加，得到下式：

C_t＝C_f+C_r+ΔC_AR (8)

其中C_f为摩擦阻力系数；C_r为剩余阻力系数；ΔC_AR为粗糙度补粘系数；C_f是由船舶阻力理论中引入的平板概念加上实际船体表面的曲率和粗糙度的影响得来，采用在西欧广泛应用的桑海公式：

C_f＝0.4631/(lgR_n)^2.6 (9)

其中R_n＝VL/ν为雷诺数，ν为流体运动粘性系数；

对于考虑横摇的船舶运动，粘性力或力矩Y_H，N_H不仅与操纵状态有关，还与横摇角有耦合关系，表示为：

船舶形状左右对称，Y_H和N_H是u的偶函数，则Y_H，N_H对u及与u耦合部分的流体动力导数为0，即Y_u＝Y_u...＝N_u＝N_u...＝0，因此Y_H，N_H与u无关；分Y_H，N_H为两部分，即

其中x_c是船舶中心在船体坐标系中的x轴坐标值；

利用井上模型给出Y_H0，N_H0的表达式：

通过大量的回转试验得到式(12)中流体粘性力系数的无量纲形式为：

式中τ′＝(d_A-d_F)/d为无因次化船尾与船艏的吃水差，其中d＝(d_A+d_F)/2，式中d_A，d_F分别为船尾与船艏的吃水深，l_v＝λ/[(π/2)λ+1.4C_bB/L]，λ＝2d/L；

采用如下公式估算：

艏摇粘性力矩的无量纲化的水动力系数估算公式：

船舶运动引起的横摇粘性力矩K_H表示为：

式中是横摇阻尼力矩，/>是横摇恢复力矩，Y_H·z_H是船体流体动力Y_H对x轴的力矩，z_H为Y_H作用点的z坐标；

横摇阻尼力矩

当横摇角较小时，即时，/>是横摇角速度的线性函数，即

螺旋桨计算模型；

以下为螺旋桨推力及转矩计算模型：

式中X_p，Y_p，K_p和N_p分别是螺旋桨的纵荡推力、横荡力、横摇力矩和艏摇力矩，t_p为螺旋桨的推力减额系数，T为螺旋桨推力，一般表示为

T＝ρn²D_p ⁴k_T(J_p) (19)

式中ρ为海水密度，n，D_p，k_T(J_p)分别代表螺旋桨每秒钟的转数，螺旋桨直径和桨的推力系数，其中k_T(J_p)是进速系数J_p的函数，而J_p的表达式形式如下：

J_p＝(1-w_p)u/nD_p (20)

式中w_p为螺旋桨的伴流系数；

舵力及力矩计算模型为：

式中X_R，Y_R和N_R分别表示舵产生的纵荡力，横荡力和艏摇力矩；其中，F_N为垂直于舵叶平面的正压力；t_R为舵力减额系数；a_H为操舵诱导船体横向力的修正因子；x_H为操舵诱导船体横向力作用中心到船舶重心的距离；x_R为作用于舵上的横向力作用点的纵向坐标，δ为舵角；

考虑螺旋桨、船体对舵的干涉时正压力的计算F_N；

在MMG模型中，安装在船、桨后的单独舵，其正压力F_N为

其中：A_R为舵叶面积；U_R为舵处来流有效流速；α_R为舵处来流有效冲角；f_α为舵的升力系数在冲角α＝0时的斜率；

综上所述，分离型船舶运动数学模型的结构为：

步骤3建立规范化的分离型船舶运动数学模型，构成分离型船舶运动模型结构；

根据式(23)，整理推导，降低模型参数辨识的复杂度和模型参数的数量，将模型的参数设为未知参数，并进行规范化，得出下式：

从式(24)中得处，辨识状态空间型的船舶运动模型，需要有航速u，横荡速度v，横摇角速度p，艏摇角速度r，螺旋桨转速n以及舵角δ的时间序列数据，这些时间数据序列是由实船航行时测量仪器测得的量测数据；式中需要辨识的模型参数有a_i(i＝1～6)，b_i(i＝1～7)，c_i(i＝1～7)，d_i(i＝1～6)，共26个未知模型参数；

步骤4离散化船舶运动模型结构；

根据式(24)对船舶运动模型进行简单的离散化，离散后的式子如下：

u(k+1)＝u(k)+T(a₁(k)u(k)²+a₂(k)v(k)²+a₃(k)r(k)²

+a₄(k)v(k)r(k)+a₅(k)n(k)²+a₆(k)sinδ(k)) (25)

v(k+1)＝v(k)+T(b₁(k)v(k)+b₂(k)r(k)+b₃(k)v(k)²+b₄(k)r(k)²

+b₅(k)u(k)r(k)+b₆(k)v(k)r(k)+b₇(k)cosδ(k)) (26)

p(k+1)＝p(k)+T(c₁(k)p(k)+c₂(k)v(k)+c₃(k)r(k)+c₄(k)v(k)²

+c₅(k)r(k)²+c₆(k)v(k)r(k)+c₇(k)cosδ(k)) (27)

r(k+1)＝r(k)+T(d₁(k)r(k)+d₂(k)v(k)+d₃(k)v(k)²+

d₄(k)r(k)²+d₅(k)v(k)r(k)+d₆(k)cosδ(k)) (28)

式中T为步长，又称为采样时间间隔；

步骤5将模型中未知的模型参数作为扩展Kalman滤波算法的状态变量；

扩展后的系统的状态方程和观测方程：

x^a(k+1)＝f^a(x^a(k),s(k),k)+w(k)

y(k)＝Hx^a(k)+e(k) (29)

式中：s(k)表示T(k)和T(k+1)两个时刻输入的采样均值，T为采样时间间隔，w(k)和e(k)是零均值白噪声序列，方差分别为Q和R；

f^a＝[f^a ₁f^a ₂...f^a ₃₀]^T，f^a ₁到f^a ₃₀的表达式：

f₁ ^a＝u(k+1)，f₂ ^a＝v(k+1)，f₃ ^a＝p(k+1)，f₄ ^a＝r(k+1)，f₅ ^a＝a₁(k)，……，f₃₀ ^a＝d₆(k)；

从f^a ₁到f^a ₃₀的表达式得，模型的参数分别扩展为从f^a ₅到f^a ₃₀的26个状态变量，扩展后系统的状态变量为30个；

步骤6对系统的观测矩阵进行增广变换；

对系统状态变量的扩展，相应的系统观测矩阵也应该进行增广变换，

步骤7对非线性船舶运动数学模型进行线性化，以求取系统的状态转移矩阵；

对于非线性系统，状态转移矩阵Φ在每个时刻都是不一样的，以每一个变量对系统的状态矩阵逐一求取偏导：

由此构成扩展Kalman滤波辨识算法：

由式(32)计算出式(27)的状态滤波值f^a，从而得出模型参数a_i、b_i，c_i和d_i的估计值；

步骤8设置扩展Kalman滤波算法的初值；

初始条件取为：

P(0)＝a I，a＝10^5，

I为30×30的单位矩阵，x^a(0)＝[u(0) v(0) p(0) r(0) 0 ... 0]^T；螺旋桨转速和舵角根据实船的量测数据决定；

步骤9按步递推辨识，获得准确的模型参数值；

步骤10应用基于扩展Kalman滤波算法的分离型船舶运动模型的预报；

应用扩展Kalman滤波算法辨识非线性船舶运动模型参数，进而建立完整的运动模型；用辨识的模型参数的最后一次递推值作为运动模型的参数值，以此模型递推计算便是船舶运动的预报，下面给出非线性船舶运动模型的预报方程：

式中N代表算法递推的总次数，这里的N＝240，航速u，横荡速度v，横摇角速度p，艏摇角速度r的初值分别取先验辨识数据最后一时刻的值，递推的步长依然是h＝1秒；

步骤四船舶运动预报；

设绕ζ轴(z轴)的旋转角度为ψ，为艏摇角；绕η轴(y轴)的旋转角度为θ，为纵摇角；绕ξ轴(x轴)的旋转角度为φ，为横摇角；

角速度与角度之间的关系式：

当不考虑纵摇时纵摇的角度θ＝0°，通过式(33)和(34)求出船舶的姿态角；

步骤五预报效果分析；

步骤六预报结果评价。

与现有技术相比，本发明的优势在于：本发明能够根据船舶航行的海况自适应修正模型的参数，适应具有干扰环境下的船舶运动预报，预报效果好。

附图说明

图1为基于模型辨识的船舶横摇角的预报曲线；

图2为几种预报算法的预报数据曲线比较；

图3为几种预报算法的预报数据曲线的局部放大图；

图4为基于模型辨识的船舶运动预报装置的模块组成框图；

图5为基于模型辨识的船舶运动预报方法的流程图。

具体实施方式

下面结合附图对本发明作出以下详细说明：

本发明公开了一种基于非线性模型辨识的船舶运动预报方法及装置，属于船舶运动控制、预报领域。本发明以分离型模型与扩展Kalman滤波算法相结合，应用扩展Kalman滤波算法辨识分离型的船舶运动非线性数学模型参数，并随环境变化自适应修正模型参数，建立船舶运动自适应预报模型，对船舶运动进行在线极短期预报。本发明的装置包括依次串行连接的数据采集模块、输入处理模块、基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块、预报模块、预报结果误差分析模块、输出模块。本发明的方法分析建立分离型船舶运动非线性模型、推导变换为包含模型参数的扩展状态方程、采用线性化方法求取系统方程的状态转移矩阵、建立扩展卡尔曼滤波器估计方程，实现在线自适应模型参数辨识算法、建立了在线自适应预报模型，该预报模型基于船舶运动非线性模型，使预报精度更进一步提高，预报时长有所增加，与以往预报方法比较物理意义明确，稳定性提高。

一种基于模型辨识的船舶运动预报方法及装置，其特征在于包括：以分离型模型的思想建立船舶运动非线性数学模型。模型是考虑横摇在内的横向船舶运动模型，即包括航速、横荡速度、横摇角速度以及艏摇角速度的四自由度船舶运动非线性模型。利用分离型模型的船舶运动模型结构，经过离散化变换，应用扩展Kalman滤波算法自适应辨识分离型模型的模型参数，进而建立完整的船舶运动非线性预报模型。需要的先验辨识数据有：航速、横荡速度、横摇角速度、艏摇角速度、舵角以及螺旋桨转速的量测数据，各量测数据的采样周期要相等。为了进行成功的递推参数辨识，采样的数据点数要足够的多。

该方法的具体步骤如下：

(1)数据采集；

(2)输入处理；

(3)基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正，下面是对这一步骤的详细分解：

3.1建立横向运动船舶模型化的动量及动量矩方程；

3.2由分离型模型的思想分析建立模型右端作用于船舶上的流体惯性力(力矩)、流体黏性力(力矩)、螺旋桨推力(力矩)、舵力(力矩)方程；

3.3建立规范化的分离型船舶运动数学模型，将模型的参数(水动力系数等)设为未知参数，构成分离型船舶运动模型结构；

3.4离散化船舶运动模型结构；

3.5将模型中未知的模型参数作为扩展Kalman滤波算法的状态变量；

3.6对系统的观测矩阵进行增广变换；

3.7对非线性船舶运动数学模型进行线性化处理，以求取系统的状态转移矩阵；

3.8设置扩展Kalman滤波算法的初值；

3.9按步递推自适应辨识，获得模型参数辨识值；

3.10基于扩展Kalman滤波算法辨识的分离型船舶运动模型，建立预报模型和预报算法。

(4)船舶运动预报；

(5)预报效果分析；

(6)预报结果评价；

(7)预报结果输出。

实现本发明的技术方案：

一、一种基于模型辨识的船舶运动预报方法，其特征在于包括以下步骤：

1、数据采集

需要的数据有：航速、横荡速度、横摇角速度、艏摇角速度、舵角以及螺旋桨转速的量测数据，各量测数据的采样周期要相等，采样周期可根据期望预报时间长度和预报数据的时间间隔选定，本发明实验用的数据采样周期是1秒。为了进行成功的递推参数辨识，采样的数据点数要足够的多。

2、输入处理

为了模型参数辨识的可靠性，需要对辨识所需数据用卡尔曼滤波算法进行数据的平滑预处理。

数据的采样周期为1秒，为了参数辨识能有足够的时间数据序列，尽量仿真时间长度放的长一些，这里用500秒，基本足够扩展卡尔曼滤波辨识算法递推辨识参数和后续对船舶运动预报的研究。为了扩展卡尔曼辨识算法能有更好的辨识效果，需要对辨识数据进行卡尔曼滤波平滑处理，这里需要写出适应的系统方程，方程如下：

X(k+1)＝ΦX(k)+w(k)

Z(k)＝HX(k)+n(k) (1)

式中，X(k)可分别取横荡速度v，横摇角速度p，艏摇角速度r等进行分批滤波平滑，也可增加状态变量的维度，进行整批的滤波平滑，Z(k)为观测矩阵，w(k)为零均值的白噪声序列，方差为Q，n(k)是零均值的白噪声序列，方差为T为采样时间间隔，T＝1秒，卡尔曼滤波递推公式：

X(k+1|k)＝ΦX(k)

P(k+1|k)＝ΦP(k)Φ^T+Q(k)

P(k+1)＝[I-K(k+1)H]P(k+1|k)

X(k+1)＝X(k+1|k)+K(k+1)·[Z(k+1)-HX(k+1|k)] (2)

式中：卡尔曼滤波递推计算的初值选取

3、基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正

3.1建立横向运动船舶模型化的动量及动量矩方程

其中m为质量，u为纵荡速度，为纵荡加速度，v为横荡速度，/>为横荡加速度，r为首摇角速度，/>为首摇角加速度，p为横摇角速度，/>为横摇角加速度，x_G为船舶重心到x轴距离，z_G为船舶重心到z轴距离，I_zz为船体对O_x轴的惯性矩，I_xx为船体对O_z轴的惯性矩，I_xz为关于zx平面的惯性积。本发明将船体坐标系建立在船舶重心处，即x_G和z_G为0。而X，Y，K和N分别代表作用于船体的水动力及力矩，它们是船舶运动变量和控制量的非线性函数。

3.2由分离型模型的思想建立方程右端作用于船舶上的流体惯性力(力矩)、流体黏性力(力矩)、螺旋桨推力(力矩)、舵力(力矩)方程

(1)船舶在流体中作变速或旋转运动时，它将迫使周围的流体介质加速或是减速，流体必然对船体施以反作用力，此即流体惯性力。

流体惯性力(力矩)表达式：

式中X_I，Y_I，K_I和N_I分别是纵荡惯性力、横荡惯性矩、横摇惯性力矩和艏摇惯性力矩。J_xx和J_zz为横摇惯性力矩和艏摇惯性力矩系数。由式(3)和(4)可知，流体惯性力(力矩)作用的结果相当于物体的质量和惯性矩均增加某一值，所以也称流体惯性力(力矩)为附加质量或附加惯性矩。

(2)船舶流体粘性力(力矩)

下面分别考虑研究涉及到的粘性力(力矩)。X_H，Y_H，K_H和N_H分别表示船舶所受的纵荡粘性力、横荡粘性力、横摇粘性力矩和艏摇粘性力矩。(第一次出现解释)

1)X_H的研究

X_H可以采用如下的估算公式来表示：

X_H＝X_uuu²+X_vvv²+X_vrvr+X_rrr² (5)

式中流体粘性力系数的无量纲估算公式为：

其中S为船体浸湿面积，船体浸湿面积可根据船体的型线图，通过微元迭加求得。但由于船体的型线图很难获得，亦可采用单桨商船的船型进行估算，估算公式如下：

S＝▽^2/3(3.432+0.305L/B+0.443B/d-0.643C_b) (7)

C_t为船舶总阻力系数，船舶的总阻力求取十分困难，需要对其进行分化求取，最后再相加，所以有下式：

C_t＝C_f+C_r+ΔC_AR (8)

其中C_f为摩擦阻力系数；C_r为剩余阻力系数；ΔC_AR为粗糙度补粘系数。C_f是由船舶阻力理论中引入的平板概念加上实际船体表面的曲率和粗糙度的影响得来，估算公式有多种，这里我们采用在西欧广泛应用的桑海公式：

C_f＝0.4631/(lgR_n)^2.6 9)

其中R_n＝VL/ν为雷诺数，ν为流体运动粘性系数。

2)Y_H，N_H的研究

对于考虑横摇的船舶运动，粘性力(力矩)Y_H，N_H不仅与操纵状态有关，还与横摇角有耦合关系，可表示为

由于船舶形状左右对称，Y_H和N_H是u的偶函数，则Y_H，N_H对u及与u耦合部分的流体动力导数为0，即Y_u＝Y_u...＝N_u＝N_u...＝0，因此可以说Y_H，N_H与u无关。为研究方便，分Y_H，N_H为两部分，即

其中x_c是船舶中心在船体坐标系中的x轴坐标值。

①利用井上模型给出Y_H0，N_H0的表达式：

通过大量的回转试验得到式(12)中流体粘性力系数的无量纲形式为：

式中τ′＝(d_A-d_F)/d为无因次化船尾与船艏的吃水差，其中d＝(d_A+d_F)/2，式中d_A，d_F分别为船尾与船艏的吃水深，l_v＝λ/[(π/2)λ+1.4C_bB/L]，λ＝2d/L。

②采用如下公式估算：

艏摇粘性力矩的无量纲化的水动力系数估算公式：

3)K_H的研究

船舶运动引起的横摇粘性力矩K_H可表示为

式中是横摇阻尼力矩，/>是横摇恢复力矩，Y_H·z_H是船体流体动力Y_H对x轴的力矩，z_H为Y_H作用点的z坐标。下面分别介绍：

①横摇阻尼力矩

当横摇角较小时，即时，/>是横摇角速度的线性函数，即/>

(3)螺旋桨计算模型

以下为螺旋桨推力及转矩计算模型：

本发明讨论的工况是船舶正向前进，即螺旋桨工作在第一象限，此时转速n＞0，u＞0，综合考虑到伴流和推力减额的螺旋桨推力和转矩模型。

式中X_p，Y_p，K_p和N_p分别是螺旋桨的纵荡推力、横荡力、横摇力矩和艏摇力矩，t_p为螺旋桨的推力减额系数，T为螺旋桨推力，一般可表示为

T＝ρn²D_p ⁴k_T(J_p) (19)

式中ρ为海水密度，n，D_p，k_T(J_p)分别代表螺旋桨每秒钟的转数，螺旋桨直径和桨的推力系数，其中k_T(J_p)是进速系数J_p的函数，而J_p的表达式形式如下：

J_p＝(1-w_p)u/nD_p (20)

式中w_p为螺旋桨的伴流系数。

(4)舵力及力矩计算模型

舵力及力矩计算模型为：

式中X_R，Y_R和N_R分别表示舵产生的纵荡力，横荡力和艏摇力矩；其中，F_N为垂直于舵叶平面的正压力；t_R为舵力减额系数；a_H为操舵诱导船体横向力的修正因子；x_H为操舵诱导船体横向力作用中心到船舶重心的距离；x_R为作用于舵上的横向力作用点的纵向坐标，δ为舵角。

1)考虑螺旋桨、船体对舵的干涉时正压力的计算F_N

在MMG模型中，安装在船、桨后的单独舵，其正压力F_N为

其中：A_R为舵叶面积；U_R为舵处来流有效流速；α_R为舵处来流有效冲角；f_α为舵的升力系数在冲角α＝0时的斜率。

综上所述，分离型船舶运动数学模型的结构为：

3.3建立规范化的分离型船舶运动数学模型，构成分离型船舶运动模型结构

根据式(23)，整理推导，为了降低模型参数辨识的复杂度和模型参数的数量，这里忽略模型中三阶以及三阶以上的项，将模型的参数(水动力系数等)设为未知参数，并进行规范化，得出下式：

从式(24)中我们知道，为了辨识状态空间型的船舶运动模型，需要有航速u，横荡速度v，横摇角速度p，艏摇角速度r，螺旋桨转速n以及舵角δ的时间序列数据，这些时间数据序列是由实船航行时测量仪器测得的量测数据。式中需要辨识的模型参数有a_i(i＝1～6)，b_i(i＝1～7)，c_i(i＝1～7)，d_i(i＝1～6)，共26个未知模型参数。

3.4离散化船舶运动模型结构

根据式(24)对船舶运动模型进行简单的离散化，离散后的式子如下：

u(k+1)＝u(k)+T(a₁(k)u(k)²+a₂(k)v(k)²+a₃(k)r(k)²

+a₄(k)v(k)r(k)+a₅(k)n(k)²+a₆(k)sinδ(k)) (25)

v(k+1)＝v(k)+T(b₁(k)v(k)+b₂(k)r(k)+b₃(k)v(k)²+b₄(k)r(k)²

+b₅(k)u(k)r(k)+b₆(k)v(k)r(k)+b₇(k)cosδ(k)) (26)

p(k+1)＝p(k)+T(c₁(k)p(k)+c₂(k)v(k)+c₃(k)r(k)+c₄(k)v(k)²

+c₅(k)r(k)²+c₆(k)v(k)r(k)+c₇(k)cosδ(k)) (27)

r(k+1)＝r(k)+T(d₁(k)r(k)+d₂(k)v(k)+d₃(k)v(k)²+

d₄(k)r(k)²+d₅(k)v(k)r(k)+d₆(k)cosδ(k)) (28)

式中T为步长，又称为采样时间间隔。

3.5将模型中未知的模型参数作为扩展Kalman滤波算法的状态变量

下面给出扩展后的系统的状态方程和观测方程：

x^a(k+1)＝f^a(x^a(k),s(k),k)+w(k)

y(k)＝Hx^a(k)+e(k) (29)

式中：s(k)表示T(k)和T(k+1)两个时刻输入的采样均值，T为采样时间间隔，w(k)和e(k)是零均值白噪声序列，方差分别为Q和R，f^a＝[f^a ₁ f^a ₂ ... f^a ₃₀]^T，下面给出f^a ₁到f^a ₃₀的表达式：

f₁ ^a＝u(k+1)，f₂ ^a＝v(k+1)，f₃ ^a＝p(k+1)，f₄ ^a＝r(k+1)，f₅ ^a＝a₁(k)，……，f₃₀ ^a＝d₆(k)。

从f^a ₁到f^a ₃₀的表达式可以看出，模型的参数分别扩展为从f^a ₅到f^a ₃₀的26个状态变量，扩展后系统的状态变量为30个。

3.6对系统的观测矩阵进行增广变换

由于对系统状态变量的扩展，相应的系统观测矩阵也应该进行增广变换，

3.7对非线性船舶运动数学模型进行线性化，以求取系统的状态转移矩阵

对于非线性系统，状态转移矩阵Φ在每个时刻都是不一样的，对状态转移矩阵的求取我们以每一个变量对系统的状态矩阵逐一求取偏导：

由此构成扩展Kalman滤波辨识算法：

由式(32)计算出式(27)的状态滤波值f^a，从而得出模型参数a_i、b_i，c_i和d_i的估计值。

3.8设置扩展Kalman滤波算法的初值

初始条件取为：P(0)＝a I，a＝10^5，I为30×30的单位矩阵，x^a(0)＝[u(0) v(0)p(0) r(0) 0 ... 0]^T。螺旋桨转速和舵角根据实船的量测数据决定。

3.9按步递推辨识，获得准确的模型参数值

经过船舶运动数学模型的变换，应用扩展Kalman滤波算法对船舶运动模型中的参数进行辨识，螺旋桨转速n一直定为1000转每分钟。扩展卡尔曼滤波辨识算法递推250次已基本看出辨识参数的收敛性和可靠性。参数辨识试验中，经过250次递推计算后，部分模型参数能不受海浪的影响，或是受海浪干扰影响较小的模型参数能较好的收敛到某一值附近，已经基本可达到真实可用的参数值，

带有海浪干扰的参数辨识，仿真时加入了三级海况下的海浪干扰，通过模型参数a₁、a₅、b₁、b₃、c₂、c₄，d₁以及d₃的辨识曲线可以看出，这些参数在扩展卡尔曼辨识算法递100次以后甚至有的参数在递推50次以后就已经收敛在某个值的附近，一是说明辨识算法的可行性和可靠性，二是说明这类参数受环境干扰小，这些参数多是与船舶自身的参数有关(水线长，吃水深，横稳心高度等)，在参数辨识中易于收敛。

通过模型参数a₄、b₅，b₆和c₅的递推辨识曲线来看，辨识曲线尽管递推了250次也依然没有收敛的迹象，恰恰说明这类模型参数与外界干扰(海浪干扰)有关，由于辨识数据由MMG非线性模型仿真得来，海浪干扰是随机的，不确定的。本发明提出的模型参数辨识可以随外部环境的变化(海况的不同)而自适应修正模型中的参数，建立船舶自适应运动非线性预报模型，提高预报精度。通过对比辨识建模的仿真数据和先验辨识数据，说明扩展Kalman滤波算法对非线性模型参数的辨识是可行的，可靠的。

3.10应用基于扩展Kalman滤波算法的分离型船舶运动模型的预报

应用扩展Kalman滤波算法辨识非线性船舶运动模型参数，进而建立完整的运动模型。用辨识的模型参数的最后一次递推值作为运动模型的参数值，以此模型递推计算便是船舶运动的预报，下面给出非线性船舶运动模型的预报方程：

式中N代表算法递推的总次数，这里的N＝240，航速u，横荡速度v，横摇角速度p，艏摇角速度r的初值分别取先验辨识数据最后一时刻的值，递推的步长依然是h＝1秒。

4、船舶运动预报

为了研究方便，我们设绕ζ轴(z轴)的旋转角度为ψ，称为艏摇角；绕η轴(y轴)的旋转角度为θ，称为纵摇角；绕ξ轴(x轴)的旋转角度为φ，称为横摇角。下面给出角速度与角度之间的关系式：

当不考虑纵摇时纵摇的角度θ＝0°，通过式(33)和(34)即可求出船舶的姿态角。

5、预报效果分析

下面给出基于模型辨识的船舶横摇角度的预报曲线图：

从图1中可看出预测值与真实值的拟合度很高，预测效果很好。为了进一步说明基于模型辨识的预报精度的优势。

6、预报结果评价

为了进一步说明基于扩展Kalman滤波辨识算法的预报精度的优势，下面与基于最小二乘算法的AR模型和基于Kalman滤波算法的AR模型预报进行比较，并通过图形和数据共同验证。

下面给出基于扩展Kalman滤波辨识算法、基于最小二乘算法的AR模型以及基于Kalman滤波算法的AR模型预报未来20秒的船舶横摇角速度的运动姿态数据曲线图：

从图2和图3中的三种预报方法的预报数据波形图可以看出，扩展Kalman滤波辨识的预报算法除了预报数据与真实数据的拟合度更好以外，相位差也较基于最小二乘算法的AR模型和基于Kalman滤波算法的AR模型的预报算法要小。

下表1展现了这三种算法预报未来20秒的数据与真实数据之间的相对误差和均方根误差：

表1三种算法预报未来20秒的相对误差和均方根误差

通过以上曲线图和误差数据可看出基于扩展Kalman的模型参数辨识预报的可靠性。

二、一种基于非线性模型辨识的船舶运动预报装置，不仅在预报方法上提出了新方法，在装置组成、结构上也提出了新的思路。本发明装置采用了模块化结构，模块化结构利于装置升级和维护；同时在线实时采集实船数据、在线建模、在线预报，是实时在线预报装置；预报效果仿真分析模块、预测结果评价模块，使应用者实时掌握预报误差，做出正确的判断、和决策。

本装置有硬件系统和软件系统组成：硬件由数据采集卡、计算机系统、数据输出接口组成；软件由各模块组成：包括数据采集模块、输入处理模块、基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块、预报模块、预报误差分析模块、输出模块，这些模块可用C/C++语言，MATLAB等语言实现。所述数据采集模块用于船舶运动的数据进行采集；所述输入处理模块可以用于对输入数据进行滤波处理；所述基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块是根据环境的变化自适应修正辨识模型的参数；所述预报模块自适应修正模型参数对船舶运动的预报；所述预报结果评价模块用于检验预报精度的各种误差指标；所述输出模块用于显示和输出船舶运动的报结果和预报误差结果。

各模块之间的信号流程为：数据采集模块、输入处理模块、基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块、预报模块、预报误差分析模块、输出模块依次串行连接。各模块必须按顺序执行，上一模块的输出是下一模块的输入。船舶运动数据采集模块后，经过输入处理模块进行数据处理；处理后的可用数据进入基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块；自适应修正辨识模型的参数作为预报模块的输入；利用预报效果仿真分析模块及结果评价模块对预报结果进行精度要求检验，如果符合要求则预报数据可用于决策、船舶的操纵、运行等操作，最后是输出预报结果和误差结果。装置可不断更新采集的数据，进行实时的在线预报和输出。

本发明的组成框图和流程图如图4、图5所示，由数据采集模块、输入处理模块、基于扩展卡尔曼滤波的在线模型参数辨识及自适应校正模块、预报模块、预报效果仿真分析模块、预报结果评价模块、输出模块。

Claims (1)

1.一种基于非线性模型辨识的船舶运动预报方法，其特征在于，包括以下步骤：

步骤1：数据采集，包括纵荡速度u、横荡速度v、横摇角速度p、艏摇角速度r、舵角δ以及螺旋桨转速n的量测数据，各量测数据的采样周期相等；

步骤2：对采集的数据进行卡尔曼滤波平滑处理；

步骤3：建立横向运动船舶模型化的动量及动量矩方程；

其中，m为质量；为纵荡加速度；/>为横荡加速度；/>为艏摇角加速度；/>为横摇角加速度；x_G为船舶重心到x轴距离，z_G为船舶重心到z轴距离，由于船体坐标系建立在船舶重心处，x_G和z_G为0；I_zz为船体对O_x轴的惯性矩；I_xx为船体对O_z轴的惯性矩；I_xz为关于zx平面的惯性积；X、Y、K和N分别代表作用于船体的水动力及力矩，它们是船舶运动变量和控制量的非线性函数；

步骤4：由分离型模型的思想建立方程右端作用于船舶上的流体惯性力和力矩、流体黏性力和力矩、螺旋桨推力、舵力和力矩方程；

步骤4.1：船舶在流体中作变速或旋转运动时，将迫使周围的流体介质加速或是减速，流体必然对船体施以反作用力，即流体惯性力；流体惯性力和力矩表达式：

式中，X_I、Y_I、K_I和N_I分别是纵荡惯性力矩、横荡惯性力矩、横摇惯性力矩和艏摇惯性力矩；J_xx和J_zz为横摇惯性力矩和艏摇惯性力矩系数；

步骤4.2：船舶所受的纵荡粘性力X_H的表达式：

X_H＝X_uuu²+X_vvv²+X_vrvr+X_rrr²

式中流体粘性力系数的无量纲估算公式为：

其中，S为船体浸湿面积，C_t为船舶总阻力系数，C_t＝C_f+C_r+ΔC_AR；C_f为摩擦阻力系数，C_f＝0.4631/(lgR_n)^2.6，R_n为雷诺数，v为流体运动粘性系数；C_r为剩余阻力系数，ΔC_AR为粗糙度补粘系数；

步骤4.3：横荡粘性力Y_H、艏摇粘性力矩N_H的表达式：

其中，式中流体粘性力系数的无量纲形式为：/>τ′＝(d_A-d_F)/d为无因次化船尾与船艏的吃水差；d＝(d_A+d_F)/2，d_A与d_F分别为船尾与船艏的吃水深；l_v＝λ/[(π/2)λ+1.4C_bB/L]；λ＝2d/L；/>式中艏摇粘性力矩的无量纲化的水动力系数为：/>x_c是船舶中心在船体坐标系中的x轴坐标值；

步骤4.4：船舶运动引起的横摇粘性力矩K_H表示为：

式中，是横摇阻尼力矩；/>是横摇恢复力矩；Y_H·z_H是船体流体动力Y_H对x轴的力矩，z_H为Y_H作用点的z坐标；

步骤4.5：螺旋桨推力T的表达式为：

T＝ρn²D_p ⁴k_T(J_p)

其中，ρ为海水密度；n和D_p分别代表螺旋桨每秒钟的转数和螺旋桨直径；k_T(J_p)表示螺旋桨的推力系数，k_T(J_p)是进速系数J_p的函数，J_p＝(1-w_p)u/nD_p，w_p为螺旋桨的伴流系数；

步骤4.6：舵力及力矩的表达式为：

式中，X_R、Y_R和N_R分别表示舵产生的纵荡力、横荡力和艏摇力矩；F_N为垂直于舵叶平面的正压力，A_R为舵叶面积，U_R为舵处来流有效流速，α_R为舵处来流有效冲角，f_α为舵的升力系数在冲角α＝0时的斜率；t_R为舵力减额系数；a_H为操舵诱导船体横向力的修正因子；x_H为操舵诱导船体横向力作用中心到船舶重心的距离；x_R为作用于舵上的横向力作用点的纵向坐标；

步骤5：建立规范化的分离型船舶运动数学模型，构成分离型船舶运动模型结构；

式中共26个未知模型参数a_i、b_i、c_i、d_i；

离散化船舶运动模型结构：

u(k+1)＝u(k)+T(a₁(k)u(k)²+a₂(k)v(k)²+a₃(k)r(k)²+a₄(k)v(k)r(k)+a₅(k)n(k)²+a₆(k)sinδ(k))

v(k+1)＝v(k)+T(b₁(k)v(k)+b₂(k)r(k)+b₃(k)v(k)²+b₄(k)r(k)²+b₅(k)u(k)r(k)+b₆(k)v(k)r(k)+b₇(k)cosδ(k))

p(k+1)＝p(k)+T(c₁(k)p(k)+c₂(k)v(k)+c₃(k)r(k)+c₄(k)v(k)²+c₅(k)r(k)²+c₆(k)v(k)r(k)+c₇(k)cosδ(k))

r(k+1)＝r(k)+T(d₁(k)r(k)+d₂(k)v(k)+d₃(k)v(k)²+d₄(k)r(k)²+d₅(k)v(k)r(k)+d₆(k)cosδ(k))

式中，T为步长，即采样时间间隔；

步骤6：将模型中未知的模型参数作为扩展Kalman滤波算法的状态变量；

扩展后的系统的状态方程和观测方程：

x^a(k+1)＝f^a(x^a(k),s(k),k)+w(k)

y(k)＝Hx^a(k)+e(k)

式中：s(k)表示T(k)和T(k+1)两个时刻输入的采样均值；w(k)和e(k)是零均值白噪声序列，方差分别为Q和R；f^a＝[f^a ₁ f^a ₂ ... f^a ₃₀]^T，f₁ ^a＝u(k+1)，f₂ ^a＝v(k+1)，f₃ ^a＝p(k+1)，f₄ ^a＝r(k+1)，f₅ ^a＝a₁(k)，……，f₃₀ ^a＝d₆(k)；模型的参数分别扩展为从f^a ₅到f^a ₃₀的26个状态变量，扩展后系统的状态变量为30个；

步骤7：对非线性船舶运动数学模型进行线性化，以求取系统的状态转移矩阵；

对于非线性系统，状态转移矩阵Φ在每个时刻都是不一样的，以每一个变量对系统的状态矩阵逐一求取偏导：

由此构成扩展Kalman滤波辨识算法：

设置扩展Kalman滤波算法的初值，计算出状态滤波值f^a，从而得出模型参数a_i、b_i，c_i和d_i的估计值；

步骤8：应用扩展Kalman滤波算法辨识非线性船舶运动模型参数，进而建立完整的运动模型；用辨识的模型参数的最后一次递推值作为运动模型的参数值，以此模型递推计算便是船舶运动的预报，非线性船舶运动模型的预报方程为：

式中，N代表算法递推的总次数；纵荡速度u、横荡速度v、横摇角速度p、艏摇角速度r的初值分别取先验辨识数据最后一时刻的值。