CN102961129B - 一种远程医疗的异常心电张量分析方法 - Google Patents

Abstract

本发明公布了一种远程医疗的异常心电张量分析方法，首先通过远程方式采集大量的标准12导联心电数据，然后通过短时傅立叶变换将心电转换为高维度的张量心电数据。然后以高维张量心电数据直接作为特征，使用直接以张量数据直接作为输入的特征抽取和特征降维的算法提取出直接用来分类的心电特征。由于这种方法是基于TTV变换法则的，所以最终可以得到基于向量存储的特征，然后使用SVM分类方法对这些向量特征进行分类。这种方法以张量心电数据直接作为输入，充分利用了心电的多导联心电的结构信息，消除了原先单导联心电单独分析带来的不精准缺陷，实现了心电分析的有效性。

Description

一种远程医疗的异常心电张量分析方法

技术领域

本发明涉及的是一种针对远程心电诊断平台，多导联张量心电数据的特征抽取和辅助分类方法。

背景技术

本系统基于物联网技术，利用远程心电诊断平台，将上海市三级甲等医院的优质医疗服务延伸至远程的县市级甚至社区级医疗中心。设计了个性化服务系统，提高了远程医疗与健康监护平台的可用性及使用效率。利用基于物联网的构建远程医疗云服务平台，并建立基层医疗单位的远程心电诊断服务示范。构建远程医疗云服务主要包括医疗信息采集、传输、数据处理、反馈诊断等关键技术。当前心电信号采集设备和医疗数据远程传输技术已经成熟，本项目核心技术主要集中在医疗数据挖掘与诊断辅助决策支持平台方面，使得能为医疗中心医生实现快速准确的诊断支持。该医疗服务平台逐步实现基层医疗单位、家庭的两种远程医疗与健康监护服务模式，形成物联网下的新型远程医疗与健康监护医疗服务模式。该项目研发医疗技术云服务平台主要创新之处包括三个方面：远程医疗诊断服务模式创新、医疗数据挖掘技术平台创新和医疗数据服务创新。

基于远程心电云服务平台，通过实现心电分析的辅助决策支持，极大提高了诊断中心医生的诊断速度和精度，使得远程心电诊断新型服务模式成为可能。在远程心电智能诊断技术方面，主要在以下几个方面提供辅助决策支持：通过对心电数据分析，提供准确的心电分析基本参数；采用结合机器学习方法和智能模式识别方法，通过对心电信号的有效表征，利用最新模式识别方法（如核方法）对心电特征进行模式分类，实现对心电数据的预诊，云服务平台能即时推荐给熟悉相应病症的医生进行诊断，提高了诊断的效率。本系统通过网格化的数据存储很好的为非结构化的海量数据提供高效存储管理和查询服务功能。这样一套完备的分析数据库，为以后的医学和工程研究提供一个大型的完备数据库

在科学计算领域，基于心电数据进行特征抽取和降维并且进行分类分析的方法比较有名的就是C.Saritha,V.Sukanya,andY.NarasimhaMurthy,“ECGSignalAnalysisUsingWaveletTransforms,”BulgarianJournalofPhysics,vol.35,pp.68-77,2008.（作者：C.Saritha,V.Sukanya,Y.NarasimhaMurthy,题目：基于小波变换的心电信号分析，杂志：保加利亚物理杂志，2008年35卷，68-77页）和Kuo-KuangJen,andYean-RenHwang,“ECGFeatureExtractionandClassificationUsingspectrumandNeuralNetworks,”JournalofMedicalandBiologicalEngineering,vol.28,no.1,2008.（作者：Kuo-KuangJen,Yean-RenHwang,题目：基于频谱和神经网络的心电特征抽取和分类，杂志：医学和生物工程杂志，2008年28卷，第一期）对于第一篇文章，首先对信号进行不同小波基底的小波分解，然后通过检测其系数来提取有效的针对不同疾病波形的特征进行分析。第二篇文章是将原来的心电转换成频谱，然后通过提取心电频谱上的特征来抽取有效的针对不同疾病的特征。最后使用神经网络来对原始的心电进行分类。需要指出的是，所有这些方法都是针对单导联心电的，而实际使用的心电诊断数据往往是多导联的。一般目前在医学上临床使用的心电是12标准导联的。这样的话针对单一导联的算法就不能完全满足需要，即使通过将多个导联心电拼接成单个导联的方法，多个导联的结构信息被破坏，有用信息有所减少并不是非常有效的方法。因此，需要研发新的将多导联心电张量心电数据作为直接输入的特征抽取和分析方法，来处理上述论文中方法无法直接处理的问题。

发明内容

本发明的目的在于克服现有分散式的心电诊断系统的非集中化带来的效率低下的不足，同时又针对原先基于单导联心电分析不够准确的缺陷，提出了一种针对多导联张量心电数据的心电分析方法。这种方法以张量心电数据直接作为输入，充分利用了心电的多导联心电的结构信息，消除了原先单导联心电单独分析带来的不精准缺陷，实现了心电分析的有效性。

本发明是通过以下技术方案实现的，其具体步骤如下：

一种远程医疗的异常心电张量分析方法，包括下述步骤：

（1）构造张量数据：

a．心电数据采集：

采集标准12导联心电数据；

b．数据预处理和去燥：

对采集的12导联心电数据首先对信号通过50hz的陷波滤波器进行滤波处理，然后对数据进行DB6小波的小波变换分解，然后去除其中最高频的和最低频的信号成份；

