航天员虚拟训练中的一维直线运动速度规划算法

王 岚，林凌杰，常 影，薛 峰

(哈尔滨工程大学机电工程学院，哈尔滨 150001)

0 引 言

航天员在太空中进行任务时，感受不到重力的存在，但惯性力依然存在[1-2]。如果物体的质量较大，则航天员很难保证其在目标点的速度为零。在作业过程中，较大的惯性力可能会使航天服手套撕裂，因此有必要在地面对航天员进行模拟训练。

本文在文献[1]所研究的基于虚拟现实技术的航天员作业训练环境下，在航天员完成搬运或安装某一目标物体的过程中，对目标物体的速度进行规划。从而对航天员进行有效的示教，提高航天员的训练速度，同时可以为航天员作业训练系统的模拟仿真和建立“虚拟人”模拟仿真系统提供参考。

速度规划算法种类繁多，最常使用的是直线型加减速控制算法、S型曲线加减速控制算法等。直线型加减速算法是最简单易行的算法，但其加速度存在突变点，导致冲击较大[3]。S型曲线保证了速度和加速度的连续，冲击较小。文献[4-6]中的方法均假设允许的最大速度和最大加速度为常量，通过相平面法，得到各阶段的运行时间的解析解。

近年来，直接方法成为时间最优速度规划的研究重点，该方法通过将时间连续模型离散化，将时间最优速度规划转化为一个静态最优化问题，得到时间最优速度规划的数值解[7]。杨亮亮等[8]将关于加速度变化时间的方程组化简为一元高次方程后，转化为单一凸型函数，进而用牛顿迭代法进行求解。Lin[9]提出了一种基于最小加加速度的关节机器人轨迹规划算法，该算法将轨迹离散化之后，采用粒子群算法和K-means聚类，求出各段轨迹所需要的时间，但每段轨迹内的加加速度均包含t-3项，当t过小或轨迹分段的数量较多时，迭代后的加加速度不易收敛，降低算法的稳定性。文献[10]和文献[11]分别提出了基于遗传算法和粒子群算法的速度规划算法，将最大加速度和最大速度简化为常量。Liang等[12]采用迭代的方式，从零开始逐步提高B样条曲线的各个控制点，直至逼近约束条件。Yuan等[13]提出了一种新型的前向和后向检查算法，通过多个解析方程，直接求出各节点的最优解。以上算法均假设在进行速度规划之前，待规划点的加速度约束为已知量或常量，而对于本文来说，只有对速度进行规划之后，才能获得加速度约束。

Bharathi等[14]将轨迹离散化之后，在每个节点处求出允许的最大速度，实现了高阶动态约束，但其算法在每个节点处都通过二分法求解，增加了求解时间。由于所研究的对象不同，上述算法均无法直接应用于本文所研究的问题中，目前还没有一种考虑人体生物力学模型的速度规划算法。

本文从运动生物力学的角度出发，建立了手臂的动力学模型及肘关节的骨骼肌模型，分析了航天员搬运物体的不同情况，并对肌肉激活度进行规划。通过粒子群算法求解出各个阶段的时间，实现了目标物体的时间最优速度规划。

1 人体上肢生物力学模型

1.1 人体上肢动力学模型

在太空中，航天员在搬运质量非常大的物体时，会在反作用力的作用下产生运动。最好的搬运方法就是让物体在航天员的矢状面内做一维直线运动，此时航天员只会发生前后的平移，而不会发生旋转，容易控制住自己的身体。航天员通过调整位置，将复杂的搬运作业分解为若干个矢状面内垂直于冠状面的一维直线运动。因此，文献[1]和文献[15]均将航天员运送物体的作业训练视为在水平面内做卧推运动，建立推拉物体模型如图1所示。图中，l5和l1分别表示左、右上臂的长度，l4和l2表示分别表示左、右前臂的长度，l6表示两个肩关节之间的距离，l3表示左右两手之间的直线距离。m1,m2,m4,m5分别表示四肢的质量，M代表目标物体的质量，T1,T2,T5,T6分别表示右肩关节、右肘关节、左肘关节、左肩关节产生的力矩，θ1表示肩关节角度，θ2表示肘关节角度。lp表示人体腕关节和肩关节在冠状面上的投影之间的距离，F′l和Fl分别为左右两臂承受的物体的惯性力。

