基于相位分辨算法的OCT弹性测量方法的研究

Study of OCT Elasticity Measurement Method based on Phase-Resolved Algorithm

摘要: 生物组织力学特性的变化对于医学研究具有重要的研究意义。为了实现对弹性信息的提取和量化，本文基于剪切波量化杨氏模量的相关理论和扫频光学相干层析成像技术提出了一种弹性测量方法，在此基础上设计并搭建了一套基于扫频OCT的光学相干弹性成像系统。系统使用压电堆栈输出双脉冲激励产生剪切波，通过相位分辨算法实现了对样品内部速度信息的提取和测量，并在此基础上实现对样品剪切波的测量和弹性模量的量化。通过对琼脂样品施加双脉冲激励进行剪切波量化杨氏模量的实验测试，测得的琼脂样品平均杨氏模量为23.29 kPa，与理论结果基本一致，从而验证了基于本系统的弹性测量方法的可行性。

Abstract: The investigation of changes in the mechanical properties of biological tissues holds paramount significance in medical research. This paper presents an elasticity measurement method aimed at extracting and quantifying elasticity information. The proposed method integrates theories related to the quantification of Young’s modulus using shear waves and utilizes swept optical coherence tomography (OCT) techniques. A novel optical coherence elasticity imaging system based on swept OCT is designed and constructed as a platform for this method. The system employs piezoelectric stack-generated double-pulse excitation to induce shear waves and utilizes a phase-resolu- tion algorithm for the extraction and measurement of the internal velocity of the sample. Furthermore, the system facilitates the quantification of sample elasticity through shear wave measurement. Experimental validation of the proposed method is conducted on agar samples utilizing double-pulse excitation to quantify the Young’s modulus of shear waves. The measured average Young’s modulus of the agar samples is found to be 23.29 kPa, demonstrating close alignment with theoretical expectations and thereby affirming the feasibility of the proposed elasticity measurement method facilitated by this system.

1. 引言

在医学研究中，生物组织力学特性的变化往往能够反映出某些病理学上的改变，因此对于组织弹性信息的提取和量化能够有效的帮助提高临床诊断的准确性，具有重要的研究价值。临床上针对生物组织弹性信息的检测，早期多以触诊法为主 [1] [2] 。然而触诊法的结果受医生主观判断的影响较大，为了能够准确、高效的定量分析生物组织的力学特性，近年来临床上在传统医学影像诊断技术的基础上发展出了弹性成像技术。弹性成像技术的原理与触诊法类似，主要步骤包括 [3] ：施加激励使被测组织产生形变或振动；观测组织在激励作用下产生的形变或振动；据此对组织的弹性模量等生物力学特性进行评估。

光学相干弹性成像技术(Optical Coherence Elastography, OCE)是一种近年来兴起的一种新型弹性成像技术，该技术将弹性测量原理与光学相干层析成像技术(Optical Coherence Tomography, OCT)相结合，在采集被测样品内部结构信息的同时，重构其在激励作用下的位移或形变，从而实现弹性测量。与超声弹性成像 [4] 、核磁共振弹性成像 [5] 等传统弹性成像技术相比，通过OCE技术采集的弹性图像空间分辨率可达到微米级，成像深度可达到毫米级，还具备了OCT的非接触、无损伤、实时性、可多模态成像的优势 [6] ，这使得OCE技术在临床上具有极大的研究和应用前景。

OCE系统主要由激励装置和OCT系统两部分组成，其中OCT的探测方式主要有散斑追踪与相敏探测两种方式。散斑追踪法就是通过跟踪OCT图像中的散斑信息来测量组织的形变，从而推断组织的弹性特性 [7] ；相敏探测法就是据干涉信号中的相位信息推断组织的形变，进而计算组织的弹性模量 [8] 。随着散斑追踪与相敏探测两种方法的相继提出，OCE技术系统中激励加载方式和对弹性量化方式的不同可被划分为压迫型OCE、共振型OCE、剪切波OCE三类，这三种方法分别是通过计算振幅位移、共振频率、剪切波波速提取杨氏模量，实现对组织生物力学特性的评估 [9] 。

