CN102334979B - 一种基于迭代收缩的双模态融合断层成像方法 - Google Patents

Abstract

本发明属于医学分子影像领域，涉及自发荧光断层成像与计算机断层成像两种模态融合的方法，尤其是一种基于迭代收缩的双模态融合断层成像方法。这种技术用于定量并定位重建目标体内部的光源强度与位置，解决由目标体表面有限的光强分布反推得到内部全部的光强分布的逆向问题。本技术方案的要点在于，利用自发荧光断层成像中获得的表面光强信息，融合计算机断层成像中获得的内部几何结构信息，采用迭代收缩法将重建中复杂的多维求优过程转化为一维并行求优的高效循环过程，从整体上得到对正则参数、l_p范数、噪声、初值鲁棒的精确重建结果。该技术能有效地应用于目标体全身生理代谢的研究，重建效率高且适用于成像系统性能较低的情况。

Description

一种基于迭代收缩的双模态融合断层成像方法

技术领域

本发明属于医学分子影像领域，涉及自发荧光断层成像与计算机断层成像两种模态融合的方法，尤其是一种基于迭代收缩的双模态融合断层成像方法。

背景技术

自发荧光断层成像是近年来新兴的一种分子影像技术，是一种可以在细胞分子水平上实现生物组织内部生物物理过程无创实时动态连续成像的一种技术。通过获取到的目标体表面出射光强分布，探测在目标体内部的光源分布则能够对目标体进行定量化研究，为荧光光源提供三维的定位信息。

目前普遍采用的断层成像技术大多是基于梯度的传统方法，由于自发荧光断层成像的信号十分微弱，导致了采集的数据集噪声较大，加之该问题的固有病态性，造成在待探测过程中对正则参数的选择、l_p范数的选择、初始值的选择、可行域的选择以及噪声的容忍度都有着很强的要求。这限制了自发荧光断层成像在实际中的应用。寻找一种高效、鲁棒、全域的探测方法一直是该领域的热点、难点问题，到目前为止仍然没有一种通用的技术能够很好地解决这个问题。

迭代收缩法是一种新兴的数值计算方法，与传统的迭代优化方法相比，它更具高效性，且近年来已被陆续应用于图像降噪以及逆问题求解中，推动了这些领域的发展。

发明内容

本发明的目的是提供一种基于迭代收缩的双模态融合断层成像技术，利用高性能探测器捕捉目标体在自发荧光断层成像中出射的表面光强分布，融合计算机断层成像中获得的目标体内部几何结构信息，使用扩散方程建立光子传输模型，采用迭代收缩法将重建中复杂的多维求优过程转化为对每个变量逐次进行一维并行求优的高效循环过程，从整体上得到对正则参数、l_p范数、噪声、初值鲁棒的精确重建结果。

为了实现上述目的，本发明涉及的基于迭代收缩的双模态融合断层成像方法具体步骤如下：

步骤a，采用波尔兹曼传输方程的一阶球谐近似得到的扩散方程建立光子在目标体内部的传输模型，即扩散方程传输模型；

步骤b，利用有限元方法将扩散方程传输模型转换成线性矩阵方程的形式；

步骤c，基于线性矩阵方程形式的扩散方程传输模型，利用正则理论，建立优化目标函数f(S)，其中，S为待求解的目标体内部光强分布；

步骤d，对f(S)进行降维处理，得到一维优化目标函数；

步骤e，基于所述一维优化目标函数，以迭代的方式求解得到最优目标体内部光强分布S。

本发明的有益效果是：通过使用基于迭代收缩的双模态融合断层成像技术，融合自发荧光断层成像与计算机断层成像两种模态的信息降低逆问题求解的病态性，由并行计算迭代收缩算子的复杂的多维求优过程转化为对每个变量逐次进行一维并行求优的高效循环过程，最终重建得到目标体内部光源分布的定量和定位信息。使用该方法的优势为，在不同正则参数、不同l_p范数、不同初值、较强噪声以及全部目标区域下均能计算得到精确的重建结果；同时，通过迭代算子的降维并行处理，使这种双模态融合的断层成像技术具有高于传统方法的重建效率。

