Journal of Earth Science  2018, Vol. 29 Issue (6): 1390-1397 PDF     0
AVA Simultaneous Inversion of Prestack Seismic Data Using Particle Swarm Optimization
Jin Zhang1,2,3, Peng Shen1, Weina Zhao1, Xubing Guo1, Wang Xing1, Song Chen1, Xiugang Xu1,2,3
1. College of Marine Geosciences, Ocean University of China, Qingdao 266100, China;
2. Key Lab of Submarine Geosciences and Prospecting Techniques, Ministry of Education, Qingdao 266100, China;
3. Qingdao National Laboratory for Marine Science and Technology, Qingdao 266200, China
ABSTRACT: A new prestack AVA simultaneous inversion using particle swarm optimization algorithm is proposed, which can obtain the elastic parameters such as P-wave and S-wave impedance from P-wave reflection data simultaneously. Compared with the conventional AVA inversion based on generalized linear technique, this method does not depend on the initial model and can reach the global minimum. In order to increase the stability of the inversion, low-frequency trends of P-wave and S-wave impedances are built into the inversion. This method has been successfully applied to synthetic and field data. The estimated P-wave and S-wave impedances can be combined to derive other elastic parameters, which are sensitive for lithology identification and fluid prediction.
KEY WORDS: AVA simultaneous inversion    particle swarm optimization (PSO)    low-frequency impedance model

0 INTRODUCTION

Nowadays, lithology and subtle reservoirs have become the main targets for oil and gas exploration. Many reservoirs have no obvious characteristics in P-wave impedance, which leads to the failure of using traditional poststack inversion. However, different combinations of lithology and fluid parameters have different AVO/AVA responses, and thus the lithology and fluid can be predicted by the variation of amplitude versus offset (Gray and Anderson, 2000; Goodway et al., 1997). Now, prestack seismic inversion combined with drilling and petrophysical analysis has become the main technique for reservoirs prediction.

Conventional prestack inversion involves two separate steps. The first step is to achieve the reflection coefficient of each common-angle gather and the second step is to calculate P-wave velocity, S-wave velocity, density or other elastic parameters using AVO approximation. In this kind of inversion, the common-angle gathers are inverted individually without fully utilizing their mutual difference as constraint, which worsens the reliability and stability of the inversion. In last decades, prestack AVA simultaneous inversion technology has made great progress. It can simultaneously invert multiple prestack gathers to acquire the elastic parameters such as P-wave velocity, S-wave velocity and density, and then gets other elastic parameters such as Poisson's ratio, Young's modulus, Lame's constant and so on, to carry out reservoir classification and fluid prediction. Prestack simultaneous inversion has been addressed by many scholars. Ma (2001) developed a method for estimating P-wave impedance (IP) and S-wave impedance (IS) simultaneously. Hampson et al. (2005) proposed a simultaneous inversion by utilizing the relationship between IP, IS and density. Bruun et al. (2005) used prestack simultaneous inversion to improve the reliability of porosity prediction. Contreras et al. (2006) successfully applied simultaneous inversion for lithology prediction of deepwater reservoirs in gulf of Mexico.

In essence, seismic inversion is an optimization problem from mathematical viewpoint, so we can choose either linear algorithms based on partial derivatives or nonlinear algorithms based on heuristic random searching algorithms. Many works have been done on prestack inversion using generalized linear inversion (GLI) technique (Pan et al., 1994; Tarantola, 1986). The advantage of the linear methods is that the convergence speed of the objective function is very fast. However, the linear methods strongly depend on the initial model and easily fall into the local minimum without priori information. The geophysical problems are nonlinear issues, therefore theoretically the best way to solve them is to use the fully nonlinear optimization algorithms, such as simulated annealing, genetic algorithm and so on. The advantage of the fully nonlinear optimization methods over GLI method is that they can find the global optimal solution and have nothing to do with the initial model. Mallick et al. (2000) and Mallick (1995) implemented a model-based AVO inversion using a genetic algorithm. Ma (2002) used simulated annealing to perform prestack AVA simultaneous inversion. Misra and Sachhi(2008, 2007) published papers on global optimization methods in AVO inversion. Contreas et al. (2005) studied AVA stochastic inversion for sand body prediction, which improves the resolution of sand body recognition.

