CN101342075B - 基于单视图的多光谱自发荧光断层成像重建方法 - Google Patents

Abstract

基于单视图的多光谱自发荧光断层成像重建方法属于光学分子影像成像领域。该方法基于扩散方程模型，并考虑了小动物体的非均匀特性，同时也考虑自发荧光光源的光谱及实际应用的特点。为达到此目的，本发明提出的基于单视图的多光谱自发荧光断层成像重建方法的实现步骤主要包括：1)数据获取；2)有限元离散化处理；3)优化可行光源区域选择；4)光源的重建。本发明克服了自发荧光断层成像重建方法重建光源不准确，重建速度比较慢，不利于实时处理以及在多光谱数据采集时可能存在误差的缺陷。

Description

基于单视图的多光谱自发荧光断层成像重建方法

技术领域

本发明涉及到光学分子影像成像模态-自发荧光断层成像(BLT)技术，尤其涉及到一种基于单视图的多光谱自发荧光断层成像重建方法。

背景技术

在过去的几年里，分子影像由于能在体揭示分子和细胞的信息而受到了越来越多的关注。因此它已成为一种疾病诊断、药物疗效评价等的重要工具。作为一种小动物分子影像成像模态，光学成像技术特别是自发荧光断层成像因为其高性能、低价格和无创伤等特性已经得到了迅速发展和广泛的应用。自发荧光断层成像是最近刚发展起来的，是一种用来在活体小动物体内重建自发荧光光源分布的光学成像技术。当荧光素与荧光素酶在有氧气和三磷酸腺苷(ATP)的条件下，就会发出荧光。而因为荧光素酶含有萤火虫素酶(firefly)，磕头虫素酶(click beetle)，所以发出的荧光有不同的光谱，波长一般在400nm-750nm。发出的荧光穿透生物体而到达生物体表面，然后用高灵敏度的液氮制冷电荷耦合器件(CCD)获得。CCD获得的动物表面的数据形成了BLT重建的基础。为采集整个生物体的表面数据，通常将生物体每隔90度旋转一次，用CCD探测器来获取生物体的前后左右四个视图。当光源在生物体的位置比较深时，采集一个视图的数据需要的时间大约是5-10分钟。而当采集时间超过一个小时之后荧光强度就会衰减，所以在多谱的情况下，如果在每个单谱段都采集四个视图的数据，那么采集时间远远超过1个小时，可能就会使得采集的数据不够准确。另外，在实际情况下存在类似像X射线(x-ray)，计算机断层成像(CT)一样的物理限制，采集数据时就受到角度限制，因此采集到的数据仅是动物体表面的一部分。多谱采集也会使得系统矩阵的维数增加，导致计算量过大。如何降低测量时间和减小系统矩阵的维数是目前的一个难点。另一方面，在进行光源重建时，因为BLT自身具有的不适定特点导致重建结果不准确。

BLT的目标是能在体实时连续的观测。因此，发展快速的重建方法就成为亟待解决的问题。目前的重建方法常用的是有限元方法。理论上有限元网格越细，得到的结果越好，但重建的时间就越长。此外，BLT是病态问题，如何降低病态性仍需要进一步的探索。

发明内容

本发明的目的在于克服了自发荧光断层成像重建方法重建光源不准确，重建速度比较慢，不利于实时处理以及在多光谱数据采集时可能存在误差的缺陷，提出了一种基于单视图的多光谱自发荧光断层成像重建方法，该方法基于扩散方程模型，并考虑了小动物体的非均匀特性，同时也考虑自发荧光光源的光谱及实际应用的特点。

为达到此目的，本发明提出的基于单视图的多光谱自发荧光断层成像重建方法的实现步骤主要包括：1)数据获取；2)有限元离散化处理；3)优化可行光源区域选择；4)光源的重建。

1)数据获取是指利用CCD探测器在被测物体表面对逃出的光子进行探测，获得边界处的光流密度；

有限元离散化处理，是指利用有限元方法将扩散方程离散为矩阵方程，并把重建区域网格剖分为{T₁，...，T_k，...}的网格序列；

优化可行光源区域的选择，主要是在对扩散方程离散化处理的基础上，利用迭代方法快速选择光源大致存在的区域，即优化可行光源区域；

光源的重建，主要是在优化可行光源区域选择的基础上，确立未知光源密度变量与被测物体表面测量值之间的线性关系，并利用Tikhonov正则化方法建立目标函数，并最终求解光源密度分布。

本发明具体包括以下步骤：

(1)数据获取：

在多光谱下，用滤波片将荧光波长λ离散为m个波段τ₁，...τ_m，其中τ_l＝[λ_l-1，λ_l]，l＝1，2，...m-1，τ_m＝[λ_m-1，λ_m]，在多光谱情况下m＝5即可提供足够的先验信息；利用CCD探测器在被测物体的表面对逃出的光子进行探测，获得每个波段τ_l的光流密度Q(x，τ_l)，其中CCD探测到的表面为γ，x表示在被测物体中的位置；根据下面公式计算被测物体表面的光子流密度Φ_meas(x，τ_l)：

$Q(x,{\mathrm{\τ}}_l)=-D(x,{\mathrm{\τ}}_l)(v\left(x\right)\·\▿{\mathrm{\Φ}}^{\mathrm{meas}}(x,{\mathrm{\τ}}_l))=\frac{{\mathrm{\Φ}}^{\mathrm{meas}}(x,{\mathrm{\τ}}_l)}{2A(x;n,n^{\′})}(x\∈\mathrm{\γ})$

D(x，τ_l)＝1/(3(μ_a(x，τ_l)+(1-g)μ_s(x，τ_l)))是生物组织的扩散系数，其中μ_a(x，τ_l)是生物组织的吸收系数，μ_s(x，τ_l)是生物组织散射系数，g是生物组织各向异性参数；v(x)是被测物体的表面

的单位法向量；A(x；n，n′)≈(1+R(x))/(1-R(x))，n′为外界媒质的折射系数；当外界媒质是空气时，R(x)可以近似为R(x)≈-1.4399n^-2+0.7099n^-1+0.6681+0.0636n，其中n为空气的折射率；

(2)有限元离散化处理：

光子在生物组织中传输的数学模型用下面的扩散方程表示：

$-\▿\·(D(x,{\mathrm{\τ}}_l)\▿\mathrm{\Φ}(x,{\mathrm{\τ}}_l))+{\mathrm{\μ}}_a(x,{\mathrm{\τ}}_l)\mathrm{\Φ}(x,{\mathrm{\τ}}_l)=S(x,{\mathrm{\τ}}_l)(x\∈\mathrm{\Ω})$

