基于諧波陷波器改進HMLD 的呼吸心跳分離算法

汪新坤，曹 樂，張文艷

（上海工程技術大學 電子電氣工程學院，上海 201620）

0 引言

基于毫米波雷達可以不受天氣，光線等一些環境因素的干擾進行非接觸式生命體征檢測，對醫療監護和智能家居健康監測具有重要意義和發展前景［1］。但是，基于毫米波雷達的生命體征信號檢測仍存在一些困難和挑戰，例如測試環境中的靜態雜波、人為肢體運動偽影、呼吸諧波等干擾因素都會影響呼吸和心跳信號的提?。?］。一般可以通過濾波處理削弱雜波和噪聲干擾，但強呼吸次諧波有可能在心跳有效的檢測頻率范圍內，且心跳信號運動幅度更小，對心跳信號提取造成一定的困難，因此強呼吸次諧波的干擾已成為毫米波雷達體征信號檢測中亟待解決的問題。

為了克服呼吸諧波對心跳信號提取的影響，Van Nguyen 團隊［3-4］前后提出一種諧波路徑算法（Harmonic Path Algorithm，HAPA）和平均諧波路徑算法（Spectrum-Averaged Harmonic Path，SHAPA），但這兩種算法存在無法找到諧波路徑的問題，并且算法容易受閾值設置的影響，算法效率較低；Zhang［5］提 出 一 種 諧 波 倍 數 循 環 檢 測 的 算 法（Harmonic Multiple Loop Detection，HMLD）算法，但該算法只計算到呼吸的二次諧波，未考慮到呼吸的三次諧波帶來的影響，若身體發生與呼吸同頻的噪聲，會增強呼吸諧波的強度，繼而影響心跳頻率的估計；張怡［6］設置的呼吸和心跳信號檢測頻段被錯開，存在頻段檢測約束問題，所以該算法沒有從真正意義上消除呼吸二次諧波及三次諧波帶來的影響。為了解決這一問題，本文提出一種基于陷波器諧波消除改進的HMLD 算法，通過陷波濾波器將鄰近心跳基波的呼吸諧波濾除，實現心跳基波頻率的精確提?。徊捎眠B續小波逆變換對已提取到的呼吸和心跳頻率重構時域信號；最后，通過靜態人體實驗進行分析，驗證了基于諧波陷波器改進的HMLD 算法能夠有效的將呼吸和心跳信號分離和提取。

1 基于毫米波雷達的生命體征檢測原理

毫米波雷達發射電磁脈沖，當接觸到人體胸廓時，由于身體的高反射率，發射的脈沖信號被反射并被接收。人體的呼吸和心跳引起胸廓周期性的擴張和收縮，雷達發射天線與胸廓的距離d（t）可由公式（1）表示［7］。

其中，d0為接收天線和人體胸廓的初始距離；dr和dh分別為呼吸信號幅度和心跳信號幅度；fr和fh分別為呼吸頻率和心跳信號頻率。

存在多個傳播信道時，接收信號可由公式（2）表示［8］。

其中，t，τ分別表示雷達回波數據中的慢時間和快時間，Akδ（τ -τk（t））和分別表示胸廓運動和靜態目標對快時間和振幅峰值的響應。τk（t）由d（t）決定，可由公式（3）表示。

其中，c表示為光速。

雷達將接收到的信號離散化為二維雷達回波矩陣，可由公式（4）表示。

其中，m，n分別以慢時間和快時間顯示采樣數；Tf是快時間采樣間隔，對應系統ADC 采樣速率；Ts是脈沖持續時間，對應系統信號采樣頻率。

雷達回波矩陣包含目標微動位移信息，Δφ目標信號前后相位差由公式（5）表示［9］。

其中，Δx為胸腔運動位移，λ為毫米波對應的波長。

根據公式（5）對雷達回波矩陣進行數據解析，提取胸廓微動信號。

2 基于諧波陷波器改進HMLD 的呼吸心跳分離

2.1 雷達回波預處理

為了準確提取人體胸廓目標微動信號，需要對雷達回波矩陣進行相應的數據預處理。預處理的算法流程圖如圖1 所示。

圖1 預處理算法流程圖Fig.1 Flow chart of preprocessing algorithm

因為雷達天線接收到物體反射的電磁波能量最大，反映在雷達回波的數據中即信號幅值最大。所以在雷達矩陣的快時間軸上檢測信號幅值最大的所在的距離單元，即目標所在的位置，可由公式（6）表示［10］。

其中，Δd表示距離單元跨度；n表示所在的距離單元；d0表示初始距離。

根據幅值最大檢測原理，在慢時間軸上依次尋找每個脈沖的目標距離單元，并提取該距離單元目標的相位信號，最后對提取的相位信號解析提取微動信號，并進行低通濾波和均值濾波，消除高頻和靜態噪聲。對原始的胸廓微動信號進行濾波處理，提取得到的信號如圖2 所示。

圖2 胸廓微動信號提取、濾波后波形Fig.2 Thoracic micro motion signal extraction and filtered waveform

2.2 基于諧波陷波器改進的HMLD 原理及方法

