本發(fā)明涉及一種光在介質(zhì)中傳播的描述方法及其應用，具體地，涉及一種基于格子波爾茲曼模型的光在介質(zhì)中傳播的描述方法及其在熒光擴散斷層成像中的應用。

背景技術：

隨著光學成像技術的發(fā)展，熒光擴散斷層成像(FDOT)技術已經(jīng)受到越來越多的關注。FDOT利用近紅外光激發(fā)組織體內(nèi)的熒光團，并根據(jù)在組織體邊界檢測到的光強值重建出熒光團濃度。該技術能夠準確定位和量化在體熒光團，可用于腫瘤的早期診斷和分子、細胞活性的監(jiān)測，具有無創(chuàng)/微創(chuàng)、相對深度深、高靈敏性和特異性等優(yōu)點。

FDOT技術是對擴散斷層成像(DOT)技術的直接擴展，其性能主要依賴于兩部分：準確描述光在組織中傳播的前向模型和重建算法。重建對于前向模型引入的噪聲敏感，且整個成像過程所需時間大部分由前向模型消耗。目前主要的前向模型求解方法有蒙特卡羅(Monte Carlo，MC)法、解析法和數(shù)值法。

MC法是20世紀40年代發(fā)展起來的一門計算科學，該方法是一種基于統(tǒng)計理論的數(shù)值方法。1983年Wilson等人將其用于組織光學領域。MC法將光在組織體中的運動看作一個隨機過程，在這個過程中光子或光子包在組織中隨機行走從而進行傳播，并利用一系列的隨機數(shù)來模擬光子在介質(zhì)中的運動和行徑，追蹤每個光子或光子包的歷史，然后對獲得的信息加以分析。

MC法是對實際物理過程的真實描述，能對任何復雜幾何形狀域進行計算，并能獲得精確的結果。但蒙特卡羅法耗時非常嚴重。

解析法又稱作格林函數(shù)法，各個位置的解可看作是格林函數(shù)和光源分布的卷積，從而易于計算。格林函數(shù)雖然能夠基本求解任意形狀下的擴散方程(Diffuse Equation,DE)，但耗時嚴重。因此，Jorge等人進一步提出了基爾霍夫近似(Kirchhoff approximation，KA)。KA的本質(zhì)是對格林函數(shù)的規(guī)范化，它是一種能模擬任意幾何體光強分布的線性方法，在求解正問題時，KA不考慮矩陣求逆問題，所以能用來生成系統(tǒng)的敏感函數(shù)或者加權值，因此，諸如代數(shù)重建技術(algebraic reconstruction techniques,ART)之類的逆問題都可用KA求解。由于KA求解過程中不含矩陣求逆問題，所以計算時間在所述的這些方法中是最少的，但是準確性問題是一個弊端。

數(shù)值法主要包括有限元法、有限差分法、有限體積法和邊界元法等。這些方法一般都是對近似的DE求解，其中有限元法(Finite Element Method,FEM)最為常用。該方法利用變分原理，建立與DE初值問題等價的積分方程，然后對整個高散射介質(zhì)進行整體離散以獲得等價的有限元方程組，最后求解該方程組獲得相應解。還有一些基于有限元的仿真軟件也常用于多種光學應用中，如COMSOL等。

光在生物組織中的傳播可以近似看作是各向異性擴散，而格子波爾茲曼(Lattice Boltzmann，LB)模型是一種常用的描述各向異性擴散的模型。加之其具有清晰的物理思想、固有的穩(wěn)定性和并行性，使之非常適合于FDOT技術中前向模型的建立。通過設計LB模型模擬光在組織體中的運動，進而利用重建算法重建出光學圖像。因此LB模型方法為實現(xiàn)FDOT的快速高效和準確性提供了實現(xiàn)途徑。

技術實現(xiàn)要素：

為了解決以上諸多問題，本發(fā)明提出了一種基于格子波爾茲曼模型的光在介質(zhì)中傳播的描述方法，其由于LB模型具有并行性結構，易于實現(xiàn)并行運算，相比于MC法，大大提高了熒光團重建效率；相對于FEM，在具有相當?shù)闹亟ň然A上，大幅提高了計算效率。

本發(fā)明的技術解決方案如下：一種基于格子波爾茲曼模型的光在介質(zhì)中傳播的描述方法，其包括以下步驟：

步驟一，初始化組織體的尺寸、獲取光源和檢測器的位置，設置組織體的組織參數(shù)，包括吸收系數(shù)μ_am、散射系數(shù)μ_sm和各向異性因子g；

步驟二，熒光團受激發(fā)出熒光；

