CN105488351A - 一种生成移动心电信号的噪声模型的方法 - Google Patents

Abstract

本发明公开了一种生成移动心电信号的噪声模型的方法，包括：使用时变系数替换传统AR模型序列中的自回归系数，形成一个多阶时变系数AR模型，即得到多阶时变系数AR模型序列，该序列包括未知的时变系数及阶数；根据得到的多阶时变系数AR模型序列，使用时变系数和序列的过去值来得到预测序列的表达式；使用基函数的方法估计多阶时变系数AR模型序列的时变系数；使用基于遗忘因子法确定非平稳随机信号的多阶时变系数AR模型序列的阶数；根据时变系数及阶数最终得到时变系数替换传统AR模型序列。本发明相对于假设序列在短时上是平稳的参数估计等方法，该发明直接使用时变AR模型来模拟心电噪声源模型，更加合理，并且在构造中误差小，致使结果更加精确。

Description

一种生成移动心电信号的噪声模型的方法

技术领域

本发明涉及信号处理的领域，具体涉及一种生成移动心电信号的噪声模型的方法。

背景技术

心电信号的噪声模型可以用于构建真实的心电信号。在心电信号研究中往往需要构建心电的仿真模型。通常心电信号的仿真是在纯净的心电信号中掺加各种噪声信号。仿真的模型越能逼真地模拟真实心电信号，相应的研究结果也越准确。

真实心电信号通常包含基线漂移、工频、肌电干扰和电极移动等各种噪声，所以要求在仿真的心电信号中掺入各种噪声，从而达到接近真实心电的目的。虽然已经有一些模型来模拟噪声，但是能够准确、高效模拟的模型却不多。

目前常用的有两种心电噪声源模型。第一种是使用不同频带的色噪声模拟各种心电噪声源，利用这些仿真噪声源生成噪声，在不同信噪比下掺入到干净仿真心电信号中形成仿真心电模型。第二种常见的噪声源是利用MIT-BIHNon-stresstestDatabase(NSTDB)数据库提供的时长3.5小时的心电信号典型噪声。

第一种模型能够根据研究需要产生足够时长且形态较多的噪声信号，但是平稳的色噪声不能够准确模拟心电信号噪声源，产生的仿真噪声与真实噪声在统计特性与波形形态上相差较大，此模型的精度较低。第二种模型提供的噪声源是真实心电噪声，应用精度高。对于一些特殊的心电信号，NSTDB数据库仅提供三种噪声，包括基线漂移(baselinewander，BW)，肌电干扰(muscleartifact，MA)，电极移动(electrodemotionartifact，EM)，不足以模拟心电信号研究中的噪声源，另外，这个模型的波形变化固定时长无法满足实时移动心电信号要求的更长时间要求。因此当前没有一个可以提供长时间并且精度较高的模型。

针对于对非平稳时间序列建模，当前有两种方法，来构建噪声源模型。传统方法是将非平稳时间序列通过差分变换等方法转换为平稳序列之后，再使用AR模型建模。但是在实际使用中，并不是所有的非平稳序列都能通过差分转换为平稳序列。另一种方法是利用非平稳信号在短时分段内信号是平稳随机的性质，用分段时不变参数AR模型法对其进行建模。这种短时持续时间对多数现象均比较适宜，但是对于心电信号，使用此方法建模时间对于此类现象来说过长，因此达不到理想结果。

发明内容

为解决现有技术存在的不足，本发明公开了一种生成移动心电信号的噪声模型的方法，该模型的精度高并且可以生成长时间的心电噪声信号。

为实现上述目的，本发明的具体方案如下：

一种生成移动心电信号的噪声模型的方法，包括以下步骤：

步骤一：使用时变系数替换传统AR模型序列中的自回归系数，形成一个多阶时变系数AR模型，即得到多阶时变系数AR模型序列，该序列包括未知的时变系数及阶数；

步骤二：根据步骤一得到的多阶时变系数AR模型序列，使用时变系数和序列的过去值来得到预测序列的表达式；

步骤三：使用基函数的方法估计多阶时变系数AR模型序列的时变系数；

步骤四：使用基于遗忘因子法确定非平稳随机信号的多阶时变系数AR模型序列的阶数；

