CN103271723B - 一种生物发光断层成像重建方法 - Google Patents

Abstract

本发明公开了一种生物发光断层成像重建方法，本发明利用生物发光断层成像数据采集平台获取小动物体内的生物发光光源在体表的信号分布情况和MicroCT系统获取小动物体的解剖结构信息，然后基于辐射传输方程的三阶简化球谐近似模型和有限元网格构建系统方程，采用基于压缩感知的稀疏正则化方法求解系统方程进而得到精确的生物体内光源的定位信息。利用本发明的生物发光断层成像的重建方法可以通过二维的小动物体表的光源分布和解剖结构信息反演出小动物体内的三维光源分布和定位信息。

Description

一种生物发光断层成像重建方法

技术领域

本发明涉及光学成像技术领域，具体涉及到一种生物发光断层成像的方法，本发明基于三阶简化球谐波近似模型和稀疏正则化方法，可应用于疾病早期诊断、基因治疗和药物研发等领域。

背景技术

近年来，由于在灵敏性、安全性、操作性和成本等方面比X光成像、磁共振成像和核素成像方式具有优势，光学成像技术得到了快速的发展，在生物研究中被广泛应用，取得了显著的科研成果。其中，生物发光断层成像(Bioluminescence tomography，BLT)作为一种领先的光学成像技术，通过高灵敏度CCD相机采集光学信号，结合生物组织光传输理论和自由空间光传输理论，采用重建算法获取生物体内光源所在的三维空间位置信息或能量信息。BLT能够克服BLI平面成像的局限，提供生物发光光源精确的空间位置和能量信息，加上本身成像信噪比相对较高，得到了国内外学者和研究人员的广泛关注和研究，目前已开始应用于基因表达、肿瘤检测和药物研发等小动物实验和预临床实验。近几年，生物发光断层成像理论得到了迅速发展和广泛关注，前向光传输模型的求解是BLT逆向光源重建的基础，快速有效地重建方法是BLT目前研究的重点和核心。

目前在光学成像领域中，辐射传输方程(Radiative transportequation，RTE)被公认为一种较为准确的生物组织光传输模型，但作为一个复杂的积分-微分方程，RTE方程的求解时间代价比较大。为此，人们提出了RTE方程的一些高阶或低阶近似模型，这些近似模型以降低求解精度为代价，大大地减少了求解所需的计算量。中国科学院自动化研究所在其专利申请文件″一种用于自发荧光成像的空间加权的有限元重建方法″(申请号200810102314.7，申请日2008.03.20，)中提出了一种基于扩散方程的自发荧光断层成像重建方法。但是作为RTE方程的低阶近似，扩散近似模型适用于生物组织高散射低吸收的情况，无法精确描述光在高吸收介质中的传输。由于生物体中很多组织对光的吸收较强，对于吸收系数较高的生物组织误差非常大，这一缺点限制了扩散近似模型在全身成像中的应用。而且扩散近似要求光源和探测器之间的相隔光子10个自由程以上，这限制了扩散近似模型在小动物成像中的应用。

生物发光断层成像的逆向重建问题是根据生物体表面获得的测量数据采用合适的重建算法反演生物体内发光光源的位置及定量信息。但是由于BLT是早期成像，在细胞和分子水平下，光源信号较弱，而且受到生物组织散射吸收的影响，在体表被CCD相机接受的信号会有衰减和噪声，再加上生物发光断层成像的逆向重建问题是一个严重病态的不适定问题，这使得开发稳定高效的重建算法成为生物发光断层成像中的充满挑战性的问题。中国科学院自动化研究所在其专利申请文件″一种基于水平集检测目标体内光源分布的方法″(申请号200910091634.1，申请日2009.08.26)中提出了一种基于/₂范数正则化的自发荧光断层成像重建方法。但是传统的/₂范数正则化方法容易产生平滑解，重建光源能量分布过于平滑，重建结果分布范围较大，远远超出光源本身大小，而且同时在BLT逆问题重建中在真实光源附近产生多个伪光源。

