中国水利水电工程进入国际先进行列，其中南水北调工程是中国的战略性工程，自建成以来已成为280多个县市的主要水利来源。由于支撑南水北调工程的关键建筑物——引水隧洞具有洞线长、洞径大、水压高、围岩地质复杂的情况，且常年受水力冲蚀、生化材料侵蚀等地质作用的影响，开始出现裂缝、塌方和露筋等典型缺陷。人工检测缺陷不仅耗时耗力，且检测的精度和时效性都很低，因此亟须取得水下机器人巡检技术的新突破，以便快速、高效地完成水下复杂环境的检测。智能水下机器人是目前研究的热点，其中水下机械臂不仅可以安装在水下航行器上，还可以选择性地安装在需要的平台上，完成清理水面、电缆敷设和修理、打捞沉物及切断绳索等任务^[1-4]。水下机械臂因自身耦合动力学及水动力学环境影响，使水下机械臂的动力学模型具有非线性瞬变特性，因此其动力学建模^[5]和运动关节的高精度控制^[6-8]成为面临的主要挑战。

本文首先提出了一种融合Newton-Euler方程、Morison方程与非线性摩擦力的水下机械臂动力学模型建模方法，并在此基础上完成了动力学参数辨识。为提高机械臂在复杂瞬变水下环境的精准控制能力，设计了一种在补偿非线性动力学模型基础上利用径向基函数(radial basis function，RBF)神经网络补偿系统未建模与建模误差的自适应滑模控制方法。通过仿真，证明了该方法比传统的比例积分微分(proportional integral differential，PID)控制和一般的RBF网络自适应滑模控制具有更高的轨迹跟踪精度，即控制精度。

1 改进的机械臂水动力学模型

水下机械臂动力学模型需要考虑重力、惯性力、离心力、Coriolis力、摩擦力、水阻力、附加质量力和浮力等环境因素，导致模型复杂难以求解。本文在现有的水下机械臂动力学模型基础上，提出一种改进的易于计算的动力学模型。建模方式为：1) 引入考虑非线性黏滞与温度影响的摩擦力模型；2) 将水阻力、附加质量力以外力的形式引入机械臂动力学模型的计算过程，简化模型的求解过程。

1.1 机械臂动力学建模

本文研究的六轴水下机械臂由k个关节与k个连杆组成，假设机械臂各个关节都是刚性连接，根据文[9-10]可建立机械臂的刚性动力学模型。各个关节坐标系速度与加速度的传递关系如下：

$ \begin{array}{c} {{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen, }}n + 1}} = _n^{n + 1}{\mathit{\boldsymbol{T}}_{{\rm{ro}}}}{{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen, }}n}} + {{\dot \theta }_{n + 1}}^{n + 1}{{\mathit{\boldsymbol{\hat Z}}}_{n + 1}}, \\ {\mathit{\boldsymbol{\theta }}_{{\rm{gen}}, 1}} = {\left[ {\begin{array}{*{20}{l}} {0, }&{0, }&0 \end{array}} \right]^{\rm{T}}}{\rm{. }} \end{array} $

(1)

其中：$\mathit{\boldsymbol{\dot \theta }}$_{gen, n}∈R^3×1，为第n(n=1, 2，…, k)个连杆相对于广义坐标系的角速度，R^a×b为a×b矩阵，R^a×1为a×1向量；_nⁿ⁺¹T_ro∈R^3×3为机械臂第n个关节移动到第n+1个关节的旋转矩阵；θ_n、$\dot \theta $_n和$\ddot \theta $_n分别表示第n个关节的角度、角速度和角加速度；ⁿ$\mathit{\boldsymbol{\hat Z}}$_n∈R^3×1，为第n个连杆相对于广义坐标系的旋转轴方向向量。

角加速度的传递关系如下：

$ \begin{array}{c} {{\mathit{\boldsymbol{\ddot \theta }}}_{{\rm{gen}}, n + 1}} = _n^{n + 1}{T_{{\rm{ro}}}}{{\mathit{\boldsymbol{\ddot \theta }}}_{{\rm{gen}}, n}} + \left( {_n^{n + 1}{T_{{\rm{ro}}}}{\mathit{\boldsymbol{{\dot \theta }}}_{{\rm{gen}}, n}}} \right) \times \\ \left( {{{\dot \theta }_{n + 1}}{{\mathit{\boldsymbol{\hat Z}}}_{n + 1}}} \right) + {{\ddot \theta }_{n + 1}}{{\mathit{\boldsymbol{\hat Z}}}_{n + 1}}{\rm{, }}\\ {{\mathit{\boldsymbol{\ddot \theta }}}_{{\rm{gen}}, 1}} = {\left[ {\begin{array}{*{20}{c}} {0, 0, 0} \end{array}} \right]^{\rm{T}}}. \end{array} $

(2)

线速度的传递关系如下：

$ {\mathit{\boldsymbol{v}}_{n + 1}} = _n^{n + 1}{\mathit{\boldsymbol{T}}_{{\rm{ro}}}}\left( {{\mathit{\boldsymbol{v}}_n} + {{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen}}, n}} \times _n^{n + 1}{\mathit{\boldsymbol{T}}_{{\rm{po}}}}} \right), {\mathit{\boldsymbol{v}}_1} = {[0, 0, 0]^{\rm{T}}}. $

(3)

其中：$\mathit{\boldsymbol{ v}}$_n∈R^3×1，为第n个连杆相对于广义坐标系的线速度; _nⁿ⁺¹T_po∈R^3×1, 为机械臂第n个关节移动到第n+1个关节的平移矩阵。

线加速度的传递关系如下：

$ \begin{array}{c} {{\mathit{\boldsymbol{\dot v}}}_{n + 1}} = _n^{n + 1}{\mathit{\boldsymbol{T}}_{{\rm{ro}}}}\left( {{{\mathit{\boldsymbol{\ddot \theta }}}_{{\rm{gen}}, n}} \times _n^{n + 1}{{\bf{T}}_{{\rm{po}}}} + {{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen}}, n}} \times } \right.\\ \left. {\left( {{{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen}}, n}} \times _n^{n + 1}{\mathit{\boldsymbol{T}}_{{\rm{po}}}}} \right) + {{\mathit{\boldsymbol{\dot v}}}_n}} \right), \\ {{\mathit{\boldsymbol{\dot v}}}_1} = {[0, 0, g]^{\rm{T}}}{\rm{. }} \end{array} $

(4)

其中：${\mathit{\boldsymbol{\dot v}}}$_n∈R^3×1，为第n个连杆相对于广义坐标系的线加速度，g为重力加速度。

质心位置线速度和线加速度为：

