基于EKF+EKS的BCG動態高斯模型濾波研究

王子民 甘智宇 劉振丙

1(桂林電子科技大學計算機與信息安全學院 廣西 桂林 541000)2(桂林電子科技大學電子工程與自動化學院 廣西 桂林 541000)3(廣西自動檢測技術與儀器重點實驗室 廣西 桂林 541000)4(廣西信息科學實驗中心 廣西 桂林 541004)

0 引 言

BCG信號是心臟在周期性地泵血過程中產生的微弱且復雜的作用力[1]。基于牛頓第三定律，在心臟規律的收縮和舒張過程中，血液流經動靜脈又回到心臟的過程會對人體產生反作用力。當血液因心房和心室的收縮流出時對人體產生作用力;當經過血管不同位置，尤其是流經大血管的彎曲處時由于沖擊力作用對身體產生反作用力；血液流動時血管的彈性運動產生的作用力。以上三種作用力的合力引起了心沖擊運動[1-2]，心沖擊圖信號源于心臟機械收縮以及大血管中的血液被突然加速而對身體產生的反作用力，因此BCG信號可以反映心臟的機械活動[3-4]，圖1所示為BCG信號波形。

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

基于微弱信號的動態模型必然關注兩個方面：濾波算法和系統的動態模型。BCG的動態高斯模型合成的信號屬于非線性周期信號，因此必須采用自適應的非線性的濾波方法。通過比較經典濾波方法、小波去噪和擴展卡爾曼濾波(Extended Kalman filter，EKF)對BCG信號進行濾波，因為BCG信號太微弱，單一的濾波方法不能滿足濾波要求。由于Reza Sameni[8]已將擴展卡爾曼濾波和擴展卡爾曼平滑運用至ECG信號的降噪，基于EKF和擴展卡爾曼平滑(Extended Kalman smoother，EKS)的特性。本文提出一種聯合擴展卡爾曼濾波和擴展卡爾曼平滑濾波方法對BCG信號進行復合濾波。由于EKF濾波是需要微弱信號真實值的，但信號的真實值通常不能確定，需要動態模型的值來作為BCG信號真實值的估計值，因此可將動態高斯模型模合成信號的數值作為濾波過程中真實值的估計值。

1 相關工作

1.1 ECG信號動態高斯模型

2003年McSharry提出一種基于高斯函數的ECG信號動態高斯模型[6]，該模型可以合成健康個體的ECG信號，ECG信號的任何一點都可以用笛卡爾坐標(x,y,z)表示，函數方程組如下所示：

(1)

表1 ECG模型參數

1.2 PCG信號動態高斯模型

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

(2)

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

表2 PCG模型的參數

通過選擇不同的模型參數來合成不同形態和不同周期的PCG波形，為評估處理醫學信號PCG的方法提供基準。本文關于BCG信號動態高斯模型的靈感來源于ECG信號和PCG信號的動態高斯模型，由于BCG信號、ECG信號和PCG信號有許多相似的地方，都屬于生物醫學微弱信號，具有信號頻率低、幅值低和噪聲強等特點。因此可將ECG信號和PCG信號動態高斯模型的思想運用到BCG信號。

1.3 擴展卡爾曼濾波和擴展卡爾曼平滑

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

(3)

(4)

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

(5)

擴展卡爾曼濾波包括兩個階段，預測和更新。預測是根據過去的狀態判斷當前狀態的值，更新階段指的是當前狀態的測試信息，用于重新精確修改之前的狀態值，之前的狀態值結果會反饋當前狀態的值。

預測方程為：

(6)

創新方程為：

(7)

更新方程為：

(8)

擴展卡爾曼平滑起源于卡爾曼平滑，EKS也是卡爾曼平滑通過泰勒級數一階展開忽略高次項的擴展方法，EKS相對擴展卡爾曼濾波前向濾波是反向平滑濾波，因篇幅有限，此處不再敘述，有興趣的讀者可以閱讀文獻[10-12]。

2 方 法

2.1 BCG信號的采集

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

圖2 BCG信號采集和模擬電路處理

2.2 BCG信號動態高斯模型

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

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

2.3 聯合擴展卡爾曼濾波與擴展卡爾曼平滑的復合濾波

