ARTÍCULO ORIGINAL

Genetic algorithm for estimating in-situ rock elastic constants by acoustic reflection records

Algoritmo genético para estimar in-situ constantes elásticas de rocas mediante registros de reflexión acústica

Luis Montes^1*, Ovidio Almanza², Alfredo Ghisays³

¹Departamento de Geociencias, Universidad Nacional de Colombia. Edificio Manuel Ancízar, Oficina 326. Carrera 30 N°. 45-03. Bogotá, Colombia.

²Departamento de Física, Universidad Nacional de Colombia. Edificio de Física, Matemáticas y Estadística, Oficina 337. Carrera 30 N°. 45-03. Bogotá, Colombia.

³Facultad de Ciencias Básicas, Universidad del Atlántico. Ciudadela Universitaria. Km 7 Vía Puerto Colombia. Barranquilla, Colombia.

^*Autor de correspondencia: teléfono: + 57 + 1 + 3165000 ext. 16539, fax: + 57 + 1 + 3165390, correo electrónico: lamontesv@unal.edu.co (L. Montes)

(Recibido el 4 de diciembre de 2012. Aceptado el 5 de agosto de 2013)

Abstract

The elastic properties of rocks can be estimated from the density (ρ), acoustic (Vp), and shear (Vs) wave velocities, whose values establish the reflected wave amplitudes. This paper presents an indirect method to estimate rocks' elastic properties, that uses Vp, Vs, and ρ values supplied by the inversion of acoustic reflection records. Because of the non-uniqueness and non-linear nature of the inversion, the solution must be sought in a search space by minimizing a cost function that measures the error between the observed datum and the inferred one. The search may converge to a local minimum, and then the true global minima may not be reached. Genetic algorithms have proven to be more efficient in finding the optimal solution for this type of search space.

A genetic algorithm coded in Matlab estimates ρ, Vp, and Vs through the inversion of the equation that relates them to the amplitudes and angles of incidence of the acoustic waves. Lamé elastic constant (λ), Poisson's ratio (v), and modulus of elasticity (E), compressibility (K) and rigidity (G) can be deemed from ρ, Vp and Vs.

The algorithm was tested in synthetic data to verify its robustness and stability, and then applied to seismic records showing a good performance in both cases. The presented method has a deeper scope than the refraction method, and is applicable to different engineering fields.

Keywords: Elastic constants, rocks, in situ, genetic algorithm, inversion

Resumen

Las propiedades elastomecánicas de las rocas se pueden estimar a partir de la densidad (ρ) y las velocidades de ondas acústica (Vp) y cizalla (Vs) cuyos valores establecen las amplitudes de las ondas reflejadas. Este artículo presenta un método indirecto para estimar propiedades elásticas de las rocas, que usa valores de Vp, Vs y ρ por la inversión de registros de reflexión acústica. Debido a la no-unicidad y a la naturaleza no-lineal de la inversión el motor de inferencia debe buscar una solución en un espacio de búsqueda, minimizando una función de costo que mide el error entre el dato observado y el inferido. La búsqueda puede converger en un mínimo local y no alcanzar el mínimo global verdadero. Los algoritmos genéticos han mostrado ser más eficientes en hallar la solución óptima en este tipo de espacios de búsqueda. Un algoritmo genético codificado en Matlab estima ρ, Vp y Vs a través de la inversión de una ecuación que las relaciona con las amplitudes y ángulos de incidencia de las ondas acústicas. La constante elástica de Lamé (λ), el coeficiente de Poisson (v), y los módulos de elasticidad (E), compresibilidad (K) y rigidez (G) se pueden estimar a partir de ρ, Vp y Vs. Para verificar la robustez y estabilidad del algoritmo, éste se probó con datos sintéticos y se aplicó a registros reales exhibiendo un buen desempeño en alcanzar las soluciones en ambos casos. El método presentado tiene una profundidad de sondeo mayor al método de refracción, siendo aplicable en distintos campos de la ingeniería.

Palabras clave: Constantes elásticas, rocas, in situ, algoritmo genético, inversión

Introduction

