Bayesian amplitude variation with angle (AVA) inversion provides a probabilistic way to estimate elastic parameters and quantify uncertainty, but conventional Metropolis–Hastings Markov chain Monte Carlo may mix slowly in multimodal posterior distributions. This study develops a nonlinear AVA inversion workflow combining an improved parallel tempering Markov chain Monte Carlo scheme with the exact Zoeppritz equations. Multiple temperature chains explore the posterior distribution, and probabilistic state exchanges improve communication among chains. The exact Zoeppritz equations serve as the forward operator, reducing errors from linearized approximations, especially at relatively large incidence angles or strong elastic contrasts. Synthetic tests, including noisy angle gathers and a Marmousi model, show that the proposed workflow recovers P-wave velocity, S-wave velocity, and density with higher correlation coefficients than conventional inversion. A field application to a sandstone–mudstone reservoir in the Tarim Basin, northwestern China, further demonstrates improved lateral continuity, with correlation coefficients between the inverted and well-log P-wave velocity, S-wave velocity, and density increasing from 0.79, 0.76, and 0.74 for conventional inversion to 0.88, 0.84, and 0.81, respectively. These results suggest that the proposed method provides a useful probabilistic framework for nonlinear elastic-parameter inversion, although its computational cost remains higher than that of deterministic inversion methods.