变拓扑多臂航天器高效统一动力学建模

王宏旭，柳子然，岳程斐，曹喜滨

(1. 哈尔滨工业大学航天学院，哈尔滨 150001； 2. 哈尔滨工业大学(深圳)空间科学与应用技术研究院，深圳 518055)

0 引 言

随着空间技术的不断发展，动态非合作目标抓捕、大型空间设施建造等航天任务对服务航天器提出了更高的要求。与传统单臂空间服务航天器相比，可变拓扑多臂航天器具备多臂协同、灵活重构的能力，是满足未来多样化在轨操作任务和需求的重要手段。变拓扑多臂航天器具有单臂作业、多臂协同、多机协同等多种任务模式，会产生开链、闭环、环链混合等多种拓扑结构，引起系统动力学模型发生相应变化。传统针对单一构型的航天器动力学建模方法效率低下，需要人工介入，极大地限制了在轨服务航天器的灵活性和工作效率。因此，快速、准确地建立其动力学模型，是变拓扑多臂航天器研究的基础与关键所在。

多臂航天器系统在不同任务形态下可能存在不同的拓扑结构，这类变拓扑机械系统也被称为变胞机构。变胞机构的概念由Dai等[1]于1996年首次提出。针对变胞机构动力学的研究主要分为拓扑构型的数学描述方法和相应的动力学建模两部分。对拓扑结构及其变化进行完整准确的描述，是深入研究变胞机构的前提和关键所在。目前常用的拓扑结构描述方法有拓扑图、Huston低序体阵列、关联矩阵和邻接矩阵等。白国超等[2]提出了用铰链邻接矩阵对机构构态变换进行描述；Dai等[3]用邻接矩阵描述了变胞机构的拓扑结构变化，给出了不同构态之间邻接矩阵的转化关系；吴艳荣等[4]将构态切换分为构件数减少、构件数增加、构件数不变三种情况，详细叙述了不同情况下各构态之间邻接矩阵的转化关系；王汝贵等[5]给出了一种既可衍生邻接矩阵又可描述变胞机构运动副空间位置的构态变换矩阵；Tian等[6]对变胞机构的各阶段进行分析，得到变胞源机构，进而得到相应的子结构。周杨等[7]基于各类变化矩阵，将拓扑结构变化过程转化为一系列的矩阵运算，并给出了统一的数学模型。刘云平等[8]基于通路矩阵和低序体矩阵对多体系统拓扑结构进行描述，建立移位算子与拓扑描述矩阵间的联系，但未建立相应的动力学模型。当前针对变胞机构拓扑构型描述的研究已经比较全面，但描述过程缺乏自主性，且拓扑构型描述与动力学建模相互独立，缺乏统一性。

在早期多臂航天器的动力学建模研究中，大多基于牛顿、欧拉、拉格朗日等经典动力学方法，以及基于达朗贝尔原理的Kane方法。但上述方法的计算效率至少为系统自由度数目的平方数量级O(n2)，在处理以多臂航天器系统为代表的大型、复杂、多自由度、多体系统时十分不利。Rodriguez[9]将卡尔曼滤波预测理论的状态方程与多体系统动力学系统进行比拟，发现了二者的内在联系，发展了O(n)阶多体系统动力学的空间算子代数法(Spatial operator algebra, SOA)，实现了多体动力学的高效率简单建模。魏承等[10]应用SOA方法针对漂浮基空间双臂机器人进行了动力学建模及仿真分析，得到兼顾高精度与高计算效率的动力学模型。胡景晨等[11]将空间算子代数与绝对节点坐标法相结合，得到了一个复杂度为O(n)，可处理非线性大变形系统的高效、高精度多柔体动力学方法。孙占庚等[12]用Kane方法建立了柔性机械臂的动力学方程，用假设模态法对方程进行离散求解，研究了不同模态截取情况对计算结果的影响。对于含有闭环的多体系统，黄文虎等[13]发展了约束的多体系统动力学Newton-Euler正交法，在继承原矢量力学基础上给出了结构完善的动力学方程解析表达式，最后在机器人系统中得到应用。陈礼等[14]通过引入约束正交补轴的概念，给出了几种典型铰/组合铰的适用于计算机编程和数值求解的约束方程生成方法，但需要给定的参数较多，输入较繁琐。Hu等[15]针对对于树系统和闭环系统，采用空间算子理论建立了系动力学模型，但本质上还是对多种拓扑构型分别建模；史加贝等[16]基于共旋坐标法建立了大变形太阳电池阵展开的多体动力学模型，研究了电池阵展开过程中的动力学响应及碰撞规律；邱雪松等[17]建立了两级往复太阳能帆板的柔性多体动力学模型，分析铰链间隙与驱动力之间的规律；游斌弟等[18]利用Lagrange和Newton方法建立了卫星太阳阵系统的多刚体动力学模型，研究了铰链副接触碰撞对太阳阵展开及卫星姿态的影响；黄泽兵等[19]基于Jourdian速度变分原理建立了空间桁架展开的动力学模型，分析了弹性变形对系统动力学特性的影响。综上，对于多臂航天器的动力学建模也已比较完善，计算效率至少为O(n)；但大部分研究仅针对单一拓扑构型，不具备变拓扑能力。一部分研究建立了变拓扑系统动力学模型，但主要应用于卫星帆板展开、天线展开等简单变拓扑，无法直接迁移至多臂航天器变拓扑动力学建模中。

