上海理工大学光电信息与计算机工程学院，上海

收稿日期：2023年4月17日；录用日期：2023年7月3日；发布日期：2023年7月11日

摘要

通常在使用ODT (Optical Doppler Tomography)系统获取组织内部速度信息时必须已知振镜扫描光束与速度方向的夹角，限制了ODT系统应用场景。若振镜的扫描速度远大于流速，则可近似认为相邻的两次扫描是对同一点的流速进行测量。本文利用MEMS (Micro Electromechanical System)实现可控的快速扫描。采用虚拟双光束测量法通过将多普勒角转化为双光束之间的夹角成功解决了ODT系统采集速度信息时必须已知扫描光束与速度方向夹角的问题。采用基于最小二乘的高斯–塞德尔迭代算法提高了系统的实时性。搭建系统后对手背皮肤进行扫描，得到了皮肤浅表层的三维结构并构建了血流信号的三维模型，验证了本系统在体测量的可行性。

关键词

光学测量，ODT系统，血流信号，速度模型

Exploration and Development of Virtual Dual-Beam ODT System

Kunyue Fu, Peigang Yang, Guanglin Li

School of Optical-Electrical and Computer Engineering, University of Shanghai for Science and Technology, Shanghai

Received: Apr. 17^th, 2023; accepted: Jul. 3^rd, 2023; published: Jul. 11^th, 2023

ABSTRACT

Usually when using ODT (Optical Doppler Tomography) system to obtain internal speed information in the organization, the angle between the galvanometer scanning beam and the speed direction must be known, which limits the application scenarios of ODT system. If the scanning velocity of the galvanometer is much greater than the flow rate, it can be approximately assumed that two adjacent scans are measuring the flow velocity at the same point. This paper uses MEMS (Micro Electromechanical System) to achieve controllable fast scanning. By using the virtual dual beam measurement method, the Doppler angle is converted into the angle between the two beams, which successfully solves the problem that the angle between the scanning beam and the velocity direction must be known when the ODT system collects velocity information. The Gaussian-Seidel iterative algorithm based on least squares is used to improve the real-time performance of the system. After building the system, the skin on the back of the hand was scanned, the three-dimensional structure of the superficial layer of the skin was obtained, and a three-dimensional model of the blood flow signal was constructed, which verified the feasibility of in vivo measurement of the system.

Keywords:Optical Measurement, ODT System, Blood Signal, Velocity Model

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

由于近年来环境的恶化、生活习惯不佳等因素的影响，癌症的患病率逐年攀升 [1] 。与此同时对于癌症病灶的早期检测，目前没有较好的防治手段，主要依靠在普通内窥镜下尽可能多的采集样品组织，之后经过固定和染色等操作后在显微镜下观察组织形态。整个过程费时费力，且需要不同专业人才互相配合，采样点也无法完全覆盖整个病灶区域，非常容易发生漏检或误检的情况 [2] 。而微血管循环作为人体内最小的循环系统，对于身体状态的反应也最为敏感，非常适合作为癌症早期病灶检测的判断标准 [3] 。由于病灶通常在身体器官内部，因此一种无创的、实时的、高分辨率检测手段亟待开发 [4] 。ODT (Optical Doppler Tomography)作为一种基于低相干光干涉原理，具有无损、非接触、高灵敏度、高分辨率、实时动态成像等特点的断层成像技术，在获得高分辨率的微脉管层析结构的同时可以获得微脉管内部的血流速度信息，非常适合用于对微米级的微血管循环进行细节呈现 [5] [6] 。ODT是常规OCT (Optical Coherence Tomography)成像技术的功能性拓展，通过对相干信号的采集，解包络，去直流等一系列数据处理，可以实现高分辨率的组织结构成像和血流动力学成像。

