蔡 雪，余锦华，汪源源，马友能

(复旦大学电子工程系，上海 200433)

0 引 言

超声成像经历半个多世纪的发展，已广泛应用于临床诊断。早在 1980年 Muir和 Carstensen就提出声波在生物组织中的传播过程存在不可忽视的非线性现象，Starritt等更于1985年通过医学仪器观察到超声波在肌肉中的非线性传播[1]。单一频率的发射波在组织中传播时，随着传播深度的变化，其传播速度在变化。由于波峰的声速大于波谷的声速，所以波形会发生畸变，产生与初始频率成倍频关系的频率分量，这些新增的频率分量叫做“谐波频率”分量[2]。在声波继续传播的过程中，将产生无限多的谐波。最后，发射的正弦波变为锯齿波。谐波产生和累积的过程就是声波传播的非线性效应。

如果声的传播过程是等熵的，则介质的物态方程可以表示成泰勒级数的形式：其中：P为声压；ρ为密度；P0、ρ0分别为介质的静态压强和静态密度。将线性项系数A、二次项系数B定义为A=(ρ∂P/∂ρ)0，B= (ρ2∂2P/∂ρ2)0，那么B和A的比值就是非线性参量B/A。它是超声波在介质传播过程中非线性效应大小的度量。与声速、声阻抗、声衰减等线性参量相比，非线性参量更能反映生物组织的结构特征以及病变状态的动态特性。作为超声医学诊断中的一个新参量，B/A近年来在超声医学诊断与生物声学中出现了可喜的应用前景。

非线性参量B/A的测量方法主要分为两大类，一类是热力学方法，另一类是有限振幅声波法[3]。虽然热力学方法测量样品的非线性参量B/A具有一定的精度，但是要求有加压装置，所以不适用于临床。与之相反，有限振幅声波法由于实验操作简单，成为了研究热点。迄今为止，研究者提出了不少基于单频发射波或多频发射波的有限振幅声波法。这些方法结合传统的图像重建方法，形成了不同的B/A成像技术。

Nakagawa等[4]于1986年提出基于参量阵和二次谐波的非线性成像方法。基于参量阵的非线性成像方法是根据参量阵理论，即两列不同频率初始波在介质中传播方向相同时，会产生差频波，且这一差频波是非线性参量沿声波传播路径的线积分。该方法在获取差频波幅值后利用传统的 CT成像算法重建出非线性声参量的层析图像。基于二次谐波的非线性成像方法以 Burges非线性传播方程为基础，即声波在介质中传播所产生的二次谐波可以看作是非线性声参量沿传播路径的线积分。通过测量声波在传播距离上产生的二次谐波幅值，以其为投影数据，利用传统的CT成像算法重建出非线性声参量的层析图像。

Fujii和 Akiyama等[7]提出了另一种测量B/A参量的方法。该方法连续发射频率分别为f0和2f0的发射信号，从发射频率为f0的接收信号中提取谐波幅值，在发射信号为2f0的接收信号中提取基波幅值，将得到的谐波幅值和基波幅值的平方相比便可得到非线性参量B/A的值。此外，龚秀芬等[6]运用该方法得到了等深度脉冲回波 C扫描B/A成像。

2011年，Varray等[8,9]考虑介质中非线性参量的非均匀分布的情况，提出了直接法(Direct Method, DM)与比较法(Comparative Method, CM)的拓展形式：扩展直接法(Extended Direct Method,EDM)和扩展比较法(Extended Comparative Method,ECM)，以提高B/A测量的精度。此外，他们通过滤波从 B型谐波图像中滤除斑点噪声，得到与二阶声场相关的图，在图上选定参考介质区和待测介质区，运用ECM即可得到B/A参量成像图。

B/A参量成像结果与B/A测量技术有直接关系,而上述B/A参量成像所用到的有限振幅声波法可分为两种：直接法和比较法。直接法指向组织发射一定频率的有限振幅超声波，从接收信号中提取所需信号幅值，计算其与发射信号幅值平方的比来得到组织B/A值。Chavrier[10]和 Wallace[11]进行不同实验验证了传统的直接法测量均一组织的B/A参量误差在 5%以内，而 Varray等[8,9]在此基础上提出了适用于B/A分布不均的扩展直接法(EDM)。Fujii和 Akiyama等[7]提出的方法也可以归于直接法，只是这一方法用2f0发射源的接收信号的基波幅值代替f0发射源的幅值。比较法则主要测量声波在两种不同介质中传播时二次谐波的变化。其中一种介质的B/A值是已知的，被称为参考介质。同样为了使比较法适用于B/A参量分布不均的组织，Varray等[8,9]提出了扩展比较法(ECM)。

