Accurate and parallelizable inverse Hankel transform for full waveform inversion

A computationally efficient way to generate an initial model for 3D full waveform inversion is to first invert for 1D velocity profiles, exploiting Hankel transform properties to speed up modeling. The inverse Hankel transform is an improper integral, and knowing at which wavenumbers to sample to numerically evaluate it with sufficient accuracy is non-trivial. In this work, we split the inverse Hankel transform into two integrals. A definite integral which is approximated with numerical integration, and a second improper integral which is treated as a convergent series using the extrema of the integral. The series can be forced to converge quicker via the method of repeated averaging. The location of the extrema correspond to the roots of the Bessel function, which are known and can be compiled into a wavenumber list. The forward modelling for each wavenumber can be then independently carried out, in a massively parallel fashion, on CPU or GPU cores. This way only the number of integration nodes and number of extrema need to be specified a-priori for accurate modeling. This is beneficial in inversion schemes where quality control of the forward modeling steps is impractical. Numerical results using only ten extrema already give low relative errors.