手臂的冗余特性导致其可达臂位形并不唯一．经过长期的自然进化，人类已形成一种潜在的臂位形选择策略使手臂具有更高的灵巧性、可操作性和舒适度．研究该策略可以有效解决虚拟人运动复现[1]、仿人机器人的运动规划[2]、康复外骨骼机器人仿人运动[3]等问题．目前国内外学者对于人臂可达运动臂位形预测的研究主要是基于指标的方法和基于学习的方法．基于指标的方法有基于单指标的方法，如关节力矩最小指标[4]、能量最小指标[5]、势能最小指标[6]、操作椭球指标[7]、快速上肢评价指标[8]和基于运动基元方法[9]等；基于混合指标的方法，如最小化手部接触力和力矩凸结合指标[10]、关节力矩最小和关节平滑度最大凸结合的指标[11]、能量、测地线、关节力矩、笛卡尔冲击、角加速度加权相加指标[12]等．文献[13]使用最小化肌肉力量指标分析人臂位形，并将其用于仿人机械的运动规划中．文献[14]使用fPCA方法从人臂日常活动中提取主要的运动模式，用于理解人臂可达运动中的臂位形．基于指标的方法通常伴随着运动生成的假设，这可能会限制运动的可变性且计算相对复杂．基于学习的方法中，文献[15]采用最小二乘法拟合手掌中心位姿与臂角的非线性关系来研究人臂运动特性．文献[16]采用多元线性回归的方法学习人臂末端位姿与特定关节角度的关系，实现人臂可达运动中的自然臂位形．文献[17]使用人工神经网络(ANN)学习人臂末端位姿与臂角非线性关系，在遥操作任务中实现仿人机械臂的仿人臂位形，但其预测性能有限．文献[18]使用深度卷积神经网络(DCNN)改进ANN学习模型提高了预测性能．现有基于学习的方法并没有考虑上下臂长度对于人臂可达运动臂位形的影响，因此普适性差、精度较低．高斯过程回归(Gaussian process regresssion，GPR)具有易实现、参数自适应推断、自带特征选择能力等优点，对处理高维、非线性、小样本等复杂问题具有更好的适应性[19]．本研究采用GPR的方式学习人臂可达运动臂位形，并将其用于虚拟人可达运动复现和机械臂的仿人运动规划，同时考虑上下臂长度对臂位形的影响，以提高预测精度和适用性．1 人臂可达运动数据采集与臂位形1.1　人臂运动学模型如图1所示，将人臂(左臂)简化为包括肩关节、肘关节、腕关节的7自由度运动学模型．关节1~3轴线相互正交为肩关节；关节4为肘关节；关节5~7轴线相互正交为腕关节．qk为第k个关节的角度；S为肩关节；E为臂平面轴肘关节；E0为参考10.13245/j.hust.210104.F001图1人臂运动学模型及臂角定义平面肘关节；W为腕关节．采用臂角作为手臂的冗余参数，分析人体利用手臂冗余性的能力．臂角φ为参考平面和臂平面之间的夹角[13]，由肩关节、肘关节、腕关节组成的平面为臂平面，关节3角度为0°时的臂平面为参考平面．1.2　人臂可达运动数据获取使用Vicon红外运动捕捉系统采集人臂运动数据．选择4名身高(168~185 cm)、体重(60~85 kg)不同的人作为受试人员．如图2(a)所示，在受试人员肩关节(S)、肘关节(E)、腕关节(W)、手掌中心(H)分别布置标记点．10.13245/j.hust.210104.F002图2人臂可达运动数据获取为了模拟人臂日常行为，以握手柄为例进行实验(如图2(b)所示)．将云台依次放在4个不同位置，在每个云台位置通过手柄刻度随机确定10个手柄位姿，共计40个不同手柄位姿．4名受试人员分别以最舒适的方式手握40组手柄位姿，使用红外运动捕捉系统采集受试人员S，E，W，H的位姿及上下臂长度，共计160组数据．受试人员以正常坐姿身体紧贴座椅左侧，通过坐垫的增减保证不同身高的受试人员肩关节位于相同位姿．臂角满足如下关系，φ=sgn φarccos((e×w)(e0×w0)||e×w||∙||e0×w0||)；sgnφ=-sgn(e0×w0∙e)，式中：e为臂平面中从肩关节到肘关节的向量；w为臂平面中从肘关节到腕关节的向量；e0为参考平面中从肩关节到肘关节的向量；w0为参考平面中从肘关节到腕关节的向量．1.3　臂长对可达运动臂位形的影响为验证臂长对可达运动臂位形的影响，计算并分析160组数据对应的臂角．在相同手柄位姿，臂角的最大差值为0.73 rad，最大标准差为0.3 rad．4个手柄位姿处3名不同受试人员的臂角位形如图3和表1所示．图中：绿色为1号人员；红色为2号人员；蓝色为3号人员．当建立人臂可达运动臂位形模型时，有必要考虑上下臂长度．10.13245/j.hust.210104.F003图3受试人员手握手柄臂位形10.13245/j.hust.210104.T001表1不同人员在相同手柄位置臂角手柄位置1号人员2号人员3号人员标准差a-1.43-1.19-1.080.15b-1.32-1.26-0.840.21c-1.12-0.92-0.850.11d-0.96-0.59-0.230.30rad2 预测模型及训练2.1　预测模型人臂可达运动中，臂位形预测的数学模型可以概括为g:m→φ，式中：m=[x,y,z,α,β,λ,d1,d2]T为人臂可达运动中手掌中心位置及上下臂长度组成的列向量，其中[x,y,z]为手掌中心位置，[α,β,λ]为手掌中心姿态的欧拉角，d1为上臂长度，d2为下臂长度；φ为该位置对应的臂角；g为可达运动中从手掌中心位姿及上下臂长度到相应臂角的映射关系，该映射关系是一种典型的高维、非线性关系．基于GPR人臂可达运动臂位形预测流程如图4所示．10.13245/j.hust.210104.F004图4基于GPR人臂可达运动臂位形预测流程从160组数据中随机选择120组数据的手掌位置及上下臂长度作为训练输入集M，相应臂角作为训练输出集ψ；剩余40组数据作为测试输入集M*，相应臂角作为测试输出集ψ*．建立GPR人臂可达运动臂位形预测模型训练集和测试集分别为D=(M,ψ),D*=(M*,ψ*)．在数据获取中通常会包含红外摄像头的噪音，将该噪音考虑到未知输入中，有φ=g(m)+ε，式中ε为独立于g(m)的高斯白噪音，且服从高斯分布．为提高计算效率，将均值函数设为0均值函数，GPR可写成φ~G(0,k(mi,mj)+σn2I)，式中：G(⋅)为高斯过程；mi和mj为训练输入集M中的样本点；k为协方差函数； σn2为噪音方差．对于测试集中一个测试点(m*,φ*)，φ*与训练数据ψ的联合分布为ψφ*~N0,  K(M,M)+σn2IK(M,m*)K(m*,M)K(m*m*)，式中：K(M,M)为训练数据的自相关协方差矩阵；K(M,m*)=K(M,m*)T为训练数据与测试点之间的互相关协方差矩阵；K(m*,m*)为测试点的自相关协方差矩阵．φ*的后验概率φ*|M,ψ,m*~N(φ¯*,cov(φ*))，式中：φ¯*为可达运动臂位形预测模型在测试点m*处的臂角预测值；cov(φ*)为方差．协方差函数k(mi,mj)将在模型训练阶段确定．2.2　预测模型训练确定适合的协方差函数的是基于GPR人臂可达运动臂位形预测的关键与难点．通过对人臂运动训练集D及不同协方差函数预测效果的对比，将人臂运动数据的先验知识归纳为：手臂末端位姿及上下臂长度8维特征输入对臂角的影响不同、总体趋势相对平滑、局部相关性强并且具备一定的全局相关性．选择各向异性Matern协方差函数、各向异性SE协方差函数、各向异性RER协方差函数和Noise协方差函数组成的复合协方差函数对人臂可达运动臂位形进行建模k(mi,mj)=k1(mi,mj)+k2(mi,mj)+k3(mi,mj)+σn2δ;    k1(mi,mj)=σf12∑n=18(1+3|mi-mj|/lnM)∙exp(-3|mi-mj|/lnM);k2(mi,mj)=σf22exp-12∑n=18((mi-mj)/lnS)2;    k3(mi,mj)=σf32∙∑n=18exp-2sin2(π(mi-mj)pn)/lnPR2,式中：σf12,σf22和σf32分别为Matern，SE，PRE协方差函数的信号方差；lnM，lnS和lnPR分别为Matern，SE，PRE协方差函数在第n个维度上的特征长度；δ为克罗内克函数，当i=j时δ=1，当i≠j时δ=0；pn为第i个维度上的周期，n=1,2,…,8为人臂运动数据的8维特征输入；(σf12,σf22,σf32,lnM,lnS,lnPR,pn,σn2)为系统超参数．目前求解最优超参数的方法有多种[19]，本研究采用最大似然估计法学习最优超参数，根据训练集D学习的最优超参数为：(σf12,σf22,σf32,σn2)=(-1.84,-0.66,-0.30,-1.15)；l1Ml1Sl1PRp1l2Ml2Sl2PRp2…l8Ml8Sl8PRp8=-2.08-1.072.893.56-0.91-0.75-1.993.03…-6.760.433.482.00．3 实验验证进行四组实验：前两组实验使用Matlab 2019b，通过与基于最小二乘法的人臂可达运动臂位形预测模型、不考虑上下臂长度的GPR人臂可达运动臂位形预测模型进行对比，验证本文方法有效性；后两组实验分别将本文方法用于虚拟人手臂可达运动复现和协作机械臂仿人运动规划，验证其优越性．3.1　人臂可达运动臂位形预测模型有效性验证使用残差(ε)、绝对平均误差(MAE)、均方根误差(RMSE)作为评价标准对模型进行评估．将测试集D*带入基于GPR的人臂可达运动臂位形预测模型，预测结果如图5所示(N为测试点序号)．10.13245/j.hust.210104.F005图5基于GPR人臂可达运动臂位形预测模型预测结果由图5可见：基于GPR预测的臂角与真实臂角整体趋势相符合，95%置信区间整体较窄且真实值都落在置信区间内，表明本模型具有较高预测精度．将测试集D*代入基于最小二乘法的人臂可达运动臂位形预测模型[15]，与本文方法进行对比，结果如图6所示．10.13245/j.hust.210104.F006图6臂长对可达运动臂位形预测模型性能影响从图6(a)中可以看出：基于最小二乘预测的臂角整体趋势与真实臂角相差较大；从图6(b)中可以看出：基于最小二乘法预测的臂角残差落在(-0.1,0.1)rad，占17.5%，落在(-0.2，-0.1)~(0.1，0.2) rad，占70%．两种方法的评估结果如表2所示，可以看出：基于GPR预测臂角的MAE和RMSE值远小于最小二乘法预测的臂角．10.13245/j.hust.210104.T002表2两种方法评估结果方法MAERMSEGPR0.053 90.004 5最小二乘法0.141 10.022 33.2　臂长对可达运动臂位形预测模型性能影响分别使用m=[x,y,z,α,β,λ,d1,d2]T和m'=[x,y,z,α,β,λ]T作为输入，训练高斯过程回归臂位形预测模型，结果如图7所示．10.13245/j.hust.210104.F007图7虚拟人手臂可达运动复现从图7(a)中可以看出：不考虑上下臂长度的预测模型整体趋势与真实臂角相差较大．从图7(b)中可以看出：不考虑上下臂长度预测模型预测的臂角残差落在(-0.1,0.1)rad，占27.5%，远低于考虑上下臂长度预测模型的预测结果．臂长对可达运动臂位形性能影响的评估结果如表3所示，可以看出：基于GPR预测臂角的MAE和RMSE值远小于最小二乘法．10.13245/j.hust.210104.T003表3臂长对可达运动臂位形性能影响的评估结果方法MAERMSE考虑臂长GPR0.053 90.004 5不考虑臂长GPR0.147 50.028 93.3　虚拟人手臂可达运动复现基于GPR人臂可达运动预测模型可用于虚拟人手臂可达运动复现．如图8所示，使用vs2017和OpenInventor搭建仿真环境，在仿真环境中建立虚拟人体模型，虚拟人身高1.8 m，上下臂长度分别为0.33和0.22 m．10.13245/j.hust.210104.F008图8仿真环境搭建及预测使用运动捕捉系统采集受试人员手掌中心相对于肩关节位姿，作为虚拟人手部可达位姿．带入虚拟人手部位姿和上下臂长度，使用本文模型预测虚拟人可达运动手臂位形，虚拟人臂与实际人臂关节角度误差如表4所示．可见，预测的虚拟人可达运动的臂位形与实际人臂的位形基本相同．10.13245/j.hust.210104.T004表4虚拟人臂与实际人臂关节角度误差位置关节1关节2关节3关节4关节5关节6关节710.51.1-3.80-2.1-1.1-2.121.61.6-3.20-3.2-1.60.0(°)-13.4　机械臂仿人运动规划基于GPR人臂可达运动预测模型也可用于冗余机械臂的仿人抓取规划．如图9(a)，在KUKA iiwa 7自由度冗余机械臂水杯抓取规划实验中，已知机械臂初始臂位形和抓取目标位姿；使用已训练的GPR人臂可达运动预测模型，代入机械臂上下臂长度和抓取目标位姿，获得抓取目标的仿人臂位形(图9(c))，随后在机械臂的位形空间中采用基于随机采样的路径规划方法(如RRT*算法)生成无碰撞的抓取路径，规划结果如图9(b)所示．基于GPR人臂可达运动预测模型可生成更为拟人的目标臂位形，提高了整个抓取规划的拟人程度．10.13245/j.hust.210104.F009图9机械臂仿人运动规划4 结语使用实验的方法验证了上下臂长度对可达运动中臂位形的影响，将人臂可达运动的手掌位姿和上下臂作为特征输入建立基于高斯过程回归人臂可达运动臂位形预测模型．实验结果表明：基于高斯过程回归的人臂可达运动臂位形预测方法的精度和有效性远高于现有人臂位形预测方法．该模型可有效解决虚拟数字人、仿人机械臂、康复外骨骼机器人的仿人运动规划等问题．