CN104146714B - 一种组织体局部血氧饱和度变化量拓扑成像装置及方法 - Google Patents

Abstract

本发明公开了一种组织体局部血氧饱和度变化量拓扑成像装置，包括光源单元、探测单元和计算机，光源单元包括激光控制器、三个不同波长的近红外激光器、波分复用器和1:9光开关；探测单元包括光纤布配器，9根源光纤、4根探测光纤、用于准直反射光的准直器、对输出光强进行调节的滤光轮、探测光子的四个光电倍增管、路由器以及多维单光子计数模块TCSPC；本发明成像方法基于最小二乘拟合算法和动态水吸收修正算法，最小二乘拟合算法具有在目标体尺寸和源探距离相近时仍能较好重建异质体形状和参数值的优势。动态水吸收修正算法能够消除水背景对光吸收的影响，进而使得氧合血红蛋白和还原红蛋白浓度的变化更加精确。

Description

一种组织体局部血氧饱和度变化量拓扑成像装置及方法

技术领域

本发明属于近红外组织光学时域拓扑成像领域，具体涉及一种基于最小二乘拟合算法和动态水吸收修正算法的三波长时域光学拓扑成像装置和方法。

背景技术

近红外光谱学(NIRS)利用组织体在“治疗窗口(600-9000nm)”的光学特性对在体组织的血氧饱和度进行探测。选择“治疗窗口(600-900nm)”范围内的波长，可以获得光在组织体中更大的穿透深度^[1]。与传统的医学诊断和成像设备(如X射线、脑电图、核磁共振)相比，近红外光成像不仅可以实现对活体组织的无损伤及实时的探测和成像，还可以直接获取组织体的生理信息^[2]。目前已有很多基于NIRS的系统被开发并应用于神经系统科学，神经心理学以及脑机接口领域。

目前脑功能成像主要采用的方法是功能磁共振成像(functional MRI,fMRI)，但该方法只能间接测量总的血红蛋白浓度(Total Haemoglobin,c_HbT)，对于更能有效反映组织体生理状态的血氧饱和度却无能为力(oxygen saturation,)。近红外光学成像方法不仅能够分别测量氧合血红蛋白浓度(Oxyhaemoglobin,c_HbO)和还原血红蛋白浓度(Deoxyhaemoglobin,c_HbR)，而且还能进一步导出血氧饱和度(SaO₂)。近年来，光学拓扑成像(Optical Topography，OT)在研究大脑皮层各调控区的应激反应过程中得到了广泛应用。光学拓扑成像即是在待测组织体表面利用反射测量(光在组织体中的漫反射路径如图1所示，其中，d表示源点和探测点之间的距离，I_O表示初始光流量，I_K表示出射光流量，scalp表示头皮组织，skull表示头骨，CSF表示脑脊液，brain表示脑髓)给出组织体表面下浅层的光学参数二维变化。例如，研究人的思维活动以及肢体运动引起的大脑皮层内血氧饱和度的变化^[2]。与其他成像方式相比，拓扑成像方法具有源探测器排布简单可变；单一源探测器对测得的数据只与其中间点处的采样点相关联，各采样点数据之间互不影响，信噪比较高；对成像目标运动的稳定性等优点^[3]。

光学拓扑成像可采用三种测量模式：连续光(Continuous Wave,CW)、频域(Frequency Domain,FD)和时域(Time Domain,TD)。CW测量方法利用光在组织中的传输的衰减来确定c_HbO和c_HbR，单一距离下无法有效区分吸收系数和散射系数的影响，测量误差较大。FD测量技术通过测量反射光相对于入射光的幅度衰减和相位延迟得到组织体吸收系数和散射系数的信息。FD测量主要是通过高频调制实现的，信噪比较低。TD测量模式用短脉冲激光作为光源，采用时间分辨测量系统对反射光进行采集，测定时间扩展函数(temporalpoint spread function，TPSF)，不但能够获得出射光强相对于入射光强的衰减信息，而且还可以获得光在组织体中传输的路径信息，从而实现吸收效应和散射效应的分离^[4]。TD模式在信息完整性、数据灵活性、系统稳定性、以及随之体现的成像质量、稳健性、多参数和多组分重建等诸多关键性能上具有其它测量模式不可比拟的综合优势，已经成为当前扩散光学成像和荧光分子成像技术的主要研究趋势^[5]。

