生物组织的δ-P1近似漫反射光学模型*

刘迎刘小君 齐贝贝田会娟

生物组织的δ-P1近似漫反射光学模型*

刘迎1)刘小君2)齐贝贝1)田会娟3)

1)(天津大学理学院，天津300072)
2)(精密仪器与光电子工程学院，光电信息技术科学教育部重点实验室，天津300072)
3)(天津工业大学，信息与通讯工程学院，天津300160)
(2010年8月2日收到;2010年9月2日收到修改稿)

本文研究了生物组织的改进的δ-P1近似漫反射光学模型，推导了含有等效光源一阶矩的双点源近似空间分辨漫反射解Rδ－P1(ρ)．研究表明，考虑等效光源一阶矩的光学模型，较好地描述了具有强的前向散射特性和较大吸收系数的生物组织散射特性;与漫射近似下的漫反射率RSDA(ρ)相比，新的光学模型能较好地描述光源附近的漫辐射强度分布，并且由于解析表达式中含有散射相函数的二阶参量γ，这对于在近光源区域由空间分辨漫反射测量数据来反演生物组织的光学参量具有重要的意义．

组织光学，δ-P1近似，等效光源，微区漫反射

PACS:42．25．Dd，42．62．Be，87．64．Cc

1.引言

对混浊介质中光子输运问题的研究为漫射光学光谱技术提供了理论基础．Boltzmann输运方程已经广泛地用于研究光在生物组织中的传播问题，它的解描述了辐射通量(fluence)、反射率(reflectance)或辐射强度(radiance)等可测量物理量与生物组织的光学参量之间的关系［1］．在输运理论中，生物组织的光学性质用吸收系数μa、散射系数μs和散射相函数p(μ)描述．在实际问题中，散射相函数p(μ)的复杂性，限制了对Boltzmann输运方程的精确求解［2］．为了得到解析的解，通常将Boltzmann输运方程以球谐函数的形式展开成为一组封闭的方程组，P1近似只取零阶和一阶项［3］．在P1近似基础上，如果在输运方程球谐函数的形式展开式中，光源项只取到零阶，这就是漫射近似．对于漫反射问题，通常要求光源与探测器之间的距离远大于一个输运平均自由程，对于生物组织，一般需要大于几个毫米．这个条件限制了漫射近似的应用，例如，当被测组织的体积比较小或组织的吸收比较大时，需要测量光源附近组织的辐射强度分布［4］;当对表层组织测量时，需要研究散射次数较少的非漫射近似光子行为［5］．

近些年对解析理论模型的研究主要是围绕克服漫射近似的这些缺点进行的．一些研究者采用不同方法来扩大漫射近似的适用范围，例如对漫射系数D进行改进［6，7］和采用双点源近似［8，9］．2005年，Ripoll等［6］将一个适合于吸收的漫射系数D用在漫射近似中，并以实验证实改进的漫射近似能够描述光在强吸收介质中的传输．Hull和Foster［7］在研究P3近似漫反射时提出了一种P3近似与漫射近似的混合模型，该理论模型将漫射近似中的漫射系数D用描述P3近似的漫射系数Dasym替代．在研究漫反射问题时，双点源近似的空间分辨漫反射解被证明能较好地描述光源附近的辐射强度分布［8，9］．

漫射近似能对接近各向同性的光辐射做出精确的描述，而光源附近的光辐射分布是各向异性的．Bevilacqua等［10］用Monte Carlo方法对光源附近约1个输运平均自由程以内的漫反射进行了研究，指出光源附近组织的光辐射与相函数的前两阶矩有很大关系，并引入了与相函数二阶矩有关的光学参量γ．Kienle等人［11］的研究表明，用Monte Carlo (MC)模型反演组织的光学参量时，仅考虑g值是不够的，还应该考虑二阶参量γ，否则产生误差会高达100%．2008年Seo等人［12］研究了δ-P1近似下的空间分辨漫反射率问题时，通过引入含有高阶参量的反照率，对光源附近的光辐射分布给出了更好的数学描述．我们近些年的研究表明，二阶参量γ的影响主要是通过光源带来的，因此考虑光源近似的高阶项的必要的．