In this paper, AVA simultaneous inversion using particle swarm optimization (PSO) has been developed. The idea of PSO comes from simulating the foraging behavior of birds. Each bird is regarded as a particle, optimizing its position and velocity through their collective cooperation in each iteration, and finally the swarm follows the best particle in the solution space. PSO is proposed by Kennedy and Eberhart (1995) and has been widely used in multi-objective optimization, data clustering, pattern recognition and other engineering fields. The facts of its fast convergence speed, simple programming, less parameter setting and especially the natural real-coded features make it suitable for dealing with real number optimization problems. PSO has been applied in geophysical inversion for many years, however, it hasn't been introduced in seismic AVA simultaneous inversion. In this paper, we use the PSO method to estimate IP and IS from the P-wave reflection data simultaneously based on Fatti's AVO approximation. In each iteration of inversion, IP and IS are randomly perturbed guided by PSO to minimize the error between the synthetic records and the observed data. In order to increase the stability of the inversion, constraints such as the low-frequency impedance model are built into inversion. The estimated IP and IS can be combined to derive other elastic parameters for lithology and fluid identification. Both the synthetic and field data have comfirmed the validity of this AVA simultaneous inversion.

1 THEORY OF PARTICLE SWARM OPTIMIZATION

In PSO algorithm, every solution to the optimization problem can be regarded as a particle in the searching space. At first, PSO generates an initial swarm, that is, initialize a group of particles randomly in a feasible solution space and each particle is evaluated by the objective function. Then each particle flies at a certain speed in the searching space, which is dynamically adjusted in accordance with the speed of its own and other particles. Usually the particles will follow the motion of the optimal particle and get the best solution after iterations. In each iteration, each particle will track two extremes, one is pbest, the optimal location corresponding to the solution for the particle itself so far, which is the particle's own flight experience, and the other is gbest, the optimal location corresponding to the solution for the whole particles, which is the whole swarm's flight experience. Each particle follows this rule and finally gets to the same position after iterations, which is the best solution to the optimization problem.

The mathematical course of PSO can be described as follows: There is an n-dimensional searching space with m particles in it. Each particle's position represents a potential solution to the objective function. That is, their position vector xi=(xi1, xi2, ···, xin) constitutes a swarm X=(x1, x2, ···, xm). If we define the speed of each particle as vi=(vi1, vi2, ···, vin), the optimal location of the i-th particle will be $p_{{\rm{best}}}^{i} = ({p_{i1}}, {p_{i2}}, \cdots, {p_{in}})$ and the global optimal location of the whole swarm will be ${g_{\rm{best}}} = ({g_1}, {g_2}, \cdots, {g_n})$ . Each particle updates its velocity and position according to the following equations

 \begin{align} {v_{id}}{\text{(}}t + 1{\text{)}} = w{v_{id}}{\text{(}}t{\text{)}} + {c_1}{r_1}\left({{p_{id}}{\text{(}}t{\text{)}} - {x_{id}}{\text{(}}t{\text{)}}} \right) + \hfill \\ \;\quad \quad \quad \quad \, {c_2}{r_2}\left({{g_{id}}{\text{(}}t{\text{)}} - {x_{{\rm{i}}d}}{\text{(}}t{\text{)}}} \right) \hfill \\ \end{align} (1)
 ${x_{id}}{\text{(}}t + 1{\text{)}} = {x_{id}}{\text{(}}t{\text{)}} + {v_{id}}{\text{(}}t + 1{\text{)}}\begin{array}{*{20}{c}} {}&{1 \leqslant i \leqslant m\begin{array}{*{20}{c}} {}&1 \end{array}} \end{array} \leqslant d \leqslant n$ (2)

