基于红外热图边界分割方法的康复机器人抗阻训练肌肉激活状态评估①

张珝， 李益斌， 徐东， 徐晗， 孙亚， 李晓龙， 傅建明

1.同济大学浙江学院 奥克兰·同济康复设备研究中心，浙江 嘉兴 314051；2.嘉兴学院 附属第二医院，浙江 嘉兴 314001

运动功能障碍是神经系统疾病或肢体损伤后的常见后遗症，在这当中，又以上肢功能障碍最为常见，严重影响了患者的日常生活[1-3]． 利用康复机器人进行训练是上肢运动功能障碍患者进行肌力训练的重要方法之一[4]，康复机器人的抗阻训练模式能够激发患者患肢的自主收缩力，是患者恢复期增强肌力的重要措施[5]． 因此，合适的评估策略能够有效地量化康复抗阻训练的效率． 临床康复评估多以评估量表为主，如Brunnstrom量表、 Fugl-Meyer量表等，但是，这些方法均为康复医师主观评定方法，且评定结果无法定量分析；康复机器人可以通过搭载各式传感器对使用者肢体的位移、 (加)速度等运动学参数进行定量测量[6-8]，但是，不同的康复机器人会对不同的运动学参数进行测量，无法形成统一的标准；等速肌力测试系统、 高速摄像机和三维测力系统等设备能够基于统一标准对肌力、 平衡、 关节运动角度等运动学参数进行精确测量[9-10]，但是，这些设备因过于庞大而不便于使用，且造价昂贵；表面肌电信号(sEMG)作为人体的生理信号，能够有效地评估不同肌肉对于康复训练过程所造成的影响以及肌肉疲劳程度等[11-12]，然而，表面肌电测量需要在目标肌肉上放置有线表面电极，这会阻碍受试者在训练期间的运动，且电极之间的距离以及电极的大小也会影响测量的结果．

红外热成像技术具有无损、 非接触和非侵入性等特点，近年来被广泛地应用到临床医学当中，通常被用作疾病治疗前后热效应的评估以及部分疾病的诊断[13-16]． 当人体处于疾病或者运动状态时，皮肤表面温度就会出现差异，早些时候的一项研究也表明，长时间的运动会增加肌肉热量的扩散[17]，因此，不少研究都在探寻运动期间肌肉的温度变化关系． González-Alonso等[18]利用热敏电阻来研究剧烈运动后肌肉局部温度的变化，但是，该方法是侵入式测量，热敏电阻需要插入受试者的静脉，对受试者造成了一定的伤害；Zagrodny等[19]通过对步态分析下肌肉的激活表明了sEMG信号和皮肤表面热特性之间的相关性；Daud等[20]通过对上肢不同运动的分析也表明了sEMG信号与热成像测量的皮肤温度之间的相关性． 上述的研究证明，热成像测量的皮肤表面温度变化与肌肉激活状态高度相关，抗阻训练会造成血液流速的加快以及肌肉的收缩从而产生大量的散热，并以此确认了红外热像技术作为评估手段的可能性，而康复抗阻训练效果的评估对于判断患者肌力的恢复以及后续的康复策略制定至关重要，因此，有必要研究出一种基于康复抗阻训练下肌肉激活状态的评估方法．

本文旨在研究红外热成像技术在上肢协调康复机器人抗阻训练模式下肌肉激活状态评估的应用． 在不同的训练条件下，获取不同训练阶段下的红外热图，并对图像进行相应的处理和分析，实现对目标肌肉区域激活状态的可视化，同时实现区域内激活状态定量分析． 这项研究的优点是，通过利用以非侵入性和非接触方式获取的热图像，可以有效地评估康复训练中肌肉的激活状态．

1 实验设备与实验步骤

1.1 红外热像设备

本研究所采用的红外热像设备是基于美国FLIR公司Boson 320长波红外热像仪机芯研发而成(图1)，该机芯搭载了9.1 mm焦距的光学镜头，分辨率为640×512，通过USB接口与电脑连接，用于对红外热像图(NIR)的捕捉(图1)．

