JUAN RODRIGUES ¹, MIGUEL BOSCH ¹, LUIS CARA ¹, ALONSO NAVARRO ² Y MANUEL DÍAZ ²

¹ Universidad Central de Venezuela, Facultad de Ingeniería, Laboratorio de Simulación e Inversión Geofísica. email: jmiguelrt@cantv.net; miguel.bosch@ucv.ve; luiscara1@yahoo.es

² Petrobras Energía Venezuela, email: alonso.navarro@petrobras.com; manuel.diaz@petrobras.com

En este trabajo presentamos un método de inversión de amplitudes sísmicas para estimar, de manera conjunta, valores de impedancia acústica y de porosidad total. Para ello empleamos información de fuentes variadas, como son la geología, la petrofísica y la geofísica. Utilizamos datos que responden al medio de propagación en escalas distintas, por lo que fue necesario desarrollar un modelo del medio a dos escalas, relacionadas entre sí mediante la fórmula de Backus para velocidades y los correspondientes promedios aritméticos para la densidad y la porosidad. El objetivo principal del trabajo es la combinación de la información de los registros de pozos, de porosidad e impedancia, con los datos de amplitudes sísmicas. Mediante esta combinación podemos mejorar la resolución vertical de las propiedades del reservorio en los sectores cercanos a los pozos, representando capas de hasta 1ms de espesor, en comparación con la resolución que obtendríamos mediante una inversión sísmica convencional que está en el orden de los 6ms de espesor en esta área. El volumen de propiedades estimado con la técnica desarrollada explica simultáneamente los datos sísmicos post-apilamiento y los registros de pozos adquiridos en el área.

Palabras clave: Inversión sísmica, Geoestadística, Optimización, Impedancia acústica, Porosidad.

In this work we present a method to jointly estimate acoustic impedance and total porosity by inverting seismic amplitudes and matching well-logged information. We combine geophysical, petrophysical and geological information in a two-scale model. The properties in the model are related across scales with the appropriate averages, i.e., the Backus average for compressional velocities and arithmetic averages for porosity and density. A major objective of the work is to combine the well-logged information of porosity and impedance with seismic information. With this combination we can improve the vertical resolution of the model in the regions of the volume influenced by the well data according to a spatial statistical model. In an example case we estimated the medium properties at 1ms interval, while the vertical resolution of a conventional seismic inversion is of the order of 6ms. The property volumes of porosity and impedance estimated with our technique jointly explain the seismic data and the well-logged data acquired in the area.

Recibido: marzo de 2007 Revisado: noviembre de 2007

INTRODUCCIÓN

El objetivo básico de la industria petrolera es hallar más y mejores yacimientos de petróleo; gracias a esto han nacido diversas técnicas de exploración, que toman enfoquesvariados sobre el tema. La técnica más importante en la búsqueda de petróleo es la sísmica. El fenómeno que estudiala sísmica es la respuesta del medio al paso de una onda sísmica, que al ser analizada permite inferir un modelo del subsuelo.

Para poder plantear un modelo mediante información sísmicadebe resolverse antes el problema inverso, que significaproponer un modelo del subsuelo que cumpla las teoríasfísicas de la propagación sísmica y los datos recolectadosen el levantamiento sísmico. Las técnicas para resolver el una técnica de optimización, que consiste en mejorar, de forma iterativa, una configuración de parámetros para maximizar la función de densidad de probabilidad.

Las propiedades que se representan en el modelo del medio son la impedancia y la porosidad logarítmica, esta última una transformada de la porosidad total. Las amplitudes sísmicas después de apilamiento tienen relación directa con la impedancia acústica, y esta última guarda relación con la porosidad total de las rocas a través de relacionespetrofísicas, como la transformada de Wyllie, entre otras. La porosidad logarítmica, por su parte, es más adecuada para los cálculos estadísticos. Más detalles sobre la selecciónde estos parámetros para problemas de inversión de amplitudes sísmicas se describen en el trabajo de Bosch (2003).

