船舶在遭遇恶劣海况时，为了避免遭受横风横浪的影响，通常会调整航向，选择纵浪或斜浪航行^[1]。即使船舶在静水中有足够的稳性，在纵浪和斜浪航行时仍有发生大角度摇摆乃至倾覆的可能。特别是当船波遭遇频率与船舶横摇固有频率比为2时，很有可能发生参数横摇甚至倾覆。参数横摇发生的基本原理是船舶在航行过程中，由于波长近似等于船长，可认为波峰分别经过船首-船中-船尾整个过程是一个完整的参数激励周期。当波峰位于船首/尾(即波谷位于船中)时，船舶首尾处的水线面要比静水中的水线面大，即此时的船舶横摇复原力矩要比静水时的偏大；而当波峰位于船中时，船舶首尾处的水线面要比静水中的水线面小，即此时的船舶横摇复原力矩要比静水时的偏小。随着船-波相对位置的周期性变化，船舶横摇复原力矩也呈现出周期性变化，从而引发显著的参数横摇运动。因此，有必要研究纵浪和斜浪中航行船舶的稳性及运动状态，分析船舶倾覆的机理，为避免船舶倾覆提供操船决策。

早在1863年，Froude^[2]就观察到在纵浪中航行的船舶受到波浪激励会发生大幅横摇的现象，只是当时的知识还无法解释该现象。随着非线性理论的发展以及计算机技术的飞速发展，各国学者针对船舶参数横摇运动进行了大量的研究并取得了一些新的成果。Matusiak^[3]对参数共振的理论进行了研究。他提出了完整的六自由度数学模型。其研究表明：产生横摇参数共振的最主要的原因是傅汝德-克雷洛夫力和复原力矩的非线性特性。Silva等^[4-5]提出了基于切片理论研究迎浪船舶动稳性的新方法，并将该方法拓展到船舶六自由度非线性运动模型，通过与试验结果对比，验证了该方法的可行性。Silva和Soares^[6]采用横摇-纵摇-垂荡耦合的数学模型，进行了数值模拟和运动稳定性分析，并用模型试验对研究结果加以验证。研究发现：如果船舶遭受波浪的波长等于其船长、波浪的频率为横摇固有频率的两倍时，将有可能出现严重的横摇参数共振现象，其横摇角可达30°。Umeda^[7]等考虑横摇恢复力矩的变化研究了规则和随机纵浪和斜浪中船舶的参数激励横摇运动，其横摇恢复力矩系数的变化根据倾覆实验求出。Vidic-Perunovic等^[8]采用一阶可靠性方法计算船舶参数横摇运动超过限定角度的概率。鲁江、卜淑霞等^[9]基于切片理论，提出了复原力的计算方法，计算过程中考虑了复原力中的辐射力和绕射力成分，并通过模型试验验证了数值计算方法的有效性。Dostal^[10]针对随机海况，应用能量包线随机平均法，理论推导了船舶参数横摇响应的概率密度函数。彭东升等^[11]针对C11集装箱船，计算了不同阻尼下船舶参数横摇响应。傅超等^[12]对C11集装箱船进行了参数敏感性分析，提出通过改进船舶设计减少参数横摇发生的概率。Wei Chai^[13]于2016年提出了一种能有效计算蒙特卡罗模拟的方法，并验证该方法适用于统计船舶在随机长峰波激励下的参数横摇响应极值。

本文在前人工作的基础上研究随机斜浪中船舶参-强激励横摇响应的概率密度函数。为实现解析求解，将数值模拟得到的复原力臂函数用近似的解析函数表示，将随机海浪谱处理为窄带随机过程以得到随机参数激励和随机强迫激励的谱密度函数，采用能量包线随机平均法计算系统响应的稳定概率密度函数。

1 船舶参数激励横摇运动方程

考虑随机斜浪，假设船舶的垂荡及纵摇运动为小量，船舶横摇运动方程如下：

