基于局部信噪比阈值调节的分光谱去相关光学相干断层扫描血管造影方法

在基于强度信息的光学相干断层扫描血管造影（OCTA）技术中，去相关映射被广泛应用，但该方法易受噪声影响，尤其是对于低信噪比区域，噪声会使帧间高相关的静态组织信号呈现出高去相关性，与血流信号的高去相关值混叠，从而会降低微血管图像的成像质量。本文提出了一种基于局部信噪比动态阈值调节的分光谱振幅去相关方法，首先探讨了图像局部信噪比与静态体素去相关值的关系及影响因素，并以分频谱处理数据的方式削弱轴向运动伪影并加强血流信息，然后依据图像局部信噪比与静态体素去相关的关系调整设定阈值，用以判定区分动静态体素，阈值以下的去相关值视为静态体素，结合Sigmoid函数予以抑制。采用该算法在高信噪比的皮肤数据和低信噪比的视网膜数据处理中均能取得较好的效果：皮肤enface图像信噪比提升约4 dB；视网膜enface图像的对比度相较其他算法优势明显，相比于未进行静态抑制的分频谱振幅去相关造影（SSADA）方法，enface图像的对比度提升了0.0182。可见该方法能够有效地抑制微血管图像中的静态体素噪声，提升图像对比度和血管可视性，有利于疾病的诊断和评估。

Abstract

Objective

In optical coherence tomography angiography (OCTA), the applications of decorrelation mapping, primarily reliant on intensity data, have caught significant attention. However, this method is particularly vulnerable to the deleterious effects of noise, especially in fields characterized by low signal-to-noise ratios (SNRs). Noise artifacts have a pronounced effect on static tissue signals, which makes them exhibit elevated decorrelation between frames and in turn tends to overlap with the high decorrelation values associated with blood flow signals. This overlap detrimentally affects the quality of microvascular image acquisition. Meanwhile, classical techniques for refining decorrelation mapping, such as frequency-domain decorrelation angiography, still struggle to yield optimal results due to this inherent challenge. In response to the spurious static voxel artifacts, some studies have resorted to employing thresholding to eliminate static voxels falling below a predefined threshold. However, the global and indiscriminate nature of such thresholding often lacks a robust theoretical foundation, making the precise suppression of static voxel artifacts a complex endeavor. To this end, we present a novel OCTA approach that incorporates considerations of SNR and dynamic threshold adjustments. This innovative method is further combined with spectral analysis principles to provide a more precise means for the identification and suppression of static voxels. The ultimate objective is to enhance the microvascular imaging quality, thereby serving as a more dependable foundation for medical diagnostics.

Methods

We introduce a method for spectral amplitude decorrelation, which features dynamic threshold adjustments based on local SNRs. The methodology commences with an in-depth exploration of the complex relationship between local image SNRs and static voxels, including a comprehensive analysis of the various factors influencing this association. Subsequently, spectral analysis techniques are employed to mitigate artifacts arising from axial motion and accentuate the visualization of blood flow data. Built upon the established connection between local image SNRs and static voxels, our approach proposes adaptive thresholds for each voxel to ensure precise differentiation between dynamic and static voxels. Voxels exhibiting decorrelation values below the established threshold are categorized as static ones and subsequently suppressed. Conversely, voxels surpassing the threshold are identified as dynamic ones and are retained. Meanwhile, we further employ a sigmoid function to apply non-linear mapping to all voxels, thereby facilitating a seamless transition at the boundary between dynamic and static voxels. After the suppression of static voxels, an averaging process is applied to the decorrelation images, which allows us to reconstruct enface microvascular images by the mean projection technique. Additionally, we have established a dedicated posterior segment SS-OCT system to collect retinal data from volunteers. The effectiveness of our algorithm is rigorously validated via the data, and we conduct comparative experiments with other classical intensity-based OCTA methods to comprehensively assess its performance.

Results and Discussions

In comparison to the conventional decorrelation mapping approach, the retinal blood flow cross-sectional images processed by our algorithm exhibit prominent blood flow signals, whereas the conventional method's results are largely submerged within the noise emanating from static tissue (Fig. 6). This disparity highlights that the SSADA algorithm affected by noise-induced interference in individual spectral amplitude decorrelation images produces lower-quality enface microvascular images after averaging. In contrast, our algorithm effectively suppresses the noise arising from static voxels within individual spectral amplitude decorrelation images, ultimately yielding high-quality enface microvascular images. Compared to other intensity-based OCTA techniques, our proposed algorithm demonstrates superior performance across both high SNR skin data and low SNR retinal data, with the same preprocessing, target extraction, and image registration protocols employed. For skin data, the enface microvascular images obtained by our algorithm exhibit an SNR enhancement of approximately 4 dB in contrast to the SSADA method without static voxel suppression (Fig. 5). In the case of retinal data, our algorithm produces enface microvascular images with significantly improved contrast ratio, achieving a contrast enhancement of 0.0182 compared to the SSADA method without static suppression (Table 1).

Conclusions