La técnica de optimización que aquí se aplicó, es una técnica novedosa que realiza la inversión en conjunto de la impedancia y la porosidad logarítmica. Además de la inversión conjunta también integra en el modelo información de otras fuentes, como la petrofísica (registros de pozo) y la geológica. Para combinar la información de registros de pozos con la de la sísmica, implementamos un modelo a doble escala relacionando ambas con los promedios apropiados para cada propiedad.

Así pues, se obtuvo un modelo de impedancia y otro de porosidad logarítmica, que cumple con los datos sísmicos cuando se está alejado del pozo y con los registros de pozo cuando se está cerca de ellos.

En la realización de este trabajo se siguieron varios pasos que a continuación se explican.

Para poder lograr relacionar la impedancia con la porosidad empleamos la transformada de Wyllie. El modelo de Wylliepara velocidades sísmicas, considera a la roca como unacomposición de su matriz y el fluido que ocupa los poros (Hilterman, 2001). Bajo esta concepción la fórmula resulta:

donde:

V_mat es la velocidad de onda P de la matriz, V_flu es la velocidad de onda P del fluido y φ y V son la porosidad y velocidad de onda P resultante de la roca respectivamente. Para la densidad se utiliza la fórmula:

donde:

ρ es la densidad de la roca, ρ_mat la densidad de la matriz, ρ_flu la densidad del fluido y φ la porosidad de la roca. Realizando el producto de la velocidad por la densidad se obtiene la fórmula para la impedancia acústica como:

donde:

La porosidad logarítmica es una propiedad de distribución normal que está definida en todos los números reales (Bosch, 2003) y su fórmula inversa es:

Sustituyendo la fórmula (5) en la (3), se obtiene la transformada de Wyllie para la porosidad logarítmica:

La resolución de la sísmica depende de la longitud de onda de la señal, que es mayor que la de los registros de pozo. Por tanto el primer paso consistió en combinar estas escalas.Para este fin utilizamos la fórmula de Backus (1962) develocidades equivalentes. El parámetro que modelamos ados escalas es la impedancia acústica, por lo que la fórmulade Backus para velocidades fue ajustada para impedancia ypara su descripción en tiempo en lugar de profundidad. La fórmula para el cambio de escala en la impedancia acústica es:

donde:

El modelo estadístico que se empleó en este proyecto fue un modelo Gausiano multivariado que se puede escribir como:

donde:

p es el arreglo de parámetros (la serie de impedancias y porosidades en alta resolución), es el valor esperado de los parámetros de acuerdo a la información previa disponible en el área, c es la constante de normalización y C^-1_M es la matriz de covarianza del modelo.

El llenado de la matriz de covarianza se logra empleando la función de covarianza que se extrae de los datos. Aquí se muestra la función de covarianza para un modelo Gausiano que es uno de los más utilizados:

donde:

h es la distancia entre los dos puntos del volumen, σ 2 es la varianza y a es el rango de covarianza, que aproximadamente corresponde a la distancia hasta donde los datos estáncorrelacionados de manera significativa. El valor de a sedefine a partir de la geología de la zona.

, en algunos puntos del

, es posible proponer un

donde:

es la media previa en el punto  x₁ , es la media de los datos conocidos (en este proyecto se considera que son iguales, porque los valores conocidos provienen de pozos y los puntos donde se desconocen los valores están en la misma área), C₁₂ es el vector de covarianza entre el punto desconocido y los puntos conocidos y  C⁻¹₂₂ es la matriz inversa de la matriz de covarianza de los puntos conocidos (Deutsch y Journel, 1998).

Aplicando la fórmula (10), para las porosidades logarítmicas, los parámetros conocidos son las porosidades en el pozo y el valor desconocido es la porosidad en un punto del volumen, que se encuentra en la posición de una traza sísmica en el volumen, entonces la fórmula resultante para la interpolación de la porosidad logarítmica es:

donde:

φ^*1_krig  es la estimación de la porosidad por Kriging en el punto de la traza.

