基于EMD-ICA的心沖擊信號(hào)降噪研究

姜 星 耿讀艷，* 張園園 付志剛

1(河北工業(yè)大學(xué)省部共建電工裝備可靠性與智能化國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300130)2(河北工業(yè)大學(xué)河北省電磁場(chǎng)與電器可靠性重點(diǎn)實(shí)驗(yàn)室，天津 300130)3(中國(guó)人民解放軍第254醫(yī)院，天津 300142)

引言

心沖擊圖(ballistocardiogram，BCG)可以展示心臟搏動(dòng)和血液在大動(dòng)脈流動(dòng)而引起的人體對(duì)外壓力或體表位移的變化，反映心臟的力學(xué)特征[1]，通過(guò)壓力、加速度、電容等傳感器來(lái)感知身體的這種微小位移變化，從而獲取心臟搏動(dòng)引起的體震信息。正常的BCG信號(hào)波形及其生理意義如圖1所示，主要包含H、I、J、K、L、M和N等7種波形，其中I、J、K反映心動(dòng)周期的射血過(guò)程，而L、M、N反映舒張過(guò)程[2]。相比Holter等心電監(jiān)測(cè)方式，BCG可實(shí)現(xiàn)更長(zhǎng)期的連續(xù)監(jiān)測(cè)，避免電極及其導(dǎo)線束縛造成對(duì)人情緒的影響，尤其適合于長(zhǎng)時(shí)間作息和睡眠監(jiān)測(cè)，比較適用于家庭應(yīng)用。然而，BCG信號(hào)自身十分微弱，且易受到呼吸、體動(dòng)、工頻噪聲的干擾[3-4]，導(dǎo)致直接測(cè)量得到的BCG信號(hào)含有太多噪聲，不能直接獲得準(zhǔn)確的生理信息特征。因此，BCG的降噪處理顯得尤為必要。

圖1 正常BCG信號(hào)波形及其生理意義Fig.1 Normal BCG signal and physiological significance

目前對(duì)于BCG信號(hào)的降噪處理主要是采用小波變換方法，一般是對(duì)周期性平穩(wěn)信號(hào)有效，而對(duì)包含尖峰或突變的非平穩(wěn)信號(hào)效果不佳[5-6]。小波降噪能同時(shí)進(jìn)行時(shí)頻分析，且能有效區(qū)分非平穩(wěn)信號(hào)的尖峰與噪聲，但需選取合適的小波基才能達(dá)到較好降噪效果，而小波降噪仍基于Fourier變換，尚不能擺脫以Fourier變換為基礎(chǔ)帶來(lái)的缺陷[7-9]。經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, EMD)是Huang在1998年提出的一種用于非線性和非平穩(wěn)時(shí)間序列信號(hào)的處理方法[10]，主要是將信號(hào)分解為若干個(gè)固有模態(tài)函數(shù)(intrinsic mode function, IMF)分量，使每個(gè)IMF包含不同的頻率成分，從而實(shí)現(xiàn)信號(hào)特征的有效分解。但是在分解過(guò)程中，由于BCG信號(hào)中有大量噪聲成分，存在較強(qiáng)的時(shí)頻耦合特征，會(huì)產(chǎn)生模態(tài)混疊現(xiàn)象[11]，造成IMF成分中包含不同的時(shí)間尺度，從而對(duì)各個(gè)模態(tài)分量的物理意義判別造成困難。獨(dú)立分量分析(independent component analysis, ICA)無(wú)需信號(hào)源先驗(yàn)知識(shí)，在時(shí)域和頻域都能將觀測(cè)得到的混合信號(hào)中的獨(dú)立成分分離出來(lái)。ICA去噪的關(guān)鍵點(diǎn)在于噪聲通道的構(gòu)建，孫云蓮等利用EMD分解得到的IMF希爾伯特時(shí)頻譜來(lái)構(gòu)建噪音通道，在信噪比較高的情況下，基于IMF時(shí)頻譜分析的ICA降噪法取得了不錯(cuò)的效果[12]。本研究在其基礎(chǔ)上，提出一種基于EMD及ICA相結(jié)合的BCG信號(hào)降噪方法。該方法對(duì)含噪BCG信號(hào)進(jìn)行EMD分解，得到頻率從高到低的多個(gè)IMF；通過(guò)引入模態(tài)相關(guān)法，找到噪聲與信號(hào)的分界IMF；對(duì)分界之前的分量構(gòu)建虛擬通道[13]，與原信號(hào)一起作為ICA的輸入，得到降噪后BCG信號(hào)。

