一种基于稀疏贝叶斯分解的抗运动干扰心率提取算法

一种基于稀疏贝叶斯分解的抗运动干扰心率提取算法

一种基于稀疏贝叶斯分解的抗运动干扰心率提取算法王群张珣刘志文北京理工大学电子信息工程学院目的提出一种棊于光电容积脉搏波(photoplcthsmography,PPG)的抗运动干扰心率检测算法。方法算法利用同步记录的一路加速度信号频谱调制离散傅里叶变换(discretefouriertransform,DET)稀疏基，在此基上对W路PPG信号进行稀疏贝叶斯分解。由于谱系数的同一稀疏约束，此方法能够识别并去除PPG信号谱中运动干扰产生的谱峰。最后从五类特征点中搜索心率谱峰。结果实验基于12组PPG数据，算法结果平均绝对误差为1.99次/分。结论本文方法在基于PPG信号的心率检测方面效果良好，具有准确度高，抗干扰性强的特点。关键词：心率;光电容积脉搏波;稀疏14叶斯分解;加速度信号;positionWangQunZhangXunLiuZhiwenSchoolofElectronicsandInformationEngineering，BeijingInstituteofTechnology;Abstract：ObjectiveToproposeamethodforheartratemonitoringduringphysicalactivitiesusingphotoplethysinography(PPG).MethodsThespectrumsof2channelPPGsignalswereestimatedbasedonaDiscreteFourierTransform(DFT)sparsebasis,whichwasmodulatedbythespectrumofone-channelsimultaneouslyrecordedaccelerationsignal.Withacommonsparsityconstraintonthespectralcoefficients,themethodcouldeasilyidentifyandremovethespectralpeaksofmotionartifact(MA)inthePPGspectrum.Finally,thespectralpeak\ncorrespondedtoheartratewassearchedfromfivekindsoffeaturepoints.ResultsExperimentswereconductedon12groupsofPPGsignals.Theaverageabsoluteestimationerrorwasi.99beats/minandthestandarddeviationwas2.44beats/min.ConclusionThemethodisofagoodperformanceinPPG-basedheartratemonitoring.Theaccuracyandthemotion-resistantabilityofthismethodaresatisfactory.Keyword：heartrate;photoplethysmography(PPG);sparseBayesiandecomposition;accelerationsignal;心率(heartrate,HR)可以反映心脏的功能状况，在人体运动过程中可作为一种衡量运动强度的指标。目前光电容积脉搏波(photoplethysmography,PPG)广泛用于心率监测，相较于电极检测心电信号(electrocardiogram,ECG)所带来的不适感，PPG的引入使心率的无创测量更加简易可行［1-2］。PPG信号描述血液容积变化，与心率基木具有一致性拉1。PPG信号虽然蕴含心率信息，但对运动干扰(motionartifact,MA)十分敏感。在运动干扰影响下，心率测量结果的可靠性下降。因此从存在运动干扰的PPG信号屮准确提取心率已成为生物医学领域的研宄热点［4-5］。典型的心率估计算法一般是将运动伪差从受干扰的PPG信号中滤除。目前己有一些文献初步实现了从PPG信号中抑制运动干扰的目的M。然而这些技术大多数针对临床应用，可以解决受试者测试过程中产生的微小运动，例如手指抖动和行走［7-9］。这种情况下运动干扰并不严重，可以获得较好的测量结果。当运动干扰剧烈时，这些算法的处理结果并不理想。其中较常用的是周期阁算法，文献［10］表明原始PPG信号在利用该方法进行谱估计时，由于频谱泄漏，心率对应的谱峰不能和摆臂造成的相近峰值相区分，这一现象促进高分辨率线性谱估计的应用。近年来，基于稀疏信号分解(sparsesignaldecomposition,SSD)的谱估计较传统的谱估计表现出较多的优势mi，该类谱估计有较高的谱分辨率，且具有鲁棒性。加速度信号能够反映人体运动信息，提供运动干扰的衡量标准。目前已有很多心率提取方法通过加速度信号的频谱特征，去除PPG信号中运动干扰，如文献［12］用加速度信号作构建卡尔曼滤波的参考信号。但由于运动的不规律性以及复杂性,使得仅依靠加速度信号进行滤波所获取心率的结果并不理想。在对PPG信号处理获取心率的研宄中，张智林等［10］提出了TR0TKA方法，该方法主要包括3个步骤:奇异值分解，用于谱估计的稀疏信号分解，带有校正的谱峰追踪。奇异值分解用以抑制原始PPG信号的部分运动干扰，但对信号稀疏度要求不高的稀疏算法，这一步可以省去以降低算法复杂度。TR01KA对一路PPG信号进行谱估计时，用稀疏信号分解以克服传统的能量谱估计存在谱叠加从而分辨率低等缺点，而且单观测向量模型可扩展为多观测向量模型以提高信号分解准确度［13-14］。在谱峰追踪方面，利用一路分解谱联合加速度信号容易出现\n谱峰追踪丢失现象，需要合理引入附加元素以增加谱峰搜索的准确度。本文对2路PPG信号进行同时分解，接着进行谱峰追踪，在稀疏解(即分解频谱)的一定范围内搜索谱峰，结合PPG的频谱、分解频谱以及归一化加速度信号的频谱，筛选心率可能出现位置，对其作最优选择。最后，将本文提出的方法应用于受强烈运动干扰影响的PPG信号，结果表现出较强的抗干扰性和较高的心率估计准确性。1方法木文提出的心率检测算法主要括2部分:基于改进稀疏基的信号稀疏分解和改进的谱峰追踪。方法流程如图1所示。信号预处理中需要对PPG信号和加速度信号进行0.4~5.0llz的带通滤波［10］，去除噪声以及运动干扰。图1本文方法流程图Fig.1Flowchartofproposedmethod1.1联合稀疏信号分解1.1.1多观测向量模型联合稀疏分解相较FFT有更高分辨率，分解出FFT不能识别的峰值，这对寻找PPG的谱峰从而确定心率尤为重要。联合稀疏分解所需的多观测向量(multiplemeasurementvectors,MNIVs)模型是单观测向量(singlemeasurementvector,SMV)模型的扩展，模型可表示为其中是由L个观测列向量构成，每个观测向量有N个元素。O为稀疏基，是未知源矩阵（或解矩阵），每一行代表1个源，VER表示观测噪声或者误差矩阵。以丽V模型为基础，许多稀疏分解算法己被提出11^1。木文采用T-MSBL算法，因为它在矩阵O的相邻列高度相关下，具有计算快速以及准确可靠的表现，同时能分解非稀疏信号。1.1.2基于T-MSBL的联合稀疏分解算法T-MSBL算法是基于一种块稀疏贝叶斯的学习（blocksparseBayesianlearning,BSBL）构架［15KBSBL构架假设所冇源（X的第i行）是互相独立的，且\n服从高斯分布，可表示为其中，丫：是非负超参数，用来控制X的行稀疏。当丫F0，对应的为零。Bi是正定矩阵，用以描述的相关结构。通过y=vec（Y）ER，D=OIL,x=vec（X）^R,v=vec（V）,将丽V模型转换成块SMV模型，表示为其中，vec（•）为矩阵向量化，（•）为矩阵转罝，表示2个矩阵的Kronecker积。假设噪声向量v的元素彼此独立，并且服从高斯分布，即p（vj~N（0，入），Vi是v的第i个元素，A是方差。已知X，则y的高斯似然为\n贝叶斯准则，依此有则X的最大后验估计为为避免过拟合，将B,统一为B。则表示矩阵对角元素为a,，…，aM。对BSBL架构中的参数Yi、B、A进行最大期望估计，得到T-SBL算法［16KT-SBL算法在分解上表现很好，但因为矩阵増大导致其运行速度缓慢。通过适当近似使计算空间降维，此时导出算法即为T-MSBL。T-MSBL在同T-SBL一样有效基础上,有着较低的计算复杂度。参数估计如下：其屮，＜1^是①的第i列，II•IIb•表示矩阵的Frobenius规范，Tr(•)表示矩阵的迹。1.2改进的稀疏基构建稀疏分解选用DFT作为稀疏基O。，并利用加速度信号的频谱对进行调制，在此为方便寻峰对加速度信号频谱进行归一化。其屮，DFT稀疏基表示为\n稀疏基调制可表示为\n其中Acc为归一化的加速度信号频谱，Ace」表示Acc的第j个元素，0。.』，＜!）.」•分别表示O。、O的第j列。由矩阵表达式可知，在观测矩阵不变情况下，稀疏棊的原子越大，对应的分解谱项越小。那么利用上述调制稀疏基对PPG信号进行分解，可使分解后PPG信号的频谱幅值在加速度信号频谱较大处得到有效抑制。因此，分解后频谱的较小幅度可认为是运动干扰产生的，若令分解系数较小项归零，可完成信号的运动干扰去除。图2为PPG信号分别在调制的DFT棊和DFT棊下分解结果，由图可知在调制DFT基下分解的结果能明显突出心率峰，相比在简单的DFT基下的分解结果更能满足心率峰值筛选需求。图2两基分解结果对比Fig.2Comparisonofdecompositionresultsoftwobasis1.3改进的谱峰搜索方法木文采用的谱峰追踪方法依次分为初始化、峰值搜索、峰值选取、心率峰值选取4个部分。原理为己知如果心率计算间隔吋间足够短，即2次吋间窗U重叠部分很多，连续2次测得的心率值变化不大则上一次心率的测量结果可成为此次谱峰搜索的中心，从而确定本次心率范围。由于PPG信号中心率成分在短时间窗内是近似周期的，而运动干扰一般是非周期的。因此心率的基波谱峰明显，随机谱峰较小。2路PPG信号的联合谱更能凸显这一现象。1.3.1初始化与峰值搜索1）初始化在进行心率搜索前需获取心率初始值，在此阶段受试者要求尽量减少手部动作。因此心率可以通过查找此段时间PPG信号频谱的最高谱峰获得。2）峰值搜索初始化后对心率进行估计，首先对峰值进行搜索，一般情况下只记录分解频谱的峰值位置即可，但考虑到存在加速度信号与PPG信号受其干扰的非一致性，导致分解谱心率峰没有突显，故同吋记录这些位置对应的频谱，距离上一心率位置集合，以及归一化的加速度信号频谱。流程如图3所示。图3峰值搜索流程图Fig.3Flowchartofpeaksearch\nprvLoc表示上次估计心率对应的峰值位罝，prvPBM表示对应的心率值。以上次心率位置prvLoc为中心，本次在分解频谱搜索范围为，其中△，为1个正整数。考虑到当前心率变化可能较大，故在内同时进行搜索峰值，其中A2*l个正整数，KA2>Ab即仏域包含R：域。若在分解频谱的R：范围内有谱峰L对应加速度信号频谱幅值过高，则将此谱峰从搜索域屮去除。若R:屮没有出现谱峰，则此次心率保留上次估计结果。为方便表达，记频谱F，分解频谱H，相对上一心率位置集合E，在上一心率位罝左侧时E为负，右侧时为正，J表示归一化的加速度信号频谱。curLoc表示当前估计心率对应谱峰的位置，curBPM表示对应的心率值。curLoc与curBPM的转化关系可表示为1.3.2峰值选取从搜索得到的分解谱峰中确定心率对应峰的位置。木文从五类特征点中选取:距上一心率位置最近点LMIN，R:屮谱峰对应频谱幅值最大的峰位置，分解谱峰最高峰的位置，R2中谱峰对应频谱幅值最大峰的位置，分解谱峰最高峰的位置。下面依次讨论依据加速度信号频谱判断PPG信号受干扰情况以及多特征点的选取。1)运动干扰强度判别。在利用加速度信号归一化频谱判断PPG谱峰受运动干扰影响程度时，考虑2种情况:一是当所选谱峰位置对应的加速度频谱幅值小于阈值Thb认为该位置谱峰受运动干扰小。二是当所选PPG谱峰位置对应的加速度频谱幅值大于阈值Th2,且大于次优峰对应的加速度谱峰幅值的2倍，认为此峰受运动干扰较大，应被次优峰替代。因此步重复使用，以AP(L,|L2)代指该过程，如图4所示，L为初步选取的最优峰位置，1、2为次优峰位置。AP(LjL2)=1表示最优峰发生替换，为0时表示最优峰保留。图4AP流程图Fig.4FlowchartofAP\n2）五特征点选取。在确定点时Lmin,有可能出现分布于上一心率两侧的等距点LMIN1,L胃2。此时设分解频谱幅值较大点L胃为最优点，分解频谱幅值较小点为LmIN2次优点，若在加速度频谱判断下不发生替代，选择LMIN1作为Lmin,否则选择LMIN2作为Lmin。分别初步选峰值在分解频谱、频谱幅值最大位置作为最优点，两谱的次高峰位置作为对应的次优点，依据AP判断是否发生替代，最终确定。流程如图5所示，多特征点可尽量列出心率可能峰值位罝。1.3.3心率峰值选取利用选取的五待征点U。，对心率峰的位置点进行综合判断，依次判断3种情况，以避免心率追踪丢失现象。1）若位于上一心率位置同侧，丑对应峰值分别大于峰值的2倍，并令分别作为相应的次优峰位置，通过加速度信号频谱判断，若不被替换，则当前心率频率位置为的均值。图5五类特征点选取流程图Fig.5Flowchartoffivekindsoffeaturepointsselection2）若位于上一心率位置两侧，则当前心率位置为的均值，以保证下次寻峰屮心偏离真值不大。3）对位于上一心率同侧情况进行判断，若对应加速度谱峰小于\nThb则选取该点作为当前心率位置;若对应加速度谱峰大于Thb本文从，U中选取心率位置，在因加速度谱峰较大发生更替或与相距大于A2时，选取作为心率位置，否则将选作当前心率位置。心率位置选取流程如图6所示。2结果与讨论2.1数据源数据源来自TROIKA所使用的12组实验数据［10］。每组数据由同步采集的1路ECG信号、2路PPG信号以及3路加速度信号组成。12位受试者为年龄在18~35岁间的健康男性。2路PPG信号分别从2部发射绿光（波长609™）的PPG信号采集装置屮获得。加速度信号从三轴加速度计获取。PPG信号采集装置和加速度计嵌入到同一个腕带内，方便穿戴。ECG信号在胸部测得。6路信号的采样频率均为125Hz。数据采集期间，要求受试者完成变速跑或走，即最初尽量减少手部动作静止2〜3s,接着依次以1〜2km/h的速度跑0.5min,再以6〜8km/h的速度跑1min,再以12〜15km/h的速度跑lmin,最后以1〜2km/h的速度跑0.5min。在运动期间要求受试者用佩戴腕带的手进行擦汗等动作，增强随机运动干扰对PPG信号的影响。图6心率位置选取流程图Fig.6Flowchartoflocationselectionofheartrate表1本文算法与TROIKA算法检测心率的平均绝对误差Table1ComparisonofproposedalgorithmandTROIKAintermsofaverageabsoluteerrors(Error1)ofheartratedetection表2本文算法与TROIKA算法检测心率的平均绝对误差率Table2ComparisonofproposedalgorithmandTROIKAintermsofaverageabsoluteerrorpercentage(Error2)ofheartratedetection2.2算法设置\n实验需要2路PPG信号以及1路加速度信号，同时从心电信号提取基准心率作为结果参照。已知3路加速度信号具有较强的相关性，选用1路即可。因运动在水平方向表现明显，x轴的加速度信号含冇足够的运动信息，所以本文选用此路加速度信号。数据处理参照所对比的TROIKA方法，即单次处理的PPG信号、加速度信号和ECG信号的时间窗口长度为8s,每次窗口向前移动2s的长度，当前处理的信号与上次信号重叠时间为6s。每个窗口计算1次心率值。稀疏分解应用MMV模型，观测矩阵Y由两路PPG信号构成，则L=2。信号由125Hz降采样到25Hz,使稀疏矩阵维度降低，减少计算量，则单次稀疏观测矩阵行数N=25X8=200。作PPG信号、加速度信号1024点谱分析，即稀疏基的列数M=1024。此时频谱的分辨率为01X60^1.46次/分。解矩阵X的列向量表示两路PPG信号在调制稀疏基下的分解谱。因为T-MSBL算法不要求迭代到函数收敛，只需要观察在几次迭代后出现的非零大系数位置即可［10］，故将迭代次数设为4。T-MSBL算法相关参数分别初始如下：A=10,B为全0矩阵，Y-l,i=l,200,X为200X2全0矩阵。谱峰追踪屮，各项参数设置为Ai=5，A2=10,Th,=O.2,Th2=0.5。2.3结果评价基准心率是通过在每个时间窗内搜索ECG信号的R峰得出。选用平均绝对误差以及平均绝对误差率对本文方法的心率提取结果作评价。平均绝对误差定义为平均绝对误差率定义为其中，Q是心率计算的总次数，HR,.re>HRm分别为利用PPG信号和ECG信号求得的心率，平均绝对误差单位为次/分。表1和表2列出本文方法以及TROIKA分别作用于12组数据结果的平均绝对误差(Error1)和平均绝对误差率(Error2)。由表可知木文提出方法结果大多优于TROTKA。本文方法相较TROTKA心率估汁误差小，同时误差的均方值小表现出算法的稳定性较强。图7为使用本方法对第3组数据的心率检测结果。可看出结果能较好地跟踪基准心率。图8是Bland-Altman图，一致性范围(limitsofagreement,LOA)为[-7.31，6.77]beats/min(标准差8=3.59beats/min),有94%的点落在一致性界限内。图9是皮尔逊系数图，皮尔逊相关系数为0.989。\n图7第3组数据估计结果Fig.7EstimationresultsonDataset3下载原图图812组数据估计结果的Bland-Altman图Fig.8Bland-Altmanplotofestimationresultson12datasets图912组数据估计结果的皮尔逊相关Fig.9Pearsoncorrelationofestimationresultson12datascts卜载原图3结论利用PPG信号可使心率在可穿戴仪器上检测更加简单易行，但在运动过程中PPG信号极易受到干扰。本文提出的算法可较好地检测运动状态下的心率。对比最近提出的TROIKA算法，本文算法具有较高的准确率和较强的抗干扰能力。[l]PengRC，ZhouXL,LinWH，etal.Extractionofheartratevariabilityfromsmartphonephotoplethysmograms[J].Computational&MathematicalMethodsinMedicine，2015，2015:1-11.[2]PetersenCL，ChenTP，Ansermino丄M，etal.Designandevaluationofalow-costsmartphonepulseoximeter[J].Sensors，2013，13(12):16882-16893.[3]MendelsonY，OchsBD.Noninvasivepulseoximetryutilizingskinreflectancephotoplethysmography[J].IEEETransactionsonBiomedicalEngineering，1988，35(35):798-805.[4]TamuraT，MaedaY，SekineM，etal.Wearablephotoplethysmographicsensors—pastandpresent[J].Electronics，2014，3(2):282-302.[5]PengF，ZhangZ，GouX，etal.Motionartifactremovalfromphotoplethysmographicsignalsbycombiningtemporallyconstrainedindependentcomponentanalysisandadaptivefilter[J].BiomedicalEngineeringOnline，2014，13(1):694-694.\n[2]PflugradtM，OrglmeisterKImprovedsignalqualityindicationforphotoplethysmographicsignalsincorporatingmotionartifactdetection[C].EngineeringinMedicineandBiologySociety，2014AnnualInternationalConferenceoftheIEEE.IEEE,2014:1872-1875.[3]KrishnanR，NatarajanBB，WarrenS.Two-stageapproachfordetectionandreductionofmotionartifactsinphotoplethysmographicda.ta.[J].IEEEtransactionsonBiomedicalEngineering,2010，57(8):1867-1876.[4]RamMR,MadhavKV，KrishnaEH，etal.AnovelapproachformotionartifactreductioninPPGsignalsbasedonAS-LMSadaptivefilter[J].TEEETransactionsonInstrumentationandMeasurement,2012，61(5):1445-1457.[5]YousefiR，NouraniM，OstadabbasS，etal.Amotion-tolerantadaptivealgorithmforwearablephotoplethysmographicbiosensors[J].IEEEJournalofBiomedicalandHealthInformatics，2014，18(2):670-681.[6]ZhangZ，PiZ，TroikaLB.八generalframeworkforheartratemonitoringusingwrist-typephotoplethysmographicsignalsduringintensivephysicalexercise[J].IEEETransactionsonBiomedicalEngineering，2015,62(2):522-531.[7]DuarteMF,BaraniukRG.Spectralcompressivesensing[J].AppliedandComputationalHarmonicAnalysis,2013,35(1):111—129.[8]LccB,HanJ,Back11J,etal.ImprovedeliminationofmotionartifactsfromaphotoplethysmographicsignalusingaKalmansmootherwithsimultaneousaccelerometry[J].PhysiologicalMeasurement,2010,31(12):1585-1603.[9]JinY,RaoBD.Supportrecoverytolinearinverseproblemswithmilltiplemeasurementvectors[J].TEEETransactionsonInformationTheory,2013,59(5):3139-3157.[10]EldarYC，RauhutH.Averagecaseanalysisofmultichannelsparserecoveryusingconvexrelaxation[J].IEEETransactionsonInformationTheory，2009，56(1):505-519.[11]ZhangZ，RaoBD.Sparsesignalrecoverywithtemporallycorrelatedsourcevectorusingsparsebaycsianlearning[J].IEEEJournalofSelectedTopicsinSignalProcessing,2011，5(5):912-926.