王子民 甘智宇 劉振丙

1(桂林電子科技大學(xué)計(jì)算機(jī)與信息安全學(xué)院 廣西 桂林 541000)2(桂林電子科技大學(xué)電子工程與自動(dòng)化學(xué)院 廣西 桂林 541000)3(廣西自動(dòng)檢測(cè)技術(shù)與儀器重點(diǎn)實(shí)驗(yàn)室 廣西 桂林 541000)4(廣西信息科學(xué)實(shí)驗(yàn)中心 廣西 桂林 541004)

0 引 言

BCG信號(hào)是心臟在周期性地泵血過程中產(chǎn)生的微弱且復(fù)雜的作用力[1]?；谂ｎD第三定律，在心臟規(guī)律的收縮和舒張過程中，血液流經(jīng)動(dòng)靜脈又回到心臟的過程會(huì)對(duì)人體產(chǎn)生反作用力。當(dāng)血液因心房和心室的收縮流出時(shí)對(duì)人體產(chǎn)生作用力;當(dāng)經(jīng)過血管不同位置，尤其是流經(jīng)大血管的彎曲處時(shí)由于沖擊力作用對(duì)身體產(chǎn)生反作用力；血液流動(dòng)時(shí)血管的彈性運(yùn)動(dòng)產(chǎn)生的作用力。以上三種作用力的合力引起了心沖擊運(yùn)動(dòng)[1-2]，心沖擊圖信號(hào)源于心臟機(jī)械收縮以及大血管中的血液被突然加速而對(duì)身體產(chǎn)生的反作用力，因此BCG信號(hào)可以反映心臟的機(jī)械活動(dòng)[3-4]，圖1所示為BCG信號(hào)波形。

心臟病已經(jīng)成為了國民健康的最大威脅，BCG信號(hào)可為心臟疾病提供有效的檢測(cè)、診斷和治療信息[5]。BCG信號(hào)單周期內(nèi)有相同的6個(gè)特征點(diǎn)，但每個(gè)周期形態(tài)有差異，同一人在不同健康狀態(tài)或者不同環(huán)境條件采集的BCG信號(hào)有差異。隨著研究的深入，科研和醫(yī)務(wù)人員均需要克服采集環(huán)境和個(gè)體差異影響的正常健康個(gè)體BCG信號(hào)，McSharry[6]提出的心電圖(Electrocardiogram，ECG)信號(hào)動(dòng)態(tài)高斯模型和Almasi[7]提出的心音(Phonocardiogram，PCG)信號(hào)的動(dòng)態(tài)高斯模型，由于BCG信號(hào)、ECG信號(hào)和PCG信號(hào)均屬于生物醫(yī)學(xué)的微弱信號(hào)，因此可以設(shè)計(jì)一個(gè)合成BCG信號(hào)的動(dòng)態(tài)高斯模型。

基于微弱信號(hào)的動(dòng)態(tài)模型必然關(guān)注兩個(gè)方面：濾波算法和系統(tǒng)的動(dòng)態(tài)模型。BCG的動(dòng)態(tài)高斯模型合成的信號(hào)屬于非線性周期信號(hào)，因此必須采用自適應(yīng)的非線性的濾波方法。通過比較經(jīng)典濾波方法、小波去噪和擴(kuò)展卡爾曼濾波(Extended Kalman filter，EKF)對(duì)BCG信號(hào)進(jìn)行濾波，因?yàn)锽CG信號(hào)太微弱，單一的濾波方法不能滿足濾波要求。由于Reza Sameni[8]已將擴(kuò)展卡爾曼濾波和擴(kuò)展卡爾曼平滑運(yùn)用至ECG信號(hào)的降噪，基于EKF和擴(kuò)展卡爾曼平滑(Extended Kalman smoother，EKS)的特性。本文提出一種聯(lián)合擴(kuò)展卡爾曼濾波和擴(kuò)展卡爾曼平滑濾波方法對(duì)BCG信號(hào)進(jìn)行復(fù)合濾波。由于EKF濾波是需要微弱信號(hào)真實(shí)值的，但信號(hào)的真實(shí)值通常不能確定，需要?jiǎng)討B(tài)模型的值來作為BCG信號(hào)真實(shí)值的估計(jì)值，因此可將動(dòng)態(tài)高斯模型模合成信號(hào)的數(shù)值作為濾波過程中真實(shí)值的估計(jì)值。

