CN103300829B - 一种基于迭代重加权的生物自发荧光断层成像方法 - Google Patents

Abstract

本发明涉及一种基于迭代重加权的生物自发荧光断层成像方法及装置。这种方案通过捕捉荧光蛋白基因的肿瘤细胞所发射出的光子信号，能够三维重建出生物体内肿瘤病灶区的大小，同时融合Micro-CT提供的生物体解剖结构信息，还能对病灶区进行定位分析。本发明结合自发荧光成像提供的功能信息与Micro-CT成像提供的结构信息，建立非匀质生物体模型和基于扩散方程的光子传输模型，利用范数正则化与迭代重加权相结合的优化策略，实现生物体内部光源的三维重建。本发明提出的方案能通过更少的观测量重建出更逼近真实解的结果，不但能有效提高求解的计算效率，还能增加重建算法的鲁棒性，适用于实际生物体内部的肿瘤三维探测与定量分析。

Description

一种基于迭代重加权的生物自发荧光断层成像方法

技术领域

本发明属于光学分子影像领域，涉及自发荧光断层成像技术，尤其是一种基于迭代重加权范数优化策略的自发荧光断层成像方法。

背景技术

自发荧光断层成像是光学分子影像的重要模态之一，它改变了传统影像学离体方法不能在体连续观测生物体的局限性，可直接无创地显示生物体内细胞分子水平的生理和病理过程。该技术主要是采用高灵敏度的超低温光学相机捕捉荧光蛋白基因的肿瘤细胞所发射出的光子信号，再结合一定的数学模型与重建算法三维探测出生物体内肿瘤病灶区的大小，同时融合微型计算机断层成像(Micro-CT)提供的生物体解剖结构信息，还能对病灶区进行定位分析。自发荧光断层成像检测灵敏度高，能对肿瘤早期的微小病灶区进行探测，已经发展成为一种理想的活体成像方法，可用于在体观测肿瘤的发生、发展和转移过程，也可应用于抗肿瘤新药研发过程中的药物筛选与疗效评价。

自发荧光信号的重建过程是采用仅从生物体表面获取的边界测量信息去估计生物体内部的肿瘤病灶区三维分布。从数学上说，这一过程中已知的信息非常有限而需求解的量却甚多，于是造成严重的病态性，不但增加了求解的难度，也限制了成像的分辨率和定量的精度。因此，发明一种准确而实用的重建算法对自发荧光断层成像技术来说是至关重要的。目前国内外的研究者都在这方面做了很多工作，例如获取多光谱的边界测量数据，也就是增加更多的独立成分来保证三维断层图像重建的准确性。另外一种常用的方法是仅使用一段混合谱段的测量数据，通过减少未知量的个数，采用较小的可行光源区域得到较准确的重建结果。然而这种方法需要在重建区域中确定一个较小的感兴趣区域，这个区域在实际应用中常常不能很好的确定。

近年来研究人员发现，由于自发荧光断层成像所探测的光子信号往往具有稀疏特征，故而可将这种稀疏性作为先验知识来提高重建图像的准确性。

发明内容

本发明的目的是提供一种基于迭代重加权范数优化策略的自发荧光断层成像方法，融合多角度自发荧光成像提供的光学功能信息与Micro-CT成像提供的生物体解剖结构信息，建立非匀质生物体模型和基于扩散方程的光子传输模型，利用范数正则化与迭代重加权相结合的优化策略，完成生物体内部光源的三维重建，可对相应的肿瘤病灶区进行定量、定位分析。

为了实现上述目的，本发明提出了一种基于迭代重加权的生物自发荧光断层成像方法，其包括如下步骤：

步骤1、将生物体中预定组织器官的Micro-CT成像数据进行合成，建立非匀质三维生物体模型，并对其进行离散化，得到离散化的非匀质三维生物体模型；其中所述非匀质三维生物体模型中各个组织器官具有相应的光学系数值；