目前对于ODT系统的功能性拓展主要集中于眼科、口腔科等可以直接进行体表观测的部位。眼科由于具有良好的光穿透性，OCT系统目前可以做到眼底视网膜网络的扫描，辅助诊疗青光眼等眼部疾病。目前在使用ODT系统对样品内部进行流速信息检测时，由于受到算法本身的限制通常需要知道扫描光束与流速方向之间的夹角才可以得到准确的流速信息。这很大程度上限制了ODT系统的进一步发展也不利于拓展其应用场景。

本文采用MEMS (Micro Electromechanical System)微镜探头搭建ODT系统的样品臂。利用MEMS微镜探头体积小，重量轻，能耗低，稳定性好的特点，通过控制微镜的翻转实现对于范围内结构信息与结构内流体速度信息的检测。利用基于MEMS的虚拟双光束测量法将多普勒角转化为双光束之间的夹角成功解决了目前主流振镜系统中必须已知光线与扫速度方向夹角的问题。采用相关性算法，利用快速傅里叶变换实现了最小二乘意义下的高斯—塞德尔迭代的快速收敛，提高系统实时性的同时消除相位噪声对于速度信息的影响，从而获得人体微血管的组织结构图、微细血管的分布信息及其中的微血流信息。

2. ODT基本原理介绍

如图1所示，扫频光源输出波长随时间变化的带宽光。带宽光经过分光器后分别进入参考臂和样品臂。参考臂上设置有固定反射镜将参考光沿原路反射回分光器；样品臂的光进入探头被探头内的反射微镜反射出探头，在样品内传播，经由样品各层面后向散射后返回分路器。参考臂的反射光与样品臂返回的光发生干涉，并被平衡探测器探测接收，完成光电转换。通过对轴向采集的点进行傅里叶变换就可以得到深度方向各层面的层析结构图，结合扫描序列通过计算机推演就可以得到样品的三维层析结构图 [7] 。

图1. ODT系统原理结构图

本文ODT系统通过测量多普勒频移实现流速场检测。对于平衡探测器检测到的干涉包络信号进行傅里叶变换，得到每一扫描点的相位信息。通过计算相邻扫描点的相位差实现流速测量。

如图2所示，当光以一定的角度θ入射到高散射介质上且介质内的物质以速度v运动时，则对于样品中的任一运动粒子A，当它以速度v在t时间内运动从位置1运动至位置2时，OCT系统的相邻两个采样点分别对粒子A采样。

图2. 多普勒测速示意图

则OCT系统接收到的位置1的相干光谱 $I_1\left(k\right)$ 为：

$I_1\left(k\right)=S\left(k\right)\left(E_R^2+2E_R{E}_s\cos \left(2knz+{\phi }_0\right)\text{d}z\right)$ (1)

其中， $S\left(k\right)$ 表示光源的光谱， $E_R$ 表示参考光的振幅， $E_s$ 表示粒子A的后向散射光的振幅，k表示波数， $nz$ 表示粒子A所在的平面对应的光程， ${\phi }_0$ 表示初始相位。

粒子A经过时间 运动至位置2，粒子A的位移 $S=tv$ ，其在纵轴的投影

$S_y=tv\cos \left(\theta \right)$ (2)

其中，θ为多普勒角度指粒子运动方向与纵轴的夹角，则OCT系统在位置2接收到的相干光谱 $I_2\left(k\right)$ 为：

$I_2\left(k\right)=S\left(k\right)\left(E_R^2+2E_R{E}_s\cos \left(2kn\left(z+tvcos\left(\theta \right)\right)+{\phi }_0\right)dz\right)$ (3)

由式1和式3可知在t时间内由粒子A的运动所产生的相位变化量 ${\Delta }_{\phi }$ 可以表示为：

$\Delta \phi =2knvt\cos \theta$ (4)

在SS-OCT系统中，干涉光谱中包含幅值和相位两部分信息。光谱做傅里叶变换后得到的复数信号可以表示为：

