1 引言

生物发光断层成像(BLT)是一种新型的光学分子成像模态,它通过定位生物体内发光光源来定量地反映生物体分子细胞水平的生理和病理变化,在恶性肿瘤早期检测、药物开发及疗效评估等的应用中极具潜力^[1-2]。

BLT通过生物体表面获取到的有限测量信息来估计生物内部自发荧光光源的分布,是一个典型的不适定逆问题,近红外光在生物组织中传播时要经历多次散射和吸收,更增加了光源重建的难度^[3]。为了降低重建的病态性,通常需要在重建过程中融合可行区域^[4]、多光谱测量 ^[5]、结构先验信息^[6]等先验信息,以获得有意义的结果。受压缩感知理论的启发,Lu等^[7]将光源的稀疏分布先验应用于BLT重建。随后,各种不同的稀疏重建算法相继被提出,如不完全变量截断共轭梯度算法(IVTCG)^[8]、分段正交匹配追踪算法(StOMP)^[9]以及各种基于非凸L_p(0<p<1)正则化方法^[10-12]。研究者们通过这些稀疏重建算法克服了传统的基于L₂正则化重建算法导致的信号过平滑的问题,有效地改善了重建结果,但在稳定性和重建精度等方面仍有提升空间。

压缩感知领域的研究表明,L_1-2正则化相对于L₁和L₂正则化方法能得到更稳定和更精确的稀疏重构^[13-14]。为进一步改善光源重建结果,本文将非凸L_1-2正则化引入到BLT的逆问题求解中,并采用一种凸差分算法来解决非凸泛函最小化问题,在每一步迭代中通过带自适应惩罚项的交替方向乘子法来高效求解,并设计了数字鼠仿体实验来验证本文方法的有效性和稳健性。

2 方 法

2.1 BLT成像模型

近红外光在生物组织中传输时,辐射传输方程的扩散近似(DA)模型可以较为准确地描述这个过程。结合Robin边界条件的扩散近似模型可以描述为^[2]

-∇·D(r)∇Φ(r)+μa(r)Φ(r)=S(r),r∈ΩΦ(ξ)+2AnD(ξ)=0,ξ∈∂ξ,(1)

式中Φ(r)为区域Ω中点r位置处的光子能量流率,D(r)=1/[3(μ_a+μ'_s)]为点r处的扩散系数,μ_a和μ'_s分别为生物组织的吸收系数和约化散射系数,S(r)为内部光源的能量密度,ξ为体表边界处的点,A为边界失配因子,n为指向边界∂Ω外侧的单位法向量。

在给定的光学参数下,基于有限元方法(FEM)^[15]求解DA模型可得到表面测量值与内部光源间的线性对应关系:

Ax=B,(2)

式中:A为m×n刚度矩阵,且m≪n;x为光源的分布,为n维向量;B为通过测量得出的边界数据,是m维向量。

由于测量得到的边界数据B是有限的,而需要求解的x是一个高维向量,因此(2)式是一个典型的病态问题。为了在有限的测量下得到稳定的重建结果,结合BLT中光源在生物组织内稀疏分布的特点,采用L_1-2正则化将BLT光源重建转化为如下的优化问题^[16]:

minx∈Rn12‖Ax-B‖22+λ(‖x‖1-‖x‖2),(3)

式中λ为正则化参数,‖x‖₁-‖x‖₂为L_1-2正则化。

2.2 基于非凸L1-2正则化的BLT成像

为求解非凸的目标函数(3)式,采用凸差分算法^[17],即利用凸差分算法函数F(x)=G(x)-H(x)将(3)式转化为

F(x)=12‖Ax-B‖22+λ‖x‖1-λ‖x‖2,(4)

式中G(x)=12‖Ax-B‖22+λ‖x‖₁,H(x)=λ‖x‖₂。

当x≠0时,‖x‖₂可以表示成 xx‖2;当x=0时,0∈∂‖x‖₂,因此凸差分算法迭代求解可表示为

xn+1=argminx∈Rn12‖Ax-B‖22+λ‖x‖1,xn=0argminx∈Rn12‖Ax-B‖22+λ‖x‖1-<x,λxn‖xn‖2>,xn≠0,(5)

式中n为外迭代次数。

在每一步迭代中,还需要解决的L₁正则化的子问题表达如下

minx∈Rn12‖Ax-B‖22+<x,v>+λ(‖x‖1),(6)

式中v∈Rⁿ表示一个常数向量。利用增广拉格朗日法可将(6)式转化为

Lδ(x,z,y)=12‖Ax-B‖22+<x,v>+λ‖z‖1+yT(x-z)+δ2‖x-z‖22,(7)

式中y为拉格朗日乘子,δ>0为惩罚因子。

采用带自适应惩罚项的交替方向乘子法^[18]对(7)式依次进行迭代求解:

xl+1=argminxLδ(x,zl,yl)zl+1=argminzLδ(xl+1,z,yl)yl+1=yl+δ(xl+1-zl+1)δl+1=min(δmax,‖x‖δl)。(8)

z的更新基于如下软阈值算子实现:

S(z,r)]i=sgn(zi)max{‖zi‖-r,0}。(9)

利用下式对自适应惩罚因子进行变换:

δl+1=min(δmax,τδl),(10)

式中δ_max为δ_l的最大值,l为内迭代次数,τ为约束权重:

τ=τ0,δzl+1-zl‖‖xl+1‖<ε1,δzl+1-zl‖‖xl+1‖≥ε,(11)

式中τ₀≥1是给定的一个较小的常数,ε为初始设定的值。

基于非凸L_1-2正则化重建算法的具体步骤总结如下:

1) 初始化设定x⁰=0,x¹≠0,ε_outer>0,ε_inner>0。

