基于蒙特卡羅法硬膜血腫厚度計(jì)算的數(shù)值仿真研究

范鑫燕，羅海軍，2，李妍妍，向 洋，羅 霞，覃 睿，郭 盼

（1.重慶師范大學(xué)光電功能材料重慶市重點(diǎn)實(shí)驗(yàn)室,重慶 401331；2.重慶國(guó)家應(yīng)用數(shù)學(xué)中心,重慶 401331）

1 引言

硬膜出血是創(chuàng)傷性腦損傷中最常見(jiàn)也是最嚴(yán)重的二次損傷病變，致死率30%～40%，致殘率接近100%［1］，其中病變位置多發(fā)生于顱骨下方腦脊液中，出血量可在幾小時(shí)內(nèi)堆積形成血腫，血腫擴(kuò)大可能引起腦中風(fēng)、神經(jīng)病理惡化［2］等后遺癥，導(dǎo)致永久性精神疾病甚至危及生命.目前醫(yī)院評(píng)估硬膜血腫是采用計(jì)算機(jī)斷層掃描（Computed Tomography，CT）和磁共振成像（Magnetic Resonance Imaging，MRI）等檢測(cè)精度高的大型精密設(shè)備，但CT 與MRI 價(jià)格高、體積龐大，不適用實(shí)時(shí)連續(xù)監(jiān)測(cè)，且往往伴隨著一定的電磁輻射，不適用重度患者、嬰幼兒和孕婦等［3］.

近紅外光譜技術(shù)（Near Infra-Red Spectroscopy，NIRS）作為監(jiān)測(cè)神經(jīng)創(chuàng)傷主流的非侵入式輔助手段，是基于生物組織結(jié)構(gòu)的光譜分析，其典型的應(yīng)用是檢測(cè)顱內(nèi)占位血腫［4］.光譜600 nm～950 nm 波段是人體組織的“光學(xué)窗口”，該波段光可穿過(guò)腦部組織幾厘米深，皮膚、骨骼等組織作為均勻介質(zhì)呈現(xiàn)出散射系數(shù)大，吸收系數(shù)小的特性，血液對(duì)于該波段內(nèi)的近紅外光敏感［5］，血腫含有豐富的血紅蛋白，呈現(xiàn)高吸收的特性，對(duì)光能量吸收強(qiáng)度是周圍正常腦組織的十倍以上，分析光的衰減信息可得到血腫層信息［6］.20 世紀(jì)末美國(guó)賓夕法尼亞大學(xué)的Chance B［7］等將NIRS 用于腦血腫檢測(cè)的研究，Salonia 等人對(duì)兒童腦出血檢測(cè)進(jìn)行了研究［8］，國(guó)內(nèi)軍事科學(xué)院張彥軍利用篩選器建立血腫吸收系數(shù)預(yù)測(cè)模型對(duì)血腫厚度進(jìn)行判斷［9］，天津工業(yè)大學(xué)王金海團(tuán)隊(duì)采用多通道差分吸光度法檢測(cè)腦異物［10］，各研究表明NIRS 可用于顱腦血腫檢測(cè).為實(shí)現(xiàn)硬膜血腫的快速檢測(cè)與厚度評(píng)估，本研究基于NIRS 與蒙特卡羅模擬（Monte Carlo，MC），結(jié)合數(shù)值分析提出一種預(yù)測(cè)硬膜血腫厚度的數(shù)學(xué)方法，即在正向測(cè)量光子信息的基礎(chǔ)上，構(gòu)建不同檢測(cè)距離下的函數(shù)模型，建立逆向函數(shù)矩陣，該方法的優(yōu)點(diǎn)是充分考慮到了檢測(cè)器的分布和硬膜血腫對(duì)光傳導(dǎo)的影響，從定量的角度分析血腫對(duì)光衰減的影響，能夠獲得一種通過(guò)顱腦表面光子信息求解硬膜血腫厚度的方法.

2 理論原理與模型建立

2.1 蒙特卡羅理論