取θ1,θ2为广义坐标，T1,T2为广义力，对右臂列拉格朗日动力学方程[16]可求得

Fl[l1cosθ1+l2cos(θ1+θ2)]

(1)

(2)

式中：

h112=h122=-h221=-m2l1lc2sinθ2。

其中，lc1,lc2分别为上臂的质心到肩关节的距离和前臂的质心到肘关节的距离，J1,J2分别为上臂和前臂的转动惯量。

图1 推拉物体模型Fig.1 Push-pull object model

如果物体只做矢状面内垂直于冠状面的一维直线运动，则左臂和右臂承受的惯性力相同，可求得

(3)

式中：M为物体的质量，a为物体的加速度。

人体腕关节和肩关节在冠状面上的投影之间的距离

(4)

将式(3)、式(4)分别代入式(1)和式(2)中可求得由肩关节和肘关节决定的物体最大加速度的绝对值分别为

(5)

(6)

物体的最大加速度

amax=min(a1max,a2max,aglove)

(7)

式中：aglove是航天服手套被撕裂的加速度。

与物体的质量相比，关节的角速度、角加速度、人体上臂的质量、转动惯量都很小，对式(5)和式(6)的影响非常小。因此，a1max和a2max的大小主要取决于|lp|和|l2cos(θ1+θ2)|的大小。图2表示出了右臂可能的三种位姿。由于航天服的限制，取θ1的范围为[30°,75°]，取θ2的范围为[10°，120°]。

图2 右臂存在的三种位姿Fig.2 Three postures of the right arm

当点C在冠状面的投影位于A点以左，由于θ1小于90°，因此在整个运动过程中，航天员的右臂只存在图2(a)一种情况，此时，由图中可以看出，lp≤0时，l2cos(θ1+θ2)<0，且|lp|一定小于|l2cos(θ1+θ2)|。

当点C在冠状面的投影位于A点右侧时，可能会出现图2(b)和图2(c)两种情况。如果在θ1取最小值时，点C在冠状面的投影也位于B点以右，那么在整个运动过程中，点C在冠状面的投影始终位于B点以右，即始终是图2(b)中的情况。此时，θ1+θ2<90°，l2cos(θ1+θ2)>0，|lp|一定大于|l2cos(θ1+θ2)|。此种情况出现的条件为lp>l1cosθ1 min，取l1=280mm，可求得lp>242.5mm。

如果在θ1取最小值时，点C在冠状面的投影也位于B点左侧，那么随着B点的运动，点C在冠状面的投影有可能变为B点右侧。如果在θ1取最大值时，点C在冠状面的投影也位于B点以左，那么在整个运动过程中，点C在冠状面的投影始终位于B点以左，即始终是图2(c)中的情况。此时，l2cos(θ1+θ2)<0，lp>0。此种情况出现的条件为lp

由以上分析可知，若lp<36.2 mm时，可以保证在运动过程中，|lp|始终小于|l2cos(θ1+θ2)|，a2max小于a1max。在实际的训练设备中，lp的值远小于36.2 mm，因此，物体的最大加速度主要受到肘关节力矩T2和航天服手套不被撕裂的限制。

1.2 肘关节及其骨骼肌模型

为了建模方便，将与肘关节有关的肱二头肌(BIC)和肱三头肌(TRI)等效为单头肌肉，通过其等效长度确定肌肉上端的等效位置点，其等效长度的计算方法采用文献[17]的计算方法。

根据Hill肌肉模型，肌纤维力Ff：

Ff=FCE+FPEE

(8)

式中：FCE为肌纤维主动力，FPEE为肌纤维被动力。

由文献[18]可知，肌纤维主动力与肌肉激活度A、肌纤维的长度lM及肌纤维收缩速度vM相关，即肌纤维主动力FCE：

FCE=F(A,lM,vM)

(9)

