可变转动惯量的两栖无人机控制

Control of amphibious unmanned aerial vehicles （UAVs） with varying moments of inertia

摘要

针对两栖无人机在飞行过程中因机械腿运动引起的转动惯量变化对飞行稳定性的影响，开发了一种适应性前馈控制器。首先构建了六旋翼无人机的动力学模型；随后基于机械腿运动角度的可变范围，确定了转动惯量的变化区间，并设计了一套能够覆盖最大、最小及变化极值的机械腿运动路径，该路径将转动惯量与时间的关系转化为非线性函数；进一步，针对两栖无人机的转动惯量变化特性，研制了一种基于前馈控制的控制策略；最后，在Simulink环境中对转动惯量动态变化条件下的前馈控制与传统的比例-积分-微分（PID）控制进行仿真比较，验证了所设计前馈控制器在处理两栖无人机转动惯量变化时的优越性能。

Abstract

An adaptive feedforward controller has been developed to address the impact on flight stability of variations in the moment of inertia caused by mechanical leg movements during the flight of amphibious unmanned aerial vehicles （UAVs）. Firstly， a dynamic model of a hexacopter UAV was established. Based on the variable range of the motion angles of the mechanical legs， the variation interval of the moment of inertia was determined， and a set of mechanical leg motion paths covering the maximum， minimum and extreme variation values of moment of inertia was designed. This converts the relationship between the moment of inertia and the time into a nonlinear function. Furthermore， given the variable moment of inertia of amphibious UAVs， a control strategy based on feedforward control was developed. Finally， a simulation in Simulink was carried out to compare feedforward control with the traditional proportional-integral-derivative （PID） control under the conditions of a dynamically changing moment of inertia. The results verified the superior performance of our feedforward controller in dealing with variations in the moment of inertia of amphibious UAVs.

Graphical abstract

关键词

两栖无人机 / 变转动惯量 / 前馈控制 / 六旋翼

Key words

amphibious UAV / variable moment of inertia / feedforward control / hexacopter

引言

随着无人机技术的迅猛发展，陆空两栖无人机因其出色的环境适应性和多样化任务执行能力而受到广泛关注。这些多功能无人机不仅能够在空中飞行，还能在陆地上移动，特别适合于湿地和森林环境监测，以及自然灾害后的搜救行动。为适应复杂多变的环境，陆空两栖无人机发展出多种结构形式，其中飞行部分以多旋翼为主，陆地移动则以轮式和履带式为常见选择。

为了提升无人机的性能，研究人员提出了多种控制策略。在传统比例-积分-微分（PID）控制方面，Canal等^［1］提出一种适应不同载荷变化的PID控制策略，以低成本和易于实施的特点通过实验验证了其在无人机上的有效性；何乃晨^［2］设计了一种两输入十输出的模糊PID控制器，以提高无人机的抗干扰能力；张长箭等^［3］采用模糊PID控制方法，针对传统PID控制响应时间长和动态性能差的问题，设计了四旋翼无人机动力学模型，有效提升了系统的稳定性和动态性能；Goje等^［4］针对无人机的稳定性问题提出了利用PID控制器实现稳定飞行的方法，通过机器学习多项式回归模型调整PID控制系统中的增益因子，以减少小的不稳定因素；他们通过在Gazebo模拟器和机器人操作系统（ROS）环境中的模拟测试，成功展示了其提出的2D路径规划算法；实验结果也验证了该方法能够在减少环境和系统干扰的同时，实现无人机的平滑无振荡飞行。在鲁棒抗扰控制技术方面，薛晶勇等^［5］针对无人机在面临外部不确定干 扰及系统参数不确定性时的控制问题，提出了一种结合不确定干扰前馈补偿与反步滑模的姿态控制方法；利用非线性干扰观测器估计干扰，设计前馈补偿反步滑模控制器，并经Lyapunov 理论验证稳定性，有效应对高阶不确定性。Yang等^［6］通过基于牛顿-欧拉公式的数学建模和二阶滑模扰动观测器（SOSMDO）重建旋转动力学的扰动，设计了一种嵌套双环控制策略，包括外环的位置误差PID控制器和内环的基于扰动补偿和快速收敛非线性微分器（RCND）的鲁棒姿态控制器，并通过Matlab/Simulink仿真验证了其有效性和鲁棒性。徐雪松等^［7］针对负载变化对四旋翼无人机飞行稳定性的影响，利用模糊控制在线调整非线性状态误差反馈控制律参数，设计了基于自适应变载荷补偿的模糊自抗扰控制系统，解决了自抗扰控制器参数自整定问题。在容错控制与场景适配方面，贺有智等^［8］提出了一种四旋翼飞行器的时延积分反演容错控制方法，旨在解决飞行过程中可能遇到的执行机构故障问题；该方法融合积分反演与时延控制技术，无需在线故障检测即可通过一步状态迭代补偿执行机构故障，提升鲁棒性。李泽宇^［9］设计了一种两栖变结构无人机，对比了互补滤波和梯度下降的姿态解算性能，以及串级PID和串级自抗扰的姿态控制性能，最终选择了互补滤波进行位置解算和串级PID进行位置控制。Giernacki等^［10］提出了一种混合准最优比例积分微分-状态相关黎卡提方程（PID⁃SDRE）控制方法，通过结合传统的PID控制和有限时间SDRE控制技术，在保持控制系统简单性的同时显著提升了控制质量。

