CN110309817A - 一种参数自适应优化vmd的脉搏波运动伪影去除方法 - Google Patents
一种参数自适应优化vmd的脉搏波运动伪影去除方法 Download PDFInfo
- Publication number
- CN110309817A CN110309817A CN201910652889.4A CN201910652889A CN110309817A CN 110309817 A CN110309817 A CN 110309817A CN 201910652889 A CN201910652889 A CN 201910652889A CN 110309817 A CN110309817 A CN 110309817A
- Authority
- CN
- China
- Prior art keywords
- decomposition
- pulse wave
- eigenmode
- mode
- frequency
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000005457 optimization Methods 0.000 title claims abstract description 28
- 230000003044 adaptive effect Effects 0.000 title claims abstract description 20
- 238000005516 engineering process Methods 0.000 title abstract description 5
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 121
- 238000000034 method Methods 0.000 claims abstract description 71
- 239000006185 dispersion Substances 0.000 claims abstract description 12
- 230000008569 process Effects 0.000 claims description 26
- 238000005314 correlation function Methods 0.000 claims description 13
- 238000005070 sampling Methods 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 5
- 239000011159 matrix material Substances 0.000 claims description 5
- 238000001914 filtration Methods 0.000 claims description 4
- 238000012216 screening Methods 0.000 claims description 4
- 230000003595 spectral effect Effects 0.000 claims description 4
- 238000001228 spectrum Methods 0.000 claims description 4
- 230000000241 respiratory effect Effects 0.000 claims description 3
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 claims description 2
- 230000008030 elimination Effects 0.000 claims description 2
- 238000003379 elimination reaction Methods 0.000 claims description 2
- 230000000694 effects Effects 0.000 abstract description 10
- 238000012545 processing Methods 0.000 abstract description 7
- 239000000203 mixture Substances 0.000 abstract 1
- 230000001133 acceleration Effects 0.000 description 5
- 230000036541 health Effects 0.000 description 5
- 238000012544 monitoring process Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 238000009499 grossing Methods 0.000 description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 2
- QVGXLLKOCUKJST-UHFFFAOYSA-N atomic oxygen Chemical compound [O] QVGXLLKOCUKJST-UHFFFAOYSA-N 0.000 description 2
- 239000008280 blood Substances 0.000 description 2
- 210000004369 blood Anatomy 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 229910052760 oxygen Inorganic materials 0.000 description 2
- 239000001301 oxygen Substances 0.000 description 2
- 208000024172 Cardiovascular disease Diseases 0.000 description 1
- 230000003190 augmentative effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 210000000748 cardiovascular system Anatomy 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 210000000624 ear auricle Anatomy 0.000 description 1
- 230000002526 effect on cardiovascular system Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000012880 independent component analysis Methods 0.000 description 1
- 238000002329 infrared spectrum Methods 0.000 description 1
- 230000001678 irradiating effect Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 230000031700 light absorption Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
- G16H40/00—ICT specially adapted for the management or administration of healthcare resources or facilities; ICT specially adapted for the management or operation of medical equipment or devices
- G16H40/60—ICT specially adapted for the management or administration of healthcare resources or facilities; ICT specially adapted for the management or operation of medical equipment or devices for the operation of medical equipment or devices
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
Landscapes
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- Biomedical Technology (AREA)
- Business, Economics & Management (AREA)
- General Business, Economics & Management (AREA)
- Epidemiology (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Primary Health Care (AREA)
- Public Health (AREA)
- Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)
Abstract
本发明涉及一种参数自适应优化VMD的脉搏波运动伪影去除方法,属于生物医学信号采集和处理技术领域;该方法通过优化变分模态分解的目标函数,并采用本征模态相关程度和中心频率离散程度相融合的相关系数对变分模态分解结果进行评价,实现对脉搏波信号的参数自适应的优化变分模态分解,得到脉搏波的本征模态{uk},根据信号不同成分对应频段的能量相对大小,设计去噪规则完成对脉搏波运动伪影的去除。对比现有技术,本发明提高了脉搏波变分模态分解的质量,提高了运动伪影去除的效果,且简化了去除运动伪影过程中的数据复杂度。
Description
技术领域
本发明涉及一种参数自适应优化VMD的脉搏波运动伪影去除方法,属于生物医学信号采集和处理技术领域。
背景技术
大量临床实测结果显示,脉搏波(PPG)包含丰富的心血管系统生理及病理信息,其波形特征与心血管疾病有着密切的联系。脉搏在临床中通常采用光电式脉搏传感器进行测量,其信号反映为光电容积脉搏波,是一种无创的测量方式。用位于可见红光光谱及红外光谱的两个光源交替照射检测部位(一般为耳垂或指尖),通过反射或透射得到脉动期间的组织的光吸收量,从而获得脉搏波的基础波形、脉率估计和血氧浓度等生理信号。
脉搏波的连续监测是当前家居健康养老、远程健康监护等系统中基础的必备测量数据,可以为用户的健康监控提供心率、血氧浓度等多种生理指标,但运动伪影对脉搏波的影响一直是高效可靠地使用脉搏血氧仪进行连续实时健康监测的重大障碍。运动中的脉搏实时采集能够通过提供心率估计帮助人们及时了解自身的健康情况并应对突发状况,从而降低运动风险,但在运动中由于脉搏波采集设备与皮肤间距的不规律变化导致漏光等现象,使得检测到的脉搏波信号存在大量的运动伪影无法去除。总之,运动伪影的存在对脉搏波的后续分析产生了不可避免的影响。
最初的运动伪影去除方法包括独立成分分析法、自适应滤波技术等算法,但这些方法普遍存在以下问题:第一,两种算法都建立在运动伪影与脉搏波信号完全不相关的基本假设基础之上,这一假设与真实状况并不相符;第二,当运动伪影与脉搏波信号发生本征模态混叠时,算法的去噪效果变差。在随后的发展中,为解决上述问题,经验本征模态分解(EMD)和集合经验本征模态分解(EEMD)等方法被应用于脉搏波信号,通过将信号分解为几个固有本征模态分量将运动伪影与脉搏波分离,达到去除运动伪影的目的,但算法缺少数学基础,可解释性差,并且在信号中存在较强运动伪影的情况下失效。
变分模态分解(VMD)是现有算法中较好的运动伪影去除算法,该算法设置合适的目标函数,通过对带约束求极值问题的变分求解,将信号分解为若干个本征模态,并在此基础上完成信号的去噪过程。但该算法现在还存在以下问题:第一,VMD算法以各个本征模态带宽最窄为目标,但其目标函数的设定没有充分考虑到分解结果中本征模态的平滑程度;第二,VMD算法中的分解本征模态数K、二次惩罚因子α一般通过经验进行提前设置,无法确认算法的最优性;第三,现有对运动伪影的去除方法一般依靠加速度信号的辅助,去噪过程需要同时采集运动信号,导致涉及到的数据量较多。
本发明针对上述存在的缺陷,提出一种参数自适应优化变分模态分解技术,致力于提高脉搏波变分模态分解的质量并简化去除运动伪影过程中的数据复杂度。
发明内容
针对现有技术中存在的不足,本发明的目的在于提出一种参数自适应优化变分模态分解技术,以完成对脉搏波信号的分解并去除运动伪影。在参数自适应选取部分,通过对变分模态分解结果的分析,本发明提出采用各个本征模态间的相关系数作为指标函数,评价分解结果的好坏,从而克服本征模态混叠问题,实现参数的自适应;在变分模态分解的具体算法实施中,本发明对目标函数进行改进,考虑变分模态的平滑程度,提出新的的优化目标函数,从而使算法针对脉搏波具有更好的效果;在运动伪影去除过程中,本发明提出基于能量分布的去噪方法,仅依靠脉搏波数据完成去噪,在不影响运动伪影去除效果的同时,相比传统方法中借助加速度相关数据去噪的方法,降低了数据的繁杂程度。
本发明旨在解决脉搏波信号中运动伪影去除的问题,提出了一种参数自适应优化VMD的脉搏波运动伪影去除方法。
一种参数自适应优化VMD的脉搏波运动伪影去除方法,包括如下步骤:
步骤1:首先对变分模态分解的目标函数进行优化设计,加入保证本征模态平滑性质的优化项,得到优化变分模态分解的目标函数;
作为优选,步骤1通过以下过程实现:
变分模态分解算法是一种以变分问题为框架的信号分解估计方法,该算法以本征模态的带宽估计函数取得最小为目标函数对输入信号进行最优求解,考虑到分解结果中本征模态的平滑程度,本发明首先设计如下的目标函数:
其中,min为取最小值,{uk}表示优化变分模态分解结果中的本征模态集合,其元素uk表示分解结果中的本征模态,{ωk}表示优化变分模态分解结果中本征模态对应的中心频率集合,其元素ωk表示分解结果中本征模态对应的中心频率,K为分解本征模态个数;
为uk的带宽估计;其中表示括号中的项对时间求梯度,δ(t)为脉冲函数,j为虚数单位,*表示卷积运算,||·||2表示二范数;
增加项的目的在于保证分解结果中各个本征模态的平滑性质;γ(1≥γ≥0)为正则加权因子,time为所处理脉搏波信号的采样点数,Δuk(t)表示本征模态uk在t处的梯度,定义为:Δuk(t)=uk(t)-uk(t-1),uk(t)为本征模态uk在采样点t处的幅值,||·||1表示1范数;
约束条件为∑uk(t)=f,其意义为分解所得所有本征模态之和等于输入信号f。
步骤2:通过引入拉格朗日乘子和二次惩罚因子,对步骤1中设计的目标函数进行去约束求解,通过对各个本征模态、本征模态对应的中心频率以及拉格朗日乘子进行循环更新,完成对脉搏波信号的优化变分模态分解,求得信号的本征模态{uk};
作为优选,步骤2通过以下过程实现:
步骤2.A:初始化变分模态分解算法相关参数;
其中,U0=[u1,0,u2,0,...,ui,0,...,uK,0]T为K个本征模态的初始化,ui,0表示第i个本征模态的初始化;ω0=[ω1,0,ω2,0,…,ωi,0,…,ωK,0]T为本征模态的中心频率初始化,ωi,0表示第i个本征模态对应的中心频率初始化,λ0为拉格朗日乘子的初始化;
步骤2.B:对所求参数进行循环更新;
首先引入拉格朗日乘子和二次惩罚因子,将公式(1)所示的带约束极值问题通过增广拉格朗日方程转化为无约束极值问题,再对信号及参数进行傅里叶变换后,运用交替方向乘数法直接在频域内进行求解;经计算,求得公式(2)~(4),用于对本征模态、中心频率以及拉格朗日乘子进行循环更新:
其中,ω为时域信号经傅里叶变换后在频域中的频率;为第k个本征模态的傅里叶变换形式在第n+1次迭代更新后的结果;为输入信号的傅里叶变换形式;为第i个本征模态的傅里叶变换形式在第n+1次迭代更新后的结果;为第i个本征模态的傅里叶变换形式在第n次迭代更新后的结果;为拉格朗日乘子的傅里叶变换形式在第n次迭代更新后的结果;ωk,n为第k个本征模态的的中心频率在第n次迭代更新后的结果;α为二次惩罚因子,其作用为改善结果的收敛性;τ为拉格朗日算子的更新参数;
步骤2.C:依据工程需要,设定阈值ε,判断筛选条件是否成立,并输出计算结果;
计算筛选条件为:
其中,ε为误差允许范围参数;
若公式(5)成立,则步骤2.B中计算出的本征模态及中心频率即为输出的uk,ωk;
若公式(5)不成立,则返回步骤2.B。
步骤3:在分解本征模态个数K和二次惩罚因子α的不同组合下,对脉搏波信号基于步骤2完成优化变分模态分解,设计参数最优选取的目标函数D,通过计算不同参数下分解结果对应的目标函数对参数进行自适应选取,得到优化变分模态分解的最优参数K和α,以及在该组参数下的分解结果,即本征模态{uk};
作为优选,步骤3通过以下过程实现:
步骤3.A:初始化相关参数,并在不同的参数组合下通过步骤2完成对脉搏波信号的变分模态分解;
初始化[Kmin,Kmax,Kstep]、[αmin,αmax,αstep],其中,Kmin为本征模态分解个数K的最小取值,Kmax为本征模态分解个数的最大取值,Kstep为本征模态分解个数的寻优步长;αmin为二次惩罚因子α的最小取值、αmax为二次惩罚因子α的最大取值,αstep为二次惩罚因子α的寻优步长;
通过初始化可以得到分解本征模态个数K和二次惩罚因子α的不同组合方式,构成维的取值矩阵,通过对矩阵遍历选择不同的参数组合通过步骤2对脉搏波信号进行变分模态分解;
步骤3.B:计算不同参数下分解结果对应的关联函数D;
在完成对脉搏波变分模态分解的过程中,可能存在欠前分解导致的模态混叠问题或过分解问题,从变分模态分解的结果考虑,期望分解结果中各个本征模态的相关程度达到最小,同时期望各个本征模态中心频率的离散程度达到最大,本发明综合考虑到以上两点要求,提出参数最优选取的关联函数如下:
其中,第一项中为各个模态间相关系数的平均值,用来衡量分解结果中模态间的相关性,其取值范围为[0,1],为分解结果中模态ui和模态uj的相关系数,分别表示模态ui和模态uj的平均值,E[·]表示·的期望;
第二项中用来衡量分解结果中各个模态中心频率的离散程度,为中心频率离散程度的度量,加入负指数的目的是保证模态相关程度和中心频率分散程度对最优解取极值的一致性,即当目标函数D取得最小值时,参数取得最优解,同时保证其取值的非负性及有界性,取值范围为[0,1];
公式(6)中的可调系数η用来调节对本征模态相关程度和中心频率离散程度的约束强度,η的参数选择范围为[0,1];
关联函数D的取值范围为[0,1],当D取0时,表示本征模态间各不相关,并且各个模态中心频率的离散程度达到无穷大(实际上不会存在这种取值),当D取1时,表明各个模态完全重叠;在对信号完成变分模态分解的过程中,期望模态间的相关程度达到最小,中心频率离散程度达到最大,即求取使关联函数D取得最小值的参数组合。
步骤3.C:选取最优参数;
对于在每一组特定的参数组合下完成的脉搏波信号优化变分模态分解结果,通过上述过程即公式(6)计算关联函数D,直到完成所有参数组合的遍历,最终选取使关联函数D取得最小值的分解本征模态个数K和二次惩罚因子α的参数组合,并保留该参数下的分解结果{uk}。
步骤4:在得到步骤3最优参数组合及该参数组合下的分解结果后,基于各个本征模态的能量谱完成运动伪影相关本征模态的去除,而后对剩余本征模态进行求和重构,得到完成去噪的脉搏波信号;
作为优选,步骤4通过以下过程实现:
步骤4.A:计算各个本征模态在不同成分频段的能量;
能量谱密度反映了信号能量与频率的分布关系,对于脉搏波信号来说,信号中脉搏波动成分对应的频率范围为0.5~4Hz,呼吸成分对应的频率范围为0.2~0.35Hz,运动伪影部分由于运动形式不确定无法准确界定其频率边界,一般对应0.1Hz左右的低频分量或更高的频段;根据不同成分的能量分布情况,将每个本征模态的能量分为三个部分并利用模态的能量谱密度计算相应频段的能量:
其中Ek1为本征模态uk中脉搏波动对应频段的能量,Ek2为本征模态uk中呼吸成分对应频段的能量,Ek1和Ek2均通过对相应频段的能量谱密度对频率积分计算,为本征模态uk的能量谱密度,由于变分模态分解的求解过程在频域中完成,故对脉搏波进行变分模态分解所得结果即为变分模态分解所得本征模态uk的频域表示Ek3为本征模态uk中运动伪影对应频段的能量,由于运动伪影成分的频率不确定,故通过本征模态总能量与其他分量能量之差计算;
步骤4.B:对运动伪影相关本征模态或频段进行去除;
通过对各个本征模态在不同成分频段能量的比较和分析,利用以下规则对各个本征模态进行更新,完成运动伪影相关模态或频段的去除:
其中,θ为设定的固定阈值,可根据具体需要在(0,1)范围内取值,用来界定运动伪影成分与脉搏波动成分的相对大小;
对公式(8)具体描述如下:
当运动伪影成分能量远大于脉搏波动成分能量时,即对应公式(8)中的Ek1/Ek3≤θ,认为该本征模态的主要构成成分为运动伪影,将该本征模态置零;当脉搏波动成分能量与运动伪影成分能量相差不大时,即对应公式(8)中的θ<Ek1/Ek3<1/θ,认为该本征模态中同时包含脉搏波有效成分及运动伪影成分,对该本征模态进行处理更新为当脉搏波动成分能量远大于运动伪影成分能量时,即对应公式(8)中的Ek1/Ek3≥1/θ,认为该本征模态的主要构成成分为脉搏波有效信息,保留该本征模态;
为满足θ<Ek1/Ek3<1/θ条件时,更新后的本征模态,其计算方法为:只保留该模态脉搏波动成分对应的频段,去除其他频段的信号,此步骤相当于对该本征模态进行带通滤波,在频域中完成,通过公式(9)即可求得
步骤4.C:将步骤4.B中完成运动伪影相关模态去除的剩余模态进行求和重构,得到完成去噪的脉搏波信号;
由于在去噪过程中对既含有运动伪影成分又含有脉搏波动成分的模态,即满足θ<Ek1/Ek3<1/θ条件的模态进行了带通滤波更新,造成其在0.5Hz和4Hz处的阶跃性质,故在重构前考虑这些本征模态在该频率处的连续性,对信号不连续处进行拟合,通过公式(10)和公式(11)完成:
其中,Δ为取值较小的频率增量,用来在不连续点前后提供拟合区间,可根据实际情况进行选取;公式(10)和公式(11)的意义为:在不连续点处通过区间端点的信号取值将该点前后小区间内的信号拟合为直线,以保证信号的连续性;
对更新后的本征模态进行求和重构,由公式(12)完成模态在频域中的求和重构:
其中,u'k(ω)为步骤4.B中通过公式(8)以及本步骤中公式(10)~(11)完成更新和不连续点拟合的本征模态,fc(ω)为去除运动伪影后的脉搏波频域信号;
最后通过公式(13)对频域重构信号进行傅里叶逆变换,即可得到去除噪声及运动伪影的脉搏波信号;
其中,fc(t)为去除运动伪影后的脉搏波信号,表示·的傅里叶逆变换,本步骤通过对去除运动伪影相关本征模态或频段后更新的本征模态进行求和运算得到去噪后的脉搏波信号。
至此,由步骤4计算得到的去除运动伪影后的脉搏波信号;由步骤1到步骤4完成了一种参数自适应优化VMD的脉搏波运动伪影去除方法。
有益效果
一种参数自适应优化VMD的脉搏波运动伪影去除方法,相比于现有技术,本方法有如下有益效果:
1.对VMD算法的目标函数及约束条件进行了更新,本发明中的目标函数增加了对脉搏波分解结果中本征模态平滑程度的约束,避免在分解结果中出现大量的带毛刺信号或突变信号,使得到的本征模态更加平滑,分解结果更具解释性,能够提高脉搏波变分模态分解的质量,进而提高运动伪影去除的效果。
2.对脉搏波信号进行变分模态分解后得到的各个本征模态对应于脉搏波信号的不同频率尺度,本征模态数K是分解过程中一个重要的参数,若K值选取过大引起过分解现象,则会导致本征模态的中心频率接近,使本征模态无法一一对应脉搏波信号的频率特征;若K值选取过小引起模态混叠现象,即某些相邻模态无法分离,可能导致无法良好去除运动伪影的现象。对VMD算法的分解本征模态数K值进行自适应选取,避免了由于K值选取不当引起的上述现象。二次惩罚因子α也是一个重要参数,若α取值过大,则不适合准确捕获模式的中心频率;若取值较小,将导致提取模式对噪声的鲁棒性的折衷。对二次惩罚因子α进行自适应选取,可以确保算法良好的收敛性,优化脉搏波变分模态分解的过程。总之,对参数的自适应选取能够消除依据经验选取参数的人为误差,使算法的结果达到最优。本发明同时考虑本征模态相关系数和模态中心频率离散程度,提出关联函数并对信号分解结果进行评价,能够使分解结果中所有模态之间的相关性尽可能小,有效地减少本征模态混叠现象的发生,确保最优的参数选择,提高了脉搏波变分模态分解的质量。
3.对运动伪影进行去除的部分不依靠加速度信号,而是对本征模态的能量的频域分布进行分析,通过计算不同成分对应频段的能量大小判断本征模态的成分属性,完成运动伪影的去除,在确保去噪效果的前提下减少了数据的繁杂程度,从而降低对处理器和寄存器的要求,提高数据实时处理的性能,能够更适用于对数据实时处理要求较高的场景。比如若在脉搏血氧仪中植入本发明中的参数自适应优化VMD的脉搏波运动伪影去除方法,则可以在测量过程脉搏波数据的同时以较快速度完成对运动伪影的去除,从而得到干净的脉搏波信号,方便对信号进一步的处理及分析,已达到更好的监测效果。
附图说明
图1为本发明“一种参数自适应优化VMD的脉搏波运动伪影去除方法”中的本方法及实施例的流程图。
图2为本发明“一种参数自适应优化VMD的脉搏波运动伪影去除方法”中的本方法及实施例中的优化变分模态分解算法的流程示意图。
图3为本发明“一种参数自适应优化VMD的脉搏波运动伪影去除方法”中的本方法及实施例中的脉搏波运动伪影去除算法流程示意图。
具体实施方式
本实施例中的数据来源于中国人民解放军总医院心血管内科的患者采样数据,采用脉搏血氧仪对患者进行脉搏波采集,采样频率为500Hz,采样时间为1分钟,获取患者在日常自由活动(非剧烈运动)情境下的脉搏波数据。
下面根据附图和实例对本发明进行详细说明,但本发明的具体实施方式不仅于此。
实施例1
本实施例阐述了将本发明“一种参数自适应优化VMD的脉搏波运动伪影去除方法”应用于脉搏波去除运动伪影场景下的流程。
图1为本方法的流程图及本实施例的流程图,从图中可以看出,本方法包含如下步骤:
步骤A:初始化相关参数,对参数取值矩阵进行遍历,在每一种参数组合下完成信号的优化变分模态分解;
图2为本步骤和本实施例中的优化变分模态分解算法的流程示意图;
步骤A.1:初始化变分模态分解算法相关参数;
初始化[Kmin,Kmax,Kstep]、[αmin,αmax,αstep],其中,Kmin为本征模态分解个数K的最小取值,Kmax为本征模态分解个数的最大取值,Kstep为本征模态分解个数的寻优步长;αmin为二次惩罚因子α的最小取值、αmax为二次惩罚因子α的最大取值,αstep为二次惩罚因子α的寻优步长;
具体到本实施例,根据脉搏波在进行变分模态分解时的一般参数,初始化的具体参数值如下:
[Kmin,Kmax,Kstep]=[3,8,1]
[αmin,αmax,αstep]=[1500,2500,50]
根据参数的初始化值,由步骤3.A可生成一个6×21的参数取值矩阵;
初始化U0、ω0、λ0,其中,U0=[u1,0,u2,0,…,uK,0]T为K个本征模态的初始化,ui,0表示第i个本征模态的初始化;ω0=[ω1,0,ω2,0,...,ωK,0]T为本征模态的中心频率初始化,ωi,0表示第i个本征模态对应的中心频率初始化,λ0为拉格朗日乘子的初始化;
具体到本实施例,采用全零法对参数进行初始化,即
ui,0=[0]1×time,i=1,2,...,K
ωi,0=0,i=1,2,...,K
λ0=0
其中time为原始信号的采样点数,具体到本实施例,time的取值为60秒×500Hz=30000点。
步骤A.2:对所求参数进行循环更新;
遍历步骤A.1中的参数矩阵,对设计的优化变分模态分解目标函数,即公式(1)进行求解,具体到本实施例,公式(1)中的正则因子γ=1;通过引入拉格朗日乘子和二次惩罚因子,对信号及参数进行傅里叶变换后,运用交替方向乘数法直接在频域内对约束极值问题进行求解;通过公式(2)~(4)对本征模态、中心频率以及拉格朗日乘子进行循环更新,循环更新的方法与步骤2.B的计算方法相同;
步骤A.3:判断停止条件并输出结果;
判断停止条件是否成立,判断方法和规则与步骤2.C中方法相同,根据公式(5)计算停止条件是否成立,若公式(5)成立,则步骤A.2中最后一次更新所得本征模态及中心频率即为输出的uk,ωk;若公式(5)不成立,则返回步骤A.2。
步骤B:对分解本征模态个数K和二次惩罚因子α进行自适应选取;
对步骤A.1中生成的参数取值矩阵进行遍历,在每一组参数取值组合下按照步骤A.2完成优化变分模态分解,并进行相应计算得到最优参数组合,具体计算步骤如下:
步骤B.1:计算不同参数下分解结果对应的关联函数D;
对于每一组参数组合下由步骤A.3中获得的uk,计算步骤3.C中提出的目标函数,以评价该参数下分解结果的好坏,通过公式(6)进行计算,计算方法与步骤3.C中的方法相同,具体到本实施例中,公式(6)中的η=0.5;
步骤B.2:选取最优参数;
对每一组特定的参数组合完成的脉搏波信号优化变分模态分解结果通过步骤B.1计算相关系数,直到完成所有参数取值矩阵的遍历,最终选取使关联函数D值取得最大的分解本征模态个数K和二次惩罚因子α的组合,并保留该参数下的分解结果;
具体到本实施例中,通过对采集到的信号的分解,得到的最优参数为K=7,α=1850。
步骤C:对脉搏波分解后的结果进行运动伪影相关本征模态的去除;
图3为本步骤和本实施例中的脉搏波运动伪影去除算法流程示意图。
步骤C.1:计算各个本征模态在不同成分频段的能量;
在步骤B.2得到的最优参数下,原始脉搏波信号由步骤A.3优化变分模态分解得到各个本征模态,根据不同成分的能量分布情况,通过模态的能量谱密度利用公式(7)计算不同成分对应频段的能量;
步骤C.2:对运动伪影相关本征模态或频段进行去除;
通过对各个本征模态在不同成分频段能量的比较和分析,利用步骤4.B中提出的规则,即公式(8)对各个本征模态进行更新,完成运动伪影相关模态或频段的去除,具体到本实施例中,公式(8)中的阈值参数选取为θ=0.6;
公式(8)中用到的更新后的本征模态通过公式(9)计算;
步骤C.3:将步骤C.2中完成运动伪影相关模态去除的剩余模态进行求和重构,得到完成去噪的脉搏波信号;
优选的,首先对满足θ<Ek1/Ek3<1/θ条件的模态在0.5Hz和4Hz处进行拟合,通过公式(10)和公式(11)完成频域信号的不连续拟合,具体到本实施例,公式(10)和公式(11)中的Δ取值为0.1Hz,处理规则及方法与步骤4.C相同;
对更新后的本征模态进行求和重构,由公式(12)完成模态在频域中的求和重构得到去噪后的脉搏波信号,最后通过公式(13)对频域信号进行傅里叶逆变换,求得去除运动伪影的脉搏波时域信号。
至此,从步骤A到步骤C完成了本实施例一种参数自适应优化VMD的脉搏波运动伪影去除方法。
在本实施例中,运用本发明中提出的一种参数自适应优化VMD的脉搏波运动伪影去除方法。相比现有方法来说具有以下优势:
在脉搏波的变分模态分解过程中,通过考虑分解结果中本征模态的平滑程度,在优化变分模态分解目标函数中增加平滑项,平滑脉搏波分解结果中的本征模态,提高了脉搏波变分模态分解的质量;通过考虑基于本征模态间相关系数和不同本征模态中心频率离散程度,提出关联函数,完成对采集到的脉搏波变分模态分解算法的分解本征模态数K和二次惩罚因子α进行自适应选取,克服过分解和模态混叠问题,实现参数的最优选择;
在去噪过程中通过考虑本征模态能量的频域分布,计算不同成分对应频段的能量大小判断本征模态的成分属性,在不借助加速度传感器的情况下去除脉搏波信号中存在的运动伪影,克服了去除运动伪影过程对加速度传感器的依赖,降低了去噪过程中数据的繁杂程度,使实施过程更具有实时性,提高了去除脉搏波运动伪影的效率。
需要说明的是,本说明书所述的仅为本发明的较佳实施例而已,以上实施例仅用于说明本发明的技术方案而非对本发明的限制。本发明不应该局限于该实施例和附图所公开的内容,凡本领域技术人员依本发明的构思通过逻辑分析、推理或者有限的实验可以得到的技术方案,都落入本发明的范围之内。
Claims (6)
1.一种参数自适应优化VMD的脉搏波运动伪影去除方法,其特征在于:包括以下步骤:
步骤A.首先对变分模态分解的目标函数进行优化设计,加入保证本征模态平滑性质的优化项,得到优化变分模态分解的目标函数;
步骤B.通过引入拉格朗日乘子和二次惩罚因子,对步骤A中设计的目标函数进行去约束求解,通过对各个本征模态、本征模态对应的中心频率以及拉格朗日乘子进行循环更新,完成对脉搏波信号的优化变分模态分解,求得信号的本征模态{uk};
步骤C.在分解本征模态个数K和二次惩罚因子α的不同组合下,对脉搏波信号基于步骤B完成优化变分模态分解,设计参数最优选取的目标函数D,通过计算不同参数下分解结果对应的目标函数对参数进行自适应选取,得到优化变分模态分解的最优参数K和α,以及在该组参数下的分解结果,即本征模态{uk};
步骤D.在得到步骤C最优参数组合及该参数组合下的分解结果后,基于各个本征模态的能量谱完成运动伪影相关本征模态的去除,而后对剩余本征模态进行求和重构,得到完成去噪的脉搏波信号。
2.根据权利要求1所述的方法,其特征在于:步骤A所述目标函数如公式(1)所示:
其中,min为取最小值,{uk}表示优化变分模态分解结果中的本征模态集合,其元素uk表示分解结果中的本征模态,{ωk}表示优化变分模态分解结果中本征模态对应的中心频率集合,其元素ωk表示分解结果中本征模态对应的中心频率,K为分解本征模态个数;为uk的带宽估计;其中表示括号中的项对时间求梯度,δ(t)为脉冲函数,j为虚数单位,*表示卷积运算,||·||2表示二范数;增加项的目的在于保证分解结果中各个本征模态的平滑性质;其中,γ为正则加权因子,1≥γ≥0,time为所处理脉搏波信号的采样点数,Δuk(t)表示本征模态uk在t处的梯度,定义为:Δuk(t)=uk(t)-uk(t-1),uk(t)为本征模态uk在采样点t处的幅值,||·||1表示1范数;约束条件为∑uk(t)=f,其意义为分解所得所有本征模态之和等于输入信号。
3.根据权利要求1所述的方法,其特征在于:所述步骤B通过以下过程实现:
步骤B.1:初始化变分模态分解算法相关参数;
其中,U0=[u1,0,u2,0,...,ui,0,…,uK,0]T为K个本征模态的初始化,ui,0表示第i个本征模态的初始化;ω0=[ω1,0,ω2,0,...,ωi,0,...,ωK,0]T为本征模态的中心频率初始化,ωi,0表示第i个本征模态对应的中心频率初始化,λ0为拉格朗日乘子的初始化;
步骤B.2:对所求参数进行循环更新;
首先引入拉格朗日乘子和二次惩罚因子,将公式(1)所示的带约束极值问题转化为无约束极值问题,再对信号及参数进行傅里叶变换后,运用交替方向乘数法直接在频域内对无约束极值问题进行求解;经计算,求得公式(2)~(4),用于对本征模态、中心频率以及拉格朗日乘子进行循环更新:
其中,ω为时域信号经傅里叶变换后在频域中的频率;为第k个本征模态的傅里叶变换形式在第n+1次迭代更新后的结果;为输入信号的傅里叶变换形式;为第i个本征模态的傅里叶变换形式在第n+1次迭代更新后的结果;为第i个本征模态的傅里叶变换形式在第n次迭代更新后的结果;为拉格朗日乘子的傅里叶变换形式在第n次迭代更新后的结果;ωk,n为第k个本征模态的傅里叶变换形式在第n次迭代更新后的结果;α为二次惩罚因子,其作用为改善结果的收敛性;τ为拉格朗日算子的更新参数;
步骤B.3:依据工程需要,设定阈值ε,判断筛选条件是否成立,并输出计算结果;
计算筛选条件为:
其中,ε为误差允许范围参数;
若公式(5)成立,则步骤B.2中计算出的本征模态及中心频率即为输出的uk,ωk;
若公式(5)不成立,则返回步骤B.2。
4.根据权利要求1所述的方法,其特征在于:所述步骤C通过以下过程实现:
步骤C.1:初始化相关参数,并在不同的参数组合下通过步骤B完成对脉搏波信号的变分模态分解;
初始化[Kmin,Kmax,Kstep]、[αmin,αmax,αstep],其中,Kmin为本征模态分解个数K的最小取值,Kmax为本征模态分解个数的最大取值,Kstep为本征模态分解个数的寻优步长;αmin为二次惩罚因子α的最小取值、αmax为二次惩罚因子α的最大取值,αstep为二次惩罚因子α的寻优步长;
通过初始化可以得到分解本征模态个数K和二次惩罚因子α的不同组合方式,构成维的取值矩阵,通过对矩阵遍历选择不同的参数组合通过步骤B对脉搏波信号进行变分模态分解;
步骤C.2:计算不同参数下分解结果对应的如下关联函数D:
其中,用来衡量分解结果中模态间的相关性,为分解结果中模态ui和模态uj的相关系数,分别表示模态ui和模态uj的平均值,E[·]表示·的期望;用来衡量分解结果中各个模态中心频率的离散程度,为中心频率离散程度的度量;可调系数η用来调节对本征模态相关程度和中心频率离散程度的约束强度;
步骤C.3:选取最优参数;
对于在每一组特定的参数组合下完成的脉搏波信号优化变分模态分解结果,通过公式(6)计算关联函数D,直到完成所有参数组合的遍历,最终选取使关联函数D取得最小值的分解本征模态个数K和二次惩罚因子α的参数组合,并保留该参数下的分解结果{uk}。
5.根据权利要求1-4任一所述的方法,其特征在于:所述步骤D通过以下过程实现:
步骤D.1:通过公式(7)计算各个本征模态在不同成分频段的能量:
其中,Ek1为本征模态uk中脉搏波动对应频段的能量,Ek2为本征模态uk中呼吸成分对应频段的能量,Ek1和Ek2均通过对相应频段的能量谱密度对频率积分计算,为本征模态uk的能量谱密度,由于变分模态分解的求解过程是在频域中完成,故对脉搏波进行变分模态分解所得结果即为变分模态分解所得本征模态uk的频域表示Ek3为本征模态uk中运动伪影对应频段的能量,由于运动伪影成分的频率不确定,故通过本征模态总能量与其他分量能量之差计算;
步骤D.2:通过公式(8)对运动伪影相关本征模态或频段进行去除:
其中,θ为设定的固定阈值,可根据具体需要在(0,1)范围内取值,用来界定运动伪影成分与脉搏波动成分的相对大小;
为满足θ<Ek1/Ek3<1/θ条件时,更新后的本征模态,其计算方法为:只保留该模态脉搏波动成分对应的频段,去除其他频段的信号,此步骤相当于对该本征模态进行带通滤波,在频域中完成,通过公式(9)即可求得
步骤D.3:将步骤D.2中完成运动伪影相关模态去除的剩余模态进行求和重构得到去除运动伪影后的脉搏波频域信号fc(ω),对fc(ω)进行傅里叶逆变换得到去噪后的脉搏波信号。
6.根据权利要求5所述的方法,其特征在于:在完成所述步骤D.2后通过如下公式对信号不连续处进行拟合:
其中,Δ为频率增量。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910652889.4A CN110309817B (zh) | 2019-07-19 | 2019-07-19 | 一种参数自适应优化vmd的脉搏波运动伪影去除方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910652889.4A CN110309817B (zh) | 2019-07-19 | 2019-07-19 | 一种参数自适应优化vmd的脉搏波运动伪影去除方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110309817A true CN110309817A (zh) | 2019-10-08 |
CN110309817B CN110309817B (zh) | 2020-10-02 |
Family
ID=68080422
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910652889.4A Active CN110309817B (zh) | 2019-07-19 | 2019-07-19 | 一种参数自适应优化vmd的脉搏波运动伪影去除方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110309817B (zh) |
Cited By (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111189639A (zh) * | 2020-01-08 | 2020-05-22 | 重庆交通大学 | 一种基于瞬时频率优化vmd的轴承故障诊断方法 |
CN112155523A (zh) * | 2020-09-27 | 2021-01-01 | 太原理工大学 | 一种基于模态能量主成分比量化的脉搏信号特征提取与分类方法 |
CN112649196A (zh) * | 2020-11-19 | 2021-04-13 | 上海交通大学烟台信息技术研究院 | 基于频域信息的信号变分模态分解预设尺度参数选取方法 |
CN113288101A (zh) * | 2021-04-13 | 2021-08-24 | 安徽通灵仿生科技有限公司 | 基于谱相减和频域eemd-cca的运动状态下icg信号的处理方法 |
CN113657268A (zh) * | 2021-08-13 | 2021-11-16 | 江苏国电南自海吉科技有限公司 | 一种应用于风电机组齿轮箱故障诊断的信号自动分解方法 |
CN115470864A (zh) * | 2022-09-28 | 2022-12-13 | 中国人民解放军总医院 | 一种基于脉冲超宽带雷达的身份识别方法 |
CN115470645A (zh) * | 2022-09-22 | 2022-12-13 | 哈尔滨工业大学(威海) | 一种变分模态分解的参数确定方法 |
CN115590489A (zh) * | 2022-09-28 | 2023-01-13 | 中国人民解放军总医院(Cn) | 一种基于调频连续波雷达的无接触血压监测方法 |
CN116763275A (zh) * | 2023-08-22 | 2023-09-19 | 中国地质大学(武汉) | 一种用于处理光电容积脉搏信号的时频分析方法及装置 |
CN117674199A (zh) * | 2024-02-01 | 2024-03-08 | 西安热工研究院有限公司 | 一种超级电容耦合锂电池的新型电力系统调频方法和系统 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20130051674A1 (en) * | 2010-05-07 | 2013-02-28 | Bart Goossens | Method and device for estimating noise in a reconstructed image |
CN109145727A (zh) * | 2018-07-11 | 2019-01-04 | 上海电力学院 | 一种基于vmd参数优化的轴承故障特征提取方法 |
-
2019
- 2019-07-19 CN CN201910652889.4A patent/CN110309817B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20130051674A1 (en) * | 2010-05-07 | 2013-02-28 | Bart Goossens | Method and device for estimating noise in a reconstructed image |
CN109145727A (zh) * | 2018-07-11 | 2019-01-04 | 上海电力学院 | 一种基于vmd参数优化的轴承故障特征提取方法 |
Non-Patent Citations (3)
Title |
---|
周旺平等: "ABC-VMD和包络谱分析在齿轮故障诊断中的应用", 《机械传动》 * |
徐海津: "抗运动伪影下基于PPG的心率估计方法研究与应用", 《中国优秀硕士学位论文全文数据库医药卫生科技辑》 * |
苏志刚等: "基于平滑先验法去除脉搏波基线漂移", 《中国医学物理学杂志》 * |
Cited By (17)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111189639A (zh) * | 2020-01-08 | 2020-05-22 | 重庆交通大学 | 一种基于瞬时频率优化vmd的轴承故障诊断方法 |
CN111189639B (zh) * | 2020-01-08 | 2021-09-14 | 重庆交通大学 | 一种基于瞬时频率优化vmd的轴承故障诊断方法 |
CN112155523A (zh) * | 2020-09-27 | 2021-01-01 | 太原理工大学 | 一种基于模态能量主成分比量化的脉搏信号特征提取与分类方法 |
CN112649196A (zh) * | 2020-11-19 | 2021-04-13 | 上海交通大学烟台信息技术研究院 | 基于频域信息的信号变分模态分解预设尺度参数选取方法 |
CN112649196B (zh) * | 2020-11-19 | 2022-09-06 | 上海交通大学烟台信息技术研究院 | 基于频域信息的信号变分模态分解预设尺度参数选取方法 |
CN113288101A (zh) * | 2021-04-13 | 2021-08-24 | 安徽通灵仿生科技有限公司 | 基于谱相减和频域eemd-cca的运动状态下icg信号的处理方法 |
CN113288101B (zh) * | 2021-04-13 | 2023-05-26 | 安徽通灵仿生科技有限公司 | 基于谱相减和频域eemd-cca的运动状态下icg信号的处理方法 |
CN113657268A (zh) * | 2021-08-13 | 2021-11-16 | 江苏国电南自海吉科技有限公司 | 一种应用于风电机组齿轮箱故障诊断的信号自动分解方法 |
CN115470645A (zh) * | 2022-09-22 | 2022-12-13 | 哈尔滨工业大学(威海) | 一种变分模态分解的参数确定方法 |
CN115590489A (zh) * | 2022-09-28 | 2023-01-13 | 中国人民解放军总医院(Cn) | 一种基于调频连续波雷达的无接触血压监测方法 |
CN115470864B (zh) * | 2022-09-28 | 2023-05-23 | 中国人民解放军总医院 | 一种基于脉冲超宽带雷达的身份识别方法 |
CN115470864A (zh) * | 2022-09-28 | 2022-12-13 | 中国人民解放军总医院 | 一种基于脉冲超宽带雷达的身份识别方法 |
CN115590489B (zh) * | 2022-09-28 | 2024-01-23 | 中国人民解放军总医院 | 一种基于调频连续波雷达的无接触血压监测方法 |
CN116763275A (zh) * | 2023-08-22 | 2023-09-19 | 中国地质大学(武汉) | 一种用于处理光电容积脉搏信号的时频分析方法及装置 |
CN116763275B (zh) * | 2023-08-22 | 2023-10-31 | 中国地质大学(武汉) | 一种用于处理光电容积脉搏信号的时频分析方法及装置 |
CN117674199A (zh) * | 2024-02-01 | 2024-03-08 | 西安热工研究院有限公司 | 一种超级电容耦合锂电池的新型电力系统调频方法和系统 |
CN117674199B (zh) * | 2024-02-01 | 2024-04-23 | 西安热工研究院有限公司 | 一种超级电容耦合锂电池的电力系统调频方法和系统 |
Also Published As
Publication number | Publication date |
---|---|
CN110309817B (zh) | 2020-10-02 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110309817B (zh) | 一种参数自适应优化vmd的脉搏波运动伪影去除方法 | |
Singh et al. | A new ECG denoising framework using generative adversarial network | |
US9294074B2 (en) | Physiological signal denoising | |
CN113349752B (zh) | 一种基于传感融合的可穿戴设备实时心率监测方法 | |
WO2019148557A1 (zh) | 基于大脑血红蛋白信息的大脑功能状态评价装置 | |
WO2013043157A2 (en) | Physiological signal denoising | |
CN114376564A (zh) | 一种基于心冲击信号的睡眠分期方法、系统、装置及介质 | |
CN113229831B (zh) | 基于肌电和肌氧信号的运动功能监测管理方法 | |
CN111358455A (zh) | 一种多数据源的血压预测方法和装置 | |
Sharma | Heart rate extraction from PPG signals using variational mode decomposition | |
CN114569096B (zh) | 一种基于视频流的非接触式连续血压测量方法及系统 | |
Mortezaee et al. | An improved SSA-based technique for EMG removal from ECG | |
WO2022237222A1 (zh) | 基于信噪比提高技术的脉搏波获取方法、系统及存储介质 | |
Wu et al. | EMGdi signal enhancement based on ICA decomposition and wavelet transform | |
Yao et al. | A new method based CEEMDAN for removal of baseline wander and powerline interference in ECG signals | |
CN115422976B (zh) | 一种基于人工网络的心肺耦合关系分析方法及监测系统 | |
Das et al. | Non-contact heart rate measurement from facial video data using a 2d-vmd scheme | |
CN114557691A (zh) | 基于多波长的ppg信号的无创血脂检测方法及系统 | |
CN116269269A (zh) | 具有数据质量评估功能的连续血压测量仪、系统及方法 | |
CN106236075B (zh) | 一种应用于便携式心电仪所测心电图的降噪方法 | |
CN111543984A (zh) | 一种基于ssda的脑电信号的眼电伪迹去除方法 | |
CN112001862A (zh) | 消除视频心冲击信号运动噪声的非接触式视心率检测方法 | |
CN114983358A (zh) | 一种基于脉搏波的中医肺藏司呼吸功能态势评估算法 | |
CN116211308A (zh) | 一种高强度运动下机体疲劳评估方法 | |
CN115374815A (zh) | 一种基于视觉Transformer的自动睡眠分期方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |