Inversión de amplitudes sísmicas a offset variable mediante técnicas de Monte Carlo

LUIS CARA Y MIGUEL BOSCH

Universidad Central de Venezuela, Facultad de Ingeniería, Laboratorio de Simulación e Inversión Geofísica. email: caraluis@gmail.com; miguel.bosch@ucv.ve

RESUMEN

Se presenta una metodología de inversión a partir de grupos de trazas sísmicas a apertura fuente-receptor variable para estimar las propiedades elásticas del medio: lentitud P, lentitud S y densidad, utilizando la técnica de Monte Carlo. Los datos observados son las amplitudes de ondas sísmicas de reflexión en el dominio del tiempo como función de la separación emisor-receptor. Para obtener sismogramas sintéticos multicomponentes a apertura variable utilizamos el método de trazado de rayos y resolvemos las ecuaciones de Zoeppritz de forma exacta. También se utiliza información estadística previa sobre los parámetros elásticos del medio, incluyendo sus correlaciones. Con este método se produce una descripción completa de las incertidumbres de los parámetros del modelo.

Palabras clave: Inversión elástica, Datos sísmicos, Apertura variable, Estadística, Velocidad sísmica, Densidad.

Inversion of variable offset seismic amplitudes using Monte Carlo techniques

ABSTRACT

We describe a method for the estimation of the elastic parameters (slowness P, slowness S and mass density) in a stratified medium by inverting variable offset seismic amplitudes with Monte Carlo techniques. The observed data are the seismic amplitudes in time for different source-receiver offsets. For the simulation of the seismic data we convolve a source waveletwith the reflectivity series obtained with ray tracing and the exact solution of the Zoeppritz equations; we also include in the inversion scheme prior information of the statistics of the elastic medium parameters. We invert the seismic data using a sampling Markov Chain Monte Carlo technique that generates a long chain of the stratified model realizations, all probable in terms of the geophysical and prior statistical information, and convergent to the posterior probability density. In this technique we combine the Gibbs and the Metropolis algorithms. We validate the method with numerical tests, using elastic models simulated on the basis of the statistics of real well-logged data, including mean property values, standard deviations and correlation across the elastic properties. With this method we produce a complete description of the uncertainties of the model parameters.

Keywords: Elastic inversion, Seismic data, Variable offset, Statistics, Seismic velocity, Density.

Recibido: febrero de 2007 Revisado: noviembre de 2007

INTRODUCCIÓN

Los métodos sísmicos son una de las fuentes principales de información en el área de la geofísica. Si bien el tratamiento de los datos para crear imágenes de reflectividades ha sidola técnica más utilizada, existe creciente interés por técnicasque producen como resultado una descripción directa delas propiedades físicas del medio terrestre como los parámetros elásticos, a través de métodos de inversión.

La inferencia de las propiedades del medio, a partir de datos observados, es lo que se conoce como un proceso de inversión (Tarantola, 1987), para el cual se deben conocer primeramente el conjunto de relaciones y leyes físicas quepredicen las observaciones. Esto es lo que se conoce comomodelado directo o problema directo. En el caso de sísmicael modelado directo es implementado a través de un algoritmoque produce un sismograma sintético (Treitel y Lines, 2001), este sismograma depende del cambio en las propiedades elásticas de las rocas y del ángulo de incidencia con que las ondas penetran el subsuelo. Las ecuaciones para los coeficientes de reflectividad y transmisividad, los cuales dependen del ángulo de incidencia, están contenidas en las ecuaciones de Zoeppritz y la solución descrita por Aki y Richards (1980), y son las que permiten construir la serie de reflectividades para así convolucionarla con una ondícula fuente adecuada.

El modelado inverso es realizado en este trabajo siguiendo el enfoque de muestreo o de Monte Carlo, utilizado porMosegaard y Tarantola (1995) y Bosch (1999). En estemétodo se produce un gran número de configuraciones del modelo a partir de las cuales es posible calcular laprobabilidad de cualquier rango de un parámetro o un conjunto de parámetros del modelo del subsuelo. En este trabajo los parámetros del modelo son las velocidades P y S, así como la densidad. La técnica está orientada a obteneruna descripción completa de la incertidumbre sobre los parámetros del modelo.

