基于經驗模態分解和形態濾波的血壓測量研究

甘永進，莫 沛，陳時東，鄭金存，*，蔣曲博

(1. 玉林師范學院 電子與通信工程學院，廣西 玉林 537000；2. 廣西科技大學鹿山學院 電氣與計算機工程系，廣西 柳州 545000；3. 桂林電子科技大學 電子工程與自動化學院，廣西 桂林 541004)

0 引言

目前，示波型電子血壓計是市場的主流[1]。然而因心率異常導致振蕩波不規則以及血管形變恢復需要一定時間，示波型電子血壓計無法實現連續測量每搏血壓的變化[2]，而連續血壓測量在家庭保健和臨床診斷中意義重大。通過光電容積脈搏波(Photoplethysmography,PPG)完成對每個心動周期血壓值的檢測，為實現血壓的無創連續測量提供一定的參考價值。但是PPG信號容易受到采集系統帶來的白噪聲、工頻干擾[3]以及由肌電干擾等因素引起的基線漂移等噪聲的干擾，對后續血壓測量的精度帶來影響。

為實現由PPG信號達到連續精準地測量血壓的目的，本文通過TI系列AFE44X0芯片結合反射式傳感器DCM03設計光電容積脈搏波采集系統，設計經驗模態分解(Empirical Mode Decomposition,EMD)結合形態學濾波的方法對原始PPG信號進行自適應濾波以及基線去除工作，對濾波后的PPG信號進行二次微分得到加速脈搏波，采用Adaboost算法對加速脈搏波進行分類，對分類后的加速脈搏波進行PPG信號傳導時間(Pulse Wave Transit Time， PWTT)的提取，最后建立不同的血壓測量模型。

1 光電容積脈搏波去噪算法

1.1 EMD去噪

經驗模態分解由數據時間尺度特征來進行信號分解以得到一系列的本征模函數(Intrinsic Mode Function，IMF)，其表征信號的不同振蕩模式，包含不同時間尺度[4-5]信號局部特征。其中，各個IMF中零點和極點個數相等或相差1，且局部極大值和極小值確定的信號包絡線的均值為0[6]，使得定義瞬時頻率的同時可避免由信號不對稱引起的瞬時頻率波動[7]。

本文采用光電反射式傳感器DCM03結合血氧模擬前端AFE44X0系列芯片進行反射式光電容積脈搏波的采集。反射式傳感器DCM03的光發射器發出光線照射到指尖后發生漫反射，DCM03的光接收器接收指尖血液的光感應信號，實現光電信號轉換，再由AFE44X0進行I-V轉換、放大濾波及A/D轉換等處理得到PPG數字信號并輸出至MCU[8]。將PPG數字信號進行EMD分解， 噪聲分布在各個IMF分量上，故將所有的IMF分量進行軟閾值去噪，可有效地對原始PPG信號的噪聲進行抑制，其步驟如下:

①對原始PPG信號進行EMD分解，得到各個IMF分量，如圖1所示；

②對各IMF分量采用Stein無偏似然估計準則進行閾值估計；

③使用wthresh進行去噪，得到去噪后各IMF分量如圖2所示；

④將去噪后的各IMF分量進行重構得到去噪后的PPG信號。

圖1 容積脈搏波的EMD分解圖
Fig.1 EMD decomposition diagram of PPG

圖2 去噪后的IMF分量
Fig.2 IMF component after denoising

1.2 基于形態學的基線漂移濾除

設采樣得到的原始一維脈搏信號為f(n)，定義域為F={0,1,…,N-1}，結構元素為g(n)，定義域為G={0,1,…,M-1}，且M

(1)

(2)

據此，由腐蝕與膨脹運算組成開、閉運算為

(f°g)(n)=(f?g⊕g)(n),

(3)

(f·g)(n)=(f⊕g?g)(n)。

(4)

一般，閉運算用于填充細小空洞，實現平滑或抑制信號波谷噪聲；開運算可用于斷開窄小的連接，消除微小的尖刺，濾除信號峰值噪聲，平滑信號邊界輪廓。常采用形態開、閉的級聯形式對信號進行處理。傳統形態的開-閉和閉-開運算以不同順序級聯開閉運算，定義為

FOC(f,g)=f°g·g,

(5)

FCO(f,g)=f·g°g。

(6)

因形態閉運算的反擴展性和形態開運算的擴展性，式(5)和式(6)定義的兩種傳統的濾波器都存在統計偏移現象，即對于開-閉濾波器而言，最后的輸出幅度偏??；但對于閉-開濾波器而言，最后的輸出幅度偏大，在一般情況下，單獨使用得到的濾波效果都不是最佳的。