1 相關(guān)工作

1.1 ECG信號(hào)動(dòng)態(tài)高斯模型

2003年McSharry提出一種基于高斯函數(shù)的ECG信號(hào)動(dòng)態(tài)高斯模型[6]，該模型可以合成健康個(gè)體的ECG信號(hào)，ECG信號(hào)的任何一點(diǎn)都可以用笛卡爾坐標(biāo)(x,y,z)表示，函數(shù)方程組如下所示：

(1)

表1 ECG模型參數(shù)

1.2 PCG信號(hào)動(dòng)態(tài)高斯模型

2011年Almasi提出一種PCG信號(hào)的動(dòng)態(tài)模型[7]，能合成正常健康個(gè)體的PCG信號(hào)。PCG信號(hào)的動(dòng)態(tài)模型是基于高斯核函數(shù)的動(dòng)態(tài)模型，與ECG信號(hào)不同的是，PCG信號(hào)模型的高斯核函數(shù)表現(xiàn)形式采用極坐標(biāo)的形式。心動(dòng)周期包括心臟收縮期和舒張期，心音信號(hào)可以分為4個(gè)階段，S1、S2、S3和S4。S3和S4必須通過特殊手段才能測(cè)聽到，相對(duì)而言S1和S2比較容易測(cè)聽到，因此PCG信號(hào)的動(dòng)態(tài)模型，通常只研究心音信號(hào)的前兩個(gè)階段S1和S2。PCG信號(hào)的高斯函數(shù)動(dòng)態(tài)模型極坐標(biāo)形式方程組如下：

(2)

式中：θ=atan2(y,x)代表信號(hào)的相位值，取值范圍-π≤θ≤π，w代表角加速度，αi、μi和σi對(duì)應(yīng)于方程組高斯核函數(shù)的幅度、寬度和相位等參數(shù)，fi代表頻率，φi代表相位的遷移。與ECG信號(hào)動(dòng)態(tài)高斯模型相同，只要模型參數(shù)確定后，PCG信號(hào)的動(dòng)態(tài)高斯模型合成的PCG信號(hào)的波形和周期也就確定。PCG信號(hào)模型參數(shù)如表2所示。

表2 PCG模型的參數(shù)

通過選擇不同的模型參數(shù)來合成不同形態(tài)和不同周期的PCG波形，為評(píng)估處理醫(yī)學(xué)信號(hào)PCG的方法提供基準(zhǔn)。本文關(guān)于BCG信號(hào)動(dòng)態(tài)高斯模型的靈感來源于ECG信號(hào)和PCG信號(hào)的動(dòng)態(tài)高斯模型，由于BCG信號(hào)、ECG信號(hào)和PCG信號(hào)有許多相似的地方，都屬于生物醫(yī)學(xué)微弱信號(hào)，具有信號(hào)頻率低、幅值低和噪聲強(qiáng)等特點(diǎn)。因此可將ECG信號(hào)和PCG信號(hào)動(dòng)態(tài)高斯模型的思想運(yùn)用到BCG信號(hào)。

1.3 擴(kuò)展卡爾曼濾波和擴(kuò)展卡爾曼平滑

卡爾曼濾波不能直接用于非線性濾波，本小節(jié)主要闡述EKF的濾波框架[9]。EKF通常適用于離散的非線性系統(tǒng)，一個(gè)離散的系統(tǒng)可以表述如下:

(3)

(4)

其中:Ak、Fk、Ck、Gk定義如下：

(5)

擴(kuò)展卡爾曼濾波包括兩個(gè)階段，預(yù)測(cè)和更新。預(yù)測(cè)是根據(jù)過去的狀態(tài)判斷當(dāng)前狀態(tài)的值，更新階段指的是當(dāng)前狀態(tài)的測(cè)試信息，用于重新精確修改之前的狀態(tài)值，之前的狀態(tài)值結(jié)果會(huì)反饋當(dāng)前狀態(tài)的值。

預(yù)測(cè)方程為：

(6)

