CN105300928B - 一种大面积获取组织光学参数及微观结构的光反射成像技术 - Google Patents

Abstract

本发明提供一种大面积获取组织光学参数及微观结构的光反射成像技术，包括在任意空间距离和在整个空间频率域中均适用的浑浊介质反射率计算方法，以及一种通过测量高和低空间频率浑浊介质反射，并以得到的光反射率反演得出介质光学参数的方法。反演方法可以是运用查表法或者公式拟合法。浑浊介质sub‑扩散和漫射光反射率的测量可用于大面积测定混浊介质的光学特性及包括相位函数等的微观结构。

Description

一种大面积获取组织光学参数及微观结构的光反射成像技术

技术领域

本发明提供了一种大面积获取组织光学参数及微观结构的光反射成像技术，可用于生物医学光学和遥感等领域。

背景技术

光的弹性散射早已被用于分析随机介质。反射率光谱和成像是一种广泛使用的非侵入性方法，用来测量随机介质(如大气、海洋和组织)的光学特性，包括吸收系数(μ_a)和约化散射系数(μ_s′)。这些参数可以提供关于介质微结构和生化成分的有价值信息，并已应用于如云遥感、监测细胞凋亡、皮肤表征和癌症检测等领域。由于辐射传输(radiativetransfer-RT)描述光在随机介质传播，使得散射光反射率本质上是一个困难的问题。特别在光源-检测器短距离情况下，对RT通常采用的漫射近似不能工作。因此从反射率测量来量化介质的相位函数仍然难以实现的，特别是获取包含有关介质微环境的基本信息的介质相位函数。在任意光源-检测器距离情况下，反射率的精确分析模型是非常需要的。这种模型将适用于光学特性，尤其是随机介质的相位函数，的快速定量评估。最近提出的在光源-探测器近距离情况下的一些反射率经验模型都具有各自的局限性。由于在光源-检测器近距离情况下，散射介质的相位函数对sub-漫反射影响显著，所以一个明确关于sub-散射反射率与相位函数关系并且可以从反射率分布反推出随机介质光学参数(包括其相位函数)的分析模型是极其需要的。

发明内容

本发明涉及一种大面积获取组织光学参数及微观结构的光反射成像技术。其核心在于结合辐射传输的小角散射近似(SAA)给出近距离sub-漫散射光反射率与散射介质相位函数的定量解析关系。

本发明提供：一种大面积获取组织光学参数及微观结构的光反射成像技术，根据

分别获得低频的反射率公式

I_snake(q)+I_diffuse(q)+μ_b/2μ_t

以及高频的反射率公式

I_SAA(q)，

其中q_c～2πβ，

并将得到的反射率依次通过低频的反射率公式以反演方法获得吸收系数μ_a以及约化散射系数μ_s′，

通过高频的反射率公式以反演方法获得各向异性因子g、传播长度尺度l_Θ以及散射介质的任意相位函数的各向同性散射本底p_b。

其反演方法采用查表法或公式拟合法。

高频的反射率公式

由SAA扩散函数

获得。

假设p_Forward(θ)满足高斯分布，使得散射介质的任意相位函数p(θ)修改成为：

以及

通过以p_iso＝2p_b，将SAA扩散函数转化为

当ql_t＞＞1和ρ＜＜l_t时，

高频反射率可简化为

，其中μ_t′≡μ_a+2p_bμ_s和传播长度尺度l_Θ≡Θ/μ_s。

利用高频公式I_SAA(q)和I_SAA(ρ)，从两个或多个介质的高频反射图谱，直接获得表征体系微观结构的参数或这些参数的组合。

低频的反射率公式

和

由蛇形光子和弥漫性光子的表达式

获得。

利用公式I_{snake，diffuse}(q)或I_{snake，diffuse}(ρ)，从一个或多个介质的低频反射图谱，直接获得其吸收系数及约化散射系数。

另外，可首先利用sub-漫散射光反射率进行拟合得到μ_b/μ_s，μ_a/μ_s和Θ/μ_s的值，确定μ_b/μ_s后设

g＝(1-2p_b)(1-Θ²/2)