本文通过仿真与仿体实验比较有限振幅声波法测量非线性参量B/A的几种经典方法。仿真和仿体实验结果均表明：扩展的比较法(ECM)测量准确性高，运用该方法结合滤波反投影仿真实现B/A参量层析成像，成像结果能很好地描述简单组织中的B/A参量分布。

1 测量非线性参量的有限振幅声波法

1.1 直接法(DM法)

超声在组织中传播时，产生的二次谐波p2与组织的非线性系数β的关系为：

其中：ρ0和c0分别表示组织的密度和声速；x为声波传播的距离，f为发射源的频率，p0为发射源的振幅。由式(2)可知，发射一定频率的有限振幅波，当其在组织中传播一定距离后，从接收到的信号中提取出二次谐波幅值，便可得到组织的非线性系数。

为了使测量结果更加准确，在计算时需要补偿声波在组织中传播的衰减。此外，还需要测量不同传播距离下产生的二次谐波幅值，绘制p2(x)/xp02对x的曲线，并在该曲线上进行线性外推，得到x→0时p2(x) /xp02的值。计算β值时只需将此值代入：

2.1 出生后24 h细胞因子比较 从生后24 h临床监测动态指标来看，非BPD早产儿入院前外周血清IL-1β、IL-6、IL-8等细胞因子水平比BPD早产儿的整体水平要高，差异有统计学意义（P＜0.05），见表1。

其中，α1和α2分别对应基波和谐波的衰减系数。

1.2 扩展的直接法(EDM法)

传统的直接法是在假设组织的非线性系数是均一的前提下进行的。对于β非均匀分布的情况，式(3)只能得到整个组织的平均值。为了准确测得非均匀组织中非线性系数的值，Varray等[8,9]将传统的直接法进行了扩展，认为组织的非线性系数不再是常数，不同传播距离上的非线性系数是不同的，于是有：

对式(4)左右两边均作求导运算，则可以得到非线性系数分布不均时的计算表达式：

1.3 Akiyama和Fujii方法

Akiyama和 Fujii[7]利用超声回波信号测量组织的B/A值。发射频率为f0的正弦波并从接收信号中提取出超声波在组织中传播所积累的二次谐波幅值。随后发射频率为2f0的正弦波，从接收信号中提取其基波幅值。于是非线性系数的计算式为

其中：p0(f0)和p0(2f0)分别为两次发射波的初始振幅；p2(x)是频率为f0的发射波产生的二次谐波幅值；pII(x)是频率为 2f0的发射波经过组织后接收到的基波幅值。虽然该方法主要针对回波信号提出，但对于透射波同样适用。

1.4 比较法(CM法)

比较法是插入取代法的特殊形式，最初是由龚秀芬等提出[5]，其核心思想是：在选定B/A参量已知的参考介质后，测量同一频率的发射源在参考介质和待测介质中传播积累的谐波幅值，从而得到待测介质的B/A值。如果参考介质和待测介质的长度d相等，两种介质的声速和密度都已知，同时考虑衰减和衍射的影响对该方法进行改善，那么，待测介质的非线性系数可以表示为

其中：下标为 0和下标为i的参量β、ρ和c分别对应参考介质和待测介质的物理量；p2i和p20分别表示待测介质和参考介质的二次声压幅值；I1和I2是衰减相关项，I1=exp(−α1d)，I2=exp(−α2d/2)，其中α1和α2分别表示基波和谐波的衰减系数；D'和D''是两个无量纲项，表示透射系数，其表达式为：

1.5 扩展的比较法(ECM法)

Varray等[8,9]指出，如式(7)所示的比较法更适用于均匀介质。为了使得比较法在非均匀介质中也同样适用，他们对原先的比较法进行了扩展。在扩展的比较法(ECM)中，依然选取非线性参量均匀分布的介质作为参考介质，那么其二次谐波幅值(p20)为

由于待测介质的非线性参量和衰减系数都是非均匀分布的，可以用章东等提出的表达形式[12]来表示其二次谐波幅值(p2i)：

