生物发光断层成像的仿真研究

刘 勇，齐树波，方勇军，陈 锋

0 引言

与传统医学成像模态，例如X线计算机断层成像（X-ray computed tomography，XCT）不同，分子成像能在特异性分子探针的帮助下，在分子水平上反映生物体生理或病理变化，在分子生物学与临床医学之间架起了相互连接的桥梁[1]。根据成像原理不同，分子影像技术主要可分为单光子发射成像（single photon emission computed tomography，SPECT）、正电子发射成像（positron emission tomography，PET）、磁共振成像（magnetic resonance imaging，MRI）、超声成像（ultrasound）以及光学成像（optical imaging，OI）等。相对于其他分子成像技术，光学分子成像技术具有无电离辐射、对生物体无害、灵敏度高且其成本相对较低等优点[1]，目前已发展成为一种理想的成像方式。

作为光学分子成像的一个重要分支，近年来，生物发光断层成像技术（bioluminescence tomography，BLT）得到了快速发展[2-5]。作为一种新的分子断层成像模式，BLT能够以非侵入的方式，在细胞和分子水平研究或监测生理和病理过程，因而有望为疾病早期诊断、基因治疗和药物研发等研究提供有力工具。然而，BLT作为一种基于扩散光的成像技术，对比XCT等传统成像技术，具有较低的成像空间分辨率。主要原因为：由于光在生物组织中的强散射特性，一方面导致了光在生物组织传播过程中的强烈衰减；另一方面，也导致了光在生物体内的传播不是沿直线方向进行。上述2个原因致使BLT的求解极具病态性（ill-posed），进一步影响了其生物医学应用。

在本文中，基于数值仿真实验，我们验证了BLT成像的可行性。同时，我们将稀疏理论应用于重建过程，以克服BLT成像的病态性。实验结果表明，基于多角度、非接触BLT成像系统（multi-view non-contactBLT imaging system），我们能够解析自发光源在成像体内的分布，定位误差小于1mm，这将有望为疾病早期诊断及药物研发等研究提供有力工具。

1 方法

1.1 BLT前向数学模型

对于BLT成像技术而言，准确描述光在生物组织中的传输过程是成像的首要问题。根据输运理论，光在生物组织中的传播可以通过辐射传输方程（radiative transferequation，RTE）来描述。直接对辐射传输方程求解虽然可以较为准确地描述光子在组织中的传播过程，但计算量较大。为了在实际中更好地进行模拟，需要对辐射传输方程进行一定的简化。对于BLT成像而言，通常可采用连续波模式下的扩散方程（diffusion equation，DE）[6]对其近似，如式（1）所示：

其中，Ω 描述了成像物体；S（r）为自发光光源；Φ（r）为光场分布；μa（r）与Φ（r）分别描述了组织的吸收和扩散系数。

根据微分方程理论，扩散方程（1）的求解需要有合适的边界条件。对于光学成像而言，常用的边界条件为外插边界条件和Robin边界条件。在本文中，我们使用Robin边界条件[6]，如式（2）所示：

考虑到有限元方法（finiteelement method，FEM）[7]可以有效地处理复杂几何形状且各向异质的成像问题，因而，在本文中，我们使用有限元方法对式（1）进行求解。当重建区被离散成小的网格后，式（1）可被改写为如下的线性方程：

式（3）中：

式中，φi（r）、φj（r）为相应的试验函数，可从一个合适的试验空间中选取，且必须满足平方可积条件。考虑到在式（3）中，矩阵A通常是正定的，于是，式（3）可进一步改写为：

式中，W为权重，用于映射成像体内未知的生物发光源S到已知的成像表面测量值Φ。

1.2 BLT逆向问题求解

BLT的逆向问题是从式（3）中求解未知的生物发光源S。然而，如前所述，考虑到：（1）光子在生物组织中的传播是高度散射的；（2）用于描述光子在成像对象体内传播的数学模型（前向模型）是近似的。因此，对比XCT或MRI成像，BLT的重建问题有更大的病态性。此外，测量数据Φ通常包含噪声，这进一步复杂了BLT的成像问题。于是，直接对式（3）进行求解不太可行。考虑到在大多数生物学应用中，生物发光源在成像体内的分布通常具有稀疏性或可稀疏性。为了提高成像的空间分辨率，在本文中我们拟基于稀疏理论进行求解[8-10]：