where w is the inertial weighting factor, which decreases linearly with the number of iterations, to make the algorithm have a good global search capability at the beginning and a good local search performance at the end. c1 and c2 are positive constants called learning factors or acceleration factors. c1 controls the step length for the particle to fly to its best position so far and c2 controls the step length for the particle to fly to the global best position. According to a number of parameter experiments by Kennedy(1998, 1997), the sum of c1 and c2 should be around 4.0, when the algorithm has better behavior. r1 and r2 are random numbers within [0, 1] and t represents the number of current iteration. In order to avoid the velocity of particles not too large for a particle flying over the optimal position, it is necessary to limit the maximum speed of particles. We define the maximum speed of each particle as Vmax, and if vid > Vmax, we set vid=Vmax; if vid < -Vmax, we set vid= -Vmax. The initial positions and velocities of the particles are generated randomly and then iterated according to the above two equations until a satisfactory solution is found.

Equation (1) consists of three parts, the first part represents the particle's previous velocity, which ensures the particle's flight; the second part is a vector from the particle's current position to its own best position; and the last part is a vector from its current position to the best position of the swarm. These three parts together determine the searching capability of particles. The first part balances the global and local searching capabilities, the second part has a strong global searching capability to avoid local minimum, and the third part reflects the information sharing between the particles. In each iteration, each particle updated its speed according to Eq. (1), then fly to the new position according to Eq. (2).

PSO is an efficient global optimization algorithm, which is easily to be implemented. Next, we will use PSO to accomplish AVA simultaneous inversion.

2 AVA SIMULTANEOUS INVERSION SCHEME

Different from the conventional AVO/AVA inversion, AVA simultaneous inversion combines the reflection coefficients acquiring and elastic parameters inversion into one step, which then becomes a global optimization problem. Now, we use PSO algorithm to implement AVA simultaneous inversion and add constraints into the objective function to enhance the efficiency and stability of the inversion.

Fatti et al. (1994) proposed the following AVO approximation with IP, IS and density as the variables

 $\begin{gathered} R\left(\theta \right) = \frac{1}{2}\left({1 + {{\tan }^2}\theta } \right)\frac{{\Delta {I_{\text{P}}}}}{{{I_{\text{P}}}}} - 4{\left({\frac{{{V_{\text{S}}}}}{{{V_{\text{P}}}}}} \right)^2}{\sin ^2}\theta \frac{{\Delta {I_{\text{S}}}}}{{{I_{\text{S}}}}} - \hfill \\ \quad \quad \;\;\, \, \, \left[ {\frac{1}{2}{{\tan }^2}\theta - 2{{\left({\frac{{{V_{\text{S}}}}}{{{V_{\text{P}}}}}} \right)}^2}{{\sin }^2}\theta } \right]\frac{{\Delta \rho }}{\rho } \hfill \\ \end{gathered}$ (3)

We replace VS/VP with IS/IP, so now the reflection coefficient only depends on IP, IS, density and incident angle. Here, IS/IP is no longer derived from the background IP and IS, but generated during each iteration. When IS/IP is close to 0.5 and the angle is small, this equation can be expressed as

 $R\left(\theta \right) = \frac{1}{2}\left({1 + {{\tan }^2}\theta } \right)\frac{{\Delta {I_{\text{P}}}}}{{{I_{\text{P}}}}} - 4{\left({\frac{{{I_{\text{S}}}}}{{{I_{\text{P}}}}}} \right)^2}{\sin ^2}\theta \frac{{\Delta {I_{\text{S}}}}}{{{I_{\text{S}}}}}$ (4)

Among them,

 $\frac{{\Delta {I_{\rm{P}}}}}{{{I_{\rm{P}}}}} = \frac{{I_{\rm{P}}^i - I_{\rm{P}}^{i - 1}}}{{I_{\rm{P}}^i + I_{\rm{P}}^{i - 1}}},\;\frac{{\Delta {I_{\rm{S}}}}}{{{I_{\rm{P}}}}} = \frac{{I_{\rm{S}}^i - I_{\rm{S}}^{i - 1}}}{{I_{\rm{S}}^i + I_{\rm{S}}^{i - 1}}}$ (5)
 $\frac{{{I_{\text{S}}}}}{{{I_{\text{P}}}}} = \frac{{I_{\text{S}}^i + I_{\text{S}}^{i - 1}}}{{I_{\text{P}}^i + I_{\text{P}}^{i - 1}}}$ (6)

