高 峰 ，吳雨桐，賈夢宇，王兵元，趙會娟

(1. 天津大學(xué)精密儀器與光電子工程學(xué)院，天津 300072；2. 天津市生物醫(yī)學(xué)檢測技術(shù)與儀器重點(diǎn)實(shí)驗(yàn)室，天津 300072)

基于拉曼光譜技術(shù)的擴(kuò)散光學(xué)層析成像技術(shù)[1]旨在將光譜技術(shù)與擴(kuò)散光學(xué)層析成像[2](diffuse opti-cal tomography，DOT)技術(shù)有機(jī)融合，在構(gòu)建復(fù)雜生物組織的三維立體空間圖像的同時(shí)，利用拉曼光譜的生化指紋特征實(shí)現(xiàn)對生物組織的成分分析和識別，進(jìn)而更加準(zhǔn)確地進(jìn)行早期病變診斷[3].而建立精確高效的拉曼光子輸運(yùn)模型是實(shí)現(xiàn)該目標(biāo)的關(guān)鍵步驟.

盡管基于有限元法(finite element method，F(xiàn)EM)實(shí)現(xiàn)的擴(kuò)散近似(diffusion approximation，DA)模型在 DOT已有廣泛應(yīng)用，計(jì)算效率高，但僅適用于遠(yuǎn)源場光子輸運(yùn)的描述[4].擴(kuò)散拉曼層析成像(diffuse Raman tomography，DRT)技術(shù)所涉及的拉曼散射主要發(fā)生在淺層生物組織中，F(xiàn)EM-DA模型并不適用[5].而蒙特卡羅(Monte-Carlo，MC)方法作為輻射傳輸方程的隨機(jī)統(tǒng)計(jì)解法，能精確地描述任意區(qū)域和光學(xué)參數(shù)分布下組織中的光子傳輸行為，因而是生物組織拉曼光譜分析和成像中適用的光子輸運(yùn)模型技術(shù)[6].

傳統(tǒng)拉曼蒙特卡羅(Raman Monte-Carlo，RMC)方法計(jì)算原理簡單、模擬準(zhǔn)確，被廣泛用于拉曼光譜技術(shù)的建模過程.Wilson等[7]應(yīng)用傳統(tǒng)RMC方法獲取并分析了給定源探分布下的多波段光子計(jì)數(shù)，從而有效區(qū)分了該模型中的特定生化物質(zhì)成分.Keller等[8]應(yīng)用傳統(tǒng) RMC方法建立了基于空間分辨的多層組織模型，通過模擬分析，實(shí)現(xiàn)了對腫瘤邊界位置的識別.但傳統(tǒng)RMC模型往往需在多個(gè)散射波段下對全體積域內(nèi)的剖分節(jié)點(diǎn)逐次激發(fā)模擬，計(jì)算量大且耗時(shí)長.Hokr等[9]對傳統(tǒng) RMC方法進(jìn)行基于GPU技術(shù)的并行化加速處理，可有效縮短模擬耗時(shí)，盡管如此，該方法仍不足以有效延伸至多波段情形；Demers等運(yùn)用NIRFAST軟件建立正向模型，實(shí)現(xiàn)了全三維活體組織的DRT重建，與RMC相比，該軟件采用 FEM-DA求解，雖計(jì)算速度快，但難以精確描述 DRT所需的近源場拉曼光子傳輸過程，且該軟件高度集成封裝，擴(kuò)展其功能和算法存在著諸多局限[10-11].因此，建立一種開源且高效的拉曼光子輸運(yùn)模擬算法，對發(fā)展多波段下復(fù)雜組織體的 DRT圖像重建方法具有重要意義.

本文提出了一種基于并行化逆向追蹤技術(shù)的淺層組織拉曼蒙特卡羅模型(GPU-Reverse RMC)，以實(shí)現(xiàn) DRT算法的快速建模.根據(jù)格林函數(shù)的互易性質(zhì)，通過調(diào)換傳統(tǒng)RMC模型中給定源探對位置，在探測器位置注入大量散射光子并逆向追蹤，從而獲取與傳統(tǒng)方法等效的探測結(jié)果[12].繼而應(yīng)用GPU技術(shù)對大量光子的逆向追蹤過程進(jìn)行單指令多線程的并行化處理，可進(jìn)一步降低該模型的計(jì)算耗時(shí)[13].本文分別建立了三維均勻和非均勻數(shù)值模型，對不同 RMC方法進(jìn)行數(shù)值模擬實(shí)驗(yàn)，又應(yīng)用 GPUReverse RMC方法快速獲取了 89、178和 356個(gè)波段下的復(fù)雜組織拉曼光譜，模擬結(jié)果表明，此方法適用于建立多波段拉曼光子輸運(yùn)模型，為 DRT圖像重建工作提供一種快速有效的建模方法.

