CN102608205A

CN102608205A - Multilayer-object nondestructive testing ultrasonic imaging method based on variable wave speed phase shifting

Info

Publication number: CN102608205A
Application number: CN2012100461511A
Authority: CN
Inventors: 秦开怀; 杨春; 孙丰
Original assignee: Tsinghua University
Current assignee: Tsinghua University
Priority date: 2012-02-24
Filing date: 2012-02-24
Publication date: 2012-07-25
Anticipated expiration: 2032-02-24
Also published as: CN102608205B

Abstract

The ultrasonic imaging method for non-destructive testing of multi-layer objects based on variable wave velocity phase migration belongs to the technical field of rapid imaging of the internal structure of multi-layer heterogeneous objects, and is characterized in that it is based on using an ultrasonic transducer to make equidistant distances in the length direction of multi-layer heterogeneous objects. Under the condition of mobile detection, in the homogeneous area, the constant wave velocity phase shift ultrasonic imaging technology and fast Fourier transform technology are used to calculate the images of each line in the homogeneous area, and in the heterogeneous area, a bounding box surrounding the boundary of the medium is established , for each heterogeneous region above and below the boundary line of the medium in the bounding box, the variable wave velocity phase shift method and the N-times single-output inverse fast Fourier transform method are used to calculate the images of each row in the heterogeneous region. The invention is suitable for rapid ultrasonic imaging of multi-layer heterogeneous objects with heterogeneous media in both the depth direction and the horizontal direction, and adopts a separate treatment method characterized by variable wave velocity near the interface of the media, with high precision.

Description

Ultrasonic imaging method for non-destructive testing of multi-layer objects based on variable wave velocity phase shift

技术领域 technical field

本发明涉及超声波无损检测技术、合成孔径聚焦成像技术和相位迁移技术，实现对复杂的非规则分层物体的内部结构快速成像。The invention relates to ultrasonic non-destructive testing technology, synthetic aperture focusing imaging technology and phase shifting technology, and realizes rapid imaging of the internal structure of complex irregularly layered objects.

背景技术 Background technique

目前，在超声波无损检测领域，对于含有分界面的分层物体的超声成像(包含接触式超声成像和液浸法超声成像，若将耦合剂层视为被测物体的一部分，则液浸法超声成像也可视为分层物体的接触式超声成像，因此本发明中将两者统一称为分层物体的超声成像)，当前主要有两种方案：一是采用时域合成孔径聚焦超声成像技术与射线跟踪(Ray Tracing)技术相结合的方法，二是采用基于相位迁移(Phase Shift Migration)技术和频域合成孔径聚焦超声成像技术的相位迁移超声成像方法。At present, in the field of ultrasonic nondestructive testing, for ultrasonic imaging of layered objects containing interfaces (including contact ultrasonic imaging and liquid immersion ultrasonic imaging, if the couplant layer is regarded as a part of the object to be tested, the liquid immersion ultrasonic imaging Imaging can also be regarded as contact ultrasonic imaging of layered objects, so in the present invention, the two are collectively referred to as ultrasonic imaging of layered objects), currently there are two main schemes: one is to use time-domain synthetic aperture focused The method combined with Ray Tracing technology, the second is the phase shift ultrasound imaging method based on phase shift (Phase Shift Migration) technology and frequency domain synthetic aperture focused ultrasound imaging technology.

合成孔径聚焦超声成像技术(SAFT)的基本思想是利用一个单一的换能器来模拟孔径阵列，对被测物体进行有序的扫描，并采用延时叠加(DAS)方法(时间延迟或相位延迟)对扫描得到的脉冲回波信号进行聚焦成像。SAFT工作模型如图1所示，其中，Z向为深度方向，X向为换能器扫描方向，当换能器位于u_n处时，目标反射物(x′，z′)至换能器的距离为r_n。SAFT技术有时域和频域之分：The basic idea of the Synthetic Aperture Focused Ultrasound Imaging Technology (SAFT) is to use a single transducer to simulate the aperture array, scan the measured object in an orderly manner, and use the delay superposition (DAS) method (time delay or phase delay ) focus imaging on the pulse-echo signal obtained by scanning. The working model of SAFT is shown in Figure 1, where the Z direction is the depth direction, and the X direction is the transducer scanning direction. When the transducer is located at u _n , the target reflector (x′, z′) to the transducer The distance is r _n . SAFT technology is divided into time domain and frequency domain:

时域SAFT基于DAS原理与多点动态聚焦技术，对目标成像区域不同深度上的聚焦点计算不同的时延。如图1，为了在目标反射物处(x′，z′)聚焦，时域SAFT技术将换能器在合成孔径长度L内的每一个扫描位置处获得的回波信号进行延时叠加处理：设s(u_n，t)为换能器在u_n处接收到的回波信号，t为采样时刻，u_n处关于目标反射物(x′，z′)的延时为Time-domain SAFT is based on the DAS principle and multi-point dynamic focusing technology, and calculates different time delays for focal points at different depths in the target imaging area. As shown in Figure 1, in order to focus on the target reflector (x′, z′), the time-domain SAFT technology performs delay superposition processing on the echo signals obtained by the transducer at each scanning position within the synthetic aperture length L: Let s(u _n , t) be the echo signal received by the transducer at u _n , t is the sampling moment, and the delay at u _n about the target reflector (x′, z′) is

${τ τ}_{n no} = = \frac{22}{c c} (({r r}_{n no} - - {z z}^{' '})) = = \frac{22}{c c} ((\sqrt{{z z}^{' ' 22} + + {(({x x}^{' '} - - {u u}_{n no}))}^{22}} - - {z z}^{' '})),, n no = = 0,1 0,1,, . . . . . .,, L L - - 11 . .$

其中，c为超声在介质中的传播速度，则在(x′，z′)处的成像为Among them, c is the propagation speed of ultrasound in the medium, then the imaging at (x', z') is

$i i (({x x}^{' '},, {z z}^{' '})) = = {Σ Σ}_{n no = = 00}^{L L - - 11} \frac{{ω ω}_{n no} s the s (({u u}_{n no},, t t - - 22 {z z}^{' '} / / c c)) δ δ ((t t - - {τ τ}_{n no}))}{{r r}_{n no}}$

其中，ω_n为变迹函数，δ函数为单位抽样序列。Among them, ω _n is the apodization function, and the δ function is the unit sampling sequence.

频域SAFT被称为波数算法(wavenumber algorithm)，主要是基于相位延迟生成波束的方法，其源自于波动方程反演理论。基本原理是对格林函数(Green function)进行傅里叶分解。算法过程为先对超声数据进行二维傅里叶变换，得到二维频谱，然后利用非线性映射(Stoltmapping)和插值实现坐标变换，最后对变换后的频谱做滤波处理，并进行二维傅里叶反变换，生成空间-时间域下的重构图像。Frequency domain SAFT is called wavenumber algorithm (wavenumber algorithm), which is mainly based on the method of phase delay beam generation, which is derived from the wave equation inversion theory. The basic principle is to perform Fourier decomposition on the Green function. The algorithm process is to first perform two-dimensional Fourier transform on the ultrasonic data to obtain a two-dimensional spectrum, then use nonlinear mapping (Stoltmapping) and interpolation to realize coordinate transformation, and finally filter the transformed spectrum and perform two-dimensional Fourier Leaf inverse transform to generate reconstructed images in the space-time domain.

相位迁移(Phase Shift Migration)超声成像方法是将反射地震学(Reflection Seismology)中的迁移技术(Migration Technique)与频域SAFT相结合，而得到的一种超声成像方法。相位迁移技术最早由Gazdag在1978年提出，一直用于地震/雷达成像中，而最近提出的相位迁移超声成像方法才将该技术运用到了超声无损检测领域。Phase Shift Migration ultrasonic imaging method is an ultrasonic imaging method obtained by combining the Migration Technique in Reflection Seismology with frequency domain SAFT. The phase migration technology was first proposed by Gazdag in 1978 and has been used in seismic/radar imaging. The recently proposed phase migration ultrasonic imaging method has applied this technology to the field of ultrasonic nondestructive testing.

相位迁移超声成像方法将超声探测模型视为爆炸物反射模型，假设待测物体中的反射物在t＝0时刻同时爆炸，每一个反射物的爆炸强度正比于其反射率，整个场强用一组换能器来测量。其主要思想是根据从水平位置(即深度方向第一行)处观测到的声场外推以计算出深度方向其他位置处的声场。具体算法包含两个主要的步骤：第一步对时域数据进行二维傅里叶变换，得到二维频谱，The phase-shift ultrasonic imaging method regards the ultrasonic detection model as an explosive reflection model, assuming that the reflectors in the object to be tested explode at the same time at t=0, and the explosion intensity of each reflector is proportional to its reflectivity. group of transducers to measure. The main idea is to extrapolate the sound field observed from the horizontal position (ie, the first line in the depth direction) to calculate the sound field at other positions in the depth direction. The specific algorithm includes two main steps: the first step is to perform two-dimensional Fourier transform on the time domain data to obtain the two-dimensional spectrum,

S(k，ω)＝2D-FT(s(u_n，t))S(k, ω) = 2D-FT(s(u _n , t))

第二步是在深度方向的循环：先对上一次循环得到的二维频谱作相移，然后作二维傅里叶反变换并取t＝0，得到一行图像，The second step is the cycle in the depth direction: firstly perform phase shift on the two-dimensional frequency spectrum obtained in the previous cycle, then perform two-dimensional inverse Fourier transform and take t=0 to obtain a row of images,

s(u_n，t＝0)＝2D-IFT(S(k，ω)·α(Δz，k，ω))s(u _n , t=0)=2D-IFT(S(k,ω)·α(Δz,k,ω))

其中，

为深度方向上Δz步长所对应的相位迁移量，c为超声在介质中的传播速度，其取值恒定，因此，该超声成像方法使用的是恒定波速的相位迁移技术。in,

is the phase shift amount corresponding to the Δz step length in the depth direction, and c is the propagation velocity of ultrasound in the medium, and its value is constant. Therefore, the ultrasonic imaging method uses a phase shift technique with constant wave velocity.

