CN101822534B - 基于相空间重构的脉搏波观测方法 - Google Patents

Abstract

本发明公开了一种基于相空间重构的脉搏波观测方法，通过对被测对象的脉搏数据的采集分析，利用非线性系统理论，对脉搏信号进行可视化处理，为疾病的诊断提供了一种辅助工具。并通过肝病分类的例子证明了本方法的可行性。本发明利用相空间重构理论，计算脉搏波样本信号的最佳时延与嵌入维，对脉搏信号进行重构，并绘制出Cao氏图，和不同嵌入维下的相空间重构图。找出二维或三维空间中的奇异吸引子。本发明为脉搏波智能诊断系统提供了一种较为简便，易于实现的分析方法。

Description

基于相空间重构的脉搏波观测方法

技术领域

本发明涉及非线性系统理论技术领域，尤其涉及一种基于相空间重构分析方法的脉搏波观测方法。

背景技术

获取脉搏数据，并通过分析脉搏数据来诊断人体疾病，是个方便和有效的途径。20世纪50年代开始至今，随着科学技术的不断发展，作为一种无创伤性诊断的手段和方法，中医脉诊得到了广泛的临床应用。很多学者在对脉象形成的生理学机制以及在生理、病理上的意义进行了研究，在此基础上尝试采用数学方程、物理模型或者波形图等形式来表述常见脉象的特征。另一方面，许多不同学科的研究者运用现代检测技术的方法和手段，将脉象的物理特征描绘、记录下来，对所得到的脉图进行定性和定量相结合的识别和分析。现有的脉搏识别和自动分析的方法，主要是基于脉象图的分析。依赖专家经验，需要建立疾病的专家诊断系统，和庞大的标准模型库支持，使得不容易做成实时分析监测的产品。

发明内容

本发明的目的在于针对现有技术的不足，提供一种基于相空间重构的脉搏波观测方法。为基于脉搏的诊断方法提供了一种新的观测角度。

本发明的目的是通过以下技术方案来实现的：一种基于相空间重构的脉搏波观测方法，包括以下步骤：

(1)采集脉搏波数据：将连续的脉搏波信息根据采样定理转换为一维时间数据序列进行分析；

(2)用互信息法求最佳时延：根据互信息法的局部最小值观点选取最佳时延；

(3)Cao方法求最小嵌入维：在最佳时延下得到Cao氏图中的曲线，从而得到最小嵌入维；

(4)脉搏波时间序列的相空间重构：在延迟时间和嵌入维数确定的情况下，可以将一维时间序列重构为m维的矩阵，其中，m为嵌入维数；

(5)将高维空间中的相空间向量映射到低维空间：将m维的矩阵投影到二维或三维空间上；

(6)在低维空间中寻找吸引子；

(7)特征提取：根据重构出的吸引子形状，利用已知样本组的分类信息，找到能够用于明显区分出两类不同样本的吸引子几何特征上的差异，提取出来作为能反映相应生理信息的指标。如果特征不是很明显，则返回步骤(4)，直到提取出用于分类的吸引子。

本发明的有益效果是：本发明通过对分组被测对象的脉搏数据进行采集分析，利用非线性系统理论，对脉搏信号进行可视化处理，为医生的诊断提供了一种辅助工具。并通过肝病分类的例子证明了本方法的可行性。

本发明利用相空间重构理论，计算脉搏波样本信号的最佳时延与嵌入维，对脉搏信号进行重构，并绘制出Cao氏图，和不同嵌入维下的相空间重构图。找出二维或三维空间中的奇异吸引子。提取出吸引子的特征。作为分类的依据。

以肝病为研究对象，采用本发明中的步骤，对正常组和患病组的脉搏波进行相空间重构。提取一些直观而简便的可视化指标用于疾病的辅助诊断与预警依据。为脉搏波智能诊断系统提供了一种较为简便，易于实现的分析方法。

附图说明

图1为本发明中用到的脉搏波样本图；

图2为本发明中的脉搏波相空间重构流程图；

图3采用cao氏法求最小嵌入维图，其中，(a)为患者最小嵌入维图；(b)为正常人最小嵌入维图；

