金 蒙，高 峰, ，馬文娟，楊 芳，趙會(huì)娟,

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

擴(kuò)散光學(xué)層析成像(diffuse optical tomography，DOT)是一種新興的無(wú)創(chuàng)功能成像模態(tài)[1-5]，其基本原理是：使用波長(zhǎng)范圍在 650～900,nm 的近紅外光“掃描”組織體，通過(guò)測(cè)量組織體表面上的光強(qiáng)分布來(lái)重建組織體內(nèi)部的光學(xué)參數(shù)圖像．DOT的成像對(duì)象是反映組織體病生理功能信息的光學(xué)參數(shù)(吸收和散射系數(shù))，具有靈敏度和特異性高、安全廉價(jià)以及實(shí)時(shí)性強(qiáng)等多重優(yōu)點(diǎn)，因此在新生兒的腦血氧監(jiān)測(cè)、乳腺腫瘤的早期檢測(cè)和基于小動(dòng)物疾病模型的分子成像[5-7]等領(lǐng)域具有廣泛的重要應(yīng)用前景．DOT研究廣泛采用基于模型的圖像重建方法(model-based image reconstruction scheme)[8]，其中擴(kuò)散方程(diffusion equation，DE)和輻射傳輸方程(radiative transfer equation，RTE)為目前應(yīng)用的 2種主要光子輸運(yùn)正向模型．?dāng)U散方程是輻射傳輸方程的球諧函數(shù)一階近似，其適用范圍有相當(dāng)?shù)木窒扌裕簝H適用于高散射、低吸收組織體的遠(yuǎn)源區(qū)光子密度場(chǎng)的描述[1]．因此，一方面，對(duì)于小尺寸組織體(如小動(dòng)物)的成像應(yīng)用，由于光源和探測(cè)器的分布較密，擴(kuò)散方程的適用性較差，影響近源場(chǎng)內(nèi)測(cè)量信息的利用和光學(xué)參數(shù)的重建；另一方面，在生物組織內(nèi)的某些區(qū)域，如高散射、高吸收的骨骼組織、低散射、低吸收的腦脊液層以及腹部大量的空腔域等，擴(kuò)散方程無(wú)法實(shí)現(xiàn)光子傳播行為的有效建模．為此，國(guó)內(nèi)外學(xué)者一直對(duì)基于輻射傳輸方程的 DOT成像技術(shù)有著廣泛的興趣．目前，基于輻射傳輸方程的光子傳播模型普遍基于零內(nèi)向流邊界條件(inflow-zero boundary condition，IZBC)，即假設(shè)邊界處向組織體內(nèi)部的輻射率為零[9-11]．應(yīng)用該邊界條件可以為求解輻射傳輸方程提供極大便利，但由于 IZBC未考慮組織體及其周圍介質(zhì)之間折射率的差異，因此，忽略了光在實(shí)際組織體邊界上的折射和反射，從而造成了正向模型和光的實(shí)際傳播物理過(guò)程之間存在較大誤差．在基于擴(kuò)散方程的擴(kuò)散層析成像中，已經(jīng)發(fā)現(xiàn)了組織體和周圍介質(zhì)的折射率差異對(duì)組織體內(nèi)部光學(xué)參數(shù)的重建有很大的影響[12]．

為了準(zhǔn)確構(gòu)建光在組織體內(nèi)傳播過(guò)程的正向模型，必須考慮采用能夠有效反映邊界反射效應(yīng)的自然邊界條件(natural boundary condition，NBC)．筆者提出了一種采用離散立體角元和有限差分方法求解基于自然邊界條件的輻射傳輸方程，并將數(shù)值求解結(jié)果和蒙特卡洛模擬結(jié)果進(jìn)行了比對(duì)，驗(yàn)證了所提方法的正確性．

1 輻射傳輸方程

在組織體區(qū)域中，精確描述光子輸運(yùn)過(guò)程的穩(wěn)態(tài)輻射傳輸方程[13]為

