脉搏波本构关系实验研究的探索*

王礼立，王 晖，丁圆圆，陈霞波，杨黎明，龚文波，浣 石，缪馥星

（1. 宁波大学冲击与安全工程教育部/浙江省重点实验室，浙江宁波 315211；2. 宁波市中医院王晖工作室，浙江宁波 315000；3. 季华实验室，广东 佛山 528200）

由连续介质力学波动理论[1]，按中医整体观思路，我们为脉搏波系统建立了一个等价的一维波力学模型[2-3]。问题的控制方程组由如下3 个守恒方程（动量守恒、质量守恒和能量守恒）加上系统本构方程（广义的状态方程）组成。

动量守恒方程为：

式(1)～(4)中：X为物质坐标，t为时间；p为压力，V为比容（与密度互为倒数），v为质点速度，e为比内能，均为(X，t)函数。本构方程p=p(V) 一般为非线性函数。

对上述脉搏波系统的控制方程组（式(1)～(4)）进行分析，并展开进一步研究[2-6]后，可知：

(1)脉搏现象包含着血液的流动和在血液中传播的携带能量（内能e+动能）及整体信息的脉搏波。前者是人眼易见的实体血液本身的物质流；后者是人眼不易见、以波形式传播的能量信息流，可以分别理解为中医的“血”和“气”。中医的“气”是什么？目前并无共识。按上述分析并联系中医气血观可知：血是气传播的物质载体（媒介），气是以波形式传播、携带能量和生命体整体信息、推动血运行的波；即脉搏波不只是反映心脏、血液和血管等循环系统本身局部实体器官的信息。血液物质流与脉搏波能量信息流两者既互有区别，又密切相关耦合。这与中医“气为血之帅，血为气之母”的观点一致。

(2)脉搏波所包含的压力波p(X,t)、质点速度波v(X,t)、比容波V(X,t)和内能波e(X,t)等不同形式的扰动，以波速C传播：

式中：K(p)= - dp/(dV/V0) 为非线性体积压缩模量。可见，波速C由非线性本构方程p=p(V) 的局部斜率（dp/dV）决定（见图1），随p而变，不是恒值。要注意区分脉搏波传播速度C和血液质点的流动速度v，前者（100～101m/s 量级）以比后者（10-1m/s 量级）快得多的速度，把心脏施加给循环血液的扰动（脉动载荷），携带着能量和信息，由近及远地传递到生命体的各部分，不到1 s 即已传遍全身。因此，心脏实际上并非像目前教科书所述那样，简单地理解成把血液以流速v泵向全身的泵，而是脉搏波发生器[6]。

图1 波速C 取决于本构方程p-V 曲线的局部斜率（dp/dV）Fig. 1 The wave velocity C depends on the local slope(dp/dV) of the p-V curve

(3)式(1)～(4)共同组成脉搏波系统的控制方程组，其中的3 个守恒方程是普适的，反映了各种脉搏波的共性方面；而非普适的系统本构方程则反映了不同脉搏波的不同特性。可见，不同脉搏波系统的传播特性主要由脉搏波系统的本构关系决定，而各种病态可以看作偏离常态本构关系的不同反映。联想到近年来引人注目的中医体质学[7-8]，就体质的重要性和地位而言，相当于前述量化的脉搏波系统中的本构关系。

需要强调，脉搏波本构关系的研究涉及对脉搏波性质的理解。目前很多研究者（例如费兆馥[9]、Luo 等[10]）把脉搏波理解为脉搏压力波（纵波），也有研究者（例如Xue 等[11]）把脉搏波理解为沿固体血管传播的径向位移波（横波），而实际上脉搏波是流-固耦合和纵波-横波耦合的复杂波[5]，其波速依赖于2 个无量纲参数：血液-血管模量比Kb(p)/Ev(p)和血管径厚比D(p)/h0，前者是本构关系的反映。对于这样的脉搏波系统的本构关系，迄今未见有关研究报导。

既然脉搏波的传播特性取决于本构关系，问题的研究核心之一就在于如何通过实验研究来确定脉搏波本构关系，以及如何采用这些方法从现有文献数据来获得脉搏波本构关系，这正是本文中研究的目的。

怎么来确定脉搏波系统本构方程呢？目前可从以下3 个途径来进行探索：(1)由实测脉搏波波速对压力的关系C(p)进行反分析；(2)对脉搏波p-V本构关系进行实测；(3)由一系列实测脉搏波波形进行Lagrange 反分析。下面分别加以探讨。

1 由实测脉搏波波速与脉压的关系 C (p) 进行反分析

