引　言

脉搏波是反映人体心血管健康状况的重要生理参数之一。目前，脉搏波测量方式主要分为侵入式、接触式和非接触式。侵入式测量主要用于临床手术，对人体有一定的损害。接触式测量需要紧贴患者皮肤，有时还需要压迫血管来获取生理参数，有一定的接触不适感。为了克服侵入式和接触式测量的不足，近年来，非接触式测量脉搏波的方法被广泛研究。

成像式光电容积描记法（imaging photo- plethysmography，IPPG）^[1]是应用广泛的非接触式测量方法之一。该方法的基本原理为：皮下浅层血管的血液流动会引起皮肤表面的颜色变化和微弱振动^[2-3]，通过获取和分析这类信号可以得到人体的某些生理参数，如心率^[4]、呼吸频率^[5]、血氧^[6]、心率变异性^[7]等。基于颜色变化获取脉搏生理参数的研究有：李江山等^[8]基于欧拉影像放大方法获取脉搏波信号，但获取的脉搏波信号的特征点（潮波和重搏波等）较少，只能简单地实现心率非接触式自动测量；Wang等^[9]通过自适应奇异谱分析法获取了保留部分特征点的脉搏波信号。通过皮肤颜色变化获取脉搏波时，由于皮肤表面的颜色变化极易受到光照、运动等因素影响，极难提取出脉搏波信号。基于振动获取脉搏生理参数的研究有：Wu等^[10]通过光学三角测量技术测量了桡动脉的皮肤振动，获取到具有一定周期性的脉搏信号；Lin等^[11]通过双目相机拍摄桡动脉表面振动，通过多点振动信号获取了与接触式方法较一致的脉搏波；Chen等^[12]在近红外光下通过颈动脉的皮肤振动获取了脉搏信号，但仅有波形上的趋势。通过皮肤振动提取的脉搏波信号，虽然周期相对稳定且有趋势性，但缺少脉搏波的特征点。这主要是由于皮肤表面振动微弱，所获取的原始脉搏信号中的特征点信号弱，被湮没在噪声中极难提取。

基于桡动脉搏动会引起表面皮肤的振动，本研究通过莫尔条纹技术放大桡动脉搏动的原始信号，使用工业相机连续拍摄在LED光源照射下的桡动脉处皮肤振动引起的莫尔条纹周期性变化的图像序列。选取条纹移动变化区域作为感兴趣区域（region of interest，ROI），对图像序列的ROI取均值来获得原始的脉搏波信号，再对其进行预处理和变分模态分解处理以获得保留多数潮波和重搏波等多个特征的完整脉搏波。

1 信号采集及原理

脉搏波是心脏的搏动沿着动脉血管和血流向外围传播而形成的波动信号。当血液流经人体皮下浅层动脉时，脉搏波会引起表层皮肤的微小振动。在临床上，桡动脉是获取脉搏波的常用部位，通过相机的快速连续拍照可以记录下桡动脉皮肤振动的变化过程，如图1所示。

脉搏波具有周期性，所以桡动脉表面的皮肤振动也具有周期性。在信号采集时，当光照在桡动脉处，大部分的光在皮肤表面发生漫反射和散射。桡动脉处皮肤的周期性振动使反射光的强度也随之发生周期性的变化，在相机记录的桡动脉搏动图像上便表现为某块区域的灰度值也具有周期性的变化。通过处理这些灰度变化的图像序列，可以获得一个周期性变化的信号，该信号能反馈脉搏波相关信息^[13]。

1.1 莫尔条纹

由于桡动脉表面皮肤振动微小，直接从连续拍摄的桡动脉搏动图像序列中很难提取脉搏波信号或特征点丢失严重。在实验中，于桡动脉表面皮肤处粘贴一片等线宽条纹光栅薄膜片，使桡动脉搏动位置处于薄膜片中心，在距离桡动脉搏动位置2～3 cm处放置一块透明玻璃片，玻璃片上粘贴相同的条纹薄膜片。调节两组条纹间的角度，可得到有较大间距的随桡动脉搏动而周期平移变化的莫尔条纹，以此来放大桡动脉处的微弱振动信号，如图2（a）所示。