$\mathrm{\Φ}(x,{\mathrm{\τ}}_l)+2A(x;n,n^{\′})D(x,{\mathrm{\τ}}_l)(v\left(x\right)\·\▿\mathrm{\Φ}(x,{\mathrm{\τ}}_l))=0(x\∈\∂\mathrm{\Ω})$

其中Ω是被测物体；

是被测物体的表面；Φ(x，τ_l)是光子流密度；S(x，τ_l)是光源密度；根据有限元理论，得到下面的弱解形式：

${\∫}_{\mathrm{\Ω}}(D(x,{\mathrm{\τ}}_l)(\▿\mathrm{\Φ}(x,{\mathrm{\τ}}_l))\·(\▿\mathrm{\Ψ}(x,{\mathrm{\τ}}_l))+{\mathrm{\μ}}_a(x,{\mathrm{\τ}}_l)\mathrm{\Φ}(x,{\mathrm{\τ}}_l)\mathrm{\Ψ}(x,{\mathrm{\τ}}_l))\mathrm{dx}$

$+{\∫}_{\∂\mathrm{\Ω}}\frac{1}{2A_n\left(x\right)}\mathrm{\Φ}(x,{\mathrm{\τ}}_l)\mathrm{\Ψ}(x,{\mathrm{\τ}}_l)\mathrm{dx}={\∫}_{\mathrm{\Ω}}S(x,{\mathrm{\τ}}_l)\mathrm{\Ψ}(x,{\mathrm{\τ}}_l)\mathrm{dx}(\∀\mathrm{\Ψ}(x,{\mathrm{\τ}}_l)\∈H^1\left(\mathrm{\Ω}\right))$

H¹(Ω)是Sobelev空间，Ψ(x，τ_l)是任意给定的试探函数。在自适应有限元分析的框架下，基于自适应的网格细分，把重建区域网格剖分表示为{T₁，...，T_k，...}的网格序列，其中网格剖分T_k包括

个离散单元。利用有限元方法，对弱解形式进行离散，得到下面矩阵方程：

(K_k(τ_l)+C_k(τ_l)+B_k(τ_l))Φ_k(τ_l)＝F_k(τ_l)S_k(τ_l)

矩阵元素的具体形式为：

$\left\{\begin{array}{c}{k}_{\mathrm{ij}}^k\left({\mathrm{\τ}}_l\right)={\∫}_{\mathrm{\Ω}}D(x,{\mathrm{\τ}}_l)(\▿{\mathrm{\Ψ}}_i^k(x,{\mathrm{\τ}}_l))\·(\▿{\mathrm{\Ψ}}_j^k(x,{\mathrm{\τ}}_l))\mathrm{dx}\\ c_{\mathrm{ij}}^k\left({\mathrm{\τ}}_l\right)={\∫}_{\mathrm{\Ω}}{\mathrm{\μ}}_a(x,{\mathrm{\τ}}_l){\mathrm{\Ψ}}_i^k(x,{\mathrm{\τ}}_l){\mathrm{\Ψ}}_j^k(x,{\mathrm{\τ}}_l)\mathrm{dx}\\ b_{\mathrm{ij}}^k\left({\mathrm{\τ}}_l\right)={\∫}_{\∂\mathrm{\Ω}}{\mathrm{\Ψ}}_i^k(x,{\mathrm{\τ}}_l){\mathrm{\Ψ}}_j^k(x,{\mathrm{\τ}}_l)/\left(2A(x;n,n^{\′})\right)\mathrm{dx}\\ f_{\mathrm{ij}}^k\left({\mathrm{\τ}}_l\right)={\∫}_{\mathrm{\Ω}}{\mathrm{\Ψ}}_i^k(x,{\mathrm{\τ}}_l){\mathrm{\Ψ}}_j^k(x,{\mathrm{\τ}}_l)\mathrm{dx}\end{array}\right.$

令K_k(τ_l)+C_k(τ_l)+B_k(τ_l)＝M_k(τ_l)，M_k(τ_l)是稀疏正定矩阵，所以得到：

${\mathrm{\Φ}}_k\left({\mathrm{\τ}}_l\right)=M_k^{-1}\left({\mathrm{\τ}}_l\right)F_k\left({\mathrm{\τ}}_l\right)S_k\left({\mathrm{\τ}}_l\right)---\left(1\right)$

通过删除[M_k ^-1(τ_l)F_k(τ_l)]中对应非测量点的行得到下面的方程：

Φ_k(τ_l)＝A_k(τ_l)S_k(τ_l)                  (2)

因此对于m个波段即得到m个单一波段的矩阵方程；

(3)优化可行光源区域选择

I)在网格剖分T₁上，将上述含有m个单一波段的矩阵方程组合成含有m个波段的一个方程：

根据光源谱信息的性质，每个谱段τ_l上所占整个谱段的能量表示为S(τ_l)＝ω(τ_l)·S，其中ω(τ_l)≥0，

ω(τ_l)可通过实验测定；S表示光源总的能量。在网格剖分T₁上，考虑光源谱信息的性质，通过方程(2)得到：

AS₁＝Φ                        (3)

其中 $\mathrm{\Φ}=\left[\begin{array}{c}{\mathrm{\Φ}}_1\left({\mathrm{\τ}}_1\right)\\ {\mathrm{\Φ}}_1\left({\mathrm{\τ}}_2\right)\\ .\\ .\\ .\\ {\mathrm{\Φ}}_1\left({\mathrm{\τ}}_m\right)\end{array}\right],$ $A=\left[\begin{array}{c}\mathrm{\ω}\left({\mathrm{\τ}}_1\right)A_1\left({\mathrm{\τ}}_1\right)\\ \mathrm{\ω}\left({\mathrm{\τ}}_2\right)A_1\left({\mathrm{\τ}}_2\right)\\ .\\ .\\ .\\ \mathrm{\ω}\left({\mathrm{\τ}}_m\right)A_1\left({\mathrm{\τ}}_m\right)\end{array}\right]$

Φ₁(τ_l)和A₁(τ_l)都由方程(2)计算得到，S₁是网格剖分T₁上的未知光源密度变量；在(3)式两边同乘A^T，得到：

A^TAS₁＝A^TΦ                    (4)

A^T是A的转置矩阵，A^T·A是一个的对称矩阵，所以(4)式是一个标准线形方程。

II)在网格剖分T₁上，确定优化光源可行区域：

