1. Introducción

En geofísica, la mejor aproximación que se tiene de la estructura interna de la tierra es por medio de la interpretación sísmica [³³]. Particularmente, se analiza la respuesta sísmica del subsuelo ante una fuerza aplicada y que es registrada como amplitudes positivas, negativas y cruces por cero (Figura 1a), las cuales son asociadas con un mapa de colores (1b). La unificación de las trazas crea un perfil o imagen sísmica (Figura 1c) y un gran número de estas imágenes forman lo que se conoce como un cubo sísmico. Este cubo es interpretado por los geocientíficos para encontrar, entre muchas otras cosas, identificadores de acumulación de hidrocarburo [¹⁶].

Fig. 1 Generación de una imagen sísmica. a) La respuesta del subsuelo a una fuerza artificial es capturada en forma de amplitudes positivas, negativas y cruces por cero. b) Asociación de las amplitudes de las trazas con un mapa de colores. c) Perfil o imagen sísmica resultante cuyos patrones coloridos representan la continuidad de las capas de roca. (1) Patrón sísmico con alta intensidad de color (2) Patrón complejo de forma elongada, irregular y de baja intensidad de color 

Los Sistemas comerciales de Interpretación Sísmica actuales (SIS) permiten tener una visión dinámica y colorida de los perfiles sísmicos, lo que resulta adecuado para apreciar más detalles que con las trazas sísmicas. Al existir una relación directa con las amplitudes, los patrones de color más intensos representan cambios drásticos en la roca (Figura 1(1)), mientras que los menos intensos significan cambios más débiles (Figura 1(2)) [⁶]. El geofísico realiza un picking o mapeo sísmico, el cual consiste en realizar diversos trazos directos sobre los patrones coloridos e identificar diferentes aspectos estratigráficos y estructurales [¹⁶]. Esta tarea depende de la calidad de las señales, por lo que puede ser una actividad trivial o bien, un problema desafiante [⁷]. El objetivo de estos trazos es el de tener una visión tridimensional de cómo y en dónde están distribuidas las capas y de qué manera son afectadas por la presencia de fracturas [¹⁰]. El especialista tiene la opción de realizar el mapeo a mano alzada o mediante los procesos automáticos (autotrackers) disponibles en los SIS [³⁴]. Sin embargo, debido al gran número de imágenes que conforman el cubo sísmico, los trazos manuales consumen demasiado tiempo y los trazos automáticos suelen no satisfacer las expectativas del experto [⁴]. En esta etapa de interpretación, los objetivos primordiales que busca el especialista son los patrones de mayor intensidad (amplitudes originales altas), no obstante, en muchas ocasiones un patrón de interés se encuentra en regiones con baja intensidad de color (amplitudes originales bajas) en donde los trazos manuales son aún más laboriosos y en donde los procesos automáticos fallan [⁷]. En la Figura 2 se muestra que el autotracking (línea amarilla) sobre patrones complejos es incorrecto. Dado un punto inicial indicado por el usuario (2a), pueden observarse brincos (2b) y truncamientos (3c) en el trazo que requieren de ajustes manuales y laboriosos que ocasionan que este tipo de zonas sean poco exploradas.

Fig. 2 Autotracking sobre patrones complejos. La línea amarilla es el resultado de un proceso automático incorrecto. a) Punto inicial del mapeo. b) Invasiones o brincos. c) Mapeo truncado 

En este artículo se propone un algoritmo alternativo de procesamiento sobre patrones de baja intensidad de color mediante técnicas de procesamiento de imágenes. Se busca que en su representación digital sea posible mejorar los trazos que realizan los autotrackers a partir de la extracción de regiones seleccionadas por un usuario. Como mapas de bits, las imágenes sísmicas pueden ser analizadas para realizar una segmentación supervisada de regiones de mayor complejidad. El algoritmo propuesto es capaz de procesar zonas de intensidad de color baja y extraer los elementos necesarios para construir geocuerpos, los cuales son una proyección tridimensional de ciertas secciones de una imagen. Los resultados muestran que los geocuerpos encontrados son tan relevantes como aquellos obtenidos a partir de patrones de intensidad de color alta.

Uno de los objetivos de los algoritmos de segmentación es el de simplificar imágenes para obtener información significativa o para la detección de objetos [²⁷]. Sin embargo, a pesar de que existen cientos de técnicas de segmentación en la literatura, no hay un solo método que funcione para todos los tipos de imágenes [²⁹]. La segmentación consiste en asociar aquellos píxeles que compartan ciertas características (como el color, iluminación o textura) dentro de una imagen [³⁶]. Dicha asociación va creando regiones que son un subconjunto de la imagen original, ya sea objetos, contornos, o incluso la imagen misma.

