本发明属于医学信号处理领域，特别涉及一种从心电信号中准确提取qrs内异常电位的方法。
背景技术：
据统计，我国每年心脏猝死(scd)的总人数高达50多万，平均每分钟就有3人因心脏原因在发病1小时内死亡，而抢救成功率却不到1％。心脏猝死出现明显的年轻化趋势，部分患者发生猝死事件往往无明显先兆症状，身体貌似健康。所以，对心脏猝死进行早期预警显得尤为重要，但目前缺乏有效的检测手段。
大量研究表明，心脏内局部区域出现的电传导延迟是心脏猝死的一个重要诱因，这种传导延迟会引起信号折返，进而会导致严重的室性心率失常。区域传导延迟可以在体表心电图(ecg)的qrs复合波中出现切迹或不易察觉的小幅波动。这些特征往往在常规心电图中无法明确反映，需要使用采样频率更高的高频心电图。目前，临床上无创检测这种心室去极化过程出现的异常特征主要包括：(1)基于信号平均心电图(saecg)的心室晚电位(vlp)，(2)qrs内异常电位(aiqps)，(3)碎裂qrs波(fqrs)。
心室晚电位检查是指体表信息叠加心电图于qrs复合波终末并延伸到st段内的高频低振幅碎裂电位，它反映了缺血区心肌迟发的电活动，是由于缺血区心肌内缓慢而不规则的拆返活动引起。在急性心肌梗死后猝死的预测中，晚电位检查具有重要价值。
临床上使用较为广泛的是时域vlp检测方法，虽然这种分析具有高的阴性预测值，但阳性预测效果不佳。vlp检测方法阳性预测率不高的主要原因是仅检测qrs复合波终末并延伸到st段内的高频低振幅碎裂电位。大量的动物心肌梗死模型和人体标测基础研究结果明确表明，这些高频低振幅碎裂电位不但存在于qrs复合波终末区域，而且也可能隐藏在qrs复合波内。在某些心肌梗死部位，产生的高频低振幅碎裂电位可能仅存在于qrs复合波内，并不反映在qrs复合波终末处，如果能准确提取qrs复合波内的这些高频低振幅碎裂电位,将可显著提高心脏猝死早期预警的可靠性。
qrs复合波内的这些高频低振幅碎裂电位又称为qrs内异常电位(aiqps),提取aiqps是一个极具挑战性的任务，因为aiqps是嵌入在qrs波内的、非常微弱、且快速变化不可预测信号。目前，提出了多种方法来应对这一挑战。离散余弦变换(dct)域的自回归移动平均(arma)模型已被用于估计aiqps，其基本思想是用一个低阶arma模型来模拟正常qrs波，从而可提取出不可预测的aiqps。也有学者提出利用小波变换分析qrs波内的高频成分，并以此为基础诊断恶性室性心律失常。很多提取方法使用的是线性模型或线性变换技术，由于人体心脏工作机制十分复杂和精细，将qrs波群建模为非线性信号可能更接近实际情况。有学者提出一种用具有相同平滑度的径向基函数非线性神经网络逼近qrs波的方法，取得了较好的效果，但主要缺点是需要调整的参数过多，神经元最佳数量因人而异，参数设置不当会高估或低估rbf神经网络逼近误差，影响提取精度。
大量研究表明，对室性心律失常高风险患者，aiqps参数可显著提高诊断准确率，但目前aiqps提取技术的提取精度仍不能满足要求，而且提取技术的鲁棒性也需要进一步提高。
在心电测量过程中采样得到的原始心电信号用x(i)表示，其中包含工频干扰p(i)、基线漂移b(i)、测量噪声n(i)和含有待提取qrs内异常电位(aiqps)的aiqp(i)；不包含上述成分的理想心电信号用xp(i)表示；
x(i)＝xp(i)+aiqp(i)+p(i)+b(i)+n(i)，(1)
其中，aiqp(i)在aiqps有效区间内为待提取qrs内异常电位，其他部分取值为0。
x(i)去除工频干扰和基线漂移后的信号用x2(i)表示。
利用公式(2)可得到包含待提取qrs内异常电位和噪声干扰的信号y(i)：
y(i)＝x2(i)-xp(i)＝aiqp(i)+n(i)，(2)
在aiqps有效区间内，当n(i)的标准差与y(i)的标准差相比较小时，可以认为y(i)≈aiqp(i)，即为待提取qrs内异常电位。问题的关键是如何尽可能准确得到理想心电信号xp(i)、如何确定aiqps有效区间以及如何定量评价提取结果是否可信。
技术实现要素：
为了克服上述现有技术的缺陷，本发明的目的在于提供一种准确提取qrs内异常电位的方法，采用非线性变换预估技术，结合样条方法，获得不包含qrs内异常电位和其他干扰成分的理想心电信号xp(i)，并最终提取qrs内异常电位。
为了达到上述目的，本发明的技术方案为：
一种准确提取qrs内异常电位的方法，包括以下步骤：
步骤一、对原始心电信号x1(i)；进行预处理，得到预处理心电信号x2(i)；当原始心电信号是测量得到的单次心搏心电信号时，对其利用低通滤波器和工频陷波器进行处理，消除基线漂移和工频干扰对后续过程的影响；当原始心电信号是测量得到的包含多个心搏的心电信号，对其利用信号平均技术进行处理，消除基线漂移、工频干扰和测量噪声对后续过程的影响。
步骤二、对预处理心电信号x2(i)进行特征点检测，确定特征点位置和qrs范围，利用非线性变换得到预估理想心电信号；首先，对预处理心电信号进行特征点检测，确定特征点位置和qrs范围；其次，分别用两个不同滤波频率的低通滤波器，对经步骤一处理后得到的预处理心电信号进行滤波；然后，对得到的这两个滤波结果相减，得到差值信号，搜索该差值信号每个特征点位置前、后的第一个过零点位置；再后，每个特征点位置前、后的第一个过零点位置包含的时间范围，用上述两个不同频率中较高滤波频率的低通滤波器滤波结果代替，其他部分用上述两个不同频率中较低频率的低通滤波器滤波结果代替，得到合成信号；最后，对得到的合成信号进行低通滤波，得到预估理想心电信号；
步骤三、根据预处理心电信号、特征点位置和预估理想心电信号，利用样条方法得到精确估计理想心电信号。首先对经步骤一处理后得到的预处理心电信号和经步骤二得到的预估理想心电信号相减，得到误差信号，搜索该误差信号过零点位置；然后，在搜索得到的误差信号过零点位置和步骤二得到的特征点位置处，取样条权重为1，其他为0；最后，根据步骤一得到的预估理想心电信号和得到的样条权重，利用三次平滑样条得到精确估计理想心电信号；
步骤四、对经步骤一处理后得到的预处理心电信号和经步骤三得到的精确估计理想心电信号相减，该相减结果经一个带通滤波器滤波后得到滤波结果；根据该滤波结果和步骤二得到的qrs范围，通过移动标准差分析技术，获取qrs内异常电位；
步骤五、对获取的qrs内异常电位进行可信度评价。利用标准差分析方法，对经步骤四获取的qrs内异常电位可信性进行评价，确认步骤四获取的qrs内异常电位是否可信，并输出评价结果。
所述的步骤二具体为：
一、利用x2(i)进行ecg特征点检测，得到qrs范围，起始位置qrsb、终止位置qrse和ecg特征点位置p(j)，特征点数量为m，j＝1，2，…，m，所述ecg特征点至少应包含qrs起始点、qrs终止点，以及q、r、s波形峰值点；
二、使用较高频率的低通滤波器对x2(i)进行滤波，得到xhi，fh为低通滤波器滤波频率，100hz≤fh≤200hz；
三、使用较低频率的低通滤波器对x2(i)进行滤波，得到xl(i)，fll为低通滤波器滤波频率，40hz≤fll≤80hz；
四、利用公式(3)计算差值信号xd(i)：
xd(i)＝xh(i)-xl(i)，(3)
其中，xh(i)是较高频率低通滤波器对x2(i)的滤波结果，xl(i)是较低频率低通滤波器对x2(i)的滤波结果。
五、根据信号xd(i)，对每个ecg特征点时间位置p(j)，j＝1，2，…，m，反向、正方向分别搜索差值信号xd(i)，得到前、后第一个过零点，分别得到对应的时间位置pb(j)和pf(j)；
六、根据pb(j)和pf(j)构成该点集合set(j)：
set(j)＝{pb(j)，pb(j)+1，…，pf(j)-1，pf(j)}，(4)
以此为基础，构成集合seth：
seth＝{set(1)，set(2)，…，set(m)}，(5)
根据集合seth，通过公式(6)合成信号xs(i)，
七、对xs(i)进行低通滤波后得到预估理想心电信号x3(i)，f3为低通滤波器滤波频率，100hz≤f3≤200hz。
所述的步骤三具体为：
一、利用公式(7)计算误差信号xe(i)：
xe(i)＝x2(i)-x3(i)，(7)
x3(i)是预估理想心电信号；
二、利用公式(8)计算样条权重w(i)：
三、根据预估理想心电信号x3(i)和样条权重w(i)，利用三次平滑样条得到精确估计的理想心电信号的x4(i)。
所述的步骤四具体为：
一、利用公式(9)计算差值信号e(i)：
e(i)＝x2(i)-x4(i)(9)
x4(i)是精确估计的理想心电信号；
二、对e(i)进行带通滤波，得到包含待提取qrs内异常电位的信号y(i)，带通滤波带宽根据具体需要选择；
三、对信号y(i)计算移动窗口方差msd(i)，利用公式(10)计算msd(i)：
其中窗口长度为2k+1，k取值范围为2ms～5ms，msd(i)计算结果是k＝2ms。
四、是计算参考msd值refmsd，把qrs起始位置qrsb前100ms到qrsb的区间定义为参考区间，先分别计算参考区间内msd(i)的均值ref_msdmean和标准差值ref_msdstd，利用公式(11)计算refmsd：
refmsd＝ref_msdmean+α*ref_msdstd，(11)
其中，α一般选择大于2；
五、根据refmsd确定aiqps的开始位置aiqpb，具体方法是，从qrs起始位置qrsb开始正向搜索，如果msd(i)＞refmsd持续时间大于等于预先设定常数m，则停止搜索，此时的位置设为tb，aiqps的开始位置aiqpb按公式(12)计算：
aiqpb＝tb-m-k，(12)
其中，m一般取值为5ms；如果搜索到qrs终止位置qrse，则aiqpb＝0并停止搜索。
六、判断是否搜索到aiqpb，如果aiqpb等于0则退出并返回失败标志，否则继续。
七、根据refmsd确定aiqps的结束位置aiqpe，具体方法是，从qrs终止位置qrse后的大约50ms处开始反向搜索，如果msd(i)＞refmsd持续时间大于等于m，则停止搜索，此时的位置设为te，aiqps的结束位置aiqpe按公式(13)计算：
aiqpe＝te+m+k(13)；
如果搜索到aiqps的开始位置aiqpb，则aiqpe＝0并停止搜索，判断是否搜索到aiqpe，如果aiqpe等于0则退出并返回失败标志，否则继续。
八、提取qrs内异常电位aiqp(i)，按公式(14)计算：
其中，aiqpb是搜索得到的aiqps开始位置，aiqpe是搜索得到的aiqps结束位置。
所述的步骤五具体为：
一、计算参考区间标准差差refstd和qrs区域标准差qrsstd，refstd为参考区间y(i)的标准差，qrsstd为qrs起始位置qrsb到qrs终止位置qrse的区间内y(i)的标准差。
提取结果的可信度判断按公式(15)计算：
其中，β＞1，具体选择可根据实际情况确定。
如果可信度等于0则返回失败标志，否则返回成功标志并同时返回提取的qrs内异常电位aiqp(i)。
有益效果
本发明提出一种利用理想心电信号二次估计技术，精确提取qrs内异常电位的方法。在理想心电信号预估阶段，采用非线性变换技术，即可对非ecg特征点区域变化趋势进行有效跟踪，又可有效消除ecg特征点对提取理想心电信号可能造成的影响。根据预估理想心电信号，利用样条方法对理想心电信号进一步估计，可精确估计出理想心电信号。与现有方法相比，本发明方法需要选择的参数少，结果更可靠。本发明还对提取的qrs内异常电位进行可信度评价，从而保证利用本发明方法进行应用开发的结果可靠性。与传统的多次叠加平均方法相比，本发明方法的一个突出特点是可对单次搏动心电信号进行qrs内异常电位提取，可大大扩展aiqps分析技术的应用范围。
本发明方法可应用的场景和范围包括：1)在心电图机中集成本发明方法，可在进行常规心电图测量时，对患者进行心脏猝死危险性评价；2)在常规多参数监护仪中集成本发明方法，可对心肌梗死患者的病情变化进行实时动态跟踪监测；3)开发便于使用的便携式或可穿戴装置，实现在家庭环境下的心脏猝死危险性预警；4)在移动设备(如手机)中集成本发明方法，可为设备使用者提供一种便捷、高效的心脏猝死危险性预警手段。
附图说明
图1是本发明的流程图。
图2是本发明的原理描述使用的模拟信号。
图3是本发明实施例的预处理流程图。
图4是本发明实施例的基线漂移消除和工频干扰去除过程示意图。
图5是本发明实施例的预估理想心电信号流程图。
图6是本发明实施例的预估理想心电信号获取过程示意图。
图7是本发明实施例的精确估计理想心电信号流程图。
图8是本发明实施例的利用三次平滑样条得到精确估计理想心电信号过程示意图。
图9是本发明实施例的qrs内异常电位提取与结果可信度评价流程图。
图10是本发明实施例的qrs内异常电位提取与结果可信度评价的过程示意图。
图11是对单次搏动ecg和多次搏动ecg叠加平均，分别用本发明方法进行qrs内异常电位提取的结果图。
图12是用本发明所述方法对一例心肌梗死患者的单次搏动ecg进行qrs内异常电位提取的结果。
图13是用本发明所述方法对一例心脏健康者的单次搏动ecg进行qrs内异常电位提取的结果。
具体实施方式
下面结合附图对本发明的原理作详细描述。
参照图1，是本发明的流程图，一种准确提取qrs内异常电位的方法，包括以下步骤：
步骤一、对原始心电信号x1(i)；进行预处理，得到预处理心电信号x2(i)；当原始心电信号是测量得到的单次心搏心电信号时，对其利用低通滤波器和工频陷波器进行处理，消除基线漂移和工频干扰对后续过程的影响；当原始心电信号是测量得到的包含多个心搏的心电信号，对其利用信号平均技术进行处理，消除基线漂移、工频干扰和测量噪声对后续过程的影响。
步骤二、对预处理心电信号x2(i)进行特征点检测，确定特征点位置和qrs范围，利用非线性变换得到预估理想心电信号；首先，对预处理心电信号进行特征点检测，确定特征点位置和qrs范围；其次，分别用两个不同滤波频率的低通滤波器，对经步骤一处理后得到的预处理心电信号进行滤波；然后，对得到的这两个滤波结果相减，得到差值信号，搜索该差值信号每个特征点位置前、后的第一个过零点位置；再后，每个特征点位置前、后的第一个过零点位置包含的时间范围，用上述两个不同频率中较高滤波频率的低通滤波器滤波结果代替，其他部分用上述两个不同频率中较低频率的低通滤波器滤波结果代替，得到合成信号；最后，对得到的合成信号进行低通滤波，得到预估理想心电信号；
步骤三、根据预处理心电信号、特征点位置和预估理想心电信号，利用样条方法得到精确估计理想心电信号。首先对经步骤一处理后得到的预处理心电信号和经步骤二得到的预估理想心电信号相减，得到误差信号，搜索该误差信号过零点位置；然后，在搜索得到的误差信号过零点位置和步骤二得到的特征点位置处，取样条权重为1，其他为0；最后，根据步骤一得到的预估理想心电信号和得到的样条权重，利用三次平滑样条得到精确估计理想心电信号；
步骤四、对经步骤一处理后得到的预处理心电信号和经步骤三得到的精确估计理想心电信号相减，该相减结果经一个带通滤波器滤波后得到滤波结果；根据该滤波结果和步骤二得到的qrs范围，通过移动标准差分析技术，获取qrs内异常电位；
步骤五、对获取的qrs内异常电位进行可信度评价。利用标准差分析方法，对经步骤四获取的qrs内异常电位可信性进行评价，确认步骤四获取的qrs内异常电位是否可信，并输出评价结果。
图2是用于本发明原理描述的模拟心电信号x(i)，包含理想心电信号用xp(i)和模拟包含待提取qrs内异常电位aiqps(i)，以及工频干扰p(i)、基线漂移b(i)和测量噪声n(i)，采样率1000hz，数据长度n＝800。
所述的步骤一具体为：参加图3和图4，图3中101是去除基线漂移，利用基线漂移消除方法去除基线漂移，得到去除后的信号x1(i)。
目前，基线漂移消除方法有多种方法，由于基线漂移对最终提取结果没有明显影响，图4中的x1(i)是用三阶、1hz的butterworth数字滤波器设计低通滤波器参数，然后对x(i)用双向零相位滤波得到。
图3中102是去除工频干扰，根据得到的信号x1(i)，利用数字工频陷波器得到去除工频干扰和基线漂移后的信号x2(i)。图4中的x2(i)是利用数字工频陷波器得到的结果。
所述的步骤二具体为：参加图5和图6，其中图6是qrs区间附近的局部图。图5中201是特征点检测算法，利用x2(i)进行ecg特征点检测，得到qrs范围(起始位置qrsb、终止位置qrse)和ecg特征点位置p(j)，特征点数量为m，j＝1，2，…，m。所述ecg特征点至少应包含qrs起始点、qrs终止点，以及q、r、s波形峰值点。
图5中202是一个低通滤波器，使用较高频率的低通滤波器对x2(i)进行滤波，得到xh(i)，fh为202低通滤波器滤波频率，100hz≤fh≤200hz。图6中的xh(i)是用三阶、fh＝150hz的butterworth数字滤波器设计低通滤波器参数，然后对x2(i)用双向零相位滤波得到。
图5中203是一个低通滤波器，使用较低频率的低通滤波器对x2(i)进行滤波，得到xl(i)，fl为203低通滤波器滤波频率，40hz≤fl≤80hz。图6中的xl(i)是用三阶、fl＝60hz的butterworth数字滤波器设计低通滤波器参数，然后对x2(i)用双向零相位滤波得到。
图5中204是信号相减运算，利用公式(3)计算差值信号xd(i)：
xd(i)＝xh(i)-xl(i)，(3)
图6中的xd(i)是信号相减后的结果。
图5中205是过零点检测与集合seth计算，根据信号xd(i)，对每个ecg特征点时间位置p(j)，j＝1，2，…，m，反向、正方向分别搜索差值信号xd(i)，得到前、后第一个过零点，分别得到对应的时间位置pb(j)和pf(j)。
图6中的xd(i)中，虚线为零值线，符号“o”代表时间位置p(j)处对应xd(i)上的点，符号“δ”和“
对每个p(j)，根据pb(j)和pf(j)构成该点集合set(j)：
set(j)＝{pb(j)，pb(j)+1，…，pf(j)-1，pf(j)}，(4)
以此为基础，构成集合seth：
seth＝{set(1)，set(2)，…set(m)}，(5)
图5中206是信号合成，根据集合seth，通过公式(6)合成信号xs(i)，
图6中的xs(i)是得到的合成信号。
图5中207是一个低通滤波器，对xs(i)进行滤波后得到预估理想心电信号x3(i)，f3为207低通滤波器滤波频率，100hz≤f3≤200hz。图6中的x3(i)是是用三阶、f3＝150hz的butterworth数字滤波器设计低通滤波器参数，然后对xs(i)用双向零相位滤波得到的预估理想心电信号
所述的步骤三具体为：参加图7和图8，其中图8是qrs区间附近的局部图。图7中301是信号相减运算，利用公式(7)计算误差信号xe(i)：
xe(i)＝x2(i)-x3(i)，(7)
图8中的xe(i)是信号相减后的结果。
图7中302是计算样条权重，利用公式(8)计算样条权重w(i)：
图8中的w(i)是根据公式(7)得到的样条权重。
图7中303是三次平滑样条运算，根据预估理想心电信号x3(i)和样条权重w(i)，利用三次平滑样条得到精确估计的理想心电信号的x4(i)，
图8中的x4(i)是得到的精确估计的理想心电信号。
所述的步骤四具体为：参加图7、图9和图8、10。
图7中304是信号相减运算，利用公式(9)计算差值信号e(i)：
e(i)＝x2(i)-x4(i)。(9)
图7中305是一个带通滤波器，通过对e(i)进行带通滤波，得到包含待提取qrs内异常电位的信号y(i)，f1和f2分别为带通滤波器的低频频率和高频低频，f1和f2可根据具体后续应用选择。图8中的y(i)是用五阶、f1＝70hz、f2＝300hz的butterworth数字滤波器设计带通滤波器参数，然后对e(i)用双向零相位滤波。
图9中401是对信号y(i)计算移动窗口方差msd(i)，利用公式(10)计算msd(i)：
其中窗口长度为2k+1，k一般取值范围为2ms～5ms。图10中msd(i)计算结果是k＝2ms。
图9中402是计算参考msd值refmsd。把qrs起始位置qrsb前大约100ms到qrsb的区间定义为参考区间，先分别计算参考区间内msd(i)的均值ref_msdmean和标准差值ref_msdstd，利用公式(11)计算refmsd：
refmsd＝ref_msdmean+α*ref_msdstd，(11)
其中，α一般选择大于2，具体选择可根据实际情况确定。
图10msd(i)中的水平虚线的幅度值代表refmsd值，α＝3。
图9中403是根据refmsd确定aiqps的开始位置aiqpb。具体方法是，从qrs起始位置qrsb开始正向搜索，如果msd(i)＞refmsd持续时间大于等于m，则停止搜索，此时的位置设为tb，aiqps的开始位置aiqpb按公式(12)计算：
aiqpb＝tb-m-k，(12)
其中，m一般取值为5ms；如果搜索到qrs终止位置qrse，则aiqpb＝0并停止搜索。
图9中404是判断是否搜索到aiqpb，如果aiqpb等于0则退出并返回失败标志，否则继续。
图9中405是根据refmsd确定aiqps的结束位置aiqpe。具体方法是，从qrs终止位置qrse后的大约50ms处开始反向搜索，如果msd(i)＞refmsd持续时间大于等于m，则停止搜索，此时的位置设为te，aiqps的结束位置aiqpe按公式(13)计算：
iqpe＝te+m+k；(13)
如果搜索到aiqps的开始位置aiqpb，则aiqpe＝0并停止搜索。
图9中406是判断是否搜索到aiqpe，如果aiqpe等于0则退出并返回失败标志，否则继续。
图9中407是提取qrs内异常电位aiqp(i)，按公式(14)计算：。
图10aiqp(i)为提取的qrs内异常电位aiqps，其中两条垂直虚线分别代表aiqpb和aiqpe。
所述的步骤五具体为：参加图9和图10。图9中408是计算参考区间标准差差refstd和qrs区域标准差qrsstd，refstd为参考区间y(i)的标准差，qrsstd为qrs起始位置qrsb到qrs终止位置qrse的区间内y(i)的标准差。
图9中409是提取结果的可信度判断，可信度按公式(15)计算：
其中，β＞1，具体选择可根据实际情况确定。
如果可信度等于0则返回失败标志，否则返回成功标志并同时返回提取的qrs内异常电位aiqp(i)。
图11是使用本发明描述中的模拟信号产生方法，分别产生1次、10次、和100次模拟心搏ecg，然后分别对多次模拟心跳ecg进行叠加平均，对得到的1次心搏ecg、以及10次和100次模拟心搏平均ecg，分别用本发明所述方法提取aiqp(i)。1次、以及10次和100次叠加平均提取的aiqp(i)，与模拟待提取aiqps(i)相比，相关系数和均方误差分别为：0.95、0.1；0.87、0.32；0.87、0.29。结果表明，本发明方法仅需要单次心搏ecg就可提取aiqp(i)，与多次叠加平均相比，在测量噪声足有小的情况下，单次心搏ecg提取的aiqp(i)准确度更高。从图11中也可看出，随着叠加平均次数增加，本发明步骤四中参考区间的干扰越来越小。所以，如果由于测量得到的ecg信号干扰较大，导致从单次心搏ecg中提取失败，可对多个心搏ecg进行叠加平均提取。
图12是用本发明所述方法对一例心肌梗死患者的单次心搏ecg进行aiqp(i)提取的结果，图13是用本发明所述方法对一例心脏健康者的单次心搏ecg进行aiqp(i)提取的结果。对比图12和图13可以发现，心肌梗死患者的aiqp(i)幅度明显大于心脏健康者的aiqp(i)幅度，且两者aiqp(i)的形态也有明显差异，这些特征可用于心脏猝死的早期预警。