CN1976628A - 信号分析方法 - Google Patents

Abstract

本发明公开了通过利用持续时间约束，对信号如心电图信号分割可靠性的改进。利用隐马尔可夫模型分析信号。持续时间约束确定了模型特定状态的最小允许持续时间。持续时间约束可以结合模型本身或结合维特比算法，用于在给定常规模型条件下，计算最概然状态序列。还公开了从模型推导出置信度，置信度可用于评估分割的质量和鲁棒性以及识别分割不可靠的任何信号，例如由于信号中噪声或异常的存在。

Description

信号分析方法

技术领域

本发明涉及一种用于信号分析的方法。通过利用隐马尔可夫建模技术结合持续时间约束，所述的分析可以包括对信号比如心电图(ECG)或其他生物医学信号进行自动分割，使其变成成份波形特征曲线。可选地，还可以结合对表示分割精确度和鲁棒性的关联置信度测量进行推导。

背景技术

ECG(也就是人们所知道的缩略词EKG)是测量心脏电功率的一种重要的非侵入信号。每个单独的心跳由许多不同的心脏病学阶段构成，这些心脏病学阶段依次引起一组ECG波形中不同的特征曲线。这些特征曲线表现了心脏特定区域肌肉细胞的去极化(放电)或复极化(充电)。图1显示了人的ECG波形和相关联的特征曲线。ECG波形的标准特征曲线是P波、QRS复合波和T波。另外小U波(在T波之后)偶尔也会出现。

心动周期从对应于心脏的心房去极化周期的P波(其起点和终点被称为P_on和P_off)开始。之后是QRS复合波，QRS复合波通常是ECG波形中最易识别的特征曲线，它对应于心室去极化周期。QRS复合波的起点和终点被称为Q点和J点。T波在QRS复合波之后对应于心室复极化周期。T波的终点被称为T_off，它代表心动周期的结束(假设没有U波的话)。

通过详细地检查ECG信号，有可能从典型的ECG波形中获得出许多含有信息的量度。这可以用于评估病人的病理状况以及探测心脏节律中存在的任何潜在异常。

一个特别重要的量度是“QT间隔”，其被定义为从QRS复合波的起点到T波终点的时间，即T_off-Q。这个定时间隔对应于心室中电功率(去极化和复极化)的总的持续时间。

QT间隔尤其重要，因为比正常间隔更长的间隔是长QT综合症(LQTS)的一个可靠指标。这是使患者易受过快反常心律(心律不齐)伤害的潜在致命条件，即所谓的“尖端扭转性室速”的。当这种心律发生时，心脏不能有效地跳动，而且通向大脑的血液循环急剧下降。结果是突然失去意识并且可能发生心脏死亡。

另一个重要的量度是PR间隔，其被定义为从P波起点到QRS复合波起点的时间，即Q-P_on。这对应于从心房去极化开始到心室去极化开始的时间，并能用于探测如房室阻滞等的状况。

除了对QT和PR间隔的测量与评估之外，把ECG波形精确而鲁棒的分割成其成分曲线也很重要，用于对诸如ST升高等状况的探测，并用于通过定位连续心博中QRS复合波顶点(即所谓的R点)对心率进行计算。

既然已经知道精确测定ECG波形中特征曲线边界的上述重要性，因此就要付出大量的努力来自动化该过程。在这个领域最初的工作集中于传统的信号处理和模式识别方法如阈值法和模板匹配法。最近，基于隐马尔可夫模型(HMMs)的统计法已经被应用于ECG的分割问题。

在隐马尔可夫建模中，存在无法被直接观测到的所关心的隐含状态序列。另外尽管如此，还存在与此状态序列概率相关的信号或观测值序列，而且从此我们可以推出给定状态序列的最概然值。

隐马尔可夫模型的核心是两个概率函数，其中一个涉及在特定时间步长且给定前一时间步长的该状态值的条件下，所关心的状态的概率，另一个涉及在给定特定状态值的条件下，信号或观测值的概率，。

隐马尔可夫建模方法特别适合于ECG分析，因为它利用了ECG波形特征曲线的独特统计特征。另外，模型还能利用波形特征曲线的连续性(即T波在QRS波之后，QRS复合波又在P波之后)，这种连续性帮助改进嘈杂的ECG信号分割。

可惜，用隐马尔可夫模型的一个重大不利因素就是他们易于产生非常不可靠的分割。这些分割以被分割的特征曲线仅仅持续少数时间采样点为特征。这个问题源于模型对非常短持续时间特征曲线的偏差，并且因此使标准的HMMs不适于某些应用，如要求高度鲁棒性的ECG分割。