The static modulus of deformation is the parameter that best represents the mechanical behaviour of a rock, and in particular when it comes to underground excavations, therefore this modulus is a cornerstone of many geomechanical analyses. Field tests to determine this parameter directly are time consuming and expensive, and the reliability of the results of these tests is sometimes questionable [1]. Because of this, the deformation modulus is often estimated indirectly from classification systems or by geophysical methods. All methods of field modulus measurement or estimation provide values that vary from laboratory values by significant amounts; deformation modulus on intact rock samples is in the order of 5 to 20 times higher than in situ values [2]. This modulus depends on the stress conditions, being higher in areas under high stresses than in rock masses under low stresses. ''It is a steadily increasing trend in the fields of engineering geology and rock mechanics to substitute geological reality by mathematical idealizations. This lack of interest in uncertainty and field observations easily leads to reducing the quality of the input parameters'' [2]. Usually, the rock's mechanical constants are measured in rock cores or witness samples in the laboratory using different methods that provides different values [3]. Laboratory tests require high-quality core samples, not always available, particularly from thinly-bedded rock masses. Therefore, the modulus of elasticity may be estimated from other rock properties using predictor equations published in literature [4]. Attempting to diminish the associated uncertainty, Artificial Intelligence methods are being applied providing low-level errors [5].

In situ estimation thereof, has been tried by measuring the velocities of acoustic waves (P) and shear (S) in seismic reflection surveys, to which the rock density must be predetermined [6], or by using well-logs in drill holes [7]. Poisson's ratio and dynamic elastic modulus are then estimated from Vp, Vs and ρ values, and therefore, the immediate deformations [8].

The parameters ρ, Vp and Vs, are related to the amplitude of the recorded waves and the distance between the source and the receiver (offset) according to a matrix system of equations [9] which are approximated to linear equations [10]. In order to estimate ρ, Vp and Vs from CMP gathers the anterior equations that are used in an inversion process named AVO analysis (Amplitude versus offset). The seismic data are essentially sensitive to both seismic wave velocities, and the density of contrasts in the subsurface rocks. Due to the significant overlap in the elastic properties between different rock types, the mapping of the elastic properties of the rock types and porosity estimation are not trivial. The inversion itself is a problem with multiple solutions (ill-posed) which, in order to be solved, should be constrained by a priori, (information-based-model), which restricts the search space [11]. Genetic algorithms have been used extensively to solve problems inherently nonlinear [12]. The basic seismic processing corrects the influence of the topography without restriction of planar horizontal and parallel layers [13].

A genetic algorithm based on AVO analysis was coded in Matlab that invert CMP gathers, and provide Vp, Vs, and ρ values, lately used to calculate the elastic constants. Tested in synthetic seismograms, the algorithm provided reliable values of the elastic parameters. The algorithm was applied in a CMP gather acquired in the vicinity of a well with available well logs, proving elastic parameter values in rocks with small errors, compared with those calculated using the well logs.

The seismic surveys

For relatively shallow surveys, less than 20 m deep, a seismic refraction survey as that shown in figure 1A can be used to measure Vp and Vs, with the limitation that each successively deeper refractor must have a higher velocity than the shallower refractor. The waves penetrate the overburden and refract along the bedrock surface. While they are traveling along this surface, they continually refract seismic waves back to the ground surface. If deeper refractors are present and are imaged by the refraction spread, they will also refract seismic waves back to the geophones on the ground surface. As a result, a shot gather that contains both acoustic and shear waves, as shown in the record of figure 1B.

In figure 1B, the direct and refracted P and S waves associated to the shallower layer must be identified separately in order to pick their first arrivals, which might be easy in the P case because it corresponds to the first low amplitude and high frequency signal that reaches the geophone. The picking of the first arrival of the converted low frequency S wave is a hard task, because it is not possible to know a priori among the different strong amplitude recorded S waves which one corresponds to the converted wave. In another expensive survey, a shearing source generates a reflected S in the record that is commonly noisy. In general, both procedures induce errors in the estimation of P and S velocities, the error being bigger in the S case. On the other hand, the rock density ρ is measured in a cored rock sample that could be obtained depending on the drilling depth.

In a configuration as of that in figure 2B, a source and a geophone are situated symmetrically around a Common Mid Point (CMP or CDP) to record one trace, then the source-geophone distance (offset) is increased and another trace is recorded, the procedure is repeated until a set of traces is obtained, as that shown in figure 2A (CDP gather).

Elastic constants of rocks