式中：Φ (r , s)表示空間位置r上、方向s上的輻射率，W/(mm3? sr)；μ和μ分別為散射系數(shù)和吸收系數(shù)，sa為散射相位函數(shù)．本文采用廣泛使用的 Henyey-Greenstein 函數(shù)[14]，即

式中g(shù)為各向異性因子．式(1)中不含源項(xiàng)，考慮邊界某一位置上垂直于邊界向內(nèi)的方向上，存在一準(zhǔn)直光源，用δ函數(shù)可以表示為

式中：0r為光源所在的邊界上的一點(diǎn)；n為位于邊界上源位置處的外法向方向單位向量．

零內(nèi)向流邊界條件為

該邊界條件的物理意義為：當(dāng)光子抵達(dá)組織區(qū)域邊界時(shí)，將從組織體表面徹底溢出，不再返回．顯然，該邊界條件僅僅是一種理想化的邊界條件，其有效性建立在組織體與其周圍媒質(zhì)光折射率相等的基礎(chǔ)上．由于實(shí)際中組織體與其周圍介質(zhì)之間折射率差異明顯，在組織體邊界處必然會(huì)發(fā)生光折射和反射．只有發(fā)生折射的光子能溢出邊界，而反射的光子則返回組織體內(nèi)繼續(xù)傳播，從而影響光子密度分布．顯然，應(yīng)用 IZBC的輻射傳輸方程作為 DOT成像的正向模型勢(shì)必引起圖像重建結(jié)果的顯著誤差．為了正確模擬組織體內(nèi)光子輸運(yùn)行為，發(fā)展基于NBC的輻射傳輸方程建模方法十分必要．

光在到達(dá)邊界處的入射角I?和折射角T?遵循斯涅耳折射定律(見(jiàn)圖1)，表示為

式中：I?為入射角；T?為折射角；R?為反射角；In為模擬組織體的折射率；Tn為組織體包圍介質(zhì)的折射率，本文中所討論的情況僅僅涉及InnΤ≥的情況．

圖1 斯涅耳折射定律Fig.1 Snell’s law of refraction

根據(jù)菲涅爾公式[15]，可以得到反射光的輻射率Φ(r , sR)和入射光輻射率Φ (r , sI)的公式為

式中

稱為菲涅爾反射系數(shù)．式(6)即為自然邊界條件．根據(jù)能量守恒定律，可以得到自然邊界條件下的表面外向光子流量率，即探測(cè)器在組織體表面測(cè)得的光子流率，其公式為

當(dāng)組織與周邊媒質(zhì)光折射率相等時(shí)， 0R= ，于是式(7)第2項(xiàng)為0，從而退化為零內(nèi)向流邊界條件的情形．

2 數(shù)值求解

在筆者以前的工作中，采用了聯(lián)合的立體角元離散法和有限差分法在零內(nèi)向流邊界條件下對(duì)輻射傳輸方程進(jìn)行了求解[16]．對(duì)于自然邊界條件下的輻射傳輸方程進(jìn)行數(shù)值求解需要在此基礎(chǔ)上進(jìn)行改進(jìn)．首先，需要對(duì)自然邊界條件進(jìn)行角度離散，本文對(duì)自然邊界條件采用了離散立體角元的方法進(jìn)行離散．根據(jù)文獻(xiàn)[16]中提出的立體角離散方法，對(duì)于入射光所在的某一立體角元IΩ肯定可以找到與入射面對(duì)稱的反射光所對(duì)應(yīng)的立體角元RΩ，且滿足關(guān)系

假設(shè)在某一離散立體角元IΩ內(nèi)，入射光輻射率為一恒定值，對(duì)式(6)在離散立體角元IΩ內(nèi)做立體角積分，可以得到

對(duì)于任意一個(gè)入射光方向所在的離散立體角元IΩ，可以得到離散的反射系數(shù)為

