仿射变换的二阶微分方程生命信号滤波算法

周游， 王鹏飞， 査富生， 郭伟， 李满天

(哈尔滨工业大学 机器人技术与系统国家重点实验室，黑龙江哈尔滨 150001)



仿射变换的二阶微分方程生命信号滤波算法

周游， 王鹏飞， 査富生， 郭伟， 李满天

(哈尔滨工业大学 机器人技术与系统国家重点实验室，黑龙江哈尔滨 150001)

由于微小型嵌入式系统尺寸和计算能力的限制，在对生命体征信号进行检查时，往往较难实现实时与精确兼顾。为了解决该问题，本文提出了一种基于二阶微分方程与仿射变换理论结合的心电信号滤波算法(ATSDEF)，并进行了与粒子滤波方法的滤波性能对比实验，证明了ATSDEF方法有较粒子滤波更好的滤波性能。通过微小型嵌入式心电测量系统对正常与异常心电信号进行在线滤波实验，验证了该算法可以有效的滤除心电信号的噪声，且能够保持心电信号中较为微弱的细节特征，从而保证了信号滤波后的精度，解决了微小型嵌入式系统的实时信号处理难题。

心电信号；微小型嵌入式系统；二阶微分方程；龙格库塔；仿射变换；滤波

生命特征信号是一种不稳定、非线性且具有概率性的信号，因此对于人体疾病监测来说，尽可能地准确滤除噪声对提取生命体征信号特征以及心脏疾病诊断显得尤为重要。特别是在人体自身以及外界随机扰动下，生命体征信号可能会出现生命信号异常、信号基线发生漂移，甚至是信号被噪声掩埋。此时如何滤除上述干扰，仍是信号处理领域的难题，特别是针对微小型、嵌入式、多任务系统，由于该类系统计算资源的限制，较难运行计算量较大的滤波算法[1-2]。

为了解决该难题，研究人员做了大量研究，并提出了众多滤波方法，如粒子滤波法[3]、变换域自适应滤波法[4]、自适应小波滤波法[5-7]、神经网络与模糊滤波法[8]等，但都没有取满意效果。尤其是近年来对微弱信号处理以及实时性要求的提高(如微纳制造、人体生命特征信息检测、超远距离通讯等等)，导致传统滤波方法解决这类情况较难获得满意的效果。本文基于二阶微分方程本身所固有的频率响应能力，结合仿射变换理论，提出了一种可处理较微弱心电信号的滤波算法(affine transformation of second order differential equation filtering，ATSDEF)，通过与粒子滤波方法的对比来说明ATSDEF方法的性能，最后，利用心电信号实时测量实验系统对正常和异常心电信号滤波的实验验证该算法的有效性。

1 二阶微分方程算法

一个典型的二阶齐次微分方程有如下频率响应形式：

(1)

式中 ω和ξ分别为二阶微分方程的固有频率和阻尼率。

由式(1)可以看出，由于二阶微分方程是众多物理现象(如机械振动、粒子振动、RLC振荡电路等等)的归一化方程，因此，其具有明显的振荡输出，即具有频率响应特性的，因此，从理论上看，式(1)应具有一定的滤波能力。

为了验证式(1)是否具有滤波能力，进行如下实验：利用Matlab模拟一组随机信号f(t)，将该信号作为式(1)的输入，可得

(2)

利用式(2)对一组随机心电信号进行滤波实验，其结果如图1所示。

图1 心电信号滤波结果Fig.1 The results of ECG filtering

从图1中可以看出，无论是式(1)的位置输出x，还是速度输出x′，都不具有信号处理能力。但同时从图1中也可以发现，速度输出x′的频率响应能力明显好于位置输出x，这主要是速度输出本质上是位置输出的微分，因此x′对频率变化更加敏感。

根据振动理论中的共振现象，振动系统的位置输出不仅具有频响能力，而且可以进行信号滤波，如目前应用广泛的机械滤波器[9-11]，而振动系统方程属于典型的二阶微分方程，因此，二阶微分方程也应具有滤波能力。然而，图1所示结果却没有发现振动系统应具有滤波能力。其主要原因是：1)振动系统参数一旦选定，其共振频率也随之固定，即其仅能对与共振频率接近的信号产生较大响应，而对其他频率的信号响应幅值迅速衰减，从而起到共振频率以外信号的抑制作用，达到滤波效果；2)图1中的输入信号为随机信号，其频率并非固定不变，而是随机变化的，而振动系统较难适应如此变化的信号，因此，其难以对该类信号进行滤波。

基于上述分析，可获得如下启发：如将式(2)的速度输出x′对随机信号的频率响应能力赋予位置输出x，将使x也具有类似x′的对随机信号的频响能力。

