二次諧波加權重構的77 GHz FMCW雷達心率監測方法

鄭春弟，李 剛，陳薈慧，王愛國

(1.佛山科學技術學院 電子信息工程學院，廣東 佛山 528225；2.清華大學 電子工程系，北京 100084)

心率是衡量人體健康狀況的關鍵生理參數，其在臨床診斷、治療效果與健康狀態評估等方面具有指標意義[1-2]。傳統上，常用接觸式傳感器進行心率監測。但在實際使用中，接觸式監測技術需要佩戴或系縛相關設備，極易引起人體的不適，特別是導聯電極貼片式監測技術，長時間監測會妨礙人體其他活動，進而造成生理與心理上的抵觸[3]，不利于24小時動態心電圖監測的順利執行。因此，隨著人們對美好生活需求的日益增長，傳統接觸式心率監測技術已無法滿足人們對智能化、無感化、舒適化健康監測設備的強烈需求。

與需要系縛的接觸式傳感器相比，心率監測雷達以電磁波為媒介，通過測量由心臟跳動引發的胸壁微小位移，能以非侵入、非接觸的方式獲取心臟跳動的頻率，可提供舒適感極高、用戶友好的心率監測方案，因此雷達技術有望成為最有應用潛力的心率監測技術[3]。雖然雷達測量心率具有非侵入、非接觸、全天時的天然優勢，但心率監測雷達發展也面臨著一系列挑戰性問題。特別地，胸壁位移由呼吸和心跳兩種疊加機械運動引起，這兩種運動的雷達回波信號在時域和頻域上相互疊加，而且心跳的位移幅度(0.01～0.5 mm)遠低于呼吸的位移幅度(約1～12 mm)[4]。因此，如何在呼吸(頻率約0.1～0.6 Hz)及其諧波強干擾條件下，測量心跳的頻率(頻率約0.8～2.5 Hz)是心率監測雷達面臨的最大挑戰之一[5-6]。

基于77 GHz調頻連續波(FMCW)雷達的心率監測技術是近幾年新涌現的生命體征監測技術[4，7-9]。總體而言，現有的77 GHz毫米波雷達心率監測技術，在處理呼吸及其諧波強干擾下的心率估計問題時，采用的方法主要有三類，即帶通濾波法[8]、最優頻率估計法[4]和多分辨率分析(MRA)法[7]。盡管通過設置合適的低端截止頻率，這些方法能夠將呼吸的基頻及二次諧波或者僅基頻濾除掉(具體情況取決于呼吸頻率的大小)。但是這些方法在濾波之后得到的信號并不是干凈的心跳信號，其中依然可能保留著呼吸的三次或者二、三次諧波，這些殘存的呼吸高次諧波將嚴重制約著心率估計的準確性。

針對呼吸及其諧波強干擾下的心率估計問題，提出了一種基于心跳二次諧波信號加權重構的心率估計方法，在77 GHz FMCW雷達上實現了心率的精準監測。由于心跳二次諧波所在頻帶遠高于呼吸及其二、三次諧波的頻帶，因此在信號重構時僅使用心跳二次諧波所在頻帶的信號，可容易地剔除呼吸及其二、三次諧波的干擾[10]。為了獲取較為干凈的心跳二次諧波信號，提出方法使用多分辨率分析對胸壁運動信號進行多頻帶分解。提出方法與文獻[7]中基于多分辨率分析的77 GHz雷達心率估計方法相比，最主要區別有以下兩點：① 在多分辨率分析基礎上，提出方法進一步將能量信息運用到心跳二次諧波信號的重構之中，而文獻[7]中的方法并未使用各頻帶能量所蘊含的信息。具體而言，提出方法依據能量大小對各頻帶進行加權重構，這非常有助于凸顯心跳二次諧波最可能出現的頻率范圍。② 提出方法是基于心跳二次諧波進行心率估計的，而文獻[7]是基于心跳基頻進行心率估計的。相比于心跳二次諧波的強度，處在心跳二次諧波范圍的呼吸諧波、呼吸心跳交調諧波的強度幾乎可以忽略[10]，因此基于心跳二次諧波進行心率估計，不僅可以避免呼吸及其諧波的強干擾，而且還可避免后續譜峰估計中模型階數選擇難題。