图2和图4显示了两个采样ECG波形包括，由专门的ECG分析员确定的波形边界(参见实心垂直线)，和那些利用标准隐马尔可夫模型推出的波形边界(参见虚线垂直线)。在两个实例中模型找到了两个QRS复合波和两个T波，当原始信号中每种波只存在一个的时候。特别地，模型分割包括一非常短的QRS复合波和T波，它们直接出现在P波的起点后和真正的QRS复合波位置之前。这样的“双搏”分割严重影响了模型的性能和可靠性。

发明内容

考虑到这个背景情况，本发明的一个目的是关于改进信号如ECG信号的分割的可靠性。

本发明的另一个目的是从模型推导出置信度。此置信度可以用于评估分割的质量和鲁棒性，以及探测任何太嘈杂或太异常的ECG波形，因为如果ECG太嘈杂或太异常，通过自动系统不能可得到可靠的分割。

根据本发明的一个方面，提供了一种用于用隐马尔可夫模型分割生物医学信号的计算机执行方法，该模型包括多个状态，该方法包括：

为至少一个状态指定一最小持续时间约束d_min；

对于具有指定最小持续时间的模型中的每个状态，用一组亚态替换该状态，亚态的总数等于最小持续时间约束d_min的值；

将该组亚态连接在一起，形成左右方向的马尔可夫链，其中，第一d_min-1亚态中的每个各有为0的自跃迁概率，为1的跃迁到紧接它们右边的状态的跃迁概率，和为0的跃迁到模型中任何其他状态的跃迁概率；以及

将该模型应用于表示生物医学信号的数据上，以获得关于将信号分割成状态的信息。

根据本发明的又一方面，提供了一种计算机执行方法，用于通过使用对维特比算法的修正，将包括一列观测值的信号分割为有限状态离散时间马尔可夫过程的一列状态，该修正包括：

对于将结合入维特比算法的有限状态离散时间马尔可夫过程的至少一个状态，定义一个持续时间变量和持续时间约束，该持续时间约束指定所述至少一个状态的最小持续时间；

将修正后的维特比算法应用于信号，以计算构成观测值序列的最概然约束持续时间的状态序列；其中，

在为有限状态离散时间马尔可夫过程的每个状态计算最概然状态序列的每个时间步长，该马尔可夫过程构成了直到那个时间步长内的观测值序列，并以所述状态结束：

对于具有持续时间约束的每个状态，利用那个状态的持续时间变量跟踪先行状态连续序列的长度，该先行状态仅由那个状态组成并在前一时间步长以那个状态结束；

如果那个状态的持续时间变量大于或等于对那个状态指定的持续时间约束，那么在给定时间步长的状态序列计算中，考虑从那个状态到马尔可夫过程中的任何其他给定状态的跃迁；

如果那个状态的持续时间约束小于对那个状态指定的持续时间约束，那么在给定时间步长的状态序列计算中，不考虑从那个状态到马尔可夫过程中任何其他状态的跃迁；并且

在计算出直到给定时间步长的那组最概然状态序列之后，更新具有持续时间约束的每个状态的持续时间变量，用以跟踪先行状态连续序列的长度，该先行状态仅由那个特定状态组成，并在所考虑的时间步长以那个状态结束。

本发明的另一方面提供了一种用于分析已根据概率分割算法分割的信号的计算机执行方法，该方法包括：

计算每个已分割信号特征曲线的置信度；

画出该置信度对各自信号特征曲线长度的坐标图；

将密度建模技术应用于测定数据空间中与高置信度特征曲线相关联的适当区域；

测定一特定信号特征曲线的置信度是否落到该区域以外。

附图说明

本发明的具体实施例将仅通过实例，参照附图来描述，附图中：

图1显示了一充分标注的ECG波形；

图2显示了一ECG波形，包括由专门的ECG分析员确定的波形边界和利用标准隐马尔可夫模型方法推出的波形边界；

图3显示了与图2中使用的同样的ECG波形，包括由专门的ECG分析员确定的波形边界和利用新的持续时间约束方法推出的波形边界；

图4显示了一第二ECG波形，包括由专门的ECG分析员确定的波形边界和利用标准隐马尔可夫模型方法推出的波形边界；

图5显示了与图4中使用的同样的ECG波形，包括由专门的ECG分析员确定的波形边界和利用新的持续时间约束方法推出的波形边界；

图6是详细描述新的持续时间约束维特比算法运算的流程图；

图7a显示了用于把ECG信号分割为P波区、基线区、QRS复合波区、T波区、可选U波区和第二基线区的标准隐马尔可夫模型的结构。另外，指出了P波的跃迁概率a和观测值模型概率分布b；

图7b显示了用于对ECG信号以类似于图7a模型的方式进行分割的标准隐马尔可夫模型的结构，除了没有可选U波；

图7c显示了用于把ECG信号分割为PR间隔区、QRS复合波区、T波区和基线区的标准隐马尔可夫模型的结构；

