廖洪海，何為，林水洋,3

(1 中國科學院上海微系統(tǒng)與信息技術研究所，上海 200050；2 中國科學院大學，北京 100049；3 隔空智能科技有限公司，上海 200050) (2019年12月20日收稿； 2020年4月23日收修改稿)

隨著物質生活水平和對自身健康狀況關注度的提高，人們對生命體征監(jiān)測技術提出了更高的需求，而非接觸生命體征監(jiān)測技術因具有體積小、功耗低、使用方便、實時性高等優(yōu)點得到眾多學者的廣泛關注[1]?，F(xiàn)階段，超聲波、WIFI和雷達等都被廣泛應用于非接觸生命體征監(jiān)測領域內，此類技術手段都是基于多普勒效應來實現(xiàn)提取心率或呼吸信號的功能[2-5]。超聲波設備存在功率高、噪聲大的缺點，WIFI設備存在干擾大和信號處理困難的缺點，由于雷達設備具備功耗低、噪聲小、信號處理技術成熟等優(yōu)點，非常適合應用到非接觸生命體征監(jiān)測技術中，所以受到眾多學者的青睞。雷達的生命體征監(jiān)測技術的研究分為3類，分別是基于單音連續(xù)波(continuous wave, CW)雷達，基于脈沖超寬帶(impulse radio ultra-wideband, IR-UWB)雷達和基于調頻連續(xù)波(frequency-modulated continuous wave, FMCW)雷達的生命體征監(jiān)測技術[6]。

CW雷達通過天線發(fā)射一個單音信號，通過剖析發(fā)射信號與反射信號的相位信息獲得距離。然而，CW雷達存在直流偏移和I / Q不平衡的問題，進而導致虛警等問題[7]。

IR-UWB雷達通過時域變換和信號處理技術分析收發(fā)機之間單極沖激脈沖信號的關系獲取目標的距離和方位等信息。再者，該技術距離識別精準度與頻率分辨率成正比，IR-UWB雷達的距離分辨率更大。但是由于脈沖信號占空比低、信噪比(signal-to-noise ratio，SNR)差，導致IR-UWB雷達存在距離分辨率低的缺點。

線性調頻連續(xù)波(linear frequency modulation continuous wave, LFMCW)雷達結合CW雷達基于相位精確測距的優(yōu)點和IR-UWB雷達距離分辨率與帶寬成正比的優(yōu)點，因此有兩大優(yōu)點[2]：

1)具有更大的帶寬以至于距離分辨率較高，因此該雷達在雜波隔離和微小運動檢測方面的能力更強。

2)基于相位實現(xiàn)精準測距，并且發(fā)射能量和信噪比較高的連續(xù)波，因此該雷達的距離分辨率更高。

近年來，在基于LFMCW雷達生命體征監(jiān)測技術的研究方面，相關學者也提出了眾多的理論模型和算法體系。劉旭陽[8]提出基于微波雷達的生命體征信號獲取與處理技術，使用集合經驗模態(tài)分解算法(ensemble empirical mode decomposition, EEMD)解決了靜止狀態(tài)下的單人呼吸心率檢測，但是因為EEMD方法不能完整解決經驗模態(tài)分解(empirical mode decomposition, EMD)中模態(tài)混淆的問題，所以該方法存在抗干擾能力不強、檢測范圍近的問題[9-10]。Wang等[11]提出基于5.8 GHz LFMCW雷達的生命體征監(jiān)測方法，由于該頻段的雷達帶寬僅有160 MHz，因此該方法僅能提供93.75 cm的距離分辨率，僅適用于呼吸檢測。Bakhtiari等[12]提出一種新的基于LM(Levenberg-Marquardt)擬合算法的心率估計算法，但是在有呼吸和身體運動的環(huán)境中該方法存在心率(heart rate, HR)信號難以檢測的問題。拜軍等[13]結合LM算法和I-Q正交技術實現(xiàn)了微弱呼吸和輕微體動環(huán)境中的心率檢測，但是該算法存在復雜度高和解決干擾不徹底的問題。

在上述問題的驅動下，本文在深入剖析載波頻率為77 GHz的LFMCW雷達技術原理的基礎上，提出基于排列熵(permutation entropy，PE)的改進集合經驗模態(tài)分解(modified ensemble empirical mode decomposition, MEEMD)篩選心率的估計算法，研究MEEMD去噪和心跳信號的PE區(qū)間選取技術，探索呼吸低頻信號和振幅大于心跳的相位雜波過濾手段。實驗表明，與快速傅里葉變換(fast Fourier transform,FFT)算法和LM擬合算法相比，本文提出的算法具有更低的算法復雜度和估計誤差，提高了估計穩(wěn)定性和準確率。因此，本文的研究可以更加準確穩(wěn)定地估計心跳，監(jiān)測人體的生命體征。

