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

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

（1.重慶醫(yī)科大學(xué) 生物醫(yī)學(xué)工程學(xué)院 超聲醫(yī)學(xué)工程國家重點(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]。有限差分或有限元方法常被用來求解控制方程，但由于每個(gè)最小波長需要大量的網(wǎng)格點(diǎn)來避免數(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方程的廣義模型來模擬超聲波在軟組織中的傳播，考慮了聲傳播的非線性、組織的不均勻性和冪律吸收。可以實(shí)現(xiàn)復(fù)雜多層生物組織介質(zhì)中的時(shí)域聲學(xué)和超聲模擬。k-wave工具箱最初在2009年由倫敦大學(xué)光聲成像小組開發(fā)，主要解決無損介質(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)的信息來求解波的傳播方程。k-wave工具箱的控制方程采用k-space偽譜方法求解，其中傅里葉配置譜方法用于計(jì)算空間梯度，k-space校正有限差分格式用于在時(shí)間上向前積分，具有快速、簡(jiǎn)便的特點(diǎn)。由于其全局性質(zhì)，譜方法可以比局部方法更精確，例如，應(yīng)用于周期問題的偽譜方法已被證明與無限階有限差分方法等價(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í)間步長[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)的基本概述和許多易于遵循的教程示例，以說明工具箱的功能；與基于時(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)生物組織的聲傳播仿真的重要工具，目前國內(nèi)外已報(bào)道基于kwave開展了相關(guān)的研究，應(yīng)用領(lǐng)域主要包括聲波在非均勻生物組織介質(zhì)中傳播的數(shù)值模擬、光聲圖像重建以及無損探傷等。

聲波在非均勻生物組織介質(zhì)中傳播的數(shù)值模擬對(duì)開展生物醫(yī)學(xué)超聲研究具有重要意義。例如，在超聲治療應(yīng)用中，數(shù)值模擬可用于研究磁共振引導(dǎo)的聚焦超聲手術(shù)中的相位像差[19-20]，并改善治療結(jié)果。對(duì)于診斷超聲，數(shù)值建模已被用作圖像重建以及理解超聲成像中圖像退化來源的重要工具[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)的方法來數(shù)值求解偏微分方程。其中基于k-space偽譜法的數(shù)值仿真具有較大優(yōu)勢(shì)[28]。有限元（FE）和有限差分（FD）建模技術(shù)可以用于建模光致聲波的傳播，但這兩種方法精確的仿真需要網(wǎng)格單元的尺寸約為十分之一波長。對(duì)于較高頻率條件下的數(shù)值仿真，這一要求加大了對(duì)計(jì)算硬件和算法的要求。COX等人[29]利用kwave進(jìn)行建模研究了二維傳感器元件記錄表面平均壓力值的理論方向性，然后針對(duì)時(shí)間反轉(zhuǎn)圖像重建，研究了傳感器方向?qū)饴晫游龀上竦挠绊憽Ｑ芯堪l(fā)現(xiàn)即使在使用完整連續(xù)的測(cè)量曲面時(shí)，傳感器指向性的加入也會(huì)帶來圖像重建的假象，研究明晰了傳感器元件陣列的使用對(duì)重建光聲圖像的結(jié)構(gòu)保真度和定量精度的影響。

除此之外，Prieur等人[30]利用k-wave工具箱模擬橫波彈性成像，將均勻和各向同性介質(zhì)中的模擬結(jié)果與解析解和有限元建模仿真的結(jié)果進(jìn)行比較。從k-wave獲得的剪切位移與解析解非常相似，在短傳播時(shí)間和大傳播時(shí)間下，聚焦區(qū)域周圍的均方根誤差分別低于9%和21%，說明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é)束語

圍繞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à)值。但目前國內(nèi)在這方面的研究尚處于起步階段，有待進(jìn)一步拓展其研究。