基于并行化逆向追踪技术的浅层组织拉曼蒙特卡罗模型

高 峰 ，吴雨桐，贾梦宇，王兵元，赵会娟

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

基于拉曼光谱技术的扩散光学层析成像技术[1]旨在将光谱技术与扩散光学层析成像[2](diffuse opti-cal tomography，DOT)技术有机融合，在构建复杂生物组织的三维立体空间图像的同时，利用拉曼光谱的生化指纹特征实现对生物组织的成分分析和识别，进而更加准确地进行早期病变诊断[3].而建立精确高效的拉曼光子输运模型是实现该目标的关键步骤.

尽管基于有限元法(finite element method，FEM)实现的扩散近似(diffusion approximation，DA)模型在 DOT已有广泛应用，计算效率高，但仅适用于远源场光子输运的描述[4].扩散拉曼层析成像(diffuse Raman tomography，DRT)技术所涉及的拉曼散射主要发生在浅层生物组织中，FEM-DA模型并不适用[5].而蒙特卡罗(Monte-Carlo，MC)方法作为辐射传输方程的随机统计解法，能精确地描述任意区域和光学参数分布下组织中的光子传输行为，因而是生物组织拉曼光谱分析和成像中适用的光子输运模型技术[6].

传统拉曼蒙特卡罗(Raman Monte-Carlo，RMC)方法计算原理简单、模拟准确，被广泛用于拉曼光谱技术的建模过程.Wilson等[7]应用传统RMC方法获取并分析了给定源探分布下的多波段光子计数，从而有效区分了该模型中的特定生化物质成分.Keller等[8]应用传统 RMC方法建立了基于空间分辨的多层组织模型，通过模拟分析，实现了对肿瘤边界位置的识别.但传统RMC模型往往需在多个散射波段下对全体积域内的剖分节点逐次激发模拟，计算量大且耗时长.Hokr等[9]对传统 RMC方法进行基于GPU技术的并行化加速处理，可有效缩短模拟耗时，尽管如此，该方法仍不足以有效延伸至多波段情形；Demers等运用NIRFAST软件建立正向模型，实现了全三维活体组织的DRT重建，与RMC相比，该软件采用 FEM-DA求解，虽计算速度快，但难以精确描述 DRT所需的近源场拉曼光子传输过程，且该软件高度集成封装，扩展其功能和算法存在着诸多局限[10-11].因此，建立一种开源且高效的拉曼光子输运模拟算法，对发展多波段下复杂组织体的 DRT图像重建方法具有重要意义.

本文提出了一种基于并行化逆向追踪技术的浅层组织拉曼蒙特卡罗模型(GPU-Reverse RMC)，以实现 DRT算法的快速建模.根据格林函数的互易性质，通过调换传统RMC模型中给定源探对位置，在探测器位置注入大量散射光子并逆向追踪，从而获取与传统方法等效的探测结果[12].继而应用GPU技术对大量光子的逆向追踪过程进行单指令多线程的并行化处理，可进一步降低该模型的计算耗时[13].本文分别建立了三维均匀和非均匀数值模型，对不同 RMC方法进行数值模拟实验，又应用 GPUReverse RMC方法快速获取了 89、178和 356个波段下的复杂组织拉曼光谱，模拟结果表明，此方法适用于建立多波段拉曼光子输运模型，为 DRT图像重建工作提供一种快速有效的建模方法.

1 模型介绍

1.1 RMC传统模型

传统RMC模拟方法分为2步.

步骤 1 在组织体表面 rs处沿法向量方向注入波长为λx的激发光 Qx(rs,λx)，采用MC法模拟获取激发光在组织体内的光子密度分布Φx(rs,λx).

步骤2 将模拟谱段离散为W个波段，在给定模拟波段的光学参数下对组织内r处的二次弥向光源进行二次模拟，获取表面rd处的拉曼边界光流量，如图1(a)所示.