具体方法为：在网格剖分T₁上，假定整个被测物体Ω为光源可行区域，给定未知光源密度的初始值S₁ ⁰，根据公式(4)，网格剖分T₁上的未知光源密度变量S₁用下面的关系式进行迭代求解：

$S_1^{n+1}=S_1^n+{\mathrm{\α}}_n{\mathrm{\β}}_n$

其中步长

搜索方向

n是迭代次数；

在迭代的过程中，用Φ₁ ^meas(τ_l)代替Φ₁(τ_l)，这是因为噪声在迭代的过程中产生的影响很小；Φ₁ ^meas(τ_l)是网格剖分T₁上的边界测量值，由Φ^meas(x，τ_l)通过最近邻域插值得到；当||β_n||≤δ或n＞N_max，迭代停止，此时便得到了未知光源密度S₁的优化解S^*，优化解S^*的大致区域为Ω^*，这个区域Ω^*我们称之为优化的光源可行区域；在每次迭代的过程中，如果光源的密度小于零，则将其置为零；δ取值范围为10^-5-10^-2，N_max为最大迭代次数，取值不超过1000。

(4)光源重建

在网格序列{T₁，...，T_k，...}上，根据选择好的优化光源可行区域Ω^*，重新将含有m个单一波段的矩阵方程组合成含有m个波段的一个方程：

考虑未知光源密度S_k(τ_l)与边界测量值Φ_k ^meas(τ_l)之间的线性关系，由方程(1)： ${\mathrm{\Φ}}_k\left({\mathrm{\τ}}_l\right)=M_k^{-1}\left({\mathrm{\τ}}_l\right)F_k\left({\mathrm{\τ}}_l\right)S_k\left({\mathrm{\τ}}_l\right),$ 得到：

${\mathrm{\Φ}}_k^{\mathrm{meas}}\left({\mathrm{\τ}}_l\right)=G_k\left({\mathrm{\τ}}_l\right)S_k\left({\mathrm{\τ}}_l\right)$

其中G_k(τ_l)通过移出[M_k(τ_l)^-1F_k(τ_l)]中对应非测量点的行而得到；考虑到光源谱信息和优化可行光源区域Ω^*，得到：

${\mathrm{\Φ}}_k^{\mathrm{meas}}=G_k{W}_k{S}_k$

其中 $G_k=\left[\begin{array}{c}w\left({\mathrm{\&tau;}}_1\right)G_k\left({\mathrm{\&tau;}}_1\right)\\ w\left({\mathrm{\&tau;}}_2\right)G_k\left({\mathrm{\&tau;}}_2\right)\\ .\\ .\\ .\\ w\left({\mathrm{\&tau;}}_m\right)G_k\left({\mathrm{\&tau;}}_m\right)\end{array}\right],$ ${\mathrm{\&Phi;}}_k=\left[\begin{array}{c}{\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_1\right)\\ {\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_2\right)\\ .\\ .\\ .\\ {\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_m\right)\end{array}\right],$ $\begin{array}{c}{W}_k=\mathrm{Diag}(w_{k\left(11\right)},w_{k\left(22\right)},...,w_{k\left(\mathrm{ii}\right)},...),\\ w_{k\left(\mathrm{ii}\right)}=\left\{\begin{array}{cc}1& s_{k\left(i\right)}\&GreaterEqual;{\mathrm{\&gamma;}}^k{S}_k^{\max }\\ 0& s_{k\left(i\right)}<{\mathrm{\&gamma;}}^k{s}_k^{\max }\end{array}\right.\end{array}$

G_k表示多谱下的边界测量值与光源位置量之间的关系，W_k为k层网格上的对角矩阵，w_k(ii)是对角矩阵W_k的第i个对角元素；Φ_k ^meas(τ_l)是波段τ_l在网格剖分T_k上的边界测量值，由Φ^meas(x，τ_l)通过最近邻域插值得到；当k＝1时，s_k(i)是通过迭代方法得到的优化解S^*，s_k ^max是优化解S^*的最大值；当k≥2时，s_k(i)和s_k ^max是网格剖分T_k上的未知光源密度的初始值和初始值的最大值，通过插值网格剖分T_k-1上的重建结果获得，即

S_k-1是在网格剖分T_k-1上重建得到的结果，I_k-1 ^k是插值因子，通过子单元继承父单元的重建结果值的方式来实现；γ_k是比例因子，是随k变化的分段常数；最终通过保留G_kW_k中对应可行光源区域的列，得到光源可行区域的未知光源密度变量S_k ^p与表面测量值Φ_k ^meas在网格剖分T_k上的关系： $A_k{S}_k^p={\mathrm{\&Phi;}}_k^{\mathrm{meas}}$

II)利用Tikhonov正则化方法建立目标函数并优化从而得到重建结果：

上一步虽然得到了未知光源密度变量S_k ^p与表面测量值Φ_k ^meas之间的关系，由于A_k矩阵的病态，不能通过直接求解的方法得到重建结果，因此利用Tikhonov正则的方法，并考虑未知光源密度变量不能为负值的特性，得到在网格剖分T_k上的目标函数：

$\underset{0\&le;s_k^p\&le;s_{\sup }^k}{\min }{\mathrm{\&Theta;}}_k\left(S_k^p\right)=\{{\left|\right|A_k{S}_k^p-{\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left|\right|}_{\mathrm{\&Lambda;}}+{\mathrm{\&lambda;}}_k\&CenterDot;{(S_k^p-S_k^{\mathrm{init}})}^T(S_k^p-S_k^{\mathrm{init}})\}$

考虑到多光谱以及非匀质模型会极大的增加计算量，因此采用谱梯度大规模优化算法对该目标函数进行优化；S_sup ^k是光源密度上界，经验取值为不大于1000，||V||_Λ＝V^TΛV，λ_k正则化系数，经验取值10^-9-10^-6，S_k ^init是网格剖分T_k上的未知光源密度初始值，取值范围为10^-7-10^-4；在判断优化过程是否中止时，采用当前优化梯度

以及优化迭代次数N_k ^te作为评判准则，即当

或

时，优化停止，得到重建结果，即光源可行区域的光源密度分布；上标t_e表示当前的迭代次数；ε_g是梯度阈值，其经验取值不超过10^-6；N_k ^max是最大优化迭代次数，经验取值不超过10⁴；

III)当优化完成后，判断重建是否停止：

利用优化求得的光源密度求解任意一个谱段τ_l上的光密度分布Φ_k ^c(τ_l)，并计算

当

