由于海洋资源的多样性和丰富程度，海洋探测技术已成为近些年的研究热点。对此，各国学者研制了多种类型的水下探测机器人。1989年，美国海洋学家Stommel首先提出水下滑翔机的概念^[1]。通常，滑翔机设计为无外加推进器的有翼水下航行器；通过调节外油囊体积产生净浮力，实现深度方向的上升与下降；通过调整内部可动质量块的位置，实现滑翔机重心的改变，从而实现姿态调节；而施加在机翼上的水动力，使滑翔机向前移动。水下滑翔机具备体积小、成本低、续航时间长等优点，受到了广泛关注。目前，典型的水下滑翔机包括Slocum^[2]、Spray^[3]、Seaglider^[4]、海翼号^[5]及海燕号^[6]等。

水下滑翔机研制过程中，采用动力学仿真方法研究其性能已成为重要环节之一。目前，各国学者已采用Newton-Euler法^[7]、Lagrange方程^[8]、Gibbs-Appell方程^[9]、微分几何^[10]等方法成功建立了多种型号的水下滑翔机动力学模型，实现了空间运动的仿真分析。基于水下滑翔机的动力学模型，Yu等^[11]建立了滑翔机的能耗模型，实现了单剖面滑翔机运动的能耗分析；Fan等^[12]研究了非均匀水流中水下滑翔机的控制方法；Yoon等^[13]以最大化前进速度为优化目标，研究了滑翔机运动轨迹的优化方法；于鹏垚等^[14]研究了水下滑翔机稳态运动的统一描述方法；顾建农等^[15]研究了水下滑翔机在恒定水流中的运动特性；Niu等^[16]确定了混合驱动水下滑翔机运动参数的吸引域，并开展了稳定性研究；Wang等^[17]研究了海洋生物附着对滑翔机运动特性的影响；Wu等^[18-20]考虑运动精度指标，开展了水下滑翔机性能分析和优化的研究。

在一些探测任务中，水下滑翔机需要下潜到水下特定位置，进行虚拟系泊，从而开展定点探测。显然，该工作模式的关键环节是水下滑翔机可以准确下潜到目标探测区域。考虑海洋环境的多变性以及滑翔机携带的电池包供电能力有限，人们希望滑翔机可以消耗尽量少的能量、尽快到达目标探测区域。为了实现上述目标，需要确定水下滑翔机的最优控制参数配置，但鲜有文献报道该项研究工作。

本文将以典型的水下滑翔机为研究对象，面向水下定点探测任务，开展滑翔机控制参数优化方法的研究，建立了滑翔机整机动力学模型和滑翔机下潜运动的能耗模型，分析了滑翔机运动轨迹随控制参数的变化规律，建立了计算性能评估参数的代理模型，结合NSGA-Ⅱ和代理模型求解最优控制参数值。

1 水下滑翔机动力学建模

本文以典型水下滑翔机为研究对象，整机布局如图 1所示。

本文不考虑水下滑翔机的转弯运动，因此仅研究可动质量块可以前后平移的情况。此外，为了更有针对性地研究控制参数对滑翔机运动的影响，不考虑水流的作用。有鉴于此，滑翔机的控制参数包括净浮力调节量c₁和可动质量块平移量c₂。

1.1 坐标系和运动学方程

为了方便后续推导，将滑翔机分为机体和可动质量块2部分进行研究，机体包含除可动质量块之外的全部结构。为描述机体和可动质量块的运动，定义了3个坐标系：惯性坐标系(I坐标系)、机体的随体坐标系(B坐标系)、可动质量块的随体坐标系(P坐标系)。I坐标系固联在惯性空间之上，由X、Y和Z轴组成。B坐标系固联在机体之上，由X_B、Y_B和Z_B轴组成。P坐标系固联在可动质量块之上，由X_P、Y_P和Z_P轴组成。3个坐标系的定义与文[20]相同，如图 2所示。B坐标系的原点O_B为c₁等于0时的整机浮心，P坐标系的原点O_P为可动质量块质心。此外，外油囊浮心和机体质心分别为O_e和O_h。

