基于多波长数字锁相检测的光学拓扑成像方法

高 峰 ，朱苹苹，易 茜，张 伟，陈 琛，张丽敏

(1.天津大学精密仪器与光电子工程学院，天津300072；2. 天津市生物医学检测技术与仪器重点实验室，天津300072)

光学拓扑成像是指在组织体表面利用反射测量给出组织体表面下浅层的光学参数二维变化．光学拓扑成像主要应用于脑功能成像研究，即监测生理或心理活动引起的大脑皮层内各调控区的应激反应过程[1]，其作为一种近红外光谱测量方法，具有完全无创性，高时间分辨率，以及可在自然条件下长时间测量等诸多优点．目前它已发展成为临床诊断和监护中的有效手段，成为现有脑功能成像模态(功能磁共振成像、脑电图等)的一种极为有益的补充．

利用反射光强信息，光学拓扑成像可以重建出组织体内部重要生理参数的变化分布图[2]，例如氧合血红蛋白和脱氧血红蛋白的浓度变化分布．虽然此成像方式不能测量确切的血容量浓度，但是可以测量分子物质的浓度变化[3]，可实时连续无创检测血氧的相对浓度变化，提供脑组织中血氧的动态变化情况，已成为临床诊断和监护的重要手段之一．

考虑到光学拓扑成像方法的优点以及应用前景，多家公司已经研制出光学拓扑成像系统，但其对探测器要求很高，成本昂贵难以普及．笔者提出的基于数字锁相检测技术的光学拓扑成像系统，采用光开关简化了系统对光源的编码过程，并且利用数字锁相检测技术提高了系统信噪比，实现了多波长同步动态测量，降低了系统对探测器的要求．

光学拓扑成像可以选择时间分辨测量、频域测量和连续波测量3 种模式．经比较，连续波测量有更实用的优点：高信噪比，高时间分辨率，动态范围大，简单性，成本低，临床方便等[2，4]．故本文选择连续波测量模式实现光学拓扑成像．假定每个测量区域的散射系数是均匀的，且在相邻的2 次生理活动前后不发生变化[5]，在此基础上以血红蛋白的吸收来测量组织中的氧合情况．

光入射到组织体后，其强度会以指数形式衰减，而数字锁相检测技术就是一种可以实现微弱光信号动态检测的有效技术．相对于其他数字检测技术，数字锁相检测技术可以实现多路调制信号的并行检测，在相同的采样时间内有更好的抗噪性能[6-7]．增大源-探测光纤间距，使得多波长耦合信号穿透表面组织到达目标位置时，依然可以利用此检测技术于噪声中提取到有用信息．

本文详细介绍了每个组成部分，设计出符合要求的源-探测光纤排布以及测量支架，并且模拟组织体光学参数进行了液态仿体的实验研究．光学拓扑成像实验过程中，系统将两路不同反射波长的激光信号耦合后作为入射光源，并以数字锁相技术同步检测反射光强．以修正的Beer-Lamber 定律为模型计算出每对邻近的源-探测光纤处生理或心理活动前后血红蛋白的浓度变化量，利用此信息大致反映出大脑皮层内部的氧合状况．

1 光学拓扑成像系统

如图1 所示，基于数字锁相检测技术的光学拓扑成像系统主要由光源系统、组织体与测量支架、检测处理系统组成．此系统可实现多路调制信号入射组织体，同步信号监测、数据采集以及并行数字锁相处理等功能．光源系统由单片机控制的两路调制信号源、激光器控制模块、激光器输出尾纤、波分复用器以及光开关组成，可以提供多路入射光源．检测处理系统由光电探测器、数据采集卡以及PC 机组成，可以并行检测、采集组织体表面四路反射光强，并将采集的信息输入到PC 机中进行锁相处理、数值计算．根据修正后Beer-Lamber 定律，可以由处理后获得的幅值信息计算出采样点处的 aμΔ ，进行线性插值后即可得到重建的图像．

图1 系统组成Fig.1 System configuration

为更好地实现本文设计的光学拓扑成像系统，选择了符合设计需求的元器件和设备，得到该系统的实物连接图，见图1(b)．下面详细介绍具体的选择情况.

1.1 光源系统

为了更好地区分氧合血红蛋白和脱氧血红蛋白[8]，本文采用波长分别为660,nm 和830,nm 的温控半导体激光器系统作为辐射光源，其由温度-电流控制模块(LTC100-B，美国Thorlabs Inc)和激光尾纤(分别为LPS-PM660-FC 以及LPS-PM830-FC，美国Thorlabs Inc)两部分组成．此外为了实验过程的简易可行，采用波分复用器(WDMs，加拿大OZ Optics Ltd)耦合两路光信号，同时采用光开关实现辐射光源入射位置的切换．