现有专利或文献中的BLT重建方法大部分采用扩散近似方程和基于/₂范数的正则化方法，最新的文献中虽然有在前向模型中采用较为精确的RTE方程的高阶球谐波近似模型或者在重建方法中采用/₁正则化算法的，但是如何结合前向和逆向问题的特点，同时兼顾前向模型和逆向重建方法的准确性和效率依然是一个挑战性的问题。

发明内容

鉴于现有技术的不足，本发明旨在于提供一种基于三阶简化球谐波近似模型和稀疏正则化的生物发光断层成像重建的方法，提高了前向模型的准确性和逆向重建的稀疏性，实现了生物发光断层成像的准确快速重建。

本发明的实现思路具体是：根据生物体的组织特异性，采用三阶简化球谐波近似模型代替扩散近似模型，保证前向模型构建的精度和效率；基于自由空间传输理论，建立探测器采集的光学数据和生物体三维表面数据的映射关系；基于hp自适应有限元方法建立内部未知光源和表面已知测量值之间的线性关系；基于/_1/2正则化的加权内点法，结合生物体内荧光光源的稀疏性，在逆向重建中同时进行定位和定量的三维重建，实现了生物发光断层成像的准确高效重建。

为了实现上述目的，本发明采用的技术方案如下：

一种生物发光断层成像重建方法，所述方法包括以下步骤：

(1)数据采集与预处理，其中，在已产生荧光光子的生物体采集多角度BLT荧光和白光数据，以及生物体的360度CT投影数据，并进行数据格式转换、除去噪声点、荧光和白光数据叠加、几何校正数据预处理；

(2)生物体组织结构信息获取，其中，设置感兴趣区域，利用Micro-CT软件进行数据重建，重建后的CT数据进行器官分割、融合形成完整的生物体数据，对分割后的数据进行生物体表面数据的提取和内部体数据的剖分，并保存两组数据；

(3)生物体表荧光数据获取，其中，对步骤(1)获取的预处理后光学数据与步骤(2)获取的生物体解剖结构，应用自由空间光传输理论，将CCD相机获取的二维光学数据映射到生物体表面，获取生物体表面的三维荧光数据分布；

(4)前向光传输模型构建，其中，利用三阶简化球谐近似模型描述生物组织内的光传输过程，推导出三阶简化球谐近似方程及其边界条件的边界测量方程，通过生物发光断层成像获取生物体的解剖结构、各生物组织的光学参数等先验信息，基于hp自适应有限元方法建立内部未知光源和表面已知测量值之间的线性关系，得到如下系统矩阵方程：

J^+，m＝GS；

其中J^+，m为生物体表面上光通量密度的测量值，G为系统矩阵，S为生物体内的光流密度分布；

(5)稀疏正则化目标函数建立，其中，由于生物发光断层成像中生物体内光源分布具有稀疏性，基于压缩感知理论，将上述线性矩阵方程表示为如下基于0范数的稀疏正则化形式：

min||GS-J^+，m||²+λ||S||₀

其中λ为正则化参数；其次，基于稀疏近似，构建BLT逆问题目标函数为如下的/_1/2正则化方程：

min||GS-J^+，m||²+λ||S||_1/2；

(6)目标函数的优化求解，其中，采用加权内点法将/_1/2正则化目标函数转化为重赋权的/₁正则化极小化问题，然后内点法求解极小化问题获取生物体内荧光光源的三维定位定量信息；

(7)三维重建结果显示。

需要说明的是，其特征在于，所述步骤(1)中的荧光数据和CT投影数据的初始采集角度相同，且至少采集四个角度间隔为90度的荧光和白光数据，其中，CT数据采集为间隔0.75度，480个投影数据。

需要说明的是，其特征在于，所述步骤(3)中的生物体表荧光数据获取，需要进行荧光图片和白光图片叠加，以及基于CT数据的标记点进行荧光数据和CT坐标系的对应。

需要进一步说明的是，其特征在于，所述目标函数的优化求解为：

(1)参数初始化，令t＝0，最大迭代次数为K，初始化S⁰＝(1，...，1)^T；

(2)目标函数转化，基于内点法求解：