针对上述问题，本文开展了变拓扑多臂航天器高效统一动力学建模研究。首先，提出了多臂航天器各连杆的构型参数，用以描述连杆间的连接关系，并提出多臂航天器拓扑描述矩阵自生成算法。其次，基于空间算子代数理论，将通路矩阵P与空间转移算子φ相结合，提出变拓扑多臂航天器统一建模方法。最后，将数值仿真结果与多体动力学仿真软件的仿真结果进行比对，验证所提方法的有效性。

1 空间算子代数建模基础理论

首先，以图1所示的n体单链系统为例，说明SOA建模理论。定义连杆标号为k,其中基座为n+1，末端连杆为1。

图1 单链系统示意图Fig.1 Diagram of a single chain system

每个连杆的速度、加速度、力均采用旋量表示如下：

(1)

式中：下角标k表示第k个关节；ωk,vk分别表示关节k的角速度和线速度矢量；Nk,Fk分别表示关节力矩和力矢量。

相邻连杆间速度、加速度递推关系为[20]

(2)

式中：φ(k+1,k)为第k根连杆到第k+1根连杆的空间转移算子；H(k)为第k个关节转轴的方向矢量或线位移单位矢量；θ(k)为绕关节k旋转的角度或线位移；a(k)为机械臂的科氏加速度。

相邻连杆间力、力矩递推关系为[20]

(3)

式中：M(k)为第k根连杆的惯性矩阵，仅与连杆自身特性有关；b(k)为机械臂的离心力。

定义机械臂的系统速度算子为V=[V(1),…,V(n-1),V(n)]T，并以同样的形式表示机械臂的加速度算子α、科氏力算子a、离心力算子b、力算子f、力矩算子T。则式(2)和式(3)可表示为

(4)

推导得到机械臂系统逆动力学方程为

(5)

式中：MG表示机械臂广义惯性矩阵；C表示机械臂的非线性力矩阵；MG表达式为

MG=HφMφTHT

(6)

以上式中φ为空间转移算子，与机械臂拓扑构型有关，对于单链结构，其表达式为

(7)

对于树状结构，可将其拆分为多个单链，每个单链为一个分支。在两个相连的分支中，更接近根部的分支为另一个分支的父分支。对于每个分支，其空间转移算子满足上述表示形式。对于整个树系统，其空间转移算子的分块元素满足如下关系：

(8)

由文献[19]可知，φj,k中的元素定义为：

φj,k(m,l)=φ(mj,lk)=φ(mj,mj-1)

φ(mj-1,mj-2)…φ(lk+1,lk)

(9)

式中：k,j分别表示第k、第j个分支；mj表示在j分支上的编号为m的连杆；lk表示在k分支上的编号为l的连杆。

在传统动力学求解过程中，对广义惯性阵MG的求逆是计算效率为O(n3)的过程，将消耗大量计算时间。在SOA理论中，可将广义惯量阵进行因式分解，简化求逆过程。具体地，采用卡尔曼滤波平滑方法将广义惯性阵表示为如下形式[21]：

