一種基于改進(jìn)多光譜的稀疏角CB-XLCT成像方法

張文元，海琳琦，劉英杰，張海波

(1.西安航空職業(yè)技術(shù)學(xué)院航空維修工程學(xué)院，陜西 西安 710089； 2.西北大學(xué)信息科學(xué)與技術(shù)學(xué)院，陜西 西安 710127)

0 引 言

X射線發(fā)光斷層成像(X-ray Luminescence Computed Tomography, XLCT)作為一種新的三維分子雙模成像技術(shù)，已成為生物光學(xué)與醫(yī)用成像領(lǐng)域的前沿?zé)狳c[1]。結(jié)合納米發(fā)光材料的光學(xué)特性，利用CT技術(shù)中的X射線源激發(fā)，發(fā)光材料受照射會發(fā)出近紅外光，進(jìn)而借助高靈敏度電耦合器件(Charge Coupled Device, CCD)采集透射出生物體表的微弱光信號，基于逆問題建模與成像算法可在分子層面對早期腫瘤進(jìn)行有效檢測[2]。

根據(jù)CT中X射線源的激發(fā)方式，XLCT技術(shù)可分為筆束型XLCT(Pencil-Beam XLCT, PB-XLCT)、扇形束XLCT(Fan-Beam XLCT, FB-XLCT)以及錐束型XLCT(Cone-Beam XLCT, CB-XLCT)[1,3-5]。其中，PB-XLCT和FB-XLCT分別利用筆束形X射線源與扇形X射線源各自的特點，利用“激發(fā)先驗”可獲得較高的成像分辨率，但由于系統(tǒng)耗時過長，不利于腫瘤的快速成像。隨后，CB-XLCT被提出，其利用錐束型射線源掃描時間短和X劑量利用率高的優(yōu)點，為XLCT應(yīng)用于早期腫瘤的實時成像提供了可行性[6-11]。

類似于傳統(tǒng)的生物發(fā)光斷層成像(BioLuminescence Tomography, BLT)[12-13]和熒光分子斷層成像(Fluorescence Molecular Tomography, FMT)[14-15]，XLCT成像問題在數(shù)學(xué)上屬于不適定逆問題[16]。BLT和FMT等經(jīng)典的重建方法大都可直接應(yīng)用于傳統(tǒng)CB-XLCT成像中。Chen等人[17]提出一種結(jié)合靶標(biāo)深度位置的自適應(yīng)Bregman算法，在算法中融入X射線衰減先驗信息，有效降低了成像逆問題病態(tài)性。Pu等人[11]將主成分分析應(yīng)用于CB-XLCT成像，通過實驗驗證了所提方法的有效性。Zhang等人[5]提出了一種基于高斯馬爾可夫隨機場的貝葉斯成像方法，用于傳統(tǒng)CB-XLCT的逆問題求解。曲璇[6]提出了一種基于多網(wǎng)格的可行區(qū)域策略，但多網(wǎng)格的多次重建間接增加了計算成本。隨后，稀疏角CB-XLCT(Sparse-view CB-XLCT)技術(shù)的出現(xiàn)，為CB-XLCT技術(shù)的實時成像轉(zhuǎn)化進(jìn)程提供了進(jìn)一步促進(jìn)。然而，由于稀疏角度數(shù)據(jù)的約束，其成像逆問題病態(tài)性相比傳統(tǒng)多角度CB-XLCT成像(一般24角度)的逆問題病態(tài)性更為嚴(yán)重，直接影響了傳統(tǒng)重建方法的成像性能[18]。

為了解決此問題，本文提出一種改進(jìn)多光譜策略的稀疏角CB-XLCT成像方法。首先，將經(jīng)典的多光譜成像策略應(yīng)用于稀疏角CB-XLCT成像中，構(gòu)建高維系統(tǒng)矩陣并建模逆問題；其次，考慮高維系統(tǒng)矩陣的規(guī)模對成像質(zhì)量的影響，基于譜回歸方法對系統(tǒng)矩陣進(jìn)行子空間學(xué)習(xí)，優(yōu)化成新的系統(tǒng)矩陣；再次，對新的系統(tǒng)矩陣進(jìn)行矩陣預(yù)處理并建模新的成像逆問題；最后，采用經(jīng)典的成像算法進(jìn)行逆問題的準(zhǔn)確求解。本文分別設(shè)計稀疏角重建實驗、多角度重建實驗以及魯棒性實驗，測試所提方法的有效性和魯棒性。

1 理論方法

1.1 CB-XLCT成像模型

CB-XLCT的前向模型分為3個部分：

1)錐束X射線源發(fā)出的X射線在生物組織中的傳輸，該過程可用下式描述：

