时间反演聚焦经颅磁声电刺激仿真与实验研究*

张 帅 黄明辉 焦立鹏 郭 梁 周振宇 徐桂芝

(1 河北工业大学 省部共建电工装备可靠性与智能化国家重点实验室 天津 300130)

(2 河北工业大学 天津市生物电工与智能健康重点实验室 天津 300130)

0 引言

神经调控技术是开展神经科学研究与神经系统疾病临床治疗的基本手段。基于电磁感应原理，通过采用电场或磁场等非侵入式方式刺激脑神经组织，来达到调节脑功能的目的，包括经颅磁刺激(Transcranial magnetic stimulation，TMS)[1]、经颅直流电刺激(Transcranial direct current stimulation，tDCS)[2]和经颅超声刺激(Transcranial ultrasound stimulation，TUS)等方法[3]。经颅磁声电刺激(Transcranial magnetic acoustic electrical stimulation，TMAES)是一种新型无创脑神经调控技术，将超声波和静磁场作用于神经组织产生电流，从而对神经组织进行刺激，刺激区域由超声波在神经组织的聚焦区域决定，同时静磁场在组织中的衰减远小于交变磁场，因此TMAES 是一种高空间分辨率的深部有效神经刺激。

Norton[4]通过计算麦克斯韦方程组，理论上证明静磁场和超声场可以产生足够强的感应电流对大脑进行刺激。杨少华[5]对静磁场中带电凝胶受超声作用而产生的电信号以及超声回波信号做了频谱分析，指出二者中心频率一致，初步验证了该方法的可行性。李慧雨等[6]使用超声换能器向静置于磁场中的导电样本打入超声，并测量沿不同方向的声压信号以及电信号，分析了施加声场和感应电场之间的一致性。Yuan 等[7]和刘世坤等[8]分别通过分析刺激后大鼠局部场电位信号的变化以及动物行为学实验，证明了磁声耦合电刺激方法可形成有效刺激。张帅等[9-10]通过建立包括刺激深度、刺激强度、聚焦面积、焦域长度以及电场强度梯度系数在内的评价指标，对经颅磁声电刺激感应电场聚焦性能进行了有效评价；并且基于H-H 神经元模型，仿真了TMAES下不同超声参数对神经元放电模式的影响，对于具体神经调控和得到最优TMAES 参数提供理论指导。

相控阵超声换能器可实现高空间分辨率、深层次的有效刺激，但是颅骨对声场的散射和衍射严重影响了焦点区域的声场[11]。骨骼的声阻抗约是软组织的10 倍[12]。Yuldashev 等[13]在阵列换能器聚焦模拟中发现了随声压变化的焦点偏移现象。为了实现更加精准的超声聚焦，Fink 等[14]将时间反演法引入了声学领域，系统理论分析的同时进行了大量的水声实验，实现了在空间和时间上的同步聚焦。高翔等[15]对时间反演法进行了深入的研究并进行改进，应用于流体、固体和分层介质，进行裂纹检测和目标定位。该方法在生物医学和磁声成像等方面也有广泛应用[16-17]。苏畅等将时间反转法用于经颅超声聚焦，证明时间反转法相比于传统的时间延时方法可以使聚焦更准确[18-20]。

现有经颅磁声电刺激相关文献大多使用超声换能器加声准直器或透镜约束方式进行神经刺激，所建立的声场和感应电场的聚焦性仍有待提高。本文采用相控阵超声换能器，通过在颅脑内设置虚拟点声源仿真模拟超声传播过程，记录各阵元接收到超声信号的时间，按照后到先发原则进行反演聚焦。相比于传统相控聚焦，可以准确地在预设点实现聚焦。通过搭建实验平台，利用水听器检测到的数值绘制聚焦点处的归一化声压图，证实了时间反演法可以解决焦点偏移，实现自适应聚焦，并对声场和感应电场分布一致性进行验证，进一步证实了时间反演用于经颅磁声电刺激的可行性。