从DAS的原理可知，时域SAFT算法需要根据换能器到反射物之间的相对距离以及波速来确定延时。而对于多层物体，声波在层与层的分界面处会发生折射、反射等声学特性，从而传播路径会发生改变。此外，物体各层的介质不一样，声波在各层的传播速度通常也不相等，因此，时延的计算不能再采用相对距离除以波速的简单计算方法，而是需要先获得波束的传播路径，然后逐段求出传播时间。其重点和难点在于快速而准确地找出波束的传播路径，而射线跟踪技术正是实现此目的之最佳方法。射线跟踪的原理主要基于Snell定理或Fermat原理，通过迭代计算找出费时最短的声波传播路径。射线跟踪技术的优点是该技术对多层物体的分界面形状并没有特殊的要求，能适用于任意形状物体的超声成像。但其严重的缺点是射线跟踪算法包含迭代运算，时间复杂度高，同时，DAS原理实际上是一种卷积运算，其计算量也比较大，成像时间比较长。如图2(b)所示，该方案生成该图像(5cm*5cm)就需要2小时。这一缺点严重限制了该技术的推广运用。It can be seen from the principle of DAS that the time-domain SAFT algorithm needs to determine the delay according to the relative distance between the transducer and the reflector and the wave speed. For multi-layer objects, the sound wave will have acoustic characteristics such as refraction and reflection at the interface between layers, so that the propagation path will change. In addition, the medium of each layer of the object is different, and the propagation speed of the sound wave in each layer is usually not equal. Therefore, the calculation of the delay can no longer use the simple calculation method of dividing the relative distance by the wave speed, but needs to obtain the propagation path of the beam first. , and then find the propagation time segment by segment. The important and difficult point is to find out the propagation path of the beam quickly and accurately, and the ray tracing technology is the best way to realize this purpose. The principle of ray tracing is mainly based on Snell's theorem or Fermat's principle, through iterative calculation to find out the shortest sound wave propagation path. The advantage of ray tracing technology is that the technology has no special requirements on the shape of the interface of multi-layer objects, and can be applied to ultrasonic imaging of objects with arbitrary shapes. But its serious disadvantage is that the ray-tracing algorithm contains iterative operations and has high time complexity. At the same time, the principle of DAS is actually a convolution operation, and its calculation load is relatively large, and the imaging time is relatively long. As shown in Figure 2(b), it takes 2 hours for this scheme to generate the image (5cm*5cm). This shortcoming seriously limits the popularization and application of this technology.

对于多层物体的成像，相位迁移超声成像方法无需求射线路径和卷积运算，成像速度可大幅提高，如生成图2(c)所示的图像仅需4分钟。并且相位迁移技术是一种基于频域的图像重构方法，比时域下的SAFT技术具有更高的距离向分辨率和更高的方位向分辨率以及更低的旁瓣，其成像质量更高。但是，由于相位迁移技术中假设声波的波速恒定，使得该方法只能适用于深度方向上存在异种介质而水平方向须为同种介质，即层与层介质之间的分界面为水平或者互相平行的直线/平面之类的规则情况，而对于深度方向和水平方向上均存在异种介质的物体，即包含非规则分界面(曲线/曲面)的物体，由于声波在各处介质中的波速不一样，而导致检测结果不够精确，如图2(c)所示，第二层分界面的成像结果已经出现严重的错误。For the imaging of multi-layer objects, the phase-shift ultrasonic imaging method does not require ray paths and convolution operations, and the imaging speed can be greatly improved. For example, it only takes 4 minutes to generate the image shown in Figure 2(c). And the phase migration technology is an image reconstruction method based on the frequency domain. Compared with the SAFT technology in the time domain, it has higher range resolution, higher azimuth resolution and lower side lobes, and its imaging quality is better. high. However, since the wave velocity of the sound wave is assumed to be constant in the phase migration technique, this method can only be applied to the presence of different media in the depth direction and the same media in the horizontal direction, that is, the interface between layers is horizontal or parallel to each other. For regular situations such as straight lines/planes, for objects with heterogeneous media in both the depth direction and the horizontal direction, that is, objects containing irregular interfaces (curves/surfaces), since the wave speeds of sound waves in various media are different , resulting in inaccurate detection results, as shown in Figure 2(c), serious errors have occurred in the imaging results of the second layer interface.

因此，在对含有非规则分界面的复杂分层物体的超声成像方面，目前还缺乏快速且准确有效的方法。Therefore, there is currently a lack of rapid, accurate and effective methods for ultrasound imaging of complex layered objects with irregular interfaces.

发明内容 Contents of the invention

本发明的目的在于提出一种超声成像方法，实现对含有非规则分界面的复杂分层物体的快速、精确、有效的成像。The object of the present invention is to propose an ultrasonic imaging method to realize rapid, accurate and effective imaging of complex layered objects with irregular interfaces.

本发明的特征在于，能够对规则分层物体或者仅含单一介质的物体成像，同时也能对含非水平且非互相平行分界面的非规则分层物体成像，且成像速度快。The present invention is characterized in that it can image regularly layered objects or objects containing only a single medium, and can also image irregularly layered objects containing non-horizontal and non-parallel interfaces at the same time, and the imaging speed is fast.

本发明的特征在于，采用相位迁移超声成像技术，并对相位迁移技术进行改进，在被测物体的深度方向和水平方向上分别考虑了存在异种介质的情况，在相位迁移公式中引入不同的波速，提出变波速的相位迁移技术，以使本超声成像方法适用于非规则分层物体。整体结构如图3所示，在成像过程中对分界线单独处理，在分界线所在的深度范围外，使用恒定波速的相位迁移公式和快速傅里叶变换技术，而在该范围内使用变波速的相位迁移技术和N次单输出快速傅里叶变换技术，以减少计算量。The present invention is characterized in that phase shift ultrasonic imaging technology is adopted, phase shift technology is improved, and the presence of different media is considered in the depth direction and horizontal direction of the measured object, and different wave velocities are introduced into the phase shift formula , a phase shift technique with variable wave velocity is proposed to make this ultrasonic imaging method suitable for irregularly layered objects. The overall structure is shown in Figure 3. In the imaging process, the boundary line is processed separately. Outside the depth range where the boundary line is located, the phase shift formula of constant wave velocity and fast Fourier transform technology are used, while variable wave velocity is used in this range. The advanced phase shift technology and N times single output fast Fourier transform technology to reduce the amount of calculation.

本发明的特征在于，依次含有以下步骤：The present invention is characterized in that it contains the following steps in sequence:

步骤(1)：构建一个由一台计算机、一个超声换能器、一套定位控制器和一个模数转换器组成的一个基于变波速相位迁移的用于对在横向和纵向均存在异种介质的三层异质物体在深度和水平两个方向形成的纵断面上作无损伤超声成像的系统，其中：Step (1): Construct a computer, an ultrasonic transducer, a set of positioning controller and an analog-to-digital converter based on variable wave velocity phase migration for detecting heterogeneous media in both horizontal and vertical directions. A system for non-invasive ultrasonic imaging of three-layer heterogeneous objects on longitudinal sections formed in depth and horizontal directions, in which:

所述超声换能器设有：与所述定位控制器的输出端相连的脉冲信号输入端，所述定位控制器的输入端与所述计算机相应的定位控制信号输出端相连，所述超声换能器还设有：与所述模数转换器的输入端相连的回波信号输出端，所述模数转换器的输出端与所述计算机的回波采样信号输入端相连。所述换能器由所述定位控制器控制，在三层异质物体表面以1步长/ms固定速率移动，所述定位控制器是控制所述超声换能器移动位置的传动装置，其参数由所述计算机输入，The ultrasonic transducer is provided with: a pulse signal input end connected to the output end of the positioning controller, the input end of the positioning controller is connected to the corresponding positioning control signal output end of the computer, and the ultrasonic transducer The energy device is also provided with an echo signal output end connected to the input end of the analog-to-digital converter, and the output end of the analog-to-digital converter is connected to the echo sampling signal input end of the computer. The transducer is controlled by the positioning controller and moves at a fixed rate of 1 step/ms on the surface of the three-layer heterogeneous object. The positioning controller is a transmission device that controls the moving position of the ultrasonic transducer. parameters are input by said computer,

三层异质物体沿x轴方向的水平长度为L，均分为L/Δx个区间，Δx为区间长度，也是所述超声换能器从x轴上的点0起向终点N_x-1止每次移动的步长，N_x＝1+L/Δx，所述超声换能器每次移动所达到的点称为探测点，共有N_x个，序号n_x＝0，1，...，N_x-1，所述定位控制器在每一个探测点处产生一个TTL脉冲，触发所述超声换能器向三层异质物体的垂直于x轴的深度方向z发射一个激励脉冲，随后换能器转为接收模式并开始计时，接收从被测物体反射的回波信号，所述模数转换器对所述换能器在探测点n_x处接收到的回波信号进行N_t次采样并存储到计算机中，采样序号n_t＝0，1，...，N_t-1，采样频率为f_s，其值为模数转换器预设，记s(z，x，t)为在横坐标为x＝n_x·Δx处接收到的回波信号在t＝n_t/f_s时刻的采样值；The horizontal length of the three-layer heterogeneous object along the x-axis direction is L, which is equally divided into L/Δx intervals, and Δx is the interval length, which is also the ultrasonic transducer from point 0 on the x-axis to the end point N _x -1 The step size of each movement, N _x =1+L/Δx, the point reached by the ultrasonic transducer each time it moves is called the detection point, there are N _x in total, and the serial number n _x =0, 1, .. ., N _x -1, the positioning controller generates a TTL pulse at each detection point, triggering the ultrasonic transducer to emit an excitation pulse to the depth direction z perpendicular to the x-axis of the three-layer heterogeneous object, Then the transducer switches to the receiving mode and starts timing to receive the echo signal reflected from the measured object, and the analog-to-digital converter performs N _t on the echo signal received by the transducer at the detection point n _x Sampling and storing in the computer, the sampling number n _t =0, 1, ..., N _t -1, the sampling frequency is f _s , and its value is the preset value of the analog-to-digital converter, denote s(z, x, t ) is the sampling value of the echo signal received at x=n _x Δx at the moment t=n _t /f _s at the abscissa;

步骤(2)：所述计算机从n_t＝0开始依序读取探测点n_x＝0处的采样值，并记录采样值第一次发生波动的采样序号

下标0表示探测点序号，然后，重复上述步骤依次读取n_x＝1，...，N_x-1各探测点处的采样值，并记录各点处采样值第一次发生波动的序号

得到N_x个探测点上总共N_x个采样值第一次发生波动的序号

从

中选取最小的序号

并计算该序号在所述三层异质物体的z方向上所对应的深度坐标值

和从

中选取最大的序号

其中v₁为第一层介质中的波速，分别通过这两个深度坐标值作两条平行于x轴的水平线，把第一层介质和第二层介质间的分界线包围起来，形成一个位于所述三层异质物体x-z纵断面上同时跨越第一、第二两层介质的第一个包围盒；Step (2): The computer reads the sampling value at the detection point n _x =0 sequentially from n _t =0, and records the sampling number at which the sampling value fluctuates for the first time

The subscript 0 represents the serial number of the detection point. Then, repeat the above steps to read the sampling values at each detection point of n _x = 1, ..., N _x -1 in sequence, and record the first fluctuation of the sampling value at each point. serial number

Get the serial number of the first fluctuation of a total of N _x sampling values on N _x detection points

from

Choose the smallest serial number

And calculate the depth coordinate value corresponding to the serial number in the z direction of the three-layer heterogeneous object

and from

Select the largest serial number from

Among them, v ₁ is the wave velocity in the first layer of medium, draw two horizontal lines parallel to the x-axis through these two depth coordinates, and enclose the boundary line between the first layer of medium and the second layer of medium, forming a The first bounding box of the three-layer heterogeneous object simultaneously spanning the first and second layers of media on the xz longitudinal section;

步骤(3)：沿着时间轴t逐个对各探测点处的采样值s(z，x，t)进行一维快速傅里叶变换1D-FFT_t(s(z，x，t))，取z的初始值z₀＝0，得到x-ω域二维信号：Step (3): Perform one-dimensional fast Fourier transform 1D-FFT _t (s(z, x, t)) to the sampling values s(z, x, t) at each detection point along the time axis t one by one, Taking the initial value z ₀ =0 of z, the two-dimensional signal in the x-ω domain is obtained:

S(z₀，x，ω)＝1D-FFT_t(s(z₀，x，t))S(z ₀ , x, ω)=1D-FFT _t (s(z ₀ , x, t))

其中，x为探测点序号为n_x处的水平向坐标值，ω为采样角频率，然后，对上述x-ω域二维信号S(z₀，x，ω)沿水平轴x依次进行一维快速傅里叶变换1D-FFT_x(S(z₀，x，ω))，得到z₀处的k-ω域二维频谱S(z₀，k，ω)：Among them, x is the horizontal coordinate value at the detection point number n _x , ω is the sampling angular frequency, and then, the above-mentioned two-dimensional signal S(z ₀ , x, ω) in the x-ω domain is sequentially processed along the horizontal axis x One-dimensional fast Fourier transform 1D-FFT _x (S(z ₀ , x, ω)), to obtain the k-ω domain two-dimensional spectrum S( _{z 0} _, k, ω) at z 0:

S(z₀，k，ω)＝1D-FFT_x(1D-FFT_t(s(z₀，x，t)))S(z ₀ , k, ω) = 1D-FFT _x (1D-FFT _t (s(z ₀ , x, t)))

其中，k为波数矢的x向分量，序号为n_k，所述波数矢表示波动相位的空间梯度，Among them, k is the x-direction component of the wave number vector, the sequence number is n _k , and the wave number vector represents the spatial gradient of the wave phase,

0≤n_k≤N_x-1，且n_k为整数；

0≤n _k ≤N _x -1, and n _k is an integer;

0≤n_ω≤N_t-1，且n_ω为整数；

0≤n _ω ≤N _t -1, and n _ω is an integer;

其中，n_ω为采样角频率ω的序号，取值在0～(N_t-1)之间，取t＝t₀＝0，一维采样数据s(z₀，x，t₀)即为三层异质物体纵断面上第z＝z₀＝0行各像素点的像素值；Among them, n _ω is the serial number of the sampling angular frequency ω, the value is between 0 and (N _t -1), taking t=t ₀ =0, the one-dimensional sampling data s(z ₀ , x, t ₀ ) is The pixel value of each pixel point in the z=z ₀ =0 row on the vertical section of the three-layer heterogeneous object;

步骤(4)：依次按以下步骤生成z方向同质区z₀＜z≤z_min内的分别由各行像素点形成的各行图像：Step (4): Follow the steps below to generate each row of images in the homogeneous area z ₀ <z≤z _min in the z direction, respectively formed by each row of pixels:

步骤(4.1)：利用下述相位迁移公式计算所述纵断面在z_min处的二维频谱S(z_min，k，ω)：Step (4.1): Use the following phase shift formula to calculate the two-dimensional spectrum S(z _min , k, ω) of the longitudinal section at z _min :

$S S (({z z}_{min min},, k k,, ω ω)) = = S S (({z z}_{00},, k k,, ω ω)) {Π Π}_{z z = = 11}^{{z z}_{min min}} α α ((Δz Δz,, k k,, ω ω,, {v v}_{z z}))$

其中，in,

其中，i为虚数单位，z₀＝0，ω为采样角频率，v_z为被测物体在深度方向上第z行介质中的声波波速，在z≤z_min时，v_z＝v₁，v₁为第一层介质中声波的波速，Δz为深度方向z上相邻两个像素之间的间隔距离，α(Δz，k，ω，v_z)为恒定波速下的相位迁移量，Among them, i is the imaginary number unit, z ₀ =0, ω is the sampling angular frequency, v _z is the sound wave velocity in the medium of the zth line of the measured object in the depth direction, when z≤z _min , v _z =v ₁ , v ₁ is the wave velocity of the sound wave in the first layer of medium, Δz is the distance between two adjacent pixels in the depth direction z, α(Δz, k, ω, v _z ) is the phase shift at a constant wave velocity,

步骤(4.2)：依次执行以下步骤：Step (4.2): Perform the following steps in sequence:

步骤(4.2.1)：取t＝t₀＝0，对所述步骤(4.1)得到的S(z_min，k，ω)作自变量ω的高斯积分，Step (4.2.1): get t=t ₀ =0, do Gaussian integral of independent variable ω to S (z _min , k, ω) that described step (4.1) obtains,

步骤(4.2.2)：对所述步骤(4.2.1)的高斯积分结果取共轭，再对k作一维快速傅里叶变换，对计算结果再取共轭后乘以1/N_x，得到第一层介质内深度方向上坐标值为z_min的一行像素点所形成的图像：Step (4.2.2): take the conjugate of the Gaussian integral result of the step (4.2.1), then perform a one-dimensional fast Fourier transform on k, take the conjugate of the calculation result and multiply it by 1/N _x , to obtain the image formed by a row of pixels with the coordinate value z _min in the depth direction in the first layer medium:

$s the s (({z z}_{min min},, x x,, {t t}_{00})) = = \frac{11}{{N N}_{x x}} {{11 D D. - - {FFT FFT}_{k k} {{[[\frac{11}{22 π π} {&Integral; &Integral;}_{- - π π}^{π π} S S (({z z}_{min min},, k k,, ω ω)) dω dω {]]}^{* *})) {}}}^{* *}$

在z方向第一层同质区z₀＜z≤z_min内，各行图像均相同：In the homogeneous area of the first layer z ₀ <z≤z _min in the z direction, the images of each row are the same:

s(z，x，t₀)＝s(z_min，x，t₀)，z₀＜z≤z_min；s(z, x, t ₀ )=s(z _min , x, t ₀ ), z ₀ <z≤z _min ;

步骤(5)：取步骤(4)中的v_z＝v₁，对z值以Δz为步长循环执行下述步骤(5.1)、步骤(5.2)，直到z≥z_max止，生成第一包围盒内深度方向z_min＜z≤z_max区间分别由各行像素点所形成的各行图像，其步骤如下：Step (5): Take v _z =v ₁ in step (4), execute the following steps (5.1) and (5.2) cyclically for the value of z with Δz as the step size, until z≥z _max , and generate the first Each row of images formed by each row of pixels in the depth direction z _min < z ≤ z _max in the bounding box, the steps are as follows:

步骤(5.1)：利用所述恒定波速下的相位迁移公式计算z+Δz处的二维频谱，z_min＜z+Δz≤z_max：Step (5.1): Calculate the two-dimensional frequency spectrum at z+Δz using the phase shift formula at a constant wave velocity, z _min <z+Δz≤z _max :

S(z+Δz，k，ω)＝S(z，k，ω)·α(Δz，k，ω，v_z)，S(z+Δz, k, ω) = S(z, k, ω) α(Δz, k, ω, v _z ),

步骤(5.2)：依次执行以下步骤：Step (5.2): Perform the following steps in sequence:

步骤(5.2.1)：取t＝t₀＝0，对所述步骤(5.1)得到的S(z+Δz，k，ω)作自变量ω的高斯积分，Step (5.2.1): get t=t ₀ =0, do Gaussian integral of independent variable ω to S (z+Δz, k, ω) that described step (5.1) obtains,

步骤(5.2.2)：对所述步骤(5.2.1)的高斯积分结果取共轭，再对k作一维快速傅里叶变换，对计算结果再取共轭后乘以1/N_x，得到深度方向上坐标值为z+Δz的一行像素点所形成的图像，z_min＜z+Δz≤z_max：Step (5.2.2): take the conjugate of the Gaussian integral result of the step (5.2.1), then perform a one-dimensional fast Fourier transform on k, take the conjugate of the calculation result and multiply it by 1/N _x , to obtain an image formed by a row of pixels whose coordinate value is z+Δz in the depth direction, z _min <z+Δz≤z _max :

$s the s ((z z + + Δz Δz,, x x,, {t t}_{00})) = = \frac{11}{{N N}_{x x}} {{11 D D. - - {FFT FFT}_{k k} {{[[\frac{11}{22 π π} {&Integral; &Integral;}_{- - π π}^{π π} S S ((z z + + Δz Δz,, k k,, ω ω)) dω dω {]]}^{* *})) {}}}^{* *};;$

步骤(6)：按以下步骤修正所述第一包围盒内深度方向z_min至z_max区间z_min＜z≤z_max中的各行图像，以消除第二层和第一层介质间的不同质而导致的误差：Step (6): according to the following steps, correct the images of each line in the interval z _min < z ≤ z _max in the depth direction z _min to z _max in the first bounding box, so as to eliminate the heterogeneity between the second layer and the first layer medium The resulting error:

步骤(6.1)：以步骤(5)得到的z_min至z_max区间图像块作为输入量，使用Canny算子边缘提取算法提取所述z_min、z_max之间的第一层介质与第二层介质的分界线c₁(x，z)，Step (6.1): Using the z _min to z _max interval image block obtained in step (5) as an input, use the Canny operator edge extraction algorithm to extract the first layer of medium and the second layer between z _min and z _max The dividing line c ₁ (x, z) of the medium,

步骤(6.2)：从z＝z_min开始，直到z≥z_max止，以所述步骤(4.1)得到的z_min处的二维频谱S(z_min，k，ω)作起始值，以Δz为步长循环执行步骤(6.2.1)和步骤(6.2.2)：Step (6.2): starting from z=z _min until z≥z _max , the two-dimensional spectrum S (z _min , k, ω) at the z _min obtained in the step (4.1) is used as the initial value, and Δz is the step size and executes step (6.2.1) and step (6.2.2) cyclically:

步骤(6.2.1)：使用N次单输出快速傅里叶反变换方法，按下式计算得到所述第一个包围盒内z+Δz处的x-ω域二维信号：Step (6.2.1): Use the N-times single-output inverse fast Fourier transform method to calculate the x-ω domain two-dimensional signal at z+Δz in the first bounding box as follows:

$S S ((z z + + Δz Δz,, x x = = {n no}_{x x} \cdot &Center Dot; Δx Δx,, ω ω)) = = \frac{11}{{N N}_{x x}} {Σ Σ}_{{n no}_{k k} = = 00}^{{N N}_{x x} - - 11} [[S S ((z z,, k k = = \frac{{22 πn πn}_{k k}}{{N N}_{x x} Δx Δx},, ω ω)) α α ((Δz Δz,, k k = = \frac{{22 πn πn}_{k k}}{{N N}_{x x} Δx Δx},, ω ω,, {v v}_{z z} ((x x = = {n no}_{x x} \cdot &Center Dot; Δx Δx))))]] {e e}^{i i \frac{22 π π}{{N N}_{x x}} {n no}_{k k} {n no}_{x x}}$

其中，N的取值等于N_x，n_x＝0，1，...，N_x-1，n_k＝0，1，...，N_x-1，α(Δz，k，ω，v_z(x))为变波速的相位迁移量：Wherein, the value of N is equal to N _x , n _x =0, 1, ..., N _x -1, n _k =0, 1, ..., N _x -1, α(Δz, k, ω, v _z (x)) is the phase shift of variable wave velocity:

其中，k、ω的取值同步骤(3)，v_z(x)是坐标(x，z)处的介质中声波的波速：Among them, the values of k and ω are the same as step (3), and v _z (x) is the wave velocity of the sound wave in the medium at the coordinates (x, z):

