This paper develops a fast numerical approach to computing spherical-wave reflection coefficients (SWRCs) for layered seabeds, which provides substantial savings in computation time when used as the forward model for geoacoustic inversion of broadband seabed reflectivity data. The approach exploits the Sommerfeld-integral representation of SWRCs as the Hankel transform of a function proportional to the plane-wave reflection coefficient (PWRC), and applies Levin integration to the rapidly oscillating integrand cast as the product of a (pre-computed) media-independent matrix and a vector involving PWRCs at a sparse sampling of integration angles. Compared to conventional Simpson's rule integration for computation of the SWRC, the Levin integration yields speed-up factors of an order of magnitude or more. Further, it results in reduced memory requirements for storage of pre-computed quantities, a desirable property when a graphics processing unit (GPU) is used for parallel computation of SWRCs. The paper applies trans-dimensional Bayesian inversion to investigate the impact of forward modeling in terms of PWRCs and SWRCs on the estimation of geoacoustic parameters and uncertainties. Model comparisons are quantified in simulated- and measured-data inversions by comparing the estimated geoacoustic parameters to the true parameters or core measurements, respectively, and by calculating the deviance information criterion for model selection.