CN101856219B - 基于频域近红外光测量的光学参数重构方法 - Google Patents

Abstract

本发明属于组织光学研究中的光学参数测量领域，涉及一种基于频域近红外光测量的光学参数重构方法：运用时域蒙特卡洛模拟生物组织模型以进行正向建模，运用快速公式进行频域信息的提取，采用基于二元多项式的快速MC模型来获得任意光学参数下的频域信息，结合L-M优化算法达到对管状组织的光学参数精确重构的目的，并引入聚类分析方法提高重构精度。本发明能同时重建组织的吸收系数与散射系数，并兼具重建速度与重建精度的优势，可以通过测量到的频域信息，即幅值和相位(A，φ)对测量样品的光学参数进行实时、准确的预测。

Description

基于频域近红外光测量的光学参数重构方法

技术领域

本发明属于组织光学中的光学参数测量领域，具体涉及一种基于频域近红外光测量的光学参数重构方法。

技术背景

近年来***发病率增高严重威胁女性健康。大量研究表明早期、及时的诊断对于***的治疗至关重要。

近红外光学法作为早期***的新的诊断方法，与传统检测方法相比，具有无创、快速、测量准确等优点。它通过对生物组织体光学参数(包括吸收系数μ_a和散射系数μ_s)的准确检测，为疾病诊断提供可靠依据。

生物组织光学参数的近红外测量方法主要分为连续光测量方法，时域测量方法⁽¹⁾和频域测量方法⁽²⁾。由于频域***的廉价，精确等特性，我们选取在频域中进行测量。在频域测量***中，入射光强度被调制到几百MHz，当此调制光通过组织体后，可通过诸如外差法获得正弦信号幅度的衰减及相位延迟.频域***因设备简单、测量速度快而在***漫射光诊断中获得重视。基于频域测量的重构过程即是由测量值(幅值信息A和相位信息Ф)获得组织体光学参量(吸收系数μ_a和散射系数μ_s)的过程。

光学参量重构(反演)建立在描述光在组织体中传输模型(正向模型)的基础上。生物组织体光学参数分布任意，边界条件复杂。采用辐射传输方程(Radiative Transfer Equation，RTE)来解决光在随机介质中的传播问题无法求出精确解析解，于是提出各种近似解法，最常用的为漫射近似⁽³⁾(也称P1近似)模型。然而，漫射近似在解决如宫颈组织之类复杂几何模型的光传输问题上，精度不高^(4，5)，不存在优势。

蒙特卡洛(Monte-Carlo，MC)方法则更加灵活易用，它是一种统计模拟随机抽样的方法，主要用于模拟各种输运现象。它根据待求解问题或物理现象本身的变化和统计规律，构造出一个合适的概率模型或随机过程，通过大量统计实验来计算所求参数。MC模型描述光在组织中的传播具有以下优势：(1)实现简单，模拟准确，其运行过程只需组织体的几个局部光学参数；(2)可模拟任意几何形状、边界条件和光学参数分布下的光传输行为；(3)能够在输出量中自然地引入本征播送噪声性能，从而有效地建立试验统计模型。

参考文献：

(1)Shay Keren，Olivier Gheysens，Craig S.Levin，et al.A Comparison Between a Time Domain andContinuous Wave Small Animal Optical Imaging System，IEEE Transactions On Medical Imaging，2008，27(1)：58-63.

(2)Anand T.N.Kumar，Scott B.Raymond，Brian J.Bacskai，et al.Comparison of frequency-domain andtime-domain fluorescence lifetime tomography，Optics Letters，2008，33(5)：470-472.

(3)谢树森，雷仕湛，光子技术，北京：科学出版社，2004：214-253.

(4)吕可诚，张春平，张光寅等，生物光子学进展，光子学报，1997，26，12：1123-1129.

(5)苏畅，丁海曙，白净，生物组织光谱学技术，国外医学生物医学工程分册，1995，18，6：316-323.

(6)Tom Collier，Dizem Arifler，Anais Malpica，Determination of Epithelial Tissue Scattering CoefficientUsing Confocal Microscopy，Journal Of Selected Topics In Quantum Electronics，2003，9(2)：307-313.

发明内容

本发明的目的是提供一种新的组织光学参数重构方法，该方法能够准确、快速提取获得任意组织光学参数下的频域信息，并在逆向求解中，具有较快的局部收敛速度和重构精度。