Rocks in depth become more compact due to the litho static pressure, which increases stiffness and elastic wave velocities. In these circumstances the behaviour of the rock changes from elastic to rigid, a procedure known as enduring by strain. To characterize a rock elastically, at least two parameters among λ (Lame's constant), v (Poisson's ratio), (Elasticity modulus), (Bulk modulus) y (shear modulus) must be known. The parameters E and v are most commonly cited in engineering. The modulus of elasticity or tensile modulus (E) relates the stress component (σ_ii) with the uniaxial deformation component (ε_ii) expressed as Hooke's Law is given by equation (1):

For Sandstones E fluctuates from 0.5 to 8x105 kg/cm2 and for Shales from 1 to 3.5x105kg/cm2 [14].

In a shear stress (σik) equation (1) becomes equation (2):

The shear deformation εik depends on the shear modulus (G) that indicates the resistance of the material to be sheared. Under uniaxial stress the material shrinks in one direction (ε_ii) and expands in the other (ε_kk), being Poisson's ratio, the rate between them, or numerically, the negative ratio of transverse to axial strain. Poisson's ratio is described by equation (3)

Under the hydrostatic pressure, the volume V of a given mass of rock will be reduced to V - ΔV when the pressure is exerted uniformly all over its surface. The change in volume divided by the original volume (ΔV/V) is named volumetric strain and is linearly related with the exerted pressure according to equation (4):

K is the Bulk modulus of the rock, and measures the resistance of the material to be deformed.

The velocities with which acoustic and shear waves travel through rocks depend on the density, and on both the shear and the incompressibility modulus, according equation (5) [14]:

and equation (6)

Consequently, and can be expressed in terms of ρ, Vs and Vp [14] according equations (7-10):

Also, the Oedometric modulus that measures the variation of rigidity in depth is estimated by equation (11):

If a variable has a true value A, then the error to estimate it is given by the ratio of the difference δA between the true and the estimated value to the true value, i.e. δA/A. Suppose the variables A,B,C, . . represent independent measurable quantities used to obtain a value for some calculated quantity U (A,B,C, . . .). According the theory of errors [15, 16] if the errors for A,B,C, . . . are independent and random, the difference between the estimated and true value of U will be given by equation (12):

Considering the Equation 7, the error percentage (δG/G) in the estimation of G due to errors in the estimation of density (δρ/ρ) and shear wave velocity (δVs/Vs) will be estimated by equation (13):

Equation 13 indicates the great sensitivity of estimating G, because it depends on the square of the shear wave velocity and on the rock density, so errors of 10% on both p and Vs propagate an error of 22.3% to G.

By a similar procedure applied to equation 8, to establish the error in the estimation of Poisson's ratio due to errors on Vp and Vs it established by equation (14):

Equation (15) points out that the error committed in the estimation of Poisson s ratio depends on the errors in Vp and Vs, and in case of a Vp/Vs ratio equals 2 and errors of 10% on both Vp and in Vs estimations, Poisson's ratio will be estimated with an error of 28%.

Considering the error analysis on equation 9, the error on the estimation of E is determined by equation (15):

Equations (13-15) indicate a high sensitivity in the estimation of the elastic parameters due to errors committed in the estimation of density and velocities of shear and acoustic waves, pointing out the necessity of a reliable and robust method to estimate them.

Amplitude analysis versus offset - AVO

When a wave strikes an interface between two media, part of its energy is reflected and part is transmitted at angles that depend on the angle of incidence (θ) and velocities V_S1, V_S2, V_P1, V_P2 of each media, while their amplitude depends, in addition, on the densities ρ₁ and ρ₂.

When ρ, V_P y V _S vary slightly from one medium to another, the complex matrix equation relating the amplitude R (θ) with all of the above parameters can be approximated by the equation (16) [10]:

The terms of the above equation are defined operationally by the equations (17-23):

Genetic algorithm