$ \left( {I + A\left( {{\omega _n}} \right)} \right)\mathit{\ddot \Phi } + {b_1}\mathit{\dot \Phi } + {b_3}{{\mathit{\dot \Phi }}^3} + g\Delta GZ\left( {\mathit{\Phi },\eta ,\mathit{\Psi }} \right) = M。$

(1)

其中：Φ为横摇角；I为横摇惯性矩；A(ω_n)为横摇附加惯性矩；b₁与b₃分别为线性阻尼系数与立方非线性阻尼系数；Δ为排水量；g为重力加速度；M为波浪强迫激励力矩；GZ(Φ, η, Ψ)为船舶复原力臂，由Φ、η和Ψ决定，其中η为长峰波的波动幅度，Ψ为船和波的相对位置，取值范围为[0，2π]。

船舶在长峰波中的复原力臂GZ(Φ, η, Ψ)可基于切片理论求解，一般采用数值的方法进行模拟^[9]。

为了能够使用随机平均理论求解运动方程(1)，将GZ(Φ, η, Ψ)展开为如下的关于横摇角Φ的多项式^[14]：

$ G{Z_a}\left( {\mathit{\Phi },\eta ,\mathit{\Psi }} \right) = {q_1}\mathit{\Phi } + {q_2}{\mathit{\Phi }^3} + {q_3}{\eta _c}\mathit{\Phi }。$

(2)

式中，η_c=ηcos(Ψ)，是由随机波浪确定的随机过程，可通过窄带海浪谱展开得到；q_i(i=1, 2, 3)为展开项的系数，使用最小二乘法确定，即使下式取得最小：

$ E = \frac{1}{{4{\rm{ \mathsf{ π} }}{\mathit{\Phi }_0}{\eta _0}}}\int_0^{2{\rm{ \mathsf{ π} }}} {\int_0^{{\eta _0}} {\int_{ - {\mathit{\Phi }_0}}^{{\mathit{\Phi }_0}} {{{\left( {GZ - G{Z_a}} \right)}^2}} } } {\rm{d}}\mathit{\Phi }{\rm{d}}\eta {\rm{d}}\mathit{\Psi }。$

(3)

2 窄带随机海浪的数学描述

对于随机海浪，一般用谐波叠加法来模拟不规则长峰波波面函数。为采用随机平均法研究船舶随机参-强激励横摇运动，以下将随机波面升高处理为窄带随机过程Z_e(x, t)，其表达式如下^[10]：

$ {Z_e}\left( {x,t} \right) = {\eta _1}(t)\cos (2{\rm{ \mathsf{ π} }}x/L) - {\eta _2}(t)\sin (2{\rm{ \mathsf{ π} }}x/L)。$

(4)

其中，η₁(t)和η₂(t)均为高斯平稳随机过程，两者互不相关、统计独立^[15-16]，表达式如下：

$ {\eta _1}(t) = \int_0^\infty {\frac{{2{\rm{ \mathsf{ π} }}\sin r}}{{{{\rm{ \mathsf{ π} }}^2} - {r^2}}}} \sin (\omega t + \zeta (\omega ))\sqrt {2S(\omega ){\rm{d}}\omega } ; $

(5)

$ {\eta _2}(t) = \int_0^\infty {\frac{{2r\sin r}}{{{{\rm{ \mathsf{ π} }}^2} - {r^2}}}} \cos (\omega t + \zeta (\omega ))\sqrt {2S(\omega ){\rm{d}}\omega } 。$

(6)

式中：ω为海浪频率；S(ω)为海浪谱；ζ(ω)为随机相位角。r=(L/2)k，其中L为特征波长，k为波数。

η₁(t)和η₂(t)对应的谱密度函数分别为：

$ {{\bar S}_1}(\omega ) = 2S(\omega ){\left( {\frac{{2{\rm{ \mathsf{ π} }}\sin r}}{{{{\rm{ \mathsf{ π} }}^2} - {r^2}}}} \right)^2}; $

(7)

