Seismic amplitude variation with offset (AVO) inversion is a cornerstone of oil and gas reservoir prediction, enabling the estimation of subsurface elastic parameters and characterization of stratigraphic interfaces. However, balancing inversion accuracy and computational efficiency remains a critical challenge. To address this, we propose a novel probabilistic AVO inversion framework integrating three key innovations. First, we derive a high-precision quadratic approximation for compressional (P-wave) reflectivity by retaining first- and second-order terms from the exact Zoeppritz equations through a perturbation strategy. This approach significantly enhances accuracy compared to conventional linear approximations, particularly in reflecting the true amplitude variation at large angles. Subsequently, to improve lateral continuity and stratigraphic resolution, we introduce an instantaneous phase constraint derived via the Hilbert transform. This constraint leverages phase sensitivity to seismic waveform coherence, ensuring geologically consistent interface characterization during stochastic inversion. Furthermore, we develop a hybrid Markov Chain Monte Carlo (MCMC) algorithm combining adaptive Gibbs sampling with the independent doubly adaptive rejection Metropolis sampling (IA2RMS) method. This framework efficiently samples high-dimensional posterior probability density functions (PDFs) of elastic parameters: Gibbs sampling generates adaptive proposal distributions, while IA2RMS accelerates Markov chain convergence through location- and scale-adjustable proposals. Numerical experiments and field seismic data demonstrate the robustness and feasibility of the proposed probabilistic AVO inversion method.