0 引　言

近年来，随着科技的发展和军民需求的推动，越来越多国家制定了海空跨域航行器的研制规划并进行工程化推进。作为一种新兴运载器，海空跨域航行器是飞行器与潜水器技术的融合，将飞行器的高速高机动、快速部署能力与潜水器的游弋能力、高隐蔽性等优势融合在一起，实现了仅用单个航行器便能完成对特定海域同时进行水面、水下等观测的任务目标。比起单一介质的水下或空中航行器，此类航行器不仅能够满足传统海洋科考等方面的需要，还可以满足往常需要多种无人海洋装备联合观测才能满足任务要求，现今已逐渐成为科学家与工程师们广泛关注的新兴焦点。

海空跨域航行器（Unmanned Aerial-Underwater Vehicle，UAUV）这一概念最早来源于1934年，苏联的鲍里斯·乌沙可夫提出的LPL（俄文“会飞的潜艇”的缩写）项目^[1]，但受限于当时的科学技术水平，该项目并未实现。近年来，随着各方面相关技术的迅速发展，UAUV才发展起来。世界上已有多个国家正在积极开展UAUV相关探索和研究，目前尚处于初步探索和样机开发阶段。欧美等发达国家借助在飞行器和潜水器等成熟技术方面的优势，首先将目光聚集在海空两用的航行器零部件探索研究之上，对相关推进原理、海空跨域航行机理与原理样机探究等方面进行开发。2011年，美国麻省理工学院的研究人员设计了一款海空两用的仿生“飞鱼”UAUV原理样机^[2]。2012年，麻省理工学院林肯实验室的研究人员参照鲣鸟收缩翅膀入水捕食的策略，设计了一款小型可重构海空跨域航行器样机^[3]。类似地，麻省理工学院的研究人员于2015年又提出了以串联支架为形式的多模式运动仿生翼^[4]。2015年，哈佛大学的研究人员通过模仿小型昆虫研制了一款仿昆虫的微小型航行器^[5]。2016年，英国帝国理工学院的研究人员通过归纳总结具备空中飞行和水下游泳能力的鸟类和鱼类等生物，设计了一款AquaMAV仿鲣鸟UAUV^{[6 − 7]}。此外很多早期研制成功的UAUV原理样机都是基于多旋翼构型，如美国奥克兰大学的Loon四旋翼UAUV^[8]和美国罗格斯大学的“领航者”八旋翼UAUV^[9]等。2011年开始，国内的中科院沈自所^{[10 − 11]}、北京航空航天大学^{[12 − 13]}等单位也对UAUV进行了一系列的研究。陈怀远^[14]设计了一种融合仿生鱼鳍和涵道螺旋桨的UAUV，并分别对其空中飞行和水下航行的流体动力特性进行分析。Bai等^[15]研制了融合固定翼、多旋翼无人机和水下滑翔机3种不同航行器特点的“哪吒”型UAUV，经过大量湖试、海试，实现了空中飞行与水下航行之间的稳定过渡。李丽云等^[16]研制了“长弓”型折叠翼UAUV，目前已公布了“长弓2号”，采用了中空机翼、辅助水下推进器等方案。

本文结合UAV和UUV两种无人航行器的优势特点，设计了一种可实现水空作业且能多次自由穿越水气界面的新型航行器——折叠翼海空跨域航行器，其兼具了空中飞行时的高机动性和水下航行时的高隐蔽性，具有十分重要的研究前景。首先，本文基于运动学与动力学方程建立UAUV入水砰击过程的数值模型，模拟其入水过程中的受力、运动加速度与速度的变化。然后，进行航行器入水的轨迹优化研究，得出不同任务下的最优入水初始条件，分析航行器入水的运动特性。

1 海空跨域航行器外形设计

飞行器的飞行速度快，可以捕捉和测量空中动态变化的大气现象；而潜水器的机动性能好，能够在水下完成小尺度海洋现象的观测与采样。但空气和水这2种流体的物理性质差别巨大。空中飞行时，需要保持较高的飞行速度，使固定翼产生足够的升力来克服UAUV的重力；而水下潜航时，由于海水密度约为空气的800倍、动力粘性系数约为空气的60倍，航行器受到的流体阻力较大，航行速度相对较低。现有的UAUV主要可分为固定翼和多旋翼两大类，考虑到固定翼飞机更快的速度与更长远的续航力，本文采用固定翼形式。另一方面，出于水下航行器减阻的需要，本文结合空中巡飞弹和水下航行器的特点，提出一款采用折叠翼形式的海空跨域航行器构型设计如图1所示。该UAUV可在空中高速飞行，在水下隐蔽潜航或悬浮待命，且具备多次跨越水气交界面的能力。

