平面波经颅超声成像相位校正及散斑跟踪∗

李倩岩 苏 畅 林伟军 郑音飞 武小晴

(1 中国科学院声学研究所 北京 100190)

(2 中国科学院大学 北京 100049)

(3 北京市海洋深部钻探测量工程技术研究中心 北京 100190)

(4 浙江大学生物医学工程与仪器科学学院 杭州 310027)

(5 之江实验室超级感知研究中心 杭州 311100)

0 引言

利用超声对脑组织进行成像，并对脑血管中的超声造影剂微泡进行速度检测，可以显示脑部血流动力学特征，对阿尔茨海默症[1]、脑卒中[2]等脑疾病的诊断具有重要作用。与目前临床脑血管成像中常用的计算机断层血管成像(Computed tomography angiography, CTA)[3]、磁共振血管成像(Magnetic resonance angiography, MRA)[4]等手段相比，超声检测具有更好的实时性，且设备便捷、成像无创无损。

传统的临床经颅超声成像主要透过颞骨、枕骨等颅骨较薄部位“声窗”进行，如经颅多普勒(Transcranial Doppler, TCD)和经颅彩色多普勒(Transcranial color-coded duplex sonography,TCCD)[5]，探测范围受限。Du等[6]提出利用超声发散波成像，可以扩大经颞窗的颅内成像范围。自Bercoff等[7]利用超快速平面波相干复合成像方法检测血流以来，颅脑血流超声成像技术得到了快速发展。Macé等[8]利用多角度平面波发射结合能量多普勒成像，实现了具有良好时空分辨率的大鼠脑功能成像，并进一步扩展至三维实时成像[9]。国内中国科学院深圳先进技术研究院等单位也在动物实验中实现了二维血流成像[10]。通过对高帧频图像中微泡的特征点定位和跟踪，Errico 等[11]提出了超声显微镜(Ultrasonic microscope, ULM)技术，可实现超高分辨率的微细脑血管成像。以上研究工作中为了增强回波能量，需要将颅骨削薄[12]，或者在血管中注入超声微泡造影剂[13]，但并未考虑颅骨对超声传播中相位变化的影响。

一方面，颅骨声阻抗与其他组织差异大、骨组织声吸收强、内部存在多孔结构等原因，使穿过颅骨的超声幅度发生衰减；另一方面，颅骨声速远高于附近其他组织，造成经颅超声的相位畸变，在对颅内目标进行超声成像和速度检测时会出现误差大、信噪比低等问题。因此，经颅超声成像中需要对颅骨导致的相位畸变加以校正，而颅骨厚度和内部声学参数的非均匀性使相位畸变随声波入射位置和角度而变化，增加了校正难度。针对经颅超声成像中的相位畸变校正，Vignon 等[14]提出一种声学自适应聚焦方法，利用颅脑两侧相对的超声换能器阵列结合时空反转滤波器，改善超声聚焦和成像效果。Guasch等[15]提出利用环绕整个颅脑的换能器阵列进行超声信号的全矩阵采集，并用全波形反演方法自适应更新颅骨声速模型，可实现颅脑高分辨率成像。这些方法需设计新的换能器阵列，而基于通用的医学超声换能器进行经颅超声成像的研究中，基本采用了利用颅脑CT 图像获取颅骨声学参数[16]并补偿相位的办法。Soulioti等[17]利用时间反转法分别在聚焦发射和接收波束形成过程中校正经颅超声相位，实现高分辨率成像。Jones 等[18]在被动声声像中、宋亚龙等[19]在平面波成像中分别对比了射线声学理论方法与时间反转法的相位补偿和成像效果，认为时间反转法精度更高，但计算耗时过长，而射线声学方法可以实现相对简单高效的成像校正。Wang等[20]、Sukhoruchkin等[21]利用数值仿真验证了射线声学校正方法可显著改善经颅聚焦扫描成像分辨率。Jiang 等在合成孔径成像中利用射线跟踪法校正超声相位[22]，又提出一种分层介质成像的全矩阵相位偏移成像方法[23]，通过数值仿真和体模实验验证了成像效果。以上研究均为对脑组织结构成像中的超声相位校正，最近Demené 等[24]采用ULM 技术进行了人体经颅血流成像，利用时间反转法校正了颅骨对超声相位的影响。

