2.1 Trazado de rayos

El problema de la prospección sísmica consiste en, dado un conjunto de observaciones en la superficie asociadas con la propagación de ondas en el subsuelo, se busca identificar cuáles son las características estructurales del subsuelo que se hacen responsables por el observable. La comparación de los resultados de las observaciones con los resultados producidos por modelos computacionales del mismo fenómeno, permiten inferir las características estructurales del subsuelo.

Para este fin se requiere entonces de un mecanismo de modelamiento de la propagación de ondas sísmicas (modelador) en un modelo del subsuelo y una estrategia de inversión (inversor) que permita, a través de la comparación entre los datos modelados y los observados, estimar las características (o parámetros) de un modelo que permitan ajustar los datos observados.

Dado que el fenómeno que se observa está asociado con la propagación de ondas sísmicas en el subsuelo, el modelo debe dar cuenta por este fenómeno físico. La propagación de ondas en el subsuelo se modela formalmente a través de la solución a la ecuación de onda

Donde P(x) es la perturbación en el campo de presión del medio y c (x) su velocidad de propagación. La solución a esta ecuación diferencial parcial es complicada, y en general debe ser resuelta numéricamente. Esto conduce a un problema de naturaleza computacional cuando se requiere explorar soluciones rápidas pero útiles del problema. Una forma de aproximar la solución a la ecuación de ondas sin requerir una solución completa del problema, es hacer uso de la aproximación a altas frecuencias que conduce a la aproximación del rayo (^{Cerveny, 2001}). En esta aproximación, se estudia el comportamiento de ondas de alta frecuencia (longitud de onda corta) que en últimas terminan siendo descritas por la dirección de propagación de la onda. Esta solución se caracteriza por la construcción de rayos, que son construcciones geométricas perpendiculares al frente de onda y que indican la dirección de propagación de las ondas. En esta aproximación, la propagación del frente de ondas asociado a un conjunto de rayos está descrita por la ecuación Eikonal (^{Margrave, 2001}).

De la cual se deduce la relación fundamental para el tiempo de propagación del frente de onda en la forma de la ecuación del rayo

donde s(x) es el recíproco de la velocidad de propagación de la onda, o slowness. En términos generales, para la tomografía de tiempos de viaje, esta ecuación es la ecuación fundamental a resolver. En esta situación, dado el observable T(x) asociado a los tiempos de propagación de los primeros arribos de los frentes de onda, la pregunta es, ¿cuál es la función del slowness del subsuelo que reproduce el observable?

El trazado de rayos se basa entonces en la idea de que la onda sísmica en las altas frecuencias sigue una trayectoria determinada por la trayectoria de rayos. Físicamente, esta aproximación describe cómo se propaga la energía de la onda considerando refracciones y reflexiones inducidas por las variaciones de la velocidad de propagación de la onda en el medio (^{Cerveny, 2001}; Udías, 2003).

La solución de la ecuación Eikonal provee los tiempos de propagación de las ondas P, ya que estos están asociados con la energía de la onda que se propaga directamente desde la fuente hasta el receptor. Además, la trayectoria encontrada numéricamente usando la ecuación Eikonal tiene una propiedad importante dada por el principio de Fermat: El camino del rayo es una curva a lo largo del espacio cuyo tiempo de propagación es estacionario. Esta propiedad puede ser explotada para motivar la construcción de un rayo entre una fuente sísmica y un receptor. En este caso su construcción se realiza a partir de una analogía entre el camino de un rayo sísmico y el camino más corto en un grid de puntos. Una estrategia como esta es denominada el método del camino más corto para calcular rayos sísmicos (^{Moser, 1991}).

Para los propósitos de este trabajo la construcción del rayo implicará diferentes etapas. Primero, una discretización del modelo del subsuelo a través del cual se estudiará la propagación de los rayos. A continuación, se describen las diferentes partes de nuestra implementación del método del camino más corto para el trazado de rayos. Finalmente, con el fin de poner en funcionamiento el modelador en el contexto de la tomografía sísmica, se presenta la estrategia de reconstrucción algebraica con la que se usó el modelador para resolver problemas de tomografía sísmica.

Construcción del grid

