Sparse Blind Deconvolution with Nonconvex Optimization for Ultrasonic NDT Application.

In the field of ultrasonic nondestructive testing (NDT), robust and accurate detection of defects is a challenging task because of the attenuation and noising of the ultrasonic wave from the structure. For determining the reflection characteristics representing the position and amplitude of ultrasonic detection signals, sparse blind deconvolution methods have been implemented to separate overlapping echoes when the ultrasonic transducer impulse response is unknown. This letter introduces the ℓ1/ℓ2 ratio regularization function to model the deconvolution as a nonconvex optimization problem. The initialization influences the accuracy of estimation and, for this purpose, the alternating direction method of multipliers (ADMM) combined with blind gain calibration is used to find the initial approximation to the real solution, given multiple observations in a joint sparsity case. The proximal alternating linearized minimization (PALM) algorithm is embedded in the iterate solution, in which the majorize-minimize (MM) approach accelerates convergence. Compared with conventional blind deconvolution algorithms, the proposed methods demonstrate the robustness and capability of separating overlapping echoes in the context of synthetic experiments.

1. Introduction

Ultrasonic nondestructive testing (NDT) is an established detection technology, which can accurately characterize multilayer or composite materials [1]. The primary ultrasonic NDT tocology mostly considered in practical systems is the pulse-echo method, which evaluates the surface and defects based on the estimation of the time of arrival (TOA), time of flight (TOF), time of flight diffraction (TOFD) and the time difference of arrival (TDOA) of received signals [2]. However, extracting the implicit information is difficult with the distortion and attenuation of the captured signals under noisy environments. Hence, the inspection aims to correctly estimate the parameter values, or separate superimposed waveforms, representing the physical properties.

The conventional method to improve the detection resolution is to increase the transducer frequency, leading to a decrease in acoustic penetration [3]. The other accepted practice estimates the waveform parameters by the Gaussian echo model (or other modified models) and prior information from the transducer [4]. These approaches have some limitations, such as reliance on hardware design and finite propagation distortion. Generally, the signal sample of ultrasonic pulse-echo testing is regarded as the convolution result between the reflection function over the propagation path and the system pulse response. Therefore, deconvolution signal processing methods can be implemented to improve temporal resolution, which significantly affects revealing any small defect object masked by superimposed echoes. The existing mainstream methods are divided into the parametric method and optimization theory method.

Parametric methods depend on prior knowledge of waveform and design specific structures to obtain the estimated parameters. With the assumption that the reflection sequence was sparse, Chang [5] applied the Gauss-Newton (GN) algorithm to parameter estimation of the basic waveform and employed the split variable augmented Lagrangian shrinkage algorithm (SALSA) in the reflection sequence solution. Jin [6] determined the reference echo signal via maximum likelihood estimation (MLE) and achieved sparse deconvolution results by the orthogonal matching pursuit (OMP) algorithm. Bossmann [7] deconvolved the superimposed echo signal by removing the additive noise with a particular matching pursuit (MP) algorithm. The variational mode decomposition (VMD) algorithm was applied to deconvolute the signal by reducing the noise level, which increased the resolution of several defects in different positions [8]. Another class of methods is the minimum entropy deconvolution (MED) technique without acquiring prior assumptions, and the main idea is to find the inverse filter to make the output as sparse as possible [9,10]. In the research conducted by Li [11], norm transformation is applied to enhance the sparsity of MED outputs and eliminate the iteration overflow caused by the long inverse filter. Morphological filtering with the sparse MED algorithm was proposed in [12], where the scale of the signal characteristics can be adaptively selected according to morphological filtering results. Under the premise of mastering the relevant information about the instrument, these algorithms were efficient [3]. Moreover, transforming the ultrasonic parameter estimation into sparse deconvolution was proved to be potentially possible.

