基于矢量场预估的二维超声弹性成像方法

周佳璐，王韬，杨军

（中国医学科学院北京协和医学院生物医学工程研究所，天津300192）

1 引 言

传统临床中对于软组织中的囊肿，脂肪瘤等存在弹性差异的病变检测一般采用触诊的检测方法，即对组织施加低频压力，通过触觉对病变组织进行定性测量与估计。由于该方法简单易用性，被广泛应用于皮肤，乳房，前列腺等部位疾病的检测中。但触诊本身的客观性与准确性则有很大的局限，特别是对于病变较深区域或弹性差异较小的病变组织很难进行测量［1］。

准静态超声弹性成像是近年来发展较快的一种新的超声成像技术。其通过对组织施加外部压力得到组织压缩前后的超声回波信号，并结合数字信号处理或者数字图像技术获取组织内部不同空间场上的应变信息，从而间接反映组织弹性模量等力学属性的差异［2］。一方面可以增强传统超声成像系统的成像精度，一方面也可以获取传统方式获取不到的病变信息。

超声弹性成像方法的主要流程一般如下：（1）施加外力获取前后对比信号；（2）求取图像位移场；（3）对位移场进行差分求取应变场获取组织应变［3］。最早的超声弹性成像法主要采用一维成像方法，即只考虑组织受压后的轴向位移（沿声波传播方向），其主要原因在于求解组织位移场时二维成像的计算量较大，且易引入较大的计算误差，无法满足算法的实时性要求，临床应用性不高［4］。然而组织在受外力发生压缩形变时，其应变特性是二维的，如果仅考虑轴向位移，计算得到的结果并不全面，精度不高。近年来出现的二维成像方法如动态规划算法（DP）［5］，快速互相关算法（FNCC）［6］，块匹配预估算法（BMAPS）［7］，光流法（OF）［8］等同时考虑了弹性形变时的横向位移，对弹性超声的临床应用有一定的推动作用。这些算法分为两类，一类是基于位移预估的后匹配算法，如块匹配预估算法（BMAPS）；一类是基于快速匹配方式的算法，如快速互相关算法等。基于快速匹配方式的算法其匹配效果较差，精度较低，临床使用较少。而现有的基于位移估计的算法在泊松比在0.5附近的情况下效果较好，即其横向位移不可过大，在泊松比小于0.45的情况下很难获得较好的结果。

本研究提出了一种基于矢量场预估弹性超声二维成像方法，采用光流法（OF）进行矢量场预估，采用相位相关算法（PS）进行位移估计，最后采用最小二乘法求取应变值并实现二维成像。实验采用有限元分析法进行数据仿真与结果对比，并通过对不同泊松比的仿真结果对比证明了该算法可实现较低泊松比情况下成像的清晰稳定。

2 方法

首先获取包络线特征峰值指纹信息，采用改进的光流跟踪法对特征指纹进行跟踪并进行矢量场预估；然后根据预估矢量场采用相位相关算法更改所跟踪的RF线对位移场进行精确计算；最后通过最小二乘法滤波求取应变场实现准静态二维弹性超声成像。从而可以获得较高精度的位移场并保证了大尺度横向应变。

2.1 跟踪特征获取

在图像跟踪中，一般采用角点跟踪法，以图像角点作为光流特征点进行光流跟踪，这种方法一般适用于分辨率清晰，图像轮廓分明的普通图像［9］。而超声获取到的射频信号（RF）归一化二维图像特征在于纹理复杂，轮廓并不清晰难以获取相关角点。因此，提出一种基于包络线峰值的特征指纹信息作为光流特征点的方法进行光流跟踪。

首先对RF信号进行归一化处理，得到RFnorm，有：

其中：min（RF）为射频信号最小值，max（RF）为射频信号最大值。

然后求取各条RF信号的包络线C，有包络线C满足：

联立方程解得上包络线信号Cup与下包络线信号Cdown，将RF信号进行分区划块，取各区块RF线上下包络线峰值与谷值作为特征跟踪点形成特征指纹信息进行光流跟踪，见图1，左图为立体图，右图为俯视图，横坐标表示第1～5条RF信号，纵坐标表示信号个数，各点颜色表示包络信号的波峰波谷位置。

2.2 改进的光流跟踪

由于RF信号具有纹理特性，如果直接采用传统光流法，其跟踪结果置信度不高，因此提出了一种将光流跟踪与指纹判别检测算法相结合的方式对光流跟踪结果进行判断以去除无效跟踪信息。

首先，采用光流法对特征点划分区域进行跟踪，设各线RF信号包络线受压后峰值A变化较小且相对位置特性不变。

设第x条RF线，第y点散射点（x，y）在t时刻的像素值（RF信号幅值）为 A（x，y，t），设压缩后该点信号幅值表示为 A（x+dx，y+dy，t+dt），o（n）为2阶无穷小，则有：

假设各特征点局部区域内包络信号值恒定不变，即有在各特征点区域内，使得下式O最小即可获得该特征点区域速度向量

图1 采用包络线峰值作为光流跟踪特征点Fig 1 The envelope peak value is used as the feature points of optical flow tracking

公式中，Ax1，Ay1，At1分别表示跟踪点（xt，yt）幅值在x，y，t轴的变化率。

然后采用之前作出的特征点矩阵作为指纹特征，对跟踪结果进行匹配判断，去除不符合矩阵排布特征的跟踪结果或置信度1/O较低的结果，其平均值为该区块的整体速度向量，则矢量跟踪场结果为：

