CN102393969A - 基于生物组织特异性的光学三维成像方法 - Google Patents

Abstract

本发明公开了一种基于生物组织特异性的光学三维成像方法，解决了现有技术中无法实现对具有不规则解剖结构和多种散射特性组织的复杂生物体进行准确快速的光学三维成像的问题。该方法基于生物组织特异性光传输混合数学模型和完全稀疏正则化方法建立目标函数，采用基于任务导向的混合优化方法进行求解，以实现体内靶向目标的光学三维成像。具体实现步骤包括：采集数据、预处理、获取解剖结构、获取表面光学数据、重建光学特性参数、建立光传输模型、建立系统方程、建立目标函数、求解目标函数和显示结果。本发明具有能够实现复杂生物体中光源的准确、快速、高分辨重建的优点，可以用于在体光学三维成像领域。

Description

基于生物组织特异性的光学三维成像方法

技术领域

本发明属于光学技术领域，更进一步涉及光学分子成像技术领域中的基于生物组织特异性的光学三维成像方法。该方法可应用于基因表达、肿瘤检测、药物研发与疗效评价等小动物实验和预临床实验中获取靶向目标的位置和强度信息。

背景技术

作为一种新兴的光学成像技术，光学三维成像是通过融合生物体体表测量的多角度光学信号、生物体的解剖结构和组织光学参数信息，基于精确的生物组织中的光传输模型重建活体生物体体内靶向目标的位置和强度分布信息，其中，生物组织中光传输过程的精确描述和靶向目标的准确快速重建是光学三维成像方法实现的基础。在光学三维成像技术中，激发荧光断层成像、扩散光学层析断层成像和自发荧光断层成像是三种主要的成像模态。

基于自发荧光断层成像模态，北京工业大学在其专利申请文件“基于单视图的多光谱自发荧光断层成像重建方法”(申请号200810116818.4，申请日2008.7.18，授权号ZL200810116818.4，授权日2010.6.2)中提出了一种基于单幅视图的多光谱自发荧光断层成像重建方法。该专利技术基于扩散近似方程，考虑生物体的非匀质特性和自发荧光光源的光谱特点，利用在单个角度测量的多个谱段荧光数据，重建生物体体内靶向目标的位置和强度分布信息。但是，由于扩散近似方程只适用于描述高散射特性组织中的光传输过程，对于低散射特性和空腔组织，它的求解精度很低。因此，该专利技术对于具有多种散射特性组织的生物体求解精度差，很难准确地获得生物体体内靶向目标的位置和强度分布信息。

为了解决基于扩散近似方程的光学三维成像方法的局限性，研究人员还提出了基于辐射传输方程的高阶近似方程的光学三维成像方法，参见Alexander D.Klose，“The inverse source problem based on the radiative transfer equation inoptical molecular imaging，”Journal of Computational Physics 202，323-345(2005)；Yujie Lu，“Spectrally resolved bioluminescence tomography with thethird-order simplified spherical harmonics approximation，”Physics inMedicine and Biology 54，6477-6493(2009)；Wengxiang Cong，“Bioluminescencetomography based on the phase approximation model，”Journal of OpticsSociety of America A 27，174-179(2010)。该类方法可以对同时存在高散射特性和低散射特性组织的生物体进行准确成像。但是，由于辐射传输方程的高阶近似方程求解难度非常大，对于具有不规则解剖结构的复杂生物体，这些方法的计算复杂度将远远超过实际应用的承受能力。此外，这些方法也不能对存在空腔组织的情况进行准确的三维成像，如动物体的胃、膀胱等。

为了对存在空腔组织的情况进行准确的三维成像，研究人员又提出了一种基于混合光传输方程的光学三维成像方法，参见Hamid Dehghani，“Optical tomographyin the presence of void regions，”Journal of Optics Society of America A17，1659-1670(2000)。该方法基于扩散近似方程和朗伯源特性方程，可以对同时存在高散射特性和空腔组织的生物体进行准确成像，但是由于扩散近似方程的局限性，这种方法对于具有多种散射特性组织的生物体仍然不能进行准确的成像。

