局部投影和离散小波变换在心音信号去噪中的应用

梁庆真，郭兴明，袁志会

(重庆大学 生物工程学院，生物流变科学与技术教育部重点实验室，重庆 400044)

心音信号是一种非线性非平稳的微振动信号，在采集的过程中不可避免的会引入干扰噪声，使得部分有用信息淹没在噪声中[1]，因此在对其进行进一步分析之前，心音去噪显得尤为重要。目前，小波去噪在心音去噪中已经得到广泛的应用，它克服了傅里叶分析不能对非平稳信号进行有效分析的缺点，但是心音信号的频带与噪声频带存在部分或全部重叠，所以采用传统的小波分解重构法对心音进行线性去噪时，会相应的去除掉与噪声频带相重叠的部分心音，在这个过程中信号本身不可避免的会损失掉部分非线性特征信息。采用非线性小波变换函数[2]的离散小波阈值去噪法(Discrete Wavelet Threshold,DWT)能较好的抑制噪声，且能得到更精确的重构信号[3]。基于相空间重构技术的局部投影降噪法(Local Projection,LP)[4-5]，它将一维的时间序列拓展到高维相空间，利用系统信号和噪声信号在高维相空间上分布特性的差异来滤除噪声信号，在滤除噪声的同时可以有效保留微弱特征信号成分[6]。本文综合两者的优点，提出噪声水平自适应估计的局部投影与离散小波阈值相结合的去噪方法，该方法既能得到精确的重构信号又能保留信号的非线性特征，采用该方法对Lorenz序列以及实测心音进行降噪研究，用信噪比和均方误差来检验该方法的去噪效果，并提取去噪前后信号的最大Lyapunov指数来说明该方法对信号非线性特征的影响。同时，为了体现本文方法的优势，还将与局部投影和离散小波阈值去噪进行对比分析。

1 基本原理

1.1 噪声水平自适应估计的局部投影降噪

Zn=[Zn,Zn+τ,Zn+2τ,…,Zn+(m-1)τ]

(1)

在理想的无噪声(Wn=0)的情况下，有如下的相空间映射关系：

(2)

当系统存在噪声(Wn≠0)时，式(2)可表示为：

(3)

式(2)在Zn的邻域Nn内可线性化展开为：

(4)

(5)

Nn,

k=[1,n)∪(n,L-(m-1)τ]}

(6)

同理，式(3)在Zn的邻域内线性化展开为：

(7)

其中:ηn为噪声总和。式(7)表明在m维的相空间里，存在一个超平面(由邻域内相点局部线性化得到)，相点与超平面的偏差代表了噪声部分。在这个m维的相空间里，存在着一个吸引子空间和零子空间(由Q个相互正交的向量组成，可记为P=[P1,P2,…,PQ])，由系统动力学行为所决定的吸引子只被局限在吸引子空间中，当系统受到噪声污染时，整个相空间中都会有分量，必定是噪声使得零子空间产生了分量，局部投影降噪算法就是通过向吸引子空间投影的方法来滤除噪声成分。投影后的修正值可表示为：

(8)

1.2 离散小波阈值降噪

带噪声的信号模型可表示为：

(9)

其中:f(i)为真实信号，y(i)为含噪信号，e(i)是噪声。一般情况下，噪声是方差为σ2的高斯噪声，服从N(0,σ2)。通常噪声不可知，由SURE程序来估计阈值，定义如下：

(10)

离散阈值小波的去噪方法[10]为：

(1) 对含噪信号y(i)做小波变换得小波系数Wj,k。

(2) 对Wj,k进行阈值处理,得到估计小波系数，

2 局部投影和离散小波阈值降噪算法

局部投影降噪是将时间序列重构的相空间划分成不同的子空间，特征值较大那部分是背景信号(相对较强的信号)，特征值较小的部分是微弱信号和噪声，故可以多次采用局部投影降噪结合小波去噪算法来提取完整信号，先取较大的m，用局部投影降噪去除随机噪声部分，保留了背景信号和微弱信号，然后取较小的m，再次用局部投影降噪算法，这时特征值较大的信号为背景信号，特征值较小的即为微弱信号，但是微弱信号中还包含一些噪声，而离散小波阈值去噪可以有效的保留信号信息，所以再用离散小波阈值降噪将此噪声去除，即可得到较完整的心音信号。算法步骤如下：

(11)

再计算Cn，Cn=(R·An)T(R·An),R为m阶对角权重矩阵，用来抑制相点两端元素畸变，保留稳定的中间元素，其中R11=Rmm=1 000。

