Nowadays, the difficulty and the cost of oil exploration and development are on the rise. Of course it is very important to find new reservoirs, but it is imperative to identify and locate reservoir fluids accurately  . Many scholars have done a lot of related work, such as some scholars combined with Biot theory and Gassmann equation, using AVO technology to identify pore fluid ; Korneev (2004) studied the influence of water or oil in the reservoir on the amplitude of seismic reflected wave and travel time under low frequency conditions, and found that when the reservoir was water filled, the influence was greater, thus finding the basis for oil and water identification under low frequency conditions; Li and Chen (2013) proposed a new reservoir fluid recognition factor based on the physical model and seismic data, which is significance for effectively reducing the multi-solution of reservoir fluid recognition and improving the success rate of reservoir drillingm. Although the above studies are of great significance to the identification and positioning of reservoir fluids   , they do not take into account the influence of reservoir heterogeneity. In fact, due to the influence of sedimentation, diagenesis and tectonics in the formation process, most of the reservoir will have uneven changes in the spatial distribution and various internal properties, that is, the reservoir heterogeneity, and due to the existence of cracks and caves, this heterogeneity will change strongly . Based on the Tesseral2D software, a sand-mud formation model containing water, oil and gas is established. The seismic data of gas filled, water filled and oil filled are obtained by setting the parameters of distance of channel, minimum offset, wavelet frequency and rock velocity and then made the difference respectively . Finally, the average frequency attribute was extracted from the difference data and its characteristics were discussed.
Basic principles of forward earthquake simulation Geometric ray method and wave equation method are two conventional seismic wave field forward modeling methods . The former can accurately calculate the kinematics characteristics of seismic wave such as ray path and travel time . But for some complicated geological structure and lithology information, it is easy to produce blind area. Wave equation simulation method focuses on seismic wave dynamics. It can simulate the seismic wave field characteristics of complex stratum more realistically   . Considering the complex occurrence state of underground rock strata, Tesseral2D, a two-dimensional full-wave field numerical simulation software based on the finite difference method, was used in this study.
2. Numerical Simulation of Sand and Mudstone Formation
The numerical simulation of this project is carried out based on Tesseral numerical simulation software, including three reservoir situations of gas filled, water filled and oil filled respectively, and post-stack time migration is carried out on the collected data to obtain the attribute profile of root-mean-square amplitude difference under the two situations of full water-oil and full gas-oil. In the process of field seismic exploration, the acquisition effect of seismic data is affected by the source frequency, channel spacing, minimum offset, source wavelet frequency, rock mass elastic parameters and other factors, whether these parameters are set reasonably determines the success or failure of seismic exploration. Before simulating and calculating the profile model of sand and mud, it is necessary to use Tesseral2D software to establish the geological model of sand and mud interbedded strata (Figure 1).
There are many kinds of seismic data processing methods, including denoising, deconvolution, dynamic and static correction, velocity analysis, stacking, migration, inversion and seismic monitoring. The main method is deconvolution and superposition. Deconvolution is the main methods to compress seismic wavelet to improve vertical seismic resolution. Superposition processing is the main processing method to enhance the reflected wave signal, suppress regular interference and random disturbance, and improve the seismic signal-to-noise ratio. Migration imaging processing is the main method to realize the spatial homing of the reflected interface. This way restores the wave field characteristics of the reflected wave to improve the seismic horizontal resolution and seismic signal fidelity. Of course, the above eight types of processing methods are closely related to each other. In addition to deconvolution, stacking and migration imaging processing, each of the other methods also has its own unique tasks in the whole process of seismic processing, and plays an irreplaceable role in other processing methods. And in some cases, other processing method could play a key role, for example, in the complicated terrain and surface conditions in the mountains and desert in seismic exploration, static correction and residual static correction processing will play a key role, etc. Therefore, according to the geological characteristics and geological tasks of the exploration area, correct and effective seismic treatment methods should be selected to form a reasonable treatment process and appropriate treatment parameters, so as to obtain high-quality treatment results.
2.1. Definition and Preprocessing of Observation System
The observation system describes the relative position of the shot point and the detection point. The geometric relation between the shot point and the detection point generally defined by the shot point. Geophone definition and source-receiver relationship definition of three parts to complete the observation system. The system include parameters wavelet detection point parameters shotpoint algorithm parameter setting, and then to sand and mud speed model and sandstone model parameter setting, get wave field snapshots and shot records, removal of
Figure 1. Sand and mud interbedded geological model.
the direct wave again. The sandstone model of argillaceous rock It is defined because the argillaceous composition is substantially different from other mineral composition in physical properties.
2.2. Normal Move out Correction
According to the seismic processing methods of major oil companies at home and abroad. The time-distance curve of the reflected wave at the common reflecting point at the horizontal interface of homogeneous medium is T = T0 +x/V when X is greater than 0, the increase of reflection time is not caused by geological factors, and reflection time cannot reflect the depth of the interface. Horizontal superposition will suppress the reflected wave, so it needs to be eliminated. In the common middle point(CMP) channel concentration, the difference between the reflection time of non-zero offset and the reflection time of zero offset is the normal time difference (dynamic time difference). NMO (normal moveout correction) is the subtraction of the normal moveout from the time of the reflected wave to flatten the time-distance curve of the reflected wave (on the same phase axis). For an in-phase axis with a common reflection point channel concentration, its T0 is fixed, and the only variable that determines the normal time difference of the in-phase axis at the position with a offset of X is velocity V, which determines whether the in-phase axis can be leveled.
Principle of horizontal superposition: a standard trace y(i) should be determined after normal moveout correction( NMO ) correction so that the minimum difference between the standard trace and each trace is the average trace of each seismic trace. The normal horizontal superposition is equal weight superposition, and the equal weight superposition may not achieve the ideal superposition effect due to the difference in the mass of each channel. The weighted coefficient of each channel is determined according to the mass of the seismic channel. The weighted coefficient of the channel with good quality is higher, and the weighted coefficient of the channel with bad quality is lower, which is helpful to improve the superposition effect. Formation of standard trace: The standard trace should be a seismic trace that better reflects the characteristics of the superimposed profile, which can be a normal superimposed trace, a superimposed trace set of adjacent common depth points, or a superimposed trace on the superimposed profile after signal enhancement processing. After a common middle point (CMP) channel set passes through the NMO, the in-phase axis of the reflected wave is leveled to achieve the in-phase superposition. After the amplitude of a CMP channel set passes through the NMO, the in-phase axis of the interfering wave cannot be leveled and there is residual time difference. The increase multiple of the stacking amplitude is far less than the N-fold of the horizontal superposition, so that the effective wave is enhanced and the interfering wave is relatively weakened.
2.4. Velocity Analysis
Velocity is a very important parameter in the digital processing and interpretation of seismic data. Velocity parameters not only affect the quality of many links in seismic data processing, but also provide important information about subsurface structure and lithology. The concept of velocity spectrum is derived from the concept of frequency spectrum, which represents the energy of seismic waves at different frequencies. We call the change of the superimposed (or related) energy of seismic waves along different velocities as velocity spectrum. The factors influencing velocity analysis include offset distribution, stacking times, signal-to-noise ratio and excision.
The core problem of migration is to improve the lateral resolution, which requires to reduce the size of the Fresnel zone, and to reduce the size of the Fresnel zone requires wave field continuation. Migration: oblique reflections return to their true underground positions and make diffracted waves converge. The seismic profiles can better demonstrate the spatial morphology and contact relationship of subsurface structures. Function of migration: Migration makes oblique reflections return to their true underground positions, and converges diffused waves, so that seismic profiles can better show the spatial morphology and contact relationship of underground structures, and improve the lateral resolution. Migration method: the earliest migration method is semi-circular scanning superposition method; then diffraction wave stacking technique; and then kirchhoff summation; another migration technique is the finite difference method migration. There is also a class of frequency-wavenumber domain migration methods.
2.6. The Seismic Profile is Directly Differentiated
Numerical simulation is carried out on the division of seismic profile. The amplitude curve is obtained according to the direct difference of seismic profile of gas-oil (Figure 2) and water-oil (Figure 3) and the correlation treatment of time-shift earthquake  . Finally, the analysis is carried out with the actual data, and the conclusion of the variation of relevant amplitude difference is obtained .
3. Results and Discussion
After obtaining the seismic data in the three states of full gas, full water and full oil, the difference has been made (Figures 4-6), and the average frequency attribute has been extracted from the difference data . Then the time-shifted seismic processing is carried out and the amplitude curve after the time-shifted seismic processing is obtained .
In this paper, Tesseral-2D software is used to study the characteristics of fluid substitution in heterogeneous reservoirs . The study shows that after oil-gas
Figure 2. Amplitude difference profile of oil-filled and gas-filled model.
Figure 3. Amplitude difference profile of oil-filled and water-filled model.
Figure 4. Average frequency property of amplitude difference between oil-filled and sand mud model.
Figure 5. Average frequency property of amplitude difference between gas-filled and sand mud model.
Figure 6. Average frequency property of amplitude difference between water-filled and sand mud model.
substitution occurs in the weak heterogeneous reservoir, the amplitude difference from the reservoir to the bottom of the reservoir will experience a mutation process, which is from a positive maximum to a negative maximum; after oil-water substitution occurs in the weak heterogeneous reservoir, the amplitude difference from the reservoir to the bottom of the reservoir will undergo a mutation process, which is from zero to a positive maximum . According to the sudden change of the amplitude difference between the two fluid substitution processes, the top and bottom of the full gas and water reservoirs can be determined, which is of great significance for identifying the nature of the fluid and we can monitoring the layer of the fluid in the process of gas or water flooding. It also helps to determine the depth of the well position and drill favorable reservoirs in the drilling work . The study also shows that the fundamental factor of the sudden change of amplitude difference value is the weak heterogeneity, and the strong scattering effect will cover up this abnormal phenomenon . This paper also finds out the fluid substitution rule of weak heterogeneous reservoir which can be widely applied in production practice. This rule can guide the fluid identification work more accurately and reduce the cost of oil and gas reservoir development  . The application of carbonate strata was not considered , and many uncertain factors, such as measurement error and environmental error, were not taken into account in the process, which should be considered in the future work.
 Castagna, J.P., Sun, S. and Siegfried, R.W. (2003) Instantaneous Spectral Analysis: Detection of Low-Frequency Shadows Associated with Hydrocarbons. The Leading Edge, 22, 129-127.
 Chuai, X.Y., Wang, S.X., Yuan, S.Y., Chen, W. and Meng, X. (2014) A Seismic Texture Coherence Algorithm and Its Application. Petroleum Science, 11, 247-257.
 Goodway, W., Chen, T. and Downton, J. (1997) Improved AVO Fluid Detection and Lithology Discrimination Using Lame Petrophysical Parameters; “λρ”, “μρ”, & “λ/μ Fluid Stack”, from P and S Inversions. 67th SEG Annual International Meeting, Dallas, 2-7 November 1997, 183-186.
 Han, J.J. and van der Baan, M. (2012) Complete Ensemble Empirical Mode Decomposition for Seismic Time-Frequency Analysis. 82th SEG Annual International Meeting, Las Vegas, 4-9 November 2012, 1-5.
 Korneev, V.A., Goloshubin, G.M., Daley, T.M. and Silin, D.B. (2004) Seismic Low-Frequency Effects in Monitoring Fluid-Saturated Reservoirs. Geophysics, 69, 522-532.
 Li, J.Y. and Dvorkin, J. (2012) Effects of Fluid Changes on Seismic Reflections: Predicting Amplitudes at Gas Reservoir Directly from Amplitudes at Wet Reservoir. Geophysics, 77, D129-D140.
 Miao, Z.Y., Chen, J.F., Wang, J., Wang, G., Zhang, C. and Li, W. (2012) Application of Butane Geochemistry of Natural Gas in Hydrocarbon Exploration. Petroleum Science, 9, 455-462.
 Sinha, S.K., Routh, P.S., Anno, P.D. and Castagna, J.P. (2003) Time Frequency Attribute of Seismic Data Using Continuous Wavelet Transform. 73th SEG Annual International Meeting, Dallas, 26th October 2003, 1481-1484.
 Stewart, R.R., Dyaur, N., Omoboya, B., de Figueiredo, J.J.S., Willis, M. and Sil, S. (2013) Physical Modeling of Anisotropic Domains: Ultrasonic Imaging of Laser-Etched Fractures in Glass. Geophysics, 78, D11-D19.
 Sun, L., Zou, C.N., Liu, X.L., Zhu, R. and Wang, X. (2014) A Static Resistance Model and the Discontinuous Pattern of Hydrocarbon Accumulation in Tight Oil Reservoirs. Petroleum Science, 11, 469-480.
 Wang, H.J., Bian, C.S., Liu, G.D., Sun, M. and Li, Y. (2014) Two Highly Efficient Accumulation Models of Large Gas Fields in China. Petroleum Science, 11, 28-38.