1 原理

1.1 TMAES原理

神经组织样本内部带电粒子受超声波作用而振动，振动的粒子在磁场B0中切割磁感线产生与超声方向垂直的洛伦兹力，导致正负电荷向相反方向移动，进而产生感应电场E，经颅磁声电刺激的原理如图1所示。超声作用下样本质点沿y方向的振动速度为v(y，t)，q为样本质点所带电荷量，电荷受到的洛伦兹力为[4]

图1 TMAES 原理图Fig.1 TMAES schematic

神经组织样本内部产生的感应电场强度E为

组织样本质点的振动速度与质点处有效声压p(y)之间的关系为

式(3)中，ρ和c分别为神经组织样本密度、超声波在神经组织样本的传播速度。将式(3)带入式(2)中，可得振动质点处有效感应电场强度为

1.2 传统相控聚焦

经颅磁声电刺激需要聚焦声束来产生感应电场以进行神经刺激，传统相控聚焦方法是计算各阵元时间延时来使声束偏转聚焦，延时时间的计算取决于中心频率、尺寸大小、焦点位置以及焦距，选取探头中心为参考点，当阵元数N为奇数时，各阵元相对参考点的时间延时为

其中，c是超声传播速度，l为焦距，d为探头阵元间距，θ为偏转角度，n即第n个阵元(n=0，±1，···，±(N -1)/1)，当阵元数为偶数时，公式(5)中n则改为(n+0.5)。

1.3 时间反演法

为解决焦点偏移现象，本文采用时间反演法进行声束聚焦，在预设焦点处设置虚拟点源，各阵元接收声波，按照后到先发的原则，使得发射时各阵元拥有准确的延迟发射时间，可以精准地聚焦在理想位置，原理如图2所示。利用互相关算法[21]可以准确地计算出各个阵元的时间延时，具体的实施步骤为：第n个阵元接收到的信号为Xn(t)，经过时间反转后，其对应的信号为Xn(T -t)。将此阵元的信号重新定义Xnr(t)，在每一个反转信号添加一个时间常量t0，各阵元时间反转信号为

图2 时间反演聚焦原理图Fig.2 Schematic diagram of time reversal focusing

其中，Kn为各阵元的补偿系数，各阵元激励信号强度不一样，一般阵元中心处权重大于边缘处阵元的权重，可提高主瓣能量，减少旁瓣数量和抑制栅瓣。

2 数值模拟过程与结果

2.1 模型建立

本文仿真所用模型数据皆从二维颅脑图像中获得。其中，颅骨的厚度为5 mm 到15 mm 不等，大脑侧面长152 mm，正面长131 mm，颅骨的声学参数密度ρ1和声速c1分别为1912 kg/m3、2300 m/s，大脑的声学参数密度ρ2和声速c2分别为1030 kg/m3、1550 m/s[22]。如图3(a)所示，设置32 个线阵元模拟超声换能器，与颅骨之间设置纯水，各阵元中心频率0.5 MHz，阵元宽度a为λ/2，阵元间距d为3λ/4，阵元长20 mm，虚拟点源置于换能器中心点正下方40 mm 处，四周设置软声场边界模拟无限大介质用来吸收出射波，防止反射的声波影响计算精度。为了防止大脑神经组织受连续脉冲波刺激产生热积累现象引发不良反应，所以发射周期性正弦波超声脉冲群，每个脉冲群含有的基波周期数Nc为100，声脉冲群重复频率(Pulse repetition frequency，PRF)为1000 Hz，脉冲重复数(Pulse repetition number，PRB)为200[23]。空间网格步长dx= dy= dz= 0.4 mm，步进时间为0.5 μs，以满足精度和稳定性的要求。

2.2 点源发射以及颅内磁场分布

