如何检测心电图波形中的模式?

83
我正在尝试从心电图中读取图像,并检测其中的主要波(P波、QRS波群和T波)。我可以读取图像并获得一个向量(例如(4.2; 4.4; 4.9; 4.7; ...))。我需要一个算法,可以遍历这个向量并检测每个波的起始点和终止点。例如:

alt text

如果它们总是具有相同的大小,或者如果我事先知道ECG中有多少个波,那将很容易。给定波:

alt text

我提取向量:
[0; 0; 20; 20; 20; 19; 18; 17; 17; 17; 17; 17; 16; 16; 16; 16; 16; 16; 16; 17; 17; 18; 19; 20; 21; 22; 23; 23; 23; 25; 25; 23; 22; 20; 19; 17; 16; 16; 14; 13; 14; 13; 13; 12; 12; 12; 12; 12; 11; 11; 10; 12; 16; 22; 31; 38; 45; 51; 47; 41; 33; 26; 21; 17; 17; 16; 16; 15; 16; 17; 17; 18; 18; 17; 18; 18; 18; 18; 18; 18; 18; 17; 17; 18; 19; 18; 18; 19; 19; 19; 19; 20; 20; 19; 20; 22; 24; 24; 25; 26; 27; 28; 29; 30; 31; 31; 31; 32; 32; 32; 31; 29; 28; 26; 24; 22; 20; 20; 19; 18; 18; 17; 17; 16; 16; 15; 15; 16; 15; 15; 15; 15; 15; 15; 15; 15; 15; 14; 15; 16; 16; 16; 16; 16; 16; 16; 16; 16; 15; 16; 15; 15; 15; 16; 16; 16; 16; 16; 16; 16; 16; 15; 16; 16; 16; 16; 16; 15; 15; 15; 15; 15; 16; 16; 17; 18; 18; 19; 19; 19; 20; 21; 22; 22; 22; 22; 21; 20; 18; 17; 17; 15; 15; 14; 14; 13; 13; 14; 13; 13; 13; 12; 12; 12; 12; 13; 18; 23; 30; 38; 47; 51; 44; 39; 31; 24; 18; 16; 15; 15; 15; 15; 15; 15; 16; 16; 16; 17; 16; 16; 17; 17; 16; 17; 17; 17; 17; 18; 18; 18; 18; 19; 19; 20; 20; 20; 20; 21; 22; 22; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 32; 33; 33; 33; 32; 30; 28; 26; 24; 23; 23; 22; 20; 19; 19; 18; 17; 17; 18; 17; 18; 18; 17; 18; 17; 18; 18; 17; 17; 17; 17; 16; 17; 17; 17; 18; 18; 17; 17; 18; 18; 18; 19; 18; 18; 17; 18; 18; 17; 17; 17; 17; 17; 18; 17; 17; 18; 17; 17; 17; 17; 17; 17; 17; 18; 17; 17; 18; 18; 18; 20; 20; 21; 21; 22; 23; 24; 23; 23; 21; 21; 20; 18; 18; 17; 16; 14; 13; 13; 13; 13; 13; 13; 13; 13; 13; 12; 12; 12; 16; 19; 28; 36; 47; 51; 46; 40; 32; 24; 20; 18; 16; 16; 16; 16; 15; 16; 16; 16; 17; 17; 17; 18; 17; 17; 18; 18; 18; 18; 19; 18; 18; 19; 20; 20; 20; 20; 20; 21; 21; 22; 22; 23; 25; 26; 27; 29; 29; 30; 31; 32; 33; 33; 33; 34; 35; 35; 35; 0; 0; 0; 0;]

我希望能够检测到以下内容:

  • [19 - 37] 中的 P 波。
  • [51 - 64] 中的 QRS 波群。
  • 等等。