However, due to distortion of the waveform in the transmission process, the received signal has a high probability that only the morphological hypothesis (e.g., the sparse prior) can satisfy. Hence, many researchers in recent years have directly abstracted ultrasonic pulse-echo deconvolution into nonconvex optimization with sparsity assumptions, where both the transmitted wavelet and the reflectivity function are unknown. Li [13] constructed a reliable two-step computing framework and applied phase lifting to initialize the blur function. Under statistical assumptions of the signal, Qu [14] studied the multichannel sparse blind deconvolution (MCS-BD) problem via the gradient descent (GD) algorithm. When multiple input signals are present, finding sparse vectors in a subspace can be directly converted into blind calibration [15]. In order to improve the stability of initialization, norm projection [16], least-squares (LS) [17], and pruned tree search [18] are proposed. Zhang [19] explored the influence of symmetric structure on the iterative solution and recovered the shift truncation over the kernel sphere. In simulation experiments, methods based on optimization theory have better antinoise performance than parametric methods.

This letter mainly focuses on sparse blind deconvolution for pulse-echo application, which is an essential branch of ultrasonic NDT. The current challenges of solving this problem with optimization theory are as follows:

Although accurate estimation parameters can be obtained through finite iterations, utilizing the optimization framework is time-consuming [20].

Without prior waveform information, the optimization process is nonconvex, which means that a robust initialization algorithm is inevitable.

We construct a two-step algorithm based on the nonconvex optimization inspired by [13] to solve these problems. The nonconvex model for multichannel ultrasonic signal deconvolution using the ratio is proposed, the solution process of which is divided into two parts. First, with multiple observations, the convolution process is transformed into a sparse blind correction problem in the sparsity case, and the initial estimations of system response and sparse sequence are obtained by the alternating direction method of multipliers (ADMM) [21]. The proximal alternating linearized minimization (PALM) algorithm [22] is applied to solve the nonconvex and nonsmooth optimization problems over each individual received sample. The majorize-minimize (MM) approximation is introduced with symmetric positive definite (SPD) matrix compensation [23] for convergence accelerating purposes. The performance of the proposed algorithms is investigated through simulation analysis and numerical experiment, which demonstrates that the superimposed signals can be significantly faster deconvoluted with reasonable initial estimations.

The remainder of this letter is organized as follows. Section 2 models the multiple blind deconvolution problem and introduces the initialization method and iteration algorithm. Section 3 evaluates the proposed method with performance analysis and compares it with classical deconvolution methods. Section 4 concludes and summarizes this paper.

2. Optimization Models and Methods

2.1. Convolution Model of Ultrasonic Inspection

The ultrasonic excitation signal is usually generated by a high-voltage pulse generator, which stimulates a sensor made of piezoelectric material. The generated wave propagates in the measured object and forms a superimposed ultrasonic echo after reflection by defects and surfaces. We integrate into one signal the measurement model from TOFD and the A-scan experiments regarding the specific ultrasonic pulse-echo application, as illustrated in Figure 1.

Figure 1

Diagram for time of flight diffraction (TOFD) testing and example of A-scan signal: (a) TOFD arrangement for defect inspection; (b) observed receive signal sample by convolution; (c) hypothetical reflectivity sequence.

To analyze the inspection and achieve separation in the time domain, one can regularly assume that the received signal convolves results by the generated ultrasonic pulse-echo with the reflectivity sequence. The latter is a small combination of the reflection coefficients. Hence, the convolution model of the multiple measurement process can be represented as where denotes the convolution operation, , is the observation sample at the receiving end, represents the impulse response (or the blurring function), is the th reflected sequence with the prior sparsity and is the additive noise. We utilized the Gaussian echo model to analyze the impulse response and assumed that the relevant parameters were stable over time in a complete measurement process. Thus: with the parameter vector , which describes amplitude, bandwidth factor, center frequency and phase, respectively. In this letter, we are only interested in the waveform of without solving the specific . Hence by ignoring the subscript in Equation (1), we propose the following blind deconvolution model to recover and from : where is the least-squares objective function with the ratio regularization function. is a regularization parameter. contains prior information about the optimization objects, such as amplitude range, and is the continuous convex function in the domain. The function is the normalized version of the function and holds a scale-invariant, which behaves correctly for deconvolution with irregular blurring functions [24] by considering the following nonconvex optimization algorithm:

Obviously, the solution to this problem is not uniquely identifiable, and we can always find a solution from the equivalence class , without noise, where is the shift symmetry matrix [13]. The convergence to the correct solution relies on the initial estimation of Equation (5). When the center frequency parameter is small, or the basic waveform is the centered Gaussian filter, the global optimal solution can be obtained by random initialization or wavelet initialization function [16]. Because of the constraint from the regularization function, the initialization of the reflected sequence is relatively relaxed. For avoiding the computational difficulty and local optimal solution of solving the model, a well-produced initialization is crucial to the success of Equation (5).

2.2. Initialization Based on Blind Gain Calibration

Without considering the noise, we modify Equation (1) as: where is the Toeplitz circulant matrix of , and can be written as

By denoting that multiple measurements can be written in matrix form, i.e., and , we rewrite Equation (6) as:

The circulant matrix has the equivalent eigenvalue decomposition, thus: where is the unitary discrete Fourier transform (DFT) matrix, and . By applying the DFT matrix to both sides of Equation (8), we have:

Let and, the measurement process can be parameterizing in a linear representation: in which is the element of the th row and the th column in , denotes the th row of . The bilinear inverse problem of is transformed into a linear inverse problem of . Even in the case of actual noise interference [25], we also can solve the problem:

The proposed initialization algorithm becomes equivalent to blind gain calibration proposed in [26], which leads to solving:

Evidently, it satisfies for the pair . To avoid this, we introduce an additional linear constraint on , thus: where . The combination of Equations (13) and (14) is a convex problem, which can be solved using the existing solvers.

However, numerical simulations indicate that the direct use of this linear equality constraint complicates the problem. Therefore, we use the ADMM algorithm, which updates local subproblems to coordinate the global problem and solve Equation (13). An equivalent operator of the sum constraint is applied in the iteration process to simplify the initialization process.

First, the augmented Lagrangian function of Equation (13) with the expansion of Equation (11) is expressed as: where is the penalty parameter and is the Lagrange multiplier matrix. The Equation (13) processes in following steps by updating parameters individually:

For update for , can be regarded as constant values, and the minimization problem is reduced to: where is the th column in . The optimization problem for is stated as follows:

Let , and thus: which means an -Lipschitzian gradient for every (see Appendix A). The Lipschitzian constants of is . Hence, the update rule of is defined based on an approximation: where:

Because is the norm function, the iterative shrinkage-thresholding (IST) method [27] can be used to find a closed-form solution: where means the th element of (resp. ). is the soft-thresholding function, which is defined as:

Mote that is separable and consists of , the update of in Equation (13) can be rewritten as:

For updating , is regarded as constant values, and thus: and when , is represented as: where , and:

Similar to , we use directly to update . Considering the constraints from formula Equation (14), we modify the update method of as: where satisfies the sum constraint.

Thus, through Equations (16), (24), and (28), the impulse response and sparse sequences from each measurement are estimated up to specific scaling. Despite the arbitrary selection of and the assumption without noise, the initialized pair is close enough to the ground truth, and guarantees Equation (5) can converge to the global optimal solution approximately.

2.3. Alternating Optimization Method Based on PALM

The basic PALM method is proposed for minimization of a sum of finite functions and requires the smooth function (or partially Lipschitz properties). However, the ratio is both nonconvex and nonsmooth, which means the smoothed ratio proposed in [16] can replace the ratio in . Hence, we rewrite Equation (4) as , and: which can be smooth approximations for sparse representation. Besides, we assume that is separable and is a proper, lower semicontinuous, convex function. Thus: where and are the indicator functions with prior knowledge. Equation (5) is represented as:

