周　浩,鄭音飛

(浙江大學 生物醫(yī)學工程教育部重點實驗室,浙江 杭州 310027)

?

非均勻組織醫(yī)學超聲非線性傳播仿真

周浩,鄭音飛

(浙江大學 生物醫(yī)學工程教育部重點實驗室,浙江 杭州 310027)

摘要:為了實現(xiàn)仿真醫(yī)學超聲波在非均勻組織中的傳播過程,建立超聲非線性傳播計算模型.由軟組織中一階非線性波動方程推導得出“聲壓-質(zhì)點振動速度”耦合超聲非線性波動方程以降低求解復雜度.采用k空間方法對非線性波動方程組求解,在保證數(shù)值計算精度的同時降低計算的內(nèi)存占用量和計算時間.通過與一維問題的解析解和二維問題的時域有限差分(FDTD)求解結(jié)果對比,驗證所述模型的精度.在空間采樣率為聲波波長的1/9、Courant-Friedrichs-Lewy(CFL)數(shù)為0.3的情況下,所述模型的平方誤差為0.012 5%,而時域有限差分方法(FDTD)的平方誤差為42.5%.利用體腹壁組織數(shù)字模型,進行醫(yī)學超聲諧波成像仿真,驗證諧波成像較基波成像能夠提高深部組織區(qū)域的圖像質(zhì)量.

關(guān)鍵詞：非線性聲學；波動方程；超聲成像；k空間方法；諧波成像

臨床醫(yī)學中應用的醫(yī)學超聲通常處于較高頻帶范圍(1～10 MHz),這導致超聲波在生物組織的傳播過程中伴隨有顯著的非線性現(xiàn)象.醫(yī)學超聲的非線性效應已被應用于醫(yī)學超聲成像、高強度聚焦超聲治療和低強度超聲理療等方面[1-4].生物組織中超聲非線性傳播仿真可應用于醫(yī)學超聲換能器設(shè)計、新型成像算法開發(fā)、聲束畸變校正研究和醫(yī)師培訓等領(lǐng)域[5-9].相較于傳統(tǒng)線性聲場的積分求解方法[10],生物醫(yī)學超聲非線性傳播仿真的關(guān)鍵在于求解超聲傳播的非線性微分方程.利用算子分裂求解非線性Khokhlov-Zabolotskaya-Kuznetsov(KZK)方程是超聲非線性傳播仿真的常用方法[11-13],但KZK方程為拋物線型方程,僅能在較小的發(fā)射孔徑內(nèi)保證仿真精度,難以獲得后向散射回波[14].利用時域有限差分方法求解非線性全波方程(full-wave equation),在仿真超聲前向傳播的同時可以獲得后向散射回波,但由于超聲波速與頻率的限制,需要極高的時間-空間采樣率以保證仿真精度[14-16].與時域有限差分(finite-difference time-domain, FDTD)方法不同,k空間方法采用傅里葉變換計算波動方程中的空間導數(shù)項,利用修正的差分算子計算時間導數(shù)項,從而可在較低的計算開銷下保證仿真精度[17-18]；然而,傳統(tǒng)的k空間方法難以獲得超聲非線性傳播時介質(zhì)內(nèi)質(zhì)點振動速度分布情況,也難以直接應用完全匹配層(perfect matched layer, PML)方法抑制聲波在仿真計算區(qū)域邊界處的周期性入射.

本文在推導出一階耦合超聲非線性波動方程組的基礎(chǔ)上,采用k空間方法模擬超聲在組織中的非線性傳播,同步獲得生物組織中聲壓與質(zhì)點速度的分布特性.此外,采用聲壓-速度耦合的波動方程組,可以直接將PML應用于k空間方法.通過超聲非線性傳播的解析解和FDTD解驗證所提出方法的準確性.最后將該方法用于醫(yī)學超聲諧波成像仿真.

1理論

生物軟組織中一階非線性波動方程[19]為

(1)

式中：ρ′為超聲傳播導致的介質(zhì)密度變化量,u為質(zhì)點速度,p為聲壓,B/A為非線性參量,c0和ρ0分別為介質(zhì)的平衡聲速和密度,μ與μB分別為切變和膨脹黏滯系數(shù),κ為熱傳導系數(shù),cv與cp分別為定容和定壓比熱，t為時間.式(1)的數(shù)值求解較為復雜,為降低計算復雜度,式(1)可進一步合并：

(2)

(3)

式中:β=1+Β/(2Α)為介質(zhì)非線性參數(shù),

δ1=κ(1/cv-1/cp)/ρ0,δ2=(μB+4μ/3)/ρ0,

