2. 兰州理工大学 电气工程与信息工程学院, 甘肃 兰州 730050;
3. 兰州理工大学 材料科学与工程学院, 甘肃 兰州 730050
2. College of Electrical and Information Engineering, Lanzhou University of Technology, Lanzhou 730050, China;
3. College of Materials Science and Engineering, Lanzhou University of Technology, Lanzhou 730050, China
脉诊是传统中医基本诊断方法之一。医者通过手指感触病人桡动脉脉搏来辨识脉象,评价其健康状况。然而,脉诊过于依赖医生的主观经验,限制了中医的传承与发展。借助现代传感器与检测技术提升脉搏检测的客观性和准确性成为亟待解决的问题。
常见的脉搏检测装置按照检测原理可分为五类:第一类是直接测量脉搏搏动力变化的压电传感器[1-2];第二类是袖带式检测装置[3-4],将袖带紧附在手腕上,检测袖带内压变化获得脉搏波形;第三类是光电脉搏传感器[5-6],根据血管透光率检测血管容积变化,以此描述脉搏波形;第四类是基于超声多普勒技术的检测装置[7];第五类是利用声学原理拾取脉搏引发的微弱振动的所谓传声器[8]。压电传感器、光电传感器、袖带式检测装置和传声器仅能获取时域信号。超声多普勒技术可探测脉管时空域运动,但技术复杂成本较高。Yoo等[9]将6个压电传感单元封装在一起,采集脉宽方向上六路脉搏信号。Chu等[10]模仿三部九侯取脉方式研制了三探头检测装置。Peng等[11]用光刻技术在柔性基板上雕刻5×5个感应电极,制成传感器阵列。这些传感器阵列可获取时空域脉搏信号,但受传感单元密度限制,其空间分辨率有不高。
本文在前期工作[12]的基础上研制了一种具有气囊式柔性探头、切脉压力可调节、采用双目视觉测量的时空域脉搏信号检测系统。研究了三维脉图构建方法,并且初步探讨了脉象要素的量化指标。
1 检测系统 1.1 总体结构检测系统主要由脉搏图像采集部分、检测参数采集与调节部分构成。图 1是检测系统总体结构框图。相机选用Basler acA1300-30gm型高速工业黑白相机,配备Computar M0814-MP定焦镜头。北京双诺公司生产的MP425数据采集模块发出同步信号,触发两部相机同时采集接触膜图像[13],图像数据通过千兆以太网卡进行高速转存。图 2是检测装置实物图。
1.2 气囊式柔性探头
气囊式柔性探头是检测装置的关键部件,图 3(a)和(b)分别显示了探头及其拆解状态。探头由缸体、环形盖、透光片、密封圈构成。丁腈橡胶制成的接触膜固定于探头底部,膜厚0.5 mm。探头内压增大促使接触膜膨胀,其外部触感和力学性质与手指指腹近似。接触膜上印有网格状结构线,见图 3(c)。结构线交点将作为双目立体视觉三维测量的标识点。
2 三维脉图构建方法三维脉图是脉搏信号空间量化形式,反映脉搏强度的指面分布状态。脉搏强度越大接触膜变形越大,利用接触膜形变描述脉搏强度是构建三维脉图的新方式。图 4是利用双目立体视觉测量接触膜形貌、构建三维脉图的流程。采用张正友标定法获得两部相机的焦距和空间位置关系[14]。一般将两部相机分别称为左相机和右相机。左右相机的位置关系用平移矩阵T和旋转矩阵R表示,左右相机焦距分别表示为fl和fr。通过图像处理技术获得网格交点的图像坐标,进而根据双目立体视觉原理得到交点的三维坐标,最终重构生成三维脉图。
2.1 网格图像交点检测微距拍摄的网格线有宽度,其交点并非具体的点,而是纵横两段网格相交区域的中心位置,属于广义交点,一般的特征点检测方法难以奏效。提出一种基于脊线的交点检测新方法。首先将原灰度图转化为二值图,提取网格图像骨架线,基于骨架线估计出网格结构参数,根据结构参数将图像划分成若干检测区域,每个区域仅包含一个待检交点。在每个区域内提取网格线的脊线,对脊线进行拟合,拟合线交点即为网格线的广义交点。
骨架线是一种与原图像具有相同连通性和拓扑结构的线性几何体,较常使用的中轴骨架线模型定义为
$ {d_8}\left( {q,\mathit{\boldsymbol{B}}} \right) < {d_8}\left( {p,\mathit{\boldsymbol{B}}} \right) + {d_8}\left( {p,q} \right) $ | (1) |
式中:d8(·) 表示八邻域棋盘距离,p(x, y) 和q(s, t) 为像素点,B表示像素点的集合。像素点之间的棋盘距离表示为
$ {d_8}\left( {p,q} \right) = \max \left( {\left| {x - s} \right|,\left| {y - t} \right|} \right) $ | (2) |
像素点与点集之间的棋盘距离表示为
$ {d_8}\left( {p,\mathit{\boldsymbol{B}}} \right) = \inf \left\{ {{d_8}\left( {p,z} \right)\left| {z \in \mathit{\boldsymbol{B}}} \right.} \right\} $ | (3) |
式中:inf表示上确界。满足式 (1) 的点p的集合构成网格图像的骨架线。根据上述定义提取的网格图像股骨架线如图 5所示。
通过下式可提取骨架线交点:
$ \mathit{\boldsymbol{I' = I}} * \mathit{\boldsymbol{k}} \cdot \mathit{\boldsymbol{I}} $ | (4) |
式中:I表示骨架的二值模式,k为3×3全1掩膜,I′中值大于3的单元组成的联通区域质心即为骨架交点。
检测区域中心可由骨架线交点确定,检测区域大小则要根据网格线宽度、方向和网眼大小确定。利用Canny算子进行边缘检测,再用Hough变换检测边缘上的直线段,见图 6。计算所有同向直线段之间的距离产生集合D,其中元素dij是第i和第j条直线段之间的距离。将D的所有元素排序后绘制于坐标系中,如图 7所示。分析发现,这些点明显聚合为若干类,其中第一类点表示的距离为网格线宽度。第二类点表示网眼大小。利用层次聚类法分别取得第一类和第二类点集,求其平均值即为网格线宽度和网眼大小估计值。
将原始灰度图像划分成若干检测区域,每个检测区域包含两段“+”形网格线。提取两段网格线的脊线,脊线交点就是要找的网格线交点。脊线是一种重要的图形图像特征,可反映形体走势。图 8显示了检测区域图像的三维形式,其脊线特征非常明显。对于灰度图像,脊线具有如下定义:
$ \frac{{\partial {k_{\max }}}}{{\partial {t_{\max }}}} = 0,\frac{{{\partial ^2}{k_{\max }}}}{{\partial t_{\max }^2}} < 0,{k_{\max }} > \left| {{k_{\min }}} \right| $ | (5) |
式中:kmax和kmin为曲面的主曲率函数,tmax和tmin为曲面的主方向函数。主曲率函数一阶偏导为零规定脊点处曲率存在极值,主曲率函数二阶偏导小于零规定曲面在脊点处应是凸的。根据定义提取脊线计算量较大,且结果中含有图像噪声引起的干扰脊线。根据网格线宽度方向上灰度近似高斯分布这一特点,提出灰度最值扫描法,设置垂直于网格线的扫描线沿网格线方向行进,提取扫面线覆盖区域灰度最值作为脊点,连接脊点得到主脊线,并对主脊线拟合进一步降低噪声影响,最终得到交点,如图 8所示。
需要指出,骨架线交点虽也在检测区域内,但它基于二值图像提取,精确度较低。以多人次手工标定均值作为参照,骨架线交点平均偏差4.49像素,而脊线交点平均偏差2.13像素,因此脊线交点更为理想。
2.2 对应交点匹配利用图像视差小、交点接近阵列分布等特点,将左右图像当中的交点集放到同一平面坐标系当中,按照欧氏距离最近原则进行匹配,匹配结果如图 9所示。
2.3 交点三维坐标计算图 10是双目立体视觉三维测量原理图。图中,Oxyz为左相机坐标系,也是世界坐标系,Orxryrzr为右相机坐标系,左右图像坐标系分别为OlXlYl和OrXrYr,空间点P(x, y, z) 在左右图像平面上的投影分别用 (Xl, Yl) 和 (Xr, Yr) 表示。
根据小孔成像原理和透视投影模型,可以得出
$ {s_1}\left[ {\begin{array}{*{20}{c}} {{X_1}}\\ {{Y_1}}\\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{f_1}} & 0 & 0 & 0\\ 0 & {{f_1}} & 0 & 0\\ 0 & 0 & 1 & 0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x\\ y\\ z\\ 1 \end{array}} \right] $ | (6) |
$ {s_{\text{r}}}\left[ {\begin{array}{*{20}{c}} {{X_{\text{r}}}} \\ {{Y_{\text{r}}}} \\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{f_{\text{r}}}}&0&0&0 \\ 0&{{f_{\text{r}}}}&0&0 \\ 0&0&1&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x_{\text{r}}}} \\ {{y_{\text{r}}}} \\ {{z_{\text{r}}}} \\ 1 \end{array}} \right] $ | (7) |
式中:sl和sr为比例常数。设左右相机之间的位置关系为
$ {\mathit{\boldsymbol{M}}_{1{\rm{r}}}} = \left[ {\mathit{\boldsymbol{R}}\left| \mathit{\boldsymbol{T}} \right.} \right] = \left[ {\begin{array}{*{20}{c}} {{r_1}} & {{r_2}} & {{r_3}} & {{t_x}}\\ {{r_4}} & {{r_5}} & {{r_6}} & {{t_y}}\\ {{r_7}} & {{r_8}} & {{r_9}} & {{t_z}} \end{array}} \right] $ | (8) |
则
$ {\left[ {\begin{array}{*{20}{c}} {{x_{\rm{r}}}} & {{y_{\rm{r}}}} & {{z_{\rm{r}}}} \end{array}} \right]^\prime } = {\mathit{\boldsymbol{M}}_{1{\rm{r}}}}{\left[ {\begin{array}{*{20}{c}} {{x_1}} & {{y_1}} & {{z_1}} & 1 \end{array}} \right]^\prime } $ | (9) |
由式 (6)~(9) 可以得出左右图像坐标系的关系:
$ {\rho _{\rm{r}}}\left[ {\begin{array}{*{20}{c}} {{X_{\rm{r}}}}\\ {{Y_{\rm{r}}}}\\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{f_{\rm{r}}}{r_1}} & {{f_{\rm{r}}}{r_2}} & {{f_{\rm{r}}}{r_3}} & {{f_{\rm{r}}}{t_x}}\\ {{f_{\rm{r}}}{r_4}} & {{f_{\rm{r}}}{r_5}} & {{f_{\rm{r}}}{r_6}} & {{f_{\rm{r}}}{t_y}}\\ {{r_7}} & {{r_8}} & {{r_9}} & {{t_z}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {z{X_1}{f_1}}\\ {z{Y_1}{f_1}}\\ z\\ 1 \end{array}} \right] $ | (10) |
式中:ρr为比例因子,于是空间点三维坐标可以表示为
$ \left\{ \begin{array}{l} x = z{X_1}{f_1}\\ y = z{Y_1}/{f_1}\\ z = \frac{{{f_1}\left( {{f_{\rm{r}}}{t_x} - {X_{\rm{r}}}{t_z}} \right)}}{{{X_{\rm{r}}}\left( {{r_7}{X_1} + {r_8}{Y_1} + {f_1}{r_9}} \right) - {f_{\rm{r}}}\left( {{r_1}{X_1} + {r_2}{Y_1} + {f_1}{r_3}} \right)}} = \\ \;\;\;\;\;\frac{{{f_1}\left( {{f_{\rm{r}}}{t_y} - {X_{\rm{r}}}{t_z}} \right)}}{{{Y_{\rm{r}}}\left( {{r_7}{X_1} + {r_8}{Y_1} + {f_1}{r_9}} \right) - {f_{\rm{r}}}\left( {{r_4}{X_1} + {r_5}{Y_1} + {f_1}{r_6}} \right)}} \end{array} \right. $ | (11) |
利用式 (11) 可获得网格线交点三维空间坐标。进一步对这些空间点集进行三次多项式曲面拟合即得到三维脉图。
3 局部脉搏有限元仿真为了验证双目视觉测量接触膜形变的结果,建立在探头作用下局部脉搏的静力有限元模型,如图 11所示。探头几何尺寸与实际相同,接触膜厚度0.5 mm,血管内直径1.4 mm,外直径2.2 mm。探头为刚体材料,接触膜采用二阶Mooney-Rivlin超弹性材料模型,材料常数根据拉伸实验获得,血管和软组织材料模型采用Holzapfel等[15]的研究成果。仿真过程分为三步完成:1) 施加探头内压,血管内壁加载舒张压;2) 给探头施加向下的压力,使接触膜与软组织充分接触;3) 血管内壁压力提升至收缩压。
4 实验及分析利用时空域脉搏信号检测系统对健康在校大学生受试者进行检测。双相机图像分辨率设置为600×600,格式为24位JPEG灰度图。由于图像数据存储量较大,受数据同步存储速率限制,双相机每秒各采集15帧图像。利用上述三维脉图构建方法对采集的图像进行三维重构。图 12(a)是双目视觉测量获得的一幅典型三维脉图,图 12(b)为有限元仿真获得三维脉图。图 12(a)中曲面的形态与 (b) 中曲面的局部十分相似,且与文献[11]的测量结果十分近似,证明本文提出的检测方法是有效的。图 12(a)显示双目视觉测量视场偏离脉搏中心,原因是相机视场较小,受试者身体细微移动导致视场中心偏离相机视场中心。另外,接触膜覆盖了桡动脉两侧的桡骨茎突和桡侧腕屈肌腱等较硬组织,使接触膜形态不规则,增加了对准脉搏中心的难度。
利用激光位移传感器实测接触膜的振动幅度,以此验证双目视觉测量结果。激光位移传感器选择Keyence IL-S065型,测量范围20 mm,精度为5 μm。图 13(a)双目视觉测量的脉搏中心振幅曲线,图 13(b)是激光位移传感器测得脉搏中心振幅曲线。
由于激光位移传感器采样与双目视觉测量并不同步,因此图 13(b)与图 13(a)波形不一致。为了获得两条曲线的振幅,首先进行高频滤波,然后利用动态差分阈值法[16]区分每个脉搏周期,最终提取各周期振幅并进行统计。双目视觉测量的平均脉搏振幅353 μm,标准差32 μm,激光位移传感器测量的平均振幅379 μm,标准差67 μm,两者无显著统计学差异,表明双目立体视觉测量具有较高的精确度。
利用三维脉图可量化描述中医脉象指标。中医研究者通常从“位数形势”四个方面描述脉象,而费兆馥等在此基础上进一步提出了7种指感因素:1) 脉位,指脉位深浅,轻取脉感明显重取脉搏变弱,是脉位较浅的表现,轻取不明显重按脉感明显,则是脉位较深的表现;2) 频率,传统记法指一呼一吸之间脉搏跳动次数,现在也指每分钟脉搏跳动的次数;3) 形态,指时域脉搏信号曲线的形态;4) 大小,手指的有脉感区域在血管径向上的宽度;5) 长短,手指的有脉感区域在血管轴向上的长度;6) 强弱,指脉搏搏动力的大小;7) 虚实,大而无力谓之虚,大而有力谓之实,虚实实际上是长短、大小和强弱的综合指标。
三维脉图中部有丘状凸起,凸起区域呈椭圆形,见图 14。图中a和b分别表示椭圆的长轴和短轴。结合三维脉图,上述7种脉象因素的量化指标如表 1所示。
表 1中f表示从中心点振幅提取的频率,Ampmax表示中心点振幅的最大值,脉搏波形态可从中心点振动曲线当中提取。根据中医理论,Ampmax与切脉压力 (即探头与手腕的接触压力) 关系曲线反映脉象沉浮,因此曲线斜率k可作为脉象沉浮指标。虚实则用Ampmax/b表示,比值越大脉象越实,比值越小脉象越虚。
5 结论1) 本文研制了一种具有气囊式仿指柔性探头的时空域脉搏信号检测装置,采用双目立体视觉获取三维脉图。
2) 与本课题组前期研制的单目装置相比,提高了三维测量精确度。与常见的压力传感器阵列检测三维脉图相比,本文方法的空间分辨率更高,且柔性探头对检测干扰更小,适合长时检测。
3) 本文方法获取的三维脉图序列当中包含节律、脉长、脉宽、脉形和强度等丰富的脉搏信息。初步提出了七种脉象因素量化指标,为准确采集、量化描述脉搏提供了一条有效途径。
4) 随着高速成像设备微型化,在本文成果基础上研制三探头检测系统,可模拟三部九侯取脉方式,为揭示中医脉诊机理奠定基础。
[1] | CIACCIO E J, DRZEWIECKI G M. Tonometric arterial pulse sensor with noise cancellation[J]. IEEE transactions on biomedical engineering, 2008, 55(10): 2388–2396. DOI:10.1109/TBME.2008.925692 |
[2] | TYAN C C, LIU S H, CHEN J Y, et al. A novel noninvasive measurement technique for analyzing the pressure pulse waveform of the radial artery[J]. IEEE transactions on biomedical engineering, 2008, 55(1): 288–297. DOI:10.1109/TBME.2007.910681 |
[3] | TANAKA S, NOGAWA M, YAMAKOSHI T, et al. Accuracy assessment of a noninvasive device for monitoring beat-by-beat blood pressure in the radial artery using the volume-compensation method[J]. IEEE transactions on biomedical engineering, 2007, 54(10): 1892–1895. DOI:10.1109/TBME.2007.894833 |
[4] |
王学民, 杨成, 陆小佐, 等. 基于柔性阵列传感器的脉象检测系统的设计[J].
传感技术学报, 2012, 25(6): 733–737.
WANG Xuemin, YANG Cheng, LU Xiaozuo, et al. The design of multi-channel pulse detection system based on flexible array sensor[J]. Chinese journal of sensors and actuators, 2012, 25(6): 733–737. |
[5] | NAM D H, LEE W B, HONG Y S, et al. Measurement of spatial pulse wave velocity by using a clip-type pulsimeter equipped with a hall sensor and photoplethysmography[J]. Sensors, 2013, 13(4): 4714–4723. DOI:10.3390/s130404714 |
[6] | WANG Dimin, ZHANG D, LU Guangming. A novel multichannel wrist pulse system with different sensor arrays[J]. IEEE transactions on instrumentation and measurement, 2015, 64(7): 2020–2034. DOI:10.1109/TIM.2014.2357599 |
[7] | VAPPOU J, LUO Jianwen, OKAJIMA K, et al. Pulse wave ultrasound manometry (PWUM):measuring central blood pressure non-invasively[C]//2011 IEEE International Ultrasonics Symposium. Orlando, FL:IEEE, 2011:2122-2125. |
[8] |
王炳和, 王海燕, 职利琴, 等. 高血压患者脉博信号的采集与频谱分析[J].
声学技术, 1999, 18(3): 135–138.
WANG Binghe, WANG Haiyan, ZHI Liqin, et al. Acquisition and spectral analysis of pulse signals of hypertension patients[J]. Technical acoustics, 1999, 18(3): 135–138. |
[9] | YOO S K, SHIN K Y, LEE T B, et al. Development of a radial pulse tonometric (RPT) sensor with a temperature compensation mechanism[J]. Sensors, 2013, 13(1): 611–625. DOI:10.3390/s130100611 |
[10] | CHU Yuwen, LUO C H, CHUNG Y F, et al. Using an array sensor to determine differences in pulse diagnosis-three positions and nine indicators[J]. European journal of integrative medicine, 2014, 6(5): 516–523. DOI:10.1016/j.eujim.2014.04.003 |
[11] | PENG J Y, LU M S C. A flexible capacitive tactile sensor array with CMOS readout circuits for pulse diagnosis[J]. IEEE sensors journal, 2015, 15(2): 1170–1177. DOI:10.1109/JSEN.2014.2360777 |
[12] |
张爱华, 李金平, 林冬梅. 基于视觉测量的脉搏图像三维重构[J].
数据采集与处理, 2012, 27(5): 570–575.
ZHANG Aihua, LI Jinping, LIN Dongmei. 3D reconstruction of pulse image based on visual measurement[J]. Journal of data acquisition & processing, 2012, 27(5): 570–575. |
[13] |
林冬梅, 张爱华, 沈蓉, 等. 双目视觉脉搏测量系统中的相机同步采集方法[J].
吉林大学学报:工学版, 2015, 45(6): 1999–2006.
LIN Dongmei, ZHANG Aihua, SHEN Rong, et al. Dual-camera synchronous acquisition method for binocular vision pulse measurement system[J]. Journal of Jilin university:engineering and technology edition, 2015, 45(6): 1999–2006. |
[14] | ZHANG Z. A flexible new technique for camera calibration[J]. IEEE transactions on pattern analysis and machine intelligence, 2000, 22(11): 1330–1334. DOI:10.1109/34.888718 |
[15] | HOLZAPFEL G A, OGDEN R W. Biomechanics of soft tissue in cardiovascular systems[M]. Vienna: Springer, 2003. |
[16] |
张爱华, 王平, 丑永新. 基于动态差分阈值的脉搏信号峰值检测算法[J].
吉林大学学报:工学版, 2014, 44(3): 847–853.
ZHANG Aihua, WANG Ping, CHOU Yongxin. Peak detection of pulse signal based on dynamic difference threshold[J]. Journal of Jilin university:engineering and technology edition, 2014, 44(3): 847–853. |