基于残缺P3近似正向模型测量混浊介质光学参数的方法

摘要

本发明公开了一种基于残缺P3近似正向模型测量混浊介质光学参数的方法，主要包括：建立辐射率的残缺P3近似正向模型，发展基于残缺P3近似正向模型光学参数反构的解析算法，从而得到被测混浊介质的吸收系数和约化散射系数。所发展的残缺P3与传统P3相比表达式更加简洁，与扩散近似相比更加准确描述光在大吸收或小散射组织中的传播。所发展的基于残缺P3解析反构算法，与传统的基于P3近似的拟合算法及基于扩散近似的解析反构算法相比，能够更加准确地反构出大吸收或小散射混浊介质的光学参数。其中的解析反构算法求解速度快有望应用于的光学参数在线监测。本发明所提出的测量方法不需要测量实验系统所用的光源强度，有利于实际应用。

说明书

技术领域

本发明属于组织光学研究中的光学参数测量领域，具体涉及一种残缺P3近似正向模型及基于此正向模型所发展的解析反构算法，获得混浊介质(组织体)的光学参数，吸收系数及约化散射系数。

背景技术

光学参数测量在生物医学光子学领域受到广泛关注。尤其是在激光治疗中，在体光学参数测量可用于计算光通量以达到在线优化光剂量的目的[1，2]。稳态测量技术由于设备简单、造价较低常用于混浊介质(组织体)的光学参数测量[3-5]。尽管表面漫反射率测量可以对皮肤等组织进行无创的光学参数测量，但由于光学穿透深度低该技术无法用于内部器官如前列腺的检测[3]。对于内部器官的检测需要通过微创测量方式来实现。目前使用的方法有空间分辨的光子密度测量方法及漫反射率光谱测量方法[4，5]。其中，光子密度测量方法需要通过积分球对系统的光源光强及探测器进行绝对矫正，不利于临床应用[6]。另外，漫反射率测量方式一般通过若干光纤耦合成的探针来测量若干源探距离(Sourcedetector separation，SDS)下的漫反射率，而有限源探距离下的测量数据限制了光学参数反构的准确性[7]。

稳态辐射率技术，在固定的SDS下测量不同方向的光强，为光学参数测量提供了一个新的自由度。该技术由于仅仅需要一至两个SDS下的角度辐射率的测量数据，减少了间质测量过程中由于探测器的平移给组织带来的损伤，因此逐渐发展成为一个新的微创测量技术[6，8，9]。目前，利用辐射率测量数据反构光学参数的算法主要有两种：基于P3近似的拟合算法及基于扩散近似(diffusion approximation，DA)的解析表达式求解方法[6，9，10]。其中，拟合算法仅能够反构出传输反照率a′＝μ′_s/(μ_a+μ′_s)大于0.9时的光学参数，这里的μ_a及μ′_s分别为吸收系数与约化散射系数。这是由于基于辐射率形状的拟合算法在求解光学参数过程中缺乏解的唯一性，即一系列的光学参数对(μ_a,μ′_s)下的辐射率的形状相同[9]。另外，尽管基于DA近似的解析表达式方法避免了拟合算法中求解的不唯一性，但该方法仅在满足DA近似条件即高a′下适用[6]。而在间质激光治疗中，会遇到a′低至0.5的情况[9]。这就需要发展一个较DA更加正确的正向模型，而且基于此正向模型能够发展一种反构解析算法，即通过解析表达式求解出混浊介质(组织体)的光学参数。

[参考文献]

1.G.Shafirstein,D.Bellnier,E.Oakley,S.Hamilton,M.Potasek,K.Beeson,andE.Parilov,“Interstitial Photodynamic Therapy-A Focused Review,”Cancers 9(2),12(2017)。

2.V.K.Nagarajan and B.Yu,“Monitoring of tissue optical propertiesduring thermal coagulation of ex vivo tissues,”Lasers Sur.Med.48(7),686–694(2016)。

