基于心電信號提取呼吸信號的算法研究

蔣蓮，陳兆學(xué)

上海理工大學(xué)醫(yī)療器械與食品學(xué)院，上海200093

前 言

呼吸率是臨床醫(yī)學(xué)中常用的監(jiān)護(hù)參數(shù)之一。目前應(yīng)用于人體呼吸信號的檢測方法主要分為兩類：一種方法是利用溫度、濕度或氣流傳感器直接檢測法，通過檢測從人體鼻部流入和流出肺部的氣體得到人體呼吸波；另一種方法是利用阻抗式或壓力傳感器的間接檢測法，通過檢測由人體呼吸引起的人體胸部形變或胸部阻抗變化而獲得呼吸波［1］。這兩種方法都有其缺點(diǎn)，前一種方法受檢位置在鼻部，會使受檢者在檢測時感到呼吸不適；后一種方法需要受檢者在安靜時進(jìn)行測量，不能夠長時間動態(tài)檢測人體呼吸。諸多的臨床資料證實(shí)，當(dāng)人體呼吸時會導(dǎo)致心電圖的波形發(fā)生變化。即心電信號波形的幅值變化與人體呼吸運(yùn)動存在一定的聯(lián)系，人體心電信號中包含呼吸信息［2］，當(dāng)人體呼吸時，胸腔部位節(jié)律性的收縮和擴(kuò)張運(yùn)動將使得描述心臟動作電位傳播方向的心電軸也出現(xiàn)周期性的變化，從而導(dǎo)致心電信號波形的QRS 波呈現(xiàn)為周期性的波動［3］。在臨床應(yīng)用中，為了減少含有各種傳感器的其它檢測設(shè)備的使用和減少受檢者的束縛，從人體心電信號中直接解析出呼吸信號信息顯得尤為重要，它不僅使動態(tài)檢測人體呼吸率成為可能，還是一種高效、低成本的呼吸率檢測方法。

許多研究者基于對心電信號的處理和分析已經(jīng)獲得了呼吸信號。Wang等［4］提出基于心電波形變化的EDR（ECG-Derived Respiratory Signal）算法，基于心電向量圖（Vector Cardiogram,VCG）提供了人體呼吸率和呼吸深度。Pinciroli等［5］和Moody等［6］提出基于心電電軸的分析法。Pilgram等［7］提出基于心率的EDR 方法，該方法需要對信號的R-R 間隔進(jìn)行奇異值分解（Singular Value Decomposition,SVD），從而得到呼吸信號的頻率。Ding 等［8］提出一種心電信號高階統(tǒng)計(jì)量的計(jì)算方法，該方法需要檢測出心電信號的R波，通過計(jì)算信號的高階統(tǒng)計(jì)量來獲取人體呼吸信息。Orphanidou 等［9］選用數(shù)據(jù)融合方法來獲得呼吸信息，將心電信號R波峰值的振幅調(diào)制和呼吸性竇性心律失常兩種方法融合在一起，從而達(dá)到了良好的呼吸信號提取效果。從心電信號中提取呼吸信號的算法雖然很多，但是各種處理算法存在優(yōu)劣，有的算法計(jì)算速度快，但準(zhǔn)確率并不高，有的算法計(jì)算準(zhǔn)確率高，但受到心電信號數(shù)據(jù)量和計(jì)算復(fù)雜度的限制。因此，本文提出一種獨(dú)立分量分析（ICA）的EDR方法。通過MATLAB 對算法進(jìn)行驗(yàn)證，并將實(shí)驗(yàn)結(jié)果與其它研究結(jié)論相比較，經(jīng)研究發(fā)現(xiàn)該算法能夠從各種干擾中直接分離出呼吸信號，避免各種復(fù)雜濾波器的設(shè)計(jì)，且計(jì)算速度快，結(jié)果準(zhǔn)確，是一種有效和可靠的EDR算法。

1 EDR算法

人體心電信號是非常微弱的信號，容易受到工頻噪聲干擾和呼吸運(yùn)動的影響。因此，通過醫(yī)用設(shè)備檢測到的人體心電信號往往混有很強(qiáng)的背景噪聲，最主要的3種干擾有：工頻噪聲干擾、肌電干擾和基線漂移［10］。其中基線漂移由于人體呼吸、胃腸蠕動以及導(dǎo)聯(lián)電極的移動等低頻干擾產(chǎn)生，表現(xiàn)為信號波形會隨著基線出現(xiàn)漂移［11］，它的頻率很低，一般在1 Hz以下。而人體呼吸信號的頻率一般為0.1～0.5 Hz，也分布在心電信號的低頻部分。因此直接通過時頻域方法提取呼吸信號會存在較大困難。ICA 方法是近年來由盲源分離技術(shù)發(fā)展而來的一種多維度信號統(tǒng)計(jì)處理方法，可以根據(jù)源信號的基本統(tǒng)計(jì)特征，由觀測數(shù)據(jù)最終恢復(fù)出源信號［12］。因此，基于ICA方法尋找源信號的思想，關(guān)鍵部分是需要從原信號中提取出包含呼吸信息的混合矩陣，再運(yùn)用ICA算法即可分離出呼吸信號。通過ICA方法的EDR算法流程如圖1所示。

