CN101739503B - 一种预测光在生物组织中传播模型的实现方法 - Google Patents

Abstract

本发明涉及一种预测光在生物组织中传播模型的实现方法，通过使用自适应多重网格模型，根据在粗网格上求解结果的误差，自适应地剖分粗网格得到细网格；然后采用多重网格上的全V循环，在细网格上进行光滑迭代消除高频残量，然后将次结果限制到较粗网格上进行残量校正，再次延拓到细网格上后进行后光滑处理，直到最终达到所要求的精度。本发明有效降低了由于网格维数的增加而带来的计算代价，提高了收敛的速度，同时由于对低频误差的直接求解也提高了计算精度，能够有效的模拟对光学分子影像中光子从体内发射，经过多次散射达到表面的这一过程。

Description

一种预测光在生物组织中传播模型的实现方法

技术领域

本发明属于光学分子影像领域，涉及预测光在生物组织中传播的规律，尤其是一种快速准确的自适应多重网格模型的实现方法。

背景技术

基因组学、蛋白组学和疾病基因组学的迅速发展，为分子影像提供了全新的途径。光学分子影像技术就是这样一种可以在分子水平上实现生物组织内部物理过程无创实时动态成像，而准确预测光子在生物组织中的传播规律对于定性、定量研究组织的性质以及对其可视化无疑具有重要意义。

在单一固定网格上对波尔兹曼辐射传输方程的一阶球谐近似方程使用各种经典迭代法如雅克比、高斯-赛德尔等方法求解。根据泛函分析的知识，该方程的刚度矩阵M的条件数是O(h^-2)，其中h代表网格的尺寸。随着网格尺寸的不断变小，刚度矩阵的条件数将会越来越小，以至于无穷。当只是使用单一固定网格上的经典迭代法时，由于各种迭代法本质上都是低通滤波器，对高频误差有着很好的迭代效果，而对低频误差收敛速度却变得很慢；而且随着所使用的网格数的增多，收敛速度将越来越慢以至于不再收敛。因此，传统在单网格上的光滑迭代求解方法既不能保证效率，也不能保证一贯的精度。

发明内容

为了解决现有技术的问题，本发明的目的是提供一种预测光在生物组织中传播的自适应多重网格模型实现方法，通过自适应细分，得到多套不同尺度的网格，在粗细网格上进行全V循环。

为了实现上述目的，本发明涉及的预测光在生物组织中传播的自适应多重网格模型实现方法的技术方案采用多重网格法，步骤如下：

步骤1：首先将生物组织离散为粗网格，在粗网格上直接通过不完全乔利斯基分解方法求解有限元生成的矩阵方程M₀Ф₀＝b₀；

步骤2：计算矩阵结果Ф₀在粗网格上每个单元上的误差，并自适应地细分误差较大的单元，获得细分网格并计算；

步骤3：在细分网格上使用高斯-赛德尔方法对有限元生成的矩阵方程M_kФ_k＝b_k进行预光滑迭代 ${\mathrm{\Φ}}_k^l=\mathrm{Smooth}({\mathrm{\Φ}}_k^{l-1},b_k),$ l＝1，2.......v₁，将计算结果Ф_k的高频误差滤除，得到预光滑结果为

步骤4：计算预光滑结果

的残量

使用限制算子

把该残量映射到第k-1层粗网格上，获得第k-1层网格上的残量r_k-1表示如下： $r_{k-1}=I_k^{k-1}(b_k-M_k{\mathrm{\Φ}}_k^{v_1});$

步骤5：在获得的所有k个网格上递归：如果k≠0时，跳转至步骤3求解M_k-1Ф_k-1＝b_k-1，否则，若k≠L时执行步骤6；

步骤6：将生物组织离散的较粗网格上求解的结果使用延拓算子

映射回该组织离散后的细网格上，并且对细网格上预光滑结果进行校正： ${\mathrm{\Φ}}_k^{v_1}={\mathrm{\Φ}}_k^{v_1}+I_{k-1}^k{\mathrm{\Φ}}_{k-1};$

步骤7：再一次在生物组织离散的细网格上使用高斯-赛德尔方法对M_kФ_k＝b_k进行后光滑迭代 ${\mathrm{\Φ}}_k^l=\mathrm{Smooth}({\mathrm{\Φ}}_k^{l-1},b_k),$ l＝v₁₊₁，...v₁+v₂，将延拓所带来的高频误差滤除，然后计算残量的相对误差： ${\left|\right|b_k-M_k{\mathrm{\Φ}}_k^{v_1+v_2}\left|\right|}_2/{\left|\right|b_k\left|\right|}_2,$ 迭代若该值小于设定的全局误差，则停止计算；否则跳转至步骤2。