，然后利用最小二乘法拟合对sub-漫射光和漫射光反射率分布进行拟合，得到所有光学参数。

本发明包括在任意空间距离和在整个空间频率域中均适用的浑浊介质反射率计算方法，以及一种通过测量高和低空间频率浑浊介质反射，并以得到的光反射率反演得出介质光学参数的方法。反演方法可以是运用查表法或者公式拟合法等。浑浊介质sub-扩散和漫射光反射率的测量可用于大面积测定混浊介质的光学特性及包括相位函数等的微观结构。散射介质的相位函数携带了关于形态学的基本信息以及单个散射体的光学特性。另外相位函数的分析可以用来预测光传播以及检测随机介质内微小结构变化或不均匀性。本发明所提出的分析模型在快速量化散射介质的所有光学特性，包括相位函数等，特别是在生物医学光学和遥感等领域，将获得广泛而重要的应用。

附图说明

图1为本发明的SAA光子的散射示意图。

图2为Monte Carlo模拟直径为1.5mm的聚苯乙烯球悬浮液在水中的光反射率对比示意图。

图3为利用SAA相位函数来近似精确米氏相位函数的示意图1。

图4为利用SAA相位函数来近似精确米氏相位函数的示意图2。

图5为利用SAA相位函数来近似精确米氏相位函数的示意图3。

图6为在原点垂直入射的平行光束，离原点径向位置r的反射率的空间分布示意图1。

图7为在原点垂直入射的平行光束，离原点径向位置r的反射率的空间分布示意图2。

图8为在原点垂直入射的平行光束，离原点径向位置r的反射率的空间分布示意图3。

图9为空间频率为q的垂直入射空间调制平面波的反射率的示意图1。

图10为空间频率为q的垂直入射空间调制平面波的反射率的示意图2。

图11为空间频率为q的垂直入射空间调制平面波的反射率的示意图3。

图12为通过公式I_SAA(q)对图2中显示的聚苯乙烯悬浮液进行拟合的结果示意图。

具体实施方式

下面结合附图对本发明实施例作进一步说明：

将一准直光束在z＝0的界面沿方向入射，通过前向峰散射介质，后向散射的方向为的光反射率。在这样的介质中，非弥漫光子只会遇到少数大角度散射，非弥漫光子根据大角度散射次数n分为n阶非弥漫光子。由于一阶非弥漫光子经历了多个小角散射和一个大角度散射，主要贡献可以描述为使用小角度散射近似到辐射传输的近距离形式。一阶非弥漫光子被称为SAA光子。二阶非弥漫光子被称为蛇形光子(“snake”光子)。

SAA光子后向散射主要是由向前方向上散射角的扩散和后向散射效率来决定。假定散射介质的相位函数p(θ)(归一化为)分为一个前向峰散射成分和各向同性散射成分，如

P_SAA(θ)＝(1-2p_b)P_Forward(θ)+(2π)^-1p_b(p_b＜＜1)。

图1.(a)经过多次小角散射(“红”散射体)和单各大角度散射(“蓝”散射体)的一阶非弥漫光子的后向散射。在第一次或最后一次散射事件，光子可以分别以概率为p_iso和1-p_iso发生各向同向散射(虚线)与前向峰散射(实线)。(b)S^eff′的列入提高了SAA光子的精度，显示了Monte Carlo模拟直径为1.5mm的聚苯乙烯球悬浮液在水中的光反射率对比。

而在分层介质中，一准直光束在原点(r＝0)以s₀方向入射到达深度z的SAA扩散函数由下式给出：

其中q是在xy平面上的空间频率，S_⊥0是S0在分界面上的投影，μ_t≡μ_s+μ_a，其中μ_s是散射系数，χ(v，z)是(1-2p_b)p_Forward(θ，z)的2维傅里叶变化。通过扩展第二指数项成具有0阶弹道项的

，可以识别出包含来自各自所有散射阶次0、1、…的S。

然后从中获得SAA光子的反射率，可以表示为