$S^{t+1}=\underset{{0<S<S}^{\sup }}{\min }{\left|\right|\mathrm{GS}-J^{+,m}\left|\right|}^2+\mathrm{\&lambda;}\left|S\right|/(\sqrt{\left|S^t\right|}+c),$

将其转化为：

$S^{t+1}=\underset{{0<S<S}^{\sup }}{\min }{\left|\right|\mathrm{GS}-J^{+,m}\left|\right|}^2+\mathrm{\&lambda;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_i-\frac{1}{w}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}[\log (u_i+s_i)+\log (u_i-s_i)];$

(3)内部参数设置，设置目标对偶距ζ＞0，w＝1/λ，S＝(1，...，1)^T；

(4)预处理共轭梯度法求解牛顿系统方程 $H\left[\begin{array}{c}\mathrm{\&Delta;S}\\ \mathrm{\&Delta;u}\end{array}\right]=-g,$ 其中H为目标函数的Hessian矩阵；

(5)采用回溯线搜索方法计算步长a；

(6)更新迭代解(S，u)＝(S，u)+a(ΔS，Δu)；

(7)计算对偶可行点υ和对偶距ζ；

(8)如果ζ/G(υ)＜1e-6则退出；否则更新w；

(9)得到S^t+1，并令t＝t+1；

(10)当t＜K时，转到步骤2；当t＝K时，输出目标解S^t。

本发明有益效果在于：

1、本发明采用三阶简化球谐波模型构建BLT前向模型能够克服扩散近似模型中存在的缺点，准确的描述光在高吸收介质中的传输过程，其精度不局限于较小的生物体，而且采用三阶近似模型在保证精度的同时兼顾了运算效率。

2、本发明不使用多光谱数据及先验可行区等先验知识，避免了采集多组重建数据和选择可行区，有效的降低了光源重建所需的时间，提高了重建效率。

3、本发明采用/_1/2正则化算法，相对于目前广泛采用的/₁正则化算法具有更好的稀疏性，更可靠的鲁棒性，能够得到更好的重建精度。

附图说明

图1为生物发光断层成像重建方法的整体流程图；

图2为/_1/2正则化方法的实施图。

具体实施方式

下面将结合附图对本发明进行描述，

如图1所示，本发明为一种生物发光断层成像重建方法，所述方法包括以下步骤：

(1)数据采集与预处理；

首先对小动物腹腔或静脉注射底物荧光素，在生物体内荧光素酶与底物荧光素发生反应产生荧光光子。固定生物体，在暗室利用BLT系统采集多角度荧光和白光等光学数据。采用CT系统扫描生物体，获取扫描数据及几何校正数据。

将采集到的BLT荧光数据进行格式转换及背景去噪处理，然后将得到的荧光数据和白光数据进行叠加处理，得到新的叠加数据。

对采集到的CT数据进行预处理。处理暗电流和空扫数据，去除坏点和去噪，然后进行几何校正，预处理后的数据用于CT数据重建。

(2)生物体组织结构信息获取；

将获取的CT数据通过专业软件进行格式转化，然后分割出骨骼肌肉等各个器官，对分割数据进行高斯平滑，还有进行各个器官的融合，形成一个完整的生物体结构数据，其中还包括整体的平滑以及去空洞等后处理操作。

重建后的CT数据进行表面数据提取和体数据剖分并保存。

(3)生物体表荧光数据获取；

将各角度光学数据的灰度值与空间位置对应，即荧光数据与CT坐标系间的对应。然后根据自由空间光传输理论，将多角度二维光学测量数据结合生物体三维组织结构信息进行数据映射(2D-3D数据映射)，得到生物体表面的三维荧光数据分布。

(4)前向光传输模型的构建；

本发明方法采用三阶简化球谐波模型(SP₃)构建BLT前向模型。SP₃方程如下：

SP₃边界条件如下：

进一步推导出SP₃的边界测量方程如下：

为了方便描述重建算法，SP₃方程的两个子方程可以采用如下的统一范式来描述