本发明的有益效果是：通过使用自适应多重网格模型，根据在粗网格上求解结果的误差，自适应地剖分粗网格得到细网格；然后采用多重网格上的全V循环，在细网格上进行光滑迭代消除高频残量，然后将次结果限制到较粗网格上进行残量校正，再次延拓到细网格上后进行后光滑处理，直到最终达到所要求的精度。有效降低了由于网格维数的增加而带来的计算代价，提高了收敛的速度，同时由于对低频误差的直接求解也提高了计算精度，能够有效的模拟对光学分子影像中光子从体内发射，经过多次散射达到表面的这一过程。克服了上述经典方法所存在的收敛速度慢甚至因不收敛得不到预测结果的问题，可以准确高效的预测光子在生物组织中的传输规律。而且对于多个感兴趣区域的情况，避免了单一固定网格上的迭代方法随着区域的增加而计算时间迅速增加的问题，本发明所涉及方法的计算时间增加不明显，这对于快速预测规律将是十分必要的。

附图说明

图1为本发明的框架流程图。

图2为的自适应多重网格方法四个主要模块的框架。

图3为多重网格算法的全V循环图。

图4为光学分子影像自适应多重网格算法的流程图。

图5为仿真实验所使用的仿体以及使用自适应多重网格方法所得到的表面的部分出射光强分布。

图6为在仿体表面上x＝13mm处一个探测圆环上的强度曲线。

图7为多个光源条件下单一网格方法和自适应多重网格方法的计算时间比较。

具体实施方式

下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

本发明是一种预测光子在生物组织中传播规律的模型。在光学分子影像中，例如通过对小动物转基因，使得小动物在病灶处产生荧光蛋白酶，进而通过静脉注射荧光底物，从而在该蛋白酶的催化作用下发射出可见或近红外光，光在生物组织中历经多次散射。最终有一部分到达小动物皮肤出射到自由空间中被探测器接收。由于哺乳动物组织对可见光及近红外光的强散射性，不像X-Ray在组织中直线传播，其传播规律极其复杂。因此，准确地预测这一过程对于定量并定位重建活体动物体内的荧光光源的逆问题是十分重要的。如图1示出本发明采用的自适应多重网格模型高效准确地预测了光在生物组织中传播的物理过程，具体实施步骤如下：

步骤1：首先将生物组织离散为粗网格，在粗网格上直接通过不完全乔利斯基分解方法求解有限元生成的矩阵方程M₀Ф₀＝b₀，其中：M₀为该粗网格所生成的刚度矩阵、Ф₀为粗网格表面上每个节点的辐射强度所构成的向量、b₀为粗网格上每个节点的光源强度所构成的向量；

步骤2：计算矩阵结果Ф₀在粗网格上每个单元上的误差，并自适应地细分误差较大的单元，获得细分网格并计算；

步骤3：在细分网格上使用高斯-赛德尔方法对有限元生成的矩阵方程M_kФ_k＝b_k进行预光滑迭代 ${\mathrm{\Φ}}_k^l=\mathrm{Smooth}({\mathrm{\Φ}}_k^{l-1},b_k),$ l＝1，2.......v₁，将计算结果Ф_k的高频误差滤除，得到预光滑结果为

M_k为第k层细分网格所生成的刚度矩阵、Ф_k为第k层细分网格表面上每个节点的辐射强度所构成的向量、b_k为第k层上每个节点的光源强度所构成的向量，Smooth表示预光滑迭代，k表示细分网格的层数，l表示第l次光滑迭代，v₁表示光滑迭代的次数；

步骤4：计算预光滑结果

的残量

使用限制算子

把该残量映射到第k-1层网格上： $r_{k-1}=I_k^{k-1}(b_k-M_k{\mathrm{\Φ}}_k^{v_1});$ r_k-1表示第k-1层网格上的残量；

步骤5：在获得的所有k个网格上递归：如果k≠0时，跳转至步骤3求解M_k-1Ф_k-1＝b_k-1，否则，若k≠L时执行步骤6；

步骤6：将生物组织离散的较粗网格上求解的结果使用延拓算子