1 方法

1.1 信號(hào)采集系統(tǒng)

在本研究中，使用的是自主研制的BCG信號(hào)采集系統(tǒng)，主要由采集人體信號(hào)的壓電薄膜式傳感器、連接傳感器對(duì)采集信號(hào)做后續(xù)處理的電路和上位機(jī)軟件組成。將傳感器置于座椅坐墊下，受試者靜止地坐在座椅上，心臟射血時(shí)引起身體微小振動(dòng)，通過(guò)傳感器獲取射血引起的重力變化信號(hào)[14]，經(jīng)過(guò)電路的多級(jí)放大和轉(zhuǎn)換，將其無(wú)線傳輸?shù)缴衔粰C(jī)軟件，最終得到測(cè)量的BCG信號(hào)。采集電路主要包括電荷放大器、低通濾波器、電壓放大器、電壓抬升電路、信號(hào)離散模塊、無(wú)線傳輸模塊和串口通信，其總體構(gòu)成如圖2所示。

圖2 電路整體設(shè)計(jì)框圖Fig.2 Overall design circuit diagram

1.2 EMD-ICA原理

經(jīng)驗(yàn)?zāi)B(tài)分解的實(shí)質(zhì)是對(duì)一個(gè)時(shí)間序列信號(hào)進(jìn)行平穩(wěn)化處理，其結(jié)果是將信號(hào)中不同尺度的波動(dòng)或趨勢(shì)逐級(jí)分解，產(chǎn)生一系列具有不同特征尺度的數(shù)據(jù)序列，每一個(gè)序列稱為IMF。分解得到的IMF必須滿足兩個(gè)條件：一是IMF中極值點(diǎn)的個(gè)數(shù)與過(guò)零點(diǎn)的個(gè)數(shù)相等或不超過(guò)1個(gè)；二是由極大值與極小值確定的包絡(luò)線均值為0。

獨(dú)立分量分析是由盲源分離計(jì)算發(fā)展起來(lái)的一種信號(hào)處理方法，要求輸入的是多維信號(hào)，并且混合信號(hào)個(gè)數(shù)要大于等于源信號(hào)個(gè)數(shù)，才能有效地分離。

由獨(dú)立分量分析的原理可知，獨(dú)立分量分析要求輸入的是多維信號(hào)，并且混合信號(hào)個(gè)數(shù)要大于等于源信號(hào)個(gè)數(shù)，才能有效地分離?；诖?，BCG信號(hào)經(jīng)EMD分解后產(chǎn)生若干個(gè)基本信號(hào)，即

(1)

式中，x(t)表示原始信號(hào)，imfi表示第i個(gè)IMF分量，rn(t)表示余項(xiàng)。

提取出噪聲所在的IMF分量進(jìn)行組合，得到噪聲的估計(jì)量；與源信號(hào)一起作為ICA多維噪聲輸入，構(gòu)建ICA噪聲通道。本研究采用基于EMD模態(tài)相關(guān)分選準(zhǔn)則[15]判定噪聲與信號(hào)的分界層，計(jì)算模態(tài)函數(shù)與原始BCG信號(hào)之間的互相關(guān)系數(shù)，有

(2)

圖3 EMD-ICA降噪原理Fig.3 EMD-ICA de-noising principle

1.3 信號(hào)采集過(guò)程