$\tilde{H}\left(z\right)=\left|\tilde{H}\left(z\right)\right|\exp \left(i\varphi \left(z\right)\right)$ (5)

其中 $\left|\tilde{H}\left(z\right)\right|$ 为复数信号的幅值， $\varphi \left(z\right)$ 为复数信号的相位。利用相邻两采样点的复数信号得到 ${\Delta }_{\phi }$ 代入式4可得：

$v=\frac{\Delta \varphi {\lambda }_0}{4\pi nt\cos \theta }$ (6)

其中 ${\lambda }_0$ 为光源的中心波长。ODT系统利用式6可以实现在对样品进行结构扫描的同时，检测样品内部粒子的运动信息，实现三维流速场的检测 [8] 。

值得注意的是，式6中探测光与速度的夹角θ需要预先已知，这很大程度限制了ODT广泛应用于各类场景的可能性。

本文利用MEMS快速扫描的优势，认为在完成相邻三个扫描点的信号采集时样品内部的速度信息维持不变，这一时间窗口的在本系统中为1 ms。利用相邻三个扫描点计算出对于中心扫描点而言的两个相位变化量 $\Delta {\varphi }_A$ 和 $\Delta {\varphi }_B$ 。由于MEMS的扫描角度可以由电压控制，我们认为这两个相位变化量是由两条夹角为γ光束在同一时间对同一点扫描所得结果，具体如图3所示。

图3. 双光束测量几何示意图

我们认为光线A，B在同一时间内对于同一点进行扫描，则光线B可以被视为光线 $B^{\prime }$ ， ${\theta }_A$ 和 ${\theta }_B$ 分别是光线A和光线B与流体运动方向的夹角。

由此，绝对流速V与由于流体流动而引起的相位变化量 $\Delta {\varphi }_A$ 和 $\Delta {\varphi }_B$ 之间的关系如下：

$\begin{array}{l}\Delta {\varphi }_A=V\Delta \frac{4\pi nT\cos {\theta }_A}{{\lambda }_0}\\ \Delta {\varphi }_B=V\Delta \frac{4\pi nT\cos {\theta }_B}{{\lambda }_0}\end{array}$ (7)

将式7中的两式相除可以消去绝对速度分量，由此多普勒角度 和 可以表示为：

${\theta }_A=\arctan \left(\frac{\Delta {\varphi }_A\cos \gamma -\Delta {\varphi }_B}{\Delta {\varphi }_A\sin \gamma }\right),\Delta {\varphi }_A>\Delta {\varphi }_B$ (8)

${\theta }_B=\arctan \left(\frac{\Delta {\varphi }_B\cos \gamma -\Delta {\varphi }_A}{\Delta {\varphi }_B\sin \gamma }\right),\Delta {\varphi }_B>\Delta {\varphi }_A$ (9)

其中γ为两光束之间的夹角。

将式8、9带入式6，可得流体绝对速度表达式：