By starting with the initial estimates obtained from the initialization, we generate the sequence via the scheme:

The PALM for solution of Equation (32) are defined respectively by: where the proximal operator is defined by the proximal map, which is associated to a specific function, i.e., and are the partial Lipschitz constants of (see Appendix A).

The solution of the problem is converted to determine the appropriate and , so that the function Equation (33) converges. From idea of gradient descent, and can be selected as: where and are the Lipschitz constants for and [22]. Similar to the algorithm in [28], can be considered constants when updating , and: where is the Toeplitz matrix operator. For , we rewrite the update step as: and is the symmetric positive definite (SPD) matrix which satisfies the MM principle [23], which can be obtained by building majorizing approximations for [29]. Thus: where is the Lipschitzian constants of . By Equation (8), we have: and . The other items in are given by the proposition established in [16].

Thus far, for each measurement sample , we can obtain the corresponding estimated by solving Equation (32) with a global initialization from Equation (16).

3. Simulation Results

3.1. Stable Initialization with Phase Transitions

Phase diagrams, which demonstrate the empirical probability of success over a range of sparsity and samples for a fixed sampling window length, are implemented to evaluate the feasibility of the proposed initialization algorithm. We classify a success initialization of if . Figure 2 presents variations of the initialization success rate by phase transitions, concerning the -sparsity of each and the number of samples . The additive noise is considered as white Gaussian noise with the variance . The length of the sequence and the recovery success rate are calculated using 50 Monte Carlo simulations on every grid point. We rescale the estimated to have the same norm as the correct .

Figure 2

Phase transitions of the proposed initialization algorithm with different noise levels: (a) noise free; (b) ; (c) ; (d) . The gray bar denotes the success rate.

As shown in Figure 2a, the initialization algorithm can estimate correctly when the sparsity and the number of samples are appropriate. To achieve successful initialization with Equation (16), must scale linearly with . However, when does not satisfy the sparse hypothesis, i.e., does not hold, blind gain calibration cannot achieve convergence to the neighborhood of . In the case of noise interference, the proposed method has a high probability of obtaining with increased samples. Because the influence of noise is not introduced in Equation (11), few samples lead to additional errors for the iterative recovery. Meanwhile, as shown in Figure 2d, the probability of successful initialization decreases with serious noise interference.

3.2. Numerical Comparison for Deconvolution Evaluation

Figure 3 illustrates the comparison of conventional parametric deconvolution methods numerically with the same noise variance (). The regularization parameters and indicator functions of [12,29] are adjusted in this comparison. In the case of , , and , offset and polarity reversal presents in the conventional fast deconvolution methods, such as MED [11] and LS [17]. The MED algorithm uses an inverse filter to solve sparse reflection sequence directly. However, the uncorrelation between the inverse filter and the blurring function results in deviation and fluctuation, as shown in Figure 3a. Besides, the MED algorithm extracts the arrival time by manually modulating the threshold, which reduces detection stability. The LS algorithm estimates the blurring function and analyzes the concentration of the error, where the relatively accurate time parameters can be obtained. The regularized coefficient of noise is based on the prior information of the signal. Therefore, the deviation appears in the estimation of amplitude parameters. The reflection sequence, estimated by LS (Figure 2b), is sparser than that estimated by MED. Although MED uses sparse constraint, the noise also affects the deconvolution results.

Figure 3

Comparison of parametric methods: (a) norm MED (the threshold is used to determine detection points); (b) LS algorithm.

Compared with parametric methods, the method based on optimization theory has better adaptability to the superimposed echo. Therefore, to evaluate the performance of the proposed algorithm, the sparsity is increased appropriately. Further, when , , and (), Figure 4 shows the initialization comparison of the Smoothed One-Over-Two (SOOT) [29] algorithm and the proposed method. The centered Gaussian filter from SOOT is different from the basic ultrasonic waveform, where converging on local optimum is inevitable. Although the blind gain calibration is implemented with time cost, the proposed algorithm provides a starting point close to the original blurring function for alternating convergence.

