基于諧波陷波器改進(jìn)HMLD 的呼吸心跳分離算法

汪新坤，曹 樂，張文艷

（上海工程技術(shù)大學(xué) 電子電氣工程學(xué)院，上海 201620）

0 引言

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

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

1 基于毫米波雷達(dá)的生命體征檢測(cè)原理

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

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

存在多個(gè)傳播信道時(shí)，接收信號(hào)可由公式（2）表示［8］。

其中，t，τ分別表示雷達(dá)回波數(shù)據(jù)中的慢時(shí)間和快時(shí)間，Akδ（τ -τk（t））和分別表示胸廓運(yùn)動(dòng)和靜態(tài)目標(biāo)對(duì)快時(shí)間和振幅峰值的響應(yīng)。τk（t）由d（t）決定，可由公式（3）表示。

其中，c表示為光速。

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

其中，m，n分別以慢時(shí)間和快時(shí)間顯示采樣數(shù)；Tf是快時(shí)間采樣間隔，對(duì)應(yīng)系統(tǒng)ADC 采樣速率；Ts是脈沖持續(xù)時(shí)間，對(duì)應(yīng)系統(tǒng)信號(hào)采樣頻率。

雷達(dá)回波矩陣包含目標(biāo)微動(dòng)位移信息，Δφ目標(biāo)信號(hào)前后相位差由公式（5）表示［9］。

其中，Δx為胸腔運(yùn)動(dòng)位移，λ為毫米波對(duì)應(yīng)的波長。

根據(jù)公式（5）對(duì)雷達(dá)回波矩陣進(jìn)行數(shù)據(jù)解析，提取胸廓微動(dòng)信號(hào)。

2 基于諧波陷波器改進(jìn)HMLD 的呼吸心跳分離

2.1 雷達(dá)回波預(yù)處理

為了準(zhǔn)確提取人體胸廓目標(biāo)微動(dòng)信號(hào)，需要對(duì)雷達(dá)回波矩陣進(jìn)行相應(yīng)的數(shù)據(jù)預(yù)處理。預(yù)處理的算法流程圖如圖1 所示。

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

因?yàn)槔走_(dá)天線接收到物體反射的電磁波能量最大，反映在雷達(dá)回波的數(shù)據(jù)中即信號(hào)幅值最大。所以在雷達(dá)矩陣的快時(shí)間軸上檢測(cè)信號(hào)幅值最大的所在的距離單元，即目標(biāo)所在的位置，可由公式（6）表示［10］。

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

根據(jù)幅值最大檢測(cè)原理，在慢時(shí)間軸上依次尋找每個(gè)脈沖的目標(biāo)距離單元，并提取該距離單元目標(biāo)的相位信號(hào)，最后對(duì)提取的相位信號(hào)解析提取微動(dòng)信號(hào)，并進(jìn)行低通濾波和均值濾波，消除高頻和靜態(tài)噪聲。對(duì)原始的胸廓微動(dòng)信號(hào)進(jìn)行濾波處理，提取得到的信號(hào)如圖2 所示。

圖2 胸廓微動(dòng)信號(hào)提取、濾波后波形Fig.2 Thoracic micro motion signal extraction and filtered waveform

2.2 基于諧波陷波器改進(jìn)的HMLD 原理及方法

HMLD 算法對(duì)胸廓微動(dòng)信號(hào)在［0.1，0.8］Hz，［0.8，2.0］Hz 呼吸和心跳檢測(cè)頻段分別進(jìn)行帶通濾波，并在對(duì)應(yīng)的頻譜中循環(huán)查找信號(hào)基頻和其他高次諧波頻率，進(jìn)而判斷呼吸心跳信號(hào)基波頻率是否存在，信號(hào)基頻與諧波頻率fN的關(guān)系可由公式（7）表示［11］。

其中，fbase為信號(hào)的基波頻率。

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

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

呼吸頻率提?。?/p>

步驟1對(duì)雷達(dá)回波信號(hào)預(yù)處理，得到胸腔運(yùn)動(dòng)信號(hào)x（t）；

步驟2在呼吸檢測(cè)頻段對(duì)x（t）進(jìn)行帶通濾波和快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT），獲取呼吸檢測(cè)頻段的頻譜信號(hào)；

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

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

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

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

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

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

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