本发明的技术方案如下：

一种基于频域近红外光测量的光学参数重构方法，包括下列步骤：

(1)、针对待检测生物组织的形状建立模型，进行时域正向时域蒙特卡洛模拟，获得组织的时域信息；

(2)、采用如下的方法从正向时域蒙特卡洛模拟结果中提取幅值和相位信息：已知每个光子包的权重W_d、平均飞行时间t_d，采用计算某一频率ω₀下的傅立叶变换，再应用式

和式

分别得出射光的幅值和相位，式中，G_R，G_I分别为

的实部和虚部，建立频域信息数据库；

(3)对已建立的数据库中的频域信息进行聚类分析，选取测量值幅度A和相位Ф以及组织体光学参量吸收系数μ_a和散射系数μ_s为聚类过程中的4个变量，采用最短距离法的原则对标准化后的变量进行聚类，聚类完成后，分别对每一大类求取其各自的幅值的平均值

和相位的平均值

吸收系数μ_a的平均值

和散射系数μ_s的平均值

将所得

及

作为类的中心点，并按照距离公式

计算待测点距各类中心的距离，，获取满足d最小的那一类，该类的类中心的

即为待测点在重构过程中的初始光学参数值，在距离公式中，

与分别代表待测点的幅值和相位信息，n是类的个数，δ_A和

分别为数据库中幅值和相位信息的标准差，α为0到1之间的常数；

(4)以步骤(3)得到的初始光学参数值作为试探值，采用L-M最优化方法，进行光学参数重构，先由快速MC模型得到初始光学参数所对应的模拟频域信息，计算模拟信息与待测点频域信息之间的误差，根据误差修改每一次迭代过程中的光学参数试探值，再由快速MC模型得到相应的模拟频域信息，并重复进行上述步骤，直至误差达到最小。

作为优选实施方式，步骤(4)中采用的快速MC模型为：根据已建立的频域信息数据库的数据，运用二元多项式模型进行曲线拟合，得到从(μ_a，μ_s)到(A，φ)的直接映射关系，从而可利用拟合函数快速获得L-M算法在重建过程中需要的任意光学参数下的模拟频域信息；二元多项式曲线拟合的阶数可以为5；步骤(4)计算误差时，通过加权的方法，削弱相位信息在重构过程中的作用。

本发明的有益效果在于：

(1)本发明采用蒙特卡洛模拟作为组织光学传输的正向模型，加入了生物组织具体的几何形状的考虑，针对空心圆柱对其边界做了处理。

(2)采用快速公式方法提取时域MC结果的频域信息，不需要建立时间柱状图，避免了由时间扩展曲线时间区间宽度带来的时间抽样造成的时间误差问题，幅值和相位的计算精度更高。

(3)本发明采用了基于二元多项式拟合的快速蒙特卡洛方法结合L-M算法的反演过程。该方法的重构效率很高，且适用于两个光学参数的同时重构。

(4)本发明通过模拟验证，选用阶数n＝5的二元多项式来对频域曲面进行回归拟合，达到较高的拟合精度。

(5)引入聚类分析方法来确定L-M算法中的迭代初值，以获得L-M算法的全局最优值，使重构精度大幅提高。

附图说明

图1本发明的重构方法的流程图。

图2(a)为幅值信息随光学参数变化曲面图；

图2(b)为相位信息随光学参数变化曲面图。

具体实施方式

本发明提供一种组织光学参数重构方法。考虑到本发明的研究对象是几何形状复杂，边界条件特殊的生物体管状组织，如空心圆柱形宫颈组织，本发明选择MC模型作为描述组织中光的传播过程的正向模型，以获得更加准确的模拟结果。

由于已知组织的光学参数描述是光在生物组织中的传输的基础，因此光在组织中的正向传输模型是组织光学参数反演问题即逆问题的基础。

正向模型的求解是光学参量反演中需多次计算的过程，其运算速度决定了光学参量反演速度由于MC的巨大计算量，基于MC正向模型的光学参量反演在临床应用中面临着难以实时给出反演结果的问题。因此，快速、准确的求得MC模拟结果的方法显得尤为重要。