2) 当‖xⁿ-x^n-1‖>ε_outer时,设σ= xn‖xn‖2,x₀=0,x₁=xⁿ,i=1,z_i=x_i,y_i=0。

3) 在满足上述条件下进行迭代求解。当‖x_l-x_l-1‖>ε_inner时,x_l+1=(A^TA+δ_lI)^-1(δ_lz_l-σ-y_l),z_l+1=shrink(x_l+1+y_l/δ_l,λ/δ_l),y_l+1=y_l+δ_l(x_l+1-z_l+1),δ_l+1=min(δ_max,τδ_l)。

4) 直到xⁿ=x_l求得一个满意的解,循环结束。

3 实验与结果

为了验证所提出的基于非凸L_1-2正则化的重建算法在BLT中的性能,本节设计了数字鼠模型上的单光源实验和双光源实验,并将其与光学分子影像中的常用稀疏重建算法——IVTCG算法^[2,8]、分离近似稀疏重构(SpaRSA)算法^[19]、StOMP算法^[9]进行对比,所有算法的正则化参数都是通过多次实验手动选择的。

数字鼠模型采用的是Digimouse^[20],在实验中截取数字鼠的躯干部分作为重建区域,躯干高34 mm,包括心脏、胃、肝脏、肾脏、肺等组织器官。单光源和双光源均设置为半径为0.5 mm、高为1 mm、初始能量密度为1 W/m³的小圆柱体。本研究采用中心定位误差(LE)评估所有算法的定位能力,采用算法耗时(T)评估算法的效率。实验中使用的光学参数如表1所示^[21]。

3.1 单光源数字鼠实验

在单光源实验中,光源放置在数字鼠的肝脏中,其中心点坐标设定为(9.2 mm,7.2 mm,19.2 mm)。为了得到模拟的表面测量数据,数字鼠模型在前向计算中被离散为一个包含14434个节点和80341个四面体单元的有限元网格。为提高算法效率和降低病态性,重建中采用了较为稀疏的逆向有限元网格,包含4600个节点和25069个四面体单元。图1分别展示了单光源数字鼠模型、前向网格以及仿真得到的表面光强分布。

图 1. (a)数字鼠模型及单光源设置;(b)前向仿真得到的表面光强分布

Fig. 1. (a) Digital mouse model and source setting in single-source case; (b) surface intensity distribution generated by forward simulation

下载图片 查看所有图片

