The Laplace transform inversion is a well-known ill-conditioned problem and many numerical schemes in literature have investigated how to solve it. In this paper, we revise the Gaver-Stehfest method by using polyharmonic splines augmented with polynomials to approximate the Laplace transform into the numerical inversion formula. Theoretical accuracy bounds for the fitting model are given. Discussions on the effectiveness of the inversion algorithm are produced and confirmed by numerical experiments about approximation errors and inversion results. Comparisons with an existing model are also presented.