为提取反射出的微弱有用信号，需要以一定频率调制激光器．如此数字锁相检测技术才能更好地去除背景噪声，分辨出多波长耦合信号中的调制信号[7]．本文采用单片机接收LABVIEW 发送的信息，控制程控放大器(LT6912，美国Linear Technologies Inc)产生正弦调制信号的特定幅度，同时控制直接数字频率合成器(direct digital synthesizer，DDS，AD9831，美国Analog Device Inc)产生正弦调制信号的特定频率．以此为基础，经过设计、制板与调试，得到多路可调频可调幅的调制信号源[9]．

1.2 组织体与测量支架

采用固态仿体和液态仿体两类进行实验．系统可行性验证阶段，选用的是固态标准仿体(吸收系数μa=0.01mm-1，约化散射系数μs′=1mm-1).在光学拓扑成像验证阶段，暂不能进行人体实验，故选用的是单独配置的液态仿体．

仿体实验时选用不同的测量支架搭配进行．以固态仿体进行实验时，采用的是平移支架，具体实现过程在后续实验中介绍．液态仿体是目前验证阶段的主体部分，在此详细介绍液态仿体的配置以及实验中用到的测量支架．

根据朗伯定律

式中：I0为初始光强；I 为测量得到的光强；μa为吸收系数；e 为常数2.718；ε(λ)为摩尔吸光系数，不随浓度C 和光程的改变而改变．由式(1)和式(2)可知，直接计算出 Δμa，即可反映出血红蛋白浓度的相对浓度变化信息．

为了相对精确地测量 aμΔ ，实验过程中使用液态仿体实现吸收系数的变化．将Intralipid 溶液与印度墨水以相应比例混合后即可配置出所需光学参数的液态仿体[10]，其中混合比例需根据所需光学参数来改变．同时，为了改善拓扑成像效果、提高空间分辨率，还要考虑源-探测光纤的分布情况以及异质体的放置位置[11]．图 2(a)展示了设定的光纤排布情况．它满足以下规则：①每对源入射光纤和探测光纤之间的距离为20,mm，此距离直接影响光程；②每对源入射光纤和探测光纤的中间点定义为采样点，且探测光纤测得的信息反映的即是中间点处的 aμΔ ；③每个采样点和相邻位置处的采样点等距．

在40,mm×40,mm 的源-探测光纤分布区域，共有 12 个采样点，相邻位置的采样点间隔为14.14,mm．为获得相对高的空间分辨率，将异质体放置在采样点位置处．图2(b)为笔者设计的测量支架，能够放置源光纤、探测光纤和塑料试管，其中塑料试管直径10,mm，用来放入异质体．小孔和大孔分别为光纤和异质体放置位置．

图2 光纤排布Fig.2 Optode array in the phantom experiment

1.3 检测处理系统

信号检测部分采用PIN Si光电二极管(PDA36A，美国Thorlabs Inc)实现微弱光信号探测[12]，数据采集卡(PCI-1712L，台湾研华Advantech Ltd)实现多通道同步数据采集．考虑到设定的源-探测光纤间距以及光电探测器自身灵敏度的限制，经探测器放大之后的信号依然比较微弱时，可以对此时的信号进行简单的滤波放大处理，以期提高数据采集的有效性．

采集到PC 机内的离散信号为

式中：A0为直流分量；Ak为2 个频率 fk(k = 1,2)的谐波幅度；fs为采样频率；fk为调制频率．它们之间满足

式中 Ns为数据采集卡设定的采样点数．本文采用的调制信号为10,kHz、15,kHz，采样频率为400,kHz，采样点 Ns为2 000．

本文需要对采集到的数据进行多波长数字锁相检测处理，此过程中需要两路离散的参考信号：同相参考信号和正交参考信号，此两路信号要有90°的相位差，其公式为

将式(3)中的 ()M n 与式(5)和式(6)中的两路参考信号 ()C n 、 ()S n 分别相乘，并通过低通滤波器滤除交流分量，此过程可以通过离散傅里叶变换的(discrete Fourier transform，DFT)矩阵来实现，即

于是，可得被测信号的幅值信号为

根据修正后的Beer-Lamber 定律，被测信号的幅值信息与血红蛋白浓度的关系为式中：Coxy和 Cdeoxy分别为含氧血红蛋白与脱氧血红蛋白的浓度；Aio(λi)为波长λi的入射光强；Ai(λi)为静止状态时波长λi的输出光强；L 为设定的源-探测器光纤间距；B 为光程因子，本文引用其经验值；ε为血红蛋白的摩尔吸光系数，在特定波长处为常数[13]．