解决逆问题的方法很多，基本的如最速下降法，牛顿法，信赖域法等。最速下降法有很好的整体收敛性，但收敛速度太慢，并不实用；牛顿法有很快的收敛速度，但只是局部收敛，而且其二阶黑森矩阵通常难以计算或花费的时间很大。针对我们的优化问题，本发明选择了将最速下降法和牛顿法结合在一起的算法：Levenberg-Marqurdt(L-M)算法⁽⁶⁾。

L-M算法是用以解决最小二乘问题的无约束最优化方法，具有全局收敛的性质及较快的局部收敛速度，为该领域内的常用反演算法。然而，在重构过程中发现，L-M算法的结果和初值选取有很大关系，即该算法在处理最小二乘问题时不具有全局最优性质，重构值往往会陷入局部最优，很大程度上影响了重构精度。因此，选取合理的迭代初值，得到全局最优解，也是重构问题面临的主要难题。

本发明中所涉及的方法进行组织光学参数重构时，首先根据模拟数据库聚类分析的结果，将待重建的频域信息归入相应的类，选取适当的的初值。再利用L-M优化算法进行逆蒙特卡洛重建。在重建过程中，采用二元多项式拟合方法结合Beer-Lambert定理来实现快速MC模型，在迭代中快速得到任意光学参数下的频域信息。重建获得相应的光学参数值。整个重构算法流程图见附图1。

(1)、模拟频域信息数据库的建立

1、针对宫颈组织的时域MC模拟：

本实施例根据宫颈组织的形状，选取了一个空心圆柱体来模拟宫颈组织；根据宫颈组织的尺寸，选取的空心圆柱体内径1cm，外径2cm。由于长度远大于半径，假设其为无限长。

蒙特卡洛的一般模拟过程：

1.根据入射条件确定起始跟踪点；

2.确定光子运动的步长；

3.确定光子行进的方向和下一次碰撞的位置；

4.确定在该位置光子的吸收和散射部分；

5.返回第二步。

如此循环计算，直到光子权重小于某一设定值，或者光子逸出生物组织表面时结束对该光子的跟踪。然后返回第一步记录另一光子，直到所设定的光子数全部跟踪完毕。

在时域蒙特卡洛模拟中，光子权重根据光子各自的平均飞行时间，被累加到不同的时间段内进行记录，进而在每个探测点形成一条时间扩展曲线。

下面重点闸释类似宫颈组织的管状组织的边界处理，给出空心圆柱模型的边界反射是如何实现的。

在边界上，非偏振光的反射率由菲涅尔反射系数给出：

$R\left({\mathrm{\α}}_i\right)=\frac{1}{2}[\frac{{\sin }^2({\mathrm{\α}}_i-{\mathrm{\α}}_t)}{{\sin }^2({\mathrm{\α}}_i+{\mathrm{\β}}_t)}+\frac{{\tan }^2({\mathrm{\α}}_i-{\mathrm{\α}}_t)}{{\tan }^2({\mathrm{\α}}_i+{\mathrm{\α}}_t)}]---\left(1\right)$

其中，α_i为入射角，α_t为出射角，n_i、n_t分别是边界两边介质的折射率。

光子包的剩余权重与碰撞点的漫反射率R_d(α_t)或漫透射率T_t(α_t)的乘积，就是此光子包的反射分量或透射分量。如果反射分量的权重值大于阈值权重，则光子包继续传输。这时我们需要计算其碰撞后的传输方向和步长。

下面我们以外圆柱面为例，详细介绍计算过程。内圆柱面与其类似。

i.反射后的光子步长

根据光子包之前的传输方向和所处的位置，计算光子包到达边界的程长Δl；剩余程长l′＝l-Δl就是反射后的光子步长

当追踪光子包发现其下一个落点在组织之外时，就要进入边界反射的程序。令边界外的落点的坐标为(x，y，z)，那上一个落点的坐标就是(x-ΔS·u_x，y-ΔS·u_y，z-ΔS·u_z)，其中ΔS为该步的步长。有了空间上的两点就可以确定穿过该两点的直线方程，如下式：

$\frac{X-x}{u_x\mathrm{\ΔS}}=\frac{Y-y}{u_y\mathrm{\ΔS}}=\frac{Z-z}{u_z\mathrm{\ΔS}}=k_1---\left(2\right)$

同时我们设该直线与外圆柱面的交点坐标为(x₁，y₁，z₁)，该点也满足以上直线方程，代入其中：