圖1 EDR算法的流程圖Fig.1 Flow chart of electrocardiogram-derived respiratory signal(EDR)algorithm

2 EDR算法的實(shí)現(xiàn)

本研究使用的是由動態(tài)心電圖［13］檢測到的不同人體的24 h 心電數(shù)據(jù)，采樣頻率為250 Hz。仿真驗(yàn)證時，將原始心電數(shù)據(jù)分成12 段，取每一段的1 min時間長度，共計(jì)15 000 個采樣點(diǎn)的心電數(shù)據(jù)信息進(jìn)行顯示。為方便說明，本文僅以其中一段數(shù)據(jù)的1 min數(shù)據(jù)為例。

2.1 QRS波群檢測

快速而準(zhǔn)確檢測出心電信號中QRS 波群的位置，是心電信號自動分析和診斷系統(tǒng)中的一個關(guān)鍵環(huán)節(jié)，各種算法層出不窮［14-16］。呼吸信號隱含在人體心電信號中，具體體現(xiàn)在：人體呼吸運(yùn)動會引起心電信號特征點(diǎn)R波的幅值變化，還對R-R波間期存在一定的影響?；谶@樣的理論基礎(chǔ)，需要從心電信號中提取出包含呼吸信號的R 波特征點(diǎn)和S 波特征點(diǎn)。本文采用Pan＆Tompkins 算法檢測信號的QRS波群，該方法具有原理簡單、計(jì)算量較小和實(shí)時性好等優(yōu)點(diǎn)［17］。其數(shù)據(jù)檢測結(jié)果如圖2所示。

圖2 經(jīng)Pan＆Tompkins算法檢測到的R波和S波特征點(diǎn)Fig.2 Detected feature points of R wave and S wave by Pan＆Tompkins algorithm

2.2 R波和S波特征點(diǎn)的三次樣條函數(shù)插值和重采樣處理

由于獲得的R 波和S 波峰值及其索引值均為離散數(shù)據(jù)點(diǎn)，若要得到連續(xù)的包含R波和S波信息的數(shù)據(jù)序列，則需要對離散點(diǎn)進(jìn)行插值。本文中選用了三次樣條函數(shù)插值法［18］。三次樣條函數(shù)插值法是基于一個特定的數(shù)據(jù)集合，如果給出n個給定點(diǎn)的數(shù)據(jù)，就可以利用這n個給定點(diǎn)將數(shù)據(jù)集合分為n-1段，用n-1 段三次多項(xiàng)式即可在每兩個連續(xù)給定數(shù)據(jù)點(diǎn)之間構(gòu)建一個三次樣條［19］。三次樣條插值算法的流程如下，設(shè)共有n個插值節(jié)點(diǎn)x1＜x2＜…＜xn-1＜xn，則經(jīng)過數(shù)據(jù)點(diǎn)(x1,y1),(x2,y2),…(xn-1,yn-1),(xn,yn)的三次樣條S(x)是一組三次多項(xiàng)式：

該多項(xiàng)式在節(jié)點(diǎn)處有連續(xù)性和在節(jié)點(diǎn)處的一階導(dǎo)和二階導(dǎo)有光滑性。根據(jù)邊界條件可以求得ai、bi、ci以及di的值，從而求得Si(x)的表達(dá)式，進(jìn)而得到擬合后的光滑曲線。經(jīng)過插值擬合的兩路信號如圖3所示。

圖3 經(jīng)樣條插值算法擬合的R波和S波序列Fig.3 R and S wave sequences fitted by spline interpolation

從心電信號的一個周期波形中可以看到，R波出現(xiàn)在S波的前面，故經(jīng)檢測到的R波和S波特征點(diǎn)在位置上會存在一定的延時。為了保持?jǐn)M合后的R 波和S波序列的一致性，還需要對該兩路序列進(jìn)行重采樣，從而獲得在同一時刻點(diǎn)采樣的兩路信號序列。經(jīng)過重采樣的兩路信號序列如圖4所示。

圖4 在同一位置重采樣得到的R波和S波序列Fig.4 R and S wave sequences re-sampled at the same positions

