CentroGeo
2010 Práctica 3
Interpolación Usando Interpolación kriging Ordinario Prof. Dr. Darío Rojas A. Ing. J.
Javier Martínez C.
Alberto ITESM 28/08/2010
CENTRO DE INVESTIGACION INVESTIGACION EN GEOGRAFIA Y GEOMATICA “ING. JORGE L. TAMAYO ” DIPLOMADO EN GEOMÁTICA MÓDULO DE GEOESTADISTICA Agosto 2009 Práctica 3. KRIGING ORDINARIO A partir de un conjunto de datos (De Wijs), hacer interpolaciones de concentraciones en puntos no medidos utilizando el método de Kriging Ordinario.
2) Introducción. En un principio el objetivo principal de la geoestadística fue el de obtener estimaciones de los valores de una variable determinada en lugares sin mediciones. El término Kriging fue introducido por Matheron en 1953 en honor D.G. Krige, un ingeniero en minas sudafricano que desarrollo métodos para estimar las concentraciones de oro y otros metales en las minas. El método Kriging provee una solución al problema de la estimación basada en un modelo continuo de variación espacial estocástica. Hace el mejor uso del conocimiento existente tomando en cuenta la forma en que la propiedad varía en el espacio a través del modelo de variograma. En su forma original una estimación de Kriging en un lugar es simplemente una combinación lineal de los datos en la vecindad. Desde entonces el método ha sido modificado para tratar con problemas más complejos en minería, ingeniería petrolera, control de la contaminación y salud. El objetivo de Kriging es el de estimar el valor de una variable aleatoria, Z, en uno o más puntos no muestreados o sobre grandes bloques, a partir de datos más o menos distribuidos en una plataforma determinada, digamos, z(x1), z(x2),…,z(xn) en x1, x2,…,xn. Los datos pueden estar distribuidos en una, dos o tres dimensiones aunque sus aplicaciones en problemas medioambientales suelen ser bidimensionales. El Kriging Ordinario es el que se aplica más comúnmente. Asume que la media es desconocida. Sí consideramos una estimación puntual, entonces se estima Z en el punto x0 con: N
Z * ( x 0 )
= ∑ λ i z ( xi ) i =1
que es un promedio ponderado de los datos (siendo λi los pesos). Para asegurar que el estimador no sea sesgado los pesos deben sumar 1 y el error esperado es: E[Z*(x0) – Z(x0)]=0. Para obtener los pesos se debe resolver el sistema de ecuaciones: Aλ=b En donde;
γ ( x1, x1)
γ ( x1, x2)
...
γ ( x1, xn)
1
...
...
...
...
1
...
...
...
...
1
γ ( xn, x1)
γ ( xn, x2)
...
γ ( xn, xn)
...
1
1
1
...
0
A=
λ 1 λ 2 λ =
. ..
λ
n
ψ ( x0 ) γ ( x1, x 0) γ ( x 2, x 0)
b=
... γ ( xn, x 0)
1
Siendo γ(xi,xj) el valor del variograma entre los puntos i y j, λi el peso i, y γ(xi,x0) el variograma entre el punto i y el punto a interpolar (x0)
2.1) Datos De Wijs. Los datos que se utilizaran en esta práctica son los datos Dewijs (Dewijs.xlsx) con muestras cada 4m, ( hoja Dewijs 4m). Los datos originales se encuentran en la Hoja Dewijs 2m y serán utilizados más adelante para validar las interpolaciones.
3) Obtención de λ. 3.1) Obtención del variograma ajustado (modelo teórico). El primer paso para hacer interpolaciones consiste en obtener un variograma ajustado a un modelo teórico. En este caso puede ajustar el variograma esférico incluido en la hoja Dewijs 4m variando los valores para el Sill y el Rango. En azul se pinta el variograma experimental (clases de distancia de 4m) y en rojo el variograma esférico.
4) Proceso de interpolación. 4.2) Obtención de la matriz A. Para construir la matriz A es necesario evaluar la función del variograma esférico ajustado para cada distancia dada entre los puntos muestreados. En este caso contamos con 59 muestras. Así, la primera muestra es x1, la segunda x2 y así sucesivamente hasta x59. Las primeras cuatro muestras son: Posición (m)
Zn 0 4 8 12
17.7 9.5 4.1 12.4
Para la matriz A sólo nos interesan los valores de la localización de las muestras, no los de concentración del Zn. Por ejemplo: γ
(x1,x1) = γ (0,0) = γ (0) Donde γ es el variograma esférico ajustado evaluado en h=0, que es la distancia de la primera muestra consigo misma. Otros ejemplos: γ
(x1,x2) = γ (0,4) = γ (4)
γ
(x2,x4) =
γ
(4,12) =
γ
(8)
Recurde que para el variograma esférico:
h h γ ( h) = c1.5 − 0.5 a a
si h<=a 3
γ ( h)
=
c
si h>a
Finalmente se requiere agregar una columna y reglón finales con 1s y un 0 en la posición de la esquina inferior derecha de la matriz, es decir A60,60 = 0. Para calcular la matriz A de manera rápida podemos hacer uso del editor de Visual Basic de Excel. Para cargarlo se deben seguir los siguientes pasos: 1) Hacer click click en el botón botón de office office de Excel Excel y elegir elegir “opciones “opciones de Excel”. Excel”.
2) En la pestaña “Más frecuente frecuentes” s” habilitar habilitar la opción opción de “Mostrar “Mostrar ficha_programador en la cinta de opciones”
En la pestaña “programador” elegir “Visual Basic”.
Escriba el siguiente código para la macro 1 sustituyendo las palabras SILL y RANGO por los valores numéricos elegidos por usted para el variograma esférico: Sub Macro1() ' ' Macro1 Macro ' For x = 2 To 60 For y = 2 To 60
If Abs(Hoja2.Cells(x, 1) - Hoja2.Cells(y, 1)) <= RANGO Then Cells(x - 1, y - 1) = SILL * (1.5 * (Abs(Hoja2.Cells(x, 1) - Hoja2.Cells(y, 1)) / RANGO) - 0.5 * (Abs(Hoja2.Cells(x, 1) - Hoja2.Cells(y, 1)) / RANGO) ^ 3) Else Cells(x - 1, y - 1) = SILL End If Next y Next x For x = 1 To 60 Cells(x, 60) = 1 Cells(60, x) = 1 Next x Cells(60, 60) = 0 End Sub
Finalmente cierre la ventana de Visual Basic, vaya a una hoja nueva de su archivo de Excel, en la pestaña de programador oprima el botón de Macros y ejecute la Macro 1.
4.3) Obtención de A-1. A λ=b implica que A-1A λ = A-1b Es decir, λ= A-1b En Excel se puede obtener una matriz inversa de A mediante la función MINVERSA. Se recomienda calcular esta matriz en una hoja nueva . En caso de que haya calculado la matriz A en la hoja 3 de Excel la formula es: =MINVERSA(Hoja3!A1:BH60) Comenzando por la celda de la fórmula. Presione F2 y, a continuación, CTRL+SHIFT+ENTER (ver ayuda de Excel).
4.4) Cálculo de la matriz b. Para construir la matriz b utilice alguno de los puntos que aparecen en los datos originales De Wijs (hoja 2m) y que no aparezcan en los datos muestreados cada 4 metros, es decir, 2,6,10,14,18,… Recuerde que lo que se usa para calcular la matriz b es la posición, no la concentración de Zn. Puede construir esta matriz en una hoja nueva de Excel. Copie la columna A (posición) de la hoja Dewijs4m y en la columna A de la hoja nueva. En la primera celda de la columna B de la nueva hoja escriba la h a la cual desea hacer la interpolación, por ejemplo 2. Posteriormente implemente en la columna la matriz matriz b con la siguiente fórmula del variograma esférico(no olvide agregar un 1 al final de la columna): =SI(ABS(C$1-$A2)<=RANGO,SILL*(1.5*ABS(C$1-$A2)/RANGO - 0.5*(ABS(C$1-$A2)/RANGO)^3),SILL)
(Nota: Los valores rango y sill tienen que establecerse)
La fórmula se puede copiar hacia la derecha de tal forma que se puedan construir otras matrices b para interpolar en otros puntos.
4.5) Obtención de la matriz de pesos λ.
El siguiente paso consiste en resolver el sistema de ecuaciones λ= A-1b que nos permite obtener los valores de la matriz λ . Como en los casos anteriores, se recomienda calcular la matriz lambda en una hoja nueva. Para multiplicar matrices utilice la función MMULT de Excel (ver ayuda de Excel). Ejemplo: =MMULT('Matriz Inversa de A'!$A1:$BH60,'Matriz b'!B2:B61)
4.5) Interpolación. Interpolación. Finalmente realice la interpolación en el punto elegido (x0).
Z(x0) = λ 1z(x1) + λ 2z(x2)+...+ λ nz(xn) Donde z(xn) es el valorde la n muestra de Zn.
4.6) Realice el proceso de interpolación para por lo menos veinte puntos más. 4.7) Gráfica de valor real vs. Valor interpolado. Realice una gráfica de dispersión XY seleccionando la posición de las interpolaciones, los datos interpolados y los datos muestreados para los mismos puntos. También calcule el RMSE entre ambas series.