图7d显示了用于把ECG信号分割为PR间隔区、QT间隔区和基线区的标准隐马尔可夫模型的结构；

图7e显示了用于把ECG信号分割为QT间隔区，和与信号的剩余部分对应的第二区(X)的标准隐马尔可夫模型的结构；

图8显示了把持续时间约束结合到隐马尔可夫模型的图解说明；

图9显示了用于ECG分割的应用嵌入持续时间约束的隐马尔可夫模型的部分结构；

图10显示了置信度对ECG波形长度的坐标图，用于许多规则ECG波形，并在99％置信下限条件下；

图11显示了包含两个受噪声(特别是“肌肉伪”噪声)干扰的ECG波形的ECG信号的一部分。还显示了两个嘈杂的波形(以A和B标注)的分割；以及

图12显示了图11中显示的两个嘈杂ECG波形的置信度。

具体实施方式

本发明的以下具体实施例是基于隐马尔可夫模型(HMM)的基础上的，该隐马尔可夫模型被训练以识别ECG信号特征曲线中的一个或多个。该模型由多个状态组成，每个状态代表了一个ECG信号的一特定区域。在图7a中显示了用于ECG分割的一类隐马尔可夫模型的结构或“拓扑”的图形描绘。该模型由六个单值的状态组成，其依次代表P波、P波的终点与QRS复合波起点间的基线部分(定义为“基线1”)、QRS复合波、T波、U波和T波(或U波，如果有一个U波的话)终点与下一心搏的P波起点间的基线部分。在图1中表示出了ECG波形的这些部分。在图7b-7e中图示了用于ECG分割的许多可选的隐马尔可夫模型结构类型。在每种情况中，模型由与被观测信号(由暗节点表示)随机相关的“隐”状态序列(由亮节点表示)组成。为了ECG分割，隐状态s_t表示在时间t起作用的特定波形特征曲线，被观测信号样本O_t表示ECG波形的关联信号样本。

模型取决于以下三个参数：i.)初始状态分布π，ii.)跃迁矩阵A，以及iii.)HMM中每个状态k的观测值概率模型b_k。训练HMM包括调整这些参数以使特定误差函数最小化。

给定一组适当的模型参数值，HMM采取ECG信号生成概率模型的形式。更确切地说，该模型定义了一用于生成ECG波形的随机过程。此过程首先选择初始状态，这通过从初始状态分布π中采样来实现。挑选出特定初始状态然后对应于ECG信号的第一波形特征曲线(即在时间采样点t＝1处起作用的特征曲线)。然后此波形特征曲线的关联ECG信号样本通过从那个特定状态的对应观测值概率模型b_k中采样产生。

该过程中的下一阶段是随即选择下一个模型状态，也就是在时间采样点t＝2处起作用的波形特征曲线。这是通过从跃迁矩阵采样来实现。该矩阵采取一简单表格的形式，该表格为模型中的每个状态定义了下一个状态发生的概率。一旦按照这样的方式选择了下一个状态，那么模型通过从适当的观测值概率模型中采样生成关联ECG信号样本。这个过程在余下的时间步长中连续叠代。

给定一特定ECG信号，将该信号分割成由模型结构定义的不同所关心特征曲线或区域(即定位每个模型状态的信号边界)包括颠倒上述“生成过程”。因而，目的是推出引起给定ECG信号的最概然隐状态序列。然后该状态序列定义ECG信号中的所关心区域的边界，和如此相应的分割。

这个分割过程可以通过使用维特比算法实现。更确切地说，维特比算法将一给定的ECG信号和一已训练的隐马尔可夫模型作为输入，并且返回最概然隐状态序列，如在L R Rabine，“隐马尔可夫模型指南及在语音识别中的精选应用”，IEEE学报，77：257-286，1989中所述。

为了模型能够成功完成分割ECG波形的任务，有必要使用具有鲁棒性的且能够提供信息的ECG编码方法。一个非常有效的代表是来自ECG信号非抽取小波变换(UWT)的一组系数，如N P Hughes、L Tarassenko和S J Roberts在“用于自动化的ECG时间间隔分析的马尔可夫模型”，神经信息处理系统的改进16，麻省理工学院出版社，2003年，中所述。

一旦已经选择了ECG信号的代表，训练HMM包括调整模型的参数以使ECG数据的概率(或可能性)即p(O|λ)最大，其中O表示ECG信号的数据集合，λ表示一组HMM参数。训练可以通过有监督或无监督的方式实现。

在有监督的学习中，训练的数据包括许多ECG波形和相应的状态标签。因而对于来自于一个给定ECG波形中的每个信号样本O_t，我们还有可利用的标签l_t，该标签指出该样本属于哪个特定状态(即波形特征曲线)。于是可以利用下面的极大似然估计计算出最初的状态分布和跃迁矩阵：

           π_i＝m_i/N