我认识一个在这个领域工作的人。你可以在这里找到他的出版物列表。如果我没记错的话,他使用隐马尔可夫模型来可靠地检测已知形状的训练集中的波浪,但你可以在论文中找到更多细节。 - Stefano Borini
您已经有了许多好的答案。我很惊讶没有人建议使用PhysioToolkit的'WFDB软件包', 特别是ecgpuwave - Amro
我对于在时间序列数据中检测模式的类似问题的回答在这里 - https://dev59.com/mGgt5IYBdhLWcg3w3xPC#11903770 - 包括了Python代码。我的方法是“切换自回归隐马尔可夫模型”(可以谷歌该短语以获取相关出版物)。 - user1149913
12个回答

55
我会做的第一件事就是看看已经有哪些研究。确实,这个特定问题已经得到了深入的研究。以下是一些非常简单的方法概述:链接
我必须回答另一个答案。我研究信号处理和音乐信息检索。表面上看,这个问题似乎与起始点检测类似,但问题背景并不相同。这种生物信号处理,即P、QRS和T波段的检测,可以利用每个波形的特定时域特征的知识。在MIR中,起始点检测并没有这样的特征。(至少不能可靠地检测)。
对于QRS检测而言,一个有效的方法是动态时间规整。当时域特征保持不变时,DTW可以工作得非常好。这是一篇使用DTW解决此问题的短小IEEE论文:链接
这是一篇很好的IEEE杂志文章,比较了许多方法:链接。您会发现已经尝试了许多常见的信号处理模型。浏览一下这篇论文,然后尝试理解其中一个基本的模型。
编辑:在浏览这些文章之后,基于小波的方法似乎最直观。DTW也可以很好地工作,而且已经有DTW模块存在,但是对我来说,小波方法最好。另外,有人通过信号的导数进行了回答。我的第一个链接研究了1990年之前的方法,但我怀疑它们不如更现代化的方法稳健。
编辑:我会尽快给出一个简单的解决方案,但我认为小波适合这里的原因是因为它们能够将各种形状进行参数化,而不考虑时间或幅度缩放。换句话说,如果你有一个具有相同重复时间形状但在不同时间尺度和幅度上的信号,小波分析仍然可以识别这些形状相似(粗略地说)。还要注意,我有点把滤波器组归为这个类别。类似的事情。


17
一部分谜题是 "起始检测",已经编写了许多复杂的算法来解决这个问题。以下是有关起始音符的更多信息。
下一步是 海明距离。这些算法允许您进行模糊比较,输入为2个数组,输出为整数"距离"或2个数据集之间的差异。数字越小,它们就越相似。这非常接近您所需要的,但不是完全准确的。我对海明距离算法进行了一些修改,以计算新的距离,它可能有一个名字,但我不知道是什么。基本上它将数组中每个元素之间的绝对距离相加并返回总和。以下是python代码。
import math

def absolute_distance(a1, a2, length):
       total_distance=0
       for x in range(0,length):
               total_distance+=math.fabs(a1[x]-a2[x])
       return total_distance

print(absolute_distance([1,3,9,10],[1,3,8,11],4))

这个脚本输出2,即这两个数组之间的距离。

现在来组合这些部分。您可以使用起点检测来查找数据集中所有波的开始位置。然后,您可以循环遍历这些位置,将每个波与样本P波进行比较。如果您命中QRS复合体,则距离将是最大的。如果您命中另一个P波,则数字不会为零,但它会小得多。任何P波和任何T波之间的距离都将非常小,但是如果您做出以下假设,则不会有问题:

任何p波和任何其他p波之间的距离都将小于任何p波和任何t波之间的距离。

系列看起来像这样:pQtpQtpQt... p波和t波紧挨着,但由于这个序列是可预测的,因此更容易阅读。

顺便说一句,这个问题可能有一个基于微积分的解决方案。但是在我看来,曲线拟合和积分使这个问题更加混乱。我编写的距离函数将找到面积差异,这非常类似于减去两个曲线的积分。

可能可以牺牲起始计算,而选择每次迭代1个点,从而执行O(n)距离计算,其中n是图中点的数量。如果你有一个所有这些距离计算的列表,并且知道有50个pQt序列,那么你就会知道50个最短距离,它们不重叠,并且所有位置都是p波。太棒了!对于简单性来说怎么样?然而,折衷方案是由于增加了距离计算的数量而导致效率下降。