1 系統模型

簡單的毫米波FMCW心率監測雷達發射與接收主要模塊如圖1所示。

圖1 毫米波FMCW心率監測雷達發射與接收示意圖

呼吸和心跳會引發胸壁機械運動，該運動將對照射的雷達波產生調制作用。首先通過對回波進行相應的解調，得到拍頻信號，然后通過適當的信號處理，雷達能夠恢復出胸壁機械運動的相關信息，進而提取出心率。

假設胸壁與雷達接收天線之間的距離為d0、發射線性調頻(Chirp)信號的調頻率為K、最大波長為λmax，則解調之后接收機正交和同相兩個通道上的拍頻信號可以表示為[4，11]

(1)

(2)

其中，Ab為拍頻信號的幅度，fb為拍頻，4πd0/λmax代表著胸壁徑向距離產生的常數相位，4πx(t)/λmax代表著由于胸壁運動引發的相位變化情況。x(t)=xh(t)+xr(t)，為心跳xh(t)和呼吸xr(t)疊加引發的胸壁準周期機械運動信號。通常，正交與同相兩路信號可以用復信號形式表示為

C(t)=I(t)+jQ(t)。

(3)

由于復信號C(t)的相位中蘊含著胸壁運動信息，因此在提取相位信息的基礎上，通過慢時間維一系列合適的信號處理方法，可估計出心跳的頻率。

由式(1)和式(2)可知，相位4πx(t)/λmax的大小與波長有關，即波長越短(工作頻率越高)，則胸壁運動引發的相位調制量越大[4，11]。例如，與24 GHz雷達相比，在同樣的胸壁運動幅度下，77 GHz雷達引發的相位調制量大約是前者的3.21倍，即77 GHz雷達的相位功率比24 GHz雷達大約高10 dB。與其他低頻段的心率監測雷達相比，除了有利于信噪比的提升之外，使用77 GHz頻段雷達還具有以下兩個優勢：① 77 GHz雷達已成為工業和汽車雷達的主流產品，大規模的商業應用有助于以極低的硬件成本實現心率監測；② 該頻段可用的帶寬大，非常有利于設計功能強大的心率監測雷達系統。鑒于以上優勢，文中選用77 GHz頻段雷達作為心率監測傳感器。

2 利用心跳二次諧波信號加權重構的心率估計方法

提出的心率估計方法流程如圖2所示。該方法主要包含3步：第1步是在估計胸壁位置的基礎上，對信號進行預處理以濾除雜波與噪聲等，從而得到較為干凈的胸壁運動信號；第2步是心跳二次諧波信號加權重構；第3步是基于多重信號分類(MUSIC)算法的頻率估計。下面分別對其進行說明。

圖2 基于心跳二次諧波信號加權重構的心率估計方法流程圖

2.1 預處理

2.1.1 胸壁位置估計

首先對時間窗內的復信號C(t)逐Chirp求快速傅里葉變換(FFT)得到距離像，然后利用文獻[12]中的標準差法計算每個距離單元的信號能量，最后將標準差最大行作為該時間窗內胸壁所在的距離單元。

2.1.2 預白化處理

在獲得胸壁所在距離單元的基礎上，對時間窗內信號的實部與虛部進行聯合預白化處理，以降低噪聲對后續相位估計的影響。

2.1.3 信號補償

在進行心率監測時，雷達元器件中的電路非理想等不利因素，會給I/Q通道引入幅度誤差、增益誤差和相位偏移等[11]。由于求取相位是一個高度的非線性過程，上述這些誤差極有可能導致較大的相位估計誤差，因此在進行信號相位提取之前，需要進行補償。文中采用文獻[13]中的方法進行I/Q通道不一致補償。

2.1.4 相位估計及解纏繞

在信號補償之后，使用四象限反正切法提取信號的相位信息。當慢時間維的采樣率設置合理時，相鄰兩次采樣獲取的相位應處在[-π，π]之內。但是當胸壁位移x(t)>λmax/4時，導致相位可能不落在[-π，π]之內，因此需要相位解纏繞處理。即當慢時間維相鄰兩次采樣的相位躍變超過π時，則需要通過加減2π的倍數來更正相位。之后再進行線性去勢操作，消除傳感器數據偏移對后續處理的影響，以便于將分析集中在相位本身的波動上。