本研究為了驗(yàn)證所得到的BCG信號(hào)的準(zhǔn)確性，使用ECG-300 A數(shù)字式三道心電圖機(jī)和所設(shè)計(jì)的BCG采集系統(tǒng)，同步采集被試者的ECG和BCG信號(hào)，并實(shí)時(shí)保存。對(duì)所采集的BCG信號(hào)分析得到心率，并與心電圖機(jī)所得心率進(jìn)行對(duì)比，分析BCG的準(zhǔn)確性。

為了驗(yàn)證裝置的有效性，并對(duì)人體心沖擊信號(hào)的特征進(jìn)行研究，本研究共采集了10位健康成年人的數(shù)據(jù)，其中有5位男性、5位女性，年齡均在23～27歲。要求受試者保持平穩(wěn)的坐姿，坐在放有壓電薄膜傳感器的座椅上，勻速呼吸，將傳感器的輸出端接在信號(hào)調(diào)理電路上，使用示波器觀察其BCG信號(hào)，等待信號(hào)穩(wěn)定后采集2 min的BCG信號(hào)，同時(shí)采集被試者的ECG信號(hào)。將采集到的數(shù)字信號(hào)經(jīng)過(guò)Matlab處理后得BCG的原始信號(hào)，圖4為某一被試者采集的原始BCG信號(hào)。

圖4 原始BCG信號(hào)Fig.4 Example of the original BCG signal

由圖4可見(jiàn)，BCG信號(hào)中可以明顯觀察到J峰，同時(shí)可以看到BCG信號(hào)的I、J、K波形是有重復(fù)性的，所以表示所采集的BCG信號(hào)是具有魯棒性的。另外，不同的被試者BCG信號(hào)的波形形狀和幅度存在著較大的差別，這些差別來(lái)源于每個(gè)人不同的心臟活動(dòng)狀態(tài)以及身體結(jié)構(gòu)。

1.4 信號(hào)處理

將采集到的BCG信號(hào)進(jìn)行EMD分解，結(jié)果如圖5所示。

圖5 EMD分解結(jié)果Fig.5 EMD decomposition results

由分解結(jié)果可以看出，原始信號(hào)經(jīng)EMD分解得到了幾個(gè)固有模態(tài)分量和一個(gè)趨勢(shì)項(xiàng)。由EMD分解過(guò)程可知，分解得到的13個(gè)IMF分量imf1，imf2，…，imf13中，噪聲對(duì)每個(gè)IMF分量的支配作用逐漸降低，信號(hào)對(duì)IMF的支配作用不斷加強(qiáng)，因此需要確定IMF分量中信號(hào)與噪聲的分界[16]。根據(jù)式(1)計(jì)算出各階固有模態(tài)分量與原始BCG信號(hào)的相關(guān)系數(shù)，得到結(jié)果后繪制相應(yīng)的曲線，如圖6所示。

圖6 相關(guān)系數(shù)曲線Fig.6 The correlation coefficient curve table

按照前述信號(hào)與噪聲的分界點(diǎn)判定方法，由圖6可以看出，imf4為分界IMF。根據(jù)噪聲與信號(hào)的分界，選擇第1～3個(gè)IMF分量作為虛擬噪聲通道。由于原始BCG信號(hào)進(jìn)行EMD分解后imf4之前的IMF分量大都是噪聲成分，所以借鑒信號(hào)時(shí)序平移思想構(gòu)造ICA的多維輸入[17]，其基本思想是：把原始BCG信號(hào)分解后得到的由第1～3個(gè)IMF分量構(gòu)建的虛擬噪聲通道時(shí)序序列向左平移l個(gè)位置，向左溢出部分拼接到該時(shí)序序列的右端，這樣就可以得到新的含噪信號(hào)作為ICA的噪聲通道輸入。這個(gè)新的含噪信號(hào)與原信號(hào)相比信噪比幾乎一致，不僅保留了原信號(hào)中的有效成分，而且信號(hào)中含噪聲的成分也沒(méi)有大的改變，僅是噪聲形態(tài)發(fā)生了改變。