We conduct a systematic examination of the intricate relationship between local SNRs and the decorrelation values of static voxels in OCT structural images. The results show that as noise levels on voxels increase, static voxels exhibit higher decorrelation values. Based on this pivotal finding, we introduce a dynamic threshold adjustment method within the context of spectral analysis. This combined approach adeptly leverages the sensitivity of decorrelation mapping to subtle differences and the efficacy of spectral analysis in mitigating artifacts stemming from axial motion. The retinal enface microvascular images produced by our algorithm adeptly differentiate capillaries in proximity to the macular region, underscoring the algorithm's competence in effectively suppressing static voxel noise within microvascular images. Furthermore, our algorithm consistently delivers favorable outcomes in retinal data characterized by low SNRs, resulting in enhanced image contrast ratio and superior vessel visibility. This enhancement has great potential in improving disease diagnosis and evaluation, contributing to more precise medical assessments.

1　引言

光学相干层析（OCT）成像基于光的低相干干涉原理，通过分析样品不同深度和不同横向位置的后向散射信息^［1］，能够对包括生物组织在内的强散射介质进行非侵入式成像^［2］。光学相干断层扫描血管造影（OCTA）作为一种OCT的功能成像技术^［3-4］，在近十年得到了快速发展。该技术通过比较不同时间重复多次采集的体扫描样品信息，利用红细胞的动态散射特性，检测信号的变化来提取血流信号，从而实现微血管的无创可视化，被广泛应用于医疗领域^［5-10］。

在OCTA领域，根据利用的信息不同而分为强度、相位、复合三类。其中，基于相位的多普勒方差^［11］、多普勒相位分辨^［12］等方法利用多普勒效应将OCT干涉相位与血液流速建立联系，能够对血管进行造影，基于复合信息的超高灵敏度光学微血管造影（UHS-OMAG）技术^［13-14］同时利用了信号相位和振幅信息做差异提取，在血管造影上呈现出较好的效果，但目前被高度关注的第三代扫频OCT（SS-OCT）系统的相位并不稳定，往往不能采用相位和复合OCTA方法。相比之下基于强度信息的OCTA无须引入复杂的光学元件和相位补偿方法^［15-16］进行相位校准，成像系统更加简单，而且只分析样品的振幅信息使其具有更快的成像速度，更适用于目前主流第三代SS-OCT系统实时成像的应用要求。

相关法（CM）作为最常用的强度OCTA方法，相较于差分法^［17-19］（DM）对体素自身灰度值大小做了归一化处理，因而在小信号区域提高了对微小差异的敏感性，在大信号区域削弱差值敏感性，将体素自身灰度值纳入考量，可以丰富血管细节信息以及抑制生物组织轮廓边缘造成的伪影。干涉测量领域常见的相关系数计算方法为Pearson相关^［20］和Berger相关^［21］两种形式，实际应用中受信噪比的影响较大，例如，在低信噪比区域，静态体素之上会附加大量噪声信息，而使其具有了较低的相关值（较高的去相关值），其性质将更倾向于动态体素，若与真正的动态信息（血流信号）混叠，就会降低血流截面以及正面微血管（enface）图像的质量。分频谱振幅去相关造影（SSADA）方法^［22］，是在相关法的基础上提出的一种经典的分频谱处理的方法，通过降低轴向分辨率以削弱轴向运动伪影，能够提高血管造影图像的质量。但对于低信噪比数据，该方法仍然具有相关法的缺陷，且当单个频谱因噪声导致血流信号与静态信号混淆在一起时，多个频谱叠加平均的结果会更差。

针对上述问题，在不考虑增加重复采集次数的情况下，研究人员提出阈值相关OCTA^［23］，设定一个确定的去相关阈值来排除噪声对静态组织的影响。通常以0.5作为全局阈值，只有超过0.5的去相关值才被认为是血流信号，而低于阈值的去相关值则被视为噪声或静态组织。但该阈值的选择对结果影响较大且缺乏依据，尤其样品的结构图像因各处介质的光学特性不同而具有明暗变化的灰度值，一个全局的阈值无法准确地处理不同灰度区域的静态体素，过高或过低的阈值可能导致血流信号的丢失或噪声干扰。

本文提出了一种结合体素局部信噪比的强度OCTA方法以实现信噪比自适应的血流提取。首先计算OCT图像体素的局部信噪比与帧间去相关值，为不同局部信噪比区域施加适当大小的动静态体素分类阈值，利用Sigmoid函数映射抑制静态体素，保留动态体素的同时过渡处理阈值附近的动静态体素。同时，还采用了分频谱处理策略，通过频谱分割来降低轴向分辨率，从而减弱轴向运动伪影的影响。该方法提升了enface图像的质量，并且在信噪比较低的视网膜数据中也能取得显著的效果，相比于未进行静态抑制的SSADA方法，enface图像的对比度提升了0.0182，具有更好的血管可视性。

2　基于局部信噪比的阈值动态调整方法

2.1　基本原理

