CN101221128A - 一种基于自适应有限元的多光谱重建方法 - Google Patents

Abstract

本发明公开一种基于自适应有限元的多光谱重建方法，包括步骤：基于自发荧光光源的多谱信息、重建目标区域的解剖结构以及光学特性参数的先验信息，结合自适应有限元理论，用后验可行光源区域选择的策略，进行基于自适应有限元的多光谱重建。本方法利用多种先验信息、多谱信息，结合自适应有限元的方法，提出后验可行光源区域的策略，该方法能有效地提高自发荧光断层成像的重建质量和速度，并大大降低了这一成像模态的重建不适定性，对现有的重建方法是一个很好的补充，在分子影像、重建算法等领域有着重要的应用价值。

Description

一种基于自适应有限元的多光谱重建方法

技术领域

本发明涉及到光学分子影像成像模态-自发荧光断层成像(BLT)技术，特别涉及到一种基于自适应有限元的多光谱重建方法。

背景技术

作为一种分子影像成像模态，自发荧光成像技术已经得到了迅速发展和广泛的应用。可是由于当前技术只能得到二维的平面图像，使当前研究人员不能三维的观测所感兴趣目标的变化，而这需要自发荧光断层成像(BLT)系统的发展来实现，由于自发荧光断层成像本身所具有的不适定特点，使确定自发光光源的位置重建在一般情况下没有唯一解，只有融入丰富的先验信息，才能够唯一定量的确定光源的信息。另外当在重建算法中融入先验信息时，如何合理利用这些信息也是需要进一步考虑的问题。

随着自发荧光断层成像的发展，当前有多种先验信息被用于光源的重建，其中包括重建区域的表面光强分布、重建区域的内部解剖结构、光学特性参数分布以及自发光光源的光谱特性和物理意义等等，它们有效的约束了光源重建。大体上可以把目前利用先验信息的重建算法分为两大类，一类是基于先验可行光源区域的方法，另一类是基于多光谱的重建方法。通过表面光强分布与区域内部解剖结构的结合，可以先验的对光源存在的区域进行大致推断，从而较好的确定一个小于整体区域的可行光源区域作为初始区域，来解决光源重建非唯一性问题。自发荧光光源的宽光谱特性也已经得到了很好的挖掘，由于不同材料对不同谱段的光有不同的吸收和散射特性，通过对不同谱段的探测，可以有利于光源的准确定位。值得一提的是目前有文献阐述如果只是考虑多光谱先验信息，而把非匀质重建目标考虑为匀质的重建目标，当进行光源重建时，特别是对于深度光源的重建是不充分的。

由于自发荧光断层成像是把整体的目标作为重建区域，这样势必增加了重建计算的负担，也不能很好的约束重建中的可行解。另外由于多光谱测量的使用以及非接触式探测方式的发展，都对当前重建算法的速度提出了严峻的挑战。在光学分子成像领域，虽然能够大大提高重建速度的解析方式已经得到发展，但是对于具有复杂内部结构的区域，这些方法还是不太合适。数值方法作为一种比较灵活的重建方法，已经被广泛的发展和应用，一些加速的方法也被进一步研究，如多网格方法，但是它需要对整体区域进行统一的网格细分，这对于自发荧光断层成像的重建产生巨大的计算负担，而且也是不必要的。

目前广泛采用的自发荧光断层成像重建算法是有限元方法，同时自发荧光断层成像重建的质量与边界测量数据的信噪比和重建目标区域的离散程度有很大的关系。所以在一定程度上，重建目标区域网格剖分得越细，自发荧光断层成像重建质量就越高。然而，如果网格剖分得过细，就会加重自发荧光断层成像问题的数值不稳定性和计算代价。相比一般的有限元方法，本发明所涉及的基于自适应有限元的多光谱重建方法利用多光谱测量数据、重建目标区域的解剖结构以及光学特性参数的先验信息来处理自发荧光断层成像问题的非唯一性和病态性；结合自适应有限元理论与后验可行光源区域选择策略，对自发荧光光源进行重建，不仅避免了由于增加多谱数据而带来的维数灾难问题，而且降低了自发荧光断层成像问题的病态性，提高了自发荧光断层成像重建质量。

发明内容

现有重建技术中，为解决自发荧光断层成像的不适定性存在利用多谱数据而带来的维数灾难问题以及自发荧光断层成像重建的速度和质量不高问题，为了解决的这些问题，本发明的核心思想是提出一种全新的基于自适应有限元多尺度多光谱重建方法。