$a_{\mathrm{ij}}=n_{\mathrm{ij}}/\underset{k}{\mathrm{\Σ}}{n}_{\mathrm{ik}}$

其中m_i是第一信号样本(遍及ECG波形)属于状态i的总次数，n_ij是在所有标签序列中从状态i到状态j的总跃迁数，N是在数据集合中单个ECG信号的数目。

通过提取属于那个特定状态的所有信号样本，然后使观测值模型直接配合那些样本，从而可以学习对于每个状态k的观测值概率模型b_k。对于高斯观测值模型，这仅仅包括计算给定样本的平均值和方差。对于高斯混合模型，有必要利用EM算法(最大期望值算法)以从给定样本中学习观测值模型参数。

在无监督的学习中，训练数据只包括ECG波形。于是可以利用Baum-Welch算法(更普遍的EM算法的一个特例)来训练隐马尔可夫模型。

为了防止模型产生不可靠的分割，本发明的具体实施例在模型上施加了一组持续时间约束。这些约束限制了模型推出在生理学上似乎不合情理的分割的能力(在关联波形特征曲线的持续时间方面)。

特别地，对于模型中的每个状态，持续时间约束表现为单数d_min的形式，它表示对于那个特定状态(在“样本”中的)的最小允许持续时间。这些持续时间约束可以用许多不同的方法来估计。如果带标签的ECG数据可用，那么对于每个波形特征曲线的持续时间约束可以被简单地估计为数据集合中存在的那个特征曲线的最小持续时间的适当部分。例如，最小T波持续时间可以被设为数据集合中存在的最小T波持续时间的80％。

如果带标签的数据不可用，并且我们只有原始的ECG信号可用，那么持续时间约束必须从数据中学习。这可以通过利用训练好的模型去分割数据集合中的ECG信号，然后估计来自于所得到的分割持续时间约束来实现(采用与先前描述的对于带标签的数据类似的方式)。在这种情况下，首先必须要排除任何不可靠的“双搏”分割，否则会影响估计的持续时间约束的质量。

一旦已经计算出了持续时间约束，那么它们可以用于改进隐马尔可夫模型的分割性能。这可以通过将持续时间约束结合维特比算法，或将持续时间约束结合HMM结构来实现。现在将详细研究这两种方法。

维特比算法是使隐马尔可夫模型可以在实践中成功运用的关键的推理过程，如在L R Rabine，“隐马尔可夫模型指南及在语音识别中的精选应用”，IEEE学报，77：257-286，1989，中详细描述的。在进行ECG分析的背景下，维特比算法使我们可以在给定ECG波形中，为P波、QRS复合波和T波定位最可能的边界位置(即开始和偏移)。

为了将持续时间约束结合入维特比算法，有必要限制算法返回的状态序列，使得每个关联波形特征曲线的持续时间服从给定的持续时间约束。

更规范地说，给定一个HMMλ(与K个状态)和一ECG信号O＝O₁O₂…O_T，我们想要推出满足持续时间约束d_min(k)，k＝1…K的最概然状态序列S＝s₁s₂…s_T。每个持续时间约束具体要求任何状态k的“运行”持续时间应该至少为d_min(k)个采样时间。这样的目的是为确保我们不会推出任何我们预先已经知道无用(在单个状态运行的持续时间方面)的状态序列。

标准维特比算法是基于一个用于计算最概然状态序列可能性的高效过程，

计算出第一个t观测值并且结束在状态i：

${\mathrm{\δ}}_t\left(i\right)=\underset{s_1{s}_2...s_{t-1}}{\max }p(s_1{s}_2...s_t=i,O_1{O}_2...O_t|\mathrm{\λ})$

该算法的核心是下面用于计算δ_t(i)的递归关系式：

${\mathrm{\δ}}_t\left(i\right)=\underset{j}{\max }\left\{{\mathrm{\δ}}_{t-1}\left(j\right)a_{\mathrm{ji}}\right\}{b}_i\left(O_t\right)$

该式利用了该状态序列是一阶马尔可夫的事实(也就是说在δ_t(i)的计算中我们只需要考虑前一状态)。于是在任何时间t，对于穿过状态格子的最大可能路径，总是有可能从其当前状态跃迁到一不同的状态，假设关联的跃迁概率非零。因此为了将持续时间约束结合入维特比算法，这条路径必须被限制，使得从一个状态到另一状态的任何跃迁只有当当前所关心状态已经被占据，并且达到大于或等于最小允许持续时间时才被允许。

