CN107713995B - 一种基于惩罚算法的生物发光断层成像光源重建方法 - Google Patents

Abstract

本发明公开了一种基于惩罚算法的生物发光断层成像光源重建方法，涉及光学分子影像技术领域，采用了保留硬约束并惩罚软约束的方法，采用的方法能够求解非凸的目标函数，运用适当的正则化参数更新方案，自适应得到最优的正则化参数，减小了算法对所设定初始正则化参数的敏感性，运用了非凸l_p范数，使得到的解更加稀疏，对光源位置的反演更加准确。

Description

一种基于惩罚算法的生物发光断层成像光源重建方法

技术领域

本发明涉及光学分子影像技术领域，特别是涉及一种基于惩罚算法的生物发光断层成像光源重建方法。

背景技术

生物发光断层成像(Bioluminescence tomography，BLT)是分子影像学的一个重要成像模态，具体过程是首先通过转基因技术将荧光素酶基因转录到靶细胞的DNA上以表达荧光素酶，再向生物体内注射底物荧光素，在生物体内环境中三磷酸腺苷、氧气及镁离子存在的情况下，荧光素酶会催化荧光素底物发生氧化反应，发射出近红外谱段的光子。产生的光子在生物体内经历吸收、散射等行为，其中有一部分光子逸出生物表面被高灵敏设备捕获，通过对于表面荧光强度数据的测量，能够定性，定量的反应生物体细胞的生理变化。

光在生物组织中的传播可由辐射传输方程(Radiative transport equation，RTE)进行准确描述，但由于辐射传输方程求解较为复杂，一般采用其一阶球谐近似模型，即扩散方程(Diffusion equation，DE)描述光在生物组织中的传播。求解扩散方程有三种方法：解析法、统计法和数值法。通常情况下，任意模型上的扩散方程的解析解很难得到，因此一般采用数值法或统计法。统计法中最常用的是蒙特卡洛(Monto Carlo，MC)方法，由于其本质上是一种随机统计方法，计算运行时间长。所以，目前最常用的求解扩散方程的方法是数值法中的有限元方法。

生物发光断层成像的光源重建问题是典型的逆问题。由于逆问题的病态性，以及CCD相机能探测到的生物发光信号极其有限，所以需要利用正则化方法才能获得近似解。这里常用的正则化算法包括：基于L₁范数的重建方法(如不完全变量截断共轭梯度算法(Incomplete variables truncated conjugate gradient，IVTCG))，基于L₂范数的重建算法(Tikhonov正则化)，基于L_p(0<p<1)范数的重建方法等。然而毫无疑问，上述正则化方法均面临着正则化参数的选择的难题，并且在光源重建的定位精度和算法的稳定性等方面仍存在不足。

发明内容

本发明实施例提供了一种基于惩罚算法的生物发光断层成像光源重建方法，可以解决现有技术中存在的问题。

本发明提供了一种基于惩罚算法的生物发光断层成像光源重建方法，所述方法包括以下步骤：

步骤100，采用CCD相机采集生物体的图像，以此获得生物体的表面荧光强度数据；

步骤200，结合步骤100中获得的表面荧光强度数据及给定的光学特异性参数，根据光在生物体中传输过程的扩散近似模型，采用有限元方法构建生物体内部光源和表面荧光强度的线性关系；

步骤300，将步骤200中生物体内部光源和表面荧光强度的线性关系转化为约束优化问题；

步骤400，对于步骤300中的约束优化问题，采用惩罚算法解出最优解x^*；

步骤500，利用MATLAB相关工具包，对成像结果进行展示。

本发明实施例中的一种基于惩罚算法的生物发光断层成像光源重建方法，具有以下效果：

1、采用了保留硬约束并惩罚软约束的方法，采用的方法能够求解非凸的目标函数，如l_p范数。

2、运用适当的正则化参数更新方案，自适应得到最优的正则化参数，减小了算法对所设定的正则化参数的敏感性。

3、运用了l_p范数，使得到的解更加稀疏。

4、对光源位置的反演更加准确。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的基于惩罚算法的生物发光断层成像光源重建方法的总体流程图；

图2(a)为包含两个相距6mm且半径均为1mm的球状光源的仿真模型及重建结果的YZ平面透视图；

图2(b)为包含两个相距6mm且半径均为1mm的球状光源的仿真模型及重建结果的XZ平面透视图；

图2(c)为包含两个相距6mm且半径均为1mm的球状光源的仿真模型及重建结果的XY平面透视图；

图3(a)为包含两个相距6mm且半径分别为1mm和2mm的球状光源的仿真模型及重建结果的YZ平面透视图；

图3(b)为包含两个相距6mm且半径分别为1mm和2mm的球状光源的仿真模型及重建结果的XZ平面透视图；

图3(c)为包含两个相距6mm且半径分别为1mm和2mm的球状光源的仿真模型及重建结果的XY平面透视图；