$ {\mathit{\boldsymbol{v}}_{{\rm{cen}}, n}} = {\mathit{\boldsymbol{\dot \theta }}_{{\rm{gen}}, n}} \times {\mathit{\boldsymbol{P}}_{{\rm{cen}}, n}} + {\mathit{\boldsymbol{v}}_n}, \quad {\mathit{\boldsymbol{v}}_{{\rm{cen}}, 0}} = {[0, 0, 0]^{\rm{T}}}; $

(5)

$ \begin{array}{c} {{\mathit{\boldsymbol{\dot v}}}_{{\rm{cen}}, n}} = {{\mathit{\boldsymbol{\ddot \theta }}}_{{\rm{gen}}, n}} \times {\mathit{\boldsymbol{P}}_{{\rm{cen}}, n}} + {{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen}}, n}} \times \left( {{{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen}}, n}} \times {\mathit{\boldsymbol{P}}_{{\rm{cen, }}n}}} \right) + \\ {{\mathit{\boldsymbol{\dot v}}}_n}, \quad {{\mathit{\boldsymbol{\dot v}}}_{{\rm{cen}}, 0}} = {[0, 0, g]^{\rm{T}}}{\rm{.}} \end{array} $

(6)

其中：v_{cen, n}∈R^3×1, ${\mathit{\boldsymbol{\dot v}}}$_{cen, n}∈R^3×1, 分别为第n个连杆相对于质心坐标系的线速度和线加速度; P_{cen, n}∈ R^3×1, 为第n个连杆相对于广义坐标系下质心位置向量。

利用Newton-Euler公式分别计算位于质心位置处的惯性力和质心力矩为：

$ {\mathit{\boldsymbol{F}}_{{\rm{gen}}, n + 1}} = {m_{n + 1}}{\mathit{\boldsymbol{\dot v}}_{{\rm{cen}}, n + 1}}{\rm{, }} $

(7)

$ {\mathit{\boldsymbol{N}}_{n + 1}} = {\mathit{\boldsymbol{I}}_{{\rm{cen}}, n + 1}}{\mathit{\boldsymbol{\ddot \theta }}_{{\rm{gen}}, n + 1}} + {\mathit{\boldsymbol{\dot \theta }}_{{\rm{gen}}, n + 1}} \times \left( {{\mathit{\boldsymbol{I}}_{{\rm{cen}}, n + 1}}{{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen}}, n + 1}}} \right). $

(8)

其中: F_{cen, n}∈R^3×1, 为第n个关节相对于质心坐标系的惯性力; N_n∈R^3×1, 为第n个关节相对于质心坐标系的力矩; m_n为第n个连杆的质量; I_{cen, n}∈ R^3×3, 为第n个连杆相对于质心坐标系的惯性。

分别根据力平衡方程和力矩平衡方程可得：

$ {\mathit{\boldsymbol{F}}_n} = _{n + 1}^n{\mathit{\boldsymbol{T}}_{{\rm{ro}}}}{\mathit{\boldsymbol{F}}_{n + 1}} + {\mathit{\boldsymbol{F}}_{{\rm{cen}}, n}}, \quad {\mathit{\boldsymbol{F}}_k} = {\left[ {\begin{array}{*{20}{c}} {0, 0, 0} \end{array}} \right]^{\rm{T}}}; $

(9)

$ \begin{array}{c} {\mathit{\boldsymbol{\tau }}_{{\rm{NE}}, n}} = {\mathit{\boldsymbol{N}}_n} + _{n + 1}^n{\mathit{\boldsymbol{T}}_{{\rm{ro}}}}{\mathit{\boldsymbol{\tau }}_{{\rm{NE}}, n}} + {\mathit{\boldsymbol{P}}_{{\rm{cen}}, n}} \times {\mathit{\boldsymbol{F}}_{{\rm{cen}}, n}} + \\ _n^{n + 1}{\mathit{\boldsymbol{T}}_{{\rm{po}}}} \times \left( {_{n + 1}^n{\mathit{\boldsymbol{T}}_{{\rm{ro}}}}{\mathit{\boldsymbol{F}}_{n + 1}}} \right). \end{array} $

(10)

其中：F_n∈R^3×1，为机械臂各个连杆的受力；τ_{NE, n}∈R^3×1，为基于Newton-Euler公式计算的第n个关节的关节力矩。

由机械臂惯性力、Coriolis力、离心力和重力组成的关节力矩τ_MCG∈R^k×1，表示如下：

$ {\mathit{\boldsymbol{\tau }}_{{\rm{MCG}}}}(n) = \mathit{\boldsymbol{\tau }}_{{\rm{NE}}{,n}}^{\rm{T}}{\mathit{\boldsymbol{\hat Z}}_n}. $

(11)

此外，τ_MCG可表示为

$ {\mathit{\boldsymbol{\tau }}_{{\rm{MCG}}}} = {\mathit{\boldsymbol{M}}_{{\rm{PDI}}}}(\mathit{\boldsymbol{q}})\mathit{\boldsymbol{\ddot q}} + \mathit{\boldsymbol{C}}(\mathit{\boldsymbol{q}}, \mathit{\boldsymbol{\dot q}})\mathit{\boldsymbol{\dot q}} + \mathit{\boldsymbol{G}}(\mathit{\boldsymbol{q}}). $

(12)

其中：$\mathit{\boldsymbol{q、\dot q、\ddot q}}$∈R^k×1，分别表示机械臂k个关节的角度、角速度和角加速度向量，其由θ_n、$\dot \theta $_n和$\ddot \theta $_n组成，q =[θ₁, θ₂, …, θ_k]^T, $\mathit{\boldsymbol{\dot q}}$=[$\dot \theta $₁, $\dot \theta $₂, …, $\dot \theta $_k]^T, $\mathit{\boldsymbol{\ddot q}}$=[$\ddot \theta $₁, $\ddot \theta $₂, …, $\ddot \theta $_k]^T；M_PDI(q)∈R^k×k, 为正定的惯性矩阵; C(q, $\mathit{\boldsymbol{\dot q}}$)∈R^k×k, 为Coriolis-向心矩阵; G(q)∈R^k+1, 为重力。

机械臂各个关节处由于具有高减速比的谐波减速器，因此电机惯性力矩也应被考虑，表示如下：

$ {\mathit{\boldsymbol{\tau }}_{{\rm{iner}}}} = {\mathit{\boldsymbol{I}}_{{\rm{motor}}}}\mathit{\boldsymbol{\ddot q}}. $

(13)

其中：I_motor∈R^k×1, 为电机的转动惯量; τ_iner∈R^k×1, 为机械臂各个关节电机的惯性力矩。

1.2 非线性摩擦力建模

非零速区域的电机摩擦模型通常可描述为：

