Nans Laroche

Sébastien Bourguignon

Jérôme Idier

Ewen Carcreff

Aroune Duclos

Fast non-stationary deconvolution of ultrasonic beamformed images for nondestructive testing

This paper addresses high-resolution ultrasonic image reconstruction from Full Matrix Capture (FMC) data in the context of nondestructive testing (NDT). In order to reduce the numerical complexity, the time-domain data and ultrasonic model are projected into the image domain through a linear beamforming procedure. The resulting model is interpreted as a shift-variant convolution process, affected by non-stationary and colored noise. An interpolation procedure is built in order to account for the spatial variations of the resulting point spread function. Under the same methodological framework, an approximate whitening filter is proposed and incorporated in the forward model. Both constructions then allow fast computations and limited memory storage. Deconvolution is performed by minimizing the least-squares data misfit error, with a penalization term favoring sparsity and spatial continuity of the output images. Results with synthetic data show that the proposed approach gives performances close to the inversion of raw FMC data, while being computationally much more efficient. The method is finally applied to laboratory data for the inspection of a stainless steel block containing closely spaced and small sidedrilled holes calibrated flaws. Successful detection and separation is achieved for flaws with diameters six times smaller than the wavelength, and distant from each other by four times less than the resolution limit given by the Rayleigh criterion.

I. INTRODUCTION

A. State-of-the-art in NDT ultrasonic imaging M ODERN ultrasonic array probes contain a large number of transducers, offering high performance in terms of flexibility, speed and quality of inspection. In the field of nondestructive testing (NDT), the main advantage of the arraybased ultrasonic modality is its ability to produce images from a single acquisition sequence through beamforming techniques [START_REF] Drinkwater | Ultrasonic arrays for non-destructive evaluation: A review[END_REF]. Beamforming has been first hardware-implemented by applying different delay laws to the transducers in order to focus the beam at different locations of the inspected region. Nowadays, software-implemented beamforming is generally preferred, being obviously more flexible. In particular, the recent growth of Graphic Processing Units (GPU) capabilities has significantly improved the speed of beamforming algorithms, which can now compute hundreds of images per second [START_REF] Sutcliffe | Realtime full matrix capture for ultrasonic non-destructive testing with acceleration of post-processing through graphic hardware[END_REF].

Until recently, software-implemented beamforming was mostly performed using the Synthetic Aperture Focusing Technique (SAFT) [START_REF] Karaman | Synthetic aperture imaging for small scale systems[END_REF], which consists in recording the pulse-echo responses of a moving single-element transducer. However, N. Laroche and E. Carcreff are with the Phased Array Company (TPAC), Nantes, France (e-mail: nans.laroche@tpac-ndt.com).

N. Laroche, S. Bourguignon, and J. Idier are with the Laboratoire des Sciences du Numérique de Nantes (LS2N), Nantes, France.

A. Duclos is with the Laboratoire d'Acoustique de l'Université du Mans (LAUM), Le Mans, France.

SAFT is becoming obsolete with the emergence of array probes, which enable more efficient data acquisition schemes. In medical ultrasound imaging, Plane Wave Imaging (PWI) [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF] is now the standard method, where all elements emit plane wave fronts at different angles and receive synchronously. This technique maximizes the transmitted energy to the tissue at a high frame rate (up to several thousand images per second), which is critical in order to capture motion information [START_REF] Tanter | Ultrafast imaging in biomedical ultrasound[END_REF].

In NDT, which is the target application area of this paper, obtaining high-resolution information is necessary in order to detect and characterize accurately possible flaws in the inspected piece [START_REF] Wilcox | Advanced Reflector Characterization with Ultrasonic Phased Arrays in NDE Applications[END_REF], [START_REF] Felice | Sizing of flaws using ultrasonic bulk wave testing: A review[END_REF]. Typical examples of interest are the control of additively manufactured parts [START_REF] Chabot | Towards defect monitoring for metallic additive manufacturing components using phased array ultrasonic testing[END_REF], [START_REF] Obaton | Efficient volumetric non-destructive testing methods for additively manufactured parts[END_REF], welds [START_REF] Zhang | Defect detection using ultrasonic arrays: The multi-mode total method[END_REF] or power generation units [START_REF] Busse | Review and discussion of the development of synthetic aperture focusing technique for ultrasonic testing (SAFT-UT)[END_REF]. In such contexts, it is a priority to collect more exhaustive information than in medical imaging, at the price of slower acquisition rates. Therefore, a beamforming method of choice is the Total Focusing Method (TFM) [START_REF] Felice | Sizing of flaws using ultrasonic bulk wave testing: A review[END_REF], [START_REF] Holmes | Post-processing of the full matrix of ultrasonic transmit-receive array data for non-destructive evaluation[END_REF] applied to Full Matrix Capture (FMC), that is, the collection of all inter-element responses from an array probe. FMC has become a reference acquisition method for various typical applications in NDT. It is highly sensitive to small defects and is used for sub-wavelength reflectors detection [START_REF] Wilcox | Advanced Reflector Characterization with Ultrasonic Phased Arrays in NDE Applications[END_REF]. It provides valuable information for the sizing and characterization of crack-like defects using diffraction and specular echoes [START_REF] Zhang | Defect detection using ultrasonic arrays: The multi-mode total method[END_REF], [START_REF] Felice | Accurate depth measurement of small surface-breaking cracks using an ultrasonic array postprocessing technique[END_REF]- [START_REF] Peng | The sizing of small surface-breaking fatigue cracks using ultrasonic arrays[END_REF]. It is also efficiently used for the separation of close scatterers [START_REF] Simonetti | Multiple scattering: The key to unravel the subwavelength world from the far-field pattern of a scattered wave[END_REF]- [START_REF] Laroche | An inverse approach for ultrasonic imaging from full matrix capture data. application to resolution enhancement in NDT[END_REF]. Let us also mention that this acquisition scheme was also used for particular applications in medical imaging, e.g., blood flow estimation [START_REF] Jensen | Synthetic aperture ultrasound imaging[END_REF].