The mathematical procedures to obtain information about the physical world on the basis of inference drawn from observations may provide multiple solutions, so the principal objective in a non-linear inverse problem is to locate a model for which a suitable defined cost function has a minimum. The optimal solution is searched in a large number of possibilities or state space that can be restricted by a priori information, e.g. a range of wave velocities or densities. The optimal model must explain the observations reasonably well within the limits of noise in the data. In presence of minima, a premature convergence to a solution may appear without considering the entire search space. Minima global optimization methods, which are generally stochastic algorithms, avoid these local minima and converge to a global optimal solution. Genetic algorithms establish an analogy between a set of solutions to a problem, called phenotype, and a set of individuals in a natural population, encoding the information of each solution in a string or chromosomes. Each of these individuals is seen as a point in the search space, and symbols that form the chain of chromosomes are called genes. Individuals evolve through iterations, called generations. In each generation, individuals are evaluated by measuring the adaptation function that measures similarities between the observed and synthetic data, then the best adapted models are selected to replace the previous generation. In each iteration, the selection, recombination and mutation generate new individuals until the full set of generations is achieved [17]. The structure of the genetic algorithm is shown in figure 3.

Each individual (model) is represented by the matrix, where I is the number of layers, M the length of each of the 3 binary strings for V_P, V_s and ρ, whose values fluctuate within pre-established ranges. An initial population is generated randomly using a uniform distribution. Amplitudes are calculated for each angle and each reflector with equation (12), to create the synthetic seismogram by convolving with the wavelet extracted from the CDP record using statistical methods. The RMS error between the synthetic seismogram, associated with each individual, and the CDP gather is calculated. By selection, combination, and mutation, new patterns are generated, iterating until the number of established generations pass and the optimal model is delivered.

Inverting synthetic seismograms

The figure 4A shows a synthetic CDP gather generated through equation 16 that simulates the seismic response in the 7- layer model depicted in figure 4B. The genetic algorithm was applied to this CDP gather doing the inversion, and providing the solution model whose numerical values are included in figure 4C. The comparison between the former and the estimated model resulted in V_P, V_S and ρ maximum errors of 0.269 km/s, 0.176 km/s and 0.089 g/cc, which are considered acceptable small errors.

Equations (7-9) together with the V_P, V_s and ρ values of both former and predicted models, were used to calculate Poisson's ratios, elasticity modules and estimation errors. The results shown in table 1 indicate errors below 20% in the estimation of these parameters. This order error is considered acceptable to design calculations in engineering that involve rocks [18].

Inverting real data

The figure 5 shows the sonic log measured in μsec/m, and the density log in g/cc. On the right, it shows a seismic image of the subsurface in the vicinity of the well, along 40m in the 1400­-1800 ms interval. The seismic in the vicinity of the well is tied by identifying lithic units, and the time-depth relationship. Then, a CDP gather located in the vicinity of the well was selected and processed by filtering and applying static corrections a plane datum, and finally, applying a dynamic correction (Normal Move out) until the section in figure 6A was achieved. For a full explanation of the anterior steps a specialized literature is recommended [13].

Figure 6A shows the 1550-1750 ms interval of a gather previously processed until the reflectors become horizontal, where the estimated incidence angles vary between 0° and 15°. As a result of applying the genetic algorithm from 1600 to 1700 ms of the CDP gather, Poisson ratio and elasticity coefficient values were obtained. The true values were calculated using well log data and equations 6, 8 and 9. Figure 6B shows on the left the predicted and the true Poisson ratio curves and on the right the true and the predicted elasticity coefficient curves, where true and predicted values look very close. The errors in the estimation of these parameters are below 20% consistent with the observed errors in the case of the synthetic model.

Conclusions

An indirect method to estimate rock's elastic properties in situ, by the inversion of reflection records, is presented. Implemented in a Matlab code, the GA was tested on synthetic data supplying reliable Poisson's ratios and elasticity coefficients. Applied to a depth interval of a CDP gather located near a well, the GA provided reliable Poisson's ratios and elasticity coefficients according to available well logs. Theoretical analysis of errors pointed out the high sensibility of elastic parameters associated to the uncertainty in shear wave velocity estimation. The method overcomes the limitations faced of the refraction method, the high uncertainty in the estimation of shear wave velocity and availability of a core rock to measure density in the laboratory. The results, achieved in the synthetic and real cases, outline the robustness and low sensibility of the GA. The level of uncertainty in the parameter estimation is considered acceptable in calculations that involve rocks. This non-expensive method has a deeper scope than the refraction method and is applicable to different engineering fields.

Acknowledgements

The authors express their gratitude to the Universidad Nacional de Colombia and in particular the Department of Geosciences, for their support during the development of the project.

