Approximation performance of model functions for blind deconvolution of ultrasonic reflections

1. Introduction

Ultrasonic inspection is offering a useful tool for material mechanical integrity testing [1-3]. Pulse-echo mode of non-destructive testing is considered in this research. Acoustic wave traveling in the test material encounters the reflection from inhomogeneities (defects) along the travel path [4]. Such echoes produced carry the information about the location of the defect as the time of flight (ToF) and the severity of the defect or even mechanical properties as the reflection amplitude. ToF can be converted to distance if propagation velocity is known and the reflection amplitude can be converted to mechanical impedance which in turn is related to density and velocity (stiffness) [3, 4]. Estimation of the aforementioned parameters for the single echo does not present significant problem, but when several reflectors are present along the wave path the reflections pile up obscuring the reliable estimation of the echoes parameters [5]. If probing signal is of indefinite bandwidth such problem does not exist, but bandwidth is always limited in practice [6, 7]. Furthermore, limited bandwidth is usually desirable, because signal acquisition system requirements can be relaxed [8], attenuation of ultrasound at high frequencies produces structural noise and signal degradation [9]. Deconvolution is used to address this problem [10-12]. It is based on assumption that reflections are sparse (based on usual defects structure which present themselves as layers or distributed cavities) and that every reflection is a time and amplitude translated copy of the original probing signal (reference). Inverse filtering [13, 14] is dividing the received signal by reference in frequency domain in order to produce the impulse response of the media being probed. Unfortunately, very high signal-to-noise ratio (SNR) is needed in order to obtain acceptable results. Wiener filtering [13] can be used to address the problem, but it does not offer the decomposition of the signal into ToF and corresponding amplitudes. Iterative deconvolution techniques allow for efficient signal decomposition either processing one reflection at a time or several reflections simultaneously [6, 10, 12, 15]. Yet, production of the reference signal constitutes a major problem in all deconvolution methods [5, 6, 16, 17]. Signal spectral content is altered along the propagation path or at reflection interface [4]. Together with noise present, this significantly degrades the results of deconvolution. Blind deconvolution [17, 18] is offering a production of reference signal using some simple mathematical model function, Gabor (harmonic oscillation with Gaussian envelope) being the most frequently used [6, 19]. Result of such deconvolution significantly improves the SNR, adapts the signal spectral content changes and significantly reduces the required amount of computer memory needed to store or transmit the acquired data [20] by compressive sensing. Analytical processing is required for the estimation of the fitting function parameters [6, 10-12, 17, 19]. Unfortunately, Gabor function goes against causality principle of ultrasonic inspection and is only suitable for narrowband signals. Attempts to split the envelope slope into two [16, 17] are not favored since analytical processing of such function is more complicated. The goal of the presented investigation was to compare the approximation performance of the model functions. New model function was proposed where the right-tailed asymmetry (causality) of the signal is easily modelled without the need for the piecewise-linear equations. Comparison of three model functions was done using four experimentally collected ultrasonic signals.

2. The ultrasonic signals

As it was mentioned above, narrowband signal is easy to approximate by Gabor function [6, 19], yet the wideband case is the more complicated. Therefore, investigation was concentrated on wideband signals. Several wideband ultrasonic transducers were used to collect the representative signals: i) spherically focused 5 MHz center frequency and 64 % bandwidth transducer IRY405 from NDT transducers PLL (further labelled 5F); ii) spherically focused 10 MHz center frequency and 71 % bandwidth transducer IRY210 from NDT transducers PLL (labelled 10F); iii) composite, unfocused 5 MHz center frequency and 86 % bandwidth transducer TF5C6N-E from Doppler Electronic Technologies (further – 5U) and iv) composite, unfocused 10 MHz center frequency and 66 % bandwidth transducer TF10C6L from Doppler Electronic Technologies (further – 10U). The polished 30 mm height and 20 mm diameter stainless steel cylinder was immersed in water tank (refer Fig. 1) and reflections collected using ultrasonic data acquisition system [8].