脑活动状态时，被测信号的幅值信息与血红蛋白浓度的关系为

式中 Ai′ (λi)为血红蛋白浓度变化后(活动状态)波长λi的输出光强．因此，血红蛋白的相对浓度变化量如式(11)所示：

由式(11)即可计算出血红蛋白浓度的相对变化量．目前系统暂时处于模拟验证阶段，本文以吸收系数的变化来反映血红蛋白浓度的变化情况，故后续实验过程中，探测的是吸收系数的相对变化量．

2 实 验

2.1 光学拓扑成像方法的可行性验证

实验过程中简化系统，以验证光学拓扑成像方法的可行性，即光源系统不采用光开关，直接将波分复用器耦合后的光入射仿体，同时检测处理系统只采用一个光电探测器以及一个采样通道．本文采用固态标准仿体模拟组织体，其光学参数：μa=0.01,mm-1，μs′=1,mm-1(800,nm 波长下测得)，并采用2 个固定支架固定源光纤和探测光纤．

实验分为以下两类．

(1) 以10,kHz 正弦信号调制发射波长为830,nm的激光器，15,kHz 正弦信号调制发射波长为660,nm的激光器．固定源光纤的支架，以0.5,mm 为步长移动探测光纤支架，实现不同的源-探测光纤间距下幅值信息的检测．将经过数字锁相处理并归一化后的信号，与外推条件下距离源光纤一定水平距离处反射光强归一化后的解析解相比较，以此来判断系统所测信号幅值的有效性．

图3中测量值与解析解并不完全一致．这是因为不同波长下，光学参数不同，基于半无限平面模型所计算出的反射光强不同，因此正常情况下测量的幅值与解析解会有一定的差异．由以上结果可见，入射光经过仿体后，系统测得的幅值随着源-探测光纤间距的改变而改变，其趋势与理论上解析解的变化趋势大体一致．

(2) 理论上，输入的调制信号幅值与多波长数字锁相检测处理后得到的信号幅值成正比关系．为相对准确地判断系统的测量结果，对于两路调制信号，固定其中一路调制信号的幅值，同时有规律地改变另一路调制信号的幅值，以此实现较准确的验证目的．

图3 实验测量值与解析解的对比Fig.3 Comparison between the phantom experimental result and the analytical solution

图4所示分别为维持15,kHz 正弦调制信号幅值不变，将10,kHz 正弦调制信号幅值从600,mV 以50,mV 为步长依次降低至100,mV 时的测量幅值曲线，以及维持10,kHz 正弦调制信号幅值不变，将15,kHz 正弦调制信号幅值从600,mV 以50,mV 为步长降低至100,mV 时的测量幅值曲线．

图4 激发光波长和调制频率分别为830,nm、10,kHz 和660,nm、15,kHz 时测量幅值与调制信号峰峰值的关系曲线 Fig.4 Normalized amplitudes measured with the input single changing using the excitation light of different wavelengths and modulation frequencies：830,nm，10,kHz；660,nm，15,kHz

由以上结果可知，调制信号的输入幅值与数字锁相检测处理后得到的幅值成正比关系．测量中可忽略手动调节的实验误差和相同调制信号峰峰值时不同波长带来的吸收系数差异．综合以上两类实验结果可知，光学拓扑成像方法是可行的，系统测量的幅值信息是有效的．

2.2 光学拓扑成像方法的实验验证

2.2.1 单异质体

背景液态仿体的光学参数设定为：aμ=0.005 mm-1，λ2＝660,nm，f2＝15,kHz，s1μ′ = mm-1，异质体的参数为：μa=0.03mm-1，μs′ = 1mm-1(以上光学参数以660,nm 波长处测得的吸光度来配置)，重建图像见图5．图5(a)、(b)分别为830,nm、660,nm 波长处Δμa的重建图像，异质体位置处求得的 Δμa分别为：Δμa(830) =0.003 9,mm-1和 Δμa(660)=0.005,9,mm-1．

在660,nm 波长处，理论吸收系数的相对变化量为： Δμ（a6 60） = 0.025mm-1，测 量 结 果 的 误 差 为76.4%．造成该误差的原因分析如下：首先，在实验过程中不可避免地会受到外界杂散光的影响，且系统检测到的反射光强度十分微弱，造成了较低的信噪比；其次，配置的液态仿体的吸收系数和散射系数分别由光谱仪测量和经验值计算所得，且液态仿体因重力影响导致纵向密度分布不均匀；此外，实验过程中还存在人为因素造成的误差．综上所述，定性的结果误差很大，但是由成像图可知，系统可以有效地分辨异质体且定位较准确，异质体对背景其他位置的影响不大，可用来定性分析、分辨吸收系数的变化．