$ {{\bar S}_2}(\omega ) = 2S(\omega ){\left( {\frac{{2r\sin r}}{{{{\rm{ \mathsf{ π} }}^2} - {r^2}}}} \right)^2}。$

(8)

如果船舶以航向角β、航速U航行，则式(4)~(8)中的海浪频率ω和海浪谱S(ω)要用对应的遭遇频率ω_e和遭遇谱S(ω_e)来替换，其表达式为：

$ {\omega _e} = \omega - kU\cos \beta ; $

(9)

$ S\left( {{\omega _e}} \right) = S(\omega )\left| {\frac{1}{{\sqrt {1 - 4U{\omega _e}\cos \beta /g} }}} \right|。$

(10)

进一步将随机过程η₂用受控自回归滑动平均模型(CARMA(2, 1))来模拟。对于CARMA(2, 1)过程，需要满足如下形式的微分方程：

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{d}}{u_1} = \left( {{u_2} - {c_2}{u_1}} \right){\rm{d}}t + {c_1}\bar \sigma {\rm{d}}W}\\ {{\rm{d}}{u_2} = - {c_3}{u_1}{\rm{d}}t} \end{array}} \right.。$

(11)

式中：u₁、u₂为状态变量；${\bar \sigma ^2}$为噪声强度；W为标准维纳过程；c₁、c₂、c₃为系统参量。对于标准维纳过程，噪声强度${\bar \sigma ^2}=1$。该CARMA(2, 1)过程的传递函数与谱密度函数分别为：

$ {H_2}(s) = \frac{{{c_1}s}}{{{s^2} + {c_2}s + {c_3}}}; $

(12)

$ {S_2}(\omega ) = {H_2}(s){H_2}( - s)\frac{{{{\bar \sigma }^2}}}{{2{\rm{ \mathsf{ π} }}}}。$

(13)

其中，s=iω。采用最小二乘法，使式(13)与式(8)的误差最小，从而确定c_i(i=1, 2, 3)。

3 船舶随机参-强激励横摇概率密度函数

采用能量包线随机平均法研究船舶随机参-强激励横摇运动概率密度函数。将船舶运动过程处理为马尔可夫过程，将响应量分成快变量与慢变量，通过对快变量的平均得到慢变量的近似方程，从而使系统的自由度缩减。

将GZ_a(Φ, η, Ψ)代入船舶参数横摇运动方程(1)，将随机波浪激励力矩近似处理为M≈Δgq₀πsin(β)η_c(t)/L，令Φ=u、η_c=ξ_t，并对方程(1)做无因次处理，有：

$ \left\{ {\begin{array}{*{20}{l}} {\dot {\bar u} = \bar v}\\ {\dot {\bar v} = - \overline {{\alpha _1}} \bar u + \overline {{\alpha _3}} {{\bar u}^3} - \overline {{\beta _1}} \bar v - \overline {{\beta _3}} {{\bar v}^3} + (\overline {{\gamma _1}} + \bar u\overline {{\gamma _2}} ){\xi _t}} \end{array}} \right.。$

(14)

其中，$\overline {{\alpha _1}} = \frac{{g\Delta {q_1}}}{{{I_{xx}} + {A_{xx}}\left( {{\omega _n}} \right)}}$，$\overline {{\alpha _3}} = - \frac{{g\Delta {q_2}}}{{{I_{xx}} + {A_{xx}}\left( {{\omega _n}} \right)}}$，$\overline {{\beta _1}} = \frac{{{b_1}}}{{{I_{xx}} + {A_{xx}}\left( {{\omega _n}} \right)}}$，$\overline {{\beta _3}} = \frac{{{b_3}}}{{{I_{xx}} + {A_{xx}}\left( {{\omega _n}} \right)}}$，$ \overline {{\gamma _1}} = \frac{\pi }{L}\frac{{\Delta g{q_0}\sin (\beta )}}{{{I_{xx}} + {A_{xx}}\left( {{\omega _n}} \right)}}$，$\overline {{\gamma _2}} = - \frac{{g\Delta {q_3}}}{{{I_{xx}} + {A_{xx}}\left( {{\omega _n}} \right)}}$。q₀为常数，本文选取q₀为船舶的初稳心高。

做如下的尺度变换：u=u、v=εv、$\overline {{\alpha _1}} = {\varepsilon ^2}{\alpha _1}$、$\overline {{\alpha _3}} = {\varepsilon ^2}{\alpha _3}$、$\overline {{\beta _1}} = {\varepsilon ^2}{\beta _1}$、$\overline {{\beta _3}} = {\beta _3}$、$\overline {{\gamma _1}} = {\varepsilon ^{5/2}}{\gamma _1}$、$\overline {{\gamma _2}} = {\varepsilon ^{5/2}}{\gamma _2}$、t_ε=εt，为了书写简便，以下用t替代t_ε，则运动微分方程(14)变为：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}u}}{{{\rm{d}}t}} = v}\\ {\frac{{{\rm{d}}v}}{{{\rm{d}}t}} = - {\alpha _1}u + {\alpha _3}{u^3} - \varepsilon \left( {{\beta _1}v + {\beta _3}{v^3}} \right) + }\\ {\;\;\;\;\;\;\sqrt \varepsilon \left( {{\gamma _1} + u{\gamma _2}} \right){\xi _t}} \end{array}} \right.。$

