二维Savitzky-Golay数字差分器剪切波弹性成像方法

潘晓畅，冯乃章，陈思平

1)医学超声关键技术国家地方联合工程实验室，广东省生物医学信息监测与超声成像重点实验室， 深圳大学生物医学工程学院，广东深圳 518060；2)深圳开立生物医疗科技股份有限公司，广东深圳 518052

剪切波弹性成像能定量反映生物组织内部的硬度，且具有非侵入、重复性好和操作方便等优点[1-2]，现已用于乳腺肿瘤[3]和肝纤维化[4]等疾病的诊断研究中. 剪切波弹性成像的基本原理是超声探头在生物组织中发射一个高强度和持续时间长的脉冲，通常称之为激励脉冲，该脉冲能在生物组织中产生声辐射力. 声辐射力在达到一定的强度后就能在生物组织中激励产生剪切波[5]. 超声探头在发射激励脉冲之后，立即发射检测脉冲来检测剪切波在生物组织中的位置，进而计算出剪切波在生物组织中的传播速度. 通过剪切波在感兴趣区域内的传播速度分布来反映组织内的杨氏模量分布[6].

在线弹性、各向同性、不可压缩及均匀的半无限大空间中，剪切波动的传播速度c、 材料密度ρ和杨氏模量E之间存在如下关系

E=3ρc2

(1)

在生物组织中也采用式(1)近似计算组织中的杨氏模量分布[6]. 剪切波传播速度估计通常采用剪切波到达时刻法(time of flight, TOF)[7-8]. 先计算出一段时间内感兴趣区域内各个点的运动速度或位移，通过每个时刻感兴趣区域内所有点的运动信息来计算剪切波到达每个点的时刻. 计算剪切波到达时刻的方法有：波形互相关(cross correlation, CC)法、剪切波峰值到达时刻(time-to-peak, TTP)和剪切波坡度峰值到达时刻(time-to-peak-slope, TTPS)[2,9]. ROUZE等[9]的研究表明，TTPS相比TTP和CC受剪切波传播中产生的回波的影响更小，计算出的剪切波速度更准确.

剪切波是一种横波，其传播方向与发射声场方向垂直，在生物组织中是从激励位置向两边传播. 在计算出感兴趣区域内各个位置的剪切波达到时刻后，通过对感兴趣区域内横向的每条线做差分运算可以计算出剪切波的传播速度. 由于直接做差分会带来噪声的放大，现有的剪切波弹性成像中普遍采用线性回归(linear regression, LR)来对一小段数据用线性函数进行拟合[2,9]. 通过拟合得到的直线斜率作为差分比直接差分鲁棒性更好，计算出的剪切波速度分布更平滑、噪声更小. 线性回归所用的数据长度越长，得到的剪切波弹性图越平滑、噪声越小，同时空间分辨率也越低[9].

估计出的剪切波到达时刻若出现错误的异常点，会影响剪切波速度估计的结果. 为提高剪切波速度估计的准确性和鲁棒性．WANG等[7]将随机抽样一致法(random sample consensus, RANSAC)应用到剪切波速度估计中，有效地去除了估计的剪切波达到时刻中的异常点. 然而，RANSAC的计算量极大，计算一帧剪切波弹性图需要186 s，不适合用于实时的剪切波弹性成像[8]. ROUZE等[8]将Radon变换用于剪切波速度估计中，该算法具有与RANSAC相当的鲁棒性，且计算速度快. 然而，使用Radon变换的前提是剪切波在各向同性的材料中进行传播. 对于肝硬化等弥散性病灶，组织是接近各向同性的，剪切波弹性成像中使用Radon变换能取得较好的结果. 但对于乳腺和甲状腺中的肿瘤性病灶，组织是各向异性的，使用Radon变换会使估计出的剪切波速度有较大的误差.

为提高剪切波弹性成像的质量，AMADOR等[10]从时空域搜索剪切波波前的峰值，相比传统的空间域搜索提高了剪切波波前峰值位置的计算精度. LIPMAN等[11]将三维方向滤波器(包含时间维度)用于方向滤波中，相比现有的二维空间滤波进一步去掉了反射波伪像，提高了剪切波弹性图的对比度. YOON等[12]通过多角度激励剪切波，对各角度分别进行剪切波速度估计后再进行空间复合，减少了剪切波弹性图中的伪像.