准确的脑血流成像对临床诊断具有重要意义，针对颅内血流超声成像中的相位校正问题，本文将脑血管中的超声造影剂微泡看作随血液流动而运动的散射点，首先采用超声平面波多角度相干复合成像和散斑跟踪方法对颅内点目标进行成像和运动速度/位移估计，再利用由CT 图像获取的颅骨声参数先验模型，采用近似的射线声学理论方法校正平面波经颅入射和反向散射波穿过颅骨时的相位畸变。分别通过数值仿真和体模实验，分析颅骨对目标点成像及位移/速度检测的影响，并验证相位校正方法对成像结果的改善效果。

1 基本原理

1.1 平面波相干复合成像

平面波成像一次发射即可覆盖较大区域，与传统聚焦扫描方法相比具有成像帧频高的优势。非聚焦发射会带来图像分辨率与对比度下降的问题，本文采用多角度相干复合方法[25]，在提高扫查帧频的同时保证成像质量。

超声换能器所有阵元以一定初始相位发射脉冲信号，使超声波按照设定的偏转角度传播，在成像目标区域内波阵面可近似看作平面波；声波在声阻抗发生变化的位置发生散射，散射点可看作被动点声源，产生的散射波以球面波形式传播，回到换能器并被各阵元接收，如图1所示。

图1 偏角平面波发射接收示意图Fig.1 Diagram of deflection plane wave transmission and reception

偏角平面波声传播路径如图2所示，假设发射平面波次数为m，θ1,θ2,···,θi,···,θm为平面波每次发射的偏转角度；假设换能器沿x轴排列，n为阵元总数，x1,x2,···,xk,···,xn为换能器阵元的横坐标，则根据成像目标区域内某一空间位置(x,z)与换能器阵元的距离关系，对各阵元的接收信号经延迟叠加可得到该位置对应像素点的灰度值：

图2 偏角平面波声传播路径示意图Fig.2 Diagram of the path of declination plane waves

其中，延迟时间τ包含平面波到达时间和散射波返回换能器阵元xk的时间，即

其中，假设介质声速均匀，声速为c时，τemit、τrec分别为

将不同角度发射时得到的图像相干叠加，即可得到最终的复合成像结果。

1.2 散斑跟踪运动目标速度矢量成像

在平面波成像得到的图像序列基础上，应用散斑跟踪方法，通过跟踪目标点回波信号特有的散射子斑点信息，可以得到目标点位移，结合成像帧频信息即可估算出运动目标的速度矢量。

假定散斑图像在短时间内除空间位置的变化偏移外稳定无变化，如图3所示，首先在超声图像上选取适当尺寸的采样窗口作为参考窗，在下一帧图像中取同一中心更大范围的窗口作为搜索窗。利用归一化的二维互相关匹配算法，在搜索窗内进行二维搜索，得出归一化互相关峰值所在位置，即参考窗与搜索窗的最佳匹配位置[26]。根据连续两帧图像之间数据窗与匹配位置的相对位移，以及获取两幅图像的时间间隔，即可计算出该处散射目标的运动速度矢量。移动采样窗口并重复以上过程，对整幅图像进行遍历，可得到整幅图像各点的运动速度。

图3 散斑跟踪原理图Fig.3 Schematic diagram of speckle tracking

1.3 近似声射线理论颅骨相位校正

颅骨声速与附近其他组织差异较大，使超声穿过颅骨时传播时间出现差异，因此超声波穿颅后相位发生畸变。且颅骨不同位置的厚度、声速分布不均匀，使超声相位变化的程度随声波入射至颅骨的位置和角度而有所不同。因此当采用平面波发射时，进入颅内的波阵面不再保持同一平面，探测目标点产生的散射波到达颅骨后也不再以球面波形式传播。若仍按前文所述在介质声速均匀的假定前提下进行成像和速度提取运算，会造成较大误差，甚至出现无法提取速度的情况。为了改善经颅成像和目标速度估计的效果，补偿颅骨造成的超声相位畸变，本文采用近似的声射线理论方法在每一帧图像的成像过程中进行补偿校正。已知图2为无颅骨时超声波入射和散射波被某一阵元接收时的声射线传播路径示意图。而当颅骨存在时，忽略超声折射造成的声波传播方向变化，用直线代替声传播路径，则入射波和散射回波分别以不同角度经过颅骨不同部位，如图4所示。计算该路径上由于颅骨声速变化带来的传播时间差异，即可进行成像校正。