根据本发明的该第一具体实施例，约束持续时间的维特比算法进行如下。在递归式的每个时间步长t，我们持续跟踪在时间t-1结束的状态i的当前运行的持续时间d(i)。更确切地说，每个状态的持续时间变量被用于跟踪先行状态连续序列的长度，这些先行状态仅包括那个特定状态，并且在前一时间步长的那个状态处结束。如果该持续时间大于或等于最小允许持续时间d_min(i)，在计算δ_t(j)时我们考虑从这个状态到其他任何状态j的跃迁(和应用标准维特比递归式一样)。但是如果当前状态占据时间小于允许的最小值，则不允许从这个状态到不同状态的跃迁，并且该状态只能进行自跃迁(直到它的最小持续时间约束被满足)。一旦计算出了所有t值的δ_t(i)，可以用标准回溯方法来产生最终状态序列。

图6的流程图示意性地说明了约束持续时间的维特比算法的全部过程。算法第一步是输入HMM参数、被观测信号和对于每个状态的最小持续时间约束。这在框10中说明。

算法下一步，如框20中所描述的，涉及所关心变量的初始化。更确切地说，对于每个状态i，我们设δ₁(i)＝π_ib_i(O_l)，ψ₁(i)＝0和d(i)＝1。和标准维特比过程一样，变量ψ_t(i)用于存储在时间t-1的最可能状态，假设该模型在时间t时处于状态i。变量d是只有约束持续时间的维特比过程才有的，并且被用于持续跟踪每个结束在时间采样点t-1时的状态的当前“运行”持续时间。更确切地说，对于每个状态i，持续时间变量d(i)用于持续跟踪先行状态连续序列的长度，该先行状态仅包括状态i，并且在时间步长t-1时结束。

在框30和40中显示了约束持续时间的维特比算法的主要部分。这两个过程在从2到T(其中T是信号的长度)的所有时间采样点处被连续叠代。框30在给定的时间步长t处，为每个状态i计算δ_t(i)和ψ_t(i)的值。在标准维特比算法中，分别用于计算δ_t(i)和ψ_t(i)的运算符“max”和“argmax”被估算，用于所有可能的先前状态(从1到K)。但是在约束持续时间的维特比算法中，这些运算符只有当被用于那些满足其最小持续时间约束d_min的状态时会被估算。另外，这些运算符当被用于当前状态i时也被估算，使得最可能的前一状态总是被允许作为我们当前正在考虑的那个状态。

假设如此，则框30的第一步是计算满足其最小持续时间约束的状态j(包括默认状态i)的集合，即那些状态k满足d(k)≥d_min(k)。然后，δ_t(i)和ψ_t(i)的计算仅需为该状态j的集合进行估计。

约束持续时间的维特比算法的下一步是更新每个状态的持续时间变量。更确切地说，如果最可能的前一状态ψ_t(i)也是状态i，那么这个状态的持续时间变量增加1，即d(i)＝d(i)+1。这对应于“自跃迁”的情况，在这种情况中，模型在连续的时间步长内保持同样的状态，并且因此使得在时间采样点t结束的状态i的当前“运行”持续时间增加1。但是如果最可能的前一状态ψ_t(i)不是状态i，那么这个状态的持续时间变量被重设为1，即d(i)＝1。这对应于从一个不同的状态跃迁到状态i的情况，并且因此结束状态i自跃迁的当前“运行”。

一旦框30和40被计算直到包括时间采样点T，标准维特比回溯过程可以用于提取最概然状态序列，通过利用先前计算出的δ和ψ变量。这个回溯过程的第一步，如框50所示，是找出使δ_T值最大的状态，即 $s_T=\underset{i}{\arg \max }\left\{{\mathrm{\δ}}_T\left(i\right)\right\}.$ 一旦计算出了这个值，在前一时间步长(即T-1)的最概然状态可以简单地从ψ_T(s_T)值中找出。利用这个递归过程，可以计算出所有剩余时间采样点的最概然状态值。

根据本发明的第二具体实施例，改进HMM分割性能的可选方法是将持续时间约束直接与模型结构相结合。一旦持续时间约束通过这种方式被结合，那么标准维特比算法可以用于推出给定ECG波形的最佳约束持续时间分割。

将持续时间约束结合入HMM结构包括将一组d_min新状态替换模型中的每个状态。这些状态以简单的左右方向的方式被连接，并且共享一个共同的观测值密度。

更确切地说，对于具有一d_min(k)的最小持续时间的每个状态k(其中d_min(k)＞1)，我们通过在初始状态k前直接加上d_min(k)-1个附加状态来增大模型。每个附加状态自跃迁概率为0，跃迁到紧接其右边的状态的跃迁概率为1，跃迁到模型中任何其他状态的跃迁概率为0。从而将这些状态放在一起形成一简单的左右方向的马尔可夫链，其中该链中的每个状态只被占据最多一个时间采样点(在该链内的任何运行期间)。

这个链最重要的特征是，每个状态的观测概率模型的参数与初始状态k的参数一样。对模型参数的这种共享通常被称为“系结”。因此，与代替初始状态k的d_min新状态相关联的观测值由与初始状态相同的参数集合支配。

