CN102034266B - 激发荧光断层成像的快速稀疏重建方法和设备 - Google Patents

Abstract

一种激发荧光断层成像的快速稀疏重建方法，包括步骤：利用有限元理论首先将扩散方程表示为线性化方程；建立未知的荧光光源分布与边界测量数据之间的线性关系；计算余量相关度向量得到最相关元素集合；将最相关元素集合与当前支撑集合并，生成新的支撑集；利用支撑集将离散后的成像空间分为允许区域和禁止区域，建立表面荧光数据与允许区域的线性关系；将最终得到的解向量中的负元素替换为0。本发明基于扩散近似模型，充分地考虑了生物组织的非匀质特性。在光源重建过程中，基于L1范数的稀疏性约束，并将TFI问题看作是压缩感知问题，利用基于支撑集的重建方法进行光源定位，有效地避免了重建结果的过平滑现象，提高了TFI成像的精度。

Description

激发荧光断层成像的快速稀疏重建方法和设备

技术领域

本发明涉及光学分子影像成像模态-激发荧光断层成像技术，特别涉及一种实现稀疏性正则化的快速荧光光源重建方法。

背景技术

作为一种光学分子影像成像模态，激发荧光成像技术已经得到了迅速的发展和广泛的应用。通过对感兴趣区域的生物组织标记光学分子探针，并对体表的荧光信号进行采集，我们可以以一种无创的方法来获取生物体分子细胞水平的信息。然而，由于在可见光和近红外光范围内，光子在生物组织内传播时会产生严重的散射现象，因此传统的二维激发荧光成像技术无法显示荧光光源的准确位置。而激发荧光断层成像(TFI)是一种三维重建技术，通过采集表面荧光数据并基于特定的逆问题模型，可以实现荧光光源的三维精确定位。

TFI是一种典型的不适定问题，这是因为测量数据仅是表面的荧光分布，而需要求解的是整个成像空间的荧光光源分布，在这种情况下，光源的位置重建是没有唯一解的，并且对噪声非常敏感。为了能得到准确而稳定的重建结果，通常的方法是在优化问题中包含正则化项，该正则化项可以看作是光源分布的先验知识。最常见的正则化是Tikhonov正则化方法，通过在优化问题中加入L2范数约束，可以使TFI问题更加稳定，然而Tikhonov正则化的缺点是结果过于平滑，造成成像精度下降。为了解决这一问题，考虑到TFI经常被用来进行肿瘤的早期检测，并且在早期阶段肿瘤通常是很小的和稀疏的，因此可以将光源的稀疏性约束融入到TFI重建过程当中。显然，稀疏性约束可以通过使用L0范数约束实现，但此时TFI问题变为了一个NP-Hard问题，难以进行高效的求解。然而当荧光光源足够稀疏的时候，可以使用L1范数代替L0范数，并且能够取得同样的效果。

发明内容

本发明的目的是提供一种激发荧光断层成像的快速稀疏重建方法和设备。

为实现上述目的，一种激发荧光断层成像的快速稀疏重建方法，包括步骤：

S1：利用有限元理论首先将扩散方程表示为线性化方程；

S2：建立未知的荧光光源分布与边界测量数据之间的线性关系；

S3：计算余量相关度向量并利用阈值的方法得到最相关元素集合；

S4：将最相关元素集合与当前支撑集合并，生成新的支撑集；

S5：利用支撑集将离散后的成像空间分为允许区域和禁止区域，建立表面荧光数据与允许区域的线性关系；

S6：将最终得到的解向量中的负元素替换为0。

本发明基于扩散近似模型，充分地考虑了生物组织的非匀质特性。在光源重建过程中，基于L1范数的稀疏性约束，并将TFI问题看作是压缩感知问题，利用基于支撑集的重建方法进行光源定位，有效地避免了重建结果的过平滑现象，提高了TFI成像的精度。另一方面，在数值计算过程中，本发明充分考虑到光源分布的稀疏性特点，避免了大量不必要的数值计算，使得光源重建的速度大大提高，此方法对激发荧光断层成像的发展具有重要的应用价值。

附图说明

图1，本发明方法的流程图。

图2，计算最相关元素集合的方法流程。

图3，具有4种不同区域的圆柱形非匀质仿体。

图4，单光源仿体和双光源仿体的三维结构图和位于z＝0平面的截面图。

图5，15组测量数据情况下的单光源重建结果。

图6，15组测量数据情况下的双光源重建结果。

图7，3组测量数据情况下的单光源重建结果。

图8，3组测量数据情况下的双光源重建结果。

具体实施方式

下面将结合附图详细描述本发明的重建方法，应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

图1是本发明的重建方法的总体流程：

步骤101：本发明采用扩散方程作为TFI的成像模型，模型包括激发过程和发射过程两个相耦合的扩散方程，通过使用有限元方法，可以将扩散方程离散化，以便于进行后续的处理；

步骤102：对于TFI逆问题而言，表面荧光分布是已知的，而内部荧光光源分布是未知的，通过对离散化后的扩散模型进行矩阵变换，可以建立表面测量数据与未知光源分布之间的线性方程。同时对于多组测量数据的情况，可以将多个线性方程组合为统一的方程；

步骤103：通过计算余量相关度向量并利用阈值方法得到最相关元素集合，如果该集合为空集，说明无法进一步减小当前余量，于是转到步骤106；

步骤104：基于步骤103中得到的最相关元素集合，更新当前支撑集，若新的支撑集不满足稀疏性约束(即支撑集的元素个数超过某一特定阈值)，则转向步骤106；

步骤105：将TFI的成像空间限制在当前支撑集所包含的范围，从而大大降低了TFI优化问题的规模。在新的成像空间中，建立表面测量数据与未知光源分布之间的线性方程，并计算最小二乘解。若循环次数尚未达到最大循环次数预设值，则转到步骤103；

步骤106：根据荧光光源的非负性特点，将重建得到的荧光光源分布向量中的负元素赋值为0。

如图2所示，步骤103包括以下子步骤：

步骤201：利用当前的解向量计算表面荧光预测值，并计算预测值与真实荧光测量值的误差，将该误差作为当前的重建余量；

步骤202：通过计算余量相关度向量，分析荧光光源分布向量中的不同元素与当前余量之间的相关程度；

步骤203：利用阈值方法将余量相关度向量中的元素分为两部分，将大于阈值部分的元素的编号组成最相关元素集合。

下面对本发明的TFI重建方法所涉及的关键步骤进行逐一详细说明，具体形式如下面所述：

步骤101：利用有限元方法将扩散方程转化为线性方程；

精确描述光在非匀质生物组织中传输的数学模型是辐射传输方程，它是一个复杂的微分-积分方程，求解十分困难。幸运的是，在可见光和近红外光谱段，光子在生物组织中传输时具有高散射、低吸收的特点，在这种情况下，扩散方程可以很好地近似辐射传输方程，而扩散方程的计算复杂度较低，可以在有限时间内进行求解，所以非常适合作为TFI成像模型。当利用连续波聚焦激光器作为激发源时，可以使用如下的两个相耦合的扩散方程对TFI的成像过程进行描述：