MG=(E+HΦK)D(E+HΦK)T

(10)

其逆矩阵为[21]

(11)

式中：K为过渡矩阵;D为对角阵，其求逆过程是算法复杂度为O(n)阶的过程，相较传统动力学算法，大大提高了计算效率。

2 拓扑构型及其描述方法

2.1 变拓扑多体系统动力学建模的问题与挑战

传统动力学建模方法大多针对于固定拓扑构型的系统，应用于变拓扑航天器动力学建模时会遇到诸多问题与挑战，具体表现为：传统动力学建模过程中，拓扑构型描述与动力学建模相互独立，没有联系在一起，建模过程繁琐；传统动力学建模中部组件编号严格与拓扑构型相关，拓扑构型改变后需重新编号，建模效率低；拓扑构型及动力学描述过程仍需要大量人工介入，建模自主性差、效率低。本章将介绍拓扑构型及其描述方法。

2.2 变拓扑多臂航天器的拓扑结构

一般的多臂航天器如图2所示，可用树形拓扑描述。树的连杆即为航天器的臂杆；连接连杆的铰即为航天器的运动关节；树的每条分支即为航天器的每条臂；树的根连杆即为航天器的中心体；根连杆与惯性空间通过六自由度虚铰连接。图中S表示连杆的首端，E表示连杆的末端，下标表示所在连杆编号。对于多臂航天器系统，所有连杆都直接或间接与根连杆相连。对于同一铰连接的两个连杆，接近根连杆一侧的称为内接体，远离根连杆的称为外接体[13]。

图2 多臂航天器的基本拓扑结构Fig.2 Topological structure of multi-arm spacecraft

当多臂航天器需要执行大范围移动或搬运的在轨任务时，基本构型下的工作空间大小无法满足任务需求，可将某个连杆甚至某个完整分支连接至另一个分支的末端，从而大幅拓展航天器的工作空间。例如，将连杆1的S1端从分支1断开，连接至分支2连杆m的Em端，形成分支2′，如图3(a)所示；如果需要更大的工作空间，可将分支1的连杆k的Sk端与根连杆断开，E1端连接至分支2的Em端，形成分支2′，如图3(b)所示。当多臂航天器执行多臂夹持任务时，某些臂的末端与被夹持物体连接，形成环状拓扑结构，如图3(c)所示。当多臂航天器执行大规模灵巧操作任务时，需要兼顾操作高灵活度和大工作空间的要求，可将多个多臂航天器按照一定规律进行连接，组成大规模的空间树网系统，形成复杂的环链混合结构，如图3(d)所示。

图3 典型多臂航天器变拓扑构型Fig.3 Typical topological configurations of multi-arm spacecraft

2.3 构型参数定义

为实现空间多臂航天器拓扑结构在轨自主描述，避免人工介入，提升在轨服务效率。本文提出空间多臂航天器构型参数的概念。

对于如图1所示的多臂航天器，其基本组成单元为臂杆，臂杆间通过关节连接。将臂杆看作树结构中的结点，连接臂杆的关节看作树结构中连接结点的连线，则可定义航天器臂杆的构型参数。具体地，定义构型参数为Cf，其结构如图4所示。臂杆序号(Link ID)为臂杆在系统内的编号，该编号在同一航天器系统内唯一。由于航天器中心体通过六自由度虚铰与惯性空间链接，臂杆编号最大值为s+6，其中s表示系统中包含的臂杆总数，数值最大的6个编号表示6个单自由度虚拟体，分别表示虚铰三个方向平动和三个方向转动。臂杆的接口数(Degree)为该臂杆允许连接的最大杆数，设计完成后根据机电接口确定。一般来说，多臂空间航天器臂杆的接口数为2。航天器中心体也可看作连杆，其接口数为多臂航天器最大臂数。接口指针(Interface)为与对应接口连接的臂杆编号，接口的编号在出厂前确定，指针的数值通过接口处实际采集接口的识别码生成，指针在构型参数中储存位置与接口编号一一对应。特别的，机械臂末端的接口指针为0，与惯性虚铰连接的接口指针为s+1。在构型参数的各个参数中，臂杆序号和允许连接数在臂杆出厂前确定并写入软件中，接口指针根据在轨实际连接的臂杆编号自主检测确定。