v₁为位于分界线c₁(x，z)之上的第一层介质中声波的传播速度，v₂为位于分界线c₁(x，z)之下的第二层介质中声波的传播速度，v ₁ is the propagation velocity of the sound wave in the first layer of medium above the dividing line c ₁ (x, z), and v ₂ is the propagation velocity of the sound wave in the second layer of medium below the dividing line c ₁ (x, z) speed,

计算步骤如下：The calculation steps are as follows:

步骤(a)：分别对不同的n_x，n_k(n_x＝0，1，...，N_x-1；n_k＝0，1，...，N_x-1)计算Step (a): calculate for different n _x , n _k (n _x =0, 1, ..., N _x -1; n _k =0, 1, ..., N _x -1) respectively

$G G (({n no}_{x x},, {n no}_{k k})) = = S S ((z z,, k k = = \frac{{22 πn πn}_{k k}}{{N N}_{x x} Δx Δx},, ω ω)) α α ((Δz Δz,, k k = = \frac{{22 πn πn}_{k k}}{{N N}_{x x} Δx Δx},, ω ω {,, v v}_{z z} ((x x = = {n no}_{x x} \cdot &Center Dot; Δx Δx)))),,$

步骤(b)：对不同的n_x(n_x＝0，1，...，N_x-1)，以G(n_x，n_k)，n_k＝0，1，...，N_x-1序列作为裁剪输出FFT算法的输入量，对n_k作G(n_x，n_k)离散傅里叶反变换，但仅取n_x这一时域点的值，作为深度方向上坐标值为z+Δz处的x-ω域二维信号：Step (b): For different n _x (n _x =0, 1, ..., N _x -1), with G(n _x , n _k ), n _k =0, 1, ..., N _{The x} -1 sequence is used as the input of the clipping output FFT algorithm, and the G(n _x , n _k ) inverse discrete Fourier transform is performed on n _k , but only the value of the time domain point n _x is taken as the coordinate value in the depth direction is the two-dimensional signal in the x-ω domain at z+Δz:

$S S ((z z + + Δz Δz,, x x = = {n no}_{x x} \cdot &Center Dot; Δx Δx,, ω ω)) = = \frac{11}{{N N}_{x x}} {Σ Σ}_{{n no}_{k k} = = 00}^{{N N}_{x x} - - 11} G G (({n no}_{x x},, {n no}_{k k})) {e e}^{i i \frac{22 π π}{{N N}_{x x}} {n no}_{k k} {n no}_{x x}},,$

步骤(6.2.2)：取t＝t₀＝0，对所述步骤(6.2.1)得到的S(z+Δz，x，ω)作自变量ω的高斯积分，得到深度方向上坐标值为z+Δz的一行像素点所形成的图像，z_min＜z+Δz＜z_max：Step (6.2.2): Take t=t ₀ =0, do the Gaussian integral of the independent variable ω for the S(z+Δz, x, ω) obtained in the step (6.2.1), and obtain the coordinate value in the depth direction It is an image formed by a row of pixels of z+Δz, z _min <z+Δz<z _max :

$s the s ((z z + + Δz Δz,, x x,, {t t}_{00})) = = \frac{11}{22 π π} {&Integral; &Integral;}_{- - π π}^{π π} S S ((z z + + Δz Δz,, x x,, ω ω)) dω dω$

同时，保存S(z+Δz，x，ω)并作一维傅里叶变换得到S(z+Δz，k，ω)，作为下一轮迭代的初始值；At the same time, save S(z+Δz, x, ω) and perform one-dimensional Fourier transform to obtain S(z+Δz, k, ω) as the initial value of the next iteration;

步骤(7)：按步骤(5)所述的方法，取v_z＝v₂，对z值以Δz为步长循环执行所述步骤(5.1)、步骤(5.2)，直到z≥z_depth止，生成深度方向z_max＜z≤z_depth区间各行像素点所形成的各行图像，v₂为第二层介质中声波的传播速度，z_depth表示三层异质物体在深度方向z上的最大坐标值，为系统预设值；Step (7): According to the method described in step (5), take v _z =v ₂ , and execute the steps (5.1) and (5.2) cyclically for the value of z with Δz as the step size until z≥z _depth , to generate images of each row formed by each row of pixels in the depth direction z _max < z ≤ z _depth interval, v ₂ is the propagation speed of the sound wave in the second layer medium, z _depth represents the maximum coordinate of the three-layer heterogeneous object in the depth direction z The value is the system default value;

步骤(8)：以步骤(7)得到的z_max至z_depth区间图像块作为输入量，使用Canny算子边缘提取算法提取所述z_max、z_depth之间的第二层介质与第三层介质的分界线c₂(x，z)，求分界线c₂(x，z)在深度方向z上的最小坐标值z′_min和最大坐标值z′_max，分别通过所述z′_min和z′_max两个深度坐标值作两条平行于x轴的水平线，把第二层介质和第三层介质间的分界线包围起来，形成一个位于所述三层异质物体x-z纵断面上同时跨越第二、第三两层介质的第二个包围盒；Step (8): Using the z _max to z _depth interval image block obtained in step (7) as an input, use the Canny operator edge extraction algorithm to extract the second layer medium and the third layer between the z _max and z _depth For the boundary line c ₂ (x, z) of the medium, find the minimum coordinate value z′ _min and the maximum coordinate value z′ _max of the boundary line c ₂ (x, z) in the depth direction z, and pass through the z′ _min and The two depth coordinates of z′ _max draw two horizontal lines parallel to the x-axis, enclosing the boundary line between the second layer of medium and the third layer of medium, forming a simultaneous xz longitudinal section of the three-layer heterogeneous object The second bounding box spanning the second and third layers of media;

步骤(9)：按照步骤(6.2)所述的方法修正第二包围盒内深度方向z′_min至z′_max区间z′_min＜z≤z′_max中的各行图像，以消除第三层和第二层介质间的不同质而导致的误差；Step (9): According to the method described in step (6.2), correct each row of images in the second bounding box in the depth direction _z′min to _z′max interval _z′min < _z≤z′max , to eliminate the third layer and The error caused by the inhomogeneity of the second layer of media;

步骤(10)：按步骤(5)所述的方法，取v_z＝v₃，对z值以z′_max为起始值，以Δz为步长循环执行所述步骤(5.1)、步骤(5.2)，直到z≥z_depth止，生成第三层介质内除了所述第二包围盒以外部分各行像素点所形成的各行图像，v₃为第三层介质中声波的传播速度。Step (10): According to the method described in step (5), take v _z =v ₃ , take z′ _max as the initial value for the z value, and execute the steps (5.1) and steps ( 5.2), until z≥z _depth , generate images of each row of pixels in the third layer of the medium except for the second bounding box, v ₃ is the propagation velocity of the sound wave in the third layer of medium.

本发明的优点在于，成像精度高，对于图2(a)所示的被测物体，若换能器直径为0.5mm，换能器移动步长为0.5mm，换能器发射的超声波中心频率为5MHz，采样频率100MHz，深度方向的步长Δz取为0.05mm，则第一层分界面误差可以控制在0.4mm以内，第二层分界面可以控制在1.0mm以内(如图4)。The present invention has the advantages of high imaging accuracy. For the measured object shown in Fig. 2(a), if the diameter of the transducer is 0.5mm, the moving step of the transducer is 0.5mm, and the center frequency of the ultrasonic waves emitted by the transducer If the sampling frequency is 5MHz, the sampling frequency is 100MHz, and the step size Δz in the depth direction is set to 0.05mm, the error of the interface of the first layer can be controlled within 0.4mm, and the interface of the second layer can be controlled within 1.0mm (as shown in Figure 4).

本发明与现有的技术相比有以下的优点：Compared with the prior art, the present invention has the following advantages:

1.能够对含有非规则分界面的物体进行快速超声成像；1. Capable of rapid ultrasonic imaging of objects with irregular interfaces;

2.成像的分辨率和精度较高。2. The imaging resolution and precision are high.

附图说明 Description of drawings

图1是SAFT的工作模型：X向为换能器扫描方向，Z向为深度方向，u_n为换能器以等步长移动时的位置，也即探测点的位置，(x’，z’)为目标反射物的坐标，r_n为目标反射物至换能器的距离。Figure 1 is the working model of SAFT: the X direction is the scanning direction of the transducer, the Z direction is the depth direction, u _n is the position of the transducer when it moves with equal steps, that is, the position of the detection point, (x', z ') is the coordinates of the target reflector, r _n is the distance from the target reflector to the transducer.

图2是非规则分层物体超声成像方法在相同实验环境下的对比图：2(a)为被测物体的剖面图；2(b)为时域SAFT与射线跟踪技术结合对2(a)生成的图像，成像时间2小时以上，而且对算法的参数非常敏感，鲁棒性不是很好；2(c)为恒定波速的相位迁移超声成像技术生成的图像，成像时间4分钟，其误差较大；2(d)为本成像方法生成的图像，成像时间6分钟。Figure 2 is a comparison of ultrasonic imaging methods for irregularly layered objects in the same experimental environment: 2(a) is the profile of the measured object; 2(b) is the combination of time-domain SAFT and ray tracing technology generated by 2(a) 2(c) is an image generated by the phase-shift ultrasonic imaging technology with constant wave velocity, and the imaging time is 4 minutes, and its error is relatively large ; 2(d) is the image generated by this imaging method, and the imaging time is 6 minutes.

图3是本成像方法总体过程示意图。FIG. 3 is a schematic diagram of the overall process of the imaging method.

图4是对图2(d)的实验数据分析图：4(a)为第一层分界线误差图，实曲线为被测物体第一层分界线的标准值，虚点曲线为本成像方法获得的第一层分界线，短划线曲线为本成像方法所获得的第一层分界线的误差曲线；4(b)为第二层分界线误差图，实曲线为被测物体第二层分界线的标准值，虚点曲线为本成像方法获得的第二层分界线，短划线曲线为本成像方法所获得的第二层分界线的误差曲线。Fig. 4 is an analysis diagram of the experimental data of Fig. 2(d): 4(a) is the error map of the first layer boundary line, the solid curve is the standard value of the first layer boundary line of the measured object, and the dotted point curve is the imaging method The obtained first layer boundary line, the dashed line curve is the error curve of the first layer boundary line obtained by this imaging method; 4(b) is the second layer boundary line error map, and the solid curve is the second layer of the measured object The standard value of the boundary line, the dashed point curve is the second layer boundary line obtained by this imaging method, and the dashed line curve is the error curve of the second layer boundary line obtained by this imaging method.

图5是本超声成像系统流程示意图。FIG. 5 is a schematic flow chart of the ultrasound imaging system.

图6是本超声成像硬件系统结构图。FIG. 6 is a structural diagram of the ultrasonic imaging hardware system.

图7是超声换能器工作示意图。Fig. 7 is a schematic diagram of the operation of the ultrasonic transducer.

图8是基于变波速相位迁移技术的多层物体无损检测超声成像算法流程图。Fig. 8 is a flowchart of an ultrasonic imaging algorithm for non-destructive testing of multi-layer objects based on variable wave velocity phase migration technology.

图9是本超声成像方法中求包围盒示意图。Fig. 9 is a schematic diagram of obtaining a bounding box in the ultrasonic imaging method.