顱腦結(jié)構(gòu)中腦組織（灰質(zhì)與白質(zhì)）為不規(guī)則形狀的起伏層，目前沒(méi)有統(tǒng)一認(rèn)可的解析式方案，復(fù)雜組織的靈敏度和穿透深度依賴于數(shù)值解決辦法［11］.MC 是一種描述光在生物組織傳播的經(jīng)典數(shù)值模擬技術(shù)，其優(yōu)勢(shì)在于不需要設(shè)定光子傳播的解析形式、復(fù)雜的邊界條件，因此不會(huì)被擴(kuò)散和其他近似方程所偏置，被認(rèn)為是光在組織傳播的“金標(biāo)準(zhǔn)”［12］.MC模擬傳播示意如圖1所示.

圖1 蒙特卡羅模擬傳播示意圖

光子以權(quán)重進(jìn)W=1 入組織，產(chǎn)生隨機(jī)步長(zhǎng)step，每產(chǎn)生一個(gè)新步長(zhǎng)都會(huì)發(fā)生吸收，散射.在傳播過(guò)程中，吸收導(dǎo)致光子權(quán)重衰減，散射會(huì)改變光子傳播方向，直至光子透射出組織表面或權(quán)重衰減至足夠小的值，光子則被定義為死亡，之后光源再發(fā)射一個(gè)新的光子進(jìn)入組織.其中R為（0～1）隨機(jī)數(shù)，（θ，φ）分別為坐標(biāo)系下極角與方位角，決定光子散射的方向；step為隨機(jī)步長(zhǎng)，由一個(gè)（0～1）隨機(jī)數(shù)和吸收、散射系數(shù)決定；ΔW為在該步長(zhǎng)內(nèi)組織對(duì)權(quán)重的吸收量，由組織的吸收和散射系數(shù)決定.記錄所有光子的傳播信息，便可模擬出光在組織里的傳播過(guò)程，如圖2 所示，圖2 表示源-檢測(cè)距離（Source Detection，SD）=1.5 cm 部分光子傳播路徑x-z軸截面圖.

圖2 x-z軸平面

其中圖例代表光子權(quán)重，可以看出，特定檢測(cè)距離下的光子傳播路徑截面呈現(xiàn)“香蕉”型，光子路徑越遠(yuǎn)，權(quán)重衰減的越多.

2.2 顱內(nèi)血腫建模

硬膜血腫為多發(fā)外傷性顱內(nèi)血腫，顱內(nèi)出血堆積形成血腫塊，血量擴(kuò)散引起血腫厚度變化，位置信息未變，故本研究選用血腫厚度作為研究對(duì)象［13］，硬膜出血多發(fā)位置位于顱骨層下方的腦脊液層，為了更好的模擬人腦的光學(xué)特性，根據(jù)人腦的結(jié)構(gòu)，本研究顱腦結(jié)構(gòu)分為五層，從外到內(nèi)依次為頭皮、顱骨、腦脊液、灰質(zhì)和白質(zhì).血腫模型中將腦脊液層替換為血腫層，模型結(jié)構(gòu)如圖3所示.840 nm波段處，血紅蛋白較其他組織對(duì)近紅外光的吸收高出10倍以上，故本研究選取840 nm為光源的波長(zhǎng)，設(shè)血腫厚度為d，空氣折射率n=1，參數(shù)［14］如表1所示.

表1 頭部模型在波長(zhǎng)為840 nm的光學(xué)參數(shù)

圖3 顱腦血腫模型圖

2.3 仿真模型驗(yàn)證

本研究采用自編蒙特卡羅程序MC自編，為驗(yàn)證程序的準(zhǔn)確性，采取兩種模型的模擬結(jié)果與其他研究者的蒙特卡羅程序模擬結(jié)果相對(duì)比.