映射回该组织离散后的细网格上，并且对细网格上预光滑结果进行校正： ${\mathrm{\Φ}}_k^{v_1}={\mathrm{\Φ}}_k^{v_1}+I_{k-1}^k{\mathrm{\Φ}}_{k-1};$

步骤7：再一次在生物组织离散的细网格上使用高斯-赛德尔方法对M_kФ_k＝b_k进行后光滑迭代 ${\mathrm{\Φ}}_k^l=\mathrm{Smooth}({\mathrm{\Φ}}_k^{l-1},b_k),$ l＝v₁₊₁，...v₁+v₂，将延拓所带来的高频误差滤除，然后计算残量的相对误差： ${\left|\right|b_k-M_k{\mathrm{\Φ}}_k^{v_1+v_2}\left|\right|}_2/{\left|\right|b_k\left|\right|}_2,$ 如果k＝L时，则判断残量的相对误差值小于设定的某一阈值，则停止计算；如果k≠L时，则执行步骤6；否则跳转至步骤2，其中v2表示后光滑的次数，L表示网格的总数。

在本发明方法中，计算引起的高频误差在细网格上处理，而低频误差则映射局限到较粗网格上处理。与传统的在单一网格上的方法相比较，既提高了计算精度，又大大减少了计算代价；而且由于在网格细分中采用自适应的策略，仅细分误差较大的区域，大大减少了问题的维数，进一步提高了计算的效率。加之使用Marshak部分反射边界条件，预测结果比在单一固定网格上的模型有了较大提高，与真实实验的结果以及蒙特卡洛统计结果吻合。

请参阅图2，使用计算机、C++和NetGen软件，实施了预测光在生物组织中传播模型的自适应多重网格方法，图2示出四个主要模块框架图，分别为多重网格的全V循环模块1、自适应网格剖分模块2、等级基函数计算模块3以及限制和延拓模块4。全V循环模块1是本发明方法的主体框架，从生物组织离散的最粗网格的计算结果开始，经过自适应网格剖分模块2，使用等级基函数计算模块3获得细分后网格上的节点基函数，进而通过限制和延拓模块4的限制算子将细网格上的结果映射到粗网格上，或使用延拓算子将粗网格上的结果映射到细网格上继续求解，直到达到所需的计算精度。

请参阅图3为多重网格算法的全V循环图。在最粗网格上直接使用不完全乔利斯基分解方法直接求解矩阵方程，而较细网格上使用光滑迭代，在粗网格与细网格之间的结果映射则使用延拓与限制。同时，细网格上的节点基函数空间由所有已经生成网格上的等级基函数空间正交形成，即 ${\mathrm{\Ω}}^k=v^1\⊕v^2\⊕...v^k,$ v^k表示第k层等级基函数所形成的空间，

为正交算子。

请参阅图4为本发明中光学分子影像自适应多重网格算法的流程图。该流程图表示一个完整V循环。从最细网格(k＝l)开始，先进行预光滑迭代v₁次： ${\mathrm{\Φ}}_k:={\mathrm{Smooth}}_k^{\left(v_1\right)}({\mathrm{\Φ}}_k,b_k),$ 再将光滑后的结果所得到的残量通过限制算子

映射到粗网格上 $b_k:=I_{k+1}^k(b_{k+1}-M_{k+1}{\mathrm{\Φ}}_{k+1}).$ 重复该操作(k:＝k-1)，直到到达最粗网格(k＝0)。在最粗网格上使用不完全乔利斯基分解直接求解 ${\mathrm{\Φ}}_0:=M_0^{-1}{b}_0.$ 然后再将该结果通过延拓算子映射到细网格上，并补偿原来在该层网格上得到的结果 $({\mathrm{\Φ}}_k:={\mathrm{\Φ}}_k+I_{k-1}^k{\mathrm{\Φ}}_{k-1})$ ，然后进行后光滑迭代 $({\mathrm{\Φ}}_k:={\mathrm{Smooth}}_k^{\left(v_2\right)}({\mathrm{\Φ}}_k,b_k)).$ 重复该操作(k:＝k+1)，直到到达最细网格(k＝l)。

请参阅图5为仿真实验所使用的仿体以及使用自适应多重网格方法所得到的表面的部分出射光强分布。仿体由六部分组成，分别模拟生物组织的心脏51、肝脏52、左肺53、右肺54、骨骼和肌肉55的非匀质特性。各个部分都提供了先验的光学特性参数。在(13，-4，-4)坐标处放置了一个半径为1mm的球形光源56，它的光强为238纳瓦/毫米³。在椭圆柱57表面上为自适应多重网格模型预测的出射光强。单位为纳瓦/毫米²。

