## 2 Problem Analysis and Theory
### Principles of PIRNN Method
#### 2.4.1 PIRNN Forward Modeling  
We performed seismic forward modeling using the first-order pressure-velocity acoustic wave equation in homogeneous media, as shown in Equation (7):  
\[\left\{ \begin{array}{l}\frac{\partial p_{x}(r,t)}{\partial t}=v^{2}(r)\frac{ \partial v_{x}(r,t)}{\partial x}\\ \frac{\partial p_{z}(r,t)}{\partial t}=v^{2}(r)\frac{\partial v_{z}(r,t)}{ \partial z}\\ \frac{\partial v_{z}(r,t)}{\partial t}=\frac{\partial p(r,t)}{\partial z}\\ \end{array} \right., \tag{7}\]  
where \(p_{x}\) and \(p_{z}\) are the pressure components in the x and z directions, respectively, \(v_{x}\) and \(v_{z}\) are the particle velocities in the x and z directions, respectively, and \(r\) represents the spatial position of the wavefield variables; and \(v(r)\) is the acoustic velocity to be inverted. When numerically solving wave equations, a time-stepping method is typically used, as it ensures that the wavefield state at each time step depends on the state of the previous time step. This characteristic aligns closely with the structure of RNNs, which can store information from previous time steps using hidden states [32]. These hidden states can be utilized to carry the wavefield information at each moment, connecting the gradient information across the entire spatiotemporal field. To further derive the time dependence of the wave equation, we discretized Equation (7) using finite differences on a second-order temporal accuracy staggered grid. Taking the x-direction as an example [33]:  
\[\left\{ \begin{array}{l}\frac{v_{x}^{k+\frac{1}{2}}[ix+\frac{1}{2},iz]-v_{x}^ {k-\frac{1}{2}}[ix-\frac{1}{2},iz]}{\Delta t}=\frac{p^{k}[ix+1,iz]-p^{k}[ix,z] }{\Delta x}\\ \frac{p_{x}^{k+1}[ix,iz]-p_{x}^{k}[ix,z]}{\Delta t}=v^{2}[ix,iz]\frac{x+\frac{ 1}{2}[ix+\frac{1}{2},iz]-v_{x}^{k+\frac{1}{2}}[ix-\frac{1}{2},iz]}{\Delta x} \end{array} \right., \tag{8}\]  
Rearranging Equation (8) to obtain the recursive formulas for \(\mathrm{v_{x}}\) and \(p_{x}\) produces:  
\[\left\{ \begin{array}{l}v_{x}^{k+\frac{1}{2}}[ix+\frac{1}{2},iz]=\frac{\Delta t [p^{k}[ix+1,iz]-p^{k}[ix,z]]}{\Delta x}+v_{x}^{k-\frac{1}{2}}[ix-\frac{1}{2}, iz]\\ p_{x}^{k+1}[ix,iz]=v^{2}[ix,iz]\Delta t\frac{v_{x}^{k+\frac{1}{2}}[ix+\frac{1}{2}, iz]-v_{x}^{k+\frac{1}{2}}[ix-\frac{1}{2},iz]}{\Delta x}+p_{x}^{k}[ix,iz] \end{array} \right., \tag{9}\]  
In each time iteration, the pressure \(\mathrm{p^{k}}\) at time \(\mathrm{k\Delta t}\) and particle velocity \(v_{x}^{k-1/2}\) at time \((\mathrm{k-1/2})\Delta t\) are first used to obtain the particle velocity \(v_{x}^{k+1/2}\) at time \((\mathrm{k+1/2})\Delta t\). Using the calculated particle velocity \(v_{x}^{k+1/2}\) at \((\mathrm{k+1/2})\Delta t\) and pressure \(p_{x}^{k}\) at \(\mathrm{k\Delta t}\), the pressure \(p_{x}^{k+1}\) is computed for the next time step. In this process, each time step of solving the wave equation utilizes wavefield data from the previous times \(\mathrm{k\Delta t}\) and \((\mathrm{k-1/2})\Delta t\). This temporal dependency enables seismic forward modeling using an RNN, where each layer of the RNN can be constructed to compute and store wavefield information at a specific time step. This structure can then be used to map the temporal evolution of the wavefield onto the hierarchical structure of the RNN. Specifically, each layer of the RNN represents the wavefield state at a particular time step, including the pressure and particle velocity. Within this framework, each RNN layer stores the current wavefield information and receives the output from the previous layer as the input. This architecture effectively simulates the temporal dependence of the wave equation, where the solution at each time step depends on the state of the previous time step.