3.S.-Y.Tzeng,J.-Y.Guo,C.-C.Yang,C.-K.Hsu,H.-J.Huang,S.-J.Chou,C.-H.Hwang,and S.-H.Tseng,“Portable handheld diffuse reflectance spectroscopysystem for clinical evaluation of skin:a pilot study in psoriasis patients,”Biomed.Opt.Express 7(2),616–628(2016)。

4.T.C.Zhu,J.C.Finlay,and S.M.Hahn,“Determination of the distributionof light,optical properties,drug concentration,and tissue oxygenation in-vivoin human prostate during motexafin lutetium-mediated photodynamic therapy,”J.Photoch.Photobio.B 79(3),231–241(2005)。

5.T.M.Baran,M.C.Fenn,T.H.Foster,“Determination of optical propertiesby interstitial white light spectroscopy using a custom fiber optic probe,”J.Biomed.Opt.18(10),107007(2013)。

6.S.Grabtchak and W.M.Whelan,“Separation of absorption and scatteringproperties of turbid media using relative spectrally resolved cw radiancemeasurements,”Biomed.Opt.Express 3(10),2371–2380(2012)。

7.T.M.Baran,J.D.Wilson,S.Mitra,J.L.Yao,E.M.Messing,D.L.Waldman,andT.H.Foster,“Optical property measurements establish the feasibility ofphotodynamic therapy as a minimally invasive intervention for tumors of thekidney,”J.Biomed.Opt.17(9),0980021(2012)。

8.D.J.Dickey,R.B.Moore,D.C.Rayner,and J.Tulip,“Light dosimetry usingthe P3 approximation,”Phys.Med.Biol.46(9),2359(2001)。

9.L.C.L.Chin,A.E.Worthington,W.M.Whelan,and I.A.Vitkin,“Determinationof the optical properties of turbid media using relative interstitialradiance measurements:Monte Carlo study,experimental validation,andsensitivity analysis,”J.Biomed.Opt.12(6),064027(2007)。

10.E.L.Hull and T.H.Foster,“Steady-state reflectance spectroscopy inthe P3 approximation,”J.Opt.Soc.Am.A 18(3),584–599(2001)。

发明内容

针对间质测量混浊介质(组织体)光学参数的需求，本发明提出了一种基于残缺P3近似正向模型的光学参数解析求解方法。

为了解决上述技术问题，本发明提出的一种基于残缺P3近似正向模型测量混浊介质光学参数的方法，包括以下步骤：

步骤一、建立辐射率的残缺P3近似正向模型，表达式如下：

式(1)中，r为源探距离，θ为探测方向相对于径向方向向量的角度，S₀为光源强度，P_l(cosθ)为勒让德多项式，Q_l(-νr)为修正的第二类球贝塞尔函数，以递推形式给出将该函数的表达式，如下：

式(1)中的ν为渐进衰减系数，其表达式为：

式(3)中，ξ＝27μ_aμ′_s+28μ_aσ₃+35σ₂σ₃，ζ＝3780μ_aμ′_sσ₂σ₃；其中，μ_a为吸收系数，μ_s为散射系数，σ_l为l阶吸收系数，定义σ_l＝μ_a+μ_s(1-g^l)，l＝2,3，其中，g为各向异性因子，μ′_s＝μ_a+μ_s(1-g)为约化散射系数；

式(1)中的h_l(ν)，l＝0,...,3，为：

式(1)中的C′为：

式(5)中，ν⁺为瞬时衰减系数，定义

步骤二、基于残缺P3近似正向模型，发展光学参数解析反构的算法，从而得到被测混浊介质的吸收系数μ_a和约化散射系数μ′_s，步骤如下：

步骤2-1、首先确定两个源探距离r和r′及分别对应的四个角度θ_i，i＝0,...,3，然后分别测量两个源探距离r和r′下四个角度θ_i，i＝0,...,3下的辐射率：L_m(r,θ_i)，L_m(r′,θ_i)；

步骤2-2、利用式(6)求解式(1)中的S₀C′h_l(ν)Q_l(-νr),l＝0,1,2，

