圆柱类构件超声相控阵聚焦模型

孔惠元 王 鉴 李仰军 韩 焱

(1 中北大学 山西省信息探测与处理重点实验室 太原 030051)

(2 中北大学 山西省现代无损检测工程技术研究中心 太原 030051)

0 引言

相控阵超声检测技术是近年来发展起来的新型超声无损检测技术，其基本原理是通过控制不同阵元的激发时间，使所有阵元在被测构件中产生的超声场“相长干涉”和“相消干涉”，从而实现合成波束在预定位置聚焦[1]，在不移动换能器的前提下实现对被测构件大范围检测[2]。与常规单阵元换能器相比，相控阵换超声检测具有声束可控、覆盖率大、焦区能量强等优势[3]。目前超声相控阵在平面构件无损检测中得到了广泛的应用。然而，圆柱类构件在工程中大量存在，如固体火箭发动机、钢管、传动轴架等，在利用相控阵换能器对此类构件进行检测时，受曲界面结构引起的入射波和回波时延的影响，扫描声束的波阵面产生弯曲。目前在该类检测应用中，往往采用迭代遍历的方式进行延迟时间的计算，从而实现波阵面的控制。然而此类算法需大量迭代离散点，计算速度慢，且计算精度受离散点取值影响，从而导致实时性较弱、检测效果差，影响了在该类构件检测中的应用。为了提高相控阵仪器针对此类结构检测时的效率，充分发挥相控阵仪器的优势，本文建立了一种具有通用性的柱类构件相控阵波束扫描模型，该模型基于换能器、耦合介质、圆柱类构件的声学、材料特性和几何关系，通过输入耦合介质和被检测构件的声速、曲率半径及换能器阵列与曲面的距离，即可计算出换能器各阵元的收发时间，利用该时间进行换能器收发时间控制即可实现柱类构件的相控阵成像。

目前针对单一介质以及双层平面结构的超声相控阵聚焦延时法则已经非常成熟，对圆柱类构件的聚焦方案研究较少，主要难点是曲面使超声波发生反射、折射等现象，已知聚焦点位置逆向求解声束路径时，声束在界面的入射点确定困难，从而无法确定延迟时间。当前对超声能量聚焦可以分为两个方向，一是物理聚焦，通过相控阵换能器控制各激励通道波形形状、时间延迟来实现声波能量在特定点的聚焦；二是后处理聚焦，通过后处理各换能器的检测信号以实现能量在指定位置处的虚拟聚焦[4]。后处理方法意味着聚焦声束不会物理地存在于被测构件中，而是通过用阵列数据的完整矩阵和相关成像算法合成其效果[5]。利用相控阵换能器控制波束在曲面构件内部聚焦已有部分学者进行了相关研究。徐娜等[6]将界面离散化，基于费马定律对每一个离散点与阵元和聚焦点之间的声程传播时间进行计算，将最短传播时间路径作为声线真实路径，该方法需要遍历界面离散点，计算量大；周正干等[7]基于折射定律对每一个入射角度进行遍历，寻找折射声线与聚焦点距离最近的声线路径进而确定声线传播时间，该方法同样需要进行大量迭代计算，计算效率低；甘勇[8]指出，上述两种求解方法的精度均取决于一个事先设定的增量∆，∆越小，解的精度越高，但求解速度越慢；高世凯等[9]将阵元垂直向下辐射声波与界面的交点视为折射点，进而确定声束路径与传播时间，该方法对折射点的确定可以理解为在曲界面上构建了多个虚拟源，只是一种假设，并不满足折射定律；Zheng[10]研究了柔性相控阵聚焦方法，但柔性阵列换能器制造工艺复杂且价格昂贵，不具备普遍性[11]。利用回波信号对曲面构件进行虚拟聚焦也有部分学者进行了相关研究。关山月等[12]基于费马定律对曲面结构全聚焦算法声束传播时间进行计算，该方法计算量同样较大；Camacho等[13]将虚拟源与合成孔径算法结合，通过阵元参数和近似几何模型求取曲界面上虚拟源坐标，将声源通过折射点到达聚焦点的路径替换为声源通过虚拟源到达聚焦点的路径，进而实现声传播时间的计算与图像重建，避免了复杂的折射点坐标计算，但要求工件曲率半径远大于阵元中心间距，否则求取虚拟源坐标误差较大。