2.2 心跳二次諧波信號的加權重構

將經過預處理獲得的相位信號記為φ(t)，首先使用Coiflet小波對其進行最大重疊離散小波變換，小波的消失矩階數選為5，層數為7，將得到的系數矩陣記為W∈(L+1)×M，其中W的第k行表示尺度為2k的小波系數，L代表層數，M代表慢時間的采樣數。然后再次使用消失矩階數為5的Coiflet小波對W進行多分辨率分析，得到矩陣WMRA。

根據多分辨率分析的含義可知：

(4)

其中行向量φ∈1×M表示時間窗內φ(t)的M個采樣組成的相位信號，表示WMRA的第k行。不妨假設心跳二次諧波所處的頻率區間[1.6 Hz，5.0 Hz]大致與WMRA的第m至第n行相對應，則此時心跳二次諧波所在頻帶的信號可以表示為

(5)

但是應當注意到，式(5)中并未體現各子頻帶的能量分布特征。

事實上，小波變換屬于線性變換，滿足能量守恒定律，W各層(對應不同的子頻帶)的小波系數具有能量量綱，它們的能量分布特征對于信號的頻率特性分析是非常有價值的。由于W各層能量的大小反映了原信號更有可能分布在哪些頻帶，因此可根據能量大小對WMRA進行加權處理，以凸顯心跳二次諧波信號更可能出現的頻帶。將W各層的能量記為[e1，…，eL+1]，則經過能量加權之后的心跳二次諧波信號可以被表示為

(6)

2.3 基于MUSIC算法的頻率估計

R=ΦΦH。

(7)

對其進行奇異值分解，可得

R=UDVH。

(8)

由于是在心跳二次諧波頻譜范圍內進行頻率估計的，此范圍內的呼吸諧波、呼吸心跳交調諧波的強度幾乎可以忽略[15]，且對各層以能量大小進行了加權處理，能量較小頻帶的貢獻進一步被削弱，因此可以將信號的模型階數定為1，即可認為僅有心跳二次諧波一個信號存在。所以此時的噪聲子空間可以表示為Un=U(∶，2∶K)，故而MUSIC譜峰可以表示為

(9)

3 實驗結果及分析

3.1 實驗設置

雷達傳感器設置：實驗使用德州儀器(TI)公司的IWR1642雷達傳感器采集胸壁運動信號。實驗中傳感器參數設置如下：起始頻率為77 GHz，調頻率為66.626 MHz/μs，Chirp脈沖寬度為60 μs，每個Chirp的ADC采樣數為256，快時間采樣率為5 MHz，每幀包含128個Chirp，幀的周期為40 ms，采用一發一收模式采集胸壁運動信號。

實驗場景設置：如圖3所示，在書房環境下采集數據，胸壁與雷達之間的距離約在0.55～0.95 m之間，門、墻等距離雷達大約1.7 m，書柜側壁距離雷達大約1 m。測量時，被測者靜坐于雷達前，并保持自然呼吸狀態。

圖3 實際測量場景

參考標準：以力康PC-3000多參數監護儀獲取的心率數據(ECG)為參考，其測量的心跳R-R間期被轉化為每分鐘跳動次數。心率監測時，監護儀和IWR1642雷達傳感器同時進行測量，其中監護儀需要在胸壁黏貼導聯電極貼片。

3.2 測量結果及其分析

在心率監測時，由于雷達參數設置為128個Chirp每幀，因此既可以每幀僅抽取1個Chirp進行心率估計，也可以每次每幀僅抽取1個Chirp，然后重復128次，再對128個測量結果取均值的方式估計心率，文中選擇后一種方式。在進行數據處理時，時間窗長為5.12 s，每隔0.72 s滑動1次。對于可能出現的野值，采用類似于文獻[16]中的方法消除。為驗證提出方法的有效性，實驗時文獻[7]中的方法僅在信號重構與譜估計階段與提出方法有所差異，其他環節兩者均采用相同的處理技術。

3.2.1 加權重構效果分析