We can estimate IP and IS from the prestack seismic gathers according to Eq. (4). The basic assumption is that the layers are approximately horizontal at common depth point, and each layer is described by IP and IS. The solution to be optimized contains n IP and n IS, where n is the number of the layers.

To implement AVA simultaneous inversion, we need to create an initial model at first, which can be random generated or obtained from smoothed sonic logging curve. For each layer expressed by IP and IS, we can calculate the reflection coefficient for each angle according to Eq. (4), and then obtain the synthetic common-angle gathers by convolving the reflection coefficients with the known wavelet. These synthetic data are compared with the observed seismic data to calculate the errors between them. Then we randomly perturb the IP and IS to create a new layer model, generate new synthetic data, and then compare with the observed data. This process iterates until the synthetic data match the observed data. We use PSO to accomplish the optimization. Once the global minimum of the objective function is found, the corresponding model parameters are considered as the resultant model.

We adopt the following objective function

 $\begin{gathered} E = \frac{{\sum\limits_{j = 1}^m {\sum\limits_{i = 1}^n {{{\left({S_{{\text{obs}}}^{ij} - S_{\bmod }^{ij}} \right)}^2}} } }}{{\sum\limits_{j = 1}^m {\sum\limits_{i = 1}^n {{{\left({S_{{\text{obs}}}^{ij}} \right)}^2}} } }} + \frac{{\sum\limits_{i = 1}^n {{{\left({I_{{\text{P}}\, {\text{pri}}}^i - I_{{\text{P}}\, \bmod }^i} \right)}^2}} }}{{\sum\limits_{i = 1}^n {{{\left({I_{{\text{P}}\, {\text{pri}}}^i} \right)}^2}} }} + \hfill \\ \quad \;\;\, \frac{{\sum\limits_{i = 1}^n {{{\left({I_{{\text{S}}\, {\text{pri}}}^i - I_{{\text{S}}\bmod }^i} \right)}^2}} }}{{\sum\limits_{i = 1}^n {{{\left({I_{{\text{S}}\, {\text{pri}}}^i} \right)}^2}} }} \hfill \\ \end{gathered}$ (7)

where $S_{{\text{obs}}}^{ij}$ is the observed seismic amplitude at time index i and channel index j, and $S_{\bmod }^{ij}$ is the synthetic seismic amplitude at time index i and channel index j. $I_{{\text{P}}\, {\text{pri}}}^i$ is the trend of the low-frequency IP at time index i. $I_{{\text{P}}\, \bmod }^i$ is IP of model at time index i, which is the smoothed result of inverted $I_{{\text{P}}\, }^i$ during each iteration. $I_{{\text{S}}\, {\text{pri}}}^i$ is the trend of the low-frequency IS at time index i. $I_{{\text{S}}\, \bmod }^i$ is IS of model at time index i, which is the smoothed result of inverted $I_{{\text{S}}\, }^i$ during each iteration. n is the number of the seismic trace samples, and m is the number of seismic traces. The objective function involves three terms. The first term is to minimize the amplitude errors between the synthetic and observed seismic data. The second and third term control the inverted IP and IS, to make their low-frequency trend match the corresponding smoothed sonic logging curve.

In this way, Eq. (7) has transformed prestack AVA inversion into global optimization problem. There are 2n parameters needed to be globally optimized, which are n IP and n IS. It should be noted that the amplitude of the seismic data must be scaled to the level matching the synthetic record obtained from the sonic logging data before inversion.

The solution of seismic inversion is band-limited, lacking of low-frequency and high-frequency information, so without suitable constraints, the inversion will become very unstable, and even lead to a great deviation from the true value. In order to reduce non-uniqueness and make the inversion stable, we apply the following constraints into the algorithm.

1. A priori IP and IS information derived from the well log data or NMO velocity is applied to the second and third term of the objective function. This makes the solution approach the true value of the low-frequency model. With this constraint, the low-frequency components will be added to the impedance solution, which are missing in the band-limited seismic data.

