Ultrasonic TOFD Imaging Enhancement Technology for Pipeline Girth Welds Testing via Time Domain Sparse Deconvolution and Frequency Domain Synthetic Aperture Focusing

1. Introduction

Ultrasonic TOFD imaging, as an important imaging method for the pipeline girth welds nondestructive testing, uses diffraction time difference as the basic parameter for quantitatively detecting internal defects . It complements with the phased array zonal ultrasonic imaging method and can provide more precise defect distribution of girth welds. TOFD has become one of the most important methods for quantitative evaluation of weld quality at present, especially for the quantitative characterization of the size and burial depth of vertical crack defects inside the weld seam, it exhibits obvious technical advantages. However, this imaging method has always been affected by signal oscillations caused by the limited frequency band of the transducer and the inherent diffraction effect of defective echo scattering, resulting in unsatisfactory longitudinal and transverse resolutions, the existence of detection blind area, and insufficient detection accuracy. In addition, the real-time limitations of its time-domain imaging method seriously affect the detection performance and online application ability of this method. For this reason, classic deconvolution techniques such as parameterized inverse filtering [1,2], high-order spectral analysis (cepstral theory) [3], and L2 regularized norm reconstruction [4,5] have been applied, effectively improving the longitudinal resolution of ultrasound TOFD imaging.Nevertheless, there are still technical pain points such as low signal-to-noise ratio and low contrast ratio, and further requirements of resolution enhancement.

For improving lateral resolution, time-domain SAFT technology is mainly used [6,7] , which calculates the wave arrival delay between the transducer and the defect, and reconstructs the defect imaging target data through the delay and sum rule. The calculation process of this method is simple, but it requires point-to-point calculation of the imaging area, with limited accuracy and efficiency. TOFD has become one of the most important methods for quantitative evaluation of weld quality at present, especially for the quantitative characterization of the size and burial depth of vertical crack defects inside the weld seam, it exhibits obvious technical advantages [8,9]. However, this numerical reconstruction method requires repeated iterative computation, which brings heavy computational burden. Recently, frequency-domain SAFT methods, also known as wavenumber-domain SAFT, have been widely studied due to their high-quality imaging result and high computational efficiency. The performance of frequency-domain SAFT have been verified in different modes, such as monostatic ultrasonic imaging [10,11], total focusing method ultrasonic imaging [12,13], compound plane wave ultrasonic imaging [14,15], etc. However, this method is established on the foundation of specific ultrasonic imaging model, which is still not be studied in TOFD ultrasonic imaging.

Based on the above background, it is proposed to conduct research on online high-resolution ultrasound TOFD imaging technology for the pipeline girth welds. The proposed method improves the TOFD image resolution in both longitudinal and transverse directions with an extreme low time costs, providing a potential to realize real-time quantitative TOFD detection. According to the simulation and experimental results, the proposed method can outline the trend of defects with very fine stripes, and the transverse resolution of image is reduced about 35%. Also, the time cost of the proposed method is 0.96 second and much less than the TOFD data acquisition time, showing the potential for real-time inspection.

2. Theoretical Expression of the General Model for Ultrasound TOFD Imaging Detection Signals

Ultrasound TOFD imaging detection mainly includes five processes from emission to reception, which are ultrasonic transducer excitation, incident propagation, defect response, reflection propagation, and ultrasonic transducer reception respectively, as shown in Figure 1. According to the classical linear system theory, the general theoretical expression for the ultrasound TOFD detection signal is obtained as:

$$y(\widehat{x},t)=e_t\left(t\right)⨂\left[\int \int g_t(t,x,z)⨂s(x,z)⨂g_b(t,x,z)dxdz\right]⨂e_r\left(t\right)$$

(1)

Among them, $y(\widehat{x},t)$ represents the recieved TOFD signal, $e_t\left(t\right)$ and $e_r\left(t\right)$ represent the exciting response function and recieving response function of the ultrasonic transducers respectively, $g_t(t,x,z)$ and $g_b(t,x,z)$ represent the time-spatial response functions of the incident propagation path and the received propagation path, $s(x,z)$ represents the spatial distribution function which used to characterize the reflection response of the defect. Where, $\widehat{x}$ represents the detection position of the TOFD signals and $x,z$ denote the spatial position of the specific point in ROI. $s(x,z)$ is also the reconstruction goal of high-resolution ultrasound TOFD imaging, and ⨂ represents the convolution operator.

Although all five processes mentioned above are correlated using convolution, their convolution dimension and method are different. $e_t\left(t\right)$ and $e_r\left(t\right)$ modulate signals in the time domain, and its convolution occurs in the time domain, mainly manifested in its effect on signal bandwidth. However, $g_t(t,x,z)$ and $g_b(t,x,z)$ are convolutions in the time spatial domain, which mainly represent the diffraction effect of ultrasound waves in physical space. Furthermore, Equation (1) can be separated into two separate parts and express it in the spatial frequency domain as:

$$y(\widehat{x},\omega )=[e_t\left(\omega \right)·e_r\left(\omega \right)]·p(\widehat{x},\omega )$$

(2)

$$p(\widehat{x},\omega )=\int \int G_t(\omega ,x-\widehat{x},z)·s(x,z)·G_b(\omega ,x-\widehat{x},z)dxdz$$

(3)

As shown in Equation (2), the reflection coefficient function $p(\widehat{x},\omega )$ is separated out as an intermediate. The above two equations are the two problems addressed in this paper, where $p(\widehat{x},\omega )$ is calculated from the measured data $y(\widehat{x},\omega )$, which mainly eliminates the problem of signal bandwidth degradation caused by the transducer response and improves the longitudinal resolution of TOFD imaging. In contrast, the use of $p(\widehat{x},\omega )$ to solve for $s(x,z)$ eliminates diffraction effects and improves the lateral resolution of TOFD imaging. Of course the exact form of the solution will be adapted to the actual situation, but it is fundamentally a temporal and spatial enhancement of the image resolution in terms of the signal bandwidth and diffraction problem.

3. High Longitudinal Resolution Ultrasound TOFD Imaging Based on Sparse Deconvolution

3.1. Sparse Deconvolution Model for Ultrasound TOFD Signals

Due to the limitation of the bandwidth of the transducer, and the transmitted and received ultrasound waves have multiple periodic oscillations,resulting in decrement in the longitudinal (i.e. depth) resolution of TOFD imaging. These periodic oscillations can be modeled as the convolution of transducer’s the emission and reception electronic pulse responses which are irrelevant to physical space, the effect of spatial dimensions can be ignored when eliminating transducer response effectiveness, at this point, $p(\widehat{x},\omega )$ and $y(\widehat{x},\omega )$ can be treated as $p\left(\omega \right)$ and $y\left(\omega \right)$. Here, we use the time-domain one-dimensional discrete convolutional model to represent Equation (2) in time domain:

$$\mathit{y}=\mathit{e}⨂\mathit{p}+\mathit{n}$$

(4)

In the formula, $\mathit{y}={[y\left(0\right)y\left(1\right),...,y(M-1)]}^T$ represent the recieved TOFD signal, and $\mathit{e}={[e\left(0\right)e\left(1\right),...,e(M-1)]}^T$ represent the overall response of the transducer’s emission and reception of electronic pulses, which can be expressed as $e\left(t\right)=e_t\left(t\right)⨂e_r\left(t\right)$. $\mathit{p}={[p\left(0\right)p\left(1\right),...,p(M-1)]}^T$ represent the reflection coefficient function of the acoustic reflector, and $\mathit{n}={[n\left(0\right)n\left(1\right),...,n(M-1)]}^T$ represent the noise.It should be noted that the above convolution operation constructs the incident reference signal with zero as the center reference point, and intercepts the first M vector dimensions after convolution as the representation of signal $\mathit{y}$.Therefore, the vector dimensions before and after convolution operation are consistent. Equation (4) is a standard linear model of ultrasound imaging. Since only the defects reflect the ultrasound waves in TOFD imaging system, the reflection coefficient function p could be used to describe the defect information. Such a linear model also used for different ultrasound imaging techniques [16,17,18]. The electronic pulse response $\mathit{e}$ can be established through the testing frequency response characteristics of various electronic devices, or a set of large plane ultrasonic echo signals can be directly selected for measurement, because the large plane echo signals already approximate the information of ultrasonic signal periodic oscillation.

Obviously, the reflection coefficient function $\mathit{p}$ contain all the defect information of the object to be tested, and is not limited by the ultrasonic signal bandwidth with extremely high longitudinal resolution. Therefore, it is necessary to obtain the reflection coefficient function $\mathit{p}$ of the object to be tested, and solve it according to Equations (4) with known $\mathit{y}$ and $\mathit{e}$. This process is called deconvolution. Due to the presence of noise in actual signals, there is not a strict convolutional relationship between TOFD signal, electronic pulse response, and reflection coefficient function. Directly performing deconvolution can easily lead to ill conditioned solution. To avoid the impact of noise on the solution results, the deconvolution process will be transformed into a linear equation solving approach. To construct the linear equation of Equation (4) , the Fourier transform matrix is introduced to perform the Fourier transform:

$$\mathit{E}=\mathit{F}\mathit{x}$$

(5)

In the formula, $\mathit{E}$ is defined as the spectrum of electronic pulse response, while the Fourier transform matrix F is defined as:

$$\mathit{F}=\left[\begin{array}{cccc}1& 1& ...& 1\\ 1& \omega & ...& 1\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& {\omega }^{M-1}& ...& {\omega }^{(M-1)(M-1)}\end{array}\right],\omega =e^{-j2\pi /M}$$

(6)

Then, by diagonalizing the frequency spectrum $\mathit{E}$ of electron pulse response, a Toeplitz cyclic matrix $\mathit{C}$ can be constructed as:

$$\mathit{C}={\displaystyle \frac{1}{M}}{\mathit{F}}^{†}diag\left(\mathit{E}\right)\mathit{F}$$

(7)

In the formula, ${\displaystyle \frac{1}{M}}$ is used to balance the energy changes during discrete Fourier calculations, and † is a conjugate operator. The Toplitz cycle matrix has generalized diagonality, which can be seen as the spatial translation of the electron pulse response $\mathit{e}$ in the matrix $\mathit{C}$. Obviously, multiplying vector $\mathit{p}$ and matrix $\mathit{C}$ yields the same result as performing truncated and equal length sum e convolution [5]. Therefore, the linear equation of Equation (4) is expressed as

$$\mathit{y}=\mathit{C}\mathit{p}+\mathit{n}$$

(8)

In the formula above, $\mathit{C}$ is a square matrix ${\mathbb{R}}^{M\times M}$. According to the least squares method, the solution in Equation (5) can be expressed as:

$$\underset{\mathit{p}}{min}\left|\right|\mathit{y}-{\mathit{C}\mathit{p}\left|\right|}_2^2$$

(9)

where $\left|\right|·{\left|\right|}_2^2$ represents the ${\ell }_2$ norm.

Theoretically speaking, when the noise level is not high or the condition number of the convolutional matrix $\mathit{C}$ is not large, using the least squares method can suppress noise partly and solve out the reflection coefficient function $\mathit{p}$ . The low noise level here is mainly manifested as signal-to-noise ratio greater than 0 dB, even higher than 12 dB. The condition number of convolutional matrix $\mathit{C}$ can be approximately understood in practical applications as the ultimate expected limit resolution of electronic pulse response. When constructing a convolutional matrix $\mathit{C}$, the size of the translation interval depends on the the ultimate limit resolution. When the translation interval increases and the correlation between row and column vectors goes down to half of its own energy, the condition number satisfies the requirement of stable solvability. However, in practical applications, if the translation interval directly takes the discrete sampling interval, it will cause the correlation between row and column vectors to be much greater than half of its own energy, leading to ill-conditioned in equation solving. Meanwhile, Equation (9) use the ${\ell }_2$ norm as the optimization objective, which is similar to a Gaussian distribution, resulting in excessively smooth solution results and limited improvement in longitudinal resolution.

Usually, in the application environment of ultrasonic non-destructive testing, the probability of defects appearing is relatively very low, and their distribution has large sparsity, that is, the proportion of defect areas in the tested object is very small, independent and discrete. Therefore, a prior assumption of sparse distribution can be introduced into the reflection coefficient function of $\mathit{p}$ the tested object, which is to add a set of sparse constraints to the above optimization objective function and rewrite it as:

$$\underset{\mathit{p}}{min}\left|\right|\mathit{y}-{\mathit{C}\mathit{p}\left|\right|}_2^2{s.t.min\left|\right|\mathit{p}\left|\right|}_0$$

(10)

The ${\ell }_0$ norm in the above formula refers to the number of non zeros in $\mathit{p}$, and this type of constraint method on the model is called regularization. Common regularization methods include ${\ell }_2$ constraints(Tikhonov regularization constraints), and${\ell }_1$ constraints etc. The ${\ell }_0$ constraint term added here has a strong sparse constraint effect on $\mathit{p}$, while trying to avoid ill-conditioned solutions caused by discomforting in the equation solving process. On the basis of correctly solving the defect position, it can effectively suppress the smoothing effect brought by the ${\ell }_2$ norm in Equation (10), and improve the longitudinal resolution of ultrasound TOFD imaging as much as possible.

3.2. Sparse deconvolution optimization method for ultrasound TOFD signals

Based on the linear equation model of Equation (10), multiple optimization methods can be used to solve it, which include common algorithms such as base pursuit (BP) [10,11,12] and match pursuit (MP) [13,14,15] , but their solution approaches are vastly different. The BP algorithm can be seen as a subtraction operation, which projects the echo signal into a new space through the inverse mapping of the matrix $\mathit{C}$. In this space, the defect energy is concentrated and the amplitude is amplified, while the noise energy is dispersed and the amplitude is reduced. Therefore, during the iterative process of solving, the main energy is retained by gradually eliminating the excess small amplitude components, thus achieving the goal of sparse reconstruction. Due to the need for the Jacobian matrix or Hessian matrix of linear equations in the BP calculation process, the solving process still relies on convex optimization theory, which has certain limitations on the model and is not compatible with ${\ell }_0$ norm constrained problems. The MP class methods are mainly based on greedy algorithms [16], similar to incremental operations. After initialization, they are set to zero vectors and continuously add non-zero terms during iteration to approximate the optimization results. Compared to other methods, for ultrasound TOFD detection, solving targets based on ${\ell }_0$ norm constraints is more suitable for optimization using MP class algorithms. Thus, the Orthogonal Matching Pursuit (OMP) algorithm [17] will be used to solve Equations (7) . Before introducing OMP, the conventional MP algorithm will be introduced, and its main process is as follows:

STEP 1: establish dictionary codebook $\mathit{D}$, stimulating atomic libraries $\mathit{B})=\{\varphi \}$, residual $\mathit{R})=\mathit{y})$, sparsity K, iterative counter $k=0$;

STEP 2: if $k⩽K$

①: calculate degree of relevance: ${\mathit{T}}^k={\mathit{D}}^k{\mathit{R}}^k$

②: select optimal atoms: ${\mathit{B}}^k=\{{\mathit{B}}^k,max\left({\mathit{T}}_i^k\right)\}$

③: update residuals: ${\mathit{R}}^k=\mathit{y}-{\mathit{D}}^T{\mathit{B}}^k$

④: next: $k=k+1$, goto STEP2;

, else goto STEP3;

STEP 3: end the loop and get the result: $\mathit{p}={\mathit{B}}^k$ .

In the initialization process, the dictionary codebook $\mathit{D}$ is a sparse spatial mapping matrix. In ultrasonic non-destructive testing, the reflection coefficient function $\mathit{p}$ itself has sparsity, at the same time, the matrix $\mathit{C}$ has cyclic traversal characteristic. Therefore, the codebook $\mathit{D}$ can be directly replaced by convolutional matrices $\mathit{C}$. In correlation calculation, the atom with the highest correlation is calculated by the inner product between the estimated residuals and the codebook. Then, the atom with the highest degree of relevance is selected and entered into the excitation atom library. Meanwhile, the residual value is updated, and this cycle is repeated until the sparsity counter is satisfied.

The dictionary codebook $\mathit{C}$ mentioned here is actually a set of basis . If this set of basis $\Phi$ is orthogonal, then:

$$\mathit{T}={\Phi }^T\mathit{y}={\Phi }^T\Phi \mathit{p}$$

(11)

It can be seen that when the codebook is an orthogonal basis, $\mathit{T}$ can be directly mapped to $\mathit{p}$. Therefore, it is reasonable to represent the main sparse components by selecting several points with the highest degree of relevance in the vector $\mathit{T}$. Typically, the codebooks $\mathit{D}$ are redundant basis, and their redundancy degree depends on the condition number of the convolutional matrix $\mathit{C}$, this leads to non-linear mapping relationship between $\mathit{T}$ and $\mathit{p}$. To eliminate the difference between $\mathit{T}$ and $\mathit{p}$ caused by the correlation between redundant basis, it is necessary to orthogonalize the matching pursuit using OMP (Orthogonal Matching Tracking) method [20].

The orthogonalization of OMP is achieved by eliminating the correlation between vectors within redundant basis, and achieves more accurate matching. The basic process of OMP is similar to MP, with slight differences in correlation calculation [21]. Assuming there are already k atoms, and the normalized basis vectors corresponding to these excited atoms in the codebook are denoted as $\{{\mathit{d}}^1,{\mathit{d}}^2,...,{\mathit{d}}^k\}$. In order to eliminate the impact of non orthogonality of each vector in the redundant basis, the new normalized basis vector is modified during the next round of correlation calculation:

$${\mathit{d}}^{k+1}={\mathit{d}}_o^{k+1}-\sum _{i=1,2,...,k}{\mathit{d}}^i<{\mathit{d}}_o^{k+1},{\mathit{d}}^i>$$

(12)

where ${\mathit{d}}_o^{k+1}$ is the normalized vector basis corresponding to the selected atom in the next iteration in the Equation (12). Obviously, Equation (12) can achieve the normalized vector to be continuously reconstructed during the iteration process, forming an effect similar to orthogonal basis to ensure the accuracy of sparse imaging as much as possible.

The entire OMP calculation has a hyperparameter, which is the sparsity K, and the selection of sparsity K needs to be determined based on the actual object of ultrasonic testing. For ultrasound TOFD, each A-scan signal will contain lateral wave and backwall reflection wave, so at least two atoms of $\mathit{B}$ need to be represented. In the case of defect presence, there will also be diffraction waves at the defect tips. Due to the independent and discrete sparse distribution of defects, it is considered that there is only one set of defects in single A-scan signal at the same time. Therefore, setting the K value to 3 or 4 is more appropriate. It should be noted that in subsequent experiments, the lateral wave is truncated and omitted after straightening, so its sparsity K is correspondingly decreased by 1.

4. High Transverse Resolution Ultrasound TOFD Imaging Based on Frequency Domain SAFT

The scanning imaging process of ultrasonic TOFD for pipeline girth welds is actually a D-scan which moves along the girth direction. The transducers are symmetrically arranged on both sides of the weld seam, and driven by a mechanical device. In an ideal situation, when the acoustic axis of the transmitting transducer passes through a defect, an extremely strong reflected echo will be generated and captured by the receiving transducer for detecting the target defect.However, in practical applications of TOFD, the directionality of sound waves is not like that of lasers, but rather generates diffraction at the defect tips, forming arc-shaped reflected waves. Therefore, even point defects can be detected by probes at multiple locations, forming a "tail" shaped curved defect image. For continuous long strip cracks, there are obvious trailing arcs at both ends of the crack, which can easily lead to inaccurate judgment of defect size and burial depth.