图1 两种模型对比Fig.1 Comparisons between the two models

忽略实际测量时探头数值孔径及其他实验因素，传统RMC方法原理为

式中：为在边界探测点处对整个空间角 Ω积分所得；c为光速；为格林函数，表示光源位于r处激发探测器于rd处可接收到的拉曼光流量；二次激发光源由给定模拟波段下组织内r处的拉曼散射系数和激发光光子密度所决定，即

式中：为 RMC第 1步模拟所得；

通过对所有离散波段进行RMC模拟，最终获取组织体表面多个源探对下的拉曼模拟光谱数据

1.2 基于GPU逆向追踪技术的RMC模型

为克服传统 RMC模型(图 1(a))计算量大且耗时长的不足，本文发展了基于格林互易性质的逆向RMC模型(reverse RMC，rRMC).已知式(1)中拉曼激发光源rQ与边界光流量Γr满足格林互异性 质[12]，即

式中：Γr(r,- n)为内部弥向光源 Qr(r,s)沿 s方向激发，在边界 rd处沿法向量方向获取的拉曼边界光流量；当调换源探位置，将光源,- n)置于 rd处沿-n方向激发，在r处可收集到沿-s方向的等效光流量数据′(r,- s).

考虑到光子的全反射效应对于边界光子探测效率的影响，为保证边界入射光源与内部弥向光源作用等效，需将光源,- n)设置为张角为ΔΩ的均匀锥束光束，ΔΩ为光子在边界处全反射角2,π立体角积分[12]. 由于此时入射光束面积发生改变，光源,- n)的作用范围缩小，为获取等效的模拟结果，应扩大光源强度至相应倍数，即

此时，将光源,- n)置于 rd处沿-n方向激发，所求等效边界光流量可看作是对整个组织体光流量的积分，即

对注入的大量拉曼光子的输运过程进行逆向追踪并记录，获取等效的散射光流量.考虑到散射光子的逆向追踪过程彼此独立，可对其进行基于 GPU的单指令多线程的并行化处理，从而进一步提高模拟计算效率，该方法的流程如图2所示.

图2 GPU-Reverse RMC模拟流程Fig.2 Flow chart of the GPU-Reverse RMC model

2 模拟实验

为验证GPU-Reverse RMC方法的有效性，进行如下模拟实验.

2.1 仿体模型

建立如图 3(a)所示 1,000,mm3的正方体数值模型，在距离上表面 1,mm的中心位置处植入体积为1,mm3的正方体异质体.以上表面中心点为三维坐标系原点，将光源置于原点位置，探测器 D1～D20沿上表面对角线均匀排布，如图3(b)所示.

选择 I型蛋白质、甘油酸脂、脱氧核糖核苷酸、β-胡萝卜素和钙化物 5种纯净物质，建立背景组织及含异质体组织的生化成分模型.

拉曼散射强度由拉曼散射系数srμ决定，约为入射光强的10-9.计算该模型的拉曼散射系数为

式中：m表示组织中的不同物质成分；(i=1,2,… ,W)表示相对拉曼散射截面，由各物质归一化拉曼基准谱给出，表示单位浓度下发生拉曼散射的概率；CRm为相对于各物质单位浓度下的相对浓度.背景及异质体中生化成分的相对浓度 CRm及拉曼散射系数srμ由表1给出[15].

图3 仿体结构Fig.3 Phantom model

表1 背景组织和异质体的生化成分的浓度占比Tab.1 RCs of the background and target constituents

实验中往往为避免荧光干扰，选择激发光波长为785,nm，拉曼散射波段范围为827～914,nm.考虑到此波段下吸收、散射光学参数的改变对于拉曼散射结果的影响甚微，模拟实验中设置散射过程为840,nm下的组织光学参数，模拟实验选择近似皮肤组织的光学参数，如表2所示[8].

表2 背景及异质体光学参数Tab.2 Optical properties of the background and the target

表2中axμ、sxμ分别代表激发光波长下吸收、散射系数；amμ、smμ分别代表散射波长下吸收、散射系数；n为光在组织体表面的折射率；g为各项异性系数，表示散射的定向性[9].

2.2 均匀仿体模拟实验

本研究计算平台所用 GPU显卡的型号是Nvidia GTX 980.设置模拟光源为 1,mW 的连续准直光源，注入光子数为 1014,s-1.模型离散为 20×20×20个立方体体元.分别使用传统 RMC(Traditional RMC，tRMC)方法和逆向 RMC(Reverse RMC，rRMC)方法对均匀背景仿体进行数值实验模拟，得到探测器 D1～D20探测点处的拉曼边界光流量，并计算其相对误差(relative error，RE)为

由图 4可知，当源探距离小于 4,mm 时，逆向RMC方法模拟结果的相对误差近似为零，说明此时的模拟结果可有效反映拉曼光的扩散效应；当源探距离大于4,mm时，获取远离光源部分的模拟数据时出现较为明显的相对误差，这是由于此时探测信号的光子量级低于 102，模拟产生误差在较弱探测信号强度下较为明显.

图4 传统 RMC方法和逆向 RMC方法模拟结果比较(D1～D20)Fig.4 Comparisons of simulating result by the tRMC method and the rRMC method(D1～D20)

逆向 RMC方法的模拟信号强度比传统方法略有偏低，这是由于该方法仅通过在边界探测点处对光源进行激发模拟并统计求和，与传统的逐点模拟方法相比，在远离表面光源的部分存在着一定程度的能量损失，因而统计求和后其结果往往略有偏低.

考虑到只有基于足够大量随机实验统计求解时，模拟结果才能有效反映真实情况，因此当探测信号的光子量级低于102时，其结果易受噪声等因素干扰，难以有效还原光子在组织中的扩散情况，参考意义较小；针对此情况，可考虑通过延长模拟积分时间以提高模拟信号的强度，从而保证信号的准确性.

对第2.1节中的数值模型分别采用tRMC方法、GPU-tRMC方法和 GPU-rRMC方法进行计算机数值模拟，注入激发光光子数为 107.将单一波段、89、178、356个波段的计算耗时估算和计算机模拟结果列于表3.

表3 3种RMC方法计算时间对比(20个探测点)Tab.3 Comparisons of the execution time among the three methods(20 detectors)

当获取单一波段下的边界光流量时，tRMC方法和GPU-tRMC方法需对组织内全部离散体元进行逐点 MC模拟，共剖分 8,000个体元，因此需要进行8,001次 MC模拟(含激发光模拟)，前者估算耗时2.7×107,min，约 450,000,h，显然无法实现，后者可通过GPU加速将计算时间有效缩短至2.5,h；本文提出的GPU-rRMC方法仅需在探测位置处共进行MC模拟，共计 21次(含激发光模拟)，计算时间可有效低至0.42,min.

当获取多波段下的模拟光谱数据时，tRMC方法由于计算量过大，只能估算其模拟耗时，89个波段下约耗时 2.4×109,min，显然无法实际应用；即使将tRMC方法经GPU加速，模拟89个波段下的光谱数据仍需耗时 1.4×104,min；采用 GPU-rRMC方法获取波段数为 89、178、356时 3种波段下的光谱数据分别耗时 50,min、85,min和 180,min，此时计算时长合理，适于实际应用.

分析上述结果的主要原因为：①tRMC方法需对组织内各离散节点逐点进行模拟，计算次数远远大于 rRMC方法，因此造成该方法下模拟耗时严重；②GPU-tRMC方法可有效提高 tRMC计算速度约105倍，但采用该方法依然需进行大量重复模拟过程，因而难以有效延伸至多波段的情形；③GPU-rRMC方法从模拟算法和实现手段两方面提高计算效率，并在模拟误差可接受范围内获取有效的模拟数据.因此该方式是获取组织拉曼光谱适用的方法.

采用GPU-rRMC方法获取3种波段下的模拟光谱，对比模拟结果以显示不同模拟波段对结果准确度的影响.获取 D10探测点的模拟结果，如图 5所示.

根据图 5可知，波段数为 89的模拟光谱在920,cm-1、1,003,cm-1、1,270～1,315,cm-1波段处出现部分光谱信息丢失现象，在 1,450,cm-1、1,521,cm-1、1,668,cm-1处的谱峰强度略低于基准谱；模拟波段数增加至178时，大部分波段下的谱峰强度得以还原，但部分特征谱峰仍未清晰反映出来；波段数为 356的光谱结果中可以清晰看到其光谱结果与基准谱较为 吻 合 ，在 845,cm-1、920,cm-1、1,245,cm-1、1,270,cm-1、1,668,cm-1和 1,298,cm-1、1,743,cm-1等位置出现了分别代表含量最高的蛋白质和脂质的特征峰；在 682,cm-1、728,cm-1等位置出现了含量为 15%,的 DNA特征峰；含量仅为 5%,的 β-胡萝卜素的特征峰在 1,154,cm-1和1,521,cm-1位置依然清晰可见.该结果说明采用GPU-rRMC方法通过增加模拟波段数量，可以更准确地获取模拟光谱数据.

图5 不同波段下均匀仿体拉曼光谱(D10)Fig.5 Raman spectrum of the homogeneous medium simulated within multi-wavelength ranges(D10)

2.3 非均匀仿体模拟实验

为进一步验证GPU-rRMC方法适用于非均匀仿体，对含的异质体仿体进行波段数为 178的数值模拟(配比浓度和模拟实验光学参数见表 1和表2).已知异质体中蛋白成分普遍高于正常组织，而DNA链由于发生大量断裂，含量相对减少，且组织中脂类和 β-胡萝卜素比正常组织含量较少，模型中异质体浓度按上述情况进行模拟配比[14].D6、D8和D10探测点的模拟结果如图6所示.

由图 6可知，在 D10探测点处获取模拟结果时，异质体几乎位于探测器正下方，此时的模拟结果在 960,cm-1出现较为明显的钙化物特征谱峰；代表蛋白质的谱线在整体光谱中明显增强，脂质、DNA和 β-胡萝卜素的特征峰明显相对减弱，这与表 1中浓度配比相吻合，说明此时的模拟结果可表现出异质体成分特性.

图6 不同探测点下的非均匀仿体拉曼光谱Fig.6 Raman spectrum of the heterogeneous medium obtained from different detectors

D8探测点的模拟结果相对于 D10处呈现整体光谱强度减弱，且此时 960,cm-1处的特征峰相对模糊；而在D6探测点处得到的光谱结果几乎与背景完全重合，说明此时异质体对模拟结果的贡献已经非常微弱.将均匀仿体和非均匀仿体在 D10探测点处的模拟光谱进行比较，结果如图7所示.

图7 均匀仿体与非均仿体模拟光谱比较(D10)Fig.7 Comparisons of the simulated Raman spectrum between the heterogeneous medium and the homogeneous medium(D10)

由图 7可知，均匀仿体的模拟光谱中以蛋白质和脂质的谱峰为主，且可看到代表 β-胡萝卜素和DNA的特征峰，未出现钙化物质特征峰；而非均匀仿体的光谱中蛋白质和DNA的特征谱峰为主，出现钙化物质的特征峰，代表脂质和 β-胡萝卜素的谱峰相对减弱.此模拟结果进一步说明了 GPU-rRMC方法对模拟非均匀仿体的拉曼光谱的有效性.

3 结 论

本文发展了一种适用于DRT算法的拉曼快速正向建模技术.通过数值模拟实验，验证了此方法对复杂生物组织模型的有效性.此方法具有如下优点：

(1)应用格林互易原理，简化了传统 RMC模型，有效地提高了计算效率.

(2)通过 GPU并行化处理，进一步提升运算速度，为获取更为精确的多波段光谱数据提供可能.

(3)适用于高密度剖分的全三维非均匀组织体，可获取任意源探分布下的光谱数据，适用于DRT计算的复杂重建算法.

(4)此方法基于MATLAB独立开发，可为DRT算法的后续开发提供开源的技术平台.

［1］Demers J L，Esmondewhite F W，Esmondewhite K A，et al. Next-generation Raman tomography instrument for non-invasive in vivo bone imaging[J].Biomedical Optics Express，2015，6(3)：793-806.

［2］Hoshi Y，Yamada Y. Overview of diffuse optical tomography and its clinical applications[J]. Journal of Biomedical Optics，2016，21：1-9.

［3］Zhao J，Zeng H，Kalia S，et al. Wavenumber selection based analysis in Raman spectroscopy improves skin cancer diagnostic specificity[J].Analyst，2016，141(3)：1034-1043.

［4］Gao F. Time-domain fluorescence molecular tomography：A FEM-diffusion-based forward model[C]//Biomedical Topical meeting. Optical Society of America，USA，2006.

［5］Uma Maheswari K，Sathiyamoorthy S，Hilda Amala J.Forward model analysis in diffuse optical tomography using RTE and FEM[J].Advances in Natural and Applied Sciences，2016，10(8)：65-74.

［6］Wang S，Zhao J，Liu H，et al. Monte Carlo simulation of in vivo Raman spectral measurements of human skin with a multi-layered tissue optical model[J].Journal of Biophotonics，2015，7(9)：703-712.

［7］Wilson R H，Dooley K A，Morris M D，et al. Monte Carlo modeling of photon transport in buried bone tissue layer for quantitive Raman spectroscopy[J].Pro SPIE，2009，7166(1)：1-10.

［8］Keller M D，Wilson R H，Mycek M A，et al. Monte Carlo model of spatially offset Raman spectroscopy for breast tumor margin analysis[J].Applied Spectroscopy，2010，64(6)：607.

［9］Hokr B H，Yakovlev V V，Scully M O. Efficient timedependent Monte Carlo simulations of stimulated Raman scattering in a turbid medium[J].ACS Photonics，2014，1(12)：1322-1329.

［10］Demers J L，Davis S C，Pogue B W，et al. Multichannel diffuse optical Raman tomography for bone characterization in vivo：A phantom study[J]. Biomedical Optics Express，2012，3(9)：2299.

［11］Dehghani H，Eames M E，Yalavarthy P K，et al. Near infrared optical tomography using NIRFAST：Algorithm for numerical model and image reconstruction[J].Communications in Numerical Methods in Engineering，2008，25(6)：711-732.

［12］Zhang Y F，Gicquel O，Taine J. Optimized emissionbased reciprocity Monte Carlo method to speed up computation in complex systems[J].International Journal of Heat & Mass Transfer，2012，55(25/26)：8172-8177.

［13］Yao R，Intes X，Fang Q. Generalized mesh-based Monte Carlo for wide-field illumination and detection via mesh retessellation[J].Biomedical Optics Express，2015，7(1)：171-184.

［14］Jermyn M，Desroches J，Aubertin K，et al. A review of Raman spectroscopy advances with an emphasis on clinical translation challenges in oncology[J].Physics in Medicine & Biology，2016，61(23)：R370-R400.

［15］于 舸，徐晓轩，吕淑华，等. 乳腺组织形态基元共焦显微拉曼光谱的研究[J]. 光谱学与光谱分析，2006，26(5)：869-873.Yu Ge，Xu Xiaoxuan，Lü Shuhua，et al. Confocal Raman microspectroscopic study of human breast morphological elements[J].Spectroscopy and Spectral Analysis，2006，26(5)：869-873(in Chinese).