第6章-船舶运动控制系统建模应用

第6章 船舶运动控制系统建模应用
6.1 引 言
数学模型化(mathematical modelling)是用数学语言(微分方程式)描述实际过程动态特性的方法。

在船舶运动控制领域，建立船舶运动数学模型大体上有两个目的：一个目的是建立船舶操纵模拟器(ship manoeuvring simulator)，为研究闭环系统性能提供一个基本的仿真平台；另一个目的是直接为设计船舶运动控制器服务。

船舶运动数学模型主要可分为非线性数学模型和线性数学模型，前者用于船舶操纵模拟器设计和神经网络控制器、模糊控制器等非线性控制器的训练和优化，后者则用于简化的闭环性能仿真研究和线性控制器(PID, LQ, LQG, H ∞鲁棒控制器)的设计。

船舶的实际运动异常复杂，在一般情况下具有6个自由度。

在附体坐标系内考察，这种运动包括跟随3个附体坐标轴的移动及围绕3个附体坐标轴的转动，前者以前进速度(surge velocity)u 、横漂速度(sway velocity)v 、起伏速度(heave velocity)w 表述，后者以艏摇角速度(yaw rate)r 、横摇角速度(rolling rate)p 及纵摇角速度(pitching rate)q 表述；在惯性坐标系内考察，船舶运动可以用它的3个空间位置000,,z y x (或3个空间运动速度
000,,z y x &&&)和3个姿态角即方位角(heading angle)ψ、横倾角(rolling angle)ϕ、纵倾角
(pitching angle)θ (或3个角速度θϕψ&&&,,)来描述，),,(θϕψ称为欧拉角[4](见图6.1.1)。

显然T ],,[w v u 和T 000],,[z y x &&&以及T
],,[r q p 和T ],,[θϕ
ψ&&&之间有确定关系[4]。

但这并不等于说，我们要把这6个自由度上的运动全部加以考虑。

数学模型是实际系统的简化，如何简化就有很大学问。

太复杂和精细的模型可能包含难于估计的参数，也不便于分析。

过于简单的模型不能描述系统的重要性能。

这就需要我们建模时在复杂和简单之间做合理的折中。

对于船舶运动控制来说，建立一个复杂程度适宜、精度满足研究要求的数学模型是至关重要的。

图6.1.1的坐标定义如下：000Z Y X O -是惯性坐标系(大地参考坐标系)，
为起始
位置，0OX 指向正北，0OY 指向正东，0OZ 指向地心；o -xyz 是附体坐标系，为船首尾之间连线的中点，ox 沿船中线指向船首，oy 指向右舷，oz 指向地心；航向角ψ以正
北为零度，沿顺时针方向取0︒～360︒；舵角δ以右舵为正。

对于大多数船舶运动及其控制问题而言，可以忽略起伏运动、纵摇运动及横摇运动，而只需讨论前进运动、横漂运动和艏摇运动，这样就简化成一种只有3个自由度的平面运动问题。

图6.1.2给出图6.1.1经简化后的船舶平面运动变量描述。

船舶平面运动模型对于像航向保持、航迹跟踪、动力定位、自动避碰等问题，具有足够的精度；但在研究像舵阻摇、大舵角操纵等问题时，则必须考虑横摇运动。

本章根据刚体动力学基本理论建立船舶平面运动基本方程，据此进一步导出状态空间型(线性和非线性)及传递函数型船舶运动数学模型，并考虑了操舵伺服系统的动态特性和风、浪、流干扰的处理方法。

这些结果将作为设计各种船舶运动控制器的基础。

计及横摇的四自由度船舶运动数学模型参见文献[5]。

ψ
图6.1.1 在惯性坐标系和附体坐标系中描述船舶的运动
Y 0
图6.1.2 船舶平面运动变量描述
6.2 船舶平面运动的运动学
(1)坐标系及运动学变量
1)惯性坐标系及与之相关的速度分量 取00Y X O -为固定于地球的大地坐标系，原点
O 设为船舶运动始点或任取，地球的曲率在此可不考虑，不过在涉及大范围航行的航线设
计问题时，需单独处理。

设船舶运动速度向量V 在0OX 方向上的分量为0u ，V 在0OY 方向上的分量为0v ，船舶当前的位置是),(00y x ，时间变量以t 表示，显有
⎪⎭
⎪
⎬⎫=-=-⎰⎰t
t
t v y t y t u x t x 00000000d )0()(d )0()( (6-2-1) 设船舶的艏摇角速度r 顺时针方向为正，有
⎰=-t
t r t 0d )0()(ψψ (6-2-2)
2)附体坐标系及与之相关的速度分量 取附体坐标系oxy 位于满载水线面内。

船舶运动速度V 在ox 方向上的分量为u ，称为前进速度，V 在oy 方向上的分量为v ，叫做横漂速度。

