基于VMD-SWT聯合算法的FMCW雷達生命體征檢測

何鵬宇,卓智海

(北京信息科技大學 信息與通信工程學院,北京 100101)

0 引言

生命體征能夠反映人體的健康狀況,呼吸率和心率是重要的生命體征,可作為診斷某些疾病的重要參數。生命體征檢測技術可分為接觸式和非接觸式兩類。傳統的生命體征檢測方法主要采用接觸式監測呼吸和心率[1-2],但接觸式設備會限制被測者的行動,使被測者感到不適。而非接觸式雷達檢測技術可以遠距離監測呼吸率和心率,符合人們對生命體征檢測技術智能化和舒適度的要求,更適合于長期監測,便于健康管理。

目前用于非接觸式生命體征檢測的雷達系統主要有三種,分別是連續波(continuous wave,CW)雷達、超寬帶(ultra-wideband,UWB)雷達和調頻連續波(frequency modulated continuous wave,FMCW)雷達[3-5],文獻[6]比較了三種雷達系統的區別。FMCW雷達的發射頻率是隨時間線性調制的,可以在探測目標距離的同時,利用多普勒效應檢測胸腔位移。通過提取中頻信號的頻率及相位,可以獲得被測對象的準確距離和位移信息[3-4]。與微波雷達相比,毫米波雷達的帶寬更寬、波長更短,因此具有更高的距離分辨率和靈敏度,適合于非接觸式生命體征檢測[7]。FMCW毫米波雷達在距離測量和微小位移測量方面都具有優異的性能[5],因此本文采用77 GHz FMCW毫米波雷達進行非接觸式生命體征檢測。

在處理生命體征信號的過程中,由于呼吸引起的胸部位移遠大于心跳引起的胸部位移,較弱的心跳信號容易受到呼吸諧波和雜波的干擾,特別是當心跳信號和呼吸諧波在頻域上接近或重疊時,呼吸諧波的頻率可能被誤認為心率。因此,如何處理呼吸信號引起的諧波干擾是心率檢測的主要問題[8-9]。文獻[10]提出了一種基于經驗模態分解(empirical mode decomposition,EMD)的生命體征檢測算法。然而,EMD的模態混合和端點效應問題導致心跳和呼吸諧波在一個模態中共存,影響了心跳信號檢測的準確性[11]。文獻[12]提出了集合經驗模態分解(ensemble empirical mode decomposition,EEMD),通過填充高斯白噪聲來提高EMD的性能,但模態分量中的殘留噪聲會降低心率估計的準確性。文獻[13]提出了基于變分模態分解(variational mode decomposition,VMD)的生命體征檢測,并通過實驗驗證了在相同的測試環境下,VMD的效果要優于EMD。但該方法不能提供心跳信號的時頻信息,當心率在短時間內出現波動時,檢測結果的誤差較大。文獻[14]提到同步壓縮變換在醫學領域中用于對腦電圖或心電圖的特征提取。其中同步壓縮小波變換(synchrosqueezed wavelet transform,SWT)是常用的同步壓縮變換之一,能夠提供高分辨的時頻圖,文獻[15-16]將其用于心電圖中的心跳信號時頻特征提取。本文將SWT引入了雷達生命體征檢測方法中,提出了VMD和SWT的聯合算法,經過實驗驗證,該方法能有效分離呼吸信號和心跳信號,降低呼吸諧波對心率估計的干擾,得到高精度的心跳信號時頻曲線。

1 FMCW雷達生命體征檢測原理

FMCW雷達發射線性調頻信號,正交接收機捕獲目標反射的回波信號并與發射信號進行正交混頻得到中頻信號。提取中頻信號的頻率和相位可以確定目標與雷達的距離和目標的胸部位移。發射的線性調頻信號可以表示為

(1)

式中:AT為信號的幅度;fc為線性調頻信號的起始頻率;B為帶寬;Tc為線性調頻信號的持續時間;φ(t)為相位噪聲。

接收信號可以認為是發射信號的延遲。若d0為目標到雷達的距離,x(t)為胸部運動引起的位移,則胸部到雷達的距離為R(t)=d0+x(t)。延遲時間為td=2R(t)/c,其中c為光速。則接收到的信號可以表示為

xR(t)=

(2)

式中:AR為接收信號的幅度。

接收信號和發射信號通過兩個正交I/Q信道進行混合,然后通過低通濾波器得到中頻信號y(t),其可以表示為

y(t)=

(3)

y(t)=ATARexp(j(2πfbt+φb(t)))

(4)

式中:fb=2Bd0/(cTc);φb(t)=4π(d0+x(t))/λc,λc為頻率fc對應的波長。