$ {\mathit{\boldsymbol{\tau }}_{{F_{{\rm{fri}}}}}} = \left( {{\mathit{\boldsymbol{\tau }}_{{F_{{\rm{coul}}}}}} + {\mathit{\boldsymbol{\tau }}_{{F_{{\rm{vis}}}}}}} \right){\rm{sign}}(\mathit{\boldsymbol{\dot q}}), $

(14)

$ {\mathit{\boldsymbol{\tau }}_{{F_{{\rm{vis}}}}}} = {\mathit{\boldsymbol{f}}_{{\rm{vis}}}}|\mathit{\boldsymbol{\dot q}}|. $

(15)

其中：τ_Ffri∈R^k×1, 为关节电机摩擦力矩；τ_Fcoul= f_coul, 为Coulomb摩擦力矩，f_coul∈R^k×1, 为Coulomb摩擦系数向量；τ_Fvis为黏性摩擦力矩；f_vis∈R^k×1, 为黏性摩擦系数向量。

然而文[11-12]中通过实验证明关节电机的黏性摩擦力通常与速度并不呈线性关系，通过引入指数参数向量δ_vis∈R^k×1，可更逼近真实的黏性摩擦力矩，即

$ {\mathit{\boldsymbol{\tau }}_{{F_{{\rm{vis}}}}}} = {\mathit{\boldsymbol{f}}_{{\rm{vis}}}}|\mathit{\boldsymbol{\dot q}}{|^{{\mathit{\boldsymbol{\delta }}_{{\rm{vis}}}}}}. $

(16)

此外，文[11]中也指出温度对于摩擦力具有明显的非线性影响，在此基础上，本文利用线性模型形式来表征温度K对摩擦模型中f_vis、δ_vis、f_coul的影响，即f_coul由z个温度影响下的Coulomb系数向量f_{coul, z}与温度K计算得到, f_vis由z个温度影响下的黏性摩擦系数向量f_{vis, z}与温度K计算得到, δ_vis由z个温度影响下的指数参数向量δ_{vis, z}与温度K计算得到，分别表示如下：

$ \begin{array}{c} {\mathit{\boldsymbol{f}}_{{\rm{coul }}}} = \sum\limits_{z = 0}^1 {{\mathit{\boldsymbol{f}}_{{\rm{coul }}, z}}} {K^z}, \\ {\mathit{\boldsymbol{f}}_{{\rm{vis}}}} = \sum\limits_{z = 0}^1 {{\mathit{\boldsymbol{f}}_{{\rm{vis}}, z}}} {K^z}, \\ {\mathit{\boldsymbol{\delta }}_{{\rm{vis}}}} = \sum\limits_{z = 0}^1 {{\mathit{\boldsymbol{\delta }}_{{\rm{vis}}, z}}} {K^z}{\rm{. }} \end{array} $

(17)

1.3 机械臂水动力学建模

在水下环境中，水下机械臂的各个关节都会受到自身运动和水流流动的影响。本文机械臂主要应用于长引水隧洞静水期检修任务，因此假设工况条件为完全静水状态。根据文[13-14]，机械臂关节n单位长度的物体与水之间的相互作用力F_{water, n}∈R^3×1, 可描述为

$ \frac{{{\rm{d}}{\mathit{\boldsymbol{F}}_{{\rm{water, }}n}}}}{{{\rm{d}}{\mathit{\boldsymbol{l}}_n}}} = \frac{{{\rm{d}}{\mathit{\boldsymbol{F}}_{{\rm{drag}}, n}}}}{{{\rm{d}}{\mathit{\boldsymbol{l}}_n}}} + \frac{{{\rm{d}}{\mathit{\boldsymbol{F}}_{{\rm{mass}}, n}}}}{{{\rm{d}}{\mathit{\boldsymbol{l}}_n}}} + \frac{{{\rm{d}}{\mathit{\boldsymbol{F}}_{{\rm{lift, }}n}}}}{{{\rm{d}}{\mathit{\boldsymbol{l}}_n}}} + \frac{{{\rm{d}}{\mathit{\boldsymbol{F}}_{{\rm{float}}, n}}}}{{{\rm{d}}{\mathit{\boldsymbol{l}}_n}}}. $

(18)

其中：F_{drag, n}∈R^3×1，为第n个关节的水阻力；F_{mass, n}∈R^3×1，为第n个关节的附加质量力；F_{lift, n}∈ R^3×1，为第n个关节的升力；F_{float, n}∈R^3×1，为第n个关节的浮力；l_n∈R^3×1，为第n个关节连杆的长度向量；且F_{drag, n}、F_{mass, n}、F_{lift, n}、F_{float, n}和l_n都相对于广义坐标系建立。由于水下机械臂没有翼型结构，因此不考虑F_{lift, n}对水下机械臂的影响。第n个连杆的F_{float, n}可表示为

$ {F_{{\rm{float, }}n}} = {m_n}g - {b_n}{\rho _{{\rm{water }}}}g = {m_n}g\left( {1 - \frac{{{\rho _{{\rm{water }}}}}}{{{\rho _{{\rm{arm }}, n}}}}} \right). $

(19)

其中：b_n为第n个连杆的体积，ρ_water为流体密度，ρ_{arm, n}为第n个机械臂的密度。

浮力对机械臂各个连杆的作用竖直向上且可认为浮力作用于各连杆质心，因此，通过修改式(6)的初始条件，即将g改为g·ρ_water/ρ_{arm, 1}, 则可将浮力代入机械臂动力学模型中，可得下式：

$ \begin{array}{c} {{\mathit{\boldsymbol{\dot v}}}_{{\rm{cen}}, n}} = {{\mathit{\boldsymbol{\ddot \theta }}}_{{\rm{gen}}, n}} \times {\mathit{\boldsymbol{P}}_{{\rm{cen}}, n}} + {{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen}}, n}} \times \left( {{{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen}}, n}} \times } \right.\\ \left. {{\mathit{\boldsymbol{P}}_{{\rm{cen}}, n}}} \right) + {{\mathit{\boldsymbol{\dot v}}}_n}, \\ {{\mathit{\boldsymbol{\dot v}}}_{{\rm{cen,}}0}} = {\left[ {0,0,g \cdot {\rho _{{\rm{water}}}}/{\rho _{{\rm{arm}},1}}} \right]^{\rm{T}}}{\rm{.}} \end{array} $

(20)

利用Morison公式可有效计算水阻力和附加质量力如下：