創(chuàng)新方程為：

(7)

更新方程為：

(8)

擴(kuò)展卡爾曼平滑起源于卡爾曼平滑，EKS也是卡爾曼平滑通過泰勒級(jí)數(shù)一階展開忽略高次項(xiàng)的擴(kuò)展方法，EKS相對(duì)擴(kuò)展卡爾曼濾波前向?yàn)V波是反向平滑濾波，因篇幅有限，此處不再敘述，有興趣的讀者可以閱讀文獻(xiàn)[10-12]。

2 方 法

2.1 BCG信號(hào)的采集

BCG信號(hào)的采集使用壓電薄膜(PVDF)傳感器[9]，安裝在凳子上，屬于無創(chuàng)傷測(cè)試方案，對(duì)用戶測(cè)試造成的心理壓力比較小。PVDF屬于模擬電壓傳感器，屬于線性傳感器，壓電薄膜所受到的壓力與輸出電壓呈正相關(guān)[14-15]。電荷放大電路主要是將電壓信號(hào)放大，通過低通濾波后放大信號(hào)。電壓跟隨器的主要作用是隔離前后一級(jí)的低通濾波放大電路[13-15]。BCG信號(hào)的前端模擬處理電路的硬件框圖如圖2所示。

圖2 BCG信號(hào)采集和模擬電路處理

2.2 BCG信號(hào)動(dòng)態(tài)高斯模型

BCG波形包含F(xiàn)波、G波、H波、I波、J波、K波、L波和M波，每一個(gè)波代表相鄰區(qū)域的一個(gè)極大值或者極小值。F波和G波在IJK波群之前，G波位置緊隨F波之后，均屬于不明顯的波形，F(xiàn)波和G波反映了心臟收縮早期活動(dòng)。心臟收縮開始期是H波，I波和J波屬于心臟收縮期的主要特征，J波位于I波之后，且J波具有BCG波形的最大幅值，K波是心臟收縮期結(jié)束的標(biāo)志。L波位置在K波之后，M波位置在L波之后，N波屬于BCG波形的最后階段。L波和M波屬于心臟舒張?jiān)缙?，N波屬于心臟舒張中期和末期。由于N波波形特征不明顯，F(xiàn)波和G波也如此，通常忽視掉這三個(gè)波形。BCG波形特征也就主要體現(xiàn)在H波、I波、J波、k波、L波和M波[4]。基于上述原理，本文的BCG信號(hào)動(dòng)態(tài)高斯模型主要包含上述6個(gè)波的完整特征。ECG信號(hào)與BCG信號(hào)的波形相似度比較高[16]，BCG波形和ECG波形均是頻率低、幅值低和噪聲強(qiáng)等特點(diǎn)，且都是通過信號(hào)的波峰和波谷表示信號(hào)特征，因此可用下式表示[6-7]：

式中：θ表示信號(hào)的相位，b表示信號(hào)的寬度，α表示信號(hào)的幅值，高斯核函數(shù)的就是將BCG信號(hào)的特征點(diǎn)H波、I波、J波、K波、L波、M波的每一個(gè)特征點(diǎn)用幅值、寬度和相位來表示。z0表示為呼吸作用，z0(t)=Asin(2πf2t)表示隨著呼吸頻率而產(chǎn)生的BCG信號(hào)的基線A=0.15 mV[6]。

2.3 聯(lián)合擴(kuò)展卡爾曼濾波與擴(kuò)展卡爾曼平滑的復(fù)合濾波

本論文對(duì)BCG信號(hào)采用的降噪濾波方法有4種，經(jīng)典濾波方法指?jìng)鹘y(tǒng)的數(shù)字濾波高通、低通和帶通濾波器有較好的濾波效果[8,17-18]。經(jīng)典的BCG信號(hào)數(shù)字濾波方式是采用高通濾波濾除呼吸信號(hào)，低通濾波濾掉因?yàn)榄h(huán)境因素、肢體活動(dòng)帶來的干擾，最后是工頻濾波去除50 Hz的工頻干擾信號(hào)[20]。小波濾波采用的是DB5小波進(jìn)行6尺度對(duì)BCG信號(hào)進(jìn)行分解?？柭鼮V波通常只適用于線性系統(tǒng)[17-18]，用于非線性的信號(hào)濾波表現(xiàn)不佳。

