CN104367315B - 一种提取心电图中j波信号的盲源分离方法 - Google Patents

Abstract

本发明涉及J波的提取方法，具体为以吉文斯旋转矩阵作为分离矩阵，基于循环平稳度准则的一种提取心电图中J波信号的盲源分离方法。J波信号是心电图（Electrocardiogram ECG)上心室复极的一种新指标，临床数据显示ECG中的J波信号形态各异且随机性很强，有时会淹没在正常ECG中观测不到，为准确提取到J波信号，本文预先采用小波包变换(Wavelet Packet Transform WPT）分解重建定位出J波所在的ST段，以减少盲源分离(Blind Source Separation BSS）的计算量；再根据ECG和J波信号的特点，使用基于高阶循环累积量的BSS方法有效的分离循环平稳的ECG和非周期的J波信号，为下一步在时频域研究J波信号，以区分其良性与高危性形态提供了保障。

Description

一种提取心电图中J波信号的盲源分离方法

技术领域

本发明涉及J波的提取方法，具体为以吉文斯旋转矩阵作为分离矩阵，基于循环平稳度准则的一种提取心电图中J波信号的盲源分离方法。

背景技术

J波属于一种正常的ECG(心电信号)变异，但当J波形态、时限及幅值发生显著变化时,如：跟随在QRS波之后的J点从基线移位持续时间＞0.03s，振幅＞0.05mv，形成尖峰状、驼峰状、圆顶状的波形或者ST段抬高，可能预示着室速、室颤等恶性心律失常事件，甚至导致心源性猝死的发生。近年来，变异的J波已经作为ECG心室复极的一种新指标得到医学界的重视，因此，区分J波的良性与高危性形态具有非常重要的临床意义。目前对J波的分类是医生凭自身经验通过观察ECG进行判断，这样仅局限于在时域研究J波信号的幅值、波形以及发生位置，尤其是在J波振幅较小，不易肉眼观察的情况下，很容易造成判断误差。因此提取J波非常重要，当J波经过分析，从中提取J波的特性，医生再结合上临床病例，可更准确诊断室速、室颤、致命性的恶性心律失常等心室疾病，但是现在没有一种可有效提取J波的方法。

发明内容

本发明为了解决现在没有一种方法能有效提取J波的问题，提供了一种提取心电图中J波信号的盲源分离方法。

本发明是采用如下的技术方案实现的：一种提取心电图中J波信号的盲源分离方法，包括以下步骤：

1)获取J波综合征患者的十二导联ECG，作为初始观测信号；

2)使用小波包变换定位初始观测信号中S点和T点，截取含有J波的ST段信号作为新观测信号：选用db3作为小波函数ψ(t)，其对应的尺度函数为φ(t)，由小波函数ψ(t)和尺度函数φ(t)计算得到滤波器H₀(z)和H₁(z)的滤波器系数h₀(k)和h₁(k)，将初始观测信号通过滤波器H₀(z)和H₁(z)在2^j个子空间中进行分解，其分解形式为二进制树结构，分解后的每个节点代表初始观测信号的一个特定频带，根据ECG中T点和S点的频谱范围，对相应的节点用下层节点进行重构，在重构信号中搜索极大值，对应于初始观测信号的T点和S点，从初始观测信号中截取J波所在的ST段作为新观测信号；

3)使用吉文斯旋转矩阵对新观测信号进行盲源分离:

建立BSS的基础模型：新观测信号 $X=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ M\\ x_m\left(t\right)\end{array}\right],$ 含有J波信号的源信号 $S=\left[\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ M\\ s_n\left(t\right)\end{array}\right],$ 观测信号X＝AS，其中A为m*n维列满秩的混合矩阵；

白化预处理：新观测信号X的协方差矩阵进行特征值分解C_X＝E[XX^H]＝DΛD^H，其中Λ为C_X的特征值矩阵，D为对应的特征向量矩阵，则白化矩阵白化处理信号Z＝VX，表示对观测信号X进行了白化处理,白化处理后的信号Z的协方差矩阵为设U＝VA，则矩阵U为酉矩阵，由Z＝VX＝VAS＝US可知，白化处理信号Z是源信号S与酉矩阵U的混合；

对白化处理信号进行吉文斯变换，得到分离信号：