本論文對BCG信號采用的降噪濾波方法有4種，經典濾波方法指傳統的數字濾波高通、低通和帶通濾波器有較好的濾波效果[8,17-18]。經典的BCG信號數字濾波方式是采用高通濾波濾除呼吸信號，低通濾波濾掉因為環境因素、肢體活動帶來的干擾，最后是工頻濾波去除50 Hz的工頻干擾信號[20]。小波濾波采用的是DB5小波進行6尺度對BCG信號進行分解。卡爾曼濾波通常只適用于線性系統[17-18]，用于非線性的信號濾波表現不佳。

由于BCG信號屬于非線性的準周期信號。若想將卡爾曼濾波用于BCG信號處理，則必須將卡爾曼濾波的線性關系轉化為非線性的關系，因此可采用擴展卡爾曼濾波方法。最常用的方法是卡爾曼濾波通過泰勒級數一階展開忽略高次項的擴展方法。在此基礎上形成擴展卡爾曼濾波[16-17]。卡爾曼平滑濾波也通常適用于線性系統的濾波。EKS也是卡爾曼平滑通過泰勒級數一階展開忽略高次項的擴展方法。EKF主要有兩大缺陷，其一是未考慮信號的分布誤差，其二是認為獨立的線性系統產生狀態誤差[21]。當微弱信號經過EKF濾波后，可再經過EKS濾波，而EKS反向平滑濾波可縮小信號的分布誤差和非線性系統產生的狀態誤差。

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

圖3 EKF+EKS在BCG信號的處理

3 實驗結果與討論

3.1 BCG信號的經典濾波處理

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

圖4 BCG的原始信號

圖5 經典濾波器處理后輸出的BCG信號

3.2 動態高斯模型合成的BCG信號

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

圖6 單周期BCG信號三維展示

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

類比推理，根據不同的高斯模型參數，可以相應地生成不同形態的BCG信號。為方便研究多周期的BCG信號，圖6的三維空間展示可轉換到二維空間展示，并加以周期延拓，形成多周期的合成BCG信號。

圖7為單周期BCG信號，該模型合成BCG信號保留了H波、I波、J波、K波、L波和M波等特征，其周期延拓結果如圖8所示。基于BCG信號的動態高斯模型，可根據需求，設置不同的模型參數，進而生成不同周期和不同形態的健康個體BCG信號。

圖7 動態高斯模型產生的單周期BCG信號

圖8 模型產生的多周期BCG信號

3.3 模型生成的BCG信號與健康個體的信號比較

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

圖9 模型合成的BCG信號添加12 dB的高斯白噪聲

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

圖10 模型合成的BCG信號特征點提取結果

健康個體采集的BCG信號如圖11所示，與圖9和圖10比較后發現，模型合成的BCG信號保留了健康個體BCG信號的有效周期特征，動態高斯模型合成的BCG信號可替代健康個體產生的BCG信號[24]。

圖11 健康個體采集的BCG信號

3.4 聯合擴展卡爾曼濾波與擴展卡爾曼平滑復合濾波

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

圖12 向動態高斯模型合成的BCG信號中人工加入SNR=10 dB的高斯白噪聲

圖13 經典濾波處理圖12信號后的BCG信號

圖14 經過小波DB5降噪處理圖12信號后的BCG信號

圖15 經過EKF處理圖12后的BCG信號

圖16 經過EKF+EKS處理圖12后的BCG信號

為了更好地比較不同方法的降噪的效果，我們對動態高斯模模型合成的BCG信號添加SNR數值范圍為6～18 dB的等信噪比間隔的高斯白噪聲，然后比較其濾波效果。

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

表4 不同濾波方法輸出信號的信噪比 dB

4 結 語

通過理解動態ECG信號和PCG信號的動態高斯模型，本文實現了一種基于高斯核函數的BCG信號模型。通過設置不同的模型參數，可以生成不同形態和不同周期的BCG信號，包含H波、I波、J波、k波、L波和M波的特征。與健康個體產生的BCG信號比較后發現，動態高斯模型可替代健康個體合成的BCG信號。在動態高斯模型合成BCG信號的基礎之上，運用了一種聯合擴展卡爾曼濾波與擴展卡爾曼平滑的方法對BCG信號進行復合濾波，通過與經典濾波、擴展卡爾曼濾波和小波濾波后的信號信噪比值比較后發現，聯合擴展卡爾曼濾波與擴展卡爾曼平滑復合濾波方法對BCG信號有更好的降噪效果。