为提高剪切波弹性成像质量，本研究对剪切波速度估计这一环节进行改进. 在获取剪切波波前位置后，将二维Savitzky-Golay(SG)差分器用于剪切波速度的差分计算中. 一维SG差分器[13-14]现已广泛用于准静态弹性成像中计算轴向应变，其本质上是线性回归的快速算法. 二维SG滤波器相比一维SG滤波器具有更高的鲁棒性，更少地受异常数据的影响[15]. 本研究采用标准的弹性仿体进行实验，比较了二维SG差分器和线性回归的剪切波弹性图的质量. 对弹性图像的质量采用超声弹性成像中普遍采用的信噪比(signal-noise ratio, SNR)和对比度噪声(contrast-to-noise ratio, CNR)进行评估．

1 实验设备

实验所用超声仿体是美国CIRS公司的弹性仿体(model 049A)，如图1. 仿体的材料特性与生物组织相近. 仿体中有不同弹性模量的内含物. 本实验所用的内含物尺寸约为10.4 mm，弹性模量为(80±12) kPa(均值±标准差). 仿体中背景的弹性模量为(25±6) kPa.

图1 本实验所用弹性仿体Fig.1 (Color online) The elastic phantom used in the experiment

实验所用超声设备是Verasonics V1系统，该系统被广泛用于剪切波弹性成像中[16-17]，所用超声探头是 L7-4线阵探头, 系统的采样频率为36 MHz. 剪切波发射的激励电压为75 V，发射800个周期的长脉冲. 剪切波的激励设置放在感兴趣区域的左侧. 超声探头采集的信号经过数字采样和波束合成后形成512条线的同相和正交(in-phase and quadrature, IQ)信号. IQ信号的横向线间距为74 μm，每条线的轴向采样点间距为77 μm. 超声探头以聚焦方式发射一次高强度、长时间的激励脉冲后，再以平面波的方式多次发射检测脉冲. 检测脉冲发射的时间间隔为 83 μs.

2 实验方法

2.1 剪切波达到时刻估计

超声探头发射激励脉冲后，会在生物组织中产生声辐射力，进而激励生成剪切波. 然后超声探头发射检测脉冲来计算剪切波的传播速度. 本研究中，剪切波的传播速度通过组织的轴向运动速度来进行计算. 组织的轴向速度由前后两次检测脉冲发射后得到的IQ数据计算出，采用KASAI等[18]提出的自相关算法为

(2)

其中，c为生物组织中的声速；f0为发射的检测脉冲的中心频率；I和Q分别为IQ信号的实部和虚部；m为IQ信号轴向上的点；M为计算一个点的运动速度所用的信号的轴向长度；n为当前检测脉冲的发射序号；v为计算出的组织的轴向运动速度.

发射一系列的检测脉冲后，用自相关方法能计算出感兴趣区域内所有点的轴向速度u(x,y,t). 这里x和y分别为感兴趣区域内的横向和轴向坐标，t为检测脉冲所发射的时刻. 对于N次发射检测脉冲，能计算出N-1组速度. 为提高剪切波速度估计的准确性，需要对组织运动速度v进行方向滤波来去除剪切波传播过程中产生的反射信号. 本研究中方向滤波采用的是DEFFIEUX等[19]提出的在频域内处理的方法. 对感兴趣区域内不同深度的横向线u(x,y0,t)分别进行滤波，这里y0指感兴趣区域内某一个深度位置.

2.2 二维SG差分器

在计算出感兴趣区域内各点的剪切波达到时刻T(x,y)后，用二维SG差分器在感兴趣区域内做二维卷积，即可得到剪切波达到时刻的横向差分Tx(x,y). 在剪切波达到时刻T(x,y)中取一小块区域 (2M+1)×(2N+1)， 用x方向p阶和y方向q阶的曲面进行拟合

(3)