TD测量模式主要采用同步条纹扫描相机(Streak camera)和时间相关单光子计数(Time-Correlated Single Photon Counting,TCSPC)^[6]两种技术来检测组织表面出射光的时间分布，利用光子飞行信息进行光学参数的测量和成像研究。前者时间分辨率很高，但体积大，成本高，动态范围有限；后者结合了光子计数和超快电子技术的综合优势，具有成本低、灵敏度高、动态范围宽以及时间分辨合理等一系列突出优点。本专利将TCSPC技术应用到光学拓扑成像系统中以实现具有纳秒持续时间量级的超微弱瞬态光信号的检测。

光学拓扑成像基本研究方法有三类^[7,8]：一是从基本的朗伯比尔定律(Lambert-Beer Law)出发，直接由动物实验或人体实验数据确定光衰减系数，此方法原则上仅适用于均一、无散射的介质，对于生物组织这样一个很强的光散射体并不严格适用；二是从修正的朗伯比尔定律(Modified Lambert-Beer Law,MLBL)出发，通过实验或模拟方法事先测定或计算光在组织体中传播的平均光路长来确定后续计算中的各项系数，此方法原则上只能获取光密度变化图像，从光密度的变化只能给出组织血氧饱和度相对变化趋势的曲线，定量信息受空间光子路径参数测量误差的影响较大；三是利用扩散方程理论(DiffusionEquation,DE)来描述光在组织体中的传输过程，然后通过某种运算重构出组织光学参数的变化量，求得组织血氧饱和度变化量,有望获得改进的定量精度[9]。本专利采用扩散方程理论实现对血氧饱和度的计算，不仅计算结果更精确，包含的信息量也更丰富。另外本专利采用动态水吸收修正算法扣除了水背景对光吸收的影响使得组织体血氧饱和度的变化更精确，在脑功能光学定位，肌肉运动损伤康复和新生儿脑发育实时监护等重要领域有着广阔的应用前景。

参考文献：

[1]覃东利，高峰，时间分辨组织体光学参数测量及扩散光学成像实验研究，2007。

[2]徐可欣，高峰，赵会娟，生物医学光子学(第二版)，科学出版社，2011。

[3]Teresa Correia,Sarah Lloyd-Fox,Nick Everdell.etc.Three-dimensionalOptical Topography of Brain Activity in Infants Watching Videos of HumanMovement.Physics in Medicine and Biology.2012.Vol.57:1135–1146。

[4]薛媛，高峰，基于平板检测和有限差分方法的时域乳房扩散光学层析成像算法研究，2007。

[5]王静怡，高峰，多通道时间相关单光子计数DOT/FMT系统集成和操作平台开发，2009。

[6]W.Becker，屈军乐，高级时间相关单光子计数技术，科学出版社，2009。

[7]Adam Liebert,Michal Kacprzak,Roman Maniewski.Time-resolvedReflectometry and Spectroscopy for Assessment of Brain Perfusion andOxygenation.Biocybernetics and Biomedical Engineering.2007.Vol.27:237-266。

[8]Davide Contini,Alessandro Torricelli，Antonio Pifferi,etal.Multichannel Time-Resolved Tissue Oximeter for Functional Imaging of theBrain.IEEE Transactions on Instrumentation And Measurement.2006.Vol.55:85-90。

[9]Feng Gao,Huijuan Zhao,Yukari Tanikawa,Yukio Yamada.OpticalTomographic Mapping of Cerebral Hemodynamic by Means of Time-domainDetection:Methodology and Phantom Validation.Physics in Medicine andBiology.2004.Vol.49:1055-1078。

发明内容

针对上述近红外光谱学在脑功能检测等领域应用的不足，本发明提供了一种基于最小二乘拟合算法和动态水吸收修正算法的组织体局部血氧饱和度变化量拓扑成像测量装置及方法。最小二乘拟合算法具有在目标体尺寸和源探距离相近时仍能较好重建异质体形状和参数值的优势。动态水吸收修正算法能够消除水背景对光吸收的影响，进而使得氧合血红蛋白和还原红蛋白浓度的变化更加精确。本发明集这些优势于一体，在探测脑部生理状态时具有更高的灵敏度。

