基于毫米波雷達(dá)的心音檢測(cè)

王 浩，張興敢

（南京大學(xué) 電子科學(xué)與工程學(xué)院，江蘇 南京 210023）

0 引 言

隨著毫米波雷達(dá)的發(fā)展，基于毫米波雷達(dá)的非接觸式生命信號(hào)檢測(cè)逐漸成為研究的熱點(diǎn)。毫米波雷達(dá)不僅可以用來連續(xù)檢測(cè)呼吸速率和心跳速率[1-2]，也可以對(duì)人體進(jìn)行睡眠監(jiān)測(cè)[3]。心音信號(hào)反映了心臟及主血管的機(jī)械運(yùn)動(dòng)，心臟的病變也會(huì)引起心音的病變，因此有效的心音信息可以輔助診斷心臟疾病[4]。心臟經(jīng)歷一次完整的收縮期和舒張期的時(shí)間稱作心動(dòng)周期，心音信號(hào)伴隨著心臟的跳動(dòng)重復(fù)出現(xiàn)。每個(gè)心動(dòng)周期出現(xiàn)的心音信號(hào)通常都包含有第一心音（S1）和第二心音（S2），但這兩種心音的產(chǎn)生原因完全不同。第一心音的產(chǎn)生是由于心室收縮時(shí)二尖瓣、三尖瓣的閉合；第二心音則是由于心室舒張時(shí)主、肺動(dòng)脈瓣的閉合。LFM（Linear Frequency Modulated）毫米波雷達(dá)可以根據(jù)瓣膜振動(dòng)引起的胸腔微小位移實(shí)現(xiàn)非接觸式的心音信號(hào)獲取。心音信號(hào)的精確分割是對(duì)S1、S2 時(shí)域上的精確定位，是心音分類[5]、異常心音識(shí)別[6]、雜音檢測(cè)[7]等心音信息分析的基礎(chǔ)，也是心音特征提取的關(guān)鍵。

