基于激光散斑衬比成像技术的脑血流分析

邬丹丹，姚 康，管凯捷，付威威，3*，孙建美，郭云宝，王 驰

(1.上海大学 机电工程与自动化学院，上海 200444；2.中国科学院 苏州生物医学工程技术研究所，江苏 苏州 215000；3. 苏州国科医工科技发展(集团)有限公司，江苏 苏州 215000；4.吉林大学 白求恩第一医院，吉林 长春 130021)

1 引 言

脑血管疾病是人类第二大致死原因，基本病因包括血管病变、血液动力学改变等。脑血管手术中，需要对术野范围内的目标血管进行监测，监测指标包括管径、流速、流量、血管吻合后血流流通和脑部血流灌注情况等。传统的血流监测方法，如核磁共振成像[1]，时空分辨率较低；如血管造影技术，需注射造影剂，存在副作用，并且会在人体内代谢，有效工作时间短。目前，脑血管手术中多采用激光多普勒方法进行单点流速测量，但其实时成像能力较低，并且无法获取血管管径来测定血流量。而传统的光谱成像方法难以精确测量血管内径，不能监测血管内的血流速度及其变化，更无法反映术野内整体血流动力学改变。

近年来，作为一种全场光学成像技术，激光散斑衬比成像(Laser Speckle Contrast Imaging， LSCI)技术逐渐发展起来[2-3]。20世纪80年代，Fercher等[4]首次应用LSCI技术测量视网膜微血管流速。Naoki等[5]通过LSCI技术发现青光眼患者的神经乳头周围的大血管流速低于常人。在皮肤血流成像方面，贾亚威等[6]应用LSCI技术进行中医理疗功效检测；张洁等[7]应用LSCI技术观测针灸对局部血流微循环的改善效果。在脑皮层血流成像方面，Nils等[8]利用LSCI技术监测人脑的局部缺血，并预测发生梗塞的脑组织范围；李苗等[9]基于LSCI技术及支持向量机获得脑皮层深穿枝血管的分布信息。

针对脑血管手术中对脑血流分析的实际需求，本文提出了基于LSCI技术的脑血流分析方法，实现了LSCI技术与血管管径测量算法的有机结合。设计并搭建了LSCI系统采集激光散斑图像，采用改进的时间衬比分析和伪彩色处理得到便于血管分割及边缘提取的血流速度图像，并反映出不同血管及同一血管在不同时刻的流速变化。手动选取两点作为感兴趣区域(Region of Interest， ROI)的基准点，以截取出管径测量所需的ROI，对ROI图像预处理后利用改进的Canny算法提取血管边缘线并计算血管管径。基于LSCI技术的脑血流分析方法具有非接触、高分辨率及可连续监测等优点，可反映出单一血管和术野内总体血流速度变化，呈现出手术过程中不同操作引起的血液动力学的动态变化。同时，可以测量得到血管内壁的有效流通直径，能够结合激光多普勒血流仪计算得到精确的血流量，监测流量的变化情况，为临床提供参考。

2 LSCI系统的工作原理及设计搭建

将LSCI技术与血管管径测量算法结合，在监测血流速度变化的同时，可测量指定血管段的有效流通直径，能够结合激光多普勒血流仪给出的速度值计算得到精确的血流量，监测血流量变化。

2.1 LSCI系统的工作原理

(1)

根据动态光散射近似模型，可将散斑衬比值K表示为相干光去相关时间τ和曝光时间T的函数[11]，如式(2)所示：

(2)

当T≫τ时，衬比值K与τ呈近似线性关系。而τ=λ/(2πν)，其中ν为散射粒子运动速度，λ为相干光波长。当T≫τ时，衬比值K与ν呈反比关系，这是利用LSCI技术进行血流成像及速度定量分析的基础。理想状态下，衬比值K在[0，1]内，当散射介质静止时，K为1；当散射介质运动时，K逐渐趋于0，且运动速度越快，K越小。

2.2 LSCI系统的设计搭建

LSCI系统结构如图1(a)所示，通常由照明光源、相机、镜头和计算机等组成。本文根据系统工作原理和结构图设计搭建了如图1(b)所示的LSCI系统，其中照明光源采用激光二极管(850 nm，100 mW)；相机采用CMOS相机(acA1920-155 μm，Basler，德国)，分辨率为1 920×1 200 pixel。

(a)LSCI系统结构示意图

(b)LSCI系统实物图

