利用k-Wave計算超聲在腹壁組織的聲輻射力?

喬玉配 宮門陽 汪海賓 李楠 劉杰惠 毛一葳 何愛軍 劉曉宙

(1 近代聲學教育部重點實驗室 南京大學聲學研究所 人工微結構科學與技術協(xié)同創(chuàng)新中心 南京 210093)

(2 江蘇科技大學理學院 鎮(zhèn)江 212003)

(3 南京大學電子科學與工程學院 南京 210093)

0 引言

聲輻射力一般定義為聲波對介質(zhì)施加的周期平均力[1]。聲輻射力的研究始于1902年Rayleigh提出“聲輻射壓”的概念[2]，已經(jīng)有百年歷史[3]。1991年，Wu[4]首先提出了聲鑷子的概念并進行了實驗，實現(xiàn)了利用聲輻射力對物體的操控。利用聲輻射力不僅可以操控生物細胞和粒子，也可以用于探索內(nèi)部解剖結構并獲取診斷信息。聲波屬于機械波，利用聲輻射力進行的操控對介質(zhì)的導電性、透光性等沒有特殊的要求，可以實現(xiàn)對粒子和生物組織的無創(chuàng)、非接觸、無標記、多功能性等操控[5?7]，且其相應的設備簡單易集成和微型化，因此，聲輻射力在精密制造、精準醫(yī)療醫(yī)學診斷，評估生物組織和液體的黏彈性特性、彈性成像等生物醫(yī)學領域具有廣泛的應用，尤其是在彈性成像領域[8?12]。彈性成像能夠獲得生物組織的彈性信息，可以檢測出生物體內(nèi)部與周圍組織彈性不同的區(qū)域，進而根據(jù)此差異判斷相應組織或器官可能發(fā)生的病理改變以及其位置、形狀和大小。聲輻射力彈性成像是一種新的無創(chuàng)彈性成像方法，是利用聲波產(chǎn)生的聲輻射力對組織施加壓，通過采集組織被激發(fā)前后的超聲圖像進行運動估計，或檢測組織在聲輻射力終止后不同時間點的應變，或檢測聲輻射力激發(fā)組織而產(chǎn)生的剪切波的傳播進行成像來反映生物組織的彈性(或硬度)[13]。因此，聲輻射力是彈性成像技術有效性和精確性的關鍵核心，對其預測和研究至關重要。

數(shù)值模擬是計算聲輻射力最直接和直觀的方法。k-Wave可以用于復雜和真實組織介質(zhì)中的時域聲學和超聲模擬，是模擬聲輻射力的一種精確而有效的工具。它結合了有限差分法的簡單性和在任何非均勻介質(zhì)中模擬的靈活性，具有快速、易用的特點[14?17]。本文基于腹壁組織圖像，利用k-Wave對超聲波在腹壁組織區(qū)域傳播時的聲場進行數(shù)值模擬，通過模擬得到其聲場分布，進而計算求得其聲輻射力的分布。

1 理論

聲波在介質(zhì)中傳播時，壓力、密度、溫度、顆粒速度等都會發(fā)生動態(tài)變化，這些變化可以用一系列基于介質(zhì)內(nèi)質(zhì)量、動量和能量守恒的一階耦合偏微分方程來描述。當能量足夠高的聲波在各向異性介質(zhì)中傳播時，聲波的傳播不再是線性的，其一階耦合方程組為[17]

其中，u是介質(zhì)質(zhì)點速度，ρ0是介質(zhì)靜態(tài)密度，p為聲壓，ρ是密度，是一個與空間有關的物理量，t為時間，c0是等熵聲速，d介質(zhì)質(zhì)點位移，B/A為非線性參數(shù)，L說明了遵循頻率冪律的聲學吸收和色散，可以表示[18]其中，τ、η為吸收和色散比例系數(shù)，y1為冪律指數(shù)。注意，當這些公式作為耦合方程被求解時，u·?ρ0和d·?ρ0被取消[19]，因此，為了提高計算效率，這些項不包括在下面給出的離散方程中。

為了減少精確模擬所需的內(nèi)存和時間步數(shù)，k-Wave使用k-space偽譜法離散化耦合聲學方程組(1)，空間導數(shù)的頻譜計算使用快速傅里葉變換，每個聲波波長只需要兩個網(wǎng)格節(jié)點即可達到可接受的精度，能夠快速高效地實現(xiàn)生物組織中非線性超聲傳播的模擬，計算得到其聲輻射力。k-Wave中基于k-space偽譜法離散化耦合聲學方程組(1)得到其離散形式如下：

式中的Nξ是在ξ方向的網(wǎng)格點數(shù)。

方程(3)中的離散方程使用時間步長?t=CFL?x/cmax迭代求解，為了能夠平衡準確性和計算效率之間的關系，CFL≤c0/cmax。在弱異構介質(zhì)，通常CFL=0.3能夠為精度和計算速度之間提供很好的平衡[19]。在每個時間步長，可以通過在計算域內(nèi)的適當網(wǎng)格點添加源值來引入質(zhì)量或力源。類似地，模擬的輸出可以通過在特定網(wǎng)格點的每個時間步長記錄聲學變量來獲得。

聲輻射力的一般表達式為[14?15,20?21]

式(4)用到擾動分析以及

其中，ρ、p、v分別表示介質(zhì)的密度、聲壓、質(zhì)點速度，n為面S的法向矢量，〈〉表示時間平均，下標0、1、2分別代表靜態(tài)、一階、二階微小量。

對于體積元(N/m3)的聲輻射力，公式(4)可以表示為

雖然公式(4)和公式(5)是無黏性流體情況得到的，嚴格地說只適用于在不存在剪切且剪切模量為零的流體中的傳播，但它考慮了由黏性衰減和不均勻性可能造成的散射而導致的波動量在封閉面內(nèi)體積沉積，而且在醫(yī)學超聲應用的實際環(huán)境中，軟組織具有相當?shù)偷奈蘸蜕⑸涮匦裕M織中的體模量遠大于剪切模量，壓縮應力比剪切應力更容易出現(xiàn)，因此，公式(4)和公式(5)計算腹壁組織的聲輻射力是合理的，對本文的問題仍然適用。基于公式(5)計算的聲輻射力所用的聲壓和質(zhì)點振速是由k-Wave運算過程中記錄的聲壓和非交錯網(wǎng)格中的質(zhì)點振速。

2 仿真模型的建立

本文的聲場是由100個陣元組成的線性相控陣換能器，如圖1(a)所示；其相應的物理參數(shù)示意圖如圖1(b)所示。仿真中設置計算區(qū)域為(Nx×Ny×Nz)，沿x方向的網(wǎng)格數(shù)Nx=88，沿y方向的網(wǎng)格數(shù)Ny=108，沿z方向的網(wǎng)格數(shù)Nz=44，陣元寬度為1個網(wǎng)格，陣元長度為12個網(wǎng)格，陣元中心距離為0，換能器的中心頻率為1MHz，換能器放置在yz平面中間位置，聲波沿x方向傳播，其聚焦點在x方向距離換能器Nx/2×dx。為了能夠保證較高的精度和計算速度取網(wǎng)格大小dx=dy=dz=1.5×10?4m和CFL=0.1。生物組織樣品(腹壁組織)放置在聲場中，其橫截面圖如圖2所示，包含結締組織、肌肉、脂肪、水等主要成分，其相應的參數(shù)值如表1所示。圖2中的Lx=Nx×dx，Ly=Ny×dy，Lz=Nz×dz。另外在仿真計算時在計算區(qū)域的周圍添加PML邊界層來模擬無限大區(qū)域。

圖1 線陣換能器示意圖Fig.1 Schematic of the linear array transducer

圖2 組織橫截面圖[22]Fig.2 Cross-sectional tissue map[22]

表1 組織中各成分對應的參數(shù)[22]Table 1 Parameters of comp onents in tissue[22]

3 數(shù)值仿真

3.1 有效性驗證

為了驗證k-Wave計算的正確性，首先利用公式(5)對放置在平面波聲場中液體柱的聲輻射力進行了仿真計算，并與已有的結果進行對比，液體柱內(nèi)聲速為930 m/s2，其密度為1000 kg/m3[23]，半徑為7.5 mm，其聲場情況如圖3所示。圖3(a)為平面波入射示意圖，圖中的圓表示柱形粒子，從圖中可以看出平面波向前傳播，遇到液體柱發(fā)生散射、折射等現(xiàn)象，從而使液體柱受到聲輻射力的作用。圖3(b)展示了時間t=40μs時其聲場情況，圖3(c)則是轉(zhuǎn)化成頻域中其聲場情況。平面入射頻率為0.5 MHz，其余參數(shù)設置和文獻[24]一致。通過k-Wave計算所得無量綱化聲輻射力為0.0039，理論計算為0.0042，誤差為7.1%，在可接受范圍內(nèi)，驗證了k-Wave計算的正確性。