心音信號(hào)的分割算法一直是心音研究的重點(diǎn)，文獻(xiàn)[8]采用基于閾值和K 類均值聚類的心電分析方法識(shí)別心音。文獻(xiàn)[9]通過心音信號(hào)的短時(shí)能量譜和自回歸參數(shù)提取信號(hào)的包絡(luò)，利用多層感知器神經(jīng)網(wǎng)絡(luò)和第一心音和第二心音的頻譜特征提取基本的心音，結(jié)合舒張和收縮間期變化完成兒童的心音分割。這些基于時(shí)域包絡(luò)特征的心音分割算法會(huì)受到提取包絡(luò)方法的影響，例如通過希爾伯特變換提取的包絡(luò)有較大的毛刺，從而影響心音分割的準(zhǔn)確度。之后又有學(xué)者把心音分割問題建模為隱馬爾可夫模型。文獻(xiàn)[10]在隱馬爾可夫模型的基礎(chǔ)上將心音信號(hào)的分割問題建模為連續(xù)時(shí)段隱馬爾可夫模型（Duration-dependent Hidden Markov Model,DHMM），利用測(cè)試數(shù)據(jù)集計(jì)算模型參數(shù)，完成心音信號(hào)的分割。文獻(xiàn)[11]將文獻(xiàn)[10]的DHMM 進(jìn)行了部分改進(jìn)，提出了基于邏輯回歸的隱半馬爾可夫模型（Logistic Regression-Hidden Semi-Markov Model, LR-HSMM），對(duì)心音信號(hào)分割的起始時(shí)間選擇更具有包容性，并采用拓展的維比特算法對(duì)LR-HSMM 進(jìn)行譯碼。文獻(xiàn)[12]提出的解碼峰值檢測(cè)算法與文獻(xiàn)[11]的LR-HSMM 相似，解決了在第二心音（S2）的干擾下，對(duì)心音信號(hào)中的第一心音（S1）的準(zhǔn)確定位問題?；陔[馬爾可夫模型的研究需要大量的數(shù)據(jù)樣本計(jì)算模型參數(shù)，導(dǎo)致算法比較復(fù)雜，分割效果受到模型參數(shù)限制。

本文以LFM 毫米波雷達(dá)為基礎(chǔ)，首先利用雷達(dá)差拍信號(hào)的相位信息獲得心臟收縮舒張時(shí)產(chǎn)生的胸腔位移，采用帶通濾波器完成對(duì)胸腔位移信息中心音信號(hào)的分離；接著利用基于雙特征的心音分割方法，以一階導(dǎo)數(shù)和頻率包絡(luò)作為心音信號(hào)特征，通過最大類間方差準(zhǔn)則計(jì)算閾值，完成心音信號(hào)中心音部分和非心音部分的初步分割；最后根據(jù)收縮期和舒張期的時(shí)間間隔的不同，對(duì)心音部分中第一心音和第二心音進(jìn)行類型判決，將一個(gè)心動(dòng)周期內(nèi)的心音信號(hào)分割為第一心音、收縮期（除去第一心音）、第二心音、舒張期（除去第二心音）四部分。本文通過毫米波雷達(dá)采集不同對(duì)象的心音信號(hào)作為實(shí)驗(yàn)數(shù)據(jù)，并與LR-HSMM 算法進(jìn)行對(duì)比。

1 心音信號(hào)檢測(cè)

1.1 LFM 毫米波雷達(dá)系統(tǒng)

整個(gè)LFM 毫米波雷達(dá)系統(tǒng)主要分為兩個(gè)子系統(tǒng)：一個(gè)是負(fù)責(zé)雷達(dá)信號(hào)發(fā)送和接收射頻模擬子系統(tǒng)；另一個(gè)是負(fù)責(zé)對(duì)差頻信號(hào)進(jìn)行信號(hào)處理的數(shù)字子系統(tǒng)。射頻模擬子系統(tǒng)主要包含發(fā)送和接收天線、模數(shù)轉(zhuǎn)換器、混頻器以及用于生成線性調(diào)頻脈沖的頻率合成器。數(shù)字子系統(tǒng)主要包含數(shù)字信號(hào)處理器、微型控制器。圖1為L(zhǎng)FM 毫米波雷達(dá)系統(tǒng)的一個(gè)簡(jiǎn)易模型，其中LNA 表示低噪聲放大器，PA 表示功率放大器。

圖1 LFM 毫米波雷達(dá)系統(tǒng)框圖

LFM 毫米波雷達(dá)的發(fā)射信號(hào)為：

式中：f0為發(fā)射信號(hào)的起始頻率；S為發(fā)射信號(hào)的調(diào)頻斜率；T為發(fā)射脈沖信號(hào)的周期。在不考慮雷達(dá)散射截面積的情況下，毫米波雷達(dá)接收到的人體反射回波信號(hào)為：

式中：td表示雷達(dá)發(fā)射信號(hào)到達(dá)接收端的時(shí)延，td= 2Rc，c 為光速，R表示人體胸部與雷達(dá)之間的距離。由于在毫米波雷達(dá)的工作過程中人體是保持靜止的，因此R=R0+ ΔR，R0為人體與毫米波雷達(dá)的距離，ΔR為胸腔運(yùn)動(dòng)位移。

發(fā)射信號(hào)sT(t)與接收信號(hào)sR(t)經(jīng)過接收端混頻器后可以得到差拍信號(hào)sb(t)，sb(t)可以表示為：

由式（3）可以得到差拍信號(hào)sb(t)的頻率fb=Std=2RSc，相位φb為：

式中：λ為毫米波雷達(dá)發(fā)射的電磁波波長(zhǎng)，λ= cf0。由式（5）可知，差拍信號(hào)sb(t)的相位變化Δφb和人體胸腔位移ΔR線性相關(guān)，差拍信號(hào)sb(t)的相位變化可以代表人體的胸腔運(yùn)動(dòng)。

1.2 心音信號(hào)提取

心音信號(hào)同心跳信號(hào)一樣都屬于微弱的生命信號(hào)，而心音信號(hào)的檢測(cè)與心跳信號(hào)檢測(cè)的最大區(qū)別在于兩者振動(dòng)的頻率不同。借助差拍信號(hào)sb(t)的相位變化確定人體胸腔的運(yùn)動(dòng)位移后，根據(jù)心音信號(hào)所在的頻率區(qū)間完成心音信號(hào)的提取，心音信號(hào)的提取流程如圖2所示。

圖2 心音信號(hào)的提取流程圖

接收端混頻后得到的差拍信號(hào)sb(t)經(jīng)過A/D 采樣和距離快速傅里葉變換后，得到一維距離像向量，代表快時(shí)間維。多個(gè)差拍信號(hào)經(jīng)過相同處理后得到的多個(gè)距離像向量構(gòu)成距離-時(shí)間二維矩陣W，矩陣的列代表每個(gè)距離單元的慢時(shí)間采樣。確定人體所在的距離單元，需要計(jì)算二維矩陣W中每個(gè)一維距離像向量的功率，再計(jì)算二維矩陣每列的平均功率，以最大功率平均值所對(duì)應(yīng)的距離單元作為人體所在距離單元，胸腔位移運(yùn)動(dòng)信息就包含在與人體距離單元數(shù)據(jù)的相位內(nèi)。如圖3 所示（圖3 中的功率均歸一化），平均功率最大值對(duì)應(yīng)的距離為0.39 m。

圖3 人體目標(biāo)檢測(cè)

在距離-時(shí)間二維矩陣W中的數(shù)據(jù)都為復(fù)數(shù)，可以通過反正切解調(diào)的方法獲得矩陣W中人體目標(biāo)所在距離列每一個(gè)數(shù)據(jù)點(diǎn)的相位。但是在LFM 毫米波雷達(dá)系統(tǒng)中存在一些非線性因素會(huì)影響反正切解調(diào)，其中最主要的是未知的直流偏移，而直流偏移的主要來源是靜態(tài)物體反射的電磁波。假設(shè)一個(gè)復(fù)距離單元數(shù)據(jù)的虛部為I(t) + DCi，實(shí)部為R(t) + DCr，那么解調(diào)所得的相位φ(t)為：

式中：DCr、DCi分別為實(shí)部和虛部的直流項(xiàng)；φ(t) =I(t)R(t)，為真實(shí)相位值。為了消除直流項(xiàng)的影響，本文采用Spath 算法[13]估計(jì)直流項(xiàng)，并通過補(bǔ)償?shù)姆绞较龜?shù)據(jù)中的直流偏移，直流偏移補(bǔ)償前后結(jié)果如圖4所示。