本文针对工业中使用一维线性相控阵换能器检测圆柱类构件聚焦法则展开研究，基于换能器、耦合介质、圆柱类构件材料特性和几何关系以及声线模型和折射定律，建立了耦合介质及被检构件的声速、曲面曲率半径、阵列与曲面间的距离等关联的延迟时间聚焦控制模型。通过计算构件内部聚焦时阵元发射声束在界面入射点坐标，得到各阵元发射声波到达聚焦点时间得到延迟时间。经有限元仿真实验，验证本文聚焦方案在效率上的优势与可行性。

1 圆柱类构件扫描波束的影响因素分析

1.1 双层结构平界面相控阵成像检测波束形成原理

双层结构平界面聚焦声束路径计算示意图如图1 所示，首层厚度h，声速c1，第二层介质声速c2，且满足c1

图1 双层结构平界面聚焦声束路径计算示意图Fig.1 Schematic diagram of path calculation of the two-layer flat interface focused sound beam

设换能器阵元数为N，中心间距为p，阵元i发出声束以入射角α入射到界面点a(m,h)，然后以折射角β到达位于第二层结构内的聚焦点F(xf,yf)，则求解点a坐标过程如下。

设阵元i位置坐标(xi,0)，由已知条件得

入射角α满足

折射角β满足

入射角α、折射角β满足折射定律

将式(2)、式(3)代入式(4)得

通过解算式(5)即可求得a(m,h)，则声束从阵元i发射到达聚焦点F(xf,yf)的总传播时间TiF为

重复上述步骤，计算所有阵元到达聚焦点F对应声时构成集合{TF}，则第i个阵元的延迟时间可以表示为式(7)，通过调整各阵元的不同延迟时间实现声场在点F波束聚焦。

1.2 圆柱类构件相控阵成像检测的影响因素分析

目前工业中对圆柱类构件的检测通常采用水浸耦合或定制楔块的方式将声场能量传递到被测构件中进行检测[14]。采用相控阵成像检测时，对某一聚焦点，计算不同阵元发射声束经耦合液(或楔块)-构件界面的入射点及时间延时至关重要[15]。不同于双层结构平界面入射点的计算，声束在耦合介质-圆柱类构件界面的入射点横纵坐标都是变量，且在界面的每一个入射点法线方向与x轴方向不垂直，随着入射点的变化，法线方向发生偏转，基于折射定律求解入射点时，阵元、入射点、聚焦点坐标和入射角、折射角难以通过直角三角形建立几何关系。当前许多学者求取折射点坐标方法可以分为基于费马定律(或折射定律)的遍历法和在界面构建虚拟源法，而基于折射定律和构件几何关系直接求取折射点坐标，进而确定延迟时间暂无相关研究。

2 圆柱类构件扫描控制波束时间

2.1 聚焦声束路径及延迟时间计算

超声波经耦合介质入射圆柱类构件时，界面形状(凸面或凹面)、耦合介质和被测构件声速相对大小都会影响透射声束的传播特性，使声束在圆柱类构件内部自动发生扩散或者聚焦现象，本节以液体为耦合剂、超声波入射凸面构件为例推导聚焦声束路径计算方法及延迟时间、声束入射临界条件。

圆柱类构件聚焦声束路径计算示意图如图2 所示，被测构件表面为圆弧，曲率为k，曲率半径为r=1/k，圆心为O，声速c2，相控阵换能器位于被测构件正上方，以阵列换能器中心为原点、水平向右为x轴且为正方向、垂直向下为y轴且为正方向建立直角坐标系，换能器中心位置与被测构件中心轴线处界面垂直距离为h，中间位置充斥着液体耦合剂，耦合剂声速c1，且满足c1

图2 曲面构件聚焦声束路径计算示意图Fig.2 Schematic diagram of focused beam path calculation for curved component

设换能器阵元数为N、中心间距为p，则聚焦点F(xf,yf)的聚焦声束路径及延迟时间可通过如图2所示的几何关系求得。

(1) 基于被检构件内聚焦点位置的曲面入射点求解

设阵元i发出声束以入射角α入射到界面点a(m,n)，然后以折射角β到达位于构件内的聚焦点F(xf,yf)。尽管从阵元i到a再到F有无数条路径，但是声波只沿遵循折射定律的直线传播，此时声线从阵元i到聚焦点F的路径是唯一的，即入射角α是唯一的，这决定了点a的坐标[16]。

设阵元i位置坐标(xi,0)，xi满足式(1)，由图2中几何关系可知入射角α为向量(-m,h+r-n)和向量(xi-m,-n)夹角的补角，即