在式(9)、式(10)中，下标为 0和i的参量分别对应参考介质和待测介质的物理量。以二次谐波幅值的表达式为基础，可以得到非均匀待测介质的非线性系数：

其中，V和W是与参考介质和待测介质的衰减系数相关的两个系数。

2 非线性参量测量实验与结果

2.1 仿真实验与结果

目前，已有几种模型来模拟声波在非均匀介质里的非线性传播过程。Pinton等[13]于 2009年提出了非线性衰减的全波(Full-wave)模型，并运用时域有限差分法(Finite Difference Time Domain,FDTD)得到了其数值解。与目前应用最为广泛的KZK抛物线方程相比，Full-wave模型不受近轴假设的限制，其对远场超声波束的描述更为准确。Full-wave模型的数学形式为

其中：p表示声压；c0为声速；ρ0表示静态密度；δ为扩散系数；β(β=1+B/ 2A)是非线性系数。式(12)左边前两项表示波动方程的线性项，第三项表示超声的热耗散，第四项为非线性项，其余两项为密度变化和弛豫衰减项。

本文在模拟超声在组织中的非线性传播时，采用 FDTD求解该模型[13]。将声波传播区域划分为若干个小的网格点，每个网格点的大小为 0.1 mm×0.1 mm，所以空间步长Δx为 0.1 mm。Fullwave模型中的各偏微分项近似为差分项后，计算每一个时间步(时间步长 Δt= Δx/c0/2)上的声场，并通过数值迭代对声场进行更新。此外，为了将仿真区域局限在有限空间，在截断边界设置完全匹配层(Perfectly Matched Layer, PML)[14]。

按照组织的非线性参量分布均匀与否，我们进行了两组仿真实验。一组是建立三种均匀的生物组织(肾、肝和肌肉)的仿真模型。由于比较法和扩展比较法都需要选定参考介质，因此我们也建立了水的仿真模型。四种组织的声学参数的设定与相关文献的值是一致的，其中水、肝、肾和肌肉的B/A参量分别设为5.2、6.8、7.2和7.5。另一组模型的组织非线性参量是变化的，具体为

仿真中设置的探头为 65个阵元的线阵探头，阵元间隔为 2个网格点。当周期为 10、振幅为 1 MPa、频率为1 MHz(或2 MHz)的高斯调制正弦脉冲波在组织中传播时，每隔5 mm接收一组数据。声波在组织中传播10 cm后停止接收。对接收到的信号进行傅里叶变换。为了提取接收信号中二次谐波量(或基波)，将其送入一个中心频率为 1 MHz(或 2 MHz)、带宽为 0.6 MHz的四阶巴特沃兹带通滤波器进行滤波。

对于第一组组织模型，五种方法测得的组织非线性参量B/A的平均值如表1所示。第二组组织模型测量结果如图1所示。

表1 均一组织B/A测量仿真结果Table 1 Simulation results of B/A measurements for homogeneous tissues

图1 非均匀组织的B/A测量仿真结果Fig.1 Simulation results of a set of inhomogeneous B/A setting measurements for inhomogeneous tissues

从结果可知：对于均匀组织，DM、EDM、CM、ECM的测量结果都很理想。Akiyama和Fujii方法结果稍差，但是误差也在 10% 以内。然而对于非线性参量分布不均的情况，Akiyama和 Fujii方法随传播深度的增加，结果越加偏离正确值，无法很好描述组织非线性参数的变化。此外，传统的 DM也无法反映组织的非线性参量变化，只能得到整个组织B/A的平均值。尽管 CM在组织参量没发生变化时测量值与设定值几乎一致，但是当组织B/A发生变化后，该方法不再有效。EDM尽管能很好地识别组织B/A的变化趋势，但是其测量值偏离设定很多。而 ECM不仅能识别组织非线性参量的变化，而且其测量值也很接近设定值，只在后端有一定振荡。

由仿真实验可以得到结论：ECM是五种方法中最有效的测量法。

2.2 仿体实验与结果

仿体实验的操作平台是基于非线性超声设备Ritec RAM-5000搭建的。该系统拥有高能量脉冲输出端和相敏超外差接收器。可使声学介质的非线性研究更加简单和精确。实验框图如图2所示。

