CN101292863A - 一种基于单谱或混合段测量的自适应有限元光源重建方法 - Google Patents

Abstract

本发明公开一种基于单谱或混合谱测量的自适应有限元光源重建方法，利用重建目标区域的光学特性参数和解剖结构信息确定出先验可行光源区域和先验禁止光源区域；根据重建目标区域边界上的单谱或混合谱光通量测量数据、结合多级自适应有限元以及后验误差估计方法，对局部网格进行优化分解，以得到重建目标区域内的荧光光源。因此，本发明在单谱或混合谱边界测量数据的基础上，采用先验可行光源区域降低了自发荧光断层成像的病态性，使其能确定出唯一解，利用自适应有限元方法对网格划分程度和算法效率以及数值稳定性进行了均衡，并根据后验误差估计对局部网格进行优化，提高了荧光光源重建的质量、效率以及鲁棒性。

Description

一种基于单谱或混合段测量的自适应有限元光源重建方法

技术领域

本发明属于分子影像领域，涉及在体荧光断层成像算法，尤其是一种基于单谱或混合谱测量的自适应有限元光源重建方法。

背景技术

自发荧光断层成像技术是近年来新兴的一种光学分子影像技术，自发荧光是利用荧光素酶标记重建目标，在ATP以及氧气存在的情况下，若荧光素酶遇到底物荧光素，荧光素酶将会催化荧光素的氧化反应并产生光子，产生光的强度与标记目标的数量成正比。在重建目标区域外，利用高灵敏度的光学检测仪器，可以直接探测到逃逸出重建目标区域的光子，然后利用有效的荧光光源重建算法，就可以得到重建目标区域内部荧光光源的位置和强度。

自发荧光断层成像属于逆向问题，具有严重病态性的特点，其原因在于重建目标区域具有强散射特性，光子在其内部的传输不再沿直线传播，而是经过大量的散射过程，从而远远偏离了原来的运行轨道。同时，利用光学检测仪器在重建目标区域边界处测量得到的光通量是各边界点的离散数值，数量有限，而要求解区域内部点的数量非常巨大，所以自发荧光断层成像是一个非适定问题，病态性非常严重，其解不唯一并且对测量误差及噪声非常敏感。因此，如何能精确的重建出目标区域内部的荧光光源，就成为自发荧光断层成像中的核心问题。在有限元分析中，自发荧光断层成像光源重建的质量不仅与边界测量数据的信噪比有关，而且还与重建目标区域的离散化程度有关。在很大程度上，网格剖分得越细，重建质量就越好。然而，过细剖分网格会加重自发荧光断层成像的病态性，增加数值的不稳定和计算的代价。

发明内容

为了解决自发荧光断层成像存在的病态性和解不唯一性，本发明利用重建目标区域的光学特性参数和解剖结构信息，将重建目标区域划分为先验可行光源区域和先验禁止光源区域，降低自发荧光断层成像的病态性并可确定出唯一解；为了解决现有重建技术过细剖分网格会加重自发荧光断层成像的病态性、增加数值的不稳定和计算的代价的问题，本发明结合光学特性参数和先验可行光源区域的先验知识并根据后验误差估计结果，对局部有限元网格进行自适应优化，均衡了网格划分程度和算法效率，不但能有效提高光源位置和强度的重建精度，而且还能提高重建算法的效率和鲁棒性。在自适应过程中，本发明采用一种改进的牛顿法和多水平、动态策略来优化算法，提供基于单谱或混合段测量的自适应有限元光源重建方法。

为了实现所述的目的，本发明基于单谱或混合段测量的自适应有限元光源重建方法的技术方案包括如下步骤：

步骤S1：利用重建目标区域的光学特性参数和解剖结构信息确定出先验可行光源区域和先验禁止光源区域；

步骤S2：根据重建目标区域边界上的单谱或混合谱光通量测量数据、结合多级自适应有限元以及后验误差估计方法，对局部网格进行优化分解，以得到重建目标区域内的荧光光源。

根据本发明的实施例，将重建目标区域划分为先验可行光源区域和先验禁止光源区域，用于降低自发荧光断层成像的病态性，使其能确定出唯一解。