请参阅图6为在仿体表面上x＝13mm处一个探测圆环(图5中的58)上的强度曲线。其中，实线为自适应多重网格模型的预测结果，虚线则为蒙特卡洛方法的结果。使用Marshak部分反射边界条件，本方法的预测结果与蒙特卡洛结果十分吻合，误差小于2％。

请参阅图7为多个光源条件下单一网格方法和自适应多重网格方法的计算时间比较。与在单一网格A上的计算时间相比较，自适应多重网格模型计算效率更高。而且，预光滑处理和后光滑处理，适合在多光源情况下对光在生物组织中传播进行数值仿真：随着光源个数的增加，自适应多重网格B模型的优势更明显：它的计算时间没有明显增加，而单一网格上的计算时间则大大增加。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (5)

1.一种预测光在生物组织中传播模型的实现方法，其特征在于：采用多重网格法，步骤如下：

步骤1：首先将生物组织离散为粗网格，在粗网格上直接通过不完全乔利斯基分解方法求解有限元生成的矩阵方程M₀Φ₀＝b₀，式中：M₀为该粗网格所生成的刚度矩阵、Φ₀为粗网格表面上每个节点的辐射强度所构成的向量、b₀为粗网格上每个节点的光源强度所构成的向量；

步骤2：计算矩阵结果Φ₀在粗网格上每个单元上的误差，并自适应地细分误差较大的单元，获得细分网格并计算；

步骤3：在细分网格上使用高斯-赛德尔方法对有限元生成的矩阵方程M_kΦ_k＝b_k进行预光滑迭代l＝1，2.......v₁，将计算结果Φ_k的高频误差滤除，得到预光滑结果为

式中：M_k为第k层细分网格所生成的刚度矩阵，Φ_k为第k层细分网格表面上每个节点的辐射强度所构成的向量，b_k为第k层上每个节点的光源强度所构成的向量，Smooth表示预光滑迭代，k表示细分网格的层数，l表示第l次光滑迭代，v₁表示光滑迭代的次数；

步骤4：计算预光滑结果

的残量

使用限制算子

把该残量映射到第k-1层粗网格上，获得第k-1层网格上的残量r_k-1表示如下： $r_{k-1}=I_k^{k-1}(b_k-M_k{\mathrm{\Φ}}_k^{v_1}),$ k：＝k-1；

步骤5：在获得的k个网格上递归：如果k≠0时，跳转至步骤3求解M_k-1Φ_k-1＝b_k-1，若k＝0时，执行步骤6；

步骤6：将生物组织离散的较粗网格上求解的结果使用延拓算子映射回该组织离散后的细分网格上，并且对细分网格上预光滑结果进行校正： ${\mathrm{\Φ}}_k^{v_1}={\mathrm{\Φ}}_k^{v_1}+I_{k-1}^k{\mathrm{\Φ}}_{k-1};$

步骤7：再一次在生物组织离散的细分网格上使用高斯-赛德尔方法对M_kΦ_k＝b_k进行后光滑迭代

l＝v₁+1，...v₁+v₂，将延拓所带来的高频误差滤除，然后计算残量的相对误差：

如果k＝L时，则判断残量的相对误差值小于设定的某一阈值，则停止计算，如果判断残量的相对误差值不小于设定的某一阈值，则跳转至步骤2；如果k≠L时，则执行步骤6；式中v₂表示后光滑的次数，L表示网格的总数。

2.根据权利要求1所述的预测光在生物组织中传播的模型的实现方法，其特征在于：采用限制算子

和延拓算子

实现中间计算结果在细分网格和粗网格之间映射。

3.根据权利要求1所述的预测光在生物组织中传播的模型的实现方法，其特征在于：所述细分网格的计算，包括使用所有已经生成网格上的等级基函数空间正交形成k层网格上的节点基函数空间Ω^k：v^k表示第k层等级基函数所形成的空间，

为正交算子。

4.根据权利要求1所述的预测光在生物组织中传播的模型的实现方法，其特征在于：所述预光滑处理和后光滑处理，适合在多光源情况下对光在生物组织中传播进行数值仿真。

5.根据权利要求1所述的预测光在生物组织中传播的模型的实现方法，其特征在于，还包括：在非匀质的生物组织中分别给定光学特性参数，并使用Marshak部分反射边界条件，预测生物组织表面出射光强与蒙特卡洛统计结果吻合。

Images