图10是N＝8时的FFT蝶形运算示意图。FIG. 10 is a schematic diagram of FFT butterfly operation when N=8.

具体实施方式 Detailed ways

本发明的具体实施过程包含三部分(如图5)：超声数据获取、成像计算和图像显示。硬件平台系统结构图如图6所示，超声成像系统由一台计算机、一个超声换能器、一套定位控制器和一个模数转换器组成，超声换能器的脉冲信号输入端与定位控制器的输出端相连，定位控制器的输入端与计算机的定位控制信号输出端相连。超声换能器的回波信号输出端与模数转换器的输入端相连，模数转换器的输出端与计算机的回波采样信号输入端相连。The specific implementation process of the present invention includes three parts (as shown in FIG. 5 ): ultrasonic data acquisition, imaging calculation and image display. The hardware platform system structure diagram is shown in Figure 6. The ultrasonic imaging system consists of a computer, an ultrasonic transducer, a positioning controller and an analog-to-digital converter. The pulse signal input terminal of the ultrasonic transducer is connected to the positioning control The output end of the controller is connected, and the input end of the positioning controller is connected with the output end of the positioning control signal of the computer. The echo signal output end of the ultrasonic transducer is connected with the input end of the analog-digital converter, and the output end of the analog-digital converter is connected with the echo sampling signal input end of the computer.

超声数据获取的方式可选择：接触式，即超声换能器直接接触待测物体的表面；或者液浸式，即待测物体置于液体中，超声换能器在液体表面进行测量。根据检测方式的不同，换能器可以使用接触式探头或者液浸式探头。系统使用单个发射/接收换能器，换能器通过定位控制器在被测物体表面或者液体中以均匀的步长Δx沿x方向(如图7)以约1步/ms固定速率移动。控制器在每一个目标位置稳定的瞬间产生一个TTL(晶体管-晶体管逻辑电平)脉冲，该脉冲用来触发超声换能器向被测物体与x方向相垂直的深度方向发射一个激励脉冲，随后换能器转为接收模式并开始计时，接收从被测物体反射的回波。换能器发射脉冲及接收回波的各个位置处为探测点。换能器接收到的回波信号由模数转换器采集并存储在存储器中。Δx需根据待测物体的实际大小和成像精度要求来确定，其值越小，生成的图像越精确，但计算时间也越长。The method of ultrasonic data acquisition can be selected: contact type, that is, the ultrasonic transducer directly contacts the surface of the object to be measured; or liquid immersion type, that is, the object to be measured is placed in the liquid, and the ultrasonic transducer measures on the liquid surface. Depending on the detection method, the transducer can use a contact probe or a liquid immersion probe. The system uses a single transmitting/receiving transducer, and the transducer moves with a uniform step size Δx along the x direction (as shown in Figure 7) at a fixed rate of about 1 step/ms through the positioning controller on the surface of the measured object or in the liquid. The controller generates a TTL (transistor-transistor logic level) pulse at the moment when the target position is stable, which is used to trigger the ultrasonic transducer to emit an excitation pulse to the depth direction of the measured object perpendicular to the x direction, and then The transducer switches to receive mode and starts timing, receiving echoes reflected from the object under test. The locations where the transducer emits pulses and receives echoes are detection points. The echo signal received by the transducer is collected by an analog-to-digital converter and stored in a memory. Δx needs to be determined according to the actual size of the object to be measured and the imaging accuracy requirements. The smaller the value, the more accurate the generated image, but the longer the calculation time.

成像计算就是以被测物体在一个纵断面上各探测点处的采样数据作为计算机输入，然后按以下成像步骤(流程图参见图8)计算被测物体的纵断面图像，其中s(z，x，t)表示在水平方向x处的换能器在t时刻接收到的信号的采样值；z表示深度方向的坐标值；Δz表示深度方向的步长(与精度有关)，即所生成的图像在纵向(亦即深度方向)上相邻两个像素点的间距；z_depth为预设的图像在纵向上的最大深度值：The imaging calculation is to take the sampling data of the measured object at each detection point on a longitudinal section as computer input, and then calculate the longitudinal section image of the measured object according to the following imaging steps (see Figure 8 for the flow chart), where s(z, x , t) represents the sampling value of the signal received by the transducer at the horizontal direction x at time t; z represents the coordinate value in the depth direction; Δz represents the step size in the depth direction (related to precision), that is, the generated image The distance between two adjacent pixels in the vertical direction (that is, the depth direction); z _depth is the maximum depth value of the preset image in the vertical direction:

步骤(1)：第一层介质为同构介质，其超声反射信号相同，而在第一层与第二层分界处，超声反射信号会产生波动。则根据超声回波信号的采样数据，搜索各探测点处第一次发生信号波动的时刻，并从这些波动时刻中找出最小时刻值t_min和最大时刻值t_max，然后根据第一层介质中的波速v₁求第一层与第二层介质间分界线的包围盒，获得包围盒在深度方向(即z向)上的两条边界线：z_min＝v₁·t_min和z_max＝v₁·t_max(如图9)；Step (1): The first layer of medium is isomorphic, and its ultrasonic reflection signal is the same, but at the boundary between the first layer and the second layer, the ultrasonic reflection signal will fluctuate. Then, according to the sampling data of the ultrasonic echo signal, search for the moment when the signal fluctuation occurs for the first time at each detection point, and find out the minimum moment value t _min and the maximum moment value t _max from these fluctuation moments, and then according to the first layer of medium Calculate the bounding box of the boundary line between the first layer and the second layer of the wave velocity v ₁ in , and obtain two boundary lines of the bounding box in the depth direction (ie z direction): z _min =v ₁ ·t _min and z _max =v ₁ ·t _max (as shown in Figure 9);

步骤(2)：取深度方向初始值z＝z₀＝0，对参数x和t进行二维傅里叶变换，将原始的超声回波信号变换为k-ω域信号：Step (2): Take the initial value z=z ₀ =0 in the depth direction, perform two-dimensional Fourier transform on the parameters x and t, and transform the original ultrasonic echo signal into a k-ω domain signal:

S(z₀，k，ω)＝∫∫s(z₀，x，t)e^-ikxe^-iωtdxdt (1)S(z ₀ , k, ω)=∫∫s(z ₀ , x, t)e ^-ikx e ^-iωt dxdt (1)

其中，ω为采样角频率，k为波数矢(wavenumber vector)的x向分量(波数矢指波动相位的空间梯度)；Wherein, ω is sampling angular frequency, and k is the x-direction component of wavenumber vector (wavenumber vector) (wavenumber vector refers to the spatial gradient of wave phase);

步骤(3)：生成深度方向z₀至z_min区间(z₀≤z≤z_min)的图像。由于被测物体在此段区间内的介质完全相同，因此此段图像也完全一样：Step (3): Generating images in the interval z ₀ to z _min (z ₀ ≤ z ≤ _{z min} ) in the depth direction. Since the medium of the measured object in this interval is exactly the same, the images in this section are also exactly the same:

步骤(3.1)：利用相位迁移公式计算z_min处的k-ω域信号：Step (3.1): Calculate the k-ω domain signal at z _min using the phase shift formula:

$S S (({z z}_{min min},, k k,, ω ω)) = = S S (({z z}_{00},, k k,, ω ω)) {Π Π}_{z z = = 11}^{{z z}_{min min}} α α ((Δz Δz,, k k,, ω ω,, {v v}_{z z})) - - - - - - ((22))$

其中，in,

i为虚数单位，v_z为被测物体在深度方向上第z行介质中的声波波速，当z≤z_min时，v_z＝v₁，v₁为第一层介质中声波的波速；i is the imaginary number unit, v _z is the sound wave velocity in the medium of the zth line of the measured object in the depth direction, when z≤z _min , v _z =v ₁ , v ₁ is the wave velocity of the sound wave in the first layer of medium;

步骤(3.2)：作二维傅里叶反变换，令t＝t₀＝0，得到纵向坐标值为z的一行图像：Step (3.2): do two-dimensional inverse Fourier transform, make t=t ₀ =0, obtain a row of images whose longitudinal coordinate value is z:

$s the s ((z z,, x x,, {t t}_{00})) = = s the s (({z z}_{min min},, x x,, {t t}_{00})) = = \frac{11}{{44 π π}^{22}} &Integral; &Integral; &Integral; &Integral; S S (({z z}_{min min},, k k,, ω ω)) {e e}^{ikx ikx} {e e}^{iωt iωt} dkdω dkdω,, (({z z}_{00} \leq \leq z z \leq \leq {z z}_{min min})) - - - - - - ((33))$

步骤(4)：生成深度方向z_min至z_max区间(z_min＜z≤z_max)的图像。取v_z＝v₁，对z值以Δz步长循环执行步骤(4.1)和步骤(4.2)直到z≥z_max：Step (4): Generating an image in the interval from z _min to z _max (z _min <z≤z _max ) in the depth direction. Take v _z =v ₁ , execute step (4.1) and step (4.2) cyclically with Δz step size for z value until z≥z _max :

步骤(4.1)：利用相位迁移公式计算z+Δz处的k-ω域信号：Step (4.1): Calculate the k-ω domain signal at z+Δz using the phase shift formula:

S(z+Δz，k，ω)＝S(z，k，ω)α(Δz，k，ω，v_z) (4)S(z+Δz,k,ω)=S(z,k,ω)α(Δz,k,ω,v _z ) (4)

步骤(4.2)：令z←z+Δz，t＝t₀＝0，作二维傅里叶反变换，得到纵向坐标值为z的一行图像：Step (4.2): Make z←z+Δz, t=t ₀ =0, do two-dimensional inverse Fourier transform, obtain a row of images whose longitudinal coordinate value is z:

$s the s ((z z,, x x,, {t t}_{00})) = = \frac{11}{{44 π π}^{22}} &Integral; &Integral; &Integral; &Integral; S S ((z z,, k k,, ω ω)) {e e}^{ikx ikx} {e e}^{iωt iωt} dkdω dkdω - - - - - - ((55))$

步骤(5)：修正z_min至z_max区间(z_min＜z≤z_max)的图像。由于步骤(4)中的计算仅使用了第一层的波速，即将该区间仍视为同种介质，得到的计算结果必然存在误差，所以此步需要修正。具体包含两步：Step (5): Correcting the image in the interval from z _min to z _max (z _min <z≤z _max ). Since the calculation in step (4) only uses the wave velocity of the first layer, that is, the interval is still regarded as the same medium, there must be errors in the calculation results obtained, so this step needs to be corrected. Specifically, it includes two steps:

步骤(5.1)：对步骤(4)得到的z_min至z_max区间的图像，使用Canny边缘提取算子提取z_min与z_max之间的分界线c(x，z)，z_min≤z≤z_max。Canny边缘提取算子的计算过程分为四步：Step (5.1): For the image in the interval from z _min to z _max obtained in step (4), use the Canny edge extraction operator to extract the dividing line c(x, z) between z _min and z _max , z _min ≤ z ≤ _zmax . The calculation process of the Canny edge extraction operator is divided into four steps:

步骤(a)：图像平滑。对原图像用二维高斯函数平滑滤波得到图像I(x，z)；Step (a): Image smoothing. The original image is smoothed and filtered with a two-dimensional Gaussian function to obtain an image I(x, z);