圖4 人體距離列數(shù)據(jù)的直流補(bǔ)償

當(dāng)相位大于π 或小于-π 時(shí)，會(huì)被相位解調(diào)函數(shù)投影在(-π,π)范圍內(nèi)，致使相位突變。通過相位解卷繞可以實(shí)現(xiàn)相位的連續(xù)變化，解卷繞后的相位變化便等同于胸腔的運(yùn)動(dòng)變化。由于心跳產(chǎn)生的心音信號(hào)非常微弱，需要對(duì)解卷繞后的相位進(jìn)行差分運(yùn)算，增強(qiáng)心音信號(hào)，即zk=φk-φk-1。心音信號(hào)中第一心音和第二心音的頻率范圍[14]大致為10～200 Hz，本文選用IIR 帶通橢圓濾波器對(duì)差分后的相位進(jìn)行濾波，分離心音信號(hào)。為避免人體低頻呼吸信號(hào)的干擾，同時(shí)考慮不同個(gè)體之間的普適性，帶通濾波器的頻帶范圍定為20～60 Hz，分離出的心音信號(hào)如圖5 所示。

圖5 LFM 毫米波雷達(dá)采集的心音信號(hào)

2 基于雙特征的心音分割

2.1 雙特征閾值分割

在時(shí)域方面，假設(shè)心音信號(hào)為s(t)，則s(t)的一階導(dǎo)數(shù)為：

當(dāng)Δt較小時(shí)，可近似為：

