基于有限元和k空間偽譜法的光聲成像仿真平臺

齊夢雨，嚴(yán)壯志, ， 趙麗麗

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

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

光聲成像技術(shù)是一種新型無創(chuàng)在體生物醫(yī)學(xué)成像技術(shù)。光聲成像結(jié)合了超聲和光學(xué)成像的優(yōu)點(diǎn)[1]，兼具超聲成像的高分辨率和純光學(xué)成像的高對比度特性，因而成為了生物醫(yī)學(xué)成像領(lǐng)域的研究前沿和熱點(diǎn)之一[2]。光聲成像的基本原理是光聲效應(yīng)[3]：使用短脈沖激光照射生物組織，組織體由于吸收能量導(dǎo)致的瞬時熱膨脹進(jìn)而輻射超聲波。通過檢測和采集這些攜帶生物組織特性的超聲信號，就可以利用一定的重建算法重建出生物組織內(nèi)部的光學(xué)吸收分布圖像。由于生物組織的光吸收特性與組織的生理特征、結(jié)構(gòu)形態(tài)特征、病變特性、代謝狀態(tài)甚至神經(jīng)活動密切相關(guān)，因此光聲成像技術(shù)已經(jīng)在腫瘤早期診斷[4-5]、生物分子成像[6]、生物組織功能成像[7]等生物醫(yī)學(xué)領(lǐng)域中取得了重大的進(jìn)展。

由于光聲成像技術(shù)是一種多模態(tài)的混合成像技術(shù)，該技術(shù)的研究常常需要昂貴的實(shí)驗(yàn)設(shè)備，以及隨之而來的繁雜的實(shí)驗(yàn)操作步驟。如果能將仿真技術(shù)應(yīng)用在光聲成像技術(shù)中，不僅可以避免真實(shí)實(shí)驗(yàn)所帶來的開銷，方便地調(diào)節(jié)光學(xué)、聲學(xué)以及各項(xiàng)系統(tǒng)參數(shù)，同時可以有選擇地忽略真實(shí)實(shí)驗(yàn)的系統(tǒng)噪聲和測量誤差，為理論研究提供有利條件。因此，光聲成像仿真技術(shù)對于光聲成像的理論研究及其在實(shí)際中的應(yīng)用有重要的指導(dǎo)意義，也是驗(yàn)證實(shí)驗(yàn)設(shè)計性能和算法性能的首選方案。

目前，計算機(jī)仿真技術(shù)已經(jīng)被廣泛應(yīng)用于光聲成像領(lǐng)域的研究中。例如，TREEBY等[8]基于Matlab平臺使用k空間偽譜法（k-space Pseudospectral Method）對光聲信號在有損介質(zhì)中的傳播進(jìn)行了模擬仿真；YUAN等[9]使用有限元方法在Matlab平臺上對光聲成像的后向過程進(jìn)行了模擬。這部分研究中的仿真主要是針對聲場仿真或重建過程仿真，并不涉及光源和光場部分的仿真，因此可擴(kuò)展性和適用性都受到了較大限制。對于包含光場和聲場的全過程光聲成像仿真，一般使用商用多物理場仿真軟件COMSOL Multiphysics或者有限元分析軟件ANSYS搭建仿真平臺進(jìn)行仿真。例如，MOON等[10]使用ANSYS對金屬元素顆粒進(jìn)行了全過程的光聲成像仿真；SOWMIYA等[11]使用COMSOL對光聲斷層成像的全過程進(jìn)行了仿真。而COMSOL或ANSYS軟件不僅價格高昂，專業(yè)性強(qiáng)，而且由于商業(yè)軟件的非開源性，內(nèi)部代碼和運(yùn)行機(jī)理無法得知。

