0 引　言

水下滑翔机是一种新型的水下航行器，其通过改变自身净浮力和质心位置来实现水下滑翔运动。与传统AUV相比，水下滑翔机具有成本低、能耗低、续航力强等优点，更适合大范围、长航程的工作任务，不足之处在于水下滑翔机的航速低、机动性较差^{[1, 2]}。水下滑翔机可通过调节尾舵舵角或改变横滚姿态实现转向，相比于传统的AUV，其在转向调节过程通常需要较大的空间和较长的时间^{[3, 4]}。尤其对于利用横滚姿态进行转向调节的水下滑翔机，所需的时间和空间尺度更大。将滑翔机的转向速度作为摄动参数，Mahmoudian等^[5]推导了滑翔机稳态转向时转向速度的近似解析表达式，Yang等^[6]推广此方法分析了装有尾舵的滑翔机转向性能。顾建农等^[7]基于柯西霍夫方程实现了水下滑翔机的动力学建模和转向运动仿真。朱时雨等^[8]分析了水柱测量中的水下滑翔机转向性能。为降低内部机构调节的能耗，在大多数工作时间内，水下滑翔机处于稳态运动。水下滑翔机的稳态螺旋运动，作为一种典型运动形式，常用于水下滑翔机航向的调节。Zhang等^[9]建立了一种稳态螺旋运动的快速求解方法，并与Seawing滑翔机的试验结果进行对比验证。赵宝强等^[10]采用李雅普诺夫方法进行了螺旋运动稳定性分析。Liu^[11]采用理论与试验方法分析了水下滑翔机的反螺旋运动。

当水下滑翔机在近海底或近壁面等限制水域运动时，水域的空间范围会限制其转向运动。水下滑翔机的转弯运动通常可视为稳态螺旋运动的一部分，因此，稳态螺旋运动在时间、空间上的运动尺度特征可以为水下滑翔机在限制水域内的转向运动提供指导。目前，在现有研究中尚未看到螺旋运动时空尺度特征的系统分析，文献[9]也仅是给出了下潜过程中回转半径与滑块转角的变化规律，并未分析其他控制参数和尺度特征。本文从水下滑翔机的运动学和动力学微分方程入手，建立了稳态螺旋运动的时间、空间尺度特征变化规律的研究思路，并以沈阳自动化研究所研制的Seawing水下滑翔机^{[9, 12]}（该水下滑翔机通过内部滑块旋转来实现转向）为例开展研究。

1 运动学建模

引入惯性坐标系 $E - \xi \eta \zeta $ 和体坐标系 $o - xyz$ 来描述水下滑翔机的空间运动。如图1所示，坐标轴 $E - \zeta $ 与重力方向一致，坐标系 $o - xyz$ 的原点位于滑翔机的浮心， $o - x$ 轴沿着滑翔机的机身轴线指向首部方向， $o - y$ 轴指向右侧机翼的方向。定义水下滑翔机在惯性坐标系下的位置和姿态分别为 ${{b}} \!= \!{\left[ {\begin{array}{*{20}{c}} {{\xi _{\rm{B}}}}&\!{{\eta _{\rm{B}}}}&\!{{\zeta _{\rm{B}}}} \end{array}} \right]^{\rm{T}}}$ 和 ${{\varOmega }} \!=\! {\left[ {\begin{array}{*{20}{c}} \varphi &\!\theta &\!\psi \end{array}} \right]^{\rm{T}}}$ ，定义水下滑翔机在体坐标系下的线速度和角速度分别为 ${{\upsilon }} = {\left[ {\begin{array}{*{20}{c}} u&v&w \end{array}} \right]^{\rm{T}}}$ 和 ${{\omega }} = {\left[ {\begin{array}{*{20}{c}} p&q&r \end{array}} \right]^{\rm{T}}}$ ，则水下滑翔机的运动学方程可表达为^[13]：

${{\dot b}} = {{{R}}_{{\rm{EB}}}}{{\upsilon }},{{\dot \varOmega }} = {{{R}}_{\Omega {\rm{B}}}}{{\omega }}\text{。}$

(1)

其中， ${{{R}}_{{\rm{EB}}}}$ 和 ${{{R}}_{{\rm{B}}\Omega }}$ 可用水下滑翔机的姿态角来表达。采用简化形式 ${\rm{c}} \cdot {\rm{ = cos(}} \cdot {\rm{)}}$ ， ${{s}} \cdot {\rm{ = sin(}} \cdot {\rm{)}}$ ， ${{t}} \cdot {\rm{ = tan(}} \cdot {\rm{)}}$ 表达三角函数，则有