其中，x取值为-M～M；y取值为-N～N；aij为曲面拟合系数．a00为曲面平滑系数，a10为曲面横向一阶导数.a10即为所需要的剪切波峰值达到时刻的差分. 用最小二乘法进行曲面拟合为

(4)

计算拟合系数aij使ε为最小值，等效于式(5)．

(5)

式(5)可化简为

(6)

式(6)可用矩阵形式表示为

AHAa=Ax

(7)

其中，AH为A的转置；a为拟合系数组成的向量，长度为(p+1)(q+1)；x为剪切波达到时刻T(x,y)组成的向量，长度为(2M+1)(2N+1)； 矩阵A为

(8)

A的大小为(2M+1)(2N+1) × (p+1)(q+1) 的矩阵，拟合系数a为

a=(AHA)-1Ax

(9)

确定二维SG滤波器尺寸后，可预先算出(AHA)-1A， 然后与感兴趣内的剪切波达到时刻T(x,y)做二维卷积，就能计得感兴趣内的横向差分.

2.3 弹性图像质量评估指标

对用二维SG差分器和线性回归分别估计出的剪切波弹性图，取4 cm×4 cm区域进行定量比较. 差分计算出的剪切波弹性图未经任何后处理操作. 对剪切波弹性图像进行质量评估所用的指标为SNR和CNR[20]. 信噪比SNR定义为

SNR=E/σ

(10)

其中，E和σ分别为所选区域内的剪切波速度的均值和标准差. 除了信噪比之外， CNR也被广泛用于弹性成像中，其定义如下

(11)

其中，EI和σI分别为所选的内含物内的剪切波速度的均值和标准差；EB和σB分别为所选背景内的剪切波速度均值和标准差.

3 实验结果

图2 剪切波在不同时刻的位置以及不同位置点的轴向运动速度随时间的变化Fig.2 (Color online) The positions of shear wave at different time and the tissue motion speed change with time at different positions

图2展示了剪切波激励后，用平面波检测组织内轴向运动速度的结果，这里轴向指探头发射声场的方向. 图2(a)为以第1次检测脉冲发射为起始时刻，第0.664 ms时计算出的探头下组织的轴向运动速度，可以看到剪切波的波前. 图2(b)为第2.075 ms时发射平面波检测脉冲，计算出的探头下组织的轴向运动速度. 比较图2(a)和(b)可见剪切波激励后向左右两边传播的趋势. 检测剪切波波前到达每个位置的时刻，就能计算出剪切波在每个位置的传播速度.

图2(c)表示图2(a)内的4个点A、B、C和D的轴向速度随时间的变化，其速度是经过方向滤波去除反射波后的结果. 这4个点所在的深度均为25 mm，横向位置分别是-2、0、2和4 mm. 剪切波右边的波前会先到达A点，然后依次到达B、C和D点. 剪切波在传播的过程中会不断衰减，如图2(c)，点A的运动速度幅度最大，然后依次是B、C和D点. 检测脉冲发射的时间间隔为83 ms，所以，图2(c)中曲线上的点的间隔也为83 ms. 组织内某点的速度达到最大值的时刻不一定位于检测脉冲的发射时刻，更多的可能是位于两次检测脉冲的发射时刻之间. 需要用插值上采样，使计算出的剪切波到达时刻的精度尽可能高. 以每一个位置的速度峰值时刻作为剪切波到达该位置的时刻，通过感兴趣区域内所有点的速度-时间曲线就可以计算出剪切波达到每个点的具体时刻.

图3(a)表示感兴趣区域内所有点的剪切波到达时刻. 由于剪切波在感兴趣区域的左侧进行激励，剪切波的传播是从感兴趣区域的左边到右边，感兴趣区域内的剪切波到达时刻从左到右的值会增大. 图3(b)和(c)分别表示用2.2和4.4 mm长度的线性回归计算出的感兴趣区域内的剪切波传播速度. 线性回归所使用的长度越短，差分计算出的剪切波速度的噪声越大，如图3(b)所示，在内含物和背景区域内都存在着较大的噪声. 线性回归的长度较长时，能起到平滑作用，明显减少噪声.

