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

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

Armando Luis Imhof*
aimhof@unsj.edu.ar

*Doctor en Ingeniería. Instituto Geofísico Sismológico Volponi, Facultad de Ciencias Exactas, Físicas y Naturales, Universidad Nacional de San Juan (UNSJ).

Carlos Adolfo Calvo**
ccalvo@unsj.edu.ar

**Magíster en Matemática Aplicada. Departamento de Matemática Aplicada, Facultad de Ingeniería Universidad Nacional de San Juan (UNSJ).

Universidad Nacional de San Juan (UNSJ) Argentina

Correspondencia: Armando Luis Imhof. Dirección: Ruta 12, km 17, Frente "Jardín de los Poetas", 5413-Rivadavia, Provincia San Juan - República Argentina. Teléfono: 11-54-264-4945015

Resumen

Las técnicas matriciales de inversión tomográfica en tiempo de viaje se basan en la construcción de matrices de distancias que contienen la información de la geometría del dispositivo, 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 por la facilidad de resolución geométrica. En el caso de este trabajo es circular por la necesidad de auscultar columnas. Las matrices de distancias se emplean para resolver el sistema lineal, que lleva a la determinación de los valores de velocidad en cada uno de los elementos del dominio. Teniendo en cuenta la discretización, existen varias posibilidades en cuanto a la construcción de las matrices de pixeles, y los mismos pueden ser definidos en un cuadrante y por simetría obtener los restantes; construir matrices con pixeles uniformes con la estrategia anterior, y por último, independizarse del cuadrante y generar cualquier disposición y tamaño individual de los mismos. Cada una de las opciones co permitirá acceder a nuevas posibilidades de investigación. o En este trabajo se desarrolló un procedimiento original para la construc-_o; g ción de matrices con elementos del último tipo.

El procedimiento se validó con datos simulados utilizando el método de SVD, lo que demostró la eficiencia del algoritmo.

Palabras clave: inversión, tomografía sísmica, trazador de rayos, e _d^e matriz de distancias, dominio circular, pixeles no uniformes.

Abstract

Tomographic matrix inversion techniques in travel time are based on the construction of distance matrices containing information of survey array geometry, the discretized domain shape and the distance traveled by rays in each area elements that compose it. The domain considered is generally rectangular for ease of geometric resolution. In this work is circular due to the need to carry out investigation in columns of this shape. Distance matrices are used to solve the linear system, which leads to the determination of velocity values in each element of the domain. Taking in mind discretization, there are several possibilities for the construction of arrays of pixels: may those be defined in one quadrant and obtaining the rest by symmetry; build arrays with uniform pixels using the previous strategy, and finally to be independent from quadrant and generate any distribution and individual size thereof. Each option will provide access to new research possibilities.

In this work an original procedure was developed for the construction of matrices with elements of the last type.

The procedure was validated with simulated data using the SVD method, which demonstrated the efficiency of the algorithm.

Keywords: inversion, ray tracer, distance matrices, seismic tomogra-phy, circular domain, non uniform pixels.

1. INTRODUCCIÓN

En medios homogéneos e isotrópicos se considera que la generación de ondas sísmicas por medios tales como golpes de martillo, explosiones o cristales piezoeléctricos (emisores) y captados en sismógrafos de varios canales (receptores) forman frentes de ondas esféricos cuyos gradientes son rectas, por lo que se puede utilizar la teoría de rayos rectos (straight-ray theory) y considerar que las transmisiones de ondas ocurren de esta forma, siempre que la longitud de onda de la señal sea mucho mayor que el tamaño de las partículas que componen el medio investigado [1]. El tiempo de viaje de los primeros arribos de estos rayos resulta una característica factible de determinar, y la alteración de este tiempo es un indicador de la presencia de anomalías. El planteo general de la tomografía es modelar (determinar forma y ubicación) dicha anomalía. Algunos de los planteos generales de solución discretizan el dominio en pequeñas celdas rectan-guiares (pixeles) planteando las ecuaciones de tiempo de viaje recorridos por los rayos, siendo las incógnitas la velocidad en cada pixel. Se llega a un sistema de ecuaciones lineales M^-s = t, donde t es el vector de tiempos medidos experimentalmente, M la matriz de distancias y s el vector de lentitudes (inversa de las velocidades).