References

1. E. Hoek, M. Diederichs. ''Empirical estimation of rock mass modulus''. International Journal of Rock Mechanics & Mining Sciences. Vol. 43. 2006. pp. 203-215.         [ Links ]

2. P. Palmstrom, S. Rajbal. ''The Deformation Modulus of Rock Masses, - Comparisons between in situ tests and indirect estimates''. Tunnelling and Underground Space Technology. Vol. 16. 2001. pp. 115-131.         [ Links ]

3. P. Santi, J. Holschen, R. Stephenson. ''Improving Elastic Modulus measurements for based on Geology''. Enviromental & Engineering Geosciences. Vol. VI. 2000. pp. 333-346.         [ Links ]

4. I. Ocak. Empirical estimation of intact rock elastic modulus. 21^th international mining congress of Turkey. Antalya, Turkey. 2009. pp. 165-172.         [ Links ]

5. I. Ocak, S. Seker. ''Estimation of Elastic Modulus of Intact Rocks by Artificial Neural Network''. Rock Mech. Rock Eng. Vol. 45. 2012. pp. 1047-1054.         [ Links ]

6. L. Briceño, R. Díaz. ''Medida de propiedades dinámicas de rocas ''in situ'' y cómputo de parámetros elastomecánicos''. Geofísica Colombiana. Vol. 3. 1995. pp. 73-79.         [ Links ]

7. T. Odumosu, C. Torres, M. Salazar, J. Ma, B. Voss, G. Wang. Estimation of Dry-Rock Elastic Moduli Based on the Simulation of Mud-Filtrate Invasion Effects on Borehole Acoustic Logs. 2007 SPE Annual Technical Conference and Exhibition. Anaheim, Ca, USA. 2007. pp. 1-16.         [ Links ]

8. A. Manilla, P. Garnica, A. Pérez. ''Evaluación indirecta de los módulos elásticos de rigidez in situ y la relación entre Vp/Vs y el ángulo de fricción interna''. Publicación Técnica. N° 225. 2003. pp. 5-25.         [ Links ]

9. K. Aki, P. Richards. Quantitative seismology. 2^nd ed. Ed. University Science Books. Sausalito, Ca, USA. 2002. pp. 119-183.         [ Links ]

10. J. Fatti, G. Smith, P. Vail, P. Strauss, P. Levitt. ''Detection of gas in sandstone reservoirs using AVO analysis: A 3-D seismic case history using Geostack technique''. Geophysics. Vol. 59. 1994. pp. 1362-1376.         [ Links ]

11. M. Sen. Seismic Inversion. 1^st ed. Edited by the Society of Petroleum Engineers. Richardson, Tx, USA. 2006. pp. 61-67.         [ Links ]

12. K. Castaño, G. Ojeda, L. Montes. ''Thin-layer detection using spectral inversion and a genetic algorithm''. Earth Sci. Res. J. Vol. 15. 2011. pp. 121-128.         [ Links ]

13. O. Yilmaz. Seismic Data Analysis: processing, inversion and Interpretation of Seismic Data. 2^nd ed. Published by the Society of Exploration Geophysicists. Tulsa, Ok, USA. 2001. pp. 25-124.         [ Links ]

14. G. Mavko, T. Mukerji, J. Dvorkin. The Rock Physics Handbook. 1^st ed. Ed. Cambridge University Press. Cambridge, UK. 2009. pp. 21-76.         [ Links ]

15. D. Cacuci. Sensitivity and uncertainty analysis: Theory. Vol. 1. Ed. Chapman & Hall/CRC. Boca Raton, USA. 2003. pp 132-141.         [ Links ]

16. A. Saltelli, S. Tarantola, F. Campolongo, M. Ratto. Sensitivity analysis in practice: a guide to assessing scientific models. 1^st ed. Ed. John Wiley & Sons. Chichester, England. 2004. pp. 109-148.         [ Links ]

17. M. Mitchell. An Introduction to Genetic Algorithms. 1^st ed. Ed. MIT Press. Boston, Mass., USA. 1998. pp. 35-81.         [ Links ]

18. I. Farmer. Engineering Property of Rocks. 1^st ed. Ed. E & F. N. Spon Ltd. London, England. 1968. pp. 59-85.         [ Links ]