$\frac{x_1-x}{u_x\mathrm{\ΔS}}=\frac{y_1-y}{u_y\mathrm{\ΔS}}=\frac{z_1-z}{u_z\mathrm{\ΔS}}=k_1---\left(3\right)$

又因(x₁，y₁，z₁)在外圆柱面上，可知下式成立：

x₁ ²+z₁ ²＝4    (4)

联立以上三式，解得k₁为：

$k_1=\frac{-(u_{x}x+u_{z}z)\±\sqrt{[{(u_{x}x+u_{z}z)}^2-(x^2+z^2-4)({u_x}^2+{u_z}^2)]}}{\mathrm{\ΔS}({u_x}^2+{u_z}^2)}---\left(5\right)$

可以证明，当±取+时，该k1对应于外圆柱面的界面反射；反之，对应于内圆柱面。于是可以得到连接前后两个落点的直线与外圆柱面的交点的坐标：

x₁＝k₁u_xΔS+x    (6)

y₁＝k₁u_yΔS+y    (7)

z₁＝k₁u_zΔS+z    (8)

所以反射后剩余的步长ΔS₁为：

$\mathrm{\Δ}{S}_1=\sqrt{{\left(k_1{u}_x\mathrm{\ΔS}\right)}^2+{\left(k_1{u}_y\mathrm{\ΔS}\right)}^2+{\left(k_1{u}_z\mathrm{\ΔS}\right)}^2}---\left(9\right)$

ii.反射后的方向

由于确定光子包反射后的方向是一个在立体空间里的柱面反射问题，因此会更为复杂一些，主要的思路是确定法线，然后由入射线和法线确定入射面，根据反射原理，反射线应位于此入射面中。具体步骤为先由两条线确定反射面。这两条线分别是组织内落点(x-ΔS·u_x，y-ΔS·u_y，z-ΔS·u_z)和组织外落点(x，y，z)所确定的入射线，以及这条线与圆柱面的交点(x₁，y₁，z₁)和该点所在的圆柱横截面的圆心(0，y₁，0)所确定的法线。这两条线可以确定一个平面。再求出通过反射点(x₁，y₁，z₁)的该平面法线，由几何知识得知该法线垂直于这一平面内的任意直线，自然也包括我们要求解的反射路径所在的那条线，由此得出式子(10)

u_yz₁u_x′-(u_xz₁-u_zx₁)u_y′-u_yx₁u_z′＝0    (10)

令反射路径的方向余弦为u_x′、u_y′和u_z′，入射角θ为：

$\mathrm{\θ}=\arccos \left(\frac{u_x{x}_1+u_z{z}_1}{r}\right)---\left(11\right)$

反射路径所在的这条直线又与确定反射平面的那两条直线分别成θ角和2θ角，因而得到式子(12)和(13)

x₁u_x′+z₁u_z′＝-r cosθ    (12)

u_xu_x′+u_yu_y′+u_zu_z′＝-cos(2θ)    (13)

联立以上两式，可以得到我们要求的出射方向余弦：

${u_z}^{\′}=\frac{\frac{u_y{z}_{1}r\cos \mathrm{\θ}}{x_1}-\frac{(u_z{x}_1-u_x{z}_1)(-\cos 2\mathrm{\θ}+\frac{u_x}{x_1}r\cos \mathrm{\θ})}{u_y}}{-(\frac{u_y{z_1}^2}{x_1}+u_y{x}_1)+\frac{(u_z{x}_1-u_x{z}_1)(\frac{u_x{z}_1}{x_1}-u_z)}{u_y}}---\left(14\right)$

其中r根据内外圆不同，可以取1cm或2cm。

再将该式代入以下两式：

${u_x}^{\′}=\frac{-r\cos \mathrm{\θ}-z_1{u_z}^{\′}}{x_1}---\left(15\right)$

${u_y}^{\′}=\frac{-\cos 2\mathrm{\θ}+(\frac{u_x{z}_{\text{1}}}{x_1}-u_z){u_z}^{\′}+\frac{u_{x}r\cos \mathrm{\θ}}{x_1}}{u_y}---\left(16\right)$

从得到反射路径的方向余弦u_x′、u_y′和u_z′，再加上前面所求得的剩余步长ΔS₁，就可以得知边界反射后的落点坐标了。