UAUV的作业模式如图2所示。该UAUV可在舰船甲板起飞，在空中巡航作业。需要执行水下任务时，UAUV可以折叠机翼直接从空中一头扎入水中，然后在水中航行。在水下潜航期间，UAUV也可悬浮在水中待命。需要执行空中任务时，UAUV即可展开机翼，调整姿态并飞出水面。执行完任务或者在能源快消耗完时，UAUV可以漂浮在水面发送无线电信号求救，等待被回收。此外，UAUV还可作为通讯中继平台，将陆、海、空、天力量统一联系起来，提高海防能力。海空跨域航行器主尺度如表1所示。

由于兼容空中飞行和水下潜航的推进系统难以设计，该UAUV采用首部前拉空气桨与尾部后推水下桨的双推进器布局结构。空气螺旋桨可采用折叠桨的形式，以减小入水时的砰击载荷。海空跨域飞行器的整体密度略小于海水，出水时首部前拉螺旋桨的牵引力可以使得UAUV在出水时高高跃出水面，并迅速转入空中飞行模式。

2 海空跨域航行器动力学建模

从力学角度对刚体进行分析时，需要首先对研究对象进行理论上的分析、建立数学模型，然后进行相关的仿真来对系统进行机理分析。为了确保航行安全，通常使航行器在水平面保持稳定状态，以减小横向运动，因此本文仅研究UAUV的垂直面运动。

2.1 运动学与动力学方程

航行器的动力学模型通常可以分为运动学和动力学2个部分。其中，航行器的垂直面运动学方程为：