附图说明

图1为基于迭代收缩的双模态融合断层成像技术框架流程图。

图2为在体非匀质小鼠的自发荧光断层成像实验的表面光强图。

图3为在体非匀质小鼠的计算机断层成像实验的解剖结构图。

图4为平滑算子θ(w)的参数w与其拟合对象的参数p的关系图。

图5为使用本技术在不同正则参数下的重建结果图。

图6为使用本技术在不同l_p范数下的重建结果图。

图7为目标优化函数值与迭代次数的演化曲线图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。应指出的是，所描述的实例仅旨在便于对本发明的理解，而对其不起任何限定作用。

本发明是一种基于迭代收缩的双模态融合断层成像技术。在自发荧光断层成像中，例如通过对小动物转基因，使得小动物在病灶处产生荧光蛋白酶，进而通过静脉注射荧光底物，从而在该蛋白酶的催化作用下发射出可见或近红外光，光在生物组织中的传输将历经多次散射。最终有一部分光子到达小动物皮肤并出射到自由空间中被探测器接收，作为自发荧光断层成像中获取到的表面光强分布信息。再融合通过计算机断层成像获取到的小动物内部解剖结构信息，借助不同组织器官的光学特性参数和分割、配准、网格剖分技术建立非匀质小动物模型，最终重建出小动物体内的光源强度与位置信息。本发明能在不同正则参数、不同l_p范数、不同初值、较强的噪声的情况下，精确重建出小动物全身的内部光源信息。

下面结合附图，详细描述本发明的断层图像重建方法。作为一种具体的实现技术方案，其框架流程图参阅图1，步骤如下：

步骤1：采用波尔兹曼传输方程的一阶球谐近似得到的扩散方程建立光子在目标体内部的传输模型，即扩散方程传输模型：

$-\▿\·[D\left(r\right)\▿\mathrm{\Φ}\left(r\right)]+{\mathrm{\μ}}_a\left(r\right)\mathrm{\Φ}\left(r\right)=S\left(r\right)$

式中，r是位置向量，D是散射系数，μ_a是吸收系数，Φ是目标体表面光强分布，S表示待求解的目标体内部光强分布。

步骤2：利用自发荧光断层成像获取目标体表面出射的光强分布测量值Φ^*；利用计算机断层成像获取目标体内部几何结构信息，通过分割出目标体内部不同的器官结构，并为不同器官设置相应的光学参数，包括吸收系数和散射系数，从而建立非匀质目标体模型。

步骤3：对非匀质目标体模型进行几何结构的网格剖分，并利用配准技术将表面光强信息与几何结构信息进行位置融合与能量融合。请参阅图2和图3为非匀质小鼠断层成像实验获取的信息。图2给出在体非匀质小鼠表面光强信息，该数据采集于自发荧光断层成像实验，并经过配准等后处理操作。图3给出在体非匀质小鼠的解剖结构信息，该数据采集于计算机断层成像实验，并经过分割等后处理操作。融合图2和图3中的数据，并作网格剖分处理，建立起非匀质小鼠的模型，包含了躯干、心脏、肺、肝脏、肌肉、骨骼等组织器官。

步骤4：基于步骤3获得的信息，比如离散化几何信息和能量信息，利用有限元方法将扩散方程传输模型转换成线性矩阵方程MS＝Φ的形式，其中，M是描述模型性质的系统矩阵。

步骤5：为了以优化的方式迭代求解目标体内部光强分布S，要基于线性矩阵方程形式的扩散方程传输模型，利用正则理论，建立优化目标函数：

$f\left(S\right)=\frac{1}{2}{\left|\right|\mathrm{MS}-{\mathrm{\Φ}}^{*}\left|\right|}_2^2+{\mathrm{\λ}\left|\right|S\left|\right|}_p^p$

式中，Φ^*为目标体表面光强分布的测量值，λ为正则化参数，

表示l_p范数。

步骤6：若直接对f(S)优化则为多维优化问题，求解十分复杂，故对其进行降维处理，得到一维优化目标函数：