2. The search boundaries are set up for each parameter. Good boundaries can reduce the non-uniqueness of the solution and minimize the computation time. The parameter boundary can be obtained from the low-frequency trend of the sonic logging curve or NMO velocity. It can neither be too wide nor too narrow. Too narrow bound may fail to achieve the true impedance model, while too wide bound will waste the computation time and easily lead to the instability of inversion. We set the width of the searching corridor as ±30% of the low-frequency model of the sonic logging curve. That is, in the inversion process, IP and IS are randomly perturbed in the range of 70%–130% of low-frequency model.

We do AVA simultaneous inversion with PSO algorithm as follows.

1. Extract at least three common-angle gathers from the forward modeling data or field data as observed seismic amplitude, that is, $S_{{\text{obs}}}^{ij}$ in Eq. (7).

2. Extract the wavelets of each common-angle gathers.

3. Initiate the particle swarm. Each particle stands for a potential solution of the model. For a n-layer model, the position vector xid of each particle represents 2n potential solutions, where the former n xid is $I_{{\text{P}}\, }^i$ and the latter n xid is $I_{{\text{S}}\, }^i.$

4. Calculate the reflection coefficient of each commonangle gather, and then convolve with the wavelets to obtain the synthetic records, that is, $S_{\bmod }^{ij}$ in Eq. (7).

5. Calculate the objective function E, where $I_{{\text{P}}\, {\text{pri}}}^i$ and $I_{{\text{S}}\, {\text{pri}}}^i$ can be obtained from the smoothed sonic logging data, which compensate the low-frequency of seismic data. $I_{{\text{P}}\, \bmod }^i$ and $I_{{\text{S}}\, \bmod }^i$ are generated from each particle, which are the smoothed results of the $I_{{\text{P}}\, }^i$ and $I_{{\text{S}}\, }^i.$ Besides, to improve the robustness and accuracy of the inversion, $I_{{\text{P}}\, }^i$ and $I_{{\text{S}}\, }^i$ are randomly perturbed in the range of 70%–130% of low-frequency model.

6. Substitute E for gid, that is, the global best solution in Eq. (1), and then update the position of each particle with Eqs. (1) and (2).

7. Iterate the process from step 4 to 6 to optimize the solutions until satisfying the end condition, for example, E is smaller than a specific error or maximum iteration number has reached. Now the local best solutions agree with the best solution, and the corresponding position vectors are final IP and IS.

In summary, the flow of AVA simultaneous inversion using PSO algorithm is show in Fig. 1.

 Download: larger image Figure 1. Flow chart of AVA simultaneous inversion using PSO algorithm.
3 SYNTHETIC DATA TEST

We use parts of a sonic logging curve (100 samples, sampling interval is 1 ms) as the model data. The reflection coefficients of 10°, 20° and 30° are calculated respectively according to Eq. (4), which then convolve with 50 Hz Ricker wavelet to get three angle traces. PSO algorithm is applied to estimate IP and IS according to Eq. (7). It should be noted that, we adopt fixed wavelet here, but in practical application, we need to extract angle wavelets of each angle from the seismic data and well log data.

Figure 2 shows the results of AVA simultaneous inversion with PSO algorithm after one iteration. To illustrate the inversion results in detail, we also show IS/IP in this figure. In the left three panels, blue curves are true values of IP, IS and IS/IP respectively, and the red curves are the inverted IP, IS and IS/IP respectively. Right two panels show the true synthetic records of 10°, 20° and 30° at a single-point CDP, and the synthetic data computed using the inverted IP and IS. We find that the inverted IP and IS are close to the true values after only one iteration, because we use the soft constraints mentioned above. It means good constraints can greatly improve the efficiency and stability of the inversion algorithm.

 Download: larger image Figure 2. The inversion results compared with the true values after one iteration. Left three panels show the comparison between the true values (blue curves) and the inverted results (red curves) of IP, IS and IS/IP. Right two panels show the true synthetic records of 10°, 20° and 30° at a single-point CDP, and the synthetic data computed using the inverted IP and IS.

Figure 3 shows the inversion results after 3 000 iterations. It can be seen that the inverted IP and IS are very close to the true values, and the calculated seismic records are also in good agreement with the true synthetic ones.

 Download: larger image Figure 3. The inversion results compared with the true values after 3 000 iterations. Left three panels show the comparison between the true values (blue curves) and the inverted result (red curves) of IP, IS and IS/IP. Right two panels show the true synthetic records of 10°, 20° and 30° at a single-point CDP, and the synthetic data computed using the inverted IP and IS.