图4 构型参数的结构Fig.4 Structure of configuration parameters

具体地，以图1中k连杆为例，其臂杆序号为k；其接口数为2；k连杆有两个接口，分别与k-1连杆、k+1连杆连接，所以其接口指针有两个，分别为k-1和k+1。k连杆的构型参数可描述为下图所示的结构。

图5 构型参数示例Fig.5 Sample of configuration parameters

2.4 拓扑描述矩阵自生成方法

多臂航天器的拓扑构型可用基于图论的关联矩阵S和通路矩阵P进行描述。关联矩阵可以说明铰hk的内接体和外接体，对于n个连杆组成的树系统，关联矩阵为(n×n)阶矩阵，其元素定义为

(12)

通路矩阵P描述由体Bi到根体Bn的通路上有哪些铰，对于n个连杆组成的树系统，其也为(n×n)阶矩阵，元素定义为

(13)

在多体系统中，通常定义Bi体的内接铰即为hi。所以通路矩阵和关联矩阵的对角线元素均为1。由通路矩阵和关联矩阵的定义，得如下关系：

S·P=E

(14)

式中：E为n阶单位方阵。

在建立构型描述矩阵的过程中，首先根据编号找到根体。再以分支为单位，从根体向末端开始遍历连杆的构型参数。根据与根体的连接关系，确定每个臂杆的内接体和外接体，从而生成关联矩阵S，对关联矩阵求逆，生成通路矩阵P。

由于多臂航天器在夹持工况下可能存在闭环拓扑结构，需要对其进行闭环检测。具体地，在某一分支的遍历过程中，如某一连杆的外接体为根体，则系统中产生了闭环，同时继续进行遍历。当系统遍历完成，对检测到的闭环进行拆分，对拆分后的系统重新遍历，生成拆分后系统的通路矩阵、关联矩阵。特别地，在夹持外部物体时，外部物体自身无法生成构型参数，需要在夹持操作完成后人为赋值。

上述算法流程如图6所示。

图6 算法流程图Fig.6 Flowchart of the algorithm

3 基于SOA与构型矩阵的动力学建模

3.1 基于构型矩阵的动力学建模方法

在空间算子代数建模理论中，只有空间转移算子φ和科氏力算子a与系统拓扑连接关系相关。而其中，科氏力算子a又可根据2.4节算法输出的内接体序列、外接体序列直接得到。所以根据拓扑构型自主快速生成空间转移算子φ是变拓扑多臂航天器动力学建模的关键所在。

从式(8)可看出，对于空间转移算子的对角线元素，其元素恒存在且不为全零矩阵；对于非对角线元素φj,k，如果j在k到根部的路径上，则元素存在，否则为全零矩阵。这与拓扑描述矩阵中的通路矩阵定义是一致的：通路矩阵的对角线元素恒存在；非对角线元素，如果i在j在到根部的路径上，则元素存在，否则元素不存在。

据此，本文提出变拓扑系统空间转移算子生成思路如下：首先，定义虚拟转移算子φo，用以表征多体系统内任意两个连杆之间均互为内接体的情况下，该系统的空间转移算子。这种多体系统实际上是不存在的，故称为虚拟转移算子。其次，将通路矩阵改写为一分块选择矩阵Ps。具体的，将通路矩阵所有1元素修改为6阶单位阵；将所有0元素修改为6阶全零方阵。最后，将选择矩阵Ps与虚拟转移算子φo作哈达玛积(Hadamard product)，即对应分块子矩阵相乘，得到当前拓扑构型下系统的转移算子φr。虚拟转移算子φo、选择矩阵Ps、真实转移算子φr的表达式如下：

(15)

(16)

(17)

式中：{i,j}表示分块矩阵第i行第j列的分块子矩阵；⊙表示哈达玛积，即对应分块子矩阵相乘。

在基于SOA的动力学方程中，除空间转移算子外，其他量均与拓扑构型无关，故将真实转移算子代入动力学方程，即可完成系统动力学建模。

传统动力学建模中，要求连杆编号满足末端到基座递增或递减的规则。但本方法采用了完整的虚拟转移算子，经选择矩阵处理后的转移算子可视为传统算子通过行列变换得到，因此对于连杆编号没有严格要求，仅要求每个连杆编号固定。这种编号性质正适合于多臂航天器拓扑构型复杂多变的特点。