运用 Ritec RAM-5000系统脉冲信号输出端经50 Ω的负载电阻与发射探头连接，接收端探头经前置放大器后进入 Ritec RAM-5000系统输入端，该信号再由系统信号显示端输入至示波器显示。实验中选用的两个平面水浸探头分别为泛美V302-SU和 V304-SU，中心频率分别为 1 MHz和 2.25 MHz。

图2 实验装置图Fig.2 Experimental layout

仿体实验主要是针对健康的猪肉组织进行的，将三种组织都分为不同厚度的小样以便得到不同传播距离下的接收信号。与仿真实验一样，为了实现 CM和 ECM，假设水为参考介质。表2是测量中用到的组织的线性声学参量汇总。表3为不同方法对三种组织的测量结果。

表2 组织的线性声学参量Table 2 Linear acoustical parameters of tissues

表3 仿体实验B/A测量结果Table 3 Phantom experiment results for B/A measurement

由于每种组织的非线性参量都是分开测量的，每种组织都可看作是均匀组织。由于EDM和ECM可测得不同传播距离下组织的B/A值。为了避免近场区域多重反射引起的干扰(实验发现测量距离大于1.8cm时干扰减小)，将仿体切割为15个不同厚度{d|d=2cm: 0.3cm: 6.5cm}的小样。表3中EDM和ECM测量值是不同距离下测得值的平均。

从测量结果可得，虽然几种方法测得组织的非线性参量B/A值的误差(DM：8%，EDM：5%，A&F：6%，CM：4%，ECM：2%)都在误差允许范围内，而 DM和 CM对应的拓展方法都比这两种方法得到的测量结果要好，但 ECM测量得到的值与相关文献给出的值最接近。

3 基于ECM的B/A成像计算机仿真

3.1 基于ECM的B/A层析成像

由以上仿真和仿体实验，可以看出ECM能较好描述组织B/A的变化。而式(9)和式(10)之比可以表示为

将式(14)中的积分因子看成一个整体，整理可得：

其中式(15)左边可看作是β'(u)沿声传播路径x的线积分。如果在某一角度上等间距求得待测介质和参考介质的二次谐波声压之比，就可以利用CT的重建算法来进行非线性参量的层析成像[15]。

设投影数据p(u,θ)为

按图3的 CT扫描系统所示，在 0°～180°间等角度地旋转组织，并在每一角度下等间隔收集投影数据，扫描完成得到投影矩阵后便可通过滤波反投影算法得到β'(u)。

重建β'(u)的关键在于斜坡滤波器|w|的选择，本文选用 Sheep-Logan滤波器在频域对投影矩阵进行滤波，其频响曲线为

图3 CT扫描系统Fig.3 CT scanning system

其中d为等间隔扫描的间隔。对滤波后的投影数据在时域进行线性内插，得到β'(u)后对其进行衰减补偿便可得到B/A的层析成像。

3.2 成像实验及结果

对图4(a)和图5(a)所示的组织模型(4×4 cm2，即 400×400网格点)进行层析成像。图4(a)所示组织模型，中心为一个圆形非线性系数β为7的脂肪类组织，外层为水，非线性系数β为 3.5。图5(a)所示组织模型为了增加组织模型的复杂度，在图4(a)所示组织模型的中心圆内又增加了一个小圆，其非线性系数β为 4.8。按图3所示的扫描方式，发射探头和接收探头之间的距离不变，待测样品每旋转 5°时，发射探头和接收探头同时沿u轴，中心阵元每隔 10个网格点移动一次进行扫描，仿真所用的探头阵元数为 65。为了减少探头阵元之间的相互干扰对成像结果的影响，每个采样步长上只取中心阵元得到的数据作为采样数据，并从采样数据中提取出二次谐波幅值。采样组织的整个扫描过程可得到不同投影区域积累的二次谐波幅值(41×37的矩阵)。选定图4(a)和图5(a)中两线之间的区域为参考介质，则将得到的二次谐波幅值矩阵除以参考区域的二次幅值，并乘以参考区域的非线性系数与其余线性参数之比便可得到投影矩阵(41×37)。对投影数据进行滤波反投影重建后，便可得到该组织的B/A层析图像，结果如图4(b)和图5(b)所示。提取其中心一条线，可得到非线性参量的一维分布，如图4(c)和图5(c)所示。

图4 组织模型I的 B/A层析成像结果Fig.4 B/A tomography results for tissue model I

