基于蒙特卡羅的腦血氧檢測光源—探測器最優(yōu)距離研究

曹 焱，邢志明，趙 斌，董祥美，史以玨，高秀敏

（1.上海理工大學(xué) 光電信息與計算機(jī)工程學(xué)院，上海 200093；2.上海交通大學(xué)醫(yī)學(xué)院附屬瑞金醫(yī)院，上海 200025）

0 引言

大腦是人體的一個重要器官，像一個中樞系統(tǒng)一樣指揮著人們的日常行為與活動。腦部的血流量約占心臟血液輸出量的15%，耗氧量約占全身耗氧量的20%［1］，所以大腦對于缺血缺氧是十分敏感的［2］。臨床上對腦血氧的監(jiān)測可以幫助醫(yī)生在手術(shù)中實(shí)時監(jiān)測患者健康狀況［3-4］，以及進(jìn)行一些腦部疾病的診斷［5-6］。因此，用于臨床監(jiān)測和診斷的血氧監(jiān)測儀器可為醫(yī)生提供很大的幫助。

利用功能性近紅外光譜技術(shù)［7］可以很容易地檢測到腦功能的血流動力學(xué)信息［8］。在實(shí)際的腦血氧測量中，被檢測者往往會因人體的生理活動（如心動周期、呼吸、血壓變化等）引入干擾，從而給腦血氧信號檢測造成了一定影響，這些生理干擾可以采用自適應(yīng)濾波法［9-10］去除。自適應(yīng)濾波法對于抑制特定的生理干擾有比較好的效果，但是需要借助額外的設(shè)備獲取參考信號。通常使用單個探測器進(jìn)行腦血氧信號檢測時會面臨以上問題，因?yàn)閮H有一個探測器無法獲取參考信號。為解決單個探測器檢測所面臨的問題，有學(xué)者提出了多距離測量方法，即使用雙探測器測量腦血氧信號［11-12］。一個探測器檢測深層組織中帶有生理活動干擾的血氧信號，另一個探測器檢測淺層信號作為參考信號，對兩個探測器檢測的信號作差分處理［13］，獲得純度更高的局部信號用于后續(xù)計算處理。此外，臨床和科研中光源—探測器（Source-Detector，SD）距離的選擇大多根據(jù)經(jīng)驗(yàn)公式［14］，即探測器能檢測到的最大組織深度是光源與探測器距離的一半［15］，或者參考市面上現(xiàn)有的腦血氧檢測產(chǎn)品來選擇光源—探測器之間的距離［16］。這些方法選擇的光源與探測器距離是固定的，但由于顱腦各部分組織結(jié)構(gòu)針對不同的年齡、性別、個體等會存在差異，如正常人和頭部有損傷的病人（如硬膜腦血腫［17］）等，經(jīng)驗(yàn)公式可能無法很好地適用于所有人群的腦血氧檢測，無法滿足光源—探測器之間的最佳距離選擇需求。

為解決上述腦血氧檢測中SD 距離選取時遇到的問題，本文根據(jù)多距離測量方法要求，給光源—探測器選擇一組恰當(dāng)?shù)木嚯x。光源放置于固定位置時，探測器需一個放置于近處，用于檢測頭皮和顱骨處的血液信號作為參考信號；另一個放置于遠(yuǎn)處，用于檢測灰質(zhì)處受到生理活動干擾的血液信號。再結(jié)合蒙特卡羅方法，從光在組織中的傳輸特性著手，同時考慮到顱腦組織的結(jié)構(gòu)差異，選定一組顱腦組織的特性參數(shù)，建立5 層顱腦組織平板模型進(jìn)行仿真。通過統(tǒng)計大量光子在組織中的傳輸情況，并對仿真得到的數(shù)據(jù)進(jìn)行分析與處理，可以很好地克服在顱腦組織存在差異時應(yīng)用經(jīng)驗(yàn)公式的缺陷，從而選擇出符合該組顱腦組織參數(shù)的最優(yōu)的光源—探測器距離。

1 蒙特卡羅模型與方法

1.1 蒙特卡羅原理

