周树然¹, 姚阳¹, 郝丽玲¹, 徐礼胜^1,2
1. 东北大学 中荷生物医学与信息工程学院, 辽宁 沈阳 110169;
2. 东北大学 教育部医学影像计算重点实验室, 辽宁 沈阳 110169
收稿日期：2016-10-20
基金项目：国家自然科学基金资助项目(61773110, 61374015, 61202258)；中央高校基本科研业务费专项资金资助项目(N161904002，N130404016)。
作者简介：周树然(1989-), 男, 内蒙古呼和浩特人, 东北大学博士研究生;
徐礼胜(1975-), 男, 安徽安庆人, 东北大学教授, 博士生导师。

摘要：针对弹性管模型适用于模拟大中型动脉血管, 不适用于模拟微小型动脉血管问题, 对人体上肢血管系统建立一种弹性管与弹性腔融合模型, 从而提高其模型的精度.首先分别在肱-桡动脉段建立弹性管模型, 在桡-指动脉段建立弹性腔模型; 然后在两模型的接触面处应用血管压力和血液流量的连续性原理, 建立二者联系, 导出融合模型; 最后对比分析弹性管模型和弹性管与弹性腔融合模型的仿真性能.在PRD, H1, H2, t1和t 2评价指标下, 弹性管与弹性腔融合模型(PRD:(4.0±2.0)%; t1:(2.5±1.7)%; t2:(3.1±1.9)%; H1:(1.9±2.1)%; H2:(1.3±1.3)%)的仿真波形和实测波形误差均低于弹性管模型(PRD:(16.5±8.3)%; t1:(3.4±2.4)%; t2:(4.5±3.7)%; H1:(22.1±14.9)%; H2:(19.7±15.3)%).结果表明, 弹性管与弹性腔融合模型能够更精确地模拟人体上肢血管系统的生理状态.
关键词：心血管系统脉搏波弹性腔模型弹性管模型融合模型
Modeling Simulation of Upper Limb Vascular System Based on Fusion of Windkessel and Tube-Load Models
ZHOU Shu-ran¹, YAO Yang¹, HAO Li-ling¹, XU Li-sheng^1,2
1. School of Sino-Dutch Biomedical & Information Engineering, Northeastern University, Shenyang 110169, China;
2. Key Laboratory of Medical Image Computing, Ministry of Education, Northeastern University, Shenyang 110169, China
Corresponding author: XU Li-sheng, E-mail: xuls@bmie.neu.edu.cn
Abstract: In response to the problem that the tube-load model is suitable to simulate large and medium-sized arteries, not suitable to simulate small arteries and arterioles, the model combining windkessel model with tube-load model was established to improve the simulation accuracy for the human upper limb vascular system. Firstly, the tube-load model was applied to the simulation from brachial artery to radial artery, and the windkessel model was established in the section from radial artery to finger artery, respectively. Secondly, based on the continuity principle of blood pressure and blood flow in the interface of two models, the relationship of the two models was built and a fusion model was derived by combining the two models. Finally, the simulation performances of the tube-load model and the fusion model were compared. With the evaluation indexes, such as PRD, H1, H2, t1 and t2, the error between the measured waveform and the simulated waveform of the fusion model (PRD: (4.0±2.0) %; t1: (2.5±1.7) %; t2: (3.1±1.9) %; H1: (1.9±2.1) %; H2: (1.3±1.3)%) is lower than that of the tube-load model (PRD: (16.5±8.3)%; t1: (3.4±2.4) %; t2: (4.5±3.7) %; H1: (22.1±14.9) %; H2: (19.7±15.3) %). The results show that the fusion model can accurately simulate the physiological state in the human upper limb vascular system.
Key Words: cardiovascular systempulse wavewindkessel modeltube-load modelfusion model
血管系统是一个复杂精细的生命系统, 维持着人体生命的正常运行.脉搏波作为人体基本的生理信号, 蕴含着丰富的心血管系统的生理病理信息, 通过采集分析脉搏波信号, 可以用来预测或诊断人体心血管系统的健康状态.
随着计算机技术的快速发展, 基于脉搏波信号的建模仿真法逐渐成为研究心血管系统的主要方法.建模仿真法分为集中参数模型和分布参数模型.集中参数模型是依据弹性腔理论, 建立心血管系统的弹性腔模型.在19世纪末, Frank^[1]首次提出了弹性腔理论, 建立了弹性腔模型, 利用一个大的弹性腔来模拟人体的主要动脉血管.之后Goldwyn等^[2]建立了双弹性腔模型, 利用两个相连的弹性腔来模拟动脉血管.何为等^[3]将心血管系统集中参数的建模细分为高度集中参数和精细集中参数两个方向.弹性腔模型从总体上刻画了心血管系统的集中特性, 却无法体现脉搏波的传播性质.分布参数模型是利用脉搏波的传播理论, 建立心血管系统的弹性管模型.Womersley^[4]提出并求解血流控制的线性方程, 研究了弹性直管、刚性直管、分叉血管以及异径直管等情况.Ling等^[5]提出非线性局部流动理论, 考虑了血管壁的大变形和非线性弹性特性以及保留了血流运动方程中对流加速度项.Alastruey等^[6]将分布参数模型扩展, 构建全身的动脉网络来探讨脉搏波的传播机制.弹性管模型虽然可以反映脉搏波的传播性质, 但不适合用来模拟微循环系统.
因此, 本文基于上肢血管系统的真实生理特性来综合两种模型的优点, 应用心血管系统的血流动力学原理并结合解剖学知识, 建立肱-指动脉的弹性管与弹性腔融合模型, 从而对人体上肢血管系统进行更加深入地研究, 为医疗设备研发和临床诊断提供新的思路与理论依据.
1 模型的建立1.1 肱-桡弹性管模型在弹性管理论中, 一般假定动脉血管中流动的血液是密度为常数的理想流体, 同时不再考虑血液的黏性特征, 也不考虑小动脉对大动脉的作用, 利用弹性管两端的纯流阻来表示^[7].通过上述的假设和简化, 最终可以得到从主动脉弓分叉出来到指端动脉为止的臂动脉是1条横截面积恒定和管壁弹性模量均匀的弹性管.因此, 可以建立描述动脉管路中血液流动的线性方程为