La discretización del problema inicia con la definición de la caja o dominio computacional que contendrá el medio modelo a través del cual se estudiará la propagación de las ondas. Definido el tamaño de la caja en la horizontal y en profundidad, se construye un grid grueso de tamaño N _x elementos en la horizontal y N _z elementos en profundidad. Cada celda en este grid tendrá como atributo un valor constante del slowness, y como tal representará un elemento del modelo. Sin embargo, estas celdas (que en lo sucesivo llamaremos celdas de tomografía) deben ser bastante grandes para asegurar contener un número suficientemente alto de rayos por celda, tal que permitan conseguir información suficiente a la hora de hacer un proceso de inversión. Esto, sin embargo, afecta el detalle con el que se puede reconstruir diferentes rayos. Para resolver este problema, lo que se hizo fue construir un subgrid que discretiza el interior de cada celda de tomografía. Sobre este subgrid el trazado de rayos se ejecuta de manera mucho más fina, permitiendo reconstruir la trayectoria del rayo con mucha más precisión sin afectar la iluminación de las celdas de tomografía.

2.2 Método del camino más corto

En el método del camino más corto, el modelo del subsuelo está representado por un grid, consistente de nodos conectados por aristas. Los grids usados para representar el modelo de velocidades son sparse; es decir, cada nodo está conectado con un número restringido de nodos en su vecindad y no con todos los nodos que están cercanos a él. Estas vecindades son denominadas listas de adyacencia de los nodos y contienen todos los nodos tal que existe una arista que los conecta, si tal arista no existe en la matriz de adyacencia hay un cero.

Para conocer los vecinos con los cuales se conectará un nodo se puede usar un esquema de vecindad rectangular en donde los vecinos se escogen a partir de un radio de propagación constante para el i -ésimo nodo y trazando la arista entre los vértices de la vecindad y el nodo en cuestión (^{Moser, 1991}; ^{Monsegny y Agudelo, 2013}). En la Figura 1 se muestran vecindades de radio 1, 2, 3 y 4. Para aumentar la cobertura de los ángulos de propagación del rayo se escoge un radio de propagación mayor como se ilustra en la Figura 1.

Figura 1 Vecinos de un nodo especifico según el radio de propagación 

Así las cosas, dado un grid y un radio de propagación se tendrá un conjunto de posibles vecinos para cada celda. Durante la propagación de un rayo que pasa por una celda se tendrá entonces que una vez el rayo sale de una celda, solo podrá propagarse a una celda vecina, esta celda estará identificada a través de la lista de vecinos.

2.2.2 Asignación de pesos

Dado que un nodo tiene un número finito de conexiones con puntos que están cercanos a él, es posible viajar desde un nodo a otro a través de estas conexiones. Para asignar el tiempo de propagación (peso de la arista) entre los puntos que se conectan, se puede hallar la distancia euclidiana y multiplicarla por un promedio de slowness entre los nodos implicados. Es decir, si el nodo i de coordenadas espaciales (x,,z) se conecta con el nodo j con coordenadas (x _j ,z _j ), entonces el tiempo de propagación t¡, entre los nodos i j, está dado por

donde y s¡ es el slowness del nodo .

2.2.3 Búsqueda del camino más corto Dijkstra y reconstrucción del rayo

Luego de tener un grid que consta de nodos, aristas y pesos, se procede a encontrar el camino más corto entre un nodo y otro a partir del algoritmo de (^{Dijkstra, 1959}). El algoritmo realiza la búsqueda de un camino mínimo desde la fuente a cualquier otro punto del grid. Para esto recorre todos los nodos del grid y asigna a cada nodo visitado un tiempo de propagación acumulado TP desde la fuente. En otras palabras, TP[j] almacena el mínimo tiempo que se puede encontrar entre el nodo j y la fuente. Este tiempo puede cambiar si al unir un vecino u del nodo j, se cumple que el tiempo acumulado en TP[u] más el tiempo de propagación desde u a j es menor que el tiempo acumulado en TP[j], es decir, si

entonces el tiempo en el nodo j es actualizado como

Cuando esto ocurre, aparece el vector de predecesores, el cual almacena el identificador del nodo predecesor que está minimizando el tiempo de propagación desde la fuente. Si la Ecuación (6) se cumple, entonces el nodo predecesor que minimiza el tiempo de propagación desde la fuente hasta el nodo j es u, y así, predecesor[j] = u. Cuando todos los nodos tengan un tiempo de propagación definitivo se cuenta con toda la información necesaria para reconstruir el rayo entre una fuente y un receptor. Una de las ventajas de este método es que encuentra el camino más corto a todos los nodos del grid y cualquier nodo diferente a la fuente puede ser escogido como un receptor.