$g\left(\mathrm{\α}\right)={\frac{1}{2}\left|\right|M_{s_k}-{\mathrm{\Φ}}^{*}-m_j{s}_k\left(j\right)+m_j\mathrm{\α}\left|\right|}_2^2+{\mathrm{\λ}\left|\right|\mathrm{\α}\left|\right|}_p^p$

式中，s_k是对S的第k次估计，m_j是M的第j列元素，s_k(j)表示s_k的第j个元素，α表示s_k+1(j)的值。这样，便能实现以一维求优的形式，结合后续步骤逐一优化s_k+1中每个元素的值，最终得到最优解S。

步骤7：为了使上式g(α)在优化中存在解析解，采用平滑算子θ(w)来拟合替代

且θ(w)＝|α|-wln(1+|α|/w)，w是平滑算子参数。请参阅图4为平滑算子θ(w)的参数w与其拟合对象

的参数p的关系曲线。每一个p值都对应一个w值，在保持这种对应的关系下，θ(w)能很好地拟合

使求解具有解析形式。

步骤8：设置初始迭代次数k＝0，并为Φ^*，M，λ，w赋值，并初始化目标体内部光强分布s₀。

步骤9：为了采用迭代收缩的策略更新元素s_k+1(j)，首先计算s_k(j)的迭代收缩算子H_k(j)：

$H_k\left(j\right)=\frac{(V_j-w-U_j)+\sqrt{{\left(\right|V_j|-w-U_j)}^2+4w\left|V_j\right|}}{2\mathrm{sign}\left(V_j\right)}$

式中，

U_j＝λ/||m_j||²，sign是符号函数。在实际运算中，为了提高计算效率，采用并行的方式计算s_k的迭代收缩算子H_k：

$H_k=\frac{(V-w-U)+\sqrt{{\left(\right|V|-w-U)}^2+4w\left|V\right|}}{2\mathrm{sign}\left(V\right)}$

式中，V＝diag^-1(M^TM)M^T(Φ^*-Ms_k)+s_k，U＝λdiag^-1(M^TM)，diag表示提取矩阵的对角元素。

步骤10：结合以上计算出的迭代收缩算子，再利用线性搜索方法，就可以得到s_k+1：

s_k+1＝s_k+ε_kH_k，

${\mathrm{\ϵ}}_k=\underset{\mathrm{\ϵ}}{\arg \min }f(s_k+\mathrm{\ϵ}{H}_k),\mathrm{\ϵ}\∈(0,\∞).$

步骤11：计算比值

判断条件R＜0.1，若为假，则执行k＝k+1，并跳至步骤9，进入下一次迭代；若为真，则重建停止，得到目标体内部光强分布S的最优解，即最后一次迭代中得到的s_k+1。

在本发明中，将复杂的多维优化目标函数转化为高效的一维优化目标函数，利用平滑算子拟合l_p范数使求解有解析形式，并采用一种基于迭代收缩的算子以并行计算的形式全域重建出整个目标体内部光源，最终得到内部光强分布的定位和定量信息。与传统的方法相比，本发明对于给定不同的初值、较强的噪声、不同l_p范数、不同正则参数均可以重建得到准确的结果，对参数不敏感、鲁棒性好。此外，在可获得的先验信息较少时，特别是对于感兴趣区域较大时，较之传统方法可行域局限的缺点，本发明能够从整体上重建得到可靠的内部光强信息。

图5为使用本技术在不同正则参数下的重建结果。图5(a)是正则参数为10^-3下的三维整体结果显示；图5(c)-(f)分别代表正则参数为10^-3，10^-5，10^-7，10^-9下的重建结果横截面显示；图5(b)为相应的计算机断层成像解剖结构横截面视图。通过测试结果发现，所有的重建偏差都在1mm之内。证明本发明对于不同的正则化参数具有良好的鲁棒性。

图6为使用本技术在不同l_p范数下的重建结果。图6(a)-(d)的l_p范数分别对应为1.3，1.5，1.71.9。通过对l_p范数的不同选择，对比重建结果的横截面视图，重建光源分布基本相似，且在所有的情况下，光源位置和分布都与实际光源吻合的很好，偏差都在1mm之内。证明本发明对于不同l_p范数具有良好的鲁棒性。