兩者分別為由熱傳導和組織黏滯導致的超聲衰減參數(shù)；δ1+δ2=δ,δ為聲衰減參數(shù),由于在生物組織中熱傳遞導致的聲衰減占總衰減比例甚少[19],在本研究中經(jīng)驗性地認為δ2=0.01δ.為了對式(3)以k空間方法求解,需要對等號右邊的非線性項進行改寫.式(3)中的非線項為二階小量,而聲壓與質(zhì)點振動速度線性關(guān)系在一階精度下成立：

(4)

可將式(4)代入式(3)中的非線性項而不影響等式推導的精度[13]，則式(3)可寫為

(5)

利用時間交錯網(wǎng)格,設(shè)Δt為時間步進長度,采用時域顯格式在時域采樣點nΔt (n=1, 2, …)上對聲壓p進行離散、在時域采樣點(n+1/2)Δt上對質(zhì)點運動速度u進行離散,則離散后的一階非線性波動方程組(2)與(5)可寫為

ξ=x, y,z.

(6)

(7)

式中：n為時間節(jié)點,ξ為質(zhì)點速度矢量方向標志.為了抑制時間導數(shù)離散引起的數(shù)值色散,利用k空間方法計算式(6)與(7)中的空間導數(shù)：

(8)

(9)

由于離散傅里葉變換具有周期性,從計算邊界一側(cè)入射的聲波會從對側(cè)邊界出射,本研究中直接利用一階完全匹配層方法使聲波在計算邊界處發(fā)生衰減[18,20].

2仿真實驗

2.1均勻介質(zhì)中的非線性傳播仿真

為了驗證所述方法對聲波非線性傳播模擬的精度,對單頻平面波在衰減介質(zhì)中的傳播進行仿真.在單頻連續(xù)平面聲波非線性傳播的情形下,式(1)轉(zhuǎn)化為Burgers方程,其解析解由文獻[21]給出.在本研究中,設(shè)定源點聲壓為

p=p0sin(2πf0t).

圖1　單頻平面波在衰減介質(zhì)中非線性傳播仿真結(jié)果Fig.1　Plane wave nonlinear propagation simulation in homogeneous lossy medium

2.2非均勻介質(zhì)中的非線性傳播仿真

圖2　非均勻介質(zhì)超聲傳播仿真Fig.2　Acoustic wave nonlinear propagation simulation in a heterogeneous medium

由于人體組織聲學參數(shù)的非均勻性,醫(yī)學超聲波在人體組織中傳播時將發(fā)生折射、散射等現(xiàn)象.為了驗證所述方法對非均勻介質(zhì)聲波非線性傳播的準確性,令聲波在如圖2(a)所示的二維非均勻介質(zhì)中傳播.此非均勻介質(zhì)由圓形異質(zhì)物和均勻背景組成,均勻背景為50mm×50mm矩形,圓形異質(zhì)物半徑為4mm,位于均勻背景中央,距圓形異質(zhì)物中央20mm處設(shè)置聲波接收點.為了與人體各類型組織間聲學參數(shù)相似,均勻背景介質(zhì)聲速設(shè)為1 500m/s,密度為1 000kg/m3,非線性參數(shù)為3.5,并設(shè)其為無衰減介質(zhì)；圓形異質(zhì)物的聲速設(shè)為1 575m/s,密度為1 050kg/m3,非線性參數(shù)為4.0,超聲衰減參數(shù)為1×10-4m2/s.令超聲發(fā)射孔徑位于均勻背景介質(zhì)的底端,以平面波方式向介質(zhì)內(nèi)發(fā)射聲波.聲波發(fā)射孔徑為20mm,發(fā)射波形為加Hanning窗的2周期正弦波,正弦波頻率為1MHz,發(fā)射波形最高幅度為0.25MPa.在使用k空間方法時,設(shè)定空間采樣間隔為非均勻介質(zhì)中聲波最小波長的1/9,CFL=0.3.采用同樣的仿真參數(shù),利用FDTD方法進行仿真以對比結(jié)果,分析在較低空間分辨率下兩種方法的精度.同時,以空間采樣間隔為非均勻介質(zhì)中聲波最小波長的1/27,CFL=0.1為仿真參數(shù),進行高時間-空間分辨率下FDTD仿真以對非均勻介質(zhì)中聲波傳播仿真提供基準.圖2(b)為3種仿真方法下,聲波接收點接收到的聲壓波形,可見k空間方法所得結(jié)果與基準結(jié)果高度符合,體現(xiàn)出k空間方法在低空間分辨率情況下的優(yōu)勢.圖2(c)為聲波接收點接收波形的頻譜A(f).k空間方法所得波形頻譜與基準波形在基頻與一次諧波附近吻合較好；由于k空間方法采用的空間采樣率較低(9倍基頻波波長),對4MHz以上的頻率成分不能有效采樣.由圖2(c)還可見,在低空間分辨率下FDTD方法幾乎不能有效包含聲波的非線性頻率成分.