$ \begin{array}{c} {{\mathit{\boldsymbol{F}}}_{{\rm{drag}}, n}} = \\ \int_0^{{l_n}} {\frac{1}{2}} {\rho _{{\rm{water}}}}{C_{{\rm{drag}}}}{S_{{\rm{drag}}, n}} \cdot {\mathit{\boldsymbol{v}}_{{\rm{flow}}, n}}(\mathit{\boldsymbol{L}})\left\| {{\mathit{\boldsymbol{v}}_{{\rm{flow, }}n}}(\mathit{\boldsymbol{L}})} \right\| \cdot {\rm{d}}\mathit{\boldsymbol{L, }}\\ {\mathit{\boldsymbol{F}}_{{\rm{mass, }}\mathit{n}}} = \int_0^{{l_n}} {{\rho _{{\rm{water }}}}} {C_{{\rm{mass }}}}{S_{{\rm{mass, }}\mathit{n}}} \cdot {{\mathit{\boldsymbol{\dot v}}}_{{\rm{flow, }}\mathit{n}}}(\mathit{\boldsymbol{L}}) \cdot {\rm{d}}\mathit{\boldsymbol{L}}. \end{array} $

(21)

其中：位置L∈R^3×1；C_drag为水阻力系数，C_mass为附加质量力系数，一般情况下，C_drag与C_mass可通过实验获得，本文根据经验设计C_drag=1.0，C_mass=2.0；S_{drag, n}为第n个关节单位长度连杆垂直于流速方向的投影面积，本文机械臂各个连杆可近似认为是圆柱体，因此S_{drag, n}=D_n×1，D_n为第n个连杆的直径；S_{mass, n}为第n个连杆的轴截面积；v_{flow, n}(L)和${\mathit{\boldsymbol{\dot v}}}$_{flow, n}(L)分别为第n个连杆在L位置所受到的水流的速度和加速度。

文[13-14]中通过分别计算各个连杆对第n个连杆的水阻力力矩和附加力矩，获得水阻力和附加质量力的计算公式，但这种方法在计算k≥6的机械臂时，计算过程烦琐，本文根据式(1)、(2)、(5)和(6)求得v_{flow, n}(L)和${\mathit{\boldsymbol{\dot v}}}$_{flow, n}(L)分别为:

$ {\mathit{\boldsymbol{v}}_{{\rm{flow}}, n}}(\mathit{\boldsymbol{L}}) = \left( {{\mathit{\boldsymbol{v}}_{{\rm{cen, }}n}} - {{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen}}, n}} \times \mathit{\boldsymbol{L}}} \right), \\ {\mathit{\boldsymbol{\dot v}}_{{\rm{flow}}, n}}(\mathit{\boldsymbol{L}}) = \left( {{\mathit{\boldsymbol{\dot v}}_{{\rm{cen, }}n}} - {{\mathit{\boldsymbol{\ddot \theta }}}_{{\rm{gen}}, n}} \times \mathit{\boldsymbol{L}}} \right). $

(22)

可进一步得到下式：

$ \begin{array}{c} {\mathit{\boldsymbol{F}}_{{\rm{drag}}, n}} = \frac{1}{2}{\rho _{{\rm{water }}}}{C_{{\rm{drag}}}}{S_{{\rm{drag}}, n}}\left[ {\left( {\int_0^{{l_n}} {{{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen, }}n}}} \times } \right.} \right.\\ \left. {\left. {\mathit{\boldsymbol{L}}\left\| {{{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{gen}}, n}} \times \mathit{\boldsymbol{L}}} \right\|{\rm{d}}\mathit{\boldsymbol{L}}} \right) + {\mathit{\boldsymbol{l}}_n}{\mathit{\boldsymbol{v}}_{{\rm{cen}}, n}}} \right], \\ {\mathit{\boldsymbol{F}}_{{\rm{mass}}, n}} = {\rho _{{\rm{water }}}}{C_{{\rm{mass }}}}{S_{{\rm{mass}}, n}}\left[ {\left( {\int_0^{{l_n}} {{{\mathit{\boldsymbol{\ddot \theta }}}_{{\rm{gen, }}n}}} \times \mathit{\boldsymbol{L}} \cdot {\rm{d}}\mathit{\boldsymbol{L}}} \right) + } \right.\\ \left. {{\mathit{\boldsymbol{l}}_n}{{\mathit{\boldsymbol{\dot v}}}_{{\rm{cen}}, n}}} \right]{\rm{. }} \end{array} $

(23)

由此可计算第n个关节的水阻力F_{drag, n}和附加质量力F_{mass, n}。不同于文[13-14]中单独计算F_{drag, n}和F_{mass, n}并代入总的力矩方程的方法，本文设计将F_{drag, n}和F_{mass, n}直接代入动力学计算公式中，将式(9)修改如下：

$ \begin{array}{c} {\mathit{\boldsymbol{F}}_n} = _{n + 1}^n{\mathit{\boldsymbol{T}}_{{\rm{ro}}}}{\mathit{\boldsymbol{F}}_{n + 1}} + {\mathit{\boldsymbol{F}}_{{\rm{cen}}, n}} + {\mathit{\boldsymbol{F}}_{{\rm{water}}}}, \\ {\mathit{\boldsymbol{F}}_{{\rm{water}}}} = {\mathit{\boldsymbol{F}}_{{\rm{drag}}}} + {\mathit{\boldsymbol{F}}_{{\rm{mass}}}}. \end{array} $

(24)

综上可计算关节力矩为

$ \mathit{\boldsymbol{\tau }} = {\mathit{\boldsymbol{\tau }}_{{\rm{MCG}}}} + {\mathit{\boldsymbol{\tau }}_{{F_{{\rm{fri}}}}}} + {\mathit{\boldsymbol{\tau }}_{{\rm{iner}}}}{\rm{.}} $

(25)

其中：τ ∈R^k×1, 为关节力矩; 此时τ_MCG为机械臂由惯性力、Coriolis力、离心力、重力、水阻力、水浮力和水下附加质量力所组成的关节力矩。

2 机械臂动力学参数辨识 2.1 激励轨迹设计

在机械臂控制中，依托动力学模型进行前馈补偿可有效提高机械臂的控制效果，因此，动力学参数辨识对于水下机械臂的控制有重要影响^[15-18]。根据式(11)、(13)、(14)、(16)和(25)对机械臂动力学模型的描述，结合本文机械臂关节没有温度传感器和水动力环境复杂多变、难以准确快速辨识的情况，拟开展不考虑水动力学环境和温度影响的动力学参数辨识。针对式(16)中摩擦力的非线性黏滞，本文设计使用二次多项式拟合黏性摩擦力矩，表示如下：

$ {\mathit{\boldsymbol{\tau }}_{{F_{{\rm{vis}}}}}} = \left( {\mathit{\boldsymbol{f}}_{{\rm{vis}}}^1|\mathit{\boldsymbol{\dot q}}| + \mathit{\boldsymbol{f}}_{{\rm{vis}}}^2|\mathit{\boldsymbol{\dot q}}{|^2}} \right){\mathop{\rm sgn}} (\mathit{\boldsymbol{\dot q}}). $

