基于块稀疏贝叶斯学习的多目标动态荧光分子重建

臧青杨，陈春晓，杨俊豪，李东升

(南京航空航天大学生物医学工程系，南京 211106)

1 引 言

动态荧光分子成像(dynamic fluorescence molecular imaging, D-FMI)是利用荧光分子探针在体内不同器官中的分布随时间变化的动态性差异进行成像的技术[1-2]。D-FMI在药代动力学、肿瘤诊断和病理信息监测等领域有重要的应用前景。动态荧光分子重建是一个严重的欠定性问题，压缩感知理论是有效解决重建欠定性的理论方法之一[3]。基于压缩感知的范数优化算法是荧光分子逆向重建的常用算法，该算法通过范数正则化项表征待求解模型稀疏度。然而，采用范数优化算法解决在体荧光目标的稀疏重建时，由于采集到的动态荧光信号帧之间相互独立、重构信号维度高及没有充分考虑数据在时空域中的相关性与块稀疏特性，导致重建稀疏度不足、定位精度低等问题[4]。

相较于传统的压缩感知算法，基于概率统计的贝叶斯学习(sparse bayesian learning, SBL)具有重建信号稀疏性好、定位精度高等优势[5]。多目标动态荧光分子重建属于多观测向量(multiple measurement vectors, MMV)模型联合稀疏重建问题。由于动态荧光信号共享相同稀疏结构，块稀疏模型能利用信号的分块稀疏特性有效挖掘多观测信号时空相关性信息[6]，因此，本研究采用基于MMV模型的块稀疏贝叶斯学习(block sparse bayesian learning, BSBL)解决多目标动态荧光分子的重建问题。

2 方法

2.1 动态荧光光子传输模型

在荧光分子断层成像中，通常使用扩散近似方程作为描述光子在生物体组织中传输过程的数学模型[7]：

(1)

在生物组织边界应用Robin边界条件，并结合有限元法求解式(1)，可以建立体表光通量密度Φm与光源能量密度S之间的关系：

Φm=AS+v

(2)

其中，S为n维列向量，n是有限元离散后的网格节点个数；Φm为m维列向量，m是体表节点个数；A为m×n的系数矩阵,主要取决于生物体的结构分布和光学特性；v为未知噪声向量。

动态荧光分子成像技术是在静态成像的基础上加入了时间维度，通过采集相继离散时间点下的荧光信号，根据式(2)建立动态荧光光子传输模型：

Φm(t)=AS(t)+v(t)，t∈[L]

(3)

其中，Φm(t)为不同观测时间下的体表光通量密度；S(t)为不同观测时间下的光源能量密度；L为动态荧光信号的采样次数，[L]表示1到L的整数集；v(t)为不同观测时间下引入的未知噪声。

2.2 BSBL算法

经典的SBL算法是假设待求解向量S的各个元素相互独立，通过给S的各个元素的权值系数赋予先验条件概率分布加以限制，并引入相关向量机(relevance vector machines, RVM)稀疏模型，在优化迭代过程中剔除权值系数趋于零的训练样本，从而获取模型的稀疏解[8-9]。

BSBL算法最初在2013年由Zhang等提出，它是一种以分解稀疏向量中每个块内的元素之间的幅值相关性作为结构先验信息的新型块稀疏重构算法[10-11]。BSBL算法首次揭示了块内相关对块稀疏信号重构性能的影响，文献[12]证明了在无噪情况下BSBL求得的全局解即为真正的最稀疏解。

BSBL算法中假定待求解稀疏向量S可以划分为g个非重叠块结构，并且S的非零元素集中分布于少数非零块中，即：

S=[S1，…，Sd1，…，Sdg-1+1，…，Sdg]T

(4)

式(2)和式(4)组成了块稀疏压缩感知模型，在BSBL算法中，各个信号块可以具备不同的大小。与之对应，系数矩阵A可分为g个块矩阵：

A=[A1，…，Ag]

(5)

BSBL算法假设第i个信号块Si服从多元高斯分布：

p(Si|γi，Bi)=N(Si；0，γiBi)

(6)

