Geofisica Colombiana N° 5 p P 1 2 - 1 7 diciembre de 2001 Bogota D.C. ISSN - 0121 - 2974 Modelacion del problema inverso en geoeh~ctrica 20 mediante elementos finitos LUIS ALBERTO BRICENO Profesor a s o c ia d o Departamento de Geociencias - Facultad de Ciencias - Universidad Na c io n a l de Colombia PEDRO MAURICIO AVELLANEDA L6PEZ Ingeniero civil - Estudiante de Maestria en RecursDs hidr aulic o s Departamento de Ingenieria Civil - Facultad de Ingenieria - Universidad NaciDnal de Colombia RESUMEN Se plantea un modelo nurnerico para la solucion del problema inverso en geoelectrica 2D utilizando el metodo de los elementos finitos. La rnetodologia de solucion se fundamenta principal mente en dos procedimientos: a) solucionar nurnericamente la ecuacion fundamental de flujo de corriente para unos valores de resistividad establecidos, para luego calcular su respuesta electrica (resistividades aparentes calculadas), y b) minimizar la diferencia entre las resistividades aparentes medidas y las calculadas utilizando un model a de inversion. El metodo es iterativo. Para el primer ciclo se establece un modelo de bloques 0 parametres hornogeneo, que se modificara sucesivamente por el proceso de inversion, procurando que el error minimo cuadratico sea menor que 10%. Luego se elabora un programa de computador que sirva como herramienta de calculo y que represente graficarnente los resultados obtenidos para facilitar el proceso interpretativo. Finalmente se present an dos casos de secciones tomograficas elaboradas en el Parque Nacional de la ciudad de Bogota D.C. PALABRAS CLAVE' RESISTIVIDAD, TOMOGRAFIA, ELEMENTOS FINITOS, MODELACI6N NUMEAICA 20. ABSTRACT A numerical model for the solution of the inverse problem in geoelectric 2D is proposed using the method of the finite elements. The solution methodology is based mainly on two procedures: a) to solve numerically the fundamental equation of current flow for established values of resistivity to calculate its electrical answer (apparent resistivities calculated), b) to diminish the difference between apparent resistivities measures and the calculated ones using an investment model. The method is iterative. For the first cycle, a model of blocks or parameters settles down homogenous, that successively will be modified by the investment process, trying that the quadratic minimum error is smaller to 10%. A computer program is elaborated that serves like calculation tool and graphically represents the obtained results to facilitate the interpretative process. Finally two cases of tomografic sections in the National Park in Bogota D.C. are presented. KEY WORDS: RESISTIVITY, TOMOGRAPHY, FINITE ELEMENT, 2D NUMERICAL MODELLING. INTRODUCCION Los model os numericos como el Metodo de los elementos tinitos (MEF) se utilizan como herramienta de soluci6n cuando se tiene una ecuaci6n de campo que gobiema el flujo, y ademas se conocen unas condiciones de frontera. Para el caso particular, en una tomograffa de resistividades electricas (TRE) tinalmente se llega a un mapa de resistividades que resulta de modelar el flujo de una corriente a traves del subsuelo. En el trabajo que aquf se presenta se analiza el modelo bidimensional, que asume que el rumbo de las estructuras geol6gicas es normal al arreglo de electrodos empleado en superficie, es decir, que la geologia perrna- nece con stante. EI procedimiento que se utiliza en este trabajo consiste en en- contrar la respuesta electrica de una primera hip6tesis que repre- sente las caracteristicas del subsuelo, y compararla con la obtenida en campo, para posteriormente efectuar una correcci6n a esta hip6- tesis inicial, y as! sucesivamente encontrar algun modelo que aproxi- me de mejor manera el comportamiento electrico registrado en cam- po. El MEF se utiliza entonces para encontrar la respuesta electrica del modelo hipotetico. 12 En general, el rnetodo consiste en discretizar el medio continuo mediante un conjunto de puntos que unidos formen una red de ele- mentos, e individualmente se comporten como un medio hornoge- neo, es decir, con resistividad constante. De esta manera se mode- Ian medios heterogeneos; basta s610 utilizar un valor de resistividad distinto para cada elemento. Las ventajas de este metodo radican en la posibilidad de simular medios heterogeneos con relativa co- modidad, ya que s610 se debe ajustar una malla apropiada de ele- mentos finitos; adernas, permite simular fronteras irregulares, tales como la topograffa de la superficie. La principal ventaja radica en que es un metoda que se puede sistematizar, y en el cual el grado de detalle depende de la densidad de elementos de la red. EI trabajo que aqui se presenta consiste en la modelaci6n del subsuelo mediante elementos tinitos utilizando un procedimiento de inversi6n de datos geoelectricos por minimos cuadrados y el rnetodo de quasi Newton (Loke y Barker, 1996). GEOFislCA COLOMBIANA, 5, DICIEMBRE DE 2001 MODELACION DEL PROBLEMA INVERSO EN GEOELECTRICA 2D MEDIANTE ELEMENTOS FINITOS MARCO TEORICO ECUACIONES CONSTITUTIVAS DEL FLUJO DE CORRIENTE ELI~CTRICA EN EL SUBSUELO Partiendo de la Ley de Ohm (Dey y Morrinson, 1979) se llega a la ecuaci6n 1, que describe matematicamente el flujo de corriente en un espacio heterogeneo tridimensional, en el cuallas fuentes de corriente son puntuales. -V.[ 1 Vl/>(x,y,z)l= Io(x )O(y)O(z )(1)p(x,y,z) 'J s s s La fuente de corriente se representa maternaticamente por la fun- ci6n Delta de Dirac 8. EI medio 10representa la resistividad p(x,y,z) en cualquier punto, y el potencial electrico la variable ¢(x,y,z). Debido a que el problema se analiza en el plano (x,z) se debe apli- car la transformada de Fourier, con el objetivo de representar todo el problema en el plano de analisis, De esta forma, la ecuaci6n (l) se transforma de la siguiente manera: En la ecuaci6n (2) la variable inc6gnita es el potencial transforma- do ¢t, y ademas se utiliza la variable de transformaci6n kyoLa selecci6n de la variable de transformaci6n se hace siguiendo la metodologfa pro- puesta por Weller y Kample (1996), la cual propone que la variable de transformaci6n 0 mimero de onda es funci6n de la abertura entre elec- trodos de corriente. k = 0.06 k YMIN a, ~ YMAX max 6 (3)<; amax=Maxima distancia entre la posici6n de un electrodo de corriente y un electrodo de potencial en superficie, para el problema en cuesti6n. amin=Mfnima distancia entre la posici6n de un electrodo de corriente y un electrodo de potencial en superficie, para el problema en cuesti6n. kymax= Maximo valor de la variable de transformaci6n. kymin= Mfnimo valor de la variable de transformaci6n. Es necesario recordar que para encontrar los potenciales reales se debe aplicar la transformada inversa de Fourier, de acuerdo con la ecua- ci6n (4): DISCRETIZACION MEDIANTE ELEMENTOS FINITOS UTILIZANDO EL METODO DE GALERKIN En general, el Metodo de los elementos finitos aplicado a un elemento de la red consiste en encontrar una soluci6n aproximada de la variable inc6gnita, que al ser remplazada en la ecuaci6n de campo (2) produce un error, ya que no es la soluci6n exacta. El Metodo de Galerkin con- siste en minimizar este error, ponderandolo por alguna funci6n, que para el caso particular seran las funciones de forma 0 de interpolaci6n GEOFlslCA COLOMBIANA, 5, OICIEMBRE DE 2001 que se utilizan para encontrar el valor del potencial dentro del elemento. Finalmente, para cada elemento se puede obtener una ecuaci6n matricial de laforma: Donde [C] representa la matriz de conductancia del elemento, (i el potencial en los nodos de la red, y S las fuentes de corriente. Para el presente caso se utilizan elementos rectangulas; por ella el subfndice cuatro en la ecuaci6n (5), pues representa el numero de nodos de elemento. Al sumar todas las contribuciones de los elementos se obtiene el sistema matricial global, que se representa de manera similar que la ecuaci6n (5), s610 que ahora involucrando la totalidad de los nodos de la malla. Solucionando el sistema de ecuaciones se encuentran los valores del potencial transformado en los nodos de la red. Posterior- mente se emplea la ecuaci6n (4) para encontrar el potencial real en los electrodos de corriente. PROCEDIMIENTO DE INVERSION Una vez calculada la respuesta electrica del modelo supuesto, se debe realizar un proceso de inversi6n de datos, en el cualla respuesta elec- trica calculada se relacione con la respuesta electrica obtenida en cam- po, para de esta forma obtener una matriz correctora del modelo inicial (parametres iniciales). La siguiente expresi6n se utiliza como ecua- ci6n de recurrencia en el proceso de inversi6n (Sasaki, 1989): ~g :Diferencia entre los datos medidos en campo gm en la iteraci6n (i) y las resistividades aparentes calculadas a partir de la modelaci6n me- diante elementos finitos. J :Matrizjacobiana que contiene las derivadas parciales de la resistividad del bloque respecto a los dernas, ~P : Vector correcci6n de los pararnetros del modelo. C: Filtro A: Factor de amortiguamiento. Cada vez que sea utilizada la ecuaci6n (6), se obtiene una matriz columna ~P correctora de los parametros iniciales que servira para la obtenci6n de un nuevo modelo. Pi + I = Pi +!1p (7) Los valores de resistividad del nuevo modelo de bloques se encuen- tra utilizando la ecuaci6n (7). ALGORITMO DE TRABAJO Los pasos por seguir son los siguientes: • Iniciar can un modelo hornogeneo; usual mente se torna el promedio de los logaritmos de las resistividades medidas (Lake and Barker, 1995). Los calculos se efecnian utilizando notaci6n logarftmica, pues de est a forma se evitan resultados negativos, y ademas se est a acorde con resultados experimentales, que rnuestran c6mo usual mente los cambios de resistividad en el subsuelo son 13 BRICENO £1 AL. l\'lALLA DE ELEMENTOS FINITOS - rt. --1 .. Nodos del elemen/ Mediante MEF se encuentra: q>1'm J I t---I--- »<: k m <1>, dentro del elemento se encuentra interpolando $. $. $k $ I' J' 'm Figura J. Malia de elementos finites que une los puntas 0 nodos que discretizan el continuo logarftrnicos. Luego se calcula el vector discrepancia entre los datos tornados y 1a respuesta del modelo. Como el primer modelo de resistividades es hornogeneo, el computo de la respuesta electric a mediante elementos finitos se puede omitir, pues te6ricamente debe ser igual al modelo supuesto. Se calcula la matriz jacobiana, luego se construye la matriz filtro C, y se termina con la soluci6n de la ecuaci6n (6). El factor de amortiguamiento se selecciona segun el nivel de perturbaci6n de los datos tornados. Generalmente se toma un maximo valor para la primera iteraci6n (A. = 0.2) Y se disminuye en las siguientes, hasta un valor minima de (A. = 0.04). El nuevo juego de parametres resulta de sumar el 'lector liP al modelo inicial. EI procedimiento asi descrito se repite tres 0 cinco iteraciones. Para controlar la convergencia del proceso se toma como indice el error minima cuadratico (EMC). Se aceptan resultados con EMC < 10%. PROGRAMA DE COMPUTADOR Geomodel es un programa que permite el analisis de datos geoe lectricos tornados en dos dimensiones 0 tomografia de resistividades electricas, cuando se emplean en campo los arreglos Wenner-Schlumberger, Wenner 0 Dipolo-dipolo (Avellaneda y Perez, 2000). El program a lee de un archivo de texto el valor de resistividad aparente tornado en cada lectura horizontal para cada nivel de da- tos, luego aproxima la profundidad de investigaci6n y efectua una interpo-laci6n de los valores de resistividad. La malla de elementos finitos se genera autornaticarnente; se esta- blecen unas dimensiones a la red y los elementos se construyen de forma tal que puedan ser contenidos por el modelo de bloques. 14 La malla esta conforrnada por una red de elemen- tos cuadrilaterales, y se consideran sus limites hasta donde aproximadamente el valor del potencial sea cero. Para efectos de la modelacion, se consideraron cuatro columnas de elementos adicionales a lado y lado de la pseudosecci6n de resistividades, y cuatro filas adicionales para los !imites inferiores. Sin em- bargo, se permite al usuario modificar estos valores para que al final pueda verificar la sensibilidad de los resultados. Seguidamente, el programa vincula de red de ele- mentos con la red de parametres 0 "bloques" que se utiliza en el proceso de inversi6n . Iterativamente se puede observar la respuesta elec- trica del modelo, y la pseudosecci6n de resistividades resultante del proceso de inversion. Se aprecia la con- vergencia del modelo cuando la pseudosecci6n de resistividades calculadas se acerca a la tomada en cam- po, y el error medio cuadratico es menor al 10%. Entrada de datos de resistividad, tipo de arreglo, niveles de datos Modelo de resistividades hornoqeneo Solucionar mediante elementos finitos k'---~ la ecuaci6n de flujo Computar la resistividad aparente calculada, con los potenciales nod ales Calcular la matriz jacobiana y resolver la ecuaci6n(6)para ~P P (i+1) = P(i) +~P 1 EMC < 10% ( Fin) Figura 2. Algoritrno de trabajo GEOFislCA COLOMBIANA, 5, DICIEMBRE DE 2001 MODELACI6N DEL PROBLEMA INVERSO EN GEOELECTRICA 2D MEDIANTE ELEMENTOS FINITOS CASO DE ESTUDIO Se realizaron dos tomografias en zona aledafia al Parque Nacional de la ciudad de Bogota D.C. Ambas se realizaron utilizando un arreglo tipo Wenner; Schlumberger, para una secci6n de 78 m, con 40 puntos de electrodos de corriente y unidad de espaciamiento entre electrodos de 2 m, 17 niveles de datos para un total de 357 datos y bloques en cada una. En las figuras 3 y 4 se muestra la pseudosecci6n de resistividades aparentes medidas para los ejemplos uno y dos; tambien, se muestra la red de elementos finitos empleada. En total se utilizaron 1892 ele- mentos con 200 I nodos en ambos casos. Adicionalmente, la red de elementos fue ampliada en cuatro elementos, a lado y lado de los Ifmites laterales del estudio, y cua- tro filas de elementos en la profundidad maxima de investigaci6n. La selecci6n de esta condici6n de frontera se obtuvo basicarnente al resolver los problemas con una red mas grande que la mostrada en las figuras 3 y 4, Y observar un cambio despreciable en los re- sultados obtenidos. ,.. 1&2 34.9 .., 53] &),1) ,.. 0.7 I! I II I ! I I! 33 '" 21.9 ,.. W D ~ " uo ~ n §~ UI lffl ~_______ = ===== "'~m Figura 3. Pseudosecci6n de resistivida- des aparentes medidas del ejernplo I. Parque Nacional, Bogota D.C. Las dis- tancias en metros .. llJ .a.s " '<> "'., ".J 53,7 72.4 I ~ II II I 11l"n I + I + I + I I I I I , I i OJ 33 ,. .. ". ... ,.. ~ In I~ ~ ~ ~ ro~ ~ I) nw ~_______ = ===== Qhm·m GEOF[SICA COLOMBIANA, 5, OICIEMBRE DE 2001 Figura 4. Pseudosecci6n de resistivida- des aparentes rnedidas del ejernplo 2. Parque Nacional, Bogota D.C. Las dis- tancias en metros 15 BRICENO fT AL. En las figuras 5 y 6 se muestran los modelos de bloques que fueron utilizados para la soluci6n del problema inverso. En ambos casos se ela- boraron tantos bloques como datos del modelo (en total 357 bloques). Adernas, la topologfa de la red de elementos finitos fue construida de forma tal que los elementos pudieran ser contenidos por algun blo- que del modelo de inversi6n. Es necesario anotar que no fue elaborado un adecuado descarte de datos an6malos, y es un factor que se debe estudiar en posteriores in- vestigaciones. Valores muy altos respecto al promedio se puntean en las figuras 5 y 6. En las figuras 7 y 8 se muestran las secciones tomograficas finales para los ejemplos respectivos. Para el primer ejemplo se realizaron Figura 5. Modelo rectangular del subsuelo para el ejernplo I Par- que Nacional, Bogota D.C. Figura 6. Modelo rectangular del subsuelo para el ejemplo 2. Par- que Nacional, Bogota D.C. cuatro iteraciones y se obtuvieron respectivamente los siguientes errores medios cuadraticos (EMC): 17.2,28.1,26.7 Y25.4%. Para el segundo ejemplo tambien se realizaron cuatro iteraciones y se obtuvieron los si- guientes errores medios cuadraticos: 8.6, 10.6,5.0 y 4.7%. En ambos ca- sos se utiliz6 un factor de amortiguamiento (= 0.1). 16 CONCLUSIONES Se presenta un algoritmo de soluci6n al problema inverso utilizando el rnetodo de los elementos finitos como herramienta para encontrar los potenciales medidos en superficie. Como semilla del procedimiento de inversi6n se supone un modelo de resistividades homogeneo, a partir del cual se obtiene la primera matriz de derivadas parciales (matriz jacobiana) que sera corregida iterativamente. GEOFf SICA COLOMBIANA, 5, QICIEMBRE DE 2001 MODELACION DEL PROBLEMA INVERSO EN GEOELECTRICA 20 MEDIANTE ELEMENTOS FINITOS ·Z.5 .. Ift.2 S3,1 ". n... 0.7 3.3 '.0 e.e 1l.J 14.0 t S 8 20 W& 24) SSS 0;9 26.,1 Hi24l-~-=~-=~-==-=-=~-==-=-=~I!l!!::lJ~'~~=-'===-~-=~-=";===-==='--'====''';===____ O",",m r:'33 f" f"113 ]4.0 ,.. fX ~,e ~.3 ~3,O 17 31 54 'n If,Z 291 43'l S