Due to its finite aperture, each transducer composing a multi-element probe has a limited bandwidth. Consequently, all linear beamforming techniques (SAFT, PWI, TFM) produce images of limited spatial resolution, in particular showing oscillations due to the shape of ultrasonic pulses. In order to improve the resolution and the contrast of ultrasonic images, previous works adopted an inverse problem approach [START_REF] Laroche | An inverse approach for ultrasonic imaging from full matrix capture data. application to resolution enhancement in NDT[END_REF], [START_REF] Tuysuzoglu | Sparsity driven ultrasound imaging[END_REF]- [START_REF] Bueno | Real-time ultrasound image reconstruction as an inverse problem on a GPU[END_REF]. In [START_REF] Laroche | An inverse approach for ultrasonic imaging from full matrix capture data. application to resolution enhancement in NDT[END_REF], we proposed a linear forward model which relates the FMC data to the reflectivity map of the inspected region, accounting for the acoustic responses of the transducers. The corresponding inversion method is able to reconstruct highly resolved images, but at the price of a computationally expensive procedure due to the large size of FMC datasets, which makes it hardly compatible with some industrial NDT applications.

This paper proposes an inversion procedure which aims to preserve the resolution quality of FMC data inversion, while reducing the computation time, by reducing the size of the dataset. More precisely, we transfer the raw FMC data model (i.e., a collection of time-domain signals related to the unknown reflectivity inside the inspected medium) into the TFM beamformed image. A new forward model between the searched reflectivity map and the beamformed image is considered. The reconstruction of the reflectivity map then amounts to an image "deconvolution" problem, where the Point Spread Function (PSF) is spatially variant.

To our best knowledge, the principle of retrieving a reflectivity map from a beamformed image has never been considered in NDT, while it has already raised interest in medical imaging. In [START_REF] Szasz | Beamforming Through Regularized Inverse Problems in Ultrasound Medical Imaging[END_REF], [START_REF] Ozkan | Inverse Problem of Ultrasound Beamforming With Sparsity Constraints and Regularization[END_REF], a regularized procedure is applied to the beamformed image, but no PSF is considered. The deconvolution of beamformed medical images is addressed in [START_REF] Michailovich | A Novel Approach to the 2-D Blind Deconvolution Problem in Medical Ultrasound[END_REF], where a shift-invariant PSF is calibrated from experimental data. In [START_REF] Florea | An Axially Variant Kernel Imaging Model Applied to Ultrasound Image Reconstruction[END_REF], an axially variant PSF is introduced by setting an analytical ad-hoc PSF model for each axial position. The recent contribution [START_REF] Besson | A Physical Model of Non-stationary Blur in Ultrasound Imaging[END_REF] also tackles nonstationary PSFs. It focuses on ultrafast ultrasound imaging, where datasets are obtained using compressed beamforming schemes. Unfortunately, it is hardly tranposable to the context of FMC imaging modality in NDT, since it would rely on repeated computations involving huge FMC datasets.

B. Contribution

The proposed method is specifically dedicated to ultrasonic imaging in NDT. We interpret the TFM beamformed image as the back-projection of the FMC data in the space domain. The resulting forward model directly relates the reflectivity map to the beamformed data through a spatially variant twodimensional (2D) convolution, which accounts for both the time-domain impulse response of the transducers and the geometrical structure of the acquisition process. However, manipulating the set of exact PSFs at all points is not numerically feasible in realistic scenarios. Akin to previous contributions to shift-variant blur approximation in astronomical imaging [START_REF] Gentile | Interpolating point spread function anisotropy[END_REF], [START_REF] Denis | Fast Approximations of Shift-Variant Blur[END_REF], we propose to accurately describe the PSF at all possible locations from the interpolation of a smaller set of reference PSFs, which are computed in two steps. First, synthetic ultrasonic data are generated considering pointwise excitations at different locations in the medium reflectivity, using the ultrasonic response model. Second, the TFM image is computed for each synthetic dataset, thereby identifying the PSF at each specific location.

A second specificity of our approach considers a forward model where the measurement noise affecting the raw FMC data is linearly transformed by the beamforming process. Such a "beamformed noise component" is not white, and its correlation structure is spatially varying. In the proposed inversion method, we account for this specific noise model, using the same methodology that was used to build the non-stationary convolution model. Finally, our inversion strategy relies on the minimization of a penalized least-squares criterion, where the penalization terms favor the presence of sparse regions of limited spatial extent.

The remainder of the paper is organized as follows. Section II introduces the data acquisition procedure, and states the corresponding model relating the FMC data to the unknown reflectivity map. In Section III, the output of the TFM procedure is related to the reflectivity image by a shift-variant 2D convolution model, for which a PSF interpolation model is built. The non-stationarity of the resulting error term is studied in Section IV, where an approximate whitening procedure is built. Then, our inversion procedure is detailed in Section V. The performance of the method in terms of image quality and of computing time is evaluated in Section VI on synthetic data. In Section VII, it is applied on NDT laboratory data acquired from a stainless steel block containing a set of close sidedrilled holes. A discussion is finally given in Section VIII.