采用文獻(xiàn)[16]中離散立體角元方法可以求得對(duì)應(yīng)于不同離散立體角元的kR．同樣，對(duì)應(yīng)于自然邊界條件的表面外向光子流率可以離散為

在文獻(xiàn)[16]中提到的采用有限差分方法求解輻射傳輸方程的迭代過(guò)程中，根據(jù)解得的表面外方向的輻射率以及式(11)，可以求解出表面內(nèi)方向上的輻射率，表面內(nèi)方向上的輻射率將影響組織體內(nèi)部光子密度的分布，具體過(guò)程如下．

步驟1 根據(jù)源的位置確定迭代初始值0Φ．

步驟2 第 1l+次迭代，迭代過(guò)程為

步驟3 根據(jù)第1l+次迭代，得到的1lΦ+及式(10)，更新lΦ對(duì)應(yīng)于介質(zhì)表面內(nèi)方向上的輻射率這里表示表面位置上向外方向上的立體角元內(nèi)的輻射率表示與介質(zhì)平面對(duì)稱的立體角元內(nèi)的輻射率．

步驟 4 根據(jù)終止判據(jù)判斷是否終止迭代，如果終止，最后一次迭代的Φ即為最終數(shù)值求解結(jié)果，根據(jù)式(11)可以獲得介質(zhì)表面的表面光子流率．如果未完成迭代則返回步驟2．

根據(jù)上述方法，可以求得有限空間內(nèi)基于自然邊界條件的輻射傳輸方程的數(shù)值解．

3 模擬驗(yàn)證

針對(duì)小動(dòng)物成像問(wèn)題，在大小為 2,cm×2,cm 的二維模擬組織體內(nèi)，對(duì)穩(wěn)態(tài)輻射傳輸方程在 2種邊界條件下進(jìn)行了數(shù)值求解．模擬組織體的光學(xué)參數(shù)為μ= 1 0 mm?1，μ =0.01mm?1，g = 0 .9，分別取3種sa折射率 nI1= 1 .56,nI2= 1 .45,nI3= 1 .38，假設(shè)模擬組織體被空氣包圍，即 nT= 1 .00．光源位于模擬組織體一邊的中心，定義在光源對(duì)面的面為透射面，在光源側(cè)面為側(cè)面，如圖 2所示．為了驗(yàn)證求解輻射傳輸方程算法的正確性，對(duì)Wang等[17]提出的無(wú)限半空間內(nèi)的光子傳播蒙特卡洛模擬程序進(jìn)行了擴(kuò)展，加入了二維邊界，當(dāng)光子溢出邊界時(shí)，記錄光子的溢出位置，并將相同位置上的光子權(quán)重進(jìn)行累加，獲得該位置上的表面光子流率．通過(guò)蒙特卡洛模擬得到的側(cè)面和透射面上的表面光子流率，如圖 3所示．同時(shí)，根據(jù)輻射傳輸方程數(shù)值求解的結(jié)果以及式(11)獲得模擬組織體的表面光子流率，然后將結(jié)果和蒙特卡洛模擬結(jié)果進(jìn)行了比對(duì)，比對(duì)結(jié)果如圖4所示．從圖3中可以看到，當(dāng)模擬組織體的折射率從 1.38增加到 1.56的過(guò)程中，在側(cè)面上和透射面上得到的表面光子流率隨著折射率的增加而增加．當(dāng)模擬組織體折射率為 1.56時(shí)，在側(cè)面上的表面光子流率比模擬組織體折射率為1.45時(shí)約高6.24%，比模擬組織體折射率為1.38時(shí)約高 11.48%，比應(yīng)用零內(nèi)向流邊界條件時(shí)約高62.55%．在透射面上的表面光子流率也有類似的分布，當(dāng)模擬組織體折射率為1.56時(shí)，透射面上的表面光子流率比模擬組織體折射率為1.45時(shí)約高7%，比模擬組織體折射率為 1.38時(shí)約高 12.37%，比應(yīng)用零內(nèi)向流邊界條件約高 62.55%．值得一提的是，在圖3(a)中，應(yīng)用零內(nèi)向流邊界條件的蒙特卡洛模擬得到的側(cè)面的表面光子流率的峰值相對(duì)于應(yīng)用自然邊界條件的蒙特卡洛模擬得到的側(cè)面的表面光子流率的峰值在y軸方向上相對(duì)于源的位置遠(yuǎn)大約1,mm．

