基于自然邊界條件求解輻射傳輸方程的方法

金 蒙，高 峰, ，馬文娟，楊 芳，趙會娟,

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

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

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

1 輻射傳輸方程

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

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

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

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

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

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

光在到達邊界處的入射角I?和折射角T?遵循斯涅耳折射定律(見圖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ù)能量守恒定律，可以得到自然邊界條件下的表面外向光子流量率，即探測器在組織體表面測得的光子流率，其公式為

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

2 數(shù)值求解

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

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

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

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

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

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

步驟2 第 1l+次迭代，迭代過程為

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

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

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

3 模擬驗證

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

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

圖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%，因此，不同的邊界條件對組織體內(nèi)部的光子密度分布影響極大．

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

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

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

圖6 非均勻組織體內(nèi)部光學參數(shù)分布以及點光源激勵下內(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在點光源激勵下得到的表面光子流率Fig.7 Outgoing photon density on side surface and transmission surface of Tissue A and B with different boundary conditions

4 結(jié) 語

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

［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]．光子學報，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.