II. FULL MATRIX CAPTURE DATA MODEL

The Full Matrix Capture acquisition modality collects the ultrasonic signals that are emitted independently by each transducer and received by all of them. For an array probe containing N el transducers, N 2 el A-scan signals y a,b (t) are recorded, where a and b respectively index the emitter and the receiver. An example of (partial) FMC data is shown in Figure 1, corresponding to the inspection of the specimen that will be described in the experimental Section VII. Typical Ascans extracted from this data set are shown in Figure 2. We can identify frontwall and backwall echoes, as well as the signature of close side-drilled holes generating mixed echoes.

Each A-scan can be modeled as the superposition of timedelayed echoes of the ultrasonic wave that propagates into the medium, where each echo corresponds to the reflection at each point of the medium, weighted by the corresponding reflectivity value [START_REF] Stepanishen | Transient Radiation from Pistons in an Infinite Planar Baffle[END_REF]. Akin to [START_REF] Laroche | An inverse approach for ultrasonic imaging from full matrix capture data. application to resolution enhancement in NDT[END_REF], we adopt the following model for the FMC data:

y = H t o + n t , (1)
where

y ∈ R N 2
el Nt is a column vector gathering all discrete-

III. FORWARD SHIFT-VARIANT CONVOLUTION MODEL ON THE BEAMFORMED IMAGE

In [START_REF] Laroche | An inverse approach for ultrasonic imaging from full matrix capture data. application to resolution enhancement in NDT[END_REF], an inversion procedure was proposed to estimate the reflectivity map o from the data y based on model (1), allowing high-resolution reconstructed images. However, dealing with the huge FMC data y makes this method impractical in many NDT situations, where relatively high-speed inspection may be necessary. In this section, we reduce the data size, by considering the image produced by the Total Focusing Method [START_REF] Holmes | Post-processing of the full matrix of ultrasonic transmit-receive array data for non-destructive evaluation[END_REF] as the new input of our reconstruction method.

A. Total Focusing Method

The Total Focusing Method is the standard beamforming procedure from FMC data. We consider the (x, z) coordinate system (see for example Figure 3), where emitter a and receiver b are respectively located at (x a , 0) and (x b , 0). In the standard configuration where the probe is in contact with the medium to inspect, the travel time from the emitter a to a given location r = (x, z), augmented with the return travel time to receiver b, reads:

τ (r, x a , x b ) = (x -x a ) 2 + z 2 + (x -x b ) 2 + z 2 c , (2)
where c is the ultrasonic velocity. The TFM image at pixel r is then defined by adding all signals delayed by the corresponding times of flight [START_REF] Holmes | Post-processing of the full matrix of ultrasonic transmit-receive array data for non-destructive evaluation[END_REF]:

O TFM (r) := Nel a=1 Nel b=1 y a,b τ (r, x a , x b ) . (3)
In discrete form, this reads:

o TFM = By (4)
where o TFM ∈ R NxNz collects the values of the TFM image in a column vector, and the binary matrix B ∈ {0; 1} NxNz×N 2 el Nt will be called the beamforming operator: for each pixel of the TFM reconstruction grid, the corresponding row in B selects the indices corresponding to the appropriate (rounded) times of flight in each A-scan according to Eq. ( 2).

B. Shift-variant convolution model

We consider the projection of the FMC data model into the spatial domain using the TFM beamforming operator. Combining Eqs. ( 1) and ( 4), we obtain:

o TFM = H s o + n s , (5)
with H s := BH t and n s := Bn t . Matrix H s ∈ R NxNz×NxNz is now a square operator converting the discretized reflectivity map into the TFM image. Considering model ( 5) instead of (1) relies on the idea that the loss of information between the FMC data and the TFM image is limited, and that it can be partly compensated by accurately modeling the resulting PSF and the noise statistical properties. Each column of H s represents the PSF at a given pixel, that is, the TFM image of a pointwise reflector at this location. Due to the acquisition geometry, the shape of the PSF is different for each pixel location, that is, the product H s o represents a shift-variant 2D convolution operation. Figure 3 shows typical PSFs obtained at four different locations in a 30 mm × 30 mm image. The time-domain transducer pulse is a symmetric Gaussian wavelet [START_REF] Demirli | Model-based estimation of ultrasonic echoes. Part I: Analysis and algorithms[END_REF], with central frequency f 0 = 5 MHz and bandwidth ratio (the frequency bandwidth over f 0 ) at -3 dB equal to 0.4. The sound velocity is 5000 m/s, such that the wavelength in the material is λ = c/f 0 = 1 mm. Synthetic data were generated with a 64-element array probe, with an inter-element spacing equal to λ/2 = 0.5 mm. The reflectivity image was composed of four point sources located at (x, z) = (7.5, 7.5), (7.5, 22.5), (22.5, 7.5), and (22.5, 22.5) mm, then yielding the four PSFs at these pixels. As expected, the PSF is wider as the depth increases, and its orientation depends on the pixel location.

C. PSF interpolation

Although the PSF can be computed at any point of the image, computing one PSF per pixel would not be reasonable nor really useful. In [START_REF] Besson | A Physical Model of Non-stationary Blur in Ultrasound Imaging[END_REF], a model accounting for an exact PSF at each location of the reconstruction grid was proposed. However, this model relies on the data propagation operator, which would be too expensive computationally with largesize FMC data. Here, we rather rely on the fact that the PSF variations between close pixels are smooth (see the examples in Figure 3). Consequently, an interpolation model can be used to approximate the PSF shape at any location from a small number of reference PSFs with controlled accuracy.

