CN106290199A - 基于稳态辐射率间质测量装置的多波长光学参数反构方法 - Google Patents

Abstract

本发明公开了一种基于稳态辐射率间质测量装置的多波长光学参数反构方法，包括：设定光源‑探测器的距离，以1°为间隔，测量水平面360°范围内辐射率的分布；选取波长，并从测量信号中，提取出该波长下的辐射率分布，利用本发明所选择的辐射传输方程在P9近似下的解析解作为光子传输正向模型拟合实验所测量的辐射率，最终重建该波长下待测样品的光学参数。本发明仅需在单个光源‑探测器距离下进行稳态辐射率间质测量，有望实现微创检测；本发明能对低约化反照率样本(如前列腺增生组织)的光学参数进行准确反构；本发明所提供方法应用前景广阔，除生物体组织方面的应用，可用于牛奶、液态仿体等光学参数的测量。

Description

基于稳态辐射率间质测量装置的多波长光学参数反构方法

技术领域

本发明属于组织光学研究中的光学参数检测领域，具体涉及一种基于稳态辐射率间质测量装置的多波长组织光学参数反构方法。

背景技术

前列腺癌是男性常见的恶性肿瘤之一，由于生活方式的改变及人口老龄化，其发病率及死亡率显著上升，已成为影响中国男性健康的重要疾病之一。近年来，光动力疗法(Photodynamic Therapy，PDT)成为前列腺治疗的新手段^[1]。在PDT治疗中，PDT的光剂量是制定治疗方案及在线优化的关键因素。光剂量是指通过测量组织吸收系数、散射系数的分布决定光通量密度的分布，由于从光学参数计算光通量密度分布的方法已较成熟，因此光剂量的关键就是获得组织的光学参数的分布。此外，由于前列腺是内部器官，因此本发明考虑的是间质光学参数测量。

从光源方面着眼，间质光学参数的测量可分为使用连续光照射的稳态方法、使用光强度调制的频域法及使用脉冲光源的时间分辨方法，其中，稳态测量系统简单、快捷、廉价，有望应用与临床医学。从测量物理量着眼，间质光学参数测量方法可分为基于光子密度的测量方法及基于辐射率的测量方法。常用的基于光子密度的测量方法需要测量多个光源-探测器距离下的光子密度分布，将会加剧内部器官的损伤^[1]。与此相对应的，基于辐射率的间质光学测量参数方法^[2]仅需要在固定的光源-探测器距离下在测量多个角度下的辐射率分布，有望实现无创或微创测量。此外，与光子密度测量方式相比，辐射率含有的角度信息对光学参数的变化更加敏感，且具有较强的抗噪能力^[3]。因此，本发明选择基于稳态辐射率的间质光学参数测量。

光学参数反构的通常思路为，依据辐射率和生物组织光学参数之间的正向模型，将不同角度下测量的辐射率代入反构算法得到生物组织光学参数的估计值，因此准确的光子传输模型显得尤为重要。在辐射率测量基本框架下，常用光子传输正向模型包括，Monte-Carlo模拟模型、辐射传输方程(radiative transfer equation，RTE)模型、扩散近似模型、PN近似模型。其中，Monte-Carlo模型能够准确地模拟光在组织中的传输过程，但其由于其耗时过长，不适用于临床中快速检测的需求；RTE被认为是描述光子在生物组织中传输的准确模型，但其解析解很难直接得到；扩散近似模型最为简单，但在间质辐射率测量方式中，当光源-探测器距离较小，组织(如前列腺增生组织)的约化返照率较低时，不满足该模型的要求^[4]；PN近似模型为RTE在N阶球谐函数展开下的解析解^[5]，其解具体为无限媒质中的各项同性点光源下的辐射率。经计算知，PN近似解的误差随阶数N的增大而减小，但当N大于9时，误差几乎不随N的增大而减小，且误差很小。与此同时，PN近似解的计算时间随阶数N的增大而增加。综合考虑，本发明选择辐射传输方程的P9近似下的解析解作为光子传输正向模型。

[参考文献]