$\begin{split} & {{{R}}_{{\rm{EB}}}} \!=\! \left(\!{\begin{array}{*{20}{c}} \!\!\!{{\rm{c}}\theta {\rm{c}}\psi }&{{\rm{s}}\varphi {\rm{s}}\theta {\rm{c}}\psi - {\rm{c}}\varphi {\rm{s}}\psi }&{{\rm{c}}\varphi {\rm{s}}\theta {\rm{c}}\psi + {\rm{s}}\varphi {\rm{s}}\psi }\\ \!\!\!{{\rm{c}}\theta {\rm{s}}\psi }&{{\rm{c}}\varphi {\rm{c}}\psi + {\rm{s}}\varphi {\rm{s}}\theta {\rm{s}}\psi }&{ - {\rm{s}}\varphi {\rm{c}}\psi + {\rm{c}}\varphi {\rm{s}}\theta {\rm{s}}\psi }\\ \!\!\!{ - {\rm{s}}\theta }&{{\rm{s}}\varphi {\rm{c}}\theta }&{{\rm{c}}\varphi {\rm{c}}\theta } \end{array}} \!\!\!\right)\text{，}\\ & {{{R}}_{{\rm{B}}\Omega }} = \left( {\begin{array}{*{20}{c}} 1&0&{ - {\rm{s}}\theta }\\ 0&{{\rm{c}}\varphi }&{{\rm{s}}\varphi {\rm{c}}\theta }\\ 0&{ - {\rm{s}}\varphi }&{{\rm{c}}\varphi {\rm{c}}\theta } \end{array}} \right)\text{。} \end{split}$

(2)

2 动力学建模

水下滑翔机可以视为不同种结构组成的多体系统，包括机身及与机身固结在一起的所有结构m_rb，可在机身内部平移和旋转的滑块结构m_p，位置不变质量可变的浮力调节机构m_b。m_rb，m_p，m_b的质心在体坐标系中的位置分别为向量 ${{{r}}_{rb}}$ ， ${{{r}}_{p}}$ ， ${{{r}}_{b}}$ 。 ${r_{rx}}$ 为m_p的质心沿着 $o - x$ 轴的位置， $\gamma $ 为m_p绕 $o - x$ 轴的转角。其中， $\gamma = 0$ 对应m_p的质心处在体坐标的 $o - xz$ 平面内，此时滑翔机将在垂直平面内做直线运动； $\gamma \ne 0$ 将导致滑翔机所受的升力偏离垂直平面，升力在水平面中的分量将作为向心力使其产生螺旋转向运动，如图2所示。

水下滑翔机在水下运动时，会带动其周围的附连水共同运动，水下滑翔机及其附连水系统的总体动能可以表达为：

$T{\rm{ = }}\frac{1}{2}{\left(\!\! {\begin{array}{*{20}{c}} {{\upsilon }} \\ {{\omega }} \end{array}}\!\! \right)^{\rm{T}}}{{M}}\left(\!\! {\begin{array}{*{20}{c}} {{\upsilon }} \\ {{\omega }} \end{array}} \!\!\right){\rm{ = }}\frac{1}{2}{\left(\!\! {\begin{array}{*{20}{c}} {{\upsilon }} \\ {{\omega }} \end{array}} \!\!\right)^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {{{{M}}_t}}&{{{{C}}_t}} \\ {{{C}}_t^T}&{{{{I}}_t}} \end{array}} \right]\left(\!\! {\begin{array}{*{20}{c}} {{\upsilon }} \\ {{\omega }} \end{array}} \!\!\right)\text{，}$

(3)

$\begin{gathered} {{{M}}_t} = {{{M}}_{{A}}} + ({m_{{\rm{rb}}}} + {m_{\rm{p}}} + {m_{\rm{b}}}){{I}} \text{，} \\ {{{C}}_t} = {{{C}}_{{A}}} - {m_{{\rm{rb}}}}{{{{\hat r}}}_{{\rm{rb}}}} - {m_{\rm{p}}}{{{{\hat r}}}_{\rm{p}}} - {m_{\rm{b}}}{{{{\hat r}}}_{\rm{b}}} \text{，}\\ {{{I}}_t} = {{{J}}_{{A}}} + {{{J}}_{{\rm{rb}}}} - {m_{\rm{p}}}{{{{\hat r}}}_{\rm{p}}}{{{{\hat r}}}_{\rm{p}}} - {m_{\rm{b}}}{{{{\hat r}}}_{\rm{b}}}{{{{\hat r}}}_{\rm{b}}} \text{。}\\ \end{gathered} $