The interpolation model proposed in this paper is inspired by previous works in Astrophysics [START_REF] Gentile | Interpolating point spread function anisotropy[END_REF], [START_REF] Denis | Fast Approximations of Shift-Variant Blur[END_REF]. In the following, we denote by r n the coordinate vector for the pixel indexed by n in a vectorized image. We first define a set of reference PSFs, say

h ref i , i ∈ [1, N ref ],
corresponding to the PSFs computed at pixels r i . Moreover, we introduce weight maps w i to modulate their contributions to PSF interpolation. That is, the PSF at any pixel r n , say h rn , is approximated as follows: Here, we use the Inverse Distance Weighting (IDW) model [START_REF] Gentile | Interpolating point spread function anisotropy[END_REF]:

h rn (r) ≈ N ref i=1 w i (r n )h ref i (r). ( 6
)
w i (r j ) = δ i,j ∀ j ∈ [1, N ref ], (7a) w i (r n ) = r n -r i -1 2 N ref j=1 r n -r j -1 2 ∀ n ∈ [1, N ref ], (7b)
where • 2 denotes the Euclidean norm.

An example is shown in Figure 4, using the four PSFs introduced in Section III-B as references. Figures 4(a)-(c) respectively show the exact PSF at the center of the image, its interpolation and the corresponding error. The main lobes are quite accurately modeled, and the relative error is approximately 30%. Finally, Figure 4(d) shows the relative PSF interpolation error at each point of the image. It is less than 30% at the center part of the image (where the acoustic field is homogeneous). Due to the inspection geometry, the PSFs at the top and on the sides of the image suffer from distortions that are more difficult to consider. These inconsistencies are related to the TFM process. Indeed, the acoustic field is not homogeneous on the sides of the grid corresponding to the edges of the probe. The PSF spatial variations are also stronger in the near surface area, therefore reducing the interpolation accuracy. However, this has a limited impact since the top of the image is strongly affected by the surface echo in contact inspection, and is usually discarded in the image analysis. Now, let O(r) denote the reflectivity value at pixel with coordinates r. From model ( 6), the spatially-variant 2D con-volution at a given pixel r n reads:

r h r (r n -r)O(r) ≈ r N ref i=1 w i (r)h ref i (r n -r)O(r) = N ref i=1 r h ref i (r n -r)O wi (r) =   N ref i=1 h ref i * O wi   (r n ), (8)
where O wi (r) := O(r)w i (r) and * denotes the 2D convolution. The forward model ( 5) can then be evaluated by summing N ref 2D convolution products, applied to the reflectivity images that have been weighted by the corresponding maps.

IV. COLORED, NON-STATIONARY, NOISE MODEL

The uncertainty term n t in the time-domain signal model (1) stands for model errors and electronic noise affecting FMC data. In NDT, it is reasonable to model them as centered Gaussian white noise if the medium is homogeneous, i.e., if the size of the material grain is small compared to the wavelength.

The covariance matrix of n t reads Cov(n t ) = E[n t n t t ] = σ 2 I, with σ 2 the noise power and I the identity matrix of size N x N z . When dealing with the projected model in the image domain proposed in Eq. ( 5), the noise term n s is filtered by the beamforming operator B, therefore its covariance matrix reads Γ s = σ 2 BB t . From the definition of matrix B in Section III-A, one can show that the element (i, j) of matrix Γ s is:

γ(r i , r j ) =σ 2 Card S ij , with S ij := (a, b) ∈ [1, N el ] 2 such that (9) [τ (r i , x a , x b )F s ] = [τ (r j , x a , x b )F s ] ,
where F s is the sampling frequency, τ (r, a, b) is the flight time defined in Eq. ( 2), and notation [•] rounds to the closest integer. S ij denotes the set of emitter-receiver pairs for which the (rounded) travel times through pixels r i and r j are equal. The noise process n s is clearly not white, but also non-stationary, because the sets S ij depend on the acquisition geometry.

A classical approach to deal with colored noise considers a "pre-whitening" procedure [START_REF] Friedlander | System identification techniques for adaptive signal processing[END_REF], [START_REF] Kessy | Optimal Whitening and Decorrelation[END_REF] in order to retrieve a simpler case, where standard techniques based on the white noise assumption can be used. That is, we seek a whitening operator represented by matrix G s ∈ R NxNz×NxNz , so that model ( 5) becomes:

G s o TFM = G s H s o + G s n s , (10)
where G s n s is a white noise process. Computing the ideal whitening operator would require the computation of the square root of Γ -1 s , which is not easily usable in our case, in particular because Γ s is a huge matrix with no specific structure.

In order to define an efficient whitening procedure that takes into account the spatial variability of the noise, we propose to build an approximate spatially-variant whitening filter, in a similar way to the PSF interpolation procedure defined in Section III-C. First, we build noise covariance functions that are defined locally, at the same pixels as the reference PSFs. To do so, we define the spatial covariance function

γ ref i , i ∈ [1,
N ref ] at the reference pixel r i , from the noise covariance between pixel r i and any other pixel-defined in Eq. ( 9)-as:

γ ref i (r) := 1 2 (γ(r i , r i + r) + γ(r i , r i -r)) for r = 0, (11a)
γ ref i (0) := γ(r i , r i ) + η. (11b)
Eq. (11a) ensures that the covariance function is symmetric, so that the power spectral density, say

