一種基于格子玻爾茲曼前向模型的熒光擴(kuò)散斷層成像系統(tǒng)

岑星星，嚴(yán)壯志,，吳煥迪

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

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

0 引言

熒光擴(kuò)散斷層成像（Fluorescence Diffuse Optical Tomography,FDOT）等光學(xué)分子影像技術(shù)，被認(rèn)為是現(xiàn)階段生物醫(yī)學(xué)取得進(jìn)一步突破的最大可能[1]。FDOT通過外部光源來激發(fā)體內(nèi)的熒光分子探針使其發(fā)出熒光，并利用體表采集到的光學(xué)信號完成體內(nèi)分子探針分布的重建工作。相較于其它傳統(tǒng)成像方式，F(xiàn)DOT具有高敏感性、高通量、低成本、安全無創(chuàng)等優(yōu)勢[2]。近年來，隨著成像系統(tǒng)、成像算法、光學(xué)分子探針等領(lǐng)域的快速發(fā)展，F(xiàn)DOT逐步走向成熟，并有望在預(yù)臨床和臨床領(lǐng)域取得廣泛應(yīng)用。

FDOT相關(guān)系統(tǒng)的開發(fā)一直是當(dāng)前生物醫(yī)學(xué)工程領(lǐng)域的研究熱點(diǎn)。獨(dú)立開源系統(tǒng)，例如SHENGHAN等[3]基于蒙特卡羅（Monte Carlo,MC）前向模型提出的用于組織體光傳播模擬的MOSE平臺；SCHWEIGER等[4]基于擴(kuò)散方程（Diffusion Equation,DE）前向模型提出的用于組織體光傳播模擬及相應(yīng)光學(xué)重建的Toast++平臺，都取得了廣泛應(yīng)用。另一方面是基于商用系統(tǒng)，如通過COMSOL、ANSYS等有限元分析工具進(jìn)行輻射傳輸方程（Radiative Transfer Equation,RTE）、DE等前向模型求解，實(shí)現(xiàn)光成像系統(tǒng)的搭建，具有更好的操作體驗(yàn)。光學(xué)成像的相關(guān)系統(tǒng)平臺已有很多，但大多數(shù)基于傳統(tǒng)前向模型進(jìn)行搭建，包括DE、RTE等前向模型。

傳統(tǒng)的前向模型有著各自的優(yōu)勢但也存在著明顯的缺陷。RTE作為麥克斯韋方程的近似，能對生物組織體中的光傳播過程進(jìn)行準(zhǔn)確描述[5]。然而，RTE同時也是結(jié)構(gòu)復(fù)雜的微積分方程，需要很高的計(jì)算成本，并且在復(fù)雜的生物組織應(yīng)用環(huán)境中，往往無法直接通過解析法進(jìn)行求解[6]。另一方面，DE作為RTE的一階近似，具有結(jié)構(gòu)簡單、易于求解、計(jì)算效率高等優(yōu)勢，是當(dāng)前FDOT應(yīng)用中使用最多的前向模型[7-8]。然而，作為最低階數(shù)的近似方程，高度的簡化處理，極大地?fù)p耗了該模型的精度，并帶給了DE局限性：①不能用于非散射介質(zhì)中的光傳播描述；②不能對光源和邊界處的光傳播過程進(jìn)行準(zhǔn)確描述[9]。RTE、DE等傳統(tǒng)模型存在的明顯缺陷，限制了相應(yīng)成像系統(tǒng)的更好發(fā)展。

格子玻爾茲曼（Lattice Boltzmann,LB）方法基于分子運(yùn)動理論，利用規(guī)則網(wǎng)格中的離散粒子運(yùn)動來對復(fù)雜的物理現(xiàn)象進(jìn)行模擬，實(shí)現(xiàn)相應(yīng)物理量的求解[10]。相比于傳統(tǒng)方法，LB方法的對應(yīng)數(shù)學(xué)模型具有計(jì)算結(jié)構(gòu)簡單、程序易于實(shí)現(xiàn)和并行化、計(jì)算穩(wěn)定性好等優(yōu)勢。因此，構(gòu)建基于LB方法的前向模型，并將其應(yīng)用于FDOT系統(tǒng)的搭建，是進(jìn)一步推動FDOT發(fā)展的極佳選擇。基于上述目的，筆者提出了基于LB前向模型的FDOT系統(tǒng)。

1 方法

1.1 LB前向模型

FDOT中的前向模型描述了光在生物組織中的傳播過程。基于LB方法的前向模型如下所示：