或k≥K_max，重建停止；Φ_k ^c(τ_l)是通过将优化求得的光源密度代入扩散方程而求得的边界处光子流密度，Φ_k ^meas(τ_l)是波段τ_l在网格剖分T_k上边界测量点的测量值，由Φ^meas(x，τ_l)通过最近邻域插值得到；ε_Φ经验取值范围10^-9-10^-7；K_max是网格最大剖分次数；如果不能满足

或k≥K_max，根据重建的结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域；通过以上方法对可行和禁止光源区域内的网格单元的选择，对区域剖分网格进行自适应的细分；并从网格剖分T_k转到T_k+1，然后转到第(4)(I)步，重复前面的步骤，直至重建停止；

本发明基于单视图测量数据，利用多光谱提供的先验信息，通过迭代方法进行可行光源区域的选择，重建光源的不准确问题得到了很好的解决，重建的速度和质量由于光源可行区域的使用而被进一步提高，此方法对自发荧光断层成像的发展具有重要的应用价值。

附图说明：

图1.非均匀仿体

图中：(a)非均匀仿体，1代表肌肉、2骨骼、3心脏、4肺、5肝和6单光源；(b)重建算法用到的初始网格；(c)包含肌肉、骨骼、肺的横截面，光源在右肺中；箭头所示的方向为CCD探测的方向，7前视图方向，8左视图方向，9后视图方向，10右视图方向。

图2.程序主流程图。

图3.确立光源可行区域的子流程图。

图4.确定光源可行区域的算法示意图。

图5.重建结果示意图。

表1.各组织在不同波长下的光学特性参数。

具体实施方式：

下面结合附图详细描述本实施例。

(1)数据获取：

为了验证本方法，我们设计了非匀质模型来近似模仿小鼠腹部内的组织，如图1(a)所示，生物组织分别为：肌肉、骨骼、心脏、肺、肝，各个组织的光学特性参数如表1所示。在此实验中，我们将重建自发光源的光源密度分布作为重建目标，在重建实验中，这个光源大小为半径1毫米，光源密度为0.238nano-Watts/mm³，位置为(-3，5，15)。它的波谱范围从500纳米到750纳米，为了进行多光谱实验，我们依据当前的探测水平，把这一波谱范围分成了5个谱段进行分别探测，分别为τ₁＝[500，550)nm，τ₂＝[550，600)nm，τ₃＝[600，650)nm，τ₄＝[650，700)nm，τ₅＝[700，750)nm。边界测量数据由蒙特卡罗方法产生，并按照图1(c)所示的方法取前视图的边界数据，并求得每个波段的光子流密度Φ^meas(x，τ_l)；

(2)有限元离散化处理：

光子在生物组织中传输的数学模型用下面的扩散方程表示：

$-\&dtri;\&CenterDot;(D(x,{\mathrm{\&tau;}}_l)\&dtri;\mathrm{\&Phi;}(x,{\mathrm{\&tau;}}_l))+{\mathrm{\&mu;}}_a(x,{\mathrm{\&tau;}}_l)\mathrm{\&Phi;}(x,{\mathrm{\&tau;}}_l)=S(x,{\mathrm{\&tau;}}_l)(x\&Element;\mathrm{\&Omega;})$

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

其中Ω是被测物体；

是被测物体的表面；Φ(x，τ_l)是光子流密度；S(x，τ_l)是光源密度；根据有限元理论，得到下面的弱解形式：

${\&Integral;}_{\mathrm{\&Omega;}}(D(x,{\mathrm{\&tau;}}_l)(\&dtri;\mathrm{\&Phi;}(x,{\mathrm{\&tau;}}_l))\&CenterDot;(\&dtri;\mathrm{\&Psi;}(x,{\mathrm{\&tau;}}_l))+{\mathrm{\&mu;}}_a(x,{\mathrm{\&tau;}}_l)\mathrm{\&Phi;}(x,{\mathrm{\&tau;}}_l)\mathrm{\&Psi;}(x,{\mathrm{\&tau;}}_l))\mathrm{dx}$

$+{\&Integral;}_{\&PartialD;\mathrm{\&Omega;}}\frac{1}{2A_n\left(x\right)}\mathrm{\&Phi;}(x,{\mathrm{\&tau;}}_l)\mathrm{\&Psi;}(x,{\mathrm{\&tau;}}_l)\mathrm{dx}={\&Integral;}_{\mathrm{\&Omega;}}S(x,{\mathrm{\&tau;}}_l)\mathrm{\&Psi;}(x,{\mathrm{\&tau;}}_l)\mathrm{dx}(\&ForAll;\mathrm{\&Psi;}(x,{\mathrm{\&tau;}}_l)\&Element;H^1\left(\mathrm{\&Omega;}\right))$

H¹(Ω)是Sobelev空间，Ψ(x，τ_l)是任意的试探函数。在自适应有限元分析的框架下，基于自适应的网格细分，把重建区域网格剖分表示为{T₁，...，T_k，...}的网格序列，其中网格剖分T_k包括

个离散单元，网格剖分T₁包含6878个离散单元，如图1(b)所示。利用有限元方法，对弱解形式进行离散，在每个波段τ_l(l＝1，2，3，4，5)上得到下面矩阵方程：

(K_k(τ_l)+C_k(τ_l)+B_k(τ_l))Φ_k(τ_l)＝F_k(τ_l)S_k(τ_l)

矩阵元素的具体形式为：

$\left\{\begin{array}{c}{k}_{\mathrm{ij}}^k\left({\mathrm{\&tau;}}_l\right)={\&Integral;}_{\mathrm{\&Omega;}}D(x,{\mathrm{\&tau;}}_l)(\&dtri;{\mathrm{\&Psi;}}_i^k(x,{\mathrm{\&tau;}}_l))\&CenterDot;(\&dtri;{\mathrm{\&Psi;}}_j^k(x,{\mathrm{\&tau;}}_l))\mathrm{dx}\\ c_{\mathrm{ij}}^k\left({\mathrm{\&tau;}}_l\right)={\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&mu;}}_a(x,{\mathrm{\&tau;}}_l){\mathrm{\&Psi;}}_i^k(x,{\mathrm{\&tau;}}_l){\mathrm{\&Psi;}}_j^k(x,{\mathrm{\&tau;}}_l)\mathrm{dx}\\ b_{\mathrm{ij}}^k\left({\mathrm{\&tau;}}_l\right)={\&Integral;}_{\&PartialD;\mathrm{\&Omega;}}{\mathrm{\&Psi;}}_i^k(x,{\mathrm{\&tau;}}_l){\mathrm{\&Psi;}}_j^k(x,{\mathrm{\&tau;}}_l)/\left(2A(x;n,n^{\&prime;})\right)\mathrm{dx}\\ f_{\mathrm{ij}}^k\left({\mathrm{\&tau;}}_l\right)={\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&Psi;}}_i^k(x,{\mathrm{\&tau;}}_l){\mathrm{\&Psi;}}_j^k(x,{\mathrm{\&tau;}}_l)\mathrm{dx}\end{array}\right.$