此外，中国科学院自动化研究所在其专利申请文件“激发荧光断层成像的快速稀疏重建方法和设备”(公开号CN102034266A，申请号201010573795.7，申请日2010.11.30)中提出了一种基于稀疏正则化策略的快速激发荧光断层成像重建方法。该方法考虑了靶向目标的稀疏分布特性，建立了基于l₁范数的稀疏正则化目标函数，能够很好地改善光学三维成像的准确性和分辨率。但是该方法仅采用了单一的优化策略对目标函数进行求解，无法实现不同问题规模下，成像准确度和成像速度的最优平衡。此外，该方法还没有考虑测量数据相对于生物体区域的稀疏特性。

综上所述，对于具有不规则解剖结构和多种散射特性组织的复杂生物体，已有技术中的基于单一近似方程或混合光传输方程的光学三维成像方法和基于单一优化策略的光学三维成像方法，都无法实现对生物体体内靶向目标的准确、快速、高分辨的重建。

发明内容

本发明的目的在于克服上述已有光学三维成像技术存在的不足，提出一种基于生物组织特异性的光学三维成像方法，以实现对具有不规则解剖结构和多种散射特性组织的复杂生物体体内靶向目标的准确、快速、高分辨率重建。

实现本发明方法的主要思路是：根据生物体在解剖结构和组织光学特性参数方面的特征差异，利用外推边界和折射率不匹配边界条件，建立生物组织特异性光传输混合数学模型。在该模型基础上，利用多级自适应有限元方法，结合生物体体内靶向目标的稀疏分布特性和测量数据相对于生物体区域的稀疏特性，建立完全稀疏正则化目标函数。采用基于任务导向的混合优化方法求解目标函数，实现生物体体内靶向目标的准确、快速、高分辨率重建。

根据上述主要思路，本发明方法的具体实现包括如下步骤：

(1)采集数据

利用多模态分子影像系统，依次采集用于光学三维成像的多角度荧光数据、用于光学特性参数重建的多角度激光数据和用于获取生物体解剖结构的计算机断层成像投影数据。

(2)预处理

2a)利用多模态分子影像系统中的预处理软件对采集的多角度激光数据和荧光数据进行去除背景噪声、提取感兴趣区域、补偿坏点预处理；

2b)利用多模态分子影像系统中的预处理软件对计算机断层成像投影数据进行补偿坏点坏线、亮场校正、几何校正预处理。

(3)获取解剖结构

利用3DMED软件对预处理后的计算机断层成像投影数据进行三维重建，获得生物体三维体数据；利用3DMED软件中的人机交互式分割方法对获得的生物体三维体数据进行器官分割，获取生物体解剖结构。

(4)获取表面光学数据

对步骤(3)获取的生物体解剖结构和步骤(2)获取的荧光和激光数据，应用非接触式光学断层成像方法中的生物体表面三维能量重建技术获取生物体表面的三维荧光和激光数据分布。

(5)重建光学特性参数

对步骤(3)获取的生物体解剖结构和步骤(4)获取的生物体表面三维激光数据分布，应用基于区域的扩散光学层析成像方法重建生物体各个组织的光学特性参数。

(6)建立光传输模型

6a)根据生物体在解剖结构和组织光学特性参数方面的差异，将生物组织划分为高散射特性组织、低散射特性组织、空腔组织三大类；

6b)根据扩散近似方程、简化球谐波近似方程、朗伯源特性方程的适用范围，采用相应的方程描述不同特性生物组织中的光传输过程；

6c)构造外推边界和折射率不匹配边界条件耦合不同特性生物组织的光传输方程，建立统一的、描述光在整个生物体中传输过程的混合数学模型。

(7)建立系统方程

7a)如果是在第一级离散网格上建立系统方程，利用Amira软件对生物体的高散射特性和低散射特性组织进行初始离散；否则，利用编写的程序对上一级离散网格进行手动调整；

7b)在第k级离散网格上，利用多级自适应有限元方法对步骤6c)建立的混合数学模型进行离散化，组装每个离散点上的子系统方程建立总的系统方程：