c．波形检测：

再次对小波进行DB6小波分解，然后找寻其中level2小波系数，然后采用过零点检测方法检测心电的QRS波尖峰，然后依次去除R波后检测P波和T波尖峰，然后采用detrend算法计算出波形的基线，通过计算各个波形与基线的交点，确定P波QRS波T波的起始和结束，也就是onset和offset；

d．逐跳切割：

一次采集心电大约包含20秒的心电数据，也就是25跳左右的心电，对于心电一次beat也就是一个P波一个QRS波一个T波进行切割；

e.R波对齐

对每一个beat的心电针对R波的峰值进行对齐，并且切割成统一的长度；

f.短时fourier变换

为了有效抽取时频域的心电诊断特征，通过使用窗口为128点的短时fourier变换对心电进行转换，最后心电被转换为128×128×12的时频空的3阶张量；这里的空指的是导联位置就是指导联轴；对于12导联（lead×time）ECG信号，s[l,n]表示离散信号在时间点n对于导联l；在时间点nΔt和频率f的短时fourier变换如下式：

$\begin{array}{c}\mathrm{STFT}\left\{s\right[l,n\left]\right\}(m,w)=S(l,m,n)\\ ={\mathrm{\Σ}}_{m=-\∞}^{\∞}w(n-m)s(l,m)e^{-j2\mathrm{\πfm}}\end{array}$

这里w[n]是一个窗函数，变换完成后数据变成一个3阶张量；

（2）基于TTV映射的张量特征抽取：

a.根据原始张量数据的判别性计算有效的投影张量：

$\begin{array}{c}{u}_k^l{|}_{l=1}^M={\arg }_{u_k^l{|}_{l=1}^M}\max {\left(\frac{1}{n}{\mathrm{\Σ}}_{i=1}^c\left((M_i^k-M_k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\mathrm{}}_{\×l}{\left(u_k^l\right)}^T\right)((M_i^k-M_k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\mathrm{}}_{\×l}{\left(u_k^l\right)}^T\right)}^T\\ -{\mathrm{\ζ}}_k^l{\mathrm{\Σ}}_{i=1}^c{\mathrm{\Σ}}_{j=1}^{n_i}\left((X_{\mathrm{ji}}^k-M_i^k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\mathrm{}}_{\×l}{\left(u_k^l\right)}^T\right)\×{\left((X_{\mathrm{ji}}^k-M_i^k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\mathrm{}}_{\×l}{\left(u_k^l\right)}^T\right)}^T)\end{array}$

b.根据优化张量可分性特点计算更为优越的投影张量：

$\underset{i=1}{\overset{C_c^2}{\mathrm{\Σ}}}\underset{j_1=1}{\overset{n_{i1}}{\mathrm{\Σ}}}\underset{j_2=1}{\overset{n_{i2}}{\mathrm{\Σ}}}\left((X_{j_1}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\×}_l{\left(u_k^l\right)}^T\right){\×\left((X_{j_1}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\×}_l{\left(u_k^l\right)}^T\right)}^T+\left((X_{j_2}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\×}_l{\left(u_k^l\right)}^T\right)\×{\left((X_{j_2}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\×}_l{\left(u_k^l\right)}^T\right)}^T$

c.对原始数据去除已抽取投影张量的维度，调整原始数据的结构：

$X_{\mathrm{ij}}^k=X_{\mathrm{ij}}^{k-1}-{\mathrm{\λ}}^{k-1}{u}_{k-1}^1\⊗u_{k-1}^2\⊗....\⊗u_{k-1}^M$

d.对原始数据计算投影张量进行加权处理：

通过调整每个张量的权重，来规避不合理张量的影响，等式如下：

$S_{\mathrm{oo}}=\underset{i,j}{\mathrm{\Σ}}{w}_{\mathrm{ij}}\underset{x{\∈A}_i,y{\∈A}_j}{\mathrm{\Σ}}w\left(d_{\mathrm{xy}}\right)S_{\mathrm{xy}}i\≠j$

最简单的方式就是取张量间距离的倒数distance(),

w(d_xy)=d_xy ^-n

或者如下式定义：