Del mismo modo se puede emplear la fórmula para estimar la desviación de la impedancia respecto a la media, que en este caso proviene de la transformada de Wyllie de la porosidad por Kriging ( φ^*1_krig). Esto se puede representar como:

donde:

Z_wk es la transformada de Wyllie de la porosidad de Kriging en el punto  x₁ , Z_p es la impedancia del pozo, Z_wp es la transformada de Wyllie de la porosidad del pozo y Z_krig es la estimación de la impedancia del Kriging.

La varianza asociada a la interpolación está definida por la expresión:

donde:

var_krig es la varianza de la estimación de mínimos cuadrados (Kriging) en el punto desconocido y σ es la desviación estándar previa de la porosidad o de la desviación de la impedancia con respecto a la transformada de Wyllie según sea el caso. Los valores de las dos desviaciones previas seextraen de los datos de registro de pozos. La varianza delvalor estimado por Kriging es menor para puntos delvolumen cercanos al pozo y mayor en la medida en que el punto se aleja del pozo.

Con base en lo antes expuesto, se interpolan los valores de la porosidad logarítmica y la impedancia acústica en todo elvolumen, a partir de los registros de pozo disponibles en elárea. Por la condición que se mencionó de la varianza, estosmodelos son idénticos a los del pozo si la traza se sitúa cerca o sobre el pozo y son iguales a la media previa si la traza se ubica más allá de la distancia de correlación establecida en el modelo estadístico.

Para simular los datos sísmicos se utilizó un método de convolución, para lo cual fue necesario estimar a partir delos datos la ondícula fuente y la serie de reflectividadcorrespondiente al modelo de impedancias acústicas. Paraobtener la serie de reflectividad, era indispensable disminuirla resolución con la que se estaba trabajando, ya que comose dijo, la señal sísmica emplea capas litológicas de mayor espesor que las correspondientes a los registros de pozo. Por tanto se aplicó la fórmula modificada de Backus (3), para promediar la impedancia acústica a la escala de la informaciónsísmica, y a partir de esta serie de capas se calcularon los correspondientes coeficientes de reflexión y transmisión en incidencia normal.

La inversión por optimización tiene por finalidad calcular un arreglo de parámetros del modelo tal que minimice la función objetivo. Para este proyecto que integra informacióngeológica, geofísica y petrofísica se presenta una función objetivo compuesta de la forma siguiente:

donde:

el primer término en la función objetivo, S₁, corresponde a la información geofísica y mide la proximidad entre los datos sísmicos observados y los calculados. En la fórmula anterior,  σ²_d es la varianza de los datos, g(Z) es la función que resuelve el problema directo de cálculo de la traza sísmica, d^obs son las amplitudes sísmicas observadas, y Z es el modelo de impedancia.

donde:

el segundo término de la función objetivo corresponde a la información petrofísica y geoestadística sobre la impedancia, y evalúa conjuntamente la proximidad entre las impedancias, y entre el valor predicho a partir de la fórmula de Wyllie en función de la porosidad y el valor interpolado a partir de los registros de los pozos. En la anterior expresión, f(φ*)es la transformada de Wyllie, φ* es el modelo deporosidad, Z_krig es el modelo de impedancia obtenido por interpolación de los datos de pozos,  φ*_krig es el modelo de porosidad obtenido por interpolación de los datos de pozos, C⁻¹z_krig es la inversa de la matriz de covarianza para la impedancia.

donde:

el último término corresponde a la información sobre la porosidad logarítmica con el valor estimado por interpolación de los datos de pozos, donde   C⁻¹_φ es la inversa de la matriz de covarianza para la porosidad. Los modelos interpolados a partir de los datos de pozos, Z_{k rig} y  φ*_krig , son independientes de la sísmica y por lo tanto son constantes durante el proceso de inversión.

El método de optimización que se aplicó fue el de Newton (Tarantola, 2005), que emplea tanto el gradiente como unaaproximación de primer orden del Hessiano de la funciónobjetivo. En los puntos minimales de la función el gradientees cero. Esta condición se usó para obtener las ecuaciones lineales para el problema de inversión. Expandiendo el gradiente de la función objetivo:

e igualando el lado izquierdo de la igualdad a cero, se obtienen las ecuaciones normales como:

En las dos últimas expresiones, Δm es el paso de actualización en el espacio de parámetros, Δm =m_n+1-m , adecuado para minimizar la función objetivo, H(S)[m_n] es el Hessiano de S evaluado en la configuración vigente delmodelo, y ∇S [m_n] es el gradiente de S evaluado en elmodelo vigente. Resolviendo Δm en este sistema deecuaciones lineales, es posible calcular el paso óptimo para minimizar la función objetivo.

El cálculo del gradiente de la función objetivo resulta como:

G=(∂g/∂z) y F=(∂f/∂φ*) son los gradientes de g(Z) y , f(φ*) respectivamente. Los operadores z ∇ y φ ∇ representan la diferenciación respecto a la impedancia y a la porosidad, respectivamente.

Si las funciones g(Z) y f(φ*) fueran lineales, la expansión del gradiente sería exacta. Pero esto no se cumple en nuestro caso, y la configuración actualizada, n+1 Δm=m −m, esuna configuración más cercana al mínimo de la funciónobjetivo, mas no la configuración mínima. Por lo tanto elprocedimiento se vuelve a repetir iterativamente hasta la convergencia al mínimo.

Multiplicando a ambos lados por la matriz de covarianza conjunta C_m , se produce el sistema:

donde:

la matriz del lado izquierdo de la ecuación es la curvatura de la función objetivo y el arreglo del lado derecho es la dirección de máxima pendiente. Después de incluir lacurvatura en la expresión del Hessiano y la matriz de covarianza conjunta, el sistema toma una estructura triangular por bloques de la manera siguiente:

Más detalles sobre la expresión anterior se describen en el trabajo de Bosch (2004). De esta forma, el sistema de ecuaciones conjuntas para las dos propiedades se reduce a la de dos sistemas de ecuaciones menores. La segunda línea en la última fórmula, es un sistema independiente de laprimera línea, que permite resolver Δm_z . Esta solución para el paso de la actualización de la impedancia, se puede usarpara resolver el paso, Δm_φ , en el sistema de la primera línea.Haciendo las sustituciones para la parte derecha del sistemacon las expresiones para la dirección de máxima pendiente yla matriz de covarianza conjunta, el paso de la optimización para el modelo de impedancia es:

y:

Para la actualización del modelo de porosidad, la ecuación que se obtiene de la primera línea del sistema lineal conjunto es:

El paso de actualización de los modelos para una iteración provienen de las expresiones (24-27). Como el problema es no lineal, g(Z), G y F son dependientes del modelo vigente, es decir, de la iteración.

El estudio en un área de ensayo en la región oriental de Venezuela se realizó a una línea sísmica de datos apilados y migrados en tiempo. La señal sísmica fue muestreada cada 2ms, el tiempo total de medición fue de 3s. También, se utilizaron los registros de velocidad de onda P, densidad y neutrónico de un pozo ubicado a 20m de distancia sobre la traza de la línea. A partir de los registros de densidad y neutrónico se derivó el registro de porosidad total.

A partir de los registros mencionados, se obtuvieron los de porosidad logarítmica y de impedancia acústica que luego se usaron para realizar la inversión. No se utilizó la totalidadde los registros del pozo, sólo se usó una sección de interésubicada entre los 1100 y 1600m de profundidad; se escogiódicha sección porque contenía a la Formación Oficina, que es la formación productora en este reservorio.

Esta división se realizó con el propósito de evitar correlacionar propiedades de puntos que estén en distintosbloques, ya que por ser una falla la causante de estaseparación, las capas a un lado y otro de la falla a la mismaaltura no deben ser las mismas y los valores de los parámetros no están correlacionados. En otras palabras, la matriz de correlación para un punto que se encuentra en el bloque A, presenta correlación cero con los valores que están en el bloque B, información que forma parte de la constitución de la matriz de covarianza del medio en las fórmulas (25-27).