为了方便后续推导，对将用到的变量进行定义，如表 1所示。

由于本文不考虑滑翔机的转弯运动以及水流的影响，机体的位置Y轴坐标y、横滚角φ、偏航角ψ、横向速度v、横滚角速度p、偏航角速度r始终等于0。根据几何关系，从B坐标系到I坐标系的变换矩阵R_BI表示为

$ \boldsymbol{R}_{\mathrm{BI}}=\left(\begin{array}{ccc} \cos \theta & 0 & \sin \theta \\ 0 & 1 & 0 \\ -\sin \theta & 0 & \cos \theta \end{array}\right) \text {. } $

(1)

在此基础上，从I坐标系到B坐标系的变换矩阵R_IB满足R_IB=R_BI^T。根据运动学关系，则有

$ \left\{\begin{array}{l} \dot{\boldsymbol{b}}=\boldsymbol{R}_{\mathrm{BI}} \boldsymbol{V}_{\mathrm{b}}, \\ \dot{\boldsymbol{\pi}}=\left(\begin{array}{ccc} 1 & 0 & \tan \theta \\ 0 & 1 & 0 \\ 0 & 0 & \sec \theta \end{array}\right) \mathit{\pmb{\Omega}}_{\mathrm{b}} . \end{array}\right. $

(2)

参考文[17]，定义运算符^，含义如下：

$ \left\{\begin{array}{l} \boldsymbol{q}=\left(\begin{array}{lrr} q_1 & q_2 & q_3 \end{array}\right)^{\mathrm{T}}, \\ \hat{\boldsymbol{q}}=\left(\begin{array}{rrr} 0 & -q_3 & q_2 \\ q_3 & 0 & -q_1 \\ -q_2 & q_1 & 0 \end{array}\right) . \end{array}\right. $

(3)

若不考虑控制参数的具体调节过程，根据运动学关系可得

$ \left\{\begin{array}{l} \boldsymbol{V}_{\mathrm{p}}=\boldsymbol{V}_{\mathrm{b}}+\mathit{\pmb{\Omega}}_{\mathrm{b}} \times \boldsymbol{r}_{\mathrm{p}}, \\ \mathit{\pmb{\Omega}}_{\mathrm{p}}=\mathit{\pmb{\Omega}}_{\mathrm{b}} . \end{array}\right. $

(4)

参考文[18]，r_p的表达式如下：

$ \boldsymbol{r}_{\mathrm{p}}=\left(0.09 \mathrm{~m}+c_2 \quad 0 \quad 0.016 \mathrm{~m}\right)^{\mathrm{T}} . $

(5)

1.2 受力分析

滑翔机在运动过程中受到的外力包括：重力、浮力、水压力、惯性水动力、黏性水动力。这里，将外力均平移至B坐标系的原点O_B。由于一些外力的作用点不是O_B，平移之后会产生附加转矩。若无特殊说明，所有外力均在B坐标系下表示。

根据重力的性质，滑翔机的重力表达式如下：

$ \left\{\begin{array}{lll} \boldsymbol{G}_{\mathrm{b}}=\boldsymbol{R}_{\mathrm{eb}}\left(\begin{array}{lll} 0 & 0 & \left.m_{\mathrm{b}} g\right)^{\mathrm{T}}, \end{array}\right. \\ \boldsymbol{G}_{\mathrm{p}}=\boldsymbol{R}_{\mathrm{eb}}\left(\begin{array}{lll} 0 & 0 & \left.m_{\mathrm{p}} g\right)^{\mathrm{T}} \end{array} .\right. \end{array}\right. $

(6)

其中：G_b和G_p分别表示机体和可动质量块的重力，g表示重力加速度。

将G_b和G_p平移至O_B后，产生的附加力矩表达式如下：

$ \left\{\begin{array}{l} \boldsymbol{M}_{\mathrm{b}}=\boldsymbol{r}_{\mathrm{h}} \times \boldsymbol{G}_{\mathrm{b}}, \\ \boldsymbol{M}_{\mathrm{p}}=\boldsymbol{r}_{\mathrm{p}} \times \boldsymbol{G}_{\mathrm{p}} \end{array}\right. $

(7)

当c₁等于0且滑翔机悬浮在水面上时，整机净浮力表示为B₀，此时整机浮力与重力大小相等、方向相反，因此可得

$ \boldsymbol{B}_0=-\left(\boldsymbol{G}_{\mathrm{b}}+\boldsymbol{G}_{\mathrm{p}}\right) . $

(8)

当c₁不等于0时，作用于外油囊浮心O_e处的净浮力表示为

$ \Delta \boldsymbol{B}=\boldsymbol{R}_{\text {IB }}\left(\begin{array}{lll} 0 & 0 & c_1 \end{array}\right)^{\mathrm{T}} . $

(9)

将ΔB平移至O_B后，产生的附加力矩表达式如下：

$ \boldsymbol{M}_{\mathrm{e}}=\boldsymbol{r}_{\mathrm{e}} \times \Delta \boldsymbol{B} \boldsymbol{.} $

(10)

随着水深的增加，水密度和水压力通常会逐渐增大，且水压力会使滑翔机的耐压壳体发生变形，从而产生浮力偏差。考虑水深z的影响，水密度ρ和耐压壳体体积V_H的表达式为

$ \left\{\begin{array}{l} \rho=\rho_0+\rho_1 z+\rho_2 z^2+\rho_3 z^3+\rho_4 z^4, \\ V_{\mathrm{H}}=V_{\mathrm{H} 0}-K_{\mathrm{H}} z . \end{array}\right. $

(11)

其中：ρ_i表示拟合系数，V_H0为在一个标准大气压(101 325 Pa)下滑翔机耐压壳体的体积，K_H为耐压壳体的压缩率。ρ_i、V_H0、K_H取值与文[21]一致。V_H随水深变化的函数如图 3所示。

从而，考虑水密度变化和水压力的影响，在O_B产生的浮力偏差可表示为^[21]

$ \boldsymbol{B}_{\mathrm{d}}=\boldsymbol{R}_{\mathrm{IB}}\left(\begin{array}{lll} 0 & 0 & \left.\left(\rho_0 V_{\mathrm{H} 0}-\rho V_{\mathrm{H}}\right) g\right)^{\mathrm{T}} . \end{array}\right. $

(12)

图 3表明，V_H的变化量较小，因此耐压壳体变形对滑翔机水动力性能的影响可以忽略不计。除耐压壳体外，与流体介质接触的结构均为非空心结构且体积较小，故水压力对其他结构的影响可以忽略不计。

从而，滑翔机所受的总浮力及附加力矩表示为

$ \left\{\begin{array}{l} \boldsymbol{F}_{\text {buо }}=\boldsymbol{B}_0+\Delta \boldsymbol{B}+\boldsymbol{B}_{\mathrm{d}}, \\ \boldsymbol{M}_{\text {net }}=\boldsymbol{r}_{\mathrm{e}} \times \Delta \boldsymbol{B} . \end{array}\right. $

(13)

当滑翔机在水中作变速运动时，将带动周围流体质点产生相应的速度变化，进而产生相互作用力，称为惯性水动力。通常将惯性水动力与滑翔机的加速度和角加速度视为线性关系，根据势流理论，并考虑滑翔机外形的对称性，可得^[18]

$ \left(\begin{array}{l} \boldsymbol{F}_{\mathrm{m}} \\ \boldsymbol{M}_{\mathrm{m}} \end{array}\right)=-\rho\left(\begin{array}{cccccc} \lambda_{11} & 0 & 0 & 0 & 0 & 0 \\ 0 & \lambda_{22} & 0 & 0 & 0 & \lambda_{26} \\ 0 & 0 & \lambda_{33} & 0 & \lambda_{35} & 0 \\ 0 & 0 & 0 & \lambda_{44} & 0 & 0 \\ 0 & 0 & \lambda_{53} & 0 & \lambda_{55} & 0 \\ 0 & \lambda_{62} & 0 & 0 & 0 & \lambda_{66} \end{array}\right)\left(\begin{array}{c} \dot{\boldsymbol{V}}_{\mathrm{b}} \\ \dot{\mathit{\pmb{\Omega}}}_{\mathrm{b}} \end{array}\right) . $

(14)

其中λ_ij表示惯性水动力系数，取值与滑翔机的水动力外形密切相关。当滑翔机的水动力外形确定后，惯性水动力系数可视为定值，并可基于文[22]提到的CFD仿真方法求解。

滑翔机运动过程中，还受到与速度和角速度相关的阻尼力，称为黏性水动力，文[20]给出了滑翔机做三维运动时的黏性水动力方程。本文不考虑水流影响及滑翔机的转弯运动，故滑翔机运动过程中的侧滑角β、p和r始终为0。因此，将文[20]中的黏性水动力方程进一步简化为

$ \left\{\begin{array}{l} \boldsymbol{F}_{\mathrm{s}}=\rho\left\|\boldsymbol{V}_{\mathrm{b}}\right\|^2 \boldsymbol{R}_\alpha\left(\begin{array}{c} -K_{\mathrm{D} 0}-K_{\mathrm{D}} \alpha^2 \\ 0 \\ -K_{\mathrm{L} 0}-K_\alpha \alpha \end{array}\right), \\ \boldsymbol{M}_{\mathrm{s}}=\rho\left\|\boldsymbol{V}_{\mathrm{b}}\right\|^2 \boldsymbol{R}_\alpha\left(\begin{array}{c} 0 \\ K_{\mathrm{M} 0}+K_{\mathrm{M}} \alpha+K_{\mathrm{q}} q \\ 0 \end{array}\right) . \end{array}\right. $

(15)

其中：α表示攻角；R_α表示与攻角相关的变换矩阵；带有下标的K表示黏性水动力系数，取值与滑翔机的水动力外形密切相关。当滑翔机的水动力外形确定后，黏性水动力系数可视为定值，并可基于文[23]提到的CFD仿真方法求解。α和R_α的表达式如下：

$ \left\{\begin{array}{l} \alpha=\arctan \left(\frac{w}{u}\right), \\ \boldsymbol{R}_\alpha=\left(\begin{array}{ccc} \cos \alpha & 0 & -\sin \alpha \\ 0 & 1 & 0 \\ \sin \alpha & 0 & \cos \alpha \end{array}\right) . \end{array}\right. $

(16)

其中u和w分别表示水下滑翔机在X_B和Z_B轴方向的速度分量。

综上，滑翔机在运动过程中受到的合外力及合外力矩为

$ \left\{\begin{array}{l} \boldsymbol{F}_{\mathrm{all}}=\boldsymbol{G}_{\mathrm{b}}+\boldsymbol{G}_{\mathrm{p}}+\boldsymbol{F}_{\mathrm{buo}}+\boldsymbol{F}_{\mathrm{m}}+\boldsymbol{F}_{\mathrm{s}}, \\ \boldsymbol{M}_{\mathrm{all}}=\boldsymbol{M}_{\mathrm{b}}+\boldsymbol{M}_{\mathrm{p}}+\boldsymbol{M}_{\mathrm{net}}+\boldsymbol{M}_{\mathrm{m}}+\boldsymbol{M}_{\mathrm{s}} . \end{array}\right. $

(17)

1.3 动力学方程

参考文[19]，直接给出水下滑翔机的整机动力学方程：

$ \left\{\begin{array}{l} \left(m_{\mathrm{b}}+m_{\mathrm{p}}\right) \dot{\boldsymbol{V}}_{\mathrm{b}}-\left(m_{\mathrm{b}} \hat{\boldsymbol{r}}_{\mathrm{h}}+m_{\mathrm{p}} \hat{\boldsymbol{r}}_{\mathrm{p}}\right) \dot{\mathit{\pmb{\Omega}}}_{\mathrm{b}}= \\ \left(m_{\mathrm{b}} \boldsymbol{V}_{\mathrm{b}}+m_{\mathrm{p}} \boldsymbol{V}_{\mathrm{p}}\right) \times \mathit{\pmb{\Omega}}_{\mathrm{b}}+m_{\mathrm{b}} \mathit{\pmb{\Omega}}_{\mathrm{b}} \times \boldsymbol{r}_{\mathrm{h}} \times \mathit{\pmb{\Omega}}_{\mathrm{b}}+\boldsymbol{F}_{\mathrm{all}}, \\ \left(\boldsymbol{J}_{\mathrm{b}}+\boldsymbol{J}_{\mathrm{p}}-m_{\mathrm{p}} \hat{\boldsymbol{r}}_{\mathrm{p}} \hat{\boldsymbol{r}}_{\mathrm{p}}\right) \hat{\mathit{\pmb{\Omega}}}_{\mathrm{b}}+\left(m_{\mathrm{b}} \hat{\boldsymbol{r}}_{\mathrm{h}}+m_{\mathrm{p}} \hat{\boldsymbol{r}}_{\mathrm{p}}\right) \dot{\boldsymbol{V}}_{\mathrm{b}}= \\ \left(\boldsymbol{J}_{\mathrm{p}} \mathit{\pmb{\Omega}}_{\mathrm{p}}\right) \times \mathit{\pmb{\Omega}}_{\mathrm{p}}+\left(\boldsymbol{J}_{\mathrm{b}} \mathit{\pmb{\Omega}}_{\mathrm{b}}+m_{\mathrm{b}} \boldsymbol{r}_{\mathrm{h}} \times \boldsymbol{V}_{\mathrm{b}}\right) \times \mathit{\pmb{\Omega}}_{\mathrm{b}}+ \\ m_{\mathrm{b}} \mathit{\pmb{\Omega}}_{\mathrm{b}} \times \boldsymbol{r}_{\mathrm{h}} \times \boldsymbol{V}_{\mathrm{b}}+m_{\mathrm{p}} \boldsymbol{r}_{\mathrm{p}} \times\left(\boldsymbol{V}_{\mathrm{p}} \times \mathit{\pmb{\Omega}}_{\mathrm{b}}\right)+\boldsymbol{M}_{\mathrm{all}} . \end{array}\right. $

(18)

联立式(2)和(18)，即可完成水下滑翔机的动力学建模。模型中，滑翔机的几何参数取值、质量参数取值、水动力系数取值与文[20]相同。

本文采用四阶显式Runge-Kutta法完成动力学仿真计算，步长设置为变步长，并利用五阶方法控制计算误差，方法的整体截断误差为步长的五次方。

2 能耗评估模型

本文主要研究滑翔机从水面下潜到水下特定位置的运动过程。首先建立滑翔机下潜运动的能耗模型。设定水下滑翔机控制参数的取值范围：c₁为2.5~5.0 N，c₂为0.01~0.04 m。

当滑翔机收到下潜指令后，将调节外油囊体积和可动质量块平移量，开始下潜运动。当下潜到目标深度后，理论上，需要将c₁和c₂均调节为0，开始虚拟系泊，开展定点探测。上述过程中，滑翔机的能量主要被净浮力调节系统、姿态调节系统、控制系统和测量系统消耗。

上述工作流程中，净浮力可视为调整2次。准备下潜时，滑翔机将净浮力从0调整到c₁，该过程是通过打开电磁阀，使液压油在大气压力作用下，从外油囊流入滑翔机内部实现。达到目标下潜深度后，滑翔机将净浮力从c₁调整到0，该过程是通过打开液压泵，将机体内部的液压油泵入外油囊实现。考虑外油囊体积变化量较小，且内部液压油视为不可压缩，因此水密度变化对外油囊的净浮力影响较小。有鉴于此，根据Archimedes定律，外油囊的体积变化量为

$ \Delta V=\frac{c_1}{\rho g} . $

(19)

设净浮力调节系统中，电磁阀的功率和流量分别表示为p_v和q_v，液压泵的功率和流量分别表示为p_b和q_b。根据文[24]，p_b和q_b是关于目标下潜深度d_T的函数，具体表达式如下：

$ \left\{\begin{array}{l} p_{\mathrm{b}}=\left(28.212+0.017\ 841 d_{\mathrm{T}} / \mathrm{m}\right) \mathrm{W}, \\ q_{\mathrm{b}}=\left(1.7 \times 10^{-6}-1.98 d_{\mathrm{T}} / \mathrm{m} \times 10^{-10}\right) \mathrm{m}^3 / \mathrm{s} . \end{array}\right. $

(20)

进一步，确定净浮力调节系统的能耗为

$ E_1=\Delta V\left(\frac{p_{\mathrm{b}}}{q_{\mathrm{b}}}+\frac{p_{\mathrm{v}}}{q_{\mathrm{v}}}\right) . $

(21)

当滑翔机准备下潜时，可动质量块平移量从0调整到c₂，设可动质量块的移动速度为v_m，驱动功率为p_m，根据文[24]，p_m与滑翔机的俯仰角θ密切相关，且满足：

$ p_{\mathrm{m}}(\theta)=(1.43+0.064|\theta|) \mathrm{W} . $

(22)

进一步，确定姿态调节系统的能耗为

$ E_2=\frac{c_2}{v_{\mathrm{m}}} p_{\mathrm{m}}(0)+\frac{c_2}{v_{\mathrm{m}}} p_{\mathrm{m}}\left(\theta_{\mathrm{T}}\right) . $

(23)

其中θ_T表示滑翔机下潜到目标深度时的俯仰角。

控制系统在下潜阶段持续工作，平均功率为p_c。设滑翔机下潜运动的总时间为t_d，则控制系统的能耗为

$ E_3=t_{\mathrm{d}} p_{\mathrm{c}}. $

(24)

水下滑翔机通常配有2种传感器：一种是持续工作传感器，平均功率表示为p_cs；另一种是深度间隔工作传感器，每隔一定下潜深度开启1次，采样深度间隔为Δd，单次开启时间为Δt，平均功率为p_is。因此，可以确定测量系统的能耗为

$ E_4=t_{\mathrm{d}} p_{\mathrm{cs}}+\frac{d_{\mathrm{T}} \Delta t}{\Delta t\left|v_{\mathrm{d}}\right|+\Delta d} p_{\text {is }} . $

(25)

其中v_d表示水下滑翔机下潜运动的平均速度。

$ v_{\mathrm{d}}=\frac{d_{\mathrm{T}}}{t_{\mathrm{d}}}. $

(26)

基于上述分析，确定滑翔机执行水下定点探测任务时，运动过程的总能耗为

$ E=\frac{\left(E_1+E_2+E_3+E_4\right)}{1000}. $

(27)

能耗模型中，关键参数取值与文[19]相同。

将c₁、c₂和d_T取值代入能耗模型，并利用动力学仿真求解θ_T和t_d，即可计算水下滑翔机的能耗。

3 下潜运动分析及代理模型 3.1 下潜运动分析

首先，在不同的c₁和c₂下开展滑翔机下潜运动的仿真分析，由于滑翔机在Y方向不发生位移，得到滑翔机运动轨迹仿真结果如图 4所示，其中d_T设为400 m。

结果表明，在不同的控制参数配置下，滑翔机下潜至水深400 m，水平位移L存在明显差异，由此可以推断，开展水下定点探测任务时，滑翔机可以采用不同的控制参数配置。考虑海洋环境的多变性以及滑翔机携带的电池包供电能力有限，本文希望滑翔机可以消耗尽量少的能量、尽快到达目标探测区域。因此，有必要确定滑翔机的最优控制参数配置。

3.2 代理模型

若采用优化迭代计算的方式确定最优控制参数值，势必要多次计算E、L和t_d。考虑本文采用数值方法求解动力学模型，需要消耗一定的计算时间，为了减少时间成本，引入代理模型，参与优化迭代计算。

这里，代理模型的输入为控制参数c₁和c₂，输出为E、L和t_d。采用优化拉丁超立方试验设计，在c₁和c₂的取值范围内抽取30组样本点，带入动力学模型开展仿真计算，后续的算例中，d_T设为400 m。仿真完成后，求得30组E、L和t_d。之后，为了确保较高的代理模型拟合精度，采用对非线性模型适应能力较强的完全四阶多项式，对所得数据进行最小二乘拟合，得到E、L和t_d的代理模型，具体表达式如下：

$ \begin{aligned} E \approx & 143-108 c_1-1\ 672 c_2+35.7 c_1^2+10\ 496 c_2^2+ \\ & 717 c_1 c_2-5.35 c_1^3-64\ 366 c_2^3-126 c_1^2 c_2- \\ & 622 c_1 c_2^2+0.308 c_1^4+1350\ 448 c_2^4+ \\ & 7.54 c_1^3 c_2-29\ 219 c_1 c_2^3+147 c_1^2 c_2^2, \end{aligned} $

(28)

$ \begin{aligned} L \approx & 2\ 388-555 c_1-76\ 060 c_2+74.1 c_1^2- \\ & 1269\ 413 c_2^2+18\ 205 c_1 c_2-5.66 c_1^3- \\ & 10\ 836\ 069 c_2^3-1\ 761 c_1^2 c_2-231\ 854 c_1 c_2^2+ \\ & 0.19 c_1^4+33\ 144\ 712 c_2^4-70.09 c_1^3 c_2+ \\ & 1\ 181\ 384 c_1 c_2^3+10\ 674 c_1^2 c_2^2, \end{aligned} $

(29)

$ \begin{aligned} t_{\mathrm{d}} \approx & 65\ 104-50\ 335 c_1-761\ 144 c_2+ \\ & 16\ 233 c_1^2+4\ 746\ 678 c_2^2+325\ 963 c_1 c_2- \\ & 2\ 431 c_1^3-29\ 162\ 978 c_2^3-57\ 121 c_1^2 c_2- \\ & 281\ 350 c_1 c_2^2+140 c_1^4+613\ 739\ 414 c_2^4+ \\ & 3\ 429 c_1^3 c_2-13\ 283\ 197 c_1 c_2^3+66\ 805 c_1^2 c_2^2 . \end{aligned} $

(30)

同时，采用决定系数R²评估代理模型的拟合精度：

$ R^2=1-\frac{\sum\limits_{i=1}^n\left(y_i-y_i^{\mathrm{A}}\right)^2}{\sum\limits_{i=1}^n\left(y_i-y_{\mathrm{ave}}\right)^2} . $

(31)

其中：n表示抽取的样本点总数，y_i表示仿真模型计算的输出值，y_I^A表示代理模型的估计值，y_ave表示所有y_i的平均值。R²越接近1，代理模型的拟合精度越高。代理模型的拟合精度评估结果如图 5—7所示。

结果表明，R²非常接近1，散点全部靠近中间对角线，说明代理模型的拟合精度较高，可以用于后续的优化迭代计算。

基于代理模型，以控制参数为自变量，以E、L和t_d为因变量绘制了函数图像，如图 8—10所示。

结果表明，随着c₁和c₂的取值变化，E、L和t_d具有非线性变化趋势，且E呈现了一定的非单调性，无法根据函数图像直接确定最优控制参数取值。因此，优化迭代计算是非常必要的。此外，通过仿真测试，在计算滑翔机的性能评估参数时，代理模型的计算时长要比动力学模型的计算时长缩短99.8%，计算效率提升明显。

4 控制参数优化计算 4.1 优化数学模型

本文的算例中，设置目标探测点的位置坐标为(700, 0, 400)m。实际探测任务中，滑翔机达到目标点附近区域即可开展定点探测。因此，当滑翔机下潜到目标深度时，L的实际值与理想值之间的差值小于许用值即可。对此，构造优化数学模型的具体表达式如下：

$ \begin{array}{cc} & \min E\left(c_1, c_2\right), \min t_{\mathrm{d}}\left(c_1, c_2\right), \\ \text { s. t. } \quad & 2.5 \mathrm{~N} \leqslant c_1 \leqslant 5 \mathrm{~N}, \\ & 0.01 \mathrm{~m} \leqslant c_2 \leqslant 0.04 \mathrm{~m}, \\ & |L-700 \mathrm{~m}| \leqslant 5 \mathrm{~m} . \end{array} $

(32)

其中：E、L和t_d采用代理模型计算。

显然，该优化问题为有约束的多目标优化问题，本文采用目前流行的多目标优化算法——第二代非劣排序遗传算法(NSGA-Ⅱ)求解上述优化问题。

综上，面向水下定点探测任务，本文提出的水下滑翔机控制参数优化方法的流程如图 11所示。

4.2 优化结果

基于3.2和4.1节，开展优化迭代计算。NSGA-Ⅱ的参数配置选用种群数16、迭代次数60、交叉指数10、变异指数20、交叉概率0.95。

迭代计算完成后，得到了由74组非劣解组成的Pareto最优解集。根据Pareto最优解集，绘制了Pareto前沿，如图 12所示。

Pareto前沿表明，2个优化目标之间具有矛盾关系。之后，从最优解集中选取3组解，开展进一步验证。3组解分别对应最小的E、最小的t_d及2个优化目标折中的解，如表 2所示。

在实际的探测任务中，若希望滑翔机能以较低的能耗到达目标探测区域，则推荐采用第1组解；若希望滑翔机能够快速到达目标探测区域，则推荐采用第2组解；一般情况，推荐采用第3组解，因为其兼顾了能耗和运动时间。之后，以第3组解为例，开展动力学仿真，计算出E、L和t_d的仿真值，并绘制了滑翔机的运动轨迹，如图 13所示。

数值算例证明了所提出方法的正确性，可用于指导实际探测任务的滑翔机控制参数配置。实际应用中，可根据具体型号滑翔机的控制器精度，对优化结果的有效数字进行取舍。

需要注意的是，本文所提出的方法对L具有一定要求，即取值需要在控制参数对应的L最小和最大值之间。虽然具有上述限制，但本文所提出的方法对滑翔机实际探测任务的顺利开展仍具有一定指导意义。

5 结论

基于水下滑翔机的动力学模型，本文面向水下定点探测任务，提出了滑翔机控制参数的优化方法，并给出了相应的数值算例。在不同的控制参数配置下，滑翔机下潜至指定水深，水平位移会存在明显差异，由此可以推断，开展水下定点探测任务时，滑翔机可以采用不同的控制参数配置。在特定的目标下潜深度下，采用完全四阶多项式可以准确描述净浮力调节量、可动质量块平移量、滑翔机能耗、下潜运动时间和水平位移之间的映射关系，且滑翔机的控制参数与性能评估参数之间的函数关系呈现明显的非线性和非单调性。结合代理模型和NSGA-Ⅱ算法完成优化迭代计算，得到了由74组控制参数值非劣解组成的Pareto最优解集。控制参数的优化结果表明，滑翔机的能耗和下潜运动时间两项性能指标之间存在矛盾关系。面向实际工程任务，进一步给出了最优解的选取准则，并采用动力学仿真对优化结果进行了验证。同时给出了所提出方法的局限性，下一步将致力于解决该局限性。