k-wave偽譜法在生物醫(yī)學(xué)超聲仿真中的研究進(jìn)展

陳思瑤，杜宇翎，周 野，李發(fā)琪，李成海*

（1.重慶醫(yī)科大學(xué) 生物醫(yī)學(xué)工程學(xué)院 超聲醫(yī)學(xué)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，重慶 400016；2.重慶市生物醫(yī)學(xué)工程重點(diǎn)實(shí)驗(yàn)室，重慶 400016）

超聲在生物組織中的非線性傳播模擬對(duì)于許多應(yīng)用具有重要作用，包括超聲波換能器的設(shè)計(jì)，聚焦超聲治療計(jì)劃的制定，診斷超聲成像等[1-2]。然而，模擬非線性波的傳播是一項(xiàng)計(jì)算困難的任務(wù)，因?yàn)榻橘|(zhì)離散必須足夠精細(xì)，以捕捉在超聲波傳播時(shí)產(chǎn)生的諧波。對(duì)于強(qiáng)激波，可能會(huì)產(chǎn)生幾百次諧波，因此需要非常密集的計(jì)算網(wǎng)格[3]。目前，多數(shù)介質(zhì)中的非線性超聲模擬主要基于Khokhlov-Zabolotskaya-Kuznetsov（KZK）方程[4]或Burgers方程[5]。雖然這些模型在許多實(shí)際情況下是準(zhǔn)確的[6]，但它們僅限于模擬定向聲束，并且通常只考慮均勻介質(zhì)中的單向傳播[7]。有限差分或有限元方法常被用來(lái)求解控制方程，但由于每個(gè)最小波長(zhǎng)需要大量的網(wǎng)格點(diǎn)來(lái)避免數(shù)值色散[8]，增加了計(jì)算的復(fù)雜程度。為了減少這種計(jì)算負(fù)擔(dān)，基于傅里葉偽譜方法的非線性超聲仿真模型被提出[9-10]。與有限差分方法相比，這種方法可以減少所需網(wǎng)格點(diǎn)的總數(shù)，最多可減小兩個(gè)數(shù)量級(jí)[11-12]。基于k-space偽譜方法的全波非線性聲學(xué)模擬作為開源k-wave聲學(xué)工具箱的一部分[11，13]。該方法可用于描述聲波在非均勻介質(zhì)中的聲傳播，考慮了冪律吸收，對(duì)波的方向性沒有限制。本文將從k-wave聲學(xué)工具箱的發(fā)展以及k-wave在超聲仿真應(yīng)用中的研究進(jìn)行綜述，為超聲在生物組織中非線性傳播的模擬研究提供參考。

1 k-wave聲學(xué)工具箱的發(fā)展

k-wave是基于MATLAB和C++編寫的開源聲學(xué)工具箱，由Bradley Treeby、Ben Cox和Jiri Jaros開發(fā)。該工具箱使用Westervelt方程的廣義模型來(lái)模擬超聲波在軟組織中的傳播，考慮了聲傳播的非線性、組織的不均勻性和冪律吸收。可以實(shí)現(xiàn)復(fù)雜多層生物組織介質(zhì)中的時(shí)域聲學(xué)和超聲模擬。k-wave工具箱最初在2009年由倫敦大學(xué)光聲成像小組開發(fā)，主要解決無(wú)損介質(zhì)[13]中光聲波場(chǎng)的模擬和重建問題。在后續(xù)的版本中擴(kuò)展了包括時(shí)變壓力和速度源、聲吸收、非線性、彈性材料和超聲波換能器模型的功能。

求解非線性Westervelt方程常見的方法為有限差分法、有限元法和譜法[14]。有限差分法和有限元法被稱為局部法，因?yàn)樗P(guān)注的波動(dòng)傳播方程僅根據(jù)附近點(diǎn)的條件在每一點(diǎn)求解。相比而言，譜方法，如k-space方法[15]和偽譜方法[16-17]，被認(rèn)為是全局的，因?yàn)樵诿恳稽c(diǎn)上都是利用整個(gè)波場(chǎng)的信息來(lái)求解波的傳播方程。k-wave工具箱的控制方程采用k-space偽譜方法求解，其中傅里葉配置譜方法用于計(jì)算空間梯度，k-space校正有限差分格式用于在時(shí)間上向前積分，具有快速、簡(jiǎn)便的特點(diǎn)。由于其全局性質(zhì)，譜方法可以比局部方法更精確，例如，應(yīng)用于周期問題的偽譜方法已被證明與無(wú)限階有限差分方法等價(jià)[18]。非線性k-space仿真的優(yōu)勢(shì)已經(jīng)由臨床超聲換能器在均勻和非均勻介質(zhì)中的波束方向圖的三維模擬進(jìn)行驗(yàn)證。通過將3D中的完整空間域離散化，并以時(shí)間步進(jìn)的方式計(jì)算域內(nèi)各處的壓力和粒子速度場(chǎng)的解。與以前基于KZK方程的超聲模擬相比，此方法對(duì)聲波的方向性或空間變化沒有任何限制。與有限差分法和有限元法相比，這種方法在相同精度下允許更粗的網(wǎng)格間距和更大的時(shí)間步長(zhǎng)[11-12]。

