基于蒙特卡罗的脑血氧检测光源—探测器最优距离研究

曹 焱，邢志明，赵 斌，董祥美，史以珏，高秀敏

（1.上海理工大学 光电信息与计算机工程学院，上海 200093；2.上海交通大学医学院附属瑞金医院，上海 200025）

0 引言

大脑是人体的一个重要器官，像一个中枢系统一样指挥着人们的日常行为与活动。脑部的血流量约占心脏血液输出量的15%，耗氧量约占全身耗氧量的20%［1］，所以大脑对于缺血缺氧是十分敏感的［2］。临床上对脑血氧的监测可以帮助医生在手术中实时监测患者健康状况［3-4］，以及进行一些脑部疾病的诊断［5-6］。因此，用于临床监测和诊断的血氧监测仪器可为医生提供很大的帮助。

利用功能性近红外光谱技术［7］可以很容易地检测到脑功能的血流动力学信息［8］。在实际的脑血氧测量中，被检测者往往会因人体的生理活动（如心动周期、呼吸、血压变化等）引入干扰，从而给脑血氧信号检测造成了一定影响，这些生理干扰可以采用自适应滤波法［9-10］去除。自适应滤波法对于抑制特定的生理干扰有比较好的效果，但是需要借助额外的设备获取参考信号。通常使用单个探测器进行脑血氧信号检测时会面临以上问题，因为仅有一个探测器无法获取参考信号。为解决单个探测器检测所面临的问题，有学者提出了多距离测量方法，即使用双探测器测量脑血氧信号［11-12］。一个探测器检测深层组织中带有生理活动干扰的血氧信号，另一个探测器检测浅层信号作为参考信号，对两个探测器检测的信号作差分处理［13］，获得纯度更高的局部信号用于后续计算处理。此外，临床和科研中光源—探测器（Source-Detector，SD）距离的选择大多根据经验公式［14］，即探测器能检测到的最大组织深度是光源与探测器距离的一半［15］，或者参考市面上现有的脑血氧检测产品来选择光源—探测器之间的距离［16］。这些方法选择的光源与探测器距离是固定的，但由于颅脑各部分组织结构针对不同的年龄、性别、个体等会存在差异，如正常人和头部有损伤的病人（如硬膜脑血肿［17］）等，经验公式可能无法很好地适用于所有人群的脑血氧检测，无法满足光源—探测器之间的最佳距离选择需求。

为解决上述脑血氧检测中SD 距离选取时遇到的问题，本文根据多距离测量方法要求，给光源—探测器选择一组恰当的距离。光源放置于固定位置时，探测器需一个放置于近处，用于检测头皮和颅骨处的血液信号作为参考信号；另一个放置于远处，用于检测灰质处受到生理活动干扰的血液信号。再结合蒙特卡罗方法，从光在组织中的传输特性着手，同时考虑到颅脑组织的结构差异，选定一组颅脑组织的特性参数，建立5 层颅脑组织平板模型进行仿真。通过统计大量光子在组织中的传输情况，并对仿真得到的数据进行分析与处理，可以很好地克服在颅脑组织存在差异时应用经验公式的缺陷，从而选择出符合该组颅脑组织参数的最优的光源—探测器距离。

1 蒙特卡罗模型与方法

1.1 蒙特卡罗原理

生物组织在光学上可近似看作混沌介质，当光照射在生物组织上时，会发生反射、吸收和散射等现象，光在生物组织中的宏观传播行为已不能用基本的朗伯—比尔定律来描述。近些年得益于计算机技术的快速发展，蒙特卡罗（Monte Carlo）方法在各个物理仿真实验中使用得越来越多，尤其在生物医学、光学等领域［18］。Monte Carlo 模拟［19-20］是一种基于概率统计的数学物理模拟方法，通过相应的概率模型或随机过程计算，模拟粒子在组织中传输时受到的影响［21］。通过对大量光子在组织中传播路径的仿真，追踪光子运动的轨迹点，可得到光在组织内部和表面的统计分布规律。

1.2 光在组织中传输的蒙特卡罗模拟方法

本文利用多层蒙特卡罗（Monte Carlo Multi-Layer，MCML）程序进行仿真，MCML 程序是Wang 等［19］用标准C 语言开发的程序包，用于模拟分层组织的光子传输［22］。蒙特卡罗程序的基本流程如下［23-24］：

（1）确定起始跟踪点和初始化光子（权重、方向、位置）后，在组织表面发射光子，并为其分配随机步长。

（2）移动步长s 后，更新光子位置，确定光子下一次的碰撞位置。

（3）光子移动步长后会随机进行吸收和散射，散射作用会改变光子的前进方向，吸收作用会降低光子权值。根据当前位置计算由于散射造成的方向偏转，更新传播方向余弦。计算光子在当前位置由于吸收造成的权重衰减，更新权重值，并对光子权重是否低于阈值进行判断。如果低于阈值，或逸出组织上下表面时，采用轮盘赌来决定光子存活。判定为死亡的，结束当前光子的追踪。

（4）对光子是否碰撞边界进行判断，根据判断结果进行相应的处理。

（5）对于未逃逸出组织的光子继续进行步骤（2）—（4），直到程序结束。

蒙特卡罗模拟流程如图1所示。

Fig.1 Flow of Monte Carlo simulation图1 蒙特卡罗模拟流程

1.3 光在脑组织中传输路径的仿真研究

研究光在颅脑组织中的传输路径，需对每个光子的传输轨迹进行追踪，但Wang 原有的程序并不能实现此功能。因此，需要在Wang 的仿真程序中进行优化，以获得光子在组织中的传输路径。由于光子在组织内传输是随机、无序的，所以采用数据链表来保存光子的传输路径。光子在组织中传输时的平均自由路程很短，在运算时保存的信息会占用很大的计算机运行内存，故用d来判断光子运动距离，优化保存的信息。

式中，d表示光子传输时相邻两个位置的距离，（xk，yk，zk）表示光子当前的位置，更新后的位置坐标为（x，y，z）。只有当d值大于某一阈值时，才将当前位置点的坐标记录到链表中［25］。蒙特卡罗仿真时需要依靠大量光子才能得到较为准确的结果，在实际仿真时，考虑到计算机的内存空间使用情况，故只关注探测器处接收到的光子，对这些被接收到的光子路径进行保存，如图2所示。

Fig.2 Schematic diagram of light source and detector placement图2 光源与探测器位置示意图

光源和探测器放置的位置如图2（a）所示，X 轴为光源和探测器所在的水平方向，Y 轴为头皮表层的竖直方向，Z轴为脑组织深度方向。图2（b）中光子入射点（0，0，0）与探测器的距离为SD，探测器与X 轴的夹角为θ。一定数量的光子从原点入射，在组织内进行随机传输运动，部分光子可到达接收面并被接收。为便于后续处理，将探测器接收面看成边长为Δs 的正方形。只有成功到达接收面的光子，才记录其在组织内的传输路径，其他光子则不予记录。

仿真程序中所使用的五层颅脑模型光学参数如表1所示［26-28］。表中各参数的定义如下：n是折射率，g是各向异性因子，μa/cm-1是吸收系数，μs/cm-1是散射系数，t/cm 是各颅脑层的组织厚度。

Table 1 Optical parameters of the five-layer cranial model表1 五层颅脑模型光学参数

2 蒙特卡罗仿真与讨论

为便于处理，本文将探测器与光子入射点设置为同一水平线，即设置图2（b）中的θ值为0°，接收面的边长Δs 为1mm，从原点（0，0，0）入射的光子数为107个，d的阈值为1mm。设置完成后，选择不同的光源—探测器（SD）距离值运行MCML 程序。值得一提的是，由于此仿真需要对大量数据进行统计分析，故对电脑内存要求较高。仿真硬件条件为：Intel Core i9-12900处理器，128GB 的RAM 内存。

2.1 光子传输路径

按如上条件设置后，可得到不同SD 距离下的光子传输路径数据，选取SD=1cm 和3.5cm 的数据进行处理，绘制光子在组织内的传输路径图，如图3所示。

Fig.3 Photon transport paths in tissue图3 光子在组织内传输路径

图3（a）、（b）分别是SD=1cm 和SD=3.5cm 的传输路径三维图，中心区域是一团聚集在一起的密集路径，边缘区域的路径逐渐稀疏。由图3（a）和（b）可明显看到，SD=3.5cm 时的传输路径核心区域比SD=1cm 时大一些。

传输路径图反映了光子在颅脑组织内传输的部分情况，而确定光源—探测器（SD）间的最优距离，还需要借助其他参数进行判断。本文通过研究光子传输时的最大可达深度和散射距离，作为选择光源—探测器最优距离的参考。

2.2 光子最大可达深度

探测深度是脑血氧检测中的一个重要参数，光子传输路径的最大可达深度一定程度上反映了光子在颅脑内可探测的深度。根据每条光子传输路径上Z 轴坐标的极大值，绘制出如图4 所示的不同SD 距离下的光子最大可达深度概率分布图。图中横坐标是光子传输的深度值，纵坐标是光子最大可达深度的概率值。可以看出，不同SD 距离对应的概率分布大致是正态分布。图4（a）稍有不同，左侧区域的概率值比右侧区域密集些，主要是由于SD 的距离过小，光子多在浅层区域传输，传输到深层组织的光子相对较少。这也说明当SD<2cm 时，只有相对较少的光子可到达深层组织，这对于选择近处探测器与光源之间的距离有很重要的参考价值。

Fig.4 Probability distribution of maximum reachable depth for 300 photonic transmission paths at different SD distances图4 不同SD距离下300条光子传输路径的最大可达深度概率分布

为便于分析，计算不同的SD 距离下光子最大可达深度在不同组织层的占比情况。如图5 所示，横坐标是不同的SD 距离，纵坐标是各颅脑组织在不同SD 距离时的最大深度占比情况。可看到随着SD 距离的增大，头皮层和颅骨层的占比不断降低；脑脊液层的占比非常小，可忽略不计；灰质层的占比基本大于15%。头皮和颅骨层占比在SD=1cm 处取得最大值23.78%。SD=2cm～4.5cm 时，光子最大可达深度在灰质层的占比于SD=3.5cm 处取得最大值19.59%

2.3 光子传输路径在不同层的散射距离

仅采用最大可达深度描述光子对深层颅脑组织探测的情况还不够客观、全面，最大可达深度只能了解光子是否到达某一深度值，其他信息不得而知。

光子在颅脑组织中传输时，散射距离也是一个重要参数。因此，计算不同SD 距离下光子传输路径在各层（取每一层深度为0.1cm）散射距离占路径总长的比例，并绘制散射距离占比变化趋势图，如图6所示。

Fig.6 Scattering distance variation of photon at different depths under different SD distances图6 不同SD距离下光子在不同深度的散射距离变化

图6 横坐标表示颅脑组织的深度，纵坐标表示光子在各层的散射距离占传输路径总长的比值。随着深度增加，光子在深层组织处的散射距离所占比重不断减少，表明光子在深层组织传输时的散射距离短。当深度达到3.5cm时，光子的散射距离所占比重基本为0，说明光子很难在大于3.5cm 的深层组织传输。通过图6 可以了解光子散射距离的整体趋势，但还不足以了解光子在不同颅脑层的情况。故绘制光子在不同颅脑层的散射距离占总散射距离的比重，如图7所示。

Fig.7 The proportion of scattering distance of different brain layers to the tatal scattering distance图7 不同颅脑层散射距离占总散射距离的比重

在图7 中，横坐标是不同的SD 距离，纵坐标是光子分别在头皮和颅骨、脑脊液、灰质层的散射距离占总散射距离的比重。可明显看到，光子的散射距离在头皮和颅骨层所占比重很大，基本都在25%左右。由于脑脊液层的特殊光学性质（低吸收、低散射），故光子在脑脊液层的散射距离比重都很小。光子在灰质层的散射距离比重不等，与SD 距离相关，当SD=2～4.5cm 时，灰质层的散射距离比重在SD=4.5cm 处取最大值10.22%。这是因为随着SD 距离的不断增大，与较小的SD 相比，有更多光子可以到达更深层的组织进行传输。故随着SD 的增大，灰质层的散射距离比重也在增大。

3 结果分析

根据多距离测量方法对SD 距离的选择要求，近距离探测器D1 检测头皮和颅骨层的血液信息，远距离探测器D2 检测脑灰质层的血液信息。在对距离进行选择时，需要同时考虑3 个方面：一是光源在传输时可达的最大深度；二是光源在组织中传输时携带信息多少；三是探测器实际接收到的信号强度。

有效深度比（Effective Depth Ratio，EDR）是探测器接收到来自有效深度层的光子与所有组织层光子的比值。将头皮和颅骨层看作第一层，将脑脊液、灰质和白质层看作第二层，则在式（2）中，第一层（最大深度）为最大可达深度在头皮和颅骨层，且并未穿过颅骨层的光子，第二层（最大深度）为穿透颅骨层的光子。

根据上文所述，计算出不同SD 距离下的EDR 值并绘制图像，如图8（a）所示。图8（a）中，横坐标为不同的SD 距离，纵坐标为EDR 值。从图8（a）中可以看出，随着SD 距离的增大，有效深度比EDR 逐渐减小，说明随着SD 的增大，传输到深层组织的光子增多。EDR 值越大，说明SD 的距离越小，光子传输时位于头皮和颅骨层的光子相对较多；反之，SD 距离越大，位于灰质层的光子相对较多。故可将EDR 值作为选择光源与探测器间距离的重要依据。

Fig.8 Effective depth ratio and grey matter layer information ratio at different SD distances图8 不同SD距离下的有效深度比与灰质层信息比

此外，光子在组织中传输时，携带的信息量大小与光子散射距离有很大关系。光子传输的EDR 值并不能完全反映探测器接收到来自灰质层的信息量大小，故对图7 中的结果作进一步处理与分析。计算在不同SD 距离下光子传输时的灰质层信息比（Gray Matter Layer Information Ratio，GMLIR），表示为：

式中，灰质层（散射距离）为光子在灰质层的散射距离与总散射距离的比；头皮、颅骨层（散射距离）为光子在头皮和颅脑层的散射距离与总散射距离的比。绘制如图8（b）所示的图像，从图中可明显看出，随着SD 距离的增大，灰质层信息不断增多。在脑血氧检测中，灰质层信息比的增大意味着探测器探测到的有效信息增多，这对于实际检测是很重要的。故GMLIR 值也作为选择光源与探测器距离的一个重要依据。

在实际选择光源与探测器距离时，还需要考虑探测器接收到的信号强度。本文在研究中将光电探测器贴近头皮表层进行必要的遮光处理后，分别使用760nm 和850nm的激光二极管（LD）在不同的SD 距离下照射头皮表层，记录探测器的电压值。如图9 所示，当不使用LD 照射时，探测器的电压均值为3.67mV（参考值），而当使用LD 时，LD和探测器间的距离与信号强度成反比。SD 距离较小时，如SD=1cm，探测器接收到760nm 和850nm 信号的电压均值大于400mV。随着SD 距离的增大，信号强度开始明显变小，说明SD 距离对于信号强度的影响非常大。探测器与LD 相距4cm 的位置处，探测器的电压均值非常接近参考值，意味着探测器接收到的信号已很微弱。在确定SD距离时，需要考虑到探测器接收信号强弱这一方面的影响，故在选择SD 时，将不再考虑SD>4cm 的情况。

Fig.9 Signal intensity received by the detector at different wavelengths with different SD distances图9 不同SD距离时不同波长下探测器接收到的信号强度

根据以上分析，EDR 与GMLIR 是选择SD 距离的重要依据。近距离探测器D1 主要检测头皮和颅骨层信息，需要EDR 值较大，且不检测灰质层信息，因此需要GMLIR 值较小；远距离探测器D2 主要检测灰质层信息，其中含有头皮和颅骨层的干扰信息少，故需要EDR 值较小，GMLIR 值较大。在选取光源与D2 的距离时，还需考虑探测器实际接收信号的强度。综上所述，本文选择的一组最优光源与探测器距离为SD1=1cm，SD2=3.5cm。

4 结语

本文通过选取一组颅脑组织的特性参数，利用蒙特卡罗方法对颅脑组织模型进行仿真研究，获得光子在头部组织的传输路径，进而分析光子传输时的最大可达深度和散射距离，以及有效深度比EDR 和灰质层信息比GMLIR。根据多距离测量方法的要求，考虑了EDR、GMLIR 和探测器实际接收到的信号强度，确定一组适合选定颅脑组织参数的最优光源—探测器距离，SD1=1cm 和SD2=3.5cm。该组距离不仅可以最大程度地检测灰质层的血氧信息，而且可以最大程度地处理来自头皮和颅骨层的干扰，从而获得纯度更高的局部脑血氧信号。利用多距离测量方法设计的脑血氧检测装置，可以选择此组最优距离进行传感器探头设计。本研究结果对脑血氧检测装置的设计，以及利用近红外光谱技术检测脑血氧的精度优化具有一定指导意义。