方 震 白忠瑞 陳賢祥 夏 攀 何征嶺 趙榮建
①(中國科學院空天信息創(chuàng)新研究院 北京100190)
②(中國科學院大學電子電氣與通信工程學院 北京100049)
心沖擊圖(BallistoCardioGram,BCG)是一種無接觸式的心臟活動監(jiān)測手段。它來源于心臟泵血過程中血液流動在人體內部產生的一系列機械沖擊力[1]。與心電圖(ElectroCatdioGram,ECG)等其他生理信息監(jiān)測技術相比,BCG避免了電極線和電極貼長時間貼在身上帶來的束縛和不適,具有無創(chuàng)、無直接接觸和檢測方便等優(yōu)勢,十分適用于居家養(yǎng)老、病房等需要長時間監(jiān)測生理信號的場景。近年來隨著傳感器和數字信號處理的發(fā)展,BCG逐漸受到研究人員的重視,研究表明,BCG可以應用于心率、心率變異性(Heart Rate Variability,HRV)監(jiān)測[2,3],也可應用于心臟收縮性以及心輸出量變化等心血管功能的評估[4]。
但BCG信號在不同測試條件(如測試設備的安置位置與角度、受試者姿勢等)下表現出較大的不一致性[5],此固有特征在一定程度上限制了分析BCG的形狀特征在評估心血管功能方面的應用潛力。相比之下,通過檢測逐拍心動周期估計心率或者計算HRV指標,是BCG可行性較高的一種應用路徑。
HRV定義為逐次心跳周期差異的變化情況,一般通過分析心電信號中R波間期序列獲得。HRV反映了心臟本身竇性心律不齊的程度以及神經體液因素與竇房結之間相互作用的平衡關系。已有的研究表明,HRV是心源性猝死、冠心病、高血壓病及慢性心力衰竭等心血管疾病及慢性阻塞性肺疾病、糖尿病等疾病預后的預測因子[6,7],還能反映出睡眠狀況、精神壓力狀況等多種信息[8,9]。因此使用無感式的BCG信號獲取心率變異性指標,對于用戶方便舒適地掌握自身健康狀況具有重要意義。
通過BCG獲取HRV的關鍵在于獲得精準可靠的逐拍心率。近年來,已有很多學者在非接觸式獲取逐拍心率方面做了研究:在BCG信號采集裝置方面,有壓電薄膜、壓電陶瓷、光纖、毫米波雷達、加速度計、壓變電阻以及各類采用MEMS工藝的微傳感器等[10–12];在逐拍心動周期的算法方面,目前有自適應閾值法、峰值探測法、倒譜分析法、多示例學習法、模板匹配法、混合判決方法[13–16]。值得注意的是,目前此方面研究更多地關注逐拍心率檢測的覆蓋率和平均心率的準確度,但事實上,逐拍心動周期計算的精度,即平均絕對誤差(Mean Absolute Error,MAE)對于心率變異性的計算也尤為重要。例如,在5 min短時HRV計算中,若BCG逐拍心動周期序列平均絕對誤差為10 ms,則由其計算得到的頻域特征HF相比ECG真實值的偏差可達15%~25%,SDNN,p NN50等時域指標偏差也可達10%左右,且此偏差具有隨機性,這使得BCG所得到的HRV特征指標的可靠性降低;而當逐拍心動周期的MAE降至4 ms時,逐拍心動周期序列所得到的HF偏差降低到8%以內,SDNN和p NN50等指標偏差降低至2%。
針對目前大多數方法在逐拍心率計算精度方面的不足,本文完成以下工作:(1)通過采用優(yōu)化的傳感前端設計,增加傳感器的靈敏度和BCG信號的時間分辨率;(2)通過對心沖擊信號進行分析,找到BCG中最適合提取逐拍心動周期的關鍵成分;(3)提出一種采用聚類的自適應模板匹配算法,以準確提取心動周期信息。最后通過實驗驗證了方法的準確性。
BCG信號采集裝置的設計影響逐拍心率提取的準確度。在各種BCG監(jiān)測方式中,壓電陶瓷傳感器由于其靈敏度高、價格低廉、設計靈活等優(yōu)勢,具有良好的應用前景。本文通過壓電陶瓷非接觸式采集BCG,計算逐拍心率。
本文所設計的傳感前端外殼為一個ABS材質塑料圓盒,壓電陶瓷片被固定在圓盒內上壁形變最大位置。外殼上蓋結構設計為一側凸起形狀,使其受到震動時響應更靈敏。實驗中將傳感前端放置在胸口下方的床墊下,采集到的身體震動信號會經放大、濾波以及模數轉換傳輸到帶有低功耗藍牙(Bluetooth Low Energy,BLE)功能的主控芯片中進行下一步的傳輸和處理。BCG信號的采樣頻率也影響逐拍心率計算的精準性,經對比測試,本文設置BCG采樣頻率為250 Hz,以保證逐拍心率計算的準確程度,同時最小化資源消耗。
傳感前端采集到的原始的壓電信號可用式(1)表示
其中,r(n)表示呼吸引起的振動成分,b(n)表示BCG成分,g(n)表示噪聲成分。心臟的跳動是一個較為復雜的機械運動過程[17],BCG中可以反映心臟機械運動和心臟泵血引起的身體機械運動的很多信息。
但對于一般所測得的BCG信號,由心臟原始的泵血引起的信號相對細微,容易被淹沒在由身體四肢震動造成的信號的主能量中,通過不同的處理方式可以凸顯出BCG中的不同成分。圖1是不同方式處理BCG信號的結果,虛線是同步采集的ECG的R波所在位置。圖1(a)是經過48 Hz低通濾波器預處理去除高頻噪聲后的壓電信號,為了求得誤差最小的逐拍心動周期誤差,我們嘗試了多種不同的方式對此信號進行進一步處理,以獲取BCG信號中最適合進行精準逐拍心率提取的“關鍵成分”,下面優(yōu)選幾種具有代表性的方法具體說明:
方法1使用3~24 Hz的巴特沃斯帶通濾波器進行濾波,得到BCG概貌波形,其結果如圖1(b)。
方法2使用與心沖擊信號最相符的db6小波基對原始壓電信號進行9級小波分解,并使用小波分解后的第3~7層細節(jié)分量重建BCG信號,其結果如圖1(c)所示。
方法3使用8~24 Hz的巴特沃斯帶通濾波器進行濾波,可以得到BCG信號中的由心臟跳動引起的細微振動成分,其結果如圖1(d)所示。
方法4對方法1處理的結果進行二次差分操作,求取其加速度信號,結果如圖1(e)所示。
對于不同的BCG處理方法,每拍心跳在信號上造成的最高點并不是同一位置,為方便起見,本文中統稱這些點為BCG的“J點”。
2.3.1模板學習
圖1 BCG不同成分提取結果
BCG形狀多變,難以提前預知其波形。我們使用了親和傳播(Affinity Propagation,AP)聚類算法[18],從10 s的BCG信號中進行模板的自學習。具體步驟如下:
(1)提取2.2節(jié)各方法處理過的10秒BCG信號段中所有的極大值點,以這些極大值點為中心,截取其前后各100采樣點,共約0.8 s長度的信號,得到數據樣本集B={b0b1···bn},其中bi代表第i個極值點處截取的信號。
(2)使用AP聚類算法對樣本集B進行聚類,并得到模板聚類。AP聚類的具體步驟如下:
步驟1初始化吸引度矩陣R和歸屬度矩陣A為0矩陣。其中R中元素r(i,k)表示數據bk適合作為數據bi的聚類中心的程度,A的元素a(i,k)表示對象bi適合選擇數據對象bk作為其聚類中心的程度。
步驟2分別按照式(2)和式(3)迭代更新吸引度矩陣與歸屬度矩陣。
其中,s(i,k)是相似度矩陣S的元素,為弱化BCG幅值變化的影響,本文沒有使用歐氏距離,而是設置s(i,k)為Pearson相關系數與常數1的差,如式(4)所示。特殊地,i=k時的相似度s(k,k)稱為參考度,影響聚類的數目,本文將s(k,k)設置為相似度的中位數M(S),可有效完成聚類。最后得到的相似度矩陣S如式(5)所示。
步驟4重復步驟2、步驟3,直到矩陣穩(wěn)定。按最大相似度規(guī)則將其他點分配到距離它最近的聚類中心點相應的聚類中,完成聚類。
(3)對各聚類,計算每個樣本的中點值的平均值,選擇此平均值最大的聚類為理想聚類,計算其算術平均作為最后的BCG模板,以BCG表示,如圖2所示。
2.3.2心跳位置探測
通過將模板在濾波后BCG數據上不斷右滑,形成相關系數函數。對于一段預處理后的BCG信號,采用以下步驟進行模板匹配,確定J點。
(1)通過以下方式構造原BCG信號等長的相關系數函數cor(x)。
(2)對前10 s BCG信號做FFT變換,然后進行頻域尋峰確定心率近似值HRe。在最初的2 s的cor(x)信號中選取最大值,作為第1個J點的位置,以此為起點,向后搜索[60/(HRe+20)s,60/(HRe–20)s]范圍內的最大值點作為下一個J點。當連續(xù)5個所選局部最大值小于閾值θ時,說明匹配質量下降,重新計算模板,本文中閾值θ設為0.75。
最后,根據式(9)和式(10)計算連續(xù)逐拍心動周期和逐拍心率。
其中,IndexJ(n)表示第n個J點的索引位置,Fs為采樣頻率。逐拍心動周期和逐拍心率之間可方便地轉換,為更清楚地比較計算結果,本文使用逐拍心動周期進行誤差分析。
按照如圖3所示方法進行實驗。本文使用壓電陶瓷設備采集BCG信號,同時為驗證BCG逐拍心率計算的準確性,在本研究中同步采集單導聯心電信號進行對比。實驗開始時按單導聯方式在受試者胸腹部佩戴一次性心電電極貼,隨后讓受試者安靜平躺于實驗床墊上,開啟設備進行BCG與ECG的同步采集。通過2.2節(jié)和2.3節(jié)所述方法處理BCG和計算BCG的逐拍心動周期,再使用PT算法[19]定位ECG的R波,計算ECG逐拍心動周期。隨后分析對于2.2節(jié)所述的不同的BCG成分,所提取出的逐拍J-J間期相對于R-R間期的覆蓋率和精準度。得到每次心跳的R-R間期與J-J間期。所有數據分析使用Py Charm與Python3.5進行。
圖4(a)和圖4(b)分別展示了受試者1的一段BCG信號在3~24 Hz帶通濾波和8~24 Hz帶通濾波的預處理方法下,使用本文第2節(jié)所述方法進行逐拍心率提取的結果,圖中使用紅色圓圈標出了所定位的J點的位置,虛線位置為ECG的R波位置。圖5則展示了兩種預處理方式所計算得到的逐拍心動周期與ECG逐拍R-R間期的對比圖??梢娊浲◣?~24 Hz的濾波器預處理之后計算得到的J點定位在原始信號上即為每次心跳所引起的壓電信號最大點位置,但如圖5(a)所示,由此所得到的逐拍心率卻與逐拍R-R間期有較大差別。而經通帶為8~24 Hz的濾波器處理后,所提取的J點位置不再是原始壓電信號每拍心跳的最大值,但卻能得到與R-R間期更為吻合的J-J間期,如圖5(b)所示。我們稱這種現象為“關鍵成分”現象,即8~24 Hz的成分是去除了較低頻的肢體振動等成分、更能反映心室泵血、血液沖擊主動脈引起的心臟原生振動的“關鍵成分”,能夠獲得與R-R間期相近的J-J間期。
圖2 BCG模板聚類與BCG模板
圖3 非接觸式逐拍心動周期檢測準確度驗證實驗示意圖
圖4 不同頻帶濾波的BCG模板匹配計算結果
表1為對受試者1的5 min BCG數據,共計351次心跳統計的結果,心跳覆蓋率定義為使用本文算法成功定位到的心跳(即定位J波與R波相隔50 ms內)的比例,平均絕對誤差定義為逐拍J-J間期與RR間期之差的絕對值的平均值,而平均誤差定義為逐拍J-J間期相對于R-R間期的平均誤差??梢妼τ谑茉囌?的5 min安靜平躺的BCG數據,在分別采用4種預處理方式后,3~24 Hz濾波以及db6小波重構都取得了100%的心跳檢出率,但卻有較高的平均絕對誤差。經8~24 Hz濾波和二次差分處理后,雖心跳檢出率略有不足,但卻取得了很低的平均絕對誤差。
圖6展示了兩種不同預處理方式下受試者1的逐拍R-R間期與J-J間期Bland-Altman圖,從中可以看出采用寬帶寬濾波器所計算的結果(圖6(a))在所有心跳間期分布范圍內誤差都明顯高于窄帶濾波器計算結果(圖6(b))。兩種預處理方式的95%置信區(qū)間分別為(–32.26 ms,32.26 ms)和(–12.34 ms,12.34 ms),且兩種處理方式所得逐拍心跳間期誤差均值都為0,說明在給定時間段內使用BCG和ECG檢測出了相同的心跳次數。此處需要說明的是,雖表1中8~24 Hz濾波器的處理方法心跳覆蓋率低于100%,但這僅說明BCG有心跳檢測位置與R波位置相差較多,總體依然是檢測到了應有的心跳次數。
本文采集了15名受試者的BCG和ECG數據,年齡23~34歲,身高156~182 cm,體重45~90 kg。表2展示了所有受試者各5 min安靜平躺狀態(tài)下采集到的數據的處理結果,并且展示了4種BCG成分對逐拍心率提取性能的比較。
由表2可見對于8~24 Hz的帶通濾波和二次差分這兩種方式處理,在檢出心跳覆蓋率方面略有不足,這是因為頻率較高時波形表現得較為復雜,影響了模板匹配的效果;而同時這兩種方式在平均絕對誤差(MAE)方面有著較大的優(yōu)勢,這是因為高頻成分更能反映心臟跳動的成分而不受整個身體震動的影響。
圖5 ECG-RR間期與不同頻率濾波的逐拍BCG-JJ間期折線圖
表1 受試者1的數據處理結果對比
圖7展示了2.1節(jié)所述各成分所計算得到的J-J間期與R-R間期誤差絕對值的分布,可以看到對于前兩種處理方式,75%的誤差分別在20 ms和16 ms以下,而后兩種處理方式的誤差有75%在4 ms以下,且中位數、下四分位數、下限均為0 ms,說明有1/2以上的J-J間期具有與R-R間期相同的采樣點數。
圖6 受試者1的逐拍RR間期與JJ間期Bland-Altman圖
表2 不同預處理方式下逐拍心率提取性能
圖7 各成分誤差分布箱型圖
此外,本文還從系統整體的角度對比了同類研究工作在靜息狀態(tài)下測量BCG提取逐拍心率的主要性能指標,如表3所示。可見,本文方法能夠在覆蓋率和逐拍心動周期的準確度方面處于領先水平,尤其是本文MAE可達到4 ms以下,這有助于提高BCG計算HRV的可靠性。
本文提出了一種基于壓電陶瓷傳感器采集BCG并提取精準逐拍心率的方法,該方法是通過優(yōu)化傳感器性能提高BCG信號質量,通過對比BCG信號的不同處理方式,找到最適合提取精準逐拍心率的成分,最后通過AP聚類和模板匹配的方法確定J點,計算逐拍心率。
表3 BCG逐拍心率提取方法性能對比
本文的逐拍心動周期平均絕對誤差為3.78 ms,心跳檢測覆蓋率為98.5%。本研究結果表明,通過設計合適的BCG采集裝置、采用合適的信號處理方法處理BCG信號以及選擇合適的特征點提取算法,可以使BCG提取逐拍心動周期(J-J間期)相對于R-R間期的誤差降至可接受的范圍,這對于提高使用BCG進行健康監(jiān)測的效果尤其是提高HRV計算的可靠性具有重要意義。對本文所提出的BCG“關鍵成分”成因進行實驗仿真和驗證,以及對BCG提取逐拍心率精度的影響因素作進一步量化分析將是下一步的研究方向。