Para conocer las coordenadas espaciales por las cuales cruza la trayectoria del rayo se usa el vector predecesor. Cada componente de este vector tiene almacenado el vértice del cual proviene el mínimo tiempo desde la fuente, por ejemplo, en el rayo de la Figura 2 el tiempo que es mínimo desde la fuente i hasta el receptor j con respecto a otros caminos es t _ia + t _ab + t _bc + t _cd + t _de + t _ej y la información contenida en el vector predecesor es; predecesor[ j] = e, predecesor[e] = d ^... predecesor[a] = i.

Figura 2 Reconstrucción del rayo 

2.3. Tomografía

2.3.1. Tomografía de tiempos de propagación

Uno de los análisis de tomografía más básicos, y que comúnmente se realiza al principio de un proceso de exploración sísmica, es la tomografía de tiempos de propagación (Traveltime Tomography), y en la cual solo se usan los primeros tiempos de llegada de una onda entre una fuente y un receptor. Para determinar la distribución de velocidades en un medio 2D dividido en celdas con velocidad constante c _i , la tomografía intenta resolver un conjunto de ecuaciones simultáneas. Si el medio está dividido en N celdas, es posible representar los tiempos de propagación como:

Donde t es el tiempo total de propagación para el -ésimo rayo, d _j es la longitud recorrida del -ésimo rayo en la j-ésima celda del modelo de velocidades y s _j =1/c es el slowness de la celda y N = N _x x N _z es la cantidad de celdas del modelo (^{Jones, 2010}). Este problema se puede escribir en notación matricial considerando M rayos como T = DS:

Muchos de los elementos de la matriz D son cero, pues un rayo no atraviesa todas las celdas del modelo, lo que convierte a la matriz D en sparse y mal condicionada. El objetivo es estimar S, la matriz que contiene el modelo de slowness, que en este caso es el vector de parámetros del modelo que se quiere invertir.

Para resolver este tipo de problemas, se puede utilizar métodos iterativos, tales como el método del gradiente conjugado (^{Scales, 1987}) y las técnicas de reconstrucción algebraica (^{Lo y Inderwiesen, 1994}; ^{Aster, Borchers y Thurber, 2018}), que funcionan bien en sistemas a gran escala (^{Jones, 2010}). En este trabajo consideraremos el uso de técnicas algebraicas como una aproximación práctica al problema de tomografía de exploración.

2.3.2 Técnica de reconstrucción algebraica iterativa simultanea (SIRT)

La técnica de reconstrucción algebraica iterativa simultanea (SIRT por sus siglas en inglés) es un algoritmo de inversión iterativo en el cual es necesario recorrer todos los hiperplanos (ecuaciones) y determinar las correcciones de cada hiperplano (ecuación) a una celda, para luego actualizar el modelo como un promedio de estas correcciones. Esta forma de actualizar el modelo convierte al método de SIRT en uno de los mejores procedimientos de tomografía cuando el ruido sigue una distribución Gaussiana, además, que es altamente resistente a datos outliers (^{Dobróka y Szegedi, 2014}). Las iteraciones se pueden apreciar geométricamente iniciando con una aproximación inicial G₀, y para encontrar la siguiente aproximación a la solución, se usa un promedio pesado entre las correcciones G/, G _t " y realizadas por los correspondientes hiperplanos respectivamente. El procedimiento continúa de esta forma generando una secuencia G₀, G _v G₂, ^... que converge a la solución X dado un criterio de tolerancia.

Para hallar la actualización de una celda del modelo es necesario calcular previamente el ajuste realizado por todos los rayos, es decir, A' s _j para /=1,--, M. y estimar un promedio pesado entre las correcciones diferentes de cero (ASj=0.0 si el rayo i no pasa por la celda j). Así la actualización para la celda j se puede determinar como:

para j=1,^…, N.

Donde W _j es el número de rayos que pasan por la j -ésima celda. El nuevo modelo estimado se actualiza a partir del promedio de correcciones As _j como:

Donde a es un parámetro que hemos introducido en este trabajo para controlar la amplitud de los cambios en los valores de los parámetros. Si As es muy grande puede pasar que la actualización del modelo salte muy lejos de la solución del problema, a se encarga de controlar la amplitud de estos cambios, entregando un mayor control sobre la forma como el problema converge a su solución a un costo computacional más alto. Después de hacer varios experimentos, encontramos que valores de a = 0.1 permiten al sistema converger suavemente a la solución a un costo computacional aceptable.

3.1 Trazado de rayos Refracciones