肌肉激活度A反映的是肌肉在神经系统的刺激下受到激活的程度，激活度越高，肌肉产生的主动收缩力也越大。在一定程度上，肌肉的激活度可以从肌电信号角度来表征，通过采集该肌肉处的皮肤表面肌电信号，并做归一化处理得到肌电信号强度，以肌电信号强度作为输入，经过肌肉激活度模型，来表征肌肉激活度。其取值范围为[0，1]。

肌纤维被动力只与肌纤维的长度有关，即肌纤维被动力：

FPEE=F(lM)

(10)

肘关节的骨骼肌模型如图3所示。图中，d1,d2分别为肱二头肌两个止点到肘关节旋转中心的距离，d3为肱三头肌上止点到肘关节旋转中心的距离，θ2表示肘关节的角度。肱三头肌是沿肘关节的外包络曲面分布的，将该曲面简化为半径为R2的球面。

图3 肘关节骨骼肌模型Fig.3 The skeletal muscle model of elbow joint

肌肉拉力线可以近似为直线，根据图3的几何关系可得，肱二头肌的肌肉长度lBIC：

因此，肌纤维长度lMBIC：

(11)

式中：lTBIC为肱二头肌的肌腱长度之和，可以认为是常量；φBIC为肱二头肌的羽状角。

肱三头肌的肌纤维长度lMTRI：

(12)

式中：θ0为肱三头肌处于最大等长收缩状态时肘关节的关节角度，φTRI为肱三头肌的羽状角。

对于肌力臂的计算，采用虚功的方法求解动态肌力臂[19]，解得肱二头肌肌力臂dBIC：

(13)

肱三头肌的肌力臂dTRI：

(14)

肘关节的总力矩T2:

T2=FfBICdBIC+FfTRIdTRI-Tf

(15)

式中：Tf为航天服的肘关节阻尼力矩，其值由文献[20]确定，如下式所示。

(16)

将式(8)-式(14)代入式(15)中，即可求得肘关节的总力矩与肘关节角度之间的函数关系式

(17)

2 基于肌肉激活度的速度规划

将式(17)代入式(6)、式(7)中可求得物体的最大加速度

(18)

由图1的几何关系可知，肩关节角度θ1与肘关节角度θ2存在一定的约束关系。可将式(18)简化为

(19)

匀加速阶段以最大加速度进行运动，根据式(19)可知，最大加速度为关节角度、角速度的函数，即在该阶段的任意t时刻时物体被移动的距离

(20)

式中：t1为匀加速阶段开始的时间，l0为匀加速阶段开始时物体到冠状面的距离。

形如式(20)的积分很难求出解析解，因此传统的相平面法无法完成此类问题的速度规划。

2.1 匀速阶段存在与否的分析

七段S型加减速曲线中，匀速段的存在条件是物体达到了最大速度。由图1和图3的几何关系可得，由肱二头肌决定的目标物体的最大速度

(21)

式中：vmaxBIC为肱二头肌最大收缩速度。

由肱三头肌决定的目标物体的最大速度

基于此，在充分研究企业管理者特质、内部控制以及企业价值关系的过程中，要建立对应的检验模型，为依据模型就能对相关参数的关系进行分析，并且要结合实际情况分别建立对应的管理者检验模型，从而合理性考量三者的关系，进一步提升模型分析的时效性。

式中：vmaxTRI为肱三头肌最大收缩速度。

目标物体的最大速度与肘关节角度的关系如图4所示。在肘关节角度为0°时，由肱三头肌决定的最大速度为0，实际上，此时规划后的速度也一定为0。因此，在搬运物体的过程中，目标物体的速度很难达到最大速度，匀速阶段不存在。

图4 物体的最大速度与肘关节角度的关系Fig.4 The relationship between the maximum velocity of the target object and the elbow joint angle

2.2 肌肉激活度规划

对于人体上肢来说，物体的最大允许加速度未知，但肌肉激活度的范围为[0，1]，因此，通过调整肌肉激活度来改变物体的加速度是较为合适的方法，可将动态的加速度约束转变为静态的肌肉激活度约束。