(1)

其中，X(r)表示生物組織中位置r處的X射線強度，X0為X射線源的初始強度，μx為是X射線穿過的組織的衰減系數(shù)[1]。

2)X射線照射生物體內(nèi)的納米靶標(biāo)，激發(fā)其產(chǎn)生近紅外光，該過程可表述為：

S(r)=εX(r)ρ(r)

(2)

其中，S(r)表示光源，ρ為目標(biāo)濃度，ε為發(fā)光產(chǎn)額。

3)近紅外光子經(jīng)過生物組織的散射、折射和吸收等過程，傳輸?shù)缴矬w表的光子信號可被高靈敏CCD探測并采集。經(jīng)典的擴散近似模型可被用來描述組織中的光傳輸過程：

(3)

其中，D(r)=(3(μa(r)+μ′s(r)))-1為散射系數(shù)，μa(r)和μ′s(r)分別為吸收系數(shù)和散射系數(shù)。Ω和?Ω分別表示成像區(qū)域和區(qū)域邊界。Φ(r)表示光子通量，ν表示單位法向量，γ的值為依賴光折射系數(shù)的邊界不匹配因子。

基于有限元理論[7-11]，式(3)可被離散為如下的矩陣形式：

Φ=M-1·F·ε·X·ρ

(4)

其中，F(xiàn)代表光源權(quán)重矩陣，M為正定矩陣。

令A(yù)=M-1·F·ε·X，得到生物體表面的發(fā)射光光子通量流率Φ與納米發(fā)光材料三維分布ρ之間的關(guān)系,其公式為：

Φ=Αρ

(5)

其中，ρ為待求解的納米靶標(biāo)，Φ為采集的生物體表光信號，A為系統(tǒng)矩陣。

1.2 改進(jìn)多光譜策略逆問題求解

經(jīng)CCD相機采集的生物體表微弱光學(xué)信號，根據(jù)經(jīng)典的多光譜策略可對生物體表微弱光學(xué)信號進(jìn)行濾波處理，將其劃分為多個不同波長的子帶信號[19]，進(jìn)而明顯增加了生物體表的邊界測量數(shù)據(jù)。因此，引入多光譜策略將式(5)表示為：

ΦM=AMρ

(6)

其中：

(7)

通過式(7)可以看出，采用多光譜策略構(gòu)建的ΦM直接增加了邊界測量數(shù)據(jù)量，有效緩減了逆問題的病態(tài)性。然而，相應(yīng)系統(tǒng)矩陣AM明顯增加了原系統(tǒng)矩陣A的矩陣規(guī)模。為了保證成像質(zhì)量同時減少高維矩陣規(guī)模的影響，本文借鑒流形學(xué)習(xí)的思想，基于譜回歸(Spectral Regression, SR)方法對系統(tǒng)矩陣AM進(jìn)行子空間學(xué)習(xí)。本質(zhì)上，SR方法來源于一個通用的譜嵌入框架[20-21]，基于該框架可概括為如下4步來獲得AM的低維投影矩陣R：

Step1構(gòu)建權(quán)重矩陣Kij與近鄰圖G。Kij表示節(jié)點之間的相鄰關(guān)系。近鄰圖G來表示AM的內(nèi)在高維幾何結(jié)構(gòu)，將高維系統(tǒng)矩陣AM中的每一行看成是一個多元統(tǒng)計變量:

(8)

其中，節(jié)點i和j分別表示原始數(shù)據(jù)集中的行向量xi和xj。如果節(jié)點i和j屬于同一個類標(biāo)簽，將其用邊連接，同時賦以權(quán)重1/lk表示邊的權(quán)值，lk為相應(yīng)標(biāo)簽內(nèi)樣本數(shù)量[20]。

Step2根據(jù)文獻(xiàn)[22]構(gòu)建拉普拉斯矩陣L=D-K,其中D為一個基于權(quán)重矩陣Kij構(gòu)造的對角矩陣且Dii=∑jKji。

Step3利用公式(9),計算特征向量y的值。

(9)

其中，T表示矩陣轉(zhuǎn)置。

Step4計算投影矩陣R=[r1,r2,…,rc-1](c

rk=(XXT+αI)-1Xyk

(10)

其中X=[x1,x2,…,xm]T，y=[y1,y2,…,ym]T。

基于上述3步計算獲得的投影矩陣R，由于是對高維系統(tǒng)矩陣AM的子空間投影，對初始逆問題式(6)進(jìn)行優(yōu)化可得：

(11)

隨后，基于稀疏感知理論[22-23]，為了進(jìn)一步保證逆問題成像的質(zhì)量，對式(11)進(jìn)行矩陣預(yù)處理(preconditioning)[10]，隨后具體步驟如下：

(12)

2)設(shè)置預(yù)處理矩陣R：