Fig. 1Experiment setup used for ultrasonic signals collection

Fig. 2A-scans of the reflections for a) 5 MHz transducers and b) 10 MHz transducers

a)

b)

Transducer excitation was done using a unipolar pulser, using 10 V rectangular pulse of 50 ns duration for 10 MHz transducers and 100 ns duration for 5 MHz transducers. Receiver was 0.5-30 MHz bandwidth with 30 dB gain. One thousand A-scans with steel slab reflection present were collected and averaged to reduce the noise. A-scans of such reflections for the corresponding transducers are presented in Fig. 2.

Spectrums of the corresponding reflections have been calculated using the discrete Fourier transform (DFT). Refer Fig. 3 for corresponding reflections spectrums and indications of –3 dB, –6 dB and –20 dB bandwidth. These bandwidth values were later used for selecting the passband parameters in additional signal filtering.

Fig. 3Reflection spectrum for a) 5 MHz transducers and b) 10 MHz transducers

a)

b)

Unfocused type transducers are composite type, and have significant bandwidth therefore their spectral temporal and spectral response is more complicated. Focused transducers are of single crystal type so their response is closer to simple resonant tank. Signals were stored for further processing.

3. Candidate functions

Ultrasonic reflection signal can be approximated as Gaussian modulated cosine or Gabor function [19]:

where $A$ is amplitude, $\tau$ is the width of the envelope of the echo, $t_0$ is the time of flight (ToF), $f_0$ is center frequency and $\phi$ is phase. These function parameters are closely related to the properties of ultrasonic signal [20]. The amplitude $A$ can be related to the size of the reflector, $f_0$ and $\tau$ are related to the used transducer, ToF represents the temporal (as well as spatial) position of the reflector, $\phi$ accounts for the distance, impedance, size, and orientation of the reflector. These parameters are expressed as vector ${\mathbf{\Gamma }}_G$:

Example of Gabor function is presented in Fig. 4. Envelope was obtained using Hilbert transform. It can be seen that such signal exhibit close similarity to Fig. 2 yet are lacking the asymmetry and causality. The decomposition of the signal using the Gaussian function can be compromised if received echoes are more complex [20]. New, Gaussian Chirplet model with additional chirp sinusoidal component was proposed in [6, 17]:

where $\psi$ is the linear chirp rate. The parameters are expressed as vector ${\mathbf{\Gamma }}_{GC}$:

Example of Gaussian Chirplet model is presented in Fig. 4.

Fig. 4a) Example of Gabor and b) Gaussian Chirplet functions

a)

b)

Unfortunately, despite this being better than simple Gaussian, the envelope of this model is still symmetric with respect to its peak position. In this study we introduce a new echo model that has the asymmetric envelope component with chirp as carrier:

where $\alpha$ is asymmetry coefficient. Envelope function was derived from lognormal distribution function with the modification of the tipping point so it corresponds to ToF position $t_0$. Corresponding parameters are presented as vector ${\mathbf{\Gamma }}_{LC}$:

Refer to Fig. 5 for representative examples of the new proposed function.

Fig. 5Example of new function with a) linear frequency modulation and b) without

a)

b)

Essential, that new function obeys the causality principle: analyzing Eq. (5) one can see that time axis and so the signal can exist only from certain temporal position in front of $t_0$. Same can be deduced from Fig. 5 analysis: there is no signal energy beyond the certain point in front of the signal. Furthermore, asymmetry of the signal is produced on uniform time axis, without the need for piecewise approximation as per [16, 17].

4. Estimation of the parameters of the approximating functions

It is important to estimate parameters of a model with least computational effort and accurately. In [19] Demirli and Saniie used Gauss-Newton algorithm to obtain the ultrasonic echo parameters. However, the drawback of this algorithm is that the results of Gauss-Newton algorithm can be local optimum rather than global optimum. The reason for this is that Gauss-Newton algorithm depends on initial solution. An Artificial Bee Colony (ABC) algorithm was introduced in 2015 to overcome those drawbacks [21, 22].