為了驗證提出的心跳二次諧波信號加權重構方法的有效性，圖4給出了加權之前、之后心跳二次諧波信號重構的時域與頻域波形圖。如圖4(a)中2個圓圈內的時域波形所示，重構信號加權處理之前，時域波形存在著一些幅度較小的快速波動，而加權處理之后，這些快速波動信號將被平滑掉。如圖4(b)所示，這些快速波動代表著高頻噪聲或干擾成分，它們會導致心跳二次諧波信號頻譜向更高的頻帶上偏移，從而使得估計的心跳二次諧波頻率與參考值對應的2倍頻率之間出現較大的偏差。但依據能量大小對各子頻帶進行加權重構之后，高頻噪聲或干擾明顯被削弱，心跳二次諧波信號將被凸顯出來，進而有利于心率估計準確性的提升。

(a) 信號時域波形圖

(a) 第1人次測量結果

3.2.2 心率測量結果分析

為評價心率測量結果的有效性，平均值、平均值相對誤差(ARE)、平均絕對誤差(AAE)和平均絕對誤差百分比(AAEP)這4個常用統計指標被用來衡量測量值與參考值之間的誤差大小[9，17]。其中后三個指標定義如下：

(10)

(11)

(12)

其中，N、hest(n)和href(n)分別表示總的心率監測次數、第n次心率的測量值和參考值。

實驗數據來自6人次，測量對象為成年男女，其中第1人次為女性，其他人次均為男性。測量時胸壁與雷達之間的距離約為0.55～0.95 m。部分測量結果隨時間變化曲線如圖5所示。

如表1所示，通過文中方法獲取的平均心率與監測儀相比，兩者非常接近。表2中給出的平均心率相對誤差進一步說明了文中方法能夠準確地監測平均心率。相比之下，文獻[7]中的方法獲取的平均心率與參考值之間存在較大誤差。

表1 各人次測量的心率平均值 次/分鐘

表2 各人次測量的平均心率相對誤差(ARE) 次/分鐘

表3和表4中分別給出了各人次測量的平均絕對誤差和平均絕對誤差百分比，這兩個統計參數再次表明文中提出的心跳二次諧波加權重構方法，能夠獲得更加準確的心率測量結果。值得注意的是，第5、6人次測量值與參考值之間差異相對較大，主要原因在于此時雷達與胸壁之間的距離相對較大，從而導致信噪比降低，進而使得測量精度有所下降。

表3 各次測量結果的平均絕對誤差(AAE) 次/分鐘

表4 各次測量結果的平均絕對誤差百分比(AAEP) %

文獻[9]使用80 GHz FMCW雷達進行心率監測，當胸壁距離雷達1 m時，其平均心率的相對誤差在1.67%～39.24%之間，誤差的中值為8.09%。相比之下，文中提出方法在胸壁與雷達之間距離0.55～0.95 m時平均心率的相對誤差為0.01%～3.14%(見表2)。此外，文中方法與文獻[7]均使用了多分辨率分析，但不同之處在于文中方法考慮了各子頻帶能量分布特征對心跳二次諧波重構的影響，通過能量加權處理，凸顯了二次諧波更可能出現的頻帶，使得心率估計精度得以明顯提升。例如，實驗中時間窗為5.12 s時，最小平均絕對誤差僅為2.77次/分鐘，最小平均絕對誤差百分比僅為4.50%，明顯優于文獻[7]中的方法。因此，與文獻[9]和文獻[7]中相同頻段雷達心率監測方法相比，文中方法有著更高的測量精度。

4 結束語

根據多分辨率分析的結果，將各頻帶能量分布信息引入到心跳二次諧波的重構之中，利用77 GHz FMCW雷達實現了短時窗高精度的心率估計。在胸壁運動信號多分辨率分析之后，提出方法依據能量大小對心跳二次諧波范圍內的各子頻帶進行加權重構，這不僅有助于凸顯心跳二次諧波最可能出現的頻率范圍，而且也有效避免了呼吸及其諧波的強干擾和模型階數選擇難題。實測數據表明，與其他77 GHz雷達測量方法相比，筆者提出的方法有著較高的心率測量精度。在現有工作基礎上，未來將著重研究時間窗更短的瞬時心率估計方法，以便于進一步分析心率變異性。