在体全眼微血流超声超分辨成像

张昌露 雷 爽 刘佳妹 商东庆 马 腾*

1（中国科学院深圳先进技术研究院 深圳 518055）2（中国科学院大学 北京 100049）

1 引 言

研究发现，多类眼科疾病都伴随有血管形态学改变[1]，如视网膜脱落、眼睛肿瘤[2]、青光眼[3]、糖尿病视网膜病变[4]和新生儿视网膜疾病[5]等。此外，眼底血管系统的相互作用在维持视神经的健康方面起着重要作用。目前，眼科疾病诊断常用的临床技术主要有光学相干断层扫描成像[6]、光学相干断层血管造影[7]、光声成像[8-9]、眼底镜、超声 B 模式成像、超声多普勒成像等。光学相干断层扫描成像、光学相干断层血管造影虽然是眼科常用的诊断技术，但是其成像深度较浅，若出现玻璃体混浊等其他影响屈光的情形，则光学成像可能无法准确探查眼底血管。除了眼睛微血流形态对维持眼睛功能至关重要外，眼睛血流动力学参数对眼部疾病的诊断也具有重要临床意义。例如：急性闭角型青光眼患者在不同眼压状态下，眼动脉、视网膜中央动脉和睫状后动脉的血流速度不同，可以根据血流速度对疾病进行分期分型。血流方向有助于了解血管网动静脉分布，可以指导眼科外科手术方案，以有效防止术中血管破裂，如眼部肿瘤及视网膜脱落手术评估等，因此，血流动力学检测对眼科疾病的诊治十分重要。受超声波衍射极限的制约，传统超声成像的成像分辨率较低，无法直接对眼底血管变化进行成像。超声多普勒成像虽然是目前常用的血管成像技术，但是其只能粗略显示血管系统，无法准确可视化微血管系统。因此，一种无创精确的眼睛微血管成像技术对眼科疾病的早期诊断和治疗评估具有重要意义。

超声成像虽然具有无辐射和性价比高等优势，但是传统超声聚集成像受成像分辨率和成像帧率的制约。超声平面波复合成像的提出打破了这一限制，将成像帧率提高至数千赫兹[10]。平面波复合成像方法的提出使得依赖于快速采集的超声功能成像成为可能。超分辨微血流成像是一种新的超声成像技术，它打破了超声波的衍射极限，在保证成像深度的同时可以可视化微米级血管[11]。与随机光学重建显微镜[12]和荧光激活定位显微镜[13]类似，超声超分辨成像主要对血管内孤立微泡进行识别、定位和追踪，形成微血管轨迹图，并估计血流速度和方向。在过去几年，国内外众多学者对超声超分辨成像算法和相关应用进行了研究，进一步推动了其临床应用进展。总的来说，主要通过改进微泡信号提取算法[14]、图像滤波算法[15-17]、微泡中心定位算法[18-20]和基于深度学习加速重建过程[21-22]。目前，超分辨微血流成像已经实现了多种器官组织的微血流标定，包括肿瘤[23]、淋巴[24]、脑[25-26]和眼底等重要组织，对于进一步研究组织器官的生理和病理血流变化具有重要意义。

临床眼科学一般采用频率为 10～20 MHz 的超声换能器，在这个频率范围内，穿透深度和成像分辨率均可达到眼底成像的要求，但传统的超声多普勒成像则无法观察到毛细血管[27]。以往研究表明，超声超分辨成像可在保持原有穿透深度的情况下，将空间分辨率提高约 10 倍[28]。因此，超声超分辨成像技术将为全眼微血管系统的可视化提供一个有力的工具[29-30]。一方面，利用二维超声超分辨成像对眼底微血流进行成像，虽然在不同的生物体组织中取得了较好的成像结果，但是二维超声超分辨成像具有高度的成像平面依赖性，若不能选取合适剖面，则可能会在临床应用中出现漏诊事故。另一方面，二维超声超分辨成像不能实现换能器阵元仰角方向的微泡追踪，会丢失仰角方向的血流信息。因此，若能进行三维空间内的微血管可视化，则对疾病的诊断和治疗评估具有重要意义。

目前，三维超声成像主要有以下方法：一是利用二维超声阵列实现三维超声成像[31-33]，这种成像方式的缺点是采集时间和帧率不足，且时间成本非常大；二是利用行列寻址进行超分辨成像[34]的信噪比较低；三是利用一维阵列通过机械扫描进行三维组织成像，通过调整合适的电机步长，实现对整个组织的三维成像[35]。

