用于经颅超声刺激的超声换能器的实时导航方法及系统

1.本发明涉及神经调控和脑功能检测领域，特别涉及一种用于经颅超声刺激的超声换能器的实时导航方法及系统。


背景技术：

2.经颅超声刺激(tus)利用低频低强度的聚焦超声，跟现阶段其他常见的非侵入式神经调控技术诸如经颅直流电刺激(tdcs)、经颅磁刺激(tms)等相比具有空间分辨率高和聚焦深度深的特点。目前，已有一批研究展示了tus对健康受试者的大脑神经活动的即时性调控效应和持续性调控效应，能引起人脑皮层的可塑性变化，也有临床研究报道了tus在癫痫、脑卒中、意识障碍等神经或精神类脑疾病的应用。这些研究结果显示tus在人类脑功能研究和脑疾病临床治疗上具有很大的潜力。然而，要充分发挥tus技术的高空间分辨率的优点，需使得超声波穿过颅骨等脑组织后精确聚焦于目标治疗靶点，这在技术上具有挑战性，需要克服不同因素带来的各种困难。
3.在目前的tus实际应用中，使用者一般是通过手动确定超声换能器的粗略位置，或者通过适用于其他神经调控技术如tms的导航系统来引导放置超声换能器。由于超声波穿过颅骨等脑组织的过程中会发生超声传播方向的偏折，从而导致在颅内聚焦的焦点相对于在没有颅骨等脑组织情况下的焦点出现位置偏移和形状畸变。现有的适用于其他神经调控技术的导航系统未能考虑到声场经颅传播对焦点位置的影响，从而导致难以引导放置超声换能器到达合适的位置而使得经颅超声精准地聚焦到目标靶区。
4.为了解决上述问题，需要开发一种适用于经颅超声刺激的导航方法和系统，该方法和系统通过融入快速的经颅声场仿真计算结果到导航系统，从而实现对tus的高空间精度的导航。首先，因为大脑颅骨对于超声波传播方向的偏折使得无法通过焦点与超声换能器的相对位置来直接确定实际经颅声场的焦点的位置；同时，又由于在神经调控过程中，低强度聚焦超声在颅内的声压或者声强难以进行实时测量，故而无法通过实时测量的数据来反馈引导调整超声换能器位置和朝向。解决这一难题的一个方法是通过头颅的影像信息来检测出头皮、颅骨等脑组织的结构信息，然后基于这些脑组织结构信息及声学特性来进行超声传播仿真计算，以确定超声换能器相对于刺激靶点的放置位置和朝向。
5.在经过仿真得到超声换能器预期的入射朝向和位置后，需要通过一种手段将现实中的换能器引导放置到该预期位置。同时因为经颅聚焦超声是非侵入式的，无法直接观测脑内部的组织结构。结合以上的要求，一种可行的导航方案是采用无框架的神经导航技术。
6.针对神经调控技术的导航系统主要分为两种：基于磁共振引导的导航技术和基于光学的导航技术。基于磁共振的导航系统的精度高，但是对于其他手术设备的要求高，技术操作难度较高；而且，用于神经调控的低强度聚焦超声引起的脑组织温度变化很小，这微小的温度变化难以通过磁共振成像技术检测出来。基于光学的导航系统操作简单，成本低，对使用环境要求低，与tus技术有很高的适配性。
7.目前，现有技术中也提出了tus导航方案，但这些方案在声场仿真方案上没有考虑
到实际聚焦超声刺激其实是正向传播的，基于时间反演方法的技术方案得到的超声换能器在摆放位置发射的超声波在颅内的聚焦焦点依然可能与目标靶点存在一定程度的空间误差；且这些方案也未能提供可视化的声场分布信息来引导超声换能器的放置。
8.因此，本领域急需开发一种用于tus的超声换能器的实时导航方法及系统，以有效降低颅骨带来的超声焦点位置偏移误差，引导超声换能器的放置，从而让经颅超声的焦点精确地落在颅内预定刺激靶点上。


技术实现要素：