人脑部血管壁厚不超过0.5 mm，而LSCI系统的采样深度在0.3～1 mm[12]，所以为了验证本文LSCI系统的可行性，设置了散斑采样深度测量实验。配制1%体积比Intralipid溶液用于模拟生物组织血流。利用0.6%体积比的Intralipid与0.4%质量比的琼脂糖配制出混合溶液，加热使之充分混合，冷却凝固后用于模拟血管周围的组织。盛放凝胶的实验模具如图2(a)所示。实物图如图2(b)所示，一根内径为0.5 mm的毛细玻璃管倾斜固定于玻璃器皿中，调节蠕动泵控制模拟血流的流动速度。在固定倾斜角度θ下，处理得到的伪彩色血流图像如图2(c)所示。根据图像中水平方向上可观察到的像素数量N、工作距离S、镜头焦距f和相机芯片像元尺寸m计算可得到采样深度Δh，即有：

(3)

其中相机芯片的像元尺寸为5.86 μm。统计得到的深度曲线如图2(d)所示，表示不同流速下计算得到的采样深度。实验结果表明，本文设计搭建的LSCI系统的采样深度约为0.65 mm，可满足应用需求。

3 改进型LSCI系统算法

改进型LSCI系统算法包括对散斑图像采用改进的时间衬比分析和伪彩色处理得到血流速度图像，然后手动选取感兴趣区域(Region of Interest，ROI)的基准点截取ROI，对ROI图像预处理后利用改进的Canny算法进行边缘提取并计算血管的有效流通直径。

(a)模具(a)Mold

(b)实物图(b)Physical map

(c)伪彩色血流图像

(d)深度曲线

3.1 改进的时间衬比分析

LSCI系统采集得到的散斑图像如图3(a)所示。比较分析常用衬比分析方法后，本文采用时间衬比分析方法，通过计算50帧散斑图像中相同位置处的像素灰度值的平均值作为衬比值，得到如图3(b)所示的衬比图像。相关实验结果表明[13]，当散斑图像帧数超过15帧时，使用时间衬比分析得到的速度值与实际速度值的相关度大于0.96。此时衬比图像整体偏暗，人眼无法很好地辨别出血流分布。人眼对彩色的分辨能力远高于对灰度的分辨能力，所以对衬比图像进行伪彩色处理，将灰度映射到伪彩色空间中，提高图像的可辨识度和细节的辨别力，增强血管边缘的清晰度，提高边缘提取的准确度。相关研究表明[14]，K2一般在[0，0.1]内，此时1/K2与运动速度ν成正比关系。为了提高图像的可视化效果，先由衬比图像计算得到1/K2图像，随后将它映射到JET颜色模型。理想状态下，散射介质静止时，其颜色为蓝色；散射介质运动时，其颜色随着运动速度的加快趋向红色。对时间衬比图像进行伪彩色处理得到如图3(c)所示的伪彩色图像，可以看出明显的血流分布，血流速度越快，颜色越红，反之越蓝。但是，此时的伪彩色图像由于采用多帧图像进行衬比计算而降低了时间分辨率。

实际应用中，为了提高时间衬比图像的信噪比，并考虑到生物及环境噪声等，需要对它进行空间平均。本文采用均值滤波对图像进行空间平均，滤波窗口大小为3×3 pixel。对时间衬比图像进行空间平均及伪彩色处理得到如图3(d)所示的血流图像，可以看到空间平均消除了部分噪声，但是也引起了空间分辨率的下降。为了平衡图像分辨率和信噪比，本文在对时间衬比图像进行空间平均前，先采用形态学闭运算进行图像增强及部分噪声消除。形态学闭运算是指图像依次经过膨胀、腐蚀处理的过程。图像被膨胀后，不连续的边缘像素点连通形成闭合边缘，增强了微小信息；接着进行腐蚀处理，去除图像噪声点，细化图像边缘。进行形态学闭运算时，需根据实际图像的结构特征选择合适的结构元素及其尺寸，常用的结构元素包括矩形、圆盘、十字叉结构等。在经过试验后，本文采用符合血管形态的圆盘结构进行形态学闭运算，其中膨胀处理的圆盘结构半径为2，腐蚀处理的圆盘结构半径为1。对散斑图像进行改进的时间衬比分析，结合形态学闭运算和空间平均，再进行伪彩色处理后得到如图3(e)所示的改进的血流速度图像。可以看出，图3(e)相比图3(c)增强了微小信息，整体改进了血流分布的显示效果，便于图像分割及边缘提取。

(a)散斑图像(a)Speckle image

(b)时间衬比图像(b)Time contrast image

(c)伪彩色图像(c)Pseudo color image

(d)血流图像(d)Blood flow image

(e)改进的血流速度图像

3.2 血管管径测量