HMLD 算法對胸廓微動信號在［0.1，0.8］Hz，［0.8，2.0］Hz 呼吸和心跳檢測頻段分別進行帶通濾波，并在對應的頻譜中循環查找信號基頻和其他高次諧波頻率，進而判斷呼吸心跳信號基波頻率是否存在，信號基頻與諧波頻率fN的關系可由公式（7）表示［11］。

其中，fbase為信號的基波頻率。

由于正常人的呼吸基頻在［0.1，0.5］Hz，心跳基頻在［0.8，2.0］Hz，若考慮到呼吸高次諧波的影響，HMLD 算法在心跳信號頻率中存在檢測效率和檢測頻段約束問題。而基于諧波陷波器改進的HMLD算法可以消除呼吸高次諧波的影響。為進一步提高HMLD 算法有效性，忽略呼吸四次諧波和心跳三次諧波的存在，將呼吸和心跳檢測頻段分別設置為［0.1，1.5］Hz，［0.8，4.0］Hz。

改進的算法主要分為兩部分：呼吸頻率提取，心跳頻率提取，其具體算法步驟如下。

呼吸頻率提?。?/p>

步驟1對雷達回波信號預處理，得到胸腔運動信號x（t）；

步驟2在呼吸檢測頻段對x（t）進行帶通濾波和快速傅里葉變換（Fast Fourier Transform，FFT），獲取呼吸檢測頻段的頻譜信號；

步驟3在頻譜圖中查找所有的峰值頻率；

步驟4在提取的所有峰值頻率中以最大峰值頻率作為呼吸信號第一條峰值路徑的待定基波頻率fR1_1；

步驟5根據公式（5），提取呼吸二次諧波頻率fR2_1和3 次諧波頻率fR3_1，式（8）和式（9）：

步驟6在頻譜圖中剔除第一條峰值頻率路徑fR1_1、fR2_1和fR3_1，并且尋找最大的峰值頻率作為呼吸信號第二條峰值路徑的待定基波頻率fR1_2；

步驟7重復步驟3～5，為提高檢測效率，防止無法找到待定基波的諧波頻率，最多查找3 條峰值頻率路徑Pn，如公式（10）。

步驟8計算每條峰值路徑的平均功率，并且以最大的平均功率的峰值路徑的待定基波頻率估計呼吸頻率fR，如公式（11）。

其中，p表示最大平均功率峰值路徑的下標。

心跳頻率提取：

步驟9在呼吸頻率提取的基礎上，胸腔運動信號x（t）通過陷波濾波器將呼吸的基波分量和諧波分量濾除得到x′（t）；

步驟10在心跳信號檢測頻段上對x′（t）進行帶通濾波和快速傅里葉變換（Fast Fourier Transform，FFT），獲取心跳檢測頻段的頻譜信號；

步驟11在心跳頻譜圖中查找所有的峰值頻率；

步驟12在提取的所有的峰值頻率中以最大峰值頻率作為心跳信號第一條峰值路徑的待定基波頻率fH1_1；

步驟13根據公式（5），提取心跳二次諧波頻率fH2_1為

步驟14在頻譜圖中剔除第一條峰值頻率路徑fH1_1和fH2_1，并且尋找最大的峰值頻率作為呼吸信號第二條峰值路徑的待定基波頻率fH1_2；

步驟15重復步驟10～12，與呼吸檢測相同，最多查找兩條峰值頻率路徑，為

步驟16計算每條峰值路徑的平均功率，并且以最大的平均功率的峰值路徑的待定基波頻率估計心跳頻率fH，即式（14）：

基于諧波陷波器改進的HMLD 算法流程圖如圖3 所示。

圖3 基于諧波陷波器改進的HMLD 算法流程圖Fig.3 Flow chart of improved HMLD algorithm based on harmonic notch filter

2.3 諧波陷波器設計

本文數據基于MATLAB 軟件平臺處理的，所以陷波濾波器的傳遞函數分子和分母系數可以通過MATLAB 的庫函數獲取。在實際的信號處理過程中，有效的數據頻率在［0.1，2］Hz，而且數據在進行FFT 后存在頻譜泄露問題，所以對陷波器的參數設置有一定的要求。經過多次的實驗驗證，本文中將陷波濾波器階數設置為8 階，濾波器的帶寬設置為0.04 Hz，并且利用MATLAB 軟件上繪制出1.0 Hz 和0.3 Hz 的陷波器幅值響應曲線如圖4 所示，可以看到該陷波器陷波深度達到了100 dB，能夠實現對單一頻率信號濾除，且對頻率周圍其他信號影響較小。

圖4 1.0 Hz 及0.3 Hz 陷波器幅值響應曲線Fig.4 Amplitude response curve of 1.0 Hz and 0.3 Hz notch filter

2.4 基于連續小波逆變換的信號重構

連續小波逆變換可以在信號的任意的時頻位置對小波系數進行修改［12］。所以可以對提取的人體呼吸心跳頻率在原始的胸腔運動信號中重構呼吸和心跳的時域信號。本文采用Morlet 小波基函數對呼吸和心跳信號進行分析，因為Morlet 小波具有良好的時頻局部化特性且形狀與心跳信號相近［13］，小波變換公式（15）：

