1.   Introduction

Transcranial focused ultrasound has the advantages of non-invasiveness, safety and large penetration depth. ^{[ 1 ]} It has great potential for development. Ultrasound can maintain a high penetration depth through the skull, so that it can directly affect the deep brain tissue and focus the thermal effect produced by ultrasound. ^{[ 2 ]} Therefore, compared with the traditional surgical treatment across the skull, transcranial focused ultrasound technology can achieve surgical treatment without craniotomy, greatly reducing the risk of infection and bleeding, and has broad application scenarios in the biomedical field. For example, at low intensity, transcranial focused ultrasound can reversibly open the blood-brain barrier without tissue ablation. ^{[ 3 , 4 ]} At high intensity, transcranial focused ultrasound can produce tissue ablation, and the outcome of tissue ablation is related to the time of tissue exposure to high intensity focused ultrasound. ^{[ 5 ]} .

However, when ultrasound passes through the skull, it will bring strong ultrasound scattering and attenuation due to the nature of the skull itself and the complexity of the structure. ^{[ 6 ]} , causing severe distortion and distortion of the transcranial sound field ^{[ 7 , 8 ]} In recent years, a technique for generating arbitrarily complex sound fields through a monolithic acoustic holographic lens has emerged. ^{[ 9 ]} Holographic technology is widely used in all kinds of energy field regulation because of its accuracy, real-time, low power consumption and other characteristics. This method stores the phase or amplitude of the required wave surface in space. ^{[ 10 ]} , so that when the holographic lens is illuminated by a suitable energy source, the wave surface is reconstructed by interference. Today, modern computer holography ^{[ 11 ]} It is no longer necessary to record the phase or amplitude of the wave surface from the real physical scene, but to calculate the phase or amplitude of the wave surface in the numerical simulation and reconstruct it based on it. In the technique of generating complex sound fields through a monolithic acoustic holographic lens, the phases or amplitudes of the ultrasound are calculated by computer numerical simulation and stored in the acoustic hologram. The acoustic holographic lens is placed on the ultrasonic propagation path, and the complex sound field can be generated by changing the phase of the ultrasound. ^{[ 12 ]} At present, this technology has been proved to be able to overcome the distortion and distortion caused by ultrasound transmission across the skull and produce focus at any position in the skull. ^{[ 13 ]} Compared with the current mainstream transcranial focusing method based on phased array. ^{[ 14 , 15 ]} This technology does not rely on complex electronic equipment, and has the advantages of low cost and convenient operation. ^{[ 16 ]} .

However, the existing transcranial focusing methods based on acoustic holographic lenses still face serious application problems. A single piece of acoustic holographic lens can only reconstruct a single focused sound field, which is lack of flexibility in practical application scenarios. Therefore, a design method of multi-frequency acoustic holographic lens for transcranial focusing is proposed in this paper. By extracting the effective information from two pieces of acoustic holographic lenses focusing at different positions at different working frequencies, it is integrated into a single piece of acoustic holographic lens to achieve the goal based on.

2.   Theoretical basis

2.1   Propagation of sound field in inhomogeneous medium

Due to the complexity of the nature and structure of the skull itself, ultrasonic waves will experience severe attenuation and scattering when they pass through the skull. ^{[ 17 ]} , resulting in severe sound field distortion ^{[ 18 ]} In order to accurately simulate the process of transcranial sound field propagation, the k-wave in MATLAB is used in this paper. ^{[ 19 ]} Open source toolbox that uses perfectly matched layers to absorb boundary reflections with the help of k Spatial pseudospectral method ^{[ 20 , 21 ]} To describe the propagation of ultrasound in a heterogeneous medium such as the skull. k The spatial pseudospectral method transforms the sound field into an integral equation in frequency-domain space to solve the sound field distribution. Compared with the traditional method, it shows higher computational efficiency in dealing with the case of inhomogeneous medium, and has been widely used in the field of acoustic simulation. This method describes the sound propagation in inhomogeneous medium by solving the following equation ^{[ 22 ]} :。