式(6)右边简写为f_l(r,θ₀,θ₁,θ₂,θ₃)，l＝0,1,2，式(6)中，α_i(r,θ₀,θ₁,θ₂)、α_i(r,θ₀,θ₁,θ₃)、β_ij(θ₀,θ₁,θ₂)和β_ij(θ₀,θ₁,θ₃)为：

步骤2-3、求解渐进衰减系数ν：

构建f₀(r,θ₀,θ₁,θ₂,θ₃)与f₀(r′,θ₀,θ₁,θ₂,θ₃)之比：

式(8)中Q₀(-νr′)＝(-νr′)^-1exp(-νr′)，求解式(8)得到渐进衰减系数ν为：

步骤2-4、求解h₁(ν)和h₂(ν)：

构建f_l(r,θ₀,θ₁,θ₂,θ₃)，l＝1,2与f₀(r,θ₀,θ₁,θ₂,θ₃)之比：

通过式(10)求得h₁(ν)和h₂(ν)为：

步骤2-5、根据式(4)得到被测混浊介质的吸收系数μ_a和约化散射系数μ′_s：

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

(1)本发明中建立的残缺P3近似正向模型，与传统的扩散近似正向模型相比能更加准确地描述光在大吸收或小散射组织中的传播。与传统的P3近似正向模型相比，其表达式更加简洁，有利用构建反构光学参数的解析表达式。

(2)本发明基于残缺P3近似正向模型所发展的光学参数解析反构算法，与传统的基于P3近似正向模型的拟合算法及基于扩散近似正向模型的解析反构算法相比，能够更加准确地反构出大吸收或小散射混浊介质的光学参数。

(3)传统的拟合算法需要大量的实验数据进行拟合求解光学参数，而本发明中所述的反构算法仅需要8个辐射率测量值，所需的实验测量值少，有望于设计为光纤探针，用于测量混浊介质(组织体)的光学参数。

(4)本发明可以反构单波长下混浊介质(组织体)的光学参数，因而在间质激光治疗中可以利用相同的激光进行光学参数测量及治疗，更加有利用临床应用。

(5)本发明所提出的解析反构算法，利用解析式从测量数据直接求解出混浊介质(组织体)的光学参数，不需要传统拟合算法中的迭代求解，因而求解速度快，有望应用于间质激光治疗中的光学参数在线监测。

(6)本发明所提出的解析反构算法，和绝对测量相比不需要测量实验系统所用的光源强度，更加有利用临床应用。

附图说明

图1示出了本发明测量方法中所确定的两个源探距离和四个角度的示意图。

具体实施方式

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

实施例：以辐射率测量值L_m(r,θ_i),L_m(r′,θ_i)i＝0,...,3，被测介质为利用印度墨水及20％Intralipid所配置的μ_a＝0.5mm^-1，μ′_s＝0.5mm^-1的液态仿体为例，其中r和r′选为5mm和10mm，θ₀，θ₁，θ₂及θ₃选为40°、80°、120°和160°，所使用的测量系统是稳态辐射率间质测量装置系统，主要结构(参见专利申请号为201610822322.3)是：包括光源模块、探测模块、位移模块及控制模块；所述光源模块包括连接至光源光纤一端的卤钨灯；所述光源光纤的另一端设有各向同性散射头，用于将所述卤钨灯发出的光转为点光源；所述探测模块包括侧照光纤和光谱仪，所述侧照光纤的一端耦合至所述光谱仪，所述侧照光纤的另一端为与光纤轴线呈45°角的探测斜面，用于探测与光纤轴线垂直的入射的光信号；所述位移模块包括手动精密位移台、电动平移台和电动旋转台，所述电动平移台的上方固定有一金属立柱，所述手动精密位移台固定在所述金属立柱的上方，所述光源光纤通过一金属固定管与所述手动精密位移台固定；所述电动旋转台固定在一金属支架上，所述侧照光纤通过一金属固定管垂直的固定在所述电动旋转台上；所述手动精密位移台用于调整光源光纤的高度，使各向同性散射头与探测斜面保持在同一水平高度；所述电动平移台用于在水平方向移动光源光纤；所述电动旋转台用于旋转侧照光纤，确定探测斜面的朝向以改变光子接收方向；所述控制模块包括步进电机控制箱和上位计算机；所述上位计算机一方面发出指令至步进电机控制箱，从而控制电动平移台及电动旋转台的运动；另一方面，用于控制光谱仪探测光信号。