本文研究了生物组织的δ-P1近似漫反射光学模型，推导了含有等效光源一阶矩的双点源近似空间分辨漫反射解Rδ-P1(ρ)．与漫射近似RSDA(ρ)仅仅是μa和μ's的函数相比，Rδ-P1(ρ)还与参量γ有关，因此新的光学模型能较好地描述光源附近的漫辐射强度分布．

2.光辐射输运方程的δ-P1近似

稳态情况下，介质中r点^s方向的辐射强度L(r，^s)满足Boltzmann方程［1］

其中，μt=μa+μs为介质的衰减系数，μa和μs分别为介质的吸收系数和散射系数．p(^s，^s')为散射相函数，表示来自^s'方向的光被散射到^s方向的概率，通常采用的实际散射相函数是Henyey-Greenstein相函数，

其中g的取值从－1到1，称为各向异性因子．g=0表示各向同性散射，g=－1表示完全后向散射，g= 1表示完全前向散射．生物组织具有前向散射的特性，一些组织的g值达到0.99—0.999［13］．如果用delta函数δ(^s－^s')描述光子的前向散射，散射相函数可表示为

辐射强度通常包括尚未被散射的约化强度Lc(r，^s)(来自光源)和漫射强度Ld(r，^s)两部分，即L(r，^s)能够表示为

将(3)式和(4)式代入(1)式后得到

其中，μ*t=μa+μ*s，μ*s=μs(1－f)．式中Sc(r，^s)，是由约化入射强度Lc(r，^s)引起的等效的源函数，又称为有效光源项．

在δ-P1近似下，(5)式中Ld(r，^s)，p*(^s，^s')和Sc(r，^s)在其球谐函数展开式中只取前两项，

式中η=cosθ=^s·^r，η0=cosγ=^s·^s'．将(6)式代入(5)式得到φ0(r)和φ1(r)满足的方程组

在δ-P1近似下，相函数pδ(^s，^s')与pHG(^s，^s')应该满足二阶相似关系．pHG(^s，^s')的Legendre矩，因此g1=g， g2=g2．根据(2)式，pδ(^s，^s')的一阶和二阶Legendre矩可表示为f+(1－f)g1*和f+(1－f)g2*．根据相似关系［7］得到

根据(6b)式，有g2*=0，因此从(8b)式得到f=g2．利用f=g2，我们得到

根据(9)式，(7)式中μ*tr=μa+μs(1－g)=μtr，因此φ0(r)和φ1(r)满足的方程组最终可表示为

漫射近似假设了q1(r)=0，因此有

对于均匀无限介质，(11a)式格林函数解为

式中D=1/3μtr是漫射系数．辐射通量φ0(r)∝ΦG(r)⊗q0(r)，符号“⊗”表示卷积;辐射强度

3.δ-P1近似下的有效光源项及其空间分辨漫反射解

从上面分析我们看到，如果不考虑有效光源的一次项q1(r)，我们并不能在输运方程的δ-P1近似理论框架下解决任何新的问题．下面我们在(10)式基础上，考虑q1(r)≠0情况下φ0(r)和φ1(r)的解．当q1(r)≠0时，由(10 b)式得到

与漫射近似的(11b)式相比，式中多了q1(r)/μtr项．利用，得到

由约化入射强度Lc(r，^s)引起的等效的源函数Sc(r，^s)是复杂的，在实际应用中往往将Lc(r，^s)作点源近似．对无限细的准直光束沿z轴入射到一半无限均匀介质的情况，考虑双点源近似［9］

式中α*=μ*s/μ*t，表示反照率;z01=2/μ*t=2 z0， z02=0，z03=－2 zb，z04=－2zb－z01，如图1所示;外推长度zb=2 AD，其中A=(1+Reff)/(1－Reff)，有效反射系数Reff与界面的相对折射率nrel有关，可近似表示为［14］

在本文中，均采用nrel=1.4．根据(15)和(16)式，(13)式中q1(r)可以表示为

q1(r)=0时的空间分辨漫反射率为［7］

其中

根据前面分析，考虑光源一次项q1(r)的δ-P1空间分辨漫反射率可以表示为

4.结果与讨论

为评估所得到的双点源δ-P1近似模型，我们采用由Wang等［15］编写的MC模拟程序得到的反射率RMC(ρ)作为参考标准，将Rδ-P1(ρ)与漫射近似漫反射解RSDA(ρ)进行了比较．对于RSDA(ρ)的计算，(19)式中φ0(ρ)和j(ρ)为［16］

在图2和图3中，图示的横坐标均为光源与探测器之间的距离;图2中的纵坐标表示空间分辨漫反射率R(ρ)，图3中的纵坐标表示漫反射率R(ρ)与检验标准RMC(ρ)之间的误差ΔR/R=|R(ρ)－RMC(ρ)|/RMC(ρ)，R(ρ)分别表示Rδ-P1(ρ)和RSDA(ρ)．此外，进行MC模拟计算时，考虑准直光束入射到半无限均匀介质的情况，模拟的光子个数为106．

在图2中，我们比较了RMC(ρ)，RSDA(ρ)和Rδ-P1(ρ)三种反射率曲线，介质的参数取值为μa= 0.02 mm－1，μs=20.0 mm－1，g=0.9．在此参数取值下，该介质的输运反照率α'=μ's/μtr=0.99，输运平均自由程l't=1/μtr=0.495 mm．为了说明(20)式中光源项q1(ρ)/μtr对反射率的影响，曲线*Rδ-P1(ρ)是(20)式中不含q1(ρ)/μtr项的反射率曲线．从图中可以看出，与漫射近似模型相比，我们所得到的双点源δ-P1近似漫反射解在一个输运平均自由程之内(ρ＜0.4 mm)与MC模拟结果仍然保持很好的一致性，并且由于q1(ρ)/μtr项的存在，使其与MC模拟结果之间的误差减小．

保持μs=20.0 mm－1，g=0.9不变，图3考虑了吸收系数μa值分别为0.01 mm－1，0.02 mm-1，0.05 mm－1，0.1 mm－1，0.2 mm－1，0.5 mm－1六种情况，相应的a'分别为0.995，0.990，0.976，0.952，0.909，0.8，而l't数值上是a'一半，因此横坐标ρ＜4 mm的范围，覆盖了大约8—10 lt．由图3所示的误差曲线可以看出，在ρ≥1.5 lt范围内，RSDA(ρ)曲线的误差均大于10%，并且随着μa的增大，误差高达20%以上．而对于Rδ-P1(ρ)曲线，当a'＜0.952时，其误差均小于10%．

图4中显示了各向异性常数g取不同值时Rδ-P1(ρ)与Monte Carlo模拟计算结果RMC(ρ)之间的误差．图中曲线具有相同的μs=20.0 mm－1和μa=0.02 mm－1，g在0.7—0.99范围内取值(对于大多数生物组织，其g值在这个范围内)．结果显示在ρ≥0.5 lt范围内，误差均小于10%．

研究表明，双点源近似下的δ-P1近似理论模型的适用范围优于漫射近似理论模型，即能够应用于吸收系数较大的散射介质，并能较好地描述光源附近的辐射强度分布．此外，由于该模型考虑了有效光源项对近场辐射强度分布的影响，使得漫反射解的表达式中含有表征生物组织的光学特性的高阶光学参量γ，而漫射近似仅含有表征生物组织的光学特性的一阶光学参量μa和μ's．这项工作将为生物组织光学特性参量的在体无损检测提供重要的理论依据．

［1］Kim A D 2004 J．Opt．Soc．Am．A 21 820

［2］Wang Q H，Zhang Y Y，Lai J C，Li Z H，He A Z 2007 Acta Phys．Sin．56 1203(in Chinese)［王清华、张颖颖、来建成、张小娟、高宗慧、田会娟2007物理学报56 1203］

［3］Liu Y，Wang L J，Guo Y F，Zhang X J，Gao Z H，Tian H J 2007 Acta Phys．Sin．56 2119(in Chinese)［刘迎、王利军、郭云峰、张小娟、高宗慧、田会娟2007物理学报56 2119］

［4］Xu T，Zhang C P，Chen G Y，Zhang X J，gao Z H，Tian H J 2005 Chin．Phys．14 1813

［5］Xu L Q，Li H，Xie S S 2008 Acta Phys．Sin．57 6024(in Chinese)［徐兰青、李晖、谢树森2008物理学报57 6024］

［6］Ripoll J，Yessayan D，Zacharakis G，Ntziachristos V 2005 J．Opt．Soc．Am．A 22 546

［7］Hull E L，Foster T H 2001 J．Opt．Soc．Am．A 18 584

［8］Tian H J，Liu Y，Wang L J，Zhang Z B，Xiao L F 2009 Acta Phys．Sin．58 243(in Chinese)［田会娟、刘迎、王利军、张智卜、肖立峰2009物理学报58 243］

［9］Wang R，Liu Y 2010 Chinese J．Lasers．37 1147(in Chinese)［王锐、刘迎2010中国激光37 1147］

［10］Bevilacqua F，Depeursinge C 1999 J．Opt．Soc．Am．A 16 2935

［11］Kienle A，Forster F K，Hibst R 2001 Opt．Lett．26 1571

［12］InSeok Seo，Carole K．Hayakawa，Vasan Venugopalan 2008 Med．Phys．35 681

［13］Zhang Z X 2008 New Technology and Appliccations of Biomedical Photonics(Beijing:Science Press)(in Chinese)［张镇西2008生物医学光子学新技术及应用(北京:科学出版社)］

［14］Groenhuis R A J，Ferwerda H A，Bosch J J T 1983 I:Theory Appl Opt．22 2456

［15］Wang L H，Jacques S L，Zheng L Q 1995 Computer Methods and Programs in Biomedicine 47 131

［16］Kienle A，Patterson M S 1997 J．Opt．Soc．Am．A 14 246

PACS:42．25．Dd，42．62．Be，87．64．Cc

*Project supported by the National Natural Science Foundation of China(Grant No．60278004)．

E-mail:liuying@tju．edu．cn

δ-P1approximation model of biologicaltissues*

Liu Ying1)Liu Xiao-Jun2)Qi Bei-Bei1)Tian Hui-Juan3)
1)(College of Science，Tianjin University，Tianjin 300072，China)
2)(Key Laboratory of Opt-electronics Information Technical Science，EMC Tianjin University，Tianjin 300072，China)
3)(School of Information and Communication Engineering，Tianjin Polytechnic University，Tianjin 300160，China)
(Received 2 August 2010;revised manuscript received 2 September 2010)

A modified delta-P1approximation is studied in this paper．An analytic solution Rδ-P1(ρ)of this model for reflectance in double source approximation is developed by taking the source series expansion terms to the first order．It is shown that when the source in extended to first order，a better description for high absorption and high scattering media will be developed．The new model can give satisfactory results in describing irradiance distribution close to source compared with the RSDA(ρ)of diffusion approximation．Moreover，Rδ-P1(ρ)provides the prediction of the second-order parameterγ，which is significant for developing the inversion algorithm to recover optical parameters from spatially-resolved diffuse reflectance by using smaller source-detector separations．

tissue optics，δ-P1approximation，effective source，diffuse reflectance close to source

*国家自然科学基金(批准号:60278004)资助的课题．

E-mail:liuying@tju．edu．cn