模型A：有界組織模型，組織厚度t=0.02 cm，折射率n=1，μa=10 cm-1，μs=90 cm-1，g=0.75，每次模擬50 000個(gè)光子，進(jìn)行5 次計(jì)算取平均值.程序計(jì)算總漫反射率Rd與總透射率Tt，計(jì)算結(jié)果與汪立宏MCML（Monte Carlo Multi-Layered）［15］、范德哈爾斯特（van de Hulst）［16］，普拉爾（Prahl）［17］的數(shù)據(jù)進(jìn)行比較，以MC 仿真結(jié)果Rd、Tt數(shù)據(jù)為基礎(chǔ)，計(jì)算與其他研究學(xué)者仿真結(jié)果的誤差率，計(jì)算如式（1）所示，仿真結(jié)果如表2所示.

表2 模型A對(duì)比數(shù)據(jù)

模型B：半無(wú)界組織模型，μs=90 cm-1，n=1.38，t=10 000 cm，g=0.75，μa取0.5 cm-1、1.0 cm-1、1.5 cm-1、2.0 cm-1、2.5 cm-1，每個(gè)吸收系數(shù)模擬50 000 個(gè)光子，進(jìn)行3 次計(jì)算取平均值與汪立宏MCML 程序計(jì)算結(jié)果進(jìn)行比較，對(duì)比結(jié)果如圖4所示.

表2 數(shù)據(jù)顯示自編程序有界介質(zhì)Rd、Tt仿真結(jié)果較其他研究者誤差均＜0.4%，圖4 顯示半無(wú)界介質(zhì)模型的數(shù)據(jù)曲線幾乎重合，最大誤差0.002 23，最小0.000 21.該對(duì)比數(shù)據(jù)表明本研究的程序數(shù)據(jù)可靠，可以正確的表達(dá)光在生物組織里的傳播.

圖4 模型B與MCML對(duì)比結(jié)果

本文采用的MC自編程序，相較于MCML，自編MC可以追蹤光子在組織里的路徑，記錄光子的實(shí)時(shí)位置及衰減權(quán)重，可以記錄特定檢測(cè)器內(nèi)光子信息，更準(zhǔn)確的分析光在組織里的傳播情況.

3 結(jié)果與分析

3.1 仿真結(jié)果正向分析

仿真模擬從坐標(biāo)原點(diǎn)（0，0，0）發(fā)射107個(gè)光子，在光源橫向位置設(shè)置9 個(gè)檢測(cè)器，SD≤5 cm，檢測(cè)器半徑為0.2 cm，收集每個(gè)檢測(cè)器的光子數(shù)量，定義源-檢測(cè)器靈敏度（Source-Detector Sensitivity，SDS）如式（2）所示：

其中N0為正常腦組織檢測(cè)的光子數(shù)，Nλ表示d=0.2 cm，0.3 cm，…，1.0 cm 情況下檢測(cè)器檢測(cè)的光子數(shù)，仿真結(jié)果如圖5、圖6所示.

圖5 模型仿真下的光子數(shù)據(jù)

圖6 不同血腫厚度下的SDS

由圖5 可得，血腫厚度、SD 與光子數(shù)具有相關(guān)性，光子數(shù)與SD 的數(shù)值關(guān)系呈現(xiàn)指數(shù)下降的趨勢(shì)，d=0 cm與d=0.2 cm，0.3 cm，…，1.0 cm有明顯的光子數(shù)差量值，因此可以有效區(qū)分是否含有血腫.由圖6 可以看出SD與SDS 存在相關(guān)增長(zhǎng)趨勢(shì)，SD＞4.0 cm，SDS 上升速率減緩且趨于平穩(wěn)，驗(yàn)證圖2 的結(jié)論，SD 增加伴隨權(quán)重與光子數(shù)的衰減，對(duì)獲得目標(biāo)層信息貢獻(xiàn)不大，數(shù)據(jù)偶然性將會(huì)增大.d＞0.8 cm，SDS 之間差值縮小，曲線基本重合，說(shuō)明入射的大部分光子已達(dá)到了組織最大的有效深度，血腫厚度增加對(duì)獲得血腫層的信息意義不大.

3.2 數(shù)據(jù)分析及求解