图1 红外热像设备

1.2 上肢协调康复机器人

上肢协调康复机器人是由同济大学浙江学院奥克兰·同济康复医疗设备研究中心自主研发、 具有完全自主知识产权的新一代康复医疗设备． 该康复机器人主要适用于脑卒中、 外伤等引起的上肢功能障碍患者，且具备多级阻力变化的抗阻训练模式[21-22]，能够有效地针对不同康复阶段的患者进行抗阻肌力训练(图2)．

图2 上肢协调康复机器人

1.3 实验步骤

在设计实验步骤时，考虑到实验中部分阻力等级较大，不易驱动，因此，实验选取了5名健康男性参与者参与实验，他们的年龄为(30±4)岁，身高为(175.20±5.36)cm，体重为(74.80±7.95)kg，每位参与者都签署了《实验知情同意书》． 参与者上肢部位肌肉或骨头没有遭受过结构性损伤，且24 h内未进行过激烈的体力运动，保证上肢部位状态的正常；同时，为了尽量保证红外热像测量的精确性，实验过程中实验室温度和湿度始终保持在(21±1)℃和(45±5)%． 图3显示了本次实验的流程： 实验主要基于上肢协调康复机器人(本文简称机器人)完成，参与者将以右上臂作为参考肌肉群，在接受红外热图采集前，将手臂暴露在空气中10 min，用于摆脱因为衣物遮挡导致的温度差异；每位参与者坐在离红外热像设备约1 m的位置，根据不同参与者身高的差异调整座椅高度，手臂自然下垂以保证被测部位完整地出现在红外热像镜头中央；为了最大限度地减小环境带来的影响，参与者手臂需要穿过带有方形孔的纸板，以此屏蔽除被测部位外的其他肢体部位． 实验开始时，采集参与者被测部位自然状态下的红外热图作为基础参照图像(Bs)，然后选择合适的位置坐于机器人前，机器人训练平台向上倾角为30°． 实验期间，参与者将基于机器人抗阻训练模式下3种不同阻力等级(L1，L5，L8，阻力等级由小至大)分别进行推、 拉两种肌力训练各5 min；单次训练结束后，参与者拥有5 min的恢复时间，考虑L8级阻力较大，在L8级训练结束后，参与者拥有10 min的休息时间以尽可能地保证目标区域状态的完全恢复；最后，分别采集参与者右上臂3种状态下： 训练刚完成(0 min)、 训练完成后1 min、 训练完成后3 min的红外热图，并进行对比分析．

图3 实验步骤

2 实验分析方法

2.1 基于归一化互相关算法的图片重新定位

重新定位在图像处理中非常重要． 在很多情况下，人们对原始图像中感兴趣区域(region of interest, ROI)的结果更感兴趣． 但是，手动选择相同的ROI较为麻烦且很难定位准确，因此，本研究采用了一种基于归一化互相关算法的重新定位技术． 当第一幅图像确认选择的ROI后，重新定位技术将自动在后续图像中找到与ROI匹配度最高的位置．

如果用矩阵R表示ROI区域，用矩阵T表示目标图像，首先对矩阵R进行归一化：

(1)

式中：Rmax和Rmin是矩阵R的最大值和最小值，RN是矩阵R归一化后形成的新的矩阵，同样，矩阵TN是矩阵T归一化后得到的． 矩阵RN和矩阵TN的互相关性为：

(2)

式中：i和j是矩阵的大小，RM是RN的平均值，TM是TN的平均值． 结果Corr可以帮助找出两个矩阵的最佳相关性，这将是原始图像R中ROI区域在目标图像T中的位置(图4)．

红色方框为基于原始图像(original)选取ROI后续图像的定位

2.2 灰度图像彩色化处理

红外热图的灰度图像没法很好地看出图像间的变化，为了使图像可视化更加清晰，需要对灰度图像彩色化处理． 彩色化处理是将灰度值矩阵中的每一个值拆分成3色(red green bule，RGB)通道，同时定义一个映射函数，用来对灰度图像值进行区分：

(3)