Técnicas de inversión elástica que persiguen el mismo objetivo de estimar los parámetros elásticos del medio a partir de sísmica a offset variable, se pueden encontrar en trabajos de Torres-Verdín et al. 1999 y Roy et al. 2004. Enestos trabajos se emplea un enfoque de optimización, en lugar del enfoque de muestreo y cálculo de probabilidadesque se usaron en este trabajo. En el trabajo de Bosch (2003) se describe una técnica similar a la aquí desarrollada, aplicada a la inversión de sísmica post-apilamiento.

METODOLOGÍA

La formulación de este problema de inversión se realiza en el marco de los métodos de inferencia probabilista. Esencialmente se describe el medio material que se exploracon la propagación de las ondas sísmicas, mediante unmodelo paramétrico constituido por los parámetros elásticos(velocidades P, S y densidad) para cada una de las capas. Alfijar los valores de estos tres parámetros en la secuencia decapas se fija la configuración particular del modelo delsubsuelo. El conjunto de todas las posibles configuraciones (usualmente infinitas) conforman el espacio de parámetros del modelo.

El conocimiento del subsuelo se describe entonces definiendo probabilidades en el espacio de parámetros, através de la definición de funciones de densidad deprobabilidad. Es necesario hacer la distinción usual en losproblemas de inferencia, entre el estado de conocimientodel medio antes y después de la realización y procesamientode las observaciones físicas. El conocimiento anterior sedescribe mediante una densidad de probabilidad previadefinida sobre el espacio de parámetros. La realización delas observaciones y sus implicaciones sobre la estructura ypropiedades del medio, según el conocimiento de lasrelaciones que existen entre ellas (leyes físicas), da lugar auna actualización de ese conocimiento, el cual se representamediante una densidad de probabilidad posterior sobre el espacio de parámetros. La densidad de probabilidad posterior, definida para este trabajo, es el resultado decombinar la información previa y los datos geofísicos, y se calcula a partir de la siguiente fórmula:

donde:

L(m) es la función de verosimilitud, la cual contiene la información relativa a la geofísica y se define usualmente mediante el desajuste entre los datos geofísicos observadosy los calculados a partir del modelo, m es el arreglo deparámetros del modelo; ρ(m) es la densidad previa, querepresenta como ya se dijo todo el conocimiento que setiene del medio antes de las observaciones; c es una constante de normalización y es la densidad de probabilidad posterior, que combina la información que se tiene del medio antes y después de las observaciones. Laexpresión anterior es común en inferencia estadística, para una mayor fundamentación de esta expresión (Tarantola, 1987).

La técnica diseñada consta de tres etapas de desarrollo, la primera de ellas es el cálculo del problema directo con la finalidad de obtener las observaciones geofísicas. La segunda etapa está vinculada a la simulación estadística de los parámetros elásticos para así obtener configuraciones de estos parámetros que sean acordes con los genuinosrangos y estadísticas del medio. Para este objetivo se utilizaron registros de pozos de una zona al oriente del país correspondiente a la Formación Oficina. La tercera etapa es la resolución del problema inverso mediante el método de Monte Carlo.

El problema directo de este trabajo se elaboró con el objetivo de obtener los sismogramas sintéticos multicomponentespre-stack, para ello se utilizó el método de trazados de rayos en un modelo de capas horizontales de igual tiempo con velocidades P, S y densidad conocidas. La trayectoria detodos los rayos es calculada a través de la ley de Snell; la amplitud de las ondas con sus respectivos cambios demodos es calculada por las ecuaciones de Zoeppritz,resueltas por Aki y Richards (1980); y los tiempos de viajes son obtenidos de la ecuación de Dix (1955).

La fase de simulación se inicia con la creación del modelo previo, el cual es el resultado de la estadística de los registros de pozos aplicada sobre los parámetros elásticos (velocidades P, S y densidad). Los valores consideradosson las medias aritméticas, desviaciones estándar, varianzasy correlaciones de cada parámetro. A continuación, seseleccionan capas al azar y se varían los tres parámetroselásticos a través del método de la raíz cuadrada de la matrizde covarianza. Esta es la matriz que vincula estadísticamente los tres parámetros elásticos por cada capa seleccionada,considerando para ello un modelo Gausiano multivariado que permite generar cambios con desviación Gausiana paracada parámetro alrededor de su media aritmética. Siguiendoeste enfoque se garantiza que se cumple la estadística de los pozos a lo largo de todo su registro. Esta es la información previa que se conoce del medio.

