基于步態信號與腦電信號的脊髓損傷大鼠運動信息分析

韓程韡，李想，石巖，王一軒

（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方法，從圖中可以看到，信號經濾波后，高能量都限定在了低頻區域，與我們的濾波結果相對應。所以短時傅里葉變換可以很好提取出大鼠腦電信號的特征，并且可以很好的對那個大鼠運動特征及步態特征，我們下一步計劃將大鼠步態信號與腦電信號進行解碼，獲取步態信息與腦電信息的對應關系。