(26)

其中：f_vis¹∈R^k×1，f_vis²∈R^k×1，分别为二次多项式非线性黏性摩擦力模型中的1次和2次性项系数向量。

根据文[19]，式(12)可改写为线性方程，此外根据式(13)和(26)也可以改写为线性方程，即：

$ \mathit{\boldsymbol{\hat \tau }} = \mathit{\boldsymbol{H}}(\mathit{\boldsymbol{q}}, \mathit{\boldsymbol{\dot q}}, \mathit{\boldsymbol{\ddot q}})\mathit{\boldsymbol{\delta }}. $

(27)

其中：$\mathit{\boldsymbol{\hat \tau }}$为观测力矩；H (q, $\mathit{\boldsymbol{\dot q}}$, $\mathit{\boldsymbol{\ddot q}}$)为观测矩阵；机械臂动力学参数向量δ ∈R^14k×1, 且δ由k个δ_n组成，δ_n∈R^14×1, 为第n个关节的动力学参数向量，具体表示如下：

$ \begin{array}{c} {\mathit{\boldsymbol{\delta }}_n} \equiv \left[ {{\mathit{\boldsymbol{I}}_{{\rm{cen}}, n}}(1, 1), {\mathit{\boldsymbol{I}}_{{\rm{cen}}, n}}(1, 2), {\mathit{\boldsymbol{I}}_{{\rm{cen}}, n}}(1, 3), } \right.\\ {\mathit{\boldsymbol{I}}_{{\rm{cen}}, n}}(2, 2), {\mathit{\boldsymbol{I}}_{{\rm{cen}}, n}}(2, 3), {\mathit{\boldsymbol{I}}_{{\rm{cen}}, n}}(3, 3), \\ {l_{{\rm{cen}}\mathit{X}{\rm{, }}\mathit{n}}}, {l_{{\rm{cen}}Y, \mathit{n}}}, {l_{{\rm{cen}}Z, n}}, {m_n}, \\ \left. {{\mathit{\boldsymbol{I}}_{{\rm{motor}}}}(n), {\mathit{\boldsymbol{f}}_{{\rm{coul}}}}(n), \mathit{\boldsymbol{f}}_{{\rm{vis}}}^1(n), \mathit{\boldsymbol{f}}_{{\rm{vis}}}^2(n)} \right]{\rm{.}} \end{array} $

(28)

其中：l_{cenX, n}、l_{cenY, n}和l_{cenZ, n}分别为第n个连杆的质心相对于广义坐标系3条坐标轴的距离。根据文[16]，δ最小参数集合可简化为10(k－1)+4+4个参数矩阵β。五次Fourier轨迹^[20]常被作为机械臂动力学参数辨识的激励轨迹：

$ \begin{array}{c} {\mathit{\boldsymbol{\theta }}_{{\rm{traj, }}, n}}(t) = \sum\limits_{j = 1}^{{N_{{\rm{freq }}}}} {\left[ {\frac{{{\mathit{\boldsymbol{A}}_{j, n}}}}{{{w_{{\rm{freq }}}}j}}\sin \left( {{w_{{\rm{freq }}}}jt} \right) - } \right.} \\ \left. {\frac{{{\mathit{\boldsymbol{B}}_{j, n}}}}{{{w_{{\rm{freq }}}}j}}\cos \left( {{w_{{\rm{freq }}}}jt} \right)} \right] + {\mathit{\boldsymbol{q}}_{{\rm{init }}}}(n){\rm{. }} \end{array} $

(29)

其中：t表示时间，t∈[0, t_freq]，t_freq为函数周期；θ_{traj, n}(t) 为第n个关节的Fourier轨迹，其基频w_freq=2π/t_freq；谐波个数N_freq=5；A_{j, n}和B_{j, n}分别表示第n个关节的j阶Fourier轨迹中余弦函数和正弦函数的振幅数值矩阵，j=[1, N_freq], 为振幅系数；q_init∈R^k×1, 为机械臂各个关节的初始角度向量。

2.2 轨迹参数设计及参数辨识

假设共有Γ_max组测试实验数据，则可得回归矩阵Φ和实际扭矩ϕ分别为：

$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{H}}_1}\left( {\mathit{\boldsymbol{q}}_1^\prime , \mathit{\boldsymbol{\dot q}}_1^\prime , \mathit{\boldsymbol{\ddot q}}_1^\prime } \right)}\\ {{\mathit{\boldsymbol{H}}_2}\left( {\mathit{\boldsymbol{q}}_2^\prime , \mathit{\boldsymbol{\dot q}}_2^\prime , \mathit{\boldsymbol{\ddot q}}_2^\prime } \right)}\\ \vdots \\ {{\mathit{\boldsymbol{H}}_\varphi }\left( {\mathit{\boldsymbol{q}}_\varphi ^\prime , \mathit{\boldsymbol{\dot q}}_\varphi ^\prime , \mathit{\boldsymbol{\ddot q}}_\varphi ^\prime } \right)}\\ \vdots \\ {{\mathit{\boldsymbol{H}}_{{\mathit{\Gamma }_{\max }}}}\left( {\mathit{\boldsymbol{q}}_{\mathit{\Gamma }_{\max }}^\prime , \mathit{\boldsymbol{\dot q}}_{{\mathit{\Gamma }_{\max }}}^\prime , \mathit{\boldsymbol{\ddot q}}_{{\mathit{\Gamma }_{\max }}}^\prime } \right)} \end{array}} \right], \quad \mathit{\boldsymbol{\phi }} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\tau }}_1^\prime }\\ {\mathit{\boldsymbol{\tau }}_2^\prime }\\ \vdots \\ {\mathit{\boldsymbol{\tau }}_\varphi ^\prime }\\ \vdots \\ {\mathit{\boldsymbol{\tau }}_{{\mathit{\Gamma }_{\max }}}^\prime } \end{array}} \right]. $

(30)

其中: H_φ(q′_φ, ${\mathit{\boldsymbol{\dot q}}_\varphi ^\prime }$, ${\mathit{\boldsymbol{\ddot q}}_\varphi ^\prime }$)和τ′_φ分别为机械臂各关节角度、角速度、角加速度和力矩的第φ组测试的采集值集合矩阵，其中φ∈[1, Γ_max]。进而可采用遗传算法通过优化Φ的条件数获得最优的A_{j, n}和B_{j, n}参数^[17]：

