CN107045728B - 生物发光断层成像复合正则化重建的自适应参数选择方法 - Google Patents

Abstract

生物发光断层成像复合正则化重建的自适应参数选择方法，属于医学图像处理领域。由于生物组织具有强散射，低吸收的特点，生物发光断层成像的重建问题在数学上是一个高度病态的问题，外界微小的测量扰动都会给重建结果带来很大的变化。为了降低BLT重建问题的病态性，在重建时可以使用正则化的求解方法将光源重建问题转变成一个非线性的最优化问题。为了能够取得更加理想的BLT重建效果，使用基于L1范数与TV范数的复合正则化方法对重建问题求解；本发明结合偏差原理，使用迭代的方式计算复合正则化参数。通过本方法，针对BLT的复合正则化重建方法可以自适应地计算得到合适的正则化参数。

Description

生物发光断层成像复合正则化重建的自适应参数选择方法

技术领域

本发明属于医学图像处理领域，涉及一种针对基于L1范数与TV范数的复合正则化生物发光断层成像重建方法的自适应正则化参数选择方法。

背景技术

随着计算机技术的快速发展及人类在生命科学领域研究的不断深入，各种医学影像技术和医学成像设备也进入了快速发展的时期。从传统的计算机断层成像(ComputedTomography，CT)、超生成像(Ultrasound Imaging，USI)、核磁共振成像(MagneticResonance Imaging，MRI)、正电子成像(Positron Emission Tomography，PET)、单光子发射断层成像(Single-Photon Emission Computed Tomography，SPECT)发展到激发荧光断层成像(Fluorescent Molecular Tomography，FMT)和生物发光断层成像(Bioluminescence Tomography，BLT)等多种成像技术，这些成像技术的发展也极大地推动了生命科学的进步。

生物发光断层成像(以下简称BLT)作为一种重要的成像技术，近年来已经成为光学分子影像的一个重要分支。BLT凭借其高灵敏度、背景噪声低等特性，通过在生物体表获得荧光信号从而重建生物体内荧光源的分布以及观测细胞分子水平变化。BLT技术不需要外在光源的激发，而是通过一种生物化学发光反应在体内发光。体内产生的荧光在生物组织内部以某种规律传播并不断地与生物组织发生相互作用，并到达体表。最后，在生物组织体表利用高灵敏度的探测器获得的荧光图像就可以重建出荧光光源在小动物体内的分布情况，从而在本质上揭示在体分子的活动规律。

由于生物组织具有强散射、弱吸收的特点，荧光光子无法在生物组织中沿直线传输，而是经历了大量的散射过程，导致BLT逆问题在数学上是一个高度病态的问题，外界微小的测量扰动都会给重建结果带来很大的变化。因此，降低BLT的病态性，如何唯一、准确地重建荧光光源是BLT研究的重点及热点。

近年来，国内外研究学者为了提高BLT成像的准确性和可靠性做了各种尝试。为了降低BLT重建问题的病态性，在重建时可以使用正则化的求解方法将光源重建问题转变成一个非线性的最优化问题。正则化方法最早在上世纪60年代由苏联科学院院士Tikhonov提出，用来求解形如Fx＝y方程的不适定问题。正则化的核心思想是在求解问题过程中引入先验信息，从而得到对原问题有意义的解。正则化函数构造主要针对正则化项的形式进行研究，目前有L2范数、总变分(total variation，TV)范数和L_p范数的正则化项形式。基于稀疏特性的正则化方法(典型的的是L1正则化)由于融入了光源的稀疏先验信息，能够提高重建图像的质量。但是该方法会导致重建的光源过于稀疏，进而降低成像的质量。TV正则化方法由于突出了光源的边界信息从而改善了成像质量，但该方法会消除重建图像中一些小的特征信息和小的物体信息(如稀疏光源)。使用联合L1正则化和TV正则化的复合正则化方法，可以融合两种正则化方法的优点，从而突破单一正则化方法的局限性，进而提高成像的质量。

在使用联合L1正则化和TV正则化的复合正则化方法进行BLT重建时，复合正则化的参数选取决定了成像效果的好坏。但目前对于复合正则化方法的参数选择问题，还没有一个完善的选择策略，只能通过人工地反复尝试以选取适合特定模型的复合正则化参数。这使得复合正则化参数选取过程过于繁琐，需要通过多次的重复实验才能得到合适的复合正则化参数。为能够自适应地选取合适的复合正则化参数，可以使用结合偏差原理(discrepancy principle)的迭代参数更新方法。使用该方法能够极大地提升参数选择的效率。