步骤2、将利用自发荧光成像得到的所述生物体的多角度二维表面光强分布信息映射到所述离散化的非匀质三维生物体模型的三维表面，得到生物体表面光强分布，并建立光子传输模型；

步骤3、采用有限元方法将上述光子传输模型转换为线性矩阵方程的形式；

步骤4、利用范数正则化与迭代重加权相结合的优化策略，对上述线性矩阵方程进行变换，得到优化目标函数；

步骤5、引入迭代收缩算子，以逐步迭代的方式对所述优化目标函数求解，得到生物体内部的光强分布。

本发明还提出了一种基于迭代重加权的生物自发荧光断层成像装置，其包括：

非匀质三维体模型建立模块，其用于将生物体中预定组织器官的Micro-CT成像数据进行合成，建立非匀质三维生物体模型，并对其进行离散化，得到离散化的非匀质三维生物体模型；其中所述非匀质三维生物体模型中各个组织器官具有相应的光学系数值；

光子传输模型建立模块，其用于将利用自发荧光成像得到的所述生物体的多角度二维表面光强分布信息映射到所述离散化的非匀质三维生物体模型的三维表面，得到生物体表面光强分布，并建立光子传输模型；

光子传输模型转换模块，其用于采用有限元方法将上述光子传输模型转换为线性矩阵方程的形式；

优化目标函数计算模块，其用于利用范数正则化与迭代重加权相结合的优化策略，对上述线性矩阵方程进行变换，得到优化目标函数；

生物体内部光强分布获取模块，其用于引入迭代收缩算子，以逐步迭代的方式对所述优化目标函数求解，得到生物体内部的光强分布。

本发明的有益效果是，通过使用基于迭代重加权范数优化策略的自发荧光断层成像方法，融合多角度自发荧光成像提供的光学功能信息与Micro-CT成像提供的生物体解剖结构信息，采用范数正则化与迭代重加权相结合的优化策略，最终实现生物体内部光源的三维重建，可对相应的肿瘤病灶区进行定量、定位分析。使用该方法的优势为，以信号的稀疏性作为先验知识，能通过更少的观测量重建出更逼近真实解的结果，同时，通过引入范数正则化与迭代重加权策略，既能平衡各系数对最优解的影响，增加了重建算法的鲁棒性，又能有效避免矩阵求逆等复杂运算，提高了求解的计算效率。

附图说明

图1为本发明优选实施例公开的自发荧光断层成像方法流程图。

图2为本发明中非匀质生物体的三维解剖结构图。

图3为本发明中非匀质生物体的三维表面光强图。

图4为使用本发明在不同正则化参数下的重建结果图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

本发明提出了一种基于迭代重加权范数优化策略的生物自发荧光断层成像方法。本发明在自发荧光成像中，底物荧光素在生物体内产生的荧光素酶催化下，与氧气、三磷酸腺苷等物质发生生物化学反应，产生的光子在组织内经历多次散射和吸收后部分穿出体表被探测器接收，从而获取多角度二维生物体表面光强分布信息，能反应出生物体内分子细胞水平的代谢变化。再融合通过Micro-CT成像获取到的生物体解剖结构信息，借助不同组织器官的光学特性参数和分割、配准、网格剖分技术，建立非匀质生物体模型和基于扩散方程的光子传输模型，利用范数正则化与迭代重加权相结合的优化策略，完成生物体内部光源的三维重建，可对相应的肿瘤病灶区进行定量、定位分析。本发明以信号的稀疏性作为先验知识，能通过更少的观测量重建出更逼近真实解的结果，同时，通过引入权函数和l₂范数正则化方法，既能平衡各系数对最优解的影响，增加了重建算法的鲁棒性，又能有效避免矩阵求逆等复杂运算，提高了求解的计算效率。

下面结合附图，详细描述本发明提出的上述断层成像方法。图1示出了本发明优选实施例公开的自发荧光断层成像方法流程图。如图1所示，该方法具体包括如下步骤：