Berger相关作为一种最常用的相关方法，是以某一体素在相邻时间的归一化振幅之差作为样本，来估计该体素在时间序列上的波动模型^［24］。以其为代表的相关映射类方法考虑了被比较的像素点的灰度值大小，相比于仅考虑差值绝对大小的方法，具有将体素差异与体素自身尺度关联从而动态调整对差值的敏感性的优势，其血流计算公式形式为

\begin{matrix} σ^{2} = \frac{1}{N - 1} \sum_{i = 1}^{N - 1} \frac{{(a_{i} - a_{i + 1})}^{2}}{a_{i}^{2} + a_{i + 1}^{2}} \end{matrix},

（1）

式中： $a_{i}$ 为同一横向位置第i帧图像中体素的灰度值；N为重复采集次数。该式的分母是对信号本身灰度值进行的归一化，在无噪声的理想状态下，归一化放大了低反射率区域的弱信号波动，能提高血管造影的灵敏度。在理想无噪声情况下，静态体素动态程度为零，计算出的去相关值为0。但实际上分子包含了动态程度与帧间噪声差两部分，需要考虑OCT结构图像中由光源和光电探测器引起的随机噪声^［25］。当静态信号自身强度远大于噪声尺度时，仍可以得到较低的去相关值，噪声的影响很小；但如果静态信号强度逐渐降低，甚至趋近于噪声信号时，静态体素的去相关值会越来越高，在血管图像中形成静态伪影，降低微血管图像的质量。

如果想要削弱静态伪影，就需要一个合适的阈值来区分静态与动态信号。可以看出，静态体素的去相关值与体素信噪比呈负相关，但OCT结构图像的信号强度与样品各处折射率相关，视网膜也因各层结构的介质折射率不同而呈现出明暗相间的多层结构，无法通过一个确定的去相关阈值来抑制整个图像的静态体素。因此可对OCT结构图像窗口化，结合体素的局部信噪比设置自适应变化的动静态分类阈值，为高信噪比区域设置一个较低的阈值，为低信噪比区域设置一个较高的阈值。将去相关在阈值以下的体素判定为静态体素，予以抑制；将去相关值高于阈值的体素判定为动态体素并保留其去相关值。

OCT在轴向上较高的分辨能力使其对轴向位移非常敏感，因而容易受到组织运动的影响而产生伪影，降低enface图像的质量。SSADA将光谱在轴向上分为几个窄带，计算窄带间的振幅去相关并叠加平均，以降低对轴向位移的敏感度，提高血管图像质量。但SSADA采用了Berger相关系数作为血流对比机制，也存在着低信噪比区域受噪声影响大的问题。分频谱之间计算得到的血流噪声淹没在静态伪影里，分频谱叠加平均后的结果甚至比分频谱的结果更差，导致发挥不出分频谱处理的作用。因此本文将提出的结合局部信噪比自适应分配阈值的方法应用在Berger相关系数上，结合分频谱处理，进一步提高enface图像的质量。

2.2　计算公式

本文将参与Berger去相关计算的空间核像素点的范围扩展到以待计算像素点为中心的矩形窗口范围内的所有像素点，该局部窗口的设置即配合了前述动静态体素分类的策略，将多帧B-scan中体素局部窗口内（时空核）的点看成在时间上具有同分布的体素，以窗口内体素的信噪比和计算所得的去相关值来判定中心像素点的分类结果。经过扩展并结合了分频谱思想的体素Berger去相关系数（D）以及对应的局部信噪比（S_NR）通过下面两个式子定义：

\begin{matrix} D = 1 - \frac{1}{M} \frac{1}{N - 1} \frac{1}{(2 P + 1) (2 Q + 1)} \sum_{m = 1}^{M} \sum_{n = 1}^{N - 1} \sum_{p = - P}^{P} \sum_{q = - Q}^{Q} \frac{A_{m, n} (x + p, z + q) A_{m, n + 1} (x + p, z + q)}{\frac{1}{2} A_{m, n} {(x + p, z + q)}^{2} + \frac{1}{2} A_{m, n + 1} {(x + p, z + q)}^{2}}, \end{matrix}

（2）

\begin{matrix} S_{N R} = 10 l g \frac{1}{M N (2 P + 1) (2 Q + 1)} \sum_{m = 1}^{M} \sum_{n = 1}^{N} \sum_{p = 1}^{P} \sum_{q = 1}^{Q} \frac{\frac{1}{2} A_{n} {(x + p, z + q)}^{2} + \frac{1}{2} A_{n + 1} {(x + p, z + q)}^{2}}{s^{2}}, \end{matrix}

（3）

式中：M表示分频谱个数；N表示重复采集次数；p为空间窗口在横向方向的索引；q为窗口深度方向的索引； $A_{m, n} (x + p, z + q)$ 表示第m个分光谱、第n个B-scan横向位置x+p、深度位置z+q处的振幅； $(2 P + 1) (2 Q + 1)$ 定义了空间核的大小； $s^{2}$ 为每一组分频谱下高斯白噪声的方差，对于SS-OCT来说，光谱干涉信号傅里叶逆变换后，结构图像的背景噪声信号 $n (x, z)$ 服从高斯白噪声的分布^［26］， $n (x, z)$ 的实部与虚部分别呈独立、均值为零、有相同方差的高斯分布。则噪声的方差即为噪声振幅平方的期望，在算法中定义为每个频谱段的多帧B-scan中背景区域像素灰度平方的均值，计算公式为