图2（b）所示是两块相同的条纹间距为W的光栅以θ角偏转叠在一起产生的莫尔条纹现象，所形成的莫尔条纹宽度为B。假设图2（c）中a、b、c为垂直光栅中相邻的3根线条，而e、f为偏转θ角光栅的某相邻的2根线条，则

$ B = \dfrac{W}{{\sin \dfrac{\theta }{2}}} $

(1)

由式（1）可知，当保持两块光栅偏转角度不变时，垂直光栅水平移动ΔW，另一块光栅保持不动，则产生的莫尔条纹将移动ΔB。在运动关系上，光栅的移动距离和莫尔条纹移动的距离呈线性相关，即可通过莫尔条纹的移动量来放大光栅的位移量，放大倍数D为

$ D = \dfrac{{\Delta B}}{{\Delta W}} = \dfrac{1}{{\sin \dfrac{\theta }{2}}} $

(2)

可见，θ角越小，产生的莫尔条纹宽度越宽，同时放大倍数也越大。

1.2 实验采集装置

图3为脉搏波信号的采集装置示意图。普通LED光源从距离受试者手腕20～30 cm处，照射桡动脉搏动区域，相机位于手腕右侧60～70 cm处采集桡动脉搏动的图像序列，视场范围大约为5 cm×5 cm。实验所使用的相机为日本Hamamatsu公司生产的工业相机（型号：C11440–36U）。在采集脉搏信号时，受试者的另一手腕同时使用一款已在国内多家医院投入临床使用的接触式动脉血压及血流动力学监测设备（浙江善时医疗器械有限公司，型号：TL–400，注册证号：豫械注准20182210225），同步测量桡动脉脉搏波及心率、血压、脉压变异度等生理参数，用于与本文方法的实验数据进行比较。

2 信号处理及分析

首先使用相机采集桡动脉搏动产生的莫尔条纹图像序列。信号的采样频率设置为40 Hz，采样时间为15 s。

2.1 选取ROI获得原始信号

通过图像灰度值的阈值分割及掩模处理，可以得到近似椭圆的光照区域。对该区域进行边缘检测，选取图像上有莫尔条纹的区域作为ROI。

对图像序列中的每一帧图像求ROI的平均灰度值，将其作为该时刻桡动脉脉搏搏动的原始信号值。即在t时刻，脉搏信号s(t)可表示为

$ s(t) = \dfrac{1}{{nm}}\displaystyle\sum\nolimits_{x = 0}^n {\displaystyle\sum\nolimits_{y = 0}^m {G(x,y)} } $

(3)

式中：$ G(x，y) $表示ROI在$ (x，y) $处的灰度值；n，m表示ROI的宽度和高度。获取的原始脉搏波信号如图4所示。

2.2 预处理

ROI会随着受试者的手腕不自主的抖动等动作产生一定程度上的偏移，从而出现误差较大的离群值。为了提高脉搏波提取的准确度，首先将采集的信号$ S=s(t)（t=1，2，3，\cdots，600） $进行归一化处理，然后将其与左右相邻两点差值大于阈值（阈值设置为0.5，正常相连的采样点间的差值不超过0.4）的异常值剔除，并用左右相邻两点的均值代替。异常值处理过程为