由于BCG信號(hào)屬于非線性的準(zhǔn)周期信號(hào)。若想將卡爾曼濾波用于BCG信號(hào)處理，則必須將卡爾曼濾波的線性關(guān)系轉(zhuǎn)化為非線性的關(guān)系，因此可采用擴(kuò)展卡爾曼濾波方法。最常用的方法是卡爾曼濾波通過泰勒級(jí)數(shù)一階展開忽略高次項(xiàng)的擴(kuò)展方法。在此基礎(chǔ)上形成擴(kuò)展卡爾曼濾波[16-17]?？柭交瑸V波也通常適用于線性系統(tǒng)的濾波。EKS也是卡爾曼平滑通過泰勒級(jí)數(shù)一階展開忽略高次項(xiàng)的擴(kuò)展方法。EKF主要有兩大缺陷，其一是未考慮信號(hào)的分布誤差，其二是認(rèn)為獨(dú)立的線性系統(tǒng)產(chǎn)生狀態(tài)誤差[21]。當(dāng)微弱信號(hào)經(jīng)過EKF濾波后，可再經(jīng)過EKS濾波，而EKS反向平滑濾波可縮小信號(hào)的分布誤差和非線性系統(tǒng)產(chǎn)生的狀態(tài)誤差。

本文采用一種基于聯(lián)合擴(kuò)展卡爾曼濾波與擴(kuò)展卡爾曼平滑的方法對(duì)BCG信號(hào)進(jìn)行復(fù)合濾波。擴(kuò)展卡爾曼濾波是正向迭代濾波，擴(kuò)展卡爾曼平滑是反向平滑濾波[7,20]。BCG信號(hào)先經(jīng)過擴(kuò)展卡爾曼正向?yàn)V波，濾除掉主要的噪聲干擾，而后BCG信號(hào)再次通過擴(kuò)展卡爾曼平滑濾波，對(duì)BCG信號(hào)反向平滑濾波，模糊其他因素的干擾。EKF+EKS復(fù)合濾波核心思想是將兩種濾波方法組合起來，BCG信號(hào)先進(jìn)行擴(kuò)展卡爾曼濾波后，再進(jìn)行擴(kuò)展卡爾曼平滑濾波，濾波效果顯著提高。聯(lián)合擴(kuò)展卡爾曼濾波與擴(kuò)展卡爾曼平滑的復(fù)合濾波原理[22]如圖3所示。

圖3 EKF+EKS在BCG信號(hào)的處理

3 實(shí)驗(yàn)結(jié)果與討論

3.1 BCG信號(hào)的經(jīng)典濾波處理

壓電薄膜傳感器采集BCG信號(hào)的原始數(shù)據(jù)后，采用MATLAB數(shù)字濾波器進(jìn)行信號(hào)經(jīng)典濾波處理。經(jīng)典濾波處理采用切比雪夫II型濾波器，切比雪夫II型3階低通濾波器，阻帶紋波30 dB，截止頻率0.05 Hz，切比雪夫II型3階高通濾波器，阻帶紋紋波30 dB，截止頻0.002 Hz，最后是采用數(shù)字陷波器去除50 Hz的工頻干擾。圖4所示為BCG采集的原始信號(hào)，圖5為經(jīng)典濾波處理后的BCG信號(hào)。

圖4 BCG的原始信號(hào)

圖5 經(jīng)典濾波器處理后輸出的BCG信號(hào)

3.2 動(dòng)態(tài)高斯模型合成的BCG信號(hào)

本文提出的BCG動(dòng)態(tài)高斯模型在笛卡爾坐標(biāo)系生成一個(gè)單周期BCG信號(hào)的三維圖像，如圖6所示。在一個(gè)三維空間中，圍繞著單位圓的軌跡顯示BCG信號(hào)的準(zhǔn)周期性質(zhì)。軌跡的每一次循環(huán)表示一個(gè)心動(dòng)周期的BCG信號(hào)。在x、y和z方向上的數(shù)值變化代表BCG波形的幅值變化。根據(jù)式(9)所描述的BCG動(dòng)態(tài)高斯模型，利用BCG信號(hào)進(jìn)曲線擬合與參數(shù)估計(jì)[23]，可以得到BCG波形中極大值和極小值點(diǎn)的動(dòng)態(tài)高斯模型參數(shù)αi、bi和θi。表3所示為BCG信號(hào)的一組動(dòng)態(tài)高斯模型參數(shù)值。