k-space計(jì)算聲場(chǎng)的物理模型是基于運(yùn)動(dòng)方程、連續(xù)性方程以及物態(tài)方程構(gòu)建的Westervelt方程，運(yùn)動(dòng)方程、連續(xù)性方程以及物態(tài)方程的描述如下：

u是聲粒子速度，p是聲壓，ρ是聲學(xué)介質(zhì)密度，ρ0是環(huán)境（或平衡）介質(zhì)密度，C0是等熵的聲波的速度，d為聲粒子位移。這里的B/A是表征有限振幅效應(yīng)對(duì)聲速的相對(duì)貢獻(xiàn)的非線性參數(shù)。與線性情況相比，質(zhì)量守恒方程包括一個(gè)附加項(xiàng)，它解釋了一個(gè)對(duì)流非線性，其中粒子速度對(duì)波速有貢獻(xiàn)。

壓力-密度關(guān)系中的L算子是一個(gè)線性積分-微分算子，它解釋了聲波的吸收和擴(kuò)散，遵循頻率冪律：

這里τ、η分別是吸收和色散比例系數(shù)：

其中α0是冪律吸收因子，單位為Np（rad/s）-ym-1，y是冪律指數(shù)。該方程考慮了二階聲學(xué)非線性、冪律吸聲和材料特性（聲速、密度、非線性和吸收系數(shù)）的非均勻分布。與傳統(tǒng)的有限差分方法相比，這提高了梯度計(jì)算的精度，從而降低了對(duì)精細(xì)計(jì)算網(wǎng)格的要求。

k-wave聲學(xué)工具箱包括了以下幾點(diǎn)優(yōu)勢(shì)：提供了較為詳細(xì)的用戶手冊(cè)，包括對(duì)控制方程和數(shù)值方法的一般介紹，還提供了軟件架構(gòu)的基本概述和許多易于遵循的教程示例，以說(shuō)明工具箱的功能；與基于時(shí)域有限差分（FDTD）格式的模型相比，k-wave數(shù)值模型的主要優(yōu)點(diǎn)是精確模擬所需的空間和時(shí)間網(wǎng)格點(diǎn)較少，這意味著這些模型運(yùn)行得更快，使用的內(nèi)存更少；能夠模擬壓力和速度源，包括光聲源，以及診斷和治療超聲換能器；可以使用定向元素指定任意檢測(cè)表面的功能，并可以選擇記錄聲壓，粒子速度和聲強(qiáng)；優(yōu)化的C++程序代碼，最大限度地提高了大型模擬的計(jì)算性能；選擇使用前向模型作為一個(gè)靈活的時(shí)間反轉(zhuǎn)圖像重建算法的光聲層析與任意測(cè)量表面；能快速進(jìn)行光聲圖像重建算法，用于記錄在線性（2D）或平面（3D）測(cè)量面上的數(shù)據(jù)。

2 k-wave在超聲仿真應(yīng)用中的研究

k-wave作為一種開源的聲學(xué)工具箱，是實(shí)現(xiàn)生物組織的聲傳播仿真的重要工具，目前國(guó)內(nèi)外已報(bào)道基于kwave開展了相關(guān)的研究，應(yīng)用領(lǐng)域主要包括聲波在非均勻生物組織介質(zhì)中傳播的數(shù)值模擬、光聲圖像重建以及無(wú)損探傷等。

聲波在非均勻生物組織介質(zhì)中傳播的數(shù)值模擬對(duì)開展生物醫(yī)學(xué)超聲研究具有重要意義。例如，在超聲治療應(yīng)用中，數(shù)值模擬可用于研究磁共振引導(dǎo)的聚焦超聲手術(shù)中的相位像差[19-20]，并改善治療結(jié)果。對(duì)于診斷超聲，數(shù)值建模已被用作圖像重建以及理解超聲成像中圖像退化來(lái)源的重要工具[21-23]。Suomi等[24]人使用的開源k-wave聲學(xué)模擬工具箱，基于三個(gè)不同患者分割的三維CT數(shù)據(jù)集預(yù)測(cè)高強(qiáng)度聚焦超聲（High Intensity Focussed Ultrasound，HIFU）治療劑量，模擬了三位患者腎臟前方組織層的衰減、反射和折射的綜合效應(yīng)如何影響超聲場(chǎng)的強(qiáng)度和形狀。模擬結(jié)果表明，超聲場(chǎng)的強(qiáng)度平均下降了11.1 dB，并且發(fā)現(xiàn)強(qiáng)度損失可以在衰減和折射之間大致平均分配，為臨床HIFU治療腎癌提供了一定的參考。Kittiphot等人[25]基于k-wave工具包模擬了高強(qiáng)度聚焦超聲熱療在腫瘤組織中療效，模擬結(jié)果可見，聚焦區(qū)域周圍聲壓的差異隨著聚焦深度的增加而增大。進(jìn)一步使用數(shù)據(jù)可視化裝置，獲得相應(yīng)的溫度分布2D圖像，以顯示乳腺癌的消融療效，觀察到除了焦點(diǎn)區(qū)域的損傷，周圍的健康組織未受影響。此外，姜翔飛等人[26]利用此工具包，對(duì)基于時(shí)間反演的經(jīng)顱二維超聲進(jìn)行仿真，研究了不同換能器陣列數(shù)和不同的超聲發(fā)射頻率對(duì)經(jīng)顱超聲聚焦效果的影響。張值豪[27]也利用k-wave工具箱對(duì)超聲斷層成像透射過程進(jìn)行仿真，在環(huán)形探頭上實(shí)現(xiàn)了對(duì)平行束以及扇形束重建效果的仿真。