步骤1、利用图像分割技术对Micro-CT成像提供的生物体解剖结构数据进行处理，采用合适的分割方法提取生物体内部主要独立组织器官的Micro-CT成像数据，如利用阈值法分割骨骼等对比度较高的器官，利用区域生长法分割肺部等体积较大的器官，而对于心脏等对比度较低且体积较小的器官则采用基于图谱的分割方法。

步骤2、将提取出的生物体内部主要独立组织器官的Micro-CT成像数据通过设置优先级的方式进行图像合成，形成仅包括主要独立组织器官的三维生物体模型；所述优先级的方式是指在图像合成时，所述主要独立组织器官重叠的区域仅保留优先级高的器官的图像数据，从而消除数据边界可能存在的数据重叠；另外，在图像合成时，所述主要独立组织器官图像的坐标位置根据所述Micro-CT成像数据中相应坐标位置确定。在合成三维生物体模型后，同时为所述主要独立组织器官赋予相应的光学系数值，包括光学吸收系数μ_a和光学扩散系数D，从而完成非匀质三维生物体模型的建立；这些光学系数值反应了不同组织器官作为不同介质，对光的不同吸收和扩散程度，具体系数请参见下表1，从而完成非匀质生物体模型的建立。

表1

图2示出了本发明中将各主要的独立组织器官合成后的三维生物体模型的解剖结构图。如图2所示，该三维生物体模型的主要组织器官包括骨骼、肺部、心脏、肝脏、胃部、肾脏、膀胱，以及皮肤。

步骤3、采用网格剖分技术将步骤2中建立的非匀质三维生物体模型进行离散化处理，并利用配准技术将自发荧光成像提供的多角度二维表面光强分布信息映射到离散化的非匀质生物体模型表面，得到融合后的生物体表面的三维光强分布Φ；同时为离散化处理后的非匀质三维生物体模型建立基于扩散方程的光子传输模型，描述光在此生物体内部的三维光强分布。其中，所述离散化处理是指用网格剖分技术对所述非匀质三维生物体模型进行离散化，即将所述三维生物体模型离散化为多个立体网格；所述自发荧光成像是指对所述生物体表面从多个角度进行自发荧光成像，并获得生物体多个角度的二维表面光强分布信息；然后将所获得的多个角度的二维表面光强分布信息映射到所述离散化的非匀质三维生物体模型的三维表面，得到融合后的生物体表面的三维光强分布Φ。其中，基于扩散方程的光子传输模型如下表示：

$\begin{array}{cc}-\▿\·[D\left(r\right)\▿\mathrm{\Φ}\left(r\right)]+{\mathrm{\μ}}_a\left(r\right)\mathrm{\Φ}\left(r\right)=S\left(r\right)& r\∈\mathrm{\Ω}\\ \mathrm{\Φ}\left(r\right)+2\mathrm{\κ}\left(n\right)D\left(r\right)[v\left(r\right)\·\▿\mathrm{\Φ}\left(r\right)]=0& r\∈\∂\mathrm{\Ω}\end{array}$

式中，是微分运算符，是偏微分运算符，Ω表示离散化处理后的非匀质三维生物体模型的整个区域，表示该三维生物体模型的边界区域，即生物体模型的表面，r是离散化的三维生物体模型上相应离散点的位置向量，包括X、Y、Z坐标；S表示待求解的生物体内部的光强分布，K是预设的生物体与外部媒介的不匹配因子函数，n是预设的生物组织的折射系数，v是的单位法向量；μ_a表示光学吸收系数；D表示光学扩散系数，Φ表示三维光强分布。

图3示出了本发明中生物体表面的三维光强分布图。如图3所示，该生物体的表面上呈现出一块明显的光斑，它是生物体内部光源经传输后在表面所形成的。