本发明一种组织体局部血氧饱和度变化量拓扑成像装置，包括光源单元、探测单元和计算机，所述光源单元包括激光控制器、三个不同波长的近红外激光器、波分复用器和1:9光开关三个激光器产生的不同波长的近红外激光由多路单模光纤耦合到所述波分复用器后由一路光纤输出；所述1:9光开关用于实现在不同的光源入射位置间进行切换；激光器为近红外皮秒脉冲半导体激光器，激光器的脉冲半高宽为40～80ps；三个激光器的波长在650～950nm、且三个激光器之间波长的差值50～150nm；所述探测单元包括光纤布配器，9根源光纤、4根探测光纤、用于准直反射光的准直器、对输出光强进行调节的滤光轮、探测光子的四个光电倍增管、路由器以及多维单光子计数模块TCSPC；其中：所述源光纤用于传导照射到待测组织体表面不同源位置的入射光；所述探测光纤用于传导由待测组织体表面不同探测位置反射的漫反射光；所述光纤布配器用于布配源光纤和探测光纤在待测组织体表面的位置；所述滤光轮内装有3～6个不同衰减系数的中值密度衰减片；所述路由器用以辨别四个光电倍增管探测到的多维漫反射光信号；所述多维单光子计数模块TCSPC用以记录多维漫反射光信号的时间点扩展函数；所述激光控制器、所述波分复用器、所述滤波轮、所述路由器和所述多维单光子计数模块TCSPC均连接至所述计算机。

本发明利用上述装置的成像方法包括以下步骤：

步骤一、数据采集，包括：

利用光纤布配器将9根源光纤和4根探测光纤固定在待测组织体表面；其中，9根源光纤形成的源点按照三行三列等距布置，相邻两个源点的距离为1～3cm；4根探测光纤形成的探测点分别位于相邻4个源点的中心位置；每个源点和与其相邻的探测点之间形成有一采样点；将待测组织体处于静息状态时记作rest状态，将待测组织体处于任务状态时记作task状态；用波分复用器耦合之后的激光束依次激励所述9个源点，并用滤光轮中的衰减片进行衰减，然后再入射到4个光电倍增管PMT的阴极，将光信号转换成电信号，4个光电倍增管的电信号经过路由器辨别之后最终送达多维单光子计数模块TCSPC进行单光子计数，在与每个源点相邻的探测点同时获取采样点在rest状态和task状态下的时间点扩展函数；

将待测组织体处于rest状态时的16个采样点的时间点扩展函数表示为：

将待测组织体处于task状态时的16个采样点的时间点扩展函数表示为：

公式(1)和公式(2)中：i＝1,2,3；j＝1,2,…,16，λ_i分别表示上述三个波长，ρ_j分别表示上述16个采样点，rest表示组织体处于静息状态，task表示组织体处于任务状态，t是时间点扩展函数中的自变量；

最终得到96个时间点扩展函数，其中包括48个rest状态下的时间点扩展函数和48个task状态下的时间点扩展函数；

步骤二、数据处理，包括：

2-1、平面半无限空间外推边界条件下，组织体表面距源点ρ处的探测点反射光流量表达式如下：

公式(3)中：ρ为相邻源和探测器之间的几何距离；c为光在真空中的传播速度；κ为扩散系数；μ_a为组织体的吸收系数，z₀为实际边界距离，z_b为外推边界距离，t为时间自变量；设本专利所采用检测系统的系统响应函数为IRF(t)，实验所测得的静息状态下时间点扩展函数为IRF(t)与^restR(ρ,t)的卷积即

对公式(4)运用最小二乘拟合算法，用在步骤一中获得的rest状态下的时间点扩展函数得出每个采样点处组织体的在三个波长下的吸收系数，进一步求出组织体中氧合血红蛋白、脱氧血红蛋白和水三种物质在rest状态下各自的浓度；

再根据下述的过定方程求出消除水吸收之后修正的氧合血红蛋白和脱氧血红蛋白的浓度；

公式(5)中：和分别表示组织体处于静息状态时在三个波长下的吸收系数，和分别表示水在三个波长下的消光系数，表示组织体处于静息状态时水的浓度，和分别表示氧合血红蛋白在三个波长下的消光系数， 和分别表示脱氧血红蛋白在三个波长下的消光系数，和分别表示消除水吸收之后修正的组织体处于静息状态时氧合血红蛋白和脱氧血红蛋白的浓度。

再利用求出rest状态下组织体的血氧饱和度；

2-2、

对公式(6)运用最小二乘拟合算法，用在步骤一中获得的rest状态下的时间点扩展函数和task状态下的时间点扩展函数得出每个采样点处组织体相对于rest状态的吸收系数的变化量，进一步求出组织体中氧合血红蛋白、脱氧血红蛋白和水三种物质相对于rest状态下各自的浓度的变化量；

再根据下述的过定方程求出消除水吸收之后修正的氧合血红蛋白和脱氧血红蛋白相对于rest状态下各自的浓度的变化量；