入射角α、折射角β满足折射定律

将sinα、sinβ向量坐标表达式代入式(12)，并与点a所在曲率圆上方程联立得方程组

通过解算式(13)即可求得a(m,n)。

(2) 解得入射点a的坐标后，声束在耦合剂中声程l1与传播时间T1为

(3) 被测构件中声程l2与传播时间T2为

(4) 阵元i发出声波到达聚集点F传播时间TiF为

重复上述步骤，计算所有阵元到达聚焦点F对应声时构成集合{TF}，则第i个阵元的延迟时间可以表示为

通过式(19)即可解出各阵元发射声束在点F聚焦的延迟时间，同样对于任一位置聚焦点坐标代入上述步骤都可求出各阵元的延时，通过调整各阵元的不同延迟时间实现声场在预定点能量增强。本节以液体耦合剂为例进行声束路径计算及延迟时间的推导，对于楔块耦合下延迟时间的计算同样适用，此外，本文以晶片阵列中心与圆柱类构件圆心处于同一垂直方向为例进行推导，而在推导过程中并不严格要求此两点处于同一垂直方向，只需获取各阵元坐标即可实现声延迟时间计算、能量聚焦，因此阵元中心与构件圆心存在偏置时也适用。

2.2 声束入射临界条件

阵元发射声线与构件左右两侧会有一点相切，此时入射角为90◦，两切线范围外的入射声束不会进入构件，此外如图2所示，入射角大于第一临界角θ时折射声束无法进入构件，设此时入射声束与界面两侧交点为p、q，则有效声束入射点位于p与q之间(这里不考虑波束角大小和波型转换现象，假设每个阵元发射声束可以覆盖整个构件，以纵波声线路径进行研究，不考虑横波影响)。设第i个阵元位置坐标(xi,0)，xi满足式(1)，设临界点p位置坐标(xp,yp)，点q坐标(xq,yq)，由折射定律有

同时入射角θ为向量(-xp,h+r-yp)和向量(xi-xp,-yp)夹角的补角，即入射角θ满足

将sinθ向量坐标表达式代入式(20)，并与点p所在曲率圆上方程联立得方程组

通过解算式(23)即可求得p(xp,yp)，在yp≤h+r范围内临界点p、q坐标同时满足式(23)方程组，因此可以解出两组实根，即为临界点p、点q坐标。

2.3 入射点方程组求解

在上述延迟时间的计算中，重点在于对非线性方程组式(13)的求解，目前求解非线性方程组的方法有很多，如直接降维法、最小二乘法、牛顿迭代法等，本文使用数据处理软件中相关工具包对方程组进行编程求解，算法思路及核心函数如图3所示。

图3 方程组求解算法流程图Fig.3 Flow of equation system solving algorithm

fsolve 函数是通过选择不同的参数来选取不同的算法，如牛顿迭代法或Levenberg-Marquardt 算法，再根据初始值和方程组的导数信息逐步接近方程组的解，本文通过选择Levenberg-Marquardt 算法寻找未知变量的值。Levenberg-Marquardt 算法是使用最广泛的非线性最小二乘法，它的关键是用模型函数对待估参数在其领域内做线性近似，忽略掉二阶以上的导数项，从而转化为线性最小二乘问题，它具有收敛速度快的优点。此外使用fsolve函数时提供一个合适的初始值非常重要，初始值选择不合理可能导致算法无法收敛。初始值的确定与阵元坐标、聚焦点坐标、检测对象参数、坐标系建立方式等有关，结合实际相控阵换能器各阵元中心间距较小(本文选用中心间距0.6 mm相控阵换能器)，则不同阵元对于初始值的选取影响较小，且本文以换能器中心建立坐标系，对于左右两侧的阵元，入射点分别位于构件中心左右两侧，因此本文选择待检测圆柱面上超声换能器阵列中心对应位置(0,5)作为初始值。

基于费马定律的迭代算法计算延迟时间时，入射点的坐标精度取决于界面的离散精度，离散精度越小，计算精度越高，但计算时间越长。分析本文的推导过程与方程组的求解过程可知，本文方法避免了费马定律遍历算法中界面离散精度的影响。然而在方程组求解过程中，所选算法不同、初始值选取不同也会影响方程组求解的精度与效率。实际检测过程中，受阵元坐标测量误差、换能器频率存在微小偏差导致近场区计算误差、构件声速与理论值存在偏差等因素影响，也会导致声束聚焦点坐标实际值与理论值存在微小偏差。

