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

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

(1 近代聲學(xué)教育部重點(diǎn)實(shí)驗(yàn)室 南京大學(xué)聲學(xué)研究所 人工微結(jié)構(gòu)科學(xué)與技術(shù)協(xié)同創(chuàng)新中心 南京 210093)

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

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

0 引言

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

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

1 理論

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

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

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

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

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

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

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

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

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

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

2 仿真模型的建立

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

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

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

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

3 數(shù)值仿真

3.1 有效性驗(yàn)證

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

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

3.2 腹壁中的聲輻射力

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖10為不同中心頻率時(shí)腹壁組織中沿x方向的歸一化聲輻射力。圖10(a)和圖10(b)分別對(duì)應(yīng)陣元間距為0、陣元寬度為0.15 mm時(shí)，陣元個(gè)數(shù)為100、頻率為2 MHz和3 MHz時(shí)對(duì)應(yīng)的軸向和橫向聲輻射力。從圖10可以看出，換能器工作頻率不同，聲輻射力的大小不同。這是因?yàn)椴煌ぷ黝l率的換能器激發(fā)的聲場(chǎng)強(qiáng)度不同。

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

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

4 結(jié)論

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