\begin{matrix} s^{2} = \frac{1}{N_{b a c k g r o u n d}} \sum_{(x, z) ϵ b g} \frac{1}{2} A_{n} {(x, z)}^{2} + \frac{1}{2} A_{n + 1} {(x, z)}^{2} \end{matrix},

（4）

按照式（2）计算得到的去相关值范围在0~1之间，其值越大，去相关度越高。以每个分频谱的背景噪声为基准，衡量每个局部窗口的灰度大小，为高灰度值的区域设立高阈值来区分动静态体素。

2.3　模拟分析

静、动态体素适应于不同局部信噪比的分类阈值通过模拟实验得到。实验采用原始光谱傅里叶逆变换之后的复信号进行，设置振幅范围为0~1的复信号模拟静态信号，并将其等分为若干份，在每一个小振幅范围内生成一定数量随机的静态复信号，这些幅值近似的信号视为同一个时空核下有着相同分布的时间序列，这样每一个小范围内的复信号数量即为模拟的时空核大小。背景噪声设置为某一振幅的高斯白噪声并使最高信噪比在25 dB左右，与现有眼底OCT结构图像的信噪比水平相当。

之后按照2.2节中式（2）~（4）计算每个时空核的去相关值D和局部信噪比（SNR）并映射在同一个坐标系下，模拟得到的静态体素分布结果如图1所示。图1（a）给出了时空核大小为25的结果。其散点分布趋近于一条曲线，符合之前提到的，随着局部信噪比的降低，静态体素的去相关值有增大的趋势。提取散点的右边界线并拟合，不同时空核的拟合结果如图1（b）所示，时空核越大，参与去相关计算的体素越多，相当于用越多的样本来估计一个模型，估计结果会越趋近于理论值，散点分布更紧凑，右边界线更陡峭。右边界线即为之前得到的结合局部信噪比的自适应去相关阈值，散点分布范围内的体素可以认为是静态或噪声体素，予以抑制。每个图片都分别计算背景噪声以及每个体素的局部信噪比，其分布理论上都符合仿真实验的结果。

图 1. 静态体素模拟分布。（a）时空核为25的静态体素分布；（b）时空核分别为10、20、25、30、40、50的分类曲线拟合结果。实线为最终选取的N=25的曲线

Fig. 1. Static voxel simulation distribution. (a) Static voxel distribution of spatiotemporal kernel of 25; (b) fitting results of classification curves with spatiotemporal kernels of 10, 20, 25, 30, 40, and 50, respectively. Solid line represents final selected curve with N=25

下载图片查看所有图片

分频谱数和窗口大小的选择要与采样品结合考虑，分频谱数越多，轴向分辨率降低的倍数越高，对轴向位移的敏感性越低，理论上可以更大程度地去除轴向运动伪影。但单个分频谱中的信息也会越少，且分辨元内可能会横跨多个红细胞与组织细胞，使去相关的计算失去意义。窗口（空间核）的选择不能过小也不能过大，窗口越大，参与计算的体素个数越多，静态体素的分布更趋近于理论值。但当窗口过大时，无法保证窗口内的体素均为同种分布，可能是由静态体素和动态体素混合组成，也可能是不同组织区域的静态体素的混合，实际和理论会产生偏差，导致分类的不准确。因此根据红细胞的大小（6~15 $μ m$ ）、眼底血管直径（几十微米）以及OCT系统的纵向分辨率（4~6 $μ m$ ），分频谱数通常设置为4~8，窗口大小通常设置为 $5 \times 5$ 比较合适。

如上所述，最终选取时空核为25的阈值曲线（实线），根据计算出的每个体素的局部信噪比确定该体素的动静态分类阈值，将去相关值低于阈值的体素判定为静态体素并将其去相关值抑制；将去相关值高于阈值的体素判定为动态体素，保留其去相关值。

2.4　通过Sigmoid函数映射过渡处理阈值附近的体素

得到动静态体素的分类曲线后，同等信噪比下，去相关值D高于和低于曲线的体素分别被认为是静态和动态体素，应对静态体素进行抑制处理。当信噪比降低到0时，说明信号强度与噪声尺度相当，理论上也已经获取不到有效信号，因此信噪比小于0的区域可认为是纯噪声。信噪比较高的静态体素受噪声影响小，去相关值较低，对微血管图像干扰较小，无须再抑制。分类曲线在去相关值大于0.1的部分可以用指数函数 $y = p e x p (q x)$ 进行拟合。