2.3超聲諧波成像仿真

在驗證了所述方法對超聲波在非均勻介質(zhì)中傳播仿真準確性的基礎(chǔ)上,對醫(yī)學超聲線陣探頭所發(fā)聚焦聲波在非均勻人體組織中的傳播進行模擬,并采集后向散射回波,進行超聲諧波成像仿真.

圖3　超聲諧波成像仿真Fig. 3　Simulation of tissue harmonic ultrasound imaging

仿真中設(shè)定空間采樣間隔為33.87μm,CFL=0.3,模擬的換能器型號為UltrasonixL9-4/34高頻線陣換能器,陣元數(shù)N=128,陣元中心間距為0.304 8mm,換能器中心頻率fc=5.0MHz.模型采用文獻[22]提供的數(shù)字化二維人體腹壁組織圖為仿真介質(zhì).該腹壁組織結(jié)構(gòu)如圖3(a)所示,包含水(W)、結(jié)締組織(C)、肌肉(M)及脂肪(F).各類介質(zhì)所對應的聲學參數(shù)已由實驗測得(見表1)[17].在此基礎(chǔ)上,對聲速和密度分別加以幅度為相應物理量平均值3%的隨機波動,以在超聲傳播過程中產(chǎn)生后向散射回波,并在22mm深度處設(shè)置充滿水的均勻圓形區(qū)域作為囊性病變.相對k空間方法仿真的空間采樣率要求,原始組織切片圖像分辨率相對較低(300dpi),仿真前須對圖像進行最鄰近點插值以提高空間采樣率.仿真時,設(shè)定區(qū)域的上邊界為換能器發(fā)射表面,通過水與腹壁組織進行聲學耦合；換能器發(fā)射信號賦值于y軸方向質(zhì)點速度采樣點,以仿真換能器縱向振動模式；換能器振動信號為加Hanning窗的2周期4MHz正弦信號,換能器振動的最高速度為0.5m/s(可以在焦點處等效產(chǎn)生0.5MPa的空間峰值聲壓).

表1　腹壁組織聲學參數(shù)