同一速度向量V 在惯性坐标系的分量),(00v u 及附体坐标系的分量),(v u 有下列明显的关
⎥⎦
⎤
⎢⎣⎡⎥⎦⎤⎢
⎣⎡-=⎥⎦⎤⎢⎣⎡v u v u ψψψψ
cos sin sin cos 00 (6-2-3)
3)两种坐标系内运动学变量之间的关系 在惯性坐标系内船舶的位置和姿态由
T 00)](),(),([t t y t x ψ确定，在附体坐标系内船舶之运动速度和角速度由[]T
)(),(),(t r t v t u 表示。

由式(6-2-1)，式(6-2-2)和式(6-2-3)知
⎪⎭
⎪⎬⎫++=-+=+=⎰⎰⎰t
t
t
t t t v t t u y t y t t t v t t u x t x t
t r t 0000000d )](cos )()(sin )([)0()(d )](sin )()(cos )([)0()(d )()0()(ψψψψψψ
(6-2-4)
可见，要确定船舶在任意时刻的位置和姿态，首先应该求出在附体坐标系内u,v,r 的变化规律，为此需要建立船舶运动的动力学方程。

(2)平面运动中船舶各点上速度之间的关系
1)刚体运动分解为移动和转动 从运动控制角度将船舶视为刚体是足够准确的，因此其运动是由移动(translation)和转动(rotation)叠加而成；可以取船上任意一点为参考点，船舶一方面整体地随该参考点平行移动，另一方面绕该参考点同时发生旋转运动；移动速度即参考点的速度，故与参考点选择有关，转动角速度则与参考点无关，即对任意的参考点均为同值，对于船舶平面运动，该转动角速度即为艏摇角速率r 。

2)船舶任意点P 处的合速度 取o 为参考点(图6.2.1)，船上任一点P 对o 点向径为j i j i ρ,,y x o +=为ox 及oy 轴上的单位向量。

以向量形式表示旋转角速度，有k ωr =,k 为沿oz 轴的单位向量，ω即为艏摇角速度向量。

由理论力学，因刚体转动而造成的速度为o r ρωV ⨯=，故P 点的合速度是
j
i ρωV V V V )()(xr v yr u o r P ++-=⨯+=+=
(6-2-5)
注意：单位向量×乘所得向量满足右手法则，如i k ⨯，右手从k 的正方向逆时针握向i 的正方向，大拇指所指方向即j 的正方向，如果方向与j 的正方向相反，结果加负号。

图6.2.1 移动与转动速度的合成
考虑船舶质心C ，其对o 点之向径为j i ρC C C y x +=，则C 点之速度为
P V
j
i j i ρωV V )()()(r x v u r x v r y u C C C C C ++=++-=⨯+=
(6-2-6)
上式最后一步是由于船舶配载对称于纵舯剖面，0=C y 。

如果取质心C 为参考点，应该从oxy 坐标系过渡到ξηC 坐标系，后者是前者沿ox 方向平行移动距离C ρ而得。

P 对C 的向径为j i d ηξ+=，于是有
j
i d ωV V )()(r v r u C C C P ξη++-=⨯+=
(6-2-7)
6.3 船舶平面运动的动力学
在推导船舶运动方程时，做下列假设：≠ 船舶是一个刚体；≡ 大地参照系是惯性参照系；≈ 水动力与频率无关，水的自由表面做刚性壁处理。

有了第一个假设就不用考虑每个质量元素之间的相互作用力的影响，而第二假设则可以消除由于地球相对于恒星参照系的运动所产生的力。

(1)平移运动方程的建立
1)刚体的动量 刚体被看做无数质量微团的集合体，各微团保持其形状及彼此之间的距离不变。

刚体动量G 为各微团动量m P d V 的积分，即
⎰⎰⎰⎰⨯+=⨯+==m m m m C C P d d d )(d d ωV d ωV V G
上式最后一项按照质心的定义应为零，设m 是刚体的总质量，则
C m V G = (6-3-1) 2)刚体动量定理 牛顿运动定律指明，刚体动量的变化率等于其所受外力之和。

以j i F Y X +=代表合外力，其中，X 是作用于ox 方向上的外力，Y 是作用于oy 方向上的外力，有
F G =t d /d (6-3-2) 利用式(6-2-6)、式(6-3-1)和式(6-3-2)，且注意到i j j i r t r t -==d /d ,d /d (因整个坐标系是建立在附体坐标系基础上的，而附体坐标系是随着船舶的移动和转动而移动和转动的，故其导数存在。

如果在惯性坐标系，则其导数为0)，参见图6.3.1，经整理得
⎭
⎬⎫
=++=--Y r x ur v m X r x vr u m C C )()(2&&& (6-3-3)
i
j i j ⋅-=t r d d
图6.3.1 单位向量微分关系
式(6-3-3)即为船舶平移的动力学基本方程，注意其形状与熟知的牛顿方程有所差异，这是由于建立船舶运动数学模型应用的oxy 是非惯性坐标系所致。

式(6-3-3)左端附加项mvr
-及mur 是船舶宏观旋转中向心惯性力分量；附加项2
r mx C -及r mx C &分别是由于质心C 对原点o 做旋转运动产生的向心惯性力及切向惯性力(离心惯性力)。

(2)旋转运动方程的建立
1)刚体的动量矩 刚体对质心C 的动量矩C H 为各微团对C 动量矩)d (m P V d ⨯的积分，即
⎰⎰⎰⎰⎰=+=+⨯⨯+=
⨯⨯+⨯=⨯=k
k j i k j i d ωd V d V d H r I m r m r m m m C P C ζζηξηξηξ]d )([d )()(d )(d )(d )(2
2
(6-3-4)
其中⎰+=m I d )(2
2ηξζζ为船舶对过C 点的垂直轴)(ςo 的惯性矩。

2)对质心C 的动量矩定理 同样由牛顿运动定律，运动着的刚体对质心C 的动量矩变化率等于其所受外力矩之和，以k M C C N =表示后者，C N 为外力矩之代数和，于是
C C t M H =d /d
即 C N r
I =&ζζ (6-3-5) 3)对于坐标系oxy 原点o 的动量矩定理 形式为式(6-3-5)的动量矩定理只适用于质心C 。

