基于時間反轉(zhuǎn)法的相控?fù)Q能器聲場的仿真研究?

韓珍珍　丁鑫　羅明艷　菅喜岐

（天津醫(yī)科大學(xué)生物醫(yī)學(xué)工程與技術(shù)學(xué)院天津300070）

基于時間反轉(zhuǎn)法的相控?fù)Q能器聲場的仿真研究?

韓珍珍丁鑫羅明艷菅喜岐?

（天津醫(yī)科大學(xué)生物醫(yī)學(xué)工程與技術(shù)學(xué)院天津300070）

相控?fù)Q能器具有焦距可調(diào)的優(yōu)勢。本文以82陣元相控?fù)Q能器建立的3D數(shù)值仿真模型為例，基于時間反轉(zhuǎn)法提取陣元的激勵信號，利用時域有限差分（FDTD）法對Westervelt聲波非線性傳播方程進(jìn)行聲場數(shù)值仿真，研究不同陣元分布、焦點偏離聲軸的距離、設(shè)定焦距大小對形成聲場的影響，可調(diào)控范圍及其消除旁瓣方法。研究結(jié)果表明，隨機分布相控陣可明顯降低聲場中的旁瓣；隨著偏離聲軸距離的增加，主瓣聲壓幅值逐漸減小，旁瓣與主瓣的最大聲強比值r逐漸增大，且沿聲軸的可調(diào)控范圍逐漸減??；隨聲軸方向上設(shè)定焦距的增加，主瓣聲壓幅值先增大后減小，r值先減小后增大；基于時間反轉(zhuǎn)法的高聲壓旁瓣消除法可在一定程度上擴大相控陣聲場的可調(diào)控范圍。

高強度聚焦超聲，相控?fù)Q能器，陣元分布，時間反轉(zhuǎn)法，高聲壓旁瓣消除法

1　引言

高強度聚焦超聲腫瘤治療具有無創(chuàng)、可重復(fù)治療等優(yōu)勢，超聲換能器是HIFU腫瘤治療系統(tǒng)中最關(guān)鍵的部件之一。目前單陣元換能器已應(yīng)用于乳腺癌、子宮肌瘤、前列腺癌等實體軟組織腫瘤的臨床治療［1］，但其焦距固定，腫瘤治療時需通過機械調(diào)整換能器與病患部位的相對位置進(jìn)行治療定位。而相控陣換能器可通過相位調(diào)控其焦距大小而實現(xiàn)動態(tài)聚焦，該研究受到眾多研究者的關(guān)注。

1996年Goss等［2］提出球面陣元隨機分布可改善相控陣聲場分布的方法。2005年Lu等［3］采用遺傳算法對256矩形陣元相控陣中每個陣元的幅值和相位進(jìn)行調(diào)控以降低聲場中旁瓣。2009年Hand等［4］通過仿真和離體實驗對254陣元隨機相控陣單焦點和多焦點的動態(tài)調(diào)控范圍進(jìn)行了測定。2012年Ji等［5］設(shè)計了用于治療深部腫瘤的90陣元大開口半隨機分布的球面相控陣。

本文根據(jù)周期、同心環(huán)和隨機分布設(shè)計的82陣元相控陣換能器建立3D數(shù)值仿真模型，并基于時間反轉(zhuǎn)法提取相位控制信號，利用Westervelt聲波非線性傳播方程進(jìn)行時域有限差分?jǐn)?shù)值仿真研究，以旁瓣與主瓣的最大聲強比值r≤0.25為判斷標(biāo)準(zhǔn)［6］討論不同聚焦位置時的聲場變化規(guī)律、焦點可調(diào)控范圍。并提出擴大焦域可調(diào)控范圍的高聲壓旁瓣消除法。

2　方法與模型

2.1聲波傳播方程式

Westervelt聲波傳播方程式［7-8］為

（1）式中：?為拉普拉斯算子；p為聲壓；c0和ρ0分別為介質(zhì)的聲速和密度；β=1+B/（2A）為聲波非線性系數(shù)，B/A為流體介質(zhì)的非線性系數(shù)；為聲波擴散系數(shù)；α為吸收系數(shù)；ω=2πf為角頻率，f為頻率。

采用時域有限差分法［9-10］對（1）式進(jìn)行中心差分，差分方程式為

i、j、k分別為直角坐標(biāo)系下x、y、z三個坐標(biāo)軸方向的坐標(biāo)，dx、dy、dz分別表示x、y、z三個坐標(biāo)軸方向的空間步長，dt為時間步長，n為計算時刻。

2.2陣元相位獲取法

基于時間反轉(zhuǎn)法［11-12］在目標(biāo)焦點S處設(shè)置正弦函數(shù)點聲源S0（t），數(shù)值仿真得到傳到相控陣編號為m陣元中心點的聲壓信號S0m（t），將該信號按時間序列進(jìn)行反轉(zhuǎn)后，得到對應(yīng)陣元m的信號S0m（T-t）。利用最小二乘函數(shù)擬合法計算一段時間內(nèi)S0m（T-t）的相對初始相位延遲Δtm，然后以同一輸入聲強對正弦信號幅值進(jìn)行調(diào)制，陣元m的激勵信號為

如果在（7）式激勵信號觸發(fā)相控陣換能器陣元形成的聲場中出現(xiàn)聲壓幅值較高的旁瓣，確定其位置R，采用同樣的方法計算聚焦在R處相控?fù)Q能器陣元m的激勵信號。對分別聚焦在S和R處的激勵信號線性疊加，調(diào)控聚焦在R處的輸入聲強使在R處聚焦時穩(wěn)態(tài)聲壓的最大幅值與聚焦在S處時R位置的穩(wěn)態(tài)聲壓幅值相等，且使聚焦在R處的相位相對聚焦在S處的相位180°延遲，從而達(dá)到降低R處高聲壓旁瓣幅值［13］的目的。

2.3數(shù)值仿真模型和參數(shù)

圖1為開口直徑100 mm、曲率半徑80 mm、曲率半徑與開口直徑之比F值為0.8的82陣元相控陣換能器和水體組成的數(shù)值仿真模型。其中數(shù)值仿真區(qū)域為邊長為100 mm的立方體，聲軸為x軸。圖2（a）、2（b）、2（c）分別為陣元按照周期、同心環(huán)和隨機分布構(gòu)成的82陣元相控?fù)Q能器，其中陣元直徑均為8 mm、陣元填充率為52.48%，周期分布相鄰陣元中心間距均為9.5 mm；同心環(huán)分布中陣元中心點所在環(huán)半徑由內(nèi)到外依次為9 mm、18 mm、27 mm、36 mm、45 mm，對應(yīng)的陣元數(shù)目分別為6、11、16、21、27；隨機分布通過Matlab隨機數(shù)生成器獲得，四個象限內(nèi)陣元數(shù)目依次為21、21、20、20，在相鄰陣元中心間距不小于8.7 mm的條件下，依次選定各象限中其余陣元的中心位置。

水體的仿真參數(shù)為ρ0=998 kg/m3，c0= 1486 m/s，β=3.5，α=0.02 Np/（m·MHz）［14］。換能器激勵函數(shù)為頻率為0.7 MHz的連續(xù)正弦波，邊界采用Mur一階邊界吸收條件?？臻g步長dx=dy=dz=0.25 mm，時間步長dt=10 ns。

圖1　數(shù)值仿真模型圖（單位：mm）Fig.1 Numerical simulation model（Unit：mm）

圖2　陣列結(jié)構(gòu)圖Fig.2 Array structure diagram

3　結(jié)果

3.1不同陣元分布形成的聲場

