年 9 月 1 引言人体组织中的精确超声聚焦在超声成像高强度聚焦超声 (High intensity focused ultrasound, HIFU) 治疗超声神经调控中都有重要作用, 对无骨骼遮挡软组织的超声聚焦已在临床中得到了广泛应用而在如 HIFU 经颅脑肿瘤治疗深部脑

Size: px

Start display at page:

Download "年 9 月 1 引言人体组织中的精确超声聚焦在超声成像高强度聚焦超声 (High intensity focused ultrasound, HIFU) 治疗超声神经调控中都有重要作用, 对无骨骼遮挡软组织的超声聚焦已在临床中得到了广泛应用而在如 HIFU 经颅脑肿瘤治疗深部脑"

皮霄王
5 years ago
Views:

1 第 35 卷第 5 期 Vol. 35, No.5 21 年 9 月 Journal of Applied Acoustics September, 年度全国检测声学会议优秀论文专栏基于虚拟源时间反转的经颅超声精确聚焦苏畅 1 彭哲凡 1,2 林伟军 1 (1 中国科学院声学研究所北京 119) (2 中国科学院大学北京 149) 摘要针对经颅超声难以精确聚焦的问题, 本文研究了基于虚拟源的时间反转方法, 建立颅骨二维数字模型, 在聚焦目标处设置虚拟声源, 利用时域有限差分法模拟超声时间反转过程, 并考察基于虚拟源时间反转方法的聚焦效果数值仿真结果表明, 基于虚拟源的时间反转方法可以实现经颅超声的精确聚焦, 聚焦强度和精度好于传统的相控聚焦 ; 换能器中心频率和数量对聚焦效果的影响规律与相控聚焦时类似 ; 该方法可以同时向多个焦点聚焦, 并自适应调节各焦点处声压幅度关键词超声聚焦, 时间反转, 经颅超声, 数值模拟中图分类号 : O42.2 文献标识码 : A 文章编号 : 1-31X(21) DOI: 1.14/j.issn.1-31X Accurate transcranial ultrasound focusing with time reversal based on virtual sources SU Chang 1 PENG Zhefan 1,2 LIN Weijun 1 (1 Institute of Acoustics, Chinese Academy of Sciences, Beijing 119, China) (2 University of Chinese Academy of Sciences, Beijing 149, China) Abstract To overcome the difficulty of ultrasound to accurately focus in the brain, a time reversal method based on virtual sources is studied in this article. A 2-D numerical model of human skull with virtual sources at the focusing targets is established. The time reversal process is simulated with finite difference time domain method, and the effectiveness of focusing with time reversal method is studied. The numerical simulation results show that time reversal based on virtual sources can be used to control the focus of the ultrasound through the skull, with better precision than using traditional phase-control focusing method. The central frequency and number of elements of the ultrasound transducer affect the focusing with similar patterns as in normal phase-control focusing. Moreover, the method allows for more than one target to be focused simultaneously, and the pressure amplitude at each focus can be adjusted adaptively. Key words Ultrasonic focusing, Time reversal, Transcranial ultrasound, Numerical simulation 收稿 ; 定稿国家重大科研仪器研制项目 ( ) 作者简介 : 苏畅 (1984- ), 女, 吉林延吉人, 博士, 副研究员, 研究方向 : 复杂介质声传播的数值模拟通讯作者 suchang@mail.ioa.ac.cn