近年来，有关脉搏波传播速度C与血压p的关联性C(p)之研究，在西医界倍受重视。

西医界把脉搏波波速PWV（pulse wave velocity，即本文式(5)所示的C）看作反映动脉弹性和硬化程度的指标。2007 年欧洲高血压指南将主动脉PWV 高于12 m/s 列为高血压的危险因素之一[12]。《中国高血压防治指南（2018 年修订版）》也将PWV 列为评估高血压患者心血管风险的重要指标[13]。进一步，人们把实测的C(p)作为无创伤、无袖带、连续监测血压的重要途径之一，加以更广泛研究。通常从Moens-Korteweg (MK)公式[14-15]和 Hughes 公式[16-17]出发，建立脉搏波波速C与血压p的关系。

MK 公式[14-15]为：

式中：E、h0和D0(=2R0)分别为血管弹性模量、初始壁厚和直径，ρb为血液密度，E0为血压为零时的血管弹性模量，ζ 为血管的材料系数。

Ma 等[18]指出，MK 公式[14-15]包含2 个基本假设：(1)血管是弹性薄壁管（h0/D0＜10）；(2)随血压的变化，动脉的厚度和半径保持不变。然而，对于人类动脉，这2 个假设并不成立；而Hughes 公式[16-17]则完全是经验的，缺乏任何理论基础。Ma 等[18]基于Fung 的超弹性模型[19]，建立了一个由以下无量纲公式表征的、无需上述假设和经验表达式的分析模型（以下简称RH 模型）：式中：α 和β 为依赖于材料和血管几何参数的常数。

其实，按照我们的整体观脉搏波模型给出的波速C表达式（式(5))，可以直接由C(p)反演系统本构关系。事实上，式(5)可改写为：

可见，RH 简化式实际上对应于指数形式的p(V)关系或对数形式的V(p)关系。具体的p-V曲线特性取决于本构参数α 和β。能不能利用文献数据获得具体的脉搏波本构关系？下面讨论2 个实例。

(1)实例1

Ma 等[18]指出，在人类血压范围（5～20 kPa）内，若取：

RH 简化式（式(7b)）能与RH 模型（式(7a)）高度一致，如图2(a)所示。

图2 RH 简化式（式(7b)）与RH 模型（式(7a)）和文献实测舒张压数据的比较Fig. 2 Comparison of the RH simplified formula (Eq.(7b)) with the RH model (Eq. (7a))and the measured diastolic blood pressure (DBP) literature data

把式(10a)代入式(9a)，取p0=5 kPa，V0=0.995×10-3m3/kg，或即ρ0=1.005×103kg/m3，可以得出对应的p-V关系式为：

式中：p的单位为kPa，V的单位为10-3m3/kg。相应的曲线如图3(a)所示。

(2)实例2

Ma 等[18]还指出，若取：

RH 简化式(7b)能与文献报道的人体舒张压实测数据较好吻合，如图2(b)所示。

把式(11a)代入式(9a)，取p0=6 kPa，V0=0.995×10-3m3/kg (ρ0=1.005×103kg/m3)，可以得出对应的p-V关系式为：

式中：p的单位为kPa，V的单位为10-3m3/kg。相应的曲线如图3(b)所示。

图3 不同情况下RH 简化式（式(7b)对应的p-V 曲线Fig. 3 p-V curves for the RH simplified formula (Eq. (7b)) under different conditions

对比图3 中两曲线后可见，图3(b)曲线陡峭得多，意味着波速C随脉压p升高而迅速增大。由此可见，不同的本构参数α 和β 对于本构关系p-V曲线有明显影响，必然会影响脉搏波波速C以及不同形式脉搏波的波形传播特征，从而影响脉搏波能量和信息的传递。这为我们下一步定量化研究中医脉诊的脉象特征与系统本构关系p-V曲线特征之间的联系，包括研究常态及病态脉象与系统本构关系p-V曲线之间的定量联系等等，提供了一个良好开端和基础。

2 对脉搏波 p-V 关系进行实测

对脉搏波p-V本构关系最直接的实验研究方法无疑是对其进行实测。为说明这一点，先回顾一下式(6)所示的MK 公式和Hughes 公式。这一在西医界获得广泛认可和采用的公式，从我们的整体观脉搏波模型的观点来重新审视的话，MK 公式只是本模型的式(5)在“血管为弹性薄壁管”特定假设下的简化特例，而Hughes 公式实质是对狗进行体内（in vivo）+体外（in vitro）的脉搏p-V关系实测后所得出的。

事实上，由弹性力学[20]知内压p下薄壁圆管的环应力σθ和环应变εθ分别为：