3 仿真验证与结果分析

3.1 仿真条件

为了验证本文方法在计算效率上的优势与有效性，本文利用数据处理软件进行数据计算，利用有限元仿真软件对如图4 所示的圆柱类构件进行仿真。模型结构的长30 mm、高20 mm，上半部分为耦合介质，下半部分为曲率半径r=50 mm的半圆，材料为钢，为节省计算时间，只截取36◦圆心角对应弧长作为耦合介质与钢的交界面。分别对耦合介质为水和有机玻璃楔块进行验证，材料参数如表1 所示，根据各向同性固体材料中纵波声速计算公式

表1 仿真使用耦合介质与钢的参数Table 1 Simulate parameters using coupling media with steel

图4 仿真模型结构图Fig.4 Simulation model structure diagram

可得纵波声速，式(24)中E为杨氏模量，ρ为密度，σ为泊松比。相控阵阵列为工程中使用比较普遍的16阵元相控阵换能器，中心频率5 MHz，阵元编号从左向右依次递增，阵列参数为阵元宽度d=0.5 mm，中心间距p=0.6 mm，阵列与钢的距离h=5 mm(指阵列中心与构件中心轴线处界面垂直距离)。

仿真中网格划分大小会直接影响到计算的精度，网格划分越小，计算结果越精确，但计算时间较长，对计算机性能要求较高，较大的网格则会导致计算结果不精确。结合本文研究内容，网格划分满足纵波波长包含12个网格即可，采用网格类型为自由三角形网格。由于超声波在两种介质中声速不同，网格划分需要分开计算，最大网格大小为

其中，fc为信号中心频率，本文设为5 MHz。时间步进的选取同样会影响计算精度，仿真中最大时间步进要求小于单个网格中超声波传播时间，本文设为T0/20，T0为激励信号周期，满足仿真要求。为减小模型中其他边界反射波对聚焦声场的影响，将耦合介质左右边界、钢中下边界设为吸收层。

由于声场只有在近场才能有效聚焦，相控阵线性阵列换能器近场大小L为

式(26)中，α为阵元宽高比系数，D为阵列孔径大小。将阵列换能器参数代入式(26)可得近场区为16 mm，因此将聚焦点设于近场区附近。

3.2 水浸耦合仿真验证

首先将耦合介质设为水，聚焦点设于被测构件中间轴线处，焦点深度距换能器下表面10 mm。根据本文方程组求解思路进行编程，同时为验证本文求解模型计算效率的优势，对本文所提方法与费马定律分别进行编程。费马定律的原理为设定某聚焦点，针对某一阵元，将该阵元与聚焦点横坐标范围内的界面离散化，依次计算局部界面各离散点与阵元和聚焦点之间距离之和并计算对应声线传播时间，所有声线路径中传播时间最小值即为声线真实传播时间，对应路径为真实传播路径，对比结果如表2所示。

表2 两种算法计算结果对比Table 2 Comparison of the calculation results of the two algorithms

由表2 可知，本文方法在计算次数上远远小于费马定律算法，且本文方法在计算效率上提升显著，对所有阵元进行计算所得相关参数如表3所示。

表3 各阵元编号、水中声程、被测构件中声程和延迟时间Table 3 The number of each array,the underwater sound path,the sound path in the measured component and the delay time

作为对比，绘制了没有施加延时的瞬态声场，仿真超声波传播过程中的声场分布，仿真结果如图5 所示。图5(a)表示所有阵元同时激发超声波，图5(b)表示随着深度增加波阵面在钢内逐渐扩散。

图5 阵元无延时瞬态声场分布Fig.5 The transient sound field distribution of the array without delay

将表3 中的延时依次加在相控阵换能器各阵元，使声波同时到达聚焦点位置，仿真结果如图6 所示。图6(a)表示换能器各阵元按照表3 的延时依次激发超声波，图6(b)表示在5.2 µs 左右在预设聚焦点位置能量增强实现聚焦。绘制聚焦声场钢内中心轴线处声压分布如图7所示(数据经过归一化处理)，由图7 可见，沿深度方向，声压逐渐增大，在10 mm 深度左右所有阵元发射声波相位一致能量叠加实现聚焦，声压达到最大值，而后随着深度增大阵元辐射声场再次发散，声压幅值逐渐减小，仿真和理论基本吻合。