本文在扫频OCT系统的基础上，提出了一种基于相位分辨算法的弹性测量方法，该方法中设计并搭建一套基于扫频OCT的剪切波光学相干弹性成像系统，通过相位分辨算法实现了对样品内部相位和速度信息的提取和测量，并在此基础上通过施加双脉冲激励实现对样品剪切波的测量，从而实现对样品弹性模量的量化。通过对该系统的理论分析和实验测试，验证了该系统结构的有效性和稳定性。针对琼脂样品的实验测试结果表明，本文提出的方法能够良好的实现对于被测样品的弹性测量。

2. 理论分析

对于扫频OCT系统，其使用的是宽带光源，光功率谱为高斯型，记作 $S\left(k\right)$ ，原理图如图1所示。

Figure 1. Schematic diagram of SS-OCT system

图1. 扫频OCT系统原理图

在扫频OCT中，扫频光源输出的光可以表示为：

$E\left(k\right)=s\left(k\right)\exp \left[-j\left(\omega t-kz\right)\right]$ (1)

式中 $k$ 为波数， $z$ 为光程， $\omega$ 为角频率。光源输出的宽带光经光纤耦合器分成两路后分别进入参考臂和样品臂，参考臂和样品臂中返回的光可分别表示为：

$E_r=\sqrt{R_r}s\left(k\right)\exp \left[-j\left(\omega t-k{z}_r\right)\right]$ (2)

$E_s={\displaystyle \underset{i}{\sum }\sqrt{R_{si}}}s\left(k\right)\exp \left[-j\left(\omega t-k{z}_{si}\right)\right]$ (3)

式中 $R_r$ 和 $R_{zi}$ 分别表示参考臂和样品臂中深度方向第i层的反射率， $Z{}_r$ 和 $Z_{si}$ 分别表示为参考臂和样品臂中样品深度方向第i层对应的光程。

参考臂和样品臂返回光经过另一光纤耦合器后发生干涉，产生的干涉信号进入平衡探测器，在探测器处接收到的干涉信号的平均光强可以表示为：

$\begin{array}{l}{I}_D\left(k\right)=\frac{\rho }{2}\langle {\left|E_r+E_s\right|}^2\rangle =\frac{\rho }{2}\cdot \left(E_r{E}_r{}^{*}+E_s{E}_s{}^{*}+2\Re \left\{E_s{E}_r^{*}\right\}\right)\\ \begin{array}{cc}& \end{array}=\frac{\rho }{2}S\left(k\right)\{R_r+{\displaystyle \underset{i}{\sum }{R}_i}+2{\displaystyle \underset{i}{\sum }\sqrt{R_r{R}_{si}}\cos \left[2k\left(z_r-z_{si}\right)\right]}+2{\displaystyle \underset{i\ne j}{\sum }\sqrt{R_{si}{R}_{sj}}\cos \left[2k\left(z_{si}-z_{sj}\right)\right]}\}\end{array}$ (4)

公式(4)是波数 $k$ 的函数，得到的真实结果即为干涉光谱图，其中 $\rho$ 表示探测器的响应率。在公式(4)中，第一项和第二项分别表示参考臂和样品臂中后向返回光的直流项，第三项表示的是样品臂中样品深度方向所有层析面的后向散射光和参考臂返回的光的互相关项，该项为所需的样品深度图像信息，第四项表示的是样品臂中样品深度方向上不同层析面上后向散射光之间的自相关项，由于生物组织中的散射和衰减较弱，自相关项相比于互相关项往往小到可以忽略不计。公式中的第一、二、四项可通过平衡探测器的差分探测方式去除 [10] ，只保留包含样品深度信息的互相关项，因而公式(4)可以改写为：

$I_D\left(k\right)=\frac{\rho }{2}S\left(k\right){\displaystyle \underset{i}{\sum }\sqrt{R_r{R}_{si}}\cos \left[2k\left(z_r-z_{si}\right)\right]}$ (5)