[1]T.C.Zhu,J.C.Finlay,S.M.Hahn,"Determination of the distribution oflight,optical properties,drug concentration,and tissue oxygenation in-vivo inhuman prostate during motexafin lutetium-mediated photodynamic therapy,"Journal of Photochemistry and Photobiology B:Biology,79(3):231-241(2005)。

[2]S.Grabtchak,L.G.Montgomery,W.M.Whelan,"Feasibility of interstitialnear-infrared radiance spectroscopy platform for ex vivo canine prostatestudies:optical properties extraction,hemoglobin and water concentration,andgold nanoparticles detection,"Journal of Biomedical Optics,19(5):057003(2014)。

[3]O.Barajas,A.M.Ballangrud,G.G.Miller et al,"Monte Carlo modeling ofangular radiance in tissue phantoms and human prostate:PDT light dosimetry,"Physics in Medicine and Biology,42(9):1675–1687(1997)。

[4]E.L.HULL,T.H.FOSTER,"Steady-state reflectance spectroscopy in theP3approximation,"Journal of the Optical Society of America A,18(3):584-599(2001)。

[5]A.Liemert,A.Kienle,"Analytical Green’s function of the radiativetransfer radiance for the infinite medium,"Physical Review,83(3):036605(2011)。

发明内容

针对多波长下组织体光学参数间质测量的需求，本发明提出一种基于稳态辐射率间质测量装置的多波长组织光学参数反构方法。

为了解决上述技术问题，本发明提出的一种基于稳态辐射率间质测量装置的多波长光学参数反构方法，其中，所采用的稳态辐射率间质测量装置的结构是，包括光源模块、探测模块、位移模块及控制模块；所述光源模块包括连接至光源光纤一端的卤钨灯；所述光源光纤的另一端设有各向同性散射头，用于将所述卤钨灯发出的光转为点光源；所述探测模块包括侧照光纤和光谱仪，所述侧照光纤的一端耦合至所述光谱仪，所述侧照光纤的另一端为与光纤轴线呈45°角的探测斜面，用于探测与光纤轴线垂直的入射的光信号；所述位移模块包括手动精密位移台、电动平移台和电动旋转台，所述电动平移台的上方固定有一金属立柱，所述手动精密位移台固定在所述金属立柱的上方，所述光源光纤通过一金属固定管与所述手动精密位移台固定；所述电动旋转台固定在一金属支架上，所述侧照光纤通过一金属固定管垂直的固定在所述电动旋转台上；所述手动精密位移台用于调整光源光纤的高度，使各向同性散射头与探测斜面保持在同一水平高度；所述电动平移台用于在水平方向移动光源光纤；所述电动旋转台用于旋转侧照光纤，确定探测斜面的朝向以改变光子接收方向；所述控制模块包括步进电机控制箱和上位计算机；所述上位计算机一方面发出指令至步进电机控制箱，从而控制电动平移台及电动旋转台的运动；另一方面，用于控制光谱仪探测光信号；

该多波长光学参数反构方法包括以下步骤：

步骤一、利用手动精密位移台、电动平移台，调整光源光纤的空间位置，从而确定各向同性散射头与探测斜面之间的距离，即光源-探测器之间的距离r；通过电动旋转台调整探测斜面的初始探测方向；设定电动旋转台的单位旋转角度为1°；

步骤二、所述卤钨灯发出白光，经过光源光纤传输至各向同性散射头，光经过各向同性散射后，以点光源的形式传输至置于样品池中的待测样品中；光经过待测样品的吸收散射后，经过探测斜面传入侧照光纤，被光谱仪接收；光谱仪测量到的信号传入上位计算机，并存储；

步骤三、以步骤一设定的单位旋转角度旋转电动旋转台，重复步骤二，直至实现在水平面内360°范围内辐射率的测量；

步骤四、波长λ下组织光学参数的反构，步骤如下：

首先，选取待反构光学参数所在波长为λ，提取出波长λ下的辐射率M＝[m₀,m₁,...,m_k,...,m₃₅₉]，其中，m_k为第k个角度下的辐射率测量值；

然后，设定最优化目标函数为：