針對上述問題，本文根據(jù)有限元法和k空間偽譜法[12]，借助開源工具箱k-Wave[13]有限元計算工具箱TOAST++[14]，設(shè)計并實(shí)現(xiàn)了光聲成像仿真平臺和圖形用戶界面(GUI) ，完成了對光聲成像全過程的模擬和仿真。光聲成像仿真平臺是基于Matlab搭建完成的。Matlab是由美國MathWorks公司出品的商業(yè)數(shù)學(xué)軟件，也是目前使用最普遍的數(shù)學(xué)軟件之一。Matlab內(nèi)部集成了極為完善的科學(xué)計算工具和圖形顯示工具，并且支持python、C/C++等其它擴(kuò)展語言。因此基于Matlab搭建的仿真平臺將更具擴(kuò)展性和普及性。

1 光聲成像理論基礎(chǔ)

在光聲成像技術(shù)中，一般使用周期性的脈沖激光照射目標(biāo)組織，吸收了激光能量的目標(biāo)組織會由于升溫而發(fā)生周期性的熱脹冷縮，從而輻射出超聲波。一般而言，為了有效產(chǎn)生超聲波，激勵源需要滿足熱學(xué)限制條件和壓力限制條件：① 脈沖激光持續(xù)時間小于組織的熱擴(kuò)散時間；② 脈沖激光的脈寬小于聲波在目標(biāo)組織中傳播所需的時間。在滿足這兩個條件的情況下，目標(biāo)組織由于激光照射而產(chǎn)生熱量的過程可以描述為：

其中，ρ是組織密度，Cp是比熱，T(r, t)表示位置r處、t時刻生物組織吸收光能量產(chǎn)生的溫度變化，λ是熱導(dǎo)率，H(r, t)是激光的能量沉積。由于熱學(xué)限制條件，這個過程可以當(dāng)作絕熱過程，所以簡化為：

如果假設(shè)目標(biāo)組織內(nèi)部的密度和聲速為常數(shù)，那么可以通過以下公式來描述聲場：

其中，矢量函數(shù)u(r, t)表示聲位移。同時，由于溫度上升導(dǎo)致的熱膨脹可以用熱膨脹方程來描述：

其中β為熱膨脹系數(shù)，c是聲速。聯(lián)立(2)～(4)式即可得到均勻生物組織中的光聲波動方程：

特別地，由于實(shí)際應(yīng)用中激光的脈寬很短，所以可以將激光的時間包絡(luò)看作是階躍函數(shù)δ(t)，此時H(r, t)=A(r)δ(t)，A(r)是光吸收系數(shù)分布函數(shù)。那么在位置r，初始時刻的t0的光聲壓強(qiáng)可以寫作：

其中，Γ=c2β/Cp稱為格日尼森系數(shù)（Grüneisen Coefficient），表示的是光能量沉積轉(zhuǎn)化為光聲波的效率。Γ是一個常數(shù)，室溫下大部分生物組織的Γ在0.11左右[13]，這說明光聲效應(yīng)產(chǎn)生的初始聲壓正比于光吸收系數(shù)分布。據(jù)此，重建生物組織光吸收分布的過程可以簡化為重建初始光聲壓強(qiáng)的過程。而重建初始光聲壓強(qiáng)的過程，則是通過采集到光聲信號后，通過特定的重建算法完成的。至此完成了光聲成像的全過程。

2 光聲成像仿真平臺的建立

仿真平臺主要可以分為四個部分：光學(xué)仿真模塊、聲學(xué)仿真模塊、圖像重建模塊以及數(shù)據(jù)輸入輸出模塊。仿真平臺的工作流程如圖1所示。

圖1 光聲仿真平臺工作流程圖Fig.1 Flow chart of the photoacoustic simulation platform

2.1 光學(xué)仿真模塊