在图8中图解说明了将持续时间约束结合入HMM结构的全部过程。图的左边是来自于初始隐马尔可夫模型结构的状态之一(如图7所示)。图的右边是这个状态通过结合持续时间约束而被变换成的新结构。特别地，d_min(k)-1个新状态(没有自跃迁循环)以简单地左右方式被连接到初始状态，并且这些状态的观测值概率模型参数“被绑定到”初始状态的相应参数。此后的过程如图中的虚线框所示。

初始的隐马尔可夫模型(即P、B1、QRS、T等等)中的每个初始状态重复上述过程。因此将持续时间约束结合入模型结构的全部作用是，用一个新的“约束持续时间的”结构取代模型中每个初始的单一状态。图9显示了作为结果的约束持续时间的HMM结构的一部分。模型中每个波形特征曲线(如P波、QRS复合波等等)此时与多个状态相关联。然后该多个状态为每个ECG波形特征曲线压缩到最小持续时间要求。

一旦根据先前描述的方式将持续时间约束结合入HMM，那么该模型可以通过使用标准维特比算法来被用于分割新的ECG波形。然后，作为结果的分割保证服从持续时间约束，因为该模型不再能生成持续时间小于允许持续时间的状态序列。

图3和图5显示了与之前在图2和图4中考虑的一样的ECG波形。在图3和图5中用实心垂直线表示由内行的ECG分析专家测定的波形边界，以及根据本发明的一个实施例如虚线垂直线所示，应用持续时间约束的隐马尔可夫模型自动测定的状态边界。不仅在模型和专家间有好的一致，而且传统模型的不可靠的“双搏”分割，如图2和4所示，已经被除去，使得作为结果的分析更鲁棒。边界可以进行进一步分析，再次自动地，如测定特定状态的持续时间，如QT间隔，以及可以对这种包括在一系列脉冲的信号之上的间隔的变化进行统计分析。以原始边界信息或进一步分析结果的形式，输出结果，例如通过在计算机系统的显示器上显示。

根据本发明的第三具体实施例，除了对ECG波形的精确分割，还可能从训练好的模型中产生“置信度”(或得分)。这个置信度对该模型在作为结果的分割中的可信度进行量化，并且可以用于评价已分割的ECG波形的质量。因此利用这个量度，有可能自动探测那些过于嘈杂或过于“异常”(在波形形态方面)的ECG波形，以提供可靠的分割及因此对QT间隔的可靠估计。该置信度可以与任何以概率方式计算的分割算法一起使用。

置信度可以定义为给定模型信号的对数概率，即logp(O₁O₂…O_T|λ)，或(给定模型)信号和最概然状态序列s₁s₂…s_T的联合对数概率，即logp(O₁O₂…O_T，s₁s₂…s_T|λ)。相对于原始概率，对数概率被用于将作为结果的值保存在计算机的动态值域里。

可以高效地计算出置信度的精确度，用于隐马尔可夫模型和这些模型的更多的改进形式，如隐半马尔可夫模型(见L.Thoraval、G.Carrault和F.Mora“用于ECG分割的连续变量持续时间HMMs”，IEEE-EMBS学报，1992年，529-530页)和阶乘隐马尔可夫模型(见Z.Ghahramani和M.I.Jordan，“阶乘隐马尔可夫模型”，机器学习，第29卷，第2期，1997年11月)。第一种形式，被定义为给定模型信号的对数概率，可以从下面给出的标准HMM前向变量中计算出来：

α_t(i)＝p(O₁O₂…O_t，s_t＝i|λ)

为了在计算机的动态值域里保存作为结果的值，通常在每个时间步长换算这个变量。这通过用下面给出的比例因子除该前向变量来实现：

$c_t=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\α}}_t\left(i\right)$

于是可以简单地计算出给定模型信号的对数概率：

$\log p(O_1{O}_2...O_T|\mathrm{\λ})=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\log c_t$

置信度的第二种形式，被定义为信号和最概然状态序列(给定模型)的联合对数概率，可以从变量增量δ_t(i)中计算出，如之前在维特比算法背景所述。通常还在每个时间步长换算这个变量，利用如下给出的比例因子：

$d_t=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\δ}}_t\left(i\right)$

于是(给定模型)信号和最概然状态序列的对数概率，可以简单地被计算出：

$\log p(O_1{O}_2...O_T,s_1{s}_2...s_T|\mathrm{\λ})=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\log d_t$

一个可替代使用上述换算过程的方法是，直接计算logα_t(i)值或logδ_t(i)值，如在L R Rabine，“隐马尔可夫模型指南及在语音识别中的精选应用”，IEEE学报，77：257-286，1989中描述的。在这种情况下，在递归过程中直接计算出置信度，且不需要进一步处理。

