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

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

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

光聲成像技術(shù)是一種新型無(wú)創(chuàng)在體生物醫(yī)學(xué)成像技術(shù)。光聲成像結(jié)合了超聲和光學(xué)成像的優(yōu)點(diǎn)[1]，兼具超聲成像的高分辨率和純光學(xué)成像的高對(duì)比度特性，因而成為了生物醫(yī)學(xué)成像領(lǐng)域的研究前沿和熱點(diǎn)之一[2]。光聲成像的基本原理是光聲效應(yīng)[3]：使用短脈沖激光照射生物組織，組織體由于吸收能量導(dǎo)致的瞬時(shí)熱膨脹進(jìn)而輻射超聲波。通過檢測(cè)和采集這些攜帶生物組織特性的超聲信號(hào)，就可以利用一定的重建算法重建出生物組織內(nèi)部的光學(xué)吸收分布圖像。由于生物組織的光吸收特性與組織的生理特征、結(jié)構(gòu)形態(tài)特征、病變特性、代謝狀態(tài)甚至神經(jīng)活動(dò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í)實(shí)驗(yàn)的系統(tǒng)噪聲和測(cè)量誤差，為理論研究提供有利條件。因此，光聲成像仿真技術(shù)對(duì)于光聲成像的理論研究及其在實(shí)際中的應(yīng)用有重要的指導(dǎo)意義，也是驗(yàn)證實(shí)驗(yàn)設(shè)計(jì)性能和算法性能的首選方案。

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

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

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

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

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

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

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

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

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

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

2 光聲成像仿真平臺(tái)的建立

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

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

2.1 光學(xué)仿真模塊

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

2.2 聲學(xué)仿真模塊

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

由于k空間偽譜法在計(jì)算時(shí)要求網(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)格中每個(gè)網(wǎng)格節(jié)點(diǎn)所對(duì)應(yīng)的多個(gè)有限元網(wǎng)格的節(jié)點(diǎn)，取它們的平均值作為均勻網(wǎng)格中該網(wǎng)格節(jié)點(diǎn)的實(shí)際值。這樣就可以將光學(xué)仿真模塊的計(jì)算結(jié)果施加在聲學(xué)仿真模塊建立的模型上。

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

2.3 圖像重建模塊

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

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

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

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

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

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

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

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

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

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

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

圖3 時(shí)域光聲信號(hào)的理論計(jì)算結(jié)果和仿真結(jié)果對(duì)比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個(gè)點(diǎn)光源均勻分布在以大圓柱體中心為球心，半徑為30 mm的球體表面。32個(gè)超聲探測(cè)器組成超聲探測(cè)陣列，均勻設(shè)置在大圓柱體的表面。為了直觀顯示目標(biāo)物的重建結(jié)果和真實(shí)位置，選取了了z=30 mm平面的重建結(jié)果，如圖5所示。

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

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

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

4 總結(jié)與展望

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