基于以上讨论，本文利用中心频率为 10 MHz的线阵超声换能器，基于机械扫描实现三维空间内全眼微血流成像，以可视化整个眼睛的微血管系统。与传统的超声图像相比，超声超分辨成像微血管图像的分辨率有了极大的提高，且可以实现全眼三维微血流成像。最终的三维超声超分辨成像结果证明本文所述方法能够可视化整个眼睛的微血管网络。

2 方法与材料

2.1 成像系统参数

在数据采集之前，需要对超声成像系统的点扩散函数（point spread function，PSF）进行标定，将直径 20 μm 的钨丝水平固定在水箱中，模拟点散射体，调整线阵换能器位置，确保成像平面垂直于导线。分别在 5 mm、10 mm 和 15 mm 深度下测试 PSF 随成像深度的变化。计算 150 帧图像的横向和轴向的平均半高宽，并将其用于微泡中心点定位。

2.2 实验装置

本研究的动物实验方案经中国科学院深圳先进技术研究院机构动物护理与使用发员会（IACUC）批准（SIAT-IRB-180213-YYS-MT-A0458）。所用的动物为健康新西兰白兔（约 2.5 kg），实验开始前用体积分数为 3% 异氟醚吸氧麻醉新西兰白兔，深度麻醉后将其置于暖垫上，维持体温37 ℃，然后将超声耦合剂凝胶均匀涂于换能器表面，使换能器与右眼表面充分耦合。在数据采集之前，每 1.5 min 耳静脉注射 0.3 mL SonoVue（Bracco Imaging SpA，意大利）微泡悬浮液，然后注射 1 mL 生理盐水冲洗。可编程超声成像平台 Verasonics Vantage 256（Kirkland，美国）和线性阵列换能器 L14-6（上海爱声生物医疗科技有限公司，中国）用于超声数据采集。超声换能器的中心频率为 10 MHz，带宽为 80%，每次传输一个信号周期，机械指数设置为 0.1。采用超快复合平面波成像，平面波偏转角度为－6°～6°等分的 5 个角度，有效脉冲重复频率（PRF）为 800 Hz。将线阵换能器沿着扫描方向安装在三轴位移平台 M8-400R（北京卓立汉光仪器有限公司，中国）上，步长为 150 μm，生成三维体图像。每个剖面共采集了 2 400 帧超声图像序列，共采集 50 个剖面，以重建三维体图像。在数据采集后，采用 delay and sum （DAS）波束合成方法将原始射频数据波束合成 IQ 数据，波束合成数据网格的大小设置为半波长。最后利用MATLAB 2019a（MathWorks，美国）保存 IQ 数据，进行后处理。

2.3 微泡信号提取

对于每个剖面的射频数据序列来说，首先应用波束合成后的 IQ 数据生成原始 B 模式超声图像，然后利用二维运动矫正算法去除因呼吸运动或者其他生理运动引起的伪影[36]。本研究采用分块时空奇异值分解（singular value decomposition，SVD）滤波器，从组织信号中提取微泡信号[16]，如图1 所示。连续采集的 IQ 序列是三维时空矩阵，空间维度（图像像素大小）为Nx和Nz，时间维度为Nt。在本研究中，每个剖面采集了2 400 帧 B 模式图像，以用于重建超分辨血流图像，每张图的大小为Nx＝316，Nz＝250。将整个数据集分割为若干个维度为Nx，Nz，Nt的子块，其中，nx≤Nx，nz≤Nz，nt≤Nt，且子块在横向和轴向方向重叠。对每个子块的 IQ 数据进行二维Casorati 矩阵重构，第一个维度表示空间，另一个维度表示时间。采用以下公式对该矩阵进行SVD：

图1 分块 SVD 滤波示意图Fig.1 Diagram of the block-wise SVD filter

其中，Sn为第n个子块的二维 Casorati 矩阵；U为Sn的空间的奇异值向量所张成的矩阵；V为Sn的时间的奇异值向量所张成的矩阵； 是奇异值矩阵，*代表矩阵的共轭转置。为了去除背景组织信号和抑制噪声，选择了两个阈值，其中，奇异值能量曲线二阶导数最大的点为低阶阈值，高阶阈值则由曲线斜率的平坦度决定。每个像素微泡信号 为：

其中，N为包含目标像素 的重叠子块的总数；Snb为第n个子块的流动微泡信号； 为剩余奇异值对应的微泡信号；Low-order 代表最小的奇异值截断阈值；High-order 代表最大的奇异值截断阈值。为了在深度上提高信噪比，减少计算时间，设定三维子块大小为 316×100×2 400，即nx＝316，nz＝100，nt＝2 400，此处，nt的选择是根据经验设定的，以分别提取眼前节和眼底微血管。其中，沿Z方向的相邻子块之间的距离（步长）为定值，如果像素网格尺寸为波长，则步长设定为 3；如果像素网格尺寸为半波长，则步长设定为 5。

