表1(Table 1)
2. 伦敦大学 国王学院, 伦敦WC2R 2LS
2. King's College London, University of London, London WC2R 2LS, UK
人体肩部的骨骼、肌肉、韧带和周围皮肤构成的肩部复合体(shoulder complex)闭环系统，是人体最为复杂的运动功能部位之一[1]。研究肩关节复合体运动特征，特别是骨骼系统的运动特征对于人体体能训练，肩部功能障碍诊断以及康复等具有重要意义。同时，正确描述肩部骨骼的姿态和运动行为是建立肩部复合体模型的关键[2]。
为了准确描述肩部骨骼的运动特征，国内外学者在肩部骨骼的机构学建模方面做了大量的研究工作，这些机构模型主要分为2类：开式链机构和闭式链机构[1]。在开式链机构方面，从1965年Dempster[3]建立了第一个简单的肩部二维串联连杆模型开始，Engin和Chen[4]于1986年又提出了一种将肱骨相对于人体躯干的运动描述成三维6自由度的刚体模型，之后Engin和Tumer[5]进一步将该肩部机构模型的自由度增加到9。Dempster[3]和Engin等[4-5]虽然建立了肩部骨骼的三维机构模型，但是这些机构模型都是从仿人机器人领域的需求出发而建立的，其目的是使模型具有与人体肱骨相近的活动度和运动范围，而非重现肩带内部骨骼的运动约束关系，即没有系统描述人体肩部复合体所有骨骼的运动本质特征。
进入20世纪90年代之后，随着光标位置检测技术[6-7]、电磁跟踪技术[8]以及三维X射线照相定位[7, 9-10]等技术的发展，肩部骨骼的姿态测量成为了可能，同时肩部机构模型也朝着更贴近于实际解剖结构的并联机构方向发展。Garner和Pandy[11]于2001年以美国国家医学图书馆的核磁共振影像资料数据库为依据建立了一个具有13自由度的人体上肢骨骼模型，其中胸锁关节(Sternoclavicular joint, SC)、肩锁关节(Acromioclavicular joint, AC)以及盂肱关节(Glenohumeral joint, GH)皆被视为球形铰链，并首次将胸廓建模为椭球，同时用肩胛骨在椭球表面的相对滑动描述肩胛胸廓关节(Scapulathoracic articulation, ST)的运动。该模型虽未从机构学角度精确定义ST几何约束，也未对该模型进行运动分析，但是却提供了一套较为完整的人体结构数据可供他人使用。基于此模型的概念，Maurel等[12-13]进一步将ST视作单个点与面接触而形成的约束，但是该研究同样没有对机构的整体进行分析，而是将闭式链拆开，以ST处的虚拟接触点的位置为输入，以此来模拟和分析肩胛骨的运动。
综上所述，在建立基于肩部骨骼解剖结构的骨骼机构模型时，肩带中的3个关节(SC、AC、ST)以及3块骨骼(锁骨、肩胛骨、胸骨)的运动必须予以充分考虑。因此，本文基于肩部骨骼的解剖结构，在Garner和Pandy[11]、Maurel等[12-13]的基础上提出了一个基于空间混联机构的描述肩带内部骨骼运动约束关系的运动模型，并建立了其位置闭环方程，得到了位置分析的闭合解，并最终通过运动实验数据验证了该模型的正确性。
1 骨骼系统机构学模型的建立人体肩部可认为是由肩带(shoulder girdle)和肱骨(humerus)2部分构成的骨骼系统，其中肩带由胸廓(thorax)、锁骨(clavicle)、肩胛骨(scapula) 3块骨骼组成，可视为一个闭环骨骼系统。与肩部运动相关的关节包括：胸锁关节、肩锁关节、肩胛胸廓关节和盂肱关节，如图 1所示。
|
| 图 1 肩关节机构模型简图 Fig. 1 Schematic of shoulder joint mechanism model |
本文使用一空间混联机构模型来模拟肩部骨骼系统，相比于传统的骨骼模型能够更真实地的反映肩部复合体的内部骨骼运动情况，如图 1所示。其中，胸锁关节a、肩锁关节b、盂肱关节e皆被视为理想3自由度球形铰链关节。与Maurel模型[12-13]和Tondu模型[14]中关于ST的定义不同的是，本文将胸廓视为一个椭球体，而肩胛骨通过其肩胛冈三角c(Trigonum spinae Scapulae, TS)与肩胛骨下角d(Angulus Inferior, AI)与胸廓椭球表面实现两点接触，其可等效于一个圆柱-平面副，具有4个自由度。其中，构件1表示锁骨，构件2表示肩胛骨，构件3表示胸廓表面，构件4表示上肢肱骨。由于本文的研究对象仅为肩部，因此假设上肢处于伸直状态，且上臂和前臂之间无相对运动，以肱骨长轴代表构件4，即GH中心e至肘关节中心f 的连线。其中f定义为外上踝(Lateral Epicondyle, EL)和内上踝(Medial Epicondyle, EM)连线的中点。关节a与b中心距离设为l1，b与c距离为l2，b与d的距离为l3，肱骨长度为l4。
2 混联机构模型的自由度分析要分析混联机构模型混联机构的运动特征，首先要分析其自由度。由Kutzbach-Grübler公式可得机构的自由度为
|
(1) |
式中：n为机构构件数目；g为运动副数目；fi为第i个运动副的相对自由度数。模型共有4个构件，4个运动副，其中SC、AC和GH为球副，自由度皆为3，ST为圆柱-平面副，自由度为4，因此肩部机构总的自由度为7。然而在肩部运动过程中，锁骨绕自身轴线方向的转动非常微小，可忽略不计[15-17]。同时从机构角度可知，对于关节a而言，其绕锁骨轴线方向的转动为虚自由度。因此整个肩部机构自由度为6，其中对应于肩带的并联部分自由度为3。
3 混联机构的运动学 3.1 坐标系的建立为了对机构进行位置分析，定义附着于单个骨骼的坐标系如图 2所示。假定胸廓固定不动，定义附着于其上的坐标系S0{a-x0y0z0}为固定系，其原点位于a点，x0、y0、z0轴分别平行于人体的解剖学主轴——冠状轴、矢状轴、垂直轴。S1{a-x1y1z1}为锁骨系，其原点位于胸锁关节a处，z1轴沿锁骨轴线方向，由a指向b，由于锁骨只有2个骨骼标记点，因此定义其x1轴位于水平面且垂直于z1，y1则由右手定则确定。S2{b-x2y2z2}为肩胛骨系，其原点位于肩锁关节b点，其z2轴为点b指向点c，x2垂直于b、c、d 3点所确定的肩胛骨平面并指向前方，且垂直于z2轴，y2轴则由右手定则确定。S3{O-x3y3z3}为胸廓系，其原点位于胸廓椭球中心，其坐标轴分别与椭球的3个主轴重合，且同时平行于人体解剖主轴，即S3//S0。S4{e-x4y4z4}为肱骨坐标系，其原点位于e点处，z4沿肱骨长轴方向且由f指向e，y4轴垂直于点e、EL、EM所确定的平面且指向内侧，x4由右手定则确定。
|
| 图 2 肩部骨骼机构坐标系定义 Fig. 2 Definition of related reference frames of shoulder skeleton mechanism |
在肩带中，胸廓、锁骨、肩胛骨互相之间通过关节SC、AC、ST相互约束，从而构成了肩部混联机构的并联部分。由于肩带机构的自由度为3，因此只需任意3个关节变量作为输入，就可确定整个肩带的姿态。本文进行运动学分析的目的在于以解析解的方式建立肩带各骨骼的运动约束关系，并验证所提出的空间混联机构模型对实际肩部运动预测的准确性。
从系S0{a-x0y0z0}到S1{a-x1y1z1}，可认为是先绕z轴旋转θ1，再绕x轴旋转θ2，再绕z轴旋转θ3而得到。从系S1到S2，可认为是先沿z轴平移l1，再绕y轴旋转θ4，再绕x轴旋转θ5，最后绕y轴旋转θ6而得到。因此有变换矩阵
|
(2) |
|
(3) |
假设椭球中心点O在固定系S0中的位置为0pO=[xO, 0 yO, 0 zO, 0 1]T，则椭球系S3相对于固定系S0的变换矩阵为
|
(4) |
设胸廓椭球曲面的半轴分别为c1、c2、c3，而点c和点d位于胸廓椭球表面上，其在系S3的位置可由椭球的参数方程表示如下:
|
(5) |
|
(6) |
式中:φc、ϕc、φd、ϕd为分别描述点c与点d位置的参数。如图 3所示，c、d同时也位于肩胛骨上。则c、d在系S2中的位置分别为
|
(7) |
|
(8) |
|
| 图 3 点c、d在肩胛骨系中的位置 Fig. 3 Position of points c and d in scapula frame |
式中：γ为
|
(9) |
|
(10) |
式(9)与式(10)即为肩带机构的关节位置求解约束方程，令其各项分量分别为0可得到6个约束方程，其中含有θ1、θ2、θ3、θ4、θ5、θ6、φc、ϕc、φd、ϕd共10个关节变量。而锁骨绕自身轴线方向的转动为虚自由度，若令θ3=0，则只需以任意3个关节变量作为输入就可求解其余6个关节变量，本文以θ1、θ2和ϕd作为输入变量。
由锁骨系的定义可知1pb=[0 0 l11]T，若已知θ1、θ2，则对于点b有
|
(11) |
图 4为肩胛骨与胸廓椭圆接触示意图。为了简化计算，首先分析肩胛骨上的b、c、d 3点，由于△bcd所代表的肩胛骨被认为是刚体，因此有
|
(12) |
|
(13) |
|
(14) |
|
| 图 4 肩胛-胸廓椭圆约束示意图 Fig. 4 Schematic of the constraint between scapula and thorax ellipsoid |
将式(5)、式(6)、式(11)代入式(12) ~式(14)即可得到3个含有变量ϕd、φd、φc、ϕc的约束方程。若ϕd为已知，则可求解得到φd、φc、ϕc。再由式(5)和式(6)可计算得到3pc和3pd的值。同时可得
|
(15) |
|
(16) |
再将0pc和0pd 代入回路方程式(9)和式(10)，则可推导得到中间变量A、B、C、D如下:
|
(17) |
|
(18) |
|
(19) |
|
(20) |
由此可得到θ4、θ5、θ6的封闭解如下：
|
(21) |
|
(22) |
|
(23) |
本文以Bolsterlee和Veeger[15]于2014年发表的上肢运动实验数据库(Shoulder Database V1.1)作为机构运动分析的数据来源。该运动数据库记录了人体上肢在完成一定的动作任务时，相关人体骨骼标记点在实验室坐标系中的空间位置坐标。所有相关数据可以从美国国家康复仿真研究中心(NCSRR)的网站中获得[18]。
该实验利用一套运动捕捉系统(optotrak system)检测分别固定于胸廓、肩胛骨、肱骨、前臂和手掌上的5个标记群的空间位置，如图 5所示[12]。所有骨骼标记点的位置按照Wu等[19]提出的ISB (International Society of Biomechanics)标准布置，如图 6所示。图中骨骼标记点已经由ISB标准详细定义[19]。
图 6中坐标系St、Sc、Ss、Sh为ISB标准中的各骨骼的局部坐标系，分别对应于计算系中的S0、S1、S2、S4。标记点的位置数据采集频率为100Hz。本文从该数据库中选取一个男性对象(年龄29岁，身高186cm，体重109kg)的运动数据作为参考。
|
| 图 6 ISB标准中的骨骼标记点与局部坐标系 Fig. 6 Skeletal landmarks and local frames defined in ISB standard |
在实验过程中，受测对象被要求完成运动范围任务，包括上肢前屈(FLEX)、外展(ABD)和肩胛平面运动(SCAP)。
4.2 机构尺寸的确定Garner和Pandy[11, 20]依据获得自美国国家医学图书馆的人体影像数据建立了具有较高尺寸精度的人体上肢的骨骼模型。本文据此确定肩部机构模型中的各杆件尺寸，如表 1和表 2所示。表 3为Garner标准模型与实验对象的骨骼特征尺寸在系St中的坐标值的对比，实验对象缩放系数为(1.20, 1.05, 1.14)，平均系数为1.13。
| 参数 | l1 /mm | l2 /mm | l3 /mm | l22 /mm | l4 /mm | γ /(°) |
| 数值 | 182.5 | 129.1 | 185.9 | 117.6 | 301.6 | 38.8 |
| 几何尺寸 | c1 | c2 | c3 | O点位置0pO | e点位置2pe |
| 数值 | 144.6 | 95.6 | 211.7 | [0-62.1 -152.1 1]T |
[28.1-43.9 -164.5 1]T |
由于Garner骨骼模型具有较高的尺寸精度且包含胸廓椭球的尺寸，因此将其作为标准模型，而实验对象的骨骼尺寸则通过对标准模型的缩放而得到。相应的缩放系数可通过对比标准模型与实验对象各骨骼标记点的位置而确定[15]。
实验对象的机构模型中各杆件尺寸可由标准模型的相应尺寸乘以平均缩放系数而确定。胸廓椭球的位置及3个半轴的长度则由标准模型的椭球尺寸乘以相应方向的缩放系数而确定。
4.3 ISB系与计算系的转换关系运动实验所采集的标记点位置及各骨骼的姿态是在Wu等[19]推荐的ISB标准系中记录和描述的，而ISB标准同时也是生物医学工程领域广泛采用的标准，本文亦将骨骼姿态的模型计算结果和实验测量结果统一放在ISB标准中来对比衡量。因此必须建立本文计算坐标系与ISB坐标系之间的转换与对应关系。
依据各局部坐标系的定义，以及表 3中各标记点的坐标，3.1节中用于计算的坐标系S0、S1、S2、S4与ISB标准坐标系St、Sc、Ss、Sh的关系可分别表示为如下旋转变换矩阵：
|
(24) |
|
(25) |
|
(26) |
|
(27) |
| 骨骼标记点 | 标准模型 | 实验对象 |
| IJ | (0, 0, 0) | (0, 0, 0) |
| PX | (31.9, -132.7, -9.8) | (56.0, -163.0, -2.4) |
| C7 | (-124.2, 54.1, 0) | (-123.1, 80.7, -4.7) |
| T8 | (-156.61, -171.5, 0) | (-179.1, -133.2, 11.2) |
| AA | (-105.6, 7.5, 182.6) | (-73.5, 9.3, 205.1) |
| TS | (-156.0, -11.7, 75.0) | (-136.3, 5.1, 121.8) |
| AI | (-156.7, -126.2, 101.9) | (-141.6, -126.5, 116.2) |
| AC | (-71.8, 26.6, 165.1) | (-37.0, 40.4, 167.2) |
| SC | (-2.8, -15.2, 1.4) | (15.3, -31.0, -6.6) |
| 椭球中心 | (-62.1, -152.1, 0) | (-61.6, -148.1, 0) |
| 椭球半轴 | (95.6, 211.7, 144.6) | (114.4, 222.7, 164.8) |
5 基于混联机构的运动模型验证 5.1 模型驱动参数的确定
通过对实验中测得的各标记点空间坐标的处理，得到关节变量θ1、θ2、ϕd，以此3个独立关节变量为输入，则可确定肩带的运动姿态。
以FLEX运动为例，图 7为测得的实验对象的锁骨姿态，并以YXZ序列下的3个欧拉角γc、βc、αc来描述[19]。
|
| 图 7 锁骨的姿态(FLEX) Fig. 7 Posture of clavicle (FLEX) |
由此可得在任意时刻，锁骨系Sc相对于实验固定系St的姿态矩阵为
|
(28) |
将式(24)和式(25)代入式(28)，可得
|
(29) |
令式(29)的矩阵的第3列元素与式(2)的矩阵的第3列前3个元素对应相等，则可得
|
(30) |
|
(31) |
由此获得了输入关节变量θ1、θ2，如图 8所示。
|
| 图 8 用于机构位置分析的输入关节变量θ1、θ2 (FLEX) Fig. 8 Joint variables θ1, θ2 inputted for mechanism position analysis (FLEX) |
另一个输入关节变量ϕd则由AI点相对于胸廓椭球系的空间位置而得到。在FLEX运动中，AI在椭球系S3中的位置3pAI(XAI, YAI, ZAI)如图 9所示。
|
| 图 9 肩胛骨下角AI在椭球系S3中的位置(FLEX) Fig. 9 Position of angulus inferior AI described in ellipsoid frame S3(FLEX) |
虽然AI非常接近于椭球表面，然而AI点始终位于肩胛骨上，与胸廓椭球表面之间存在一定间隙，与机构模型中所构想的点d并不完全一致。本文通过将肩胛下角AI沿椭球中心线方向的投影点视为等效的d点。如图 10所示，AI与椭球中心O的连线交椭球表面于d。
|
| 图 10 AI投影于胸廓椭圆表面示意图 Fig. 10 Schematic of point AI projected on thorax elliptical surface |
由此得到了3pd=[xd, 3 yd, 3 zd, 31]T，再依据式(6)，有
|
(32) |
由此则可得到ϕd的值，如图 11所示。
|
| 图 11 关节变量ϕd(FLEX) Fig. 11 Joint variable ϕd (FLEX) |
至此，完整地得到了FLEX运动相应的3个输入关节变量θ1、θ2、ϕd。同样地可以对ABD和SCAP运动进行分析，在此不再赘述。
5.2 肩胛骨姿态验证确定了输入关节变量θ1、θ2、ϕd后，将其代入式(12)~式(23)中，则可求解得到其余关节变量φd、φc、ϕc与θ4、θ5、θ6的值，进而确定整个肩带的运动。将各关节变量代入式(1)和式(2)，则可得
|
(33) |
由此可得系S2相对于系S0的姿态矩阵20R，再依据式(24)和式(26)可得
|
(34) |
同时，在肩带机构中假设以肩胛骨为末端执行件。在YXZ序列下其姿态可用3个欧拉角γs、βs、αs用以描述[19]。则肩胛骨Ss相对于固定系St的姿态可表示为
|
(35) |
因此，将姿态矩阵stR按照式(35)的YXZ旋转序列进行分解，可得到γs、βs、αs 3个关节角度值。
图 12为3组运动相应的肩胛骨姿态的模型预测结果(γmod、βmod、αmod)与实验测量结果(γmeas、βmeas、αmeas)的对比，从其中可以看出利用本文所提出的机构模型对于肩带运动姿态的描述与实验测量结果吻合较好，其基本趋势一致。
对于每一个运动任务，每一运动时刻，肩胛骨姿态的模型预测值与实验测量值之间的平均偏差为
|
(36) |
|
| 图 12 肩胛骨相对于胸廓的运动姿态(FLEX、ABD、SCAP) Fig. 12 Movement posture of scapula with respect to thorax (FLEX, ABD, SCAP) |
R值越小，则肩胛骨姿态的模型预测值与实验测量值的相近程度越好。平均偏差的计算结果如图 13所示，R为R的平均值。
|
| 图 13 肩胛骨姿态的模型预测值与实验测量值平均偏差 Fig. 13 Average deviations between scapula postures obtained from model prediction and experiment measurement |
从图 13中可以看出，对于3组运动而言，模型预测的肩胛骨姿态欧拉角与实际测量值之间的偏差较小。最小的R存在于SCAP运动中，为3.56°，最大的R存在于ABD运动中，仅为5.71°。
6 结论1) 本文通过将肩胛骨与胸廓的相对运动关系定义为类似于圆柱-平面副的线面接触运动约束，从而将包含胸骨、锁骨、肩胛骨和肱骨在内的肩部骨骼系统抽象成一个4杆6自由度的混联空间机构模型，并完成了肩带机构的位置建模和封闭解求解。
2) 通过骨骼尺寸的缩放，建立了基于个体特征的肩部骨骼机构学模型，为个性化肩部模型的建立和运动分析提供了基础。
3) 通过获得自肩部运动实验的数据，验证了所提出的混联机构模型对于人体肩部各骨骼之间运动约束关系描述的正确性。
| [1] | YANG J, FENG X, KIM J H, et al. Review of biomechanical models for human shoulder complex[J]. International Journal of Human Factors Modelling & Simulation, 2010, 1 (3): 271–293. |
| [2] | YANG J J, FENG X, XIANG Y, et al. Determining the three-dimensional relation between the skeletal elements of the human shoulder complex[J]. Journal of Biomechanics, 2009, 42 (11): 1762–1767. DOI:10.1016/j.jbiomech.2009.04.048 |
| [3] | DEMPSTER W T. Mechanisms of shoulder movement[J]. Archives of Physical Medicine & Rehabilitation, 1965, 46 : 49–70. |
| [4] | ENGIN A E, CHEN S M. Statistical data base for the biomechanical properties of the human shoulder complex-Ⅰ:Kinematics of the shoulder complex[J]. Journal of Biomechanical Engineering, 1986, 108 (3): 215–221. DOI:10.1115/1.3138605 |
| [5] | ENGIN A E, TUMER S T. Three-dimensional kinematic modelling of the human shoulder complex-Part Ⅰ:Physical model and determination of joint sinus cones[J]. Journal of Biomechanical Engineering, 1989, 111 (2): 107–112. DOI:10.1115/1.3168351 |
| [6] | XU X, DICKERSON C R, LIN J H, et al. Evaluation of regression-based 3-D shoulder rhythms[J]. Journal of Electromyography and Kinesiology, 2015, 17 : 28–33. |
| [7] | MASO F D, BLACHE Y, RAISON M, et al. Glenohumeral joint kinematics measured by intracortical pins, reflective markers, and computed tomography:A novel technique to assess acromiohumeral distance[J]. Journal of Electromyography and Kinesiology, 2015, 29 : 4–11. |
| [8] | HARISH V, BAKSH A, UNGI T, et al.Measurement of electromagnetic tracking error in a navigated breast surgery setup[C]//Proceedings of the SPIE Medical Imaging.Bellingham:SPIE Press, 2016, 9786:1-8. |
| [9] | ZHANG C, SKALLI W, LAGAC P Y, et al. Investigation of 3D glenohumeral displacements from 3D reconstruction using biplane X-ray images:Accuracy and reproducibility of the technique and preliminary analysis in rotator cuff tear patients[J]. Landscape Design Journal of the Landscape Institute, 2015, 29 (10): 12–20. |
| [10] | DAL M F, RAISON M, LUNDBERG A, et al. Glenohumeral translations during range-of-motion movements, activities of daily living, and sports activities in healthy participants[J]. Clinical Biomechanics, 2015, 30 (9): 1002–1007. DOI:10.1016/j.clinbiomech.2015.06.016 |
| [11] | GARNER B A, PANDY M G. Musculoskeletal model of the upper limb based on the visible human male dataset[J]. Computer Methods in Biomechanics & Biomedical Engineering, 2001, 4 (2): 93–126. |
| [12] | MAUREL W, THALMANN D. Human sholder modeling including scapulo-thoracic constraint and joint sinus cones[J]. Computers & Graphics, 2000, 24 (2): 203–218. |
| [13] | MAUREL W, THALMANN D, HOFFMEYER P, et al.A biomechanical musculoskeletal model of human upper limb for dynamic simulation[C]//Proceedings of the Eurographics Workshop on Computer Animation and Simulation.Vienna:Springer, 1996:121-136. |
| [14] | TONDU B.Modelling of the shoulder complex and application the design of upper extremities for humanoid robots[C]//Proceedings of the IEEE-Ras International Conference on Humanoid Robots.Piscataway, NJ:IEEE Press, 2005:313-320. |
| [15] | BOLSTERLEE B, VEEGER H E J. Modelling clavicular and scapular kinematics:From measurement to simulation[J]. Medical & Biological Engineering & Computing, 2014, 52 (3): 283–291. |
| [16] | SHOLUKHA V, JAN S V S. Combined motions of the shoulder joint complex for model-based simulation:Modeling of the shoulder rhythm(ShRm)[M]. London: Springer, 2014: 205-232. |
| [17] | WARNER M B, CHAPPELL P H, STOKES M J. Measurement of dynamic scapular kinematics using an acromion marker cluster to minimize skin movement artifact[J]. Journal of Visualized Experiments Jove, 2015, 96 . |
| [18] | The National Center for Simulation in Rehabilitation Research(NCSRR).Shoulder movements database:The supplement of Delft Shoulder and Elbow Model Project[DS/OL]. (2015-03-20)[2016-05-30]. https://simtk.org/projects/dsem. |
| [19] | WU G, HELM F C T V D, VEEGER H E J, et al. ISB recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion-Part Ⅱ:Shoulder, elbow, wrist and hand[J]. Journal of Biomechanics, 2005, 38 (5): 981–992. DOI:10.1016/j.jbiomech.2004.05.042 |
| [20] | GARNER B A, PANDY M G. A kinematic model of the upper limb based on the visible human project (VHP) image dataset[J]. Computer Methods in Biomechanics & Biomedical Engineering, 1999, 2 (2): 107–124. |