Deconvolution can be divided into two minimization problems: i) “pure” deconvolution where signal is treated as convolution of test object reflectivity function with probing (reference) signal [10-12] and ii) reference signal production [18, 19, 23]. Deconvolution is the optimization problem that uses sparse deconvolution methods such as matching pursuit, Prony model or orthogonal matching pursuit in order to deconvolve the mixed data [5, 6, 15, 16]. In general, reference signal can be treated as time-invariant here and can be pre-recorded from some interface exhibiting strictly step response. Unfortunately, sometimes the reference signal is not available or is time-variant. In such case second optimization problem arises which was named blind deconvolution [17-20, 22, 23]. Here signal echo is simulated as simple mathematical function whose model parameters can be iteratively adjusted. In current research we only concentrate on reference signal derivation, adjusting the model parameters. Adjustment can be carried out using iterative Newton method. We suggest different approach with following improvements: initially $t_0$ (ToF), $A$ (amplitude) and $\phi$ (phase) are directly derived using the correlation processing:

where $\tilde{x}$ is the complex correlation function:

where $s$ is the original signal, ${\widehat{s}}_{\mathrm{c}\mathrm{o}\mathrm{s}LNCh}$ is the model function Eq. (5) without $\phi$ and $\psi$ components and ${\widehat{s}}_{\mathrm{s}\mathrm{i}\mathrm{n}LNCh}$ is the model function Eq. (5) without $\phi$ and $\psi$ components and sin instead cos. Both model and the original signals are discrete versions. ToF ($t_0$) is interpolated using cosine technique described in [24] and $A$ and $\phi$ are using linear interpolation.

Parameters $\alpha$, $f_0$, $\tau$ and $\psi$ are adjusted using Nelder-Mead direct search (simplex search [25]) which does not require derivatives but only the objective function e.g. remainder of subtraction:

Such approach slightly differs from commonly used ${\mathcal{l}}_2$ norm: instead of analyzing just the remainder ${‖s\left(t\right)-{\widehat{s}}_{\mathrm{\Gamma }}\left(t\right)‖}_2$ we analyze how much of the reduction was obtained relative the ${\mathcal{l}}_2$ norm of the original signal.

Since both numerator and denominator at Eq. (14) can be attributed to the energy, this parameter was termed remainder energy reduction. The goal of the deconvolution is the complete removal of the reflection from the signal so it is not interfering with remaining reflections. The stronger is the remainder reduction the better is performance of the succeeding deconvolution steps. It might be presumed that remainder energy (or ${\mathcal{l}}_2$ norm) should be a sufficient metric for approximation (so the deconvolution) performance evaluation, but reflection amplitude can be small, so even slight reduction can give acceptable results (refer Fig. 6 for the example of the subtraction remainder representation in time domain).

Fig. 6a) Original A-scan (thin solid line) of the 5 MHz focused transducer signal and its approximation by Gabor (solid) and new function (dash-dot) and b) the corresponding remainders in time domain

a)

b)

Such metric does not represent the actual approximation performance. Ideal approximation should produce very small remainder assuming signal is free from noise (which was accomplished by significant averaging). Low quality of approximation will not be able to offer full stripping of the reflection and the remainder will be high. Then the remainder energy reduction should be the metric to be used and it was further used as approximation performance evaluation criterion. This optimization criterion also satisfies a general deconvolution goal: a reference signal that produces minimal remainder after subtraction.

Initial guess of the $f_0$, was produced using center of gravity of signal spectrum $S\left(f\right)$:

Initial guess of the parameter $\tau$ was produced using half value full width criteria of Hilbert envelope of signal temporal representation:

Initial guess of the parameter $\alpha$ was produced using:

Initial guess of the parameter $\psi$ was produced using:

where $\beta$ is the envelope bandwidth:

However, a new problem arises with these algorithms. Depending on initial conditions, the algorithm may remain on the local minimum.

Iterative Gauss Newton optimization method was applied [6] for its high convergence speed. Minimizing the sum of squares of the residual is a common problem. Taking advantage of the echo model, the Gauss Newton iteration formula for estimation of the parameter vector with an initial guess ${\mathrm{\Gamma }}^{\left(k-1\right)}$ can be written as:

where $D_{\mathrm{\Gamma }}{\widehat{s}}_{{\mathrm{\Gamma }}^{\left(k-1\right)}}\left(t\right)$ denotes the Jacobian at initial solution ${\mathrm{\Gamma }}^{\left(k-1\right)}$, i.e.:

Hence, the update vector $d\mathrm{\Gamma }={\left(d\phi ,d\psi ,d\alpha \right)}^T$ is obtained from the Eq. (27):

This problem can be solved directly. Then we obtain a new update:

and after fixed $r$ of Newton steps, the new estimate is obtained:

Performance investigation of the aforementioned algorithm is presented in next section.

5. Experimental results

Quality of the convergence was inspected, by artificially varying the model parameters around the final optimization result in order to confirm that global minimum was reached. Refer Fig. 7 for convergence analysis of the parameters $\alpha$, $\phi$ and $\psi$ for focused 5 MHz transducer.

Fig. 72D representation of the remainder energy reduction versus a) φ and ψ and b) α and ψ parameters with simplex (star) and Gauss Newton (circle) optimization for focused 5 MHz transducer

a)

b)

Fig. 8Approximation performance of candidate functions for ultrasonic reflections collected using different transducers when original signal was filtered (horizontal lines – unfiltered case)

a)

b)

c)

d)

Initial guess was out of presented range. Nelder-Mead simplex search produced near optimal result (star) and Gauss Newton optimization (circle) provided fine adjustment. It can be concluded that proposed optimization algorithm locates the global minimum and can be used for the comparison of the model functions.

Further investigation was aimed on deciding whether additional signal filtering is required. Approximation of all four representative signal types was analyzed for the remainder energy reduction in case of raw and filtered signals. Signal cut-off frequencies were chosen basing the –3 dB, –6 dB and –20 dB points established in Section 2. Zero-phase filtering by fourth order Butterworth IIR filter was used. Refer to Fig. 8 for the approximation performance comparison.

It could have been expected that approximation performance will not be high due to the complex nature of the reflections (refer Fig. 9 for remainders comparison of the two most differing reflections, note the new scale compared to Fig. 6). Filtering reduces the complexity of the signal and further improves the SNR. Filters used do not alter the main passband of the signal (–3 dB worst case), albeit improve the approximation performance.

Fig. 9Remainder comparison for a) complex and b) more simple reflections

a)

b)

It’s essential to point out that zero-phase filtering has to be used because otherwise signal spectral components are spread in both directions (due to high and low pass filer transition band). It can also be noted that new asymmetrical model function offers the best performance compared to Gabor and Gaussian Chirplet functions.

6. Conclusions

It was demonstrated that simple model functions like Gabor or Gaussian Chirplet demand for envelope asymmetry. New model function, obeying the causality of the ultrasonic signals was proposed. Function possess smooth envelope so can be easily differentiated.

Approximation parameters optimization technique was developed. Technique is based on combination of Nelder-Mead simplex search and iterative Gauss Newton optimization method. Remainder energy reduction of the model function subtraction from the original signal was used as optimization criteria. Grid plots of the remainder reduction versus optimization parameters confirm that technique locates the global maximum.

Four wideband ultrasonic signals were collected as candidates for approximation. Signals were deliberately chosen to significantly deviate from simple model function. Approximation performance investigation was carried out on Gabor, Gaussian Chirplet and newly proposed model functions. It was expected that approximation quality will be low, albeit after some filtering of the signal remainder energy reduction down to 29 dB was achieved.