圖3 平面波入射到液體柱時的聲場分布Fig.3 Sound field distribution when a plane wave incident on a liquid cylinder

3.2 腹壁中的聲輻射力

為了計算生物組織(腹壁)的聲輻射力，首先用k-Wave對其聲場進行仿真，記錄整個空間中的聲壓以及聲速，之后用來計算其聲輻射力。仿真模型以及計算中所需信息如第2節(jié)所述。仿真得到該組織的歸一化聲壓振幅空間分布如圖4(a)～圖4(c)所示。圖4(a)給出了換能器中心頻率分別為1 MHz、3 MHz時，其歸一化聲壓振幅隨聲波主傳播軸x的變化；圖4(b)～圖4(c)為換能器中心頻率為1 MHz、3 MHz時，xy平面內(nèi)歸一化聲壓振幅分布。另外，為了分析腹壁組織的非均勻性對聲場的影響，在圖4(d)～圖4(f)中給出了均勻介質(zhì)的歸一化聲壓振幅空間分布作為對比。圖4(d)為換能器中心頻率為1 MHz、3 MHz時，均勻介質(zhì)中歸一化聲壓振幅隨聲波主傳播軸x的變化；圖4(e)～圖4(f)為換能器中心頻率為1 MHz、3 MHz時，均勻介質(zhì)中xy平面內(nèi)歸一化聲壓振幅分布。從圖4中可以看出，聲壓振幅在聚焦點的位置達到最大值。對比腹壁組織和均勻介質(zhì)中聲場情況可以看出，聲波在非均勻介質(zhì)中傳播，由于介質(zhì)密度和聲速等的變化而產(chǎn)生散射、反射等現(xiàn)象。對于同一種組織、不同的頻率激發(fā)下，產(chǎn)生的聲場不同。

圖4 腹壁和均勻介質(zhì)中歸一化聲壓振幅空間分布Fig.4 Spatial distribution of the normalized pressure amplitude in abdominal wall tissue and homogeneous media

在計算聲輻射力時，得到了在三維空間網(wǎng)格每個點上的力的3個時變分量，然而，在實際的應用中，僅使用零仰角平面(xy平面)內(nèi)的力作為輸入。該平面是力最大的平面，也是用于實際建立彈性圖[15]的平面。因此，基于公式(5)計算得到其聲輻射力如圖5所示，圖中聲輻射力采用最大聲輻射力進行了歸一化。

圖5 腹壁組織中歸一化聲輻射力Fig.5 The normalized acoustic radiation force in abdominal tissue

圖5分別為歸一化的軸向聲輻射力、橫向聲輻射力以及總聲輻射力在xy平面的分布。從此圖中可以清楚得知聲輻射力的分布和大小，在實際的應用中可以更加有效、快捷地使用聲輻射力對生物組織施壓，以探測生物體的彈性，從而判斷生物體是否發(fā)生病變。聲輻射力的研究為彈性成像技術奠定了基礎，使得彈性成像技術更加精準、高效。

另外也可以通過聲輻射力的等值面了解到聲輻射力的分布，圖6給出了在?20 dB、?10 dB處的聲輻射力等值面。由聲輻射力等值面可以識別具有相同大小的力的所在位置。通過查看相鄰等值線的間距可以大致了解力的大小的分布層次。等值面使聲壓的分布一面了然。

圖6 腹壁組織中歸一化聲輻射力等值面Fig.6 The normalized acoustic radiation force isosurface in abdominal wall tissue

最后對面陣換能器的陣元寬度、間距、陣元個數(shù)以及工作頻率等參量對聲輻射力的影響進行了計算與分析，如圖7～10所示。

圖7為不同陣元寬度時腹壁組織中沿x方向的歸一化聲輻射力。圖7(a)和圖7(b)分別對應陣元間距為0時，不同陣元寬度的軸向和橫向聲輻射力。為了避免換能器寬度帶來的影響，陣元寬度為0.15 mm和0.30 mm時的陣元個數(shù)分別取100和50個以保證換能器的寬度相等，此時換能器的寬度都為15 mm。圖7(c)和圖7(d)分別對應陣元間距為0.15 mm時，不同陣元寬度的軸向和橫向聲輻射力。同樣的，對陣元寬度為0.15 mm和0.30 mm時的陣元個數(shù)分別取30和20個以保證換能器的寬度相等，對應的換能器的寬度都為9 mm。從圖7可以看出，換能器寬度不變的情況下，陣元間距為0時，陣元寬度對聲輻射力幾乎沒有影響。陣元間距為0.15 mm時，隨著陣元寬度的增大，聲輻射力變大。這是因為，在陣元間距為0時，不管怎么改變陣元的寬度，面陣換能器實際起作用的面的范圍沒有改變，其激勵下產(chǎn)生的聲場沒有改變。當陣元之間有間距時，隨著陣元寬度的增大，換能器實際起作用的面的范圍變大，其激勵下產(chǎn)生的聲場強度變大，從而使得組織所受的聲輻射力增大。