$V=\{\begin{array}{c}\frac{{\lambda }_0\Delta {\varphi }_B}{4\pi nT\cos \left(\arctan \left(\frac{\Delta {\varphi }_B\cos \gamma -\Delta {\varphi }_A}{\Delta {\varphi }_B\sin \gamma }\right)\right)}\Delta {\varphi }_B>\Delta {\varphi }_A\\ \frac{{\lambda }_0\Delta {\varphi }_A}{4\pi nT\cos \left(\arctan \left(\frac{\Delta {\varphi }_A\cos \gamma -\Delta {\varphi }_B}{\Delta {\varphi }_A\sin \gamma }\right)\right)}\Delta {\varphi }_B<\Delta {\varphi }_A\end{array}$ (10)

3. 系统设计

基于扫频光源的血流OCT成像系统由光学系统模块、信号采集模块、系统控制模块和数据处理模块构成。光学系统模块中返回的样品干涉光信号，经由信号采集模块传输至计算机，经过数据处理模块实现图像的显示与储存。系统控制模块用于统一控制相干信号的扫描、传输与采集的时序同步。基于扫频光源的血流OCT成像系统的原理图如图4所示。

图4. 基于扫频光源的血流OCT成像系统原理图

图4中黑色框内为本文实验设备的光学硬件部分原理图。光源选用AXSUN公司中心波长为1310 nm，3bD带宽为80 nm，扫频速率为50 K赫兹的扫频光源。扫频光源输出的瞬时窄带宽光，通过一个分光比为98:2的光纤耦合器分别进入参考臂和样品臂。参考臂经由一个扩束准直器后被固定反射镜沿原光路反射回参考臂，形成参考光。样品臂的光经过MEMS探头的微镜反射后照射在样品上，并在样品内迅速传播。样品不同深度层面的后向散射光沿原光路返回探头形成样品光。样品光经过三端环形器后进入第二个光纤耦合器后和到达此处的参考光发生干涉形成干涉信号。干涉信号通过平衡探测器进行差分探测后转换为电信号。高速采集卡采集到干涉电信号后传输给电脑进行信号处理，得到扫描范围内组织的三维结构图像和微血管血流图像 [7] 。

图4中平衡探测器和DAQ数据采集卡组成了系统的信号采集模块。平衡探测器选用ThorlabsPDB 570C差分式平衡探测器实现干涉信号到相应电信号的转化；高速采集卡ATS9350通过与光源的kclock信号及trigger信号直连的方式与MEMS探头匹配相应的扫描速率实现对于信号的可控采集。

样品臂采用无锡微奥公司的MEMS探头，其结构如图5所示。利用GRIN透镜实现光束聚集，经由反射微镜反射出探头窗口实现侧向扫描。探头长8 mm，直径为2 mm，有效工作距离在2~3 cm [8] 。

图5. MEMS探头(a为探头结构图，b为探头实物图)

系统控制模块由计算机、DAQ高速数据采集卡、扫频光源和同步控制卡组成，互相之间配合实现系统时钟同步并控制系统扫描与采集。同步控制卡输出4路MEMS控制信号，控制MEMS微镜沿XZ平面和YZ平面在一定角度范围内旋转，同时输出一路同步信号至DAQ数据采集卡。DAQ接收光源输出的扫频信号和k-clock信号并将其作为外部时钟源，与控制采集卡的同步信号相结合，实现对干涉信号的帧同步采集。

计算机收到被测样品组织的相干干涉信号后，经过数据处理模块得到样品内部的层析结构图与三维流速分布信息，其中数据处理的过程如下：首先，对采集到的干涉信号进行去直流、高斯整流、色散补偿以提高系统的轴向分辨力；其次，对数据进行光谱整形便于后续数据处理；接下来对处理后的数据进行离散傅里叶逆变换，得到信号的幅值信息和相位信息：对于幅值信息，进行背景噪声减除、对数压缩后得到样品的层析结构图；对于相位信息，进行相关性处理、相位展开、背景去噪、峰值检测后得到相应的流速图像。

一般而言，OCT的扫描方式根据扫描范围不同分为轴向扫描方式(Aline)和横向扫描方式(Bline)两种。轴向扫描通过轴向移动，进而得到组织内部的探测深度信息。横向扫描方式指在连续轴向扫描的同时，移动样品臂使得探测光束扫过整个样品表面，从而得到三维组织信息。本文中，轴向扫描通过扫频光源改变光波波长的方式实现，横向扫描通过MEMS控制微镜探头自身转动的方式实现。

此外，为了保证采集到的信号中包含可提取可用的流速信息，信号必须满足如下两个条件 [9] ：

1) 相邻两个Aline扫描之间存在相位变化，且相位变化完全由粒子的移动所产生，因此扫描过程中不能引入额外的相位变化。

2) 为了保证数据存在足够强的相关性，在横向位置上相邻的Aline扫描光斑必须存在一定的重叠。