图4 基于近似射线声学理论的相位校正示意图Fig.4 Schematic diagram of phase correction based on approximate ray acoustics theory

若无颅骨存在，根据前文分析，当换能器阵列以角度θ发射平面波，对成像区域内像素点(x,z)进行延迟叠加计算时，位于xk处的换能器接收波形的延迟时间可由公式(3)～(4)计算得到。而有颅骨存在时，入射波和散射回波的传播时间分别变为

其中，cskull(l)代表声波在颅骨内传播的声速，大小随位置发生变化；c是颅内脑实质的平均声速。用校正后的延迟时间代入公式(1)～(2)进行成像，对每个平面波发射角度下的回波信号均重复以上过程，再进行常规的相干复合运算，可得到校正后的多角度相干复合图像。校正后的图像序列可以直接利用1.2节中的散斑跟踪方法，进行速度矢量提取。

2 模型与方法

2.1 数值仿真

为了分析颅骨对平面波和散斑跟踪对点目标成像及速度估计的影响，并验证相位补偿方法的效果，建立如图5所示的介质模型，图中上方黄色线段标出了换能器阵列所在位置，中间非均匀灰色部分为颅骨，下方绿色圆点为用于成像和运动速度检测的散射目标点，红色箭头标出了目标点的运动方向。

图5 平面波颅内运动目标速度估计模型Fig.5 Estimation model of plane wave intracranial motion target velocity

换能器阵列为具有64 个阵元的平面线阵，位于距颅骨外缘约5 mm 处，中心频率5 MHz，阵元间距0.2 mm(约0.65波长)。颅骨模型是从颅脑CT图像中提取的，CT 扫查空间分辨率为0.48 mm，对CT 图像经旋转处理后截取了58×76 个网格的区域，再经插值处理使模型的网格间距符合数值模拟算法的要求，最终形成的模型区域网格点数为913×1201，网格点间距0.03 mm，模型区域对应的空间范围为27.4 mm×36.0 mm。将颅内外其他软组织看作均匀介质，设置为人体软组织平均参数；颅骨部分则根据CT 值(反映组织对X射线的吸收率，单位为HU)与颅骨声学参数的线性变换关系[16]，逐点计算对应的声速和密度，建立非均匀颅骨声参数模型，仿真所用具体介质参数见表1[27−29]。注意此处横波速度和密度参数用于声场的数值仿真计算，而在成像和相位校正中只考虑反射波的首波相位，只需用到介质纵波声速。

表1 数值仿真中用到的声参数Table 1 Acoustic parameters used in numerical simulation

为了比较和分析不同速度、不同方向运动目标的成像效果，设置了规则排列的散射目标转动模型，如图5中绿色点所示。每个目标为直径0.32 mm 的实心圆，声速和密度略大于背景介质声速与密度，见表1。目标大小约为成像所用超声波的一个波长，可将其看作点散射目标。目标点共40 个每5 个一组，以成像区域的中心点(0,20 mm)为圆心，沿周向8个方向均匀排列，即以x轴正方向为起始角，每组目标点对应的周向角度依次为0、π/4、π/2、3π/4、π、5π/4、3π/2、7π/4；组内5 个点所在半径依次为1 mm、2 mm、3 mm、4 mm、5 mm。这些点组成的圆盘绕圆心顺时针旋转，角速度为100 rad/s，则从圆心向外目标点的线速度依次为0.1 m/s、0.2 m/s、0.3 m/s、0.4 m/s、0.5 m/s。