Dentro del área de procesamiento de patrones, se identifican cuatro técnicas de segmentación: las que binarizan la imagen basadas en un umbral (divide los píxeles en dos grupos: blanco para el fondo y negro como primer plano); aquellas que detectan límites o bordes (buscan los cambios drásticos de la intensidad de los píxeles en los bordes u orillas de los objetos); las que obtienen regiones (encuentra regiones coherentes formadas por píxeles que tienen características similares) y las técnicas híbridas (combina crecimiento de región y detección de bordes) [¹⁵].

Por lo general, estos métodos analizan la imagen completa para llevar a cabo la segmentación de acuerdo con distintos objetivos. Debido a lo anterior, no tienen un buen desempeño si se pretende obtener una segmentación más específica y al aplicarlos en imágenes sísmicas se generan demasiados segmentos basura. El algoritmo propuesto en este artículo efectúa un procesamiento más controlado que reproduce el trazo manual que realizaría un usuario humano sobre la continuidad de los patrones sísmicos.

El presente artículo está organizado de la siguiente manera: en la sección 2 se muestran los trabajos relacionados y la aportación principal de este trabajo; la sección 3 describe el IMP-2DMA, se explica en qué consiste la selección de semillas y guías de direccionamiento; los resultados experimentales se presentan en la sección 4 y finalmente la sección 5 es de conclusiones y perspectivas.

2. Trabajos relacionados

La automatización de tareas sobre imágenes sísmicas ha tomado gran relevancia en los últimos años [⁴²]. La simplificación de las imágenes como esqueletos sísmicos (líneas), tuvo una aceptación dividida dentro de la industria [²², ²³], sin embargo también ha sido objeto de análisis en el dominio del procesamiento de imágenes. La idea principal consiste en aplicar algoritmos de detección de contornos para simplificar la imagen pero conservando su estructura original (Figura 3a) [²⁰, ³]. No obstante, los resultados suelen ser confusos y requieren de edición e interpretación por parte del experto. A pesar de eso, la simplificación de las imágenes sísmicas dio lugar al planteamiento de nuevos objetivos, como la identificación de secciones sísmicas (facies) [¹⁹, ⁴⁴, ¹³], el reconocimiento de fracturas [¹², ³², ⁴³, ⁴¹, ¹⁴] y la identificación de curvas asociadas a cuerpos salinos (3b, línea amarilla) [⁴², ³⁸, ²⁴, ¹] o círculos pequeños asociados a dolinas kársticas (3c) [³⁰] entre muchos otros.

Fig. 3 Procesamiento de imágenes sísmicas. a) Detección de contornos. b) La detección de curvas favorece la identificación de cuerpos salinos. c) Los círculos representan estructuras geológicas 

La opciones de color que surgieron en los SIS amplió el panorama de interpretación y fue posible apreciar más detalles presentes en los cubos sísmicos [⁴]. La visualización de geocuerpos, que son la conjunción tridimensional de las imágenes 2D, fue posible al proyectar automáticamente las capas individuales de ciertos patrones seleccionados. Sin embargo, como ya se ha comentado, los patrones mostrados en las Figuras 1(1), 1(2) y Figura 2 no pueden ser procesados de forma trivial ni por los expertos ni por los autotrackers. Las distintas técnicas de segmentación basadas en umbrales [⁴⁰, ³⁹], detección de contornos [⁵, ¹⁸], crecimiento de región [³⁵, ¹⁷] y las técnicas híbridas [¹¹, ²¹] tampoco son una opción viable debido a que aplicados en este tipo de imágenes generan demasiados segmentos basura. El algoritmo propuesto realiza un mejor procesamiento de los patrones sísmicos, independientemente de su forma irregular, elongada y las intensidades de color de la que estén formados.

3. El algoritmo desarrollado: el IMP-2DMA

En este artículo se propone un algoritmo de segmentación llamado IMP-2DMapping Algorithm o IMP-2DMA que actúa en imágenes sísmicas a color dentro del espacio CIELAB. El algoritmo está divido en tres partes: la primera es la selección de semillas, rangos y variables. La segunda parte consiste en la creación de máscaras binarias que muestran la segmentación realizada. En esta etapa se muestra el algoritmo principal de todo el procesamiento. La última parte muestra cómo se construyen y visualizan los geocuerpos a partir de las máscaras obtenidas en la etapa anterior.