然而，现有研究多关注于抵抗外部干扰，对于转动惯量大幅变化的无人机控制研究较少。本文设计中的无人机在陆地移动部分采用六足结构，显著增强了越障能力；在飞行模式下，其机械腿可执行多种动作以完成不同任务。针对两栖无人机飞行时机械腿运动引起的转动惯量非线性变化进行研究，通过平移轴定理建立了机械腿关节角度与转动惯量之间的关系，并计算出对应转动惯量最大值和最小值的关节角度。基于此，设计了覆盖转动惯量变化范围的机械腿运动轨迹，进一步针对转动惯量大幅变化的特性，设计了前馈控制策略，有效减少了转动惯量变化对位置环和姿态环的影响。与传统PID控制相比，本设计展现出更优的动态性能。

1 六旋翼建模

本文所研究的两栖无人机由六旋翼无人机和6条机械腿组成，两栖无人机机体为正六边形结构，6条机械腿均匀地布置在正六边形机体的6个角上，结构如图1所示。然后建立坐标系，o_B 点位于六旋翼飞行器的质心，x_B 指向旋翼5，z_B 轴垂直于旋翼平面向下，y_B 轴垂直于o_B，x_B 、z_B 平面向前，符合右手定则。（o_B，x_B，y_B，z_B ）组成的坐标系表示机体坐标系｛B｝，用来描述无人机的滚转角φ，俯仰角\theta，偏航角\psi，和无人机绕o_Bx_B、o_By_B、o_Bz_B 旋转的角速度\dot{\phi }、\dot{\theta }、\dot{\psi }。（o_E，x_E，y_E，z_E ）组成的坐标系表示地面坐标系｛E｝，用来说明无人机的位置，如图2所示。

无人机6个旋翼的转速分别定义为{\omega }_{\mathrm{1}}、{\omega }_{\mathrm{2}}、{\omega }_{\mathrm{3}}、{\omega }_{\mathrm{4}}、{\omega }_{\mathrm{5}}、{\omega }_{\mathrm{6}}，6个旋翼同时转动会产生向上的拉力U_L，无人机开始产生向上的运动趋势。当旋翼1、2、3转速增加或减少，旋翼4、5、6转速减少或增加时，无人机产生绕y_B轴旋转的力矩{\tau }_y，滚转角开始变化；当旋翼2、5的转速不变，旋翼1、6和旋翼3、4的转速分别增大和减小时，产生绕x_B轴旋转的力矩{\tau }_x，俯仰角开始变化；当旋翼1、3、5转速增大，旋翼2、4、6转速减少，无人机开始产生绕z_B轴旋转的力矩{\tau }_z，偏航角开始改变^［3］，则可以得到