(1)

式中:L为血管单位长度流感; C为血管单位长度流容; p为弹性管任一点的瞬时血压; q为弹性管任一点的瞬时血流量.
对式(1)进行拉普拉斯变换, 即可得到此波动方程的一般解.由于脉搏波信号是一个不规则的连续信号, 无法用一个已知的函数来表示, 所以需要先将其离散化, 然后用已知函数近似代替.本文利用多个小阶梯来逼近真实脉搏波信号, 当弹性管近端的输入压力为脉搏波p_A0(t)时, 可以将它考虑为N个小压力阶跃p_0i组成, 其表达式为

(2)

最后, 按照线性叠加理论, 可得到弹性管模型的血压和血流量传输方程:

(3)

其中:β₀为弹性管近端反射系数; β_l为弹性管远端反射系数; C₀为脉搏波波速; l为臂动脉长度; Z_c为弹性管特性阻抗; U(t)为单位阶跃函数.
肱-桡弹性管模型是对肱动脉到桡动脉的动脉血管段建立弹性管模型.将肱动脉处的实测压力脉搏波作为弹性管模型的输入p_A0(t), 并将不同的力学参数β₀, β_l和C₀代入上述的传输方程, 则可在桡动脉处得到对应的压力脉搏波波形^[6].
1.2 桡-指弹性腔模型在弹性腔理论中, 一般将动脉看作是一个弹性腔, 假定对于某一特定时刻, 腔内的各种流动参数如血压和血流等处处相同, 它们只是时间的函数, 与位置无关, 并且假定弹性腔内的体积与弹性腔内的压力呈正比例关系^[8].
桡-指弹性腔模型是对桡动脉到指端动脉的动脉血管段建立弹性腔模型.桡-指弹性腔模型由集中参数R, L, C二阶系统构成.根据等效电路原理, 并应用质量守恒与动量守恒定律, 可得到桡-指弹性腔容积脉搏波的方程式如下:

(4)

从而能够导出其微分方程:

(5)

其中:将桡动脉处的压力脉搏波作为模型输入p_in, 其由压力脉搏传感器测得; 将指端处的容积脉搏波作为模型输出q_out, 其由光电容积脉搏传感器测得.根据上述数学表达式, 若将表示不同生理条件的微循环模型参数R, L, C代入模型中, 并给定桡动脉处压力脉搏波p_in, 则可以在指端处得到对应模型参数的容积脉搏血流波形q_out.
1.3 肱-指弹性管与弹性腔融合模型弹性管模型是偏微分方程, 它的脉搏波波形是随时间t和位置x两个变量一起变化的, 而弹性腔模型是常微分方程, 它的脉搏波波形只是时间t的函数, 与位置x无关.通过两个模型(弹性腔模型和弹性管模型)的融合, 消去它们在空间维度上的差异, 实现了两种尺度模型的建立^[9].在模型融合的过程中, 应用血管压力和血液流量在弹性管模型与弹性腔模型接触面处的连续性原理, 建立二者的联系, 具体关系表达式如下:

(6)

(7)

(8)

式(6), 式(7)中:上角标代表两个模型的界面, 即(1D→0D)代表弹性管模型到弹性腔模型的界面处; 下角标代表不同的模型类型, 即0代表弹性腔模型, 1代表弹性管模型; p代表血管压力; q代表血液流量; l代表血管段长度.在弹性管模型与弹性腔模型的界面处, 把弹性管模型终端的血管压力p和血液流量q的输出波形作为弹性腔模型始端的输入波形.式(8)中:上角标代表不同的血管段, 即b代表肱动脉，r代表桡动脉; l_br代表肱动脉到桡动脉的距离.
弹性管模型和弹性腔模型的结合发生在两个模型相连的血管段处, 将肱动脉作为模型的起始端, 把指端动脉作为模型的终端, 中间经过的是距离肱动脉始端处长度l_br的桡动脉.在弹性管模型的始端处输入实测的肱动脉血压波形, 然后压力信号沿着弹性管传输到弹性管末端即桡动脉处, 将此时弹性管输出的压力波形作为桡-指弹性腔模型的输入信号, 经过弹性腔的作用, 最终在弹性腔的输出端可得到指端的脉搏波信号, 见图 1.
图 1(Fig. 1)

图 1 肱-指弹性管与弹性腔融合模型Fig.1 Fusion model from brachial artery to finger artery

2 模型的评价2.1 数据采集本实验以东北大学在校大学生为研究对象, 共采集24人, 其中男生12人, 女生12人, 年龄范围为23~28岁, 均为健康的年轻者, 无不良过往病史.脉搏数据采集的主要设备为成都泰盟科技有限公司研发生产的BL-420多通道生理信号记录仪.在采集脉搏信号的同时并记录了每位受试者的年龄、身高、体重、血压等主要生理信息.
本实验的采集过程共分为三个步骤:第一步, 测量血压.受试者在平静状态下处于坐姿, 取三次血压平均值为最终血压值.第二步, 采集指端动脉处容积脉搏波信号.采样频率为1 kHz, 采集时间为3 min.第三步, 同步采集肱动脉、桡动脉、指端动脉三处压力脉搏波信号.采样频率为1 kHz, 采集时间为3 min.第四步, 测量距离.利用软尺分别测量三处动脉之间的距离.
2.2 模型的评价方法本文利用波形相似度的方法对模型的输出波形给予评价.在进行波形相似度的比较时, 选择的是均方根百分误差^[10]和特征参数^[11]两个指标.其中, 均方根百分误差是从全局考虑, 反映的是波形整体的性质; 而特征参数指标, 则从局部出发, 表明的是波形的细节特征.它们各从不同的角度刻画了波形的相似度, 能全面反映波形性能, 以下为两种方法的具体计算过程:
1) 均方根百分误差评价指标.均方根百分误差表示的是实测波形与模型输出波形二者误差的百分比, 具体计算方法为