图3(d)和(e)分别表示用2.2和4.4 mm长度的二维SG数字差分器计算出的感兴趣区域内的剪切波传播速度. 同样随着差分器长度的增加，计算出的剪切波速度的噪声越小. 比较图3(b)和(d)，图3(c)和(e)，二维SG数字差分器计算出的剪切波速度分布图比线性回归的要更加光滑，噪声更小，尤其是在背景上对噪声去除更明显. 同时，二维SG滤波器计算出的剪切波传播速度图上，内含物的边界更加清晰和光滑，有利于在实际诊断中分割出异物的位置.

图3 感兴趣区域内剪切波到达的时刻和不同长度的二维SG差分器和线性回归计算出的剪切波速度 Fig.3 (Color online) The shear wave arrival time in the interesting region and the shear wave speed obtained by 2D SG differentiator and linear regression with different lengths

在内含物的中心深度29 mm的剪切波速度分布如图4(a)和(b)所示. 图4(a)对应的是2.2 mm长度的二维SG差分器和线性回归计算出的剪切波速度，二维SG差分器计算出的值在内含物区域内上下抖动更小. 在内含物与背景的边界的过渡带斜率越高，剪切波弹性图内的内含物越明显，图像的横向分辨率越高. 二维SG差分器和线性回归的边界过渡带斜率一致，说明这两种方法得到的剪切波弹性图的分辨率相同. 图4(b)为长度为4.4 mm差分的结果，二维SG差分器和线性回归在内含物的边界过渡带斜率也几乎相同. 比较图4(a)和图4(b)可见，计算差分时所用的长度越短，内含物与背景的过渡带斜率越高，对应的剪切波弹性图的分辨率越高. 计算差分时所用的长度越长，计算出的内含物内的剪切波速度的值方差越小，对应的剪切波弹性图的信噪比越高.

图4 用不同长度的二维SG差分器和线性回归计算出的内含物中心处的剪切波速度分布Fig.4 (Color online) The shear wave speed profile of inclusion obtained by 2D SG differentiator and linear regression

表1比较了用不同长度的数据进行差分所得到的剪切波速度图像的信噪比和对比度噪声. 二维SG差分器得到的剪切波弹性图像比线性回归法在对比度噪声上要高5.1～6.5 dB，在内含物的信噪比上要高2.9～3.1 dB，在背景的信噪比上要高0.5～2.6 dB. 二维SG差分器得到的剪切波弹性图像质量要优于线性回归法得到的剪切波弹性图像. 随着差分器长度的增加，剪切波弹性图的信噪比和对比度噪声都会增加，但付出的代价是图像分辨率降低.

表1 不同差分长度下，二维SG差分器和线性回归计算出的剪切波速度图的信噪比和对比度噪声Table 1 The comparsion of CNR and SNR of shear waveelastogram obtained by 2D SG differentiator andlinear regression

结 语

本研究基于超声剪切波弹性成像的原理，将二维SG滤波器用于计算剪切波的传播速度中，提高了剪切波弹性图的图像质量. 在超声弹性仿体上进行了实验，对比了本研究方法和传统剪切波弹性成像方法的图像质量，以信噪比和对比度噪声作为评估指标. 与现有剪切波速度估计所用的线性回归相比，本研究提出的方法能将剪切波弹性图背景的信噪比提高0.5～2.6 dB，内含物的信噪比提高2.9～3.1 dB，对比度噪声提高5.1～6.5 dB，且不降低剪切波弹性图的横向分辨率. 通常2D SG差分器的长度越长，剪切波弹性图的信噪比和对比度会提高，同时分辨率会降低. 在实际中，需要根据病灶的尺寸选择最佳长度的2D SG差分器. 后续将着手将该方法用于临床的乳腺剪切波弹性成像中，期望能提高乳腺肿瘤诊断的准确性.

引文：潘晓畅，冯乃章，陈思平. 二维Savitzky-Golay数字差分器剪切波弹性成像方法[J]. 深圳大学学报理工版，2018，35(1)：22-28.

/