MC 仿真的正問(wèn)題是已知顱腦組織的參數(shù)、血腫位置及厚度信息，研究顱腦表面場(chǎng)域檢測(cè)器內(nèi)的光子分布特性；逆問(wèn)題即對(duì)顱腦中是否含有血腫以及血腫厚度進(jìn)行有效預(yù)測(cè).MC 正問(wèn)題仿真得到的數(shù)據(jù)是理想與無(wú)測(cè)量噪聲的，因此預(yù)測(cè)模型采用仿真數(shù)據(jù).逆問(wèn)題分析中，為了得到低誤差的預(yù)測(cè)值，采用控制變量法可通過(guò)對(duì)已知量的了解來(lái)減少對(duì)未知量估計(jì)的誤差.

由SDS 可知檢測(cè)器得到的光子數(shù)據(jù)，能夠有效反映血腫改變的差異，之間的三維關(guān)系如圖7 所示，整合檢測(cè)器光子數(shù)與SD、血腫厚度的函數(shù)關(guān)系，是評(píng)估顱內(nèi)血腫厚度的關(guān)鍵.

圖7 血腫厚度、SD與光子數(shù)三維關(guān)系圖

圖7顯示，光子數(shù)與SD、血腫厚度具有相關(guān)性，設(shè)D為光源橫向放置的檢測(cè)器，光子數(shù)矩陣函數(shù)N與檢測(cè)器的關(guān)系如式（3）所示：

血腫厚度增大，光子數(shù)有明顯的下降趨勢(shì)，為了更直觀，準(zhǔn)確得到血腫層對(duì)光的衰減情況，設(shè)血腫厚度為變量x，h（x）為血腫厚度與檢測(cè)器有關(guān)光子數(shù)的函數(shù)關(guān)系式，f（Di）與h（x）存在以下關(guān)系：

構(gòu)建h（x），即構(gòu)建正向函數(shù)方程組，本研究根據(jù)h（x）的趨勢(shì)選取兩種函數(shù)模型.

函數(shù)模型1：冪律衰減模型，可以通過(guò)異速關(guān)系從光與血腫組織的相互作用模型中導(dǎo)出［18］，表達(dá)如式（5）所示：

其中，η為比例系數(shù)，u為標(biāo)度指數(shù)，由自變數(shù)P的數(shù)值決定，利用最小二乘法得到最優(yōu)擬合變量［19］，求出殘差平方

Sr對(duì)η，u求導(dǎo)，使導(dǎo)數(shù)等于0，解方程組，可求解η，u.

函數(shù)模型2：指數(shù)律衰減模型，可通過(guò)血腫厚度的變化和光衰減作用關(guān)系得出［20］，表達(dá)如式（7）所示：

其中，P0為補(bǔ)償系數(shù)，A為前置因子，R為衰變速率，采用LM（Levenberg-Marquardt）優(yōu)化算法進(jìn)行非線性曲線擬合［21］，求出殘差平方

將P0、A，R作為變量，使得Sr最小可求得P0、A，R.

將h（x）代入兩種函數(shù)模型中，以Pearson's R、RSquare（COD）作為函數(shù)構(gòu)建的指標(biāo)，Pearson's R 表示變量間線性相關(guān)的程度，R-Square（COD）表示依變數(shù)未知量的變異中占據(jù)的百分比，由控制的自變數(shù)來(lái)解釋.上述SDS可知d＞0.8 cm，大部分入射光子已經(jīng)達(dá)到了最大的有效深度，透射到檢測(cè)器光子數(shù)之間的差值很小，為提高模型的精確度，在構(gòu)建模型時(shí)將不考慮d＞0.9 cm的數(shù)據(jù)，兩種函數(shù)模型構(gòu)建結(jié)果如表3所示.

表3 正向函數(shù)方程組

表3 可得，SD＜5.0 cm，兩種函數(shù)模型Pearson's R、R-Square（COD）均＞0.99，SD=5.0 cm，冪律衰減函數(shù)模型的Pearson's R、R-Square（COD）＜0.99，到達(dá)顱腦頭皮的光子數(shù)變少使數(shù)據(jù)具有較大的偶然性.Pearson's R、RSquare 顯示，函數(shù)矩陣構(gòu)建有效，SD≤5.0 cm，兩種函數(shù)組模型可以表達(dá)光在顱腦組織里傳播.