为了实现将式(2)的速度输出x′对随机信号的频率响应能力赋予位置输出x，则需将x′的变化特性赋予x，因此，将式(2)写成一阶方程组的形式：

(3)

式中：v为速度输出。考虑到存在一随机输入，则式(3)可写为

(4)

式中：λ为表征随机输入强度的系数。由于式(4)中含有随机输入f(t)，根据随机微分方程定义，式(4)是具有典型随机特征的随机微分方程。对于随机微分方程的解法，Milstein等相继给出了基于二阶、四阶龙格库塔的数值解法。然而，与常微分方程不同，随机微分方程中含有不确定的且与时间有关的随机项或随机过程，使得二阶、四阶龙格库塔法所获得的数值解精度不高。针对该问题，王鹏等提出了三阶半隐式随机龙格库塔法，该方法与二阶、四阶龙格库塔法相比具有更高的稳定性和精度[12]。因此，本文中选用该方法对式(4)进行数值求解。

根据文献[12]研究，构造式(4)的三阶半隐式随机龙格库塔数值解表达式为

(5)

其中

(6)

式中：s为截取的随机输入长度，本文中取随机输入的50个连续采样点。令

(7)

则，可得

(8)

将式(5)代入式(7)，可得

(9)

其中

(10)

(11)

将式(9)代入式(8)，可得

(12)

令

(13)

则(11)式变为

(14)

为了将速度输出vn的频率响应能力置换给位置输出xn，根据仿射变换理论，令矩阵D为仿射矩阵:

(15)

利用式(13)右乘N，可得

(16)

对比N与式(16)，可以发现二者的列向量是互换的。

令

(17)

利用式(17)右乘M，可得

(18)

则式(13)可以写为

Vn+1=PVn+Rfn

(19)

其中

(20)

将式(19)写成迭代形式：

Vn+1=PVn+Rfn

Vn+2=PVn+1+Rfn+1=P2Vn+PRfn+Rfn+1

Vn+3=PVn+2+Rfn+2=P3Vn+P2Rfn+PRfn+1+

Rfn+2

⋮

Vn+m=PVn+2+Rfn+2=PmVn+Pm-1Rfn+

Pm-2Rfn+1+…+PRfn+m-2+Rfn+m-1

(21)

式(21)可以写为

(22)

2 实验验证

粒子滤波(partical filter,PF)方法是较为广泛使用的信号滤波方法[13-16]，为了验证所建立的滤波算法的正确性，本文进行了与粒子滤波对比实验，根据文献[13-16]，选择粒子数为200。对于本文所建立的ATSDEF方法，其参数选择为p=0.6，q=1.1，h=0.01，λ=0.8。对比实验结果如图2所示。

图2 ATSDEF与粒子滤波对比实验Fig.2 Compared with particle filter

从图2中可以直观看出，ATSDEF处理结果比粒子滤波处理结果更加接近原始信号。为了进一步比较二者滤波的结果，对图2中的数据进行信噪比、

均方根误差和信噪比增益三个方面计算，如表1所示。

表1 ATSDEF与粒子滤波性能对比

从表1中的数据可以看出，ATSDEF滤波方法在信噪比、均方根误差和信噪比增益三个方面均明显优于粒子滤波。

为了更好的观察和分析所建立的滤波算法，分别对正常(简称Z)和异常(简称Y)心电信号进行滤波实验，实验系统如图3所示。

图3 心电信号实时测量系统Fig.3 ECG real-time measurement system

图3所示系统主要包括电源、手机无线接口、信号处理电路和测量点，从图3中可以看出，实验系统中的信号处理电路很小，面积约为2.5×2.5 cm2。

2.1 正常心电信号滤波

从图4中可以直观的看出，本文所建立的滤波算法对6个心脏功能正常的人的心电信号的滤波效果非常明显，如图4中每个子图的第二行所示。且将原始信号减去滤波后信号，即可将噪声有效的分离，如图4中每个子图的第三行所示。

图4 正常心电信号降噪结果Fig.4 Normal ECG signal filtering results

为了进一步观察滤波效果，将图4中部分滤波结果进行放大，以便观察滤波后的局部细节，如图5所示。图5中横坐标与纵坐标的定义与图4中相同。

图5 图4中部分心电信号滤波结果的局部放大Fig.5 Partial amplification of partial ECG signal in Fig.4

从图5中可以看出，滤波后的心电信号各个波形被完整的保留了，如图5中各个箭头所指示的P波、QRS波群以及T波。即使被噪声几乎完全掩盖的心电信号波形(如图4中(e)Z5号心电信号记录所示)，也能被本文所建立的滤波算法有效滤出，如图5中(c)所示，其中P波被有效的滤出来，这为下一步的心电信号特征提取提供了非常有利的基础。