目前，常采用的扫描方式主要有逐点扫描和均速扫描两种。逐点扫描是指对于每一个横向位置微镜都静止一段时间，采集N条光谱，然后微镜旋转到下一角度进行扫描，将每个横向位置所采集到的N条光谱作为一组数据求平均后分别求出相邻的两条光谱之间的相位差。这种扫描方式保证了扫描过程中不会引入额外的相位但是对于流速测量而言有着一个致命的缺点：扫描速度慢。因为血管内的血流流速不固定，跟随着心脏脉搏的周期不断变化，扫描速度过慢导致同一幅扫描图像中的流速并非同一时间的速度信息，这会直接让流速测量失去意义。同时由于生物颤动等原因，会额外引入随机的相位噪声，且长时间让生物样品保持静止不动是很困难的，这同样会导致测量结果不准确。

匀速扫描方式相比于逐点扫描而言扫描速度大大提升，对于血流的测量而言可行性较高。均速扫描的问题在于扫描本身会引入额外的相位差，同时相邻横向位置扫描之间的相关性会下降。对于上述问题可以通过在后期的数据处理过程中加入消除扫描机构所引入的相位噪声的算法来解决，对于相关性的问题可以通过增加每幅图像的横向扫描点数或减少横向扫描范围的方式解决。基于上述原因，本文采用匀速栅格扫描的方式对样品进行扫描。

4. 实验验证

4.1. 扫描参数设置

如前文所述，本实验采用匀速栅格扫描方式进行。选用这一方法的核心就是在扫描速度和相邻Aline的聚焦光斑之间的覆盖率之间找到一个平衡点。当光束照射在样品上时，其单个光斑的实际扫描区域很小，扫描范围x与物镜焦距f和微镜偏转角度α成正比，具体比例关系为： $x=2f\tan \alpha$ 。物镜焦距f由探头本身决定，而微镜的偏转角度α则与实验设置的扫描电压幅值成对应关系，比例系数ρ由微镜探头本身的电压/度数比例因子决定。当扫描范围在毫米级别时，为了使得光斑覆盖率足够大，横向位置上的采样点必须足够多。若采样点密度太低，则相邻Aline之间的光斑覆盖率变低，此时将极大影响多普勒测速的准确性甚至导致测量失败。因为如果光斑覆盖率太低，样品的反射面存在一个因为扫描结构转动所产生的微小变化且变化本身所导致的相位变化远大于由血流流动所产生的相位变化。此时，稀疏采样后所得到的相位变化大部分由这一变化产生，导致流速测量准确性下降甚至失败。通常而言，我们认为当光斑覆盖率达到95%时，可以认为相邻Aline之间的相位差主要由血流变化引起，而不是扫描机构的结构所导致的。

4.2. 皮肤组织扫描实验

应用本文所设计搭建的扫描系统对手背皮肤进行检测试验，从而验证系统对于在体检测的可行性。对如图6中手背静脉血管(黑框位置)进行扫描，其中Aline设置1024个采样点，Bline设置512个采样点，扫描频率50 kHz，扫描1 mm × 1 mm的区域，微镜在Aline方向上的最大偏转角度为11.5度，Bline方向上的最大偏转角度为8.35度。

得到如图7所示的二维组织结构图，从图中可以清晰分辨皮肤的三层结构：表皮、真皮、皮下组织。成像深度最高可达1.5 mm。图中左侧箭头所指处为微血流系统，右侧箭头所指处为较大血流系统。

图8(a)所示为将所有二维扫描图像沿Bline方向叠加后的侧视图。图中可以清晰分辨出皮肤的三层结构。同时如图中方框所示，因为生物颤动的影响图像中存在纵向形变，由此可见扫描时间对于在体生物成像系统而言是一个关键影响因素。由于生物颤动属于随机误差，本文通过对同一位置的多次扫描取平均值消除此误差，得到如图8(b)所示的侧视图。从图中可以看出由生物颤动所造成的图像不连续现象得到了极大程度的优化。