图4为本发明中采用的两组样本的脉搏波二维重构图，其中，(a)为健康组脉搏波重构图；(b)为肝病组脉搏波重构图；

图5为本发明中采用的两组样本的脉搏波三维重构图，其中，(a)为健康组脉搏波重构图；(b)为肝病组脉搏波重构图；

图6为本发明中采用的两组样本在最佳时延下的脉搏波二维重构图，其中，(a)为健康组脉搏波重构图；(b)为肝病组脉搏波重构图；(c)为肝病组畸变脉搏波重构图；

图7为互信息发求最佳时延示意图。

具体实施方式

1、本发明利用非线性系统分析的理论背景如下：

1.1相空间重构

相空间重构是通过一系列具体算法将隐含在部分变量时间序列中的原系统整体提取出来的方法，也称为系统的动力学重构(Reconstruction of Dynamics)。

对于单变量的时间序列信号x1，x2，...，xN，重构后的相空间为：

Xi＝[xi，xi+τ，...，xi+(m-1)τ]T(1)

其中：i＝1，2，...，L，L＝N-(m-1)τ；

Xi-重构后的相空间矢量；

τ-延迟时间；

m-嵌入维数；

N-原始时间序列点数；

L-重构后相空间矢量个数。

由(1)式可以得到重构后的相空间轨道矩阵：

X1＝[x1，x1+τ，...，x1+(m-1)τ]T

X2＝[x2，x2+τ，...，x2+(m-1)τ]T

XL＝[xL，xL+τ，...，xL+(m-1)τ]T(2)

上述重构相空间的过程相当于将时间序列映射到m维的欧氏空间中，在相空间重构的过程中，嵌入维数m和延迟时间τ是两个重要的参数，它们与重构后的相空间的质量密切相关

1.2吸引子

吸引子指耗散系统经过长时间演化后达到平衡状态，其轨道归于状态空间中的一个有限区域，该区域即称为吸引子。吸引子的几何特性，如维数等，为动力学系统提供了重要的信息。若动态系统的吸引子为分形，则称该吸引子为奇怪吸引子(Strange Attractor)。

1.3互信息法求最佳时延

考虑两离散信息时间序列{s₁，s₂，...，s_n}和{q₁，q₂，...，q_n}构成的系统S和Q。则根据信息论，从两系统测量中所获得的平均信息量，即信息熵分别为：

