1. Introducción
La exploración sísmica busca generar imágenes del subsuelo para estudiar sus características estructurales a través del análisis de la propagación de ondas sísmicas. Esta técnica es de uso común en situaciones como la exploración petrolera, minería, exploración de acuíferos, obras de ingeniería civil, etc. (Hernández et al., 2006). Específicamente, en áreas estructuralmente complejas, las imágenes en escala de profundidad revelan las características del subsuelo, lo que convierte la técnica de imágenes sísmicas en una herramienta bastante utilizada por las compañías petroleras para la ubicación de nuevas reservas en regiones con geología compleja.
Existen diferentes mecanismos por medio de los cuales se puede hacer exploración sísmica. Una de esas aproximaciones es la tomografía sísmica. La tomografía es la obtención de imágenes de cortes o secciones de algún objeto a partir de información recolectada en su superficie y hace parte de un conjunto de problemas denominado en matemáticas como Problema Inverso, el cual consiste en la obtención de un conjunto de parámetros de un modelo a través de la interpretación de datos observados que en su mayoría son inexactos y poseen ruido (Jones, 2010).
En el transcurso del siglo XX, episodios notables marcaron avances en los métodos de prospección sísmica. Si bien muchas de estas tecnologías tardaron tiempo para pasar del formalismo a la práctica, cada una generó finalmente nuevas oportunidades de exploración. En 1920, se introdujeron los disparos analógicos de cobertura simple para detectar capas inclinadas del subsuelo (Albertin et al., 2002), mostrando esto la necesidad de conocer una aproximación del terreno antes de perforar el suelo. Desde esta época, el problema de encontrar yacimientos de hidrocarburos se ha tornado cada vez más complejo y los costos de realizar excavaciones en busca de estos son altos. Por ende, diferentes técnicas computacionales de tomografía se han utilizado para tener una vista previa del subsuelo a explorar y mejorar los índices de éxito en la detección de hidrocarburos. El trazado de rayos y la tomografía de tiempos de propagación son utilizados para este propósito. El uso de técnicas como el trazado de rayos para aproximar el fenómeno de la propagación de ondas en el subsuelo se puede atacar de diferentes maneras. La solución de la ecuación Eikonal por separación de variables y su correspondiente solución numérica usando el método de diferencias finitas (Vidale, 1988) ha resultado muy común. Una alternativa menos restrictiva al método de solución directa, es el trazado de rayos en un grid usando la técnica del camino más corto (Moser, 1991). Estos métodos de trazados de rayos son usados en el modelado sísmico, para encontrar los tiempos de propagación de la onda.
Las primeras referencias relacionadas al trazado de rayos en un grid son de (Nakanishi y Yamaguchi, 1986) quienes proponen un trazador de rayos en su estudio de inversión no lineal de los primeros tiempos de llegada, para encontrar el camino del rayo y los tiempos de propagación de un punto fuente a un número de puntos receptores. (Moser, 1991, 1992) utilizó esta idea e investigó métodos para mejorar la eficiencia de la búsqueda de la trayectoria de un rayo en el tiempo mínimo y además mostró cómo es posible tratar los tiempos de llegada de los rayos reflejados. Su idea consistía de un grid de puntos que son conectados con sus puntos vecinos por una arista y la longitud de dicha arista es tan larga como el tiempo que tarda la onda en ir de un punto a otro. La idea de Moser fue brillante, pero la implemen-tación computacional no era muy eficiente. Durante los últimos años, la tomografía basada en rayos ha tenido gran evolución gracias a las nuevas herramientas de cómputo que día a día procesan gran cantidad de datos en un tiempo cada vez menor (Fajardo y Castillo, 2013; Meléndez et al., 2013; Malony, McCumsey, et al., 2017; Zhang et al., 2017, 2018; Monil et al., 2018).
Uno de los aspectos computacionales más complejos en el problema de trazado de rayos a través de la aproximación del camino más corto es la identificación de los segmentos de línea que dan forma a cada rayo. En teoría de grafos, un algoritmo muy eficiente para determinar el camino más corto en un grid de puntos fue propuesto por (Dijkstra, 1959). Varias modificaciones de este algoritmo han sido usadas generalmente en el trazado de rayos sobre grids. (Gallo y Pallottino, 1986) describen versiones modernas del algoritmo de Dijkstra que son lo suficientemente eficaces para hacer que en la actualidad el método del camino más corto se vuelva competitivo con otros métodos para el modelado sísmico.
En aras de mejorar cada vez más los tiempos de cómputo, diferentes mejoras a la técnica del camino más corto fueron propuestas por: (Fischer y Lees, 1993; Cheng y House, 1996; Zhang, Chen y Xu, 2004; Zhang et al., 2013). La aparición de herramientas de cómputo mucho más poderosa impulsó el desarrollo del trazado de rayos ofreciendo la facilidad de trabajar con grids mucho más grandes y obtener resultados en menor tiempo. De esta manera los algoritmos proporcionan en la actualidad un marco atractivo para la tomografía basada en rayos, algunos de estos algoritmos se encuentran trabajos recientes (Piçta y Dwornik, 2012; Giroux y Larouche, 2013; Malony, McCumsey, et al., 2017; Malony, Monil, et al., 2017; Monil et al., 2018), lo que hace a la técnica del trazado de rayos para tomografía un mecanismo que aún es atractivo para la industria y la academia.
En este trabajo se propone la implementación de un modelador de la propagación de ondas sísmicas en el subsuelo a través de la aproximación de rayos, específicamente, a través del uso de la aproximación del método del camino más corto. El trazador implementado que permite calcular simultáneamente trayectorias de rayos por refracción, reflexión en la última capa del modelo y de reflexiones en puntos de interés en el medio. Estas estrategias han sido usadas individualmente en otros trabajos (Meléndez et al., 2015) y (Zhang et al., 2018). El modelador se usó en el contexto de la tomografía sísmica acoplado a un método de reconstrucción algebraica. Esta es la primera vez que en la literatura se reporta el uso de este tipo de técnicas de reconstrucción acopladas a un modelador de trazado de rayos por medio del método del camino más corto y su aplicación se convierte en un aspecto novedoso en el contexto de la tomografía sísmica.
2. Métodos
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.
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
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.
2.2.4. Múltiples fuentes y reflexiones flotantes
Para mejorar la iluminación del medio con más rayos, se propone realizar varios disparos desde fuentes distintas. La forma como se realiza el trazado de los rayos es absolutamente independiente de la posición de su origen o fin, así que esto hace que solo sea necesario aplicar el algoritmo de Dijkstra tantas veces como fuentes se desee ubicar.
La solución clásica de la ecuación Eikonal convencionalmente falla cuando se tienen medios con fuertes variaciones laterales o anisotropías bien localizadas espacialmente en el medio. Con el fin de atacar este problema es posible generar rayos que se reflejen en algunos puntos específicos, con el fin de mejorar la iluminación en ciertas regiones y de modelar de forma más directa los primeros arribos de las reflexiones en las anomalías o interfaces de interés. Para este fin se escoge, en el grid de nodos por donde se trazan los rayos, el id del nodo en el cual se quiere reflejar un rayo. A estos puntos los llamamos puntos de reflexiones flotantes, ya que, a diferencia de las reflexiones en el fondo del medio, corresponderán a lugares que pueden variar de modelo a modelo. Para trazar los rayos se escoge tanto la fuente como el receptor, equidistantes al nodo que se encuentra en el punto de corte entre la superficie y la perpendicular que pasa por el punto de reflexión. Se aplica el algoritmo de Dijkstra, usando como fuente el reflector elegido, y trazando las trayectorias desde este punto hasta la fuente y el receptor. Al igual que se mencionó anteriormente, a causa de la discretización del medio esta forma de hallar la reflexión no garantiza con precisión la ley de reflexión en el reflector, pero es una variante para mejorar la iluminación y encontrar una trayectoria mínima entre una fuente y receptor fijo para cualquier otro punto de interés en el medio.
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 G0, 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 G0, G v G2, ... 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. Resultados
Esta sección presenta los resultados del trazador de rayos implementado usando datos de refracción y reflexión, además del uso del trazador en la obtención de imágenes tomográficas del subsuelo a partir de datos sintéticos.
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.
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.
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.
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.
3.2.2. Tomografía SIRT: Inversión basada en modelos.
La Figura 5 corresponde a las soluciones de un sistema de ecuaciones lineales sparse y mal condicionado de tamaño m x n, en donde la matriz de distancias D no cambia durante las iteraciones y en la cual se verificó que las actualizaciones realizadas al modelo inicial convergen al modelo verdadero. En los problemas reales de tomografía sísmica se tiene que DS = T es decir, la matriz de distancias depende del modelo inicial elegido y por lo tanto debe cambiar en cada actualización del modelo de velocidades. En este caso, el modelo inicial no puede ser un vector de ceros como en el caso anterior, pues los rayos se propagarían todos por la superficie del medio desde las fuentes hasta los receptores. Por esta razón y para garantizar una iluminación similar a la del modelo sintético, el modelo inicial debe tener una configuración cercana al modelo verdadero (una restricción que no es exclusiva de los problemas de tomografía de tiempo de propagación, sino asociada a cualquier problema de inversión). El objetivo es encontrar un modelo que minimice la diferencia entre los tiempos de propagación calculados en el modelo y los tiempos de propagación observados en la superficie.
Una simple iteración del método de SIRT consiste en realizar una actualización al modelo inicial escogido. Luego de la primera actualización del modelo, los rayos son trazados nuevamente usando el modelo ajustado, y se determinan las nuevas trayectorias de los rayos con la misma configuración de fuentes y receptores. Esto implica que el operador de distancias D debe ser calculado nuevamente según las nuevas trayectorias en el modelo actualizado. Este proceso se denomina inversión basada en modelos y se repitió hasta que la norma euclidiana entre el tiempo observado T obs y el tiempo calculado en cada iteración T cal fuera menor a 0.001s.
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.
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.
4. Discusión y conclusiones
En este trabajo se ha presentado una estrategia para el modelado de propagación de ondas en el subsuelo a través de la aproximación del método del camino más corto con un enfoque particularmente útil para la tomografía sísmica de tiempos de propagación. Nuestro modelador considera los efectos sobre la trayectoria del rayo inducidos por refracción, reflexión en el fondo del modelo y reflexión en puntos flotantes en el interior del modelo.
La estrategia usa una discretización de dos niveles que permite reconstruir de manera fina la trayectoria del rayo mientras que a su vez permite recolectar información suficiente sobre la iluminación del medio en cada celda para los propósitos de la inversión necesaria para la tomografía. La aproximación propuesta además hace uso del algoritmo de Dijkstra para adelantar el trazado de los rayos de manera eficiente, lo que hace del algoritmo una opción práctica desde el punto de vista computacional.
Es ciertamente un hecho que las soluciones que se consiguen a través de la tomografía de tiempos de propagación como las propuestas en este trabajo, especialmente de primeros arribos, son restringidas en términos de la cantidad de información que se puede recuperar del subsuelo. Sin embargo, sigue siendo una aproximación a primer orden extremadamente valiosa para hacer primeras exploraciones de la estructura del subsuelo. Más importante aún, resultan ser muy valiosas para proveer la información requerida en la construcción de modelos iniciales para estrategias mucho más complejas, pero más sensibles a las características del modelo inicial escogido, como por ejemplo migración de tiempo reverso (RTM) o inversión de onda completa (FWI), en los que la escogencia del modelo inicial es un ingrediente crucial.
Sigue habiendo limitaciones asociadas con la estrategia de trazado de rayos, y que son aún limitaciones en el modelador que se presenta en este trabajo. Para propósitos de tomografía, se tienen limitaciones en términos de la resolución espacial implicada, toda vez que la discretización de las celdas de tomografía implican la necesidad de usar celdas suficientemente grandes de tal manera que se pueda asegurar una iluminación por celda para que la inversión sea robusta. Este efecto hace que el modelo sufra de pixelación y tenga una resolución espacial limitada.
Una novedad de este trabajo es la combinación de la estrategia de trazado de rayos con el método del camino más corto acoplado a una estrategia de inversión por medio de una técnica de reconstrucción algebraica. Se ha mostrado que el método de reconstrucción algebraica se acopla bien al esquema discretizado del método del camino más corto, ofreciendo resultados satisfactorios en lo que tiene que ver con la identificación de estructura en el subsuelo. Además, el uso combinado de estas técnicas resulta ser computacionalmente económico, lo que las hace aún bastante atractivas para experimentos de exploración de modelos, en los que se requiere probar muchas realizaciones de diferentes modelos o diferentes modelos iniciales