式中，R(x，y)代表RGB的任意一个通道矩阵，X(x，y)代表灰度矩阵于，Rc1和Rc2表示自定义RGB单一通道矩阵下颜色强弱区分的阈值，Gc1和Gc2表示自定义灰度矩阵区分的阈值，灰度矩阵根据(3)式进行变换，结合另外两个通道形成新的RGB矩阵，完成灰度图像的彩色化处理．

2.3 基于边界轮廓的图像分割方法

在本实验中，不同肌肉区域内目标位置基于灰度值变化的对比度非常重要． 为了突出这项变化，本研究采用了一种基于边界轮廓的图像分割方法，对图像进行一种模糊增强，从而达到增强灰度值对比度的目的．

首先，利用图像二值化对目标图像进行处理：

(4)

式中，G(x，y)表示目标图像灰度矩阵，Tg表示分割阈值． 将小于Tg的像素变为黑色(0)，将大于Tg的像素变为白色(255)，完成图像的黑白转化，从而突出图像的边界，考虑图像不同位置Tg不同，本研究将目标灰度图像切成4份，通过不同的Tg来对图像边界精确划分，最后通过将图像合并回原始图像：

(5)

式中，Gi为G切割中的一部分．

黑白图像边缘轮廓可视为一条二维封闭的边缘曲线L，且可以表示为有限个点的集合P={P1，P2，…，Pk}，P∈[0，L]，集合的坐标表达式为：

(6)

式中：pi为边界Pi(pi，qi)的x坐标；qi为该点的y坐标；

接下来，先将图像在纵向进行分割，之后以纵向分割为基础对图像进行横向分割，以此得到新的包含边界的网格交点坐标：

P′=(p′s，q′s，v)

(7)

式中：

(8)

式中：s为新分割图像横坐标，s∈[0，m]，m为纵向网格数量，v为基于s的纵坐标，v∈[0，n]，n为横向网格数量．

由此，目标图像被分割成新[m×n]的图像，然后通过计算每个网格的平均灰度值来生成新的图像．

2.4 目标图像处理

原始目标图像如图5a所示，利用2.1的重新定位方法修剪右上臂ROI区域(图5b)． 应用移动平均滤波器后，采用大津的全局阈值法[23]突出ROI区域(图5c)． 然后将ROI区域利用2.3所述的图像分割方法将图片平均分为20个水平部分和45个垂直部分(图5d)． 通过计算ROI区域中每个分割的网格平均灰度值生成模拟热图(图5e)，并使用彩色图显示(图5f)． 为了定量评估肌肉不同条件下的激活状态，将以图5e为基准提取相对灰度值进行分析．

图5 图像处理过程

2.5 相对灰度值

在红外热像技术中，红外热图的灰度值与目标温度成正相关线性关系[24]，被测目标温度的高低可以使用红外热图的灰度值大小来表示，同样，基于灰度值变换而来的彩色热图也能通过不同颜色的区分来表示温度的变化． 本研究采用的红外热像设备能够很好地显示被测目标温度的高低，但是，由于该设备并未进行过校准，被测目标的灰度值会随着测试范围内温度上限的变化而发生变化(目标温度并未发生变化)，灰度值能够表示目标温度的变化却不能表示绝对温度． 为了解决这个问题，本研究采用了一种相对灰度值的表示方法(图6)，即提取目标两个相邻位置的灰度值之差作为该位置对比度变化的相对灰度值：

rGS=GS1-GS2

(9)

式中：GS1表示相邻位置区域1，GS2表示相邻位置区域2．

为了有效地评估肌肉激活状态，本实验基于上臂运动时主要参与的肌肉并结合生理解剖学将目标部位划分为3块区域，即三角肌区域、 肱二头肌区域和肱三头肌区域(图6a)，每个区域提取一个相对灰度值(rGS)作为肌肉激活状态参数(图6b中1，2，3)，并利用2.1所述方法重新定位其他图像的相同区域(图6c)． 为了更好地识别相对灰度值，本研究选取了以各个训练状态下训练刚完成时(0 min)的图像作为目标图像进行定位．

图6 相对灰度值选取