$f({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)=\sqrt{{\mathrm{\Σ}}_{k=0}^{359}{(\frac{m_k}{m_0}-\frac{l_k({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)}{l_0({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)})}^2}---\left(1\right)$

式(1)中；分别为样品池中待测样品的吸收系数、散射系数的估计值；是光源-探测器距离为r、探测角度为k/180×π处P9近似下辐射传输方程的解，即P9近似下的辐射率，其表达式为：

$l_k({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)=\underset{n=0}{\overset{N}{\Σ}}\frac{2n+1}{4\mathrm{\π}}{\mathrm{\φ}}_n({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s,r)P_n\left(cos\right(\frac{k\·\mathrm{\π}}{180}\left)\right)---\left(2\right)$

式(2)中：n＝0,1,...,N，N＝9为近似阶数，k＝0,1,...,359为第k个探测角度；P_n为n阶勒让德多项式；φ_n(μ_a,μ_s,r)为第n阶光子密度计算值，其表达式为：

${\mathrm{\φ}}_n({\mathrm{\μ}}_a,{\mathrm{\μ}}_s,r)=\frac{1}{4\mathrm{\π}}\underset{i=1}{\overset{\frac{N+1}{2}}{\Σ}}{\mathrm{\υ}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)G_i\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)l_n\right({\mathrm{\υ}}_i\left({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\right)){\mathrm{\κ}}_n\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)r\right)---\left(3\right)$

其中，κ_n(x)的表达式如下：

$\left\{\begin{array}{c}{\mathrm{\κ}}_0\left(x\right)=\frac{\exp \left(x\right)}{x}\\ {\mathrm{\κ}}_1\left(x\right)=(1-\frac{1}{x})\frac{\exp \left(x\right)}{x}\\ {\mathrm{\κ}}_n\left(x\right)={\mathrm{\κ}}_{n-2}\left(x\right)-\frac{2n-1}{x}{\mathrm{\κ}}_{n-1}\left(x\right)\end{array}\right.---\left(4\right)$

定义递推式：B_n+1(ω(μ_a,μ_s))＝(2n+1)σ_nB_n(ω(μ_a,μ_s))+ω(μ_a,μ_s)n²B_n-1(ω(μ_a,μ_s))，其中，σ_n≡μ_a+μ_s(1-gⁿ)，g为组织体的各向异性因子，g取值为0.9；

当，B₁(ω(μ_a,μ_s))＝1，B₂(ω(μ_a,μ_s))＝3σ₁，时

$B_{N+1}\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=R\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\underset{l=0}{\overset{\frac{N-1}{2}}{\Σ}}{a}_l{\mathrm{\ω}}^l({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)---\left(5\right)$

当，B₀(ω(μ_a,μ_s))＝1，B₁(ω(μ_a,μ_s))＝μ_a，时

$B_{N+1}\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=Q\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\underset{l=0}{\overset{\frac{N+1}{2}}{\Σ}}{b}_l{\mathrm{\ω}}^l({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)---\left(6\right)$

ω_i(μ_a,μ_s)为Q(ω(μ_a,μ_s))＝0的根，

${\mathrm{\υ}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)=\sqrt{-{\mathrm{\ω}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)}---\left(7\right)$

$G_i\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\frac{1}{\frac{b_{N+1}}{2}}\frac{R\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)}{{\mathrm{\Π}}_{n=1,n\≠i}^{\frac{N-1}{2}}\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)-{\mathrm{\ω}}_n({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)}---\left(8\right)$

$\left\{\begin{array}{c}{l}_N\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\frac{v_i^{N}N!}{R\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s))}\\ l_{N+1}\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=0\\ l_{n-1}\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\frac{1}{n}\[\frac{2n+1}{{\mathrm{\υ}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)}{\mathrm{\σ}}_n{l}_n\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)-(n+1)l_{n+1}\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)\]\end{array}\right.---\left(9\right)$

最后，通过Powell方向加速法计算式(1)的最小值，与该最小值对应处的光学参数即为波长λ下组织光学参数的反构值。

与现有技术相比，本发明的有益效果是：

1.本发明所述光学参数反构方法基于双光纤稳态测量装置，所述装置与传统的间质测量装置相比，以旋转测量模式代替移动测量模式，减少了对组织体的损伤程度，可实现对组织体进行无创或微创的间质辐射率测量，从而用于光学参数的在体获取，此外，所述装置设备简单，造价低；

2.本发明所述光学参数反构方法基于辐射率测量模式，所述测量模式有较强的抗噪能力，且对光学参数更加敏感，有利于对光学参数进行更准确的反构；

3.本发明所述光学参数反构方法使用的正向模型为RTE在P9近似下的解，此正向模型在小光源-探测器距离、低约化返照率下具备更高的准确性，能对低约化反照率样本(如前列腺增生组织)的光学参数进行准确反构；

4.本发明提出的方法可以对多波长下的组织体的光学参数进行测量；

5.本发明所提供方法应用前景广阔，除生物体组织方面的应用，可用于牛奶、液态仿体等光学参数的测量。

附图说明

图1为本发明中采用的稳态辐射率间质测量装置原理图。

具体实施方式

下面结合附图和具体实施例对本发明技术方案作进一步详细描述，所描述的具体实施例仅对本发明进行解释说明，并不用以限制本发明。

如图1所示，本发明中作采用的稳态辐射率间质测量装置的结构是：包括光源模块、探测模块、位移模块及控制模块。

所述光源模块包括连接至光源光纤2一端的卤钨灯1；所述光源光纤2的另一端设有各向同性散射头3，用于将所述卤钨灯1发出的光转为点光源，所述光源光纤2的下侧包覆有金属管，便于固定。

所述探测模块包括侧照光纤5和光谱仪6，所述侧照光纤5的一端耦合至所述光谱仪6，所述侧照光纤5的另一端为与光纤轴线呈45°角的探测斜面4，用于探测与光纤轴线垂直的入射的光信号，所述侧照光纤5的下侧包覆有金属管，便于固定。

所述位移模块包括手动精密位移台7、电动平移台8和电动旋转台9，所述电动平移台8的上方固定有一金属立柱，所述手动精密位移台7固定在所述金属立柱的上方，所述光源光纤2通过一金属固定管与所述手动精密位移台7固定；所述电动旋转台9固定在一金属支架上，所述侧照光纤5通过一金属固定管垂直的固定在所述电动旋转台9上；所述手动精密位移台7用于调整光源光纤2的高度，使各向同性散射头3与探测斜面4保持在同一水平高度；所述电动平移台8用于在水平方向移动光源光纤2，以实现特定的光源-探测器距离(各向同性散射头3与探测斜面4之间的距离)；所述电动旋转台9用于旋转侧照光纤5，使探测斜面4改变光子接收方向。

所述控制模块包括步进电机控制箱10和上位计算机11；所述上位计算机11一方面发出指令至步进电机控制箱10，从而控制电动平移台8及电动旋转台9的运动；另一方面，用于控制光谱仪6探测光信号；

使用上述稳态辐射率间质测量装置实现多波长光学参数反构方法，步骤如下：

步骤一、利用手动精密位移台7、电动平移台8，调整光源光纤2的空间位置，从而确定各向同性散射头3与探测斜面4之间的距离，即光源-探测器之间的距离r；通过电动旋转台9调整侧照光纤5以得到探测斜面4的初始探测方向；设定电动旋转台9的单位旋转角度为1°；

步骤二、所述卤钨灯1发出白光，经过光源光纤2传输至各向同性散射头3，光经过各向同性散射后，以点光源的形式传输至置于样品池12中的待测样品中；光经过待测样品的吸收散射后，经过探测斜面4传入侧照光纤5，被光谱仪6接收；光谱仪6测量到的信号传入上位计算机11，并存储；

步骤三、以步骤一设定的单位旋转角度旋转电动旋转台9，重复步骤二，直至实现在水平面内360°范围内辐射率的测量；

步骤四、波长λ下组织光学参数的反构，步骤如下：

首先，选取待反构光学参数所在波长为λ，提取出波长λ下的辐射率M＝[m₀,m₁,...,m_k,...,m₃₅₉]，其中，m_k为第k个角度下的辐射率测量值；

然后，设定最优化目标函数为：

$f({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)=\sqrt{{\mathrm{\Σ}}_{k=0}^{359}{(\frac{m_k}{m_0}-\frac{l_k({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)}{l_0({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)})}^2}---\left(1\right)$

式(1)中；分别为样品池中待测样品的吸收系数、散射系数的估计值；是光源-探测器距离为r、探测角度为k/180×π处P9近似下辐射传输方程的解，即P9近似下的辐射率，其表达式为：

$l_k({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)=\underset{n=0}{\overset{N}{\Σ}}\frac{2n+1}{4\mathrm{\π}}{\mathrm{\φ}}_n({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s,r)P_n\left(cos\right(\frac{k\·\mathrm{\π}}{180}\left)\right)---\left(2\right)$

式(2)中：n＝0,1,...,N，N＝9为近似阶数，k＝0,1,...,359为第k个探测角度；P_n为n阶勒让德多项式；φ_n(μ_a,μ_s,r)为第n阶光子密度计算值，其表达式为：

${\mathrm{\φ}}_n({\mathrm{\μ}}_a,{\mathrm{\μ}}_s,r)=\frac{1}{4\mathrm{\π}}\underset{i=1}{\overset{\frac{N+1}{2}}{\Σ}}{\mathrm{\υ}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)G_i\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)l_n\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right){\mathrm{\κ}}_n\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)r\right)---\left(3\right)$

其中，κ_n(x)的表达式如下：

$\left\{\begin{array}{c}{\mathrm{\κ}}_0\left(x\right)=\frac{\exp \left(x\right)}{x}\\ {\mathrm{\κ}}_1\left(x\right)=(1-\frac{1}{x})\frac{\exp \left(x\right)}{x}\\ {\mathrm{\κ}}_n\left(x\right)={\mathrm{\κ}}_{n-2}\left(x\right)-\frac{2n-1}{x}{\mathrm{\κ}}_{n-1}\left(x\right)\end{array}\right.---\left(4\right)$

定义递推式：B_n+1(ω(μ_a,μ_s))＝(2n+1)σ_nB_n(ω(μ_a,μ_s))+ω(μ_a,μ_s)n²B_n-1(ω(μ_a,μ_s))，其中，σ_n≡μ_a+μ_s(1-gⁿ)，g为组织体的各向异性因子，g取值为0.9；

当，B₁(ω(μ_a,μ_s))＝1，B₂(ω(μ_a,μ_s))＝3σ₁，时

$B_{N+1}\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=R\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\underset{l=0}{\overset{\frac{N-1}{2}}{\Σ}}{a}_l{\mathrm{\ω}}^l({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)---\left(5\right)$

当，B₀(ω(μ_a,μ_s))＝1，B₁(ω(μ_a,μ_s))＝μ_a，时

$B_{N+1}\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=Q\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\underset{l=0}{\overset{\frac{N+1}{2}}{\Σ}}{b}_l{\mathrm{\ω}}^l({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)---\left(6\right)$

ω_i(μ_a,μ_s)为Q(ω(μ_a,μ_s))＝0的根，

${\mathrm{\υ}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)=\sqrt{-{\mathrm{\ω}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)}---\left(7\right)$

$G_i\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\frac{1}{\frac{b_{N+1}}{2}}\frac{R\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)}{{\mathrm{\Π}}_{n=1,n\≠i}^{\frac{N-1}{2}}\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)-{\mathrm{\ω}}_n({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)}---\left(8\right)$

$\left\{\begin{array}{c}{l}_N\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\frac{v_i^{N}N!}{R\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s))}\\ l_{N+1}\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=0\\ l_{n-1}\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\frac{1}{n}\[\frac{2n+1}{{\mathrm{\υ}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)}{\mathrm{\σ}}_n{l}_n\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)-(n+1)l_{n+1}\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)\]\end{array}\right.---\left(9\right)$

最后，通过Powell方向加速法计算式(1)的最小值，与该最小值对应处的光学参数即为波长λ下组织光学参数的反构值。

实施例：

卤钨灯1选择美国OCEAN OPTICS公司的HL-2000卤钨灯，发光波长范围为360nm～2000nm；光源光纤2的纤芯直径为400μm，数值孔径为0.37；各向同性散射头选择瑞典Medlight公司的IP85isotropic probe，散射头的直径为850μm；侧照光纤5的纤芯直径为600μm，数值孔径为0.22；光谱仪选用海洋光学HR2000光谱仪，其接收光波长范围为200nm～1100nm；利用Intralipid、India Ink及蒸馏水配制出特定光学参数的待测液态仿体，分别测量所述液态仿体在660nm、780nm和830nm三个波长下的光学参数为例。

按照上述方法，并从测量信号中，提取出660nm波长下的辐射率分布，并存储于上位计算机中，利用上位计算机按照步骤四提供的具体算法对所存储的信号进行处理，重建660nm波长下待测液态仿体的光学参数。

同理，分别重建780nm和830nm波长下待测液态仿体的光学参数。

尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。

Claims (1)

1.一种基于稳态辐射率间质测量装置的多波长光学参数反构方法，其中，所采用的稳态辐射率间质测量装置的结构是，包括光源模块、探测模块、位移模块及控制模块；所述光源模块包括连接至光源光纤(2)一端的卤钨灯(1)；所述光源光纤(2)的另一端设有各向同性散射头(3)，用于将所述卤钨灯(1)发出的光转为点光源；所述探测模块包括侧照光纤(5)和光谱仪(6)，所述侧照光纤(5)的一端耦合至所述光谱仪(6)，所述侧照光纤(5)的另一端为与光纤轴线呈45°角的探测斜面(4)，用于探测与光纤轴线垂直的入射的光信号；所述位移模块包括手动精密位移台(7)、电动平移台(8)和电动旋转台(9)，所述电动平移台(8)的上方固定有一金属立柱，所述手动精密位移台(7)固定在所述金属立柱的上方，所述光源光纤(2)通过一金属固定管与所述手动精密位移台(7)固定；所述电动旋转台(9)固定在一金属支架上，所述侧照光纤(5)通过一金属固定管垂直的固定在所述电动旋转台(9)上；所述手动精密位移台(7)用于调整光源光纤(2)的高度，使各向同性散射头(3)与探测斜面(4)保持在同一水平高度；所述电动平移台8用于在水平方向移动光源光纤(2)；所述电动旋转台(9)用于旋转侧照光纤(5)，探测斜面(4)的朝向以改变光子接收方向；所述控制模块包括步进电机控制箱(10)和上位计算机(11)；所述上位计算机|(11)一方面发出指令至步进电机控制箱(10)，从而控制电动平移台(8)及电动旋转台(9)的运动；另一方面，用于控制光谱仪(6)探测光信号；

其特征在于，该多波长光学参数反构方法包括以下步骤：

步骤一、利用手动精密位移台(7)、电动平移台(8)，调整光源光纤(20的空间位置，从而确定各向同性散射头(3)与探测斜面(4)之间的距离，即光源-探测器之间的距离r；通过电动旋转台(9)调整探测斜面(4)的初始探测方向；设定电动旋转台(9)的单位旋转角度为1°；

步骤二、所述卤钨灯(1)发出白光，经过光源光纤(2)传输至各向同性散射头(3)，光经过各向同性散射后，以点光源的形式传输至置于样品池(12)中的待测样品中；光经过待测样品的吸收散射后，经过探测斜面(4)传入侧照光纤(5)，被光谱仪(6)接收；光谱仪(6)测量到的信号传入上位计算机(11)，并存储；

步骤三、以步骤一设定的单位旋转角度旋转电动旋转台9，重复步骤二，直至实现在水平面内360°范围内辐射率的测量；

步骤四、波长λ下组织光学参数的反构，步骤如下：

首先，选取待反构光学参数所在波长为λ，提取出波长λ下的辐射率M＝[m₀,m₁,...,m_k,...,m₃₅₉]，其中，m_k为第k个角度下的辐射率测量值；

然后，设定最优化目标函数为：

$f({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)=\sqrt{{\Σ}_{k=0}^{359}{(\frac{m_k}{m_0}-\frac{l_k({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)}{l_0({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)})}^2}---\left(1\right)$

式(1)中；分别为样品池中待测样品的吸收系数、散射系数的估计值；是光源-探测器距离为r、探测角度为k/180×π处P9近似下辐射传输方程的解，即P9近似下的辐射率，其表达式为：

$l_k({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s)=\underset{n=0}{\overset{N}{\Σ}}\frac{2n+1}{4\mathrm{\π}}{\mathrm{\φ}}_n({\widetilde{\mathrm{\μ}}}_a,{\widetilde{\mathrm{\μ}}}_s,r)P_n\left(cos\right(\frac{k\·\mathrm{\π}}{180}\left)\right)---\left(2\right)$

式(2)中：n＝0,1,...,N，N＝9为近似阶数，k＝0,1,...,359为第k个探测角度；P_n为n阶勒让德多项式；φ_n(μ_a,μ_s,r)为第n阶光子密度计算值，其表达式为：

${\mathrm{\φ}}_n({\mathrm{\μ}}_a,{\mathrm{\μ}}_s,r)=\frac{1}{4\mathrm{\π}}\underset{i=1}{\overset{\frac{N+1}{2}}{\Σ}}{\mathrm{\υ}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)G_i\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)l_n\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right){\mathrm{\κ}}_n\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)r\right)---\left(3\right)$

其中，κ_n(x)的表达式如下：

$\left\{\begin{array}{c}{\mathrm{\κ}}_0\left(x\right)=\frac{\exp \left(x\right)}{x}\\ {\mathrm{\κ}}_1\left(x\right)=(1-\frac{1}{x})\frac{\exp \left(x\right)}{x}\\ {\mathrm{\κ}}_n\left(x\right)={\mathrm{\κ}}_{n-2}\left(x\right)-\frac{2n-1}{x}{\mathrm{\κ}}_{n-1}\left(x\right)\end{array}\right.---\left(4\right)$

定义递推式：B_n+1(ω(μ_a,μ_s))＝(2n+1)σ_nB_n(ω(μ_a,μ_s))+ω(μ_a,μ_s)n²B_n-1(ω(μ_a,μ_s))，其中，σ_n≡μ_a+μ_s(1-gⁿ)，g为组织体的各向异性因子，g取值为0.9；

当，B₁(ω(μ_a,μ_s))＝1，B₂(ω(μ_a,μ_s))＝3σ₁，时

$B_{N+1}\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=R\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\underset{l=0}{\overset{\frac{N-1}{2}}{\Σ}}{a}_l{\mathrm{\ω}}^l({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)---\left(5\right)$

当，B₀(ω(μ_a,μ_s))＝1，B₁(ω(μ_a,μ_s))＝μ_a，时

$B_{N+1}\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=Q\left(\mathrm{\ω}\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\underset{l=0}{\overset{\frac{N+1}{2}}{\Σ}}{b}_l{\mathrm{\ω}}^l({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)---\left(6\right)$

ω_i(μ_a,μ_s)为Q(ω(μ_a,μ_s))＝0的根，

${\mathrm{\υ}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)=\sqrt{-{\mathrm{\ω}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)}---\left(7\right)$

$G_i\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\frac{1}{b_{\frac{N+1}{2}}}\frac{R\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)}{{\Π}_{n=1,n\≠i}^{\frac{N-1}{2}}\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)-{\mathrm{\ω}}_n({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)}---\left(8\right)$

$\left\{\begin{array}{c}{l}_N\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\frac{v_i^{N}N!}{R\left({\mathrm{\ω}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)}\\ l_{N+1}\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=0\\ l_{n-1}\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)=\frac{1}{n}\[\frac{2n+1}{{\mathrm{\υ}}_i({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)}{\mathrm{\σ}}_n{l}_n\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)-(n+1)l_{n+1}\left({\mathrm{\υ}}_i\right({\mathrm{\μ}}_a,{\mathrm{\μ}}_s\left)\right)\]\end{array}\right.---\left(9\right)$

最后，通过Powell方向加速法计算式(1)的最小值，与该最小值对应处的光学参数即为波长λ下组织光学参数的反构值。