9.本技术的目的在于提供一种用于经颅超声刺激的超声换能器的实时导航方法及系统，能够有效降低颅骨带来的超声焦点位置偏移误差，引导超声换能器放置，从而让经颅超声的焦点精确地落在颅内预定靶点上。
10.本技术提供了一种用于经颅超声刺激的超声换能器的实时导航方法，包括以下步骤：
11.(a)建立超声穿过颅骨的声场计算模型，基于颅骨形态和刺激需求预先标注颅内预定刺激靶点，并预先确定超声换能器的初始位置和角度；
12.(b)基于超声换能器的初始位置和角度放置代表聚焦型的超声换能器的凹面声源，进行正向经颅超声声场仿真计算，进而得到所述超声换能器产生的经颅超声声场的仿真声学焦点(超声换能器仿真声学焦点)的位置；
13.(c)计算此时所述超声换能器与所述颅内预定刺激靶点的距离、以及所述超声换能器与经颅超声声场的仿真声学焦点的位置之间的距离，然后利用梯度下降法，计算得到中间轮次的超声换能器的位置和角度；
14.(d)基于步骤(c)中计算得到的中间轮次的超声换能器的位置和角度放置代表聚焦型的超声换能器的凹面声源，进行正向经颅超声声场仿真计算，从而得到此轮中的超声换能器产生的经颅超声声场的仿真声学焦点；
15.(e)重复迭代步骤(c)和(d)，当达到迭代停止条件时，停止迭代，进而得到用于导航的最佳超声换能器的位置和角度；其中，每一轮迭代计算使用的超声换能器的位置和角度是基于梯度下降法由上一轮迭代计算的结果得到的。
16.在另一优选例中，在获得一中间轮次的超声换能器的位置和角度后进行一次正向经颅超声声场仿真计算，从而得到在中间轮次的超声换能器的位置和角度下的此时的超声换能器产生的经颅超声声场的仿真声学实际焦点。
17.在另一优选例中，在使用梯度下降法进行下一中间轮次的超声换能器的位置和角度计算时，需要确定在上一中间轮次的超声换能器的位置和角度下进行正向超声声场仿真计算得到的声压声场在预定刺激靶点p0处的梯度。
18.在另一优选例中，考虑到经颅超声的仿真声场可能存在多个局部极小值点或多个局部极大值点，故先使用高斯平滑滤波器算法对仿真计算得到经颅超声声场进行平滑，以使得平滑后的经颅超声声场的半峰全宽区域内不存在局部极小值点，在最大声压值点外也不存在局部极大值点，然后使用梯度下降法进行下一中间轮次的超声换能器的位置和角度计算时，需要确定平滑后的超声仿真声场在预定刺激靶点p0处的梯度。
19.在另一优选例中，在步骤(c)中还包括：在利用梯度下降法计算得到中间轮次的超
声换能器的位置和角度之前，包括对超声声场仿真计算得到的此时的经颅超声声场进行平滑的步骤，从而使得平滑后的经颅超声声场的半峰全宽区域内不存在局部极小值点，在最大声压值点外也不存在局部极大值点。
20.在另一优选例中，所述对超声声场仿真计算得到的经颅超声声场进行平滑的步骤发生在步骤(c)中的计算此时所述超声换能器与所述颅内预定刺激靶点的距离以及所述超声换能器与经颅超声声场的仿真声学焦点的位置之间的距离之后。
21.在另一优选例中，通过使用高斯平滑滤波器算法对仿真计算得到经颅超声声场进行平滑，从而在使用梯度下降法进行下一中间轮次的超声换能器的位置和角度计算时，需要确定平滑后的经颅超声声场在预定刺激靶点p0处的梯度
22.在另一优选例中，在步骤(b)中，在凹面声源各位置设置大量的点声源，通过进行经颅声场仿真计算可以得到经颅超声声场的半峰全宽区域，并将该区域的几何中心点或峰值点作为所述超声换能器产生的经颅超声声场的仿真声学焦点。
23.在另一优选例中，在步骤(c)中，在第一轮迭代计算中，通过以下公式计算超声换能器与颅内预定刺激靶点的距离d
1,0
，以及超声换能器与超声换能器仿真声学焦点的位置之间的距离d
1,b
：
24.d
1,0
＝|p
0-p
1,a
|；
25.d
1,b
＝|p
1,b-p
1,a
|；
26.其中，p0为颅内预定刺激靶点的位置，p
1,a
为超声换能器的初始位置，p
1,b
为步骤(b)中超声声场仿真计算得出的超声换能器的经颅仿真声学焦点的位置。
27.在另一优选例中，在步骤(c)中利用梯度下降法，计算中间轮次的超声换能器的位置和角度包括：在第一轮迭代计算中，通过来计算第一中间轮次的超声换能器的角度和位置，其中，
28.第一中间轮次的超声换能器的角度按照以下计算公式确定：
[0029][0030][0031]
其中l为设定的迭代步长；
[0032]
第一中间轮次的超声换能器的位置p
2,a
按照以下计算公式确定：
[0033][0034]
在另一优选例中，所述迭代步长l根据操作者的经验进行设置，例如在在另一优选例中，所述迭代步长l根据操作者的经验进行设置，例如在之间。
[0035]
在另一优选例中，所述迭代停止条件为：当迭代轮数大于预设迭代轮数阈值n时，停止迭代。
[0036]
在另一优选例中，所述迭代停止条件为：当迭代轮数为n时，p
n,b
与p0之间的距离误差小于预设的导航空间误差阈值e时，或当两轮相邻的迭代所得到的导航空间误差值的下降量小于某一预设阈值时，停止迭代。
[0037]
在另一优选例中，步骤(c)和(d)中的超声换能器的角度和位置的具体迭代计算公
式如下：
[0038][0039][0040][0041]
其中m＝1,2,
…
为迭代次数，p
m,a
为第m轮经颅超声声场仿真计算中的超声换能器的位置，p
m,b
为第m轮经颅超声声场仿真计算中得到的超声换能器产生的经颅超声声场的仿真声学焦点的位置，是第m轮声场仿真计算得到的经颅超声声场在预定刺激靶点p0处的仿真声场梯度，为第m中间轮次的超声换能器的角度，p
m 1,a
为第m中间轮次的超声换能器的位置。
[0042]
在另一优选例中，是通过使用高斯平滑滤波器算法对第m轮声场仿真计算得到的经颅超声声场进行平滑后，平滑后的经颅超声声场在预定刺激靶点p0处的仿真声场梯度。
[0043]
在另一优选例中，所述经颅超声声场仿真计算的方法采用的是k空间伪谱法。
[0044]
在另一优选例中，所述超声换能器为单阵超声换能器或多阵元超声换能器组成的阵列。
[0045]
在另一优选例中，在步骤(a)之前还包括步骤(a0)：对受试者头颅的影像数据进行预处理，从而获得受试者头颅结构模型以及声场仿真介质数据，其中所述声场仿真介质数据用于建立超声经过颅骨的声场计算模型。
[0046]
在另一优选例中，所述方法还包括注册配准步骤，所述注册配准步骤通过基于受试者头颅结构模型，将在现实坐标系下头颅的最佳超声换能器的位置和角度转换成在图形界面的图像坐标系下头颅结构模型的对应位置和角度，从而能在图形界面的图像坐标系下引导现实坐标系下的超声换能器到达指定位置。
[0047]
在另一优选例中，所述注册配准步骤还包括在受试者头颅结构模型上选取人体面部特征点在图像坐标系中的对应点而形成配对点集，并使用刚体获取所述图像坐标系中选定的所述人体面部特征点在现实坐标系中的位置。
[0048]
在另一优选例中，所述人体面部特征点为3个或3个以上。
[0049]
在另一优选例中，所述体面部特征点选自下组：眼眶外端、耳朵的远端、耳垂根部、鼻尖、鼻根、或唇角。
[0050]
在另一优选例中，所述注册配准步骤通过奇异值分解法对特征点的进行配对计算，实现注册配准。
[0051]
在另一优选例中，假设所述现实坐标系是第一刚性系统，所述图像坐标系是第二刚性系统，所述现实坐标系和图像坐标系通过旋转和平移进行转换，从而得出现实坐标系和图像坐标系的转换关系，进而得到现实坐标系上的任意一个点在图像坐标系上的对应位置。
[0052]
在另一优选例中，所述注册配准步骤发生在声场仿真计算步骤之前或之后或同时进行；为避免受试者等待，注册配准步骤可在声场仿真计算步骤之后进行。
[0053]
本技术还提供了一种用于经颅超声刺激的超声换能器的实时导航系统，包括：
[0054]
图像预处理模块，所述图像预处理模块被配置为对受试者头颅的图像数据进行预处理，从而获得受试者头颅结构模型以及声场仿真介质数据；
[0055]
声场仿真计算模块，所述声场仿真计算模块被配置为基于声场仿真介质数据建立超声声场仿真计算模型，并通过梯度下降法调整换能器的位置和角度并进行声场计算，进而获得最佳超声换能器的朝向和位置；
[0056]
注册配准模块，所述注册配准模块被配置为基于受试者头颅结构模型，将在现实坐标系下头颅的最佳超声换能器的位置和角度和在图像坐标系下头颅结构模型中的对应位置和角度的互相转换。
[0057]
在另一优选例中，所述系统还包括空间追踪定位设备，所述空间追踪定位设备被配置为利用红外光学追踪定位仪，获取定位刚体在现实空间坐标系下的坐标数据。
[0058]
在另一优选例中，所述空间追踪定位设备包括设置第一刚体、第二刚体以及第三刚体，所述第一刚体在现实空间坐标系下指向头颅表面的某一特征点，用于注册配准时获取特征点在现实坐标系下的坐标数据，所述第二刚体绑定在超声换能器上，通过与超声换能器进行标定，用于获取超声换能器在现实坐标系的坐标数据，所述第三刚体绑定在受试者的头部，用于获取头颅在现实坐标系下的坐标数据，以追踪头颅的空间位置。
[0059]
在另一优选例中，所述系统还包括设备通讯模块，所述设备通讯模块被配置为实时获取第一刚体、第二刚体以及第三刚体在现实坐标系下的位置和角度数据，并将获得的数据传递给所述注册配准模块。
[0060]
在另一优选例中，所述系统还包括实时导航模块，所述实时导航模块被配置为实时显示超声换能器在图像坐标系内的具体位置和角度以及引导超声换能器至最佳位置和角度。
[0061]
在另一优选例中，所述图像预处理模块还被配置为可视化展示头颅的三维模型、并展示所述超声换能器和颅内预定刺激靶点以及预期(最佳)超声换能器位置彼此之间的相对位置和角度关系，并实时在头颅结构模型上叠加显示在预期(最佳)超声换能器位置和角度上的超声换能器所产生的经颅超声声场分布。
[0062]
应理解，在本发明范围内中，本发明的上述各技术特征和在下文(如实施例)中具体描述的各技术特征之间都可以互相组合，从而构成新的或优选的技术方案。限于篇幅，在此不再一一累述。
附图说明
[0063]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。应理解，下面描述中的附图仅仅是本发明的一些实施实例，本领域的普通技术人员，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的实施实例。
[0064]
图1为根据本技术的一个实施例的用于经颅超声刺激的超声换能器的实时导航系统的结构框图；
[0065]
图2为根据本技术的一个实施例的经过注册配准，通过图形界面引导超声换能器达到指定位置的示意图；
[0066]
图3为根据本技术的一个实施例的基于梯度下降的原理寻找超声换能器最佳入射位置和朝向的迭代算法示意图(其中a：超声换能器的初始位置和角度；b和c：基于梯度下降的迭代算法的中间过程；d：迭代算法的终止条件)；
[0067]
图4为根据本技术的一个实施例的用于经颅超声刺激的超声换能器的实时导航方法的流程图。
[0068]
各附图中，各标识如下：
[0069]
1-超声换能器
[0070]
2-空间追踪定位设备
[0071]
3-第一刚体(指针刚体)
[0072]
4-第二刚体(换能器刚体)
[0073]
5-第三刚体(头动刚体)
具体实施方式
[0074]
本发明人通过广泛而深入的研究，首次开发了一种用于经颅超声刺激的超声换能器的实时导航方法及系统，本技术的实时导航系统和方法可以通过经颅超声声场仿真寻找单阵元超声换能器或多阵元超声换能器组成的阵列的最佳放置位置及朝向，并通过注册配准将在现实坐标系下头颅的最佳超声换能器的位置和角度和在图像坐标系下头颅结构模型中的对应位置和角度的互相转换，从而能在图形界面的图像坐标系下引导现实坐标系下的超声换能器到达指定位置，进而实现超声换能器或换能器阵列发出的超声波经过颅骨等脑组织后精确地聚焦于颅内预先规划的刺激靶点。
[0075]
术语
[0076]
如本文所使用的，“超声换能器的朝向”和“超声换能器的角度”可互换使用；
[0077]
如本文所使用的，“预设刺激靶点”、“预定刺激靶点”、“颅内预定刺激靶点”以及“目标靶点”可互换使用；
[0078]
如本文所使用的，“经颅超声声场的仿真声学焦点”即为“超声换能器产生的经颅超声声场的仿真声学焦点”，可以称为“超声换能器仿真声学焦点”或者“超声换能器的经颅仿真声学焦点”。
[0079]
经颅超声刺激
[0080]
经颅超声刺激(transcranial ultrasound stimulation,tus)是近年快速发展的一种非侵入式功能性神经调控技术，通过透过头皮和颅骨等脑组织的低频低强度脉冲超声波来调控(增强或抑制)大脑神经活动。
[0081]
神经导航技术
[0082]
神经导航技术又称图像引导神经外科手术(image-guided neurosurgery，igns)，是基于计算机断层扫描成像(computed tomography，ct)、磁共振成像(magnetic resonance image，mri)等术前影像数据建立图像引导空间，借助光学或者其他方式追踪设备实时跟踪显示手术仪器相对于脑组织和病变部位的位置关系，从而指导医生进行手术操作的技术。
[0083]
超声换能器：
[0084]
超声换能器根据压电效应，可将电能转化为机械(声)能，也可将机械(声)能转化
为电能，可以负责产生超声波，也可以接收并检测超声波。一般由五个主要部件组成：具有压电特性的晶体/陶瓷元件、元件表面的正电极和地电极、阻尼(背衬)块、匹配层和外壳
[0085]
配准：
[0086]
图像配准是将不同来源的数据转换到一个坐标系的过程。数据可能来自于不同传感器、时间、深度、视角。常常用于计算机视觉，医学成像，军事自动化目标识别等等；配准是比较和整合不同测量方式获得的数据的必要过程。
[0087]
术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。本专利的申请文件中，如果提到根据某要素执行某行为，则是指至少根据该要素执行该行为的意思，其中包括了两种情况：仅根据该要素执行该行为、和根据该要素和其它要素执行该行为。多个、多次、多种等表达包括2个、2次、2种以及2个以上、2次以上、2种以上。需要说明的是，在本专利的申请文件中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。
[0088]
在本发明中，所有方向性指示(诸如上、下、左、右、前、后等)仅用于解释在某一特定姿态(如附图所示)下各部件之间的相对位置关系、运动情况等，如果该特定姿态发生改变时，则该方向性指示也相应地随之改变。
[0089]
本发明的主要优点：
[0090]
(a)本技术的用于经颅超声刺激的超声换能器的实时导航方法及系统通过基于梯度下降原理，进行经颅超声声场迭代仿真计算寻找出经颅超声刺激中的超声换能器最佳放置位置和朝向，能够有效降低颅骨带来的经颅超声焦点位置的偏移误差，引导超声换能器的放置，从而让经颅超声的焦点精确地落在颅内预设的刺激靶点上，并提供可视化的声场分布信息来引导超声换能器的放置；
[0091]
(b)在本技术的用于经颅超声刺激的超声换能器的实时导航方法及系统中，根据使用需要，可采用多种思路在头皮处设置超声换能器的初始位置和朝向，然后经过正向经颅声场仿真得到超声换能器此时的实际焦点以及声场数据，并结合声场焦点和目标靶点之间的关系，基于梯度下降的原理推算出下一步的超声换能器的放置位置和朝向，然后再进行下一轮的正向经颅声场仿真计算，不断迭代直到满足终止条件得到所求(最佳)超声换能器的位置和朝向；
[0092]
(c)现有导航技术没有考虑到颅骨等脑组织对经颅超声传播的影响，或是基于其他方法确定超声换能器的位置和朝向来进行导航，这都存在一定程度的空间误差，且未能提供正向经颅声场分布信息引导的导航；本发明进行正向的经颅超声声场计算，通过基于梯度下降的原理进行正向迭代声场仿真计算，不断修正超声换能器位置和朝向，最终使得超声换能器发出的超声波经过颅骨传播后的声学焦点与预定刺激靶点的空间位置误差更小，从而实现高空间精度的经颅超声刺激导航；
[0093]
(d)本技术的用于经颅超声刺激的超声换能器的实时导航方法及系统分别通过软件交互和光学追踪设备(空间追踪定位设备)获取头部多个解剖特征点在图像坐标系中和
现实坐标系中的位置，得到处于两个坐标系中的配对点集，通过使用奇异值分解计算得到两个点集之间的旋转矩阵和向量，得到两个坐标系之间的转换关系，用于软件图形界面上显示超声换能器在图像坐标系的位置；
[0094]
(e)本技术中的导航系统采用奇异值分解求解刚体旋转平移矩阵来进行注册配准，相较于现有技术，本技术提出的注册配准方法通过增加参与注册配准的配对点数目来降低噪声影响和获取配对点坐标时的操作误差，实现高精度的注册配准，从而得到更高精度的经颅超声刺激导航；
[0095]
(f)本技术设计得到tus导航方法将有力地推动tus的发展和应用，为有关基于tus的脑科学研究以及脑疾病临床治疗提供一种新的工具，充分发挥tus的高空间精度优点。
[0096]
在以下的叙述中，为了使读者更好地理解本技术而提出了许多技术细节。但是，本领域的普通技术人员可以理解，即使没有这些技术细节和基于以下各实施方式的种种变化和修改，也可以实现本技术所要求保护的技术方案。
[0097]
用于经颅超声刺激的超声换能器的实时导航方法
[0098]
在本技术中提供了一种用于经颅超声刺激的超声换能器的实时导航方法，该方法包括以下步骤：
[0099]
(a0)对受试者头颅的图像数据进行预处理，从而获得受试者头颅结构模型以及声场仿真介质数据；
[0100]
(a)基于声场仿真介质数据建立超声穿过颅骨的声场计算模型，基于颅内预定的刺激靶点确定超声换能器的初始位置和角度；
[0101]
(b)基于超声换能器的初始位置和角度放置代表聚焦型的超声换能器的凹面声源，进行正向经颅超声声场仿真计算，进而得到所述超声换能器产生的经颅超声声场的仿真声学焦点的位置；
[0102]
(c)计算此时所述超声换能器与所述颅内预定刺激靶点的距离以及所述超声换能器与经颅超声声场的仿真声学焦点的位置之间的距离，然后利用梯度下降法，计算中间轮次的超声换能器的位置和角度；
[0103]
(d)基于步骤(c)中计算得到的中间轮次的超声换能器的位置和角度，进行正向经颅超声声场仿真计算，从而得到新的超声换能器产生的经颅超声声场的仿真声学焦点；
[0104]
(e)重复迭代步骤(c)和(d)，当达到迭代停止条件时，停止迭代，进而得到用于导航的最佳超声换能器的位置和角度；其中，每一轮迭代计算使用的超声换能器的位置和角度是基于梯度下降法由上一轮迭代计算的结果得到。
[0105]
其中，在步骤(b)中，在代表超声换能器的凹面声源各位置设置大量的点声源，通过进行声场仿真计算之后可以得到经颅超声声场的半峰全宽区域，并将该区域的几何中心点作为所述超声换能器产生的经颅超声声场的仿真声学焦点；也可将半峰全宽区域的峰值点作为所述超声换能器产生的经颅超声声场的仿真声学焦点。
[0106]
另外，在步骤(c)中，考虑到经颅超声的仿真声场可能存在多个局部极小值点或多个局部极大值点，故先使用高斯平滑滤波器算法对仿真计算得到经颅超声声场进行平滑，以使得平滑后的经颅超声声场的半峰全宽区域内不存在局部极小值点，在最大声压值点外也不存在局部极大值点，然后使用梯度下降法进行下一中间轮次的超声换能器的位置和角度计算时，需要确定平滑后的超声仿真声场在预定刺激靶点p0处的梯度。
[0107]
换句话说，在步骤(c)中还包括：在利用梯度下降法计算得到中间轮次的超声换能器的位置和角度之前，且在步骤(c)中的计算此时所述超声换能器与所述颅内预定刺激靶点的距离以及所述超声换能器与超声换能器仿真声学焦点的位置之间的距离之后，包括对超声声场仿真计算得到的此时的经颅超声声场进行平滑的步骤，从而使得平滑后的经颅超声声场的半峰全宽区域内不存在局部极小值点，在最大声压值点外也不存在局部极大值点。
[0108]
优选地，通过使用高斯平滑滤波器算法对仿真计算得到经颅超声声场进行平滑，从而在使用梯度下降法进行下一中间轮次的超声换能器的位置和角度计算时，需要确定平滑后的经颅超声声场在预定刺激靶点p0处的梯度需要注意的是，在每次迭代中，在利用梯度下降法，计算得到中间轮次的超声换能器的位置和角度之前，需要对此时超声声场仿真计算得到的此时的经颅超声声场进行平滑，以使得平滑后的经颅超声声场的半峰全宽区域内不存在局部极小值点，在最大声压值点外也不存在局部极大值点，然后使用梯度下降法进行下一中间轮次的超声换能器的位置和角度计算时，需要确定平滑后的超声仿真声场在预定刺激靶点p0处的梯度。
[0109]
该方法还包括注册配准步骤，通过基于受试者头颅结构模型，将在现实坐标系下头颅的最佳超声换能器的位置和角度转换成在图像坐标系下头颅结构模型的对应位置和角度，从而实现能在图形界面的图像坐标系下引导现实坐标系下的超声换能器到达指定位置。
[0110]
注册配准步骤还包括在受试者头颅结构模型上选取人体面部特征点在图像坐标系中的对应点而形成配对点集，并使用刚体获取所述图像坐标系中选定的所述人体面部特征点在现实坐标系中的位置。
[0111]
注册配准步骤通过奇异值分解法进行特征点的配对。在一实施例中，假设所述现实坐标系是第一刚性系统，所述图像坐标系是第二刚性系统，所述现实坐标系和图像坐标系通过旋转和平移进行转换，从而得出现实坐标系和图像坐标系的转换关系，进而得到现实坐标系上的任意一个点在图像坐标系上的对应位置。在一实施例中，注册配准步骤发生在声场仿真计算步骤之前或之后或同时进行；为避免受试者等待，注册配准步骤可在声场仿真计算步骤之后进行。
[0112]
其中注册配准步骤通过注册配准模块、空间追踪定位设备以及设备通讯模块实现。
[0113]
用于经颅超声刺激的超声换能器的实时导航系统
[0114]
如图1所示，本技术还提供了一种用于经颅超声刺激的超声换能器的实时导航系统，包括：图像预处理模块、声场仿真计算模块、空间追踪定位设备、设备通信模块、注册配准模块和实时导航模块。
[0115]
其中，实时导航模块提供了实时显示超声换能器在图像坐标系内的具体位置以及引导超声换能器至目标位置(预定刺激靶点)的辅助指示。设备通讯模块和声场仿真计算模块为导航系统提供了获取超声换能器现实坐标系下的坐标和图像坐标系的目标位置，注册配准模块为导航系统提供两个坐标系的转换关系，图像数据预处理模块则实现头颅结构模型的三维可视化展示的功能，更为直观地展示超声换能器和颅骨内的预设靶点以及目标换能器位置的相对位置关系，并在头颅结构模型上叠加显示在最终确定的超声换能器位置和
朝向上的超声换能器所产生的经颅超声声场的分布。
[0116]
图像预处理模块
[0117]
图像预处理模块首先将从受试者采集到的ct或者mri图像数据进行预处理，然后从图像数据构建注册配准模块中需要的头颅结构模型，以及声场仿真计算模块中需要的声场介质数据。mri对于大脑的灰质和白质等软组织成像较为清晰，而一般结构像mri对于骨组织的成像信号较弱；由于骨组织吸收射线能力较强，ct对骨组织的成像与其他组织相比有明显的对比度，是获取声场仿真介质数据的主要获取方式。本发明将主要阐述对于ct图像数据的预处理，具体步骤如下所述：
[0118]
首先，通过对受试者头部进行ct扫描得到大小为n
x
×
ny×
nz三维图像，像素间隔为d
x
×dy
×dz
，同时三维图像数据在图2所示的图像坐标系的原点坐标为o
x
×
oy×
oz。图像索引为(i,j,k)的体素通过以下公式转化为在图像坐标系的坐标(x,y,z)
[0119][0120]
其中i＝1，2，...，n
x
，j＝1，2，...，ny，k＝1，2，...，nz。然后，设定颅骨的ct值阈值t，对图像像素数据进行遍历，若像素强度大于t则保留，否则将去除，最后得到只保留颅骨信息的三维数据。由于颅骨不同位置上孔隙度不同以及密度不同，不是均匀介质，本发明利用下述ct值与密度以及声速的经验公式(2)、(3)和(4)，并通过方程(1)映射到图像坐标系中，得到头部声速以及密度在图像不同位置上的分布。ct值与密度以及声速的经验公式为
[0121][0122]c(i，j，k)
＝c
水
ψ
(i，j，k)
c
骨
(1-ψ
(i，j，k)
)，
ꢀꢀꢀꢀꢀ
(3)
[0123]
ρ
(i，j，k)
＝ρ
水
ψ
(i，j，k)
ρ
骨
(1-ψ
(i，j，k)
)，
ꢀꢀꢀꢀꢀꢀꢀ
(4)
[0124]
其中h是图像范围内索引为(i，j，k)的ct值，ψ是孔隙度，c
水
和c
骨
是超声波在水中的声速和在该位置颅骨中的声速，ρ
水
和ρ
骨
是水的密度和该位置的颅骨的密度。
[0125]
对于图像数据处理：在本发明实施例中，我们主要是以头颅的ct影像数据作为输入，来提取颅骨等组织信息作为仿真声场介质参数输入到仿真程序中，使用方程(2)-(4)进行ct值到介质密度和声速的转换。在其他实施例中，另外也可以由头部的磁共振成像(mri)数据来合成头部的ct图像，然后由合成的ct图像来提取颅骨等组织信息作为仿真声场介质参数输入到仿真程序中。也有方法直接使用mri图像作为输入，从mri图像中分割出颅骨并取颅骨密度和声速的代表值。
[0126]
声场仿真计算模块
[0127]
声场仿真计算模块用于建立超声声场计算模型，确定超声换能器的放置位置以及朝向。具体的步骤如图4所示。步骤(1)：确定超声换能器的初始位置和朝向。本发明先通过图像预处理模块处理受试者原始图像后得到声场参数模型，然后在图像空间中标定预先规划的刺激靶点(预定刺激靶点)，并在图像空间上标定超声换能器出射面轴心点和头皮接触的位置，作为超声换能器的初始位置。该初始位置可以根据需要通过多种方法确定和调整，如可由操作人员根据经验确定，或通过计算各位置头皮与靶点的距离并选定距离最近的头皮位置作为超声换能器的初始位置，或通过计算找到头皮上的某一点作为超声换能器的初
始位置，该点的颅骨切平面与该点和靶点之间的连线垂直。在初步设定超声换能器的初始位置点后，进一步计算该初始位置点在颅骨表面的向内的法向量，如该初始位置点与靶点之间的连线与该法向量夹角过大，则该初步设定的位置点不宜作为超声换能器的初始位置点，应进一步调整初始位置点，以使得该夹角较小为宜。这是因为超声波沿颅骨的法向量入射，颅骨对超声波传播方向的偏折小；超声波入射方向偏离颅骨的法向量越大，则颅骨引起超声传播方向的偏折越大。经调整并确定超声换能器的初始位置后，把该位置点向靶点的连线方向作为超声换能器的中轴线初始朝向。
[0128]
步骤(2)：进行声场计算，基于本发明提出的梯度下降法调整超声换能器的朝向和位置。根据所使用的超声换能器的具体尺寸和结构设计，设置用于仿真计算的凹面声源，凹面声源的位置和朝向按步骤(1)中确定的超声换能器的初始位置和朝向来设定。需要说明的是该凹面声源是由多个点声源组成，具有相同的超声声源发射参数，这些点声源通过正向经颅传播后得到了整个凹面声源的仿真声学焦点，该焦点与目标靶点进行空间位置比较后开始利用本发明提出的梯度下降法调整超声换能器的朝向和位置，以减少代表超声换能器的凹面声源的聚焦焦点与预先设定的刺激靶点之间的空间误差。本发明每一轮经颅声场计算后均采用梯度下降法调整超声换能器的朝向和位置，然后按照调整后的朝向和位置进行下一轮的经颅声场计算，并按新的经颅声场计算结果采用梯度下降法进一步调整超声换能器的朝向。如此循环进行多轮经颅声场计算和朝向及位置调整，直到某一轮朝向调整后经颅声场计算出来的焦点和预先设定的靶点之间的空间距离小于使用者所设定的某一可接受的阈值或者迭代轮次达到使用者设定的迭代次数阈值，则停止迭代计算，最后一次调整的超声换能器的位置和朝向就是本发明所确定超声换能器的位置和朝向。使用该方法计算出来超声换能器的位置和头皮之间可能存在一定间距，则需在超声换能器和头皮之间填充超声耦合介质(如除气水囊、特制的凝胶垫片等)。
[0129]
本发明提出的梯度下降法的原理是向经颅超声仿真声场分布上当前点的对应梯度或者近似梯度的反方向进行迭代搜索，到达函数的局部极小值点。具体来说，本发明中要优化的目标函数是仿真得到的经颅声场聚焦焦点到目标靶点的距离，利用经颅超声仿真声场上靶点位置的梯度方向作为超声换能器的调整方向，并且给定步长更新超声换能器的位置以及朝向。
[0130]
超声声场仿真计算方法
[0131]
本发明中所述的超声声场仿真计算，为简化仿真步骤以及计算量和减少仿真计算对计算机内存的要求，并且不影响超声仿真结果精度的要求下，可假设超声在整个声场范围内是线性传播的，并且由于经颅聚焦超声的频率较低，可不考虑声场中的非线性衰减因子。同时为了降低对原始图像数据的空间分辨率要求以及提高仿真计算速度，本发明使用伪谱k空间(k-space)法进行仿真计算，该方法的计算步骤如下：
[0132]
对于均匀、无损耗的流体介质，线性声学波动方程可以通过以下三个方程推出
[0133][0134]
[0135][0136]
其中u是声学质点速度，p是声压，ρ是声学密度，c0是声速，ρ0是净场密度。该等式假设是在热粘性流体中，不存在弛豫机制。将以上的公式进行结合得到二阶波动方程
[0137][0138]
但为了计算的便利以及便于后续消除误差，本发明的声场仿真都将通过上述的一阶方程(5)、(6)和(7)进行计算，而非使用更为熟悉的二阶波动方程(8)。
[0139]
如果假设存在声波的吸收以及是非均质的情况，则需要在进一步在方程(6)和(7)上加入修正项，得到
[0140][0141][0142]
其中d是质点位移，l是计算声吸收和散射的线性整微分算子。
[0143]
以上公式描述了声场如何传播，可以通过添加初始条件来描述介质中如何产生声波
[0144][0145][0146]
其中sf是初始外部作用力源，sm则是初始质量输入源。
[0147]
本发明引入了伪谱k空间法求解偏微分方程。通过这种方法，波动方程中的偏微分可以通过傅里叶方程算出，而不是通过传统的有限差分近似来求解，以避免产生比较大的线性相位误差。同时，伪谱k空间法可以被视作无穷阶精确有限差分法的极限情况，理论上只需要每波长2个节点就能得到精确平滑变化的系数。对于一般的二阶波动方程，其具体形式为
[0148][0149]
其中是标量波数。
[0150]
对于一般三维离散的k空间法，可由以下方程进行计算。
[0151][0152]
[0153][0154][0155][0156]
其中，方程(14)和(16)表示基于伪谱k空间法的梯度计算公式，方程(15)、(17)和(18)则是对应连续变量情况下的公式(11)、(12)和(10)的k空间修正后的一阶精确前向差分；ξ＝x,y,z是三维线性空间的坐标轴方向；n和n 1是对应变量的当前时刻和下一时刻；和分别是傅里叶变换和傅里叶反变换，可以通过快速傅里叶变换求解；δt是时间间隔；是k空间运算符，c
ref
是背景声速；k
ξ
代表在ξ方向上的波数，离散形式如下
[0157][0158]
其中，n
ξ
是在ξ方向上的网格数；ld是幂律吸收的离散形式，具体公式为
[0159][0160]
其中分别代表吸收和散射比例系数，α0是幂律前置系数，y是幂律指数。
[0161]
同时整个声场的声强可以结合以下公式算出
[0162][0163]
为了能让离散仿真准确，需要保持收敛条件其中cfl数是
[0164]
最后，伪谱k空间法再引入完美匹配层(perfect matching layer,pml)来吸收仿真范围边界外声波，声波传播至边界处将会急剧衰减并且不会反射声波，防止傅里叶空间变换在边界处产生的震荡效应影响仿真的结果。伪谱k空间法通过迭代计算声场的声压梯度来得到下一时刻的声速分布，又通过计算声速的梯度来得到下一时刻声压的分布，最终完成整个仿真过程。
[0165]
根据以上的仿真公式，在给定超声换能器的结构信息和声传播介质信息后，就可以进行超声换能器凹面声源的正向传播仿真计算。
[0166]
需要说明的是，相比于其他计算方法，k空间伪谱法在计算精度和计算速度上面有
优势。
[0167]
基于梯度下降法调整超声换能器的朝向(角度)和位置
[0168]
在以上声场计算的基础上，本发明提出的用于确定超声换能器的朝向和位置的梯度下降法的步骤流程如图3和图4所示。具体的仿真步骤如下：
[0169]
步骤(1)：预设的经颅超声刺激靶点p0的具体位置可表示为(x0,y0,z0)，如图3(a)所示。首先，将在图像预处理模块中已将原始的图像数据转换成的声场仿真的介质信息，例如介质的声速、介质的密度、介质的形状和边界等，作为仿真参数输入。按上文所述的方法确定超声换能器的初始位置p
1,a
及超声换能器的中轴线初始朝向并在该初始位置上计算颅骨表面上的法向量判断初始放置位置和朝向的合理性。其中，法向量的计算可采用如下方法：对该位置点进行k近邻搜索，对邻域点构建协方差矩阵，求解协方差矩阵的特征向量及特征值，其中最小特征值对应的特征向量就是法向量。
[0170]
步骤(2)：为了模拟超声换能器的正向传播，需要在p
1,a
设置如图3所示的超声换能器的凹面声源，其中凹面曲率半径，换能器口径为，同时凹面声源的位置和朝向已通过步骤(1)确定。在凹面声源表面各位置设置大量的点声源，声场中设定外部声源的输入，在点声源处声压随着时间变化的方程为
[0171]
p
source
＝a
·
sin(2πft),
ꢀꢀ
(22)
[0172]
其中a是超声的声强振幅，f是超声的振荡频率，t是超声波发射的时刻。
[0173]
之后通过进行正向声场仿真计算之后可以得到放置于当前位置的超声换能器的经颅声压分布，凹面声源经过颅骨传播后的实际焦点处声场为长轴在声场传播方向的一长椭球形，为了后续计算方便，取半高全宽区的中心作为聚焦焦点p
1,b
。
[0174]
对于凹面单阵元换能器，其在自由声场中的焦距为一固定值，换而言之，如果实际焦点到凹面声源的距离和目标靶点到凹面声源的距离存在较大差距时，仅调整超声换能器的朝向无法使得目标靶点与实际焦点重合，还需要调整超声换能器到头皮的距离。首先计算此时目标靶点和超声换能器放置位置的距离d
1，0
＝|p
0-p
1，a
|以及经颅仿真声场焦点和超声换能器的距离d
1，b
＝|p
1，b-p
1，a
|。
[0175]
步骤(3)：确定正向经颅声场仿真计算出的声压声场在目标靶点p0处的梯度为本发明采用梯度下降的策略来调整超声换能器的朝向，使得超声换能器的朝向调整依据梯度的方向来确定，如图3(b)所示。记向量的方向来确定，如图3(b)所示。记向量本发明按经过梯度下降的策略确定的下一轮声场仿真计算的目标焦点位置为其中l为设定的迭代步长。则下一轮经颅声场仿真计算中，超声换能器的朝向按向量来设定。
[0176]
步骤(4)：利用步骤(2)计算得到的距离d
1，0
和d
1，b
来更新下一轮的经颅声场仿真计算中超声换能器的位置超声换能器的朝向和位置更新后，新一轮的正向经颅声场仿真计算示意图如图3(c)所示，可计算得到此时的经颅声场焦点为p
2，b
。
[0177]
本发明然后按照步骤(2)到步骤(4)重复迭代进行声场计算，具体迭代计算的方程如下
[0178][0179][0180][0181][0182]
其中m＝1，2，...，为迭代轮数，p
m，a
为第m轮经颅声场仿真计算中的超声换能器的位置，p
m，b
为第m轮经颅声场仿真计算中得到的经颅超声声场的焦点位置，p
m，c
为由第m轮经颅声场仿真计算中得到的声场按梯度下降的策略确定的第m 1轮声场轮仿真计算的目标焦点位置，是第m轮声场仿真计算得到的目标靶点处的声场梯度。本发明按上述迭代方程不断迭代调整超声换能器的朝向和位置，进行正向传播仿真，得到下一轮经颅超声换能器的声场分布，直至达到迭代停止条件时停止迭代计算。此时得到超声换能器的朝向和位置用于导航系统引导超声换能器放置的目标朝向和位置。迭代停止条件可根据需要来采用不同的方法设置，如当迭代轮数大于预设的迭代轮数阈值n时，或者当迭代轮数为n时，p
n，b
与p0之间的距离误差小于预设的导航空间误差阈值e时，或当两次相邻的迭代所得到的导航空间误差值的下降量小于某一预设阈值时，如图3(d)所示。
[0183]
优选地，考虑到经颅超声的仿真声场可能存在多个局部极小值点或多个局部极大值点，是通过使用高斯平滑滤波器算法对第m轮声场仿真计算得到的经颅超声声场进行平滑后，平滑后的经颅超声声场在预定刺激靶点p0处的仿真声场梯度，其中平滑后的经颅超声声场的半峰全宽区域内不存在局部极小值点，在最大声压值点外也不存在局部极大值点。
[0184]
空间追踪定位设备
[0185]
空间追踪定位设备2是一种通过跟踪附着在目标上的主动或被动红外标记的位置来实时确定目标位置的装置，而红外反射标记的位置是用一个空间追踪定位设备2(例如摄像系统)确定的(图2)。本实施例中的空间追踪定位设备使用的是加拿大ndi公司的polaris vicra红外光学定位仪。红外标记刚体如图2中的第一刚体3以及图2中第二刚体4、第三刚体5所示。第一刚体3通过旋转校准标定来确定针尖相对于刚体主体之间的相对距离，从而得到在现实坐标系下的位置，具体操作是：将第一刚体3的针尖插入在一圆孔凹槽内，然后让第一刚体3绕着圆孔旋转，同时不断获取第一刚体3工具坐标系相对于现实坐标系的刚体变换矩阵。由于针尖在现实坐标系和工具坐标系中位置都没有变化，但工具坐标系相对于现实坐标系的位置则会因为不断旋转第一刚体3而不断变化，具体是分布在现实坐标系下以针尖为球心、半径为工具坐标系原点到针尖的距离的球体表面上，根据以上关系可以得到相关方程，求解得到针尖在第一刚体3工具坐标系的位置，从而确定针尖在工具坐标系的坐标。
[0186]
本实施例同样可以设定超声换能器1在第二刚体4工具坐标系下的坐标，来确定超声换能器的位置和朝向。
[0187]
由上可知，空间追踪定位设备可以包括第一刚体3、第二刚体4以及第三刚体5，第一刚体3的针尖在现实坐标系下指向头颅表面的某一特征点，用于注册配准时获取特征点在现实坐标系下的坐标数据，第二刚体4绑定在超声换能器1上，通过与超声换能器1进行标
定，用于获取超声换能器在现实坐标系的坐标数据，第三刚体5绑定在受试者的头部，用于获取头颅在现实坐标系下的坐标数据，以追踪头颅的空间位置。
[0188]
设备通讯模块
[0189]
设备通讯模块则是利用空间追踪定位设备与计算机的连接例如通用串口、有线网线或者无线网络进行通信以及设备商提供的软件开发包的api接口来实时获取刚体3、4和5在现实坐标系下的位置和朝向，并将获得的数据传递给注册配准模块。
[0190]
注册配准模块以及注册配准方法步骤
[0191]
注册配准模块提供了一个选择图像坐标系和现实坐标系中对应点的交互图像界面以及配对点之间的配准算法。本实施例使用的是解剖学特征配对点注册配准的方式。具体来说，使用者需要通过鼠标在包含由图像预处理模块重建的三维头颅模型的交互界面上，选取人体面部特征如眼眶外端、耳朵的远端、鼻尖、鼻根、唇角等在图像坐标系中位置确定的点形成配对点集{di}，然后使用第一刚体的针尖点选图像坐标系中选定的各特征点在现实坐标系中对应的位置，并通过设备通讯模块与空间定位追踪设备通信，获取得到配对点集{si}的在现实坐标系中的位置。这样的配对点至少选择四对，且不处于一个平面上，而在本实施案例中则需要选择五对以上为佳。分别获取配对点集在图像坐标系和现实坐标系中的位置之后，实施案例使用奇异值分解法进行配对点的配准。具体来说，假设两个坐标系是刚性系统，两个坐标系可以通过旋转和平移进行互相转换，配对点集在在现实坐标系中的坐标{si}和在图像坐标系中的坐标{di}存在以下转换关系
[0192]di
＝rsi d，
ꢀꢀꢀ
(27)
[0193]
其中r是旋转矩阵，d是位移向量。r和d可以根据以下步骤解出：
[0194]
步骤(1)：首先分别求两个点云的质心
[0195][0196][0197]
步骤(2)：求各点相对于质心的位移向量，即d
′i＝d
i-d，s
′i＝s
i-s。
[0198]
步骤(3)：利用质心位移向量，计算h矩阵
[0199][0200]
步骤(4)：对h矩阵进行svd分解
[0201]
h＝uλv
t
.
ꢀꢀꢀ
(31)
[0202]
步骤(5)：基于矩阵u和v，计算旋转矩阵r
[0203]
r＝vu
t
.
ꢀꢀꢀ
(32)
[0204]
步骤(6)：计算det(r)，如果det(r)＝1，结果保留，否则应重新进行配准步骤；当两点集{di}和{si}不是分布在一个平面上且误差较小时，不会出现der(r)＜0的情况。
[0205]
步骤(7)：最后计算两个点集间的位移向量d
[0206]
d＝d-rs.
ꢀꢀꢀ
(33)
[0207]
此时就得出了两个坐标系之间的转换关系。利用这样的转换关系就可以得到现实坐标系上的任意一个点在图像坐标系上的对应位置，从而可以在实时导航模块的基于图像坐标系的图形界面中显示第二刚体的坐标，指示超声换能器位置，然后引导至预期超声换能器位置。
[0208]
实时导航模块
[0209]
实时导航模块提供了实时显示超声换能器在图像坐标系内的具体位置以及引导超声换能器至基于梯度下降法迭代计算得到的预期位置和朝向的辅助指示。
[0210]
实时导航模块通过设备通讯模块获取到空间追踪定位设备中第二刚体以及第三刚体在现实坐标系的坐标，进一步使用注册配准模块得到的刚体变换矩阵，得到上述刚体在图像坐标系下的坐标，从而指示了超声换能器以及头颅在图像坐标系的坐标。通过变换实际超声换能器的位置和角度使得其在图像坐标系下与基于梯度下降法迭代计算得到的预期(最佳)超声换能器的位置和角度重合。
[0211]
在本技术提及的所有文献都被认为是整体性地包括在本技术的公开内容中，以便在必要时可以作为修改的依据。此外应理解，在阅读了本技术的上述公开内容之后，本领域技术人员可以对本技术作各种改动或修改，这些等价形式同样落于本技术所要求保护的范围。