图4(a)为包含两个水平相距6mm，垂直相距2mm，且半径均为1mm的球状光源的仿真模型及重建结果的YZ平面透视图；

图4(b)为包含两个水平相距6mm，垂直相距2mm，且半径均为1mm的球状光源的仿真模型及重建结果的XZ平面透视图；

图4(c)为包含两个水平相距6mm，垂直相距2mm，且半径均为1mm的球状光源的仿真模型及重建结果的XY平面透视图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

参照图1，本发明提供了一种基于惩罚算法的生物发光断层成像光源重建方法，该方法包括：

步骤100，采用CCD相机采集生物体的图像，以此获得生物体的表面荧光强度数据，将包括吸收系数、散射系数和各向异性因子在内的光学特异性参数作为先验信息；

定义成像域为一个高25mm，半径20mm的均匀圆柱状仿体，成像目标内部光源发射出光子，穿透仿体并被高灵敏度的CCD相机采集到。

步骤200，结合步骤100中获得的表面荧光强度数据及给定的光学特异性参数，根据光在生物体中传输过程的扩散近似模型，采用有限元方法构建生物体内部光源和表面荧光强度的线性关系。具体方法如下：

步骤210：利用光在生物体中传播的模型以描述光在生物体中的传输过程，所采用的模型为辐射传输方程的一阶球谐近似，即扩散方程；

步骤220：利用网格化软件amira对目标离散化，将成像域剖分为多个四面体单元构成的有限元网格，得到网格化后各节点的坐标，节点和四面体之间的对应关系，节点和四面体的面之间的对应关系；

步骤230：对步骤210中的扩散方程进行离散化，结合给定的光学特异性参数构建表面荧光强度数据和内部光源分布之间的线性关系：

b＝Ax

其中，b为表面荧光强度数据，A为系统矩阵，表示表面荧光强度数据和内部光源分布之间的线性关系，x是要求解的内部光源能量分布。

步骤300，将步骤200中生物体内部光源和表面荧光强度的线性关系转化为约束优化问题。具体实现方式如下：

将步骤230所得的线性关系转化为以下约束优化问题：

s.t.||Ax-b||≤σ

其中，

是非凸目标函数，σ＝1e^-10代表噪声容限。

步骤400，对于步骤300中的约束优化问题，采用惩罚算法解出最优解x^*。具体实现方式如下：

步骤410：构建步骤300中目标函数更加简化的形式：

其中

f_λ,μ(x):＝h_λ,μ(||Ax-b||²-σ²)

步骤420：设定步骤410中目标函数的初始可行解x^feas＝A⁺b，其中A⁺＝A^T(AA^T)^-1为A的伪逆，设定初始参数λ₀＝1,μ₀＝1,ε₀＝1,ρ＝2,σ＝1e^-10,σ＝1e^-10,θ＝1/ρ,τ＝2，令x⁰＝0作为算法的起始点；

步骤430：若

则将x^k,0的值替换为x^feas的值，用非单调近端梯度(NPG)方法解出步骤410中目标函数的近似驻点，也就是函数的一阶导数与0的距离小于ε_k的点，即满足：

的步骤410中目标函数的近似驻点x^k，其中ε_k为近似程度；

步骤440：判断得出结果的残差和近似程度是否都小于设定容限，即：

max{(||Ax^k-b||²-σ²)₊,0.01ε_k}≤tol

其中tol是人工设定的终止条件，若满足则x^k就是目标函数的最优解x^*，若不满足则执行参数更新，使λ_k+1＝ρλ_k,μ_k+1＝θμ_k,ε_k+1＝θε_k,x^k+1＝x^k,k＝k+1并返回步骤430；

具体的，步骤430中应用的NPG方法为：

步骤431：选择参数L_max＝10⁸,L_min＝1,τ＝2,c＝10^-4,M＝4,k＝0；

步骤432：选择初始

步骤433：求解：

步骤434：若u的目标函数值稍小于或等于前M+1次迭代中最大的目标函数值即若：

则结束搜寻并用Barzilai-Borwein(BB)方法选择下一次内部迭代时非单调近端梯度方法的参数

并使x^k+1＝μ,k＝k+1，若大于前M+1次迭代中最大的目标函数的值，则更新L_k＝τL_k并返回步骤433；

步骤435：若x^k,l满足

的近似一阶最优条件，即：

和

得出x^k＝x^k,l，若不满足，则返回步骤433，其中x^k,l代表NPG中内迭代l次，主迭代k次得出的解，Diag(x^k,l)代表第i个对角元素是

的对角矩阵，|x^k,l|^p代表第i个元素是

的矩阵。

步骤500，利用MATLAB相关工具包，对光源重建结果进行展示。

实验结果分析

结合附图2、3、4，对本发明基于惩罚算法的生物发光断层成像光源重建方法进一步进行说明：

附图2(a)、(b)、(c)中，在半径为20mm，高为24mm的圆柱状仿体中，两光源被放置在高18mm，相距6mm的地方。星型标识点为真实光源位置，不规则部分表示重建的荧光光源。两光源都是半径为1mm的球形光源，可以看出，本方法得到的重建结果对光源中心点的定位准确。