步骤4、为了便于后续步骤中的优化处理，采用有限元方法将步骤3中的光子传输模型转换为线性矩阵方程MS＝Φ的形式，其中M为描述非匀质三维生物体模型光学特性的系统矩阵；再利用范数正则化策略，对上述线性矩阵方程进行变换，得到优化目标函数：

$f\left(S\right)=\frac{1}{2}{\left|\right|\mathrm{MS}-\mathrm{\Φ}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|S\left|\right|}_p^p$

式中，λ为正则化参数，表示l₂范数，表示l_p范数，p是介于0与2间的任意数，M中的任一元素如下表示：

其中，是所述有限元方法的点基函数，Ω_x表示离散化处理后的非匀质三维生物体模型的整个区域上的第x个元素，表示该三维生物体模型的边界区域上的第x个元素。

步骤5、为了将优化目标函数中复杂的l_p范数计算转化为仅需要计算l₂范数的形式，引入S的权函数W(S)，它是一个主对角线上的元素值为|S[k]|^2-p的对角矩阵，表达式如下，其中k是大于等于0的整数：

步骤6、结合步骤5中引入的权函数W(S)，可将步骤4中的优化目标函数简化为如下仅含l₂范数的形式：

$f\left(S\right)=\frac{1}{2}{\left|\right|\mathrm{MS}-\mathrm{\Φ}\left|\right|}_2^2+\mathrm{\λ}{S}^T{W}^{-1}\left(S\right)S$

步骤7、为了提高步骤6中优化目标函数的计算效率，引入一个迭代收缩算子P，可避免既费时又占存储空间的矩阵W(S)求逆过程，P的表达式如下：

$P={(\frac{2\mathrm{\λ}}{\mathrm{\α}}I+W\left(S_j\right))}^{-1}W\left(S_j\right)$

式中，I为单位矩阵，α是不小于1的常数，j为迭代次数。

步骤8、结合步骤7中引入的迭代收缩算子P，再利用迭代的方式求解优化目标函数，就可以得到S_j+1，该步骤具体包括：

步骤81，设置初始化生物体内部的光强分布S和迭代次数j；其中，S_j＝S₀，j＝0，S₀为预设的生物体内部的光强分布初始值，可为全零值或其他值；

步骤82，计算权函数W(S_j)：

步骤83，计算迭代收缩算子P，并利用迭代的方式求解优化目标函数，即建立迭代方程，逐步求解生物体内部光强分布S_j+1：

$S_{j+1}=P(\frac{1}{\mathrm{\α}}M(\mathrm{\Φ}-M{S}_j)+S_j);$

步骤84：计算比值判断条件R＜0.1，若为假，则执行j＝j+1，并跳至步骤e2，进入下一次迭代；若为真，则重建停止，得到生物体内部光强分布S的最优解，即最后一次迭代中得到的S_j+1。

在本发明中，融合多角度自发荧光成像提供的光学功能信息与Micro-CT成像提供的生物体解剖结构信息，采用范数正则化与迭代重加权相结合的优化策略，最终实现生物体内部光源的三维重建，可对相应的肿瘤病灶区进行定量、定位分析。与传统方法相比，本发明以信号的稀疏性作为先验知识，能通过更少的观测量重建出更逼近真实解的结果，同时，通过引入权函数和l₂范数正则化方法，既能平衡各系数对最优解的影响，增加了重建算法的鲁棒性，又能有效避免矩阵求逆等复杂运算，提高了求解的计算效率。

图4为使用本发明在不同正则化参数下的重建结果图。图4(a)为重建得到的生物体内部光源在三维非匀质模型中的显示效果，可直观地看出该光源位于生物体的肝脏部位；图4(b)-(g)为重建所得光源在同一横截面上的二维显示效果，其中黑色圆圈表示真实光源的位置；图4(c)-(g)代表当正则化参数分别为10^-3，10^-5，10^-7，10^-9，10^-11时基于非匀质生物体模型的重建结果，可见重建所得光源的中心与真实光源的中心基本吻合，且随着正则化参数的改变，重建结果依然准确而稳定；图4(b)则代表当正则化参数为10^-7时基于匀质生物体模型的重建结果，即不考虑生物体内部组织器官对光子传输的影响，将该结果与4(e)比较可发现，即使是在同一正则化参数下，基于匀质模型的重建结果明显不如基于非匀质模型的重建结果准确。