本文使用fortran 语言编写的空间四阶精度与时间二阶精度的交错网格时域有限差分程序，求解线性弹性波方程，模型外边界处利用卷积完美匹配层(Convolution perfectly matched layer, CPML)处理以消除反射。计算中时间步长取0.005 µs，空间步长为0.03 mm(约等于声波在背景介质中传播波长的1/10)，满足数值计算的收敛性和稳定性条件，总迭代步数为10001，计算总时长50 µs。本文重点考察颅骨对声波相位畸变的影响和校正，而没有对幅度衰减做处理，因此仿真算法中没有考虑颅骨对声波能量的吸收作用。实际中颅骨对声波幅度的衰减是影响成像的重要因素，作者将在后续工作中考虑幅度衰减的补偿办法。

利用以上仿真软件，分别在有颅骨和无颅骨的模型中，使换能器阵列以−8° ～8°、间隔2°、共计9 个不同偏转角度发射平面波，仿真计算声传播过程，记录换能器各阵元处的回波数据。完成一组计算后，根据模型中各目标点不同的速度和方向，计算目标点移动后的位置，重复多角度平面波发射的声场仿真计算，记录回波数据。对每组仿真计算得到的回波数据，利用前文所述的相干复合成像方法和相位校正方法分别成像，成像区域空间范围为横向大小为12 mm(−6～6 mm)，纵向大小为16 mm(12～28 mm)。对目标点位置发生变化的前后两帧图像，设成像帧频为1 kHz，即两幅图像成像时间间隔为1 ms，利用散斑跟踪法进行速度矢量成像。散斑跟踪计算过程中，采用的参考窗大小为1 mm×0.3 mm(横向×纵向)，搜索窗口大小为2 mm×1.3 mm(横向×纵向)。

2.2 体模实验

体模实验平面如图6(a)所示，主要装置包括计算机及超声成像系统、线阵换能器、三轴扫描装置及固定夹具、仿真颅骨及成像线靶等。超声成像系统(深圳华声，四叶草)具备64 个独立通道，可控制线阵换能器(L15-4)的64个阵元发射超声并接收回波信号。实验时由计算机向成像系统发送指令，设置发射和采集参数，成像系统采集的回波信号经正交调制(IQ)后传输至计算机存储，用于后续的成像处理。线阵换能器及仿真颅骨利用夹具固定于水槽中，成像线靶与步进电机驱动的三轴扫描装置相连，由计算机控制扫描装置移动线靶位置，以模拟颅内的运动目标。

换能器、仿真颅骨和线靶的位置关系如图6(b)所示，换能器与仿真颅骨距离约1 mm 左右，与线靶距离约25 mm。仿真颅骨为环氧树脂制成的平板，厚度7.06 mm，密度为1.25×103kg/m3，声速为2700 m/s。为了增强成像目标反射能量，成像线靶由平行排列的铜丝组成，每条靶线直径约1 mm，延伸方向与换能器阵扫查平面垂直，沿平行于换能器表面方向排列，间距约2～3 mm。实验过程中，成像线靶由扫描装置带动向靠近换能器方向移动，以模拟颅内的运动目标。

图6 实验装置示意图Fig.6 Schematic diagram of the experimental device

由于成像系统和扫描装置配合的实时采集和瞬时运动速度难以控制，本文采用了非实时的定点测量方式，即固定线靶位置进行扫描成像，使线靶移动一定距离后再进行成像，从两幅图像中计算目标点位移，来代替速度估算。由于散斑跟踪法本质上是跟踪前后两帧图像上特征散射点的位移，再除以成像帧之间的时间间隔来计算速度，对位移矢量的估算精度即决定了对速度的估算效果，因此实验中用颅内目标点的位移矢量成像代替速度成像具有一定合理性。

换能器阵列的64 个阵元以−8° ～8°等间隔递增的9 个角度发射平面波，脉冲中心频率为5 MHz，相邻阵元中心距离为0.2 mm。采集的回波数据利用多角度平面波相干复合方法进行成像处理，成像区域范围内检测到3个目标点，每得到一帧图像，控制成像线靶向靠近换能器方向移动0.2 mm，再重复多角度发射和采集成像过程。得到两帧目标点处于不同位置的图像序列后，对两帧图像进行散斑跟踪处理，得到最终的目标位移矢量分布图。散斑跟踪计算过程中采用的参考窗大小仍为1 mm×0.3 mm(横向×纵向)，搜索窗口大小为2 mm×1.3 mm(横向×纵向)。分别在无颅骨、有仿真颅骨条件下采集数据并进行成像和位移计算，对有仿真颅骨的数据进行相位校正处理并对比结果。