F γ ref i (the Fourier transform of γ ref i ), is real-valued. Then, if F γ ref i (ν) takes negative values, the addition of η := -min ν F γ ref i (ν) in (11b) ensures that F γ ref i (ν)
≥ 0 for any spatial frequency ν. Finally, we now define the local whitening filter at pixel r i , whose transfer function is obtained by Wiener filtering:

Fg ref i (ν) := F γ ref i (ν) F γ ref i (ν) + µ , µ > 0, (12)
where g ref i denotes the corresponding impulse response. In order to account for non-stationary noise, we finally propose to replace the spatially-variant PSF model in Section III-C by the whitening PSF model, where reference PSFs h ref i are substituted for

f ref i := g ref i * h ref i . ( 13
)
Figure 5 shows the reference whitening filters g ref i associated to the four PSFs h ref i introduced in Section III-C and represented in Figure 3. Such whitening filters show similar characteristics to those of the PSFs h ref i (similar orientations, slow spatial variations), although with a smaller spatial extension.

V. INVERSION PROCEDURE

We now design an inversion method to retrieve the reflectivity map o from the TFM image o TFM . Such an inverse problem is ill-posed since the TFM image lacks information (in particular at high frequencies), due to the limited bandwidth of the ultrasonic transducers. We adopt a standard regularization strategy, where the ultrasonic image is obtained by minimizing the penalized least-squares criterion:

o s = arg min o J LS (o) + φ reg (o). (14)
The two terms in Eq. ( 14) are discussed hereafter. From the whitening procedure described in Section IV and with similar notations, the whitened forward model ( 10) can be approximated by:

o W TFM = H W s o + n W s , (15)
where:

• o W TFM is the vectorized form of the approximate whitened TFM image O W TFM :=

N ref i=1 g ref i * O wi TFM , where O wi TFM (r) := O TFM (r)w i (r); • H W s o
is the vectorized form of the approximate whitened forward model operator applied to the reflectivity image, that is,

N ref i=1 f ref i * O wi .
By construction, the noise term n W s can now be considered as stationary and white, with covariance matrix equal to I. Then, considering the standard least-squares criterion:

J LS (o) := o W TFM -H W s o 2 2 (16)
is statistically founded.

In NDT, sparsity-enhancing regularization is often considered [START_REF] Laroche | An inverse approach for ultrasonic imaging from full matrix capture data. application to resolution enhancement in NDT[END_REF], [START_REF] Tuysuzoglu | Sparsity driven ultrasound imaging[END_REF], [START_REF] Guarneri | A Sparse Reconstruction Algorithm for Ultrasonic Images in Nondestructive Testing[END_REF]. It is adapted to the detection of few reflectors in homogeneous media, and it improves the resolution by introducing high-frequency information in the solution. Here, we also introduce spatial smoothness in the solution in order to enhance the reconstruction of reflectors with some spatial extension. Hence, we define the penalization function φ reg (o) as:

φ reg (o) = µ 1 o 1 + µ 2 Do 2 2 , µ 1 , µ 2 > 0, ( 17
)
where Do is the image gradient. The 1 -norm term favors sparsity, and the second term enhances spatial smoothness. It can be shown that for any µ 1 ≥ µ max 1 with

µ max 1 := 2 (H W s ) t o W TFM ∞ , (18)
the reconstructed image is identically zero, whatever the value of µ 2 [START_REF] Laroche | An inverse approach for ultrasonic imaging from full matrix capture data. application to resolution enhancement in NDT[END_REF]. In practice, we set µ 1 to a fraction of µ max 1 and µ 2 is set to a small positive value (see the experimental Sections VI and VII).

The minimization of the cost function in Eq ( 14) resorts to convex, non-smooth, optimization, which is performed in this paper with the FISTA algorithm (Fast Iterative Shrinkage Thresholding Algorithm) [START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF], which is a popular iterative algorithm for such kind of problems. It requires the evaluation of the forward H W s • and adjoint (H W s ) t • operators at each iteration. Following the interpolation procedures described in Sections IV and V, both operations essentially resort to computing N ref 2D convolution products. In this work, we exploit the fact that the practical range of each PSF is limited to a small area, so that evaluation of convolution products in the space domain turns out to be more efficient than computations based on two-dimensional Fast Fourier Transform (FFT). The FISTA algorithm is implemented in Matlab and the convolution products are computed using CPU-based (C++) and GPU-based (CUDA) implementation. Indeed, since the reconstructed images are very sparse, forward operations are more efficiently computed with CPU, whereas non-sparse adjoint operations advantageously benefit from GPU implementation.

VI. RESULTS WITH SYNTHETIC DATA

In this section, the proposed method is evaluated on synthetic data. Section VI-A first describes the simulation setup. Then the relevance of the proposed model (PSF interpolation, noise prewhitening) is studied in Subsection VI-B. Its impact on the computation time is addressed in Subsection VI-C.

A. Simulation setting