表2为使用本发明与另两种传统方法的重建速度比较，包括经典牛顿法和共轭梯度法。为使重建结果的比较更公平，不同方法均针对相同大小的生物体网格数据量，该网格为3397个节点和18344个四面体单元。从这三种方法的重建结果可见，基于本发明的重建时间9.28秒明显短于另两种传统方法，这是因为本发明所引入的优化策略可有效避免费时的矩阵求逆等过程，提高了重建的效率。

表2

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于迭代重加权的生物自发荧光断层成像方法，其包括如下步骤：

步骤1、将生物体中预定组织器官的Micro-CT成像数据进行合成，建立非匀质三维生物体模型，并对其进行离散化，得到离散化的非匀质三维生物体模型；其中所述非匀质三维生物体模型中各个组织器官具有相应的光学系数值；

步骤2、将利用自发荧光成像得到的所述生物体的多角度二维表面光强分布信息映射到所述离散化的非匀质三维生物体模型的三维表面，得到生物体表面光强分布，并建立光子传输模型；

步骤3、采用有限元方法将上述光子传输模型转换为线性矩阵方程的形式；

步骤4、利用范数正则化与迭代重加权相结合的优化策略，对上述线性矩阵方程进行变换，得到优化目标函数；

步骤5、引入迭代收缩算子，以逐步迭代的方式对所述优化目标函数求解，得到生物体内部的光强分布。

2.如权利要求1所述的自发荧光断层成像方法，其特征在于，步骤1具体包括：

步骤11、利用图像分割技术将生物体的Micro-CT成像数据进行分割，提取出生物体内部预定组织器官的Micro-CT成像数据；

步骤12、将所提取的预定组织器官的Micro-CT成像数据进行图像合成，建立仅包括所述预定组织器官的三维生物体模型；

步骤13、对于所建立的三维生物体模型中的各个组织器官赋予相应的光学系数值，以形成非匀质的三维生物体模型；

步骤14、对所述非匀质的三维生物体模型进行离散化处理，得到离散化的非匀质三维生物体模型。

3.如权利要求1所述的自发荧光断层成像方法，其特征在于，步骤2具体包括：

步骤21、利用自发荧光成像技术从多角度获得所述生物体表面的二维表面光强分布信息；

步骤22、利用配准技术将从多角度获得的生物体的二维表面光强分布信息映射至所述离散化的非匀质三维生物体模型表面，得到生物体表面的三维光强分布信息；

步骤23、建立所述离散化的非匀质三维生物体模型的所述三维光强分布信息、各个组织器官的光学系数值与生物体内部光强分布信息的光子传输模型。

4.如权利要求3所述的自发荧光断层成像方法，其特征在于，所述光子传输模型如下表示：