步驟三，設計擴散系數(shù)，與距離光源點的距離有關，當距離光源點遠時，擴散系數(shù)??；反之，當距離光源點近時，擴散系數(shù)大，這樣能實現(xiàn)光的強度隨距離光源點的距離增大而衰減；

步驟四，建立格子波爾茲曼(Lattice Boltzmann，LB)模型的D3Q7模型(表示3維空間7個離散速度方向)，初始化平衡態(tài)分布函數(shù)；

步驟五，計算LB模型(指“格子波爾茲曼模型”)中的碰撞過程；

步驟六，判斷光顆粒是否運動到邊界，若是則對邊界處的光顆粒進行邊界處理，否則直接進入步驟七；

步驟七，計算LB模型中的遷移過程；

步驟八，更新各節(jié)點的光強值，即對各節(jié)點處強度值求和；

步驟九，判斷是否到達最大的迭代次數(shù)，如果是則顯示輸出，否則返回步驟三；

步驟十，更新輸出結果。

優(yōu)選地，所述擴散系數(shù)為利用朗伯-比爾定律設計的擴散系數(shù)，表示為下式：

其中，L為熒光射入組織體的厚度，該擴散系數(shù)能添加與初始熒光光強有關的系數(shù)；表示為未加入組織體厚度因素的擴散系數(shù)。

對于散射作用遠大于吸收作用的生物組織，用約化散射系數(shù)(reduced scattering coefficient)μ_sm′＝μ_sm(1-g)代替散射系數(shù)μ_sm，則上式能改寫為下式：

優(yōu)選地，所述步驟四中，提出對LB模型演化方程的吸收項的設計，輻射傳輸方程(RTE)、擴散方程DE、LB方法的宏觀方程分別寫成以下三式：

其中，式中為組織體吸收和散射造成的輻射率的減??；μ_a(r)為組織體位于空間節(jié)點r上的吸收系數(shù)；μ_s(r)為組織體位于空間節(jié)點r上的散射系數(shù)；為位于空間節(jié)點r處方向上的源項；式中的μ_aΦ(r，t)+Q(r，t)表示組織體吸收造成的輻射率的減小和源項；因此，對于LB方法，設計吸收項F表示吸收造成的節(jié)點處輻射值的減少和源項，表示為：

F＝-μ_af_α(r，t)+Q(r，t)；

其中，μ_a表示為吸收系數(shù)，f_α(r，t)表示為在t時刻、節(jié)點r處的光顆粒數(shù)，Q(r，t)表示為在t時刻、節(jié)點r處的源項光子密度。

優(yōu)選地，所述吸收項保證了能量守恒，由于光在組織體擴散時不斷被組織體吸收，該部分光子不參與后續(xù)的擴散，本發(fā)明設計吸收項以確保能量守恒。

優(yōu)選地，所述步驟六中，提出采用羅賓(Robin)邊界條件，能準確求解LB模型，該邊界條件考慮了組織體的折射率n_i大于周圍介質(zhì)的折射率n_t時邊界出射光的反射和折射效應，反射概率能夠通過反射系數(shù)R表示，反射效應和折射效應的邊界條件上的輻射率表示為：

其中，r表示位于組織邊界上的空間點，時，表示由方向反射到方向上的輻射率；時，表示由方向折射到方向上的輻射率，組織體與周圍介質(zhì)的折射率相差不大時，折射角相對于入射的偏移量不大，并且出射的輻射率不會再對組織體內(nèi)產(chǎn)生影響。

對于由反射回到組織體內(nèi)的光顆粒仍然遵守RTE描述，在組織體內(nèi)繼續(xù)傳播，而折射出組織體的光顆粒能被探測器測量得到，表示為：

在LB方法中，光顆粒始終沿著格線在節(jié)點間運動，即與e_α相隔θ_i，即此時，表示為：

光顆粒并非位于節(jié)點時，利用LB方法的插值格式處理。

優(yōu)選地，所述LB模型采用D3Q7模型模擬光子在組織體中的擴散，所述模型包括所有DmQn類模型(表示m維空間n個離散速度方向)，常用的LB模型包括D2Q5、D3Q7、D3Q15、D3Q19等，LB模型的步驟主要包括：

步驟十一，初始化組織體尺寸，設置組織體的參數(shù)，包括吸收系數(shù)、散射系數(shù)和各向異性因子，并獲取光源和探測器位置；

步驟十二，使用近紅外光激發(fā)熒光團發(fā)出熒光；

步驟十三，分別使用提出的LB模型模擬激發(fā)光和熒光的擴散過程；

步驟十四，根據(jù)步驟十三得到的結果，采用重建算法重建組織體內(nèi)的熒光團濃度。