$ \begin{array}{*{20}{c}} {\min }&{{\rm{cond}}(\mathit{\boldsymbol{ \boldsymbol{\varPhi} }})}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}}&{{\theta _{\min , n}} \le {\mathit{\boldsymbol{\theta }}_{{\rm{traj}}, n}} \le {\theta _{\max , n}}, }\\ {}&{{{\dot \theta }_{\min , n}} \le {{\mathit{\boldsymbol{\dot \theta }}}_{{\rm{traj, }}n}} \le {{\dot \theta }_{\max , n}}, }\\ {}&{{{\ddot \theta }_{\min , n}} \le {{\mathit{\boldsymbol{\ddot \theta }}}_{{\rm{traj}}, n}} \le {{\ddot \theta }_{\max , n}}.} \end{array} $

(31)

其中：θ_{min, n}和θ_{max, n}分别为关节n的最小角度值和最大角度值，${\dot \theta }$_{min, n}和${\dot \theta }$_{max, n}分别为关节n的最小角速度值和最大角速度值，$\ddot \theta $_{min, n}和$\ddot \theta $_{max, n}分别为关节n的最小角加速度值和最大角加速度值。

最后通过最小二乘可求解得出动力学参数β的辨识参数：

$ \mathit{\boldsymbol{\hat \beta }} = {\left( {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}} \right)^{ - 1}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} \phi}} . $

(32)

3 基于RBF神经网络自适应滑模控制 3.1 基于RBF神经网络逼近的自适应滑模控制

式(25)中τ_iner的I_motor可与τ_MCG的M_PDI(q)合并，可得

$ \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{q}}) = {\mathit{\boldsymbol{M}}_{{\rm{PDI}}}}(\mathit{\boldsymbol{q}}) + {\rm{diag}}\left( {{\mathit{\boldsymbol{I}}_{{\rm{motor}}}}} \right). $

(33)

因此式(25)可改写为

$ \mathit{\boldsymbol{\tau }} = \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{q}})\mathit{\boldsymbol{\ddot q}} + \mathit{\boldsymbol{C}}(\mathit{\boldsymbol{q}}, \mathit{\boldsymbol{\dot q}})\mathit{\boldsymbol{\dot q}} + \mathit{\boldsymbol{G}}(\mathit{\boldsymbol{q}}) + {\mathit{\boldsymbol{\tau }}_{{F_{{\rm{fri}}}}}}. $

(34)

以$\mathit{\boldsymbol{x}} = \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{x}}_1} = \mathit{\boldsymbol{q}}}\\ {{\mathit{\boldsymbol{x}}_2} = \mathit{\boldsymbol{\dot q}}} \end{array}} \right]$为系统输出，其中x₁∈R^k×1，x₂∈R^k×1，x ∈R^2k×1，系统输入u = τ ∈R^k×1，构建控制系统为

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot x}}}_1} = {\mathit{\boldsymbol{x}}_2}, }\\ {{{\mathit{\boldsymbol{\dot x}}}_2} = {\mathit{\boldsymbol{M}}_1}\mathit{\boldsymbol{u}} + {\mathit{\boldsymbol{M}}^{ - 1}}\left( { - \mathit{\boldsymbol{C\dot q}} - \mathit{\boldsymbol{G}} - {\mathit{\boldsymbol{\tau }}_{{F_{{\rm{fri}}}}}}} \right) + \mathit{\boldsymbol{d}}(t).} \end{array}} \right. $

(35)

其中：M₁为由向量M_RBF构造成的角矩阵，M_{RBF, n}为由M逆的对角元组成的向量，即M_{RBF, n}= M^－1(n, n)。令M₁=diag[M_RBF，1, M_RBF，2, …, M_{RBF, k}]，M₁∈R^k×k；d(t)∈R^k×1, 为干扰向量且有界，即d(t)≤ D_max，D_max∈R^k×1, 为最大干扰误差向量。

令fun₁(x)=M₁，fun₂(x)= M^－1(－C${\mathit{\boldsymbol{\dot q}}}$－ G － τ_Ffri)，且fun₁(x)∈R^k×k，fun₂(x)∈R^k×1，控制系统可改写为

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot x}}}_1} = {\mathit{\boldsymbol{x}}_2}, }\\ {{{\mathit{\boldsymbol{\dot x}}}_2} = {\rm{fu}}{{\rm{n}}_1}(\mathit{\boldsymbol{x}})\mathit{\boldsymbol{u}} + {\rm{fu}}{{\rm{n}}_2}(\mathit{\boldsymbol{x}}) + \mathit{\boldsymbol{d}}(t).} \end{array}} \right. $

(36)

根据动力学辨识参数$\mathit{\boldsymbol{\hat \beta }}$可获得辨识结果$\mathit{\boldsymbol{\hat M }}$₁、$\mathit{\boldsymbol{\hat M }}$、$\mathit{\boldsymbol{\hat C }}$、$\mathit{\boldsymbol{\hat G }}$和$\mathit{\boldsymbol{\hat \tau }}$_Ffri，然而系统存在未建模量，如水动力学和摩擦力温度等，因此控制系统可表示为

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{fu}}{{\rm{n}}_2}(\mathit{\boldsymbol{x}}) = [ - (\mathit{\boldsymbol{\hat C}} + \mathit{\boldsymbol{\tilde C}})\mathit{\boldsymbol{\dot q}} - (\mathit{\boldsymbol{\hat G}} + \mathit{\boldsymbol{\tilde G}}) - }\\ {\left. {\quad \;\;\;\;\left( {{{\mathit{\boldsymbol{\hat \tau }}}_{{F_{{\rm{fri}}}}}} + {{\mathit{\boldsymbol{\tilde \tau }}}_{{F_{{\rm{fri}}}}}}} \right)} \right]{{(\mathit{\boldsymbol{\hat M}} + \mathit{\boldsymbol{\tilde M}})}^{ - 1}}, }\\ {{\rm{fu}}{{\rm{n}}_1}(\mathit{\boldsymbol{x}}) = \left( {{{\mathit{\boldsymbol{\hat M}}}_1} + {{\mathit{\boldsymbol{\tilde M}}}_1}} \right).} \end{array}} \right. $

(37)

其中${\mathit{\boldsymbol{\tilde C}}}$、${\mathit{\boldsymbol{\tilde G}}}$、${\mathit{\boldsymbol{\tilde \tau }}}$_Ffri、${\mathit{\boldsymbol{\tilde M}}}$和${\mathit{\boldsymbol{\tilde M}}}$₁为误差量。

神经网络可有效拟合系统误差项^[21-22], 本文采用RBF网络分别逼近fun₁(x)与fun₂(x)，RBF网络的输入输出算法可表示如下：