En la figura 1 se muestran los gráficos cruzados (1a, 1b, y 1c) que se obtienen para los tres parámetros considerados siguiendo la estadística previa para un conjunto de realizaciones del modelo tomada al azar, y en los gráficos 1d, 1e y 1f se muestra el histograma de frecuencias para cadaparámetro. En las simulaciones Gausianas correlacionadasutilizadas en la figura 1, se utilizó una correlación positivaentre las lentitudes P y S, y una correlación negativa de éstas respecto a la densidad. Estas correlaciones, lascorrespondientes desviaciones estandar y valores mediosse obtuvieron a partir del análisis de los registros de pozo.

Figura 1. Gráficos cruzados e histogramas de las propiedades elásticas de los estratos simulados estadísticamente. (a) Gráfico cruzado de lentitud P contra lentitud S. (b) Gráfico cruzado de lentitud P contra densidad. (c) Gráfico cruzado de lentitud S contra densidad. (d) Gráfico de porcentaje de frecuencia para la lentitud P. (e) Gráfico de porcentaje de frecuencia para la lentitud S. (f) Gráfico de porcentaje de frecuencia para la densidad.

Para completar el modelo estadístico de la densidad de probabilidad posterior, hace falta aún definir la función deverosimilitud geofísica, basada en el ajuste entre los datoscalculados a partir del modelo y los datos observados. En este trabajo se define como:

donde:

d^cal y d^obs son las amplitudes sísmicas calculadas y observadas, respectivamente, para todas las muestras de tiempo y aperturas (offsets). Los índices i se refieren a las muestras de tiempo, y j a las aperturas emisor-receptor, σ²_d es la varianza de los datos.

Finalmente para resolver el problema de inferencia, siguiendo el enfoque de muestreo o de Monte Carlo, se requiere generar una colección grande de configuraciones del modelo. Con este objetivo, se dispone de un gran arsenal de técnicas y algoritmos estudiados y desarrollados en el campode la estadística. La mayor parte de estas técnicas estánsoportadas en métodos de cadenas, como las cadenas deMarkov, que consisten en generar una serie de realizaciones del modelo a partir de la modificación de la anterior, y respetando ciertas reglas estadísticas en cada transición o paso de la cadena. Al repetirse un gran número de veces elconjunto de configuraciones visitadas por la cadena, éstaconverge a una muestra estadística de la densidad de probabilidad objetivo. Una de estas técnicas es el algoritmo de Metrópolis, el cual permite producir una configuracióncandidata generada a partir de una distribución de muestreofuente, que no es la distribución objetivo. La distribución de muestreo fuente puede ser cualquiera, pero es deseablepara eficiencia del algoritmo que sea cercana de alguna manera a la distribución objetivo que se desea muestrear.

El algoritmo se basa en comparar la configuración candidata y la configuración corriente, para decidir si la candidata se acepta como el siguiente paso de la cadena o si se rechaza repitiéndose la configuración corriente como nuevo eslabón.

En síntesis, para generar nuestra cadena de Markov se:

- Escogen al azar varias capas a perturbar.

- Se generan los tres parámetros elásticos de esas capas según una distribución Gaussiana multivariada.

- Se evalúa la verosimilitud de las amplitudes sísmicas.

- Se aplica el algoritmo de Metrópolis para aceptar o rechazar la modificación.

Finalmente, de la colección de realizaciones obtenidas se calculan las estadísticas finales: valores medios ehistogramas de probabilidad. Se eliminan para estos cálculosuna primera fase de la cadena que está influenciada por la configuración inicial.

EJEMPLO CON DATOS SINTÉTICOS

A continuación validamos la técnica de Monte Carlo desarrollada para la estimación de parámetros elásticos, con la realización de un grupo de pruebas que denominamospruebas sintéticas. Estas pruebas no provienen de datos de campo sino de datos simulados en el computador. Para ellose produce una simulación geoestadística a partir de registros de pozos, discutida anteriormente, al cual se le calculó la correspondiente respuesta sísmica en reflexión a apertura variable, siguiendo para ello la metodología delproblema directo. El conjunto de parámetros elásticos obtenidos en una realización tomada al azar se considera como la configuración «verdadera» y los datos sísmicos calculados a partir de ella como los datos «observados». Utilizando estas observaciones sintéticas se procede a aplicar la metodología propuesta para inferir el perfil de parámetros elásticos verdaderos.

