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

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

1 引 言

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

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

2 方法

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

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

(1)

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

Φm=AS+v

(2)

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

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

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

(3)

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

2.2 BSBL算法

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

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

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

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

(4)

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

A=[A1，…，Ag]

(5)

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

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

(6)

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

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

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

(7)

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

同時假設(shè)觀測噪聲v服從p(v；λ)=N(0，λI)分布。根據(jù)貝葉斯準則可得出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)優(yōu)化算法對式進行求解，得到所有參量的學(xué)習(xí)規(guī)則為：

(10)

(11)

(12)

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

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

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

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

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

(13)

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

同時假設(shè)S(t)的所有列向量滿足同一相關(guān)模型B，在式(3)模型下，根據(jù)貝葉斯準則對觀測信號Φm(t)進行建模可得：

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

(14)

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

(15)

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

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

(16)

(17)

3 實驗結(jié)果

3.1 數(shù)值仿真實驗

本研究數(shù)值仿真實驗設(shè)置見圖2。圖2(a)為數(shù)值仿真幾何模型，大圓柱體作為生物組織，半徑為1.5 cm，高度為3 cm，內(nèi)置三個小圓柱體作為多目標(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)熒光光源的動態(tài)濃度設(shè)置見圖2(b)。

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

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

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

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

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

圖3 t=4時,數(shù)值仿真重建結(jié)果(a)、(b).真實光源位置和濃度；(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è)計了在豬肉組織內(nèi)部植入雙腫瘤源的實驗對重建算法進行驗證，動態(tài)腫瘤源由不同濃度的吲哚菁綠(indocyanine green, ICG)溶液密封在透明導(dǎo)管構(gòu)成。豬肉組織實驗設(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)，雙腫瘤源動態(tài)濃度設(shè)置見圖5(b)。

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

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

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

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

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

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

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

圖7t=4時,豬肉組織實驗重建結(jié)果
(a)、(b).真實腫瘤位置和濃度；(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 豬肉組織實驗動態(tài)熒光分離結(jié)果(a)、(b)分別為tube1、tube2的動態(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.

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

表1 重建效果對比

4 結(jié)論

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