CN105286779B - 一种在体视网膜血流动力学的成像与绝对流速测量方法 - Google Patents

Abstract

一种在体视网膜血流动力学的成像与绝对流速测量方法，本发明提出基于常规FD·OCT，整合OCT血管造影技术和多普勒OCT成像技术，实现视网膜血管绝对流速和血流动力学的测量。该发明通过OCT血管造影技术重建三维血管几何结构分布，根据三维血管的几何结构，计算血管与光束的夹角，即多普勒角度，采用多普勒OCT技术获得血流的多普勒频移，并根据计算的多普勒角度进行矫正，从而获取绝对流速，解决视网膜血管多普勒角度的难题；该发明通过OCT血管造影技术重建的三维血管结构，结合特定位置的绝对流速值，通过模拟血液在血管中的流动，可实现在体视网膜血流动力学的测量。

Description

一种在体视网膜血流动力学的成像与绝对流速测量方法

技术领域

本发明涉及眼科影像学领域，特别涉及一种在体视网膜血流动力学的成像与绝对流速测量方法。

背景技术

视网膜灌注或流速的改变在青光眼、糖尿病视网膜病变等眼科疾病的发病机制中起着重要的作用，因此在在体视网膜血流动力学的测量是眼科影像学领域一个非常重要的研究内容。视网膜血流动力学的定量评估可以为眼病的病理变化以及早期诊治提供非常有用的信息。

目前，有多种成像技术可以用来评估视网膜血管流速，如视频荧光素眼底血管造影、激光散斑流速成像仪、激光多普勒测速仪、激光扫描多普勒流速仪、超声彩色多普勒成像技术等。然而，现有的这些临床方法各有其局限性。视频荧光素眼底血管造影和激光散斑流速成像仪虽然可通过动静脉通过时间和平均模糊率等间接指标表达视网膜血流变化情况，但无法直接测量流速；激光扫描多普勒血流仪(例如，海德堡视网膜血流仪)可以在二维扫描区域沿探头波束方向形成流速分量的地形图；然而，获取的流速分量只是一个流量的相对变化，无法获取绝对流速数值；彩色多普勒超声可以对血管流速进行横断面和三维体成像，但～200μm的空间分辨率，限制了其对视网膜中央血管微细结构的显像能力，从而降低了其对血流测量的准确性。因此，目前客观测量在在体视网膜血流动力学仍面临极大的挑战。

傅里叶域光学相干层析成像(FD-OCT)采用相位分辨可实现多普勒成像，以更简单和更为可靠地方法测量具有深度分辨信息的血管结构和流速，但早期报道的多普勒OCT技术只能探测到沿光束方向的流速分量，在探测光束与流速的夹角(多普勒夹角)未知的情况下，无法测得流体的绝对流速，只能获得流速的相对变化。另外，在视网膜流体测量上，因视网膜大多血管流速方向与光束方向垂直，基本无法获取多普勒频移的信息。为解决这些技术困难，近些年，有研究者提出双光束方向的多普勒OCT技术用于实现视网膜绝对流速的测量。这种方法利用两束不同入射光线方向的光束测量各自的多普勒流速分量，然后通过三角关系计算绝对流速，该方法解决了多普勒角度测量的困难，但是因为双光束扫描需要复杂的硬件改装，这种方法不易移植到现有商业化的OCT设备上。

发明内容

为了解决现有技术存在的问题，本发明提供了一种在体视网膜血流动力学的成像与绝对流速测量方法，包括以下步骤：

(1)OCT血管图像造影；

(2)采用FD-OCT获取视网膜视乳头区的三维结构数据；

(3)二维血管信息提取；

(4)三维血管几何结构的重建：根据二维血管的几何特征，初步估计血管的位置，及中心线的位置，进而对根据二维血管的位置对三维体数据进行重采样，获取单根血管的三维体数据，通过对体数据进行展开，并应用基于扩展的动态规划的三维分层算法获取血管边界，最终获得眼底血管的三维形态结构数据，实现视网膜ONH区的3-D血管几何结构的重建；

(5)获取OCT的多普勒图像：采用环形扫描或水平扫描方式获取OCT的多普勒图像；

(6)流速校正：根据步骤(3)获取的血管的几何边界数据，使用最小二乘法拟合血管边界，从而获得血管中心位置，可以获得沿着血管中心线方向需要计算的血管位置的方向，从而获得入射光线与血管之间的夹角，即多普勒角度(θ)，对步骤(4)获取的血流分量进行角度矫正：

从而获取视网膜ONH区血管特定位置的绝对流速。