圖6 單周期BCG信號(hào)三維展示

Index(i)HIJKLMTime/s-0.1220-0.115000.06210.10230.1280θi/rad-π2-3π100π42π5π2αi0.5193-3.8544.741-2.178-1.160-1.611bi0.24310.2089-0.14010.37290.56150.4935

類比推理，根據(jù)不同的高斯模型參數(shù)，可以相應(yīng)地生成不同形態(tài)的BCG信號(hào)。為方便研究多周期的BCG信號(hào)，圖6的三維空間展示可轉(zhuǎn)換到二維空間展示，并加以周期延拓，形成多周期的合成BCG信號(hào)。

圖7為單周期BCG信號(hào)，該模型合成BCG信號(hào)保留了H波、I波、J波、K波、L波和M波等特征，其周期延拓結(jié)果如圖8所示?；贐CG信號(hào)的動(dòng)態(tài)高斯模型，可根據(jù)需求，設(shè)置不同的模型參數(shù)，進(jìn)而生成不同周期和不同形態(tài)的健康個(gè)體BCG信號(hào)。

圖7 動(dòng)態(tài)高斯模型產(chǎn)生的單周期BCG信號(hào)

圖8 模型產(chǎn)生的多周期BCG信號(hào)

3.3 模型生成的BCG信號(hào)與健康個(gè)體的信號(hào)比較

動(dòng)態(tài)高斯模型合成的BCG信號(hào)屬于無干擾、無噪聲和幅值穩(wěn)定的BCG信號(hào)。但通過包含壓電薄膜的模擬電路采集的信號(hào)不能完全消除環(huán)境、設(shè)備和人為因素的干擾，因此動(dòng)態(tài)高斯模型合成的BCG信號(hào)需人工添加一定比例的高斯白噪聲，目的是為模擬健康個(gè)體產(chǎn)生的BCG信號(hào)動(dòng)態(tài)高斯模型合成的BCG信號(hào)添加SNR=12 dB的噪聲后的信號(hào)如圖9所示。

圖9 模型合成的BCG信號(hào)添加12 dB的高斯白噪聲

雖然BCG信號(hào)屬于非線性信號(hào)，但是每個(gè)周期的波形有規(guī)律性，BCG信號(hào)一個(gè)周期包含H波、I波、J波、K波、L波和M波等6個(gè)特征點(diǎn)，且特征點(diǎn)出現(xiàn)的先后順序和幅值大小順序保持不變。通過判斷特征點(diǎn)的出現(xiàn)先后順序、特征點(diǎn)的幅值大小順序和信號(hào)走向趨勢(shì)判斷是否屬于健康個(gè)體的BCG信號(hào)。由于模型的參數(shù)可以調(diào)節(jié)，可以產(chǎn)生不同形態(tài)不同周期的BCG信號(hào)。但無論何種形態(tài)的健康個(gè)體BCG信號(hào)的J波幅值最大和H波、I波、J波、K波、L波和M波出現(xiàn)先后不變。因此通過特征點(diǎn)的數(shù)量及其順序和特征點(diǎn)幅值可判斷動(dòng)態(tài)高斯模型合成的BCG信號(hào)可替代健康個(gè)體產(chǎn)生的信號(hào)。圖10為提取了該BCG信號(hào)的所有的特征點(diǎn)，并標(biāo)識(shí)了所有周期BCG波形包含的6個(gè)特征，分別是H波、I波、J波、K波、L波、M波，且先后順序和幅值大小一致。

圖10 模型合成的BCG信號(hào)特征點(diǎn)提取結(jié)果

健康個(gè)體采集的BCG信號(hào)如圖11所示，與圖9和圖10比較后發(fā)現(xiàn)，模型合成的BCG信號(hào)保留了健康個(gè)體BCG信號(hào)的有效周期特征，動(dòng)態(tài)高斯模型合成的BCG信號(hào)可替代健康個(gè)體產(chǎn)生的BCG信號(hào)[24]。