光學(xué)仿真模塊不僅可以模擬和計算光在組織中的擴(kuò)散，同時可以設(shè)置光學(xué)仿真過程中的各項(xiàng)參數(shù)，包括光源的屬性、目標(biāo)組織和介質(zhì)的光學(xué)參數(shù)等。由于光在組織中的傳輸和擴(kuò)散屬于光聲成像的前向過程，因此正確計算并得到吸收光的光場分布至關(guān)重要。光學(xué)仿真模塊使用有限元方法（FEM）對光擴(kuò)散方程進(jìn)行求解從而得光場分布。所以在使用光學(xué)仿真模塊進(jìn)行計算前，需要先對目標(biāo)組織進(jìn)行離散化和網(wǎng)格劃分，得到標(biāo)準(zhǔn)的網(wǎng)格文件（*.msh），這一過程可以利用AutoCAD、Gmsh等專業(yè)制圖軟件來完成。值得注意的是，在有限元方法中，網(wǎng)格劃分的好壞是影響求解精度和速度的關(guān)鍵因素。網(wǎng)格劃分應(yīng)首先考慮網(wǎng)格的數(shù)目，一般來說，網(wǎng)格數(shù)量越多計算精度越高，但計算速度也就越慢，所以在確定網(wǎng)格的數(shù)目時應(yīng)該有所權(quán)衡綜合考慮。此外，網(wǎng)格劃分還有一些其它原則，例如：在硬件條件能達(dá)到要求的情況下，最好使用形狀規(guī)則的網(wǎng)格單元；在計算時數(shù)據(jù)變化梯度較大的部位（比如組織的邊界）需要采用比較密集的網(wǎng)格；對于研究中需要重點(diǎn)關(guān)注的部位（例如待重建的目標(biāo)組織）使用更為精細(xì)的網(wǎng)格；當(dāng)幾何模型的結(jié)構(gòu)和形狀對稱時，網(wǎng)格也應(yīng)劃分成對稱網(wǎng)格等。

2.2 聲學(xué)仿真模塊

光聲成像理論中，在滿足壓力和熱力學(xué)限制的條件下，利用格日尼森系數(shù)可以將組織的光吸收分布轉(zhuǎn)化為初始聲壓分布，進(jìn)而對聲在組織體內(nèi)的傳播進(jìn)行建模從而得到聲場分布。為了對聲場進(jìn)行較精確的仿真，如果繼續(xù)使用有限元方法則需要進(jìn)一步提高網(wǎng)格的數(shù)目和精度，除了會提高仿真平臺的運(yùn)行時間外，也會使仿真平臺占用更多的硬件資源。因此，聲場仿真的數(shù)值方法采用k空間偽譜法。k空間偽譜法用于模擬復(fù)雜介質(zhì)中光聲信號的傳播，和傳統(tǒng)的仿真方法，例如有限元和有限差分法相比，不僅速度快而且計算結(jié)果也較為精確[15]。利用開源的Matlab聲學(xué)仿真工具箱k-Wave，可以較方便地在Matlab上使用k空間偽譜法對光聲信號的傳播進(jìn)行計算。

由于k空間偽譜法在計算時要求網(wǎng)格必須是均勻網(wǎng)格，而光學(xué)仿真模塊中使用的有限元網(wǎng)格為非均勻網(wǎng)格，這就導(dǎo)致了兩種網(wǎng)格在網(wǎng)格形狀和節(jié)點(diǎn)數(shù)目上并不一致。所以想要在聲學(xué)仿真模塊中使用光學(xué)仿真模塊的結(jié)果，首先要進(jìn)行網(wǎng)格的統(tǒng)一[16]。相比于k空間偽譜法使用的均勻網(wǎng)格，有限元網(wǎng)格更加精細(xì)，因此需要選取均勻網(wǎng)格中每個網(wǎng)格節(jié)點(diǎn)所對應(yīng)的多個有限元網(wǎng)格的節(jié)點(diǎn)，取它們的平均值作為均勻網(wǎng)格中該網(wǎng)格節(jié)點(diǎn)的實(shí)際值。這樣就可以將光學(xué)仿真模塊的計算結(jié)果施加在聲學(xué)仿真模塊建立的模型上。

此外，聲學(xué)仿真模塊的另一個功能是定義聲場仿真中需要的各種參數(shù)，例如組織介質(zhì)和目標(biāo)組織的聲學(xué)相關(guān)屬性，不同組織區(qū)域的聲速和密度等。