心跳頻率提?。?/p>

步驟9在呼吸頻率提取的基礎(chǔ)上，胸腔運(yùn)動(dòng)信號(hào)x（t）通過陷波濾波器將呼吸的基波分量和諧波分量濾除得到x′（t）；

步驟10在心跳信號(hào)檢測(cè)頻段上對(duì)x′（t）進(jìn)行帶通濾波和快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT），獲取心跳檢測(cè)頻段的頻譜信號(hào)；

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

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

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

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

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

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

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

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

2.3 諧波陷波器設(shè)計(jì)

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

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

2.4 基于連續(xù)小波逆變換的信號(hào)重構(gòu)

連續(xù)小波逆變換可以在信號(hào)的任意的時(shí)頻位置對(duì)小波系數(shù)進(jìn)行修改［12］。所以可以對(duì)提取的人體呼吸心跳頻率在原始的胸腔運(yùn)動(dòng)信號(hào)中重構(gòu)呼吸和心跳的時(shí)域信號(hào)。本文采用Morlet 小波基函數(shù)對(duì)呼吸和心跳信號(hào)進(jìn)行分析，因?yàn)镸orlet 小波具有良好的時(shí)頻局部化特性且形狀與心跳信號(hào)相近［13］，小波變換公式（15）：

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

根據(jù)Morlet 小波轉(zhuǎn)換公式進(jìn)行呼吸和心跳信號(hào)重構(gòu)，具體步驟如下：

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

（2）根據(jù)頻率轉(zhuǎn)換尺度的關(guān)系，計(jì)算尺度參數(shù)a，公式（16）：

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

（3）根據(jù)不同尺度的Morlet 小波函數(shù)對(duì)胸腔運(yùn)動(dòng)信號(hào)x（t）進(jìn)行小波變換，得到小波系數(shù)矩陣C；

（4）通過系數(shù)矩陣C進(jìn)行小波逆變換，獲得重構(gòu)信號(hào)r（t）。

3 實(shí)驗(yàn)及結(jié)果分析

3.1 實(shí)驗(yàn)設(shè)計(jì)

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

圖5 實(shí)驗(yàn)場(chǎng)景圖Fig.5 Experimental scene

3.2 呼吸心跳分離

選取測(cè)試對(duì)象1 的雷達(dá)回波數(shù)據(jù)對(duì)提取到的胸腔運(yùn)動(dòng)信號(hào)進(jìn)行呼吸心跳信號(hào)分離。根據(jù)諧波陷波器改進(jìn)的HMLD 算法，先對(duì)呼吸的檢測(cè)頻段進(jìn)行帶通濾波，并進(jìn)行FFT 得到呼吸檢測(cè)頻譜，如圖6 所示。

圖6 呼吸檢測(cè)頻譜Fig.6 Respiratory detection spectrum

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

圖7 心跳檢測(cè)頻譜Fig.7 Heartbeat detection spectrum

根據(jù)提取的呼吸頻率fR和心跳頻率fH，通過連續(xù)小波逆變換進(jìn)行呼吸和心跳時(shí)域信號(hào)重構(gòu)，結(jié)果如圖8 所示，可以觀察到重構(gòu)后的呼吸和心跳信號(hào)整體趨勢(shì)平穩(wěn)、光滑。

圖8 呼吸心跳時(shí)域信號(hào)重構(gòu)Fig.8 Time domain signal reconstruction of respiration and heartbeat

3.3 結(jié)果分析

本文采用變分經(jīng)驗(yàn)?zāi)B(tài)分解（Variational Mode Decomposition，VMD）算法對(duì)每名測(cè)試對(duì)象進(jìn)行呼吸心跳分離，并對(duì)提取的結(jié)果進(jìn)行信噪比（Single Nosie Ratio，SNR）比較，公式（17）：

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

基于諧波陷波器改進(jìn)的HMLD 算法與VMD 算法提取心跳的結(jié)果見表1。

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

表1 基于陷波器改進(jìn)HMLD 與VMD 算法的實(shí)驗(yàn)結(jié)果Tab.1 Experimental results of improved HMLD and VMD algorithm based on notch filter

4 結(jié)束語

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