基于近红外血管成像的浅静脉三维定位研究

赵会梦，李臣鸿

（中南民族大学生物医学工程学院膜离子通道与药物研发实验室，武汉 430070）

0 引言

随着新型冠状病毒肺炎疫情的暴发，医护人员所必需的厚重防护服无疑会导致静脉穿刺失败率增加，这可能会延误抢救患者的时机，甚至会危及患者的生命，目前没有特定的医疗辅助工具可以解决这一问题。因此，研究静脉图像的识别及其处理，定位血管的三维位置对于辅助静脉穿刺等医疗手段具有重要意义。

血液中血红蛋白对近红外光的吸收特性使得近红外光在血管相关的临床应用范畴中具有独特优势，国内外均积极探索近红外光在血管显影相关领域的应用[1]。中国科学院西安光学精密机械研究所以其精湛的近红外成像技术成功研制出国内首台高速投影式血管显像仪，成像速率达25帧/s[2]。该研究对血管进行了二维平面定位，但缺乏血管的深度信息[2]。拉脱维亚电子与计算机科学研究所的学者们构建了用于血管结构分析的近红外成像系统，该系统能够测试适合掌静脉成像的近红外波长，但是同样缺乏血管的深度信息[3]。

提供二维信息的静脉显像仪的发展到目前为止已经颇为成熟，能提供三维信息的设备却仍然处在各大高校及研究机构的研发之中，并没有足以碾压其他产品的优势产品。本研究希望为血管的三维定位研究提供新的方法和思路，进而为医护人员静脉穿刺提供参考，达到辅助医疗的目的。因此，计算血管的皮下深度是本研究将静脉信息从二维升级为三维的重点，通过探寻不同波长下静脉图像中静脉灰度值和静脉皮下深度之间的关系，建立模拟方程，以此实现提供血管的三维位置的目的。

1 近红外静脉图像采集系统的搭建

近红外静脉图像采集系统分为2个部分：第一部分为图像采集设备。使用无畸变低照度近红外定制摄像头，加入700～1 100 nm的宽动态滤光片，能够滤除可见光的干扰。第二部分为自行设计制作的多波长近红外光源，该光源需要满足的条件如下：（1）多波长近红外光源在同一块电路板上；（2）不同波长光源位置分布均匀，以免造成误差；（3）不同波长的光源自动切换，以免手动切换导致误差。基于以上设计考量，本研究中多波长近红外光源由单片机控制，每隔3 s切换下一个波长，单片机型号为STM32F103c8T6。光源的灯珠购自深圳市超越光电有限公司，选取了780、805、825、850、880和940 nm 6种波长[4]以用于后续多波长近红外图像处理。多波长近红外光源灯珠摆放位置示意图及其实物图如图1所示，近红外静脉图像采集系统实物图如图2所示。

图1 多波长近红外光源灯珠摆放位置示意图及其实物图

图2 近红外静脉图像采集系统实物图

2 二维血管增强与提取

由于近红外静脉图像不够清晰，会影响到对血管边界的提取，本研究先对850 nm波长近红外光下拍摄的静脉图像进行图像增强[5]，以利于下一步将静脉与背景进行分割，明确静脉的二维位置，然后细化出静脉中心线，更好地辅助医护人员进行静脉穿刺。本研究进行图像处理的平台为MATLAB。

2.1 图像增强

对比度受限自适应直方图均衡化（contrast limited adaptive histogram equalization，CLAHE）处理对于医学领域的图像，尤其是医学红外图像的增强效果尤为明显。本研究中CLAHE程序运行时间（11～16 s）由tic、toc（MATLAB内的计时程序，可以定量计算完成制定程序所消耗的时间资源）显示。手背静脉原图如图3（a）所示，经过CLAHE处理后如3（b）所示。

2.2 血管与背景的分割