Among, $u$ Stands for the acoustic particle velocity, $d$ Is the acoustic particle displacement, $p$ Representing sound pressure, $\rho$ Is the acoustic density, ${\rho _0}$ Is the density of the medium, ${c_0}$ Is the isentropic sound speed, through the linear operator L Frequency-dependent acoustic absorption and dispersion are introduced ^{[ 23 ]} :。

Among, $\tau = - 2{\alpha _0}c_0^{\gamma - 1}$ , $\eta = - 2{\alpha _0}c_0^\gamma \tan ({\text{π}}\gamma /2)$ Absorption and dispersion proportional coefficients, respectively, ${\alpha _0}$ Is the power-law prefactor, $\gamma$ Is the power law absorption index.

2.2   Design principle of multifrequency acoustic holographic lens

At present, the time reversal method is a mainstream design method of transcranial focusing acoustic holographic lens, which uses the time reversal invariance of acoustic propagation ^{[ 24 ]} The counter-propagating field is obtained by numerical simulation by setting a series of virtual sound sources at the focusing position. ^{[ 25 ]} Firstly, a virtual sound source is set inside the skull, and the ultrasonic signal from the virtual source is recorded on the holographic surface outside the skull model. Then, the signal recorded by the sensor is Fourier transformed, and the phase at the working frequency is extracted as the acoustic hologram. In order to realize multi-frequency sound field control based on a single acoustic holographic lens, this paper attempts to integrate the effective information from two groups of acoustic holograms driven at different frequencies and focused at different positions. Specifically, by using the time reversal method, ${p_{\text{f}}}(x, y)$ , followed by Fourier transforming the received signal and extracting the phase at the operating frequency as an acoustic hologram, denoted ${H_1}(x, y)$ And ${H_2}(x, y)$ , define the cost function as follows ^{[ 26 ]} :。

Among, $E(|{p_{\text{f}}}(x, y)|)$ A pixel on that holographic plane driven at a certain frequency $(x, y)$ The average sound pressure value received at, ${\text{std(|}}{p_{\text{f}}}(x, y)|)$ A certain pixel point on that holographic plane driven at a certain frequency $(x, y)$ The standard deviation of the sound pressure value received at. $S$ Is the smoothing term found with a second order accurate finite difference approximation of the Laplace operator, $\alpha$ And $\beta$ Is the weight parameter. First, the multi-frequency hologram $H(x, y)$ Initialize with zero and then traverse one by one $H(x, y)$ For each pixel in, set it to ${H_1}(x, y)$ Or ${H_2}(x, y)$ The value of the corresponding pixel in. Specifically, ${H_1}(x, y)$ And ${H_2}(x, y)$ Corresponding ${p_{\text{f}}}(x, y)$ Substituting (resp. 5 ), calculate and compare the cost functions of the two, and assign the higher cost function to $H(x, y)$ Repeat this operation until all pixels are updated. Repeat the iteration and recalculate after each iteration. $S$ The iteration is terminated when the number of pixels changed from the previous iteration is less than 0.5% of the total pixels of the surface. The initial iteration time is $\beta$ = 0, each computation updates after iteration $\beta$ , make $E(|{p_{\text{f}}}(x, y)|) = - 0.5\beta S$ , variable $\alpha$ Always set to 0.3. After obtaining the multi-frequency acoustic hologram, refer to He et al. ^{[ 27 ]} Through 6 The extracted phase information is converted into the height value corresponding to each pixel on the surface of the acoustic holographic lens. $\xi (x, y)$ :。

Among, ${c_{{\text{water}}}}$ And ${c_{{\text{lens}}}}$ The sound velocities in water and in the holographic lens material are 2495 m/s and 1500 m/s, respectively. $f$ Is the average of the two frequencies corresponding to the multi-frequency acoustic holographic lens, which is 1.5 MHz.

