快速熒光擴散斷層成像的格子玻爾茲曼前向模型

岑星星，嚴壯志,

1 上海大學(xué) 通信與信息工程學(xué)院，上海市，200444

2 上海大學(xué) 生物醫(yī)學(xué)工程研究所，上海市，200444

0 引言

熒光擴散斷層成像（Fluorescence Diffuse Optical Tomography,FDOT）是一種新型的醫(yī)學(xué)成像模態(tài)，通過體表采集的熒光信號，完成體內(nèi)分子探針的三維成像功能，能在細胞、分子水平實現(xiàn)相應(yīng)病灶的定位及監(jiān)測工作[1]。FDOT所具有的高靈敏度、安全無創(chuàng)、低成本等優(yōu)勢，使得該成像技術(shù)在近幾年得到迅速發(fā)展，并在臨床、預(yù)臨床等醫(yī)學(xué)領(lǐng)域取得廣泛應(yīng)用，如：癌癥檢查、藥物研究等[2]。對于FDOT的發(fā)展，關(guān)于前向模型和重建算法的研究必不可少，前向模型描述了光在特定介質(zhì)中的傳播過程；重建算法則通過前向模型和邊界采集的光學(xué)信號，對體內(nèi)標記物的光場分布進行重建計算[3]。由于生物組織的高散射性和探測器的測量不充分，相應(yīng)的重建計算屬于不適定問題，使得所用前向模型的微小偏離會造成重建結(jié)果的巨大誤差。因此，建立一個準確且計算可行的前向模型對于FDOT的實現(xiàn)非常重要。

輻射傳輸方程（Radiative Transfer Equation,RTE）和擴散方程（Diffusion Equation，DE）是FDOT中的兩個經(jīng)典前向模型。RTE是一種被廣泛接受的高精度模型，但受制于昂貴的計算成本，使其僅出現(xiàn)于簡單的FDOT運用場景中[4]。相較于RTE，RTE的簡化模型有著更好的適用性，其中的DE作為RTE的一階球諧近似，具有更高的計算效率，并在FDOT的相關(guān)研究中有著非常廣泛的應(yīng)用。然而，DE也存在著一系列的問題，無法對靠近光源和邊界區(qū)域的光傳播過程進行準確描述[5]。此外，DE經(jīng)有限元等數(shù)值方法求解的計算效率依然無法滿足我們的需求，尤其是在增加計算網(wǎng)格、探測器數(shù)量以保證FDOT成像質(zhì)量的情況下。

格子玻爾茲曼（Lattice Boltzmann,LB）方法具有物理含義清晰、計算程序簡單、計算結(jié)構(gòu)并行性高等優(yōu)點[6]，開發(fā)基于LB的前向模型，對于FDOT的性能提升有著巨大潛能。LB方法源于格子氣自動機（LGA）方法，是一種特殊的離散分子動力學(xué)模型，已成功應(yīng)用于各種物理學(xué)研究，包括流體力學(xué)、熱力學(xué)等[7]。近年來，LB方法也被發(fā)展用于輻射傳輸描述，涉及源項、介質(zhì)和邊界上的輻射傳輸描述[8-10]。研究表明，用于輻射描述的LB能提供精確的計算結(jié)果，且有著較高的計算效率。然而，當前所涉及的大多數(shù)LB輻射傳輸研究都是處于封閉長方體容器中的熱輻射問題，極少有完整關(guān)于FDOT前向模型的研究報道。

基于我們最初的想法[11]，本文提出了基于LB的FDOT前向模型，通過RTE離散向原有LB方程引入光傳播描述[12-13]。進一步，本文給還出了基于LB前向模型的FDOT算法流程，選擇代數(shù)重建技術(shù)（Algebra Reconstruction Technique，ART）作為重建算法，尋找成像標記物的相應(yīng)光源分布，將其代入LB前向模型，使計算得到的模擬測量能與實際的測量達到最佳擬合。

1 基于格子玻爾茲曼方法的前向模型

1.1 光在生物組織內(nèi)的傳播模型

LB前向模型的光傳播描述，通過RTE引入。類比于傳統(tǒng)LB模型的離散化推導(dǎo)，得到經(jīng)RTE離散后的LB前向模型。

傳統(tǒng)的LB基本模型可由玻爾茲曼方程經(jīng)時間、速度、相空間的離散化得到[13]。玻爾茲曼方程描述了粒子在介觀尺度上的統(tǒng)計行為，相應(yīng)的Bhatnagar-Gross-Krook（BGK）近似如下所示[14]：

