基于迭代测量的光学断层重建方法与流程

本发明属于生物光学成像技术领域，具体地涉及一种基于迭代测量的光学断层重建方法。

背景技术

光学断层成像是生物医疗研究中一种非常重要的成像方式，是一种高灵敏度的无损伤检测技术。利用该技术可以对生物组织各部分的吸收系数、散射系数等光学参数进行成像，从而提供生物组织体的重要功能信息，具有广阔的发展前景和研究价值。其显著特点包括：无损伤性、便携性、连续性、功能性测量。目前，该技术已成功应用于血氧检测，脑出血的早期诊断与定位、检测乳腺癌等领域。

光学断层图像重建包括正向问题和逆向问题。正向问题即根据一定的光子传输模型和光学参数得出组织边界测量的理论估计值的过程。根据一组边界测量值和特定的光子传输模型重建表征组织体内在特性的光学参数的空间分布称为逆向问题，或图像重建问题。

由于在重建过程中需要重复计算正向问题，因此，正向模型成为了决定重建效率的关键因素之一。众所周之，可以利用辐射传输方程对光在生物组织中传输过程进行建模。然而，辐射传输方程计算量非常大。因此，可采用扩散近似作为正向模型。但是扩散近似具有内在的局限性：首先，对于介质是低散射或者非散射的区域，扩散方程不适用。其次，扩散近似无法准确地对平行光源附近的光传输进行建模。鉴于以上局限性，扩散近似只有在高散射区域才能提供相对精确的传输模型。此外，在重建时，需要利用测量数据来恢复光学参数的分布。因此，探测器的数量和位置会对重建结果产生影响。实际上，可以通过增加边界探测器数量来提高重建精度。然而，这会导致矩阵规模增大，从而增加重建的计算代价。利用简化模型的方法，采用降维的子空间来表示未知数，此方法需要利用基向量来构建变换矩阵，增加了重建过程的计算负担。此外，利用数据子集的方法，可以提高成像速度，但在目标区域的成像误差较大。光学断层图像重建的两个重要的性能指标是重建速度和重建精度。为了提高重建效率，需要寻找一个有效的解决方法。

技术实现要素：

针对上述技术问题，本发明目的是：提供一种基于迭代测量的光学断层重建方法，可以弥补由扩散近似造成的正向模型误差，从而提高对光传输过程的建模精度。与基于辐射传输方程的重建方法相比，该方法基于扩散近似进行计算，因此进一步降低了计算代价。此外，可以减少灵敏度矩阵的计算量，从而加速整个重建过程，提高重建效率。

本发明的技术方案是：

一种基于迭代测量的光学断层重建方法，包括以下步骤：

s01：建立正向模型误差；

s02：以经过成像中心的直线将均匀分布于边界的测量数据划分为数量相等的两组测量数据，即ya、yb；

s03：利用第一组测量数据计算迭代更新矢量δxa，根据δxa对待重建的吸收系数x进行迭代计算，即x＝x+δxa，判断建立的目标函数值，若目标函数值小于预设的阈值η，则结束迭代计算，得到最终的待重建的吸收系数x，反之进入步骤s04；

s04：利用第二组测量数据计算迭代更新矢量δxb，根据δxb对待重建的吸收系数x进行迭代计算，即x＝x+δxb，判断建立的目标函数值，若目标函数值小于预设的阈值η，则结束迭代计算，得到最终的待重建的吸收系数x，反之进入步骤s03。

优选的技术方案中，所述步骤s01中，利用高斯模型表示正向模型误差ε，其中，n表示高斯分布，γε和分别表示正向模型误差ε的协方差和均值，利用高斯模型表示待重建的吸收系数x，其中，γx和分别表示x的协方差和均值。

优选的技术方案中，所述步骤s01还包括：

s11：利用高斯模型产生样本x^l(l＝1,2,…,r)，其中，r表示样本数量，利用辐射传输方程以及获得样本的精确光子密度φ，其中，φ表示辐射量，在有限元框架下，辐射传输方程为矩阵方程其中，q表示光源，表示辐射向量，ar为刚度矩阵；

s12：利用扩散近似获得样本的扩散近似光子密度φd，在有限元框架下，扩散近似为矩阵方程adφd＝q，其中，q表示光源，ad为刚度矩阵；

s13：根据ε^l＝φ(x^l)-φd(x^l)计算正向模型误差，根据和分别计算正向模型误差的协方差和均值。

优选的技术方案中，所述步骤s03中，其中，矩阵和表示cholesky因子，j为灵敏度矩阵。

优选的技术方案中，所述步骤s03中建立的目标函数ψ(x)＝||ya-φd(x)||。

优选的技术方案中，所述步骤s04中，其中，矩阵和表示cholesky因子，j为灵敏度矩阵。

优选的技术方案中，所述步骤s04中建立的目标函数ψ(x)＝||yb-φd(x)||。

与现有技术相比，本发明的优点是：

将辐射传输方程与扩散近似之间的正向模型误差融入到成像过程中。该方法的特点是可以弥补由扩散近似造成的正向模型误差，从而提高对光传输过程的建模精度。此外，与基于辐射传输方程的重建方法相比，该方法基于扩散近似进行计算，进一步降低了计算代价。在重建时，对边界测量数据进行分组，在此基础上进行图像重建。该方法可以获得很高的重建精度和重建速度。