1 原理與建模

通常情況下，如果雷達被放在一個振蕩物體(例如近似靜止的身體中周期運動的心臟)前面，并發(fā)送一系列等間隔的chirp信號，該雷達的接收天線Tx接收到的每一個chirp信號和發(fā)射的chirp信號混頻之后將得到range-FFT頻譜中的峰值。根據(jù)LFMCW雷達測距原理中距離與頻率成正比的性質，那么在身體近似靜止的狀態(tài)下，上述峰值對應的頻率將不會發(fā)生太大變化，但是峰值對應的相位將隨著人體微弱的振蕩運動變化，即微動信號包含心臟運動的幅度和周期。但是由于考慮到在同一距離的身體其他部位微弱動作的疊加，實際上的微動信號包含諸多可被視為心率的噪聲，后文對這些噪聲進行了數(shù)學建模分析。MEEMD算法通過成對添加高斯白噪聲和均值策略，完整地保留了分解得到的模態(tài)分量的物理意義，解決了模態(tài)混疊的問題，具有理論完備、算法復雜度低、穩(wěn)定等優(yōu)點，所以可以通過MEEMD算法對微動信號進行分解。微動信號被MEEMD算法分解之后的結果中必定存在心臟的周期信號和其他的噪聲干擾信號，可以通過本文提出的PE閾值準確地篩選出心臟的周期信號，最后通過簡單的峰值檢測算法便能夠獲得心率HR。

1.1 基于PE的MEEMD篩選的LFMCW雷達心率估計算法

本文提出的算法是在對LFMCW雷達心率估計算法深入剖析的基礎上，研究基于PE的MEEMD濾波器去除呼吸低頻信號和振幅大于心跳的相位雜波等噪聲和干擾，并以精度更高的77 GHz毫米波雷達為研究載體，提出基于PE的MEEMD篩選的LFMCW雷達心率估計算法。通過設計更優(yōu)良的參數(shù)，達到更加精確的心率估計。算法流程圖如圖1所示，具體解釋如下：

圖1 基于PE的MEEMD心率估計算法的流程圖Fig.1 Flow chart of PE-based MEEMD heartrate evaluation algorithm

1)連續(xù)的N個中頻信號sIF(t)構成N×M的原始數(shù)據(jù)矩陣M[nslow,mfast]，其中nslow=1,2,…,N，mfast=1,2,…,M，N表示發(fā)送的chirp數(shù)，M是每個chirp的采樣點數(shù)。

2)對原始數(shù)據(jù)矩陣M[nslow,mfast]的快時間維度進行傅里葉變換得到Rp[nslow,mfast][14]。

5)通過峰值檢測算法獲得心率估值HR[15]。

1.2 雷達信號建模

雷達發(fā)送信號[16]為

sTx(t)=exp(j(2πfct+πγt2+φ)).

(1)

(2)

其中：σ表示接收信號的幅度，由反射物體的雷達散射截面(radar cross section, RCS)和傳播損耗共同決定。sRx(t)經過下變頻得到中頻信號sIF(t)

(3)

(4)

(5)

對微動信號建模

(6)

1.3 基于PE的MEEMD濾波器

1.3.1 心跳信號的PE區(qū)間選取

PE是一種一維時間序列的復雜度的平均熵度量方式[17]。較小的噪聲本質上不會改變含噪信號的復雜度，所以可以對任意現(xiàn)實世界時間序列計算排列熵。由于PE具有快速和健壯的特點，因此在存在大量數(shù)據(jù)集且沒有時間進行預處理和參數(shù)微調時是可取的[18]。PE的計算方法如下。對序列{x(i),i=1,2,…,N}進行相空間重構得到序列X

(7)

其中:m是嵌入維數(shù)，λ是時間延遲,X(k)={x(k),x(k+λ),…,x(k+(m-1)λ)},k=1,2,…N-(m-1)λ。將X(k)按升序排列，如下所示

Xsort(k)={x(i1),…，x(iq)，…，x(im)}.

(8)

其中：iq∈{k,k+λ,…,k+(m-1)λ}為下標，x(i1)<…

(9)

歸一化后得到

(10)

其中：Pj為全排列矩陣中的第j行出現(xiàn)的頻率，lnm!為排列熵為Hp(m)的最大值。Hp_norm(m)越大說明序列越隨機。