對于身體健康的成年人,心跳的位移幅度為0.1～0.5 mm,呼吸的位移幅度為1～12 mm,心跳頻率為0.8～2 Hz,呼吸頻率為0.1～0.5 Hz。在相位φb中,胸部位移x(t)和波長λc都處于毫米范圍,因此微小的胸部位移也會引起相位φb(t)的明顯變化,通過持續發射線性調頻信號并提取中頻信號相位可以實現對胸部位移x(t)采樣,進而估計呼吸率和心率。

2 VMD-SWT聯合算法

2.1 概述

基于VMD-SWT聯合算法的FMCW雷達生命體征檢測流程如圖1所示,主要包括目標定位、相位提取、信號分離和頻率估計4個步驟。首先,對中頻信號的采樣數據進行距離快速傅里葉變換(fast Fourier transform,FFT)以獲得距離特征(range profile)。然后,根據距離特征的極值確定目標對應的距離單元(range bin)。提取目標所在的距離單元的相位值進行相位展開和相位差分運算。之后對相位差分信號使用VMD算法進行信號分解,得到呼吸信號和心跳信號,通過搜索這兩個信號的頻譜峰值可以初步估計呼吸率和心率。但心跳信號易受呼吸諧波的干擾,且該方法不能得到心跳信號的時頻信息,因此需要使用 SWT進一步處理心跳信號,并采用模極大值法提取小波脊線得到心跳信號的時頻曲線,提高對心跳信號的檢測精度。

2.2 基于變分模態分解的信號分解方法

變分模態分解算法[17]假設任何信號都是由一系列具有特定中心頻率和有限帶寬的子信號組成,即本征模態函數(intrinsic mode function,IMF)。IMF定義為調幅調頻信號,可以表示為

uk(t)=Ak(t)cos(φk(t))

(5)

式中:Ak(t)和φk(t)分別為第k個模態uk(t)的幅值和相位,幅值Ak(t)和瞬時相位φk′(t)的變化都比較小,且模態uk(t)具有中心頻率ωk。模態uk(t)和其中心頻率ωk可通過變分問題求解,求解過程需要滿足兩個條件:①模態帶寬之和最小;②模態之和等于原始信號。假設原始信號f(t)被分解成k(正整數)個模態,根據上述條件,約束變分表達式為

(6)

式中:δ(t)為單位沖激函數;*表示卷積。引入拉格朗日乘數λ和懲罰因子α,將約束問題變為非約束問題,使約束最小化問題轉為求增廣拉格朗日函數L的鞍點,其表達式為

L({uk},{ωk},λ)=

(7)

uk(t)和ωk可以通過交替方向乘子法迭代得到,迭代公式如下:

(8)

(9)

(10)

(11)

(12)

式中:εr和εa分別為相對誤差和絕對誤差。當兩式均滿足時,停止迭代。

對相位信號使用VMD方法分解,可獲得一組模態分量,根據模態分量的中心頻率ωk確定所對應的信號,實現呼吸信號和心跳信號的分離。

2.3 基于同步壓縮小波變換的時頻分析

同步壓縮小波變換以連續小波變換為基礎,給定母小波函數ψ(t),對信號s(t)進行連續小波變換,其小波變換系數Ws(a,b)定義為

(13)

根據Plancherel定理,信號s(t)以ψ(t)為母小波函數的連續小波變換Ws(a,b)可以重寫為

(14)

(15)

(16)

通過模極大值法可提取同步壓縮小波變換系數Ts(ωl,b)的小波脊線。

使用 SWT對VMD分解結果中的心跳信號分量進行處理,并提取小波脊線即可得到心跳信號的時頻曲線。

3 實驗結果與分析

本文采用德州儀器的毫米波雷達IWR1443,一個幀周期發射一個啁啾(Chirp)信號,將接收信號與發射信號混頻得到中頻信號,對中頻信號采樣并進行數據處理。參數配置如表1所示。

表1 參數設置Table 1 Parameter configuration

被測者坐在雷達傳感器前,胸部與雷達處于同一水平高度且沒有雜物干擾,距離雷達0.3～2 m范圍內。實驗設備及實驗場景如圖2所示。同時被測者使用指夾式脈搏血氧儀測量心跳作為參考值。

圖2 實驗設備及實驗場景Fig.2 Experimental equipment and experimental setup

根據實測數據繪制的距離特征如圖3所示。圖3顯示了1 min內測量的1 200幀數據,目標的信息反映在圖中較亮的部分。從圖中可以看出,目標大約在距離雷達0.7 m的位置。

圖3 距離特征Fig.3 Range profile