令K_k(τ_l)+C_k(τ_l)+B_k(τ_l)＝M_k(τ_l)，M_k(τ_l)是稀疏正定矩阵，所以得到：

${\mathrm{\&Phi;}}_k\left({\mathrm{\&tau;}}_l\right)=M_k^{-1}\left({\mathrm{\&tau;}}_l\right)F_k\left({\mathrm{\&tau;}}_l\right)S_k\left({\mathrm{\&tau;}}_l\right)$ (1)通过删除

[M_k ^-1(τ_l)F_k(τ_l)]中对应非测量点的行得到下面的方程：

Φ_k(τ_l)＝A_k(τ_l)S_k(τ_l)                     (2)

对于5个波段得到5个单一波段的矩阵方程；

(3)优化可行光源区域选择

I)在网格剖分T₁上，将上述含有5个单一波段的矩阵方程组合成含有5个波段的一个方程：

根据光源谱信息的性质，每个谱段τ_l上所占整个谱段的能量表示为

S(τ_l)＝ω(τ_l)·S，ω(τ_l)≥0，

其中S表示光源总的光源能量，值为1；ω(τ_l)通过实验测定分别为：ω(τ₁)＝0.222，ω(τ₂)＝0.167，ω(τ₃)＝0.222，ω(τ₄)＝0.167，ω(τ₅)＝0.222。

在网格剖分T₁上，考虑光源谱信息的性质，通过方程(2)得到：

AS₁＝Φ                                     (3)

其中 $\mathrm{\&Phi;}=\left[\begin{array}{c}{\mathrm{\&Phi;}}_1\left({\mathrm{\&tau;}}_1\right)\\ {\mathrm{\&Phi;}}_1\left({\mathrm{\&tau;}}_2\right)\\ {\mathrm{\&Phi;}}_1\left({\mathrm{\&tau;}}_3\right)\\ {\mathrm{\&Phi;}}_1\left({\mathrm{\&tau;}}_4\right)\\ {\mathrm{\&Phi;}}_1\left({\mathrm{\&tau;}}_5\right)\end{array}\right],$ $A=\left[\begin{array}{c}0.222A_1\left({\mathrm{\&tau;}}_1\right)\\ 0.167A_1\left({\mathrm{\&tau;}}_2\right)\\ 0.222A_1\left({\mathrm{\&tau;}}_3\right)\\ 0.167A_1\left({\mathrm{\&tau;}}_4\right)\\ 0.222A_1\left({\mathrm{\&tau;}}_5\right)\end{array}\right]$

Φ₁(τ_l)和A₁(τ_l)都由方程(2)计算得到，S₁是网格剖分T₁上的未知光源密度变量。在(3)式两边同乘A^T，得到：

A^TAS₁＝A^TΦ                                  (4)

A^T是A的转置矩阵，A^T·A是一个34390×34390的对称矩阵，所以(4)式是一个标准线形方程。

III)在网格剖分T₁上，确定优化光源可行区域：

具体方法为：在网格剖分T₁上，假定整个被测物体Ω为光源可行区域，给定未知光源密度的初始值

根据公式(4)，网格剖分T₁上的未知光源密度变量S₁根据下面的关系式进行迭代求解：

$S_1^{n+1}=S_1^n+{\mathrm{\&alpha;}}_n{\mathrm{\&beta;}}_n$

其中步长

搜索方向

n是迭代次数。

在迭代的过程中，用Φ₁ ^meas(τ_l)代替Φ₁(τ_l)，Φ₁ ^meas(τ_l)是τ_l在网格剖分T₁上的边界测量值，由Φ^meas(x，τ_l)通过最近邻域插值得到；参数取值分别为：δ＝5×10^-3，N_max＝500。经过5次迭代后，||β_n||＝4.9e-5，||β_n||≤δ，所以迭代停止，此时便得到了未知光源密度S₁的优化解S^*，优化解S^*的大致区域为Ω^*，这个区域Ω^*我们称之为优化的光源可行区域。在每次迭代的过程中，如果光源的密度小于零，则将其置为零。

(4)光源重建

(I)在被测物体Ω的网格剖分{T₁，...，T_k，...}的网格序列上，根据选择好的优化光源可行区域Ω^*，重新将含有5个波段的矩阵方程组合成含有5个波段的一个方程：

考虑未知光源密度S_k(τ_l)与边界测量值Φ_k ^meas(τ_l)之间的线性关系，由方程(1)： ${\mathrm{\&Phi;}}_k\left({\mathrm{\&tau;}}_l\right)=M_k^{-1}\left({\mathrm{\&tau;}}_l\right)F_k\left({\mathrm{\&tau;}}_l\right)S_k\left({\mathrm{\&tau;}}_l\right),$ 得到：

${\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_l\right)=G_k\left({\mathrm{\&tau;}}_l\right)S_k\left({\mathrm{\&tau;}}_l\right)$

其中G_k(τ_l)通过移出[M_k(τ_l)^-1F_k(τ_l)]中对应非测量点的行而得到。考虑到光源谱信息和优化可行光源区域Ω^*，得到：

${\mathrm{\&Phi;}}_k^{\mathrm{meas}}=G_k{W}_k{S}_k$