其中，γi为一个非负尺度参数，用于控制信号块Si的稀疏度。当γi=0时，对应信号块Si则为0。在学习过程中，受自动相关性确定机制影响，大部分γi将趋于0，通过设定合适的阈值将趋于0的γi置为0，从而促成了解的块稀疏；Bi为一个未知的正定矩阵，表征着该块内的元素之间的相关结构模型信息。

在假设信号块与块之间相互独立的前提下，待求解向量S的先验分布可以建模为：

p(S|{γi，Bi})=N(S；0，Γ)

(7)

其中，Γ=diag-1(γ1B1，…，γgBg)。

同时假设观测噪声v服从p(v；λ)=N(0，λI)分布。根据贝叶斯准则可得出S的后验分布为：

p(S|Φm，{γi，Bi}，λ)=N(S；μ，∑)

(8)

其中，∑-1=Γ-1+ATλ-1A，μ=∑ATλ-1Φm。

+ΦmT(λI+AΓAT)-1Φm

(9)

采用快速边缘似然最大化(fast marginal likelihood maximization, FMLM)优化算法对式进行求解，得到所有参量的学习规则为：

(10)

(11)

(12)

2.3 基于BSBL的多目标动态荧光分子重建

现有的BSBL算法多针对于单观测向量(single measurement vector, SMV)模型，但在实际应用中常需从MMV模型中重构稀疏信号，本研究的多目标动态荧光分子重建即属于MMV模型联合稀疏重构问题。研究表明，相比于SMV模型，MMV模型更能准确估计出稀疏解向量非零行的位置，重构性能更优[13-14]。

MMV模型联合稀疏重构问题需假定重构源信号是K稀疏的，而且不同源信号之间共享相同的稀疏结构[15]。MMV稀疏重构模型见图1，源信号s(t)的每一个列向量都是块稀疏的，同时不同列向量之间共享相同的稀疏结构且具备时域相关性，因此信号块si的划分可以利用到信号的时空相关性和块稀疏特性。

在假设信号块与块之间相互独立的前提下，信号块Si符合矩阵正态分布，则动态荧光信号S(t)的先验分布可以建模为：

图1 MMV稀疏重构模型Fig 1 MMV sparse reconstruction model p(S(t)|{Γ，B})=MN(S(t)；0，Γ，B)

(13)

其中，Γ=diag-1(Di)，Di为正定对称矩阵，用来描述稀疏信号块Si的相关结构信息。

同时假设S(t)的所有列向量满足同一相关模型B，在式(3)模型下，根据贝叶斯准则对观测信号Φm(t)进行建模可得：

p(Φm(t)|λ，B)
=MN(Φm(t)；AS(t)，λIM，B)

(14)

为了估计参量Di与B，采用第二类最大似然估计方法得到代价函数L：

(15)

其中，CA=λIM+AΓAT。

采用优化算法对上式代价函数进行求解，可得待估计参量的学习规则为：

(16)

(17)

3 实验结果

3.1 数值仿真实验

本研究数值仿真实验设置见图2。图2(a)为数值仿真几何模型，大圆柱体作为生物组织，半径为1.5 cm，高度为3 cm，内置三个小圆柱体作为多目标荧光光源，小圆柱半径为0.2 cm，高度为0.4 cm，中心坐标分别为(-0.7,0,1.5)、(0.5,0.5,1.5)和(0.5,-0.5,1.5)，多目标荧光光源的动态浓度设置见图2(b)。

为获取体表测量数据，将几何模型进行有限元四面体网格剖分，通过COMSOL前向仿真获取体表光通量密度Φm(t)。M-FOCUSS(multiple fOCal underdetermined system solver)算法是基于MMV模型的稀疏信号重构算法，它利用Lp范数模型对待求解列向量进行约束，相较于L0和L1范数优化模型有着更优的信号重构性能和更高的计算效率。使用M-FOCUSS算法和本研究块稀疏贝叶斯方法进行重建，重建结果见图3。

从图3的对比结果可以看出，M-FOCUSS算法虽然可以重建出光源的大致位置与强度，但存在重建精确度和稀疏度不足的问题，本研究方法可以获取更精确、更稀疏的重建结果。