$\left\{\begin{array}{cc}-\▿\·[D\left(r\right)\▿\mathrm{\Φ}\left(r\right)]+{\mathrm{\μ}}_a\left(r\right)\mathrm{\Φ}\left(r\right)=S\left(r\right)& r\∈\mathrm{\Ω}\\ \mathrm{\Φ}\left(r\right)+2\mathrm{\κ}\left(n\right)D\left(r\right)[v\left(r\right)\·\▿\mathrm{\Φ}\left(r\right)]=0& r\∈\∂\mathrm{\Ω}\end{array}\right.$

式中，是微分运算符，是偏微分运算符，Ω表示离散化的非匀质三维生物体模型的整个区域，表示该三维生物体模型的边界区域，r是离散化的非匀质三维生物体模型上相应离散点的位置向量；S表示生物体内部的光强分布，κ是预设的生物体与外部媒介的不匹配因子函数，n是预设的生物组织的折射系数，v是的单位法向量；μ_a和D分别为各个组织器官上的光学系数值，分别表示光学吸收系数和光学扩散系数，Φ表示生物体表面的三维光强分布信息。

5.如权利要求4所述的自发荧光断层成像方法，其特征在于，步骤3中所述线性矩阵方程如下表示：

MS＝Φ

其中，M为通过有限元方法得到的描述非匀质三维生物体模型的光学特性矩阵。

6.如权利要求5所述的自发荧光断层成像方法，其特征在于，步骤4中利用范数正则化优化策略所得到的目标函数如下表示：

$f\left(S\right)=\frac{1}{2}{\left|\right|\mathrm{MS}-\mathrm{\Φ}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|S\left|\right|}_p^p$

其中，λ为正则化参数，表示l₂范数，表示l_p范数，p是介于0与2间的任意数。

7.如权利要求6所述的自发荧光断层成像方法，其特征在于，步骤4中将利用范数正则化优化策略所得到的所述目标函数，再利用迭代重加权优化策略进一步优化得到如下表示的目标函数：

$f\left(S\right)=\frac{1}{2}{\left|\right|\mathrm{MS}-\mathrm{\Φ}\left|\right|}_2^2+{\mathrm{\λS}}^T{W}^{-1}\left(S\right)S$

其中，W(S)为S的权函数，其是一个主对角线上的元素值为|S[k]|^2-p的对角矩阵，其中k是大于等于0的整数，表达式如下所示：

8.如权利要求7所述的自发荧光断层成像方法，其特征在于，步骤5中所述的迭代收缩算子表达式为：

$P={(\frac{2\mathrm{\λ}}{\mathrm{\α}}I+W\left(S_j\right))}^{-1}W\left(S_j\right)$

其中，I为单位矩阵，α是不小于1的常数，j为迭代次数。

9.如权利要求8所述的自发荧光断层成像方法，其特征在于，所述步骤5具体包括以下步骤：

步骤51、设置初始化生物体内部的光强分布S和迭代次数j；其中，S_j＝S₀，j＝0，S₀为预设的生物体内部的光强分布初始值；

步骤52、计算权函数W(S_j)：

步骤53、计算迭代收缩算子P，并利用迭代的方式求解优化目标函数，逐步求解生物体内部的光强分布S_j+1：

$S_{j+1}=P(\frac{1}{\mathrm{\α}}M(\mathrm{\Φ}-{\mathrm{MS}}_j)+S_j);$

其中，α是不小于1的常数；

步骤54、计算比值判断条件R＜0.1是否为真，若为假，则执行j＝j+1，并跳至步骤52，进入下一次迭代；若为真，则重建停止，得到生物体内部光强分布S的最优解，即最后一次迭代中得到的S_j+1。

10.一种基于迭代重加权的生物自发荧光断层成像装置，其包括：

非匀质三维体模型建立模块，其用于将生物体中预定组织器官的Micro-CT成像数据进行合成，建立非匀质三维生物体模型，并对其进行离散化，得到离散化的非匀质三维生物体模型；其中所述非匀质三维生物体模型中各个组织器官具有相应的光学系数值；

光子传输模型建立模块，其用于将利用自发荧光成像得到的所述生物体的多角度二维表面光强分布信息映射到所述离散化的非匀质三维生物体模型的三维表面，得到生物体表面光强分布，并建立光子传输模型；

光子传输模型转换模块，其用于采用有限元方法将上述光子传输模型转换为线性矩阵方程的形式；

优化目标函数计算模块，其用于利用范数正则化与迭代重加权相结合的优化策略，对上述线性矩阵方程进行变换，得到优化目标函数；

生物体内部光强分布获取模块，其用于引入迭代收缩算子，以逐步迭代的方式对所述优化目标函数求解，得到生物体内部的光强分布。