臧青楊，陳春曉，楊俊豪，李東升

(南京航空航天大學(xué)生物醫(yī)學(xué)工程系，南京 211106)

1 引 言

動(dòng)態(tài)熒光分子成像(dynamic fluorescence molecular imaging, D-FMI)是利用熒光分子探針在體內(nèi)不同器官中的分布隨時(shí)間變化的動(dòng)態(tài)性差異進(jìn)行成像的技術(shù)[1-2]。D-FMI在藥代動(dòng)力學(xué)、腫瘤診斷和病理信息監(jiān)測等領(lǐng)域有重要的應(yīng)用前景。動(dòng)態(tài)熒光分子重建是一個(gè)嚴(yán)重的欠定性問題，壓縮感知理論是有效解決重建欠定性的理論方法之一[3]。基于壓縮感知的范數(shù)優(yōu)化算法是熒光分子逆向重建的常用算法，該算法通過范數(shù)正則化項(xiàng)表征待求解模型稀疏度。然而，采用范數(shù)優(yōu)化算法解決在體熒光目標(biāo)的稀疏重建時(shí)，由于采集到的動(dòng)態(tài)熒光信號(hào)幀之間相互獨(dú)立、重構(gòu)信號(hào)維度高及沒有充分考慮數(shù)據(jù)在時(shí)空域中的相關(guān)性與塊稀疏特性，導(dǎo)致重建稀疏度不足、定位精度低等問題[4]。

相較于傳統(tǒng)的壓縮感知算法，基于概率統(tǒng)計(jì)的貝葉斯學(xué)習(xí)(sparse bayesian learning, SBL)具有重建信號(hào)稀疏性好、定位精度高等優(yōu)勢[5]。多目標(biāo)動(dòng)態(tài)熒光分子重建屬于多觀測向量(multiple measurement vectors, MMV)模型聯(lián)合稀疏重建問題。由于動(dòng)態(tài)熒光信號(hào)共享相同稀疏結(jié)構(gòu)，塊稀疏模型能利用信號(hào)的分塊稀疏特性有效挖掘多觀測信號(hào)時(shí)空相關(guān)性信息[6]，因此，本研究采用基于MMV模型的塊稀疏貝葉斯學(xué)習(xí)(block sparse bayesian learning, BSBL)解決多目標(biāo)動(dòng)態(tài)熒光分子的重建問題。

2 方法

2.1 動(dòng)態(tài)熒光光子傳輸模型

在熒光分子斷層成像中，通常使用擴(kuò)散近似方程作為描述光子在生物體組織中傳輸過程的數(shù)學(xué)模型[7]：

(1)

在生物組織邊界應(yīng)用Robin邊界條件，并結(jié)合有限元法求解式(1)，可以建立體表光通量密度Φm與光源能量密度S之間的關(guān)系：

Φm=AS+v

(2)

其中，S為n維列向量，n是有限元離散后的網(wǎng)格節(jié)點(diǎn)個(gè)數(shù)；Φm為m維列向量，m是體表節(jié)點(diǎn)個(gè)數(shù)；A為m×n的系數(shù)矩陣,主要取決于生物體的結(jié)構(gòu)分布和光學(xué)特性；v為未知噪聲向量。

動(dòng)態(tài)熒光分子成像技術(shù)是在靜態(tài)成像的基礎(chǔ)上加入了時(shí)間維度，通過采集相繼離散時(shí)間點(diǎn)下的熒光信號(hào)，根據(jù)式(2)建立動(dòng)態(tài)熒光光子傳輸模型：

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

(3)

其中，Φm(t)為不同觀測時(shí)間下的體表光通量密度；S(t)為不同觀測時(shí)間下的光源能量密度；L為動(dòng)態(tài)熒光信號(hào)的采樣次數(shù)，[L]表示1到L的整數(shù)集；v(t)為不同觀測時(shí)間下引入的未知噪聲。

2.2 BSBL算法