2.3 應(yīng)用小波理論重構(gòu)一路呼吸信號

小波變換即是通過小波技術(shù)對原信號進(jìn)行分解重構(gòu)，主要包括空間（時間）和頻率的局部變換，應(yīng)用平移和伸縮運(yùn)動等運(yùn)算功能可實(shí)現(xiàn)信號或函數(shù)的多尺度的細(xì)化分析。當(dāng)對信號進(jìn)行一階小波分解后，信號被解為兩部分：近似信號A和細(xì)節(jié)信號D。其中，近似信號包含信號的低頻成分，細(xì)節(jié)信號包含信號的高頻成分［20］。當(dāng)運(yùn)用小波理論重構(gòu)呼吸信號時，首先需要確定信號的分解層數(shù)。相關(guān)研究表明，在基于小波變換的心電信號頻率分解和重建中，coif4小波基函數(shù)是心電信號數(shù)字濾波的最佳選擇［21］，其中，小波coif4函數(shù)的頻率分解規(guī)則如表1所示。

本文使用的心電信號頻率為250 Hz，而正常人體每小時的呼吸速率約20 次/min，換算成頻率為0.2～0.5 Hz。參照表1中coif4 小波函數(shù)的頻率分解規(guī)則，則需要進(jìn)行9層頻率分解才可得到呼吸信息所在頻段。為了從心電信號中重構(gòu)出一路呼吸信號，需要原信號頻率作十次分解，從而得到十個細(xì)節(jié)分量值，通過保留第九層的小波分量，并把分解后信號的其它層次細(xì)節(jié)分量均置為零，再采用waverec小波函數(shù)對第九層細(xì)節(jié)分量進(jìn)行重構(gòu)就可得到一路呼吸信號序列，重構(gòu)的呼吸信號序列如圖5所示。

2.4 ICA算法理論

表1 coif4函數(shù)的頻率分解規(guī)則Tab.1 Frequency decomposition rule of coif4 function

圖5 重構(gòu)的呼吸信號序列Fig.5 Reconstructed respiratory signal sequence

ICA是盲源分離領(lǐng)域的一項(xiàng)關(guān)鍵技術(shù)，它在沒有任何先驗(yàn)知識的前提下就能從多個混合信號中提取或分離出源信號［22］。ICA問題的數(shù)學(xué)表達(dá)式為：

式中，s是n維源信號向量，s=[s1,s2,…sn] ；A是混合矩陣；x是m維觀測的混合信號向量，x=[x1,x2,…xm]T；觀測信號x為信源向量s的瞬時線性組合。

ICA 算法的實(shí)質(zhì)是依據(jù)觀測信號x和源信號s分量之間的統(tǒng)計(jì)獨(dú)立性假設(shè)，根據(jù)源信號s在概率分布上的某些先驗(yàn)條件，最終估算出一個混合矩陣A。即通過求解出一個分離矩陣W，使得y=Wx的各分量盡可能保持統(tǒng)計(jì)獨(dú)立，并把y看作為源信號s的估計(jì)。ICA的原理框圖如圖6所示。

圖6 ICA的原理框圖Fig.6 Principle diagram of independent component analysis

運(yùn)用ICA模型進(jìn)行信號分析時，需滿足以下假設(shè)條件：（1）源信號s(t)各分量之間是統(tǒng)計(jì)獨(dú)立的，即分量間是線性無關(guān)的，并且所有分量中最多一個信號滿足高斯分布；（2）觀測的混合信號個數(shù)要大于或等于源信號個數(shù)，即m ＞n；（3）混合矩陣A是列滿秩的可逆矩陣。

目前，常用的ICA 算法大致可分為兩大類：一類是基于隨機(jī)梯度法的自適應(yīng)算法，該算法可以通過神經(jīng)網(wǎng)絡(luò)實(shí)現(xiàn)，但是受選取的學(xué)習(xí)速率參數(shù)的影響，其運(yùn)算收斂速度慢。另一類是基于批運(yùn)算相關(guān)目標(biāo)函數(shù)的快速算法，即FastICA［23］，該算法主要是尋找一個方向使得WTx具有最大高斯性，主要采用了定點(diǎn)迭代的思想，其計(jì)算簡單且具有較好的收斂性。因此，本文中選用了基于負(fù)熵的FastICA算法獲取源信號序列。將經(jīng)過插值和重采樣分別得到的R 波序列和S波序列，利用小波函數(shù)分解和重構(gòu)的一路呼吸信號以及現(xiàn)有的連續(xù)24 h 心電原始信號序列4 路信號混合，構(gòu)成一個觀測矩陣。選用基于負(fù)熵的FastICA 算法分離得到不同的獨(dú)立源分量序列。經(jīng)過ICA得到的3路源信號序列如圖7所示。