利用本发明方法对上述实施例中被测介质进行光学参数测量的过程如下：

步骤一、建立辐射率的残缺P3近似正向模型：

P3近似用于描述无限媒质下点光源激励下的辐射率测量[10]。由于P3近似表达式的复杂，无法构建出反构光学参数的解析表达式。本发明提出了一种表达式较为简单的残缺P3(incomplete P3approximation)，简写为P3_in，正向模型。P3_in为P3近似的渐近衰减解部分，表达式如下：

式(1)中，r为源探距离，θ为探测方向相对于径向方向向量的角度，S₀为光源强度，P_l(cosθ)为勒让德多项式，Q_l(-νr)为修正的第二类球贝塞尔函数，以递推形式给出将该函数的表达式，如下：

式(1)中的ν为渐进衰减系数，其表达式为：

式(3)中，ξ＝27μ_aμ′_s+28μ_aσ₃+35σ₂σ₃，ζ＝3780μ_aμ′_sσ₂σ₃；其中，μ_a为吸收系数，μ_s为散射系数，σ_l为l阶吸收系数，定义σ_l＝μ_a+μ_s(1-g^l)，l＝2,3，其中，g为各向异性因子，μ′_s＝μ_a+μ_s(1-g)为约化散射系数；

式(1)中的h_l(ν)，l＝0,...,3，为：

式(1)中的C′为：

式(5)中，ν⁺为瞬时衰减系数，定义

步骤二、基于残缺P3近似正向模型，发展光学参数解析反构的算法，从而得到被测混浊介质的吸收系数μ_a和约化散射系数μ′_s，步骤如下：

步骤2-1、首先确定r′＝5mm和r＝10mm两个源探距离及分别对应的θ₀＝40°、θ₁＝80°、θ₂＝120°和θ₃＝160°四个角度，如图1所示，然后分别测量两个源探距离r和r′下四个角度θ_i，i＝0,...,3下的辐射率：L_m(r,θ_i)，L_m(r′,θ_i)；

步骤2-2、利用式(6)求解式(1)中包含ν,h₁(ν)和h₂(ν)的S₀C′h_l(ν)Q_l(-νr),l＝0,1,2，

式(6)右边简写为f_l(r,θ₀,θ₁,θ₂,θ₃)，l＝0,1,2，式(6)中，α_i(r,θ₀,θ₁,θ₂)、α_i(r,θ₀,θ₁,θ₃)、β_ij(θ₀,θ₁,θ₂)和β_ij(θ₀,θ₁,θ₃)为：

步骤2-3、求解渐进衰减系数ν：

为了求解ν且同时避免绝对测量即消除光强S₀，构建f₀(r,θ₀,θ₁,θ₂,θ₃)与f₀(r′,θ₀,θ₁,θ₂,θ₃)之比：

式(8)中Q₀(-νr′)＝(-νr′)^-1exp(-νr′)，求解式(8)得到渐进衰减系数ν为：

步骤2-4、求解h₁(ν)和h₂(ν)：

为了求解h₁(ν)和h₂(ν)，构建f_l(r,θ₀,θ₁,θ₂,θ₃)，l＝1,2与f₀(r,θ₀,θ₁,θ₂,θ₃)之比：

公式(9)已经求得ν，因此Q_l(-νr)(l＝0,1,2)为已知，代入公式(10)求得h₁(ν)和h₂(ν)为：

步骤2-5、根据式(4)得到被测混浊介质的吸收系数μ_a和约化散射系数μ′_s：

最终求得液态仿体的吸收系数及约化散射系数反构结果为：μ_a＝0.4852mm^-1，μ′_s＝0.5096mm^-1，相对误差为：2.96％及1.92％。

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