經(jīng)典的SBL算法是假設(shè)待求解向量S的各個(gè)元素相互獨(dú)立，通過給S的各個(gè)元素的權(quán)值系數(shù)賦予先驗(yàn)條件概率分布加以限制，并引入相關(guān)向量機(jī)(relevance vector machines, RVM)稀疏模型，在優(yōu)化迭代過程中剔除權(quán)值系數(shù)趨于零的訓(xùn)練樣本，從而獲取模型的稀疏解[8-9]。

BSBL算法最初在2013年由Zhang等提出，它是一種以分解稀疏向量中每個(gè)塊內(nèi)的元素之間的幅值相關(guān)性作為結(jié)構(gòu)先驗(yàn)信息的新型塊稀疏重構(gòu)算法[10-11]。BSBL算法首次揭示了塊內(nèi)相關(guān)對(duì)塊稀疏信號(hào)重構(gòu)性能的影響，文獻(xiàn)[12]證明了在無噪情況下BSBL求得的全局解即為真正的最稀疏解。

BSBL算法中假定待求解稀疏向量S可以劃分為g個(gè)非重疊塊結(jié)構(gòu)，并且S的非零元素集中分布于少數(shù)非零塊中，即：

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

(4)

式(2)和式(4)組成了塊稀疏壓縮感知模型，在BSBL算法中，各個(gè)信號(hào)塊可以具備不同的大小。與之對(duì)應(yīng)，系數(shù)矩陣A可分為g個(gè)塊矩陣：

A=[A1，…，Ag]

(5)

BSBL算法假設(shè)第i個(gè)信號(hào)塊Si服從多元高斯分布：

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

(6)

其中，γi為一個(gè)非負(fù)尺度參數(shù)，用于控制信號(hào)塊Si的稀疏度。當(dāng)γi=0時(shí)，對(duì)應(yīng)信號(hào)塊Si則為0。在學(xué)習(xí)過程中，受自動(dòng)相關(guān)性確定機(jī)制影響，大部分γi將趨于0，通過設(shè)定合適的閾值將趨于0的γi置為0，從而促成了解的塊稀疏；Bi為一個(gè)未知的正定矩陣，表征著該塊內(nèi)的元素之間的相關(guān)結(jié)構(gòu)模型信息。

在假設(shè)信號(hào)塊與塊之間相互獨(dú)立的前提下，待求解向量S的先驗(yàn)分布可以建模為：

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

(7)

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

同時(shí)假設(shè)觀測噪聲v服從p(v；λ)=N(0，λI)分布。根據(jù)貝葉斯準(zhǔn)則可得出S的后驗(yàn)分布為：

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)優(yōu)化算法對(duì)式進(jìn)行求解，得到所有參量的學(xué)習(xí)規(guī)則為：

(10)

(11)

(12)

2.3 基于BSBL的多目標(biāo)動(dòng)態(tài)熒光分子重建

現(xiàn)有的BSBL算法多針對(duì)于單觀測向量(single measurement vector, SMV)模型，但在實(shí)際應(yīng)用中常需從MMV模型中重構(gòu)稀疏信號(hào)，本研究的多目標(biāo)動(dòng)態(tài)熒光分子重建即屬于MMV模型聯(lián)合稀疏重構(gòu)問題。研究表明，相比于SMV模型，MMV模型更能準(zhǔn)確估計(jì)出稀疏解向量非零行的位置，重構(gòu)性能更優(yōu)[13-14]。

MMV模型聯(lián)合稀疏重構(gòu)問題需假定重構(gòu)源信號(hào)是K稀疏的，而且不同源信號(hào)之間共享相同的稀疏結(jié)構(gòu)[15]。MMV稀疏重構(gòu)模型見圖1，源信號(hào)s(t)的每一個(gè)列向量都是塊稀疏的，同時(shí)不同列向量之間共享相同的稀疏結(jié)構(gòu)且具備時(shí)域相關(guān)性，因此信號(hào)塊si的劃分可以利用到信號(hào)的時(shí)空相關(guān)性和塊稀疏特性。