(9)

式中:u(i)为实测波形; v(i)为仿真波形.
2) 特征参数评价指标.脉搏波波形特征参数主要分为时间类特征参数和幅值类特征参数两类.本文采用了4个主要的特征参数来进行分析, 前两个参数为时间类参数, 后两个参数为幅值类参数, 其计算方法见表 1.
表 1(Table 1)

表 1 主要的特征参数Table 1 Main characteristic parameters 符号 波形特征参数 计算方法 t1 主波起始点到降中峡最低点时间 降中峡最低点时间与主波起始时间的差值 t2 主波起始点到重搏波最高点时间 重搏波最高点时间与主波起始时间的差值 H1 降中峡高度 降中峡幅值与降支最小值的差值 H2 重搏波高度 重搏波幅值与降支最小值的差值

表 1 主要的特征参数 Table 1 Main characteristic parameters

3 结果与分析3.1 仿真波形结果在弹性管模型和融合模型的始端, 输入同一实测肱动脉脉搏波, 经过两种模型的传输, 则可在各自模型终端得到相应的指端动脉仿真波形.由图 2可见, 输出指端动脉波形与实测指端动脉波形, 从主波到切迹部分存在较大偏差, 从切迹到重搏波部分存在较小偏差, 从重搏波到波形终点有一定波动, 其余部分拟合良好.由图 3可见, 输出指端动脉波形与实测指端动脉波形, 只在切迹到重搏波部分存在较小偏差, 其余部分拟合良好.
图 2(Fig. 2)

图 2 指端动脉的实测波形与弹性管模型仿真波形对比Fig.2 Comparison of the measured waveform and the waveform simulated by the tube-load model at finger artery

图 3(Fig. 3)

图 3 指端动脉的实测波形与融合模型仿真波形对比Fig.3 Comparison of the measured waveform and the waveform simulated by the fusion model at finger artery

3.2 统计分析结果对弹性管模型和融合模型, 在PRD, t1, t2, H1, H2指标下相应实测与仿真波形的误差进行统计分析.由图 4可知, 两种模型的PRD, H1, H2误差值相差较大, 而t1, t2误差值相差较小.由表 2可知, 融合模型指标的误差值均低于弹性管模型对应指标的误差值.对PRD, t1, t2, H1和H2指标, 融合模型比弹性管模型的误差值分别降低了12.5%, 0.9%, 1.4%, 20.2%和18.4%.
图 4(Fig. 4)

图 4 弹性管模型和融合模型的不同指标的误差分布对比Fig.4 Comparison of error distribution of different indexes between the tube-load model and the fusion model (a)—均方根百分误差; (b) —主波起始点到降中峡最低点时间; (c) —主波起始点到重搏波最高点时间; (d) —降中峡高度; (e) —重搏波高度.

表 2(Table 2)

表 2 弹性管模型和融合模型的不同指标误差的统计结果Table 2 Statistical results of the error of different indexes between the tube-load model and the fusion model % 评价指标 弹性管模型(肱-指) 融合模型(肱-指) PRD 16.5±8.3 4.0±2.0 t1 3.4±2.4 2.5±1.7 t2 4.5±3.7 3.1±1.9 H1 22.1±14.9 1.9±2.1 H2 19.7±15.3 1.3±1.3 ??注：表中数据均采用mean±SD(均值±标准差)形式进行表示.

表 2 弹性管模型和融合模型的不同指标误差的统计结果 Table 2 Statistical results of the error of different indexes between the tube-load model and the fusion model