更简单的技术,比如检测时域振幅增加,通常会导致过高的误报或漏报数量。这正是我担心的。在我的(不太理想的)解决方案中,我也提出了同样的建议。 - Vivin Paliath
3
你的算法很有意思,可能会有一定的成功。这是一个非常复杂的问题,没有完美的解决方案。 - rook

8
你可以使用互相关。对每个模式取一个样本,并将它们与信号进行相关分析。当相关性高时,您会得到峰值。我期望这种技术可以很好地提取qrs和t波。之后,您可以通过在相关信号中寻找qrs前的峰值来提取p波。
互相关是一种相当易于实现的算法。基本上:
x is array with your signal of length Lx
y is an array containing a sample of the signal you want to recognize of length Ly
r is the resulting correlation

for (i=0; i<Lx - Ly; i++){
  r[i] = 0;
  for (j=0; j<Ly ; j++){
    r[i] += x[i+j]*y[j];
  }
}

寻找r中的峰值(例如超过阈值的值)。

这是一个不错的初步尝试,因为波形总是遵循某种模式。但对于这个问题,时间缩放和幅度缩放都可能有所不同,因此最终,这种方法在不同对象中将无法稳健运作。 - Steve Tjoa
是的,这只是一个初步的尝试。不够健壮,但足够简单的编码来尝试一下。模式匹配通常是最简单的技术,仍然可以提供一些结果。当然,小波变换肯定更好。 - nacmartin

7

我会先简化数据。

不要分析绝对数据,而是分析一个数据点到下一个数据点的变化量。

这里有一个快速的一行代码,可以将以;分隔的数据作为输入,并输出该数据的差值。

perl -0x3b -ple'( $last, $_ ) = ( $_, $_-$last )' < test.in > test.out

运行您提供的数据,这是输出结果:

0;0;20;0;0;-1;-1;-1;0;0;0;0;-1;0;0;0;0;0;0;1;0;1;1;1;1;1;1;0;0;2;0;-2;-1;-2;-1;-2;-1;0;-2;-1;1;-1;0;-1;0;0;0; 0;-1;0;-1;2;4;6;9;7;7;6;-4;-6;-8;-7;-5;-4;0;-1;0;-1;1;1;0;1;0;-1;1;0;0;0;0;0;0;-1;0;1;1;-1;0;1;0;0;0;1;0;-1;1; 2;2;0;1;1;1;1;1;1;1;0;0;1;0;0;-1;-2;-1;-2;-2;-2;-2;0;-1;-1;0;-1;0;-1;0;-1;0;1;-1;0;0;0;0;0;0;0;0;-1;1;1;0;0;0; 0;0;0;0;0;-1;1;-1;0;0;1;0;0;0;0;0;0;0;-1;1;0;0;0;0;-1;0;0;0;0;1;0;1;1;0;1;0;0;1;1;1;0;0;0;-1;-1;-2;-1;0;-2;0; -1;0;-1;0;1;-1;0;0;-1;0;0;0;1;5;5;7;8;9;4;-7;-5;-8;-7;-6;-2;-1;0;0;0;0;0;1;0;0;1;-1;0;1;0;-1;1;0;0;0;1;0;0;0; 1;0;1;0;0;0;1;1;0;2;1;1;1;1;1;1;1;1;1;-1;1;0;0;-1;-2;-2;-2;-2;-1;0;-1;-2;-1;0;-1;-1;0;1;-1;1;0;-1;1;-1;1;0;-1; 0;0;0;-1;1;0;0;1;0;-1;0;1;0;0;1;-1;0;-1;1;0;-1;0;0;0;0;1;-1;0;1;-1;0;0;0;0;0;0;1;-1;0;1;0;0;2;0;1;0;1;1;1;-1; 0;-2;0;-1;-2;0;-1;-1;-2;-1;0;0;0;0;0;0;0;0;-1;0;0;4;3;9;8;11;4;-5;-6;-8;-8;-4;-2;-2;0;0;0;-1;1;0;0;1;0;0;1;-1; 0;1;0;0;0;1;-1;0;1;1;0;0;0;0;1;0;1;0;1;2;1;1;2;0;1;1;1;1;0;0;1;1;0;0;-35;0;0;0;