其中，a表示尺度參數；b為平移參數；w0為Morlet 母小波中心頻率。

根據Morlet 小波轉換公式進行呼吸和心跳信號重構，具體步驟如下：

（1）根據諧波消除改進的HMLD 算法提取的呼吸或心跳頻率f，設置小波變換的頻率范圍［f -BW，f ＋BW］，BW為陷波濾波器的帶寬；

（2）根據頻率轉換尺度的關系，計算尺度參數a，公式（16）：

其中，w0為Morlet 母小波中心頻率，fs為信號采樣頻率，f∈［f -BW，f ＋BW］，將頻率范圍等間隔離散化成［f1，f2，…，fn］，即可獲得尺度序列S ＝［a1，a2，…，an］。

（3）根據不同尺度的Morlet 小波函數對胸腔運動信號x（t）進行小波變換，得到小波系數矩陣C；

（4）通過系數矩陣C進行小波逆變換，獲得重構信號r（t）。

3 實驗及結果分析

3.1 實驗設計

本文采用自主設計的單通道的雷達信號采集系統進行實驗信號采集。系統由模擬前端、主控模塊及信號調理電路組成。其中模擬前端采用脈沖相干式A111 毫米波雷達芯，該芯片工作頻率為60 GHz，波長λ＝5 mm，檢測范圍為［0.2，0.8］m，距離單元跨度Δx ＝0.48 mm，系統的采樣率為50 Hz。實驗選取3 名平均年齡為24 歲的成年人作為實驗對象，在初始位置d0＝0.6 m 采集原始雷達回波信號，并用小米手環的實時監測心率值作為實驗的心跳參考信號，每次采集時間為2 min。實驗場景如圖5 所示。

圖5 實驗場景圖Fig.5 Experimental scene

3.2 呼吸心跳分離

選取測試對象1 的雷達回波數據對提取到的胸腔運動信號進行呼吸心跳信號分離。根據諧波陷波器改進的HMLD 算法，先對呼吸的檢測頻段進行帶通濾波，并進行FFT 得到呼吸檢測頻譜，如圖6 所示。

圖6 呼吸檢測頻譜Fig.6 Respiratory detection spectrum

在呼吸檢測頻譜中查找平均功率最大的呼吸峰值路徑為P ＝［0.300 5，0.608 40，0.908 5］，因此此時呼吸頻率為峰值路徑下的基波頻率fR ＝0.300 5 Hz。對于心跳頻率的提取，先將呼吸的基波頻率、二次諧波和三次諧波消除掉，然后在心跳檢測頻段進行帶通濾波和FFT，得到心跳信號檢測頻譜圖如圖7 所示。

圖7 心跳檢測頻譜Fig.7 Heartbeat detection spectrum

根據提取的呼吸頻率fR和心跳頻率fH，通過連續小波逆變換進行呼吸和心跳時域信號重構，結果如圖8 所示，可以觀察到重構后的呼吸和心跳信號整體趨勢平穩、光滑。

圖8 呼吸心跳時域信號重構Fig.8 Time domain signal reconstruction of respiration and heartbeat

3.3 結果分析

本文采用變分經驗模態分解（Variational Mode Decomposition，VMD）算法對每名測試對象進行呼吸心跳分離，并對提取的結果進行信噪比（Single Nosie Ratio，SNR）比較，公式（17）：

其中，s2（l）表示提取信號頻譜峰值的功率，∑s2（f）為信號頻譜的功率和。

基于諧波陷波器改進的HMLD 算法與VMD 算法提取心跳的結果見表1。

通過表1 可知，3 名測試對象采用改進的HMLD 算法提取心跳頻率誤差率分別為4.61%、11.5%、7.78%，提取心跳信號信噪比相比于VMD 平均提高了2.66 dB。，因此采用諧波陷波器改進的HMLD 算法提取的呼吸和心跳信號更為準確。

表1 基于陷波器改進HMLD 與VMD 算法的實驗結果Tab.1 Experimental results of improved HMLD and VMD algorithm based on notch filter

4 結束語

本文針對基于毫米波雷達生命體征信號檢測中存在呼吸諧波導致呼吸心跳難以分離的問題，提出了一種基于諧波陷波器改進的HMLD 算法在信號頻譜中提取呼吸和心跳頻率。首先，對3 名測試對象進行雷達回波預處理提取胸廓微動信號；其次，采用HMLD 算法提取呼吸頻率，采用諧波陷波器消除呼吸諧波，再根據HMLD 算法提取心跳頻率；最后，對提取的呼吸和心跳頻率進行連續小波逆變換，對時域信號進行重構。結果表明，提取的心跳頻率誤差率在11.5%以內，心跳信噪比相比于VMD 平均提高了2.66 dB，說明基于諧波陷波器改進HMLD 算法能夠準確的提取呼吸心跳信號，為基于毫米波雷達的高精度生命體征信號檢測提供了一種有效的研究方法。