所述的步骤(2)中FD-OCT视网膜视乳头区图像纵向分辨率为3μm，采集速率不小于50kHz。三维血管造影数据由512×128×8条A扫描组成。

所述的步骤(2)中三维结构数据获得的方法为：采用方块扫描方式获取三维结构数据，每个三维数据由128幅断层图组成，每幅断层图像有在同一个位置8次重复B扫描组成，然后采用基于强度信息的多普勒方差方法通过在同一个位置做重复扫描，利用相邻两次B-scan之间强度的变化，实现对微血管的成像。所述的基于强度信息的多普勒方差方法的算法如下：

其中M代表平均扫描深度的像素点，N代表同一位置扫描的次数。

所述的步骤(3)二维血管信息提取方法为：

a、首先通过基于扩展的动态规划的三维分层算法，获取视网膜内界膜的表面，基于三维表面边界自动探测的过程主要分为三维数据预处理和利用最短路算法寻找边界，在读入视网膜三维结构数据后，首先进行简单的平滑滤波，计算方法如下式所示：

其中v(x，y，z)表示在横向x，纵向y及深度方向z的强度值，v′(x，y，z)为滤波后的强度值，(a，b，c)分别为在三个方向上滤波核的大小；

对于三维OCT数据v(x，y，z)，基于三维扩展的状态转移方程可以表示为：

其中C₁和C₂分别表示在x-z平面和y-z平面上的路径的累计费用，参数d₁，d₂和a₁，a₂用于控制探测的目标曲面的平滑度；w(x+i，x)则表示点(x+i，y-i，z)到点(x，y，z)的路径的权重，可以由下式获得：

w(x+i，x)＝2-r(x+i，y-1，z)-r(x，y，z)

其中r(x，y，z)为点(x，y，z)处的梯度值，并归一化0至1区间；

通过对体数据做最大值投影，获取视乳头部位二维血管形态结构：

p(x，y)＝max{v(x，y，k)|k＝z₁，z₂...z_h}

其中p(x，y)为投影后的二维血管形态结构(图3b)，在深度方向上投影范围为从内界膜起点向下300个像素点范围；

b、对二维的血管采用半自动的基于迪克斯特拉最短路径算法获取二维血管的几何结构，根据相邻像素点之间的亮度值不同，给每条边赋予一定的权重，权重公式如下式所示：

w_ab＝2-w_a-w_b

其中w_ab表示节点a与节点之间的权重，w_a表示节点a的亮度值，w_b表示节点b的亮度值，通过输入血管的起点(s)和终点(e)的位置，采用戴克斯特拉算法获取起点与终点之间血管的位置，初始时，起点(s)的路径长度为0，若存在能达到边(s，m)，则把起点(s)到节点(m)的距离设置d_s，m，不能到达则设置为无穷大，不断对外接变进行拓展，若存在节点(n)使得d_s，n+w_n，m＜d_s，m，则更新d_s，m，且设置节点(m)上一节点的为节点(n)。从而可以获得起点到终点之间的最短路径，即为需要探测的血管的位置。

所述的步骤(5)的OCT的多普勒图像由1024个像素点组成，采集图像的时间＞4秒，保证能测量在不少于4个心动周期内的视网膜血流的搏动，通过分析相邻两个A扫描直接的相位差，如下式所示：

其中为计算所得的相位差，Im(A_j，z)表示在第j条线，深度为z位置处虚部，Re(A_j，z)则为实部。J和N分别为横向A扫描和深度方向Z的平均数，沿着深度方向的流速分量v_z为：

其中且λ_C为光源中心波长，n为视网膜的折射率，ΔT为扫描速度的倒数。

本发明的有益效果是：本发明提供了一种在体视网膜血流动力学的成像与绝对流速测量方法，本发明提出基于常规FD-OCT，整合OCT血管造影技术和多普勒OCT成像技术，实现视网膜血管绝对流速和血流动力学的测量。该发明通过OCT血管造影技术重建三维血管几何结构分布，根据三维血管的几何结构，计算血管与光束的夹角，即多普勒角度，采用多普勒OCT技术获得血流的多普勒频移，并根据计算的多普勒角度进行矫正，从而获取绝对流速，解决视网膜血管多普勒角度的难题；该发明通过OCT血管造影技术重建的三维血管结构，结合特定位置的绝对流速值，通过模拟血液在血管中的流动，可实现在在体视网膜血流动力学的测量。

附图说明

图1为本发明功能实现流程图。

图2为本发明步骤(3)三维血管提取方法流程图。

