劉璐瑤,張 森,肖文棟
1) 北京科技大學自動化學院,北京 100083 2) 北京市工業(yè)波譜成像工程技術研究中心,北京 100083 3) 北京科技大學順德研究生院,廣東 528399
生理信號蘊含著許多有價值的人體健康信息,可用于疾病的診斷和預防. 一般情況下,在心臟驟停等不良事件發(fā)生6~24 h前,心率呼吸等生命體征會出現(xiàn)異常[1]. 近年來,許多可穿戴傳感器被用來檢測生命體征信號,例如心電圖(ECG)、光容積描記(PPG)、呼吸帶等[2?5]. 雖然這些傳感器的測量結(jié)果相對準確,但往往給目標對象帶來不適及額外的負擔,特別是對一些特殊的人群如嬰兒和燒傷患者等[6]. 近年來,非接觸式生理信號檢測在睡眠呼吸暫停監(jiān)測、嬰兒猝死綜合征(SIDS)監(jiān)測、疲勞監(jiān)測、臨床醫(yī)療、家庭保健等方面受到了越來越多的關注[7?8]. 不同于傳統(tǒng)的測量方法,非接觸式生理信號檢測不需要人體佩戴任何傳感器,不會造成額外的負擔. 非接觸式生理信號檢測方法包括紅外、視頻成像、靜電場、超聲波、電磁波等等. 例如,成像式光電容積描記技術[9](Imaging photo-plethysmography, IPPG)及遠程光電描記技術[10](Remote photo-plethysmography, RPPG)使用攝像頭等電子成像設備采集人體體表皮膚視頻信息,經(jīng)處理提取人體生理參數(shù). 文獻[11]~[12]已經(jīng)嘗試解決光照環(huán)境變化及運動干擾對生理信號檢測的影響,但這對基于IPPG或者RPPG技術的生理信號檢測仍然是一個挑戰(zhàn). 紅外成像技術[13]也被用于非接觸式生理信號檢測,但易受溫度、天氣等環(huán)境因素的影響. 與上述技術相比,毫米波雷達不受光照、溫度等因素影響,穿透能力較強,不侵犯個人隱私,可實現(xiàn)毫米級高精度檢測.
近年來,非接觸式生理信號檢測使用的雷達主要有三種類型:連續(xù)波(Continuous wave, CW)多普勒雷達[14?15]、超寬帶(Ultra-Wideband, UWB)脈沖 雷 達[16?17]和 調(diào) 頻 連 續(xù) 波 (Frequency modulated continuous wave, FMCW)雷達[18?19]. 連續(xù)波多普勒雷達具有結(jié)構(gòu)簡單、功耗低的優(yōu)點,但沒有距離分辨率[20],因此其生理信號檢測容易受環(huán)境中其他物體或人體反射信號的干擾. 超寬帶雷達系統(tǒng)具有穿透力強、距離分辨率高等特點,但信號易受脈沖寬度和峰值信號強度的控制[21]. FMCW雷達不僅具有超寬帶雷達的測距能力,而且具有連續(xù)波多普勒雷達的靈敏度和魯棒性[22]. 此外,F(xiàn)MCW雷達具有體積小、重量輕、功耗低等優(yōu)點[23]. 在本文中,我們選擇工作在76.4 GHz頻段的毫米波FMCW雷達檢測生理信號.
目前,生理信號非接觸式檢測的信號處理方法主要有三種:基于快速傅里葉變換(Fast fourier transform, FFT)的 方法[24?25],基 于連續(xù)小 波變換(Continuous wave transform, CWT)的方法[26?27]以及基于時域信號處理的方法[16, 28]. FFT可以獲得生命體征速率,但不能跟蹤生命體征隨時間的變化[29].短時傅里葉變換(Short time fourier transform, STFT)已用于雷達心跳呼吸檢測,但其窗口長度不能隨頻率和時間變化,限制了算法的頻率分辨率[25]. 與FFT相比,CWT具有更靈活的時頻分辨率,可以提高生命體征的檢測精度. 文獻[26]使用具有高分辨率的時間頻率譜(Time frequency spectrum, TFS)獲得更準確的呼吸和心率. 除了傳統(tǒng)的頻域技術,如FFT和CWT,許多時域信號處理技術也得到了廣泛的發(fā)展. 文獻[28]提出峰值檢測方法計算生命體征速率,但對時域信號波形有較高的要求,除非獲得的呼吸心率信號有明顯的尖峰,否則檢測效果較差. 近些年,壓縮感知(Compressed sensing,CS)被用于計算生命信號隨時間變換的頻率[30],但如果生命體征速率在短時間內(nèi)發(fā)生顯著變化,其檢測精度也會大大降低.
以上介紹的方法雖然已經(jīng)被證明能夠?qū)崿F(xiàn)非接觸式生命體征檢測,由于心率與呼吸諧波存在頻率范圍重疊,可能會將呼吸諧波頻率誤判為心率. 同時心率信號本身比較微弱,容易被噪聲淹沒,導致心率檢測錯誤. 上述兩個問題是目前生理信號非接觸式檢測精度低的主要原因. 因此,本文結(jié)合小波分析和自相關計算減小呼吸諧波、環(huán)境等雜波對生理信號的影響,提高檢測精度.
圖1為基于FMCW雷達的非接觸式生理信號檢測模型,該雷達系統(tǒng)主要包括信號發(fā)生器、放大器、低通濾波器、ADC模塊等. 假定FMCW雷達與目標對象的距離為,表示目標對象生理運動引起的胸壁運動位移,表示目標對象胸壁運動與雷達天線之間的距離變化關系,其發(fā)射信號可近似表示為
圖1 基于 FMCW雷達的非接觸式生理信號檢測模型Fig.1 Noncontact vital signs detection model based on FMCW radar
發(fā)射信號從發(fā)射機天線(TX)向目標對象胸部發(fā)射. 經(jīng)過胸部反射之后接收機天線(RX)獲得的反射信號可近似表示為:
本文提出的基于FMCW雷達的非接觸式生理信號檢測方法總體流程圖如圖2所示,包括信號預處理和生理信號提取. 信號預處理由目標檢測、直流偏置去除及相位解纏構(gòu)成,旨在從采集的數(shù)據(jù)中準確提取出對應于目標對象的相位信息.生理信號提取由小波包分解、自相關計算及連續(xù)小波變換構(gòu)成. 小波包分解被用于分離重構(gòu)呼吸心跳信號,自相關計算被用于降低雜波對心跳信號的影響,連續(xù)小波變換被用于對呼吸心跳信號進行時頻分析,提取生理信號速率.
圖2 基于FMCW雷達的非接觸式生理信號檢測方法流程圖Fig.2 Noncontact vital signs detection processing procedure based on FMCW radar
2.1.1 目標檢測
雷達信號的發(fā)射是面向整個空間環(huán)境,其中除了目標物體還包含其他物體,因此反射回來的信號存在很多雜波信息. 為了準確地提取對應于目標對象的相位信息,在雷達視角范圍內(nèi)識別出目標對象的位置是十分必要的. 目標位置的識別通過FMCW雷達測距實現(xiàn),首先對采集的ADC數(shù)據(jù)進行快速傅立葉變換(FFT)得到和環(huán)境中各物體一一對應的距離信息,然后選擇與目標對象相對應的距離. 最后,沿著選定的距離提取對應于目標對象的相位信息.
2.1.2 直流偏置去除
對于FMCW雷達,在相位解調(diào)前必須消除復信號虛分量和實分量的直流偏置,否則會影響相
2.1.3 相位解纏
用下面arctan函數(shù)計算雷達信號的相位
2.2.1 小波包分解
經(jīng)信號預處理獲取的相位信號為目標對象的呼吸、心跳及其他雜波的混合. 為了提取呼吸及心跳信號,本文對相位信號進行小波包分解(WPD). WPD適用于非平穩(wěn)信號的時頻局部分析,與小波分解相比,小波包分析具有更高的時頻分辨率.
本文對得到的相位信號進行6級小波包分解,如圖3所示. 在第6層,可以得到包含64個節(jié)點的小波系數(shù),這些節(jié)點間的頻率差為0.15625 Hz.第1至第3個節(jié)點的低頻分量用于重構(gòu)呼吸信號,第6至第12個節(jié)點的高頻分量用于重構(gòu)心跳信號.
圖3 小波包分解圖Fig.3 Wavelet packet decomposition diagram
2.2.2 自相關計算
通過WPD獲取的心跳信號仍然存在一些干擾,可能會影響心率檢測的準確性,例如呼吸信號高次諧波的頻率可能會被誤判為心率. 為了解決這一問題,本文利用自相關計算從快時間軸和慢時間軸兩個方面對心跳信號進行分析. 自相關通過計算序列信號在不同時間與自身的相似性,提取被噪聲掩沒的周期信號,被看作是時間間隔的函數(shù). 快時間軸上的信號為原始心跳信號,慢時間軸上的信號為心跳信號經(jīng)過平移之后與原始心跳信號的自相關系數(shù). 自相關計算分析過程如圖4所示. 圖中 H R(t)表示原始心跳信號, H RA(t)表示經(jīng)過自相關計算后的信號,T表示原始心跳信號慢時間軸上開始的時刻,表示原始心跳信號沿著慢時間軸平移的時間. 因為心跳信號具有較好的周期性,所以心跳信號的自相關系數(shù)相比于雜波更明顯,可以有效突出心跳信號,提高心跳信號在頻譜分析中的分辨率. 具體過程如下.
圖4 自相關計算示意圖Fig.4 Autocorrelation computation diagram
(2)基于自相關計算公式獲取平移后的心跳信號與原始心跳信號的自相關系數(shù).
2.2.3 連續(xù)小波變換
連 續(xù) 小 波 變 換 (Continuous wavelet transform,CWT)是一種處理非平穩(wěn)信號的方法. 它一方面通過尺度變換提供信號的頻譜信息,另一方面通過小波變換保留信號的時域信息. 一個常見的小波變換計算方法如下:
心跳信號經(jīng)過自相關計算得到與其相對應的自相關系數(shù),利用連續(xù)小波變換將自相關系數(shù)變換為包含“時間?頻率”的頻域信息,通過時間頻譜(Time frequency spectrum, TFS)得到和時間相對應的頻率,最后利用滑動窗口算法得到隨時間變化的心率. 主要步驟如下:
(1)對自相關系數(shù)進行連續(xù)小波變換,通過改變縮放因子e及時移因子v覆蓋整個信號得到一系列小波系數(shù). 該小波系數(shù)可看作一個二維矩陣,矩陣的行表示時間信息,列表示頻率信息.
(2)以矩陣行表示的時間為單位,沿著行選擇每一個時間點對應頻率的最大值,得到TFS.
(3)利用滑動窗口算法進行頻率計算. 對于心率,窗口長度設置為 24 s,滑動窗口長度設置為 6 s;對于呼吸,窗口長度設置為60 s,滑動窗口長度設置為 6 s.
本文采用德州儀器(Texas Instruments, TI)毫米波IWR1642雷達傳感器,工作頻率范圍為77~81 GHz,有兩個發(fā)射天線和四個接收天線. 采集數(shù)據(jù)時,數(shù)據(jù)通過USB接口傳輸?shù)絇C機,后續(xù)信號處理在matlab上進行.
表1總結(jié)了FMCW雷達配置參數(shù). 每個frame包含128個chirps,每個chirp采集68個數(shù)據(jù)點. Frame發(fā)射間隔為 50 ms,即慢采樣率fslow為 20 Hz. 此外,ADC 采樣速率(快采樣率ffast)為 3.2 MHz. 雷達起始工作頻率為76.4 GHz. 每個chirp循環(huán)時間為64 μs,前一個chirp結(jié)束到下一個chirp開始之間的空閑時間16 μs,即一個chirp的持續(xù)時間Td=48 μs. 每個 chirp的斜率S為 20 MHz·μs?1,掃頻帶寬為 960 MHz.
表1 雷達參數(shù)配置Table 1 Radar configuration parameters
所有實驗均在普通室內(nèi)實驗室進行,如圖5所示. 受試者分別被要求坐在距離雷達0.5、1.0、1.5、2.0、2.5和 3.0 m 的椅子上,面向雷達. 利用可穿戴傳感器CM19(生產(chǎn)廠家:江蘇智海電子技術有限公司;型號:DEG-01-S;精度:采樣率 300 Hz)同時測量ECG信號和呼吸信號作為參考.
圖5 基于FMCW雷達的非接觸式生理信號檢測實驗場景Fig.5 Scenario of noncontact vital signs detection based on FMCW radar.
為了評價本文方法的實驗效果,本文引用平均絕對誤差(Average absolute error, AAE)和平均絕對誤差率(Average absolute error percentage, AAEP)作為評價指標[31],分別定義如下:
經(jīng)預處理獲取的相位信號為心率、呼吸及噪聲的混合. 圖 6分別展示了 50 s的相位信號(a)以及相位信號的頻譜(b). 從圖 6(a)可以發(fā)現(xiàn),呼吸信號具有較大的振幅且變化明顯,疊加在呼吸波形上的微小波動為心跳信號. 圖6(b)中表示呼吸頻率,表示心跳頻率,心跳頻率相比呼吸頻率及呼吸的2次諧波頻率是微弱的. 從相位圖及其頻譜中可以發(fā)現(xiàn)心跳信號相比于呼吸信號是微弱的,更容易受呼吸諧波及噪聲的干擾.
圖6 雷達相位信號(a)以及雷達相位信號頻譜(b)Fig.6 Radar phase signal (a) and radar phase frequency spectrum (b)
經(jīng)小波包分解和移動平均濾波獲取的心跳及呼吸信號如圖7所示. 120 s的雷達呼吸信號和參考信號比較如圖 7(a)所示. 圖 7(b)展示了 30 s的心跳信號和ECG參考信號的比較結(jié)果. 該方法獲得的心跳信號具有清晰可見的峰值并與參考ECG信號的R峰一一對應. 對于呼吸,結(jié)果也表明雷達的測量信號和參考信號之間具有較高的一致性.
圖7 雷達和參考傳感器的呼吸和心跳信號比較. (a)呼吸信號;(b)心跳信號Fig.7 Time domain respiration and heartbeat signals from the radar system and reference sensor: (a) respiration signal; (b) heartbeat signal
圖8(a)和(b)顯示了 180 s的呼吸速率(Breath rate, BR)和心跳速率(Heartbeat ratee, HR). 對于呼吸,每分鐘的次數(shù)(Beats per minute, bpm)為 17,本文方法提取的BR與參考傳感器的測量結(jié)果幾乎是完全重合的. 心跳的平均速率約為每分鐘75次(bpm). 雖然雷達測量值與參考結(jié)果之間存在較小偏差,但總體變化趨勢保持一致,所以本文方法獲得的心率結(jié)果也是值得信賴的.
圖8 雷達和參考信號的呼吸速率及心跳速率. (a)呼吸速率;(b)心跳速率Fig.8 Instantaneous BR and HR from the radar system and reference signal: (a) instantaneous BR; (b) instantaneous HR
為進一步評價本文方法的性能,選取10名不同受試者(6名男性,4名女性)進行實驗. 實驗對象沒有已知的心臟、呼吸系統(tǒng)或任何其他疾病.如前所述,受試者被要求坐在椅子上面對雷達,每個受試者在六個距離(0.5, 1.0, 1.5,2.0,2.5和 3.0 m)進行180 s的生理信號檢測實驗. 實驗過程中受試者進行正常呼吸,測試結(jié)果如表2和表3所示. 在距離雷達0.5 m,10名受試者呼吸率和心率檢測平均絕對誤差率的平均值分別小于1.65%和1.83%.受試者距離雷達小于2.5 m時,心率及呼吸檢測平均絕對誤差率的平均值分別小于5.01%和3.44%;平均絕對誤差的平均值分別小于4.39(bpm)和0.69 (bpm). 當距離大于 2.5 m,生理信號檢測平均絕對錯誤率及平均絕對誤差會明顯增加.
表2 10名受試者在不同距離生理特征速率測量平均絕對誤差率Table 2 AAEP of radar instantaneous vital sign rates detection from ten subjects at six different distances
表3 10名受試者在不同距離生理特征速率測量平均絕對誤差Table 3 AAE of radar instantaneous vital sign rates detection from ten subjects at six different distances
本文提出了一種基于小波分析和自相關計算的FMCW雷達非接觸式生理信號檢測方法. 該方法 基 于 小 波 包 分 解 (Wavelet packet decomposition,WPD)從原始信號中分解重構(gòu)心跳和呼吸信號,利用自相關計算減小雜波對心跳信號的影響,提高檢測精度. 從雷達原始信號中提取的呼吸和心跳信號與參考傳感器的信號相比具有很高的一致性. 本文選取不同受試者距離雷達不同距離進行非接觸式呼吸心跳檢測實驗. 隨著檢測距離的增加,信號信噪比的降低,非接觸式生理信號檢測難度漸漸增大. 本文方法主要考慮了環(huán)境噪聲及呼吸諧波對生理信號檢測的影響,而由運動引起的干擾強度較大且具有局部性、多發(fā)性及隨機性,會對生理信號檢測產(chǎn)生漂移、不規(guī)則、不連續(xù)等非線性變化,文中尚未考慮運動干擾對生理信號檢測的影響,我們將在未來的工作中進行更加深入的研究.