$ \begin{array}{*{20}{c}} {{h_j} = \exp \left( {\frac{{{{\left\| {\mathit{\boldsymbol{x - }}{\mathit{\boldsymbol{c}}_{{\rm{neu}}, j}}} \right\|}^2}}}{{2b_{{\rm{RBF}}}^2}}} \right), }\\ {\mathit{\boldsymbol{h}}(\mathit{\boldsymbol{x}}) = {{\left[ {{h_1}, {h_2}, \cdots , {h_j}, \cdots , {h_J}} \right]}^{\rm{T}}}, }\\ {{\rm{fu}}{{\rm{n}}_2}(\mathit{\boldsymbol{x}}) = {\mathit{\boldsymbol{Q}}^ * }\mathit{\boldsymbol{h}}(x) + {\mathit{\boldsymbol{\varepsilon }} _2}, }\\ {{\rm{fun}}_1^\prime (\mathit{\boldsymbol{x}}) = {\mathit{\boldsymbol{V}}^ * }h(\mathit{\boldsymbol{x}}) + {\mathit{\boldsymbol{\varepsilon }}_1}, }\\ {{\rm{fu}}{{\rm{n}}_1}(\mathit{\boldsymbol{x}}) = {\rm{diag}}\left[ {{\rm{fun}}_1^\prime (\mathit{\boldsymbol{x}})} \right].} \end{array} $

(38)

其中：h_j为第j个神经元的Gauss基函数; c_{neu, j}∈R^2k×1，为第j个神经元中心向量; b_RBF为神经元Gauss核宽度; J为神经元输出个数; x为网络输入；h (x)∈R^J×1, 为Gauss基函数的输出向量；Q^*∈R^k×J, V^*∈R^k×J, 分别为逼近fun′₁(x)(fun′₁(x)∈R^k×1)和fun₂(x)的理想网络权值矩阵，fun₁(x)由fun′₁(x)重组而成；ε₁∈R^k×1和ε₂∈R^k×1均为网络逼近误差向量，其最大值分别为ε_M1(ε_M1∈R^k×1)和ε_M2(ε_M2∈R^k×1)，即|ε₁|≤ ε_M1，| ε₂|≤ ε_M2。则RBF的输出为

$ {\widehat {{\rm{fun}}}_1}(\mathit{\boldsymbol{x}}) = \mathit{\boldsymbol{\hat Vh}}(\mathit{\boldsymbol{x}})\quad {\widehat {{\rm{fun}}}_2}(\mathit{\boldsymbol{x}}) = \mathit{\boldsymbol{\hat Qh}}(\mathit{\boldsymbol{x}}){\rm{.}} $

(39)

其中: $\widehat {{\rm{fun}}}$₁(x)∈R^k×k和$\widehat {{\rm{fun}}}$₂(x)∈R^k×1均为RBF网络每次的迭代输出矩阵；$\mathit{\boldsymbol{\hat Q}}$∈R^k×J和$\mathit{\boldsymbol{\hat V}}$∈R^k×J均分别为RBF网络的迭代网络权值矩阵。设计误差函数为e = q_de－ q_re；q_de∈R^k×1，为机械臂各个关节角度期望轨迹向量；q_re∈R^k×1，为机械臂各个关节的实际轨迹向量；e ∈R^k×1，为误差向量。设计滑模函数如下：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{s}} = \mathit{\boldsymbol{\dot e}} + \mathit{\boldsymbol{ce}}, }\\ {\mathit{\boldsymbol{\dot s}} = \mathit{\boldsymbol{\ddot e}} + \mathit{\boldsymbol{c\dot e}}.} \end{array} $

(40)

其中: c ∈R^k×1，为滑模函数参数向量; ${\mathit{\boldsymbol{\dot e}}}$和${\mathit{\boldsymbol{\ddot e}}}$分别为误差的导数和误差的二次导函数；${\mathit{\boldsymbol{\dot s}}}$为滑膜函数的导函数。基于此设计控制率，即系统输入u为

$ \mathit{\boldsymbol{u}} = \frac{1}{{{{\widehat {{\rm{fun}}}}_1}(\mathit{\boldsymbol{x}})}}\left[ { - {{\widehat {{\rm{fun}}}}_2}(\mathit{\boldsymbol{x}}) + {{\mathit{\boldsymbol{\ddot q}}}_{{\rm{de}}}} + \mathit{\boldsymbol{c\dot e}} + \mathit{\boldsymbol{\eta }}{\mathop{\rm sgn}} (\mathit{\boldsymbol{s}})} \right]. $

(41)

其中：${\mathit{\boldsymbol{\ddot q}}}$_de∈R^k×1, 为机械臂各个关节的期望角加速度向量；η ∈R^k×1, 为干扰上界向量，η ≥ D_max。

将(41)代入(40)可得

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot s}} = \left( {{{\widehat {{\rm{fun}}}}_2}(\mathit{\boldsymbol{x}}) - {{\widehat {{\rm{fun}}}}_2}(\mathit{\boldsymbol{x}})} \right) - \mathit{\boldsymbol{\eta }}{\mathop{\rm sgn}} (\mathit{\boldsymbol{s}}) + }\\ {\left. {\widehat {\rm{fun}}{_1}(\mathit{\boldsymbol{x}}) - {\rm{fu}}{{\rm{n}}_1}(\mathit{\boldsymbol{x}})} \right)\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{d}}(t){{\widetilde {{\rm{fun}}}}_2}(\mathit{\boldsymbol{x}}) - }\\ {\mathit{\boldsymbol{\eta }}{\mathop{\rm sgn}} (\mathit{\boldsymbol{s}}) + {{\widetilde {{\rm{fun}}}}_1}(\mathit{\boldsymbol{x}})\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{d}}(t) = \mathit{\boldsymbol{\tilde Qh}}(\mathit{\boldsymbol{x}}) - {\mathit{\boldsymbol{\varepsilon }}_2} - }\\ {\mathit{\boldsymbol{\eta }}{\mathop{\rm sgn}} (\mathit{\boldsymbol{s}}) + \left[ {\mathit{\boldsymbol{\tilde Vh}}(\mathit{\boldsymbol{x}}) - {\mathit{\boldsymbol{\varepsilon }}_1}} \right]\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{d}}(t).} \end{array} $

(42)

其中${\widetilde {{\rm{fun}}}}$₁(x)、${\widetilde {{\rm{fun}}}}$₂(x)、${\mathit{\boldsymbol{\tilde Q}}}$和${\mathit{\boldsymbol{\tilde V}}}$为误差量。

3.2 控制器稳定性证明

设计Lyapunnov函数：

$ \mathit{\boldsymbol{\mu }} = \frac{1}{2}{\mathit{\boldsymbol{s}}^2} + \frac{1}{{2{\mathit{\boldsymbol{\gamma}} _1}}}\mathit{\boldsymbol{\tilde Q}}{\mathit{\boldsymbol{\tilde Q}}^{\rm{T}}} + \frac{1}{{2{\mathit{\boldsymbol{\gamma}} _2}}}\mathit{\boldsymbol{\tilde V}}{\mathit{\boldsymbol{\tilde V}}^{\rm{T}}}. $