式（7）中，ε用于描述测量数据包含的噪声水平；S是待重建的纳米荧光团的密度。

2 仿真实验设置及实验结果

2.1 BLT成像系统

在本文中，我们所模拟的成像系统是多角度、非接触的BLT成像系统（multi-view non-contact BLT imagingsystem）。该系统主要包括光子探测器（CCD）、旋转台及计算机等设备，如图1所示。

图1 多角度、非接触BLT成像系统

2.2 仿真实验设置

在数字仿真实验中，一个直径3.0 cm、高5.0 cm的圆柱体被使用作为仿体，如图2所示。仿体内的光学参数被配置为 μa=0.3 cm-1、μ′s=10.0 cm-1。这与小鼠组织的平均光学参数类似。此外，1个直径0.4 cm、高0.6 cm的小圆柱体被浸入在仿体中，用以模拟生物发光源。

图2 仿真实验参数设置示意图

直径3.0 cm、高5.0 cm的圆柱体被使用作为仿体；仿体内的光学参数为 μa=0.3 cm-1、μs′=10.0 cm-1；直径0.4 cm、高0.6 cm的圆柱体被使用作为生物发光源。

为了提高BLT成像的空间分辨率，在BLT成像过程中，我们总计获取了6张自发光投影图像，等间隔 60°。

2.3 实验结果

对于BLT重建，首先，图2中的圆柱仿体被离散为5 675个节点及29 317个四面体单元，如图3所示。其次，重建被执行。在重建过程中，我们使用稀疏重建算法求解式（7）。对于每个投影，检测点分布在圆柱仿体140°×4.6 cm的范围内。BLT的重建结果如图4所示。

图3 圆柱形仿体的有限元离散结果

图4 BLT重建结果

为了定量评估BLT重建性能，我们计算了重建的定位误差，见表1。这里，定位误差被定义为：实际生物发光源的中心与BLT重建发光源的中心最大值之间的距离。

表1 BLT重建的定位误差

实验结果显示，基于多角度、非接触的BLT成像系统，结合稀疏重建算法，我们能够在较少（6个）的角度下，获取较好的BLT重建结果，如图4及表1所示。

3 结语

基于多角度、非接触的BLT成像系统，结合稀疏重建算法，本文解析了生物发光源在成像体内的三维分布信息。数值仿真实验结果显示，定位误差小于1mm。在今后的研究中，我们将重点搭建BLT成像系统，进行物体仿体及在体小动物实验。

[1] Ntziachristos V，Ripoll J，WANG Lyu，etal.Looking and listening to light：the evolution of whole-body photonic imaging[J].Nature Biotechnology，2005，23（3）：313-320.

[2] WANG Ge，CONG Wen-xiang，Durairaj K，et al.In vivo mouse studies with bioluminescence tomography [J].Optics Express，2006，14（17）：7 801-7 809.

[3] Lv Y，Tian J，Cong W，et al.A multilevel adaptive finite element algorithm for bioluminescence tomography[J].Optics Express，2006，14（18）： 8 211-8 223.

[4] Gu X，ZhangQ，Larcom L，etal.Three-dimensional bioluminescence tomography with model-based reconstruction[J].Optics Express，2004，12（17）： 3 996-4 000.

[5] Lv Y，Tian J，Cong W，et al.Spectrally resolved bioluminescence tomography with adaptive finite element analysis：methodology and simulation[J].Physics inmedicine and biology，2007，52（15）：4 497.

[6] Schweiger M，Arridge S R，Hiraoka M，et al.The finite element method for the propagation of light in scattering media：boundary and source conditions[J].Medical Physics，1995，22（11）：1 779-1 792.

[7] Joshi A，Bangerth W，Sevick-Muraca E.Adaptive finite element based tomography for fluorescence optical imaging in tissue[J].Optics Express，2004，12（22）：5 402-5 417.

[8] Candès E J，Wakin M B.An introduction to compressive sampling[J].IEEE Signal Processing Magazine，2008，25（2）：21-30.

[9] Elad M，Aharon M.Image denoising via sparse and redundant representationsover learned dictionaries[J].IEEE Transactionsan Image Processing，2006，15（12）：3 736-3 745.

[10] Tropp JA，Gilbert A C.Signal recovery from random measurements via orthogonal matching pursuit[J].IEEE Transactionsan Information Theory，2007，53（12）：4 655-4 666.