GEOESTADISTICA
Teoría y Práctica Aplicación en la Caracterización de Yacimientos Bogotá, Febrero 2010
GEOESTADISTICA Aplicación en la Caracterización de Yacimientos
Ramón Giraldo H. PhD Estadística Universidad Nacional de Colombia
2
Contenido Introducción 1. Análisis Exploratorio de Datos Espaciales 1.1 1.2 1.3 1.4
Datos Espaciales Tipos de Variables y Escalas de Medida Medidas Descriptivas Gráficos Exploratorios
2. Definiciones Básicas de Geoestadística 2.1. 2.2. 2.3. 2.4. 2.5.
Definición y Origen de la Geoestadística Variable Regionalizada Estacionariedad Fuerte e Intrínseca Isotropía Ejemplo.
3. Dependencia o Correlación Espacial 3.1. Funciones de Correlación Espacial 3.1.1. Variograma y Semivariograma 3.1.2. Covariograma 3.1.3. Correlograma 3.1.4. Semivariograma experimental 3.2. Modelos Teóricos de Semivarianza 3.3. Ejemplo.
3
4. Predicción Espacial 4.1. 4.2. 4.2. 4.3. 4.4. 4.5.
Predicción Espacial Optima Métodos Kriging Kriging Ordinario Kriging Indicador Cokriging Ordinario Ejemplo.
5. Apéndice 5.1 Álgebra de Matrices 5.2. Conceptos Estadísticos Básicos 5.2 Regresión y Mínimos Cuadrados.
6. Aplicaciones a Datos Reales 6.1. Aplicación en sísmica 6.2 Aplicación con información de pozos
7. Bibliografía.
4
Introducción La necesidad de acudir a herramientas estadísticas para el análisis de datos en todas las áreas del conocimiento, ha hecho que aparezcan con el correr de los años nuevas metodologías que, no obstante se centran en fundamentos probabilísticos comunes, son específicas para cada una de las diversas disciplinas del saber. Algunos ejemplos son, entre otros, minería, geología y geofísica. La gran relevancia que tiene a nivel mundial el tema de la caracterización de reservorios ha hecho que los profesionales en estadística encaminen esfuerzos en el desarrollo de nuevas técnicas apropiadas para el análisis de información enmarcada dentro de este contexto. Dentro de este los métodos geoestadísticos juegan un papel preponderante. La geoestadística permite cuantificar la incertidumbre y especificar la forma en que ésta varía en el espacio-tiempo. Uno de sus campos de aplicación es la caracterización de reservorios, que involucra un conjunto de métodos probabilísticos, cuyo objetivo es definir el modelo más probable de un reservorio, con sus formas de cuerpos, heterogeneidades petrofísicas, geometría estructural y caracterización paleoambiental. Los yacimientos poseen pozos irregularmente distribuidos en función de cómo haya sido la historia de su desarrollo. Cuando una empresa decide llevar adelante una tarea de perforación necesita conocer qué chances va a tener de encontrar crudo y eso implica minimizar las incertezas que se desprenden de la falta de homogeneidad de los cuerpos. De esta forma, las posibilidades de hallar el recurso buscado aumentan o disminuyen según cuáles sean las condiciones de porosidad y permeabilidad, entre otros factores. Ahí es donde entra la geoestadística, por ser una herramienta que permite predecir en un punto qué valor aproximado se va a tener de una determinada propiedad, y qué incertidumbre asociada se tiene a esa predición, que combinada con la geofísica de reservorio permite integrar la información de pozos y el dato sísmico a fin de determinar nuevas locaciones para drenar las zonas saturadas. En el documento se presentan las definiciones y los conceptos teóricos y en el desarrollo del curso se harán aplicaciones con datos geofísicos reales (capítulo 6 no incluido en el docuemnto). Para el seguimiento completo de la teoría descrita se requiere tener conocimientos básicos de álgebra de matrices y de estadística matemática. Sin embargo aquellas personas que estén poco familiarizadas con estos temas, podrán obviar la lectura de algunas secciones en las que se hacen desarrollos teóricos y centrar su atención en la filosofía de los métodos presentados y en las aplicaciones mostradas en cada uno de los capítulos del documento. Un resumen no exhaustivo de conceptos de álgebra lineal y de estadística es hecho al final en el apéndice. Aunque en el texto se cubren a manera de introducción diversos temas geoestadísticos y se hacen aplicaciones de métodos recientes, es necesario acudir a la lectura de artículos científicos y textos formales para lograr un buen dominio de esta metodología.
5
Capitulo Uno Análisis Exploratorio de Datos Espaciales 1.1. Datos Espaciales. Las mediciones de las características de interés en un estudio regionalizado tienen implícitamente asociadas las coordenadas geográficas de los sitios en donde estas fueron tomadas. Generalmente cuando el área de estudio es considerablemente grande se usa un geoposicionador para establecer dichas coordenadas. En otros casos es suficiente con hacer asignaciones según planos cartesianos. Un esquema general de datos georreferenciados (datos espaciales) es el siguiente: Sitio
1 2 3 4 . . . n
Latitud Norte
Longitud Este
Var. 1
Var. 2
.
X11 X21 X31 X41 . . . Xn1
X12 X22 X32 X42 . . . Xn2
. . . . . . . .
.
.
. . . . . . . .
Var. p
. . . . . . . .
X1p X2p X3p X4p . . . Xnp
En la tabla anterior n es el número de sitios muestreados y p el de variables medidas en cada uno de ellos. Cada Xij corresponde a una medida de una variable (variable j) que puede ser numérica (discreta o continua) o categórica (ordinal o nominal). En general la metodología geoestadística trabaja con datos correspondientes a variables numéricas. Algunas de las variables pueden estar más intensamente muestreadas que las otras (Xij faltantes). Las coordenadas pueden ser planas, geográficas (grados, minutos y segundos) o cartesianas. Sin embargo la posible utilización de unas u otras depende del software empleado para los análisis. 1.2. Tipos de Variables y Escalas de Medida Existen distintas categorías para las variables y se han propuesto numerosas clasificaciones para expresar su variabilidad. Las dos más comunes son: 1) De doble estado o binarias, que son aquellos que pueden tomar sólo dos valores, p.ej. los datos de presencia o ausencia de un mineral. 2) De Multiestado, que son aquellas en que las medidas pueden tomar tres o
6
más valores, éstas pueden ser Cualitativas o Cuantitativas, y presentan diferentes escalas de medida (Digby & Kempton, 1992). 1.2.1 Variables Cualitativas Son aquellos que expresan cualidades no mensurables y se dividen en cualitativas sin secuencia lógica y cualitativas con secuencia lógica dependiendo de su escala de medida. Las posibles escalas de medida de las variables cualitativas son: a) Escala Nominal: Se presenta cuando las observaciones de la variable no pueden ser ordenadas en una secuencia de grados del atributo. Esta es la escala más simple de medida, p. ej. épocas climáticas (seco, lluvioso), sitios geográficos, etc. b) Escala Ordinal: Se presenta cuando las mediciones pueden ser ordenadas de menor a mayor o viceversa, pero las distancias entre los elementos ordenados no tienen ningún sentido físico y si lo tienen, no son iguales a todo lo largo de la escala, p.ej. tipo de grano (arcilla, limo), dureza de un mineral . 1.2.2. Variables Cuantitativas Expresan magnitudes o cantidades, que son el resultado de mediciones de algún instrumento, conteos de eventos o de operaciones matemáticas simples. Estos pueden ser: 1.2.2.1. Variables Discretas: Son aquellas que representan cantidades expresables sólo por un número finito de valores en la escala real, generalmente las que sólo pueden tomar valores enteros, sin fracciones. 1.2.2.2. Variables Continuas: Son aquellas en los que existe potencialmente un número infinito de valores entre dos puntos de la escala. Pueden ser datos enteros o fraccionarios, p. ej. Características físicas o químicas. 1.2.2.3. Variables Derivadas: Son aquellas en que los datos son generados a partir de cálculos simples entre medidas de variables cuantitativas o cualitativas, p. ej. índices, tasas, proporciones, etc. Las variables continuas tienen dos posibles escalas de medida, estas son : a) Escala de Intervalo: Es una escala ordinal en donde las distancias tienen un sentido físico igual a todo lo largo de la escala, pero el punto de valor cero es fijado arbitrariamente, p. ej. el tiempo (el tiempo inicial (t0) puede ser cualquier momento), la altitud (cero se refiere al nivel del mar) y la temperatura en grados Celsius en la que por ejemplo el valor de cero grados no indica ausencia de calor o de agitación de moléculas y por consiguiente no es posible afirmar que un cuerpo de 20 grados tiene el doble de calor que uno de 10. b) Escala de Razón: Es una escala de intervalo en la que no hay que fijar un cero arbitrario, entonces el resultado de dividir o multiplicar un valor de la escala por otro tiene un sentido físico, p. ej. Las variables químicas.
7
1.3.Medidas Descriptivas. Siempre que se va a realizar un análisis estadístico, es conveniente realizar un estudio exploratorio de los datos. Esto implica establecer si los datos están muy agrupados o muy dispersos, cual es el punto representativo de la agrupación y si hay observaciones muy alejadas de las restantes. Estos aspectos se tratan a continuación: 1.3.1. Medidas de Localización Estas medidas indican alrededor de que valor se agrupan los datos, generando valores representativos de las observaciones (tabla 1). a) Media: Es el promedio aritmético de las observaciones y una medida representativa cuando no hay valores muy extremos en los datos, porque en esos casos es afectada por ellos, corriéndose hacia un lado de la distribución. b) Mediana: Se define como el valor de la variable que supera el 50 % de las observaciones y es superado por el restante 50 %. Esta medida tiene en cuenta sólo el orden de los datos más no su magnitud, por esto no se deja afectar por los valores atípicos (extremos) y puede ser más representativa que la media en muchos casos. c) Cuantilas: Particionan en intervalos de igual amplitud la distribución. Particularmente los cuartiles (Qi , i=1,2,3) son de gran utilidad, como se verá más adelante, en la detección de observaciones atípicas
Tabla 1. Medidas de localización con sus respectivos cálculos muestrales Medida
Cálculo.
Media
n
x Mediana
xi i 1
n
Estadísticas de orden:
x( 0) x(1) x( 2 ) ... x( n) x x( n 1 ) / 2 Si n impar ~
Si n par
xn / 2 x( n / 2 1 ) ~ x 2
Similar a la mediana pero se divide Cuantilas (Cuartiles, Deciles, sobre 4 en el caso de cuartíles o sobre 10 en los Deciles. etc)
8
1.3.2. Medidas de Variabilidad Indican cuanto se alejan o dispersan los datos con respecto a las medidas de localización o el grado de homogeneidad de los mismos (tabla 2). a) Varianza: Es una medida de la dispersión en la distribución de probabilidad de una variable aleatoria, expresada en unidades cuadradas. b) Desviación Estándar: Indica en promedio cuánto se alejan las observaciones de la media aritmética; está dada en las mismas unidades de la variable, a diferencia de la varianza. c) Coeficiente de Variación (C.V.): Es una medida relativa de variabilidad y, en general, se acepta que un conjunto de datos es relativamente homogéneo si el C.V. es menor del 30%, aunque algunos autores refutan este concepto. d) Rango y Rango Intercuartílico: Representan el recorrido de la variable y la distancia entre los cuartíles, respectivamente. Son útiles cuando se comparan dos o más distribuciones. Tabla 2. Medidas de variabilidad con sus respectivos cálculos muestrales Medida
Cálculo
Varianza
n
S2 Desviación Estándar Error Estándar
xi X i 1
2
n1
S S2 E .E
S2 n
Coeficiente de Variación Rango Rango entre Cuartíles
S ~x Xmáx - Xmin
Q3 Q1
1.4. Gráficos Exploratorios A continuación se presentan algunos gráficos que resumen la información de un conjunto de datos, indicando aparte de las medidas de localización y variabilidad, aspectos importantes como la detección de observaciones atípicas.
9
a). Histogramas Un histograma es un conjunto de rectángulos, cada uno de los cuales representa un intervalo de agrupación o clase. La base de cada rectángulo es igual a la amplitud del intervalo, y la altura es proporcional a la frecuencia (absoluta o relativa) de cada clase. Para obtener una buena representación de las frecuencias de las observaciones, se recomienda calcular 2 n como el número aproximado de intervalos, sin embargo se debe tener en cuenta que esto depende de la variable estudiada y no es una regla que se debe seguir siempre (Peña, 1987). b). Diagrama de Caja Este diagrama resume un conjunto de observaciones univariadas, suministrando un análisis exploratorio de los datos, útil para estudiar simetría, supuestos distribucionales y detectar observaciones atípicas (Hoaglin et al., 1983). El gráfico (Fig. 1) divide los datos en cuatro áreas de igual frecuencia. La caja central encierra el 50%, tomando la línea vertical como la mediana. La línea horizontal va desde el primer cuartíl hasta menos 1.5 veces el rango intercuartílico del primer cuartíl, y desde el tercer cuartíl hasta mas 1.5 veces el rango intercuartílico del tercer cuartíl. Los puntos que estén por fuera de la línea horizontal se consideran puntos afuera y en algunos casos cuando están a más de tres veces del rango entre cuartiles se consideran como puntos muy afuera o muy alejados
Q1 - 1.5R.I
Q1
Me
Q3
Q3 + 1.5 R.I
Figura 1. Representación de un diagrama de caja. c). Diagrama de Tallos y Hojas El diagrama de tallos y hojas de Tukey es una forma semi-gráfica de presentar la información para variables cuantitativas, especialmente cuando el número total de datos es pequeño (menor que 50). Los principios para construirlo son: Redondear los datos a dos o tres cifras significativas, expresándolos en unidades convenientes. Disponerlos en una tabla con dos columnas separadas con una línea como sigue: Para datos con dos dígitos, escribir a la izquierda de la línea los dígitos de las decenas, formando así el tallo, y a la derecha las unidades que serán las ramas. Por ejemplo 87 se escribirá 87.
10
Para datos con tres dígitos el tallo estará formado por los dígitos de las decenas y centenas, que se escribirán a la izquierda, separados de las unidades. Por ejemplo, 127 será 127. Cada tallo define una clase y sólo se escribe una vez. El número de hojas representa la frecuencia de dicha clase (Figura a) de abajo). En algunas ocasiones, hay muchas observaciones en cada fila y conviene abrir cada tallo en dos, y en algunas otras es posible abrir cada fila en 5 clases (Figura b) de abajo). Las observaciones afuera y muy afuera del diagrama de caja son indicadas en el diagrama de caja con las expresiones bajo y alto (ver figura 4, sección 1.6).
a) 11 12 13 14 15 16
34 24577 345 27 2 1
b) 11* 11º 12* 12º 13* 13º 14*
0011223333444 55567777788888999 1112333344 55567789 00123 6678 22
d). Gráfico de Datos Clasificados en Intervalos Una clasificación de los datos en intervalos de clase, definidos por símbolos, dentro del área de estudio, es útil en la identificación de posibles tendencias en los valores de la variable, de zonas con mayor o menor magnitud o de observaciones extremas (ver figura 5 en la sección 1.6) 1.5. Relación entre variables a). Covarianza y Correlación. La covarianza mide la variabilidad conjunta de dos variables. Es una extensión de la varianza al caso bidimensional: n
( X i X )( Yi Y )
COV ( X ,Y ) i 1
n
El coeficiente de correlación mide el grado de asociación lineal que existe entre dos variables X y Y. Se calcula mediante: n
r
COV ( X ,Y ) Sx S y
( xi x )( yi y ) i 1
n
( xi x )2 i 1
n
( yi y )2 i 1
11
El coeficiente de correlación es un número en el intervalo [-1, 1]. Un valor de r = -1 indica una relación lineal negativa perfecta entre X y Y , mientras que una valor de r = 1 señalará una asociación positiva perfecta de X y Y. Si r = 0, entonces se concluirá que no existe ninguna relación lineal entre X y Y. b) Gráficos de Dispersión. Son muy útiles tanto para la detección de relaciones entre las variables como para la identificación de tendencias en el valor promedio de la variable en la región (relación entre la variable medida y las coordenadas geográficas). Un supuesto fundamental en el análisis geoestadístico es que el fenómeno sea estacionario, para lo cual, entre otros aspectos, el nivel promedio de la variable de estudio debe ser constante en todos los puntos del área de estudio. Una detección de tendencia en el gráfico de dispersión puede ser una muestra de que no se satisface dicho supuesto. El gráfico se construye tomando como eje de las abcisas la variable que representa la coordenada geográfica (latitud o longitud) y en el eje de las ordenadas la variable cuantitativa de estudio. La observación de la nube de puntos resultante o incluso el ajuste de una línea de regresión (Fox, 1984), permiten establecer si existe dicha tendencia.
12
Capitulo Dos Definiciones Básicas de Geoestadística 2.1. Origen y Definición de Geoestadística Los orígenes de la geoestadística se encuentran en el campo de la minería. Como antecedentes suelen mencionarse trabajos de Sichel (1947; 1949) (citado en Samper & Carrera, 1990) y Krige (1951). El primero observó la naturaleza asimétrica de la distribución del contenido de oro en las minas surafricanas, la equiparó a una distribución de probabilidad lognormal y desarrolló las fórmulas básicas para esta distribución. Ello permitió una primera estimación de las reservas, pero bajo el supuesto de que las mediciones eran independientes, en clara contradicción con la experiencia de que existen “zonas” más ricas que otras. Una primera aproximación a la solución de este problema fue dada por geólogo G. Krige que propuso una variante del método de medias móviles, el cual puede considerarse como el equivalente al krigeado simple que, como se verá más adelante, es uno de los métodos de estimación lineal en el espacio con mayores cualidades teóricas. La formulación rigurosa y la solución al problema de estimación vino de la mano de Matheron (1962). En los años sucesivos la teoría se fue depurando, ampliando su campo de validez y reduciendo las hipótesis necesarias (Samper & Carrera, 1990). De la minería las técnicas geoestadísticas, se han "exportado" a muchos otros campos como hidrología, física del suelo, ciencias de la tierra y más recientemente al monitoreo ambiental y al procesamiento de imágenes de satélite. La geoestadística es una rama de la estadística que trata fenómenos espaciales (Journel & Huijbregts, 1978). Su interés primordial es la estimación, predicción y simulación de dichos fenómenos (Myers, 1987). Esta herramienta ofrece una manera de describir la continuidad espacial, que es un rasgo distintivo esencial de muchos fenómenos naturales, y proporciona adaptaciones de las técnicas clásicas de regresión para tomar ventajas de esta continuidad (Isaaks & Srivastava, 1989). Petitgas (1996), la define como una aplicación de la teoría de probabilidades a la estimación estadística de variables espaciales. La modelación espacial es la adición más reciente a la literatura estadística. Geología, ciencias del suelo, agronomía, ingeniería forestal, astronomía, o cualquier disciplina que trabaja con datos colectados en diferentes locaciones espaciales necesita desarrollar modelos que indiquen cuando hay dependencia entre las medidas de los diferentes sitios. Usualmente dicha modelación concierne con la predicción espacial, pero hay otras áreas importantes como la simulación, el diseño muestral y los modelos en enmallados (lattices) (Cressie, 1989). Cuando el objetivo es hacer predicción, la geoestadística opera básicamente en dos etapas. La primera es el análisis estructural, en la cual se describe la correlación entre puntos en el espacio. En la segunda fase se hace predicción en sitios de la región no muestreados por medio de la técnica kriging (capitulo 4). Este es un proceso que calcula un promedio ponderado de las observaciones muestrales. Los pesos asignados a los valores muestrales 13
son apropiadamente determinados por la estructura espacial de correlación establecida en la primera etapa y por la configuración de muestreo (Petitgas, 1996). Los fundamentos básicos de estas etapas son presentados a continuación. Se realiza también una revisión del caso en que se miden simultáneamente varias variables en cada sitio de muestreo y se desea hacer predicción de una de ellas con base en información de las otras. En este caso la técnica de predicción es conocida como cokriging (capitulo 4). Algunos temas especiales como el diseño de una red de muestreo óptima, en términos de varianza de predicción y costos, y el análisis de componentes principales regionalizado también serán estudiados (capitulo 5). 2.2. Variable Regionalizada. Una variable distribuida en el espacio de forma que presente una estructura de correlación, se dice que es una variable regionalizada. De manera más formal se puede definir como un proceso estocástico con dominio contenido en un espacio euclidiano m-dimensional Rm, {Z(x) : x D Rm }. Si m = 2, Z(x) puede asociarse a una variable medida en un punto x del plano (Díaz-Francés, 1993). En términos prácticos Z(x) puede verse como una medición de una variable aleatoria (por ejemplo concentración de un contaminante) en un punto x de una región de estudio. Recuérdese que un proceso estocástico es una colección de variables aleatorias indexadas; esto es, para cada x en el conjunto de índices D, Z(x) es una variable aleatoria. En el caso de que las mediciones sean hechas en una superficie (sólo se tengan longitud y latitud como coordenadas), entonces Z(x) puede interpretarse como la variable aleatoria asociada a ese punto del plano (x representa las coordenadas, planas o geográficas, y Z la medición de la variable en cada una de ellas). Estas variables aleatorias pueden representar la magnitud de una variable ambiental medida en un conjunto de coordenadas de la región de estudio. 2.3. Estacionariedad 2.3.1. Estacionariedad de Segundo Orden Sea {Z(x) : x D Rm } una variable regionalizada definida en un dominio D contenido en Rm (generalmente una variable medida en la superficie de una región) se dice que Z(x) es estacionario de segundo orden si cumple (Díaz-Francés, 1993): a.
E [ Z(x)] = k, k R, x D Rm . El valor esperado de la variable aleatoria es finito y es una constante para todo punto en el dominio (el valor promedio es igual en todo punto de la región).
b. COV [ Z(x) , Z(x+h)] = C(h) < Z (x) tiene covarianza finita y es una función única del vector de separación h entre cada pareja de puntos. 14
Obviando la dirección de variación, el supuesto de estacionariedad en la media puede ser estudiado a través de un gráfico del promedio en función de la distancia (Fig. 5).
Valor Promedio.
19 18 17 16 0
10000 20000 Distancia (m)
30000
Figura 5. Gráfico de dispersión de valores promedios de una variable simulada en función de la distancia entre puntos de muestreo La figura anterior fue elaborada con base en simulación. Se generaron datos de una variable hipotética con valores uniformemente distribuidos entre 16 y 19. Luego estos fueron asignados aleatoriamente a las coordenadas dadas en la matriz del anexo. En el proceso de elaboración de la figura se tomaron vecindades, se calculó la media y la distancia euclidiana promedio entre las correspondientes coordenadas dentro de la vecindad. El gráfico corresponde a la nube de puntos de las dos variables. Se puede establecer que el valor promedio está fluctuando alrededor de un valor entre 17 y 18 y que por consiguiente la media no tiene ninguna tendencia de cambio en función de la distancia existente entre los puntos de muestreo. 2.3.2. Estacionariedad Débil Generalmente se trabaja sólo con la hipótesis que pide que los incrementos [Z(x+h)- Z(x)] sean estacionarios, esto es (Clark, 1979): a. Z(x) tiene esperanza finita para todo punto en el dominio. Lo que implica que la esperanza de los incrementos es cero. E [ Z(x+h) - Z(x)] = 0 b. Para cualquier vector h, la varianza del incremento está definida y es una función única de la distancia. V [ Z(x+h) - Z(x)] = 2 (h)
15
El concepto de estacionariedad es muy útil en la modelación de series temporales (Box & Jenkins, 1976). En este contexto es más fácil la identificación, puesto que sólo hay una dirección de variación (el tiempo). En el campo espacial existen múltiples direcciones y por lo tanto se debe asumir que en todas el fenómeno es estacionario. Cuando el nivel promedio de la variable no es el mismo en todas las direcciones o cuando la covarianza o correlación dependan del sentido en que se determinan, no habrá estacionariedad. Si la correlación entre los datos no depende de la dirección en la que esta se calcule se dice que el fenómeno es isotrópico, en caso contrario se hablará de anisotropía. En Journel & Huijbregts (1978), se trata el caso de la anisotropía y se proponen algunas soluciones. Cressie (1986) discute cual debe ser el tratamiento en caso de que la media no sea constante.
16
Capitulo Tres Dependencia o Correlación Espacial 3.1. Funciones de Correlación Espacial La primera etapa en el desarrollo de un análisis geoestadístico es la determinación de la dependencia espacial entre los datos medidos de una variable. Esta etapa es también conocida como análisis estructural. Para llevarla a cabo, con base en la información muestral, se usan tres funciones: El semivariograma, el covariograma y el correlograma experimental. A continuación se hace una revisión de los conceptos asociados a cada una de ellas y se describen sus bondades y limitaciones. 3.1.1. Variograma y Semivariograma. Representa la varianza de los incrementos de la variable regionalizada y se denota por 2(h). De acuerdo con lo anterior utilizando la definición teórica de la varianza en términos del valor esperado de una variable aleatoria, tenemos:
2 ( h) V Z( x h) Z( x) 2 2 E Z( x h) Z( x) E Z( x h) Z( x) 0 2 E Z( x h) Z( x) .
La mitad del variograma ( (h)), se conoce como la función de semivarianza y caracteriza las propiedades de dependencia espacial del proceso. Dada una realización del fenómeno, la función de semivarianza es estimada por medio del semivariograma experimental, que se calcula mediante (Wackernagel, 1995):
( h)
Z( x h) Z( x)
2
2 n ( h)
donde Z (x) es el valor de la variable en un sitio x, Z (x+h) es otro valor muestral separado del anterior por una distancia h y n(h) es el número de parejas que se encuentran separadas por dicha distancia.
17
3.1.2. Covariograma y Correlograma. De acuerdo con la fórmula de la covarianza dada en el capitulo uno, la función de covarianza espacial entre parejas de observaciones que se encuentran a una distancia h está dada, con base en los datos muestrales, por: n
( Z( x h) m)( Z( x) m)
COV Z( x h), Z( x) i 1
n
n
Z( x h) Z( x)
i 1
n
m 2 C( h )
donde m, si el proceso es estacionario de segundo orden, representa el valor promedio en todo punto de la región de estudio y n es el número de parejas de puntos que se encuentran a una distancia h. De otro lado para determinar la fórmula del correlograma, dado que se conoce el covariograma, sólo faltaría establecer las desviaciones estándar en cada punto del dominio. Sin embargo, si el fenómeno es estacionario, éstas al igual que la media son constantes y por consiguiente: r ( h)
COV( Z( x h), Z( x)) C( h) C( h) C( h) Sx h Sx V( Z( x)) C(0) S2x
Cualquiera de las tres funciones de dependencia espacial mencionadas, es decir semivariograma, covariograma o correlograma, puede ser usada en la determinación de la relación espacial entre los datos. Sin embargo como se puede observar en las fórmulas, la única que no requiere que la media del proceso (m) sea conocida, es la función de semivarianza. Por esta razón, fundamentalmente, en la práctica se emplea el semivariograma y no las otras dos funciones. A continuación se presenta un ejemplo ilustrativo del cálculo de la función de semivarianza experimental Ejemplo. Suponga que se tienen medidas sobre una variable hipotética cuyos valores están comprendidos entre 28 y 44 unidades y su configuración en la región de estudio es como se presenta en el esquema de la siguiente página. Como se indica en la representación, la distancia entre cada par de puntos contiguos es de 100 unidades. Luego si existe un punto faltante la distancia entre los dos valores ubicados a cada lado de éste será de 200 unidades. Veamos como calcular bajo esta situación el semivariograma experimental. Por simplicidad se calcularán sólo los semivariogramas en sentido “occidente-oriente” (izquierda-derecha) y “sur-norte”(inferior-superior), debido a que para obtener un semivariograma experimental en el que sólo se tenga en cuenta la distancia y no la orientación, se requeriría calcular la distancia euclidiana entre todas las parejas de puntos.
18
44
40
42
40
39
37
36
42
43
42
39
39
41
40
38
37
35
38
37
37
33
34
35
37
36
36
35
35
34
33
32
29
28
29
30
32
37
37
35
38
36
35
36
38
37
35
30
200
100
En primer lugar en sentido izquierda-derecha se encuentran todas las parejas de puntos que están a una distancia de 100 unidades. Una vez detectados estos puntos se aplica la fórmula del semivariograma experimental. De igual forma se procede para las distancias de 200, 300, 400 y 500 unidades. Específicamente en el caso de las distancias de 100 y 200 unidades se realiza la siguiente operación: (100) = (38 - 37)2 + (37 - 35)2 + (29 - 30)2 + ... + (37 - 36)2 /2* 36 = 1.458 (200) = (40 - 44)2 + (40 - 40)2 + (42 - 39)2 + ... + (29 - 32)2 /2* 36 = 3.303 Similarmente procedemos para las otras distancias y para el sentido inferior-superior. Los resultados se muestran en la siguiente tabla. Tabla 5. Valores de la función de semivarianza experimental en dos direcciones para el conjunto de datos hipotéticos de la configuración de datos dada arriba. Distancia 100 200 300 400
Semivarianza Sentido Este - Oeste 1.45 3.30 4.31 6.69
Semivarianza Sentido Norte - Sur 5.34 9.87 18.88 27.53
Al graficar los valores de la función de semivarianza experimental dados en la tabla anterior (Fig. 8) se observa que en sentido inferior-superior el semivariograma es mayor que en sentido izquierda-derecha, luego la conclusión más relevante para este conjunto de datos es que la estructura de correlación espacial no sólo depende de la distancia entre las observaciones, sino de su orientación. En otras palabras el fenómeno es anisotrópico (debido a que el sill en las dos direcciones es distinto existe anisotropía zonal; Samper & Carrera, 1990)
19
Semivariogramas Experimentales
Semivarianza
30 25 20 15 10 5 0 100
Norte-Sur 200
300
400
Este-Oeste
Distancia
Figura 8. Función de semivarianza experimental en dos direcciones para el conjunto de datos hipotéticos del ejemplo de esta sección. 3.2. Modelos Teóricos de Semivarianza. Existen diversos modelos teóricos de semivarianza que pueden ajustarse al semivariograma experimental (función de semivarianza calculada con los datos muestrales). En Samper & Carrera (1990) se presenta una discusión respecto a las características y condiciones que éstos deben cumplir. En general dichos modelos pueden dividirse en no acotados (lineal, logarítmico, potencial) y acotados (esférico, exponencial, gaussiano) (Warrick et al., 1986). Los del segundo grupo garantizan que la covarianza de los incrementos es finita, por lo cual son ampliamente usados cuando hay evidencia de que presentan buen ajuste. Los parámetros básicos de estos modelos son el efecto pepita, la meseta y el rango (David, 1977).Antes de estudiar los modelos usados para ajustar los semivariogramas experimentales se definirán dichos parámetros: Efecto Pepita Se denota por C0 y representa una discontinuidad puntual del semivariograma en el origen. Puede ser debido a errores de medición en la variable o a la escala de la misma. En algunas ocasiones puede ser indicativo de que parte de la estructura espacial se concentra a distancias inferiores a las observadas. Meseta Es la cota superior del semivariograma. También puede definirse como el limite del semivariograma cuando la distancia h tiende a infinito. La meseta puede ser o no finita. Los semivariogramas que tienen meseta finita cumplen con la hipótesis de estacionariedad fuerte; mientras que cuando ocurre lo contrario, el semivariograma define un fenómeno natural que cumple sólo con la hipótesis intrínseca. La meseta se denota por C1 o por (C0 + C1) cuando la pepita es diferente de cero. Si se interpreta la pepita como un error en las mediciones, esto explica porque se sugiere que en un modelo que explique bien la realidad, la pepita no debe representar mas del 50% de la meseta. Si el ruido espacial en las
20
mediciones explica en mayor proporción la variabilidad que la correlación del fenómeno, las predicciones que se obtengan pueden ser muy imprecisas. Rango Es la distancia a partir de la cual dos observaciones son independientes. El rango se interpreta como la zona de influencia. Existen algunos modelos de semivariograma en los que no existe una distancia finita para la cual dos observaciones sean independientes; por ello se llama rango efectivo a la distancia para la cual el semivariograma alcanza el 95% de la meseta. Entre más pequeño sea el rango, más cerca se esta del modelo de independencia espacial. El rango no siempre aparece de manera explícita en la fórmula del semivariograma. En el caso del modelo esférico (3.2.1), el rango coincide con el parámetro a, que se utilizará en las ecuaciones más adelante. Sin embargo, en el modelo exponencial (3.2.2), el rango efectivo es a/3 y en el modelo gaussiano (3.2.3) es a/√3. 3.2.1. Modelo Esférico Tiene un crecimiento rápido cerca al origen (Fig. 9), pero los incrementos marginales van decreciendo para distancias grandes, hasta que para distancias superiores al rango los incrementos son nulos. Su expresión matemática es la siguiente:
3 h 1 h 3 C ( h) 1 2 a 2 a C1
ha ha
En donde C1 representa la meseta, a el rango y h la distancia. 3.2.2. Modelo Exponencial Este modelo se aplica cuando la dependencia espacial tiene un crecimiento exponencial respecto a la distancia entre las observaciones. El valor del rango es igual a la distancia para la cual el semivariograma toma un valor igual al 95% de la meseta (Fig. 9). Este modelo es ampliamente usado. Su expresión matemática es la siguiente:
3h ( h) C1 1 exp a 3.2.3. Modelo Gaussiano Al igual que en el modelo exponencial, la dependencia espacial se desvanece solo en una distancia que tiende a infinito. El principal distintivo de este modelo es su forma parabólica cerca al origen (Fig.9). Su expresión matemática es:
h2 ( h) C1 1 exp 2 a 21
30
Semivariograma
25 20 Esférico Exponencial
15
Gaussiano 10 5 0 0
50
100
150
200
250
300
Distancia(h)
Figura 9. Comparación de los modelos exponencial, esférico y Gaussiano. La línea punteada vertical representa el rango en el caso del modelo esférico y el rango efectivo en el de los modelos exponencial y gaussiano. Este tiene un valor de 210, respecto a una escala simulada entre 0 y 300. El valor de la meseta es 30 y el de la pepita 0. El 95% de la meseta es igual a 28.5. 3.2.4. Modelo Monómicos. Corresponden a los modelos que no alcanzan la meseta (Fig. 10). Su uso puede ser delicado debido a que en algunos casos indican la presencia de no estacionariedad en alguna dirección. Su fórmula matemática es la siguiente:
( h) kh 02 Obviamente cuando el parámetro es igual a uno el modelo es lineal y k representa la pendiente de la ecuación de regresión con intercepto cero. Gráficamente se pueden representar así: (h)
1< <2 =1 0< < 1
h Figura 10. Comportamiento típico de los modelos de semivarianza monómicos. 22
3.2.5. Modelo Pepita Puro. Es indicativo de carencia de correlación espacial entre las observaciones de una variable (Fig. 11). Es común sumar este modelo a otro modelo teórico de semivarianza, para obtener lo que se conoce como semivariograma anidado. Lo anterior se sustenta en una propiedad de los semivariogramas que dice que cualquier combinación lineal de semivariogramas con coeficientes positivos es un semivariograma. Su expresión matemática es:
0 ( h) C 0
h0 , donde C0 >0 h0
Su representación gráfica es la siguiente: (h)
C0
h
Figura 11. Modelo de semivarianza teórico para variables sin correlación espacial.
23
Capitulo Cuatro Predicción Espacial 4.1. Predicción Espacial Optima. De la teoría de la decisión se conoce que si Z0 es una cantidad aleatoria y Z*0 es su predictor 1, entonces L( Z 0 ; Z 0* ) representa la pérdida en que se incurre cuando se predice
Z 0 con Z 0*
y el mejor predictor será el que minimice E L Z 0 ; Z 0* / Z con Z Z 1 , Z 2 , , Z n , es decir el predictor óptimo es el que minimice la esperanza
condicional de la función de pérdida. Si LZ 0 ; Z 0* Z 0 Z 0* Z *0 E Z 0 / Z . La expresión anterior indica que para encontrar el predictor óptimo se requiere conocer la distribución conjunta de la n+1 variables aleatorias. 2
4.2. Definición de Kriging. La palabra kriging procede del nombre del geólogo sudafricano D. G. Krige, cuyos trabajos en la predicción de reservas de oro, realizados en la década del cincuenta, suelen considerarse como pioneros en los métodos de interpolación espacial. El kriging es un conjunto de métodos de predicción espacial que se fundamentan en la minimización del error cuadrático medio de predicción. En la tabla 6 se mencionan los tipos de kriging y algunas de sus propiedades. En la secciones 4.3 y 4.4, se hace una presentación detallada de ellos. Tabla 6. Tipos de predictores kriging y sus propiedades. TIPO DE NOMBRE PREDICTOR LINEAL
NO LINEAL
PROPIEDADES
Simple Ordinario Universal
Indicador Probabilístico Log Normal, Trans-Gaussiano Disyuntivo
Son óptimos si hay normalidad multivariada. Independiente de la distribución son los mejores predictores linealmente insesgados. Son los mejores predictores independientemente de la distribución.
1
La palabra estimación es utilizada exclusivamente para inferir sobre parámetros fijos pero desconocidos; predicción es reservada para inferencia sobre cantidades aleatorias. 24
Los métodos kriging se aplican con frecuencia con el propósito de predicción, sin embargo estas metodologías tienen diversas aplicaciones, dentro de las cuales se destacan la simulación y el diseño de redes óptimas de muestreo (capitulo 5). 4.3. Kriging Ordinario Suponga que se hacen mediciones de la variable de interés Z en los puntos xi de la región de estudio, es decir se tienen las observaciones Z(x1), . . . , Z(xn), y se desea predecir Z(xo), en el punto xo donde no hubo medición. Lo anterior puede representado en el siguiente esquema: Y
Z(x1) Z(x4) Z(xj)
Z(x2) * Z(x0)
Z(xi)
Z(x3) Z(x5) Z(xn) X
Los puntos negros representan las coordenadas de la región donde se hizo medición de la variable de interés. El asterisco indica la ubicación del punto donde se requiere predecir la variable. Asociado a cada punto hay una correspondiente coordenada X, Y. En esta circunstancia, el método kriging ordinario propone que el valor de la variable puede predecirse como una combinación lineal de los valores medidos así: Z*(x0) = 1 Z(x1) + 2 Z(x2) + 3 Z(x3) + 4 Z(x4) + 5 Z(x5) + . . . + n Z(xn) n
=
i Z ( xi ) i 1
En donde los i representan los pesos o ponderaciones de los valores originales. Dichos pesos se calculan en función de la distancia entre los puntos muestreados y el punto donde se va a hacer la correspondiente predicción. La suma de los pesos debe ser igual a uno para que los errores de predicción tengan promedio cero. Esto último se conoce como el requisito de insesgamiento. Matemáticamente la propiedad de insesgamiento se expresa a través de: E [ Z * ( x0 )Z ( x )] 0
Asumiendo que el proceso es estacionario de media k y utilizando las propiedades del valor esperado, se demuestra que la suma de las ponderaciones es igual a uno: E [ Z * ( x0 )Z ( x0 )] 0 E [ Z * ( x0 )Z ( x0 )] E [ Z * ( x0 )] E [ Z ( x0 )]
25
n
E [ i Z ( xi )] k 0 i 1
1 E( Z ( x1 ) 2 E( Z ( x 2 ) ... n E( Z ( x n ) k k k k
k( 1 2 ... n )k n
n
i 1
i 1
k( i )k i 1
Se dice que Z*(x0) es el mejor predictor porque los pesos se obtienen de tal manera que minimicen la varianza del error de predicción, es decir que minimicen la expresión: VAR [ Z* ( x0 )Z ( x0 )]
Esta última es la característica distintiva del kriging , ya que existen otros métodos de interpolación como el de distancias inversas o el poligonal que no garantizan varianza mínima de predicción (Samper & Carrera, 1990). La estimación de los pesos se obtiene minimizando V [Z* (x0 )Z(x0 )] sujeto a
n
i
1 .
i 1
Se tiene que V [ Z* ( x0 )Z ( x0 )] V [ Z* ( x0 )] 2COV [ Z* ( x0 ),Z ( x0 )] V [ Z ( x0 )] Desagregando las componentes de la ecuación anterior se obtiene los siguiente: n n n V [ Z * ( x0 )] V i Z ( xi ) i j COV [ Z ( xi ),Z ( x j )] i 1 i 1 j 1
Nota: En adelante COV [ Z ( xi ),Z ( x j )] se notara por Cij. V [ Z ( x0 )] 2 n COV [ Z * ( x0 ),Z ( x0 )] COV i Z ( xi ),Z ( x0 ) i 1 n
n
i 1
i 1
iCOV [ Z ( xi ), Z ( x0 )] iCi0
Entonces reemplazando, tenemos que: n
n
n
V [ Z * ( x0 ) Z ( x0 )] i j Cij 2 iCi0 2 i 1 j 1
i 1
n
Luego se debe minimizar la función anterior sujeta a la restricción
i 1 . Este problema i 1
de minimización con restricciones se resuelve mediante el método de multiplicadores de Lagrange: 26
n
n
n
k2 i jCij 2 iCi0 2 i 1 j 1
i 1
n i 1 i 1 Multiplicador deLagrange 0 2
Siguiendo el procedimiento acostumbrado para obtener valores extremos de una función, se deriva e iguala a cero, en este caso con respecto a y :
n
n
n
n
n
( 12C11 21 j C1 j i j Cij ) 2 i Ci0 2 2 i 1
( k2 ) 1
j 2
i 2 j 1
i 1
i 1
1
n 21C11 2 j C1 j 2C10 2 2 j
n
n
j 1
j 1
2 j C1 j 2C10 2 0 j C1 j C10 ( 1 )
De manera análoga se determinan las derivadas con respecto a 2 , ..., n : n n ( k2 ) 2 j C 2 j 2C 20 2 0 j C 2 j C 20 ( 2 ) 2 j 1 j 1
. . . n n ( k2 ) 2 j C nj 2Cn0 2 0 j Cnj Cn0 ( 3 ) n j 1 j 1
por último derivamos con respecto a : n ( k2 ) n 2 i 20 1 1( 4 ) i 1 i 1
De (1), (2), (3), (4) resulta un sistema de (n + 1) ecuaciones con (n + 1) incógnitas, que matricialmente puede ser escrito como: C11 . . . C n1 1
. .
. .
. .
. .
. C 1n . . . . . C nn . 1
1 1 C10 . . . . . . . . . 1 n C n0 0 1
Cij Ci0 por lo cual los pesos que minimizan el error de predicción se determinan mediante la función de covariograma a través de: 27
= Cij -1 Ci0. Encontrando los pesos se calcula la predicción en el punto xo. De forma análoga se procede para cada punto donde se quiera hacer predicción. Los pesos también pueden ser estimados utilizando la función de semivarianza, expresando la función de covariograma en términos de la función de semivarianza, mediante la siguiente relación: Notación: ij = (h), donde h es la distancia entre los puntos i y j, análogamente Cij = C(h) , además 2 = V(Z(x)). 1 2
ij E ( Z ( x j )Z ( xi ))2
1 E ( Z ( x j )) 2 2( Z ( x j )Z ( xi )( Z ( xi )) 2 2
1 1 E ( Z ( x j )2 E Z ( x j )Z ( xi ) E ( Z ( xi ))2 2 2
1 1 E [( Z ( x j ))2 ] k 2 E [( Z ( xi ))2 ] k 2 E Z ( x j )Z ( xi ) k 2 2 2
1 V ( Z ( x )) 1 V ( Z ( x ))COV Z ( x j )Z ( xi ) 2 2
V Z ( x )COV Z ( x j )Z ( xi ) 2 Cij Cij 2 ij
(5)
Reemplazando (5) en (1) ,(2) y (3) se determinan los pesos óptimos en términos de la función de semivarianza: n ( k2 ) n j C1 j C10 j ( 2 1 j ) ( 2 10 ) 1 j 1 j 1
1
n
2 j j 1 j 2 10 j 1
j 1
n
n
j 1
j 1
2 j 1 j 2 10 j 1 j 10
Similarmente, ( k2 ) n j 2 j 20 2 j 1
( k2 ) n j nj n0 n j 1
28
El sistema de ecuaciones se completa con (4). De acuerdo con lo anterior los pesos se obtienen en términos del semivariograma a través del sistema de ecuaciones: .
11 . . . n1 1
. .
. .
. .
. .
.
1n . . .
. . nn . 1
1 . . . n
1 . . . 1 0
10 . . . n0 1
Los pesos de kriging también pueden ser estimados mediante el uso del correlograma aplicando la siguiente relación: ij
Cij
2
.
La varianza de predicción en cada punto es calculada por (Cressie, 1993): n
2o i io i 1
en donde io , y (i) son interpretados de igual forma a como fueron descritos anteriormente. 4.3.1. Validación del kriging. Existen diferentes métodos para evaluar la bondad de ajuste del modelo de semivariograma elegido con respecto a los datos muestrales. El más empleado es el de validación cruzada, que consiste en excluir la observación de uno de los n puntos muestrales y con los n-1 valores restantes y el modelo de semivariograma escogido, predecir vía kriging el valor de la variable en estudio en la ubicación del punto que se excluyó. Se piensa que si el modelo de semivarianza elegido describe bien la estructura de autocorrelación espacial, entonces la diferencia entre el valor observado y el valor predicho debe ser pequeña. Este procedimiento se realiza en forma secuencial con cada uno de los puntos muestrales y así se obtiene un conjunto de n “errores de predicción” . Lo usual es calcular medidas que involucren a estos errores de predicción para diferentes modelos de semivarianza y seleccionar aquel que optimice algún criterio como por ejemplo el del mínimo error cuadrático medio (MECM). Este procedimiento es similar al método Jacknife, una técnica de re-muestreo, empleado en diversos contextos estadísticos para estimar varianzas de estimadores, entre otros aspectos. Una forma fácil de hacer la validación cruzada es mediante un gráfico de dispersión de los valores observados contra los valores predichos. En la medida en que la nube de puntos se ajuste más a una línea recta que pase por el origen, mejor será el modelo de semivariograma utilizado para realizar el kriging.
29
4.3.2. Representación de las predicciones Una vez se ha hecho la predicción en un conjunto de puntos diferentes de los muestrales vía kriging, se debe elaborar un mapa que dé una representación de global del comportamiento de la variable de interés en la zona estudiada. Los más empleados son los mapas de contornos, los mapas de residuos y los gráficos tridimensionales. En el caso de los mapas de contornos, en primer lugar se divide el área de estudio en un enmallado y se hace la predicción en cada uno de los nodos de éste mismo. Posteriormente interpolando se unen los valores predichos con igual valor, generando así las líneas de contorno (isolíneas de distribución). 4.4. Otros Métodos Kriging A continuación se mencionan algunos aspectos generales de otros métodos de predicción espacial. Un estudio riguroso de ellos puede hacerse en Cressie (1993), Deutsch & Journel (1998) y Samper & Carrera (1990) Kriging Simple Su expresión es similar a la del kriging ordinario. Es menos usado en la práctica porque requiere conocer la media del proceso (). El procedimiento para estimar los pesos se resuelve de manera análoga al del sistema kriging Ordinario. n
Z i ( Z i ) * 0
i 1
Kriging Universal Se aplica cuando el proceso estocástico de estudio no es estacionario en la media. Suponga que la tendencia en la media puede ser modelada a través de la siguiente ecuación de regresión:
Z ( si ) m( si ) R( si ) donde Z ( si ) es el proceso estocástico de interés, m( si ) representa la tendencia, que es modelada como una función determinística de las coordenadas geográficas y R( si ) hace referencia al error de estimación . El predictor en este caso tiene la forma
Z 0* m( s0 ) R0* ,
n
R0* i Ri i 1
* * donde Z 0 es la predicción en el sitio de interés y m( s0 ) y R0 corresponden a la tendencia ajustada y la predicción del residual, llevada a cabo a través de kriging ordinario, en este mismo.
30
Kriging Indicador Suponga que se tiene un proceso estocástico espacial Zi.. Con base en los valores observados se construye la siguiente variable indicadora:
1 Si Z i z Ii 0 Otro caso entonces: n
E I 0 / I PI 0 1 / I I 0* i I i , donde I I 1 , I 2 , , I n . i 1
El kriging indicador consiste en hacer una transformación de los valores observados a una variable indicadora (utilizando por ejemplo la mediana o los cuartiles) y posteriormente aplicar kriging ordinario o simple para predecir en sitios de la región de estudio no muestreados probabilidades de que la función indicadora tome el valor 1. Este procedimiento también es válido para procesos estocásticos en los que la variable estudiada en cada sitio es de tipo doble estado (por ejemplo cuando se mide presencia -ausencia de una especie). Kriging Probabilístico Es un predictor basado en cokriging (sección 4.5) que utiliza como variables predictoras una variable indicadora y una variable generada a través de la transformación uniforme. Sea Zi la variable observada, transformaciones:
i = 1,2, . . ., n entonces se definen las siguientes
1 Si Z i z Ii 0 Otro caso R( Z i ) para todo i, i = 1,2,. . . , n Ui n
con R( Z i ) igual al rango (posición que ocupa dentro de los datos ordenados de menor a mayor) la i-ésima observación muestral. La predicción de probabilidad de éxito en el sitios de interés está dada por: n
n
i 1
i 1
I i I i v iU i * 0
Los pesos i y (sección 4.5).
vi se estiman mediante el sistema de ecuaciones del método cokriging 31
Kriging Log-Normal y Multi-Gaussiano Estos dos procedimientos asumen que las variable regionalizada considerada sigue distribución normal en cada punto del dominio. El primero de estos consiste en aplicar kriging simple u ordinario a la transformación logarítmica de los datos. En el segundo se asume que el proceso estocástico sigue distribución normal con igual media y varianza y a cada valor observado le asigna su "score" normal (probabilidad acumulada, hasta el correspondiente valor, bajo el supuesto de normalidad). Posteriormente se realiza kriging simple u ordinario para hacer predicción en sitios no muestreados de las correspondientes probabilidades acumuladas. Estos dos métodos, aunque fáciles de implementar, no son muy realistas porque están sumiendo conocida la distribución de probabilidad y los parámetros de la misma. Kriging Disyuntivo Kriging de transformaciones polinomiales, f i Z i , específicas de los datos. n
Z 0* f i ( Z i ) . i 1
En la derivación del sistema de ecuaciones correspondiente, se emplean conceptos referentes a espacios de Hilbert y polinomios de Legendre, Jacobi y Hermite (Samper & Carrera, 1990). 4.5.Geoestadística Bivariada y Cokriging Ordinario Si se tienen dos variables regionalizadas Zv1 y Zv2 tomadas en cada uno de los puntos. El semivariograma cruzado de estas dos, se estima por: n (6) 1 h v v ( h) Zv1 (x h) Zv1 (x) Zv2 (x h) Zv2 (x) 1 2 2n h
Donde nh es el número de parejas de datos que se encuentran a una distancia h (Bogaert et al., 1995).
4.5.1. Modelo Lineal de Corregionalización (MLC) El MLC asume que todos los semivariogramas simples y el semivariograma cruzado pueden expresarse como una suma de modelos básicos (exponencial, esférico, gaussiano, etc.) idénticos. Para el caso de dos variables:
32
v ( h )
h )... m m ( h )
v ( h )
h )... m m ( h )
0 0(
1
0 0(
2
(7)
v v ( h )
0 0(
1 2
h )... m m ( h )
donde: v1 y v2 son los semivariogramas simples, v1v2 es el semivariograma cruzado. 0(h), 1(h), . . ., m(h) son los modelos básicos de semivariograma y , y son constantes. Matricialmente:
v ( h ) v1v2 ( h ) ( h ) 1 ( h ) v v ( h ) v 2 12
m
B ( h ) , donde s s
s 0
Bs s s
s s
0 s ( h ) 0 ( h ) s
s ( h )
(8)
A (h) se le conoce como matriz de corregionalización. Esta puede ser también calculada con base en covarianzas cruzadas y correlaciones cruzadas, para lo cual se aplican las formulas dadas en la sección 3.1.2. 4.5.2. Cokriging El método de predicción espacial cokriging consiste en hacer predicción espacial de una variable con base en su información y en la de algunas variables auxiliares que este correlacionadas espacialmente con ella. El predictor cokriging tiene la siguiente expresión en el caso en el que se considera una sola variable auxiliar: n1
Zˆ v*1 ( xo )
i 1
n2
b Z
ai Z v1 ( xi )
j
j 1
v2
( xj )
(9)
El lado izquierdo de la igualdad representa la predicción de la variable de interés en el punto x0 no muestreado. Z v1 ( xi ) con i=1, 2 , ... , n1, representa los valores observados de la variable primaria. Así mismo, Z v2 ( x j ) con j=1, 2, . . ., n2, representa los valores observados de la variable auxiliar. ai y bj, con i=1, 2 , ... , n1 y j=1, 2, . . ., n2 respectivamente, representan los pesos o ponderaciones de las observaciones de las variables primaria y auxiliar y se estiman con base en el MLC ajustado a los variogramas simples y cruzados. Los pesos ai y bj se estiman de manera análoga al proceso descrito para el método kriging, es decir estos serán los que minimizan la varianza de predicción sujeta a la restricción de que el predictor sea 33
insesgado. La estimación de los parámetros se obtiene resolviendo el siguiente sistema de ecuaciones (Isaaks & Srivastava, 1989):
v1 ( 1,1 ) ( 1, n ) v1 v1v 2 ( 1,1 ) v1v 2 ( 1, m ) 1 0
v1 ( n ,1 ) v1 ( n , n ) v1v 2 ( n ,1 ) v 1v 2 ( n , m ) 1 0
v1v 2 ( 1,1 ) v1v 2 ( 1, n ) v 2 ( 1,1 ) v 2 ( 1, m ) 0 1
v1v 2 ( m ,1 ) v 1v 2 ( m , n ) v 2 ( m ,1 ) v 2 ( m,m ) 0 1
1 1 0 0 0 0
0 a1 0 an 1 b1 1 bm 0 1 0 2
v1 ( 0 ,1 ) ( 0,n ) v1 v1v 2 ( 0 ,1 ) v 1v 2 ( 0 , m ) 1 0
(10)
La matriz del lado izquierdo contiene los valores de las funciones de semivarianza y de semivarianza cruzada calculadas para todas las distancias entre las parejas de puntos consideradas. Las dos ultimas filas de dicha matriz son las correspondientes a la restricción de insesgamiento del predictor. ai y bj con i = 1, 2, ..., n y j = 1, 2, ...., m, son los parámetros a estimar, 1 y 2 son los multiplicadores de Lagrange empleados para la restricción de insesgamiento y el vector del lado derecho contiene los valores de la funciones de semivarianza y semivarianza cruzada evaluados para las distancia entre los sitios de muestreo (de ambas variables) y el sitio donde se desea hacer la predicción. Las dos últimas filas del vector están asociadas a la condición de insesgamiento. La correspondiente varianza de predicción del método cokriging se calcula como (Bogaert et al, 1995): k2 CovZ v1 x0 , Z v1 x0 1 ai CovZ v1 xi , Z v1 x0 b j CovZ v 2 x j , Z v 2 x0 n
m
i 1
j 1
(11)
donde 1 es el multiplicador de Lagrange empleado para la restricción dado por la n condición de insesgamiento ai 1 . CovZ vi xk , Z vi xl vi2 vivi k ,l es la i 1 función de covarianza espacial de la variable i, i=1,2, evaluada para la distancia entre los sitios de muestreo k, l. La ventaja del método cokriging frente al kriging radica en el hecho de que cuando la variable auxiliar está ampliamente correlacionada con la variable de interés se puede obtener un disminución en la varianza de predicción, no obstante dicha variable tenga menor densidad de muestreo. En situaciones en las que la variable objetivo tiene costos altos de muestreo se recomienda la aplicación de esta metodología (Bogaert et al., 1995). El método cokriging se ha presentado en esta sección en función de dos variables, sin embargo este puede ser extendido de manera natural al caso en que se tengan más de dos variables, sin embargo en estos casos puede ser dispendioso el ajuste del MLC.
34
Apéndice 6.1. Álgebra de Matrices.
La gran mayoría de métodos estadísticos, incluyendo la geoestadística, pueden ser tratados de forma mucho más sencilla a través del uso del álgebra de matrices. Por ésta razón es útil, si no esencial, tener un cierto conocimiento mínimo de ésta área de las matemáticas. Lo anterior es cierto siempre y cuando el interés sea usar los métodos como una herramienta. La notación del álgebra matricial algunas veces puede resultar desanimante. Sin embargo, no es difícil entender sus principios básicos. 6.1.1. Matriz Una matriz A de tamaño (mxn) es un arreglo rectangular de m filas con n columnas.
a 11 a 12 a 21 a 22 . . A . . . . a m1 a m2
. .
. .
. .
.
.
.
a 1n a 2n . . . a mn
6.1.2. Suma y Producto de Matrices El procesos aritmético de adición, sustracción, multiplicación y división tiene sus contraparte con matrices. Si A y D son dos matrices de orden 3x2, entonces su suma se define como:
a11 a12 d 11 A D a 21 a 22 d 21 a 31 a 32 d 32
d 12 d 22 d 23
a11 d 11 a12 d 12 a 21 d 21 a 22 d 22 a 31 d 31 a 32 d 32
En el caso de la multiplicación se debe cumplir que el número de columnas de la primera matriz sea igual ala número de filas de la segunda.
a11 a12 b11 A B a 21 a 22 b 21 a 31 a 32
b12 b 22
a1i b i1 b13 a 2i b i1 b 23 a 3i b i1
a1i b i2 a1i b i3 a 2i b i 2 a 2i b i 3 a 3i b i2 a 3i b i32
6.1.3. Inversa y Determinante de una Matriz.
35
Si k es un número, es cierto que k x k-1 = 1. De forma similar si A es una matriz cuadrada (número de filas igual al número de columnas) su inversa es A-1, donde AA-1 = A1A = I, con I igual a la matriz idéntica (matriz de unos en la diagonal y cero por fuera de ella). Un ejemplo de matriz inversa es: 1 2 1 2 / 3 1 / 3 1 2 1 / 3 2 / 3
Esto puede comprobarse observando que: 2 1 2 / 3 1 / 3 1 0 1 2 1 / 3 2 / 3 0 1
la inversa de una matriz 2x2, si existe, puede determinarse fácilmente por medio del siguiente cálculo:
a11 a12 a 21 a 22
1
a 22 a 21
a12 a11
Donde = a11a22 - a12a21. La cantidad es llamada el determinante de la matriz. Claramente la inversa no está definida si el determinante es igual a cero. Con matrices grandes el cálculo de la inversa es tedioso y se debe usar un programa de computo para realizarlo. 6.1.4. Valores y Vectores Propios. Dada una matriz A de orden (n x n), si existe un vector x (n x 1) y un número tal que Ax = x. ó (A - I)x = 0 donde I es la matriz idéntica de orden (n x n) y 0 es un vector (n x 1), entonces se llama a y x, respectivamente, valor y vector propio de la matriz A. Pueden encontrarse hasta n valores propios y hay tantos vectores propios como valores propios se encuentren. Los valores de deben satisfacer que el determinante de A - I = 0. Los vectores propios se calculan después de reemplazar los valores propios encontrados en la expresión Ax = x. Al igual que con la inversa, para matrices grandes se debe emplear un software especializado para su obtención. A continuación, a manera de ilustración, se realiza el cálculo de los vectores y valores propios de una matriz de orden 2 x 2.
6 3 Sea A , entonces 3 4
6 3 1 0 A I 0 0 3 4 0 1
36
6 3 0 0 3 4 0 (6 ) 3 3 (4 )
0
(6 )(4 ) 9 0 2 10 15 0
b b 2 4ac 2a
( 10) 100 4(15)
2 81623 . , 18377 .
10 40 2
Para cada valor propio existe un vector propio, el cual se obtiene reemplazando el valor propio correspondiente en la primera expresión de la página anterior y usando la condición de que los respectivos vectores propios estén normalizados.
x1 Un vector x se dice que está normalizado si satisface que x2
x12 x 22 1 .
Teniendo en cuenta lo anterior se calculan los vectores propios de la siguiente forma: (A - I)x = 0 3 x1 0 (6 ) (4 ) x 2 0 3
(6 ) x1 3x 2 0 3x1 (4 ) x 2 0 Restando las dos ecuaciones anteriores y factorizando, obtenemos:
x1 (6 3) x 2 (3 4 ) 0 x1 (3 ) x 2 ( 1 ) 0 x1
(1 ) x 2 (3 )
Entonces para = 8.1623 y = 1.8377 se tiene respectivamente: 37
x1 = 1.3847x2 y x1 = -0.7207x2 . Ahora utilizando la restricción de que los vectores estén normalizados se obtiene:
x12 13847 . 2 1 x12
x12 (13847 . ) 2 x12 (13847 . )2
2 (13847 x12 1 13847 . . )2
x12
(13847 . )2
2 1 13847 .
x1
13847 .
2 1 13847 .
0.8107
x1 0.8107 0.5855 . 13847 . 13847 . x1 0.8107 Luego el vector propio asociado al valor propio = 8.1623 es x 2 0.5855 Efectuando un procedimiento similar se puede comprobar que el vector propio asociado al x1 0.5847 valor propio = 1.8377 es x 2 0.8113 Reemplazando el valor de x1, obtenemos que x 2
En resumen dada la matriz del ejemplo entonces se puede comprobar que:
. 0 0.8107 0 6 3 81623 81623 . 0.5855 0 3 4 0 y, con el segundo valor y vector propio, que
. 0 0.5847 0 6 3 18377 18377 . 0.8113 0 3 4 0 6.1.5. Teorema del Valor Singular Sea Xnxp una matriz real, existen V y U ortogonales (VTV = I y UT U=1) y una matriz D diagonal tal que:
X VDU T donde : 1 i. D 0
0 , q
1 2 q
38
ii. Las columnas de V son los vectores propios de XTX y las filas de U son los vectores propios de XTX. Si se calculan los vectores propios de XTX y se multiplica a la derecha por U en la ecuación
X VDU T , entonces:
XU VD . Ahora multiplicando por D-1 a la derecha,
se tiene:
V XUD 1 6.2. Conceptos de Probabilidad A continuación se presenta una revisión no exhaustiva y a manera introductoria de conceptos básicos de la teoría de probabilidades. Un estudio profundo y formal de estos se puede hacer en Mood, Graybill & Boes (1974)
6.2.1. Variable Aleatoria Si X es una función que le asigna a cada uno de los resultados de un experimento aleatorio (aquel cuya respuesta no puede ser establecida de antemano) un número real, entonces X se llama una Variable Aleatoria. Estas pueden ser discretas o continuas. 6.2.2. Función de Probabilidad Si X es una variable aleatoria discreta. Se llamará a f(x) = P (X = x) función de probabilidad de la variable aleatoria X, si satisface las siguientes propiedades: i. f x 0 x RX ii.
f x 1 . x
Si existe una función f(x) tal que: i. f x 0 , x ii.
f x dx 1
iii. Pa X b a f x dx para cualquier a y b, entonces f(x) es la función de densidad de b
probabilidad de la variable aleatoria continua X. La función de probabilidad acumulada, notada como F(x), es igual a P X x y se evalúa a través de una sumatoria o de una integral dependiendo de si X es discreta o continua. 39
6.2.2.1. Valor Esperado y Varianza Si X es una variable aleatoria, el valor esperado de una función de la variable aleatoria X, g X está dado por: g x f x E g X x g x f x dx
X discreta X continua
como caso particular, xf x E X x xf x dx
X discreta X continua
La varianza de la variable aleatoria X está definida como:
V ( X ) 2 E X 2
x 2 f x x x 2 f x dx
X discreta X continua
La raíz cuadrada de la varianza se denomina desviación estándar y se denota por . Se cumple que: 1. 2. 3. 4.
EaX aE X , con a constante
EaX b aE X b , con a y b constantes
V aX a 2V X y a constante
V X E X 2 E X 2
6.2.2.2. Función de Probabilidad Binomial y Normal. Modelo Binomial Suponga que hay un experimento que consiste en examinar n individuos y evaluar o medir en cada uno de ellos si tienen o no una característica dada (sólo hay dos posibles resultados).Sea p la probabilidad de ¨éxito¨ y q = 1-p la de ¨fracaso¨ en cada uno de los n ensayos. Se asume que esta probabilidad es constante en cada uno de ellos.
40
Sea X= Número de éxitos en los n ensayos, entonces asumiendo conocido p entonces es posible establecer las probabilidades de ocurrencia de cada evento mediante la siguiente ecuación, denominada modelo de probabilidad binomial: n P( X x) p x (1 p) n x x 0, 1, 2, . . ., n x En este modelo: E ( X ) np 2 V ( X ) np(1 p)
Modelo Normal El modelo de probabilidad normal (Gaussiano) es útil para encontrar las probabilidades asociadas a eventos de variables aleatorias cuyas distribuciones de frecuencias son simétricas alrededor del valor promedio. Algunos ejemplos de este tipo de variables aleatorias son los siguientes: Sea el valor promedio de la variable (E(X)) y 2 su correspondiente varianza (V(X)), entonces las probabilidades de ocurrencia de eventos asociados a los posibles resultados de la variable estudiada pueden ser encontrados usando la siguiente expresión, llamada modelo de probabilidad normal: b
Pa X b a
x 2
1 / 2 1 e 2
dx .
Obviamente resultaría muy dispendioso tener que calcular estas integrales para cada valor de a, b, y . Por esta razón se acude a un procedimiento llamado estandarización, el cuál consiste en hacer la transformación Z X . La variable anterior tendrá (si la distribución
de frecuencias de X se ajusta a un modelo de probabilidad normal con media y varianza 2) una distribución de frecuencias que se ajusta a un modelo de probabilidad normal con media cero y varianza uno, es decir que: z
1
2 a 1 2z b P( a X b ) e Z z1 Z z2 z 1 2
2
La ecuación anterior también puede resultar difícil de evaluar, sin embargo para cualquier valor de a, b, y las correspondientes probabilidades pueden hallarse, sin necesidad de resolver la integral, empleando la tabla de distribución acumulada normal estándar que aparece en los textos de estadística. 6.2.3. Función de Probabilidad Bivariada. Si X y Y son dos variables aleatorias discretas. La probabilidad de X = x y Y = y está determinada por la función de probabilidad bivariada f x, y PX x, Y y donde : i. f x, y 0, x, y RX , RY ii.
f x , y 1 x
y
41
Si existe una función f x , y tal que la probabilidad conjunta: Pa X b, c Y d
f x, y dydx
b d
a c
para cualquier valor de a, b, c y d en donde f x , y 0 , x, y y entonces f x , y es la función de probabilidad bivariada de X y Y.
f x, y dydx 1
,
La función de probabilidad acumulada F x , y es igual a PX x,Y y y se evalúa a través de una doble sumatoria o de una doble integral dependiendo de si las variables aleatorias son discretas o continuas, respectivamente. 6.2.3.1. Función de Probabilidad Marginal Si X y Y son dos variables aleatorias con función de probabilidad conjunta f x , y . Las funciones de probabilidad marginales de Y y Y están dadas por f x f x , y y
si X y Y son variables aleatorias discretas
f y f x , y x
ó por f x
f y
f x , y dy f x , ydx
si X y Y son variables aleatorias continuas
6.2.3.2. Función de Probabilidad Condicional Sean X y Y dos variables aleatorias con función de densidad conjunta f x , y . La función de probabilidad condicional de la variable aleatoria X, denotada por f x / y , para un valor fijo y de Y, está definida por: f x / y
f y 0 .
f x , y , donde f y es la función de probabilidad marginal de Y de manera tal que f y
De manera análoga, la función de probabilidad condicional de Y para un valor fijo x de X se define como:
42
f y / x
f x 0 .
f x , y , donde f x es la función de probabilidad marginal de X de manera tal que f x
6.2.3.3. Independencia Estadística. Sean X y Y dos variables aleatorias con función de densidad conjunta f x , y . X y Y son independientes si y sólo si: f x , y f x f y
donde f x y f y son las funciones de probabilidad marginales.
6.2.3.4. Valor Esperado, Varianza y Covarianza Sean X y Y dos variables aleatorias que se distribuyen conjuntamente. El valor esperado de una función de X y Y, g x , y , se define como: g x , y f x , y E g X ,Y x y g x , y f x , y dydx
si X y Y son discretas si X y Y son continuas
La covarianza entre X y Y, denotada por Cov (X, Y), se define como: E X X Y Y E XY XY Y X X Y E XY E X EY
donde X y Y representan los valores esperados de X y Y respectivamente. Si la covarianza de X y Y se divide por el producto de las desviaciones estándar de X y Y, el resultado es una cantidad sin dimensiones que recibe el nombre de coeficiente de correlación y se denota por X ,Y . X ,Y
Cov X ,Y
XY
6.2.3.5. Propiedades del Valor Esperado y la Varianza. Si X y Y son dos variables aleatorias con densidad conjunta, entonces se cumple que: 1. E X Y E X EY 43
2. V X Y V X V Y 2CovX ,Y
n
3. V ai X i i 1
ai a j CovX i , X j . n
n
i 1 j 1
Observación: CovX i , X j CovX j , X i y
Cov X i , X i V X i
Como caso particular:
V a1 X 1 a2 X 2 a12V X 1 a22V X 2 2Cov X 1 , X 2 3. Si E X EY , entonces
1 1 1 E X Y 2 V X V Y Cov X ,Y . 2 2 2
6.3. Algunos Métodos Estadísticos. 6.3.1. Regresión Simple En el modelo de regresión simple se establece una relación lineal entre la esperanza condicional de una variable aleatoria Y dados unos valores fijos de una variable X. Modelo Poblacional
Yi 0 1 xi i EY / X i Yˆi 0 1 xi Yi : i-ésimo valor de la variable respuesta o dependiente en la población xi :i-ésimo valor de la variable predictora o independiente en la población 0 y 1 son parámetros poblacionales que representan el intercepto y la pendiente, respectivamente i : i-ésimo error aleatorio en la población. Supuestos del Modelo. 1. E i 0
2 2. V i 3. Cov i , j 0
4. i N 0 , 2
Modelo Muestral
yi ˆ 0 ˆ 1 xi ei yi ˆyi ei
44
yi : i-ésimo valor de la variable respuesta en la muestra, xi : i-ésimo valor de la variable predictora, ˆyi : Estimación del promedio de Y dado el i-ésimo valor de X en la muestra,
ˆ 0 y ˆ 1 son las estimaciones de los parámetros con base en la información muestral,
ei : i-ésimo error muestral. Estimación de 0 y 1 Uno de los métodos de estimación de los parámetros es el de mínimos cuadrados, que consiste en encontrar los estimadores que hacen mínima la suma de cuadrados de los errores, es decir aquellos valores que hacen más pequeña
i2 Yi Yˆi n
n
i 1 n
i 1
i2 i 1
0
2
i2 Yi Yˆi n
n
i 1
i 1
2
.
n
Yi 0 1 xi 2 . Derivando e igualando a cero se obtiene: i 1
n
n
i2
i 1
1
2 Yi 0 1 x i = 0 y
i 1
n
2 X i Yi 0 1 x i = 0. i 1
Al simplificar las dos ecuaciones anteriores y distribuir las sumas se tiene: n
Yi i 1
n
n 0 1 xi i 1
n
n
n
i 1
i 1
i 1
xi Yi 0 xi 1 xi2 Las dos ecuaciones anteriores se conocen como ecuaciones normales. Dadas las realizaciones y1 , y 2 , , y n las ecuaciones pueden resolverse para encontrar los estimados de los parámetros: n
n
i 1
i 1
yi nˆ 0 ˆ 1 xi y ˆ 0 ˆ 1 x ˆ y ˆ x 0
1
xi yi y ˆ 1 x xi ˆ 1 xi2 n
n
n
i 1
i 1
i 1
45
n n yi xi n i 1 xi yi n ˆ 1 i 1n i 1
n n 2 ˆ x i 1 x i i 1 i 1 2
n xi y x i i n n i 1 i 1 i 1 ˆ 1 ˆ 1 xi2 xi y i n n i 1 i 1 n
n
n
n
y i xi
ˆ 1
x i y i i 1
n
i 1
n
i 1
xi n i 1 2 x i n i 1 n
2
xi x y i y i 1
n
x i x 2 i 1
Se puede demostrar que los errores estándar estimados de los estimadores de los parámetros corresponden a:
s ˆ 1
s n
x i x 2 i 1
y s ˆ 0
n x i2 i 1 s n 2 n x i x i 1
, con s
n
ei2 i 1
n2
46
Capitulo Seis Aplicaciones
Esta sección será desarrollada durante la realización del curso con aplicaciones a datos geofísicos reales. Se emplea el software de distribución libre R en el análisis de los datos. Por reserva con la información considerada no se incluyen en el texto ni los datos ni los resultados obtenidos.
47
REFERENCIAS Biau, G., E. Zorita, H. von Storch & H. Wackernagel. 1997. Estimation of precipitation by kriging in EOF space. GKSS, 97, E45. Box, G. E. P. y G.M. Jenkins. (1976). Time Series Analysis Forecasting and Control. Holden -Day, San Francisco, 575pp. Bogaert, P., P. Mahau & F. Beckers. 1995. The Spatial Interpolation of Agroclimatic Data. Cokriging Software and Source Code. FAO, Rome. Carr, J. , D. Myers y Ch. Glass.(1985). Cokriging - A Computer Program. Computers & Geosciences. 11(2), 111-127. Caselton , W. F. & J. V. Zidek. (1984). Optimal Monitoring Desings. Statistics & Probability Letters. 2, 223(27). Clark, I. 1979. Practical geostatistics. Elsevier Publishing , New York. Cressie, N. 1989. Geostatistics. The American Statistician. 43(4), 611(23). Cressie, N. 1993. Statistical for spatial data. John Wiley & Sons, New York. Curia, D. 1998. Geoestadística Aplicada a la Geología del Petróleo. Larriestra, Curia y Asociados, Buenos Aires. David, M. 1977. Geostatistical ore reserve estimation. Elsevier, Amsterdam. Davis J.C., Statistics and Data Analysis in Geology, Wiley, 3rd Edition, 2002, 638 p. Deutsh , C. V. & A. G. Journel. 1992. GSLIB: Geostatistical Software Library and User´s Guide. Oxford University Press, New York. Díaz- Francés, E. (1993). Introducción a Conceptos Básicos de Geoestadística. Memorias Seminario Estadística y Medio Ambiente. Centro de Investigación en Matemáticas, CIMAT. Guanajuato, Mexico. Diggle, P., L. Harper y S. Simon. (1995). Geoestatistical Analysis of Residual Contamination from Nuclear Weapons Testing. Programme Abstracts of the third SPRUCE International Conference. Merida, Mexico. Digby, P & R. Kempton. 1992. Multivariate Análisis of Ecological Communities. Chapman and Hall, Londres Englund, E. & A. Sparks. 1988. GeoEAS, User´s Guide. EPA, Las Vegas. Evangelos A. & G. T. Flatman. 1988. On sampling nonstationary spatial autocorrelated data. Computers and Geosciences, 14, 5, 667-686. Fox, J. 1984. Linear statistical models and related methods with applications to social research. John Wiley & Sons, New York. Gamma Desing. 1995. GS+. Geostatistical software for the agronomic and biological science, version 2.3. Plainwell, Michigan. Garmin International, Inc. 1993. Garmin communication & navigation. GPS 100 SRVY II personal surveyor. Owner´s manual. Lenexa, U. S. A. Giraldo, R., W. Troncoso, J. E. Mancera & N. Méndez. 2000. Geoestadística: Una Herramienta para la Modelación en Estuarios. Rev. Acad. Col. Ciencias. 24(90):59-72.
48
Hair, J. F., R. E. Anderson, R. L. Tatham y W. C. Black. (1995). Multivariate Data Analysis With Readings. Fourth edition. Prentice Hall, New Jersey. Hoaglin, D. F., F. Mosteller & J. Tukey. 1983. Understanding Robust and Exploraory Data Análisis. John Willey & Sons, New York. Isaaks, E. & R. M. Srivastava. 1989. Applied geostatistics. Oxford University Press, New York. Journel, A.G. y Ch. J. Huijbregts. 1978. Mining Geostatistics, Academics Press, New York. 590 pp. Krige, D. G. 1951. A statistical Approoach to Some Basic Mine Valuation Problems on the Witwatersrand. Journal of the Chemical, Metallurgical and Mining Society of South África, 52: 119-139. Link, P. 1992. Basic Petroleum Geology. OGCI Publications, Oil & Gas Consultants Inc Matheron, G. 1962. Traite de Geostatistique Apliquee, Tome I. Memoires bureau de Recherches Geologiques et Minieres, N 24. Editions Bureau de Recherche et Minieres, Paris. Mardia, K. y C. Goodall, C. 1990. Spatial - Temporal Analysis of Multivariate Environmental Monitoring Data. LATEX documents for the North - Holland Series in Statistics and Probability. McBratney, A. B., Webster, R. and Burgess, T. M. „The design of optimal sampling schemes for local estimation and mapping of regionalized variables I‟, Computers and Geosciences, 7(4), 331-334, 1981 Mood, A., F. A. Graybill & D. C. Boes. 1974. Introduction to the Theoty of Statistics. McGraw-Hill, New York. Morton, T, Thompson, A. and M. Woods. 1997. Development Geology Reference Manual, AAPG Methods in Exploration No. 10. Myers, D. E. 1987. Optimization of Sampling Locations for Variogram Calculations. Water Resources Research. 23(3), 283(93). Petitgas, P. 1996. Geostatistics and their applications to fisheries survey data 5: 114-142. In: B. A. Megrey & E. Mosknes, (E). Computers and fisheries research. ChapmanHall, Londres. Samper, F.J. & J. Carrera 1990. Geoestadística. Aplicaciones a la Hidrogeología Subterránea. Centro Internacional de Métodos Numéricos en Ingeniería. Universitat Politécnica de Catalunya. Barcelona, España. 484 pp. Stanley; T. Practical Statistics for Petroleum Engineers. Petroleum Publishing Company, Tulsa 1973. Wackernagel. H. 1995. Multivariate geostatistics. An introduction with applications. Springer-Verlag, Berlín. Warrick, A. W., D. E. Myers & D. R. Nielsen. 1986. Geostatistical methods applied to soil science. Methods of soil analysis. Part 1. Physical and mineralogical methodsagronomy monograph 9 : 53 - 81.
49