为建立生物体内发光光源分布与生物体表测量信息之间的线性关系，将三阶简化球谐近似模型方程转化为线性矩阵方程，建立系统方程。

根据边界测量离散点，在hp自适应有限元框架下，得到SP₃的边界测量方程如下：

其中β₁和β₂可以通过SP₃的边界测量方程得到，G₁和G₂为移除相应行列后的矩阵。由于存在噪声的影响，所以使用边界测量数据J^+，m取代J^+，b，得到系统方程：

J^+，m＝GS。

(5)稀疏正则化目标函数建立；

基于压缩感知理论，结合生物体内发光光源稀疏分布的特性，建立稀疏正则化的目标函数如下：

$\underset{{0<S<S}^{\sup }}{\min }{\left|\right|\mathrm{GS}-J^{+,m}\left|\right|}^2+\mathrm{\&lambda;}{\left|\right|S\left|\right|}_0$

其中S^sup为光源能量的上界，λ为正则化参数。

由于基于0范数的稀疏正则化目标函数是一个NP组合优化问题，因此需要对其进行转化求解，基于稀疏近似，构建BLT逆问题目标函数为如下的/_1/2正则化方程：

$\underset{0<S<S^{\sup }}{\min }{\left|\right|\mathrm{GS}-J^{+,m}\left|\right|}^2+\mathrm{\&lambda;}{\left|\right|S\left|\right|}_{1/2}.$

(6)目标函数的优化求解；

采用加权内点法将/_1/2正则化目标函数转化为重赋权的/₁正则化极小化问题，然后内点法求解极小化问题获取生物体内荧光光源的三维定位定量信息。

(7)重建结果的三维显示；

基于步骤6的重建结果，我们得到了生物体内发光光源的位置和能量信息，为了更好的展现重建效果，结合步骤2获取的生物体组织结构信息，我们使用专业软件对重建结果进行三维立体显示，获取生物体各个组织、内部光源及重建光源的立体效果图。

如图2所述，所述目标函数的优化求解为：

(1)参数初始化，令t＝0，最大迭代次数为K，初始化S⁰＝(1，...，1)^T；

(2)目标函数转化，基于内点法求解：

$S^{t+1}=\underset{{0<S<S}^{\sup }}{\min }{\left|\right|\mathrm{GS}-J^{+,m}\left|\right|}^2+\mathrm{\&lambda;}\left|S\right|/(\sqrt{\left|S^t\right|}+c),$

将其转化为：

$S^{t+1}=\underset{{0<S<S}^{\sup }}{\min }{\left|\right|\mathrm{GS}-J^{+,m}\left|\right|}^2+\mathrm{\&lambda;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_i-\frac{1}{w}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}[\log (u_i+s_i)+\log (u_i-s_i)];$

(3)内部参数设置，设置目标对偶距ζ＞0，w＝1/λ，S＝(1，...，1)^T；

(4)预处理共轭梯度法求解牛顿系统方程 $H\left[\begin{array}{c}\mathrm{\&Delta;S}\\ \mathrm{\&Delta;u}\end{array}\right]=-g,$ 其中H为目标函数的Hessian矩阵；

(5)采用回溯线搜索方法计算步长a；

(6)更新迭代解(S，u)＝(S，u)+a(ΔS，Δu)；

(7)计算对偶可行点υ和对偶距ζ；

(8)如果ζ/G(υ)＜1e-6则退出；否则更新w；

(9)得到S^t+1，并令t＝t+1；

(10)当t＜K时，转到步骤2；当t＝K时，输出目标解S^t。

对于本领域的技术人员来说，可根据以上描述的技术方案以及构思，做出其它各种相应的改变以及变形，而所有的这些改变以及变形都应该属于本发明权利要求的保护范围之内。

Claims (4)

1.一种生物发光断层成像重建方法，其特征在于，所述方法包括以下步骤：

(1)数据采集与预处理，其中，在已产生荧光光子的生物体采集多角度BLT荧光和白光数据，以及生物体的360度CT投影数据，并进行数据格式转换、除去噪声点、荧光和白光数据叠加、几何校正数据预处理；