2.3 圖像重建模塊

利用聲學(xué)仿真模塊采集到光聲信號后，使用圖像重建模塊及其圖形用戶界面，可以非常便捷地對采集到的時域光聲信號進(jìn)行解析法重建，例如快速傅里葉重建算法和時間反演重建 (Time Reversal Reconstruction, TR)等。此外，圖像模塊集成的重建算法還有：用于求解大型線性方程組的經(jīng)典迭代算法LSQR算法；基于L1范數(shù)正則化的稀疏重建算法YALL1算法[17]。特別地，針對實(shí)際應(yīng)用中使用廣泛的其它迭代算法，該模塊雖然不能直接進(jìn)行重建，但是可以生成并導(dǎo)出這些迭代算法在重建時所需要的系統(tǒng)測量矩陣和探測器數(shù)據(jù)以供使用。

2.4 數(shù)據(jù)輸入輸出模塊和圖形用戶界面

數(shù)據(jù)的輸入輸出也是仿真平臺中非常重要的一個部分。這部分主要包括模型和網(wǎng)格的導(dǎo)入，仿真數(shù)據(jù)的導(dǎo)出等，這些操作都可以通過圖形用戶界面來完成。與此同時，用戶圖形界面還負(fù)責(zé)對各個模塊進(jìn)行控制以及進(jìn)行交互顯示等重要工作。圖形用戶界面不僅使仿真平臺的工作過程更加直觀和簡潔，同時也極大地提高了仿真平臺的易用性和交互性，有利于平臺的進(jìn)一步開發(fā)和推廣。借助于Matlab的圖形界面工具集GUIDE，完成了對仿真參數(shù)設(shè)置、仿真過程控制、仿真結(jié)果展示以及數(shù)據(jù)導(dǎo)入導(dǎo)出等操作的圖形化和窗口化。圖2所示為仿真平臺的主界面以及光學(xué)仿真模塊的操作界面。

圖2 仿真平臺主界面以及光學(xué)仿真模塊控制界面Fig.2 Interfaces of main program and the optical module

3 仿真平臺的可靠性驗(yàn)證和重建結(jié)果分析

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

通過設(shè)計仿真實(shí)驗(yàn)，分別從前向過程、重建過程兩方面對仿真平臺的可靠性進(jìn)行驗(yàn)證。

前向過程的可靠性驗(yàn)證。使用仿真平臺的光學(xué)和聲學(xué)仿真模塊進(jìn)行光聲成像的前向過程仿真，得到光聲信號時域波形。通過對比仿真波形和理論計算波形的形狀結(jié)構(gòu)來驗(yàn)證仿真平臺前向過程的可靠性。

重建過程的可靠性驗(yàn)證。使用光聲成像仿真平臺的圖像重建模塊進(jìn)行圖像重建。通過對比重建結(jié)果和目標(biāo)物的實(shí)際位置來驗(yàn)證重建結(jié)果是否準(zhǔn)確。

3.2 前向過程的可靠性驗(yàn)證

通過對比理論計算結(jié)果和仿真結(jié)果的時域波形結(jié)構(gòu)，可以驗(yàn)證光聲仿真平臺的前向過程的仿真是否可靠。理論研究表明，使用短波脈沖光源照射球狀目標(biāo)組織產(chǎn)生的時域光聲信號具有“N”型結(jié)構(gòu)[18]。因此實(shí)驗(yàn)中首先建立球狀目標(biāo)組織的幾何模型，進(jìn)行離散化和網(wǎng)格劃分后導(dǎo)入仿真平臺進(jìn)行仿真，并將采集得到的光聲信號進(jìn)行對比。圖3分別展示了光聲波形的理論計算結(jié)果和仿真平臺的仿真結(jié)果。

圖3 時域光聲信號的理論計算結(jié)果和仿真結(jié)果對比Fig.3 Comparison of theoretical temporal waveform and simulation result of photoacoustic signal