内圆柱面反射的方法与外圆柱面相同，只是入射角要取π-θ，选r值时取1cm，且求k₁的公式中±号要选减。

2、频域信息的提取：

本发明使用下述快速公式方法从时域MC模拟结果中提取幅值和相位信息。

当已知每个光子包的权重W_d、平均飞行时间t_d，可采用式(17)计算某一频率ω₀时的傅立叶变换，再应用(18)和(19)式即可精确得出射光的幅值和相位。

${\tilde{u}}_{\mathrm{sc}}\left({\mathrm{\ω}}_0\right)=\underset{d}{\mathrm{\Σ}}{W}_d\exp \left(i{\mathrm{\ω}}_0{t}_d\right)---\left(17\right)$

$A=\sqrt{{G_R}^2+{G_I}^2}---\left(18\right)$

$\mathrm{\Φ}=\arctan \left(\frac{G_I}{G_R}\right)---\left(19\right)$

其中G_R，G_I分别为的实部和虚部。

本发明利用上述方法获得了μa＝0.10，0.15，0.20，......0.80，μs＝30，35，40，45，......100共225组光学参数的频域信息，建立了一定光学参数范围下的频域信息数据库。由于不同光学参数下对应的的频域信息(幅值A和相位Ф)均成连续、平滑的曲面(见附图2)，所以为以下使用数学拟合方法提供了合理依据。

(2)、聚类分析方法与初值选取

1、聚类分析方法的引入

由于在重建过程中采用的一般优化算法均为局部优化算法，本发明中所使用的L-M优化方法也不例外。该方法只能保证求得局部最优点，其优化结果与初值选取有关。迭代设置的初始值不同，光学参数重构精度也发生变化。待重构的光学参数值距离迭代的初始值越近，重构精度越高，反之距离迭代的初始值越远，重构精度越低。

由此，我们设想如果能找到一种方法使得根据不同的光学参数，设置与之距离较近的迭代初始值从而提高重构精度。基于这一思想选择了聚类分析这种方法。本发明的聚类方法运用的是***聚类方法。***聚类方法主要包括最短距离法、最长距离法、中间距离法、类平均法、重心法以及离差平方和法(Ward方法)这六大类。为保证尽量将距离较近的幅值、相位所对应的多组光学参数划分在一类而且划分类的数量不能过多，本实施例选取了最短距离法。

以下用d_ij表示第i个样本与第j个样本的距离，G₁，G₂，......表示类，D_KL表示G_K与G_L之间的距离。一开始开始每个样本自成一类，类与类之间的距离与样本之间的距离相同，即D_KL＝d_KL，所以最初的距离矩阵全部相同，记为D(0)＝(d_ij)。

2、聚类的步骤与分类中心的确定：

首先，提取μ_a(0.1～0.8)步长为0.05；μ_s(30～100)步长为5的对应的225组光学参数的幅值和相位。

其次，为了保证频域信息相近的点的光学参数也尽可能相近，我们选取A，Ф，μ_a与μ_s作为分类的4个依据，即聚类过程中的4个变量，这四个变量在聚类时先要进行标准化，再按照最短距离法的原则进行聚类。类的个数可视待聚类的数据量而定。

再次，聚类完成后，分别对每一大类求取其各自的幅值的平均值和相位的平均值

求取μ_a的平均值

和μ_s的平均值

将所得及

作为类的中心点。这样各类的中心就确定下来。

接下来，就是判别待测点的频域信息所属的类。待测点距各类中心的距离用下式计算：

$d=\underset{1\≤i\≤n}{\min }\sqrt{{(\tilde{A}-\stackrel{\‾}{A})}^2+{(\widetilde{\mathrm{\φ}}-\stackrel{\‾}{\mathrm{\φ}})}^2}---\left(20\right)$

式中

与

分别代表待测点的幅值和相位信息，n是类的个数。

我们发现，在聚类过程中幅值信息的影响较相位信息大，为此，在本发明中，我们将上式(20)的距离公式改进为：

$d=\underset{1\≤i\≤n}{\min }\sqrt{{(\tilde{A}-\stackrel{\‾}{A})}^2/+{\mathrm{\δ}}_A+\mathrm{\α}*{(\widetilde{\mathrm{\φ}}-\stackrel{\‾}{\mathrm{\φ}})}^2/{\mathrm{\δ}}_{\mathrm{\φ}}}---\left(21\right)$