這樣就構(gòu)造了多個(gè)輸入信號(hào)，包括兩個(gè)虛擬噪聲通道和原始BCG信號(hào)。在獨(dú)立分量分析中，常用的算法是Fast ICA和Infomax算法，從計(jì)算速度方面考慮，本研究采用一種基于負(fù)熵的判據(jù)和一種非常高效的FastI CA[18-19]算法進(jìn)行分析，其框圖如圖7所示。負(fù)熵定義如下：

(3)

式中，E為均值運(yùn)算，g為非線性函數(shù)，Ygauss是高斯變量。Fast ICA其步驟如下：

1) 對(duì)觀測(cè)信號(hào)進(jìn)行中心化處理，使其均值為0。

2) 利用特征值分解法對(duì)中心化后的數(shù)據(jù)進(jìn)行白化，得到標(biāo)準(zhǔn)化數(shù)據(jù)Z(t)。

3) 選擇要估計(jì)的獨(dú)立成分個(gè)數(shù)m，迭代次數(shù)為p，初始值為p=1。

4) 選擇一個(gè)隨機(jī)的初始分離矩陣W。

5) 迭代更新分離矩陣，令

(4)

(5)

Wp=Wp/‖Wp‖

(6)

式(4)中,g(y)=tanh(a1,y),a1=1。

6) 判斷Wp是否收斂，如果不收斂，返回步驟5)。

圖7 Fast ICA算法框圖Fig.7 Fast ICA algorithm block diagram

7) 令p=p+1，假如p≤m，返回步驟4)，重復(fù)執(zhí)行步驟4)～7)，直到將所有需要的獨(dú)立分量提取出來(lái)為止，算法結(jié)束。

1.5 算法性能評(píng)價(jià)

為了定量說(shuō)明本方法的降噪效果，同時(shí)采用EMD方法和小波去噪方法[20]，分別對(duì)本研究的含噪BCG信號(hào)進(jìn)行降噪處理，比較不同方法的降噪效果。為了定量分析3種方法的降噪效果，計(jì)算其降噪后信號(hào)與原始信號(hào)的相關(guān)系數(shù)(R值)、信噪比(Rsn值)和能量比(Esn值)。相關(guān)系數(shù)值說(shuō)明了降噪后的信號(hào)與原始信號(hào)之間整體波形的相似度，信噪比說(shuō)明了信號(hào)的去噪情況，而能量比值反映了降噪后信號(hào)占原信號(hào)的能量百分比。其中，信噪比Rsn定義[21]為

(7)

若降噪后信號(hào)具有較高的信噪比，說(shuō)明降噪效果好。

信號(hào)能量計(jì)算公式[16]如下：

(8)

式中，x(t)為信號(hào)幅值，T為采樣周期，m為采樣點(diǎn)數(shù)。

降噪后的信號(hào)占原信號(hào)的能量百分比Esn定義為

Esn=E0/E

(9)

式中，E表示原始信號(hào)能量，E0表示降噪后信號(hào)能量。

Esn越大，說(shuō)明降噪后的信號(hào)越高，越能保持原信號(hào)特征，也越接近原信號(hào)。

2 結(jié)果

2.1 采集系統(tǒng)

為了驗(yàn)證信號(hào)采集系統(tǒng)的準(zhǔn)確性，經(jīng)過(guò)心沖擊圖和心電圖方法計(jì)算得出心率[22]，如表1所示。對(duì)其用配對(duì)t檢驗(yàn)方法進(jìn)行檢驗(yàn)，結(jié)果顯示兩種方法計(jì)算得出的心率不具有顯著性差異(P>0.05)，說(shuō)明本研究的BCG采集系統(tǒng)具有一定的準(zhǔn)確性。

表1 被試者的心率計(jì)算結(jié)果

圖8 BCG信號(hào)降噪前后對(duì)比。(a)降噪前BCG信號(hào)；(b)降噪后BCG信號(hào)Fig.8 BCG signal de-noising before and after the comparison. (a) Noise reduction before BCG signal；(b) Noise reduction after BCG signal