其中，后向散射系数μ_b≡μ_s(z)p(π，z)，s_⊥＝s_⊥in+s_⊥out，S^eff是对具有相同相位函数和两倍吸收与散射的有效介质的扩散函数，S^eff′是对相同相位函数和散射而吸收系数被修改为2μ_a+2p_bμ_s的第二有效介质的扩散函数。在上述公式(2)第二项说明当光子进入或逃逸介质时，进行各向同性散射路径而不是前向峰散射路径，他们的概率分别为p_iso和1-p_iso(p_iso＜＜1)(见图1)。第二项的列入提高了SAA的准确性(见图2)。注意在S^eff′中应去除弹道项，以避免重复计算。

假设p_Forward(θ)满足高斯分布，使得散射介质的任意相位函数p(θ)修改成为：

由中得出

以及

P_SAA的n阶矩由得到。具体地，各向异性因子(n＝1)是

g＝(1-2p_b)[1-Θ²/2]。 (4)

通过将公式(2)中的p_iso＝2p_b，将SAA扩散函数即公式(1)简化为

其中erf是误差函数。相比其它辐射传输的近似解截断相位函数矩的顺序，SAA光子的反射率(2)包含来自相位函数所有矩的贡献。

P_SAA的表达与由组织和细胞引起的光散射的统一米氏和分形模型是一致的。高斯项捕获来自米氏散射的贡献，而各向同性散射项与背景的折射率波动相关联。均方根散射角Θ随着组织和细胞中的米氏散射体(大结构)尺寸减小。

当ql_t＞＞1和ρ＜＜l_t时，

高频反射率可简化为

，其中μ_t′≡μ_a+2p_bμ_s和传播长度尺度l_Θ≡Θ/μ_s。

利用高频公式I_SAA(q)和I_SAA(ρ)，从两个或多个介质的高频反射图谱，可以直接反演获得表征体系微观结构的参数或这些参数的组合。

而低频的反射率公式

和

由蛇形光子和弥漫性光子的表达式

获得。

其中β≡μ_a+μ_s′，μ_s′≡μ_s(1-g)，g为各向异性因子，

G^{(snake，diffuse)}分别是蛇形光子和弥漫性光子的格林函数，

G^(snake)(r，r′)＝exp(-β|r-r′|)/4π|r-r′|²为蛇形光子格林函数，在各向同性散射混浊介质内是各向同性源弹道传播。在傅立叶域，蛇形光子和弥漫性光子的反射率可以简化为：

和

其中z_e是依赖于界面处折射率不匹配的外推长度。

利用公式I_{snake，diffuse}(q)或I_{snake，diffuse}(ρ)，从一个或多个介质的低频反射图谱，可以直接反演获得其吸收系数及约化散射系数。

图3、图4和图5比较了由蒙特卡罗模拟空间域和傅里叶域下水中直径为1.5μm(没有吸收和有吸收，μ_a/μ_s′＝0.16)和0.49μm(没有吸收)聚苯乙烯球悬浮液得到的结果与SAA光子、蛇形光子和弥漫性光子光反射率。

在蒙特卡罗模拟中，入射总光子数设为106，假定半无限介质和周围的折射率是匹配的。通过对聚苯乙烯粒子的折射率添加非零虚部实现散射介质的吸收。相对应于SAA光子的参数分别是p_b＝0.0177和Θ＝0.451，p_b＝0.0171和Θ＝0.447，p_b＝0.0194和Θ＝0.587。

图6、图7、图8比较SAA光子、蛇形和弥漫性光子反射率与蒙特卡罗模拟水中直径为1.5μm(图6：没有吸收，图7：μ_a/μ_s′＝0.16；g＝0.92)和0.49μm(图8：没有吸收，g＝0.86)聚苯乙烯球悬浮液得到的结果。入射光波长是0.515μm。

图3-图5显示了利用SAA相位函数来近似精确米氏相位函数。图6-图8显示了在原点垂直入射的平行光束，离原点径向位置r的反射率的空间分布。图9、图10和图11显示了空间频率为q的垂直入射空间调制平面波的反射率。

表1.聚苯乙烯球悬浮液模型(直径为d，吸收系数为μ_a＝0)和蒙特卡罗模拟之间匹配关系。R²＝1表示模型与真实值相近程度。