2 年 9 月 1 引言人体组织中的精确超声聚焦在超声成像高强度聚焦超声 (High intensity focused ultrasound, HIFU) 治疗超声神经调控中都有重要作用, 对无骨骼遮挡软组织的超声聚焦已在临床中得到了广泛应用而在如 HIFU 经颅脑肿瘤治疗深部脑神经超声调控等情况下, 超声需穿透颅骨并在颅脑内聚焦由于颅骨不同位置具有不同的结构厚度密度声速和吸收系数, 特别是颅骨与其他组织声阻抗差异较大, 导致声波相位变化和能量吸收, 采用传统的聚焦控制方法会发生聚焦点偏移能量降低颅骨和皮肤升温等问题, 难以实现精确聚焦 [1] 因此必须提高对经颅超声声场的认识, 对传统的聚焦方法作出改进 Hynynen 等首先通过电子计算机断层扫描 (Computed tomography,ct) 了解颅骨厚度和密度方面的信息, 利用超声分层的衍射折射和反射的数字化模型计算时间延迟以实现相控聚焦 [2 3], 这种相位控制方法实现较为困难 199 年 Thomas [4] 等提出利用时间反转 (Time reversal,tr) 方法自适应地补偿超声经颅传播时的相位偏移以实现 [5] 经颅聚焦, 具有较高的聚焦精度 Fink 等对时间反转方法进行了全面的理论分析和实验要实现经颅超声时间反转需要在聚焦目标点处设置声源, 具体办法包括在颅脑内植入水听器作为主动声源, 以及植入反射体或利用高强度聚焦超声激发空泡作为被动声源等, 这些方法具有侵入性或对生物体有潜在危害 [] 而随着计算技术的发展, 利用高精度声场数值模拟方法, 可以在计算机上基于虚拟源模拟时间反转过程, 实现真正无损的经颅超声精确 [7] [8] 聚焦 Aubry 等和 Marquet 等利用三维 CT 扫描图像建模, 将基于虚拟源的数值模拟与动物实验结合, 实现了猴颅脑的经颅超声聚焦, 能够在聚焦点处获得较高能量 Leduc 等 [9] [1] 王意喆等利用基于虚拟源和数值模拟的时间反转方法消除伪焦点, 减轻颅骨升温, 进一步证明了时间反转方法在超声经颅聚焦中的良好应用前景为了进一步研究时间反转方法在经颅超声精确聚焦中的作用, 为超声经颅神经调控仪器的研发作理论指导, 本文在二维颅骨模型上, 利用有限差分法进行了基于虚拟源的时间反转过程模拟, 与传统的相控聚焦方法比较, 考察聚焦效果, 分析了超声频率换能器数量等参数对聚焦效果的影响, 并探索了利用时间反转实现多焦点同时聚焦的方法和聚焦效果 2 基于虚拟源的时间反转方法基于虚拟源的超声经颅时间反转聚焦过程如图 1 所示, 主要可分三个步骤 :(a) 从颅骨 CT/MRI 图像出发, 建立包含颅脑和换能器阵的物理模型 ; (b) 在聚焦目标处设置虚拟的点声源, 由该点发射超声脉冲或脉冲串, 利用数值模拟方法仿真计算该点声源激发的超声波经过颅骨传播到达换能器阵列的过程, 并记录换能器各阵元接收到的声压波形信号 ;(c) 将换能器处接收的波形在时间轴上反转, 反转后的波形再通过换能器阵列发出, 此过程可通过数值模拟或实验进行, 最终在目标点处再次聚焦图 1 中所采用的颅骨模型是从二维颅骨图像中获取的, 颅骨半径约为.1 m, 厚度从 5 mm 到 2 mm 不等颅骨白色区域声速为 33 m/s, 密度为 158 kg/m 3, 灰色区域声速 1 m/s, 密度 15 kg/m 3, 模型中其他区域均假设为水, 声速为 15 m/s, 密度 1 kg/m 3 换能器采用半圆形阵列, 半径为. m, 聚焦目标点选在颅内半圆形的圆心处 (a) (b) (c) 图 1 基于虚拟源的超声经颅时间反转聚焦过程 Fig. 1 Steps of transcranial ultrasound focusing with time reversal based on a virtual source