步骤(b)：计算图像I(x，z)每个像素的梯度M和方向Q。采用2×2模板作为对x方向和z方向的偏微分的一阶近似，即Step (b): Calculate the gradient M and direction Q of each pixel of the image I(x, z). A 2×2 template is used as a first-order approximation to the partial differential in the x-direction and z-direction, namely

$p p = = \frac{11}{22} (\begin{matrix} - - 11 & 11 \\ - - 11 & 11 \end{matrix}),,$ $q q = = \frac{11}{22} (\begin{matrix} 11 & 11 \\ - - 11 & - - 11 \end{matrix})$

则梯度大小M和方向Q为Then the gradient size M and direction Q are

$M m = = \sqrt{p p \times \times p p + + q q \times \times q q},,$ $Q Q = = arctan arctan ((\frac{q q}{p p}));;$

步骤(c)：梯度图像的非极大值抑制。对每个像素点I(x，z)，若其梯度值M(x，z)小于沿着其梯度方向Q(x，z)的两个相邻点的梯度幅值，说明该点不是边缘点，则将I(x，z)的灰度设为0；Step (c): Non-maximum suppression of the gradient image. For each pixel point I(x, z), if its gradient value M(x, z) is smaller than the gradient magnitude of two adjacent points along its gradient direction Q(x, z), it means that the point is not an edge point, set the gray level of I(x, z) to 0;

步骤(d)：双阈值化处理。设定好双阈值方法检测和连接边缘需要的低阈值Low和高阈值High，对梯度图像进行双阈值化处理。梯度幅值大于高阈值High的是边缘；梯度幅值小于低阈值Low的不是边缘；梯度幅值介于两者之间的，则判断该像素的八个相邻像素中是否存在大于高阈值High的边缘像素，若存在则是边缘像素，否则不是；Step (d): Double thresholding. Set the low threshold Low and high threshold High required by the dual-threshold method to detect and connect edges, and perform dual-thresholding processing on the gradient image. If the gradient magnitude is greater than the high threshold High, it is an edge; if the gradient magnitude is smaller than the low threshold Low, it is not an edge; if the gradient magnitude is between the two, it is judged whether there is an edge greater than the high threshold High among the eight adjacent pixels of the pixel. The edge pixel of , if it exists, it is an edge pixel, otherwise it is not;

步骤(5.2)：对z值从z＝z_min开始，以Δz步长循环执行步骤(5.2.1)和步骤(5.2.2)直到z≥z_max：Step (5.2): Starting from z=z _min for z value, execute step (5.2.1) and step (5.2.2) cyclically with Δz step size until z≥z _max :

步骤(5.2.1)：利用变波速的相位迁移公式计算z+Δz处的k-ω域信号：Step (5.2.1): Calculate the k-ω domain signal at z+Δz using the phase shift formula with variable wave velocity:

S(z+Δz，k，ω)＝S(z，k，ω)α(Δz，k，ω，v_z(x)) (6)S(z+Δz,k,ω)=S(z,k,ω)α(Δz,k,ω, _vz (x)) (6)

其中，α(Δz，k，ω，v_z(x))为变波速的相位迁移量：Among them, α(Δz, k, ω, v _z (x)) is the phase shift of variable wave velocity:

vz(x)的取值如下：The value of vz(x) is as follows:

v₁，v₂分别等于分界线c(x，z)所在的上下两层介质中声波的传播速度；v ₁ and v ₂ are respectively equal to the propagation speed of sound waves in the upper and lower layers of the medium where the boundary line c(x, z) is located;

步骤(5.2.2)：令z←z+Δz，t＝t₀＝0，作二维傅里叶反变换，得到纵向坐标值为z的一行图像：Step (5.2.2): make z←z+Δz, t=t ₀ =0, do two-dimensional Fourier inverse transform, obtain a row of images whose longitudinal coordinate value is z:

$s the s ((z z,, x x,, {t t}_{00})) = = \frac{11}{{44 π π}^{22}} &Integral; &Integral; &Integral; &Integral; S S ((z z,, k k,, ω ω)) {e e}^{ikx ikx} {e e}^{iωt iωt} dkdω dkdω - - - - - - ((77))$

步骤(6)：生成深度方向z_max至z_depth区间(z_max＜z≤z_depth)的图像。令v_z＝v₂，按步骤(4)的方法计算z_max至z_depth区间的各行图像；Step (6): Generating an image in the interval from z _max to z _depth (z _max <z≤z _depth ) in the depth direction. Make v _z =v ₂ , calculate each line image of z _max to z _depth interval by the method for step (4);

步骤(7)：处理后续分界线。对步骤(6)所生成的图像部分(z_max＜z≤z_depth)使用Canny边缘提取算子提取新分界线c(x，z)，令v₁，v₂分别为分界线c(x，z)所在的上下两层介质中声波的传播速度，并令z_min和z_max分别为新分界线c(x，z)在深度方向上的上下边界，形成新的包围盒，如果z_min＝＝z_max＝＝z_depth，则结束；否则，执行步骤(5.2)修正z_min至z_max区间的图像；Step (7): Process the subsequent dividing line. Use the Canny edge extraction operator to extract a new boundary line c(x, z) from the image part (z _max <z≤z _depth ) generated in step (6), let v ₁ and v ₂ be the boundary line c(x, The speed of sound wave propagation in the upper and lower layers of the medium where z) is located, and let z _min and z _max be the upper and lower boundaries of the new dividing line c(x, z) in the depth direction to form a new bounding box, if z _min = =z _max ==z _depth , then end; otherwise, perform step (5.2) to modify the image in the interval from z _min to z _max ;

步骤(8)：如果z_max＝＝z_depth，则结束；否则，重复执行步骤(6)和步骤(7)。Step (8): If z _max ==z _depth , then end; otherwise, repeat step (6) and step (7).

其中，为实现步骤(2)中的二维傅里叶变换(如公式(1))，本发明采用对参数x和t按顺序进行两次一维快速傅里叶变换(1D-FFT)的方法，以提高计算速度，即公式(1)实现为：Wherein, in order to realize the two-dimensional Fourier transform (as formula (1)) in the step (2), the present invention adopts the method of carrying out twice one-dimensional fast Fourier transform (1D-FFT) to parameter x and t in order method, to improve the calculation speed, that is, the formula (1) is realized as:

步骤(3.2)和步骤(4.2)中需作二维傅里叶反变换并取t＝t₀＝0的值作为图像数据，则可以先令t＝t₀＝0，并对ω作积分，然后再对k作积分，以减少第二维积分运算的数据量。那么公式(3)(5)化为：In step (3.2) and step (4.2), it is necessary to do two-dimensional inverse Fourier transform and get the value of t=t ₀ =0 as the image data, then t=t ₀ =0 can be ordered first, and ω is integrated, Then k is integrated to reduce the amount of data in the second dimension integral operation. Then formula (3)(5) turns into:

$s the s ((z z,, x x,, {t t}_{00})) = = \frac{11}{22 π π} &Integral; &Integral; {{\frac{11}{22 π π} &Integral; &Integral; S S ((z z,, k k,, ω ω)) dω dω}} {e e}^{ikx ikx} dk dk$

其中，对ω的积分

可以利用数值积分方法(如高斯Gaussian积分)快速求得。而对k的傅里叶反变换仍通过一维快速傅里叶变换(1D-FFT)来实现，只是先将数值积分得到的数据取共轭([Integration_ω(S(z，k，ω))]*)，然后直接利用1D-FFT子程序进行运算，再将运算结果取一次共轭，并乘以1/N_x，即可得到。N_x为换能器移动步数，即N_x＝1+L_x/Δx，其中L_x为待测物体在x方向的长度。因此，公式(3)(5)的具体实现过程表示为：where the integral over ω

It can be quickly obtained by numerical integration methods (such as Gaussian integral). The inverse Fourier transform of k is still realized by one-dimensional fast Fourier transform (1D-FFT), but the data obtained by numerical integration is firstly conjugated ([Integration _ω (S(z, k, ω) )]*), and then directly use the 1D-FFT subroutine to perform the operation, and then take the conjugate of the operation result once, and multiply it by 1/N _x to get it. N _x is the number of moving steps of the transducer, that is, N _x =1+L _x /Δx, wherein L _x is the length of the object to be measured in the x direction. Therefore, the specific implementation process of formula (3)(5) is expressed as:

$s the s ((z z,, x x,, {t t}_{00})) = = \frac{11}{{N N}_{x x}} {{11 D D. - - {FFT FFT}_{k k} (([[Integratio Integration {n no}_{ω ω} ((S S ((z z,, k k,, ω ω)))) {]]}^{* *})) {}}}^{* *}$

相位迁移量α(Δz，k，ω，v_z)需要针对每一个速度v_z，分别对不同的k、ω作计算，其结果预先保存在三维数组中。其中，k、ω的取值分别为：The phase shift amount α(Δz, k, ω, v _z ) needs to be calculated for each velocity v _z for different k and ω, and the results are stored in a three-dimensional array in advance. Among them, the values of k and ω are respectively:

0≤n_k≤N_x-1，且n_k为整数；

0≤n _k ≤N _x -1, and n _k is an integer;

0≤n_ω≤N_t-1，且n_ω为整数；

0≤n _ω ≤N _t -1, and n _ω is an integer;

其中，Δx为换能器沿x方向移动的步长，f_s为采样频率(其值为系统预设)，N_x＝1+L_x/Δx为换能器沿x方向的探测点数，L_x/Δx为换能器移动步数，L_x为被测物体的长度(x方向)，N_t为换能器在每一探测点处对回波信号的采样次数。由于相位迁移超声成像技术将其工作模型视为爆炸物反射模型，仅考虑声波从反射物处向上(沿z轴负方向)传播的部分，所以在相位迁移量α(Δz，k，ω，v_z)的计算中仅取ω＜0部分的值，而令ω≥0部分的相位迁移量为0，同时若b＜0，则取

由ω取值的周期性知，当要求ω＜0时，计算中可取N_t/2≤n_ω＜N_t部分。Among them, Δx is the step size of the transducer moving along the x direction, f _s is the sampling frequency (its value is the system preset), N _x =1+L _x /Δx is the number of detection points of the transducer along the x direction, L _x /Δx is the number of moving steps of the transducer, L _x is the length of the object to be measured (in the x direction), and N _t is the number of times the transducer samples the echo signal at each detection point. Since the phase-shift ultrasonic imaging technology regards its working model as the reflection model of explosives, and only considers the part of the sound wave propagating upward from the reflector (along the negative direction of the z-axis), the phase shift α(Δz, k, ω, v In the calculation of _z ), only the value of the part of ω<0 is taken, and the phase shift of the part of ω≥0 is 0, and if b<0, then take

From the periodicity of the value of ω, when ω<0 is required, the part N _t /2≤n _ω <N _t can be used in the calculation.

步骤(5.2)中，变波速相位迁移量α(Δz，k，ω，v_z(x))是x的函数，波速v_z非常数，所以步骤(5.2.1)不能再直接计算，而需要合并到步骤(5.2.2)的傅里叶反变换中，即In step (5.2), the variable wave velocity phase shift α(Δz, k, ω, v _z (x)) is a function of x, and the wave velocity v _z is non-constant, so step (5.2.1) can no longer be directly calculated, but needs Incorporated into the inverse Fourier transform of step (5.2.2), namely

