基于步态信号与脑电信号的脊髓损伤大鼠运动信息分析

韩程韡，李想，石岩，王一轩

（1.北京航空航天大学自动化科学与电气工程学院，北京 100191；2.中国传媒大学，媒体融合与传播国家重点实验室，北京 100024）

1 引言

医学信号是人体所表现出来的物理信息和化学信息。正确处理医学信号无论对认识人体自身还是应用到工程技术当中，都至关重要。在工程中，脑机接口技术是医学信号处理的重要应用场所，而数据处理与分析系统则是脑机技术中起决定性作用的环节。在智能假肢等应用中，正确分析脑电数据和步态数据的信号，探求二者间的关联，是实现控制系统的关键。

步态是人或动物在行走过程中表现出来的相应特征和模式，观察并研究人或动物的行走过程并从中提取出相应的特征进行分析则被称为步态分析。因行走是一种较为复杂的运动，由生物体大脑、神经、骨骼、关节、肌肉等多环节配合才能完成，且受生物身高、体重、日常习惯等多方面因素影响，故步态分析在生物医学、信息科学等多领域有着广泛应用。在一般的研究中，研究者多通过计算机视觉或压力分析来研究步态特征，如在山东大学纪阳阳的硕士毕业论文中选用了运动视觉的方式来获得人体行走时的步态信息［1‑2］，安徽大学的李彦琳在硕士学位论文中则是根据足底触觉和压力分布分析出人的步态特征［3‑4］。这些方法在应用时所需的数据量大，计算复杂。若数据不经处理直接应用到算法当中，则可能会使计算量加大，运算时间过长。此外因未经处理的数据间相关性高、信息冗余也可能会造成过拟合现象，从而导致对步态分析的结果出现偏差。此时主成分分析法可以很好地提取步态数据中最主要的信息，将数据投影到方差最大的空间中，减少数据间的相关性，更好地分析和提取出步态特征。

人的行走过程，离不开大脑活动的控制。当大脑在活动的时候，脑部的神经元放电产生电位，形成了可以测得的脑电信号。脑电信号作为一种机理复杂的非平稳信号，有着信号强度弱，采集过程中引入噪声大、频率范围低、随机性很强的特点。研究发现，与运动感知相关的信号多集中于低频区域［5］。因此，选择正确的滤波方法可以更好分析出运动相关脑电信息。例如在Hammer等人的研究中通过研究脑电信号的低频成分，很好地分析出位置与速度等相关的运动信息［6］。经滤波处理后的脑电信号在经过不同方法的特征提取后，可以更突出信号运动特征。不同的特征提取方法可针对信号不同层面的特点进行关注，现在有的主要特征提取方式可分为频率域、时间域、时间‐频率域和非线性动力学等方法。例如在Gotman的文章中，他通过分析信号的频率谱，提取信号在频率域上的特征。现在应用的滤波器中，因选择滤波器本身的特性，会因其幅值特性的起伏使滤波后信号变形，产生幅值偏差［7］。我们在研究中选择巴特沃斯滤波器，其通频带内的频率响应曲线最大限度平坦的特性，让过滤后的信号最大程度保留原有特性。在完成信号的特征提取时，传统的单纯时间域和频率域上的特征提取不能保证信号在该维度上拥有最清晰的特征，提取出最大的信息，而时间‐频率域的特征提取方式结合了两种信号分析方式的优点，最大程度体现信号的信息。

2 步态数据处理

为了减少数据相关度和冗余的信息，我们需对数据进行降维处理。在这里我们选用主成分分析法对数据进行降维。在减少数据维度的同时尽最大可能保留数据原有信息。

在降维处理数据之前，我们需要进行数据读取工作。原始数据均存在txt文本中：第一列是时间，后续15列分别为骨盆、臀部、膝盖、脚踝、脚掌这五个点在空间中的x、y、z坐标，每列数据以制表符相隔。我们用pandas包中的read_csv函数将数据读取为dataframe的表格形式后，再将时间和五个点的三维坐标分别存在六个列表。