El IMP-2DMA fue evaluado con imágenes que corresponden a cubos sísmicos reales y fueron obtenidas directamente del sofware Petrel R [³⁷]. Las imágenes tienen una resolución de 1680x600 bit/píxeles. Se seleccionó el espacio de color CIELAB (creado por la Commission Internationale de l’Eclairage en 1976 [⁸]) por ser ampliamente recomendado para realizar comparaciones entre colores, porque puede mostrar todos los colores visibles al ojo humano y además de ser independiente de dispositivo [⁹]. La ventaja de CIELAB sobre el estandar RGB (espacio que define los colores en términos de los colores primarios rojo, verde y azul) es que se logran identificar diferencias muy pequeñas y significativas entre colores [²⁵].

3.1. Selección de semillas, rangos y variables

Una semilla se define como un punto de coordenadas (x, y) dentro de una imagen (indicado por el usuario) que se convierte en el punto inicial del procesamiento, ya que define el color predominante del patrón seleccionado. En el espacio CIELAB, todo imagen consta de tres coordenadas (L, a, b), que corresponden a la iluminación L y a las coordenadas cromáticas a (coordenadas rojo/verde) y b (coordenadas amarillo/azul). El usuario selecciona un patrón sísmico y define rangos y variables necesarias para el funcionamiento del algoritmo.

En la Figura 4 se muestran los valores iniciales que selecciona el usuario. La semilla central s define la posición y el color predominante del patrón inicial. Las flechas alrededor indican un vecindario de 4 u 8 píxeles (a elección del usuario) que se convierten, a su vez, en un conjunto de semillas base. El usuario determina también el vecindario V , que sirve de límite al crecimiento vertical a partir de un punto p. A la diferencia de color entre dos muestras se le conoce como ΔE o error delta y cuando se aplica en valores CIELAB se escribe ΔE*. Este valor permite conocer la diferencia entre dos colores e involucra la iluminación y los valores de croma, definidos como: ΔL*, Δa* y Δb*.

Fig. 4 Rangos y variables iniciales. La semilla central s define el color predominante en el patrón. Los valores ΔE*, el umbral δ y el vecindario vertical V son definidos por el usuario de acuerdo a las características del patrón seleccionado 

Para calcular la diferencia entre dos píxeles p(x, y) y q(x, y), se calcula este valor con la fórmula CIE76 [³¹]. Sean (L1*,a1*,b1*) y (L2*,a2*,b2*) dos colores en el espacio CIELAB, ΔE* se calcula con:

(1)

siendo ΔE* ≈ 2,3 un valor apropiado para notar lo que se conoce como una diferencia apenas notable (Just Noticeable Difference). Si este valor supera un umbral δ, indica que la diferencia de los colores es evidente. Dada la distribución casi horizontal de los patrones sísmicos, se propusieron dos guías de direccionamiento para conseguir un mejor control del recorrido a través de las imágenes sísmicas. Los puntos de la guía G ₁ = {p ₁, p ₂, . . . , p _n } deben de seguir aproximadamente la trayectoria del patrón que seleccionó el usuario (Figura 5). Dados estos puntos, se realiza un ajuste de curva de interpolación para obtener los puntos intermedios. El objetivo de las guías de dirección es el de resolver las discontinuidades existentes sobre los patrones sísmicos seleccionados.

Fig. 5 Guías de direccionamiento. El usuario define una guía G ₁ = {p ₁, p ₂, . . . , p _n } en donde cada punto de la guía debe seguir de manera aproximada al patrón sísmico seleccionado 

Debido a que los patrones que conforman una imagen sísmica van cambiando gradualmente de forma y posición (desplazamientos verticales), la guía G ₂ = {p ₁, p ₂, . . . , p _n } determina los cambios de altura que tiene un patrón en particular que haya seleccionado el usuario. A manera de ejemplo, en el perfil sísmico de la Figura 6 se muestra un patrón de color azul (6a, vista frontal) que tiene un desplazamiento hacia abajo hasta llegar a una posición final (6a, vista en perspectiva). El usuario marca los puntos de la guía G ₂ únicamente sobre la vista frontal acorde al desplazamiento vertical del patrón. Al igual que con G ₁, se realiza el cálculo de los valores intermedios para tener una altura inicial en cada imagen. Las dos guías permiten que este patrón sea extraído en todas las imágenes desde su posición inicial hasta la final.