表2，对比蒙特卡罗模拟聚苯乙烯球体悬浮液(d＝1.5μm，μ_a/μ_s′＝0.16)光反射率拟合的参数与理论值。

SAA光子和合并蛇形与弥散射的光子分别很好地描述实际空间中光源与检测器在短距离和长距离的反射率(图6-图8)，以及频率域中高和低空间频率的反射率(图9-图11)。

当距离大于β^-1时，对于后向散射光这种划分与独立于相位函数具体形式的通用径向分布观测相一致。公式(8)完美描述了具有低中度吸收(μ_a/μ_s′＜1)的前向峰散射介质的光反射率。表1显示了在傅里叶空间模型和蒙特卡罗模拟计算在0≤q/μ_s′≤100内的结果相近程度。模型与蒙特卡罗模拟之间的相近程度定义如下

其中SSR＝∑(logI_model-logI_MC)²是模型的均方误差，

是真实值偏差，是logI_MC的平均值。R²意味着模型与真实值完美重合匹配。相近程度值是在傅里叶空间0≤q/μ_s′≤100范围内计算的。当吸收增大(μ_a/μ_s′＞≈1)时，光子在其传播方向不再完全随机化或在介质中光散射不再前向峰散射，公式(8)精度降低。

光反射率的封闭形式：

光反射率的封闭形式(8)和极限形式(6)，(7)表明，除了众所周知的在低空间频率(远距离)作用于光漫射和反射率的传输平均自由程l_t＝1/(1-g)μ_s外，有一新的表征高空间频率(近距离)sub-漫射光反射率的传播长度尺度l_Θ≡Θ/μ_s，它仅由相位函数前向峰部分的角度宽度决定。对于前向峰散射介质，上述两尺度可以通过公式(4)进一步建立内在联系，近似满足μ_s ³l_tl_Θ ²＝2。通过公式(6，7)可以方便测量高频(近距离)光反射率。图12显示了通过公式(6)对图3-图11中显示的聚苯乙烯悬浮液进行拟合的结果。散射光的扩散角度Θ分别由2.4％、3.6％和5.6％的精度决定。

通过拟合出低和高空间频率的光反射率确定混浊介质的完整一组光学特性，该组光学特性包括完全由散射光扩散角度Θ定义的SAA相位函数、散射介质的任意相位函数的各向同性散射本底p_b、后向散射系数μ_b、散射系数μ_s、吸收系数μ_a和各项异性因子g。

举一例子，表2显示了对蒙特卡罗模拟的聚苯乙烯悬浮液(d＝1.50μm，μ_a/μ_s′＝0.16)的光反射率进行拟合，参数结果如表2所示。

另外，可以首先利用sub-漫散射光反射率进行拟合得到μ_b/μ_s，μ_a/μ_s和Θ/μ_s的值，确定μ_b/μ_s后设g＝(1-2p_b)(1-Θ²/2)，然后利用最小二乘法拟合对sub-漫射光和漫射光反射率分布进行拟合，得到所有光学参数。

由光反射率得到的这组光学参数与理论值很好地吻合。尤其，所提取的散射特性和相位函数参数(μ_s，g，Θ)精度是非常好的，这部分原因是由g和Θ之间的制约因素引起的。g和Θ分别表征由前向峰散射介质引起的漫散射和sub-漫散射光反射率。

总之，本发明提出了一种描述在任意光源与探测器距离下前向峰散射介质的反射率分析模型，及其在分析前向峰散射介质的光学参数及相位函数的应用。当吸收处于中低度水平时，该分析模型在整个空间尺度显示出优异的性能，并已成功证明了应用该模型从介质sub-扩散和漫射光反射率精确测定混浊介质的光学特性和相位函数。散射介质的相位函数携带了关于形态学的基本信息及其单个散射体的光学特性。另外相位函数的分析也可以用来预测光传播以及检测随机介质内微小结构变化或不均匀性。随着近年来如空间调制照明的发展，快速测量宽视场中sub-散射和漫散射光的反射率已成为可能。因此，本发明所提出的分析模型在快速量化散射介质的所有光学特性，包括相位函数等，特别是在生物医学光学和遥感等领域，将获得广泛而重要的应用。

Claims (9)