公式(7)中：和分别表示组织体相比于静息状态时在三个波长下的吸收系数变化量，和分别表示水在三个波长下的消光系数，表示组织体相对于静息状态时水的浓度变化量，和分别表示氧合血红蛋白在三个波长下的消光系数，和分别表示脱氧血红蛋白在三个波长下的消光系数，_correctΔc_HbO和_correctΔc_HbR分别表示消除水吸收后修正的组织体氧合血红蛋白和脱氧血红蛋白的浓度相比于静息状态时的各自的变化量。

2-3、利用rest状态消除水吸收之后修正的的氧合血红蛋白和脱氧血红蛋白的浓度与task状态消除水吸收之后修正的氧合血红蛋白和脱氧血红蛋白相对于rest状态下各自的浓度的变化量求出task状态下的氧合血红蛋白和脱氧血红蛋白的浓度；再利用求出task状态下组织体的血氧饱和度；

2-4、求出每个采样点在task状态下相比于rest状态的血氧饱和度的变化量ΔSaO₂＝^taskSaO₂-^restSaO₂；

步骤三、根据步骤二获得的每个采样点在task状态下相比于rest状态的血氧饱和度的变化量，用matlab软件画出16个采样点血氧饱和度变化量的二维拓扑成像。

与现有技术相比，本发明的有益效果是：

1.当异质体直径大于源和探测器距离时，基于时域扩散方程解析解的最小二乘拟合算法可以明显改善基于修正朗伯比尔定律(MLBL)方法在量化度方面出现的过估计问题。

2.当异质体直径大于源和探测器距离时，基于时域扩散方程解析解的最小二乘拟合算法可以明显改善基于修正朗伯比尔定律(MLBL)方法出现的异质体中间值估计过高，而边缘值估计过低，导致对目标体形状的重建不准确的问题。

3.本发明直接用血氧饱和度成像并将之用作组织体生理状态的指示器。相比于氧合血红蛋白或者脱氧血红蛋白浓度的变化更能有效反映组织体的生理状态的变化。

4.本发明采用三个波长，并运用了动态水吸收修正算法，通过求解超定方程能够消除水背景对光吸收的影响，能更加精确地求解氧合血红蛋白和脱氧血红蛋白浓度及其变化量，进而可以精确求解血氧饱和度及其变化量。

5.本发明采用反射扩散光测量方式，并根据组织体的外形来合理布局测量点，有效克服了近红外光在实现组织体特别是脑组织测量时透射能力较低的缺点。

附图说明

图1是光在脑组织中的漫反射路径示意图；

图2是本发明三波长时域光学拓扑成像系统结构图；

图3是本发明中光纤帽上源‐探测光纤布配平面图；

图4是所有16个采样点血氧饱和度变化量的二维拓扑成像；

图5是三个波长的时序图以及相应的时间点扩展函数示意图；

具体实施方式

下面结合附图和具体实例对本发明技术方案作进一步详细描述。

如图2所示，本发明是一种基于最小二乘拟合的三波长时域光学拓扑成像测量装置，主要包括光源单元、探测单元和计算机。

所述光源单元包括激光控制器1、三个不同波长的近红外激光器2、波分复用器3和1:9光开关4。三个激光器2为近红外皮秒脉冲半导体激光器，激光器2的脉冲半高宽为40～80ps；三个激光器2的波长在650～950nm、且三个激光器2之间波长的差值50～150nm；选择三个激光器2的波长的优选方案是分别为660nm、780nm和830nm。激光控制器1可以调节三个激光器的强度、脉冲频率以及组合时序，用激光控制器1控制三个不同波长的激光器2发射如图5所示的时序脉冲，三个激光器2产生的不同波长的近红外激光由多路单模光纤耦合到所述波分复用器3后并由一路光纤输出；所述1:9光开关4用于实现在不同的光源入射位置间进行切换。

所述探测单元包括9根源光纤5、光纤布配器6、4根探测光纤7、用于准直反射光的准直器8、对输出光强进行调节的滤光轮9、探测光子的四个光电倍增管PMT10、路由器11以及多维单光子计数模块TCSPC12；其中：所述源光纤5用于传导照射到组织体表面不同源位置的入射光，9根源光纤芯径为50～100μm，数值孔径为0.22～0.27，一端与光开关4相连，另一端与光纤布配器6的9个不同的源点相连；所述光纤布配器6用于布配源光纤5和探测光纤7在待测组织体表面的位置；所述探测光纤7用于传导由待测组织体表面不同探测位置反射的漫反射光，4根探测光纤芯径为500～1000μm,数值孔径为0.22～0.27，四根探测光纤7的一端与光纤帽四个探测点相连，另一端与四个准直器8分别相连；所述滤光轮9内装有3～6个不同衰减系数的中值密度衰减片；所述路由器11用以区分四个光电倍增管10探测到的多维漫反射光信号；所述多维单光子计数模块TCSPC12用以记录多维漫反射光信号的时间点扩展函数。