本发明提出一种针对基于L1范数与TV范数的复合正则化BLT重建方法的自适应正则化参数选择方法。该方法通过结合偏差原理，使用迭代的方式计算出合适的复合正则化参数。

发明内容

由于生物组织具有强散射，低吸收的特点，生物发光断层成像(bioluminescencetomography,BLT)的重建问题在数学上是一个高度病态的问题，外界微小的测量扰动都会给重建结果带来很大的变化。为了降低BLT重建问题的病态性，在重建时可以使用正则化的求解方法将光源重建问题转变成一个非线性的最优化问题。

为了能够取得更加理想的BLT重建效果，使用基于L1范数与TV范数的复合正则化方法对重建问题求解，其表达式如下：

在上式中，求解生物发光的光源分布x应使BLT正则化方程中的目标函数f(x)取最小值。

代表数据拟合项。A为包含了生物组织的组织结构与光学参数信息的系数矩阵。在BLT光源重建的过程中是根据得到的表面光子通量流率去寻找光源的分布，y为边界得到的光子通量流率测量值。α||x||_TV与β||x||₁分别代表BLT重建正则化方程的TV正则化项与L1正则化项，其中α与β分别为TV正则化项和L1正则化项的正则化参数。

在有限元框架下，||x||_TV表示为形如||x||_TV＝||Bx||₁的形式。此时，式(1)被改写为如下的形式：

正则化参数α与β对BLT重建结果往往具有较大的影响，若正则化参数选取不当，会造成求得的解与问题的真实解存在明显偏差，如图2所示。

BLT重建的复合正则化参数选择只能通过人工反复实验确定，为提升参数的选择效率，本方法提出一种针对TV与L1复合正则化的自适应参数选择方法，其描述如下：

对于式(2)所示的BLT重建复合正则化方法，为了选择合适的正则化参数α与β，应使其满足复合正则化参数偏差原理(multi-parameter discrepancy principle)，如式(3)所示。

||Ax^δ(α,β)-y_δ||₂＝cδ (3)

其中x^δ(α,β)表示在正则化参数取值为(α,β)时所求得的BLT重建光源结果，y_δ为上文提到的光子通量流率边界测量值。cδ用于表示偏差项，c表示调节系数，取值范围为：c＞1；δ表示偏差值，取值范围为：0＜δ＜1。

定义：

对式(4)分别求α与β的偏导，那么得到：

因此，式(3)所示的BLT重建复合正则化参数偏差原理表示为如下的形式：

F(α,β)-αF_α′(α,β)-βF_β′(α,β)＝c²δ² (7)

基于变分原理，得到如下等式：

此时，

如果将式(9)中的

近似地表示为T||x^δ||₁的形式，其中T是一个正数常量，那么结合式(5)与式(6)得：

F(α,β)≈||y_δ||²-αF_α′(α,β)-(β+T)F_β′(α,β) (10)

使用模型函数m(α,β)表示F(α,β)，将式(10)改写为以下的形式：

那么m(α,β)表示为如下的形式：

其中，C、D与T为迭代过程中的第一中间变量、第二中间变量、第三中间变量。

此时，得到：

其中，α_k与β_k为第k次迭代时得到的正则化参数值。由此得到第k次迭代时第一中间变量C_k、第二中间变量D_k、第三中间变量T_k的值，如下式所示：

此时结合式(11)得到k+1次迭代时，α_k+1的表达式：

结合启发式算法(heuristic algorithm)，得到β_k+1存在如下关系：

选取收缩系数(contraction factor)ω(0＜ω＜1)，则式(16)表示为如下形式：

综上所述，本方法的流程如表1所示。

表1算法流程

附图说明

图1为仿体初始的光源分布示意图，图中面积较大的黑色圆形区域为仿体区域，黑色区域中面积较小的两个白色区域为光源分布区域；

图2为不同正则化参数下的BLT重建结果；

图3迭代过程中α、β与||Ax-y||₂的结果变化曲线；

图4为使用本方法计算得到的正则化参数用于计算时得到的光源重建结果。

具体实施方式

下面根据具体实施示例与附图对本发明进行说明。