$s the s ((z z + + Δz Δz,, x x,, {t t}_{00})) = = \frac{11}{22 π π} &Integral; &Integral; {{\frac{11}{22 π π} &Integral; &Integral; [[S S ((z z,, k k,, ω ω)) α α ((Δz Δz,, k k,, ω ω,, {v v}_{z z} ((x x))))]] {e e}^{ikx ikx} dk dk}} {e e}^{iωt iωt} dω dω - - - - - - ((88))$

公式(8)中的外层积分可以利用数值积分方法(如高斯Gaussian积分)求得，而内层积分运算(如式(9))则需通过修改离散傅里叶反变换计算公式求得(如式(10))，The outer integral in formula (8) can be obtained by numerical integration method (such as Gaussian integral), while the inner integral operation (such as formula (9)) needs to be obtained by modifying the calculation formula of inverse discrete Fourier transform ( Such as formula (10)),

$s the s ((z z + + Δz Δz,, x x,, ω ω)) = = \frac{11}{22 π π} &Integral; &Integral; [[S S ((z z,, k k,, ω ω)) α α ((Δz Δz,, k k,, ω ω,, {v v}_{z z} ((x x))))]] {e e}^{ikx ikx} dk dk - - - - - - ((99))$

$S S ((z z + + Δz Δz,, x x = = {n no}_{x x} \cdot &Center Dot; Δx Δx,, ω ω)) = = \frac{11}{{N N}_{x x}} {Σ Σ}_{{n no}_{k k} = = 00}^{{N N}_{x x} - - 11} [[S S ((z z,, k k = = \frac{{22 πn πn}_{k k}}{{N N}_{x x} Δx Δx},, ω ω)) α α ((Δz Δz,, k k = = \frac{{22 πn πn}_{k k}}{{N N}_{x x} Δx Δx},, ω ω,, {v v}_{z z} (({n no}_{x x}))))]] {e e}^{i i \frac{22 π π}{{N N}_{x x}} {n no}_{k k} {n no}_{x x}} - - - - - - ((1010))$

式(10)中，n_x为换能器所在的探测点在x方向上的序号，取值范围同于n_k，而α(Δz，k，ω，v(_zx))是自变量为n_k和n_x的函数，不能直接使用离散傅里叶反变换(IDFT)或者快速傅里叶反变换(IFFT)方法。因此，本发明中使用N次单输出快速傅里叶反变换方法(NsIFFTPruning)。该方法基于快速傅里叶变换中蝶形运算的思想，但每一次变换仅输出一个变换域中的值。NsIFFTPruning具体算法如下：In formula (10), n _x is the serial number of the detection point where the transducer is located in the x direction, and the value range is the same as n _k , and α(Δz, k, ω, v( _z x)) is the independent variable The functions of n _k and n _x cannot directly use the inverse discrete Fourier transform (IDFT) or inverse fast Fourier transform (IFFT) methods. Therefore, in the present invention, an N-times single-output inverse fast Fourier transform method (NsIFFTrunning) is used. This method is based on the idea of butterfly operation in fast Fourier transform, but each transform only outputs a value in the transform domain. The specific algorithm of NsIFFTrunning is as follows:

For n_x＝0 to N_x-1For n _x ＝0 to N _x -1

For n_k＝0 to N_x-1For n _k ＝0 to N _x -1

$G G (({n no}_{x x},, {n no}_{k k})) = = S S ((z z,, k k = = \frac{{22 πn πn}_{k k}}{{N N}_{x x} Δx Δx},, ω ω)) α α ((Δz Δz,, k k = = \frac{{22 πn πn}_{k k}}{{N N}_{x x} Δx Δx},, ω ω,, {v v}_{z z} (({n no}_{x x}))))$

endend

S(z+Δz，x＝n_x·Δx，ω)＝IFFTPruning(G，N_x，n_x)S(z+Δz, x=n _x Δx, ω) = IFFTPruning(G, N _x , n _x )

endend

其中，当n_x固定时，函数IFFTPruning(G，N_x，n_x)计算G(n_x，n_k)，n_k＝0，1，…，N_x-1序列的离散傅里叶反变换，但只取n_x这一时域点的值，即Wherein, when n _x is fixed, the function IFFTPruning (G, N _x , n _x ) calculates G(n _x , n _k ), n _k =0, 1,..., the inverse discrete Fourier transform of N _x -1 sequence , but only take the value of n _x this time domain point, namely

$S S ((z z + + Δz Δz,, x x = = {n no}_{x x} \cdot &Center Dot; Δx Δx,, ω ω)) = = \frac{11}{{N N}_{x x}} {Σ Σ}_{{n no}_{k k} = = 00}^{{N N}_{x x} - - 11} G G (({n no}_{x x},, {n no}_{k k})) {e e}^{i i \frac{22 π π}{{N N}_{x x}} {n no}_{k k} {n no}_{x x}} - - - - - - ((1111))$

式(11)可以使用IDFT或者IFFT算法来计算，但与一般IDFT问题不同的是：式(11)仅需计算G(n_x，n_k)在n_x＝h这一时域点处的值，如果直接使用IDFT或者IFFT算法，则必然存在大量的冗余计算，例如，在图10中，若只需计算n_x＝2的点，则图中虚线部分的蝶形运算是冗余的。为了减少计算量，本发明使用裁剪输出的FFT算法(IFFTPruning)。该算法仍基于FFT算法，采用蝶形运算方法，其原理如下：根据FFT算法中蝶形分组方法，将计算分为H个步骤(Stage)，其中H满足2^H＝N_x。在步骤h(h＝1，...，H)中，将N_x点输入序列G(n_x，n_k)，n_k＝0，1，…，N_x-1平均分为2^H-h个组，每个组中分别从0开始递增1编号b_h(b_h＝0，1，...，2^h-1)，在每组蝶形运算中除了编号为b_h＝b_h+1mod 2^h-1的点外，其余的点都为冗余的。在每一个计算步骤中，仅计算各组蝶形运算中编号为b_h的点。如图10所示，N_x＝8，n_x＝2时，步骤(1)中计算每组蝶形运算中编号为b₁＝0的点，步骤(2)中计算每组蝶形运算中编号为b₂＝2的点，步骤(3)中计算每组蝶形运算中编号为b₃＝2的点.Equation (11) can be calculated using IDFT or IFFT algorithm, but different from the general IDFT problem: Equation (11) only needs to calculate the value of G(n _x , n _k ) at the time domain point of n _x =h, If the IDFT or IFFT algorithm is directly used, there must be a large amount of redundant calculations. For example, in FIG. 10, if only the points of n _x =2 need to be calculated, the butterfly operation in the dotted line in the figure is redundant. In order to reduce the amount of calculation, the present invention uses the FFT algorithm (IFFTrunning) of clipped output. The algorithm is still based on the FFT algorithm and adopts the butterfly operation method. The principle is as follows: According to the butterfly grouping method in the FFT algorithm, the calculation is divided into H steps (Stage), where H satisfies 2 ^H =N _x . In step h (h=1,...,H), the N _x point input sequence G(n _x , n _k ), n _k =0, 1,..., N _x -1 is equally divided into 2 ^Hh group, each group starts from 0 and increments number b _h by 1 (b _h =0, 1, ..., 2 ^h -1), except for the number b _h = b _h+1 in each group of butterfly operations Except for the point of mod 2 ^h-1 , the remaining points are redundant. In each calculation step, only the points numbered b _h in each group of butterfly operations are calculated. As shown in Figure 10, when N _x =8, n _x =2, the point numbered b 1 =0 in each group of butterfly operations is calculated in step (1), and the point numbered b ₁ =0 in each group of butterfly operations is calculated in step (2). For the point numbered b ₂ =2, calculate the point numbered b ₃ =2 in each group of butterfly operations in step (3).

根据前述成像步骤计算得到的图像宽度为N_x，而非待测物体的宽L_x，因此，需要在成像计算阶段的最后一步进行x向插值运算。具体过程为：在x向每相邻两个点间线性插入inst＝Δx/Δz-1个点。The width of the image calculated according to the aforementioned imaging steps is N _x , not the width L _x of the object to be measured. Therefore, an x-direction interpolation operation needs to be performed in the last step of the imaging calculation stage. The specific process is: linearly insert inst=Δx/Δz-1 points between every two adjacent points in the x direction.

图像显示即将成像计算阶段得到的二维图像数据显示在显示器上，可根据需要显示灰度图像或彩色图像。Image display is to display the two-dimensional image data obtained in the imaging calculation stage on the monitor, and grayscale images or color images can be displayed as required.

对于待测对象是规则分层物体的情况，成像步骤(5.1)中使用Canny边缘提取算子提取的分界线c(x，z)为直线，成像计算步骤不变。For the case where the object to be measured is a regular layered object, the boundary line c(x, z) extracted by using the Canny edge extraction operator in the imaging step (5.1) is a straight line, and the imaging calculation step remains unchanged.

对于待测对象中仅含单一介质的情况(说明为接触式超声测量方式，否则应视为规则分层物体)，则z_min＝＝z_max＝＝z_depth，成像计算执行完步骤(3)后即停止。For the case where the object to be measured contains only a single medium (it is described as a contact ultrasonic measurement method, otherwise it should be regarded as a regular layered object), then z _min ==z _max ==z _depth , and the imaging calculation completes step (3) Then stop.

Claims

1. The multilayer object nondestructive testing ultrasonic imaging method based on variable wave speed phase shift is characterized by sequentially comprising the following steps of:

step (1): constructing a system which is composed of a computer, an ultrasonic transducer, a set of positioning controllers and an analog-to-digital converter and is used for carrying out nondestructive ultrasonic imaging on a vertical section formed in two directions of depth and horizontal of a three-layer heterogeneous object with heterogeneous media in the transverse direction and the longitudinal direction based on variable wave speed phase shift, wherein:

the ultrasonic transducer is provided with: the pulse signal input end is connected with the output end of the positioning controller, the input end of the positioning controller is connected with the corresponding positioning control signal output end of the computer, and the ultrasonic transducer is also provided with: and the output end of the analog-to-digital converter is connected with the echo sampling signal input end of the computer. The transducer is controlled by the positioning controller to move at a fixed speed of 1 step/ms on the surface of a three-layer heterogeneous object, the positioning controller is a transmission device for controlling the moving position of the ultrasonic transducer, and parameters of the transmission device are input by the computer,

the horizontal length of the three-layer heterogeneous object along the x-axis direction is L, the three-layer heterogeneous object is divided into L/delta x intervals, delta x is the interval length, and the interval length is from a point 0 to an end point N of the ultrasonic transducer on the x-axis_x1 step size per move, N_xThe point reached by the ultrasonic transducer in each movement is called probe point, and the total number of the probe points is N_xNumber n_x＝0，1，...，N_x-1, said positioning controller generating a TTL pulse at each detection point triggering said ultrasound transducer to transmit an excitation pulse in a depth direction z perpendicular to the x-axis of the three-layered heterogeneous object, then the transducer switches to a receiving mode and starts timing to receive the echo signal reflected from the object to be measured, said analog-to-digital converter couples said transducer at a detection point n_xProcessing the received echo signal by N_tSub-sampling and storing in computer, sampling number n_t＝0，1，...，N_t-1, sampling frequency f_sThe value is predetermined for the analog-to-digital converter, and s (z, x, t) is represented by x ═ n on the abscissa_xThe echo signal received at Δ x is n at t_t/f_sSampling values of the moments;