图3为交互式二维血管几何结构提取方法，其中图3a为二维血管结构提取流程图，图3b和3c分别为二维血管造影投影图像和经过滤波后图像，图3d和3e)示血管在垂直和水平情况下自动探测算法能良好的探测血管的位置。图3f为眼底主要大血管二维结构探测结果。

图4为本发明原理示意图，其中图4a和4d分别为OCT血管造影和多普勒OCT扫描模式，图4b为OCT血管造影在视乳头部位血管造影图，图4c为眼底视网膜主要大血管图像分割结果。图4e和4f为多普勒OCT结果，颜色流速的大小，图4g为多普勒OCT测量血管血流随着心脏的跳动发生波动情况，红线为矫正多普勒角度前，蓝色为矫正多普勒角度后。

图5为FD-OCT获取的多普勒OCT图像和OCT血管造影图。

图6为基于OCT血管造影数据重建的视网膜三维血管结构。

具体实施方式

本发明整合OCT血管造影与多普勒OCT流速检测技术，提出一种在在体视网膜血流动力学的成像与测量方法。其包括以下步骤：

OCT血管图像造影；

1、采用FD-OCT获取视网膜视乳头区(ONH)的三维结构数据。涉及的FD-OCT，其图像纵向分辨率为3μm，采集速率不小于50kHz。三维血管造影数据由512×128×8条A扫描组成，采用方块扫描方式获取。每个三维数据由128幅断层图组成，每幅断层图像有在同一个位置8次重复B扫描组成。基于强度信息的多普勒方差方法通过在同一个位置做重复扫描，利用相邻两次B-scan之间强度的变化，实现对微血管的成像。基于强度信息的多普勒方差方法的算法如下图式所式：

其中M代表平均扫描深度的像素点，N代表同一位置扫描的次数。当M和N增加时，性噪比(SNR)增加。在发明中，我们将M设为2，将N设为8。A_n，m代表OCT数据的幅度值。最后采用基于平均强度直方图分析的阈值来消除噪声。

2：基于获取的血管造影体数据，如图2所示，首先通过基于扩展的动态规划的三维分层算法，获取视网膜内界膜的表面。基于三维表面边界自动探测的过程主要分为三维数据预处理和利用最短路算法寻找边界。在读入视网膜三维结构数据后，首先进行简单的平滑滤波，如下式所示：

其中v(x，y，z)表示在横向x，纵向y及深度方向z的强度值，v′(x，y，z)为滤波后的强度值，(a，b，c)分别为在三个方向上滤波核的大小，本发明中均设置为2。进一步通过在深度方向上做应用[-1：1]模板对三维体数据做卷积计算，增强视网膜内界膜的边界。

基于三维扩展的动态规划算法的基本原理是将边界探测的问题转换成寻找最小路径费用的问题，动态规划算法可以高效寻找最优路径，并保持较好的鲁棒性。对于三维OCT数据v(x，y，z)，基于三维扩展的状态转移方程可以表示为：

其中C₁和C₂分别表示在x-z平面和y-z平面上的路径的累计费用，参数d₁，d₂和a₁，a₂用于控制探测的目标曲面的平滑度。w(x+i，x)则表示点(x+i，y-i，z)到点(x，y，z)的路径的权重，可以由下式获得：

w(x+i，x)＝2-r(x+i，y-1，z)-r(x，y，z)

其中r(x，y，z)为点(x，y，z)处的梯度值，并归一化0至1区间。

通过对体数据做最大值投影，获取视乳头部位二维血管形态结构：

p(x，y)＝max{v(x，y，k)|k＝z₁，z₂...z_h}

其中p(x，y)为投影后的二维血管形态结构(图3b)，在深度方向上投影范围为从内界膜起点向下300个像素点范围。

进而对二维的血管采用半自动的基于迪克斯特拉最短路径算法获取二维血管的几何结构，如图3a流程图所示，与基于三维扩展的动态规划算法相似，将每一幅梯度图像的像素点看成是图的一个节点，相邻的像素之间的相互连接为边，根据相邻像素点之间的亮度值不同，给每条边赋予一定的权重，权重公式如下式所示

w_ab＝2-w_a-w_b

其中w_ab表示节点a与节点之间的权重，w_a表示节点a的亮度值，w_b表示节点b的亮度值。通过把通过输入血管的起点(s)和终点(e)的位置，采用戴克斯特拉算法获取起点与终点之间血管的位置。初始时，起点(s)的路径长度为0，若存在能达到边(s，m)，则把起点(s)到节点(m)的距离设置d_s，m，不能到达则设置为无穷大。不断对外接变进行拓展，若存在节点(n)使得d_s，n+w_n，m＜d_s，m，则更新d_s，m，且设置节点(m)上一节点的为节点(n)。从而可以获得起点到终点之间的最短路径，即为需要探测的血管的位置。

