基于77 GHz 毫米波雷達的非接觸式精確心跳檢測*

潘海鵬鄒永洋顧敏明

(浙江理工大學機械與自動控制學院，浙江 杭州 310018)

心跳檢測對于日常健康監(jiān)測和臨床醫(yī)療診斷有著重大作用。 一方面，心率能反映基本的健康狀態(tài)[1]；另一方面，心跳的節(jié)律和強度在很大程度上能體現(xiàn)內臟器官的病理變化[2]。 相比傳統(tǒng)的接觸式心跳檢測，如ECG、PPG、腕帶式血氧儀等，基于毫米波雷達的非接觸式心跳檢測可實現(xiàn)遠程監(jiān)護，在某些特殊場合有著顯著優(yōu)勢，如災后生命探測、嬰兒監(jiān)護、疲勞監(jiān)測等。

毫米波雷達具有體積小、空間分辨率高、大帶寬等優(yōu)點，毫米波波長較短，對胸壁位移這類微小運動具有更高靈敏度[3-4]。 Obeid 等人[5]通過2.4 GHz、5.8 GHz 和60 GHz 三種多普勒雷達對比，驗證了一個事實:雷達工作頻率越高，對小位移檢測靈敏度越高。 因此，77 GHz 毫米波雷達已被廣泛應用于車載雷達[6]、人體姿態(tài)識別[7]、體征檢測[8]等領域。

雷達心跳檢測一大挑戰(zhàn)在于呼吸及其諧波強干擾。 由于呼吸引起的胸壁位移(1 mm～12 mm)遠大于心跳引起的胸壁位移(0.1 mm ～0.5 mm)，導致心跳信號被淹沒。 Hu 等人[9]提出連續(xù)小波變換(CWT)與經(jīng)驗模態(tài)分解(EMD)相結合的方法，實現(xiàn)了心肺信號分離與提取，但該方法對信噪比要求較高，需預先基帶信號校準。 隨后，變分模態(tài)分解(VMD)[10]算法被提出，在變分框架中自適應分解信號，解決了EMD 模態(tài)混疊問題。 Duan 等人[11]將VMD 算法應用于雷達心跳檢測領域，實現(xiàn)了呼吸、心跳和噪聲的分離，但VMD 分解參數(shù)的選擇仍是一個亟待解決的問題。

雷達心跳檢測另一大挑戰(zhàn)在于頻譜分辨率不足。傳統(tǒng)心率估計通過采用超過10 s 的時間窗口來提高頻率分辨率[12-13]。 然而，時間窗口越長，丟失的心跳細節(jié)越多，實時性越差。 對于短時窗下的心率估計，Nosrati等人[14]首次將頻率時間相位回歸技術(FTPR)用于估計心率，顯著提高了頻率分辨率，但短時窗下會出現(xiàn)嚴重主頻估計偏差。 然而，實時心率估計往往只需要對心率范圍內的窄帶頻譜區(qū)間進行細微觀察與分析，頻譜細化算法[15-16]能對信號頻譜中感興趣的頻帶進行局部放大，有利于心率的準確估計。

本文提出了一種基于毫米波雷達的準確且快速的心跳檢測方法。 首先，對原始雷達回波信號進行預處理，實現(xiàn)心肺信號的初步分離。 然后，提出基于自適應變分模態(tài)分解(AVMD)的心跳信號提取算法，準確地恢復心跳信號。 進一步引入基于線性調頻Z 變換(CZT)的頻譜細化算法，用于估計實時心率。 最后，通過模擬駕駛行為實驗驗證本文提出的心跳檢測算法的有效性和魯棒性。

1 雷達體征檢測原理

連續(xù)波多普勒雷達體征檢測基本原理是捕捉微動目標引起的頻移，圖1 為77 GHz 毫米波雷達體征檢測系統(tǒng)模型。 本文將雷達放置于受試者正前方，距離受試者d0處，提取由呼吸和心跳引起的胸壁位移信號x(t)。 雷達發(fā)射信號可表示為:

圖1 77 GHz 毫米波雷達心跳檢測系統(tǒng)模型

式中:f為雷達工作頻率，φ(t)為初始相位。 首先雷達向受試者通過發(fā)射端發(fā)射T(t)。 隨后T(t)信號的相位被受試者胸壁運動x(t)所調制，反射信號被雷達接收端所捕獲，反射信號可表示為:

式中:c為無線電的波速，λ＝c/f表示雷達波長。 一般而言，R(t)可以近似看作經(jīng)過2d0/c時間延遲并相位調制后的T(t)。R(t)經(jīng)過前端的低噪聲放大器后，得到基帶信號B(t):

式中:θ＝4πd0/λ+θ0為常數(shù)相位偏移，與接收機本身參數(shù)以及d0相關；殘余相位噪聲簡化為Δφ(t)＝φ(t)-φ(t-2d0/c)。 然后，由正交混頻器解調基帶信號B(t)，同相和正交信號分別表示為:

經(jīng)過增益放大器后，由呼吸和心跳引起的胸壁信號x(t)由反正切解調和解卷繞技術計算可得，表達式如下:

最后，模擬信號x(t)經(jīng)過模數(shù)轉換器得到數(shù)字信號x(t)，用于數(shù)字信號處理以及心跳檢測算法的開發(fā)。

2 心跳檢測算法

為實現(xiàn)復雜環(huán)境下多普勒雷達對心跳信號的精確估計，本文提出了一種雷達體征數(shù)據(jù)處理的具體實現(xiàn)方案，主要包括三部分:雷達信號預處理、心跳信號提取、心跳頻譜估計。 雷達捕獲的受試者胸壁位移信號x(t)作為本文算法的輸入。

2.1 雷達信號預處理

預處理步驟包括去除背景噪聲和多項式趨勢、帶通濾波。 由于雷達測量范圍內除受試者胸壁運動外還存在大量靜止物體，通過減去原始雷達回波數(shù)據(jù)的平均值來去除靜止背景噪聲。 另外，由于熱噪聲和時間漂移引起的發(fā)射單元幅值的不穩(wěn)定導致回波數(shù)據(jù)呈現(xiàn)嚴重的多項式趨勢，通過減去最小二乘法擬合的曲線來去除多項式趨勢。

此外，基于呼吸頻率典型范圍為0.1 Hz～0.3 Hz、心率范圍為1 Hz～3 Hz 這一事實，帶通濾波器選取巴特沃斯濾波器，其通帶范圍設置為0.8 Hz～4.0 Hz，并應用于去除噪聲后的x(t)，得到預處理后的信號~x(t)。

2.2 心跳信號提取

2.2.1 變分模態(tài)分解

變分模態(tài)分解(VMD)作為一種完全非遞歸的自適應信號處理算法，通過迭代計算每個本征模態(tài)分量的中心頻率和帶寬，以每個模態(tài)估計帶寬之和最小為準測來確定變分模態(tài)模型的最優(yōu)解，自適應分解信號。

VMD 算法中，本征模函數(shù)(IMF)被定義為一個調幅-調頻信號，其表達式為:

式中:Ak(t)、φk(t)分別為瞬時幅值和相位，則uk(t)的瞬時頻率wk(t)＝dφ(t)/dt。

為獲取IMF 分量，VMD 摒棄EMD 的循環(huán)篩選的方式，而在變分框架下分解。 在變分模型的迭代求解過程中，更新各個IMF 分量的中心頻率和帶寬，根據(jù)信號頻域特征自適應分解信號頻帶，最終得到窄帶IMF 分量。 假設將原始信號~x(t)分解為K個IMF 分量，相應的約束變分模型可表示為:

式中:uk表示VMD 分解的K個IMF 分量；wk表示IMF 分量的中心頻率。

為求解上述約束變分問題最優(yōu)解，引入二次懲罰因子α和拉格朗日算子λ構造增廣拉格朗日函數(shù):