Fig. 6 Desplazamiento vertical de un patrón sísmico desde una posición inicial (a) a una final. b) Los puntos de G ₂ definen la altura del recorrido sobre cada uno de los puntos de G ₁ 

3.2. Creación de máscaras binarias

El IMP-2DMA recibe como entrada un conjunto de imágenes I y devuelve un conjunto M de máscaras 2D binarias. Los pasos del algoritmo IMP-2DMA se muestran en la Tabla 1. Todas las variables son almacenadas en una colección de objetos O _j . El recorrido de la imagen va de izquierda a derecha sobre los puntos de G ₁ y G ₂. Se lleva a cabo la conversión de la i-ésima imagen al formato CIELAB con I _c = C(i) (punto 3). Por cada píxel p _i se realizan dos validaciones verticales (V ,′ ↑′,′ ↓′) de dicho píxel mediante la función α(s, V , I _c , δ), donde s ← γ(p, g) es un vector de semillas que incluye a los vecinos inmediatos, g = 4 o g = 8, de p. Cada nuevo píxel q ∈ V es evaluado contra todos los elementos de p ∈ s.

Tabla 1 Algoritmo principal del IMP-2DMA. Cada objeto recibe un conjunto de imágenes a ser procesadas y sus respectivas variables. Por cada punto de las guías G ₁ y G ₂, se realiza la validación del vecindario V para crear las máscaras binarias 

La función α es la que realiza la comparación entre píxeles cuyo algoritmo se muestra en la Tabla 2: recibe el vector de semillas s, V y δ y realiza la comparación entre cada píxel q ∈ V contra los elementos de p ∈ s. Si se cumple que 2,3 ≤ ΔE*(p, q) ≤ δ, se considera un píxel aceptado y se marca con un valor binario en M mediante set(M, p(x, y),′ 0′), donde M es la imagen de salida o máscara binaria.

Tabla 2 Algoritmo de la función α. Dado un píxel q ∈ V , se evalúa contra los píxeles p ∈ s con la ecuación 1. Todo punto que cumpla que ΔE*(p, q) < δ, es un píxel aceptado y será agregado como parte de la misma región 

El algoritmo termina cuando se han evaluado todos los píxeles de G ₁. En la Figura 7 se muestra un resultado preliminar de la máscara binaria que se obtuvo de un perfil sísmico. La guía G ₁ está distribuida en un patrón de grises a lo largo de la horizontal de la imagen (7a). La máscara obtenida (Figura 7b) es una extracción del patrón seleccionado mediante el algoritmo IMP-2DMA. El proceso se repite en cada una de las imágenes del cubo sísmico, obteniendo un conjunto de máscaras que servirán para la visualización tridimensional.

Fig. 7 a) Selección de un patrón sísmico de color grisáceo y de estructura alongada. Se muestra la serie de puntos de la guía G ₁ que siguen la trayectoria del patrón sísmico y el alcance del vecindario vertical V . b) Máscara binaria 2D obtenida 

3.3. Visualización de geocuerpos

Para visualización de los resultados se necesita construir una proyección tridimensional de las máscaras 2D obtenidas. Se utilizaron funciones ya diseñadas en Matlav v9. Como primer paso, se construye un espacio (meshgrid) para en donde se alojarán cada una de las máscaras m _i ∈ M. Se procesa cada máscara para obtener su respectivo contorno (Figura 8a), se realiza un enlace punto a punto con los contornos de la siguiente imagen (8b) y finalmente, se agregan efectos de luz y sombras para dar un efecto de volumen. La utilidad de estos geocuerpos, que aún deben de ser interpretados, queda fuera del alcance de este artículo, sin embargo, sí es posible validar de una manera cuantitativa las máscaras binarias que se obtuvieron.

Fig. 8 Creación y visualización de geocuerpos. a) Se crea un espacio tridimensional para alojar cada una de las máscadas 2D obtenidas. b) y c) El geocuerpo correspondiente se obtiene al entrelazar las máscaras y agregar efectos de luz y sombra para dar un efecto de volumen 

4. Resultados experimentales

Para evaluar el desempeño del algoritmo, se realizaron diversas pruebas en las imagenes de un cubo sísmico real y se analizaron las máscaras binarias obtenidas. Se seleccionaron tres patrones distintos y visualmente atractivos de forma irregular y colores de baja intensidad (9, columna (a)). Las distintas técnicas de segmentación que existen dividen una imagen de tal modo que pueden contener demasiadas (oversegmented) o pocas regiones (undersegmented) correspondientes a objetos dentro de la imagen [²⁸].