在假設(shè)信號(hào)塊與塊之間相互獨(dú)立的前提下，信號(hào)塊Si符合矩陣正態(tài)分布，則動(dòng)態(tài)熒光信號(hào)S(t)的先驗(yàn)分布可以建模為：

圖1 MMV稀疏重構(gòu)模型Fig 1 MMV sparse reconstruction model p(S(t)|{Γ，B})=MN(S(t)；0，Γ，B)

(13)

其中，Γ=diag-1(Di)，Di為正定對(duì)稱矩陣，用來描述稀疏信號(hào)塊Si的相關(guān)結(jié)構(gòu)信息。

同時(shí)假設(shè)S(t)的所有列向量滿足同一相關(guān)模型B，在式(3)模型下，根據(jù)貝葉斯準(zhǔn)則對(duì)觀測信號(hào)Φm(t)進(jìn)行建?？傻茫?/p>

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

(14)

為了估計(jì)參量Di與B，采用第二類最大似然估計(jì)方法得到代價(jià)函數(shù)L：

(15)

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

采用優(yōu)化算法對(duì)上式代價(jià)函數(shù)進(jìn)行求解，可得待估計(jì)參量的學(xué)習(xí)規(guī)則為：

(16)

(17)

3 實(shí)驗(yàn)結(jié)果

3.1 數(shù)值仿真實(shí)驗(yàn)

本研究數(shù)值仿真實(shí)驗(yàn)設(shè)置見圖2。圖2(a)為數(shù)值仿真幾何模型，大圓柱體作為生物組織，半徑為1.5 cm，高度為3 cm，內(nèi)置三個(gè)小圓柱體作為多目標(biāo)熒光光源，小圓柱半徑為0.2 cm，高度為0.4 cm，中心坐標(biāo)分別為(-0.7,0,1.5)、(0.5,0.5,1.5)和(0.5,-0.5,1.5)，多目標(biāo)熒光光源的動(dòng)態(tài)濃度設(shè)置見圖2(b)。

為獲取體表測量數(shù)據(jù)，將幾何模型進(jìn)行有限元四面體網(wǎng)格剖分，通過COMSOL前向仿真獲取體表光通量密度Φm(t)。M-FOCUSS(multiple fOCal underdetermined system solver)算法是基于MMV模型的稀疏信號(hào)重構(gòu)算法，它利用Lp范數(shù)模型對(duì)待求解列向量進(jìn)行約束，相較于L0和L1范數(shù)優(yōu)化模型有著更優(yōu)的信號(hào)重構(gòu)性能和更高的計(jì)算效率。使用M-FOCUSS算法和本研究塊稀疏貝葉斯方法進(jìn)行重建，重建結(jié)果見圖3。

從圖3的對(duì)比結(jié)果可以看出，M-FOCUSS算法雖然可以重建出光源的大致位置與強(qiáng)度，但存在重建精確度和稀疏度不足的問題，本研究方法可以獲取更精確、更稀疏的重建結(jié)果。

通過非負(fù)矩陣分解(non-negative matrix factorization, NMF)算法對(duì)多目標(biāo)動(dòng)態(tài)熒光分子重建結(jié)果進(jìn)行分離，獲取多目標(biāo)光源的動(dòng)態(tài)濃度分布?xì)w一化曲線，見圖4。從圖3和圖4可以看出BSBL算法對(duì)動(dòng)態(tài)熒光濃度的重構(gòu)稀疏度更好、精度更高。

圖2數(shù)值仿真實(shí)驗(yàn)設(shè)置
(a).仿真幾何模型；(b).動(dòng)態(tài)熒光濃度設(shè)置

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

圖3 t=4時(shí),數(shù)值仿真重建結(jié)果(a)、(b).真實(shí)光源位置和濃度；(c)、(d).M-FOCUSS算法重建結(jié)果；(e)、(f).塊稀疏貝葉斯重建結(jié)果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 生物組織實(shí)驗(yàn)