優(yōu)選地，所述步驟十一中的組織體能夠模擬所有組織的光學參數(shù)，所述的熒光團能夠設置為任意數(shù)量、尺寸大小及位置。

優(yōu)選地，所述組織體包括第一組織體、第二組織體、第三組織體，組織體是填充有脂肪乳的參數(shù)確定的仿體。

優(yōu)選地，所述擴散系數(shù)控制LB模型的擴散速度，本發(fā)明對常用擴散系數(shù)進行優(yōu)化，利用朗伯-比爾定律設計出與組織體厚度和熒光光強相關的擴散系數(shù)。

優(yōu)選地，所述組織體是基于FDOT(熒光擴散斷層成像)這一應用使用的具體介質(zhì)。該發(fā)明廣泛適用于多種介質(zhì)。所述的模型為一種新型的描述光在介質(zhì)中傳播的模型，并基于該模型構建了新型FDOT系統(tǒng)。所述的方法可廣泛應用于諸如光學成像、光聲成像、電磁源成像、微波成像、電阻抗成像等圖像重建中前向模型的建立和求解。

本發(fā)明的積極進步效果在于：本發(fā)明由于LB模型具有并行性結構，易于實現(xiàn)并行運算，相比于MC法，大幅提高了熒光團重建效率；相對于FEM，在具有相當?shù)闹亟ň然A上，大幅提高了計算效率。本發(fā)明所述的模型為一種新型的描述光在介質(zhì)中傳播的模型，并基于該模型構建了新型FDOT系統(tǒng)。所述的模型可廣泛應用于圖像重建中前向模型的建立和求解。

附圖說明

圖1是本發(fā)明基于格子波爾茲曼模型的光在介質(zhì)中傳播的描述方法中基于D3Q7的LB模型的各矢量方向圖。

圖2是本發(fā)明基于格子波爾茲曼模型的光在介質(zhì)中傳播的描述方法的流程示意圖。

圖3是本發(fā)明中基于LB方法的FDOT重建算法流程圖。

圖4是本發(fā)明中Robin邊界條件下光顆粒在格子邊界的反射與折射圖。

圖5是本發(fā)明仿真模型示意圖。

圖6a至6h是第一組織體基于本發(fā)明和COMSOL方法得到的三維擴散結果圖。

圖7a至7h是第二組織體基于本發(fā)明和COMSOL方法得到的三維擴散結果圖。

圖8a至8h是第三組織體基于本發(fā)明和COMSOL方法得到的三維擴散結果圖。

圖9a至9d是本發(fā)明中基于LB方法得到的三維重建結果圖。

具體實施方式

下面結合附圖給出本發(fā)明實施例，以具體說明本發(fā)明的技術方案。

如圖2、圖4所示，本發(fā)明公開了一種基于格子波爾茲曼模型的光在介質(zhì)中傳播的描述方法，其主要包括以下步驟：

步驟一，初始化組織體的尺寸、獲取光源和檢測器的位置，設置組織體的組織參數(shù)，包括吸收系數(shù)μ_am、散射系數(shù)μ_sm和各向異性因子g；

步驟二，熒光團受激發(fā)出熒光；

步驟三，設計擴散系數(shù)，與距離光源點的距離有關，當距離光源點遠時，擴散系數(shù)??；反之，當距離光源點近時，擴散系數(shù)大，這樣能實現(xiàn)光的強度隨距離光源點的距離增大而衰減；

步驟四，建立LB模型的D3Q7模型(表示3維空間7個離散速度方向)，初始化平衡態(tài)分布函數(shù)；

步驟五，計算LB模型中的碰撞過程；

步驟六，判斷光顆粒是否運動到邊界，若是則對邊界處的光顆粒進行邊界處理，否則直接進入步驟七；

步驟七，計算LB模型中的遷移過程；

步驟八，更新各節(jié)點的光強值，即對各節(jié)點處強度值求和；

步驟九，判斷是否到達最大的迭代次數(shù)100，如果是則顯示輸出，否則返回步驟三；

步驟十，更新輸出結果。

所述步驟三中的擴散系數(shù)為利用朗伯-比爾定律設計的擴散系數(shù)，表示為式(1)：

μ_sm表示為熒光波長λ_m下的散射系數(shù)，μ_am表示為熒光波長λ_m下的吸收系數(shù)，L表示為熒光射入組織體的厚度，g(κ，L)表示為利用朗伯-比爾定律設計的擴散系數(shù)。該擴散系數(shù)能添加與初始熒光光強有關的系數(shù)。