实际成像中会将参考臂的等光层面调整在样品表面，即有 $z_{si}=z_r+z_i$ ，其中 $z_i$ 为样品深度方向上每层所对应的光程。公式(5)中余弦项来自于样品臂与参考臂之间的干涉，式中 $2k\left(z_r-z_{si}\right)$ 可以改写为 $2kz+\phi$ ， $\phi$ 即为该干涉光谱的相位 [11] 。通过对探测器接受到的干涉信号 $I_D\left(k\right)$ 进行逆傅里叶变换，可以得到样品的深度信息。

当样品臂中扫描的样品发生运动时，公式(5)可以改写为 [12] ：

$I^{\prime }\left(k\right)=\frac{\rho }{2}S\left(k\right){\displaystyle \underset{i}{\sum }\sqrt{R_r{R}_{si}}}\cdot {\displaystyle {\int }_0^{\tau }\cos \left[2k\left(\Delta z+vt\right)\right]\text{d}t}$ (6)

根据多普勒效应，运动样品的多普勒频移 $f_D$ 与其运动方向的速度 $v$ 的关系可表示为：

$v=\frac{{\lambda }_0{f}_D}{2n\cos \theta }$ (7)

式中 ${\lambda }_0$ 为光源中心波长， $n$ 为样品折射率， $\theta$ 为扫描光束与样品运动方向之间的夹角，即多普勒角。

样品的运动使干涉光谱产生的多普勒频移可表示为：

$f_D=\frac{\Delta \phi }{2\pi \cdot \Delta t}$ (8)

式中 $\Delta \phi$ 和 $\Delta t$ 分别表示相邻两次扫描的干涉信号的相位差和时间间隔，根据 Kasai自相关算法可将 $\Delta \phi$ 表示为 [13] ：

$\Delta \phi =\arctan \left[\frac{\mathrm{Im}\left(F_m\cdot F^{*}{}_{m+1}\right)}{\mathrm{Re}\left(F_m\cdot F^{*}{}_{m+1}\right)}\right]$ (9)

式中 $F_m$ 和 $F_{m+1}$ 分别表示相邻两次扫描的数据经FFT后的复数信号， $\mathrm{Re}\left(\right)$ 和 $\mathrm{Im}\left(\right)$ 分别表示该信号的实部与虚部， $F_m^{*}$ 是 $F_m$ 的共轭。其原理是利用公式(4)中的互相关项提取经过傅里叶变换的干涉信号的相位项，计算相邻OCT信号的相位项差值。

根据式(7)、式(8)、式(9)，被测样品的运动速度 $v$ 表示为：

$v=\frac{{\lambda }_0}{4\pi n\cos \theta }\cdot \frac{\Delta \phi }{\Delta t}$ (10)

可以看出相位差 $\Delta \phi$ 与样品的运动速度 $v$ 成正相关，依据式(10)就可以通过采集的OCT干涉信号里的相位信息来获得待测样品的速度信息和位移信息。剪切波OCE便是在此基础上获得剪切波的传播信息，通过量化剪切波波速就可以直接得到剪切模量。

剪切波波速 $V_S$ 与剪切模量G的关系可以表示为：

$G=\rho \cdot V_S^2$ (11)

式中 $\rho$ 为被测样品的密度。剪切模量G与杨氏模量E的关系可以表示为：

$E=2G\cdot \left(1+\upsilon \right)$ (12)

式中 $\upsilon$ 为被测样品的泊松比，对于生物组织泊松比一般默认取0.5，则由式(11)和式(12)可得：

$E=2\cdot \rho \cdot V_S{}^2\cdot \left(1+\upsilon \right)=3\rho \cdot V_S{}^2$ (13)

3. 实验过程与结果分析

根据上述理论分析搭建的OCE系统主要由OCT系统与激励装置两部分组成，利用OCT作为探测装置，配合激励装置实现对样品结构信息与弹性信息的检测。其中OCT系统主要由光学模块、扫描控制模块、数据采集模块等部分组成，其示意图如图2所示。光学模块包括扫频激光光源(AXP50124-8，AXSUN)、马赫曾德干涉仪、平衡探测器(PDB500C，Thorlabs)、参考臂、样品臂，其中马赫曾德干涉仪是系统的核心。扫频光源的中心波长为1310 nm，带宽为110 nm，扫频速率为100 kHz，光源输出的瞬时窄带光通过分光比为10:90的光纤耦合器分别进入参考臂和样品臂。参考臂内的光经由平面反射镜返回；样品臂内的光经由GRIN透镜和MEMS微镜组成的内窥扫描探头射向待测样品，其在不同样品深度的反射或散射光再沿原光路返回。参考臂和样品臂返回的光经过另一分光比为50:50的光纤耦合器发生干涉，产生的干涉信号进入差分平衡探测器中进行光电转换，转换成的电信号经由高速数据采集卡(ATS9350，AlazarTech)传输到计算机中，进行后续的数据处理与分析。