(2)生物体组织结构信息获取，其中，设置感兴趣区域，利用Micro-CT软件进行数据重建，重建后的CT数据进行器官分割、融合形成完整的生物体数据，对分割后的数据进行生物体表面数据的提取和内部体数据的剖分，并保存两组数据；

(3)生物体表荧光数据获取，其中，对步骤(1)获取的预处理后光学数据与步骤(2)获取的生物体解剖结构，应用自由空间光传输理论，将CCD相机获取的二维光学数据映射到生物体表面，获取生物体表面的三维荧光数据分布；

(4)前向光传输模型构建，其中，利用三阶简化球谐近似模型描述生物组织内的光传输过程，推导出三阶简化球谐近似方程及其边界条件的边界测量方程，通过生物发光断层成像获取生物体的解剖结构、各生物组织的光学参数先验信息，基于hp自适应有限元方法建立内部未知光源和表面已知测量值之间的线性关系，得到如下系统矩阵方程：

J^+,m＝GS

其中J^+,m为生物体表面上光通量密度的测量值，G为系统矩阵，S为生物体内的光流密度分布。

(5)稀疏正则化目标函数建立，其中，由于生物发光断层成像中生物体内光源分布具有稀疏性，基于压缩感知理论，将上述线性矩阵方程表示为如下基于0范数的稀疏正则化形式：

min||GS-J^+,m||²+λ||S||₀

其中λ为正则化参数；其次，基于稀疏近似，构建BLT逆问题目标函数为如下的l_1/2正则化方程：

min||GS-J^+,m||²+λ||S||_1/2；

(6)目标函数的优化求解，其中，采用加权内点法将l_1/2正则化目标函数转化为重赋权的l₁正则化极小化问题，然后内点法求解极小化问题获取生物体内荧光光源的三维定位定量信息；

(7)三维重建结果显示。

2.根据权利要求1所述的生物发光断层成像重建方法，其特征在于，所述步骤(1)中的荧光数据和CT投影数据的初始采集角度相同，且至少采集四个角度间隔的荧光和白光数据，其中，CT数据采集为间隔0.75度，480个投影数据，角度为90度。

3.根据权利要求1所述的生物发光断层成像重建方法，其特征在于，所述步骤(3)中的生物体表荧光数据获取，需要进行荧光图片和白光图片叠加，以及基于CT数据的标记点进行荧光数据和CT坐标系的对应。

4.根据权利要求1所述的生物发光断层成像重建方法，其特征在于，所述目标函数的优化求解为：

(1)参数初始化，令t＝0,最大迭代次数为K，初始化S⁰＝(1,...,1)^Τ；

(2)目标函数转化，基于内点法求解：

$S^{t+1}=\underset{0<S<S^{\sup }}{\min }{\left|\right|\mathrm{GS}-J^{+,m}\left|\right|}^2+\mathrm{\&lambda;}\left|S\right|/(\sqrt{\left|S^t\right|}+c);$

将其转化为：

$S^{t+1}=\underset{0<S<S^{\sup }}{\mathrm{mi}}{\left|\right|\mathrm{GS}-J^{+,m}\left|\right|}^2+\mathrm{\&lambda;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_i-\frac{1}{w}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}[\log (u_i+s_i)+\log (u_i-s_i)];$

(3)内部参数设置，设置目标对偶距ζ＞0，w＝1/λ，S＝(1,...,1)^Τ；

(4)预处理共轭梯度法求解牛顿系统方程 $H\left[\begin{array}{c}\mathrm{\&Delta;S}\\ \mathrm{\&Delta;u}\end{array}\right]=-g,$ 其中H为目标函数的Hessian矩阵；

(5)采用回溯线搜索方法计算步长a；

(6)更新迭代解(S,u)＝(S,u)+a(ΔS,Δu)；

(7)计算对偶可行点υ和对偶距ζ；

(8)如果ζ/G(υ)＜1e-6则退出；否则更新w；

(9)得到S^t+1，并令t＝t+1；

(10)当t<K时，转到步骤2；当t＝K时，输出目标解S^t。