To solve the problem of reduced transverse resolution mentioned above, it is often necessary to use Synthetic Aperture Focusing (SAFT) [16,22] technology to reconstruct images. However, the D-scan image of TOFD is different from conventional two-dimensional B-scan images, its essence is to project three-dimensional spatial information into two-dimensional images, which will inevitably lead to certain information lossness. The most direct explanation example is that D-scan only collects one A-scan line on each axial section and estimates the defect position based on the time delay of reference echo. However, for the same echo time delay, there are multiple possible defect positions on its axial section, and its specific position information cannot be determined. Usually, the SAFT technique for TOFD imaging approximates the central imaging assumption of defects, assuming that the defects are located at the centerline position of the transducer pair. When the defect is assumed to fall on the centerline of the transducer pair, the emission and reception processes of sound waves naturally satisfy the condition of symmetrical sound path, that is, the emission path and reception path are symmetrical about the centerline.

According to Equation (1), the spatial diffraction propagation process of ultrasound TOFD can also be constructed by a convolutional model, which is usually defined in the analytical model as a finite amplitude Delta signal with amplitude information $s(x,z)$. As shown in Equation (3), the convolution model can be simplified by multiplying the transmitting frequency-spatial function $G_t(\omega ,x-\widehat{x},z)$, the amplitude information $s(x,z)$ and the receiving frequency-spatial fucntion $G_b(\omega ,x-\widehat{x},z)$ in frequency domain. In order to associate the angular frequency $\omega$ with the propagating direction of the ultrasonic waves, $\omega$ is replaced with the wavenumber $k=\omega /c$. Thus Equation (3) is rewritten as

$$p(\widehat{x},k)=\int \int G_t(k,x-\widehat{x},z)·s(x,z)·G_b(k,x-\widehat{x},z)dxdz$$

(13)

The imaging reconstruction is the inverse computation process of $s(x,z)$ from $p(x,k)$. Neglecting the multiple reflections of sound waves in propagating process, the transceiver path is symmetry. Thus, the frequency-spatial functions for transmitting and receiving are equalvalent without differences

$$G(k,x-\widehat{x},z)=G_t(k,x-\widehat{x},z)=G_b(k,x-\widehat{x},z)$$

(14)

Since $G(k,x-\widehat{x},z)$ can be regarded as the 2-D free-space Green’s function $G(k,x-\widehat{x},z)=exp\left(ik\sqrt{\left(x-\widehat{x}\right)+z^2}\right)/\sqrt{\left(x-\widehat{x}\right)+z^2}$, Equation (13) is then obtained as

$$\begin{array}{cc}\hfill & p(\widehat{x},k)=\int \int {\displaystyle \frac{exp\left(2·i{k}^{\prime }\sqrt{{\left(x-\widehat{x}\right)}^2+z^2}\right)}{\sqrt{{\left(x-\widehat{x}\right)}^2+z^2}}}·s(x,z)dxdz\hfill \\ \hfill & k^{\prime }=2k={\displaystyle \frac{2\omega }{c}}={\displaystyle \frac{\omega }{c_{erm}}}\hfill \end{array}$$

(15)

In TOFD, the transmitting and receiving probes are symmetrically distributed on both sides of the weld. Therefore, the processes of wave incident onto central weld defects, and wave scattering from the defects and received by receiver are also symmetrical. To simplify the calculation of symmetrical bidirectional wave field propagation, the Explosion Reflection Model (ERM) is introduced. Under ERM, we assume that the sound speed of the propagating medium is half of its real value (cerm=c/2), and the propagation process can be equivalently modeled as ultrasonic waves being autonomously emitted from each pixel point, which are then received by the ultrasonic transducer. The phase information of the sound waves is consistent with the self-emission and reception scenario, and the waveform is equivalent along the time axis. In this model, all scatter points are set to emit sound waves at the same time—referred to as the “explosion” moment—and their intensity is directly proportional to the reflection coefficient of each point. The sound fields produced by these points are emitted simultaneously and received together by the ultrasonic transducer. This simplifies the bidirectional wave propagation and only considers the sound field received from the upward returning probe.

According to Weyl’s integral, the Green’s function in the 2-D free space can be expanded and expressed as:

$${\displaystyle \frac{exp\left(2·i{k}^{\prime }\sqrt{{\left(x-\widehat{x}\right)}^2+z^2}\right)}{\sqrt{{\left(x-\widehat{x}\right)}^2+z^2}}}={\displaystyle \frac{-i}{4\pi }}\int {\displaystyle \frac{exp\left(i{k}_x\left(x-\widehat{x}\right)-i\left|z\right|\sqrt{k^{\prime 2}-k_x^2}\right)}{\sqrt{k^{\prime 2}-k_x^2}}}d{k}_x$$

(16)