欲有效地抑制采集到的脈搏信號中的不同噪聲，減小最終輸出的單向偏移，由兩種濾波器的平均組合形式[9]，將形態開閉-閉開濾波器定義為

FOC-CO(f,g)=(f°g·g+f·g°g)/2。

(7)

欲將統計偏移現象更深一層地減小，根據不同尺寸的結構元素，定義廣義開-閉和閉-開濾波器為

FOC(f,g1,g2)=f°g1·g2,

(8)

FCO(f,g1,g2)=f·g1°g2。

(9)

相比于傳統的形態濾波器，廣義形態濾波器能夠對信號中的各種噪聲進行有效的抑制[10]。為校正基線漂移[11]，本文通過尺寸不同的結構元素的開-閉和閉-開運算組合進行處理。首先，設原始脈搏血氧信號為x(n)，對受基線漂移干擾的x(n)進行廣義形態閉-開運算處理；然后，將x(n)進行廣義形態開-閉運算處理；之后，再把以上兩個步驟的結果進行求和平均，得到基線分量；最后，將x(n)與基線作差，得到校正基線后的信號y(n)，形態學去除基線漂移[12]框圖如圖3所示。 由表1小波及本文算法對信號去噪后的信噪比對比知，該算法去除基線漂移效果較佳。經過EMD分解以及形態學濾波后的波形和頻譜如圖4、圖5所示，由圖4和圖5可見，信號中的高頻噪聲以及基線漂移得到一定程度的抑制，說明了設計的算法對光電容積脈搏波消噪效果的可行性。

圖3 去除基線漂移框圖
Fig.3 Block diagram of baseline wander removal

表1 幾種不同去基線算法的信噪比
Tab.1 SNR of different algorithms

序號方法SNR 1db2小波38.607 22db4小波38.216 33db5小波38.987 34db6小波38.527 55sym4小波38.587 76sym8小波38.497 37形態學濾波41.027 6

圖4 各個階段波形的頻譜圖
Fig.4 Spectrogram of each stage

圖5 各個階段波形
Fig.5 Wave form of each stage

2 血壓計算方法

2.1 脈搏波傳導時間計算

PPG信號經EMD和形態學濾波后，再提取PPG信號的傳導時間來建立PPG信號傳導時間與血壓的數學關系，以此建立血壓計算模型。其中，PPG信號及相應的加速脈搏波波形如6所示。

圖6 PPG信號及其加速脈搏波波形
Fig.6 PPG and accelerating pulse wave

圖6中，由血液微循環機理知加速脈搏波A-C段的時間tAC可以反映血液由心臟傳播到指尖末端毛細血管并反射匯合的時間，而這個時間可準確的表示脈搏波傳導時間。脈搏波傳導時間的提取的關鍵是找到特征點A和C的位置。其中，特征點A是波形的極大值點，準確地識別該位置較為容易。故特征點C的位置識別成為脈搏波傳導時間提取的關鍵。特征點C位置出現有兩種情況：一種是特征點C位于B和Z之間，如圖7(a)所示；另一種是特征點C位于Z和M之間，如圖7(b)所示。

圖7 特征點C位于不同位置的加速脈搏波
Fig.7 Acceleration pulse wave with different feature points

本文首先對經數字信號處理后的PPG信號的加速脈搏波進行分類[13]，再提取兩種不同類型的PPG信號的傳導時間，進而分別建立兩種PPG信號(特征點C位于B和Z之間或Z和M之間)的血壓計算模型。

2.1.1 加速脈搏波波形識別

為識別出如圖7中的兩種不同的加速脈搏波，本文利用迭代算法Adaboost進行實驗。 Adaboost在訓練過程中自適應調整訓練樣本權重，強化學習難以區分的訓練樣本，由加權投票機制聯合弱分類器，使分類效果好的弱分類器的權重得到提高，而分類效果差的弱分類器的權重降低，形成高性能的強分類器。算法描述如下：

①輸入訓練集D：

D={(x1,y1),(x2,y2),…,(xn,yn)},

其中，xi∈X,yi∈Y,X為訓練樣本集，Y為分類類別標志，Y={-1,+1}。

②初始化權值：D1(x)=1/n。

③設置訓練輪數T：fort=1,2，…，T。

④將弱學習算法在權值Dt下訓練，得到預測函數ht：X～{-1,+1}。