[1] GARRA B S. Elastography: history, principles, and technique comparison[J].Abdominal Imaging,2015,40(4):680-697.

[2] PALMERI M L, WANG M H, DAHL J J, et al. Quantifying hepatic shear modulus in vivo using acoustic radiation force[J].Ultrasound in Medicine & Biology,2008,34(4):546-558.

[3] LEE E J, JUNG H K, KO K H, et al. Diagnostic performances of shear wave elastography: which parameter to use in differential diagnosis of solid breast masses?[J].European Radiology,2013,23(7):1803-1811.

[4] FERRAIOLI G, TINELLI C, DAL BELLO B, et al. Accuracy of real-time shear wave elastography for assessing liver fibrosis in chronic hepatitis C: a pilot study[J]. Hepatology, 2012, 56(6):2125-2133.

[5] SARVAZYAN A P, RUDENKO O V, SWANSON S D, et al. Shear wave elasticity imaging: a new ultrasonic technology of medical diagnostics[J]. Ultrasound in Medicine & Biology,1998,24(9):1419-1435.

[6] BERCOFF J, TANTER M, FINK M.Supersonic shear imaging: a new technique for soft tissue elasticity mapping[J].IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control,2004,51(4):396-409.

[7] WANG M H, PALMERI M L, ROTEMBERG V M, et al.Improving the robustness of time-of-flight based shear wave speed reconstruction methods using RANSAC in human liver in vivo[J].Ultrasound in Medicine & Biology,2010,36(5):802-813.

[8] ROUZE N C, WANG M H, PALMERI M L,et al.Robust estimation of time-of-flight shear wave speed using a radon sum transformation[J]. IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control,2010,57(12):2662-2670.

[9] ROUZE N C, WANG M H, PALMERI M L,et al. Parameters affecting the resolution and accuracy of 2-D quantitative shear wave images[J]. IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control,2012,59(8):1729-1740.

[10] AMADOR CARRASCAL C, CHEN Shigao, MANDUCA A,et al.Improved shear wave group velocity estimation method based on spatiotemporal peak and thresholding motion search[J]. IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control,2017,64(4):660-668.

[11] LIPMAN S L, ROUZE N C, PALMERI M L,et al. Evaluating the improvement in shear wave speed image quality using multidimensional directional filters in the presence of reflection artifacts[J].IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control,2016,63(8):1049-1063.

[12] YOON H, AGLYAMOV S R, EMELIANOV S Y. Dual-phase transmit focusing for multi-angle compound shear-wave elasticity imaging[J].IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control,2017,64(10): 1439-1449.

[13] LUO Jianwen, YING Kui, HE Ping, et al. Properties of Savitzky—Golay digital differentiators[J].Digital Signal Processing,2005,15(2):122-136.

[14] LUO Jianwen, BAI Jing, HE Ping,et al. Axial strain calculation using a low-pass digital differentiator in ultrasound elastography[J]. IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control,2004,51(9):1119-1127.

[15] GOWRI B G, HARIHARAN V, THARA S, et al. 2D image data approximation using Savitzky Golay filter-smoothing and differencing[J]. IEEE. 2013, 7903: 365-371.

[16] DENG Yufeng, ROUZE N C, PALMERI M L, et al. Ultrasonic shear wave elasticity imaging sequencing and data processing using a verasonics research scanner[J]. IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control,2017,64(1):164-176.

[17] DENG Yufeng, ROUZE N C, PALMERI M L,et al. On system-dependent sources of uncertainty and bias in ultrasonic quantitative shear-wave imaging[J]. IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control,2016,63(3):381-393.

[18] KASAI C, NAMEKAWA K, KOYANO A, et al. Real-time two-dimensional blood-flow imaging using an auto-correlation technique[J]. IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control,1985,32(3):458-464.

[19] DEFFIEUX T, GENNISSON J, BERCOFF J,et al. On the effects of reflected waves in transient shear wave elastography[J]. IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control,2011,58(10):2032-2035.

[20] VARGHESE T, OPHIR J. An analysis of elastographic contrast-to-noise ratio[J]. Ultrasound in Medicine & Biology,1998,24(6):915-924.