$ s\left( t \right) = \left\{ {\begin{array}{*{20}{c}} {s\left( {t + 1} \right)}&{t = 1} \\ {\left[ {s\left( {t - 1} \right) + s\left( {t + 1} \right)} \right]/2}&{1 < t < 600} \\ {s\left( {t - 1} \right)}&{t = 600} \end{array}} \right. $

(4)

采样频率为40 Hz。由采样定理可知，通过采样点重建的信号频率为0～20 Hz。

通常人的脉搏频率为0.6～4 Hz，一般不会超过5 Hz^[14]。本文选择0.6～5 Hz的带通滤波对数据进行滤波处理，去除趋势和高频噪声，再归一化处理。预处理前后信号的频谱图如图5（a）和（b）所示，经预处理后的脉搏波如图5（c）所示。

在图5（b）中，幅值最大的频率对应的就是心脏搏动的频率，即心率（heart rate）R_h

$ {R_{\text{h}}} = \max (freq)60 $

(5)

式中，$ \max (freq) $表示在0.6～5 Hz信号频谱中最大幅值对应的频率值。

图5（c）中脉搏波特征点不够明显。在预处理后，尽管经带通滤波后的脉搏波信号是该频段内的主要信号，但混杂在该频段内的其余噪声并没有被去除，还需要对其进一步分解，做去噪处理。

2.3 变分模态分解

变分模态分解（variational mode decomposi- tion，VMD）^[15]是一种自适应、完全非递归的模态变分和信号处理方法，适用于非线性时间序列信号。该方法主要是用求解变分问题的思想，对信号进行提取，在不丢失原始信号特征的情况下，把一个原始信号分解成多个不同中心频率的信号。其原理是将非平稳信号f分解为k个模态分量子信号$ {u_k} $，且每个分量有一个确定的有限带宽和中心频率$ {\omega _k} $，同时要求各个模态的估计带宽之和最小。约束条件为所有模态之和与原始信号相同。重构时，首先通过对各个模态函数分量$ {u_k}(t) $进行Hilbert变换，获得其解析信号；然后将各解析信号与$ \text{e}^{-\mathrm{j}\omega_kt} $进行混合，将相应频谱调制到预估基频带，计算出平移后信号梯度的L₂范数平方，用于估计带宽^[16]。相应表达式为

$ \left\{ \begin{gathered} \mathop {\min }\limits_{\left\{ {{u_k}} \right\},\left\{ {{\omega _k}} \right\}} \left\{ {\displaystyle\sum\limits_k {\left. {\left\| {{\partial _t}\left[ {\left( {{\delta(t)} + \dfrac{{\text{j}}}{{{\text{π }}t}}} \right){u_k}(t)} \right]{{\text{e}}^{ - {\text{j}}{\omega _k}t}}} \right.} \right\|_2^2} } \right\} \\ {\mathrm{s.t.}}\displaystyle\sum\limits_k {{u_k}(t)} = f(t) \\ \end{gathered} \right. $

(6)

式中：{u_k}={u₁,u₂,u₃,···,u_k}，{ω_k}={ω₁,ω₂,ω₃,···,ω_k}分别为所有模态及其中心频率；k为分解的模态个数；∂_t为对t求偏导；δ（t）为狄拉克分布。

为了求解上式，引入Lagrange乘法算子$ \lambda $，将其转换为非约束变分问题，得到增广Lagrange函数，表达式为

$\begin{split} L(\{ {u_k}\} ,\{ {\omega _k}\} ,\lambda ) = & \alpha \displaystyle\sum\limits_k {\left. {\left\| {{\partial _t}\left[ {\left( {{\delta}(t) + \dfrac{{\text{j}}}{{{\text{π }}t}}} \right){u_k}(t)} \right]{{\text{e}}^{ - {\text{j}}{\omega _k}t}}} \right.} \right\|_2^2} + \\ & \left\| {\left. {f(t) - \displaystyle\sum\limits_k {{u_k}(t)} } \right\|} \right._2^2 + \\ & \left\langle {\left. {\lambda (t),f(t) - \displaystyle\sum\limits_k {{u_k}(t)} } \right\rangle } \right. \end{split} $

(7)

式中，α为二次惩罚因子，其作用是降低高斯噪声干扰，保证信号的重构精度。利用交替方向上乘子算法迭代出$ \widehat u_k^{n + 1} $、$ \widehat \omega _k^{n + 1} $和$ {\widehat \lambda ^{n + 1}} $，求解最优解。求解的表达式为

$ \widehat u_k^{n + 1}(\omega ) = \dfrac{{\widehat f(\omega ) - \displaystyle\sum\limits_{l < k} {\widehat u_l^{n + 1}(\omega )} - \displaystyle\sum\limits_{l > k} {\widehat u_l^n(\omega )} + \dfrac{{{{\widehat \lambda }^n}(\omega )}}{2}}}{{1 + 2\alpha {{(\omega - \omega _k^n)}^2}}} $

(8)

$ \omega _k^{n + 1} = \dfrac{{\displaystyle\int_0^\infty {\omega {{\left| {\widehat u_k^{n + 1}(\omega )} \right|}^2}{\text{d}}\omega } }}{{\displaystyle\int_0^\infty {{{\left| {\widehat u_k^{n + 1}(\omega )} \right|}^2}{\text{d}}\omega } }} $

(9)

$ {\widehat \lambda ^{n + 1}}(\omega ) = {\widehat \lambda ^n}(\omega ) + \tau \left[ {\widehat f(\omega ) - \displaystyle\sum\limits_k {\widehat u_k^{n + 1}(\omega )} } \right] $

(10)

式中：$ \widehat f(\omega ) $，$ \widehat u_l^n(\omega ) $，$ {\widehat \lambda ^n}(\omega ) $和$ \widehat u_k^{n + 1}(\omega ) $分别为$ f(t) $、$ u_l^n(t) $、$ {\lambda ^n}(t) $、$ u_k^{n + 1}(t) $的傅里叶变换；$ \tau $为噪声容忍度；n为迭代次数。

VMD算法的具体过程如下：

1）设置最大模态分量k，n=0，初始化$ \widehat u_k^{n + 1} $，$ \omega _k^{n + 1} $，$ {\widehat \lambda ^{n + 1}} $为$ \widehat u_k^1 $，$ \omega _k^1 $，$ {\widehat \lambda ^1} $；

2） 根据式（3）～（5）迭代更新$ \widehat u_k^{n + 1} $，$ \omega _k^{n + 1} $，$ {\widehat \lambda ^{n + 1}} $；

3） 判断$ {\widehat u_k} $是否满足迭代停止条件，满足条件则停止迭代，得到最终模态分量$ {\widehat u_k} $及中心频率$ {\omega _k} $，迭代停止条件为

$ \dfrac{{\displaystyle\sum\limits_k {\left\| {\widehat u_k^{n + 1} - \widehat u_k^n} \right\|} _2^2}}{{\left\| {\widehat u_k^n} \right\|_2^2}} < {\varepsilon _1} $

(11)

式中，ε₁为预设的收敛误差。

对滤波预处理后的脉搏信号进行层数K值为5的VMD分解，分解后各频带的本征模态函数（intrinsic mode function，IMF）分量及频率如图6（a）和图6（b）所示。心跳频率及其谐波频率是脉搏波的主频率^[14]。IMF1～IMF5最大幅值对应的频率分别是0.80 Hz、1.13 Hz、2.20 Hz、3.33 Hz和4.53 Hz。对图6（b）分析可知，IMF1为低频噪声，IMF2为心跳频率，IMF3为二次谐波频率，IMF4为三次谐波频率，IMF5为高频噪声。所以IMF2、IMF3和IMF4为脉搏波主频率，将其进行重构可获得脉搏波信号。重构后的脉搏波信号如图6（c）所示。

2.4 脉搏波的特征点

一个完整周期的典型脉搏波一般包含主波、潮波和重搏波^[17-18]。如图7（a）所示，其波峰、波谷的特征点有明确的生理意义：A点表示脉搏周期开始；B点表示主波峰；C点表示潮波波谷；D点表示潮波波峰；E点表示降中峡；F点表示重搏波峰；G点即为下一个脉搏波的起点。这些特征点可用于血压、血氧等生理参数的预测。特征点的幅值大小能反映出脉搏的强弱，这对中医的诊断具有重要意义。T表示相邻脉搏波主波峰间的时间差，可用于计算心率变异性^[19]，计算式为

$ {V_{{\text{hr}}}} = std(T) $

(12)

式中：V_hr表示心率变异性（heart rate variablity，HRV）；$ std(T) $表示采集时间内每个相邻周期T的标准差。

本文方法所获取的脉搏波及提取的特征点如图7（b）所示。

3 实验结果及评价

本文设计了两个实验：实验一，邀请20名受试者（年龄20～65岁，其中男性13名，女性7名），用本文方法采集其桡动脉脉搏波信号，以验证本文方法的有效性；实验二，用血压及血流动力学监测设备TL–400同步测量脉搏波波形及心率等生理参数，并与本文方法的测量结果进行对比，以验证本文方法的可靠性和准确性。

3.1 实验一

采集受试者的脉搏波时，将照明光源调节到合适角度，受试者保持身体放松，呼吸平稳。随着脉搏搏动，粘贴在一侧桡动脉皮肤处的条纹薄膜片与玻璃片上的条纹薄膜片产生莫尔条纹。同时，在受试者另一手腕使用TL–400测量桡动脉脉搏波及其他生理参数。在一个测量周期内，实验测得的脉搏波特征点个数及计算的心率、心率变异性和特征周期占比如表1所示，对应获取的脉搏波波形如图8所示。

以实验组1为例，在采样时间内受试者的心率为68，即每分钟心脏搏动68次。本实验单次测量时长为15 s，理论上在一个测量时长内应有17个脉搏周期。从实验组1测得的脉搏波中提取特征点：A点出现17次，B点出现17次，C点出现17次，D点出现17次，E点出现13次，F点出现13次。R_AB、R_CD、R_EF分别表示主波、潮波和重搏波出现的周期次数占测量时间内的脉搏周期次数的百分比。实验组1的R_AB、R_CD、R_EF分别为100%，100%，76%。特征点周期次数占比越高，表示获取的脉搏波中的特征点信号越明显，脉搏波中特征点信号的保留程度越完整。

由表1和图8可知，在本实验20名受试者的脉搏波信号中，实验组1～18的脉搏波周期稳定，主波周期占比为100%，即在测量时间内可以完全测得每一个脉搏周期的主波信号。获取的这18组脉搏波中的特征点信号十分明显。潮波周期占比平均为88.1%，重搏波周期占比平均为78.9%，即潮波和重搏波的保留完整程度分别为88.1%和78.9%，可以认为绝大多数的细节特征得到保留。实验组19的潮波周期占比为18%，重搏波周期占比为65%，可判定为潮波缺失，而重搏波表达较为明显。初步认为，这是由于信号采集时光源强度过强，导致采集的图像产生过曝光，使得脉搏波特征缺失。实验组20所获取的脉搏波波形杂糅紊乱，仅测得几个完整的脉搏周期，无法计算出准确的心率变异性。波形杂糅紊乱的主要原因是受试者的腕部体脂过高，导致其桡动脉处皮肤表面振动十分微弱。

3.2 实验二

为验证本文方法的准确性，将获取到的波形、心率等与TL–400测得的数据作对比，并用Bland–Altman法^[20]进行一致性评估。其原理为：用两种方法对同一批受试对象同时进行测量时，一般不会获得完全相同结果，会存在一定趋势的差异。这种差异被称为偏倚，可以用两种方法的测量结果之差的平均数d进行估算。平均数d的变异情况可用差值的标准差S_d来描述。如果差值的分布服从正态分布，则 95% 的差值应该落在 [d−1.96S_d，d+1.96S_d]。该区间称为 95%一致性界限（95% limits of agreement，95% LoA）。当绝大多数差值都位于该区间内，从统计学意义上就可认为这两种方法具有较好的一致性，可以互相代替。

Bland–Altman法对心率的分析结果如图9所示。本文方法与TL–400设备测得的心率误差均值为0.1，置信区间为[−1.41, 1.61]，两者具有一致性。由此可以认为，从采集到的波形数据中得到了真实的心率参数。

为了进一步验证本文方法的准确性，分别取2名受试者通过两种方法获得的脉搏波作比较，详见图10。两种方法获得的脉搏波波形整体特征一致，表明用本文方法提取到的脉搏波与用接触式测量设备得到的脉搏波具有相同的特征点信号。

4 结　论

人体桡动脉的搏动会引起表面皮肤的振动。利用在桡动脉皮肤处粘贴的条纹光栅薄膜片与外置的玻璃片上相同的条纹薄膜片间产生的莫尔条纹，可有效放大桡动脉皮肤表面振动的原始信号。使用相机连续帧照相方式，在LED光源照射下，拍摄由桡动脉处皮肤振动产生的莫尔条纹图像序列，然后使用带通滤波和变分模态分解与重构算法获取人体脉搏波信号。测试结果及对比分析表明，本文方法能获得周期稳定，并保留大多数潮波和重搏波等特征的完整脉搏波信号。本文方法可用于非接触、无感状态下的连续脉搏波监测及与脉搏波相关的生理参数的测量。