$w\left(d_{\mathrm{xy}}\right)=\left\{\begin{array}{cc}=1& \mathrm{if}{d}_{\mathrm{xy}}\∈N_{\%}~M_{\%}\\ =0& \mathrm{if}{d}_{\mathrm{xy}}\∉N_{\%}~M_{\%}\end{array}\right.$

组合两种形式得到：

$w\left(d_{\mathrm{xy}}\right)=\left\{\begin{array}{cc}={d_{\mathrm{xy}}}^{-n}& \mathrm{if}{d}_{\mathrm{xy}}\∈N_{\%}~M_{\%}\\ =0& \mathrm{if}{d}_{\mathrm{xy}}\∉N_{\%}~M_{\%}\end{array}\right.$

e.算法整体迭代收敛过程如下：

整个算法过程，是一个严格的单调收敛过程，逐次迭代目标值关系如下式所示：

$\begin{array}{c}{a}_k=g(u_k^1,1)\≤g(u_k^2,1)\≤...\≤g(u_k^M,1)\≤g(u_k^1,2)\≤g(u_k^2,2)\\ \≤...\≤g(u_k^1,t)\≤g(u_k^2,t)\≤...\≤g(u_k^1,T)\≤g(u_k^2,T)\\ \≤...\≤g(u_k^M,T)\\ =b_k\end{array}$

T-＞无穷时，算法收敛于最终目标极限收敛值；

f.判断计算过程结束终止条件：

使用如下方法来判断算法是否收敛，并且判断算法什么时候应该停止；误差值小于一定阈值则判断算法停止：

||F_k-F_k-1||_Fro≤ε

采用这种方法来判断收敛与否和终止整个算法；

（3）选择合理初值：

求解一个最为近似的张量：

$\min f(a^{\left(1\right)},...,a^{\left(N\right)})\≡\frac{1}{2}{\left|\right|Z-\left[\right[a^{\left(1\right)},...,a^{\left(N\right)}\left]\right]\left|\right|}^2$

a．无约束情况张量算法：

交替最小二乘方法求解目标等式如下：

其展开等式如下式:

$={\min }_{a^{\left(n\right)}}{\left|\right|Z_{\left(n\right)}-a^{\left(n\right)}{(a^{\left(N\right)}\⊗...\⊗a^{(n-1)}\⊗a^{(n+1)}\⊗...\⊗a^{\left(1\right)})}^T\left|\right|}^2$

这里表示kroneckerproduct，而Z_(n)表示按第nmode转换张量Z一个矩阵；这个问题的解犹如下式：

张量梯度下降方法：

可以将上述目标函数进行展开，写成如下的形式：

第一项没有涉及变量，所以：

$\frac{\∂f_1}{\∂a^{\left(n\right)}}=0$

这里0表示一个0向量，长度为I_n，第二项犹如如下式子：

$f_2\left(x\right)=Z{\×}_{m=1}^N{a}_r^{\left(n\right)}$

$={\left(Z{\×}_{m=1,m\≠n}^N{a}_r^{\left(m\right)}\right)}^T{a}_r^{\left(n\right)}$

第二项求导之后得到如下式:

$\frac{\∂f_2}{{\∂a}^{\left(n\right)}}=\left(Z{\×}_{m=1,m\≠n}^N{a}_r^{\left(m\right)}\right)$

第三项如下：

$f_3\left(x\right)=\underset{m=1}{\overset{N}{\mathrm{\Π}}}{a^{\left(m\right)}}^T{a}^{\left(m\right)}$

因此

$\frac{\∂f_3}{\∂a^{\left(n\right)}}=2\left(\underset{m=1,m\≠n}{\overset{N}{\mathrm{\Π}}}{a^{\left(m\right)}}^T{a}^{\left(m\right)}\right)a^{\left(n\right)}$

综合以上三项就可以得到最终结果；

b.有约束张量情况:

带约束非线性最小二乘：

带约束优化问题，通过转换对求解的约束条件有所放松，再去求解就可以得到更加优越的计算结果，如下式：

$\begin{array}{cc}\min \frac{1}{2}{\mathrm{\Σ}}_{i=1}^l{f_i\left(x\right)}^2& \min \frac{1}{2}{z}^{T}z\\ g_j\left(x\right)=0j=1,...,m_e& f_i\left(x\right)-z_i=0i=1,...,l\\ g_j\left(x\right)\≥0j=m_e+1,...,m& g_j\left(x\right)=0j=1,...,m_e\\ x_l\≤x\≤x_u& g_j\left(x\right)\≥0j=m_e+1,...,m\\ & x_l\≤x{\≤x}_u\end{array}$

将左边的优化问题转换为右边的优化问题，然后将结果代入，使用标准SQP算法求解，来求得最终的优化结果；

（4）分类比对：

最后使用SVM来对抽取出来的以vector方式存储的向量特征进行分类，

这直接通过求解如下优化主问题：

$\underset{W,b,\mathrm{\ξ}}{\min }C{\mathrm{\Σ}}_{n=1}^N{\mathrm{\ξ}}_n+\frac{1}{2}{\left|\right|w\left|\right|}^2$

Subjecttoy_i(w^Tφ(x_i)+b)≥1-ξ_n,ξ_n≥0,i=1,2…,n

这里的参数C>0在松弛变量和惩罚因子间的平衡，而他的Lagrangian乘子变换成如下等式：

$L(w,b,a)=\frac{1}{2}{\left|\right|w\left|\right|}^2+C{\mathrm{\Σ}}_{n=1}^N{\mathrm{\ξ}}_n-{\mathrm{\Σ}}_{n=1}^N{a}_n\{t_{n}y\left(x_n\right)-1+{\mathrm{\ξ}}_n\}-{\mathrm{\Σ}}_{n=1}^N{\mathrm{\μ}}_n{\mathrm{\ξ}}_n$

这里{a_n≥0}而且{μ_n≥0}是Lagrangian乘子，而对偶Lagrangian问题如下式：

$\tilde{L}\left(a\right)={\mathrm{\Σ}}_{i=1}^N{a}_n-\frac{1}{2}{\mathrm{\Σ}}_{n=1}^N{\mathrm{\Σ}}_{n=1}^N{a}_n{a}_m{t}_n{t}_{m}k(x_n,x_m)$

它有约束项0≤a_n≤C和.而且k(x,x′)=φ(x)^Tφ(x′)是核函数。

本发明的有益效果是：

(1)充分利用远程心电集中式的优点，充分利用云技术和物联网技术。有效提升心电诊断的高度集中化统一调度诊断。使得处理速度非常快，计算效率高，可以做到实时处理。

(2)系统算法针对12导联医用标准数据库，比起原始但导联的算法具有更大的优势。

(3)提出直接针对张量心电的特征抽取和特征降维，使得高阶张量心电可以保留完整的结构信息不丢失，使诊断分类可以更加有效的进行。

(4)采用ECG多导联信号的时频域特性作为诊断系数，比起单纯时域的特征具有更好的性能。

附图说明

图1为本发明硬件架构图。

图2为本发明数据流程图。

图3为本发明方法实现流程图。

图4为本发明GTR1DA方法进行特征抽取后降维到3维和的分类效果图。

图5为ica方法进行特征抽取后降维到3维和的分类效果图。

图6为lda方法进行特征抽取后降维到3维和的分类效果图。

图7为pca方法进行特征抽取后降维到3维和的分类效果图。

图8为TR1DA方法进行特征抽取后降维到3维和的分类效果图。

图9为UMPCA方法进行特征抽取后降维到3维和的分类效果图。

具体实施方式

我们首先通过远程方式采集大量的标准12导联心电数据，然后通过短时傅立叶变换（Short-timeFourierTransform）将心电转换为高维度的张量心电数据（通常128×128×12）。然后以高维张量心电数据直接作为特征，使用直接以张量数据直接作为输入的特征抽取和特征降维的GeneralizedTensorRankOneDiscriminantAnalysis算法提取出直接用来分类的心电特征。由于这种方法是基于TTV变换法则的，所以最终可以得到基于向量存储的特征，然后使用SVM分类方法对这些向量特征进行分类。提出一种以张量作为直接输入的特征降维和特征抽取算法来直接对于张量的心电数据进行处理是本发明的核心创新点。

为了更好地描述本发明的内容，首先描述一下远程心电系统的架构。

心电的远程医疗系统主要实现心电信号的远程采集、传输、诊断功能，框架上结合了物联网和分布式的优势特点，总共分为采集端、服务器端、诊断端三大模块，如图1所示。

患有心脏疾病的患者不需要前往大型医院排队，而是由社区医院的医生携带轻便的采集设备前往居民区，在现场对病人进行心电数据的采集；采集结束后，采集医生通过蓝牙、wifi等方式将心电数据从采集设备传输至数据接收前置机，该前置机可以是手机、平板电脑或PC等；前置机将心电数据和病人基本信息按照特定的网络传输协议进行打包，传输至服务器端的数据格式转换服务器。

数据格式转换服务器在接收到数据包之后，按照约定的传输协议进行解包，生成格式化的心电数据类，发送至数据库服务器，将心电数据保存至数据仓库，其中，由于各采集中心的配套设备存在差异，因此打包传输的网络传输协议可能有所不同，该数据格式转换服务器的作用及时形成统一接口，产生统一的数据类结构；数据库服务器对数据进行保存和判断，如果是新采集的心电数据，则立即推送至诊断中心服务器。

诊断中心服务器接收到新的推送数据（诊断请求），会根据各诊断医生的工作量动态地给各个诊断端推送诊断数据；诊断端是一款带有GUI界面的心电在线诊断软件，诊断医生根据波形显示和统计参数给出诊断结论并生成诊断报告，并将诊断报告发送至数据库服务器。数据库服务器一旦接收到新的诊断报告过来，会把该报告显示给该心电数据所属的病人患者及其采集医生。

整个数据流程图如图2所示。

为了更好地说明本系统的核心算法，给出如下定义：

如图3所示本发明包括以下步骤：

1.构造张量数据：

我们的系统使用的标准12导联心电数据，而不是单一导联的心电数据。而心电的有效特征又往往不是仅仅存在于时域信号中，在频域上往往也存在着有效的分类特征。所以，我们的方法首先通过预处理，去除心电波形中的各种干扰和噪声，然后通过短时Fourier变换，将心电转换为张量时频域的表征形式。

具体来说，包含如下步骤：

(1)数据预处理和去燥

对于心电来说，必须要去除包括工频噪声，肌肉电干扰，基线漂移等多类噪声。首先对信号做工频频率的滤波处理，使信号通过50hz的陷波滤波器。然后对数据进行DB6小波的小波变换分解，然后去除其中最高频的和最低频的信号成份。因为低频部分包含了基线漂移，而最高频部分包含了肌肉点干扰。

(2)波形检测

再次对小波进行DB6小波分解，然后找寻其中level2小波系数，然后采用过零点检测方法检测心电的QRS波尖峰。然后依次去除R波后检测P波和T波尖峰。然后采用detrend算法计算出波形的基线，通过计算各个波形与基线的交点，确定P波QRS波T波的起始和结束。也就是onset和offset。

(3)逐跳切割

一次采集心电大约包含20秒的心电数据，也就是25跳左右的心电。我们对于心电一次beat也就是一个P波一个QRS波一个T波进行切割。

(4)R波对齐

出去诊断有效性的考虑，我们对波形进行了对齐，也就是对每一个beat的心电针对R波的峰值进行对齐，并且切割成统一的长度。

(5)短时fourier变换

为了有效抽取时频域的心电诊断特征，我们通过使用窗口为128点的短时fourier变换对心电进行转换，最后心电被转换为128×128×12的时频空的3阶张量。这里的空指的是导联位置就是指导联轴。对于12导联(lead×time)ECG信号，s[l,n]表示离散信号在时间点n对于导联l。在时间点nΔt和频率f的短时fourier变换如下式：

$\begin{array}{c}\mathrm{STFT}\left\{s\right[l,n\left]\right\}(m,w)=S(l,m,n)\\ ={\mathrm{\Σ}}_{m=-\∞}^{\∞}w(n-m)s(l,m)e^{-j2\mathrm{\πfm}}\end{array}$

这里w[n]是一个窗函数，变换完成后数据变成一个3阶张量。

2.基于TTV映射的张量特征抽取

(1)根据原始张量数据的判别性计算有效的投影张量

$\begin{array}{c}{u}_k^l{|}_{l=1}^M={\arg }_{u_k^l{|}_{l=1}^M}\max {\left(\frac{1}{n}{\mathrm{\Σ}}_{i=1}^c\left((M_i^k-M_k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\mathrm{}}_{\×l}{\left(u_k^l\right)}^T\right)((M_i^k-M_k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\mathrm{}}_{\×l}{\left(u_k^l\right)}^T\right)}^T\\ -{\mathrm{\ζ}}_k^l{\mathrm{\Σ}}_{i=1}^c{\mathrm{\Σ}}_{j=1}^{n_i}\left((X_{\mathrm{ji}}^k-M_i^k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\mathrm{}}_{\×l}{\left(u_k^l\right)}^T\right)\×{\left((X_{\mathrm{ji}}^k-M_i^k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\mathrm{}}_{\×l}{\left(u_k^l\right)}^T\right)}^T)\end{array}$

(2)根据优化张量可分性特点计算更为优越的投影张量

$\underset{i=1}{\overset{C_c^2}{\mathrm{\Σ}}}\underset{j_1=1}{\overset{n_{i1}}{\mathrm{\Σ}}}\underset{j_2=1}{\overset{n_{i2}}{\mathrm{\Σ}}}\left((X_{j_1}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\×}_l{\left(u_k^l\right)}^T\right){\×\left((X_{j_1}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\×}_l{\left(u_k^l\right)}^T\right)}^T+\left((X_{j_2}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\×}_l{\left(u_k^l\right)}^T\right)\×{\left((X_{j_2}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}{\×}_l{\left(u_k^l\right)}^T\right)}^T$

(3)对原始数据去除已抽取投影张量的维度，调整原始数据的结构

$X_{\mathrm{ij}}^k=X_{\mathrm{ij}}^{k-1}-{\mathrm{\λ}}^{k-1}{u}_{k-1}^1\⊗u_{k-1}^2\⊗....\⊗u_{k-1}^M$

通过这个处理之后，后续计算过程就不用考虑，直接计算出来的投影张量的影响，完全可以独立进行运算。

(4)对原始数据计算投影张量进行加权处理。

通过调整每个张量的权重，来规避不合理张量的影响，等式如下：

$S_{\mathrm{oo}}=\underset{i,j}{\mathrm{\Σ}}{w}_{\mathrm{ij}}\underset{x{\∈A}_i,y{\∈A}_j}{\mathrm{\Σ}}w\left(d_{\mathrm{xy}}\right)S_{\mathrm{xy}}i\≠j$

最简单的方式就是取张量间距离的倒数distance(),

w(d_xy)=d_xy ^-n

或者如下式定义:

$w\left(d_{\mathrm{xy}}\right)=\left\{\begin{array}{cc}=1& \mathrm{if}{d}_{\mathrm{xy}}\∈N_{\%}~M_{\%}\\ =0& \mathrm{if}{d}_{\mathrm{xy}}\∉N_{\%}~M_{\%}\end{array}\right.$

组合两种形式得到:

$w\left(d_{\mathrm{xy}}\right)=\left\{\begin{array}{cc}={d_{\mathrm{xy}}}^{-n}& \mathrm{if}{d}_{\mathrm{xy}}\∈N_{\%}~M_{\%}\\ =0& \mathrm{if}{d}_{\mathrm{xy}}\∉N_{\%}~M_{\%}\end{array}\right.$

(5)算法整体迭代收敛过程如下：

整个算法过程，是一个严格的单调收敛过程，逐次迭代目标值关系如下式所示：

$\begin{array}{c}{a}_k=g(u_k^1,1)\≤g(u_k^2,1)\≤...\≤g(u_k^M,1)\≤g(u_k^1,2)\≤g(u_k^2,2)\\ \≤...\≤g(u_k^1,t)\≤g(u_k^2,t)\≤...\≤g(u_k^1,T)\≤g(u_k^2,T)\\ \≤...\≤g(u_k^M,T)\\ =b_k\end{array}$

T->无穷时，算法收敛于最终目标极限收敛值.

(6)判断计算过程结束终止条件

(7)我们使用了如下方法来算法是否收敛，并且判断算法什麽时候应该停止。误差值小于一定阈值：

||F_k-F_k-1||_Fro≤ε

采用这种方法来判断收敛与否和终止整个算法。

3.选择合理初值

我们方法初值选择是否合理，严重影响结果的好坏，和算法的收敛情况。所以，这里使用了一种张量逼近算法，来提供的一个有效的初始张量来提高张量算法的计算结果。基本思路就是求解一个最为近似的张量：

$\min f(a^{\left(1\right)},...,a^{\left(N\right)})\≡\frac{1}{2}{\left|\right|Z-\left[\right[a^{\left(1\right)},...,a^{\left(N\right)}\left]\right]\left|\right|}^2$

(1)无约束情况张量算法

a.交替最小二乘方法

这种方法的求解目标等式如下：

我们可以展开这个等式如下式:

$={\min }_{a^{\left(n\right)}}{\left|\right|Z_{\left(n\right)}-a^{\left(n\right)}{(a^{\left(N\right)}\⊗...\⊗a^{(n-1)}\⊗a^{(n+1)}\⊗...\⊗a^{\left(1\right)})}^T\left|\right|}^2$

这里表示kroneckerproduct，而Z_(n)表示按第nmode转换张量Z一个矩阵。这个问题的解犹如下式：

b.张量梯度下降方法

可以将上述目标函数进行展开，写成如下的形式：

第一项没有涉及变量，所以：

$\frac{\∂f_1}{\∂a^{\left(n\right)}}=0$

这里0表示一个0向量，长度为I_n，第二项犹如如下式子：

$f_2\left(x\right)=Z{\×}_{m=1}^N{a}_r^{\left(n\right)}$

$={\left(Z{\×}_{m=1,m\≠n}^N{a}_r^{\left(m\right)}\right)}^T{a}_r^{\left(n\right)}$

第二项求导之后得到如下式:

$\frac{\∂f_2}{{\∂a}^{\left(n\right)}}=\left(Z{\×}_{m=1,m\≠n}^N{a}_r^{\left(m\right)}\right)$

第三项如下：

$f_3\left(x\right)=\underset{m=1}{\overset{N}{\mathrm{\Π}}}{a^{\left(m\right)}}^T{a}^{\left(m\right)}$

因此,

$\frac{\∂f_3}{\∂a^{\left(n\right)}}=2\left(\underset{m=1,m\≠n}{\overset{N}{\mathrm{\Π}}}{a^{\left(m\right)}}^T{a}^{\left(m\right)}\right)a^{\left(n\right)}$

综合以上三项就可以得到最终结果。

(2)有约束张量情况

a)带约束非线性最小二乘

带约束优化问题，通过转换对求解的约束条件有所放松，再去求解就可以得到更加优越的计算结果。如下式：

$\begin{array}{cc}\min \frac{1}{2}{\mathrm{\Σ}}_{i=1}^l{f_i\left(x\right)}^2& \min \frac{1}{2}{z}^{T}z\\ g_j\left(x\right)=0j=1,...,m_e& f_i\left(x\right)-z_i=0i=1,...,l\\ g_j\left(x\right)\≥0j=m_e+1,...,m& g_j\left(x\right)=0j=1,...,m_e\\ x_l\≤x\≤x_u& g_j\left(x\right)\≥0j=m_e+1,...,m\\ & x_l\≤x{\≤x}_u\end{array}$

将左边的优化问题转换为右边的优化问题，然后将结果代入，使用标准SQP算法求解，来求得最终的优化结果。最终结果可以得到明显改善。

4.分类比对

我们的方法最后使用SVM来对抽取出来的以vector方式存储的向量特征进行分类。

这直接通过求解如下优化主问题：

$\underset{W,b,\mathrm{\ξ}}{\min }C{\mathrm{\Σ}}_{n=1}^N{\mathrm{\ξ}}_n+\frac{1}{2}{\left|\right|w\left|\right|}^2---\left(\mathrm{1.4.1}\right)$

Subjecttoy_i(w^Tφ(x_i)+b)≥1-ξ_n,ξ_n,≥0,i=1,2…,n

这里的参数C>0在松弛变量和惩罚因子间的平衡。而他的Lagrangian乘子变换成如下等式：

$L(w,b,a)=\frac{1}{2}{\left|\right|w\left|\right|}^2+C{\mathrm{\Σ}}_{n=1}^N{\mathrm{\ξ}}_n-{\mathrm{\Σ}}_{n=1}^N{a}_n\{t_{n}y\left(x_n\right)-1+{\mathrm{\ξ}}_n\}-{\mathrm{\Σ}}_{n=1}^N{\mathrm{\μ}}_n{\mathrm{\ξ}}_n---\left(\mathrm{1.4.2}\right)$

这里{a_n≥0}而且{μ_n≥0}是Lagrangian乘子.而对偶Lagrangian问题如下式：

$\tilde{L}\left(a\right)={\mathrm{\Σ}}_{i=1}^N{a}_n-\frac{1}{2}{\mathrm{\Σ}}_{n=1}^N{\mathrm{\Σ}}_{n=1}^N{a}_n{a}_m{t}_n{t}_{m}k(x_n,x_m)---\left(\mathrm{1.4.3}\right)$

它有约束项0≤a_n≤C和.而且k(x,x′)=φ(x)^Tφ(x′)是核函数。

在实际应用中采用了医院方面提供的大约10万条心电作为样本，对本发明的张量心电的方法做了实验评测。该方法和基于向量的方法PCA，ICA和LDA做了比较，同时又和基于张量的方法UMPCA，TR1DA方法做了实验比较。图4-9为不同方法进行特征抽取后降维到3维和的分类效果。可以发现，我们的方法所展示的效果可分性明显优于其它方法。实验分类精度如下所示：图7：PCA(83.4841%)，图5：ICA(81.0724%)，图6：LDA(77.3902%)，图9：UMPCA(86.6387%)，图8：TR1DA(98.4388%)，图4：GTR1DA(99.2140%).实验结果解释我们的方法GTR1DA比其它基于向量的方法明显优越，而对于其它几种基于张量的算法也有不同程度的提升。

Claims (1)

1.一种远程医疗的异常心电张量分析方法，包括下述步骤：

(1)构造张量数据：

a.心电数据采集：

采集标准12导联心电数据；

b.数据预处理和去噪：

对采集的12导联心电数据首先对信号通过50hz的陷波滤波器进行滤波处理，然后对数据进行DB6小波的小波变换分解，然后去除其中最高频的和最低频的信号成份；

c.波形检测：

再次对小波进行DB6小波分解，然后找寻其中第二级的小波系数，然后采用过零点检测方法检测心电的QRS波尖峰，然后去除R波后检测P波和T波尖峰，然后采用detrend算法计算出波形的基线，通过计算各个波形与基线的交点，确定P波QRS波T波的起始和结束，也就是onset和offset；

d.逐跳切割：

一次采集心电大约包含20秒的心电数据，也就是25个完整的由P波QRS波和T波组成的完整心跳波形，对于心电一次完整心跳波形也就是一个P波一个QRS波一个T波进行切割；

e.R波对齐：

对每一个完整心跳波形的心电针对R波的峰值进行对齐，并且切割成统一的长度；

f.短时fourier变换

为了有效抽取时频域的心电诊断特征，通过使用窗口为128点的短时fourier变换对心电进行转换，最后心电被转换为128×128×12的时频空的3阶张量；这里的空指的是导联位置就是指导联轴；对于12导联(lead×time)ECG信号，s[l,n]表示在时间点n，对于第l个导联轴的离散信号值；在时间点nΔt和频率f的短时fourier变换如下式：

$\begin{array}{c}\mathrm{STFT}\left\{s\right[l,n\left]\right\}(m,\mathrm{\ω})\≡S(l,m,m)\\ ={\mathrm{\Σ}}_{m=-\∞}^{\∞}\mathrm{\ω}(n-m)s(l,m)e^{-j2\mathrm{\πfm}}\end{array}$

这里ω[n]是一个窗函数，变换完成后数据变成一个3阶张量；

(2)基于TTV映射的张量特征抽取：

a.根据原始张量数据的判别性计算有效的投影张量：

$\begin{array}{c}{u}_k^l{|}_{l=1}^M={\arg }_{u_k^l{|}_{l=1}^M}\max (\frac{1}{n}{\mathrm{\Σ}}_{i=1}^c((M_i^k-M_k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}\×l{\left(u_k^l\right)}^T)\×{((M_i^k-M_k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}\×l{\left(u_k^l\right)}^T)}^T\\ -{\mathrm{\ζ}}_k^l{\mathrm{\Σ}}_{i=1}^c{\mathrm{\Σ}}_{j=1}^{n_i}((X_{\mathrm{ji}}^k-M_i^k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}\×l{\left(u_k^l\right)}^T)\×((X_{\mathrm{ji}}^k-M_i^k)\underset{l=1}{\overset{M}{\mathrm{\Π}}}\×l{\left(u_k^l\right)}^T{)}^T)\end{array}k$

b.根据优化张量可分性特点计算更为优越的投影张量：

$\underset{i=1}{\overset{C_c^2}{\mathrm{\Σ}}}\underset{j_1=1}{\overset{n_{i1}}{\mathrm{\Σ}}}\underset{j_2=1}{\overset{n_{i2}}{\mathrm{\Σ}}}{((X_{j_1}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}\×l{\left(u_k^l\right)}^T\×{((X_{j_1}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}\×l{\left(u_k^l\right)}^T)}^T+\left({(X_{j_2}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}\×l\left(u_k^l\right)}^T\right){\×((X_{j_2}-M_{j_1{j}_2})\underset{l=1}{\overset{M}{\mathrm{\Π}}}\×l\left(u_k^l\right))}^T)}^T$

对原始数据去除已抽取投影张量的维度，调整原始数据的结构：

$X_{\mathrm{ij}}^k=X_{\mathrm{ij}}^{k-1}-{\mathrm{\λ}}^{k-1}{u}_{k-1}^1\⊗u_{k-1}^2\⊗....\⊗u_{k-1}^M$

c.对原始数据计算投影张量进行加权处理：

通过调整每个张量的权重，来规避不合理张量的影响，等式如下：

$S_{\mathrm{oo}}=\underset{i,j}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{ij}}\underset{x\∈A_i,y\∈A_j}{\mathrm{\Σ}}\mathrm{\ω}\left(d_{\mathrm{xy}}\right)S_{\mathrm{xy}},i\≠j$

或者取张量范数距离的倒数

w(d_xy)＝d_xy ^-n

或者如下式定义：

$w\left(d_{\mathrm{xy}}\right)=\left\{\begin{array}{cc}=1& \mathrm{if}{d}_{\mathrm{xy}}\∈N_{\%}~M_{\%}\\ =0& \mathrm{if}{d}_{\mathrm{xy}}\∉N_{\%}~M_{\%}\end{array}\right.$

或者组合两种形式：

$w\left(d_{\mathrm{xy}}\right)=\left\{\begin{array}{cc}={d_{\mathrm{xy}}}^{-n}& \mathrm{if}{d}_{\mathrm{xy}}\∈N_{\%}~M_{\%}\\ =0& \mathrm{if}{d}_{\mathrm{xy}}\∉N_{\%}~M_{\%}\end{array}\right.$

d.算法整体迭代收敛过程如下：

整个算法过程，是一个严格的单调收敛过程，逐次迭代目标值关系如下式所示：

$\begin{array}{c}{a}_k=g(u_k^1,1)\≤g(u_k^2,1)\≤...\≤g(u_k^M,1)\≤g(u_k^1,2)\≤g(u_k^1,2)\≤(u_k^2,2)\\ \≤...\≤g(u_k^1,t)\≤g(u_k^2,t)\≤...\≤g(u_k^1,T)\≤g(u_k^2,T)\\ \≤...\≤g(u_k^M,T)\\ =b_k\end{array}$

T趋向于无穷时，算法收敛于最终目标极限收敛值；

e.判断计算过程结束终止条件：

使用如下方法来判断算法是否收敛，并且判断算法什么时候应该停止；误差值小于一定阈值则判断算法停止：

||F_k-F_k-1||_Fro≤ε

采用这种方法来判断收敛与否和终止整个算法；

(3)选择合理初值：

求解一个最为近似的张量：

$\min f(a^{\left(1\right)},...,a^{\left(N\right)})\≡\frac{1}{2}{\left|\right|Z-\left[\right[a^{\left(1\right)},...,a^{\left(N\right)}\left]\right]\left|\right|}^{}$

a.无约束情况张量算法：

交替最小二乘方法求解目标等式如下：

其展开等式如下式:

$={\min }_{a^{\left(n\right)}}{\left|\right|Z_{\left(n\right)}-a^{\left(n\right)}{(a^{\left(N\right)}\⊗...\⊗a^{(n-1)}\⊗a^{(n+1)}\⊗...\⊗a^{\left(1\right)})}^T\left|\right|}^2$

这里表示克罗内克积，而Z_(n)表示按模式n展开转换张量Z为一个矩阵，即张量矩阵化；这个问题的解犹如下式：

张量梯度下降方法：

可以将上述目标函数进行展开，写成如下的形式：

第一项没有涉及变量，所以：

$\frac{\∂f_1}{\∂a^{\left(n\right)}}=0$

这里0表示一个0向量，长度为I_n，第二项犹如如下式子：

$\begin{array}{c}{f}_2\left(x\right)=Z{\×}_{m=1}^N{a}_r^{\left(n\right)}\\ ={\left(Z{\×}_{m=1,m\≠n}^N{a}_r^{\left(m\right)}\right)}^T{a}_r^{\left(n\right)}\end{array}$

第二项求导之后得到如下式:

$\frac{\∂f_2}{\∂a^{\left(n\right)}}=\left(Z{\×}_{m=1,m\≠n}^N{a}_r^{\left(m\right)}\right)$

第三项如下：

$f_3\left(x\right)=\underset{m=1}{\overset{N}{\mathrm{\Π}}}{a}^{{\left(m\right)}^T}{a}^{\left(m\right)}$

因此

$\frac{{\∂f}_3}{\∂a^{\left(n\right)}}=2\left(\underset{m=1,m\≠n}{\overset{N}{\mathrm{\Π}}}{a}^{{\left(m\right)}^T}{a}^{\left(m\right)}\right)a^{\left(n\right)}$

综合以上三项就可以得到最终结果；

b.有约束张量情况:

带约束非线性最小二乘：

带约束优化问题，通过转换对求解的约束条件有所放松，再去求解就可以得到更加优越的计算结果，如下式：

$\begin{array}{cc}\min \frac{1}{2}{\mathrm{\Σ}}_{i=1}^l{f}_i{\left(x\right)}^2& \min \frac{1}{2}{z}^{T}z\end{array}$

g_j(x)＝0j＝1，...，m_ef_i(x)-z_i＝0i＝1，...，l

g_j(x)≥0j＝m_e+1，...，mg_j(x)＝0_j＝1，...，m_e

x_l≤x≤x_ug_j(x)≥0j＝m_e+1，...，m

x_l≤x≤x_u

将左边的优化问题转换为右边的优化问题，然后将结果代入，使用标准SQP算法求解，来求得最终的优化结果；

(4)分类比对：

最后使用支持向量机SVM来对抽取出来的以向量方式存储的向量特征进行分类，这直接通过求解如下优化主问题：

$\underset{W,b,\mathrm{\ξ}}{\min }C{\mathrm{\Σ}}_{n=1}^N{\mathrm{\ξ}}_n+\frac{1}{2}{\left|\right|\mathrm{\ω}\left|\right|}^2$

按照相应约束

y_i(ω^Tφ(x_i)+b)≥1-ξ_n，ξ_n≥0，i＝1,2...，n

这里的参数C＞0在松弛变量和惩罚因子间的平衡，而他的拉格朗日乘子变换成如下等式：

$L(\mathrm{\ω},b,a)=\frac{1}{2}{\left|\right|\mathrm{\ω}\left|\right|}^2+C{\mathrm{\Σ}}_{n=1}^N{\mathrm{\ξ}}_n-{\mathrm{\Σ}}_{n=1}^N{a}_n\{t_{n}y\left(x_n\right)-1+{\mathrm{\ξ}}_n\}-{\mathrm{\Σ}}_{n=1}^N{\mathrm{\μ}}_n{\mathrm{\ξ}}_n$

这里{a_n≥0}而且{μ_n≥0}是拉格朗日乘子，而对偶拉格朗日问题如下式：

$\tilde{L}\left(a\right)={\mathrm{\Σ}}_{i=1}^N{a}_n-\frac{1}{2}{\mathrm{\Σ}}_{n=1}^N{\mathrm{\Σ}}_{n=1}^N{a}_n{a}_m{t}_n{t}_{m}k(x_n,x_m)$

它有约束项0≤a_n≤C和而且k(x,x′)＝φ(x)^Tφ(x′)是核函数。