4 FIELD DATA EXAMPLE

We apply our method to a 2D seismic data of Shengli Oilfield in China. The previous seismic data in this area show that the reservoirs are sand-shale interbedded, and each single sand body is very small and changes fast laterally. Therefore the conventional poststack inversion cannot recognize the sand body. P-wave velocity of sandstone is 2 600–3 600 m/s and the mudstone is 2 000–2 600 m/s according to the well log data. Prior to the AVA inversion, we have done the amplitude preserved processes such as spherical compensation, surface consistent deconvolution, random noise attenuation and so on, to make sure the final prestack amplitudes represent the reflection strength of the subsurface interfaces.

In order to reduce the cost of computation, we applied our approach to parts of three seismic common-angle gathers at 5°, 15°, and 25° (time range is 2 200–2 600 ms, 101 samples in total; CDP range is 290–510, 211 traces in total) with well gs1 located at CDP 376.

In the inversion process, we use low-frequency models of IP and IS as the initial model, which are created through the corresponding sonic logging curve smoothing and extrapolation. Figures 4 and 5 show the IP and IS profiles obtained from the AVA simultaneous inversion. The oil sand has been drilled by gs1 well at around 2 400 ms as white ellipse denoted. From the well data analysis, we find that VP of the oil sand is lower than the dry sand, but still higher than that of the surrounding mudstone, while VS of the oil sand will not be influenced with the fluid, so VS of oil sand is still high. Therefore, IP and IS have the same trend of VP and VS, which has agreement with the results in Figs. 4 and 5.

 Download: larger image Figure 4. Inverted IP section with gs1 well located at CDP 376. White ellipse denotes the reservoir.
 Download: larger image Figure 5. Inverted IS section with gs1 well located at CDP 376. White ellipse denotes the reservoir.

Since we have got IP and IS from the AVA simultaneous inversion, other elastic parameters more sensitive to the lithology and fluid, such as λρ, μρ and VP/VS can be easily achieved, because according to petrophysics knowledge, λρ=IP2–2IS2, μρ=IS2, and VP/VS=IP/IS. Figures 6 to 8 show the λρ, μρ and VP/VS profiles calculated from the inverted IP and IS profiles.

 Download: larger image Figure 6. λρ section with gs1 well located at CDP 376. White ellipse denotes the reservoir.
 Download: larger image Figure 7. μρ section with gs1 well located at CDP 376. White ellipse denotes the reservoir.
 Download: larger image Figure 8. VP/VS section with gs1 well located at CDP 376. White ellipse denotes the reservoir.

Same as Figs. 4 and 5 shown, our target reservoir is illustrated by white ellipse. The reservoir has relatively lower λρ, higher μρ and lower VP/VS, which also agree with the well log.

5 CONCLUSIONS

Prestack AVA simultaneous inversion can obtain the elastic parameters such as IP and IS at the same time, and then achieve the other elastic parameters characterizing the reservoir. Conventional prestack AVA simultaneous inversion methods often use generalized linear inversion algorithm, which has high computational efficiency, but strongly depends on the initial model and easily falls into the local optimum. In this paper, PSO algorithm is applied to prestack AVA simultaneous inversion, which does not depend on the initial model and can reach global minimum. Low-frequency impedance models of P-wave and S-wave are added into the objective function as constraints, which greatly improve the speed and stability of the inversion. After the synthetic data test, we apply this method to a 2D seismic data, and extract a variety of seismic attribute profiles, which are consistent with the logging results. Therefore, AVA simultaneous inversion using PSO algorithm has good potential in lithology identification and fluid prediction.

ACKNOWLEDGMENTS

This study was supported by the National Natural Science Foundation of China (Nos. 41004096 and 41230318). We thank Shengli Oilfield for providing the field data. Finally we thank for the anonymous reviewers for their constructive suggestion. The final publication is available at Springer via https://doi.org/10.1007/s12583-017-0809-6.

