基于Kinect的 帕金森病步態(tài)不對稱性識別方法①

張幼安，侯振杰，坎標(biāo)，姚恩，張家瑋

1.常州大學(xué)華羅庚學(xué)院，江蘇常州市213164；2.常州大學(xué)信息數(shù)理學(xué)院，江蘇常州市213164；3.常州大學(xué)機(jī)械工程學(xué)院，江蘇常州市213164

帕金森病(Parkinson's disease,PD)是一種特征性運動障礙，臨床存在多種表現(xiàn)，如慌張步態(tài)、凍結(jié)步態(tài)、靜止性震顫、肌僵直等。作為老年人多發(fā)病，國內(nèi)外學(xué)者做了很多努力預(yù)防或治療此病及其并發(fā)癥[1-9]。Gray等[10]研究顯示，帕金森病患者軀體運動功能下降，尤其是步態(tài)障礙與高死亡風(fēng)險相關(guān)。

臨床對帕金森病患者病情的評估主要采用統(tǒng)一帕金森病分級量表(Unified Parkinson's Disease Rating Scale,UPDRS)，但該方法參數(shù)較多，較為繁瑣[11-13]。實驗室一般通過穿戴設(shè)備或測量足底壓力，需要鋪設(shè)專業(yè)行道或在腳底墊測力墊，設(shè)備價格昂貴，也影響行走步態(tài)[14]。有很多學(xué)者通過視頻圖像提取數(shù)據(jù)，但建模過程一般較為繁瑣，無法推廣[15-17]。

鑒于步態(tài)障礙的主要特征為步態(tài)不對稱性，本研究提出一種基于Kinect骨骼數(shù)據(jù)的帕金森病步態(tài)不對稱性識別方法；針對帕金森病“N”字形發(fā)病規(guī)律，結(jié)合Kinect采集患者骨骼點數(shù)據(jù)，通過聚類算法得到行走加速度分布，使用字符串相似度的矩陣算法[18]，比較左右腳加速度分布，量化步態(tài)對稱性指標(biāo)，最后建立對稱步態(tài)隱馬爾科夫模型(Hidden Markov Model,HMM)完成常規(guī)識別工作。

1 對象與方法

1.1 研究對象

常州市武進(jìn)區(qū)南夏墅敬老院及湖塘鎮(zhèn)南夏墅村志愿者18例，采集時間為2016年7月至8月，調(diào)研過程為常州大學(xué)重點暑期社會實踐項目。

健康組納入標(biāo)準(zhǔn)：①年齡≥40周歲；②無神經(jīng)內(nèi)科疾病史；③無腰痛、關(guān)節(jié)痛等導(dǎo)致下肢活動異常的癥狀；④無明顯殘疾或步態(tài)障礙。

患者組納入標(biāo)準(zhǔn)：①年齡≥40周歲；②經(jīng)確診有帕金森病或有明顯帕金森綜合征癥狀，符合中國帕金森病的診斷標(biāo)準(zhǔn)；③可順利步行。

健康組10例，其中男性7例，女性3例；年齡40～65歲；身高158～176 cm。患者組8例，均為男性；年齡55～70歲；身高168～177 cm。

所有研究對象均知情同意。

1.2 方法

采用微軟公司Kinect V2.0，包含RGB彩色攝像頭、紅外線發(fā)射器和紅外線CMOS攝像頭；每秒最多獲取30幀圖像；可識別距離0.5～4.5 m。

Kinect骨骼跟蹤技術(shù)通過處理景深數(shù)據(jù)，建立人體各個關(guān)節(jié)的坐標(biāo)。骨骼追蹤能確定人體各個部分，如頭、腿等。骨骼跟蹤產(chǎn)生X、Y、Z軸坐標(biāo)，確定這些骨骼點的位置。因此，Kinect提供的骨骼數(shù)據(jù)為行走過程中身體各關(guān)節(jié)點的空間坐標(biāo)序列。使用時通過調(diào)用自帶SDK進(jìn)行實驗對象行走過程中的全身骨骼點采集。

設(shè)立步道長7 m、寬2 m，排除外界干擾。將Kinect平放在步道一端。受試者以正常舒服姿勢站立于Kinect正前方6 m處，操作者對Kinect進(jìn)行距離和角度等參數(shù)校準(zhǔn)。當(dāng)Kinect指示燈亮起，受試者以最舒適狀態(tài)面向Kinect走去，直線行走6 m或至無法識別人體關(guān)節(jié)點坐標(biāo)后停止。每例受試者測試6次。

采集數(shù)據(jù)為全身25個關(guān)節(jié)點的三維坐標(biāo)序列，通過該骨骼點序列提取步態(tài)周期，繼而提取步態(tài)參數(shù)，進(jìn)行對比實驗，驗證設(shè)備可靠性。計算行走過程中位移和加速度序列，聚類分析患者左右腳步態(tài)對稱性，最后完成常規(guī)識別。

1.3 觀察指標(biāo)

采用MATLAB軟件編程提取步態(tài)特征。提取步態(tài)參數(shù)包括步長、步幅、步寬和聚類后的加速度序列。

1.3.1 步態(tài)特征提取

先由空間坐標(biāo)序列的周期性變化計算步態(tài)周期，再通過步態(tài)周期計算各項步態(tài)參數(shù)。

1.3.1.1 步態(tài)周期計算

使用Kinect采集的人體骨骼點空間變化坐標(biāo)序列。由于Kinect攝像頭擺放無法使X軸與地面完全平行，也無法以Y軸的0值視為落地點，所以不能根據(jù)返回坐標(biāo)(如Y軸的0值)直接確定雙腳落地點來確定周期。

行走過程中，人體中心點呈現(xiàn)左右、上下擺動的周期性變化。本文將每次擺動的最低點作為特殊點，出現(xiàn)特殊點的幀作為特殊幀。通過求特殊幀位置即可計算步態(tài)周期。步態(tài)周期計算流程圖如圖1。

1.3.1.1.1選取中心點

行走過程中，上半身產(chǎn)生的彈性變形接近0，彎腰程度接近0，上半身可近似剛體；而下半身由于運動過于復(fù)雜，無法準(zhǔn)確描繪中心點周期性變化；結(jié)合Kinect V2.0的采集點，取脊柱底部點作為中心點(圖2)。

1.3.1.1.2數(shù)據(jù)低通濾波

采集取中心點Y軸坐標(biāo)序列。由于Kinect采樣頻率默認(rèn)為30幀/s，因此時間域和頻域一致。做該點Y軸坐標(biāo)頻域變化曲線圖，手動提取步態(tài)穩(wěn)定區(qū)間，一般為前80%點。

第1次濾波：做相鄰數(shù)據(jù)差值，若該值大于前n組差值和，則將該值置為前n組差值的均值。由于Kinect采集性能較穩(wěn)定，n置為5。

第2次濾波：對坐標(biāo)序列行RLOWESS濾波，得到新的Y坐標(biāo)變化序列。

1.3.1.1.3求特殊幀

根據(jù)濾波后Y坐標(biāo)序列，做變化曲線圖。計算曲線上每一點的瞬時變化率，得到瞬時變化序列。在13～17之間依次取數(shù)作為周期基數(shù)，套取每個周期內(nèi)的瞬時變化率最小值。根據(jù)相鄰套取結(jié)果的差值判斷：若與其他相鄰套取結(jié)果差之間存在倍數(shù)差關(guān)系，則該次套取失敗，重新選循環(huán)基數(shù)，直致篩選出正確特殊點，并求得特殊點序列。特殊點所在幀數(shù)即特殊幀位置。

人正常行走速度 V=4～4.7 km/h，即 1.1～1.3 m/s，以S=0.6 m作為正?？绮介L，計算步頻為1.85～2.20/s。

Kinect采集頻率30幀/s，所以n=13.64～16.22幀/周期，因此Kinect采集的步態(tài)周期數(shù)在13～17之間，所以將13～17作為周期基數(shù)。

1.3.1.1.4計算步態(tài)周期

步態(tài)周期(T)為行走過程中一側(cè)足跟著地至該側(cè)足跟再次著地時所經(jīng)過的時間，即連續(xù)3幀特殊幀出現(xiàn)所經(jīng)歷的時間。

其中，tk為第k個特殊幀出現(xiàn)的時間，tk+2為第k+2個特殊幀出現(xiàn)的時間。

圖1 步態(tài)周期計算流程圖

圖2 Kinect V2.0的采集點

1.3.1.2 步態(tài)參數(shù)

步態(tài)參數(shù)定義見圖3，步長、步幅、步寬用來驗證Kinect采集步態(tài)的可靠性，加速度序列用來進(jìn)行對稱性檢驗和步態(tài)識別。

圖3 步態(tài)參數(shù)定義

1.3.1.2.1步長

行走時一側(cè)足跟著地到緊接著的對側(cè)足跟著地的距離。使用空間點坐標(biāo)距離公式計算：

其中，spi為第i次跨步步長，Xk、Yk、Zk為特殊點序號為k的空間坐標(biāo)，Xk+1、Yk+1、Zk+1為特殊點序號為k+1的空間坐標(biāo)。

1.3.1.2.2步幅

行走時由一側(cè)足跟著地到該側(cè)足跟二次著地所行進(jìn)的距離。

1.3.1.2.3步寬

先以一側(cè)足(本文選右足)的2個落地點和1個懸空點，通過三點式方程構(gòu)造平面空間；然后計算對側(cè)足(本文選左足)落地點至該平面的距離。

兩個右腳落地點R1(X1,Y1,Z1)和R2(X2,Y2,Z2)，懸空點R3(X3,Y3,Z3)建立的空間平面方程：

對側(cè)足落地點到平面的距離，即步寬：

其中，(X0,Y0,Z0)為對側(cè)落地點坐標(biāo)。

1.3.1.2.4位移

位移的計算使用下一時刻點坐標(biāo)至該點坐標(biāo)的直線距離。

1.3.1.2.5速度

每次采集捕獲骨骼點的瞬時速度。由于采集頻率為30幀/s，所以相鄰幀時間間隔為0.033 s，用相鄰3個點的平均速度代替中點的瞬時速度，中點的瞬時速度：

其中，fc為 Kinect采集頻率，δt-1為點 Bt-1到點 Bt的距離，δt為點 Bt到Bt+1的距離。

1.3.1.2.6加速度

相鄰幀速度的變化率：

根據(jù)采集到的m+1個點的坐標(biāo)，可以確定時間參數(shù)序列，即位移序列X1,X2,X3,……Xm；速度序列V1,V2,V3,……Vm-1；加速度序列a1,a2,a3,……am-3。

1.4 對比實驗

在地面鋪設(shè)白紙，使用米尺測距、秒表測時，采集一次步態(tài)周期實際步長、步幅、步寬，并計算速度。

1.5 對稱性檢驗

臨床研究表明[19]，出現(xiàn)凍結(jié)步態(tài)的帕金森病患者姿勢控制更差。帕金森病患病呈“N”字型進(jìn)展：首先單側(cè)腿出現(xiàn)病癥，隨后進(jìn)展至同側(cè)上肢，再發(fā)展至對側(cè)腿部，再至對側(cè)上肢。推測初期單腿患病時，肌肉控制力下降，患者會多使用健康腿，導(dǎo)致步行周期左右腿加速度不同；對比兩腿加速度序列便可以比較兩腿步態(tài)對稱性。

本研究通過比較左右腳加速度分布字符串的相似度，確定步態(tài)的對稱性，對比參數(shù)為左右膝的加速度序列。在確保加速度序列誤差足夠小的基礎(chǔ)上，分別對左、右序列進(jìn)行K-means聚類，得到左、右一次聚類字符串；取一側(cè)聚類結(jié)果不變，用另一側(cè)所有可能的聚類結(jié)果字符串與前一側(cè)聚類字符串用字符串相似度矩陣算法，求相似度，取相似度最大值作為實際對稱性指標(biāo)。

1.6 常規(guī)識別

建立對稱步態(tài)HMM模型進(jìn)行常規(guī)識別。根據(jù)中心點波形分析，左右腳加速度分布會表現(xiàn)在中心點的加速度分布上。

收集健康組中心點加速度序列，使用Baum-Welch算法訓(xùn)練HMM模型參數(shù)；通過前向算法統(tǒng)計不同程度帕金森病患者步態(tài)序列的發(fā)生概率，提取最小概率作為閾值。計算受試者步態(tài)序列相對于健康人HMM模型的發(fā)生概率，若在患者閾值范圍內(nèi)，則記為患者。

在應(yīng)用HMM模型預(yù)測步態(tài)加速度分布時，我們將已知加速度序列看作是可見狀態(tài)鏈，將人體腿部各動作激勵的肌肉或神經(jīng)信號作為隱含狀態(tài)鏈。

2 結(jié)果

2.1 可靠性檢驗

使用Kinect采集的人體步行1次骨骼數(shù)據(jù)，其中心點Y軸坐標(biāo)頻域如圖4。

圖4 中心點Y軸坐標(biāo)頻域變化

取前80%點作為穩(wěn)態(tài)序列，對穩(wěn)態(tài)序列兩次濾波，作出穩(wěn)態(tài)序列圖像，并在同一圖上作出濾波后的曲線，如圖5。

對曲線中每一點求斜率，并作斜率圖，如圖6。

圖5 穩(wěn)態(tài)數(shù)據(jù)及濾波圖

圖6 斜率圖

一名測試者Kinect采集的步態(tài)參數(shù)，與手工測量步態(tài)參數(shù)比較見表1。Kinect采集的三種步態(tài)參數(shù)誤差均在5%以內(nèi)，提取步態(tài)參數(shù)的方法可以信任。見表1。

表1 兩種測量方式比較(cm)

2.2 步態(tài)對稱性檢驗

展示1例健康人數(shù)據(jù)和3例最有區(qū)別度的患者數(shù)據(jù)。由于每次實驗采集骨骼點超過100條序列，僅顯示前7次數(shù)據(jù)，如表2。其中患者1為拖沓步態(tài)，患者2為輕度跛行，患者3為單腿萎縮。

表2 膝關(guān)節(jié)加速度原始數(shù)據(jù)(m/s2)

加速度狀態(tài)可分為5類：加速度=0、加速度減小的加速運動、加速度增大的加速運動、加速度減小的減速運動、加速度增加的減速運動。因此，K-means聚類時設(shè)K=5，結(jié)果見表3。

表3 膝關(guān)節(jié)加速度聚類結(jié)果

左右腳狀態(tài)相同的速度聚為同一類，將一側(cè)聚類后的結(jié)果進(jìn)行循環(huán)種類號替換，聚類結(jié)果不同的可以任意替換(如1234可任意替換為1324、4123等)，聚類結(jié)果相同的只能替換成相同的聚類效果(如1233只能替換為2144、4311等)。一側(cè)加速度序列有120種相似度結(jié)果，提取相似度最大值作為該受試者的實際相似度。結(jié)果見表4。

表4 健康人與帕金森病患者左右腳加速度分布相似度比較

健康人左右相似度最大，患者1、患者2、患者3病情依次加重，左右相似度依次減小。

2.3 常規(guī)性識別

采用Baum-Welch算法學(xué)習(xí)HMM模型，加速度序列通過DBSCAN聚類，結(jié)果有6類，即HMM模型中Markov鏈每個狀態(tài)對應(yīng)可能的觀察值數(shù)目為6；預(yù)設(shè)Markov鏈狀態(tài)數(shù)目N，初步訓(xùn)練得到模型參數(shù)，再運用前向算法計算由該初步模型產(chǎn)生加速度序列的概率P；不斷重設(shè)N值，重復(fù)該步驟至概率收斂，得最佳N值為5。

結(jié)果：M=6，N=5

其中，M為每個狀態(tài)對應(yīng)的可能的觀察值數(shù)；N為模型中Markov鏈狀態(tài)數(shù)；A為狀態(tài)轉(zhuǎn)移概率矩陣；B為觀察值概率矩陣；Pi為初始狀態(tài)概率矢量。

采集計算多組帕金森病癥患者的加速度序列，聚類標(biāo)準(zhǔn)化，前向算法多次計算后

某受試者經(jīng)加速度序列聚類標(biāo)準(zhǔn)化后，前向算法計算得logP=-4.2×10-2，高于閾值，則判斷該受試者有可能患帕金森病。

對8例患者重新采集數(shù)據(jù)，共8組26條數(shù)據(jù)(一些患者由于身體原因無法重復(fù)多次采集)，并提取出能代表其特征的7組14條數(shù)據(jù)，成功識別出12條，成功率85.71%。提取了10組60條健康人數(shù)據(jù)，在有效的46條數(shù)據(jù)中，成功區(qū)別出35條，區(qū)別度76.09%。

3 討論

根據(jù)帕金森病不對稱步態(tài)進(jìn)行患病識別一直受到關(guān)注。王人成等[20]采用基于FITTS定律提出的二維極坐標(biāo)下的步態(tài)對稱性指標(biāo)；Shorter等[21]采用偏差閾值分析，比較不對稱行走的屈伸行為與未裝支具組的標(biāo)準(zhǔn)，使用該方法可量化雙側(cè)關(guān)節(jié)的對稱性和不對稱性；鄒亞等[22]根據(jù)可穿戴慣性傳感器觀察參數(shù)提出對稱性計算式。

多數(shù)學(xué)者主要以左右腳參數(shù)比值做主要因式[23-24]，但對于不對稱步態(tài)的識別依據(jù)一直缺乏獨到且可執(zhí)行的解釋，尤其是帕金森病患者不對稱性步態(tài)。廖福元等[25]將患者左、右腿支撐時間和跨步時間排成序列對，以多個尺度逼近序列對的相似度來量化步態(tài)的對稱性，但沒有結(jié)合病理提出識別原因。

本研究基于帕金森病“N”字形發(fā)病規(guī)律，將使患者患病初期產(chǎn)生較為明顯的不對稱步態(tài)，通過Ki-nect體感設(shè)備建立人體連桿模型，采集關(guān)節(jié)點，使設(shè)備對正常行走步態(tài)的影響降到最低。

本研究通過骨骼點三維坐標(biāo)序列類正弦曲線的周期性特征，結(jié)合曲線斜率大小求波峰、波谷的方法提取周期，然后提取步態(tài)特征參數(shù)。對比足跡測量，證明該提取方法可行。周期提取還可以通過計算斜率方差的方式提取波峰、波谷，再提取周期。在提取周期后，除可計算本文所用的特征參數(shù)之外，還可以通過空間坐標(biāo)的時間變化曲線提取其他步態(tài)參數(shù)，如步頻、步速、支撐相時間等。

對于帕金森病的步態(tài)，本文基于帕金森病患者左右腳加速度分布不同的假設(shè)，通過字符串相似度的矩陣算法識別不對稱步態(tài)。由于K-means聚類方法對噪聲不敏感，誤差可能較大；也可選擇使用DBSCAN聚類，但需人工多次干預(yù)，直到聚類到所需的聚類書目。在使用字符串相似度的矩陣算法時，由于需要對一側(cè)聚類結(jié)果進(jìn)行窮舉，再與另一側(cè)進(jìn)行匹配算法，較為復(fù)雜，有待改進(jìn)。最后，我們使用HMM預(yù)測模型進(jìn)行患者和非患者的識別，證明有一定效果。但HMM模型的訓(xùn)練常需要較大數(shù)據(jù)，因此需要插值算法，或需要將幾次采集結(jié)果組合，會產(chǎn)生一定誤差，還需要進(jìn)一步研究。