其中，f=f(r,ξ,t)為分子密度函數(shù)，代表t時刻，位置r處，速度為ξ的分子數(shù)量；λ為松弛時間，代表f趨近于平衡態(tài)feq的速度快慢。假設(shè)離散的時間間隔為△t，離散后的空間為Ωk，離散速度模型為DmQn（m代表空間維數(shù)，n代表離散速度的數(shù)量），經(jīng)玻爾茲曼BGK離散得到的LB基本模型如下所示[13]：

其中，f(r,ξ,t)為離散后的分子密度函數(shù)，代表t時刻，位置r∈Ωk處，以離散速度ξi進行運動的分子數(shù)量（下標i為離散后速度方向的索引，介于1到n之間）；f(r+ξ△t,ξ,t+△t)為下一時刻t+△t上的分子密度函數(shù)；τ≡λ/△t為無量綱的松弛因子；feq(r,ξi,t)為離散后的平衡態(tài)分布函數(shù)。

同時，RTE可以改寫為玻爾茲曼方程的類似形式，公式及物理描述的相似性提供了條件支持。對于物理空間Ω中的單色光傳播，RTE可給出精確描述[4]。進一步，RTR可重寫為如下形式：

其中，w(r)=μs(r)/β(r)代表反照率；為時刻沿方向運動的光源項；為散射項函數(shù)，描述了光子經(jīng)散射后由方向運動到方向的可能性。對比式（1）和式（3）可以看出，RTE和玻爾茲曼BGK具有相似的數(shù)學(xué)形式，并且前者的對應(yīng)后者的平衡態(tài)分布函數(shù)feq，cβ(r)對應(yīng)松弛時間λ的倒數(shù)。此外，對兩者描述的物理量進行觀察，發(fā)現(xiàn)RTE的輻射率φ和玻爾茲曼BGK的分子密度函數(shù)f，都是關(guān)于時間、物理空間和速度的函數(shù)。雖然，輻射率φ從能量角度進行描述，不同于分子密度函數(shù)f中的粒子數(shù)量，但光傳播過程中的輻射率φ和相應(yīng)的光子密度分布函數(shù)有著直接聯(lián)系[12]：

其中，h為普朗克常量；v為光的頻率大小。

通過上述過程建立了RTE和玻爾茲曼BGK的聯(lián)系。接下來本文類比于LB基本模型的離散化推導(dǎo)，通過RTE離散得到相應(yīng)的LB前向模型。對應(yīng)式（3）、式（4），RTE經(jīng)時間、速度和相空間離散后的LB形式如下所示：

其中，Θ(si,sj)為離散后的散射相函數(shù)；wi為權(quán)重系數(shù)。下標i、j指向不同的離散速度方向，當取D3Q6離散速度模型（圖1）時，為對應(yīng)箭頭方向，且介于1到6之間。

圖1 LB方法中D3Q6離散模型Fig.1 The D3Q6 for discrete model in the LB

1.2 格子玻爾茲曼前向模型的建立

1.2.1 平衡態(tài)權(quán)重系數(shù)權(quán)重系數(shù)決定了LB前向模型中平衡分布函數(shù)S(r,t)的具體表達，本文通過局部輻射場的物理守衡進行權(quán)重系數(shù)求解。假設(shè)輻射率在達到平衡態(tài)時保持各向同性，且大小為(r,t)，同時假設(shè)光在均勻介質(zhì)中以恒定速度大小c進行傳播，k階輻射強度的矩方程經(jīng)過高斯求積后可近似為：

對式（8）～（9）進行整理，可以得到如下矩陣方程：

選擇DmQn并確定離散方向后，相應(yīng)的權(quán)重系數(shù)通過式（10）求解得到，不同離散模型對應(yīng)的權(quán)重系數(shù)如下所示：

1.2.2 平衡態(tài)的散射相函數(shù)

1.2.3 邊界條件

為了能讓LB前向模型的邊界計算變得簡便，同時又能維持較好的精度，在邊界節(jié)點r∈?Ωk（k為離散空間節(jié)點的索引）進行如下光子密度φ(r,t)的計算：

并與DE中的Robin邊界條件進行關(guān)聯(lián)，得到如下邊界計算公式[18]：