至于Hughes 公式，回顾其原文[16-17]可知，Hughes 等在狗的活体体内实验中，用三轴超声导管测量了主动脉内外径随血压的变化，如图4(a)所示。由式(12)的分析可知，由此实际上已得出比容随压力的变化V(p)，即相应的活体p-V本构关系。而Hughes 等在死狗体外实验中，则直接对p-V本构曲线进行了测量，如图4(b)所示。

图4(b) Hughes 等[16-17]在死狗体外直接测得的p-V 曲线Fig. 4(b) The p-V curve directly measured in the in-vitro measurements for dogs by Hughes et al[16-17]

图4(a) Hughes 等[16-17]在狗的活体体内实验中，用三轴超声导管测得的主动脉内外径随血压的变化Fig. 4(a) Aortic inner and outer diameters along with aortic pressure recorded simultaneously and continuously with a triaxial ultrasonic catheter in the in-vivo measurements for dogs by Hughes et al[16-17]

Hughes 等[16-17]的实验研究结果当时以E(p)的形式（式6(b)）表述，而没有以更具实质意义的p-V本构关系的形式表述，显然是受到当时已被广泛认可的MK 公式的影响。其实，由Hughes公式（式(6b)）和式(12c)可导出：

由此可见，MK-Hughes 公式（式(6)）实际上对应于式(14)所示的指数函数形式的V(p)，或对数函数形式的p(V)本构关系，式中的η0和η1依赖于p0、V0、D0/h0、E0和本构参数ζ。

Hughes 等的贡献常常被归结于提出了Hughes 公式（式(6b)）。其实，更有价值的贡献在于对脉搏波p-V本构关系进行体内+体外的直接实测。Hughes 等是开创这方面研究的先驱，虽然尚限于狗的实验。

把上述脉搏波p-V关系实测法与C(p)关系反分析法两者相比较，前者属于直接实测法，其误差一般可小于后者的间接反演法；但前者属于有创法，而后者属于无创法，因而更便于推广应用。不过，不论是基于RH 简化式反演的p-V本构关系（式(9)和图3），还是基于Hughes 等实测结果建立的p-V本构关系（式(14)），都只是提供了一个基于大量实测数据的统计平均结果。从中医脉诊学的观点来看，脉象因人、因时、因地而异，换句话说，脉搏波的本构关系必定因人、因时、因地而异；甚至于同一个人因测量部位、测量时辰、病况变化而异。原则上，如同中医体质分类那样，本构关系也有不同类型，需要因人而定。

如何能够针对个人、无创地由实测数据来确定其脉搏波p-V关系？本文中建议采用脉搏波反分析法，即通过实测一系列波形进行Lagrange 反分析[2-3,21]，下面作进一步具体讨论。

3 由系列实测波形进行Lagrange 反分析

由连续介质力学波动理论[1-3]，Lagrange 反分析基于动量守恒方程（式(1)）和质量守恒方程（式(2)）。前者建立了压力p(X,t)的偏导数和质点速度v(X,t)的偏导数之间的关系，而后者建立了比容V(X,t)的偏导数和质点速度v(X,t)的偏导数之间的关系。通过求解v(X,t)作为过渡，可以求取脉搏波的p-V本构方程（式(4)）。具体步骤如下：

(1)在设定的测点设置由n(n≥3)个相邻的测压元件组成的阵列式压力传感器。一旦实测得到p(Xi,t)，i=1, 2 , ··· ,n，就可由式(1)求得 ∂v/∂t；应用初始条件v|t=0=v0，即可积分求得v(X,t)。

(2)求得v(X,t)后，就可由式(2)求得 ∂V/∂t；应用初始条件V|t=0=V0，即可积分求得V(X,t) 。

(3)由p(X,t)和V(X,t)消去t，即得p-V本构方程。

根据实测的一系列脉搏波信息去反求系统本构方程，力学上称为解“反问题”（反演）。其实，中医的切脉相当于根据脉诊信息对病况做出诊断，以便对症下药，予以诊治，以数理语言来表述，同样是在解“反问题”。两者在思维方式上是相通的。

为正确而成功地施行脉搏波的反分析，有必要强调以下2 点：(1)正确选取测点；(2)满足必需的测量分辨率和精度。

