劉艷秋，胡先功，張 衡，郭紅波，賀小偉，3*

（1.西北大學(xué) 信息科學(xué)與技術(shù)學(xué)院，陜西 西安 710127；2.西北大學(xué) 西安市影像組學(xué)與智能感知重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710127；3.西北大學(xué) 網(wǎng)絡(luò)和數(shù)據(jù)中心，陜西 西安 710127）

1 引 言

生物發(fā)光成像（Bioluminescence Imaging，BLI）是一種非侵入性、非電離的分子成像技術(shù)，能夠在活體中顯示熒光素標(biāo)記的腫瘤細(xì)胞。由于腫瘤細(xì)胞具有較高的敏感性和特異性，BLI在小動(dòng)物的預(yù)臨床研究中發(fā)揮著越來越重要的作用[1-3］。與BLI相比，生物發(fā)光斷層成像（Bioluminescence Tomography，BLT）是采用生物發(fā)光源三維重建，能更準(zhǔn)確地定位和量化腫瘤信息[4］，已廣泛應(yīng)用于預(yù)臨床研究。然而在實(shí)際生物應(yīng)用中，由于光子在生物組織中的復(fù)雜傳播和探測(cè)到的表面光子的限制，BLT的病態(tài)性限制了其光源重建精度[5-6］。

在生物發(fā)光成像中，光從生物體內(nèi)部光源發(fā)出，通過復(fù)雜生物組織結(jié)構(gòu)到達(dá)生物體表面。描述光在生物體內(nèi)的傳輸過程的一個(gè)標(biāo)準(zhǔn)模型為輻射傳輸方程（Radiative Transfer Equation，RTE）[7］。然而，在實(shí)際的生物發(fā)光成像中，RTE的計(jì)算成本很高[8］。最典型的方法是使用RTE的擴(kuò)散近似方程（Diffusion approximation Equation，DE）作為前向模型。DE基本上是RTE的一階球諧波近似的特例，是描述光在組織中傳播速度最快的模型，但其精度有限，特別是在低散射區(qū)域、高吸收區(qū)域和邊界處[9］。為了克服DE模型的局限性，需要對(duì)RTE進(jìn)行高階近似。相比于其它的高階近似模型（如SN、PN模型），簡化球諧近似（Simplified spherical harmonic approximation，SPN）的計(jì)算負(fù)載更低，三階簡化球諧近似（3rd Simplified spherical harmonic approximation equation，SP3）已被廣泛應(yīng)用于生物發(fā)光成像[10-12］。但由于SP3也是對(duì)RTE的近似，模型誤差仍然存在，而且高階近似需要求解的未知量是DE的幾倍，帶來了更高的計(jì)算負(fù)載。混合SP3和DE光傳輸模型（Hybrid Simplified spherical harmonics with Diffusion Equation，HSDE）的提出，一定程度上平衡了光傳輸精度和速度[13-14］，然而簡化模型對(duì)于RTE的近似誤差仍然存在，SP3模型導(dǎo)致的計(jì)算負(fù)載也不能完全得到緩解。

另外，由于將成像物體表面檢測(cè)到的二維發(fā)光能量轉(zhuǎn)化為物體內(nèi)重建的三維生物發(fā)光源的BLT逆問題是一個(gè)數(shù)學(xué)不適定問題[4］。因此，需要融合足夠的先驗(yàn)知識(shí)來保證解的唯一性和穩(wěn)定性。可行域策略和多光譜方法分別通過減少未知變量和增加已知數(shù)據(jù)量減少光源重建的病態(tài)性，是BLT重建中常用的兩種策略[15］。此外，基于壓縮感知（Compressed Sensing，CS）理論的各種稀疏算法，如采用稀疏正則化項(xiàng)[16-17］（如L0和L1），使用貪婪策略[18］或采用稀疏貝葉斯方法[5，19］等算法，用來減少噪聲對(duì)重建結(jié)果的影響，使光學(xué)分子斷層成像的重建精度不斷提高。

在我們之前的研究中，分析了RTE的擴(kuò)散近似產(chǎn)生的光傳輸誤差、數(shù)據(jù)采集過程中成像系統(tǒng)存在的噪聲以及逆問題固有的病態(tài)性對(duì)BLT光源重建精度的影響，提出了光譜差分策略，減弱DE光傳輸模型和成像系統(tǒng)的系統(tǒng)誤差，并緩解重建問題的病態(tài)性[20］。在此基礎(chǔ)上，本文將光譜差分策略分別應(yīng)用到DE和SP3光傳輸模型，首先對(duì)這兩種RTE近似產(chǎn)生的誤差進(jìn)行分析，對(duì)比光譜差分策略對(duì)兩種光傳輸誤差的衰減作用，進(jìn)而將光譜差分策略分別應(yīng)用到DE和SP3光傳輸模型進(jìn)行光源重建，旨在有效提高光源重建的準(zhǔn)確性和效率。

2 理論方法

2.1 光傳輸過程及誤差分析

穩(wěn)態(tài)的RTE為：

（1）DE光傳輸模型誤差

穩(wěn)態(tài)擴(kuò)散方程表示為：

通過有限元近似表示為：

A(DE)為DE近似構(gòu)造的系統(tǒng)矩陣，S為光源能量密度，Φm為表面節(jié)點(diǎn)上的光通量。

DE是RTE的一階球諧展開近似，在我們之前的研究中，分析了該近似過程對(duì)給定波長λ對(duì)應(yīng)的系統(tǒng)矩陣A(DE)λ產(chǎn)生的誤差[20］，現(xiàn)將其簡單表示為：

其中：受波長λ影響的量用Xλ表示，其它與波長無關(guān)的量用Y表示。

對(duì)于兩個(gè)波長λj和λk對(duì)應(yīng)的誤差做差：

（2）SP3光傳輸模型誤差

高階簡化球諧近似方程可以表示為：

其中：μa，n=μa+μs(1-gn)，g為各向異性參數(shù)，n為階數(shù)；φ為勒讓德多項(xiàng)式。N=3時(shí)，SP3模型的方程為：

φi是輻射度的勒讓德矩φi的線性組合。通過有限元近似表示為：

其中，A(SP3)為SP3近似構(gòu)造的系統(tǒng)矩陣，類似于DE近似。SP3近似對(duì)給定波長λ對(duì)應(yīng)的系統(tǒng)矩陣A(SP3)λ產(chǎn)生的誤差可以表示為：

其中：受波長λ影響的量用Xλ表示，其它與波長無關(guān)的量用Y表示。

對(duì)于兩個(gè)波長λj和λk之間的差值：

表明對(duì)于DE和SP3光傳輸模型，由于波長接近的光在通過相同光學(xué)路徑時(shí)遇到相近的系統(tǒng)響應(yīng)，對(duì)相近波長的測(cè)量數(shù)據(jù)做差可以一定程度上消除光傳輸模型產(chǎn)生的誤差。

2.2 光譜差分策略在光源重建中的應(yīng)用

由于成像問題的非唯一性以及光源重建的病態(tài)性，已證實(shí)多光譜方法可以緩解其病態(tài)性[21，22］。假設(shè)測(cè)量數(shù)據(jù) 包含3組波長，對(duì)于DE和SP3光傳輸模型，將各自系統(tǒng)矩陣組合：

其中：Aλi是給定波長λi對(duì)應(yīng)的系統(tǒng)矩陣，Φmλi是 相同波長的表面光能量，Amulti和Φmmulti是組合的系統(tǒng)矩陣和表面光能量。

在多光譜方法的基礎(chǔ)上，光譜差分策略利用每個(gè)測(cè)量波長之間的數(shù)據(jù)差，將式（13）和（14）變換為：

其中，Adiffe和Φmdiffe分別為差分的系統(tǒng)矩陣和差分的表面光能量。

由于光源在生物體中通常是稀疏分布的，且表面測(cè)量數(shù)據(jù)嚴(yán)重不足，基于壓縮感知理論，采用Lp正則化方法將式（16）的重建模型轉(zhuǎn)化為以下最小化問題：

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

3.1 實(shí)驗(yàn)設(shè)置

采用數(shù)值仿真實(shí)驗(yàn)驗(yàn)證光譜差分對(duì)于DE和SP3光傳輸模型誤差的緩解作用，以及光譜差分策略在BLT光源重建中的性能。前向仿真實(shí)驗(yàn)和光源重建實(shí)驗(yàn)均采用常用的數(shù)字鼠模型，只研究數(shù)字鼠的軀干部分，高度為35 mm，包括肌肉、心臟、胃、肝臟、腎臟和肺六個(gè)主要組織，如圖1（a）所示。在肝臟中設(shè)置一個(gè)半徑為1 mm的球形光源，其中心點(diǎn)坐標(biāo)為（18 mm，8 mm，14.8 mm），位置如圖1（b）。在前向仿真實(shí)驗(yàn)中，三維數(shù)字鼠軀干離散網(wǎng)格如圖1（c）所示，包含20 263個(gè)節(jié)點(diǎn)。光源重建實(shí)驗(yàn)網(wǎng)格如圖1（d）所示，包含10 139個(gè)節(jié)點(diǎn)。本文選用BLT中常用波段范圍的三個(gè)波長630 nm、650 nm和670 nm進(jìn)行實(shí)驗(yàn)，由文獻(xiàn)[25］中總結(jié)的公式，計(jì)算得到各生物組織在各個(gè)波長下對(duì)應(yīng)的光學(xué)參數(shù)如表1所示。

在前向仿真實(shí)驗(yàn)中，以蒙特卡羅方法（Monte Carlo extreme，MC）[26］獲得的表面光能量作為對(duì)比標(biāo)準(zhǔn)。以平均誤差（Average error）和余弦相似度（Cosine similarity）作為定量評(píng)價(jià)指標(biāo)，分析DE和SP3引起的光傳輸模型誤差。

測(cè)量的表面光能量差值的計(jì)算方法如下：

計(jì)算630 nm、650 nm和670 nm波長下，分別通過MC和DE（SP3）獲得的表面光能量的差異，以及MC和DE（SP3）獲得的不同波長的表面光能量任意兩組波長組合，差分之后的能量誤差。進(jìn)而采用平均誤差定量評(píng)價(jià)上述能量誤差：

其中：EEnengydifferencei代表第i個(gè)表面節(jié)點(diǎn)對(duì)應(yīng)的能量誤差，N為表面節(jié)點(diǎn)總數(shù)目。

圖1 數(shù)值仿真實(shí)驗(yàn)設(shè)置Fig.1 Numerical simulation experiment settings

表1 不同波長下數(shù)字鼠組織的光學(xué)參數(shù)Tab.1 Optical parameters of the mouse tissues at different wavelengths.

同時(shí)采用余弦相似度評(píng)估表面光能量的差異以及差分?jǐn)?shù)據(jù)的差異：

其中：Ai、Bi分別表示被比較的向量A、B的分量。其值越接近1，表明兩個(gè)向量越接近。

在光源重建實(shí)驗(yàn)中，采用常用的定量分析指標(biāo)：定 位 誤 差（Location Error，ILE）、Dice系 數(shù)（IDice）和對(duì)比噪聲比（Contrast-to-Noise Ratio，ICNR），定量分析各方法的目標(biāo)定位、形狀恢復(fù)和圖像對(duì)比度性能。這些指標(biāo)的計(jì)算方法如下：

其中：ILE表示重建光源與真實(shí)光源之間的位置偏差，（x，y，z）和（x0，y0，z0）分別為重建能量加權(quán)中心點(diǎn)的坐標(biāo)和真實(shí)光源中心的坐標(biāo)。

Dice系數(shù)用來評(píng)價(jià)形狀恢復(fù)，其中，X、Y分別為重建光源區(qū)域和實(shí)際光源區(qū)域。

對(duì)比噪聲比用來評(píng)價(jià)圖像對(duì)比度，其中：下標(biāo)ROI和BCK分別表示被成像物體的目標(biāo)區(qū)域和背景區(qū)域，μ、ω、σ分別表示平均強(qiáng)度值、加權(quán)因子和方差。

3.2 前向仿真實(shí)驗(yàn)結(jié)果

為了分析比較DE和SP3模型在光傳輸過程中產(chǎn)生的誤差，進(jìn)行了前向仿真實(shí)驗(yàn)，進(jìn)而驗(yàn)證光譜差分對(duì)于光傳輸誤差消除的有效性。

首先，采用MC方法獲得數(shù)字鼠模型對(duì)應(yīng)630 nm、650 nm和670 nm波長的表面光能量，為了便于展示光能量的細(xì)節(jié)，將其投影到X-Z平面，如圖2（a）所示。相應(yīng)的，分別采用DE和SP3模型獲得數(shù)字鼠模型在各波長的表面光能量，分別如圖2（b）、（c）所示。對(duì)于各個(gè)光傳輸模型獲得的表面光能量，都表現(xiàn)為波長越長，對(duì)應(yīng)的表面光能量越強(qiáng)。這是因?yàn)樵谳^長波長下，組織散射效應(yīng)減弱，光穿透能力增強(qiáng)。對(duì)比圖2（a）、（b）和（c），直觀上，DE模型在630 nm獲得的表面光能量和MC與SP3差異較明顯，而在650 nm和670 nm對(duì)應(yīng)的表面光能量較為接近。

圖2 MC、DE和SP3模型獲得的表面光能量在X-Z平面投影Fig.2 X-Z plane projection of surface luminescence fluxes obtained by the MC，DE and SP3，respectively.

為了定量分析DE和SP3模型的光傳輸誤差，采用公式（18）計(jì)算各特定光譜MC和DE（SP3）獲得的表面光能量的誤差，再利用公式（19）計(jì)算平均誤差。進(jìn)而應(yīng)用光譜差分，任意兩組波長組合，利用式（18）計(jì)算各組測(cè)量波長之間的能量差異，獲得對(duì)應(yīng)的差分后的能量誤差，再利用式（19）計(jì)算出平均誤差。DE和SP3模型以及光譜差分后DE和SP3對(duì)應(yīng)的平均誤差值如圖3所示，對(duì)于DE模型，630 nm對(duì)應(yīng)的能量平均誤差較大，650 nm、670 nm對(duì)應(yīng)的能量誤差相應(yīng)減小。SP3模型在各波長對(duì)應(yīng)的能量平均誤差和DE模型呈現(xiàn)相同的變化趨勢(shì)，且各能量平均誤差值都相對(duì)DE模型減小，表明SP3模型比DE模型具有更高的傳輸精度。分別對(duì)DE模型和SP3模型獲得的表面光能量進(jìn)行光譜差分，能量平均誤差相對(duì)減小，尤其對(duì)于630 nm對(duì)應(yīng)的能量平均誤差改善較為明顯，而且應(yīng)用光譜差分策略后SP3模型比DE模型獲得了更小的能量平均誤差。

圖3 MC和DE（SP3）獲得的表面光能量的平均誤差，以及MC和DE（SP3）獲得的表面光能量通過光譜差分后的平均誤差Fig.3 Average errors of surface light energy obtained by MC and DE（SP3）and the average errors after applying the spectral differential

為了進(jìn)一步驗(yàn)證前向仿真實(shí)驗(yàn)的結(jié)果，由公式（20）計(jì)算了MC和DE（SP3）在同一波長下獲得的表面光能量的余弦相似度，如圖4所示，余弦相似度越接近1，表明能量值越接近。由圖4可以看出，SP3模型的余弦相似度大于DE模型。進(jìn)而通過MC與DE（SP3）獲得的各波長的能量差計(jì)算余弦相似度，從圖4可以看出，光譜差分后DE和SP3對(duì)應(yīng)的余弦相似度顯著增加。計(jì)算各個(gè)特定光譜和光譜差分對(duì)應(yīng)的余弦相似度的平均值和方差，DE和SP3模型對(duì)應(yīng)的三組波長余弦相似度的平均值和方差分別為0.985 8±0.006 0和0.988 1±0.003 0，任意兩組波長組合光譜差分后DE和SP3模型對(duì)應(yīng)的余弦相似度的平均值和方差分別為0.989 3±6.127 6e-04和0.990 8±1.042 7e-04。表明光譜差分能夠有效減小DE模型的傳輸誤差，使其接近SP3模型的傳輸精度，并且對(duì)SP3模型應(yīng)用光譜差分能夠獲得更高的傳輸精度。

圖4 MC和DE（SP3）獲得的表面光能量的余弦相似度，以及通過MC與DE（SP3）獲得的各波長的能量差計(jì)算的余弦相似度Fig.4 Cosine similarity of surface light energy obtained by MC and DE（SP3）and the cosine similarity calculated by the difference of energy between each measured wavelength obtained by MC and DE（SP3）

進(jìn)而對(duì)各方法的運(yùn)算效率進(jìn)行論證，在獲得各波長DE和SP3模型對(duì)應(yīng)的表面光能量的同時(shí)記錄下了其所耗時(shí)間，DE模型在630 nm、650 nm和670 nm波長下耗時(shí)分別為：212.90 s、208.41 s和211.73 s，而SP3模型在三個(gè)波長耗時(shí)分別為：992.55 s、990.35s和994.26 s，DE模型耗時(shí)遠(yuǎn)低于SP3。對(duì)三組波長對(duì)應(yīng)的表面光能量進(jìn)行光譜差分，DE模型的光能量差分總耗時(shí)為0.031 s，SP3模型的光能量差分總耗時(shí)為0.030 s。即采用DE模型獲得三組波長的光能量并進(jìn)行光譜差分，總耗時(shí)仍低于SP3模型獲得一組波長光能量的耗時(shí)。

3.3 光源重建實(shí)驗(yàn)結(jié)果

為了驗(yàn)證光譜差分策略在光源重建中的可行性與適用性，進(jìn)行逆向光源重建實(shí)驗(yàn)。采用MC方法模擬 的630 nm、650 nm和670 nm波長的表面光能量作為測(cè)量的表面數(shù)據(jù)。分別對(duì)DE和SP3模型獲得的三組波長對(duì)應(yīng)的系統(tǒng)矩陣進(jìn)行光譜差分，獲得式（16）的表達(dá)，采用ISTA方法重建，并將MC獲得的630 nm波長的表面光能量分別采用DE和SP3模型進(jìn)行光源重建作為對(duì)比。

圖5 光源重建實(shí)驗(yàn)結(jié)果Fig.5 Experimental results of source reconstruction

DE和SP3模型在630 nm對(duì)應(yīng)的重建結(jié)果分別如圖5（a）和（b）所示，采用光譜差分策略的DE（DE-spectral differential）、SP3方 法（SP3-spectral differential）對(duì)應(yīng)的重建結(jié)果分別如圖5（c）和（d）所示。三維視圖中的紅色球體表示真實(shí)光源的實(shí)際位置，綠色不規(guī)則形狀為重建光源。相應(yīng)的矢狀面、冠狀面和橫切面根據(jù)真實(shí)光源的中心位置來確定，截面圖中的紅色圓圈代表真實(shí)光源的位置。圖5（a）、（b）中重建光源距離真實(shí)光源有一定的偏離，且由于重建光源位置比真實(shí)光源偏上，在Z=14.8 mm橫切面圖中觀測(cè)不到重建光源信息。圖5（c）、（d）中DE-spectral differential和SP3-spectral differential方法都能夠較為準(zhǔn)確地定位真實(shí)光源，然而，相比于DE-spectral differential方法，SP3-spectral differential方法 的重建 光源與真實(shí)光源重疊更好。

表2中的評(píng)價(jià)指標(biāo)定量分析驗(yàn)證了圖5中的結(jié)論，DE和SP3模型在630 nm對(duì)應(yīng)的重建結(jié)果LE均大于1 mm，且Dice較小。這和前向仿真實(shí)驗(yàn)中，在630 nm波長DE和SP3模型獲得的表面光能量和MC的誤差較明顯的結(jié)論一致。DEspectral differential方法獲得的LE小于1 mm，且Dice和CNR值相對(duì)增大，提高了光源重建的準(zhǔn)確性。而SP3-spectral differential方法重建光源的LE最?。?.724 mm），Dice和CNR值最大，在目標(biāo)定位、形狀恢復(fù)和圖像對(duì)比度等方面有更好的效果。這是由于光譜差分策略在光源重建過程中，有效減少了DE和SP3模型的光傳輸誤差，提高了模型精度，并且緩解了BLT逆問題的病態(tài)性，使光源重建的位置誤差小于1 mm，獲得了更準(zhǔn)確的重建精度。

表2 光源重建實(shí)驗(yàn)的定量結(jié)果Tab.2 Quantitative results in source reconstruction experiment

此 外，DE模型 獲 得630 nm、650 nm和670 nm波長對(duì)應(yīng)的系統(tǒng)矩陣耗時(shí)分別為33.61 s、33.66 s和34.81 s，SP3模型獲得三組波長對(duì)應(yīng)的系統(tǒng)矩陣耗時(shí)分別為：1 530.68 s、1 520.63 s和1 522.52 s，DE模型構(gòu)建系統(tǒng)矩陣耗時(shí)遠(yuǎn)低于SP3。對(duì)三組波長對(duì)應(yīng)的系統(tǒng)矩陣以及表面光能量進(jìn)行光譜差分，獲得式（16）中差分的系統(tǒng)矩陣和差分的表面光能量，DE模型差分總耗時(shí)為0.401 s，SP3模型差分總耗時(shí)為0.414 s。即采用DE模型獲得三組波長的系統(tǒng)矩陣，并對(duì)系統(tǒng)矩陣和表面光能量進(jìn)行光譜差分，總耗時(shí)仍低于SP3模型獲得一組波長數(shù)據(jù)的耗時(shí)。

4 結(jié) 論

在BLT研究中，高階光傳輸模型可提高光在生物組織中的傳輸精度，多光譜方法可降低逆問題病態(tài)性，從不同角度提升光源信息的重建精度。在此基礎(chǔ)上，本文對(duì)DE和SP3光傳輸模型獲得的表面光能量進(jìn)行光譜差分，對(duì)比了光譜差分策略對(duì)兩種光傳輸誤差的衰減作用，前向仿真實(shí)驗(yàn)結(jié)果表明光譜差分策略能有效地減少DE和SP3的模型誤差，對(duì)DE模型采用光譜差分，能夠獲得接近SP3模型的傳輸精度，并且降低高階近似對(duì)運(yùn)算時(shí)間和存儲(chǔ)空間的高要求。進(jìn)而將光譜差分策略分別應(yīng)用到DE和SP3光傳輸模型進(jìn)行光源重建，實(shí)驗(yàn)結(jié)果表明光譜差分策略在提高兩種光傳輸模型精度的同時(shí)，緩解了BLT中逆問題的病態(tài)性，使光源重建在目標(biāo)定位、形狀恢復(fù)和圖像對(duì)比度等方面具有更準(zhǔn)確的效果。其中DE模型結(jié)合光譜差分策略，在重建精度和速率上達(dá)到了較好的平衡。