當(dāng)輸入聲強為1 W/cm2，設(shè)定焦點為幾何焦點（80，0，0）mm時，周期分布、同心環(huán)分布和隨機分布相控陣在x軸上形成的聲壓變化曲線如圖3所示，其中黑色虛線為周期分布、灰色實線為同心環(huán)分布、黑色實線為隨機分布。由圖3可知，三種分布陣列形成的最大聲壓（主瓣）幅值分別為4.800 MPa、4.716 MPa、4.777 MPa，實際焦距分別為78.25 mm、77.75 mm、77.75 mm，均小于設(shè)定的80 mm幾何焦距，周期分布和同心環(huán)分布陣列在主瓣左側(cè)旁瓣多且幅值較大，隨機分布相控陣僅在靠近主瓣左側(cè)有一個較小旁瓣。圖4為三種不同分布陣列在幾何焦點聚焦時實際形成的焦平面聲場分布圖，由圖4可知周期分布陣列焦域周圍有六個明顯的旁瓣，同心環(huán)分布陣列焦域周圍有兩層旁瓣，隨機分布陣列僅在焦域附近有分布無規(guī)則的較小旁瓣。對上述三種不同分布相控陣在幾何焦點聚焦時實際形成的三維聲場進(jìn)行分析，周期、同心環(huán)、隨機分布陣列中旁瓣與主瓣的最大聲強比值分別為0.1400、0.2358、0.1247，隨機分布陣列的比值最低。

圖3　不同陣元分布相控陣在x軸上形成的聲壓變化曲線圖Fig.3The acoustic pressure curve along the xaxis of phased arrays with different element distribution

圖4　在幾何焦點聚焦時的焦平面聲場分布圖Fig.4 The acoustic field distribution of the focal plane when focusing at geometric focus

3.2偏離聲軸聚焦的聲場分布

當(dāng)輸入聲強為1 W/cm2，設(shè)定焦點位置為x軸80 mm、z軸0 mm、沿y軸正方向以1 mm為步長偏離x軸聚焦時，形成焦點位置處最大聲壓的變化曲線如圖5（a）所示，其中灰色實線為周期分布、黑色虛線為同心環(huán)分布、黑色實線為隨機分布，由圖5（a）可知隨著偏離x軸距離的增加三種分布相控陣最大聲壓值均逐漸減小，且變化趨勢基本一致。圖5（b）為與圖5（a）相應(yīng)的旁瓣與主瓣最大聲強比值r的變化曲線，由圖5（b）可知隨著偏離x軸距離的增加該比值均逐漸增大，同心環(huán)分布和周期分布陣列的r值均明顯高于隨機分布陣列。由圖3～5可知，與周期和同心環(huán)分布陣列相比，隨機分布陣列可以明顯降低聲場中的旁瓣，為此選用隨機分布相控?fù)Q能器進(jìn)行聲場討論分析。

圖5　三種不同分布陣列沿y軸聚焦時的聲場變化Fig.5The acoustic variation as a function of focal spot location along the y-axis for three phased arrays with different element distribution

以圖2（c）隨機分布相控陣為例，當(dāng)設(shè)定焦點為（80，5，0）、（80，10，0）、（80，15，0）時實際形成的焦平面聲場分布如圖6所示，由圖6可知隨著偏離x軸距離的增加，分布在焦域周圍的旁瓣逐漸移至焦域下方，并且旁瓣數(shù)量逐漸增多、幅值逐漸增大。由圖5（b）可知，隨機分布相控陣在偏離x軸8 mm聚焦形成的聲場中r=0.2664，超過了0.25，其xy平面聲場分布如圖7（a）所示，對圖7（a）所示的靠近焦域左下方的旁瓣進(jìn)行高聲壓旁瓣消除后，xy平面的聲場分布如圖7（b）所示，使其旁瓣幅值明顯降低，且r=0.1534。同時高聲壓旁瓣消除后，二次聲壓旁瓣幅值由1.963 MPa下降為1.192 MPa，下降率為39.28%，而主瓣聲壓幅值由3.803 MPa變?yōu)?.403 MPa，下降率僅為10.52%。當(dāng)對隨機分布相控陣進(jìn)行高聲壓旁瓣消除后，可使其實現(xiàn)最大偏離x軸12 mm的調(diào)控聚焦。