步骤五：通过比较多阶时变系数AR模型序列的序列值与步骤二中的序列的预测值之间的差值来确定基函数的维数，根据该维数进一步确定多阶时变系数AR模型序列的时变系数，根据时变系数及阶数最终得到时变系数替换传统AR模型序列。

进一步的，步骤一中，对于一个时间序列x(n)，阶数为p的传统AR模型是：

x(n)＝-a₁x(n-1)-a₂x(n-2)-...-a_px(n-p)+v(n)(1)

(1)式中v(n)是均值为0，方差为σ²的平稳白噪声残差序列，p是模型的阶数，a_i是AR模型的自回归系数，i＝1,…,p，将上式中的系数a_i替换为时变系数a_i(n),i＝1,…,p，则形成p阶时变系数AR模型，替换过程为：

$\begin{array}{c}x\left(n\right)=-a_1\left(n\right)x(n-1)-a_2\left(n\right)x(n-2)-\mathrm{...}-a_p\left(n\right)x(n-p)+v\left(n\right)\\ =-\[x(n-1),x(n-2),\mathrm{...},x(n-p)\]\left[\begin{array}{c}{a}_1\left(n\right)\\ a_2\left(n\right)\\ .\\ .\\ .\\ a_p\left(n\right)\end{array}\right]+v\left(n\right)\\ =-\underset{i=1}{\overset{p}{\Σ}}{a}_i\left(n\right)x(n-i)+v\left(n\right)\end{array}---\left(2\right)$

v(n)仍然是均值为0，方差为σ²的平稳白噪声，p是模型的阶数，a_i(n),i＝1,…,p是模型的时变系数，表示的意义是第n时刻的模型参数。

进一步的，步骤二中，使用模型的时变系数和序列的过去值预测序列表达式为：

$\hat{x}\left(n\right)=-\underset{i=1}{\overset{p}{\Σ}}{a}_i\left(n\right)x(n-i)---\left(3\right)$

进一步的，步骤三中，假设时变系数a_i(n)，i＝1,…,p是一组基函数的线性组合，表达式如下：

$a_i\left(n\right)=\underset{j=0}{\overset{m}{\Σ}}{a}_{ij}{g}_j\left(n\right)---\left(4\right)$

式中，a_ij为使用最小二乘估计方法估计由时变系数变换得到的常系数，为矩阵系数A的元素，{g_j(n)，j＝0,…,m}是一组基时间函数，m为基函数的维数的总数，并将这一组基函数代入序列的预测表达式(3)中，并且令

A＝[a₁₀,a₁₁,...,a_1m|...|a_p0,a_p1,...,a_pm]^T(5)

$X_n^T=\[x(n-1)g_0\left(t\right),...,x(n-1)g_m\left(n\right)|...|x(n-p)g_0\left(n\right),...,x(n-p)g_m\left(n\right)\]---\left(6\right)$

其中，A为系数矩阵，序列表达式变为：

$x\left(n\right)=-X_n^{T}A+v_n---\left(7\right)$

进一步的，使用傅里叶基数来表示时变信号，其中，N是观测时间内信号采样数，ω＝1/N，j表示基函数的当前的维数或个数，j＝0,1,2,…,m，k＝1,2,…,N；

进一步的，模型阶数为12，即p＝12。

进一步的，确定基函数维数时，选择维数m的取值范围从1到8，原则是通过真实信号与拟合信号之间误差最小的方法来确定基函数维数，即根据实际值与估计值的最小二乘误差进行确定。

进一步的，使用最小二乘估计方法估计由时变系数变换得到的常系数矩阵A，并且利用估计残差的方差最小的优化准则，由式(7)可得到A的最小二乘估计式：

${\hat{A}}_N={\left({\mathrm{\φ}}^T\mathrm{\φ}\right)}^{-1}{\mathrm{\φ}}^{T}Y---\left(9\right)$

式中Y^T＝[x₁…x_N]，1,…,N是采样点时刻，即采样数，估计出A后，再根据式(4)即可求得时变系数a_i(n)在各个时刻的值。

本发明的有益效果：