\left\{\begin{array}{l}{U}_L=c_T({\omega }_{\mathrm{1}}^{\mathrm{2}}+{\omega }_{\mathrm{2}}^{\mathrm{2}}+{\omega }_{\mathrm{3}}^{\mathrm{2}}+{\omega }_{\mathrm{4}}^{\mathrm{2}}+{\omega }_{\mathrm{5}}^{\mathrm{2}}+{\omega }_{\mathrm{6}}^{\mathrm{2}})\\ {\tau }_x=\frac{\sqrt[]{\mathrm{3}}}{\mathrm{2}}d{c}_T(-{\omega }_{\mathrm{1}}^{\mathrm{2}}+{\omega }_{\mathrm{3}}^{\mathrm{2}}+{\omega }_{\mathrm{4}}^{\mathrm{2}}-{\omega }_{\mathrm{6}}^{\mathrm{2}})\\ {\tau }_y=\frac{\mathrm{1}}{\mathrm{2}}d{c}_T(-{\omega }_{\mathrm{1}}^{\mathrm{2}}-\mathrm{2}{\omega }_{\mathrm{2}}^{\mathrm{2}}-{\omega }_{\mathrm{3}}^{\mathrm{2}}+{\omega }_{\mathrm{4}}^{\mathrm{2}}+\mathrm{2}{\omega }_{\mathrm{5}}^{\mathrm{2}}+{\omega }_{\mathrm{6}}^{\mathrm{2}})\\ {\tau }_z=c_M({\omega }_{\mathrm{1}}^{\mathrm{2}}-{\omega }_{\mathrm{2}}^{\mathrm{2}}+{\omega }_{\mathrm{3}}^{\mathrm{2}}-{\omega }_{\mathrm{4}}^{\mathrm{2}}+{\omega }_{\mathrm{5}}^{\mathrm{2}}-{\omega }_{\mathrm{6}}^{\mathrm{2}})\end{array}\right.（1）

式中，c_T为升力系数，c_M为反扭矩系数，d为轴距的一半。假设六旋翼为刚体，则可以得到六自由度动力学模型为

\left\{\begin{array}{l}\ddot{x}=\frac{\mathrm{1}}{m}(a_{ux}-k_x\dot{x}+f_x)\\ \ddot{y}=\frac{\mathrm{1}}{m}(a_{uy}-k_y\dot{y}+f_y)\\ \ddot{z}=\frac{\mathrm{1}}{m}(a_{uz}-k_z\dot{z}+f_z)-g\\ \ddot{\phi }=\frac{\mathrm{1}}{J_{xx}}\left((J_{yy}-J_{zz})\dot{\theta }\dot{\psi }+{\tau }_x-k_{\phi }\dot{\phi }+{\tau }_{f\phi }\right)\\ \ddot{\theta }=\frac{\mathrm{1}}{J_{yy}}\left((J_{zz}-J_{xx})\dot{\phi }\dot{\psi }+{\tau }_y-k_{\theta }\dot{\theta }+{\tau }_{f\theta }\right)\\ \ddot{\psi }=\frac{\mathrm{1}}{J_{zz}}\left((J_{xx}-J_{yy})\dot{\phi }\dot{\theta }+{\tau }_z-k_{\psi }\dot{\psi }+{\tau }_{f\psi }\right)\end{array}\right.（2）

式中k_x、k_y、k_z、k_{\phi }、k_{\theta }、k_{\psi }为阻力系数，f_x、f_y、f_z为外部干扰力，{\tau }_{f\theta }、{\tau }_{f\phi }、{\tau }_{f\psi }为外部干扰力矩，a_{ux}、a_{uy}、a_{uz}为位置回路的虚拟控制量，计算式为

\left\{\begin{array}{l}{a}_{ux}=U_L(\mathrm{c}\mathrm{o}\mathrm{s}\phi \mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{c}\mathrm{o}\mathrm{s}\psi +\mathrm{s}\mathrm{i}\mathrm{n}\phi \mathrm{s}\mathrm{i}\mathrm{n}\psi )\\ a_{uy}=U_L(\mathrm{c}\mathrm{o}\mathrm{s}\phi \mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{s}\mathrm{i}\mathrm{n}\psi +\mathrm{s}\mathrm{i}\mathrm{n}\phi \mathrm{c}\mathrm{o}\mathrm{s}\psi )\\ a_{uz}=U_L\mathrm{c}\mathrm{o}\mathrm{s}\phi \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\right.（3）

2 机械腿运动轨迹设计

两栖无人机的机械腿由根关节、髋关节、膝关节3个关节和基节、股节、胫节3个腿部关节连杆组件组成。根关节将机体和基节连接在一起，根关节的转轴平行于z_B轴，使基节能在x_Bo_By_B 组成的水平面内转动，髋关节将基节和股节连接在一起，膝关节将股节和胫节连接在一起，髋关节和膝关节的转轴相互平行且垂直于根关节转轴，使股节和胫节能在竖直面内旋转^［11-12］。以机械腿的基节、股节、胫节均在x_Bo_By_B 平面内，且都平行于x_B 轴所对应的关节角为各关节的0°位置，根关节的转角范围为［-30°，30°］，髋关节的转角范围为［-90°，0°］，膝关节的转角范围为［-90°，90°］，如图3所示。

将根关节电机与机体视为一个整体，将髋关节电机和膝关节电机与股节视为一个整体，则可得到机体和机械腿各部分在各自质心处的转动惯量，如表1所示。机械腿各部分的长度和质量如表2所示。

根据平移轴定理可推算得到机械腿各部分在坐标系｛B｝下的转动惯量公式为

\left\{\begin{array}{l}{J}_{ix}^{\mathrm{\text{'}}}=J_{ix}+m_{\mathrm{c}}\left(\left(l_B^{}\left|\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{2}i\mathrm{\pi }-\mathrm{\pi }}{\mathrm{6}}\right)\right|\right.\right.+\\ l_{\mathrm{c}z}{\left.\left|\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{c}}\right|\right)}^{\mathrm{2}}+\left.l_h^{\mathrm{2}}\right)+\\ m_{\mathrm{t}}\left(\left(l_B^{}\left|\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{2}i\mathrm{\pi }-\mathrm{\pi }}{\mathrm{6}}\right)\right|\right.\right.+(l_{\mathrm{c}}+\\ l_{\mathrm{t}z}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{t}}){\left.\left|\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{c}}\right|\right)}^{\mathrm{2}}+\left.l_h^{\mathrm{2}}\right)+\\ m_{\mathrm{s}}\left(\left(l_B\left|\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{2}i\mathrm{\pi }-\mathrm{\pi }}{\mathrm{6}}\right)\right|\right.\right.+(l_{\mathrm{c}}+(l_{\mathrm{t}}+\\ l_{\mathrm{s}z}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{t}}){\left.\left|\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{c}}\right|\right)}^{\mathrm{2}}+\left.l_h^{\mathrm{2}}\right)\\ J_{iy}^{\mathrm{\text{'}}}=J_{iy}+m_{\mathrm{c}}\left(\left(l_B\left|\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\mathrm{2}i\mathrm{\pi }-\mathrm{\pi }}{\mathrm{6}}\right)\right|\right.\right.+\\ l_{\mathrm{c}z}{\left.\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{c}}\right)}^{\mathrm{2}}+\left.l_h^{\mathrm{2}}\right)+m_{\mathrm{t}}\left(\left(l_B^{}\left|\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\mathrm{2}i\mathrm{\pi }-\mathrm{\pi }}{\mathrm{6}}\right)\right|\right.\right.+\\ (l_{\mathrm{c}}+l_{\mathrm{t}z}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{t}}){\left.\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{c}}\right)}^{\mathrm{2}}+\left.l_h^{\mathrm{2}}\right)+\\ m_{\mathrm{s}}\left(\left(l_B\left|\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\mathrm{2}i\mathrm{\pi }-\mathrm{\pi }}{\mathrm{6}}\right)\right|\right.\right.+(l_{\mathrm{c}}+(l_{\mathrm{t}}+\\ l_{\mathrm{s}z}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{t}}){\left.\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{c}}\right)}^{\mathrm{2}}+\left.l_h^{\mathrm{2}}\right)\\ J_{iz}^{\mathrm{\text{'}}}=J_{iz}+J_{ix}^{\mathrm{\text{'}}}+J_{iy}^{\mathrm{\text{'}}}-J_{ix}-J_{iy}-\mathrm{2}{l}_h^{\mathrm{2}}(m_{\mathrm{c}}+m_{\mathrm{t}}+m_{\mathrm{s}})\end{array}\right.（4）

式中，i表示第i（i=1、2、3、4、5、6）条机械腿，位置与旋翼相同，J_{ix}、J_{iy}、J_{iz}表示第i条机械腿的各部分在自身质心处的x方向、y方向、z方向的转动惯量总和，J_{ix}^{\mathrm{\text{'}}}、J_{iy}^{\mathrm{\text{'}}}、J_{iz}^{\mathrm{\text{'}}}表示第i条腿在坐标系｛B｝中的转动惯量，m_{\mathrm{c}}、m_{\mathrm{t}}、m_{\mathrm{s}}分别表示基节、股节、胫节的质量，l_{\mathrm{c}}、l_{\mathrm{t}}、l_{\mathrm{s}}分别表示基节、股节、胫节的长度，l_B=0.295 m表示z_B轴到根关节的距离，6条机械腿的l_B相同，l_{\mathrm{c}z}、l_{\mathrm{t}z}、l_{\mathrm{s}z}分别表示基节、股节、胫节的质心到前一关节的距离，{\alpha }_{\mathrm{c}}、{\alpha }_{\mathrm{t}}、{\alpha }_{\mathrm{s}}分别表示根关节、髋关节、膝关节的旋转角度，l_h=0.085 m表示根关节到x_Bo_By_B 平面的距离。

整个两栖无人机在坐标系｛B｝下的转动惯量表示为

\left\{\begin{array}{l}{J}_{xx}=J_{Bx}+{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{6}}}{J}_{ix}^{\mathrm{\text{'}}}\\ J_{yy}=J_{By}+{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{6}}}{J}_{iy}^{\mathrm{\text{'}}}\\ J_{zz}=J_{Bz}+{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{6}}}{J}_{iz}^{\mathrm{\text{'}}}\end{array}\right.（5）

式中，J_{xx}、J_{yy}、J_{zz}分别表示两栖无人机在x_B、y_B、z_B 轴的转动惯量，J_{Bx}、J_{By}、J_{Bz}分别表示六旋翼机体在x_B、y_B、z_B 轴的转动惯量。

为找出两栖无人机转动惯量的最大值和最小值，设定5条机械腿的运动轨迹：第一条运动轨迹为髋关节保持90°，膝关节保持0°，根关节以60（°）/s的角速度从-30°转到30°；第二条运动轨迹为髋关节保持0°，膝关节保持-90°，根关节以60（°）/s的角速度从-30°转到30°；第三条运动轨迹为髋关节保持0°，膝关节保持0°，根关节以60（°）/s的角速度从-30°转到30°；第四条运动轨迹为髋关节保持0°，膝关节以90（°）/s的角速度从-90°转到0°，根关节以60（°）/s的角速度从-30°转到30°；第五条运动轨迹为膝关节保持0°，髋关节以90（°）/s的角速度从-90°转到0°，根关节以60（°）/s的角速度从-30°转到30°。将运动轨迹代入到式（4）、（5中计算出六旋翼无人机转动惯量的变化情况，计算结果如图4所示。将每一条轨迹中的转动惯量最大值和最小值及其对应的关节角统计出来，可得到转动惯量的最大值和最小值对应的关节角，如表3所示。

由表3可以看出，5条轨迹中最高的J_{xx}的最大值对应的关节角为（±30°，0°，0°）（其中第一项为α_c，第二项为α_t，第三项为α_s，下同），5条轨迹中最低的J_xx 最小值对应的关节角为{\alpha }_{\mathrm{c}}=\mathrm{0}°、{\alpha }_{\mathrm{t}}和{\alpha }_{\mathrm{s}}取范围内任意值；5条轨迹中最高的J_{yy}的最大值对应的关节角为（0°，0°，0°），最低的J_yy 最小值对应的关节角为（±30°，-90°，0°）；5条轨迹中最高的J_{zz}的最大值对应的关节角为（±30°，0°，0°），最低的J_zz 最小值对应的关节角为（0°，90°，0°）。依此所设计的供测试控制策略使用的运动轨迹分为4段：第一段髋关节以90（°）/s的角速度从（-30°，-90°，0°）位置转到（-30°，0°，0°）；第二段根关节以30（°）/s的角速度从（-30°，0°，0°）转到（0°，0°，0°）；第三段根关节以30（°）/s、膝关节以90（°）/s的速度从（0°，0°，0°）转到（30°，-90°，0°）；第四段根关节以60（°）/s的角速度从（30°，-90°，0°）转回起点（-30°，-90°，0°）。运动轨迹如图5所示。

其中第一段J_{zz}从最小值变为最大值，第二段J_{xx}从最大值变为最小值，第三段J_{yy}从最大值变为最小值，第四段回到起点。每段内关节匀速转动，4 s为一个周期。将这条运动轨迹中的关节转角写成时间的函数，表示为

（6）

将式（6）代入式（4））、（5）中，得到两栖无人机转动惯量随时间变化的函数，将其在Simulink中进行20 s的仿真，得到这段轨迹中无人机转动惯量随时间变化的曲线，如图6所示。

3 前馈控制设计

由于两栖六旋翼无人机实际是一种非线性、耦合和欠驱动系统，具有六自由度，但只有4个实际输入的系统，同时该系统的转动惯量也随时间不断变化，因此需要有一种既能跟踪位置和姿态角，又能适应时变的转动惯量的控制器^［10］。

控制器通过控制电机的转速来控制两栖无人机的位姿，由式（1）可知，控制U_L、{\tau }_x、{\tau }_y、{\tau }_z的值就能够控制电机转速，又由式（2）和（3）可知，控制U_L、{\tau }_x、{\tau }_y、{\tau }_z能够控制无人机的位姿，因此本设计通过调控U_L、{\tau }_x、{\tau }_y、{\tau }_z以控制电机转速，进而控制无人机位姿来实现对六旋翼无人机的控制，控制器结构如图7所示。

3.1 位姿控制律设计

将U_L、{\tau }_x、{\tau }_y、{\tau }_z作为控制律，针对3个位置子系统和3个姿态子系统设计控制器，实现x\to x_d、y\to y_d、z\to z_d、\phi \to {\phi }_d、\theta \to {\theta }_d、\psi \to {\psi }_d。由式（1）可得位置控制模型为

\left\{\begin{array}{l}\ddot{x}=\frac{a_{ux}}{m}-\frac{k_x\dot{x}}{m}\\ \ddot{y}=\frac{a_{uy}}{m}-\frac{k_y\dot{y}}{m}\\ \ddot{z}=\frac{a_{uz}}{m}-g-\frac{k_z\dot{z}}{m}\end{array}\right.（7）

姿态控制模型为

\left\{\begin{array}{l}\ddot{\phi }=\frac{-k_{\phi }\dot{\phi }}{J_{xx}}+{\tau }_x\\ \ddot{\theta }=\frac{-k_{\theta }\dot{\theta }}{J_{yy}}+{\tau }_y\\ \ddot{\psi }=\frac{-k_{\psi }\dot{\psi }}{J_{zz}}+{\tau }_z\end{array}\right.（8）

设x_d、y_d、z_d、{\phi }_d、{\theta }_d、{\psi }_d为期望值，x、y、z、\phi、\theta、\psi为实际值，x_e、y_e、z_e、{\phi }_e、{\theta }_e、{\psi }_e为实际值减去期望值得到的误差值，采用比例-微分（PD）控制方法得到控制律为

\left\{\begin{array}{l}{a}_{ux}=-K_{px}{x}_e-K_{dx}{\dot{x}}_e-K_{ix}{\int }_{\mathrm{0}}^t{x}_e\mathrm{d}t\\ a_{uy}=-K_{py}{y}_e-K_{dy}{\dot{y}}_e-K_{iy}{\int }_{\mathrm{0}}^t{y}_e\mathrm{d}t\\ a_{uz}=-K_{pz}{z}_e-K_{dz}{\dot{z}}_e-K_{iz}{\int }_{\mathrm{0}}^t{z}_e\mathrm{d}t+g+{\ddot{z}}_d+\frac{k_z{\dot{z}}_d}{m}\\ {\tau }_x=-K_{p\phi }{\phi }_e-K_{d\phi }{\dot{\phi }}_e+{\ddot{\phi }}_d+\frac{k_{\phi }\dot{\phi }}{J_{xx}}\\ {\tau }_y=-K_{p\theta }{\theta }_e-K_{d\theta }{\dot{\theta }}_e+{\ddot{\theta }}_d+\frac{k_{\theta }{\dot{\theta }}_d}{J_{yy}}\\ {\tau }_z=-K_{p\psi }{\psi }_e-K_{d\psi }{\dot{\psi }}_e+{\ddot{\psi }}_d+\frac{k_{\psi }{\dot{\psi }}_d}{J_{zz}}\end{array}\right.（9）

式中，K_px、K_dx、K_py、K_dy、K_pz、K_dz 为位置环参数，K_pψ 、K_dψ 、K_pθ 、K_dθ 、K_pφ 、K_dφ 为姿态环参数。根据Hurwitz判据，二阶系统的比例系数、微分系数和积分系数均大于0即可保证位置控制律和姿态控制律系统稳定^［4］。

3.2 期望姿态角控制律设计

为了使\theta \to {\theta }_d，\phi \to {\phi }_d，需要求解这两个姿态角。由式（3）中的第三式可得

U_L=\frac{a_{uz}}{\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_d\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_d}（10）

将式（10）代入式（3）的前两式可得

{\phi }_d=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\mathrm{s}\mathrm{i}\mathrm{n}\psi \mathrm{c}\mathrm{o}\mathrm{s}\psi a_{ux}-\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\psi a_{uy}}{a_{uz}}\right)（11）

{\theta }_d=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\mathrm{c}\mathrm{o}\mathrm{s}\psi (\mathrm{s}\mathrm{i}\mathrm{n}\psi a_{ux}-\mathrm{c}\mathrm{o}\mathrm{s}\psi a_{uy})}{a_{uz}}\right)（12）

令E=\frac{\mathrm{c}\mathrm{o}\mathrm{s}\psi (\mathrm{s}\mathrm{i}\mathrm{n}\psi a_{ux}-\mathrm{c}\mathrm{o}\mathrm{s}\psi a_{uy})}{a_{uz}}，由式（12）可知，只有当E的值在［-1，1］内{\theta }_d才存在，则当E>1时，取{\theta }_d=\frac{\mathrm{\pi }}{\mathrm{2}}；当E<-1时，取{\theta }_d=-\frac{\mathrm{\pi }}{\mathrm{2}}；当\left|E\right|\le \mathrm{1}时，取{\theta }_d=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}E。将得到的{\phi }_d和{\theta }_d代入式（10）可计算出位置控制律U_L。

3.3 期望姿态角速度解算器

本设计使用有限时间三阶微分器进行角速度求解。传统微分器的工作原理基于微分的数学定义，即对函数进行局部线性化和差分运算。然而这种直接微分的方法十分容易受到干扰，通过引入低通滤波环节，可有效抑制干扰的放大效应，并提高微分信号的质量，将有限时间三阶微分器安排在系统的过渡过程中，有效减少系统的超调和振荡，从而提高系统响应的品质。解算器为

\left\{\begin{array}{l}{\dot{x}}_{\mathrm{1}}=x_{\mathrm{2}}\\ {\dot{x}}_{\mathrm{2}}=x_{\mathrm{3}}\\ {\xi }^{\mathrm{3}}{\dot{x}}_{\mathrm{3}}=-{\mathrm{2}}^{\frac{\mathrm{3}}{\mathrm{5}}}\cdot \mathrm{4}(x_{\mathrm{1}}-v(t)+(\xi x_{\mathrm{2}}{)}^{\frac{\mathrm{9}}{\mathrm{7}}}{)}^{\frac{\mathrm{1}}{\mathrm{3}}}-\mathrm{4}({\xi }^{\mathrm{2}}{x}_{\mathrm{3}}{)}^{\frac{\mathrm{3}}{\mathrm{5}}}\\ y_{\mathrm{1}}=x_{\mathrm{2}}\\ y_{\mathrm{2}}=x_{\mathrm{3}}\end{array}\right.（13）

式中，v(t)为输入信号，俯仰角的缩放系数\xi=0.1，滚转角的缩放系数\xi=0.04。上述三阶微分器解决了非连续函数的求导问题，使得为达到期望轨迹而计算出的期望俯仰角和期望滚转角指令信号{\theta }_d和{\phi }_d可以不连续。在位置控制中引入切换函数，三阶微分器的最后一层在接收到带有干扰的信号时，通过连续的积分可以有效减少干扰信号的影响^［13］。

为了证明三阶微分器在有限时间内的稳定性，定义系统为

V(t)=\frac{\mathrm{1}}{\mathrm{2}}{x}_{\mathrm{2}}^{\mathrm{2}}(t)+\frac{\mathrm{1}}{\mathrm{14}}{\left[x_{\mathrm{1}}(t)-v(t)+\frac{\zeta x_{\mathrm{2}}(t)}{\mathrm{9}}\right]}^{\mathrm{14}/\mathrm{3}}（14）

其中的所有项都是平方项或更高次幂项，因此这个函数是正定的。V(t)对时间求导得到

\dot{V}(t)=x_{\mathrm{2}}(t){\dot{x}}_{\mathrm{2}}(t)+x_{\mathrm{3}}(t){\dot{x}}_{\mathrm{3}}(t)+

\frac{\mathrm{7}}{\mathrm{6}}{\left(x_{\mathrm{1}}(t)-v(t)+\frac{\zeta x_{\mathrm{2}}(t)}{\mathrm{9}}\right)}^{\mathrm{15}/\mathrm{3}}\left(x_{\mathrm{2}}(t)+\frac{\zeta x_{\mathrm{3}}(t)}{\mathrm{9}}\right)（15）

将式（13）中的{\dot{x}}_{\mathrm{1}}、{\dot{x}}_{\mathrm{2}}和{\dot{x}}_{\mathrm{3}}代入式（15）中可得到

\begin{array}{l}\dot{V}(t)=-\frac{\mathrm{2}}{\mathrm{3}}{x}_{\mathrm{3}}(t){\left(x_{\mathrm{1}}(t)-v(t)+\frac{\zeta x_{\mathrm{2}}(t)}{\mathrm{9}}\right)}^{\mathrm{7}/\mathrm{3}}-\\ \frac{\mathrm{4}}{\mathrm{35}}{\zeta }^{\mathrm{2}}{x}_{\mathrm{3}}^{\mathrm{2}}(t)\end{array}（16）

注意到\dot{V}(t)中的-\frac{\mathrm{4}}{\mathrm{35}}{\zeta }^{\mathrm{2}}{x}_{\mathrm{3}}^{\mathrm{2}}(t)项总是非正的，而-\frac{\mathrm{2}}{\mathrm{3}}{x}_{\mathrm{3}}(t){\left(x_{\mathrm{1}}(t)-v(t)+\frac{\zeta x_{\mathrm{2}}(t)}{\mathrm{9}}\right)}^{\mathrm{7}/\mathrm{3}}的符号取决于x_{\mathrm{3}}(t)和{\left(x_{\mathrm{1}}(t)-v(t)+\frac{\zeta x_{\mathrm{2}}(t)}{\mathrm{9}}\right)}^{\mathrm{7}/\mathrm{3}}的符号，又由于{\left(x_{\mathrm{1}}(t)-v(t)+\frac{\zeta x_{\mathrm{2}}(t)}{\mathrm{9}}\right)}^{\mathrm{7}/\mathrm{3}}是一个奇数次幂，它将与x_{\mathrm{3}}(t)有相同的符号，因此\dot{V}(t)是负定的，保证了V(t)会随时间下降。为了证明有限时间稳定性，需要证明V(t)在有限时间内达到零。由于\dot{V}(t)是负定的，V(t)是严格递减的，可以进一步分析\dot{V}(t)的下降速率：当x_{\mathrm{1}}(t)-v(t)+\frac{\zeta x_{\mathrm{2}}(t)}{\mathrm{9}}较大时，非线性项{\left(x_{\mathrm{1}}(t)-v(t)+\frac{\zeta x_{\mathrm{2}}(t)}{\mathrm{9}}\right)}^{\mathrm{7}/\mathrm{3}}的值迅速增大，导致\dot{V}(t)的值迅速减小，即V(t)的值迅速下降，这种迅速减小意味着系统在某个时间T内就能达到稳态，函数V(t)将减小到零，此时系统状态x_{\mathrm{1}}、x_{\mathrm{2}}、x_{\mathrm{3}}达到稳态。线性阻尼项-\frac{\mathrm{4}}{\mathrm{35}}{\zeta }^{\mathrm{2}}{x}_{\mathrm{3}}^{\mathrm{2}}(t)确保了V(t)下降不会无限快，从而避免了系统的不稳定^［14］。因此通过分析函数及其导数，证明了给定的三阶微分器在有限时间内收敛到零，从而证明了其有限时间稳定性。

4 仿真分析

4.1 仿真参数设定

为了验证所设计的前馈控制器的有效性，基于六旋翼模型和转动惯量变化模型进行仿真分析。两栖无人机初始位置和姿态设定为

x(\mathrm{0})=\mathrm{0}，y(\mathrm{0})=\mathrm{0}，z(\mathrm{0})=\mathrm{0}

\phi (\mathrm{0})=\mathrm{0}，\theta (\mathrm{0})=\mathrm{0}，\psi (\mathrm{0})=\mathrm{0}

两栖无人机的运动轨迹设定为螺旋曲线，随时间变化的方程为

x_d=\mathrm{3}\mathrm{c}\mathrm{o}\mathrm{s}t，y_d=\mathrm{3}\mathrm{s}\mathrm{i}\mathrm{n}t，z_d=\mathrm{2}+\mathrm{0.5}t

设定两栖无人机受到的干扰为

f_x=-\mathrm{2}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{0.3}\mathrm{\pi }t)，f_y=\mathrm{2.4}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{0.3}\mathrm{\pi }t)，f_z=\mathrm{2.4}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{0.3}\mathrm{\pi }t)，{\tau }_{f\phi }=-I_x\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{0.3}\mathrm{\pi }t)，{\tau }_{f\theta }=I_y\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{0.3}\mathrm{\pi }t)， {\tau }_{f\psi }=I_z\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{0.3}\mathrm{\pi }t)

两栖无人机结构参数如表4所示。PID控制参数如表5所示。

4.2 仿真结果分析

前馈控制器的仿真结果如图8~11所示。从图中可以看出，无人机的位置跟踪效果较好，滚转角和俯仰角跟踪时，先出现了大幅度的振荡，后又快速地稳定了下来；偏航角跟踪时出现了轻微振荡，从图11的振荡部分放大图可以看出，无人机偏航角的振荡幅度在0.01 rad左右，小于0.05 rad，处在较小范围。

为了验证前馈控制器的性能，本文将其与传统PID控制器进行了同步对比仿真。图12（a）显示，在x方向上，前馈控制器的响应速度略慢于PID控制器。图12（a）和图12（b）展示了x和y方向上的位置随时间变化的情况。在这两幅图中，前馈控制轨迹与期望轨迹几乎完全重合，显示出前馈控制在跟踪期望位置方面的高精度。而PID控制轨迹虽然在大部分时间内能够跟随期望轨迹，但在波峰和波谷处出现较大的偏差，这可能导致系统在实际应用中出现过度振荡或响应迟缓的问题。特别是在图12（c）中，当存在时变的转动惯量和干扰时，PID控制的无人机在z轴方向出现了明显的振荡，而前馈控制的无人机则保持了稳定。此外，图12（d）期望轨迹与前馈控制轨迹非常接近，几乎重合，这表明前馈控制在跟踪期望轨迹方面表现出色。相比之下，PID控制轨迹虽然也试图跟随期望轨迹，但在某些区域出现了明显的偏差，尤其是在轨迹的转弯处，这表明PID控制在处理复杂轨迹时不如前馈控制稳定和精确，进一步证明了前馈控制器的优越性。这些仿真结果表明，前馈控制器在处理复杂动态和不确定性方面具有明显的优势。

综合以上分析可以得出，前馈控制在期望轨迹跟踪方面明显优于传统的PID控制。前馈控制能够在转动惯量不断变化的情况下更精确地跟随复杂的动态轨迹，减少偏差和振荡，从而提高系统的稳定性和响应速度。

5 结束语

本文提出了一种针对陆空两栖无人机的机械腿在空中运动时转动惯量发生变化，从而影响无人机飞行稳定问题的处理方法，设计出覆盖转动惯量最大值、最小值和最大变化值的机械腿运动轨迹，然后将该轨迹对应的转动惯量转化为转动惯量与时间之间的关系，改变了传统的无人机转动惯量是定值的设定，创新性地将非线性时变的转动惯量引入到无人机控制模型中，并且提出了一种转动惯量不断变化的两栖无人机的前馈控制策略。仿真结果表明，本文设计的前馈控制策略相比于传统的PID控制，在转动惯量和干扰都随时间变化的情况下能够更准确地实现位置跟踪。

参考文献

[1] CANAL I， REIMBOLD M， CAMPOS M. Drone use to combat COVID⁃19： adaptive tuning proposal of the control system under variable load［J］. IEEE Latin America Transactions， 2021， 19（6）： 901-908.

[2] 何乃晨. 基于模糊自适应PID控制的四旋翼无人机飞控系统研究［D］. 天津：河北工业大学，2017.

[3] HE N C. Research on flight control system of four⁃rotor UAV based on fuzzy adaptive PID control［D］. Tianjin： Hebei University of Technology， 2017. （in Chinese）

[4] 张长箭，张帆，朱正吼，等. 四旋翼无人机建模及模糊PID控制［J］. 南昌航空大学学报（自然科学版），2024，38（2）：39-48.

[5] ZHANG C J， ZHANG F， ZHU Z H， et al. Modeling and fuzzy PID control of a quadcopter［J］. Journal of Nanchang Hangkong University（Natural Sciences）， 2024， 38（2）： 39-48. （in Chinese）

[6] GOJE S， BEDMUTHA P， SHAH P，et al. Cascaded PID control system for UAV with gain factor prediction using ML［J］. International Journal for Research in Applied Science and Engineering Technology， 2021， 9（11）： 665-674.

[7] 薛晶勇，杨波，朱永乐，等. 四旋翼不确定干扰前馈补偿与反步滑模姿态控制［J］. 计算机集成制造系统，2022，28（2）：345-354.

[8] XUE J Y， YANG B， ZHU Y L， et al. Feedforward compensation and backstepping sliding mode control for attitude of quadrotor under uncertain disturbance［J］. Computer Integrated Manufacturing Systems， 2022， 28（2）： 345-354.（in Chinese）

[9] YANG C S， YANG Z， HUANG X N， et al. Modeling and robust trajectory tracking control for a novel six⁃rotor unmanned aerial vehicle［J］. Mathematical Problems in Engineering， 2013， 2013： 673525.

[10] 徐雪松，吴仪政，倪赟. 变载荷四旋翼无人机的自适应模糊自抗扰控制［J］. 传感器与微系统，2022，41（7）：101-105， 109.

[11] XU X S， WU Y Z， NI Y. Adaptive fuzzy auto disturbance rejection control of variable load quadrotor UAV［J］. Transducer and Microsystem Technologies， 2022， 41（7）： 101-105， 109. （in Chinese）

[12] 贺有智，刘同其. 四旋翼飞行器时延积分反演容错控制［J］. 系统工程与电子技术，2015，37（10）：2341-2346.

[13] HE Y Z， LIU T Q. Time delay integral backstepping based fault tolerant control of quadrotor aircraft［J］. Systems Engineering and Electronics， 2015， 37（10）： 2341-2346. （in Chinese）

[14] 李泽宇. 陆空两栖变结构无人机控制［D］. 广州：广东工业大学，2022.

[15] LI Z Y. Control of land-air amphibious variable structure UAV［D］. Guangzhou： Guangdong University of Technology， 2022. （in Chinese）

[16] GIERNACKI W， STĘPIEŃ S， CHODNICKI M， et al. Hybrid quasi⁃optimal PID⁃SDRE quadrotor control［J］. Energies， 2022， 15（12）： 4312.

[17] 熊有伦. 机器人学： 建模、控制与视觉［M］. 武汉：华中科技大学出版社，2020.

[18] XIONG Y L. Robotics modeling， control and vision［M］. Wuhan： Huazhong University of Science and Technology Press， 2020. （in Chinese）

[19] 庄红超. 电驱动大负重比六足机器人结构设计及其移动特性研究［D］. 哈尔滨：哈尔滨工业大学，2014.

[20] ZHUANG H C. Electrically driven large⁃load⁃ratio six⁃legged robot structural design and its mobile characteristics research［D］. Harbin： Harbin Institute of Technology， 2014. （in Chinese）

[21] SUN C R， AGHA S A， MOHAMED Z， et al. Optimised sliding mode control of a hexacopter： simulation and experiments［J］. Electronics， 2022， 11（16）： 2519.

[22] RECALDE L F， GUEVARA B S， CARVAJAL C P， et al. System identification and nonlinear model predictive control with collision avoidance applied in hexacopters UAVs［J］. Sensors， 2022， 22（13）： 4712.