Carbonate cave reservoirs have good hydrocarbon-bearing properties and are an important oil carrier. How to accurately identify the nature of fluids in carbonate cave reservoirs is the key to oil exploration and development. In response to this problem, geophysicists have proposed a variety of new technologies for directly identifying fluids. Tang Wenbang  found that the application of frequency difference analysis technology can effectively distinguish the nature of cave filling. Cai Rui  proposed a method for fluid identification of carbonate reservoir based on spectral decomposition technique. This method is based on the sensitivity of seismic waves to the porosity of carbonate reservoirs. Goodway et al.   combined the Biot theory with the Gassmann equation to identify fluids using AVO techniques. Han Gehua  established a carbonate cave reservoir imaging technique based on the characteristics of the Ordovician reservoirs in the Tahe Oilfield, with the pre-stack time migration and the target layer fine processing technology as the core. This technique provides detailed basic information for carbonate cave reservoir prediction.
In the processing and interpretation of seismic data, it is very important to carry out forward modeling research on deterministic underground geological bodies. On the one hand, forward seismic records can be used as the basis for interpretation of seismic data, and the results of data processing can be tested. On the other hand, forward modeling can be used as the basis for inversion studies. Wandler  used physical model techniques to prove that the AVO response of carbonate cave can be used as a discriminating factor for fluid identification. The numerical simulation method has been widely used in forward modeling because of its advantages in efficiency. Yao Yao  applied a non-homogeneous elastic wave equation to the forward model of the carbonate cave reservoir, and analyzed the characteristics of the seismic wave field in detail. Wu Junfeng et al.  used the non-homogeneous elastic wave equation to carry out forward modeling of the collapsed cavity and multi-hole composite cavity model, and systematically analyzed and summarized the forward modeling results. Min Xiaogang et al.  carried out a forward modeling of several different types of carbonate cave reservoirs in Tahe Oilfield, and summarized the reflection characteristics of the caves on the stack and migration sections.
In this paper, the numerical simulation of cross-grid high-order finite difference method is carried out, and the seismic response characteristics of carbonate cave reservoirs filled with different fluids are analyzed by AVO (amplitude versus offset). Three sets of models were designed: the caves varied in width, the caves filled with different solids, and the oil-gas-water model. By using the numerical simulation technology to obtain the seismic records of the three groups of models and analyzing the AVO response characteristics, the fluid properties in the carbonate caverns are identified.
According to the reflection and transmission theory of seismic waves  , when a plane primary wave is incident on the interface of two media with different elastic properties, a reflected primary wave and a reflected secondary wave are generated in the first medium, and a transmitted primary wave and a transmitted secondary wave are generated in the second medium. If the rock reservoir contains different fluids, the seismic wave velocity will be different. This change causes the amplitude of the reflected wave to change regularly with the offset (AVO). The AVO is directly related to the elastic properties on both sides of the interface.
The reflection and transmission coefficient can be expressed by the Zoeppritz equation. Which is a fourth-order matrix equation system composed of medium density, primary wave velocity VP and secondary wave velocity VS, incident angle and reflection angle and transmission angle. The complete Zoeppritz equation is very complicated and the physical meaning is not clear  . In order to be able to visually observe the relationship between amplitude and physical parameters, we used the Shuey approximation of the Zoeppritz equation in the AVO analysis  .
In formula (1):
where R0(θ) is the reflection coefficient when the incident angle is θ, VP is the primary wave velocities difference between the upper and lower layer, is the density difference between the upper and lower layer, VP is the primary wave velocity, is the layer density, is the Poisson’s ratio. When the angle of incidence is small, the Shuey formula can be simplified to:
where P is the zero-offset primary wave reflection coefficient; G represents the reflection coefficient relative to the offset, which is a parameter directly related to the elastic parameter. This approximation is less accurate than the Zoeppritz equation, but is often effective in describing AVO anomalies when the incident angle is small.
3. Model Design
In this paper, three sets of models are designed to study the AVO response characteristics of caves with different fluids.
The first set of models simulates the AVO response characteristics of single cave with different width, including seven sub-models (1-1; 1-2; 1-3; 1-4; 1-5; 1-6; 1-7). The elastic parameters of the underground geological body are the same, and the width changes. The specific parameters of seven sub-models are shown in Table 1.
The second set of model simulates the AVO response characteristics of the cave filled with different solid, this model includes six sub-models (2-1; 2-2; 2-3; 2-4; 2-5; 2-6). The width of the underground geological bodies is the same, and the elastic parameters are changed. The specific parameters are shown in Table 2.
The third model simulates the AVO response characteristics of the caves filled with oil and gas in the same layer, including three horizontal layer and six cave filled oil, gas, water and solid (3-1; 3-2; 3-3; 3 -4; 3-5; 3-6). The specific parameters are shown in Table 3.
Table 1. The parameter of the different horizontal scale model.
Table 2. The parameter of the cave lithology change.
Table 3. The parameter of the oil-gas-water.
Figure 1 shows the schematic diagram of the first and second sets of models and the seismic acquisition geometry. The underground geological bodies of the two models have identical surrounding rock, buried depth and acquisition parameters. The specific acquisition parameters are as follows: single-sided shooting; the receive channels are 120, the gun spacing is 40 m, the track spacing is 20 m, the number of coverage is 60, and the sampling rate is 1 ms.
Figure 2 is a schematic diagram of the third group (oil-gas-water) model and s the seismic acquisition geometry. The model is based on reality, and consists of three layers. The oil, gas, water caves and other three caves filled solid for comparison are located in the third layer. The size, depth and spacing of the six
Figure 1. The 1st and 2nd groups of models and acquisition schematic diagram.
Figure 2. The oil-gas-water model and acquisition schematic diagram.
geological bodies are the same. The specific acquisition parameters are as follows: single-sided shooting; the receive channels are 120, the shot spacing is 40 m, the track spacing is 20 m, the number of coverage is 60, and the sampling rate is 1 ms.
4. Numerical Simulations
The staggered-grid high-order finite difference method of acoustic wave equation was used to carry out numerical simulation  . This method has faster calculation speed and takes up less memory, and is suitable for seismic wave field simulation of simple models. The time difference and spatial difference precision of numerical calculation are 2nd order and 10th order respectively. When performing cross-grid finite difference numerical simulation, the stress and velocity components are located at t +∆t/2 and t. The staggered grid node distribution is shown in Figure 3. Perfect matching layer (PML) boundary absorption conditions are applied during the simulation. The grid spacing is 5 m, the sampling rate is 1 ms, and the acoustic frequency is 30 HZ.
The two-dimensional acoustic wave equation is shown in Equation (3):
where u is pressure, v and vz represent particle velocity in the x and y direction, is medium density, vp is primary velocity.
5. Results and Discussion
Figure 4 shows the single shot record corresponding to seven sub-models in the first set of models. When the width of the underground geological body is 50 m (Figure 4(a)), the seismic response is a weak hyperbolic reflection axis. As the width of the geological body increases (Figure 4(b)-(f)), the diffraction tails produced at both ends of the geological body become more and more obvious. When the width of underground geological body is 5000 m (Figure 4(g)), the seismic response becomes a strong hyperbolic reflection axis. The single shot record of the numerical simulation is consistent with the actual situation.
The CDP data at the center of the seven sub-model geological bodies is extracted to form an angular gathers, according to the angular set we obtain the AVO curve (Figure 5). In Figure 5 the dotted line is the simulation data, the solid line is the theoretical data calculated by the Shuey approximate equation. The simulated data is basically consistent with the theoretical data which prove that the numerical simulation is accuracy.
The seven AVO curves of the first set of models were overlapped in Figure 6, from which we can see that when the width of the geological bodies is 50 m and 100 m, the curve is gentle and almost no change. When the width is greater than 200 m, the characteristic of AVO curve shows a change. When the width is more than 600 m, the variation characteristics of the curve tend to be consistent.
The six AVO curves of the second set of models were overlapped in Figure 7, from which we can see that each curve has the characteristics of negative intercept and positive slope. The absolute value of the amplitude decreases with the increase of the offset, and the phase inversion occurs after a certain angle, which can be classified into the fourth type of AVO response characteristics.
Previous studies have shown that a small amount of gas can significantly reduce the primary wave velocity  , but has no effect on the shear wave velocity,
Figure 3. The staggered grid distribution map of 2D acoustic equation velocity-pressure.
Figure 4. The shot records of the different horizontal scale model.
so the gas content will cause the Vp/Vs to decrease. As can be seen from Figure 7, as the ratio of Vp/Vs decreasing, the absolute value of the intercept of the AVO curve increases.
The third model of oil-gas-water is in accordance with the actual geological conditions. The CDP data at the center point of each cave is also extracted to form an angle gather, and the AVO curves are obtained (Figure 8). Compared with other curves, the AVO curve of Cavity 1 (gas) has obvious differences in response characteristics, so AVO technology can effectively identify gas cave. The difference in AVO characteristics of caves containing water, oil or other solids (Cavity 2-6) is not obvious, and oil and water cannot be identified by AVO alone.
For AVO fluid identification technology, we have the following conclusions: 1) When the width of the geological bodies is 50 m, the AVO curve is smooth and
Figure 5. AVO response of the different horizontal scale model.
Figure 6. AVO response of the different horizontal scale model.
Figure 7. AVO response of the different lithology rock model.
Figure 8. AVO response of the oil-gas-water model.
its characteristics are almost unchanged; When the width is greater than 200 m, the AVO curve shows a significant change with the increase of the incident angle; When the width is larger than 600 m, the variation characteristics of the AVO curve are basically not affected by the lateral scale. 2) When the surrounding rock is constant, the absolute value of the intercept of the AVO curve increases as the ratio of the Vp/Vs decreases. 3) The inclusion of gas in carbonate cave reservoir will greatly affect the seismic wave velocity, and the AVO technique can effectively identify the gas reservoir.
AVO technology has the following limitations. First, it relies heavily on true amplitude, and second, the solution is not unique. So, it is impossible to judge whether the cave filling is water or oil by AVO technology alone.
 Bortfeld, R. (1961) Approximations to the Reflection and Transmission Coefficients of Plane Longitudinal and Transverse Waves. Geophysical Prospecting, 9, 485-502.
 Pei, Z.L. and Mu, Y.G. (2003) High-Order Finite Difference Method Simulation of Staggered Grid for Seismic Wave Propagation in Uneven Media. Journal of China University of Petroleum (Edition of Natural Science), 27, 17-21.