其中 $G_k=\left[\begin{array}{c}w\left({\mathrm{\&tau;}}_1\right)G_k\left({\mathrm{\&tau;}}_1\right)\\ w\left({\mathrm{\&tau;}}_2\right)G_k\left({\mathrm{\&tau;}}_2\right)\\ w\left({\mathrm{\&tau;}}_3\right)G_k\left({\mathrm{\&tau;}}_3\right)\\ w\left({\mathrm{\&tau;}}_4\right)G_k\left({\mathrm{\&tau;}}_4\right)\\ w\left({\mathrm{\&tau;}}_5\right)G_k\left({\mathrm{\&tau;}}_5\right)\end{array}\right],$ ${\mathrm{\&Phi;}}_k=\left[\begin{array}{c}{\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_1\right)\\ {\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_2\right)\\ {\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_3\right)\\ {\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_4\right)\\ {\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_5\right)\end{array}\right],$ $\begin{array}{c}{W}_k=\mathrm{Diag}(w_{k\left(11\right)},w_{k\left(22\right)},...,w_{k\left(\mathrm{ii}\right)},...),\\ w_{k\left(\mathrm{ii}\right)}=\left\{\begin{array}{cc}1& s_{k\left(i\right)}\&GreaterEqual;{\mathrm{\&gamma;}}^k{S}_k^{\max }\\ 0& s_{k\left(i\right)}<{\mathrm{\&gamma;}}^k{s}_k^{\max }\end{array}\right.\end{array}$

G_k表示多谱下的边界测量值与光源位置量之间的关系，W_k为k层网格上的对角矩阵，Φ_k ^meas(τ_l)是波段τ_l在网格剖分T_k上的边界测量值。当k＝1时，s_k(i)是通过迭代方法得到的优化解S^*，s_k ^max是优化解S^*的最大值；当k≥2时，s_k(i)和s_k ^max是网格剖分T_k上的初始值和初始值的最大值，通过插值网格剖分T_k-1上的重建结果得到，即

S_k-1是在网格剖分T_k-1上重建得到的结果，I_k-1 ^k是插值因子，通过子单元继承父单元的重建结果值的方式来实现；γ^k是比例因子，是随k变化的分段常数；k＝1时，参数γ^k＝0；k＞1时，γ^k初始为0.1，并随着k增加，有γ^k+1＝10γ^k。最终通过保留G_kW_k中对应可行光源区域的列，得到光源可行区域的未知光源密度变量S_k ^p与表面测量值Φ_k ^meas在网格剖分T_k上的关系：

$A_k{S}_k^p={\mathrm{\&Phi;}}_k^{\mathrm{meas}}$

(II)利用Tikhonov正则化方法建立目标函数并优化从而得到重建结果：

上一步虽然得到了未知光源密度变量S_k ^p与表面测量值Φ_k ^meas之间的关系，由于A_k矩阵的病态，不能通过直接求解的方法得到重建结果，利用Tikhonov正则的方法，并考虑未知光源密度变量不能为负值的特性，得到在网格剖分T_k上的目标函数：

$\underset{0\&le;s_k^p\&le;s_{\sup }^k}{\min }{\mathrm{\&Theta;}}_k\left(S_k^p\right)=\{{\left|\right|A_k{S}_k^p-{\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left|\right|}_{\mathrm{\&Lambda;}}+{\mathrm{\&lambda;}}_k\&CenterDot;{(S_k^p-S_k^{\mathrm{init}})}^T(S_k^p-S_k^{\mathrm{init}})\}$

考虑到多光谱以及非匀质模型会极大的增加计算量，因此采用谱梯度大规模优化算法对该目标函数进行优化。S_sup ^k是光源密度上界，取值为100，||V||_Λ＝V^TΛV，λ_k正则化系数，取值1.0×10^-8，S_k ^init是在网格剖分T_k上未知光源密度变量的初始值，取值1.0×10^-5。在判断优化过程是否中止时，采用了当前优化梯度

以及优化迭代次数N_k ⁱ作为评判准则，即当

或时，优化停止，得到重建结果，即光源可行区域的光源密度分布，其中ε_g是梯度阈值，取值1.0×10^-8，N_k ^max是最大优化迭代次数，取值1.0×10⁴。

(III)当优化完成后，判断重建是否停止：

利用优化得到的光源密度分布求解任意一个谱段τ_l上的光子流密度Φ_k ^c(τ_l)，并计算

当

(ε_Φ为正数)或k≥K_max(K_max是网格最大剖分次数)，重建停止。Φ_k ^c(τ_l)是通过将优化求得的光源密度代入扩散方程而求得的边界处光子流密度，Φ_k ^meas(τ_l)是波段τ_l在第k层网格上边界测量点的测量值；如果不能满足

或k≥K_max，利用重建的结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域。通过以上方法对可行和禁止光源区域内的网格单元的选择，并对区域剖分网格进行自适应的细分，从网格剖分T_k转到T_k+1，然后转到第(4)(I)步，重复前面的步骤，直至重建停止。

停止阈值ε_Φ和网格最大细分次数K_max分别是1×10^-8和3。当在网格剖分T₁上重建完成后，计算在τ₁＝[500-550)nm的光密度分布Φ_k ^c(τ_l)，求得

此时

(k＝1)＜(K_max＝3)所以重建未停止；利用网格剖分T₁上的重建结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域。通过以上方法对可行和禁止光源区域内的网格单元的选择，并对区域剖分网格进行自适应的细分，并从网格剖分T₁转到网格剖分T₂，然后转到第(4)(I)步，重复前面的步骤；优化完成后计算在τ₁＝[500-550)nm的光密度分布Φ_k ^c(τ_l)，求得

因此重建停止，重建的结果如图5所示，重建的位置为(-1.53，5.3，14.98)，此时重建的最大光源密度为0.224nano-Watts/mm³。

表1

Claims (1)

1.一种基于单视图的多光谱自发荧光断层成像重建方法，其特征在于，包括以下步骤：

1)数据获取：

在多光谱下，用滤波片将荧光波长λ离散为m个波段τ₁，...，τ_m，其中τ_l＝[λ_l-1，λ_l]，l＝1，2，...，m-1，τ_m＝[λ_m-1，λ_m]；利用CCD探测器在被测物体的表面

对逃出的光子进行探测，获得每个波段τ_l的光流密度Q(x，τ_l)，其中CCD探测到的表面为γ，x表示在被测物体中的位置；根据下面公式计算被测物体表面的光子流密度Φ^meas(x，τ_l)：

D(x，τ_l)＝1/(3(μ_a(x，τ_l)+(1-g)μ_s(x，τ_l)))是生物组织的扩散系数，其中μ_a(x，τ_l)是生物组织的吸收系数，μ_s(x，τ_l)是生物组织散射系数，g是生物组织各向异性参数；v(x)是被测物体的表面

的单位法向量；A(x；n，n′)≈(1+R(x))/(1-R(x))，n′为外界媒质的折射系数；当外界媒质是空气时，R(x)近似为R(x)≈-1.4399n^-2+0.7099n^-1+0.6681+0.0636n，其中n为空气的折射率；

2)有限元离散化处理：

光子在生物组织中传输的数学模型用下面的扩散方程表示：

$-\&dtri;\&CenterDot;(D(x,{\mathrm{\&tau;}}_l)\&dtri;\mathrm{\&Phi;}(x,{\mathrm{\&tau;}}_l))+{\mathrm{\&mu;}}_a(x,{\mathrm{\&tau;}}_l)\mathrm{\&Phi;}(x,{\mathrm{\&tau;}}_l)=S(x,{\mathrm{\&tau;}}_l)(x\&Element;\mathrm{\&Omega;})$

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