(15)

将响应分量进一步写为快变量和慢变量，即系统由状态空间(u, v)变换为状态空间(u, H)，其中H为系统总能量，其表达式为：

$ H(u,v) = \frac{{{v^2}}}{2} - \frac{1}{4}{\alpha _3}{u^4} + \frac{1}{2}{\alpha _1}{u^2}。$

(16)

式中，每个能量值H(u, v)对应于特定的横摇运动。对于确定的α₁、α₃值，可得到对应于式(16)的能量等值线，用b(H)表示在特定能量值H下，系统所能达到的最大横摇角，以下用b(H)来度量船舶的横摇运动。

定义船舶横摇的稳性消失角如下：

$ {B_1} = \left( {\sqrt {{\alpha _1}/{\alpha _3}} ,0} \right)。$

(17)

以该稳性消失角为临界横摇角，其所在的闭合曲线能量值为临界能量H_c，有H_c=α₁²/(4_α3)。

使用如下的转换关系：

$ {\left( {\frac{{{\rm{d}}u}}{{{\rm{d}}t}}} \right)^2} = 2H - 2\left( { - \frac{1}{4}{\alpha _3}{u^4} + \frac{1}{2}{\alpha _1}{u^2}} \right)。$

(18)

则运动方程(15)进一步写为：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}u}}{{{\rm{d}}t}} = \sqrt {2H - {\alpha _1}{u^2} + \frac{1}{2}{\alpha _3}{u^4}} }\\ {\frac{{{\rm{d}}H}}{{{\rm{d}}t}} = \left[ { - {\beta _1}{{\left( {\frac{{{\rm{d}}u}}{{{\rm{d}}t}}} \right)}^2} - {\beta _3}{{\left( {\frac{{{\rm{d}}u}}{{{\rm{d}}t}}} \right)}^4}} \right]\varepsilon + \left( {\frac{{{\rm{d}}u}}{{{\rm{d}}t}}} \right)\left( {{\gamma _1} + u{\gamma _2}} \right)\sqrt \varepsilon {\xi _t}} \end{array}} \right.。$

(19)

引入中间变量Q(u, H)=v²，有：

$ Q(u,H) = 2H - {\alpha _1}{u^2} + \frac{1}{2}{\alpha _3}{u^4}。$

(20)

则式(19)可进一步写为如下的随机微分方程式：