所述激光控制器1、所述波分复用器3、所述滤波轮9、所述路由器11和所述多维单光子计数模块TCSPC12均连接至所述计算机13；所述计算机13用以对系统各个组成部分进行集成控制和成像重建。控制三个激光器2按照如图5所示所示时序发射脉冲串，并根据反射光的强弱调节滤光轮9的位置以便对反射光进行不同程度的衰减，使光子计数率低于TCSPC 12同步计数率的二十分之一，且尽量大，确保TCSPC 12有效计数，避免堆积效应引起TCSPC 12的计数损失，根据所探测到的时间点扩展函数进行拟合运算以及其后的成像重建。

采用上述测量装置获取二维拓扑成像包括以下步骤：

步骤一、数据采集，包括：

利用光纤布配器6将9根源光纤5和4根探测光纤7固定在待测组织体表面，可以使脑组织也可以使肌肉组织；其中，9根源光纤5形成的源点(图3中三角所示)按照三行三列等距布置，源点依次编号为1-9，相邻两个源点的距离为1～3cm；4根探测光纤)形成的探测点(如图3中圆圈所示)分别位于相邻4个源点的中心位置，探测点编号为1-4；每个源点和与其相邻的探测点之间形成有一采样点，共有16个源-探测对即采样点(图3中五角星所示)，如图3所示；从探测点反射的漫射光经过准直器8准直后射入滤光轮9，经滤光轮9衰减后射入光电倍增管10，光电倍增管10将接收到的微弱光信号转换成电脉冲信号经路由器送至TCSPC 12计数模块。

四根探测光纤7的一端与光纤帽四个探测点相连，另一端与四个准直器8分别相连，用于传导由组织体表面不同探测位置出射的漫反射光；源-探测器配置如图3所示，9个源点(图3中三角所示)按3行3列均匀分布，源点依次编号为1-9。四个探测点(如图3中圆圈所示)按照2行2列均匀分布，探测点编号为1-4。共有16个源-探测对即采样点(图3中五角星所示)。从探测点反射的漫射光经过准直器8准直后射入滤光轮9，经滤光轮9衰减后射入光电倍增管PMT 10，光电倍增管PMT 10将接收到的微弱光信号转换成电脉冲信号经路由器11送至TCSPC 12计数模块。

本发明方法中，将待测组织体处于静息状态时记作rest状态，将待测组织体处于任务状态时记作task状态；

用波分复用器3耦合之后的激光束依次激励所述9个源点，并用滤光轮9中的衰减片进行衰减，然后再入射到4个光电倍增管PMT10的阴极，将光信号转换成电信号，4个光电倍增管PMT10的电信号经过路由器11辨别之后最终送达多维单光子计数模块TCSPC12进行单光子计数，因为组织体不同地方血红蛋白浓度分布不均，进而对光的吸收也不同，导致不同源-探测对的反射光强也不一样。为了提高信噪比并减少堆积效应引起的计数损失，这就需要切换滤光轮9到相应的衰减片对反射光进行不同程度的衰减，控制光子计数率在同步计数率的二十分之一以内，并尽量大。根据调节之后的光强选择合适的积分之间进行积分，每切换一次光开关，就在与每个源点相邻的探测点同时获取采样点在rest状态和task状态下的时间点扩展函数；

将待测组织体处于rest状态时的16个采样点的时间点扩展函数表示为：

将待测组织体处于task状态时的16个采样点的时间点扩展函数表示为：

公式(1)和公式(2)中：i＝1,2,3；j＝1,2,…,16，λ_i分别表示上述三个波长，ρ_j分别表示上述16个采样点，rest表示组织体处于静息状态，task表示组织体处于任务状态，t是时间点扩展函数中的自变量；

最终得到96个时间点扩展函数，其中包括48个rest状态下的时间点扩展函数和48个task状态下的时间点扩展函数；

步骤二、数据处理，包括：

2-1、我们以其中一个采样点为例介绍数据处理的基本原理如下：假设氧合血红蛋白、脱氧血红蛋白和水三种物质的浓度在源点和探测点之间的局部区域是均匀的，且其他物质对光的吸收可忽略不计。

反射测量方式下，源和探测器距离相对大脑外表尺寸较小，故适用于平面半无限空间外推边界条件。此时，组织体表面距源点ρ处的探测点反射光流量表达式如下：