图7为目标优化函数值与迭代次数的演化曲线，该图对应p＝1.3。在40步迭代中，目标函数的值单调下降，直到最后变化甚微。对于其他p下的演化曲线也与上述曲线类似，最终收敛性良好。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于迭代收缩的双模态融合断层成像方法，其特征在于，该方法包括以下步骤：

步骤a，采用波尔兹曼传输方程的一阶球谐近似得到的扩散方程建立光子在目标体内部的传输模型，即扩散方程传输模型，所述扩散方程传输模型表示为：

-▽·[D(r)▽Φ(r)]+μ_a(r)Φ(r)＝S(r)，

其中，r是位置向量，D是散射系数，μ_a是吸收系数，Φ是目标体表面光强分布，S表示待求解的目标体内部光强分布；

步骤b，利用有限元方法将扩散方程传输模型转换成线性矩阵方程的形式，所述线性矩阵方程表示为MS＝Φ，其中，M是描述模型性质的系统矩阵；

步骤c，基于所述步骤b转换得到的线性矩阵方程形式的扩散方程传输模型，利用正则理论，建立优化目标函数f(S)，其中，S为待求解的目标体内部光强分布，所述优化目标函数f(S)表示为：

$f\left(S\right)=\frac{1}{2}{\left|\right|\mathrm{MS}-{\mathrm{\Φ}}^{*}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|S\left|\right|}_p^p,$

其中，Φ^*为目标体表面光强分布的测量值，λ为正则化参数，

表示l_p范数；

步骤d，对f(S)进行降维处理，得到一维优化目标函数，所述一维优化目标函数g(α)表示为：

$g\left(\mathrm{\α}\right)=\frac{1}{2}{\left|\right|{\mathrm{Ms}}_k-{\mathrm{\Φ}}^{*}-m_j{s}_k\left(j\right)+m_j\mathrm{\α}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|\mathrm{\α}\left|\right|}_p^p,$

其中，s_k是对S的第k次估计，m_j是M的第j列元素，s_k(j)表示s_k的第j个元素，α表示s_k+1(j)的值；

其中，采用平滑算子θ(w)来拟合替代

以使所述一维优化目标函数g(α)在优化中存在解析解，其中，θ(w)＝|α|-wln(1+|α|/w)，且w是平滑算子参数；

步骤e，基于所述一维优化目标函数，以迭代的方式求解得到最优目标体内部光强分布S，所述步骤e包括以下步骤：

步骤e1，设置初始迭代次数k＝0，并为Φ^*，M，λ，w赋值，并初始化目标体内部光强分布s₀；

步骤e2，采用并行的方式计算s_k的迭代收缩算子H_k，所述H_k的计算公式为：

$H_k=\frac{(V-w-U)+\sqrt{{\left(\right|V|-w-U)}^2+4w\left|V\right|}}{2\mathrm{sign}\left(V\right)},$

其中，V＝diag^-1(M^TM)M^T(Φ^*-Ms_k)+s_k，U＝λdiag^-1(M^TM)，diag表示提取矩阵的对角元素；

步骤e3，利用线性搜索方法得到s_k+1；

步骤e4，计算比值

判断条件R＜0.1，若为假，则执行k＝k+1，并跳至步骤e2，进入下一次迭代，若为真，则重建停止，得到目标体内部光强分布S的最优解，即最后一次迭代中得到的s_k+1；

所述步骤b之前还包括建立非匀质目标体模型、对非匀质目标体模型进行几何结构的网格剖分，并利用配准技术将表面光强信息与几何结构信息进行位置融合与能量融合的步骤，所述建立非匀质目标体模型的步骤具体为：

利用自发荧光断层成像获取目标体表面光强分布的测量值Φ^*；利用计算机断层成像获取目标体内部几何结构信息，通过分割出目标体内部不同的器官结构，并为不同器官设置相应的光学参数，包括吸收系数和散射系数，从而建立非匀质目标体模型。

Images