We consider a synthetic array probe made of 64 elements with inter-element distance equal to λ/2, with wavelength λ = 1 mm. The medium is inspected in contact and contains 25 pairs of close pointwise reflectors, separated horizontally from d F = λ/2 = 0.5 mm. Each pair of reflectors is centered on a grid with coordinates (x, z) ∈ [START_REF] Tanter | Ultrafast imaging in biomedical ultrasound[END_REF][START_REF] Zhang | Defect detection using ultrasonic arrays: The multi-mode total method[END_REF][START_REF] Peng | The sizing of small surface-breaking fatigue cracks using ultrasonic arrays[END_REF][START_REF] Jensen | Synthetic aperture ultrasound imaging[END_REF][START_REF] Ozkan | Inverse Problem of Ultrasound Beamforming With Sparsity Constraints and Regularization[END_REF] mm × [START_REF] Tanter | Ultrafast imaging in biomedical ultrasound[END_REF][START_REF] Zhang | Defect detection using ultrasonic arrays: The multi-mode total method[END_REF][START_REF] Peng | The sizing of small surface-breaking fatigue cracks using ultrasonic arrays[END_REF][START_REF] Jensen | Synthetic aperture ultrasound imaging[END_REF][START_REF] Ozkan | Inverse Problem of Ultrasound Beamforming With Sparsity Constraints and Regularization[END_REF] mm. The pixel size is 0.1 × 0.1 mm 2 . The synthetic reflectivity map equals 1 at the reflectors location and 0 elsewhere. The FMC data are generated with model [START_REF] Drinkwater | Ultrasonic arrays for non-destructive evaluation: A review[END_REF], where the transducer impulse response is a symmetric Gaussian wavelet [START_REF] Demirli | Model-based estimation of ultrasonic echoes. Part I: Analysis and algorithms[END_REF], with central frequency 5 MHz and bandwidth ratio equal to 0.4. White Gaussian noise is added, with signal-to-noise ratio SNR := 10 log 10 H t o 2 2 /N 2 el N t σ 2 = 10 dB. We compare the following deconvolution methods of the TFM image obtained by:

• using one single PSF computed at the center of the grid and assuming white noise in the TFM data model ( 5), named Deconv 1 ; • using four reference PSFs computed on a 3 mm × 3 mm map (see Figure 3) and the PSF interpolation model defined in Section III-C, still with the white noise assumption for model ( 5), named Deconv 4 ; • using the former PSF interpolation model, but now with the colored noise assumption and the whitened model in [START_REF] Peng | The sizing of small surface-breaking fatigue cracks using ultrasonic arrays[END_REF], named Deconv W 4 . In all experiments in this section, the regularization parameter µ 1 is set to 0.8 H t o TFM ∞ , where H denotes the appropriate forward operator in each case, that is, using on single centered PSF for Deconv 1 , using four reference PSFs for Deconv 4 and Deconv W 4 . In the latter case, the whitened model is used and the whitened TFM image o W TFM is considered instead of o TFM . We also consider the inversion procedure proposed in [START_REF] Laroche | An inverse approach for ultrasonic imaging from full matrix capture data. application to resolution enhancement in NDT[END_REF] using the full FMC data set y and based on model (1)the same penalization framework is used-which can be viewed as the reference inverse solution, since there is no information loss in the model. In that case, we similarly set µ 1 = 0.8 H t t y ∞ . Parameter µ 2 is always set to 0, since only pointwise reflectors are considered in this section.

Optimization of criterion ( 14) is run with the FISTA procedure until the relative norm between successive iterates becomes smaller than 0.1%.

B. Impact of PSF variability and colored noise models

The TFM image is shown in Figure 6, and the reconstructed image lines at the five exact depths of the reflectors are shown in Figure 7 (for better visibility, results are split in two panels).

All methods improve the resolution of the TFM image and are able to separate some of the pairs of reflectors. Deconv 1 achieves good results at the center of the imaging area, where the (invariant) PSF model is the most accurate, but misses detections on both edges and at the top of the image. Including the PSF interpolation model with Deconv 4 improves the detection and separation of the reflectors, in particular at low depth where the PSF varies more strongly, but some detections are still missed at the edges of the scene. Finally, the inversion that includes both the PSF interpolation and the whitening procedure is able to correctly reconstruct all the reflectors pairs. This experiment was repeated ten times, with different noise realizations with SNR = 10 dB. The following metrics were computed, where quantities were averaged over the data realizations and over the 25 pairs of reflectors:

• The position error measures the distance between the detected and the true reflectors. We consider that a pointwise reflector is detected if a local maximum in the reflectivity map is distant from the true reflector by less than d F /4, that is, no farther than one of the pixels neighboring the true location. If no reflector is detected, a maximum distance is set to d F /2. • The Peak to Center Intensity Difference (PCID) [START_REF] Fan | A comparison between ultrasonic array beamforming and super resolution imaging algorithms for non-destructive evaluation[END_REF] computes the difference in amplitude between the maximum of the estimated reflectivity due to the two reflectors and the minimum in the area between them, in dB. If the two reflectors are not separated, then the PCID is set to 0. Since the image is sparse, an infinite PCID may be reached. In this case, the value is set to -50 dB. • The resolution capability indicator measures the ratio of resolved pairs of reflectors, considering that the pair is resolved if the corresponding PCID is below -6 dB [START_REF] Fan | A comparison between ultrasonic array beamforming and super resolution imaging algorithms for non-destructive evaluation[END_REF]. • The quadratic error corresponds to the 2 -norm reconstruction error, expressed as a percentage of the 2 norm of the true reflectivity map. • Finally, the amplitude of the estimated reflectivity at the reflectors locations is computed, as a percentage of the true reflectivity value.