對于散射作用遠大于吸收作用的生物組織，用約化散射系數(shù)(reduced scattering coefficient)μ_sm′＝μ_sm(1-g)代替散射系數(shù)μ_sm，則式(1)能改寫為式(2)：

μ_sm表示為熒光波長λ_m下的散射系數(shù)，μ_am表示為熒光波長λ_m下的吸收系數(shù)，L表示為熒光射入組織體的厚度，g表示為各向異性因子，g(κ，L)表示為利用朗伯-比爾定律設計的擴散系數(shù)。

所述步驟四中，提出對LB模型演化方程的吸收項的設計，RTE、DE、LB方法的宏觀方程分別寫成式(3)、式(4)、式(5)：

式(3)中，c表示為光在組織體中的傳播速度，表示為組織體中r點，某個確定方向在某時間點t的輻射率，μ_α(r)表示為該組織體位于空間節(jié)點r上的吸收系數(shù)；μ_s(r)表示為該組織體位于空間節(jié)點r上的散射系數(shù)，表示為散射相位函數(shù)，表示為位于空間節(jié)點r處，某個確定方向上的源項輻射率。

式(3)中表示為組織體吸收和散射造成的輻射率的減??；表示為源項帶來的輻射率增加。

式(4)中，c表示為光在組織體中的傳播速度，Φ(r，t)表示為節(jié)點r處，在某時間點t的光子密度，κ表示為擴散系數(shù)，μ_a表示為吸收系數(shù)，Q(r，t)表示為節(jié)點r處，在某時間點t的光源項的光子密度。

式(5)中，τ表示為松弛時間，e_α表示為光顆粒速度矢量，表示為局部平衡態(tài)分布函數(shù)。

因此，對于LB方法，設計吸收項F表示吸收造成的節(jié)點處輻射值的減少和源項，表示為式(6)：

F＝-μ_af_α(r，t)+Q(r，t) (6)；

μ_a表示為吸收系數(shù)，f_α(r，t)表示為在t時刻、節(jié)點r處的光顆粒數(shù)，Q(r，t)表示為在t時刻、節(jié)點r處的源項光子密度。

所述吸收項保證了能量守恒，由于光在組織體擴散時不斷被組織體吸收，該部分光子不參與后續(xù)的擴散，本發(fā)明設計吸收項以確保能量守恒。

所述步驟六中，提出采用Robin邊界條件，該邊界條件考慮了組織體的折射率n_i大于周圍介質(zhì)的折射率n_t時邊界出射光的反射和折射效應，反射概率能夠通過反射系數(shù)R表示，反射效應和折射效應的邊界條件上的輻射率表示為式(7)：

r表示為邊界上的空間點，表示為出射光方向，表示為該邊界點上的外向單位法向向量。時，表示由方向反射到方向上的輻射率；時，表示由方向折射到方向上的輻射率，組織體與周圍介質(zhì)的折射率相差不大時，折射角相對于入射的偏移量不大，并且出射的輻射率不會再對組織體內(nèi)產(chǎn)生影響。

對于由反射回到組織體內(nèi)的光顆粒仍然遵守RTE描述，在組織體內(nèi)繼續(xù)傳播，而折射出組織體的光顆粒能被探測器測量得到，表示為式(8)：

r表示為邊界上的空間點，表示為該邊界點上的外向單位法向向量，表示為在某空間點r上，某時刻t時在方向上的輻射率。

在LB方法中，光顆粒始終沿著格線在節(jié)點間運動，即與e_α相隔θ_i，即此時，表示為式(9)：

e_α表示為光顆粒速度矢量。

光顆粒并非位于節(jié)點時，利用LB方法的插值格式處理。

所述模型通過利用LB模型模擬光子在組織體中的擴散，所述模型包括所有DmQn類模型(表示m維空間n個離散速度方向)，常用的LB模型包括D2Q5、D3Q7、D3Q15、D3Q19等，LB模型的步驟主要包括：

步驟十一，初始化組織體尺寸，設置組織體的參數(shù)，包括吸收系數(shù)、散射系數(shù)和各向異性因子，并獲取光源和探測器位置；

步驟十二，使用近紅外光激發(fā)熒光團發(fā)出熒光；

步驟十三，分別使用提出的LB模型模擬激發(fā)光和熒光的擴散過程；

步驟十四，根據(jù)步驟十三得到的結果，采用重建算法重建組織體內(nèi)的熒光團濃度。

所述步驟十一中的組織體能夠模擬所有組織的光學參數(shù)，所述的熒光團能夠設置為任意數(shù)量、尺寸大小及位置。