在采樣率確定的情況下，范圍為50～120 beat/min的心跳信號的歸一化排列熵Hp_norm(m)與心率呈正比，這與PE的定義相符。與此同時，采樣率降低，正比關系整體趨勢不變，雖然在局部會有波動，但是排列熵的區(qū)間用于篩選信號仍是有效的。如圖2(a)所示，當采樣率為20 Sa/s(sample per second, 每秒采樣次數(shù))，PE值線性分布在0.31～0.44。如圖2(b)所示，當采樣率升高到4 000 Sa/s，周期現(xiàn)象變強，PE值降低，線性分布在0.107～0.111?？梢愿鶕?jù)這個性質篩選MEEMD分解之后的分量，獲得心跳信號。

圖2 不同采樣率下50～120心跳的PE區(qū)間Fig.2 PE interval of 50-120 beat/min heartbeat at different sampling rates

1.3.2 基于PE區(qū)間的MEEMD心跳信號篩選方法

為消除或減弱靜態(tài)雜波、線性趨勢等干擾和噪聲，本文基于PE理論提出心跳信號的MEEMD篩選算法，具體實現(xiàn)步驟如圖3所示。

圖3 基于PE區(qū)間的心跳信號濾波算法流程圖Fig.3 Flow chart of heartbeat signal filteralgorithm based on PE interval

1) 分別在待處理信號S(t)中添加零均值白噪聲ni(t)和-ni(t)，如下

(11)

其中：ai為噪聲幅值，一般選擇信號的0.1～0.2的SD，i=1,2,…,Ne，Ne表示添加白噪聲對數(shù)，一般Ne≤100[13]。

(12)

4) 如果Im(t)是噪聲或者干擾，在原始信號S(t)中去掉Im(t)得到新的S(t)=S(t)-Im(t)，執(zhí)行步驟1)～4)，否則執(zhí)行步驟5)。

2 實驗與對比

2.1 實驗平臺

本實驗使用DCA1000EVM數(shù)據(jù)采集模塊和IWR1642毫米波雷達，使用mmWave Studio軟件設計雷達信號，使用MATLAB 2018b進行算法仿真。

2.2 LFMCW信號設計

本文設計雷達的chirp周期Tc為59 μs，瞬時頻率斜率α為67.495 MHz/μs，雷達的工作頻率從77 GHz線性變化到81 GHz，帶寬達到3.982 21 GHz。根據(jù)式(13)設計距離分辨率ΔR為3.75 cm，更小的距離分辨率能更加有效地隔離雜波。LFMCW雷達的具體參數(shù)如表1所示。

ΔR=c/2B.

(13)

其中：ΔR代表距離分辨率，c代表光速，B代表雷達帶寬。

由表1可知，本文設計的雷達的起始頻率fstart為77 GHz，根據(jù)式(14)可以求得微動信號的測量精度ΔRmic為0.302 2Δφmicmm,即該雷達的測量精度達到了毫米級。因為心跳引起的胸腔起伏的幅度為1～2 mm，所以本文設計的雷達從理論上做到了測量心跳所需的精度。與Wang等提出的基于5.8 GHz LFMCW雷達的呼吸測量算法[16]相比，本文設計的雷達有效地解決了微動信號精準度不足的問題，可以實現(xiàn)心率高精度監(jiān)測。

表1 77 GHz線性調頻連續(xù)波雷達參數(shù)Table 1 Parameters of 77 GHz LFMCW radar

ΔRmic=Δφmic×c/4πfc.

(14)

其中：Δφmic為相位，fc代表雷達的中心頻率，為78.991 105 GHz，由下式確定：

fc=fstart+B/2=78.991 105 GHz.

(15)

由上分析可知，本文還保證了每一個chirp的相干性，如圖4(a)所示,即在chirp周期Tc間預留一段Tu，并且保證每一個chirp起始相位相同，本文的初始相位為0°。發(fā)射波形如圖4(b)所示。接收信號是包含了環(huán)境信息的回波信號，如圖4(c)所示。

圖4 LFMCW雷達頻率曲線和收發(fā)波形Fig.4 LFMCW radar frequency curve and transceiver waveform

2.3 目標檢測

為更好地實現(xiàn)人體目標的心率信號監(jiān)測，首先需要使用一種有效的人體目標檢測算法確定待測人體目標的具體位置，本文采用結合以式(4)為理論基礎的距離譜圖和基于2D-FFT的速度-距離譜圖這2項綜合指標確定人體位置。最終以式(5)為理論基礎從距離譜圖中人地位置所對應的下標獲取含有心率信號的微動信號，為2.4節(jié)的心率估計做準備。本節(jié)布置了站姿、坐姿和平躺3種姿態(tài)的人體目標檢測環(huán)境，驗證目標檢測算法的有效性。