利用图3(a)仿真模型进行仿真，声波向四周发射，图3(b)为脑组织内虚拟点声源发射单正弦脉冲传播过程中的二维截面图。超声经大脑、颅骨和水域到达换能器各个阵元。由图3(b)可见，颅骨对超声的能量衰减、传播路径有较大的影响。由于颅骨声阻抗较高，只有小部分的超声能量穿过颅骨，并因折射使传播路径发生变化，部分超声能量通过反射回到颅内。

图3 仿真模型与点源发射声波Fig.3 The simulation model and point source emits sound waves

记录超声波到达换能器各阵元时的声压幅值以及时间并绘制归一化声压图，图4为17号～32号阵元接收到的归一化声压信号图(将图2中换能器阵列右边第一个定义为1 号阵元)。图4中左边出现第一个峰值的曲线为17号阵元的归一化声压图，最后一条为32号阵元，因为在非均匀介质中传播路径的不同，距离虚拟点源近的阵元相比于较远的阵元接收到的声压幅值较高，声波到达时间较短。记录每个阵元的归一化声压峰值所对应的时间，并将32号阵元作为参考阵元，计算其余阵元与参考阵元的时间差，即为各个阵元的延迟发射时间。

图4 部分阵元接收信号Fig.4 Some elements receive signal

永磁体提供经颅磁声电刺激所需的静磁场，图5为基于颅脑结构的TMAES 系统模型，颅骨两边各放置一块直径80 mm、厚20 mm 的柱形永磁体。图6为预设点源处沿永磁体中轴线上磁通密度分布曲线，曲线对称分布，z轴0 mm 处为极小值0.5 T。将磁通密度降至最大值95%的区域近似看成磁场均匀分布区域，则分布区域为-12～12 mm，超声聚焦区域的大小为毫米级别的，所以在聚焦点处的磁场可近似看成0.5 T 的均匀磁场，满足磁声电刺激的基本要求。

图5 TMAES 系统模型Fig.5 TMAES system model

图6 磁通密度模Fig.6 Norm of magnetic flux density

2.3 仿真结果分析

超声换能器发射的声波经水域、颅骨和大脑在预设点上聚焦，图7(a)为传统相控聚焦声场切面图，图中红色部分为聚焦区域，蓝色为非聚焦区域，单位为声压Pa。焦点形状狭长不规则，选用焦点横向和纵向声压峰值分别衰减到0.707倍时的区域为焦域，所得焦点纵向分辨率28.4 mm，横向分辨率5.7 mm，实际焦点位置相比于预设焦点径向偏移4 mm，轴向偏移2.5 mm。因为颅骨的厚度不均匀，使得在聚焦区域左边有长约8 mm、宽约2 mm 的焦域，可对非聚焦区域形成刺激。图7(b)为运用时间反演法的聚焦声场示意图，焦点形状规则分布，呈梭形，焦点纵向分辨率23.2 mm，横向分辨率4.6 mm，相比于传统相控聚焦，聚焦区域减小，实际焦点位置相比于预设焦点径向偏移1 mm，轴向偏移0.5 mm，实际焦点与预设焦点吻合度较高，说明时间反演法在分层的非均匀介质中具有良好的聚焦效果。

图7 传统相控聚焦和时间反演聚焦Fig.7 The traditional phased array focusing and time reversal focusing

由公式(4)可得，当磁感应强度恒定为0.5 T时，神经组织密度和超声传播速度为定值，感应电场强度数值变化规律与声压大小呈线性关系，声场分布可近似为电场分布。图8中，实线为运用时间反演法焦点径向声压最大处轴向截线电场分布，虚线为传统相控聚焦焦点径向声压最大处轴向截线电场分布，两种聚焦方式下测得的声压最大值分别为1.44 MPa 和1.32 MPa，经公式(4)计算得，运用时间反演法刺激靶点处的感应电场峰值为0.48 V/m，大于传统相控聚焦的感应电场峰值0.44 V/m，聚焦强度略有增加。传统相控聚焦的聚焦区域与所设虚拟点源位置有一定的偏差，难以实现精准聚焦，而运用时间反演法的聚焦区域与虚拟点源的位置基本吻合，这是由于颅骨折射声波使传播路径改变，进而导致焦点偏移，同时提高声束聚焦精度和提高焦点处感应电场峰值。