其中，D(r)1/[3(μa(r)+(1-g)μs(r))]為光擴散系數(shù)；是邊界處的外法線向量；當組織周圍填充空氣，組織折射率為n時，A(r,n,n')≈[1+R(r)]/[1-R(r)]，其中的R(r)≈-1.439 9n-2+0.709 9n-1+0.668 1+0.063 6n。

2 熒光分子斷層成像實現(xiàn)

2.1 LB前向模型計算流程

本文使用迭代算法完成LB前向模型的求解。假設(shè)模擬光源為穩(wěn)定光源，通過迭代計算直到輸出結(jié)果穩(wěn)定為止。LB前向模型的具體計算流程如下：

（1）選擇離散速度模型DmQn，并用適當數(shù)量的格子對連續(xù)空間Ω進行網(wǎng)格劃分；

（2）設(shè)置初始參數(shù)μa(r)、μs(r)、g和q(r,t)，并初始化φ(r,t)；

（3）根據(jù)式（10）和（12），分別計算權(quán)重系數(shù)和離散后的散射相函數(shù)；

根據(jù)新經(jīng)濟地理學(xué)和空間集聚等理論，在產(chǎn)業(yè)集聚過程中，由于生產(chǎn)要素和知識技術(shù)的集聚使得資源得到合理利用，基礎(chǔ)設(shè)施得到共享，進而使邊際成本降低；隨著大量產(chǎn)業(yè)集聚于同一區(qū)域，由于能源、勞動力等是有限的，過度集聚會造成邊際成本上升。因此，假設(shè)產(chǎn)業(yè)集聚和邊際成本可能具有非線性關(guān)系：

（4）使用式（7），計算平衡態(tài)分布函數(shù)S(r,t)；

（6）根據(jù)式（13）和（14），在邊界節(jié)點使用羅賓（Robin）邊界條件進行計算；

（7）通過式（15）計算相對差異，如最大差異值低于閾值（本文取10-2），則終止迭代過程，否則返回步驟4；

（8）計算格子節(jié)點r∈Ωk處的光子密度φ(r,t)：

2.2 基于LB前向模型的FDOT重建

熒光分子斷層成像方程由下式給出：

（1）初始化待重建光場分布S0；

（2）采集邊界測量值φmeas；

（3）通過LB前向模型計算生成靈敏度矩陣W；

（4）計算當前估計值WS'和真實測量值φmeas的差值φmeas-WS'（t代表計算時刻）；

（5）用當前差值矯正下一時刻的光場分布：

（6）重復(fù)步驟（4）和（5），使得相應(yīng)差值φmeas-WSt足夠小，從而實現(xiàn)待重建光場分布的求解。

3 實驗與結(jié)果分析

本文分別進行了數(shù)值仿體和物理仿體實驗，通過對比不同重建結(jié)果，評估了LB前向模型用于FDOT成像時的性能表現(xiàn)。

數(shù)值仿體實驗以直徑3 cm、高5 cm的圓柱作為成像物體，兩個高、寬都為0.2 cm的圓柱標記物被放置于成像物體中，使用不同位置的標記物，如表1所示，完成不同實驗?zāi)Ｐ偷脑O(shè)置，具體如圖2(a)所示。為簡化問題，圓柱體被填充均勻介質(zhì)，相應(yīng)的光學(xué)參數(shù)μa(r)、μs(r)、g分別為0.3 cm-1、10 cm-1和0.75。

在實驗計算中，LB模型被用于光傳播模擬，通過D3Q6進行速度離散，并將物理空間離散為15×15×25的節(jié)點網(wǎng)格，在離散節(jié)點進行LB迭代計算，并在最大差異值小于10-2時終止，具體計算流程見2.1。其中，實驗以DE作為對比，進行相似程度的離散，對應(yīng)5 831個節(jié)點，并在離散節(jié)點使用限元方法進行求解，相應(yīng)計算在Toast++上完成[19]。

表1 數(shù)值仿體實驗及標記物位置Tab.1 Numerical simulation experiment and targets positions

圖2 數(shù)值仿體模型和不同前向模型的重建結(jié)果Fig.2 Numerical phantom model and reconstruction results from different forward models

對應(yīng)不同的實驗組，圖2給出了基于LB和基于DE前向模型的FDOT重建結(jié)果，左右兩邊分別對應(yīng)重建結(jié)果的三維、二維顯示，二維結(jié)果對應(yīng)三維空間光源所在層（對應(yīng)a-c的黑色圓圈）上的重建光場分布。通過對比可以看出，基于LB前向模型的重建結(jié)果相似于基于DE模型的重建結(jié)果，和正確標記物位置都有著良好重疊；相對而言，基于LB前向模型的FDOT重建結(jié)果在正確區(qū)域上顯示得更加飽滿，而基于DE的重建結(jié)果有著更少的偽影表現(xiàn)。