Substituting Equation (16) into Equation (15) yields

$$p\left(\widehat{x},k^{\prime }\right)={\displaystyle \frac{-i}{4\pi }}\int {\displaystyle \frac{exp\left(i{k}_{\widehat{x}}\widehat{x}\right)}{\sqrt{k^{\prime 2}-k_{\widehat{x}}^2}}}\left(\int \int exp\left(-i{k}_{\widehat{x}}x-i\left|z\right|\sqrt{k^{\prime 2}-k_{\widehat{x}}^2}\right)·s\left(x,z\right)dxdz\right)d{k}_{\widehat{x}}$$

(17)

It can be seen that the relationship between $p(\widehat{x},k)$ and $s(x,z)$ established by Equation (17) has a complex Fourier decomposition relationship, and here we perform an inverse Fourier transform over $k^{\prime }$. Equation (17) can be rewritten as:

$$p\left(k_{\widehat{x}},k^{\prime }\right)={\displaystyle \frac{-i}{4\pi }}·{\displaystyle \frac{1}{\sqrt{k^{\prime 2}-k_{\widehat{x}}^2}}}\left(\int \int exp\left(-i{k}_{\widehat{x}}x-i\left|z\right|\sqrt{k^{\prime 2}-k_{\widehat{x}}^2}\right)·s\left(x,z\right)dxdz\right)$$

(18)

For convenience, we assume that $k_x=k_{\widehat{x}}$ when the size of ROI remains unchanged in the depth direction. Similarly, the amplitude information $s(x,z)$ can be replaced by the spectrum $S\left(k_x,z\right)$ through a Fourier transform over x

$$P\left(k_{\widehat{x}},k^{\prime }\right)={\displaystyle \frac{-i}{4\pi }}·{\displaystyle \frac{1}{\sqrt{k^{\prime 2}-k_{\widehat{x}}^2}}}\int exp\left(-i\left|z\right|\sqrt{k^{\prime 2}-k_x^2}\right)·S\left(k_x,z\right)dz$$

(19)

According to the analysis based on plane wave model, the acoustic phase $exp\left(-i\left|z\right|\sqrt{k^{\prime 2}-k_x^2}\right)$ has a much larger modulation effect on the signal than $\sqrt{k^{\prime 2}-k_{\widehat{x}}^2}$, and thus $\sqrt{k^{\prime 2}-k_{\widehat{x}}^2}$ can be neglected here. Since the final image can be normalized, the constant coefficient $-i/\left(4\pi \right)$ is also discarded. Equation (19) is simplified as

$$P\left(k_x,k^{\prime }\right)=\int exp\left(-i\left|z\right|\sqrt{k^{\prime 2}-k_x^2}\right)·S\left(k_x,z\right)dz$$

(20)

Equation (20) can be replaced with recursive formula

$$\begin{array}{cc}\hfill & S\left(k_x,z\right)=P\left(k_x,k^{\prime }\right)exp\left(i{k}_{z}z\right)\hfill \\ \hfill & k_z=\sqrt{k^{\prime 2}-k_x^2}\hfill \end{array}$$

(21)

where $k_z$ is the wavenumber along the z-axis. The whole process of wavefield derivation is to calculate $S\left(k_x,z\right)$ and thus estimate the defect signal $s(x,z)$ with the given wavefield $P\left(k_x,k^{\prime }\right)$. However, unlike the conventional B-scan, the deduced depth information z of D-scan is not directly equal to the defect depth, and because of the lateral distance PCS between the transducer’s scanning plane and the imaging plane, the deduced depth should be corrected as follows:

$$z^{\prime }=\sqrt{PC{S}^2+z^2}$$

(22)

Its corresponding derivation formula should be written as

$$S\left(k_x,z\right)=P\left(k_x,k^{\prime }\right)exp\left(i{k}_z\sqrt{PC{S}^2+z^2}\right)$$

(23)

The amplitude information $s(x,z)$ is calculated by performing an inverse Fourier transform over $k_x$. Since the actual signals is band limited, then the final image is obtained by superimosing all the $\left(k_x,z\right)$ values over the wavenumber $k^{\prime }$ as

$$s(x,z)=\int \int S\left(k_x,z\right)exp\left(i{k}_{x}x\right)d{k}_{x}d{k}^{\prime }$$

(24)

In summary, the high lateral resolution is obtained with the following calculation steps

The longitudinal resolution optimized data $\mathit{p}$ is written as $[{\mathit{p}}_1,{\mathit{p}}_2,\cdots ,{\mathit{p}}_K]$, where ${\mathit{p}}_k$ is equal to $\mathit{p}(k,\omega )={[{\mathit{p}}_k\left(0\right),{\mathit{p}}_k\left(0\right),\cdots ,{\mathit{p}}_k(M-1)]}^T$. ${\mathit{p}}_k$ is calculated with the sparse deconvolution for each position as