最后利用交替方向乘子算法(ADMM)求解上述增廣拉格朗日函數(shù)的鞍點，得到約束模型最優(yōu)解，然后得到~(t)的分解結果。

2.2.2 自適應變分模態(tài)分解

分解的模態(tài)數(shù)K和二次懲罰因子α的選擇控制原始信號和分解信號的重構誤差，參數(shù)的經(jīng)驗選擇可能導致較高的重構誤差，甚至丟失原始信號的某些重要信息。K過小，會出現(xiàn)模態(tài)混疊的現(xiàn)象；K過大，會出現(xiàn)過分解，產(chǎn)生虛假分量。α越小，IMF分量帶寬越大，越容易模態(tài)混疊；α越大，IMF 分量帶寬越小，避免模態(tài)混疊，但計算量增大[17]。 實際應用中為采樣頻率，一般可取α＝fs。

基于上述分析，本文提出自適應變分模態(tài)分解(AVMD)，避免了VMD 分解參數(shù)經(jīng)驗選擇導致嚴重重構誤差。 AVMD 算法提取心跳過程如下:首先根據(jù)采樣頻率設定二次懲罰因子α，隨后對~x(t)進行VMD 分解并計算每個IMF 分量與~x(t)的相關系數(shù)，設置相關系數(shù)閾值δ，一般可取0.1。 若相關系數(shù)最小值小于δ則停止分解，得到IMF 分量，其中最大的相關系數(shù)對應的IMF 分量即為心跳信號xh(t)。 算法1總結了AVMD 算法提取心跳信號的處理流程。

算法1 心跳信號提取(AVMD)

2.3 心跳頻譜估計

2.3.1 線性調頻Z 變換

線性調頻Z 變換(CZT)由Rabiner 等人[18]提出，相比傳統(tǒng)的FFT 只能得到粗略的“全體頻譜”，CZT 能實現(xiàn)信號的窄帶分析，顯著提高頻率估計準確性。 CZT 本質是計算沿Z 平面上一段螺旋線周線等角度間隔采樣的有限采樣點的Z 變換值，其定義可表示為:

式中:x(n)為采樣序列，n＝1，2，…，N-1；A0和θ0分別表示起始采樣點的矢量半徑和相角；W0表示螺旋線伸展趨勢；φ0表示相鄰采樣點間的夾角；M表示目標頻帶的采樣點數(shù)，m＝0，1，…，M-1。

算法2 詳細闡述CZT 算法心跳頻譜估計流程。

算法2 心跳頻譜估計(CZT)

2.3.2 心率估計性能對比

為分析CZT 算法心率估計的準確性，與傳統(tǒng)的FFT 進行對比分析。 以上述雷達信號處理過程得到的心跳信號xh(t)為例，ECG 測量的心率作為參考，信號采樣頻率fs＝32 Hz，時間窗口T＝3 s，靜坐心率頻帶為1 Hz～2 Hz，其頻譜計算結果如圖2 所示，實線和虛線分別表示CZT 和FFT 算法的局部頻譜。由圖可知，F(xiàn)FT 和CZT 估計的心率分別為1.333 Hz和1.427 Hz，參考心率為1.418 Hz。 對比心率估計誤差，F(xiàn)FT 是CZT 的10 倍左右。 此外，表1 列舉了四組不同時間段下兩種算法的心率估計誤差，CZT心率估計誤差明顯低于FFT，驗證了CZT 頻譜細化算法心率估計的準確性。

圖2 心跳信號局部頻譜

表1 不同時間段下心率估計誤差 單位:Hz

3 實驗過程及結果分析

為驗證本文提出的心跳檢測算法的有效性與魯棒性，論文在不同的實驗條件下測量了多組實驗數(shù)據(jù)，并采用多個評價指標對本文方法的性能進行評估。

3.1 實驗參數(shù)介紹