其中Ω是被测物体；

是被测物体的表面；Φ(x，τ_l)是光子流密度；S(x，τ_l)是光源密度；根据有限元理论，得到下面的弱解形式：

${\&Integral;}_{\mathrm{\&Omega;}}(D(x,{\mathrm{\&tau;}}_l)(\&dtri;\mathrm{\&Phi;}(x,{\mathrm{\&tau;}}_l))\&CenterDot;(\&dtri;\mathrm{\&Psi;}(x,{\mathrm{\&tau;}}_l))+{\mathrm{\&mu;}}_a(x,{\mathrm{\&tau;}}_l)\mathrm{\&Phi;}(x,{\mathrm{\&tau;}}_l)\mathrm{\&Psi;}(x,{\mathrm{\&tau;}}_l))\mathrm{dx}$

$+{\&Integral;}_{\&PartialD;\mathrm{\&Omega;}}\frac{1}{2A_n\left(x\right)}\mathrm{\&Phi;}(x,{\mathrm{\&tau;}}_l)\mathrm{\&Psi;}(x,{\mathrm{\&tau;}}_l)\mathrm{dx}={\&Integral;}_{\mathrm{\&Omega;}}S(x,{\mathrm{\&tau;}}_l)\mathrm{\&Psi;}(x,{\mathrm{\&tau;}}_l)\mathrm{dx}(\&ForAll;\mathrm{\&Psi;}(x,{\mathrm{\&tau;}}_l)\&Element;H^1\left(\mathrm{\&Omega;}\right))$

H¹(Ω)是Sobelev空间，Ψ(x，τ_l)是任意给定的试探函数；在自适应有限元分析的框架下，基于自适应的网格细分，把重建区域网格剖分表示为{T₁，...，T_k，...}的网格序列，其中网格剖分T_k包括

个离散单元；利用有限元方法，对弱解形式进行离散，得到下面矩阵方程：

(K_k(τ_l)+C_k(τ_l)+B_k(τ_l))Φ_k(τ_l)＝F_k(τ_l)S_k(τ_l)

矩阵元素的具体形式为：

$\left\{\begin{array}{c}{k}_{\mathrm{ij}}^k\left({\mathrm{\&tau;}}_l\right)={\&Integral;}_{\mathrm{\&Omega;}}D(x,{\mathrm{\&tau;}}_l)(\&dtri;{\mathrm{\&Psi;}}_i^k(x,{\mathrm{\&tau;}}_l))\&CenterDot;(\&dtri;{\mathrm{\&Psi;}}_j^k(x,{\mathrm{\&tau;}}_j))\mathrm{dx}\\ c_{\mathrm{ij}}^k\left({\mathrm{\&tau;}}_l\right)={\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&mu;}}_a(x,{\mathrm{\&tau;}}_l){\mathrm{\&Psi;}}_i^k(x,{\mathrm{\&tau;}}_l){\mathrm{\&Psi;}}_j^k(x,{\mathrm{\&tau;}}_l)\mathrm{dx}\\ b_{\mathrm{ij}}^k\left({\mathrm{\&tau;}}_l\right)={\&Integral;}_{\&PartialD;\mathrm{\&Omega;}}{\mathrm{\&Psi;}}_i^k(x,{\mathrm{\&tau;}}_l){\mathrm{\&Psi;}}_j^k(x,{\mathrm{\&tau;}}_l)/\left(2A(x,n,n^{\&prime;})\right)\mathrm{dx}\\ f_{\mathrm{ij}}^k\left({\mathrm{\&tau;}}_l\right)={\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&Psi;}}_i^k(x,{\mathrm{\&tau;}}_l){\mathrm{\&Psi;}}_j^k(x,{\mathrm{\&tau;}}_l)\mathrm{dx}\end{array}\right.$

令K_k(τ_l)+C_k(τ_l)+B_k(τ_l)＝M_k(τ_l)，M_k(τ_l)是稀疏正定矩阵，所以得到：

${\mathrm{\&Phi;}}_k\left({\mathrm{\&tau;}}_l\right)=M_k^{-1}\left({\mathrm{\&tau;}}_l\right)F_k\left({\mathrm{\&tau;}}_l\right)S_k\left({\mathrm{\&tau;}}_l\right)---\left(1\right)$ 通过删除

[M_k ^-1(τ_l)F_k(τ_l)]中对应非测量点的行得到下面的方程：

Φ_k(τ_l)＝A_k(τ_l)S_k(τ_l)           (2)

因此对于m个波段即得到m个单一波段的矩阵方程；

3)优化可行光源区域选择

I)在网格剖分T₁上，将上述含有m个单一波段的矩阵方程组合成含有m个波段的一个方程：

根据光源谱信息的性质，每个谱段τ_l上所占整个谱段的能量表示为S(τ_l)＝ω(τ_l)·S，其中ω(τ_l)≥0，ω(τ_l)通过实验测定；S表示光源总的能量；在网格剖分T₁上，考虑光源谱信息的性质，通过方程(2)得到：

AS₁＝Φ             (3)

其中 $\mathrm{\&Phi;}=\left[\begin{array}{c}{\mathrm{\&Phi;}}_1\left({\mathrm{\&tau;}}_1\right)\\ {\mathrm{\&Phi;}}_1\left({\mathrm{\&tau;}}_2\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Phi;}}_1\left({\mathrm{\&tau;}}_m\right)\end{array}\right],$ $A=\left[\begin{array}{c}\mathrm{\&omega;}\left({\mathrm{\&tau;}}_1\right)A_1\left({\mathrm{\&tau;}}_1\right)\\ \mathrm{\&omega;}\left({\mathrm{\&tau;}}_2\right)A_2\left({\mathrm{\&tau;}}_2\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \mathrm{\&omega;}\left({\mathrm{\&tau;}}_m\right)A_1\left({\mathrm{\&tau;}}_m\right)\end{array}\right]$

Φ₁(τ_l)和A₁(τ_l)都由方程(2)计算得到，S₁是网格剖分τ₁上的未知光源密度变量；在(3)式两边同乘A^T，得到：

A^TAS₁＝A^TΦ            (4)

A^T是A的转置矩阵，A^T·A是一个

的对称矩阵，所以(4)式是一个标准线形方程；

II)在网格剖分T₁上，确定优化光源可行区域：

具体方法为：在网格剖分T₁上，假定整个被测物体Ω为光源可行区域，给定未知光源密度的初始值S₁ ⁰，根据公式(4)，网格剖分T₁上的未知光源密度变量S₁用下面的关系式进行迭代求解：