3 实验结果与讨论

3.1 实验结果

表1显示了5名参与者在进行推、 拉力抗阻模式训练下3种阻力等级下3块肌肉区域平均相对灰度值在Bs—0 min—1 min—3 min时的变化．

表1 抗阻训练模式下参与者平均相对灰度值

可以看到的是，在推、 拉力抗阻训练中，所有训练刚结束时(0 min)与基础参照图像(Bs)相比，3个肌肉区域的相对灰度值都(rGS)有了较明显的提升，之后训练完成1 min和训练完成后3 min阶段，rGS逐渐下降．

为了了解抗阻模式训练下目标区域热图的变化，图7选取了其中一名参与者推力抗阻训练模式L8阻力等级下的完整训练状态彩色热图．

图7a-图7d表示目标肌肉区域Bs—0 min—1 min—3 min的彩色热图变化． 以2.5中的区域划分为参照，可以看到的是，0 min(图7b)彩色热图中，3个肌肉区域相对灰度值对比度较之Bs(图7a)相同区域有明显的提升，而在1 min((图7c))时，3个肌肉区域对比度有所降低，最后在3 min((图7d))时，三角肌区域和肱二头肌区域对比度继续下降，而肱三头肌区域对比度有所提升．

图7 一名参与者L8推力抗阻训练模式下彩色热图集

不同的阻力等级可能产生不同的热图变化． 图8展示了同一名参与者在推力和拉力抗阻训练模式下阻力等级L1和L8 3个肌肉区域在0 min时颜色对比度的变化．

图8 抗阻训练模式下阻力等级L1和L8彩色热图对比

图8a，b两图中均可以看出，阻力等级L8在3个肌肉区域的颜色对比度变化均比阻力等级L1更强；而从所有参与者在肱三头肌区域平均相对灰度值(rGS)的比较来看，在推力抗阻训练模式下，阻力等级L8在0 min的rGS变化比阻力等级L1更高(图9a)；然而在拉力抗阻训练模式下，L1在0 min的rGS变化比L8更高(图9b)．

图9 抗阻训练模式下阻力等级L1和L8的相对灰度值在肱三头肌区域的比较

为了对比推力和拉力抗阻训练模式下，3个肌肉区域彩色热图的变化，图10选取了同一参与者在推、 拉力抗阻训练模式下阻力等级L8 3个肌肉区域在0 min时颜色对比度的变化．

图10 推/拉力抗阻训练模式下阻力等级L8彩色热图对比

从图11中可以看出，推力和拉力抗阻训练模式下，3个肌肉区域的颜色对比度都有所增强，在推力抗阻训练中，肱三头肌区域对比度增强更为明显，拉力抗阻训练模式下，肱二头肌区域对比度增强更为明显；在从所有参与者3个肌肉区域平均相对灰度值(rGS)的比较来看，0 min下的rGS均有增加，推力抗阻训练模式下，肱三头肌区域rGS增加最大(图11a)，而拉力抗阻训练模式下，肱二头肌区域rGS增加最大(图11b)．

图11 推/拉力抗阻训练模式下的rGS在三肌肉区域的比较

3.2 讨论

肌力训练是康复运动疗法中的一种重要方法，是指在肌肉收缩时给予阻力可以使肌肉随时间产生适应而增加肌力[25]． 在Brunnstrom肌力分级理论中，在3级肌力到4级肌力进展时，应该进行阻力等级自小而大的抗阻训练． 本研究中的上肢协调康复机器人能够提供不同等级的抗阻训练，能够有效地对上肢肌肉进行肌力训练． 因此，本研究也将上臂三块主要肌肉所在的区域选为目标观察区域；同时，选取3个由小至大的阻力等级进行实验也符合肌力康复训练循序渐进的康复理论．