现由该式出发对力矩和动量矩进行变换以导出适用于o 点的动量矩定理表达式。

以k M N o =表示外力矩之和，其中N 是作用于船舶的绕z 轴的外力矩，以zz I 表示船舶对oz 轴的惯性矩，由理论力学的力矩和惯性矩移轴公式，有F ρM M ⨯+=C C o 及
2
C zz m I I ρζζ+=，这样由式(6-3-4)和式(6-3-5)可推出
k k j i k k )(r x ur v m x r I Y x r I N c c c &&&&+++=⨯+=ςςςς
N ur v mx r
I C zz =++)(&& (6-3-6)
式(6-3-6)即为船舶转动的动力学基本方程，其形状与式(6-3-5)的区别在于，左端的附加
项v mx C &及ur mx C 分别代表由于质心C 对原点o 做旋转运动所产生的离心惯性力矩和向心
惯性力矩。

6.4 船舶平面运动的线性化数学模型
综合式(6-3-3)和式(6-3-6)，得下列形式的船舶平面运动基本方程
⎪⎭
⎪
⎬⎫
=++=++=--N ur v mx r I Y r x ur v m X r x vr u m C zz C C )()()(2&&&&& (6-4-1)
当附体坐标系原点取在质心C 时，0=C x ，可得最简形式的船舶平面运动基本方程
⎪⎭
⎪
⎬⎫
==+=-N r I Y ur v
m X vr u m zz &&&)()( (6-4-2)
式(6-4-1)代表着3种力的平衡关系：左端是船体本身的惯性力和力矩，右端是流体对船体运动的反作用力，实际上包含了流体惯性力和力矩及黏性力和力矩。

式(6-4-1)本质是非线性的，其左端显式地出现vr ur ,等非线性项，尤其右端的N Y X ,,将是运动变量和控制变量的多元非线性函数，结构异常复杂。

(1)船舶平面运动的非线性模型和线性模型
船舶运动数学模型分线性化数学模型和非线性数学模型两大类。

研究船舶数学模型通常有两种目的：一种目的是建立精密程度不同的船舶运动仿真器(又称船舶运动模拟器)，用于通过仿真对船舶操纵特性进行研究，对船舶运动闭环控制系统进行研究，对船舶运动控制器性能进行评价。

这种模型必须是非线性的，以包含尽可能多的机理细节；另一种模型目的是用于船舶运动控制器设计，这种模型主要是线性的，因为迄今为止，线性反馈控制理论仍是能够提供各种控制器设计系统性方法的惟一控制论分支。

当引用神经网络控制或模糊控制时，非线性船舶运动数学模型可以提供训练和学习的数据。

1)船舶平面运动非线性数学模型 为应用方程式(6-4-1)求解船舶平面运动的基本变量r v u ,,，必须具体讨论流体动力X,Y 和力矩N 的结构形式。