肌肉激活度的变化有三种情况：增加、减小和不变。前两种情况显而易见是存在的。当肌肉激活度不变时，物体的加速度近似不变。若物体的加速度为0时，则物体处于匀速状态，此时肘关节力矩近似始终等于航天服阻尼力矩，肌肉激活度近似保持不变。由于肘关节的最大力矩大于航天服阻尼力矩，因此此时肌肉激活度一定不为1。由于匀速阶段不存在，因此这种情况不存在；若物体的加速度为非零定值时，考虑到时间最短的原则，该加速度一定为最大加速度，即达到了最大加速度限制。对于肌肉激活度规划来说，最大加速度限制就是肌肉激活度达到了1时的加速度，即匀加速阶段肌肉激活度一定近似为1。因为肱二头肌和肱三头肌是一对拮抗肌，其中某根肌肉收缩时，另外一根肌肉一定放松，因此其肌肉激活度可能会始终为0。

由上述分析可知，在肌肉激活度规划中，肌肉激活度的变化只存在四种情况：增加、减小、始终为0、始终为1。

根据式(19)可知，物体的加速度与时间的关系为

(22)

(23)

式中：y0为初始时刻物体到冠状面的距离。

图5 A、B、C、D点的位置Fig.5 Location of points A, B, C and D

根据起始位置、目标位置和航天服在充气后肘关节自然弯曲的角度θT0的不同情况，可将搬运过程分为表1所示的6种情况，其中的A、B、C、D点的位置如图5所示，此时θ2=θT0=30°。

6种情况下的起始位置和终止位置的肌肉激活度见表2，“+”表示肌肉激活度不为0。

表1 不同起止点的分类情况Table 1 Classification of different starting and ending points

表2 各情况下的起止点肌肉激活度Table 2 Muscle activation at start and stop points under different cases

通过分析起始位置和终止位置，就可得到各肌肉起始时刻和终止时刻的肌肉激活度，以Case 5为例，肱二头肌起始时刻的肌肉激活度

(24)

肱三头肌终止时刻的肌肉激活度

(25)

式中：θ0和θd分别为起始时刻和终止时刻的肘关节角度。

在实际的运动过程中，由于运动距离的不同，Case 1、Case 3、Case 4、Case 6又可各分为两种情况，共计十种情况。本文仅以Case 4为例，其肌肉激活度规划如图6所示。

图6 Case 4情况下的肌肉激活度规划Fig.6 Muscle activation planning in Case 4

Case 4a和Case 4b两种情况下肱二头肌和肱三头肌的肌肉激活度与时间的关系如式(26)和式(27)所示。

Case 4a:

(26)

Case 4b:

(27)

式中：J为肌肉激活度的变化率。

将式(26)和式(27)代入式(22)即可求得各个阶段的加速度。以Case 4a为例，在终止时刻，肱二头肌的肌肉激活度为aBICd，因此

aBIC0-Jt1+J(t2-t1)-J(t3-t2)=aBICd

(28)

令t1=Tt1，t2-t1=Tt2，t3-t2=Tt3，代入式(28)可得

因此，给定Tt1和Tt2，即可求得在该种情况下肱二头肌和肱三头肌的肌肉激活度与时间的关系。

2.3 粒子群算法求解

采用Matlab/Simulink搭建上肢模型如图7所示。

图7 上肢Simulink模型Fig.7 Upperlimb model by Simulink

根据约束条件可知，物体在最终时刻的速度为0，位移为目标物体的搬运距离Δy，即

(29)

(30)

式中：Δy=yd-y0，yd为目标物体的给定位置到冠状面的距离。

通过联立式(29)和式(30)，可以唯一确定Tt1和Tt2的值。在本文中，通过标准粒子群算法求解Tt1和Tt2的值。

将Tt1和Tt2的值输入到图7所示的模型中，获得最终时刻的速度v及位置L，以此构建多目标优化算法的适应度函数f(x)：

f(x)=q1|v|+q2|L-yd|

(31)

式中：q1和q2分别为速度误差和位置误差的权重。

根据适应度，更新种群的位置，反复迭代之后，即可求得最优解。