使用吉文斯旋转矩阵 $W_{\mathrm{ij}}=\left[\begin{array}{ccccccc}1& L& 0& L& 0& L& 0\\ M& O& M& & M& & M\\ 0& L& \cos \mathrm{\θ}& L& \sin \mathrm{\θ}& L& 0\\ M& & M& O& M& & M\\ 0& L& -\sin \mathrm{\θ}& L& \cos \mathrm{\θ}& L& 0\\ M& & M& & M& O& M\\ 0& L& 0& L& 0& L& 1\end{array}\right]$ 作为分离矩阵变换白化处理信号Z，得到分离信号Y＝W_ij ^HZ＝W_ij ^H(US)＝W_ij ^HUS＝[y₁(t)；L y_i(t)；L y_j(t)；L y_n(t)]，所得的分离信号Y即为源信号X的估计；

4)选用循环平稳度DCS作为分离准则选取参数θ，从而得到分离矩阵W：若使任意分离信号y_i(t)的循环平稳度值g(θ)＝DCSα在心电信号的循环频率α处取得最大值，则存在θ_k+1＝θ_k+△θ使得W_ij(θ_k+1)＝U，此时Y＝I_nS，即Y为S的最优估计；

5)S＝U^HZ，由源信号S可得J波的波形图，即提取到了J波。

Daubechies小波简称db小波，dbN表示小波的阶次是N次，是本领域的技术人员公知的函数。

本发明在盲源分离算法之前，使用小波包分解，利用分解系数重构T波和S波的波峰能量所在频段的信号，以避免R点在奇异值检测中的干扰，然后选择合适的阈值，在重构出来的信号上搜索极大值对应为初始观测信号上的T点和S点，可以相对准确的定位到心电图上T点和S点，截取到J波所在的ST段，减少了盲源分离的计算量，并且提高了分离准确度。在盲源分离的过程中先进行白化预处理来简化过程，然后使用吉文斯旋转矩阵作为分离矩阵，以综合考虑十二导联的观测信号，分离的结果选取最优作为源信号的估计，提高了分离精度，在算法选择参数的时候，充分利用心电信号的循环平稳特性，以循环平稳度作为分离准则，并且使用变步长的学习速率，最终不仅可以得到源信号的最优估计，有效的提取到J波信号，而且提高了算法的收敛速度。

附图说明

图1为本发明的流程图。

具体实施方式

一种提取心电图中J波信号的盲源分离方法，包括以下步骤：

1)获取J波综合征患者的十二导联ECG，作为初始观测信号；

2)从初始观测信号中提取含有J波的新观测信号：选用db3作为小波函数ψ(t)，其对应的尺度函数为φ(t)，由小波函数ψ(t)和尺度函数φ(t)计算得到滤波器H₀(z)和H₁(z)的滤波器系数h₀(k)和h₁(k)，