附图说明

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

图1为本发明基于迭代测量的光学断层重建方法的流程图；

图2为单目标体仿真模型；

图3为基于扩散近似的重建结果与本发明重建结果的单目标体重建结果对比图；

图4为双目标体仿真模型；

图5为基于扩散近似的重建结果与本发明重建结果的双目标体重建结果对比图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。

实施例：

如图1所示，基于迭代测量的光学断层重建方法，包括：

正向模型误差的融入，用于弥补由扩散近似造成的正向模型误差，从而提高对光传输过程的建模精度；

测量数据的分组迭代，用于减少重建计算量，提高重建速度。

步骤1、正向模型误差的融入，包括：

1)利用高斯模型表示正向模型误差ε，其中，n表示高斯分布，γε和分别表示正向模型误差ε的协方差和均值。利用高斯模型表示待重建的吸收系数x，其中，γx和分别表示x的协方差和均值。

2)根据高斯模型产生一组样本x^l(l＝1,2,…,r)，其中，r表示样本数量，利用辐射传输方程对上述样本进行正向问题的计算，即利用辐射传输方程以及获得样本的精确光子密度φ，其中，φ表示辐射量。在有限元框架下，辐射传输方程可表示为矩阵方程其中，q表示光源，表示辐射向量，根据公式

计算刚度矩阵ar，其中，ar表示单元刚度矩阵中第i行第j列的元素，ω表示角频率，ω是成像区域，是成像区域的边界，μs表示散射系数，μa表示吸收系数，c表示光速，u＝[u1,u2,…,um]^t为空间节点的基函数，m为节点基函数的数量，v＝[v1,v2,…,vp]^t为角度基函数，p为角度基函数的数量，表示边界上的单位法向量，θ′和θ分别表示光子的入射角度和散射角度，(cosθ,sinθ)表示单位圆上的点，

3)利用扩散近似对高斯模型所产生的样本x^l(l＝1,2,…,r)进行正向问题的计算。即利用扩散近似获得扩散近似光子密度φd。在有限元框架下，扩散近似可以表示为矩阵方程adφd＝q，其中，q表示光源。根据计算刚度矩阵ad，其中，ad表示单元刚度矩阵中第i行第j列的元素，▽表示梯度，表示扩散系数，g为各向异性因子。

4)根据ε^l＝φ(x^l)-φd(x^l)计算正向模型误差；根据和分别计算正向模型误差的协方差和均值。

步骤2、测量数据的分组迭代，将边界测量数据分组，轮流用于迭代重建，直到目标函数值小于预设的阈值，从而减少重建计算量，提高重建速度。测量数据的分组迭代按照如下步骤进行：

1)将待重建的吸收系数x初始化，即以经过成像中心的直线将均匀分布于边界的测量数据y划分为数量相等的两组测量数据，即ya和yb。

2)利用第一组测量数据ya计算迭代更新矢量

δxa用于迭代更新重建未知量x，即x＝x+δxa，对矩阵和进行cholesky分解，即矩阵和表示cholesky因子。按照计算灵敏度矩阵j；按照ψ(x)＝||ya-φd(x)||计算目标函数ψ。若目标函数值小于预设的阈值η，则结束迭代计算，获得最终的重建结果x；否则进入第3)步，即利用第二组测量数据yb再进行迭代计算。

3)利用第二组测量数据yb计算迭代更新矢量

δxb用于迭代更新重建未知量x，即x＝x+δxb，其中，按照计算灵敏度矩阵j；按照ψ(x)＝||yb-φd(x)||计算目标函数ψ。若目标函数值小于预设的阈值η，则结束迭代计算，获得最终的重建结果x；否则回到第2)步，即利用第一组测量数据ya再进行迭代计算。

使用实施例公开的基于迭代测量的光学断层重建方法分别重建两个模型，并将其性能与基于扩散近似的重建方法进行比较。

首先对包含单个目标体的仿真模型进行重建，仿真模型如图2所示。采用30个探测器和4个光源均匀分布于模型四周，有限元离散网格包含212个单元和122个节点。目标区域μa为0.3mm^-1，μs为0.03mm^-1。背景区域μa为0.06mm^-1，μs为0.03mm^-1。g为0.1，η为0.001。图3(a)给出了基于扩散近似的重建结果，图3(b)给出了基于本发明的重建结果。为了进一步对重建结果进行定量评估，引入均方误差，计算如下：其中n表示成像区域的节点数量，和表示重建结果和真实值。表1列出了不同重建方法的性能指标。

表1

接下来，对包含双目标体的仿真模型进行重建，如图4所示。采用30个探测器和4个光源均匀分布在重建模型四周。有限元离散网格包含264个单元和148个节点。目标区域μa分别为0.2mm^-1和0.3mm^-1，μs为0.03mm^-1。背景区域μa为0.06mm^-1，μs为0.03mm^-1。g为0.1，η为0.001。图5(a)给出了基于扩散近似的重建结果，图5(b)给出了基于本发明的重建结果。表2列出了不同方法的重建性能。

表2

由实验结果可知：

本发明公开的基于迭代测量的光学断层重建方法可以更好地重建出目标体的结构，有效提高图像重建的速度和精度。

应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。