(4)

其中： ${{M}}$ 为系统的广义质量矩阵； ${{{M}}_{{A}}}$ ， ${{{J}}_{{A}}}$ 和 ${{{C}}_{{A}}}$ 分别为体坐标系下附连水的质量矩阵，转动惯量矩阵和耦合项矩阵； ${{{J}}_{{\rm{rb}}}}$ 为体坐标系固定质量的转动惯量矩阵； ${{I}}$ 为单位对角矩阵。

则系统在体坐标系下的线动量 ${{P}}$ 和角动量 ${{\varPi }}$ 可利用其动能表达为

${{P}} = \frac{{\partial T}}{{\partial {{\upsilon }}}},{{\varPi }} = \frac{{\partial T}}{{\partial {{\omega }}}}\text{，}$

(5)

进一步表达有

$\left( {\begin{array}{*{20}{c}} {{P}} \\ {{\varPi }} \end{array}} \right){\rm{ = }}{{M}}\left( {\begin{array}{*{20}{c}} {{\upsilon }} \\ {{\omega }} \end{array}} \right)\text{。}$

(6)

在稳态运动时，系统的广义质量矩阵将保持不变，即 ${{\dot M}} = 0$ 。对式(6)两端求导，可得系统的动量变化率为

$\left( {\begin{array}{*{20}{c}} {{{\dot P}}} \\ {{{\dot \Pi }}} \end{array}} \right){\rm{ = }}{{M}}\left( {\begin{array}{*{20}{c}} {{{\dot \upsilon }}} \\ {{{\dot \omega }}} \end{array}} \right)\text{。}$

(7)

基于动量定理和坐标变换，可得系统的动量变化率 ${{\dot P}}$ ， ${{\dot \Pi }}$ 与其所受外力的关系

${{\dot P}} = {{P}} \times {{\omega }} + {{R}}_{{EB}}^{\rm{T}}({{G}} + {{B}}) + {{{F}}_{\rm{B}}}\text{，}$

(8)

${{\dot \Pi }} = {{\Pi }} \times {{\omega }} + {{P}} \times {{\upsilon }} + {{{r}}_{{\rm{CG}}}} \times {{R}}_{{EB}}^{\rm{T}}{{G}} + {{{M}}_{\rm{B}}}\text{。}$

(9)

其中： ${{G}}$ 和B分别为水下滑翔机的重力和浮力， ${{{r}}_{{\rm{CG}}}}$ 为体坐标系下的重心位置向量， ${{{F}}_{{B}}}$ 和 ${{{M}}_{{B}}}$ 分别为体坐标系下水下滑翔机的粘性水动力和水动力矩。对于 ${{{F}}_{{B}}}$ 和 ${{{M}}_{{B}}}$ ，可结合滑翔机运动速度和水动力系数来计算^{[9, 13]}。

将式（6）、式（7）代入式（8）、式（9），即得到水下滑翔机稳态运动的动力学方程：

$\left(\!\! {\begin{array}{*{20}{c}} {{{\dot \upsilon }}} \\ {{{\dot \omega }}} \end{array}} \!\!\right){\rm{ = }}{{{M}}^{ - 1}}\left(\!\! {\begin{array}{*{20}{c}} \begin{gathered} ({{{M}}_t}{{\upsilon }} + {{{C}}_t}{{\omega }}) \times {{\omega }}+ \\ {{R}}_{{\rm{EB}}}^{\rm{T}}({{G}} + {{B}}) + {{{F}}_{\rm{B}}} \\ \end{gathered} \\ \begin{gathered} ({{C}}_t^T{{\upsilon }} + {{{I}}_t}{{\omega }}) \times {{\omega }} + ({{{M}}_t}{{\upsilon }} + {{{C}}_t}{{\omega }}) \times {{\upsilon }}+ \\ \begin{array}{*{20}{c}} {}&{}&{} \end{array} {{{r}}_{{\rm{CG}}}} \times {{R}}_{{\rm{EB}}}^{\rm{T}}{{G}} + {{{M}}_{\rm{B}}} \\ \end{gathered} \end{array}} \!\!\right)\text{。}\!\!\!\!\!\!\!\!$