3 结果与讨论

3.1 数值仿真结果及误差分析

数值仿真结果如图7所示，从左至右依次为无颅骨、有颅骨未做校正、有颅骨并利用近似射线法校正后的成像结果，其中灰度图为平面波相干复合成像法得到的一帧图像，红色箭头为对两帧图像利用散斑跟踪法得到的速度矢量分布图。

对比图7(a)、图7(b)的散斑成像和速度矢量分布图可以直观地看到，颅骨的存在使图像质量骤降，对速度的估计结果产生巨大影响：(1)散斑图像伪像增多、图像分辨率与对比度降低，使得散斑跟踪过程中无法正确分辨具有有效速度的散射点，出现速度缺失。(2)由于颅骨内声速传播较快，速度分布整体向换能器方向发生偏移，偏移距离约为1.4 mm。由图7(c)可以看出，进行近似射线法的相位补偿后：(1)图像质量提升，能够区分具有有效速度的散射点。(2)通过校正将散射点的偏移距离减小至约0.3 mm，同时速度大小与速度方位的估计值更加准确。

图7 多角度相干复合平面波散斑跟踪速度矢量估计分布图Fig.7 Multi-angle coherent composite plane wave speckle tracking velocity vector estimation distribution

为了定量对比速度大小的估计结果，对以上3种情况下得到的速度矢量分布图，在各散射点实际位置附近1.4 mm×1.4 mm 范围内取估算的平均速度，并对速度矢量取模后，与目标点真实运动速度大小做对比，如图8所示，其中图8(a)～(c)分别为无颅骨、有颅骨未校正、有颅骨并做校正3 种情况下的速度矢量模的估计值与真实值的对比。根据模型设置，在圆盘模型周向8 个方向上，每条半径方向有5 个目标点，分别具有0.1 m/s、0.2 m/s、0.3 m/s、0.4 m/s、0.5 m/s 共5 种线速度，在图8中用星号标出，对每个目标点的速度估计值用蓝色点表示。图8(d)为3 种情况下速度矢量模估计值与真实值之间的误差分布对比。

图8 3 种情况下不同半径速度估计幅值大小及误差对比结果图Fig.8 Comparison results of estimated amplitudes and errors of velocity at different radii in three cases

从图8(a)～(c)中可以看到：(1)颅骨的存在使部分有效速度点丢失，出现估计值为0 即未检测到的散射点，如图8(b)中红色箭头所示，这些点在校正后均得到恢复。(2)颅骨的存在使速度大小估计值发生明显偏差，就可检测到的散射点而言，最大误差已达60%以上，难以进行有效分辨，校正后，估计值与真实值偏差明显减小。(3)近似射线声学校正方法可提升估计精度，经计算，无颅骨时速度大小估计相对误差为5.07%，存在颅骨时速度估计相对误差为54.98%，近似射线法校正后速度估计相对误差减小至11.94%。

根据图8(d)，对比目标点以不同线速度大小运动时，速度估计值的偏差情况，可以看到：(1)无颅骨存在时，速度估计精度与散射点实际速度大小无明显关系，速度大小估计基本准确。(2)当颅骨存在时，目标点移动速度越快，速度估计误差越大，这可能与目标点所在位置有关，由于模型设置，线速度大的目标点均位于远离成像区域中心处，可能对成像精度有一定影响。(3)利用近似射线法校正以后，虽未完全恢复至无颅骨时的速度估计精度，但基本克服了颅骨的影响，速度估计误差减小，且与目标移动速度的大小不再相关。