$ \left\{ \begin{array}{r@{\;}l} \dot{x} &= u \cos \theta - w \sin \theta ，\\ \dot{z} &= -u \sin \theta + w \cos \theta ，\\ \dot{\theta} &= q 。\end{array} \right. $

(1)

式中：x和z分别为惯性坐标系下航行器的纵向、垂向位移；u和w分别为载体坐标系下航行器的纵向、垂向速度；θ 为航行器的俯仰角；q为航行器的俯仰角速度。惯性坐标系与航行器的载体坐标系如图3 (a)所示。

UAUV在运动时会受到诸多力的作用，具有非线性、多自由度耦合性等特点。同时，由于力的作用点与方向不同，也将分别产生各自的力矩。本文假设UAUV受到的力和力矩作用是可以线性叠加的，其垂直面动力学方程为：

$ \left\{ \begin{aligned} & m(\dot{u} + wq) = -(G-B)\sin\theta + F_{p0} + X_{u}\dot{u} +\\ & \quad X_{|u|u}|u|u + X_{|u|w}|u|w + X_{|u|q}|u|q ，\\ & m(\dot{w} - uq) = (G-B)\cos\theta + Z_{w}\dot{w} + Z_{q}\dot{q} +\\ & \quad Z_{0} + Z_{|u|u}|u|u + Z_{|u|w}|u|w + Z_{|u|q}|u|q ，\\ & I_{yy}\ddot{z} = -(z_{G}G - z_{B}B)\sin\theta + M_{w}\dot{w} + \\ & \quad M_{q}\dot{q} + M_{0} + M_{|u|u}|u|u + M_{|u|w}|u|w + M_{|u|q}|u|q 。\end{aligned} \right. $

(2)

式中：m为航行器质量；I_yy为航行器绕y轴转动惯量；G和B分别为航行器所受的重力和浮力；z_G和z_B分别为载体坐标下重心和浮心的垂向坐标位置；$ X_{{a}} $和$ \boldsymbol{X}_{|\mu| \mu} $等分别为作用于航行器的惯性与粘性水动力系数。

在UAUV入水的过程中，航行器部分在空中、部分浸在水中，所受流体动力与水下航行有很大差别。因此，UAUV的完整入水过程可分为3个阶段：1）俯冲阶段，航行器完全处于空气中，不存在浮力；2）砰击阶段，航行器部分处于水中，其重力大于浮力；3）完全入水，航行器完全浸没于水下，其重力等于浮力。

目前，针对出入水问题的受力时变性与强非线性，其仿真建模方法主要有阶跃法、密度变化法、线性假设法与函数假设法等。整个入水过程中视为无控制输入阶段，UAUV仅受到重力、浮力和流体动力的作用，不存在推进力和操纵力。对于不同成因与不同性质的力，应当采用不同的建模方法以描述各种改变量。由于UAUV的型长远大于其型宽、且本文研究的UAUV入水姿态角小于60°，因此忽略吃水不对称性的影响，假设这些力都仅随入水浸深s的增大而变化。由于UAUV入水时受到的空气动力远小于水动力，可以认为其入水过程中的受的流体作用力主要取决于流体和结构物之间的沾湿面。本节采用流体作用力矩阵的变化率，来研究入水过程中的流体动力的变化速率。由于入水时UAUV浸没在水中的体积排挤液体，会自由产生液面隆起，本章采用等效后的水线WL₁来代替实际的静水面WL。

重力G的大小和定系下的重心位置在整个入水过程中保持不变为定值，而浮力B受自由水面隆起的影响较小。因此，忽略等效自由液面与UAUV纵轴夹角的影响，假设浮力的大小和浮心位置仅取决于入水浸深s变化，如图3 (a)所示。

2.2 惯性类水动力

航行器水下的附加质量矩阵可通过BEM方法进行计算。空气对于航行器的惯性作用力相对于水的惯性作用力较小，可忽略不计。由文献[17]可知，物体在入水过程中的附加质量可认为仅与其入水部分的体积有关，而UAUV入水部分的外形随着入水浸深s的变化而改变。由于时变附加质量矩阵的复杂性，使得求解析解的变化率非常困难，因此本节采用插值的方法来进行计算。

为了确定哪部分航行器外表面与液体接触，也可以将其离散成若干切片。对于三维面元法，UAUV的表面可以用许多块小的四边形面元来近似代替真实的几何形状，所形成的面元是斜高很小的圆台侧面。UAUV的切片计算方法与四边形面元划分如图3 (b)所示。

UAUV在水下的附加质量矩阵如表2所示，其中λ_ij表示附加质量矩阵中的元素。可以看出，UAUV的附加质量矩阵对称性良好，基本呈现对称阵的形式，而且可以发现航行器在z轴方向上的主附加质量与相应回转体平动的理论结果较为接近。因此，可以认为BEM计算得到的附加质量结果均比较精确，能够满足工程需求。

2.3 粘性类水动力

本文基于商业软件STAR-CCM+计算UAUV的粘性类水动力的导数，选取在实尺度下将UAUV在Solidworks软件中建立的三维模型并进行适当简化，保留航行器主体、导流罩、机翼和尾翼等主要结构，仅去除其空气螺旋桨、水下螺旋桨以及舵机等复杂的机构（前文已对螺旋桨的影响单独进行考虑），且忽略重力场。

参考风洞试验方法，直航与斜航的计算域为长方体流体区域，全拘束的UAUV模型在计算域横剖面的中心位置，如图4(a)所示。计算域的大小根据UAUV的主尺度来确定：沿UAUV载体左右两侧分别取5l，上下两侧分别取4l，纵向取16l，其中计算域入口取在首向6l，出口取在尾向8l。因此，建立上述的虚拟数值水池计算域为−8.0l ≤ x ≤6.0l，−5.0l ≤ y ≤ 5.0l和−4.0l ≤ z ≤ 4.0l。

旋转的计算域为绕旋转中心的扇形柱体流体区域，如图4 (b)所示。以模拟UAUV的垂直面旋臂运动为例，令其运动速度为旋转线速度U = qR，其中R为载体动系原点到旋转中心的距离，如图5所示。旋转运动中，保持U不变，通过改变R的大小来改变p。计算域扇形的圆心角为150°，来流截面距离UAUV约60°，周围边界及其他计算设置与直航/斜航相同。

系列的改变计算域速度入口处的来流速度U，来模拟UAUV不同速度的直航工况，分别监测UAUV的流体动力（矩），速度设置分别为0 ~ 5 m/s。对比图6的计算结果可以看出，UAUV载体直航时，所受的阻力、零升力与零升力矩的大小与航速大小的平方近似成正比，与经验计算的公式以及流体动力学理论相符。这是由于粘性流体动力的成因主要是层流边界层的线性摩擦、湍流边界层的二次摩擦和涡流造成的二次阻力。

拘束模型斜航试验是潜水器操纵性试验中不可或缺的重要组成部分，为了便于UAUV模型姿态角的调整，采用重叠网格方法旋转模型，此外设置与直航的计算过程相同。根据斜航试验的经验，在弱机动运动假设下，航行器的运动参数变化δu和δw等都较小，因此可假设u为几乎不变，且通常情况下飞行器或潜水器流体动力的导数与其航速的关系不大^[18]。统一取水下速航为u = 3 m/s。图7 (a)展示了UAUV不同自由度的力与力矩随α 的变化。当α 持续增大时，UAUV受到的垂向力和俯仰力矩分别线性变化。根据流体力学知识，在每个平移自由度上，UAUV匀速运动时受到该运动方向的流体动力为阻碍运动的阻尼力，因此总为负值^[19]。此外，由于UAUV载体外形的前部相对肥大、后部较为瘦削，所以总有M_w<0，这可以从图中得到印证。

保持UAUV的水下巡航速度不变，系列地改变UAUV模型的俯仰角速度q，分别以0 ~ 0.2 rad/s的角速度在计算域中绕其y轴定常回转，监测UAUV所受的力（矩）。其中，旋转角速度为0即为UAUV以相同航速定常直航的工况。图7 (b)展示了UAUV不同自由度的力与力矩随q的变化，可以看到旋臂计算结果的线性程度也比较高，可以认为UAUV在计算范围内符合小机动假设。对以上计算所得的数据进行最小二乘法处理后，可以拟合求得关于线速度的一系列流体动力系数，结果如表3所示。

考虑到入水过程中粘性流体的复杂受力，此处综合采用密度变化法与线性假设法对入水的粘性流体动力系数进行建模。初步假设在近水面区域的一个过渡区：当UAUV全长的l /2部分进入水中之后（即等效水线WL₁开始浸没至细长碳管机身），粘性流体作用开始产生；当UAUV浸没水中的深度为其全长的平方l ²时，粘性流体作用发展充分，此时水动力系数完全等于UUV模式。当s = l/2 ~ l ²时，粘性水动力系数随s线性变化。以$ X_{|u|u} $为例，计算公式如下：

$ X_{|u|u}(s)=\left\{\begin{array}{l} 0, s < l / 2，\\ \dfrac{2 s-l}{2 l^{2}-l} X_{|u|u} , \dfrac{l}{2} \leqslant s < l^{2}，\\ X_{|u|u}, s \geqslant l^{2}。\end{array}\right. $

(3)

2.4 入水砰击力

从数学的角度考虑，入水砰击力的数值具有强非线性与衰减性。通常在数学建模中，常用的衰减拟合函数有指数函数等。根据文献[17]中对NACA机翼入水的砰击力结果与文献[20]中对UAUV砰击力系数进行分析，本节选用倒数函数来模拟入水砰击力的衰减。以入水时刻（位置）为0，认为入水过程中砰击力系数的倒数随t（即s）呈线性变化，如图8所示。假设UAUV在浸深为s₀处的砰击力达到峰值，并假设入水砰击力的方向总与UAUV的运动速度方向一致。航行器入水过程中的砰击力系数C_p0(s)估算式为：

$ C_{p 0}(s)=\frac{1}{\left(\dfrac{1}{C_{p 0}^{\max }}+\gamma \left(s-s_{0}\right)\right)} 。$

(4)

式中：$C_{p 0}^{\max } $为砰击力峰值系数；s₀为砰击力达到峰值处对应的浸深。

因此，航行器入水过程中的砰击力F_p0为：

$ F_{p 0}=\frac{1}{2} \rho C_{p 0}(s) u^{2} l^{2} 。$

(5)

3 入水轨迹优化结果

为了避免复杂控制策略的切换造成的不稳定，将UAUV的自由入水过程视作初值问题，即入水轨迹完全由入水前初始时刻的运动状态确定。本节以入水攻角和俯仰姿态为优化变量，对不同任务要求下的入水轨迹进行优化。计算某一任务对应的指标函数值时，需要将初始条件代入UAUV入水仿真的运动学与动力学方程中进行计算。对UAUV初始变量即初始姿态角θ _s与攻角α _s（为了改变速度V_s）进行优化，在可行域内找到最优解使得指标函数值最小。为了解决这一轨迹优化问题，利用Matlab中用于解决非确定性多项式（Nondeterministic Polynomial，NP）问题的fmincon函数来完成优化求解。其中，本文涉及到该函数返回值ef（Exit Flag）含义为：ef = 1表示满足一阶最优性条件，ef = 2表示变量小于函数的容许范围。

以最短时间完成入水与最接近拟定目标点2个任务为例，采用变积分步长的4阶Runge-Kutta方法进行求解，取初始时间步长dt = 0.0005 s。仿真过程中，dt的计算原则是在每个时间步长之内位移的变化大小不超过基础尺寸dl。

3.1 最短时间完成入水

任务1：UAUV最速完成入水，即在最短时间内到达水下一定深度处，此时的优化指标函数J₁为：

$ J_{1}=\int_{t_{0}}^{t_{1}} {\mathrm{d}} t 。$

(6)

式中：t_s和t_f分别为航行器入水过程起始和结束的时间。由于UAUV的深度信息由深度计获取，以航行器在空中从某一固定高度h、以一定初速度V_s自由释放点作为入水任务起始时刻t_s，以航行器入水后深度传感器测得的深度h大于目标深度值ζ_f作为结束时刻t_f。通常来说，入水时间越短，越有利于航行器的隐身性安全与其跨介质运动的顺利切换。

由于初始速度越大则入水过程越快，因此任务1中保持初始速度的大小V_s = 12 m/s不变。其中，V_s = |(u_s, v_s)|。取初始高度h = 0.5 m，搜索空间为UAUV的初始俯仰角θ_s与攻角α_s，其他初始量皆为0，以J₁（即t_f）作为评价指标，找到θ _s与α _s的最优组合，使其满足任务一的需求。

系列的改变目标深度ζ₁，得到任务1的优化结果如表4所示，其中θ _s与α _s为最优初始变量。任务1的优化指标函数图像如图9 (a)所示，可见当UAUV入水的初始姿态角在可行目标域的范围内变化时，此处所求得的最优解可认为是全局最优解而非局部最优。从表4可知，在各目标深度下，最优的α _s均为0、θ _s均为−85.5°。由此可知，不同的ζ₁对任务1的优化结果几乎没有影响。

3.2 最接近拟定目标点

任务2：UAUV最接近拟定目标点，即入水过程完成时UAUV的位置距离目标点的位置误差最小，此时的优化指标J₂为：

$ J_2=\sqrt{\left(\xi\left(t_f\right)-\xi_2\right)^2+\left(\zeta\left(t_f\right)-\zeta_2\right)^2}。$

(7)

式中：(ξ₂, ζ₂)为拟定目标点；(ξ(t_f), ζ(t_f))为最后时刻UAUV的位置。

取拟定目标点(ξ₂, ζ₂)为(0, 4)，即要求UAUV入水路径终点的水平位移最小。保持初始高度h = 0.5 m不变，搜索空间仍为UAUV的初始俯仰角θ _s与攻角α_s，其他初始量皆为0，以J₂（即|ξ(t_f)|）作为评价指标，找到θ _s与α _s的最优组合，使其满足任务2的需求。

系列的改变V_s，得到任务2的优化结果如表5所示，其中θ _s与α _s为最优初始变量。任务2的优化指标函数图像如图9 (b)所示，同样可以认为所求得的最优解是全局最优。可知，θ_s与α_s均对J₂有影响，且θ_s对J₂的影响更大。可知在各初始条件V_s下，最优的θ_s均接近垂直入水，而α_s始终约为2.04°。由于本文UAUV在动系下的重心略低于浮心且外形关于水线面并不是完全对称，所以存在横向受力，因此需要具备一定的攻角初始值α_s来弥补其入水运动的不对称性。

任务2所得的部分仿真时历曲线如图10所示。由图10 (a)可知，在UAUV达到4 m水深处时，其攻角减小至−0.4°以下，纵向航速减小至2 m/s左右，随后可以开启螺旋桨实现螺旋下潜。从图10 (b)中还可以发现，不同工况的俯仰角变化曲线都在约s = 3.5 m，θ = 90°点处相交，且交点左右曲线与横轴所围的面积之比近似相等，这同样可见UAUV入水运动的不对称性。

在任务2所要求的达到满足最接近指定目标点的基础上，后续还可加入入水时间约束，只需在优化指标函数中加入时间项考虑权重系数，即可兼顾上述2个任务的要求。另外，在以上2个仿真过程中发现，由于ef条件的不同，求解器完成任务一的计算量明显小于任务2。这可能是因为在UAUV的入水运动中，纵向速度u越大时入水时间越短，由此任务1等价于只需对俯仰角θ 单变量进行优化、搜索空间维数为一维，故收敛速度较快。综上所述，任务1与任务2是分别针对UAUV几何外形不对称程度与运动不对称程度的轨迹优化问题。

4 结　语

本文提出一种可实现水空两栖作业且能多次自由穿越水气界面的航行器—海空跨域航行器。对UAUV入水过程可控变量进行优化，得到了满足不同任务需求的入水初始状态值，保证其精准、高效完成入水过程。通过计算分析得出如下结论：

1）UAUV入水轨迹取决于入水开始时的初始状态，不同目标深度所对应最速入水的最优初始攻角均为0、最优初始俯仰角也相同；

2）不同速度所对应最接近拟定目标点的最优初始攻角不为0，而最优初始俯仰角接近垂直状态。