图2是4种算法得到的重建结果的3D视图和在光源中心位置处的XY平面(Z=19.2 mm)的截面视图。表2给出了4种算法在单光源实验中的定量重建结果。从图2、表2中可以看出,尽管非凸L_1-2正则化算法比SpaRSA算法和StOMP算法的速度稍慢,但其总体中心定位误差最小。

图 2. 单光源重建结果对比。(a) L1-2、(b) IVTCG、(c) SpaRSA、(d) StOMP重建结果的3D视图;(e) L1-2、(f) IVTCG、(g) SpaRSA、(h) StOMP重建结果在真实光源中心Z=19.2 mm处的 XY平面的截面视图,其中红色圆圈代表真实光源

Fig. 2. Comparison of reconstruction results in single-source case. 3D views of the reconstruction results by (a) L1-2, (b) IVTCG, (c) SpaRSA, (d) StOMP; transverse views at Z=19.2 mm of the reconstruction results by (e) L1-2, (f) IVTCG, (g) SpaRSA, (h) StOMP, where the red circle represents the real source

下载图片 查看所有图片

3.2 双光源实验

双光源实验评估的不仅是算法的中心定位能力,还有算法对双光源的空间分辨能力。实验中,双光源放置在肾脏中,中心坐标分别设定为(9.6 mm,15.3 mm,27.9 mm)和(12.5 mm,15.3 mm,27.9 mm)。前向计算时,有限元网格包含14752个节点和82148个四面体单元,而在逆向重建中数字鼠模型被离散为仅包含3878个节点和21199个四面体单元的网格。图3展示了双光源数字鼠模型以及仿真得到的表面光强分布。

图4是双光源实验重建结果的比较,包含3D重建视图和在光源中心位置Z=27.9 mm处的XY平面视图。表3总结了4种算法在双光源实验中的定量重建结果。从图4和表3可以看出,在双光源实验中,非凸L_1-2正则化算法比其他几种算法的中心定位误差小,时间比SpaRSA算法和StOMP算法稍长,但重建的位置偏差更小。

3.3 稳健性测试

在实验测量中,噪声是不可避免的,对于不适定问题,算法的稳健性也是算法性能的重要评价指标。因此本节评估了不同噪声条件下,所提算法在单光源实验中的重建稳定性。设置的高斯白噪声水平分别为5%、10%、15%,鉴于噪声的随机性,在每个噪声水平下,重建算法独立运行10次,表4给出了不同噪声水平下单光源重建结果的平均值。由表4可以看出,随着噪声水平的增加,中心定位误差稍有增大,但总体重建结果比较稳定。

图 3. (a)数字鼠模及双光源位置;(b)前向仿真得到的表面光强分布

Fig. 3. (a) Digital mouse model and sources setting in double-source case; (b) surface intensity distribution generated by forward simulation

下载图片 查看所有图片

图 4. 双光源重建结果对比。(a) L1-2、(b) IVTCG、(c) SpaRSA、(d) StOMP重建结果的3D视图;(e) L1-2、(f) IVTCG、(g) SpaRSA、(h) StOMP重建结果在真实光源中心Z=27.9 mm处的XY平面的截面视图,其中红色圆圈代表真实光源

Fig. 4. Comparison of reconstruction results in double-source case. 3D views of the reconstruction results by (a) L1-2, (b) IVTCG, (c) SpaRSA, (d) StOMP; transverse views at Z=27.9 mm of the reconstruction results by (e) L1-2, (f) IVTCG, (g) SpaRSA, (h) StOMP, where the red circle represents the real source

下载图片 查看所有图片

4 结论

将一种新型的非凸L_1-2正则化方法引入到BLT光源重建中,并采用凸差分算法来解决非凸泛函最小化问题,在每一步迭代中采用带自适应惩罚项的交替方向乘子法进行高效求解。非匀质数字鼠模型的仿真结果表明,在单光源和双光源实验中,基于非凸L_1-2正则化算法的中心定位误差都控制在1 mm以内,优于其他几种对比算法,在不同的噪声条件下也获得了相对稳定的重建结果。但从仿真结果中也可以看出,该算法得到的重建能量密度偏小,下一步研究将结合自适应有限元方法对算法进行改进,以降低重建能量的偏差,进一步改善重建结果。