一种用于近红外光脑功能研究的自引导扩散光层析成像方法

摘要

本发明公开了一种用于近红外光脑功能研究的自引导扩散光层析成像方法，包括：首先获得每个源点对应探测点位置光强；对探测点位置光强按照MLBL-OT获得探测区域吸收系数变化的二维拓扑图像并进行图像分割，实现吸收系数变化区域有效定位，生成该区域定位模板矩阵K；按照DOT重建方式采用Newton-Raphson迭代法求解，计算节点的Jacobi矩阵J，获得OT-DOT重建方程；计算得出预定探测内吸收系数变化，绘制OT-DOT重建图像。本发明可以在近红外光学模态测量方式下根据MLBL-OT定位脑功能变化区域改善DOT逆问题欠定性，从而提高重建速度的近红外脑功能自引导成像方法。

说明书

技术领域

本发明属于生物医学中近红外脑功能成像方法，尤其涉及一种近红外光脑功能研究中 自引导扩散光层析成像方法。

背景技术

人脑活动可导致灰质中血液流动和血氧代谢活动变化，改变大脑灰质中对应区域的组 织光学参数。由于人体生物组织对波长在650-900nm范围的入射光吸收较小，因此可以穿过 头皮和头骨到达大脑皮层。近红外脑功能研究利用该波段的近红外光入射光，对人体静息 状态和脑活动状态下头部预定探测区域进行反射式测量，获得两个状态下的头皮处的光场 分布。通过对灰质组织中主要吸收体的吸收系数变化进行有效重建，获得组织中氧合血红 蛋白、脱氧血红蛋白、血氧饱和度等重要信息，进而研究脑功能变化引起的大脑相关区域 的血液流动及血氧代谢情况。近红外光脑功能成像具有无创、无辐射、较高的时间分辨率、 价廉可移动等特点，为新生儿及成人脑思维或意识研究提供客观依据。

目前，在近红外脑功能成像中主要采用基于修正朗伯-比尔的拓扑成像方法 (MLBL-OT)。MLBL-OT假设源-探点之间为光学性质均匀的组织体，采样点假设位于源- 探点间的中心区域，能够获得脑血氧变化的二维拓扑图像。MLBL-OT具有计算简单、强鲁 棒性、时间分辨率高、适合实时成像的特点。但是MLBL-OT的空间分辨率受制于源-探距离 和源-探分布密度，算法中头部差分路径因子(DPF)存在个体差异和不确定性从而影响定 量性。MLBL-OT虽然获取图像空间分辨率和量化度较低，但具有对脑血氧变化区域变化实 现快速定位的优势。

为了提高近红外脑功成像的空间分辨率和定量能力，基于扩散层析成像(DOT)理论 的近红外脑成像系统获得更多关注。DOT方式在正问题中以光子输运理论为基础计算预定 探测区域光场分布，在逆问题中采用多次迭代策略对组织光学参数进行重建。与MLBL-OT 方式相比，DOT方式具有更加准确的数学模型，无需考虑个体或不同位置的DPF，并且实现 了重建区域的三维成像。DOT方式虽然提高了图像的空间分辨率和量化率，但离散化导致 测量数远小于重建参数个数，致使逆问题欠定性严重影响重建的速度。

近年来，随着大脑研究领域内多个学科进步和发展，多模态技术在脑研究中被逐步采 用并取得了一定成果。一方面，在脑结构研究中头部CT/MRI图像的分割提取已经相当成熟， 可以为脑功能研究提供大脑生理组织层面的脑结构几何信息。另一方面，在脑功能研究中 脑电、脑磁、功能核磁、近红外功能成像等多模态间相互引导、佐证和补充。但是脑功能 多模态技术也存在多测量系统同时采用所带来的仪器费用和测量操作难度增加的明显不 足。同时，脑功能多模态测量结果之间的配准和融合相当难度给算法方面带来挑战。

发明内容

考虑现有技术中近红外脑功能研究中的优势和不足，鉴于MLBL-OT的变化定位和DOT 量化度高的特点，本发明采用近红外脑功能模态的同一测量数据研究MLBL-OT引导DOT， 提供了一种用于近红外光脑功能研究的自引导扩散光层析成像方法，可以在同一光学模态 下根据MLBL-OT定位脑功能变化区域改善DOT逆问题欠定性，从而提高重建速度的自引导 成像方法。

为了解决上述技术问题，本发明提出的一种用于近红外光脑功能研究的自引导扩散光 层析成像方法，包括下列步骤：

步骤一、在头部预定探测区域按扩散层析成像探测方式、以相邻最近两点之间的距离 为10mm并按照矩形阵列排布源-探位置点，在各位置点放置光纤；探测时，其中一个位置 点为源点，其余各点为探点；首先，指定其中一位置点为源点，获得该源点对应的探点位 置的光强，依此类推，直至以每一位置点均作为源点为止，分别测量得到预定探测区域脑 活动状态的探测器光强测量结果向量Ma和人体静息状态光强测量结果向量Mr；

步骤二、头皮和头骨层的厚度为h，选择与源点距离大于等于2h的探点位置光强，按 照修正朗伯-比尔拓扑成像方法获得头部预定探测区域吸收系数变化△μ_a的二维拓扑图像；

步骤三、对二维拓扑图像进行图像分割，实现吸收系数变化区域有效定位，并生成所 述预定探测区域的由修正朗伯-比尔拓扑成像引导扩散层析成像的模板矩阵K，包括以下步 骤：

对修正朗伯-比尔拓扑成像结果以吸收系数变化△μ_a中的最小值作为阈值进行区域划分， 得到二维拓扑图像下的脑功能变化区域ROI和未发生脑功能变化区域n-ROI；

对头部预定探测区域数学模型采用有限元方法计算中三维节点坐标(x,y,z)和静息状 态模拟光强分布向量Fr；

根据头部预定探测区域内三维节点坐标(x,y,z)的(x,y)是否属于脑功能变化区域 ROI或者未发生脑功能变化区域n-ROI，生成一个模板矩阵K计算公式为：

式(1)中：k——为模板矩阵K内的元素；

nr^ALL——为模板矩阵K的行号，表示脑功能变化区域ROI内三维节点在预定探测区域 内节点中的索引号；

nr^ROI——为模板矩阵K的列号，表示脑功能变化区域ROI内三维节点在脑功能变化区 域ROI中的索引号；

步骤四、按照扩散层析成像重建方式采用Newton-Raphson迭代法求解，计算预定探测 区域内节点的Jacobi矩阵J，根据矩阵K和J生成由修正朗伯-比尔拓扑成像引导扩散层析 成像OT-DOT重建方程；

根据求解组织体内三维光学参数分布时散层析成像逆问题数学方程表示OT-DOT方式 计算公式为：

$\frac{\mathrm{Ma}·\mathrm{Fr}}{\mathrm{Mr}}-F\left({\mu }_a\right)=J·\mathrm{K\Delta }{{\mu }_a}^{\mathrm{ROI}}---\left(2\right)$

式(2)中：F(μ_a)——为通过正问题求解的边界处探测点光强向量；△μ_a^ROI——为经过 模板矩阵筛选后ROI区域内吸收系数变化；

步骤五、对OT-DOT重建方程进行计算得出预定探测区域内吸收系数变化，并且利用 MATLAB软件绘制OT-DOT重建图像。

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

本发明成像方法是在近红外光模态测量方式下采用同一测量数据，在DOT重建之前采 用MLBL-OT生成模板矩阵，将模板矩阵应用于DOT重建过程实现近红外功能成像的自引 导，从而降低DOT逆问题的欠定性，提高了重建速度。

附图说明

图1(a)是预定探测区域静息状态下平板模型，其中，黑点所示为源-探点；

图1(b)是预定探测区域脑活动状态下平板模型，其中，圆柱为吸收系数变化区域；

图2是本发明中MLBL-OT方式成像；

图3是本发明中MLBL-OT图像分割后吸收区域定位图；

图4(a)是采用传统的DOT重建图像；

图4(b)是利用本发明MLBL-OT引导DOT重建图像，其中，虚线为吸收变化区域真实 位置。

具体实施方式

下面结合附图和具体实施例对本发明技术方案作进一步详细描述，所描述的具体实施 例仅仅对本发明进行解释说明，并不用以限制本发明。

本发明用于近红外光脑功能研究的自引导扩散光层析成像方法，具体步骤如下：

对于近红外脑功能成像的预定探测区域，可以合理简化为二层平板模型上层为头皮和 头骨层10mm，下层为灰质层20mm，如图1(a)所示，在头部预定探测区域按DOT探测 方式、源-探分布采用DOT高密探测方式两点间最短距离为10mm，并按照矩形阵列排布源 -探位置点，本实施例是4×4矩形阵列，在各位置点放置光纤，探测时，当1个位置布置光 源时，其他15个位置均布置探测器，获得该源点对应的探点位置的光强，依此类推，直至 以每一位置点均作为源点为止，然后分别测量得到预定探测区域脑活动状态的探测器光强 测量结果向量Ma和人体静息状态光强测量结果向量Mr；本实施例中是采用数学模型模拟脑 活动状态引起的脑功能变化区域(ROI)，如图1(b)中的圆柱区域所示，设圆柱区域半径 为5mm高为10mm，上表面圆心位于平板下层上表面中心。通过对图1(a)和图1(b)模 型进行空间剖分获得所有参与重建节点坐标，并且对DOT正问题进行有限元计算(FEM)， 得到Mr，Ma。在计算过程中静息状态下模型各剖分节上层吸收和约化散射系数向量为μ_a^T和 μ'_s^T，下层吸收和约化散射系数向量为μ_a^B和μ'_s^B。通过参考前人文献中人脑切片光学测量值 得到(见：文献Yaroslavsky,A.N.,Schulze,P.C.,Yaroslavsky,I.V.,Schober,R.,Ulrich,F.,and  Schwarzmaier,H.J,“Optical properties of selected native and coagulated human brain tissues in  vitro in the visible and near infrared spectral range,”Physics in medicine and biology,47(12), 2059(2002).)。

选择与源点距离大于等于20mm的探点位置光强，按照MLBL-OT方法获得头部预定 探测区域吸收系数变化△μ_a的二维拓扑图像；具体过程如下：

基于MLBL-OT作为光在组织内扩散的一阶近似，用来描述光能在高散射组织内的衰减 过程，计算公式为：

$\mathrm{\Delta OD}=-\ln \frac{\mathrm{Ma}}{\mathrm{Mr}}=\Delta {\mu }_a\mathrm{BL}$

式中：△OD——为由静息状态到脑活动状态探测位置光强变化向量；

△μ_a——为吸收系数变化向量；

B,L——分别为差分路径因子(DPF)和光源与探测器之间距离；

当采用MLBL-OT计算预定探测区域内吸收系数变化时，采样点认为处于源点和探测点 中间的区域，并假设源-探间的组织是均匀的，获得的二维拓扑图像，图2示出的是采样点 下各层吸收系数变化的平均效应。因此得到的吸收系数变化结果会远小于真实值；

对二维拓扑图像进行图像分割，实现吸收系数变化区域有效定位，并生成所述预定探 测区域的由OT-DOT的模板矩阵K，包括以下步骤：

对MLBL-OT结果以吸收系数变化△μ_a中的最小值作为阈值进行区域划分，得到二维拓 扑图像下的ROI和未发生脑功能变化区域(n-ROI)；将ROI区域记为“区域1”区域内的 所有像素点值记为1，n-ROI区域记为“区域0”所有像素点值记为0，X-Y平面上的二维 阈值区域划分结果，如图3所示。对头部预定探测区域数学模型采用有限元方法计算中三 维节点坐标(x,y,z)和静息状态模拟光强分布向量Fr；根据头部预定探测区域内三维节点 坐标(x,y,z)的(x,y)是否属于ROI或者n-ROI，生成一个模板矩阵K计算公式为：

式(1)中：k——为模板矩阵K内的元素；

nr^ALL——为模板矩阵K的行号，表示ROI内三维节点在预定探测区域内节点中的索引 号；

nr^ROI——为模板矩阵K的列号，表示ROI内三维节点在ROI中的索引号；

除行列索引号为(nr^ALL,nr^ROI)外，K中的其他元素值都为0。模板矩阵K按照二维ROI和 n-ROI区域实现了对三维预定探测区域内所有FEM节点的区域划分。

由于MLBL-OT方式成像的分辨率较低，因此按照MLBL-OT得到的ROI范围会略大 于真实ROI，并且模板矩阵K是由二维平面区域到三维立体区域的转换，三维ROI区域内 必然包括一些未发生吸收系数变化的节点。但是，模板矩阵K可以达到按照MLBL-OT成 像的分辨率筛选出预定探测区域内未发生吸收系数变化的节点。

按照传统三维DOT重建方式，预定探测区域内所有节点都将参与重建。但DOT重建 中正问题和逆问题的计算中需要对成像空间进行节点离散，三维重建将导致测量数据远远 少于重建参数个数，因此逆问题的欠定性将十分严重，影响重建速度。为了有效的减少参 与重建的节点个数从而改善欠定性，在OT-DOT方式下：将模板矩阵K与DOT方式重建中 经过Newton-Raphson迭代法求解得到的Jacobi矩阵做矩阵乘法运算。这样所得到的矩阵中 包括真实ROI区域和小部分n-ROI区域的重建节点，对传统三维DOT重建中所有节点均参 与重建方式进行优化，达到减少n-ROI区域的重建节点数的目的；根据求解组织体内三维 光学参数分布时DOT逆问题数学方程表示OT-DOT方式计算公式为：

$\frac{\mathrm{Ma}·\mathrm{Fr}}{\mathrm{Mr}}-F\left({\mu }_a\right)=J·\mathrm{K\Delta }{{\mu }_a}^{\mathrm{ROI}}$

式中：F(μ_a)——为通过正问题求解的边界处探测点光强向量；

J——为传统DOT方式下经过Newton-Raphson迭代法计算得到的Jacobi矩阵；

△μ_a^ROI——为经过模板矩阵筛选后ROI区域内吸收系数变化；

虽然OT-DOT方式逆问题中n-ROI区域重建节点较传统DOT方式大幅度减少，但该方 程仍然是欠定的，因此需要对该方程采取代数重建技术(ART)的正规化方法进行求解。这 里分别给出的DOT和OT-DOT方式对ROI重建模拟成像结果测量值由正问题的FEM计算 产生，背景光学μ_a^T设为0.015mm^-1，μ'_s^T设为2.0mm^-1，μ_a^B设为0.02mm^-1，μ'_s^B设为1.2mm^-1， ROI吸收系数设为0.03mm^-1。MLBL-OT计算中B＝6.3，L＝20。由于FEM计算得到离散点 的重建值，在做图过程中采用了MATLAB的插值函数，绘制z＝12.5mm处X-Y截面重建图， y＝30mm处X-Z截面重建图如图4(a)和图4(b)所示。在迭代过程中，迭代终止条件设 为迭代误差变化小于0.01。在达到相同的迭代终止条件情况下DOT方式迭代次数为28次， OT-DOT迭代次数为9次，表明OT-DOT方式较DOT方式明显提高了重建速度。

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