其中，i11，i12，imn为各区块中有效值的个数。

2.3 相位相关与应变计算

根据光流跟踪预估偏移结果（u，v）进行偏移。由于u=dx/dt，压缩前后只有两帧数据，则有dt=1，可得u=dx即光流跟踪时横向像素偏移量；同理v=dv/dt，则有v=dy即光流跟踪时纵向像素偏移量。在光流跟踪处理中，横向像素数等同于RF线总数，纵向像素数目则是超声RF信号总采样点个数，设采样频率为fs，则有纵向像素偏移量转换为时间轴偏移量如下dt=dy/fs=v/s：。设原信号为RF，压缩后信号为RF。设压缩前某段信号为RFi（t1，t2），代表第i条RF信号t1到t2之间的一段信号，则压缩后，该段信号对应的进行相位相关精确匹配的压缩后的信号段应为 RF（i+u）（t1+v）/f2，t2+v/fs）。

图2 基于特征指纹的光流矢量跟踪场示意图Fig 2 Sketch of optical flow vector tracking field based on feature fingerprint

通过希尔伯特变换对压缩前后RF信号进行处理，获取含有相位信息的解析信号，有：

然后，利用跟踪矢量场作为相关计算的先验信息，设声速为s，取预估矢量场横向位移m=Mu（n，t），纵向时延 tP=Mv（n+m，t）/s，则有压缩前后 RF解析信号对应的预估偏移信号为：

超声频率ω0已知，则通过相关函数相移角度计算两端信号之间的时延：

则有最终计算得到各时间点t横向位移为△x=Mu（n，t），纵向位移为△y=s△t。

采用最小二乘法进行应变计算［10］，采用分段线性拟合算法对位移场求差分进行应变拟合即可求得应变场。

3 实验

3.1 仿真实验设计

采用有限元仿真软件COMSOL对模型应变情况进行仿真，采用FIELD II声场仿真软件对超声数据进行仿真处理［11］。

首先用COMSOL建立一个边长为80 mm的正方形二维模型，其内部中心包含一个半径为10 mm的圆。设置正方形模型材料与人类组织弹性纤维参数一致：杨氏模量为2MPa，泊松比为0.45。设置圆形材料为病变组织，杨氏模量设为正方形模型材料的5倍：杨氏模量为10MPa，泊松比为0.43［12］。设计受压模型即正方形底部为有限元分析边界固定条件，正方形上部为刚性连接，受向下的均衡力为0.1 N/m2。

有限元分析采用Delaunay方法对模型进行自由三角剖分，共产生40000个自由度扩散点，其仿真结果见图3：

图3 （a）.有限元分析模型；（b）.模型受压后Y轴位移场；（c）.模型受压后X轴位移场Fig 3 （a）.Finite element analysis model；（b）.Y axial displacement field after compression；（c）.X axial displacement field after compression

采用FIELD II声场仿真软件设计超声探头，其中心频率为3.5 MHz，采样率为100 MHz，线阵超声探头共计192个阵元，每组发射64个采样有效阵元，阵元宽度为声波波长，阵元间隙为波长1/20，（当阵元间隙与阵元宽度之和大于1/2波长时即可避免阵元间相互作用）。阵元长度设置为5 mm，并将每个物理阵元划分为10个计算阵元。设计声场扫描模型内部分布40 000个随机散射点，并根据模型有限元分析结果建立声场仿真模型见图4。

图4 （a）.压缩前声场模型；（b）.压缩后声场模型Fig 4 （a）.sound field modelbeforecompression；（b）.sound field model after compression

3.2 仿真实验结果与对比

将矢量场预估算法计算得到的纵向应变结果与有限元分析结果进行对比，比较其偏差值作为算法评价依据，见图5。

对于传统一维弹性超声成像算法，BMAPS块匹配预估算法，矢量场预估二维成像算法在0.4～0.48不同泊松比情况下进行对比，结果见表1。

由结果可看出，传统一维弹性超声算法在较低泊松比情况下其计算误差成指数增长，经过数据观察，其只能在中心RF线上保证较好的计算结果，而本算法在较低泊松比条件下的计算误差也低于BMAPS块匹配预估算法。实验证明，本算法是一种清晰度更高，鲁棒性更好的超声弹性成像方法。

表1 不同泊松比情况下各算法的偏差性能对比Table 1 Comparison of the deviation performance of each algorithm under Different Poisson's ratio

图5 （a）应变计算结果（b）有限元分析计算结果Fig 5 （a）strain calculation results（b）finite element analysis and calculation results

4 总结

本研究旨在提出一种高清晰度高鲁棒性的超声弹性成像方法，采用光流跟踪法对包络线峰值指纹特征进行跟踪，以跟踪矢量场结果作为先验信息对超声RF信号进行相位相关计算，一方面实现横向位移的准确估计，一方面也避免了相位相关计算时的相位混叠问题。与传统一维相位相关法相比，本算法在计算相关函数之前将跟踪矢量场位移作为预估量加入计算，避免在存在横向位移时依旧同一条RF线上进行信号匹配，与近年来的二维相位相关法相比，采用光流法预估矢量场大大提高了横向位移预估的精确性，可保证在软组织泊松比较小的情况下依旧可以进行准确追踪。实验结果最终也证明了该算法的准确有效性，对于超声弹性成像技术的系统研发具有一定的推动作用。