相对于假设序列在短时上是平稳的参数估计等方法，本发明直接使用时变AR模型来模拟心电噪声源模型，更加合理，结果更加精确，生成的噪声源的数学统计特性和波形变化与NSTDB数据库中提供的有限样本噪声序列基本一致，具有很高的精确度。本发明不再局限于模拟高频噪声、工频噪声、低频噪声等平稳噪声序列，有很大的应用范围。

NSTDB数据库可以提供精确的噪声源，但是时间只有3.5小时，并且对于移动心电等信号，数据库中三个噪声源无法模拟，该发明解决了数据库中样本噪声序列不足的问题，可以提供长时间、复杂的心电噪声源。

附图说明

图1是构造时变AR模型的流程图。

图2a-图2b分别是原始BM与时变AR模型预测BM。

图3a-图3b分别是原始MA与时变AR模型预测MA。

图4a-图4b分别是原始EM与时变AR模型预测EM。

具体实施方式：

下面结合附图对本发明进行详细说明：

如图1所示，为构造时变AR模型的流程图，具体步骤为：

1、从MIT数据库中下载BW、MA、EM三种噪声模型。

2、利用AR模型建立跟真实心电信号相似的心电噪声源模型。

3、对比使用时变模型模拟出来的噪声源模型和MIT数据库中真实心电信号噪声源两者的误差，来判断模型的效果。

其中步骤2以传统AR模型为基础建立心电噪声源模型主要包括以下几个步骤：

1)使用时变系数替换传统AR模型序列中的自回归系数，形成一个多阶时变系数AR模型。在此过程中，要使用序列过去值的线性组合来拟合序列的当前值。

2)使用时变系数和序列的过去值来预测序列的表达式。此时，将对非平稳随机信号建立时变AR模型的过程转化成对该时变模型进行定阶和求时变系数的过程。

3)使用基函数的方法估计模型的时变系数。当时变系数被看作是一个时间函数时，就将时变系数的估计问题就变成常参数的估计。模型建立之后，使用最小二乘法估计由时变系数变换得到的常系数矩阵。

4)使用基于遗忘因子法确定非平稳随机信号的时变AR模型的阶数。

5)通过比较序列值与序列的预测值之间的差值来确定基函数的维数，使模型具有最好的精度和效率。

算例说明：

为了更好的描述本发明，下面给出了更为详细的实施例子：

1、从MIT数据库中下载BW、MA、EM三种噪声模型。

2、利用AR模型建立跟真实心电信号相似的心电噪声源模型，主要包括以下几个步骤：

1)对于一个时间序列x(n)，阶数为p的传统AR模型是：

x(n)＝-a₁x(n-1)-a₂x(n-2)-...-a_px(n-p)+v(n)(1)

(1)式中v(n)是均值为0，方差为σ²的平稳白噪声残差序列，p是模型的阶数，a_i(i＝1,…,p)是AR模型的自回归系数。将上式中的系数a_i替换为时变系数a_i(n),i＝1,…,p，则形成p阶时变系数AR模型，替换过程为：

$\begin{array}{c}x\left(n\right)=-a_1\left(n\right)x(n-1)-a_2\left(n\right)x(n-2)-\mathrm{...}-a_p\left(n\right)x(n-p)+v\left(n\right)\\ =-\[x(n-1),x(n-2),\mathrm{...},x(n-p)\]\left[\begin{array}{c}{a}_1\left(n\right)\\ a_2\left(n\right)\\ .\\ .\\ .\\ a_p\left(n\right)\end{array}\right]+v\left(n\right)\\ =-\underset{i=1}{\overset{p}{\Σ}}{a}_i\left(n\right)x(n-i)+v\left(n\right)\end{array}---\left(2\right)$

v(n)仍然是均值为0，方差为σ²的平稳白噪声，p是模型的阶数。a_i(n),i＝1,…,p是模型的时变系数，表示的意义是第n时刻的模型参数。

2)根据(2)式，使用模型的时变系数和序列的过去值预测序列表达式为：

$\hat{x}\left(n\right)=-\underset{i=1}{\overset{p}{\Σ}}{a}_i\left(n\right)x(n-i)---\left(3\right)$