其中，ω(r)=μs(r)+β(r)代表反照率，wi為權(quán)重系數(shù)，可通過局部輻射場的守恒進(jìn)行求解，如下所示：

其中，ψk(s)為s的多項(xiàng)式。不同離散速度模型所對應(yīng)的權(quán)重系數(shù)如表1所示。

表1 LB方法中不同離散速度模型對應(yīng)的權(quán)重系數(shù)Tab.1 The weights of different discrete models in LB method

LB前向模型的邊界節(jié)點(diǎn)計(jì)算由羅賓（Robin）邊界條件給出。求出邊界節(jié)點(diǎn)r∈eΩk（k為離散空間節(jié)點(diǎn)的索引）上的光子密度Φ(r,t)：

并和羅賓邊界條件進(jìn)行關(guān)聯(lián)，得到如下邊界計(jì)算公式[18]：

其中，D(r)=1/[3(μa(r)+(1-g)μs(r))]為光擴(kuò)散系數(shù)，g為各向異性因子；是邊界處的外法線向量；當(dāng)組織周圍填充空氣，組織折射率為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。

設(shè)置不同位置的光源并進(jìn)行LB光傳播模擬求解，構(gòu)建出靈敏度矩陣W，并得到FDOT重建計(jì)算的相關(guān)線性方程，具體方程如下：

其中，Φmeas為邊界測量值；S為待重建的光場分布。

1.2 ART重建算法

ART可用于線性方程的求解，對式（6）的線性方程WS=Φmeas進(jìn)行ART求解，具體公式如下：

其中，n是當(dāng)前光源的索引值；t是迭代索引值；wn代表靈敏度矩陣W的第n行元素；φm是邊界測量向量Φmeas的第m個分量；λ'n代表松弛參數(shù)。其基本思想是通過當(dāng)前邊界估計(jì)值wn[st]T和真實(shí)邊界測量值Φm的差值Φm-wn[st]T，去校正下一時刻的未知光源分布估計(jì)值st+1。通過多次迭代計(jì)算，使得估計(jì)值st+1能滿足上述線性方程，從而實(shí)現(xiàn)相應(yīng)方程的求解。

2 系統(tǒng)設(shè)計(jì)

FDOT系統(tǒng)主要用于實(shí)現(xiàn)光傳播模擬和FDOT重建，系統(tǒng)的具體框架如圖1所示。對應(yīng)不同功能，將用戶界面劃分成兩個部分，包括光傳播模擬和FDOT重建，并單獨(dú)進(jìn)行操作。雖然兩個界面所要實(shí)現(xiàn)的功能不同，但相應(yīng)的程序操作流程有著極大的重合性，包括光學(xué)參數(shù)設(shè)置、LB前向模型計(jì)算等。因此，對兩個不同的界面使用一樣的核心處理模塊，包括：參數(shù)模塊、算法模型、結(jié)果顯示模塊和數(shù)據(jù)交互模塊，并進(jìn)行了詳細(xì)介紹。

圖1 FDOT系統(tǒng)框架圖Fig.1 FDOT system framework

2.1 參數(shù)模塊

參數(shù)模塊涉及系統(tǒng)中的各種輸入?yún)?shù)設(shè)置，具體的結(jié)構(gòu)如圖2所示。

圖2 FDOT系統(tǒng)參數(shù)結(jié)構(gòu)圖Fig.2 Parameter structure diagram in FDOT system

由圖2可見，參數(shù)模塊的主要工作內(nèi)容分為四塊：完成實(shí)驗(yàn)?zāi)Ｐ蛥?shù)、光源參數(shù)、組織體的光學(xué)參數(shù)和算法參數(shù)。實(shí)驗(yàn)?zāi)Ｐ蛥?shù)由ExperimentModel類進(jìn)行管理，包含模型形狀及相應(yīng)的長、寬、高。光源參數(shù)由Source類進(jìn)行管理，包含模擬光源的空間位置。組織體的光學(xué)參數(shù)由Optical類進(jìn)行管理，包含光吸收系數(shù)、散射系數(shù)和g。算法參數(shù)由Algorithm類文件進(jìn)行管理，包括LB前向模型計(jì)算時的迭代次數(shù)、DmQn模型、空間離散尺度、邊界條件和ART重建計(jì)算時的松弛參數(shù)等。

2.2 算法模塊

算法模塊與系統(tǒng)所用算法的計(jì)算實(shí)現(xiàn)相關(guān)，包括LB前向模型和ART重建算法的計(jì)算流程，算法相關(guān)結(jié)構(gòu)如圖3所示。