然而，动静态体素分界线附近会出现动静态体素交叠的现象，交叠的多少与样品的种类、局部窗口的大小、图像整体信噪比等因素有关，如果简单地将其两侧体素以0和1进行阈值分割，则很可能损失一部分血流信息，并误判一部分静态体素为动态血流信号。为了克服阈值分割法非黑即白的缺陷，算法利用Sigmoid函数为每一个体素生成对应的模板因子 $β$ ，从而对分界线附近的像素点进行过渡处理，它是一种常见的S型函数^［27］，相比于阶跃函数更加平滑地将变量映射到区间（0，1），常用于二分类问题。将分界线置于Sigmoid函数自变量为0的位置，则分界线上的体素最终的去相关值为0.5，意味着难以倾向静、动态体素中的任意一类，使映射更加合理。Sigmoid函数定义为

S (x) = \frac{1}{1 + e^{- x}}

。（5）

为利用Sigmoid函数在动静态分界线处起到平滑的作用，可逐行以分界线的D值作为原点构建Sigmoid函数。因此结合式（5），关注某个在I-D空间下坐标为（d，s_nr）的体素，以该点到分界线的横向距离作为Sigmoid自变量即可得到如图2所示的模板。图中黄色区域为动态体素，模板因子取1；蓝色区域为静态体素，模板因子取0；分界线附近体素的模板因子按照Sigmoid函数的趋势从0过渡至1。体素点（d，s_nr）在D坐标方向上与分界线的距离可通过下式计算：

图 2. 模板在I-D空间的示意图

Fig. 2. Schematic diagram of template in I-D space

下载图片查看所有图片

Δ D = d - \frac{l n (\frac{s_{n r}}{p})}{q}

。（6）

将 $Δ D$ 作为自变量代入式（5），像素点（d，s_nr）处的模板因子 $β$ 即为

\begin{matrix} β = S (Δ D) = \frac{1}{1 + e x p [\frac{l n (\frac{s_{n r}}{p})}{q} - d]}, \end{matrix}

（7）

式中：s_nr和d分别表示体素的局部信噪比和去相关值；p、q为拟合边界线参数。

2.5　算法流程

算法整体流程如图3所示，主要包含目标提取、分频谱、计算局部信噪比与去相关值、静态体素抑制、enface生成5个部分。

图 3. 算法流程图

Fig. 3. Algorithm flow chart

下载图片查看所有图片

1）叠加平均同一横向位置连续采集的N张B-scan，提取目标区域。

2）采用主瓣宽度尽量窄、旁瓣小、能量集中于主瓣的汉明窗作为窗函数将原始光谱分割为M个分光谱，并使相邻两个频谱有部分重叠，以避免信号截断导致的频谱泄漏。

3）计算各分频谱的背景强度平方的均值，之后计算 $M \times N$ 个结构图像的局部Berger去相关系数与局部信噪比。

4）计算对应分类曲线的Sigmoid平滑模板，将步骤3）中计算出的去相关值D与信噪比映射到S-D空间中，为每个像素点附加模板因子 $β$ ，以实现对动静态体素更准确的分类。

5）对逐个横向位置的N张B-scan进行步骤1）~4）的处理，并将处理结果沿轴向均值投影，生成最终的enface图像。

3　实验与分析

3.1　系统搭建

本文OCTA数据是在如图4所示的眼后节SS-OCT系统下采集的。系统采用1060 nm±50 nm、扫频速度为200 kHz的扫频激光光源（AXSUN），光源发出的光经9∶1耦合器分别进入样品臂和参考臂。样品光被经过准直器、二维扫描振镜（GVSK2-EC，Thorlabs）和光学透镜组后射入人眼并聚焦在视网膜上，人眼返回光从9∶1耦合器另一端口输出并进入50∶50耦合器；参考光束样品臂中的光经过发射准直器、色散补偿玻璃块后直接被接收准直器接收，在偏振控制器中调整偏振态后也进入50∶50耦合器中。参考光与样品光在50∶50耦合器中发生干涉，干涉信号被平衡探测器（PDB481C-AC，Thorlabs）接收并转化为电信号，干涉电信号由采样率为1 GS/s的高速采集卡（ATS9371，AlazarTech）采集并传入计算机处理。

图 4. SS-OCT系统结构示意图

Fig. 4. Structure diagram of SS-OCT system

下载图片查看所有图片

光源输出A-trigger信号触发采集卡并输出K-clock信号作为采集卡采样的外部时钟信号。采集卡Aux I/O口的输出信号触发波形发生卡（PCIe-6321，NI），使其输出不同的波形信号驱动二维扫描振镜，实现断层扫描和体扫描以采集OCTA数据。

3.2　实验结果

3.2.1　皮肤血管造影