3.3 重建過程的可靠性驗(yàn)證

圖4 所示為驗(yàn)證實(shí)驗(yàn)中使用的仿體幾何模型和離散后的網(wǎng)格模型。仿體分為兩部分，大圓柱體為組織介質(zhì)，其底面半徑為15 mm，高為60 mm；內(nèi)部的小圓柱體為目標(biāo)物，底面半徑為4 mm，高為6 mm。目標(biāo)物和介質(zhì)的光吸收系數(shù)的比值設(shè)置為20:1。共16個點(diǎn)光源均勻分布在以大圓柱體中心為球心，半徑為30 mm的球體表面。32個超聲探測器組成超聲探測陣列，均勻設(shè)置在大圓柱體的表面。為了直觀顯示目標(biāo)物的重建結(jié)果和真實(shí)位置，選取了了z=30 mm平面的重建結(jié)果，如圖5所示。

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

圖3展示了一個周期內(nèi)光聲信號時域波形的仿真值和理論值。其中圖中虛線為理論計算結(jié)果，實(shí)線為仿真平臺得到的仿真值。當(dāng)目標(biāo)物體為實(shí)心球狀物時，仿真得到的光聲信號波形和理論計算值基本重合，不但呈現(xiàn)出了“N”型的結(jié)構(gòu)，而且從信號出現(xiàn)到消失的持續(xù)時間也均為4 μs左右。這驗(yàn)證了仿真平臺前向過程仿真的可靠性。

圖4 實(shí)驗(yàn)?zāi)Ｐ虵ig.4 Experimental model

圖5 重建結(jié)果和目標(biāo)物的真實(shí)位置（z=30 mm平面）Fig.5 Reconstructed results and real position of target (plane z=30mm)

圖5 展示的是重建結(jié)果和目標(biāo)物的真實(shí)位置。圖中的亮團(tuán)是仿真平臺的重建結(jié)果，目標(biāo)物的真實(shí)位置用實(shí)線進(jìn)行了標(biāo)注。重建過程將大圓柱體離散成為長方體，所以仿體的真實(shí)邊界在圖中也以實(shí)線標(biāo)出。由圖5可知，重建出的亮團(tuán)雖然存在著一定程度的偽影，但是整體呈現(xiàn)出較為規(guī)則的圓形形狀，這和目標(biāo)物本身的形狀相符。此外，亮團(tuán)最亮的部位和目標(biāo)物的真實(shí)位置基本重合，這說明仿真平臺的重建結(jié)果是準(zhǔn)確的，因此證明了仿真平臺后向過程的可靠性。

4 總結(jié)與展望

本文闡述了基于Matlab的光聲成像仿真平臺的設(shè)計和實(shí)現(xiàn)。對圖形用戶界面以及各模塊的工作原理進(jìn)行了介紹。通過仿真實(shí)驗(yàn)，從前向過程和后向過程兩個方面證明了仿真平臺的可靠性，同時提供了利用該平臺進(jìn)行光聲成像仿真和重建的簡單演示。與使用COMSOL，ANSYS等商業(yè)軟件相比，該仿真平臺可以大大降低仿真的成本，在非復(fù)雜實(shí)驗(yàn)的場景下提高仿真效率。同時，該仿真平臺的開源特性也使其具有更好的可移植性和可維護(hù)性。為光聲成像技術(shù)的理論研究提供了一個較為可靠和高效的研究工具。但是距離成為一個成熟的仿真工具，該仿真平臺仍需進(jìn)一步地改進(jìn)和完善。針對該仿真平臺的改進(jìn)工作將會主要集中在以下幾個方面：①集成幾何建模和網(wǎng)格離散的功能，使仿真平臺更具獨(dú)立性；②提供更多的可設(shè)置系統(tǒng)參數(shù)，進(jìn)一步細(xì)化仿真過程以應(yīng)對更加復(fù)雜的仿真場景；③為不同類型的真實(shí)實(shí)驗(yàn)提供更加豐富的重建算法并提高重建精度。