进一步，根据二维血管的几何特征，初步估计血管的位置，及中心线的位置，进而对根据二维血管的位置对三维体数据进行重采样，获取单根血管的三维体数据。通过对体数据进行展开，并应用基于扩展的动态规划的三维分层算法获取血管边界，最终获得眼底血管的三维形态结构数据，实现视网膜ONH区的3-D血管几何结构的重建。3：采用环形扫描或水平扫描方式获取OCT的多普勒图像，涉及的环形扫描图像由1024个像素点组成，采集图像的时间＞4秒，保证能测量在不少于4个心动周期内的视网膜血流的搏动。通过分析相邻两个A扫描直接的相位差，如下式所示：

其中为计算所得的相位差，Im(A_j，z)表示在第j条线，深度为z位置处虚部，Re(A_j，z)则为实部。J和N分别为横向A扫描和深度方向Z的平均数。则沿着深度方向的流速分量v_z为：

其中λ_C为光源中心波长，n为视网膜的折射率，ΔT为扫描速度的倒数。

4：根据第1和2步骤获取的3D视网膜血管的几何结构，根据血管的几何边界数据，使用最小二乘法拟合血管边界，从而获得血管中心位置，可以获得沿着血管中心线方向需要计算的血管位置的方向，从而获得入射光线与血管之间的夹角，即多普勒角度(θ)，对步骤3获取的血流分量进行角度矫正：

从而获取视网膜ONH区血管特定位置的绝对流速。

5：结合获取的三维几何结构数据，并结合多普勒OCT测量获得绝对流速，将其导入流体模拟软件(如fluent软件)，评价三维血管形态对血流动力学的影响，如灌注压下降的速度，根据流速的情况预测血流的粘滞度等血流动力学。

如说明书附图所示，图中4a和4d分别为OCT血管造影和多普勒OCT扫描模式，其中OCT血管造影采用方块状扫描模式，多普勒OCT采用了横向位置重复扫描模式。图4b为OCT血管造影在视乳头部位血管造影图，图4c为眼底视网膜主要大血管图像分割结果。图4e和4f为多普勒OCT结果，颜色流速的大小，4g为多普勒OCT测量血管血流随着心脏的跳动发生波动情况，红线为矫正多普勒角度前，蓝色为矫正多普勒角度后。

Claims (5)

1.一种在体视网膜血流动力学的成像与绝对流速测量方法，其特征在于，包括以下步骤：

(1)OCT血管图像造影；

(2)采用FD-OCT获取视网膜视乳头区的三维结构数据；

(3)二维血管信息提取，方法为：

a、首先通过基于扩展的动态规划的三维分层算法，获取视网膜内界膜的表面，基于三维表面边界自动探测的过程主要分为三维数据预处理和利用最短路算法寻找边界，在读入视网膜三维结构数据后，首先进行简单的平滑滤波，计算方法如下式所示：

$I^{\′}(x,y,z)=\frac{1}{(2a+1)(2b+1)(2c+1)}\underset{i=-a}{\overset{a}{\Σ}}\underset{j=-b}{\overset{b}{\Σ}}\underset{k=-c}{\overset{c}{\Σ}}I(x+i,y+j,z+k)$

其中I(x，y，z)表示在横向x，纵向y及深度方向z的强度值，I′(x，y，z)为滤波后的强度值，(a，b，c)分别为在三个方向上滤波核的大小；

对于三维OCT数据I(x，y，z)，基于三维扩展的状态转移方程可以表示为：

$C_1(x,y,z)=\underset{-d_1\≤i\≤d_1}{min}\{C_1(x+i,y-1,z)+{\mathrm{\α}}_1|i|+w(x+i,x)\}$