皮肤组织相对于视网膜后向散射信号更强，采集到的OCT结构图像信噪比更高，噪声的影响更小，所以相关映射类方法能够得到一个较好的结果，图5（a）、5（b）、5（d）中的3种方法属于相关映射类方法，图像质量明显高于采用差分法的结果［图5（c）］。图5（c′）与5（d′）分别是图5（c）与5（d）中两种方法某一帧B-scan提取的血流图对比，通过观察可以发现，在上表面位置（虚线框），差分法残留了一些生物组织的轮廓，导致最终投影重建的微血管图像质量较差，而相关法在边缘部分提取得更加准确。这是因为样品浅层的信号强度更高，OCT结构图像在边缘轮廓位置灰度值较高，对于同样的差值，相关法考虑了体素自身的灰度大小，从而抑制了生物组织轮廓的伪影。

图 5. 皮肤数据成像结果。（a）提出算法的结果；（b）SSADA算法的结果；（c）SVOCT算法的结果；（d）Pearson相关法的结果；（c′）SVOCT算法的血流截面图；（d′）Pearson相关法的血流截面图。（a）~（d）中方框为选取的静态信号区域

Fig. 5. Skin data imaging results. (a) Results of proposed algorithm; (b) results of SSADA algorithm; (c) results of SVOCT algorithm; (d) results of Pearson correlation method; (c′) blood flow cross-section of SVOCT algorithm; (d′) blood flow cross-section of Pearson correlation method. Boxes in (a)-(d) represent selected static signal areas

下载图片查看所有图片

图5的结果显示，在高信噪比的皮肤数据中，本文算法同样具有和相关法同水平的成像效果，并在此基础上对静态体素进行了抑制。为了定量比较本文算法的性能，将图5（a）~（d）中的大血管标记为动态信号，虚线框区域标记为静态信号，计算了enface图中动态信号和静态组织的信噪比（S_NR）：

S_{N R} = 20 l g (\frac{I_{d y n a m i c}}{σ_{s t a t i c}}),

（8）

式中： $I_{d y n a m i c}$ 为动态信号平均值； $σ_{s t a t i c}$ 为静态信号标准差。由于实验无法获得微血管的真实图像，信噪比的结果与主观选择的静态区域有关，质量评估会带有一定的主观性，但也能在一定程度上体现图像的质量。在以正方形框区域作为静态区域时，本文算法与Pearson相关映射和散斑方差（SVOCT）算法相比具有最高的信噪比，相比没有静态体素抑制的SSADA算法，信噪比从15.01 dB提升到了19.83 dB，能得到更高质量的enface结果。

3.2.2　视网膜血管造影

算法的核心思想是考虑了体素局部信噪比分情况处理而非全局设定阈值，因此数据的信噪比越低，本文算法的提升效果越明显。对于散射信号微弱、信噪比较低的视网膜数据来说，其受噪声影响很大，相关法效果往往比差分类方法差，而采用本文算法能够在摒除噪声信号影响的情况下发挥出相关法对微小差异敏感的优势，显著提升enface图像的质量。本文对5位志愿者的视网膜进行了500 Ascan×2×500 Bscan的体扫描，扫描范围为黄斑区附近3 mm×3 mm的方形区域，每个断层图像尺寸为500 pixel×500 pixel，逐个横向位置重复扫描两次，采集时间间隔为2.5 ms。

实验中视网膜区域下边界选取为色素上皮层（RPE）约5个像素处，根据工程经验，分频谱个数设置为4，局部窗口大小设置为5×5。图6对比了一个B-scan在同样采取分频谱思想的情况下，有静态体素抑制和无静态体素抑制的血流截面图。在低信噪比的眼底图像中，直接进行去相关计算得到的除了血流信息还有大量噪声与静态体素相叠加引起的高去相关值。导致血流信号淹没在静态伪影中，如图6（a）所示。采取本文算法对静态体素进行抑制，能清晰地看到血管的分布，如图6（b）中箭头所标识之处。

图 6. 进行静态抑制前后的血流截面图。（a）未进行静态体素抑制；（b）静态体素抑制后；（c）OCT原始结构图像

Fig. 6. Cross-section of blood flow before and after static suppression. (a) No static voxel suppression; (b) after static voxel suppression;(c) OCT original structural image

下载图片查看所有图片

图7分别展示了5组数据采用不同方法得到的成像结果，从最后一组数据的放大图中可以看出，普通的Pearson相关法效果最差，而采用了分频谱思想的SSADA算法对图像质量有一定的提升，但仍然逊色于SVOCT算法。同样在分频谱的情况下，本文算法对图像质量提升明显，大部分静态体素被有效地抑制，黄斑区的毛细血管更加清晰，图像可视性和对比度显著提高。

图 7. 本文提出算法与其他强度OCTA的结果对比。（a）提出算法的结果；（b）SSADA算法的结果；（c）SVOCT算法的结果；（d）Pearson相关法的结果。实线框为感兴趣区域

Fig. 7. Result comparison of proposed algorithm and other strength OCTA. (a) Results of proposed algorithm; (b) results of SSADA algorithm; (c) results of SVOCT; (d) results of Pearson correlation method. Solid boxes are regions of interest

下载图片查看所有图片

由于眼底毛细血管较为密集，难以划定静态组织区域，故采用式（7）计算四种算法的原始图像对比度来衡量算法的效果：