$\left\{\begin{array}{c}H\left(S\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{P}_s\left(s_i\right){\log }_2{P}_s\left(s_i\right)\\ H\left(Q\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{P}_q\left(q_j\right){\log }_2{P}_q\left(q_j\right)\end{array}\right.$

其中，P_s(s_i)和P_q(q_j)分别为S和Q中事件s_i和q_j的概率。

在给定S的情况下，能获取的系统Q的的信息，即系统S和Q的互信息为：

I(Q，S)＝H(Q)-H(Q|S)

$H\left(Q\right|S_i)=-\underset{j}{\mathrm{\Σ}}\left[\frac{P_{\mathrm{sq}}(s_i,q_j)}{P_s\left(s_i\right)}\right]{\log }_2\left[\frac{P_{\mathrm{sq}}(s_i,q_j)}{P_s\left(s_i\right)}\right]$

$I(Q,S)=\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}{P}_{\mathrm{sq}}(s_i,q_j){\log }_2\left[\frac{P_{\mathrm{sq}}(s_i,q_j)}{P_s\left(s_i\right)P\left(q_j\right)}\right]$

其中，P_sq(s_i，q_j)为事件s_i和事件q_j的联合分布概率。

接着定义[s，q]＝[X(t)，X(t+1)]，其中s代表时间序列X(t)，q为其延迟时间为τ的时间序列X(t+1)，则I(Q，S)显然是与时间延迟τ有关的函数，不妨记为I(t)。I(t)的大小代表了在已知系统S即X(t)，的情况下，系统Q也就是X(t+1)，的确定性的大小。I(t)＝0，表示X(t+1)完全不相关；而I(t)的极小值，则表示x(t)与X(t+1)是最大可能的不相关。采用I(t)的第一个极小值点作为最优时间延迟。

1.4Cao方法求最小嵌入维

将序列X(t_i)(i＝1，2，...，N)构造的m维相空间矢量，即式4-2-1，记为X_i(m)，构造的m+1维相空间矢量记为X_i(m+1)。定义

$a(i,m)=\frac{\left|\right|X_{n(i,m)}(m+1)-X_i(m+1)\left|\right|}{\left|\right|X_{n(i,m)}\left(m\right)-X_i\left(m\right)\left|\right|}$

式中，i＝1，2，...，N-m；X_n(i，m)(m)是离轨线；X_i(m)最近的轨线；n(i，m)是满足条件1≤n(i，m)≤N-m的正整数且依赖变量i和m；‖‖表示欧氏距离下的最大值范数，即

‖X_k(m)-X_i(m)‖＝max_{0≤j≤(m-1)}|X_k+jt-X_i+jt|

记a(i，m)的均值为

$E\left(m\right)=\frac{1}{N-\mathrm{mt}}\underset{i=1}{\overset{N-\mathrm{mt}}{\mathrm{\Σ}}}a(i,m)$

这里的E(m)独立于变量嵌入维数m和时间延迟τ，为了找到从m到m+1变化的最佳嵌入维数，定义

E₁(m)＝E(m+1)/E(m)

如果时间序列所描述的是动力系统的混沌现象，当E₁(m)自某个m₀开始停止变化，则m₀+1即为所寻找的最佳嵌入维数。

1.5计算关联维数

本发明中采用G-P算法计算关联维数。该算法解决了从有限长离散时间序列及有限长度离散距离向量中估计关联维大小的问题。在已知延时τ，嵌入维数m情况下，从原数据中得到重构相空间中的向量点：

X_i(d)＝(x_i，x_i+τ，...，x_i+(m-1)τ)，i＝1，2，...，N-(d-1)τ

关联积分定义为如下函数，其中θ(x)为Heaviside函数。

$C(m,N,r,t)=\frac{2}{M(M-1)}\underset{1\≤i\≤j\≤M}{\mathrm{\Σ}}\mathrm{\θ}(r-d_{\mathrm{ij}});r>0,d_{\mathrm{ij}}=\left|\right|X_i-X_j\left|\right|$

数学上，关联维定义为：

$D(m,t)=\underset{r\→0}{\lim }\frac{\log C(m,r,t)}{\log r}$

式中： $C(m,r,t)=\underset{N\→\∞}{\lim }C(m,N,r,t)$

鉴于工程中无限长序列无法得到且精度上r不能无限小，G-P算法提出用一个关联积分曲线中某个线性区域的斜率来近似估计这个关联维；原因是当r过大时，C恒等于1，表现为斜率恒为0，当r过小时，积分曲线对系统中微小变化十分敏感，导致噪声影响明显，表现为临近点变化率大幅跳变，线性区域则表示r大小合适，符合系统内在特性要求[2]。故关联维求取式变为：

$D\left(m\right)\≈\frac{\mathrm{\Δ}\log C(N,r)}{\mathrm{\Δ}\log r}$

本发明利用非线性系统理论，对脉搏信号进行分析与处理，包括以下步骤：

(1)采集脉搏波数据：

将连续的脉搏波信息根据采样定理转换为一维时间数据序列进行分析。

(2)用互信息法求最佳时延：

如(图7)中，根据互信息法的局部最小值观点，在0.5-3之间时，互信息函数呈递减，3-6之间时，互信息函数呈递增，局部极值为3，选取最佳时延τ＝3。

(3)Cao方法求最小嵌入维：

如(图3)所示：在最佳时延下，根据1.4中的算法，得到Cao氏图中的曲线。以(图3b)为例，E2随嵌入维数的变化发生明显的变化，这说明数据具备非线性特性。根据E1曲线的变化情况，取最小嵌入维m＝20。

(4)脉搏波时间序列的相空间重构

根据1.1中的算法，在延迟时间和嵌入维数确定的情况下，可以将一维时间序列重构为m维的矩阵。以上面最小嵌入维为20的数据序列为例，一个长度为3000的序列被重构为20维的矩阵，该矩阵的每行是相空间中的一个点向量。结果说明该脉搏信号重构所需维数高于正常的三维。

(5)将高维空间中的相空间向量映射到低维空间

由于维数已经超过了便于人眼识别的三维理解范围，为了将这个空间形象化，同时也更加利于分析和观察，统一将这个空间投影到二维或三维空间上，为了方便结果的统一比较和分析，重构时，横轴表示第一个向量点，纵轴表示第二个向量点(类推)。据此作出相空间中的投影图。如(图6)所示。

(6)在低维空间中寻找吸引子

在相空间重构中，m值取值是否合理对吸引子的获取起到至关重要的作用，当嵌入维m取值不合理的时候，反映在投影图上，就是相空间中投影图的扭曲变形，是对原混沌吸引子不真实的体现。为了消除在实验进行中的一些误差，需要对投影在二维或三维空间上的吸引子进行检验。如果不能将混沌吸引子充分打开，并且体现出一类脉搏波的自相似性。那么对相空间重构的参数进行调整，直到能正确反映系统特征的吸引子出现。

(7)特征提取

根据重构出的吸引子形状，提取出能反映相应生理信息的指标。总结出脉搏波相空间重构的实际意义。

下面根据附图和实施例详细说明本发明，本发明的目的和效果将变得更加明显。

实施例1

针对具体脉搏波序列对上述步骤适当裁剪，以肝病为研究对象，采集两组脉搏波数据进行相空间重构。步骤如下：

(1)选择测试样本。采集二组实验数据作对照分析，以下是这二组数据的具体情况：第一组是健康正常的年轻男性和女性脉搏数据。第二组是采自浙江省第一人民医院的，二十例肝病患者的脉搏数据，包括男女老少和不同疾病。

利用HK-2000C\D型集成化脉搏传感器进行稳态下的人体左手桡动脉脉搏信号的提取，即200Hz采样的脉搏波时间序列。

(2)提取用于分析的脉搏波时间序列。分别选取正常组9人，和患有不同程度肝病患者6人的脉搏数据。由于测量时一些外界因素的影响，有些测量数据因为被测量者的轻微的动作例如说话，数据失真，所以要将数据先复原为波形，然后从中截取没有失真的波段的数据进行分析，为了方便数据的处理和对比，所有数据都截取3000个。一共有15列3000个脉搏波时间序列。

采用Matlab(The MathWorks，Inc.)作为数据处理软件。版本为R2009a。将脉搏采集数据导入工作空间(Workspace)中，进行后续分析。

(3)用互信息法求最佳时延τ_d

调用根据互信息法编写的MATLAB求时延函数，分别求15列脉搏波数据的最佳时延。结果见表1中的最佳时延纵列。

(4)Cao氏方法求最小嵌入维m

调用根据Cao氏法编写的MATLAB求最小嵌入维函数，分别计算15列脉搏波数据的最小嵌入维并画出Cao氏图。最小嵌入维结果见表1中的最小嵌入维纵列。为了便于分析说明，这里选取了两组数据中较为典型的Cao氏图作为对照。(见图3)正常人的Cao氏图中在达到最小嵌入维后，E1和E2都趋于稳定，但是肝病患者的多波折，变化幅度比较大。如图3。其中(a)图为典型患者最小嵌入维图，(b)图为典型正常人最小嵌入维图。

(5)画出时延1的二维和三维相空间重构图

5.1取τ＝Δt，m＝2，重构出一个2×(n-Δt-1)的2维向量。n为原始一维时间序列的长度。Δt为一个采样间隔，m为嵌入维。画出二维相空间重构图。(如图4)a组为肝病组的脉搏波重构图，b组为正常组的男女典型脉搏波重构图。b组中的二维脉搏波重构图有很强的相似性，故分别挑选男女两个典型重构图代表其整体特征。

5.2取τ＝Δt，m＝3，重构出一个3×(n-Δt-1)的3维向量。n为原始一维时间序列的长度。Δt为一个脉搏信号采样间隔，m为嵌入维。画出三维相空间重构图。(如图5)a组为肝病组的脉搏波重构图，b组为正常组的男女典型脉搏波重构图。b组中的三维脉搏波重构图有很强的相似性，故分别挑选男女两个典型重构图代表其整体特征。

(6)画最佳时延τ_d下的二维相空间重构图

取τ＝τ_d，m＝2，重构出一个2×(n-τ_d-1)的2维向量。n为原始一维时间序列的长度。τ_d为最佳时延，m为嵌入维。画出二维相空间重构图。(如图6)a组图为正常组的男女典型脉搏波重构图，b组为典型肝病组的脉搏波重构图。c组为肝病组中有一定程度畸变的二维重构图。造成分类上的困难，故提取出来。

(7)观测吸引子差异

在时延为τ＝Δt时，男、女性样本的脉搏数据的重构图均为棒球棒型。但，男性样本的曲线聚合度明显高于女性样本，即重构图的棒型的内部浅色区域的面积大于女性的重构图的浅色区域面积。男性脉搏的重构图为“平行”状，女性样本脉搏重构图的特征则与男性样本有明显不同，为“交叉”状。”(图4b)中左图为女性典型二维重构图；右图为男性典型二维重构图。

两组脉搏波数据的二维相空间重构图(如图6)。其中在健康组脉搏波重构图中随机抽取了4组二维重构图(图6a)作为对照。对比可以看出：与正常人二维重构图相比，这些肝病患者(图6b)的重构图，所有的线都相对集中，而这些患者几乎都是肝病的重病者，比如肝癌、乙肝。

例外的是：另外两个患有肝硬化的患者没有出现这样的现象：如图6(c)。

(8)特征提取

8.1罹患肝病的患者出现了非常明显的最小嵌入维变小的特点。

8.2检测个体如果有肝脏疾病，那他的脉搏数据测量结果的二维相空间重构图可能反映出棒状中心空白面积变小，棒形变细的特点。

本发明中的提出的推论并不适用来用来肝病的诊断，这与本发明中的试验为探索性而非证明性试验有关。即并不能通过相空间重构的结果来确认究竟是否罹患肝病。而是通过该病例的分析证明相空间重构的方法可以作为医生或智能设备用来辅助诊断的一种工具。

表1：对两组样本分析后得到的数据表

                健康人的脉搏参数

序号	状况	心率(次/分)	最佳时延	最小嵌入维
					1	正常	64	24	19
2	正常	74	26	16
					3	正常	72	27	15
4	正常	64	34	16
					5	正常	72	26	17
6	正常	76	31	15
					7	正常	66	34	15
8	正常	90	24	15
					9	正常	100	28	20

                肝病患者的脉搏信号参数

编号	疾病	心率(次/分)	最佳时延	最小嵌入维
					10	肝硬化	68	16	14
11	肝硬化	88	23	15
					12	乙肝	76	25	15
13	肝病	48	27	大于30
					14	血栓、肝癌	94	26	10
15	脾肿大、乙肝	66	27	17

上述实施例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。

Claims (1)

1.一种基于相空间重构的脉搏波观测方法，其特征在于，包括以下步骤：

(1)采集脉搏波数据：将连续的脉搏波信息根据采样定理转换为一维时间数据序列进行分析；

(2)用互信息法求最佳时延：根据互信息法的局部最小值观点选取最佳时延；

(3)Cao方法求最小嵌入维：在最佳时延下得到Cao氏图中的曲线，从而得到最小嵌入维；

(4)脉搏波时间序列的相空间重构：在延迟时间和嵌入维数确定的情况下，可以将一维时间序列重构为m维的矩阵，其中，m为嵌入维数；

(5)将高维空间中的相空间向量映射到低维空间：将m维的矩阵投影到二维或三维空间上；

(6)在低维空间中寻找吸引子；

(7)特征提取：根据重构出的吸引子形状，利用已知样本组的分类信息，找到能够用于明显区分出两类不同样本的吸引子几何特征上的差异，提取出来作为能反映相应生理信息的指标。如果特征不是很明显，则返回步骤(4)，直到提取出用于分类的吸引子。

Images