图6. 手背扫描位置示意图

图7. 二维组织扫描图

图8. 三维组织侧视图

通过计算机三维重构后得到图9所示的手背皮肤三维结构图，图中可以清晰看出皮肤组织是三层结构，通过三维数据的旋转和切片可以清晰的从不同角度观察扫描位置及其深度信息。实验结果显示，本文所设计系统可以实现皮肤结构的二维扫描及三维结构重建，具有一定的临床应用价值。

图9. 手背皮肤三维结构图

4.3. 皮肤血流信息提取

应用上节中设置的扫描参数，计算可得相邻Aline之间的时间窗口为1 ms，位移为3.1 um，成像距离设置为2 mm时，光斑直径为13.16 um。相邻Aline之间的位移远小于光斑直径，我们认为相邻Aline的扫描是对于同一点的重复测量，因此可以应用前文中所述的基于MEMS的双光束测量法对样品进行速度信息的提取。

首先对原始信号做离散傅里叶变换，得到如图10所示原始相位图。

图10. 扫描位置的原始相位图

图中明暗表示了速度信号的强弱，即暗部表征速度信号较弱即实际样品的静态部分，亮部表征速度较强即实际样品的动态部分。从图中可以看出动态部分的相位信息连续性较弱，主要是由于相位信息由各数据点独立计算导致算法对于各种扰动极为敏感所致。

因此，对于离散傅里叶变化后的相位信息我们采用Yang等人提出的Kasai算法 [10] ，计算出相邻5个Aline间的平均相位变化率。得到如图11(a)所示的均值化相位变化量示意图，图中可以看到速度信息的连续性得到了大幅改善且减少了一部分背景噪声。

图11. kasai算法相位信息示意图

如图11(a)中箭头所指处所示，相位信息中存在跳变现象。这主要是由于Kasai算法中的相位变化量存在正负，其符号表示速度值的变化方向；另外由于MEMS扫描机构本身以及外部噪声信号的干扰也可能导致相位不连续 [11] 。基于上述原因，对图11(a)中各点的相位变换量取绝对值，同时对图像采用自适应滤波和背景去噪后得到如图11(b)所示的相位信息示意图。

从图11(b)中可以看出相位变化量不连续的现象得到了较大程度的缓解，但仍存在一些不足。如图11(b)中箭头所指处所示，依旧可以看到一些相位跳变的情况。这是由于Kasai算法本身的数学限制，由公式11中计算得出的相位变化量如若没有落在[−π, π]之中，则会被强制折叠进入[−π, π]范围内，其对应的速度值也会因此局限于一个较小范围内 [12] 。由此造成了速度信息的突变。

$\Delta {\varphi }_{\left(i,j\right)}=\arctan \left\{\frac{\mathrm{Im}\left[{{\displaystyle \sum }}_j^{j+N-1}{\tilde{p}}_{\left(i,j\right)\cdot {\tilde{p}}^{*}{}_{\left(i,j+1\right)}}\right]\}}{\mathrm{Re}\left[{{\displaystyle \sum }}_j^{j+N-1}{\tilde{p}}_{\left(i,j\right)\cdot {\tilde{p}}^{*}{}_{\left(i,j+1\right)}}\right]}\right\}$ (11)

因此我们需要对相位信息进行相位展开。主流的相位展开方法主要有路径跟踪算法和最小范数算法。路径跟踪算法高度依赖路径的选择，相位图中的噪点、断点、阴影等区域会造成沿不同路径所得的积分结果不同。最小范数法是一种全局相位图解包算法，其实质是一种目标函数优化算法，目的是寻找包裹相位差与真实相位的偏导数的差分值最小的情况，其中最为典型的算法是最小二乘法。基于高斯—塞德尔迭代法的最小二乘算法是离散状况下最为经典的算法，但其收敛速度较慢不适用于实时解包处理。为了提高该迭代算法的运算速度，本文利用基于快速傅里叶变换的最小二乘解包算法以提高高斯—塞德尔迭代法的收敛速度。具体算法过程如下：