2.3   量化指标

为了评估基于本文提出方法设计的多频声全息透镜的声场聚焦性能, 将生成的声全息透镜分别在仿真中和实验中进行测试. 用不同频率的平面波进行经颅聚焦实验, 并测量聚焦平面上的声压值, 并通过对比聚焦面上的焦点与预设的目标焦点之间的距离来验证基于该方法生成的多频声全息透镜的声场聚焦准确性. 除此之外, 分别计算在两种频率下重建的声场分布与预设的目标声场分布之间的峰值信噪比(PSNR), 从而评估基于该方法的声场聚焦质量, 验证所提出的方法的性能. PSNR的计算公式如下^[28]:

其中MAX代表声压图像中像素点的最大值, MSE代表均方误差. PSNR值越高, 说明声场聚焦质量越好.

3.   数值仿真与实验

3.1   多频声全息透镜的设计

在数值仿真中, 设置仿真网格大小为0.25 mm, 库朗数为0.2, 数值仿真通过MATLAB中的k-wave工具包实现. 具体来说, 利用时间反演法, 在聚焦平面目标焦点的位置设置点源, 在声全息透镜所在的平面接收时域信号并提取在工作频率下的相位来获取声全息图. 分别获取两组目标焦点在1 MHz和2 MHz频率下的声全息图, 如图1所示.

图 1  目标焦点、时间反演法生成的声全息图以及多频声全息图

Fig. 1.  Target foci, acoustic holograms generated by the time inversion method, and multi-frequency acoustic holograms.

下载: 全尺寸图片 幻灯片

3.2   仿真和实验平台设置

在获得多频声全息图后, 基于(6)式可以获得多频声全息透镜每个像素的高度, (6)式的频率参数设置为两种频率的均值, 即1.5 MHz, 制成的声全息透镜如图2所示.

图 2  用于经颅聚焦的多频声全息透镜　(a)多频声全息透镜的3D视图; (b)多频声全息透镜俯视图; (c)多频声全息透镜的正视图

Fig. 2.  Multi-frequency acoustic holographic lens for transcranial focusing: (a) 3D view of the multi-frequency acoustic holographic lens; (b) top view of the multi-frequency acoustic holographic lens; (c) front view of the multi-frequency acoustic holographic lens.

下载: 全尺寸图片 幻灯片

为了验证该全息透镜的聚焦性能, 建立了经颅聚焦的仿真和实验平台. 在仿真和实验中, 颅骨、全息面、聚焦平面之间的相对位置设置如图3所示. 仿真平台利用MATLAB中的k-wave工具包搭建, 仿真参数与3.1节一致. 实验平台以针式水听器(ONDA-HNP-0400)为核心, 如图3所示. 多频声全息透镜和颅骨仿体由3D打印技术制造得到. 材料选用Verowhite, 其中声速为2495 m/s^[29], 密度为1175 kg/m³. 颅骨模型的制造分辨率为750 μm. 颅骨模型上配有支架, 可与声全息透镜一同安装在超声换能器上. 使用信号发生器(KEYSIGHT MXG N5182 B), 发射1 MHz或2 MHz的正弦脉冲, 并通过线性射频放大器连接到超声换能器. 整套实验装置被浸没在充满水的水箱中. 利用针式水听器扫描虚拟源所在的平面, 扫描区域大小为50 mm×50 mm, 步长为0.5 mm. 通过对获得的信号进行滤波处理, 采集到时变信号并计算其能量, 以绘制声压图.

图 3  经颅聚焦实验平台

Fig. 3.  Experimental platform for transcranial focusing.

下载: 全尺寸图片 幻灯片

4.   结果分析与讨论