Results are gathered in Figure 8, confirming the behavior observed in Figure 7. In general, all metrics are improved by deconvolution methods compared to the TFM-except for amplitude estimation, because regularization introduces some bias in the estimated reflectivity. Deconv 1 , however, performs quite poorly, since less than one half of the pairs of reflectors are resolved, with many positioning errors. Including the PSF variability strongly improve all criteria, which are still improved by accounting for the colored noise model. In particular, the resolving capability is 100%, with less positioning errors and better amplitude recovery than with Deconv 4 . This method performs almost as well as the inversion on the FMC data-although it is much more efficient computationally, as will be shown in the next subsection.

C. Computation times

Figure 9 shows the computation times for the different methods, as a function of the number of transducers and of the image size. Note that the TFM image, the reference PSFs, whitening filters and interpolating weights can be computed prior to optimization. Therefore, the corresponding computation time (typically a few seconds) is not considered here. All computation times have been evaluated on a laptop computer with 16.0 GB RAM, equipped with an Intel Core i7-7820 HK @ 2.90 GHz quad-core processor and a NVIDIA GeForce GTX 1080 graphics card.

As expected, the computational complexity of the different deconvolution methods does not depend on the number of transducers N el , since all computations have been moved to the image space. This is obviously not the case for FMC inversion, where the data size is proportional to N 2 el . For a 251 × 251pixel image, the computation time for Deconv W 4 is always about 2 seconds, while the time ratio between Deconv W 4 and FMC inversion varies from 3 (with 32 transducers) to 400 (with 256 transducers).