A_kS_k＝Φ_k

其中，A_k是第k级离散网格上的系统矩阵，依赖于生物体内三类特性生物组织的分布和生物组织的光学特性参数；

S_k是第k级离散网格上的靶向目标能量密度分布；

Φ_k是第k级离散网格边界节点上的光通量密度。

(8)建立目标函数

根据第k级离散网格边界节点上的光通量密度计算值和测量值之间的误差，结合靶向目标的稀疏分布约束，建立下列目标函数：

$\min \underset{S_k^{\inf }\≤S_k\≤S_k^{\sup }}{\mathrm{\Θ}\left(S_k\right)}=\{\frac{1}{2}{\left|\right|A_k{S}_k-{\mathrm{\Φ}}_k^m\left|\right|}_1+{\mathrm{\λ}}_k{\left|\right|S_k\left|\right|}_1\}$

其中，Θ(S_k)是第k级离散网格上的目标函数；

是第k级离散网格上的靶向目标能量密度的下限；

是第k级离散网格上的靶向目标能量密度的上限；

是第k级离散网格边界节点上的光通量密度的测量值；

‖F‖₁定义为求解矩阵F的l₁范数；

λ_k是第k级离散网格上的正则化因子。

(9)求解目标函数

9a)采用基于任务导向的混合优化方法求解建立的目标函数，根据第k级离散网格上形成的系统矩阵的大小选择合适的优化方法，获得第k级离散网格上的靶向目标能量密度分布；

9b)利用靶向目标能量密度分布计算第k级离散网格边界节点上的光通量密度；

9c)判断第k级离散网格边界节点上的光通量密度的测量值和计算值之差，如果小于给定阈值，则目标函数求解过程结束，获得光学三维成像的靶向目标重建结果，转向步骤(10)；否则，转向步骤9d)；

9d)根据步骤9a)获得的靶向目标能量密度分布和步骤9b)获得的边界节点上的光通量密度计算值，调整第k+1级离散网格，转向步骤(7)。

(10)显示结果。

本发明与现有技术相比具有如下优点：

第一，本发明由于同时考虑生物体在解剖结构和组织光学特性参数方面的差异建立光传输混合数学模型，克服了现有技术中基于单一近似方程或混合光传输方程的光学三维成像方法的在重建精度和效率方面的局限性，能够对具有不规则解剖结构和多种散射特性组织的复杂生物体的靶向目标进行准确、快速成像。

第二，本发明由于同时考虑靶向目标的稀疏分布特性和测量数据相对于生物体区域的稀疏特性建立目标函数，克服了现有技术中仅考虑靶向目标的稀疏分布特性带来的成像分辨率不足的问题，能够有效地改善光学三维成像的分辨率。

第三，本发明由于采用基于任务导向的混合优化方法对目标函数进行求解，克服了现有技术中采用单一优化策略求解目标函数过程中存在的成像准确度和成像效率不能兼顾的问题，能够实现复杂生物体内靶向目标的准确、快速重建。

附图说明

图1为本发明的流程图；

图2为本发明实施例中的生物组织特异性光传输混合数学模型的建立过程示意图；

图3为本发明实施例中耦合不同特性生物组织中光传输方程的一种方式示意图；

图4为本发明实施例中耦合不同特性生物组织中光传输方程的另一种方式示意图；

图5为本发明实施例中耦合不同特性生物组织中光传输方程的另一种方式示意图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照图1，本发明的具体步骤如下：

步骤1，采集数据

利用多模态分子影像系统，依次采集用于光学三维成像的多角度荧光数据、用于光学特性参数重建的多角度激光数据和用于获取动物体解剖结构的计算机断层成像投影数据。

多角度荧光数据的采集，控制成像体等间隔旋转一定角度，一般不大于90°(本例中选90°)，利用多模态分子影像系统中的CCD相机采集不少于四个角度的荧光数据(本例中为四个角度)。继续旋转，使成像体转回到荧光数据初始采集的位置。