圖6　隨機分布相控陣在不同焦點位置聚焦時歸一化后焦平面聲場分布圖Fig.6 The acoustic field distribution of the focal plane for the quasi-random phased array with different focal location

圖7　隨機分布相控陣在（80，8，0）聚焦時歸一化后xy平面聲場分布圖Fig.7The acoustic field distribution in the xy plane for the quasi-random phased array when focusing at（80，8，0）

3.3隨機分布相控陣沿聲軸聚焦

當(dāng)輸入聲強為1 W/cm2，設(shè)定目標(biāo)焦點在z=0平面，y方向上分別沿x軸、偏離x軸5 mm、10 mm、12 mm，沿x軸方向的焦距范圍均在55～90 mm調(diào)控的條件下，形成聲場中最大聲壓值隨設(shè)定焦點位置的變化曲線如圖8（a）所示，其中灰色虛線為沿x軸、灰色實線為偏離x軸5 mm、黑色虛線為偏離x軸10 mm、黑色實線為偏離x軸12 mm，由圖8（a）可知在y軸方向偏離x軸距離相同的條件下，隨x軸方向上設(shè)定焦距的增加，形成的最大聲壓值先增大后減小，且偏離x軸的距離不同形成最大聲壓的位置也不同，分別為x=70 mm、75 mm、75 mm、80 mm。在x軸方向上調(diào)控位置相同時，隨著y軸方向上偏離x軸距離的增加，形成的最大聲壓值逐漸降低。圖8（b）為與圖8（a）相同條件下隨設(shè)定焦點位置的r值變化曲線，由圖8（b）可知隨著y方向上偏離x軸距離的增加，r值逐漸增大，隨著x軸方向上設(shè)定焦距值的增加r值先減小后增大，形成最小r值的位置有所不同，分別為x=65 mm、65 mm、75 mm、80 mm。在r≤0.25的判斷條件下，x軸上的可調(diào)控范圍為55～90 mm，y方向上偏離x軸5 mm、10 mm沿x軸方向聚焦的可調(diào)控范圍分別為55～90 mm、60～80 mm，y方向上偏離x軸12 mm僅有x=80 mm和x=90 mm兩點。隨著y方向上偏離x軸距離的增加，沿x軸方向的聲場可調(diào)控范圍逐漸減小。同時由圖8（b）可知，在偏離x軸10 mm、12 mm聚焦時，分別有x=55、60、85、90 mm和x=55、60、65、70、75、85 mm位置處r值均大于0.25，且在x=85 mm調(diào)控聚焦時均有一波峰值，分別進(jìn)行高聲壓旁瓣消除后的r值變化曲線如圖8（c）所示，由8（c）可知偏離x軸10 mm、12 mm分別對應(yīng)x軸方向可調(diào)控范圍分別不小于60～90 mm、80～90 mm。焦點在（80，12，0）位置調(diào)控聚焦時形成的xy和yz平面聲場分布如圖9所示。由圖9可知，在焦域附近沒有較大的聲壓旁瓣出現(xiàn)。

圖8　隨機分布相控陣y方向上偏離x軸不同距離時沿x軸以5 mm步長調(diào)控聚焦時的聲場變化規(guī)律Fig.8 The acoustic variation as a function of focal spot location along the x-axis and off y-axis different distance from x-axis for the quasi-random phased array

圖9　隨機分布相控陣在（80，12，0）位置聚焦時的歸一化聲場分布圖Fig.9 The acoustic field distribution when focusing at（80，12，0）for the quasi-random phased array

4　討論與結(jié)論

本文設(shè)計了周期、同心環(huán)、隨機三種不同分布的相控陣，基于時間反轉(zhuǎn)法進(jìn)行相位調(diào)控，分析比較其聲場分布特點和變化規(guī)律，并結(jié)合基于時間反轉(zhuǎn)法的高聲壓旁瓣消除法增大隨機分布相控陣聲場的可調(diào)控范圍。結(jié)論如下：

