CN110309817B - 一种参数自适应优化vmd的脉搏波运动伪影去除方法 - Google Patents

一种参数自适应优化vmd的脉搏波运动伪影去除方法 Download PDF

Info

Publication number
CN110309817B
CN110309817B CN201910652889.4A CN201910652889A CN110309817B CN 110309817 B CN110309817 B CN 110309817B CN 201910652889 A CN201910652889 A CN 201910652889A CN 110309817 B CN110309817 B CN 110309817B
Authority
CN
China
Prior art keywords
decomposition
eigenmode
pulse wave
mode
motion artifact
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.)
Active
Application number
CN201910652889.4A
Other languages
English (en)
Other versions
CN110309817A (zh
Inventor
郭树理
李亚静
戴忠健
何昆仑
韩丽娜
范利
刘宏斌
王春喜
曹丰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chinese PLA General Hospital
Beijing Institute of Technology BIT
Original Assignee
Chinese PLA General Hospital
Beijing Institute of Technology BIT
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Chinese PLA General Hospital, Beijing Institute of Technology BIT filed Critical Chinese PLA General Hospital
Priority to CN201910652889.4A priority Critical patent/CN110309817B/zh
Publication of CN110309817A publication Critical patent/CN110309817A/zh
Application granted granted Critical
Publication of CN110309817B publication Critical patent/CN110309817B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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的脉搏波运动伪影去除方法
技术领域
本发明涉及一种参数自适应优化VMD的脉搏波运动伪影去除方法，属于生物医学信号采集和处理技术领域。
背景技术
大量临床实测结果显示，脉搏波(PPG)包含丰富的心血管系统生理及病理信息，其波形特征与心血管疾病有着密切的联系。脉搏在临床中通常采用光电式脉搏传感器进行测量，其信号反映为光电容积脉搏波，是一种无创的测量方式。用位于可见红光光谱及红外光谱的两个光源交替照射检测部位(一般为耳垂或指尖)，通过反射或透射得到脉动期间的组织的光吸收量，从而获得脉搏波的基础波形、脉率估计和血氧浓度等生理信号。
脉搏波的连续监测是当前家居健康养老、远程健康监护等系统中基础的必备测量数据，可以为用户的健康监控提供心率、血氧浓度等多种生理指标，但运动伪影对脉搏波的影响一直是高效可靠地使用脉搏血氧仪进行连续实时健康监测的重大障碍。运动中的脉搏实时采集能够通过提供心率估计帮助人们及时了解自身的健康情况并应对突发状况，从而降低运动风险，但在运动中由于脉搏波采集设备与皮肤间距的不规律变化导致漏光等现象，使得检测到的脉搏波信号存在大量的运动伪影无法去除。总之，运动伪影的存在对脉搏波的后续分析产生了不可避免的影响。
最初的运动伪影去除方法包括独立成分分析法、自适应滤波技术等算法，但这些方法普遍存在以下问题：第一，两种算法都建立在运动伪影与脉搏波信号完全不相关的基本假设基础之上，这一假设与真实状况并不相符；第二，当运动伪影与脉搏波信号发生本征模态混叠时，算法的去噪效果变差。在随后的发展中，为解决上述问题，经验本征模态分解(EMD)和集合经验本征模态分解(EEMD)等方法被应用于脉搏波信号，通过将信号分解为几个固有本征模态分量将运动伪影与脉搏波分离，达到去除运动伪影的目的，但算法缺少数学基础，可解释性差，并且在信号中存在较强运动伪影的情况下失效。
变分模态分解(VMD)是现有算法中较好的运动伪影去除算法，该算法设置合适的目标函数，通过对带约束求极值问题的变分求解，将信号分解为若干个本征模态，并在此基础上完成信号的去噪过程。但该算法现在还存在以下问题：第一，VMD算法以各个本征模态带宽最窄为目标，但其目标函数的设定没有充分考虑到分解结果中本征模态的平滑程度；第二，VMD算法中的分解本征模态数K、二次惩罚因子α一般通过经验进行提前设置，无法确认算法的最优性；第三，现有对运动伪影的去除方法一般依靠加速度信号的辅助，去噪过程需要同时采集运动信号，导致涉及到的数据量较多。
本发明针对上述存在的缺陷，提出一种参数自适应优化变分模态分解技术，致力于提高脉搏波变分模态分解的质量并简化去除运动伪影过程中的数据复杂度。
发明内容
针对现有技术中存在的不足，本发明的目的在于提出一种参数自适应优化变分模态分解技术，以完成对脉搏波信号的分解并去除运动伪影。在参数自适应选取部分，通过对变分模态分解结果的分析，本发明提出采用各个本征模态间的相关系数作为指标函数，评价分解结果的好坏，从而克服本征模态混叠问题，实现参数的自适应；在变分模态分解的具体算法实施中，本发明对目标函数进行改进，考虑变分模态的平滑程度，提出新的优化目标函数，从而使算法针对脉搏波具有更好的效果；在运动伪影去除过程中，本发明提出基于能量分布的去噪方法，仅依靠脉搏波数据完成去噪，在不影响运动伪影去除效果的同时，相比传统方法中借助加速度相关数据去噪的方法，降低了数据的繁杂程度。
本发明旨在解决脉搏波信号中运动伪影去除的问题，提出了一种参数自适应优化VMD的脉搏波运动伪影去除方法。
一种参数自适应优化VMD的脉搏波运动伪影去除方法，包括如下步骤：
步骤1：首先对变分模态分解的目标函数进行优化设计，加入保证本征模态平滑性质的优化项，得到优化变分模态分解的目标函数；
作为优选，步骤1通过以下过程实现：
变分模态分解算法是一种以变分问题为框架的信号分解估计方法，该算法以本征模态的带宽估计函数取得最小为目标函数对输入信号进行最优求解，考虑到分解结果中本征模态的平滑程度，本发明首先设计如下的目标函数：
Figure GDA0002634852720000021
s.t.∑uk(t)＝f
其中，min为取最小值，{uk}表示优化变分模态分解结果中的本征模态集合，其元素uk表示分解结果中的本征模态，{ωk}表示优化变分模态分解结果中本征模态对应的中心频率集合，其元素ωk表示分解结果中本征模态对应的中心频率，K为分解本征模态个数；
Figure GDA0002634852720000031
为uk的带宽估计；其中
Figure GDA0002634852720000032
表示括号中的项对时间求梯度，δ(t)为脉冲函数，j为虚数单位，*表示卷积运算，||·||2表示二范数；
增加
Figure GDA0002634852720000033
项的目的在于保证分解结果中各个本征模态的平滑性质；γ(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)，用于对本征模态、中心频率以及拉格朗日乘子进行循环更新：
Figure GDA0002634852720000041
Figure GDA0002634852720000042
Figure GDA0002634852720000043
其中，ω为时域信号经傅里叶变换后在频域中的频率；
Figure GDA0002634852720000044
为第k个本征模态的傅里叶变换形式在第n、n+1次迭代更新后的结果；
Figure GDA0002634852720000045
为输入信号的傅里叶变换形式；
Figure GDA0002634852720000046
为第i个本征模态的傅里叶变换形式在第n+1次迭代更新后的结果；
Figure GDA0002634852720000047
为第i个本征模态的傅里叶变换形式在第n次迭代更新后的结果；
Figure GDA0002634852720000048
为拉格朗日乘子的傅里叶变换形式在第n、n+1次迭代更新后的结果；ωk,n、ωk,n+1为第k个本征模态的中心频率在第n、n+1次迭代更新后的结果；α为二次惩罚因子，其作用为改善结果的收敛性；τ为拉格朗日算子的更新参数；
步骤2.C：依据工程需要，设定阈值ε，判断筛选条件是否成立，并输出计算结果；
计算筛选条件为：
Figure GDA0002634852720000049
其中，ε为误差允许范围参数；
Figure GDA00026348527200000410
为第k个本征模态的傅里叶变换形式在第n次迭代更新后的结果；
若公式(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和二次惩罚因子α的不同组合方式，构成
Figure GDA0002634852720000051
维的取值矩阵，通过对矩阵遍历选择不同的参数组合通过步骤2对脉搏波信号进行变分模态分解；
步骤3.B：计算不同参数下分解结果对应的关联函数D；
在完成对脉搏波变分模态分解的过程中，可能存在欠前分解导致的模态混叠问题或过分解问题，从变分模态分解的结果考虑，期望分解结果中各个本征模态的相关程度达到最小，同时期望各个本征模态中心频率的离散程度达到最大，本发明综合考虑到以上两点要求，提出参数最优选取的关联函数如下：
Figure GDA0002634852720000052
其中，第一项中
Figure GDA0002634852720000053
为各个模态间相关系数的平均值，用来衡量分解结果中模态间的相关性，其取值范围为[0,1]，
Figure GDA0002634852720000054
为分解结果中模态ui和模态uj的相关系数，
Figure GDA0002634852720000055
分别表示模态ui和模态uj的平均值，E[·]表示·的期望；
第二项中
Figure GDA0002634852720000056
用来衡量分解结果中各个模态中心频率的离散程度，
Figure GDA0002634852720000057
为中心频率离散程度的度量，加入负指数的目的是保证模态相关程度和中心频率分散程度对最优解取极值的一致性，即当关联函数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左右的低频分量或更高的频段；根据不同成分的能量分布情况，将每个本征模态的能量分为三个部分并利用模态的能量谱密度计算相应频段的能量：
Figure GDA0002634852720000061
其中Ek1为本征模态uk中脉搏波动对应频段的能量，Ek2为本征模态uk中呼吸成分对应频段的能量，Ek1和Ek2均通过对相应频段的能量谱密度对频率积分计算，
Figure GDA0002634852720000071
为本征模态uk的能量谱密度，由于变分模态分解的求解过程在频域中完成，故对脉搏波进行变分模态分解所得结果即为变分模态分解所得本征模态uk的频域表示
Figure GDA0002634852720000072
Ek3为本征模态uk中运动伪影对应频段的能量，由于运动伪影成分的频率不确定，故通过本征模态总能量与其他分量能量之差计算；
步骤4.B：对运动伪影相关本征模态或频段进行去除；
通过对各个本征模态在不同成分频段能量的比较和分析，利用以下规则对各个本征模态进行更新，完成运动伪影相关模态或频段的去除：
Figure GDA0002634852720000073
其中，θ为设定的固定阈值，可根据具体需要在(0,1)范围内取值，用来界定运动伪影成分与脉搏波动成分的相对大小；uk(ω)为变分模态分解所得本征模态，u'k(ω)为完成运动伪影相关模态或频段的去除的各个本征模态；
对公式(8)具体描述如下：
当运动伪影成分能量远大于脉搏波动成分能量时，即对应公式(8)中的Ek1/Ek3≤θ，认为该本征模态的主要构成成分为运动伪影，将该本征模态置零；当脉搏波动成分能量与运动伪影成分能量相差不大时，即对应公式(8)中的θ＜Ek1/Ek3＜1/θ，认为该本征模态中同时包含脉搏波有效成分及运动伪影成分，对该本征模态进行处理更新为
Figure GDA0002634852720000074
当脉搏波动成分能量远大于运动伪影成分能量时，即对应公式(8)中的Ek1/Ek3≥1/θ，认为该本征模态的主要构成成分为脉搏波有效信息，保留该本征模态；
Figure GDA0002634852720000075
为满足θ＜Ek1/Ek3＜1/θ条件时，更新后的本征模态，其计算方法为：只保留该本征模态的脉搏波动成分对应的频段，去除其他频段的信号，此步骤相当于对该本征模态进行带通滤波，在频域中完成，通过公式(9)即可求得
Figure GDA0002634852720000076
Figure GDA0002634852720000081
步骤4.C：将步骤4.B中完成运动伪影相关模态去除的剩余模态进行求和重构，得到完成去噪的脉搏波信号；
由于在去噪过程中对既含有运动伪影成分又含有脉搏波动成分的模态，即满足θ＜Ek1/Ek3＜1/θ条件的模态进行了带通滤波更新，造成其在0.5Hz和4Hz处的阶跃性质，故在重构前考虑这些本征模态在该频率处的连续性，对信号不连续处进行拟合，通过公式(10)和公式(11)完成：
Figure GDA0002634852720000082
Figure GDA0002634852720000083
其中，Δ为取值较小的频率增量，用来在不连续点前后提供拟合区间，可根据实际情况进行选取；公式(10)和公式(11)的意义为：在不连续点处通过区间端点的信号取值将该点前后小区间内的信号拟合为直线，以保证信号的连续性；
对更新后的本征模态进行求和重构，由公式(12)完成模态在频域中的求和重构：
Figure GDA0002634852720000084
其中，u'k(ω)为步骤4.B中通过公式(8)以及本步骤中公式(10)～(11)完成更新和不连续点拟合的本征模态，fc(ω)为去除运动伪影后的脉搏波频域信号；
最后通过公式(13)对频域重构信号进行傅里叶逆变换，即可得到去除噪声及运动伪影的脉搏波信号；
fc(t)＝Φ-1[fc(ω)] (13)
其中，fc(t)为去除运动伪影后的脉搏波信号，Φ-1[·]表示·的傅里叶逆变换，本步骤通过对去除运动伪影相关本征模态或频段后更新的本征模态进行求和运算得到去噪后的脉搏波信号。
至此，由步骤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)中用到的更新后的本征模态
Figure GDA0002634852720000121
通过公式(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 (5)

1.一种参数自适应优化VMD的脉搏波运动伪影去除方法，其特征在于：包括以下步骤：
步骤A.首先对变分模态分解的目标函数进行优化设计，加入保证本征模态平滑性质的优化项，得到优化变分模态分解的目标函数；
步骤B.通过引入拉格朗日乘子和二次惩罚因子，对步骤A中设计的目标函数进行去约束求解，通过对各个本征模态、本征模态对应的中心频率以及拉格朗日乘子进行循环更新，完成对脉搏波信号的优化变分模态分解，求得信号的本征模态{uk}；
步骤C.在分解本征模态个数K和二次惩罚因子α的不同组合下，对脉搏波信号基于步骤B完成优化变分模态分解，设计参数最优选取的关联函数D，通过计算不同参数下分解结果对应的目标函数对参数进行自适应选取，得到优化变分模态分解的最优参数K和α，以及在所述最优参数K和α下的分解结果，即本征模态{uk}；
步骤D.在得到步骤C最优参数组合及该参数组合下的分解结果后，基于各个本征模态的能量谱完成运动伪影相关本征模态的去除，而后对剩余本征模态进行求和重构，得到完成去噪的脉搏波信号；
步骤A所述目标函数如公式(1)所示：
Figure FDA0002634852710000011
s.t.∑uk(t)＝f
其中，min为取最小值，{uk}表示优化变分模态分解结果中的本征模态集合，其元素uk表示分解结果中的本征模态，{ωk}表示优化变分模态分解结果中本征模态对应的中心频率集合，其元素ωk表示分解结果中本征模态对应的中心频率，K为分解本征模态个数；
Figure FDA0002634852710000012
为uk的带宽估计；其中
Figure FDA0002634852710000013
表示括号中的项对时间求梯度，δ(t)为脉冲函数，j为虚数单位，*表示卷积运算，||·||2表示二范数；增加
Figure FDA0002634852710000014
项的目的在于保证分解结果中各个本征模态的平滑性质；其中，γ为正则加权因子，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所述的方法，其特征在于：所述步骤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)，用于对本征模态、中心频率以及拉格朗日乘子进行循环更新：
Figure FDA0002634852710000021
Figure FDA0002634852710000022
Figure FDA0002634852710000023
其中，ω为时域信号经傅里叶变换后在频域中的频率；
Figure FDA0002634852710000024
为第k个本征模态的傅里叶变换形式在第n、n+1次迭代更新后的结果；
Figure FDA0002634852710000025
为输入信号的傅里叶变换形式；
Figure FDA0002634852710000026
为第i个本征模态的傅里叶变换形式在第n+1次迭代更新后的结果；
Figure FDA0002634852710000027
为第i个本征模态的傅里叶变换形式在第n次迭代更新后的结果；
Figure FDA0002634852710000028
为拉格朗日乘子的傅里叶变换形式在第n、n+1次迭代更新后的结果；ωk,n、ωk,n+1为第k个本征模态的傅里叶变换形式在第n、n+1次迭代更新后的结果；α为二次惩罚因子，其作用为改善结果的收敛性；τ为拉格朗日算子的更新参数；
步骤B.3：依据工程需要，设定阈值ε，判断筛选条件是否成立，并输出计算结果；
计算筛选条件为：
Figure FDA0002634852710000029
其中，ε为误差允许范围参数；
若公式(5)成立，则步骤B.2中计算出的本征模态及中心频率即为输出的uk,ωk；
若公式(5)不成立，则返回步骤B.2。
3.根据权利要求1所述的方法，其特征在于：所述步骤C通过以下过程实现：
步骤C.1：初始化相关参数，并在不同的参数组合下通过步骤B完成对脉搏波信号的变分模态分解；
初始化[Kmin,Kmax,Kstep]、[αmin,αmax,αstep]，其中，Kmin为本征模态分解个数K的最小取值，Kmax为本征模态分解个数的最大取值，Kstep为本征模态分解个数的寻优步长；αmin为二次惩罚因子α的最小取值、αmax为二次惩罚因子α的最大取值，αstep为二次惩罚因子α的寻优步长；
通过初始化可以得到分解本征模态个数K和二次惩罚因子α的不同组合方式，构成
Figure FDA0002634852710000031
维的取值矩阵，通过对矩阵遍历选择不同的参数组合通过步骤B对脉搏波信号进行变分模态分解；
步骤C.2：计算不同参数下分解结果对应的如下关联函数D：
Figure FDA0002634852710000032
其中，
Figure FDA0002634852710000033
用来衡量分解结果中模态间的相关性，
Figure FDA0002634852710000034
为分解结果中模态ui和模态uj的相关系数，
Figure FDA0002634852710000035
分别表示模态ui和模态uj的平均值，E[·]表示·的期望；
Figure FDA0002634852710000036
用来衡量分解结果中各个模态中心频率的离散程度，
Figure FDA0002634852710000037
为中心频率离散程度的度量；可调系数η用来调节对本征模态相关程度和中心频率离散程度的约束强度；ωi、ωi+1为模态ui、ui+1的中心频率；
步骤C.3：选取最优参数；
对于在每一组特定的参数组合下完成的脉搏波信号优化变分模态分解结果，通过公式(6)计算关联函数D，直到完成所有参数组合的遍历，最终选取使关联函数D取得最小值的分解本征模态个数K和二次惩罚因子α的参数组合，并保留该参数组合下的分解结果{uk}。
4.根据权利要求1-3任一所述的方法，其特征在于：所述步骤D通过以下过程实现：
步骤D.1：通过公式(7)计算各个本征模态在不同成分频段的能量：
Figure FDA0002634852710000041
其中，Ek1为本征模态uk中脉搏波动对应频段的能量，Ek2为本征模态uk中呼吸成分对应频段的能量，Ek1和Ek2均通过对相应频段的能量谱密度对频率积分计算，
Figure FDA0002634852710000042
为本征模态uk的能量谱密度，由于变分模态分解的求解过程是在频域中完成，故对脉搏波进行变分模态分解所得结果即为变分模态分解所得本征模态uk的频域表示
Figure FDA0002634852710000043
Ek3为本征模态uk中运动伪影对应频段的能量，由于运动伪影成分的频率不确定，故通过本征模态总能量与其他分量能量之差计算；
步骤D.2：通过公式(8)对运动伪影相关本征模态或频段进行去除：
Figure FDA0002634852710000044
其中，θ为设定的固定阈值，可根据具体需要在(0,1)范围内取值，用来界定运动伪影成分与脉搏波动成分的相对大小；uk(ω)为变分模态分解所得本征模态，u'k(ω)为完成运动伪影相关模态或频段的去除的各个本征模态，ω为本征模态的中心频率；
Figure FDA0002634852710000045
为满足θ＜Ek1/Ek3＜1/θ条件时，更新后的本征模态，其计算方法为：只保留该本征模态的脉搏波动成分对应的频段，去除其他频段的信号，此步骤相当于对该本征模态进行带通滤波，在频域中完成，通过公式(9)即可求得
Figure FDA0002634852710000046
Figure FDA0002634852710000051
步骤D.3：将步骤D.2中完成运动伪影相关模态去除的剩余模态进行求和重构得到去除运动伪影后的脉搏波频域信号fc(ω)，对fc(ω)进行傅里叶逆变换得到去噪后的脉搏波信号。
5.根据权利要求4所述的方法，其特征在于：在完成所述步骤D.2后通过如下公式对信号不连续处进行拟合：
Figure FDA0002634852710000052
Figure FDA0002634852710000053
其中，Δ为频率增量。
CN201910652889.4A 2019-07-19 2019-07-19 一种参数自适应优化vmd的脉搏波运动伪影去除方法 Active CN110309817B (zh)

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 CN110309817A (zh) 2019-10-08
CN110309817B true 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)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111189639B (zh) * 2020-01-08 2021-09-14 重庆交通大学 一种基于瞬时频率优化vmd的轴承故障诊断方法
CN112155523B (zh) * 2020-09-27 2022-09-16 太原理工大学 一种基于模态能量主成分比量化的脉搏信号特征提取与分类方法
CN112649196B (zh) * 2020-11-19 2022-09-06 上海交通大学烟台信息技术研究院 基于频域信息的信号变分模态分解预设尺度参数选取方法
CN113288101B (zh) * 2021-04-13 2023-05-26 安徽通灵仿生科技有限公司 基于谱相减和频域eemd-cca的运动状态下icg信号的处理方法
CN113657268B (zh) * 2021-08-13 2023-01-31 江苏国电南自海吉科技有限公司 一种应用于风电机组齿轮箱故障诊断的信号自动分解方法
CN115470864B (zh) * 2022-09-28 2023-05-23 中国人民解放军总医院 一种基于脉冲超宽带雷达的身份识别方法
CN115590489B (zh) * 2022-09-28 2024-01-23 中国人民解放军总医院 一种基于调频连续波雷达的无接触血压监测方法
CN116763275B (zh) * 2023-08-22 2023-10-31 中国地质大学（武汉） 一种用于处理光电容积脉搏信号的时频分析方法及装置
CN117674199B (zh) * 2024-02-01 2024-04-23 西安热工研究院有限公司 一种超级电容耦合锂电池的电力系统调频方法和系统

