TRABAJOS DE INVESTIGACIÓN

Diseño e implementación de un algoritmo trazador de rayos en un dominio circular para tomografía sísmica 2D

Design and implementation of a 2D circular domain seismic tomography ray tracer algorithm

 

Carlos Adolfo Calvo¹, Armando Luis Imhof², Adriana Martin³

¹ Departamento Matemática Aplicada – FI – UNSJ – San Juan – Argentina.
² Instituto Geofísico S. Volponi – FCEFN – UNSJ – San Juan – Argentina
³ Departamento Informática – FCEFN – UNSJ – San Juan – Argentina.
E-mail: ccalvo@unsj.edu.ar

---

RESUMEN

Existen numerosos métodos de inversión por la técnica de tomografía en tiempo de viaje. Varios de ellos se basan en métodos matriciales que requieren la construcción de matrices de distancias que contienen la información de la geometría del dispositivo (ubicación de sensores: emisores y receptores), así como de la forma del dominio discretizado y la distancia que recorren los rayos en el mismo en cada uno de sus elementos de área componentes. Los dominios considerados por lo general son rectangulares, primero por la facilidad de resolución geométrica, y además por el uso extensivo del arreglo cross-hole con los emisores y receptores ubicados vertical y paralelamente en sendas perforaciones de exploración petrolera. En este trabajo se desarrolló un procedimiento original para la construcción de estas matrices, necesarias para resolver el sistema lineal, que lleva a la determinación de los valores de velocidad o lentitud en cada uno de los elementos, pero en un dominio circular. El procedimiento se validó utilizando el método de mínima longitud del vector solución (Minimum Length Solution). Se simularon datos sísmicos sintéticos a partir de dos modelos teóricos propuestos, y los resultados obtenidos a partir de la inversión de los mismos demostraron la eficiencia del algoritmo para construir matrices de distancias.

Palabras Clave: Trazador de rayos; Matriz de distancias; Tomografía sísmica; Dominio circular; Pixeles no uniformes.

ABSTRACT

There are several investigation methods concerning travel time tomography. Some are based on matrix methods requiring the construction of distance matrices that contain geometric information of the survey configuration (location of sensors, both sources and receivers), the shape of the discretized domain, and the distance the rays travel for each individual area element. The domains are generally rectangular for ease of geometric resolution and also, in petroleum exploration perforations, for the extended use of cross-hole arrays with the sources and receivers located vertically and parallel. In this work an original procedure is developed for the construction of these matrices, which lead to linear systems for the determination of velocity in each of the elements of a circular domain. The procedure was validated using the method of Minimum Length Solution. Synthetic seismic data were simulated from two theoretical models, and the results from the inversion of the data demonstrate the efficacy of the algorithm to constitute distance matrices.

Keywords: Ray tracer; Distance matrices; Seismic tomography; Circular domain; Non uniform pixels.

---

 

INTRODUCCIÓN

Desde hace años se investigan debilidades y/o fracturas en diversos tipos de estructuras utilizando auscultaciones geofísicas a tal fin, por ejemplo Kanli, 2008; Binda et al., 2003; Binda et al., 2001. Estos métodos indirectos se basan en determinar propiedades físicas específicas del medio investigado, cuyas variaciones a su vez se relacionan con fallas o anomalías en dichas estructuras. Diferentes metodologías se aplican a fin de abordar el problema desde diferentes ángulos (i.e.se estudian diferentes tipos de propiedades físicas con instrumentos específicamente adaptados para cada caso). Una de ellas, extensamente utilizada es la tomografía sísmica, por ejemplo en Gheslaghi et al., 1995 y Santamarina, 1994. La misma se basa en colocar en un sector de la periferia del dominio estudiado emisores de ondas mecánicas; y en otro receptores a fin de recibir y almacenar las ondas provenientes de aquellos. Por lo general cada receptor almacenará por vez los eventos provocados por cada uno de los emisores (Imhof, 2007).
El problema a resolver luego, se divide en dos partes: a) Pre-procesar las señales a fin de establecer tendencias en los datos, depurar los mismos y llevar a cabo una primer interpretación cualitativa; b) Proceder a la inversión de los datos, a fin de poder calcular los valores de parámetros interiores del dominio a partir de las mediciones efectuadas en el contorno, que en tomografía sísmica en tiempo de propagación por lo general es la lentitud, inversa de la velocidad.
Para poder realizar la inversión existen varias alternativas, comenzando por plantear y establecer hipótesis sobre el dominio que contiene la/s anomalía/s que se desea investigar. La gran pregunta es ¿Qué tipo de teoría justifica el medio estudiado?.
Como parte de la hipótesis de partida, se debe tener en cuenta la longitud de onda de la señal de emisión (que dependerá también del tipo de material atravesado) en relación con el tamaño de la anomalía; si es posible utilizar la hipótesis del rayo y, finalmente entre estos dos aspectos, si el rayo considerado puede ser recto o curvo.
En este trabajo se considera que el tamaño de la anomalía es suficientemente mayor al de la longitud de onda de la señal, por lo que se descarta fenómenos de difracción (Santamarina y Fratta, 2005). Se considera al dominio como medio homogéneo e isotrópico, por lo que las ondas sísmicas se asumen como trayectorias rectilíneas, y se utiliza la hipótesis del rayo recto, linearizando el problema (Tarantola, 1987). Por último se consideran que los rayos recorren distancias pequeñas en las anomalías, por lo que no se consideran los fenómenos de refracción en las interfases. En muchos casos normales esto no produce un error perceptible (Tarantola, 2005).
Se determinan los tiempos de primeros arribos t de los rayos sísmicos; la alteración de los mismos respecto del normal es un indicador de la presencia de anomalía/s. Resumiendo, el planteo general de la tomografía es modelar (i.e. determinar forma y ubicación) las mismas (Menke, 1989).
De los numerosos planteos generales de solución de este problema de inversión, varios discretizan el dominio en pequeñas celdas rectangulares (píxeles) planteando las ecuaciones de tiempo de viaje recorridos por los rayos, siendo las incógnitas la velocidad en cada píxel. Se llega a un sistema de ecuaciones lineales M.s = t , donde t es el vector de tiempos medidos, M la matriz de distancias (construida a partir de un algoritmo trazador de rayos) y s es el vector de lentitudes (inversa de la velocidad). Normalmente en lo relativo a software de modelado geofísico, se han implementado algoritmos trazadores de rayos con el objetivo de construir matrices de distancias para geometrías rectangulares, debido a su extensiva utilización en perforaciones petroleras verticales enfrentadas (dispositivo cross-hole). Sin embargo no se pueden utilizar estos trazadores para auscultar columnas y vigas de geometría diferente. Además los mismos emplean pixeles de igual tamaño. La construcción de un algoritmo trazador de rayos para dominios circulares y con pixeles de tamaños generales es el objeto de este trabajo. Imhof et al. (2010) demostraron que aumentar la densidad de pixeles por sectores del dominio de tal forma de lograr iguales longitudes de propagación en todas las celdas (i.e. contenido de información) mejora el condicionamiento de la matriz M, tornando más robusta la inversión. Esto justifica ampliamente la división del dominio en pixeles de diferente tamaño.
Por último, se pueden proyectar varias alternativas en la división del dominio; dos de ellas ya han sido publicadas (Imhof y Calvo, 2014a; Imhof y Calvo, 2014b); en este trabajo se decidió dividir un cuadrante en tamaños generales de celdas y luego, por simetría trasladar esa división a los otros tres.

DEFINICIÓN DEL PROBLEMA

Para un modelo sintético, se construye un dominio circular (Figura 1). Este dominio se discretiza con líneas verticales y horizontales que forman celdas (pixeles) rectangulares (excepto en los bordes). Éstas celdas son de dimensión arbitraria en cada cuadrante, pero simétricas según los ejes coordenados en los otros tres. Sobre el borde se colocan los emisores ei y los receptores rj en forma libre cuyos rayos cortan a la grilla formando trazos en las celdas interceptadas que constituyen los elementos de la matriz de distancias

Figura 1. Dominio circular con 3 emisores (e1, e2 y e3), 2 receptores (r1 y r2) y 6 rayos rectos.
Figure 1. Circular Domain with 3 sources (e1, e2 y e3), 2 receivers (r1 and r2) and 6 straight rays

A cada fila de la matriz M le corresponde un rayo k. Y a cada columna tiene asignada una celda c. En su camino
el rayo k-ésimo al atravesar la celda c, recorre una distancia dkc siendo esta el elemento de la matriz M. Para
aquellas celdas en que el rayo no las toca dkc = 0.
Los datos del problema son los parámetros de la grilla y la ubicación de los emisores y receptores.

CONSTRUCCIÓN DE LA GRILLA

Se introducen N valores positivos de ordenadas crecientes, con pasos variables; siendo radio del dominio (Figura 2).

Figura 2. Grilla N=3. En trazo continuo yg = [3 4 5]. En línea de puntos la grilla completa de 12 celdas (pixeles)
Figure 2. N=3 grid. In solid line yg = [3 4 5]. Complete 12 grid cells (pixels) in dashed line.

Se obtienen las respectivas abscisas
Se completan con cero y valores negativos (ya que todo se construye en forma simétrica a un cuadrante)
resultando el vector de ordenadas y el vector de abscisas:

La grilla se construye con líneas horizontales de ordenadas Yg y líneas verticales de abscisas Xg; resultando 2 celdas en la primera fila (fila superior), 4 en la segunda, 6 en la tercera y así creciendo en incrementos de 2 hasta llegar a 2N en la fila N. Por simetría se completa la mitad inferior de la grilla.
La grilla obtenida es simétrica respecto de ambos ejes con 2N(N+1) celdas.

Numeración de la grilla
La grilla se numera por fila desde arriba hacia abajo c=1:2N(N+1) (Figura 3).

Figura 3. Celdas numeradas de 1 a 24. Los pares (i,j) ubican las celdas con 2 coordenadas.
Figure 3. Numbered cells from 1 to 24. Pairs (i,j) locate them with 2 coordinates.

Es importante determinar el número de la celda en función de la ubicación para relacionarla con el camino del rayo. Para ello se introduce los ejes (i,j) que ubican filas y columnas de la grilla, donde i=1:2N y j=1:2N .
Para obtener la relación c=f(i,j), primero se observa que sobre una fila, c se obtiene sumando o restando el valor de i a partir de una columna fija (i.e. i=N), esto es:

En segundo lugar, se observa que sobre la columna i=N, j=2N:-1:N+1 (grilla superior), c toma los valores 1 ,4 , 9 ,…N2, para j=2N,2N-1,2N-2…N+1.
Resulta:

Figura 4. La trayectoria del rayo k-ésimo interseca la grilla en 8 puntos determinando 7 elementos dkc de la matriz M.
Figure 4. k-th raypath intersects grid in 8 points determining 7 elements dkc from matrix M.

Numeración y camino de los rayos
Tanto los n emisores como los m receptores están ubicados en forma genérica sobre el contorno de un dominio circular y su número es arbitrario. Cada par (ei,rj) define un rayo k ; estos se numeran desde 1 a n.m fijando por vez un emisor y recorriendo todos los índices de los receptores resultando:

FORMACIÓN DE LA MATRIZ DE DISTANCIA

VALIDACIÓN DEL ALGORITMO

Problemas de Inversión por Tomografía
Se diseñó un par de modelos sintéticos, con distribución de sensores ei y rj; con rayos rectos normales respectivos, representados en Figura 5-a.

Figura 5. (a) Distribución de emisores ei y receptores rj y rayos en el dominio circular. Radio del dominio R=0.3m. (b) Discretización del dominio circular en pixeles. nc=2N(N+1)=2x28x29=1624 elementos.
Figure 5. (a) Sources ei and receivers rj distribution in circular domain. Domain Radius R=0.3m; b) Pixel circular domain discretization. nc=2N(N+1)=2x28x29=1624 elements.

Con los datos definidos del dominio base se procedió a la construcción de la matriz de distancias con el algoritmo trazador generado y de acuerdo a una grilla de 1624 pixeles rectangulares irregulares, según se aprecia en la Figura 5-b.
Una vez definidos el número de sensores y la densidad de la grilla de elementos, se diseñaron 2 modelos con anomalías de diverso tipo y tamaño (Figura 6 y Tabla 1), construyendo la matriz de datos para cada uno de ellos con un programa de modelado implementado para tal fin (Figura 7 a-b).

Figura 6. Representación gráfica de los modelos sintéticos. Columna cilíndrica de radio R=0.3m. Sensores indicados con círculos en los bordes. (a) Modelo A. (b) Modelo B (Tabla 1).
Figure 6. Synthetic models representation. Cylindrical column with radius R=0.3m. Edge sensors located as small circles. (a) Model A; (b) Model B. (Table 1),

Figura 7. Modelos A y B con ubicación de anomalías. Compárese los resultados con modelos teóricos de Tabla 1.
Figure 7. Models A and B with anomalies location. Compare results with theoretical models from Table 1.

Tabla 1. Dos Modelos diferentes para validar la inversión (Mod. A y Mod. B). Los ítems ai denotan anomalía i-ésima, con valores de parámetros indicados.
Table 1. Two distinct models to validate inversion (Mod. A and Mod. B). Items ai means i-th anomaly, with indication of parameter’s values.

Inversión
Los métodos de tomografía sísmica que utilizan como hipótesis restrictiva la utilización de rayos rectos que implican la determinación de valores de lentitud/velocidad en dominios discretizados en forma de pixeles, conllevan la resolución de sistemas lineales que necesariamente deben ser sub-determinados, ya que cada pixel es una incógnita a resolver, según la forma:

Donde Mn.m,nc es la matriz de distancias, xn.m,1 el vector de lentitudes/velocidades y tn.m el vector de tiempos de llegada de los eventos sísmicos.

Para la grilla seleccionada n=m=17; por lo que el número de filas de M es 289; y el número de columnas (i.e. cantidad de de celdas), nc=2380 pixeles. Esto representa un sistema lineal muy sub-determinado.
El problema es el siguiente: para lograr un sistema lineal al menos equi-determinado (ni qué decir sobredeterminado) que posibilite la aplicación del método de mínimos cuadrados y a la vez obtener imágenes lo más resolutivas posibles (lo que implica reducir el tamaño de los pixeles), será imprescindible incrementar enormemente el número de sensores de tal forma de conseguir igual número de filas (rayos sísmicos; nxm) que columnas (número de pixeles) de M. Esto desde luego es inviable desde el punto de vista práctico, tanto por el costo de instrumental, el volumen de trabajo necesario y la falta de espacio para colocar los sensores.
Una de las opciones (no la única) es aplicar el método MLS (eng. Minimum Length Solution) que se utiliza en sistemas sub-determinados (Santamarina y Fratta, 2005) cuya expresión es:

Donde el término h-g es la inversa generalizada de Penrose (Penrose, 1955). Es conveniente añadir alguna información adicional a fin de favorecer la convergencia. Usualmente se estima una lentitud/velocidad para el medio base. La expresión queda:

Donde x0nc,1 es el vector de estimación inicial de lentitud/velocidad.
Los resultados de la inversión se muestran en la Figuras 7 a-b.

CONCLUSIONES

