Trabajo ganador del área Gestión y Economía Minera de CONAMIN 2023.

Por: Alfredo Marín Suárez, Universidad Nacional de Ingeniería.

Resumen

Este trabajo obtiene resultados algo similares a los que históricamente se lograron con la Simulación Condicional Multivariable, usando el método Bandas Rotantes, a la gran minería en diciembre 1978, en un yacimiento de hierro de la región Pilbara al noroeste del continente australiano, y plasmada en mi tesis de doctor ingeniero, desarrollada bajo la dirección del doctor George Matheron.

Para esto se transforma las leyes con Anamorfosis Inversa (nscore), se construye variogramas en todas las direcciones y luego se aplica la Simulación Secuencial Gaussiano. Finalmente se hace una transformación por Anamorfosis Directa para así poder emitir una serie de reportes de resultados para la operación minera, a partir de 30 a100 puntos simulados por cada bloque unitario (smu), tales como:

Mineral a concentradora, mineral a lixiviación, mineral a pad rom, mineral a stock y desmonte a botaderos. Mineral con variabilidad real controlada para la planta de tratamiento metalúrgico.  

Planos con ley simulada y variabilidad en cada bloque unitario (smu).

Planos de la probabilidad y ley media encima de un cut off requerido.

Planos del intervalo de confianza de la ley simulada por cada bloque unitario, planos de “errores”, es decir de la precisión de la ley simulada por cada bloque unitario, que permiten también clasificar los recursos recuperables teniendo en cuenta la selectividad minera en medidos e indicados, con su grado de confiabilidad porcentual. a diferencia de lo tradicional.

En suma, estos planos van a permitir una mejor definición de los polígonos con su tonelaje y ley media, así como realizar un Ore Control con mejor blending en el planeamiento de minado a corto plazo. Mostrado inicialmente en pruebas piloto en Cerro Verde, Las Bambas, Tantahuatay (2022-2023).

Una versión de la simulación presentada también se está aplicado a variables de geotecnia de la mina Antamina (2019-2023).

Estos resultados probabilísticos obtenidos conllevan al análisis de riesgo de proyectos mineros.

Se usó C++, Fortran, Python, y software auxiliar libre Sgems.

Introducción

Durante la planificación de la explotación de una mina, es necesario prever las dispersiones de las reservas a la salida de diversos procesos de extracción, almacenaje, etc. Si se conociera perfectamente el yacimiento minero esta tarea sería llevada sin mucha dificultad, pero esto no ocurre en la realidad, solo se conoce el yacimiento cuando este ya ha sido explotado.

Entonces, al no disponer de esta información se puede optar por simular el yacimiento, condicionando a los datos observados para darle mayor robustez a los resultados.

Para tal efecto usamos el modelo probabilístico propuesto por la geoestadística, que consiste en interpretar la repartición espacial de la ley de un mineral z(x) como una realización de una función aleatoria Z(x), x ∈ Rⁿ, caracterizado por sus dos primeros momentos. La simulación condicional consistirá en obtener varias realizaciones de la función aleatoria Z(x) que tengan el mismo variograma, y en el caso multivariable, los mismos variogramas cruzados y correlaciones de la realidad, condicionados por los datos experimentales z₀(X), es decir, que los valores simulados z(x) deben ser iguales a los valores experimentes z₀(X), en los puntos de observación.

El número de realizaciones z(x) del yacimiento simulado estará limitado solamente por las necesidades del usuario y por los recursos computacionales disponibles.

En este trabajo se ilustra la aplicación en un yacimiento tipo pórfido de Cu andino, cuya ley estudiada es el Cu total, yacimiento representativo de los grandes depósitos tipo pórfido, tales como Warintza en Ecuador; Cerro Verde, Toquepala y Cuajone en Perú, y Chuquicamata y Escondida en Chile, entre otros.

Estimación o simulación

Si nosotros representamos las leyes de los datos en la primera gráfica y luego en la segunda las leyes estimadas en los puntos X, tendremos lo mostrado en la Figura 1.