在采集卡接受到的干涉信息中，包裹相位图可以表示为

${\varphi }_{i,j}={\phi }_{i,j}+2\pi k_{i,j}$ (12)

其中， ${\phi }_{i,j}\in \left[-\pi ,\pi \right]$ 为包裹相位函数， ${\varphi }_{i,j}$ 是与 ${\phi }_{i,j}$ 一一对应的真实相位值。 $i,j$ 对应M*N个数据点。在x和y方向上分别定义包裹相位差 ${\Delta }_{i,j}^x$ 和 ${\Delta }_{i,j}^y$ ，则

$\{\begin{array}{c}{\Delta }_{i,j}^x=W\left({\phi }_{i+1,j}-{\phi }_{i,j}\right)\\ {\Delta }_{i,j}^y=W\left({\phi }_{i,j+1}-{\phi }_{i,j}\right)\end{array}$ (13)

其中，W为包裹算子，其作用是保证包裹相位差 ${\Delta }_{i,j}^x$ 和 ${\Delta }_{i,j}^y$ 落在区间[−π, π]之中。

根据包裹相位差与真实相位的偏导数的差分值最小这一最终目的可得如下数学表达式

$\min \left(J\right)=\underset{i=0}{\overset{M-2}{{\displaystyle \sum }}}\underset{j=0}{\overset{N-1}{{\displaystyle \sum }}}{\left({\varphi }_{i+1,j}-{\varphi }_{i,j}-{\Delta }_{i,j}^x\right)}^p+\underset{i=0}{\overset{M-2}{{\displaystyle \sum }}}\underset{j=0}{\overset{N-1}{{\displaystyle \sum }}}{\left({\varphi }_{i,j+1}-{\varphi }_{i,j}-{\Delta }_{i,j}^y\right)}^p$ (14)

其中p为范数，本文采用最小二乘法因此p取2。

当p = 2时，将式14对 ${\phi }_{i,j}$ 求导并令其等于0，可得：

${\varphi }_{i+1,j}+{\varphi }_{i-1,j}+{\varphi }_{i,j+1}+{\varphi }_{i,j-1}-4{\varphi }_{i,j}={\rho }_{i,j}$ (15)

其中

${\rho }_{i,j}=\left({\Delta }_{i,j}^x-{\Delta }_{i-1,j}^x\right)+\left({\Delta }_{i,j}^y-{\Delta }_{i,j-1}^y\right)$ (16)

根据高斯–塞德尔迭代法

${\varphi }_{i,j}^{n+1}=\frac{\left({\varphi }_{i+1,j}^n+{\varphi }_{i-1,j}^n+{\varphi }_{i,j+1}^n+{\varphi }_{i,j-1}^n\right)-{\rho }_{i,j}}{4}$ (17)

对式17进行快速傅里叶变换后可得

${\hat{\varphi }}_{i,j}=\frac{{\hat{\rho }}_{i,j}}{2\cos \left(\pi i/M\right)+2\cos \left(\pi j/N\right)-4}$ (18)

其中 ${\hat{\rho }}_{i,j}$ 是对 ${\rho }_{i,j}$ 做快速傅里叶变换后的结果。对 ${\hat{\varphi }}_{i,j}$ 进行快速傅里叶逆变换就可以得到最小二乘意义下解包相位 ${\varphi }_{i,j}$ 。

本文的数据处理均基于仿真平台进行，图12展示了对于使用最小二乘解包算法前后程序对于单幅图像的处理速度，可以看到在使用本文提出的基于快速傅里叶变换的最小二乘解包算法后，对于图像的处理速度提高了约15%。

图12. 程序对于单幅图像的处理速度

对图11中的相位信息经过上述算法运算后得到如图13所示的相位展开图。因为最小二乘法作为一种全局解包算法，没有将相位点质量的高低纳入考虑范畴，故而存在部分误差由低质量区域传导至高质量区域的现象，由此导致图13中“拉线”情况的出现，且“拉线”主要集中于图片的左侧，由此推测MEMS本身在起始位置时的精度不高，在后续处理时可以考虑去除这一部分的图像。

图13. 相位展开后的相位信息图

我们结合前一节中得到的三维结构模型，对于相位信息进行阈值匹配计算。根据前一节得到的样品层析结构图，认为只有在层析边界内的部分才存在实际有意义的速度信息，在相位信息解算时加入了边界条件，得到如图14所示的相位信息图。

图14. 加入边界条件后的相位信息图

对处理后的相位信息，按照式10进行速度信息转换。为了提高对比度，对得到的速度信息进行峰值检测和直方图量化后得到如图15所示的速度信息图像。

系统的灵敏度取决于相邻Aline之间的相关性即当系统的总体相位方差达到临界值时系统测得的速度，也可以认为是系统可测的的最小流速。系统可测的最大流速则取决于探测器的占空比。当探测器的占空比固定时，可测的最大流速与扫描频率成正相关。本系统中测得的血流速度最大值为0.6916 mm/s，最小值为0.4120 mm/s，最小的流速变化量为0.0042 mm/s。

图15. 样品内部速度信息图像

图16. 扫描区域的三维速度模型示意图

对所有扫描数据做速度信息转换后，通过计算机三维重构建立了如图16(a)所示的速度信息三维模型。从速度信息的三维图中可以看到，皮肤内部实际存在这丰富的流体运动。从XY平面观察速度三维模型可以得到如图16(b)所示的二维平面图。从俯视图中可以清晰分辨出皮肤表层毛细血管的分布信息以及血流速度的分部信息。图中亮部展示了皮肤内部大量毛细血管中的血流信号且强度较高，相对而言分布较少的灰部可以看到都集中于亮部的周围，其表示的是强度较小的血流信号，这也符合血液沿血管流动时中心部分速度最快，靠近血管壁的部分速度较慢的特点。由此表明本文设计的基于MEMS的ODT系统可以实现在体的微血管测量及血流信号的检测。为ODT技术在临床微血管系统的检测提供了理论上和技术上的准备。

5. 结论

近年来，癌症的患病率逐年上升，对于癌症的早期病灶检查却仍然依赖于病理学的活体取样检测，费时费力的同时非常容易发生误检，给患者的身心带来巨大痛苦。医学界一直诉求一种对于早期癌症病灶的快速检测手段。本文设计的ODT系统在获得高分辨率的层析结构的同时建立速度模型，对于微血管系统的结构及速度信息都可以高效的捕捉。通过基于MEMS的双光束测量法解决了现有ODT系统中速度采集必须已知扫描夹角的问题，同时利用基于最小二乘的高斯—塞德尔迭代法提高了系统的实时性。经实验验证，本文所设计的系统对于微血管系统可以实现高精度的扫描，对于临床癌症病灶的早期检测的研究提供了一种可行的路径选择。

本文也确实存在一定的局限性，由于存在生物颤动，呼吸，心跳等因素的影响，在ODT系统的最终成像中表现为图像模糊，在组织的静态区域会产生虚假的血流信号。未来可以加入运动误差矫正算法与图像匹配算法，以减小被测组织的不自主运动对于图像质量的影响。

致谢

感谢所有在本文所设计的基于MEMS的ODT系统中提供理论指导和技术支持的各位老师同学。

文章引用

傅坤跃,杨沛刚,李广林. 虚拟双光束ODT系统探索与研制
Exploration and Development of Virtual Du-al-Beam ODT System[J]. 建模与仿真, 2023, 12(04): 3427-3440. https://doi.org/10.12677/MOS.2023.124314

参考文献