具体实施方式
下面结合具体的实施例及附图对呼吸信号检测方法和装置的技术方案进行详细的描述,以使其更加清楚。
如图1所示,在一个实施例中,一种呼吸信号检测方法,包括以下步骤:
步骤S110,采集原始信号。
具体的,可通过雷达探测仪在灾后救援现场或反恐场景中采集原始信号,该原始信号中可能包括呼吸信号。该原始信号组成雷达回波矩阵,该雷达回波矩阵中通常包含大量噪声,杂波和干扰。雷达回波矩阵可由式(1)表示:
式(1)中,表示雷达回波矩阵;r[m,n]表示由于胸腔表面运动产生的信号;c[n]表示由背景环境产生的杂波;w[m,n]表示系统噪声,通常符合高斯白噪声模型;d[m]表示由采样引入的不稳定的直流基线;l[m,n]表示由于采样幅值不稳定性引入的线性趋势,在慢时间域上通常显示为基线漂移。
步骤S120,对该原始信号进行滤波处理,且该滤波处理中采用的滤波参数为根据预先分析得到的慢时间域呼吸信号的谐波结构确定的。
对原始信号进行滤波处理的目的是为了抑制c[n]、w[m,n]、d[m]和l[m,n],保留r[m,n]。
在一个实施例中,在步骤S110之前,还包括步骤:预先分析得到慢时间域呼吸信号的谐波结构。
如图2所示,在一个实施例中,预先分析得到慢时间域呼吸信号的谐波结构包括如下步骤:
步骤S210,获取并根据脉冲回波延迟时间、慢时间域采样周期和快时间域采样周期构建离散时间回波矩阵的仿真模型,并设置该离散时间回波矩阵的仿真参数。
由于UWB脉冲雷达回波数据存在两个维度,即快时间维度和慢时间维度,因而通常被称为回波矩阵。UWB雷达对人体呼吸运动的检测主要是对人体胸腔周期起伏运动的检测。当UWB雷达对人体探测时,人体胸腔的起伏运动会导致胸腔表面与雷达之间的径向距离的周期性变化。这一变化导致回波脉冲的延迟时间的周期波动。当人体吸气时,胸腔膨胀,胸腔表面与雷达间的距离变短,经胸腔表面反射的脉冲回波的延迟时间缩短;反之,当人体呼气,胸腔收缩,胸腔表面与雷达间的距离变长,经胸腔表面反射的脉冲回波的延迟时间延长。回波脉冲延迟时间的周期波动反应在二维回波数据(快时间域和慢时间域)上是在慢时间域上信号具有周期性。为此,呼吸信号是指因胸腔运动产生回波脉冲延迟时间的周期波动,再由该周期波动导致二维回波数据在慢时间域上形成的周期信号。其中,快时间域是指单个脉冲发射及回波时间形成的时间域,慢时间域是指每个脉冲发射时间点形成的时间域。
本实施例中,采用高阶余弦信号对人体呼吸导致胸腔表面的起伏运动进行建模。设被测目标与雷达间的距离为d0,人体胸腔起伏运动的幅度为B,人体的呼吸率为fr,则人体胸腔表面与雷达间的距离d(t)随时间t的变化可用式(2)表示:
d(t)=d0-B×(cosπfrt)u (2)
式(2)中,u表示高阶余弦信号的阶数。通常u取6时,(2)式可以较好的近似人体胸腔表面与雷达间距离的变化过程。图3给出了当d0=4m,B=0.01m,fr=0.25Hz,u=6的情况下,根据该模型给出的胸腔表面与雷达间距离随时间变化的示意图,其中,m为米,Hz为赫兹。图3中横坐标为时间,单位为秒,纵坐标为人体胸腔表面与雷达间的距离,单位为米。
根据d(t)计算反射的脉冲回波的延迟时间td(t),如式(3)所示:
式(3)中,c表示光速。
采用p(t)表示雷达发射机产生的脉冲波形。设雷达慢时间域采样周期为Tst,雷达快时间域采样周期为Tft。根据脉冲回波延迟时间td(t)、慢时间域采样周期Tst和快时间域采样周期Tft构建离散时间回波矩阵的仿真模型,可得到离散时间回波矩阵的仿真模块如式(4)所示:
r[m,n]=p(nTft-td(mTst)) (4)
本实施案例中,雷达发射机采用的是一阶高斯脉冲,脉冲波形p(t)表达式如下:
式(5)中,A为脉冲幅度因子,σ为时间因子。因脉冲幅度因子的取值并不影响后续分析,故在此不作特殊设置。时间因子σ在本实施案例中取值128ps,ps为皮秒,即10-12秒。把式(5)代入式(4)可得在本实施案例下离散时间回波矩阵表达式:
其中,Tst表示慢时间域采样周期,Tft表示雷达快时间域采样周期。
设置该离散时间回波矩阵的仿真参数,Tst和Tft的取值应该确保采样在慢时间域和快时间域上均满足采样定理。为了确保离散时间回波矩阵在慢时间域上是周期信号,在仿真参数设置时应确保为fr的整数倍。通常,为fr的100倍。该值较为合理,一方面既保证为fr的整数倍,另一方面,较大,可确保在慢时间域上采样满足采样定理。Tft的取值可通过脉冲波形p(t)计算脉冲频谱,根据频谱确定脉冲的频带宽度,根据频带宽度确定Tft,保证在快时间域上采样满足采样定理。
步骤S220,根据该仿真模型和仿真参数获取离散时间回波矩阵的仿真结果。
具体的,采用matlab软件作为仿真平台。根据仿真模型,设置仿真参数,经matlab进行仿真,以获得图4所示的回波矩阵,图4中,d0=4m,B=10mm,fr=0.25Hz,A=10-9,σ=128ps,Tst=0.04s,Tft=10ps,其中,m为米,mm为毫米,Hz为赫兹,ps为皮秒,横坐标为快时间域变量n,纵坐标为慢时间域变量m,且右边为与左边对应的灰度刻度条(即图中的-6至6)。
步骤S230,根据该仿真结果提取慢时间域呼吸信号。
具体的,对于回波矩阵r[m,n],rm[n]用于表示单路快时间域信号,其中,快时间域离散时间变量n可看作该函数的自变量,慢时间域离散时间变量m可看作固定参数;类似的,rn[m]用于表示单路慢时间域信号,其中,慢时间域离散时间变量m可看作该函数的自变量,快时间域离散时间变量n可看作固定参数。
在雷达信号处理中,为了抑制杂波,任何信号的直流分量都会被移除,因此,对各路rn[m]取均值,获得不含直流分量的慢时间域信号计算均值的方法是:由于rn[m]的周期已知(在本实施案例中为100),首先仿真出一个周期的rn[m],然后对一个周期的仿真数据求和再除以周期即得到rn[m]的均值。典型的慢时间域呼吸信号指的是在回波矩阵中具有最大功率的单路慢时间域呼吸信号其中,各路慢时间域信号功率可通过对信号平方再取平均即可获得。
在本实施案例中,按照上述方法进行计算,得到本案例典型的慢时间域呼吸信号为n=2663。
步骤S240,利用离散傅立叶级数分析该慢时间域呼吸信号,得到该慢时间域呼吸信号的谐波结构向量。
首先定义呼吸信号的谐波结构为一向量,该向量的第n个元素对应第n次谐波,该向量的第n个元素的取值为第n次谐波的功率在信号总功率中的比重。由于典型的慢时间域呼吸信号是离散周期信号,可以采用离散傅立叶级数对其进行分析,获得其谐波结构。
在本实施案例中,由图4中所给参数设置知,fst=100×fr,fst为慢时间域采样频率,且这意味着,这些呼吸信号经过慢时间域采样后都是周期的,且离散周期TM=100。
对于离散周期信号,采用离散傅立叶级数进行分析。离散傅立叶级数的计算公式为:
k=0,1,...,TM-1 (6)
在本实施案例中,对做离散傅立叶级数的计算,结果如图5所示,图5中离散傅立叶级数k=0,1,...,9。
由于直流分量已被移除, 和对应第k次谐波,且基于此,式(7)给出了计算第k次谐波的功率在总功率中所占比重Ratio(k)的计算方法:
根据图5给出的离散傅立叶级数计算呼吸信号的谐波结构向量。图6用画图的方式显示了典型的慢时间域呼吸信号的谐波结构向量,横坐标为谐波阶次,纵坐标为功率比重。从图6可以看出,在本实施案例中,1次谐波和2次谐波占据了信号的绝大部分功率,因此,在谐波结构向量中只保留1次和2次谐波,得到典型的慢时间域呼吸信号的谐波结构向量:
上述通过离散傅立叶级数分析得到了呼吸信号的谐波结构,方便结合呼吸信号的谐波结构设置相应的滤波参数,以提高滤波的质量,避免滤除存在的呼吸信号。
进一步的,在一个实施例中,步骤S120包括以下步骤:
(1)对原始信号中每一路信号进行快时间域滤波处理得到相应路的快时间域滤波输出信号,并组成快时间域滤波输出矩阵。
具体的,快时间域滤波处理采用快时间域滤波器进行滤波处理。该快时间域滤波器为带通线性相位有限脉冲响应数字滤波器。原始信号为包含噪声的回波矩阵,回波矩阵的每一路快时间信号经过该带通线性相位有限脉冲响应数字滤波器获得相应的一路输出信号,多路输出信号组成快时间域滤波输出矩阵。该带通线性相位有限脉冲响应数字滤波器的频率响应的带宽应与雷达脉冲频带相匹配。例如,雷达脉冲的带宽是从1GHz到3GHz,则为了确保雷达脉冲在幅度谱上不出现失真,滤波器的频率响应在1GHz到3GHz上应尽可能保持近似为1,而在其它频段,为了抑制噪声,应尽可能接近为0。此外,快时间域滤波处理采用带通线性相位有限脉冲响应数字滤波器是为了保持回波脉冲的时域波形。
快时间域滤波输出矩阵rftF[m,n]的表达式如下:
式(8)中,为雷达回波矩阵(即原始信号组成的回波矩阵),hftF[n]为带通线性相位有限脉冲响应数字滤波器的单位脉冲响应;*n表示快时间域的时域卷积计算。
在本实施案例中,雷达脉冲的带宽为0.45GHz到3.555GHz,而由于使用的天线的带宽为0.9GHz到5GHz,综合考虑,回波脉冲的带宽为0.9GHz到3.555GHz。采用经典的线性相位有限脉冲响应数字滤波器设计技术,得到本实施案例中所使用的快时间域滤波参数如图7A和图7B所示,图7A表示快时间域滤波器的幅度谱响应,图7A中横坐标为频率,纵坐标为幅值;图7B表示快时间域滤波器的冲击响应,图7B中横坐标为时间,纵坐标为幅值。该快时间域滤波器的阶数为83,群延时为42。
通过该快时间域滤波,回波矩阵r中的d[m]能被有效抑制。
在本实施案例中,通过实验,采集了一套实验数据,如图8所示,为存在呼吸信号的雷达回波实验数据示意图,且右边为与左边对应的灰度刻度条(即图8中的50.3至51)。在图8中,各点的数值都在50.3至51.1之间变化,d[m]>50,其值远远大于其它信号的数值,为了抑制d[m],对实验数据在快时间域上滤波,输出如图9所示。图9为快时间域滤波后的输出矩阵图,且右边为与左边对应的灰度刻度条(图9中的-0.15至0.2)。从图9中可以看到,各点数值在-0.2至0.2之间变化,d[m]得到有效抑制。
(2)对该快时间域滤波输出矩形的每一路快时间域滤波输出信号进行带通慢时间域滤波处理得到相应路的带通慢时间域滤波输出信号,并组成慢时间域滤波输出矩阵,其中,该带通慢时间域滤波处理中采用的滤波参数为根据该慢时间域呼吸信号的谐波结构确定的。
具体的,慢时间域滤波处理包括带通慢时间域滤波处理和高通慢时间域滤波处理,慢时间域滤波处理采用慢时间域滤波器进行,该慢时间域滤波器由2个线性相位有限脉冲响应数字滤波器组成,即带通慢时间域滤波器和高通慢时间域滤波器。快时间域滤波输出矩阵的每一路快时间域滤波输出信号经过慢时间域滤波器获得相应的输出信号。
带通慢时间域滤波器作用是保留回波矩阵中的r[m,n],抑制c[n]、w[m,n]和l[m,n]。根据呼吸信号的谐波结构确定带通慢时间域滤波器的参数,即根据呼吸信号的谐波结构确定带通慢时间域滤波处理中采用的滤波参数。在本实施案例中,在呼吸信号的谐波结构分析中,已知呼吸信号主要包含1次谐波和2次谐波。基于此,考虑到人体呼吸运动的频率主要在0.2Hz到0.5Hz之间浮动,呼吸信号的频带主要为0.2Hz至1Hz。保留r[m,n]对于慢时间域滤波来说就是保留呼吸信号,因而,带通慢时间域滤波处理中的滤波参数包括滤波频率,该滤波频率为0.2赫兹到1赫兹。
为此,带通慢时间域滤波器的频率响应在0.2Hz到1Hz上应尽可能保持近似为1,而在其它频段,为了抑制噪声,应尽可能接近为0。这里需要指出,c[n]和l[m,n]在慢时间域上通常集中在极低频段,因而通过该慢时间域滤波后均会被有效抑制。
此外,慢时间域滤波处理采用带通线性相位有限脉冲响应数字滤波器是为了保持呼吸信号的时域波形。
带通慢时间域滤波输出矩阵rstF[m,n]的计算如下式所示:
rstF[m,n]=rftF[m,n]*mhstF[m] (9)
式(9)中,hstF[m]为带通慢时间域滤波器的单位脉冲响应;*m表示慢时间域的时域卷积计算。
本实施案例中所使用的带通慢时间域滤波器参数如图10A和图10B所示。图10A表示带通慢时间滤波器的幅度谱,图10A中横坐标为频率,纵坐标为幅值;图10B表示带通慢时间域滤波器的冲击响应,图10B中横坐标为时间,纵坐标为幅值。该带通慢时间域滤波器的阶数为401,群延时为200。
本实施案例中,经带通慢时间域滤波器处理后的回波矩阵(即实验输出数据)如图11所示,且右边为与左边对应的灰度刻度条(即图11中的-0.015至0.015)。
(3)对该快时间域滤波输出矩阵的每一路快时间域滤波输出信号进行高通慢时间域滤波处理得到相应路的高通慢时间域滤波输出信号,并组成高通慢时间域滤波输出矩阵。
高通慢时间域滤波器主要是为了获取一个只包含w[m,n],用于接下来估计噪声功率。高通慢时间域滤波器的性能要求就是下截止频率足够高,能够有效抑制r[m,n]、c[n]和l[m,n]。因此根据前述分析,在本实施案例中,高通慢时间域滤波器的下截止频率至少应大于1Hz。
高通慢时间域滤波输出矩阵的计算如下式所示:
式(10)中,为高通慢时间域滤波器的单位脉冲响应;*m表示慢时间域的时域卷积计算。
本实施案例中所使用的带通慢时间域滤波器参数如图12A和图12B所示。图12A表示高通慢时间滤波器的幅度谱,图12A中横坐标为频率,纵坐标为幅值;图12B表示高通慢时间域滤波器的冲击响应,图12B中横坐标为时间,纵坐标为幅值。该高通慢时间域滤波器的阶数为201,群延时为100。
本实施案例中,经高通慢时间域滤波器处理后的回波矩阵(即实验输出数据)如图13所示,且右边为与左边对应的灰度刻度条(即图13中的-1.5至1.5)。
步骤S130,计算滤波处理后的原始信号的慢时间域信号的周期图。
具体的,步骤130具体包括:根据带通慢时间域滤波输出矩阵的每一路带通慢时间域滤波输出信号进行周期图计算,得到各路慢时间域信号的周期图。
在本实施案例中,基于带通慢时间域滤波器实验输出矩阵,计算各路慢时间域信号的周期图In[k],计算公式如下:
k=0,1,2,...,M-1 (11)
式(11)中,M为带通慢时间域滤波器输出矩阵的单路慢时间域信号的数据长度。在本实施案例中M=2000,为慢时间域信号的快速傅立叶变换(FFT)。图14中给出了带通慢时间域滤波器输出矩阵的周期图,且右边为与左边对应的灰度刻度条(即图14中的0.005至0.02)。图14中,在n=175附近,周期图In[k]存在较显著的一次谐波和二次谐波信号。
步骤S140,估算滤波处理后的原始信号的慢时间域信号的噪声功率。
具体的,步骤S140具体包括:根据该高通慢时间域滤波输出矩阵的每一路高通慢时间域滤波输出信号估算各路慢时间域信号的噪声功率。
根据高通慢时间域滤波输出矩阵的每一路慢时间域信号估计各路慢时间域信号的噪声功率,计算公式如下:
式(12)中,Pn表示第n路慢时间域信号的噪声功率估计值,为高通慢时间域滤波输出矩阵,为高通慢时间域滤波器的单位脉冲响应。
步骤S150,根据该慢时间域呼吸信号的谐波结构、慢时间域信号的周期图及慢时间域信号的噪声功率计算谐波图。
具体的,谐波图Hn[k]计算公式如下:
其中,Jh表示周期信号所含谐波阶次的最大值;aj等于谐波结构向量的第j个分量的取值,In[k]为各路慢时间域信号的周期图。
在本实施案例中,由谐波结构分析知,Jh=2,a1=0.8562,a2=0.1404。
根据谐波图定义,由实验数据计算所得谐波图如图15所示,横坐标为快时间域,纵坐标为频率,且右边为与左边对应的灰度刻度条(即图15中的5至40)。
步骤S160,判断该谐波图是否大于预设门限,若是,则执行步骤S170,若否,则执行步骤S180。
基于谐波图,由呼吸信号的二元假设判决,判决公式如下:
式(14)中,H1表示存在呼吸信号,H0表示不存在呼吸信号。当Hn[k]大于预设门限r时,H1为判决结果,反之,Hn[k]小于等于预设门限r时,H0为判决结果。当预设门限r过高,会降低检测算法的检测概率,而预设门限r过低,又会增大检测算法的虚警率。
在图15所示的本实施案例的谐波图中,通过适当的选取门限值,例如r=10,将能成功检测到呼吸信号。
步骤S170,表示存在呼吸信号。
步骤S180,表示不存在呼吸信号。
上述呼吸信号检测方法,采用慢时间域呼吸信号的谐波结构确定滤波参数,根据确定后的滤波参数滤波器对采集的原始信号进行滤波处理,若原始信号中存在呼吸信号,则通过此滤波处理可较好的滤除环境噪声且尽可能的保留呼吸信号,根据呼吸信号谐波结构及计算得到的慢时间域信号周期图和噪声功率计算谐波图,通过谐波图与预设门限进行判断,提高了检测的性能,更加有效的检测呼吸信号是否存在。
进一步的,在一个实施例中,当存在呼吸信号时,上述呼吸信号检测方法还包括:根据该谐波图的最大值估算目标距离和/或目标呼吸率。
具体的,在实际系统中,由于各种电路元器件和数字滤波器的非理想特性,脉冲往往会被展宽,同时,考虑到在应用场景中不可避免的多经效应,单个人体的呼吸运动在回波矩阵的快时间域中会被展宽和延伸。这一点在图15所示的谐波图中很好的得到了验证。但是,考虑到人体胸腔表面的直接反射信号往往强于其它多径信号,当判决存在呼吸信号时,对目标距离和呼吸率进行估算的步骤包括:
(1)获取Hn[k]的最大值,其对应的快时间域变量n用表示,其对应的离散时间频率用表示。
(2)将作为呼吸率的估计值。
(3)估计目标距离d。
具体的,按传统的雷达脉冲信号计算距离的方式结合快时间域滤波器对距离估计的影响可得,目标距离计算公式如下:
式(15)中,delayftF为快时间域滤波器的群延时,c表示光速,Tft为雷达快时间域采样周期,表示快时间域变量n。目标距离是指被测目标与检测装置之间的距离。
在本实施案例中,图15中小箭头“”所指向的Hn[k]即为谐波图的最大值,根据该最大值即可按式(15)估计呼吸率和目标距离。
上述呼吸信号检测方法对呼吸信号的检测性能明显优于传统的UWB脉冲探测方法(即将呼吸信号视作单一频点的正弦信号的检测方法)对呼吸信号检测。如图16所示为两种检测方法的检测性能对比示意图,横坐标为信噪比(单位dB),纵坐标为检测概率。图16所示,在高斯白噪声信道模型下,通过MonteCarlo仿真,获得两种方法在虚警率Pfa为0.01和0.001的情况下对于不同的信噪比(SNR)所能达到的检测概率的对比图,图中黑色圆点(即图16中“●”)连线为传统UWB探测方法在Pfa为0.01时不同信噪比所达到的检测概率,黑色矩形点(即图16中“■”)连线为本发明中呼吸信号检测方法在Pfa为0.01时不同信噪比所达到的检测概率,无色圆点(即图16中“○”)连线为传统UWB探测方法在Pfa为0.001时不同信噪比所达到的检测概率,无色矩形点(即图16中“□”)连线为本发明中呼吸信号检测方法在Pfa为0.001时不同信噪比所达到的检测概率,从图16中可以看出本发明中呼吸信号检测方法的检测性能明显优于传统的UWB探测方法。
如图17所示,在一个实施例中,一种呼吸信号检测装置,包括采集模块110、噪声抑制模块120、周期图计算模块130、噪声功率估算模块140、谐波图计算模块150和判断模块160。其中:
采集模块110用于采集原始信号。
具体的,采样模块110为雷达探测仪或其他信号检测装置,本实施例中,通过雷达探测仪在灾后救援现场或反恐场景中采集原始信号,该原始信号中可能包括呼吸信号。该原始信号组成雷达回波矩阵,该雷达回波矩阵中通常包含大量噪声,杂波和干扰。雷达回波矩阵可由式(1)表示:
式(1)中,表示雷达回波矩阵;r[m,n]表示由于胸腔表面运动产生的信号;c[n]表示由背景环境产生的杂波;w[m,n]表示系统噪声,通常符合高斯白噪声模型;d[m]表示由采样引入的不稳定的直流基线;l[m,n]表示由于采样幅值不稳定性引入的线性趋势,在慢时间域上通常显示为基线漂移。
噪声抑制模块120,用于对该原始信号进行滤波处理,且该滤波处理中采用的滤波参数为根据预先分析得到的慢时间域呼吸信号的谐波结构确定的。
对原始信号进行滤波处理的目的是为了抑制c[n]、w[m,n]、d[m]和l[m,n],保留r[m,n]。
进一步的,上述呼吸信号检测装置还包括呼吸信号谐波结构分析模块170,用于预先分析得到慢时间域呼吸信号的谐波结构。
如图17中,呼吸信号谐波结构分析模块170包括建模单元171、仿真单元173、提取单元175和分析单元177。其中:
建模单元171,用于获取并根据脉冲回波延迟时间、慢时间域采样周期和快时间域采样周期构建离散时间回波矩阵的仿真模型,并设置该离散时间回波矩阵的仿真参数。
本实施例中,采用高阶余弦信号对人体呼吸导致胸腔表面的起伏运动进行建模。设被测目标与雷达间的距离为d0,人体胸腔起伏运动的幅度为B,人体的呼吸率为fr,则人体胸腔表面与雷达间的距离d(t)随时间t的变化可用式(2)表示:
d(t)=d0-B×(cosπfrt)u (2)
式(2)中,u表示高阶余弦信号的阶数。通常u取6时,(2)式可以较好的近似人体胸腔表面与雷达间距离的变化过程。图3给出了当d0=4m,B=0.01m,fr=0.25Hz,u=6的情况下,根据该模型给出的胸腔表面与雷达间距离随时间变化的示意图,其中,m为米,Hz为赫兹。图3中横坐标为时间,单位为秒,纵坐标为人体胸腔表面与雷达间的距离,单位为米。
根据d(t)计算反射的脉冲回波的延迟时间td(t),如式(3)所示:
式(3)中,c表示光速。
采用p(t)表示雷达发射机产生的脉冲波形。设雷达慢时间域采样周期为Tst,雷达快时间域采样周期为Tft。根据脉冲回波延迟时间td(t)、慢时间域采样周期Tst和快时间域采样周期Tft构建离散时间回波矩阵的仿真模型,可得到离散时间回波矩阵的仿真模块如式(4)所示:
r[m,n]=p(nTft-td(mTst)) (4)
本实施案例中,雷达发射机采用的是一阶高斯脉冲,脉冲波形p(t)表达式如下:
式(5)中,A为脉冲幅度因子,σ为时间因子。因脉冲幅度因子的取值并不影响后续分析,故在此不作特殊设置。时间因子σ在本实施案例中取值128ps,ps为皮秒,即10-12秒。把式(5)代入式(4)可得在本实施案例下离散时间回波矩阵表达式:
其中,Tst表示慢时间域采样周期,Tft表示雷达快时间域采样周期。
设置该离散时间回波矩阵的仿真参数,Tst和Tft的取值应该确保采样在慢时间域和快时间域上均满足采样定理。为了确保离散时间回波矩阵在慢时间域上是周期信号,在仿真参数设置时应确保为fr的整数倍。通常,为fr的100倍。该值较为合理,一方面既保证为fr的整数倍,另一方面,较大,可确保在慢时间域上采样满足采样定理。Tft的取值可通过脉冲波形p(t)计算脉冲频谱,根据频谱确定脉冲的频带宽度,根据频带宽度确定Tft,保证在快时间域上采样满足采样定理。
仿真单元173,用于根据该仿真模型和仿真参数获取离散时间回波矩阵的仿真结果。
具体的,采用matlab软件作为仿真平台。根据仿真模型,设置仿真参数,经matlab进行仿真,以获得图4所示的回波矩阵,图4中,d0=4m,B=10mm,fr=0.25Hz,A=10-9,σ=128ps,Tst=0.04s,Tft=10ps,其中,m为米,mm为毫米,Hz为赫兹,ps为皮秒,横坐标为快时间域变量n,纵坐标为慢时间域变量m,且右边为与左边对应的灰度刻度条(即图中的-6至6)。
提取单元175,用于根据该仿真结果提取慢时间域呼吸信号。
具体的,对于回波矩阵r[m,n],rm[n]用于表示单路快时间域信号,其中,快时间域离散时间变量n可看作该函数的自变量,慢时间域离散时间变量m可看作固定参数;类似的,rn[m]用于表示单路慢时间域信号,其中,慢时间域离散时间变量m可看作该函数的自变量,快时间域离散时间变量n可看作固定参数。
为了抑制杂波,任何信号的直流分量都会被移除,因此,对各路rn[m]取均值,获得不含直流分量的慢时间域信号计算均值的方法是:由于rn[m]的周期已知(在本实施案例中为100),首先仿真出一个周期的rn[m],然后对一个周期的仿真数据求和再除以周期即得到rn[m]的均值。典型的慢时间域呼吸信号指的是在回波矩阵中具有最大功率的单路慢时间域呼吸信号其中,各路慢时间域信号功率可通过对信号平方再取平均即可获得。
在本实施案例中,按照上述方法进行计算,得到本案例典型的慢时间域呼吸信号为n=2663。
分析单元177,用于利用离散傅立叶级数分析该慢时间域呼吸信号,得到该慢时间域呼吸信号的谐波结构向量。
首先定义呼吸信号的谐波结构为一向量,该向量的第n个元素对应第n次谐波,该向量的第n个元素的取值为第n次谐波的功率在信号总功率中的比重。由于典型的慢时间域呼吸信号是离散周期信号,可以采用离散傅立叶级数对其进行分析,获得其谐波结构。
在本实施案例中,由图4中所给参数设置知,fst=100×fr,fst为慢时间域采样频率,且这意味着,这些呼吸信号经过慢时间域采样后都是周期的,且离散周期TM=100。
对于离散周期信号,采用离散傅立叶级数进行分析。离散傅立叶级数的计算公式为:
k=0,1,...,TM-1 (6)
在本实施案例中,对做离散傅立叶级数的计算,结果如图5所示,图5中离散傅立叶级数k=0,1,...,9。
由于直流分量已被移除, 和对应第k次谐波,且基于此,式(7)给出了计算第k次谐波的功率在总功率中所占比重Ratio(k)的计算方法:
根据图5给出的离散傅立叶级数计算呼吸信号的谐波结构向量。图6用画图的方式显示了典型的慢时间域呼吸信号的谐波结构向量,横坐标为谐波阶次,纵坐标为功率比重。从图6可以看出,在本实施案例中,1次谐波和2次谐波占据了信号的绝大部分功率,因此,在谐波结构向量中只保留1次和2次谐波,得到典型的慢时间域呼吸信号的谐波结构向量:
上述通过离散傅立叶级数分析得到了呼吸信号的谐波结构,方便结合呼吸信号的谐波结构设置相应的滤波参数,以提高滤波的质量,避免滤除存在的呼吸信号。
进一步的,噪声抑制模块120包括:
快时间域滤波模块121,用于对该原始信号中每一路信号进行快时间域滤波处理得到相应路的快时间域滤波输出信号,并组成快时间域滤波输出矩阵。
具体的,快时间域滤波模块121为快时间域滤波器。该快时间域滤波器为带通线性相位有限脉冲响应数字滤波器。原始信号为包含噪声的回波矩阵,回波矩阵的每一路快时间信号经过该带通线性相位有限脉冲响应数字滤波器获得相应的一路输出信号,多路输出信号组成快时间域滤波输出矩阵。该带通线性相位有限脉冲响应数字滤波器的频率响应的带宽应与雷达脉冲频带相匹配。例如,雷达脉冲的带宽是从1GHz到3GHz,则为了确保雷达脉冲在幅度谱上不出现失真,滤波器的频率响应在1GHz到3GHz上应尽可能保持近似为1,而在其它频段,为了抑制噪声,应尽可能接近为0。此外,快时间域滤波器采用带通线性相位的滤波器是为了保持回波脉冲的时域波形。
快时间域滤波输出矩阵rftF[m,n]的表达式如下:
式(8)中,为雷达回波矩阵(即原始信号组成的回波矩阵),hftF[n]为带通线性相位有限脉冲响应数字滤波器的单位脉冲响应;*n表示快时间域的时域卷积计算。
在本实施案例中,雷达脉冲的带宽为0.45GHz到3.555GHz,而由于使用的天线的带宽为0.9GHz到5GHz,综合考虑,回波脉冲的带宽为0.9GHz到3.555GHz。采用经典的线性相位有限脉冲响应数字滤波器设计技术,得到本实施案例中所使用的快时间域滤波参数如图7A和图7B所示,图7A表示快时间域滤波器的幅度谱响应,图7A中横坐标为频率,纵坐标为幅值;图7B表示快时间域滤波器的冲击响应,图7B中横坐标为时间,纵坐标为幅值。该快时间域滤波器的阶数为83,群延时为42。
通过该快时间域滤波模块121滤波后回波矩阵中的d[m]能被有效抑制。
在本实施案例中,通过实验,采集了一套实验数据,如图8所示,为存在人体呼吸运动的雷达回波实验数据示意图。在图8中,各点的数值都在50.3至51.1之间变化,d[m]>50,其值远远大于其它信号的数值,为了抑制d[m],对实验数据在快时间域上滤波,输出如图9所示。图9为快时间域滤波后的实验输出数据示意图。从图9中可以看到,各点数值在-0.2至0.2之间变化,d[m]得到有效抑制。
慢时间域滤波模块123,用于对快时间域滤波输出矩阵的每一路快时间域滤波输出信号滤波获得相应的输出信号。慢时间域滤波模块123由2个线性相位有限脉冲响应数字滤波器组成,即带通慢时间域滤波器1231和高通慢时间域滤波器1233。
带通慢时间域滤波器1231,用于对快时间域滤波输出矩形的每一路快时间域滤波输出信号进行带通慢时间域滤波处理得到相应路的带通慢时间域滤波输出信号,并组成慢时间域滤波输出矩阵,其中,该带通慢时间域滤波处理中采用的滤波参数为根据该慢时间域呼吸信号的谐波结构确定的。
带通慢时间域滤波器作用是保留回波矩阵中的r[m,n],抑制c[n]、w[m,n]和l[m,n]。根据呼吸信号的谐波结构确定带通慢时间域滤波器的参数,即根据呼吸信号的谐波结构确定带通慢时间域滤波处理中采用的滤波参数。在本实施案例中,在呼吸信号的谐波结构分析中,已知呼吸信号主要包含1次谐波和2次谐波。基于此,考虑到人体呼吸运动的频率主要在0.2Hz到0.5Hz之间浮动,呼吸信号的频带主要为0.2Hz至1Hz。保留r[m,n]对于慢时间域滤波来说就是保留呼吸信号,因而,带通慢时间域滤波处理中的滤波参数包括滤波频率,该滤波频率为0.2赫兹到1赫兹。
为此,带通慢时间域滤波器的频率响应在0.2Hz到1Hz上应尽可能保持近似为1,而在其它频段,为了抑制噪声,应尽可能接近为0。这里需要指出,c[n]和l[m,n]在慢时间域上通常集中在极低频段,因而通过该慢时间域滤波后均会被有效抑制。
此外,慢时间域滤波处理采用带通线性相位有限脉冲响应数字滤波器是为了保持呼吸信号的时域波形。
带通慢时间域滤波输出矩阵rstF[m,n]的计算如下式所示:
rstF[m,n]=rftF[m,n]*mhstF[m] (9)
式(9)中,hstF[m]为带通慢时间域滤波器的单位脉冲响应;*m表示慢时间域的时域卷积计算。
本实施案例中所使用的带通慢时间域滤波器参数如图10A和图10B所示。图10A表示带通慢时间滤波器的幅度谱,图10A中横坐标为频率,纵坐标为幅值;图10B表示带通慢时间域滤波器的冲击响应,图10B中横坐标为时间,纵坐标为幅值。该带通慢时间域滤波器的阶数为401,群延时为200。
本实施案例中,经带通慢时间域滤波器处理后的回波矩阵如图11所示。
高通慢时间域滤波器1233,用于对该快时间域滤波输出矩阵的每一路快时间域滤波输出信号进行高通慢时间域滤波处理得到相应路的高通慢时间域滤波输出信号,并组成高通慢时间域滤波输出矩阵。
高通慢时间域滤波器主要是为了获取一个只包含w[m,n],用于接下来估计噪声功率。高通慢时间域滤波器的性能要求就是下截止频率足够高,能够有效抑制r[m,n]、c[n]和l[m,n]。因此根据前述分析,在本实施案例中,高通慢时间域滤波器的下截止频率至少应大于1Hz。
高通慢时间域滤波输出矩阵的计算如下式所示:
式(10)中,为高通慢时间域滤波器的单位脉冲响应;*m表示慢时间域的时域卷积计算。
本实施案例中所使用的带通慢时间域滤波器参数如图12A和图12B所示。图12A表示高通慢时间滤波器的幅度谱,图12A中横坐标为频率,纵坐标为幅值;图12B表示高通慢时间域滤波器的冲击响应,图12B中横坐标为时间,纵坐标为幅值。该高通慢时间域滤波器的阶数为201,群延时为100。
本实施案例中,经高通慢时间域滤波器处理后的回波矩阵如图13所示。
周期图计算模块130,用于计算滤波处理后的原始信号的慢时间域信号的周期图。
具体的,周期图计算模块130还用于根据带通慢时间域滤波输出矩阵的每一路带通慢时间域滤波输出信号进行周期图计算,得到各路慢时间域信号的周期图。
在本实施案例中,基于带通慢时间域滤波器实验输出矩阵,计算各路慢时间域信号的周期图In[k],计算公式如下:
k=0,1,2,...,M-1 (11)
式(11)中,M为带通慢时间域滤波器实验输出矩阵的单路慢时间域信号的数据长度。在本实施案例中M=2000,为慢时间域信号的快速傅立叶变换(FFT)。图14中给出了带通慢时间域滤波器输出矩阵的周期图。图14中,在n=175附近,周期图In[k]存在较显著的一次谐波和二次谐波信号。
噪声功率估算模块140,用于计算滤波处理后的原始信号的慢时间域信号的噪声功率。
具体的,噪声功率估算模块140还用于根据该高通慢时间域滤波输出矩阵的每一路高通慢时间域滤波输出信号估算各路慢时间域信号的噪声功率。
根据高通慢时间域滤波输出矩阵的每一路慢时间域信号估计各路慢时间域信号的噪声功率,计算公式如下:
式(12)中,Pn表示第n路慢时间域信号的噪声功率估计值,为高通慢时间域滤波输出矩阵,为高通慢时间域滤波器的单位脉冲响应。
谐波图计算模块150,用于根据所述慢时间域信号的周期图及慢时间域信号的噪声功率计算谐波图。
具体的,谐波图Hn[k]计算公式如下:
其中,Jh表示周期信号所含谐波阶次的最大值;aj等于谐波结构向量的第j个分量的取值,In[k]为各路慢时间域信号的周期图。
在本实施案例中,由谐波结构分析知,Jh=2,a1=0.8562,a2=0.1404。
根据谐波图定义,由实验数据计算所得谐波图如图15所示。
判断模块160,用于判断所述谐波图是否大于预设门限,若是,则表示存在呼吸信号,若否,则表示不存在呼吸信号。
基于谐波图,由呼吸信号的二元假设判决,判决公式如下:
式(14)中,H1表示存在呼吸信号,H0表示不存在呼吸信号。当Hn[k]大于预设门限r时,H1为判决结果,反之,Hn[k]小于等于预设门限r时,H0为判决结果。当预设门限r过高,会降低检测算法的检测概率,而预设门限r过低,又会增大检测算法的虚警率。
在图15所示的本实施案例的谐波图中,通过适当的选取门限值,例如r=10,将能成功检测到呼吸信号。
上述呼吸信号检测装置,采用慢时间域呼吸信号的谐波结构确定滤波参数,根据确定后的滤波参数滤波器对采集的原始信号进行滤波处理,若原始信号中存在呼吸信号,则通过此滤波处理可较好的滤除环境噪声且尽可能的保留呼吸信号,根据呼吸信号谐波结构及计算得到的慢时间域信号周期图和噪声功率计算谐波图,通过谐波图与预设门限进行判断,提高了检测的性能,更加有效的检测呼吸信号是否存在。
进一步的,如图18所示,在一个实施例中,上述呼吸信号检测装置,除了包括采集模块110、噪声抑制模块120、周期图计算模块130、噪声功率估算模块140、谐波图计算模块150和判断模块160,还包括距离和呼吸率估算模块180。其中:
距离和呼吸率估算模块180,用于在判断存在呼吸信号时,根据所述谐波图的最大值估算目标距离和/或目标呼吸率。
具体的,当判决存在呼吸信号时,距离和呼吸率估算模块180,对目标距离和呼吸率进行估算包括以下过程:
(1)获取Hn[k]的最大值,其对应的快时间域变量n用表示,其对应的离散时间频率用表示。
(2)将作为呼吸率的估计值。
(3)估计目标距离d。
具体的,按传统的雷达脉冲信号计算距离的方式结合快时间域滤波器对距离估计的影响可得,目标距离计算公式如下:
式(15)中,,delayftF为快时间域滤波器的群延时,c表示光速,Tft为雷达快时间域采样周期,表示快时间域变量n。目标距离是指被测目标与检测装置之间的距离。
在本实施案例中,图15中小箭头“”所指向的Hn[k]即为谐波图的最大值,根据该最大值即可按式(15)估计呼吸率和目标距离。
以上所述实施例仅表达了本发明的几种实施方式,其描述较为具体和详细,但并不能因此而理解为对本发明专利范围的限制。应当指出的是,对于本领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干变形和改进,这些都属于本发明的保护范围。因此,本发明专利的保护范围应以所附权利要求为准。