（1）與周期和同心環(huán)分布陣列相比，隨機分布陣列可明顯降低聲場中的旁瓣，改善聲場分布，這與Gavrilov等［6］的研究結(jié)果一致。

（2）隨聲軸方向上設(shè)定焦點和換能器間距離的增加，主瓣聲壓幅值先增大后減小，旁瓣與主瓣的最大聲強比值先減小后增大。同時隨著偏離聲軸距離的增加，主瓣聲壓幅值逐漸減小，旁瓣與主瓣的最大聲強比值逐漸增大，且沿聲軸的聲場可調(diào)控范圍逐漸減小，該結(jié)果與Raju等［15］的研究結(jié)果相近。

（3）基于時間反轉(zhuǎn)法的高聲壓旁瓣消除法可在一定程度上增大相控陣聲場的可調(diào)控范圍。

上述結(jié)論是在高強度聚焦聲場的仿真研究中得出的，對于聲壓幅值較小的超聲波而言，一般呈線性傳播，但對HIFU治療腫瘤而言，其焦點聲壓一般較高，因此需要考慮超聲傳播的非線性特性。

在基于時間反轉(zhuǎn)法獲取陣元的激勵信號時，聲波的不同傳遞過程中具有一定非線性差異，2010年P(guān)inton等［16］的研究結(jié)果表明對相控陣的相位差計算時其影響可忽略不計。為此本文在獲取陣元激勵信號時未考慮其非線性的影響。

本文僅以水體內(nèi)相控?fù)Q能器形成的聲壓場進(jìn)行分析討論。對于非均質(zhì)人體組織的仿真研究中，可利用2003年Aubry等［17］提出的基于CT數(shù)據(jù)獲取組織聲學(xué)參數(shù)的方法，建立結(jié)合臨床實際的數(shù)值仿真模型，研究其形成的聲場特性。基于人體CT數(shù)據(jù)的數(shù)值仿真正在進(jìn)行中。

［1］KENNEDY J E.High-intensity focused ultrasound in the treatment of solid tumours［J］.Nature Reviews Cancer，2005，5（4）：321-327.

［2］GOSS S A，F(xiàn)RIZZELL L A，KOUZMANOFF J T，et al.Sparse random ultrasound phased array for focal surgery［J］.IEEE Trans.Ultrason.Ferroelectr.Freq. Control，1996，43（6）：1111-1121.

［3］LU Mingzhu，WAN Mingxi，XU Feng，et al.Focused beam control for ultrasound surgery with spherical-section phased array：sound field calculation and genetic optimization algorithm［J］.IEEE Trans.Ultrason.Ferroelectr. Freq.Control，2005，52（8）：1270-1290.

［4］HAND J W，SHAW A，SADHOO N，et al.A random phased array device for delivery of high intensity focused ultrasound［J］.Phys.Med.Biol.，2009，54（19）：5675-5693.

［5］JI Xiang，SHEN Guofeng，BAI Jingfeng，et al.The characterization of an ultrasound spherical phased array for the ablation of deep-seated tissue［J］.Appl.Acoust.，2012，73（5）：529-534.

［6］GAVRILOV L R，HAND J W.A theoretical assessment of the relative performance of spherical phased arrays for ultrasound surgery［J］.IEEE Trans.Ultrason.Ferroelectr. Freq.Control，2000，47（1）：125-139.

［7］HAMILTON M F，BLOCKSTOCK D T.Nonlinear Acoustics［M］.Boston：Academic Press，1998.

［8］BAILEY M，KHOKHLOVA V，SAPOZHNIKOV O，et al. Physical mechanism of the therapeutic effect of ultrasound（A review）［J］.Acoustic Physics，2003，49（4）：369-388.