图5 激发光波长分别为830 和660,nm 时单异质体的Δμa重建图像 Fig.5 Reconstructed-images of the single inhomogeneous target using the excitation light of different wavelengths(830,nm，660,nm)

2.2.2 不同吸收系数的双异质体

背景液态仿体的光学参数设定为：μa=0 .02mm-1，约化散射系数μs′ = 1mm-1，异质体光学参数为：μa= 0.04mm-1和μa= 0.06mm-1，μs′ = 1mm-1(以上光学参数是以波长660,nm 处的吸光度来配置)，重建图像如图6 所示．图6(a)为830,nm 波长处 Δμa的重建图像，上、下侧异质体位置处测得的 Δμa分别为：Δμa′(8 30) =0.006,1,mm-1和 Δμa(8 30) =0.003,1,mm-1． 图6(b)为660,nm 波长处 Δμa的重建图像，上、下侧异质体位置处测得的 Δμa分别为：Δμa′(6 60)= 0.007,2,mm-1和 Δμa(6 60)= 0.003,8,mm-1．

图6 激发光波长分别为830 nm 和660 nm 时不同吸收系数的双异质体的Δμa 重建图像 Fig.6 Reconstructed Δμa images of two inhomogeneous tar- gets with different concentrations using the excitation light of different wavelengths(830,nm，660,nm)

在660,nm 波长处，理论上吸收系数的相对变化量分别为：Δμa(6 60)= 0.02,mm-1、Δμa′ (660) = 0.04mm-1，故误差分别为81%、82%，双异质体的吸收系数比例的测量误差为5.25%．测量误差比较大，其原因与单异质体的情况相同，可以通过优化液态仿体光学参数的配置以及改善遮光条件等方法降低测量误差．另外由吸收系数变化的成像图可知，系统可以有效地分辨异质体且定位较准确，双异质体的吸收系数比例与理论值基本一致，即系统可以对吸收系数的变化趋势正确成像．

［5］Gao F，Zhao H J，Tanikawa Y，et al. Optical tomographic mapping of cerebral haemodynamics by means of time-domain detection：Methodology and phantom validation[J]. Physics in Medicine and Biology，2004，49：1055-1078.

［6］Lasker J M，Masciotti J M，Li Y，et al. Dynamic optical tomographic imager with optimized digital lock-in filtering[J]. SPIE-OSA Biomedical Optics ， 2007 ，6629(3)：1-11.

［7］Masciotti J M，Lasker J M，Hielscher A H. Digital lock-in detection for discriminating multiple modulation frequencies with high accuracy and computational efficiency[J]. IEEE Transactions on Instrumentation and Measurement，2008，57(1)：182-189.

［8］Everdell N L，Gibson A P，Tullis I D C. A frequency multiplexed near-infrared topography system of imaging functional activation in the brain[J]. Review of Scientific Instruments，2005，76(9)：093705-1-093705-5.

［9］焦玉婷. 动态扩散光学层析成像中多波长数字锁相光电检测技术的研究[D]. 天津：天津大学精密仪器与光电子工程学院，2010.Jiao Yuting. Research on Multi-Wavelength Digital Lock-In Photoelectric Detection in Dynamic Diffuse Optical Tomography[D]. Tianjin：School of Precision Instrument and Opto-Electronics Engineering，Tianjin University，2010(in Chinese).

［10］van Staveren H J，Moes C J M，van Marie，et al. Light scattering in Intralipid-10% in the wavelength range of 400-1100,nm[J]. Applied Optics ， 1991 ， 30(31)：4507-4514.

［11］Yamamoto T，Maki A，Kadoya T，et al. Arranging optical fibres for the spatial resolution improvement of topographical images[J]. Physics in Medicine and Biology，2002，47：3429-3440.

［12］安毓英，曾晓东. 光电探测原理[M]. 西安：西安电子科技大学出版社，2009.An Yuying，Zeng Xiaodong. Photoelectrical Detection Theory[M]. Xi’an：Xi’an University of Electronic Science and Technology Press，2009(in Chinese).

［13］Matcher S J，Cooper C E，Cope M，et al. Performance comparison of several published tissue near-infrared spectroscopy algorithms[J]. Analytical Biochemistry，1995，227：54-68.