多角度激光数据的采集，打开多模态分子影像系统中的激光器光源，在对侧利用CCD相机采集成像体表面透射的激光数据。控制成像体等间隔旋转一定角度，一般不大于90°(本例中选90°)，采集不少于四个角度的激光数据(本例中为四个角度)。最终，关闭激光器光源，旋转成像体使之转回到激光数据初始采集的位置。

多角度计算机断层成像投影数据的采集，控制成像体等间隔旋转一定小角度，一般不大于1°(本例中选1°)，采集不少于360个角度的投影数据(本例中为360个角度)。

步骤2，预处理

2a)利用多模态分子影像系统中的预处理软件对采集的多角度激光数据和荧光数据进行去除背景噪声、提取感兴趣区域、补偿坏点预处理；

2b)利用多模态分子影像系统中的预处理软件对计算机断层成像投影数据进行补偿坏点坏线、亮场校正、几何校正预处理。

步骤3，获取解剖结构

采用3DMED软件对预处理后的计算机断层成像投影数据进行三维重建，获得动物体三维体数据；应用3DMED软件中的人机交互式分割方法对获得的动物体三维体数据进行器官分割，获取动物体解剖结构。

步骤4，获取表面光学数据

利用步骤3获取的动物体解剖结构和步骤2获取的荧光和激光数据，应用非接触式光学断层成像方法(申请号200910024292.1，申请日2009.10.13，授权号ZL200910024292.1，授权日2011.4.6)中描述的生物体表面三维能量重建技术获取动物体表面的三维荧光和激光数据分布。

步骤5，重建光学特性参数

利用步骤3获取的动物体解剖结构和步骤4获取的动物体表面三维激光数据分布，应用基于区域的扩散光学层析成像算法重建动物体内各个组织的光学特性参数。

基于区域的扩散光学层析成像的实现包含如下两层内容。首先，对微型计算机断层成像系统获取的动物体解剖结构进行器官的区域抽取和数字化表示，采用基于区域标记的体元网格表示不同的器官，并假设器官内的组织光学特性参数上均匀的，将大规模的体元光学特性参数重建变为多器官内部均匀光学特性参数的重建。第二，在上一步获得的目标动物体基于区域标记的体元网格的基础上，通过对区域内同类体元重建的光学特性参数的合并，实现不同区域组织光学特性参数的高精度同时重建。由于面向动物体全身测量，复杂的光学特性参数分布特征突破了基于扩散近似方程的生物组织中光传输模型的有效范围，因此，基于区域的扩散光学层析成像算法是采用辐射传输方程的简化球谐波近似方程描述生物组织中的光传输过程。

步骤6，建立光传输模型

下面结合图2详细说明生物组织特异性光传输混合数学模型的建立过程。

6a)动物体内生物组织的特异性分类

根据步骤3获取的动物体解剖结构和步骤5获取的各个生物组织的光学特性参数，对动物体内的生物组织进行特异性分类。生物组织的特异性分类包含两层含义，即根据解剖结构差异分类和根据组织光学特性参数差异分类。首先，根据动物体在解剖结构方面的特征差异，划分不同器官的生物组织区域。图2(a)是基于动物体解剖结构差异的生物组织分类示意图，其中不同颜色的区域代表不同类别的生物组织，在本实施例中包括骨骼、心脏、肾脏、肝脏、肺、膀胱、胃和脂肪八种组织。

其次，根据不同生物组织在光学特性参数方面的特征差异，应用下式和下面的准则将生物组织划分为具有不同散射特性的区域：

$\▿=\frac{{\mathrm{\μ}}_s^{\′}}{{\mathrm{\μ}}_a}$

式中，▽是不同散射特性组织的划分标准因子，