(43)

其中：γ₁, γ₂＞0，γ₁, γ₂∈R^k×1，为自适应率参数向量。

对式(43)求导可得

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot \mu }} = \mathit{\boldsymbol{s\dot s}} + \frac{1}{{{\gamma _1}}}\mathit{\boldsymbol{\tilde Q}}{{\mathit{\boldsymbol{\dot {\tilde Q}}}}^{\rm{T}}} + \frac{1}{{{\gamma _2}}}\mathit{\boldsymbol{\tilde V}}{{\mathit{\boldsymbol{\dot {\tilde V}}}}^{\rm{T}}} = }\\ {\mathit{\boldsymbol{\tilde Q}}\left( {\mathit{\boldsymbol{sh}}(\mathit{\boldsymbol{x}}) - \frac{1}{{{\mathit{\boldsymbol{\gamma }}_1}}}{{\mathit{\boldsymbol{\dot {\hat Q}}}}^{\rm{T}}}} \right) + \mathit{\boldsymbol{\tilde V}}\left( {\mathit{\boldsymbol{sh}}(\mathit{\boldsymbol{x}})\mathit{\boldsymbol{u}} - \frac{1}{{{\mathit{\boldsymbol{\gamma }}_2}}}{{\mathit{\boldsymbol{\dot {\hat V}}}}^{\rm{T}}}} \right) + }\\ {\mathit{\boldsymbol{s}}\left( { - {\mathit{\boldsymbol{\varepsilon }}_2} - \mathit{\boldsymbol{\eta }}{\mathop{\rm sgn}} (\mathit{\boldsymbol{s}}) - {\mathit{\boldsymbol{\varepsilon }}_1}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{d}}(t)} \right).} \end{array} $

(44)

设计自适应控制率为

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\dot {\hat Q}}} = {\mathit{\boldsymbol{\gamma }}_1}\mathit{\boldsymbol{sh}}(\mathit{\boldsymbol{x}}), }\\ {\mathit{\boldsymbol{\dot {\hat V}}} = {\mathit{\boldsymbol{\gamma }}_2}\mathit{\boldsymbol{sh}}(\mathit{\boldsymbol{x}})\mathit{\boldsymbol{u}}.} \end{array}} \right. $

(45)

令η > ε_M2+ ε_M1u + d_max+ η₀, η₀∈R^k×1, η₀>0，则有$\mathit{\boldsymbol{\dot \mu }}$≤- η₀|s|≤0，且当$\mathit{\boldsymbol{\dot \mu }}$≡0时, s ≡0，根据LaSalle不变集原理，t→∞时，s →0。综上可得Lyapunnov函数是(43)正定的，其导数是负定的，从而证明了控制系统的稳定性。

4 机械臂动力学建模及参数辨识仿真 4.1 机械臂动力学建模

机械臂的动力学模型参数通过数模或根据经验设计的，参数如表 1所示，并假设该参数为真实的动力学模型参数。此外设定水密度ρ_water=1 000 kg/m³，重力加速度g =9.798 8 m/s²。

根据表 1中摩擦力的参数和式(14)、(16)及(17)建立机械臂各个关节的非线性摩擦模型。

关节2的非线性摩擦模型如图 1所示。

4.2 机械臂动力学参数辨识仿真

通过遗传算法优化，即式(31)，可得如图 2和3所示的机械臂最优激励轨迹图。

各机械臂的最优激励参数如表 2所示。

设计仿真系统的采样频率为1 000 Hz，完成数据采样，基于式(30)和(32)完成机械臂仿真环境下的动力学参数辨识，参数辨识结果如表 3所示。

将表 3中的动力学辨识参数代入式(27)中，可得图 4所示的辨识结果，由图 4可知辨识力矩对真实力矩具有较好的辨识效果。

将表 3中的动力学辨识参数代入式(27)中，并设计不同幅度和不同起始位置的Fourier验证轨迹，验证结果如图 5所示，其中6个关节的辨识误差分别为：5.66%、7.74%、9.97%、1.79%、5.16%和1.29%。

4.3 水动力学环境对机械臂的影响

基于表 1设定的水动力学参数，依据水动力学公式，即式(18)、(21)、(23)、(24)和(25)建立六自由度水下机械臂的水动力学公式，并完成仿真建模，如图 6所示。

由图 6可知，水动力学对于机械臂动力学的稳定性具有重要的影响，其中各个关节的最大水动力对力矩的影响分别为：23.21、101.86、35.49、13.84、1.28和0 N·m。

5 自适应滑模轨迹跟踪控制仿真

通过第4章的动力学建模，可获得一个较为详尽的水下机械臂动力学仿真环境。针对各个关节设立幅值分别为：5、20、15、10、8和5 N·m的Gauss随机误差，以表示噪声干扰，使用第4章中的验证轨迹进行验证，系统温度随时间变化而升高，设定仿真开始时系统温度为15 ℃，仿真结束时系统温度为35 ℃，其间呈线性增长。

在此基础上，本章将开展基于动力学模型的机械臂仿真控制。首先设计了实验1(P1)，使用双环PID控制器进行控制仿真；其次，为提高控制精度还开展了实验2(P2)，基于RBF神经网络补偿的滑模控制仿真；最后，为进一步减低滑模控制的振动，在控制系统中加入了动力学模型补偿完成实验3(P3)。

P1使用了双环PID控制，分别对速度环和位置环进行PID控制，6个关节轨迹跟踪误差如图 7所示，可知6个关节轨迹的跟踪误差最大可达0.022 rad。

P2使用RBF神经网络对系统建模误差和未建模项进行拟合补偿，P3则在P2的基础上添加了动力学模型补偿，其跟踪误差如图 8所示。由图 8可知：P2与P3的轨迹跟踪效果明显优于P1的效果，P3在对动力学模型进行补偿之后，跟踪效果也优于P2。

6 结论

针对水下机械臂精准建模及控制技术，本文提出了一种融合Newton-Euler方程、Morison方程与非线性摩擦力的水下机械臂动力学模型建模方法。通过引入多项式函数建立非线性黏滞摩擦力模型，并在此基础上设计了激励轨迹优化方法和动力学参数辨识方法。通过实验证明了动力学参数辨识有效。

为提高机械臂在复杂多变水环境下的精准控制能力，本文设计了一种在补偿非线性动力学模型的基础上利用RBF神经网络补偿系统未建模与建模误差的自适应滑模控制方法。通过仿真，本文证明了该方法具有比传统PID控制和一般的RBF神经网络自适应滑模控制精度更突出的优势。

下一步将进一步开展关于非线性摩擦力自适补偿方法的研究，并在实物平台上开展实验测试。