图5 组织模型II的B/A层析成像结果Fig.5 B/A tomography results for tissue mode II

从成像结果可以看出，基于 ECM的非线性成像结果与组织模型较吻合，能够很好地反应组织的非线性参量分布。

4 总结与展望

为了获得较好的B/A成像结果，本文比较了生物组织非线性B/A参量测量的几种经典算法，希望从中选择最优的方法进行成像。仿真和仿体实验均表明：ECM 法不但测量精确，而且其测量结果能很好地反映组织中B/A的变化。运用 ECM进行B/A成像的仿真模拟，结合ECM与滤波反投影得到的B/A层析成像结果与设定的组织模型吻合程度较高。

在本文的仿真和仿体实验中，忽略了组织中散射点对结果的影响。而在实际中，散射点对超声波在组织中传播是有影响的，以后对于B/A成像的研究也应考虑散射点的影响。

参考文献

[1] 余锦华, 汪源源. 医学超声成像的模拟研究[J]. 声学技术, 2011,30(1): 33-40.YU Jinghua, WANG Yuanuan. Simulation study of medical ultrasound imaging[J]. Technical Acoustics, 2011, 30(1): 33-40.

[2] Choudhry S, Gorman M, Charboneau W, et al. Comparison of tissue harmonic imaging with conventional US in abdominal disease[J]. Radiographics, 2000, 20(4): 1127-1135.

[3] BjØrnØ L. Forty years of nonlinear ultrasound[J]. Ultrasonics,2002, 40(1-8): 11-17.

[4] Nakagawa Y, Nakagawa M, Yoneyama M, et al. Nonlinear parameter imaging computed tomography by parametric acoustic array[C]// Proceeding of IEEE International Ultrasonics Symposium(IUS), 1984: 673-676.

[5] Gong X, Feng R, Zhu C, et al. Ultrasonic investigation of the nonlinearity parameter B/A in biological media[J]. J. Acoust. Soc.Am., 1984, 76(3): 949-950.

[6] Gong X, Zhang D, Liu J, et al. Study of acoustic nonlinearity parameter imaging methods in reflection mode for biological tissues[J]. J. Acoust. Soc. Am., 2004, 116(3): 1819-1825.

[7] Fujii Y, Taniguchi N, Akiyama I. A new system for in vivo assessment of the degree of nonlinear generation using the second harmonic component in echo signals[J]. Ultrasound Med. Biol., 2004,30(11): 1511-1516.

[8] Varray F, Pasovic M, Cachard C, et al. Acoustic nonlinearity parameter of tissue on echo mode: review and evaluation of the different approaches for B/A imaging[C]// Proceeding of IEEE International Ultrasonics Symposium (IUS), 2009: 41-44.

[9] Varray F, Basset O, Tortoli P, et al. Extensions of nonlinear B/A parameter imaging methods for echo mode[J]. IEEE Trans. Ultrasonics. Ferr., Freq. Control, 2011, 58(6): 1232-1244.

[10] Chavrier F, Lafon C, Birer A, et al. Determination of the nonlinear parameter by propagating and modeling finite amplitude plane waves[J]. J. Acoust. Soc. Am., 2006, 119(5): 2639-2644.

[11] Wallace K, Lloyd C, Holland MR, et al. Finite amplitude measurements of the nonlinear parameter B/A for liquid mixtures spanning a range relevant to tissue harmonic mode[J]. Ultrasound Med. Biol., 2007, 33(4): 620–629.

[12] Zhang D, Gong X, Ye S. Acoustic nonlinearity parameter tomography for biological specimens via measurements of the second harmonics[J]. J. Acoust. Soc. Am., 1996, 99(4): 2397-2402.

[13] Pinton G, Dahl J, Rosenzweig S, Trahey G E. A heterogeneous nonlinear attenuating full-wave model of ultrasound[J]. IEEE Trans. Ultrasonics. Ferr. Freq. Control, 2009, 56(3): 474-488.

[14] Li X. PML condition for the numerical simulation of acoustic wave[C]// Proceeding of International Conference on Computing Control and Industrial Engineering, 2010: 129-132.

[15] Zhang D, Gong X. Experimental investigation of the acoustic nonlinearity parameter tomography for excised pathological biological issues[J]. Ultrasound Med. Biol., 1999, 25(4): 593-599.