血管管径测量方法主要分为两类：(1)基于血管两侧边界，通过提取血管网络和中心线计算管径值。例如，Li等[15]采用基于信息融合的边缘检测器获取血管边缘，生成血管中心线的垂直线计算管径。该类方法通过血管边缘线和中心线计算管径，边界识别不准确及中心线位置的偏差会导致管径的测量误差。(2)基于血管横截面的灰度分布特性统计曲线，采用不同曲线模型对它进行拟合，计算曲线半高度等宽得到血管管径。例如，Lowell等[16]采用改进的高斯微分滤波器拟合血管横截面的灰度分布曲线。该类方法与选定的血管横截面位置及选取的拟合曲线模型等有着密切关系，参数调节复杂，不利于临床操作。综合考虑上述血管管径测量方法存在的问题，本研究基于改进的血流速度图像进行管径测量，算法流程如图4所示。

图4 血管管径测量算法流程

3.2.1 图像预处理

图像预处理包括在改进的血流速度图像的基础上手动选取ROI，灰度化后利用自适应对比度增强(Adaptive Contrast Enhancement， ACE)算法进行局部增强，最后进行OTSU阈值分割得到二值图像。具体步骤如下：

(1)ROI选取

实际应用中，医生需要得到某一特定血管段内的有效流通直径，所以需要进行ROI选取。如图5(a)所示，本文通过在血管段两侧手动选取A，B两点作为ROI的基准点，截取出指定大小的ROI，如图5(b)所示。通过手动选取基准点截取ROI，将管径测量区域局限在较小范围内，血管边缘线可视为两条平行直线，避免了背景噪声及其他血管干扰等带来的管径测量误差，极大地减少了图像处理的数据量，便于临床操作。

(2)ACE增强

ROI的灰度图像如图5(c)所示。为了提高图像分割的准确度，本文在进行阈值分割前先采用ACE算法进行图像的局部增强，通过差分计算目标点与周围像素点的相对阴暗关系来校正最终像素值，拉升对比度来调整图像的动态范围。处理得到的ACE图像如图5(d)所示，可以看出，ACE算法调整了图像的对比度，充分考虑了局部信息，并且限制了噪声的放大，避免了低频背景的干扰，对待测量的血管段具有较高的增强效果。

(3)OTSU阈值分割

OTSU算法是一种自适应的阈值确定方法，又称作最大类间方差法。它是按照图像的灰度特性，用阈值将图像分为背景和目标两部分，并计算背景和目标之间的类间方差，遍历图像得到使类间方差最大的阈值T0，取T0对图像进行分割得到二值图像，如图5(e)所示。可以看出，图像分割效果较好，从背景中提取出了待测量的血管段，消除了噪声。

(a)选取示意图

(b)感兴趣区域(b)ROI

(c)灰度图(c)Grayscale image

(d)增强图像(d)Enhancing image

(e)二值图像(e)Binary image

3.2.2 改进的Canny算法

得到二值图像后，需进行边缘提取得到血管边缘线。在对常用边缘提取算子进行尝试后，本文采用不易受噪声干扰，且能够检测到弱边缘的Canny算子提取血管边缘，其简要步骤如下：

(1)利用高斯滤波器对图像进行平滑，去除噪声；

(2)采用一阶偏导的有限差分计算图像中每个像素点的梯度幅值和方向；

(3)对梯度幅值进行非极大值抑制，将当前像素与沿正负梯度方向上的两个像素进行梯度幅值比较，如果当前像素的梯度幅值最大，则保留为边缘点，否则抑制该点；

(4)采用双阈值算法检测和连接边缘，设置高、低阈值将边缘像素标记为强、弱边缘像素，抑制低于低阈值的边缘像素，其次查看弱边缘像素的8个邻域像素，如果其中存在强边缘像素，则保留该像素点。

Canny算子边缘提取结果如图6(a)所示，与二值图像的边缘像素点逐个对比发现，此时提取得到的边缘对二值图像边缘进行了局部放大，导致测量得到的血管管径值略微偏大，并且Canny算子会导致边界部分的边缘缺失。针对以上问题，本文在Canny算子边缘提取前对图像进行了腐蚀操作，并将Canny算子与Sobel算子进行有机结合，生成改进的边缘提取图，如图6(b)所示。可以看出，改进的Canny算法补齐了边界部分的边缘，并且抑制了局部放大，提高了边缘提取的准确度。然后根据血管边缘方向逐列或逐行查找边缘线上两端点并标记，如图6(c)所示，在某一列标记了C，D两端点。遍历图像重复这一过程标记出边缘线上全部端点，即由红色十字标记的上端点和绿色三角形标记的下端点，如图6(d)所示(彩图见期刊电子版)。逐列计算上、下端点的间距得到多个距离，对其求均值作为端点间平均距离。由于手动选取ROI，管径测量区域局限在较小范围内，血管边缘线可近似为平行直线，所以分别对上、下端点数据进行最小二乘法拟合得到红色和绿色直线，如图6(e)所示，并计算出相应的斜率，利用两个斜率的均值和端点间平均距离计算血管管径。

(a)Canny算子(a)Canny operator

(b)改进的Canny算法(b)Improved Canny algorithm

(c)标记示意(c)Marking diagram