是生物组织的约化散射系数，μ_a是生物组织的吸收系数。划分的准则：如果▽≥10，那么生物组织划分为高散射特性组织；如果0＜▽＜10，那么生物组织划分为低散射特性组织；如果▽＝0，那么生物组织划分为空腔组织。图2(b)是基于组织光学特性参数差异的生物组织分类示意图，其中三个矩形框分别代表具有不同散射特性的组织类。在本实施例中，选用中心波长在650nm的组织光学参数来说明基于组织光学特性参数差异的生物组织分类过程。表1列出了上述几种常见生物组织的组织光学特性参数，并计算了相应的划分标准因子。由于胃和膀胱的内部是一个囊腔状区域，在里面只发生光的吸收，几乎不发生光的散射，因此这两种生物组织划分为空腔组织，如图2(b)中的底部的矩形框中所示。骨骼、心脏、肾脏和脂肪，由于它们的划分标准因子▽≥10，因此这四种生物组织划分为高散射特性组织，如图2(b)的顶部的矩形框中所示。肺和肝脏，由于它们的划分标准因子▽＜10，因此这两种生物组织划分为低散射特性组织，如图2(b)的中间层的矩形框中所示。

表1不同生物组织的组织光学特性参数(650nm)

6b)不同特性生物组织的特异性表述

光在生物组织中的传输过程可以用辐射传输方程精确地描述。辐射传输方程是一个复杂的积分-微分方程，在复杂的生物组织中很难直接求解，并且求解的时间代价很高。因此，通常采用辐射传输方程的近似方程描述光在生物组织中的传输过程，例如扩散近似方程、球谐波近似方程、简化球谐波近似方程、离散坐标近似方程和相位近似方程。这些近似方程都有其各自的优缺点和适用范围：扩散近似方程求解速度快，计算复杂度低，但只能准确描述高散射特性生物组织中的光传输过程；球谐波、简化球谐波、离散坐标和相位近似方程等辐射传输方程的高阶近似形式，在较高阶数情况下，能够准确地描述光在任意散射特性生物组织中的传输过程，但是求解难度非常大；尤其是球谐波近似、离散坐标近似和相位近似方程，计算复杂度将远远超过实际应用的承受能力。研究表明，五阶简化球谐波近似能够有效的求解，并能够达到实际应用所能承受的求解精度。然而，所有的这些辐射传输方程的近似方程，都无法准确描述光在空腔组织中的传输过程。空腔是一种散射特性为零，只存在吸收的组织。光在空腔组织中沿直线传输，可以用辐射度学理论或辐射传输方程准确描述。考虑辐射传输方程求解的时间代价，一般采用基于辐射度学理论的漫射光的朗伯源特性方程描述光在空腔组织中的传输过程。

根据上述辐射传输方程的各种近似方程和漫射光的朗伯源特性方程的适用范围和优劣势，采用不同方程描述不同特性生物组织中的光传输过程，完成不同特性生物组织的特异性表述：对于高散射特性组织，采用求解速度快的扩散近似方程；对于低散射特性组织，采用简化球谐波近似方程；对于空腔组织，采用漫射光的朗伯源特性方程。图2(c)是采用不同光传输方程描述不同特性生物组织中光传输过程的分类示意图，图2(d)是建立的不同特性生物组织的特异性表述结果示意图。

6c)光传输混合数学模型的建立

步骤6b)已经建立了不同特性生物组织的特异性表述，但是描述不同特性生物组织的光传输方程之间还是相互独立的，需要构造合适的边界条件对不同光传输方程进行耦合。根据不同散射特性生物组织的区域划分，同时结合不同光传输方程的特点，将不同光传输方程或不同类区域之间的耦合划分为三类情况，包括空腔组织向散射特性组织的耦合、散射特性组织向空腔组织的耦合以及散射特性组织之间的耦合。图3描述了空腔组织向散射特性组织的耦合过程，图4是散射特性组织向空腔组织耦合的过程示意图，图5是散射特性组织之间耦合的过程示意图。下面结合图3、图4和图5详细说明三类情况下的耦合过程。

