随机激励下非线性时滞船舶横摇模型的随机分岔

Random Bifurcation of Nonlinear Time-Delay Ship Rolling Motion under Random Excitation

摘要: 本文研究了同时带有随机激励与时滞反馈控制的船舶横摇动力学模型。首先，分析了系统发生Hopf分岔的临界条件；其次，借助中心流形理论和最大Lyapunov指数，考察了系统的局部稳定性及随机D分岔的产生条件，并基于奇异边界理论对系统的全局稳定性进行了讨论；最后，通过数值模拟揭示了噪声强度与时滞系数对系统动力学行为的影响机制，从而验证了相关理论结果的正确性。

Abstract: This paper investigates a ship rolling dynamics model with both random excitation and time-delay feedback control. First, the critical conditions for the occurrence of Hopf bifurcations are analyzed. Second, using central manifold theory and the maximum Lyapunov exponent, the local stability of the system and the conditions for the occurrence of random D-bifurcation are investigated. The global stability of the system is discussed based on singular boundary theory. Finally, numerical simulations reveal the mechanisms by which noise intensity and time-delay coefficients influence the system’s dynamic behavior, validating the theoretical results.

1. 引言

船舶是海上运输的核心载体，其载货能力与航行安全性一直备受关注[1]。在航行过程中，船舶常受到风浪等环境因素的干扰，引发复杂的非线性运动[2]-[4]。一般而言，船舶运动可分解为艏摇、纵摇和横摇三种角运动，以及升沉、横荡和纵荡三种平移运动。其中，由风浪引发的横摇运动尤为显著，强烈的横摇不仅会影响船体稳定性、中断正常航行、造成船员不适，还可能严重威胁船舶与人员安全，甚至引发倾覆事故[5]。因此，深入研究船舶横摇运动的动力学特性、稳定性[6]及其控制方法，对提升航行安全具有重要实际意义。

近年来，国内外不少学者都对船舶横摇运动进行了研究。Kannan等人[7]提出了一种高效的勒让德小波谱方法，用于研究船舶横摇运动模型中的非线性阻尼系数，并用数值算例验证了所提小波方法的有效性和适用性。朱等人[8]提出了一个具有随机振幅的级数展开模型，基于该模型，研究了两个相关波浪力作用下的船舶横摇运动，重点讨论了两个随机波浪力相关性的影响。Jensen等人[9]在一阶可靠性方法的框架内，考虑船舶横摇运动的极值预测，重点关注泊松假设对峰值统计的适用性。

实际上，海洋环境天然具有显著的不规则性和随机特征，因此研究船舶在随机波浪激励下的动态响应行为更具有实际意义。近年来，随机动力系统理论受到广泛关注，众多学者对其进行了深入研究，涵盖了随机稳定性[10]、随机分岔[11]、随机共振[12]、反馈控制[13]等多个重要方向。随着该领域的不断发展，一系列近似与数值方法被广泛应用于非线性随机系统的分析中，如随机平均法[14]，路径积分法[15]等。

具体到随机激励下船舶横摇运动动力学这个课题，陈等人[16]将随机波浪激励描述为谐波激励和高斯白噪声激励的组合，利用有限元法和Crank-Nicolson方法对相应的瞬态概率密度函数进行了数值求解。魏等人[17]针对一类船舶在非光滑扰动和随机激励即风浪作用下在纵波中横摇运动的随机分岔与混沌，通过随机平均法得到了横摇运动的稳态概率密度函数，并对其定性分析了随机P分岔。周等人[18]提出了一种新的极值统计策略，通过极少的模型试验观测数据来预报船舶参数横摇角的极值。Narayanan等人[19]研究了包括非线性在内的随机波浪激励下光滑船舶和撞击船舶的横摇动力学，采用FD方法求解四维马尔可夫系统相应的FPK方程来研究船舶的随机响应。

目前，在船舶横摇动力学研究中，同时纳入时滞效应与随机激励因素的文献尚不多见。本文在现有非线性船舶横摇运动模型的基础上，将实际航行中遇到的随机干扰近似描述为乘性噪声，并结合时滞反馈控制，构建了一类兼具随机激励与时滞特性的船舶横摇模型。文章结构安排如下：第二部分建立了带有时滞的非线性随机船舶横摇模型；第三部分对该系统进行理论化简，推导出相应的Itô随机微分方程，并从概率角度分析了系统的随机分岔行为；第四部分通过数值模拟对理论结果进行了验证。

2. 数学模型与稳定性分析

考虑如下单自由度非线性船舶横摇运动微分方程：

$\ddot{\phi }+2\mu \dot{\phi }+{\mu }_1{\dot{\phi }}^3+{\omega }_0^2\phi +{\alpha }_1{\phi }^3+{\alpha }_2{\phi }^5=F\cos \Lambda t$ (1)

其中$\phi$ 表示船舶横摇角，$\mu$ 和${\mu }_1$ 为无量纲阻尼扭矩系数，${\omega }_0$ 为船舶横摇运动固有圆频率，${\alpha }_1$ 和${\alpha }_2$ 为无量纲恢复扭矩系数，$F$ 和$\Lambda$ 分别为波浪激励幅值和圆频率。