首先，测点取何处？原则上凡能测量到脉搏波的地方都可以，当然以脉搏波信号越强越好。联系到至今仍指导着中医临床实践的“寸口三部九候诊法”， 研究者们显然会更多关注桡动脉切脉部位“寸、关、尺”3 处。初步容易设想，能不能直接对寸、关、尺3 部的脉搏波进行Lagrange 反分析[2-3]。但实践发现并不可行，这一方面由于寸、关、尺3 部的间距约为10 mm，达不到反分析差分数值分析的精度要求；另一方面，也是更重要的，由于按照中医寸口脉分候脏腑学说，左手寸、关、尺分候心、肝、肾，右手寸、关、尺分候肺、脾、命门（右肾），则显然在寸、关、尺3 处分别确定的本构关系既有基本共性，又分别包含不同的分候特性，不能等同视之。事实上，脉诊实践中常常出现关部脉搏波最强，寸部次之，尺部最弱的情况。由于Lagrange 反分析适用于波强单调变化的简单波，这样的强弱起伏的脉搏波分布显然不符合这一要求。所以就寸口脉而言，应该对寸、关、尺分别测量，在每个测点采用阵列式传感器测取系列波形后进行反分析。这样才能更好结合中医脉诊特点，以促进传统中医脉诊的定量化发展。

其次，在进行Lagrange 反分析时都要经历数值微分和数值积分运算，从而存在一个如何满足对于测量分辨率和计算精度等要求的问题。下面通过一个实例来加以讨论。

Hu 等[22]采用图5(a) 所示的电容型阵列传感器测量了脉搏压力波的时空分布。阵列传感器长10 mm，宽7.5 mm，含4 个长度方向和3 个宽度方向共12 个微传感器。典型的“关”部的测量结果如图5(b)所示，图中纵坐标是无量纲压力波幅值，横坐标是采样数。把图5(b)的中间一列前3 行（即图中Ch2、Ch5、Ch8）的实测无量纲压力波形换算为p-t图，则如图5(c)所示。图中波形沿脉搏波传播方向依次递减，可满足Lagrange 反分析的要求。然而，实验中采用的扫描速率为100 Hz，相当于采样时间即时间分辨率为10 ms，其精度对于Lagrange 反分析是远远不够的。由于微传感器中心间距为2.5 mm（参看图5(a)），如果波速按照10 m/s 计，对应的时间差为0.25 ms。因此，传感器的时间测量精度至少应该是25 µs，才能保证两位有效数的测量精度。

图5(a) Hu 等[22]的实验中采用的电容型阵列传感器Fig. 5(a) The capacitor array sensor used in the experiment by Hu et al[22]

图5(b) Hu 等[22]的实验中“关”部各通道的实测波形Fig. 5(b) The waveform of each channel at Guan pulse taking position in the experiment by Hu et al[22]

图5(c) 对应于Hu 等[22]的实验中电容型阵列传感器通道2、5、8 的p-t 图的电信号Fig. 5(c) Electrical signals corresponding to the p-t waveforms measured at Ch2, Ch5 and Ch8 of the capacitor array sensor in the experiment by Hu et al[22]

Hu 等[22]的实验虽然提供了有益的信息，但因其扫描速率不够高，尚难以据此进行Lagrange 反分析，由此足以说明测量精度的重要性。今后需要研制具有微秒级别时间分辨率的毫米级阵列传感器，以推动Lagrange 反分析法在脉搏波本构关系研究中的应用。

这是目前尚无由脉搏波Lagrange 反分析来确定本构关系之实例的主要原因，也是我们正在进行中的研究工作，将另文讨论。

4 结 论

脉搏波控制方程组的本构关系决定着脉搏波的传播特征，但迄今缺乏这方面的有关研究报导。如何通过实验研究来确定脉搏波本构关系，以及如何通过这些方法从现有文献数据来获得脉搏波本构关系，是当前研究的核心之一。

本文中探索了3 个可行途径：(1)由实测脉搏波波速C对压力p的关系C(p)进行反分析（无创法)；(2)直接对脉搏波p-V关系进行实测（有创法）；(3)由一系列实测脉搏波波形进行Lagrange 反分析（无创法）。

采用上述方法，根据现有文献数据，成功地获得了不同类型的脉搏波本构关系。由C(p) 关系的RH 简化式实际上可推得指数型p(V)本构关系。由MK-Hughes 式实际上可推得对数型p(V)本构关系。脉搏波传播特性随非线性本构参数发生显著变化。这一研究还为今后从大量已有文献数据探索其背后更具实质意义的本构关系提供了一个新思路、新途径。

按照中医体质分类观点，脉搏波本构关系必定也有不同类型，具体因人、因时、因地而异；并且对个人还因测量部位、测量时辰、病况变化而异，原则上不存在单一普适的关系。在这个意义上，脉搏波的Lagrange 反分析具有广阔发展前景，但它对正确选择测点和提高测量敏感度和精度等方面提出了更高要求。