上述文本中插入了原本不在输出中的换行。


完成这个任务后,查找qrs复合物就很容易了。

perl -F';' -ane'@F = map { abs($_) > 2 and $_ } @F; print join ";", @F'< test.out

;;20;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;4;6;9;7;7;6;-4;-6;-8;-7;-5;-4;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;5;5;7;8;9;4;-7;-5;-8;-7;-6
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;4;3;9;8;11;4;-5;-6;-8;-8;-4;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;-35;;;

这里的数据点20-35是由于原始数据以0开头和结尾而产生的。

要找到其他数据点,您需要依靠模式匹配。


如果您查看第一个p波,您可以清楚地看到一种模式。

0;0;0;0;0;0;1;0;1;1;1;1;1;1;0;0;2;0;-2;-1;-2;-1;-2;-1;0;-2;-1;1;-1;0;-1;0;0;0;0;
#           \________ up _______/   \________ down _________/

不过,第二个P波上的模式并不容易看出来。这是因为第二个P波扩散得更远。

0;0;0;1;0;1;1;0;1;0;0;1;1;1;0;0;0;-1;-1;-2;-1;0;-2;0;-1;0;-1;0;1;-1;0;0;-1;0;0;0;
#     \________ up _______/       \________________ down ________________/

第三个p波比其他两个波略微不规则。
0;0;0;0;0;1;-1;0;1;0;0;2;0;1;0;1;1;1;-1;0;-2;0;-1;-2;0;-1;-1;-2;-1;0;0;0;0;0;
#                \_______ up ______/  \__________ down __________/

您需要以类似的方式找到t波,与p波相比,它们发生的时间不同。


这应该足够的信息让您开始了解了。

这两个一行代码仅为示例,不建议日常使用。


1
我只看到手动解决问题的方法,用户可以通过绘制数据并自行选择起始点来完成。 - Hannes Ovrén

4
那两个陡峭的山峰和山谷也是QRS波群吗?
我想你需要做的是计算每个点上的图形斜率。然后,您还需要查看斜率变化的速度(二阶导数???)。如果出现突然变化,则知道已经到达某种陡峭的山峰。当然,您希望限制检测的变化,因此可能要执行类似于“如果时间间隔T内斜率变化了X”,以便不捕获图形中的微小颠簸。
我已经很久没有做过任何数学了...而且这似乎是一个数学问题 ;)哦,我也没有进行任何信号分析:)。
只是添加另一点。我认为您还可以尝试信号平均。例如,对最后3或4个数据点进行平均。我认为您也可以通过这种方式检测到突然变化。

有趣的算法加1分。但我认为这个问题有点更复杂。 - rook
是的,另外两个峰和谷是QRS复合波。实际上,该图像有3个P波、3个QRS复合波和3个T波。这是一种有趣的方法,但如果我没有一个函数,我不知道如何计算第二导数。我想你是说给变化的值打分,并选择那些得分高的变化,比如某些事情的开始和结束,是吗?我会尝试一下,等我有了一些结果后再发布更新。感谢您的回答。 - Alaor
是的,差不多了。你有点在评分,但是你是通过计算斜率或者观察振幅随时间变化来做到这一点的。 - Vivin Paliath

4

我并不是这个具体问题的专家,但是根据我的一般知识,让我们假设您知道QRS波群(或其他特征之一,但我将在此示例中使用QRS波群)大致发生在长度为L的某个固定时间段内。我想知道您是否可以按照以下方式将其视为分类问题:

  1. 将信号分成长度为L的重叠窗口。每个窗口要么完全包含QRS波群,要么不包含。
  2. 对每个窗口进行傅里叶变换。您的特征是每个频率上的信号强度。
  3. 在一些手动注释的数据上训练决策树、支持向量机等模型。

