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.


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 0 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, 2 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 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.

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.
Sensors 2020, 20, x FOR PEER REVIEW 3 of 14 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.

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. 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, 1 ≤ ≤ , y ∈ ℝ is the observation sample at the receiving end, h ∈ ℝ represents the impulse response (or the blurring function), ̅ ∈ ℝ is the th reflected sequence with the prior sparsity and w 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: 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, 1 ≤ l ≤ L, y l ∈ R n is the observation sample at the receiving end, h l ∈ R s represents the impulse response (or the blurring function), x l ∈ R n is the ith reflected sequence with the prior sparsity and w l 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 h(θ, t) without solving the specific θ. Hence by ignoring the subscript l in Equation (1), we propose the following blind deconvolution model to recover h and x from y: where is the least-squares objective function f 0 (x, h) = 1 2 h * x − y 2 2 with 1 / 2 the ratio regularization function. λ is a regularization parameter. g(x, h) contains prior information about the optimization objects, such as amplitude range, and is the continuous convex function in the domain. The 1 / 2 function is the normalized version of the 1 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 y l = R + l h * −1 R − l x l , 0 without noise, where R is the shift symmetry matrix [13]. The convergence to the correct solution relies on the initial estimation (x 0 , h 0 ) 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 h 0 is crucial to the success of Equation (5).

of 14
The circulant matrix C(g) has the equivalent eigenvalue decomposition, thus: where F is the n × n unitary discrete Fourier transform (DFT) matrix, andg = √ nFg ∈ C n . By applying the DFT matrix F to both sides of Equation (8), we have: LetY = FY and, the measurement process can be parameterizing in a linear representation: in whichy i,l is the element of the ith row and the lth column inY, F H i denotes the ith row of F. The bilinear inverse problem of (g, X) is transformed into a linear inverse problem of (τ, X). 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 (0, 0). To avoid this, we introduce an additional linear constraint on τ, thus: where 1 T = [1, 1, . . . , 1]. 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 ρ > 0 is the penalty parameter and λ is the Lagrange multiplier matrix. The Equation (13) processes in following steps by updating parameters individually: Sensors 2020, 20, 6946 6 of 14 For update for X, τ (k) can be regarded as constant values, and the minimization problem is reduced to: whereb (k) l is the lth column in diag τ (k+1) Y . The optimization problem for x l is stated as follows: 2 , and thus: which means an L-Lipschitzian gradient for every x ∈ R n (see Appendix A). The Lipschitzian constants L 0 of ζ 0 (x) is ρ F 2 2 . Hence, the update rule ofx l is defined based on an approximation: where: Because x l 1 is the 1 norm function, the iterative shrinkage-thresholding (IST) method [27] can be used to find a closed-form solution: where x (k+1) l|j means the jth element of x (k) l (resp. z (k) j ). soft(·) is the soft-thresholding function, which is defined as: Mote that X (k+1) is separable and consists of x (k+1) l| j , the update of X (k+1) in Equation (13) can be rewritten as: For updating τ (k) , X is regarded as constant values, and thus: and when ∇L ρ X (k+1) , τ i , λ (k) = 0, τ (k+1) i is represented as: where U, V ∈ C n , and: Similar to X (k+1) , we use U (k+1) /V directly to update τ (k+1) . Considering the constraints from formula Equation (14), we modify the update method of X (k+1) as: where trace diag τ (k+1) = c 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 c and the assumption without noise, the initialized pair (x 0 , h 0 ) is close enough to the ground truth, and guarantees Equation (5) can converge to the global optimal solution approximately.

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 1 / 2 ratio is both nonconvex and nonsmooth, which means the smoothed 1 / 2 ratio proposed in [16] can replace the 1 / 2 ratio in f (x, h). Hence, we rewrite Equation (4) which can be smooth approximations for sparse representation. Besides, we assume that g(x, h) is separable and is a proper, lower semicontinuous, convex function. Thus: where g 1 and g 2 are the indicator functions with prior knowledge. Equation (5) is represented as: By starting with the initial estimates x (0) , h (0) 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., prox σ,t (z) = argmin c (k) and d (k) are the partial Lipschitz constants of ∇ f (see Appendix A). The solution of the problem is converted to determine the appropriate c (k) and d (k) , so that the function Equation (33) converges. From idea of gradient descent, c (k) and d (k) can be selected as: where L 1 and L 2 are the Lipschitz constants for ∇ x f and ∇ h f [22]. Similar to the algorithm in [28], x (k+1) can be considered constants when updating h (k+1) , and: where T x : H n → H n×n is the Toeplitz matrix operator. For x (k+1) , we rewrite the update step as: and A x (k) , h (k) is the symmetric positive definite (SPD) matrix which satisfies the MM principle [23], which can be obtained by building majorizing approximations for f (x, h) [29]. Thus: where L 1 is the Lipschitzian constants of ∇ x f 0 x, h (k) . By Equation (8), we have: and L 1 = M 2 2 . The other items in A x (k) , h (k) are given by the proposition established in [16]. Thus far, for each measurement sample y i , we can obtain the corresponding estimated x i ,ĥ i by solving Equation (32) with a global initialization (X o , h 0 ) from Equation (16).

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 Figure 2 presents variations of the initialization success rate by phase transitions, concerning the K-sparsity of each x l and the number of samples L. The additive noise w l is considered as white Gaussian noise with the variance σ 2 . The length of the sequence n = 100 and the recovery success rate are calculated using 50 Monte Carlo simulations on every grid point. We rescale the estimated h 0 to have the same norm as the correct h.
proposed initialization algorithm. We classify a success initialization ℎ of ℎ if ℎ − ℎ / ℎ ≤ 0.3. 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 = 100 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 ℎ. 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 As shown in Figure 2a, the initialization algorithm can estimate h correctly when the sparsity and the number of samples are appropriate. To achieve successful initialization with Equation (16), L must scale linearly with K. However, when K does not satisfy the sparse hypothesis, i.e., K n does not hold, blind gain calibration cannot achieve convergence to the neighborhood of h. In the case of noise interference, the proposed method has a high probability of obtaining h 0 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. Figure 3 illustrates the comparison of conventional parametric deconvolution methods numerically with the same noise variance (σ = 0.02). The regularization parameters and indicator functions of [12,29] are adjusted in this comparison. In the case of K = 3, L = 20, and n = 1000, 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.

Numerical Comparison for Deconvolution Evaluation
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 0 constraint, the noise also affects the deconvolution results.
functions of [12,29] are adjusted in this comparison. In the case of = 3, = 20, and = 1000, 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. 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 = 20, = 20, and = 1000 ( = 0.02), 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. 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 K = 20, L = 20, and n = 1000 (σ = 0.02), 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.       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 ( ) − ( ) ≤ 10 . 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 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 2,h and 1,x 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 x (k+1) − x (k) 1 ≤ 10 −3 . 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.

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 1 / 2 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.