3)假设时变系数a_i(n)，i＝1,…,p是一组基函数的线性组合，表达式如下：

$a_i\left(n\right)=\underset{j=0}{\overset{m}{\Σ}}{a}_{ij}{g}_j\left(n\right)---\left(4\right)$

式中{g_j(n)，j＝0,…,m}是一组基时间函数，m为基函数的维数的总数。并将这一组基函数代入序列的预测表达式中，并且令

A＝[a₁₀,a₁₁,...,a_1m|...|a_p0,a_p1,...,a_pm]^T(5)

$X_n^T=\[x(n-1)g_0\left(t\right),...,x(n-1)g_m\left(n\right)|...|x(n-p)g_0\left(n\right),...,x(n-p)g_m\left(n\right)\]---\left(6\right)$

其中，A为系数矩阵，序列表达式变为：

$x\left(n\right)=-X_n^{T}A+v_n---\left(7\right)$

在本发明中使用傅里叶基数来表示时变信号，其中，N是观测时间内信号采样数，ω＝1/N,j表示基函数的当前维数(个数)，j＝0,1,2,…,m，k＝1,2,…,N。

4)参考EURASIPJournalonAppliedSignalProcessing2007中第94到104页里的文章MultichannelECGandnoisemodeling:applicationtomaternalandfetalECGsignals中的结论，当模型的阶数在12-16之间时，足以满足其构建心电信号噪声的时变AR模型，因此本发明确定模型阶数为12，即p＝12。

5)确定基函数维数时，选择维数m的取值范围从1到8，原则是通过真实信号与拟合信号之间误差最小的方法来确定基函数维数。根据实际值与估计值的最小二乘误差进行确定，其结果如表1所示。从表中可以看出，当基函数的维数从1到8变化时，随着维数的增大，误差逐渐减小。当m的值达到一定程度时，误差减少的幅度越来越少，预测精度逐渐趋于稳定，因此，选择的维数m没必要过大。选择维数同时也要考虑计算效率的影响。因此本发明在此选择维数为5，即m＝5。

表1

将确定的维数和阶数代入模型中，此时的AR模型转换成可以反映心电信号噪声的模型。

3、对比使用时变模型模拟出来的噪声源模型和MIT数据库中真实心电信号噪声源两者的误差，来判断模型的效果。

使用模型模拟的BW、MA、EM与数据库中真实的心电噪声源对比如图2a-图2b、图3a-图3b、图4a图4b所示。图2a、图3a及图4a分别是原有的MIT数数据库中真实的噪声源，图2b、图3b及图4b分别是用本发明模型拟合的噪声源，可以看出，使用时变模型构造出来的噪声源模型其波形和动态特征与MIT数据库中真实的信号差异不大，精度很高。因此时变AR模型对于非平稳信号具有很好的拟合性，能贴近实际信号拟合出仿真信号，具有很好的实际应用价值。

使用以基函数为基础建立的时变信号AR模型，不仅可以很好的反映人体心电信号的噪声的动态特征和波形变化，大大提高了心电噪声源模型的精确度，而且可以提供多种长时间噪声，不再局限于数据库中。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (8)

1.一种生成移动心电信号的噪声模型的方法，其特征是，包括以下步骤：

步骤一：使用时变系数替换传统AR模型序列中的自回归系数，形成一个多阶时变系数AR模型，即得到多阶时变系数AR模型序列，该序列包括未知的时变系数及阶数；

步骤二：根据步骤一得到的多阶时变系数AR模型序列，使用时变系数和序列的过去值来得到预测序列的表达式；

步骤三：使用基函数的方法估计多阶时变系数AR模型序列的时变系数；

步骤四：使用基于遗忘因子法确定非平稳随机信号的多阶时变系数AR模型序列的阶数；

步骤五：通过比较多阶时变系数AR模型序列的序列值与步骤二中的序列的预测值之间的差值来确定基函数的维数，根据该维数进一步确定多阶时变系数AR模型序列的时变系数，根据时变系数及阶数最终得到时变系数替换传统AR模型序列。

2.如权利要求1所述的一种生成移动心电信号的噪声模型的方法，其特征是，步骤一中，对于一个时间序列x(n)，阶数为p的传统AR模型是：