Dadas las condiciones del problema que se estudia, lo primero es estudiar el comportamiento de la estrategia propuesta para la realización del trazado de rayos en un medio cuya velocidad aumenta con la profundidad de forma lineal, y que puede ser expresada como una función lineal v = a + bz. Hay que aclarar que se pueden usar otros modelos que posean algún cambio en profundidad, aunque para algunas situaciones el gradiente lineal puede ser una aproximación aceptable a la realidad. En este tipo de experimentos se evidencia el comportamiento del trazador de rayos en un problema completamente dominado por procesos de refracción.

En la Figura 3 se pueden observar los rayos trazados cuyas trayectorias están dominadas por la refracción que experimentan en un medio con variación lineal en la velocidad de propagación de la onda de la forma v = 1800+0.9z, usando 1 fuente (3a) y una resolución de grid de 100x50 celdas, y 5 fuentes (5b) con 240 rayos por cada fuente en un modelo con resolución de 360x160 celdas.

En la Figura 3 se pueden observar los rayos trazados cuyas trayectorias están dominadas por la refracción que experimentan en un medio con variación lineal en la velocidad de propagación de la onda de la forma v = 1800+0.9z, usando 1 fuente (3a) y una resolución de grid de 100x50 celdas, y 5 fuentes (5b) con 240 rayos por cada fuente en un modelo con resolución de 360x160 celdas.

Figura 3 Trazado de rayos en un medio con gradiente lineal. (a) una fuente y 80 rayos en un grid de (100x50). (b) 5 fuentes y 240 rayos en un grid de (320x160). 

Reflexiones en el fondo

Como se describió previamente, en este modelado es posible introducir reflexiones sobre una interfaz plana en el "fondo" del modelo inspirado en un proceso análogo al del software de la empresa NORSAR, en el cual se dispara rayos desde la fuente hacia un reflector ubicado en el fondo de la región de interés. A partir de una relación trigonométrica que permite asegurar el cumplimiento de las leyes de reflexión, es posible determinar el ángulo de incidencia de un rayo disparado desde la fuente y con llegada en el reflector. Conocido este ángulo se aplica nuevamente el algoritmo de Dijk-stra, donde la nueva fuente es el reflector (basado en el principio de Huygens) y realizando una búsqueda sobre los receptores del modelo, para determinar en cuál de ellos el ángulo de incidencia desde la fuente coincide con el ángulo de salida del reflector hacia uno de los receptores; cumpliéndose así la Ley de reflexión.

Uno de los problemas de la discretización en el método del camino más corto es que los ángulos barridos por las trayectorias de los rayos comprenden un conjunto finito de valores no muy fino, por ende, la búsqueda de un reflector cuyo ángulo de incidencia desde la fuente sea igual al ángulo reflejado hacia un receptor no es estrictamente satisfecha. Para resolver este problema, se usó un criterio en el cual se acepta un punto como un reflector si la diferencia entre el ángulo de incidencia y el reflejado es menor que 0.05rad « 2.8°. Esta condición puede variar significativamente según el radio r de propagación de los rayos, ya que para radios pequeños se barren menor cantidad de ángulos entre 0 rad y n/2 rad. La Figura 4a muestra un medio con refracciones y reflexiones en el fondo simultáneas.

Reflexiones flotantes

Como se mencionó en la sección anterior, la estrategia de modelado que se propone en este trabajo permite ubicar en el modelo lugares flotantes en los cuales se da la reflexión de los rayos, esto con el fin de aproximar de mejor manera la reproducción de características estructurales localizadas y con potenciales fuertes de variaciones laterales. La Figura 4b muestra el resultado de dichas reflexiones en un par de regiones planas horizontales que pueden representar por ejemplo los bordes superior e inferior de una capa de contraste o un elemento difractor. La oportunidad de identificar tales zonas en el trazado de rayos y correspondiente inversión, representa una posibilidad valiosa para identificar potenciales contrastes de impedancia en el medio durante una tomografía.

Figura 4 Trazado de rayos. a) Refracciones y reflexiones simultáneas en un medio con variación lineal en profundidad. b) Refracciones y reflexiones flotantes en nodos puntos específicos del modelo 

3.2. Resultados tomografía