1 模型介紹

1.1 RMC傳統(tǒng)模型

傳統(tǒng)RMC模擬方法分為2步.

步驟 1 在組織體表面 rs處沿法向量方向注入波長為λx的激發(fā)光 Qx(rs,λx)，采用MC法模擬獲取激發(fā)光在組織體內(nèi)的光子密度分布Φx(rs,λx).

步驟2 將模擬譜段離散為W個(gè)波段，在給定模擬波段的光學(xué)參數(shù)下對組織內(nèi)r處的二次彌向光源進(jìn)行二次模擬，獲取表面rd處的拉曼邊界光流量，如圖1(a)所示.

圖1 兩種模型對比Fig.1 Comparisons between the two models

忽略實(shí)際測量時(shí)探頭數(shù)值孔徑及其他實(shí)驗(yàn)因素，傳統(tǒng)RMC方法原理為

式中：為在邊界探測點(diǎn)處對整個(gè)空間角 ?積分所得；c為光速；為格林函數(shù)，表示光源位于r處激發(fā)探測器于rd處可接收到的拉曼光流量；二次激發(fā)光源由給定模擬波段下組織內(nèi)r處的拉曼散射系數(shù)和激發(fā)光光子密度所決定，即

式中：為 RMC第 1步模擬所得；

通過對所有離散波段進(jìn)行RMC模擬，最終獲取組織體表面多個(gè)源探對下的拉曼模擬光譜數(shù)據(jù)

1.2 基于GPU逆向追蹤技術(shù)的RMC模型

為克服傳統(tǒng) RMC模型(圖 1(a))計(jì)算量大且耗時(shí)長的不足，本文發(fā)展了基于格林互易性質(zhì)的逆向RMC模型(reverse RMC，rRMC).已知式(1)中拉曼激發(fā)光源rQ與邊界光流量Γr滿足格林互異性 質(zhì)[12]，即

式中：Γr(r,- n)為內(nèi)部彌向光源 Qr(r,s)沿 s方向激發(fā)，在邊界 rd處沿法向量方向獲取的拉曼邊界光流量；當(dāng)調(diào)換源探位置，將光源,- n)置于 rd處沿-n方向激發(fā)，在r處可收集到沿-s方向的等效光流量數(shù)據(jù)′(r,- s).

考慮到光子的全反射效應(yīng)對于邊界光子探測效率的影響，為保證邊界入射光源與內(nèi)部彌向光源作用等效，需將光源,- n)設(shè)置為張角為ΔΩ的均勻錐束光束，ΔΩ為光子在邊界處全反射角2,π立體角積分[12]. 由于此時(shí)入射光束面積發(fā)生改變，光源,- n)的作用范圍縮小，為獲取等效的模擬結(jié)果，應(yīng)擴(kuò)大光源強(qiáng)度至相應(yīng)倍數(shù)，即

此時(shí)，將光源,- n)置于 rd處沿-n方向激發(fā)，所求等效邊界光流量可看作是對整個(gè)組織體光流量的積分，即

對注入的大量拉曼光子的輸運(yùn)過程進(jìn)行逆向追蹤并記錄，獲取等效的散射光流量.考慮到散射光子的逆向追蹤過程彼此獨(dú)立，可對其進(jìn)行基于 GPU的單指令多線程的并行化處理，從而進(jìn)一步提高模擬計(jì)算效率，該方法的流程如圖2所示.

圖2 GPU-Reverse RMC模擬流程Fig.2 Flow chart of the GPU-Reverse RMC model

2 模擬實(shí)驗(yàn)

為驗(yàn)證GPU-Reverse RMC方法的有效性，進(jìn)行如下模擬實(shí)驗(yàn).

2.1 仿體模型

建立如圖 3(a)所示 1,000,mm3的正方體數(shù)值模型，在距離上表面 1,mm的中心位置處植入體積為1,mm3的正方體異質(zhì)體.以上表面中心點(diǎn)為三維坐標(biāo)系原點(diǎn)，將光源置于原點(diǎn)位置，探測器 D1～D20沿上表面對角線均勻排布，如圖3(b)所示.

選擇 I型蛋白質(zhì)、甘油酸脂、脫氧核糖核苷酸、β-胡蘿卜素和鈣化物 5種純凈物質(zhì)，建立背景組織及含異質(zhì)體組織的生化成分模型.

拉曼散射強(qiáng)度由拉曼散射系數(shù)srμ決定，約為入射光強(qiáng)的10-9.計(jì)算該模型的拉曼散射系數(shù)為

式中：m表示組織中的不同物質(zhì)成分；(i=1,2,… ,W)表示相對拉曼散射截面，由各物質(zhì)歸一化拉曼基準(zhǔn)譜給出，表示單位濃度下發(fā)生拉曼散射的概率；CRm為相對于各物質(zhì)單位濃度下的相對濃度.背景及異質(zhì)體中生化成分的相對濃度 CRm及拉曼散射系數(shù)srμ由表1給出[15].

圖3 仿體結(jié)構(gòu)Fig.3 Phantom model

表1 背景組織和異質(zhì)體的生化成分的濃度占比Tab.1 RCs of the background and target constituents

實(shí)驗(yàn)中往往為避免熒光干擾，選擇激發(fā)光波長為785,nm，拉曼散射波段范圍為827～914,nm.考慮到此波段下吸收、散射光學(xué)參數(shù)的改變對于拉曼散射結(jié)果的影響甚微，模擬實(shí)驗(yàn)中設(shè)置散射過程為840,nm下的組織光學(xué)參數(shù)，模擬實(shí)驗(yàn)選擇近似皮膚組織的光學(xué)參數(shù)，如表2所示[8].

表2 背景及異質(zhì)體光學(xué)參數(shù)Tab.2 Optical properties of the background and the target

表2中axμ、sxμ分別代表激發(fā)光波長下吸收、散射系數(shù)；amμ、smμ分別代表散射波長下吸收、散射系數(shù)；n為光在組織體表面的折射率；g為各項(xiàng)異性系數(shù)，表示散射的定向性[9].

2.2 均勻仿體模擬實(shí)驗(yàn)

本研究計(jì)算平臺所用 GPU顯卡的型號是Nvidia GTX 980.設(shè)置模擬光源為 1,mW 的連續(xù)準(zhǔn)直光源，注入光子數(shù)為 1014,s-1.模型離散為 20×20×20個(gè)立方體體元.分別使用傳統(tǒng) RMC(Traditional RMC，tRMC)方法和逆向 RMC(Reverse RMC，rRMC)方法對均勻背景仿體進(jìn)行數(shù)值實(shí)驗(yàn)?zāi)M，得到探測器 D1～D20探測點(diǎn)處的拉曼邊界光流量，并計(jì)算其相對誤差(relative error，RE)為

由圖 4可知，當(dāng)源探距離小于 4,mm 時(shí)，逆向RMC方法模擬結(jié)果的相對誤差近似為零，說明此時(shí)的模擬結(jié)果可有效反映拉曼光的擴(kuò)散效應(yīng)；當(dāng)源探距離大于4,mm時(shí)，獲取遠(yuǎn)離光源部分的模擬數(shù)據(jù)時(shí)出現(xiàn)較為明顯的相對誤差，這是由于此時(shí)探測信號的光子量級低于 102，模擬產(chǎn)生誤差在較弱探測信號強(qiáng)度下較為明顯.

圖4 傳統(tǒng) RMC方法和逆向 RMC方法模擬結(jié)果比較(D1～D20)Fig.4 Comparisons of simulating result by the tRMC method and the rRMC method(D1～D20)

逆向 RMC方法的模擬信號強(qiáng)度比傳統(tǒng)方法略有偏低，這是由于該方法僅通過在邊界探測點(diǎn)處對光源進(jìn)行激發(fā)模擬并統(tǒng)計(jì)求和，與傳統(tǒng)的逐點(diǎn)模擬方法相比，在遠(yuǎn)離表面光源的部分存在著一定程度的能量損失，因而統(tǒng)計(jì)求和后其結(jié)果往往略有偏低.

考慮到只有基于足夠大量隨機(jī)實(shí)驗(yàn)統(tǒng)計(jì)求解時(shí)，模擬結(jié)果才能有效反映真實(shí)情況，因此當(dāng)探測信號的光子量級低于102時(shí)，其結(jié)果易受噪聲等因素干擾，難以有效還原光子在組織中的擴(kuò)散情況，參考意義較??；針對此情況，可考慮通過延長模擬積分時(shí)間以提高模擬信號的強(qiáng)度，從而保證信號的準(zhǔn)確性.

對第2.1節(jié)中的數(shù)值模型分別采用tRMC方法、GPU-tRMC方法和 GPU-rRMC方法進(jìn)行計(jì)算機(jī)數(shù)值模擬，注入激發(fā)光光子數(shù)為 107.將單一波段、89、178、356個(gè)波段的計(jì)算耗時(shí)估算和計(jì)算機(jī)模擬結(jié)果列于表3.

表3 3種RMC方法計(jì)算時(shí)間對比(20個(gè)探測點(diǎn))Tab.3 Comparisons of the execution time among the three methods(20 detectors)

當(dāng)獲取單一波段下的邊界光流量時(shí)，tRMC方法和GPU-tRMC方法需對組織內(nèi)全部離散體元進(jìn)行逐點(diǎn) MC模擬，共剖分 8,000個(gè)體元，因此需要進(jìn)行8,001次 MC模擬(含激發(fā)光模擬)，前者估算耗時(shí)2.7×107,min，約 450,000,h，顯然無法實(shí)現(xiàn)，后者可通過GPU加速將計(jì)算時(shí)間有效縮短至2.5,h；本文提出的GPU-rRMC方法僅需在探測位置處共進(jìn)行MC模擬，共計(jì) 21次(含激發(fā)光模擬)，計(jì)算時(shí)間可有效低至0.42,min.

當(dāng)獲取多波段下的模擬光譜數(shù)據(jù)時(shí)，tRMC方法由于計(jì)算量過大，只能估算其模擬耗時(shí)，89個(gè)波段下約耗時(shí) 2.4×109,min，顯然無法實(shí)際應(yīng)用；即使將tRMC方法經(jīng)GPU加速，模擬89個(gè)波段下的光譜數(shù)據(jù)仍需耗時(shí) 1.4×104,min；采用 GPU-rRMC方法獲取波段數(shù)為 89、178、356時(shí) 3種波段下的光譜數(shù)據(jù)分別耗時(shí) 50,min、85,min和 180,min，此時(shí)計(jì)算時(shí)長合理，適于實(shí)際應(yīng)用.

分析上述結(jié)果的主要原因?yàn)椋孩賢RMC方法需對組織內(nèi)各離散節(jié)點(diǎn)逐點(diǎn)進(jìn)行模擬，計(jì)算次數(shù)遠(yuǎn)遠(yuǎn)大于 rRMC方法，因此造成該方法下模擬耗時(shí)嚴(yán)重；②GPU-tRMC方法可有效提高 tRMC計(jì)算速度約105倍，但采用該方法依然需進(jìn)行大量重復(fù)模擬過程，因而難以有效延伸至多波段的情形；③GPU-rRMC方法從模擬算法和實(shí)現(xiàn)手段兩方面提高計(jì)算效率，并在模擬誤差可接受范圍內(nèi)獲取有效的模擬數(shù)據(jù).因此該方式是獲取組織拉曼光譜適用的方法.

采用GPU-rRMC方法獲取3種波段下的模擬光譜，對比模擬結(jié)果以顯示不同模擬波段對結(jié)果準(zhǔn)確度的影響.獲取 D10探測點(diǎn)的模擬結(jié)果，如圖 5所示.

根據(jù)圖 5可知，波段數(shù)為 89的模擬光譜在920,cm-1、1,003,cm-1、1,270～1,315,cm-1波段處出現(xiàn)部分光譜信息丟失現(xiàn)象，在 1,450,cm-1、1,521,cm-1、1,668,cm-1處的譜峰強(qiáng)度略低于基準(zhǔn)譜；模擬波段數(shù)增加至178時(shí)，大部分波段下的譜峰強(qiáng)度得以還原，但部分特征譜峰仍未清晰反映出來；波段數(shù)為 356的光譜結(jié)果中可以清晰看到其光譜結(jié)果與基準(zhǔn)譜較為 吻 合 ，在 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等位置出現(xiàn)了分別代表含量最高的蛋白質(zhì)和脂質(zhì)的特征峰；在 682,cm-1、728,cm-1等位置出現(xiàn)了含量為 15%,的 DNA特征峰；含量僅為 5%,的 β-胡蘿卜素的特征峰在 1,154,cm-1和1,521,cm-1位置依然清晰可見.該結(jié)果說明采用GPU-rRMC方法通過增加模擬波段數(shù)量，可以更準(zhǔn)確地獲取模擬光譜數(shù)據(jù).

圖5 不同波段下均勻仿體拉曼光譜(D10)Fig.5 Raman spectrum of the homogeneous medium simulated within multi-wavelength ranges(D10)

2.3 非均勻仿體模擬實(shí)驗(yàn)

為進(jìn)一步驗(yàn)證GPU-rRMC方法適用于非均勻仿體，對含的異質(zhì)體仿體進(jìn)行波段數(shù)為 178的數(shù)值模擬(配比濃度和模擬實(shí)驗(yàn)光學(xué)參數(shù)見表 1和表2).已知異質(zhì)體中蛋白成分普遍高于正常組織，而DNA鏈由于發(fā)生大量斷裂，含量相對減少，且組織中脂類和 β-胡蘿卜素比正常組織含量較少，模型中異質(zhì)體濃度按上述情況進(jìn)行模擬配比[14].D6、D8和D10探測點(diǎn)的模擬結(jié)果如圖6所示.

由圖 6可知，在 D10探測點(diǎn)處獲取模擬結(jié)果時(shí)，異質(zhì)體幾乎位于探測器正下方，此時(shí)的模擬結(jié)果在 960,cm-1出現(xiàn)較為明顯的鈣化物特征譜峰；代表蛋白質(zhì)的譜線在整體光譜中明顯增強(qiáng)，脂質(zhì)、DNA和 β-胡蘿卜素的特征峰明顯相對減弱，這與表 1中濃度配比相吻合，說明此時(shí)的模擬結(jié)果可表現(xiàn)出異質(zhì)體成分特性.

圖6 不同探測點(diǎn)下的非均勻仿體拉曼光譜Fig.6 Raman spectrum of the heterogeneous medium obtained from different detectors

D8探測點(diǎn)的模擬結(jié)果相對于 D10處呈現(xiàn)整體光譜強(qiáng)度減弱，且此時(shí) 960,cm-1處的特征峰相對模糊；而在D6探測點(diǎn)處得到的光譜結(jié)果幾乎與背景完全重合，說明此時(shí)異質(zhì)體對模擬結(jié)果的貢獻(xiàn)已經(jīng)非常微弱.將均勻仿體和非均勻仿體在 D10探測點(diǎn)處的模擬光譜進(jìn)行比較，結(jié)果如圖7所示.

圖7 均勻仿體與非均仿體模擬光譜比較(D10)Fig.7 Comparisons of the simulated Raman spectrum between the heterogeneous medium and the homogeneous medium(D10)

由圖 7可知，均勻仿體的模擬光譜中以蛋白質(zhì)和脂質(zhì)的譜峰為主，且可看到代表 β-胡蘿卜素和DNA的特征峰，未出現(xiàn)鈣化物質(zhì)特征峰；而非均勻仿體的光譜中蛋白質(zhì)和DNA的特征譜峰為主，出現(xiàn)鈣化物質(zhì)的特征峰，代表脂質(zhì)和 β-胡蘿卜素的譜峰相對減弱.此模擬結(jié)果進(jìn)一步說明了 GPU-rRMC方法對模擬非均勻仿體的拉曼光譜的有效性.

3 結(jié) 論

本文發(fā)展了一種適用于DRT算法的拉曼快速正向建模技術(shù).通過數(shù)值模擬實(shí)驗(yàn)，驗(yàn)證了此方法對復(fù)雜生物組織模型的有效性.此方法具有如下優(yōu)點(diǎn)：

(1)應(yīng)用格林互易原理，簡化了傳統(tǒng) RMC模型，有效地提高了計(jì)算效率.

(2)通過 GPU并行化處理，進(jìn)一步提升運(yùn)算速度，為獲取更為精確的多波段光譜數(shù)據(jù)提供可能.

(3)適用于高密度剖分的全三維非均勻組織體，可獲取任意源探分布下的光譜數(shù)據(jù)，適用于DRT計(jì)算的復(fù)雜重建算法.

(4)此方法基于MATLAB獨(dú)立開發(fā)，可為DRT算法的后續(xù)開發(fā)提供開源的技術(shù)平臺.

［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，F(xiàn)ang 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］于 舸，徐曉軒，呂淑華，等. 乳腺組織形態(tài)基元共焦顯微拉曼光譜的研究[J]. 光譜學(xué)與光譜分析，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).