式中Ts為心音信號(hào)s(t)的采樣周期。s′(t)取絕對(duì)值之后得到心音信號(hào)s(t)的時(shí)域特征sp(t)，從圖5 心音信號(hào)波形圖中可以發(fā)現(xiàn)，由振動(dòng)產(chǎn)生的第一心音和第二心音信號(hào)波形會(huì)產(chǎn)生明顯的上下波動(dòng)，且第一和第二心音的幅度明顯高于非心音部分。根據(jù)文獻(xiàn)[10]，第一心音和第二心音的持續(xù)時(shí)間小于收縮期（除去第一心音)和舒張期（除去第二心音），對(duì)應(yīng)心音信號(hào)中S1 和S2 部分的時(shí)域特征sp(t)的值會(huì)明顯高于非心音部分。在時(shí)頻方面，通過短時(shí)傅里葉變換（Short Time Fourier Transform,STFT）對(duì)心音信號(hào)s(t)進(jìn)行時(shí)頻分析：

式中h(t)為窗函數(shù)，窗長(zhǎng)的選擇需要滿足兩個(gè)條件：一是要保證具有良好的時(shí)間分辨率；二是盡量不超過第一心音或是第二心音的持續(xù)時(shí)間，本文將窗長(zhǎng)設(shè)定為100 ms。通過心音信號(hào)s(t)的短時(shí)傅里葉變換，得到心音信號(hào)頻帶范圍為[f1,f2]的頻率包絡(luò)FEnv(t)為：

式中：頻率f的間隔為Δf；頻率包絡(luò)FEnv(t)代表[f1,f2]內(nèi)頻率能量之和隨時(shí)間的變化。心音信號(hào)中屬于S1 和S2 時(shí)間段的頻率能量會(huì)高于非心音時(shí)間段。本文設(shè)定了三種心音信號(hào)的頻率包絡(luò)特征：FEnv1、FEnv2 和FEnv3，它們代表了三種不同頻帶范圍的頻率包絡(luò)，F(xiàn)Env1 的 頻 率 范 圍 為20～30 Hz，F(xiàn)Env2 的 頻 率 范 圍 為30～40 Hz，F(xiàn)Env3 的頻率范圍為40～60 Hz。

心音信號(hào)的時(shí)域特征sp(t)和時(shí)頻特征FEnv(t)構(gòu)成了心音信號(hào)的一組特征向量，即gt=[sp,FEnv]。通過設(shè)定閾值，將心音信號(hào)在時(shí)間維度上初步分為兩個(gè)部分：第一心音或第二心音部分以及非心音部分。令θ表示閾值向量，即θ=[θ1,θ2]，若sp(t) ＞θ1且FEnv(t) ＞θ2，即gt＞θ，則時(shí)刻t歸為第一心音或第二心音部分；否則，時(shí)刻t歸為非心音部分，其中θ1為特征向量sp對(duì)應(yīng)的閾值，θ2為特征向量FEnv 對(duì)應(yīng)的閾值。

2.2 計(jì)算閾值向量θ

計(jì)算閾值向量θ等價(jià)于計(jì)算閾值θ1和θ2。本文利用最大類間方差法[15]求解閾值θ1和閾值θ2（閾值θ2計(jì)算步驟與θ1相同）。θ1的計(jì)算過程如下：

1）對(duì)向量sp中的元素進(jìn)行量化，所有元素值根據(jù)四舍五入原則只取小數(shù)點(diǎn)后三位，量化后的向量sp表示為sp=[sp,1,sp,1,…,sp,L]，其 中L為向量sp的元素個(gè)數(shù)。向量sp中元素值為sp,i的元素個(gè)數(shù)為ni，計(jì)算元素值為sp,i的元素出現(xiàn)的概率pi為：

初始化閾值θth，令θth= mean(sp)，其中函數(shù)mean(·)為求均值函數(shù)。

2）根據(jù)閾值θth可以將向量sp中的元素分為兩部分：一部分為元素值大于θth的S1 或S2 部分；另一部分為元素值小于θth的非心音部分。計(jì)算sp中元素屬于S1或S2 部分的概率P1以及非心音的概率P2：

3）計(jì)算sp中屬于S1 或S2 部分的元素平均值μ1以及非心音部分的元素平均值μ2：

閾值為θth時(shí)，類間方差為：

4）令θth=θth+ 0.001，若θth大于sp的最大元素值則停止，否則，轉(zhuǎn)至步驟2）。

經(jīng)過以上步驟得到閾值θ1：

2.3 心音類型識(shí)別