第一类情况，空腔组织向散射特性组织的耦合。研究表明，入射到散射介质的准直笔形光束可以等效为一个各项同性光源，位于距离散射介质表面一个平均自由程的位置。也就是说，光子传输的一个平均自由程是任意光子转化为完全漫射光的长度尺度。漫射光经过空腔组织作用后，会转变为非漫射光。如前所述，这种非漫射光进入散射特性组织后，需要继续传输一个光子传输平均自由程的距离才能转化为一个各项同性的光源。因此，在处理空腔组织向散射特性组织耦合的情况下，需要将空腔组织与散射特性组织之间的边界向散射特性组织区域外推一个光子传输平均自由程的距离，形成外推边界，并将非漫射光转化为位于外推边界上的各项同性光源。图3描述了空腔组织向散射特性组织耦合情况下的处理过程示意图。在外推边界上形成的各项同性光源可以通过如下公式确定：

$S_0\left(r_s\right)={\∫}_S{\∫}_S\mathrm{dSd}{S}^{\′}\frac{\mathrm{\ϵ}\left(r^{\′}\right)\mathrm{\φ}\left(r^{\′}\right)}{\mathrm{\π}}G(r,r^{\′})\frac{1}{{|r_s-r|}^2}\mathrm{\δ}\left(\right|r_s-r|-\frac{1}{{\mathrm{\μ}}_t^{\′}})\mathrm{\δ}({\mathrm{\υ}}_{r-r^{\′}}-\frac{r_s-r}{|r_s-r|})$

式中，S₀(r_s)是在外推边界的r_s点形成的各项同性光源；S是空腔组织和散射特性组织之间的边界；ε(r′)是空腔组织和散射特性组织之间的折射率不匹配因子，可以通过式ε(r′)＝1/2A_n(r′)计算，其中A_n(r′)＝-1.4399n^-2+0.7099n^-1+0.6681+0.0636n；φ(r′)是有散射特性组织入射到空腔组织的光流率；G(r，r′)是漫射光在空腔组织中传输的传输函数；

是一个光子传输平均自由程；定义为从点r′到点r的单位方向向量。

第二类情况，散射特性组织向空腔组织的耦合。任意光经在散射特性组织中的充分传输后，会转化成完全漫射光，这种完全漫射光会在散射特性组织和空腔组织的边界形成指向空腔组织的光通量，而这个光通量就是形成空腔组织中光传输的面光源。考虑散射特性组织和空腔组织之间的折射率不匹配条件，经散射特性组织充分散射后形成的光流率可以通过下式转化为指向空腔组织的光通量：

J⁺(r′)＝ε(r′)φ(r′)

式中，J⁺(r′)是在散射特性组织与空腔组织边界上的r′点处的光通量，其方向指向空腔组织。图4描述了散射特性组织向空腔组织耦合情况下的处理过程示意图。

第三类情况，散射特性组织之间的耦合。图5描述了散射特性组织之间耦合情况下的处理过程示意图。在这种情况下，需要将在散射特性组织边界上形成的光流率转化为小体光源，具体实现过程如下所述。首先，考虑不同散射特性组织之间的折射率不匹配条件，将在边界上形成的光流率转化为光通量：

J_n(r′)＝ε(r′)φ(r′)

式中，J_n(r′)是不同散射特性组织边界上的r′点处的光通量，其方向指向出射组织。

其次，应用下述公式将形成的光通量转化为小体光源：

$q_0\left(r^{\′}\right)=\frac{J_n\left(r^{\′}\right)\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i}{\underset{j=1}{\overset{Q}{\mathrm{\Σ}}}{V}_j}$

式中，q₀(r′)是在不同散射特性组织边界点r′处形成的小体光源；M是与点r′相连接的所有面片的数量，S_i是其中第i个面片的面积；Q是与点r′相连接的所有四面体的体积，V_j是其中第j个四面体的体积。

将不同特性生物组织中的光传输方程与上面建立的边界耦合条件进行联立，同时采用数值方法对建立的联立方程进行数值离散，建立描述光在多种散射特性生物组织中传输的混合数学方程。

步骤7，建立系统方程

7a)如果是在第一级离散网格上建立系统方程，利用Amira软件对生物体的高散射特性和低散射特性组织进行初始离散；否则，利用编写的程序对上一级离散网格进行手动调整；