⑦訓練完畢，輸出函數：

為驗證Adaboost算法識別不同脈搏波的可行性，本文由預測靈敏度SE、陽性率TPR以及準確率AC這3個指標來評價該算法區分兩種脈搏波的性能。其中，

(10)

(11)

(12)

其中，TP和TN分別表示正確分類了第一類和第二類脈搏事件，FP和FN分別表示錯誤分類了第一類和第二類脈搏事件。

對250個波形信號以Adaboost算法進行識別，其中第一類(特征點C位于B和Z之間)130個，第二類(特征點C位于Z和M之間)120個，得到測量數據如表2所示。由表1知，以特征點A-Z的時間段tAZ作為分類依據，采用Adaboost算法可以準確的區分出兩種不同的加速脈搏波。

表2 分類評估參數
Tab.2 Evaluation parameters

TPTNFPFN靈敏度陽性率準確率11910411160.8810.9150.892

2.1.2 特征點提取

特征點A是加速脈搏波的波峰，在波形任意位置開始，在一個周期范圍內判斷出幅值最大的位置即為特征點A的位置。

特征點C的位置提取分兩種情況：①當C介于B與Z之間時，如B與Z間存在極大值，則極大值坐標即為C的位置；如沒有極大值，則以B與Z之間斜率最小的點的坐標作為C的位置；②當C介于Z與M之間時，若B、M兩點所在的直線斜率為k，以Z、M間所有點的斜率與k差值最小時的坐標作為C的位置。準確判斷特征點A和C的位置后，A與C之間的采樣點數與采樣時間的乘積即脈搏波傳導時間(PWTT)。

2.2 血壓模型的建立

準確識別出兩種不同的脈搏波并獲取脈搏波傳導時間后，由脈搏波傳導時間與血壓間的線性關系，進行回歸分析得出回歸系數及常數，建立血壓測量回歸方程式。經回歸分析以及雙彈性腔模型理論得出C點位于B和Z之間的血壓方程如式(13)、(14)所示，C點位于Z與M之間血壓方程如式(15)、(16)所示，其中，Ps為收縮壓，Pd為舒張壓。

Ps=348.852 1-1.389 2*PWTT，

(13)

Pd=20.883 6+0.402 6*Ps，

(14)

Ps=359.532 6-1.523 1*PWTT，

(15)

Pd=5.235 5+0.602 1*Ps。

(16)

3 實驗分析

為驗證所建立的血壓估算模型的可靠性，本研究采用水銀血壓計與本系統對測試者同時進行檢測以分析該設計的準確性。其中，以水銀血壓計測得數據作為真值，本設計測量數據為測值。在室溫25 ℃下，對年齡(18～65)及性別各不相同的150位測試者(包括110名血壓正常的測試者和40名高血壓患者)進行實驗，測量得到的數據及相應的誤差圖如圖8、圖9所示。部分測量數據如表3所示(表中測值和真值單位為mmHg)。

由圖8、圖9及表3知，本系統測得血壓值與真值之間均保持較好的一致性，數據誤差均保持在±5以內，優于AAMI標準。兩種加速脈搏波的分類具有較高準確性，具有一定的臨床意義。

表3 實驗測量數據
Tab.3 Measured data

序號性別/年齡 PS測值PS真值誤差PD測值PD真值誤差分類結果真實類別1M/209899-16063-3112M/25102107-56263-1113F/199695+16265-3224M/29110107+36564+1115F/27116119-37473+1226M/33118116+27580-5227F/33120125-57781-4228F/35123121+27877+1219M/27106110-46265-31110F/28108109-16972-32111M/39136138-27574+11112M/45144139+59188+32213M/43139142-37681-51114F/40145142+39288+42215F/42147146+19395-222

圖8 采集數據分布圖
Fig.8 Distribution graph of collected data

圖9 測值和真值誤差圖
Fig.9 Error between measure data and true data

4 結論

結合EMD及形態學方法對光電容積脈搏波進行降噪處理，對比不同濾波方法的信噪比，發現本文方法濾波效果優于小波變換。為了驗證Adaboost識別不同脈搏波的有效性，本文基于預測靈敏度、陽性率及準確率3個指標進行評價，實驗表明Adaboost識別不同脈搏波具有一定的可靠性。為驗證血壓模型的有效性，對150位測試者進行測量，實驗結果證明，論文建立的血壓計算模型誤差在±5 mmHg以內，具有一定的準確性和可靠性，達到了較好的效果。