圖11 健康個(gè)體采集的BCG信號(hào)

3.4 聯(lián)合擴(kuò)展卡爾曼濾波與擴(kuò)展卡爾曼平滑復(fù)合濾波

本節(jié)的內(nèi)容是對(duì)比不同濾波方法的降噪效果，信號(hào)的信噪比數(shù)值越大，說明信號(hào)的噪聲越低，該方法的降噪效果越好。圖12所示為向動(dòng)態(tài)高斯模型合成的BCG信號(hào)中人工加入SNR=10 dB的高斯白噪聲，圖13所示為圖12的信號(hào)經(jīng)過經(jīng)典濾波后的BCG信號(hào)，經(jīng)典濾波由傳統(tǒng)數(shù)字高通、低通和數(shù)字陷波器濾波組成。圖13中BCG信號(hào)的信噪比為17 dB。圖14所示為圖12信號(hào)經(jīng)過DB5小波進(jìn)行6層分解處理后的結(jié)果，信噪比為16 dB。圖15所示為圖12信號(hào)經(jīng)過EKF處理后的結(jié)果，信噪比為21 dB。圖16所示為圖12信號(hào)經(jīng)過EKF+EKS處理后的結(jié)果，信噪比為22 dB。

圖12 向動(dòng)態(tài)高斯模型合成的BCG信號(hào)中人工加入SNR=10 dB的高斯白噪聲

圖13 經(jīng)典濾波處理圖12信號(hào)后的BCG信號(hào)

圖14 經(jīng)過小波DB5降噪處理圖12信號(hào)后的BCG信號(hào)

圖15 經(jīng)過EKF處理圖12后的BCG信號(hào)

圖16 經(jīng)過EKF+EKS處理圖12后的BCG信號(hào)

為了更好地比較不同方法的降噪的效果，我們對(duì)動(dòng)態(tài)高斯模模型合成的BCG信號(hào)添加SNR數(shù)值范圍為6～18 dB的等信噪比間隔的高斯白噪聲，然后比較其濾波效果。

表4表明，具有相同信噪比的BCG信號(hào)通過代表4種不同方法的濾波器,分別是經(jīng)典濾波器、DB5小波進(jìn)行6尺度分解、EKF濾波和EKF+EKS復(fù)合濾波。通過濾波器后BCG信號(hào)的信噪比相比較,通過EKF+EKS的BCG信號(hào)的信噪比數(shù)值明顯大于通過其他三種濾波器BCG信號(hào)的信噪比,證明針對(duì)BCG信號(hào)的濾波,基于EKS+EKF復(fù)合濾波效果比傳統(tǒng)數(shù)字濾波降噪效果更好。

表4 不同濾波方法輸出信號(hào)的信噪比 dB

4 結(jié) 語

通過理解動(dòng)態(tài)ECG信號(hào)和PCG信號(hào)的動(dòng)態(tài)高斯模型，本文實(shí)現(xiàn)了一種基于高斯核函數(shù)的BCG信號(hào)模型。通過設(shè)置不同的模型參數(shù)，可以生成不同形態(tài)和不同周期的BCG信號(hào)，包含H波、I波、J波、k波、L波和M波的特征。與健康個(gè)體產(chǎn)生的BCG信號(hào)比較后發(fā)現(xiàn)，動(dòng)態(tài)高斯模型可替代健康個(gè)體合成的BCG信號(hào)。在動(dòng)態(tài)高斯模型合成BCG信號(hào)的基礎(chǔ)之上，運(yùn)用了一種聯(lián)合擴(kuò)展卡爾曼濾波與擴(kuò)展卡爾曼平滑的方法對(duì)BCG信號(hào)進(jìn)行復(fù)合濾波，通過與經(jīng)典濾波、擴(kuò)展卡爾曼濾波和小波濾波后的信號(hào)信噪比值比較后發(fā)現(xiàn)，聯(lián)合擴(kuò)展卡爾曼濾波與擴(kuò)展卡爾曼平滑復(fù)合濾波方法對(duì)BCG信號(hào)有更好的降噪效果。