圖3 LB和ART算法結(jié)構(gòu)圖Fig.3 Algorithm structure diagram of LB and ART

由圖3可見，算法模塊的主要工作內(nèi)容分為兩個部分：LB前向模型的計(jì)算及相應(yīng)參數(shù)的讀取和結(jié)果存儲；ART的重建計(jì)算及相應(yīng)參數(shù)的讀取及結(jié)果存儲。LB前向模型的計(jì)算對應(yīng)方法RhoCalculation( )，ART重建算法的計(jì)算對應(yīng)方法ImDataCalculation( )，具體的計(jì)算流程如圖4所示。

圖4 兩種算法流程圖Fig.4 Two kinds of algorithm flow chart

經(jīng)LB計(jì)算后的結(jié)果存儲在變量Rho中，經(jīng)ART重建計(jì)算后的結(jié)果存儲于ImData中。

2.3 結(jié)果顯示模塊

結(jié)果顯示模塊，用于顯示計(jì)算后的結(jié)果，包括光傳播模擬時，經(jīng)LB前向模型計(jì)算后的光場分布；FDOT重建時，經(jīng)LB、ART計(jì)算后的重建結(jié)果。

觀察圖5可以知道，結(jié)果顯示模塊的主要功能分為兩個部分：對LB計(jì)算后的光場分布進(jìn)行三維、兩維和一維的顯示（二維、一維顯示內(nèi)容以光源位置為中心），分別對應(yīng)Show3D( )、Show2D( )和Show1D( )方法；對LB、ART重建計(jì)算后的重建結(jié)果進(jìn)行多個維度的顯示。

圖5 LB和ART重建結(jié)果顯示結(jié)構(gòu)圖Fig.5 Reconstruction result display diagram of LB and ART

2.4 數(shù)據(jù)交互模塊

數(shù)據(jù)交互模型主要用于完成系統(tǒng)計(jì)算結(jié)果的數(shù)據(jù)導(dǎo)出及FDOT重建時的數(shù)據(jù)導(dǎo)入工作。導(dǎo)出內(nèi)容包括：光傳播模擬時，經(jīng)LB前向模型計(jì)算后的光場分布數(shù)據(jù)；FDOT重建時，經(jīng)LB、ART計(jì)算后的重建結(jié)果數(shù)據(jù)。導(dǎo)入內(nèi)容為FDOT的待重建數(shù)據(jù)，包括FDOT重建時的激發(fā)光源位置、探測器位置及探測器上的測量值。

2.5 用戶圖形界面

用戶圖形界面分為光傳播模擬界面和FDOT重建界面。光傳播模擬界面分為參數(shù)設(shè)置區(qū)、操作區(qū)、模型顯示區(qū)和模擬結(jié)果顯示區(qū)，具體如圖6所示。在參數(shù)設(shè)置區(qū)設(shè)置相應(yīng)參數(shù)，得到如模型顯示區(qū)所顯示的實(shí)驗(yàn)?zāi)Ｐ?；在操作區(qū)進(jìn)行LB光傳播模擬計(jì)算后顯示在結(jié)果顯示區(qū)內(nèi)。

圖6 光傳播模擬界面Fig.6 Optical propagation simulation interface

FDOT重建的界面設(shè)計(jì)風(fēng)格和光傳播模擬界面保持一致，也分為參數(shù)設(shè)置區(qū)、操作區(qū)、模型顯示區(qū)和模擬結(jié)果顯示區(qū)，具體如圖7所示。操作流程類似于光傳播模擬界面的使用流程。

圖7 FDOT重建界面Fig.7 FDOT reconstruction interface

通過圖形界面的使用，可以使得FDOT系統(tǒng)的操作變得更加直觀、簡潔，并提升了該系統(tǒng)的可操作性和交互性。

3 系統(tǒng)驗(yàn)證和結(jié)果

3.1 光傳播模擬的可靠性驗(yàn)證

本研究通過對比MCxyz的MC光傳播模擬結(jié)果及Toast++的DE光傳播計(jì)算結(jié)果，來驗(yàn)證本系統(tǒng)是否能可靠用于光傳播模擬計(jì)算。實(shí)驗(yàn)?zāi)Ｐ筒捎弥睆綖? cm、高為5 cm的圓柱，在圓柱中心放置初始能量密度為1 nw/sr且持續(xù)發(fā)光的點(diǎn)光源，如圖8所示。以MC的計(jì)算結(jié)果作為標(biāo)準(zhǔn)。

圖8 FDO可靠性模擬實(shí)驗(yàn)?zāi)Ｐ虵ig.8 The simulation model of FDOT to test reliability

