MÉTODOS ESTADÍSTICOS EN PROBLEMAS ESPACIALES por
HECTOR NICOLAS FUNES Profesor en Matemática y Física
Tesis presentada en parcial cumplimiento para la obtención del grado de
MAGISTER EN ESTADISTICA APLICADA
Comité de Supervisión Dr. Raúl Pedro Mentz (Director) Dr. Aldo José Viollaz Mg. Santiago Mario Di Lullo
INSTITUTO DE INVESTIGACIONES ESTADISTICAS FACULTAD DE CIENCIAS ECONOMICAS UNIVERSIDAD NACIONAL DE TUCUMAN
San Miguel de Tucumán, Mayo 2004
Índice Índice ................................................................................................................................ i Agradecimientos ............................................................................................................ iv Resumen .......................................................................................................................... v CAPITULO 1: INTRODUCCIÓN ............................................................................... 1 1.1 Estadística para datos espaciales. ........................................................................... 1 1.1.1 Algunos ejemplos de datos espaciales............................................................... 1 1.2 Modelo espacial general. .......................................................................................... 2 1.2.1 Ejemplos de datos en Geoestadística. ............................................................... 3 1.2.2 Ejemplos de datos lattice. .................................................................................. 3 1.3 Geoestadística: definición y alcance........................................................................ 4 CAPITULO 2: ESTADÍSTICA DESCRIPTIVA PARA DATOS GEOESTADÍSTICOS.................................................................................................... 6 2.1 Primeras herramientas exploratorias..................................................................... 6 2.1.1 Presentación de los datos. Mapas de datos. ..................................................... 6 2.1.2 Box Plot y/o Diagrama de tallos y hojas........................................................... 7 2.1.3 Scatter plot tridimensional. ............................................................................... 8 2.1.4 Post plot............................................................................................................... 8 2.2 Otras herramientas exploratorias........................................................................... 9 2.2.1 Diagrama de puntos de Variable versus Variable índice. ............................ 10 2.2.2 Medias y medianas por filas y columnas....................................................... 10 2.2.3 El Estadístico u. .............................................................................................. 12 2.2.4 h- Scatter Plot .................................................................................................. 13 2.2.5 Nube de variogramas. ..................................................................................... 16 2.2.6 Nubes de las diferencias de las raíces cuadradas. ........................................ 18 2.2.7 El Pocket plot................................................................................................... 19 CAPITULO 3: CAMPOS ALEATORIOS................................................................. 23 3.1 Campos Aleatorios.................................................................................................. 23 3.2 Función de distribución y momentos de un campo aleatorio. ............................ 23 3.3 Campos aleatorios estacionarios. .......................................................................... 24 3.4 Isotropía................................................................................................................... 26 3.5 Meseta y Alcance. ................................................................................................... 27 3.6 Campos Aleatorios Intrínsecamente Estacionarios............................................. 27 3.7 Variograma vs. Covariograma.............................................................................. 29 3.7.1 Propiedades del Variograma........................................................................... 29 3.7.2 Covariograma y Correlograma....................................................................... 30 CAPÍTULO 4: ANÁLISIS ESTRUCTURAL............................................................ 32 4.1 Estimación del variograma. ................................................................................... 32 4.1.1 Método de los momentos. ................................................................................ 32 4.1.2 Variograma para datos irregularmente espaciados...................................... 34 4.1.3 Estimación robusta del variograma. .............................................................. 38 4.2 Modelos de variograma.......................................................................................... 39 4.2.1 Modelos de semivariogramas isotrópicos....................................................... 39 Construcción de otros modelos de variogramas en Rd. ...................................... 44 4.2.2 Anisotropía........................................................................................................ 44
i
4.3 Ajuste a modelos de variogramas. ........................................................................ 46 4.3.1 Máxima Verosimilitud (ML). .......................................................................... 47 4.3.2 Mínimos cuadrados. ......................................................................................... 47 Ajuste mediante mínimos cuadrados generalizados........................................... 48 4.3.3 Ajuste a sentimiento. ........................................................................................ 49 4.4 Validación cruzada del variograma ajustado. ..................................................... 49 CAPITULO 5: KRIGING. .......................................................................................... 51 5.1 Prediccion Espacial y Kriging. .............................................................................. 51 5.2 El mejor predictor lineal - Kriging simple. .......................................................... 53 5.3 Kriging ordinario.................................................................................................... 55 Obtención de las ecuaciones de kriging ordinario. ............................................. 56 5.4 Aspectos prácticos. ................................................................................................. 61 5.4.1 Efectos de distribución de los datos. ............................................................... 61 Efecto pantalla........................................................................................................ 63 Efecto de agrupamiento......................................................................................... 65 Distribución angular de los datos. ........................................................................ 65 5.4.2 Dibujo de curvas de nivel. ................................................................................................66 5.5 Kriging lognormal. ................................................................................................ 68 CAPÍTULO 6: KRIGING UNIVERSAL. .................................................................. 70 6.1 Kriging Universal. .................................................................................................. 70 6.1.1 Modelo supuesto. .............................................................................................. 70 6.1.2 Predictor Supuesto. .......................................................................................... 71 6.1.3 Predicción espacial óptima del proceso Z. ..................................................... 72 6.1.4 Ecuaciones de kriging universal. .................................................................... 72 6.2 Estimación del variograma para el kriging Universal. ....................................... 74 6.3 Kriging mediana polish.......................................................................................... 76 Datos grillados........................................................................................................ 77 Datos no grillados................................................................................................... 77 6.3.1 Mediana polish. ................................................................................................ 78 6.3.2 Superficie mediana polish................................................................................ 82 6.3.3 Kriging basado en los residuos de la mediana polish.................................... 82 CAPITULO 7: APLICACIÓN A LA HIDROGEOLOGÍA, "ACUIFERO DE TUCSON” ..................................................................................................................... 84 7.1 Análisis descriptivo de los datos. ........................................................................... 84 7.1.1 Otras herramientas exploratorias. ................................................................. 87 7.2 Análisis Estructural................................................................................................ 89 7.2.1 Estimación del variograma.............................................................................. 89 7.2.2 Ajuste a modelos de semivariogramas. .......................................................... 92 Modelando con el Geoeas y el Variowin. ............................................................. 93 7.2.3 Isotropía. ........................................................................................................... 99 7. 3 Kriging.................................................................................................................. 101 7.4 Logaritmo del contenido de calcio. ..................................................................... 103 7.4.1 Presentación y análisis descriptivo de los datos. ......................................... 103 7.4.2 Semivariograma para la variable logaritmo del contenido de Calcio. ...... 104 7.4.3 Kriging de la variable logaritmo del contenido de calcio. .......................... 106 7.5 Comparación de resultados. ................................................................................ 108
ii
CAPITULO 8: APLICACIÓN A LA HIDROGEOLOGIA: “DATOS DEL ACUÍFERO DE LA CALDERA” ............................................................................. 110 8.1 Introducción. ......................................................................................................... 110 8.2 Ubicación de los pozos.......................................................................................... 111 8.3 Concentración de cloruros................................................................................... 111 8.3.1 Análisis descriptivo. ....................................................................................... 111 8.3.2 Análisis estructural. ....................................................................................... 114 8.3.2.1 Estimación del variograma. ................................................................... 114 8.3.2.2 Ajuste a un modelo de semivariograma............................................... 115 8.3.3 Isotropía. ........................................................................................................ 116 8.3.4 Kriging............................................................................................................ 117 8.4 Nitratos. ................................................................................................................. 119 8.4.1 Análisis descriptivo. ....................................................................................... 119 8.4.2 Análisis estructural. ....................................................................................... 122 8.4.2.1 Estimación del variograma. .................................................................... 122 8.4.2.2 Ajuste a un modelo de semivariograma................................................ 122 8.5 Mapas generados por el programa Surfer. ........................................................ 125 8.6 Conclusiones.......................................................................................................... 128 CAPITULO 9: UNA APLICACIÓN A LAS CIENCIAS DEL MEDIO AMBIENTE ................................................................................................................ 129 9.1 Introducción. ......................................................................................................... 129 9.2 Presentación de los datos. .................................................................................... 129 9.3 Modelado de la tendencia. ................................................................................... 133 9.4 Kriging de los residuos. ........................................................................................ 134 9.4.1 Estimación del variograma............................................................................ 134 9.4.2 Ajuste a un modelo de semivariograma. ...................................................... 135 9.4.3 Kriging............................................................................................................. 138 9.4 Kriging Mediana Polish. ...................................................................................... 139 CAPÍTULO 10: CONCLUSIONES.......................................................................... 142 CAPITULO 11: BIBLIOGRAFIA............................................................................ 144 APÉNDICE C5 ........................................................................................................... 146 1. Obtención de las ecuaciones del Kriging simple................................................. 146 2. Obtención de las ecuaciones del kriging ordinario............................................. 147 APÉNDICE C6 ........................................................................................................... 152 1. Obtención de las ecuaciones del Kriging Universal en términos de variogramas. ................................................................................................................................... 152 2. Obtención de las ecuaciones del Kriging Universal en términos de covariogramas. ......................................................................................................... 155
iii
Agradecimientos Deseo expresar mi agradecimiento al Dr. Raúl Pedro Mentz, director de esta tesis, mi guía durante estos largos años de vida donde estuvo inserta la realización de este trabajo. Él junto al equipo de docentes, investigadores y personal de apoyo del INIE influyeron notablemente en mi formación académica y humana. Deseo agradecer a todos mis compañeros del Magister y en especial a mi amiga María Cristina Ahumada quien constantemente apoyó y alentó la empresa que emprendí con mucho esfuerzo. Mi gratitud para con mis profesores del profesorado, en especial a la Lic. Elda Canterle y al Estadístico Di Veltz, quienes despertaron en mi el interés por la Estadística en mis épocas de estudiante de grado en la Universidad Nacional de Salta. A mis compañeros de trabajo mi reconocimiento por sus aportes y apoyo. Al igual que a ambas universidades, la de Salta y Tucumán, por el apoyo económico que hicieron posible la realización de este programa de Magister. A los grandes forjadores de las ideas matemáticas y estadísticas por permitirme compartir y disfrutar su ciencia. millón de gracias
millón de gracias millón de gracias a mi familia Finalmente, les doy un por el sacrificio, espera, aguante, y el compartir un papá con la que día tras día crecía más y más: la tesis. Quiero dedicar de todo corazón esta tesis a toda mi familia, mis padres, mis hermanos. En especial a Marisa mi compañera de toda la vida, a mis retoños: Beatriz Analía, Nicolás Alberto, Héctor Raúl y Juan Pablo.
iv
Resumen En este trabajo se presentan once capítulos dedicados al tema "Métodos Estadísticos en Problemas Espaciales". De entre dichos métodos se centraliza en el enfoque para tratar datos espaciales que se denomina Geoestadística. En los primeros ocho capítulos se pone énfasis en los resúmenes de lecturas de la bibliografía, parte del desarrollo del trabajo de tesis plasmado en el plan de trabajo. En el primer capítulo se presenta un marco general para el análisis y modelado de datos espaciales. Mediante distintos ejemplos provenientes de distintas ciencias se muestra la necesidad de ese nuevo enfoque. Luego se define e indica los alcances del enfoque que se trata en este trabajo: el enfoque Geoestadístico. El segundo capítulo esta dedicado a la presentación de las herramientas descriptivas y exploratorias del enfoque Geoestadístico. A través de un ejemplo se muestra como ellas son utilizadas para tal fin. A los gráficos y resumen clásicos se deben agregar otros instrumentos que enfaticen la distribución de los valores de los datos y las posiciones de los mismos. Los datos espaciales deben ser pensados como una realización de una colección espacial de variables aleatorias dependientes, cuya dependencia está fuertemente ligada a las ubicaciones espaciales. Es necesario investigar en forma exploratoria los supuestos que se hagan sobre la distribución de los valores de los datos y las estructuras de dependencia. Para ello se presentan nuevas herramientas. En el tercer capítulo se presentan los conceptos teóricos que sustentan la Geoestadística. Se definen los conceptos de campos aleatorios, los momentos de primer y segundo orden de un campo aleatorio, destacando la herramienta fundamental de la Geoestadística: el variograma. Como así también las definiciones necesarias de campo aleatorio estacionario, intrínsecamente estacionario, y el carácter isótropico de los mismos. Se presentan las propiedades que debe tener el variograma como así también la función de covarianza estacionaria a los efectos de comparación. En el cuarto capítulo se trata la estimación del variograma. Se trabaja en presencia del supuesto de isotropía, y se expone la corrección en caso de anisotropía geométrica. Para realizar una predicción en una determinada posición es preciso conocer el valor del variograma empírico, que no siempre esta definido en dicha posición, y se necesita ajustar algún modelo de variograma teórico. Por eso en este capítulo se presentan modelos de variogramas teóricos y criterios de ajuste de los mismos a los variogramas empíricos. Por último se trata la validación cruzada que es una forma de medir el ajuste y diagnosticar algunos problemas con el mismo. En el quinto capítulo se presenta los fundamentos teóricos básicos de la predicción espacial que se usará en el desarrollo de este trabajo, poniendo énfasis en el método de kriging. Se fundamenta que el kriging es sinónimo de predicción óptima en algún sentido, y en base a los supuestos necesarios se deducen las ecuaciones del kriging simple y kriging ordinario. En una segunda parte se presentan características prácticas de esta metodología de predicción, poniendo énfasis en como contribuyen los
v
datos en la predicción de acuerdo a sus posiciones. Por último se introducen los resultados teóricos referentes al kriging lognormal. En el capítulo cinco se presentaron los métodos de predicción bajo el supuesto de estacionarierad intrínseca. En el sexto capítulo se presentará el método de predicción general cuando la variable de interés presente algún tipo de tendencia. Así se presentan los supuestos básicos de lo que se denomina "kriging universal”, tanto como las ecuaciones y las soluciones de las mismas. Este método fue el primero que planteó el problema de la predicción de funciones no intrínsecas de una forma global. También como un caso particular de este tipo de predicción se presentará el método denominado kriging mediana polish. En los capítulos siete a nueve se presentan aplicaciones, porque la Geoestadística está muy orientada a la aplicación. Así además de sus conceptos y técnicas, se requiere entender conceptualmente el fenómeno en estudio y trabajar interdisciplinariamente con los científicos "dueños de los datos". En el capítulo siete y en el siguiente se presentan problemas similares pero con diferentes fuentes que involucran datos cuyo tratamiento requiere de las herramientas que han sido presentadas a lo largo de los capítulos anteriores. Se trabaja con datos correspondientes a variables geoquímicas de aguas subterráneas. Estos lotes de datos espaciales provienen del campo de investigación de la Hidrogeología. Así, en primer lugar se trabaja con un conjunto de datos obtenidos a través de Internet de la página http:/www.u.arizona.edu/. De acuerdo a la información suministrada en dicho sitio se pretendió caracterizar el acuífero de Tucson. En el octavo capítulo se trabaja con datos del sistema acuífero de La Caldera. Para el desarrollo de la investigación, se contó con la información proporcionada por la cátedra de Hidrogeología en el marco del proyecto de investigación Hidrogeología del Sistema Acuífero La Caldera. Se trató la concentración de cloruros y la concentración de nitratos de las aguas subterráneas, por que de acuerdo a los especialistas son las variables químicas de importancia para la caracterización del sistema acuífero de La Caldera, a los efectos de determinar posibles focos de contaminación antrópica. Es de destacar que para llevar a cabo este trabajo se realizaron las distintas etapas de un estudio estadístico, así la primer tarea fue la elaboración de la base de datos. En ella principalmente se registró la información de la ubicación de los pozos donde se realizaron las observaciones de las variables de interés. En el noveno capítulo se aplicó las herramientas proporcionadas por el enfoque estadístico presentado en los 6 primeros capítulos a un problema proveniente de las denominadas ciencias del medio del ambiente. Para ello se utilizó información suministrada por el equipo de investigación del proyecto 577: “Determinación de Dióxido de Nitrógeno en la Atmósfera de Salta (capital)” de la Facultad de Ciencias de Exactas. Este grupo de investigación indagó acerca de la calidad del aire que se respira en la ciudad de Salta y por lo tanto se interesan en las concentraciones de los contaminantes mayoritarios en la troposfera en especial la concentración de dióxido de nitrógeno. Finalmente, en el capítulo décimo figuran las conclusiones y en el capítulo undécimo la bibliografía consultada.
vi
CAPITULO 1: INTRODUCCIÓN A través de este trabajo se pretende presentar un enfoque para el análisis y modelado de datos espaciales. Previo a ello en este capítulo se presenta un marco general para luego en su último apartado definir e indicar los alcances del enfoque propuesto.
1.1 Estadística para datos espaciales. Muchos de los modelos estadísticos simples consideran muestras aleatorias. Estas presuponen variables aleatorias independientes e idénticamente distribuidas. La independencia es un supuesto conveniente que hace que la teoría estadística sea más tratable. Sin embargo los modelos que involucran dependencia estadística son más reales. Así por ejemplo, los modelos de series de tiempo son basados en observaciones de una muestra de variables aleatorias idénticamente distribuidas que son dependientes y ocurren generalmente en tiempos igualmente espaciados. Los datos espaciales son otro ejemplo de la necesidad de crear modelos que involucren dependencia entre las variables. Existen muchas disciplinas que trabajan con datos recolectados desde diferentes ubicaciones en el espacio, estos son los datos espaciales. La noción de que ellos (al igual que los temporales) pueden estar muy juntos posibilita la correlación de los mismos, o sea no pueden ser modelados como estadísticamente independientes. A diferencia de los datos temporales, en los datos espaciales la dependencia está presente en todas las direcciones, y en general se vuelve más débil cuando las localizaciones de los datos están más alejadas. 1.1.1 Algunos ejemplos de datos espaciales. En geología por ejemplo, constituyen datos espaciales las medidas de contenido mineral en cada punto del espacio. En ciencias de la atmósfera: el registro mensual de lluvias en distintas localidades. Observemos en este caso que tienen tanto componentes espaciales como temporales. Así, para cada una de las localidades, los datos forman una serie de tiempo, y en un determinado tiempo, los datos de las distintas localidades son espaciales. En los últimos tiempos científicos recolectaron conjuntos de datos espacio temporales meteorológicos para estudiar los efectos de la polución atmosférica, en particular el fenómeno denominado lluvia ácida. En ecología, el tamaño, la forma de los organismos, y la densidad de la distribución local, deben tener en cuenta la ubicación de los datos. Así, estas y otras áreas del conocimiento como procesamiento de imágenes, geología, epidemiología, agronomía, silvicultura, astronomía, etc., necesitan desarrollar modelos que tengan en cuenta la existencia de dependencia entre las medidas en diferentes ubicaciones, es decir se precisan modelos para trabajar con los datos espaciales. Estos modelos deben ser más flexibles que sus contrapartes temporales, porque el pasado, presente, y el futuro no tienen analogía en el espacio. Las primeras manifestaciones de la estadística sobre datos espaciales aparecen en la forma de mapas de datos. Por ejemplo Halley (1686) superpuso sobre un mapa, las direcciones de cambio de los vientos alisos y monzones entre y en las cercanías de los trópicos, y procuró asignarles causas físicas. Mucho tiempo después aparecen los modelos espaciales. Así en 1907, Student estudió la distribución de partículas a través de los líquidos.
1
Fisher conocía claramente la dependencia en experimentos agrícolas. Entre los años 1920 y 1930, en la estación experimental Rothamsted en Inglaterra, estableció los principios de aleatorización, bloqueo y replicación a los efectos de suprimirla. Además de controlar el sesgo no deseado, la aleatorización también neutraliza el efecto de correlación espacial, aunque no la remueve (Yates, 1938). Sin embargo, se debe tener en cuenta que la aleatorización no neutraliza la correlación espacial en escalas espaciales más grandes o más pequeñas que las dimensiones de las parcelas. Algunas veces se puede esperar que las observaciones cercanas tiendan a ser disímiles. La competencia entre plantas por la luz y los nutrientes del suelo podría inducir la existencia de plantas más sanas rodeadas por plantas menos robustas. Mead (1967) analizó los datos sobre coliflores e investigó el efecto de competición. Los métodos de los vecinos más próximos para el análisis de los experimentos agrícolas intentan tener en cuenta la dependencia espacial, indirectamente, por el uso de los residuos desde las parcelas vecinas como covariantes o por diferenciación. En la actualidad, la ciencia Estadística puede encontrarse en cualquier disciplina científica cuantitativa, de esta manera se producen interesantes variaciones en los temas que fueron tradicionales. En áreas tales como geología, ecología y ciencias del medio ambiente, no es a menudo posible o no es apropiado aleatorizar, bloquear y replicar los datos. Existe una necesidad de nuevos modelos y aproximaciones estadísticas que se dirijan a nuevos problemas surgidos desde la tecnología. Muchos de los problemas resultantes, tales como valuación de recursos, monitoreo ambiental (por ej. el calentamiento global), e imágenes en medicina, son espaciales en su propia naturaleza.
1.2 Modelo espacial general. Supongamos que se mide el grado de mineralización en una mina de cobre. Se puede pensar que las mediciones del contenido de cobre son realizaciones de una variable aleatoria Z(s) (Z ∈ R1), en cada punto s del espacio considerado: un segmento de recta, un subconjunto del plano ó un subconjunto del espacio tridimensional. Muchas veces, son varias las propiedades o características observadas para cada s. Por ejemplo si se considera el grado de mineralización en una mina de cobre, el cobre puede ser de principal interés pero otros metales como el molibdeno, la plata, el oro y el zinc pueden también ser de interés crítico. Así para cada s, las mediciones se pueden pensar como realizaciones de un vector aleatorio Z(s). Además s varía en un conjunto D ⊂ Rd (en este caso d=1, 2 ó 3) a fin de generar el campo aleatorio multivariado: Z(s) : s ∈ D ⊂ R d (1.2.1) Esto permite que los datos espaciales puedan ser pensados como los resultados de un campo aleatorio multivariado {Z(s): s ∈ D⊂ Rd} donde Z ∈ Rm, m ≥ 1, y la ubicación s puede ser un punto de un espacio Euclideano d- dimensional con d ≥1. Una realización de (1.2.1) se denota por {z(s): s ∈ D}. Cressie(1991) consideró que la mayoría de los problemas de la estadística espacial están identificados en una de cuatro categorías según D sea fijo ó aleatorio y D sea continuo ó discreto. En este trabajo de tesis nos interesa trabajar con problemas de estadística espacial que se pueden pensar como una realización parcial del campo aleatorio (1.2.1) con D fijo y continuo. Presentaremos estos problemas y los correspondientes a D fijo
{
}
2
discreto, a los efectos de comparar. Además consideraremos el caso más simple, el caso univariado. Si los datos se pueden suponer provenientes de una realización parcial de un campo aleatorio donde D se considere fijo y continuo se dice que los datos son geoestadísticos. En cambio si se pueden suponer provenientes de una realización parcial de un campo aleatorio donde D se considere fijo y discreto se dice que son datos lattice. Un lattice de ubicaciones evoca una idea de puntos espaciados regularmente en Rd, pero este será denominado lattice regular, dando la posibilidad de que sea irregular. 1.2.1 Ejemplos de datos en Geoestadística. Además del grado de mineralización en un bloque de mineral, que motivó a Matheron pionero de la Geoestadística, para crear el método de predicción kriging que estudiaremos, se pueden citar otros ejemplos de datos geoestadísticos. El pH del suelo en agua, la conductividad eléctrica en el suelo, la tensión sueloagua, la infiltración agua – suelo y la erosión del suelo debido al agua, son algunas de las variables que son de interés para los científicos de la tierra, quienes realizan mapas de las propiedades del suelo de un campo a partir de un número pequeño de observaciones en ubicaciones conocidas, suponiendo que varían continuamente a lo largo del campo. La distribución espacial de la leucemia, y los niveles de ruidos acústicos, polvos y químicos en un ambiente de trabajo, son algunos de los problemas abordados por los científicos de la Salud. La cartografía, fue uno de los primeros dominios de aplicación de los métodos propuestos por Matheron para modelar numéricamente las características de los terrenos. El precio de mercado de los cerdos, o de las frutas, son variables que generan datos del tipo geoestadísticos. En hidrogeología, las variables que están relacionadas con el flujo (nivel piezométrico y trasmisividades) y la calidad química e isotópica de las aguas subterráneas, son fuentes de datos de interés. 1.2.2 Ejemplos de datos lattice. Los sensores remotos de satélites ofrecen un medio eficiente de recolección de datos. Por ejemplo estos permiten que se reúna información rápidamente sobre los patrones del estado atmosférico, la distribución mineral y las extensiones de terreno sembrado sin tener una labor intensiva de reconocimiento tradicional. Los satélites orbitan la tierra y reciben los datos en forma de ondas electromagnéticas reflectantes con una cantidad de frecuencias, incluyendo aquellas de la parte visible del espectro. Para diversos muestreos y métodos de integración, la tierra está dividida en pequeños rectángulos llamados pixeles. Supóngase una zona agrícola de interés (alrededor de 34000 km2) que tiene ciertas proporciones dedicadas al trigo, maíz, soja, y así sucesivamente, que necesitan ser estimadas. Estos cultivos tienen sus propias propiedades de reflectancia que, junto con ruidos, son remotamente sensorizados. Así, los datos son recibidos como un lattice regular en R2 (dejando de lado la curvatura de la tierra) y son identificados con los centros de sus pixeles respectivos. Existe una gran superposición entre técnicas de sensores remotos y técnicas de imágenes en medicina, aunque las escalas espaciales son diferentes, la forma de los datos y las cuestiones de existencia preguntadas son a menudo semejantes. Los modelos
3
estadísticos para tales datos necesitan expresar el hecho que observaciones cercanas tienden a ser semejantes. En contraste a los problemas geoestadísticos, los datos desde problemas de lattices pueden ser exhautivos del fenómeno. Por supuesto, es posible realizar muestreo; por ejemplo utilizando sólo una ventana pequeña de todos los datos. Los problemas geoestadísticos se distinguen claramente de los problemas de datos lattice por la capacidad del índice espacial s de variar continuamente sobre un subconjunto de Rd. Esto no afirma que los métodos de una clase de problemas no puedan ser apropiados para ser aplicados a la otra clase de problemas.
1.3 Geoestadística: definición y alcance. Previamente se verá qué se entiende etimológicamente por Geoestadística y algunas acepciones que fueron dadas por distintos estadísticos. El prefijo “geo” indica la relación estrecha entre el tema en estudio y todo lo perteneciente a la tierra, es decir en este caso la estadística relativa a la tierra, y verdaderamente este fue su significado original. Hart en 1954 le dió un contexto geográfico para denotar las técnicas estadísticas que enfatizan la localización dentro de las distribuciones zonales. Matheron usó el término en un contexto geológico para denotar la teoría y los métodos para inferir las reservas de mineral en bruto desde datos espacialmente distribuidos en un bloque (volumen de tierra y roca que puede ser minado si es suficientemente rico). En Francia bajo el impulso de Matheron principalmente, emergió antes de 1980 la disciplina Geoestadística como una mezcla de ingeniería en minas, geología, matemática y estadística, que a diferencia de otros enfoques más clásicos, tiene en cuenta la tendencia espacial y la correlación espacial que en terminología minera corresponden a variabilidad espacial en gran escala y en pequeña escala respectivamente. Las mismas ideas también fueron desarrolladas en forma independiente por Gandin en la Unión Soviética, pero aplicadas a la Meteorología. Sea {Z(s): s ∈ D} una función aleatoria (proceso aleatorio), donde D es un subconjunto fijo de Rd llamado conjunto de índices, y s varía continuamente en D. Matheron, padre de la Geoestadística en su forma actual, la definió como “la aplicación del formalismo de las funciones aleatorias al reconocimiento y estimación de fenómenos naturales”. Cressie(1991) considera a la Geoestadística como una rama de la estadística que abarca teorías y aplicaciones para procesos aleatorios con índices espaciales continuos. Así en ambos casos, el término Geoestadística pierde su significado etimológico. La definición dada por Cressie es más general que la de Matheron porque deja abierta la posibilidad de que el fenómeno en estudio sea proveniente de la naturaleza o no, y así responder a nuevos problemas que surjan de la tecnología. Para este trabajo, el término Geoestadística tendrá la última acepción. El fin último de la Geoestadística es la caracterización del fenómeno, lo que conduce a varios tipos de aplicaciones. El primero es la predicción (estimación), esto es, conocida la variable en una serie de puntos, predecir su valor en otros. La innovación de la Geoestadística es que permite obtener no sólo la predicción sino también una medida de la incertidumbre asociada a ella. La predicción suele producir mapas que son mucho más “suaves” que la realidad. Por ello, en los casos en que la variabilidad espacial sea de interés es necesario recurrir a técnicas de simulación, segundo grupo de aplicaciones, a fin de obtener
4
realizaciones plausibles de la variable estudiada. Otros tipos de aplicaciones son las que resultan del hecho de que al proporcionar medidas sobre la incertidumbre de la predicción, la Geoestadística constituye el marco ideal para seleccionar la ubicación de puntos de muestreo de forma que se minimice la incertidumbre de la predicción. Se puede advertir que los objetivos del Análisis de Series Temporales son muy similares a los de la Geoestadística. De hecho, las diferencias no son tanto de tipo conceptual como de campo de aplicación. Sin embargo, como el desarrollo de ambos ha sido independiente, sus nomenclaturas son bastantes diferentes, de forma que el intercambio de técnicas no es inmediato pese a no revestir dificultades fundamentales. Por otro lado, los tipos de problemas que se plantean son lo suficientemente distintos como para que dicho intercambio no sea todo lo provechoso que cabría esperar. Así, la metodología Geoestadística está pensada para datos distribuidos de forma arbitraria en el espacio, por lo que sus técnicas son más generales pero menos potentes que las del Análisis de Series Temporales. Como se deduce de la definición, la Geoestadística está muy orientada a la aplicación por lo que se requiere no sólo el conocimiento de las técnicas y metodología de la misma, sino también entender conceptualmente el fenómeno en estudio. Así a lo largo de los próximos cinco capítulos se presentará el método Geoestadístico y en los últimos se aplicarán a distintos campos del conocimiento científico.
5
CAPITULO 2: ESTADÍSTICA DESCRIPTIVA PARA DATOS GEOESTADÍSTICOS. Antes de dar una presentación más formal de los modelos geoestadísticos, se presentará un ejemplo y se analizarán sus datos en una forma exploratoria a los efectos de mostrar las distintas herramientas utilizadas para tal fin. Los datos espaciales son a menudo observados y resumidos sólo en las formas clásicas, como si fueran resultados de una muestra aleatoria o quizás de una serie temporal, en vez de ser pensados como una realización de una colección espacial de variables aleatorias dependientes, cuya dependencia está fuertemente ligada a las ubicaciones espaciales. Veremos cómo podemos hacer una exploración descriptiva de los datos.
2.1 Primeras herramientas exploratorias. 2.1.1 Presentación de los datos. Mapas de datos. Los siguientes datos se obtuvieron de Internet y representan el análisis de 60 observaciones de suelos en una región contaminado con arsénico, cadmio y plomo. La Figura 2.1 muestra las ubicaciones donde fueron tomados. La distribución de las coordenadas define aproximadamente un rectángulo de 240 pies en la dirección oeste este por 210 pies en la dirección sur – norte. Las ubicaciones de las observaciones están irregularmente espaciadas, aunque existen algunos claros y algunas agrupaciones. El mapa de las ubicaciones de las mediciones sirve además para identificar posibles errores en las coordenadas.
Figura 2.1: Posiciones de los datos correspondientes al estudio de la contaminación del suelo.
Para ejemplificar en las siguientes secciones, analizaremos las mediciones de la variable “contenido de cadmio”. Las medidas de la misma en p.p.m. se presentan con sus respectivas posiciones en la Figura 2.2.
6
Norte
Este Figura 2.2: Medidas de cadmio en sus respectivas posiciones
2.1.2 Box Plot y/o Diagrama de tallos y hojas. Por su propia naturaleza, el análisis exploratorio de datos no considera las observaciones de la misma manera, ya que aísla para una lectura cuidadosa las observaciones atípicas. La existencia de algún modelo subyacente, será investigado por el análisis exploratorio. En el contexto espacial, un modelo estocástico fundamental es aquel en el que todos los datos provienen de una distribución normal conjunta, cuya estructura de correlación depende de las ubicaciones espaciales. Claro está, que esta hipótesis es muy difícil de investigar, pero se puede indagar si el conjunto de datos independientemente de sus posiciones se puede suponer proveniente de una distribución normal. Al tratar a los datos como un conjunto de números ignorando su condición espacial, para investigar la existencia de valores alejados podríamos realizar un diagrama de tallos y hojas y/o un box-plot. Además dichos gráficos nos brindan alguna idea de la forma de la distribución.
Figura 2.3:Box-Plot para el contenido de Cadmio.
De la Figura 2.3 se observa que la distribución es prácticamente simétrica, y sin la existencia de valores alejados. Claro está que un histograma sería la herramienta exploratoria clásica más usual. Ambas herramientas, se presentan simultáneamente al realizar una exploración de los datos con el programa GEOEAS.
7
2.1.3 Scatter plot tridimensional. Los mismos datos se pueden representar en un scatter plot tridimensional, como el de la Figura 2.4. Un scatter plot tridimensional es una representación gráfica donde en los ejes x e y se ubican las coordenadas de los puntos de ubicación (Norte - Este) y en el tercer eje se ubican los valores de la variable de interés (Cadmio). La utilidad del mismo como una herramienta exploratoria es bastante limitada, porque los rasgos interesantes en los datos se pueden ocultar fácilmente, por ejemplo la tendencia que puede existir en una dirección particular. En el scatter plot se espera poder detectar valores alejados y observar posibles tendencias en direcciones determinadas. 250 300
Este
350 400 450
15 10 5 0 150 200 Norte
250 300
Figura 2.4: Scatter plot para las medidas de la variable cadmio.
Para la variable contenido de cadmio el scatter plot correspondiente muestra valores distribuidos en forma casi uniforme. El scatter plot es una herramienta útil cuando se lo programa en forma iterativa y se pueden rotar los ejes para visualizar posibles tendencias, valores muy alejados, etc. Esto, se puede realizar por ejemplo con el soft Mathematica. Para el conjunto de datos que se está tratando, los distintos gráficos generados por el programa, no muestran más que el de la Figura 2.4. 2.1.4 Post plot. El análisis exploratorio de datos espaciales requiere algunas nuevas propuestas, sugeridas por las técnicas de análisis univariadas y multivariadas. Uno de los más usadas es el gráfico Post- plot, que es similar al gráfico de la Figura 2.2, pero que en vez de ubicar los valores de la variable en su posición coloca símbolos y/o colores de acuerdo a una codificación que se aplica a los valores de la variable. Para codificar la variable se procede a dividir su recorrido de acuerdo a los cuartiles. En la figura 2.5 se representa el post plot para la cantidad de cadmio en cada observación. Con el triángulo y el color celeste se representan los valores correspondientes a los valores menores o iguales al primer cuartil, o sea a 5.3. Con la estrella y el color verde se representan los valores comprendidos entre 5.3 y el valor mediano 7.95. Con el diamante y el color azul se representan los valores comprendidos entre 7.95 y el tercer cuartil 10.8 y con cuadrados rojos a los valores mayores que el tercer cuartil.
8
Figura 2.5: Post Plot para la variable cadmio.
A partir de este gráfico se observa una tendencia general, los valores más altos ocurren en la banda central de la dirección oeste - este. Mientras los valores menores se presentan en bandas paralelas en el sur y en el norte. Además se detecta un posible valor alejado espacial: el valor 11.5 correspondiente a las coordenadas (290, 310), porque en las posiciones vecinas tienen valores pequeños. Este valor debería ser investigado. Otros gráficos similares pueden ser considerados al realizar la codificación utilizando mas percentiles, esto es cuando la cantidad de datos es muy grande para que tenga sentido la codificación. A estos últimos gráficos algunos autores los denominan Percentage Plots.
2.2 Otras herramientas exploratorias. Para tratar con datos geoestadísticos es necesario agregar supuestos al modelo estadístico subyacente. Existen distintas herramientas exploratorias para validar esos supuestos. Los métodos exploratorios necesitan ser resistentes a observaciones atípicas para un modelo subyacente. En un contexto no espacial, los valores alejados significan una desviación del modelo Gaussiano. En un contexto espacial, además de estas desviaciones pueden esperarse desviaciones de los supuestos de estacionariedad de la media y/o de la medida de dependencia, informalmente hablando esta estacionariedad se refiere a que las propiedades variacionales no cambian a través de la región de interés. Supongamos que el contenido de Cadmio en un punto s es la realización de un proceso aleatorio {Z (s) : s ∈ D} y está observado en los puntos {s i : i = 1, 2, K , n} de la región de interés. Entonces, un proceso intrínsecamente estacionario se define como aquel que cumple con los siguientes supuestos sobre las diferencias:
E (Z (s + h) − Z (s) ) = 0 Var (Z (s + h) − Z (s) ) = 2 γ (h)
(2.2.1) (2.2.2)
La cantidad 2γ (h ) es conocida como el variograma y γ (h) se denomina semivariograma. Estos conceptos serán ampliados en el próximo capítulo. 9
Las herramientas exploratorias para investigar si la información disponible a través de los datos no contradice al supuesto sobre la esperanza, o sea sobre (2.2.1), son entre otras: • Diagrama de puntos: Variable versus Variable índice. • Medias y medianas por filas y columnas. Esta última, sirve además para la detección de los posibles valores alejados, a través de un estadístico denominado u. En cuanto a las herramientas para investigar acerca de la influencia de los posibles valores extremos espaciales sobre el variograma, se utiliza: • El h-scatter-plot. • Nubes de variogramas en una dirección dada. • Nubes de las raíces cuadradas de las diferencias en una dirección dada. Estos dos últimos instrumentos se basan en los valores de las estimaciones del variograma. • Pocket plot Es una técnica para identificar un área localizada que sea atípica con respecto al modelo estacionario impuesto a la medida de dependencia, y también se basa en los valores de las estimaciones del variograma. Todas estas herramientas se detallarán en los siguientes apartados. 2.2.1 Diagrama de puntos de Variable versus Variable índice. Una herramienta exploratoria muy empleada es el clásico diagrama de puntos de la variable de interés versus cada una de las variables que definen la posición de los puntos. A los efectos de determinar la posible existencia de dependencia entre ellas. Por ejemplo en la Figura 2.6, en nuestro caso se trata de observar si a medida que aumenta la distancia oeste- este la cantidad de cadmio aumenta o disminuye o no se presenta ninguna tendencia.
a) Figura 2.6: Diagrama de puntos de: a) cadmio vs. Coordenada este.
b) b) cadmio vs. Coordenada norte.
En ambos diagramas se observa que no existe tendencia lineal entre las coordenadas y la variable cadmio. En la 2.6 b) se observa una leve tendencia cuadrática de la variable cadmio versus la variable que mide las coordenadas en la dirección norte. 2.2.2 Medias y medianas por filas y columnas. Si los datos están regularmente ubicados en una grilla (Figura 2.7), puedo identificar filas y columnas (direcciones O-E y S-N) y calcular medias y medianas muestrales a lo largo de ellas.
10
fila
columna Figura 2.7: Grilla de ubicación de datos que definen filas y columnas.
Como este conjunto de datos no está regularmente espaciado, es decir los datos no están ubicados en una grilla regular, se realiza un agrupamiento de baja resolución de las observaciones como en una tabla de doble entrada, como lo indica la Figura 2.8. Las filas (azul) y las columnas (negro) indican las posiciones donde se fijan los valores cercanos.
Figura 2.8: Grilla de baja resolución para realizar el cálculo de las medias muestrales.
Otra posibilidad es la de considerar a franjas en la dirección Este – Oeste (Sur – Norte) como filas (columnas), para el agrupamiento (Figura 2.9).
Figura 2.9: Definición de columnas por franjas en la dirección Sur- norte.
La Figura 2.10 es un intento de resumir la posible no estacionariedad en la media a lo largo de filas y columnas, es decir en las direcciones oeste- este y norte- sur, usando la media muestral y la mediana muestral a lo largo de las filas y las columnas respectivamente. Entonces, el aspecto espacial de los datos participa por primera vez. 11
a)
b)
Figura 2.10: En ambas figuras con azul se representan los valores de la media y con rojo los valores de la mediana de los valores de contenido de cadmio. a) Según columnas. b) Según filas.
En la Figura 2.10 a) se observa que no existe ninguna tendencia en los valores medios y medianos a través de las columnas. En cambio en la 2.10 b) se observa una tendencia cuadrática de los valores medios y medianos a través de las filas. Esta tendencia es más notable con las medias y medianas que con los datos (ver Figura 2.6 b)). 2.2.3 El Estadístico u. Otro propósito del cálculo de los valores medios y medianos por filas y columnas es el siguiente: desde el punto de vista del Análisis Exploratorio de Datos, la mediana es un resumen estadístico resistente, la comparación de la mediana con la media (resumen estadístico no resistente) tiene la función adicional de destacar las filas y las columnas que posiblemente presenten outliers espaciales. Si la diferencia media − mediana es bastante grande, entonces la fila o columna debe ser investigada por los posibles outliers. ¿Pero qué significa bastante grande? Veamos los siguientes resultados teóricos, extraídos de Cressie(1991), para determinar la magnitud admisible de dichas diferencias. Supongamos que Y1 , Y2 , K , Yn es una muestra aleatoria (es decir i.i.d.) de una población con función de densidad simétrica f, con valor esperado µ y varianza σ 2 . La mediana puede ser aproximada por: n signo(Yi − µ ) ~ 1 Y ≅µ+ (con f ( µ ) ≠ 0 ) 2 f (µ ) n i =1 y en forma trivial la media puede ser expresada como: n 1 (Yi − µ ) Y =µ + n i =1 Por lo tanto la diferencia: n signo(ε i ) ~ 1 ε i − donde ε i = Yi − µ i = 1, 2, K , n Y −Y ≅ 2 f ( µ ) n i =1
∑
∑
∑
Se supone que ε i = Yi − µ i = 1, 2, K , n se distribuyen normalmente. El siguiente
~ σ2 Var (Y − Y ) = n Así, la diferencia media − mediana estandarizada es: resultado es válido:
12
π − 1 2
(1)
u ≡n
1
2
~ Y −Y 0.7555 σˆ
Según Cressie se debe prestar atención a los valores de u que estén cercanos a 3 o sean más grandes que 3. En el contexto del análisis exploratorio de datos se usa como un estimador de σ a: rango intercuartil de las Y σˆ = 2 × 0.6745 La tabla 2.1 presenta las diferencias media − mediana estandarizadas, para los valores que estamos tratando. De la tabla 2.1 se observa que el valor u señala que en la columna 1 podría existir un outliers, así el valor sospechoso es z(254.4, 216.0)= 14.9 p.p.m.(ver Figura 2.5) que es muy grande comparados con los demás valores de la columna que no superan el valor del primer cuartil 5.3 p.p.m.. Cuando existe covarianza positiva entre las Y, (1) es una subestimación de la ~ verdadera varianza de la diferencia Y − Y y por lo tanto u indica la existencia de más filas y columnas atípicas que las que debería. Fila Coordenada Norte media u Fila
1 118.6
2 138
3 150
4 163
5 172
6 181.3
7 193
8 204.3
9 215.5
10 227
1.03 11
0 12
0 13
0.71 14
0 15
0.76 16
-0.18 17
0.59
0.27
0
Coordenada Norte media u
246
253.0
268.5
271.5
285
295
313
0.75
0.63
------
0
0.04
0.42
- 0.07
Columna 1 2 3 4 5 6 7 8 9 10 11 Coordenada 254.4 275.4 284.8 333.6 345.6 354.2 364.8 412.8 434.8 443.6 492.0 Este media u 3.35 0.04 -0.35 -0.87 0 -1.64 1.08 -0.7 2.02 -1.08 -1.21 Tabla 2.1: Diferencias entre las medias y las medianas estandarizadas según filas y según columnas.
2.2.4 h- Scatter Plot El h-scatter plot es un diagrama de puntos que representa los valores de la variable en la posición s versus el valor de la misma variable en la posición s + h, donde h es un vector separación. Si los datos están regularmente espaciados, los h-scatter plots son construidos de acuerdo al procedimiento ilustrado en la Figura 2.11.
Figura 2. 11: Construcción de h-scatter plot con datos irregularmente espaciados.
13
Si los datos están irregularmente espaciados, todos los pares que tengan una distancia de separación cercana a la del vector separación h serán retenidos para construir un h-scatter plot. En la práctica, se establece una tolerancia sobre el vector separación h. Esta tolerancia se puede establecer de diversas maneras, las dos más usadas son las siguientes: • Con centro en el punto final del vector separación h se considera un rectángulo de lados como se indica en la Figura 2.12 a). Así por ejemplo el vector h1 es retenido para la construcción del h-scatter plot en la dirección h. •
Se establece tolerancias sobre la distancia o sea sobre la magnitud del vector separación y también para su dirección, como se lo indica en la Figura 2.12 b). El rectángulo ahora es curvilíneo.
a)
b)
Figura 2.12: a) La tolerancia es aplicada en el punto final del vector h. b) La tolerancia es aplicada tanto en la dirección como en la magnitud del vector h.
El variograma es una medida de dependencia en la dirección h y Matheron propone como estimador al promedio de las diferencias cuadradas de los valores de la variable de interés, que se encuentran ubicadas en puntos que estén separados una distancia h = h y en la dirección de h en el caso más simple, esto y la situación en que los puntos estén irregularmente espaciados será tratado ampliamente en el capítulo 4. Fijada una dirección y h, el h-scatter plot nos muestra los pares de valores que están involucrados en el cálculo del semivariograma en la dirección fijada. Los puntos de un h-scatter plot que estén cerca de la recta y = x no influyen en gran medida en el valor del semivariograma. En cambio los puntos que se encuentran muy alejados de la recta son correspondientes a pares de valores cuyos aportes a la estimación del semivariograma serán grandes. Así el h-scatter plot tiene como principal función identificar pares que tienen una fuerte influencia en la estimación de la medida de la dependencia del proceso. Cuando el h-scatter plot muestre una forma de “mariposa volando”, no puede resumirse adecuadamente con un único número, es posible que una transformación de los datos, tales como la transformación logarítmica, la transformación de score normal, o la transformación de rangos, puedan corregir dicha situación. Cuando la nube de puntos del h-scatter plot es elíptica o casi elíptica; el conjunto de pares puede ser adecuadamente resumido con un solo estadístico, por ejemplo el coeficiente de correlación, la covarianza o el momento de inercia alrededor de la línea diagonal. Este último es precisamente el estimador propuesto por Matheron. Además, los h-scatter plots son frecuentemente útiles para chequear diversos supuestos que se realizan en geoestadística: • Los h-scatter plots de una variable univariada Gaussiana proveen una manera visual rápida de verificación para la hipótesis de multi-Gaussianidad. Si la función
14
•
aleatoria que representa el fenómeno es multi-Gaussiana, todos los h-scatter plots son nubes elípticas alrededor de la línea diagonal con una alta densidad de puntos en el centro de la nube. Una sucesión de h-scatter plots, calculado para incrementos de la magnitud del vector de separación h, provee una herramienta para investigar la estacionariedad. Si tal sucesión de h-scatter plots muestra que el centro de la nube de pares se desvía de la línea diagonal, entonces la función aleatoria que representa el fenómeno estudiado no puede ser considerada estacionaria.
Las figuras 2.13 y 2.14 muestran distintos h-scatter plots para la variable cadmio, para la dirección Oeste-Este, para distintos magnitudes promedios del vector de separación. Los h-scatter plot presentados a modo de ejemplo no muestran la existencia de posibles extremos. a)
b)
Figura 2.13: H-scatter plots para la dirección Oeste-Este con una tolerancia de 150. a)h=33.6 b)h=62.4. c)
d)
Figura 2.14: H-scatter plots para la dirección Oeste-este con una tolerancia de 150.c)h=91.03 d)h=125.1
Se puede ir obteniendo un conjunto de observaciones cada vez más grande, aunque una explicación de su “falta de comportamiento” puede ser la falta de estacionariedad en la media del proceso Gaussiano subyacente. Si ellos necesitan
15
especial atención será determinado a partir de su influencia en la estimación del variograma y en el método de predicción espacial: el kriging. 2.2.5 Nube de variogramas. El estimador para el variograma propuesto por Matheron es: 1 (Z (s i ) − Z (s j ) )2 2 γˆ (h) = N (h) N (h )
∑
donde la suma se realiza sobre todos los pares que contribuyen al cálculo del variograma empírico en la dirección del vector h. Por ser un promedio de diferencias elevadas al cuadrado, está afectado por observaciones atípicas. Una manera de investigar la influencia de posibles valores alejados sobre el variograma es a través de la herramienta exploratoria denominada nube de variogramas. Este instrumento consiste en graficar para una dirección dada, por ejemplo Oeste- Este, y para cada distancia de separación h = h los valores de las diferencias al cuadrado de los valores de la variable que contribuyen al cálculo del variograma o del semivariograma empírico. En forma más simple la nube de variogramas en la dirección e es simplemente un gráfico x - y de puntos de los valores (Z (s j ) − Z (s i ) )2 :s i + h e = s j , (i, j )∈ N (he) versus los valores de h. La Figura 2.15 muestra la nube de variogramas para la variable que estamos analizando en la dirección Oeste- Este, para 0 < h ≤ 150 pies .
{
}
Figura 2.15: Nube de semivariogramas en la dirección Oeste- Este.
Se debe adoptar un criterio para determinar que tan grande debe ser el valor que precise atención. Es por ello que se recurre al box plot de la nube de variogramas para cada h si es que los datos están distribuidos en una grilla regular o para h promedios correspondientes a grupos de h cuando los datos están irregularmente espaciados. La Figura 2.16 muestra los “box- plot” de la nube de variogramas en la dirección Oeste-Este, donde los h están agrupados en intervalos de amplitud 10. La Figura 2.17 muestra el mapa en donde se identifican los pares cuyos valores de diferencias de cuadrados son “muy grandes”, donde además se muestra que para los cálculos en la dirección en cuestión se tuvo en cuenta una tolerancia angular de 150.
16
Figura 2.16: Box Plot de la nube de variogramas correspondientes a la figura 2.15. Las unidades en el eje vertical son (p. p. m. de cadmio)2.Por ejemplo: el valor h =2 corresponde al intervalo 20 ≤ h < 30 .
Figura 2.17: Mapa en donde se indican los pares de posiciones cuyos cuadrados de las diferencias de contenido de cadmio son valores grandes con respectos a los demás valores para cada h.
La tabla 2.2 muestra los valores de la mitad del cuadrado de la diferencia de los contenidos de cadmio que son valores alejados. Una gran cantidad de pares de datos se debería considerar para un análisis exhaustivo. Pero, es importante destacar que la nube de variogramas puede llevar a confundir la asimetría de la distribución de los valores de γˆ (h) para cada h con los valores alejados debido a que las diferencias están elevadas al cuadrado. Esta última afirmación está sustentada en el siguiente resultado extraído de Cressie(1991): Si Z es un proceso Gaussiano, (Z(s + h) − Z(s))2 se distribuye como 2γ(h)χ21 donde χ21 es una variable aleatoria chi cuadrado con 1 grado de libertad; por lo tanto, 2γ(h) es el primer momento de una variable aleatoria altamente asimétrica.
17
h
Identificación
3 4 5 5 6 6 6 8 8 8 8 10 10 10 10 12 13
32 1 57 28 5 16 14 48 48 48 16 37 33 53 50 14 52
9 34 59 45 53 56 24 57 46 27 57 29 12 57 48 43 47
Distancia
Z(+)
39.44 46.68 52.80 51.382 61.846 63.739 64.830 81.937 84.148 85.00 89.746 101.052 101.117 105.60 107.6678 124.804 132.00
1.20 11.50 16.70 12.10 11.20 8.30 11.5 6.80 680 6.8 8.30 8.7 1.7 6.5 14.90 11.50 11.60
[Z ( +) − Z ( −)]2
Z(-) 5.2 1.2 6.90 1.60 6.50 11.0 5.3 16.7 15.0 14.5 16.70 0.9 9.5 16.7 6.8 4.40 3.40
2
8 53.045 48.020 55.125 11.045 3.649 19.22 49.005 33.66 29.645 35.28 30.42 30.42 52.020 32.80 25.105 33.60
Tabla 2.2: Valores de los identificadores para cada vector, h, los valores de la variable en el origen del vector y en el extremo final del vector, y por último el valor de la mitad del cuadrado de la diferencia de los contenidos de cadmio para cada vector.
2.2.6 Nubes de las diferencias de las raíces cuadradas. Para evitar el problema planteado en la última parte del apartado anterior Cressie y Hawkins propusieron un estimador del variograma basándose en las raíces cuadradas de los módulos de las diferencias. El estimador es el siguiente: 4
1 1 0.457 + 0.494 2 γ (h ) = ∑ Z ( s i ) − Z (s j ) 2 N (h ) N (h) N (h) Asociado a este estimador, propusieron además una herramienta exploratoria similar a la nube de las raíces cuadradas de las diferencias. Así, la nube de las raíces cuadradas de las diferencias en la dirección e es simplemente un gráfico x - y de puntos 1 de los valores (Z (s i + h e) − Z (s i ) )2 : s i + h e = s j , (i, j ) ∈ N ( he) versus los valores de h. En la Figura 2.18 se representan los boxs plots correspondientes a la nube de las raíces cuadradas de los módulos de las diferencias. A diferencia de la nube de variogramas representada en la Figura 2.14 sólo presenta valores alejados para la fila 6. Esos valores se presentan en la tabla 2.3. h 6 6 6 6
Identificación 51 5 16 14
60 53 56 24
Distancia 64.800 61.846 63.739 64.830
Z(+)
Z(-)
9.9 11.2 8.3 11.5
9.9 6.5 11.0 5.3
(1 2 )[Z ( +) − Z (−)]1 2 0 1.084 0.821 1.245
Tabla 2.3: Valores de los identificadores para cada vector, h, los valores de la variable en el origen del vector y en el extremo final del vector, y por último el valor de la mitad de la raíz cuadrada de la diferencia de los contenidos de cadmio para cada vector.
18
Figura 2.18: Box Plot de la nube de raíces cuadradas de las diferencias. Las unidades en el eje vertical son (p. p. m. de cadmio)1/2.Por ejemplo: el valor h =2 corresponde al intervalo 20 ≤ h < 30 .
2.2.7 El Pocket plot. Las técnicas de análisis de datos presentadas hasta ahora han sido útiles para detectar tendencias, valores alejados y sus influencias sobre la estimación del variograma. Ahora se presentará una técnica para identificar un área localizada que sea atípica con respecto al modelo estacionario impuesto a la medida de dependencia. Ésta explota la naturaleza espacial de los datos a través de las coordenadas de las filas y de las columnas. Aunque se mostró un enfoque robusto para la estimación del variograma, persiste la preocupación de que alguna fracción importante de diferencias {Z (s i + h e) − Z (s i )} sea inapropiada en la estimación de 2γ (he) . Las ubicaciones de la grilla que exhiban medidas diferentes desde el resto necesitan ser identificadas. Estos focos de no - estacionariedad, una vez descubiertos, pueden ser removidos de la estimación del variograma, por supuesto eventualmente deben ser modelados e incorporados en la valuación de los recursos finales. El pocket plot, llamado así porque sirve para detectar focos de noestacionariedad, es una simple idea que se ilustrará a través de un conjunto de datos que se ubican en una grilla como la de la Figura 2.19. Concentrémonos en una fila j de la grilla. Para cualquier otra fila, digamos k, existe un cierto número N jk de diferencias, cuyas ubicaciones están a la distancia
h= k− j
en la dirección Sur - Norte. Si con Y jk denotamos la media de estas 1
diferencias 2 promediadas sobre los N jk términos, y definimos
Yh =
1 N ( h e)
∑
1
Z (s i + he) − Z (s i ) 2
N ( h e)
donde e es el vector unitario en la dirección Sur - Norte. En definitiva Yh es una media ponderada de las Y jk tales que h = k − j .
19
Figura 2.19: Posiciones y datos para ejemplificar el pocket plot.
Se define Pjk = Y jk − Yh El conjunto de estos valores constituye la contribución residual de la fila j al estimador del variograma en los diferentes retardos. Idealmente, estos puntos se distribuirán alrededor del cero, pero si existe algún valor inusual en la fila j, entonces dará una contribución inusual a todos los retardos y típicamente mostrará una distribución de los puntos arriba del cero. Al variar la fila j y ubicando los scatterplots juntos unos de otros se forma el gráfico denominado pocket plot. La tabla 2.4 (pág. 22) presenta los resultados para los datos de la grilla de la Figura 2.19. Los resultados de la última columna son representados en la Figura 2.20. Claramente las filas 5 y 8 son atípicas, los valores {P5k } y {P8 k } están por encima del cero.
Figura 2.20 : Pocket plots en la dirección Sur-Norte
20
En la Figura 2.21 se presenta los box-plots de los pocket plot, que muestra claramente que las filas anteriormente nombradas son los focos de no estacionariedad.
Figura 2.21: box-plots de los pocket plot en la dirección Sur-Norte
Como vemos, las herramientas exploratorias empleadas en Geoestadística son las clásicas para un conjunto de datos cualesquiera, pero se agregan otras que tienen en cuenta la distribución espacial de los datos. Además se incorporan herramientas exploratorias propias de la geoestadística, que sirven para la detección de tendencias, influencia de valores alejados que afectan la estructura de dependencia o la falta de comportamiento estacionario. Éstas últimas explotan la naturaleza espacial de los datos a través de las coordenadas de las filas y de las columnas. Es de destacar que al intentar describir o explorar las estructuras de dependencia espacial se necesita más información que la que se necesita cuando se propone o supone que el modelo subyacente se corresponde con una muestra aleatoria. Los paquetes o software geoestadísticos proveen herramientas exploratorias, como el GEOEAS y Variowin. Las provistas por el GEOEAS, no siempre pueden aprovecharse en su real medida, ya que no tienen una interfaz muy amigable para conectarse con otros programas. Además no presentan todas las técnicas y se hace necesario muchas veces recurrir a la programación para poder usarlas. En este trabajo se ha programado en el MATHEMATICA, generando algunos gráficos como el Scatter plot, el Post plot, el Box-plot de la Nube de Variogramas, el Pocket plot. El VARIOWIN es un soft dedicado solo a una parte de un estudio geoestadístico y sólo presenta herramientas descriptivas relacionadas con el Variograma, y también fue usado aquí, con el fue generada la nube de variogramas y los h-scatter plots. Los software profesionales no específicos para solucionar problemas de índole geoestadístico como el SPPS, STATISTICA, ect. aportan sus instrumentos dedicados a la descripción y la exploración de datos.
21
Tabla 2.4:Valores de las medias de las raíces cuadradas de las diferencias, de las medias ponderadas y de los residuos para las distintas filas j separadas de la fila k en un retardo h.
22
CAPITULO 3: CAMPOS ALEATORIOS. En este capítulo se presentan conceptos teóricos que sustentan la Geoestadística. Se definen los conceptos de campos aleatorios, los momentos de primer y segundo orden de un campo aleatorio, destacando la herramienta fundamental de la Geoestadística: el variograma. Como así también las definiciones necesarias de campo aleatorio estacionario, intrínsecamente estacionario, y el carácter isótropico de los mismos. A los efectos de comparación con el variograma se presenta la función de covarianza estacionaria como así también sus propiedades.
3.1 Campos Aleatorios. En teoría de probabilidad un conjunto de k variables aleatorias Z1 , Z 2 ,..., Z k
define un vector aleatorio Z = ( Z1 , Z 2 ,..., Z k ) con k componentes. De esta manera,
Z(s) = (Z 1 (s), Z 2 (s), L , Z k (s) ) es un vector aleatorio definido en s, con s perteneciente a
un subconjunto D del espacio R d de dimensión d. A cada punto s0 del espacio le corresponde un vector de variables aleatorias Z(s0). Se tiene así una familia de vectores de variables aleatorias que se denomina función de variable aleatoria. El caso más simple es cuando k=1, es decir Z escalar, será tratado en este trabajo. La función aleatoria también es conocida como campo aleatorio o proceso aleatorio. En símbolos un campo aleatorio: { Z(s): s ∈ D}
(3.1.1)
donde D es un subconjunto de Rd, el espacio euclídeo d-dimensional, y s varía continuamente en la región D. Una gran variedad de problemas, tales como los presentados en el capítulo1, pueden ser resueltos usando los métodos geoestadísticos. Lo común a los mismos es que los datos espaciales pueden ser pensados como una realización de la función aleatoria (3.1.1), cuando D se supone no aleatorio y continuo. Matheron denominó variable regionalizada a la realización de un campo aleatorio.
3.2 Función de distribución y momentos de un campo aleatorio. Considérese un campo aleatorio Z(s) definido en D ⊂ Rd. Para m puntos cualesquiera: s1 , s 2 , K , s m ; el vector aleatorio ( Z (s 1 ) , Z (s 2 ), ..., Z (s m )) se caracteriza por su función de distribución m-dimensional: FS
1
, . . . ,Sm (z1
, . . ., zm ) ≡ P[ Z(s1 ) ≤ z1 , . . . , Z(sm ) ≤ zm ]
(3.2.1)
El conjunto de todas estas distribuciones para todo valor de m y para cualquier selección de puntos en D constituye la “ley espacial de probabilidad” del campo aleatorio Z(s). En Geoestadística son suficientes los dos primeros momentos de la distribución de Z(s). De hecho, en la mayoría de las aplicaciones prácticas la información disponible no permite inferir momentos de mayor orden.
23
El momento de primer orden es la esperanza matemática definida como E(Z(s)) = µ(s)
(3.2.2)
Si existe para todo s∈D se denomina la tendencia del proceso aleatorio Z. Algunos autores lo denominan deriva, y otros tendencia y deriva. Los tres momentos de segundo orden considerados en Geoestadística son: n La varianza o momento de segundo orden de Z(s) respecto de µ (s):
σ 2 (s) = Var [ Z (s)] = E {[ Z (s) − µ (s)]2 }
(3.2.3)
La Var[Z(s)] es una función de s.
n La covarianza de dos variables aleatorias Z (s i ) y Z (s j ) , C (s i , s j ) definida como:
{
}
C (s i , s j ) = E [ Z (s i ) − µ (s i )][ Z (s j ) − µ (s j ) ]
(3.2.4)
La covarianza es una función solamente de s i ys j .
n El variograma 2 γ (s i , s j ) que se define como: 2 γ (s i , s j ) = Var [ Z (s i ) − Z (s j )]
(3.2.5)
El semivariograma es por tanto γ (s i , s j ) . Sin embargo, hay autores que usan indistintamente ambos términos para referirse a la función γ (s i , s j ) . En general, no es posible la inferencia estadística a partir de una sola realización, de la misma manera que no es posible determinar la función de distribución de una variable aleatoria (por ejemplo el resultado de tirar un dado) a partir de una sola observación (se ha obtenido 3 en el experimento del dado). Para hacer posible la inferencia estadística en este enfoque, se hace imprescindible introducir hipótesis adicionales acerca de Z(s) para poder reducir el número de parámetros desconocidos a fin de investigar la función de distribución. Estas hipótesis tienen que ver con la homogeneidad espacial de la función aleatoria. Por ejemplo, suponer que la función aleatoria se “repite” en el espacio y que esta “repetición” proporciona la información equivalente a muchas realizaciones de la misma función aleatoria, permitiendo de esta forma la posibilidad de la inferencia estadística.
3.3 Campos aleatorios estacionarios. Se dice que un campo aleatorio es estrictamente estacionario si su función de distribución (3.2.1) es invariante respecto a traslaciones de vector h, o lo que es lo mismo, la función de distribución del vector aleatorio ( Z (s 1 ) , Z (s 2 ), ..., Z (s m )) es
idéntica a la del vector ( Z (s 1 + h) , Z (s 2 + h), ..., Z (s m + h) ) para cualquier h y todo m y todas las ubicaciones posibles s1 , s 2 ,..., s m . Sin embargo, puesto que la geoestadística se basa en los dos primeros momentos del campo aleatorio, es suficiente suponer que estos dos momentos existen y limitar la
24
hipótesis de estacionariedad a los dos primeros momentos. Se dice que un campo aleatorio Z(s) es estacionario de orden 2 o de segundo orden si: a)
E(Z(s)) = µ
para todo s
(3.3.1)
b) Para toda pareja de variables aleatorias Z(s + h), Z(s) su covarianza existe y sólo depende del vector separación h, es decir, C ( s + h , s) = E { Z ( s + h ) Z ( s) } − µ 2 = C ( h )
(3.3.2)
O sea los momentos de segundo orden no dependen de las ubicaciones s, s + h sino del vector de separación entre las ubicaciones: h, comúnmente denominado vector retardo. La función C(h) se denomina covariograma o función de covarianza estacionaria. La estacionariedad de la covarianza implica que la varianza Var [ Z(s)] existe, es
finita y no depende de s, es decir Var [ Z (s)] = C (0).
Además, bajo esta hipótesis el variograma también es estacionario, o sea que el variograma no depende de s ni de s + h, sino únicamente de h, en símbolos:
{
2γ ( s + h, s) = 2γ (h) = E [ Z (s + h) − Z (s)]2
}
(3.3.3)
Podría considerarse que el variograma es repetitivo, redundante e innecesario ya que mide la variabilidad espacial del fenómeno de forma similar a la más conocida función de covarianza. Esto es debido a que a partir de la relación:
Var (Z (s1 ) − Z (s 2 ) ) = Var (Z (s1 ) ) + Var (Z (s 2 ) ) − 2 Cov(Z (s1 ) , Z (s 2 ) ) y el supuesto que Z(⋅) es estacionario de segundo orden, por lo tanto valen (3.3.1) y (3.3.2). Entonces Var (Z (s1 ) − Z (s 2 ) ) se expresa como: Var(Z(s1) − Z(s2)) = 2{ C(0) − C (s 1 − s 2)} haciendo h = s 1 − s 2 ; se observa que el variograma: 2γ(h) = 2{C(0) − C(h)}
(3.3.4)
Es decir, bajo la hipótesis de estacionariedad el semivariograma resulta ser igual a la varianza menos la covarianza.
25
3.4 Isotropía. Una gran simplificación se obtiene al suponer que las estructuras de primer orden y de segundo orden son funciones sólo de la distancia; ello es debido a que en la práctica en la mayoría de los casos el número de datos no es suficiente para caracterizar el comportamiento del campo en las distintas direcciones. La estacionariedad puede ser pensada como una propiedad de invariancia bajo el grupo de transformaciones de las traslaciones de las coordenadas. Para un campo aleatorio en R d podemos considerar también la invariancia bajo rotaciones y reflexiones. Stein(1999) define a un campo aleatorio Z en R d como estrictamente isotrópico si y sólo si sus distribuciones conjuntas finitas son invariantes bajo todos los movimientos rígidos. Esto es, para cualquier matriz ortogonal H d × d y cualquier s ∈ R d Pr{Z (H s1 + s) ≤ t1 , K , Z (H s n + s) ≤ t n }= Pr{Z (s1 ) ≤ t1 , K , Z (s n ) ≤ t n } para todo n finito, para todas las posiciones s1 , K , s n ∈ R d . La condición de isotropía equivale a suponer que no existe razón para distinguir una dirección de otra para el estudio del campo aleatorio bajo consideración. Un campo aleatorio Z(s) en R d es débilmente isotrópico ó isotrópico de orden 2 o de segundo orden si: tal que E(Z(s)) = µ para todo s (3.4.1) a) existe una constante µ b) Para toda pareja de variables aleatorias Z(s + h), Z(s) su covarianza existe y es una función C no negativa tal que depende de la magnitud del vector separación h, es decir, que C (s + h,s)= E{Z (s+h)Z (s)}− µ 2 =C ( h ) para todo s , s + h ∈ R d (3.4.2) La función C ( h ) recibe el nombre de función de autocovarianza isótropica para Z. c) Cuando el covariograma es isotrópico entonces el variograma es isotrópico y vale que: 2γ( h ) = 2{C(0) − C( h )} (3.4.3) Notación: A fin de simplificar la notación simbolizaré h con h. La Figura 3.1 muestra los gráficos donde se observa la equivalencia total entre el semivariograma y el covariograma isotrópicos.
Figura 3.1: Relación entre el semivariograma y la función de covarianza.
26
Un resultado que se deriva fácilmente de las definiciones dadas es el siguiente: un campo aleatorio (estrictamente/débilmente) isotrópico es siempre un campo aleatorio (estrictamente/débilmente) estacionario. Si un campo aleatorio que no sea isotrópico, a través de una transformación lineal de las coordenadas se convierte en isotrópico entonces decimos que el campo aleatorio es geométricamente anisotrópico. Mas detalles, acerca de los tipos de anisotropías serán presentados en el capítulo 4 cuando se dé la manera de detectarlas en forma empírica.
3.5 Meseta y Alcance. Si C (h) → 0 cuando h→∞ en la relación (3.3.4) entonces 2γ (h) → 2C (0) . La cantidad C (0) = γ (∞) es llamada la meseta (the “sill”). Esto significa que C(0) es el valor del semivariograma para cuando el proceso Z(s) es tal que en dos ubicaciones cualesquiera las variables no están correlacionadas. Sí h0 es el menor valor de h para el cual 2γ( h ( 1 + ε) ) = 2C(0) para cualquier ε > 0 se denomina el alcance del variograma: en la dirección h0 / h 0 . En otras palabras la distancia en la cual el variograma alcanza su meseta es llamada el alcance, y marca la zona de influencia en torno a un punto mas allá de la cual la autocorrelación es nula. De (3.3.4) γ (0) = 0 , aunque con frecuencia el semivariograma es discontinuo en el origen, con un salto finito que se llama pepita o efecto pepita (efecto nugget). El alcance, la meseta y el efecto pepita definen las características del variograma. En la Figura 3.2 se representa los mismos para un variograma isotrópico. La razón del efecto pepita a la meseta es a menudo llamada como el efecto pepita relativo y se expresa usualmente en porcentaje.
ALCANCE
M E S E T A
EFECTO PEPITA
Figura 3.2 Parámetros de un semivariograma.
3.6 Campos Aleatorios Intrínsecamente Estacionarios. Como se acaba de ver, para un campo aleatorio estacionario de segundo orden existen la varianza y la covarianza. Sin embargo existen campos aleatorios y fenómenos físicos reales que muestran una capacidad casi ilimitada de variación. Para estos campos no están definidas ni la varianza ni la función de covarianza estacionaria. Por ejemplo en Hidrogeología variables con este comportamiento son el nivel piezométrico en un
27
acuífero con pronunciados gradientes hidráulicos (Fennessy,1982; Neuman y Jacobson, 1984), y la concentración de ciertas especies químicas disueltas en el agua (Myers et al., 1982). Un ejemplo sencillo de función aleatoria no estacionaria es el definido por un proceso de Wiener-Levy en la recta. Este proceso sólo toma valores en puntos discretos 1,2,3, ...k, k+1, ... etc, de forma que los valores en dos puntos vecinos Zk y Zk+1 están relacionados mediante: Z k +1 = Z k + ε k (3.6.1) donde los ε k son variables aleatorias independientes e independientes de los Zk con un valor esperado nulo y varianza unidad. De esta manera, Z k + h = Z k + ε k + ε k +1 + K + ε k + h −1 La varianza de Zk+h viene dada por Var ( Z k + h ) = Var ( Z k ) + h Es decir la varianza del proceso no es estacionaria porque depende de k.
(3.6.2)
(3.6.3)
Si se consideran los incrementos de la variable se tiene que: k +h −1 E ( Z k +h − Z k ) = E ∑ ε i = 0 i =1
(3.6.4) k + h −1 Var ( Z k +h − Z k ) = Var ∑ ε i = h i =1 El valor esperado y la varianza de los incrementos no dependen de k, además el variograma de este proceso es estacionario porque viene dado por h, o sea 2γ (h) = Var ( Z k +h − Z k ) = h
Acabamos de ver que existen campos aleatorios cuya varianza no es estacionaria y sin embargo sus incrementos tienen una varianza que depende de h, pero no de k, es decir el variograma es estacionario. Esta es la motivación para definir el concepto de campos aleatorios intrínsecamente estacionarios o simplemente funciones aleatorias intrínsecas. Z(⋅) es un campo intrínsecamente estacionario si para todo s y s+h perteneciente a D: E[Z(s + h) − Z(s)] = 0
(3.6.5)
Var [Z(s + h) − Z(s)] = 2γ(h)
(3.6.6)
En la literatura habitual (3.6.5) y (3.6.6) definen un proceso de incrementos estacionarios. La condición (3.6.5) es equivalente a decir que la media (esperanza matemática) es constante. Es evidente que una función aleatoria estacionaria de segundo orden es siempre intrínsecamente estacionario, pero como hemos visto, el recíproco no siempre es cierto. La clase de todos los procesos estacionarios de segundo orden está estrictamente contenido en la clase de todos los procesos intrínsecamente estacionarios. Un ejemplo
28
más general de un proceso para el cual 2γ (⋅) está definido pero C (⋅) no lo está, es el movimiento Browniano d- dimensional isotrópico. Si {W(s): s∈ Rd } es tal proceso, entonces Var(W(s+h) − W(s)) = h pero sin embargo: cov (W(u), W(v)) = ½ (u+v−u −v), u , v ∈ Rd no es sólo una función de u − v. Un campo aleatorio se dice intrínsecamente isotrópico sí para todo s y s+h perteneciente a D: E[Z(s + h) − Z(s)] = 0
(3.6.7)
Var [Z (s + h) − Z (s)] = 2 γ ( h )
(3.6.8)
3.7 Variograma vs. Covariograma. Según se suponga que el proceso aleatorio {Z(s): s ∈ D} D⊂ Rd es estacionario de segundo orden o intrínsecamente estacionario, la estructura de dependencia del proceso estocástico quedará especificada por el covariograma o por el variograma en el primer caso y por el variograma en el segundo caso. En los siguientes apartados se mostrarán propiedades generales que deben cumplir ambas medidas de la estructura de la dependencia.
3.7.1 Propiedades del Variograma. El variograma para procesos intrínsecamente estacionarios cumplen con las siguientes propiedades: • 2 γ(h) = 2 γ(−h)
(3.7.1)
• 2 γ(0) = 0
(3.7.2)
El comportamiento del variograma en una vecindad del origen es muy informativo acerca de las propiedades de continuidad del proceso aleatorio Z (⋅) . Los tipos más comunes fueron categorizados por Matheron(1971) según Cressie(1991) como: •
Si 2 γ (⋅) es continua en el origen, entonces Z (⋅) es L2 continua.
•
Si 2 γ (⋅) no se aproxima a 0 cuando h se aproxima al origen, entonces Z (⋅) no es L2 continua y es altamente irregular. (3.7.4)
(3.7.3)
Esta discontinuidad de 2 γ (⋅) en el origen es lo que Matheron denominó efecto nugget. Más precisamente si 2γ(h)→ 2c0 > 0, c0 es el efecto nugget (pepita). Esto es porque pensó que la variación en microescala (pequeñas pepitas) eran las causantes de la discontinuidad en el origen.
29
Nótese que la continuidad L2 de Z(⋅) no significa que las realizaciones son seguramente continuas. El variograma debe satisfacer una propiedad llamada condicional definida negativa. Esto es: •
m
m
∑∑
ai aj 2γ( si − sj ) ≤ 0 para cualquier número finito de ubicaciones espaciales
i =1 j =1
m
{ si : i = 1,2, ..., m } y números reales { ai : i = 1, ... ,m } que satisfacen que
∑a
i
=0.
i =1
(3.7.5) • El variograma debe crecer más lentamente que h . Esto es que: 2 γ (h) lim h → ∞ =0 2 h 2
(3.7.6)
Observación: Algunos autores definen erróneamente el variograma como E(Z(s1 ) − Z(s2))2. Esta definición coincide con la definición 2γ(h)= Var (Z(s1 ) − Z(s2)) con h = s1 − s2 si el proceso Z (⋅) es intrínsecamente estacionario, pero si el proceso Z(⋅) es representado por: Z(s ) = µ(s) + δ(s ) donde δ(s) es un proceso estocástico intrínsecamente estacionario con variograma 2γ(⋅) y la media µ(s) no es constante, entonces E(Z(s1 ) − Z(s2))2 = 2γ(s1 − s2) + (µ(s1 ) −µ (s2))2 que no es en general una función de s1 − s2. Ni necesariamente satisfará la última propiedad citada anteriormente que todos los variogramas deben satisfacer. 3.7.2 Covariograma y Correlograma. El covariograma o función de covarianza estacionaria es relativa al proceso estocástico Z (⋅) estacionario de segundo orden. El covariograma tiene las siguientes propiedades: • C (h) = C (−h)
(3.7.7)
• C(0) = Var[Z(s)] ≥ 0 • C(h) ≤ C(0)
∀ s∈ D ⊂ Rd
(Desigualdad de Cauchy-Schwarz)
• C(h) debe ser definida positiva.
(3.7.8)
(3.7.9) (3.7.10)
Es también verdadero que cualquier función definida positiva corresponde al covariograma de un proceso estocástico estacionario de segundo orden. Este resultado es una consecuencia de la teoría espectral. • En el apartado 3.3 se probó la relación existente entre el variograma y el covariograma cuando el proceso es estacionario:
30
2γ(h) = 2{C(0) − C(h)} Además del variograma y el covariograma se podrían definir el correlograma, C (h) Si C(0) > 0 se define como correlograma a ρ (h) = C ( 0) γ (h) Además de (3.3.4) se obtiene que: ρ (h) = 1 − C ( 0) Es fácil verificar que ρ(h) = ρ(-h) y ρ( 0) = 1. En series de tiempo, estimaciones del correlograma tradicionalmente son usadas por los analistas para diagnosticar la no estacionariedad, la determinación del tipo de dependencia estacionaria, el ajuste del modelo, etc.. En geoestadística no constituye un instrumento más importante que el covariograma o el variograma según corresponda. Anteriormente se vió que el variograma está definido en algunos casos donde el covariograma no lo está, en aquellos casos cuando, en base a los datos se estime ésta última función, se estimaría un parámetro inexistente. Cressie(1991) aporta justificativos teóricos de la preferencia de la estimación del variograma a la del covariograma. Teniendo en cuenta que las funciones aleatorias intrínsecas son más generales que las estacionarias sumado a los resultados citados en el párrafo anterior hablan a favor del uso en Geoestadística del variograma (semivariograma) en lugar del covariograma, aunque la razón principal es la costumbre entre los que habitualmente recurren a las herramientas geoestadísticas para enfocar sus problemas. Por lo tanto en esta tesis se trabajará con el variograma (semivariograma).
31
CAPÍTULO 4: ANÁLISIS ESTRUCTURAL. Se entiende por análisis estructural al proceso de selección del modelo geoestadístico, en el marco de los conceptos definidos en el capítulo anterior. Así, el análisis estructural implica especificar el tipo de hipótesis que se van a hacer sobre la variabilidad del fenómeno en estudio. Es decir, implica decidir si la variable se puede considerar estacionaria, o no; si se requiere la especificación de la tendencia y, en caso de requerirla, la forma que tendrá dicha tendencia; si es suficiente suponer que la variable es intrínseca, etc.. Además de lo anterior se incluye dentro del análisis estructural la estimación del variograma. De hecho, con frecuencia el término análisis estructural se reserva para esta tarea. Sin embargo, la estimación del variograma está tan ligada a la hipótesis sobre el tipo de variable (estacionaria, intrínseca, no intrínseca), que la separación entre ambos procesos resulta un tanto artificial. En este capítulo, se tratará la estimación del variograma, es decir la determinación del variograma empírico. En particular se trabaja en presencia del supuesto de isotropía, y se presenta la corrección en caso de anisotropía geométrica. El variograma empírico es una función que no está definida para cualquier distancia y no da ninguna garantía de cumplir con las propiedades. Es por ello que se presentan los modelos de variogramas teóricos y los criterios de ajuste de los mismos a los variogramas empíricos. Por último se trata la validación cruzada, una forma de medir el ajuste y diagnosticar algunos problemas con el mismo.
4.1 Estimación del variograma. 4.1.1 Método de los momentos. Se supondrá que los datos {Z(si): i = 1, ... , n} pueden modelarse con un proceso estacionario intrínseco. Bajo el supuesto de media constante (3.4.7) un estimador natural basado en el método de los momentos, dado por Matheron es: 2 γ$ (h) ≡
[∑ (Z(si)
− Z(sj))2 ]
/ N(h)
h ∈ Rd
(4.1.1)
donde la suma es sobre N(h) = { (si , sj ): si − sj = h , i,j = 1, ..., n}
(4.1.2)
N(h) es el número de pares distintos contados en N(h). Nótese que N(−h) = {(si , sj): si − sj =− − h , i, j = 1, ..., n} y entonces N (−h)
= N (h)
y
porque N(−h) es coordinable con N(h), y 2 γ$ (h) = 2 γ$ (−h), preservando una propiedad del variograma teórico. Este estimador es insesgado, o sea E(2 γ$ (h)) = 2 γ(h). La prueba es la siguiente, E [[∑ (Z(si ) − Z(sj))2 ]
/ N(h)] = [∑ E(Z(si ) − Z(sj))2 ] / N(h) = [∑ Var(Z(si ) − Z(sj))] / N(h),
debido a que Z
es intrínsecamente estacionario. Por lo tanto: E [[∑ (Z(si ) − Z(sj))2 ]/ N(h)] =[∑2 γ(h)]/ N(h)= 2γ(h)N(h)/N(h)=2 γ(h) 32
El estimador por el método de los momentos en lo sucesivo se denominará estimador clásico. Ejemplo 4.1: Sea los datos sobre un segmento de recta en R1: 2 + 1
3 + 2
4 + 3
3 + 4
5 + 5
4 + 6
6 + 7
5 + 8
6 + 9
6 Z observado + 10
Para h = 1
N(1)= 9
N(1) = {(2,1); (3,2); (4,3); (5,4); (6,5); (7,6); (8,7); (9,8);(10,9)}
2 γ$ (1) = {12 + 12 + 12 + 12 + 12 + 12 + 22 + 12 + 12}/ 9 = 12 / 9 ≅ 1.33 Para h = 2
N(2)= 8
N(2)= {(1,3); (2,4); (3,5); (4,6); (5,7); (6,8); (7,9); (8,10)}
2 γ$ (2) = {22 + 02 + 02 + 22 + 02 + 12 + 12 + 02}/ 8 = 10 / 8 ≅ 1.2 De esta manera se completa la siguiente tabla: h N(h) γ (h)
1 2 9 8 1.33 1.2
3 7 2.1
4 6 3.3
5 5 3.8
6 4 5.7
7 3 8
8 2 9
9 1 16
De esta tabla se puede observar lo siguiente: el número de parejas disminuye al aumentar la distancia h, si bien esto no tiene porque ser siempre así, es común que el número de parejas se reduzca a partir de una cierta distancia. Esto hace que para grandes distancias la estimación del variograma sea poco fiable y limita el máximo valor de h para el que se puede estimar el variograma. Ejemplo 4.2: Sea la grilla en R2:
Si h = (1,0) N((1,0)) = {(s2, s1); (s3,s2); (s5 ,s4); (s6,s5); (s8,s7); (s9,s8)}
33
N((1,0))= 6
La representación esquemática de los vectores que aportan al cálculo del variograma para h = (1,0) es la siguiente:
Con zi representamos el valor de Z(si ) observado en si para i =1,2, ... ,9. Así: 2 γ$ ((1,0)) ={(z2 − z1)2 + (z3 − z2)2 + (z5 − z4)2 + (z6 − z5)2 + (z8 − z7)2 + (z9 − z8)2 }/ 6
Los vectores retardos en las direcciones donde se puede calcular el variograma clásico 2 γ$ (h) son:
h
0 0 0
1 0 1
0 1 1
1 1 2
1 −1 2
0 2 2
2 0 2
2 1 5
1 2 5
1 − 2 5
2 −1 5
2 2 8
2 −2 8
N ( h)
9
6
6
4
4
3
3
2
2
2
2
1
1
h
4.1.2 Variograma para datos irregularmente espaciados. En la práctica, y especialmente cuando se trabaja en dos o tres dimensiones, se presenta el inconveniente de que los datos están distribuidos de forma irregular. En este caso habrá valores de h para los que N(h)será un número muy pequeño con lo cual la varianza de 2 γ$ (h) será muy grande. Esto es una consecuencia de que la varianza del estadístico X (el promedio muestral de una muestra aleatoria de la variable aleatoria X) es inversamente proporcional al tamaño de la muestra y 2 γ$ (h) es un promedio. En este sentido es deseable poder aumentar el número de parejas N(h)empleadas en la estimación de 2 γ$ (h). Para el caso unidimensional en el que h no es más que un escalar se consideran una serie de intervalos de distancia (hj, hj+1) de longitud Lj = hj+1 − hj, generalmente con igual longitud L. Para cada intervalo j se consideran todas las posibles parejas de puntos tales que su distancia xi − yi esté comprendida entre los límites hj y hj+1, es decir,
h j ≤ xi − yi < h j +1 De esta forma se aumenta el número de parejas en desmedro de “discretizar” el variograma (o semivariograma), ya que se obtiene un solo valor del variograma muestral para cada intervalo 2γ * (h *j ) . La distancia h *j puede tomarse como la media de las distancias de todas las parejas de puntos (posiciones) que se emplearon para calcular 2γ * (h *j ) .
34
En dos o tres dimensiones es necesario definir también un ángulo de tolerancia alrededor de la dirección definida por el vector h ya que la consideración de intervalos de distancia es en general insuficiente para garantizar un númeroN(h)grande debido a que los puntos si generalmente no están alineados. La Figura 4.1 nos muestra esquemáticamente como se procede en el plano, para calcular el variograma experimental. El vector separación h viene definido por su ángulo trigonométrico θ y su módulo h= h . A lo largo de la dirección definida por h se determinan una serie de intervalos de longitud L y se considera una zona de tolerancia definida entre las dos direcciones con los ángulos θ − ∆ϕ y θ + ∆ϕ . De esta forma cada intervalo define ahora un cuadrilátero curvilíneo, por ejemplo ABCD, para el que el variograma muestral se calcula tomando todas aquellas parejas de puntos cuyo vector de separación x i − y i cae dentro de dicho cuadrilátero.
C D
A
B
Figura 4.1 Representación esquemática de los intervalos y de la tolerancia para el cálculo del variograma muestral en dos dimensiones.
Puesto que 2γ (h) es una función de la distancia h y de la dirección θ , conviene π π 3π en general calcular 2γˆ (h) en varias direcciones (por ejemplo: θ = 0, , , , en dos 4 2 4 dimensiones) para comprobar si el variograma depende de θ , es decir para contrastar sí el variograma es o no isotrópico.
En forma similar se procede en el espacio tridimensional para la estimación del variograma. En general cuando los datos están irregularmente espaciados en Rd, el estimador clásico del variograma es usualmente suavizado usando en su lugar: 2 γ + (h( k)) = promedio {(Z(si) − Z(sj))2
: (i,j) ∈ N(h) ; h ∈ T(h(k)) }
donde la región T(h(k)) es una región de tolerancia específica en Rd alrededor de h(k), k = 1, 2, ... , r y el “promedio {} ⋅ ” denota un posible promedio ponderado sobre los elementos en {} ⋅ .
35
Las regiones de tolerancia deberán ser tan pequeñas como sea posible para retener la resolución espacial, sin embargo bastante grandes, de modo que el estimador 2γ + (⋅) sea estable. Journel y Huijbregts(1978) (citado por Cressie 1991) recomiendan que el número de pares distintos U {N(h); h ∈ T(h(k))} en T(h(k)) sea al menos 30; así las regiones de tolerancia se deberían elegir de modo que la mayoría de ellas satisfagan esta condición. A menudo las regiones {T(h(k)): k =1,2, ...,r} se eligen disjuntas y exhaustivas, en forma análoga a la elección de los intervalos para construir el histograma de un conjunto de datos univariados. Ejemplo 4.3: Los siguientes gráficos presentan los valores del semivariograma a partir de las observaciones de la variable “contenido de Cadmio” correspondiente al conjunto de datos presentado en el capítulo 2. Los tres primeros gráficos fueron obtenidos con el programa VARIOWIN. En la figura 4.2 se presenta el semivariograma clásico omnidireccional que es una herramienta muy utilizada en geoestadística para obtener estimaciones de las características del variograma: el alcance, la meseta y el efecto pepita. Para obtener el semivariograma clásico omnidireccional se considera la dirección de 00 con una tolerancia angular de 900 a ambos lados de la dirección especificada, de esta manera permite que se incluyan todas las parejas de puntos independientemente de la dirección. Esto, maximiza el número de parejas en cada clase de distancia, pero produce un suavizado del semivariograma.
Figura 4.2:Semivariograma omnidireccional para la variable Cadmio. En este caso se eligió 20 intervalos de clases de distancia con una longitud de 15.00 pies
A pesar de haber elegido 20 intervalos de clase de distancia, en la gráfica se observan 19 valores del semivariograma empírico en vez de los 21 que debería generar el programa VARIOWIN, porque VARIOWIN agrega otra clase (denominada retardo 0) que retiene las parejas con una separación menor que la mitad de la longitud de un intervalo. En esta situación la clase de retardo 0, al igual que la última clase, no está representada porque sólo 2 parejas aportarían para el cálculo del semivariograma empírico, entonces el programa no lo representa. En la figura 4.3 se presenta el mismo semivariograma empírico que en la figura anterior pero además se indican la cantidad de parejas de puntos que permiten el cálculo
36
del valor correspondiente a la distancia promedio representante de cada clase de distancia.
Figura 4.3: Semivariograma omnidireccional para la variable “Cadmio”. Los valores numéricos indican el número de parejas que intervienen en el cálculo.
El soporte de los gráficos anteriores es de 0 hasta la distancia máxima entre las parejas de ubicaciones, es decir hasta 302.362 pies. En general los semivariogramas no son válidos más allá de la mitad de dicha distancia. La Figura 4.4 presenta el semivariograma empírico correspondiente a 10 intervalos de clases de distancia con una longitud de 15.00 pies cada uno.
Figura 4.4:Semivariograma omnidireccional para la variable Cadmio.
En la figura 4.5 se presenta un gráfico del semivariograma empírico producido por el programa VARIO del GEOEAS. A pesar de tomar los recaudos necesarios en cuanto a la elección de los parámetros para este gráfico es bastante similar al de la Figura 4.4 pero no son iguales debido a la manera distinta de tomar los intervalos de tolerancia de las distancias de separación. En este programa los intervalos son de igual tamaño: 15 pies (es una selección realizada), en cambio en el VARIOWIN el correspondiente al retardo cero es de 7.5 pies y el resto de 15 pies.
37
Figura 4.5:Semivariograma omnidireccional para la variable Cadmio producido por el GEOEAS.
En el gráfico de la Figura 4.6 se eligió 5 clases de intervalos de distancia con una longitud de 31.00 pies, se observa una mayor suavización del semivariograma.
Figura 4.6:Semivariograma omnidireccional para la variable Cadmio.
4.1.3 Estimación robusta del variograma. Recordemos desde (4.1.1) que el estimador clásico del variograma es: 2 γ$ (h) ≡ [∑(Z(si ) − Z(sj))2 ]
/ N(h)
h ∈ Rd
donde N(h) = { (si , sj ): si − sj = h , i = 1, ..., n} y N(h) es el número de pares distintos de N(h). Este estimador, está afectado por observaciones atípicas debido al término cuadrático presente en el sumando del estimador clásico. Cressie y Hawkins propusieron un método más robusto para la estimación del variograma. Si Z es un proceso Gaussiano, (Z(s + h) − Z(s))2 se distribuye como
38
2γ(h)χ21 donde χ21 es una variable aleatoria chi cuadrado con 1 grado de libertad. Por lo tanto 2γ(h) es el primer momento de una variable aleatoria altamente asimétrica. Usando las transformaciones de Box y Cox, Cressie y Hawkins encontraron que la raíz cuarta de χ21 tiene una asimetría de 0.08 y una kurtosis de 2.48 (comparada con 0 y 3 de la distribución Gaussiana). De esta forma varios estimadores de posición pueden ser aplicados a: {Z(si) − Z(sj) 1/2: (si , sj ) ∈ N(h) } Finalmente, estos estimadores son elevados a la cuarta potencia, para llevarlos a una escala correcta, y se ajustan por sesgo. Por ejemplo si los estimadores de posición son la media o la mediana, resultan los estimadores: 2 γ (h) ≡
{[∑Z(si) − Z(sj) 1/2 ] / N(h) }4 /(0.457 + 0.494/ N(h))
2 γ~ (h) ≡ [Med{Z(si) − Z(sj) 1/2: (si , sj ) ∈ N(h) }]4/ B(h)
(4.1.3) (4.1.4)
donde Med {} ⋅ denota la mediana de la sucesión {} ⋅ y B(h) corrige el sesgo [el valor asintótico de B(h) = 0.45]. Hay otra ventaja en usar Z(si) − Z(sj) 1/2 en vez de (Z(si ) − Z(sj))2. Nótese que los sumandos en el estimador clásico y en (4.1.3) no son independientes, y tanto más dependientes ellos son, menos eficiente es su promedio para la estimación del variograma. En Cressie (1991) se muestra que los sumandos {Z(si ) − Z(sj) 1/2} en (4.1.3) son menos correlacionados que los sumandos {(Z(si ) − Z(sj)) 2} en (4.1.1). Un estudio comparativo realizado por Omre muestra que el semivariograma muestral así como el estimador de Cressie y Hawkins son bastante sensibles a la hipótesis de que la variable Z (⋅) tenga una distribución normal.
4.2 Modelos de variograma. Los métodos descriptos en la sección 4.1 son útiles para estimar el valor del variograma en un retardo h dado. En la mayoría de las aplicaciones prácticas se requiere sin embargo conocer todos los valores de la función variograma o semivariograma. Como dichas funciones difieren en un factor 2, y los programas trabajan con el semivariograma nos referiremos a este último. Ahora bien, los semivariogramas y las funciones de covarianzas deben satisfacer las propiedades ya citadas. En especial, no cualquier función puede servir de semivariograma ya que ha de ser condicionalmente definida negativa y tener valor nulo en el origen. A las funciones que cumplen estas condiciones se les suele denominar modelos válidos del semivariograma. El semivariograma empírico es una función que no está definida para todo h y no da ninguna garantía de cumplir con las propiedades. En la práctica, lo que se hace es calcularlo y ajustarlo a algún modelo de semivariograma.
4.2.1 Modelos de semivariogramas isotrópicos. Dado que h es un vector y γ es una función escalar, en general γ puede depender h de la distancia h = h como de la orientación u = . Es decir γ(h) es anisótropa. Como h se vió en la sección 3.4 una gran simplificación se obtiene al suponer que las estructuras
39
de dependencia son funciones sólo de la distancia; ello es debido a que en la práctica en la mayoría de los casos, el número de datos no es suficiente para estimar la anisotropía y se suele realizar la hipótesis de que γ (o C) es independiente de la orientación de h, es decir γ es isotrópico. Los tres modelos isotrópicos básicos dados en términos del semivariograma son: el lineal, el esférico, y el exponencial. Además, los modelos de semivariogramas pueden ser convenientemente clasificados en dos tipos; aquellos que alcanzan una meseta y aquellos que no. Los modelos del primer tipo son a menudo denominados modelos de transición. Como se definió en la sección 3.5 la distancia a la que un variograma alcanza la meseta es el alcance. Existen modelos de transición que no alcanzan una meseta efectivamente, pero si se aproximan asintóticamente a un valor constante, que se puede considerar como meseta. En estos modelos se habla de un alcance "práctico" cuando el 95% de ese límite sea alcanzado. En todos los modelos el parámetro c0 representa el efecto pepita. Si c0 = 0 significa que el semivariograma es continuo en el origen. • Modelo lineal ( válido en Rd, d ≥ 1):
h=0
0 γ(h;θ θ)= c + b h 0 l
(4.2.1)
h≠0
donde el vector θ = (c0 , bl )′ con c0 ≥ 0 y bl ≥ 0 . Si el parámetro bl = 0, estamos en presencia de un modelo indicativo de un fenómeno sin ninguna autocorrelación espacial, y se denomina modelo de efecto pepita puro. El modelo lineal es lineal y diferenciable en sus parámetros, los demás modelos no son lineales en sus parámetros aunque sí diferenciables.
a)
b)
Figura 4.7: a) Efecto pepita pura
b)Modelo de semivariograma lineal
• Modelo esférico(válido en Rd, d =1, 2, 3):
40
0 si h= 0 3 3 h 1 h γ(h;θ θ)= c0 + ce − si 0 < h ≤ a e 2 a e 2 a e c0 + ce si h ≥ ae donde el vector θ = ( c0 ; ce ; a e ) l con c0 ≥ 0 , ce ≥ 0 y ae > 0.
(4.2.2)
Este modelo se caracteriza porque alcanza la meseta c0 + ce para una distancia
. finita h = a e el alcance. La pendiente en el origen es igual a 15
c0 + ce ae
.
Este modelo es indicativo de fenómenos continuos aunque no derivables, es decir, fenómenos cuya representación pueden presentar quiebros.
Figura 4.8: Modelo de semivariograma esférico.
• Modelo exponencial (válido en Rd, d ≥ 1): 0 si h=0 γ(h;θ θ)= h si h≠0 c0 +c p 1− exp − a p donde el vector θ = (c0 , c p , a p )′ con c0 ≥ 0 , cp ≥ 0 y ap > 0. El modelo exponencial alcanza la meseta en forma asíntotica:
Como
alcance
h lim c0 + c p 1 − exp − = c0 + c p h →∞ a e práctico a ′p suele tomarse la
distancia
(4.2.3)
a
la
cual γ ( h ) = 0.95 (c0 + c p ) , que es aproximadamente a ′p ≅ 3a p . La pendiente en el origen es
c0 + c p ap
que es menor que la de un semivariograma esférico con el mismo parámetro de
alcance, pero que es mayor a igualdad de alcance práctico ( a ′ ). Es decir, a igualdad de alcance práctico ( a ′ ) el semivariograma exponencial “eleva” y se aproxima a la meseta más rápidamente que el esférico.
41
Figura 4.9: Modelo de semivariograma exponencial.
Otros modelos de semivariogramas usados son: • modelo Gaussiano (válido en Rd, d ≥ 1): 2 h γ(h;θ θ)= c0 1− exp( − 2 ) a donde el vector θ = (c0 , a)′ con c0 ≥ 0 y a > 0 .
(4.2.4)
También alcanza su meseta asintóticamente, lim γ (h; θ ) = c0 , por lo tanto h →∞
hablando estrictamente, no tiene alcance, pero su alcance práctico es a ′ ≅ 3 a , valor para el cual el semivariograma es igual a 0.95 c0. Su comportamiento en el origen de tipo parabólico con pendiente nula es indicativo de una gran regularidad.
Figura 4.10: Modelo de semivariograma Gaussiano.
• modelo racional cuadrático (válido en Rd, d ≥ 1):
0 si γ(h;θ)= c + c h 2 / (1 + h 2 / a ) si r 0 r
42
h= 0 (4.2.5)
h ≠0
Figura 4.11: Modelo de semivariograma racional cuadrático.
Este modelo también alcanza su meseta asintóticamente lim γ (h;θ )=c0 + c r ⋅ a r , por lo h →∞
tanto hablando estrictamente, no tiene alcance, pero su alcance práctico es 17 a ′r = a r , valor para el cual el semivariograma es igual a 0.95 (c0 + c r ⋅ a r ) . 3 • modelo potencial (válido en Rd, d ≥ 1):
0 γ(h;θ)= c + b h 0 t
si h = 0 (4.2.6) λ
si h ≠ 0
donde el vector θ = (c 0 , bt , λ )′ con c0 ≥ 0 , bt≥ 0 y 0 ≤ λ < 2. Nótese que si λ = 0 corresponde a un efecto pepita puro, y este parámetro no puede ser igual o mayor que 2, pues no se satisfaría la propiedad (3.5.6). Estos semivariogramas no tienen meseta, porque tienden a infinito cuando lo hace h = h . Este modelo es de interés especial ya que posee un amplio rango de comportamientos en el origen dependiendo del valor de λ .
Figura 4.12: Modelo de semivariograma potencial para distintos valores de λ .
Los estudios de Yaglon y Cristakos validan esta familia de modelos de semivariogramas, los cuales corresponden a un movimiento Browniano isotrópico fraccional.
43
Construcción de otros modelos de variogramas en Rd. A partir de los modelos presentados, se pueden obtener otros modelos sobre la base de las propiedades siguientes: • Todo modelo de variograma isotrópico en Rd es también válido en Rn mientras d sea mayor a n. • Si 2γ 1 (⋅) y 2γ 2 (⋅) son variogramas válidos en Rd , para b> 0 entonces: n 2 γ (⋅) ≡ 2 γ 1 (⋅) + 2 γ 2 (⋅) es un variograma válido en Rd n b[2γ 1(⋅)] también es un variograma válido en Rd . Por lo tanto la familia de variogramas válidos constituye un cono convexo. En forma más simple las dos últimas propiedades se pueden expresar de la siguiente forma: Los modelos de variogramas se pueden combinar linealmente (con coeficientes positivos) para obtener otras funciones que son variogramas válidos. El modelo resultante, denominado modelo anidado, se puede expresar como
2 γ (h) =
n
∑ b 2 γ (h) i
i
i =1
con bi > 0 y 2 γ i un modelo de variograma válido para todo i. Los modelos presentados en esta sección son modelos isotrópicos, que dependen solamente de la magnitud del vector separación h, en la sección siguiente veremos como se procede cuando no se puede suponer que el campo es isotrópico.
4.2.2 Anisotropía. Cuando la dependencia entre Z(s) y Z(s + h) es una función de la magnitud y la dirección de h, el variograma no es más una función solamente de la distancia entre dos ubicaciones espaciales, entonces el proceso Z se denomina anisotrópico. Las anisotropías son causadas por algún proceso físico subyacente que se desarrolla diferencialmente en el espacio. Por ejemplo la presencia del campo gravitatorio provoca que el proceso en la dirección vertical sea diferente de aquel en las direcciones horizontales. En geología, se tiene otro ejemplo si una mineralización en las rocas ocurre en lentes rectangulares, entonces el variograma será diferente en varias direcciones horizontales. Para estudiar la presencia de anisotropía es necesario calcular el semivariograma en varias direcciones, lo cual suele requerir una cantidad de datos muy superior a lo normalmente disponible. Si esto es posible, puede dibujarse cada semivariograma separadamente (Ver figura 4.13), si los semivariogramas son marcadamente distintos, hay que pensar en la presencia de anisotropía.
Figura 4.13: Semivariogramas en distintas direcciones.
44
Sin embargo, antes de estudiar la anisotropía es necesario tener en cuenta la posibilidad de que el comportamiento direccional puede ser consecuencia de que el proceso Z no es estacionario o ni siquiera intrínseco. La anisotropía puede ser de dos tipos genéricos: elíptica y zonal. La anisotropía elíptica se manifiesta en que el alcance varía con la dirección (figura 4.14). La anisotropía zonal se presenta cuando existe una fuente de variabilidad adicional en una de las direcciones, lo que se manifiesta en que la meseta depende de la dirección.
Figura 4.14: Semivariogramas que indican la existencia de anisotropía elíptica.
El caso de anisotropía elíptica también denominada anisotropía geométrica, se puede tratar mediante una transformación lineal de coordenadas. Para este caso se denomina razón de anisotropía al cociente entre el alcance mínimo y el máximo y las direcciones correspondientes se denominan principales. En el caso de un proceso definido en el plano, la transformación lineal entre las coordenadas originales (x, y) y las nuevas ( x ′, y ′) , en las que el semivariograma es isótropo, viene dado por x ′ λ 0 cos φ = y ′ 0 1 − sen φ
sen φ x cos φ y
(4.2.7)
donde λ es la razón de anisotropía y φ es el ángulo formado por los ejes de la elipse de anisotropía (direcciones principales) con los de las coordenadas originales (figura 4.15). Realizando el cambio de referencias indicado en (4.2.7) se pueden emplear todas las ecuaciones con el semivariograma isótropo.
Figura 4.15: Transformación de coordenadas para corregir la anisotropía elíptica.
45
En general, algunas veces la anisotropía puede ser corregida por una transformación lineal del vector de retardo h. Esto es, el variograma anisotrópico de Z es geométricamente anisotrópico. Es decir: 2γ(h) = 2γο (Ah ) h ∈ Rd (4.2.8) donde A es una matriz d x d y 2γο es una función de una variable real. En la práctica, aunque los semivariogramas direccionales se calculen con frecuencia, no es tan común tener en cuenta la anisotropía en la predicción porque requiere más datos de los habitualmente disponibles. Por ello, si los semivariogramas muestrales direccionales muestran anisotropías, lo adecuado es examinarlas a nivel conceptual, lo que puede sugerir las causas que las expliquen.
4.3 Ajuste a modelos de variogramas. Los distintos estimadores de variogramas, como los presentados en este capítulo, 2 γ$ (⋅), 2 γ (⋅) y 2 γ~ (⋅), no pueden ser usados para la predicción espacial. Ellos no son necesariamente condicionalmente definidos negativos; la ausencia de esta propiedad puede resultar en desconcertantes errores de predicción en media cuadrática negativos. La idea fundamental es buscar un variograma válido que, como una medida de la dependencia espacial, represente más ajustadamente a la dependencia espacial presente en los datos: Z = (Z(s1), ..., Z(sn) )’. El espacio de todos los variogramas válidos es un gran conjunto, usualmente se elige una familia paramétrica de variogramas. Por ejemplo, supongamos que se elige la familia de variogramas isotrópicos lineales: {2γ : γ(h) = c0 + bl h ; c0 ≥ 0, bl ≥ 0 } (4.3.1) Entonces se busca un elemento del conjunto (4.3.1) que sea el que mejor ajuste a los datos, en algún sentido. En general, sí: (4.3.2) P = {2γ : 2γ(⋅) = 2γ(⋅ ;θ ) ; θ ∈ Θ } es el subconjunto de variogramas válidos requerido, entonces se busca un elemento de dicho subconjunto que mejor ajuste al variograma muestral. Diversos criterios de bondad del ajuste para obtener el mejor elemento de P se han propuesto. El criterio de Máxima Verosimilitud requiere supuestos acerca de la distribución de Z, en cambio los criterios basados en Mínimos Cuadrados no requieren de dichos supuestos para la estimación de θ . Otro método de ajuste es el denominado ajuste a sentimiento, el cúal no garantiza un modelo de variograma único ya que se basa en apreciaciones subjetivas y en la experiencia del usuario. Además no permite conocer la calidad y el grado de fiabilidad del modelo ajustado. No obstante se presentarán criterios estadísticos que ayuden a comparar el grado de bondad de un modelo de variograma frente a otro.
46
4.3.1 Máxima Verosimilitud (ML). Supongamos que los datos Z son multivariados Gaussianos; y teniendo en cuenta el modelo lineal: Z=Xβ + δ donde X es una matriz nxq no estocástica, de rango q < n, β es un vector de parámetros fijos (qx1), δ es un vector aleatorio (nx1) El vector δ es multivariado Gaussiano con E(δ δ) = 0 y matriz de varianzas n x n ∑(θ θ ), donde ∑(θ θ ) = (Cov[δ δ(si), δ(sj)]) = (Cov[Z(si), Z(sj)] ) depende de θ como consecuencia de que 2γ depende de θ. Por lo tanto Z es Gaussiano con esperanza E(Z) = X β y matriz de varianzas y covarianzas ∑(θ θ). La verosimilitud de Z es: 1 ( 2π )− n / 2 ∑ (θ θ ) − 1/2 exp { − ( Z − X β )’ [∑(θ θ )]−1 (Z − X β ) } 2 y el negativo del logaritmo de la verosimilitud es: L(β β ,θ θ ) = (n/2) log(2π) + (1/2) log∑ (θ θ ) + (1/2) ( Z − X β )’ [∑(θ θ )]−1 (Z − X β ) β ∈ Rq , θ∈Θ Θ. Los estimadores M.L. β* y θ* satisfacen: L(β β *,θ θ* ) = inf { L(β β ,θ θ ) : β ∈ Rq , θ∈Θ Θ} La obtención de las estimaciones se realiza mediante procedimientos numéricos computacionales.(Consultar Cressie 1991, pág. 472). La estimación por máxima verosimilitud presenta el problema de que los estimadores son altamente sesgados. Otra opción, es filtrar los datos de tal manera que la distribución conjunta ya no depende de β , esto es lo que se conoce con el procedimiento de máxima verosimilitud restringida.
4.3.2 Mínimos cuadrados. El método de ajuste de variogramas de Máxima Verosimilitud ignora la apariencia visual del gráfico del variograma muestral, digamos, {(h, 2 γ$ (he) ) : h = h(1),h(2), ... ,h(K)}, y el ajuste a una curva de variograma teórico que esté cercano a él. Independientemente del método de ajuste, la comparación del gráfico del variograma experimental y los valores del variograma teórico a ajustar es una herramienta de diagnóstico invaluable y altamente recomendada. Para medir la cercanía entre el variograma muestral y el variograma teórico ha sido propuesto, la suma de cuadrados de las diferencias entre un estimador de variograma genérico 2γ#(he) y un modelo 2γ(he; θ). El método de mínimos cuadrados ordinarios especifica que θ se estime por minimización de: K
∑ { 2γ#(h(j) e) − 2γ(h(j) e; θ) }2
(4.3.5)
j =1
en alguna dirección e. A pesar de tener una muy buena interpretación geométrica θo# el estimador de θ mediante (4.4.5) no tiene en cuenta la variación y covariación distribucional del estimador genérico 2γ#.
47
Ajuste mediante mínimos cuadrados generalizados. El método de mínimos cuadrados ordinarios es puramente un procedimiento numérico que tiene una interpretación geométrica actrativa. Para retener la geometría e introducir además el concepto de covariación en el procedimiento, consideramos el criterio de mínimos cuadrados generalizados. Supongamos que un estimador de variograma 2γ# se obtiene en K retardos h(1), h(2), ..., h(K) , donde K es fijo y la cantidad de datos que contribuyen a la estimación en cada retardo es grande,(al menos 30 pares de acuerdo a lo sugerido por Journel y Huijbregts). Además, sea 2γ(h;θ θ) un modelo de variograma cuya forma exacta es conocida excepto por el vector de parámetros desconocidos θ. Sea 2γγ# un vector aleatorio Kx1, 2γγ# = (2γ# (h(1)), 2γ# (h(2)) , ..., 2γ# (h(K)))’; cuya matriz de varianzas y covarianzas es Var(2γγ# ) = V, que puede depender de θ. Entonces se elige el valor de θ que minimiza: θ))’V−1 (2γγ# − 2γγ(θ θ)) (2γγ# − 2γγ(θ
(4.3.6)
donde 2γγ(θ θ )≡ ( 2γ(h(1); θ) , 2γ(h(2); θ) , ... , 2γ(h(K); θ) )’ es el modelo teórico evaluado en los retardos h(1), h(2), ..., h(K). Llamemos al estimador θv#. Además del estimador por mínimos cuadrados ordinarios θo# y el estimador por mínimos cuadrados generalizados θv#, se define el estimador mínimos cuadrados ponderados θ∆#, donde ∆ ≡ diag { var ( 2γ#(h(1))), var ( 2γ#(h(2))), ... , var ( 2γ#(h(1))) } (4.3.7) es una matriz diagonal con las varianzas especificadas a lo largo de la diagonal. Es decir θ∆# es el estimador que se obtiene al minimizar (2γγ# − 2γγ(θ θ))’ ∆ −1 (2γγ# − 2γγ(θ θ)) Al aplicar mínimos cuadrados generalizados no se hace ningún supuesto acerca de cual es la distribución de los datos. Carroll y Rupert mostraron que poseen mejores propiedades de robustez que el estimador máximo verosimil cuando la distribución de Z está mal especificada. La determinación de V en (4.3.6) no es siempre fácil. Cressie(1991) propone como encontrar V en el caso del estimador clásico.
48
4.3.3 Ajuste a sentimiento. El método de los mínimos cuadrados produce un ajuste basado unícamente en el número Ni de parejas sin tener en cuenta ciertos aspectos cualitativos del semivariograma. Es decir, posiblemente deje de lado el aspecto crucial de la representación adecuada del semivariograma cerca del origen. El método de ajuste a sentimiento consiste en seleccionar los parámetros del semivariograma teniendo en cuenta una serie de consideraciones de tipo cualitativo tales como (Clark, 1979): 1. Basta con que el modelo ajustado refleje los principales aspectos del semivariograma muestral (la anisotropía, varianza, etc.). No se deben intentar ajustar los mínimos detalles ya que en general éstos no son una característica del verdadero semivariograma sino más bien fluctuaciones muestrales. 2. El comportamiento de γ * (h) a grandes distancias junto al conocimiento de la varianza muestral s2 determinaran la presencia o no de una meseta S. 3. El valor del efecto pepita puede ser obtenido extrapolando los primeros puntos del semivariograma muestral hasta cortar el eje de ordenadas. 4. En general el ajuste del modelo al semivariograma muestral puede mejorarse considerando modelos compuesto del tipo: γ (h) = ∑ γ i (h) donde cada uno de los i
sumandos son modelos básicos (exponencial, esférico, etc.). Este tipo de semivariogramas puede presentarse cuando la variabilidad aleatoria de Z responde al efecto combinado de varios mecanismos que actúan a diferentes escalas. 5. El sentido común y el conocimiento físico del fenómeno o de la variable que se estudia, son esenciales a lo largo de todo el proceso de estimación del semivariograma. Las apreciaciones subjetivas y la experiencia del usuario no garantizan un modelo de variograma único. Además no permite conocer la calidad y el grado de fiabilidad del modelo ajustado. No obstante en el próximo apartado se verá que es posible establecer criterios estadísticos que ayuden a comparar el grado de bondad de un modelo de variograma frente a otro.
4.4 Validación cruzada del variograma ajustado. Supóngase que el modelo de variograma 2 γ (h ; θ$ ), h ∈R d se ajustó en base de los datos { Z (s i ): i =1, 2,... n } . Una forma para diagnosticar algunos problemas con el ajuste obtenido es mediante la validación cruzada. La idea básica es borrar algún dato y usar los restantes datos para predecir las observaciones borradas. Por lo tanto el error de predicción puede ser obtenido como el valor predicho menos el valor actual. La repetición de este procedimiento sobre muchos subconjuntos permite el conocimiento de la variabilidad del error de predicción. En un contexto de estimación, el borrado de observaciones para mejorar la inferencia de un parámetro estimable fue llamado Jackknifing por Tukey. En el Jackknifing la elaboración de valores falsos lo hace diferente del enfoque de validación cruzada tomada aquí para predicción espacial. En el capítulo 5 se presenta la técnica de predicción espacial conocida como kriging. A los fines de esta sección, se supone que se conoce la predicción Z$ (s 0 ) del
49
valor Z(s0) en la ubicación s0 ∈ D, junto con una medida de su error de predicción en media cuadrática: σ 2k (s 0 ) . Si el modelo de variograma describe adecuadamente la dependencia espacial implícita en el conjunto de datos, entonces el valor predicho Z$ (s 0 ) debería ser cercano al verdadero valor Z(s0). Idealmente, observaciones adicionales de Z(⋅) podrán ser tomadas para chequearlas, o inicialmente algunos de los datos pueden ser reservados para validar el predictor espacial. Pero probablemente, todos los datos se usan para ajustar el variograma y construir el predictor espacial, y no existe la posibilidad de tomar más observaciones. En este caso el enfoque de validación cruzada puede ser usado. Se obtiene a partir de todos los datos el modelo de variograma ajustado 2 γ (h; θ$ ) ; luego se deja de lado un dato Z(sj) y se lo predice con Z$ (s ) basado en 2 γ (h; θ$ ) y los datos −j
j
restantes. El error de predicción en media cuadrática asociado es σ 2− j (s j ) el cual depende inter alia del modelo de variograma ajustado. La cercanía de los valores predichos a los verdaderos valores puede ser caracterizada en diversas formas; por ejemplo: 1 n Z (s j ) − Z$ − j (s j ) ∑ σ (s ) n j =1 −j j
n Z (s ) − Z$ (s ) 2 j −j j 1 ∑ n j =1 σ − j (s j )
(4.6.1)
1
2
Z (s j ) − Z$ − j (s j ) diagrama de tallos y hojas de: : j = 1,K , n σ − j (s j )
(4.6.2)
(4.6.3)
En todos estos resúmenes, se utilizan los residuos de predicción estandarizadas para obtener una herramienta de diagnóstico del ajuste del modelo de variograma 2 γ (h; θ$ ) . La media en (4.6.1) debería ser aproximadamente 0, la raíz cuadrada de la media de los cuadrados en (4.6.2) debería ser aproximadamente 1, y el histograma en (4.6.3) puede ser examinado por la posible presencia de outliers. Samper y Neuman(1989) suponen que los errores de predicción son Gaussianos con correlaciones despreciables; esto permite la construcción de una verosimilitud, basada en los errores de predicción, para ser minimizado con respecto a θ . Sin embargo, aún si las correlaciones entre los errores de predicción son reconocidas en la distribución conjunta, es insensato pensarlo como una verosimilitud porque los parámetros son usados para definir los “datos”(es decir los errores de predicción).
50
CAPITULO 5: KRIGING. En este capítulo se presentará los fundamentos teóricos mínimos de la predicción espacial que se usará en el desarrollo de este trabajo, poniendo énfasis en el método de kriging. Se fundamenta que el kriging es sinónimo de predicción óptima en algún sentido, y en base a los supuestos necesarios se deducen las ecuaciones del kriging simple y kriging ordinario. En otro apartado se presentan características prácticas de esta metodología de predicción. Por último se introducen los resultados teóricos referentes al kriging lognormal.
5.1 Prediccion Espacial y Kriging. Sea {Z(s): s∈D⊂ Rd} una función aleatoria o proceso aleatorio, como se definió en la sección 3.1, desde la cual n datos Z(s1), Z(s2),..., Z(sn) son recolectados. Los datos se utilizan para realizar inferencias sobre el proceso, y así predecir alguna funcional conocida g({Z(s): s∈D⊂ Rd}) [o, más simplemente, g (Z (⋅) ) ] de la función aleatoria Z (⋅) . Los siguientes son dos ejemplos de funcionales: * La predicción puntual supone que g(Z(⋅)) = Z(s0), donde s0 es una ubicación espacial conocida. * La predicción promedio en un bloque supone que g (Z (.) ) = Z ( B) =
1 Z (s) ds B B B ⊂ D , donde B es un bloque cuya ubicación y geometría son conocidas y cuyo volumen d-dimensional es B.
∫
Con predicción espacial se quiere decir predecir g(Z(⋅)) a partir de los datos Z(s1), Z(s2), ..., Z(sn) observados en ubicaciones espaciales conocidas s1,s2, ..., sn. Ésta terminología abarca las nociones temporales de suavizado (o interpolación), filtrado, y predicción, las cuales cuentan con el orden del tiempo para su distinción. Si se dispone de datos temporales del pasado y del presente, el suavizado se refiere a la predicción de g (Z (⋅) ) en el pasado, el filtrado se refiere a la predicción de g (Z (⋅) ) en el tiempo presente; y predicción se refiere a la predicción de g (Z (⋅) ) en puntos del tiempo del futuro. Kriging es un método de predicción espacial que minimiza el error cuadrático medio esperado, él que depende de las propiedades de segundo orden del proceso Z (⋅) . La palabra kriging es sinónima de “predicción óptima”. En otras palabras, se aplica para hacer inferencias de manera óptima sobre valores no observados del proceso ′ aleatorio Z (⋅) desde los datos Z = (Z (s1 ),K, Z (s n ) ) ; observados en el conjunto de ubicaciones espaciales conocidas {s1,...,sn }. D.G.Krige, un ingeniero en minas de Sudáfrica, en el año 1950, desarrolló métodos empíricos para la determinación exacta de la distribución del grado de mineralización desde distribuciones basadas en el muestreo de grados de mineralización. Sin embargo, la formulación de la predicción espacial lineal óptima no proviene del trabajo de Krige. Ésta, se debe a Matheron, quien en honor a Krige puso el
51
nombre “kriging” al método de predicción que desarrollo en sus trabajos de geoestadística. Al mismo tiempo que la geoestadística fue desarrollada en Ingeniería en Minas bajo G. Matheron en Francia, las mismas ideas eran desarrolladas en Meteorología por L.S.Gandin en la Unión Soviética. La original, y simultánea, contribución de estos autores fue poner la predicción lineal óptima en términos de variogramas, es decir en un contexto espacial. El nombre que Gandin usó para su enfoque fue Análisis Objetivo, y la terminología interpolación óptima en vez de kriging. Notación Denotamos con p (Z; g ) al predictor genérico de g (Z (⋅) ) . Cuando g (Z (⋅) ) = Z ( B) se escribe como p (Z; B ) y en el caso particular de B = {s 0 } como p (Z; s 0 ) Cuando existen muchas superficies posibles desde las cuales elegir, una táctica estándar es tratar con éstas estadísticamente. A la luz de las observaciones Z de Z(⋅), se puede considerar la distribución de Z(⋅) condicionada al vector de datos Z. Esto es, las inferencias sobre el proceso [por ejemplo la predicción de Z(s0)] podrían involucrar esta distribución condicional. Denotamos con L(Z(s0), p(Z, s0)) a la pérdida incurrida cuando Z(s0) se predice con p (Z, s 0 ) . Un predictor óptimo p es aquel que minimiza el valor esperado de la función perdida: E { L( Z (s 0 ) ; p( Z , s 0 ))} donde E (⋅) denota la esperanza con respecto a la distribución conjunta de Z(s0) y Z. Es un resultado bien conocido que el mejor predictor minimiza E{L(Z(s0), p(Z,s0) ) / Z}, donde E{ . /Z} denota la esperanza (posterior) con respecto a la distribución condicional de Z(s0) / Z. Así, un predictor óptimo de las partes no observadas de Z (⋅) se obtendrán condicionalmente sobre las partes observadas Z, reforzando la discusión anterior. Una medida condicional de pérdida de predicción es E{L(Z(s0), p)/ Z}, a diferencia de la medida no condicional E{L(Z(s0), p)}. Regiones de predicción Para una función de pérdida L dada y predictor p(Z,s0), se definen las regiones de predicción para Z(s0). La región { Z (s 0 ): L( Z (s 0 ) ; p( Z , s 0 )) < k α } es una región de predicción del 100 (1−α )% si se puede elegir una constante kα (la cual podría depender de s0) de modo que P{ L( Z (s 0 ); p( Z , s 0 )) < k α } = 1 − α 0≤α≤1 En muchos problemas de predicción, se usa como función de pérdida el error cuadrático:
L( Z ( s 0 ) ; p ( Z , s 0 ) ) = ( Z ( s 0 ) − p ( Z , s 0 ) )
2
El predictor óptimo resultante, que minimiza E{(Z(s0) − p(Z,s0) )2/ Z} es la esperanza condicional de Z(s0) dado Z: ° p (Z, s0) = E(Z(s0) / Z)
52
La región de predicción del 100 (1 − α )% es el intervalo simétrico: ( p(Z,s0) − kα1/2, p(Z,s0) + kα1/2 ) La discusión precedente demuestra que se necesita la distribución condicional de Z(s0)/ Z. Ésta es calculada desde la distribución conjunta (n+1)-dimensional de (Z (s0), Z’); sin embargo en la práctica, no es posible la estimación de tal distribución desde los n datos disponibles, a menos que se hagan algunos supuestos simplificadores del modelo. El supuesto más simple a hacer es que Z (⋅) sea un proceso aleatorio Gaussiano, porque entonces la distribución conjunta de (Z(s0), Z’) es Gaussiana y E(Z(s0)/ Z) es lineal en Z, dependiendo solo de: µ (s i ) = E (Z (s i ) ); i = 0, K, n , y C(si,sj) =Cov( Z(si), Z(sj) ) 0 ≤ i ≤ j ≤ n. Aún así, hay potencialmente (n + 1) + (1/2) (n + 1) (n + 2) parámetros a estimar desde sólo n datos. Más supuestos simplificadores, tales como: E(Z(s)) ≡µ y C(si,sj) = C*( si − sj ) donde C * (.) es una función definida positiva sobre Rd, permiten realizar procedimientos de inferencia.
5.2 El mejor predictor lineal - Kriging simple. Para la pérdida dada por el error cuadrático, el mejor predictor es E (Z (s 0 ) / Z ) , el cual no siempre es lineal en Z. En vez de preguntar por el mejor predictor, uno podría preguntarse por el mejor predictor lineal, esto es; obtener l1,l2 , ... , ln , k en n
p(Z,s0) =
∑ li Z(si) + k, i=1
tal que minimice E(Z(s0) − p(Z, s0))2. 2
n Es decir, se debe minimizar E Z (s 0 ) − ∑ li Z (s i ) − k con respecto a l1, l2, ... ,ln, k. i =1 Se puede expresar: 2
n n n E Z (s 0 ) − ∑ li Z (s i ) − k = var Z (s 0 ) − ∑ li Z (s i ) + µ (s 0 ) − ∑ li Z (s i ) − k i =1 i =1 i =1 donde µ(s) =E(Z(s)) s∈ D.
2
n
Si se elige k 0 = µ (s 0 ) − ∑ li Z (s i ) el segundo sumando alcanza su valor mínimo: 0. i =1
Además eligiendo l ′ = c′ Σ −1 se minimiza el primer sumando, donde l ′ = (l1 , l 2 , K , l n )′ , ′ c ≡ (C (s 0 , s1 ),K, C (s 0 , s n ) ) y ∑ es una matriz n x n cuyo elemento (i, j ) es C (s i , s j ) .
[
]
Σ −1 existe con seguridad si el proceso Z es tal que Cov Z (s i ), Z (s j ) = C (s i − s j ) porque al ser esta función definida positiva, el determinante de la matriz del sistema es no nulo entonces l ′ = c′ Σ −1 . Nótese que esta matriz no depende de s por lo que es fácil calcular los coeficientes l i para distintos puntos s simplemente cambiando el vector de términos independientes c. 53
Por lo tanto el predictor lineal óptimo Z* (s0) es: Z * (s 0 ) = p * (Z; s 0 ) = c′ Σ −1 (Z − µ ) + µ (s 0 ) ′ donde µ ≡ (µ (s1 ),L, µ (s n ) ) .
(5.2.1)
El error cuadrático medio de predicción minimizado, a menudo denominado varianza de la predicción es: σsk2(s0) ≡ E(Z(s0) − p* (Z,s0) )2 cuya expresión es σsk2(s0) = Var (Z(s0)) − c’ ∑−1 c = C(s0, s0) − c’ ∑−1c
(5.2.2)
Por lo tanto la varianza de predicción es menor que la varianza de la variable a predecir. Las justificaciones matemáticas de la obtención de las expresiones (5.2.1) y (5.2.2) paso por paso se encuentran en el apéndice C5. Observaciones: • Matheron llamó a tal predicción espacial kriging simple porque cuenta con el conocimiento de la función media µ (.) . Si µ (.) es desconocida (5.2.1) ya no es un predictor.
• El predictor (5.2.1) tiene un interés fundamentalmente teórico, porque en general no se conocen los n+1 valores esperados µ (s 0 ) ; µ ( s1 ) ;K; µ (s n ) ni la matriz de autocovarianzas. • La predicción en s0 requiere de la inversión de una matriz n x n. • Debido a que los datos aparecen linealmente en (5.2.1), el predictor p no es resistente a los outliers. • El predictor p es óptimo entre todos los predictores lineales insesgados. Cuando Z proviene de un proceso Gaussiano es el mejor entre todos los predictores (lineales o no). Está última afirmación se justifica en el siguiente apartado. Datos Gaussianos y datos no Gaussianos. Sea u una variable aleatoria y v un vector aleatorio distribuidos normalmente, con valores esperados E(u) y E(v), respectivamente. Sea Var(u) la varianza de u, CVV la matriz de covarianzas de v y CuV el vector de covarianzas cruzadas entre u y v, el valor esperado y la varianza de u condicionada a v, vienen dadas por: −1 E ( u / v ) = E (u) + Cu′V CVV [ v − E ( v )] −1 Var ( u / v ) = Var (u) − Cu′V CVV Cu V Si u = Z (s 0 ) y v = Z
E ( Z ( s 0 ) / Z ) = p 0 ( Z , s 0 ) = µ ( s 0 ) + c ′ ∑ −1 [ Z − µ ]
54
que coincide con el predictor obtenido en (5.2.1). Es decir si Z (⋅) es un proceso °
Gaussiano entonces el predictor óptimo p coincide con el predictor lineal óptimo p* (bajo la función de pérdida del error cuadrático). Además la varianza de la predicción coincide con la Var ( Z (s 0 ) / Z ) . Sin embargo, no debe olvidarse que la predicción lineal óptima permite predictores tratables que pueden comportarse bastante mal cuando Z (⋅) está lejos de la normalidad. El supuesto Gaussiano tiene otra característica, la homocedasticidad condicional, esto es, la Var (Z (s 0 ) / Z ) no depende de Z. Intuitivamente, el error cuadrático medio de predicción condicional E{(Z(s0) − p(Z,s0))2/ Z}es una medida más apropiada de la variación del predictor que el error cuadrático medio de predicción no condicional E{(Z(s0) − p(Z,s0) )2 }. Sin embargo, para Z (⋅) Gaussiano, p * (Z; s 0 ) = p o (Z; s 0 ) = E ( Z (s 0 ) / Z) , y la medida condicional es indistinguible de la incondicional. No todos los datos se comportan como realizaciones de un proceso Gaussiano. Algunas veces es necesario transformarlos para que eso ocurra. El predictor óptimo E(Z(s0) / Z) bajo el supuesto de modelo no Gaussiano es típicamente no lineal. Puesto que el error cuadrático medio de predicción condicional se puede expresar como: E{(Z(s0) − p(Z,s0) )2/ Z} = Var(Z(s0)/ Z) + {E(Z(s0)/ Z) − p(Z,s0) }2
Así, si los dos primeros momentos de la distribución condicional son conocidos o pueden ser estimados, entonces el error cuadrático medio condicional puede ser ° obtenido [también el predictor óptimo p (Z, s0) = E(Z(s0) / Z)]. °
Para un proceso Gaussiano Z (⋅) : p* (Z; s0) = p (Z, s0) y E{(Z(s0) − p*(Z,s0) )2/ Z}= Var(Z(s0)/ Z) = σsk2(s0); donde σsk2(s0) es el error cuadrático medio de predicción (incondicional). Por otra parte,
{(
)
}
{
}
E Z (s 0 ) − p ∗ (Z, s 0 ) / Z = Var (Z (s 0 ) / Z ) + p o (Z, s 0 ) − p ∗ (Z, s 0 ) 2
2
En situaciones no Gaussianas, donde se use un predictor de kriging simple, la expresión precedente es el error cuadrático medio de predicción condicional y su esperanza produce σsk2(s0). Los procesos que son casi Gaussianos deberían dar origen a predicciones casi lineales.
5.3 Kriging ordinario. Con kriging ordinario se referirá a predicción espacial bajo los supuestos que a continuación se detallan.
Modelo • Z(s) = µ + δ(s) s∈ D, µ ∈ R, y µ es desconocido. • El proceso δ (s) es tal que E (δ (s) ) = 0 para todo s. • El proceso Z (⋅) tiene variograma 2γ (h) = Var (Z (s + h) − Z (h) )
55
(5.3.1) (5.3.2) (5.3.3)
Predictor n
p (Z; B) =
∑
n
λ i Z (s i )
con
∑λ
i
=1
(5.3.4)
i =1
i =1
Esta última condición sobre los coeficientes del predictor lineal garantiza la insesgadez uniforme, es decir: n n n E ( p(Z; B) ) = E λ i Z (s i ) = ∑ λ i E ( µ + δ ( s i ) ) = ∑ λ i µ = µ i =1 i =1 i =1
∑
Obtención de las ecuaciones de kriging ordinario. El predictor óptimo p (Z; B ) se obtiene minimizando el error cuadrático medio de predicción: σk2(s0) ≡ E(Z(s0) − p (Z,s0) )2 (5.3.5) n
sobre la clase de los predictores lineales
∑ λ i Z ( si ) que satisfacen i =1
n
∑λ
i
= 1.
i =1
Se debe minimizar: n n E Z (s 0 ) − λ i Z (s i ) − 2m λi − 1 (5.3.6) i =1 i =1 con respecto a λ 1 , λ 2 , ..., λ n y m, donde m es el multiplicador de Lagrange asociado
∑
n
∑λ
con la restricción
i
∑
= 1.
i =1
n
−
En el apéndice C5 se muestra que minimizar (5.3.6) equivale minimizar (5.3.7) n n (5.3.7) λ i λ j γ (s i −s j )+ 2 λ i γ (s 0 −s i ) − 2m λi − 1 j =1 i =1 i =1 n
∑∑ i =1
∑
∑
Después de derivar (5.3.7) con respecto a λ 1 , λ 2 , ..., λ n y m, e igualar los resultados a cero, se obtiene que los pesos óptimos satisfacen las ecuaciones: −
n
∑ j =1
λ j γ ( s i − s j ) + γ (s 0 − s i ) − m = 0
n
∑λ
i
i = 1, 2, ... , n.
=1
i =1
Esto es, el óptimo λ 1 , λ 2 , ..., λ
λ 0 = Γ0 −1 γ
n
puede ser obtenido desde: (5.3.8)
0
donde
λ 0 ≡ ( λ 1 ,λ 2 , L ,λ n ,m)′ γ 0 ≡ (γ ( s 0 − s 1 ) ,L , γ (s 0 − s n ) ,1 )
56
′
γ (s i − s j ) i = 1,..., n Γ0 ≡ 1 i = n +1 0 i =n +1 Γ0 es una matriz simétrica (n+1) x (n+1).
j = 1,..., n j = 1,..., n j =n+1
Desde (5.3.8) el vector de coeficientes λ*nx1 ≡ ( λ*1 ,K , λ*n ) ′ esta dado por
λ
* nx1
′
′ ( 1 − 1′ Γ −1 γ ) −1 = γ + 1 Γ −1 1 Γ 1 ′
(5.3.9)
y 1 − 1′ Γ −1 γ m= − 1′ Γ −1 1
(5.3.10)
′ ′ donde γ ≡ (γ ( s 0 − s 1 ) ,L , γ (s 0 − s n ) ) , 1′ ≡ (1,L ,1) y Γ es una matriz n x n cuyo elemento genérico es γ (si − sj ) . Las demostraciones de la validez de (5.3.9) y (5.3.10) se la presentan en el apéndice C5. El predictor óptimo (kriging ordinario) es p$ (Z ; s 0 ) = λ*’Z = Z$ (s 0 ) . El error cuadrático medio de predicción mínimo es llamado a veces la varianza del kriging (o de la predicción), a saber: n
n
(
)
σk2 = − ∑ ∑ λ *i λ*j γ s i − s j + 2 i =1 j =1
= 2λ
n
∑λ
* i
γ (s 0 − s i )
(5.3.11)
i =1
γ nx1 − λ *nx1 Γ λ *nx1 = λ *nx1 γ nx1 + λ *nx1 − λ *nx1 Γ λ *nx1 = λ *nx1 γ nx1 + m = λ ′0γ 0 * nx 1
(5.3.12)
Además σk se puede expresar como: 2
σk2 = λ
* nx 1
γ nx1 + m
′ 1 − 1′ Γ −1γ nx1 −1 1 − 1′ Γ −1γ nx1 = γ nx1 + 1 Γ − 1′ Γ −1 1 1′ Γ −1 1 (1′ Γ −1γ nx1 − 1 ) 2 = γ ′nx1 Γ γ nx1 − 1′ Γ −1 1 −1
(5.3.13)
Con estos resultados se puede construir los intervalos de predicción, bajo el supuesto que el proceso sea Gaussiano, el intervalo de predicción de Z(s0) al 95% es:
(
. σ k (s 0 ) , Z$ (s 0 ) + 196 . σ k (s 0 ) I ≡ Z$ ( s 0 ) − 196
57
)
Observaciones: • Se necesita el conocimiento del variograma 2γ (⋅) . • La predicción en s0 requiere la inversión de una matriz n x n. • El predictor p es óptimo entre todos los predictores lineales homogéneos. • p(z,si) = z(si) para i =1,2, ... ,n es decir p es un interpolador exacto. Nótese que la media constante µ no necesita ser estimada. Pero en realidad se puede suponer que la media es una combinación lineal desconocida de funciones conocidas.
Ejemplo 5.3.1: Dado un proceso estocástico en un segmento de recta. Se quiere predecir ′ Z (n + 1) desde el vector de datos Z = (Z (1), Z (2), L , Z (n) ) . Modelo supuesto: • Z(i) = µ + δ(i) i = 1, 2, K , n, K µ es desconocida. • δ(i) es tal que E (δ (i ) ) = 0 ∀ i = 1, 2, K , n, K y ∑δ es la matriz de varianzas y covarianzas cuyo elemento genérico es Cov( Z (i ) , Z ( j ) ) = σ 0 ρ i − j
(
• El proceso Z (⋅) tiene como variograma a 2 γ (i − j ) = σ 0 1 − ρ i − j
)
0 < ρ <1 (5.3.14)
Predictor supuesto: n
p(Z; n + 1) =
∑
λ i Z (si )
n
con
∑λ
i
= 1.
i =1
i =1
Cuyo óptimo es el que minimiza a σk2(n + 1) ≡ E(Z(n + 1) − p (Z, n + 1))2. El predictor óptimo es p$ ( Z ; n + 1 ) = λ * Z donde λ * esta dado por: 1 − 1 ′ Γ −1 γ (5.3.15) 1 ′ Γ −1 1 Para obtener los pesos óptimos se debe encontrar la inversa de la matriz Γ(n × n) . Esta es de la forma:
λ *' = γ ′ Γ −1 + 1′ Γ −1
0 1− ρ 1− ρ 2 ... 1− ρ n −1 0 1− ρ ... 1− ρ n −2 1− ρ 2 Γ= σ 0 ... ... ... ... ... n−2 1− ρ n−3 ... 0 1− ρ 1− ρ 1− ρ n −1 1− ρ n −2 ... 1− ρ 0 que se puede expresar como:
58
(5.3.16)
1 ρ 2 Γ=σ 0 − ... ρ n−2 ρ n −1
ρ 1
ρ 2 ... ρ n −1 1 ρ ... ρ n −2 1
...
...
...
ρ n −3 ρ n−2
... ...
1
1 1
... + ... ... ρ 1 1 1 1 1
ρ
1 ... 1 1 ... 1 ... ... ... o sea ... 1 1 ... 1 1
Γ=σ 0 (−Ω+11′) 2
(5.3.17)
cuya inversa es:
1 −1 Ω −1 11′Ω −1 Γ = 2 −Ω − 1−1′Ω −1 1 σ0 −1
(5.3.18)
Se puede comprobar que: −ρ 1 − ρ 1 + ρ 2 1 ... ... Ω −1 = 1 − ρ2 0 0 0 0
0 −ρ
... ...
...
...
...
1+ ρ
...
−ρ
2
0 0 ... − ρ 1
(5.3.19)
Reemplazando (5.3.19) en (5.3.18), se obtiene que: 1 1 Γ −1 = − 2 1− ρ2 σ0
−ρ 1 − ρ 1+ ρ 2 M M 0 0 0 0
...
0
...
0
M
M
... 1+ ρ 2 ...
−ρ
1− ρ ... 1− ρ 1 1 1− ρ (1− ρ ) 2 ... (1− ρ ) 2 1− ρ 1 M + M M M M M 2 − ρ (1− ρ )(n −1) 1− ρ (1− ρ ) 2 ... (1− ρ ) 2 1− ρ 1 1− ρ ... 1− ρ 1 1 0 0
Reemplazando en (5.3.15) se obtiene el vector de pesos óptimos: Cuya primer componente es 1− ρ λ1 = n−(n−2)ρ
(5.3.20)
La componente i-ésima con i = 2, 3, K , n − 1 es
λi =
(1 − ρ ) 2
(5.3.21)
n − (n − 2 ) ρ
y la componente n -ésima:
λn = ρ +
1− ρ n − ( n − 2) ρ
Entonces el predictor óptimo es p$ ( Z ; n + 1 ) = λ * Z , es decir: 59
(5.3.22)
n −1
p$ ( Z , n + 1 ) = ρ Z (n) + ( 1 − ρ )
Z (1) + ( 1 − ρ ) ∑i = 2 Z (i ) + Z (n)
(5.3.23)
[ n − ( n − 2) ρ ]
Observación: Nótese que si ρ = 0, se obtiene como predictor óptimo a Z . Reemplazando Γ −1 en (5.3.10) se obtiene la siguiente expresión para m; 1− ρ 2 2 m=σ 0 n−(n−2) ρ n
Por lo tanto la varianza del kriging σ 2k = ∑ λ i γ (s 0 − s i ) + m se puede expresar como: i =1
1 (n−nρ − ρ 2 + 2 ρ + 1)+ ρ n−(n−2) ρ
σ k2 (n+1)=σ 0 (1− ρ ) 2
=σ 0
2
(1− ρ 2 )(1+ ρ ) 2 (1− ρ )+ n−(n−2) ρ
¿Qué sucede si se emplea Z en vez del predictor (5.3.23)? A los fines de la comparación se determina el error cuadrático medio de predicción cuando se utiliza Z para predecir Z (n +1)
E ( Z (n + 1) − Z ) = Var ( Z (n + 1) − Z ) = 2
1 ρ2 ρ n 1 = σ 20 1 + 1 + 2 1 − ρ n −1 ρ − −2 2 1− ρ n n(1 − ρ ) n
(
)
A los fines de la comparación se presenta las siguientes tablas para distintos valores de n y de ρ . n
ρ Tabla 5.1: Valores de la varianza del kriging para distintos valores del tamaño de la muestra (n) y del coeficiente de correlación ρ . n
ρ Tabla 5.2: Valores del error cuadrático medio de predicción cuando se utiliza Z (n +1) en vez del predictor de kriging.
Z para predecir
De la comparación de ambas tablas se observa que si ρ = 0.1 , es decir si se presenta correlación baja, independientemente del tamaño muestral ambas medidas no
60
difieren significativamente. Si se considera una correlación media (0.5) los valores son comparables pero los de la varianza del kriging son siempre menores que los correspondientes valores del error cuadrático medio. En cambio en presencia de correlación alta, independientemente del tamaño muestral, ambas medidas son muy distintas. Siendo la varianza del kriging muy menor que el correspondiente error cuadrático medio cuando se utiliza como predictor a la media muestral. Esto habla a favor del uso del estimador de kriging cuando se presenta correlación.
5.4 Aspectos prácticos. 5.4.1 Efectos de distribución de los datos. Una de las ventajas del kriging sobre la mayoría de los interpoladores es la forma en que tiene en cuenta la distribución de los datos. Se esperaría que cualquier predictor asignase pesos máximos a los puntos situados más cerca y que éstos fuesen disminuyendo a medida que aumenta la distancia entre puntos donde se quiere predecir y la ubicación de los datos. Igualmente, parecería razonable exigir que si dos datos están ubicados muy próximos, sus pesos serían menores que si están alejados. Estas propiedades se cumplen en forma natural cuando el kriging es el método de predicción elegido. A los fines de ejemplificar este efecto en las predicciones, se emplearan tres semivariogramas: uno potencial, uno esférico sin efecto pepita y el tercero esférico con efecto pepita, es decir tres modelos que representan diferentes estructuras de dependencia. Así el primer modelo corresponde a una fuerte estructura de dependencia, el segundo a una moderada y la tercera a una débil. En la Figura 5.1 se presentan gráficamente dichos modelos.
Figura 5.1:De izquierda a derecha: modelo potencial, modelo esférico sin efecto pepita, y modelo esférico con efecto pepita 1.
En la figura 5.2 se presentan distintas disposiciones de los puntos para las predicciones con las tres diferentes estructuras de dependencia citadas anteriormente. En todas ellas el punto donde se quiere predecir se encuentra rodeado de los puntos donde se dispone información. En el caso A) estos se encuentran distribuidos simétricamente respecto del punto de predicción y cada uno a una distancia de una unidad de dicho punto. En el caso B) los puntos ya no se distribuyen simétricamente, sino el punto que se encuentra al norte se encuentra a una distancia 1.5 del punto de predicción. En el caso C) la distribución de los puntos para realizar la predicción es simétrica, pero a diferencia de A) la cantidad de los mismos se ha duplicado.
61
A)
B)
C)
Figura 5.2: Distintas disposiciones de los puntos para realizar la predicción en el punto coloreado de azul.
Los pesos obtenidos mediante la solución de las ecuaciones de kriging empleando cada uno de los 3 semivariogramas para el caso A) son iguales a 0.25 para todos los puntos debido a la simetría y a la distancia que cada uno de ellos se encuentra del punto de predicción. En la figura 5.3 se presentan esquemáticamente los pesos para este caso.
1)
2)
3)
Figura 5.3: Pesos de los puntos para el caso A usando para el semivariograma: 1) el modelo potencial. 2) el modelo esférico sin efecto pepita. 3) el modelo esférico con efecto pepita.
Los pesos obtenidos para el caso B) son como los muestra la figura 5.4, para todos los semivariogramas, los pesos correspondientes a los puntos simétricos respecto al punto de predicción tienen valores iguales disminuyendo sus valores a medida que disminuye el grado de dependencia. En cuanto a los correspondientes a los puntos asimétricos se puede observar que presentan valores menores los puntos que están más alejados.
1)
2)
3)
Figura 5.4: Pesos de los puntos para el caso B usando el semivariograma: 1) potencial. 2) esférico sin efecto pepita. 3) esférico con efecto pepita.
62
Para el caso C) lo que llama la atención es la existencia de pesos negativos, cuando se está en presencia de las dos primeras estructuras de dependencia. Esto se debe a que los puntos situados más próximos al de predicción reciben los mayores pesos y reducen, llegando hasta en algunos casos a hacerlos negativos, los de los puntos que están más lejos. ¿Pero que sucede con la tercer estructura, la que presenta el efecto pepita? Los pesos asignados a todos los puntos de predicción no son negativos. Esto responde a que la presencia de efecto pepita tiende a uniformizar los pesos, Cressie (1991). 1)
2)
3)
Figura 5.5: Pesos de los puntos para el caso C usando el semivariograma: 1) potencial. 2) esférico sin efecto pepita. 3) esférico con efecto pepita.
Efecto pantalla. En la figura 5.6 se presentan otras disposiciones de los puntos. En todas ellas el punto donde se quiere predecir se encuentra rodeado de los puntos donde se dispone información. En el caso D) cuatro de estos se encuentran distribuidos simétricamente respecto del punto de predicción y cada uno a una distancia de una unidad de dicho punto mientras que un quinto se encuentra al Este a una distancia de 2 unidades. En el caso E) 4 puntos presentan una distribución similar a la disposición B) con un quinto punto que se encuentra al norte a una distancia 2 del punto de predicción. En el caso F) la distribución de los puntos es similar a la del caso C) pero se agrega un punto en la dirección Noroeste. D)
E)
F)
Figura 5.6: Distintas disposiciones de los puntos para realizar la predicción en el punto coloreado de azul.
Las Figuras 5.7, 5.8 y 5.9 también indican la existencia de pesos negativos. Como se puede observar, siempre que un punto de observación tiene un peso negativo existe otro situado entre aquel y el punto de predicción. Este efecto se denomina efecto pantalla y consiste en que los puntos situados más próximos al de predicción reciben los mayores pesos y reducen, llegando hasta en algunos casos a hacerlos negativos, los 63
de los puntos que están detrás. Así, se dice que los puntos más próximos “apantallan” a los que quedan detrás.
1)
2)
3)
Figura 5.7: Pesos de los puntos para el caso D usando el semivariograma: 1) potencial. 2) esférico sin efecto pepita. 3) esférico con efecto pepita.
1)
2)
3)
Figura 5.8: Pesos de los puntos para el caso E usando el semivariograma: 1) potencial. 2) esférico sin efecto pepita. 3) esférico con efecto pepita. 1)
2)
3)
Figura 5.9: Pesos de los puntos para el caso F usando el semivariograma: 1) potencial. 2) esférico sin efecto pepita. 3) esférico con efecto pepita.
La existencia de pesos negativos tiene su razón de ser, pero puede conducir a resultados extraños, como dar como predicción un valor negativo en casos de que la variable ha de ser necesariamente positiva. Para evitarlo, Svidarovszky (1985) propone un método de kriging en el que se impone que los pesos sean no negativos. Otro argumento es el de Journel (1986) que indica que si el problema se reforma adecuadamente, no es preciso realizar la restricción antes mencionada en forma explícita. Por ejemplo se puede suponer que la variable en cuestión no es intrínseca, o hacer una transformación logarítmica de los datos etc., pero esto requiere hipótesis adicionales. En todos los gráficos de los casos 3, se muestra que la presencia de efecto pepita tiende a uniformizar los pesos y reduce considerablemente el apantallamiento. Si la variable no presenta ninguna estructura (efecto pepita puro), todos los pesos son iguales a 1/n, siendo n el número de datos. En resumen, los pesos serán más uniformes (independientes de la distribución de los puntos de observación) cuanto menor sea la estructura de la variable. Es decir, dependerán mucho de la distribución si
64
el semivariograma tiene pendiente nula en el origen (todos los ejemplos 1 de las figuras), algo menos si tiene comportamiento lineal en el origen (ejemplos 2 de todas las figuras) y mucho menos si existe pepita (ejemplos 3 de todas las figuras). En este contexto, Journel y Huigbregts (1978) hablan del efecto de suavización causado por la pepita. Efecto de agrupamiento. Otro efecto importante de la distribución de los datos es el llamado efecto de agrupamiento, según el cual puntos muy próximos tienden a comportarse como si se agrupasen en uno solo. En la Figura 5.10 existe una pareja de puntos muy próximos. Así, los pesos de los dos puntos muy próximos ubicados al Oeste suman 0.254 para el modelo 1) y los pesos de cada uno de los restantes puntos rondan en 0.25. Para el modelo 2) la suma es 0.264 que difiere del peso de su simétrico en 0.017.En cuanto al 3) la suma es 0.322, mientras que los pesos de cada uno de los restantes puntos rondan en 0.23.
1)
2)
3)
Figura 5.10: Pesos de los puntos usando el semivariograma: 1) potencial. 2) esférico sin efecto pepita. 3) esférico con efecto pepita.
Al igual que en el párrafo anterior, el efecto de agrupamiento es tanto mayor cuanto más estructurada esté la variable. Distribución angular de los datos. Los pesos son también sensibles a la distribución angular de los datos. En la figura 5.11 se presentan distintas disposiciones de los puntos para mostrar la dependencia de los pesos del kriging con la distribución de los datos y los modelos de semivariogramas. G)
H)
Figura 5.11: Distintas disposiciones de los puntos para realizar la predicción en el punto coloreado de azul.
En ambas disposiciones las distancias entre los puntos de observación y el de predicción son iguales. Pero en el primer caso (caso G) el punto de predicción se encuentra rodeado por los de observación. En cambio en el segundo (caso H) el punto de predicción se encuentra en una esquina. 65
1)
2)
3)
Figura 5.12: Pesos de los puntos para el caso G usando el semivariograma: 1) potencial. 2) esférico sin efecto pepita. 3) esférico con efecto pepita.
1)
2)
3)
Figura 5.13: Pesos de los puntos para el caso H usando el semivariograma: 1) potencial. 2) esférico sin efecto pepita. 3) esférico con efecto pepita. La distribución de los pesos es mucho más uniforme en los casos presentados en la Figura 5.12 que en los de la Figura 5.13. Esto es consecuencia de la diferencia en la distribución de los puntos de observación con respecto al de predicción. Las conclusiones del párrafo anterior en relación con la uniformidad de los pesos en función de la estructura siguen siendo válidas. Así, los pesos del caso H 3) (mínima estructura) son mucho más uniformes que los de H 2) y los de éste mucho más que los del H 1) (máxima estructura), en el que un punto tiene un peso de 1.09. Además las varianzas de predicción son mucho más pequeñas en el caso en el que los datos rodean el punto de predicción que en el caso en que éste queda en una esquina. Es evidente, que los valores de los pesos son demasiados sensibles a la posición de los puntos de observación, como para que sean predecibles. Lo importante de la discusión anterior es tener presente que conviene que los datos estén distribuidos lo más uniformemente posible. Ello lleva al campo del diseño de redes de observación, que consiste en la selección de la posición de los puntos donde se mide a fin de obtener una buena predicción. Un criterio es elegirlo de forma que se minimice la varianza de predicción, la cual depende de la distribución de los puntos.
5.4.2 Dibujo de curvas de nivel. Luego de resolver las ecuaciones de kriging, es fácil obtener los valores predichos en cualquier punto s, simplemente cambiando el término independiente. Así, para dibujar curvas de nivel, lo más apropiado es hacer que s vaya recorriendo los nodos de una malla regular. A partir de los valores predichos de Z en la malla, se pueden dibujar las curvas de nivel mediante algún programa de computación, por ejemplo el programa SURFER, que generalmente requiere conocer los valores de Z sobre alguna malla regular o el aún más conocido programa para la investigación de datos geoestadísticos el GEOEAS. En forma similar que se obtienen las curvas de nivel de las predicciones para la variable se obtienen simultáneamente las curvas de nivel de las estimaciones de los desvíos estándares de la predicción. Así, de esta manera se permite visualizar la incertidumbre de la predicción.
66
Los resultados geoestadísticos se presentan frecuentemente en forma cartográfica. En la figura 5.14 a) se muestra el mapa con los resultados del kriging para el problema de la concentración de Calcio presentado en el capítulo 2. Para la generación del mismo se considero una malla o grilla regular con origen en el punto (260, 120) y un incremento de 20 pies en ambas direcciones de tal manera de abarcar el rectángulo definido por las posiciones de los datos. Mientras que en la figura 5.14 b) se presenta el mapa con los desvíos estándares de la predicción. a)
b)
Figura 5.14:
a) Mapa del kriging del contenido de Cadmio. b) Mapa de los desvíos estándares de predicción del contenido de Cadmio.
67
5.5 Kriging lognormal. El kriging es óptimo cuando el proceso Z tiene una distribución normal o gaussiana, claro que su optimabilidad puede degradarse considerablemente al tratar con otras distribuciones. Para superar esta limitación es conveniente realizar una transformación sobre Z de forma que los transformados tengan una distribución Gaussiana. De esta manera por ejemplo, se define el proceso aleatorio lognormal. Un proceso aleatorio {Z (s): s ∈D} de valores positivos se dice lognormal sí y solo sí la transformación Y (s) ≡ log Z (s) s ∈D es un proceso Gaussiano. (5.5.1) La meta del kriging log-normal, al igual que la del kriging ordinario, es predecir Z(s 0 ) desde el vector de observaciones Z ≡ ( Z (s1 ),K , Z (s n ) ) ′ o en forma más general predecir Z(B). La idea es transformar el problema desde la escala Z a la escala Y. En un principio se supondrá que el proceso aleatorio Y es intrínsecamente estacionario. El predictor de Y(s 0 ) es n
n
i =1
i =1
p$ Y (Z ; s 0 ) ≡ ∑ λ i log Z (s i ) = ∑ λ i Y (s i )
(5.5.2)
donde λ 1,K, λ n se obtienen resolviendo (5.3.8) en la escala Y; es decir, el variograma usado en (5.3.8) es 2γ Y (h)≡ var(Y (s+ h)−Y (s) ), h ∈R d . El predictor exp( p$ Y (Z ; s 0 )) es sesgado para Z(s 0 ) . Bajo los supuestos siguientes: el proceso Y (⋅) es Gaussiano e intrínsecamente estacionario con media µ Y y
variograma 2γ Y (⋅) y con varianza finita σ Y2 (s) ≡ var (Y (s)), s ∈D ; entonces un predictor insesgado para Z(s 0 ) es 1 1 ( p Z ( Z ; s 0 ) ≡ exp p$ Y (Z ; s 0 ) + σ Y2 (s 0 ) − var p$ Y (Z ; s 0 ) 2 2 (5.5.3) 1 2 $ = exp pY (Z ; s 0 ) + σ Y ,k (s 0 ) − mY 2 2 donde σ Y ,k (s 0 ) y mY son, desde (5.3.9) y (5.3.10), la varianza del kriging y el multiplicador de Lagrange en la escala de Y . El error cuadrática medio de predicción es: 2 E (Z (s 0 )− p( Z (Z;s 0 ) ) = exp 2µ Y +σ Y2 (s 0 ) × exp σ Y2 (s 0 ) +exp(var( pˆ Y (Z;s 0 )) ) −
{ (
)} { (
)
−2 exp(cov(Y (s 0 ), pˆ Y (Z;s 0 ) )) Todos estos resultados están presentados en Cressie(1991).
}
(5.5.4)
La principal característica de (5.5.4) es que se necesita conocer el variograma 2 γ Y y los momentos µ Y y σ Y2 (⋅) o tienen que ser estimados; provocando problemas de inferencia difíciles de resolver. A los fines de simplificarlos se puede imponer la condición más fuerte sobre el proceso Y, a saber, estacionario de segundo orden. Entonces σ Y2 (s 0 ) = CY (0) que puede ser estimado desde la relación del variograma con
el covariograma: 2 γ Y (h) = 2( CY (0) − CY (h)) . Además µ Y puede ser estimado por mínimos cuadrados ponderados basados en los datos transformados Y ≡ ( Y (s 1 ),K , Y (s n ) ) ′ . 68
El predictor kriging - lognormal (5.5.3) es óptimo en la clase de predictores que satisfacen la siguiente condición: Minimizar E (Y (s 0 ) − log p(Z ; s 0 ) n
∑λ i =1
i
2
sujeto a p(Z ; s 0 ) = e
∑ in = 1 λ
log Z (s ) + k i i
,
= 1 y E( Z(s 0 )) = E( p(Z ; s 0 )) . Es importante notar que (5.5.3) no coincide con el
mejor predictor insesgado E ( Z(s 0 ) / Z ) .
69
CAPÍTULO 6: KRIGING UNIVERSAL. En el capítulo anterior se presentaron métodos de kriging bajo el supuesto que la variable a estimar es intrínseca. En muchos casos, la variable no satisface estas condiciones, y el fenómeno se caracteriza por presentar una tendencia. Así por ejemplo, en hidrología subterránea, los niveles piezométricos muestran una tendencia global en la dirección del flujo (Samper Calvete – Carrera Ramirez, 1996). Una vez que se ha detectado que Z no es intrínseca, caben entre otras, las siguientes posibilidades: 1) Suponer que la variable es localmente estacionaria. Esto es, aplicar el método de kriging ordinario pero limitándolo a los puntos de observación que se encuentran a una distancia menor a una prefijada (Kriging en un entorno). Esto es posible llevarlo a cabo con un programa como el GEOEAS. 2) Suponer que estas variables se descomponen como la suma de la tendencia, µ(s), tratada como función determinista, y una componente estocástica δ(s) que se puede tratar como función aleatoria intrínseca con valor esperado nulo. 3) Otra alternativa, propuesta por Matheron, se basa en la idea de filtrar linealmente los datos, para ello definió las funciones aleatorias intrínsecas de orden k. En este capítulo se presenta la segunda alternativa, que se denomina "kriging universal”. Este método fue el primero que planteó el problema de la predicción de funciones no intrínsecas de una forma global. También como un caso particular de este tipo de predicción se tratará el método denominado kriging mediana polish.
6.1 Kriging Universal. En esta sección se supondrá que E(Z(s)) ya no es más constante sino una combinación lineal desconocida de funciones conocidas f 0 (s) ,L , f p (s) . Aunque
{
}
cada f i (s) se ha expresado como una función de s, cualesquiera de ellas podrían ser una constante por ejemplo 1, o un valor de una variable explicatoria asociada con el dato en s, s ∈ D. 6.1.1 Modelo supuesto. En esta sección se adopta el modelo: p +1
Z ( s ) = ∑ f j −1 ( s ) β j −1 + δ ( s )
s∈D
(6.1.1)
j =1
donde: p +1
• µ (s) = ∑ f j −1 (s) β j −1 es denominada tendencia o deriva que tiene carácter j =1
determinista. • β ≡ ( β 0 ,L , β p ) ′ ∈R p +1 es un vector de parámetros desconocidos.
70
• δ(⋅) es un proceso aleatorio intrínsecamente estacionario con E(δ(⋅)) = 0 y variograma 2 γ(⋅). El vector de datos Z se puede expresar como:
Z (s1 ) f 0 (s 1 ) Z ( s 2 ) = f 0 ( s 2 ) M M Z ( s n ) f 0 ( s n )
f 1 (s1 ) L f 1 (s 2 ) L
f p −1 (s 1 ) f p −1 (s 2 )
M L f 1 (s n ) L
M f p −1 (s n )
f p (s1 ) f p ( s 2 ) M f p ( s n )
β 0 δ (s1 ) β 1 + δ ( s 2 ) M M β p δ ( s n )
o en forma más compacta: Z = X β +δ
(6.1.2)
Donde X es una matriz n x (p+1) cuya componente (i, j ) es f j −1 (s i ) ; β es el vector anteriormente citado y δ es un vector n x 1 tal que el i- ésimo elemento es δ(si). Si se quiere predecir Z en una determinada ubicación s0, Z(s0) debe satisfacer el modelo, o sea: Z(s 0 ) = x ′ β + δ (s 0 )
(
x ≡ f 0 (s 0 ) , f 1 (s 0 ),L , f p (s 0 )
donde
(6.1.3)
)
′
6.1.2 Predictor Supuesto. A partir del vector de datos Z se desea predecir linealmente Z(s0) usando un predictor uniformemente insesgado. Esto es, el predictor es de la forma: n
p( Z ; s 0 ) = ∑ λ i Z ( s i ) = λ ′ Z
(6.1.4)
i =1
donde λ es un vector de pesos n x 1. La condición necesaria y suficiente para que sea un predictor uniformemente insesgado es: λ ′ X = x′ (6.1.5) ya que E [ p( Z ; s 0 ) ] = λ ′ E ( Z ) = λ ′ E ( X β + δ ) = λ ′ X β (6.1.6) es igual a: E ( Z(s 0 )) = x ′ β
(6.1.7)
para todo β ∈R p + 1 si y solo si λ ′ X = x ′ . Nótese que si p=0 y f0(s) ≡ 1 se obtiene el kriging ordinario y en este caso λ ′ X = x ′ se n
reduce a
∑λ
i
= 1.
i =1
71
6.1.3 Predicción espacial óptima del proceso Z. En el kriging universal, el predictor lineal insesgado óptimo, se expresará como p$ ( Z ; s 0 ) y es aquel que minimiza el error cuadrático medio de predicción:
σ 2e = E [ Z (s 0 ) − p( Z ; s 0 )] sobre λ 1,K , λ n sujeto a λ ′ X = x ′ .
2
(6.1.8)
El adjetivo “universal” fue utilizado por Matheron (1969) para referirse a la insesgadez del predictor cuando la tendencia es una combinación desconocida de funciones conocidas. El problema de optimización puede ser expresado equivalentemente usando los multiplicadores de Lagrange como; Se debe minimizar:
E [ Z (s 0 ) − λ ′ Z ] − 2 m .( λ ′ X − x ′ ) Con respecto a los vectores λ ′ y m; donde m ≡ ( m0 ,K , m p ) ′ . Es decir, se debe minimizar 2
(6.1.9)
2
p +1 n n E Z (s 0 ) − ∑ λ i Z (s i ) − 2 ∑ m j −1 ∑ λ i f j −1 (s i ) − f j −1 (s 0 ) i =1 j =1 i =1
(6.1.10)
con respecto a λ 1 ,K , λ n , m0 ,K , m p . Suponiendo que: 2 γ (h) = Var ( Z (s + h) − Z (s)) = Var ( δ (s + h) − δ (s) ) = E ( δ (s + h) − δ (s) ) 2
(6.1.11)
(6.1.10) se transforma en: n
n
− ∑ ∑ λ i λ j γ (s i − s j ) + 2 i =1 j =1
n
∑λ i =1
p +1
j =1
i =1
n
i γ ( s 0 − s i ) − 2 ∑ m j −1 ∑ λ i f j −1 ( s i ) − f j −1 ( s 0 )
(6.1.12)
Es importante resaltar que minimizar (6.1.12) equivale minimizar (6.1.10) si alguna función del conjunto f j −1 (s): j = 1,K , p + 1 es idénticamente 1. Caso contrario
{
}
las ecuaciones apropiadas para el kriging universal se escriben en términos de las funciones de covarianzas. 6.1.4 Ecuaciones de kriging universal. Derivando con respecto a λ 1,K , λ n , m0 K , m p e igualando a cero, los pesos óptimos son obtenidos desde
λ u = Γu−1 γ u donde
λ u ≡ ( λ 1 ,K , λ n , m0 ,K , m p ) ′
o
λ u ≡ ( λ ′ , m ′)
′
72
(6.1.13) (6.1.14)
o
γ u ≡ (γ (s 0 − s1 ),L , γ (s 0 − s n ) ,1, f 1 (s 0 ),L , f p (s 0 )) ′ γ u ≡ ( γ ′ , x ′) ′
(6.1.15)
y Γu es una matriz simétrica (n + p + 1) x (n + p + 1) definida como sigue:
i = 1,K , n γ (s i − s j ) Γu ≡ f j −1− n (s i ) i = 1,L , n 0 i = n + 1,L , n + p + 1
j = 1,L , n j = n + 1,L , n + p + 1
(6.1.16)
j = n + 1,L , n + p + 1
que se puede expresar en términos de submatrices por: Γ Γu = X ′ donde: • Γ es la matriz n x n , Γ= γ ( s i − s j )
X 0
i , j = 1, 2 ,K , n .
• X es la matriz n x ( p+1) X = f j −1− n (s i )
i = 1,2,K , n j = n + 1,K , n + p + 1
• 0 es la matriz nula (p+1) x (p+1). Para encontrar la matriz inversa de Γu se aplica el algoritmo de Gauss generalizado. Las soluciones del sistema (6.1.13) son: El vector de pesos óptimos,
′
λ ′ = (γ + X ( X ′ Γ −1 X ) −1 ( x − X ′ Γ −1 γ )) Γ −1
(6.1.17)
y el vector de los multiplicadores de Lagrange: ′ m ′ = − x − X ′ Γ −1 γ ( X ′ Γ −1 X ) −1
(
)
(6.1.18)
Las justificaciones de las expresiones (6.1.17) y (6.1.18) se encuentran en el apéndice C6. La varianza del kriging se puede expresar como: σ 2k (s 0 ) = − λ*′ Γ −1 λ* + 2 λ*′ γ = λ*′ γ + ( λ*′ γ − λ*′ Γ −1λ* ) donde la expresión entre paréntesis luego de trabajo algebraico se transforma en m ′ x . Por lo tanto; σ 2k (s 0 ) = λ*′ γ + m ′ x = λ*u′ γ u (6.1.19) Otra expresión de la varianza del kriging es:
σ 2k (s 0 ) = γ ′ Γ −1 γ − ( x − X ′ Γ −1 γ ) ′ ( X ′ Γ −1 X ) −1 ( x − X ′ Γ −1 γ )
(6.1.20)
Los pesos de ponderación óptimos y la varianza del kriging obtenidos a través de las ecuaciones (6.1.17) y (6.1.20), ó, de las obtenidas con la función de covarianzas, permiten construir el intervalo de predicción nominal del 95% para Z (s 0 ) es: A ≡ ( Z$ (s ) − 196 . σ (s ) , Z$ (s ) + 196 . σ (s ) ) 0
k
0
0
k
0
Bajo el supuesto que Z(⋅) es Gaussiano, la Prob{ Z(s0) ∈ A }= 95%.
73
Observaciones: • Tendencia polinómica. Si s ∈ R2 a menudo E ( Z(s)) = µ (s) se expresa como una combinación lineal de polinomios en las coordenadas espaciales s = ( x , y ) ′ . Una superficie de tendencia de grado r es: µ (s) = ∑ ∑ a kl x k y l (6.2.7) 123 0 ≤ k +l ≤ r
Por ejemplo, una superficie de tendencia cuadrática es: µ (s) = a 00 + a10 x + a 01 y + a 20 x 2 + a11 x y + a 02 y 2 La expresión s = ( x , y )′ y
(6.2.7) es un caso particular de la tendencia en (6.1.1), con
(r + 1)(r + 2) −1 . 2 • Las ecuaciones de kriging universal pueden considerarse una generalización de las de kriging ordinario. Si se supone que Z es una variable intrínseca, con esperanza µ, entonces puede formularse como no intrínseca con p=1, β l = µ y f l (s) = 1. • Para predecir óptimamente un valor no conocido Z(s0), solo se necesita conocer: var Z (s i ) − Z (s j ) ; 0 ≤ i ≤ j ≤ n
f 0 ( s) = 1
f 1 ( s) = x
K
{ (
f p ( s) = y r
donde p =
}
)
No es necesario conocer los coeficientes β l de la combinación lineal desconocida de funciones conocidas
{f
0
}
( s) , L , f p ( s) .
• La descomposición (6.1.1) no puede ser obtenida fácilmente. Dentro de cada disciplina, los científicos tienen a menudo una buena idea acerca de que parte de Z es debida a factores controlables y a determinadas variables exógenas. Pero aún así no existe unanimidad, porque la descomposición Z(s) = µ (s) + δ (s) en definitiva depende de preferencias y gustos personales.
6.2 Estimación del variograma para el kriging Universal. En las ecuaciones del kriging universal se supone que el variograma es conocido. En la práctica, debe ser estimado, y la estimación del mismo no es sencilla. El uso de los estimadores citados en el capítulo 3 sería inapropiado porque: 2 2 E Z (si )− Z (s j ) =Var Z (si )− Z (s j ) + E Z (si )− Z (s j )
(
)
(
)[ (
)]
2
p +1 = 2γ (s i −s j ) + ∑ β k −1 ( f k −1 (si )− f k −1 (s j ) ) k =1 es decir, los estimadores propuestos en el capítulo 3 serían sesgados. Si β fuera conocido, un estimador para el variograma podría estar basado en p +1
δ (⋅) ≡ Z (⋅) − ∑ β k −1 f k −1 (⋅) porque E (δ (s i ) − δ (s j ) ) = 2 γ (s i − s j ) . 2
k =1
El vector de parámetros β es desconocido, pero se puede estimar fácilmente.
74
Como Z satisface el modelo lineal general, donde E (Z ) = Xβ y Var(Z ) = Σ , el estimador obtenido por mínimos cuadrados generalizados de β es: −1 β$ = X ′ Σ −1 X X ′ Σ −1Z GLS
(
)
(
)
que necesita el conocimiento de la matriz de varianzas y covarianzas Var (Z ) = Σ , o equivalentemente basta conocer el variograma. Pero, 2 γ (⋅) es desconocido, llevando la discusión a foja cero. Esta circularidad sumada a la última observación del apartado anterior provocan algunas disconformidades con el kriging universal. En forma general, considérese los residuos g.l.s. en el contexto general del modelo (6.1.1). Desde (6.1.2) Z = X β +δ el estimador g.l.s. de β es −1 β$ = X ′ Σ −1 X X ′ Σ −1Z (6.2.1) GLS
y los residuos correspondientes son W = Z − X β$ gls = Σ − X ′ Σ −1 X
(
(
)
−1
)
X ′ Σ −1 Z
(6.2.2)
Los estimadores del variograma basados en W son sesgados, además de la dependencia estadística (expresada en la matriz de varianzas Σ ), W satisface (p + 1)
(
(
restricciones lineales, la matriz de proyección Σ − X ′ Σ −1 X
)
−1
)
X ′ Σ −1 en (6.2.2) es de
rango n− p− 1. Intuitivamente, tales restricciones algebraicas inducen a los residuos W que exhiban más correlaciones negativas que aquellas de los errores δ . Argumentos similares se aplican a los residuos o.l.s. Un gran número de autores propuso soluciones para los problemas de sesgo que se presentan cuando se usan los residuos para estimar el variograma Una de estas propuesta es la de Neuman y Jacobson. El enfoque de Neuman y Jacobson es un proceso iterativo cuyas etapas son: Etapa 1: obtener una estimación de β por mínimos cuadrados ordinarios, es decir ajustar µ(s) por mínimos cuadrados ordinarios. Etapa 2: realizar una estimación del variograma desde los residuos y ajustar un modelo de variograma. Etapa 3: obtener una estimación de β por mínimos cuadrados generalizados basándose en el modelo de variograma ajustado en la etapa anterior. El proceso interactivo consiste en repetir las dos últimas etapas hasta lograr convergencia. Pero resultados empíricos como los presentados en Cressie (1990) muestra que este método no soluciona el problema del sesgo. Por otro lado, se sabe que para un estimador basado en n observaciones ocurre que ( sesgo) 2 = (O( 1 n )) y var ianza = O( 1 n ) . Es decir, cuando n tiende a infinito el 2
( sesgo) 2 tiende a cero más rápido que la varianza y se considera que el sesgo puede ser ignorado. Sin embargo, en muchos problemas de Geoestadística n no es usualmente grande y por lo tanto no es adecuado el resguardo asintótico. Es generalmente verdadero que el sesgo de un estimador de variograma basado en los residuos es pequeño en retardos cercanos al origen pero más sustancial en retardos distantes. Cuando un modelo de variograma es ajustado por mínimos cuadrados
75
generalizados o por mínimos cuadrados ponderados automáticamente pone más peso al estimador en los retardos pequeños, y por lo tanto el efecto del sesgo sería pequeño. Además, si el kriging es llevado a cabo en un entorno local, el variograma ajustado es solamente evaluado en retardos pequeños, precisamente donde el variograma ha sido bien ajustado. Esto es, en definitiva, el sesgo en el variograma no influiría en gran medida en los valores de las predicciones. La estimación de la varianza del kriging es la más propicia para ser afectada por el sesgo en el estimador del variograma. Por simplicidad, supóngase que el variograma tenga una meseta, esto es lim 2 γ (h) = 2 σ 2 . h →∞
W′ W . Según los resultados de n − p −1 Cressie (1990), bajo el supuesto que el covariograma es positivo, el sesgo de σ$ 2 es 2 O(1 n) y negativo, y como la varianza del kriging σ k es directamente proporcional a Un estimador de σ 2 comúnmente usado es σ$ 2 =
σ 2 , se obtiene usualmente una estimación sesgada de σ k 2 . Los resultados de Cressie y Grondona (1992) reafirman que el sesgo de la estimación del variograma es importante para n pequeño y h (retardo) grande. Además concluyen que la presencia de sesgo es causada por los contrastes lineales en los residuos R (o W). En conclusión, aunque el predictor de kriging universal puede ser poco influenciando por el sesgo, existe evidencia tanto experimental como teórica que la varianza del kriging estimada puede ser más pequeña que lo que debería ser.
6.3 Kriging mediana polish. Los datos espaciales se pueden pensar como un muestreo parcial de una realización de un proceso aleatorio {Z(s): s ∈ D}. Además el proceso Z se lo modela de acuerdo a la descomposición (6.1.1) Z(s) = µ (s) + δ (s) (6.3.1) donde µ (⋅) ≡ E ( Z (⋅)) es la estructura de medias y δ (⋅) es la estructura del error. En la realidad µ (⋅) no es conocida y una de las formas de modelarla ya fue presentada en los apartados anteriores. En los espacios de dimensión dos o más, es natural suponer que µ (⋅) se descompone aditivamente en componentes direccionales. Por ejemplo en R2,
µ ( s) = a + c ( x ) + r ( y )
s ∈D ⊂ R 2 (6.3.2) Además, si { s i : i =1,2,K , n} están ubicados sobre una grilla de p filas y q columnas:
{( x l , y k ) ′: k = 1,K , p; l = 1,K , q}
Entonces, por simplicidad (6.3.2) se puede expresar como:
µ (s i ) = a + rk + cl donde s i = ( x l , y k ) ′
(6.3.3) El efecto fila rk puede ser estimado mediante la replicación en la otra dimensión; es decir, rk puede ser estimado desde:
{Z (s i ): la segunda coordenada de s i sea
y k ; i = 1,K , n} donde k = 1,K , p . En forma similar se puede estimar el efecto columna cl , l = 1,K , q .
76
Datos grillados. Los datos espaciales grillados en R2 pueden ser considerados como una tabla de doble entrada (o de más vías en Rd). Es importante tener en cuenta que los espaciados en ambas direcciones, horizontal y vertical, no tienen porque ser iguales. Miller y Kahn (1962) propusieron un análisis de la varianza formal de dos vías para testar la no estacionariedad a través de las filas y a través de las columnas mediante el uso del estadístico F. Desafortunadamente, los tests basados en valores críticos obtenidos de una tabla F son incorrectos porque los datos son correlacionados. Datos no grillados. Para tratar con datos no grillados se realiza el trazado de un mapa de baja resolución de las ubicaciones espaciales. La resolución de las coordenadas espaciales es elegida a menudo en una manera ad hoc de modo que cada combinación ( x l , y k ) tenga aproximadamente una observación Z ( x l , y k ) en ( x l , y k ) . En la práctica, esto se hace por la superposición de una grilla sobre el mapa de alta resolución y asignando a las ubicaciones {s i : i = 1,K , n} de los datos a los nodos
más cercanos de la grilla {( x l , y k ) ′ : k = 1,K , p ; l = 1,K q} . De acuerdo a esta resolución, Z (s i ) es escrito como Z ( x l , y k ) . Entonces, se aplica el método de la
mediana polish al conjunto de datos {Z ( x l , y k )} , pero antes se debe elegir la resolución y la orientación de la grilla.
Análisis por las medias. Bajo el supuesto de una observación en cada nodo de una grilla, los estimadores mínimos cuadrados ordinarios de las componentes aditivas de µ (⋅) son: n
a$ =
∑ Z (s ) i
i =1
,
n
∑ Z (s ) i
r$k =
N ( yk )
c$l =
M ( xl )
q
− a$
k = 1,K , p
− a$
l = 1,K , q
∑ Z (s ) i
p
N ( y k ) ≡ {i : s i = (⋅ , y k ) ′; i = 1,K , n }
(6.3.4)
M ( x l ) ≡ {i : s i = ( x l , ⋅ ) ′; i = 1,K , n }
Y el estimador mínimo cuadrados ordinarios de µ (s i ) es: µ$ (s i ) = a$ + r$k + c$l donde s i = ( x l , y k ) ′ (6.3.5) Para s = ( x , y ) ′ ubicados en la región acotada por las líneas que unen los cuatros nodos, ( x l , y k ) ′ , ( x l +1 , y k ) ′, ( x l , y k +1 ) ′, y ( x l +1 , y k +1 ) ′ , donde x l < x l +1 y y k < y k +1 se define el interpolador plano (6.3.6): y − yk x − xl µ$ (s) ≡ a$ + r$k + (r$k +1 − r$k ) + c$l + (c$l +1 − c$l ) k = 1,K, p − 1; l = 1,K, q − 1 yk +1 − yk xl +1 − xl Este estimador es no paramétrico, porque por ejemplo no se ha impuesto forma de tendencia lineal o cuadrática, y presenta una desventaja que consiste en que los
77
residuos
{Z (s i ) − µ$ (si ): i = 1,K, n} permiten
estimadores sesgados de la dependencia
espacial desconocida del proceso del error δ (⋅) . El enfoque tomado en esta sección consiste de dos etapas. Primero, la estructura de la media es estimada y removida. Entonces, la estructura de la dependencia espacial es estimada. El procedimiento podría repetirse, en donde una estimación más eficiente de la estructura de la media es obtenida la segunda vez, basada en la estimación de la dependencia espacial en la primera vez. En efecto, las dos etapas podrían ser iteradas hasta que existan pocos cambios en los resultados. A los efectos de evitar el problema del sesgo en los estimadores de la dependencia espacial, se propusieron estimadores no lineales de a, { rk } y { cl }. Cuando la distribución del proceso del error es simétrica, entonces E ( promedio{Z (s i ): i ∈ A}) = E (mediana{Z (s i ): i ∈ A}) . Además la mediana de un conjunto de datos tiene la propiedad de ser resistente a los outliers. Como los datos generalmente son transformados para que su distribución sea aproximadamente simétrica, (Cressie, 1985), las medianas análogas a las medias aritméticas de (6.3.4) pueden permitir residuos menos sesgados, en el sentido que los estimadores de la dependencia espacial desconocida del proceso del error sean menos sesgados. Este algoritmo es denominado de la mediana polish. 6.3.1 Mediana polish. A través de un ejemplo se presentará el algoritmo de la mediana polish. Bajo el supuesto natural que µ (⋅) se descompone aditivamente en componentes direccionales, este algoritmo produce el efecto mediano total a~ , los efectos medianos fila { r~k : k = 1,K , p} y los efectos medianos columnas { c~l : l = 1,K , q} desde un arreglo de números de tamaño p x q {Ykl : k = 1,K , p; l = 1,K , q } . En el contexto espacial, los datos grillados { Z (s i ): i = 1,2,K , n} juegan el papel de las Y. Ejemplo: Se desea aplicar el algoritmo de la mediana polish a los datos grillados que se presentan en la siguiente tabla: 8 9 11 13 7
3 11 8 4 2
6 12 4 7 3
6 7 2 6 5
5 3 3 1 0
Etapa inicial o etapa 0. Para iniciar el algoritmo, a las 25 celdas se le agregarán 11 celdas con ceros, es decir se añadirá una fila y una columna de ceros. 8 9 11 13 7 0
3 11 8 4 2 0
6 12 4 7 3 0
78
6 7 2 6 5 0
5 3 3 1 0 0
0 0 0 0 0 0
A los efectos de la explicación se agregará otra fila y otra columna. En la última columna se cargan las medianas de cada fila. 8 9 11 13 7 0
3 11 8 4 2 0
6 12 4 7 3 0
6 7 2 6 5 0
5 3 3 1 0 0
0 0 0 0 0 0
6 9 4 6 3 0
Etapa 1: a cada fila se le resta la mediana de la de la etapa anterior a los datos y se acumulan las medianas en la penúltima columna. 2 0 7 7 4 0
-3 2 4 -2 -1 0
0 3 0 1 0 0
0 -2 -2 0 2 0
-1 -6 -1 -5 -3 0
6 9 4 6 3 0
Como en la próxima etapa se necesitan las medianas por las columnas, se las agrega en la última fila. 2 0 7 7 4 0 4
-3 2 4 -2 -1 0 -1
0 3 0 1 0 0 0
0 -2 -2 0 2 0 0
-1 -6 -1 -5 -3 0 -1
6 9 4 6 3 0 6
Obsérvese que el número 6 de la última celda de la última fila no es más que la mediana de las medianas por fila. Etapa 2: se remueve las medianas por columnas de la etapa anterior desde los datos y se acumulan en la penúltima fila. -2 -2 0 0 0 0 -4 3 3 -2 -5 3 3 5 0 -2 0 -2 3 -1 1 0 -4 0 0 0 0 2 -2 -3 4 -1 0 0 -1 6
Como en la próxima etapa se necesitan las medianas de los datos por las filas, se las agrega en la última columna.
79
-2 -4 3 3 0 4
-2 3 5 -1 0 -1
0 3 0 1 0 0
0 -2 -2 0 2 0
0 -5 0 -4 -2 -1
0 3 -2 0 -3 6
0 -2 0 0 0 0
Etapa 3: se procede al igual que en la etapa 1, es decir se remueve las medianas por filas de la etapa anterior desde los datos y se acumulan en la penúltima columna. -2 -2 3 3 0 4
-2 5 5 -1 0 -1
0 5 0 1 0 0
0 0 -2 0 2 0
0 -3 0 -4 -2 -1
0 1 -2 0 -3 6
Como en la próxima etapa se necesitan las medianas por las columnas, se las agrega en la última fila. -2 -2 0 0 0 0 -2 5 5 0 -3 1 3 5 0 -2 0 -2 3 -1 1 0 -4 0 0 0 0 2 -2 -3 4 -1 0 0 -1 6 0 0 0 0 -2 0 Etapa 4: se procede al igual que en la etapa 2, obteniéndose la siguiente tabla: -2 -2 3 3 0 4
-2 5 5 -1 0 -1
0 5 0 1 0 0
0 0 -2 0 2 0
2 -1 2 -2 0 -3
0 1 -2 0 -3 6
0 0 -4 0 2 0 0
2 -1 0 -2 0 -3 0
0 1 0 0 -3 6 0
Etapa 5: se procede al igual que en la etapa 3. -2 -2 1 3 0 4 0
-2 5 3 -1 0 -1 0
0 5 -2 1 0 0 0
80
0 0 2 0 0 0
En este ejemplo el método iterativo converge, puesto que la próxima tabla no se modificará ya que se deberán desagregar ceros y agregar ceros según corresponda. La tabla original 5x5 es reemplazada por una tabla 5x5 de residuos, más el efecto mediano total que se encuentra en la celda (6,6), más los efectos filas medianos ubicados en la columna 6 y los efectos columnas medianos que se encuentran ubicados en la fila 6. En general, el algoritmo Mediana polish es el siguiente: Para i =1,3,5, ..., se define
{ + med {Y
} : l = 1,K , q}
Ykl( i ) ≡ Ykl(i −1) − med Ykl( i −1) : l = 1,K , q Yk(,iq) +1 ≡ Yk(,iq−+11)
( i −1) kl
k = 1,K , p + 1 l = 1,K , q k = 1,K , p + 1
(6.3.7)
Para i =2,4,6, ..., se define
{ + med {Y
} : k = 1,K , p}
Ykl(i ) ≡ Ykl(i −1) − med Ykl( i −1) : k = 1,K , p Yp(+i )1, l ≡ Yp(+i −11, l)
( i −1) kl
k = 1,K , p l = 1,K , q + 1 l = 1,K , q + 1
(6.3.8)
donde med { y1 , y 2 ,K , y n } es la mediana de { y1 , y 2 ,K , y n } . Para iniciar el algoritmo se supone que Ykl k = 1,K , p l = 1,K , q ( 0) Ykl = 0 en otro caso Es decir se inicia con las pxq celdas correspondientes a los datos a los que se les agregan p+q+1 celdas con ceros. Se usa (6.3.6) para remover las medianas por filas desde los datos y acumular las cantidades removidas en las p celdas filas extras (que forman la columna q + 1). De la misma manera en las columnas de la tabla, removiendo las medianas por columnas desde no solo los datos sino también de la columna de las remociones acumuladas de las filas. Esta última cantidad removida es la entrada de la celda extra (p+1, q+1). Este proceso se debe repetir hasta lograr convergencia. Suponiéndose convergencia, los efectos estimados son a~ ≡ Yp(+∞1), q +1 r~ ≡ Y ( ∞ )
k = 1,K , p
c~l ≡ Yp(+∞1), l
l = 1,K , q
k
k , q +1
(6.3.9)
con la propiedad que Ykl = a~ + r~k + c~l + Ykl( ∞ )
(6.3.10)
Esto es, la tabla original p x q es reemplazada por la tabla p x q de residuos {Y : k = 1,K , p ; l = 1,K q } , y las p+q+1 celdas extras contienen los efectos medianos por filas { r~k : k = 1,K , p} , los efectos medianos por columnas { c~l : l = 1,K , q} y el efecto mediano total. (∞) kl
81
6.3.2 Superficie mediana polish. En un contexto espacial Z se puede expresar como Z (s i ) = a~ + r~k + c~l + R (s i )
si = ( xl , y k ) ′
(6.3.11)
~ (s ) ≡ a~ + r~ + c~ es un estimador flexible de la superficie media que es Donde µ i k l resistente a los outliers. Para los pares ordenados que representan posición: s = ( x , y ) ′ ubicados en la región acotada por las líneas que unen los cuatros nodos, ( x l , y k ) ′ , ( x l +1 , y k ) ′, ( x l , y k +1 ) ′, y ( x l +1 , y k +1 ) ′ , donde x l < x l +1 y y k < y k +1 se define el interpolador plano (6.3.12):
y− y k y k +1 − y k
µ~(s)≡a~ +r~k +
~ ~ ~ x − xl (rk +1 −rk )+cl + xl +1 − xl
~ ~ (cl +1 −cl )
k = 1,K, p − 1;
l = 1,K, q − 1 También es posible extrapolar más allá de la grilla de observaciones. Supóngase que x < x1 pero y1 ≤ y k ≤ y ≤ y k +1 ≤ y p entonces, para s = ( x , y ) ′ se define y − yk ~ x − x1 ~ ~ ( rk +1 − r~k ) + c~1 + ( c2 − c1 ) y k +1 − y k x 2 − x1
~ (s) ≡ a~ + r~ + µ k
k = 1,K , p − 1
Una fórmula similar se obtiene cuando y está fuera de rango o cuando ambas x y y lo están. Por lo tanto, a través de la interpolación y la extrapolación se define la ~ (s)): s ∈R 2 para todo el plano. superficie mediana polish (s , µ
{
}
Observaciones. • Si en cada nodo de la grilla el número de observaciones es distinto, la única modificación necesaria es notacional. El algoritmo no se altera en su esencia, sucesivamente remueve las medianas de las filas y las medianas de las columnas desde las entradas de la tabla. Si una fila (columna) entera de nodos de la grilla no tiene observaciones, entonces la fila (columna) es ignorada. • El algoritmo se inicia con la remoción por las filas, pero sin problema alguno el proceso podría iniciarse por las columnas. • En la práctica, se necesita un criterio para detener el algoritmo; por ejemplo, cuando otra iteración deja inalterable cada entrada de la tabla dentro de una tolerancia ε preestablecida, la mediana polish termina y el algoritmo se dice que converge. 6.3.3 Kriging basado en los residuos de la mediana polish. Los residuos { R (s i ): i =1,K , n } obtenidos mediante mediana polish se piensan como un nuevo conjunto de datos espaciales, a los cuales se les aplica kriging ordinario. Modelo supuesto • R(s) ≈ δ(s) s∈ D. • El proceso δ(s) es tal que E(δ(s)) = 0 para todo s. • El proceso R(.) tiene variograma 2γ(h) = Var (R( s + h) − R(s ))
82
Predictor supuesto R$ (s 0 ) =
n
∑λ
n
i
∑λ
R( s i ) con
i =1
i
=1
i =1
Esta última condición sobre los coeficientes del predictor lineal garantiza insesgadez uniforme. El vector de coeficientes λ*nx1 ≡ ( λ*1 ,K , λ*n ) ′ esta dado por (5.3.9). El kriging basado en los residuos se comporta como una proxi para el kriging basado en los errores desconocidos {δ (s i ): i = 1 ,K , n } . ~ (s ) puede ser extendida para s ∈R 2 y el Como la estimación mediana polish µ 0
0
kriging basado en los residuos permite R$ (s 0 ) para todo s 0 ∈R 2 . De esta manera, el predictor kriging mediana polish de Z(s0) es definido como ~ ~ (s ) + R$ (s ) Z (s 0 ) ≡ µ s0 ∈ R 2 (6.3.13) 0 0 ~ ~ Nótese que Z (⋅) es un interpolador exacto; esto es, Z (s i ) = Z (s i ) i = 1,K , n . ~ La varianza del kriging asociado con el predictor kriging mediana polish Z (s 0 ) es definida como la varianza del kriging ordinario basados en los residuos obtenidos con la mediana polish. Se la simboliza como σ 2m (s 0 ) y su expresión es n
σ 2m (s 0 ) ≈ ∑ λ i γ (s 0 − s i ) + m
(6.3.14)
i =1
El procedimiento para el caso de que los datos no estén grillados es el siguiente: Primero se define la grilla de baja resolución, luego se aplica el método de la mediana polish para extraer la tendencia, y por último los residuos son pensados como un nuevo conjunto para realizar el kriging. Para el caso grillado vale la propiedad de interpolación exacta, Z~ (s i )=Z (s i ), i = 1K , , n . En el caso no grillado, la preservación de esta propiedad es considerada importante.
83
CAPITULO 7: APLICACIÓN A LA HIDROGEOLOGÍA, "ACUIFERO DE TUCSON” En este capítulo y en el siguiente se presentarán problemas que involucran datos espaciales provenientes del campo de investigación de la Hidrogeología, cuyo tratamiento requiere de las herramientas que han sido presentadas a lo largo de los capítulos anteriores. Se trabajará con datos correspondientes a variables geoquímicas de aguas subterráneas. En esta primera parte se trabajó con un conjunto de datos obtenidos en Internet de la página http:/www.u.arizona.edu/. De acuerdo a la información suministrada en dicho sitio, provienen de la tesis de Hidrología UA MS 19886, con la cual se pretendió caracterizar el acuífero de Tucson. En el capítulo 8 se trabajará con datos del sistema acuífero de La Caldera. En una primera sección se analizan los datos en forma descriptiva. En la segunda, se realiza el análisis estructural, como en todo estudio geoestadístico esta etapa juega un papel preponderante. En la tercera sección se realiza el kriging. Como en la mayoría de los trabajos con variables hidrogeológicas se transforman las variables, más precisamente se trabaja con el logaritmo de las mismas, en la cuarta sección se presenta un estudio considerando la transformación de la variable de interés. En la última sección se comparan los resultados obtenidos.
7.1 Análisis descriptivo de los datos. Se dispone de valores de las variables; concentración de: Calcio, Magnesio, Sodio, Potasio, HCO3, SO4, Cloro y NO4. Como así también los valores de las coordenadas Este y Norte de 114 lugares donde se midieron las variables de interés. La Figura 7.1 muestra las ubicaciones donde fueron tomadas las 114 observaciones. La distribución de las coordenadas define aproximadamente un rectángulo de 10.8 pies en la dirección oeste - este por 5.7 pies en la dirección sur – norte. Las ubicaciones de las observaciones están irregularmente espaciadas, aunque existen algunos claros y grupos, destacándose los claros en la zona del vértice superior izquierdo, en la zona del vértice inferior derecho y en la zona del vértice superior derecho.
Figura 7.1: Posiciones de los datos correspondientes al estudio del acuífero de Tucson.
84
De todas las variables químicas disponibles se analizará el comportamiento de la variable concentración de Calcio en la región de interés. Las medidas de la misma en p.p.m. se presentan con sus respectivas posiciones en la Figura 7.2.
6
1.85
1.22 1.45 2.87
0.93 0.79 0.95 2.15 0.8 1.84 0.55 0.3 7.53
5 2.5511.46 1.98
1.32
2.5 1.8 2
4
3
2
1
1.5 0.51 0.45 0.340.48 0.640.76 0.34
3.19
1.91 2.54
2.2 1.72 1.5 1.65 1.68 1.08 0.54 3.21 1.81 1.77 2.21 1.96 3.34 3.19 1.23 3.69 2.6 3.74 2.99 2.99 1.31 1.9 0.85 0.94 2.253.3 0.851.21 3.99 1.1 2.2 2.59 0.88 1.36 0.55 1.83 1.39 2.4 2.2 0.650.51 0.870.58 1.7 1.39 0.79 0.8 1.55 1.43 1.72 1.5 2.54 1.15 1.12 2.02 2.291.34 2.04 1.08 1.9 0.94 2.2 0.79 2.45 1.85 1.3 0.88 1.07 1.12 1.25 3.08 1.11 0.6 1.47 2.28 2.42 1.19
2
4
6
8
0.63
0.55 3.09 1.59
2.15 0.74
10
Figura 7.2: Medidas de calcio en sus respectivas posiciones correspondientes al estudio del acuífero de Tucson.
En el mapa de los datos se observa en rojo los puntos de coordenadas aproximadas (4.5, 5.5) y (5.0, 4.8) donde la variable presenta valores alejados con respecto a los otros. Esto queda evidenciado a través del scatter plot tridimensional, como el presentado en la Figura 7.3. 3D Scatterplot (Tucwate.STA 10v*114c)
Figura 7.3: Scatter -plot para las medidas de la variable calcio.
85
En la Figura 7.4 se observa que la distribución de la concentración de Calcio observada es asimétrica y presenta valores muy alejados 7.53 y 11.46, los cuales coinciden con los anteriormente observados.
0
5
10
15
CALCIO Figura7.4: Box-Plot para el contenido de Calcio.
En la figura 7.5 se representa el post plot para la cantidad de calcio en cada observación. Con el triángulo y el color celeste se representan los valores correspondientes a los valores menores o iguales al primer cuartil, o sea a 0.93. Con la estrella y el color verde se representan los valores comprendidos entre 0.93 y el valor mediano 1.5. Con el diamante y el color azul se representan los valores comprendidos entre 1.5 y el tercer cuartil 2.2 y con cuadrados rojos a los valores mayores que el tercer cuartil.
norte 6 5 4 3 2 1 este 2
4
6
8
10
Figura 7.5: Post Plot para la variable calcio.
Si bien todos los valores están dispersos en la región, los menores se presentan casi todos por debajo de la coordenada 3 del Norte. Los valores de cada categoría determinada por los cuartiles tienden a alinearse de alguna manera. La mayoría de los valores más grandes se encuentran entre las coordenadas 4 y 7 del Este.
86
7.1.1 Otras herramientas exploratorias. Los scatter plot de la variable de interés versus cada una de las variables que definen la posición de los puntos, mostrados en la Figura 7.6, señalan que no se presenta tendencia en los valores de la variable calcio ni en la dirección Oeste-Este ni en la dirección Sur-Norte.
15
10
10
CALCIO
CALCIO
15
5
5
0 0
2
4
6 ESTE
8
10
0 0
12
1
2
a)
3 4 NORTE
5
6
7
b)
Figura 7.6: Diagrama de puntos de: a) calcio vs. Coordenada este. b) calcio vs. Coordenada norte.
A los efectos de detectar posibles tendencias y calcular el estadístico u se procede a calcular medias y medianas por filas y columnas. Con franjas de amplitud 0.2 se definieron las filas en la dirección Sur - Norte y las columnas en la dirección Oeste Este. La Figura 7.7 es un intento de resumir la posible no estacionariedad en la media a lo largo de filas y columnas, es decir en las direcciones oeste- este y sur- norte, usando la media muestral y la mediana muestral a lo largo de las filas y las columnas respectivamente.
a)
b)
Figura 7.7: En ambas figuras con azul se representan los valores de la media y con rojo los valores de la mediana de los valores de contenido de calcio. a) Según columnas. b) Según filas.
Cuando sólo se presenta el color rojo por banda es porque en aquellos puntos se superponen media y mediana. Además, se observa que no existe ninguna tendencia en los valores medios y medianos a través de las columnas ni de las filas.
87
En las tablas 7.1 y 7.2 se presentan las diferencias estandarizadas, para los valores que estamos tratando, es decir los valores del estadístico u presentados en la sección 2.2.3.
Tabla 7.1: Valores del estadístico u según columnas. Tabla 7.2: Valores del estadístico u según filas.
88
En la tabla 7.1 el valor del estadístico u señala que en las columnas correspondientes a los valores medios 4.894 y 6.14 de la variable este, la variable en estudio podría presentar outliers. Así para la primer columna sospechosa, el valor es Calcio(4.86, 4.77)= 11.46 p.p.m.(ver Figura 7.5), es decir en esta franja se ubica uno de los outliers ya detectado anteriormente. A la lista de valores alejados se le agrega Calcio(6.23, 2.83)= 3.69 p.p.m. aportado por la franja cuyo valor medio de “este” es 6.14. Este es un valor alejado con respecto a los valores de la franja cuyas ubicaciones tienen las siguientes características: estemedia
6.14
este 6.05 6.1 6.11 6.21 6.23
norte 0.8 1.67 2.34 2.11 2.83
calcio 1.21 1.15 1.1 0.65 3.69
El efecto del otro valor alejado detectado anteriormente sobre el valor del estadístico u quedó neutralizado por la composición de la franja correspondiente a un valor medio de “este” igual a 4.325. En la tabla 7.2 se presentan los valores de u correspondientes a las filas. El valor de u=2.79486 es producto de la franja cuya coordenada norte media es 1.05, ésta recibe el aporte de 5 valores. El nuevo valor alejado espacial es Calcio(3.95, 1.05)= 3.08 que es más grande que el tercer cuartil: 2.2, mientras los otros valores son menores que la mediana: 1.5 (las medidas de posición son las correspondientes al conjunto total de datos). En la franja cuya coordenada norte media es 4.5175 el valor Calcio(4.85, 4.59)= 1.98 es más grande que la mediana y los tres valores restantes son inferiores al primer cuartil: 0.93. El valor alejado Calcio(4.86, 4.77)= 11.46 es detectado por la franja cuya coordenada norte media es 4.69. En cambio Calcio(4.29, 5.4)= 7.53 es neutralizado por los otros dos valores que determinan la franja con coordenada norte media: 5.49667.
7.2 Análisis Estructural. 7.2.1 Estimación del variograma. En este apartado, se intenta sobre la base de los datos disponibles determinar la estructura de correlación espacial del proceso aleatorio Z(x): “contenido de calcio”. Bajo el supuesto que el proceso Z(x) es intrínsecamente aleatorio, la función semivariograma γ(h) cuantifica dicha estructura. A partir de una realización se debe estimar el semivariograma. Con el programa Vario que forma parte del Geoeas se genera el semivariograma experimental “omnidireccional”, el cual permite utilizar todos los pares de datos independientemente de la dirección. Se considera la dirección de 0° y una tolerancia angular de 90° y para completar la definición de los intervalos planos las siguientes condiciones para las distancias: Mínimo =0, Máximo = 11.156. Es decir se toma en cuenta la máxima distancia entre los puntos, y el incremento se fija en 0.561, de esta manera se generan 20 intervalos de distancia. En la figura 7.8 se presenta el semivariograma experimental “omnidireccional” correspondiente. 89
Figura 7.8: Semivariograma experimental “omnidireccional”.
En la figura 7.9 se presenta el mismo semivariograma, en a) se conectaron los puntos del semivariograma a los efectos de la visualización, y en b) se indican la cantidad de pares de valores que intervinieron para el cálculo del semivariograma correspondiente al valor medio de los valores de separación que definen el intervalo plano. Así se observa que a partir del retardo 14 la cantidad de pares disminuye en forma rápida pasando desde 64 hasta llegar al valor 4 que es la cantidad de pares que aportan al valor del semivariograma 1.694 correspondiente a una distancia promedio de 10.945.
Figura 7.9: Semivariograma experimental “omnidireccional”.
Para observar que las elecciones de los intervalos de distancias no influyen en la “forma” del semivariograma muestral se realizan sucesivos gráficos con distintas longitudes de los intervalos de distancia como el presentado en la figura 7.10. Resultados comparables se obtienen si se trabaja con el programa Variowin, para la obtención del semivariograma empírico. Dichos resultados no son iguales debido a la manera diferente de considerar la construcción de los intervalos de distancia. La figura 7.11 muestra un semivariograma empírico generado por este programa.
90
Figura 7.10: Semivariograma experimental “omnidireccional”, con las siguientes condiciones para las distancias: Mínimo =0, Máximo= 11.156, y el incremento se fija en 1.2.
Figura 7.11: Semivariograma experimental “omnidireccional” generado con Variowin, con las siguientes condiciones para las distancias:10 intervalos de amplitud 1.12 más el intervalo coorespondiente al retardo 0 de amplitud de 0.56.
En esta última figura se observa que la cantidad de pares de valores que intervinieron para el cálculo del semivariograma en los distintos retardos es aproximadamente el doble que el que utiliza el otro programa. Este es por la manera de la elección de los intervalos de retardo, y porque Variowin usa los datos en forma doble para los cálculos. Una regla empírica afirma que los variogramas no son generalmente válidos más allá de la mitad de la máxima distancia entre muestras (en esta situación 11.156). Esto y lo observado en la figura 7.9 b); en cuanto a la confianza de los valores dado por la cantidad de pares que intervienen en el cálculo del semivariograma empírico. Se considera la dirección de 0° y una tolerancia angular de 90° y para completar la definición de los intervalos planos las siguientes condiciones para las distancias: 91
Mínimo = 0, Máximo = 7.4, y el incremento se fija en 0.561. Así se tiene 13 clases de distancias iguales.
Figura 7.12: Ssemivariograma experimental “omnidireccional”, con las siguientes condiciones para las distancias: Mínimo =0, Máximo= 7, y el incremento se fija en 0.561.
A los efectos de comparar más adelante los ajustes de los semivariograma empíricos a semivarigramas teóricos realizados con ambos programas, se considera la generación de los intervalos de distancia de la misma manera en ambos programas.
Figura 7.13: Semivariogramas experimentales, los de la derecha fueron generados por el Geoeas y los de la izquierda por el Variowin.
En la Figura 7.13, en la primera fila, los intervalos de distancia se generan con centro en 1.12 k con k=0,1,2,3,4,5,6 y radio 0.56. En la segunda fila, intervalos con centro en 0.86 k con k=0,1,2,3,4,5,6 ,7,8 y radio 0.43. 7.2.2 Ajuste a modelos de semivariogramas. La idea es buscar un semivariograma válido que represente más ajustadamente a la dependencia espacial presente en los datos: Z = (Z(s1), ..., Z(sn) )’. El espacio de 92
todos los semivariogramas válidos es un gran conjunto, usualmente se elige una familia paramétrica de semivariogramas. Para la variable Calcio de acuerdo a las publicaciones y a la forma de los semivariogramas empíricos se elige la familia de semivariogramas esféricos. Modelando con el Geoeas y el Variowin. Con el programa Geoeas se realiza un ajuste a sentimiento. El cual no garantiza un modelo de variograma único ya que se basa en apreciaciones subjetivas y en la experiencia del usuario. Luego a través de la validación cruzada se trata de determinar la calidad y el grado de fiabilidad del modelo ajustado. No se debe intentar ajustar los mínimos detalles ya que en general éstos no son una característica del verdadero semivariograma sino más bien fluctuaciones muestrales. El valor del efecto pepita puede ser obtenido extrapolando los primeros puntos del semivariograma muestral hasta cortar el eje de ordenadas. Modelo 1. Por ensayo de prueba y error se ajusta el modelo esférico: 0 si h= 0 3 3 h 1 h γ(h;θ θ)= c0 + ce − si 0 < h ≤ a e Donde el vector de parámetros: 2 a 2 ae e θˆ = (cˆ0 , cˆ e , aˆ e )' = (0.9, 0.9, 2 )' c + c si h ≥ a 0 e e
Figura 7.14: Ajuste del semivariograma esférico (línea sólida) al semivariograma omnidireccional experimental.
En la figura 7.14 se muestra dicho ajuste. Es decir se ajusta un modelo esférico isotrópico con un efecto pepita igual a 0.9, y una meseta de 1.8 para un alcance de 2. A los efectos de determinar la calidad y el grado de fiabilidad del modelo ajustado se realizará la validación cruzada (ver capítulo 4). Así se considera un nuevo conjunto de datos: las diferencias entre las predicciones y los valores actuales correspondientes; es decir los residuos de la predicción.
93
La Figura 7.15 muestra el mapa de los residuos de la predicción cuando se realiza el ajuste antes mencionado. Se observa en la misma que se producen grandes residuos (los símbolos son proporcionales a los residuos) cuando se predicen los valores extremos 7.53 y 11.46. Otros residuos grandes pero de menor magnitud se presentan cuando la predicción (kriging) es realizada en zonas donde la información es escasa.
Figura 7.15: Mapa de los residuos de predicción cuando se ajusta el modelo 1.
Scatterplot with Box Plot (ajuste12345.STA 23v*114c)
Scatterplot with Box Plot (ajuste12345.STA 23v*114c)
6
7
4
6
2 0 RESIDUO1
CALCIOE1
5 4 3 2
-2 -4 -6 -8
1
-10 -12
0 0
2
4
6
8
10
12
0
14
1
2
3
4
5
6
CALCIOE1
CALCIO
Figura 7. 16: a) Scatter plot de los valores de la variable calcio vs. los valores predichos. b) Scatter plot de los valores de los residuos de predicción vs. los valores predichos.
En la Figura 7.16 a) se observa la influencia de los valores alejados de la variable Calcio, los cuales son predichos por valores pequeños en comparación a sus valores actuales. Pero hay otros valores que son predichos por valores más grandes que los valores actuales. Se esperaría que los puntos se distribuyan alrededor de la recta y = x , y estén muy próximos a ella.
94
7
En la Figura 7.16 b) se observa una leve tendencia lineal entre los residuos de predicción y los valores predichos, dicha estructura es consecuencia de la existencia de valores alejados. En la validación cruzada la herramienta de diagnóstico usual del ajuste es el estudio de la distribución de los residuos de predicción estandarizados. La Figura 7.17 muestra la distribución de los residuos (RESIDUO1) y los residuos de predicción estandarizados (RESEST1). Box Plot (ajuste12345.STA 23v*114c)
RESEST1
RESIDUO1
-12
-10
-8
-6
-4
-2
0
2
4
6
Non-Outlier Max Non-Outlier Min Median; 75% 25% Outliers Extremes Outliers Extremes
Figura 7.17: Box plot de los residuos de predicción y de los residuos de predicción estandarizados para el ajuste1.
La media de los residuos de predicción estandarizados es 0.012 con un desvío estándar de 1.15 y la mediana es de 0.2. Por otro lado ajustando el modelo 1 propuesto con el programa Model de Variowin da un índice de ajuste del 0.13249. Este índice de ajuste es una propuesta de Pannatier (1996) para medir el grado de ajuste a través de una cantidad en vez del ajuste a sentimiento. Modelo 2. Con el programa Model de Variowin se ajusta el semivariograma empírico a un modelo esférico que produce el mejor valor del índice de ajuste propuesto en dicho programa. El valor del indicador es 0.058295. En la figura 7.18 se muestra dicho ajuste. El modelo es esférico isotrópico con un efecto pepita igual a 0.6, y una meseta de 1.95 para un alcance de 1.
Figura 7.18: Ajuste del modelo 2 (línea sólida) al semivariograma omnidireccional experimental.
95
A los efectos de comparar los modelos se trabaja con el programa XVALID del GEOEAS para realizar la validación cruzada. La Figura 7.19 muestra el mapa de los residuos de la predicción cuando se realiza el ajuste del modelo 2. Al igual que al trabajar con el modelo 1 se observa en la misma que se producen grandes residuos cuando se predicen los valores extremos 7.53 y 11.46. Otros residuos grandes pero de menor magnitud se presentan cuando la predicción (kriging) es realizada en zonas donde la información es escasa.
Figura 7.19: Mapa de los residuos de predicción cuando se ajusta el modelo 2.
La Figura 7.20 muestra los box plots de los residuos (RESIDUO3) y de los residuos estandarizados (RESEST3) de predicción cuando se ajusta el modelo 2. La media de los residuos de predicción estandarizados es 0.011 con un desvío estándar de 1.16 y la mediana es de 0.18. Box Plot (ajuste12345.STA 23v*114c)
RESEST3
RESIDUO3
-12
-10
-8
-6
-4
-2
0
2
4
6
Non-Outlier Max Non-Outlier Min Median; 75% 25% Outliers Extremes Outliers Extremes
Figura 7.20: Box plot de los residuos de predicción y de los residuos de predicción estandarizados para el ajuste 2.
96
Modelo 3. En general el ajuste del modelo al semivariograma muestral puede mejorarse considerando modelos compuesto del tipo: γ (h) = ∑ γ i (h) donde cada uno de los i
sumandos son modelos básicos (exponencial, esférico, etc.). Una propuesta de este tipo es considerar la suma de un efecto pepita de 0.6, el modelo esférico de meseta 0.8 para un alcance de 2 y el modelo esférico de meseta 0.6 para un alcance de 4.
Figura 7.21: Ajuste del modelo3 (línea sólida) al semivariograma omnidireccional experimental.
La Figura 7.22 muestra el mapa de los residuos de la predicción cuando se realiza el ajuste del modelo3.
Figura 7. 22: Mapa de los residuos de predicción cuando se ajusta el modelo 3.
La Figura 7.23 muestra los box plots de los residuos (RESIDUO4) y de los residuos estandarizados (RESEST4) de predicción cuando se ajusta el modelo 3. La media de los residuos de predicción estandarizados es 0.048 con un desvío estándar de 1.33 y la mediana es de 0.24.
97
RESEST4
Non-Outlier Max Non-Outlier Min Median; 75% 25% Outliers Extremes Outliers Extremes
RESIDUO4
-12
-10
-8
-6
-4
-2
0
2
4
6
Figura 7.23: Box plot de los residuos de predicción y de los residuos de predicción estandarizados para el ajuste 3.
El valor del índice de ajuste proporcionado por el programa Variowin es 0.167. Comparación de los modelos. Estos tres ajustes son ejemplos de los infinitos modelos que se podrían ajustar al semivariograma empírico, por el procedimiento de ajuste a sentimiento o cuantificando el índice propuesto por Pannatier en Variowin. El sentido común y el conocimiento físico del fenómeno o de la variable que se estudia, son esenciales a lo largo de todo el proceso de estimación del semivariograma. Ésta, en conjunto con la experiencia del modelador, es una herramienta fundamental para la elección del semivariograma teórico a ajustar. Los tres modelos ajustados dan resultados similares, como se observa al comparar los mapas de residuos de predicción correspondientes y la comparación de la distribución de los residuos estandarizados, que se muestra en la Figura 7.24.
c) RESEST4
b) RESEST3
a) RESEST1
-12
-10
-8
-6
-4
-2
0
2
4
6
Non-Outlier Max Non-Outlier Min Median; 75% 25% Outliers Extremes Outliers Extremes Outliers Extremes
Figura 7.24: Comparación de las distribuciones de los residuos de predicción estandarizados para los tres modelos a) modelo 1, b) modelo 2, c) modelo 3.
Si se tiene en cuenta el índice de ajuste (0.058) proporcionado por el Variowin el modelo 2 sería el más adecuado. Pero este es un modelo demasiado empinado
98
comparado con los otros que alcanzan la meseta en forma más suave. El modelo 1 es el ajuste que considero el más adecuado, porque a falta de información precisa acerca del fenómeno en estudio, la solución debe ser la más sencilla posible. 7.2.3 Isotropía. A esta altura del análisis estructural, es importante chequear la hipótesis de isotropía. Se verá a través de los semivariogramas direccionales si se presenta alguna tendencia de los valores de la variable a lo largo de alguna dirección particular. Todos los pares de puntos que intervienen en el cálculo del semivariograma empírico omnidireccional se dividen en cuatro grupos que aportan al cálculo de los cuatro semivariogramas direccionales en las direcciones de 00, 450, 900, y 1350 con una tolerancia de 22.50. La gráfica de estos cuatro semivariogramas direccionales en las direcciones de 00, 450, 900, y 1350 con una tolerancia de 22.50 se presentan en las Figuras 7.25 - 7.28 respectivamente
Figura 7.25: Semivariograma direccional 00 con una tolerancia de 22.50. Modelo ajustado: modelo1.
Figura 7.26: Semivariograma direccional 450 con una tolerancia de 22.50. Modelo ajustado: modelo1.
99
Figura 7.27: Semivariograma direccional 900 con una tolerancia de 22.50. Modelo ajustado: modelo1.
Figura 7.28: Semivariograma direccional 1350 con una tolerancia de 22.50. Modelo ajustado: modelo1.
Estos semivariogramas direccionales, como todos los semivariogramas direccionales, esconden la mayor parte de la forma, distorsionan los valores del efecto pepita y de la meseta. De esta manera se justifica porque usualmente se modela primero el semivariograma omnidireccional. Sin embargo, a partir de estas cuatro figuras, se puede confirmar el supuesto acerca de la isotropía. Todos los puntos correspondientes a los distintos retardos del semivariograma empírico en la dirección de 00, a excepción del segundo, se encuentran por debajo del modelo 1 ajustado (ver figura 7.25). De esta manera sugiere que en esta dirección el rango debería ser más grande. En las direcciones 450 y 1350, el semivariograma omnidireccional ajustado parece ser adecuado. En la dirección de 900 la distribución de los puntos correspondientes a los distintos retardos del semivariograma empírico, con respecto al modelo1 ajustado, sugiere que el alcance debe ser menor que
100
2. Por lo tanto solo bastaría considerar un rango de alcances, que por ensayo y error se establece entre 1 y 3. Esta información es de utilidad porque cuando se realice la predicción, es decir el kriging, se considera que el alcance del modelo de semivariograma direccional proviene de un patrón elíptico.
7. 3 Kriging. Con el programa Kriging del Geoeas se procede a realizar la predicción. Éste, produce una grilla regular de puntos predichos, usando las ecuaciones del kriging ordinario de la sección 5.3 del capítulo 5. La resolución de dichas ecuaciones permiten encontrar los pesos correspondientes a los valores de la variable que intervienen en el promedio ponderado. A los efectos de disminuir el tiempo de computación las predicciones puntuales se realizan a partir de los valores muestrales que se encuentran en una elipse con centro en el punto a predecir, además se puede especificar la cantidad mínima y máxima de puntos que intervienen en la predicción. El kriging se realiza sobre puntos de una grilla cuyo origen es el punto (0.3, 0.8), es decir 0.3 de la dirección Este y 0.8 del Norte. La separación en la dirección Este es de 1 unidad y en la dirección Norte es de 0.5. Para realizar el kriging es fundamental el conocimiento del semivariograma válido, a través de la etapa de análisis estructural se decidió que el modelo de semivariograma válido es el modelo1 propuesto en la sección 7.2 A partir de los valores predichos de la variable Calcio en la grilla, distintos programas posibilitan dibujar las curvas de nivel, es decir las curvas que tienen igual contenido de calcio. Dichas curvas de nivel se presentan en las Figuras 7.29 y 7.30, la primera es una salida del programa Conrec (componente del Geoeas) y la segunda es producida por el Statistica. Ambos gráficos, tienen en general la misma distribución de las curvas de nivel, claro que el segundo goza de la belleza del color y se notan más claramente las referencias.
N o r t e
Figura 7.29: Curvas de nivel para los valores de kriging de la variable Calcio, producido con el programa Geoeas.
101
Figura 7.30: Curvas de nivel para los valores de kriging de la variable Calcio, producida con el programa Statistica.
Para visualizar la incertidumbre de la predicción se presenta el gráfico de las curvas de nivel para los errores de predicción en las Figuras 7.31 y 7.32, que muestran los resultados razonables: los errores de kriging son más altos en zonas con menor densidad muestral.
Figura 7.31: Curvas de nivel para los valores de los residuos estandarizados del kriging de la variable Calcio, producida con el programa Geoeas
Figura 7.32: Curvas de nivel para los valores de los residuos estandarizados del kriging de la variable Calcio, producida con el programa Statistica.
102
7.4 Logaritmo del contenido de calcio. 7.4.1 Presentación y análisis descriptivo de los datos. En la figura 7.4 se observó que la distribución del contenido de calcio es asimétrica. Si se supone normalidad, las ecuaciones de kriging permiten realizar predicciones con propiedades óptimas. Entonces es apropiado trabajar con una variable transformada para corregir la asimetría y proporcionar a la distribución empírica algunas de las características que posee la distribución teórica. En Hidrogeología, en general a las concentraciones de los iones mayoritarios se les aplica una transformación logarítmica. La figura 7.33 muestra el box- plot de los logaritmos de la variable contenido de calcio. Se observa que presenta una distribución casi simétrica, con un valor alejado. El valor del coeficiente de asimetría es –0.031.
26
LOGCALCI
-1.0
-.5
0.0
.5
1.0
1.5
Figura 7.33: Box-Plot de la variable logaritmo del contenido de Calcio.
Una visualización del buen comportamiento normal de la variable logaritmo del contenido de Calcio es permitida con el Q-Q plot de la Figura 7.34.
Quantile-Quantile Plot of LNCALCIO (lncalci.STA 7v*114c) Distribution: Normal y=0.365+0.638*x+eps .01
.05 .1
.25
.5
.75
.9 .95
.99
3
Observed Value
2
1
0
-1
-2 -4
-3
-2
-1
0
1
2
3
Theoretical Quantile
Figura 7.34: Q-Q plot de la variable logaritmo del contenido de Calcio.
103
4
7.4.2 Semivariograma para la variable logaritmo del contenido de Calcio. El gráfico del semivariograma empírico de la variable transformada presenta menos oscilaciones que el correspondiente a la variable sin transformar. En la figura 7.35 se presenta dicho semivariograma.
Figura 7.35: Semivariograma omnidireccional para la variable logaritmo del Calcio.
Por ensayo de prueba y error se ajusta un modelo esférico isotrópico con un efecto pepita igual a 0.28, y una meseta de 0.14 para un alcance de 4. En la figura 7.36 se muestra dicho ajuste.
Figura 7.36: Ajuste del semivariograma esférico (línea sólida) al semivariograma omnidireccional experimental de la variable logaritmo del Calcio.
La Figura 7.37 muestra el mapa de los residuos de la predicción cuando se realiza el ajuste antes mencionado. Se observa en la misma que se producen grandes residuos cuando se predicen los valores correspondientes a los valores extremos sin
104
transformar: 7.53 y 11.46. Otros residuos grandes pero de menor magnitud se presentan cuando la predicción es realizada en zonas donde la información es escasa.
Figura 7.37: Mapa de los residuos de predicción cuando se ajusta el modelo esférico antes mencionado cuando se trabaja con la variable Logaritmo de contenido de Calcio.
Las distribuciones de los residuos y los residuos estandarizados cuando se realiza la validación cruzada se presentan en la figura 7.38. La misma es un intento de cuantificar el ajuste propuesto para el semivariograma de la variable transformada. Box Plot (lncalci.STA 7v*114c)
RESIDEST
RESIDUO
-5
-4
-3
-2
-1
0
1
2
3
4
Non-Outlier Max Non-Outlier Min Median; 75% 25% Outliers Outliers
Figura 7.38: Distribuciones de los residuos de predicción y de los residuos de predicción estandarizados para el modelo ajustado al semivariograma de la variable transformada.
Las distribuciones de los residuos y de los residuos estandarizados cuando se realiza la transformación inversa se presentan en la figura 7.39.
105
Box Plot (lncalci.STA 12v*114c)
RESEST
RES
-4
-2
0
2
4
6
8
10
12
Non-Outlier Max Non-Outlier Min Median; 75% 25% Extremes Outliers Extremes
Figura 7.39: Distribuciones de los residuos de predicción y de los residuos de predicción estandarizados para el modelo ajustado al semivariograma de la variable transformada.
En la figura 7.40 se comparan las distribuciones de los residuos de predicción estandarizadas cuando se ajusta el modelo1, es decir cuando no se transformó la variable, y cuando se realiza el ajuste del semivariograma de la variable transformada.
Box Plot (juntos.STA 35v*114c)
RESEST
RESEST1
-10
-6
-2
2
6
10
14
Non-Outlier Max Non-Outlier Min Median; 75% 25% Outliers Extremes Outliers Extremes
Figura 7.40: Comparación de las distribuciones de los residuos de predicción estandarizados para el modelo 1, y para cuando se modeló el semivariograma del logaritmo del contenido de calcio.
7.4.3 Kriging de la variable logaritmo del contenido de calcio. A los fines de comparar con los resultados obtenidos en el apartado 7.3, el kriging se realiza sobre puntos de una grilla cuyo origen es el punto (0.3, 0.8), es decir
106
0.3 de la dirección Este y 0.8 del Norte. La separación en la dirección Este es de 1 unidad y en la dirección Norte es de 0.5. Para realizar el kriging es fundamental el conocimiento del semivariograma válido, a través de la etapa del análisis estructural se decidió que el modelo de semivariograma válido es el modelo esférico isotrópico con un efecto pepita igual a 0.28, y una meseta de 0.14 para un alcance de 4. La figura 7.41 muestra las curvas de nivel para el kriging de la variable logaritmo del contenido de calcio. Spline 7 6 5
NORTE
-0.277 -0.212 -0.146 -0.081 -0.016 0.05 0.115 0.181 0.246 0.311 0.377 0.442 0.508 0.573 0.638
4 3 2 1 0 -2
0
2
4
6
8
10
12
14
ESTE
Figura 7.41: curvas de nivel para el kriging de la variable logaritmo del contenido de calcio, obtenido a partir de la grilla indicada.
Las curvas de nivel de las desviaciones estándares del kriging muestran que los valores krigeados en la parte central de la grilla, la que coincide con la presencia de mayor información, se corresponden con los valores más pequeños de las desviaciones estándares. La figura 7.42 presenta dichos resultados. Spline 7 6 5
NORTE
0.578 0.587 0.597 0.606 0.616 0.625 0.635 0.644 0.654 0.663 0.673 0.682 0.692 0.701 0.711
4 3 2 1 0 -2
0
2
4
6
8
10
12
14
ESTE
Figura 7.42: curvas de nivel de las desviaciones estándares del kriging de la variable logaritmo del contenido de calcio, obtenido a partir de la grilla indicada.
107
7.5 Comparación de resultados. Para poder comparar los resultados obtenidos en el apartado anterior con los obtenidos en el apartado 7.3 a partir de la predicción de la variable contenido de calcio sin transformar, se realiza las transformaciones inversas adecuadas de los valores predichos con la variable transformada. La figura 7.43 muestra las curvas de nivel correspondientes a dichos valores. Spline 7 6 5
NORTE
0.878 0.981 1.085 1.188 1.292 1.395 1.498 1.602 1.705 1.809 1.912 2.016 2.119 2.223 2.326
4 3 2 1 0 -2
0
2
4
6
8
10
12
14
ESTE
Figura 7.43: curvas de nivel para los valores transformados de los valores predichos con la variable transformada.
Las curvas de nivel para los valores transformados de las desviaciones estándares del kriging del logaritmo del contenido de calcio se muestran en la figura 7.44. 3D Contour Plot (tuc2.STA 6v*132c) Spline 7 6 5
NORTE
0.552 0.622 0.692 0.762 0.832 0.902 0.971 1.041 1.111 1.181 1.251 1.321 1.39 1.46 1.53
4 3 2 1 0 -2
0
2
4
6
8
10
12
14
ESTE
Figura 7.44: curvas de nivel de los desvíos estándares al aplicar la transformación inversa adecuada.
En todos los mapas presentados los resultados tienen validez en el rectángulo definido por la grilla de puntos definida para realizar el kriging es más, dada la irregularidad de la distribución de los puntos donde se observó la variable, “la región 108
válida de predicción” sería un polígono irregular. Esta grilla tiene sus límites basados en el rango de las coordenadas Este y Norte de los puntos de observación de la variable de interés. Para obtener otras predicciones se puede variar los valores que definen la grilla. En el siguiente mapa se presentan la grilla de puntos donde se realiza las predicciones y la ubicación de los puntos donde se observó la variable de interés. También se marca las regiones donde las predicciones no tendrían ningún valor, porque los valores extrapolados no tienen mucha confianza. De acuerdo a los resultados presentados en el apartado 5.5.1 del capítulo 5, se observa que en la zona donde se dispone más información los desvíos estándares de la predicción son menores, dando más confianza a los valores predichos. Scatterplot (tuc2.STA 6v*132c) 7
6
NORTE
5
4
3
2
1
0 0
2
4
6
8
10
12
ESTE
Figura 7.45: mapa donde se muestran la grilla de predicción y las ubicaciones de los puntos donde se observó la variable de interés.
Por supuesto que el conocimiento geológico del problema tratado ayudaría a tratar de determinar la significación de los resultados. Es ahí cuando la interdisciplina alcanza su punto máximo. De la comparación de las curvas de nivel presentadas en las figuras 7.30 y 7.43 se observa que los valores predichos con la modelación del semivariograma de la variable transformada son altos en zonas o regiones que se esperarían valores más bajos, esta suposición es más respetada cuando la predicción es llevada a cabo sin la transformación de la variable. Así, los resultados del kriging obtenidos en el apartado 7.3 con la variable sin transformar son más razonables que los obtenidos con la variable transformada. Estos resultados no hablan en contra de la transformación de las variables, pero muestra que es conveniente realizar ambos tratamientos a los fines de decidir por cual herramienta es más útil. Aparte que la interpretación de los resultados a partir de la variable transformada son más difíciles de realizar.
109
CAPITULO 8: APLICACIÓN A LA HIDROGEOLOGIA: “DATOS DEL ACUÍFERO DE LA CALDERA”
8.1 Introducción. El principal objetivo de este capítulo, es aplicar la metodología Geoestadística en el modelado de las variables hidrogeológicas de interés, que están relacionadas con la calidad química de las aguas subterráneas. Para el desarrollo de la investigación, se contó con la información proporcionada por la cátedra de Hidrogeología en el marco del proyecto de investigación Hidrogeología del Sistema Acuífero La Caldera. De acuerdo a los especialistas, las variables químicas consideradas de importancia para la caracterización del sistema acuífero de La Caldera, a los efectos de determinar posibles focos de contaminación antrópica, son la concentración de cloruros y la concentración de nitratos de las aguas subterráneas. En la región de interés, es decir lo que se denomina sistema acuífero de La Caldera, y para la época (1991-1995) considerada como unidad temporal de análisis, se consideran las planillas de análisis físico- químico de 47 pozos de los que estaban en funcionamiento, para alimentar el sistema de agua potable del sector norte y el centro de la ciudad de Salta. Una de las primeras tareas durante la investigación fue la elaboración de la base de datos. En ella principalmente se registró la información de la ubicación de los pozos donde se realizaron las observaciones de las variables de interés. Éstas, también se clasificaron de acuerdo a los meses en que fueron observadas. Al registrar la información de primera mano, proveniente de las fichas de los análisis químicos de las muestras de agua, tomadas de los diferentes pozos junto a los legajos de perforaciones realizadas por la Dirección General de Obras Sanitarias, se encontró dificultades para determinar la ubicación exacta de los pozos. Pero, con la colaboración de miembros del equipo de investigación, fueron subsanadas. Aquellos casos en que un pozo presentó en un mismo día o dentro del mismo mes más de una observación, se tomó como representativa el mayor valor de concentración para cada ión. Las épocas para las que se dispone de información son los meses de julio, agosto y setiembre de 1991; agosto, setiembre y octubre de 1993. Algunas pocas fichas aportan datos muy aislados en los meses de octubre y noviembre de 1992 (un total de 8 pozos). Otros muy pocos datos se presentan en los primeros meses de 1991, 1992 y 1993, y solo dos observaciones se encuentran disponibles en el año 1995. Así, a los fines de intentar realizar el análisis geoestadístico se consideró a toda la información como proveniente de un solo momento geológico. Por otro lado, el pozo ASP1232 presenta datos como para realizar un estudio transversal, esto es, se podría estudiar el comportamiento de las variables de interés a través del tiempo, el cual no es el motivo de este trabajo. Este pozo era seguido con los análisis físico-químico mes a mes, por corresponder a una planta de alimentos balanceados.
110
8.2 Ubicación de los pozos. El mapa que se muestra en la figura 8.1, presenta las ubicaciones de los 46 pozos donde se observaron los valores de las variables de interés: concentración de cloruros y concentración de nitratos. La distribución espacial de los mismos es irregular.
N O R T E
Figura 8.1:Mapa que muestra la posición de los pozos de agua donde se observaron las variables de interés.
ESTE
A los efectos de poder trabajar con los programas específicos, las coordenadas Gauss Kruger que determinan la posición de los pozos fueron modificadas restando a las coordenadas este 3500 y a las del norte 7200. De los 46 pozos cuyas fichas aportaron datos de concentración de cloruros, 3 fichas (correspondientes a los pozos: AS0129, AS393 y AS0548) no aportan datos de la concentración de nitratos.
8.3 Concentración de cloruros. 8.3.1 Análisis descriptivo. El siguiente diagrama de puntos presenta los valores de la concentración de cloruros observados.
Figura 8. 2:Scatter plot de la concentración de cloruros, las unidades del eje vertical es mg/lt.
111
En el gráfico de la figura 8.2 se nota claramente el valor alejado aportado por el pozo ASP1232. También se observa una leve tendencia lineal, en la dirección Sur- Norte. Las posibles tendencias lineales se investigarán a través de los scatter-plot de la concentración de cloruros versus cada una de las variables coordenadas.
Figura 8.3: Diagrama de puntos de la variable concentración de Cloruros vs. la variable Este.
En la figura 8.3, se observa que no existe dependencia lineal entre la concentración de cloruros y la variable Este, cuyos valores son las coordenadas este de los puntos de posición. Lo cuál queda confirmado con un p-valor muy alto (0.9667) que no permite rechazar la hipótesis que la pendiente sea nula (a=-.07075). En cambio en la figura 8.4 se observa una leve tendencia, a medida que aumenta los valores de la coordenada norte disminuye la concentración de cloruros. Pero es sólo efecto de las escalas porque a través de un test de hipótesis (a=-1.6941, p-value=0.0119), se rechaza la existencia de tendencia.
Figura 8.4: Diagrama de puntos de la variable concentración de Cloruros vs. la variable Norte.
112
Es también importante destacar que en la figura 8.5 se observa una fuerte dependencia lineal entre la concentración de nitratos y la concentración de cloruros, que permitiría aumentar la cantidad de datos de la primer variable prediciéndola a partir de los valores de la segunda en los pozos ya citados.
Figura 8.5: Diagrama de puntos de la variable concentración de Nitratos vs. la variable concentración de Cloruros.
Para este ajuste se obtuvo un coeficiente de correlación igual a 0.744 y un test de hipótesis para testar si la pendiente es nula es rechazado con un p-valor= 0.00000 para un valor estimado de la pendiente igual a 0.601848. Al dejar de lado la posición y solo considerar la variable concentración de cloruros se observa que su distribución es asimétrica positiva y presenta un valor alejado: 71.84 mg/l, el correspondiente al pozo ASP1232. Esto, se observa en el box plot siguiente.
Figura 8.6: Box plot de la variable concentración de cloruros.
113
Para algunos modelistas es fundamental realizar una transformación de variables a los efectos de lograr simetría en la distribución, característica fundamental de la distribución teórica normal. Porque, si las medidas de concentración de cloruros provienen de una distribución normal, las predicciones obtenidas con el kriging ordinario son óptimas (Cressie, 1991). En este apartado, se trabajará con los datos sin transformar. 8.3.2 Análisis estructural. 8.3.2.1 Estimación del variograma. En este apartado, se intenta sobre la base de los datos disponibles determinar la estructura de correlación espacial del proceso aleatorio Z(x): “concentración de cloruros”. Bajo el supuesto que el proceso Z(x) es intrínsecamente estacionario, la función semivariograma γ(h) cuantifica dicha estructura. A partir de una realización se debe estimar el semivariograma. Con el programa Vario que forma parte del Geoeas se genera el semivariograma experimental “omnidireccional”, el cual permite utilizar todos los pares de datos (761 pares) independientemente de la dirección. Se considera la dirección de 0° y una tolerancia angular de 90°. Las siguientes condiciones para las distancias: Mínimo =0, Máximo = 5 (km), es decir se toma en cuenta la mitad de la máxima distancia entre los puntos, y un incremento de 0.5 (km), permiten generar 11 intervalos planos de distancia. En la figura 8.7 se presenta el semivariograma experimental “omnidireccional”.
Figura 8.7: Semivariograma experimental “omnidireccional” para la concentración de cloruros.
Para observar que las elecciones de los intervalos de distancias no influyen en la “forma” del semivariograma muestral se realizaron sucesivos gráficos con distintas longitudes de los intervalos de distancia, pero la “forma” del mismo no varió sustancialmente.
114
8.3.2.2 Ajuste a un modelo de semivariograma. La idea es buscar un semivariograma válido que represente más ajustadamente a la dependencia espacial presente en los datos: Z = (Z(s1), ..., Z(sn) )’. El espacio de todos los semivariogramas válidos es un gran conjunto, usualmente se elige una familia paramétrica de semivariogramas. Para la variable en cuestión de acuerdo a la forma de los semivariogramas empíricos se elige la familia de semivariogramas esféricos. Con el programa Geoeas se realiza un ajuste a sentimiento. El cual no garantiza un modelo de semivariograma único ya que se basa en apreciaciones subjetivas y en la experiencia del usuario. Es de destacar que no se debe intentar ajustar los mínimos detalles ya que en general éstos no son una característica del verdadero semivariograma sino más bien fluctuaciones muestrales. Por sucesivas pruebas de ensayos de prueba y error se ajustó el semivariograma esférico: 0 si h= 0 3 3 h h 1 γˆ (h;θˆ) = 65 + 120 * − si 0 < h ≤ 1.6 2 1 .6 2 1 .6 185 si h ≥1.6 En la figura 8.8 se muestra dicho ajuste. Es decir se ajusta un modelo esférico isotrópico con un efecto pepita igual a 65 (unidades de concentración al cuadrado), y una meseta de 185 (unidades de concentración al cuadrado) para un alcance de 1.6 km.
Figura 8.8: Ajuste del semivariograma esférico (línea sólida) al semivariograma omnidireccional experimental.
A los efectos de determinar la calidad y el grado de fiabilidad del modelo ajustado se realizará la validación cruzada. Así se considera un nuevo conjunto de datos: las diferencias entre las predicciones y los valores actuales correspondientes a las posiciones donde se realizaron las observaciones; es decir los residuos de la predicción. La Figura 8.9 muestra el mapa de los residuos de la predicción cuando se realiza el ajuste antes mencionado. Se observa que se produce un residuo elevado (los símbolos son proporcionales a los residuos) cuando se predice el valor correspondiente al pozo 115
ASP1232. Otros residuos grandes pero de menor magnitud se presentan cuando la predicción (kriging) es realizada en zonas donde la información es escasa.
Figura 8.9: Mapa de los residuos de predicción cuando se ajusta el modelo esférico propuesto.
8.3.3 Isotropía. A esta altura del análisis estructural, es importante chequear la hipótesis de isotropía. A través de los semivariogramas direccionales se investiga si se presenta alguna tendencia de los valores de la variable a lo largo de alguna dirección particular. Todos los pares de puntos que intervienen en el cálculo del semivariograma empírico omnidireccional se dividen en dos grupos que aportan al cálculo de dos semivariogramas direccionales en las direcciones de 00 y 900, con una tolerancia de 450. Las gráficas de estos semivariogramas direccionales se presentan en las figuras 8.10 y 8.11 respectivamente.
Figura 8.10: Semivariograma direccional 00 con una tolerancia de 450.
116
Figura 8.11: Semivariograma direccional 900 con una tolerancia de 450
Como todos los semivariogramas direccionales, éstos esconden la mayor parte de la forma del semivariograma, distorsionando los valores del efecto pepita y la meseta. Sin embargo, a partir de estas dos figuras, se puede confirmar el supuesto acerca de la isotropía. En ambas direcciones, el semivariograma omnidireccional ajustado parece ser adecuado. En la dirección de 00 la distribución de los puntos correspondientes a los distintos retardos del semivariograma empírico, con respecto al modelo ajustado, sugiere que el alcance debe ser mayor que 1.6. Por lo tanto solo bastaría considerar un rango de alcances, que por ensayo y error se establece entre 1.2 y 2.5. Esta información es de utilidad porque cuando se realice la predicción, se considera que el alcance del modelo de semivariograma direccional proviene de un patrón elíptico. 8.3.4 Kriging. Con el programa Kriging del Geoeas se procede a realizar el Kriging (la predicción). Éste, produce una grilla regular de puntos predichos, usando las ecuaciones del kriging ordinario para encontrar los pesos correspondientes a los valores de la variable que intervienen en el promedio ponderado. El kriging se realiza sobre puntos de una grilla cuyo origen es el punto (56, 57), es decir 56 de la dirección Este y 57 del Norte. La separación de los puntos, que forman la grilla, en ambas direcciones es de 0.5 km. Para realizar el kriging es fundamental el conocimiento del semivariograma válido, a través de la etapa de análisis estructural se decidió que el modelo de semivariograma válido es el modelo esférico propuesto en la sección 8.3.2.2. A partir de los valores predichos de la variable Concentración de cloruros en la grilla, distintos programas posibilitan dibujar las curvas de nivel, es decir las curvas cuyos puntos tienen igual concentración de cloruros. Dichas curvas de nivel se presentan en la figura 8.12 que es una salida del programa Conrec (parte del Geoeas).
117
N O R T E
ESTE Figura 8.12: Curvas de nivel para los valores de kriging de la variable concentración de cloruros, producido con el programa Geoeas.
Figura 8.13: Curvas de nivel para los valores de kriging de la variable concentración de cloruros, producido con el programa Geoeas, en una grilla de más resolución.
118
La Figura 8.13 muestra las curvas de nivel pero para otra grilla con mayor resolución (2.5 km. en ambas direcciones) donde se realiza el kriging. Los errores de predicción se presentan en el mapa de la figura 8.14.
Figura 8.14: Curvas de nivel para los valores de los errores estándares del kriging de la variable concentración de cloruros, producido con el programa Geoeas.
Los mapas de curvas de nivel de las predicciones de concentración de cloruros muestran que el pozo ASP1232 es un posible foco de contaminación. Porque se predicen valores altos en cercanías de dicho pozo, cuyas coordenadas en este mapa son (59.31, 57.38). Son valores altos, pero si se tuviese más información en posiciones cercanas a dicho pozo, las predicciones tendrían más precisión, es decir las curvas de nivel de los errores de predicción en esta zona llegarían a tomar valores más bajos.
8.4 Nitratos. 8.4.1 Análisis descriptivo. El siguiente post- plot muestra el comportamiento de la variable concentración de nitratos a través de los 43 pozos que aportaron información. Se observa un comportamiento similar a la de la concentración de los cloruros, los valores más bajos de concentración de nitratos se presentan en la parte norte del acuífero. En cambio los valores más altos se encuentran diseminados en la parte sur. El valor correspondiente al pozo ASP1232 no tiene muchas observaciones cercanas, así que cuando se realice la validación cruzada es muy probable que se observe un residuo muy grande.
119
Figura 8.15: Ubicación de los pozos donde se observó la variable concentración de nitratos(mg/l) y el valor de la misma categorizada según los cuartiles.
El box - plot siguiente muestra la distribución asimétrica de las observaciones de la variable. En él se destaca el valor alejado: 61.92 mg/l correspondiente al pozo hasta aquí cuestionado, los demás valores no superan el valor 39 mg/l.
Figura 8.16: Box- plot de la variable concentración de nitratos.
A los efectos de detectar posibles tendencias en las direcciones pre-especificadas se grafican los scatter – plot: Nitratos vs Este y Nitratos vs Norte.
120
Figura 8.17: Diagrama de puntos de la variable concentración de Nitratos vs. la variable Este.
En la figura 8.17, se observa que no existe dependencia lineal entre la concentración de nitratos y la variable Este, cuyos valores son las coordenadas este de los puntos de posición. Lo cuál queda confirmado con un p-valor muy alto (0.683768) que no permite rechazar la hipótesis que la pendiente sea nula (a=0.155854). En cambio en la figura 8.18 se observa una leve tendencia, a medida que aumenta los valores de la coordenada norte disminuye la concentración de nitratos. Pero es sólo efecto de las escalas porque a través de un test de hipótesis (a=-1.6941, p-value=0.0119), se rechaza la existencia de tendencia.
Figura 8.18: Diagrama de puntos de la variable concentración de Nitratos vs. la variable Norte.
121
8.4.2 Análisis estructural. 8.4.2.1 Estimación del variograma. Con el programa Vario se genera el semivariograma experimental “omnidireccional”. Se utilizan los mismos parámetros para la generación de los intervalos planos que en la definición del semivariograma muestral de la variable concentración de cloruros. De esta manera se generan 11 intervalos de distancia. En la figura 7.19 se presenta el semivariograma experimental “omnidireccional”.
Figura 7.19: Semivariograma experimental “omnidireccional”.
Para observar que las elecciones de los intervalos de distancias no influyen en la “forma” del semivariograma muestral se realizaron sucesivos gráficos con distintas longitudes de los intervalos de distancia, pero la “forma” del mismo no varió sustancialmente. 8.4.2.2 Ajuste a un modelo de semivariograma. Con el programa Geoeas se realiza un ajuste a sentimiento. Por sucesivas pruebas de ensayos de prueba y error se ajustó el semivariograma esférico: 0 si h= 0 3 3 h 1 h si 0 < h ≤ 2.5 γˆ (h;θˆ) = 5 0+ 125 * − 2 2 . 5 2 2 . 5 175 si h ≥ 2 .5 En la figura 8.20 se muestra dicho ajuste. Es decir se ajusta un modelo esférico isotrópico con un efecto pepita igual a 50 (unidades de concentración al cuadrado), y una meseta de 175 (unidades de concentración al cuadrado) para un alcance de 2.6 km.
122
Figura 8.20: Ajuste del semivariograma esférico (línea sólida) al semivariograma omnidireccional experimental.
A los efectos de determinar la calidad y el grado de fiabilidad del modelo ajustado se realiza la validación cruzada. La Figura 8.21 muestra el mapa de los residuos de la predicción cuando se realiza el ajuste antes mencionado. Se observa en la misma que se produce un residuo elevado cuando se predice el valor correspondiente al pozo ASP1232. Otros residuos elevados se presentan cuando la predicción es realizada en zonas donde la información es escasa.
Figura 8.21: Mapa de los residuos de predicción cuando se ajusta el modelo esférico propuesto.
123
A pesar de trabajar con muy poca información, a través de los semivariogramas direccionales se investigó si se presenta alguna tendencia de los valores de la variable a lo largo de alguna dirección particular, sin obtener resultados positivos. 8.4.3 Kriging. Con el programa Kriging del Geoeas se procede a realizar la predicción. El kriging se realiza sobre puntos de una grilla cuyo origen es el punto (56, 57), es decir 56 de la dirección Este y 57 del Norte. La separación de los puntos en ambas direcciones es de 0.25 km. A partir de los valores predichos de la variable Concentración de nitratos, se dibujan las curvas de nivel. Dichas curvas de nivel se presentan en la figura 8.22.
Figura 8. 22: Curvas de nivel para los valores predichos de la variable concentración de nitratos.
El mapa de curvas de nivel de las predicciones de concentración de nitratos muestran que el pozo ASP1232 es un posible foco de contaminación.
124
Para visualizar la incertidumbre de la predicción se presenta el gráfico de las curvas de nivel para los errores de predicción en la Figura 8.23, que muestra los resultados razonables: los errores de kriging son más altos en zonas con menor densidad muestral.
Figura 8. 23: Curvas de nivel para los valores de los residuos estandarizados del kriging de la variable concentración de nitratos, producida con el programa Geoeas.
8.5 Mapas generados por el programa Surfer. A los efectos de presentar mapas con mejores definiciones que los producidos por el programa Geoeas, a partir del modelado del semivariograma para los cloruros y nitratos se realizó el Kriging con el programa Surfer, muy utilizado para realizar mapas, y se procedió a representar las curvas de nivel. En las figuras 8.24 y 8.25 se presentan los mapas de iso-contenidos de cloruros y nitratos generados por el programa Surfer. A estos mapas hay que recortarles las curvas de nivel en la zona del noroeste indicada con rojo en ambos mapas, porque en esa zona no se presenta ninguna observación y por lo tanto las predicciones en dichas regiones no tienen validez.
125
Figura 8.24: Mapa de isolíneas para la concentración de cloruros.
126
7266
7265
7264
7263
7262
7261
7260
7259
7258 3557
3558
3559
3560
Figura 8. 25: Mapa de isolíneas para la concentración de nitratos.
127
3561
8.6 Conclusiones. A pesar de contar con una cantidad de observaciones inferior a lo recomendado para realizar un estudio geoestadístico, los mapas de curvas de nivel logrados a través de este procedimiento presentan características fundamentales que los especialistas intuían “acerca de la realidad”. Aunque, la realidad nunca llega a ser conocida. Las curvas de igual contenido de cloruros presentadas en el mapa correspondiente, al igual que las de igual concentración de nitratos, muestran que el pozo ASP1232 es un posible foco de contaminación. Así los pozos cercanos están influenciados por sus valores altos de concentración de cloruros y de nitratos. Si se perforasen otros pozos en la vecindad, en un radio de 1km, se obtendrían valores similares de los indicadores. Otros focos posibles de contaminación, para las aguas subterráneas del acuífero de La Caldera, pero de menor grado se encuentran en la zona norte. A juicio de los especialistas esto es un llamado de atención para la actual empresa encargada de proveer el agua a los habitantes del sector norte y el centro de la ciudad de Salta a los efectos de la perforación de nuevos pozos, pero es muy poco determinante en la calidad de las aguas que llegan al consumidor porque la pureza del agua provista por el acuífero se ve afecta en gran medida por el sistema de distribución obsoleto. Las curvas de nivel de los desvíos estándares de las predicciones muestran valores altos en las zonas cercanas al valor alejado correspondiente al pozo cuestionado, como así también, en zonas donde la información es escasa. Dichos valores también dependen del semivariograma elegido, pero por ensayo y error, modelos distintos a los adoptados produjeron valores más elevados.
128
CAPITULO 9: UNA APLICACIÓN A LAS CIENCIAS DEL MEDIO AMBIENTE 9.1 Introducción. En este Capítulo se pretende aplicar las herramientas proporcionadas por el enfoque estadístico presentado en los 6 primeros capítulos a un problema proveniente de las denominadas ciencias del medio del ambiente. Para ello se utilizará información suministrada por el equipo de investigación del proyecto 577: “Determinación de Dióxido de Nitrógeno en la Atmósfera de Salta (capital)” de la Facultad de Ciencias de Exactas realizados entre los años 1996 y 1998. Este grupo de investigación indaga acerca de la calidad del aire que se respira en la ciudad de Salta y por lo tanto se interesan en las concentraciones de los contaminantes mayoritarios en la troposfera. El dióxido de nitrógeno es uno de los componentes mayoritarios de la troposfera de las ciudades modernas debido, principalmente, a los procesos de combustión, como los que ocurren en los motores de combustión interna de los vehículos automotores. Estos producen óxidos de nitrógeno que son descargados directamente a la atmósfera. La producción de NO X en los motores de los automotores es inevitable bajo las condiciones de baja temperatura y presión de funcionamiento. En un automóvil cuya mezcla aire/combustible es la adecuada, la producción de CO e hidrocarburos sin quemar es mínima pero la de NO es máxima. (Informe final del proyecto 567). Los óxidos de nitrógeno, y particularmente el dióxido de nitrógeno pueden impactar negativamente en la salud humana y contribuir a la degradación del medio ambiente por generación de lluvia ácida y formación del nebluno fotoquímico. El tránsito intenso de vehículos no debidamente controlados en cuanto a las condiciones de emisión de los gases de escape, y las características topográficas de la Ciudad de Salta, plantearon la necesidad de determinar la concentración de NO X en el estrato inferior de la atmósfera. La contribución de este trabajo consiste en proporcionar mapas de la concentración de dióxido de nitrógeno que permita predecirla en lugares donde no se realizó observación por falta de instrumentos debido a dos causas fundamentales: ∧ el costo de los dispositivos y de los componentes químicos para la detección de la sustancia química de interés. ∧ por la destrucción de los mismos ocasionados por la ignorancia acerca de la utilidad de los dispositivos instalados para tal fin.
9.2 Presentación de los datos. Las observaciones fueron realizadas en determinados puntos estratégicos del micro, del macro centro, y la periferia de la ciudad de Salta. A los efectos de tener como referencia o medida basal local se tomaron mediciones de la concentración de dióxido de nitrógeno en un ambiente netamente rural en el sitio denominado La Choza. Además se monitorearon cuatro vías de acceso a la ciudad. Para la determinación de NO2 en la atmósfera los investigadores utilizaron la técnica del muestreador pasivo, basado en la difusión del gas en tubos colectores. La descripción de los mismos se encuentran en el informe final del proyecto de investigación n°577. Los tubos colectores fueron ubicados en las veredas de las calles
129
cerca de las esquinas siguiendo las dos direcciones Norte-Sur y Este-Oeste favorecidos por la distribución de las calles en el casco principal de la ciudad de Salta. En el siguiente mapa se presentan las 22 posiciones de los puntos donde se logró obtener información de la concentración de dióxido de nitrógeno en febrero de 1997 (uno de los periodos de recolección de 15 días de duración). A los efectos de determinar las coordenadas de los puntos de observación en el mapa se consideró un sistema de coordenadas ortogonales con centro en un punto referencia de la ciudad: el vértice nortdeste de la plaza principal 9 de Julio. El eje de las abscisas en la dirección OesteEste y el eje de las ordenadas en la dirección Sur- Norte. 1600
1100
NORTE
600
100
-400
-900
-1400 -1600
-1100
-600
-100
400
900
1400
OESTE
Figura 9.1: Posiciones donde se realizaron las observaciones.
Las distancias están medidas en metros. Así, en metros, están expresadas las variables que definen ambos ejes del gráfico de la figura 9.1. Los puntos de observación se distribuyen en un rectángulo de aproximadamente 2.5 kilómetros en la dirección Oeste- Este por 2.7 kilómetros en la dirección Sur- Norte. Algunos puntos fueron cambiados de sus posiciones originales a los efectos de que formen parte de una grilla irregular de 6 filas por 7 columnas. En la figura 9.2, los puntos correspondientes al microcentro están indicados con triángulos, con cuadrados los del macrocentro y los restantes corresponden a lugares periféricos. De acuerdo a la cantidad promedio de vehículos que circulan por las calles de la ciudad, se observa que existe más concentración de puntos de observación en la zona denominada centro de la ciudad. Aproximadamente el 67% cumple con dicho requisito. Claro que la densidad de puntos es más alta cuanto más nos acercamos al centro comercial de la ciudad (microcentro). También en la Figura 9.2, con la letra T queda indicada la posición de un punto de observación en la zona de la Terminal de ómnibus, con M la del Mercado municipal y con A posiciones de puntos de observación en principales avenidas de acceso al centro comercial de la ciudad.
130
A
A
A
A
M T A A
Figura 9.2: Mapa donde se indican las posiciones de puntos estratégicos de observación.
La variable de interés: concentración de dióxido de nitrógeno, está medida en
µg de NO2 por cada m3. En la figura 9.3 se presenta las concentraciones del NO2 observadas en febrero de 1997 en sus respectivas posiciones. 1600
1.2
.8
2.4
3.4
1100
600
8.5
11.5
7.4
2.3
2.7
11.0
100 .4
4.3
5.5 4.0
NORTE
-400
.8
-900 .2
3.3
8.7
2.6
2.3
1.9
.3
-1400 -1600
-1100
-600
-100
400
900
1400
OESTE
Figura 9.3: Medidas de la concentración de dióxido de nitrógeno en sus respectivas posiciones correspondientes a la ciudad de Salta.
El punto que representa a la zona de la terminal de ómnibus al igual que el que representa al mercado municipal presentan los valores más elevados 11 µg / m 3 y 11.5 µg / m 3 respectivamente. Sin considerar la posición de las observaciones, en la Figura 9.4 se observa que la distribución de los valores de la variable de interés es asimétrica, no presenta valores alejados. El valor mínimo de las concentraciones de dióxido de nitrógeno es .20 µg / m3 , el valor máximo es 11.5 µg /m3 , el valor medio es 4.07 µg / m 3 y el valor mediano es 2.65 µg / m 3 .
131
Figura 9.4: Box Plot de la concentración de dióxido de Nitrógeno.
Los scatter plot de la variable de interés versus cada una de las variables que definen la posición de los puntos, mostrados en la Figura 9.5, señalan que no se presenta tendencia en los valores de la variable concentración de nitrógeno ni en la dirección Oeste-Este ni en la dirección Sur-Norte. b)
14
14
12
12
10
10
8
8 NOX297
NOX297
a)
6 4
4 2
2
0
0 -2 -2000
6
-1500
-1000
-500
0
500
1000
-2 -2000
1500
OESTE
-1500
-1000
-500
0
500
1000
1500
NORTE
Figura 9.5: Diagramas de puntos de a) concentración de nitrógeno vs. Coordenada oeste- este. b) concentración de nitrógeno vs. Coordenada sur-norte.
A los efectos de detectar posibles tendencias y calcular el estadístico u se procede a calcular las medias y medianas por filas y columnas.
Figura 9.6: En ambas figuras con azul se representan los valores de la media y con rojo los valores de la mediana de los valores de concentración de nitrógeno. a) Según columnas. b) Según filas.
La Figura 9.6 es un intento de resumir la posible no estacionariedad en la media a lo largo de filas y columnas, es decir en las direcciones sur- norte y oeste- este, usando 132
2000
la media muestral y la mediana muestral a lo largo de las filas y las columnas respectivamente. En la Figura 9.6 a) se observa que la columna correspondiente a la coordenada oeste-este –175 m se destaca de las otras que presentan una pequeña tendencia lineal. En cuanto a la Figura 9.7 b) muestra que existe una leve tendencia no lineal sino cuadrática. Luego, por el método de la mediana polish se tratará de modelar dicha tendencia. En las tablas 9.1 y 9.2 se presentan los valores del estadístico u definidos en la sección 2.2.3.
Tabla 9.1: Valores del estadístico u según columnas.
Tabla 9.2: Valores del estadístico u según filas.
En la tabla 9.1 el valor del estadístico u señala que en las columnas correspondientes a los valores de la variable oeste-este la variable en estudio no presenta outliers. De forma similar ocurre en la tabla 9.2 para las filas. Así las diferencias estandarizadas entre medias y medianas muestran que por fila (o columna) la variable no presenta valores alejados.
9.3 Modelado de la tendencia. Mediante el método de la mediana polish se procederá a modelar la tendencia que presenta la variable concentración de dióxido de nitrógeno en la región de estudio. Luego de 20 iteraciones del proceso de extracción de los efectos medianos según fila y según columna, de acuerdo a la rutina indicada en el apartado 6.3.1, se obtuvieron los resultados indicados en la tabla 9.3.
Efectos medianos Filas
Residuos
Efectos medianos columnas
Efecto Mediano Total
Tabla 9.3: Valores de los efectos medianos según filas, columnas; efecto mediano total y residuos luego de aplicar Mediana polish a las concentraciones de dióxido de nitrógeno.
133
La superficie mediana polish formada por la unión de planos definidos sobre los rectángulos que determinan la grilla con ecuaciones 6.4.11 se muestra en la figura 9.8. En dicha representación se observa como se modelo la tendencia que se presenta particularmente en la dirección Sur- Norte.
Figura 9.7: Superficie de planos interpolados usando mediana polish para los datos de concentración de dióxido de nitrógeno.
9.4 Kriging de los residuos. Los residuos obtenidos mediante la aplicación del procedimiento de mediana polish se consideran como un nuevo conjunto de datos espaciales sin presencia de tendencia. Por lo tanto, se puede llevar acabo la predicción de los residuos en posiciones determinadas mediante el procedimiento de kriging ordinario. Luego ésta se sumará a la tendencia para obtener la predicción de la concentración de dióxido de nitrógeno en los puntos de interés, de acuerdo a la descomposición aditiva presentada en 6.4.10. Para ello se deben seguir las distintas etapas que se desarrollan a continuación. 9.4.1 Estimación del variograma. En este apartado, se intenta sobre la base de los datos disponibles determinar la estructura de correlación espacial del proceso aleatorio R(x):“residuos obtenidos al aplicar mediana polish”. Con el programa Vario (Geoeas) se genera el semivariograma experimental “omnidireccional”, el cual permite utilizar todos los pares de datos (284 pares) independientemente de la dirección. Se considera la dirección de 0° y una tolerancia angular de 90° y para completar la definición de los intervalos planos las siguientes condiciones para las distancias: Mínimo 245, Máximo = 1800 (m) es decir la mitad de la máxima distancia entre los puntos, y el incremento se fija en 300 m, de esta manera se generan 6 intervalos de distancia. En la figura 9.8 se presenta el semivariograma experimental “omnidireccional” y en la tabla 9.4 se muestran las cantidades de pares que determinan el valor del semivariograma correspondiente a los intervalos de 134
distancias. La cantidad de los mismos en los dos primeros intervalos planos es menor que 30 pares.
Figura 9.8:Semivariograma experimental “omnidireccional” para los residuos usando como parámetros para las distancias: mínima=245 m; máxima=1800 m y un incremento 300 m.
1 2 3 4 5 6
Pairs 10 22 39 34 30 33
Avg Distance 418.671 682.636 986.583 1264.834 1585.817 1884.789
Estimate 6.508 7.196 8.884 7.965 8.202 7.213
Tabla 9.4: Valores de la cantidad de pares, distancias promedios y la estimación del semivariograma correspondiente.
Los semivariogramas direccionales no aportan mucha información debido a que la cantidad de observaciones disponible hace que la cantidad de pares de puntos que intervienen en el cálculo de la estimación del semivariograma en cada intervalo plano de distancia sea pequeña. 9.4.2 Ajuste a un modelo de semivariograma. Con el programa Geoeas se realiza un ajuste a sentimiento. Por sucesivas pruebas de ensayos de prueba y error se ajustó el semivariograma esférico: 0 si h= 0 3 3 h ) 1 h ) si 0 < h ≤ 800 γ (h;θ ) = 2.4+ 5.6 * − 2 800 2 800 8 si h ≥ 800 En la figura 9.9 se muestra dicho ajuste. Es decir se ajusta un modelo esférico isotrópico con un efecto pepita igual a 5.6 (unidades de concentración al cuadrado), y una meseta de 8 (unidades de concentración al cuadrado) para un alcance de 800 m. A los efectos de determinar la calidad y el grado de fiabilidad del modelo ajustado se realiza la validación cruzada. Para ello se considera un nuevo conjunto de datos: las diferencias entre las predicciones y los valores actuales correspondientes.
135
Figura 9.9: Ajuste del semivariograma esférico (línea sólida) al semivariograma omnidireccional experimental.
La Figura 9.10 muestra el mapa de los errores de la predicción cuando se realiza el ajuste propuesto. Se observa en la misma que se producen dos diferencias elevadas (los símbolos son proporcionales a la diferencia entre el valor predicho y el valor de la variable residuos). Una, cuando se predice el valor correspondiente a la posición (-175, 175) que presenta un valor alejado de la variable residuo igual a –6.03; la otra diferencia alta positiva se corresponde al punto (805, 420) en donde la variable presenta el valor alejado 6.25. Es de destacar que el programa no puede predecir en las posiciones (-420, 1400); (-1400, 1400) ; (1050, -1330). Las dos primeras son el vértice superior izquierdo y su vecino más cercano, que no presentan vecinos cercanos para realizar la predicción. Lo mismo sucede con el último, pero su posición es el vértice inferior derecho (del rectángulo que definen los datos).
Figura 9.10: Mapa de los residuos de predicción cuando se ajusta el modelo esférico propuesto.
136
Seguidamente a los efectos de comparar se realizó un ajuste usando el criterio de mínimos cuadrados ponderados. Al aplicar un ajuste mediante mínimos cuadrados ponderados propuesto en 4.3.2 mediante rutinas de estimación no lineal de los paquetes SPSS o STATISTICA se obtienen resultados similares al anterior El modelo elegido para ajustar es el esférico: VAR=(c+b*((3/2)*(DIST/a)-(1/2)*(DIST/a)**3))*(DIST
=a) La variable dependiente es VAR, la independiente es DIST y los parámetros a, b y c a optimizar la función de pérdida: ∑(ENES/PRED**2)*(OBS-PRED)**2, donde ENES se refiere a la cantidad de pares que intervienen en el cálculo de cada uno de los seis valores del semivariograma experimental. STATISTICA Método Cuasi - Newton Final loss: .821909587 R=.75514 Variance explained: 57.024% Estimate A=1017.404 B=4.196404 C =3.9132 STATISTICA Método Hooke - Jeeves pattern moves Final loss: 0.821909590 R=.75514 Variance explained: 57.023% Estimate A= 1017.400 B=4.196639 C=3.9130 SPSS Método Sequencial Quadratic Programming Loss function value: 0.8219095869 Estimate A=1017.404 B=4.196403 C=3.913228
En los tres casos el valor estimado del alcance es 1017.4, es decir que en 1017.4m el semivariograma crece desde un efecto pepita de 3.91 hasta alcanzar una meseta de 8.1. Por lo tanto el semivariograma que mejor ajusta, según el criterio de mínimos cuadrados ponderados, es un esférico con las características dadas. En al figura 9.11 se presenta dicho ajuste. C:3
9
C:5
C:4
8
C:6
C:2
7
C:1
SEMIVAR
6 5 4 3 2 1 0 0
500
1000
1500
h
Figura 9.11: Ajuste mínimos cuadrados ponderados de un semivariograma esférico (línea sólida) al semivariograma omnidireccional experimental.
La distribución de los residuos de predicción utilizando este último modelo no difieren del propuesto intuitivamente, como lo muestra el mapa de residuos de predicción de la figura 9.12.
137
Figura 9.12: Mapa de los residuos de predicción cuando se ajusta el modelo esférico con un alcance de 1017.4m, un efecto pepita de 3.91 y una meseta de 8.1.
9.4.3 Kriging. Basado en el último semivariograma ajustado y en el conjunto de los datos de los residuos, las ecuaciones de kriging ordinario (5.3.8) permiten la predicción de los residuos en los puntos de interés. Para dibujar las curvas de nivel de los valores predichos se hace que la posición s vaya recorriendo los nodos de una malla regular con origen en el punto de coordenadas (-1400, -1330), con una separación de 100m entre líneas de la malla que cubre el rectángulo definido por las observaciones. A partir de los valores predichos de los residuos (R) en la malla, se pueden dibujar las curvas de nivel.
Figura 9.13: Mapa que presenta las curvas de nivel de predicción de la variable residuos cuando se realiza un kriging con el modelo esférico con un alcance de 1017.4m, un efecto pepita de 3.91 y una meseta de 8.1.
138
En forma similar que se obtienen las curvas de nivel de las predicciones para la variable se obtienen simultáneamente las curvas de nivel de las estimaciones de las desviaciones estándar de predicción. Así, de esta manera se permite visualizar la incertidumbre de la predicción.
Figura 9.14: Mapa que presenta las curvas de nivel de las desviaciones estándar de predicción de la variable residuos.
9.4 Kriging Mediana Polish. Desde las ecuaciones de interpolación y extrapolación (6.4.11), la estimación mediana polish µ~(s 0 ) puede ser definida para todo s0 del plano, y además, el kriging de los residuos brinda Rˆ (s 0 ) para cualquier posición del plano entonces la predicción ) mediana polish de Z (s ) está dada por: Z~ (s ) = µ~ (s ) + R(s ) 0
0
0
0
En particular se estima la tendencia mediana polish en los puntos de la grilla definida en el kriging; es decir en la malla de puntos cuyo origen es (-1400, -1330) con una separación de 100m entre líneas. Así en los 624 puntos de la malla se realiza la suma de los valores de la tendencia y la de los residuos predichos. Entre dichos puntos se encuentran los de interés para los investigadores, los valores de sus posiciones, dióxido de nitrógeno predichos y sus desviaciones estándares se presentan en la tabla 9.5. oeste norte concentraciones de nitrógeno predichas -1200 -830 0.36 -400 -30 7.08 400 870 4.24 800 -830 3.97
desviaciones estándar de la predicción 2.76 2.58 2.58 2.72
Tabla 9.5: Valores predichos de la concentración de dióxido de nitrógeno y sus desviaciones estándares en los puntos de interés para los investigadores.
139
Bajo el supuesto que Z: "concentración de dióxido de Nitrógeno" es normal, los intervalos de predicción del 95% para cada punto de interés serían: oeste norte concentraciones de nitrógeno predichas -1200 -830 (-5.05, 5.77) -400 -30 (2.02, 12.14) 400 870 (0, 9.30) 800 -830 (0, 9.30) Tabla 9.6: Intervalos de predicción de la concentración de dióxido de nitrógeno en los puntos de interés para los investigadores.
En la figura 9.15 se presenta un mapa con las curvas de nivel de las predicciones mediana polish de las concentraciones de dióxido de nitrógeno. El mapa de los errores estándares de la predicción correspondiente se presentó en la figura 9.14.
Figura 9.15: Mapa que presenta las curvas de nivel de las predicciones mediana polish de la variable concentración de dióxido de nitrógeno.
A modo de conclusión. Como en todo estudio geoestadístico, el análisis estructural tuvo un peso importante en su desarrollo. Para el problema planteado se decidió ante la presencia de tendencia modelarla mediante el procedimiento de mediana polish, y de esa manera trabajar para la predicción con el kriging correspondiente. Otros modelistas probablemente elegirían otro forma para la tendencia, y se producirían diferentes resultados. Pero de acuerdo a los resultados teóricos con respecto al sesgo y ante la presencia de muy poca información, se inclinó por este procedimiento. La elección del ajuste del semivariograma mediante mínimos cuadrados ponderados también se fundamenta en los resultados teóricos que fundamentaron el párrafo anterior. La cantidad de puntos de observación tanto como sus ubicaciones son en general adecuadas debido a los altos costos de la implementación de los sistemas de detección y
140
medición de la concentración de dióxido de nitrógeno, pero se considera que es necesario agregar otros en el acceso de camiones de carga a la ciudad en la zona Este. Se sugiere, para próximos estudios, que los datos estén distribuidos lo más uniformemente posible, de esa manera los pesos que se asignan a los datos no estén tan influenciados por la posición de los puntos de observación a los efectos de mejorar las predicciones. Por último, el mapa de los valores predichos de la concentración de dióxido de nitrógeno junto con el de los desvios estándar del kriging son el resultado de este trabajo, así el procedimiento geoestadístico permitió dar respuesta a la inquietud de los investigadores.
141
CAPÍTULO 10: CONCLUSIONES. Como bien lo dijo Matheron la Geostadística es una mezcla de matemática, estadística y computación, esta tesis es una muestra de esta afirmación. A lo largo de los primeros capítulos de este trabajo se presentaron las bases del enfoque Geoestadístico en una forma más ordenada y ampliada. Esto es el fruto de la conjunción de la lectura de la bibliografía disponible con las informaciones de distintas páginas de Internet acerca de la temática. Para realizar las aplicaciones, se debe recurrir a paquetes de computación que están compuestos por distintos programas que abarcan las diferentes etapas de un estudio geoestadístico. Se trabajó especialmente con el paquete GEOEAS (Englund, E. & Sparks, A..1991) muy recomendado por modelistas en distintas páginas visitadas de Internet. Existen otros software libres al igual que el GEOEAS, quizás con un poco más de herramientas como el BLUEPACK, pero la plataforma y el manual del usuario del GEOEAS hacen sencilla la tarea de cualquier persona que se inicie en el estudio de esta metodología. La desventaja del paquete GEOEAS es su entorno DOS. A los efectos de comparación también se utilizó el programa VARIOWIN, que como su nombre lo indica funciona bajo WINDOWS, pero sólo esta diseñado para tratar las etapas de estimación y modelado del variograma. Las herramientas descriptivas clásicas fueron trabajadas con software de estadística general como ser SPPS y STATISTICA. Las otras, es decir las propias de la Geoestadística, fueron proporcionadas por el GEOEAS o por el VARIOWIN. También se programó en lenguaje MATHEMATICA algunos instrumentos descriptivos no disponibles en los softwares antes mencionados. Se usó además el software MATHEMATICA para comprobar las virtudes del kriging según la disposición de los datos. Si bien el programa SURFER no es específico de Geoestadística se lo usó en alguna aplicación para la generación de los mapas de las curvas de "iso- contenidos". Es importante destacar que para realizar mapas de las predicciones de la variable de interés mediante kriging (una de las alternativas que presenta el programa SURFER) se debe introducir el tipo de variograma adecuado y las características que lo definen, no se debe usar los valores que presenta por defecto el programa porque no se estaría representando en forma adecuada. Es necesario un estudio geoestadístico previo que puede ser llevado a cabo con el GEOEAS. En versiones más actuales de SURFER se pueden realizar algunas etapas de dicho estudio. Se consideró importante en cada una de las aplicaciones realizar el análisis estructural poniendo énfasis en la estimación del variograma, porque es la herramienta fundamental en el kriging de la variable de interés. Se trabajó con datos provistos por las páginas de Internet a los efectos de ganar confianza en el tratamiento de este tipo de información comparando con los resultados o propuestas de otros investigadores, uno de dichos estudios es el presentado en el capitulo 7. En las distintas etapas de un estudio geoestadístico es necesario la iteracción con los especialistas proveedores de los datos de esa manera los resultados y la interpretación de los mismos son más provechosas, esto ocurrió fundamentalmente con la gente de Hidrogeología cuando se realizó el estudio del acuífero de la Caldera. En el estudio de las concentraciones de dióxido de nitrógeno sólo se logró una reunión con los investigadores y la provisión de los datos.
142
Los estudios realizados para completar esta tesis, han permitido detectar varios temas directamente relacionados con el enfoque propuesto que serían muy interesante analizar en el futuro, entre ellos los siguientes: a) Generación de curvas de nivel. b) Simulación Geoestádistica. c) Muestreo espacial. d) Métodos Geoestadísticos multivariantes Generación de curvas de nivel. Como los resultados geoestadísticos generalmente se presentan en mapas bidimensionales a través de las curvas de "iso contenido" es interesante investigar como se produce la interpolación para dibujar dichas curvas que puede ser tipo Splain como lo realiza el Geoeas, o distintas alternativas que presentan otros programas como el Statistica. La interpolación fractal sería otra alternativa interesante. Vale aclarar que para esta tesis los mapas con curvas de niveles fueron generados según la propuesta del programa Geoeas. Simulación Geoestadística. El fin último de la Geoestadística es la caracterización del fenómeno, lo que conduce a varios tipos de aplicaciones. La primera de ella es la predicción (estimación), la que fue tratada en esta tesís. La predicción suele producir mapas que son mucho más “suaves” que la realidad. Por ello, en los casos en que la variabilidad espacial sea de interés es necesario recurrir a técnicas de simulación, segundo grupo de aplicaciones, a fin de obtener realizaciones plausibles de la variable estudiada. Muestreo espacial Otro tipo de aplicación de este enfoque que puede ser estudiado, es el que resulta de que la Geoestadística permite obtener no sólo la predicción sino también una medida de la incertidumbre asociada a ella. Así la Geoestadística constituye el marco ideal para seleccionar la ubicación de puntos de muestreo de forma que se minimice la incertidumbre de la predicción. Métodos Geoestadísticos multivariantes. En esta tesis se puso énfasis en el tratamiento de datos considerados como una realización de un proceso aleatorio univariado, pero se puede estudiar el caso de datos que sean una realización de un proceso aleatorio multivariado. En estos procesos además de la correlación espacial cabe esperar la correlación entre las variables.
Se destacan distintas líneas de estudio relacionadas con el enfoque propuesto en esta tesis pero por supuesto que los otros enfoques para tratar datos espaciales como el Método Lattice, propuesto en Cressie (1991), son importantes de explorar.
143
CAPITULO 11: BIBLIOGRAFIA • Candela,L., Olea, R. & Custodio, E..1988. “Lognormal kriging for the assessment of reliability in groundwater quality control observation networks”. Journal of Hydrology 103, 67-84. • Cox Lawrence H. & Piegorsch Walter W..1996. “Combining environmental information. I: Environmental monitoring, measurement and assessment”. Enviorometrics 7, 299-308. • Cressie Noel A.C..1991. Statistics for Spatial Data. John Wiley & Sons, Inc. •
Dirección General de Obras Sanitarias, 1995. Legajos de Perforaciones efectuadas por la D.G.O.S. y por empresas particulares. Inédita. Salta.
• Englund, E. & Sparks, A..1991. Geo- EAS 1.2.1 User’s Guide, US-EPA Report #600/8-91/008, EPA-EMSL, Las Vegas, NV. • Erickson B.H., Nosanchuk T.A.. 1992. “Understanding Data”. Open University Press. Buckingham. •
Foster, S. y Hirata, R. 1991. Determinación del riesgo de contaminación de aguas subterráneas. CEPIS (Centro de Planeamiento De ingeniería Sanitaria y Ciencias del Medio Ambiente). Organización mundial de la Salud; Organización Panamericana de la Salud (Programa de Salud Ambiental, HPE), Lima, Perú.
•
Fuertes, A. et.al 1999. Plan de trabajo del proyecto de investigación N0 577: "Hidrogeología del Sistema Acuifero La Caldera"
•
Fuertes, A. et.al 1999."Informe de avance del proyecto de investigación N0 577:"Hidrogeología del Sistema Acuifero La Caldera"
• Gantmacher F.R. 1989.“The theory of matrices”. Volume one. Chelsea Publishing Company. • Hoaglin D. C., Mosteller F., Tukey J. W.. 1983. “Understanding Robust and Exploratory Data Analysis”. John Wiley & Sons, Inc • Isaaks E. H., Srivastava R.M.,1989. “An Introduction to Applied Geostatistics”. Oxford University Press. • Lea Cox B. & Wang, J.S. Y. 1993. "Fractal surfaces: measurement and applications in the earth sciences". Fractals, Vol. 1, 87- 115. • Manual de SPSS for Windows. Release 10.01. Standard Version. • Manual de Statistica '99 Edition. Kernel Release 5.5
144
• Musso de Dip, H, E. et al 1999."Informe Final del proyecto de investigación N0 577: Determinación de dióxido de Nitrógeno en la atmósfera de Salta (Capital)" Consejo de Investigación de la Universidad de Salta. • Notas del curso Applied Geostatistical Workshop, dictado por los Drs. David B. Mark, Carol A. Gotway y Gary W. Hergert de la University of Nebraska-Lincoln, en Santa Marta Colombia. • Notas del curso Geoestadística. 1997, dictado por el Dr. Javier Sánchez Vila, en Santa Rosa de la Pampa. Argentina. • Oliver M.A. & Webster R. 1991.“How geostatistics can help you”. Soil, use and management 7, 216-217. • Oliver M.A.& Webster R. & McGrath S. P.1996. “Disjuntive kriging for environmental management”. Enviorometrics 7, 333-357. • Pannatier, Y.. 1996. "Variowin,Software for Spatial Data". Springer Verlag. • Peraudin, J. J. et al. 1998. "La Feulle, The Universe of Geostatistics" Geovariances Newsletter,#10. • Ripley, Brian D.. 1981. "Spatial Statistics". John Wiley & Sons, Inc. • Samper Calvete, F.J. y Carrera Ramírez, J.. 1996. "Geoestadística, aplicaciones a la hidrología subterránea". Centro Internacional de Métodos Numéricos en Ingeniería. Universitat Politécnica de Catalunya. Barcelona. • Stein M. L., 1999. “Interpolation of Spatial Data: Some Theory for Kriging”. Springer Verlag. • Wolfram S., 1991. “The Mathematica book” — 3rd ed. Mathematica Version 3. Sitios de Internet visitados en búsqueda de información: •
http://www.geoavariances.fr
•
http://www.curie.ei.jrc.it/ai-gostats.htm
•
http://www-sst.unil.ch/geostatistics.html
•
http://www.u.arizona.edu/~donaldm/
•
http://www.springer-ny.com/supplements/variowin.html
•
http://www.statsoft.com
•
http://www.goldensoftware.com/products/surfer/surfer.shtml
145
APÉNDICE C5 1. Obtención de las ecuaciones del Kriging simple. Para la pérdida dada por el error cuadrático, el mejor predictor es E (Z (s 0 ) / Z ) , el cual no siempre es lineal en Z. En vez de preguntar por el mejor predictor, uno podría preguntarse por el mejor predictor lineal, esto es; obtener l1,l2 , ... , ln , k en n
p(Z,s0) =
∑
li Z(si) + k,
i =1
tal que minimice E(Z(s0) − p(Z,s0))2. Es decir, se debe minimizar 2
n E Z (s 0 ) − ∑ li Z (s i ) − k , con respecto a l1,l2, ... ,ln, k . i =1
Se puede expresar: 2
n n n E Z (s 0 ) − ∑ li Z (s i ) − k = var Z (s 0 ) − ∑ li Z (s i ) + µ (s 0 ) − ∑ li Z (s i ) − k i =1 i =1 i =1 donde µ(s) =E(Z(s)) s∈ D.
2
Al ser suma de dos expresiones positivas o nulas , el valor mínimo de: n E Z ( s 0 ) − ∑ li Z ( s i ) − k i =1 es la suma de los valores mínimos de ambas expresiones.
2
n
Si se elige k 0 = µ (s 0 ) − ∑ li Z (s i ) el segundo sumando alcanza su valor mínimo: 0. i =1 n Minimizar var Z (s 0 ) − ∑ li Z (s i ) respecto a los li i = 1,2, ... ,n equivale minimizar i =1 :
n
var ( Z (s 0 )) + ∑ j =1
(
n
)
n
∑ l j li cov Z (s i ) ; Z (s j ) − 2 ∑ li cov( Z (si ) ; Z (s 0 )) = i =1
i =1
= var ( Z (s 0 )) + l ′ Σ l − 2 l ′ c
donde l’ = ( l1 , l2 , ... , ln )’ , c ≡ (C(s0,s1), ... , C(s0,sn))’ y ∑ es una matriz n x n cuyo elemento (i,j) es C(si,sj). Como: ∂ Var Z (s 0 ) −
∂l
n
∑l
i
i =1
Z (s i )
= 2 ∑ l −2c = 0 por lo tanto ∑ l = c.
146
[
]
∑−1 existe con seguridad si el proceso Z es tal que Cov Z (s i ), Z (s j ) = C (s i − s j ) porque al ser esta función definida positiva, el determinante de la matriz del sistema es no nulo entonces lo’ = c’ ∑−1. Nótese que esta matriz no depende de s por lo que es fácil calcular los coeficientes li para distintos puntos s simplemente cambiando el vector de términos independientes c. Por lo tanto el predictor lineal óptimo p* (Z; s0) [o, más simplemente, Z* (s0)] es: Z*(s0) = p* (Z; s0) = lo’ Z + ko = c’ ∑−1Z + µ(s0) − lo’µ µ= −1 c’ ∑ (Z −µ µ) + µ(s0) donde µ ≡ (µ(s1), ... , µ(sn) )’. El error cuadrático medio de predicción minimizado, a menudo denominado varianza de la predicción es:
σsk2(s0) ≡ E(Z(s0) − p* (Z,s0) )2 = E(Z(s0) − c’ ∑−1( Z −µ µ ) − µ(s0))2 = E(Z(s0) − c’ ∑−1 Z + c’ ∑−1 µ − µ(s0))2 Recordando que: E(X − a)2 = (E(X) − a)2 + Var(X) donde a es una constante. Sea X = Z(s0) − c’ ∑−1 Z y a = c’ ∑−1 µ − µ(s0) E(X) = µ(s0) − c’ ∑−1 µ Var(X) = Var (Z(s0))+ Var (c’ ∑−1 Z) − 2 Cov (Z(s0), c’ ∑−1 Z) = Var (Z(s0) )+ c’ ∑−1 ∑ ∑−1 c − 2 c’ ∑−1 c = Var (Z(s0)) − c’ ∑−1 c Reemplazando en la expresión de σsk2(s0) se obtiene: σsk2(s0) = Var (Z(s0)) − c’ ∑−1 c = C(s0, s0) − c’ ∑−1c .
2. Obtención de las ecuaciones del kriging ordinario. El predictor óptimo p (Z; B ) se obtiene minimizando el error cuadrático medio de predicción: σk2(s0) ≡ E(Z(s0) − p (Z,s0) )2 (1) n
sobre la clase de los predictores lineales
n
∑ λ i Z ( si ) que satisfacen
∑λ
i =1
i
= 1.
i =1
Se debe minimizar: n n E Z (s 0 ) − λ i Z (s i ) − 2m λi − 1 (2) i =1 i =1 con respecto a λ 1 , λ 2 , ..., λ n y m, donde m es el multiplicador de Lagrange asociado
∑
n
con la restricción
∑λ
i
∑
= 1.
i =1 n
La condición
∑λ
i
= 1 permite que la siguiente identidad sea válida:
i =1 n
[Z(s0) − ∑ λ i i =1
Z ( s i ) ]2 = −
1 2
n
n
∑∑λ i λ i =1 j =1
j
(
Z ( s i ) − Z (s j )
147
)
2
+
n
∑λ i =1
i
( Z (s 0 ) − Z (s i )) 2
Aplicando el operador esperanza miembro a miembro: n
E[Z(s0)− ∑ λ i Z ( s i ) ]2= i =1
−
1 2
∑ ∑ λ i λ j E ( Z ( s i ) − Z (s j )) + ∑ λ i E ( Z (s 0 ) − Z (s i )) 2 n
n
2
i =1 j =1
n
i =1
se obtiene: n
E[Z(s0)− ∑ λ i Z ( s i ) ]2 = − i =1
∑ ∑ λ i λ j γ ( s i − s j ) + 2 ∑ λ i γ (s 0 − s i ) n
n
n
i =1 j =1
i =1
Por lo tanto minimizar (2) equivale minimizar: n
n
(
)
− ∑ ∑ λ i λ jγ si − s j + 2 i =1 j =1
n
n
i =1
i =1
∑ λ i γ (s 0 − s i ) − 2m [ ∑ λ i − 1 ]
(3)
Después de derivar (3) con respecto a λ 1 , λ 2 , ..., λ n y m, e igualar los resultados a cero, se obtiene que los pesos óptimos satisfacen las ecuaciones: −
n
∑ j =1
λ j γ ( s i − s j ) + γ (s 0 − s i ) − m = 0
n
∑λ
i
i = 1, 2, ... , n.
=1
i =1
Esto es, el óptimo λ 1 , λ 2 , ..., λ
n
λ 0 = Γ0 −1 γ donde
puede ser obtenido desde: (4)
0
λ 0 ≡ ( λ 1 ,λ 2 , L ,λ n ,m)′ γ 0 ≡ (γ ( s 0 − s 1 ) ,L , γ (s 0 − s n ) ,1 )
γ (s i − s j ) i = 1,..., n Γ0 ≡ 1 i = n +1 0 i =n +1 Γ0 es una matriz simétrica (n+1) x (n+1).
′
j = 1,..., n j = 1,..., n j =n+1
Desde (4) el vector de coeficientes λ*nx1 ≡ ( λ*1 ,K , λ*n ) ′ esta dado por ′ ( 1 − 1′ Γ −1 γ ) −1 * ′ λ nx1 = γ + 1 Γ −1 1′ Γ 1 y
148
(5)
m= −
1 − 1′ Γ −1 γ 1′ Γ −1 1
(6)
′ ′ donde γ ≡ (γ ( s 0 − s 1 ) ,L , γ (s 0 − s n ) ) , 1′ ≡ (1,L ,1) y Γ es una matriz n x n cuyo elemento genérico es γ (si − sj ) . A continuación se demuestran la validez de (5) y (6) . Γ0 se puede expresar como: Γnxn 1nx1 Γ0 = donde Γnxn ≡ γ (s i − s j ) i , j = 1,K , n 0 11′ xn Para encontrar los pesos óptimos, se necesita determinar la inversa de Γ0. Para ello se procederá de la siguiente manera: Γ0 se puede expresar de la siguiente forma: Γnxn − 11′ 0 Γ0 = + 11′ = Ω + 11′ −1 0′ con el tamaño del vector de unos adecuado en cada caso para que se puedan realizar las operaciones matriciales. 0 es un vector de ceros de tamaño n. La inversa de Γ0 es: Γ = ( Ω + 11′) −1 0
−1
Ω −1 11′Ω −1 = Ω − 1 + 1′Ω −1 1 −1
(7)
donde: −1
Ω
−1
( Γ − 11′) −1 Γ − 11′ 0 = = −1 0′ 0′
0 −1
(8)
entonces, 1′Ω −1 1 = 1′( Γ − 11′) −1 1 − 1
(9)
y el escalar del denominador del segundo término en (7) es: 1 + 1′ Ω −1 1 = 1′ ( Γ − 11′) −1 1
(10)
La matriz Ω −1 11′Ω −1 se expresa como:
(Γ − 11′) −1 (11′)( Γ − 11′) −1 Ω 11′Ω = −1′( Γ − 11′) −1 −1
−1
− (Γ − 11′ ) −1 1 1
(11)
La inversa de Γ − 11′ es: (Γ − 11′) −1 = Γ −1 +
Γ −1 11′ Γ −1 1 = Γ −1 + Γ −1 11′ Γ −1 −1 α 1 − 1′ Γ 1
donde α = 1 − 1′ Γ −1 1 Reemplazando (12) en (9) y en (11)
(12) (13)
149
1 + 1 ′ Ω −1 1 =
1− α
(14)
α
1 −1 1 −1 ( Γ + Γ 11′ Γ −1 ) (11′ )( Γ −1 + Γ −1 11′ Γ −1 ) α α Ω −1 11′Ω −1 = 1 −1′( Γ −1 + Γ −1 11′ Γ −1 ) α
− ( Γ −1 +
Γ −1 11′ Γ −1 )1 α 1 1
(15) Donde la matriz n x n (Γ −1 +
1
α
Γ −1 11′ Γ −1 )(11′)( Γ −1 +
1
α
Γ −1 11′ Γ −1 ) se puede expresar
como: ( Γ −1 +
1
α
Γ −1 11′ Γ −1 )(11′)( Γ −1 +
Γ −1 11′ Γ −1 +
2
α
1
α
Γ −1 11′ Γ −1 ) =
Γ −1 1(1′ Γ −1 1)1′ Γ −1 +
1
α
2
Γ −1 1(1′ Γ −1 1)(1′ Γ −1 1)1′ Γ −1 =
1 2 −1 −1 −1 −1 1 + ( 1′ Γ 1) + 2 (1′ Γ 1) Γ 11′ Γ = α α 2
1 ′ Γ −1 1 −1 1 −1 −1 −1 1 + Γ 11′ Γ = 2 Γ 11′ Γ α α Por lo tanto:
1 −1 1 Γ 11′ Γ −1 − (α Γ −1 + Γ −1 11′ Γ −1 )1 2 α α Ω −1 11′Ω −1 = 1 −1 −1 −1 − 1′( α Γ + Γ 11′ Γ ) 1 α
1 Γ −1 11′ Γ −1 Ω 11′Ω α (1 − α ) = −1 1 − 1′ Γ 1 − 1 1′ Γ −1 1− α −1
−1
−
1 Γ −1 1 1− α α 1− α
(16)
Así reemplazando (16) y (8) en (7) se obtiene:
−1 1 −1 Γ + Γ 11′ Γ −1 Γ = α 0′ −1 0
1 Γ −1 11′ Γ −1 0 α (1 − α ) − 1 ′ Γ −1 −1 − 1− α
150
Γ −1 1 − 1 − α = α1 1 − α
1 −1 −1 −1 Γ − 1 − α Γ 11′ Γ = 1 1′ Γ −1 1− α
1 Γ −1 1 1− α 1 − 1 − α
En consecuencia: 1 1 −1 −1 −1 λ* nx1 Γ ( Ι − 1′ Γ −1 1 11′ Γ ) 1′ Γ −1 1 Γ 1 γ nx1 λ0 = = 1 1 1 −1 m 1′ Γ − −1 −1 1′ Γ 1 1′ Γ 1 De donde:
λ* nx1 = Γ −1 γ nx1 + 1
1 − 1′ Γ −1γ nx1 1′ Γ −1 1
y m= −
151
1 − 1′ Γ −1γ nx1 como se quería demostrar. 1′ Γ −1 1
APÉNDICE C6 1. Obtención de las ecuaciones del Kriging Universal en términos de variogramas. En el kriging universal, el predictor lineal insesgado óptimo, se simboliza como p$ ( Z ; s 0 ) y es aquel que minimiza el error cuadrático medio de predicción:
σ 2e = E [ Z (s 0 ) − p( Z ; s 0 )] sobre λ 1,K , λ n sujeto a λ ′ X = x ′ .
2
(1)
El problema de optimización puede ser expresado equivalentemente usando los multiplicadores de Lagrange como; Se debe minimizar:
E [ Z (s 0 ) − λ ′ Z ] − 2 m .( λ ′ X − x ′ ) Con respecto a los vectores λ ′ y m; donde m ≡ ( m0 ,K , m p ) ′ . Es decir, se debe minimizar 2
(2)
2
p +1 n n E Z (s 0 ) − ∑ λ i Z (s i ) − 2 ∑ m j −1 ∑ λ i f j −1 (s i ) − f j −1 (s 0 ) i =1 j =1 i =1
(3)
con respecto a λ 1 ,K , λ n , m0 ,K , m p . Suponiendo que f 0 (s) ≡ 1 lo cuál garantiza que
n
∑λ
i
=1
(4)
i =1
Ahora, se obtiene que: 2
n 2 Z ( s ) − λ i Z ( s i ) = [x ′ β + δ ( s 0 ) − λ ′ Z ] ∑ 0 i =1
= [x ′ β + δ (s 0 ) − λ ′ X β − λ ′ δ ] n = δ (s 0 ) − ∑ λ i δ (s i ) i =1
n
2
(δ (s ) − δ ( s ))
n
= − ∑ ∑ λ iλ
i
j
i =1 j =1
n
2
+ ∑ λ i (δ (s 0 ) − δ (s i ) )
2
j
2 2
(5)
i =1
Suponiendo que: 2 γ (h) = Var ( Z (s + h) − Z (s)) = Var ( δ (s + h) − δ (s) ) = E ( δ (s + h) − δ (s) ) 2
152
(6)
(3) se transforma en: n
n
− ∑ ∑ λ i λ j γ (s i − s j ) + 2 i =1 j =1
n
p +1
i =1
j =1
i =1
n
∑ λ i γ (s 0 − s i ) − 2 ∑ m j −1 ∑ λ i f j −1 (si ) − f j −1 (s 0 ) (7)
Ecuaciones de kriging universal. Derivando con respecto a λ 1,K , λ n , m0 K , m p e igualando a cero, los pesos óptimos son obtenidos desde
λ u = Γu−1 γ u donde o
o
(8)
λ u ≡ ( λ 1 ,K , λ n , m0 ,K , m p ) ′ λ u ≡ ( λ ′ , m ′)
(9)
′
γ u ≡ (γ (s 0 − s1 ),L , γ (s 0 − s n ) ,1, f 1 (s 0 ),L , f p (s 0 )) ′ γ u ≡ ( γ ′ , x ′) ′
(10)
y Γu es una matriz simétrica (n + p + 1) x (n + p + 1) definida como sigue:
i = 1,K , n γ (s i − s j ) Γu ≡ f j −1− n (s i ) i = 1,L , n 0 i = n + 1,L , n + p + 1
j = 1,L , n j = n + 1,L , n + p + 1
(11)
j = n + 1,L , n + p + 1
que se puede expresar en términos de submatrices por: Γ Γu = X ′ donde: • Γ es la matriz n x n , Γ= γ ( s i − s j )
X 0
i , j = 1, 2 ,K , n .
• X es la matriz n x ( p+1) X = f j −1− n (s i )
i = 1,2,K , n j = n + 1,K , n + p + 1
• 0 es la matriz nula (p+1) x( p+1). Para encontrar la matriz inversa de Γu se aplica el algoritmo de Gauss generalizado, procedimiento que se fundamenta en el apéndice B.
I n + p +1 Γu Sea la matriz formada por los bloques, la matriz identidad de orden − I n + p +1 0 n + p +1 n+p+1: I n + p +1 , su opuesta: − I n + p +1 y la matriz nula de orden n+p+1: 0 n + p +1 . Que se expresa como,
153
Γ X′ −In 0 ( p +1) xn
X 0
0 nx ( p +1) I p +1 0 nx ( p +1) 0 p +1
In 0 ( p +1) xn
0 nx ( p +1) − I p +1
0n 0 ( p +1) xn
Aplicando el algoritmo de escalonamiento generalizado de Gauss: En una primera etapa suponiendo la existencia de Γ −1 se obtiene:
Γ 0 nx ( p +1) 0n 0 ( p +1) xn
X − X ′ Γ −1 X
In − X ′ Γ −1
Γ −1 X − I p +1
Γ −1 0 ( p +1) xn
0 nx ( p +1) I p +1 0 nx ( p +1) 0 p +1
Para una segunda etapa se supone la existencia de ( X ′ Γ −1 X ) −1 ;
Γ 0 nx ( p +1) 0n 0 ( p +1) xn
X − X ′ Γ −1 X
In − X ′ Γ −1
0 nx ( p +1) 0 p +1
− Γ −1 X ( X ′ Γ −1 X ) −1 X ′ Γ −1 + Γ −1 ( X ′ Γ −1 X ) −1 X ′ Γ −1
Γ −1 X ( X ′ Γ −1 X ) −1 − ( X ′ Γ −1 X ) −1 0 nx ( p +1) I p +1
(12)
La matriz inversa de Γu es:
− Γ −1 X ( X ′ Γ −1 X ) −1 X ′ Γ −1 + Γ −1 Γ = ( X ′ Γ −1 X ) −1 X ′ Γ −1 −1 u
Γ −1 X ( X ′ Γ −1 X ) −1 − ( X ′ Γ −1 X ) −1
Por lo tanto se obtiene:
′
λ ′ = (γ + X ( X ′ Γ −1 X ) −1 ( x − X ′ Γ −1 γ )) Γ −1 y
(
m ′ = − x − X ′ Γ −1 γ
)
′
(13)
(14)
( X ′ Γ −1 X ) −1
La varianza del kriging se puede expresar como:
σ 2k (s 0 ) = − λ*′ Γ −1 λ* + 2 λ*′ γ = λ*′ γ + ( λ*′ γ − λ*′ Γ −1λ* ) donde la expresión entre paréntesis luego de trabajo algebraico se transforma en m ′ x . Por lo tanto; σ 2k (s 0 ) = λ*′ γ + m ′ x = λ*u′ γ u (15) Además, otra expresión de la varianza del kriging es:
σ 2k (s 0 ) = γ ′ Γ −1 γ − ( x − X ′ Γ −1 γ ) ′ ( X ′ Γ −1 X ) −1 ( x − X ′ Γ −1 γ )
154
(16)
2. Obtención de las ecuaciones del Kriging Universal en términos de covariogramas. Ahora, se establece el supuesto que el proceso Z(⋅)es estacionario de segundo orden. Entonces, el covariograma C (h) = cov( Z (s + h); Z (s)) está bien definido. En consecuencia el primer término de (5.2.9) se transforma en : 2
n n E Z ( s 0 ) − ∑ λ i Z ( s i ) = E δ ( s 0 ) − ∑ λ i δ ( s i ) i =1 i =1
2
n = Var δ (s 0 ) − ∑ λ i δ (s i ) i =1
(17)
= Cov(δ (s 0 ), δ (s 0 )) + λ ′ Σ λ − 2 λ ′ c
(
donde Σ = Cov δ (s i ); δ (s j )
)
i , j = 1,2, ..., n y c = ( Cov (s 0 − s1 ) ,L , Cov (s 0 − s n ) )
′
Sin la suposición adicional f 0 (⋅) = 1 , se debe minimizar: Cov(δ (s 0 ); δ (s 0 )) + λ ′ Σ λ − 2 λ ′ c − 2 λ ′X m − 2 x ′ m
(18)
Las ecuaciones obtenidas por la diferenciación y la igualación a cero se expresan matricialmente como: Σ U λ U = cU
(19)
o
X λ c = O − m x
Σ X ′
(20)
La inversa de Σ U se obtiene en forma similar a la inversa de Γu , así el vector de pesos óptimos y el vector de multiplicadores de Lagrange están dados por:
′
λ ′={c+ X ( X ′Σ −1 X ) −1 (x− X ′Σ −1c)} Σ −1 m ′=(x− X ′Σ −1 c) ′( X ′Σ −1 X ) −1
(21)
En esta situación la varianza del kriging es: n
n
n
σ 2k (s 0 ) = C (0) − 2∑ λ i C (s 0 − s i ) + ∑ ∑ λ i λ j C (s i − s j ) i =1
i =1 j =1
n
p +1
i =1
j =1
= C (0) − ∑ λ i C (s 0 − s i ) + ∑ m j −1 f j −1 (s 0 ) = C (0) − c ′ Σ −1c + (x − X ′ Σ −1 c) ′ ( X ′ Σ −1 X ) −1 (x − X ′ Σ −1 c)
155
(22)