通过非负矩阵分解(non-negative matrix factorization, NMF)算法对多目标动态荧光分子重建结果进行分离，获取多目标光源的动态浓度分布归一化曲线，见图4。从图3和图4可以看出BSBL算法对动态荧光浓度的重构稀疏度更好、精度更高。

图2数值仿真实验设置
(a).仿真几何模型；(b).动态荧光浓度设置

Fig2Numericalsimulationexperimentsetup
(a).geometricmodel；(b).dynamicfluorescenceconcentrationsettings

图3 t=4时,数值仿真重建结果(a)、(b).真实光源位置和浓度；(c)、(d).M-FOCUSS算法重建结果；(e)、(f).块稀疏贝叶斯重建结果Fig 3 Results of numerical simulation reconstruction, when t=4 (a)、(b).the position and intensity of true light source；(c)、(d).the reconstruction results of M-FOCUSS algorithm；(e)、(f).the reconstruction results of BSBL

3.2 生物组织实验

本研究设计了在猪肉组织内部植入双肿瘤源的实验对重建算法进行验证，动态肿瘤源由不同浓度的吲哚菁绿(indocyanine green, ICG)溶液密封在透明导管构成。猪肉组织实验设置见图5，其中猪肉组织的长为5.5 cm，宽为3 cm，高度为3 cm，内置肿瘤源半径为0.1 cm，高度为0.2 cm，中心坐标分别为(1.7,1.5,2.1)和(3.8,1.5,2.1)，双肿瘤源动态浓度设置见图5(b)。

体表光学图像由课题组搭建的小动物光学成像系统AOIS获取，采集的二维动态荧光图像见图6。

图4动态荧光分离结果
(a)、(b)、(c)分别为tube1、tube2和tube3的动态荧光浓度归一化曲线

Fig4Resultsofdynamicfluorescenceseparation
(a)，(b)，(c)，arenormalizedcurvesofdynamicfluorescenceconcentrationfortube1,tube2andtube3respectively

图5 猪肉组织实验设置(a).猪肉组织几何模型；(b).双肿瘤源动态荧光浓度设置Fig 5 Pork tissue experiment setup(a).pork tissue geometric model；(b).dual tumor dynamic fluorescence concentration setting

图6 二维动态荧光图像Fig 6 Two dimensional dynamic fluorescence image

使用M-FOCUSS算法和本研究方法对猪肉实验数据进行重建，重建结果见图7。

对猪肉组织实验重建结果进行动态浓度分离，获取动态肿瘤源的动态浓度分布曲线，如图8所示。

图7t=4时,猪肉组织实验重建结果
(a)、(b).真实肿瘤位置和浓度；(c)、(d).M-FOCUSS算法重建结果；(e)、(f).块稀疏贝叶斯重建结果

Fig7Reconstructionresultsofporktissue,whent=4
(a)and(b)arethepositionandintensityoftruetumour；
(c)and(d)arethereconstructionresultsofM-FOCUSSalgorithm；(e)and(f)arethereconstructionresultsofBSBL

图8 猪肉组织实验动态荧光分离结果(a)、(b)分别为tube1、tube2的动态荧光浓度归一化曲线Fig 8 Dynamic fluorescence separation results of pork tissue experiment(a)，(b) are normalized curves of dynamic fluorescence concentration for tube1 and tube2 respectively.

本研究采用定位误差和浓度偏差对M-FOCUSS算法和本研究算法重建结果进行评价，见表1。定位误差为重建光源中心与真实光源中心差的欧式距离。从图7和表1可以看出本研究算法在肿瘤定位精度和浓度偏差方面占据优势，重建结果收敛性也比较好。

表1 重建效果对比

4 结论

本研究针对多目标动态荧光分子逆向重建存在稀疏性不足、定位精度低的问题，提出了基于块稀疏贝叶斯学习的重建方法。该方法充分挖掘了信号的时空相关性信息和分块稀疏特性，实现多目标动态荧光分子的稀疏重建。通过数值仿真和生物组织实验验证，本研究重建算法相较于传统压缩感知重构算法具有更高的定位精度和更好的鲁棒性。