在许多情况下，给定的ECG信号会包含许多所关心的ECG特征曲线(如多个ECG搏动或多个QT间隔)，并且因此最好为信号中每个所关心的单独特征曲线计算各自的置信度(而不是全部信号的总度量)。对于在时间采样点tl开始且在时间采样点t2结束的特定特征曲线，该置信度可以简单地通过从直到时间t2内信号部分的置信度中减去直到(但不包括)时间t1内信号部分的置信度来计算，即：

logp(O₁O₂…O_t2|λ)-logp(O₁O₂…O_t1-1|λ)

或使用第二种定义的置信度：

logp(O₁O₂…O_t2，s₁s₂…s_t2|λ)-logp(O₁O₂…O_t1-1，s₁s₂s_t1-1|λ)

当在实践中使用这个置信度来评估模型在其特定特征曲线分割上所具有的置信程度时，在分析中必须把所讨论的特征曲线的长度考虑进去。这是必要的，因为两种形式的置信度都随特定ECG特征曲线的长度而改变。

将ECG特征曲线长度结合入置信度分析的最有效的方法是，制作一个置信度(y轴)对特征曲线长度(x轴)的特征曲线坐标图，用于很多规则的ECG特征曲线(即这些特征曲线不被认为是异常的，也不被噪声干扰)。图10显示了这样一个特征曲线坐标图，其中所关心的特定特征曲线被看成一个单独的ECG波形(即ECG信号的该部分对应于一个完整的心搏)。

给定置信度对ECG特征曲线长度的特征曲线坐标图，那么通过使用标准密度建模技术，有可能识别出该二维空间中具有高置信度的区域。特别地，对应于规则(即高置信度)ECG特征曲线的二维数据点集可以通过任何适当的密度函数来建模，如高斯分布或概率线性回归模型。

该给定密度函数会自然地为这个二维数据空间中的高置信度区域定义一个边界。在高斯分布的情况下，这个边界可以被看成位于平均数的两标准偏差间的数据空间区域。在概率线性回归模型的情况下，这个边界可以被看成，在一个合适的统计意义水平上，位于平均回归线的置信度下限之上的数据空间区域(见R.J.Freund和W.J，Wilson“回归分析”，学术出版社，1998年)。图10中用虚线表示了这个下限(在99％的意义水平上)。

给定许多新的ECG特征曲线分割的置信度，和每个被分割特征曲线的长度，那么测定是否有任何置信度落到高置信度区域边界以外就是件很简单的事了。然后可以将那些置信度落到该区域以外的特征曲线的分割结果从最终分析中除去，或可选地将这些特征曲线突出显示留待训练有素的ECG分析员做进一步的人工分析。

图11显示了一部分ECG信号，它包含有两个被肌肉伪噪声干扰的ECG波形。该图还显示了每个通过带有嵌入持续时间约束的隐马尔可夫模型波形的分割。这两个波形在x轴的下面用“A”和“B”标记。

不幸地的是，肌肉伪噪声具有与标准ECG信号类似的带宽，并且因此不能通过简单的过滤技术将其除去。噪声掩盖了下面的ECG信号波形特征曲线，并且使得精确分割波形非常困难。由于这个原因，能够探测到这样的不可靠分割是很重要的。

图12显示了置信度对波形长度的特征曲线坐标图，对于两个嘈杂的ECG波形。还显示出了99％置信度下限(虚线)。两个嘈杂波形的置信度得分都低于这条线，因此使得这些波形能在信号自动分析时被自动探测到。

已经在ECG信号分析背景中描述了本发明的上述具体实施例。尽管如此，本发明不限于这种应用，并且可以被用于分析其他生物医学信号。适合的生物医学信号的例子包括：脑电图(EEG)信号，其中该分析可以用于窒息探测；肌电图(EMG)信号，测量肌肉中的电功率；侵入获得的心脏信号，如在心脏电生理学中，将导管插入心脏，每个导管包含许多电极；以及心脏内电描记图信号，在植入心脏起博器之前使用，包括通过将电极插入身体获得的来自心脏的信号。约束持续时间的维特比算法和置信度测量技术不限于和隐马尔可夫模型一起使用。例如，约束持续时间的维特比算法可以与有限状态离散时间马尔可夫过程一起使用，并且置信度测量可以与任何以概率方式计算的分割算法一起使用。约束持续时间的维特比算法可以与非生物医学信号一起使用，如在检/纠错码领域。

本发明的具体实施例由计算机系统上执行的计算机程序来完成。该计算机系统可以是任何类型的计算机系统，但典型地是传统的个人电脑，它可以执行用任何适当语言编写的计算机程序。该计算机程序可以存储在一计算机可读介质中，该介质可以是任何类型，例如：记录介质，如可插入计算机系统驱动器的圆盘状介质，它可以采用磁性、光学或磁光方式存储信息；计算机系统的固定记录介质如硬盘驱动器；或者一固态计算机存储器。要分析的信号可以直接输入计算机系统，例如从一个ECG装置，或者计算机系统可以从之前获得信号的存储器中读取表示信号的信息。