圖6 為經(jīng)過2.1 節(jié)的方法處理后得到的簡(jiǎn)易心音信號(hào)時(shí)域分布模型。模型中的心音信號(hào)在時(shí)間維度被分為了不同長(zhǎng)度的信號(hào)片段，其中黑色的心音信號(hào)片段代表第一心音或第二心音，白色部分代表收縮期（除去第一心音）或舒張期（除去第二心音）的非心音部分。只有確定黑色心音片段的心音類型之后才能完成心音信號(hào)的最終分割。根據(jù)文獻(xiàn)[16]，一次完整的心動(dòng)周期中心室收縮期持續(xù)時(shí)間一般小于舒張期，且第一心音和第二心音的起始分別對(duì)應(yīng)著心室收縮期和舒張期的起始。根據(jù)圖6 的心音信號(hào)時(shí)域分布模型，心音類型識(shí)別方法為：若tstart,i-tstart,i-1＜tstart,i+1-tstart,i，則shs,i為第一心音；若tstart,i-tstart,i-1＞tstart,i+1-tstart,i，則shs,i為第二心音，其中shs,i為第i個(gè)S1 或S2 部分信號(hào)片段，tstart,i為第i個(gè)S1 或S2 部分信號(hào)片段的起始時(shí)刻。心音類型識(shí)別完成后，一個(gè)心動(dòng)周期內(nèi)心音信號(hào)可以分割為4 種連續(xù)的時(shí)間階段：S1、Sys、S2、Dia，它們分別對(duì)應(yīng)第一心音、除去第一心音之外的收縮期、第二心音、除去第二心音之外的舒張期。

圖6 心音信號(hào)時(shí)域分布模型

3 實(shí)驗(yàn)結(jié)果及分析

本文使用的毫米波雷達(dá)為德州儀器公司的AWR1642，工作頻段為77～81 GHz。在PC 端完成配置雷達(dá)參數(shù)后，毫米波雷達(dá)差拍信號(hào)采樣數(shù)據(jù)會(huì)通過DCA1000 采集卡傳回電腦，在PC 端通過Matlab R2020b進(jìn)行信號(hào)處理和心音分割。實(shí)驗(yàn)通過6 位不同的測(cè)試者采集心音信號(hào)數(shù)據(jù)，在實(shí)驗(yàn)過程中，測(cè)試者在近距離面對(duì)雷達(dá)而坐，同時(shí)毫米波雷達(dá)高度與測(cè)試者的心臟高度保持一致。毫米波雷達(dá)的參數(shù)配置為：調(diào)頻斜率S為60.012 MHz·μs-1，脈沖起始頻率f0為77 GHz，ADC 采樣率fs為2 MHz·μs-1，脈沖持續(xù)時(shí)間T為50 μs，每條數(shù)據(jù)的采集總時(shí)長(zhǎng)為10 s。

3.1 心音分割結(jié)果

對(duì)于每一條心音信號(hào)數(shù)據(jù)，根據(jù)2.1 節(jié)的內(nèi)容，通過式（8）和式（10）分別提取心音信號(hào)s(t)的一階導(dǎo)數(shù)s′(t)和頻率包絡(luò)FEnv(t)，之后對(duì)s′(t)取絕對(duì)值得到特征向量sp和FEnv，心音信號(hào)的時(shí)域特征sp(t)和三種頻率包絡(luò)如圖7 所示（三種頻率包絡(luò)均歸一化）。從圖7b）中可以看出，與原始的心音信號(hào)相比，sp(t)放大了S1 和S2 部分與非心音部分之間的幅度差距，圖7c）中頻率包絡(luò)FEnv(t)的變化也與心音信號(hào)的幅度包絡(luò)變化一致，即S1 和S2 部分的頻率能量高于非心音部分，將兩種特征結(jié)合可以使得S1 和S2 部分與非心音部分之間的分割更加可靠。

圖7 心音信號(hào)的時(shí)域特征和三種頻率包絡(luò)