3 第 35 卷第 5 期苏畅等 : 基于虚拟源时间反转的经颅超声精确聚焦 379 二维直角坐标系下, 将模型中介质参数看作随空间坐标变化的函数 f(x, y), 弹性波在介质中的传播规律可以写成关于应力速度的一阶偏微分方程组 : [ ρ t t v x v y σ xx σ yy σ xy ] [ ] σ = xx / x + σ xy / y, σ xy / x + σ yy / y λ + 2µ λ v = x / x λ λ + 2µ v y / y + µ v x / y + v y / x g xx g yy, (1) 其中 v 为质点振动速度,σ 为应力张量,g 为外力源,ρ 为介质密度,λ,µ 为介质拉梅系数, 可表示为 λ = ρvp 2 2ρvs, 2 µ = ρvs, 2 其中 v p, v s 为别为介质纵横波波速在本文的经颅时间反转模型中, 仅考虑介质中的纵波, 只需设置横波速度 v s = 即可计算中声源由外力源 g ii 给出, 可由随时间变化的部分与随空间变化的部分组合而成, 例如空间中点 (x s, y s ) 处的点声源可写作 g ii = s (t) δ (x x s, y y s ), (2) 本文模拟中采用的声源波形 s(t) 为具有一定中心频率的正弦脉冲串, 中心频率为 5 khz 或 1 MHz 在换能器阵中, 由于单个阵元尺寸和阵元之间的间距与介质中声波的波长相比都比较小, 计算时可将每个阵元看作一个点声源将偏微分方程组 (1) 用中心差分格式近似, 利用空间 8 阶精度, 时间 2 阶精度的二维交错网格时域有限差分法求解当声源中心频率为 5 khz 时, 空间网格步长取.4 mm, 而声源中心频率为 1 MHz 时, 空间网格加密至.2 mm, 相当于水中一个声波波长内有 7 个网格点, 时间步长取. µs, 满足算法对稳定性和精度的要求为了防止计算区域截断边界反射影响声场计算的精度, 在边界处采用了非分裂的卷积完全匹配层 (NPML) 以模拟无限大介质 [11] 3 数值模拟结果和讨论 3.1 时间反转方法聚焦效果验证为了验证时间反转方法聚焦的效果, 以图 1 所示的数值仿真模型为例, 利用半圆形换能器阵列经颅骨向换能器阵列的几何中心聚焦换能器阵列共有 5 个阵元, 则相邻阵元间距约 1.1 mm, 小于水中声波波长的 1/2 首先由聚焦目标处作为声源, 激发的超声波经过颅骨由换能器各阵元接收, 将换能器各阵元接收的波形经时间反转处理后再次激发, 产生的声场快照如图 2 所示由于各阵元激发的波形相位不同, 在颅骨以外的波阵面形状并不规则, 而在穿透颅骨后波阵面为半圆形, 在目标点处实现同相位聚焦从声场幅度上看, 图 2(b) 中可看到大部分超声能量能够穿过颅骨, 少部分能量被颅骨反射作为对照, 我们计算了传统相控聚焦方法产生的声场假设每个阵元发出的超声在颅骨处垂直入射并穿透颅骨, 根据到达目标点路径上经过的介质声速计算所需的时间延迟, 使各阵元以一定的时间延迟发出相同幅度的超声信号计算得到几个时刻的声场快照如图 3 所示从模拟结果中可见穿过颅骨后的波阵面形状不规则, 各路信号到达目标点的时间不同, 会影响聚焦强度和精度同时, 从图 3(b) 中可看出, 被颅骨反射的超声能量所占比例较大进一步地, 取声场中各点声强峰值分布, 将两种聚焦方式作比较, 如图 4 所示时间反转方法得到的聚焦点更为集中, 而相控聚焦方式受颅骨非均匀性影响较大, 在焦点以外的位置产生了较强的噪声, 形成了形状不规则的较大焦斑 3.2 换能器频率阵元数对聚焦效果的影响根据以上的模型, 我们改变声源频率和换能器阵中阵元数量, 分别计算了基于虚拟源的时间反转方法和简单的相控聚焦方法向目标点聚焦的声场, 并比较聚焦效果模拟结果表明, 不同参数条件下时间反转聚焦的焦点位置均比较准确, 而换能器各阵元激发强度颅骨外表面反射的声压以及焦点处声压则有较明显的变化在经颅聚焦超声中希望

21 年 9 月 38 1 1 2 2 3 3 1 2 3 4 1 (a) t= µs 2 3 4 (b) t=9 µs 1 1 2 2 3 3 1 2 3 4 1 (c) t= µs 2 3 4 (d) t=15 µs 图2 时间反转聚焦法声场快照 Fig.

4 21 年 9 月 (a) t= µs (b) t=9 µs (c) t= µs (d) t=15 µs 图2 时间反转聚焦法声场快照 Fig. 2 Sound ﬁeld snapshots of ultrasound focusing with time reversal method (a) t= µs (b) t=9 µs (c) t= µs 1 2 (d) t=15 µs 图3 相控聚焦法声场快照 Fig. 3 Sound ﬁeld snapshots of ultrasound focusing with phase-control method