(3) 计算矩阵Cn的特征值Λ和特征向量P，其中Λ=diag(λ1,λ2,…,λm)，且(λ1≥λ2≥…≥λm)，P=[P1,P2,…,Pm]，对于特征向量P，[P1,P2,…,Pm0]且(m0

(7) 修正去噪后的相点，并更新所有相点，对更新后的相点进行加权反重构处理。

算法分离噪声的流程图如图1。

图1 分离噪声的流程图

3 实验及结果分析

3.1 数值仿真分析

为验证本文算法的有效性，以Lorenz系统产生的混沌吸引子为例进行分析。

Lorenz系统的混沌动力学方程为：

x′=σ(y-x)y′=x(R-z)-yz′=xy-bz(12)

其中:σ=10,R=28,b=8/3。设定积分步长h=0.003，取x方向的12 000个点作为研究对象，添加SNR为4.142 6高斯白噪声，重构嵌入维数m=11，时延τ=20。这时特征值较大的表征吸引子整体流行的信号，特征值较小的表征精细结构的信号，二者结合白噪声构成了该Lorenz信号。分别用离散小波阈值(DWT)、局部投影(LP)以及本文算法对其进行去噪。DWT法选取了去噪效果较好的db5小波，并采用了SURE阈值估计。

图2 Lorenz系统降噪效果图

图2中的(a)、(b)、(c)、(d)、(e)分别为原始信号、加噪信号、DWT去噪信号、LP去噪信号和本文方法去噪信号。可以看出，直接用DWT去噪效果不是很理想，LP法已基本保留了特征值较大和较小的信号成分，已将大的随机噪声去除，但是去噪后的曲线在细微处依然不是很平滑，说明特征值较小的微弱信号中还含有部分随机噪声，影响了Lorenz的精细结构；本文方法又用DWT法将微弱信号中随机噪声去除，去噪后的吸引子曲线较为平滑，有着较为精细的结构，与原始信号吸引子大致相同，对噪声有明显的抑制作用。

为了检验该算法的去噪效果，选择了两个评价指标——信噪比(SNR)和均方误差(MSE)。去噪后SNR越高，且MSE越小，表明越接近原始信号，去噪效果就越好。在不同信噪比的噪声下，分别计算三种方法降噪后的信噪比和均方误差，如表1所示。

表1 不同方法降噪后的信噪比和均方误差

由表1可知，在SNR不同的噪声中，本文方法较DWT法和LP算法都有较高的SNR和较低的MSE，说明该方法能有效的滤除噪声，且能较大程度的逼近原始信号。

最大Lyapunov指数是诊断和描述非线性混沌系统的重要参数[11]，最能表征信号的非线性混沌特征。本文采用了改进的Rosenstein方法[12-13]计算最大Lyapunov指数，它对重构时延，维数和噪声水平具有一定的鲁棒性，算出的最大Lyapunov指数较为准确。提取加噪前后以及去噪后该系统的最大Lyapunov指数如表2所示。

表2 不同方法降噪后的最大Lyapunov指数

由表2可知，加噪后的最大Lyapunov指数明显偏离实际值，导致有效信息丢失，说明噪声对最大Lyapunov指数这一特性参数有影响。DWT和LP算法虽然能使这一参数朝着原始信号逼近，但是存在较大的误差，不能作为原始信号特征值的近似，本文方法使该参数在较小的误差范围内等同于未加噪时的原始信号，证实该去噪方法能有效保留信号的非线性特征。

3.2 实测信号分析

心音信号由运动心力检测仪(ECCM，专利号01256971.2)采集，采样频率为11 025 Hz。正常心音(Normal Heart Sound，NHS)、二尖瓣狭窄(Mitral Stenosis，MS)、主动脉关闭不全(Aortic Incompetence，AI)三类心音的采样的时间序列长度分别为7 300、6 000、4 800。LP法中由互信息法和CAO算法求得三类心音的延迟时间τ和嵌入维数m见表3。

表3 不同类型心音信号的延迟时间τ和嵌入维数m

限于篇幅，本文仅列出正常心音信号的相空间降噪图。对一例含噪的正常心音进行相空间重构，分别用DWT、LP和本文方法进行降噪的效果如图3。

图4为不同降噪算法的去噪时序图，图中的ⓐ、ⓑ、ⓒ、ⓓ分别代表原始信号、DWT去噪、LP去噪和本文方法去噪后的信号。

图3 各种降噪算法的相空间去噪效果图

图4 三种类型心音的降噪时序图

图3～4中可以看出，DWT法去噪后相空间曲线流形依然杂乱，LP法降噪效果优于DWT法，已经趋于平滑；从时序图来看，在病理信号上两者都保留了体现信号病理特征的高频部分，保证了心音信号的完整性，但是LP法去噪后的信号还有些许毛刺；本文方法降噪后的心音信号相空间图曲线较为平滑，也比较清晰，去噪后的时序图几乎没有毛刺，也保留了心音信号的精细结构，说明该方法既能抑制噪声又能得到较为精确的信号，从而表明该方法的有效性。

提取去噪后信号的最大Lyapunov指数如表4所示(表中不含噪信号的最大Lyapunov指数是相关学者研究的数据[14])。可以看出，三种类型心音的最大Lyapunov指数均为正值，表明心音信号中存在混沌现象[15]，其中正常心音的最大Lyapunov指数最小，说明病理心音较正常心音混沌性强，正常心音具有较强的周期性。并且DWT法和LP法得出的最大Lyapunov指数已经明显偏离了不含噪信号的最大Lyapunov指数值，不能作为不含噪信号非线性特征的近似，而本文算法得出的最大Lyapunov指数在较小的误差范围内等同于不含噪信号，表明该方法有效保留了信号的非线性特征。

表4 不同方法降噪后的最大Lyapunov指数

4 结 论

综合局部投影和离散小波去噪的优点，本文提出了噪声自适应估计的局部投影与离散小波阈值相结合的降噪方法，并提取信号的最大Lyapunov指数这一表征信号非线性的特征量，通过Lorenz序列以及实测心音这种微振动生理信号的实验分析，证实了该方法能有效滤除噪声且能在误差范围内逼近不含噪信号，得到精确的重构信号，保留了更多的非线性特征信息，对比局部投影和离散小波阈值去噪算法，表明了该方法的优越性，为心音信号的后续处理奠定了基础。

[1] 韩阳.基于心音信号分析的心脏病诊断方法研究[D].兰州：兰州理工大学,2012.

[2] 熊新兵,陈亚光.非线性小波变换及其应用[J].中南民族大学学报(自然科学版),2005,24(2):43-45.

XIONG Xin-bing, CHEN Ya-guang. Nonlinear wavelet transform and its application[J]. Journal of south-central university for nationalities(Natural Science Edition),2005,24(2): 43-45.

[3] Emeterio J L S, RodriguezHernandez M A. Wavelet denoising of ultrasonic a-scans for detection of weak signals[C]. //IEEE on 19th International Conference on Systems, Signals and Image Processing, 2012:48-51.

[4] 张来斌,李峰,段礼祥.基于噪声水平自适应估计的往复压缩机振动信号的局部投影降噪方法[J].振动与冲击,2010,29(1):53-57.

ZHANG Lai-bin, LI Feng, DUAN Li-xiang. Local projection denoise method for reciprocating compressor vibration signals based on adaptive estimation of noise level[J].Journal of Vibration and Shock, 2010,29(1):53-57.

[5] 阳建宏,徐金梧,杨德斌,等.邻域自适应选取的局部投影非线性降噪方法[J].振动与冲击,2006,25(4):64-67.

YANG Jian-hong,XU Jin-wu, YANG De-bin et al. Nonlinear denoise method by local projection with Adaptive neighborhood selection [J].Journal of Vibration and Shock, 2006,25(4):64-67.

[6] Miller N B. Reflections from gradual transition sound absorbers[J]. Journal of the Acoustical Society of America，1958,30(10).

[7] Cao L. Practical method for determining the minimum embedding dimension of a scalar time series. Physica D: Nonlinear Phenomena[J]. 1997,110(1-2): 43-50.

[8] Fraser A M, Swinney H L. Independent coordinates for strange attractors from mutual information[J]. Physical Review A. 1986,33(2): 1134-1140.

[9] 杨志安,王光瑞,陈式刚. 用等间距分格子法计算互信息函数确定延迟时间[J]. 计算物理,1995,12(4):442-447.

YANG Zhi-an, WANG Guang-rui, CHEN Shi-gang. Spacing divided grid method to calculate the mutual information function to determine the delay time[J].Chinese Journal of Computational Physics, 1995, 12(4):442-447.

[10] 刘希灵,李夕兵,洪亮,等.基于离散小波变换的岩石SHPB测试[J].爆炸与冲击，2009,29(1):67-71.

LIU Xi-ling, LI Xi-bing, HONG Liang, et al. SHPB test of rock based on discrete wavelet transform[J]. Journal of Explosion and Shock Waves, 2009, 29(1): 67-71.

[11] 吕金虎,陆君安,陈士华. 混沌时间序列分析及其应用[M]. 武汉:武汉大学出版社,2005.

[12] 梁勇,孟桥,陆佶人. Lyapunov指数的算法改进与加权预测[J].声学技术,2006,25(5):463-467.

LIANG Yong, MENG Qiao, LU Ji-ren. Improved algorithm of Lyapunov exponent and weighted prediction[J]. Technical Acoustics, 2006,25(5):463-467.

[13] Kantz H, Schreiber T. Nonlinear time series analysis[M]. Cambridge: Cambridge Univ. Press, 1977：65-74.

[14] 丁晓蓉.基于混沌理论的心音信号非线性动力学分析[D].重庆:重庆大学,2012.

[15] 姜桂仁. 混沌时序的特征量分析及相空间重构[D]. 镇江:江苏大学,2005.