For all methods, the computation time increases similarly as a function the image size. As expected, Deconv FMC data have been acquired using the Pioneer platform from TPAC (West Chester, Ohio, USA). The probe is a 128-element probe designed by Imasonic (Voray-sur-l'Ognon, France), with central frequency equal to 3 MHz, and with an inter-element spacing equal to 0.8 mm. At this frequency, the wavelength is λ = 1.88 mm. The Rayleigh criterion, measuring the resolution limit of standard ultrasonic imaging methods [START_REF] Felice | Sizing of flaws using ultrasonic bulk wave testing: A review[END_REF], [START_REF] Simonetti | Multiple scattering: The key to unravel the subwavelength world from the far-field pattern of a scattered wave[END_REF]- [START_REF] Fan | A comparison between ultrasonic array beamforming and super resolution imaging algorithms for non-destructive evaluation[END_REF], is defined as: R := 0.61λ/ sin θ, with tan θ = D/(2z SDH ) and D the aperture of the probe. In this configuration setup, R = 1.64 mm, which corresponds to the distance between the SDHs located at x ≈ 56 mm. The TFM image of the full sample is shown in Figure 11, confirming this resolution limit.

Image reconstruction methods are applied on the area framed in red in Figure 11. The pixel size is 0.05 × 0.05 mm 2 , so that the image dimension is N x × N z = 1373 × 173 pixels. The elementary waveform used in model [START_REF] Drinkwater | Ultrasonic arrays for non-destructive evaluation: A review[END_REF] and in the generation of the reference PSFs is a Gaussian wavelet [START_REF] Demirli | Model-based estimation of ultrasonic echoes. Part I: Analysis and algorithms[END_REF], whose parameters have been estimated by fitting the backwall echoes in the FMC data. It is represented in red on the Ascan signals in µ 2 = 10 -4 H t o TFM 2 , with the forward model operator and the data adapted to the method (see the end of Section VI-A for more details).

The reconstructed images are displayed in Figure 12. For more precise analysis, intensity lines along the x direction are represented in Figure 13. For each x, the intensity is averaged over z SDH ± 0.6 mm. The TFM image shows many artifacts and, in particular, the right side of the image is hardly interpretable: the closest resolved SDHs (with PCID below -6 dB) are the 22 nd and 23 rd ones, located around x ≈ 65 mm and distant from 1.4 mm. The image obtained by Deconv 1 achieves satisfactory detection and separation of the most distant SDHs, up to x ≈ 60 mm. Then, some SDHs are missed. The last resolved SDH pair is located around x ≈ 73 mm (29 th and 30 th ones, distant from 1.0 mm). Better separation is achieved by considering the shift-variant PSF interpolation model, in particular for x ∈ [65, 80] mm. For Deconv 3 , the last resolved SDHs pair is located around x ≈ 76 mm (32 nd and 33 rd , distant from mm). Both Deconv 1 and Deconv 3 , however, show many artifacts below most flawssee in particular for x ≤ 45 mm and for x ∈ [55, 65] mm in Figure 12. both the PSF interpolation model and the whitening procedure, Deconv W 3 strongly reduces these artifacts, which are only visible at the two extremities of the reconstructed image, and still improves the resolution. The last resolved SDH with PCID below -6 dB is also the 33 rd one, but the PCID is around -5.5 dB between the 33 rd and 34 th SDHs, and between the 34 rd and 35 th ones. The closest resolved flaws with the proposed method are separated by 0.7 mm center to center (and by 0.4 mm edge to edge), which is more than four times smaller than the Rayleigh criterion.

VIII. CONCLUSION

Our contribution interprets the TFM image as a linear projection of the FMC data in the space domain. More precisely, it relies on a slowly varying, non-stationary blur model to mathematically describe the relation between the reflectivity map of the media and the TFM image. The spatially varying point spread function depends on both the geometrical and the acoustical configurations of the setup. We proposed a computationally efficient interpolation method to account for these spatial variations. We also proposed a whitening scheme to account for the fact that the TFM image contains spatially correlated, non-stationary noise. The latter is also based on an interpolation model, and its overhead computing cost is negligible.

The resulting non-stationary deconvolution method also relies on a sparsity assumption on the reflectivity map. In the synthetic data case, it produced results of similar quality as the reference inversion of FMC data. In the experimental data case, the method was able to separate close SDHs with a diameter six times smaller than the wavelength and distant from a quarter of the Rayleigh criterion.

In our interpolation model, offline computations (which are necessary to determine the reference PSFs, the whitening filters and the associated weights, for instance) can be performed in advance as far as the geometrical and acoustical configuration is known, which is often the case in industrial contexts. Under our current implementation, the proposed method cannot perform real time reconstructions. However, for a standard configuration of 64 transducers and 251 × 251 pixels, the online computation time is around one second using a rather standard laptop computer, which is still reasonable in some practical applications.

In practice, for more computational efficiency, we could imagine the following inspection strategy: first, the TFM image is computed. Second, if something suspicious is detected on the TFM image, the proposed method is applied to increase the quality of the reconstructed image. Similarly to what was performed in the experimental section of the paper, the reconstruction procedure could then be run only on parts of the inspected piece where high-resolution imaging is required.

Finally, let us mention that the current model assumes the shape invariance of the ultrasonic pulse during the propagation.

Significant distortions could occur in attenuative and dispersive materials [START_REF] Kak | Signal Processing of Broadband Pulsed Ultrasound: Measurement of Attenuation of Soft Biological Tissues[END_REF]- [START_REF] Carcreff | A linear model approach for ultrasonic inverse problems with attenuation and dispersion[END_REF]. Our ongoing research consists in incorporating such distortions in our spatially-varying PSF model, in order to adapt our non-stationary deconvolution method to such complex materials.

Fig. 1 .

 1 Fig. 1. Typical FMC data in NDT, obtained with a 128-element probe (partial). Five blocks are represented, corresponding to the 128 collected signals after emission by transducer 1, 32, 64, 96 and 128.

Fig. 2 .

 2 Fig.2. Some time-domain signals from the FMC data set in Figure1. In red, the estimated elementary waveform that will be used in the model (see Section VII).

Fig. 3 .Fig. 4 .

 34 Fig. 3. Variability of the PSF corresponding to model (5). (a): four PSFs computed at several pixels of the reconstruction grid (centered on the black squares). (b)-(e): zooms on the corresponding zones.

Fig. 5 .

 5 Fig. 5. Impulse responses of the whitening filters g ref i at the same locations as the four PSFs in Figure 3. The absolute value is displayed in logarithmic scale for a better visibility.

Fig. 6 .

 6 Fig. 6. TFM image for a synthetic scene composed of 25 pairs of close reflectors.

Fig. 7 .Fig. 8 .

 78 Fig.7. Reconstruction results on synthetic data. Intensity lines (in dB) of reflectivity profiles at the depths of the five reflectors: TFM image (red line), Deconv 1 (blue stem), Deconv 4 (green stem), Deconv W 4 (magenta stem), and FMC inversion (black stem). The true locations are represented by the cyan crosses.

Fig. 9 . 4 (

 94 Fig. 9. Computation times for Deconv 1 (blue), Deconv 4 (green), Deconv W 4 (magenta), and FMC inversion (black). Top: as a function of the number of transducers, for a 251 2 -pixel image. Bottom: as a function of the image size, with 64 transducers. Both ordinate axes are in logarithmic scale.

1

 1 Figure9shows the computation times for the different methods, as a function of the number of transducers and of the image size. Note that the TFM image, the reference PSFs, whitening filters and interpolating weights can be computed prior to optimization. Therefore, the corresponding computation time (typically a few seconds) is not considered here. All computation times have been evaluated on a laptop computer with 16.0 GB RAM, equipped with an Intel Core i7-7820 HK @ 2.90 GHz quad-core processor and a NVIDIA GeForce GTX 1080 graphics card.As expected, the computational complexity of the different deconvolution methods does not depend on the number of transducers N el , since all computations have been moved to the image space. This is obviously not the case for FMC inversion, where the data size is proportional to N 2 el . For a 251 × 251pixel image, the computation time for Deconv W 4 is always about 2 seconds, while the time ratio between Deconv W 4 and FMC inversion varies from 3 (with 32 transducers) to 400 (with 256 transducers).For all methods, the computation time increases similarly as a function the image size. As expected, Deconv 1 is fastest one. With 64 transducers, accounting for the PSF interpolation and the whitening procedure increases the computation time by a factor of 1.4 for the 101 × 101-pixel image (0.15 s vs. 0.23 s), and by a factor of 3.7 for the biggest, 1-million-pixel, one (32.2 s vs. 119.2 s). Finally, with 64 transducers, Deconv W 4 is approximately 30 times faster than FMC inversion, whatever the image size.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Stainless steel block containing a set of 41 close aligned SDHs (∅ 0.3 mm), with decreasing distance from left to right.

Figure 2 .Fig. 12 .Fig. 13 .

 21213 Fig. 12. Ultrasonic imaging of the stainless steel piece showed in Figure 10. (a): the three reference PSFs used in the interpolation method (the black frames show the considered spatial extension of the kernels). (b): TFM reconstruction. (c): reconstruction by Deconv 1 . (d): reconstruction by Deconv 3 . (e): reconstruction by Deconv W 3 .

ACKNOWLEDGMENT

We would like to thank George D. Connolly from Electric Power Research Institute (EPRI, Charlotte, USA) who kindly lent us the experimental specimen.

This work was partially funded by the French ANRT (Association Nationale Recherche Technologie), project 2017/1083.