Full Waveform Inversion using the Wasserstein metric for ultrasound transducer array based NDT

Abstract

Ultrasonic imaging methods often assume linear direct models, while in reality, many nonlinear phenomena are present, e.g. multiple reflections. A family of imaging methods called Full Waveform Inversion (FWI), which has been developed in the field of seismic imaging, uses full acoustic wave simulations as direct models, taking into account virtually all nonlinearities, which can ultimately enhance the accuracy of ultrasonic imaging. However, the problem of cycle skipping -- the existence of many local minima of the Least Squares (L2) misfit function due to the oscillatory nature of the signals -- is worsened when FWI is applied to ultrasound data because of a lack of low-frequency components. In this paper, we explore the use of the squared Wasserstein (W2) Optimal Transport Distance as the metric for the misfit between the acquired and the synthetic data, applying the method to Nondestructive Evaluation with ultrasonic phased arrays. An analytical continuous time-domain derivation of the adjoint acoustic field related to the W2 misfit is presented and used for the computation of the gradients. To cope with the computational burden of FWI, we apply a low-memory strategy that allows for the computation of the gradients without the storage of the full simulated fields. The GPU implementation of the method (in CUDA language) is detailed, and the source code is made available. Six prototypical cases are presented, and the corresponding sound speed maps are reconstructed with FWI using both the L2 and the W2 misfit functionals. In five of the six cases, the pixel-wise sum of squared errors obtained with W2 was at least one order of magnitude lower than that obtained with W2, with an increase in the gradient computation time not exceeding 2\%. The results highlight both the adequacy of the W2 misfit for ultrasonic FWI with phased arrays and its computational feasibility.

Figures (14)

• Figure 1: FWI algorithm.
• Figure 2: Cycle-skipping phenomena in $L^2_2$ misfit for time-shifted oscillatory signals.
• Figure 3: Division of the simulation domain in subdomains to fully utilize CUDA thread blocks. The purple area represents one subdomain, which have its pressure values calculated by a $32 \times 32$ thread block. The pink area represents the memory containing the pressure values necessary for this calculation, which is $2N$ pixels wider and higher than the subdomain, where $N$ is the spatial Laplacian precision. This pink area is copied only once to the SM shared memory, minimizing RAM access.
• Figure 4: Fourth-order 2-D Laplacian approximation cross stencil. To approximate the Laplacian in the point $(x,y)$ represented as the red circle, its neighbours in each dimension are used. The purple positions are used to approximate $\partdev{^2u}{y^2}$, and the yellow positions are used to approximate $\partdev{^2u}{x^2}$.
• Figure 5: Simulating fields backward to parallelize gradient calculation. Both schemes first calculate the direct field (1a and 1b) to obtain the metric value and adjoint source signals. At the top, after the first simulation (1a), the adjoint field (2a) is simulated to obtain the starting conditions for the reversed adjoint field (3a), which is run in parallel with the direct field again (4a). At the bottom, the final conditions of the first simulation (1b) are used as the starting conditions for the second simulation (2b), which runs in parallel with the third (3b). On both, the parallelized fields are used to obtain the gradient. Noriega2018.