所述組織體包括第一組織體、第二組織體、第三組織體，是填充有脂肪乳的參數(shù)確定的仿體，本發(fā)明中用到3種參數(shù)不同的第一組織體、第二組織體、第三組織體。如下表1所示；所述光源為產(chǎn)生發(fā)射波長為λ_x激發(fā)光的激光器，所述探測器為CCD相機；所述的LB模型：通過利用LB模型模擬光子在組織體中的擴散，常用的LB模型包括D2Q5—二維5個方向擴散，D3Q7—三維7個方向擴散。所述的基于LB模型的FDOT系統(tǒng)在包括PC或服務器在內(nèi)的硬件計算平臺中運行。

表1：圓和圓柱體仿真實驗的光學參數(shù)設置

所述擴散系數(shù)(即松弛因子，與松弛時間τ成反比)控制LB模型的擴散速度，本發(fā)明對常用擴散系數(shù)進行優(yōu)化，利用朗伯-比爾定律設計出與組織體厚度和熒光光強相關的擴散系數(shù)。

如圖5所示，本發(fā)明在第一組織體、第二組織體、第三組織體中設計4個case(實驗例)，利用LB模型的D3Q7模型模擬光子擴散過程，

如圖1所示，LB模型的D3Q7模型是由均勻網(wǎng)格組成，其中網(wǎng)格節(jié)點對應于組織體離散后的節(jié)點，即光顆粒沿著網(wǎng)格線在相鄰的網(wǎng)格點間運動。4個實驗例的仿真模型如圖5，所述D3Q7模型如圖1，包含有7個方向的離散速度，分別表示為式(10)：

α表示為光顆粒運動方向下標，e_αD3Q7表示為光顆粒速度矢量。

D3Q7模型擁有7個離散速度，模型中微觀粒子擴散演化方程為式(11)：

f_α(r，t)表示為在t時刻、節(jié)點r處的光顆粒數(shù)，表示為依賴于局部強度和速度的局部平衡分布函數(shù)，ρ為光顆粒密度，e_α表示為光顆粒速度矢量，τ表示為松弛時間，δt表示為時間步長，F(xiàn)_α表示為吸收項。

其中，τ為松弛時間，與前面所述擴散系數(shù)的關系為式(12)：

g(κ，L)表示為利用朗伯-比爾定律設計的擴散系數(shù)，τ表示為松弛時間。

在D3Q7模型中，A_α表示為式(13)：

α表示為光顆粒運動方向下標，A_α表示為平衡態(tài)加權系數(shù)。

FDOT的完整過程包含兩次DOT的過程，第一次模擬光在組織體中的傳播，第二次模擬熒光團發(fā)出的熒光在組織體中的傳播，兩次過程極為相似。因此本發(fā)明假設熒光團已經(jīng)被激發(fā)，第一組織體中的三維仿真方法直接采用如圖2所示的流程，待擴散完成后，輸出結果。如圖3所示，F(xiàn)DOT重建算法包括以下步驟:

步驟二十，初始化；

步驟二十一，實驗得測量值；

步驟二十二，求解正問題，得值；

步驟二十三，ART求解方程；

步驟二十四，判斷條件是否成立，是則轉步驟二十五，否則轉步驟二十一；

步驟二十五，結束。

第一組織體擴散結果如圖6a至6h所示，其重建結果如圖9a至9d所示，上述實施的硬件條件為處理器Xeon W3550(英特爾至強W3550處理器)，主頻3.07GHz，內(nèi)存8G；軟件配置為Window7系統(tǒng)和matlab2012b。對于3中組織體的四個實驗例，仿真所需時間與COMSOL(多物理場仿真軟件)結果對比如下表2所示。

表2：LB方法和COMSOL進行三維仿真實驗所需時間

第二組織體實施步驟與上述一致，利用上述4種實驗例進行基于LB模型的光子擴散，擴散結果如圖7a至7h、9a至9d所示，仿真所需時間與COMSOL結果對比如表2所示。

第三組織體實施步驟與上述一致，利用上述4種實驗例進行基于LB模型的光子擴散，擴散結果如圖8a至8h、9a至9d所示，仿真所需時間與COMSOL結果對比如表2所示。

實驗顯示，針對FDOT，采用基于LB模型前向模型能夠得到較好的重建效果，并且大幅提高了計算效率，該方法有望用到實際應用中。

以上所述的具體實施例，對本發(fā)明的解決的技術問題、技術方案和有益效果進行了進一步詳細說明，所應理解的是，以上所述僅為本發(fā)明的具體實施例而已，并不用于限制本發(fā)明，凡在本發(fā)明的精神和原則之內(nèi)，所做的任何修改、等同替換、改進等，均應包含在本發(fā)明的保護范圍之內(nèi)。