基于惯性传感器的上肢康复训练评估
神经肌肉系统疾病患者，往往需要进行长期的康复训练恢复自由活动的能力或保持肌肉活性。在肢体康复训练中，临床医生给患者开一系列康复训练的处方，无论是在患者面前演示还是移动患者肢体做动作，都需要沿着所需路径移动。经常是医生在患者面前重复几遍训练动作，患者通过观察加强了对训练动作的认知，包括动作的幅度、速度、角度或者其他体位姿势。然后根据指导方案，患者在家中持续训练。在家庭环境下进行康复训练，由于没有医生或理疗师在场指导，患者不能及时掌握自身的训练情况，医生也不能提供有效指导意见，在这种现实条件下可使用计算机系统辅助评估患者的康复训练。本研究提出一种基于惯性传感器的上肢康复训练评估框架，可自动评估上肢康复训练的质量，输出训练动作的分数。该框架的主要部分包括：训练动作的量化指标，将性能指标映射为动作质量分数的评分函数，通过监督学习回归动作质量分数的神经网络模型。框架结构图如下所示。
本研究选取《Effective Functional Progressions in Sport Rehabilitation》一书中推荐的一种训练患者上肢肌肉和关节的动作作为示例动作。此动作分为四个步骤，起始时将小臂平举于胸前，手掌向下；第二步向身体后方肘击，肩关节活动应达到最大角度；第三步手臂归为；最后将小臂打开，整条手臂向身体后方摆动，肩关节活动达到最大角度。
本研究采集一名健康人以及两名不同患病程度的神经肌肉系统疾病患者的上肢动作数据，用采集数据完成拟议框架。两名患者中，一人为轻度脑卒中患者，可自行完成训练动作；一人为小儿脑瘫患者，需要护士辅助完成动作。
4.1 运动测量系统设计与实现
（1）硬件设计
为完整的描述上肢训练动作，选择在右臂腕关节、肘关节以及肩关节处佩戴惯性传感器模块MPU6050，采集腕关节和肘关节动作产生的三轴加速度和角速度数据，肩关节仅采集三轴角速度数据。示例动作起始状态时对三个传感器执行初始化，腕关节处传感器模块水平安放，X轴指向手指，Y轴指向人体正前方，Z轴指向上方；肘关节和肩关节传感器坐标轴设定与腕关节一致。
如图所示，MPU6050模块使用IIC方式与单片机连接，其中XDA与XCL为次IIC接口，可接磁力传感器组成九轴传感器模块。AD0为选址引脚，AD0接高电平时，MPU6050地址为0X69，AD0接低电平时，MPU6050地址为0X68。本研究需要使用三个MPU6050模块，选择了三组通用IO口模拟三条IIC总线分别控制三个传感器模块采集数据。
图3.MPU6050模块原理图
系统硬件设计完整原理图及PCB附于附录。
（2）软件设计
系统软件实现了对三个传感器节点数据的采集与上传。MPU6050模块初始化流程如图所示。
图4. MPU6050初始化流程图
MPU6050模块中的加速度传感器的量程可设置为±2g、±4g、±8g、±16g，为取得更高精度，本次研究将量程设为±2g，寄存器中的读数换算成加速度公式为：
公式（1）中指的是加速度寄存器读数的转换结果，。陀螺仪的可选倍率有：250°/S、500°/S、1000°/S、2000°/S，为取得更高的测量范围，本次研究将陀螺仪倍率设置为2000°/S，换算公式为：
公式中指的是角速度寄存器读数的转换，。加速度计和陀螺仪的采样频率设置为50Hz。框图中设置其他参数指关闭中断、关闭AUX_IIC接口、禁止FIFO等步骤。
三个MPU6050模块依次初始化，间隔时间100us。初始化完成后读取相关寄存器得到加速度、角速度数据，读取的数据由JTAG发送给PC端，波特率设置为115200。
（3）数据展示
本次研究通过自主设计的惯性数据采集装置，可采集试验者上肢运动时手腕、肘部以及肩膀三处关节节点的加速度和角速度数据，现将部分采集数据作简单的介绍和分析。
图5.手腕三轴加速度数据
图5为30S内手腕处三轴加速度数据波形，其中y轴波动最大，x轴其次，当试验者右臂肘击到最大角度时，x轴和y轴加速度数据达到第一个峰值，右臂回收时，加速度值由正变负，当试验者右臂全伸展到最大角度时，x轴和y轴达到第二个峰值，因此时手腕沿y轴运动距离最大，所以y轴的第二个峰值明显大于第一个峰值。因为拟定动作应尽量保持水平执行，所以z轴波动最小，只有2g以内的波动，而重力加速度分量的存在，使得z轴数据最大。Z轴数据有小幅度范围内的波动，是因为上肢活动时人体的肌肉收缩或摆动等原因所致，无法避免，但通过去噪加窗数据等预处理手段，可将噪声信号干扰降至最低。在30S内，出现10次左右的
峰值，说明此项动作执行一次大约需要3秒的时间。
图6.三节点Z轴角速度数据度数据
图6为15S内，腕关节、肘关节和肩关节的Z轴角速度图示。由图可知，在一个动作周期内，手腕部和肘部的Z轴角速度数据曲线起落一致，几乎重叠，而肩关节在本示例动作中，只做往返的摆动，所以在一个动作周期内，不像手腕和肘部数据曲线的峰值有明显差异。
以上两图中，曲线情况基本符合预期，呈现出规律性，说明采集到的数据有效，可以进行下一步处理。
目前采集的数据有一组健康试验者所做的示例动作，以及两组患者所做的动作。两组患者动作中，一组由患者自主完成，一组由护理人员帮助患者完成动作。要求试验者重复做示例动作两分二十秒，去除首尾各十秒钟，保留的有效数据为两分钟。回传速率设置为50Hz，即每秒可采集到每一轴数据量为50个，每次采集到的样本组成一个6050×15的矩阵，本文中将采集到的健康试验者动作数据集合标记为，将采集到的患者动作数据集标记为和。
图7为健康试验者所做动作和模拟患者动作对比，显示的曲线为示例动作手腕部位y轴加速度数据。上边的图数据来自，下边来自。由于此患者无法自主完成动作，所以由护理人员帮助其完成训练动作，从数据曲线中可以直观的看出，在短时间内，患者的数据曲线显得杂乱，虽然也能从中区分出动作周期，但很明显，每一个动作周期的峰值、周期时间长短都不一致，这说明患者对自身肢体控制能力偏弱。由图可知，健康试验者的动作数据和患者的动作数据间存在着明显差别，这为评估患者训练质量提供了可能性，可把健康试验者动作数据作为参照，评价患者动作质量。
图7.正常动作数据与患者动作数据对比图
4.2 数据预处理
本次研究的数据预处理过程包括去噪、加窗、归一化。
（1）去噪
初步采集到的加速度和角速度数据中，由于设备的抖动、杜邦线连接的不稳定等原因，夹杂着许多干扰信号，从图中可以看到有许多毛刺，在处理这些数据之前，有必要先进行滤波处理。本文采用移动平均滤波的方法对数据滤波处理，移动平均滤波是一个低通滤波器，可以有效减少随机干扰的影响，适合用在本次研究当中。移动平均滤波基于统计规律，将连续的采样数据看成一个长度固定为N，在新的一次测量后，上述队列的首数据去掉，其余N-1个数据依次前移，并将新的采样数
图8.滤波前后对比
据插入，作为新队列的尾；然后对这个队列进行算术运算，并将其作为本次测量结果。图7为截取部分手腕部位X轴加速度信号的滤波前后对比，N设置为20，可以看出信号平滑了很多。
（2）加窗
肢体动作数据的采集往往需要较长的时间，一段时间内可能会采集到多组动作。而在动作评估中，需要评估的不是一段时间内的动作，而是每一个完整的训练动作。所以需要对数据流进行窗口分割，每个或者相邻几个窗口描述一个完整的动作。本研究采用固定长度的滑动窗口分割技术对传感器数据流进行划分，窗口的长度设置为2秒，每两个相邻窗口间有50%的重叠。对于每个时间窗口，它所包含的数据可以用一个N×1维的向量来表示。由上文中对数据曲线的分析可知，正常人做一个示例动作大约需要三秒钟，所以每三个邻近的特征值可描述一个动作。每个窗口提取的特征类型有三轴加速度信号的均值、方差、最大值、最小值，三轴角速度信号的均值，以及三轴旋转角的均值。其中由角速度信号计算得到，它能有效地反映手臂运动转向的大小。旋转角计算公式如下：
（3）归一化
在机器学习算法中，将特征数据归一化第一可加快计算速度，第二有可能提高精度。在本研究中，为了方便之后的建立数学模型及训练神经网络等步骤，也需要对特征数据进行归一化处理。采用下列公式将向量归一化到区间：
（4）
把采集到的三组数据预处理之后，得到三个120×42的矩阵。为了方便下文中高斯混合模型的建立，将三个矩阵合并，组成一个新的360×42的矩阵，记为。此矩阵囊括了健康者和两种患者的动作数据，下文将对此矩阵进行降维和建模操作。
4.3 降维
在本研究中，一共有三个传感器节点采集动作信号，其中腕关节和肘关节采集加速度和角速度信号，肩关节采集加速度信号，去噪加窗后提取这些数据的特征值，包括均值、最大（小）值、方差，所以获得的特征向量共有42维（18×2+2×3）。当提取的特征数量过多时，直接对高维特征进行分类会遇到很多问题，包括产生较大的计算量、数据不够直观以及数据可视性较差等。而在特征集合中，这些特征之间有些是相互关联的，并且在一定程度上存在着信息的重叠，比如图5所示，手腕部位和肘部的z轴角速度数据几乎重叠。
本研究中，采用主成分分析的方法对特征向量进行降维，去除数据冗余，以找到特征的最优线性组合，并排列出它们的重要性。主成分分析的主要思想是用一个新的低维数的特征集合来取代原始高维数据，让其尽可能的反映原来数据的信息。
对于新的特征集
中的每个主成分，都能够表示为原特征集
中所有特征的一个线性组合：
。其中a由原始特征空间的特征值分解来确定。每个主成分间相关系数为0。
采用主成分分析法进行特征降维的步骤如下：
1、计算训练样本集合的均值向量，并将所有样本中心化；
2、求出其协方差矩阵和m个特征值，并按照特征值的大小进行排列；
3、通过设置累计方差贡献率的阈值，选取前l个特征值所对应的特征向量，组成投影矩阵；
4、将原始m维特征向量投影到低维的向量空间中，得到降维后的特征集合。
本研究通过调用MATLAB中princomp函数将训练样本降至5维，测试样本降维时，对测试样本减去训练样本均值，再乘以系数矩阵矩阵，得到测试样本降维后的数据，保证了测试样本和训练样本转换到相同的样本空间中。
4.4 建立数学模型
为对患者训练动作进行评估，得到动作质量分数，首先要用动作数据建立数学模型。本次研究选择高斯混合模型（Gaussian Mixture Model,GMM）对示例动作建模。
高斯混合模型是一种用高斯概率密度函数来表达数据的参数概率模型，常用于人体运动建模。对于一个多维数据集，有个高斯分量的高斯混合模型通常为如下形式：
（5）
公式⑤表达的是一个概率密度函数，由个高斯分量线性叠加而成，其中为每个高斯分量的系数或权重
，
，
， 表示第个分模型的高斯分布密度。
作为一种生成模型，GMM提供了一种确定数据集最优成分数量的方法。由于生成模型本身就是数据集的概率分布，因此可以利用模型来评估数据的似然估计。
在本研究中，我们对矩阵使用高斯混合模型进行概率建模，矩阵中的数据，可以看做是从数个高斯分量中生成出来的。已知矩阵中有三种数据，所以将值设定为3，既本模型由、、三个高斯分量组成。那么必然有的数据服从，有的数据服从或者，不同于k-means等其他聚类算法，高斯混合模型不会直接输出每个数据属于哪一类，它会告诉我们某个数据属于每个高斯分量的概率值，这个概率值可以转化为本研究需要的“分数”。本研究利用高斯混合模型这一特点，把健康人所做动作的数据集和患者（模拟）所做动作的数据集放到同一个矩阵中，生成高斯混合模型，由此得到正常人和患者（模拟）所做动作的“分数”，这个分数被用来训练神经网络。
如公式⑤所示，建立高斯混合模型需要估计的参数为每个高斯分量的系数、均值以及协方差矩阵，估计这几个参数的过程被称作“参数估计”。如前文所述，数据点被组成一个360×5的矩阵，我们希望找到一组参数，它所确定的概率分布生成这些数据的概率最大，这个概率实际上为，这个乘积被称作似然函数。为避免浮点数下溢，取似然函数的对数，将乘积变成加和，得到对数似然函数。接下来只要将这个函数最大化，即找到了一组参数值使似然函数取得了最大值，完成了参数估计的过程。由于GMM的对数似然函数中有加和，不能直接用求导解方程的方法求最大值，所以采取EM（expectation maximization）算法求得GMM的参数值。EM算法分两步，第一步先求出要估计参数的粗略值，第二步使用第一步的值最大化似然函数。EM算法流程如下：
① 取参数的初始值开始迭代
② E步：根据当前模型参数，计算后验概率
， （6）
③ M步：根据E步中计算的后验概率计算新的参数
（7）
（8）
（9）
④ 计算（5）式对数似然函数
（10）
⑤ 重复②②步，直到对数似然函数收敛。
在E步中，公式③表达的意思是，对于每个数据来说，它由第个高斯分量生成的概率。式中和也是需要估计的值，所以采用迭代法，假设和已知，通过重复步骤迭代所得的值。假设E步中得到的就是正确的“数据由高斯分量生成的概率”，集中考虑所有的数据点，可以看做数据点由高斯分量生成，由于每个告诉分量都是一个标准的高斯分布，所以参数和可由公式 ⑦ ⑧求出。
本研究使用MATLAB实现了对降维后的数据建模，根据以上介绍的方法步骤求解高斯混合模型的参数。将高斯分量设置为6，既用经过63次迭代后，求得值如下（保留小数点后三位）：
0.144
0.234
0.162
0.123
0.146
0.191
的值组成一个6行5列的矩阵，表示第个高斯分量、第维数据的均值。求得值如下所示：
-0.149 0.587 0.436 0.316 0.008
0.253 0.287 -0.317 -0.052 -0.358
-0.333 0.373 0.225 -0.244 0.368
0.490 -0.802 0.239 0.151 -0.183
-0.560 -0.368 0.290 -0.178 -0.215
0.203 -0.321 -0.504 0.072 0.400
的值为6个5行5列的协方差矩阵，主对角线上的元素为方差，非主对角线上的元素为两两特征间的相关系数。6组的值如下：
0.5483 -0.1015 0.0474 -0.0292 0.0181
-0.1015 0.1149 -0.0208 -0.0575 0.0070
0.0474 -0.0208 0.0470 -0.0086 -0.0242
-0.0292 -0.0575 -0.0086 0.1654 0.0117
0.0181 0.0070 -0.0242 0.0117 0.0553
0.3162 -0.0092 0.0732 -0.0594 0.0571
-0.0092 0.2099 -0.0498 -0.0660 0.0723
0.0732 -0.0498 0.0610 0.0038 -0.0108
-0.0594 -0.0660 0.0038 0.1933 -0.0574
0.0571 0.0723 -0.0108 -0.0574 0.0995
0.6174 0.1803 0.0016 -0.0136 0.0375
0.1803 0.3940 0.0052 -0.0080 -0.0552
0.0016 0.0052 0.0135 0.0232 0.0113
-0.0136 -0.0080 0.0232 0.0889 0.0262
0.0375 -0.0552 0.0113 0.0262 0.0680
0.0672 0.0764 0.0574 0.0388 0.0011
0.0764 0.2044 -0.0055 0.1388 -0.0250
0.0574 -0.0055 0.1153 -0.0447 0.0369
0.0388 0.1388 -0.0447 0.2119 0.0089
0.0011 -0.0250 0.0369 0.0089 0.0758
0.6119 0.2256 0.1256 -0.1073 -0.0787
0.2256 0.1324 0.0490 -0.0453 -0.0603
0.1256 0.0490 0.1566 0.0011 -0.0329
-0.1073 -0.0453 0.0011 0.0814 0.0383
-0.0787 -0.0603 -0.0329 0.0383 0.0551
0.2748 -0.0466 0.1046 0.0086 0.0435
-0.0466 0.2118 -0.0734 0.0831 -0.0044
0.1046 -0.0734 0.1896 -0.0146 0.0893
0.0086 0.0831 -0.0146 0.1315 0.0154
0.0435 -0.0044 0.0893 0.0154 0.0918
N维高斯函数，其数据的维度也是N维的，此时均值向量是由N个均值向量组成，其物理意义仍然是高斯函数的中心，协方差矩阵依旧表示的是高斯函数的形状。
4.5 评分函数
本研究使用高斯混合模型完成了对动作的量化，度量公式如下：
（11）
（12）
公式（11）为高斯混合模型密度函数的对数似然形式，公式（12）被称为度量公式。此公式中，指一个数据集合中动作重复次数，根据前文分析，示例动作执行一次大概3秒钟，而选取的数据集中一共包含了6分钟的数据，既数据集中包含动作个数的期望值为120，所以取120。
之后采用组测试集的数据，代入到度量公式中，可计算出每个动作的度量值。计算结果如
图9.三组动作不同度量值
图9所示，*代表正常人动作度量值，☆代表病患1动作度量值，△代表病患2的动作度量值。很明显，正常人的度量值更为集中，病患的度量值有不同程度的分散，这是上肢运动不稳定的表现，度量值越分散，说明患者对自身的上肢控制能力越弱。
为将度量值映射为0到1之间的质量分数，针对正常人动作和患者动作，构建两个评分函数。由此产生的分数值在拟议框架中有着重要作用。第一，在现实环境的训练评估中，质量分数可更为直观的展现出性能度量，例如，提供一个88分的质量分数更便于患者理解，它还可以使患者通过一段时间内接受到的质量分数，方便的自我监督训练进展；第二，本研究通过质量分数对神经网络模型进行有监督的训练。对于正常人和患者的训练动作性能指标，提出以下两个评分函数：
（13）
（14）
在上述公式中
，，
，a1与a2为待定参数，通过演算，当
时可将得分值确定在所需范围；的作用是缩小度量值范围，确保所得分数接近1。通过公式（13）（14），计算结果如图10所示。同样的，*代表正常人动作得分值，△代表患者1动作得分值，☆代表患者2的动作得分值。计算结果满足预期，
图10.三组动作不同分值
正常人动作得分分布在0.9到1之间，患者动作分数呈现出不稳定、分值偏低的特点。
4.6 神经网络模型
在众多神经网络算法中，BP算法是最为常见的一种，它也作为本次研究中神经网络算法的根本，以BP算法为根本的前反馈神经网络，又被人们叫做BP人工神经网络，即BP网络。
图11.BP神经网络结构图
BP网络可以理解成一种有单一传播途径的复杂前反馈网络，这个网络中存在输入层、输出层和中间层（即隐含层）这三部分。在同一层级的神经节点里没有任何关联，各层之间可以做到全关联，经典的BP神经网络如图11所示。
该网络为（n，p，q）结构。此模型将输入其中的信号按一定的顺序通过输入层节点向每个隐藏层节点传递过去。接下来传到输出层。每个级别节点的输出只能影响它的旁边节点的输出。各个相邻的层级中的每个节点通过一定的连接权值实现连接。BP网络可以被当成一个由输入通向输出的高度非线性对应关系。
BP算法可理解为运用实际输出和期望输出的差值对网络的每层权值从后到前分部地完成修正的一个计算方式。根据对网络权值与阈值的适当修正，让网络性能函数，既学习目标函数，沿着梯度方向下降。
采集到的加速度、角速度数据作为神经网络的输入，根据高斯混合模型得到的分数值作为输出。希望得到的输出相比于实际输出的差值返回给神经网络。
前文提到，预处理之后的数据组成一个360×42的矩阵，为提高神经网络训练的速度，对这个矩阵做出调整，只保留1-3、13-15、19-21、31-33、37-39列数据，既三个传感器节点加速度和角速度的窗口均值，组成360×15的矩阵作为神经网络的输入层。构建的神经网络模型中囊括了输入层、两层中间层还有输出层，并且在中间层里面，第一层包含了20个神经元且第二次仅有单个神经元。设从输入层到中间层的连接权为，而从中间层再到输出层的连接权为；中间层内囊括的各个单元的阈值为；输出层当中的各个单元的阈值为。模型的期望输出表示为，实际输出为。模型中的隐含层里的激活函数都使用S型函数：
（15）
激活函数的导数为：
（16）
其误差函数为：
（17）
此模型能够依照误差的负梯度来相应的对进行调整，让其保持输出误差单调递减。
BP网络的学习，大体上可以分为四个过程：模式顺序传播，误差反向传播，记忆过程，学习收敛步骤。算法流程图如下：
图12.BP神经网络算法流程图
（学校时写的报告复制来的，有些符号没有显示，不要在意细节，大概能看到流程和思路。）
