The classical pre-stack seismic inversion technique uses the Zoeppritz equation and its simplified versions to calculate the PP and PS reflection coefficients at different incidence angles, aiding in inverting the subsurface velocity and density parameters. Despite its widespread application, the amplitude variation with angle (AVA) inversion based on the Zoeppritz equation has limitations regarding the accuracy. The AVA neglects transmission losses and the effects of multiple reflections during seismic wave propagation, resulting in reduced resolution. In contrast, the propagation matrix theory offers a comprehensive range of reflection coefficients for P- and S-waves in multilayered media at arbitrary incidence angles, thereby theoretically enhancing the inversion accuracy. However, the seismic responses obtained using this method exist in the slowness–frequency domain and require constant slowness for consistency along a profile. This assumption is violated when variations in the P-wave velocity occur within the subsurface, affecting the incidence angle of propagating seismic waves. This study modifies the propagation matrix theory to compute AVA seismic responses and applies it to pre-stack multiparameter inversion. The effectiveness of the modified method was validated by deriving theoretical AVA seismic responses and comparing them to solutions from a typical layered media model. The modified theory was also employed for seismic pre-stack inversion. Numerical simulations and field data tests demonstrated that the new propagation matrix method offers a high accuracy and stability.