En la figura 2 se muestran 5 realizaciones que se seleccionaron al azar de un conjunto de dos millones que se obtuvieron en la fase de inversión, estas realizaciones son consistentes con el modelo estadístico considerado como previo y además todos estos modelos ajustan los datos de sísmica considerados como observados.

Figura 2. Muestra de (líneas en gris) cinco modelos generados por la cadena de Markov que toma muestras según la densidad de probabilidad posterior. En negro el modelo “verdadero”. (a) Lentitud P. (b) Lentitud S. (c) Densidad.

En la figura 3 se muestra el ajuste que hay entre los datos calculados y observados, es decir, correspondientes al modelo verdadero, para uno de los modelos mostrados anteriormente en la figura 2. Podemos hacer notar el buenajuste entre los datos observados y los calculados. Para lasprimeras capas el error en el ajuste de los datos es mayor, debido a los altos ángulos de incidencia que se están considerando en ese caso, los cuales producen altas amplitudes.

Figura 3. Sismogramas calculados, observados y residuos separados en sus componentes vertical y horizontal para una de las configuraciones mostradas en la figura 2. (a) Sismograma calculado vertical. (b) Sismograma observado vertical. (c) Diferencia entre los sismogramas: calculado vertical y observado vertical. (d) Sismograma calculado horizontal. (e) Sismograma observado horizontal. (f) Diferencia entre los sismogramas: calculado horizontal y observado horizontal.

La figura 4 muestra la evolución del ajuste de los datos calculados y observados pero a través de todo el númerode iteraciones que se consideraron en la fase de inversión. Para medir el ajuste de los datos se utilizó la estadística Chicuadrado.Para el cálculo de esta estadística se divide lasumatoria de la diferencia cuadrática entre los datoscalculados y observados (comúnmente se le llama residuocuadrático) entre la varianza de los datos y luego entre el número de observaciones.

Figura 4. Gráficos del progreso en la disminución del logaritmo en base 10 de los residuos de los datos, medidos mediante la estadística Chi-cuadrado. (a) Componente vertical. (b) Componente horizontal.

En general los modelos obtenidos a partir de la cadena como realizaciones de la densidad de probabilidad combinada, son ejemplos de modelos que cumplen con la información geofísica y geoestadística que se han dado. El conjunto formado por dos millones de realizaciones del modelo es una muestra significativa de la densidad de probabilidad, y los aspectos variables entre estas realizaciones están indicando la incertidumbre en la estimación resultante de la inversión de amplitudes sísmicas. En la figura 5 se muestra la probabilidad acumulada para las tres propiedades elásticas (lentitud-P, lentitud-S y densidad) en función del tiempo de reflexión, calculada a partir del conjunto de realizaciones generadas. Los mismos describen de manera completa la incertidumbre en la estimación de estos parámetros.

Figura 5. Histograma de frecuencias de los 2 millones de realizaciones calculadas con la cadena convergente a la densidad de probabilidad posterior. La línea gris clara representa el modelo “verdadero”, la línea negra el promedio calculado a partir de la muestra. (a) Lentitud P. (b) Lentitud S. (c) Densidad.

Estos gráficos muestran que para rangos reales de propiedades elásticas y aperturas emisor-receptor, la señal sísmica contiene información sobre la lentitud P, lentitud S y densidad. De esta manera es posible estimar estos parámetros elásticos a partir de la sísmica mediante las técnicas de inversión. Las pruebas efectuadas, en particular la correspondiente a los datos de parámetros elásticos de la Formación Oficina, mostraron que la certidumbre en la estimación de estos parámetros es mayor para la lentitud P,seguida de la lentitud S y la densidad, esta última reflejando mayor incertidumbre.

Como se observa en la figura 3, utilizamos para estas pruebas un rango de aperturas entre 0 y 4000 m para una zona deestudio que su ubica entre 0.6 y 2 s de tiempo de viaje en iday vuelta para reflexiones verticales. Este rango en tiempocorresponde aproximadamente a una profundidad entre 1300 y 5300 m para la región modelada, según las velocidadessísmicas consideradas. Los ángulos de incidencia varían enrango según la profundidad, encontrándoseaproximadamente entre 0° y 56° para la parte superior de laregión modelada (a 1300 m de profundidad) y entreaproximadamente 0° y 20 para las capas más profundas (a5300 m de profundidad). Como es de esperar, en la medida en que se incrementa la profundidad disminuye el rango de ángulos de incidencia; tanto profundidad como ángulos de incidencia influyen en la estimación. Vemos por lo tanto en la figura 5, que para las tres propiedades la incertidumbre se incrementa con la profundidad, con mayor acento para la densidad.