更多的細節(jié)對比由圖3給出，對應(yīng)圖2(e)和2(f)白色虛線上的重建結(jié)果。通過對比可以看出，基于LB前向模型的重建結(jié)果能和標記物的正確位置有著更好的重疊，非重疊區(qū)域的面積更小。另外，在中間非標記物的區(qū)域，基于LB前向模型的重建結(jié)果值為0，說明相應(yīng)的重建能將兩個標記物的位置進行完全分離。其次，實驗還對比了不同前向模型重建結(jié)果的對比信噪比（Contrast-to-noise ratio，CNR），將真實標記存在和不存在空間分別定義為感興趣區(qū)域（ROI）和背景區(qū)域（BCK），相應(yīng)的計算公式如下：

圖3 不同前向模型對應(yīng)重建結(jié)果的一維對比Fig.3 The 1D comparison of reconstruction results from different forward models

其中，μROI和μBCK分別為ROI和BCK區(qū)域的平均重建強度；σ2ROI和σ2BCK分別為相應(yīng)區(qū)域重建強度的方差；wROI和wBCK分別代表相應(yīng)區(qū)域的容積大小。另外，實驗還記錄了不同前向模型的FDOT計算時間。上述實驗的CNR值和計算時間由表2給出。通過觀察可以發(fā)現(xiàn)，LB和DE模型所對應(yīng)重建結(jié)果的CNR值非常接近；基于LB的FDOT重建比基于DE的FDOT重建快上將近5倍。

表2 不同前向模型重建結(jié)果的CNR值和時間Tab.2 The CNR and calculation time from different models

本文除數(shù)值仿體實驗外，還進行了物理仿體實驗。使用FDOT/XCT混合的成像系統(tǒng)，完成成像物體表面光學(xué)信號的采集。實驗以直徑為3 cm、高為4.4cm的圓柱為成像物體，圓柱體內(nèi)填充均勻介質(zhì)，相應(yīng)的光學(xué)參數(shù)μ a(r)、μ s(r)、g分別為0.3 cm-1、10 cm-1和0.75，一個填充了熒光標記物（ICG）、直徑為0.4 cm的管子被放置于仿體模型中，如圖4（a）所示。

在實驗計算過程中，使用了類似于數(shù)值仿體實驗的LB、DE計算。相應(yīng)的重建結(jié)果由圖4給出，左右兩邊分別對應(yīng)重建結(jié)果的三維、二維顯示，其中的二維顯示對應(yīng)三維圓柱黑色圓圈所在層的內(nèi)容。通過三維顯示對比可以看出，基于LB的重建結(jié)果呈現(xiàn)出了類似于真實標記物形狀的圓柱形；而基于DE的重建結(jié)果呈現(xiàn)出了和真實標記物形狀差異較大的球型；稍顯不足的是，基于LB的FDOT重建結(jié)果有著較高的偽影。通過重建結(jié)果的二維對比，可以看出基于LB的重建結(jié)果在正確區(qū)域顯得更亮，能和真實的重建區(qū)域進行更好擬合。

圖4 物理仿體模型和不同前向模型的重建結(jié)果Fig.4 Physical phantom model and reconstruction results from different forward models

更多的細節(jié)對比由圖5給出，分別對應(yīng)圖4（e）和2（f）白色虛線上的重建結(jié)果。通過對比可以看出，基于LB和DE的FDOT重建結(jié)果都能和標記物所在位置的輪廓進行很好重合。

圖5 不同前向模型的一維重建結(jié)果對比Fig.5 The 1D comparison of reconstruction results from different forward models

其次，實驗還記錄了上述物理仿體實驗中不同前向模型FDOT重建結(jié)果的CNR值和相應(yīng)的計算時間，由表3給出。經(jīng)對比表發(fā)現(xiàn)，LB前向模型的使用可以使得FDOT的重建速度提升5倍左右，CNR提升2倍左右。

表3 不同前向模型重建結(jié)果的CNR值和時間Tab.3 The CNR and calculation time from different models

4 結(jié)論

FDOT作為一種新型成像技術(shù)，具有很高的應(yīng)用價值。其中的前向模型對于FDOT的成像性能有著直接影響。為了進一步提升FDOT的性能，有必要建立一個計算快速且準確的前向模型。本文，我們提出了一種基于LB方法的前向模型，并將其用于FDOT成像。本文通過數(shù)字、物理仿體實驗，驗證了LB前向模型對于FDOT的性能提升。實驗證明，LB前向模型的使用可以讓FDOT的計算時間得到大量縮短，且具有較好的精度表現(xiàn)。