公式(3)中：式中，ρ为相邻源和探测器之间的几何距离；c为光在真空中的传播速度；κ为扩散系数；μ_a为组织体的吸收系数，z₀为实际边界距离，z_b为外推边界距离，t为时间自变量；设本专利所采用检测系统的系统响应函数为IRF(t)，实验所测得的静息状态下时间点扩展函数为IRF(t)与^restR(ρ,t)的卷积即

用在步骤一中获得的rest状态下的时间点扩展函数，分别进行最小二乘拟合获取组织体在三个波长下的吸收系数列如下矩阵方程

公式(5)中：和分别表示组织体在三个波长下的吸收系数； 和分别表示氧合血红蛋白在三个波长下的消光系数，和分别表示脱氧血红蛋白在三个波长下的消光系数，和分别表示水在三个波长下的消光系数，^restc_HbO，^restc_HbR和分别表示氧合血红蛋白、脱氧血红蛋白和水三种物质在rest状态下的浓度。

解式(5)可以求出^restc_HbO，^restc_HbR和

再根据下述的过定方程求出消除水吸收之后修正的氧合血红蛋白和脱氧血红蛋白的浓度；

公式(6)中：和分别表示组织体处于静息状态时在三个波长下的吸收系数，和分别表示水在三个波长下的消光系数，表示组织体处于静息状态时水的浓度，和分别表示氧合血红蛋白在三个波长下的消光系数， 和分别表示脱氧血红蛋白在三个波长下的消光系数，和分别表示消除水吸收之后组织体处于静息状态时修正的氧合血红蛋白和脱氧血红蛋白的浓度。

再利用求出rest状态下组织体的血氧饱和度；

2-2、实际临床应用中，研究者更关心光学参数的相对变化量。待测组织从rest状态变化到task状态期间，约化散射系数变化量△μ_s'远小于吸收系数改变量△μ_a，可忽略不计。结合式(4)，有如下关系式：

根据公式(8)运用最小二乘拟合算法，用在步骤一中获得的task状态下的时间点扩展函数和rest状态下的时间点扩展函数得出每个采样点处组织体相对于rest状态的吸收系数的变化量，列式如下

公式(9)中：和分别表示组织体在三个波长下task状态相比于rest状态吸收系数的变化量；△c_HbO，△c_HbR和分别表示三种物质task状态相比于rest状态浓度的变化量，解式(9)可以求出△c_HbO，△c_HbR和

再根据下述的过定方程求出消除水吸收之后修正的氧合血红蛋白和脱氧血红蛋白相对于rest状态下各自的浓度的变化量；

公式(10)中：和分别表示组织体相比于静息状态时在三个波长下的吸收系数变化量，和分别表示水在三个波长下的消光系数，表示组织体相对于静息状态时水的浓度变化量，和分别表示氧合血红蛋白在三个波长下的消光系数，和分别表示脱氧血红蛋白在三个波长下的消光系数，_correctΔc_HbO和_correctΔc_HbR分别表示消除水吸收后修正的组织体氧合血红蛋白和脱氧血红蛋白的浓度相比于静息状态时的各自的变化量。

2-3、根据如下关系式

利用rest状态消除水吸收之后修正的的氧合血红蛋白和脱氧血红蛋白的浓度和task状态消除水吸收之后修正的氧合血红蛋白和脱氧血红蛋白相对于rest状态下各自的浓度的变化量求出task状态下的氧合血红蛋白和脱氧血红蛋白的浓度；

再利用求出task状态下组织体的血氧饱和度；

2-4、求出每个采样点在task状态下相比于rest状态的血氧饱和度的变化量ΔSaO₂＝^taskSaO₂-^restSaO₂；

步骤三、根据步骤二获得的每个采样点在task状态下相比于rest状态的血氧饱和度的变化量，并作出如图4所示二维图像，用matlab软件画出16个采样点血氧饱和度变化量的二维拓扑成像。

尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。

Claims (1)

1.一种局部血氧饱和度变化量拓扑成像方法，所采用的组织体局部血氧饱和度变化量拓扑成像装置，包括光源单元、探测单元和计算机；

所述光源单元包括激光控制器(1)、三个不同波长的近红外激光器(2)、波分复用器(3)和1:9光开关(4)；三个激光器(2)产生的不同波长的近红外激光由多路单模光纤耦合到所述波分复用器(3)后由一路光纤输出；所述1:9光开关(4)用于实现在不同的光源入射位置间进行切换；

激光器(2)为近红外皮秒脉冲半导体激光器，激光器(2)的脉冲半高宽为40～80ps；三个激光器(2)的波长在650～950nm,且三个激光器(2)之间波长的差值为50～150nm；

所述探测单元包括光纤布配器(6)、9根源光纤(5)、4根探测光纤(7)、用于准直反射光的准直器(8)、对输出光强进行调节的滤光轮(9)、探测光子的四个光电倍增管PMT(10)、路由器(11)以及多维单光子计数模块TCSPC(12)；其中：

所述源光纤(5)用于传导照射到待测组织体表面不同源位置的入射光；

所述光纤布配器(6)用于布配源光纤(5)和探测光纤(7)在待测组织体表面的位置；

所述探测光纤(7)用于传导由待测组织体表面不同探测位置反射的漫反射光；

所述滤光轮(9)内装有3～6个不同衰减系数的中值密度衰减片；

所述路由器(11)用以辨别四个光电倍增管(10)探测到的多维漫反射光信号；

所述多维单光子计数模块TCSPC(12)用以记录多维漫反射光信号的时间点扩展函数；

所述激光控制器(1)、所述波分复用器(3)、所述滤光轮(9)、所述路由器(11)和所述多维单光子计数模块TCSPC(12)均连接至所述计算机(13)；

其特征在于，该方法包括以下步骤：

步骤一、数据采集，包括：

利用光纤布配器(6)将9根源光纤(5)和4根探测光纤(7)固定在待测组织体表面；其中，9根源光纤(5)形成的源点按照三行三列等距布置，相邻两个源点的距离为1～3cm；4根探测光纤(7)形成的探测点分别位于相邻4个源点的中心位置；每个源点和与其相邻的探测点之间形成有一采样点；

将待测组织体处于静息状态时记作rest状态，将待测组织体处于任务状态时记作task状态；

用波分复用器(3)耦合之后的激光束依次激励所述9根源光纤(5)形成的源点，并用滤光轮(9)中的衰减片进行衰减，然后再入射到4个光电倍增管(10)的阴极，将光信号转换成电信号，4个光电倍增管(10)的电信号经过路由器(11)辨别之后最终送达多维单光子计数模块TCSPC(12)进行单光子计数，在与每个源点相邻的探测点同时获取采样点在rest状态和task状态下的时间点扩展函数；

将待测组织体处于rest状态时的16个采样点的时间点扩展函数表示为：

${\hat{R}}_{rest}^{{\mathrm{\λ}}_i}({\mathrm{\ρ}}_j,t)---\left(1\right)$

将待测组织体处于task状态时的16个采样点的时间点扩展函数表示为：

${\hat{R}}_{task}^{{\mathrm{\λ}}_i}({\mathrm{\ρ}}_j,t)---\left(2\right)$

公式(1)和公式(2)中：i＝1,2,3；j＝1,2,…,16，λ_i分别表示上述三个波长，ρ_j分别表示上述16个采样点，rest表示组织体处于静息状态，task表示组织体处于任务状态，t为时间自变量；

最终得到96个时间点扩展函数，其中包括48个rest状态下的时间点扩展函数和48个task状态下的时间点扩展函数；

步骤二、数据处理，包括：

2-1、平面半无限空间外推边界条件下，组织体表面距源点ρ处的探测点反射光流量表达式如下：

$R(\mathrm{\ρ},t)=\frac{1}{2}{\left(4\mathrm{\π}\mathrm{\κ}\right)}^{-\frac{3}{2}}{t}^{-\frac{5}{2}}{e}^{-\left({\mathrm{\μ}}_{a}ct\right)}\[z_0{e}^{\frac{-r_1^2}{4\mathrm{\κ}t}}+(z_0+2z_b)e^{\frac{-r_2^2}{4\mathrm{\κ}t}}\]---\left(3\right)$

公式(3)中：ρ为相邻源和探测点之间的几何距离；c为光在真空中的传播速度；κ为扩散系数；μ_a为组织体的吸收系数，z₀为实际边界距离，z_b为外推边界距离，t为时间自变量；设所采用的组织体局部血氧饱和度变化量拓扑成像装置的系统响应函数为IRF(t)，实验所测得的静息状态下时间点扩展函数为IRF(t)与^restR(ρ,t)的卷积即

${\hat{R}}_{rest}^{{\mathrm{\λ}}_i}({\mathrm{\ρ}}_j,t)=IRF\left(t\right)*R_{rest}(\mathrm{\ρ},t)---\left(4\right)$

其中，^restR(ρ,t)是在rest状态下，组织体表面距原点ρ处的探测点反射光流量，对公式(4)运用最小二乘拟合算法，用在步骤一中获得的rest状态下的时间点扩展函数得出每个采样点处组织体在三个波长下的吸收系数，进一步求出组织体中氧合血红蛋白、脱氧血红蛋白和水三种物质在rest状态下各自的浓度；

