Solución Analítica-Numérica para la Trayectoria de Rayos Sísmicos en Medios con Heterogeneidad Vertical

Armando Luis Imhof,

Carlos Adolfo Calvo,

Manuel Sánchez

Resumen

La propagación de ondas sísmicas generadas en emisores y detectadas en receptores en medios heterogéneos y anisotrópicos sigue una trayectoria curva determinada por el principio de tiempo mínimo de Fermat. La aplicación del cálculo de variaciones conduce a una ecuación diferencial ordinaria. Debido a la compactación del terreno normalmente la velocidad aumenta con la profundidad, estando esta variación sujeta a leyes experimentales para cada suelo, lo que lleva a resolución numérica en los casos generales. En este trabajo se estudia el ajuste de datos experimentales de velocidad por una función exponencial; la integración analítica de la ecuación diferencial y la determinación numérica de las constantes de integración. Los datos experimentales suelen determinarse por algún método geofísico tal como up-hole o down-hole. Su aplicación principal se centra en el diseño de modelos numéricos de trayectorias curvas y por ende luego tiempos de primeros arribos a través de algoritmos de inversión tomográficos para detección y modelado de anomalías en los primeros 12m. de profundidad.

Palabras Clave: Anisotropía, Rayos curvos, principio de fermat, teoría del rayo, heterogeneidad vertical

Abstract

The propagation of seismic waves generated by transmitters and receivers acquired in heterogeneous and anisotropic media follows a curved path determined by the principle of Fermat. The application of the calculus of variations leads to an ordinary differential equation. Because soil compaction wave velocity typically increases with depth, being this variation subject to experimental laws for each subsoil. That leads to numerical resolution in the general case. In this work, the adjustment, by an exponential function, of experimental velocity data set is studied; besides analytical integration of the differential equation and the numerical determination of the constants of integration. Experimental data are usually determined by a geophysical method as up-hole and down-hole. Main application is focused on the design of numerical models of curved rays and therefore first picks for wave arrivals for the construction of algorithms for inversion of seismic data for detection and modeling of velocity anomalies up to the first 12m. depth.

Keywords: Anisotropy, Curved Rays, fermat principle, ray theory, vertical heterogeneity

Recibido: 19 de noviembre 2014 Aprobado: 20 de febrero 2015

1. Introducción

La propagación de ondas sísmicas generadas por medios tales como golpes de martillo, explosiones o cristales piezoeléctricos (emisores); detectadas por geófonos (receptores) y registradas en sismógrafos, está regida por la ecuación de onda, cuya solución completa implica resolver la distribución de frentes de ondas que se propagan en el medio considerado, que puede ser homogéneo, heterogéneo, isótropo, anisótropo o cualquier combinación entre heterogeneidad y anisotropía.

Sin embargo, muchos casos suelen utilizar una simplificación considerando a los frentes como rayos sísmicos (Cerveny, 2001) que viajan en medios tanto homogéneos como heterogéneos y/o anisotrópicos. En medios homogéneos los frentes de ondas serán esféricos y los rayos normales a los mismos y rectos. Esto no sucederá en los medios heterogéneos donde los rayos serán curvos, a causa de la refracción producida por el cambio de velocidades sísmicas de acuerdo a la ley de Snell. Además, en los medios anisótropos los rayos no serán normales a los frentes de ondas. Por otra parte, según se demostró (Santamarina & Cascante, 1996) la anisotropía considerada únicamente, no produce curvatura del rayo, sino que la acentúa. Es decir, en síntesis, los rayos que atraviesen un medio homogéneo y anisótropo serán rectos, pero no normales a los frentes de ondas.

Para resolver las trayectorias de los rayos se han desarrollado métodos y algoritmos numéricos denominados trazadores de rayos (eng. ray tracing) basados en la teoría del rayo sísmico (eng. ray theory). Debido a la dificultad que los planteos analíticos proponen, hasta ahora siempre se han llevado a cabo aproximaciones numéricas (Menke, 1989).

En este trabajo se desarrolla una metodología analítica-numérica con el propósito de resolver la trayectoria de rayos sísmicos que se transmiten en medios de diferentes heterogeneidades y anisotropías verticales, no aplicándose para variaciones laterales de velocidad, lo que llevaría otro tipo de desarrollo.