2.2 异常心电信号滤波

异常心电信号作为人体心脏的功能异常及其活动过程在体表的反映，其正确处理对心血管疾病的诊断和治疗具有重要作用，同时也能为预防心血管疾病提供依据。对6个心脏功能异常的人的心电信号的滤波结果如图6所示。

对比图4可以发现，图6所示的心电信号峰值不具有规律性(如图6中(a)、(d)、(e)和(f)所示)，或心电信号波形异常(如图6中(b)和(c)所示)。但本文所建立的滤波算法都能对上述异常心电信号进行有效的滤波，如图6中各个子图的第二行所示。其具体细节如图7和图8。

图6 异常心电信号的滤波Fig.6 Filtering of abnormal ECG signals

图7 节拍异常的心电信号细节Fig.7 Abnormal heart beat details

图8 幅值异常的心电信号细节Fig.8 Amplitude anomaly ECG signal details

在图7中，每个子图的第一行为原始心电信号，第二行为自适应阈值窗函数降噪方法的滤波结果，横坐标与纵坐标的意义与图6相同。从图6中可以看出，节拍异常的心电信号在某个时间段会出现频率异常，这种异常有时发生在QRS波群或者R波上，如图6中(a)和(d)所示，由于QRS波群或者R波的峰值较大，因此也并不会对滤波噪声很多影响。然而，如果节拍异常发生在Q波时，由于Q波本身幅值相对较小，特别是在某个时间段内连续出现多个Q波时，造成Q波不明显，如图7中各个子图的第一行中虚线箭头所示，将会对滤波算法带来很大的负面影响。从图7中第二行点线箭头所示结果上看，所建立的自适应阈值窗函数降噪方法能有效滤除Q波节拍异常的心电信号噪声，还原心电信号的本来波形。

在图8中，图8所示的为另一种异常心电信号，即波形峰值异常。此类心电信号R波峰值过小，甚至没有R波峰值，一旦遇到噪声，很有可能被噪声淹没，如图8中虚线箭头所示，这使得对此类信号进行滤波变得更加困难。但从图8中点线箭头所示的结果上看，所建立的滤波算法也能有效滤除噪声。

3 结论

1)结合仿射变换理论，将二阶微分方程中速度输出对输入信号频率的敏感性通过仿射变换赋予位置输出，使得二阶微分方程的位置输出因具有了频率响应能力，从而可以获得较好的滤波效果。

2)建立了心电信号滤波算法(ATSDEF)，并进行了与粒子滤波的对比实验，通过对比实验可以看出，在信噪比、均方根误差和信噪比增益等重要的信号处理指标上，本文所建立的滤波算法相较于粒子滤波都有非常明显的提高。

3)利用嵌入式心电信号检测系统对ATSDEF算法进行在线实时心电信号滤波实验，通过心电信号滤波实验验证了该算法可以在保持微弱的细节特征的同时有效的滤除心电信号的噪声，解决了高精度处理难得问题，为后续特征识别与提取奠定了基础。

[1]ZHANG X， HUANG J C， SONG E L， et al. Design of small MEMS microphone array systems for direction finding of outdoors moving vehicles [J]. Sensors， 2014， 14: 4384-4398.

[2]刘富， 于鹏， 刘坤. 采用独立分量分析Zernike矩的遥感图像飞机目标识别[J]. 智能系统学报， 2011， 6(1)： 51-56.

LIU Fu， YU Peng， LIU Kun. Research concerning aircraft recognition of remote sensing images based on ICA Zernike invariantmoments[J]. CAAI Transactions on intelligent systems， 2011， 6(1)： 51-56.

[3]PARIS S， DURAND F. A fast approximation of the bilateral filter using a signal processing approach[J]. Int J comput Vis, 2007， 81: 24-52.

[4]齐林， 周丽晓. 变换域自适应滤波算法的研究[J]. 郑州大学学报：理学版， 2007， 39(1)： 61-66. QI Lin， ZHOU Lixiao. Study on transform-domain adaptive filtering algorithm[J]. J of Zhengzhou Univ Nat.Sci.Ed， 2007， 39(1)： 61-66.

[5]XU Y S， WEAVER J B， HEALY D M， et al. Wavelet transform domain filters: a spatially selective noise filtration technique[J]. IEEE trans. image process， 1994， 3: 747-758.

[6]USMAN I， KHAN A， ALI A， et al. Reversible watermarking based on intelligent coefficient selection and integer wavelet transform[J]. Int J innova-tive comput Inf control， 2009， 5: 4675-4682.

[7]XU R， CHEN Y W. Wavelet-based multiresolution medical image registration strategy combining mutual information with spatial information[J]. Int J Innovative comput Inf control， 2007， 3: 285-296.