式中α为0到1之间的常数。δ_A和

分别为数据库中幅值和相位信息的标准差。待测点频域信息所属的类，即为使式(21)中的d达到最小的那一类，该类中心的

即为重构过程中的初始光学参数值。

将测量数据带入算法中进行光学重构，程序会自动将数据分入与初始值距离最近的一组，在迭代过程中初始值的频域信息与被测数据较为接近，保证边缘数据重构的精确性提高，整体重构精度被大大增高。

本发明经过多次模拟验证，考虑了重构误差大的某些区域，最终确定类的个数为42个。(3)、使用L-M算法进行快速光学参数重构：

所谓光学参数的反演问题就是利用特定的光传播模型F，求解p使得在光传输模型下模拟测量结果

与实际测量结果y之间的误差达到最小。这一问题可以采用最优化方法来求解。

Levenberg-Marquardt算法(简称L-M算法)是一种将最速下降法和高斯牛顿法结合在一起的优化算法，克服了最速下降法收敛速度过慢和高斯牛顿法局部收敛并且黑森矩阵难以计算的缺点，具有较快速的全局收敛性，适用于本发明的光学参数重构。

L-M算法的基本形式如下：

p_k+1＝p_k-(J(p_k)^TJ(p_k)+λdiag[J(p_k)^TJ(p_k)])^-1J(p_k)^Tf(p_k)    (22)

针对本发明的频域测量问题，光学参数的重构问题即通过设置试探值(μ_a ⁽⁰⁾，μ_s ⁽⁰⁾)^T，计算其对应的模拟结果(A⁽⁰⁾，Ф⁽⁰⁾)^T，与测量值(A_mea，Ф_mea)^T之间的误差，并根据误差修改试探值，观察修改后模拟结果与测量值的误差，并重复进行上述步骤，直到误差达到最小。具体过程如下：

设初始光学参数为修正后的光学参数应为：

$p_{k+1}=\left(\begin{array}{c}{\mathrm{\μ}}_a^{\left(1\right)}\\ {\mathrm{\μ}}_s^{\left(1\right)}\end{array}\right)=\left(\begin{array}{c}{\mathrm{\μ}}_a^{\left(0\right)}\\ {\mathrm{\μ}}_s^{\left(0\right)}\end{array}\right)+\left(\begin{array}{c}\mathrm{\Δ}{\mathrm{\μ}}_a\\ \mathrm{\Δ}{\mathrm{\μ}}_s\end{array}\right)---\left(23\right)$

其中修正量

由L-M算法计算得出。根据L-M算法式，

雅克比矩阵为

其中的偏导数由定义求出

$\frac{\∂A}{\∂{\mathrm{\μ}}_a}=\frac{A({\mathrm{\μ}}_a+\mathrm{\Δ}{\mathrm{\μ}}_a,{\mathrm{\μ}}_s)-A({\mathrm{\μ}}_a,{\mathrm{\μ}}_s)}{\mathrm{\Δ}{\mathrm{\μ}}_a}---\left(26\right)$

类似的，也可以求出

等量。

由此计算出黑森矩阵为

代入式(27)，得

根据克莱姆法则，则两根为：

$p_{k+1}-p_k=\left(\begin{array}{c}\mathrm{\Δ}{\mathrm{\μ}}_a\\ \mathrm{\Δ}{\mathrm{\μ}}_s\end{array}\right)=\left(\begin{array}{c}\frac{\left|D_1\right|}{\left|D\right|}\\ \frac{\left|D_2\right|}{\left|D\right|}\end{array}\right)---\left(30\right)$

其中：

重复以上过程，直到重构结果收敛或满足迭代的终止条件为止。式(22)中λ为L-M算法中的有效下降因子，本发明中取λ＝0.001

本发明在解决光学参数的重构问题中，采用了替代传统MC模拟方法的快速求得MC模拟的近似值的方法，以实现对光学参数实时重构的目的，解决逆蒙特卡洛模拟问题中的关键问题。本发明中，采用二元多项式回归方法，快速获得L-M算法在重建过程中需要的任意光学参数下的频域信息，代替耗时的正向MC模拟。下面对该方法做详细说明：

该方法利用形式如下的二元多项式函数建立数学模型，对已建立的频域信息数据库进行曲面拟合。

