## 2 Methodology of the Forward Modelling
### Approximation of AVA Seismic Data Simulation
\tag{5}\]The new layer propagation operator can be expressed as:  
\[\mathbf{B}_{iz}^{new}=\bigg{(}\mathbf{F}_{iz,1}\mathbf{F}_{iz,1}\bigg{(}\frac{1}{2}\bigg{)} \bigg{)}\backslash\bigg{(}\mathbf{F}_{iz+1,2}\big{/}\mathbf{F}_{iz+1,2}\bigg{(}\frac{1}{ 2}\bigg{/}\bigg{)}\bigg{)}. \tag{6}\]  
where matrix \(F_{iz,k}\) and \(E_{iz,k}\) can be expressed in the following forms:  
\[\mathbf{F}_{iz,k}=\mathrm{i}\omega \begin{pmatrix}\beta_{\mathrm{P}}^{iz,k}&\beta_{ \mathrm{S}}^{iz,k}&\beta_{\mathrm{P}}^{iz,k}&\beta_{\mathrm{S}}^{iz,k}\\ -\xi_{\mathrm{P}}^{iz,k}&-\xi_{\mathrm{S}}^{iz,k}&\xi_{\mathrm{P}}^{iz,k}&\xi_ {\mathrm{S}}^{iz,k}\\ -Z_{\mathrm{P}}^{iz,k}&-Z_{\mathrm{S}}^{iz,k}&-Z_{\mathrm{P}}^{iz,k}&-Z_{ \mathrm{S}}^{iz,k}\\ W_{\mathrm{P}}^{iz,k}&W_{\mathrm{S}}^{iz,k}&-W_{\mathrm{P}}^{iz,k}&-W_{ \mathrm{S}}^{iz,k}\end{pmatrix} , \tag{7}\]  
\[\mathbf{E}_{iz,k}=\mathrm{diag}\Big{[}\mathrm{e}^{i\omega\varsigma_{\mathrm{3}}^{ iz,k}}\cdot\mathrm{e}^{i\omega\varsigma_{\mathrm{3}}^{iz,k}z}\cdot\mathrm{e}^{-i \omega\varsigma_{\mathrm{3}}^{iz,k}z}\cdot\mathrm{e}^{-i\omega\varsigma_{ \mathrm{3}}^{iz,k}z}\cdot\mathrm{e}^{-i\omega\varsigma_{\mathrm{3}}^{iz,k}z} \Big{]}. \tag{8}\]  
When \(k=1\), the horizontal slowness is represented by \(s_{1}=\sin\theta/V_{\mathrm{P}_{\mathrm{P}_{\mathrm{z}}}}\). When \(k=2\), the horizontal slowness is represented by \(s_{1}=\sin\theta/V_{\mathrm{P}_{\mathrm{z}-1}}\). \(\beta\), \(\bar{\zeta}\), \(W\), and \(Z\) can be calculated using the following equations:  
\[\beta=\sqrt{\frac{c_{55}\mathrm{s}_{1}^{2}+c_{33}s_{3}^{2}-\rho}{c_{11}s_{1}^{ 2}+c_{33}s_{3}^{2}+c_{55}\big{(}s_{1}^{2}+s_{3}^{2}\big{)}-2\rho}}, \tag{9}\]  
\[\xi=\sqrt{\frac{c_{11}s_{1}^{2}+c_{55}s_{3}^{2}-\rho}{c_{11}s_{1}^{2}+c_{33}s_ {3}^{2}+c_{55}\big{(}s_{1}^{2}+s_{3}^{2}\big{)}-2\rho}}. \tag{10}\]  
In Equations (9) and (10), when \(s_{3}=s_{3\mathrm{P}}\) is calculated, \(\xi=\xi_{\mathrm{P}}\) and \(\beta=\beta_{\mathrm{P}}\) are computed, and when \(s_{3}=s_{3\mathrm{S}}\) is calculated, \(\xi=\xi_{\mathrm{S}}\) and \(\beta=\beta_{\mathrm{S}}\) are computed. \(\theta\) represents the angle of incidence. \[W=c_{55}(\xi s_{1}+\beta s_{3}),Z=\beta c_{13}s_{1}+\xi c_{33}s_{3}. \tag{11}\]  
In Equation (11), when \(s_{3}=s_{3\mathrm{P}},\xi=\xi_{P},\beta=\beta_{P}\) are calculated, \(W_{\mathrm{P}}\) and \(Z_{\mathrm{P}}\) are computed, and when \(s_{3}=s_{3\mathrm{S}}\), \(\xi=\xi_{\mathrm{S}}\), \(\beta=\beta_{\mathrm{S}}\) are calculated, \(W_{\mathrm{S}}\) and \(Z_{\mathrm{S}}\) are computed. \(s_{3}\) and \(s_{1}\) represent the horizontal slowness, and can be expressed as  
\[s_{1}=\frac{\sin(\theta)}{V_{\mathrm{P}}}, \tag{12}\]  
\[s_{3}=\frac{1}{\sqrt{2}}\sqrt{K_{1}\pm\sqrt{K_{1}^{2}-4K_{2}K_{3}}}, \tag{13}\]  
where \(V_{\mathrm{P}}\) represents the P-wave velocity. In Equation (13), \(+\) denotes the vertical slowness of the P-wave \(s_{3\mathrm{P}}\), and \(-\) denotes the vertical slowness of the S-wave \(s_{3\mathrm{S}}\). \(K_{1}\), \(K_{2}\), and \(K_{3}\) are represented as follows:  
\[ \begin{array}{l}K_{1}=\rho\Big{(}\frac{1}{c_{55}}+\frac{1}{c_{33}}\Big{)}+ \frac{1}{c_{55}}\Big{[}\frac{c_{13}}{c_{33}}(c_{13}+2c_{55})-c_{11}\Big{]}s_{1 }^{2},\\ K_{2}=\frac{1}{c_{33}}(c_{11}s_{1}^{2}-\rho),\ K_{3}=s_{1}^{2}-\frac{\rho}{c_ {35}},\end{array}  \tag{14}\]  
where \(c_{ij}\) represents the elastic modulus. According to Equation (1), the relationship between the particle vibration velocity and stress of the first layer and the amplitude of the seismic wave can be expressed as follows:  
\[ \begin{pmatrix}v_{1}\\ v_{3}\\ \sigma_{33}\\ \sigma_{13}\end{pmatrix} =\textbf{F}_{1,1} \begin{pmatrix}U_{\text{P}}^{-}\\ U_{\text{S}}^{-}\\ U_{\text{P}}^{+}\\ U_{\text{S}}^{+}\end{pmatrix} . \tag{15}\]  
The relationship between the particle vibration velocity and stress of the bottom layer, and the reflection and transmission coefficients can be expressed as follows:  
\[A_{2,n}\textbf{r}= \begin{pmatrix}v_{1}\\ v_{3}\\ \sigma_{33}\\ \sigma_{13}\end{pmatrix} =F_{nz,1}