REFERENCES CITED
 Bruun, A., Schiett, R. C., Marsden, G., et al., 2005. Pre-Stack Simultaneous AVO Inversion on the South Arne Chalk Field, Danish North Sea. Ex-panded Abstracts of 67th EAGE Conference and Exhibition. Contreras, A., Torres-Verdín, C., Fasnacht, T., 2006. AVA Simultaneous Inversion of Partially Stacked Seismic Amplitude Data for the Spatial Delineation of Lithology and Fluid Units of Deepwater Hydrocarbon Reservoirs in the Central Gulf of Mexico. Geophysics, 71(4): E41-E48. DOI:10.1190/1.2212276 Contreras, A., Torres-Verdin, C., Kvien, K., et al., 2005. AVA Stochastic Inversion of Pre-stack Seismic Data and Well Logs for 3D Reservoir Modeling. Expanded Abstracts of 67th EAGE Conference and Exhibition. Fatti, J. L., Smith, G. C., Vail, P. J., et al., 1994. Detection of Gas in Sandstone Reservoirs Using AVO Analysis:A 3-D Seismic Case History Using the Geostack Technique. Geophysics, 59(9): 1362-1376. DOI:10.1190/1.1443695 Goodway, W., Chen, T., Downton, J., 1997. Improved AVO Fluid Detection and Lithology Discrimination Using Lamé Parameters; Ir, Mr and 1/m Fluid Stack from P and S Inversions:Ann. Mtg., Can. Soc. Expl. Ge-ophys.: 148-151. Gray, D., Andersen, E., 2000. Application of AVO and Inversion to For-mation Properties. World Oil, 221(7): 85-90. Hampson, D. P., Russell, B. H., Bankhead, B., 2005. Simultaneous Inver-sion of Pre-stack Seismic Data. Ann. Mtg. Abstracts, SEG.: 1633-1637. Kennedy, J., 1997. The Particle Swarm:Social Adaptation of Knowledge. IEEE International Conference on Evolutionary Computation (Indian-apolis, Indiana). IEEE Service Center, Piscataway.: 303-308. Kennedy, J., 1998. The Behavior of Particles. 7th Annual Conference on Evolutionary Programming, San Diego: 581-589. Kennedy, J., Eberhat, R. C., 1995. Particle Swarm Optimization. Proceedings of the 4th International Conference on Neural Networks. IEEE Service Center, Piscataway.: 1942-1948. Ma, X. Q., 2002. Simultaneous Inversion of Prestack Seismic Data for Rock Properties Using Simulated Annealing. Geophysics, 67(6): 1877-1885. DOI:10.1190/1.1527087 Ma, X. Q., 2001. Technical Article:Global Joint Inversion for the Estimation of Acoustic and Shear Impedances from AVO Derived P-and S-Wave Reflectivity Data. First Break, 19(10): 557-566. DOI:10.1046/j.1365-2397.2001.00211.x Mallick, S., Huang, X. R., Lauve, J., et al., 2000. Hybrid Seismic Inversion:A Reconnaissance Tool for Deepwater Exploration. The Leading Edge, 19(11): 1230-1237. DOI:10.1190/1.1438512 Mallick, S., 1995. Model-Based Inversion of Ampli-tude-Variations-with-Offset Data Using a Genetic Algorithm. Geophysics, 60(4): 939-954. DOI:10.1190/1.1443860 Misra, S., Sacchi, M. D., 2008. Global Optimization with Model-Space Preconditioning:Application to AVO Inversion. Geophysics, 73(5): R71-R82. DOI:10.1190/1.2958008 Misra, S., Sacchi, M. D., 2007. Nonlinear One Dimensional Pre-Stack Seismic Inversion with Edge Preserving Smoothing Filter. Extended Abstract, 69th EAGE conference and exhibition, London.: 26-30. Pan, G. S., Young, C. Y., Castagna, J. P., 1994. An Integrated Tar-get-Oriented Prestack Elastic Waveform Inversion:Sensitivity, Cali-bration, and Application. Geophysics, 59(9): 1392-1404. DOI:10.1190/1.1443697 Tarantola, A., 1986. A Strategy for Nonlinear Elastic Inversion of Seismic Reflection Data. Geophysics, 51(10): 1893-1903. DOI:10.1190/1.1442046