1.一种大面积获取组织光学参数及微观结构的光反射成像方法，

其特征在于：根据

分别获得低频的反射率公式

I_snake(q)+I_diffuse(q)+μ_b/2μ_t

以及高频的反射率公式

I_SAA(q)，

其中q_c＝2πβ,

并将得到的反射率依次通过低频的反射率公式以反演方法获得吸收系数μ_a以及约化散射系数μ_s'，μ_b/2μ_t来源于SAA光子的弹道贡献，q_c为两个条件的交点，μ_t≡μ_s+μ_a，μ_s是散射系数，

通过高频的反射率公式以反演方法获得各向异性因子g、传播长度尺度l_Θ以及散射介质的任意相位函数的各向同性散射本底p_b，

I(q)为光反射率表达式，I_snake(q)为snake光子的光反射率表达式，I_SAA(q)为SAA光子的光反射率表达式，I_diffuse(q)为diffuse光子的光反射率表达式；q是在xy平面上的空间频率，μ_b是后向散射系数，β≡μ_a+μ_s'，μ_s'≡μ_s(1-g)，z是深度，S^eff是对具有相同相位函数和两倍吸收与散射的有效介质的扩散函数，s_⊥是光传播方向s在分界面上的投影，p_iso是进入介质时进行各向同性散射路径的概率，S^eff'是对相同相位函数和散射而吸收系数被修改为2μ_a+2p_bμ_s的第二有效介质的扩散函数；z_e是依赖于界面处折射率不匹配的外推长度，q_c为两个条件的交点。

2.根据权利要求1所述的一种大面积获取组织光学参数及微观结构的光反射成像方法，其特征在于：其反演方法采用查表法或公式拟合法。

3.根据权利要求1所述的一种大面积获取组织光学参数及微观结构的光反射成像方法，其特征在于：高频的反射率公式

由SAA扩散函数

获得，ζ是积分变量，χ是傅里叶空间的光散射相位函数，s_⊥0是是光传播方向s₀在分界面上的投影。

4.根据权利要求3所述的一种大面积获取组织光学参数及微观结构的光反射成像方法，其特征在于：假设p_Forward(θ)满足高斯分布，使得散射介质的任意相位函数p(θ)修改成为：

由中得出

以及

通过以p_iso＝2p_b，将SAA扩散函数转化为

p_SAA(θ)是SAA光子的相位函数，θ为散射角度，Θ为均方根散射角，

p(θ)为散射介质的任意相位函数，

erf为误差函数。

5.根据权利要求3所述的一种大面积获取组织光学参数及微观结构的光反射成像方法，其特征在于：当ql_t＞＞1和ρ＜＜l_t时，高频反射率可简化为

和

，

其中μ_t'≡μ_a+2p_bμ_s和传播长度尺度l_Θ≡Θ/μ_s，Γ为伽马函数。

6.根据权利要求5所述的一种大面积获取组织光学参数及微观结构的光反射成像方法，其特征在于：利用高频公式I_SAA(q)和I_SAA(ρ)，从两个或多个介质的高频反射图谱，直接获得表征体系微观结构的参数或这些参数的组合。

7.根据权利要求1所述的一种大面积获取组织光学参数及微观结构的光反射成像方法，其特征在于：低频的反射率公式

和

由蛇形光子和弥漫性光子的表达式

获得，r为原点径向位置。

8.根据权利要求7所述的一种大面积获取组织光学参数及微观结构的光反射成像方法，其特征在于：利用公式I_{snake,diffuse}(q)或I_{snake,diffuse}(ρ)，从一个或多个介质的低频反射图谱，直接获得其吸收系数及约化散射系数。

9.根据权利要求1所述的一种大面积获取组织光学参数及微观结构的光反射成像方法，其特征在于：另外，可首先利用子表面(sub)-漫散射光反射率进行拟合得到μ_b/μ_s，μ_a/μ_s和Θ/μ_s的值，确定μ_b/μ_s后设

g＝(1-2p_b)(1-Θ²/2)，

然后利用最小二乘法拟合对子表面(sub)-漫射光和漫射光反射率分布进行拟合，得到所有光学参数，Θ为均方根散射角。