7b)系统方程的建立，在第k级离散网格上，利用多级自适应有限元方法对步骤6c)建立的混合数学方程进行离散化，通过组装每个离散点上的子系统方程建立总的系统方程：

A_kS_k＝Φ_k

其中，A_k是第k级离散网格上的系统矩阵，依赖于动物体内三类特性生物组织的分布和生物组织的光学特性参数；

S_k是第k级离散网格上的靶向目标能量密度分布；

Φ_k是第k级离散网格上的边界节点上的光通量密度。

步骤8，建立目标函数

根据第k级离散网格边界节点上的光通量密度计算值和测量值之间的误差，结合靶向目标的稀疏分布约束，建立下列目标函数：

$\min \underset{S_k^{\inf }\≤S_k\≤S_k^{\sup }}{\mathrm{\Θ}\left(S_k\right)}=\{\frac{1}{2}{\left|\right|A_k{S}_k-{\mathrm{\Φ}}_k^m\left|\right|}_1+{\mathrm{\λ}}_k{\left|\right|S_k\left|\right|}_1\}$

其中，Θ(S_k)是建立的目标函数；

是第k级离散网格上的靶向目标能量密度的下限；

是第k级离散网格上的靶向目标能量密度的上限；

是第k级离散网格边界节点上的光通量密度的测量值；

‖F‖₁定义为求解矩阵F的l₁范数；

λ_k是第k级离散网格上的正则化参数。

步骤9，求解目标函数

9a)采用基于任务导向的混合优化方法求解建立的目标函数，根据第k级离散网格上形成的系统矩阵的大小选择合适的优化方法，获得第k级离散网格上的靶向目标能量密度分布；其中，对于小系统矩阵，采用求解速度较快的改进奇异值分解法进行求解；对于大系统矩阵，采用基于预处理系统矩阵的混合迭代法进行求解；

9b)利用靶向目标能量密度分布计算第k级离散网格边界节点上的光通量密度；

9c)如果第k级离散网格边界节点上的光通量密度的测量值和计算值之差小于给定阈值，则目标函数求解过程结束，获得光学三维成像的靶向目标重建结果，转向步骤10；否则，转向步骤9d)；

9d)根据步骤9a)获得的靶向目标能量密度分布和步骤9b)获得的边界节点上的光通量密度计算值，调整第k+1级离散网格，转向步骤7。

步骤10，显示结果

对步骤9c)获得的靶向目标重建结果和步骤3获取的动物体解剖结构进行图像融合，将重建的靶向目标在动物体中进行三维显示。

Claims (5)

1.基于生物组织特异性的光学三维成像方法，包括：

(1)采集数据

利用多模态分子影像系统，依次采集用于光学三维成像的多角度荧光数据、用于光学特性参数重建的多角度激光数据和用于获取生物体解剖结构的计算机断层成像投影数据；

(2)预处理

2a)利用多模态分子影像系统中的预处理软件对采集的多角度激光数据和荧光数据进行去除背景噪声、提取感兴趣区域、补偿坏点预处理；

2b)利用多模态分子影像系统中的预处理软件对计算机断层成像投影数据进行补偿坏点坏线、亮场校正、几何校正预处理；

(3)获取解剖结构

利用3DMED软件对预处理后的计算机断层成像投影数据进行三维重建，获得生物体三维体数据；利用3DMED软件中的人机交互式分割方法对获得的生物体三维体数据进行器官分割，获取生物体解剖结构；

(4)获取表面光学数据

对步骤(3)获取的生物体解剖结构和步骤(2)获取的荧光和激光数据，应用非接触式光学断层成像方法中的生物体表面三维能量重建技术获取生物体表面的三维荧光和激光数据分布；

(5)重建光学特性参数

对步骤(3)获取的生物体解剖结构和步骤(4)获取的生物体表面三维激光数据分布，应用基于区域的扩散光学层析成像方法重建生物体各个组织的光学特性参数；

(6)建立光传输模型

6a)根据生物体在解剖结构和组织光学特性参数方面的差异，将生物组织划分为高散射特性组织、低散射特性组织、空腔组织三大类；