(d)标记(d)Marking

(e)拟合(e)Fitting

为验证改进的Canny算法对血管管径测量误差的降低程度，通过公式计算相对误差(Relative Error，Er)，即有：

Er=(M-T)/T，

(4)

其中：M为图6(d)中上、下边缘端点间的平均距离，T为图5(e)二值图像中上、下边缘端点间的平均距离。表1是针对采用改进的时间衬比分析和伪彩色处理对散斑图像处理得到的50幅血流速度图像，对每幅图像手动选取不同位置的10个ROI，在CPU为i7-8700 3.2 GHz，内存为8 GB的计算机上运行不同算法得到的Er对比结果。可以看出，改进的Canny算法降低了血管管径测量的相对误差，其相对误差可达5.1%，低于其他算法。

表1 血管边缘提取方法比较

4 动物实验及结果分析

人脑部血管直径不超过5 mm，兔子耳背中央大血管直径最粗在5 mm左右，并考虑到人脑部血管壁厚(≤0.5 mm)，两者具有极大的相似性，所以设置了动物实验采集兔子耳背的散斑图像。实验过程中所有操作均得到中国科学院苏州生物医学工程技术研究所实验动物伦理委员会的批准，并且符合《江苏省实验动物管理方法》的要求。

4.1 实验对象和方法

实验对象为购自江苏省苏州市镇湖实验动物科技的健康成年实验兔，依据体型将实验兔分为特大、正常和小型实验兔，每种体型的实验兔各两只。进行实验前需适应性喂养3天，实验环境温度为18.0～28.0 ℃。为避免对血管造成不必要的伤害且便于操作，利用剃毛膏在兔子耳背上清理出观测窗口。随后将兔子静置半小时左右，等待其生理状态稳定。实验过程中，兔子始终固定在兔盒中，并放置于防震光学平台上，如图7所示。拍摄散斑图像时，可将少量生理盐水滴到兔子耳背的观测窗口处以减少镜面反射的影响，提高图像的散斑效果。

图7 动物实验示意图

4.2 实验结果及分析

4.2.1 血流速度变化监测

按住兔子耳背中间大血管后，整个伪彩色图像的颜色趋向蓝色，表示血流速度大幅度降低，出现了血流阻断现象，如图8(a)所示。放松对血管的按压后，伪彩色图像的颜色又趋向红色，如图8(b)所示(彩图见期刊电子版)。可以看出，图8(b)与图3(e)所显示的血流分布情况基本一致，即出现了血流灌注现象。图3(e)、图8(a)及图8(b)的动态循环过程验证了LSCI系统监测实际血流速度变化的可靠性。实验结果表明，此时的伪彩色血流速度图像处理及显示帧率不超过每秒12帧。

(a)阻断图像(a)Blocking image

(b)再灌注图像(b)Reperfusion image

4.2.2 血管管径测量

本文采用改进的时间衬比分析和伪彩色处理对散斑图像处理得到血流速度图像，手动选取不同位置的ROI进行管径测量。选取一幅具有代表性的改进的血流速度图像，在该图中手动选取3个不同位置的ROI，如图9(a)所示，截取出的ROI如图9(b)所示，提取出的血管边缘线如图9(c)所示。

如表2所示，对3种体型的实验兔分别采用本文提出的管径测量法、人工确定血管直径法及LSCI软件进行测量，得到的血管管径值依次记为DO，DR，DL。其中，人工确定血管直径法是在软件中对正常光照下拍摄的兔子耳背图像进行血管边界标识后计算管径，并通过软管实验测得该方法的平均误差为0.19 μm；LSCI软件为RFLSI Ⅲ型激光散斑血流成像仪(深圳市瑞沃德生命科技)所用软件。从表2中数据可以看出，3种方法测量得到的血管管径值具有较高的一致性，本文算法与人工确定血管直径法的平均误差为2.2 μm，LSCI软件与人工确定血管直径法的平均误差为5.3 μm，所以本文算法测得的管径值更接近人工确定血管直径法测得的管径值，优于LSCI软件测得的管径值。

(a)ROI选取(a)ROI selection

(b)ROI(b)ROI

(c)血管边缘(c)Vessel edge

表2 管径测量结果比较

5 结 论

本文提出了基于LSCI技术的脑血流分析方法，实现了LSCI技术与血管管径测量算法的结合，通过改进的时间衬比分析和Canny算法对血管管径测量方法进行了改进。实验结果表明，本文研究的脑血流分析方法可有效监测血流速度变化，同时测量血管的有效流通直径，结合激光多普勒血流仪能够实现目标血管的流速及流量的精确监测，为临床提供参考，具有较高的实用价值。下一步研究的重点是利用现场可编程门阵列(FPGA)图像采集卡实现术中脑血流的实时监测和多个ROI的实时管径测量。