$S_1^{n+1}=S_1^n+{\mathrm{\&alpha;}}_n{\mathrm{\&beta;}}_n$

其中步长搜索方向

n是迭代次数；

在迭代的过程中，用Φ₁ ^meas(τ_l)代替Φ₁(τ_l)，Φ₁ ^meas(τ_l)是网格剖分T₁上的边界测量值，由Φ^meas(x，τ_l)通过最近邻域插值得到；当||β_n||≤δ或n＞N_max，迭代停止，此时便得到了未知光源密度S₁的优化解S^*，优化解S^*的大致区域为Ω^*，这个区域Ω^*我们称之为优化的光源可行区域；在每次迭代的过程中，如果光源的密度小于零，则将其置为零；δ取值范围为10^-5-10^-2，N_max为最大迭代次数，取值不超过1000；

4)光源重建

I)在网格序列{T₁，...，T_k，...}上，根据选择好的优化光源可行区域Ω^*，重新将含有m个单一波段的矩阵方程组合成含有m个波段的一个方程：

考虑未知光源密度S_k(τ_l)与边界测量值Φ_k ^meas(τ_l)之间的线性关系，由方程(1)：

得到：

${\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_l\right)=G_k\left({\mathrm{\&tau;}}_l\right)S_k\left({\mathrm{\&tau;}}_l\right)$

其中G_k(τl)通过移出[M_k(τ_l)^-1F_k(τ_l)]中对应非测量点的行而得到；考虑到光源谱信息和优化可行光源区域Ω^*，得到：

${\mathrm{\&Phi;}}_k^{\mathrm{meas}}=G_k{W}_k{S}_k$

其中 $G_k=\left[\begin{array}{c}w\left({\mathrm{\&tau;}}_1\right)G_k\left({\mathrm{\&tau;}}_1\right)\\ w\left({\mathrm{\&tau;}}_2\right)G_k\left({\mathrm{\&tau;}}_2\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ w\left({\mathrm{\&tau;}}_m\right)G_k\left({\mathrm{\&tau;}}_m\right)\end{array}\right],$ ${\mathrm{\&Phi;}}_k^{\mathrm{meas}}=\left[\begin{array}{c}{\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_1\right)\\ {\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_2\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left({\mathrm{\&tau;}}_m\right)\end{array}\right],$ W_k＝Diag(w_k(11)，w_k(22)，…w_k(ii)，…)，

$w_{k\left(\mathrm{ii}\right)}=\left\{\begin{array}{cc}1& s_{k\left(i\right)}\&GreaterEqual;{\mathrm{\&gamma;}}^k{s}_k^{\max }\\ 0& s_{k\left(i\right)}<{\mathrm{\&gamma;}}^k{s}_k^{\max }\end{array}\right.$

G_k表示多谱下的边界测量值与光源位置量之间的关系，W_k为k层网格上的对角矩阵，w_k(ii)是对角矩阵W_k的第i个对角元素；Φ_k ^meas(τ_l)是波段τ_l在网格剖分T_k上的边界测量值，由Φ^meas(x，τ_l)通过最近邻域插值得到；当k＝1时，s_k(i)是通过迭代方法得到的优化解S^*，s_k ^max是优化解S^*的最大值；当k≥2时，s_k(i)和s_k ^max是网格剖分T_k上的未知光源密度的初始值和初始值的最大值，通过插值网格剖分T_k-1上的重建结果获得，即

S_k-1是在网格剖分T_k-1上重建得到的结果，I_k-1 ^k是插值因子，通过子单元继承父单元的重建结果值的方式来实现；γ^k是比例因子，是随k变化的分段常数；最终通过保留G_kW_k中对应可行光源区域的列，得到光源可行区域的未知光源密度变量S_k ^p与表面测量值Φ_k ^meas在网格剖分T_k上的关系：

$A_k{S}_k^p={\mathrm{\&Phi;}}_k^{\mathrm{meas}}$

II)利用Tikhonov正则化方法建立目标函数并优化从而得到重建结果：

上一步虽然得到了未知光源密度变量S_k ^p与表面测量值Φ_k ^meas之间的关系，由于A_k矩阵的病态，不能通过直接求解的方法得到重建结果，因此利用Tikhonov正则的方法，并考虑未知光源密度变量不能为负值的特性，得到在网格剖分T_k上的目标函数：

$\underset{0{\&le;s}_k^p\&le;s_{\sup }^k}{\min }{\mathrm{\&Theta;}}_k\left(S_k^p\right)=\{{\left|\right|A_k{S}_k^p-{\mathrm{\&Phi;}}_k^{\mathrm{meas}}\left|\right|}_{\mathrm{\&Lambda;}}+{\mathrm{\&lambda;}}_k\&CenterDot;{(S_k^p-S_k^{\mathrm{init}})}^T(S_k^p-S_k^{\mathrm{init}})\}$

考虑到多光谱以及非匀质模型会极大的增加计算量，因此采用谱梯度大规模优化算法对该目标函数进行优化；S_sup ^k是光源密度上界，经验取值为不大于1000，||V||_Λ＝V^TΛV，λ_k正则化系数，经验取值10^-9-10^-6，S_k ^init是网格剖分T_k上的未知光源密度初始值，取值范围为10^-7-10^-4；在判断优化过程是否中止时，采用当前优化梯度以及优化迭代次数N_k ^te作为评判准则，即当

或时，优化停止，得到重建结果，即光源可行区域的光源密度分布；上标te表示当前的迭代次数；ε_g是梯度阈值，其经验取值不超过10^-6；N_k ^max是最大优化迭代次数，经验取值不超过10⁴；

III)当优化完成后，判断重建是否停止：

利用优化得到的光源密度求解任意一个谱段τ_l上的光密度分布Φ_k ^c(τ_l)，并计算||Φ_k ^c(τ_l)-Φ_k ^meas(τ_l)||，当或k≥K_max，重建停止；Ф_k ^c(τ_l)是通过将优化求得的光源密度代入扩散方程而求得的边界处光子流密度，Φ_k ^meas(τ_l)是波段τ_l在网格剖分T_k上边界测量点的测量值，由Φ^meas(x，τ_l)通过最近邻域插值得到；ε_Ф经验取值范围10^-7-10^-9；K_max是网格最大剖分次数；如果不能满足

或k≥K_max，根据重建的结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域；通过以上方法对可行和禁止光源区域内的网格单元的选择，对区域剖分网格进行自适应的细分；并从网格剖分T_k转到T_k+1，然后转到第4)(I)步，直至重建停止。

Images