首先，通过matlab的nirfast建立仿体。本实验使用圆形区域仿体，其包含两个光源区域，如图1所示。

为了避免逆行为(inverse crime)，前向仿真的有限元网格和重建的网格节点数不同。一般情况下，为了保证重建时计算量不会过大，重建模型的有限元节点数往往小于前向的有限元节点数。在本实验中，前向网格的有限元节点数设置为3508个，共包含6807个面元；重建网格的有限元节点数设置为1309个，共包含2491个面元。

由于BLT逆问题求解中基于单光谱的重建结果可能存在不唯一性，因此实验中采用了两个不同的波段。不同光波段在仿体中拥有不同的光吸收系数与光散射系数，实验选择的两个谱段分别为600nm和630nm，此时仿体中的光学特性参数如表2所示。

表2仿体的光学特性参数

分别对不同光谱下添加光源的仿体区域进行前向仿真，合并后得到边界测量值y，并可以计算得到重建需要的系数矩阵A。

使用表1所示的算法流程进行实验。输入变量中y_δ、A与B基于实验模型求解得到，y_δ的值为边界测量值y，A和B为根据实验模型求解得到的矩阵。ε的值设定为10^-6，c与δ的值分别取2和0.1。实验中α与β的初值分别设为0.1与1。

本发明结合偏差原理，使用迭代的方式计算复合正则化参数。通过本方法计算，α、β与||Ax-y||₂的变化曲线如图3所示。最终得到的光源重建结果如图4所示，可以看出使用本方法计算得到的正则化参数用于计算时能够得到有较好的BLT重建结果。实验结果表明，本发明针对BLT的复合正则化重建方法可以自适应地计算得到合适的正则化参数。

Claims (1)

1.生物发光断层成像复合正则化重建的自适应参数选择方法，其特征在于：

使用基于L1范数与TV范数的复合正则化方法对重建问题求解，其表达式如下：

在上式中，求解生物发光的光源分布x应使BLT正则化方程中的目标函数f(x)取最小值；

代表数据拟合项；A为包含了生物组织的组织结构与光学参数信息的系数矩阵；在BLT光源重建的过程中是根据得到的表面光子通量流率去寻找光源的分布，y为边界得到的光子通量流率测量值；α||x||_TV与β||x||₁分别代表BLT重建正则化方程的TV正则化项与L1正则化项，其中α与β分别为TV正则化项和L1正则化项的正则化参数；

在有限元框架下，||x||_TV表示为形如||x||_TV＝||Bx||₁的形式；此时，式(1)被改写为如下的形式：

对于式(2)所示的BLT重建复合正则化方法，为了选择合适的正则化参数α与β，应使其满足复合正则化参数偏差原理，如式(3)所示；

||Ax^δ(α,β)-y_δ||₂＝cδ (3)

其中x^δ(α,β)表示在正则化参数取值为(α,β)时所求得的BLT重建光源结果，y_δ为上文提到的光子通量流率边界测量值；cδ用于表示偏差项，c表示调节系数，取值范围为：c＞1；δ表示偏差值，取值范围为：0＜δ＜1；

定义：

对式(4)分别求α与β的偏导，那么得到：

因此，式(3)所示的BLT重建复合正则化参数偏差原理表示为如下的形式：

F(α,β)-αF′_α(α,β)-βF′_β(α,β)＝c²δ² (7)

基于变分原理：

此时，

如果将式(9)中的

近似地表示为T||x^δ||₁的形式，其中T是一个正数常量，那么结合式(5)与式(6)得：

F(α,β)≈||y_δ||²-αF′_α(α,β)-(β+T)F′_β(α,β) (10)

使用模型函数m(α,β)表示F(α,β)，将式(10)改写为以下的形式：

那么m(α,β)表示为如下的形式：

其中，C、D与T为迭代过程中的第一中间变量、第二中间变量、第三中间变量；

此时，得到：

其中，α_k与β_k为第k次迭代时得到的正则化参数值；由此得到第k次迭代时第一中间变量C_k、第二中间变量D_k、第三中间变量T_k的值，如下式所示：

此时结合式(11)得到k+1次迭代时，α_k+1的表达式：

结合启发式算法(heuristic algorithm)，得到β_k+1存在如下关系：

选取收缩系数(contraction factor)ω(0＜ω＜1)，则式(16)表示为如下形式：

Images