5 第 35 卷第 5 期苏畅等 : 基于虚拟源时间反转的经颅超声精确聚焦 (a) (b) 图 4 半圆形换能器阵超声经颅聚焦数值模拟得到的声强分布 Fig. 4 Sound intensity distribution from numerical simulation of ultrasound transcranial focusing with half circular transducer array 尽量减小换能器激发幅度, 提高目标焦点处的声压, 因此, 我们提出焦点声压增益来定量衡量经颅超声聚焦的效果, 即焦点处声压峰值与颅骨以外部分声压峰值的比值比较两种聚焦方式不同频率和阵元数情况下的焦点声压增益, 结果如表 1 所示表中数据表明 : (1) 各种参数条件下, 时间反转方法所得焦点声压增益均高于相控聚焦方法 (2) 相同聚焦方式, 同样换能器阵元数情况下, 随声源频率增高, 焦点声压增益降低, 即高频超声穿透颅骨的能力下降但低频超声形成的焦斑较大, 分辨率较低, 因此应结合实际需求合理选择超声频率 (3) 同样聚焦方式, 声源频率相同情况下, 增加阵元数目, 可提高焦点声压增益这是由于增加阵元数目后, 要在焦点达到同样的声压, 每个阵元所需的激发强度降低但要实现时间反转方法要求换能器各阵元具有独立的控制电路系统, 增加阵元数意味着大大增加系统的复杂程度, 因此实际中应根据需要作出取舍, 合理选择换能器阵列的阵元数量表 1 超声频率和换能器阵元数对焦点声压增益的影响 Table 1 Effect on focal pressure gain of ultrasound frequency and element number of the transducer 时间反转方法相控聚焦方法 5 khz 1 MHz 5 khz 1 MHz 21 阵元阵元多目标点聚焦在深部脑神经超声调控中, 往往需要对多个目标点同时进行超声刺激, 需要对发射的超声相位和幅度精确控制, 使超声同时聚焦至不同位置的多个目标点, 并尽量保证各目标点处声压的一致性用传统的相控聚焦方法很难实现这种多焦点同时聚焦, 只能采用由部分阵元分别控制各焦点的方式, 例如在图 5 的模型中, 要采用四组换能器线阵分别向 4 个目标点聚焦, 由下方上方换能器阵分别向 3 号目标点聚焦, 由左右两侧换能器阵分别向 1 2 号目标点聚焦此时对于每一个焦点, 所用的换能器阵孔径较小, 阵元数目较少, 因此必然导致聚焦精度和焦点声压受到影响数值模拟相控聚焦得到的声强分布如图 (a) 所示, 可见在各焦点均存在聚焦位置不集中的问题, 特别是受换能器孔径限制, 沿换能器阵法线方向焦斑较长 y/cm 图超声经颅多目标点聚焦模型 Fig. 5 Model of transcranial ultrasound focusing on multiple targets 24

382 21 年 9 月 y/cm 1..8..4 (a) 24 1 2 3 图 7 各焦点处的声压幅度比较 Fig. 7 Comparison of sound pressure amplitude at each focus 图 y/cm (b) 多目标点超声聚焦数值模拟得到的声强分布 Fig.

6 年 9 月 y/cm (a) 图 7 各焦点处的声压幅度比较 Fig. 7 Comparison of sound pressure amplitude at each focus 图 y/cm (b) 多目标点超声聚焦数值模拟得到的声强分布 Fig. Sound intensity distribution from numerical simulation of ultrasound focusing on multiple targets 采用基于虚拟源的时间反转方法, 可以很方便地在多个目标点处设置虚拟源, 由这些虚拟源同时激发声场, 所有换能器阵列同时接收并作时间反转处理时间反转法聚焦时, 每个焦点都由所有换能器阵列共同作用, 与上述相控聚焦法相比增大了起作用的孔径, 模拟计算得到图 (b) 的声强分布图中可以看到, 时间反转方法可以明显改善聚焦效果,4 个焦点能量均比较集中进一步探讨各目标点声压幅度一致性的问题在相控聚焦过程中, 由于颅骨不同位置厚度不同, 从不同方向聚焦时穿过颅骨到达目标处的超声能量差别较大, 导致各焦点之间声压幅度不同比较模拟计算所得 4 个焦点处的声压幅度, 如图 7 中星号所示, 最弱的 1 号焦点声压仅为号焦点声压的一半而在时间反转聚焦方法中, 也存在不同方向穿透颅骨超声能量的差异, 我们在时间反转过程中并没有对各阵元的波形幅度进行补偿, 因此不同位置的焦点处声压幅度也存在一定差别比较模拟计算所得 4 个焦点处的声压幅度, 结果如图 7 方块所示, 最强与最弱的焦点声压存在约 8% 的差距 24 利用基于虚拟源的时间反转方法, 可以迭代地改变虚拟源的强度, 自适应地调整焦点声压幅度具体过程是 : 首先设置多个强度相同的虚拟源, 数值模拟时间反转过程, 得到聚焦后目标点处的声压幅度 ; 将目标点处得到的声压幅度归一化并取倒数, 作为虚拟源激发波形的幅度 ; 由新的幅度不同的虚拟源同时激发声场, 重新进行时间反转针对图 5 的模型, 经过以上处理后, 得到 4 个焦点处的声压幅度如图 7 三角号所示, 具有较好的一致性 4 结论本文针对经颅超声聚焦中由于颅骨影响, 难以精确聚焦的问题, 采用基于虚拟源的时间反转方法实现超声的精确聚焦利用二维时域有限差分法进行数值模拟, 与传统相控聚焦方法作比较, 并分析了聚焦时各种参数对聚焦效果的影响结果表明, 基于虚拟源的时间反转方法可以克服颅骨非均匀性对超声聚焦的影响, 在颅内深部形成精确的聚焦, 聚焦强度和精度好于传统的相控聚焦使用时间反转方法时, 超声频率升高, 焦点声压增益降低, 而换能器阵元数目增加对聚焦效果有改善作用, 这些规律与相控聚焦方法中的规律一致利用基于虚拟源的时间反转方法可以同时向多个目标点聚焦, 并自适应地调整各目标点的虚拟源激发强度, 以保证各焦点处声压幅度的一致性这些结果可为经颅治疗和神经调控设备的研发提供支持和参考本研究仅采用了二维模型和线性声波方程, 基于三维模型的数值模拟研究正在进行中, 下一步研究工作还包括结合非线性声传播方程和生物组织热传导方程, 综合考虑生物体组织非线性和温升效