2.4 微泡中心定位

在进行时空滤波后得到了微泡信号，采用三维时空块匹配去噪滤波器[15]，这样可以在保持微泡图像特征的同时抑制噪声，然后对每个像素的归一化强度应用硬阈值，去除信噪比低的微泡信号。采用二维样条插值方法对每帧流动微气泡进行上采样，像素尺寸为 8 μm×8 μm，而单个微泡的半径范围为 1～3 μm，远小于超声波半波长，因此可以认为是一个散射点。对于微泡定位，首先将校准后的平均 PSF 拟合为二维高斯核模型，然后与插值后的微泡信号进行二维归一化互相关，具体数学表达如式（3）所示，然后将相关系数低于 0.6 的像素点置为零，以识别潜在的微气泡区域，最后计算互相关稀疏的局部极大值，即可看作为微气泡中心点。

2.5 微泡追踪

在 2.4 小节中，完成了对微泡信号的定位。为了生成血流的速度和方向图，需要对每一帧的微泡进行追踪，本文利用库恩-曼克尔斯分配算法实现微泡的多帧追踪[37]。首先，计算所有微泡在连续两帧之间的距离，然后，根据总距离最小化的原则对连续两帧的微气泡进行配对，并剔除总长度小于 5 帧的轨迹，进一步抑制噪声的影响。血流速度计算是将平均距离除以配对微泡之间的时间间隔，然后将各帧中成对微泡的位置连接起来，生成血流速度幅值和方向图。

3 实验结果

3.1 成像系统参数测量

按照 2.1 小节中的方法对成像系统的点扩散函数进行测量，测量结果如图2 所示。图2（a）举例说明了深度为 5 mm、10 mm、15 mm 处的钨丝导线横截面的 B 模式超声图像，各个散射点的PSF 的横向和轴向强度分布如图2（b）～（c）所示。平均横向分辨率和轴向分辨率分别为 251 μm 和155 μm，理论计算的超分辨成像的定位精度约为8 μm[11]，可用于后续微泡中心点定位处理。

图2 超声成像系统的校准——配备 L14-6 探头的 Verasonics 256 研究平台Fig.2 Calibration of ultrasound imaging system—Verasonics 256 research platform equipped with an L14-6 probe

3.2 二维超声超分辨微血流成像结果

按照 2.2 小节所述的数据采集方式，图3 展示了在瞳孔正中间位置采集的超声成像平面超分辨血流图像。图3（a）是超声成像平面的示意图，成像深度约为 22 mm。在图3（b）的 B 型图像中，可以清楚地看到眼前段的角膜、晶状体和虹膜。如图3（c）所示，应用 SVD 滤波将静止组织和背景噪声信号进行滤除，经过多帧叠加后得到超声多普勒血管图像。图3（d）为重建后的超分辨微血管图像，成像剖面中眼眶和虹膜的微血管清晰可见。图3（e）表示微泡流在血管中的灌注情况，不同的颜色编码分别表示不同的血流速度。在图3（f）中，不同的颜色编码分别表示不同的血流方向。

图3 第一个成像剖面Fig.3 An image section of the rabbit eye

图4 为兔眼的另一成像剖面，在这个剖面内，可见巩膜和球后微血管，眼前血管沿巩膜呈圆形分布，并且可以观察到，视网膜和脉络膜有密集的脉管系统，血管的数量和密度沿视网膜到脉络膜的深度增多。血流速度和血流方向如图4（e）～（f）所示，从图4（e）中可以看出，超分辨成像可以检测到低至 4.1 mm/s 的微血流。

图4 另一成像剖面Fig.4 Image of another section of the rabbit eye

图5（a）为图4（c）中白色框内的放大视图，图5（b）为图4（d）中白色框内的放大视图。从图5（a）～（b）中可以看出，微血管的细小分支不能通过功率多普勒超声显示，但超声超分辨成像却可以可视化微血管。局部放大超分辨图像标记位置的截面轮廓如图5（c）所示，可以看出，超分辨成像可以清楚地区分两个相距 61 μm 的微血管，并且可以可视化直径小于 50 μm 的单个微血管。

图5 多普勒与超分辨成像分辨率对比分析Fig.5 The contrast of Doppler and super-resolution imaging

3.3 二维超声超分辨微血流成像结果