$$\begin{array}{c}\hfill \underset{{\mathit{p}}_k}{min}{∥{\mathit{y}}_k-{\mathit{C}}_k{\mathit{p}}_k∥}_2^{2}s.t.min{∥{\mathit{p}}_k∥}_0\end{array}$$

(25)

The Fourier transform matrix ${\mathit{F}}_{k\widehat{x}}$ over $\widehat{x}$ is defined as

$${\mathit{F}}_{k\widehat{x}}=\left[\begin{array}{cccc}1& 1& \dots & 1\\ 1& k_{\widehat{x}}& \dots & 1\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& {k_{\widehat{x}}}^{K-1}& \dots & {k_{\widehat{x}}}^{(K-1)(K-1)}\end{array}\right],k_{\widehat{x}}=e^{-i2\pi /K}$$

(26)

The following recursive wavefield derivation at depth z is to left multiply the dispersion relation matrix $\mathit{Z}$

$$\begin{array}{cc}\hfill & \mathit{Z}=\left[\begin{array}{cccc}{k}_z(0,0)& k_z(0,1)& \dots & k_z(0,M-1)\\ k_z(1,0)& k_z(1,1)& \dots & k_z(1,M-1)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ k_z(K-1,0)& k_z(K-1,1)& \dots & k_z(K-1,M-1)\end{array}\right],\hfill \\ \hfill & k_z(k,m)=exp\left\{i\sqrt{{\displaystyle \frac{{\left[\frac{2\pi }{M\Delta t}·\left(-\frac{M}{2}+m\right)\right]}^2}{c_{erm}^2}}-{\left[{\displaystyle \frac{2\pi }{K\Delta x}}·\left(-{\displaystyle \frac{K}{2}}+k\right)\right]}^2}\sqrt{PC{S}^2+z^2}\right\}\hfill \end{array}$$

(27)

The inverse Fourier transform matrix ${\mathit{F}}_x$ over $k_x$ is defined similarly as

$${\mathit{F}}_x=\left[\begin{array}{cccc}1& 1& \dots & 1\\ 1& x& \dots & 1\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& x^{K-1}& \dots & x^{(K-1)(K-1)}\end{array}\right],x=e^{-i2\pi /K}$$

(28)

The superimposition over $\omega$ is to right multiply the constant matrix $\mathit{E}$ as

$$\mathit{E}=\left[\begin{array}{c}1\\ 1\\ .\\ .\\ .\\ 1\end{array}\right]$$

(29)

The calculation of the final image ${\mathit{I}}_z$ at depth z is calculated with the aforementioned matrixed as

$${\mathit{I}}_z=\left\{{\mathit{F}}_x\left[\mathit{Z}\odot \left({\mathit{F}}_{k\widehat{x}}\mathit{p}\right)\right]\right\}\mathit{E}$$

(30)

where ⊙ is the symbol of Hadamard matrix multiplication. The wavefield calculation for ${\mathit{I}}_z$ is then performed for each depth.

The entire computational workflow, depicted in Figure 2, primarily consists of three steps: frequency domain spatial mapping, phase migration, and imaging spatial inverse mapping. In the frequency domain spatial mapping step, a two-dimensional Fourier transform is employed to map the received signal $p(\widehat{x},t)$ into $P\left(k_{\widehat{x}},k^{\prime }\right)$, with a computational complexity of approximately $O(N^{2}log{N}^2)$. The phase migration step involves extrapolating the mapped wavefield layer by layer to gather acoustic field information at different depths. The core computational step of this process is accomplished using the Hadamard product, resulting in a complexity of $O\left(N^2\right)$. The overall computation complexity is $O\left(N^3\right)$ for the case of N imaging depth layers. The final step, imaging spatial inverse mapping, is conducted through the imaging condition formula, with a complexity of $O(N^{2}log{N}^2)$. When combining these three steps, the computation for the entire imaging reconstruction is $O\left(N^3\right)$, which is evidently more efficient than the $O\left(N^4\right)$ required by traditional SAFT.

5. Ultrasound TOFD Imaging with Both High Longitudinal and Lateral Resolutions

Since the process of improving the longitudinal and lateral resolutions is holistic, we unify the two procedures as a whole. The acquired raw data $\mathit{y}$ is denoted as $[{\mathit{y}}_1,{\mathit{y}}_2,\cdots ,{\mathit{y}}_K]$, where ${\mathit{y}}_k$ is equal to $\mathit{y}(k,\omega )={[{\mathit{y}}_k\left(0\right),{\mathit{y}}_k\left(0\right),\cdots ,{\mathit{y}}_k(M-1)]}^T$. First, the optimized data $\mathit{p}$ containing all the scanning positions is calculated in order to improve the longitudinal resolution. Second, the lateral resolution is improved with the frequency domain SAFT and then the final image is achieved. The whole process can be written as