7 第 35 卷第 5 期苏畅等 : 基于虚拟源时间反转的经颅超声精确聚焦 383 应对经颅超声聚焦的影响等参考文献 [1] TER HAAR G. Acoustic surgery[j]. Phys. Today, 21, 54(): [2] SUN J, HYNYNEN K. Focusing of therapeutic ultrasound through a human skull: A numerical study[j]. J. Acoust. Soc. Am., 1998, 14(3): [3] CLEMENT G T, HYNYNEN K. A non-invasive method for focusing ultrasound through the human skull[j]. Phys. Med. Biol., 22, 47(8): [4] 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, 199, 43(): [5] FINK M A, MONTALDO G, TANTER M. Time reversal acoustics in biomedical and engineering[j]. Annu. Rev. Biomed. Engrg., 23, 5(1): [] KYRIAKOU A, NEUFELD E, WERNER B, et al. A review of numerical and experimental compensation techniques for skull-induced phase aberrations in transcranial focused ultrasound[j]. Int. J. Hyperther., 213, 3(1): 3 4. [7] AUBRY J F, TANTER M, PERNOT M, et al. Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans[j]. J. Acoust. Soc. Am., 23, 113(1): [8] MARQUET F, PERNOT M, AUBRY J F, et al. Noninvasive transcranial ultrasound therapy based on a 3D CT scan: Protocol validation and in vitro results[j]. Phys. Med. Biol., 29, 54(9): [9] 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]. J. Biomech. Sci. Eng., 2, 7(1): [1] 张千, 王意喆, 周文征, 等. 高强度聚焦超声经颅治疗温度场的数值仿真研究 [J]. 生物医学工程学杂志, 214, 31(): ZHANG Qian, WANG Yizhe, ZHOU Wenzheng, et al. Research on numerical simulation of temperature distribution during transcranial tumor therapy with high intensity focused ultrasound[j]. J. Biomed. Eng., 214, 31(): [11] 李佩笑, 林伟军, 张秀梅, 等. 常规分裂和非分裂完全匹配层吸收边界比较研究 [J]. 声学学报, 215, 4(1): LI Peixiao, LIN Weijun, ZHANG Xiumei, et al. Comparisons for regular splitting and non-splitting perfectly matched layer absorbing boundary conditions[j]. Acta Acustica, 215, 4(1):

/MPa / kg m - 3 /MPa /MPa 2. 1E ~ 56 ANSYS 6 Hz (a) 一阶垂向弯曲 (b) 一阶侧向弯曲 (c) 一阶扭转 (d) 二阶侧向弯曲 (e) 二阶垂向弯曲 (f) 弯扭组合 2 6 Hz

/MPa / kg m - 3 /MPa /MPa 2. 1E ~ 56 ANSYS 6 Hz (a) 一阶垂向弯曲 (b) 一阶侧向弯曲 (c) 一阶扭转 (d) 二阶侧向弯曲 (e) 二阶垂向弯曲 (f) 弯扭组合 2 6 Hz 31 3 Vol. 31 No. 3 218 9 Journal of Shijiazhuang Tiedao University Natural Science Edition Sep. 218 1 1 2 1 2 1 1. 543 2. 543 U462. 3 217-2 - 16 A 295-373 218 3-63 - 6 1-4 5-7 8-11 1 11 11 398 mm 86 mm