生物組織在光學(xué)上可近似看作混沌介質(zhì)，當(dāng)光照射在生物組織上時，會發(fā)生反射、吸收和散射等現(xiàn)象，光在生物組織中的宏觀傳播行為已不能用基本的朗伯—比爾定律來描述。近些年得益于計算機(jī)技術(shù)的快速發(fā)展，蒙特卡羅（Monte Carlo）方法在各個物理仿真實(shí)驗(yàn)中使用得越來越多，尤其在生物醫(yī)學(xué)、光學(xué)等領(lǐng)域［18］。Monte Carlo 模擬［19-20］是一種基于概率統(tǒng)計的數(shù)學(xué)物理模擬方法，通過相應(yīng)的概率模型或隨機(jī)過程計算，模擬粒子在組織中傳輸時受到的影響［21］。通過對大量光子在組織中傳播路徑的仿真，追蹤光子運(yùn)動的軌跡點(diǎn)，可得到光在組織內(nèi)部和表面的統(tǒng)計分布規(guī)律。

1.2 光在組織中傳輸?shù)拿商乜_模擬方法

本文利用多層蒙特卡羅（Monte Carlo Multi-Layer，MCML）程序進(jìn)行仿真，MCML 程序是Wang 等［19］用標(biāo)準(zhǔn)C 語言開發(fā)的程序包，用于模擬分層組織的光子傳輸［22］。蒙特卡羅程序的基本流程如下［23-24］：

（1）確定起始跟蹤點(diǎn)和初始化光子（權(quán)重、方向、位置）后，在組織表面發(fā)射光子，并為其分配隨機(jī)步長。

（2）移動步長s 后，更新光子位置，確定光子下一次的碰撞位置。

（3）光子移動步長后會隨機(jī)進(jìn)行吸收和散射，散射作用會改變光子的前進(jìn)方向，吸收作用會降低光子權(quán)值。根據(jù)當(dāng)前位置計算由于散射造成的方向偏轉(zhuǎn)，更新傳播方向余弦。計算光子在當(dāng)前位置由于吸收造成的權(quán)重衰減，更新權(quán)重值，并對光子權(quán)重是否低于閾值進(jìn)行判斷。如果低于閾值，或逸出組織上下表面時，采用輪盤賭來決定光子存活。判定為死亡的，結(jié)束當(dāng)前光子的追蹤。

（4）對光子是否碰撞邊界進(jìn)行判斷，根據(jù)判斷結(jié)果進(jìn)行相應(yīng)的處理。

（5）對于未逃逸出組織的光子繼續(xù)進(jìn)行步驟（2）—（4），直到程序結(jié)束。

蒙特卡羅模擬流程如圖1所示。

Fig.1 Flow of Monte Carlo simulation圖1 蒙特卡羅模擬流程

1.3 光在腦組織中傳輸路徑的仿真研究

研究光在顱腦組織中的傳輸路徑，需對每個光子的傳輸軌跡進(jìn)行追蹤，但Wang 原有的程序并不能實(shí)現(xiàn)此功能。因此，需要在Wang 的仿真程序中進(jìn)行優(yōu)化，以獲得光子在組織中的傳輸路徑。由于光子在組織內(nèi)傳輸是隨機(jī)、無序的，所以采用數(shù)據(jù)鏈表來保存光子的傳輸路徑。光子在組織中傳輸時的平均自由路程很短，在運(yùn)算時保存的信息會占用很大的計算機(jī)運(yùn)行內(nèi)存，故用d來判斷光子運(yùn)動距離，優(yōu)化保存的信息。

式中，d表示光子傳輸時相鄰兩個位置的距離，（xk，yk，zk）表示光子當(dāng)前的位置，更新后的位置坐標(biāo)為（x，y，z）。只有當(dāng)d值大于某一閾值時，才將當(dāng)前位置點(diǎn)的坐標(biāo)記錄到鏈表中［25］。蒙特卡羅仿真時需要依靠大量光子才能得到較為準(zhǔn)確的結(jié)果，在實(shí)際仿真時，考慮到計算機(jī)的內(nèi)存空間使用情況，故只關(guān)注探測器處接收到的光子，對這些被接收到的光子路徑進(jìn)行保存，如圖2所示。

Fig.2 Schematic diagram of light source and detector placement圖2 光源與探測器位置示意圖

光源和探測器放置的位置如圖2（a）所示，X 軸為光源和探測器所在的水平方向，Y 軸為頭皮表層的豎直方向，Z軸為腦組織深度方向。圖2（b）中光子入射點(diǎn)（0，0，0）與探測器的距離為SD，探測器與X 軸的夾角為θ。一定數(shù)量的光子從原點(diǎn)入射，在組織內(nèi)進(jìn)行隨機(jī)傳輸運(yùn)動，部分光子可到達(dá)接收面并被接收。為便于后續(xù)處理，將探測器接收面看成邊長為Δs 的正方形。只有成功到達(dá)接收面的光子，才記錄其在組織內(nèi)的傳輸路徑，其他光子則不予記錄。

仿真程序中所使用的五層顱腦模型光學(xué)參數(shù)如表1所示［26-28］。表中各參數(shù)的定義如下：n是折射率，g是各向異性因子，μa/cm-1是吸收系數(shù)，μs/cm-1是散射系數(shù)，t/cm 是各顱腦層的組織厚度。

Table 1 Optical parameters of the five-layer cranial model表1 五層顱腦模型光學(xué)參數(shù)

2 蒙特卡羅仿真與討論

為便于處理，本文將探測器與光子入射點(diǎn)設(shè)置為同一水平線，即設(shè)置圖2（b）中的θ值為0°，接收面的邊長Δs 為1mm，從原點(diǎn)（0，0，0）入射的光子數(shù)為107個，d的閾值為1mm。設(shè)置完成后，選擇不同的光源—探測器（SD）距離值運(yùn)行MCML 程序。值得一提的是，由于此仿真需要對大量數(shù)據(jù)進(jìn)行統(tǒng)計分析，故對電腦內(nèi)存要求較高。仿真硬件條件為：Intel Core i9-12900處理器，128GB 的RAM 內(nèi)存。

2.1 光子傳輸路徑

按如上條件設(shè)置后，可得到不同SD 距離下的光子傳輸路徑數(shù)據(jù)，選取SD=1cm 和3.5cm 的數(shù)據(jù)進(jìn)行處理，繪制光子在組織內(nèi)的傳輸路徑圖，如圖3所示。

Fig.3 Photon transport paths in tissue圖3 光子在組織內(nèi)傳輸路徑

圖3（a）、（b）分別是SD=1cm 和SD=3.5cm 的傳輸路徑三維圖，中心區(qū)域是一團(tuán)聚集在一起的密集路徑，邊緣區(qū)域的路徑逐漸稀疏。由圖3（a）和（b）可明顯看到，SD=3.5cm 時的傳輸路徑核心區(qū)域比SD=1cm 時大一些。

傳輸路徑圖反映了光子在顱腦組織內(nèi)傳輸?shù)牟糠智闆r，而確定光源—探測器（SD）間的最優(yōu)距離，還需要借助其他參數(shù)進(jìn)行判斷。本文通過研究光子傳輸時的最大可達(dá)深度和散射距離，作為選擇光源—探測器最優(yōu)距離的參考。

2.2 光子最大可達(dá)深度

探測深度是腦血氧檢測中的一個重要參數(shù)，光子傳輸路徑的最大可達(dá)深度一定程度上反映了光子在顱腦內(nèi)可探測的深度。根據(jù)每條光子傳輸路徑上Z 軸坐標(biāo)的極大值，繪制出如圖4 所示的不同SD 距離下的光子最大可達(dá)深度概率分布圖。圖中橫坐標(biāo)是光子傳輸?shù)纳疃戎担v坐標(biāo)是光子最大可達(dá)深度的概率值。可以看出，不同SD 距離對應(yīng)的概率分布大致是正態(tài)分布。圖4（a）稍有不同，左側(cè)區(qū)域的概率值比右側(cè)區(qū)域密集些，主要是由于SD 的距離過小，光子多在淺層區(qū)域傳輸，傳輸?shù)缴顚咏M織的光子相對較少。這也說明當(dāng)SD<2cm 時，只有相對較少的光子可到達(dá)深層組織，這對于選擇近處探測器與光源之間的距離有很重要的參考價值。

Fig.4 Probability distribution of maximum reachable depth for 300 photonic transmission paths at different SD distances圖4 不同SD距離下300條光子傳輸路徑的最大可達(dá)深度概率分布

為便于分析，計算不同的SD 距離下光子最大可達(dá)深度在不同組織層的占比情況。如圖5 所示，橫坐標(biāo)是不同的SD 距離，縱坐標(biāo)是各顱腦組織在不同SD 距離時的最大深度占比情況。可看到隨著SD 距離的增大，頭皮層和顱骨層的占比不斷降低；腦脊液層的占比非常小，可忽略不計；灰質(zhì)層的占比基本大于15%。頭皮和顱骨層占比在SD=1cm 處取得最大值23.78%。SD=2cm～4.5cm 時，光子最大可達(dá)深度在灰質(zhì)層的占比于SD=3.5cm 處取得最大值19.59%

2.3 光子傳輸路徑在不同層的散射距離

僅采用最大可達(dá)深度描述光子對深層顱腦組織探測的情況還不夠客觀、全面，最大可達(dá)深度只能了解光子是否到達(dá)某一深度值，其他信息不得而知。

光子在顱腦組織中傳輸時，散射距離也是一個重要參數(shù)。因此，計算不同SD 距離下光子傳輸路徑在各層（取每一層深度為0.1cm）散射距離占路徑總長的比例，并繪制散射距離占比變化趨勢圖，如圖6所示。

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

圖6 橫坐標(biāo)表示顱腦組織的深度，縱坐標(biāo)表示光子在各層的散射距離占傳輸路徑總長的比值。隨著深度增加，光子在深層組織處的散射距離所占比重不斷減少，表明光子在深層組織傳輸時的散射距離短。當(dāng)深度達(dá)到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 中，橫坐標(biāo)是不同的SD 距離，縱坐標(biāo)是光子分別在頭皮和顱骨、腦脊液、灰質(zhì)層的散射距離占總散射距離的比重?？擅黠@看到，光子的散射距離在頭皮和顱骨層所占比重很大，基本都在25%左右。由于腦脊液層的特殊光學(xué)性質(zhì)（低吸收、低散射），故光子在腦脊液層的散射距離比重都很小。光子在灰質(zhì)層的散射距離比重不等，與SD 距離相關(guān)，當(dāng)SD=2～4.5cm 時，灰質(zhì)層的散射距離比重在SD=4.5cm 處取最大值10.22%。這是因?yàn)殡S著SD 距離的不斷增大，與較小的SD 相比，有更多光子可以到達(dá)更深層的組織進(jìn)行傳輸。故隨著SD 的增大，灰質(zhì)層的散射距離比重也在增大。

3 結(jié)果分析

根據(jù)多距離測量方法對SD 距離的選擇要求，近距離探測器D1 檢測頭皮和顱骨層的血液信息，遠(yuǎn)距離探測器D2 檢測腦灰質(zhì)層的血液信息。在對距離進(jìn)行選擇時，需要同時考慮3 個方面：一是光源在傳輸時可達(dá)的最大深度；二是光源在組織中傳輸時攜帶信息多少；三是探測器實(shí)際接收到的信號強(qiáng)度。

有效深度比（Effective Depth Ratio，EDR）是探測器接收到來自有效深度層的光子與所有組織層光子的比值。將頭皮和顱骨層看作第一層，將腦脊液、灰質(zhì)和白質(zhì)層看作第二層，則在式（2）中，第一層（最大深度）為最大可達(dá)深度在頭皮和顱骨層，且并未穿過顱骨層的光子，第二層（最大深度）為穿透顱骨層的光子。

根據(jù)上文所述，計算出不同SD 距離下的EDR 值并繪制圖像，如圖8（a）所示。圖8（a）中，橫坐標(biāo)為不同的SD 距離，縱坐標(biāo)為EDR 值。從圖8（a）中可以看出，隨著SD 距離的增大，有效深度比EDR 逐漸減小，說明隨著SD 的增大，傳輸?shù)缴顚咏M織的光子增多。EDR 值越大，說明SD 的距離越小，光子傳輸時位于頭皮和顱骨層的光子相對較多；反之，SD 距離越大，位于灰質(zhì)層的光子相對較多。故可將EDR 值作為選擇光源與探測器間距離的重要依據(jù)。

Fig.8 Effective depth ratio and grey matter layer information ratio at different SD distances圖8 不同SD距離下的有效深度比與灰質(zhì)層信息比

此外，光子在組織中傳輸時，攜帶的信息量大小與光子散射距離有很大關(guān)系。光子傳輸?shù)腅DR 值并不能完全反映探測器接收到來自灰質(zhì)層的信息量大小，故對圖7 中的結(jié)果作進(jìn)一步處理與分析。計算在不同SD 距離下光子傳輸時的灰質(zhì)層信息比（Gray Matter Layer Information Ratio，GMLIR），表示為：

式中，灰質(zhì)層（散射距離）為光子在灰質(zhì)層的散射距離與總散射距離的比；頭皮、顱骨層（散射距離）為光子在頭皮和顱腦層的散射距離與總散射距離的比。繪制如圖8（b）所示的圖像，從圖中可明顯看出，隨著SD 距離的增大，灰質(zhì)層信息不斷增多。在腦血氧檢測中，灰質(zhì)層信息比的增大意味著探測器探測到的有效信息增多，這對于實(shí)際檢測是很重要的。故GMLIR 值也作為選擇光源與探測器距離的一個重要依據(jù)。

在實(shí)際選擇光源與探測器距離時，還需要考慮探測器接收到的信號強(qiáng)度。本文在研究中將光電探測器貼近頭皮表層進(jìn)行必要的遮光處理后，分別使用760nm 和850nm的激光二極管（LD）在不同的SD 距離下照射頭皮表層，記錄探測器的電壓值。如圖9 所示，當(dāng)不使用LD 照射時，探測器的電壓均值為3.67mV（參考值），而當(dāng)使用LD 時，LD和探測器間的距離與信號強(qiáng)度成反比。SD 距離較小時，如SD=1cm，探測器接收到760nm 和850nm 信號的電壓均值大于400mV。隨著SD 距離的增大，信號強(qiáng)度開始明顯變小，說明SD 距離對于信號強(qiáng)度的影響非常大。探測器與LD 相距4cm 的位置處，探測器的電壓均值非常接近參考值，意味著探測器接收到的信號已很微弱。在確定SD距離時，需要考慮到探測器接收信號強(qiáng)弱這一方面的影響，故在選擇SD 時，將不再考慮SD>4cm 的情況。

Fig.9 Signal intensity received by the detector at different wavelengths with different SD distances圖9 不同SD距離時不同波長下探測器接收到的信號強(qiáng)度

根據(jù)以上分析，EDR 與GMLIR 是選擇SD 距離的重要依據(jù)。近距離探測器D1 主要檢測頭皮和顱骨層信息，需要EDR 值較大，且不檢測灰質(zhì)層信息，因此需要GMLIR 值較小；遠(yuǎn)距離探測器D2 主要檢測灰質(zhì)層信息，其中含有頭皮和顱骨層的干擾信息少，故需要EDR 值較小，GMLIR 值較大。在選取光源與D2 的距離時，還需考慮探測器實(shí)際接收信號的強(qiáng)度。綜上所述，本文選擇的一組最優(yōu)光源與探測器距離為SD1=1cm，SD2=3.5cm。

4 結(jié)語

本文通過選取一組顱腦組織的特性參數(shù)，利用蒙特卡羅方法對顱腦組織模型進(jìn)行仿真研究，獲得光子在頭部組織的傳輸路徑，進(jìn)而分析光子傳輸時的最大可達(dá)深度和散射距離，以及有效深度比EDR 和灰質(zhì)層信息比GMLIR。根據(jù)多距離測量方法的要求，考慮了EDR、GMLIR 和探測器實(shí)際接收到的信號強(qiáng)度，確定一組適合選定顱腦組織參數(shù)的最優(yōu)光源—探測器距離，SD1=1cm 和SD2=3.5cm。該組距離不僅可以最大程度地檢測灰質(zhì)層的血氧信息，而且可以最大程度地處理來自頭皮和顱骨層的干擾，從而獲得純度更高的局部腦血氧信號。利用多距離測量方法設(shè)計的腦血氧檢測裝置，可以選擇此組最優(yōu)距離進(jìn)行傳感器探頭設(shè)計。本研究結(jié)果對腦血氧檢測裝置的設(shè)計，以及利用近紅外光譜技術(shù)檢測腦血氧的精度優(yōu)化具有一定指導(dǎo)意義。