Se estimaron los parámetros para la transformada petrofísica, fórmula (6), es decir las velocidades Vmat y Vflu, y las densidades ρ_mat y ρ_flu que caracterizan la matriz y el fluido. Para ello se desarrollaron los cálculos correspondientes utilizando un método de regresión no lineal. Los valoresque se obtuvieron después de la regresión aplicada a los datos de los registros de pozo fueron, V_mat=5728 m/s,V_flu=1622 m/s, ρ_mat=2953kg/m³, y ρ_flu=1285kg/m³. La figura 1 muestra el conjunto de datos de registros de pozo usados para la regresión y la transformada petrofísica de la impedancia en función de la porosidad utilizada en la inversión.

Figura 1. Gráfico cruzado de la impedancia acústica y la porosidad logarítmica. Línea negra representa la transformada de Wyllie con parámetros calibrados al pozo por regresión.

Con esta caracterización concluida, se procedió a generar el modelo de porosidad e impedancia, a alta resolución vertical, que corresponde a la interpolación de los datos de registrosde pozo. Se utilizó una distancia de correlación lateral de 600m en la función de covarianza espacial. La función decovarianza vertical para la porosidad logarítmica y para la desviación de la impedancia respecto a la impedancia predicha por la transformada petrofísica se modeló a partir de los datos del registro de pozos.

A partir del pozo, se extrajo la media de la porosidad logarítmica que resultó igual a -0.465, con una desviación estándar de 0.409; el otro valor que se calculó fue la desviación de la diferencia entre la impedancia observada yla impedancia predicha por la fórmula de Wyllie, cuyo valorfue 1.037 MPa. En la figura 2 se presentan los modelos de porosidad total y de impedancia, respectivamente, generados a partir de esta interpolación geoestadística; en estos modelos se puede notar que cuando la traza está cerca del pozo, el modelo es muy parecido al del pozo, mientras que más lejos tiende a la media previa de cada propiedad.

Figura 2. Secciones de porosidad total y de impedancia acústica (en tonos de grises) resultantes de la interpolación por Kriging de los datos de registros de pozo. Los datos del registro de pozo se superponen sobre el CDP 45 con la misma escala de grises que la propiedad correspondiente calculada con la inversión; la línea punteada blanca indica la posición del pozo. La línea negra gruesa indica la intersección de la sección con la superficie de falla. Las líneas de sombra muestran los datos sísmicos.

Tomando como modelo previo el resultado de la extrapolación geoestadística de los datos del pozo, el siguiente paso consistió en realizar la simulación de los datos,para lo cual se promediaron verticalmente las capas muestreadas a 1ms, mediante la fórmula (7), y se obtuvo la serie de impedancias equivalentes para las capas gruesas (6ms de espesor), con las que se calcularon los coeficientes de reflexión.

A partir del pozo, se extrajo la media de la porosidadlogarítmica que resultó igual a -0.465, con una desviación estándar de 0.409; el otro valor que se calculó fue la desviación de la diferencia entre la impedancia observada y la impedancia predicha por la fórmula de Wyllie, cuyo valor fue 1.037 MPa. En la figura 2 se presentan los modelos de porosidad total y de impedancia, respectivamente, generados a partir de esta interpolación geoestadística; en estos modelos se puede notar que cuando la traza estácerca del pozo, el modelo es muy parecido al del pozo,mientras que más lejos tiende a la media previa de cada propiedad.

Tomando como modelo previo el resultado de laextrapolación geoestadística de los datos del pozo, el siguiente paso consistió en realizar la simulación de los datos, para lo cual se promediaron verticalmente las capas muestreadas a 1ms, mediante la fórmula (7), y se obtuvo laserie de impedancias equivalentes para las capas gruesas (6ms de espesor), con las que se calcularon los coeficientes de reflexión.