抗阻训练过程中，是训练时和训练后肌肉即时变化为负荷、 恢复的过程，在负荷阶段能量物质被大量消耗，代谢产物被蓄积起来，人体功能下降，产生疲劳；停止运动后，到了恢复阶段，机体内环境(热、 酸碱和水)开始恢复平衡，肌肉内的消耗得到了补充[26]． 如表1所示，所有肌肉区域在训练状态下0 min时的相对灰度值(rGS)均出现明显的提升，因为当抗阻训练刚结束时，整个上肢肌肉区域的肌肉产生热量打破了体内平衡，进入恢复状态后产生散热而导致温度的提升，表明了训练时3个区域的肌肉均处于激活状态；到1 min和3 min状态时，rGS逐渐下降，肌肉从激活状态恢复到正常状态，图7也很好地说明了这一点；但是，部分训练状态下(拉力抗阻训练L5和L8)，出现了rGS上升的情况，这可能是因为肌肉产生热量并缓慢地通过皮肤灌注引起激活肌肉附近血管扩张，这种热传递会导致皮肤表面温度变化的延迟，而L5和L8级的抗阻训练负荷更大，造成产热更大，因而导致在3 min这段监测时间内，热量一直持续散出．

相对灰度值(rGS)结合彩色红外热图能够更好地表现出肌肉的激活状态． 结合上面的理论，当阻力越大，肌肉运动产生的负荷越大，散热量也越大，肌肉的激活程度越高，图8在推、 拉力抗阻训练下阻力等级L8比L1下肌肉的激活程度更高，图9(a)肱三头肌区域rGS的变化也说明了这一点；但是在图9(b)中，出现了肱三头肌区域阻力等级L1的肌肉激活程度高于L8的情况，这可能是因为部分参与者在L8训练状态下选取rGS区域时，并未选取到对比度最大的位置，导致rGS在取平均值的过程中，出现了偏差；图10反映了在两种相反的抗阻训练(阻力等级L8)下3个肌肉区域的激活情况，从图10上能够看出，三角肌的激活状态相差不大，在推力抗阻训练时，肱三头肌区域激活程度更高，而在拉力抗阻训练时，肱二头肌区域激活程度更高，这是因为推力抗阻训练时，更加侧重于伸肌(肱三头肌)的发力，而拉力抗阻训练时，更加侧重于屈肌(肱二头肌)的发力，二者互为主动/拮抗肌，在做两组相反的运动时，产生了不同的激活效果，图11a和图11b也反映了相同的结果，这也是图9分别选择这两个区域作为对比的原因；同时还需注意到的是，由于上肢协调康复机器人调高了设备倾角，通常作为协同肌出现的三角肌在所有的训练状态下均被激活，在上肢协调康复机器人临床实验中，康复治疗师提出了调高设备倾角以达到提高患者上肢肌肉群肌力的康复效率，本实验结果也验证了这一理论．

在本次研究中，受试者的数量有限，且受试者的肌肉质量、 皮肤颜色和年龄等均会影响本研究的结果，这也是表1中部分数据标准差偏大的原因之一；为了实现测量的可重复性，有必要将相机和被摄体之间的距离以及其他变化降至最低，因此，可以尝试目标区域的运动跟踪；在实验中，虽然利用了基于相邻区域对比度最高的训练刚结束时(0 min)图像提取的相对灰度值(rGS)参数实现了参与者训练前后温度变化的量化分析，但是，该方法对相邻区域的选取要求较高，容易出现偏差，且不能更加直观地观测到目标区域温度的变化，因此，可以对本次实验使用的红外热像镜头进行校准，或更进一步完成温度的标定；本次实验全部选取了健康人作为参与目标，所有肌肉区域都处于正常可激活状态，这与上肢协调康复机器人实际应用中针对上肢功能障碍患者的情况有所偏离，未来工作是在改进方法的同时进一步验证该方法在上肢功能障碍者上的准确性和有效性．

4 结论

本研究通过红外热成像技术评估上肢3个肌肉区域在上肢协调康复机器人抗阻训练模式下的激活状态，并使用一种边界轮廓分割的图像处理方法将红外灰度图简化并彩色化处理，实现了红外热图可视化；同时，通过相对灰度值这一参数实现了对上肢肌肉激活状态的定量分析，之后有必要进行进一步的研究，将该方法改善后扩展到临床应用中去．