再根据下述的过定方程求出消除水吸收之后修正的氧合血红蛋白和脱氧血红蛋白的浓度；

公式(5)中：和分别表示组织体处于静息状态时在三个波长下的吸收系数，和分别表示水在三个波长下的消光系数，表示组织体处于静息状态时水的浓度，和分别表示氧合血红蛋白在三个波长下的消光系数， 和分别表示脱氧血红蛋白在三个波长下的消光系数，和分别表示消除水吸收之后修正的组织体处于静息状态时氧合血红蛋白和脱氧血红蛋白的浓度；

再利用求出rest状态下组织体的血氧饱和度；

2-2、待测组织从rest状态变化到task状态期间，约化散射系数变化量△μ_s'远小于吸收系数改变量△μ_a，忽略不计，结合式(4)，有如下关系式：

${\hat{R}}_{rest}(\mathrm{\ρ},t)\·e^{-{\mathrm{\Δ\μ}}_a\·ct}={\hat{R}}_{task}(\mathrm{\ρ},t)---\left(6\right)$

对公式(6)运用最小二乘拟合算法，用在步骤一中获得的rest状态下的时间点扩展函数和task状态下的时间点扩展函数得出每个采样点处组织体在三个波长下相对于rest状态的吸收系数的变化量，进一步求出组织体中氧合血红蛋白、脱氧血红蛋白和水三种物质相对于rest状态下各自的浓度的变化量；

再根据下述的过定方程求出消除水吸收之后修正的氧合血红蛋白和脱氧血红蛋白相对于rest状态下各自的浓度的变化量；

$\left[\begin{array}{c}{\mathrm{\Δ\μ}}_a^{{\mathrm{\λ}}_1}-{\mathrm{\ϵ}}_{H_{2}O}^{{\mathrm{\λ}}_1}\·{\mathrm{\Δc}}_{H_{2}O}\\ {\mathrm{\Δ\μ}}_a^{{\mathrm{\λ}}_2}-{\mathrm{\ϵ}}_{H_{2}O}^{{\mathrm{\λ}}_2}\·{\mathrm{\Δc}}_{H_{2}O}\\ {\mathrm{\Δ\μ}}_a^{{\mathrm{\λ}}_3}-{\mathrm{\ϵ}}_{H_{2}O}^{{\mathrm{\λ}}_3}\·{\mathrm{\Δc}}_{H_{2}O}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\ϵ}}_{HbO}^{{\mathrm{\λ}}_1}{\mathrm{\ϵ}}_{HbR}^{{\mathrm{\λ}}_1}\\ {\mathrm{\ϵ}}_{HbO}^{{\mathrm{\λ}}_2}{\mathrm{\ϵ}}_{HbR}^{{\mathrm{\λ}}_2}\\ {\mathrm{\ϵ}}_{HbO}^{{\mathrm{\λ}}_3}{\mathrm{\ϵ}}_{HbR}^{{\mathrm{\λ}}_3}\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}_{correct}{c}_{HbO}\\ \mathrm{\Δ}_{correct}{c}_{HbR}\end{array}\right]---\left(7\right)$

公式(7)中：和分别表示组织体相比于静息状态时在三个波长下的吸收系数变化量，和分别表示水在三个波长下的消光系数，表示组织体相对于静息状态时水的浓度变化量，和分别表示氧合血红蛋白在三个波长下的消光系数，和分别表示脱氧血红蛋白在三个波长下的消光系数，_correctΔc_HbO和_correctΔc_HbR分别表示消除水吸收后修正的组织体氧合血红蛋白和脱氧血红蛋白的浓度相比于静息状态时的各自的变化量；

2-3、利用rest状态消除水吸收之后修正的氧合血红蛋白和脱氧血红蛋白的浓度与task状态消除水吸收之后修正的氧合血红蛋白和脱氧血红蛋白相对于rest状态下各自的浓度的变化量求出task状态下的氧合血红蛋白和脱氧血红蛋白的浓度；再利用求出task状态下组织体的血氧饱和度；

其中，表示消除水吸收之后修正的组织体处于task状态时氧合血红蛋白的浓度，表示消除水吸收之后修正的组织体处于task状态时脱氧血红蛋白的浓度；

2-4、求出每个采样点在task状态下相比于rest状态的血氧饱和度的变化量ΔSaO₂＝^taskSaO₂-^restSaO₂；

步骤三、根据步骤二获得的每个采样点在task状态下相比于rest状态的血氧饱和度的变化量，用计算机上的matlab软件画出16个采样点血氧饱和度变化量的二维拓扑成像。