x(n)＝-a₁x(n-1)-a₂x(n-2)-...-a_px(n-p)+v(n)(1)(1)式中v(n)是均值为0，方差为σ²的平稳白噪声残差序列，p是模型的阶数，a_i是AR模型的自回归系数，i＝1,…,p，将上式中的系数a_i替换为时变系数a_i(n),i＝1,…,p，则形成p阶时变系数AR模型，替换过程为：

$\begin{array}{c}x\left(n\right)=-a_1\left(n\right)x(n-1)-a_2\left(n\right)x(n-2)-\mathrm{...}-a_p\left(n\right)x(n-p)+v\left(n\right)\\ =-\[x(n-1),x(n-2),\mathrm{...},x(n-p)\]\left[\begin{array}{c}{a}_1\left(n\right)\\ a_2\left(n\right)\\ .\\ .\\ .\\ a_p\left(p\right)\end{array}\right]+v\left(n\right)\\ =-\underset{i=1}{\overset{p}{\Σ}}{a}_i\left(n\right)x(n-i)+v\left(n\right)\end{array}---\left(2\right)$

v(n)仍然是均值为0，方差为σ²的平稳白噪声，p是模型的阶数，a_i(n),i＝1,…,p是模型的时变系数，表示的意义是第n时刻的模型参数。

3.如权利要求1所述的一种生成移动心电信号的噪声模型的方法，其特征是，步骤二中，使用模型的时变系数和序列的过去值预测序列表达式为：

$\hat{x}\left(n\right)=-\underset{i=1}{\overset{p}{\Σ}}{a}_i\left(n\right)x(n-i)---\left(3\right).$

4.如权利要求1所述的一种生成移动心电信号的噪声模型的方法，其特征是，步骤三中，假设时变系数a_i(n)，i＝1,…,p是一组基函数的线性组合，表达式如下：

$a_i\left(n\right)=\underset{j=0}{\overset{m}{\Σ}}{a}_{ij}{g}_j\left(n\right)---\left(4\right)$

式中，a_ij为使用最小二乘估计方法估计由时变系数变换得到的常系数，为矩阵系数A的元素，{g_j(n)，j＝0,…,m}是一组基时间函数，m为基函数的维数的总数，并将这一组基函数代入序列的预测表达式(3)中，并且令

A＝[a₁₀,a₁₁,...,a_1m|...|a_p0,a_p1,...,a_pm]^T(5)

$X_n^T=\[x(n-1)g_0\left(t\right),\mathrm{...},x(n-1)g_m\left(n\right)\left|\mathrm{...}\right|x(n-p)g_0\left(n\right),\mathrm{...},x(n-p)g_m\left(n\right)\]---\left(6\right)$ 其中，A为系数矩阵，序列表达式变为：

$x\left(n\right)=-X_n^{T}A+v_n---\left(7\right).$

5.如权利要求4所述的一种生成移动心电信号的噪声模型的方法，其特征是，使用傅里叶基数来表示时变信号，其中，N是观测时间内信号采样数，ω＝1/N，j表示基函数的当前的维数或个数，k＝1,2,…,N；

6.如权利要求1所述的一种生成移动心电信号的噪声模型的方法，其特征是，模型阶数为12，即p＝12。

7.如权利要求1所述的一种生成移动心电信号的噪声模型的方法，其特征是，确定基函数维数时，选择维数m的取值范围从1到8，原则是通过真实信号与拟合信号之间误差最小的方法来确定基函数维数，即根据实际值与估计值的最小二乘误差进行确定。

8.如权利要求4所述的一种生成移动心电信号的噪声模型的方法，其特征是，使用最小二乘估计方法估计由时变系数变换得到的常系数矩阵A，并且利用估计残差的方差最小的优化准则，由式(7)可得到A的最小二乘估计式：

${\hat{A}}_N={\left({\mathrm{\φ}}^T\mathrm{\φ}\right)}^{-1}{\mathrm{\φ}}^{T}Y---\left(9\right)$

式中1,…,N是采样点时刻，即采样数，估计出A后，再根据式(4)即可求得时变系数a_i(n)在各个时刻的值。