La motivación principal ha sido la de contar con una herramienta para poder diseñar algoritmos computacionales para resolver problemas de tomografía sísmica en suelos muy superficiales (0-12m) en los cuales se demostró que se puede aproximar la variación de velocidad incremental con una función exponencial. Una vez conocida la variación de velocidades en profundidad en algún punto del sector de trabajo (por medio de up-hole, o down hole superficial). Esto se concluyó por la inspección visual de los resultados experimentales obtenidos (Imhof, 2007). Por otro lado resulta también posible estimar la forma de incremento de la velocidad con la profundidad, a partir de aproximaciones sucesivas para resolver, por ejemplo una tomografía sísmica, en caso de utilizar esta metodología para tal fin.

En su propagación por el subsuelo, los rayos siguen una trayectoria curva determinada por el Principio de Fermat, que establece que los rayos seguirán trayectorias de tiempo mínimo y no de distancia mínima. El desarrollo analítico de la trayectoria para un caso de heterogeneidad vertical fue llevado a cabo por Santamarina, Klein & Fam (2002). Esto produce que en caso de medios heterogéneos los rayos se curven.

2. Formulación del problema y modelo matemático

En un medio plano heterogéneo y anisótropo verticalmente, se ubica el emisor de onda E y el receptor R.

La distribución de velocidades V(θ ) según la inclinación θ, sigue una ley elíptica (Calvo, Imhof & Jofre, 2008). El coeficiente de anisotropía d=Vv/Vh se considera como constante, lo que impli ca considerar que las variaciones verticales de velocidad van acompañadas del mismo orden en las horizontales, siendo Vv , Vh las velocidades correspondientes a las direcciones principales (vertical y horizontal), ver Figura 1.

Figura 1. Distribución elíptica de velocidades

El vector es tangente a la curva que representa la trayectoria y=y(x) siendo

dy /dx=y'=tag(θ)=Vy /Vx

V(θ) = √Vx2+Vy2 representa el modulo de

La ecuación de la elipse resulta:

V(θ ) = Vv√1+y'2 /d2+y'2 (1)

Para Vv = cte el medio es homogéneo. Para d=1 el medio es isotrópico. Debido a la compactación diferencial del terreno se produce una heterogeneidad donde Vv=Vv(y), ver Figura 2.

Figura 2. Variación de velocidad con la profundidad

Siguiendo el principio de Fermat de tiempo mínimo, debe ser minimizado el funcional tiempo.

Para ello se utiliza la fórmula de Euler del cálculo de variaciones (Elsgoltz, 1969), que para este caso es F-y'Fy' = C(cte) , lo que conduce a una ecuación diferencial ordinaria de primer orden:

Vv (y) √d2+y'2 = C (3)

Las coordenadas del emisor y receptor constituyen las condiciones de contorno:

y(xe) = ye y(xr)=yr (4)

La constante C es característica de cada rayo. La integración de la Ecuación (3), depende de la expresión Vv (y) la que se determina a partir de datos experimentales, en general se resuelve numéricamente (Santamarina & Fratta, 2005). Para algunas expresiones Vv (y) tiene solución analítica. Las condiciones de contorno permiten encontrar y (x) y a la constante C.

3. Ajuste de la velocidad con una función exponencial

Imhof (2007) demostró que la variación de velocidad de onda compresional en los primeros metros bajo la superficie sigue la forma presentada en Tabla 1, para la mayoría de los suelos. Esta agrupación de datos se puede ajustar mediante una función exponencial del tipo:

Vv (y) = a-be-ay (5)

Para la construcción de la Figura 2 se dispone de una tabla (yi ,Vv( yi )) con i=1,2,3...n obtenida experimentalmente (Tabla 1). Se observa que la función Vv(y) crece a medida que aumenta la profundidad hasta un nivel en que se estabiliza. Este carácter asintótico es descrito adecuadamente por la función exponencial Vv = a-be-ay , donde:

a: Velocidad (máxima) asintótica (para y→∞) en el subsuelo; esto ocurre (aproximadamente) para profundidades iguales o mayores a 12m: yc ≅ 12m es la profundidad a partir de la cual se estabiliza Vv . Pudiéndose considerar como constante (en la escala de estudio). Si se considera que el término exponencial negativo es despreciable para argumentos iguales o superiores a 5 (e−5 ≅ 0), entonces αyc ≅ 5

b: es el cambio de velocidad debido a la heterogeneidad, resultando:

a-b: Velocidad (mínima) a nivel del piso

α: Coeficiente de atenuación de la velocidad según el medio (aproximadamente α = 5/yc)

Para el ajuste se determina el error local ei (en cada valor yi) y el error total de ajuste E:

ei =Vv (yi) - Vi = (a-be-αyi - Vi (6)

(7)

Donde (a, b, α) representan las variables de ajuste. Se buscan aquellas que hacen mínimo al error E, para su determinación se usa el algoritmo de Nelder - Mead (Shampine & Reichelt, 1997). Éste es un proceso iterativo que está implementado en MatLAB. Es necesario introducir la expresión de E y los valores iniciales (a(0), b(0), α(0))para que arranque el proceso iterativo. Los valores usados son 8:

a(0) = Vmax b(0) = Vmax - Vmin α(0) = 5 / yc (8)

Donde el supra índice señala la etapa iterativa (cero para la inicial). Luego de finalizar este proceso se obtiene el valor óptimo (a, b, α) con el que se construye Vv (y) = a-be-αy

4. Solución de la ecuación diferencial con la función exponencial Vv = a-be-αy

La Ecuación (1) a resolver es (a-be-αy) √d2+y'2 = C

Para las ecuaciones diferenciales de la forma h (y, y') = 0 en que falta la variable independiente x, existen métodos de resolución que consisten en la sustitución de y o y' por una expresión paramétrica (Elsgoltz, 1969).

Para este caso la sustitución propuesta es: (ver solución (9)

y' = d.tg(t) por lo que √d2+y'2 = d/cos(t)

4.1 Aplicación de las Condiciones de Contorno

El parámetro t ∈ al intervalo [t0, t1], correspondiendo t0 al emisor y t1 al receptor, reemplazando en Ecuación (9) e imponiendo las condiciones de contorno, se llega a un sistema de 4 ecuaciones no lineales con 4 incógnitas: (t0, t1, C, C1)

Para resolver el sistema (10), se puede aplicar:

El algoritmo iterativo de Newton:

Pasando a notación vectorial:

u = [u1, u2, u3, u4]T = [t0, t1, C, C1]T

f= [f 1, f 2, f 3, f 4]T J(u) = ∂fi / ∂uj

u(n+1) = u(n) + ∆u n= 0, 1, 2, 3,....

Donde ∆u surge de resolver el sistema lineal en cada etapa:

J(u(n)) • ∆u = − f(u (n))

La elección del valor inicial u(0) condiciona la convergencia del algoritmo que no siempre se logra, por ello se ha desarrollado un método específico para este caso, sin uso de derivada.

Método sin uso de derivada

Consiste en elegir C en el intervalo [Cmin, Cmax] que posibilite valores reales (x, y) siendo Cmin= d(a−be−αM). M es el máximo entre ye e yr. Ello surge de (ver ecuación 11):

Y de la expresión de f1 (o de f3), surge: C ≤ ad= Cmax

Luego se hace C = λCmin+(1−λ)Cmax eligiéndose 0≤λ≤1

Se calcula t0 y t1 de las expresiones f2 (t0,C) y f4(t1, C):

t0 = arcos (da−dbe−αye/C)

t1 = arcos (da−dbe−αyr/C)

Luego se elimina C1 entre f1 y f3 haciendo f1 (t0, C, C1) − f3 (t1, C, C1) = f (t0, t1, C) = 0 reemplazando en esta última ecuación las expresiones de t0, t1 y de C se llega a una expresión sólo función de

λ : f(λ) = 0

Para el cálculo de la raíz λ , se usa un método iterativo, estimando 2 valores iniciales:

λ1,λ2 e iterando según método de la secante:

Calculado λ, se determina C, t0, t1, C1

5. Aplicación a un caso determinado:

Datos: (xe, ye) = (0,1) ; (xr, yr) = (3,3); d = 1.3

Los valores de Vv según Tabla 1

Nota: Basado en Imhof, A. L. (2007). Caracterización de Arenas y Gravas con Ondas Elásticas por Tomografía Sísmica en Cross-Hole (Tesis doctoral en Ingeniería). Universidad Nacional de Cuyo. Mendoza. Argentina.

Ajuste de Vv=a−be−αy: Algoritmo de Nelder-Mead

Mínimo de

Los valores iniciales usados:

a(0)= 950 b(0)=660 α(0)=5/12

El ajuste arroja los valores

a=947.9536; b=654.9423; α=0.4428

Condiciones de contorno

M = yr=3 Cmin=1006.8 Cmax=1232.3 λ1= 0.8

Por iteración (método de la secante) resulta:

λ= 0.9066

Luego se calcula:

C= 1027.8630, t0=0.8407, t1= 0.2028, C1=4.7022

La trayectoria del rayo queda (ver 12):

El código MatLab propuesto para el ajuste de la función exponencial es el siguiente:

function E=faminimizar(c)

%E=faminimizar(c) es la función a minimizar

% c=[C1 C3 C3] =[a b alfa] es el vector que contiene los valores iniciales de los parámetros de ajuste

% DATOS

y=[0 2 4 6 8 10 12]; %profundidad

V=[290 690 840 900 930 940 950]; % velocidad

Va= c(1)-c(2)*exp(-c(3)*y);% valores ajustados

E=norm(V-Va,1);%devuelve el error de minimización

end

Para llamar la rutina se usa:

>>co=[950 660 5/12]; % valores iniciales Vmax, VVmin, alfa/yc

>>c=fminsearch('faminimizar',co); % minimiza a la función escrita en faminimizar

6. Extensión a otras expresiones Vv = f (y)

El método desarrollado es aplicable a numerosas expresiones Vv = f (y). Sólo requieren poder explicitar y de Vv (y) √d 2+y' 2 = C y luego poder integrar dx , además para adaptarse al ajuste deben tener y''<0, por ejemplo para Vv(y) = √a+by se obtiene como solución una cicloide. TambiénVv(y) = 4√a+by conduce a una solución adecuada, pero por su carácter asintótico la exponencial es la que mejor se adapta.

7. Conclusiones

Se ha desarrollado un método que permite encontrar expresiones analíticas para la trayectoria del rayo curvo en medios verticalmente heterogéneos y anisótropos. La imposición de las condiciones de contorno se logra con un método numérico desarrollado para este caso. La aplicación de este trabajo se centra en la validación de modelos numéricos.

Las conclusiones son importantes, ya que conociendo la ley de variación de velocidades verticales en los primeros 12 metros, se calculan las trayectorias y con las mismas los tiempos de arribo de los rayos. Con estos tiempos se pueden construir modelos de tiempos de primeros arribos teóricos de tomografía sísmica y compararlos con aquellos reales medidos y desarrollar por ende algoritmos numéricos de inversión por aproximaciones sucesivas, los que posibilitarán modelar anomalías irregulares de mayor o menor velocidad a la del medio que las contiene, insertas en medios heterogéneos y anisótropos, cuales son todos en los primeros metros del subsuelo. Este es el paso siguiente de la investigación.

Cabe destacar que la presencia de niveles acuíferos no impide la aplicación de este método, puesto que siempre de antemano hay que hacer una estimación previa de valores de velocidades según Tabla 1, que luego se ajustarán con el algoritmo numérico propuesto.

Referencias bibliográficas

Calvo, C.; Imhof, A. & Jofre, S. (2008). Determinación de la Matriz de Distancias en un Medio Heterogéneo y Anisótropo. En: IV Congreso Internacional de Matemática Aplicada. Universidad de Buenos Aires (UBA). Argentina

Cerveny, V. (2001) Seismic Ray Theory. New York: Cambridge University Press.

Elsgoltz, L. (1969) Ecuaciones diferenciales y cálculo variacional. Rusia: Ed. MIR.

Imhof, A. L. (2007). Caracterización de Arenas y Gravas con Ondas Elásticas por Tomografía Sísmica en Cross-Hole (Tesis doctoral en Ingeniería). Universidad Nacional de Cuyo. Mendoza. Argentina.

Menke, W. (1989). Geophysical Data Analysis: Discrete Inverse Theory. NY. USA: Academic Press, Inc.

Santamarina, J.C. & Cascante, G. (1996) Stress anisotropy and wave propagation – A micromechanical view. Canadian Geotechnical Journal, 33, 770-782.

Santamarina, J.C. & Fratta, D. (2005) Discrete Signals and Inverse Problems. An Introduction for Engineers and Scientists. UK: Wiley & Sons.

Santamarina, J.C.; Klein, K.A.. & Fam, M.A. (2002). Soils and Waves. Ed. Wiley & Sons

Shampine, L. F & Reichelt, M. W. (1997) The MATLAB ODE Suite, SIAM Journal on Scientific Computing, 18-1

Sobre los autores

Carlos Adolfo Calvo. Departamento de Matemática, Facultad de Ingeniería. UNSJ , San Juan, Argentina.

Correo electrónico: ccalvo@unsj.edu.ar.

Armando Luis Imhof. Departamento de Geofísica y Astronomía – Inst. Geofísico Volponi, Fac. de Cs Exactas, Fcas y Nat., UNSJ – Argentina.

Correo electrónico: aimhof@unsj.edu.ar.

Manuel Sánchez. Instituto de Mecánica Aplicada, Facultad de Ingeniería. UNSJ , San Juan, Argentina.

Correo electrónico: msanchez@unsj.edu.ar

Analia Moyano. Departamento de Matemática, Facultad de Ingeniería. UNSJ , San Juan, Argentina.

Correo electrónico: anamyo75@yahoo.com.ar