本发明通过多种先验信息的融入，确保了光源重建的唯一性，自适应有限元方法的引入，在获得局部网格细分的基础上，也带来了后验可行光源区域选择的策略，这些不仅降低了重建的不适定性，而且改善了重建结果。

为了实现所述的目的，本发明基于自适应有限元多光谱重建方法的技术方案包括：基于自发荧光光源的多谱信息、重建目标区域的解剖结构以及光学特性参数的先验信息，结合自适应有限元理论，用后验可行光源区域选择的策略，进行基于自适应有限元的多光谱重建。

根据本发明的实施例，具体步骤包括以下：

步骤S1：在重建目标区域的第l-1(l＞1)层剖分网格上，利用有限元理论把单谱段扩散方程离散为线性方程；

步骤S2：基于后验可行光源区域的选择，在单谱段上建立未知光源变量与边界测量数据之间的线性关系，并利用自发光光源光谱性质，把单谱段组合成含有多个谱段的统一方程；

步骤S3：利用正则理论确立优化的目标函数，然后利用谱投影梯度基的大尺度优化算法对目标函数进行优化，以获得l-1层上的重建结果；

步骤S4：利用重建结果求解单谱段上的光密度分布，用评判准则评判重建是否停止，如果满足准则中的任意一个，重建停止；

步骤S5：如果不能满足评判准则的任意一个，利用后验的测度算法，对可行光源区域进行选择，并对可行光源区域和禁止光源区域进行自适应的网格细分；

步骤S6：从第l-1层向第l层转换，完成必要的参数调整，然后转到第一步，继续重建。

根据本发明的实施例，利用光源的谱特性和后验可行光源区域的策略，组合成含有多个谱段的统一方程包括：

${\mathrm{\Φ}}_l^{\mathrm{meas}}={F_{l}W}_l^S{S}_l$

其中Φ_l ^meas和S_l表示多谱下第l离散层上的边界测量值和光源未知量；F_l表示多谱下的边界测量值与光源未知量之间的关系；W_l ^S用于后验可行光源区域的选择；它们的具体形式如下：

${\mathrm{\Φ}}_l^{\mathrm{meas}}=\left[\begin{array}{c}{\mathrm{\Φ}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_1\right)\\ {\mathrm{\Φ}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_2\right)\\ .\\ .\\ .\\ {\mathrm{\Φ}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_K\right)\end{array}\right],$ $F_l=\left[\begin{array}{c}\mathrm{\ω}\left({\mathrm{wb}}_1\right)f_l\left(w{b}_1\right)\\ \mathrm{\ω}\left({\mathrm{wb}}_2\right)f_l\left({\mathrm{wb}}_2\right)\\ .\\ .\\ .\\ \mathrm{\ω}\left({\mathrm{wb}}_K\right)f_l\left({\mathrm{wb}}_K\right)\end{array}\right]$

