史振志，趙會娟，陳文亮，徐可欣

(天津大學精密儀器與光電子工程學院，天津 300072)

在體生物組織光學參數(shù)測量以及光子在生物組織中的傳輸問題是生物醫(yī)學光子學領域的重要基礎研究課題，其在光動力療法[1]、無創(chuàng)血糖檢測[2-3]、癌癥診斷[4-5]等方面都有著潛在的應用．對漫反射光的測量可采取時間分辨、頻域和穩(wěn)態(tài)(或稱連續(xù)光，continuous wave，CW)空間分辨測量方法．基于在體穩(wěn)態(tài)測量的光學參數(shù)反構的一般原理是依據漫反射光和生物組織光學參數(shù)之間的正向模型，將空間多點測量得到的漫反射光代入反構算法得到生物組織光學參數(shù)的估計值[6-9]，因此發(fā)展快速準確的光子傳輸模型和反構方法顯得尤為重要．

通常，輻射傳輸方程(radiative transfer equation，RTE)被認為是描述光子在生物組織中傳輸?shù)淖顪蚀_的模型，但一般來說，很難直接得到 RTE的解析解，所以研究者們發(fā)展了離散坐標法(SN)和球諧函數(shù)法(PN)等一些近似的解析解法．漫射近似(diffusion approximation，DA)作為一階 PN 近似，是一種成功的近似模型，在生物醫(yī)學領域得到了廣泛的應用，然而此模型基于一些限制條件，只適用于大約化反照率( 'a)和大檢測距離(source-detector separations，SDSs)的情況．隨著近紅外光譜在生物醫(yī)學中的應用，需要新的正向模型來描述小約化反照率近光源(小檢測距離)情況下的光子傳輸模型，例如在1,000～2,000,nm 波長的無創(chuàng)血糖檢測的應用中，由于該波段水和皮膚組織的高吸收特性，約化反照率比較?。徊⑶覟榱颂岣邫z測信噪比通常在滿足探測深度的條件下采用近光源檢測，所以DA不適合這樣的應用．為了解決該問題，Hull和 Foster用球諧函數(shù)法推導了輻射傳輸方程，給出了 P3近似的格林方程，并證明了在高吸收介質近光源檢測的情況下，P3近似比DA更加準確．為了修正DA和簡化P3近似以適用于高吸收的介質，他們在DA的基礎上結合P3近似發(fā)展了一種混合漫射近似(hybrid diffuse approximation，DAH)模型[10]．Tian 等[11]曾使用 DAH 模型來反構光學參數(shù)．另外，Klose等[12]提出了一種簡化球諧函數(shù)(SPN)模型，采用有限元方法近似求解復雜的 RTE，并且證明了 SPN近似能夠顯著提高光子在高吸收介質小體積物體上的傳輸精度．

本文在大約化反照率范圍(0.50～0.99)以及小檢測距離(0.4～8.0,mm)下，研究了用不同的正向模型反構吸收系數(shù)(aμ)和約化散射系數(shù)('sμ)的準確性．在介紹了 6種光子傳輸模型(單點源和雙點源模型下的 DA、DAH和 P3近似)以及反構光學參數(shù)的方法后，筆者對6種正向模型的適用范圍進行數(shù)值模擬和分析，提出了一種使用不同正向模型聯(lián)合反構光學參數(shù)的方法，給出了在5個不同約化反照率區(qū)間反構aμ和sμ′的最佳正向模型以及相應的反構誤差．

1 方 法

1.1 正向模型

當忽略波動效應時，光子在生物組織中的傳輸可用RTE描述為

式中：L ( r,s?)為r處s?方向的平均光功率密度；q(r,s?)為光源分布；μs為散射系數(shù)；p(s? ?s?′ )為生物組織的散射相函數(shù)，描述粒子從s?′方向到s?方向散射的概率．p ( s? ?s?′ )通常用 Henyey-Greenstein相位函數(shù)[13]來描述為

式中：g為各向異性因子；θ為從?s′方向到?s方向的角度．

球諧函數(shù)法是用球諧函數(shù)作為基函數(shù)對 RTE進行展開，再對展開的 PN方程組進行求解，得到 PN近似．DA一般也稱為P1近似，是在1階PN近似方程的基礎上進行近似簡化求解得到．在無限介質以及各向同性穩(wěn)態(tài)光源的情況下，DA的格林函數(shù)解可以表示為

式中：Φ為輻射通量；ρ為探測距離；effμ為有效衰減系數(shù)；D為漫射系數(shù)；sμ′為約化散射系數(shù)．

P3近似是由RTE展開的PN方程組的3階近似解，理論上考慮了相函數(shù)的2階矩和3階矩，能夠更加準確地描述組織的光輻射分布．因為篇幅的限制，這里不給出 P3近似具體的推導過程，具體可以參照文獻[10]．

DAH的格林函數(shù)解是在DA格林函數(shù)解的基礎上，分別用P3近似里面的v?和asymD 系數(shù)代替DA中的effμ和D[10]．對應方程(3)，DAH的格林函數(shù)解可表示為

其中

式中 g1、g2和 g3分別為相位函數(shù)的 1階矩、2階矩和 3階矩．根據 Henyey-Greenstein相位函數(shù)，可得γ=1 + g ，δ =1 + g + g2．

在半無限介質的情況下，DA和 DAH的推導主要涉及光源的特征項以及合適的邊界條件．

對于光源項，由于實際光源很復雜，將其直接用到 RTE中去求解非常困難，因此通常將光源進行簡化處理．在半無限散射介質中，通常將光源看成一無限細的筆形光源，考慮該光束由z軸入射到半無限的散射介質中，有效的光源分布[14]可表示為

式中a′為約化反照率.

為了簡化無限細光源的分布特征，通常采用簡單的各向同性點源來代替無限細筆形光源的分布．采用一個與空氣-組織界面具有相同偶極矩的光源來近似有效光源項，稱之為單點源，如圖 1(a)所示；而采用同時滿足具有相同偶極矩和四極矩的 2種離散點源分布來近似有效光源項，稱之為雙點源[10-11]，如圖1(b)所示．

對于單點源近似，考慮位于0z處單點源的偶極矩與式(6)表示的有效光源的偶極矩相等，可得

式(7)表示無限細光源可以用位于 z0= 1 /μt′的各向同性點源來表示，其中光源強度為約化反照率a′．

對于雙點源近似，將式(7)的有效光源處理為 z01和 z02處的2個強度都為 a ′/2 的點源．考慮無限細光源和雙點源的四極矩和偶極矩相等，得

式中有效的光源深度為 z01= 2 /μt′和z02= 0 ．

對于半無限介質的邊界的處理，本文采用外推邊界條件來描述光在空氣-組織界面的輻射通量．外推邊界條件假設輻射通量在物理界面上不為零，而是在介質表面外 zb處為零(如圖 1所示)[15-16]，外推邊界可表示為

圖1 點源近似以及外推邊界條件示意Fig.1 Approximation for point sources and the extrapolated boundary condition

其中

式中effR 為有效反射系數(shù)．

為了滿足外推邊界條件，必須引入鏡像光源，其強度與真實的光源相同，但處于外推邊界的另一側．單點源的情況下，鏡像光源的位置為z=? ( 2 zb+ z0)；雙點源的情況下，鏡像光源的位置分別為 z1=? ( 2 zb+ z01)和 z2=? ( 2 zb+ z02)．因此，半無限介質在外推邊界條件下的單點源光源項 q1(z)和雙點源光源項 q2(z)可以用真實光源和鏡像光源之和表示，即

根據以上的討論，在外推邊界條件下 DA和DAH的輻射通量可以表示[17]為

界面的漫反射率[15]可以表示為

式中fres()R θ為菲涅耳反射系數(shù)．

1.2 逆問題

為了研究在空間分辨漫反射測量條件下 RTE的不同解對光學參數(shù)反構的影響，用蒙特卡羅模擬產生空間分辨漫反射測量數(shù)據，而蒙特卡洛模擬中輸入的吸收系數(shù)、散射系數(shù)、各向異性因子和折射系數(shù)作為光學參數(shù)的“真實值”．蒙特卡羅模擬采用了目前最廣泛使用的 MCML程序[18]，用于多層半無限介質穩(wěn)態(tài)光的模擬傳輸，程序采用了無限細的筆形光束作為光源，用 Henyey-Greenstein相位函數(shù)來計算散射角度．

反構中，分別使用單點源和兩點源模型下的DA、DAH、P3近似的解通過 Levenberg-Marquadt(LM)非線性擬合對“測量值”進行擬合，反構結果與“真實值”之間的相對誤差作為評價模型適應性的標準．

2 結 果

在蒙特卡羅模擬中采用了 0.1,mm的空間分辨率，最大的記錄半徑為 30,mm，使用了 108的光子進行運算．P3近似的推導基于 Hull等[10]的研究結果，程序代碼由Finlay教授提供．

為了驗證各種模型對不同約化反照率下光學參數(shù)的重構精度，在蒙特卡羅模擬中保持 μs'不變但使吸收系數(shù)從小到大變化，采用的吸收系數(shù)為 μa＝0.01,mm-1、0.1,mm-1、0.5,mm-1、1.0,mm-1，其他光學參數(shù)為sμ′＝1,mm-1，g＝0.9，n＝1.4．蒙特卡羅模擬中半無限介質的空間分辨探測如圖 2(a)所示，探測的范圍 ρ為 ρmin～8.0,mm，間隔 0.1,mm．ρmin代表在進行空間分辨測量時所采用的最小探測距離，ρmin的范圍為 0.4～4.0,mm，間隔是 0.1,mm．對同一組光學參數(shù)在一系列 ρ下蒙特卡羅模擬得出的“測量值”與使用單點源和兩點源模型下的DA、DAH和P3近似的解進行非線性擬合，得到光學參數(shù)的重構值，光學參數(shù)的反構原理如圖2(b)所示．

圖 2 蒙特卡羅模擬中半無限介質的空間分辨探測及組織光學參數(shù)反構示意Fig.2 Scheme of steady-state spatially resolved diffuse reflectance of Monte Carlo simulation under semiinfinite medium and the derivation of optical properties

圖 3～圖 6給出了使用上述 6種模型在不同 a '下反構μa和μs′的結果．圖3～圖6的結果表明：光學參數(shù)的準確計算取決于正向模型的選擇，并且不同的ρmin影響μa和μs′的反構結果．當a'= 0 .99時(如圖 3所示)，反構μa最佳的正向模型為雙點源的DA和DAH，在最小探測距離 ρmin為 0.6～4.0,mm 的情況下，反構誤差＜10%；而反構μs′最佳的正向模型為單點源的 DA 和 DAH，在 ρmin為 0.4～4.0,mm 的情況下，反構誤差在0.6%～7.0%之間．當 a '= 0 .91時(見圖4)，反構μa和μs′最佳的正向模型都為單點源的DAH，在 ρmin為0.7～4,mm的情況下，μa的反構誤差＜10%；在 ρmin為1.7～4.0,mm的情況下，μs′的反構誤差＜10%．當 a'= 0 .67時(圖 5)，雙點源的 DAH模型優(yōu)于其他的模型，在ρmin為0.4～2.6,mm的情況下，μa的反構誤差＜5%；在ρmin為0.5～2.6,mm的情況下，μs′的反構誤差＜5%．當a'=0.5時(圖6)，單點源DA適合反構μa，在 ρmin為0.7～2.4,mm的情況下，反構誤差＜6%；對于μs′，這些模型的反構準確度都不高，雙點源的 DAH 模型的誤差最小，在 ρmin為0.4～1.2,mm 的情況下，反構誤差僅為 20%左右．但是，無論是單點源的 P3近似還是雙點源的 P3近似在光學參數(shù)的計算中都沒有得到好的結果．這是因為雖然 P3是高階近似，并且已經證明了在高吸收介質近光源檢測的情況下P3近似比DA更加準確，但P3近似也不是任何時候都比DA更加準確．

圖3 a′=0.99時吸收系數(shù)和約化散射系數(shù)的反構結果Fig.3 Derived absorption coefficient and reduced scattering coefficient at a′=0.99

圖4 a′=0.91時吸收系數(shù)和約化散射系數(shù)的反構結果Fig.4 Derived absorption coefficient and reduced scattering coefficient at a′=0.91

圖5 a′=0.67時吸收系數(shù)和約化散射系數(shù)的反構結果Fig.5 Derived absorption coefficient and reduced scattering coefficient at a′=0.67

圖6 a′=0.50時吸收系數(shù)和約化散射系數(shù)的反構結果Fig.6 Derived absorption coefficient and reduced scattering coefficient at a′=0.50

為了更深入地研究以上模型的適用性，將 a '為0.50～0.99的區(qū)間進一步細分，使用 19組不同的吸收 系 數(shù) (μa＝ 0.01,mm-1，0.02,mm-1，0.03,mm-1，0.04,mm-1，0.05,mm-1，0.06,mm-1，0.07,mm-1，0.08,mm-1，0.09,mm-1，0.10,mm-1，0.20,mm-1，0.30,mm-1，0.40,mm-1，0.50,mm-1，0.60,mm-1，0.70,mm-1，0.80,mm-1，0.90,mm-1，1.00,mm-1)，設μs′=1,mm-1，g＝0.9，n＝1.4，形成了 5 個約化反照率的區(qū)間 ，即0.96 ≤ a '≤ 0 .99的 “ 樣 品 ” 有 4個(μa＝0.01,mm-1，0.02,mm-1，0.03,mm-1，0.04,mm-1)，0.91≤a'＜0.96的樣品有7個(μa＝0.04,mm-1，0.05,mm-1，0.06,mm-1，0.07,mm-1，0.08,mm-1，0.09,mm-1，0.10,mm-1)，0.71 ≤ a '＜0 . 91的樣品有4個(μa＝0.10,mm-1，0.20,mm-1，0.30,mm-1，0.40,mm-1)，0.62≤a'＜0.71的樣品有3個(μa＝0.40,mm-1，0.50,mm-1，0.60,mm-1)，0.5≤a'＜0.62的樣品有5個(μa＝0.60,mm-1，0.70,mm-1，0.80,mm-1，0.90,mm-1，1.00,mm-1)．

對各個模型反構結果進行比較，得到了在不同約化反照率范圍條件下最佳的正向模型，如表1所示．

表1 在不同約化反照率情況下反構吸收系數(shù)和約化散射系數(shù)的最佳模型Tab.1 The best model to derive absorption coefficient and reduced scattering coefficient together at different ranges of reduced albedo

0.96 ≤ a'≤ 0.99下反構μa的最佳正向模型為雙點源的混合漫射近似，反構μs′的最佳正向模型為單點源的混合漫射近似；在0.91≤a'＜0.96下，反構μa和μs′的正向模型同為單點源的混合漫射近似；在0.71≤a'＜0.91下，反構μa的最佳正向模型為單點源的漫射近似，反構μs′的最佳正向模型為雙點源的漫射近似；在0.62≤a'＜0.71下，反構μa和μs′的正向模型同為雙點源的混合漫射近似；在0.50'0.62a≤＜下，反構aμ的最佳正向模型為單點源的漫射近似，反構sμ′的最佳正向模型為雙點源的混合漫射近似．

根據以上結果，提出一種計算光學參數(shù)的聯(lián)合反構方法，和以往反構aμ和sμ′時使用同一種正向模型不同的是，聯(lián)合反構算法不僅在不同約化反照率時使用不同的正向模型，而且反構aμ和sμ′時也使用不同的正向模型．例如當μa＝0.7,mm-1，sμ′＝1,mm-1時，如果只是使用傳統(tǒng)的單一的正向模型(如單點源的漫射近似)，吸收系數(shù)的反構誤差為 0.14%，約化散射系數(shù)的反構誤差僅為 49.1%；使用了聯(lián)合反構算法，吸收系數(shù)的反構誤差因為同為單點源漫射近似所以同為 0.14%，而約化散射系數(shù)的反構誤差為 10.4%，有了很大的提高．

在光學參數(shù)的反構問題中，最小二乘問題是求方程的平方和的最小化，而正向問題的準確性是反構問題的關鍵．從圖 7可以看出，各個理論模型的漫反射率與模擬結果的偏差各不相同，尤其在本文所選擇的最小檢測距離的范圍(0.4,mm≤ρmin≤4.0,mm)下，變化的浮動較大，而隨著光源探測的距離減小，漫反射率的值也越大．所以選擇不同的最小檢測距離 ρmin，使得方程的平方和的值變化浮動增大，從而在最小化的擬合過程中，使得反構結果的誤差增大，從而表明了合理選擇ρmin的重要性．為了將這一結論對實際的光學參數(shù)計算具有指導意義，研究了最小探測距離ρmin的選取范圍．使用表 1給出的最佳正向模型，表2給出了對19組光學參數(shù)進行反構時最佳的ρmin范圍，并給出了相應的反構誤差范圍．結果表明，在光學參數(shù)計算的過程中，選擇合適的最小探測距離 ρmin對提高光學參數(shù)的反構精度有很大的幫助．

圖7 理論模型與模擬結果的漫反射率的偏差Fig.7 Variety of diffuse reflectance of theoretical models and simulation results

表2 19組光學參數(shù)條件下最優(yōu)的最小探測距離ρmin范圍以及相應的反構誤差Tab.2 ρmin and the corresponding derivation error range in case of 19 kinds of e derived absorption coefficient and reduced scattering coefficient

3 結 語

研究了在近光源(光源檢測器范圍0.4～8.0,mm)和大的約化反照率范圍(0.50～0.99)的情況下吸收系數(shù)和約化散射系數(shù)的反構問題，使用了6種光學傳輸模型對光學參數(shù)的計算進行了討論．為了細致地研究各個正向模型的適用性，在約化反照率范圍為0.50～0.99的情況下，選取了 19組光學參數(shù)通過對模擬數(shù)據的擬合來進行光學參數(shù)的計算．由于篇幅的原因，本文沒有給出所有的計算結果，只給出了在5個約化反照率范圍用于反構aμ和sμ′的最佳模型．從結果中可以看出光學參數(shù)的準確計算取決于正向模型的正確選擇，并且aμ和sμ′的計算都有各自的適用模型，在此 'a范圍內未必使用高階的PN近似就能取得好的結果，而DA和DAH適用于這些 'a范圍內的aμ和sμ′重建．另外還研究了最小探測距離ρmin的選取范圍并給出了相應的反構誤差，得出不同的ρmin值的選擇影響aμ和sμ′的反構精度．

根據計算和分析，在近光源(光源檢測器范圍0.4～8.0,mm)和大的約化反照率范圍(0.50～0.99)的情況下，提出了計算光學參數(shù)的聯(lián)合反構方法，即在不同的約化反照率下，使用不同的正向模型來分別反構吸收系數(shù)和約化散射系數(shù)，如表 1所示．如果在重構之前或重構過程中已知 'a的范圍，則可根據 'a選擇最佳的正向模型，從而達到降低重構誤差的目的．但在實際的光學參數(shù)測量過程中，光學參數(shù)的值是未知的，所以無法判斷 'a的值．故將聯(lián)合反構算法應用于實際之前，需要構建一個正確選擇正向模型的預判體系．本文主要介紹光學參數(shù)聯(lián)合反構的方法，而這個正確選擇正向模型的預判體系會在結合具體實際應用的過程中進一步去進行研究．

本文只討論了在sμ′固定(同為1,mm-1)的情況下各個正向模型在不同 'a下的適用性，未來會研究由不同的sμ′所導致的不同'a范圍內正向模型的適用性，從而得到更加完善的光學參數(shù)聯(lián)合反構方法．

［1］ Wilson B C，Patterson M S. The physics，biophysics and technology of photodynamic therapy[J]. Physics in Medicine and Biology，2008，53(9)：R61-R109.

［2］ Bykov A V，Kirillin M Y，Priezzhev A V，et al. Simulations of a spatially resolved reflectometry signal from a highly scattering three-layer medium applied to the problem of glucose sensing in human skin[J]. Quantum Electronics，2006，36(12)：1125-1130.

［3］ Yang Yue，Chen Wenliang，Shi Zhenzhi，et al. The reference point of floating-reference method for blood glucose sensing[J]. Chinese Optical Letter，2010，8(4)：421-424.

［4］ Zhu C，Palmer G M，Breslin T M，et al. Use of a multiseparation fiber optic probe for the optical diagnosis of breast cancer[J]. Journal of Biomedical Optics，2005，10(2)：024032.

［5］ Pavlova I，Weber C R，Schwarz R A，et al. Monte Carlo model to describe depth selective fluorescence spectra of epithelial tissue：Applications for diagnosis of oral precancer[J]. Journal of Biomedical Optics，2008，13(6)：064012.

［6］ Doornbos R M P，Lang R，Aalders M C，et al. The determination of in vivo human tissue optical properties and absolute chromophore concentrations using spatially resolved steady-state diffuse reflectance spectroscopy[J].Physics in Medicine and Biology，1999，44(4)：967-981.

［7］ Arifler D. Sensitivity of spatially resolved reflectance signals to coincident variations in tissue optical properties[J]. Applied Optics，2010，49(20)：4310-4320.

［8］ Cen H，Lu R. Quantification of the optical properties of two-layer turbid materials using a hyperspectral imagingbased spatially-resolved technique[J]. Applied Optics，2009，48(29)：5612-5623.

［9］ Pilz M，Honold S，Kienle A. Determination of the optical properties of turbid media by measurements of the spatially resolved reflectance considering the pointspread function of the camera system[J]. Journal of Biomedical Optics，2008，13(5)：054047.

［10］ Hull E L，F(xiàn)oster T H. Steady-state reflectance spectroscopy in the P3 approximation[J]. Journal of the Optical Society of America A，2001，18(3)：584-599.

［11］ Tian Huijuan，Liu Ying，Wang Lijun，et al. Hybrid diffusion approximation in highly absorbing media and its effects of source approximation[J]. Chinese Optical Letter，2009，7(6)：515-518.

［12］ Klose A D，Larsen E W. Light transport in biological tissue based on the simplified spherical harmonics equations[J]. Journal of Computational Physics，2006，220：441-470.

［13］ Henyey L G，F(xiàn)reyer J L. Diffuse radiation in the galaxy[J]. The Astrophysical Journal，1941，93：70-83.

［14］ Farrell T J，Patterson M S，Wilson B C. A diffusion theory model of spatially resolved，steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo[J]. Medical Physics，1992，19(4)：879-888.

［15］ Kienle A，Patterson M S. Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium[J].Journal of the Optical Society of America A，1997，14(1)：246-254.

［16］ Haskell R C，Svaasand L O，Tsay T，et al. Boundary conditions for the diffusion equation in radiative transfer[J]. Journal of the Optical Society of America A，1994，11(10)：2727-2741.

［17］ Tualle J M，Prat J，Tinet E，et al. Real-space Green's function calculation for the solution of the diffusion equation in stratified turbid media[J]. Journal of the Optical Society of America A，2000，17(11)：2046-2055.

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