Citations (2)

* Cited by examiner, † Cited by third party
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参数优化的轴承故障特征提取方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
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 (2)

* Cited by examiner, † Cited by third party
Title
基于平滑先验法去除脉搏波基线漂移;苏志刚等;《中国医学物理学杂志》;20181031;第35卷(第10期);第1197-1202页 *
抗运动伪影下基于PPG的心率估计方法研究与应用;徐海津;《中国优秀硕士学位论文全文数据库医药卫生科技辑》;20180915(第09期);第15-16页、第30-33页 *

Also Published As

Publication number Publication date
CN110309817A (zh) 2019-10-08

Similar Documents

Publication Publication Date Title
CN110309817B (zh) 一种参数自适应优化vmd的脉搏波运动伪影去除方法
US9294074B2 (en) Physiological signal denoising
CN111358455B (zh) 一种多数据源的血压预测方法和装置
CN110680308B (zh) 基于改进emd与阈值法融合的心电信号去噪方法
WO2019148557A1 (zh) 基于大脑血红蛋白信息的大脑功能状态评价装置
CN113229831B (zh) 基于肌电和肌氧信号的运动功能监测管理方法
Motin et al. PPG derived heart rate estimation during intensive physical exercise
Sharma Heart rate extraction from PPG signals using variational mode decomposition
CN111310570A (zh) 一种基于vmd和wpd的脑电信号情感识别方法及系统
WO2022237222A1 (zh) 基于信噪比提高技术的脉搏波获取方法、系统及存储介质
CN112263253B (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
CN111543984A (zh) 一种基于ssda的脑电信号的眼电伪迹去除方法
CN111493821B (zh) 一种基于modwt及中值滤波的ppg信号实时去噪方法
CN117357080A (zh) 近红外光谱信号去噪方法及装置、终端设备、存储介质
CN114557691B (zh) 基于多波长的ppg信号的无创血脂检测方法及系统
CN116211308A (zh) 一种高强度运动下机体疲劳评估方法
CN114569096B (zh) 一种基于视频流的非接触式连续血压测量方法及系统
Klein Optimizing real-time fNIRS in BCI and neurofeedback: A comprehensive overview of strategies to improve reliability, spatial specificity, and signal quality
Mohamed et al. Diagnosis of Cardiac Disorders in ECG Signal using Artificial Neural Network Algorithm
CN114983358A (zh) 一种基于脉搏波的中医肺藏司呼吸功能态势评估算法
CN114569096A (zh) 一种基于视频流的非接触式连续血压测量方法及系统
CN116269285B (zh) 一种非接触式常态化心率变异性估计系统
CN115422976B (zh) 一种基于人工网络的心肺耦合关系分析方法及监测系统

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