本发明属于信号处理技术领域，具体涉及一种生命信号的特征提取方法。
背景技术：
呼吸与心跳是生命体征信息的重要指标。一方面，心肺体征信息可用于判断有无生命体存在，以及生命体的基本状态；另一方面，心肺活动参数的异常往往伴随着医学上的突发紧急事件。因此实时地进行心肺体征监测在诸多场合有非常重要的实用价值。常见的生命体征信号检测方法中，用于人体呼吸检测的方法主要有：压力传感器法、温度传感器法、电阻抗式呼吸测量法、呼吸感应体积描记法；与心跳相关的检测方法有：心电图、心音、光电式脉搏波测量法等。
以上方法大多基于接触式，需要与人体皮肤相接触才能测量出生命体征参数，限制了其在一些特殊场合的应用。如在对老人的监护中，长时间佩戴仪器设备会引起监护对象的不适；对于大面积烧伤的患者，电极或传感器可能会造成二次伤害；在发生一些灾难后的搜救过程中，通过接触式的方法来检测生命体征显然是不切实际的，因此非接触式生命信号特征提取技术的发展有着非常重要的价值。
目前，利用毫米波雷达来检测生命信号特征是一种常见的方法，这是由于其具有穿透能力强、分辨率高、波长较短的特点。距离信息的微弱的变化能够使得毫米波雷达回波相位发生较大的变化，通过提取相位来进行对生命体征信号的检测是一种非常有效的方法。
技术实现要素：
发明目的：本发明旨在提供一种提取生命信号特征的方法，该方法能够在非接触的情况下，精确提取目标的呼吸频率和心跳频率。
技术方案：本发明采用如下技术方案：
一种基于线谱跟踪的生命信号特征提取方法，所述生命信号为呼吸信号和心跳信号，包括如下步骤：
(1)对目标的生命信号引起的胸腔运动进行建模，得到目标与雷达的距离为：
r(t)＝r0-a(t)
其中，r0为目标的初始距离，a(t)为生命信号引起的胸腔运动的模型；
(2)lfmcw毫米波雷达连续发射n个调频信号，对每一个雷达回波信号进行deramp混频处理，获得中频信号if；
利用傅里叶变换从中频信号中获取目标的位置和相位信息，n个调频信号的回波得到n个相位信息，构成长度为n的一维向量φ＝(φ1,φ2,…,φn)；其中φk为第k个调频信号回波得到的相位；
(3)设置长度为m的窗函数，重叠点数为l，对φ进行短时傅里叶变换，获得m×k维的时频谱矩阵ps，其中k＝fix([n-(m-l)])/l，fix表示取整操作；
(4)选取频率范围为[0.1,0.5]hz频谱数据，构建第一隐马尔科夫线谱跟踪模型γ1＝(ψ,ω,ξ)；其中ψ、ω、ξ分别为第一隐马尔科夫线谱跟踪模型的状态转移矩阵、观测概率矩阵和初始概率；
(5)将k列频谱ps均匀划分成q组，每组数据块为p列，将每组的p列数据作为观测量序列，利用viterbi算法计算第一隐马尔科夫线谱跟踪模型每组观测量序列对应的最优状态序列，将q组最优状态序列拼接起来，即得到各时刻的呼吸频率所在的频率单元；
(6)选取频率范围为[0.5,3]hz频谱数据，并且将时频谱矩阵ps中对应呼吸频率的多次谐波分量的元素置为零，得到第二时频谱矩阵ps′，构建第二隐马尔科夫线谱跟踪模型γ2＝(ψ,ω′,ξ)；其中ω′为第二隐马尔科夫线谱跟踪模型的观测概率矩阵；
(7)将k列频谱p′s均匀划分成q组，每组数据块为p列，将每组的p列数据作为观测量序列，利用viterbi算法计算第二隐马尔科夫线谱跟踪模型每组观测量序列对应的最优状态序列，将q组最优状态序列拼接起来，即得到各时刻的心跳频率所在的频率单元。
所述生命信号引起的胸腔运动的模型为：
其中，rai表示呼吸波形的i次谐波的幅度，fr表示呼吸频率，ha表示心跳波形的幅度，fh表示心跳频率，num为最大谐波次数。
所述步骤(2)具体包括以下步骤：
(2.1)lfmcw毫米波雷达系统发射的线性调频信号为：
s(t)＝σexp(j2πfct+jπγt2)
其中，σ为信号幅度，fc为载波频率，γ为线性调频斜率；
回波信号x(t)为：
其中ρk为第k个回波的散射系数，τ＝2r(t)/c为目标的回波时延，r(t)表示t时刻雷达与目标之间的距离，c为电磁波传播速度；利用发射信号对回波信号进行混频处理，即可获得if信号：
其中，f(t)表示中频信号，*为共轭处理；考虑到时延τ较小，忽略二次时延相位项，if信号可近似表示为：
(2.2)由于τ＝2r(t)/c，式(22)可以改写为：
其中λ＝1/fc，从中频信号中提取相位φb；
目标的距离与所提取出的相位的关系为：
(2.3)假设车载雷达系统采样频率为fs，调频连续波时宽为t，则一个时宽内采样点数s＝fst。对n个if信号进行采样，并将采样点按列存储，形成s*n大小的帧信号，其中第n个if信号的第s个采样点为f(s,n)；s＝0,1,2,…,s-1，n＝0,1,2,…,n-1；
将帧信号通过一维傅里叶变换获得一帧的距离fft图，可表示为：
其中，w(s)为高斯窗函数，u＝1,2,…,s。
(2.4)对n个回波信号做完距离fft之后，找出每个回波信号中目标所对应的距离单元rk，根据式(24)提取出n回波的相位信息并做相位解缠绕处理，构成长度为n的一维向量φ＝(φ1,φ2,…,φn)。
所述步骤(4)具体包括以下步骤：
(4.1)设置每个频率单元的初始概率默认为均匀分布的，即：
m表示窗函数长度，也即频率单元的总数；
(4.2)设状态转移概率满足均值为0方差为
在短时傅里叶变换得到的时频图中，第i个频率单元表示为：
[fi,fi+1]i＝1,…,m(28)
其中心频率为：
频率从第i个单元转移到第j个单元的概率为：
其中ε预设的状态最大转移范围，超过该范围的状态转移概率统一置为零；对状态转移概率进行归一化处理，即：
以
(4.3)根据时频谱矩阵ps计算观测概率矩阵ω，ω中的元素为：
其中zp为p时刻的观测量，zp＝1,…,m；xp为p时刻的状态量，即频率所在的频率单元序号，p＝1,2,…,k；pr(zp|xp＝i)表示p时刻的频率在第i个频率单元时，观测结果为zp的概率；
以ωi(zp)为第i行p列元素的m×k维矩阵为观测概率矩阵ω；
(4.4)第一隐马尔科夫线谱跟踪模型γ1表示为：
γ1＝(ψ,ω,ξ)(26)。
所述步骤(5)中计算ps第q个数据块作为观测量序列所对应的最优状态序列，q＝1,2,…,q，具体包括以下步骤：
(5.1)定义q(v)为第q个数据块第v列数据所对应的时刻；v＝1,2,3,...,p；
令v＝1，初始化状态空间中每个状态对应的局部状态δq(v)(i)和状态索引θq(v)(i)：
δq(1)(i)＝ξi·ωi(zq(1))；
θq(1)(i)＝0；
i＝1,2,...,m；
其中ξi为m×1维向量ξ的第i个元素，本发明中，ξi均为
(5.2)向后递推，对于v＝2,3,...,p时q(v)时刻的每一个状态sj，计算对应的局部状态δq(v)(i)，即从q(v)-1时刻状态为si到q(p)时刻状态为sj所对应的最大概率，并将q(v)-1时刻的状态索引存储到θq(v)(i)：
j＝1,2,...,m；
其中
(5.3)通过最大化δq(p)(j)找到在时刻q(p)的状态估计
(5.4)最优路径回溯：通过存储在θq(p)中的后向指针将时刻p状态的路径追溯到初始时刻：
得到最优状态序列
所述步骤(6)具体包括以下步骤：
(6.1)第二隐马尔科夫线谱跟踪模型的状态转移矩阵和初始概率与第一隐马尔科夫线谱跟踪模型的状态转移矩阵和初始概率相同，为ψ和ξ；
(6.2)根据第二时频谱矩阵ps′计算观测概率矩阵ω′，ω′中的元素为：
其中z′p为p时刻的观测量，z′p＝1,…,m；x′p为p时刻的状态量，即频率所在的频率单元序号，p＝1,2,…,k；pr(z′p|x′p＝i)表示p时刻的频率在第i个频率单元时，观测结果为z′p的概率；
以ω′i(z′p)为第i行p列元素的m×k维矩阵为观测概率矩阵ω′；
(6.3)第二隐马尔科夫线谱跟踪模型γ2表示为：
γ2＝(ψ,ω′,ξ)。
所述步骤(7)中计算p′s第q个数据块作为观测量序列所对应的最优状态序列，q＝1,2,…,q，具体包括以下步骤：
(7.1)定义q(v)为第q个数据块第v列数据所对应的时刻；v＝1,2,3,...,p；
令v＝1，初始化状态空间中每个状态对应的局部状态δq(v)(i)和状态索引θq(v)(i)：
δq(1)(i)＝ξi·ω′i(z′q(1))；
θq(1)(i)＝0；
i＝1,2,...,m；
其中ξi为m×1维向量ξ的第i个元素，本发明中，ξi均为
(7.2)向后递推，对于v＝2,3,...,p时q(v)时刻的每一个状态sj，计算对应的局部状态δq(v)(i)，即从q(v)-1时刻状态为si到q(p)时刻状态为sj所对应的最大概率，并将q(v)-1时刻的状态索引存储到θq(v)(i)：
j＝1,2,...,m；
其中
(7.3)通过最大化δq(p)(j)找到在时刻q(p)的状态估计
(7.4)最优路径回溯：通过存储在θq(p)中的后向指针将时刻p状态的路径追溯到初始时刻：
得到最优状态序列
为了提高抗干扰性，所述生命信号引起的胸腔运动的模型为：
其中ra表示呼吸信号基波的幅度，fr表示呼吸频率，ha表示心跳波形的幅度，fh表示心跳频率，num为最大谐波次数，mb(t)为人体抖动信号：
其中n表示有身体抖动的时间单元的数目，t1,...,tn分别为每个时间单元的时长，a1,...,an分别为每个时间单元的最大抖动幅度。
有益效果：与现有技术相比，本发明公开的基于线谱跟踪的生命信号特征提取方法利用时间的相关性来跟踪生命体征信号的线谱并提取出相对应的频率，可以在非接触的情况下精确提取到目标的呼吸和心跳频率。
附图说明
图1为本发明公开的生命信号特征提取方法的流程图；
图2为生命体征信号的波形图；
图3为对回波信号做距离fft之后的结果图；
图4为提取到的相位信息做解缠绕之后的处理图；
图5为生命体征信号的线谱跟踪曲线图；
图6为加入模拟人体抖动信号后的生命体征信号波形图；
图7为加入模拟人体抖动信号后的线谱跟踪曲线图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面结合附图对本发明的具体实施案例做说明。
实施例1：
如图1所示，本发明公开了一种基于线谱跟踪的生命信号特征提取方法，本发明中的生命信号为呼吸信号和心跳信号，提取的是呼吸信号和心跳信号的频率。
该方法包括如下步骤：
步骤1、对目标的生命信号引起的胸腔运动进行建模，得到目标与雷达的距离为：
r(t)＝r0-a(t)
其中，r0为目标的初始距离，a(t)为生命信号引起的胸腔运动的模型；
本实施例中，生命信号引起的胸腔运动的模型为：
其中，rai表示呼吸波形的i次谐波的幅度，fr表示呼吸频率，本实施例中为20次/分钟，变化幅度大小为2次，即20±2次/分钟；ha表示心跳波形的幅度，fh表示心跳频率，本实施例中为72次/分钟，变化幅度大小为5次；num为最大谐波次数。初始的呼吸波形ra最大值为6mm，变化幅度大小为3mm，n次呼吸谐波的最大值为ra/2n(n>＝2)，初始的心跳波形最大值ha为1mm，变化幅度大小为0.1mm。由此模拟出的胸腔运动波形如图2所示。
目标距离雷达的初始距离r0为0.9m，则目标与雷达的距离为：
本实施例中num的取值为4。
步骤2、lfmcw毫米波雷达连续发射n个调频信号，对每一个雷达回波信号进行deramp混频处理，获得中频信号；
利用傅里叶变换从中频信号中获取目标的位置和相位信息，n个调频信号的回波得到n个相位信息，构成长度为n的一维向量φ＝(φ1,φ2,…,φn)；其中φk为第k个调频信号回波得到的相位；具体步骤为：
(2.1)lfmcw毫米波雷达系统发射的线性调频信号为：
s(t)＝σexp(j2πfct+jπγt2)
其中，σ为信号幅度，fc为载波频率，本实施例中为77ghz，γ为线性调频斜率；γ＝b/t，本实施例中带宽b取值为2000mhz，调频连续波时宽t为50μs；
回波信号x(t)为：
其中ρk为第k个回波的散射系数，τ＝2r(t)/c为目标的回波时延，r(t)表示t时刻雷达与目标之间的距离，c为电磁波传播速度；利用发射信号对回波信号进行混频处理，即可获得if信号：
其中，f(t)表示中频信号，*为共轭处理；考虑到时延τ较小，忽略二次时延相位项，if信号可近似表示为：
其中，多普勒频率
(2.2)由于τ＝2r(t)/c，式(22)可以改写为：
其中λ＝1/fc，从中频信号中提取相位φb；
目标的距离与所提取出的相位的关系为：
(2.3)假设车载雷达系统采样频率为fs，调频连续波时宽为t，则一个时宽内采样点数s＝fst。对n个if信号进行采样，并将采样点按列存储，形成s*n大小的帧信号，其中第n个if信号的第s个采样点为f(s,n)；s＝0,1,2,…,s-1，n＝0,1,2,…,n-1；本实施例中fs为5mhz，n为10240；
将帧信号通过一维傅里叶变换获得一帧的距离fft图，可表示为：
其中，w(s)为高斯窗函数，u＝1,2,…,s。
(2.4)对n个回波信号做完距离fft之后，找出每个回波信号中目标所对应的距离单元rk，本实施例中选取幅度最大的距离单元作为目标单元，如图3所示，图像中灰度最大的距离单元即为目标单元。根据式(24)提取出n个回波的相位信息并做相位解缠绕处理，如图4所示，构成长度为n的一维向量φ＝(φ1,φ2,…,φn)。
步骤3、设置长度为m的窗函数，重叠点数为l，对φ进行短时傅里叶变换(shorttimefouriertransform,stft)，获得m×k维的时频谱矩阵ps，其中k＝fix([n-(m-l)])/l，fix表示取整操作；本实施例中，m＝1024，重叠点数为l＝60；计算得到k＝154；
步骤4、选取频率范围为[0.1,0.5]hz频谱数据，构建第一隐马尔科夫线谱跟踪模型γ1＝(ψ,ω,ξ)；其中ψ、ω、ξ分别为第一隐马尔科夫线谱跟踪模型的状态转移矩阵、观测概率矩阵和初始概率；具体包括以下步骤：
(4.1)设置每个频率单元的初始概率默认为均匀分布的，即：
m表示窗函数长度，也即频率单元的总数；
(4.2)设状态转移概率满足均值为0方差为
在短时傅里叶变换得到的时频图中，第i个频率单元表示为：
[fi,fi+1]i＝1,…,m(28)
其中心频率为：
频率从第i个单元转移到第j个单元的概率为：
其中ε预设的状态最大转移范围，超过该范围的状态转移概率统一置为零；对状态转移概率进行归一化处理，即：
以
(4.3)根据时频谱矩阵ps计算观测概率矩阵ω，ω中的元素为：
其中zp为p时刻的观测量，zp＝1,…,m；xp为p时刻的状态量，即频率所在的频率单元序号，p＝1,2,…,k；pr(zp|xp＝i)表示p时刻的频率在第i个频率单元时，观测结果为zp的概率；
以ωi(zp)为第i行p列元素的m×k维矩阵为观测概率矩阵ω；
(4.4)第一隐马尔科夫线谱跟踪模型γ1表示为：
γ1＝(ψ,ω,ξ)(26)。
步骤5、将k列频谱ps均匀划分成q组，每组数据块为p列，将每组的p列数据作为观测量序列，利用viterbi算法计算第一隐马尔科夫线谱跟踪模型每组观测量序列对应的最优状态序列，将q组最优状态序列拼接起来，即得到各时刻的呼吸频率所在的频率单元；
计算ps第q个数据块作为观测量序列所对应的最优状态序列，q＝1,2,…,q，具体包括以下步骤：
(5.1)定义q(v)为第q个数据块第v列数据所对应的时刻；v＝1,2,3,...,p；
令v＝1，初始化状态空间中每个状态对应的局部状态δq(v)(i)和状态索引θq(v)(i)：
δq(1)(i)＝ξi·ωi(zq(1))；
θq(1)(i)＝0；
i＝1,2,...,m；
其中ξi为m×1维向量ξ的第i个元素，本发明中，ξi均为
(5.2)向后递推，对于v＝2,3,...,p时q(v)时刻的每一个状态sj，计算对应的局部状态δq(v)(i)，即从q(v)-1时刻状态为si到q(p)时刻状态为sj所对应的最大概率，并将q(v)-1时刻的状态索引存储到θq(v)(i)：
j＝1,2,...,m；
其中
(5.3)通过最大化δq(p)(j)找到在时刻q(p)的状态估计
(5.4)最优路径回溯：通过存储在θq(p)中的后向指针将时刻p状态的路径追溯到初始时刻：
得到最优状态序列
步骤6、选取频率范围为[0.5,3]hz频谱数据，并且将时频谱矩阵ps中对应呼吸频率的多次谐波分量的元素置为零，得到第二时频谱矩阵p′s，构建第二隐马尔科夫线谱跟踪模型γ2＝(ψ,ω′,ξ)；其中ω′为第二隐马尔科夫线谱跟踪模型的观测概率矩阵；
具体包括以下步骤：
(6.1)第二隐马尔科夫线谱跟踪模型的状态转移矩阵和初始概率与第一隐马尔科夫线谱跟踪模型的状态转移矩阵和初始概率相同，为ψ和ξ；
(6.2)根据第二时频谱矩阵ps′计算观测概率矩阵ω′，ω′中的元素为：
其中z′p为p时刻的观测量，z′p＝1,…,m；x′p为p时刻的状态量，即频率所在的频率单元序号，p＝1,2,…,k；pr(z′p|x′p＝i)表示p时刻的频率在第i个频率单元时，观测结果为z′p的概率；
以ω′i(z′p)为第i行p列元素的m×k维矩阵为观测概率矩阵ω′；
(6.4)第二隐马尔科夫线谱跟踪模型γ2表示为：
γ2＝(ψ,ω′,ξ)。
步骤7、将k列频谱p′s均匀划分成q组，每组数据块为p列，将每组的p列数据作为观测量序列，利用viterbi算法计算第二隐马尔科夫线谱跟踪模型每组观测量序列对应的最优状态序列，将q组最优状态序列拼接起来，即得到各时刻的心跳频率所在的频率单元。
计算p′s第q个数据块作为观测量序列所对应的最优状态序列的具体步骤与步骤5类似。
本实施例中提取的呼吸频率和心跳频率结果如图5所示。
与真实的呼吸频率和心跳频率进行对比，利用均方误差mse来评估效果，得到的结果为：
其中
实施例2
为了检验本发明的抗干扰性能，在实施例1的基础上加入了对人体抖动的模拟，具体步骤如下：
假设人体的抖动信号为三角波，表达形式如下：
其中n表示有身体抖动的时间单元的数目，本实施例取总时间单元的20％，每个时间单元为4s；t1,...,tn分别为每个单元的时长，本例中均设置为4s；a1,...,an分别为每个单元的最大抖动幅度，本例中设置为幅度在[010]mm之间的随机值，抖动信号如图6所示。
此时生命信号引起的胸腔运动的模型为：
目标与雷达的距离为：
经过处理后得到的线谱跟踪结果如图7所示，计算可得呼吸频率与心跳频率的均方误差mse分别为0.0048、0.0450，因此本发明在干扰条件下，依旧能够对生命体征信号具有很好的检测性能。