El objetivo de este trabajo es construir un algoritmo para generar matrices de distancia utilizando distribuciones generales (ubicación, tamaño y distribución) de pixeles para dominios circulares.

2. DEFINICIóN DEL PROBLEMA

Para un ensayo de laboratorio 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). Sobre el borde se colocan los emisores (e_i) y los receptores (r) en forma libre, cuyos rayos cortan a la grilla formando trazos en las celdas interceptadas que constituyen los elementos de la matriz de distancia. A cada fila de la matriz le corresponde un rayo k y a cada columna le corresponde una celda c. En su camino, el rayo k al atravesar la celda c recorre una distancia d_w constituyendo este un elemento de la matriz M. Para aquellas celdas que no son tocadas por el rayo, d_ic = 0.

Los datos del problema son los parámetros de la grilla (coordenadas x y y) y la ubicación de los emisores y receptores.

3. CONSTRUCCIóN DE LA GRILLA

El centro del dominio circularde radio R (figura 2) es el centro de coordenadas (x, y). La grilla se genera por la intersección de dos vectores: el de abscisas xg = [-R, x_v x₂, x₃,... R] con n_x=length (xg) divisiones y el de ordenadas yg=[-R, yl, y2, y3,...R]con n_y=length (yg) divisiones. El número de celdas o pixeles será, entonces, de

Cabe destacar que npix representará el número total de pixeles de una grilla rectangular cuyos límites son precisamente [-x^ x_g] y [-y, y_g], que para el caso de figura 2 es de 30. Dentro de esta grilla queda inscripto el dominio circular. Se deberá determinar el número de elementos que quedan dentro del círculo (nc), que para este ejemplo coinciden (nc= npix=30).

La grilla se construye con líneas horizontales de ordenadas y_g y líneas verticales de abscisas x_g; resultando celdas rectangulares que forman una grilla no simétrica (en general) respecto de los ejes coordenados. Esta constituye la diferencia más notable con otro trabajo similar [2].

3.1 Numeración de la grilla

Una de las claves del proceso es la numeración de los elementos. Se parte de una grilla rectangular en la cual está incluido el dominio circular. Por otra parte, se deben numerar solo los pixeles que están dentro del círculo, lo cual hará que en ciertas ocasiones existan elementos fuera del mismo que no se numeran. Cada pixel numerado tendrá asignado un número y una ubicación.

La figura 3 representa un grillado con sus elementos numerados por fila de arriba hacia abajo (c = 1, 2, 3,... nc). Es importante notar que c representa a los elementos dentro del dominio circular.

Se debe determinar el número de la celda c en función de la ubicación para relacionarla con el camino del rayo. Para ello se introduce los ejes que ubican filas y columnas de la grilla. El origen de este sistema está en la esquina superior izquierda, siguiendo la notación habitual de los índices de las matrices.

3.1. a Relación entre los sistemas c = f (i, j)

Para obtener la relación c = f (i, j) se deben obtener los índices i,j extremos (del círculo).

índices extremos jpy jf:

Se debe determinar para cada fila i cuáles son los índices jp y jf correspondientes a la primera y última celda de la grilla.

La fila i comprendida entre la ordenadas y_k-1 e y_k (siendo k = n+l-i ), determina sobre la circunferencia los puntos de abscisas extremas A y B de valores *„ = ±Jj¡^j -yi.jp resulta entonces el máximo número de abscisas de x^que no supera a -x_y jf el máximo número de abscisas de x_g que es inferior a x_f

En notación de MatLab se indica como

Donde la función find(...) da los índices del vector x_g que cumplan la condición dada dentro del paréntesis.

3.1.b Número de celdas y matriz de índices

En cada fila el número de celdas es jf-jp+1; sumando para todas las filas se obtiene el número total de celdas de la grilla circular:

En la matriz de índices, el número de cada fila corresponde al número de la celda c, siendo sus elementos El siguiente algoritmo en lenguaje MatLab la genera:

Esta matriz índice permite pasar del sistema al c, y viceversa.

3.2 Numeración y camino de los rayos

Tanto los n emisores como los m receptores están ubicados en forma libre sobre el contorno de un dominio circular y su número es arbitrario. Cada par (e_i, r) define un rayo k; estos se numeran desde 1 a n.m, fijando un emisor y recorriendo todos los índices del receptor, resultando

Este rayo k (figura 4) recorre un camino recto de ecuación y = ax+b, donde la pendiente a y la ordenada b se determinan por

Los puntos intersecciones del ravo con la grilla tienen por abscisas a Xg (líneas verticales) y a las (provienen de líneas horizontales, por ejemplo, coordenada Xh de P4 en figura 4).

El conjunto de abscisas ampliadas son las de los puntos emisor y receptor {x, x, Xg, Xh} ordenadas de menor a mayor, que restringidas al segmento interno emisor/receptor definen el vector X de abscisas de los puntos P, y el correspondiente vector Y de ordenadas.

Siendo p el número de celdas atravesadas por el rayo c.

Los puntos P, P₂,... P, P_p+1 determinan las distancias d_fc y puntos P_mc (x_m, y_mc) interiores en cada celda, que permiten conocer (i,j) y posteriormente el número c de esta.

Donde la función de MatLab find( ) calcula los índices del vector X_g (o Y_g) que cumplan la condición entre paréntesis.Luego, mediante el empleo de la matriz índices se determina c.

4. FORMACIóN DE LA MATRIZ DE DISTANCIA

La matriz M = {d_fc}tiene tantas filas (n.m) como número de rayos, y tantas columnas como número de celdas nc.

Al sumar por filas:

Se determina la distancia D_k que recorre el rayo k entre el emisor i y el receptor j a través de todas las celdas, siendo d_fc = 0 para aquellas en que el rayo no las toca.

El proceso de construcción de la matriz D se realiza por filas (o sea, fijando el rayo k). Se calcula la pendiente a y la ordenada b; se realiza la intersección (X,Y) con la grilla, se hallan los valores d_k, los puntos interiores p_mc (y_m¿ y_mc)a la celda interceptada, y con estos se determinan los valores c con los cuales se redesigna d_{v c}. Se completa con d_{v c} = 0 para aquellas en que el rayo no las toca. Luego se repite para los demás rayos.

5. VALIDACIóN DEL ALGORITMO

5.1 Modelo para Inversión por Tomografía

Se diseñó un modelo sintético equivalente a uno real, construyendo una planilla de tiempos de primeros arribos a partir de la geometría descrita (figura 5), ubicando sensores distribuidos a lo largo del contorno circular de un tanque lleno con agua (V_base=1400 m/s) en cuyo interior se encontraba una anomalía cilíndrica de concreto (V_anom=3500 m/s) según la descripción. La disposición de los sensores se considera normal al eje del tanque y ubicados en su sector medio. Además se representaron los rayos que parten de cada emisor y arriban a los receptores.

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 con una grilla de 420 pixeles rectangulares irregulares, según se aprecia en la figura 6.

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 lentidud/velocidad en dominios discretizados en forma de pixeles conllevan la resolución de sistemas lineales que necesariamente deben ser subdeterminados [3], [4], ya que cada pixel es una incógnita que se debe resolver, según la forma

M . x = t (6)

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

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 el número de sensores (ei, rj) 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 como por el volumen de trabajo necesario [5].

El método de inversión por SVD es una poderosa herramienta para diagnosticar problemas de inversión y estimar la información disponible. Permite además la determinación de la pseudoinversa [6], [7]. Cualquier matriz genérica M puede descomponerse en 3 matrices ortogonales:

M = U • A • V^T (7)

Y a partir de la ortogonalidad de matrices se puede calcular la pseudoinversa generalizada de Penrose [8]:

M^-g = V^<p> • A^-1 • (U^<P>)^T (8)

por lo que la ecuación general quedará

x^(est) = M^-g • t (9)

Los elementos de la matriz diagonal A en la ecuación. (8) son los valores singulares X de M.M^T en orden descendente. La figura 7 muestra los mismos.

Usualmente en sistemas subdeterminados este método suele converger muy bien; excepto cuando la subdeterminación es muy grande, como es el caso presente; por lo que es necesario añadir alguna información adicional a fin de favorecer la convergencia. Usualmente se estima una velocidad para el medio base. Modificando la ecuación (9), la misma asume la forma

x^(est) = x0 + M^-g • (t - M • x0) (10)

Donde x0_nc1 es el vector de estimación inicial de lentitud/velocidad.

Parámetro p y su relación con la convergencia

Otra forma de expresar la ecuación (9) se indica a continuación:

El número máximo de valores singulares de M (p_mXx) será el producto de emisores por receptores: m.n; es decir, igual a su número de filas. Esta ecuación indica que valores singulares pequeños li en problemas mal condicionados (como lo son los de tomografía sísmica) magnificarán los errores en el modelo en M, y además errores de medición en t. Es posible disminuir los errores restringiendo el número de valores singulares descartando los muy pequeños y considerando solo aquellos suficientemente grandes. En este trabajo, luego de la inspección de la figura 7, se probó realizar la inversión de datos utilizando 3 posibles valores de p: 40, 46 y 49. Los resultados se muestran en las figuras 8 a 10.

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 necesita 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, lo que constituye una nueva herramienta computacional para resolver la inversión de datos de tomografíasísmica en columnas cilíndricas mediante métodos matriciales.

Además, el algoritmo posibilita 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 yla distribución de sensores en la periferia y de pixeles en el interior del dominio [9], lo cual abre nuevas posibilidades de investigación.

La validación demostró el adecuado funcionamiento del algoritmo, resultando la inversión en la localización precisa de las anomalías existentes en el dominio.

Por las conclusiones anteriores será posible diseñar experimentos reales con óptimo número de sensores a fin de detectar a muy bajo costo inclusiones en medios aproximadamente homogéneos por ejemplo, oquedades en columnas cilíndricas de hormigón, o bien el volumen de agua transportada por una cañería, etc.

Reconocimientos

Este trabajo de investigación forma parte de los proyectos "Microtomografía sísmica y determinación de módulos dinámicos con ondas P y S en medios elásticos con arreglos genéricos" y "Desarrollos de métodos numéricos para detección de fallas en medios homogéneos y heterogéneos con técnicas tomográficas", subvencionados por la Secretaría de Ciencia y Técnica (CICITCA) de la Universidad Nacional de San Juan, República Argentina.

REFERENCIAS

[1] J.C. Santamarina, K. A. Klein and M. A. Fam, Soils and Waves. New York: Wiley & Sons, 2002.         [ Links ]

[2] C. A. Calvo and A. L. Imhof, "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. vol. 35, n°. 1, 78-87, 2014.         [ Links ]

[3] A. L. Imhof, 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, 2008.         [ Links ]

[4] W. Menke, Geophysical Data Analysis: Discrete Inverse Theory. New York: Academic Press, 1989.         [ Links ]

[5] J. C. Santamarina and D. Fratta, Discrete Signals and Inverse Problems: An Introduction for Engineers and Scientists. New York: Wiley, 2005.         [ Links ]

[6] A. Tarantola, Inverse Problem Theory. Amsterdam: Elsevier, 1987, p. 613.         [ Links ]

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

[8] R. A. Penrose, "A Generalized Inverse for Matrices", Proceedings of Cambridge Phil. Society, n° 51, pp. 406-413, 1955.         [ Links ]

[9] A. L. Imhof, C.A. Calvo, and J. C. Santamarina,"Seismic Data Inversion by Cross Hole Tomography Using Geometrically Uniform Spatial Coverage", Brazilian Journal of Geophysics, vol. 28, n° 1, pp. 79-88, 2010.         [ Links ]