研究中把船体、螺旋桨和舵视为
一个整体，此时X,Y,N 将是移动速度),(v u 、转动角速度)(r 、它们的时间导数),,(r v u
&&&、舵角)(δ以及螺旋桨转速)(n 的非线性函数
⎪⎭
⎪
⎬⎫
==),,,,,,,(),,,,,,,(),,,,,,,(n r v u r v u N N n r v u r v u Y Y n r v u r v u X X δδδ&&&&&&&&&＝
(6-4-3)
完全从理论上确定式(6-4-3)的函数关系极为困难，迫使研究者不得不转向半理论半经验的方法或多元数据回归方法。

Abkowitz 提出一种小扰动和Taylor 展开研究X,Y,N 的表示式的方法，其主要思路是，
考虑船舶等速直线运动这一平衡状态：0,0,0,0========r v u
r v V u u &&&δ，这时在式(6-4-3)中的自变量n 将不出现；从该点出发，研究偏离平衡状态不远的运动：
δδδ∆∆∆∆=∆=∆=∆=∆=∆=∆+=,,,,,,,,,,0Λ&&&&&&v u r r v v u u r r v v u u u 是小量；将
X,Y,N 在平衡点附近展成Taylor 级数时，在展开式中将仅出现v u &&,和r &的一次项，因为流体对船舶的惯性反作用力只取决于平移加速度v u &&,以及转动角加速度r &
本身，而与它们的各阶导数无关；至于和r v u ,,∆有关的黏性力各项及与δ有关的舵力各项，则取至3阶为止，更高阶的项全部略去。

将式(6-4-3)的展开式代回式(6-4-1)并进行移项整理，可得到Abkowitz
非线性船舶运动方程[6]。

Norrbin 发展一种非线性船舶运动数学模型[7]，该模型有两个特点，一是适用于运动变量),,(r v u 的整个变化范围；二是它不像Abkowitz 模型那样，完全按数学方式处理流体动力，以至其Taylor 展开式的某些项缺乏物理意义，而是在更深的层次上依赖于流体动力学的基本原理，构成一种半理论半经验的模型格局。

以上所述的Abkowitz 和Norrbin 船舶运动非线性数学模型属于“整体式”模型，本节将做较详细的介绍。

与此相对应，日本船舶操纵数学模型小组(Manoeuvring Model Group, MMG)提出了一种分离式船舶运动数学模型[8]，后者是在单独考虑船体、螺旋桨、舵的流体动力学特性的基础上再研究在它们构成一个推进和操纵系统时，各部分之间的相互干扰。

这种分离式模型的优势是具有完整的理论支持，易于进行实验研究从而获得较为通用的数据回归结果，对于希望建立自己的复杂程度不同的船舶操纵模拟器的各类研究人员均有裨益。

有关MMG 模型的结构和细节，有兴趣的读者可参考文献[9]。

2)船舶平面运动线性数学模型 沿用Abkowitz 的研究方案，在把流体动力X,Y,N 展开成Taylor 级数时只保留一阶小量[6]，同时在船舶运动基本方程左端也进行线性化处理，从而得到平面运动线性数学模型，有
δδδδN r N v N r N v N r u mx v mx r I Y r Y v Y r Y v Y r mx r mu v m u X u X u m r v r v c c zz r v r v c u u ++++=++++++=++∆+∆=∆&&&&&&&&&&&&
&&&00 以矩阵形式表示之，有
δδδ⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡+⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡∆⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡∆⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-----N Y r v u u mx N N mu Y Y X r v u N I N mx Y mx Y m X m C r v
r v u r zz v C r C v u 0)(0
)(000)()(0)()(000)
(00&&&&&&&&(6-4-4
)
式(6-4-4)对研究平面运动稳定性有用。

3)前进运动与横漂、转首运动的解耦 式(6-4-4)表明，在线性化前提下，前进运动与其他两个自由度上的运动互相独立，从航速控制的角度，该自由度的运动可以单独考虑；横漂及转首运动之间存在着强耦合，这两个自由度上的运动与船舶航向、航迹控制密切相关，是
本章研究的重点，故而将式(6-4-4)重新写为
u X u
X m u u ∆=-&&)( (6-4-5)
δδδ
⎥⎦
⎤
⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡--=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡----N Y r v u mx N N mu Y Y r v N I N mx Y mx Y m C r v r v r zz v C r C v )()()()()()
(00&&&&&&
(6-4-6)
4)流体动力导数的无量纲化 船舶线性化数学模型的进一步推演主要涉及10个流体动
力导数δδN Y N N Y Y N N Y Y r v r v r v r v ,,,,,,,,,&&&&，前4个称为“速度导数”，第5～第8个称为“加速度导数”，最后两个称为“舵力和舵力矩导数”。

由于船舶(包括桨、舵)几何形状的复杂性，应用理论流体动力学方法计算这些流体动力导数是不可能的，因此它们的确定必须依赖于船模试验。

为了数据处理的科学性以及使用的方便性，根据相似原理和量纲分析方法，
应该采用无量纲的流体动力导数。