Se observa que con los valores estimados se suavizan las fluctuaciones de la variable ley, es decir, que con ellos no podemos cuantificar la variabilidad de la ley, perdiendo una propiedad tan valiosa. En cambio, la Simulación Condicional si reproduce las fluctuaciones con las características estadísticas y estructurales de la realidad, haciendo posible mejorar ampliamente la calidad del blending, planeamiento de minado, entre otras aplicaciones^{[1], [2], [4]}.

Base de la simulación geoestadística

Serie de números al azar azar (ver Figura 4).

Metodología de la simulación condicional

Análisis estadístico de datos

Se realiza el análisis estadístico como parte de la llamada ingeniería de datos que consiste en estudiar las leyes que caracterizan al tipo de yacimiento, información que debe encontrarse regularizada (compositada) proveniente principalmente de los huecos de perforación de la voladura (blast hole), aparte de la información de las leyes de los taladros de exploración mediante correlaciones, histogramas, detección de poblaciones anómalas y de grado de lognormalidad principalmente.

Decluster

Que consiste en desagrupar las leyes, afectando a las leyes un peso que depende del área o volumen que representa. Es denotar que también el Kriging hace esta función. 

 Proceso de Anamorfosis Inversa

ν Anamorfosis Inversa (nscore).

ν Consiste en transformar el valor de las leyes observadas en valores gaussianos.

ν La función de transformación en rigor es aproximada con una suma de polinomios de Hermite. Para verificar si la transformación ha sido correcta, construimos el histograma de valores transformados y saldrá gauss con media cero y desviación típica 1.

A continuación, se presenta la Anamorfosis Inversa de las leyes de mineral, como puede verse a cada ley real se resigna un valor gaussiano^{[1], [3], [4]}, como se observa en la Figura 3.

Variables regionalizadas y sus momentos

Se interpreta a las variables regionalizadas z(x) como realizaciones particulares de una función aleatoria (F.A.) estacionaria, aplicando la hipótesis estacionaria de orden 2.

a) m=E [Z(x)] media constante (momento de orden 1)

b) C(h)=E{Z(x).Z(x+h)} función covarianza

c) 2γ(h)=E{Z(x)-Z(X+h)}² función variograma

Que está relacionado con (b) mediante:

d) γ(h)=C(o)-C(h)  

Ejemplo de un variograma de las leyes transformadas por anamorfosis gaussiano, con una pepita y dos estructuras esféricas^{[1], [3], [4], [5]} (ver Figura 5).

Variografía 

Se realiza con las leyes que resultan de la aplicación de la Anamorfosis Inversa Gaussiana, en un dominio geológico del yacimiento. Consiste en construir los variogramas promedio en las diferentes direcciones del yacimiento con el fin de detectar la dirección de mayor alcance, es decir, la dirección del eje a₁ del elipsoide de variogramas, seguido de la definición del eje a₂ y a₃ del elipsoide. Cuando la modelación de los variogramas presenta dos o más estructuras se considerarán elipsoides adicionales por cada estructura encontrada.

Luego se determina los radios de un elipsoide de búsqueda, que permiten encontrar información para condicionar la simulación.

Estos permiten hacer Maching Learning, en el dominio geológico estudiado, a partir de los parámetros de los variogramas en las diferentes direcciones, es decir, de sus elipsoides de influencia^{[1], [3], [4], [5], [6], [7], [8]}.

Generación de valores simulados

A este nivel se puede usar la Simulación Secuencial Gaussiana, bastante simplificada y que básicamente consiste en lo siguientes pasos:

ν Se configura un camino aleatorio para cada nodo. 

ν Se usa un Kriging simple de estimación para definir una distribución gaussiana en el primer nodo, a partir de la data existente.

Se repite la estimación y muestreo Montecarlo para generar valor simulado secuencialmente para cada nodo en su camino aleatorio. De 30 a 100 valores por bloque unitario^[3].