Figure 4

Comparison of initialization algorithms.

Figure 5 reveals the results of the initialization algorithm and alternating solution. Although the SOOT algorithm reduces the influence of initialization by variable metric strategies and increasing the number of subiterations, the spike interference still affects the solution where the initialization of the centered Gaussian filter is not applicable to the true ultrasonic response. Despite limited missing solutions and amplitude deviation, the proposed algorithm can accurately estimate the position and amplitude of the reflected sequence in the case of relative sparsity.

Figure 5

Deconvolution results of different methods with optimization theory: (a) SOOT; (b) proposed algorithm.

For a fair comparison, we used the same noise conditions and parameter settings as illustrated in [29], where the fixed noise standard variances (0.01, 0.02, and 0.03) and regularization parameters are considered. Table 1 compares the residual deconvolution error, where and represent the norm distances between the estimated values and the correct values for different noise levels, respectively. The running time is obtained by taking the average value from experiments, and the stopping criterion for iterative algorithms is set as . Although the deconvolution based on LS [17] is superior in speed, the cost is a large estimation error. MED-based methods [9,10,11] achieve relatively good sparse sequence estimation in less computing time. However, because of the characteristics of the inverse filter, the basic response of the signal is inaccessible. A solution using alternate optimization can obtain a better solution at the cost of computing time. The CVX toolbox without improvement causes additional time consumption due to the introduction of random initialization. Besides, because the optimization conditions do not make specific constraints on the ultrasonic waveform, the estimation error of the blurring function is larger than that of the reflection sequence. Meanwhile, the proposed algorithm is significantly faster than the other algorithms based on nonconvex optimization. The proposed initialization is closer to the ground truth, where the algorithm can achieve the stopping criterion more quickly.

Table 1

Comparison of algorithm deconvolution results. The interpretation platform used the Intel Core i5-9400f processor with 16GB memory.

Algorithm	σ=0.01		σ=0.02		σ=0.03		Running Time (s)
	ℓ_2,h	ℓ_1,x	ℓ_2,h	ℓ_1,x	ℓ_2,h	ℓ_1,x
LS	2.3981	1.1092	2.5495	1.5212	3.4781	2.5008	3.89
ℓ_0 MED	\	0.8121	\	1.1002	\	1.3063	5.62
M-S-MED	\	0.7203	\	0.8212	\	1.0304	8.20
OMED	\	0.5548	\	0.6167	\	0.9033	10.34
ℓ_1/ℓ_2 CVX	1.5301	0.1023	1.6623	0.4180	2.0301	0.8902	65.29
SOOT	0.8923	0.0221	1.2340	0.0936	1.6721	0.1645	34.50
proposed	0.3622	0.0133	0.4528	0.0584	0.8973	0.1045	22.32

4. Conclusions

In this letter, a deconvolution method based on nonconvex optimization is proposed. By simplifying the problem to blind gain calibration, a two-step iterative initialization can effectively obtain the initial guesses for ultrasonic pulse-echo response and reflection sequence with the ADMM algorithm. Then, using ratio regularization, an exact algorithm based on the PALM framework is applied. Simulation and numerical analysis illustrate that the proposed method algorithm has the probability of successful initialization. In the algorithm comparison experiment, deconvolution results are obtained with high accuracy and acceptable time cost.

The proposed method has potential gains in ultrasonic NDT applications for defect detection and sensor parameter estimation. When the specific instrument parameters are ambiguous or the propagation is distorted, the parameters representing the measured object can be obtained accurately through sparse blind deconvolution. However, the time cost limits the application of the algorithm in other areas (such as ultrasonic imaging) unless the processing speed is further improved.