Figure 2. Schematic structure of SS-OCT system

图2. 扫频OCT系统结构示意图

激励装置由信号发生器(AFG3022B，Tektronix)、压电堆栈、探针构成，结构示意图如图3所示。信号发生器发出驱动信号驱动压电堆栈开始激励，再通过探针将激励传导到被测样品上诱导其产生振动，剪切波将以激励焦点为中心向四周传播。在本文的剪切波OCE系统实验中，信号发生器输出一个双脉冲信号激励压电堆栈，使被测样品内部产生剪切波，同时以OCT的扫描同步信号作为外部参考输入，从而实现激励与扫描的时序同步。

Figure 3. Schematic structure of the excitation device

图3. 激励装置结构示意图

当OCE系统中信号发生器驱动压电堆栈发生振动后，被测样品内部会对应产生剪切波，剪切波以探针与样品接触处为中心沿x轴方向向四周传播。然而，样品组织内部的剪切波传播速度可达到几米至几十米每秒，传统的OCT结构成像方式中A-scan采集速度小于剪切波的传播速度，无法在一次B-scan扫描周期内检测到剪切波的传播。为了解决这一问题，实现对组织中剪切波的量化与测量，可采用M-B扫描模式采集样品的相位信息。

M-B扫描模式的原理图如图4所示，其原理是使用OCT在样品同一空间位置上重复进行M次A-scan，定义为M-scan，通过时间上的连续扫描得到该点处深度方向上随时间变化的相位分布图，根据多普勒OCT原理该相位图能够反映该位置处样品产生振动的群速度。在完成一次M-scan后，调节位移台横向移动样品到下一个位置进行M-scan，重复该过程最终可到包含时空相位信息的M-B OCT图像。通过计算M-B OCT图像中两检测点之间的距离 $\Delta r$ 和剪切波传播 $\Delta r$ 距离所需要的时间延迟 $\Delta t$ 可将两检测点之间剪切波传播的平均速度 $V_S$ 表示为：

$V_S=\frac{\Delta r}{\Delta t}$ (14)

将式(14)计算得到的剪切波速度 $V_S$ 代入式(13)即可得到被测样品的杨氏模量。

Figure 4. Diagram of the M-B scan mode

图4. M-B扫描模式原理图

Figure 5. Diagram of the agar sample

图5. 琼脂样品实物图

为验证本文搭建的OCE系统的可行性，我们对琼脂样品进行了弹性测试实验，测量琼脂样品上两检测点之间在一定深度范围处剪切波传播的平均波速，从而测得该区域内的平均杨氏模量。如图5所示，琼脂样品主要由琼脂颗粒和0.8% v/v的脂肪乳剂混合制作而成，与脂肪乳剂混合的目的是为了增加样品的光散射系数，使OCT扫描的效果更好，经测试本实验中制备的琼脂样品密度约为1003 kg/m³。

本实验是对样品受到激励所引起的剪切波传播距离上一定距离的两个点进行检测，具体实验方法如图6所示。O点为振源点，A、B为两个检测点，这三个点在同一条直线上。将OCT扫描模式调整为M-B模式，分别对A、B两点进行扫描，可以得到这两个点在深度方向上随时间变化的相位图像。由于样品在受到激励时内部产生剪切波会导致相位发生变化，因而由这两幅时空相位图像就可以得到A、B两点间剪切波的传播情况。其中，激励装置与扫频OCT系统通过信号控制卡输出的帧同步信号进行同步控制。

Figure 6. Diagram of the elasticity test experiment

图6. 弹性测试实验原理图

