Estimación de Recursos Mineros
Estimación de Recursos Mineros
Estimación de Recursos Mineros
por
Mayo, 2007
Prefacio
Todo lo malo que encuentre en este texto se debe a mi persona y todo lo bueno se
debe a mis dos maestros, profesores de la Escuela Nacional Superior de Minas de
París: Georges Matheron, creador de la geoestadística, quién me enseñó esta
apasionante disciplina y Phillipe Formery quien me inculcó el “gusto” por las
probabilidades.
1
Tabla de Materias
Prefacio 1
Tabla de materias 2
Capítulo 1: Los métodos tradicionales de estimación de recursos 3
La media aritmética 4
Los polígonos 6
El método del inverso de la distancia 7
Crítica general de los métodos tradicionales de estimación de leyes 10
2
Capítulo I
Se tiene un conjunto de leyes z1, z2, . . . , zN de mineral localizadas en los puntos x1,
x2, . . . , xN
S, con el fin de localizar las zonas ricas y pobres dentro de esta zona S.
3
Figura I.2: Estimación local con bloque unitario o unidad básica de cálculo.
1 + 1 + 1 + 3 + 2 + 2 + 1 11
zˆS = = = 1.57
7 7
4
La fórmula general es:
N
1
zˆS =
N
∑z
i =1
i
la figura anterior existe una agrupación de datos en la zona de alta ley: el valor
Figura I.4: Una planta en el depósito de MM. En rojo la intersección de los sondajes.
5
I.2 Los polígonos.
El método de los polígonos se basa en lo siguiente: asignar a cada punto del espacio
la ley del dato más próximo. Para estimar una zona S se ponderan las leyes de los
Figura I.5: Método de los polígonos. Hay que calcular el área de 7 polígonos.
zˆ S = 1 .3 6
1 N
zˆS = ∑ Si zi (S = S1 + S 2 + ... + S N )
S i =1
Comentarios:
6
• En general no es adecuado en estimaciones locales porque asigna la misma ley
a todos los bloques que están dentro de un mismo polígono. Produce problemas
Figura I.6: Salar de Atacama (se trata de salmueras con contenidos de litio, boro, potasio, sulfatos).
Puntos de muestreo y de bombeo. Observar la agrupación (cluster) de los puntos en la zona central.
El método del inverso de la distancia se basa en lo siguiente: asignar mayor peso a las
7
α
ponderar las leyes por 1/ di , (α = 1, 2, . . . ; di = distancia entre la muestra i y el
Figura I.7: Método del inverso de la distancia a la potencia alfa. Hay que calcular las distancias entre los
datos y el centro del bloque
N
zi
∑ dα
i =1
zˆS = N
i
(α > 0)
1
∑ dα
i =1 i
8
Comentarios:
A veces, para evitar el problema de las agrupaciones de datos, se utiliza una búsqueda
Figura I.8: Búsqueda octogonal. Se podrían utilizar las k (¿cómo definir k?) muestras más cercanas al
centro del bloque, dentro de cada octante.
9
En este ejemplo solo las muestras 1, 2, 3 intervienen en la estimación del bloque S.
métodos estudiados:
• Son empíricos...
• Demasiado geométricos.
leyes son erráticas y otros más favorables en los cuales las leyes son
regulares.
Figura I.9: Arriba, sondaje vertical en mina de hierro La Japonesa. Abajo, sondaje vertical en mina de
hierro Boquerón Chañar. Ambos tienen la misma ley media, 32% Fe. ¿Cuál es más continuo?
10
ii) la posible presencia de anisotropías, es decir direcciones en las cuales la
Figura I.10: Anisotropías de las leyes de cobre en una planta del depósito RT. dimensiones en metros.
Ejemplo: en el caso de la figura I.11, las leyes z1 = 1.25 y z2 = 1.75 son simétricas con
Figura I.11: Bloque e isoleyes anisotrópicas. En Chuquicamata la situación es similar, las leyes decrecen
sistemáticamente de oeste a este, ver figura I.12.
11
Figura I.12: Proyección horizontal de los sondajes de la mina Chuquicamata. Se observa la degradación
de leyes en la zona central, de oeste a este.
z1 + z2 1.25 + 1.75
zˆS = = = 1.5
2 2
12
Sea zs la ley verdadera desconocida de S. Sería interesante poder escribir una
zS = zˆS ± error
un yacimiento todo o nada: o bien hay mineral con ley 1.0 o bien hay estéril con
ley 0.0. El mineral se representa de color gris en la figura siguiente (el valor
sondajes positivos en rojo, con ley 1.0 y negativos en blanco, con ley 0.0):
Figura I.13: Oficina Ossa. Depósito de nitratos-yodo. Malla inicial de reconocimiento de 400mx400m.
13
Debido a que los sondajes caen en el centro del bloque, no se puede aplicar inverso de
por la ley del sondaje central. En este caso se puede calcular el porcentaje
Observamos que el promedio de las leyes verdaderas es 0.55 y que el promedio de las
Sin embargo, existe el sesgo condicional: para leyes altas ocurre que siempre la ley
estimada es superior a la ley real, y para leyes bajas, la ley estimada es siempre infe-
Al aplicar una ley de corte sobre las estimaciones, ocurre entonces que la ley mina es
El sesgo condicional se puede comprobar en una mina a cielo abierto, al comparar las
leyes estimadas de los bloques (modelo de largo plazo) con el promedio de los pozos
14
Figura I.14: Pozos de tronadura en mina Lomas Bayas. Leyes de cobre. Dimensiones en metros. Se
observa un bloque de 25mx25mx15m.
La figura I.15 muestra la información de largo plazo y de corto plazo en la mina Los
Bronces.
Figura I.15: Información de largo plazo (sondajes, en negro) y de corto plazo (pozos de tronadura, en
colores). Los sondajes están agrupados en las zonas de alta ley.
15
Capítulo II
Sea x un punto del espacio. Se designa la variable regionalizada por la notación z(x).
por la letra x. Por ejemplo la ley en el punto x se representa por z(x). Por consiguiente,
16
Se observa que existen problemas de notación: se acostumbra a designar una variable
regionalizada con la letra z, lo cual coincide con la notación utilizada para la cota o
elevación.
galería:
17
Figura II.3: Leyes de canaletas entre A y A’. Mina Carola.
Fig. II.4: Precio del cobre (promedio mensual 1987-2005) en centavos de dólar / libra.
18
Figura II.5 (a): Depósito de nitratos-yodo: la zona mineralizada, de color rojo en la figura, se llama
caliche.
Figura II.5 (b): Isópacas en Oficina Sur Lagunas (depósito de yodo). Las potencias están en metros. Los
límites corresponden a propiedad minera.
19
Ejemplo 4: En el espacio de tres dimensiones, sea z(x1, x2, x3) = z(x) = Ley de Cu en
Figura II.6: Caso típico de depósito de óxidos-sulfuros. La capa superior corresponde a grava o coluvio.
Figura II.7. Planta en mina RT. Leyes de bloques de 25mx25mx15m. Zona de óxidos.
20
En un depósito de este tipo se puede comprobar que la ley de cobre se comporta de
Ejemplo 6: En el espacio de tres dimensiones, sea z(x1, x2, x3) = z(x) = densidad de
21
Figura II.9: Densidades superficiales en ton/m3 en la mina Chuquicamata. La falla (llamada falla oeste) la
cual delimita mineral de estéril queda representada por el cambio de densidades.
Ejemplo 7: en el espacio de tres dimensiones sea z(x1, x2, x3) = z(x) = ley de arsénico
en un depósito de cobre. Esta variable regionalizada corresponde a una impureza
(figura II.10)
22
Figura II.10: Leyes de arsénico en un banco de mina Los Pelambres. Los datos corresponden a pozos
de tronadura en bancos de 15 metros. Se trata de una variable muy irregular, concentrada en vetas.
Observación:
Los ejemplos anteriores nos muestran que una variable regionalizada es simplemente
una función z(x) del punto x. Sin embargo, esta función no se comporta como las
adecuado, por ejemplo, en la figura II.6 se podrían distinguir dos campos disjuntos, los
23
Entonces en un mismo depósito minero D pueden haber varios campos o unidades D1,
Figura II.11: Unidades D1, D2, D3, D4 en una sección del depósito de cobre porfídico de Inca de Oro.
24
Figura II.12: Un testigo. Tiene un cierto largo l.
En el caso en que los testigos que constituyen el sondaje son de tamaño irregular, es
Figura II.13: Regularización de un sondaje a un largo constante b. Esta operación produce errores.
25
La figura II.14a muestra una sección transversal en un depósito de óxidos de cobre.
Las líneas representan los sondajes de exploración. El punto rojo se denomina collar
del sondaje. El collar está caracterizado por las coordenadas x0, y0, z0 y por dos
ángulos: (θ, φ) azimuth e inclinación.
Figura II.14a: Sección en el depósito de cobre de RT. Se observan las unidades grava ( estéril), lixiviado,
óxidos y sulfuros. Un compósito está caracterizado por sus coordenadas x, y, z, las leyes de cobre total,
de cobre soluble, un código que indica la unidad, además del nombre del sondaje que contiene al
compósito.
Cada compósito está caracterizado por sus coordenadas x, y, z, sus leyes, un código
que indica el dominio o unidad geológica y la identificación del sondaje, eventualmente
otra información. Se tiene así la base de datos de sondajes del depósito, la cual, en
formato de texto, puede ser incorporada a cualquier paquete computacional.
Para tratar las desviaciones de los sondajes, se divide el sondaje en tramos rectilíneos
L1, L2, …, Lr.
26
Figura II.14b: Azimuth θ (se mide desde el norte) e inclinación φ (se mide desde la horizontal) de un
sondaje.
estimaciones.
Ejemplo: la figura II.15 siguiente representa el caso de una variable regionalizada z(x) =
ley de cobre definida en un soporte cuadrado de lado axa:
27
La ley de corte es c = 0.5
B=T(m–c)
Figura II.15 (a): La zona explotada varía según el tamaño del bloque unitario.
28
Figura II.15 (b): Importancia económica del soporte y la anisotropía. A medida que aumenta el soporte,
se diluyen las leyes. Observar que la ley de corte es mayor que la ley media. Repetir los cálculos para
una ley de corte de 0.40
29
II.5 El modelo matemático de la geoestadística: las funciones aleatorias.
Una función aleatoria es una función Z(x) que asigna a cada punto x del espacio un
Al hacer un experimento sobre la función aleatoria se obtiene una función ordinaria (no
Figura II.16: Función aleatoria y variable regionalizada. Los colores indican rangos de la variable.
30
Un ejemplo de función aleatoria es la siguiente:
Figura II.17: Leyes de litio en el Salar de Atacama para dos campañas de reconocimiento en el tiempo.
Se observan los puntos de bombeo y el decrecimiento de las leyes en el tiempo. No se indican las
escalas de las leyes.
Observaciones:
Esto tendría el mismo sentido que decir “el número 6 es una variable aleatoria”:
31
El enunciado correcto de la hipótesis probabilística de la geoestadística es: “z(x)
b) Para que esta hipótesis probabilística tenga un sentido real, es necesario poder
cual supone que la inferencia estadística (es decir el cálculo de parámetros que
USACH, 2005.
32
Capítulo III
El variograma
geoestadística.
siguiente:
G G
γ (h ) =
1 ⎡
( G 2
E Z (x + h ) − Z (x) ⎤
2 ⎢⎣ ⎥⎦ )
33
Esta ecuación es la que hay que adaptar en cada situación práctica (mallas regulares e
γ (0) = 0
γ ( h) ≥ 0
G G
γ (−h ) = γ (h )
La última relación proviene del hecho que si dos leyes z1 y z2 están a la distancia h,
34
Figura III.3: Línea recta con muestras regulares.
( z2 − z1 ) 2 + ( z3 − z2 ) 2 + ... + ( z N − z N −1 ) 2
γ (b) =
2( N − 1)
( z3 − z1 ) 2 + ( z4 − z2 ) 2 + ... + ( z N − z N − 2 ) 2
γ (2b) =
2( N − 2)
( z4 − z1 ) 2 + ( z5 − z2 ) 2 + ... + ( z N − z N − 3 ) 2
γ (3b) =
2( N − 3)
35
N −k
1
γ (kb) = ∑
2( N − k ) i =1
( zi + k − zi ) 2
Posteriormente los valores γ (0), γ (b), γ (2b), γ (3b), ... se llevan a un gráfico:
G
Estudiaremos el comportamiento de la función γ ( h) para | h | pequeño, para lo cual
36
Caso 1: Leyes muy regulares y continuas.
Para una distancia b pequeña, las dos leyes de la figura son casi iguales, lo que implica
que para ⏐h⏐ pequeño, γ(h) es próximo a cero; luego el gráfico de γ(h) en una
37
Caso 2: Continuidad y regularidad promedio
Figura III.7: Leyes con continuidad promedio. La variable es continua pero no es derivable.
En este caso, para una distancia pequeña, la diferencia de leyes es significativa; luego
el gráfico de γ(h) en una vecindad del origen será:
38
Caso 3: Existencia de micro variaciones:
Figura III.9: Presencia de una estructura a menor escala. La variable es más discontinua.
un crecimiento más moderado (debido a la variación a gran escala): se dice que existe
39
En la práctica la equidistancia b es mayor que d y se tendrá un gráfico del tipo:
El nombre efecto de pepita proviene del estudio de los depósitos de oro. Consideremos
40
Figura III.12: Efecto de pepita en un testigo de una mina de oro.
41
Por muy pequeña que sea la distancia b, las leyes de dos puntos a esta distancia son
prácticamente independientes. El gráfico de γ(h) será:
γ(0) = 0, γ(h) = C si h ≠ 0.
III.15:
42
III.2 Comportamiento del variograma para grandes distancias.
Estudiaremos ahora el comportamiento de la función γ(h) para |h| grande, para lo cual
Se dice que existe una deriva o tendencia. Al hacer el cálculo se observará que γ(h)
siempre crece:
43
Caso 2: Leyes con pseudo-periodicidades:
44
Se dice que el variograma presenta efecto de hoyo o de agujero. En la figura, d = 9
unidades proporciona una medida del pseudo-período ; ∆ es una medida de la
45
Se observa que a partir de una cierta distancia, del orden a = 6 unidades, la función
γ(h) permanece aproximadamente constante:
Esto quiere decir que da lo mismo que la distancia que separa los puntos sea 6, 7, 8 o
más unidades; en otras palabras, dos puntos cuya distancia sea superior a a = 6
El alcance proporciona una medida de la zona de influencia de una muestra porque dos
46
Figura III.23: Zona de influencia de una muestra localizada en el punto x0.
Dos muestras cuya distancia sea inferior al alcance a están correlacionadas entre sí.
C =σ2
47
Figura III.24: Sondaje vertical en una mina de carbón. Los números superiores representan la potencia
real de los mantos.
⎧1 si x ∈ carbón
z ( x) = ⎨
⎩0 si x ∉ carbón
i) Cálculo de la media:
N
1 40
m=
N
∑z
i =1
i =
80
= 0.5
Este resultado nos indica que la mitad del sondaje es carbón y la otra mitad estéril.
1 N
40(1 − 0.5) 2 + 40(0 − 0.5)2
σ =2
N
∑
i =1
( zi − m) =2
80
= 0.25
48
iii) Cálculo de γ(h):
Se obtiene fácilmente:
γ(5) = 0.12
γ(10) = 0.20
γ(15) = 0.24
γ(20) = 0.25
γ(25) = 0.25
49
El alcance tiene una significación geológica en este caso: corresponde al espesor
50
Figura III.27: componentes del vector h
i) Fijemos la dirección θ del vector h; que sea por ejemplo θ = 90º, es decir, la
dirección NS. El vector h sólo puede ser:
51
Calculemos γ(h1) = γNS(10). Al aplicar el algoritmo hay que considerar las
(diferencias)2 posibles: (zi - zj)2 cuando ambos datos zi y zj están definidos. La figura
Figura III.29: Parejas posibles para calcular gama de 10 metros en la dirección NS (hay 36 vectores).
Luego:
52
ii) Sea ahora la dirección θ = 0º, es decir la dirección EW. El vector h sólo puede
ser:
Figura III.31: Parejas posibles para calcular gama de 10 metros en la dirección EW (hay 36 vectores).
Se obtiene entonces:
53
γ(h1) = 0.0146 (36 parejas)
γ(h2) = 0.0330 (33 parejas)
γ(h3) = 0.0431 (27 parejas)
iii) La práctica demuestra que para estudiar las estructuras basta con calcular γ(h)
Figura III.32: Cálculo de gama de h en la dirección de 45°. La distancia entre parejas es ahora 14.41
metros.
Figura III.33: Cálculo de gama de h en la dirección de 135°. La distancia entre parejas es 14.41 metros.
54
Gráfico de γ(h):
Figura III.34: Variograma anisótropo. La variación de las leyes es más regular en la dirección EW que en
la NS.
Se observa una clara anisotropía que nos indica que el fenómeno es más regular en la
dirección EW que en la NS. (Esto se puede comprobar al mirar como varían las leyes
Ejemplo: Los datos que se dan a continuación provienen de un banco en una mina de
fierro:
55
Figura III. 35: Datos de leyes de fierro.
56
Figura III.37: Variograma EW
57
Observamos que γ(h) es casi el mismo según las direcciones: podemos concluir que el
fenómeno es isótropo.
ponderado de los valores del variograma (ponderación por el número de parejas N'):
⏐h⏐ γ(h)
10.00 4.17
14.41 5.73
20.00 8.31
28.28 11.60
30.00 11.94
58
Ejemplo: la figura III.41 muestra la fotografía aérea de un bosque de coníferas. La
59
Figura III.43: Variogramas bosque de coníferas.
periodicidades).
⎧1 si x ∈ granos (negro)
z ( x) = ⎨
⎩0 si x ∉ granos (blanco)
variogramas!
En este caso de dos estados (negro y blanco), si se considera que z(x) es la realización
de una función aleatoria Z(x), entonces se dice que Z(x) es un conjunto aleatorio.
60
Figura III.44: Aproximadamente isótropo.
Figura III.45: Aproximadamente isótropo, con fronteras difusas. Presencia de efecto de pepita
61
Figura III.47: Presencia de dos estructuras, a pequeña escala y a gran escala.
Figura III.48: Fallas. La figura se obtuvo dibujando rectas al azar, se obtienen así polígonos. Cada
polígono se pinta de negro si al tirar una moneda sale cara, en caso contrario se deja igual..
conjunto aleatorio), se puede tomar una imagen y asignar a cada pixel su código
62
Figura III.49: Variograma de una imagen (agua). Observar que la variable regionalizada es estacionaria
en la dirección este-oeste, mientras que en la norte-sur tiene una deriva, la cual se manifiesta solo para
grandes distancias: Se dice que la realización es “localmente estacionaria” en la dirección NS.
promedio; este promedio es bueno cuando el número N' de parejas es grande. Sin
embargo, a medida que ⏐h⏐ crece, N' decrece; la práctica justifica entonces la regla
siguiente:
esta dirección hasta una distancia máxima del orden de 200 metros:
63
Figura III.50: Compósitos proyectados en una planta de la Extensión Norte de Mina Sur. El ancho del
cuerpo es del orden de 400 metros en la dirección EW.
Ejercicio:
Los datos de la figura corresponden a una zona dentro del yacimiento de nitratos-yodo
de Lagunas. Encontrar la función γ(h) hasta ⏐h⏐ = 250 m. según las direcciones: θ
64
Figura III.51: Leyes de NaNO3 Oficina Lagunas. Profundidad del sondaje = 2 metros.
65
c) Cálculo del variograma para mallas irregulares.
66
Lo más probable es que no encontremos ningún o muy pocos pares de datos que estén
aproximadamente a la distancia h.
67
Este método de aproximación presenta problemas:
diferencias en el cálculo.
Algunos paquetes computacionales definen otro tipo de zona para evitar este problema
68
En este caso hay que definir tres parámetros: θ, ε y d (d se llama a veces ancho
de banda).
Hay que tener presente que es necesario conocer bien el variograma en una vecindad
Figura III.57: Compósitos en el espacio de tres dimensiones, con leyes de cobre en mina Chuquicamata.
Comparar con la figuras I.11 y II.9.
69
Figura III.58a: Aproximación en el espacio de 3 dimensiones: Una especie de cono.
70
i) El variograma experimental, calculado a partir de los datos.
variograma experimental:
teórico.
Para tener un buen ajuste, hay que considerar que uno de los objetivos finales es la
bloque:
71
Figura II.60a: El modelo de bloques es fundamental para la planificación minera.
Figura III.60b: Modelo tridimensional de bloques. Existen tres unidades geológicas (a bloque completo)
Figura III.61: Vecindad de estimación. Para estimar el bloque de la figura solo se utilizan los compósitos
que están dentro del círculo de radio R.
72
Si la vecindad de búsqueda es circular o esférica, sólo se utilizará la función γ(h) hasta
una distancia máxima de ⏐h⏐ = 2R; luego conviene ajustar γ(h) hasta ⏐h⏐ = 2R.
darán los mismos resultados; pero el modelo 2 es más simple y más fácil de ajustar:
ajusta un modelo a un histograma de datos (existen tests de “bondad del ajuste” que
ajustado.
73
Figura III.63: Histograma experimental y modelo lognormal ajustado.
74
III.1.1 Los modelos de variograma
γ (0) = 0
γ ( h) ≥ 0
G G
γ (−h ) = γ (h )
Sin embargo, la teoría demuestra que estas condiciones no son suficientes: En efecto,
se puede probar que γ(h) pertenece a una familia Ω de funciones (llamadas funciones
de tipo negativo condicional) y la elección del modelo debe quedar restringida a esta
Por consiguiente, sólo hay que utilizar los modelos que describiremos a continuación (o
bien combinaciones de ellos, las cuales consisten en sumar dos o más modelos). No
a) El modelo potencia
75
γ(h) = ω⏐h⏐α
76
(h = ⏐h⏐)
El alcance es a y la meseta es C.
77
Figura III.66: A la izquierda la base de datos en una planta (75000 pozos de tronadura), a la derecha el
variograma esférico (en azul) y el experimental (puntos rojos). El alcance es del orden de 100 metros.
C(2(h/a) - (h/a)2) si h ≤ a
γ(h) =
C si h > a
78
d) El modelo exponencial: Crece más lentamente que el esférico o cuadrático y
γ(h) = C (1 - exp(- h / ω) )
Figura III.68: Modelo exponencial o modelo de Formery. Este modelo se presenta a veces en leyes que
están asociadas a fallas (ver la figura III.48, ver también la figura III.25). El alcance práctico es
aproximadamente iguala 3ω.
79
Figura III.69: Variograma vertical para el cobre, material quebrado remanente de la explotación por
block-caving en mina Salvador. En azul el modelo exponencial, en rojo, los puntos experimentales.
80
f) El modelo gaussiano: Tiene un comportamiento parabólico en el origen; su
ecuación es:
81
Figura III.72: Modelo cúbico. Es derivable en el origen.
82
a) Caso isótropo:
omnidireccional:
γ(h) = γomni(⏐h⏐)
En esta notación:
h = (hx , hy ) h = ( hx , hy , hz )
| h |= hx 2 + hy 2 h = hx 2 + hy 2 + hz 2
Ejemplo: En el caso estudiado en las figura III.34, el modelo isótropo ajustado es:
γ (h) = 0.417 | h |
83
γ (hx , hy ) = 0.417 hx 2 + hy 2
b) Caso anisótropo:
Anisotropía geométrica:
isótropo mediante una transformación lineal de las coordenadas. El caso más común
diferentes alcances:
84
En la figura se ha representado una anisotropía geométrica (en el caso isótropo lo
Sea k = a1/a2 > 1 la razón entre el alcance mayor y menor. Las fórmulas de
γ(h) = γ1(⏐h'⏐)
γ(h) = ω(θ)⏐h⏐
85
EW: γ1(h) = 0.0015|h| = ⏐h⏐ / 666.67
γ (h) = γ 1 (h ') = γ 1 ( hx 2 + k 2 hy 2 )
γ (h) = 0.0015 hx 2 + 12.96hy 2
Anisotropía zonal
En este caso la anisotropía no puede ser reducida por una transformación lineal simple
de las coordenadas.
86
Se define entonces el modelo de anisotropía zonal como un modelo anidado (o
imbricado); es decir:
G G
γ (hx , hy , hz ) = γ 1 (| h |) + γ 2 (| hz |)
En que:
87
G
| h |= hx 2 + hy 2 + hz2
88
Capítulo IV
El error de estimación
N = Nº de datos
En la práctica, se estima la ley media desconocida por una fórmula lineal del tipo:
89
zˆV = α1 z ( x1 ) + α 2 z ( x2 ) + ... + α N z ( xN )
α1 + α 2 + ... + α N = 1
αi = 1 / N ← media aritmética
αi = Si / S ← polígonos
ε = zˆV − zV
en que zˆV es la ley estimada (conocida) y zv es la ley real (desconocida).
Mencionemos que es equivalente definir el error como:
ε = zV − zˆV
90
Debido a que zv es desconocido, entonces ε es desconocido. Renunciamos entonces
Asumimos entonces que ε es una magnitud aleatoria; es decir, una variable aleatoria.
Esta magnitud aleatoria tiene una cierta ley de probabilidad caracterizada por una
esperanza matemática mE y una varianza σE2. Respecto de la ley de probabilidad del
La ley de probabilidad del error es la ley normal o de Gauss (en Geoestadística esta
aproximación es razonable).
91
Figura IV.3: Ley normal. σ es la desviación estándar.
Figura IV.4: Ley normal, regla de los dos sigma. Es el caso más utilizado en la práctica.
92
En otras palabras, utilizando probabilidades:
P(m - σ ≤ ε ≤ m + σ) = 0.68
P(m - 3σ ≤ ε ≤ m + σ) = 0.997
Estudio de mE:
En términos teóricos:
mE = E(ε)
Este valor es nulo porque los errores se compensan (siempre que el método de
estimación verifique la condición de insesgado Σα = 1).
Luego:
mE = 0
Estudio de σE2
En términos teóricos:
σ E2 = E ⎡⎣(ε − mE ) 2 ⎤⎦ = E ⎡⎣ε 2 ⎤⎦
σ E2 = E ⎡⎣ε 2 ⎤⎦
93
P[ - σE ≤ ε ≤ σE ] = 0.68 (i)
o bien, se puede afirmar, con 95% de confianza que la ley verdadera es igual a la ley
estándar σ E = σ E2 .
Cálculo de σE2 :
Sabemos que:
σ E2 = E (ε 2 ) = E ⎡⎣( zˆV − zV ) 2 ⎤⎦
94
Por otra parte:
N
zˆV = ∑α
i =1
i z ( xi )
1
zV =
V ∫ z ( x)dx
V
(La ley real desconocida se calcula, en el caso de ser posible, por el promedio de las
leyes de todos los puntos x dentro de V). La integral anterior puede ser simple, doble o
N
1
ε = ∑ α i z ( xi ) − ∫ z ( x)dx
i =1 VV
siguiente:
N
1 1 N N
σ 2
E = 2∑ α i ∫ i
γ ( x , x ) d x − ∫ ∫ γ ( x, y )dxdy − ∑ ∑ α α i j γ ( xi , x j )
i =1 V V V2 V V i =1 i =1
95
Figura IV.6: Vector que une los puntos p y q.
Antes de explicar cómo se calculan los términos en forma numérica, observamos que
σE2 depende de:
i) Término γ(xi, xj): Representa el valor de γ(h) siendo h el vector que une los
puntos xi y xj:
96
Figura IV.7: Vector que une el dato en xi con el dato en xj.
1
i) Término V ∫ γ ( xi , x) dx : Representa el valor medio de la función γ(h) (h es el
V
97
Entonces, la aproximación es:
1 1 k G
∫ γ ( xi , x)dx ≈ ∑ γ (h j )
VV k j =1
i) Si V es bidimensional: k ≥ 36 ptos.
1
iii) Término V2 ∫ ∫ γ ( x, y)dxdy
VV
: Representa el valor medio de la función γ(h)
98
Esta integral también se calcula por discretización de V en k puntos.
1 1 k Gk
V2 ∫V V∫ γ ( x , y ) dxdy ≈
k2
∑∑ γ (hij )
i =1 j =1
99
a) Si se estima zV por zV = 0.5 Z(x1) + 0.5 Z(x2) = 38, un calculo computacional
proporciona: 2 σE = 3.02
Luego:
zV = 38 ± 3.02
Luego:
zV = 35 ± 2.84
En este caso se utilizó un solo dato (contra 2 del caso anterior) y sin embargo se
utilizó mal la información de los sondajes, atribuyéndoles a ambos el mismo peso 0.5.
Es evidente que el dato z(x1) debe tener mayor peso que el dato z(x2).
Luego:
zV = 36.5 ± 2.38
Como este error (con confianza del 95%) es menor que los dos anteriores, se elimina la
paradoja. Este ejemplo nos lleva a la conclusión que no basta con tener datos sino que
hay que utilizarlos bien (es decir, ponderarlos en forma adecuada). Sin embargo, nos
100
queda una inquietud (que responderemos en el capítulo siguiente): ¿Existe otra
Observaciones:
N
1 1 N N
σ E2 = 2 ∑ α i ∫ γ ( xi , x ) d x − ∫ ∫ γ ( x , y ) d xd y − ∑∑α α i j γ ( xi , x j )
i =1 V V V2 V V i =1 i =1
101
Capítulo V
El Krigeado
V.1. Introducción
102
El nombre krigeado proviene de los trabajos de Daniel Krige en las minas de oro
sudafricanas de la Rand Corporation, en los años 50. La teoría fue formalizada una
Figura V.2: Dos bloques contiguos en la oficina Ossa (ver figura I.13). El color amarillo representa
mineral (con ley 1) y el blanco estéril (conley 0). Se observa que en ambos casos es necesario utilizar
datos externos al bloque.
El interés práctico más importante del krigeado, proviene, no del hecho que asegura la
mejor precisión posible, sino más bien porque permite evitar un error sistemático. En la
mayoría de los depósitos mineros, se deben seleccionar, para la explotación, un cierto
número de bloques, considerados como rentables y se deben abandonar otros bloques
considerados no-explotables. Daniel Krige demostró que, si esta selección se realizara
considerando exclusivamente las muestras interiores a cada bloque, resultaría
necesariamente (en promedio) una sobre-estimación de los bloques seleccionados. La
razón de este problema es que el histograma de las leyes reales de los bloques tiene
menos leyes extremas, ricas o pobres, luego tiene más leyes intermedias que el
histograma calculado con las muestras interiores, y, si se calcula el efecto de una
selección sobre este último histograma, los paneles eliminados serán en realidad
menos pobres que lo que se había previsto, y los paneles conservados menos ricos
(figura V.2).
103
Figura V.3: Histograma de bloques y de muestras.
zˆ K = λ1 z ( x1 ) + λ 2 z ( x 2 ) + ... + λ N z ( x N )
λ 1 + λ 2 + ... + λ N = 1
Como es natural, el krigeado atribuye pesos altos a las muestras cercanas a V y pesos
débiles a las alejadas. Sin embargo, esta regla intuitiva puede fallar en ciertas
104
situaciones en las cuales se habla de efecto de pantalla o de transferencia de
influencia.
i) λ1 ≥ λi (i = 2, 3, . . . , 8)
ii) λ 2, λ 3 ≥ λ i (i = 4, 5, . . . , 8)
iii) λ6 ≅ 0
iv) λ4 ≅ λ5 + λ6 + λ7 + λ8
105
Se dice que hay una transferencia de influencia; es decir, la influencia de la
muestra S4 se reparte entre las muestras S5, S6, S7 y S8. Se puede afirmar que
Para obtener las ecuaciones de krigeado hay que minimizar la expresión de σE2
N N N
1 1
σ 2
E = 2 ∑ λi ∫ i
γ ( x , x ) dx − ∫ ∫ γ ( x , y ) dxdy − ∑ ∑ λ λ γ ( x , x
i j i j )
i =1 V V V2 V V i =1 i =1
∑i =1
λi = 1
El método clásico para minimizar la expresión de σE2 (igualar a cero las derivadas
parciales de σE2 respecto de λ1, λ2, . . . , λN) no asegura que la suma de los λi sea 1.
En este caso hay que utilizar el método de Lagrange, el cual explicaremos con un
ejemplo:
Ejemplo:
106
A' = x2 + y2 - 2µ(x + y - 1)
1, entonces A' = A.
∂A' / ∂x = 2x - 2µ = 0
∂A' / ∂y = 2y - 2µ = 0
∂A' / ∂µ = -2(x + y - 1) = 0
x = 0.5
y = 0.5
µ = 0.5
107
Se demuestra que al realizar N + 1 derivaciones se obtiene el sistema de ecuaciones
siguiente:
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
λ1γ(xN, x1) + λ2γ(xN,x2) + ... + λNγ(xN, xN) + µ = (1/V)∫vγ(xN, x)dx
λ1 + λ2 + ... + λN = 1
que es un sistema lineal de N+1 ecuaciones con N+1 incógnitas (λ1, λ2, . . . , λN, µ)
Se demuestra que el sistema siempre tiene solución (se supone que el modelo de γ(h)
es autorizado), salvo el caso en el cual existen dos (o más) muestras diferentes con las
estimación.
108
N
1 1
σ 2
E = ∑λ
i =1
i
V ∫ γ ( x , x ) dx − V ∫ ∫ γ ( x , y ) dxdy + µ
i 2
V V V
Ejemplo 1:
Figura: V.5: Bloque a krigear, en negro la ley de la muestra, en azul el número de la muestra.
zˆV = λ1 z ( x1 ) + λ2 z ( x2 )
se obtiene:
λ1 = 0.77
λ2 = 0.23
2 σE = 2.37
109
Ejemplo 2:
zˆV = λ1 z ( x1 ) + λ 2 z ( x 2 ) + λ 3 z ( x3 ) + λ 4 z ( x 4 ) + λ 5 z ( x5 )
λ1 = 0.47
λ2 = 0.20
λ3 = 0.20
λ5 = 0.08
λ4 = 0.08
110
2σE = 0.12
datos 2 y 3 que a 4 y 5.
Krigeado Puntual
sistema siguiente:
con la varianza:
N
σ 2
E = ∑ λ γ (x , x
i =1
i i 0 )+ µ
111
Figura V.7: Krigeado del punto x0. Se puede generar una grilla de valores interpolados al hacer variar x0.
Esta técnica tiene aplicación en la cartografía automática y en la simulación de leyes.
que si se desea estimar la ley en un punto conocido (por ejemplo el punto A de la figura
V.7), el krigeado proporciona la ley del dato con una varianza σE2 = 0. Se dice que el
Figura V.8: Krigeado puntual en una dimensión versus interpolador de mínimos cuadrados.
112
Esta propiedad no la tienen otros interpoladores tales como los mínimos cuadrados.
a. Propiedad de simetría
Si γ(h) es isótropo, entonces datos que son simétricos respecto de V y con respecto a
En el ejemplo de la figura:
λ1 = λ1
λ2 = λ4 = λ6 = λ8
λ3 = λ5 = λ7 = λ9
Esta propiedad era útil cuando se resolvían los sistemas de krigeado “a mano”
113
b. Composición de krigeados
respectivos:
V1 z1 + V2 z2
zˆ =
V1 + V2
114
V.4. Vecindad de Estimación.
Figura V.11: Vecindad de estimación. Para el krigeado no importa la agrupación de datos al lado
izquierdo del bloque (el krigeado desagrupa la información). ¿Cuál radio tomar?
datos disponibles (krigeado completo). Sin embargo, esta situación implica cálculos
muy largos; por otra parte, las muestras alejadas tendrían un peso casi nulo. Por esta
razón la práctica recomienda restringirse a una vecindad de estimación que puede ser
una esfera o círculo, o bien un elipsoide o elipse (3D y 2D). Como recomendación
práctica, el radio de búsqueda en una cierta dirección no debe ser inferior al alcance en
esa dirección.
que contenga un promedio del orden de 8 muestras, los resultados son buenos. En el
115
Figura V.12: En el espacio 3D hay que elegir los parámetros de búsqueda de manera de que se
produzca “interpolación” entre los sondajes. En esta mina se observa que quedarán bloques mejor
estimados que otros y habrá que categorizarlos en medidos, indicados e inferidos.
Esta estrategia establece los parámetros que hay que utilizar para la búsqueda de
• Radios de búsqueda (Rx, Ry, Rz). En primera aproximación se pueden utilizar los
alcances del variograma en las direcciones (x, y, z), en una vecindad con forma
de elipsoide.
Figura V.13: Elipsoide de búsqueda. En algunas situaciones este elipsoide puede estar inclinado.
El centroide es el centro de gravedad del bloque.
116
• Mínimo k de muestras para krigear. Sirve para controlar el caso en que solo una
muestra cae en la vecindad. Si, por ejemplo, se pone k = 2, solo se krigearán los
I.8)
bloque.
Los parámetros l y s deben ser utilizados con precaución. Para no introducir artefactos,
es recomendable que estos valores sean altos, lo que hace que su uso no sea
117
Figura V.14: Existencia de sondaje inclinado: Al usar octantes con l = 1, se utilizan exclusivamente los
compósitos 1 y 5. Luego, hay que usar un l mayor.
Figura V.15: Al usar un máximo s = 1 de compósitos por octante, sólo se usan los compósitos 8 y 3
(casi a la misma cota que el bloque), sin tomar en cuenta la variación en la vertical. Luego, hay que usar
un s mayor.
118
VI. Bibliografía.
U. de Chile, 2007.
1996.
de Chile, 2003.
2005.
Ch. Huijbregts
P. Delfiner
www.cg.ensmp.fr/bibliotheque
119