$\left\{\begin{array}{c}-\▿\·[D_x\left(r\right)\▿{\mathrm{\Φ}}_x\left(r\right)]+{\mathrm{\μ}}_{\mathrm{ax}}\left(r\right){\mathrm{\Φ}}_x\left(r\right)=\mathrm{\Θ\δ}(r-r_l)\\ -\▿\·[D_m\left(r\right)\▿{\mathrm{\Φ}}_m\left(r\right)]+{\mathrm{\μ}}_{\mathrm{am}}\left(r\right){\mathrm{\Φ}}_m\left(r\right)={\mathrm{\Φ}}_x\left(r\right)\mathrm{\η}{\mathrm{\μ}}_{\mathrm{af}}\left(r\right)\end{array}\right.(r\∈\mathrm{\Ω})$

上述扩散方程的边界条件可以描述为：

${\mathrm{\Φ}}_{x,m}\left(r\right)+2q{D}_{x,m}\left(r\right)[\stackrel{\→}{v}\left(r\right)\·{\mathrm{\Φ}}_{x,m}\left(r\right)]=0(r\∈\∂\mathrm{\Ω})$

其中，Ω表示整个成像空间，

表示空间边界，下标x和m分别表示激发光和发射光，μ_ax，am是吸收系数，

是扩散系数，

是约化散射系数，Φ_x，m表示光子密度，Θδ(r-r_l)表示各向同性的点状激发光源，Θ表示光源的强度，

是边界的外向单位法向量，q是取决于边界光学反射系数偏差的特定常量，ημ_af表示待重建的荧光量子产额。

在有限元理论框架下，上述扩散方程及其边界条件首先被表示为如下的弱解形式：

$\left\{\begin{array}{c}{\∫}_{\mathrm{\Ω}}(D_x\left(r\right)\▿{\mathrm{\Φ}}_x\left(r\right)\·\▿\mathrm{\Ψ}\left(r\right)+{\mathrm{\μ}}_{\mathrm{ax}}\left(r\right){\mathrm{\Φ}}_x\left(r\right)\mathrm{\Ψ}\left(r\right))\mathrm{dr}\\ +{\∫}_{\∂\mathrm{\Ω}}\frac{1}{2A(r;n,n^{\′})}{\mathrm{\Φ}}_x\left(r\right)\mathrm{\Ψ}\left(r\right)\mathrm{dr}={\∫}_{\mathrm{\Ω}}\mathrm{\Θ\δ}(r-r_l)\mathrm{\Ψ}\left(r\right)\mathrm{dr}\\ {\∫}_{\mathrm{\Ω}}(D_m\left(r\right)\▿{\mathrm{\Φ}}_m\left(r\right)\·\▿\mathrm{\Ψ}\left(r\right)+{\mathrm{\μ}}_{\mathrm{am}}\left(r\right){\mathrm{\Φ}}_m\left(r\right)\mathrm{\Ψ}\left(r\right))\mathrm{dr}\\ +{\∫}_{\∂\mathrm{\Ω}}\frac{1}{2A(r;n,n^{\′})}{\mathrm{\Φ}}_m\left(r\right)\mathrm{\Ψ}\left(r\right)\mathrm{dr}={\∫}_{\mathrm{\Ω}}{\mathrm{\Φ}}_x\left(r\right){\mathrm{\η\μ}}_{\mathrm{af}}\left(r\right)\mathrm{\Ψ}\left(r\right)\mathrm{dr}\end{array}\right.$

其中，Ψ(r)表示任意测试函数。然后，利用四面体对整个成像空间进行剖分，并将相应的基函数作为测试函数Ψ(r)，从而上述弱解形式可以被离散化为如下的两个等式：

K_xΦ_x＝S_x

K_mΦ_m＝FX

其中，K_x，m是系统矩阵，矩阵F是通过对未知的荧光光源量子产额离散化得到的，向量X表示需要重建的光源量子产额。

步骤102：建立未知的荧光光源分布与表面荧光测量值之间的线性关系；

对于激发过程而言，光子密度Φ_x可以通过求解K_xΦ_x＝S_x而直接获得，得到的Φ_x将作为发射过程的能量源。对于TFI逆问题而言，由于K_m为对称正定矩阵，因此可以得到如下的矩阵方程：

${\mathrm{\Φ}}_{m,l}=K_{m,l}^{-1}\mathrm{FX}=B_{l}X$

去除向量Φ_m，l中的非表面测量元素以及矩阵B_l中对应的行，上述等式可以进一步转化为：

${\mathrm{\Φ}}_{m,l}^{\mathrm{meas}}=A_{l}X$

对于TFI而言，通常使用多个点状激发光源来得到多组荧光测量数据，每一组数据都可以得到上述的矩阵方程，将多个方程组装起来可以得到如下的方程：

Φ＝AX

其中，

$\mathrm{\Φ}=\left\{\begin{array}{c}{\mathrm{\Φ}}_{m,1}^{\mathrm{meas}}\\ {\mathrm{\Φ}}_{m,2}^{\mathrm{meas}}\\ .\\ .\\ .\\ {\mathrm{\Φ}}_{m,L}^{\mathrm{meas}}\end{array}\right\},$ $A=\left[\begin{array}{c}{A}_1\\ A_2\\ .\\ .\\ .\\ A_L\end{array}\right]$

由于TFI是一个不适定的逆问题，因此需要融入一定的正则化方法来使TFI问题更加稳定，本发明方法采用的是L1范数稀疏性约束，在该约束条件下，TFI转化为具有如下形式的优化问题：

min‖X‖₁subject to AX＝Φ

这里，如果将矩阵A的所有列看作是一个基函数集合，并将向量X中的元素看作是对应于每个基函数的系数，则TFI逆问题可以看作是压缩感知问题，在这种情况下，最优解X_opt仅包含有限个较大的系数。

步骤103：通过计算余量相关度向量并采用阈值方法得到最相关元素集合S_n，若S_n为空集，则转向步骤106；

为了进行未知光源重建，这里采用的方法是首先计算最相关元素集合S_n，集合S_n中包含的是某些元素的编号，相比那些不在S_n中的元素，通过改变位于S_n中的元素的数值，能够在更大程度上减小当前的重建误差。当集合S_n为空集时，说明改变当前重建结果的值无法在很大程度上减小当前重建误差，即当前重建结果的精度已经达到预期要求，因此算法会直接转向步骤106。

步骤104：将集合S_n与当前支撑集I_n-1合并，生成新的支撑集I_n，根据I_n中元素的个数决定算法是否停止；

在步骤103中计算得到的集合S_n指明了通过改变哪些元素的值可以更好地减小当前重建余量r_n-1，因此需要将集合S_n与当前支撑集I_n-1合并，从而构成新的支撑集I_n：

I_n＝I_n-1∪S_n

在本重建方法中，荧光光源的稀疏性特征被作为先验知识融入到重建过程中。为了实现稀疏性成像，需要指定支撑集I_n中元素的最大个数P_max，对于新的支撑集I_n，若|I_n|＞P_max，则表明下一次迭代得到的结果将不满足稀疏性约束，因此算法会直接转向步骤106。

步骤105：将TFI的重建空间限制在I_n所指出的范围，并建立表面测量数据Φ与荧光光源分布

之间的线性方程，利用共轭梯度法计算其最小二乘解。若当前循环次数n尚未达到最大循环次数N_max，则转到步骤103；

在步骤104中已经计算得到了新的支撑集I_n，并且假设荧光光源仅存在于I_n所指出的空间范围，从而TFI优化问题的规模大大降低了，因此重建的不适定性也得到了极大的缓解。这里，首先将矩阵A的列向量中编号不在I_n中的向量去掉，同时也将向量X中对应的元素去掉，从而可以得到表面荧光测量数据与X中的可允许元素之间的线性关系：

$A_{I_n}{X}_{I_n}=\mathrm{\Φ}$

因为|I_n|远小于向量X的维数，因此上述方程通常是超定的(Over-Determined)，由于噪声等原因产生的误差，使得上述方程往往矛盾的，无法得到精确的解。因此，本方法将计算其的最小二乘解，即求解如下形式的方程：

$A_{I_n}^T{A}_{I_n}{X}_{I_n}=A_{I_n}^T\mathrm{\Φ}$

在本方法中，将使用共轭梯度法对上述方程进行迭代求解，迭代的次数被设定为|I_n|，然后X_n可以根据

构建出来。如果算法当前的迭代次数尚未达到预设值，则转到步骤103。

步骤106：将重建结果中的负元素替换为零元素。

本发明方法在重建过程中并未附加非负性约束，在实际应用中，由于荧光探针量子产额具有非负性，因此在本重建方法结束时，会将重建结果中的负元素置为0。

其中，步骤103包括以下3个子步骤：

步骤201：利用当前的迭代结果计算重建余量；

在重建方法的第n步迭代中，需要对当前重建结果的准确性进行评价，对于TFI成像而言，重建得到的结果是荧光光源的分布，而在实际应用中无法预先知道真实的光源位置，因此可以基于当前重建得到的光源分布来计算相应的表面荧光数据，将该计算得到的表面数据作为预测值，与真实的测量数据进行比较，从而评价当前的重建结果。具体来说，需要计算如下的重建余量：

r_n-1＝Φ-AX_n-1

假设A是一个m×k的稠密矩阵，于是计算r_n-1需要m×k个乘法运算，如果考虑了X_n-1的稀疏性，则相应的计算量可以大幅度削减。假设X_n-1中有k′个非零元素，则计算r_n-1实际上仅需要m×k′个乘法运算，因此基于稀疏性的数值运算可以大大提高算法的计算效率。

步骤202：计算余量相关度向量，分析荧光光源分布向量中的不同元素与当前余量之间的相关程度；

在步骤201中计算得到的重建余量可以反映目前重建的荧光光源分布的精确程度，重建余量r_n-1在尺度上越小，表明当前结果越精确。为了能够进一步减小重建余量r_n-1，可以首先考察荧光光源分布X中的每个元素与r_n-1的相关程度，相关程度高表明通过改变该元素的值能够很好地减小r_n-1的尺度。余量相关度向量可以通过如下公式计算：

c_n＝A^Tr_n-1

相关度向量c_n中的正值表示正相关，负值表示负相关。然而，在实际应用中考察的是相关度的绝对值。

步骤203：通过使用阈值方法，将相关度向量中的元素分为两部分，在重建过程中只考察大于阈值的部分。

在步骤202中计算得到的相关度向量c_n表征了荧光光源分布X中的每个元素与余量的相关程度，通过合理地改变相关程度高的元素的数值，可以更好地减小当前余量，从而使重建结果更加精确。

在本步骤中，首先设定相关度阈值，对于当前迭代而言，相关度阈值设定为t_n＝α‖c_n‖_∞，其中α∈(0，1)为预先设定的阈值参数。这种设定阈值的好处是根据每次迭代重建的结果不同，t_n会自适应地变化，从而能够获得更加理想的效果。在实际应用中，最优的α通常在0.7到0.9之间变化。

计算得到当前相关度阈值t_n后，可以将相关度向量c_n中的元素分为两个部分，将大于阈值t_n的元素的序号组成集合S_n。这里需要指出的是，在生成集合S_n时，需要去掉已经包含在当前支撑集中的元素，于是可以得到如下的集合S_n的表达式：

运行结果

为了验证本发明方法的精确性，我们利用三维非匀质仿真仿体进行了TFI重建实验。

在本实验中，采用的非匀质仿体为圆柱形，其中包含了4个区域，分别代表肌肉(M)、肺(L)、心脏(H)和骨骼(B)，表1给出了不同组织对应的光学系数。该仿体的直径和高度均为20mm，图3显示了该仿体的三维结构图和位于z＝0平面的截面图。在本实验中，我们分别进行了单光源和双光源的重建实验，如图4所示，所有光源均为球形，直径为2mm，中心位于z＝0平面内，光源的量子产额设定为0.5，图3的截面图中的黑色圆点表示点状激发光源的位置，对于每个点光源，表面荧光分布是在对面的160°视野内的圆柱表面测量的。为了对光源进行重建，该非匀质仿体被离散化为2710个点和13316个四面体。

在第一组实验中，我们使用了15个点状激发源生成荧光测量数据，图5和图6显示了单光源和双光源的重建结果，重建时间分别为0.062秒和0.073秒，从结果中可以看出，本方法的精度是比较高的。在第二组实验中，我们仅采用了由编号1、6、11指明的点光源来生成荧光测量数据，从而来模拟一种严重的不适定性情况，图7和图8显示了对应的重建结果，重建时间分别为0.023秒和0.032秒，结果表明，本重建方法依然能够得到十分准确的结果，这表明本发明方法对严重不适定问题依然是适用的。

表1仿体中不同区域的光学系数(单位：mm^-1)

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (10)

1.一种激发荧光断层成像的快速稀疏重建方法，包括步骤：

S1：利用有限元理论首先将扩散方程表示为线性化方程；

S2：建立未知的荧光光源分布与边界测量数据之间的线性关系；

S3：计算余量相关度向量并利用阈值的方法得到最相关元素集合；

S4：将最相关元素集合与当前支撑集合并，生成新的支撑集；

S5：利用支撑集将离散后的成像空间分为允许区域和禁止区域，建立表面荧光数据与允许区域的线性关系；

S6：将最终得到的解向量中的负元素替换为0。

2.按权利要求1所述的方法，其特征在于所述扩散方程表示为如下的弱解形式：

$\left\{\begin{array}{c}{\∫}_{\mathrm{\Ω}}(D_x\left(r\right)\▿{\mathrm{\Φ}}_x\left(r\right)\·\▿\mathrm{\Ψ}\left(r\right)+{\mathrm{\μ}}_{\mathrm{ax}}\left(r\right){\mathrm{\Φ}}_x\left(r\right)\mathrm{\Ψ}\left(r\right))\mathrm{dr}\\ +{\∫}_{\∂\mathrm{\Ω}}\frac{1}{2A(r;n,n^{\′})}{\mathrm{\Φ}}_x\left(r\right)\mathrm{\Ψ}\left(r\right)\mathrm{dr}={\∫}_{\mathrm{\Ω}}\mathrm{\Θ\δ}(r-r_l)\mathrm{\Ψ}\left(r\right)\mathrm{dr}\\ {\∫}_{\mathrm{\Ω}}(D_m\left(r\right)\▿{\mathrm{\Φ}}_m\left(r\right)\·\▿\mathrm{\Ψ}\left(r\right)+{\mathrm{\μ}}_{\mathrm{am}}\left(r\right){\mathrm{\Φ}}_m\left(r\right)\mathrm{\Ψ}\left(r\right))\mathrm{dr}\\ +{\∫}_{\∂\mathrm{\Ω}}\frac{1}{2A(r;n,n^{\′})}{\mathrm{\Φ}}_m\left(r\right)\mathrm{\Ψ}\left(r\right)\mathrm{dr}={\∫}_{\mathrm{\Ω}}{\mathrm{\Φ}}_x\left(r\right)\mathrm{\η}{\mathrm{\μ}}_{\mathrm{af}}\left(r\right)\mathrm{\Ψ}\left(r\right)\mathrm{dr}\end{array}\right.$

其中，Ψ(r)表示任意测试函数，Ω表示整个成像空间，表示空间边界，下标x和m分别表示激发光和发射光，μ_ax，am是吸收系数，

是扩散系数，

是约化散射系数，Φ_x，m表示光子密度，Θδ(r-r_l)表示各向同性的点状激发光源，Θ表示光源的强度，ημ_af表示待重建的荧光量子产额。

3.按权利要求1所述的方法，其特征在于所述步骤S2包括：

将由多组测量数据生成的线性方程组装成一个统一的线性方程。

4.按权利要求1所述的方法，其特征在于，所述步骤S3包括步骤：

S3-1：计算当前余量r_n-1＝Φ-AX_n-1；

S3-2：计算余量相关度向量c_n＝A^Tr_n-1；

S3-3：将余量相关度向量中大于特定阈值的元素的编号组成最相关元素集合S_n：

其中，Φ表示光子密度，A表示系统矩阵，X_n-1表示当前荧光光源分布，t_n表示当前相关阈值。

5.按权利要求1所述的方法，其特征在于，基于当前支撑集I_n，建立表面荧光测量数据Φ与支撑集范围内的荧光光源分布

之间的线性关系：

$A_{I_n}{X}_{I_n}=\mathrm{\Φ}.$

6.按权利要求5所述的方法，其特征在于利用共轭梯度法计算所述线性关系的最小二乘解：

$A_{I_n}^T{A}_{I_n}{X}_{I_n}=A_{I_n}^T\mathrm{\Φ}.$

7.按权利要求1所述的方法，其特征在于步骤S3还包括：如果计算得到的集合S_n为空集

时，则转到步骤S6。

8.按权利要求1所述的方法，其特征在于所述步骤S4还包括：

当目前的支撑集中的元素个数大于可允许的最大数目P_max时，则转到步骤S6。

9.按权利要求1所述的方法，其特征在于所述步骤S5还包括：如果迭代次数未达到预设的最大数目N_max，则转到步骤S3。

10.一种激发荧光断层成像的快速稀疏重建设备，包括：

利用有限元理论首先将扩散方程表示为线性化方程的装置；

建立未知的荧光光源分布与边界测量数据之间的线性关系的装置；

计算余量相关度向量并利用阈值的方法得到最相关元素集合的装置；

将最相关元素集合与当前支撑集合并，生成新的支撑集的装置；

利用支撑集将离散后的成像空间分为允许区域和禁止区域，建立表面荧光数据与允许区域的线性关系的装置；

将最终得到的解向量中的负元素替换为0的装置。

Images