根据本发明的实施例，利用光学检测仪器在重建目标区域边界上，测量光通量Q(r)并求取光通量密度Φ(r)。

根据本发明的实施例，对重建目标区域Ω进行第k层网格剖分，利用有限元理论方法将扩散方程转换成线性矩阵方程。

根据本发明的实施例，利用正则理论形成第k层网格剖分的目标函数Θ_k(S_k ^p)：

$\min S_{\inf }^k\≤S_k^p\≤S_{\sup }^k\{{\left|\right|A_k{S}_k^p-{\mathrm{\Φ}}_k^m\left|\right|}_{\mathrm{\Λ}}+{\mathrm{\λ}}_k{\mathrm{\η}}_k\left(S_k^p\right)\}$

其中，S_inf ^k和S_sup ^k分别表示第k层网格剖分的光源最小和最大边界；Λ为加权矩阵；λ_k为正则化参数；η_k(S_k ^p)为补偿函数。

根据本发明的实施例，优化目标函数Θ_k(S_k ^p)进行迭代求解，得到第k层网格剖分的光源重建结果。

根据本发明的实施例，局部网格细分步骤如下：

(1)利用第k网格剖分的光源重建结果计算光通量密度Φ_k ^c：

(2)验证光源重建停止评判准则，即利用区域边界