R=(ΛΛT+λI)-1/2UT

(13)

其中，I為單位矩陣，λ為正則參數(shù)(λ>0),T表示對矩陣進(jìn)行轉(zhuǎn)置操作。

3)對公式(11)進(jìn)行矩陣預(yù)處理：

(14)

4)可獲得優(yōu)化后的新的逆問題如下所示：

(15)

其中，σ為正則化參數(shù)。對于式(15)，基于貪婪方法可采用經(jīng)典的正交匹配追蹤(Orthogonal Matching Pursuit, OMP)算法進(jìn)行求解[24]。

2 實驗與結(jié)果

2.1 實驗設(shè)置

為了驗證本文所提方法的有效性和魯棒性，借鑒文獻(xiàn)[9]的工作，分別仿真模擬了稀疏角度4投影和全角度24投影的成像對比實驗。仿真數(shù)據(jù)采用的是非勻質(zhì)仿體模型，仿體半徑為10 mm，高為30 mm。該仿體模性包含了6個器官，分別是心臟、左肺、右肺、肝臟、骨骼以及肌肉。在肺部置入1個半徑為0.5 mm、高為1 mm的光源來模擬早期腫瘤目標(biāo)，光源距離X軸方向最近表面深度為4 mm。

在具體對比中，分別引入經(jīng)典的主成分分析(Principal Components Analysis, PCA)以及不采用本文方法直接重建進(jìn)行成像對比，最終的重建算法均為OMP算法。在結(jié)果量化分析中，引入中心定位誤差(Location Error, LE)[2]、恢復(fù)濃度(Relative Quantity, RQ)[7]以及相對定量誤差(Relative Quantity Error, RQE)[7]來作為量化評價指標(biāo)。在重建中，為了避免經(jīng)驗選取正則參數(shù)帶來的主觀影響，本文采用經(jīng)典L曲線方法自適應(yīng)選取正則參數(shù)[1]。所有仿真實驗都是基于Matlab平臺，在配置為3.60 GHz Intel? Xeon? Processor i7-7820X處理器和64 GB內(nèi)存的臺式工作站上進(jìn)行計算。

2.2 稀疏角度與全角度重建實驗對比

在前向仿真中，本節(jié)分別模擬稀疏角度4次X射線激發(fā)與全角度24次X射線激發(fā)，采用XLCT領(lǐng)域中5種常見的子帶波長(585 nm, 627 nm, 667 nm, 720 nm 和 802 nm)進(jìn)行逆向重建[1-2,7-10,25-29]，不同子帶波長下的各個器官的光學(xué)參數(shù)根據(jù)Alexandrakis等人[30]的工作進(jìn)行計算。仿真實驗納米發(fā)光材料模擬NaGdF4: Eu3，相應(yīng)的X射線的衰減系數(shù)為0.054 mm-1[31]。目標(biāo)的納米發(fā)光產(chǎn)額設(shè)置為0.15 cm3/mg[7]。納米靶標(biāo)顆粒濃度設(shè)置0.3183μg/mm3[7]。經(jīng)過前向仿真計算，對于單光源數(shù)字仿體實驗的5種子帶波長(585 nm, 627 nm, 667 nm, 720 nm 和 802 nm)對應(yīng)各自的能量權(quán)重系數(shù)分別為0.06, 0.15, 0.22, 0.26和0.31[29]。圖1為采用的數(shù)字仿體以及2種激發(fā)方案下各自前向仿真結(jié)果展示。基于有限元方法，前向仿真采用的仿體被離散為167663個四面體單元和29597個結(jié)點。

圖1 實驗采用的仿體以及2種激發(fā)方案下各自前向仿真結(jié)果展示

在逆向多光譜重建中，分別采用本文所提方法(SR+preconditioning)、結(jié)合PCA的方法(PCA+preconditioning)以及直接重建(with no processing)進(jìn)行對比，采用的仿體包含47873個四面體單元和8674個結(jié)點。圖2和圖3分別是稀疏4角度投影和全24角度投影各自的重建結(jié)果展示。表1和表2分別為相應(yīng)的量化重建結(jié)果。

表1 稀疏4角度投影下的量化重建結(jié)果

表2 全24角度投影下的量化重建結(jié)果

圖2 稀疏4角度投影下，重建結(jié)果在Z=15 mm平面處的二維截面圖與重建結(jié)果的三維展示圖