$\left\{\begin{array}{c}{h}_0\left(k\right)=<{\mathrm{\&phi;}}_{\mathrm{1,0}}\left(t\right),{\mathrm{\&phi;}}_{0,k}\left(t\right)>=\&Integral;\mathrm{\&phi;}\left(\frac{t}{2}\right){\mathrm{\&phi;}}^{*}(t-k)\mathrm{dt}\\ h_1\left(k\right)=<{\mathrm{\&psi;}}_{\mathrm{1,0}}\left(t\right),{\mathrm{\&phi;}}_{0,k}\left(t\right)>=\&Integral;\mathrm{\&psi;}\left(\frac{t}{2}\right){\mathrm{\&phi;}}^{*}(t-k)\mathrm{dt}\end{array}\right.---\left(1\right),$ 将初始观测信号用滤波器H₀(z)和H₁(z)在2^j个子空间中进行分解，其分解形式为二进制树结构，结构中各个子空间的基函数和的递推公式为式： $\left\{\begin{array}{c}{\mathrm{\&psi;}}^{2p}\left(\frac{t}{2^{j+1}}\right)=\sqrt{2}\underset{k=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{h}_0\left(k\right){\mathrm{\&psi;}}^p(\frac{t}{2^j}-k)\\ {\mathrm{\&psi;}}^{2p+1}\left(\frac{t}{2^{j+1}}\right)=\sqrt{2}\underset{k=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{h}_1\left(k\right){\mathrm{\&psi;}}^p(\frac{t}{2^j}-k)\end{array}\right.---\left(2\right),$ 结构中每个子空间结点(j，p)代表初始观测信号的一个特定频带，其分解系数可表示为： $d_j^p\left(n\right)=<x\left(t\right),{\mathrm{\&psi;}}_j^p\left(t\right)>=\&Integral;x\left(t\right){\mathrm{\&psi;}}^p(\frac{t}{2^j}-k)\mathrm{dt}---\left(3\right),$ 根据ECG中T点和S点的频谱范围，对相应的节点用下层节点的分解系数进行重构，在重构信号中搜索极大值，对应于初始观测信号的T点和S点，从初始观测信号中截取J波所在的ST段作为新观测信号；

3)使用吉文斯旋转矩阵对新观测信号进行盲源分离：

建立BSS的基础模型：新观测信号 $X=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ M\\ x_m\left(t\right)\end{array}\right],$ 含有J波信号的源信号 $S=\left[\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ M\\ s_n\left(t\right)\end{array}\right],$ 观测信号X＝AS，其中A为m*n维列满秩的混合矩阵；

白化预处理：新观测信号X的协方差矩阵进行特征值分解C_X＝E[XX^H]＝DΛD^H，其中Λ为C_X的特征值矩阵，D为对应的特征向量矩阵，则白化矩阵白化处理信号Z＝VX，表示对观测信号X进行了白化处理，白化处理信号Z的协方差矩阵设U＝VA，则矩阵U为酉矩阵，由Z＝VX＝VAS＝US可知，白化处理信号Z是源信号S与酉矩阵U的混合；

对白化处理信号进行吉文斯变换，得到分离信号：

使用吉文斯旋转矩阵 $W_{\mathrm{ij}}=\left[\begin{array}{ccccccc}1& L& 0& L& 0& L& 0\\ M& O& M& & M& & M\\ 0& L& \cos \mathrm{\&theta;}& L& \sin \mathrm{\&theta;}& L& 0\\ M& & M& O& M& & M\\ 0& L& -\sin \mathrm{\&theta;}& L& \cos \mathrm{\&theta;}& L& 0\\ M& & M& & M& O& M\\ 0& L& 0& L& 0& L& 1\end{array}\right]$ 作为分离矩阵变换白化处理信号Z，得到分离信号Y＝W_ij ^HZ＝W_ij ^H(US)＝W_ij ^HUS＝[y₁(t)；L y_i(t)；L y_j(t)；L y_n(t)]，其中所得的分离信号即为源信号的估计；

4)选用循环平稳度作为分离准则选取参数θ，从而得到分离矩阵W：分离信号的四阶循环累积量的循环平稳度定义为式：

${\mathrm{DCS}}^{\mathrm{\&alpha;}}=\frac{\&Integral;\&Integral;\&Integral;{\left|C_{4y}^{\mathrm{\&alpha;}}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3)\right|}^2{\mathrm{d\&tau;}}_1{\mathrm{d\&tau;}}_2{\mathrm{d\&tau;}}_3}{\&Integral;\&Integral;\&Integral;{\left|C_{4y}^0({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3)\right|}^2{\mathrm{d\&tau;}}_1{\mathrm{d\&tau;}}_2{\mathrm{d\&tau;}}_3}---\left(4\right)$ 其中，τ₁，τ₂，τ₃表示对于固定的时间延迟，对任意分离信号若s_q(t)(1≤q≤n)是非周期的J波信号，则其循环累积量应为零，即又因为各个源信号相互独立，所以互累积量均为零，综上可得， $C_{{4y}_i}^{\mathrm{\&alpha;}}=\underset{p=1;p\&NotEqual;q}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{ip}}{\left(\mathrm{\&theta;}\right)}^4{C}_{{4s}_p}^{\mathrm{\&alpha;}}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3);C_{{4y}_i}^0=\underset{p=1;}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{ip}}{\left(\mathrm{\&theta;}\right)}^4{C}_{{4s}_p}^0({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3),$ 代入(4)式，得到一个关于θ的函数 $g\left(\mathrm{\&theta;}\right)={\mathrm{DCS}}^{\mathrm{\&alpha;}}=\frac{\&Integral;\&Integral;\&Integral;{\left|\underset{p=1;p\&NotEqual;q}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{ip}}{\left(\mathrm{\&theta;}\right)}^4{C}_{{4s}_p}^{\mathrm{\&alpha;}}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3)\right|}^2{\mathrm{d\&tau;}}_1{\mathrm{d\&tau;}}_2{\mathrm{d\&tau;}}_3}{\&Integral;\&Integral;\&Integral;{|\underset{p=1;p\&NotEqual;q}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{ip}}{\left(\mathrm{\&theta;}\right)}^4{C}_{4\mathrm{sp}}^0({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3)+f_{\mathrm{iq}}{\left(\mathrm{\&theta;}\right)}^4{C}_{{4s}_q}^0({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3)|}^2{\mathrm{d\&tau;}}_1{\mathrm{d\&tau;}}_2{\mathrm{d\&tau;}}_3}---\left(5\right)$

分离信号y_i(t)为ECG信号，其循环平稳度值g(θ)在心电信号的循环频率α处取得最大值，设定初试参数θ₀＝0，迭代步长 $\mathrm{\&Delta;\&theta;}=\left\{\begin{array}{cc}\frac{1}{\mathrm{\&pi;}}& R\&GreaterEqual;0.3,k\&le;k_0\\ \frac{1}{\mathrm{\&pi;}}{e}^{-k_0(k-k_0)}& R<0.3,k>k_0\end{array}\right.,$ 其中 $R_{\mathrm{ij}}=\frac{\left|\mathrm{\&Sigma;}{y}_i\left(t\right)y_j\left(t\right)\right|}{\sqrt{\mathrm{\&Sigma;}{y_i}^2\left(t\right){y_j}^2\left(t\right)}},$ 通过迭代θ_k+1＝θ_k+△θ使得g(θ_k+1)取得最大值，此时W_ij(θ_k+1)＝U，Y＝W^HUS＝I_nS，即Y为S的最优估计，

5)S＝U^HZ，MATLAB仿真得到源信号S中各个分量的波形图，即可得到J波的波形图。

Claims (1)

1.一种提取心电图中J波信号的盲源分离方法，其特征在于包括以下步骤：

1)获取J波综合征患者的十二导联ECG，作为初始观测信号；

