k-wave偽譜法在生物醫學超聲仿真中的研究進展

陳思瑤，杜宇翎，周 野，李發琪，李成海*

（1.重慶醫科大學 生物醫學工程學院 超聲醫學工程國家重點實驗室，重慶 400016；2.重慶市生物醫學工程重點實驗室，重慶 400016）

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

1 k-wave聲學工具箱的發展

k-wave是基于MATLAB和C++編寫的開源聲學工具箱，由Bradley Treeby、Ben Cox和Jiri Jaros開發。該工具箱使用Westervelt方程的廣義模型來模擬超聲波在軟組織中的傳播，考慮了聲傳播的非線性、組織的不均勻性和冪律吸收?？梢詫崿F復雜多層生物組織介質中的時域聲學和超聲模擬。k-wave工具箱最初在2009年由倫敦大學光聲成像小組開發，主要解決無損介質[13]中光聲波場的模擬和重建問題。在后續的版本中擴展了包括時變壓力和速度源、聲吸收、非線性、彈性材料和超聲波換能器模型的功能。

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

k-space計算聲場的物理模型是基于運動方程、連續性方程以及物態方程構建的Westervelt方程，運動方程、連續性方程以及物態方程的描述如下：

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

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

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

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

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

2 k-wave在超聲仿真應用中的研究

k-wave作為一種開源的聲學工具箱，是實現生物組織的聲傳播仿真的重要工具，目前國內外已報道基于kwave開展了相關的研究，應用領域主要包括聲波在非均勻生物組織介質中傳播的數值模擬、光聲圖像重建以及無損探傷等。

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

為了精確求解光聲反演問題，需要建立光致聲波在生物軟組織中產生和傳播的數值模型。有大量的文獻描述了超聲在生物軟組織中的傳播，部分理論可直接適用于光聲學，同時開發了相應的方法來數值求解偏微分方程。其中基于k-space偽譜法的數值仿真具有較大優勢[28]。有限元（FE）和有限差分（FD）建模技術可以用于建模光致聲波的傳播，但這兩種方法精確的仿真需要網格單元的尺寸約為十分之一波長。對于較高頻率條件下的數值仿真，這一要求加大了對計算硬件和算法的要求。COX等人[29]利用kwave進行建模研究了二維傳感器元件記錄表面平均壓力值的理論方向性，然后針對時間反轉圖像重建，研究了傳感器方向對光聲層析成像的影響。研究發現即使在使用完整連續的測量曲面時，傳感器指向性的加入也會帶來圖像重建的假象，研究明晰了傳感器元件陣列的使用對重建光聲圖像的結構保真度和定量精度的影響。

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

3 結束語

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