进一步地，为了分析对速度方向的估计准确度，将速度矢量分解为横向、纵向分量分别比较。图9(a)～(c)显示了横向速度分量估计值与真实值的差，随目标点的真实横向速度分量大小的分布；图9(d)～(f)为纵向速度分量估计值与真实值的差，随目标点的真实纵向速度分量大小的分布；真实速度分量大小在图中用星号标出，对每个目标点的速度分量估计误差用粉色点表示。图9(a)、图9(d)为无颅骨情况，图9(b)、图9(e)为有颅骨未校正的结果，图9(c)、图9(f)为有颅骨且做校正后的结果。需要注意的是，由于有颅骨未做校正时误差较大，为了清楚看到每种情况下的误差分布情况，图9(b)、图9(e)坐标纵轴范围是另两种情况的6倍。

图9 横纵速度分量估计值误差对比图Fig.9 Comparison chart of the estimated errors of horizontal and vertical velocity components

观察图9可以发现：(1)比较无颅骨时的横向和纵向速度分量误差分布情况，二者变化范围比较接近，说明散斑跟踪方法估算速度的精度不受运动方向影响，这也是散斑跟踪法与多普勒方法的主要差别之一。(2)有颅骨存在时，速度横向和纵向分量的估计误差均增大，校正后两个方向速度分量的估计精度都得到提升。(3)校正后，速度横向分量误差总体上大于纵向分量的误差，且大部分向x轴正方向偏移，目前尚不清楚这种偏移趋势是否与所选取的颅骨模型声学参数和厚度分布特点相关，后续工作中还需对更多不同部位颅骨的模型进行计算分析。出现这种结果的另一个可能原因是在近似射线法中忽略了声波折射角度，造成相位校正的误差。

3.2 实验结果及误差分析

实验得到无颅骨、有仿真颅骨未做校正、有仿真颅骨并经近似射线法校正处理后的成像及目标位移矢量分布结果如图10 所示，将各图内3个运动目标点从左到右分别记为点 1○2○3○。对比图10(a)～(c)，由于仿真颅骨的存在，目标点位置及移动位移估计值的准确度均出现偏差。

对比图10(a)～(c)，可以发现：(1)图10(b)中目标对应的亮斑位置向换能器一侧靠近，与无颅骨时偏差约为3 mm，使用近似射线法校正后，亮斑回到与无颅骨时基本相同的位置。(2)图10(b)图像分辨率下降，横向方向更为明显，1○2○两点在图中几乎显示为连续界面，校正后图像分辨率与对比度得到提升。(3)仿真颅骨的存在导致对目标点位移的大小和方向估计均出现了偏差，近似射线校正方法可以改善目标点位移估计的效果。3 个目标点的实际位移均为向正上方移动0.1 mm；无颅骨时对3 个目标点位移大小和角度估计值的平均误差分别为0%和0%；有仿真颅骨时，位移大小和角度估计值平均误差为16.47%和27.67 %；校正后，误差降低为1.13%和1.41%，检测精度明显提升。3 种情况下3个运动目标的具体数据见表2。

表2 位移估计结果Table 2 Displacement estimation results

图10 经颅实验速度估计结果图Fig.10 Estimated results of transcranial experimental velocity.

对比实验与数值仿真结果，实验结果中对位移估计的误差小于数值仿真结果，其原因之一可能是条件所限，实验中仅设置了朝单一方向运动的3 个靶点；另外，实验所用仿真颅骨为均匀介质且为厚度固定的平板，而非内部声速、厚度非均匀的真实颅骨。

4 结论

本文利用超声平面波多角度相干复合方法对颅内目标进行成像，在高帧频图像序列基础上利用散斑跟踪法估算目标的位移或速度矢量。针对颅骨造成的超声相位畸变，利用由CT 图像获取的颅骨声速先验模型，采用近似射线声学理论方法进行校正。数值仿真结果表明，颅骨的存在会使目标点对应的图像散斑偏离实际位置，同时分辨率和对比度降低，图像质量的下降导致对运动目标位移/速度的估算误差高达54.98%。基于准确的颅骨声速先验模型，利用近似的射线声学理论方法可以改善成像效果，在校正目标点图像位置偏差的同时，目标点运动位移/速度的大小和方向估计精度提升近40%。体模实验中对目标位移大小和角度估计的平均误差在校正前分别为16.5%和27.7%，校正后降至1.1%和1.4%。进一步验证了相位校正的效果。

致谢文中数值仿真所用颅脑CT 图像由浙江大学医学院附属第四医院提供。