圖4為數(shù)值求解輻射傳輸方程所得之透射面和反射面外向光子流率和蒙特卡洛模擬結(jié)果的比對(duì)結(jié)果，從圖4中可以看到，在考慮組織內(nèi)外折射率差異的自然邊界條件及零內(nèi)向流邊界條件下數(shù)值求解結(jié)果分別與其相應(yīng)蒙特卡洛模擬結(jié)果相吻合，從而證明了本文所提出的應(yīng)用離散立體角元和有限差分方法求解自然邊界條件下的輻射傳輸方程的算法的正確性．

圖2 模擬仿體形狀和光源位置Fig.2 Shape of simulated phantom and source position

圖3 根據(jù)蒙特卡洛模擬得到的表面光子流率Fig.3 Outgoing photon density acquired from MC simulations

圖5所示為在零內(nèi)向流邊界條件和自然邊界條件 ( nI= 1 .56,nT= 1 .00)下求解輻射傳輸方程所得的模擬組織體內(nèi)部光子密度分布．可以看出，自然邊界條件下所得到的組織體內(nèi)部的光子密度比零內(nèi)向流邊界條件下的值平均高出 178.26%，因此，不同的邊界條件對(duì)組織體內(nèi)部的光子密度分布影響極大．

圖4 蒙特卡洛模擬和輻射傳輸方程數(shù)值解得到的歸一化表面光子流率比對(duì)Fig.4 Comparisons of normalized outgoing photon density acquired from MC simulations and FDM solution to RTE with NBC and IZBC

圖5 不同邊界條件下點(diǎn)光源激勵(lì)下均勻組織體內(nèi)的光子密度分布Fig.5 Comparisons of photon density inside tissue acquired from FDM solution to RTE with IZBC

本文中模擬了 2個(gè)非均勻模擬組織體在如圖 2所示的準(zhǔn)直光源激勵(lì)的情況下，對(duì)輻射傳輸方程在零內(nèi)向流邊界條件下和自然邊界條件下進(jìn)行了數(shù)值求解．非均勻模擬組織體內(nèi)部的吸收系數(shù)和散射系數(shù)分布分別如圖 6(a)～6(d)所示．零內(nèi)向流條件下獲得的模擬組織體內(nèi)部的光子密度分布分別如圖 6(f)、(h)所示，自然邊界條件下得到的模擬組織體內(nèi)部的光子密度分布如圖 6(e)、(g)所示．在 2種不同邊界條件下得到的 2個(gè)模擬組織體表面上的表面光子流率分布如圖7所示．

圖6 非均勻組織體內(nèi)部光學(xué)參數(shù)分布以及點(diǎn)光源激勵(lì)下內(nèi)部光子密度分布Fig.6 Optical coefficients distribution of Tissue A and B,and photon density inside Tissue A and B with different boundary conditions

圖7 不同邊界條件下非均勻組織體 A和 B在點(diǎn)光源激勵(lì)下得到的表面光子流率Fig.7 Outgoing photon density on side surface and transmission surface of Tissue A and B with different boundary conditions

4 結(jié) 語(yǔ)

提出了一種在自然邊界條件下求解輻射傳輸方程的離散立體角元-有限差分方法．通過(guò)對(duì)小動(dòng)物尺寸二維組織域進(jìn)行求解，并與相應(yīng)的蒙特卡洛模擬結(jié)果進(jìn)行比較，證明了該算法的可行性和有效性．在此基礎(chǔ)上，進(jìn)一步探討了不同邊界條件對(duì)基于輻射傳輸方程的 DOT成像正向模型所產(chǎn)生的影響．研究表明，不同邊界條件下的正向模型數(shù)值解差異明顯，因此，在今后的基于輻射傳輸方程的 DOT逆問(wèn)題研究中應(yīng)采用更符合實(shí)際的自然邊界條件．