A parte de este método de Generación de Valores Simulados, para obtener reportes en esta publicación, hay un primigenio histórico, con otro algoritmo de cálculo denominado Simulación Condicional Multivariable Bandas Rotantes, presentado en el anexo 1.

Anamorfosis Directa 

Finalmente transformamos la ley gaussiana simulada para cada bloque unitario a su valor real con la Anamorfosis Directa, tal como se muestra en la Figura 8^{[1], [3], [4], [6], [7], [8]}.

Aplicaciones de la simulación

Reportes de Ore Control a obtener: 

1. Mineral a concentradora ◊ Chancadora millchan.

2. Mineral a lixiviación ◊ Chancadora hidrochan.

3. Mineral a pad rom.

4. Mineral a stock.

5. Desmonte a botaderos.

En las Figuras 10, 11 y 12, presentamos un plano con los blasthole de la operación minera que permite la definición del contorno poligonal, para calcular su ley promedio y tonelaje de los minerales mencionados.

Planos de trabajo obtenidos a partir del número de puntos (de 30 a 100) simulados por bloque unitario

Los valores que aparecen en cada bloque unitario de los planos de las Figuras 13 a la 15, son producto de estadísticas de la elección de 30 hasta 100 valores simuladas condicionadas.

Caso de verificación visual de la simulacion condicional

Este caso se grafica en la Figuras 16 y 17. 

Intervalo de confianza de la simulación

Se calcula el intervalo de confianza de la ley simulada, es decir, la precisión de la ley por cada bloque, cálculo que va a permitir también clasificar los recursos medidos e indicados, con su grado de confiabilidad porcentual a diferencia de lo tradicional. 

En suma, estos planos van a permitir una mejor definición del perímetro de la explotación local con su tonelaje y ley, así como realizar un Ore Control con mejor blending para el planeamiento de minado a corto plazo^{[2], [4], [9], [10], [11]} (ver Figuras 18 y 19).

Conclusiones

1. Como se puede observar, los reportes a corto plazo del día a día de Ore Control obtenidos por la simulación condicional, dan más información que cualquier técnica de estimación incluyendo el mejor estimador como lo es la técnica de Kriging, por lo tanto, la simulación condicional da resultado novedosos y más eficientes para la minería a cielo abierto, dando recursos recuperables teniendo en cuenta la selectividad minera en medidos e indicados.

2. Existen más aplicaciones, así los resultados probabilísticos obtenidos conllevan al análisis de riesgo de proyectos mineros. 

3. Se ha iniciado la aplicación con una prueba piloto (2022-2023), a través de un curso taller de geoestadística con Simulación Condicional, aplicada al Ore Control con Blending en el planeamiento de minado a corto plazo, en las minas de cobre de Cerro Verde y Las Bambas, igualmente se ha iniciado su aplicación en el yacimiento de oro de Tantahuatay-Coimolache de Compañía de Minas Buenaventura.

4. Una versión de la simulación presentada también se está aplicado a variables de geotecnia de la Mina Antamina (2019-2023).

5. Es de notar que la Simulación Condicional presentada es una técnica de Maching Learning por excelencia.

6. El software fue construido en Fortran y Python, con la ayuda de software libre Sgems, para los cálculos complementarios y presentaciones gráficas.

Referencias

1. Matheron G. 1962, 1963. Traité de Géostatistique Appliquée. Ed. Technip, Paris VOL. 1; VOL. 2

2. Guibal D. 1972. Simulation de Schémas Intrinsèques. N-291 E.N.S.M.P.

3. Maréchal A., Deraisme J., Journel A., Matheron G. 1978. Cours de Géostatistique non Linéaire. C-74 E.N.S.M.P.

4. Marín Suárez A. 1978. Méthodologie de L'estimation et Simulation Multivariable des Grands Gisements Tridimensionnels. Thèse présentée à I'école Nationale Supérieure des Mines de Paris.

5. Marín Suárez A. 1986. Modelo Geoestadístico de Filones de Almadén. Ed. Minas de Almadén S.A., Almadén (España).

6. Marín Suárez A. 2004. “Estimación de reservas de la mina Toromocho”, para David Lowell Exploration”. Perú.

