基于多任务贝叶斯压缩感知方法的生物发光断层成像重建算法

摘要

一种基于多任务贝叶斯压缩感知方法的生物发光断层成像重建算法，属于医学图像处理领域，该方法首先采用高阶近似模型模型光在生物组织中的传输规律，并基于多任务学习方法探索多光谱之间的内在相关性，将多光谱之间的相关性作为先验信息融入重建算法中以降低BLT重建的病态性，最后在此基础上实现了荧光光源的三维重建。本发明方法与其他生物发光断层成像重建算法相比较，进一步地融入了多光谱相关性，降低了BLT重建病态性，不仅能对荧光光源进行准确重建定位，而且可以极大地提高了计算效率。

说明书

技术领域

本发明属于医学图像处理领域，涉及一种基于多任务贝叶斯压缩感知方法的生物发光 断层成像重建算法。

背景技术

光学分子影像是一种快速发展的分子影像技术，它通过将光学过程与一定的分子性质 相结合，对目标体内的生物荧光或激发荧光进行分析和处理，并进行定性和定量研究。

光学分子影像技术中最具代表性的成像方式是激发荧光成像技术(Fluorescence Imagi  ng)和生物发光成像技术(Bioluminescence Imaging,BLI)。它们均是二维的生物发光成像 技术，虽然这种成像技术应用起来方便、简单，但这种二维成像方式在应用过程中存在局 限性，尤其是对于成像深度的限制，二维的荧光图像不能反映光源深度信息和且难以定量 化。它们仅可以反映生物体内的荧光探针在某一角度的投影信息，而这种投影信息是多个 深度的信号叠加而成的。因此，二维成像方法具有很低的分辨率。

在二维生物发光成像技术基础上发展起来的生物发光断层成像技术(Bioluminescence  Tomography,BLT)由于能够反映信号深度的信息，已经成为光学分子影像的一个重要分支。 生物发光断层成像技术不需要外在光源的激发，而是通过一种生物化学发光反应在体内发 光。体内产生的荧光在生物组织内部以某种规律传播并不断地与生物组织发生相互作用， 并到达体表。最后，在生物组织体表利用高灵敏度的探测器获得的荧光图像就可以重建出 荧光光源在小动物体内的分布情况，从而在本质上揭示在体分子的活动规律。

但光子在生物组织中不沿直线传输，而是经历了大量的散射过程，导致BLT逆问题在 数学上是一个高度病态的问题，外界微小的测量扰动都会给重建结果带来很大的变化。国 内外研究人员为降低其病态性做了很多工作，一般都是从不同角度为该问题提供各种不同 的先验知识。

现有的降低其病态性的方法大多是基于多光谱信息和光源的稀疏特性展开BLT重建方 法的研究，但是，这些方法并没有考虑多光谱之间的相关性。在数学上，多光谱方法增加 了待求解方程的已知量，因此可以改善求解的结果。对于BLT而言，基于多光谱的方法是 在光源保持不变，利用多个滤波片来获取不同波段的数据，从而进行光源的重建。实际上， 多个光谱之间的测量数据不是孤立的，而是相关关联的，它们的共同任务是重建荧光光源。 因此，如果在重建时考虑到多光谱之间的相关性，将有可能改善BLT重建图像质量。基于 此，本发明提出了一种基于多光谱内在相关性的BLT重建方法。

本发明基于多任务学习方法，首先采用高阶近似模型模型光在生物组织中的传输规律， 并基于多任务学习方法探索多光谱之间的内在相关性，最后在此基础上实现了荧光光源的三 维重建。

发明内容

针对现有荧光断层成像重建算法中存在的上述问题，本发明提出了一种基于多任务学习 的方法，它采用高阶近似模型模型光在生物组织中的传输规律，并基于多任务学习方法探索 多光谱之间的内在相关性，最后在此基础上实现了荧光光源的三维重建。

本发明采用的技术方案为：首先采用高阶近似模型模型光在生物组织中的传输规律，并 基于多任务学习方法探索多光谱之间的内在相关性，将多光谱之间的相关性作为先验信息融 入重建算法中以降低BLT重建的病态性，最后在此基础上实现了荧光光源的三维重建。具体 包括下述步骤：

步骤一，问题定义及初始化——获得P个谱段在M个测量点上的值，并设置α₀服从的伽 马先验分布的形状参数a和尺度参数b，其中为出射光测量值Φ(τ_i)服从的高斯分布的方差；

在模拟光在生物组织中传输规律时可以采用两种模型：辐射传输方程和扩散方程。辐 射传输方程是一种精确描述光传播规律的数学模型，但是即使在非常简化的条件下RTE 的求解仍然十分困难，且计算量非常巨大。扩散方程是辐射传输方程的低阶近似，因此其 求解简单，但求解精度低。进一步，扩散方程在应用时必须满足生物组织强散射、低吸收 的条件，而且不适用于光源位于生物组织边界时的情况，限制了扩散方程在生物医学领域 中的应用。基于此，本发明选择高阶球谐近似(Simplified Spherical Harmonics Approxima  tion，SPN)方程。该模型能够较为精确的描述光在强吸收低散射的生物组织中传输，同时， 它还保持着适中的计算复杂度。用该模型替代扩散近似传输方程以准确地模拟光在生物组 织中的传输规律，以提高前向求解的精度。

考虑到光谱影响，在波段为τ_i时，位置为r时，SP7模型共有四个方程：

其中，S代表光源函数；μ_a是光吸收系数；μ_s是光散射系数；μ_ai＝μ_a+μ_s(1-gⁱ)，gⁱ为各 向异性因子；代表辐射度的勒让德矩(Legendre Moments)φ_i的线性组合，可表示为：

在SP7方程的基础上，令φ₆＝φ₄＝0，并保留公式(1)和(2)，即可得到SP3方程：

在SP7的边界条件的基础上，令φ₆＝φ₄＝0，得到SP3的边界条件：

其中，代表垂直边界向外的单位法向量；A_i、B_i和C_i为一系列常数，与边界反射率的 角度距有关。对应的体表出射光Φ(r,τ_i)：

其中，J₀，…，J₃为一系列常数。

步骤二，搭建线性关系——基于有限元离散的SP3扩散方程方法，针对每一个波长搭建边 界测量值与未知光源的关系模型；

在公式(6)和(7)中，S(r,τ_i)表示的是荧光光源分布，是需要求解的未知量，为了使 用有限元法对SP3方程及其边界条件进行求解，首先将SP3的两个方程写成对应的弱解形 式：

其中，Ψ(r,τ_i)是测试函数，公式(8)和(9)两个边界条件也要写为弱解形式：

然后利用格林公式(Green’s Formula)将边界条件的弱解形式融入到SP3方程的弱解形 式中。

其中e₁、e₂、e₃和e₄分别为：

$e_1=\frac{1}{W}(\frac{D_1(1+8C_2)}{8{\mu }_{a3}}-\frac{(1+B_2)(1+2A_1)}{14{\mu }_{a3}})$

$e_2=\frac{1}{W}(\frac{(1+B_2)(1+8C_2)}{56{\mu }_{a3}}-\frac{D_1(7+24A_2)}{24{\mu }_{a3}})$

$e_3=\frac{1}{W}(\frac{(1+B_1)(1+8C_2)}{24{\mu }_{a1}}-\frac{D_2(7+2A_1)}{2{\mu }_{a1}})$

$e_4=\frac{1}{W}(\frac{D_2(1+8C_1)}{8{\mu }_{a1}}-\frac{(1+B_1)(7+24A_2)}{72{\mu }_{a1}})$

上式中，W具有如下形式：

$W=\frac{(1+B_1)(1+B_2)}{21{\mu }_{a1}{\mu }_{a3}}-\frac{D_1{D}_2}{{\mu }_{a1}{\mu }_{a3}}$

进行网格剖分，将光源分布S(r,τ_i)写成只与谱段相关的S(τ_i)，选择相应基函数，合成整 体的刚度矩阵M_ij，将SP3方程离散为如下形式的矩阵方程：

其中B为N×N大小的矩阵，和是是两个分块子矩阵，可表示为：

这里需要建立出表面出射光Φ(τ_i)和未知荧光光源分布S(τ_i)之间的线性关系。在实验过 程中，由于只能采集到生物体表面的光学信号，因此，只保留Φ(τ_i)中边界上节点值，去掉矩 阵中不包含边界节点的行，得到方程：

上述方程建立了未知荧光光源分布和边界出射光之间的关系，即：

Φ(τ_i)＝A(τ_i)S(τ_i)

在多光谱情况下，谱段的能量比ω(τ_i)通过预先的谱分析测得，令S代表所有谱段的光源 的总能量，即p为谱段个数。将所有谱段的光源和边界测量值进行整合， 多光谱情况下未知荧光光源分布和边界出射光之间的关系如下：

A_mulS＝Φ_mul

此时，

$A_{mul}=\left(\begin{array}{c}\omega \left({\tau }_1\right)A\left({\tau }_1\right)\\ \omega \left({\tau }_2\right)A\left({\tau }_2\right)\\ .\\ .\\ .\\ \omega \left({\tau }_P\right)A\left({\tau }_P\right)\end{array}\right),{\Phi }_{mul}=\left(\begin{array}{c}\Phi \left({\tau }_1\right)\\ \Phi \left({\tau }_2\right)\\ .\\ .\\ .\\ \Phi \left({\tau }_p\right)\end{array}\right)$

步骤三，共享先验估计——基于经验贝叶斯极大似然函数，推断表征多光谱相关的参数 α；

在实际问题中，出射光测量值Φ(τ_i)通常含有噪声，若噪声服从均值为零，方差为的 高斯随机分布，则Φ(τ_i)关于荧光光源分布S(τ_i)和α₀的极大似然函数表示为：

$p\left(\Phi \right({\tau }_i\left)\right|S\left({\tau }_i\right),{\alpha }_0)={\left(\frac{{\alpha }_0}{2\pi }\right)}^{\frac{M}{2}}\exp \{-\frac{{\alpha }_0}{2}||\Phi \left({\tau }_i\right)-A\left({\tau }_i\right)S\left({\tau }_i\right)|{|}_2^2\}$

为了表示谱段间关系，引入一种多任务的思想，将整个光源作为一个整体任务T。其中， S(τ_i)是T中的第i个任务模型，用S_j(τ_i)表示为S(τ_i)的第j个组成部分。若S(τ_i)服从均值为零的 高斯先验分布，参数α₀服从Gamma先验分布，即有：

$p\left(S\right({\tau }_i\left)\right|\alpha ,{\alpha }_0)={\Pi }_{j=1}^{N}N\left(S_j\right({\tau }_i\left)\right|0,{\alpha }_0^{-1}{\alpha }_j^{-1})$

p(α₀|a,b)＝Ga(α₀|a,b)

其中，α＝(α₁,…α_j,…α_N)^T是超参数，用来表示先验信息。那么，在所有的谱段中，α_j的 取值相同，即用α表征各个谱段之间的相关性。为了估计超参数α，使用经验贝叶斯的方法， 得到α的边界最大似然函数，即有：

$\left(\begin{array}{c}L=\underset{i=1}{\overset{p}{\Sigma }}\log p\left(\Phi \right({\tau }_i\left)\right|\alpha )\\ ={\Sigma }_{i=1}^p\log (\int P\left(\Phi \left({\tau }_i\right)|S\left({\tau }_i\right),{\alpha }_0\right)P\left(S\left({\tau }_i\right)|\alpha ,{\alpha }_0\right)P\left({\alpha }_0|a,b\right)){\mathrm{dS}}_{\tau i}{\mathrm{d\alpha }}_0\\ =-\frac{1}{2}{\Sigma }_{i=1}^p\left[\left(n_i+2a\right)\log \left(\Phi {\left({\tau }_i\right)}^T{B}_i^{-1}\Phi \left({\tau }_i\right)+b\right)+\frac{\log \left|B_i\right|}{2}\right]+const\end{array}\right)$

其中，B_i＝I+A(τ_i)Λ^-1A(τ_i)^T，Λ^-1＝diag(α₁,α₁,…,α_N)。

步骤四，未知光源重建——根据已知的超参数的估计和边界测量值Φ_mul，利用极大似 然函数重建荧光光源。

得到估计的后，再将所有谱段的光源和边界测量值进行整合，则未知光源S的大似然函 数可表示为：

$\left(\begin{array}{c}P\left(S\right|{\Phi }_{mul},\hat{\alpha })=\int P\left(S\right|{\Phi }_{mul},\hat{\alpha },{\alpha }_0)P\left({\alpha }_0\right|a,b){\mathrm{d\alpha }}_0\\ =\frac{\Gamma \left(a+\frac{N}{2}\right){\left[1+\frac{2}{2b}{\left(S-{\mu }_{mul}\right)}^T{\Sigma }_{mul}^{-1}\left(S-{\mu }_{mul}\right)\right]}^{-\left(a+\frac{N}{2}\right)}}{\Gamma \left(a\right){\left(2\pi b\right)}^{N/2}{\left|{\Sigma }_{mul}\right|}^{1/2}}\end{array}\right)$

这里，

${\mu }_{mul}={\Sigma }_{mul}{A}_{mul}^T{\Phi }_{mul}$

$\Sigma mul={(A_{mul}^T{A}_{mul}+\Xi )}^{-1}$

式中，多任务模型框架如图1所示。

本发明基于多任务学习方法，首先采用高阶近似模型模型光在生物组织中的传输规律， 并基于多任务学习方法探索多光谱之间的内在相关性，最后在此基础上实现了荧光光源的 三维重建。从重建结果来看，本发明不但将小鼠体内的荧光光源进行了准确重建定位，而 且计算效率也得到了极大提高。

附图说明

图1为多任务模型框架结构图；

图2为数字鼠胸部测量值采集区示意图；

图3为肺部中光源的重建结果图(a)为肺部中光源的重建位置，(b)-(d)表示对应XCT的 BLT重建结果的冠状切片。在(b)-(d)中，用圈标记示意真实光源的位置；

图4为基于L1正则化重建方法的光源重建结果(a)-(c)表示对应XCT的BLT重建结果的 冠状切片，用圈标记示意真实光源的位置；

图5为肝脏中光源重建结果图(a)和(b)是利用micro-CT的重建切片。(c)-(f)对应基于本章 算法重建的光源分布，(g)-(j)对应使用L1范数正则化方法的重建结果。(c)和(i)是y＝51mm时 的冠状切片，(f)和(j)为z＝12.97mm的横截面。图中用圈标记真实光源位置。

图6为两种算法的重建结果。(a)和(c)对应本算法重建的结果，(b)和(d)对应基于L1正则 化算法的重建结果。(a)和(b)为y＝51.05mm的冠状面，(c)和(d)为z＝15.97的横截面。真实光 源位置用圈标记。

图7为肝脏中光源的重建结果(a)是micro-CT数据的冠状切片。(b)和(d)分别为IRIRM 和L1-RMRM重建的结果。(c)和(e)为两种算法在XCT上重建结果。用圈标记真实光源位置。

图8光源的重建结果(a)为基于贝叶斯压缩感知方法的重建结果，(b)为本算法重建结果。 用圈标记真实光源位置。

具体实施方式

下面结合附图和实施例对本发明做更详细的说明。

为了检验本发明提出方法的有效性，采用了非匀质的数字鼠仿体进行仿真实验，共进行 了两项实验。在实验中，认为生物组织的光学特性参数已知，详细参数如表1所示。本发明 采用有限元方法产生前向边界测量数据。为有效地避免“逆行为”，在产生前向数据时，采 用SP3方程产生前向数据，前向数据的网格含有42342个节点和208966个四面体单元，重建 时所需要的网格含有20196个节点和108086个四面体。考虑到光谱分布，实验中选择的两个 谱段分别为620nm和700nm。为真实模拟荧光数据采集，仅使用了数字鼠胸部的测量数据(共 1311个边界测量点)，如图2所示。

表1生物组织光学特性参数

实验一：

第一组：第一组试验中，将半径大小为1.5mm的球形光源放置在小鼠的肺中，球心位 置为(14.5,44.0,12.8)mm，光源位置距生物体表面距离为4.8mm。

使用本发明提出的重建算法进行荧光光源重建，并对重建结果进行分析，结果如图3 所示。可以看到，利用本发明提出的方法可以准确的重建荧光光源，重建光源的中心位置 为(14.69,44.17,12.91)mm，与真实光源的绝对距离为0.28mm。

进一步，为了验证本章算法的有效性，与基于L1正则化的重建方法进行了对比，基 于L1正则化的重建方法重建结果如图4所示。利用L1正则化重建方法重建的光源中心 位置为(14.59,43.63,12.88)mm，与真实光源中心距离为0.39mm。

另外，利用L1正则化重建方法，重建耗费的计算时间为1773.9秒，而利用本章算法 重建，耗费的计算时间仅为0.67秒。

第二组：在第二组实验，在实验中，保持光源的大小不变，将光源放置于肝脏中，光 源中心位置为(22.0,51.0,13.0)mm，距离生物体表面7.06mm。

然后，用本发明的方法进行光源重建，结果如图5所示。可以看到，利用本章提出的 方法可以准确的重建荧光光源，重建光源的中心位置为(22.37,50.68,13.15)mm，与真实光 源的绝对距离为0.5612mm。

为了验证算法有效性，依然与基于L1正则化的重建方法进行了对比，结果如图6所 示。利用L1正则化重建方法重建的光源中心位置为(23.01,50.06,14.87)mm，与真实光源中 心距离为2.3239mm。

因此可以看出，基于L1正则化的重建方法在重建的荧光光源位于生物组织中较深位置 时，重建效果较差，误差较大；而本章提出的方法可以准确地定位荧光光源。同时，基于 L1正则化方法重建的计算时间为1738.8秒，而本文的算法的计算时间仅为0.936秒。

将以上两组实验进行总结性量化比较，将本章提出的基于多任务压缩感知方法的重建 算法简写为IRIRM，将基于L1范数正则化的重建算法简写为L1-RMRM。对比结果如表2 所示。

表2两种算法的量化比较

实验二：

为了进一步测试本算法的有效性，开展了第二项实验：双光源的实验。

第一组：在第一组实验中，首先在肝脏和肺中各放置一个半径为1.5mm的球形光源， 其中心位置分别为(14,51,16)mm和(21,51,16)mm，光源到生物体体表的距离分别为 5.56mm和5.82mm。利用本发明的方法和基于L1正则化的方法分别进行重建，重建结果 如图7所示。从图7中可以看到，两种算法都可以准确的重建荧光光源。

第二组：在第二组实验中，将两个光源均放置于肝脏中，利用上述两种算法进行光源 重建，重建结果如图8所示。从图8可以看到，基于L1正则化的重建方法仅重建出了一 个光源，而本发明提出的算法可以准确的重建两个光源，这进一步证明了本方法的有效性。

将以上两组双光源的实验进行总结性量化比较，对比结果如表3所示。

表3双光源两种算法的量化比较

在本发明中，利用小鼠形状的仿体，先后在肺部和肝脏中设置了单光源进行仿真实验， 并将本发明提出的算法与基于L1正则化的重建算法进行对比。进一步地，又开展了双光源 实验，分别用本发明的重建方法和基于L1正则化的重建算法进行光源重建，并对比实验结 果。实验结果表明，本算法不但可以在小鼠体内对荧光光源进行准确重建定位，而且可以 极大地提高计算效率，进一步验证了本算法的有效性。