本文模擬駕駛行為的體征數(shù)據(jù)采集實驗如圖3所示。 實驗的參考心跳數(shù)據(jù)由傳統(tǒng)的接觸式ECG模塊測量(圖3 中右下角電路開發(fā)板)，ECG 電極貼于受試者心臟附近，采樣頻率設為500 Hz。 實驗采用德州儀器公司的IWR1443Boost 毫米波雷達平臺(圖3 中右上角電路開發(fā)板)來測量體征信息，其相關實驗參數(shù)如表2 所示。

圖3 模擬駕駛行為的體征數(shù)據(jù)采集實驗

表2 實驗參數(shù)

實驗開始時，受試者位于雷達正前方0.5 m 左右距離，初始時刻調節(jié)雷達支架的位置來確保受試者胸腔的位置與雷達基本處于同一水平面，此后保持正常呼吸，記錄實驗數(shù)據(jù)。 本文通過兩種不同狀態(tài)分別模擬駕駛的不同情形:靜坐狀態(tài)用于模擬駕駛過程中直線行駛等常規(guī)駕駛情形；往復擺臂運動用于模擬駕駛過程中轉向、人體扭動等復雜情形。

時間窗口的大小會影響檢測的精度，一般而言，時間窗口越短，時間分辨率越好，呈現(xiàn)更多心跳細節(jié)，但同時頻率分辨率越低，檢測精度越差。 為驗證本文算法的普適性，選擇3 s 和10 s 兩種不同時窗長度用于實時心率檢測。

3.2 實驗結果分析

通過分析3s 和10s 兩個不同時窗下心率和心率變異性相關指標，對比了FTPR、CZT 算法的心率檢測性能。

3.2.1 HR 分析

心率準確度(HRA)通常定義為多普勒雷達檢測到的心率與ECG 檢測到的參考心率偏差在一定偏差以內的時間百分比[19]，本文選取±2%、±3%這兩個偏差作為心率準確度指標。

為直觀地展示本文提出的算法在不同時窗下的性能優(yōu)勢，圖4 展示了靜坐時不同時窗下不同算法的心率估計曲線，以ECG 作為參考。 如圖4(a)所示，取心率偏差在±2%誤差范圍，F(xiàn)TPR、CZT 算法的心率準確率高達95.5%、97.3%。 然而，當心率出現(xiàn)較大突變時FTPR 心率曲線中存在嚴重的偏離點。如圖4(b)所示，取心率偏差在±2%誤差范圍，由于3s 時窗較短，頻譜分辨率較低，F(xiàn)TPR 算法的心率準確率僅為45.76%，但CZT 算法的心率準確率仍超過了91%，且呈現(xiàn)了更多的心跳細節(jié)。 此外，當短時心率發(fā)生較大突變時CZT 算法估計的心率值與參考心率一致性仍較好，進一步驗證了本文算法的有效性與魯棒性。

圖4 靜坐時不同時窗下心率估計曲線

圖5 展示了往復擺臂時不同時窗下不同算法的心率估計曲線。 往復擺臂時出現(xiàn)與胸壁運動相當?shù)纳眢w運動，會降低心率檢測精度。 如圖5(a)所示，取心率偏差在±2%或±3%誤差范圍，兩種算法心率準確率均超過95%。 如圖5(b)所示，取心率偏差在±2%誤差范圍時，F(xiàn)TPR、CZT 算法的心率準確率分別為64.41%、86.44%，兩者準確率均有所降低。 此外，取心率偏差在±3%誤差范圍時，兩者算法的心率準確率分別為88.98%、96.61%，這表明CZT 算法在模擬復雜駕駛的往復擺臂情形下仍能保證較高的心率檢測準確度。

圖5 往復擺臂時不同時窗下心率估計曲線

選取均方根誤差(RMSE)作為心率誤差評估指標:

式中:NBPM表示實驗測量時間內時間窗口的個數(shù)，BPMest(i)表示雷達估計的第i 個心率，BPMref(i)表示ECG 估計的第i個參考心率。