7. Lowell Mineral Exploration. 2005. “Estimación geoestadísticos de reservas de Warintza”, Ecuador.

8. Lowell Mineral Exploration. 2006. “Estimación de recursos del yacimiento de Titanio en Lateritas” Paraguay.

9. Lowell Mineral Exploration. 2007. “Estimación geoestadística de recursos del proyecto Toromocho” Perú (2007). 

10. “Estimación de Recursos del Proyecto Minero Michiquillay”, Perú. 2017. Para Activos Mineros del Ministerio de Energía y Minas, Proinversión y el Fonafe.

11. “Simulación de Variables Geotécnicas por Simulación Condicional en la mina Antamina”, Perú. 2020, 2021 y 2022.

Anexo 1. Fundamentos de la Simulación Bandas Rotantes

a) Simulación del valor de una ley en R³ a partir de R¹

Previamente hay que definir las dimensiones del bloque unitario escogido. 

Sea: {D_i,i=1 a n} rectas en el espacio R³,cada recta D_i queda definido por el vector unitario u_i

Se genera sobre cada recta y de manera independiente, realizaciones de la F.A. Y_i(u), a una dimensión con:

E[Y_i (u)]=0

E[Y_i (u) Y_i (u+h)]=C₁ (h)

Luego si:

x=(X₁, X₂,…,X_n) es un punto de Rⁿ y estudiamos la F.A. Y(x) este podrá ser expresado a partir de los Y_i (u) obtenidos anteriormente.

Como se ve en cada recta tenemos una banda que rotará (1), (4).

b) Visualización de los resultados simulados no condicionales, mostrados en tres etapas.

ν La primera etapa es mostrar en el espacio los valores simulados Y(x) con el método de las Bandas Rotantes.

ν “Ocultar” el valor simulado Y(x) en el punto x.

ν Estimar por Kriging el valor de Y*(x) en el punto x, a partir de los valores simulados, que se encuentran en su elipsoide de influencia determinados por los variogramas.

c) Condicionado

Siendo el objetivo central reponer la variabilidad de la variable ley, estimado por Kriging, calculamos el incremento de ley ε, afectando a cada ley estimada este incremento, con lo cual lograremos reponer la pérdida de variabilidad.

En términos matemáticos lo expresamos así: 

Tenemos que calcular el error ε en la expresión

Y_Sc (x)= Y_K^* (x) + ε

Siendo: ε = [Y_S(x) - Y_Sk^*(x)]

Con lo cual tenemos:

Y_Sc (x)= Y_K^*(x) + [Y_S(x) - Y_Sk^*(x)]

Definidas para la variable ley en un pequeño bloque unitario v

Donde:

ν Y_Sc(x) representa a la variable ley simulada en forma condicional.

ν Y_S(x) es la variable simulada (no condicional).

ν Y_K^*(x) es el valor estimado por Kriging de Y(x), sobre un pequeño bloque unitario centrado en x.

ν Y_Sk^*(x) es el valor estimado por Kriging de la variable simulada Y_S(x), en el punto ¨ocultado¨, sobre el mismo bloque unitario.

Este procedimiento está explicado en la Figura 4A1.

La relación anterior podemos escribirla así:

Y_Sc(x) = Y_S(x) + [Y_K^*(x) - Y_Sk^*(x)]

Los ponderadores de las estimaciones por Kriging a media conocida nula, de Y_K*(x) y de Y_Sk^*(x) son las mismas con una configuración geométrica de los datos también igual. Entonces es suficiente realizar solo un Kriging de diferencias entre Y_α y Y_Sα: valor real y valores simulados gausianos. Entonces la expresión anterior se transforma en:

Y_Sc(x) = YS(x) + λ^α [Y_α - Y_Sα]

Donde λ^α son los ponderadores del sistema de ecuaciones del Kriging aplicados en cada punto alfa.

Siendo el resultado obtenido Y_Sc(x) un valor simulado de la variable ley para un bloque unitario.