根據距離特征的極值確定目標所在的距離單元并進行相位提取、相位展開和相位差分,結果如圖4所示。圖4顯示了30 s內目標的胸部位移信息,圖中頻率低幅度高的分量由呼吸引起,而頻率高幅度低的分量主要由心跳引起,同時也存在呼吸信號的諧波以及其他干擾,相位信號的頻譜如圖5所示。從圖5可以看出在0.25 Hz和0.50 Hz處有明顯的峰值,分別對應呼吸率和呼吸諧波頻率。而心跳頻率的峰值不明顯,通過求取正常心率范圍0.8～2.0 Hz之間的極值,得到峰值為1.17 Hz。

圖4 相位信號Fig.4 Phase signal

圖5 相位信號頻譜Fig.5 Phase signal spectrum

對圖4相位信號使用EMD算法進行分解,得到5個分量IMF1～IMF5以及殘余分量res,其分解結果如圖6左列所示,對每個分量進行傅里葉變換,得到其對應的頻譜,結果如圖6右列所示。從圖中可以看出,IMF3為呼吸信號,IMF2為呼吸諧波信號與心跳信號。EMD的分解模態數不能預先設定,會一直迭代到滿足終止條件,模態數不確定導致呼吸信號與心跳信號所處的分量也不確定,對分解后的處理十分不便,且分解結果也不夠理想,心跳和呼吸諧波在一個模態中共存。

而VMD算法可以提前設定模態數量,將其預先設定為4,其分解結果如圖7左列所示,每個分量對應的頻譜如圖7右列所示。從圖中可以看出,IMF4對應呼吸信號,IMF3對應呼吸諧波信號,IMF2對應心跳信號。與EMD相比,VMD的分解結果更接近正弦波形,其頻譜峰值明顯,在提取心跳信號方面明顯優于EMD。在IMF2中可以看到2個較明顯的峰值分別為1.00 Hz和1.17 Hz,出現兩個峰值的原因可能是呼吸信號的諧波干擾,或者是心跳信號的頻率在測量時間內有波動,其在兩個峰值頻率所保持的時間較長。由于VMD的分解結果不能提供時頻信息,無法得知心跳信號在測試時間內的變化趨勢,因此需要對心跳信號使用SWT進一步處理。

圖6 EMD模態分量及頻譜Fig.6 EMD mode components and spectrum

圖7 VMD模態分量及頻譜Fig.7 VMD mode components and spectrum

對VMD的模態IMF2使用SWT進行時頻分析,其結果如圖8所示。通過模極大值法提取小波脊線,可得到心跳信號的時頻曲線,其結果如圖9所示。

圖8 心跳信號的SWT時頻面表示Fig.8 SWT time-frequency representation of heartbeat signal

圖9 心跳信號的時頻曲線Fig.9 Time-frequency curve of heartbeat signal

將圖9的心跳頻率轉換為每分鐘心跳次數,指夾式脈搏血氧儀測量的心率作為參考值ref,兩者對比的結果如圖10所示。從圖中可以看出,兩條曲線雖然不能完全重合,但其變化趨勢基本一致,且SWT的結果與參考值偏差也較小。因此VMD-SWT聯合算法可以獲得準確的心跳信號時頻曲線。

圖10 SWT時頻曲線與指夾式脈搏血氧儀檢測結果對比Fig.10 Comparison of SWT time-frequency curves with finger pulse oximeter detection results

為了進一步驗證所提算法的準確性和有效性,對10名測試者進行了實驗測試。每位測試者采集3 min的數據,并從中提取出1 min的數據長度進行呼吸率和心率估計,并使用均方根誤差(root mean square error,RMSE)對心率估計進行性能比較。實驗結果如表2所示,其中EMD和VMD的結果為心跳信號的頻譜峰值所對應的心率,指壓式脈搏血氧儀和VMD-SWT的聯合算法為1 min內的心率平均值。由表2可以看出,VMD-SWT聯合算法的準確性高于其他方法。

表2 VMD-SWT聯合算法與其他方法的心率檢測結果對比Table 2 Comparison of heart rate detection results between VMD-SWT joint algorithm and other methods 次/min

4 結束語

本文基于FMCW毫米波雷達實現生命體征監測,針對雷達生命體征檢測中心跳信號易受到呼吸諧波干擾的現象,提出了基于VMD-SWT的聯合算法,對VMD分解結果中的心跳信號使用SWT做進一步處理,獲取心跳信號的時頻信息,得到更準確的心率估計,并通過實驗將該方法與EMD和VMD方法進行對比。實驗結果表明,基于VMD-SWT的聯合算法能夠準確有效地分離呼吸信號和心跳信號,與EMD和VMD方法相比,該方法可以提供準確的時頻信息,降低呼吸諧波對檢測結果的干擾,提高雷達生命體征檢測精度和穩定性。