step (2): the computer starts from n_tStarting to read the probe point n sequentially when the probe point n is 0_xThe sampling value at 0 and the sampling sequence number of the first fluctuation of the sampling value is recorded

Subscript 0 denotes the probe number, and then repeatsThe above steps read n in sequence_x＝1，...，N_x-1 sampling values at the probe points and recording the number of first fluctuations of the sampling values at the points

To obtain N_xN in total on each probe point_xNumber of samples fluctuating for the first time

From

Select the smallest serial number

And calculating the depth coordinate value corresponding to the sequence number in the z direction of the three-layer heterogeneous object

And from

Select the largest serial number

Wherein v is₁Respectively making two horizontal lines parallel to the x axis according to the two depth coordinate values to surround a boundary between the first layer of medium and the second layer of medium to form a first bounding box which is positioned on the x-z longitudinal section of the three-layer heterogeneous object and simultaneously spans the first layer of medium and the second layer of medium;

and (3): the sampled values s (z, x, t) at the detection points are entered along the time axis t one by oneRow one-dimensional fast Fourier transform 1D-FFT_t(s (z, x, t)), and taking the initial value z of z₀Obtaining an x-omega domain two-dimensional signal:

S(z₀，x，ω)＝1D-FFT_t(s(z₀，x，t))

wherein x is the serial number n of the detection point_xAt horizontal coordinate value, where ω is sampling angular frequency, and then, for the above x- ω domain two-dimensional signal S (z)₀X, ω) are sequentially subjected to one-dimensional fast fourier transform 1D-FFT along a horizontal axis x_x(S(z₀X, ω)) to yield z₀K-omega domain two-dimensional spectrum S (z)₀，k，ω)：

S(z₀，k，ω)＝1D-FFT_x(1D-FFT_t(s(z₀，x，t)))

Wherein k is the x-component of the wavenumber vector with the sequence number n_kThe wavenumber vector represents a spatial gradient of the wave phase,

0≤n_k≤N_x-1, and n_kIs an integer;

0≤n_ω≤N_t-1, and n_ωIs an integer;

wherein n is_ωThe value of the sampling angular frequency omega is 0 to (N)_t-1), take t ═ t₀One-dimensional sampled data s (z) of 0₀，x，t₀) I.e. the z-th ═ z on the vertical section of the three-layer heterogeneous object₀Pixel values of pixel points in 0 row;

and (4): generating a z-direction homogeneous region z by the following steps₀＜z≤z_minThe inner images of all rows formed by pixel points of all rows respectively:

step (4.1): calculating the profile in z using the following phase shift equation_minTwo-dimensional spectrum S (z) of_min，k，ω)：

Wherein,

wherein i is an imaginary unit, z₀0, ω is the sampling angular frequency, v_zThe wave velocity of the sound wave in the z-th row of medium of the measured object in the depth direction is equal to or less than z_minWhen, v_z＝v₁，v₁Is the wave velocity of the acoustic wave in the first layer of medium, and Δ z is the separation distance between two adjacent pixels in the depth direction z, and α (Δ z, k, ω, v)_z) Is the amount of phase migration at a constant wave velocity,

step (4.2): the following steps are performed in sequence:

step (4.2.1): get t ═ t₀(iv) 0, for S (z) obtained in the step (4.1)_minK, ω) is taken as the gaussian integral of the independent variable ω,

step (4.2.2): to pairConjugating the Gaussian integral result of the step (4.2.1), performing one-dimensional fast Fourier transform on k, conjugating the calculation result, and multiplying the result by 1/N_xObtaining the coordinate value z in the depth direction in the first layer medium_minThe image formed by the row of pixel points:

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>min</mi> </msub> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msub> <mi>N</mi> <mi>x</mi> </msub> </mfrac> <mo>{</mo> <mn>1</mn> <mi>D</mi> <mo>-</mo> <msub> <mi>FFT</mi> <mi>k</mi> </msub> <mo>{</mo> <mo>[</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>π</mi> </mrow> </mfrac> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>π</mi> </mrow> <mi>π</mi> </msubsup> <mi>S</mi> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>min</mi> </msub> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <mi>dω</mi> <msup> <mo>]</mo> <mo>*</mo> </msup> <mo>)</mo> <msup> <mo>}</mo> <mo>*</mo> </msup> </mrow> </math>

in the z-direction the first layer homogeneous zone z₀＜z≤z_minAnd in the inner step, the images of all rows are the same:

s(z，x，t₀)＝s(z_min，x，t₀)，z₀＜z≤z_min；

and (5): taking v in step (4)_z＝v₁And circularly performing the following steps (5.1) and (5.2) by taking the step delta z as the step size for the z value until z is more than or equal to z_maxGenerating a first enclosureDirection of internal depth z_min＜z≤z_maxEach line of image formed by each line of pixel points in the interval respectively comprises the following steps:

step (5.1): calculating a two-dimensional frequency spectrum at z + Δ z using the phase shift formula at the constant wave velocity, z_min＜z+Δz≤z_max：

S(z+Δz，k，ω)＝S(z，k，ω)·α(Δz，k，ω，v_z)，

Step (5.2): the following steps are performed in sequence:

step (5.2.1): get t ═ t₀(ii) the gaussian integral of the argument ω is taken for S (z + Δ z, k, ω) obtained in the step (5.1) described above as 0,

step (5.2.2): conjugating the Gaussian integral result of the step (5.2.1), performing one-dimensional fast Fourier transform on k, conjugating the calculation result and multiplying by 1/N_xObtaining an image formed by a row of pixel points with the coordinate value of z + delta z in the depth direction, z_min＜z+Δz≤z_max：

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>+</mo> <mi>Δz</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msub> <mi>N</mi> <mi>x</mi> </msub> </mfrac> <mo>{</mo> <mn>1</mn> <mi>D</mi> <mo>-</mo> <msub> <mi>FFT</mi> <mi>k</mi> </msub> <mo>{</mo> <mo>[</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>π</mi> </mrow> </mfrac> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>π</mi> </mrow> <mi>π</mi> </msubsup> <mi>S</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>+</mo> <mi>Δz</mi> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <mi>dω</mi> <msup> <mo>]</mo> <mo>*</mo> </msup> <mo>)</mo> <msup> <mo>}</mo> <mo>*</mo> </msup> <mo>;</mo> </mrow> </math>

And (6): correcting the depth direction z in the first bounding box as follows_minTo z_maxInterval z_min＜z≤z_maxTo eliminate errors caused by inhomogeneities between the second layer and the first layer of media:

step (6.1): z obtained in step (5)_minTo z_maxTaking an interval image block as an input quantity, and extracting z by using a Canny operator edge extraction algorithm_min、z_maxThe boundary c between the first layer medium and the second layer medium₁(x，z)，

Step (6.2): from z to z_minInitially until z ≧ z_maxZ obtained in said step (4.1)_minTwo-dimensional spectrum S (z) of_minK, ω) as a starting value, and step (6.2.1) and step (6.2.2) are cyclically executed with Δ z as a step size:

step (6.2.1): and calculating to obtain an x-omega domain two-dimensional signal at a z + delta z position in the first bounding box by using an N-time single-output inverse fast Fourier transform method according to the following formula:

wherein the value of N is equal to N_x，n_x＝0，1，...，N_x-1，n_k＝0，1，...，N_x-1，α(Δz，k，ω，v_z(x) Phase shift amount as a wave velocity:

wherein, the values of k and omega are the same as the step (3), v_z(x) Is the wave velocity of the acoustic wave in the medium at coordinates (x, z):

v₁to be located at the boundary line c₁Speed of propagation, v, of acoustic waves in a first layer of medium above (x, z)₂To be located at the boundary line c₁(x, z) the propagation velocity of the acoustic wave in the second layer medium below (x, z),

the calculation steps are as follows:

step (a): respectively for different n_x，n_k(n_x＝0，1，...，N_x-1；n_k＝0，1，...，N_x-1) calculating

Step (b): for different n_x(n_x＝0，1，...，N_x-1) with G (n)_x，n_k)，n_k＝0，1，...，N_x-1 sequence as input to the clipped output FFT algorithm, for n_kAs G (n)_x，n_k) Inverse discrete Fourier transform, but taking only n_xThe value of this time-domain point is taken as an x- ω domain two-dimensional signal at a coordinate value z + Δ z in the depth direction:

step (6.2.2): get t ═ t₀(iv) when S (z + Δ z, x, ω) obtained in the step (6.2.1) is subjected to gaussian integration of the argument ω to obtain an image formed by one line of pixel points whose coordinate value in the depth direction is z + Δ z, z being equal to 0_min＜z+Δz＜z_max：

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>+</mo> <mi>Δz</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>π</mi> </mrow> </mfrac> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>π</mi> </mrow> <mi>π</mi> </msubsup> <mi>S</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>+</mo> <mi>Δz</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <mi>dω</mi> </mrow> </math>

Meanwhile, storing S (z + delta z, x, omega) and performing one-dimensional Fourier transform to obtain S (z + delta z, k, omega) which is used as an initial value of the next iteration;

and (7): according to the method described in step (5), v is taken_z＝v₂And circularly executing the step (5.1) and the step (5.2) by taking the step deltaz as a step size for the z value until z is more than or equal to z_depthTo generate a depth direction z_max＜z≤z_depthEach row image v formed by each row of pixel points in interval₂Is the propagation velocity of the acoustic wave in the second layer medium, z_depthThe maximum coordinate value of the three-layer heterogeneous object in the depth direction z is represented as a system preset value;

and (8): z obtained in step (7)_maxTo z_depthTaking an interval image block as an input quantity, and extracting z by using a Canny operator edge extraction algorithm_max、z_depthThe boundary c between the second layer medium and the third layer medium₂(x, z) determining the boundary line c₂(x, z) minimum coordinate value z 'in depth direction z'_minAnd maximum coordinate value z'_maxRespectively through said z'_minAnd z'_maxTwo depth coordinate values are used for making two horizontal lines parallel to the x axis, and a boundary between the second layer of medium and the third layer of medium is surrounded to form a second surrounding box which is positioned on the x-z vertical section of the three-layer heterogeneous object and simultaneously spans the second layer of medium and the third layer of medium;

and (9): correcting the depth direction z 'in the second surrounding box according to the method of the step (6.2)'_minTo z'_maxZone z'_min＜z≤z′_maxThe images of the rows in the image table are used for eliminating errors caused by the heterogeneity between the third layer medium and the second layer medium;

step (10): according to the method described in step (5), v is taken_z＝v₃Z is z'_maxThe step (5.1) and the step (5.2) are circularly executed by taking the step delta z as a step size until the step z is more than or equal to z_depthGenerating the parts except the second bounding box in the third layer of mediumEach row image formed by row pixel points, v₃Is the propagation speed of the acoustic wave in the third layer medium.