Aplicadas en una imagen sísmica, ocurre que cualquiera que sea la técnica utilizada, el resultado es un exceso de regiones que dificultan la extracción y análisis de una sección específica. Debido a que la salida a evaluar son imágenes binarias, se seleccionó un algoritmo de binarización basado en umbral por ser uno de los más sencillo de implementar y de ajustar. En este tipo de pruebas, el trazo realizado manualmente se toma como el más acertado. En la Figura 9 columna (b) se observa la solución ideal para cada muestra y en las columnas (c) y (d) los resultados que devuelven ambos algoritmos. Puede observarse que los resultados del IMP-2DMA (columna (c)) tienen una mayor precisión y no generan tantos segmentos ajenos a la solución ideal como los del algoritmo basado en umbrales (columna (d)).

Fig. 9 Resultados de la segmentación. a) Imagen original. b) Segmentación manual. c) Segmentación obtenida por el IMP-2DMA. d) Segmentación basada en umbral. El rectángulo rojo acota el área de comparación de este algoritmo contra la segmentación manual 

Para realizar un análisis cuantitativo, se utilizó la prueba de suma de rangos de Wilcoxon [²⁶] en lugar de otro análisis como la prueba-t [²]. Lo anterior es debido a que el tamaño de una de las muestras puede llegar a ser muy grande respecto a la otra, lo cual es mejor manejado por la prueba de Wilcoxon.

Dadas las máscaras obtenidas por los distintos métodos, M _HU (segmentación manual), M _T H _R (segmentación por umbral) y M _IMP (segmentación del IMP-2DMA), se plantea la hipótesis nula de que las máscaras sean iguales, simbolizada como:

(2)

Siendo M ₁ la segmentación manual y M _i , con i = 1, 2 donde 1 es el M _T H _R y 2 es el M _IMP (Figura 9). Es evidente que el algoritmo basado en umbral, que actúa sobre toda la imagen, genera demasiados segmentos basura que afectan su comparación con los otros dos métodos (siendo el rectángulo rojo el que acota el área de evaluación). En la Tabla 3 se reportan los resultados experimentales. Se muestran los valores de p obtenidos con α = 0,05. Es claro que los valores por p > 0,05 concluyen que la hipótesis nula H ₀ no puede ser rechazada, por lo que el resultado del M _IMP es superior al del algoritmo basado en umbral.

Tabla 3 Resultados experimentales. Se muestran los valores de p al comparar los resultados del M _THR y el M _IMP contra la segmentación manual M _HU  

5. Conclusiones y perspectivas

En este artículo se ha propuesto un método de segmentación de imágenes a color basado en el espacio CIELAB: el IMP-2DMA. El algoritmo está enfocado en la extracción de patrones complejos visibles en las imágenes sísmicas. La proyección de las máscaras 2D obtenidas permite visualizar los geocuerpos que surgen a partir de patrones de mayor complejidad. El algoritmo fue evaluado con diversos patrones sísmicos con información de cubos sísmicos reales.

Este trabajo es una aportación al área de la interpretación sísmica, ya que cumple con un objetivo más amplio de lo que existe en la literatura. Dentro del área de reconocimiento de patrones, se encontró que los esfuerzos actuales se han enfocado en extraer exclusivamente cuerpos de sal o domos salinos [²⁴, ¹, ¹³], dejando relegado el mejoramiento de técnicas para una mayor y completa exploración de las imágenes sísmicas. Desde luego que no se trata de competir contra lo que hacen los paquetes especializados, sino de mostrar que es posible realizar procesos exploratorios más profundos con un mínimo de información disponible.

Una de las ventajas del IMP-2DMA es que es posible el procesamiento y proyección tridimensional de cualquier patrón, independiente de su forma y su intensidad de color, superando algunos de los resultados que se obtienen con los SIS. Sin embargo, los geocuerpos extraídos aún deben de ser interpretados por los especialistas para darles el adecuado contexto geológico. El algoritmo aún debe de mejorarse para aumentar la velocidad de procesamiento de las imágenes. Como trabajo futuro se contempla aplicar el IMP-2DMA a otro tipo de imágenes para la recuperación de fracturas en fotografías de rocas, ya que se puede extraer su trayectoria completa, conservando su espesor mediante las técnicas propuestas en este trabajo.