Posteriormente se realizó la inversión de los datos, resolviendo iterativamente el sistema lineal descrito en lasección de metodología. En la inversión utilizamos una desviación tipo de los datos calculados respecto a los observados del uno por ciento. En la figura 3 mostramos los datos observados y los calculados a partir del modelo conjunto impedancia-porosidad resultante de la inversión, evidenciando que este modelo explica las observaciones sísmicas dentro del margen de error establecido.

Figura 3. (a) Trazas sísmicas resultantes del apilamiento de los datos sísmicos registrados en campo. (b) Trazas sísmicas calculadas a partir del modelo de impedancias estimado mediante el proceso de inversión. (c) Diferencia entre las trazas sísmicas en los gráficos (a) y (b).

En la figura 4a se presenta el resultado de la inversión para la porosidad total. En este gráfico se observa que existebuena correspondencia entre los valores predichos y losvalores reales del pozo y que mantiene hasta la distancia decorrelación una alta resolución vertical, superior a la correspondiente a la inversión sísmica lejos del pozo.

Figura 4. Secciones de porosidad total e impedancia acústica resultantes de la inversión geoestadística de las amplitudes sísmicas. Los datos del registro de pozo se superponen sobre el CDP 45 con la misma escala de grises que la propiedad correspondiente calculada con la inversión; la línea punteada blanca indica la posición del pozo. La línea negra gruesa indica la intersección de la sección con la superficie de falla. Las líneas de sombra muestran los datos sísmicos.

También se aprecia que la resolución alrededor de la falla varía según la posición del pozo respecto al bloque. Esto sedebe a que la correlación lateral a través de bloques distintoses cero. Otra característica importante, es el predominio deltono de gris correspondiente a la media de los valores y quees más evidente en las zonas donde no existe correlacióncon el pozo; mientras que en las zonas cercanas al pozo, los valores son muchos más dispersos, gracias a que se resuelve una estructura más fina en la estratigrafía, generando así capas más contrastantes en el gráfico.

Comparando ambas imágenes se observa que presentan valores extremos en los mismos puntos, que al ser comparados a su vez con la sísmica, reflejan una buenacorrespondencia, lo que evidencia el buen funcionamiento de la técnica.

En este trabajo se desarrolló una metodología para la inversión de amplitudes sísmicas que teniendo como base las técnicas de optimización genera, de forma iterativa, unmodelo del subsuelo que cumple tanto con los datos sísmicos como con los registros de pozo. Con la técnica se estiman simultáneamente la porosidad y la impedancia del medio,descritas a una resolución vertical superior a la que proporciona una inversión sísmica convencional.

Los autores agradecen a Petrobras Energía de Venezuela por la cooperación en este trabajo mediante el suministro de datos del área y participación de su personal. Igualmente agradecemos el apoyo del CDCH-UCV, a través del proyecto PG-08-00-5631-2004. Este trabajo se desarrolló en el marco del convenio de cooperación UCV-Petrobras.

1. BACKUS, M. (1962). Long-wave elastic anisotropy produced by horizontal layering, Journal of Geophysical Research, 67, 4427.        [ Links ]

2. BOSCH, M. (2003). Inferencia estadística de porosidad e impedancia a partir de ondas sísmicas vía métodos de Monte Carlo y optimización. Trabajo de ascenso. Universidad Central de Venezuela, Caracas.        [ Links ]

3. BOSCH, M. (2004). The optimization approach to lithological tomografy: Combining seismic data and petrophysics for porosity prediction, Geophysics, 69, 1272-1282.        [ Links ]

4. DEUTSCH, C. V. and A. G. JOURNEL (1998). Geostatistical Software Library and User’s Guide, New York, USA, 363pp.        [ Links ]

5. HILTERMAN, F. J. (2001., Seismic amplitude interpretation: Geophysical Development Corporation.         [ Links ]

6. TARANTOLA, A. (2005). Inverse Problem Theory and Methods for Model Parameter Estimation: SIAM publications, Philadelphia.        [ Links ] ]]> 1 1962 67 4427 2 2003 3 2004 69 1272-1282 4 1998 363pp 5 2001 6 2005