经过仿真和实验, 基于本文设计的多频声全息透镜得到一系列经颅聚焦声场. 基于此, 通过2.3节所提出的量化指标对聚焦结果进行评估, 以验证该方法的可靠性. 图4为对应的仿真与实验结果. 图4(a), (c)分别为在1 MHz和2 MHz的平面波激发下, 在聚焦平面附近的纵向声场, 图4(b), (d)分别为在1 MHz和2 MHz的平面波激发下, 目标焦点的位置与在仿真和实验中在聚焦平面上的横向声场. 图4(e), (f)分别为在1 MHz和2 MHz的平面波激发下, 仿真与实验中聚焦平面上的横向声场在x = 7.3 mm的虚线处的声压对比. 根据图4结果可知, 通过本文所提出的方法生成的多频声全息透镜在不同频率的超声激发下, 能够在仿真中在聚焦平面的不同位置处生成精准的声聚焦点, 并且通过实验验证该多频声全息透镜在实际应用场景下依旧具备良好的性能, 具有灵活的应用性.

图 4  多频声全息透镜的仿真与实验结果　(a) 1 MHz仿真中, 焦点附近的纵向声场; (b)1 MHz激发下的仿真与实验重建的聚焦平面与目标焦点; (c) 2 MHz仿真中, 焦点附近的纵向声场; (d) 2 MHz激发下的仿真与实验重建的聚焦平面与目标焦点; (e) 1 MHz激发的情况下, 在仿真和实验中, 沿着x = 7.3 mm的虚线上的声压对比; (f) 2 MHz激发的情况下, 在仿真和实验中, 沿着x = 7.3 mm的虚线上的声压对比

Fig. 4.  Simulation and experimental results of multi-frequency acoustic holographic lenses: (a) Longitudinal acoustic field near the focal point in 1 MHz simulation; (b) focusing plane and target foci reconstructed by simulation and experiment under 1 MHz excitation; (c) longitudinal acoustic field near the focal point in 2 MHz simulation; (d) focusing plane versus target foci reconstructed from simulation and experiment under 2 MHz excitation; (e) sound pressure comparison along the dashed line at x = 7.3 mm in simulation and experiment for the case of 1 MHz excitation; (f) sound pressure comparison along the dashed line at x = 7.3 mm in simulation and experiment for the case of 2 MHz excitation.

下载: 全尺寸图片 幻灯片

此外, 为了评估基于本文提出方法的超声聚焦质量, 在数值仿真和实验中, 皆对聚焦声场与目标焦点之间的PSNR进行了计算, 如表1所示. 相对而言, 在仿真和实验中, 较高的PSNR证明该声全息透镜可以在不同频率的超声激励下在聚焦平面重建出不同的高质量焦点.

表 1  聚焦平面的声场与目标图像之间的PSNR

Table 1.  PSNR between the sound field at the focal plane and the target image

频率/MHz	PSNR(仿真)	PSNR(实验)
1	32.16	27.48
2	40.18	32.33

下载: 导出CSV 

| 显示表格

为了验证通过该方法制作的声全息透镜在多焦点情况下, 仍旧可以生成高质量的焦点, 在1 MHz场景下的聚焦平面中设置多个焦点, 设计了第2片用于经颅聚焦的声全息透镜来进行重复实验以验证该方法在多焦点的情况下仍旧具备可行性. 仿真和实验结果如图5所示. 图5(a), (c)分别为在仿真中, 在1 MHz和2 MHz激发下焦点平面的纵向声场. 图5(b), (d)分别为在1 MHz和2 MHz的激发下, 仿真和实验中在聚焦平面处形成的焦点和目标焦点. 图5(e), (f)分别为在1 MHz和2 MHz的激发下, 仿真与实验中聚焦平面上的横向声场在y = 7.3 mm的虚线处的声压对比. 由结果可知, 即使在多点情况下, 用不同频率的超声激发多频声全息透镜, 仍能在不同频率超声激发下, 生成准确的经颅聚焦点.