为了把血管轮廓完全从背景中分割出来，本研究对图像进行了二值化处理。二值化处理目的是将一定范围内的灰度突出成灰度级内最大值，将其他所有灰度降低至0。二值化处理可以简化图像，减小图像所包含的数据量，能够很方便地突出感兴趣目标的轮廓[6]。本研究通过随机采集感兴趣区域的灰度值，取其平均值作为阈值来实现二值化。二值化处理后的近红外静脉图像如图3（c）所示。

观察图3（c）可以看出，血管的检测情况良好，但是血管边缘呈锯齿状，毛刺颇多，所以对该图做了一些形态学处理[7]。形态学处理中，膨胀和腐蚀是常用的方法[8-9]。本研究使用5×5的矩形算子对二值化处理结果图先进行开运算再进行闭运算，处理结果如图3（d）所示。由结果图比对可以看出，形态学处理后血管边缘细小毛刺得到平滑，而且血管形态流畅，完整度高。

2.3 细化处理

图像的细化通常是指将一幅二值化图像骨架化的操作。本研究中对图3（d）做细化处理，即提取血管的中心线条轨迹，将血管部分细化成只有一个像素宽度的线条，如图3（e）所示。细化结果可以显示血管走向，有助于标识出静脉穿刺部位。

图3 二维处理前后的近红外静脉图像

3 多波长近红外静脉深度恢复数学模型的建立与计算

手背静脉三维定位的关键在于血管深度的获取，以下详细介绍通过构建多波长近红外静脉深度恢复数学模型计算血管深度并进行验证的过程。

3.1 手背皮肤光学模型的建立

比尔-朗伯定律指出，电磁波穿过物体时被吸收的量与光程中存在的吸光分子数目成正比[10]。具体吸收量与被穿透物体的吸光系数、厚度等因素有关[11]。吸光系数即单位光程中存在的吸光分子数目的具象表示。比尔-朗伯定律适用于所有的电磁辐射和所有的吸光物质，包括气体、固体、液体、分子、原子和离子[12]。因此，本研究在开发血管深度计算算法时，参考了该定律，在此基础上结合实际情况，建立了不同波长近红外光下的手背皮肤光学模型。

皮肤位于人体表面，人身体各处皮肤形态颇有差异，本研究主要以手背静脉网上方皮肤为例。本研究中手背皮肤光学模型的分层主要可以分成3层：上皮层、真皮层和皮下脂肪层，透过脂肪再向下就是本研究的待测目标——手背静脉。在该手背皮肤光学模型中，血管位于皮下的深度即上皮层、真皮层和皮下脂肪层厚度之和。由皮肤分层可以看到，该深度的波动主要由皮下脂肪层厚度的改变导致，所以本研究默认上皮层和真皮层的厚度固定不变，皮下脂肪层厚度上下波动。

在同一波长的近红外光垂直照射下，上皮层、真皮层和皮下脂肪层的吸收系数各有不同，在不同波长的近红外光垂直照射下，各层吸收系数也不相同。近红外光于皮肤各层的吸收系数是很多研究方向都需要用到的数据，所以在有规律可循的层面已经有学者给出了计算方法，没有规律的层面则通过测量得到。如福建省激光精密加工工程技术研究中心的研究人员在进行老化皮肤光分布的蒙特卡洛模拟时，总结了上皮层、真皮层和血管吸收系数的计算方法[13]。2018年，李博等[14]在与肥胖相关的近红外光谱研究中测量了脂肪层在近红外波段的数据，计算中用到的氧合血红蛋白、水等数据来源于上海化工数据共享服务平台。

在各波段的近红外静脉图像上取血管中间的灰度值，取其平均值得到数据如下：780 nm灰度值为130，805 nm灰度值为127，825 nm灰度值为121，850 nm灰度值为100，880 nm灰度值为113，900 nm灰度值为94。根据图像的对比度和近红外光在静脉中的传播特性，从中挑选805、850、940 nm 3个波长为计算模型的建立提供数据。

3.2 多波长近红外静脉深度恢复数学模型的建立