图6 中心轴线处聚焦瞬态声场分布Fig.6 Focus transient sound field distribution at the central axis

图7 钢内中心轴线声压分布Fig.7 Sound pressure distribution in the central axis of steel

将聚焦点设于换能器中心位置右侧5 mm、深度10 mm (偏转加聚焦)，根据本文原理计算得各阵元激发延迟时间，同时对本文方法与费马定律算法分别进行编程，对比结果如表4所示。

表4 两种算法计算结果对比Table 4 Comparison of the calculation results of the two algorithms

由表4 可知本文方法迭代次数较少，且计算效率有效提升，对各阵元进行计算所得参数如表5所示。

表5 各阵元编号、水中声程、被测构件中声程和延迟时间Table 5 The number of each array,the underwater sound path,the sound path in the measured component and the delay time

将表5 中的延时数据依次加到相控阵换能器各阵元，使声波在预定点聚焦，仿真结果如图8 所示。图8(a)表示所有换能器各阵元按照表5的延时依次激发超声波，图8(b)表示在5.5 µs 左右在预设聚焦点位置能量增强实现聚焦。

图8 偏转聚焦瞬态声场声布Fig.8 Deflection focusing transient sound field sound distribution

由焦点坐标可得聚焦点位于换能器中心位置右侧向深度方向偏转63◦，绘制该方向声压分布，如图9 所示(数据经过归一化处理)，由图9 可知随着深度增大，声压逐渐增大，在10 mm 深度左右所有阵元发射的超声波相位一致能量叠加实现聚焦，声压达到最大值，之后随着深度增加声场再次发散声压逐渐减小，仿真结果与理论基本吻合。

图9 偏转63°方向聚焦声场声压分布Fig.9 Deflection 63°direction focused sound field sound pressure distribution

此外计算了其他位置点聚焦的延时时间并进行声场仿真验证，与上述聚焦结果类似，声场可以在指定位置实现聚焦。

3.3 楔块耦合仿真验证

将耦合剂换为有机玻璃制楔块，由于偏转聚焦情况更为普遍，对于只聚焦的情况这里不再阐述。将聚焦点设于换能器中心位置右侧4 mm、深度10 mm，根据上述原理得各阵元延迟时间，同时对本文方法与费马定律算法分别进行编程，对比结果如表6所示。

表6 两种算法计算结果对比Table 6 Comparison of the calculation results of the two algorithms

由表6 可知本文方法迭代次数较少，且计算效率有效提高，对各阵元进行计算所得参数如表7所示。

表7 各阵元编号、有机玻璃中声程、被测构件中声程和延迟时间Table 7 The number of each array,the plexiglass middle sound path,the measured component middle sound path and the delay time

将表7中的延时数据加到换能器各阵元仿真声场分布，如图10 所示。图10(a)所有换能器各阵元按照表7 的延时依次激发超声波，图10(b)表示在3.9 µs 左右所有阵元辐射超声波在预设焦点实现聚焦。

由焦点坐标可得聚焦点位于换能器中心位置右侧向深度方向偏转68◦，绘制该方向声压分布，如图11所示(数据经过归一化处理)，在10 mm深度左右声压值达到最大，之后随着深度增加声压减小，与预期结果基本吻合。此外计算了其他聚焦点的延时法则并进行声场仿真验证，与上述聚焦结果类似，声束都可以在指定位置实现聚焦。

图11 偏转68°方向聚焦声场声压分布Fig.11 Deflection direction 68° focused sound field sound pressure distribution

4 结论

针对工业中相控阵超声检测圆柱类结构时延迟时间计算量大的问题，利用解析法对声束聚焦延迟控制时间进行了理论推导并进行仿真验证，建立了基于圆柱类构件特征和耦合介质特性的超声相控阵扫描成像的聚焦模型。基于换能器、耦合介质、曲面的声学、材料特性和几何关系，利用声线模型和折射定律建立了耦合介质及被检构件的声速、构件表面曲率半径、阵列与曲面间的距离等关联的延迟时间聚焦模型。利用该模型计算出延迟时间，对各阵元激发超声波时间进行控制，进而实现扫描声束的聚焦。本文以水和有机玻璃制楔块作为耦合介质，对钢曲面的相控阵声束聚焦进行了仿真实验，结果表明本文方案在计算效率上提升显著，同时声束可以较理想地在指定位置实现聚焦，验证了本模型在计算效率上的优势与有效性。