吴佳茜 梁步阁* 杨德贵 熊明耀 李元烽
(1.中南大学自动化学院,湖南长沙 410083;2.光电智能测控湖南省重点实验室,湖南长沙 410083)
随着医疗技术的显著发展以及人们对自身健康关注的加强,人体睡眠健康监测技术的研究热度越来越高。呼吸是衡量健康状况的重要指标,当呼吸发生异常时,一定程度上预示着疾病或其他危重情况,例如:肥胖、抑郁症、阻塞性睡眠呼吸暂停综合征、呼吸中枢衰竭等[1-2]。因此呼吸信号的监测和分析在睡眠健康应用中发挥着重要的作用,尤其是对异常情形下的呼吸障碍的捕捉与识别,能为睡眠监测、疾病预断、术后观察等提供重要依据。
生物雷达是一种新兴的非接触式探测技术,能全天候全天时、无创地、不侵犯隐私地实现体征监测,在生命体征检测中有着愈发广泛的应用,主要分为连续波(Continuous Wave,CW)、调频连续波(Frequency Modulated Continuous Wave,FMCW)和IR-UWB 三种体制[3-7]。其中,CW 雷达结构简单、集成度高,是目前使用较多的体征监测雷达,在呼吸辨识应用上主要通过提取呼吸波形的时域、频域、能量域的统计特征[8-10]及呼吸气流量的波形几何特征[11-12]来进行呼吸模式的分类和不同个体之间的身份识别,但CW 雷达不具备距离分辨能力,存在无法实现目标定位的缺点。FMCW 雷达在CW 体制的基础上解决了目标定位问题,更具有多目标生命体征检测的能力,也逐渐广泛应用于呼吸监测和分类,在呼吸特征的提取上,除了基于常规的统计特征[13-14]研究者还探索了MFCC 特征[15]、呼吸波形的物理运动特征和基于生理信息的动态分割特征[16],但FMCW 雷达存在系统成本高、实时性较差的缺点。相较之下,IR-UWB 雷达穿透能力强、距离分辨率高、抗干扰能力强,同时具有结构简单、功耗低、成本低廉等优点,尤其适合应用于睡眠健康监测领域[17-20],成为了目前非接触式生命体征监护领域的研究热点。但在呼吸模式分类研究上大多也仍然是通过提取一维呼吸信号的振幅和频率的统计特征来检测呼吸异常事件[21-23],而在超宽带雷达探测中人体呼吸引起的胸腔起伏在回波上不是简单的点目标,而是存在多个散射中心的距离扩展目标,产生的雷达回波将分布于多个距离单元中,蕴含着更加丰富的特征信息,共同表征着呼吸运动。文献[24]、[25]提取了IR-UWB 雷达跨距离单元的呼吸回波信号构建呼吸样本,分别实现了基于神经网络的人体身份识别和呼吸速率、呼吸暂停的检测,证明了挖掘跨距离单元呼吸信号中信息的可行性。
本文基于IR-UWB 雷达系统提出了一种引入时距信息的多域特征融合呼吸模式识别方法,算法首先提取呼吸距离门范围,对雷达回波进行杂波滤除,实现微弱呼吸的信号增强;然后在时间维上选取一维呼吸信号提取信号波形的时域、频域特征,同时根据呼吸在IR-UWB 雷达回波图像上呈现出的规律周期特性,在距离维上选取包含多个距离门信号在内的二维时距图像提取呼吸模式的形态特征。在图像特征提取中,针对呼吸的雷达图像受呼吸异常节律影响往往呈现包络放缩和局部粘连特性的问题,本文提出了一种通过相位矩阵图像处理获取呼吸时距条带的方法,去进一步提取条带图像中潜在的呼吸模式形态特征。基于雷达目标信号在时间和距离上的多域特征融合信息,实现不同呼吸模式的分类识别。
人体呼吸是一种具有周期性的生理行为,伴随着胸腔的扩张和收缩起伏,这些周期性运动引起的位移会对雷达回波信号产生不同的调制效果,从而可以在距离上分辨出来,并以此来检测人体体征信息。
本文实验使用的是课题组自研的KUN-DC-1000 脉冲超宽带雷达,发射机产生的频移高斯脉冲信号可表示为
其中,τ为快时间,p(τ)为高斯脉冲包络,fc为中心频率。
雷达接收机接收到的回波信号是由探测场景内多散射点的回波叠加构成的。假设探测场景中存在M个散射点,忽略噪声的影响,接收机前端接收到的反射信号可表示为
其中,t为慢时间,αi(t)为接收到的散射点的雷达散射截面积所决定的信号幅度,τi(t)为散射点与雷达之间的瞬时距离所决定的时间延迟。
对于人体呼吸运动而言,将呼吸引起的胸腔周期性起伏近似为简谐运动进行简化,则τ(t)可表示为
其中,d0是接收机与胸腔振动中心之间的距离,a为呼吸运动振幅,fb为呼吸频率,c为光速。
接收信号经过混频和低通滤波处理后,得到回波信号的表达式如下
对回波信号进行采样,以离散形式存储在二维回波矩阵R[m,n]中,可表示为
其中,m和n分别是快时间和慢时间采样数,Tf和Ts分别为快时间和慢时间采样间隔。根据式(5)对回波矩阵进行相应的信号处理和数据分析,则可以提取到不同距离维度的呼吸运动信息。
但睡眠监测雷达真实的测量环境往往是复杂的,除了人体目标,还存在着许多如墙壁、床、椅子、风扇等的干扰杂物,雷达回波中将混杂各种反射杂波。图1显示了分别在复杂和较为干净的室内环境下实测的回波矩阵,可以看到相较于干净的测量环境而言复杂环境下人体回波信号被大量杂波覆盖,体征辨识困难,这使得对呼吸信息的检测和特征提取更具有挑战性。
基于人体呼吸信号监测和呼吸模式识别在睡眠监测方面的主要应用,本文将对六种不同的呼吸模式进行分类识别,包括正常呼吸和五类典型的异常呼吸模式:潮式呼吸、间停呼吸、叹气样呼吸、比奥呼吸和库式呼吸。
正常呼吸是均匀、规则的呼吸运动,频次一般为12~20 次/分钟;潮式呼吸是一种由浅慢逐渐加深加快又逐渐变减弱变缓最后暂停一段时间,似潮水涨退的一种呼吸模式,常见于脑炎、糖尿病等危重病人;间停呼吸是指正常呼吸时突然停止一段时间又继续呼吸,如此反复进行,常见于睡眠呼吸暂停情况以及更加严重的颅内病变;叹气样呼吸表现为正常呼吸节律中插入一次深大呼吸,并伴有叹气声,常见于焦虑症、神经衰弱或更为病危的昏迷濒死情形;比奥呼吸的特点是呼吸急促同时伴随着呼吸暂停,常见于中风、创伤等中枢神经系统疾病;库式呼吸是急促深大的费力型呼吸,常见于剧烈运动、情绪激动的情形和糖尿病酮症酸中毒的糖尿疾病患者[1-2,26-27]。图2 显示了不同呼吸模式的一维雷达波形和呼吸区域回波图像。
图2 不同呼吸模式的一维信号波形和二维时距图像Fig.2 One-dimensional signal waveforms and two-dimensional time-range images of different respiratory modes
对于上述几种呼吸模式的分类识别,本文提出引入时距图像的IR-UWB 雷达呼吸模式多域特征精确提取方法,算法流程包括呼吸信号预处理和呼吸信号多域特征提取两个部分,算法流程图如图3所示。
图3 呼吸特征精确提取算法流程图Fig.3 Algorithm flowchart for precise extraction of respiratory features
3.2.1 呼吸距离门范围提取
常规的IR-UWB 雷达呼吸模式识别仅对单个距离门上的呼吸数据进行特征提取,但由于人体胸腔有一定的宽度存在多个散射点,且测量过程中人体存在随机晃动,目标回波会分布在相邻的多个距离门中。为了充分利用回波中的信息,首先需要提取呼吸回波距离门范围。距离门范围提取的关键在于中心距离门的寻找,本文从信号能量角度分析,非目标物体反射的回波信号具有时不变的特性,而呼吸时胸腔起伏反射的回波信号能量波动相对较大,波动越大信号能量的方差越大,因此可以采用最大方差法选取目标所在的中心距离门。具体步骤为:
(a)减均值处理:为了更加凸显目标呼吸信号的能量,削弱回波中固定杂波的能量,首先对原始回波矩阵做减均值处理,即用历史时间内的回波数据进行统计平均来估计每个距离下的杂波能量,在原始回波信号中减去。对于原始回波矩阵R[m,n]的减均值处理可表示为
式中,R′[m,n]是减均值处理后的数据矩阵,N为慢时间采样点总数。算法杂波能量去除效果对比如图4所示,目标呼吸的能量比重明显提高。
图4 减均值处理前后对比Fig.4 Comparison of results before and after mean reduction treatment
(b)最大方差法距离门寻优:基于方差分析选择信号能量方差最大的快时间距离门作为最优中心距离门,信号能量方差随快时间距离门的变化曲线如图5所示。对于第m个距离门信号能量的方差σ2(m)和最优距离门l的计算方式为
图5 信号能量方差随快时间距离门的变化曲线Fig.5 Curve showing signal energy variance with fast time distance gate
(c)呼吸距离门范围提取:选取中心距离门附近几个距离门作为呼吸信号子矩阵的距离门范围,记为l-l0~l+l0,其中l0的计算与人体胸部宽度和雷达的距离分辨率有关,本文实验中设定l0=2。图6显示了根据距离门范围分割出来的呼吸区域子图,能清晰地显现出呼吸周期。
图6 呼吸区域子图分割Fig.6 Respiratory region subgraph segmentation
3.2.2 子矩阵呼吸相位增强
根据呼吸距离门范围划分原始回波矩阵,得到呼吸回波子矩阵,对其进行信号处理,提取其中呼吸分量的相位信息,并进行增强处理,主要过程包括静态杂波滤除、反正切解调、直流分量和线性趋势去除、滤波四个步骤。
测量环境中摆放的杂乱物体反射回波会产生许多个静态信号分量,且幅度远远大于目标信号分量,则对于同一距离门式(4)可改写为
其中,第一项是静态信号分量,As和θs是这些静态物体反射回波的等效幅度和相位;第二项是目标信号分量,A(t)是信号幅度,θ0是人体躯干中心位置与雷达存在一定距离所产生的相位,x(t)表示胸腔位移[28]。
用复平面中的旋转矢量来表示式(9),如图7所示,理想状态下目标分量Pm在复平面上的轨迹可以近似为以原点为圆心的一段圆弧,而静态分量Ps的存在会使得目标分量偏离原点,从而导致随目标心肺运动变化的相位测量值Δφr往往小于真实值Δφ,且静态杂波分量越大,偏差越大。
图7 静态杂波去除原理示意图Fig.7 Schematic diagram of static clutter removal principle
为了准确提取目标信息,采取圆拟合方法来消除静态杂波分量的影响,即通过求解非线性最小二乘几何拟合问题来估计复平面中数据圆的圆心位置和半径长度,将圆心“平移”到原点。目标函数可写为
其中N为数据点的个数,z=(z1,z2)和r为估计的圆心坐标和半径。
经过静态杂波去除后,利用反正切解调提取目标相位信息
其中Pπ为π 的倍数,代表相位解缠绕过程,目的是消除信号轨迹越过复平面中两个相邻象限边界时发生的相位不连续性。
一般而言,θ0随时间几乎不变化,可以通过减均值消除相位曲线的直流分量
但睡眠状态发生变化的情形下,胸腔起伏的中心位置有时会随时间发生缓变的轻微移动,使得θ0呈现随时间逐渐变大或减小的线性趋势,对此采用线性趋势消除算法[29]通过线性最小二乘拟合估计出呼吸相位在慢时间方向上的线性趋势并减去,即
其中,β=[β1,β2],β1=[0,1,…,N-1]Τ,β2=[1,1,…,1]Τ。
最后,考虑到人体呼吸频率的合理范围为0.1~0.8 Hz,对所提取的相位数据进行带通滤波,消除高频噪声和不规则心跳的干扰以获取尽可能干净的呼吸信息。
呼吸相位信息提取效果如图8所示,其中图(a)显示了静态杂波滤除前后的数据IQ星座图,可见该静态杂波去除方法可以正确地拟合IQ 平面中目标运动形成的圆弧轨迹,较好地去除静态杂波;图(b)显示了增强处理前后的相位曲线,可见该相位提取和增强处理方法能够较好地消除静态杂波和线性趋势的影响,完整地保留并还原了呼吸的相位信息。
图8 相位增强处理前后对比Fig.8 Comparison of results before and after phase-enhancement processing
3.3.1 基于一维信号的波形特征提取
在呼吸相位提取预处理的基础上,针对一维呼吸信号分别提取时域和频域特征两类典型特征。
(1)时域特征提取
不同呼吸模式的呼吸深度及其变化规律存在很大差异,正常呼吸的深度保持稳定,异常呼吸却往往伴随着不同程度的加深或减弱,而从时域波形中提取出来的波峰波谷信息能够反映每个呼吸周期的潮气水平和起伏强度,时域波形如图9 所示。针对时域波形提取波峰、波谷信息,通过提炼数据的统计学特性,采用平均值、最大值、最小值、方差、峰谷值和波峰数量10 个特征参数构成时域特征向量。
图9 呼吸波形的波峰波谷搜寻Fig.9 Peak and valley search results for respiratory waveforms
(2)频域特征提取
呼吸频率反映了呼吸节律,是区分呼吸模式的重要特征。正常呼吸的频率较低且单一稳定,而异常呼吸模式则往往存在着不同的节律改变。因此,可以提取呼吸信号的频域信息,以频域频带、频率变化来表征呼吸的节律,作为判断呼吸模式的依据之一。本文从频域频带分布和时频变化两方面分析频域特征。
a.全局信号的频带特征
传统方法往往是单一地提取频谱峰值作为频率特征,而频带曲线的形态反映了呼吸信号在整个探测时长内的频率分量及其分布(平均频域特征),不同呼吸模式的频谱形态变化、频率峰值的位置均有所差异。例如:正常呼吸和库式呼吸的频谱更加干净,往往只有一个单一且尖锐的峰;潮式呼吸、间停呼吸和比奥呼吸中存在呼吸暂停,则频谱会在零频附近存在较高幅值。
因此,为了更好地体现不同呼吸模式的频率成分,如图10所示,本文基于频带曲线提取1/4频带面积和3/4 频带面积处所对应的频率值分别作为频带的频率下限flower和频率上限fupper,并提取其对应的频带宽度fband,以这3 个特征共同构成频带分布特征,其中
图10 某潮式呼吸频谱Fig.10 Τidal breathing spectrum
b.局部信号的时频特征
呼吸障碍相较于正常呼吸存在显著的呼吸频率变化,而呼吸的雷达回波信号属于典型的非平稳信号,傅里叶分析得到的频谱是全局信号分析,无法表现出信号局部的频率。因此,采用时频分析方法获取呼吸信号的时频图像,可进一步提取描述呼吸频率变化特性的时变特征。
短时傅里叶变换(Short-Τime Fourier Τransform,SΤFΤ)是常用的一种时频分析方法,其本质是选取一定长度的窗函数,移动窗函数来截取局部信号进行傅里叶变换[30],表达式为
其中,y(t)为呼吸雷达回波,w(t)为窗函数。对呼吸信号做SΤFΤ 分析,则可以获得呼吸的频率变化特征。
图11展示了不同呼吸模式下的信号时频图,可以看出不同呼吸模式的时频图中频率随时间的变化规律存在显著差异。本文利用重心法提取每段局部信号的重心频率用以估计该短时间内的瞬时呼吸频率,从中提取频率时变特征。对于局部呼吸信号作SΤFΤ 计算得到的离散功率谱,重心频率的计算公式可表示为
图11 不同呼吸模式的时频图Fig.11 Τime-frequency plots of different respiratory patterns
其中,S(f)为功率谱,f1~f2为频率范围。对于得到的重心频率数据G[fg1,fg2,…,fgk],其中k为窗函数移动获取的局部呼吸频率的个数,求取其均值、方差、最大值、最小值、范围5 个特征作为时频特征,同上述3个频带分布特征共同构成频率特征向量。
3.3.2 基于二维图像的形态特征提取
呼吸信号的雷达回波时距图像与单行呼吸信号波形类似,也呈周期性变化,但受目标呼吸节律影响,图像往往呈现包络放缩和局部粘连特性。图12为分割出来的呼吸子图,明暗区域表示人体呼吸的起伏。文献[25]提出通过边缘检测、膨胀腐蚀的图像处理来提取呼吸图像中的“斑点”用以显示呼吸周期,但这往往只适用于图12(a)所示的“斑点”轮廓清晰的呼吸图像中,而当人体呼吸较为微弱或者起伏较大时,“斑点”轮廓变得模糊且发生粘连,如图12(b)所示,直接对图像进行图像处理效果不佳。
图12 不同情形下的呼吸图像Fig.12 Respiratory images in different situations
为此,本文提出了结合相位矩阵图像处理提取子图时距呼吸条带来获取呼吸周期的方法。具体步骤如下:
步骤1将由呼吸回波子矩阵提取出来的相位矩阵转化为归一化灰度图。
步骤2对灰度矩阵进行图像的滑窗阈值分割处理。考虑不同呼吸周期中强度不同导致不同区域图像灰度相差很大,这里采用基于滑窗的局部阈值方法,对相位灰度矩阵选取适宜的滑窗长度,沿慢时间方向滑动滑窗,在每个滑窗区域内对灰度矩阵做二值化处理得到二维检测条带图。
步骤3用二值化矩阵提取呼吸子图中的周期性呼吸时距条带,用于后期特征提取。整个流程的处理结果如图13所示。
图13 呼吸时距条带提取Fig.13 Extraction of breath time-interval bands
对处理得到的呼吸时距条带分别提取以下图像特征:
a.条带宽度:宽度特征表征着呼吸的缓急程度,宽度越小代表呼吸越急促,宽度越大代表呼吸越平缓。
b.条带间距:时距呼吸条带的间距也包含着呼吸的快慢信息,间距越小代表呼吸越快,间距越大代表呼吸越慢。
c.条带能量:能量特征表征着呼吸的强度,尤其是异常呼吸模式下的时距呼吸条带能量具有相应的变化规律。条带内能量越大代表当前呼吸强度越大,引起的胸腔起伏越大,能量越小代表呼吸越微弱,引起的胸腔起伏越小。
d.条带峰宽比:不同呼吸状态下时距呼吸条带的能量空间形状不同,用峰宽比表示,即能量峰值高度与呼吸条带的宽度之比。峰宽比越大表示呼吸条带空间形状越尖锐,呼吸越剧烈。图14所示为时距呼吸条带峰宽比示意图。
图14 呼吸时距条带峰宽比示意图Fig.14 Schematic diagram of peak-to-width ratio of respiratory time-interval bands
求取以上特征的最大值、最小值、均值和方差共16个呼吸图像的形态特征参数,共同构成图像特征向量。
实验使用课题组自研的KUN-DC-1000 脉冲超宽带雷达进行数据采集,雷达系统参数如表1所示。实验邀请了6 名体型不同、身体健康、年龄在22~25 岁之间的志愿者进行呼吸样本采集,包括2 名女生和4 名男生。采集设置在空旷的走廊、有少数桌椅摆放的会议室和较多杂物堆放的实验室这三个场景内进行,以模拟各类真实环境。志愿者卧躺在小床上,根据正常呼吸、潮式呼吸、间停呼吸、叹气样呼吸、比奥呼吸和库式呼吸六种呼吸模式的特点对每种呼吸进行模仿,单个测试时间为30 s,每个模式重复15 次,总共得到540 个呼吸数据,图15 为数据采集过程。对采集到的雷达回波数据进行信号预处理和相应的呼吸特征提取,并添加对应的类别标签1~6,完成特征数据集的构建。
表1 IR-UWB雷达系统参数Tab.1 IR-UWB radar system parameters
图15 数据采集场景Fig.15 Data collection scenarios
为验证本文提出的二维时距图像特征及其与一维时域、频域特征构成的融合特征在IR-UWB 雷达呼吸模式识别问题上的有效性,本文对传统一维特征、二维时距图像特征及多域融合特征的分类性能进行了对比,并分析引入所提图像特征后对各种呼吸模式的分类效果。分类性能通过分类准确率AC、精确率PPV 和召回率ΤPR 三种评价指标衡量,评价指标定义如下:
其中,ΤP表示正确预测的正样本数,ΤN表示正确预测的负样本数,FP表示错误预测的正样本数,FN 表示错误预测的负样本数。实验采用5折交叉验证法来对分类模型进行训练和验证,即将所有数据分成5 部分,其中4/5 的数据用于训练,1/5 的数据用于测试,重复此过程5 次,用5 次交叉训练验证得到的准确率平均值对分类模型性能进行评价。
首先,针对决策树(DΤ)、线性判别(LDA)、K-最近邻法(KNN)、支持向量机(SVM)四种现有研究工作中常见的机器学习模型对六种呼吸模式的不同特征数据进行分类识别,几种典型分类器的分类准确率如表2所示。
表2 不同分类器对不同特征集的分类准确率Tab.2 Classification accuracies of different classifiers for different feature sets
从表2中可以看到,对于不同的机器学习模型,二维图像特征与一维时、频域特征对六种呼吸模式的分类效果相当,体现了所提二维图像特征相较于传统一维时、频域特征的有效性,而将图像特征与时、频域特征融合,可以获得更加高于一维特征的分类准确率。总体上看,采用图像特征及融合特征的LDA、KNN、SVM 三种分类器都达到了90%以上的准确率,本文所提特征及特征融合处理表现出较好的呼吸模式区分度。其中SVM 分类器显著优于其他分类器,可以实现96.3%的最佳准确率。
为评估时距图像特征及融合特征在不同呼吸模式识别中的有效性,分析不同特征下采用SVM分类器对六种呼吸模式的分类结果,其准确率和性能对比曲线分别如表3 和图16 所示。根据实验结果,单独使用时域特征或频域特征时呼吸模式的识别效果欠佳;当同时使用时域特征和频域特征进行分类时,识别率提高到了93.1%;使用本文提出的图像特征识别准确率同样大幅提高到了93.5%,分类效果明显;而利用时域、频域、图像的融合特征进行呼吸模式分类,识别率可进一步提高到96.3%,各类呼吸模式的ΤPR、PPV 指数都有所提高,绝大部分呼吸模式能正确识别。可见,本文提出的引入时距信息的IR-UWB 雷达多域特征提取和呼吸模式识别方法能够弥补单个距离门上提取的呼吸信号特征的不足,更好地表征呼吸,提高呼吸模式的识别准确率。
表3 不同特征集的SVM分类准确率Tab.3 SVM classification accuracies for different feature sets
图16 不同特征集分类性能对比Fig.16 Comparison of classification performances with different feature sets
图17显示了采用SVM 分类器对六种呼吸模式分类的混淆矩阵。如图所示,潮式呼吸、间停呼吸和叹气样呼吸基本都能正确识别,比奥呼吸的少部分样本数据易被识别为潮式呼吸和间停呼吸;而正常呼吸和库式呼吸由于特征相近容易互相混淆,这是由于个体差异因素导致不同被测对象的正常呼吸和库式呼吸存在临界情形下频率、幅度等相近的问题。
图17 融合特征的SVM分类混淆矩阵Fig.17 SVM classification confusion matrix based on fusion features
针对不同情形的呼吸模式分类,本文提出了引入时距信息的IR-UWB 雷达多域特征融合呼吸模式识别方法。区别于传统方法仅在单个距离门的一维呼吸信号上提取特征,本文针对IR-UWB 雷达目标回波信息距离扩展的特性,在提取了一维呼吸信号时域、频域特征的基础上,还选取了多个距离门的呼吸数据,从时距图像角度挖掘了呼吸模式在雷达回波图像上呈现出来的形态特征。最后在实验中验证了融合多域特征识别的有效性,并基于SVM 分类器实现了正常呼吸和5 种呼吸障碍的识别,达到了96.3%的准确率。实验结果表明本方法在呼吸模式的辨识上相较于传统方法具有更优的识别性能。