船舶在波浪中航行时，所受的波浪激励力不是确定的，而是随机的。对于船舶横摇运动，当船舶横摇角较大时，波浪冲击船体的位置和角度随横摇角改变，导致波浪力矩的幅度依赖于当前横摇状态，可以将乘性噪声作为随机因素引入非线性船舶横摇模型。此外，还考虑引入控制函数，该函数表示由具有时间延迟的延时状态反馈产生的控制信号[20] [21]。就一类主动减摇系统的工作原理而言，其抗横摇力矩由偏心配重的运动引起，该重物的运动由一个受延时控制器控制的执行器完成，延时控制器接收来自安装在船舶上的传感器的输入信号。因此，在上述模型(1)中分别引入乘性噪声激励和时滞反馈控制，具体模型如下：

$\ddot{\phi }+2\mu \dot{\phi }+{\mu }_1{\dot{\phi }}^3+{\omega }_0\phi +{\alpha }_1{\phi }^3+{\alpha }_2{\phi }^5=F\cos \Lambda t+\delta \phi \xi \left(t\right)+U$ (2)

其中$U={\beta }_1\phi \left(t-\tau \right)+{\beta }_2\dot{\phi }\left(t-\tau \right)$ ，${\beta }_1$ 和${\beta }_2$ 分别表示时滞位移反馈和速度反馈的控制增益，$\tau$ 为时滞，$\xi \left(t\right)$ 为均值为零的高斯白噪声，其强度用$\delta$ 表示。

船舶横摇运动模型(2)的线性部分为：

$\ddot{\phi }+2\mu \dot{\phi }+{\omega }_0\phi -{\beta }_1\phi \left(t-\tau \right)-{\beta }_2\dot{\phi }\left(t-\tau \right)=0$ (3)

系统(2)在平衡点$\left(0,0\right)$ 处的超越特征方程如下所示：

$\Delta =$ ${\lambda }^2+2\mu \lambda +{\omega }_0^2-{\beta }_1{\text{e}}^{-\lambda \tau }-{\beta }_2\lambda {\text{e}}^{-\lambda \tau }=0$ (4)

当$\tau =0$ 时，特征方程的解为：

${\lambda }_{1,2}=\frac{-\left(2\mu -{\beta }_2\right)\pm \sqrt{{\left(2\mu -{\beta }_2\right)}^2-4\left({\omega }_0^2-{\beta }_1\right)}}{2}$

当且仅当$\mu >\frac{{\beta }_2}{2}$ 时，系统平衡点$\left(0,0\right)$ 处是渐近稳定的。当$\tau \ne 0$ 时，设存在一对共轭特征根${\lambda }_{1,2}=v\pm i\omega$ ，把${\lambda }_1=v+i\omega$ 代入到特征方程(4)中，得到：$\Delta =H+iS=0$

$H=v^2-{\omega }^2+2\mu v+{\omega }_0^2-{\text{e}}^{-v\tau }\left({\beta }_1\cos \omega \tau +{\beta }_{2}v\cos \omega \tau +{\beta }_2\omega \sin \omega \tau \right)=0$ ，

$S=2\omega \left(v+\mu \right)+{\text{e}}^{-v\tau }\left({\beta }_1\sin \omega \tau -{\beta }_2\omega \cos \omega \tau +{\beta }_{2}v\sin \omega \tau \right)=0$ ，

通过使$v=0$ 可简化为：

