Obtención de imágenes de la distribución de impedancia eléctrica del subsuelo. Aplicación a la detección de objetos locales

Problema inverso 5-1 5. Problema inverso Menke (1989) dice que el problema inverso es simplemente el conjunto de métodos usados para extraer información útil de nuestro entorno a partir de medidas físicas o datos. La información útil vendrá especificada como valores numéricos de alguna propiedad de este entorno. Estas propiedades también se referirán como parámetros del modelo. Se presupone que hay algún método específico (normalmente una teoría matemática o modelo) que relaciona los parámetros con los datos. El problema inverso contrasta con el problema directo, donde se predicen los datos a partir de los parámetros y de un modelo. Normalmente el problema inverso es más difícil de resolver que su correspondiente problema directo. La teoría del problema inverso en su sentido más amplio ha sido desarrollada por los investigadores que trabajan con métodos geofísicos. La razón es que dichos investigadores tratan de entender el interior de la Tierra sólo a partir de datos obtenidos desde la superficie. Sin embargo, el problema inverso aparece en muchas otras ramas de las ciencias físicas, como pueden ser la tomografía médica, el procesado de imagen o el ajuste de curvas. En nuestro caso los parámetros serán las resistividades o conductividades del suelo, los datos serán las tensiones medidas en la superficie y el modelo queda aún por determinar. 5.1. Antecedentes La interpretación de los sondeos eléctricos a fin de determinar las resistividades y espesores de las capas en un medio estratificado ha sido un tema de investigación desde principios de siglo. Hasta la disponibilidad de ordenadores, el intérprete se basaba en los procedimientos de ajuste de curvas. Desde que el problema directo para medios estratificados fue resuelto por medio de la teoría lineal de filtros (Gosh, 1971a, 1971b), han aparecido muchos trabajos que tratan sobre la interpretación automática y numérica (Inman, 1975; Koefoed, 1979; Pous, Marcuello y Queralt, 1987; Zohdy, 1989). En los últimos años ha habido un incremento en el uso de imágenes en dos dimensiones (imágenes que representan una sección transversal del subsuelo). Smith (1986) y Lowry y Shive (1990) se basan en el método de Bristow (Bristow, 1966), una técnica gráfica simple para detectar cavidades en el subsuelo. La intersección de las líneas equipotenciales de las medidas con valores anómalos delimita de una forma cualitativa las cavidades. Los métodos basados en la retroproyección ponderada utilizan un concepto similar y provienen en gran parte de la tomografía médica (Barber y Brown, 1984; Kotre, 1994). Su implementación en aplicaciones geoeléctricas se ha realizado con un cierto éxito (Noel y Xu, 1991; Tsourlos et al., 1993). Estas técnicas intentan reconstruir una sección de resistividad usando una suma ponderada de los potenciales medidos. Generalmente estos métodos son de un solo paso, aunque también se han utilizado versiones iterativas (Yorkey y Webster, 1987; Tsourlos et al., 1993). Problema inverso 5-2 Sin embargo, la mayoría de técnicas utilizan métodos basados en el criterio de mínimos cuadrados. En estos métodos, debido a que el problema está mal condicionado, es necesario aplicar técnicas de regularización (Tikhonov y Arsenin, 1977). Matemáticamente, un problema mal condicionado es aquel en el que pequeñas variaciones (o errores) en los datos provocan grandes variaciones en los parámetros. Ello es debido a que los problemas geofísicos son indeterminados por dos razones: la falta intrínseca de datos y los errores en los datos y en el modelo (Tarantola, 1987). Como resultado, el problema inverso no tiene una solución única, es decir, puede haber más de una solución (conjunto de parámetros) que satisfagan los datos con un error prescrito. La regularización consiste básicamente en introducir alguna clase de información a priori (por ejemplo minimizando la norma de los parámetros) para “estabilizar” el problema. Hua, Webster y Tompkins (1988) comparan tres formas de penalización para regularizar el método de reconstrucción. Sasaki (1992) utiliza el método de Occam (S. C. Constable, R. L. Parker y C.G. Constable, 1987). Loke y Barker (1995) utilizan esta técnica con solo una iteración, a fin de mejorar la velocidad del algoritmo. Los mismos autores (Loke y Barker, 1996a) implementan una versión iterativa basada en un método quasi-Newton que reduce el tiempo de cálculo. Avis y Barber (1994) implementan un método de regularización basado en la técnica SVD (Singular Value Decomposition). Lines y Treitel (1984) revisan diversos métodos de inversión basados en los mínimos cuadrados y comentan la relación entre el método de Marquardt y la técnica SVD. Ellis y Oldenburg (1994) dicen que todas estas técnicas de regularización intentan resolver algo imposible: obtener unos parámetros únicos (aquí los parámetros son las conductividades) a partir de un problema inverso que no tiene una solución única. Por ese motivo aducen que se ha de incluir en el algoritmo toda la información a priori posible. Aunque todas las estructuras geológicas son de tres dimensiones (3D), en realidad las medidas y los algoritmos que tengan en cuenta esta tercera dimensión han sido poco utilizados hasta la fecha. Esto es debido sobre todo a la gran cantidad de medidas que se han de realizar y procesar. Por ello se requiere una instrumentación adecuada y ordenadores con gran capacidad de cálculo. Últimamente, la aparición de sistemas de medida automáticos y de algoritmos más eficientes ha comportado una mayor utilización de imágenes 3D. La mayoría de los algoritmos vistos para el caso 2D se pueden adaptar al caso 3D. Kotre (1996a) utiliza el método de retroproyección ponderada. Oldenburg, McGillivray y Ellis (1993) utilizan una técnica de subespacios de vectores de conductividades para reducir las dimensiones de las matrices. Escogiendo un número reducido de vectores base logra reducir considerablemente el tiempo de cálculo sin deteriorar excesivamente las imágenes. Sasaki (1994) describe dos métodos (uno completo y otro aproximado) que utilizan el método de Occam. Loke y Barker (1996b) adaptan su método (Loke y Barker, 1996a) para el caso 3D. Los métodos iterativos tienen que obtener en cada iteración el modelo (la matriz de sensibilidad o Jacobiana) y los datos usando algún método numérico (elementos finitos, diferencias finitas). A pesar de que algunos autores proponen técnicas especiales (Oldenburg, McGillivray, y Ellis, 1993; Loke y Barker, 1996b), el tiempo de cálculo puede ser considerable, sobre todo en el caso 3D. Además, los algoritmos iterativos pueden tener problemas de convergencia. Los algoritmos de un Problema inverso 5-3 solo paso son un caso particular donde sólo se realiza la primera iteración. Este hecho, junto con la elección de un suelo homogéneo como modelo de partida, reduce considerablemente el tiempo de cálculo (Loke y Barker, 1995). Los mismos autores comentan que las imágenes obtenidas tienen que ser tomadas como una primera estimación de la verdadera distribución de conductividad del subsuelo, que puede ser mejorada utilizando técnicas iterativas. Gasulla, Jordana y Pallás (1999) muestran que esta primera estimación es suficientemente buena para la detección de objetos locales, que es el caso que nos ocupa. Además se aplica un factor de corrección en las imágenes que tiene en cuenta la no-linealidad del problema original. Por lo tanto, todos los algoritmos implementados en este trabajo son de un sólo paso. 5.2. Determinación del modelo El capítulo 2 define cómo realizar medidas de la resistividad aparente del subsuelo. Si el suelo es homogéneo podemos determinar su conductividad σ a partir de una sola medida. El apartado 3.4 describe una primera aproximación al problema inverso al determinar la profundidad y el radio de objetos esféricos y cilíndricos con la calicata Schlumberger cuando el objeto equidista de los electrodos inyectores. Este capítulo aborda el problema inverso de forma más general, al obtener imágenes 2D y 3D de la distribución de conductividad del subsuelo. En los siguientes apartados describimos la teoría matemática en la que se basan los algoritmos de reconstrucción implementados. 5.2.1. Teorema de la Sensibilidad El primer paso para resolver el problema inverso consiste en determinar el modelo (teoría matemática) que relaciona la tensión diferencial en la superficie (datos) con la conductividad del subsuelo (parámetros). El teorema de la sensibilidad proporciona esta relación matemática. Consideremos el medio V limitado por la superficie S de conductividad σ(x,y,z) mostrado en la Figura 5.1, donde se supone que no hay fuentes internas. Sea φ(x,y,z) una distribución de potencial en el medio y Jφ(x,y,z) la densidad de corriente asociada con él. Sea también ψ(x,y,z) una segunda distribución de potencial en el medio y Jψ(x,y,z) la densidad de corriente asociada. Aplicando el teorema de la divergencia tenemos ( ) ( ) ∫ ∫∫ =∇+∇=∇vsvdsJdvJJdvJφφφφ ψψψψ (5.1) donde V es el volumen y S la superficie que encierra al volumen. Como no hay fuentes internas ∇Jφ es cero y ψ es finito en el volumen. Por tanto, la ecuación anterior se reduce a ∫∫∫ ∇∇−=∇=vvsdvdvJdsJψφσψψ φφ (5.2) Similarmente, para Jψ y φ Problema inverso 5-4 ∫∫∫ ∇∇−=∇=vvsdvdvJdsJψφσφφ ψψ (5.3) Vemos que las dos integrales son iguales, por tanto ∫∫∫ ∇∇−== VSS dVdsJdsJψφσφψψφ(5.4) Supongamos que el potencial φ es originado por la aplicación de la corriente Iφ entre los electrodos AB y el potencial ψ es originado por la corriente Iψ aplicada entre los terminales MN (Figura 5.1). En la expresión anterior, la densidad de corriente Jφ será cero excepto en A y B. Similarmente, Jψ será diferente de cero en M y N, por tanto, ∫∇∇== V MNAB dVII ψφσφψψφ (5.5) donde ψAB y φMN son diferencias de tensión. Definimos la impedancia mutua z como ∫ ∇ ∇ === V MNAB dV IIII z ψφφψ σ ψ σ φ σ φ ψ )()( (5.6) expresión que relaciona las diferencias de tensión en la superficie con la conductividad del medio. Como era de esperar se cumple el principio de reciprocidad. MV σ (x,y,z) N AB Iψ Iφ ds S Figura 5.1. Volumen de conductividad σσ(x,y,z) El problema no es lineal debido a que ∇φ y ∇ψ son funciones de la conductividad σ, y en general la estimación de la distribución de conductividad (problema inverso) a partir de la expresión anterior no es posible. Sin embargo, se puede emplear algún método iterativo para resolver el problema (Menke, 1989; Tarantola, 1987). En los métodos iterativos se realiza una hipótesis inicial sobre la distribución de conductividad (normalmente se supone el medio homogéneo), se resuelve el problema directo, se calculan las diferencias respecto a las medidas Problema inverso 5-5 obtenidas y se formula una nueva hipótesis de la distribución de conductividad. El proceso se repite hasta que la diferencia entre las tensiones calculadas (problema directo) y las medidas sea menor que una cota prefijada. Cuando la distribución de conductividad cambia de σ0(x,y,z) a σ0(x,y,z)+∆σ0(x,y,z), el cambio en la impedancia mutua ∆z para los pares de electrodos de corriente (AB) y tensión (MN) es ∫∆+∇∇ ∆−=∆VdV II z ψφ σσψσφ σ)()( 000 0(5.7) Geselowitz (1971) y Lehr (1972) demuestran este resultado conocido como teorema de la sensibilidad. Si expandimos el término ∇ψ(σ0+∆σ0) respecto a ∆σ0, la ecuación anterior se puede expresar como (Murai y Kagawa, 1985) ( ) ∫     ∆+ ∇∇ ∆−=∆V o oo oOdV II z2)()( σ σψσφ σ ψφ (5.8) donde O((∆σ0)2) indica el término de orden más elevado respecto a ∆σ0. Si ∆σ0 es suficientemente pequeño este término puede ser despreciado. Se observa que ahora los gradientes de potencial ∇φ y ∇ψ son independientes de ∆σ0, y por tanto ∆z es lineal respecto a ∆σ0. Si definimos σ como la conductividad actual del medio, la impedancia mutua medida z es ∫∇∇ ∆−=∆+≈∆+== V AB dV II zzzzz I z ψφψ σψσφ σ ψ)()( 00 00000 (5.9) donde z0 es la impedancia mutua debida a una distribución de conductividad σ0. Si no se dispone de información a priori, se suele escoger una distribución de conductividad homogénea como σ0. El proceso iterativo para determinar σ es (Murai y Kagawa, 1985) nnn v nn nnn dV II zzz σσσ σψσφ σ ψφ ∆+= ∇∇ ∆−=−=∆ + ∫ 1 )()( (5.10) donde zn es la impedancia mutua para la distribución de conductividad σn y n = 0, 1, 2,.. es el número de la iteración. En cada iteración n se calcula el cambio de conductividad ∆σn y se obtiene la nueva distribución de conductividad σn+1. El proceso continúa hasta que la diferencia ∆zn es menor que una cierta cota prefijada. El cálculo de zn y de los gradientes ∇φ y ∇ψ en cada iteración se suele realizar mediante métodos numéricos (por ejemplo el Método de los Elementos Finitos). Para calcular ∆σn normalmente se discretiza el medio en celdas o pixeles de conductividad Problema inverso 5-6 constante, si bien también es posible tratar la conductividad como una variable continua (Menke, 1989). 5.2.2. La matriz de sensibilidad La obtención de imágenes de la distribución de conductividad del subsuelo consiste en la determinación de la conductividad de cada una de las celdas en las que se divide. Si el subsuelo está formado por capas horizontales de diferente conductividad y espesor, el medio se discretiza en capas más que en celdas. Este problema lo han tratado diferentes autores (Inman, 1975; Koefoed, 1979; Pous, Marcuello y Queralt, 1987; Zohdy, 1989) y no va a ser estudiado aquí. En el caso 2D se discretiza en celdas cúbicas (o en forma de paralepípedos) la sección transversal a la superficie que está justo debajo de la agrupación de electrodos. La Figura 5.2 muestra un ejemplo con 80 celdas (16 en horizontal por 5 en vertical). En el caso 2D ½ las celdas tienen la forma de una barra infinita en la dirección perpendicular a la sección transversal. Esta forma tan especial de discretizar el subsuelo se utiliza principalmente para la detección de estructuras alargadas en el subsuelo (p.e. tuberías), colocando los electrodos perpendicularmente a dichas anomalías. En el caso 3D se discretizan varias secciones transversales, por lo que los resultados pueden ser más realistas que en el caso 2D. El modelado del subsuelo en 2D tiene la ventaja respecto al 3D de reducir el número de celdas y, por tanto, el tiempo de cálculo de la matriz de sensibilidad. Por el contrario, un modelado 3D da una visión más realista de la distribución de conductividad del suelo. Si sabemos que las estructuras a detectar son extensas en una dirección, de sección constante y paralelas a la superficie, el caso 2D½ puede ser la mejor opción, ya que incorpora las ventajas de los dos modelados anteriores. x y z Figura 5.2. Discretización 2D del subsuelo. El corte vertical se ha dividido en 16 ×× 5 celdas. El incremento de la impedancia mutua en la iteración n se puede expresar ahora como un sumatorio ∑∫ = ∆−=∆Q j j nnn j ndvLLz 1j celda )()( ψφσ(5.11) Problema inverso 5-7 donde Q es el número de celdas, y para simplificar la nomenclatura se han definido los vectores ψ φ σψ ψ σφ φ I L I L n n n n )( )( )( )( ∇ = ∇ = (5.12) La expresión (5.11) se refiere a una posición de los 4 electrodos ABMN. En la práctica se realizan diferentes medidas, cada una de ellas con diferentes situaciones de los electrodos sobre la superficie. El conjunto de todas las colocaciones par inyector - par detector que utilizamos es lo que se denomina configuración electródica (ver capítulo 2), y para referirnos a una medida en concreto (posición de 4 electrodos) utilizamos un único subíndice, en este caso el i: PidvLLzQ j j n i n i n j n i..1 )()( 1j celda =∆=∆∑∫ = ψφσ(5.13) donde P es el número de medidas. El sistema de ecuaciones, en forma matricial es el siguiente: nnn SZ∆Σ=∆(5.14) donde ∆Zn y ∆Σn son vectores de dimensiones P × 1 y Q × 1 respectivamente, y Sn es una matriz de dimensiones P × Q. Ésta última recibe el nombre de matriz de sensibilidad porque contiene las sensibilidades de todas las medidas con respecto a todas las celdas. Los elementos de Sn son ∫ −= j celda )()( j n i n i n ij dvLLsψφ(5.15) La estimación de la distribución de conductividad se ha convertido ahora en encontrar la conductividad asignada a cada elemento en el que hemos dividido el subsuelo. El proceso iterativo dado por (5.10) queda ahora como ...Q,j Piszzz n j n j n j Q j n jij n ii n i 21 ...2,1 1 1 =∆+= =∆=−=∆ + = ∑ σσσ σ(5.16) o en forma matricial como nnn nnn SZ ∆Σ + Σ = Σ ∆Σ=∆ +1(5.17) Problema inverso 5-8 La diferencia de los distintos métodos existentes para estimar la distribución de conductividad radica básicamente en la manera de calcular ∆Σn, como veremos al describir los diferentes algoritmos implementados. 5.3. Métodos de un solo paso Un caso particular de los métodos iterativos son los algoritmos de un solo paso, en los que sólo se realiza la primera iteración. El algoritmo expresado en forma matricial se reduce en este caso a ∆Σ+Σ=Σ ∆Σ = ∆ 0est SZ(5.18) donde Σest es el vector de conductividad estimada. Si se escoge la conductividad inicial Σ0 como una distribución homogénea, la impedancia mutua correspondiente Z0 se puede calcular analíticamente (capítulo 2). Definimos r1 a r4 como las distancias de los electrodos ABMN a un punto arbitrario P (Figura 5.3). r3r4 r1 P r 2 MNB A Figura 5.3. Distancias de los electrodos ABMN a un punto P del subsuelo El vector L(φ) debido a los electrodos de inyección A y B es           −= 3 2 2 3 1 1 2 )( r r r r Lr r r r π ρ φ(5.19) y el debido a los electrodos detectores M y N es           −= 3 4 4 3 3 3 2 )( r r r r Lr r r r π ρ ψ(5.20) Para calcular el elemento sij de la matriz de sensibilidad debemos integrar en el pixel j el producto escalar de estos vectores debido a la disposición de los electrodos ABMN cuando realizamos la medida i. Sin embargo, esta integral no tiene en general expresión analítica. Por lo Problema inverso 5-9 tanto será necesario integrar numéricamente esta función escalar. Entre los muchos métodos de integración numérica existentes, se ha elegido el método de Gauss. Para la mayoría de funciones este método da resultados más exactos que otros métodos comúnmente usados, como el método de los rectángulos o el de Romberg (Loke y Barker, 1995). La exactitud y el tiempo de cálculo de la integral crecen con el número de puntos que se escojan en cada celda para evaluar la integral. Las celdas más cercanas a los electrodos necesitan de mayor número de puntos ya que los campos eléctricos varían más rápidamente. Si solo se escoge un punto por celda para evaluar la integral, es equivalente a multiplicar el producto escalar de los campos normalizados por el volumen de la celda. En el caso 2D½ antes de la integración numérica se realiza una integración analítica en la dirección en la que la celda es de dimensión infinita (Loke y Barker, 1995). Cano (1999) describe con detalle cómo se realiza la integración numérica en los casos 2D, 3D y 2D ½. Para acelerar el proceso de reconstrucción la matriz de sensibilidad puede ser calculada con anterioridad y almacenada en un fichero. Los algoritmos implementados que se describen en los siguientes apartados tienen en común que son de un solo paso y que estiman la conductividad ponderando las medidas obtenidas por los coeficientes de sensibilidad. La diferencia entre ellos es como se realiza esta ponderación. 5.3.1. Métodos de inversión basados en el criterio de mínimos cuadrados Puesto en forma matricial, el problema es determinar la distribución de conductividades a partir del sistema de ecuaciones ∆Σ = ∆ SZ(5.21) Si la matriz S es cuadrada (P = Q) la solución se obtiene simplemente invirtiendo la matriz de sensibilidad ZS∆=∆Σ −1(5.22) Pero en general no tiene por qué haber el mismo número de medidas que celdas de conductividad. Si P > Q el sistema es sobredeterminado y no tiene por qué existir una solución exacta. Si definimos el vector error de predicción o residuo como ∆Σ − ∆ = SZe(5.23) y su norma L2 (también denominada norma cuadrática o euclídea) como )()( ∆Σ−∆∆Σ−∆== SZSZeeETT (5.24) Problema inverso 5-16 A B M N Figura 5.5. Líneas equipotenciales formadas por los electrodos MN La conductividad se determina según (5.40) y los coeficientes de sensibilidad se calculan como    =contrario casoen 0 ialesequipotenc lineas las entre está celda la si ij ij s s Este método será denominado como retroproyección equipotencial. Noel y Xu (1991) proponen un método similar, pero además ponderan la sensibilidad por la parte proporcional de la celda que está entre las líneas equipotenciales. Tsourlos et al. (1993) proponen una solución iterativa. La idea de retroproyectar entre líneas equipotenciales fue sugerida originalmente por Barber, Brown y Freston (1983) para el caso de la tomografía médica en analogía al caso de la tomografía computerizada por rayos X. Barber y Brown (1984) presentan imágenes in vivo de la resistividad de los tejidos. El algoritmo que implementan, que recibe el nombre de Applied potential tomography, utiliza unos pesos calculados geométricamente y realiza un filtrado frecuencioespacial en un espacio transformado para mejorar la apariencia de las imágenes (Barber y Brown, 1986; Barber y Seagar, 1987a). Barber y Seagar (1987b) describen el sistema de medida utilizado. Powell, Barber y Freeston (1987) adaptan el método utilizando una agrupación lineal de electrodos. 5.4. Imágenes con datos sintéticos Los datos utilizados en los algoritmos de reconstrucción pueden obtenerse de tres fuentes diferentes: a) a partir de medidas experimentales, b) a partir de soluciones numéricas (por ejemplo, el método de los elementos finitos o el método de los elementos de contorno), c) a partir de una solución analítica. La primera alternativa requiere realizar medidas experimentales de campo o sobre un modelo analógico (por ejemplo una cubeta de plástico llena de agua y donde se introducen objetos de conductividad diferente a la del agua). Además, se debe disponer de la instrumentación adecuada. Otros inconvenientes son que el tiempo necesario para realizar las medidas puede que sea largo, y que las medidas tienen asociado un error (debido a la instrumentación, a la posición de los electrodos, etc.). La segunda alternativa supone emplear alguno de los métodos numéricos Problema inverso 5-17 existentes para resolver el problema directo. El problema con estos métodos radica en que el tiempo de cálculo puede ser elevado, sobre todos si se requiere que la solución sea precisa. Por el contrario, una solución analítica (tercera alternativa) puede ser rápida y suficientemente precisa. Sin embargo, sólo disponemos de solución analítica para ciertos casos. En el capítulo 3 se describieron diferentes soluciones para el caso de una esfera y un cilindro inmersos en un suelo homogéneo. Para obtener los datos utilizaremos la expresión (3.3) para el caso de la esfera y la expresión (3.29) para el caso del cilindro. En el capítulo 6 se realizarán medidas experimentales para validar los algoritmos utilizados. Las configuraciones electródicas utilizadas son la Schlumberger y la doble dipolo (Capítulo 2). La configuración Schlumberger presenta un buen compromiso entre la visibilidad y el margen dinámico necesario del detector (apartado 3.3). La configuración doble dipolo es la que tiene mayor visibilidad. Sin embargo, como se comentó en el Capítulo 4, se ve más afectada por los errores en las medidas. Las imágenes obtenidas por los diferentes métodos muestran la distribución de conductividad absoluta. En todos los casos se supone un suelo homogéneo de conductividad 1 S/m y una inyección de corriente de 1 A. Los coeficientes de sensibilidad utilizados por los algoritmos se obtienen a partir de (5.15). A no ser que se diga lo contrario, la integral se calcula en 64 puntos equiespaciados. En el apartado 5.4.1 se comparan los diferentes algoritmos de reconstrucción utilizando los datos sintéticos para el caso de la esfera. El apartado 5.4.2 estudia las características de la matriz de sensibilidad a partir de la descomposición SVD. El apartado 5.4.3 muestra el efecto del parámetro λ en los métodos Marquardt y Occam. En el apartado 5.4.4 se obtiene un factor de escalado de las imágenes que tiene en cuenta la no-linealidad del problema inverso. El apartado 5.4.5 y 5.4.6 estudian el efecto del modelo utilizado y el error en las medidas sobre las imágenes obtenidas. El apartado 5.4.7 considera el caso de un objeto cilíndrico perpendicular a la agrupación de electrodos. Por último, el apartado 5.4.8 muestra que si el objeto no se encuentra justo debajo de la agrupación de electrodos es necesario utilizar imágenes 3D. 5.4.1. Comparación de los diferentes algoritmos de reconstrucción En este apartado comparamos cualitativamente y cuantitativamente las imágenes 2D obtenidas por cinco de los algoritmos descritos en el apartado 5.3: Marquardt (5.33), Occam (5.32), TSVD, retroproyección total y retroproyección equipotencial. En el método de Occam se ha utilizado como matriz L una aproximación discreta de un operador derivada segunda (Sasaki, 1992). El apartado 5.4.3 muestra que un valor adecuado para λ es de diez a cien veces el obtenido con el método de la curva L. En el método TSVD el rango r de la matriz inversa en (5.30) se determina a partir de la curva L discreta. En los métodos de retroproyección no se aplica ningún tipo de filtrado en el dominio de la frecuencia y se escoge experimentalmente un factor de amplificación k = 10. Las imágenes 2D corresponden a un corte transversal a la superficie justo debajo de la agrupación de electrodos. La Figura 5.6 muestra un corte dividido en 85 celdas (17 en horizontal y 5 en vertical) Problema inverso 5-18 de dimensión una unidad en las tres direcciones del espacio. La unidad espacial se escoge de 1 m, aunque su valor puede ser cualquiera sin afectar a las imágenes que se presentarán en este apartado. El número de electrodos es de 16 separados entre ellos una unidad. La posición del electrodo 1 es x = -7,5, y = 0, z = 0 y la del electrodo 16 es x = 7,5, y = 0, z = 0. El origen de coordenadas está situado entre los electrodos 8 y 9. x yz 1 2 3 5 6 7 8 9 10 11 12 1314 15 164 1 2 17 18 34 51 68 85 35 52 69 Figura 5.6. Discretización 2D del subsuelo en 85 celdas (17 en horizontal y 5 en vertical) del corte transversal justo debajo de la agrupación de electrodos. La Figura 5.7 muestra las imágenes obtenidas con los diferentes métodos utilizando la configuración Schlumberger y doble dipolo en el caso de una esfera conductora de radio 0,5 unidades y profundidad 1,5 situada en x = y = 0 (entre los electrodos 8 y 9). La Figura 5.8 muestra la curva L para los métodos Marquardt, TSVD y Occam. El valor de λ escogido para el método de Marquardt y Occam es λ = 10c. Para el método TSVD se representa la variación discreta de la curva L en función del rango r de la matriz inversa. En este caso se ha escogido r = 35. La Figura 5.9 muestra las imágenes obtenidas para la misma esfera conductora situada en x = y = 0, z = 2,5 (λ = 100c, r = 50). a) b) Problema inverso 5-19 c) d) e) Figura 5.7. Imagen 2D obtenida con la configuración Schlumberger (izquierda) y doble dipolo (derecha) a partir de datos sintéticos para una esfera conductora de radio 0,5 unidades situada en x = y = 0, z = 1,5. El modelo utilizado es de 17 ×× 1 ×× 5 celdas de lado unidad. a) Marquardt (λλ = 1,50××10-8; 1,83××10-8), b) TSVD (r = 35), c) Occam (λλ=1,51××10-9; 2,21××10-9), d) retroproyección total, e) retroproyección equipotencial. a) Problema inverso 5-20 b) c) Figura 5.8. Curva L obtenida con la configuración Schlumberger (izquierda) y doble dipolo (derecha) a partir de datos sintéticos para una esfera conductora de radio 0,5 unidades situada en x = y = 0, z = 1,5. a) Marquardt, b) TSVD, c) Occam. a) b) Problema inverso 5-21 c) d) e) Figura 5.9. Imagen 2D obtenida con la configuración Schlumberger (izquierda) y doble dipolo (derecha) a partir de datos sintéticos para una esfera conductora de radio 0,5 unidades situada en x = y = 0, z = 2,5. El modelo utilizado es de 17 ×× 1 ×× 5 celdas de lado unidad. a) Marquardt (λλ = 1,86××10-12; 2,85××10-12), b) TSVD (r = 50), c) Occam (λλ=1,86××10-12; 3,07××10-12), d) retroproyección total, e) retroproyección equipotencial. Los métodos Marquardt, TSVD y Occam (basados en la minimización por mínimos cuadrados) son claramente superiores a los métodos de retroproyección. En éstos la imagen obtenida presenta poca resolución y una gran dependencia de la configuración electródica utilizada. Además, el cambio de conductividad estimada decrece con la profundidad. Gasulla, Jordana y Pallás (1998b) corroboran con medidas experimentales la superioridad del método de Marquardt frente al de retroproyección total. Kotre (1994,1996a) utiliza filtros en el dominio de la frecuencia para mejorar las imágenes, pero éstos también podrían ser utilizados en los métodos de mínimos cuadrados. El método de Occam “suaviza” la imagen de la esfera alargándola en sentido vertical. En cambio los métodos Marquardt y TSVD tienen un comportamiento muy similar y localizan correctamente la esfera. El valor de λ utilizado en los métodos Marquardt y Occam disminuyen con la profundidad de la esfera. En el método TSVD el rango r aumenta. Problema inverso 5-22 Para estimar la bondad de las imágenes en términos cuantitativos definimos el error de conductividad normalizado como ( ) ( ) ∑ =       −= Q iideal i ideal i est i est i maxmax Q ECN 1 2 1 σ σ σ σ(5.43) donde σiest es la conductividad estimada por los algoritmos. σiideal = 1 para i = 26 cuando z = 1,5; i = 43 cuando z = 2,5; y σiideal = 0 para el resto de celdas. La Tabla 5.1 muestra el valor de ECN para los casos analizados anteriormente. Los resultados confirman las conclusiones anteriores. Los métodos de retroproyección tienen errores elevados en todos los casos. Los métodos Marquardt y TSVD presentan los errores menores. El método de Marquardt será el utilizado preferentemente. El método TSVD presenta el inconveniente de que la determinación del rango r de la matriz inversa se ha de realizar por prueba y error, al no disponer de la Toolbox Spline de Matlab necesaria para determinarlo de forma automática (Hansen, 1992). Schlumberger Doble dipolo x = y = 0, z = 1,5 x = y = 0, z = 2,5 x = y = 0, z = 1,5 x = y = 0, z = 2,5 Marquardt 2,40×10-2 2,91×10-2 2,08×10-2 3,00×10-2 TSVD 2,27×10-2 3,58×10-2 1,43×10-2 2,80×10-2 Occam 0,105 6,12×10-2 7,39×10-2 6,25×10-2 Retro total 0,226 0,373 0,248 0,308 Retro equipot 0,302 0,349 0,234 0,282 Tabla 5.1. ECN para los diferentes métodos de reconstrucción. La Figura 5.10 muestra la imagen 2D obtenida con la configuración Schlumberger y el método Marquardt para una esfera conductora situada en x = y = 0, z = 3,5 y en x = 4, y = 0, z = 1,5. La localización de la esfera es correcta en ambos casos. a) b) Figura 5.10. Imagen 2D obtenida con la configuración Schlumberger y el método de Marquardt para una esfera conductora de radio 0,5 unidades situada en a) x = y = 0, z = 3,5 (λλ = 1,92××10-16), b) x = 4, y = 0, z = 1,5 (λλ = 2,45××10-8). Problema inverso 5-23 5.4.2. Estudio de las matrices de sensibilidad La descomposición en valores propios (SVD) dada por (5.27), revela las dificultades asociadas con el mal acondicionamiento de la matriz S. La Figura 5.11 muestra el logaritmo de los valores propios normalizados al valor mayor para la configuración Schlumberger y doble dipolo con un modelo de 17 × 1 × 5 celdas de lado unidad (Figura 5.6). Los valores propios decaen gradualmente a cero, lo que es una característica de los problemas mal definidos. La configuración doble dipolo está mejor condicionada (número de condicionamiento = 8,27×109) que la configuración Schlumberger (número de condicionamiento = 3,63×1010). Valores propios SVD matriz S Valores propios λi 0 20 40 60 80 Valores propios normalizados log10(λj/λ1) -10 -8 -6 -4 -2 0 doble diplo Schlumberger Figura 5.11. Representación de los valores propios normalizados de la matriz S con un modelo de 17 ×× 1 ×× 5 celdas de lado unidad para la configuración Schlumberger y doble dipolo. Las columnas de la matriz V en (5.27) son los vectores propios del espacio de parámetros (Menke, 1989). La Figura 5.12 muestra los vectores asociados con los valores propios λ5, λ20, λ40, λ60, λ80, para las configuraciones Schlumberger (izquierda) y doble dipolo (derecha). Las celdas más superficiales están asociadas con los valores propios mayores mientras que las celdas inferiores se corresponden con los valores propios menores. Por lo tanto cuando regularizamos el problema, al eliminar o atenuar los valores propios menores estamos eliminando información acerca de las celdas inferiores. Estas celdas contribuyen muy poco a la diferencia de potencial medida en la superficie y, por tanto, el ruido presente en las medidas puede ser comparable o superior a la contribución de estas celdas. Es lógico, pues, que si intentamos obtener información acerca de las celdas más profundas (regularizando poco el problema) el ruido domine la solución. Problema inverso 5-24 a) b) c) d) e) Figura 5.12. Imagen de los vectores propios de la matriz V (espacio de parámetros) para la configuración Schlumberger (izquierda) y doble dipolo (derecha) con el modelo de la Figura 5.6. Imagen asociada con el valor propio a) λλ5, b) λλ20, c) λλ40, d) λλ60, e) λλ80. Problema inverso 5-25 Otra forma de evaluar la matriz de sensibilidad es mediante la matriz de resolución de los parámetros R, que caracteriza cómo los parámetros estimados se aproximan a la solución verdadera. Imaginemos que existe un conjunto de parámetros ∆Σreal que solucionan el problema ∆Z = S∆Σreal. Sustituyendo esta expresión en (5.30), la conductividad estimada es ( ) realrealest RSSZS∆Σ=∆Σ=∆=∆Σ ++ (5.44) Si R = I, entonces ∆Σest = ∆Σreal y todos los parámetros se obtienen de forma única. Por ejemplo, en un problema sobredeterminado S+ = (STS)-1ST, (5.25), y R = I. Si R no es la matriz identidad, entonces los parámetros estimados son una contribución ponderada de los parámetros reales. Por ejemplo, en un problema indeterminado S+ = S T(SST)-1, (5.26), y R = ST(SST)-1S. La matriz R depende de la configuración electródica, del método de regularización y del modelo utilizado. La Figura 5.13 muestra R utilizando la configuración doble dipolo con el modelo de la Figura 5.6 y S+ = (STS)-1ST. En este caso R = I. Figura 5.13. Matriz de resolución R para la configuración doble dipolo con el modelo de la Figura 5.6 y S+ = (STS)-1ST. En este caso R = I. La covarianza de los parámetros depende de la covarianza de los datos y del modo en que los errores en los datos se propagan a los parámetros. Menke (1989) define la matriz de covarianza de los parámetros estimados, que caracteriza el grado de amplificación del error, como [ ] T n est SCS++ =∆Σcov (5.45) donde Cn es una matriz de covarianza de los datos. Si consideramos ICnn 2 σ=, [ ] ( ) 1 2 cov − =ΣSST n est σ y sustituimos S por (5.28), [ ] T rprn est VV 2 2 cov − Λ=Σσ. Así pues, la covarianza de los parámetros estimados es muy sensible a los valores propios menores y por tanto Problema inverso 5-32 un contraste resistivo (σesf/σh) menor que 0,01 o mayor que 100 equivale a que el objeto sea en la práctica totalmente aislante o conductor respectivamente. En lo que sigue se aplicará el factor de corrección en las imágenes de conductividad obtenidas. 5.4.5. Influencia del modelo Cuando el modelo no se adapta bien al objeto pueden producirse resultados extraños. La Figura 5.22 muestra la imagen obtenida con la configuración Schlumberger y el método de Marquardt para una esfera aislante de radio 0,7 a una profundidad de 1,5 situada en x = y = 0. Los coeficientes de sensibilidad se han calculado con 1 punto por celda. El resultado es irreal ya que la conductividad estimada de la esfera es de –0,2. Esto se debe a que la esfera tiene un volumen mayor que la celda y el algoritmo aumenta el cambio de la conductividad para compensar este efecto. Si el modelo es de celdas de lado 1,4 unidades y el centro de la primera fila coincide con el centro de la esfera los resultados con coherentes (Figura 5.23). Figura 5.22. Imagen 2D obtenida con la configuración Schlumberger y el método de Marquardt para una esfera aislante de radio 0,7 situada en x = y = 0, z = 1,5. La conductividad de la esfera es negativa. Se ha utilizando 1 punto por celda para obtener los coeficientes de sensibilidad. Figura 5.23. Imagen 2D obtenida con la configuración Schlumberger y el método de Marquardt para una esfera aislante de radio 0,7 situada en x = y = 0, z = 1,5. Las celdas son de lado 1,4 unidades y la primera fila está centrada en z = 1,5. Se ha utilizado 1 punto por celda para obtener los coeficientes de sensibilidad. Problema inverso 5-33 La Figura 5.24 muestra la imagen obtenida con un modelo de 16 × 1 × 5 celdas de lado unidad. La esfera se atribuye ahora a dos celdas. En este caso, se han utilizado 8 puntos por celda para calcular los coeficientes de sensibilidad. Figura 5.24. Imagen 2D obtenida con la configuración Schlumberger para una esfera aislante de radio 0,7 situada en x = y = 0, z = 1,5. Se han utilizando 8 puntos por celda para obtener los coeficientes de sensibilidad. La Figura 5.25 representa las imágenes obtenidas para una esfera aislante de radio unidad situada en x = y = 0, z = 2, utilizando dos modelos diferentes: uno de 17 × 1 × 5 y otro de 16 × 1 × 5 celdas de lado unidad. Este segundo modelo delimita con más precisión la anomalía. En ambos casos aparecen conductividades negativas debido a que las celdas son de 1 unidad en la dirección y mientras que el diámetro de la esfera es de 2 unidades. La sensibilidad de cada celda se ha obtenido utilizando 1 punto. a) b) Figura 5.25. Imagen 2D obtenida con la configuración Schlumberger para una esfera aislante de radio unidad situada en x = y = 0, z = 2. a) Modelo 17 ×× 1 ×× 5 celdas. b) Modelo 16 ×× 1 ×× 5 celdas. Las celdas son de lado unidad. Se ha utilizado 1 punto por celda para obtener los coeficientes de sensibilidad. La Figura 5.26 muestra las imágenes obtenidas con celdas de dimensión 1 × 2 × 1. Los valores negativos de conductividad desaparecen debido a que el modelo se adapta mejor a las dimensiones de la esfera. Problema inverso 5-34 a) b) Figura 5.26. Imagen 2D obtenida con la configuración Schlumberger para una esfera aislante de radio unidad situada en x = y = 0, z = 2. a) Modelo 17 ×× 1 ×× 5 celdas. b) Modelo 16 ×× 1 ×× 5 celdas (λλ = 1,51××10-8). Las celdas son de dimensión 1 ×× 2 ×× 1 unidades. La Figura 5.27 muestra la imagen 2D obtenida con la configuración Schlumberger con un modelo de 9 × 1 × 3 celdas de lado 2 unidades. La sensibilidad de cada celda se ha obtenido utilizando 8 puntos. La esfera se localiza correctamente. Figura 5.27. Imagen 2D obtenida con la configuración Schlumberger para una esfera aislante de radio unidad situada en x = y = 0, z = 2. El modelo es de 9 ×× 1 ×× 3 celdas de lado 2 unidades. La sensibilidad de las celdas se ha obtenido utilizando 8 puntos. La Figura 5.28 muestra la imagen obtenida y la curva L con la configuración doble dipolo si utilizamos un modelo de 32 × 1 × 10 celdas de dimensión 0,5 × 2 × 0,5 unidades. La anomalía no queda bien definida debido a que el valor de λ encontrado no es el adecuado. La causa es que el número de celdas es muy superior al número de medidas disponibles. Con un valor de λ = 10-17 los resultados mejoran ostensiblemente (Figura 5.29). Este valor se ha obtenido por prueba y error comparando la imagen obtenida con la imagen “ideal”. Problema inverso 5-35 a) b) Figura 5.28. Distribución de conductividad y curva L con la configuración doble dipolo para una esfera aislante de radio unidad a una profundidad de 2 unidades usando el método Marquardt (λλ = 4,18××10-6). El modelo es de 32 ×× 1 ×× 10 celdas de lado 0,5 unidades en las direcciones x, z y de 2 unidades en la dirección y. La esfera está situada en x = 0, y = 0. Figura 5.29. Imagen 2D obtenida con la configuración doble dipolo para una esfera aislante de radio unidad situada en x = y = 0, z = 2 (λλ =10-17). El valor de λλ se ha encontrado por prueba y error. El modelo es de 32 ×× 1 ×× 10 celdas de dimensión 0,5 ×× 2 ×× 0,5 unidades. Sin embargo, en un caso real no sabremos cuál es la imagen ideal (si lo supiéramos no haría falta obtenerla), por lo que el procedimiento anterior no es adecuado. Lo que sí podemos tener es cierta información a priori de la distribución de conductividad. Esta información previa puede tener origen en un conocimiento previo de la zona o a través de datos obtenidos con otros métodos geofísicos. En nuestro caso también puede obtenerse con los algoritmos de retroproyección, que no precisan estimar ningún parámetro λ. La Figura 5.30 muestra la imagen obtenida con el método de retroproyección entre líneas equipotenciales. La imagen no es muy buena pero puede servir como una primera estimación para delimitar la zona donde se encuentra el objeto. Problema inverso 5-36 Figura 5.30. Imagen 2D obtenida con la configuración doble dipolo y el método retroproyección equipotencial para una esfera aislante de radio unidad situada en x = y = 0, z = 2. El modelo es de 32 ×× 1 ×× 10 celdas de dimensión 0,5 ×× 2 ×× 0,5 unidades. La Figura 5.30 indica que el objeto se encuentra entre x = -2,5 y x = 2,5. La Figura 5.31 muestra la curva L y la imagen obtenida con el método de Marquardt cuando el modelo sólo incluye las celdas probables de contener el objeto. Esto sería equivalente a considerar la matriz de covarianza ICmm 2 σ=en (5.37), con 1= m σ para las celdas situadas entre x = -2,5 y x = 2,5; y 0= m σ(no permitimos que varíe su valor) para el resto de celdas. La curva L queda bien definida en este caso y se utiliza un valor de λ = 100c. Figura 5.31. Imagen 2D obtenida con la configuración doble dipolo para una esfera aislante de radio unidad situada en x = y = 0, z = 2 (λλ =7,27××10-16). El modelo es de 10 ×× 1 ×× 10 celdas de dimensión 0,5 ×× 2 ×× 0,5 unidades. La curva L queda bien definida. Una vez localizado el objeto podemos discretizar con más detalle la zona donde se encuentra (de hecho volvemos a introducir más información previa). La Figura 5.32 muestra la imagen obtenida y la curva L con la configuración doble dipolo usando un modelo de 10 × 1 × 10 celdas de dimensión 0,25 × 2 × 0,25 centrado en x = y = 0, z = 2. Problema inverso 5-37 Figura 5.32. Imagen 2D obtenida con la configuración doble dipolo para una esfera aislante de radio unidad situada en x = y = 0, z = 2 (λλ =1,23××10-21). El modelo es de 10 ×× 1 ×× 10 celdas de dimensión 0,25 ×× 2 ×× 0,25 unidades. La curva L queda bien definida. Otra alternativa para obtener el valor de λ adecuado en el caso de la Figura 5.28 es incrementar el número de medidas (o bien reducir el número de celdas). La Figura 5.29 muestra la imagen 2D obtenida utilizando 21 electrodos separados 0,75 unidades (189 medidas independientes). La posición del electrodo 1 es x = -7,5 y la del electrodo 16 es x = 7,5. La curva L queda bien definida. Figura 5.33. Imagen 2D y curva L obtenida utilizando 21 electrodos separados 0,75 unidades con la configuración doble dipolo para una esfera aislante de radio unidad situada en x = y = 0, z = 2 (λλ = 5,62××10-15), resultando 189 medidas. El modelo es de 32 ×× 1 ×× 10 celdas de dimensión 0,5 ×× 2 ×× 0,5 unidades. Si no es posible incrementar el número de medidas reales, podemos interpolar las medidas disponibles. La Figura 5.34 muestra la imagen obtenida con la configuración doble dipolo si interpolamos las 104 medidas (obtenidas utilizando 16 electrodos) con la función spline de Matlab, resultando 194 medidas. La curva L asociada queda bien definida. Problema inverso 5-38 Figura 5.34. Distribución de conductividad y curva L con la configuración doble dipolo para una esfera aislante de radio unidad situada en x = y = 0, z =2 (λλ = 7,07××10-9). El modelo es de 32 ×× 1 ×× 10 celdas de dimensión 0,5 ×× 2 ×× 0,5 unidades. Las medidas se han obtenido interpolando las 104 medidas existentes con la función spline de Matlab, resultando 194 medidas 5.4.6. Errores en las medidas Los errores en las medidas pueden ser debidos, entre otros factores, a la inexactitud en la posición de los electrodos, al ruido geoeléctrico y a la inexactitud del sistema de medida. En el laboratorio, el ruido “geoeléctrico” vendrá dado por la variación de la resistividad del agua con la temperatura ambiente y con las pequeñas vibraciones mecánicas de la cubeta (al vibrar la estructura del edificio). Las medidas de laboratorio se realizan con 16 electrodos separados 2 cm (1 unidad). La Figura 5.35 muestra las imágenes 2D obtenidas para una esfera aislante situada en x = y = 0, z = 2, 3, 4 unidades, con la configuración Schlumberger (izquierda) y doble dipolo (derecha). Se han utilizado datos sintéticos sin ruido. El modelo es de 16 × 1 × 5 celdas de dimensión 1 × 2 × 1 unidades. La esfera se localiza correctamente en todos los casos. El parámetro λ disminuye al aumentar la profundidad de la esfera. a) b) Problema inverso 5-39 c) d) e) f) Figura 5.35. Imagen 2D obtenida con la configuración Schlumberger (izquierda) y doble dipolo (derecha) a partir de datos sintéticos sin ruido para una esfera conductora de radio unidad situada en x = y = 0, a) z = 2, λλ = 1,18××10-4; b) z = 2, λλ = 5,75××10-5; c) z = 3, λλ = 2,28××10-8; d) z = 3, λλ = 7,73××10-9; e) z = 4, λλ = 7,68××10-12; f) z = 4, λλ = 9,65××10-12. El modelo utilizado es de 16 ×× 1 ×× 5 celdas de dimensión 1 ×× 2 ×× 1 unidades. Añadimos ruido gaussiano a las medidas. Definimos la relación S/N para una medida j como         =       j j j z N S σ 0 (5.51) donde zj0 es la impedancia mutua de la medida j y j σ es la desviación típica del ruido gaussiano para la medida j. En el apartado 4.3.1 vimos que para la configuración Schlumberger Vdmin = 35,7 mVpp y para la configuración doble dipolo Vdmin = 0,777 mVpp. Basados en las medidas experimentales realizadas en el laboratorio, la amplitud del ruido añadido se ajusta para que (S/N)min = 70 dB en la configuración Schlumberger y (S/N)min = 37 dB para la configuración doble dipolo, que equivale a unos 10 µV de ruido. La Figura 5.36 muestra las imágenes obtenidas con la esfera situada a una profundidad de 4 unidades. Las imágenes se deterioran respecto a las de la Figura 5.35, pero la esfera aún se localiza correctamente con las dos configuraciones. El parámetro λ aumenta para atenuar el efecto del ruido. Las imágenes a profundidad 2 y 3 unidades no se muestran ya que prácticamente no cambian. Problema inverso 5-40 a) b) Figura 5.36. Imagen 2D obtenida con la configuración a) Schlumberger (λλ = 1,68××10-7) y b) doble dipolo (λλ = 1,04××10-7) a partir de datos sintéticos a los que se ha añadido ruido gaussiano. La esfera aislante está situada en x = y = 0, z = 4. La Figura 5.37 muestra las imágenes obtenidas con la esfera en z = 3, 4, cuando la relación S/N es del 0,1 % para todas las medidas y en las dos configuraciones. Las imágenes se degradan bastante para z = 4, y la localización de la esfera es incorrecta. a) b) c) d) Figura 5.37. Imagen 2D obtenida con la configuración Schlumberger (izquierda) y doble dipolo (derecha) obtenida a partir de datos sintéticos a los que se ha añadido ruido gaussiano con una relación S/N = 0,1 % para todas las medidas. La esfera está situada en x = y = 0. a) z = 3 (λλ = 3,04××10-5), b) z = 3 (λλ = 9,66××10-7), c) z = 4 (λλ = 7,83××10-5), d) z = 4 (λλ = 1,51××10-5). Si disponemos de información a priori acerca del error en las medidas conviene tenerla en cuenta. Partiendo de (5.36) Problema inverso 5-41 ( ) ZCSCSCSn T mn T∆+=∆Σ − − −− 1 1 11 consideramos ICmm 2 σ=, y Cn = ( ) 22 1,.., P diag σσ con P el número de medidas. En el ejemplo de la Figura 5.36 la desviación típica es constante para todas las medidas y (5.36) se reduce al método de Marquardt, con ( ) 2 mnσσλ=. En el ejemplo de la Figura 5.37 la desviación típica es proporcional a las medidas y (5.36) se puede reescribir como ( ) rel T relrel T relrel ZSISS ∆+=∆Σ −1 λ(5.52) donde el subíndice rel indica relativo. El parámetro λ y los coeficientes de los vectores ∆Σrel, ∆Zrel y de la matriz Srel son 2 2 0 0 0 , 0 , 0 , ..1 1 −               = = = ∆ =∆ = ∆ =∆ N S s z s Qi z z z ..P j m ij i relij j i reli j relj σ σ λ σ σ σ σ (5.53) donde σ0 es la conductividad del medio y λ es inversamente proporcional al cuadrado de la relación S/N. Es decir, cuanto menor sea el ruido menor será el valor de λ utilizado. La expresión (5.52) es el método de Marquardt con la matriz de sensibilidad, los datos y los parámetros expresados de forma relativa. La Figura 5.38 muestra las imágenes obtenidas en el caso de la Figura 5.37 utilizando la expresión (5.52). La configuración doble dipolo permite, ahora, la correcta localización de la esfera situada a profundidad 4 unidades, mientras que la configuración Schlumberger sigue atribuyéndola a una profundidad de 3 unidades. La Figura 5.39 muestra los valores propios de las matrices de sensibilidad relativas (al descomponerlas con la técnica SVD) para las dos configuraciones electródicas. La configuración doble dipolo está mejor condicionada, por lo que será menos vulnerable al ruido. Es decir, para un factor λ similar (al ser la relación S/N la misma para las dos configuraciones), el número de vectores propios no “atenuados” será mayor para la configuración doble dipolo y, por tanto, contendrá más información sobre las celdas más profundas (que como vimos estaban asociadas a los valores propios menores), lo que explica su mejor comportamiento. Problema inverso 5-48 La resolución del problema inverso consiste básicamente en invertir la matriz de sensibilidad para obtener la distribución de conductividades. Sin embargo el problema está mal definido, lo que implica que la matriz de sensibilidad está mal condicionada y, por tanto, pequeños errores en los datos provocan grandes variaciones en la solución. Los métodos de mínimos cuadrados tratan de resolver esta dificultad regularizando el problema, que consiste básicamente en añadir información a priori. La regularización de Tikhonov es uno de los métodos más conocidos. Los algoritmos de Marquardt-Levenberg y de Occam son casos particulares. El método TSVD deriva de la descomposición SVD de la matriz de sensibilidad. Tarantola y Valette (1982) abordan el problema inverso desde un punto de vista probabilístico. Su método es general y algunos de los algoritmos clásicos pueden ser obtenidos como casos particulares. Los métodos de retroproyección, que provienen en su mayoría de la tomografía médica, regularizan el problema utilizando la traspuesta de una matriz de sensibilidad ponderada. Los algoritmos de reconstrucción (de la conductividad) se comparan a partir de una estimación cualitativa y cuantitativa de las imágenes 2D obtenidas. La mayoría de autores obtienen los datos a partir de métodos numéricos (elementos finitos, diferencias finitas). En este trabajo los datos se obtienen a partir de la solución analítica para una esfera (Capítulo 3) conductora y aislante inmersa en el subsuelo. Se utilizan dos configuraciones electródicas: la doble dipolo y la Schlumberger. Los métodos de mínimos cuadrados son claramente superiores a los métodos de retroproyección. El estudio de la matriz de resolución de los parámetros explica estos resultados. Los métodos de Marquardt y TSVD obtienen los menores errores. El método de Occam suaviza la imagen alargando verticalmente el objeto por lo que su localización es más imprecisa. El parámetro λ se puede determinar automáticamente para el método de Marquardt, por lo que éste es el método escogido. Un valor de λ demasiado pequeño inestabiliza la solución y la imagen es muy “ruidosa”. Por el contrario, un valor de λ demasiado grande estabiliza demasiado la solución y el objeto se localiza incorrectamente en las capas superiores. La no-linealidad del problema original se resuelve normalmente mediante el uso de métodos iterativos. En este trabajo se tiene en cuenta aplicando un factor de corrección a las imágenes obtenidas. Si el modelo escogido (discretización del subsuelo) no es adecuado las imágenes pueden ser irreales, apareciendo, por ejemplo, conductividades negativas. Si el número de celdas es muy superior al número de medidas el problema es muy indeterminado, lo que hace difícil encontrar una solución adecuada. Esto se solventa, aparte de reduciendo el número de celdas, añadiendo información a priori (utilizando los métodos de retroproyección) o aumentando el número de medidas (reales o interpoladas). Los errores en las medidas afectan la localización de los objetos más profundos. La configuración doble dipolo obtiene mejores resultados que la Schlumberger cuando en aquélla se eliminan las medidas más erróneas. Además, al ser el error dominante proporcional a la medida, es mejor utilizar una matriz de sensibilidad relativa. En el caso de tener objetos cuya resistividad no varía a lo largo del eje y se realiza una integración previa a lo largo de dicho eje para el cálculo de los coeficientes de sensibilidad. Además, el factor de corrección de las imágenes es diferente. Si el objeto no se encuentra en el mismo plano vertical que contiene la Problema inverso 5-49 agrupación de electrodos la posición estimada del objeto no es correcta. La utilización de imágenes y medidas 3D permite una localización correcta.