如圖5(a)所示，當環(huán)境中沒有待測目標時，中頻信號包含Rx-Tx耦合信號，距離較近的桌子邊緣和較遠靜態(tài)物體的反射波等多種成分組成，波形較復雜。圖5(b)能夠找到各個分量的值。同理，圖6(a)的距離譜圖也能夠找到0.133 3 m處的桌子邊緣等分量。

圖5 坐姿情況下的中頻信號和距離譜圖Fig.5 IF signal and distance spectrum in sitting style

如圖6(g)所示，當環(huán)境中放置1把椅子，圖5(d)中0.133 3 m處的桌子邊緣和較遠的靜態(tài)物體的能量相對較小，0.888 4 m處的椅子靠背和0.666 3 m處的椅子把手的能量較大，圖6(c)的距離譜圖也能夠找到椅子等分量。如圖5(c)所示，時域波形是椅子靠背的反射波疊加椅子把手等反射波，仍然比較復雜。

圖6 站立、平躺和坐姿3種姿態(tài)下的距離譜圖和速度距離譜圖Fig.6 Distance spectrum and speed-distance spectrum in three styles: standing, lying, and sitting

如圖6(h)所示，估計志愿者心率時，如圖5(f)所示，由于被0.577 5 m處的身體遮擋，椅子的反射能量急劇下降，Rx-Tx耦合，桌子邊緣和較遠的靜態(tài)物體的能量更小。如圖5(e)所示，由于距離分辨率小，能量集中，時域信號近視正弦/余弦信號。

本文使用2D-FFT計算速度-距離譜圖準確定位志愿者，有助于獲取微動信號。如圖6(f)所示，0.577 5 m處的志愿者存在0～8 m/s的不同速度，這表明該位置的物體存在瞬時運動，如圖6(d)所示，0.888 4 m處的椅子只有小于0.1 m/s的速度，這有可能是地板晃動導致的。根據(jù)人體瞬時運動的特性，使用速度-距離譜圖能識別人體，有助于提取微動信號，進行準確的心率估計[19]。

由于待測目標的整體平移運動，包含心跳的微動數(shù)據(jù)中可能引入間隙脈沖干擾，可以通過全局距離校準算法校正，避免影響微動數(shù)據(jù)的采集[20]。

為了驗證本文提出的方法適用于不同姿態(tài)的心率估計，本文還對站立姿態(tài)和平躺姿態(tài)進行實驗。如圖6(i)所示，能夠在距離雷達0.86 m處檢測到站立姿態(tài)的志愿者與其他位置的雜物。如圖6(j)所示，能夠在距離雷達0.64 m處檢測到平躺姿態(tài)的志愿者與其他位置的雜物。通過圖6(k)和6(l)所示，能夠輕松的檢測出待測目標的微動信號，這與圖6(m)和6(n)的實驗環(huán)境一致，充分證明了本文使用的檢測方法的有效性。

2.4 基于PE的MEEMD篩選算法的心跳估計

根據(jù)前文的分析可知，微動信號中包含身體的晃動，手臂、衣物、皮膚的晃動，呼吸引起的胸腔起伏等噪聲。如圖7(a)所示，身體的晃動幅度比心跳的幅度更大，心跳信號疊加在約為±0.1 m身體晃動信號之上。環(huán)境噪聲和身體的各種微動也疊加在微動信號中，因此從原始微動信號無法直接估計心跳。

圖7 FFT法，LM法和本文提出的方法的實驗結果對比Fig.7 Comparison of experimental results between FFT method, LM method, and the method proposed in this paper

由圖7(b)可知，現(xiàn)有的FFT方法通過在幅度譜中尋找心跳頻率，當0.8～2 Hz通帶內存在干擾時，這種方法因噪聲問題導致準確率降低。為更加形象地闡述上述問題，圖7(c)為放大后的FFT法頻率在0～2 Hz的信號特征，在該實驗中志愿者的真實心跳為88，然而FFT法的估計結果卻為79，準確率達到89.77%。而且這種方法的性能需要較高的采樣率和更大數(shù)據(jù)樣本空間支撐。例如，當采樣率為20 Sa/s時，該算法需要采樣12 000個樣本點，即10 min數(shù)據(jù)才能實現(xiàn)心跳分辨率為1跳。