3.2 臂杆不对称情况下的动力学建模

在空间算子代数建模理论中，使用到了每个臂杆的惯性矩阵M(k)，其表达式为

(18)

但在实际工程实践中，机械臂的臂杆往往为非对称臂，其质心通常不在机械臂中心，不同本体坐标系下Ik和Pk(k)均不一致。在SOA动力学建模理论中，每个连杆的本体系坐标原点定义在内接铰中心，坐标系的定义与拓扑构型相关，由此导致系统惯量参数变化。例如，对如图7所示的臂杆，当铰k为该臂杆的内接铰时，其本体坐标系为k系，质心向量为b；当铰k+1为该臂杆的内接铰时，其本体坐标系为k+1系，质心向量为a。惯性张量矩阵同理。

图7 臂杆坐标系及质心示意Fig.7 The coordinate and the center of mass of the link

针对上述因拓扑构型变化，导致的动力学参数不一致的问题。本文提出一种基于构型参数的动力学参数选择方法。具体地，预先分别将不同坐标系定义下的质心向量、惯性张量矩阵储存至计算机中。在动力学建模过程中，根据图6所示算法输出的内接体序列，确定与每个臂杆的内接铰，确定本体坐标系定义，再根据坐标系选择对应的惯量参数。

3.3 含闭环拓扑系统的动力学建模

对于含闭环拓扑系统，将系统在闭环处进行分割，将系统拆分为树系统。对于分割得到的树系统，采用SOA进行动力学建模，再对动力学方程附加闭环约束方程，得到含闭环系统动力学方程[22]。

具体地，对于闭环结构，定义闭环约束为

(19)

对式(19)应用虚功原理，得到闭环约束力为

(20)

将式(19)和式(20)合并至系统动力学方程，可得

(21)

式中：

(22)

(23)

在动力学求解过程中，首先通过式(23)的第二行求解约束力，再通过第一行求解系统动力学。

对于多个多臂航天器组成的环链混合系统，同样地，将系统在闭环处进行分割，将系统拆分为树系统，进行树系统动力学建模，再附加闭环约束。与单个航天器的闭环系统相比，多航天器的环链系统存在多个中心体，需要将其中一个中心体视作整个树系统的根，其他中心体视作树系统的连杆，对整个树系统进行动力学建模。

3.4 对多柔体系统的推广

(24)

式中：φ(tk,k+1)为k体经有限元离散后与k+1铰连接的元素到k+1体的转移算子；Km(k)为柔性体的刚度矩阵；ϑ(k)表示k体模态坐标和铰坐标组成的广义坐标；下角标m表示柔性体。式(24)写成算子形式为

(25)

(26)

对式(25)、式(26)与式(2)～式(6)，显然相比于多刚体系统，多柔体系统的动力学模型结构并没有发生改变，式(10)～式(12)的关系仍然成立。因此，对于多柔体系统，本文在3.1节和3.2节提出的动力学建模方法仍适用。

4 仿真验证

4.1 仿真条件

为验证上述动力学建模方法的准确性和有效性，本文分别使用数值仿真软件和多体动力学仿真软件建立了双臂航天器的动力学模型。基本构型(构型1)下双臂航天器如图8所示，双臂航天器由一个中心体和两个三自由度机械臂组成，航天器参数如表1所示。为验证双臂航天器变拓扑动力学模型，本文设计了几种航天器构型变化，其拓扑构型示意如图9所示。为简化研究，该双臂航天器仅可进行二维运动，且基座固定。SOA理论在漂浮基机器臂及三维空间内运动领域已有大量研究，可认为SOA理论本身是准确有效的。且本文研究重点为系统的变拓扑特性，并没有对系统的变拓扑特征进行简化。因此，本文对空间双臂航天器所作的简化是合理的。

图8 双臂航天器示意Fig.8 Diagram of a dual-arm spacecraft

图9 双臂航天器的拓扑构型变化示意图Fig.9 Diagram of topology configuration changes of a dual-arm spacecraft

表1 航天器参数表Table 1 Parameters of the spacecraft

4.2 树状系统动力学仿真