為了精確求解光聲反演問題，需要建立光致聲波在生物軟組織中產(chǎn)生和傳播的數(shù)值模型。有大量的文獻(xiàn)描述了超聲在生物軟組織中的傳播，部分理論可直接適用于光聲學(xué)，同時(shí)開發(fā)了相應(yīng)的方法來(lái)數(shù)值求解偏微分方程。其中基于k-space偽譜法的數(shù)值仿真具有較大優(yōu)勢(shì)[28]。有限元（FE）和有限差分（FD）建模技術(shù)可以用于建模光致聲波的傳播，但這兩種方法精確的仿真需要網(wǎng)格單元的尺寸約為十分之一波長(zhǎng)。對(duì)于較高頻率條件下的數(shù)值仿真，這一要求加大了對(duì)計(jì)算硬件和算法的要求。COX等人[29]利用kwave進(jìn)行建模研究了二維傳感器元件記錄表面平均壓力值的理論方向性，然后針對(duì)時(shí)間反轉(zhuǎn)圖像重建，研究了傳感器方向?qū)饴晫游龀上竦挠绊?。研究發(fā)現(xiàn)即使在使用完整連續(xù)的測(cè)量曲面時(shí)，傳感器指向性的加入也會(huì)帶來(lái)圖像重建的假象，研究明晰了傳感器元件陣列的使用對(duì)重建光聲圖像的結(jié)構(gòu)保真度和定量精度的影響。

除此之外，Prieur等人[30]利用k-wave工具箱模擬橫波彈性成像，將均勻和各向同性介質(zhì)中的模擬結(jié)果與解析解和有限元建模仿真的結(jié)果進(jìn)行比較。從k-wave獲得的剪切位移與解析解非常相似，在短傳播時(shí)間和大傳播時(shí)間下，聚焦區(qū)域周圍的均方根誤差分別低于9%和21%，說(shuō)明k-wave是模擬聲輻射力和剪切波傳播的一種有效工具。另外，Gu等人[31]模擬聲波在強(qiáng)非均勻介質(zhì)中傳播時(shí)，使用k-wave的結(jié)果用作比較和驗(yàn)證的基準(zhǔn)，系統(tǒng)地評(píng)估了改進(jìn)混合域方法。吳其洲等人[32]在圓柱體缺陷構(gòu)件的超聲探傷重構(gòu)研究中，利用k-wave工具箱建立仿真平臺(tái)，在圓柱體周圍設(shè)置360個(gè)傳感器，同時(shí)接受散射信號(hào)，實(shí)現(xiàn)缺陷信號(hào)的采集，為缺陷重構(gòu)提供數(shù)據(jù)支撐。孫明健等人[33]在基于光聲信號(hào)的高鐵鋼軌表面缺陷檢測(cè)方法中，使用有限元及k-wave方法對(duì)鋼軌表面的光聲圖像進(jìn)行了重建，結(jié)果表明，采用這種方法能較好的檢測(cè)到表面微裂紋，此結(jié)果在鋼軌探傷領(lǐng)域有較大的可行性及發(fā)展?jié)摿?。張思思等人[34]基于k-wave建立了表面缺陷的超聲檢測(cè)模型，結(jié)果為表面缺陷檢測(cè)、數(shù)值模擬、超聲回波信號(hào)處理以及特征提取等技術(shù)提供了理論依據(jù)。綜上可見，基于k-wave進(jìn)行的聲傳播模擬在較多的領(lǐng)域都具有重要的應(yīng)用價(jià)值。

3 結(jié)束語(yǔ)

圍繞k-wave偽譜法在生物醫(yī)學(xué)超聲仿真中的應(yīng)用，本文首先介紹了k-wave聲學(xué)工具箱的發(fā)展，進(jìn)一步綜述了k-wave在聲學(xué)仿真應(yīng)用中的相關(guān)研究。通過綜述可見，基于k-wave偽譜法對(duì)于模擬不均勻多層組織中的聲傳播具有重要應(yīng)用價(jià)值。但目前國(guó)內(nèi)在這方面的研究尚處于起步階段，有待進(jìn)一步拓展其研究。