［9］YEE K S.Numerical solution of initial boundary value problems involving Maxwell's equations［J］.IEEE Trans. Antennas Propag，1966，14（3）：302-307.

［10］RADZEVICIUS S J，CHEN C C，ETERS J L，et al.Nearfield dipole radiation dynamics through FDTD modeling［J］.Journal of Applied Geophysics，2003，52（2）：75-91.

［11］THOMAS J L，F(xiàn)INK M A.Ultrasonic beam focusing through tissue inhomogeneities with a time reversal mirror：application to transskull therapy［J］.IEEE Trans. Ultrason.Ferroelectr.Freq.Control，1996，43（6）：1122-1129.

［12］王意喆，張千，周文征，等.HIFU經(jīng)顱治療中降低顱骨聲壓的仿真研究［J］.計算機仿真，2014，31（5）：209-213. WANG Yizhe，ZHANG Qian，ZHOU Wenzheng，et al. Simulation study on reducing acoustic pressure on skull during HIFU transcranial therapy［J］.Computer Simulation，2014，31（5）：209-213.

［13］LEDUC N，OKITA K，SUGIYAMA K，et al.Focus control in HIFU therapy assisted by time-reversal simulation with an iterative procedure for hot spot elimination［J］. Journal of Biomechanical Science and Engineering，2012，7（1）：43-56.

［14］GINTERS.Numericalsimulationofultrasoundthermotherapy combining nonlinear wave propagation with broadband soft-tissue absorption［J］.Ultrasonics，2000，37（10）：693-696.

［15］RAJU B I，HALL C S，SEIP R.Ultrasound therapy transducers with space-filling non-periodic arrays［J］.IEEE Trans.Ultrason.Ferroelectr.Freq.Control，2011，58（5）：944-954.

［16］PINTON G，AUBRY J F，F(xiàn)INK M，et al.Effects of nonlinear ultrasound propagation on high intensity brain therapy［J］.Medical Physics，2011，38（3）：1207-1216.

［17］AUBRY J F，TANTER M，PERNOT M，et al.Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans［J］. The Journal of the Acoustical Society of America，2003，113（1）：84-93.

The simulation study of acoustic pressure field about HIFU phased array transducer based on time reversal method

HAN ZhenzhenDING XinLUO MingyanJIAN Xiqi
（Department of Biomedical Engineering and Technology，Tianjin Medical University，Tianjin 300070，China）

Phased array transducer has the advantage of adjustable focal length.In this paper，the incentive signal of each element is abstracted by time reversal method and the acoustic pressure numerical simulation is undergone through finite difference time domain（FDTD）of Westervelt acoustic nonlinear transmission formula based on the 3D numerical simulation model of 82-element phased array transducer.The effect of different element distribution，the distance off the acoustic axis and the setting focal length on the forming acoustic pressure field is discussed.The adjustable range of acoustic pressure field and the method which eliminate sidelobe is also analyzed.The results show that the phased array transducer with quasi-random distribution can suppress the sidelobe significantly.As the distance off the acoustic axis increasing，the acoustic pressure amplitude of the main lobe decreases，the acoustic intensity ratio between the maximum sidelobe and the main lobe（r value）increases，and the adjustable range along the acoustic axis decreases.The acoustic pressure amplitude of the main lobe increases firstly and then decreases，while r value decreases firstly and then increases with the increase of setting focal length along the acoustic axis.To some extent，the acoustic adjustable range can be extended by the high acoustic pressure sidelobe elimination method based on time reversal method.

High-intensity focused ultrasound，Phased array transducer，Element distribution，Time reversal method，High acoustic pressure sidelobe elimination method

R318.5

A

1000-310X（2015）04-0344-07

10.11684/j.issn.1000-310X.2015.04.009

2014-10-29收稿；2015-02-16定稿

?國家自然科學(xué)基金項目（81272495）

韓珍珍（1989-），女，山東濰坊人，碩士研究生，研究方向：超聲醫(yī)學(xué)。?

E-mail：jianxiqi@tmu.edu.cn