图 5  多频声全息透镜在多焦点情况下的仿真与实验结果　(a) 1 MHz仿真中, 焦点附近的纵向声场; (b)1 MHz激发下的仿真与实验重建的聚焦平面与目标焦点; (c)在2 MHz仿真中, 焦点附近的纵向声场; (d) 2 MHz激发下的仿真与实验重建的聚焦平面与目标焦点; (e)1 MHz激发的情况下, 在仿真和实验中, 沿着y = 7.3 mm的虚线上的声压对比; (f)2 MHz激发的情况下, 在仿真和实验中, 沿着y = 7.3 mm的虚线上的声压对比

Fig. 5.  Simulation and experimental results of multi-frequency acoustic holographic lenses in the case of multiple focal points: (a) Longitudinal acoustic field near the focal point in the 1 MHz simulation; (b) focused plane and target foci reconstructed by simulation and experiment for the 1 MHz excitation; (c) longitudinal acoustic field near the focal point in the 2 MHz simulation; (d) focused plane and target foci reconstructed by simulation and experiment for the 2 MHz excitation; (e) sound pressure comparison in the 1 MHz excitation case in the simulation and experiment along the y = 7.3 mm dashed line; (f) sound pressure comparison along the dashed line at y = 7.3 mm in simulation and experiment for the case of 2 MHz excitation.

下载: 全尺寸图片 幻灯片

同样, 计算了在数值仿真和实验中在聚焦平面重建的聚焦声场与目标图像之间的PSNR, 如表2所示. 在仿真中和实验中, 较高的PSNR证明该多频声全息透镜即使在多焦点的情况下, 仍旧能在不同频率的超声激发下生成高质量的焦点.

表 2  聚焦平面的声场与目标图像之间的PSNR

Table 2.  PSNR between the sound field at the focal plane and the target image.

频率/MHz	PSNR(仿真)	PSNR(实验)
1	29.39	23.30
2	39.89	32.17

下载: 导出CSV 

| 显示表格

综上所述, 本文所设计的多频声全息透镜, 在单焦点和多焦点情况下均能在相应频率的激发下生成高质量焦点, 其物理原因在于通过对全息透镜相位分布的相干叠加设计实现对复合声波束干涉和衍射效应的特殊控制. 在1 MHz的超声激发下进行的数值仿真中, 两个多频声全息透镜在聚焦平面上分别形成了单个焦点和多个焦点. 相较于单焦点, 多焦点中的能量较高的焦点和较低的焦点分别为单焦点能量的0.63倍和0.44倍, 主要原因在于多焦点聚焦需要通过控制不同区域的相位信息来使多个焦点处的波前相干叠加, 这导致能量在多个焦点之间分散, 从而引发能量损失. 此外, 随着激发频率从1 MHz增大到2 MHz, 生成的焦点的分辨率明显提高, 同时PSNR也得到了显著提升, 这一现象可归因于频率的增加导致波长缩短, 同时全息透镜相对于波长而言的数值孔径变大, 从而使声波的聚焦点更加精确, 体现了较短的波长在提高聚焦分辨率和聚焦质量方面的优势.

5.   结　论

针对目前用于经颅聚焦的声全息透镜只生成单一声场, 缺乏应用灵活性这一问题, 在本文中, 我们提出了一种用于经颅聚焦的多频声全息透镜的设计方法, 通过提取不同频率下, 聚焦在不同焦点位置的声全息透镜的有效信息, 生成用于经颅聚焦的多频声全息透镜. 数值仿真与实验结果证明, 无论在单焦点还是多焦点任务下, 该透镜均能在不同频率的平面波激励下, 克服颅骨造成的声波散射问题, 能够跨骨形成对应的精确且高质量的焦点. 在一定程度上提高了声全息透镜在经颅聚焦这一领域的应用灵活性. 但是值得一提的是, 该工作目前只应用于颅骨仿体中, 在真实的颅骨中需要考量更复杂的情况, 比如松质骨和皮质骨对聚焦精度产生的影响, 这些将在之后的工作进行研究.

下载: 全尺寸图片 幻灯片