P(μ_a，μ_s，n)＝(a₀+a₁μ_a+a₂μ_a ²+......+a_nμ_a ⁿ)×(b₀+b₁μ_s+b₂μ_s ²+......+b_nμ_s ⁿ)    (34)

上式为两个同阶数的一元多项式相乘而得，表示变量P是μ_a，μ_s的二元函数。式(34)中，系数a₀，a₁，a₂，......a_n和b₀，b₁，b₂，......b_n满足最小二乘回归，n为二元多项式的阶数。

针对幅值和相位信息，上述二元多项式函数可以写为：

A(μ_a，μ_s，n)＝(a₀+a₁μ_a+a₂μ_a ²+......+a_nμ_a ⁿ)×(b₀+b₁μ_s+b₂μ_s ²+......+b_nμ_s ⁿ)    (35)

φ(μ_a，μ_s，n)＝(c₀+c₁μ_a+c₂μ_a ²+......+c_nμ_a ⁿ)×(d₀+d₁μ_s+d₂μ_s ²+......+d_nμ_s ⁿ)    (36)

其中，a₀，a₁，a₂，......a_n与b₀，b₁，b₂，......b_n由幅值决定，c₀，c₁，c₂，......c_n与d₀，d₁，d₂，......d_n由相位决定。

利用已建立的离散频域信息数据库，选择这225组光学参数，及其各点的幅值相位信息作为已知信息。进行二元多项式曲面拟合，返回形式如下的幅值系数矩阵p_A与相位系数矩阵p_φ：

$p_A=\left[\begin{array}{ccccccc}{a}_n{b}_n& a_n{b}_{n-1}& ...& a_n{b}_3& a_n{b}_2& a_n{b}_1& a_n{b}_0\\ a_{n-1}{b}_n& a_{n-1}{b}_{n-1}& ...& a_{n-1}{b}_3& a_{n-1}{b}_2& a_{n-1}{b}_1& a_{n-1}{b}_0\\ .& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ a_3{b}_n& a_3{b}_{n-1}& ...& a_3{b}_3& a_3{b}_2& a_3{b}_1& a_3{b}_0\\ a_2{b}_n& a_2{b}_{n-1}& ...& a_2{b}_3& a_2{b}_2& a_2{b}_1& a_2{b}_0\\ a_1{b}_n& a_1{b}_{n-1}& ...& a_1{b}_3& a_1{b}_2& a_1{b}_1& a_1{b}_0\\ a_0{b}_n& a_0{b}_{n-1}& ...& a_0{b}_3& a_0{b}_2& a_0{b}_1& a_0{b}_0\end{array}\right]---\left(37\right)$

$p_{\mathrm{\φ}}=\left[\begin{array}{ccccccc}{c}_n{d}_n& c_n{\text{d}}_{n-1}& ...& c_n{d}_3& c_n{d}_2& c_n{d}_1& c_n{d}_0\\ c_{n-1}{d}_n& c_{n-1}{d}_{n-1}& ...& c_{n-1}{d}_3& c_{n-1}{d}_2& c_{n-1}{d}_1& c_{n-1}{d}_0\\ .& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ c_3{d}_n& c_3{d}_{n-1}& ...& c_3{d}_3& c_3{d}_2& c_3{d}_1& c_3{d}_0\\ c_2{d}_n& c_2{d}_{n-1}& ...& c_2{d}_3& c_2{d}_2& c_2{d}_1& c_2{d}_0\\ c_1{d}_n& c_1{d}_{n-1}& ...& c_1{d}_3& c_1{d}_2& c_1{d}_1& c_1{d}_0\\ c_0{d}_n& c_0{d}_{n-1}& ...& c_0{d}_3& c_0{d}_2& c_0{d}_1& c_0{d}_0\end{array}\right]---\left(38\right)$

该系数矩阵的产生原则由最小二乘法确定。

各未知系数确定后，随光学参数变化所对应的频域信息曲面便可用已建立的二元多项式函数关系替代。将μ_a，μ_s直接代入上述数学函数中，根据矩阵乘法，任意光学参数点下的频域信息即可由下式快速计算而得：

A＝[μ_a ⁿ...μ_a ³ μ_a ² μ_a ¹ μ_a ⁰]×p_A×[μ_s ⁿ...μ_s ³ μ_s ² μ_s ¹ μ_s ⁰]    (39)