为此选择一些基本的度量单位：长度0L ──L (船长)，速度0V ──V (航速)，时间0t ──L/V ，质量0m ──3)2/1(L ρ，力0F ──22)2/1(L V ρ，力矩0M ──32)2/1(L V ρ，其中ρ──水密度。

这样将得到各量的无量纲值：
质量：)2
1
/(3'L m m ρ= 长度：L x x C C /'=
速度：V v v /'= 转首角速度：V rL r /'= 力：)2
1/(22'L V F F ρ= 力矩：)2
1
/(23'V L N N ρ=
惯性矩：)21/(5'
L I I zz zz
ρ= 16
2
mL I zz =
以此为基础，将得到无量纲速度导数v v v Y VL Y F V V v F Y v Y Y 2
0000'
''2
)/()/(/ρ==∂∂=
∂∂=,
'''/r N N r ∂∂=Λ,))/(2(4r N VL ρ=，以此类推。

以上介绍的无量纲化流体动力导数称为“一
撇”系统(prime system)，由美国造船与轮机工程师协会(SNAME)于1950年提出；此后Norrbin
又提出了“两撇”(bis system)[7]，其独到之处是采用与上述不同的基本度量单位，如：长度──L ，速度──gL ，时间──g L /，质量──∇∇,ρ为排水体积，力──∇g ρ，力
矩──L g ∇ρ。

由此得出的无量纲流体动力导数以Λ,,"
"r
v N Y 表示。

5)线性流体动力导数的估算公式 Clarke 整理大量船模试验数据，给出关于10个线性
流体动力导数的回归公式[10]，汇集如下：
2
2 ')/(π])/(1.5/16.01[L T L B T B C Y b v ⋅-+-=& 22
')/(π])/(0033.0/67.0[L T T B L B Y r ⋅--=&
2
')/(π]/041.0/1.1[L T T B L B N v ⋅--=&
2
')/(π]/33.0/017.012/1[L T L B T B C N b r ⋅-+-=&
2 ')/(π]/40.01[L T T B C Y b v ⋅+-= (6-4-7)
2')/(π]/080.0/2.22/1[L T T B L B Y r ⋅-+--=
2
')/(π]/4.22/1[L T L T N v
⋅+-= 2' )/(π]/56.0/039.04/1[L T L B T B N r ⋅-+-= 2 '/0.3L A Y δδ= '' )2/1(δδY N -=
上式中B,T,b C ,A δ分别指船宽、吃水、方形系数、舵叶面积。

上式中的''
'',,,r v
r v N N Y Y 是船体本身的流体动力导数，在实际应用时应考虑舵对船体流体动力的干扰，尚需对这些流体动力
导数做一定的修正，需修改的增量按下式确定[10] ' 'δγY Y v -=∆
' '21
v r Y Y ∆-
=∆ '' 21
v v Y N ∆-=∆ (6-4-8)
' '4
1
v r Y N ∆=∆
30.0=γ
(2)状态空间型船舶平面运动数学模型
状态空间型的船舶运动数学模型是船舶运动控制器设计的基础，它可以有多层次的模型化方案，不同维数的模型用于不同的设计目的和精度要求，详见文献[9]。

1)二自由度状态空间型船舶线性数学模型 在式(6-4-6)的第一行两端除以32
1L ρ，第二
行除以
42
1L ρ，并转化成无量纲流体动力导数，则有
δδδ⎥⎥⎥⎥
⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎥⎥⎦⎤
⎢⎢⎢⎣⎡--=⎥⎥⎦⎤⎢⎢⎣⎡⎥⎦⎤⎢⎣⎡----'2'2''''
'
'
'
.
.
''''
'''''
')()()()(N L V Y L V r v x m N V N L
V m Y V Y L V r v N I L N x m Y x m L Y m C r v
r v
r zz v
C r C v
&&&&
(6-4-9)
上式可简记为
U Q X P X I ')2()2(')2()2(.
')2(+=
(6-4-10)
其中
⎥⎦⎤⎢⎣⎡----=)()(''''
''
''''')
2(r zz v
C r C v N I L N x m Y x m L Y m &&&&
I
⎥
⎥⎥
⎦
⎤
⎢⎢⎢⎣⎡--=)()('
''''''')
2(C r v r v x m N V N L V m Y V Y L V P
⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡='2'2 ')
2(δδN L V Y L V Q 分别是惯性力导数矩阵、黏性力导数矩阵及舵力导数矩阵，T )2(][r v =X 是状态向量，δ=U 是控制输入。

将式(6-4-10)化成标准的状态空间形式，得
δ)2()2()2()2(.
B X A X += (6-4-11)
其中
⎥⎦
⎤
⎢⎣⎡==-22211211
')2(1 ')2()2()(a a a a P I A ⎥⎦
⎤⎢⎣⎡==-2111'
)2(1' )2()2()(b b Q I B
并且
L
Y x m N x m Y m N I S S L V N Y m Y N x m b S V N Y x m Y N I b S V x m N Y m m Y N x m a S L V N Y m Y N x m a S LV x m N Y x m m Y N I a S V N Y x m Y N I a r C v C v r zz v v C r C r zz C r v r v C v v v v C C r r C r r zz v r C v r zz )])(())([(//])()([/])()[(/)])(())(([//])()([/)])(())([(/])()[(''''''''''112'''''''2112'
''''''111'
'''''''''221
'''''''211''''''''''121
''
'''''11&&&&&&&&&&&&&&&&-----=-+--=---=--+---=-+--=-----=---=δδδδ
(6-4-12)
2)三自由度状态空间型船舶线性数学模型 在式(6-4-11)的基础上，增加一个便于研究问题的状态变量ψ∆(航向偏差)，且r r ψψψψ,-=∆为设定航向，使状态向量成为
T
)3(][ψ∆=r v X 。

因r =∆.
ψ，可得
δ)3()3()3()3(.
B X A X += (6-4-13) 其中 ⎥⎥⎦
⎤⎢⎢⎣⎡
=01000)2()
3(A
A ⎥⎦⎤⎢⎣⎡=0)2()3(
B B
3阶模型是最基本的，由此可演化成其他更高阶的模型形式，直接利用3阶模型可进行线性二次型(Linear Quadratic, LQ)最优控制器设计。

3)四自由度状态空间型船舶线性数学模型 在式(6-4-13)基础上再叠加以舵机伺服系统的模型，后者一般被视为一个1阶惯性环节，其时间常数为T r ，则有 r r
r T T δδδ1
1.
+-
= (6-4-14) 其中：δr 为命令舵角，则状态变量成为T )4(][δψ∆=r v X ，可得到
r δ)4()4()4()4(.
B X A X += (6-4-15)
其中 ⎥⎥⎥⎦
⎤
⎢⎢⎢
⎣⎡-=r T 1000)3()3()
4(B A A T )4(1
|000⎥⎦⎤⎢⎣⎡=r T B 4)考虑随机干扰时的线性船舶数学模型 考虑海上环境干扰对船舶的影响，并把这种干扰简化为一种白噪声[]T
21)2(w w =w ，
则船舶运动数学模型将从确定性系统变为随机系统，这样有
[]⎪⎭
⎪⎬⎫
=++=T
21)2()2()2()2()2()2(.
w w w w B X A X δ (6-4-16)
[]⎪⎭⎪
⎬⎫=++=T
321)3()3()3()3()3()3(.
w w w w w B X A X δ (6-4-17) []⎪⎭
⎪
⎬⎫=++=T
4321)4()4()4()4()4()4(.
w w w w r w w B X A X δ (6-4-18)
其中白噪声3w 代表航向角ψ受到的高频噪声，4w 代表海浪对舵叶驱动伺服系统的干扰作用。

(3)传递函数型的船舶运动数学模型
传递函数型数学模型在经典控制论以至智能控制范畴内用于分析船舶运动的动态行为，
并且可作为设计航向、航迹控制器的基础。

1)3阶传递函数模型 对于船舶航向控制来说，采用3个自由度的状态空间数学模型式(6-4-12)，加上输出方程
)3(CX =m ψ (6-4-19) 其中m ψ为量测航向，[]100=C ，将此状态空间模型转换成传递函数形式为
)
1)(1()
1(][)(21301+++=
-=-s T s T s s T K s s G B A I C ψδ
(6-4-20)
这是一个3阶系统，具有两个非零极点和一个零点，且有
211222112
11
a a a a T T -= 21
2
13
01121211121322112
12
1)(11)(b T T T K a b a b b T a a T T T T =-=+-=+ 由此不难解得3个时间常数321,,T T T 以及一个系统增益系数0K 。

2)2阶传递函数模型(Nomoto 模型) 野本(Nomoto)对3阶船舶模型式(6-4-20)做了一项出色的简化工作，使之降为2阶[11]。

论证的出发点在于，对于船舶这种大惯性的运载工具来说，其动态特性只在低频段是重要的，故在式(6-4-20)中令0j →=ωs ，且利用一个熟知的近似关系：当0→x 时有)1/(1)1(x x +≈-，并忽略2阶和3阶小量，由此导出著名的Nomoto 模型
)
1()(00
+=
s T s K s G ψδ (6-4-21)
其中增益0K 与3阶模型相同，时间常数3210T T T T -+=，或直接由下式求出 21
122211112121110a a a a a b a b K --=
1121211121
2112221122110a b a b b a a a a a a T --
-+-= 式(6-4-21)广泛应用于船舶自动舵的控制器设计中。

用Nomoto 模型进行船舶运动控制器设计有两个好处：一是在低频范围，其频谱与高阶模型的频谱非常相近；二是设计出的控
制器阶次低，易于实现。

求解本节所述船舶运动数学模型需要已知8个船舶参数，即航速V ，两柱间长L ，船宽
B ，满载吃水T ，方形系数b
C ，排水量∇，重心距中心距离C x ，舵叶面积δA 。

首先将这8个已知参数代入式(6-4-7)，求出10个流体动力导数，并用式(6-4-8)修正，然后代入式(6-4-12)，即可求出各种自由度的数学模型。

6.5 船舶平面运动的一种简洁非线性数学模型
(1)用于船舶运动闭环控制系统仿真的六自由度非线性模型
各种线性船舶数学模型只用于在不同情况下进行控制器设计，当用于船舶闭环控制系统仿真研究时，必须以非线性模型表述被控过程的动态特性，并且还需考虑风、浪、流造成的
环境干扰。

从式(6-4-10)出发，在其右端加上非线性流体动力项NON F 、风力项WIND F 、浪力项WAVE F ，则无量纲的二自由度非线性船舶运动数学模型将呈下列形式 '
WAVE 'WIND 'NON ')2()2(')2()2(.
')2(F F F U Q X P X I ++++=
(6-5-1) 其中
⎥⎥⎦
⎤⎢⎢⎣⎡=421
NON 3
2
1NON ' NON
//L N L Y ρρF
⎥⎥⎦
⎤⎢⎢⎣⎡=421
WIND 32
1WIND '
WIND
//L N L Y ρρF
⎥⎥⎦
⎤⎢⎢⎣⎡=421
WAVE 32
1WAVE '
WAVE
//L N L Y ρρF NON Y ,WIND Y ,WAVE Y 及NON N ,WIND N ,WAVE N 分别是非线性力、风力、浪力在y 方向的合力及
在绕z 轴方向的合力矩。

由式(6-5-1)和式(6-4-14)不难看出在T )4(][δψ∆=r v X 的4个自由度上有非线性状态方程
⎥⎥⎥
⎦
⎤⎢⎢⎢⎣⎡++++=-00][)'('
WAVE 'WIND 'NON 1 )2()4()
4()4()4(.
F F F I B X A X r δ (6-5-2)
考虑到船舶位置T 00],[y x 的两个自由度上的运动学关系
⎭
⎬⎫
+=-=ψψψψcos sin sin cos 00v u y v u x && (6-5-3)
则式(6-5-2)与式(6-5-3)构成了六自由度的船舶运动非线性数学模型的基本框架，状态变量变为 T 00)6(][y x r v δψ∆=X
各研究者关于式(6-5-2)中非线性流体动力'
NON F 的取法不同是区别到目前为止形形色色的非
线性船舶运动数学模型的主要标志。

(2)Norrbin 关于非线性力的简化表示式
Norrbin 在研究船舶参数辨识问题时提出了一种关于非线性流体动力'NON F 的简洁表示
式[7,12]，如下所示：
⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡=),(),('NON 'NON 'NON
r v Cf r v Cf N Y N
Y F
(6-5-4) 其中
⎪⎪⎪
⎩
⎪
⎪⎪⎨⎧
∞<-<+⋅≤-≤---⋅-<-<∞---⋅=r v L r v L r r T r v L r v L r v L r r T r v L r v L r r T r v f Y 121])(1121[21121])(132121[211])(1121[),(2233
22 (6-5-5)
⎪⎪⎪
⎩
⎪
⎪⎪⎨⎧
∞
<-<⋅≤-≤-+--⋅-<-<∞--⋅=r v L r v L r r T r v L r v L r v L r r T r v L r v L r r T r v f N 121)](161[21121])(161)(141321[211)](161[),(44
22 (6-5-6)
式(6-5-4)中的比例系数C 为无量纲横流系数，其值通常在0.3～0.8范围内。

Norrbin 关于'
NON
F 的横流模型式(6-5-5)、式(6-5-6)的优越之处在于其表示式在各类非线性模型中最为简单，它的导出具有比较明确的理论基础，并且公式中除了船舶吃水和船长之外，不需任何关于船体结构的数据，应用甚为方便。

据笔者的经验，由式(6-5-2)～式(6-5-6)组成简化的非线性船舶运动数学模型用于在自动舵控制下的闭环系统的仿真研究，结果是可信的[13~15]。

应指出，对式(6-5-5)和式(6-5-6)中同时出现0=v 和0=r 的情况应做专门处理。

(3)风力干扰
在式(6-5-2)中，风力'WIND 'WIND ,N Y 分成平均风力'
WIND 'WIND ,N Y 及脉动风力
'WIND 'WIND ~,~N Y [16]。

平均风力计算见图6.5.1。

图6.5.1 平均风力计算
平均风力的表示式如下：
⎪⎭
⎪⎬⎫==421221'WIND 321221'
WIND /)(/)(L L A V C N L A V C Y L R A R N L R A R Y ρργρργ
00X O
(6-5-7)
式(6-5-7)中)(),(R N R Y C C γγ为无量纲的风力和风力矩系数，文献[17]给出了这两个系数的一系列图谱可资利用，其估算公式参见文献[16]；L A 为船舶水线以上侧投影面积，A ρ为空气密度；R V 为相对风速，R γ为相对风速与首向间的夹角，称为风舷角，由绝对风速T V 、绝对风向WIND α以及航速V (u,v )按下式计算
222WIND WIND )sin()cos(R
R R T R T R v u V v V v u
V u +=---=---=ψαψα
(6-5-8)
⎪⎩
⎪⎨⎧>>-<>=+=其他
00,00,0arctan
R R R R R
R
R v u v u u v π
πυυγ 上式中ψα,WIND 变动范围为0︒～360︒，R γ变动范围为0︒～±180︒，相对风从右舷来时R γ>0。

脉动风力'
WIND 'WIND ~,~N Y 是由大气的湍流所造成的，按文献[16]，它们被认为是某种白噪
声的实现，该白噪声的标准差Y σ,N σ与绝对风速T V 的平方成正比
3
2
22)(2.0)(2.0L
C V L C V R N T
A N R Y T A Y γρ
σγρσ== (6-5-9)
(4)浪力干扰
浪力'
WAVE 'WAVE ,N Y 分为两个组成部分：高频的一次力，它是与波浪宏观振荡运动同步
的周期力，幅值可较船舶的推进力或因运动而产生的流体动力高一个数量级，但由于大惯性
船舶本体的滤波作用，一次力产生的振荡运动(艏摇、横荡等)被限制在允许范围之内；低频的二次力，数量级较小，数值变动缓慢，产生船位的漂移。

1)一次力的计算 采用文献[16]的结果，把波浪看成规则波，这种波浪只有一个频率ω、一个周期ωT 和一个波高ωh ；而把船舶看成一个简单的六面体；在小扰动假设下压力由波形抬高按Bernoulli 公式求出，浪力是在船体水下表面上把压力积分而得，并表成封闭的解析形式。

更准确地可采用不规则波概念，把不同风力下的波谱分解成一系列波谱段(例如10段)，每一段波谱对应着一定的频率和波高，这样不规则波就由一系列规则波叠加而成；船体也被分解成一系列六面体分段(例如20段)；分别计算各种波浪分量在每一分段上的力，最后按频率和船长进行二维求和可得到总的浪力，但计算量大为增加[18]，未予采用。

规则波对于船的传播方向称为浪向角，以χ表之，参见图6.5.2，有
)(πWIND ψαχ--= (6-5-10)
χ=0为顺浪，χ＝π为顶浪，2π±=χ为横浪(2
π
+
表示浪从右舷来)；船对波浪的遭遇频率是 ψχωωsin cos kv ku e +-=
(6-5-11)
图6.5.2 浪向角χ
式(6-5-11)中ω为规则波本身的圆频率，k 为波数，有
ωωωL g gT k π
2π422
2=== (6-5-12) 式中：ωL 为波长，T ω为波浪周期，与风速有关，其具体的依存关系视考察的海区而有所不同。

Kallstrom 根据海上观测数据进行最小二乘回归[16]，得到)(T V T ω和波高)(T V h ω公式如下：
⎭
⎬⎫
+=++-=5.1015.0)(6.5042.00014.0)(223T T T T T V V h V V V T ωω (6-5-13)
注意式(6-5-13)只适用于,m/s 20≤T V 对于m/s 20>T V 的情况应谨慎进行外推处理；并且在T V ＝0时仍给出1.5 m 的波高和5.6 s 的波浪周期，这对大西洋上的情况可能是合适的，但用于我国近海海域，可能稍有误差。

Kallstrom 还对Zuidweg 的工作[19]略加修改，给出浪
力表示式如下：
)()sin cos sin sin cos sin ()(sin sin 22
2
22WAVE WAVE t b
b b b
c L c c c c b B ak N t s b
c
b aL
Y ξ-⋅--⋅=⋅-=
(6-5-14)
0X O
其中 ⎪⎪
⎪⎪
⎭
⎪
⎪⎪⎪
⎬⎫
==⋅=⋅=-=-)cos()2()()sin()2()(sin 2/cos 2//)e 1(2t h t t kh t s kB c kL b k g a e e kT ωξωχχρω
ω
(6-5-15) 式(6-5-14)中)(t ξ代表在附体坐标系原点处波面的振荡，)(t s 则表明沿波浪传播方向上的波面ξ的斜率在原点处的值。

在进行仿真时应对式(6-5-15)的e ω进行适当的滤波[16]，滤波的方法如下：
)())(()1())(()(k k S B k k S A k e ef ef ωωω++--= (6-5-16) )
1()
1()(-+-⋅=
k S A k S A k S (6-5-17)
式(6-5-16)、式(6-5-17)表明所采用的是一个时变滤波器。

ef ω为经过滤波的遭遇频率；e ω为未经滤波的遭遇频率；A ,B 是两个常数，取A =0.999，B =0.001(即A +B =1)；S (0)取为0.999，随着递推次数k 的增加，由式(6-5-17)知S (k )下降，当k →∞时S (k )→0，此时式(6-5-16)趋于一个定常滤波器
)(001.0)1(999.0)(k k k e ef ef ωωω+-= (6-5-18)
滤波的结果是：在连续的采样周期ef ω也连续变化。

2)二次力的计算 目前尚无简捷而可靠的方法。

海浪干扰的另一种简单模拟方法是用白噪声驱动一个典型的2阶振荡环节(相当于2阶低通滤波器)[20]。

其中白噪声的带宽为0.5 Hz ，2阶滤波器具有低阻尼比，参数为0.05，自然频率
g U n /cos 2
00γωωω-=
其中0ω为海浪频谱的峰值频率，U 为船速，γ为航向与海浪方向之间的夹角，g 为重力加速度。

例如，如果模拟的海况为5级风，中浪，参数可取0ω为0.15 rad/s ，船速为7 m/s(约14 kn)，γ为60︒。

(5)流干扰
仿真时通常假定流是恒定并且均匀的，它只改变船舶运动的位置和速度，而不改变船舶
的航向，有下列速度平衡方程：
c
c c
c V v u y V v u x γψψγψψsin cos sin cos sin cos 00++=+-=&&
(6-5-19)。