在数值仿真软件中，采用脚本文件的形式，基于本文提出的动力学理论，分别针对基本构型、构型2、构型3建立其数值仿真模型。相应地，在多体动力学仿真软件中，分别建立不同构型的多臂航天器动力学模型。假定仿真前，所有关节均处于零位，初始加速度为0。在仿真过程中，分别对6个关节施加q=45°·sin(π·t/2000)的关节角运动，其中t为仿真时间，单位为s，比较各关节所需输出的关节力矩。以多体动力学仿真软件的仿真结果为基准，验证算法的准确性。算法生成的通路矩阵如式(27)～式(29)所示，式(27)对应图8所示拓扑构型，动力学仿真结果如图10所示；式(28)对应图9(a)所示拓扑构型，动力学仿真结果如图11所示；式(29)对应图9(b)所示拓扑构型，动力学仿真结果如图12所示。

(27)

图10 构型1动力学仿真结果Fig.10 Dynamics simulation results of Configuration 1

(28)

图11 构型2动力学仿真结果Fig.11 Dynamics simulation results of Configuration 2

(29)

图12 构型3动力学仿真结果Fig.12 Dynamics simulation results of Configuration 3

从式(27)～式(29)中可看出，对于不同拓扑构型，本文所提方法均能准确生成相应的通路矩阵。从仿真结果可看出，数值仿真软件的仿真结果能准确地描述系统的动力学响应。对比多体动力学仿真软件的仿真结果，数值仿真软件的仿真误差与实际结果相差6个数量级，几乎可以忽略不计。因此，本文所提出的变拓扑系统动力学建模方法可应用于多种树状拓扑系统构型，建模方法准确有效。

4.3 闭环系统动力学仿真

为验证上述建模方法对于闭环系统动力学建模的有效性和准确性，分别建立了如图13所示的闭环动力学系统的数值仿真模型和多体动力学仿真模型，航天器参数如表1和表2所示。同样地，为简化模型，验证建模方法对不同拓扑构型的有效性，将航天器基座固定。

图13 闭环构型示意Fig.13 Diagram of the close-loop configuration

表2 操作目标参数Table 2 Parameters of the target

仿真前，所有关节角均为-45°，初始速度、初始加速度均为0。在仿真过程中，对其中5个关节施加0.1 N·m的恒定力矩，其余关节自由，比较关节的角加速度。以多体动力学仿真软件的仿真结果为基准，验证算法的准确性。仿真结果如图14所示，生成的通路矩阵如式(30)所示。

(30)

图14 闭环构型动力学仿真结果Fig.14 Dynamics simulation results of the close-loop configuration

从式(30)中可看出，对于如图13所示的闭环系统，本文提出的算法将系统在连杆3与连杆8中间断开，变成与图9(b)类似的单链系统，完成了由闭环系统向树状系统的转换与描述。从仿真结果可看出，数值仿真软件的仿真结果能准确地描述系统的动力学响应。对比多体动力学仿真软件的仿真结果，数值仿真软件的角加速度误差曲线为规则曲线，这可能是积分误差导致的。但误差与实际结果相差3个数量级，几乎可以忽略不计。因此，本文所提出的变拓扑系统动力学建模方法可应用于闭环拓扑系统，建模方法准确有效。

5 结 论

针对多臂航天器可能存在的开链、闭环、环链混合等多种拓扑构型及相应动力学模型变化的问题，本文提出了基于图论和空间算子代数理论的变拓扑多体系统动力学统一建模方法，且具有O(n)阶的高计算效率。本文提出多臂航天器构型参数的概念，用于描述航天器臂杆间的连接关系，并生成用于描述拓扑构型的关联矩阵和通路矩阵；将通路矩阵与空间转移算子结合，考量多臂航天器实际特性，基于空间算子代数方法建立树系统和闭环系统的动力学模型，该建模方法适用于开链、闭环、环链混合等多种拓扑构型，并可扩展到多柔体系统，是一种具有统一形式的高效的动力学建模方法。数值仿真软件的仿真与多体动力学仿真软件的仿真结果对比表明：本文所提方法对于不同拓扑构型的多臂航天器均能实现较高精度的建模，建模有效，仿真结果准确。本文所提出的多提系统变拓扑统一建模方法对未来多臂航天器的设计与应用具有一定的理论参考价值与工程指导意义。