在實(shí)驗(yàn)過程中，實(shí)驗(yàn)?zāi)Ｐ捅灰来翁畛洳煌鈱W(xué)參數(shù)的介質(zhì)，模擬常規(guī)和低吸收組織環(huán)境，分別對應(yīng)Case1和Case2，如表2所示。

表2 實(shí)驗(yàn)介質(zhì)的光學(xué)參數(shù)Tab.2 Optical parameters of mediums in experiment

本系統(tǒng)中的LB前向模型經(jīng)D3Q26離散化，離散的網(wǎng)格尺寸為0.1 cm，具體的計(jì)算流程由章節(jié)2.2的圖4給出。

3.2 FDOT重建的可靠性驗(yàn)證

通過與待重建標(biāo)記物的真實(shí)位置進(jìn)行對比，來驗(yàn)證本系統(tǒng)是否能可靠用于FDOT重建。實(shí)驗(yàn)以直徑為3 cm、高為5 cm的圓柱模型作為成像物體，在成像物體內(nèi)部放置兩個高、直徑都為0.2 cm的圓柱形標(biāo)記物，標(biāo)記物的中心位置分別為(－0.4 cm,0.4 cm,2.7 cm)、(0.4cm,－0.4 cm,2.7 cm)，并在成像物體周圍放置探測器數(shù)據(jù)采集點(diǎn)。成像物體被填充均勻介質(zhì)，介質(zhì)的光學(xué)參數(shù)μa(r)、μs(r)、g分別為0.3 cm-1、10 cm-1和0.75。

3.3 驗(yàn)證結(jié)果分析

不同系統(tǒng)的光傳播模擬計(jì)算結(jié)果由圖9給出，顯示內(nèi)容為圓柱實(shí)驗(yàn)?zāi)Ｐ椭虚g虛線圓圈所在層上的光場分布對數(shù)值。通過對比可以發(fā)現(xiàn)，經(jīng)本系統(tǒng)LB模型計(jì)算的光場分布和經(jīng)MCxyz系統(tǒng)MC計(jì)算的光場分布更加接近，特別是在低吸收組織環(huán)境中，兩者的結(jié)果幾乎一致，如圖9（a2）～（c2）所示；經(jīng)Toast++系統(tǒng)DE計(jì)算的光場分布要明顯亮于其它系統(tǒng)的模型計(jì)算結(jié)果，圖9（b1）～（b2）所示。通過上述實(shí)驗(yàn)，驗(yàn)證了本系統(tǒng)用于光傳播模擬時的可靠性。

實(shí)驗(yàn)?zāi)Ｐ腿鐖D10（a）所示。經(jīng)本系統(tǒng)計(jì)算得到的FDOT重建結(jié)果，如圖10（b）（c）所示，分別對應(yīng)三維、二維顯示。為了進(jìn)行更直觀的對比，成像標(biāo)記物的真實(shí)位置被放置于重建結(jié)果中，如圖10（c）中的白色圓圈所示。通過對比可以看出，經(jīng)本系統(tǒng)計(jì)算得到的重建光源位置能和真實(shí)標(biāo)記物的正確位置進(jìn)行良好重合，無論是在三維顯示中，圖10（a）（b）所示，還是在二維顯示中，圖10（c）所示。通過上述實(shí)驗(yàn)，驗(yàn)證了本系統(tǒng)用于FDOT重建的可靠性。

4 總結(jié)

圖9 不同前向模型的二維計(jì)算結(jié)果對比Fig.9 The comparison of 2D results from different forward models

圖10 實(shí)驗(yàn)?zāi)Ｐ秃椭亟ńY(jié)果的三維顯示Fig.10 Experimental model and 3D display of reconstruction results

本研究以格子玻爾茲曼前向模型作為核心，搭建了新型的FDOT系統(tǒng)，經(jīng)光傳播模擬實(shí)驗(yàn)和FDOT重建實(shí)驗(yàn)，評估了所搭建系統(tǒng)的可靠性。主要內(nèi)容包括系統(tǒng)的核心算法、框架、核心模塊及用戶可視化界面。對該系統(tǒng)的實(shí)用性進(jìn)行分析，和其它商用系統(tǒng)相比，該系統(tǒng)的開源性和簡易性，能極大地降低FDOT成像的使用成本；和其它開源系統(tǒng)相比，由于新型前向模型的使用，使得該系統(tǒng)具有更佳的并行化及速度提升優(yōu)勢。因此，提出的FDOT系統(tǒng)具有很高的推廣價值。