2.2 降噪效果

根據(jù)上述的分析方法，由Fast ICA最終得到降噪后的BCG信號(hào)。圖8所示為BCG信號(hào)降噪前后的波形，可以觀察到降噪處理后BCG的信號(hào)特征得到了明顯的保留，信號(hào)波形較之前也趨于平滑。為了進(jìn)一步證明降噪效果，計(jì)算了降噪前后BCG信號(hào)的功率譜密度(power spectral density，PSD)，如圖9所示。可以看出，在高頻部分BCG信號(hào)功率譜明顯趨于平滑，且在低頻部分沒(méi)有明顯的衰減，進(jìn)一步說(shuō)明本研究采用的方法在去除噪聲的同時(shí)，很好地保留了BCG信號(hào)的低頻特性。

利用算法評(píng)價(jià)指標(biāo)，分別計(jì)算10名被試的原始BCG信號(hào)降噪前后的相關(guān)系數(shù)、信噪比、能量百分比，得到的結(jié)果如圖10所示。對(duì)其進(jìn)行單因素ANOVA檢驗(yàn)，結(jié)果顯示3種降噪方法計(jì)算得出的相關(guān)系數(shù)、信噪比、能量比具有顯著性差異(P<0.05)；并計(jì)算其平均值，結(jié)果如表2所示。

表2 不同降噪算法效果對(duì)比Tab.2 Comparison of different noise reduction algorithms

圖9 降噪前后BCG信號(hào)的PSD。(a)降噪前；(b)降噪后Fig.9 PSD before and after de-noising BCG signal. (a) Before de-noising；(b) After de-noising

圖10 不同降噪方法處理前后的對(duì)比。(a)相關(guān)系數(shù)對(duì)比；(b)信噪比對(duì)比；(c)能量比對(duì)比Fig.10 Comparison of quantitative figure before and after different noise reduction methods. (a) Comparison of correlation coefficient；(b) Comparison of signal-to-noise ratio； (c) Comparison of energy percentage

由圖10可以看出，EMD-ICA方法得到的相關(guān)系數(shù)、信噪比和能量百分比的各項(xiàng)指標(biāo)均優(yōu)于小波和EMD方法的指標(biāo)。如表2所示，對(duì)其進(jìn)行配對(duì)t檢驗(yàn)，EMD-ICA相比EMD方法信噪比顯著提高(5.19±1.29vs14.87±3.04，P<0.001)，相比小波方法能量百分比顯著提高 (88.81±2.81vs96.64±2.92，P<0.001)。雖然小波降噪方法的信噪比較高，但是相關(guān)系數(shù)和能量百分比相對(duì)于其他兩種方法最低。EMD降噪方法雖然相關(guān)系數(shù)和能量百分比相對(duì)較高，但是信噪比相比其他方法最低。3種方法比較來(lái)看，小波降噪方法的降噪效果雖較好，但也損失了部分原信號(hào)的信號(hào)特征；EMD降噪方法能較好地保持原信號(hào)特征，接近原信號(hào)，但降噪效果不明顯；EMD-ICA降噪方法的降噪效果最明顯，最大限度地保留了原信號(hào)的信號(hào)特征。

3 討論

消除噪聲是BCG信號(hào)預(yù)處理的重要內(nèi)容，由于采集系統(tǒng)和外部環(huán)境等因素的影響，采集的BCG信號(hào)數(shù)據(jù)含有不同類別的噪聲，具有非線性、非平穩(wěn)性，嚴(yán)重影響心率等生理信息測(cè)量的準(zhǔn)確性。綜合考慮降噪效果和計(jì)算量，本方法避免了文獻(xiàn)[23]中需要大量實(shí)驗(yàn)確定小波基與分解尺度的問(wèn)題。從信號(hào)分解角度，小波變換中基函數(shù)是固定的，不能匹配信號(hào)的全部，用自適應(yīng)EMD分解方法代替復(fù)雜的小波變化，大大降低了算法的復(fù)雜度。文獻(xiàn)[14]中EMD分解BCG信號(hào)確定IMF階數(shù)后，直接對(duì)各分量進(jìn)行重構(gòu)，不能有效消除IMF中由于高頻噪聲的干擾以及時(shí)間尺度信號(hào)變化等因素引起的模態(tài)混疊噪聲。本研究提出信號(hào)時(shí)序平移方法構(gòu)造多個(gè)觀測(cè)信號(hào)，通過(guò)ICA可以實(shí)現(xiàn)混疊頻率的有效分離，從而估計(jì)出IMF中的原始信號(hào)，能有效消除EMD分解過(guò)程中模態(tài)混疊現(xiàn)象對(duì)降噪結(jié)果的影響，提高降噪效果。

