## 1 Introduction
### Theory and Algorithm  
In an isotropic acoustic medium, the wave equation is as follows [8; 10]:  
\[ \begin{cases}\left(\frac{1}{v_{0}^{2}}\frac{\partial^{2}}{\partial t^{2}}-\rho _{0}\
abla\frac{1}{\rho_{0}}\cdot\
abla\right)p_{0}(\mathbf{x};t;\mathbf{x}_{s })=\delta(\mathbf{x}-\mathbf{x}_{s})\delta(t)\\ d(\mathbf{x}_{r};t;\mathbf{x}_{s})=p_{0}(\mathbf{x}=\mathbf{x}_{r};t;\mathbf{x }_{s})\end{cases}  \tag{1}\]  
where \(v_{0}\) and \(\rho_{0}\) are velocity and density, respectively; \(p_{0}(\mathbf{x};t;\mathbf{x}_{s})\) is the pressure wavefield at any location of \(\mathbf{x}\) due to a source located at \(\mathbf{x}_{s}\); and \(d(\mathbf{x}_{r};\mathbf{x}_{s};t)\) is the recorded data, i.e., the measured value of the pressure wavefield at a receiver position \(\mathbf{x}=\mathbf{x}_{r}\). For a second medium with small perturbations in velocity (\(\delta v\)) and density (\(\delta\rho\)) compared to the previous medium, the velocity and density are \(v_{0}+\delta v\) and \(\rho_{0}+\delta\rho\), respectively. The pressure wavefield in the perturbed medium \(p_{0}+\delta p\) satisfied the same acoustic wave equation as follows:  
\[ \begin{cases}\left(\frac{1}{(v_{0}+\delta v)^{2}}\frac{\partial^{2}}{\partial t ^{2}}-(\rho_{0}+\delta\rho)\
abla\frac{1}{\rho_{0}+\delta\rho}\cdot\
abla \right)(p_{0}+\delta p)=\delta(\mathbf{x}-\mathbf{x}_{s})\delta(t)\\ \delta d(\mathbf{x}_{r};t;\mathbf{x}_{s})=\delta p(\mathbf{x}=\mathbf{x}_{r}; t;\mathbf{x}_{s})\end{cases}  \tag{2}\]  
where \(\delta p\) is the wavefield perturbation, and the observed data perturbation is represented by \(\delta d\). Based on the Born approximation, the following equation for \(\delta p\) can be derived by subtracting Equation (1) from Equation (2):  
\[\left(\frac{1}{v_{0}^{2}}\frac{\partial^{2}}{\partial t^{2}}-\rho_{0}\
abla \frac{1}{\rho_{0}}\cdot\
abla\right)\delta p(\mathbf{x};t;\mathbf{x}_{s}) \approx\left(\frac{2\delta v}{v_{0}^{3}}\frac{\partial^{2}}{\partial t^{2}}- \left(\
abla\frac{\delta\rho}{\rho_{0}}\right)\cdot\
abla\right)p_{0}( \mathbf{x};t;\mathbf{x}_{s}) \tag{3}\]  
By following [PERSON] et al. [10] and using the asymptotic approximation, we obtained the ray-based relationship between \(\delta d\), \(\delta v\) and \(\delta\rho\):  
\[\delta d(\mathbf{x}_{r};\omega;\mathbf{x}_{s})=-\int\frac{2\omega^{2}}{v_{0}( \mathbf{x})^{2}}\bigg{(}\frac{\delta v(\mathbf{x})}{v_{0}(\mathbf{x})}+\text{ cos}^{2}\theta\frac{\delta p(\mathbf{x})}{\rho_{0}(\mathbf{x})}\bigg{)} \times A(\mathbf{x}_{s};\mathbf{x}_{r};\mathbf{x})e^{i\omega T(\mathbf{x}_{s}; \mathbf{x}_{r})}d\mathbf{x} \tag{4}\]  
where \(T(\mathbf{x}_{s};\mathbf{x}_{r};\mathbf{x})=\tau(\mathbf{x};\mathbf{x}_{s})+ \tau(\mathbf{x}_{r};\mathbf{x})\) and \(A(\mathbf{x}_{s};\mathbf{x}_{r};\mathbf{x})=A(\mathbf{x};\mathbf{x}_{s})A( \mathbf{x}_{r};\mathbf{x})\) represent the traverting summation and the amplitude product of the Green's function from the source location \(\mathbf{x}_{s}\) to the image point \(\mathbf{x}\) and reflected back to the receiver location \(\mathbf{x}_{r}\), respectively, and \(\theta\) is the subsurface reflection angle. Equation (4) represents the forward modeling formulation under a high-frequency assumption. The detailed derivations of Equations (3) and (4) can be found in [10]. Following a similar method to [5] and [8], we can invert Equation (4) for the composite model parameter perturbation, \(\frac{\delta v(\mathbf{x})}{v_{0}(\mathbf{x})}+\cos^{2}\!\theta\frac{\delta p( \mathbf{x})}{\theta v_{0}(\mathbf{x})}\). In 2D, the composite model parameter perturbation can be described in terms of the perturbed wavefield as follows (detailed derivation is given in Appendix A):  
\[