本研究中多波长近红外静脉深度恢复数学模型搭建步骤如下：（1）对同一个待测血管搭建不同波长的近红外光照射环境就会产生一个参数矩阵，因为皮肤模型分为3层，所以从多波长近红外图像中挑选3个波长的图像用来产生3×3的参数矩阵。（2）在每个波长的近红外静脉图像中分别获取血管中间点的灰度值，经处理后与上述参数矩阵组成增广矩阵。（3）对该增广矩阵进行求解，得到各层的深度计算值，将超声血管成像设备提供的血管深度及血管直径作为参考进行辅助运算和误差修正。（4）仍然使用该待测血管，另取其他部位样本点进行验证运算，计算取样部位的血管深度，同样使用超声图像进行验证。多波长近红外静脉深度恢复数学模型如下：

式中，n为灰度级最大值（本文为255）；m为灰度值；d1、d2、d3、d4分别表示上皮层、真皮层、皮下脂肪层的厚度和血管直径，单位均为分别表示波长为λ时上皮层、真皮层、皮下脂肪层和血液的吸收系数。

多波长近红外静脉深度恢复数学模型的超声血管图及灰度值取样示意图如图4所示，其中图4（a）是用来作为标准的超声血管图，由图可得：d4=0.141 cm，因此建立的方程组如下：

在d1+d2+d3=0.15 cm的约束条件下解该增广矩阵并进行误差修正，对结果保留2位小数，可得d1=0.01 cm、d2=0.09 cm、d3=0.05 cm。

3.3 深度恢复计算结果

在MATLAB中确定805 nm波长的灰度图像左上角坐标为（338，755）、右下角坐标为（352，821），根据图4（b）中红色矩形标记所示的矩形样本读取样本灰度值，使用该值计算脂肪层厚度dcalculate。脂肪层厚度dcalculate与d1、d2的和组成血管深度模拟值dsimulation，使用血管深度模拟值dsimulation与其在原图中的坐标数据（即位置信息）绘制三维散点图（如图5所示）。计算血管深度模拟值dsimulation与超声显示血管参考深度的误差，误差数据柱状图如图6所示，误差数据三维散点图如图7所示。

图4 多波长近红外静脉深度恢复数学模型的灰度值取样示意图及参照标准（超声血管图）

图5 血管位于皮下深度三维散点图

图6 误差数据柱状图

从图5中可以看出，整体结果分布较为集中，第338行至339行波动较大。从图6中可以看出，有90.85%的数据误差≤0.015 cm，有96.92%的数据误差≤0.021 cm。从图7中可以看出，在1 005个验证点中，误差较大的点多位于样本边界处，推测取样时可能超出了血管边界，导致该行附近数据误差较大。去除第338行和339行样本数据，使用（340，755）到（352，821）共计871个样本点重新计算误差，得到的误差修正柱状图如图8所示。

图7 误差数据三维散点图

图8 误差修正柱状图

从图8中可以看出，去除第338行和339行样本数据后，误差最大值已经由0.041 cm降至0.033 cm。误差修正后，有91.50%的数据误差≤0.015 cm，有97.82%的数据误差≤0.021 cm，较修正前有少许提升。由误差分布规律以及计算原理可以得出，越靠近血管中心轨迹，计算误差越小。

4 结语

本研究通过建立多波长近红外静脉深度恢复数学模型，恢复了二维场景中丢失的深度信息，从而对手背静脉进行了三维定位。经871个样本点的计算和验证，有97.82%的把握将血管深度测量的误差控制在0.021 cm以内。本研究为现行市场上的二维投影式静脉显像仪升级为三维提供了参考[15]，将来可能为疫情下的医护人员和患者提供更多帮助。

本研究仍存在不足之处，例如并未完全脱离医护人员实现智能自动静脉穿刺，尚需人为参与等。因此在后续研究中，我们会致力于开发自动计算并判断最佳穿刺位置的医疗器械，为其配备相应的硬件，使其能够按照静脉定位数据做到自动配药[16]、自动静脉穿刺以及自动给药，为医护人员提供更大的帮助。