(10)

3 稳态螺旋运动的时空描述

基于水下滑翔机的运动学方程和动力学方程，当给定其控制输入量和 ${{b}}$ ， ${{\varOmega }}$ ， ${{\upsilon }}$ ， ${{\omega }}$ 的初值后，即可实现水下滑翔机的三维空间运动仿真。对于依靠质心移动实现转向的水下滑翔机，其运动姿态的控制输入量包括移动滑块沿 $o - x$ 轴的位置 ${r_{{{rx}}}}$ 、移动滑块绕 $o - x$ 轴的转角 $\gamma $ 和净浮力质量m_b。本文以Seawing水下滑翔机为例进行分析，水下滑翔机自身的质量参数和水动力参数见文献[9]。将运动状态量 ${{b}}$ ， ${{\varOmega }}$ ， ${{\upsilon }}$ ， ${{\omega }}$ 的初值取为零向量，并将控制参数取值为 $\gamma {\rm{ = 2}}{{\rm{0}}{\text{°}}}$ ， ${r_{{{rx}}}}{\rm{ = 0}}{\rm{.42\;m}}$ ，m_b=0.3 kg，得到水下滑翔机的运动轨迹、运动速度变化规律分别如图3和图4所示。从图3可以看出，可采用回转半径 $R$ ，螺距 $H$ 和回转周期 $T$ 来表达水下滑翔机的空间和时间尺度特征。其中，回转半径 $R$ 的方向与滑翔机绕 $o - \zeta $ 轴的角速度方向一致，螺距 $H$ 的方向与滑翔机沿 $o - \zeta $ 轴的线速度方向一致，回转周期 $T$ 为滑翔机一个周期回转运动的时间。若水下滑翔机需要转向某一角度，可以利用转向角数值和螺旋运动的时空尺度特征，换算得到对应转向运动所需的时间、空间量。

从图4可以看出，随着时间的变化，水下滑翔机的运动速度 ${{\upsilon }}$ ， ${{\omega }}$ 将逐渐稳定，相应的水下滑翔机运动也将逐渐达到稳态。为进一步定量地验证本文稳态运动仿真方法的正确性，将本文回转半径的分析结果与文献[9]中的结果进行对比，如图5所示，其中 ${r_{{{rx}}}}{\rm{ = 0}}{\rm{.421\;6\;m}}$ , m_b=0.3 kg。可以看出，2种方法结果的变化趋势一致，但量值上存在较小的差别，原因为文献[9]中为采用递归方法实现稳态螺旋运动的快速分析，在推导过程中忽略了质量矩阵中影响较小的耦合项，而本文方法从滑翔机的运动微分方程入手，更为详细地计入了质量矩阵。总体看来，本文的运动仿真方法可用于水下滑翔机稳态螺旋运动时空尺度的分析。

4 时空尺度的分析

根据水下滑翔机沿垂向的运动方向进行分类，水下滑翔机的稳态螺旋运动可分为上浮螺旋运动和下潜螺旋运动。当重力大于浮力，即m_b>0时，滑翔机将进行下潜螺旋运动，此时需调节 ${r_{{{rx}}}}$ 使水下滑翔机的重心在浮心的前方，从而保证滑翔机沿着机头的方向滑行。而当重力小于浮力，即m_b<0时，滑翔机将进行上浮螺旋运动，此时需调节 ${r_{{{rx}}}}$ 使滑翔机的重心在浮心的后方。以Seawing水下滑翔机为例进行研究，通过改变滑翔机下潜和上浮运动中的控制参数，系统地设计一系列稳态运动计算工况，进而分析水下滑翔机的时间、空间尺度特征随控制输入参数间的变化规律。