融合模型相比弹性管模型误差减少, 是由于在忽略尺动脉的情况下, 肱动脉到桡动脉部分主要由1根中型动脉构成, 细的分支较少, 主要集中1根血管, 可以用弹性管模型来近似; 而桡动脉到指端动脉部分处于大量微动脉与微静脉的交汇处, 聚集着大量的毛细血管, 此部分属于手掌的微循环, 不能再用1根弹性管来近似, 这样带来较大误差, 也不符合实际情况^[7].在人体上肢动脉血管段处, 弹性管模型适合模拟肱动脉到桡动脉部分, 不适合模拟肱动脉到指端动脉部分, 所以应选用能反映微循环生理特性的弹性腔模型来仿真.
4 结语通过上述指标分析, 融合模型的均方根百分误差是最小的, 而且切迹幅值与重搏波幅值也与实测波形最接近, 很好地模拟了波的反射情况.因此, 对于人体上肢血管系统而言, 融合模型优于弹性管模型, 更加接近上肢血管系统的真实状态.
未来建模的发展趋势是多尺度模型的结合, 在接下来的研究中, 将2维或3维模型与0维模型和1维模型结合, 构成三种或四种尺度模型, 充分整合各种模型的优点, 进一步提高模型的精度.同时将建模的范围扩展到更多的人体血管甚至全身各处, 建立包括心脏在内的完整循环系统.
参考文献

[1]	Frank O. Die grundform des arteriellen pulses[J].Zeitschrift für Biologie, 1899, 37(19): 483–526.
[2]	Goldwyn R M, Watt T B. Arterial pressure pulse contour analysis via a mathematical model for the clinical quantification of human vascular properties[J]. IEEE Transactions on Biomedical Engineering, 1967, 14(1): 11–17.
[3]	何为, 余传祥. 心血管动力学参数测量原理与临床应用[M]. 北京: 科学出版社, 2010: 276-304. ( He Wei, Yu Chuan-xiang. The measuring principle and clinical application of cardiovascular dynamics parameters[M]. Beijing: Science Press, 2010: 276-304.)
[4]	Womersley J R. Oscillatory flow in arteries:the constrained elastic tube as a model of arterial flow and pulse transmission[J].Physics in Medicine & Biology, 1957, 2(2): 178–187.
[5]	Ling S C, Atabek H B. A nonlinear analysis of pulsatile flow in arteries[J].Fluid Mechanics, 1972, 55(3): 493–511.DOI:10.1017/S0022112072001971
[6]	Alastruey J, Khir A W, Matthys K S, et al. Pulse wave propagation in a model human arterial network:assessment of 1-D visco-elastic simulations against in vitro measurements[J].Journal of Biomechanics, 2011, 44(12): 2250–2258.DOI:10.1016/j.jbiomech.2011.05.041
[7]	罗志昌, 张松, 杨益民. 脉搏波的工程分析与临床应用[M]. 北京: 科学出版社, 2006: 8-155. ( Luo Zhi-chang, Zhang Song, Yang Yi-min. Engineering analysis and clinical applications of pulse wave[M]. Beijing: Science Press, 2006: 8-155.)
[8]	王昊, 伏全海, 徐礼胜, 等. 基于血管弹性腔模型的脉搏波发生器的现场可编程门阵列设计与实现[J].生物医学工程学杂志, 2014, 31(5): 989–993. ( Wang Hao, Fu Quan-hai, Xu Li-sheng, et al. Design and implementation of the pulse wave generator with field programmable gate array based on windkessel model[J].Journal of Biomedical Engineering, 2014, 31(5): 989–993.DOI:10.7507/1001-5515.20140186)
[9]	Liu H, Liang F, Wong J, et al. Multi-scale modeling of hemodynamics in the cardiovascular system[J].Acta Mechanica Sinica, 2015, 31(4): 446–464.DOI:10.1007/s10409-015-0416-7
[10]	Nguyen P, Kim J M. Adaptive ECG denoising using genetic algorithm-based thresholding and ensemble empirical mode decomposition[J].Information Sciences, 2016, 373(1): 499–511.
[11]	Huotari M, M??tt? K, R?ning J. Photoplethysmographic measurements of arterial and aortic pulse waveform characteristics[J].Finnish Journal of Ehealth & Ewelfare, 2015, 7(2/3): 83–87.