\begin{matrix} C = \sum δ {(i, j)}^{2} P_{δ} (i, j), \end{matrix}

（9）

式中： $δ (i, j) = |i - j|$ 表示相邻像素灰度差； $P_{δ} (i, j)$ 表示灰度差为 $δ$ 的像素分布概率。在采用同样的预处理、目标提取、图像配准、投影等操作时，本文提出算法与SSADA、SVOCT、Pearson相关映射法相比具有最高的对比度（如表1所示），能提供较好的视网膜微血管图像。

表 1. 四种算法的对比度

Table 1. Contrast ratios of four algorithms

Algorithm	Contrast ratio
Proposed	0.0201
SSADA	0.0018
SVOCT	0.0021
Pearson corellation	0.0003

查看所有表

无论是高信噪比的皮肤组织，还是低信噪比的眼底视网膜，本文算法都能取得较好的效果，尤其是对于低信噪比样品，本文算法能提供较大的改善，显著提高了血管可视性，为医疗领域提供了更可靠的诊断依据。

4　结论

提出了一种基于强度信息的分频谱并结合局部信噪比去除静态体素的OCTA算法，以分频谱的思想削弱轴向运动伪影的同时增强了血流信息，拓展了一组振幅信号的去相关与逆信噪比计算公式并模拟了静态和噪声体素去相关值和局部信噪比之间的分布规律，为不同局部信噪比区域分配不同的动静态阈值提供了依据。由于包含了局部窗口内的相关运算，本文提出算法也存在细节模糊的问题，模糊程度与参与运算的时空核大小有关。在重复采样次数为2的情况下，综合考虑模拟实验分布结果和眼底红细胞尺寸等因素，选取了5×5的窗口大小，但当重复采集次数增加的情况下，空间核大小可以随之减小，在血管造影图像质量不变的情况下减轻模糊现象。空间核的大小实际上也需要结合样品种类和扫描空间间隔进行考量，当慢轴方向扫描的实际距离小、图像更细腻时，空间核也可以随之增加，只要保证空间核范围内的体素大致服从同种分布即可。

此外，背景噪声方差参与了每一个局部逆信噪比的计算，具有实际意义，故其值必须准确，因此如何自适应地提取一个纯背景区域是一个难点。在提取准确的情况下，相比于其他强度OCTA方法，本文提出的方法在高、低信噪比数据中均能获得更好的微血管图像质量，能够为疾病诊断提供更可靠的依据。

参考文献

[1] 屈慧, 汪毅, 娄世良, 等. 纯随机相位板散斑去相关光学相干层析成像[J]. 光学学报, 2023, 43(1): 0111002.

Qu H, Wang Y, Lou S L, et al. Speckle decorrelation optical coherence tomography using pure random phase plate[J]. Acta Optica Sinica, 2023, 43(1): 0111002.

[2] Swanson E A, Izatt J A, Hee M R, et al. In vivo retinal imaging by optical coherencetomo graphy[J]. Optics Letters, 1993, 18(21): 1864-1866.

[3] Huang D, Swanson E A, Lin C P, et al. Optical coherence tomography[J]. Science, 1991, 254(5035): 1178-1181.

[4] 时群, 冯锦平, 郑烨, 等. 基于微激励光学相干弹性成像的角膜固频在体测量[J]. 光学学报, 2022, 42(10): 1012005.

Shi Q, Feng J P, Zheng Y, et al. Micro-Force Optical Coherence Elastography for in vivo Corneal Natural Frequency Measurement[J]. Acta Optica Sinica, 2022, 42(10): 1012005.

[5] Kim H, Eom T J, Kim J G. Vascular morphometric changes during tumor growth and chemotherapy in a murine mammary tumor model using OCT angiography: a preliminary study[J]. Current Optics and Photonics, 2019, 3(1): 54-65.

[6] Tan B Y, Chua J, Barathi V A, et al. Quantitative analysis of choriocapillaris in non-human primates using swept-source optical coherence tomography angiography (SS-OCTA)[J]. Biomedical Optics Express, 2018, 10(1): 356-371.

[7] Zhang J, Nguyen T, Potsaid B, et al. Multi-MHz MEMS-VCSEL swept-source optical coherence tomography for endoscopic structural and angiographic imaging with miniaturized brushless motor probes[J]. Biomedical Optics Express, 2021, 12(4): 2384-2403.

[8] 马艳红, 张鹏飞. 基于OCT高精度定量的视网膜光生理功能成像技术[J]. 中国激光, 2023, 50(9): 0907105.

Ma Y H, Zhang P F. High-accuracy quantitative functional imaging technique for retinal opto-physiology with optical coherence tomography[J]. Chinese Journal of Lasers, 2023, 50(9): 0907105.

[9] 苗清然, 王海霞, 于洋, 等. 光学相干层析成像在指尖生物特征识别中的应用[J]. 激光与光电子学进展, 2023, 60(8): 0811012.

Miao Q R, Wang H X, Yu Y, et al. Application of optical coherence tomography in fingertip biometrics[J]. Laser & Optoelectronics Progress, 2023, 60(8): 0811012.