图8 电场径向分布曲线Fig.8 Radial electric field distribution curve

时间反演法同样适用于多点聚焦，通过控制两组或者多组阵元可以实现两个或者多个焦点。多数神经类疾病的产生并不是孤立的大脑区域异常的结果，因此研究不同区域的多焦点刺激同样具有重要应用价值[24]。图9为运用时间反演生成的两点聚焦声场分布图，将32 个阵元分成两组，预设的两个虚拟点源先后发射声波，分别记录两组阵元对应预设点源接收到的信号，各阵元接收到的信号幅值不同，通过调节各个阵元发射声压的幅值，可以改变各焦点处声压能量，因为颅骨厚度不均匀，在超声传播过程中对超声能量的阻挡程度不同，传播路径也不同，导致两个焦点的形状和大小不一样。

图9 两点聚焦Fig.9 Two-point focus

3 实验验证

3.1 实验设计

本文设计仿体实验来对比传统相控聚焦和时间反演聚焦的效果，利用本实验室Verasonics 开放式多通道超声研究平台实现焦点可视化，可以对刺激位置准确定位。如图10所示，在水箱中放置由CT模型3D 打印而成的仿体颅骨，选择与真实颅骨声学特性相近的聚酯材料。在矢状缝与冠状缝交点正下方20 mm 处固定64 阵元相控阵超声探头(P4-2v探头，Verasonics，美国)。在颅骨仿体内部距离超声探头中心正上方40 mm 处放置水听器(HNR-0500，ONDA，美国)，水听器稳定粘合于千分尺，用于检测颅骨内部的声压信号。因为水听器的输出信号十分微弱，所以需要辅助使用声压放大器(5660B PREAMP，OLYMPUS，日本)放大超声声压信号，将60 dB 前置声压放大器连接在示波器(TDS3014，Tektronix，美国)上，探头由Verasonics 系统控制各阵元时间延时发射连续超声脉冲波，以1 mm 为单位在水平方向移动水听器，并检测声压。

图10 实验平台Fig.10 Experiment platform

此外，又设计检测实验来考察声场与感应电场之间的关系。如图11 所示，在静磁场中放入与超声换能器处于同一水平面的导电样本，将金属导线放置于导电样本中，并与超声换能器探头表面相距40 mm，在相对位置不变的情况下，控制Verasonics系统超声相控阵的焦距以2 mm 为单位从20 mm递增至60 mm，金属导线连接差分放大器放大测量到的电信号，使用万用表记录数值，并使用水听器检测与金属导线处于同一区域的超声信号，将测得的两种信号归一化处理并对比。

3.2 实验结果

本实验使用的超声相控阵参数如下：阵元数N为64，阵元宽度a为0.25 mm，阵元自身长度h为14 mm，阵元间距d为0.3 mm，相控阵总长l为18.9 mm，如图12 所示。当阵元间距d远小于超声波长λ时，波束指向性较差[25]，为得到形状规则的焦点，设置各阵元基波频率为1 MHz。预先设定这些参数进行仿真，以获得阵列延迟传输所需时间。图13(a)为仿真模式下使用时间反演法得到的阵元中心处纵向截面声场分布图，焦域纵向分辨率约为21.3 mm，横向分辨率约为1.5 mm。图13(b)为聚焦实验中计算机屏幕显示的GUI 控制界面，其中，弧形区域为颅骨轮廓，图中的一个亮斑即为水听器尖端，移动水听器，使示波器显示的电压值最大，即可视为聚焦区域的中心点，此外在探头发射声波聚焦的过程中不考虑剪切波的影响[26]。

图12 探头组成示意图Fig.12 Schematic diagram of probe composition