3

一个很可能会产生良好结果的方法是曲线拟合:

  • Divide the continuous wave into intervals (probably it's best to have the interval borders about half way between the sharp peaks of the qrs complexes). Only consider a single interval at a time.
  • Define a model function that can be used to approximate all possible variations of electrocardiographic curves. This is not as difficult as it seems first. The model function can be constructed as a sum of three functions with parameters for the origin (t_), amplitude (a_) and width (w_) of each wave.

       f_model(t) = a_p   *  f_p  ((t-t_p  )/w_p) + 
                    a_qrs *  f_qrs((t-t_qrs)/w_qrs) +
                    a_t   *  f_t  ((t-t_t  )/w_t)
    

    The functions f_p(t), f_qrs(t), f_t(t) are some simple function that can be uses to model each of the three waves.

  • Use a fitting algorithm (e.g. the Levenberg-Marquardt-Algorithm http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm) to determine the fitting parameters a_p, t_p, w_p, a_qrs, t_qrs, w_qrs, a_t, t_t, w_t for the dataset of each intervall.

    The parameters t_p, t_qrs and t_p are the ones you are interested in.


3
这是一个很好的问题!我有一些想法: 动态时间规整 可能是一个有趣的工具。您可以为三个类别建立“模板”,然后使用DTW来查看信号“块”(将信号分成,例如,0-0.5、0.1-0.6、0.2-0.7秒等)与您的模板之间的相关性。我曾经用过类似的东西来进行加速度计数据的步态分析,效果还不错。
另一种选择是结合信号处理和机器学习算法。将信号再次分成“块”。再次制作“模板”(每个类别需要十几个)。对每个块/模板进行FFT,然后使用朴素贝叶斯分类器(或其他ML分类器,但NB应该够用)对三个类别进行分类。我也尝试过在步态数据上使用此方法,并能够获得相对复杂信号的98%以上的精确度和召回率。让我知道这个方法如何运作,这是一个非常令人兴奋的问题。

1
"

小波变换”可能是一个相关的关键词。我曾经参加过一个演讲,演讲者使用这种技术来检测嘈杂的心电图中不同的心跳阶段。

就我有限的理解而言,它有点像傅里叶变换,但是使用的是(缩放后的)脉冲副本,例如您的心跳形状。

"

1

首先,标准心电图波的各种组成部分可能会在任何给定的绘图中缺失。这样的绘图通常是不正常的,通常表明存在某种问题,但不能保证它们一定存在。

其次,识别它们与艺术和科学同样重要,特别是在出现问题的情况下。

我的方法可能是尝试训练神经网络来识别这些组件。您将向其提供前30秒的数据,归一化使最低点为0,最高点为1.0,它将有11个输出。不是异常评级的输出将是过去10秒钟的加权值。 0.0表示距离现在10秒钟之前,而1.0表示现在。输出将是:

  1. 最近一次P波开始的位置
  2. 最近一次P波结束的位置
  3. 最近一次P波的异常评级,其中一个极端是“缺失”。
  4. 最近一次QRS复合波开始的位置
  5. 最近一次QRS复合波中Q部分转变为R部分的位置。
  6. 最近一次QRS复合波中R部分转变为S部分的位置。
  7. 最近一次QRS复合波结束的位置。
  8. 最近一次QRS复合波的异常评级,其中一个极端是“缺失”。
  9. 最近一次T波开始的位置。
  10. 最近一次T波结束的位置。
  11. 最近一次T波的异常评级,其中一个极端是“缺失”。

我可能会与其他人建议的其他类型的分析进行双重检查,或者将这些其他类型的分析与神经网络的输出一起使用,以便为您提供答案。

当然,这个神经网络的详细描述不应该被视为规定性的。我相信我并没有选择最优的输出,例如,我只是随意地想出了一些可能的想法。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接