在下潜和上浮运动中保持 ${r_{{{rx}}}}$ 和m_b的取值不变，通过改变 $\gamma $ 进行稳态运动的计算，可得螺旋运动的时空尺度特征随转角 $\gamma $ 的变化规律，如图6所示。其中， ${r_{{{rx}}}}{\rm{ = 0}}{\rm{.42\;m}}$ , m_b=0.3 kg，对应下潜螺旋运动； ${r_{{{rx}}}}{\rm{ = 0}}{\rm{.42\;m}}$ , m_b=0.3 kg，对应上浮螺旋运动。可以看出，回转半径 $R$ 关于坐标原点中心对称，其方向与转角 $\gamma $ 的方向一致；螺距 $H$ 和回转周期 $T$ 关于 $\gamma {\rm{ = 0}}$ 位置左右对称。随着转角 $\gamma $ 绝对值的增大，上浮和下潜螺旋运动中回转半径 $R$ 、螺距 $H$ 和回转周期 $T$ 的绝对值均减小。这表明在有限尺度的空间中转向时，可考虑增大 $\gamma $ 来降低转向过程中所需的时间、空间尺度。

在下潜和上浮运动中保持 $\gamma $ 和m_b的取值不变，通过改变 ${r_{{{rx}}}}$ 进行稳态运动的计算，可得到时空尺度特征随 ${r_{{{rx}}}}$ 的变化规律，如图7所示。其中，m_b=0.3 kg, $\gamma {\rm{ = 9}}{{\rm{0}}{\text{°}}}$ 和m_b=–0.3 kg, $\gamma {\rm{ = 9}}{{\rm{0}}{\text{°}}}$ ，分别对应下潜螺旋运动和上浮螺旋运动。可以看出，在下潜螺旋运动中，回转半径 $R$ 、螺距 $H$ 和回转周期 $T$ 均存在最大值；在上浮螺旋运动中，回转半径 $R$ 存在最大值，而且上浮过程中的回转周期 $T$ 要小于下潜过程。 ${r_{{{rx}}}}$ 在0.4~0.405 m附近时，水下滑翔机的重心和浮心的纵向位置比较接近。尽管在最值位置两侧的单调范围中调节 ${r_{{{rx}}}}$ ，均可降低所需的时空尺度，但总体来看，当调节 ${r_{{{rx}}}}$ 来减小水下滑翔机的重心和浮心的纵向距离（使 ${r_{{{rx}}}}$ 接近0.4~0.405 m范围），更有助于快速地降低螺旋转向过程中的时间、空间尺度。

在下潜和上浮运动中保持 $\gamma $ 和 ${r_{{{rx}}}}$ 的取值不变，通过改变m_b进行稳态运动的计算，可得到螺旋运动的时空尺度随m_b的变化规律，如图7所示。其中， ${r_{{{rx}}}}{\rm{ = 0}}{\rm{.42\;m}}$ ， $\gamma {\rm{ = 9}}{{\rm{0}} ^\circ }$ 和 ${r_{{\rm{rx}}}}{\rm{ = 0}}{\rm{.38\;m}}$ ， $\gamma {\rm{ = 9}}{{\rm{0}}^\circ}$ ，分别对应下潜螺旋运动和上浮螺旋运动。可以看出，回转半径 $R$ 、螺距 $H$ 和回转周期 $T$ 在下潜运动中的变化规律与其在上浮运动中的变化规律完全不同；在下潜螺旋运动中，随着m_b绝对值的增大，回转半径 $R$ 、螺距 $H$ 和回转周期 $T$ 单调递增；在上浮螺旋运动中，随着m_b绝对值的增大，回转半径 $R$ 先增大后减小，回转周期 $T$ 单调递减，螺距 $H$ 呈波动变化。因此，当通过调节m_b来降低螺旋转向运动的时间、空间尺度时，可考虑在下潜螺旋运动中减小m_b绝对值，而由于上浮螺旋运动中的变化规律复杂，应根据转向时刻m_b的具体量值来判断。

5 结　语

本文通过对水下滑翔机稳态螺旋运动的时间、空间尺度进行分析，得到如下结论：1）建立水下滑翔机的稳态螺旋运动仿真方法，并验证了该方法的正确。2）系统地分析了稳态螺旋运动的时空尺度特征随控制输入参数间的变化规律，发现了上浮和下潜运动中变化规律的不同之处。3）当通过调节移动滑块的转角或纵向位置来降低螺旋运动中的时间、空间尺度时，可考虑增大滑块转角、减小水下滑翔机的重心和浮心的纵向距离来实现。4）当通过调节净浮力来降低螺旋运动中的时间、空间尺度时，在下潜过程中可考虑降低净浮力，在上浮过程中应根据转向时刻m_b的具体量值来调节。