本文设计的算法程序流程图如图8所示。

图8 算法流程图Fig.8 The flow chart of the algorithm

3 仿真算例

在Matlab2018编程环境下验证本算法。取目标物体的质量M=200 kg，y0=0.3 m，yd=0.55 m。航天服手套被撕裂的加速度aglove=2m/s2，J=3，速度误差的权重q1=19，位置误差的权重q2=115，人体及肌肉的各项参数设定见表3，肱二头肌和肱三头肌的其他结构参数按文献[21]总结的数据设定，见表4。

可通过以上算法可解得结果如下：

Tt1=0.3838 s，Tt2=0.5666 s，Tt3=0.3696 s。

运动过程中，肱二头肌和肱三头肌的肌肉激活度Tf如图9所示，目标物体的加速度、速度、目标物体到冠状面的距离如图10所示。

表3 人体及肌肉参数Table 3 The parameters of the human body and muscle

表4 肱二头肌和肱三头肌的结构参数Table 4 Structual parameters of the biceps and triceps

图9 肌肉激活度曲线Fig.9 Muscle activation profile

图10 规划结果Fig.10 Planning results of the algorithm

根据图9和图10可知，经过本文算法的规划之后，肌肉激活度进行了匀速的增加或减小，在末端时刻，物体的速度、加速度均为0，位置达到了给定位置0.55 m。肘关节的输出力矩连续变化，没有发生突变。在整个运动过程中，物体的加速度连续变化，实现了较好的规划效果。

将此过程中的加速度和加速度限制amin,amax绘在同一张图中，如图11所示。由图可知，实际加速度可以较为贴合的靠近加速度限制，并没有超出允许的加速度范围。此过程中的肱二头肌和肱三头肌的肌肉力FBIC和FTRI如图12所示。

采用文献[11]中的方法进行求解，将其对加速度的约束条件变为本文中的加速度约束条件，限制最大加加速度Jmax=5m/s3，解得规划时间为1.324 s，而本文算法求得的时间为1.322 s。

图11 加速度限制曲线Fig.11 Acceration constraints

图12 肌肉力曲线Fig.12 Muscle force profile

文献[11]规划结果的加速度和加速度限制如图13所示。对比图13和图11可发现，本文所采用的方法能达到加速度限制的边界并维持在加速度边界一段时间，而文献[11]所采用的方法只能在加速度的两个极值点达到加速度边界，而没有加速度基本维持在加速度边界的阶段，因此求解出的规划时间略长。

图13 加速度限制曲线Fig.13 Acceration constraints

根据文献[11]的规划结果反解出肱二头肌和肱三头肌的肌肉激活度如图14所示。从图中可知，该方法无法对肌肉激活度进行约束，求解出的肱二头肌肌肉激活度在0.6 s左右出现转折点，即肌肉激活度的变化率出现变化。对于航天员来说，肌肉激活度的变化率维持不变是比较符合人体生物力学的做法。本文求解的肌肉激活度的变化率只在肌肉激活度为0或者1时发生变化，即肌肉激活度匀速增加，对于航天员来说是较为合适的搬运方法。

图14 肌肉激活度曲线Fig.14 Muscle activation profile

4 结 论

本文建立了上臂和前臂的动力学模型及肘关节的骨骼肌模型，分析了航天员在身着航天服时搬运物体的十种情况，并对肌肉激活度进行规划。将动态的加速度约束转变为静态的肌肉激活度约束，通过粒子群算法求解出各个阶段的时间，从而实现了目标物体的在动态加速度约束下的时间最优速度规划。在多人协同虚拟训练和建立“虚拟人”模拟仿真系统中，提供φ0机器人和虚拟人一个相对符合人体生物力学的给定速度，具有实际应用可行性。

下一步的研究包括：

1)对较复杂的空间六自由度物体的运动过程进行分析；

2)将作业时间和作业过程中消耗的能量作为优化目标进行多目标优化；

3)根据被训练者的感官体验，确定肌肉激活度的变化率J的取值；采集每个被训练者的肌电信号，对结果进行动态调整，以增加算法的适用性。