图13 单点聚焦和实验界面Fig.13 Single point focus and experimental interface

图14中实线为仿真模式下y轴等于0 mm处横向归一化声压图，由图14 可见声压最大值位于z轴0 mm 处，基本没有发生焦点偏移现象，两种聚焦方式出现声压峰值的位置不同，反演聚焦时在z= 0 mm 处达到峰值，相控聚焦则在z=-1 mm处达到峰值，而用时间反演法根据各采样点绘制的归一化声压图与仿真模式下绘制的声压图分布大致相似，相关系数为0.823，证明了时间反演法能够解决焦点偏移现象，实现精准聚焦。

图14 两种聚焦方式声场对比Fig.14 Sound field comparison of two focusing methods

将实验检测到的声场与电场数据归一化处理并绘制如图15所示的归一化强度图，两条曲线分布基本重合，随着焦距的增加，声压和电场的强度均有所增强，并在40 mm 处达到峰值，焦距超过40 mm后，强度均有所下降。通过计算得到两条曲线的相关系数为0.913，据此可得，换能器声场和磁声耦合产生的感应电场分布具有较高的一致性。

图15 电场与声场分布Fig.15 Electric field and sound field distribution

为了便于进行比较，将水听器经前置放大器连接到数字示波器显示界面的波形(图16)记录下来，得到时域响应图，并将实验数据导入到Matlab设置采样频率进行频谱分析，图17(a)为超声信号频谱分析图，横轴为频率，纵轴表示能量大小；同理将导电凝胶两端的电流连接示波器显示波形并记录，利用Matlab 进行频谱分析，得到图17(b)的电流信号的频谱分析图。电流波形频谱分析要比超声波形的频谱分析效果要好，这是因为测量超声时噪声影响较大，使得采集到的信号频带较宽，但仍可得到大部分集中在1 MHz 左右，而测量电流信号时导电凝胶杂质较少，使得信号较为纯净，相对频带窄得多，超声基波频率为1 MHz，所测得的电流信号是频率为1 MHz 的交流信号。两种波形频谱分析得到的中心频率一致的结论，可证明磁场和声场耦合确实可以产生电流进行刺激。

图16 测量超声波形和测量电流波形Fig.16 Measuring ultrasonic wave shape and measuring current waveform

图17 超声波形频谱分析和电流波形频谱分析Fig.17 Ultrasonic shape spectrum analysis and current waveform spectrum analysis

此外，为进一步验证检测到的电信号确为声场与磁场共同作用下产生，将导电样本分别放置于静磁场和声场中，均未检测到电信号的产生，此实验联合声电场检测实验综合证明了经颅磁声电刺激可行性。

4 结论

通过建立基于真实颅骨和大脑的物理场模型，分别运用传统相控聚焦方法和时间反演法对大脑进行经颅磁声电刺激，并对结果进行分析。由于传统相控阵聚焦方法只考虑聚焦点和各阵元的相对位置，并未对真实经颅超声的传播过程进行计算，颅骨的形状、厚度、轮廓引起的相位失真、能量衰减和折射效应等未得到考量，运用时间反演法能够在非均匀介质中实现自适应聚焦。仿真结果证明基于时间反演的相位调控能够较好地解决因传播介质的非均匀性而导致的焦区偏移问题，提高了聚焦精度。

时间反演聚焦时各阵元的延迟信号可以通过计算机模拟来获得，从而避免在人体组织植入声源造成损伤。通过搭建实验平台，对两种聚焦方式的声场归一化处理，验证了时间反演法可以补偿焦点偏移，并通过对磁声耦合产生的电场与声场强度归一化处理，分析得出了它们之间存在较高的一致性。但本实验仅考虑了颅骨对超声的影响，后续需进一步考虑头皮、白质、灰质对超声聚焦的影响。基于时间反演的经颅磁声电刺激拥有高空间分辨率以及有效精准刺激的优点，有助推动精准神经调控的发展。