附图3(a)、(b)、(c)中，在半径为20mm，高为24mm的圆柱状仿体中，两光源被放置在高18m，相距6mm的地方，星型标识点为真实光源位置，不规则部分表示重建的荧光光源。两光源的半径分别为1mm和2mm。

附图4(a)、(b)、(c)中，在半径为20mm，高为24mm的圆柱状仿体中，两光源其中一个被放置在高16mm的地方，另一个被放置在高18mm的地方。水平距离为6mm，星型标识点为真实光源位置，不规则部分表示重建的荧光光源。

由附图2、3、4可以看出，对于不同的光源设定，基于本发明的生物发光断层成像光源重建方法所得到的重建结果具有光源中心定位准确、光源形状近似的优点。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (2)

1.一种基于惩罚算法的生物发光断层成像光源重建方法，其特征在于，所述方法包括以下步骤：

步骤100，采用CCD相机采集生物体的图像，以此获得生物体的表面荧光强度数据；

步骤200，结合步骤100中获得的表面荧光强度数据及给定的光学特异性参数，根据光在生物体中传输过程的扩散近似模型，采用有限元方法构建生物体内部光源和表面荧光强度的线性关系；

步骤300，将步骤200中生物体内部光源和表面荧光强度的线性关系转化为约束优化问题；

步骤400，对于步骤300中的约束优化问题，采用惩罚算法解出最优解x^*；

步骤500，利用MATLAB相关工具包，对光源重建结果进行展示；

步骤200具体包括：

步骤210：利用光在生物体中传播的模型以描述光在生物体中的传输过程，所采用的模型为辐射传输方程的一阶球谐近似，即扩散方程；

步骤220：利用网格化软件amira对目标离散化，将成像域剖分为多个四面体单元构成的有限元网格，得到网格化后各节点的坐标，节点和四面体之间的对应关系，节点和四面体的面之间的对应关系；

步骤230：对步骤210中的扩散方程进行离散化，结合给定的光学特异性参数构建表面荧光强度数据和内部光源分布之间的线性关系：

b＝Ax

其中，b为表面荧光强度数据，A为系统矩阵，表示表面荧光强度数据和内部光源分布之间的线性关系，x是要求解的内部光源能量分布；

步骤300具体包括：

将步骤230所得的线性关系转化为以下约束优化问题：

s.t.||Ax-b||≤σ

其中，

是非凸目标函数，σ＝1e^-10代表噪声容限；

步骤400具体包括：

步骤410：构建步骤300中目标函数更加简化的形式：

其中

f_λ,μ(x):＝h_λ,μ(||Ax-b||²-σ²)

步骤420：设定步骤410中目标函数的初始可行解x^feas＝A⁺b，其中A⁺＝A^T(AA^T)^-1为A的伪逆，设定初始参数λ₀＝1,μ₀＝1,ε₀＝1,ρ＝2,σ＝1e^-10,σ＝1e^-10,θ＝1/ρ,τ＝2，令x⁰＝0作为算法的起始点；

步骤430：若

则将x^k,0的值替换为x^feas的值，用非单调近端梯度方法解出步骤410中目标函数的近似驻点，也就是函数的一阶导数与0的距离小于ε_k的点，即满足：

的步骤410中目标函数的近似驻点x^k，其中ε_k为近似程度；

步骤440：判断得出结果的残差和近似程度是否都小于设定容限，即：

max{(||Ax^k-b||²-σ²)₊,0.01ε_k}≤tol

其中tol是人工设定的终止条件，若满足则x^k就是目标函数的最优解x^*，若不满足则执行参数更新，使λ_k+1＝ρλ_k,μ_k+1＝θμ_k,ε_k+1＝θε_k,x^k+1＝x^k,k＝k+1并返回步骤430。

2.如权利要求1所述的基于惩罚算法的生物发光断层成像光源重建方法，其特征在于，步骤430中应用的NPG方法为：

步骤431：选择参数L_max＝10⁸,L_min＝1,τ＝2,c＝10^-4,M＝4,k＝0；

步骤432：选择初始

步骤433：求解：

步骤434：若u的目标函数值小于或等于前M+1次迭代中最大的目标函数值即若：

则结束搜寻并用Barzilai-Borwein方法选择下一次内部迭代时非单调近端梯度法的参数

并使x^k+1＝μ,k＝k+1，若大于前M+1次迭代中最大的目标函数的值，则更新L_k＝τL_k并返回步骤433；

步骤435：若x^k,l满足

的近似一阶最优条件，即：

和

得出x^k＝x^k,l，若不满足，则返回步骤433，其中x^k,l代表NPG方法中内迭代l次，主迭代k次得出的解，Diag(x^k,l)代表第i个对角元素是

的对角矩阵，|x^k,l|^p代表第i个元素是

的矩阵。

Images