Acoustic radiation force (ARF) imaging techniques have been developed to characterize the viscoelasticity of biological soft materials by measuring the motion excited by the ARF noninvasively. Due to the unknown stress distribution induced by the ARF, it is hard to accurately solve the inverse problem, and simplified single degree-of-freedom (SDF) models have been developed to analyze it. In this study, an inverse finite element (FE) procedure based on a Bayesian formulation is developed. In the presence of measurement noise and prior parameter uncertainties, the Bayesian approach aims to estimate a distribution of the mechanical properties rather than a best-fit value that is generally obtained with an optimization procedure. In this approach, prior parameters, such as ARF distribution, local heterogeneity profile and known mechanical property estimations, are represented as probability distributions. Gaussian process metamodels based on a small number of FE model runs are used in the uncertainty quantification. An analytical approximation technique based on Taylor expansion is applied to solve the integral for the computation of the posterior distribution of the material parameter. Simulation results demonstrate the computational efficiency of the proposed approach and provide practical probability distributions for the estimated material properties. In a comparison study with the SDF models, the Bayesian approach improves the estimation results even when the prior parameters have large uncertainty levels.