2)使用小波包变换定位初始观测信号中S点和T点，截取含有J波的ST段信号作为新观测信号：选用db3作为小波函数ψ(t)，其对应的尺度函数为φ(t)，由小波函数ψ(t)和尺度函数φ(t)计算得到滤波器H₀(z)和H₁(z)的滤波器系数h₀(k)和h₁(k)，将初始观测信号通过滤波器H₀(z)和H₁(z)在2^j个子空间中进行分解，其分解形式为二进制树结构，分解后的每个节点代表初始观测信号的一个特定频带，根据ECG中T点和S点的频谱范围，对相应的节点用下层节点进行重构，在重构信号中搜索极大值，对应于初始观测信号的T点和S点，从初始观测信号中截取J波所在的ST段作为新观测信号；

3)使用吉文斯旋转矩阵对新观测信号进行盲源分离:

建立BSS的基础模型：新观测信号 $X=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ .\\ .\\ .\\ x_m\left(t\right)\end{array}\right],$ 含有J波信号的源信号 $S=\left[\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ .\\ .\\ .\\ s_n\left(t\right)\end{array}\right],$ 观测信号X＝AS，其中A为m*n维列满秩的混合矩阵；

白化预处理：新观测信号X的协方差矩阵进行特征值分解C_X＝E[XX^H]＝DΛD^H，其中Λ为C_X的特征值矩阵，D为对应的特征向量矩阵，则白化矩阵白化处理信号Z＝VX，表示对观测信号X进行了白化处理,白化处理后的信号Z的协方差矩阵为C_Z＝E[ZZ^H]＝E[VX(VX)^H]＝VE[XX^H]V^H＝VAA^HV^H＝VA(VA)^H＝I，设U＝VA，则矩阵U为酉矩阵，由Z＝VX＝VAS＝US可知，白化处理信号Z是源信号S与酉矩阵U的混合；

对白化处理信号进行吉文斯变换，得到分离信号：

使用吉文斯旋转矩阵作为分离矩阵变换白化处理信号Z，得到分离信号 $Y={W_{ij}}^{H}Z={W_{ij}}^H\left(US\right)={W_{ij}}^{H}US=\&lsqb;y_1\left(t\right);...y_i\left(t\right);...y_j\left(t\right);...y_n\left(t\right)\&rsqb;,$ 所得的分离信号Y即为源信号X的估计；

4)选用循环平稳度DCS作为分离准则选取参数θ，从而得到分离矩阵W：若使任意分离信号y_i(t)的循环平稳度值g(θ)＝DCS^α在心电信号的循环频率α处取得最大值，则存在θ_k+1＝θ_k+△θ使得W_ij(θ_k+1)＝U，此时Y＝I_nS，即Y为S的最优估计；

5)S＝U^HZ，由源信号S可得J波的波形图，即提取到了J波。