读取数据完成后，我们即对数据进行降维处理。我们用主成分分析（PCA）法，将数据从三维降至二维。在操作上我们先将945组数据中有效数据筛选出来，即去除掉数据中含有NAN的数据，并将缩短后的数据存至新的列表中。如左腿骨盆的945组数据，在去除无效数据后剩余674组数据，每组数据是三维的坐标点以［x,y,z］形式组成的列表。我们将缩短后的有效数据看做三维的列向量，在去均值处理后求出其协方差矩阵。此矩阵为3*3的对称矩阵，求出其特征值和对应特征向量后，保留数值最大的两个特征值所对应的特征向量并进行归一化处理，组成坐标变换矩阵（矩阵大小为3*2）。在此变换矩阵的作用下我们可以得到一个674行2列的矩阵，可以看做674组2维数据。至此我们将数据从三维降至二维。

处理完数据后我们将降维后的二维数据在之前删掉的NAN的地方补充上0，使其变成和时间相同长度的列表。在三维图上，我们将二维数据在时间轴的配合下在一张三维图上进行呈现，并在同一时刻的不同点中间用线段链接，以直观地体现出大鼠腿部运动随时间的角度变化。同样思路，我们也可以将数据在二维图像上进行展示：将x轴和时间轴融合。将每组数据的x坐标转换至时间步长范围内，在结合数据对应时间点，可将数据在二维图上展示出来。

因大鼠腿部共有五个标注点对应的五组坐标，其中位于中间的臀部、膝盖和脚踝这三个点可对应三个关节来求出大鼠运动过程中三个关节的角度变化。具体思路为：将一个关节点和其临近的两个点分别组成两个向量，通过计算两个向量间夹角可以求出这些关节间实时角度变化，并通过二维折线图进行输出，如图1所示。

一共要进行两组实验数据结果分析，第一组为正常大鼠，第二组为进行脊髓损伤打击后的大鼠。对该两组实验视频进行提取后获得相应的实验坐标数据。之后分别对两组进行降维并获取角速度然后进行对比，来获得相应的实验结果。获得实验及数据预处理结果之后，做了一个仿真图像，将所获数据用matlab图象化出来。图2所示为进行数据预处理之后获得的步态数据。在图2中，最上面的小图是健康状况下的大鼠二维数据步态运动模型，中间的图为做完脊髓损伤后，最下面的图为脊髓损伤手术后恢复一段时间的大鼠运动模型。从图中可以看出，在健康大鼠脊髓损伤手术后，大鼠步态幅度减小明显，几乎没有太大变化，大腿根部及膝盖部分保持位置几乎不变。在健康状况下,大鼠能保持一定频率的步态行走，脊髓损伤之后，只有轻微的行走过程，且并不明显。在恢复一段时间后，基本的步态幅度保持不变，但是可以看出，步态的规律性变差，经常会有不同幅度的不规律动作，其最大幅度远大于正常状况下，小幅度也会小于正常大鼠的变化。

图2 正常大鼠与脊髓损伤大鼠腿部位置随时间变化对比图（从上到下依次为健康、脊髓损伤、脊损恢复）

除了图2，我们还通过了matlab获得了每个坐标点的速度变化图像。而速度获得公式为

其中，x和y分别为降维后相邻两坐标之间的差值。从上述图中，我们可以了解到，在刚刚做完脊髓损伤手术后的大鼠的步态数据速度变化并不剧烈明显。在刚刚做完手术后的大鼠步态单点速度图没有参考意义。在康复一段时间后，我们获得以下数据：

在图3中，最上面的图为健康大鼠各个点速度变化，中间的图为刚刚完成脊髓损伤手术后的大鼠速度图，下面的图为做完脊髓损伤手术后并康复一段时间后的大鼠各点速度变化。很明显可以看出，大鼠在健康状态下，速度变化比较规律，有一定的周期性，在脊髓损伤情况下，大鼠步态行走产生了不规律的变化，且变化范围增大。正常情况下，速度在0‑0.7之间变化，而在术后恢复后，变化范围扩大了0.5.且并不规律。

图3 正常大鼠与脊髓损伤术后恢复后大鼠步速对比图（从上到下依次为健康、脊髓损伤、脊损恢复）

通过以上数据可以看出。大鼠脊髓损伤之后，其步态行走受到了深远的影响。原本能够行走的有规律，且变化程度小。在脊髓损伤之后，步态行走变得没有规律，且变化较大。这说明脊髓与大鼠的步态行走密不可分，对其行走的规律性有着很重要的影响。

如图4所示，可以看出，健康情况下，角度变化比较规律，刚刚术后，幅度变化频率小。