Aunque el trazado de rayos per-se es un resultado interesante que resuelve el problema de modelado en un proceso de tomografía o inversión sísmica, su uso toma valor precisamente cuando se usa en conjunto con un procedimiento de inversión. Dado que el objetivo principal de este trabajo es explorar la estrategia de modelado más que la inversión, concentramos toda nuestra atención en discutir el uso del modelador en situaciones sintéticas simples con los resultados obtenidos para estimar velocidades usando la técnica de reconstrucción algebraica. Vale la pena anotar que una aproximación como estas no se encuentra reportada en la literatura y resulta interesante explorar dado que el comportamiento robusto ofrecido por un método de reconstrucción algebraica puede ofrecer ventajas a la hora de enfrentar situaciones afectadas por discretización, como es el método del camino más corto. Además, se espera que, si la estrategia de modelado funciona bien con una estrategia de inversión simple como esta, con certeza deberá funcionar en escenarios mucho más complejos como los ofrecidos por estrategias de inversión con gradiente. Para esto entonces primero se muestran los resultados del método en la solución de sistemas de ecuaciones lineales 3.2.1, común en los procesos de tomografía sísmica, y segundo los resultados obtenidos en inversión basada en modelos 3.2.2.

3.2.1 Solución de matrices sparse usando técnicas de reconstrucción algebraica

El objetivo de estos experimentos es verificar la efectividad del método de inversión para resolver el sistema de ecuaciones lineales DS = T. En la Figura 5 se pueden ver los resultados de la inversión, donde se trazaron rayos por refracción (como en la Figura 3) en el modelo sintético de la Figura 5a y se hallaron los tiempos de propagación T y la matriz de distancias D.

Figura 5 Recuperación de un modelo con gradiente lineal a partir de un modelo con slowness cero. Solo rayos refractados 

Para resolver el sistema DS = T, se usa como modelo inicial el vector S (Figura 5c) con entradas 0.0 en todas sus componentes. Esto significa que c -»oo, lo cual no tiene sentido físico, pero muestra la solidez del método para encontrar la solución a pesar de tomar dicha condición inicial. El modelo S se va actualizando a partir de la Ecuación (10) y dejando fijos la matriz D y el vector T. En la parte inferior y límites del trazado de rayos se puede observar los efectos en la reconstrucción del modelo sintético debidos a la poca iluminación sobre las celdas. A mayor cantidad de rayos atravesando una celda mayor será la velocidad de convergencia del método.

Modelo con gradiente lineal

La Figura 6 muestra la inversión en un medio donde la velocidad varía linealmente con la profundidad. En la Figura 6b, se muestra los resultados obtenidos usando inversión basada en modelos. El modelo sintético tiene una velocidad determinada por la función lineal v = 1800 + 1.1z y fue usado para encontrar los tiempos observados T ^obs . Para la tomografía se eligió un modelo inicial con velocidad v = 1800 + 1.4z, este modelo es diferente en todas las capas al modelo sintético, a excepción de la primera capa (z = 0.0m) donde la velocidad es de 1800 m/s para ambos.

Figura 6 Recuperación de un modelo con gradiente lineal partiendo de un modelo con igual cantidad de capas y un gradiente mayor 

La mejor forma de interpretar los resultados obtenidos, es a partir de la gráfica Figura 6d, que muestra la diferencia entre el modelo sintético (6a) y la tomografía final (6b). Allí se puede apreciar que en la región por donde pasan los rayos refractados las celdas del modelo se acercan al modelo inicial (entre más oscuro es el azul, mejor es la solución de tomografía).

Modelos con gradiente lineal y anomalía rectangular

En este experimento se tomó un modelo sintético con gradiente en profundidad dado por v = 1800+0.2 z y una anomalía rectangular (difractor) ubicada a mitad de profundidad con una velocidad dada por v = 1600+0.1z. Para este modelo, se trazaron rayos por refracción y se combinaron con reflexiones flotantes en las caras superior e inferior de la anomalía rectangular, como se muestra en la Figura 7. El modelo inicial tomado fue uno con gradiente lineal dado por v = 1800+0.2z.

Para este experimento las reflexiones flotantes se ubicaron en posiciones correspondientes a las caras superior e inferior de la zona del difractor, de forma análoga a como se muestra en la Figura 4b. Es interesante ver como solo usando la información de las posiciones de las reflexiones flotantes y un modelo inicial solo con un gradiente en profundidad, nuestra estrategia de trazado de rayos+tomografía es capaz de recuperar información en profundidad del obstáculo difractor. Hay que aclarar aquí que las posiciones de las reflexiones flotante no deben coincidir con la posición del difractor, basta con que los rayos reflejados pasen a través del obstáculo para que el mecanismo de tomografía lo identifique.