圖7 不同陣元寬度時腹壁組織中沿x方向的歸一化聲輻射力Fig.7 The normalized acoustic radiation force on abdominal wall tissue with different element widths along x direction

圖8為不同陣元間距時腹壁組織中沿x方向的歸一化聲輻射力。圖8(a)和圖8(b)分別對應陣元寬度為0.15 mm時，不同陣元間距的軸向和橫向聲輻射力。陣元間距為0和0.15 mm時的陣元個數(shù)分別取100和50個以保證換能器的寬度相等，此時換能器的寬度都為15 mm。從圖8可以看出，隨著陣元間距的變大，軸向聲輻射力變小，橫向聲輻射力變大。

圖8 不同陣元間距時腹壁組織中沿x方向的歸一化聲輻射力Fig.8 The normalized acoustic radiation force on abdominal tissue with different element spacing along x direction

圖9為不同陣元個數(shù)時腹壁組織中沿x方向的歸一化聲輻射力。圖9(a)和圖9(b)分別對應陣元間距為0、陣元寬度為0.15 mm時，陣元個數(shù)取100和50個時對應的軸向和橫向聲輻射力。從圖9可以看出，陣元個數(shù)的改變使得聲輻射力出現(xiàn)峰值的位置發(fā)生了偏移。這是由于在陣元間距和寬度不變的情況下，隨著陣元個數(shù)的改變，換能器寬度發(fā)生了改變，使得聲場分布發(fā)生改變，從而導致聲輻射力分布發(fā)生改變。

圖9 不同陣元個數(shù)時腹壁組織中沿x方向的歸一化聲輻射力Fig.9 The normalized acoustic radiation force on abdominal wall tissue with different number of elements along x direction

圖10為不同中心頻率時腹壁組織中沿x方向的歸一化聲輻射力。圖10(a)和圖10(b)分別對應陣元間距為0、陣元寬度為0.15 mm時，陣元個數(shù)為100、頻率為2 MHz和3 MHz時對應的軸向和橫向聲輻射力。從圖10可以看出，換能器工作頻率不同，聲輻射力的大小不同。這是因為不同工作頻率的換能器激發(fā)的聲場強度不同。

圖10 不同中心頻率時腹壁組織中沿x方向的歸一化聲輻射力Fig.10 Thenormalized acoustic radiation forceon abdominal tissuewith different center frequencies along x direction

從以上分析可知，在聲輻射力彈性成像中，可以根據(jù)需求對陣換能器的陣元寬度、間距、陣元個數(shù)以及工作頻率等參量進行調(diào)節(jié)，從而使該彈性成像技術更加精確、高效。

4 結論

利用聲輻射力對粒子和生物組織的無創(chuàng)、非接觸、多功能性等的操控是一項具有實際應用前景的技術，尤其是在彈性成像領域。分析操控需求，預測相應的聲輻射力大小和分布使其應用技術更高效、精確。數(shù)值模擬是計算聲輻射力最直接和直觀的方法。而k-Wave可以用于復雜和真實組織介質(zhì)中的時域聲學和超聲模擬，是模擬聲輻射力的一種精確而有效的工具。它結合了有限差分法的簡單性和在任何非均勻介質(zhì)中模擬的靈活性，具有快速、易用的特點。本文基于腹壁組織圖像，利用k-Wave對超聲波在腹壁組織區(qū)域傳播時的聲場進行數(shù)值模擬，得到其聲場和聲輻射力的分布。對面陣換能器的陣元寬度、間距、陣元個數(shù)以及工作頻率等參量對聲輻射力的影響進行了計算與分析，在聲輻射力彈性成像中可以按實際需要進行選擇和調(diào)整這些參量，使該彈性成像技術更加精確、高效。本文研究為聲輻射力在彈性成像技術中的應用奠定了基礎，為實現(xiàn)精密制造和精準醫(yī)療提供核心技術支持。