結合表3 中靜坐時心率估計結果，可進一步評估不同算法的實時心率估計性能。 由表3(a)可知，F(xiàn)TPR、CZT 算法的平均心率準確率(2%)分別為95.95%、99.25%，兩者的平均RMSE 均小于1BPM，表明估計心率與參考心率一致性較好。 如表3(b)所示，由于時間窗口變短，頻率分辨率降低，導致心率檢測精度降低，F(xiàn)TPR、CZT 算法的平均心率準確率(2%)分別為42.09%、89.27%，F(xiàn)TPR 算法檢測精度顯著降低，但CZT 算法檢測精度仍保持90%左右，算法魯棒性更好。 此外，CZT 算法的平均RMSE小于1BPM，誤差較小，短時窗下呈現(xiàn)出更多的心跳細節(jié)。

表3 靜坐時HR 分析結果

表4 為往復擺臂時心率估計結果。 由表4(a)可知，取±3%誤差范圍，F(xiàn)TPR、CZT 算法的平均心率準確率均超過98%。 由表4(b)可知，取±2%誤差范圍時，F(xiàn)TPR、CZT 算法的平均心率準確率分別為62.57%、78.67%。 取±3%誤差范圍時，F(xiàn)TPR、CZT算法的平均心率準確率分別為82.77%、92.23%，CZT 算法顯著優(yōu)于FTPR 算法。 另外，CZT 算法的平均RMSE 為1.52BPM，相對于靜坐時精度降低。

表4 擺臂時HR 分析結果

3.2.2 HRV 分析

心率變異性(HRV)是指連續(xù)心跳之間的時間間隔，揭示了每次心跳的規(guī)律性。 心率變異性分析對心臟疾病的診斷和治療具有重要意義。 不同于在時域中通過檢測相鄰心跳的時間間隔來確定心率，本文采用在頻域中搜索心跳頻率的心率檢測方法[20]，其中連續(xù)心跳時間間隔BBIs 可用BPM 來表示:

另外，平均相對誤差(MRE)能被表示為:

HRV 分析評價指標:RR 間期標準差(SDNN)、連續(xù)RR 間期的均方根差(RMSSD)，其中，NBBI為連續(xù)心跳間隔數(shù)，可分別表示為:

靜坐時心率變異性分析結果如表5 所示。

表5 靜坐時HRV 分析結果

如表5(a)所示，CZT 算法的MRE 誤差范圍為0.41%～0.60%，表明雷達與ECG 間一致性很好。 另外，所有的RMSSD 都小于3 ms，這表明長時窗時心率細節(jié)嚴重丟失，與心率分析結果一致。 如表5(b)所示，F(xiàn)TPR、CZT 算法的平均MRE 分別為2.34%、1.01%。 兩種算法與ECG 間心率變異性評價指標差值的平均值分別為SDNN(4.60 ms、1.68 ms)、RMSSD(4.93 ms、1.53 ms)。 對比上述誤差，CZT 檢測性能明顯優(yōu)于FTPR，驗證了短時窗下本文算法的優(yōu)越性。

往復擺臂時心率變異性分析結果如表6 所示。由表6(a) 可知，CZT 算法的MRE 誤差范圍為0.47%～0.85%，一致性仍較好。 由表6(b)可知，F(xiàn)TPR、CZT 算法的平均MRE 分別為1.87%、1.35%。往復擺臂運動時誤差同樣較小，進一步驗證了本文心跳檢測算法的魯棒性和有效性。

表6 擺臂時HRV 分析結果

4 結束語

本文通過理論分析和人體實驗，驗證了我們提出的基于毫米波雷達的實時心跳檢測算法的有效性和魯棒性。 在呼吸的強干擾下本文提出的心跳提取算法實現(xiàn)了心跳信號的精確恢復。 此外，將頻譜細化算法引入心跳頻譜估計，該方法實現(xiàn)了快速、準確的心率估計，并解決了短時窗下頻譜分辨率不足的難題。 模擬駕駛行為的實驗結果表明，該方法具有較高的準確率和較小的平均誤差，驗證了雷達與ECG 間心率特征的高度一致性。 因此，本文提出的實時心跳檢測方法對于非接觸式生命體征檢測領域存在一定研究價值。