本研究設(shè)計(jì)了在豬肉組織內(nèi)部植入雙腫瘤源的實(shí)驗(yàn)對(duì)重建算法進(jìn)行驗(yàn)證，動(dòng)態(tài)腫瘤源由不同濃度的吲哚菁綠(indocyanine green, ICG)溶液密封在透明導(dǎo)管構(gòu)成。豬肉組織實(shí)驗(yàn)設(shè)置見圖5，其中豬肉組織的長為5.5 cm，寬為3 cm，高度為3 cm，內(nèi)置腫瘤源半徑為0.1 cm，高度為0.2 cm，中心坐標(biāo)分別為(1.7,1.5,2.1)和(3.8,1.5,2.1)，雙腫瘤源動(dòng)態(tài)濃度設(shè)置見圖5(b)。

體表光學(xué)圖像由課題組搭建的小動(dòng)物光學(xué)成像系統(tǒng)AOIS獲取，采集的二維動(dòng)態(tài)熒光圖像見圖6。

圖4動(dòng)態(tài)熒光分離結(jié)果
(a)、(b)、(c)分別為tube1、tube2和tube3的動(dòng)態(tài)熒光濃度歸一化曲線

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

圖5 豬肉組織實(shí)驗(yàn)設(shè)置(a).豬肉組織幾何模型；(b).雙腫瘤源動(dòng)態(tài)熒光濃度設(shè)置Fig 5 Pork tissue experiment setup(a).pork tissue geometric model；(b).dual tumor dynamic fluorescence concentration setting

圖6 二維動(dòng)態(tài)熒光圖像Fig 6 Two dimensional dynamic fluorescence image

使用M-FOCUSS算法和本研究方法對(duì)豬肉實(shí)驗(yàn)數(shù)據(jù)進(jìn)行重建，重建結(jié)果見圖7。

對(duì)豬肉組織實(shí)驗(yàn)重建結(jié)果進(jìn)行動(dòng)態(tài)濃度分離，獲取動(dòng)態(tài)腫瘤源的動(dòng)態(tài)濃度分布曲線，如圖8所示。

圖7t=4時(shí),豬肉組織實(shí)驗(yàn)重建結(jié)果
(a)、(b).真實(shí)腫瘤位置和濃度；(c)、(d).M-FOCUSS算法重建結(jié)果；(e)、(f).塊稀疏貝葉斯重建結(jié)果

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

圖8 豬肉組織實(shí)驗(yàn)動(dòng)態(tài)熒光分離結(jié)果(a)、(b)分別為tube1、tube2的動(dòng)態(tài)熒光濃度歸一化曲線Fig 8 Dynamic fluorescence separation results of pork tissue experiment(a)，(b) are normalized curves of dynamic fluorescence concentration for tube1 and tube2 respectively.

本研究采用定位誤差和濃度偏差對(duì)M-FOCUSS算法和本研究算法重建結(jié)果進(jìn)行評(píng)價(jià)，見表1。定位誤差為重建光源中心與真實(shí)光源中心差的歐式距離。從圖7和表1可以看出本研究算法在腫瘤定位精度和濃度偏差方面占據(jù)優(yōu)勢，重建結(jié)果收斂性也比較好。

表1 重建效果對(duì)比

4 結(jié)論

本研究針對(duì)多目標(biāo)動(dòng)態(tài)熒光分子逆向重建存在稀疏性不足、定位精度低的問題，提出了基于塊稀疏貝葉斯學(xué)習(xí)的重建方法。該方法充分挖掘了信號(hào)的時(shí)空相關(guān)性信息和分塊稀疏特性，實(shí)現(xiàn)多目標(biāo)動(dòng)態(tài)熒光分子的稀疏重建。通過數(shù)值仿真和生物組織實(shí)驗(yàn)驗(yàn)證，本研究重建算法相較于傳統(tǒng)壓縮感知重構(gòu)算法具有更高的定位精度和更好的魯棒性。