$C_2(x,y,z)=\underset{-d_2\≤i\≤d_2}{min}\{C_2(x-1,y+j,z)+{\mathrm{\α}}_2|j|+C_1(x,y,z)$

其中C₁和C₂分别表示在x-z平面和y-z平面上的路径的累计费用，参数d₁，d₂和α₁，α₂用于控制探测的目标曲面的平滑度；w(x+i，x)则表示点(x+i，y-i，z)到点(x，y，z)的路径的权重，可以由下式获得：

w(x+i，x)＝2-r(x+i，y-1，z)-r(x，y，z)

其中r(x，y，z)为点(x，y，z)处的梯度值，并归一化0至1区间；

通过对体数据做最大值投影，获取视乳头部位二维血管形态结构：

p(x，y)＝max{I(x，y，k)|k＝z₁，z₂...z_h}

其中p(x，y)为投影后的二维血管形态结构，在深度方向上投影范围为从内界膜起点向下300个像素点范围；

b、对二维的血管采用半自动的基于迪克斯特拉最短路径算法获取二维血管的几何结构，根据相邻像素点之间的亮度值不同，给每条边赋予一定的权重，权重公式如下式所示：

w_ab＝2-w_a-w_b

其中w_ab表示节点a与节点之间的权重，w_a表示节点a的亮度值，w_b表示节点b的亮度值，通过输入血管的起点(s)和终点(e)的位置，采用戴克斯特拉算法获取起点与终点之间血管的位置，初始时，起点(s)的路径长度为0，若存在能达到边(s，m)，则把起点(s)到节点(m)的距离设置d_s，m，不能到达则设置为无穷大，不断对外接变进行拓展，若存在节点(n)使得d_s，n+w_n，m＜d_s，m，则更新d_s，m，且设置节点(m)上一节点的为节点(n)。从而可以获得起点到终点之间的最短路径，即为需要探测的血管的位置；

(4)三维血管几何结构的重建：根据二维血管的几何特征，初步估计血管的位置，及中心线的位置，进而对根据二维血管的位置对三维体数据进行重采样，获取单根血管的三维体数据，通过对体数据进行展开，并应用基于扩展的动态规划的三维分层算法获取血管边界，最终获得眼底血管的三维形态结构数据，实现视网膜视乳头区的3-D血管几何结构的重建；

(5)获取OCT的多普勒图像：采用环形扫描或水平扫描方式获取OCT的多普勒图像；

(6)流速校正：根据步骤(3)获取的血管的几何边界数据，使用最小二乘法拟合血管边界，从而获得血管中心位置，可以获得沿着血管中心线方向需要计算的血管位置的方向，从而获得入射光线与血管之间的夹角，即多普勒角度θ，对步骤(4)获取的血流分量进行角度矫正：

$v=\frac{v_z}{cos\mathrm{\θ}}$ ，

其中v为绝对流速，v_z为沿着深度方向的流速分量，从而获取视网膜视乳头区血管特定位置的绝对流速。

2.根据权利要求1所述的一种在体视网膜血流动力学的成像与绝对流速测量方法，其特征在于，所述的步骤(2)中FD-OCT视网膜视乳头区图像纵向分辨率为3μm，采集速率不小于50kHz，三维血管造影数据由512×128×8条A扫描组成。

3.根据权利要求1所述的一种在体视网膜血流动力学的成像与绝对流速测量方法，其特征在于，所述的步骤(2)中三维结构数据获得的方法为：采用方块扫描方式获取三维结构数据，每个三维数据由128幅断层图组成，每幅断层图像有在同一个位置8次重复B扫描组成，然后采用基于强度信息的多普勒方差方法通过在同一个位置做重复扫描，利用相邻两次B-scan之间强度的变化，实现对微血管的成像。

4.根据权利要求3所述的一种在体视网膜血流动力学的成像与绝对流速测量方法，其特征在于，所述的基于强度信息的多普勒方差方法的算法如下：

${\mathrm{\σ}}^2=\frac{1}{T^2}\[1-\frac{{\mathrm{\Σ}}_{m=1}^M{\mathrm{\Σ}}_{n-1}^{N-1}\left(\right|A_{n,m}\left|\right|A_{n+1,m}\left|\right)}{{\mathrm{\Σ}}_{m=1}^M{\mathrm{\Σ}}_{n-1}^{N-1}\frac{1}{2}\left(\right|A_{n,m}{|}^2+\left|A_{n+1,m}{|}^2\right)}\]$

其中M代表平均扫描深度的像素点，N代表同一位置扫描的次数，A_n，m代表OCT数据的幅度值。

5.根据权利要求1所述的一种在体视网膜血流动力学的成像与绝对流速测量方法，其特征在于，所述的步骤(5)的OCT的多普勒图像由1024个像素点组成，采集图像的时间＞4秒，保证能测量在不少于4个心动周期内的视网膜血流的搏动，通过分析相邻两个A扫描直接的相位差，如下式所示：

其中为计算所得的相位差，Im(A_j，z)表示在第J条线，深度为z位置处虚部，Re(A_j，z)则为实部。J和N分别为横向A扫描和深度方向Z的平均数，沿着深度方向的流速分量v_z为：

其中λ_c为光源中心波长，n为视网膜的折射率，ΔT为扫描速度的倒数。