[10] 张璇晔, 朱疆. 内窥光学相干层析成像的研究进展与应用[J]. 中国激光, 2023, 50(21): 2107103.

Zhang X Y, Zhu J. Research progress and applications of endoscopic optical coherence tomography[J]. Chinese Journal of Lasers, 2023, 50(21): 2107103.

[11] Chen Z P, Zhao Y H, Srinivas S M, et al. Optical Doppler tomography[J]. IEEE Journal of Selected Topics in Quantum Electronics, 1999, 5(4): 1134-1142.

[12] Chen Z P, Zhao Y H, Nelson M D J S. Phase-resolved optical coherence tomography and optical Doppler tomography: technology and applications[J]. Proceedings of SPIE, 2001, 4241: 182-185.

[13] An L, Qin J, Wang R K. Ultrahigh sensitive optical microangiography for in vivo imaging of microcirculations within human skin tissue beds[J]. Optics Express, 2010, 18(8): 8220-8228.

[14] Wang R K, An L. Doppler optical micro-angiography for volumetric imaging of vascular perfusion in vivo[J]. Optics Express, 2009, 17(11): 8926-8940.

[15] Wei X, Hormel T T, Jia Y L. Phase-stabilized complex-decorrelation angiography[J]. Biomedical Optics Express, 2021, 12(4): 2419-2431.

[16] Miao Y S, Siadati M, Song J, et al. Phase-corrected buffer averaging for enhanced OCT angiography using FDML laser[J]. Optics Letters, 2021, 46(16): 3833-3836.

[17] Blatter C, Weingast J, Alex A, et al. In situ structural and microangiographic assessment of human skin lesions with high-speed OCT[J]. Biomedical Optics Express, 2012, 3(10): 2636-2646.

[18] Huang Y P, Zhang Q Q, Thorell M R, et al. Swept-source OCT angiography of the retinal vasculature using intensity differentiation-based optical microangiography algorithms[J]. Ophthalmic Surgery, Lasers & Imaging Retina, 2014, 45(5): 382-389.

[19] Mariampillai A, Standish B A, Moriyama E H, et al. Speckle variance detection of microvasculature using swept-source optical coherence tomography[J]. Optics Letters, 2008, 33(13): 1530-1532.

[20] Jonathan E, Enfield J, Leahy M J. Correlation mapping method for generating microcirculation morphology from optical coherence tomography (OCT) intensity images[J]. Journal of Biophotonics, 2011, 4(9): 583-587.

[21] Cha M, Phillips R D, Wolfe P J, et al. Two-stage change detection for synthetic aperture radar[J]. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53(12): 6547-6560.

[22] Jia Y L, Tan O, Tokayer J, et al. Split-spectrum amplitude-decorrelation angiography with optical coherence tomography[J]. Optics Express, 2012, 20(4): 4710-4725.

[23] Enfield J, Jonathan E, Leahy M. In vivo imaging of the microcirculation of the volar forearm using correlation mapping optical coherence tomography (cmOCT)[J]. Biomedical Optics Express, 2011, 2(5): 1184-1193.

[24] Husvogt L, Ploner S B, Chen S Y, et al. Maximum a posteriori signal recovery for optical coherence tomography angiography image generation and denoising[J]. Biomedical Optics Express, 2020, 12(1): 55-68.

[25] 周晖. 扫频OCT系统中相对强度噪声来源与抑制技术研究[D]. 南京: 东南大学, 2015.

ZhouH. Source and suppression of relative intensity noise in swept-source OCT system[D]. Nanjing: Southeast University, 2015.

[26] 樊金宇. 高灵敏扫频光学相干层析在体实时血管成像技术的研究[D]. 合肥: 中国科学技术大学, 2021.

FanJ Y. Research on high-sensitivity swept-source optical coherence tomography in vivo real-time vasulature imaging[D]. Hefei: University of Science and Technology of China, 2021.

[27] 张舞杰, 李迪, 叶峰. 基于Sigmoid函数拟合的亚像素边缘检测方法[J]. 华南理工大学学报(自然科学版), 2009, 37(10): 39-43.

Zhang W J, Li D, Ye F. Sub-pixel edge detection method based on sigmoid function fitting[J]. Journal of South China University of Technology (Natural Science Edition), 2009, 37(10): 39-43.

2.4　通过Sigmoid函数映射过渡处理阈值附近的体素

王露桐, 汪毅, 徐玉帅, 娄世良, 蔡怀宇, 陈晓冬. 基于局部信噪比阈值调节的分光谱去相关光学相干断层扫描血管造影方法[J]. 光学学报, 2024, 44(5): 0517001. Lutong Wang, Yi Wang, Yushuai Xu, Shiliang Lou, Huaiyu Cai, Xiaodong Chen. Split-Spectrum Threshold Decorrelation Optical Coherence Tomography Angiography Method Based on Local Signal-to-Noise Ratio[J]. Acta Optica Sinica, 2024, 44(5): 0517001.