［1］ Arridge S R. Optical tomography in medical imaging[J]．Inv Probl，1999，15(2)：41-93.

［2］ Yodh A，Chance B. Spectroscopy and imaging with diffusing light[J]. Phys Today，1995，48：34-40.

［3］ Arridge S R，Hebden J C. Optical imaging in medicine(Part 2)：Modelling and reconstruction[J]. Phys Med Biol，1997，42(5)：841-853.

［4］ Gibson A，Hebden J C，Arridge S R. Recent advances in diffuse optical imaging[J]. Phys Med Biol，2005，50(4)：1-43.

［5］ Zhao Huijuan，Gao Feng，Yukari Tanikawa，et al.Time-resolved diffuse optical tomographic imaging for the provision of both anatomical and functional information about biological tissue[J]．Appl Opt，2005，44(10)：1905-1916.

［6］ Gao Feng，Zhao Huijuan，Tanikawa Yukari，et al.Optical tomographic mapping of cerebral haemodynamics by means of time-domain detection：Methodology and phantom validation[J]．Phys Med Biol，2004，49(6)：1055-1078.

［7］ Zhao H，Gao F，Tanikawa Y，et al. Time-resolved Diffuse optical tomography and its application to in vitro and vivo imaging[J]．Journal of Biomedical Optics，2007，12(6)：No. 062107.

［8］ Klose A D，Uwe N，Jürgen B，et al. Optical tomography using the time-independent equation of radiative transfer(Part 1)：Forward model[J]．Journal of Quantitative Spectroscopy and Radiative Transfer，2002，72(5)：691-713.

［9］ Kim Hyun Keol，Charette A. Frequency domain optical tomography using a conjugate gradient method without line search[J]．Journal of Quantitative Spectroscopy and Radiative Transfer，2007，104(2)：248-256.

［10］ Balima O，Pierre T，Charette A，et al. A least square finite element formulation of the collimated irradiation in frequency domain for optical tomography applications[J]．Journal of Quantitative Spectroscopy and Radiative Transfer，2010，111(2)：280-286.

［11］ Tarvainen Tanja，Vauhkonen Marko，Kolehmainen Ville，et al. Coupled radiative transfer equation and diffusion approximation model for photon migration in turbid medium with low-scattering and non-scattering regions[J]．Phys Med Biol，2005，50(20)：4913-4930.

［12］ Chen Bingquan，Stamnes K，Jakob S J. Validity of the diffusion approximation in bio-optical imaging [J]．App Opt，2001，40(34)：6356-6366.

［13］ Duderstadt J J，Martin W R. Transport Theory[M].New York：Wiley，1979.

［14］ Henyey L G，Greenstein J L. Diffuse radiation in the galaxy[J]．Annales d'Astrophysique，1940，3：117-137.

［15］ Hecht E. Optics [M]．2nd ed. USA：Addison and Wesley Publishing Company，Inc，1987.

［16］ 金 蒙，高 峰，李 嬌，等. 二維穩(wěn)態(tài)輻射傳輸方程的有限差分求解法[J]．光子學(xué)報(bào)，2010，39(9)：1594-1601.Jin Meng，Gao Feng，Li Jiao，et al. A finitedifference-method solution to two-dimensional steadystate radiative transfer equation[J]． Acta Photonica Sinica，2010，39(9)：1594-1601(in Chinese).

［17］ Wang Lihong，Jacques S L，Zheng Liqiong. MLMCMonte Carlo modeling of light transport in multi-layered tissues[J]．Computer Methods and Programs in Biomedicine，1995，47(2)：131-146.