基于時間反轉法的相控換能器聲場的仿真研究?

韓珍珍　丁鑫　羅明艷　菅喜岐

（天津醫科大學生物醫學工程與技術學院天津300070）

基于時間反轉法的相控換能器聲場的仿真研究?

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

（天津醫科大學生物醫學工程與技術學院天津300070）

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

高強度聚焦超聲，相控換能器，陣元分布，時間反轉法，高聲壓旁瓣消除法

1　引言

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

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

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

2　方法與模型

2.1聲波傳播方程式

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

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

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

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

2.2陣元相位獲取法

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

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

2.3數值仿真模型和參數

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

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

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

圖2　陣列結構圖Fig.2 Array structure diagram

3　結果

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

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

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

圖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）隨機分布相控陣為例，當設定焦點為（80，5，0）、（80，10，0）、（80，15，0）時實際形成的焦平面聲場分布如圖6所示，由圖6可知隨著偏離x軸距離的增加，分布在焦域周圍的旁瓣逐漸移至焦域下方，并且旁瓣數量逐漸增多、幅值逐漸增大。由圖5（b）可知，隨機分布相控陣在偏離x軸8 mm聚焦形成的聲場中r=0.2664，超過了0.25，其xy平面聲場分布如圖7（a）所示，對圖7（a）所示的靠近焦域左下方的旁瓣進行高聲壓旁瓣消除后，xy平面的聲場分布如圖7（b）所示，使其旁瓣幅值明顯降低，且r=0.1534。同時高聲壓旁瓣消除后，二次聲壓旁瓣幅值由1.963 MPa下降為1.192 MPa，下降率為39.28%，而主瓣聲壓幅值由3.803 MPa變為3.403 MPa，下降率僅為10.52%。當對隨機分布相控陣進行高聲壓旁瓣消除后，可使其實現最大偏離x軸12 mm的調控聚焦。

圖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隨機分布相控陣沿聲軸聚焦

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

圖8　隨機分布相控陣y方向上偏離x軸不同距離時沿x軸以5 mm步長調控聚焦時的聲場變化規律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　討論與結論

本文設計了周期、同心環、隨機三種不同分布的相控陣，基于時間反轉法進行相位調控，分析比較其聲場分布特點和變化規律，并結合基于時間反轉法的高聲壓旁瓣消除法增大隨機分布相控陣聲場的可調控范圍。結論如下：

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

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

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

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

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

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

［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，FRIZZELL 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，FINK 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］.計算機仿真，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，FINK 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定稿

?國家自然科學基金項目（81272495）

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

E-mail：jianxiqi@tmu.edu.cn