2016, Vol. 37
脉搏波中包含着丰富的人体生理健康等方面的信息，历来是医学、健康等领域的研究热点，目前的脉搏波分析手段局限在已经提出很久的时域及频域的分析手段上，许多研究者正尝试新的方法对脉搏波的复杂性进行研究.本文将分析复杂系统的复杂网络分析方法引入到了脉搏波的分析之中，从复杂系统的角度对脉搏波进行分析，为脉搏波的分析手段提供了一种新的视角.
1 脉搏波分析方法作为人体的重要生理信号，脉搏波一直是中外医学界的重要研究对象.中医对脉搏波的研究虽有丰富的体系结构，却受医生的经验影响较大，不易于传播与继承.西方医学对脉搏波的研究着眼于使用其他领域的理论方法来分析，其研究大致经历了从理论描述到模型分析以及线性化理论到非线性化理论的发展过程[1-2].
脉搏波分析方法主要有时域分析法、频域分析法、时频联合分析法等方法.现有的脉搏波分析方法研究，都是建立在已有的基层分析之上的应用研究，即使用已有的脉搏波特征提取方法获得一些特征量，寻找与这些特征量有关的病理信息以及健康信息；而在脉搏波特征量的提取上，通过新的脉搏波分析方法，从不同的角度分析脉搏波，发现新维度上的脉搏波特征量，是很有意义的.
复杂网络分析方法是对复杂系统的一种分析方法，将复杂网络的分析方法引入到脉搏波的分析之中，从复杂网络的角度分析获得脉搏波的新的特征参量，将丰富脉搏波的分析手段.
2 脉搏波网络化本文在对复杂网络进行一定的研究之后，确定了使用基于水平可见性(horizontal visibility graphs, HVG)的时间序列网络化方法[3]对脉搏波进行分析.
西班牙马德里理工大学的Lacasa等于2008年将复杂网络的分析方法引入到时间序列的分析过程中, 提出了基于可见性(visibility graph, VG)的时间序列网络化方法[4]，并于2009年对该方法进行了简化与改进，提出了HVG网络化方法.近几年内，该方法被引入到了许多新领域的问题的复杂网络分析之中，并取得了较好的研究成果[5-8].已经有人将HVG网络化方法引入到了人体生理信号的分析中，并取得了一些新的成果[9]，而在脉搏波的研究领域，还没有人关注该方法的应用.
VG网络化方法是将时间序列与垂直柱状图以及网络中的节点对应起来，如图 1所示，数据或节点间的关系为柱状图之间的可视关系.可视化的建网规则：任意两个数据A(ta, ya)和B(tb, yb)，具有彼此可视关系的条件是，居于这两个数据之间的任意一个数据C(tc, yc)满足以下条件：
|
(1) |
|
图 1 VG网络化过程示意图 Fig.1 VG network process schematic diagram |
即数据点A(ta, ya)和B(tb, yb)间的任意一个数据C(tc, yc)处于数据点A和B的连线的下方.满足式(1)的数据点间具有可视的特性，相对应的网络中的节点间就相应有一条连边.
HVG网络化方法是在VG网络构建方法上简化而来的一种保持了网络特性的时间序列网络化方法.该方法将VG中的可见性的标准加以限制，将可见性规则限制在垂直柱状条顶端的水平方向上，直观的视图描述如图 2所示，数据与网络中节点间的对应关系继承VG网络构建方法的对应关系.
|
图 2 HVG网络化过程示意图 Fig.2 HVG network process schematic diagram |
在对脉搏波进行分析的基础上，使用脉搏波的周期序列作为网络化的输入数据，即首先使用新的周期自动识别算法识别出脉搏波的波峰位置，如图 3所示，图中绿色标志为极小值位置，红色标志为识别出的极大值位置，黑色标志为在极大值中筛选出的脉搏波波峰位置.
|
图 3 脉搏波周期识别效果 Fig.3 Recognition results of pulse wave period |
以波峰为分割点对脉搏波进行划分，将脉搏波以周期为区间划分为序列，如图 4所示，即可得到脉搏波的周期长度序列.由此使用HVG网络化方法建立脉搏波的复杂网络.
|
图 4 脉搏波周期序列 Fig.4 Pulse wave period sequence |
使用上述方法将脉搏波进行网络化之后，就要对其网络特性进行分析，以获得脉搏波的网络方面的特性.
在实验中，共采集了包括不同性别、不同年龄段的112名志愿者的脉搏波数据，对采集到的数据进行了相应的分析，分析结果显示，脉搏波数据具有相同的规律性.因此，在后续的结果分析中，使用随机选择的一位志愿者的脉搏波数据，作为分析的示例数据.
3.1 度分布度分布是反映网络性质的一个重要标准，它指的是一个网络中节点的度值满足的概率分布，反映出的是节点在网络中的重要程度.使用HVG网络构建方法构建的脉搏波网络的度分布如图 5所示.
|
图 5 HVG网络度分布 Fig.5 Distribution of HVG network |
由于时间序列的开始与结束位置缺少前面或后面的数据，计算出的相应的节点的度值会变小，应该将脉搏波的HVG网络中度值太小的几个节点的度值在度分布中去除.由图 5可知，脉搏波的HVG网络的度分布是一个长拖尾分布.
将以HVG方法构建的脉搏波网络的度分布放在双对数坐标系下，并对脉搏波网络的度分布进行线性拟合，结果如图 6所示.
|
图 6 HVG网络度分布拟合图 Fig.6 Fitting figure of HVG network distribution |
拟合曲线的表达式为
|
(2) |
拟合效果的评价参数如表 1所示.
|
表 1 HVG网络度分布拟合效果 Table 1 Fitting results of HVG network distribution |
由表中参数可知，拟合效果较好，拟合优度达到0.931 5, 均方误差较小，拟合直线能较好反映因变量lgP(k)随自变量lgk的变化关系.
由脉搏波的HVG网络的度分布及其线性拟合的结果可知，构建出的脉搏波网络是度分布符合幂律分布的网络，即网络都是无标度网络.少数节点与前后的较多节点都有影响关系，而大部分节点却与较少的其他节点具有影响关系.从人体状态的角度来看这一特性，可以理解为某一异常脉搏的出现是较长时间的影响累积的结果，并将影响之后的较长时间，而通常脉搏则仅受此前较短时间的影响，并影响此后的较短时间，二者具有一致性.
由脉搏波网络的无标度特性可知，脉搏波在网络化之后即可以使用现有的无标度网络方面的复杂网络研究方法与结论对其进行类似的分析，并使用无标度网络的演化模型，对脉搏波在人体状态发生变化的过程中的变化进行分析，以发现其动态演化规律.
3.2 网络直径为分析对比由脉搏波构建的网络的直径随网络规模的变化而变化的情况，实验中采集了某一志愿者不同时间长度的脉搏波数据，以分析在网络规模不断增加的情况下脉搏波网络直径的变化情况.志愿者的脉搏波采集以一分钟为单位，每个实验样本比前一个样本的采集时长延长一分钟，即实验样本是每次递增一分钟长度的脉搏波样本，上述不同规模脉搏波样本序列共39个样本.由这些脉搏波数据样本进行网络化分析，并计算其网络直径，其网络直径与网络节点数的关系如图 7所示.
|
图 7 网络直径与网络规模关系 Fig.7 Relationship between diameter and network size |
由脉搏波的HVG网络的网络规模与网络直径之间关系的散点图大致可以看出，二者基本上是呈现线性关系的，为验证这种关系，对其进行了线性拟合，拟合结果如图 8所示.
|
图 8 网络直径与网络规模关系拟合图 Fig.8 Fitting figure of the relationship between diameter and network size |
拟合曲线的表达式为
|
(3) |
拟合效果评价参数如表 2所示.
|
|
表 2 网络直径与网络规模关系拟合效果 Table 2 Fitting results of the relationship between diameter and network size |
拟合优度为0.874 1，即两者有较好的线性关系.由拟合效果可见，随着网络规模的增长，网络直径基本上是随之线性增长，因此，由脉搏波构建的该网络不具有小世界的特性.
4 结语本文在对脉搏波的研究现状进行分析的基础上，将HVG网络化方法引入到了脉搏波的分析之中，使用HVG网络化方法构建的脉搏波网络具有无标度的特性，却不具有小世界的特性.复杂网络的分析方法扩展了脉搏波的研究视角，使研究者获得了更多维度的脉搏波特征参量.
| [1] |
王炳和, 相敬林.
脉搏系数建模与脉象信息分析的研究进展[J]. 生物医学工程学杂志 , 2002, 19 (2) : 329–333.
( Wang Bing-he, Xiang Jing-lin. Research progress in analysis of pulse coefficient modeling and pulse condition information[J]. Journal of Biomedical Engineering , 2002, 19 (2) : 329–333. ) |
| [2] |
罗志昌, 张松, 杨益民.
脉搏波的工程分析与应用[M]. 北京: 科学出版社, 2006 : 3 -185.
( Luo Zhi-chang, Zhang Song, Yang Yi-min. Engineering analysis and application of pulse wave[M]. Beijing: Science Press, 2006 : 3 -185. ) |
| [3] | Luque B, Lacasa L, Ballesteros F, et al. Horizontal visibility graphs:exact results for random time series[J]. Physical Review E , 2009, 80 (4) : 046103. DOI:10.1103/PhysRevE.80.046103 |
| [4] | Lacasa L, Luque B, Ballesteros F, et al. From time series to complex networks:the visibility graph[J]. Proceedings of the National Academy of Sciences , 2008, 105 (13) : 4972–4975. DOI:10.1073/pnas.0709247105 |
| [5] | Schleussner C F, Divine D V, Donges J F, et al. Indications for a North Atlantic ocean circulation regime shift at the onset of the Little Ice Age[J]. Climate Dynamics , 2015, 45 (11) : 3623–3633. |
| [6] | Zhang B, Wang J, Fang W. Volatility behavior of visibility graph EMD financial time series from Ising interacting system[J]. Physica A:Statistical Mechanics and Its Applications , 2015, 432 : 301–314. DOI:10.1016/j.physa.2015.03.057 |
| [7] | Aragoneses A, Carpi L, Tarasov N, et al. Unveiling temporal correlations in the intensity of a fibre laser during the transition to optical wave turbulence[J]. Physical Review Letters , 2015, 116 (13) : 1329–1336. |
| [8] | Liu H L, Tang Q, Liu J.Nonlinear analysis of typical hydrological time series under the viewpoint of visibility graph methods and its two modification versions[C]//Proceedings of the 2015 Chinese Intelligent Automation Conference.Berlin:Springer, 2015:327-337. |
| [9] | Zhu G, Li Y, Wen P P. Epileptic seizure detection in EEGs signals using a fast weighted horizontal visibility algorithm[J]. Computer Methods and Programs in Biomedicine , 2014, 115 (2) : 64–75. DOI:10.1016/j.cmpb.2014.04.001 |