Claims (16)

1、一种采用隐马尔可夫模型分割生物医学信号的计算机执行方法，该模型包括多个状态，该方法包括如下步骤：

为至少一个状态指定一最小持续时间约束d_min；

对于具有指定最小持续时间的模型中的每个状态，用一组亚态替换该状态，亚态的总数等于最小持续时间约束d_min的值；

将该组亚态连接在一起，形成左右方向的马尔可夫链，其中第一d_min-1亚态中的每个各有为0的自跃迁概率，为1的跃迁到紧接它们右边的状态的跃迁概率，和为0的跃迁到模型中任何其他状态的跃迁概率；以及

将该模型应用于表示生物医学信号的数据上，以获得关于将信号分割成状态的信息。

2、根据权利要求1所述的方法，其中该链中最后的亚态具有与已经被替换的初始状态一样的跃迁概率。

3、根据权利要求1或2所述的方法，其中所有d_min亚态的观测值概率模型被连在一起，使得马尔可夫链中的d_min状态中任何一个的观测值模型的参数与链中任何其他状态的相应参数相等。

4、根据权利要求3所述的方法，其中马尔可夫链中每个亚态的观测值模型的参数被设置成与已经被替换的初始状态的观测值模型的参数相等。

5、根据前述权利要求任一项所述的方法，其中应用模型的步骤包括将模型和要分割的信号提供给维特比算法。

6、一种计算机执行方法，用于通过使用对维特比算法的修正，将包括一列观测值的信号分割为有限状态离散时间马尔可夫过程的一列状态，该修正包括：

对于将结合入维特比算法的有限状态离散时间马尔可夫过程的至少一个状态，定义一个持续时间变量和持续时间约束，该持续时间约束指定所述至少一个状态的最小持续时间；

将修正后的维特比算法应用于信号，以计算构成观测值序列的最概然约束持续时间的状态序列；其中，

在为有限状态离散时间马尔可夫过程的每个状态计算最概然状态序列的每个时间步长，该马尔可夫过程构成了直到那个时间步长内的观测值序列，并以所述状态结束：

对于具有持续时间约束的每个状态，利用那个状态的持续时间变量跟踪先行状态连续序列的长度，该先行状态仅由那个状态组成并在前一时间步长以那个状态结束；

如果那个状态的持续时间变量大于或等于对那个状态指定的持续时间约束，那么在给定时间步长的状态序列计算中，考虑从那个状态到马尔可夫过程中的任何其他给定状态的跃迁；

如果那个状态的持续时间约束小于对那个状态指定的持续时间约束，那么在给定时间步长的状态序列计算中，不考虑从那个状态到马尔可夫过程中任何其他状态的跃迁；并且

在计算出直到给定时间步长的那组最概然状态序列之后，更新具有持续时间约束的每个状态的持续时间变量，用以跟踪先行状态连续序列的长度，该先行状态仅由那个特定状态组成，并在所考虑的时间步长以那个状态结束。

7、根据权利要求6所述的方法，其中为有限状态离散时间马尔可夫过程的每个状态提供了一持续时间变量和持续时间约束。

8、根据权利要求6或7所述的方法，其中该有限状态离散时间马尔可夫过程包括一隐马尔可夫模型。

9、一种用于分析已根据概率分割算法分割的信号的计算机执行方法，该方法包括：

计算每个已分割信号特征曲线的置信度；

画出该置信度对各自信号特征曲线长度的坐标图；

将密度建模技术应用于测定数据空间中与高置信度特征曲线相关联的适当区域；

测定一特定信号特征曲线的置信度是否落到该区域以外。

10、根据权利要求9所述的方法，其中该信号特征曲线包括一更长信号的一部分。

11、根据权利要求9或10所述的方法，其中该置信度是从以下组中选取的一个，该组包括：给定模型的信号特征曲线的对数概率；以及给定模型的信号特征曲线和最概然状态序列的联合对数概率。

12、根据前述权利要求任一项所述的方法，其中该信号是从以下组中选择的生物信号，该组包括：心电图信号、脑电图信号、肌电图信号、心脏电生理学信号和心脏内电描记图。

13、根据前述权利要求任一项所述的方法，其中进一步包括显示有关该信号分析结果的信息。

14、根据前述权利要求任一项所述的方法，其中该信号是心电图信号，并且从该信号分割中获得有关QT间隔的信息。

15、一种计算机程序，能够在计算机系统上运行，并且当这样运行时能够使计算机系统执行根据前述权利要求任一项所述的方法。

16、一种计算机可读介质，存储根据权利要求15所述的计算机程序。

Images