实验开始时，OCT系统收到信号控制卡输出的帧同步信号，系统以M-B扫描模式在A点处进行连续的A-scan；与此同时信号发生器也收到了帧同步信号，开始向压电堆栈发出频率为1 kHz的方波双脉冲信号以驱动样品在O处沿深度方向振动。振动产生的剪切波沿着垂直于深度的方向以O点为中心向四周传播，经过一定时间到达检测点A，此时A点的相位发生变化，对应的A点处M-scan图像记录了相位变化的过程。在M-B扫描模式中，每次M-scan包含1000条A-scan，由于扫频频率为100 kHz，因此每条A-scan之间相差10 μs。对于M-scan生成的图像，图像的横坐标为时间，纵坐标为不同深度的相位差。B点处的数据采集与处理方法与A点处完全一致，通过采集两个检测点处相位差随时间的变化，就可以得到O点处激励引起的剪切波到达A、B两检测点的时间差。

对相距2 mm的两个采样点A、B在相同激励下按照上述方法进行扫描，实验结果如图7所示。其中，图7(a)和图7(c)经过M-B扫描得到的OCT相位差图，图中红色部分表明样品此时因受到脉冲激励作用相位差发生变化，颜色越深表示相位差变化的幅度越大。对比两幅图可以发现距离振源更远的B点剪切波传播到达的时间要晚于A点，且相位差变化的幅度发生了衰减。分析图中每一行的相位差分布，将相位差分布曲线上第一个波峰或波谷对应的时间视作剪切波到达该点的时间，从而就能够精确检测剪切波的传播时间。在本文的实验中，通过提取A、B两检测点深度方向上每一行的首波到达时刻，再分别取众数以作为A、B两检测点处的首波到达时刻。选取图7(a)和图7(c)中同一深度处(图中白线标记处)进行分析，对应的相位差分布曲线分别如图7(b)和图7(d)所示，可以发现A、B两点剪切波首波波峰的到达时刻分别为1.91 ms和2.63 ms，再根据式(13)和式(14)就可以计算出此处对应的平均杨氏模量为23.21 kPa。

Figure 7. Phase difference plots obtained from the measurement of the detection points: (a) and (c) are the phase difference plots of the two detection points A and B, respectively; (b) and (d) are the phase difference distribution curves at the white line markers in (a) and (c), respectively

图7. 检测点测量得到的相位差图：(a)和(c)分别为A、B两检测点的相位差图；(b)和(d)分别为(a)和(c)中白线标记处的相位差分布曲线

按照上述方法对琼脂样品上的不同位置进行多次测量实验，经计算得到样品的平均剪切波波速和平均杨氏模量的结果如表1所示。

Table 1. Experimental results of OCE measurements on agar samples

表1. 琼脂样品OCE测量实验结果

由表1可见，通过本方法测得的琼脂样品平均杨氏模量为23.29 kPa，这个结果与文献中报道的同比例配置的琼脂样品的杨氏模量相差不大 [14] 。琼脂样品的弹性成像测量实验的结果表明，本文中采用的基于扫频OCT和双脉冲激励的弹性成像方法能够通过量化样品中的剪切波波速以实现对样品杨氏模量的测量，具有较好的可行性。

4. 总结

光学相干弹性成像技术是一种近年来兴起的一种新型弹性成像技术，该技术将弹性测量原理与光学相干层析成像技术相结合，在采集被测样品内部结构信息的同时，通过重构其在激励作用下的位移或形变实现弹性测量。本文基于剪切波量化杨氏模量的相关理论和扫频光学相干层析成像技术，设计并搭建了一套基于扫频OCT的光学相干弹性成像系统，系统使用压电堆栈输出双脉冲激励产生剪切波，通过相位分辨算法实现了对样品内部速度信息的提取和测量，并在此基础上实现对样品剪切波的测量和弹性模量的量化。通过基于双脉冲激励的剪切波量化杨氏模量的实验测试，经计算得到了被测琼脂样品的杨氏模量为23.29 kPa，实验结果与目前文献中已报道的研究结果相比差距不大，从而验证了基于本系统的弹性测量方法的可行性。在后续的研究中可以考虑使用空气脉冲、声辐射力等非接触激励与本文的弹性测量方法相结合，做出进一步的改进。