$$\begin{array}{cc}\hfill & {\mathit{I}}_z=\left\{{\mathit{F}}_x\left\{\mathit{Z}\odot \left[{\mathit{F}}_{k\widehat{x}}\left(\underset{\mathit{p}}{min}{∥\mathit{y}-\mathit{Cp}∥}_2^{2}s.t.min{∥\mathit{p}∥}_0\right)\right]\right\}\right\}\mathit{E}\hfill \\ \hfill & \mathit{I}=\left[{\mathit{I}}_1,{\mathit{I}}_2,\cdots ,{\mathit{I}}_M\right]\hfill \end{array}$$

(31)

where M represents the maximum depth in the final image.

6. Simulation and Experimental Research

In order to verify the effectiveness of the above method, simulation and experimental studies are carried out, and the design drawings and physical photographs of the specimen used are shown in Figure 3. The specimen is a flat plate butt-welded, and the inside of the weld contains five types of defects such as open transverse cracks, unfused bevels, root under-welding, longitudinal cracks and slag entrapment, and the thickness of the plate is 21 mm.

6.1. Simulation

The data generation for the simulation is based on the simulation software CIVA. The transceiver TOFD transducer with a diameter of 6 mm, center frequency of 5 MHz, bandwidth of 70% (-6 dB) is selected, wedge refracting angle is ${60}^{\circ }$, Probe Center Separation(PCS) is 48.5 mm, scanning interval is 0.5mm, and transmitting transducer and receiving transducer are symmetrically arranged on both sides of the weld. At each scanning position, $10\mu s$ data is recorded with sampling rate $50MHz$, and as a result, each A-scan line data has a 500 points. The whole scanning step is 80 steps. The collected simulation data is a 80 × 500 matrix, and forms the D-scan TOFD image shown in Figure 4.

Due to the diffraction effect of ultrasound waves, the tips of defects presents arc-shaped artifacts, which can easily cause errors in the judgment of defect size and burial depth. For this reason, the above proposed frequency-domain synthetic aperture and sparse deconvolution methods are fused to realize the acoustic reconstruction of the D-scan TOFD imaging, results shown as Figure 5. The figure shows that the spatial resolution of ultrasonic TOFD imaging has been greatly improved, which lays a good foundation for the quantitative detection of annular weld defects in terms of size and embedment depth, but it should be noted that the weaker transverse crack signals are easily rejected after the above processing.

6.2. Experimental Studies

The specimen used in the experiment is consistent with those used in the simulation, with a V-type groove welded by argon-arc welding and manual electric-arc welding modes, and 5 sets of artificial defects are embedded same as the simulation. TOFD transducers are same to used in the simulation too, PCS is 48.5mm. These two transducers are tightly pressed against the steel specimen through a mechanical spring-compression device, and water is used to couple between the wedge and the steel plate. The scanning is controlled by a robotic arm, with a step accuracy of 0.5mm. The original imaging result of the scanning is shown in Figure 6.

The time-domain and frequency- domain synthetic aperture focusing imaging are shown in Figure 7 and Figure 8 respectively, which shows that in the unwelded portion of the root, the frequency domain imaging method have a higher imaging contrast and can separate the defects from the bottom echo. At the same time, the combination of sparse deconvolution and frequency domain synthesized aperture can further improve the longitudinal resolution, and as shown in Figure 9, except for the transverse defects are difficult to identify, the longitudinal resolution of the other four defects can be greatly improved, so that the depths of defects can be clearly located as show in Figure 10. In Figure 11, the maximum intensity projection of each defects shows the improvements of transverse resolutions. Especially, in the localized zoomed in image, it can be seen that the full width at half maximum (FWHM) of the defects in raw data image, image reconstructed by time-domain SAFT, image reconstructed by frequency-domain SAFT and image processed by prosed method are 7mm, 5mm, 5mm and 4.5mm respectively. The proposed method can reduce transverse resolution about 35% comparing with raw data and better than 28% for traditional method. Moreover, the noise level of image reconstructed by proposed method is nearly decreased to zero. Thanks to the reduction of the computational volume of the frequency domain synthetic aperture technology, for the 300 ×1000-point TOFD signal data, the synthetic aperture reconstruction process takes about 0.1s, even if the frequency domain synthetic aperture and sparse inverse convolution fusion reconstruction takes about 0.96s, the reconstruction time is less than the sweeping data acquisition time, which meets the needs of on-line imaging inspection applications.

7. Conclusions

In this paper, an ultrasonic TOFD imaging enhancement technology for pipeline grith welds is proposed. This enhancement aims to use the reflection coefficients to image the defects. The reflection coefficients could be solved by eliminating the electrical and spatial pulse responses with sparse deconvolution and frequency-domain synthetic aperture focusing, respectively. Benefiting from the Fast-Fourier Transform, the enhancement process in frequency domain extremely reduces the computational time consumptions. The simulations and experiments demonstrated that using reflection coefficient representing the defects in image could significantly improve both the image longitudinal and transverse resolutions, which is advantage to point out the defect positions accurately. As a result, the proposed method could image the defects with high accuracy and efficiency.

This Work is Supported by Science Foundation of Donghai Laboratory (Grant No. DH-2022ZY0011).