DISCUSIÓN Y CONCLUSIONES

En este trabajo se ha desarrollado una metodología de inversión de amplitudes sísmicas para apertura emisorreceptorvariable, mediante la cual se genera un númerogrande de configuraciones de modelos del subsuelo, queresultaron probables en términos de su ajuste con lasobservaciones (de datos sintéticos) y con la estadística quese obtuvo de pozos localizados en una zona al oriente delpaís. Dicha metodología está resumida en lo que se conocecomo técnicas de muestreo, o Monte Carlo, utilizandocadenas de Markov. Éstas garantizan la convergencia a una muestra de la densidad de probabilidad posterior.

La metodología empleada para el cálculo del problema directo es suficientemente eficiente para aplicar el método de Monte Carlo, ya que permite calcular los sismogramas multicomponentes en tiempos razonables.

La técnica se validó en sus diversas componentes: lasimulación de los tiempos de viaje y coeficientes de reflexión, la simulación estadística de las propiedades elásticas, y las pruebas sintéticas de inversión. Uno de los aspectos adestacar de esta técnica de inversión, que utiliza métodos de Monte Carlo, es que la misma permite una descripción completa de la incertidumbre en los parámetros elásticos estimados, como resultado de la combinación de los errores de los datos y la relación no lineal entre los datos y los parámetros elásticos del medio.

La investigación sigue en proceso y prevemos hacer uso dedatos de sísmica pre-stack suministrados por Petrobras para combinarlos con la información geoestadística caracterizada a partir de registros de pozos, para así obtener una inversión que no sólo contenga a los parámetros elásticos sino también descripción del reservorio en términos de la porosidad, litología y fluidos.

AGRADECIMIENTOS

Los autores agradecen a PETROBRAS, y en particular a Alonso Navarro, Manuel Díaz y Luz Mery Rodríguez, porhaber suministrado valiosa información como los registros de pozos utilizados para realizar los ejemplos con datos sintéticos de este proyecto, así como datos de sísmica quese utilizaran en la segunda fase del mismo. Se agradece elapoyo para esta investigación del CDCH-UCV a través de los proyectos: PG-08-00-5631-2004 y PI-08-00-5627-2004.

REFERENCIAS

1. AKI, K., AND RICHARDS, P.G. (1980). Quantitative seismology: Theory and methods.W.H. Freeman and Co, New York.        [ Links ]

2. BOSCH, M. (1999). Lithologic tomography: From plural geophysical data to lithology estimation. Journal ofGeophysical Research, 104, 749-766.        [ Links ]

3. BOSCH, M. (2003). Inferencia Estadística de Porosidad e Impedancia a partir de ondas Sísmicas Vía Métodos de Monte Carlo. Trabajo de Ascenso, Caracas, Universidad Central de Venezuela, 68 p.        [ Links ]

4. DIX, C. H. (1955). Seismic velocities from surface measurements, Geophysics, 20, 68-86.        [ Links ]

5. ROY, I. G., SEN, M. K., TORRES-VERDÍN, C., AND VARELA, O. J. (2004). Prestack inversion of a Gulf of Thailand oceanbottom cable (OBC) data set. Geophysics, 69, 1470-1477.        [ Links ]

6. MOSEGAARD, K. AND TARANTOLA, A. (1995). Monte Carlo Sampling of Solutions to Inverse Problems: Journal of Geophysical Research. 100, B7, 12.431-12.447.        [ Links ]

7. TARANTOLA, A. (1987). Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation. New York: Elsevier. 613 p.        [ Links ]

8. TREITEL, S., AND L. LINES, L. (2001). Past, present and future of geophysical inversion- A new millennium analysis. Geophysics, 66, 1, 21-24.        [ Links ]

9. TORRES-VERDÍN, C., VICTORIA, M., MERLETTI, G., AND PENDREL, J. (1999). Trace-based and geostatistical inversion of 3- D seismic data for thin-sand delineation: An application in San Jorge Basin, Argentina. The Leading edge,18, 1070- 1077.        [ Links ]