${W_l}^S=\mathrm{Diag}(w_{l\left(11\right)}^S,w_{l\left(22\right)}^S,...,w_{l\left(\mathrm{ii}\right)}^S,...,w_{l\left(N_{\mathrm{Pl}}{N}_{\mathrm{Pl}}\right)}^S),$ $w_{l\left(\mathrm{ii}\right)}^S=\left\{\begin{array}{cc}0& \{{S^{\&prime;}}_{l\left(i\right)}<{\mathrm{\&gamma;}}_{\mathrm{sp}}^{\left(l\right)}{S_l}^{\&prime;\max },l>1\}\\ 1& \{l=1\}\mathrm{or}\{{S^{\&prime;}}_{l\left(i\right)}\&GreaterEqual;{\mathrm{\&gamma;}}_{\mathrm{sp}}^{\left(l\right)}{S_l}^{\&prime;\max },l>1\}\end{array}\right.$

其中，ω(wb_k)是离散的第k谱段在整个自发荧光光源谱段上占有的能量比，满足ω(wb_k)≥0和 $\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\mathrm{\&omega;}\left({\mathrm{wb}}_k\right)\&ap;1;$ Φ_l ^meas(wb_k)和f_l(wb_k)表示在第l层上第k谱段上的边界测量值和光源未知量；s′_l(i)和s_l′^max是l-1层重建结果在网格细分的基础上插值得到的l层上的初始重建结果以及结果中的最大值，即：

${S^{\&prime;}}_l=I_{l-1}^l{S}_{l-1}^r,(l>1)$

γ_sp ^(l)是比例因子，用来确定所要选择在l层上的可行光源区域的大小。

根据本发明的实施例，采用三种测度准则来评判重建是否停止：单谱段上边界测量值和计算值之间的误差、每一层优化中的梯度g_Θl ^(l)的大小和进行重建层数l的上限，具体表示如下：

$\left|\right|{\mathrm{\&Phi;}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_k\right)-{\mathrm{\&Phi;}}_l^c\left({\mathrm{wb}}_k\right)\left|\right|<{\mathrm{\&epsiv;}}_{\mathrm{\&Phi;}}$ 、 $\left|\right|g_{{\mathrm{\&Theta;}}_l}^{\left(l\right)}\left|\right|<{\mathrm{\&epsiv;}}_g$ 和l≥L_max，其中Φ_l ^meas(wb_k)是第l离散层第k谱段上的边界测量值；Φ_l ^c(wb_k)是相应的计算值；ε_Φ、ε_g和L_max分别是相应的测度阈值。

根据本发明的实施例，在禁止光源区域和可行光源区域，采用等级值亏修正基的后验估计技术和直接最大值选择技术对要进行细分的网格进行选取。

根据本发明的实施例，在禁止光源区域，仅利用单谱段上的误差估计来实现单元选择，用于提高重建速度。

本发明的技术效果：随着自发荧光断层成像模态的发展，有许多问题需要被解决，本发明提出了基于自适应有限元的多光谱重建方法，基于扩散近似模型，该算法考虑了重建区域的非均匀特性，也考虑自发荧光光源的谱特点。通过多光谱测量和自适应网格的进化，不仅有效的降低了自发荧光断层成像的不适定性，而且也避免了由于利用多谱数据而带来的维数灾难问题。另外自发荧光断层成像重建的速度和质量进一步从后验可行光源区域的选择中进一步提高。

本发明利用自适应有限元方法综合了当前主流的两种自发荧光断层成像重建方法，发展了后验可行光源区域的策略，光源的唯一性得到了很好的解决，重建的速度和质量由于多尺度方法的引进被进一步提高，此方法对自发荧光断层成像的发展具有重要的应用价值。

附图说明

图1本发明方法的流程图。

图2通过微核磁共振扫描的小鼠模型。

图3基于单光谱和多光谱的重建结果，(a)和(b)是单光谱的重建结果，光源分别在小鼠模型的半径中心位置和小鼠模型中心位置；(c)和(d)是相应多光谱的重建结果。

图4重建结果进化示意图。

图5考虑光源特性参数误差的重建结果，(a)和(b)是光源在小鼠模型的半径中心位置，光学特性误差分别为正、负50％时的重建结果；(c)和(d)是光源在小鼠模型中心位置，相应的重建结果。

具体实施方式

下面将结合附图详细描述本发明的重建方法，应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

本发明的重建方法包括以下步骤：

(1)在l-1(l＞1)层上，对重建区域进行网格剖分，利用有限元方法在单谱段上把扩散方程离散为线性方程，以便于进行后续的处理；

(2)考虑到未知光源与边界测量数据的线性关系，在单谱段上，通过矩阵变换建立两者的关系，然后考虑自发荧光光源的光谱特点和后验可行光源区域选择，把以上多个单谱上的关系式组合为一个表征多谱特点的方程；

(3)利用正则的方法来确立优化的目标函数，然后利用谱投影梯度基的大尺度优化算法对目标函数进行优化，来获得l-1层上的重建结果；

(4)通过最大迭代次数和优化梯度降低程度判断优化过程是否完成，当优化完成后，利用上一步的重建结果求解单谱段上的光密度分布，然后利用多种评判准则评判重建是否停止，如果满足准则中的任意一个，重建停止；

(5)如果所有的评判准则中的任意一个均不能被满足，利用后验的测度算法，对可行光源区域进行选择，并利用后验的误差估计方法对可行和禁止光源区域进行自适应的网格细分；

(6)从第l-1层向第l层转换，完成必要的参数调整，然后转到第一步，继续重建。

下面对本发明方法涉及的关键步骤进行逐一详细说明，具体形式如下面所述：

图1示出了本发明方法的流程图。

步骤1：利用有限元方法转化扩散方程成为线性方程

精确描述光子在物质中传输的数学模型是辐射传输方程，但是由于此方程是一个积分-微分方程，因此，利用此方程进行光源重建将需要花费许多时间，幸运的是，大部分物质都具有强散射的特性，在这种情况下，扩散方程作为一种近似可以很好的对光子传输进行模拟。当自发荧光成像实验在一个全黑的环境下进行时，光源发光可以暂时认为是稳定的，稳态的扩散方程就可以满足模拟条件，进一步考虑到不同谱段对重建物质的光学特性参数的影响，可以得到以下扩散方程形式：

-·(D(x，λ)Φ(x，λ))+μ_a(x，λ)Φ(x，λ)＝S(x，λ)(x∈Ω)

这里Ω是目标重建区域；Φ(x，λ)是光子流密度；S(x，λ)是光源密度；μ_a(x，λ)是生物组织的吸收系数；D(x，λ)＝1/(3(μ_a(x，λ)+(1-g)μ_s(x，λ)))是扩散系数，其中μ_s(x，λ)是散射系数，g是各向异性参数。进一步考虑到重建区域折射系数n和外部媒介折射系数n′的非匹配性，扩散方程的边界条件可以被表达为：

Φ(x，λ)+2A(x；n，n′)D(x，λ)(v(x)·Φ(x，λ))＝0(x∈Ω)

其中v(x)是重建区域边界Ω的单位法向量；A(x；n，n′)≈(1+R(x))/(1-R(x))，当外界媒质是空气时，R(x)可以近似为R(x)≈-1.4399n^-2+0.7099n^-1+0.6681+0.0636n。而表面能够测量到的量是流出光流密度Q(x，λ)，即为：

$Q(x,\mathrm{\&lambda;})=-D(x,\mathrm{\&lambda;})(v\&CenterDot;\&dtri;\mathrm{\&Phi;}(x,\mathrm{\&lambda;}))=\frac{\mathrm{\&Phi;}(x,\mathrm{\&lambda;})}{2A(x;n,n^{\&prime;})}(x\&Element;\&PartialD;\mathrm{\&Omega;})$

在实际的实验中，需要通过滤波片对探测光信号进行谱段分离，因此在自发荧光光源谱段范围被离散为几个谱段，即wb_k∈[λ_k，λ_k+1]，k＝1，2，...，K-1

在自适应有限元分析的框架下，基于自适应的网格细分，可以把重建区域网格剖分表示为{T₁，...T_l，...}的网格序列，其中网格剖分T_l包括N_Tl个离散单元和N_Pl个离散点。利用有限元方法，对以上扩散方程进行离散，可以得到以下线性方程：

(K_l(wb_k)+C_l(wb_k)+B_l(wb_k))Φ_l(wb_k)＝F_l(wb_k)S_l(wb_k)

步骤2：建立源变量与表面测量值之间的线性关系，利用光源谱信息重组单谱段方程；

让M_l(wb_k)＝K_l(wb_k)+C_l(wb_k)+B_l(wb_k)，M_l(wb_k)是一个对称正定矩阵。考虑未知源变量S_l(wb_k)与边界测量值Φ_l ^meas(wb_k)之间的线性关系，有

${\mathrm{\&Phi;}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_k\right)=f_l\left({\mathrm{wb}}_k\right)S_l\left({\mathrm{wb}}_k\right)$

其中f_l(wb_k)通过移出M_l(wb_k)^-1F_l(wb_k)中对应非测量点的行而得到。进一步考虑光源的谱特性，在谱段wb_k上所占整个谱段的能量可以表示为S(wb_k)＝ω(wb_k)S，其中ω(wb_k)≥0， ${\mathrm{\&Sigma;}}_{k=1}^K\mathrm{\&omega;}\left({\mathrm{wb}}_k\right)\&ap;1;$ S表示总的光源密度。考虑能量谱分布和可行光源区域，可以得到：

${\mathrm{\&Phi;}}_l^{\mathrm{meas}}=F_l{W}_l^S{S}_l$

其中Φ_l ^meas和S_l表示多谱下的边界测量值和光源未知量；F_l表示多谱下的边界测量值与光源位置量之间的关系；W_l ^S用于后验的选择可行光源区域，它们其中的部分具体形式如下：

${\mathrm{\&Phi;}}_l^{\mathrm{meas}}=\left[\begin{array}{c}{\mathrm{\&Phi;}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_1\right)\\ {\mathrm{\&Phi;}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_2\right)\\ .\\ .\\ .\\ {\mathrm{\&Phi;}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_K\right)\end{array}\right],$ $F_l=\left[\begin{array}{c}\mathrm{\&omega;}\left({\mathrm{wb}}_1\right)f_l\left({\mathrm{wb}}_1\right)\\ \mathrm{\&omega;}\left({\mathrm{wb}}_2\right)f_l\left({\mathrm{wb}}_2\right)\\ .\\ .\\ .\\ \mathrm{\&omega;}\left({\mathrm{wb}}_K\right)f_l\left({\mathrm{wb}}_K\right)\end{array}\right]$

${W_l}^S=\mathrm{Diag}(w_{l\left(11\right)}^S,w_{l\left(22\right)}^S,...,w_{l\left(\mathrm{ii}\right)}^S,...,w_{l\left(N_{\mathrm{Pl}}{N}_{\mathrm{Pl}}\right)}^S),$ $w_{l\left(\mathrm{ii}\right)}^S=\left\{\begin{array}{cc}0& \{{S^{\&prime;}}_{l\left(i\right)}<{\mathrm{\&gamma;}}_{\mathrm{sp}}^{\left(l\right)}{S_l}^{\&prime;\mathrm{msx}},l>1\}\\ 1& \{l=1\}\mathrm{or}\{{S^{\&prime;}}_{l\left(i\right)}\&GreaterEqual;{\mathrm{\&gamma;}}_{\mathrm{sp}}^{\left(l\right)}{S_l}^{\&prime;\max },l>1\}\end{array}\right.$

其中s′_l(i)和s_l′^max是l-1层重建结果在网格细分的基础上插值得到的l层的初始重建结果以及结果中的最大值，即：

${S^{\&prime;}}_l=I_{l-1}^l{S}_{l-1}^r,(l>1)$

γ_sp ^(l)是比例因子，用来确定所要选择在l层上的可行光源区域的大小。通过保留F_lW_l ^S中对应可行光源区域的列，可以得到最后的未知源变量与边界条件在第l层的关系为：

$A_l{S}_l^p={\mathrm{\&Phi;}}_l^{\mathrm{meas}}$

步骤3：利用正则理论得到优化目标函数，用优化方法进行优化，得到重建结果；

步骤2虽然得到了未知源变量与边界测量值之间的关系，由于A_l矩阵的病态，不能通过直接求解的方法得到重建结果，利用正则理论，并考虑未知源变量的物理意义，可以得到第l层的优化目标函数Θ_l(S_l ^p)为：

$\underset{0\&le;S_l^p\&le;S_l^{\sup }}{\min }{\mathrm{\&Theta;}}_l\left(S_l^p\right)=\{{\left|\right|A}_l{S}_l^p-{\mathrm{\&Phi;}}_l^{\mathrm{meas}}{\left|\right|}_{\mathrm{\&Lambda;}}+{\mathrm{\&lambda;}}_1{\mathrm{\&eta;}}_1\left(S_l^p\right)\}$

其中S_l ^sup是第l层源密度的上界；Λ是权重矩阵，‖V‖_Λ＝V^TΛV；λ_l是正则参数；η_l(·)是惩罚函数，通过选择有效的优化方法，可以获得较好的自发荧光断层重建。这里利用的是修正的谱梯度级的大尺度优化算法，在判断优化过程是否中止时，采用了当前优化梯度g_Θl ⁱ与初始优化梯度g_Θl ⁰的比例γ_gi ^(l)以及优化迭代次数N_i ^(l)作为评判准则，即当 ${\mathrm{\&gamma;}}_{g_i}^{\left(l\right)}<{\mathrm{\&epsiv;}}_g^{\left(l\right)}$ 或 $N_i^{\left(l\right)}>N_{\max }^{\left(l\right)}$ 时，优化停止，得到需要的重建结果，其中 ${\mathrm{\&gamma;}}_{g_i}^{\left(l\right)}=\left|\right|g_{{\mathrm{\&Theta;}}_l}^i\left|\right|/\left|\right|g_{{\mathrm{\&Theta;}}_l}^0\left|\right|,$ N_max ^(l)是最大允许迭代次数。

步骤4：基于上一步重建结果，计算光在可测量点上的分布，判断重建是否停止，若不停止，利用后验测度选取网格，然后对其细分；

当第l层的优化停止后，利用优化的重建结果，求解单谱段上光在可测量点上的分布Φ_l ^c(wb_k)，然后采用了三种测度准则来评判重建是否停止，分别是单谱段上边界测量值与计算值之间的误差、每一层优化中的梯度g_Θl ^(l)的大小和进行重建层数l的上限，具体表示如下：

$\left|\right|{\mathrm{\&Phi;}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_k\right)-{\mathrm{\&Phi;}}_l^c\left({\mathrm{wb}}_k\right)\left|\right|<{\mathrm{\&epsiv;}}_{\mathrm{\&Phi;}}$ 、 $\left|\right|g_{{\mathrm{\&Theta;}}_l}^{\left(l\right)}\left|\right|<{\mathrm{\&epsiv;}}_g$ 和l≥L_max，其中Φ_l ^meas(wb_k)是单谱段上的边界测量值；Φ_l ^c(wb_k)是在l层上的计算值；ε_Φ、ε_g和L_max分别是相应的测度阈值。当重建不停止时，利用重建的结果选取新的可行光源区域，然后在可行和禁止光源区域，通过两种不同的方法选取要细分的单元，分别是直接最大值选择方法和等级值亏修正基的误差估计方法。由于在可行光源区域，重建结果值相对大的单元表达了真实的光源位置，选择它们作为细分的对象来进一步改善重建的质量。考虑到计算量问题，虽然光源多谱特性作为先验信息被利用，可是由于每个谱段都是在同一生物组织内的光子传输，因此它们的传输模型是一致的，计算误差分布也是一致的，利用单谱段上的等级值亏修正基的估计技术在禁止光源区域来选取要细分的单元，从而减少了计算量，即：

${\tilde{e}}_{V_l}^{W_l}\left({\mathrm{wb}}_k\right)=D_{E_l}^{-1}\left({\mathrm{wb}}_k\right)r_{E_l}\left({\mathrm{wb}}_k\right)$

是近似的误差指示器；E_l是通过分解相应于V_l的二次元空间W_l而得到的一个子空间；D_El(wb_k)是对M_El(wb_k)的近似；r_El(wb_k)表不在E_l空间迭代求解的残留值。

步骤5：对离散区域进行局部网格细分；

相比于全局的网格细分方法，进行局部的网格细分有它自己的更加复杂的方式，这关系到区域离散基本单元的选择和细分规则的确立。由于三角形和四面体是较常用的描述复杂区域的基本单元形状，这样更适于实际重建目标的执行，因此选择选择四面体作为基本的区域离散单元。局部网格细分技术中关键一点是对在细分中产生的影响有限元计算的悬点进行处理，在这里选择“红-绿”结合的细分规则，首先通过“红”的方法对选定的四面体单元进行细分，获得八个子四面体单元。然后通过“绿”的封闭方法对在“红”规则细分中产生的悬点，通过对邻近单元细分的方法予以处理，从而最终达到局部网格细分的目的。

步骤6：进行参数转换，转入第1步，进入下一层重建；

通过以上的步骤，就完成了一层上的重建，为了转入下一层的重建，需要对依赖于层数的参数进行调整，然后转入第1步继续进行重建。

运行结果

为了验证本发明的方法，我们利用微核磁共振成像(MicroMRI)建立的整体数字鼠(如图2所示)和蒙特卡罗方法生成的边界测量数据进行了实验。

在此实验中，我们把标记PC3M-Luc细胞的荧光素酶催化发出的自发光源为重建目标，在重建实验中，这个光源大小为半径1mm，光源能量为1.0 nano-Watts。它的波谱范围从500nm到750nm，为了进行多光谱实验，我们依据当前的探测水平，把这一波谱范围分成了5个谱段进行分别探测。至于数字鼠模型和边界测量数据的生成方面，由于脑部光学参数的缺少，我们只是取整体小鼠颈部一直到尾部的部分作为生成测量数据的模型，利用我们自己研发的基于蒙特卡罗方法的仿真平台MOSE来产生表面测量数据。

为了对方法进行评估，多谱的MOSE平台利用离散为包含大约14600个三角面片的小鼠去获得表面测量数据。在仿真中，半径为1.0mm，能量为1.0 nano-Watts总能量的实球体光源的每一个谱段被采样1.0×10⁶个光子。另外，上界S_l ^sup为一个充分大的正数、权重矩阵Λ为单位阵和惩罚函数η_l(X)＝X^TX。梯度容忍比率ε_g ^(l)和最大迭代次数N_max ^(l)在每一层都是一样的，即为1.0×10^-5和5000。重建终止梯度范数ε_g、停止域值ε_Φ和最大细分层数L_max设为7.0×10^-9、1.0×10^-7和3。在最粗糙的网格层上，整个的重建区域为可行光源区域，当网格细分后，选择因子γ_sp ^(l)被初始设置为10^-4。禁止和可行光源区域的网格细分比例为1.5％和50％重建过程设置未知变量S_init ^(l)的初始猜想为1.0×10^-6。所有的重建过程在奔腾四2.8GHz和1GB内存的个人电脑上执行。

下面的附表是：

重建中利用的小鼠不同组织的在不同谱段的光学特性参数

波段	500～550nm		550～600nm		600～650nm		650～700nm		700～750nm
											μa *	μs’*	μa	μs’	μa	μs’	μa	μs’	μa	μs’
	骨头	9.0e-3	3.34	2.2e-3	2.93	6.0e-4	2.61	3.6e-4	2.34	5.9e-4											2.11
脂肪											6.2e-4	1.34	2.5e-4	1.28	2.1e-4	1.22	3.4e-4	1.18	1.2e-3	1.14
	心脏	9.1e-3	1.28	2.2e-3	1.13	6.9e-4	1.00	5.8e-4	0.91	1.4e-3											0.82
肠											1.7e-3	1.47	5.3e-4	1.32	2.8e-4	1.19	3.8e-4	1.09	1.3e-3	0.99
	肾	1.0e-2	3.04	2.6e-3	2.66	8.6e-4	2.36	8.0e-4	2.11	2.2e-3											1.90
肝脏											5.4e-2	0.83	1.2e-2	0.76	3.4e-3	0.70	2.0e-3	0.65	3.0e-3	0.60
	肺	2.7e-2	2.41	7.1e-3	2.30	2.0e-3	2.21	1.4e-3	2.12	2.8e-3											2.04
胰腺											1.0e-2	3.04	2.6e-3	2.66	8.6e-4	2.36	8.0e-4	2.11	2.2e-3	1.90
	脾	5.4e-2	0.83	1.3e-2	0.76	3.4e-3	0.70	2.0e-3	0.65	3.0e-3											0.60
胃											1.8e-3	1.74	6.0e-4	1.60	3.8e-4	1.48	5.7e-4	1.38	2.0e-3	1.29

*单位：mm^-1

重建结果如图3所示，在重建中，我们首先基于提出的方法和单光谱进行重建实验，由于没有利用多光谱的先验信息，对在小鼠模型半径中心位置以及小鼠模型中心位置的光源位置进行重建时，不能得到正确的重建结果，而且随着光源深度的增加，光源重建的不正确性越来越大，单光谱的重建结果，光源分别在小鼠模型的半径中心位置和小鼠模型中心位置，如图3(a)和3(b)所示。当把多光谱信息融入到整个重建中去以后，对相同位置的光源重建都得到了较好的结果，重建结果被显示在图3(c)和3(d)中。在重建中，共进行了两次网格的细分。光源位于小鼠模型半径中心位置时重建的网格进化显示在图4(a)、4(b)和4(c)中。

为了进一步对我们的方法进行验证，我们考虑光学特性参数的变化进一步进行了重建实验。由于在实际的活体小鼠的实验中，需要对光学特性参数进行实时的测量，这样不可避免会引入测量误差。在这个实验中，我们对光学特性参数扰动正负50％的误差。利用提出本发明的方法进行重建实验，当光源位置在小鼠模型半径中心位置上时，不论是高估还是低估光学特性参数，所得重建结果与精确的光学特性参数时的重建基本一致的。光源在小鼠模型半径中心位置时，光学特性误差分别为正负50％时的重建结果被显示在图5(a)和5(b)。随着光源深度的增加，光学特性参数的误差越来越影响光源的重建，图5(c)和5(d)显示了光源在小鼠模型中心位置时的重建结果。从图中可以看出，虽然在对小鼠模型中心位置的光源进行重建时有少许的位置误差，但是基本也对光源进行了较好的重建，这进一步显示了本发明所提方法的鲁棒性。

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

Claims (6)

1.一种基于自适应有限元的多光谱重建方法，其特征在于：

基于自发荧光光源的多谱信息、重建目标区域的解剖结构以及光学特性参数的先验信息，结合自适应有限元理论，用后验可行光源区域选择的策略，进行基于自适应有限元的多光谱重建。

2.根据权利要求1所述的方法，其特征在于，包括以下步骤：

步骤S1：在重建目标区域的第l-1(l＞1)层剖分网格上，利用有限元理论把单谱段扩散方程离散为线性方程；

步骤S2：基于后验可行光源区域的选择，在单谱段上建立未知光源变量与边界测量数据之间的线性关系，并利用自发光光源光谱性质，把单谱段组合成含有多个谱段的统一方程；

步骤S3：利用正则理论确立优化的目标函数，然后利用谱投影梯度基的大尺度优化算法对目标函数进行优化，以获得l-1层上的重建结果；

步骤S4：利用重建结果求解单谱段上的光密度分布，用评判准则评判重建是否停止，如果满足准则中的任意一个，重建停止；

步骤S5：如果不能满足评判准则的任意一个，利用后验的测度算法，对可行光源区域进行选择，并对可行光源区域和禁止光源区域进行自适应的网格细分；

步骤S6：从第l-1层向第l层转换，完成必要的参数调整，然后转到第一步，继续重建。

3.如权利要求2所述的方法，其特征在于，利用光源的谱特性和后验可行光源区域的策略，组合成含有多个谱段的统一方程包括：

${\mathrm{\&Phi;}}_l^{\mathrm{meas}}=F_l{W_l}^S{S}_l$

其中Φ_l ^meas和S_l表示多谱下第l离散层上的边界测量值和光源未知量；F_l表示多谱下的边界测量值与光源未知量之间的关系；W_l ^S用于后验可行光源区域的选择；它们的具体形式如下：

${\mathrm{\&Phi;}}_l^{\mathrm{meas}}=\left[\begin{array}{c}{\mathrm{\&Phi;}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_1\right)\\ {\mathrm{\&Phi;}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_2\right)\\ .\\ .\\ .\\ {\mathrm{\&Phi;}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_K\right)\end{array}\right],$ $F_l=\left[\begin{array}{c}\mathrm{\&omega;}\left(w{b}_1\right)f_l\left({\mathrm{wb}}_1\right)\\ \mathrm{\&omega;}\left({\mathrm{wb}}_2\right)f_l\left({\mathrm{wb}}_2\right)\\ .\\ .\\ .\\ \mathrm{\&omega;}\left({\mathrm{wb}}_K\right)f_l\left({\mathrm{wb}}_K\right)\end{array}\right]$

${W_l}^S=\mathrm{Diag}(w_{l\left(11\right)}^S,w_{l\left(22\right)}^S,...,w_{l\left(\mathrm{ii}\right)}^S,...,w_{l\left(N_{\mathrm{Pl}}{N}_{\mathrm{Pl}}\right)}^S),$ $w_{l\left(\mathrm{ii}\right)}^S\left\{\begin{array}{cc}0& \{{S^{\&prime;}}_{l\left(i\right)}<{\mathrm{\&gamma;}}_{\mathrm{sp}}^{\left(l\right)}{S_l}^{\&prime;\max },l>1\}\\ 1& \{l=1\}\mathrm{or}\{{S^{\&prime;}}_{l\left(i\right)}\&GreaterEqual;{\mathrm{\&gamma;}}_{\mathrm{sp}}^{\left(l\right)}{S_l}^{\&prime;\max },l>1\}\end{array}\right.$

其中，ω(wb_k)是离散的第k谱段在整个自发荧光光源谱段上占有的能量比，满足ω(wb_k)≥0和 $\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\mathrm{\&omega;}\left({\mathrm{wb}}_k\right)\&ap;1;$ Φ_l ^meas(wb_k)和f_l(wb_k)表示在第l层上第k谱段上的边界测量值和光源未知量；s′_l(i)和s_l′^max是l-1层重建结果在网格细分的基础上插值得到的l层上的初始重建结果以及结果中的最大值，即：

${S^{\&prime;}}_l=I_{l-1}^l{S}_{l-1}^r,(l>1)$

γ_sp ^(l)是比例因子，用来确定所要选择在l层上的可行光源区域的大小。

4.如权利要求2所述的方法，其特征在于，采用三种测度准则来评判重建是否停止：单谱段上边界测量值和计算值之间的误差、每一层优化中的梯度g_Θl ^(l)的大小和进行重建层数l的上限，具体表示如下：

$\left|\right|{\mathrm{\&Phi;}}_l^{\mathrm{meas}}\left({\mathrm{wb}}_k\right)-{\mathrm{\&Phi;}}_l^c\left({\mathrm{wb}}_k\right)\left|\right|<{\mathrm{\&epsiv;}}_{\mathrm{\&Phi;}}$ 、 $\left|\right|g_{{\mathrm{\&Theta;}}_l}^{\left(l\right)}\left|\right|<{\mathrm{\&epsiv;}}_g$ 和l≥L_max，其中Φ_l ^meas(wb_k)是第l离散层第k谱段上的边界测量值；Φ_l ^c(wb_k)是相应的计算值；ε_Φ、ε_g和L_max分别是相应的测度阈值。

5.如权利要求2所述的方法，其特征在于，在禁止光源区域和可行光源区域，采用等级值亏修正基的后验估计技术和直接最大值选择技术对要进行细分的网格进行选取。

6.如权利要求2所述的方法，其特征在于，在禁止光源区域，仅利用单谱段上的误差估计来实现单元选择，用于提高重建速度。

Images