$\{\begin{array}{l}-{\omega }^2+{\omega }_0^2={\beta }_1\cos \omega \tau +{\beta }_2\omega \sin \omega \tau ,\hfill \\ 2\mu \omega =-{\beta }_1\sin \omega \tau +{\beta }_2\omega \cos \omega \tau .\hfill \end{array}$ (5)

对两个等式的两边求平方并把结果相加可以得到：

${\omega }^4-\left(2{\omega }_0^2-4{\mu }^2+{\beta }_2^2\right){\omega }^2+{\omega }_0^4-{\beta }_1^2=0$ (6)

求解方程(6)得：

$\omega =\sqrt{\frac{\left(2{\omega }_0^2-4{\mu }^2+{\beta }_2^2\right)\pm \sqrt{{\left(2{\omega }_0^2-4{\mu }^2+{\beta }_2^2\right)}^2-4\left({\omega }_0^4-{\beta }_1^2\right)}}{2}}$

通过联立式(5)，计算临界时滞为：

${\tau }_c=\frac{1}{\omega }\arctan \frac{\omega \left(2\mu {\beta }_1-{\omega }_0^2{\beta }_2+\omega {\beta }_2\right)}{{\omega }_0^2{\beta }_1+{\omega }^2\left({\beta }_1-2\mu {\beta }_2\right)}$

对特征方程两端作关于$\tau$ 的微分得：

$\frac{\text{d}\lambda }{\text{d}\tau }=\frac{-\lambda {\text{e}}^{-\lambda \tau }\left({\beta }_1+{\beta }_2\lambda \right)}{{\text{e}}^{-\lambda \tau }\left(\tau {\beta }_1+\tau {\beta }_2\lambda -{\beta }_2\right)+2\left(\mu +\lambda \right)}$ ，${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=\frac{{\text{e}}^{-\lambda \tau }\left(\tau {\beta }_1+\tau {\beta }_2\lambda -{\beta }_2\right)+2\left(\mu +\lambda \right)}{-\lambda {\text{e}}^{-\lambda \tau }\left({\beta }_1+{\beta }_2\lambda \right)}$

根据Hopf分岔定理，若$\mathrm{Re}\left[{{\left(\text{d}\lambda /\text{d}\tau \right)}^{-1}|}_{\lambda =i\omega }\right]\ne 0$ ，则系统(2)存在周期解。

3. 随机分岔分析

本节，我们将使用随机中心流形定理[22]简化船舶横摇运动模型(2)，并通过奇异边界理论和最大 Lyapunov指数分析系统的随机稳定性[23]。

定义双线性算子，令${\Phi }_k\in C\left(\left[-\tau ,0\right],R^2\right)$ ，${\Psi }_j\in \widehat{C}\left(\left[0,\tau \right],R^2\right)$ ，$j,k=1,2$ ，则有：

$\left({\Psi }_j\left(s\right),{\Phi }_k\left(\theta \right)\right)=\left({\Psi }_j\left(0\right),{\Phi }_k\left(0\right)\right)+\left({\beta }_1+{\beta }_2\right){\displaystyle {\int }_{-\tau }^0\left({\Psi }_j\left(\zeta +\tau \right),{\Phi }_k\left(\zeta \right)\right)\text{d}\zeta }$ . (7)

若特征方程(4)有一对纯虚根$\lambda =\pm i\omega$ 时，其所有特征值的解空间$C$ 可以表示为子空间$P$ 和子空间$Q$ 的直和，即$C=P\oplus Q$ 。其中$P$ 表示由特征值$\lambda =\pm i\omega$ 生成的二维子空间，$Q$ 表示由$\Delta =0$ 的其余特征值生成的无限维子空间。系统(3)的子空间$P\in C$ 的基函数$\Phi \left(\theta \right)$ 和对应的伴随子空间$\widehat{P}\in \widehat{C}$ 的基函数$\Psi \left(s\right)$ 为：

$\Phi \left(\theta \right)=\left(\begin{array}{cc}\cos \omega \theta & -\sin \omega \theta \\ \omega \sin \omega \theta & \omega \cos \omega \theta \end{array}\right)=\left({\Phi }_1,{\Phi }_2\right),$ (8)

$\Psi \left(s\right)=\left(\begin{array}{cc}\cos \omega s& -\sin \omega s\\ \omega \sin \omega s& \omega \cos \omega s\end{array}\right)={\left({\Psi }_1,{\Psi }_2\right)}^{\text{T}}$ . (9)

定义内积矩阵$\left(\Psi \left(s\right),\Phi \left(\theta \right)\right)=\left({\Psi }_j,{\Phi }_k\right)$ ，$j,k=1,2$ ，得到：

$\left(\Psi \left(s\right),\Phi \left(\theta \right)\right)=\left(\begin{array}{cc}\left({\Psi }_1,{\Phi }_1\right)& \left({\Psi }_1,{\Phi }_2\right)\\ \left({\Psi }_2,{\Phi }_1\right)& \left({\Psi }_2,{\Phi }_2\right)\end{array}\right)$ ， (10)

其中

$\left({\Psi }_1,{\Phi }_1\right)=\cos \omega s\cos \omega \theta +{\omega }^2\sin \omega s\sin \omega \theta ,\left({\Psi }_1,{\Phi }_2\right)=-\cos \omega s\sin \omega \theta +{\omega }^2\sin \omega s\cos \omega \theta$

$\left({\Psi }_2,{\Phi }_1\right)=-\sin \omega s\cos \omega \theta +{\omega }^2\cos \omega s\sin \omega \theta ,\left({\Psi }_2,{\Phi }_2\right)=\sin \omega s\cos \omega \theta +{\omega }^2\cos \omega s\cos \omega \theta$

把内积矩阵$\left(\Psi \left(s\right),\Phi \left(\theta \right)\right)$ 代入双线性算子(7)中，可以得到非奇异矩阵：

${\left(\Psi ,\Phi \right)}_{nsg}=\left(\begin{array}{cc}{\psi }_{11}& {\psi }_{12}\\ {\psi }_{21}& {\psi }_{22}\end{array}\right)$ ， (11)

其中

${\psi }_{11}=1+\frac{{\beta }_1+{\beta }_2}{2}\left[\left(\frac{1}{\omega }-\omega \right)\sin \omega \tau +\tau \left({\omega }^2+1\right)\cos \omega \tau \right],{\psi }_{12}=\frac{{\beta }_1+{\beta }_2}{2}\tau \left({\omega }^2+1\right)\sin \omega \tau ,$

${\psi }_{21}=\frac{{\beta }_1+{\beta }_2}{2}\tau \left(-{\omega }^2-1\right)\sin \omega \tau =-{\psi }_{12},{\psi }_{22}={\omega }^2+\frac{{\beta }_1+{\beta }_2}{2}\left[\left(\omega -\frac{1}{\omega }\right)\sin \omega \tau +\tau \left({\omega }^2+1\right)\cos \omega \tau \right].$

通过使用新的基函数$\overline{\Psi }\left(s\right)={\left(\Psi ,\Phi \right)}_{nsg}^{-1}\Psi \left(s\right)$ ，把基函数$\Psi \left(s\right)$ 归一化，从而，

$\overline{\Psi }\left(s\right)=\left(\begin{array}{cc}{\overline{\Psi }}_{11}\left(s\right)& {\overline{\Psi }}_{12}\left(s\right)\\ {\overline{\Psi }}_{21}\left(s\right)& {\overline{\Psi }}_{22}\left(s\right)\end{array}\right)$ ， (12)

其中

$\begin{array}{l}{\overline{\Psi }}_{11}\left(s\right)={\left({\psi }_{11}{\psi }_{22}+{\psi }_{12}^2\right)}^{-1}\left({\psi }_{22}\cos \omega s+{\psi }_{12}\sin \omega s\right),\\ {\overline{\Psi }}_{12}\left(s\right)=\omega {\left({\psi }_{11}{\psi }_{22}+{\psi }_{12}^2\right)}^{-1}\left({\psi }_{22}\sin \omega s-{\psi }_{12}\cos \omega s\right),\\ {\overline{\Psi }}_{21}\left(s\right)=-{\left({\psi }_{11}{\psi }_{22}+{\psi }_{12}^2\right)}^{-1}\left({\psi }_{21}\cos \omega s+{\psi }_{11}\sin \omega s\right),\\ {\overline{\Psi }}_{22}\left(s\right)=-\omega {\left({\psi }_{11}{\psi }_{22}+{\psi }_{12}^2\right)}^{-1}\left({\psi }_{21}\sin \omega s-{\psi }_{11}\cos \omega s\right).\end{array}$

设$y\left(t\right)=\left(\overline{\Psi }\left(s\right),{\varphi }^P\left(\theta \right)\right)\in R^2$ ，若方程(2)有唯一解，且该唯一解可以表示为${\phi }_t\left(\varphi \left(\theta \right),\mu \right)={\phi }_t^P\left(\varphi \left(\theta \right),\mu \right)+{\phi }_t^Q\left(\varphi \left(\theta \right),\mu \right)$ ，在这当中$\varphi \left(\theta \right)={\varphi }^P\left(\theta \right)+{\varphi }^Q\left(\theta \right)$ ，${\varphi }^P\left(\theta \right)\in P$ ，${\varphi }^Q\left(\theta \right)\in Q$ 。${\phi }_t^P\left(\varphi \left(\theta \right),\mu \right)$ 和${\phi }_t^Q\left(\varphi \left(\theta \right),\mu \right)$ 分别是${\phi }_t\left(\varphi \left(\theta \right),\mu \right)$ 在半群$J\left(t,\tau \right)$ 和无穷小生成元$A\left(\theta ,\tau \right)$ 下投影到子空间$P$ ，$Q$ 上的解，它们保持不变。纯虚特征值$\lambda =\pm i\omega$ 张成的特征子空间$P$ ，与中心流形$M_{\mu }\subseteq C\left(\left[-\tau ,0\right],R^2\right)$ 相切。根据变量$x_t^P\left(\theta \right)=\Phi \left(\theta \right)y\left(t\right)+x_t^Q\left(\theta \right)$ ，$\theta =-\tau$ 的一阶近似为：

$\{\begin{array}{l}\phi \left(t-\tau \right)=y_1\left(t\right)\cos \omega \tau +y_2\left(t\right)\sin \omega \tau \hfill \\ \dot{\phi }\left(t-\tau \right)=-\omega y_1\left(t\right)\sin \omega \tau +\omega y_2\left(t\right)\cos \omega \tau \hfill \end{array}$ (13)

对于系统(2)，若满足Hopf分岔的条件，相应的随机延迟微分系统中心流形形式为：

$\{\begin{array}{l}\begin{array}{c}{\dot{y}}_1\left(t\right)=-\omega y_2\left(t\right)+{\overline{\Psi }}_{12}[-{\mu }_1{y}_2^3\left(t\right)-{\alpha }_1{y}_1^3\left(t\right)-{\alpha }_2{y}_1^5\left(t\right)+F\cos \Lambda t+\sigma y_1\left(t\right)\xi \left(t\right)\\ +{\beta }_2\left(-y_1\left(t\right)\omega \sin \omega \tau +y_2\left(t\right)\omega \cos \omega \tau \right)+{\beta }_1\left(y_1\left(t\right)\cos \omega \tau +y_2\left(t\right)\sin \omega \tau \right)]\end{array}\hfill \\ \begin{array}{c}{\dot{y}}_2\left(t\right)=\omega y_1\left(t\right)+{\Psi }_{22}[-{\mu }_1{y}_2^3\left(t\right)-{\alpha }_1{y}_1^3\left(t\right)-{\alpha }_2{y}_1^5\left(t\right)+F\cos \Lambda t+\sigma y_1\left(t\right)\xi \left(t\right)\\ +{\beta }_2\left(-y_1\left(t\right)\omega \sin \omega \tau +y_2\left(t\right)\omega \cos \omega \tau \right)+{\beta }_1\left(y_1\left(t\right)\cos \omega \tau +y_2\left(t\right)\sin \omega \tau \right)]\end{array}\hfill \end{array}$， (14)

其中

${\overline{\Psi }}_{12}={\overline{\Psi }}_{12}\left(0\right)=\frac{-\omega {\psi }_{12}}{{\psi }_{11}{\psi }_{22}+{\psi }_{12}^2},{\overline{\Psi }}_{22}={\overline{\Psi }}_{22}\left(0\right)=\frac{\omega {\psi }_{11}}{{\psi }_{11}{\psi }_{22}+{\psi }_{12}^2}$

引入变量代换，即$\Gamma :\left(y_1,y_2\right)\to \left(m_1,m_2\right)$ ，如下所示为：

$\Gamma :\{\begin{array}{l}{m}_1=y_1-\frac{{\overline{\Psi }}_{12}}{{\overline{\Psi }}_{22}}{y}_2,\hfill \\ m_2=y_2+\frac{{\overline{\Psi }}_{12}}{{\overline{\Psi }}_{22}}{y}_1.\hfill \end{array}$ (15)

再令$l=\frac{{\overline{\Psi }}_{12}}{{\overline{\Psi }}_{22}}$ ，得到$\Gamma$ 的逆变换：

${\Gamma }^{-1}:\{\begin{array}{l}{y}_1=\frac{1}{1+l^2}{m}_1+\frac{l}{1+l^2}{m}_2,\hfill \\ y_2=-\frac{l}{1+l^2}{m}_1+\frac{1}{1+l^2}{m}_2.\hfill \end{array}$ (16)

把(16)式代入(14)式，化简得到中心流形：

$\{\begin{array}{l}{\dot{m}}_1=-\omega m_2\hfill \\ \begin{array}{c}{\dot{m}}_2=\omega m_1+p_0(p_1{p}_9{m}_1+p_1{p}_8{m}_2+p_3{p}_7{m}_1^3+p_3{p}_6{m}_2^3+p_3{p}_4{m}_1^2{m}_2+p_3{p}_5{m}_1{m}_2^2\\ \begin{array}{l}+p_{10}{m}_1^5+p_{10}{p}_{11}{m}_1^4{m}_2+p_{10}{p}_{12}{m}_1^3{m}_2^2+p_{10}{p}_{13}{m}_1^2{m}_2^3+p_{10}{p}_{14}{m}_1{m}_2^4+p_{10}{p}_{15}{m}_2^5)\\ +p_{0}F\cos \Lambda t+\delta \left(p_0{p}_1{m}_1+p_0{p}_2{m}_2\right)\xi \left(t\right),\end{array}\end{array}\hfill \end{array}$ (17)

其中

$\begin{array}{l}{p}_0={\overline{\Psi }}_{12}l+{\overline{\Psi }}_{22},p_1=\frac{1}{1+l^2},p_2=\frac{l}{1+l^2},p_3=\frac{1}{{\left(1+l^2\right)}^3},p_4=-{\mu }_1-l^3{\alpha }_1,p_5=l^3{\mu }_1-{\alpha }_1,\\ p_6=-3l^2{\mu }_1-3l{\alpha }_1,p_7=3l{\mu }_1-3l^2{\alpha }_1,p_8=\left({\beta }_2-l\omega {\beta }_1\right)\cos \omega \tau -\left(l{\beta }_2+{\beta }_1\omega \right)\sin \omega \tau ,\\ p_9=\left(l{\beta }_2+\omega {\beta }_1\right)\cos \omega \tau +\left({\beta }_2-l{\beta }_1\omega \right)\sin \omega \tau ,p_{10}=-\frac{{\alpha }_2}{{\left(1+l^2\right)}^5},p_{11}=5l,p_{12}=10l^2,\\ p_{13}=10l^3,p_{14}=5l^4,p_{15}=l^5.\end{array}$

引入变量代换$m_1=R\left(t\right)\cos \vartheta \left(t\right),m_2=R\left(t\right)\sin \vartheta \left(t\right),\vartheta \left(t\right)=\omega t+\phi \left(t\right)$ ，其中，$R\left(t\right)$ 和$\phi \left(t\right)$ 分别为解的振幅和相位。通过极坐标变换，可以得到关于$R\left(t\right)$ 和$\phi \left(t\right)$ 的随机微分方程，如下所示：

$\{\begin{array}{l}\dot{R}\left(t\right)=f_1\left(R,\vartheta \right)+g_1\left(R,\vartheta \right)\xi \left(t\right),\hfill \\ \dot{\phi }\left(t\right)=f_2\left(R,\vartheta \right)+g_2\left(R,\vartheta \right)\xi \left(t\right),\hfill \end{array}$ (18)

其中

$\begin{array}{c}{f}_1\left(R,\vartheta \right)=p_0(p_1{p}_{9}R\sin \vartheta \cos \vartheta +p_1{p}_{8}R{\sin }^2\vartheta +p_3{p}_7{R}^3\sin \vartheta {\cos }^3\vartheta +p_3{p}_6{R}^3{\sin }^4\vartheta \\ +p_{10}{p}_{11}{R}^5{\sin }^2\vartheta {\cos }^4\vartheta +p_{10}{p}_{12}{R}^5{\sin }^3\vartheta {\cos }^3\vartheta +p_{10}{p}_{13}{R}^5{\sin }^4\vartheta {\cos }^2\vartheta \\ +p_3{p}_4{R}^3{\sin }^2\vartheta {\cos }^2\vartheta +p_3{p}_5{R}^3{\sin }^3\vartheta \cos \vartheta +p_{10}{R}^5\sin \vartheta {\cos }^5\vartheta \\ +p_{10}{p}_{14}{R}^5{\sin }^5\vartheta \cos \vartheta +p_{10}{p}_{15}{R}^5{\sin }^6\vartheta +F\sin \vartheta \cos \Lambda t),\end{array}$

$\begin{array}{c}{f}_2\left(R,\vartheta \right)=p_0(p_1{p}_9{\cos }^2\vartheta +p_1{p}_8\sin \vartheta \cos \vartheta +p_3{p}_7{R}^2{\cos }^4\vartheta +p_3{p}_6{R}^2{\sin }^3\vartheta \cos \vartheta \\ +p_{10}{p}_{11}{R}^4\sin \vartheta {\cos }^5\vartheta +p_{10}{p}_{12}{R}^4{\sin }^2\vartheta {\cos }^4\vartheta +p_{10}{p}_{13}{R}^4{\sin }^3\vartheta {\cos }^3\vartheta \\ +p_{10}{p}_{14}{R}^4{\sin }^4\vartheta {\cos }^2\vartheta +p_{10}{p}_{15}{R}^4{\sin }^5\vartheta \cos \vartheta +F{R}^{-1}\cos \vartheta \cos \Lambda t\\ +p_3{p}_4{R}^2\sin \vartheta {\cos }^3\vartheta +p_3{p}_5{R}^2{\sin }^2\vartheta {\cos }^2\vartheta +p_{10}{R}^4{\cos }^6\vartheta ),\end{array}$

$g_1\left(R,\vartheta \right)=\delta \left(p_0{p}_{1}R\cos \vartheta \sin \vartheta +p_0{p}_{2}R{\sin }^2\vartheta \right),$

$g_2\left(R,\vartheta \right)=\delta \left(p_0{p}_1{\cos }^2\vartheta +p_0{p}_2\sin \vartheta \cos \vartheta \right).$

一般来说，船舶横摇模型中的非线性项较弱，但非线性效应不可忽略，根据随机平均法[24]，第一步，在随机时滞船舶横摇模型中，我们使用了高斯白噪声作为随机激励项，使系统的响应成为Markov扩散过程；第二步，船舶横摇运动固有频率与波浪激励频率不同，作时间平均，消除快变量，从而降低系统维数。此时得到受一维概率密度函数主导的Markov过程，该过程由如下Itô方程控制：

$\text{d}R=m\left(R\right)\text{d}t+\sigma \left(R\right)\text{d}W\left(t\right)$ (19)

$m\left(R\right)$ 和$\sigma \left(R\right)$ 分别为式(18)的漂移系数和扩散系数，$W\left(t\right)$ 为标准Wiener过程，其中：

$m\left(R\right)=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }\left[f_1\left(R,\vartheta \right)+\frac{1}{2}\frac{\partial g_1\left(R,\vartheta \right)}{\partial R}{g}_1\left(R,\vartheta \right)+\frac{1}{2}\frac{\partial g_1\left(R,\vartheta \right)}{\partial \vartheta }{g}_2\left(R,\vartheta \right)\right]\text{d}\vartheta }$

${\sigma }^2\left(R\right)=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }{\left[g_1\left(R,\vartheta \right)\right]}^2\text{d}\vartheta }$

把$f_1\left(R,\vartheta \right)$ ,$g_1\left(R,\vartheta \right)$ ,$g_2\left(R,\vartheta \right)$ 代入，计算三角函数积分时注意其正交性，可以大大简化计算，得到：

$m\left(R\right)=\frac{1}{16}{u}_{1}R+\frac{1}{8}{u}_2{R}^3+\frac{1}{16}{u}_3{R}^5$ ，$\sigma \left(R\right)=\sqrt{\frac{1}{8}{u}_4{R}^2}$ (20)

$u_1=8p_0{p}_1{p}_8+3{\delta }^2{p}_0^2{p}_1^2+5{\delta }^2{p}_0^2{p}_2^2,u_2=3p_0{p}_3{p}_6+p_0{p}_3{p}_4$

$u_3=5p_0{p}_{10}{p}_{15}+p_0{p}_{10}{p}_{11}+p_0{p}_{10}{p}_{13},u_4={\delta }^2{p}_0^2{p}_1^2+3{\delta }^2{p}_0^2{p}_2^2$

接下来，本文将对噪声激励下采用时滞位移反馈控制的船舶横摇模型的随机动力学行为作进一步研究。根据Lyapunov指数的定义，可推导出其近似表达式：

$\rho =\underset{t\to \infty }{\lim }\frac{1}{t}\ln \left(\Vert R\left(t\right)\Vert \right)=\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t\left(\frac{1}{16}{u}_1-\frac{1}{16}{u}_4\right)\text{d}t}=\frac{u_1-u_4}{16}$

根据最大Lyapunov指数法，可以得出：当$\rho <0$ ，即$u_1<u_4$ 时，系统(2)的平凡解局部渐近稳定；当$\rho >0$ ，即$u_1>u_4$ 时，系统(2)的平凡解，是局部渐近不稳定的；当$\rho =0$ ，即$u_1=u_4$ 时，系统处于临界状态，可能会发生随机分岔。Lyapunov指数仅能用于判断船舶横摇模型在平衡点附近的局部稳定性，无法反映系统在全局范围内的稳定性。对于随机系统的全局稳定性分析，需借助奇异边界理论。通过随机平均法，可将船舶横摇模型依概率收敛于一维扩散过程。在该过程中，系统在边界处的行为决定了其全局动力学特性，即扩散过程的边界性态可用于推断系统的全局概率稳定性。这一行为由漂移系数$m\left(R\right)$ 和扩散系数$\sigma \left(R\right)$ 共同决定。当$R\to 0$ 时，$m\left(R\right)=0$ ，$\sigma \left(R\right)=0$ 。当$R\to \infty$ 时，$m\left(R\right)=\infty$ 。根据奇异边界理论的分类可知，左边界$R=0$ 为第一类奇异边界，右边界$R=+\infty$ 为第二类奇异边界。

对于左边界$R=0$ ，可以得到对应的扩散系数$\gamma =2$ ，漂移系数$\varsigma =1$ ，和特征值$c=\frac{u_1}{u_4}$ ：(i) 若$c<1$ ，即$u_1<u_4$ ，那么边界$R=0$ 是吸引自然的，这时船舶横摇模型平凡解是随机稳定的；(ii) 若$c=1$ ，即$u_1=u_4$ ，那么边界$R=0$ 是严格自然的，这时系统处于临界状态；(iii) 若$c>1$ ，即$u_1>u_4$ ，那么边界$R=0$ 是排斥自然的，这时船舶横摇模型的平凡解是随机不稳定的。

对于右边界$R=+\infty$ ，可以得到对应的扩散系数$\gamma =2$ ，漂移系数$\varsigma =5$ ，和特征值$c=-\frac{u_3}{u_4}$ ：(i) 若$c<-1$ ，即$u_3>u_4$ ，那么边界$R=+\infty$ 为越出边界；(ii) 若$c=-1$ ，即$u_3=u_4$ ，那么边界$R=+\infty$ 为规则边界；(iii) 若$c>-1$ ，即$u_3<u_4$ ，那么边界$R=+\infty$ 为进入边界。

根据三指数法，当右边界$R=\infty$ 为入口边界时，系统的状态仅依赖于左边界$R=0$ 的情况。通过前面的分析，当$c=-\frac{u_3}{u_4}>-1$ 时，边界$R=+\infty$ 为进入边界，此时只需考虑左边界的情况，若左边界$R=0$ 是第一类奇异边界，平稳概率密度函数可以表示为：

$p\left(R\right)=\frac{C}{u_4{R}^2}\exp \left[{\displaystyle \int \frac{2\left\{\frac{1}{16}{u}_{1}R+\frac{1}{8}{u}_2{R}^3+\frac{1}{16}{u}_3{R}^5\right\}}{\frac{1}{8}{u}_4{R}^2}}\right],$ (21)

其中$C$ 是归一化常数，且${\displaystyle {\int }_0^{+\infty }p\left(R\right)\text{d}R=1}$ 。将$R^2=m_1^2+m_2^2$ 代入式(21)即可得到联合概率密度函数表达式。

由于$\gamma =2,\varsigma =1$ ，平稳概率密度函数可以简化为：

$p\left(R\right)=O\left(R^{c-\gamma }\right)=O\left\{R^{\frac{u_1}{u_4}-2}\right\}$

根据文献[25]，船舶横摇模型在$u_1=u_4$ 时发生随机D分岔，在$\frac{u_1}{u_4}=2$ 发生随机P分岔。

4. 数值模拟

通过平稳概率密度函数图形的变化，可以判断是否发生P分岔，本节，我们将进行数值模拟实验分析平稳概率密度函数形状变化。模型参数选取如下[23] [26]：$\mu =0.21$ ，${\mu }_1=0.01$ ，${\omega }_0=1$ ，${\alpha }_1=-1$ ，${\alpha }_2=-0.5$ ，${\beta }_2=0.1$ ，将分别分析噪声强度$\delta$ 和时滞控制增益系数${\beta }_1$ 对平稳概率密度函数形状的影响。

首先，固定时滞系数${\beta }_1=1.000$ ，改变噪声强度$\delta$ 的大小，观察随机船舶横摇模型的分岔行为，实验结果如图1~图3。

(a) (b)

(c)

Figure 1. When noise intensity $\delta =0.434$ ,$\frac{u_1}{u_4}=1.97<2$ : (a) Joint probability density function of the ship roll model; (b) Cross-section of figure (a); (c) Steady probability density function

图1. 噪声强度$\delta =0.434$ ，$\frac{u_1}{u_4}=1.97<2$ ：(a) 船舶横摇模型的联合概率密度函数图；(b) 图(a)的截面图；(c) 平稳概率密度函数图

(a) (b)

(c)

Figure 2. When noise intensity $\delta =0.452,\frac{u_1}{u_4}=2.00$ : (a) Joint probability density function of the ship roll model; (b) Cross-section of figure (a); (c) Steady probability density function

图2. 噪声强度$\delta =0.452,\frac{u_1}{u_4}=2.00$ ：(a) 船舶横摇模型的联合概率密度函数图；(b) 图(a)的截面图；(c) 平稳概率密度函数图

(a) (b)

(c)

Figure 3. When noise intensity $\delta =0.513,\frac{u_1}{u_4}=2.19>2$ : (a) Joint probability density function of the ship roll model; (b) Cross-section of figure (a); (c) Steady probability density function

图3. 噪声强度$\delta =0.513,\frac{u_1}{u_4}=2.19>2$ ：(a) 船舶横摇模型的联合概率密度函数图；(b) 图(a)的截面图；(c) 平稳概率密度函数图

(a) (b)

(c)

Figure 4. When time delay coefficient${\beta }_1=0.998,\frac{u_1}{u_4}=1.93<2$ : (a) Joint probability density function of the ship roll model; (b) Cross-section of figure (a); (c) Steady probability density function

图4. 时滞系数${\beta }_1=0.998,\frac{u_1}{u_4}=1.93<2$ ：(a) 船舶横摇模型的联合概率密度函数图；(b) 图(a)的截面图；(c) 平稳概率密度函数图

(a) (b)

(c)

Figure 5. When time delay coefficient ${\beta }_1=1.0012,\frac{u_1}{u_4}=2.00$ : (a) Joint probability density function of the ship roll model; (b) Cross-section of figure (a); (c) Steady probability density function

图5. 时滞系数${\beta }_1=1.0012,\frac{u_1}{u_4}=2.00$ ：(a) 船舶横摇模型的联合概率密度函数图；(b) 图(a)的截面图；(c) 平稳概率密度函数图

(a) (b)

(c)

Figure 6. When time delay coefficient ${\beta }_1=1.013,\frac{u_1}{u_4}=2.09>2$ : (a) Joint probability density function of the ship roll model; (b) Cross-section of figure (a); (c) Steady probability density function

图6. 时滞系数${\beta }_1=1.013,\frac{u_1}{u_4}=2.09>2$ ：(a) 船舶横摇模型的联合概率密度函数图；(b) 图(a)的截面图；(c) 平稳概率密度函数图

当噪声强度$\delta =0.434$ 时，从图1(a)与1(b)可以看出，联合概率密度函数呈现尖锐的峰形；而在图1(c)中，平稳概率密度函数随R变化先单调递增后递减，并于$R=0$ 处取得最大值。这表明当$R=0$ 时，船舶横摇模型在概率意义下是渐近稳定的，即系统随时间推移将逐渐趋于平衡状态。当噪声强度$\delta =0.452$ 时，如图2(a)~(b)所示，船舶横摇模型的联合概率密度函数呈现一个较为平缓的单峰形态。图2(c)与图1(c)相似，此处的函数图像同样随$R$ 值增大先上升后下降，并在$R=0$ 处取得峰值，不过其变化斜率更为平缓，此时系统处于临界状态。当噪声强度$\delta =0.513$ 时，与图2(a)和图2(b)相比，图3中的联合概率密度函数由单峰形态转变为火山口状。从图3(c)可见，平稳概率密度函数在$R=0$ 处不再取得最大值，表明系统在概率意义上已处于不稳定状态。这意味着船舶横摇模型逐渐由稳定过渡到失稳，即发生了随机分岔。进一步对比图1(c)至图3(c)可看出，随着噪声强度增大，平稳概率密度函数的峰值持续降低，其形态也逐渐趋于平缓。

在固定噪声强度的条件下，改变时滞系数取值，船舶横摇模型的随机分岔行为如图4~图6所示。当时滞系数由0.998增大至1.0012时，从图4(a)至图6(a)可看出，系统的联合概率密度函数由尖锐峰形逐渐过渡为较平缓的单峰，最终演变为火山口形态，表明系统逐渐失稳，处于平衡状态的概率持续减小。同时，从图4(c)至图6(c)可见，随着时滞系数变大，平稳概率密度函数的峰值持续减小，其斜率也逐渐趋于平缓。

5. 结论

本文基于已有非线性船舶横摇模型，引入时滞反馈与随机噪声激励，建立了带有时滞控制与随机激励的非线性船舶横摇动力学模型。首先，分析了线性系统的稳定性及其Hopf分岔条件；随后，应用中心流形理论与随机平均法，将随机时滞船舶横摇模型转化为一般随机方程，并推导出相应的随机微分方程。在此基础上，借助最大Lyapunov指数的符号判别系统局部稳定性，并通过平稳概率密度函数图像考察其全局稳定性。最后，通过数值实验研究了噪声强度和时滞系数对系统随机动力学行为的影响。