楊國成,余慧敏
(湖南師范大學 信息科學與工程學院,湖南 長沙 410006)
超寬帶雷達通過測量人體身體微動(如呼吸和心跳運動)導致的雷達接收信號脈沖時延的變化檢測人體生命信息[1-3]。關于生命信息檢測,國內外研究小組一直在開發(fā)雷達原型,并提出了許多不同的呼吸和心跳檢測算法[4-5]。文獻[6]闡述一些信號去噪方法,如平均對消法來消除靜態(tài)雜波,然后提取已濾波矩陣的最高頻譜分量作為呼吸頻率。文獻[7]將檢測到的特征時間位置從胸部表面的最大回聲處移動到接近肺后部的組織層,以避免錯誤提取其他反射信號。盡管通過這些方法可以獲得呼吸頻率,但仍然難以準確檢測到心跳信號。另外,存在一些研究利用EMD 或EEMD 分解信號以分離系統(tǒng)噪聲和生命信號[8-10]。文獻[11]將EMD 用于處理回波數(shù)據(jù),并提前預設采樣時間為攜帶呼吸和心跳信息的14 ~16 ns,但是由于測量環(huán)境會發(fā)生變化,僅僅憑借先前的經驗難以確定呼吸和心跳所在的位置。文獻[12]將EEMD 和希爾伯特黃變換用于處理回波數(shù)據(jù),能夠識別呼吸信號但識別不了心跳信號。由于呼吸運動的強度更大,且呼吸和心跳的頻帶范圍接近,心跳信號總是被呼吸諧波和其他雜波所覆蓋,僅憑現(xiàn)有方法很難從低信噪比的環(huán)境中檢測到心跳信息。為了解決上述問題,本文提出了一種基于區(qū)域谷值雙層EEMD 的超寬帶雷達生命信號檢測方法。
1.1.1 人體胸部運動模型
人體的呼吸和心跳運動可近似為簡諧振動,則雷達天線到人體胸部的瞬時距離為:
其中,t 為慢時間,d0為天線到胸部振動中心的平均距離,Ar和Ah分別為呼吸和心跳的位移幅度,fr和fh分別為呼吸頻率和心跳頻率。天線到人體胸部的距離變化導致接收信號的延遲變化,接收信號的時間延遲為:
其中,V0為光速。
1.1.2 信道模型
對于理想的多徑傳播信道,脈沖響應h(t,τ)為人體生命體征信號的響應和其他靜止目標的響應之和:
1.1.3 接收機模型
接收信號一般為發(fā)射脈沖和信道脈沖響應的卷積加上系統(tǒng)噪聲,表示為:
其中,R(t,τ)為接收信號,s(t,τ)為發(fā)射信號,n(t,τ)為系統(tǒng)噪聲。
偽二維聚類經驗模態(tài)分解首先使用聚類經驗模態(tài)分解沿特定時間軸分解數(shù)據(jù)矩陣的各個列或行信號,然后通過相應的級別將分解的本征模態(tài)函數(shù)重新排列成不同的二維信號矩陣[13]。
已知y(t,τ)是二維原始數(shù)據(jù)矩陣,同時是一個快慢時間場,τ 代表快時間的N 個時間位置,t 代表慢時間的M 個時間位置:
首先,PBDEEMD 分解將對y(t=m,τ)的數(shù)據(jù)矩陣(y 的第個m 行向量)進行EEMD 分解,以獲得一系列分解結果cj,如下:
其中,cj=(cm,1,jcm,2,j…cm,N,j),1 ≤j ≤J,cj是第j 個IMF 的行向量,J 是IMF 的數(shù)量。
其次,將在慢時間坐標t=m+1 處的數(shù)據(jù)矩陣進行EEMD 分解,將相應第j 個比例級別的結果重新排列為M 行N 列的cj(t,τ)矩陣。矩陣cj(t,τ)就是第j個偽二維本征模態(tài)函數(shù)PBDIMF,表示為:
PBDEEMD 分解也可以依次對列向量進行EEMD 分解,例如y(t,τ=n),按照類似的步驟進行。
基于區(qū)域谷值雙層EEMD 的信號檢測方法可分為3 步。第一,對回波數(shù)據(jù)進行PBDEEMD 分解,獲得攜帶生命信息的目標矩陣;第二,計算目標矩陣關于快時間的能量函數(shù),選擇能量函數(shù)中目標時間區(qū)域的谷值作為特征時間指數(shù)(Feature Time Index,F(xiàn)TI);第三,提取特征時間指數(shù)所對應的慢時間信號,并對慢時間信號進行EEMD 分解,分離出呼吸和心跳信號。該方法可以有效濾除回波數(shù)據(jù)中的靜態(tài)雜波和其他噪聲,并實現(xiàn)呼吸信號和心跳信號的有效分離,具體流程如圖1 所示。
圖1 具體流程
(1)對回波數(shù)據(jù)進行PBDEEMD 分解,以濾除系統(tǒng)靜態(tài)雜波。然后,選擇第一級PBDIMF 矩陣即cj(t,τ)作為攜帶生命信息的目標矩陣,并將第2級到第8 級PBDIMF 矩陣作為靜態(tài)雜波進行濾除。
(2)計算目標矩陣PBDIMF 在快時間上的能量函數(shù),選擇能量函數(shù)第一峰值和第二峰值之間區(qū)域作為目標時間區(qū)域,并選擇目標時間區(qū)域內的谷值作為特征時間指數(shù)。
(3)提取目標矩陣PBDIMF 在FTI 處的慢時間信號,使用EEMD 將信號分解為若干個IMF,并對所有IMF 進行快速傅里葉變換獲得頻譜圖,依據(jù)呼吸和心跳頻段分離出高頻雜波和低頻呼吸心跳信號。
能量函數(shù)ej(τ)是將目標矩陣PBDIMF 的能量按慢時間累加獲得,表示為:
為簡單起見,在下文中使用e(τ)表示ej(τ)。
根據(jù)無線電反射原理,在信道模型的傳播路徑上,回波波形的最大值具有最大能量即能量函數(shù)第一峰值,該處時間指數(shù)為τpeak1;回波波形的第二大值具有第二大能量即能量函數(shù)第二峰值,該處時間指數(shù)為τpeak2。能量函數(shù)第一峰值和第二峰值之間的時間區(qū)域定義為目標時間區(qū)域。由于回波信號在該區(qū)域中的幅度較小,即在該區(qū)域內的能量強度也較弱。目標時間區(qū)域的谷值即目標時間區(qū)域中的能量最小值,表示如下:
其中,時間指數(shù)τvalley是目標時間區(qū)域中的能量最小值即特征時間指數(shù),此處的能量強度最小,適合提取信號中的生命信息?;夭〝?shù)據(jù)的能量函數(shù)如圖2 所示,第一峰值的時間指數(shù)τpeak1為99,第二峰值的時間指數(shù)τpeak2為101,兩峰值間的谷值(τvalley=100)就是特征時間指數(shù)FTI。
圖2 回波矩陣的能量函數(shù)
雷達回波的脈沖波形具有的能量越多,檢測到呼吸信號的機會就越大。然而,在這樣的條件下仍然很難找到心跳信號,可以認為系統(tǒng)噪聲或靜態(tài)雜波掩蓋了微小的心跳信號。實際上,心跳總是存在但很難準確測量。基于區(qū)域谷值的雙層EEMD 方法的優(yōu)勢在于,與τFTI之外的時間比,在τFTI處的慢時間信號的波形幅度很小,而生命活動的幅度相對較大。因此,利用雜波和生命信號之間的幅度差異進行分離,導致呼吸和心跳信號可以通過不同的頻率特性區(qū)分。
本節(jié)將通過MATLAB 仿真驗證該方法的有效性,此處只考慮空氣胸腔界面的一次反射,并添加靜態(tài)雜波和高斯噪聲以模擬復雜的現(xiàn)實環(huán)境。其中,呼吸和心跳的頻率為0.3 Hz、1.1 Hz;呼吸和心跳的位移幅度為2.3 mm、0.5 mm,快時間采樣間隔為0.05 ns,快時間采樣點數(shù)為300,脈沖發(fā)射重復頻率為100 Hz,總發(fā)射時間為20 s,模擬信號SNR 范圍為2 ~20 dB,仿真參數(shù)詳具體如表1 所示?;夭〝?shù)據(jù)模型的建立過程如圖3 所示。
表1 模擬參數(shù)
圖3 回波信號的傳播模型
UWB 脈沖是一個非常短的時間脈沖,接收器獲得的噪聲主要來自雷達系統(tǒng)本身,因此在本節(jié)驗證該方法能在不同SNR 下提取呼吸和心跳信號。在仿真中添加不同SNR(2 dB、5 dB、10 dB、20 dB)的高斯白噪聲,并用如上所述的區(qū)域谷值方法確定特征時間指數(shù)獲得τvalley=100。仿真結果數(shù)據(jù)如圖4所示,對于SNR=10 dB 和20 dB,IMF5 和IMF6 的主要頻率峰值分別等于預先設定的呼吸和心跳速率,分別為1.1 Hz 和0.3 Hz,測量值與設定值一致。圖4(a)和圖4(b)中,當SNR 等于2 dB 和5 dB時,由于強噪聲增加了提取心跳信號的難度,IMF5中的心跳頻率漂移到1.15 Hz,而IMF6 中的呼吸頻率仍保持在0.3 Hz。仿真結果表明,在不同信噪比的環(huán)境下,該方法都可以準確測量檢測目標的呼吸和心跳頻率,具有很強的抗干擾能力。
圖4 各種SNR 下的IMF 以及相應的呼吸和心跳速率
為了解決接收信號中心跳信息難以準確測量的問題,本文提出一種基于區(qū)域谷值雙層EEMD 的生命信號檢測方法。相比于已有的濾波矩陣中最高頻譜分量判作呼吸速率,該方法使用目標時間區(qū)域的谷值來確定特征時間指數(shù),并且將慢時間信號分解成本征模態(tài)函數(shù),進而有效地檢測過去被忽略具有小能量的信號,實現(xiàn)了呼吸與心跳信號的有效分離。仿真結果表明,在不同信噪比環(huán)境下,該方法都可以準確檢測出呼吸心跳頻率。