6b)根据扩散近似方程、简化球谐波近似方程、朗伯源特性方程的适用范围，采用相应的方程描述不同特性生物组织中的光传输过程；

6c)构造外推边界和折射率不匹配边界条件耦合不同特性生物组织的光传输方程，建立统一的、描述光在整个生物体中传输过程的混合数学模型；

(7)建立系统方程

7a)如果是在第一级离散网格上建立系统方程，利用Amira软件对生物体的高散射特性和低散射特性组织进行初始离散；否则，利用编写的程序对上一级离散网格进行手动调整；

7b)在第k级离散网格上，利用多级自适应有限元方法对步骤6c)建立的混合数学模型进行离散化，组装每个离散点上的子系统方程建立总的系统方程：

A_kS_k＝Φ_k

其中，A_k是第k级离散网格上的系统矩阵，依赖于生物体内三类特性生物组织的分布和生物组织的光学特性参数；

S_k是第k级离散网格上的靶向目标能量密度分布；

Φ_k是第k级离散网格边界节点上的光通量密度；

(8)建立目标函数

根据第k级离散网格边界节点上的光通量密度计算值和测量值之间的误差，结合靶向目标的稀疏分布约束，建立下列目标函数：

$\min \underset{S_k^{\inf }\≤S_k\≤S_k^{\sup }}{\mathrm{\Θ}\left(S_k\right)}=\{\frac{1}{2}{\left|\right|A_k{S}_k-{\mathrm{\Φ}}_k^m\left|\right|}_1+{\mathrm{\λ}}_k{\left|\right|S_k\left|\right|}_1\}$

其中，Θ(S_k)是第k级离散网格上的目标函数；

是第k级离散网格上的靶向目标能量密度的下限；

是第k级离散网格上的靶向目标能量密度的上限；

是第k级离散网格边界节点上的光通量密度的测量值；

‖F‖₁定义为求解矩阵F的l₁范数；

λ_k是第k级离散网格上的正则化因子；

(9)求解目标函数

9a)采用基于任务导向的混合优化方法求解建立的目标函数，根据第k级离散网格上形成的系统矩阵的大小选择合适的优化方法，获得第k级离散网格上的靶向目标能量密度分布；

9b)利用靶向目标能量密度分布计算第k级离散网格边界节点上的光通量密度；

9c)判断第k级离散网格边界节点上的光通量密度的测量值和计算值之差，如果小于给定阈值，则目标函数求解过程结束，获得光学三维成像的靶向目标重建结果，转向步骤(10)；否则，转向步骤9d)；

9d)根据步骤9a)获得的靶向目标能量密度分布和步骤9b)获得的边界节点上的光通量密度计算值，调整第k+1级离散网格，转向步骤(7)；

(10)显示结果。

2.根据权利要求1所述的基于生物组织特异性的光学三维成像方法，其特征在于：所述步骤(1)中数据采集的实现是利用多模态分子影像系统等间隔采集不少于四个角度的荧光数据、不少于四个角度的激光数据、不少于360个角度的计算机断层成像投影数据。

3.根据权利要求1所述的基于生物组织特异性的光学三维成像方法，其特征在于：所述步骤(1)中的荧光数据、激光数据和计算机断层成像投影数据的初始采集角度相同。

4.根据权利要求1所述的基于生物组织特异性的光学三维成像方法，其特征在于：所述步骤6b)中光传输过程是采用扩散近似方程描述光在高散射特性生物组织中的传输过程，采用简化球谐波近似方程描述光在低散射特性生物组织中的传输过程，以及采用朗伯源特性方程描述光在空腔组织中的传输过程。

5.根据权利要求1所述的基于生物组织特异性的光学三维成像方法，其特征在于：所述步骤9a)中根据系统矩阵的大小选择合适的优化方法是通过采用如下规则实现的：对于小系统矩阵，采用求解速度较快的奇异值分解法进行求解；对于大系统矩阵，采用基于预处理系统矩阵的混合迭代法进行求解。

Images