Modelar sísmicamente por métodos matriciales dominios diferentes al rectangular resulta complejo por tener que distribuir pixeles rectangulares en dominios curvos. Pero en muchos casos se necesitan este tipo de algoritmos a fin de resolver situaciones reales.
El algoritmo trazador de rayos sísmicos rectos diseñado en este trabajo permite obtener matrices de distancias para geometrías circulares, utilizando celdas de tamaños diferentes, lo que constituye una nueva herramienta computacional para resolver la inversión de datos de tomografía sísmica en columnas cilíndricas mediante métodos matriciales.
Además el algoritmo posibilita una forma particular de dividir el dominio en pixeles de diferente tamaño, para poder de este modo estudiar la relación entre el número de condición y la distribución de sensores en la periferia y de pixeles en el interior del dominio, lo cual abre nuevas posibilidades de investigación ya que permite estudiar la relación entre el número de condición, el tamaño relativo de los pixeles y la geometría de los sensores.
La validación llevada a cabo demostró el adecuado funcionamiento del algoritmo, resultando la inversión en la localización precisa de las anomalías existentes en el dominio.

Agradecimientos

El presente trabajo de investigación forma parte de proyectos subvencionados por la Secretaría de Ciencia y Técnica (CICITCA) de la Universidad Nacional de San Juan, República Argentina.

REFERENCIAS

1. Binda, L., A. Saisi, C. Tiraboschi, (2001). Application of sonic tests to the diagnosis of damaged and repaired structures. NDT&E International, 34: 123-138.         [ Links ]

2. Binda, L., A. Saisi, L. Zanzi, (2003). Sonic tomography and flat-jack tests as complementary investigation procedures for the stone pillars of the temple of S. Nicolò l’Arena (Italy); NDT&E International, 36: 215-227.

3. Gheslaghi, F., J.C. Santamarina, D. Wiese, M. Thomas, M. Polak, G. Caratin, (1995). Tomographic imaging concrete structures. Proceedings of the International Symposium on NDT in Civil Engineering (NDT-CE). Ed. G. Schickert and H. Wiggenhauser, 1: 297-302.         [ Links ]

4. Imhof, A.L., (2007). Caracterización de arenas y gravas con ondas elásticas: tomografía sísmica en cross-hole". Tesis Doctoral. Universidad Nacional de Cuyo. Mendoza. República Argentina.         [ Links ]

5. Imhof, A.L., C. Calvo, (2014). Diseño e implementación de un algoritmo trazador de rayos general en un dominio circular para tomografía sísmica 2D. Revista Ingeniería y Desarrollo. Universidad del Norte, 32(2): 242-259. Colombia.         [ Links ]

6. Imhof, A.L., C. Calvo, (2014). Diseño y Validación de un algoritmo trazador de rayos para celdas uniformes en un dominio circular para tomografía sísmica 2D. Revista Investigación Operacional. La Habana, Cuba. 35(1): 78- 87.         [ Links ]

7. Imhof, A.L., C.A. Calvo, J.C. Santamarina, (2010). Seismic data inversion by cross-hole tomography using geometrically uniform spatial coverage. Revista Brasileira de Geofísica, 28(1): 79-88.         [ Links ]

8. Kanli, A.I., (2008). Image reconstruction in seismic and medical tomography, J. of Engineering and Environmental Geophysics (JEEG), 13(2): 85-97.         [ Links ]

9. Menke, W., (1989). Geophysical Data Analysis: Discrete Inverse Theory. Ed. Academic Press Inc., New York, USA.         [ Links ]

10. Penrose, R.A., (1955). A Generalized Inverse for Matrices, Proceedings of Cambridge Phil. Society, 51: 406-413.         [ Links ]

11. Santamarina, J.C. y D. Fratta, 2005. Discrete signals and inverse problems: an introduction for engineers and scientists, first ed., Ed. Wiley.         [ Links ]

12. Santamarina, J.C., (1994). An introduction to geotomography, in Geophysical Characterization of Sites, R.D. Woods (Eds), New Hampshire: International Science Publisher, pp.35-44.         [ Links ]

13. Tarantola, A., (1987). Inverse problem theory, Ed. Elsevier, Amsterdam.         [ Links ]

14. Tarantola, A., (2005). Inverse Problem Theory and Model Parameter Estimation, Ed. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, USA, 2005.         [ Links ]

Recibido: 06-04-2015
Aceptado: 15-07-2015