$ \left\{ \begin{array}{l} {\rm{d}}u = \sqrt {Q(u,H)} {\rm{d}}t\\ {\rm{d}}H = Q(u,H){\left( { - {\beta _1} - {\beta _3}Q(u,H)} \right)_\varepsilon }{\rm{d}}t + \\ \;\;\;\;\;\;\sqrt {Q(u,H)} \left( {{\gamma _1} + u{\gamma _2}} \right)\sqrt \varepsilon {\xi _t}{\rm{d}}t。\end{array} \right. $

(21)

根据随机平均法，针对0≤H≤H_c，写出上式的漂移与扩散系数分别为^[17]：

$ \begin{array}{l} m\left( H \right) = \frac{1}{{T\left( H \right)}}\int_0^{T(H)} Q (u(t),H)\left[ { - {\beta _1} - {\beta _3}Q(u(t)),} \right.\\ H)]{\rm{d}}t + \frac{1}{{T(H)q}}\int_{ - \infty }^0 {{R_{\xi t\xi t}}} (\tau )。\\ \left\{ {\int_0^{4K(\bar k)} {\frac{{c{n_{t + \tau }} \cdot {\rm{d}}{n_{t + \tau }} \cdot \left[ {b(H){\gamma _2} \cdot s{n_t} + {\gamma _1}} \right]\left[ {b(H){\gamma _2} \cdot s{n_{t + \tau }} + {\gamma _1}} \right]}}{{c{n_t} \cdot d{n_t}}}} \cdot } \right.\\ {\rm{d}}(qt)\} {\rm{d}}\tau 。\end{array} $

(22)

$ \begin{array}{l} {\sigma ^2}(H) = \frac{{b{{(H)}^2}q}}{{T(H)}}\int_{ - \infty }^\infty {{R_{{\xi _t}{\xi _t}}}} (\tau )\left\{ {\int_0^{4K(\bar k)} c {n_t} \cdot {\rm{d}}{n_t} \cdot } \right.\\ c{n_{t + \tau }} \cdot {\rm{d}}{n_{t + \tau }} \cdot \left[ {b(H){\gamma _2} \cdot s{n_t} + {\gamma _1}} \right]\left[ {b(H){\gamma _2} \cdot s{n_{t + \tau }} + } \right.\\ \left. {\left. {{\gamma _1}} \right]{\rm{d}}(qt)} \right\}{\rm{d}}\tau 。\end{array} $

(23)

其中，T(H)为时间区间，其表达式为：

$ T(H) = 2\int_{ - b(H)}^{b(H)} {\frac{{{\rm{d}}u}}{{\sqrt {Q(u,H)} }}} = \frac{{4K(\bar k)}}{q}。$

(24)

上式中，K(k)为第一类完全椭圆积分, 且有：

$ \bar k = \frac{{b{{(H)}^2}}}{2}\sqrt {\frac{{{\alpha _3}}}{H}} ; $

(25)

$ q = \frac{{\sqrt {2H} }}{{b(H)}}; $

(26)

$ b(H) = \sqrt { - \frac{{ - {\alpha _1} + \sqrt {\alpha _1^2 - 4{\alpha _3}H} }}{{{\alpha _3}}}} 。$

(27)

sn、cn、dn为雅克比椭圆函数，且有如下定义：

$ s{n_t} = sn(qt,\bar k); $

(28)

$ c{n_t} = cn(qt,\bar k); $

(29)

$ {\rm{d}}{n_t} = {\rm{d}}n(qt,\bar k); $

(30)

$ s{n_{t + \tau }} = sn(q(t + \tau ),\bar k); $

(31)

$ c{n_{t + \tau }} = \mathit{cn}(q(t + \tau ),\bar k); $

(32)

$ {\rm{d}}{n_{t + \tau }} = dn(q(t + \tau ),\bar k)。$

(33)

系统的能量H满足如下的伊藤随机微分方程：

$ {\rm{d}}H = m(H){\rm{d}}t + \sigma (H){\rm{d}}{W_t}。$

(34)

对于上述伊藤随机微分方程，其响应为扩散过程，对应的转移概率密度p(H, t|H₀, t₀)满足如下的FPK方程：

$ \frac{{\partial p}}{{\partial t}} = m(H)\frac{{\partial p}}{{\partial H}} + \frac{1}{2}{\sigma ^2}(H)\frac{{{\partial ^2}p}}{{\partial {H^2}}}。$