以上結(jié)果顯示，已知血腫層厚度信息，可通過(guò)函數(shù)矩陣正向求解檢測(cè)器的光子信息，分析光在血腫層的傳播情況.為了通過(guò)檢測(cè)判斷顱腦中的出血情況，本研究基于正向函數(shù)矩陣模型構(gòu)建逆向函數(shù)矩陣，為驗(yàn)證兩種函數(shù)模型的正確性，本研究選取d=0.25 cm，0.35 cm，0.45 cm，0.55 cm，0.65 cm，0.75 cm 進(jìn)行仿真驗(yàn)證，光子數(shù)為107.將不同血腫厚度檢測(cè)到的光子數(shù)代入函數(shù)矩陣中，求解x，x平均值與平均絕對(duì)誤差如表4、表5和圖8所示.

圖8（a）顯示預(yù)測(cè)值與參照值十分接近，血腫厚度＜0.7 cm，冪律衰減模型與指數(shù)律衰減模型的預(yù)測(cè)值與參照值曲線幾乎重合，兩種函數(shù)模型平均絕對(duì)誤差＜3.6%，血腫厚度＞0.7 cm，兩種函數(shù)模型與參照值誤差增大.表4、表5 得，冪律衰減模型計(jì)算值最大誤差0.010 752，最小誤差0.001 160，最大平均絕對(duì)偏差6.473%，最小為1.13%；指數(shù)律衰減模型計(jì)算值最大誤差0.046 886，最小誤差0.000 316，最大平均絕對(duì)偏差10.777%，最小為0.981%，兩種函數(shù)對(duì)比結(jié)果在有效誤差范圍內(nèi).

圖8 結(jié)果對(duì)比分析

表4 冪律衰減函數(shù)驗(yàn)證結(jié)果

表5 指數(shù)律衰減函數(shù)驗(yàn)證結(jié)果

圖8（b）顯示，冪律衰減模型在誤差，平均絕對(duì)誤差以及標(biāo)準(zhǔn)差都要比指數(shù)律衰減模型數(shù)值小，血腫厚度=0.75 cm，冪律衰減模型平均絕對(duì)誤差比指數(shù)律衰減模型＜4.3%，誤差＜6.3886%，預(yù)測(cè)平均值更接參照值，構(gòu)建效果更好，選擇冪率的衰減模型可以更好表達(dá)血腫層的衰減信息，可作為顱腦血腫厚度求解的函數(shù)模型.血腫厚度與誤差具有相關(guān)性，這是因?yàn)樵谙嗤Ｐ拖鹿庾釉诮M織里的路徑增加，能量衰減逐漸增多，透射出頭皮組織表面的光子數(shù)量減少，導(dǎo)致誤差變大.

4 總結(jié)

NIRS 與MC 可模擬光在正常顱腦組織與含有血腫的顱腦組織中的傳播過(guò)程，采用控制變量法改變血腫厚度，對(duì)模型正向仿真，提出血腫檢測(cè)及厚度預(yù)測(cè)的逆問(wèn)題可行性分析方法，分析血腫厚度與SDS，結(jié)果顯示該方法可以有效檢測(cè)顱腦是否含有血腫，檢測(cè)靈敏度與檢測(cè)距離存在對(duì)數(shù)相關(guān)性.利用仿真數(shù)據(jù)建立冪律衰減模型與指數(shù)律衰減模型的正向方程矩陣，可以有效表達(dá)光在血腫層的傳播情況.逆向函數(shù)矩陣的預(yù)測(cè)值與參照值的對(duì)比顯示，在標(biāo)準(zhǔn)差、平均絕對(duì)誤差以及預(yù)測(cè)值等方面，冪律衰減模型比指數(shù)律衰減模型更精確，可作為預(yù)測(cè)并評(píng)估血腫層厚度有效數(shù)學(xué)工具.