φ＝[μ_a ⁿ...μ_a ³ μ_a ² μ_a ¹ μ_a ⁰]×p_φ×[μ_s ⁿ...μ_s ³ μ_s ² μ_s ¹ μ_s ⁰]    (40)

由于使用不同阶数的多项式拟合精度不同，因此需要对阶数的选取进行研究，一般情况下，拟合精度随阶数增加而增加，但并不是阶数越高越好，相反在高于某个阶数时，会产生过拟合现象，造成更大的拟合误差。本发明选用阶数n＝5的二元多项式来对频域曲面进行回归拟合，以达到较高的拟合精度。

发明人发现，在反演过程中，幅值信息的影响同样大于相位信息的影响，所以针对最优化问题中最小化误差的表达式，可以将误差ε的求取公式，从表达式(41)改动为表达式(42)，以削弱相位的信息在重构过程中的作用，增强幅值信息在重构过程中的作用。从而进一步的提高了重建精度。

ε＝min(A_待测-A_理论)²+(φ_待测-φ_理论)²    (41)

ε＝min[(A_待测-A_理论)²+α(φ_待测-φ_理论)²]    (42)

式(42)中的α为0到1之间的权值常数，本发明中取α＝0.8。

本发明公开和揭示的所有组合和方法可通过借鉴本发明公开内容产生，尽管本发明的组合和方法已通过详细实施过程进行了描述，但是本领域技术人员明显能在不脱离本发明内容、精神和范围内对本发明所述的方法进行改动，更具体地说，所有相类似的改动对本领域技术人员来说是显而易见的，他们都被视为包括在本发明精神、范围和内容之中。

Claims (4)

1.一种基于频域近红外光测量的光学参数重构方法，其特征在于，包括下列步骤：

（1）、针对待检测生物组织的形状建立模型，进行时域正向时域蒙特卡洛模拟，获得组织的时域信息；

（2）、采用如下的方法从正向时域蒙特卡洛模拟结果中提取幅值和相位信息：已知每个光子包的权重W_d、平均飞行时间t_d，采用

计算某一频率ω₀下的傅立叶变换，再应用式

和式

分别得到出射光的幅值A和相位Φ，式中，G_R，G_I分别为

的实部和虚部，建立频域信息数据库；

（3）对已建立的数据库中的频域信息进行聚类分析，选取出射光的幅值A和相位Φ以及组织体光学参量吸收系数μ_a和散射系数μ_s为聚类过程中的4个变量，采用最短距离法的原则对标准化后的变量进行聚类，聚类完成后，分别对每一大类求取其各自的幅值的平均值和相位的平均值

吸收系数μ_a的平均值

和散射系数μ_s的平均值

将所得

及

作为类的中心点，并按照距离公式计算待测点距各类中心的距离，获取满足d最小的那一类，该满足d最小的类的类中心的

即为待测点在重构过程中的初始光学参数值，在距离公式中，

与

分别代表待测点的幅值和相位信息，n是类的个数，δ_A和

分别为数据库中幅值和相位信息的标准差，α为0到1之间的常数；

（4）以步骤（3）得到的初始光学参数值作为试探值，采用L-M最优化方法，进行光学参数重构，先由快速MC模型得到初始光学参数所对应的模拟频域信息，计算模拟频域信息与待测点频域信息之间的误差，根据误差修改每一次迭代过程中的光学参数试探值，再由快速MC模型得到相应的模拟频域信息，并重复进行上述步骤，直至误差达到最小。

2.根据权利要求1所述的基于频域近红外光测量的光学参数重构方法，其特征在于，步骤（4）中采用的快速MC模型为：根据已建立的频域信息数据库的数据，运用二元多项式模型进行曲线拟合，得到从(μ_a,μ_s)到(A，φ)的直接映射关系，从而可利用拟合函数快速获得L-M算法在重建过程中需要的任意光学参数下的模拟频域信息。

3.根据权利要求2所述的基于频域近红外光测量的光学参数重构方法，其特征在于，二元多项式曲线拟合的阶数为5。

4.根据权利要求1所述的基于频域近红外光测量的光学参数重构方法，其特征在于，步骤（4）计算误差时，通过加权的方法，削弱相位信息在重构过程中的作用。

Images