(35)

上式中，(H₀, t₀)表示初始状态。

式(35)可进一步写成如下的算子形式：

$ \frac{{\partial p}}{{\partial t}} - {L^*}p = 0。$

(36)

其中，L^*为椭圆算子的伴随算子，且有：

$ {L^*}p = m(H)\frac{{\partial p}}{{\partial H}} + \frac{1}{2}{\sigma ^2}(H)\frac{{{\partial ^2}p}}{{\partial {H^2}}}。$

(37)

当随机激励输入系统的能量与系统阻尼耗散的能量达到统计平衡时，有∂p/∂t=0，对应的概率密度函数为平稳概率密度函数p_st(H)，其满足如下的平稳FPK方程：

$ {L^*}{p_{st}}(H) = 0。$

(38)

引入尺度函数l(u)与速度函数υ(u)：

$ l(u) = \int_{u0}^u {{\varphi ^{ - 1}}} (U){\rm{d}}U; $

(39)

$ v(u) = \int_{u0}^u {\frac{{\varphi (U)}}{{{\sigma ^2}(U)}}} {\rm{d}}U。$

(40)

其中，$\varphi (u) = \exp \left[ {\int {\frac{{2m(u)}}{{{\sigma ^2}(u)}}} {\rm{d}}u} \right]$。式(38)可进一步写为：

$ {L^*}{p_{st}} = \frac{1}{2}\frac{\partial }{{\partial v}}\left( {\frac{{\partial {p_{st}}}}{{\partial l}}} \right) = 0。$

(41)

求得式(41)的平稳概率密度函数为：

$ {p_{st}}(H) = \frac{{\varphi (H)}}{{{\sigma ^2}(H)}}\left[ {{e_1}l(H) + {e_2}} \right]。$

(42)

其中，参数e₁、e₂由边界条件确定。

对于低强度噪声激励有e₁=0^[13]，则式(42)

可进一步写为：

$ {p_{st}}(H) = \frac{{{e_2}}}{{{\sigma ^2}(H)}}\exp \left( {2\int_0^H {\frac{{m(H)}}{{{\sigma ^2}(H)}}} } \right)。$

(43)

应用传递函数关系p_st(b)=p_st(H) ((d/dH)b)^-1，得到基于变量b的平稳概率密度函数为：

$ {p_{st}}(b) = {p_{st}}(H)\sqrt {\alpha _1^2 - 4{\alpha _3}H} \sqrt {\frac{{{\alpha _1} - \sqrt {\alpha _1^2 - 4{\alpha _3}H} }}{{{\alpha _3}}}} 。$

(44)

对平稳概率密度函数式(44)在某一角度范围[φ₁, φ₂]内积分，可进一步得到b(H)在该角度范围内的概率，从而对船舶的横摇稳定性进行判断。

4 计算结果与分析 4.1 船舶横摇响应

本文针对C11集装箱船展开研究，该船主尺度如表 1所示。

采用文献[9]中所述的方法，在Φ∈[-0.6，0.6] rad、η∈[-10, 10] m、Ψ∈[0，2π] rad范围内数值模拟C11集装箱船的复原力臂GZ(Φ, η, Ψ)，得到不同Φ、η和Ψ时的复原力臂。采用GZ_a(Φ, η, Ψ)近似GZ(Φ, η, Ψ)，得到对应式(2)的各项拟合系数为q₁=2.164 3、q₂=-1.775 3、q₃=-0.106 1。

本文研究采用ITTC海浪谱。对于特征波高3 m、特征波长262 m、航向角150°、航速1.43 m/s的情况，得到对应于式(14)的船舶横摇运动方程为：

$ \left\{ {\begin{array}{*{20}{l}} {\dot {\bar u} = \bar v}\\ {\dot {\bar v} = - 0.2373\bar u + 0.1947{{\bar u}^3} - 0.0269\bar v - }\\ {16.8258{{\bar v}^3} + (0.0013 + 0.0116\bar u){\xi _t}} \end{array}} \right.。$

(45)

取小参量取参量ε=0.1，得到对应于式(15)的无因次横摇运动方程为：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}u}}{{{\rm{d}}{t_\varepsilon }}} = v}\\ {\frac{{{\rm{d}}v}}{{{\rm{d}}{t_\varepsilon }}} = - 23.73u + 19.47{u^3} - 0.1 \times (2.69v + }\\ {\left. {16.8258{v^3}} \right) + \sqrt {0.1} \times (0.4057 + 3.6777u){\xi _{{t_\varepsilon }}}} \end{array}} \right.。$

(46)

采用龙格库塔法数值求解式(45)和式(46)，得到船舶横摇角响应时间历程，如图 1所示。

图 1表明，时间尺度变换前后，船舶横摇响应完全相同，验证了对船舶横摇方程进行尺度变换的正确性。

4.2 概率密度函数的蒙特卡洛验证

蒙特卡洛方法通过数值方法在计算机上模拟实际的响应过程，然后进行统计。该方法在系统可靠度评估、风险评估及结构非线性随机动力响应分析方面具有广泛的应用^[19]。以下采用蒙特卡洛方法验证本文理论计算概率密度函数的正确性。

蒙特卡洛法计算参数如下：随机种子数n₁=25、迭代步数n₂=5 000 000、时间步长Δt=0.01 s，针对4.1中给出的海况，计算关于系统能量H的概率密度函数，并与采用随机平均法计算得到结果进行对比，如图 2所示。

图 2表明，采用随机平均法计算的概率密度函数与用蒙特卡洛法得到的概率密度函数吻合较好，验证了本文理论推导及半解析计算方法的正确性。在第2节中，将随机海况的波面升高近似为CARMA(2，1)过程，从而得到海浪的自相关函数，该函数是能量包线随机平均法中漂移系数和扩散系数的组成部分。由于近似拟合存在一定的误差，因此在种子数足够多的情况下，随机平均法与蒙特卡洛法计算得到的概率密度函数仍会存在一定差别。

4.3 概率密度函数和概率

考虑航速U=1.43 m/s、航向角β=45°、特征波高H=10 m，特征波长L分别为100、150、200、250、300 m，计算船舶机参-强激励横摇响应的概率密度函数和概率。

针对以上工况，由式(4)~(10)，可得到遭遇频率下随机过程η₂(t)的谱密度函数，如图 3(a)所示。根据式(11)~(13)，应用CARMA(2, 1)过程拟合随机过程η₂，得到拟合的谱密度函数，如图 3(b)所示。得到拟合系数c_i(i=1, 2, 3)和γ₁如表 2所示，计算得到不同工况关于b的概率密度函数，如图 4所示。

对图 4中概率密度函数进一步在[0，0.05]、[0.05，0.10]、[0.10，0.15]、[0.15，0.20]、[0.20，0.25]、[0.25，0.30]及[0.30，0.6] rad的区间内积分，得到b在各区间的概率，结果如表 3示。表 3表明，当遭遇波长(L/cosβ)远离船长时，b在小角度区间的概率较大；当遭遇波长接近船长时，b在大角度区间的概率增加，该结论与已有研究工作的结论相同^[1]。以b在[0.25，0.30] rad区间的概率为例，L=200 m(遭遇波长为282.8 m)时为0.018 1，是L=150 m(遭遇波长为212.1 m)时的3.3倍，是L=100 m(遭遇波长为141.4 m)时的16.5倍。

5 结语

本文基于能量包线随机平均法，以C11船为例，半解析的研究了斜浪中船舶随机参-强激励横摇响应的概率密度函数及响应概率。结果表明，当遭遇频率远离船长时，b(H)在大角度区间的概率较小；当遭遇频率接近船长时，b(H)在大角度区间的概率增加。由b(H)可定量分析斜浪中船舶随机参-强激励横摇运动概率。