根據(jù)2.2 節(jié)的步驟計(jì)算特征向量sp和FEnv 對(duì)應(yīng)的閾值θ1和θ2，若sp(t) ＞θ1且FEnv(t) ＞θ2，則時(shí)刻t歸為第一心音或第二心音部分；否則，時(shí)刻t歸為非心音部分。利用閾值向量θ完成對(duì)心音信號(hào)s(t)中S1 或S2 部分和非心音部分在時(shí)間軸上的分割，如圖8 所示。此時(shí)，心音信號(hào)被劃分為持續(xù)時(shí)間不等的S1 或S2 部分以及非心音部分的信號(hào)片段，由于臨床上心室的收縮期一般小于心室的舒張期，所以根據(jù)2.3 節(jié)心音類型的識(shí)別方法，將S1 或S2 部分的信號(hào)片段進(jìn)行第一心音或第二心音的判決，最終將整個(gè)心音信號(hào)分割為4 種時(shí)間階段：S1、Sys、S2、Dia。頻率包絡(luò)選定為FEnv1 時(shí)，心音分割結(jié)果如圖9 所示。從圖9 可以看出，本文所提出的算法對(duì)于心音信號(hào)的分割定位取得了很好的效果。

圖8 雙特征閾值分割

圖9 心音信號(hào)分割

3.2 心音分割效果對(duì)比

為了準(zhǔn)確評(píng)估所提算法對(duì)于心音信號(hào)的分割效果，以手動(dòng)分割心音信號(hào)的結(jié)果作為參考標(biāo)準(zhǔn)，將本文所提算法與文獻(xiàn)[11]中的LR-HSMM 算法進(jìn)行對(duì)比，選取評(píng)價(jià)指標(biāo)F1分?jǐn)?shù)來評(píng)估兩種算法，F(xiàn)1分?jǐn)?shù)是衡量測(cè)試準(zhǔn)確性的常用指標(biāo)，它同時(shí)包含了精確率p和召回率r，計(jì)算公式為：

根據(jù)本文方法，6 位測(cè)試者心音分割的p、r和F1分?jǐn)?shù)結(jié)果如表1 所示。兩種算法的指標(biāo)對(duì)比結(jié)果如表2 所示。表1 和表2 中結(jié)果為所有數(shù)據(jù)指標(biāo)的平均值，黑體字表示最優(yōu)F1分?jǐn)?shù)，其中LR-HSMM 算法中的特征提取部分選用同態(tài)包絡(luò)。從表1 中6 位不同測(cè)試者的平均F1分?jǐn)?shù)可知，當(dāng)頻率包絡(luò)選取FEnv1 時(shí)，本文方法對(duì)于心音信號(hào)的分割效果明顯優(yōu)于FEnv2 和FEnv3。由表2 可知，當(dāng)頻率包絡(luò)為FEnv1時(shí)，無論是第一心音（S1）還是第二心音（S2），本文方法在p、r以及F1分?jǐn)?shù)三項(xiàng)指標(biāo)上均高于文獻(xiàn)[11]的LR-HSMM 算法。其中對(duì)于第一心音的F1分?jǐn)?shù)，本文方法達(dá)到了（92.29±0.94）%，比LR-HSMM 算法提高了7.07%；對(duì)于第二心音的F1分?jǐn)?shù)，本文方法達(dá)到了（88.71±4.61）%，比LR-HSMM 算法提高了7.76%。

表1 心音分割的F1 分?jǐn)?shù) %

表2 兩種算法的F1 分?jǐn)?shù)對(duì)比 %

4 結(jié) 論

本文利用LFM 毫米波雷達(dá)獲得人體胸腔位移信息，通過帶通濾波提取人體的心音信號(hào)。利用心音信號(hào)的一階導(dǎo)數(shù)和頻率包絡(luò)特征，同時(shí)結(jié)合心音間隔時(shí)間的生理特征完成對(duì)心音信號(hào)的分割定位。實(shí)驗(yàn)結(jié)果表明，本文提出的算法對(duì)于心音信號(hào)的分割定位具有較高的精確度，為后續(xù)的心音記錄自動(dòng)分析奠定了基礎(chǔ)。下一步的研究將包括心音在內(nèi)諸多人體身體特征，通過毫米波雷達(dá)應(yīng)用于個(gè)人身份的識(shí)別，從而進(jìn)一步擴(kuò)大心音的用途。

注：本文通訊作者為張興敢。