圖7(d)為LM擬合法降噪之后的頻譜圖，圖7(e)為0.8～2 Hz的頻段放大，心跳為88的志愿者的估計結果也為79，準確率達到89.77%。該方法雖然能夠在一定程度上解決噪聲問題，但準確率仍然不高，微動干擾下性能提升不大。

如圖7(f)所示，本文提出的基于PE的MEEMD心跳篩選算法，能夠清晰識別微動信號中含有87個周期，由此可見，本文所提出的算法準確率達到98.86%。如圖7(g)和圖7(h)所示，為同一志愿者分別在平躺和站立姿態(tài)使用本文提出的算法得到的微動信號。由于距離雷達更近，信號幅度更大；平躺姿態(tài)下心率更緩。這2個特點都與實際相符。除此之外，采樣率為20 Sa/s時，僅需采樣1 200點便可估計心跳，與FFT方法相比具有更小的計算復雜度，提高了心跳測量的實時性和準確性。

2.5 同類算法比較

為比較本文提出的基于PE的MEEMD心跳篩選算法的穩(wěn)定性和準確性，定義準確率和穩(wěn)定性指標，實現(xiàn)上述算法的性能比較。

估計準確率A為

(16)

其中:HRestimate為估計心跳值，HRreal為心跳實際值，本文中以接觸式心率儀作為基準。

平均估計準確率Aaver為

(17)

其中：Ai為第i次估計的準確率，其中i=1,2,…,N，N為估計總次數(shù)。

算法準確率和穩(wěn)定性的綜合評價指標fHR為

fHR=α×mseHR+(1-α)stdHR.

(18)

其中：mseHR為均方誤差，如式(19)所示，代表估計心跳與基準心跳的累計誤差，表示估計的準確性，越小說明算法越準確。stdHR為心跳估計的總體偏差，如式(20)所示，表示估計的穩(wěn)定性，越小說明算法越穩(wěn)定。

(19)

(20)

實驗采集了30名志愿者共30組數(shù)據(jù)。表2記錄了采用3種算法與接觸式心率儀的心跳估計結果。為了進一步驗證本文提出的算法更加準確和穩(wěn)定，采用前文提出的2個評價指標Aaver和fHR進行算法的評估。圖8(a)展示了30名志愿者分別使用接觸式心率儀，F(xiàn)FT算法、LM擬合法和本文提出的基于PE的MEEMD心跳篩選方法的估計結果，本文提出的方法與接觸式心率儀的結果曲線更逼近[10,21]。圖8(b)展示了30名志愿者使用3種算法進行的心跳估計的準確率，證明本文提出的方法準確率更高。圖8(d)為30名志愿者使用3種算法進行心跳估計的fHR指標，證明本文提出的方法的fHR指標最小，即估計誤差最小和最穩(wěn)定，具體指標值見表3。圖8(c)展示了3種算法的平均準確率，其中FFT法為94.18%，LM擬合法為96.40%，本文方法為98.30%。綜上所述，本文提出的基于PE的MEEMD篩選方法的穩(wěn)定性和準確率更優(yōu)。

表2 3種算法估計的心跳和心率儀測量的心跳的數(shù)據(jù)Table 2 The data estimated by the three algorithms and heart rate measured by cardiotachometer

圖8 3種心跳估計算法的性能對比Fig.8 Performance comparison of three heartbeat estimation algorithms

表3 3種算法的綜合評價指標對比Table 3 Comparison of comprehensive evaluationindicators of the three algorithms

3 結論

本文提出基于PE的MEEMD篩選的LFMCW雷達心率估計算法，通過設計77 GHz的LFMCW雷達，并在對身體微動信號屬性進行充分剖析，對其提取手段深入研究的基礎上，提出基于PE的MEEMD濾波器去除呼吸低頻雜波、靜態(tài)雜波、線性趨勢等干擾和噪聲，提高了心跳估計的穩(wěn)定性和準確性。最后為了驗證和評估本文提出的算法的合理性和科學性，通過構建物理實驗平臺分別對現(xiàn)有的FFT方法、LM估計方法和本文提出的方法進行評估，實驗數(shù)據(jù)表明本文提出的算法降低了算法復雜度和估計誤差，提高了估計穩(wěn)定性和準確率，可以更加準確穩(wěn)定地估計心跳以實現(xiàn)人體生命體征的實時監(jiān)測。但是實時性的生命體征監(jiān)測技術的廣泛應用仍然有待提高，未來的研究工作會圍繞如何更加準確且快速地估計心率進行研究。