3.2 小节展示了两只眼睛剖面的超分辨图像，利用机械扫描方法一共采集了 50 个剖面，本小节将 50 个剖面重建为三维微血流图像。图6（a）为三维功率多普勒图像，图6（c）为超分辨成像微血管图像。图6（b）为与原 B 模式图像融合后的多普勒血图，图6（d）为与原 B 模式图像融合后的微血管图像。由图6（a）～（d）可知，眼底有致密的微血管，且视网膜和脉络膜的血管包裹眼球，这与已知眼睛的解剖结构是一致的。

图6 兔全眼机械扫描三维微血流图像Fig.6 The 3D volume images of the whole eye of the rabbit

4 讨 论

本文应用超声超分辨微血流成像技术对全眼微血管进行了在体研究。实验结果表明，距离为 61 μm 的两个微血管可以被清晰区分，成像算法可以可视化直径小于 50 μm 的单个微血管，其中，理论分辨率极限约为 96 μm[11]。总的来说，在超过 15 mm 的成像深度下，三维超声超分辨成像成功地实现了对眼睛后部微血管的检测，分辨率约为 50 μm。

眼睛血液微循环的改变是各种眼部疾病的早期征兆，如青光眼、黄斑变性、早产儿视网膜病变和糖尿病视网膜病变等。光学相干断层血管造影、眼底镜等光学成像技术是临床眼科诊断中常用的高分辨率眼血管可视化技术[7-8]。然而，由于其穿透深度较浅，不能有效穿透不透明组织，因此，这些光学成像方式很难可视化球后血管[9]。超声成像具有无创、不受屈光介质的影响等优点，已成为眼科检查不可缺少的工具。然而，由于衍射的限制，传统的超声多普勒成像的分辨率被限制在半波长以下，很难检测到几十微米或更小的微血管[28]。因此，本研究采用基于超快复合平面波成像技术的超分辨成像对全眼微血管进行可视化，对眼科疾病的诊断具有重要价值。与已有超声超分辨眼底微血流成像相比[29-30]，本文所提出的方法可以对全眼微血流系统进行成像，并且可以通过机械扫描对整个眼部区域的微血流进行三维成像，可以进一步提高眼科疾病的早期诊断效率。

然而，本研究存在一定的局限性。

第一个局限性是成像技术采用线阵换能器，无法直接在三维空间中定位微气泡，因此，在阵元仰角方向上的分辨率受到了限制。此外，由于换能器通过机械扫描重建三维超声图像，因此无法消除机械扫描方向上的运动伪影。这是使用一维线性阵列换能器进行三维超分辨成像所存在的一个缺点。在未来的研究中，使用二维矩形阵列换能器或其他方法可能会解决这个问题[38]。

第二个局限性是超分辨成像所需要的采集时间比较长，特别是对于血流速度较低的小血管而言，超分辨率微血管图的重建需要多帧微泡中心点的累积。在本研究中，每个扫描平面在 800 Hz的帧率下获得 2 400 帧。因此，每个切片的成像大约需要 8 s（数据采集 3 s，数据保存 5 s），在机械扫描步长为 150 μm 的情况下，对整个眼睛进行完整的 3D 扫描大约需要 7 min 的成像时间。然而，由于眼睛的高敏感性，在临床眼科诊断应用中，获得大量数据具有挑战性。为了减少采集时间，可以在未来的工作中使用基于深度学习的定位方法[21-22]来定位高密度的微气泡。

第三个局限性是微泡检测方法。本文利用预先测定的平均 PSF 与提取的微泡信号进行二维归一化互相关。然而，不同半径、形状的流动微泡的 PSF 可能不同，如果使用高互相关阈值，则可能会将微泡信号滤除，不利于超分辨图像重建。如果使用低阈值，则可能将其余噪声信号当作微泡信号，从而引入误差信息。由于无法准确估计微气泡的 PSF，因此很难在微气泡对应的信号水平和噪声信号水平之间设定一个精确的折衷。在后面的研究中，需要设计更加精确、快速的微泡定位算法，以重建超分辨图像。

5 结 论

综上所述，本研究成功利用三维超声超分辨成像显示了整个眼睛的三维微血管系统。与传统的光学成像和超声多普勒图像相比，超声超分辨成像可以区分两个直径小于 50 μm 且距离为61 μm 的血管，成像分辨率远高于传统多普勒成像，且成像深度也远高于光学成像，除了可以精确提供血管的形态学参数外，还可以得到血管的动力学参数。体内实验结果表明，三维超声超分辨成像在眼部疾病的早期诊断方面具有临床应用潜力。