設(shè)定換能器發(fā)射焦點距離F=22mm,發(fā)射陣元數(shù)為32(發(fā)射孔徑數(shù)F/#=2.3),發(fā)射線數(shù)為32.在仿真聚焦聲波傳播過程中,同步記錄計算區(qū)域上邊界的聲壓值作為陣列換能器接收到的回波信號.對回波信號進行波束合成、解調(diào)、抽取、對數(shù)壓縮、雙線性插值,獲得基頻B型超聲圖像[23],如圖3(b)所示.利用Butterworth帶通濾波器濾取波束合成信號的二次諧波成分,可獲得諧波圖像[1],如圖3(c)所示.對比圖3(b)與(c)可見,諧波成像所獲得的囊性暗區(qū)邊界比基頻成像所獲囊性暗區(qū)邊界清晰,驗證了諧波成像的優(yōu)勢.由圖3(c)可見,諧波成像模式不能對5mm以內(nèi)的淺層組織良好顯像,這是由于超聲諧波須在組織中傳輸一定距離后才能生成,較淺深度范圍內(nèi)回波中的諧波能量仍較弱.

3討論

本文所述方法使用基于傅里葉變換的k空間算子計算超聲非線性方程組中的空間導數(shù)(式(8)與(9)),計算精度要高于FDTD中的空間、時間差分運算,因而允許仿真時采用較FDTD方法更低的時間-空間采樣率.本文方法還可精確獲取超聲傳播時介質(zhì)內(nèi)質(zhì)點的振動速度,從而可用于超聲傳播時組織內(nèi)聲強矢量分布和組織升溫的仿真計算[24].

(10)

(11)

除與離散傅里葉變換相關(guān)運算過程外,式(6)與(7)所描述的迭代步進計算可以在空間各采樣點上并行完成.為了展示并行運算對仿真的加速效果,改變2.2節(jié)中所述問題的采樣點數(shù)目,在計算機操作系統(tǒng)中將CPU分別設(shè)置為單核與雙核運行條件下,測量k空間方法在時間域上單步步進所需的計算時間.測量結(jié)果如圖4所示,可見隨著網(wǎng)格數(shù)量(Gsize)的上升,雙核對仿真計算的加速效果逐漸明顯.尤其是可利用在通用圖形處理器(generalpurposegraphicprocessingunit,GPGPU)中運行的并行快速傅里葉變換程序[25],將本文所述方法移植到GPGPU中進一步提高運算速度.

圖4　不同矩陣規(guī)模下單、雙線程程序單步計算時間對比Fig.4　Comparison of compute time between single and double-thread programs under different computation grid sizes

4結(jié)語

本文在推導出一階超聲非線性波動方程的基礎(chǔ)上,提出了利用k空間方法仿真超聲在非均勻組織中傳播的模型,同步獲得介質(zhì)內(nèi)聲壓與質(zhì)點振動速度分布情況.與FDTD方法相比,k空間方法能夠顯著降低仿真時的內(nèi)存占用量和計算用時.通過仿真聲波在一維衰減組織與二維非均勻組織中的傳播,驗證了k空間方法的準確性,并在此基礎(chǔ)上進行了超聲諧波成像的仿真.超聲諧波成像的仿真結(jié)果表明,諧波成像較基波成像能夠明顯提高深部組織區(qū)域的圖像質(zhì)量.

參考文獻(References)：

[1]SZABOTL. Diagnostic ultrasound imaging [M]. Burlington: Elsevier, 2004:381-428.[2] 周浩, 王友釗, 鄭音飛. 醫(yī)學超聲編碼激勵技術(shù)研究進展 [J]. 浙江大學學報:工學版, 2011, 45(2):387-391.

ZHOU Hao, WANG You-zhao, ZHENG Yin-fei. The advance of medical ultrasonic coded excitation [J]. Journal of Zhejiang University:Engineering Science, 2011, 45(2): 387-391.

[3] JI X, BAI J F, SHEN G F, et al. High-intensity focused ultrasound with large scale spherical phased array for the ablation of deep tumors [J]. Journal of Zhejiang University SCIENCE B, 2009, 10(9):639-647.

[4] YING Z M, LIN T, YAN S G. Low-intensity pulsed ultrasound therapy: a potential strategy to stimulate tendon-bone junction healing [J]. Journal of Zhejiang University SCIENCE B, 2012, 13(12): 955-963.

[5] DHANALIWALA A H, HOSSACK J A, MAULDIN F W. Assessing and improving acoustic radiation force image quality using a 1.5-D transducer design [J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2012, 59(7): 1602-1608.

[6] VARRAY, F, BASSET O, TORTOLI P, et al. Extensions of nonlinear B/A parameter imaging methods for echo mode [J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2011, 58(6): 1232-1244.

[7] 鄭音飛, 邱飛岳, 李鵬, 等. 基于發(fā)射束迭代算法的超聲波色差抑制方法研究 [J]. 傳感技術(shù)學報, 2009, 22(2): 240-243.

ZHENG Yin-fei, QIU Fei-yue, LI Peng, et al. Study of reduction aberration of ultrasound image based on iteration of transmit [J]. Chinese Journal of Sensors and Actuators, 2009, 22(2): 240-243.

[8] GJERALD S U, BREKKEN R, HERGUM T, et al. Real-time ultrasound simulation using the GPU [J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2012, 59(5): 885-892.

[9] H?ILUND-KAUPANG H, M?S?Y S E. Transmit beamforming for optimal second harmonic generation [J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2011, 58(8): 1559-1569.

[10] 吳施偉, 吳海騰, 金浩然, 等. 聚焦探頭水浸檢測下的頻域合成孔徑聚焦技術(shù) [J]. 浙江大學學報:工學版, 2015, 49(1):110-115.

WU Shi-wei, WU Hai-teng, JIN Hao-ran, et al. Frequency-domain synthetic aperture focusing technique for immersion testing using focused transducer [J]. Journal of Zhejiang University:Engineering Science, 2015, 49(1):110-115.

[11] VARSLOT T, TARALDSEN G. Computer simulation of forward wave propagation in soft tissue [J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2005, 52(9): 1473-1482.

[12] JING Y, CLEVELAND R O. Modeling the propagation of nonlinear three-dimensional acoustic beams in inhomogeneous media [J]. Journal of the Acoustical Society of America, 2007, 122(3): 1352-1364.

[13] 劉偉, 楊軍, 田靜. 有限振幅波的三維時域建模[J]. 聲學學報, 2010,36(4): 349-357.

LIU Wei, YANG Jun, TIAN Jing. Time-domain modeling of finite-amplitude sound beams in three-dimensional Cartesian coordinate system [J]. Acta Acoustica, 2010, 36(4): 349-357.

[14] 解卓麗, 周浩, 鄭音飛. 非均勻組織醫(yī)學超聲發(fā)射聲場仿真. 聲學學報, 2013, 38(6):657-662.

XIE Zhuo-li, ZHOU Hao, ZHENG Yin-fei. Simulation of the acoustic field emitted from medical linear transducer in a heterogeneous tissue [J]. Acta Acoustica, 2013, 38(6):657-662.

[15] HALLAJ O, CLEVELAND R O. FDTD simulation of finite-amplitude pressure and temperature fields for biomedical ultrasound [J]. Journal of the Acoustical Society of America, 1999, 105(5):L7-L12.

[16] PINTON G F, DAHL J, ROSENZWEIG S, et al. A heterogeneous nonlinear attenuating full-wave model of ultrasound [J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2009, 56(3): 474-488.

[17] JING Y, WANG T, CLEMENT G T. A k-space method for moderately nonlinear wave propagation [J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2012, 59(8): 1664-1673.

[18] TABEI M, MAST T D, WAAG R C. A k-space method for coupled first-order acoustic propagation equations [J]. Journal of the Acoustical Society of America, 2002, 111(1): 53-63.

[19] HAMILTON M F, MORFEY C L. Model equations [M]∥ Nonlinear Acoustics. Melville: Acoustical Society of America, 2008: 41-62.

[20] YUAN X, BORUP D, WISKIN J W, et al. Formulation and validation of Berenger’s PML absorbing boundary for the FDTD simulation of acoustic scattering [J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 1997, 44(4): 816-822.

[21] BLACKSTOCK D T, HAMILTON M F, PIERCE A D. Progressive waves in lossless and lossy fluids [M]∥ Nonlinear Acoustics. Melville: Acoustical Society of America, 2008: 65-147.

[22] MAST T D, HINKELMAN L M, ORR M J, et al.. Simulation of ultrasonic pulse propagation through the abdominal wall [J]. Journal of the Acoustical Society of America, 1997, 102(2): 1177-1190.

[23] ZHOU H, ZHENG Y F. An efficient quadrature demodulator for medical ultrasound imaging [J]. Frontiers of Information Technology and Electronic Engineering, 2015, 16(4): 301-310.

[24] COX B. Acoustics for ultrasound imaging [EB/OL]. [2015-3-10]. http:∥www.ucl.ac.uk/medphys/staff/people/bcox/bens_lecture_notes.

[25] NVIDIA Corporation. cuFFT [CP/OL]. [2015-3-10]. https://developer.nvidia.com/cufft.

DOI:10.3785/j.issn.1008-973X.2016.03.023

收稿日期：2015-01-04.

基金項目：中央高?；究蒲袠I(yè)務費專項資金資助項目(2015FZA5019,2016FZA5015).

作者簡介：周浩(1984-)，男，博士生，從事醫(yī)學超聲成像技術(shù)研究. ORCID：0000-0001-6894-1139. E-mail:bmezhou@zju.edu.cn 通信聯(lián)系人：鄭音飛，男，副教授. ORCID：0000-0001-6837-2634. E-mail:zyfnjupt@126.com

中圖分類號：TN 98

文獻標志碼：A

文章編號：1008-973X(2016)03-06-0574

Simulation of nonlinear ultrasound propagation in heterogeneous tissue

ZHOU Hao, ZHENG Yin-fei

(KeyLaboratoryforBiomedicalEngineeringofMinistryofEducation,ZhejiangUniversity,Hangzhou310027,China)

Abstract:A numerical model was proposed for the simulation of the nonlinear ultrasound propagation in heterogeneous tissue. First, the coupled nonlinear wave equations for pressure and velocity were obtained based on 1-st order nonlinear wave equations in soft tissue to reduce the complexity of numerical computation. Then, k-space method was used to solve the derived nonlinear wave equations to reduce the memory usage and the computation time of the simulation, while preserving the computation accuracy. Compared with the analytic solution of a 1-dimensinal problem and the finite-different time-domain (FDTD) results of a 2-dimensinal problem, and the accuracy of the proposed model was validated. With grid size of 1/9 of the wavelength and Courant-Friedrichs-Lewy (CFL) of 0.3, the square errors of the proposed model and the FDTD method are 0.0125% and 42.5%, respectively. Medical harmonic ultrasound imaging was simulated using the proposed method based on a digital human abdominal map. The results show that image quality can be improved in the deeper tissue by using the harmonic signal．

Key words:nonlinear acoustics; wave equation; medical ultrasound imaging; k-space method; harmonic medical ultrasound imaging