本发明属于超声检测及超声成像技术领域，具体涉及超声二维平面换能器阵列，以及经颅低频超分辨率容积三维成像方法。
背景技术：
目前应用超声对颅骨内的监测手段主要是经颅多普勒技术，即通过对回波信号的频移信息反推超声源与散射或反射目标间的相对运动，用来检测颅内脑底主要动脉的血流动力学及血流生理参数。然而经颅多普勒技术只能显示出血管血流的二维信息。
与传统经颅二维成像技术相比，经颅三维成像技术有着诸多优势：首先，三维成像更加直观，能直接显示出颅内结构的三维信息；其次，对三维图像进行断层分层分析，可从传统成像方法无法实现的角度进行观察，方便临床全面了解检测部位的解剖结构和疾病情况；再次，三维成像技术还可以提供更多的医学信息，例如组织表面积、组织容积，为临床诊断治疗提供更加准确可靠的信息。更重要的，脑血管的分布成像及相关疾病诊断需要的是颅内的三维信息，因此，经颅三维成像技术是经颅成像发展的必然趋势。
三维超声成像也称为超声体成像，根据超声三维体数据获取方式的不同，可以分为两种：一种是静态成像，即通过机械控制探头扫描得到空间截面的二维信号，再将二维信号进行空间重建得到三维图像，此方法有着较高的分辨率，但一定程度上丢失了时间信息，且需要精准的定位系统；第二种是动态成像，即容积成像技术，利用面阵容积探头发射空间超声波束直接获得三维超声体数据，不仅具有二维彩超换能器全部功能，还具有立体成像、图像切割、图像旋转及高平面图像分析功能。容积成像技术的优点在于图像获取简单快速，但缺点在于由于现有超声设备通道数的限制，为了兼容现有超声设备，因此面阵探头尺寸较小，造成视野有限，图像分辨率也比二维超声图像分辨率低。目前的容积成像主要应用于对心脏、子宫、膀胱等器官的三维成像，还未有对颅内组织的检测和观察。
综上所述，如何在保证实时性的同时，有效提高经颅超声成像的分辨率也是一个非常关键的难点。国际上，法国郎之万实验室基于随机光学重建显微镜成像技术，开发出联合造影微泡的超声经颅超分辨成像方法，其有三个特点：第一、实验目标为颅骨较薄的大鼠，因而采用高频超声成像而非低频；第二、采用低浓度、单个微泡进行定位追踪；第三、使用上万帧数据合成，获取时间较大。这些因素使得该项技术在临床使用中面临诸多问题。而以色列理工学院生物医学工程系研究小组在基于超分辨光学波动成像方法上，开发出可用于高浓度、低获取时间的联合微泡超声超分辨成像方法，但是其分辨率相较于单微泡追踪方法有明显不足，且其不适合进行微泡/血流的速度成像。更重要的，这些技术都是针对二维成像，而对于超声经颅容积三维成像，还未有明确的超分辨率成像技术。
因此，如何能在现有超声成像设备基础上，并结合射频数据采集装置，针对联合造影经颅超声脑血管分布，将动态三维超声成像技术应用于颅脑，在保留实时性的同时，提出一种小区域经颅血管容积三维成像算法来改善经颅容积成像的分辨率，是亟待解决的问题。
技术实现要素：
本发明的目的在于提出一种超声低频经颅容积超分辨率三维造影成像方法及系统，其将容积成像技术应用于经颅血管成像领域，实现对颅骨内部结构的实时三维成像，同时，利用体素最近邻方法(vnn)和三维马尔科夫链蒙特卡洛多目标追踪算法实现经颅血管的小区域三维超分辨率成像。
为了实现上述目的，本发明采用了以下技术方案：
一种超声低频经颅容积超分辨率三维造影成像方法，包括以下步骤：
1)利用二维平面换能器向颅内发射超声平面波，所述超声平面波于相互垂直的两个方向上进行扫描，通过接收所述扫描的回声信息，获得连续n帧超声造影回波射频数据；
2)对所述回波射频数据进行二维波束合成，得到波束合成后的射频数据bn(n＝1,...,n)；
3)依据波束合成后的射频数据，确定每一帧成像的空间规则体积像素网络wa×b×c，其中，a为x轴方向像素数目，b为z轴方向像素数目，c为y轴方向像素数目；
4)对每一体积像素点wa,b,c(a＝1,...,a；b＝1,...,b；c＝1,...,c)，简称体素点，遍历每一帧成像中所有扫描图像上所有像素点，计算扫描图像每个像素点距离该空间规则体积像素网络wa×b×c的距离，用距离体积像素点wa,b,c最近的像素值对该体积像素点wa,b,c进行插值；插值后得到与空间规则体积像素网络wa×b×c体积相同的三维体积像素网络w′a×b×c；
5)将得到的三维体积像素网络w′a×b×c沿y轴分割成n1幅二维图像bm(m＝1,...,n1)，n1为所述二维平面换能器次轴阵元数，对bm(m＝1,...,n1)中每一幅图像的n帧图像bm,n(n＝1,...,n)依次选取其中连续的nmed帧得到bm,％n(％n＝n-nmed+1,...,n)，对bm,％n用时域中值滤波函数进行处理，得到背景图像gm,n：
其中，distance(·)表示像素间距离，i表示像素点横坐标，j表示像素点纵坐标；
6)从原始图像bm,n中减去背景图像gm,n得到前景图像fm,n；
7)对前景图像fm,n进行伪迹去除及局部最大值处理后，将n1幅前景图像按照原体积像素网络wa×b×c的体积和分布拼接合成为三维数据，得到第n帧图像，然后将得到的每一帧中μn个像素点作为造影微泡的位置，记录造影微泡的位置坐标
8)将第n帧图像中所有μn个像素点的坐标组合，得到
9)在y中设定距离范围(≤1mm)，以各个造影微泡位置为中心，将小于距离范围的其他造影微泡位置预设为中心造影微泡位置可能产生轨迹的链接造影微泡位置；
10)利用马尔科夫链蒙特卡洛多目标追踪对造影微泡轨迹进行推演；
11)对推演的最终轨迹集进行成像。
优选的，所述二维平面换能器采用128或256阵元，主轴阵元中心间距为1mm～1.4mm，主轴有效孔径为16mm～22.4mm，次轴阵元中心间距为1mm～1.4mm，次轴有效孔径为16mm～22.4mm，中心频率为1.8mhz～2.2mhz，所述换能器的脉冲宽度小于等于10μs，带宽大于等于60％，相邻阵元串扰小于等于-30db，背材衰减小于等于-60db。
优选的，所述步骤1)中，成像帧率frt为0.8khz～1.2khz，以满足造影微泡的瞬态物理信息变化要求。
优选的，所述二维波束合成采用基于特征空间的最小方差自适应波束合成算法(例如，超声二维面阵的三维宽波束小区域快速空化成像方法)。
优选的，所述步骤4)中，所述最近的像素值应为距离体积像素点wa,b,c最近的扫描平面(扫描平面为每一帧三维扫描图像内的一幅二维图像平面)上该体素点(体素点即体积像素点)的投影的像素。在实现过程中可以添加距离的阈值量hmax，如果该体素点与对应最近的扫描平面的距离大于等于预设的阈值量，则说明该体素点离自由扫描图像较远，不需要对该体积像素点进行插值(体积像素点取值为0)，如此可以降低算法的时间复杂度。
优选的，所述伪迹去除包括以下步骤：
7.1)构建二维高斯核函数，通过二维卷积对前景图像fm,n滤波，得到滤波后的前景图像f’m,n：
其中，na表示二维高斯核函数的横向(x轴)尺寸，nb表示二维高斯核函数的纵向(z轴)尺寸，二维高斯核函数为：
其中，δz和δx分别为图像轴向和横向的像素尺寸，
7.2)设定阈值fthr，去除f'm,n中小于fthr的像素值，后经局部最大值处理后，以μn个像素点作为造影微泡的位置，记录其坐标，并可根据位置间距和时间差得到对应的速度信息；
优选的，所述阈值fthr为-20db～-15db，可通过滤波进一步去除尾迹。
优选的，所述步骤10)具体包括以下步骤：
10.1)根据抽取的推演状态，改变y中的造影微泡位置，并推演得到轨迹集ω＝[ω1ω2…ωm]，其中m为轨迹总数；
10.2)构造时域迭代状态方程：
其中a为过程矩阵：
δt＝1/frt，frt为成像帧率，
10.3)根据卡尔曼滤波，预测的微泡位置表述为：
其中，
10.4)统计计算轨迹终点数目zt，当前轨迹总数nt，新生轨迹起点个数at，被选中为轨迹的目标点数dt，虚警目标总数ft，并计算从上次迭代中衍生的新轨迹数目ct及未被选中为轨迹的目标点数gt：
ct＝nt-zt-at
gt＝nt-dt
10.5)计算当前推演状态下的先验概率值：
其中，pz及pd分别表示目标终止或消失概率及目标被选中为轨迹点的概率，λb及λf分别表示理想轨迹中目标点新生的概率值及虚警概率值；
10.6)假设虚警数目和推演路径均服从均匀分布，计算基于当前推演轨迹集和各轨迹目标点的似然值：
其中，m为推演轨迹总数，τk为当前推演轨迹中目标点总数，
10.7)计算当前推演状态下的后验概率值：
其中，p(y)表示探测造影微泡位置的先验概率；
10.8)以当前推演轨迹集ω为马尔科夫链的某状态，依据最大后验概率方法计算接受当前推演概率：
其中，ω′为前次推演轨迹集，ω为本次推演轨迹集(即当前推演轨迹集)；
10.9)产生服从均匀分布的随机数ξ，0<ξ<1，若ξ<a(ω′,ω)，则接受本次推演，将新轨迹集ω替换老轨迹集ω′；否则，拒绝本次推演，轨迹集ω′保持不变；
10.10)重复步骤10.1)-10.9)，直至完成迭代次数。
优选的，所述n的取值为80～200，推演迭代次数l为1500～3000，以期尽量达到稳态过程。
优选的，所述推演状态包括新生、消亡、更新、断裂、融合、延拓、减少和转换；推演状态的抽取包括以下步骤：计算迭代过程下轨迹总数m，若m＝0，则选中新生状态；若m＝1，则在除转换状态外的其他推演状态中随机抽取一种；若m>1，则在各推演状态中随机抽取一种。
一种超声低频经颅容积超分辨率三维造影成像系统，包括低频经颅专用超声平面换能器以及回波信号接收模块、脑血管超声造影成像模块、造影微泡提取模块、伪迹去除模块、造影微泡轨迹推演模块和超分辨率成像模块；其中：
所述回波信号接收模块用于对所述超声平面换能器接收的超声回波信号进行采样并转化为射频数据；
所述脑血管超声造影成像模块用于根据射频数据，采用二维波束合成及像素的体素最近邻方法，生成脑血管超声造影图像，并转化为空间规则体积像素网络的形式；
所述造影微泡提取模块用于去除脑血管超声造影图像中通过时域中值滤波生成的背景信息，并保留造影微泡的图像信息；
所述伪迹去除模块用于滤除保留造影微泡的图像中的非造影微泡的移动轨迹；
所述造影微泡轨迹推演模块用于通过马尔科夫链蒙特卡洛多目标追踪对造影微泡轨迹进行推演；
所述超分辨率成像模块用于根据推演得到的轨迹集进行脑血管造影超分辨率成像，其可以极大提高脑血管造影成像的分辨率，并提供可靠的血流信息。
优选的，所述低频经颅专用超声平面换能器采用上述二维平面换能器。
本发明的有益效果体现在：
本发明将vnn重建算法与基于马尔科夫链蒙特卡洛多目标追踪相结合，可以有效的提高经颅超声脑血管造影三维成像分辨率，可达到微米级或像素级分辨率的超分辨率成像，因此可以有效反映颅内毛细血管的分布情况。
进一步的，本发明所述二维平面换能器可以在颞骨的小区域范围内通过进行相互垂直的两个方向的扫描，进而在进行观测时能够在一定范围内收集一次扫描中的全部回声信息，实现实时的三维成像。其具有双向扫描、实时、低频、高分辨率及小区域成像等特点。
进一步的，本发明采用的三维图像恢复算法为像素的体素最近邻方法(vnn)和三维基于马尔科夫链蒙特卡洛多目标追踪的经颅超声脑血管造影超分辨率成像方法，可有效提高经颅容积成像的分辨率。
附图说明
图1为实施例列举的低频经颅超声面阵换能器的阵元排列示意图。
图2为容积换能器扫描示意图。
图3为经颅超声成像时颅骨作用位置示意图。
图4为经颅超声容积成像实验系统及流程图。
图5为容积换能器三维成像流程图。
图6为三维推演状态示意图。
图7为二维基于马尔科夫链蒙特卡洛多目标追踪算法利用轨迹集得到的最终的超分辨率成像结果；(a)为二维超分辨率成像结果；(b)为超分辨率成像与原始造影成像复合叠加结果。
具体实施方式
下面结合附图和实施例对本发明做进一步的详细说明。
本发明提出了一种用于经颅专用小区域容积成像换能器阵列设计，以及在此基础上的经颅容积超分辨率成像方法。
参见图1，经颅专用小区域容积成像换能器采用超声小区域二维平面换能器阵列设计：由n1×n2个阵元组成，每个阵元尺寸为l×l，阵元中心距离为d，焦点距离面阵几何中心为f。参见图2，所述容积成像换能器为经颅专用低频超声容积换能器，其采用二维面阵换能器，使用换能器时可以进行两个方向的扫描，能够在一定范围内收集一次扫描中的全部回声信息，然后进行联机重建三维图像。
以下述二维面阵为例：经颅专用小区域容积成像换能器为32×8共256阵元，主轴阵元中心间距为1mm～1.4mm，主轴有效孔径为16mm～22.4mm，次轴阵元中心间距为1mm～1.4mm，次轴有效孔径为16mm～22.4mm，中心频率为1.8mhz～2.2mhz，换能器的脉冲宽度小于等于10μs，带宽大于等于60％，相邻阵元串扰小于等于-30db，背材衰减小于等于-60db。符合实际设计要求，且提供多样化的工作模式。
本发明所使用的经颅专用低频超声容积换能器，此换能器在次轴方向上可以加大有效孔径，提高成像分辨率和灵敏度，并且可以在两个方向上进行扫描成像。同时，主轴有效孔径较小，可在区域较小的部位成像，能够很好的贴合颞骨。该换能器可以在颞骨的小区域范围内进行相互垂直的两个方向的扫描，进而在进行观测时能够在一定范围内收集一次扫描中的全部回声信息，实现实时的三维成像。其具有双向扫描、实时、低频、高分辨率及小区域成像等特点。
参见图3，经颅组织成像的观测位点位于整个颅骨最薄的位置，即图3中圆点标注的颅骨两侧的颞骨位置，可以看出颞骨的区域很小，而且并不平整，因此在颞骨处进行颅骨内超声成像会对换能器尺寸的要求较高，本发明的超声容积换能器可以较好的贴合在颞骨的位置进行成像观察。
本发明提出的超分辨率经颅容积成像方法中：容积换能器发射空间超声波束直接获得三维超声体数据。三维数据是通过容积换能器扫描平面的移动而获得的大量连续的二维断面图，然后每一断面的二维图像信息连同其位置信息一起被数字化后输入计算机。每一时刻，计算机获得的数据是由许多断面组成的三维数据信息，三维数据信息包含一系列的体积像素，每一像素既包括灰度值，也包括亮度值。最后计算机对这些三维数据进行图像的重建和数据处理，最终能够在屏幕上显示出来。容积成像的关键就在于将采集到的数据重建，然后进一步的进行图像处理，提高成像分辨率，从而得到分辨率较高的三维图像。本发明涉及的三维图像恢复算法为像素的体素最近邻方法(vnn)和三维基于马尔科夫链蒙特卡洛多目标追踪的经颅超声脑血管造影超分辨率成像方法，能够有效提高经颅容积成像分辨率。
参见图4，利用本发明设计的经颅专用低频超声容积换能器透过实验室制作的颅骨模型对血管模型进行微泡造影成像实验，由于真实颅骨难以得到并且不易于操作，因此实际实验中利用3d打印技术打印出与真实人体组织声衰减较为接近的模型(c.bai,m.ji,j.zong,etal.,“a3d-printedskullmodelwithcorrespondingacousticcharacteristicofhumanskullforultrasoundbrainimaginganddiagnosis”,proceedingsofinternationalsocietyfortherapeuticultrasound(istu)17thannualsymposium,nanjing,china,jun.,2017)，模型内含有脑组织仿体，造影微泡按照血管仿体的管道以流动泵推动持续流经脑组织仿体。低频超声换能器在开放式超声成像平台主机的控制下发射低频超声信号、接收回波信号并将回波信号传送至主机；主机将回波信号传输至开放式超声成像平台接收模块；接收模块对收到的回波信号采样、存储为射频采样数据并传送给计算机；计算机对回波射频采样数据采用容积换能器三维成像算法进行处理，随后依照标准的成像过程显示得到超声图像。
参见图5，容积换能器三维成像算法具体流程如下：
1)经颅专用小区域容积成像换能器放置于一侧颞骨附近(此区域对超声波的声衰减最弱)，向颅内发射超声平面波，获得连续n＝100帧超声造影回波射频数据，成像帧率frt为1khz；
2)对回波射频数据进行二维波束合成，得到所有的波束合成后的数据；
3)依据波束合成后的射频数据，确定每一帧成像的空间规则体积像素网络wa×b×c，其中，a为x轴方向像素数目，b为z轴方向像素数目，c为y轴方向像素数目。规则体素网络的构建需要考虑两个方面，体素网络的大小和单个体素点的大小，体素网络的大小需要能够包围所有的自由扫描图像；单个体素体积大小的设置不能过小，过小会导致重建图像中“分界线”的出现，同时加大算法计算量，但对重建图像误差与信噪比并没有提升。同时单个体素体积设置不能过大，体积过大的单个体素会造成重建图像由于体素数目较少存在明显的“锯齿”现象。体素体积大小的设置应该以图像采集间隔为标准，使体素体积大小近似保持在扫描图像的采集间隔附近；
4)对每一体积像素点wa,b,m(a＝1,...,a；b＝1,...,b；m＝1,...,c)，遍历每一帧成像所有扫描图像上所有像素点，计算每个像素点距离该空间规则体素网络wa×b×c的距离，用距离该体素点wa,b,c最近的像素值对该体素点进行插值；
5)距离该体素点wa,b,c最近的像素值应为距离该体素点wa,b,c最近的扫描平面上该体素点wa,b,c的投影的像素w′a,b,c。在实现过程中需要添加距离的阈值量hmax，如果距离该体素点最近的像素距离大于hmax，则说明该体素点离自由扫描图像较远，则不需要用空间自由扫描图像进行插值(即该体素点取值为0)，如此可以降低算法的时间复杂度；
对vnn算法，插值距离阈值hmax决定了距离体素点的最近像素点是否可以用于对该体素点赋值，hmax过小，则会导致重建图像中存在很多不能被插值的空白体素出现，影响三维重建结果，随着hmax增大，空白体素逐渐减小到扫描区域完全被插值。同时hmax设置过大，则可能对某些远离扫描图像的区域进行赋值，重建结果中引入噪声，所以hmax的数值的确定是在保证可以对连续扫描区域内体素赋值的情况下，hmax较小；
6)将得到的三维体积像素网络w′a×b×c沿y轴分割成n1(n1＝8，次轴阵元数目)幅二维图像bm(m＝1,...,n1)，对每一幅图像的n帧图像bm,n(n＝1,...,n)依次选取其中连续的nmed帧得到bm,％n(％n＝n-nmed+1,...,n)，对bm,％n用时域中值滤波函数进行处理，得到背景图像gm,n：
其中，distance(·)表示像素间距离，i表示像素点横坐标，j表示像素点纵坐标；
7)从原始图像bm,n中减去背景图像gm,n得到前景图像fm,n，微泡以及其他移动目标的图像都包含在前景图像中；
8)构建二维高斯核函数，通过二维卷积对前景图像fm,n滤波，以减少来自非微泡的
移动目标的伪迹对图像带来的影响，得到滤波后的前景图像f'm,n：
其中，na表示二维高斯核函数的横向尺寸，nb表示二维高斯核函数的纵向尺寸，二维高斯核函数为：
其中，δz和δx分别为图像轴向和横向的像素尺寸，
9)对f'm,n设定阈值fthr＝-20db，像素值小于fthr的信息去除，并进行局部最大值处理；随后将n1幅前景图像按照原三维空间分布合成为三维数据，然后得到每一帧μn个像素点作为造影微泡的位置，记录其坐标为
10)第n帧图像中所有μn个像素点的坐标组合，得到
11)在y中设定距离范围，以各个目标点(造影微泡位置)为中心，将小于距离范围的其他目标点预设为中心目标点可能产生轨迹的链接目标点；
12)参见图6，为三维基于马尔科夫链蒙特卡洛多目标追踪的经颅超声脑血管造影超分辨率成像方法过程中的不同推演状态，初始化设定路径推演迭代次数l＝2000，其推演状态包含八种，分别为：新生、消亡、更新、断裂、融合、延拓、减少和转换。
13)计算迭代过程下轨迹总数m，若m＝0，则选中新生状态(即新生状态概率1)；若m＝1，则在除转换状态外的其他推演状态中随机抽取一种，即各状态概率1/7；若m>1，则在八种推演状态中随机抽取一种，即各状态概率1/8；
14)根据抽取的状态，改变y中的目标点，以推演得到轨迹集ω＝[ω1ω2…ωm]，其中m为轨迹总数；
15)设
其中a为过程矩阵：
δt＝1/frt，frt为成像帧率，
16)根据卡尔曼滤波，预测的微泡位置可以表述为：
其中，
17)统计计算轨迹终点数目zt，当前轨迹总数nt，新生轨迹起点个数at，被选中为轨迹的目标点数dt，虚警目标总数ft，并计算从上次迭代中衍生的新轨迹数目ct及未被选中为轨迹的目标点数gt：
ct＝nt-zt-at
gt＝nt-dt
18)计算当前推演状态下的先验概率值：
其中，pz及pd分别表示目标终止或消失概率及目标被选中为轨迹点的概率，λb及λf分别表示理想轨迹中目标点新生的概率值及虚警概率值；
19)假设虚警数目和推演路径均服从均匀分布，计算基于当前推演轨迹集和各轨迹目标点的似然值：
其中，m为推演轨迹总数，τk为当前推演轨迹中目标点总数，
20)计算当前推演状态下的后验概率值：
21)当前推演轨迹集ω可被认为是马尔科夫链的某状态，依据最大后验概率方法计算接受当前推演概率：
其中，ω′为前次推演轨迹集，ω为本次推演轨迹集；
22)产生服从均匀分布的随机数ξ(0<ξ<1)，若ξ<a(ω′,ω)，则接受本次推演，新轨迹集ω将替换老轨迹集ω′；否则，拒绝本次推演，轨迹集ω′保持不变；
23)重复步骤13)-22)，直至完成迭代次数。
24)对推演的最终轨迹集进行成像。
在此基础上，结合脉冲逆转技术、编码发射技术、微泡追踪超分辨率成像技术、参量成像技术等，可实现不同功能的经颅血管三维成像。
参见图7，(a)为二维超分辨率成像结果，其中血管1内径1mm，血管2内径0.7mm；(b)为(a)中超分辨率成像与原始造影成像复合叠加结果，根据成像结果对比，表明利用基于马尔科夫链蒙特卡洛多目标追踪的经颅超声脑血管造影算法进行颅内超分辨率成像可提高图像分辨率。
三维基于马尔科夫链蒙特卡洛多目标追踪的经颅超声脑血管造影超分辨率成像方法为在二维的基础上将运算维度改为三维，基本原理保持不变，在选用vnn对扫描图像进行处理基础上，三维基于马尔科夫链蒙特卡洛多目标追踪的成像方法可以实现三维经颅脑血管的超分辨率成像。