[8]PARISS A M. A recurrent fuzzy filter for the analysis of lung sounds[J]. BioSystems， 2007， 87: 164-171.

[9]ADAM M H， SCOTT G， JLHN P H， et al. Miniature low-pass mechanical filter for improved frequency response with MEMS microphones & low-pressure transducers[J]. Sensors and actuators A: Physical， 2014， 210: 51-58.

[10]LUTSENKO V I， KOSAREV O I. A highly sensitive precision piezo accelerometer with a built-in mechanical filter[J]. Journal of machinery manufacture and reliability， 2009， 38(1)： 91-93.

[11]MATTHEW J M， JAMES A S， SIETSE M N. Mechanicalltering by the boundary layer anduid-structure interaction in the supercial neuromast of thesh lateral line system[J]. J comp physiol， A 2008， 194: 795-810.

[12]王鹏， 吕显瑞， 张伸煦. 求随机微分方程的三阶半隐式随机龙哥库塔方法[J]. 吉林大学学报：理学版， 2008， 46(2)： 119-223.

WANG Peng， LU Xianrui， ZHANG Shenxuan. Three-stage semi-implicit stochastic runge-kutta methods for stochastic differential equations[J]. Journalof Jilin University： science edition， 2008， 46(2)： 119-223.

[13]王法胜， 鲁明羽， 赵清杰，等. 粒子滤波算法[J]. 计算机学报， 2014， 37(8)： 1679-1694.

WANG Fasheng， LU Mingyu， ZHAO Qingjie， et al. Particle filtering algorithm[J]. Chinese journal of computers， 2014， 37(8)： 1679-1694.

[14]付巍， 郑宾. 基于改进旳粒子滤波的静电目标跟踪算法. 高技术通讯， 2014， 24(2)： 138-143.

FU Wei， ZHENG Bin. An electrostatic target tracking algorithm based on improved particle filtering[J]. High technology letters， 2014， 24(2)： 138-143.

[15]周帆， 江维， 李树全，等. 基于粒子滤波的移动物体定位和追踪算法[J]. 软件学报， 2013， 24(9)： 2196-2213.

ZHOU Fan， JIANG Wei， LI Shuquan， et al. Moving target localization and tracking algorithms: A particle filter based method[J]. Journal of software， 2013， 24(9)： 2196-2213.

[16]张琪， 乔玉坤， 孔祥玉,等. 随机摄动强跟踪粒子滤波算法[J]. 物理学报， 2014， 63(11)： 110505: 1-9.

ZHANG Qi， QIAO Yukun， KONG Xiangyu, et al. Study on stochastic perturbation strong tracking particlelter[J]. Acta phys sin， 2014， 63(11)： 110505: 1-9.

本文引用格式：

周游， 王鹏飞， 査富生，等. 仿射变换的二阶微分方程生命信号滤波算法[J]. 哈尔滨工程大学学报， 2017， 38(5)： 752-758.

ZHOU You， WANG Pengfei， ZHA Fusheng， et al. Affine transformation-based life signal filtering algorithm for two-order differential equation[J]. Journal of Harbin engineering university， 2017， 38(5)： 752-758.

Affine transformation-based life signal filtering algorithm
for two-order differential equation

ZHOU You， WANG Pengfei， ZHA Fusheng， GUO Wei， LI Mantian

(State Key Laboratory of Robotics and Systems， Harbin Institute of Technology， Harbin 150001， China)

Due to limitations in the size and computational power of miniature embedded systems， it is often difficult to achieve real-time and accurate reconciliation when testing vital signs. To solve this problem， we propose an ECG signal filtering algorithm based on a second-order differential equation and the affine transformation theory (ATSDEF). We compared the filtering performance of our proposed algorithm with that of the particle filter method and the results confirm the ATSDEF to be more effective. Lastly， we validated the algorithm in normal and abnormal ECG filtering experiments， the results of which demonstrate its ability to not only effectively filter ECG signal noise， but also to retain detailed weak characteristics in the ECG signal to ensure the accuracy of signal filtering， thereby solving the longstanding problems of the real-time signal processing of miniature embedded systems.

ECG signal; miniature embedded system; second-order differential equation; Longge Kuta; affine transformation; filter

2016-03-11.

日期：2017-04-26.

国家自然基金项目(61375097；61473105);黑龙江省自然基金项目(F2015008);机器人技术与系统国家重点实验室基金项目(SKLRS201620B；SKLRS201603C).

周游(1983-)，男，博士研究生; 査富生(1975-)，男，副教授，博士生导师.

査富生，E-mail: zfsh751228@163.com.

10.11990/jheu.201603036

TN911.72

A

1006-7043(2017)05-0752-06

网络出版地址：http://www.cnki.net/kcms/detail/23.1390.u.20170426.1237.056.html