上光通量密度的测量值Φ_k ^m与计算值Φ_k ^c之差(‖Φ_k ^m-Φ_k ^c‖或最大网格优化数L_max判断光源重建程序是否结束；

(3)若 $\left|\right|{\mathrm{\&Phi;}}_k^m-{\mathrm{\&Phi;}}_k^c\left|\right|<{\mathrm{\&epsiv;}}_{\mathrm{\&Phi;}}$ 或k≥L_max，则结束重建程序，得到重建光源S_k ^r；否则，对先验可行光源区域S_k ^p和先验禁止光源区域S_k ^*分别采用分级误差修正技术和直接最大值选择法进行后验误差估计；

(4)根据后验误差估计结果对局部网格进行细分，对局部网格进行优化分解采用“红”优化和“绿”闭合相结合的优化方法。

根据本发明的实施例，从第k层网格剖分转换至第(k+1)层网格剖分，即令k＝k+1，在新网格上重新计算边界测量数据，完成必要的参数调整，再次利用有限元理论方法将扩散方程转换成线性矩阵方程，继续进行光源迭代重建，直至满足光源重建停止评判准则。

根据本发明的实施例，所采用单谱光通量测量数据或混合谱光通量数据，用于缩短数据采集时间、提高光源重建速度。

根据本发明的实施例，所述的后验误差估计，分别对光源可行区域和光源禁止区域采取两种不同的后验误差估计方法，对光源可行区域采取直接最大选择方法，对光源禁止区域则采取分级误差修正技术。

本发明的有益效果是：利用光学特性参数与解剖结构信息的先验知识，将重建目标区域划分为先验可行光源区域和先验禁止光源区域，解决了自发荧光断层成像的病态性问题和解不唯一性问题；采用多级自适应有限元方法对重建结果进行优化，解决了过细划分网格会加重自发荧光断层成像的病态性、增加数值的不稳定和计算的代价的问题，从而使得采用本发明所述方法的自发荧光断层成像不仅提高了重建光源的位置精度及强度精度，而且还能提高重建算法的效率和鲁棒性。

附图说明

附图1为自发荧光断层成像的实验流程图。

附图2为基于本发明的自发荧光断层成像光源重建程序流程图。

附图3为用于仿真实验的非匀质生物组织仿体。其中，图3(a)表示的是带有一个光源的非均质仿体，它包括肌肉、肺、心脏、骨骼、肝脏以及光源六部分组成；图3(b)为光源重建过程中的第1层网格剖分。

具体实施方式

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

本发明的具体实施例如图1为自发荧光断层成像的实验流程图所示，具体如下：

1.利用重建目标区域的光源特性参数与解剖结构信息，将重建目标区域划分为先验可行光源区域和先验禁止光源区域。

2.在重建目标区域边界，利用高灵敏度的光学检测仪器(如液氮制冷CCD探测器)对逃逸出的光子进行单谱或混合谱探测，获得边界逸出光通量Q(r)。

3.利用扩散方程 $(-\&dtri;\&CenterDot;[D\left(r\right)\&dtri;\mathrm{\&Phi;}\left(r\right)]+{\mathrm{\&mu;}}_a\left(r\right)\mathrm{\&Phi;}\left(r\right)=S\left(r\right),r\&Element;\mathrm{\&Omega;})$ 和Robin边界条件 $(\mathrm{\&Phi;}\left(r\right)+2A(r;n,n^{/})D\left(r\right)(v\left(r\right)\&CenterDot;\&dtri;\mathrm{\&Phi;}\left(r\right))=0,r\&Element;\&PartialD;\mathrm{\&Omega;}),$ 结合多级自适应有限元算法和后验误差估计方法，对重建目标区域内部的光源进行重建。

其中，Ω和

分别是重建目标区域及其边界；Φ(r)是在点r处的光子通量密度；D(r)是与位置相关的光子扩散系数 $(D\left(r\right)=\frac{1}{3[{\mathrm{\&mu;}}_a\left(r\right)+(1-g){\mathrm{\&mu;}}_s\left(r\right)]}),{\mathrm{\&mu;}}_s\left(r\right)$ 为散射系数，g是各向异性参数；S(r)是荧光光源密度；v(r)是组织边界的外法线单位向量；A(r；n，n^/)是一个表示生物组织和周边介质折射系数不匹配的常数，它可以近似的表示为： $A(r;n,n^{/})\&ap;\frac{1+R\left(r\right)}{1-R\left(r\right)};$ n和n^/分别表示为生物组织和周边介质的折射系数，如若周边介质为空气，那么n^/≈1.0，R(r)可以近似表示为R(r)≈-1.4399n^-2+0.7099n^-1+0.6681+0.0636n。

下面对上述的步骤3进行详细说明，具体形式如下面所述：

3.1重建目标区域对光子的作用主要有散射和吸收两种，辐射传输方程是能够精确描述光子在重建目标区域中传输的数学模型，但是由于辐射传输方程是一个积分-微分方程，利用该方程进行光源重建将非常耗时。然而重建目标区域对于光子是高散射低吸收的媒质，所以扩散方程能够对光子传输进行近似模拟。当自发荧光成像实验在一个全黑的环境下进行时，光源可以认为是稳定的，于是就可以得到稳态的扩散方程：

$-\&dtri;\&CenterDot;[D\left(r\right)\&dtri;\mathrm{\&Phi;}\left(r\right)]+{\mathrm{\&mu;}}_a\left(r\right)\mathrm{\&Phi;}\left(r\right)=S\left(r\right),(r\&Element;\mathrm{\&Omega;})$

进一步考虑到重建目标区域折射系数n和外部媒介折射系数n^/的非匹配性，上述扩散方程的边界条件可以被表达为：

$\mathrm{\&Phi;}\left(r\right)+2A(r;n,n^{/})D\left(r\right)(v\left(r\right)\&CenterDot;\&dtri;\mathrm{\&Phi;}\left(r\right))=0,(r\&Element;\&PartialD;\mathrm{\&Omega;})$

3.2利用光学检测仪器探测(如液氮制冷CCD探测器)在重建目标区域边界测量得到的是光通量Q(r)，然后结合下面的公式： $Q\left(r\right)=-D\left(r\right)(v\&CenterDot;\&dtri;\mathrm{\&Phi;}\left(r\right))=\frac{\mathrm{\&Phi;}\left(r\right)}{2A(r;n,n^{\&prime;})},$ 求取边界上的光通量密度Φ(r)。

3.3自发荧光断层成像光源重建程序读入边界测量数据，初始化全局参数，设置层数为第一层，即k＝1，初始化层参数。自发荧光断层成像光源重建程序流程图如图2所示。

3.4对重建目标区域Ω进行第k层网格剖分，利用Galerkin方法可以得到Φ(r)的弱解形式：

$\underset{\mathrm{\&Omega;}}{\&Integral;}(D\left(r\right)(\&dtri;\mathrm{\&Phi;}\left(r\right))\&CenterDot;(\&dtri;\mathrm{\&Psi;}\left(r\right))+{\mathrm{\&mu;}}_a\left(r\right)\mathrm{\&Phi;}\left(r\right)\mathrm{\&Psi;}\left(r\right))\mathrm{dr}+\underset{\&PartialD;\mathrm{\&Omega;}}{\&Integral;}\frac{1}{2A_n\left(r\right)}\mathrm{\&Phi;}\left(r\right)\mathrm{\&Psi;}\left(r\right)\mathrm{dr}=\underset{\mathrm{\&Omega;}}{\&Integral;}S\left(r\right)\mathrm{\&Psi;}\left(r\right)\mathrm{dr}$

其中，Ψ(r)∈H¹(Ω)。

3.5设{J₁，J₂，…，J_k，J_k+1，…}为基于自适应有限元方法对区域Ω进行不同层网格剖分而得到的单元集合序列。其中，J_k为第k层网格剖分而得到的单元集合。假设对区域Ω进行第k层网格剖分，将会得到

个单元和

个节点，光通量密度Φ(r)和光源密度S(r)的近似表达式如下：

${\mathrm{\&Phi;}}_k\left(r\right)=\underset{i=1}{\overset{N_{P_k}}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}_i^k\left(r\right){\mathrm{\&psi;}}_i^k\left(r\right)$ $S_k\left(r\right)=\underset{i=1}{\overset{N_{S_k}}{\mathrm{\&Sigma;}}}{s}_i^k\left(r\right){\mathrm{\&gamma;}}_i^k\left(r\right)$

其中，

为节点基函数；φ_i ^k(r)为第i个节点上的光通量密度；

和s_i ^k分别为插值基函数的个数和插值点上的光源密度；

为插值基函数(k＝1，2，3…)。

3.6把光通量密度Φ(r)和光源密度S(r)的近似表达式代入到Φ(r)的弱解形式中，得到：

$\left(\underset{i=1}{\overset{N_{P_k}}{\mathrm{\&Sigma;}}}(\underset{{\mathrm{\&tau;}}_i}{\&Integral;}(D(\&dtri;{\mathrm{\&psi;}}_m^k)\&CenterDot;(\&dtri;{\mathrm{\&psi;}}_n^k)+{\mathrm{\&mu;}}_a{\mathrm{\&psi;}}_m^k\&CenterDot;{\mathrm{\&psi;}}_n^k)\mathrm{dr}+\underset{\&PartialD;{\mathrm{\&tau;}}_i}{\&Integral;}\frac{{\text{\&psi;}}_m^k\&CenterDot;{\mathrm{\&psi;}}_n^k}{2A_n}\mathrm{dr})\right)\&CenterDot;{\mathrm{\&Phi;}}_k=(\underset{i=1}{\overset{N_{\mathrm{pk}}}{\mathrm{\&Sigma;}}}\underset{{\mathrm{\&tau;}}_i}{\&Integral;}{\mathrm{\&psi;}}_m^k\&CenterDot;{\mathrm{\&gamma;}}_n^k\mathrm{dr})\&CenterDot;S_k$

其中，τ_i是第i个离散单元；

是第i个离散单元的边界。根据上式，可以得到线形矩阵方程，如下：M_kΦ_k＝F_kS_k，其中M_k为稀疏对称正定矩阵，进而可表示为： ${\mathrm{\&Phi;}}_k=M_k^{-1}{F}_k{S}_k.$

3.7删除Φ_k中的区域内部离散节点值，便得到在区域边界上测得的光通量值Φ_k ^m。

3.8利用重建目标区域光学特性参数与解剖结构信息的先验知识，S_k可划分为先验可行光源区域S_k ^p和先验禁止光源区域S_k ^*，将M_k ^-1F_k中对应于S_k ^*的列删掉而保留对应于Φ_k ^m的行。从而，我们可以得到下式： $A_k{S}_k^p={\mathrm{\&Phi;}}_k^m.$

3.9基于Tikhonov正则理论对方程 $A_k{S}_k^p={\mathrm{\&Phi;}}_k^m$ 进行求解，可以得到第k层网格剖分的目标函数Θ_k(S_k ^p)：

$\min S_{\inf }^k\&le;S_k^p\&le;S_{\sup }^k\{{\left|\right|A_k{S}_k^p-{\mathrm{\&Phi;}}_k^m\left|\right|}_{\mathrm{\&Lambda;}}+{\mathrm{\&lambda;}}_k{\mathrm{\&eta;}}_k\left(S_k^p\right)\}$

其中，S_inf ^k和S_sup ^k。分别表示第k层网格剖分的荧光光源最小和最大边界；Λ为加权矩阵；λ_k为正则化参数；η_k(S_k ^p)为补偿函数。采用带有效集策略的改进牛顿法求解上述最优化问题，得到第k层最优化结果S_k ^r。

3.10在第k层网格剖分时，如果梯度范数 $\left|\right|g_{{\mathrm{\&Theta;}}_k}^i\left|\right|=\left|\right|\&dtri;{\mathrm{\&Theta;}}_k\left(S_k^i\right)\left|\right|$ 小于阈值ε_g ^(k)或优化最后两步的间距 $\left|\right|d_{S_k}^i\left|\right|=\left|\right|S_k^i-S_k^{i-1}\left|\right|$ 小于阈值ε_d ^(k)，则计算边界光通量密度Φ_k ^c；否则继续迭代求解3.9中的最优化问题。

3.11验证光源重建停止评判准则，即利用区域边界

上光通量密度的测量值Φ_k ^m与计算值Φ_k ^c之差(‖Φ_k ^m-Φ_k ^c‖或最大网格优化数L_max判断光源重建程序是否结束。若 $\left|\right|{\mathrm{\&Phi;}}_k^m-{\mathrm{\&Phi;}}_k^c\left|\right|<{\mathrm{\&epsiv;}}_{\mathrm{\&Phi;}}$ 或k≥L_max，则结束重建程序，得到重建光源S_k ^r；否则，对先验可行光源区域S_k ^p和先验禁止光源区域S_k ^*分别采用分级误差修正技术和直接最大值选择法进行后验误差估计。

3.11根据后验误差估计结果对局部网格进行细分，网格的细分采用“红”优化和“绿”闭合相结合的优化方法。首先通过“红”优化的方法对选定的四面体单元进行细分，获得八个子四面体单元。然后通过“绿”闭合的方法对在“红”优化细分中产生的悬点，通过对邻近单元细分的方法予以处理，从而最终达到局部网格细分的目的。

3.12从第k层网格剖分转换至第(k+1)层网格剖分，即令k＝k+1，从粗糙层转换至精细层，在新网格上重新计算边界测量数据，完成必要的参数调整，转至3.4步，继续进行光源迭代重建，直至满足光源重建停止评判准则。

运行结果

为了验证本发明的方法，我们利用非匀质生物组织仿体进行了仿真实验：

第一个实验利用非匀质生物组织仿体(如图3(a)所示)对基于单谱或混合谱测量的自适应有限元光源重建方法进行了验证。该非匀质生物组织仿体是一个直径为20mm、高度为30mm的圆柱。其中，肌肉组织为圆柱形，其光学特性参数为：μ_a＝0.01mm^-1，μ_s＝4.0mm^-1，g＝0.9；两个肺组织均为椭圆形，其光学特性参数为：μ_a＝0.35mm^-1，μ_s＝23.0mm^-1，g＝0.94；心脏组织为椭圆形，其光学特性参数为：μ_a＝0.2mm^-1，μ_s＝16.0mm^-1，g＝0.85；骨骼组织为圆柱形，其光学特性参数为：μ_a＝0.002mm^-1，μ_s＝20.0mm^-1，g＝0.9；肝脏组织为椭圆形，其光学特性参数为：μ_a＝0.035mm^-1，μ_s＝6.0mm^-1，g＝0.9。根据本发明所提出的这种基于单谱或混合谱测量的自适应有限元光源重建方法，验证主要包含下列步骤：

1.在右侧肺组织中嵌入一个球体荧光光源，光源位于(-3.0，5.0，15.0)，其直径为2.0mm，密度为0.238nano-Watts/mm³。对该生物组织仿体进行网格剖分，第1层网格剖分如图3(b)所示。

2.利用分子光学仿真平台MOSE对仿体进行网格剖分，并求出仿体边界节点上的光通量密度；通过插值法求出仿体边界上其他点的光通量密度。

3.设第k层网格剖分的光源最小边界 $S_{\inf }^k=0.0$ ，最大边界S_sup ^k为一足够大的正数，加权矩阵Λ为单位矩阵，补偿函数η_k(X)＝X^TX，正则化参数λ_k＝1.16×10^-10，梯度阈值 ${\mathrm{\&epsiv;}}_g^{\left(k\right)}=1.0\&times;{10}^{-17}$ ，距离阈值 ${\mathrm{\&epsiv;}}_d^{\left(k\right)}=1.0\&times;{10}^{-16}$ ，终止阈值ε_Φ＝2.1×10^-9，最大优化网格数L_max＝4，根据不同的后验误差估计方法，在每一次局部网格细分过程中，先验可行光源区域S_k ^p中将有50％的网格被细分，而在先验禁止光源区域S_k ^*中将有1.5％的网格被细分。

4.根据从仿体前侧、右侧、后侧以及左侧视角观测仿体得到的光强分布信息，可以粗略地得到光源的位置和先验可行光源区域。利用仿体的解剖结构信息与各组织光学特性参数的先验知识以及光源深度和逸出光半高全宽值(FWHM)的关系，可以进一步指定先验可行光源区域为：

P_S＝{(x，y，z)|13.0＜z＜17.0，(x，y，z)∈Right lung}。

5.设定光源密度初始估计值 $S_1^0=1.0\&times;{10}^{-6}\mathrm{nano}-\mathrm{Watts}/{\mathrm{mm}}^3$ ，经过两步迭代， $\left|\right|{\mathrm{\&Phi;}}_k^m-{\mathrm{\&Phi;}}_k^c\left|\right|=2.05\&times;{10}^{-9},$ 满足光源重建停止评判准则，光源重建程序结束，最终得到的最大光源密度为0.219nano-Watts/mm³，相对误差为RE＝|S_recons-S_real|/S_real≈8.0％，重建光源的中心位于(-3.20，5.77，15.48)。

第二个实验将验证基于本发明的多光源重建情况。在非匀质生物组织仿体的左侧肺和右侧肺中分别嵌入两个光源，这四个光源的密度均为0.238nano-Watts/mm³，初始先验可行光源区域为：

P_S＝{(x，y，z)|11.0＜z＜19.0，(x，y，z)∈Lungs}

利用本发明所提出的这种基于单谱或混合谱测量的自适应有限元光源重建方法，对嵌入生物组织仿体中的四个光源进行重建，结果如下表：

光源序号	真实光源位置	重建光源位置	重建光源密度	相对误差
					1	(-3.0，5.0，17.0)	(-3.72，5.46，15.99)	0.209	12.1％
2	(-3.0，5.0，13.0)	(-3.64，4.41，14.25)	0.175	26.5％
					3	(-3.0，-5.0，17.0)	(-2.80，-5.98，17.05)	0.134	43.7％
4	(-3.0，-5.0，13.0)	(-3.78，-5.90，13.29)	0.127	46.6％

第三个实验用于验证基于本发明重建光源的空间分辨率。在生物组织仿体中嵌入两个球形光源，其中心的间距分别为3.5mm和3.0mm，球形光源的直径均为2.0mm，球形光源的密度均为0.238nano-Watts/mm³，初始先验可行光源区域为：

P_S＝{(x，y，z)|11.5＜z＜18.5，(x，y，z)∈Right lung}

利用本发明所提出的这种基于单谱或混合谱测量的自适应有限元光源重建方法，对嵌入生物组织仿体中的光源进行重建，结果如下表：

前面已经具体描述了本发明的实施方案，应当理解，对于一个具有本技术领域的普通技能的人，在不背离本发明的范围的情况下，在上述的和在附加的权利要求中特别提出的本发明的范围内进行变化和调整能同样达到本发明的目的。

Claims (10)

1、一种基于单谱或混合谱测量的自适应有限元光源重建方法，其特征在于：

步骤S1：利用重建目标区域的光学特性参数和解剖结构信息确定出先验可行光源区域和先验禁止光源区域；

步骤S2：根据重建目标区域边界上的单谱或混合谱光通量测量数据、结合多级自适应有限元以及后验误差估计方法，对局部网格进行优化分解，以得到重建目标区域内的荧光光源。

2、根据权利要求1所述的方法，其特征在于：将重建目标区域划分为先验可行光源区域和先验禁止光源区域，用于降低自发荧光断层成像的病态性，使其能确定出唯一解。

3、根据权利要求1所述的方法，其特征在于：利用光学检测仪器在重建目标区域边界上，测量光通量Q(r)并求取光通量密度Φ(r)。

4、根据权利要求1所述的方法，其特征在于：对重建目标区域Ω进行第k层网格剖分，利用有限元理论方法将扩散方程转换成线性矩阵方程。

5、根据权利要求4所述的方法，其特征在于：利用正则理论形成第k层网格剖分的目标函数Θ_k(S_k ^p)：

${\min }_{S_{\inf }^k\&le;S_k^p\&le;S_{\sup }^k}{\{\left|\right|A_k{S}_k^p-{\mathrm{\&Phi;}}_k^m\left|\right|}_{\mathrm{\&Lambda;}}+{\mathrm{\&lambda;}}_k{\mathrm{\&eta;}}_k\left(S_k^p\right)\}$

其中，S_inf ^k和S_sup ^k分别表示第k层网格剖分的光源最小和最大边界；Λ为加权矩阵；λ_k为正则化参数；η_k(S_k ^p)为补偿函数。

6、根据权利要求5所述的方法，其特征在于：优化目标函数Θ_k(S_k ^p)进行迭代求解，得到第k层网格剖分的光源重建结果。

7、根据权利要求6所述的方法，其特征在于：局部网格细分步骤如下：

(1)利用第k层网格剖分的光源重建结果计算光通量密度Φ_k ^c；

(2)验证光源重建停止评判准则，即利用区域边界

上光通量密度的测量值Φ_k ^m与计算值Φ_k ^c之差(‖Φ_k ^m-Φ_k ^c‖)或最大网格优化数L_max判断光源重建程序是否结束；

(3)若 $\left|\right|{\mathrm{\&Phi;}}_k^m-{\mathrm{\&Phi;}}_k^c\left|\right|<{\mathrm{\&epsiv;}}_{\mathrm{\&Phi;}}$ 或k≥L_max，则结束重建程序，得到重建光源S_k ^r；否则，对先验可行光源区域S_k ^p和先验禁止光源区域S_k ^*分别采用分级误差修正技术和直接最大选择法进行后验误差估计；

(4)根据后验误差估计结果对局部网格进行细分，对局部网格进行优化分解采用“红”优化和“绿”闭合相结合的优化方法。

8、根据权利要求6所述的方法，其特征在于：从第k层网格剖分转换至第(k+1)层网格剖分，即令k＝k+1，在新网格上重新计算边界测量数据，完成必要的参数调整，再次利用有限元理论方法将扩散方程转换成线性矩阵方程，继续进行光源迭代重建，直至满足光源重建停止评判准则。

9、根据权利要求1所述的方法，其特征在于：所采用单谱光通量测量数据或混合谱光通量数据，用于缩短数据采集时间、提高光源重建速度。

10、根据权利要求1所述的方法，其特征在于：所述的后验误差估计，分别对光源可行区域和光源禁止区域采取两种不同的后验误差估计方法，对光源可行区域采取直接最大值选择方法，对光源禁止区域则采取分级误差修正技术。

Images