图4 正常大鼠与脊髓损伤术后恢复后大鼠腿部角度对比图（从上到下依次为健康、脊髓损伤、脊损恢复）

3 脑电数据处理

3.1 脑电滤波

巴特沃斯滤波器最大的特点是满足其幅值平方要求。即

此中滤波器在性能上有着显著的优势。如下图5所示，相比较别的IIR滤波器，它原理简单、设计便利，具有通频带内频率响应曲线最大限度平坦的优点。

图5 巴特沃斯滤波器和其他经典IIR滤波器的相频特性比较

我们采用的脑电信号为成年大鼠脑部植入电极后所采集到的信号。信号长度为20s 左右，采集频率为40000 Hz。在采集信号过程中，大鼠被安置在走步机上，可随履带的运动速度前进。在20s左右的信号录入过程中，发生过三次很明显的运动。分别在视频中4.5s‑6s,10s‑11s,16s‑19s.

在信号处理之前，信号原始图片随下图6所示。可见信号有大概三段比较剧烈的抖动，整体信号不尽平稳，信号起伏相差较大，无法探求信号的特征所体现的频率范围，且信号的运动特征在频率和幅值中隐藏的关系无法明晰体现。

图6 脑电信号处理前图像

在处理信号的过程中，我们采用了4阶截止频率为8‑30 Hz的巴特沃斯带通滤波器。滤波后的信号如下图7所示。

图7 经过8⁃30 Hz 4阶巴特沃斯滤波器处理后的脑电信号图像

在滤波中，我们在这个区间段主要关注alpha波和beta波段的脑电信号。可以看到在处理后的信号基准线被拉平，很好过滤了原始信号中的高频和低频的噪音。为后续进行脑电信号的解码奠定了很好的基础，减去了大面积的杂音干扰。

此外，经过滤波后的信号波动区间界定更加清晰，将频率限定在一定区间后，统一了原始信号中幅值和频率蕴含信息模糊的问题，将运动信息更加聚集在信号幅值当中。

3.2 脑电特征提取

短时傅里叶变换（STFT）是傅里叶变换的一个变形。傅里叶变换可以很好的将时域信号变换到频率域，以探求信号的频率空间上的特征。但是傅里叶存在的前提是所出来的信号必须是平稳信号，而对于脑电这样的典型非平稳信号来讲，傅里叶变换是不存在的。

而短时傅里叶变换的原理是从傅里叶变换的思路出发，将非平稳随机信号截取成一系列平稳信号，用窗函数在原始函数上的平移，完成一系列的傅里叶变换，最后得到函数的时频域特征。

如下图8所示，我们对滤波后的函数进行了窗函数长度为100的短时傅里叶变换，得到信号的时频域特征结果，并将结果通过频谱图表示出来。

图8 短时傅里叶变换进行特征提取后的信号图像

在图像中横坐标是时间轴，纵坐标是频率轴，而信号的颜色则表示在这一点上信号的“能量”，即反映了幅值的大小（dB）。从图中可以看到，信号经滤波后，高能量都限定在了低频区域，与我们的滤波结果相对应。

同时，除了频率域信息外，还可以从横向读出时间域上的信息，可以明显看到在对应大鼠运动的三个时间段上，图像颜色更黄，信号能量很大，与大鼠运动相对应。所以短时傅里叶变换可以很好提取出大鼠脑电信号的特征。

我们可以看到，相比较单独的频率域特征和时间域特征来看，时频域特征可以更好地结合两者的优势，提取出信号更详尽的信息。但同时在分析图像中，可以看到色块较大，此方法的分辨率并不算很高。并且，从其原理出发，并不一定窗函数所画区间就和平稳随机函数区间相对应，这也是影响此方法准确率的潜在威胁之一。

4 结论及展望

获得了大鼠的步态数据，并对数据进行了降维处理，通过数据我们了解到了，脊髓损伤对大鼠步态运动的规律性有了很重要的影响。在脑电数据中，通过处理后的内容可以看到，通过STFT方法，从图中可以看到，信号经滤波后，高能量都限定在了低频区域，与我们的滤波结果相对应。所以短时傅里叶变换可以很好提取出大鼠脑电信号的特征，并且可以很好的对那个大鼠运动特征及步态特征，我们下一步计划将大鼠步态信号与脑电信号进行解码，获取步态信息与脑电信息的对应关系。