通過FastICA 算法，將4 路信號分離后共得到3路獨(dú)立的源信號序列，分別為源信號Z1序列、源信號Z2序列和源信號Z3序列。觀察源序列波形圖，可以看出前兩路源信號與重構(gòu)的呼吸信號有著良好的相關(guān)性，而源信號Z3 序列則保留了更多的基線漂移信號特征。由于本章主要需要獲得人體呼吸率，故將對前兩路源信號序列的結(jié)果著重進(jìn)行對比分析。

3 EDR算法驗(yàn)證

由于人體呼吸率的頻率低，頻譜中存在較大其他頻譜干擾，使得參數(shù)檢測誤差較大。因而采用時域的方法檢測該參數(shù)更為有利。在時域內(nèi)，采用下采樣和滑動平均的方法對提取出的源信號序列進(jìn)行濾波處理，得到更為準(zhǔn)確的代表呼吸信號的序列［24］，再采用峰值查找法標(biāo)記出各序列的峰值。將基于FastICA算法得到的Z1序列和Z2序列峰值與基于R波幅度調(diào)制和基于S波幅度調(diào)制提取的呼吸波形峰值相比較，各信號序列的峰值標(biāo)記波形圖如圖8所示。

由圖8可以得知，應(yīng)用FastICA算法得出的兩路源信號Z1和Z2序列與基于R波或S波幅度調(diào)制方法提出的呼吸信號波形均有著良好的相關(guān)性，將該4路信號的幅值進(jìn)行歸一化，得到4路信號的波形如圖9所示。

圖7 ICA得到的3路源信號序列Fig.7 Three source signal sequences obtained with independent component analysis

圖8 各信號序列的峰值標(biāo)記波形圖Fig.8 Waveforms with labeled peak signals of each signal sequences

圖9 歸一化幅值的4路信號波形序列Fig.9 Four signal waveform sequences with normalized amplitude

統(tǒng)計(jì)數(shù)學(xué)上，常用內(nèi)積來反映各數(shù)據(jù)間的關(guān)聯(lián)性。經(jīng)計(jì)算得到4 路信號的內(nèi)積值如表2所示。內(nèi)積值越接近1，表明兩數(shù)據(jù)的關(guān)聯(lián)性越大。內(nèi)積值越接近0，說明兩列數(shù)據(jù)不相關(guān)。從內(nèi)積表中，我們可以得出經(jīng)過ICA 得到的Z1 序列和Z2 序列均可以提取出呼吸信號。因此，將Z1序列和Z2序列二者數(shù)據(jù)取平均得到新的代表呼吸信號的序列。

通過滑動平均法檢測序列的峰值，各序列統(tǒng)計(jì)得到的峰值數(shù)值即反映了人體每分鐘呼吸次數(shù)。將現(xiàn)有的24 h 連續(xù)信號以每2 h 為間隔，共分為12 段，對每段內(nèi)的1 min時長信號序列進(jìn)行峰值統(tǒng)計(jì)，共累積統(tǒng)計(jì)120 min 內(nèi)的呼吸總數(shù)，最后，通過求平均得到的每兩小時內(nèi)呼吸率結(jié)果如表3所示。

4 結(jié)論

為了減少額外醫(yī)療設(shè)備的使用，降低受測者的不適感，本文提出一種EDR 算法。首先選用幾路相互獨(dú)立的涵有呼吸信號的數(shù)據(jù)序列，其次運(yùn)用ICA分析方法分離出源信號序列，然后將得到的源信號序列與其他方法提取的呼吸信號序列進(jìn)行比較，最后確定以ICA 分離出來的兩路信號的均值作為提取的呼吸信號序列。從驗(yàn)證的結(jié)果來看，每分鐘內(nèi)呼吸率檢測結(jié)果在人體正常情況下每分鐘呼吸結(jié)果范圍內(nèi)，證明了上述各種算法的正確性，雖然在時域內(nèi)峰值的檢測存在一定的誤差，但是通過該方法提取出的呼吸信號波形與其他方法提取出的呼吸信號存在良好的相關(guān)性，其平均相似度達(dá)到95.94%以上，說明該方法提取呼吸信號的有效性。并且，使用該算法避免了各種復(fù)雜濾波器的設(shè)計(jì)，算法實(shí)現(xiàn)簡單，滿足了一般病人呼吸參數(shù)檢測的需求，具有相當(dāng)?shù)膶?shí)用價值和意義。

表2 4路信號序列的夾角余弦值Tab.2 Cosine coefficient of 4 signal sequences

表3 EDR算法仿真結(jié)果比較Tab.3 Simulation results of EDR algorithm