EMD方法因在處理非線性、非平穩(wěn)信號(hào)上的優(yōu)越性而得到了廣泛的應(yīng)用，能夠自適應(yīng)地分解信號(hào)，但其分解過(guò)程中面臨的端點(diǎn)效應(yīng)問(wèn)題阻礙了其發(fā)展。對(duì)于該問(wèn)題，目前提出的很多方法都是從算法上入手來(lái)抑制端點(diǎn)效應(yīng)的影響，忽略了噪聲信號(hào)在EMD分解過(guò)程中對(duì)端點(diǎn)效應(yīng)的影響。本研究選擇增加BCG信號(hào)長(zhǎng)度，以減少端點(diǎn)效應(yīng)對(duì)信號(hào)降噪的影響。此外，分離信號(hào)的幅度和順序上存在不確定性，但是BCG信號(hào)的特征存在于相對(duì)幅頻特征中，因此對(duì)BCG信號(hào)的研究影響不大。

由于壓電薄膜傳感器的測(cè)量原理，BCG信號(hào)極易受到身體晃動(dòng)的干擾，咳嗽、扭動(dòng)身體等引起的輕微運(yùn)動(dòng)都會(huì)影響信號(hào)的質(zhì)量，使得特征難以辨認(rèn)，不能有效地提取信號(hào)。本系統(tǒng)目前適用于安靜狀態(tài)下獲取BCG信號(hào)，比較適合睡眠時(shí)檢測(cè)，特別適用于家庭應(yīng)用。因此，未來(lái)工作首先要解決BCG信號(hào)在各種條件下的采集與降噪問(wèn)題，以保證BCG信號(hào)的質(zhì)量。其次，隨著對(duì)BCG信號(hào)生理意義研究的不斷深入，希望探求一種既能代替心電圖實(shí)現(xiàn)多參數(shù)無(wú)感檢測(cè)又可以將BCG信號(hào)包含的心血管系統(tǒng)生理信息用于輔助醫(yī)生的方法，從而對(duì)各種心血管疾病進(jìn)行診斷與分析。這是今后研究的方向，其實(shí)現(xiàn)將有力于促進(jìn)人類健康水平的提高。

4 結(jié)論

本研究綜合考慮EMD與ICA的技術(shù)優(yōu)點(diǎn)，提出一種基于EMD-ICA的BCG信號(hào)降噪方法。該方法利用EMD將BCG信號(hào)分解，根據(jù)模態(tài)相關(guān)準(zhǔn)則將分解信號(hào)進(jìn)行信號(hào)層與噪聲層的判定，得到重構(gòu)信號(hào)并作為ICA的輸入，最后去除了原始信號(hào)中的噪聲成分。該方法克服了EMD分解過(guò)程中產(chǎn)生的模態(tài)混疊問(wèn)題，同時(shí)解決了ICA方法中多維輸入的要求，實(shí)現(xiàn)了多周期BCG信號(hào)的降噪。通過(guò)與小波和EMD降噪方法進(jìn)行對(duì)比實(shí)驗(yàn)研究，表明該方法能夠更有效地識(shí)別心臟的動(dòng)力學(xué)信號(hào)，這對(duì)后期研究日常作息條件下與人體健康狀態(tài)相關(guān)的生理病理信息有重要意義。