通過上述結(jié)果可以看出，對于全24角度投影重建實驗，所有方法的重建誤差LE均在1 mm以內(nèi)，恢復(fù)濃度均在0.265 μg/mm3～0.295 μg/mm3，相應(yīng)的相對量化誤差RQE均在20%以內(nèi)。對于稀疏4角度投影重建實驗，由于測量數(shù)據(jù)的明顯不足，使得逆問題病態(tài)性明顯加劇，直接重建的誤差LE已大于1 mm，恢復(fù)濃度已降至0.235 μg/mm3，相應(yīng)的相對量化誤差RQE已增至26.26%。而對于本文提出的SR+preconditioning方法與PCA+preconditioning方法，重建的誤差LE仍然在1 mm以內(nèi)，恢復(fù)濃度相對于全24角度投影重建結(jié)果略有降低，但相應(yīng)的相對量化誤差RQE仍在20%以內(nèi)。特別地，不論是重建誤差LE還是恢復(fù)濃度RQ及相應(yīng)的相對量化誤差RQE，本文所提方法的成像質(zhì)量最好。

2.3 魯棒性測試

在上一節(jié)重建實驗結(jié)果的基礎(chǔ)上，本節(jié)針對稀疏4角度投影，進(jìn)一步設(shè)計了2組魯棒性實驗?？紤]影響成像質(zhì)量常見的因素[6-7,19]，分別設(shè)計了靶標(biāo)深度水平測試與噪聲水平測試。其中，靶標(biāo)深度水平測試中設(shè)置光源距X軸方向最近表面深度分別為8 mm、7 mm、6 mm、 5 mm和4 mm，噪聲水平測試中分別設(shè)置噪聲為高斯白噪聲，該噪聲服從均值為0的正態(tài)分布，方差以表面測量值均值的5%為間隔在10%到30%之間變換。實驗結(jié)果如圖4和表3、表4所示。

圖4 稀疏4角度投影下深度水平測試與噪聲水平測試中所提方法各自的重建結(jié)果展示

表3 稀疏4角度投影下深度水平測試中所提方法量化重建結(jié)果

表4 稀疏4角度投影下噪聲水平測試中所提方法量化重建結(jié)果

通過重建結(jié)果可以看出，盡管目標(biāo)深度的增加或者噪聲水平的增加都會不同程度增加逆問題重建的難度，進(jìn)而影響成像質(zhì)量。但本文所提方法在上述因素的測試中，重建誤差LE均在1 mm以內(nèi)。對于深度水平測試，隨著目標(biāo)深度的減小，恢復(fù)濃度有所改善，相應(yīng)的相對量化誤差也逐漸減小。對于噪聲水平測試，重建誤差LE在噪聲水平的變化下整體較穩(wěn)定，而恢復(fù)濃度RQ及相應(yīng)的相對量化誤差RQE相對于重建誤差LE更易受噪聲的干擾。其中，當(dāng)噪聲水平設(shè)置為30%時，相對量化誤差RQE達(dá)到了30.62%?？傮w上，本文所提方法對于深度因素和噪聲因素的影響，成像質(zhì)量是可接受的。

3 結(jié)束語

稀疏角CB-XLCT成像作為一種新型的光學(xué)斷層成像技術(shù)，具有對早期腫瘤進(jìn)行實時監(jiān)測的應(yīng)用潛力。但受限于測量數(shù)據(jù)的嚴(yán)重不足，相對于傳統(tǒng)多角度成像，逆問題病態(tài)性明顯加劇。本文提出了一種改進(jìn)多光譜的稀疏角CB-XLCT成像方法，將譜回歸與矩陣預(yù)處理先后對多光譜策略生成的高維系統(tǒng)矩陣進(jìn)行子空間學(xué)習(xí)以及進(jìn)一步優(yōu)化重建逆問題，最終采用經(jīng)典OMP算法對優(yōu)化后的逆問題進(jìn)行準(zhǔn)確求解。本文設(shè)計了稀疏4角度投影仿真實驗、全24角度投影仿真實驗、深度水平測試以及噪聲水平測試。實驗結(jié)果表明，本文所提方法不僅可以準(zhǔn)確定位目標(biāo)和恢復(fù)濃度，同時對于常見干擾因素也具有較好的魯棒性。

總之，作為醫(yī)學(xué)分子影像的研究熱點，CB-XLCT應(yīng)用于早期腫瘤的實時監(jiān)測還有很多問題需要解決，比如在成像中減少投影數(shù)量進(jìn)行單視圖重建以及在重建中對腫瘤形狀進(jìn)行準(zhǔn)確勾勒；此外，對于核譜回歸方法在CB-XLCT實時成像中的進(jìn)一步探索也是筆者未來的研究重點。