CN102940482A - 自适应的荧光断层成像重建方法 - Google Patents

Abstract

一种自适应的荧光断层成像重建方法，利用有限元方法将扩散方程转化为线性方程；建立未知的荧光光源分布与表面荧光测量值之间的线性关系；计算当前正则化参数，并将残差相关系数中绝对值最大的元素选入支撑集合I；将更新后的支撑集合I中当前所有元素在矩阵A中对应的列取出，组成一个矩阵，获得下一步的搜索方向；计算下一步的步长并更新支撑集合；根据求解的搜索方向及步长迭代获取下一步结果，并更新正则化参数；判断是否达到停止条件，若达到，则重建过程结束，否则转到步骤S4。本发明不用***正则化参数，而是在重建过程中自适应地决定正则化参数。本发明通过自适应的正则化参数选取策略提高了重建的鲁棒性，大大提升了重建效率。

Description

自适应的荧光断层成像重建方法

技术领域

本发明涉及光学分子影像成像模态激发荧光断层成像(TomographicFluorescence Imaging，TFI)技术，特别涉及到一种自适应的荧光断层成像重建方法。

背景技术

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

TFI是一种典型的不适定问题，这是因为测量数据仅是表面的荧光分布，而需要求解的是整个成像空间的荧光光源分布，在这种情况下，光源的位置重建是没有唯一解的，并且对噪声非常敏感。为了能得到准确而稳定的重建结果，通常的方法是在优化问题中包含正则化项，该正则化项可以看作是光源分布的先验知识。最常见的正则化是Tikhonov正则化方法，通过在优化问题中加入L2范数约束，可以使TFI问题更加稳定。传统TFI重建方法根据经验值预先估计正则化参数的值，并使正则化参数在整个迭代重建过程保持不变。但不同的重建问题所要求的正则化参数通常不同，传统重建方法很难准确预测合适的正则化参数从而完成精确重建。因此合适的正则化参数选择仍是一项挑战性问题。

发明内容

本发明针对TFI重建技术中的正则化选取问题，提出了一种自适应正则化参数的快速鲁棒TFI重建方法。

为实现上述目的，一种自适应的荧光断层成像重建方法，包括：

S1利用有限元方法将扩散方程转化为线性方程；

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

S3计算当前正则化参数，并将残差相关系数中绝对值最大的元素选入支撑集合I；

S4将更新后的支撑集合I中当前所有元素在矩阵A中对应的列取出，组成一个矩阵，获得下一步的搜索方向；

S5计算下一步的步长并更新支撑集合；

S6根据求解的搜索方向及步长迭代获取下一步结果，并更新正则化参数；

S7判断是否达到停止条件，若达到，则重建过程结束，否则转到步骤S4。

本发明不用***正则化参数，而是在重建过程中自适应地决定正则化参数。本发明不仅通过自适应的正则化参数选取策略提高了重建技术的鲁棒性，而且大大提升了重建效率。

附图说明

图1是本发明方法的流程图；

图2是多光源仿体的三维结构图和位于z＝0平面截面图；

图3是4组测量数据情况下的单光源重建结果；

图4是4组测量数据情况下的双光源重建结果；

图5是4组测量数据情况下的三光源重建结果。

具体实施方式

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

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

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

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

步骤103：计算当前正则化参数，并将残差相关系数中绝对值最大的元素选入支撑集合；

步骤104：通过求解方程获得下一步的搜索方向；

步骤105：计算下一步的步长并更新支撑集合；

步骤106：根据求解的搜索方向及步长迭代获取下一步结果，并更新正则化参数；

步骤107：判断是否达到停止条件。若达到，则重建过程结束；否则转到步骤104。

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

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

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

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

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

${\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是吸收系数，D_x，m＝1/3(μ_ax，am+μ′_sx，sm)是扩散系数，μ′_sx，sm是约化散射系数，Φ_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转化为具有如下形式的优化问题：

$\underset{X}{\min }{F}_{\mathrm{\λ}}\left(X\right)=\frac{1}{2}{\left|\right|\mathrm{AX}-\mathrm{\Φ}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|X\left|\right|}_1$

λ为正则化参数，F_λ(X)为待优化的目标函数。

步骤103：为实现对目标函数的最小化，首先对F_λ(X)求导得到：

${\∂}_{X_{\mathrm{\λ}}}{F}_{\mathrm{\λ}}\left(X_{\mathrm{\λ}}\right)=A^T({\mathrm{AX}}_{\mathrm{\λ}}-\mathrm{\Φ})+{\mathrm{\λ}\∂}_{X_{\mathrm{\λ}}}{\left|\right|X_{\mathrm{\λ}}\left|\right|}_1$

其中

是||X_λ||₁的偏微分，且由下式确定：

${\∂}_{X_{\mathrm{\λ}}}{\left|\right|X_{\mathrm{\λ}}\left|\right|}_1=\{u\∈R^n|\left.\begin{array}{cc}{u}_i=\mathrm{sgn}\left(X_{\mathrm{\λ},i}\right),& X_{\mathrm{\λ},i}\≠0\\ u_i\∈[-\mathrm{1,1}],& X_{\mathrm{\λ},i}\≠0\end{array}\right\}$

集合I＝{i：X_λ(i)≠0}用来表示未知向量X_λ的支撑集合。当前残差的相关性系数c由下式确定：

c＝A^T(Φ-AX_λ)

获取残差的相关系数c后，再确定初始的正则化参数λ：

I＝{j：|c_k(j)|＝||c_k||_∞＝λ}

其中j为向量c中绝对值最大的元素，将第j元素加入支撑集合I中，并使正则化参数λ赋值为c中绝对值最大的元素的值。

步骤104：将更新后的支撑集合I中当前所有元素在矩阵A中对应的列取出，组成一个矩阵A_I。接下来，通过解以下等式，获得下一步的搜索方向p_k：

$A_I^T{A}_I{p}_k\left(I\right)=\mathrm{sgn}\left(c_k\left(I\right)\right)$

在获得p_k的值后，将p_k向量中不在支撑集合I中的分量元素的值置为0，从而得到p_k(I)。

步骤105：计算下一步搜索方向的步长γ_k，并更新支撑集合。

步骤106：在步骤104获取了下一步搜索方向p_k，以及步骤105获取了下一步搜索方向上的步长γ_k后，可以更新目前的未知数X_k：

X_k＝X_k-1+γ_kp_k

获取新的迭代结果以后，接下来将自动更新正则化参数λ：

λ＝λ-γ_k

步骤107：该步骤判断迭代停止的条件是否达到。计算更新后目标函数的值F_λ(X_k)，而更新前目标函数值为F_λ(X_k-1)。若|_Fλ(X_k)-F_λ(X_k-1)|/F_λ(X_k)＜T，则停止迭代，否则转到步骤104开始下一轮迭代。其中阈值T的取值为0.0001～0.01。

其中，步骤105包括以下5个子步骤：

步骤201：计算步长

${\mathrm{\γ}}_k^{+}=\underset{i\∈I^c}{\min }\left\{\frac{\mathrm{\λ}-c_k\left(i\right)}{1-a_i^T{A}_I{p}_k\left(I\right)}\right\}$

其中I^c是I的补集。a_i是矩阵A的第i列。获得的同时纪录最小值出现的位置pos+。

步骤202：计算步长

${\mathrm{\γ}}_k^{-}=\underset{i\∈I}{\min }\{-X_k\left(i\right)/p_k\left(i\right)\}$

获得的同时纪录最小值出现的位置pos-。

步骤203：比较

与

的大小。如果

转到步骤204；否则转到步骤S205；

步骤204：根据纪录的位置pos+找到相应的元素，将其从补集I^c中移到支撑集合I中，从而支撑集合元素的个数增加，并使

步骤205：根据纪录的位置pos-找到相应的元素，将其从支撑集合I中移除到补集I^c中，从而支撑集合元素的个数减少，并使

运行结果

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

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

我们使用了4个点状激发源生成荧光测量数据，图4、图5及图6分别显示了单光源、双光源及三光源的重建结果。重建时间分别为0.39秒，0.42秒及0.48秒。结果表明，本重建方法可以不用***正则化参数，而是在重建过程中自适应的决定正则化参数的值，能够得到十分准确的结果，且具有较快的重建速度。

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

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

Claims (3)

1.一种自适应的荧光断层成像重建方法，包括： 

S1利用有限元方法将扩散方程转化为线性方程； 

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

S3计算当前正则化参数，并将残差相关系数中绝对值最大的元素选入支撑集合I； 

S4将更新后的支撑集合I中当前所有元素在矩阵A中对应的列取出，组成一个矩阵，获得下一步的搜索方向； 

S5计算下一步的步长并更新支撑集合； 

S6根据求解的搜索方向及步长迭代获取下一步结果，并更新正则化参数； 

S7判断是否达到停止条件，若达到，则重建过程结束，否则转到步骤S4。 

2.如权利要求1所述的方法，其特征在于所述步骤S5包括： 

计算步长

并记录最小值出现的位置pos+； 

计算步长

并记录最小值出现的位置pos-； 

比较

与的大小。如果

转到步骤S4；否则转到步骤S5； 

将补集I^c中的pos+位置的元素加入到支撑集合I中，并使

将支撑集合I中的pos-位置的元素移除到其补集I^c，并使

3.如权利要求1所述的方法，其特征在于判断是否达到停止条件的阈值为0.0001-0.01。 

Images