ANAL/SIS DE SENSIBIL/DAD EN LA ESTIMACION DE LA FUNCION CARACTERISTICA DE PUNTa DE HAMILTON DE LA SISMICA GEOMETRICA LUIS ALFREDO MONTES VIDES Profesor Asociado, Departamento de Geociencias - Universidad Nacional de Colombia Montes, L (2000): Analisis de sensibilidad en la estirnacion de la funcicin caracteristica de Punto de Hamilton de la sismica geometrica. Geofis. Colomb. 4:31-36. ISSN 0121-2974 Santa Fe de Bogota, D.C. - Colombia RESUMEN Usando la aproxirnacion de segundo orden de los tiempos de Transite en modelos sismicos constituidos por capas hornoqeneas e isotr6picas separadas por interfaces suaves, los tiernpos de los eventos de reflexi6n registrados en vecindad del rayo central pueden estimarse a partir de una funci6n hiperb61ica que corresponde a la funci6n caracteristica de punto de Hamilton, la cual no es conocida nurnericarnente en la optica geometrica. En la sismica, esta se puede determinar a traves de una regresi6n lineal entre las coordenadas fuente-detector y los tiempos de Transite de por 10 menos nueve rayos apropiadamente seleccionados y registrados en configuraciones de cubrirniento multiple. Un anal isis de sensibilidad permite conocer la estabilidad y robustez del rnetodo de regresi6n que estima la funci6n caracteristica de los tiempos de Transite. ABSTRACT Using the second-order approximations of the traveltimes in a seismic model which consists of a pile of homogeneous layers separated by smooth interfaces, the traveltimes of reflection events registered in the vicinity of the central ray can be estimated from an hyperbolic function that corresponds to Hamilton's point characteristic in optics, where it is not numerically known. In seismic, it can be determined through a linear regression between the source-detector coordinates and the traveltimes of at least nine rays properly selected and registered in multiple fold configurations. A sensibility analysis allows to know the stability and robustness of the linear regression method to estimate the characteristic traveltimes function. 1. INTRODUCCION Las leyes generales para representar propiedades globales de los rayos transmitidos y reflejados en un sistema sismico fueron establecidas usando las ideas de la optica qeometrica (Bortfeld, 1989). En un modelo de la tierra constituido por capas hornoqeneas e isotropicas separadas por interfaces suaves con velocidades y densidades arbitrarias, el sistema queda limitado por la superficie anterior considerada plana (de adquisici6n), y la posterior arbitrariamente curvada (de reflexi6n). EI tiempo de Transite de cualquier rayo cercano a un rave de referencia denominado rayo central se expresa por una sencilla funci6n rnatricial, que es equivalente a la funci6n caracteristica de punta de Hamilton de la optica. Maternaticamente, la teoria de aproximacion paraxial es la contra parte sismica de la 6ptica qausstana. Problemas, tales como fenornenos de enfocamiento (Bortfeld & Kemper, 1991) Y miqracion de amplitudes (Bortfeld & Kiehn, 1992) han sido tratados exitosamente sin conocer la estructura detallada del sistema sismico. La teoria provee formulas mas simples y compactas, deseables en problemas de inversion. EI proposito de este articulo es evaluar el proceso que estima la funci6n que suministra el tiempo de un evento reflexion asociado a una superficie del sistema, y analizar como los errores cornetidos en la identificacion de los primeros arribos puede influir en el estimative de dicha funcion: error que se transfiere a los resultados de cualquier proceso de inversion que haga uso de esta teoria. 2. RAVOS V TIEMPOS DE TRANS ITO El sistema se considera compuesto de capas homoqeneas e isotropicas separadas por interfaces suaves, esta limitado por la superficie sobre la cual se posicionan las fuentes y los geofonos lIamada superficie anterior, y por la superficie donde tiene lugar la reflexion, denominada superficie posterior. En la aproxirnacion geometrica 0 paraxial, un rayo transmitido que incida normalmente sobre el reflector se usa como rave central, a este se referencian los rayos vecinos 0 paraxiales. EI sistema de coordenadas se fija al rayo central, haciendo coincidir el origen del sistema cartesiano con el punto donde el rayo intersecta la interface, asi el plano xyes tangente a la superficie anterior en el punta de partida del rayo central. El plano xy' es tangente al punto de interseccion del rayo central con la superficie posterior. Con la ayuda de estos sistemas cualquier rayo paraxial puede representarse por los vectores posicion y lentitud (slowness) proyectados sobre los pianos tangentes (xy ,xy') en las superficies anterior (X, P) y posterior (x' , p'), ver Fig.1. Para diferenciar los vectores asociados al rayo central de aquellos asociados a los rayos paraxiales, la pareja de vectores asociados al rayo central se escribe con el suscrito "0", asi en la superficie anterior la pareja sera (x 0 ' Po)' En la aproximacion de segundo orden de los tiempos de Transite las parejas de vectores en la superficie de reflexion se relacionan con los vectores de la superficie anterior a traves de la ecuacion: [X'] [Ao Bo][ X J p' = Co Do P-Po (1) con Ao,Bo,Cr,,12o matrices [acobianas 2x2. Un rayo trazado desde X en la superficies anterior a x' en la posterior y de all[ hasta x" sobre la anterior representa una onda que atravesando el sistema es reflejada con un tiempo de Transite: I(X,Y) =~ +Po(X+Y)+~(X+Y)D~ICo~(X+Y) +~(Y -x)s;;'Ao~(Y -x) (2) donde 7;, es el tiempo de Transite del rayo central de ida y vueIta (two way time). EI tiempo de Transito zero-offset se obtiene cuando x = x' . Superficie posterior FIGURA 1. SISTEMA SISMICO MOSTRANDO EL RAVO CENTRAL CON LOS SISTEMA CARTESIANOS DEFINIDOS EN LAS SUPERFICIES ANTERIOR V POSTERIOR. EL RAVO PARAXIAL ATRAVIESA EL SISTEMA HAST A EL REFLECTOR. La ecuacion (2) corresponde a la caracteristica de punta de Hamilton en optica, lIamada asi por dar el tiempo de propaqacion entre dos puntos. Aun usada en el disefio de instrumentos opticos la funcion caracteristica de punto de Hamilton no se conoce explicita ni numencamente. Sin embargo en sismologia el sistema es la Tierra y fa forma de la funcion se expresa en (2), pudiendo estirnarse nurnericarnente. Las matrices B~JAo,/2~Jf:o son sirnetricas, en consecuencia la funcion (2) posee nueve incognitas: seis en B~I Ao,D~ICo, una en 7;, y dos en Po' Utilizando al menos nueve tiempos de transite apropiadamente registrados las incognitas de fa expresion (2) pueden resolverse y por tanto estimar numericamente la funcion caracteristica de los tiempos de transite. 3. CONFIGURACION V TIEMPO DE ARRIBO Para estimar la funcion caracteristica se requiere de por 10 menos nueve mediciones (picking) de tiempos de arribo de las ondas reflejadas, ya que la ecuacion (2) posee nueve incognitas. Para simplificar las ecuaciones denotaremos los parametros as;: Una configuraci6n con disparo en el origen del sistema cartesiano da lugar a una familia de trazos fuente comun (common-shot). AI situar los geofonos en las coordenadas (x gi , Y gi) Yal usar la ecuaci6n (2) se tiene un sistema de ecuaciones que matricialmente tiene la forma: Ix y 2 X Y y2 2 X Y 2g9 g9 XK9 1'9 g9 1'9 xg9 g9 g9 Yg9 EI sistema presenta tres pares de columnas iguales 0 linea/mente dependientes, por tanto no admite soluci6n. Este resultado indica que es imposible conocer la funci6n caracteristica usando una familia de trazos con fuente comun. Un analisis similar en trazos obtenidos en receptor comun CRP /leva a un sistema con seis grados de libertad, imposibilitando conocer los parametres a partir de una familia de trazos CRP. AI considerar trazos registrados en configuraci6n CMP se obtiene el sistema de ecuaciones: = Por el hecho de contar con seis grados de libertad se concluye que en un unico arreglo CMP no es posible conocer la funci6n caracteristica. -p -q u~ Por ultimo en una familia de trazos common-offset nuevamente se obtiene un sistema con seis grados de libertad encontrando /a misma imposibi/idad hallada en los tres casos anteriores. Lo anterior permite afirmar que para determinar /a funci6n caracteristica se requiere utilizar tiempos de arribo registrados en distintas configuraciones. 4. APROXIMACION HIPERBOLICA En la exploraci6n sismica se sabe hace tiempo que en medios con capas planas y horizontales, los tiempos de transite de las reflexiones aproximadamente verticales se ajustan mejor por funciones hiperb61icas que por funciones de tipo parab61icas (Ursin, 1982). AI considerar el cuadrado de la funci6n de tiempos de transite dada en (2) y despreciando los terrninos mayores a la segunda potencia se obtiene: t(x,t"l =[J: -2])0 -~(x+xjf +2J:G(x+xj- D~'Co~(X+~)+~(~ -x)-B,;l&i(~ -x)] (3: La ecuaci6n (3) representa los tiempos de transite en la aproximaci6n hiperb6lica. EI conjunto de N tiempos registrados en distintas configuraciones da lugar a un sistema matricial de N ecuaciones. 5. ANALISIS DE SENSIBILIDAD La estabilidad se define como la propiedad que posee un estimador central a ser insensible a la presencia de pequefios errores aleatorios en los datos, al suministrar informaci6n "a priori" sobre el modelo se puede garantizar la estabilidad en /a soluci6n. La robustez, indica la sensibilidad del pararnetro a un numero reducido de errores apreciables en los datos. Estimar la funci6n caracteristica usando tiempos de transito de eventos de reflexi6n es un problema de regresi6n lineal de un sistema de ecuaciones sobredeterminado que se expresa de la forma: Minimizar (4) La soluci6n para (1) seran los va/ores de tfl que generen el menor error If , que se obtiene a traves de la matriz inversa generalizada, ver Apendice A. Para estudiar la sensibilidad en la determinaci6n de la funci6n caracteristica, la ecuaci6n (3) se perturba introduciendo errores aleatorios en los tiempos de transite. Estas secuencias de nurneros aleatorios se generan con una distribuci6n normal con media fl = 0 y desviaci6n estandar 0" * 0 . Se lIevaron a cabo nueve ensayos, con gaussianas de 0", ={2,4, ...,18ms}, el valor mayor representa el tiempo usado por un frente de onda con frecuencia de 100 Hz para viajar una longitud de onda A a una velocidad de 1000 m/s. Con este intervalo se indica que el error se relaciona con la incertidumbre para identificar en una longitud de onda el tiempo de arribo (picking). Con fuentes y deteclores en configuraci6n CMP la funci6n caracteristica viene dada por la expresi6n: t(O, X) = To + xUX, (5) donde U = B~I A 0 y To se pueden estimar a traves del metoda de la inversa generalizada. Por otro lade cuando fuente y detector estan en configuraci6n zero-offset la funci6n caracleristica de tiempos de transite toma la forma: lex) = To + J\x + xVx (6) donde E =D~JCo Y Pose obtienen usando la matriz de Moore-Penrose. ver Apendice A. El tiempo para una configuraci6n zero-offset se puede aproximar al tiempo de arribo al detector mas proximo a la fuente. Ligeras perturbaciones se introdujeron en los tiempos de transite del sistema de ecuaciones (3) y se estudi6 como los errores en esos tiempos de arribo influyen en los parametres, para ello se cre6 un programa en Matlab 4.0 que calcula la funci6n caracteristica exacta a partir de tiempos de transite obtenidos de un modelo de capas plano horizontales, utilizando el programa de trazamiento de rayos SEIS88 desarrollado por Cerveny et al. (1980) en la Universidad de Charles en Praga. Posteriormente se adicionaron a los tiempos de arribo cadenas de errores estocasticos generados por una funci6n de distribuci6n normal, con media cera y 0", * 0; los resultados arrojados en el proceso se muestran en la Tabla 1. TABLA 1 a(s) ()(s) p, r, ~I 1~2 Vn VI1 V12 Vn 0 0.5 0 0 1 0 1 0 0 0 0.002 0.4997 -0.0006 -0.0009 0.9937 0.0109 0.9922 0.0011 0.001 0.0021 0.004 0.4993 -0,0011 -0.0017 0.9874 0.0219 0,9844 0.0026 0.0011 0.0041 0.006 0.499 -0.0015 -0.0026 0,981 0.033 0.9766 0.0045 0.0002 0.0062 0.008 0.4987 -0.002 -0.0034 0.9746 0.0443 0.9687 0.0068 -0.0016 0.0084 0.01 0.4984 -0.0024 -0.0042 0,9681 0,0556 0.9608 0.0095 -0.0043 0.0106 0.012 0.4981 0.0027 -0.005 0.9615 0.067 0,9528 0.0126 -0.008 0.0128 0.14 0.4978 -0.003 -0.0057 0.9548 0,0785 0.9448 0,0161 -0.0126 0.015 0.16 0.4975' -0.0033 -0,0065 0,9481 0.0901 0.9368 0.0199 -0,0181 0.0173 0.18 0.4973 -0.0035 -0,0072 0.9413 0,1019 0.9287 0.0242 -0.0246 0.0196 Los resultados muestran que en el caso extrema de maximo nivel de ruido en los tiempos de arriba (0"=18ms), los parametros J3;,IAo (U) y To se obtuvieron con errores bajos (menores a 3%), indicando la buena estabilidad de los parametres. AI perturbar a nivel maximo los tiempos en zero-offset se hall6 que el Po estimado no se desvi6 significativamente de su valor verdadero, mientras que D;,I£ () se calcul6 con un error de 10%. Para disminuir ese nivel de incertidumbre los tiempos zero-offset deben ser calculados a partir del ajuste de datos en configuraci6n CMP. En general los resultados indican que a mayor error en estimar los tiempos (picking) mayor error en la determinaci6n de los parametres. CONCLUSIONES En el presente trabajo se muestra el procedimiento para estimar la funci6n caracteristica de los tiempos de transite en la aproximaci6n hiperb61icausando para ello un conjunto de tiempos de reflexi6n registrados en cobertura multiple. Un anallsis de sensibilidad sobre el procedimiento de estimaci6n de la funci6n caracteristica indica que se pueden conocer los parametres de esta funci6n con un pequeno margen de error, el cual esta ligado a la calidad de los datos colectados, La funci6n estimada representa los tiernpos de eventos de reflexion en la vecindad del rayo central en el rango de offsets donde es valida la correcci6n NMO. AGRADECIMIENTOS Este trabajo es un resultado parcial obtenido durante mis estudios conducentes a Ph.D. en la Universidad Federal de Para en Brasil, bajo la direcci6n de los Doctores Walter Sollner y Peter Hubral a quienes manifiesto un gran agradecimiento por sus valiosas ideas y apoyo durante el desarrollo del traba]o. as! como en la sustentaci6n del mismo. REFERENCIAS BIBLIOGRAFICAS Arthanari, T. S. & Y. Dodge. (1981): Mathematical Programming in Statistics. J, Wiley and Sons. Bortfeld, R. (1989): Geometrical Ray Theory: Rays and Traveltimes in Seismic Systems. Geophysics 54:342-349. Bortfeld, R., & R.M. Kemper. (1991): Geometrical Ray Theory: Line Foci and Point Foci in the Anterior Surface of Seismic Systems: Geophysics 56:806-811. Bortfeld, R.. & M. Kiehn. (1992): Reflection Amplitudes and Migration Amplitudes, Geophys. Prospect. 40:873- 884. Cerveny, V., L. Klimes & V. Psencik. (1988): Complete Seismic-Ray Tracing in 3-D Structures. D. Doornbos Ed. Seismological Algorithms. Academic Press. Dantzig, G. B. (1963): Linear Programming and Extensions. Princeton University Press. Menke, W. (1973): Geophysical Data Analysis: Discrete Inverse Theory. Academic Press. Ursin, B, (1982): Quadratic Wavefront and Traveltime Approximations in Homogeneous Layered Media with Curved Interfaces. Geophysics, 47:1012-1021. APENDICEA SOLUCION DEL PROBLEMA INVERSO Sea el sistema algebraico lineal de n-inc6gnitas (mt, m2 , ... , mn) relacionadas a traves del sistema acoplado de r-ecuaciones dado por: gttmt + gt2m2 + + gtnmn = d, g21ml + g22m2 + + g2nmn = d2 (A-1) Donde gIl Y b, tienen valores conocidos, Cuando el numero de ecuaciones es mayor que el nurnero de inc6gnitas el problema se dice que es sobredeterminado. Matricialmente el sistema de ecuaci6nes A1 se representa asi: (A-2) La soluci6n consiste en estimar los parametres del modelos dado par m'" I de modo que los datos previsibles a pre se aproximen 10 mejor posible a los datos observados a ohs . Sus diferencias se expresan por el vector de error If = a obs - J pre, siendo la mejor solucion aquella que genere el menor error medido a traves de la norma del error IIL( If )l!. Las normas mas usadas son la de error absoluto minimo y la de error cuadratico medio minimo definidas respectivamente como: (If)11= 2::)&, y (----""-~~ (If)II=/''I& 2 "JL...A ' Para cada norma se supone que los datos se comportan sequn un tipo de estadistica, as! para la norma 2 los datos obedecen a una distribuci6n de Gauss y en este caso la soluci6n al sistema A-2 viene dado par el metoda de la matriz inversa generaJizada de Moore-Penrose (Menke, 1987): (A-3) el superindice t indica traspuesta. Suponiendo que existe la matriz inversa (0' arl se obtiene: (A-4) En el caso de la norma 1, para el sistema A-2 la soluci6n se establece a partir de su representaci6n como un problema de programaci6n lineal (Arthanari & Dodge, 1981). La representaci6n sequn la teoria de investigaci6n de operaciones viene dado par: Minimizar 2:: Sujeto a: Om - If = a , (A-5) El cual fue resuelto par el denominado rnetodo Simplex 0 Simplex revisado de la programaci6n lineal (Dantzig, 1963). Generalmente los datos contienen ruido que causa error al estimar los parametros del modelo, si los datos tiene asociada una funci6n de distribuci6n caracterizada par una matriz de covarianza dada por [cov(ffi)]=MaM , donde _M es la matriz inversa generalizada de Moore- Penrose. Si los datos no estan correlacionados y tienen una varianza a~ entonces la solucion que genera el menor error cuadratico tendra una matriz de covarianza dada por: que resulta tinalmente en: (A-6) La expresion A-6 indica que fa matriz de covarianza de los parametres esta controlada por la rnatriz de covarianza de los datos y por la matriz inversa generalizada. La introduccion de informacion "a priori" garantiza la obtencion de la solucion en el problema de minimos cuadrados, esa informacion es una restriccion sobre los valores que tomaran los parametres, Esto se presenta como el caso donde un parametro debe tomar un valor constante conocido 0 multiplicado de una matriz de pesos ponderados, dando un mayor valor a aquellos que presentan mayor estabilidad. En el caso particular de ajuste a una parabola la ecuacion A-2 se expresa en forma matricial asl: (A-4) EI calculo de ff1 por minimos cuadrados de A-4 segun la inversa generalizada A-4, sera: