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
Nowadays, lithology and subtle reservoirs have become the main targets for oil and gas exploration. Many reservoirs have no obvious characteristics in Pwave 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 commonangle gather and the second step is to calculate Pwave velocity, Swave velocity, density or other elastic parameters using AVO approximation. In this kind of inversion, the commonangle 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 Pwave velocity, Swave 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 Pwave impedance (I_{P}) and Swave impedance (I_{S}) simultaneously. Hampson et al. (2005) proposed a simultaneous inversion by utilizing the relationship between I_{P}, I_{S} 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 modelbased 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 multiobjective 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 realcoded 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 I_{P} and I_{S} from the Pwave reflection data simultaneously based on Fatti's AVO approximation. In each iteration of inversion, I_{P} and I_{S} 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 lowfrequency impedance model are built into inversion. The estimated I_{P} and I_{S} 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 OPTIMIZATIONIn 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 p_{best}, 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 g_{best}, 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 ndimensional searching space with m particles in it. Each particle's position represents a potential solution to the objective function. That is, their position vector x_{i}=(x_{i1}, x_{i2}, ···, x_{in}) constitutes a swarm X=(x_{1}, x_{2}, ···, x_{m}). If we define the speed of each particle as v_{i}=(v_{i1}, v_{i2}, ···, v_{in}), the optimal location of the ith particle will be
$ \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. c_{1} and c_{2} are positive constants called learning factors or acceleration factors. c_{1} controls the step length for the particle to fly to its best position so far and c_{2} 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 c_{1} and c_{2} should be around 4.0, when the algorithm has better behavior. r_{1} and r_{2} 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 V_{max}, and if v_{id} > V_{max}, we set v_{id}=V_{max}; if v_{id} < V_{max}, we set v_{id}= V_{max}. 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 SCHEMEDifferent 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 I_{P}, I_{S} 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 V_{S}/V_{P} with I_{S}/I_{P}, so now the reflection coefficient only depends on I_{P}, I_{S}, density and incident angle. Here, I_{S}/I_{P} is no longer derived from the background I_{P} and I_{S}, but generated during each iteration. When I_{S}/I_{P} 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 I_{P} and I_{S} 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 I_{P} and I_{S}. The solution to be optimized contains n I_{P} and n I_{S}, 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 I_{P} and I_{S}, we can calculate the reflection coefficient for each angle according to Eq. (4), and then obtain the synthetic commonangle 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 I_{P} and I_{S} 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
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 I_{P} and n I_{S}. 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 bandlimited, lacking of lowfrequency and highfrequency 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 nonuniqueness and make the inversion stable, we apply the following constraints into the algorithm.
1. A priori I_{P} and I_{S} 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 lowfrequency model. With this constraint, the lowfrequency components will be added to the impedance solution, which are missing in the bandlimited seismic data.
2. The search boundaries are set up for each parameter. Good boundaries can reduce the nonuniqueness of the solution and minimize the computation time. The parameter boundary can be obtained from the lowfrequency 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 lowfrequency model of the sonic logging curve. That is, in the inversion process, I_{P} and I_{S} are randomly perturbed in the range of 70%–130% of lowfrequency model.
We do AVA simultaneous inversion with PSO algorithm as follows.
1. Extract at least three commonangle gathers from the forward modeling data or field data as observed seismic amplitude, that is,
2. Extract the wavelets of each commonangle gathers.
3. Initiate the particle swarm. Each particle stands for a potential solution of the model. For a nlayer model, the position vector x_{id} of each particle represents 2n potential solutions, where the former n x_{id} is
4. Calculate the reflection coefficient of each commonangle gather, and then convolve with the wavelets to obtain the synthetic records, that is,
5. Calculate the objective function E, where
6. Substitute E for g_{id}, 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 I_{P} and I_{S}.
In summary, the flow of AVA simultaneous inversion using PSO algorithm is show in Fig. 1.
Download:


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 I_{P} and I_{S} 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 I_{S}/I_{P} in this figure. In the left three panels, blue curves are true values of I_{P}, I_{S} and I_{S}/I_{P} respectively, and the red curves are the inverted I_{P}, I_{S} and I_{S}/I_{P} respectively. Right two panels show the true synthetic records of 10°, 20° and 30° at a singlepoint CDP, and the synthetic data computed using the inverted I_{P} and I_{S}. We find that the inverted I_{P} and I_{S} 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:


Figure 3 shows the inversion results after 3 000 iterations. It can be seen that the inverted I_{P} and I_{S} are very close to the true values, and the calculated seismic records are also in good agreement with the true synthetic ones.
Download:


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 sandshale interbedded, and each single sand body is very small and changes fast laterally. Therefore the conventional poststack inversion cannot recognize the sand body. Pwave 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 commonangle 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 lowfrequency models of I_{P} and I_{S} as the initial model, which are created through the corresponding sonic logging curve smoothing and extrapolation. Figures 4 and 5 show the I_{P} and I_{S} 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 V_{P} of the oil sand is lower than the dry sand, but still higher than that of the surrounding mudstone, while V_{S} of the oil sand will not be influenced with the fluid, so V_{S} of oil sand is still high. Therefore, I_{P} and I_{S} have the same trend of V_{P} and V_{S}, which has agreement with the results in Figs. 4 and 5.
Download:


Download:


Since we have got I_{P} and I_{S} from the AVA simultaneous inversion, other elastic parameters more sensitive to the lithology and fluid, such as λρ, μρ and V_{P}/V_{S} can be easily achieved, because according to petrophysics knowledge, λρ=I_{P}^{2}–2I_{S}^{2}, μρ=I_{S}^{2}, and V_{P}/V_{S}=I_{P}/I_{S}. Figures 6 to 8 show the λρ, μρ and V_{P}/V_{S} profiles calculated from the inverted I_{P} and I_{S} profiles.
Download:


Download:


Download:


Same as Figs. 4 and 5 shown, our target reservoir is illustrated by white ellipse. The reservoir has relatively lower λρ, higher μρ and lower V_{P}/V_{S}, which also agree with the well log.
5 CONCLUSIONSPrestack AVA simultaneous inversion can obtain the elastic parameters such as I_{P} and I_{S} 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. Lowfrequency impedance models of Pwave and Swave 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.
ACKNOWLEDGMENTSThis 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/s1258301708096.
REFERENCES CITED
Bruun, A., Schiett, R. C., Marsden, G., et al., 2005. PreStack Simultaneous AVO Inversion on the South Arne Chalk Field, Danish North Sea. Expanded Abstracts of 67th EAGE Conference and Exhibition. 
Contreras, A., TorresVerdí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): E41E48. DOI:10.1190/1.2212276 
Contreras, A., TorresVerdin, C., Kvien, K., et al., 2005. AVA Stochastic Inversion of Prestack 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 3D Seismic Case History Using the Geostack Technique. Geophysics, 59(9): 13621376. 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. Geophys.: 148151. 
Gray, D., Andersen, E., 2000. Application of AVO and Inversion to Formation Properties. World Oil, 221(7): 8590. 
Hampson, D. P., Russell, B. H., Bankhead, B., 2005. Simultaneous Inversion of Prestack Seismic Data. Ann. Mtg. Abstracts, SEG.: 16331637. 
Kennedy, J., 1997. The Particle Swarm:Social Adaptation of Knowledge. IEEE International Conference on Evolutionary Computation (Indianapolis, Indiana). IEEE Service Center, Piscataway.: 303308. 
Kennedy, J., 1998. The Behavior of Particles. 7th Annual Conference on Evolutionary Programming, San Diego: 581589. 
Kennedy, J., Eberhat, R. C., 1995. Particle Swarm Optimization. Proceedings of the 4th International Conference on Neural Networks. IEEE Service Center, Piscataway.: 19421948. 
Ma, X. Q., 2002. Simultaneous Inversion of Prestack Seismic Data for Rock Properties Using Simulated Annealing. Geophysics, 67(6): 18771885. 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 Pand SWave Reflectivity Data. First Break, 19(10): 557566. DOI:10.1046/j.13652397.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): 12301237. DOI:10.1190/1.1438512 
Mallick, S., 1995. ModelBased Inversion of AmplitudeVariationswithOffset Data Using a Genetic Algorithm. Geophysics, 60(4): 939954. DOI:10.1190/1.1443860 
Misra, S., Sacchi, M. D., 2008. Global Optimization with ModelSpace Preconditioning:Application to AVO Inversion. Geophysics, 73(5): R71R82. DOI:10.1190/1.2958008 
Misra, S., Sacchi, M. D., 2007. Nonlinear One Dimensional PreStack Seismic Inversion with Edge Preserving Smoothing Filter. Extended Abstract, 69th EAGE conference and exhibition, London.: 2630. 
Pan, G. S., Young, C. Y., Castagna, J. P., 1994. An Integrated TargetOriented Prestack Elastic Waveform Inversion:Sensitivity, Calibration, and Application. Geophysics, 59(9): 13921404. DOI:10.1190/1.1443697 
Tarantola, A., 1986. A Strategy for Nonlinear Elastic Inversion of Seismic Reflection Data. Geophysics, 51(10): 18931903. DOI:10.1190/1.1442046 