DOCTORA L T H E S I S

Computational Methods and Strategies

for Geometallurgy

Pierre-Henri Koch

Mineral Processing
Computational Methods and Strategies
for Geometallurgy

Pierre-Henri Koch

May 13, 2019

Luleå University of Technology

Department of Civil, Environmental & Natural Resources Engineering
Division of Minerals and Metallurgical Engineering
Printed by Luleå University of Technology, Graphic Production 2019

ISSN 1402-1544
ISBN 978-91-7790-346-8 (print)
ISBN 978-91-7790-347-5 (pdf)
Luleå 2019
www.ltu.se
Abstract

At the interface of geology and mineral processing, geometallurgy is a powerful tool for en-
hancing resource efficiency. A spatial model that represents the ore body in terms of miner-
alogy and physical properties is combined with a process model that describes the concen-
tration process. The performance of a given ore in the process is computed in terms of grade
and recovery of the mineral of interest in the concentrate, but also the presence of potential
penalty elements and energy costs. The inclusion of ore performance indicators in a block
model yields a geometallurgical model that considers the variations in an ore body.

Progress has been made in recent years to list and study different processing options in
terms of data requirements and implementation costs. While providing useful data, little
advance was made to guide decision-making and to handle uncertainty. The objective has,
therefore, been to develop, choose and validate computational methods that suggest optimal
decisions in the scope of geometallurgical strategies for an iron ore and a porphyry copper
deposit.

The selected approach is based on an analysis of structure and regularity from the ore block
down to the mineral grains. By selecting the appropriate mathematical tool for each scale,
the dimension of the data is reduced and the different scales are then taken into account in
making decisions. Methods introduced for dimension reduction include machine learning
models, statistical models and spectral descriptors. The decision models rely on stochastic
multi-armed bandits which are a form of reinforcement learning. The presentation of the
different models proceeds by zooming in from coarse scale to fine scale then taking a step
back and analyze the implications. Data that was collected during sampling campaigns and
industrial plant surveys is used to design and verify the proposed models.

i
 With regard to the dimension reduction problem, results showed the method’s ability to
classify mineral textures and identify mineral phases with more than 90 percent accuracy on
the selected data sets of optical images and incorporate different physical properties into a
geometallurgical ore type classification. Decision results showed that strategies in the case of
a feed grade control or when different ore types were identified, resulted in a twofold increase
of a reward function which is either Boolean (the product fulfills quality requirements or
not), or continuous (an economic objective). The cumulative value of the reward function
measured the optimality of a processing strategy. Quantitative methods were introduced to
evaluate ore classification as well as geometallurgical strategies.

The achieved results suggest the introduction of these computational methods in the prac-
tice of geometallurgy. The increased knowledge of different ore type performances and appro-
priate models lead to optimal decisions for improved resource efficiency along the ore value
chain. This is achieved by both maximizing profit and decreasing environmental impact, for
example by choosing processing routes that minimize energy consumption.

ii
Key words

Geometallurgy, process simulation, minerals, machine learning, textures, decision-making

iii
Acknowledgments

I wish to thank my supervisors for their support in both practical and technical aspects. Pro-
fessor Jan Rosenkranz for keeping this work on track by interesting discussions and for letting
me the freedom to explore computational methods from various horizons. Pertti Lamberg
for his ideas and previous work on liberation and simulation. Cecilia Lund for her friendship
and support, expertise in geology, and daily discussions.
This study is part of PREP (Primary resource efficiency by enhanced prediction), a project
within the Swedish Mining Research program SIP-STRIM and partially funded by VINNOVA. I
am grateful to our partners LKAB, Boliden, Outotec, Lundin Mining and Chalmers University
of Technology. The help of Antti Remes and Antti Roine from Outotec, Therese Lindberg
and Mattias Gustafsson from LKAB and Adam McElroy, Iris Wunderlich from Boliden was
particularly appreciated.
All my colleagues from the MiMeR department from Luleå University of Technology were also
a powerful support. I would like to thank in particular Viktor Lishchuk for his ideas about
PREP and his work in spatial modeling as well as his help in characterizing the samples we
collected and Mehdi Parian for providing interesting discussions about textures and provid-
ing data. Bertil Pålsson for his experience in mineral processing and ideas about simulation is
also kindly thanked. I also thank Malin Johansson for her full support and Thursday’s lunches.

I would like to thank my Mum for her language skills and availability, and my Dad for his
organization and sense of humor. Mira is thanked for her careful review of my thesis.

I would like to thank my wife Coralie for her incredible patience and constant support.

iv
Contents

Abstract i

Key words iii

Acknowledgments iv

List of papers and tools vii

Abbreviations x

1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Aim and objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Structure of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2 Selection of the case studies 10

2.1 Aitik copper ore deposit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Leveäniemi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3 Methods, concepts and tools 18

3.1 Drill core scanning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 Geometallurgical tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.3 Methods for texture quantification and classification . . . . . . . . . . . . . . . . 19
3.4 Sequential decisions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

4 Numerical results and discussion 33

4.1 Geometallurgical classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

v
 4.1.1 Mineralogical characterization . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.1.2 Grindability tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.2 Texture quantification and classification . . . . . . . . . . . . . . . . . . . . . . . 38
4.3 Data fusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.4 Sequential decisions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

5 Conclusion 56

Bibliography 60

Papers 69

vi
List of papers and tools

List of papers included in this thesis

• PAPER I : Koch, P.-H., Lund, C., and Rosenkranz, J. (2019a). Automated drill core
mineralogical characterization method for texture classification and modal mineralogy
estimation for geometallurgy. Minerals Engineering, 136:99–109

• PAPER II :Koch, P.-H., Tiu, G., Semsari, P., Parian, M., and Lishchuk, V. (2019c). Textures
in geometallurgy : from characterization to applications. pages 1–28. Submitted to
Minerals Engineering

• PAPER III :Koch, P.-H., Lund, C., and Rosenkranz, J. (2019b). Data fusion for enhanced
prediction of geometallurgical variables. Submitted to Applied computing & Geosciences

• PAPER IV :Koch, P.-H. and Rosenkranz, J. (2019). Sequential decision-making in

geometallurgy. pages 1–22. Submitted to Minerals Engineering

List of papers not included in this thesis

• Koch, P.-H., Lamberg, P., and Rosenkranz, J. (2015). How to build a process model in a
geometallurgical program? In Andre-Mayer, A., Cathelineau, M., Muchez, P., Pirard, E.,
and Sindern, S., editors, Proceedings of the 13th SGA Biennial meeting, page 85, Nancy,
France

• Koch, P.-H. and Rosenkranz, J. (2017). Texture-based liberation models for comminu-
tion. In Proceedings of the Conference in Minerals Engineering, pages 83–96, Luleå,
Sweden

vii
 • Koch, P.-H. (2018). A numerical study of the effects of microwave pre-treatment on
value liberation from a zinc ore. In Proceedings of the 29th International Mineral Pro-
cessing Congress, Moscow, Russian Federation

• Koch, P.-H., Lund, C., Lishchuk, V., Lamberg, P., and Rosenkranz, J. (2017). Automated
classification of drill cores textures for geometallurgy. In Proceedings of Process Miner-
alogy 2017, Cape Town, South Africa

• Lishchuk, V., Koch, P.-H., Lund, C., and Lamberg, P. (2015a). The geometallurgical
framework. Malmberget and Mikheevskoye case studies. In Proceedings of the XVth
Conference of PhD students and Young Scientists, pages 57–66, Szlarska Poreba (Poland)

• Lishchuk, V., Lund, C., Koch, P.-H., Pålsson, B. I., and Gustafsson, M. (2018b). Geomet-
allurgical characterisation of Leveäniemi iron ore – Unlocking the patterns. Minerals
Engineering, 131:325–335

• Guntoro, P. I., Ghorbani, Y., Koch, P.-H., and Rosenkranz, J. (2019). X-ray Microcom-
puted Tomography (µCT) for Mineral Characterization: A Review of Data Analysis Meth-
ods. Minerals, 9(3):183

Contribution of the author to the content of each paper :

• PAPER I : Concept, code for acquisition, texture and pixel classification. Sampling,
characterization and verification of the case study. The paper was mainly written by
the author, reviewed and complemented by the co-authors.

• PAPER II : Concept, parts on ergodicity and stationarity including figures, coordina-

tion of the different cases. The paper was mainly written by the author, reviewed and
complemented by the co-authors.

• PAPER III : Concept, design of the networks’ architecture in TensorFlow and Keras,
sampling and characterization of the samples to the data set, training and testing the
results in Python. Hyperparameters optimization tests with Talos.The paper was mainly
written by the author, reviewed and complemented by the co-authors.

• PAPER IV : Implementation of all the algorithms, sampling and characterization of the

samples for the case studies. The paper was mainly written by the author, reviewed and
complemented by the co-authors.

viii
Tools and models produced during this thesis :

• A MATLAB framework to analyze and simulate stationary and non-stationary textures

based on Potts model.

• Process unit models in .NET implemented in releases of HSC Sim for blasting and
grinding

• Segmentation models in WEKA for drill core scanning and mineral identification

• Deep convolutional neural networks in Python with TensorFlow/Keras, Talos, scikit and
sklearn for classification and regression of geometallurgical variables

• A MATLAB framework for Bandits algorithms in the context of geometallurgy

• Process models for Aitik’s flotation circuit and Malmberget in HSC 9

ix
Abbreviations

Type of deposits

IOCG Iron Oxide Copper Gold

AIO Apatite Iron ore

Elements

Au Gold
Cu Copper
Fe Iron
Ti Titanium
Mo Molybdenum
Pb Lead
V Vanadium
Zn Zinc

Minerals

Ap Apatite
Amph Amphibole
Bt Biotite
Chl Chlorite
Cpy Chalcopyrite
Ep Epidote
Fsp Feldspar
Grt Garnet
Hbl Hornblende
Hem Hematite
Ilm Ilmenite
K-Fsp Potassium Feldspar

x
Mgt Magnetite
Ms Muscovite
Na-Pl Sodium Felsdpar
Pl Plagioclase
Po Pyrrhotite
Px Pyroxene
Py Pyrite
Qtz Quartz
Ti-Mgt Titano-Magnetite
Ttn Titanite
Tur Tourmaline
Zrn Zircon

Methodology

BWi Bond Work index

GCT Geometallurgical Comminution Test
(S)AG (Semi)Autogeneous Mill
EMC Element to mineral conversion
SEM Scanning Electron Microscope
BSE Back-scattered electrons
EDS Electron Dispersive Spectroscopy
XRD X-ray Diffraction
XRM X-ray Microscope
FFT Fast Fourier Transform
µXCT X-ray Computed Micro-tomography

Algorithms and data

LS Least-squares
kNN k-nearest neighbors
ANN Artificial Neural networks
CNN, ConvNet Convolutional Networks
SGD Stochastic gradient descent
i.i.d Independent identically distributed
S- Stationary

xi
NS- Non-stationary
σ Standard deviation
CI Confidence interval
MAB Multi-armed bandits
UCB Upper confidence bounds
KL Küllback-Leibler divergence
Adam Adaptive moment (stochastic gradient descent)
SW Sliding window

Image processing

px Pixels
vx Voxels
[] Unitless
BW Black and White (Images)
RGB Red, Green and Blue Channels (Visible-light color Images)
LBP Local Binary Patterns
GLCM Gray-level co-occurrence matrix

xii
Chapter 1

Introduction

1.1 Background

Resources Emissions

Mine

Value Capital Ore Data

Value
Society Processes
Capital

Emissions Resources Emissions Resources

ENVIRONMENT

Figure 1.1: Ore value chain

Geology and metallurgy combined, form a multi-disciplinary approach that links economic
geology, process mineralogy and metallurgy to establish and validate a predictive model of the

1
ore value chain performance. The proposed definition distinguishes characterization from
modelling and accounts for current research activity in flexibility, sustainability, and variabil-
ity of mineral processing circuits (Koch, 2017, Lamberg, 2011, McQuiston and Bechaud, 1968,
Vann et al., 2011).This way of approaching mining and mineral processing is not new but the
combination of these disciplines into a unified framework, that can be used in practice to im-
prove the creation of value and minimize the impact on environment, requires appropriate
instruments, data, software and interpretation of the results. The industrial implementation
of geometallurgy in a raw materials’ context is called a geometallurgical program. Geomet-
allurgy reflects past and current research at universities, research centers and in the mineral
industry. Predictive models are preferred to explanatory modelling since predictions offer a
possibility of action before the change occurs. The final output of such a program is called a
geometallurgical block model (Aasly and Ellefmo, 2014, Alruiz et al., 2009, Lund et al., 2013,
Schouwstra et al., 2013) . In order to integrate different scales and aspects of the ore value
chain, three sub-models exist.

1. Spatial model: 3D block model in which geometallurgical data is available. This model
is used for data management, visualization, for production planning and to test differ-
ent extraction scenarios.

2. Process model: a flow sheet or arithmetic model in which geometallurgical data is

used to forecast the products. Depending on the approach, this model can predict
the recovery rate, the element grade or mineral grade in the concentrate and tailings,
particle size distribution, water, and energy consumption.

3. Economic model: a set of equations that take into account the OPEX and CAPEX for
both the mine and the concentrator, the commodity price, the revenue and uncertainty
in order to calculate the available value contained in a block of the spatial model at a
given time (Runge, 1998). This model relies on the previous models to obtain its initial
data. Once the value has been calculated, it can be placed back to the spatial model to
help planning and production.

Several approaches exist in geometallurgy for characterizing an ore within the different sub-
models. The most common ones found in literature can be classified into three main types.

2
 • The traditional approach, where both the spatial and process models are based on
elemental grades from drill cores assays.

• The proxy approach, where the process model is based on many small-scale tests that
provide proxy variables to describe the metallurgical response of the ore in the process.

• The mineralogical approach, where the spatial and process models are based on quan-
titative mineral characterization.

In this work, the mineralogical approach was selected (Koch, 2017). Minerals, or more pre-
cisely mineral particles encapsulate enough information to generate the elemental composi-
tion and study beneficiation processes (Lamberg, 2011). Moreover, it is possible to generate
a population of mineral particles by simulation and use it as a stream for process simulation.
This is made possible by the stream structure available in some simulation systems, for in-
stance HSC Sim (Outotec, 2019). Depending on the level of detail required, one can choose
to work with either

• One-dimensional data : mineral phases, chemical and modal composition but no size
information

• Two-dimensional data : chemical and modal composition by size fractions

• Three-dimensional data : chemical and modal composition and liberation distribution

by particles size fractions

The value of geometallurgy resides in its capacity to cover the entire value chain of min-
erals with error estimates for every step and its ability to give more reliable predictions. The
first step is then to acquire information by either measuring the physical properties of the
intact ore (on drill cores for example), particles obtained from crushing and grinding the ore,
or samples from the plant. Once a sufficient amount of data has been collected, it can be used
to improve the grade or recovery of the mineral of interest, optimize the throughput or reduce
the environmental footprint of the process. In a geometallurgical program, a large amount of
data from different sources is acquired, structured and made available as a basis for decision-
making for mining and process planning. As such, the last step is to take decisions that will
either prevent a degradation of the performances or maximize them.

3
1.2 Motivation

Geometallurgy as a field that combines geology, mineralogy and metallurgy is limited by

the available data and the numerical models that are used. Conceptually, the links between
the spatial models, the process models and the appropriate decisions are currently being
investigated, see Figure 1.2.

Building geometallurgical block models has been a challenge since the information avail-
able in a geometallurgical program is not always additive. Moreover, the data itself can be
numeric or categorical. The validity of the theory of kriging and the degree of validation
of multiple-point geostatistics or machine learning methods has recently been studied by
Lishchuk et al. (2018a). Drill core information is commonly used to populate block models.
However, few automatic scanning systems are available. The detection and classification of
mineral phases is often done by geologists who need to rely on visual inspection to classify
the ore. An automated system based on visible light and machine learning allows supporting
human expertise but can also be trained and improve over time.

Process models have been developed similarly to the simulators used in chemical engineer-
ing. As such, they were not initially designed for particulate solid material. Here, no chemical
or phase equilibria are used. In addition to that, the solids properties are distributed, e.g.
size distribution, size dependent composition etc., thus requiring a stream data structure
adapted to particle flows. This resulted in having empirical or semi-empirical models for
comminution and concentration processes. Blasting, crushing and grinding aim at reduc-
ing the particle size and as a result increase the degree of mineral liberation of the particle
population. Depending on the ore properties, a high degree of liberation can be difficult
to reach, due to the texture of the ore. The question of ore texture has recently received
more attention due to the increased complexity of ores being mined but also to technological
progress (Bonnici, 2012, Donskoi et al., 2008, Pérez-Barnuevo et al., 2018, Tungpalan et al.,
2015). While qualitative descriptions of textures have been widely used in geology, quantita-
tive descriptions are still not commonly used for predicting production. It has been observed
that textures have an impact on the two main operations of mineral processing operations, i.e.
in comminution (Barbery, 1991, Barbery and Leroux, 1988, Gay, 1999, 2004, King, 1979, Koch
and Rosenkranz, 2017, Pérez-Barnuevo et al., 2013) and concentration (Parian et al., 2016,

4
Tungpalan et al., 2015, Vos, 2017). It can be assumed that when more data is available, pro-
cess simulations will better reproduce the complex behavior observed in mineral processing
plants. Many methods are now available to collect chemical, mineralogical but also geomet-
ric parameters from both intact ore and particles. Table 1.1 shows the different techniques
available to obtain textural information from a sample. One can see that each method oper-
ates at a different level of detail but also at different scales. The question of extracting relevant
information and combine the data from different sources remains a challenge (Lawson, 2017,
van den Boogaart and Tolosana-Delgado, 2018).

Table 1.1: Ore characterization methods

Properties Sub-category Data collected

Chemical List of chemical elements

Geological Mineralogical List of minerals

Physical e.g., hardness, density
Breakages Energy, PSD, liberation
Concentrations Recovery, mass pull
Processing
Qualities Grades, penalties
Environmental Acid rock drainage
Qualities Recoveries, mass pull
Metallurgical
Environmental CO2 emission
Visual e.g., color, luster, fluorescence
Reflective
Multispectral Mineral phases
Spatial Location of origin Coordinates x, y, z
Mechanical Rock quality designation RQD
Conductivity Conductivity/Resistivity
Geophysical
Susceptibility Susceptibility

Despite a growing literature corpus on geometallurgy, few studies have focused on decision-
making based on the data gathered within a geometallurgical program. While human-based
decisions can always be taken depending on an ore classification or process mineralogy
knowledge, no automated tool offers suggestions for action using the results. Given the costs

5
that a geometallurgical program represents as well as the change of culture it implies, one
would like to maximize the impact of both the spatial and process model. As pointed out
earlier in Koch (2017), the real value of such programs lies in the ability to predict and take
appropriate actions. Decision models move geometallurgy from observations and model
based predictions to actions. So far, the use of comprehensive data to take optimal decisions
is not yet supported by automatic decision support systems but is instead still done manually
outside the plant automation system.

Ore body Ores Processing plant Products Market

Predictions Predictions

Predictions

Samples

Ore Particles Metallurgical Market

characterization characterization characterization analysis

? ?
Paper 1 : Drillcores Paper 3 : Data fusion Paper 4 : Decisions

Spatial model Blocks Process models Performance indicators Economic model

Predictions
Samples

Information

Main objects

Figure 1.2: Full workflow in geometallurgy. Question marks indicate challenging areas

Even though characterization methods exist for each step, the links between characteriza-
tion and the models used for simulation are still largely missing, compare Figure 1.2. In partic-
ular, the difference of scale between measurements makes it difficult to validate one method
with another. One challenge with the models themselves is the delay between sampling and
characterization which makes the calibration of highly detailed models complicated in the
time frame required for short and mid-term decisions. The value of the data could be con-
siderably increased by collecting relevant data in large amounts, combining different sources
and use these to guide decisions.

6
1.3 Aim and objectives

The main hypothesis underlying this thesis is that texture characterization, results from ge-
ometallurgical tests, together with mineralogical analyses, provide quantitative data that can
be used to take optimal decisions with regard to processing options in geometallurgy. Based
on that, the objectives and research questions of this study can be derived.
The first objective is to identify, test and apply the appropriate tools to quantify mineral
textures. In particular, the following questions are addressed

• What is a texture in the context of mineral processing and geometallurgy?

• Which mathematical tools are available and relevant to quantify textures and discrimi-
nating between them?

• How to use textural information?

The second objective is to design a recommendation system that uses the available infor-
mation to choose from different processing routes in a sequential context. In particular, the
following questions are addressed:

• When a mass of ore is to be processed, how is variability distributed between the spatial
model and the process model?

• How can the available information be used to support decisions?

• How can decisions and strategies in geometallurgy be quantitatively compared?

7
1.4 Structure of the thesis

Ore body Ores Processing plant Products Market

Predictions Predictions

Predictions

Samples

Ore Particles Metallurgical Market

characterization characterization characterization analysis

Paper 1 : Drillcores Paper 2 : Textures Paper 3 : Data fusion Paper 4 : Decisions

Spatial model Blocks Process models Performance indicators Economic model

Predictions
Samples

Information

Main objects

Figure 1.3: Entire workflow in geometallurgy - Filling the gaps

This thesis addresses the challenges identified by both academics and the industry. Results
are organized in four main areas presented in Figure 1.3.

1. Drill core scanning Paper 1 offers a novel method to identify and estimate mineral
grades from visible light drill core scanning. The method is tested on a porphyry cop-
per deposit. Results indicate that such a system can help reduce costs and improve
accuracy. The data can then be added to the spatial model or help the process model.

2. Textures: Paper 2 generalizes the concept of texture and introduces fundamental con-
cepts for studying and modelling of a mineral map. Examples of relevant industrial
applications in copper, iron, and lithium industry are presented. These tools help se-
lect relevant information from textures.

3. Data fusion: Paper 3 introduces neural networks with an appropriate architecture to

accept various kinds of data, among these also textural information, and make predic-
tions regarding the metallurgical performance of the ores. This class of models allows
both prediction but also regression in presence of noise in the measurements. The
models also offer flexibility and dimension reduction.

8
 4. Optimal decisions: Paper 4 presents an approach that helps operators take optimal
decisions in different settings when variability is present, supported by the available
data or even in absence of data. By interacting with an economic model, the impact of
decisions on the economic performance of an operation is quantified.

The aforementioned research questions are approached by following a multi-scale ap-

proach. First, the available data is presented, then the relevant tools for quantitative charac-
terization are introduced at the mineral map level. Once the finest level of detail has been
reached, the subsequent chapters will broaden to the decisions one can take based on the
textural information. Finally, challenges and future opportunities that come with the new
tools and concepts developed are discussed.

9
Chapter 2

Selection of the case studies

Within this work, two industrial case studies have been studied. The first one is a porphyry
copper deposit located in the north of Sweden from which copper and gold are concentrated.
The second one is an iron oxide ore deposit from which magnetite is extracted. This section
presents a summary of the geological features as well as the main ore mineralogy.

2.1 Aitik copper ore deposit

Aitik is a porphyry copper-gold deposit with an overprint that exhibits IOCG-like features
(Wanhainen et al., 2012). Its structure is characterized by strong deformations and metamor-
phic features. This Norrbotten Cu-Au-Fe oxide district shares a number of characteristics in
terms of intrusions, mineral assemblages and alteration (Wanhainen et al., 2006, Wanhainen
and Martinsson, 1999).
The Aitik mine is currently operated by New Boliden AB, resources and reserves are pub-
licly reported on an annual basis, see Table 2.1.

10
 Table 2.1: Resource and reserves at Aitik (Boliden, 2018)

Reserves Quantity 2018 [kt] Au [g/t] Ag [g/t] Cu [wt%]

Proven 787000 0.15 1.2 0.22

Probable 361000 0.13 1.2 0.23
Resources

Measured 204000 0.08 0.8 0.15

Indicated 1127000 0.09 0.8 0.17
Inferred 175000 0.11 0.5 0.14
Production

Milled ore 38472 0.13 1.98 0.28

Concentrate 404 24.58

Figure 2.1: A : General location of both Aitik and Leveäniemi, B : Northern Norrbotten ore
province, Leveäniemi is located in the Gruvberget area, C : Geological map of Aitik (Wan-
hainen et al., 2003)

11
 The current genetic model has been developed by Wanhainen (2005) and consists of the
following steps:

1. The Haparanda suite formed in a context of oceanic crust subduction under the Ar-
chaean craton ( 1.9 Ga).

2. Early mineralization occurred as a result of the alteration of a monzodiorite intrusion

around 1.88 Ga. The alteration pattern and grade distribution follow the porphyry
copper style (low-grade inner zone to high-grade outer zone + stockwerk).

3. From 1.8 Ga to at least 1.73 Ga, IOCG-type hydrothermal fluids circulated through the
system (overprinting). This is related to presence of magnetite, garnet, quartz, REE,
epidote, tourmaline, scapolite, K-feldspar. This later event seems to have a limited
economic impact on the deposit.

Figure 2.2: Cross-section of Aitik deposit, A : Intrusion, B : Current location with the open pit
(Wanhainen et al., 2003)

12
The footwall consists mainly of quartz/micro-quartz monzodiorite and feldspar-biotite-amphibole
gneiss. The ore zone consists of garnet-biotite schist (footwall side) and quartz-muscovite-
sericite schist (hanging wall side). The hanging wall consists of feldspar-biotite-amphibole
gneiss. Throughout the deposit, pegmatite dykes can be found. Mineralization occurs mostly
as disseminated sulphides or veinlets.

2.2 Leveäniemi

The massive magnetite mineralisation at Leveäniemi is hosted by metasediments and metavol-

canic rocks. Magnetite occurs as massive veins and replacements of these rocks. The main
alterations present are tremolite, biotite, potassium feldspar, scapolite and magnetite.

Figure 2.3: Left : Geological map of Leveäniemi. 1:Magnetite ore, 2:Calcite-rich magnetite ore,
3:Hematite-altered ore, 4:Ore breccia, 5:Leptite, 6:Sericite schist, 7:Metabasite, 8:Lina granite,
9:Skarn (Frietsch, 1966). Right: Profile of the ore body, the solid line indicates the boundary
of the mineralization. (Eriksson and Hallgren, 1975)

13
 Magnetite is the dominant ore mineral with minor amounts of hematite in both massive
and brecciated magnetite mineralization. The deposit is highly amphibolized and biotized.
Hornblende and biotite are also common with up to 20%. Ti and Mn in the system can be
carried by ilmenite (FeTiO3 ). It has been observed that the V is inside the magnetite lattice
because very few other minerals contain vanadium. A source of P could be apatite, but it was
rarely observed in the collected samples. The waste rock can carry up to 25% Fe in silicates.
The mine is operated by LKAB, resources and reserves are publicly reported on a yearly basis,
see Table 2.2

Table 2.2: Resource and reserves at Leveäniemi (LKAB, 2018)

Reserves Quantity 2018 [kt] Fe [wt%]

Proven 87000 48.1

Probable 9000 37.5
Resources

Measured 33000 46.6

Indicated 38000 42.5
Inferred 78000 41.8
Concentration

Milled ore from Leveäniemi 7000 44

Concentrate at Malmberget* 94000 67
*includes feed from other locations

2.3 Sampling

Two main types of samples have been collected in several campaigns, deposit samples and
process samples. Deposit samples have been used for mineralogical analyses, grindability
and concentration tests. Process samples have been used for mineral liberation analyses,
mass balancing and calibrating the process simulations. The difference between metallurgi-
cal and geometallurgical samples lies in their masses, i.e. metallurgical samples are larger.

14
 Table 2.3: Samples collected in the study

Metallurgical Geometallurgical Process samples

samples samples

Objectives mineral particles behavior liberation data

characterization information
crushing mass balancing
characterization
calibration of
simulations

Typical mass 60 kg/sample 10 kg/sample 2 kg/sample

Mass of 12 kg/sample 2 kg/sample 500 g/sample

duplicates

Procedure drill core drill core stable process

open pit samples open-pit samples 34 streams
one sample/45 min

Another difference between geometallurgical and metallurgical samples is that geomet-

allurgical ore behavior in the process is linked to the geological properties. Metallurgical
samples are usually more general and include bigger volumes, than geometallurgical ore
types.

Sampling from Aitik

Both geological and processing ore types have been already predefined at Aitik. Geometallur-
gical ore types have been established within this work based on the samples collected during
the campaigns and laboratory tests, see Table 2.3. The main hypothesis is that particles from
different geometallurgical ore types will behave differently in the process. This approach,
based on particles characterized by size and mineralogy, links the available geological infor-
mation and ore behavior in the process plant. Samples collected in Aitik were selected based
on several geological parameters.
Only laboratory tests define which of these geological parameters play a significant role in
processing and how geological variability can be linked to the processing variability. In order

15
to satisfy requirements given in Table 2.3 and provide good representation, the following
geological parameters were considered in sampling:

• Rock type: Hornblende banded gneiss, muscovite schist, biotite schist, biotite gneiss,
pegmatite,
amphibole-biotite gneiss, diorite.

• Mineralization: chalcopyrite, pyrite, pyrrhotite, magnetite, malachite.

• Alteration: k-feldspar, epidote, silicification, muscovite, albite, bleaching, biotite.

• Structure: foliation, veining, weathering, cavities

Drill core logging has been done in collaboration with mine geologists. Among all the prop-
erties, mineralization, alteration and structure are given on a scale between 0 and 5, where 0
represents no occurrence and 5 is the maximum possible occurrence. Additional comments
were also recorded whenever a specific feature was present. In total 22 samples were collected
from the drill cores and 10 samples from the pits. The total weight of the pit samples is 900
kg.

Table 2.4: Drill cores sampled in Aitik (the location is given in the local coordinate system)

Hole ID North East Relative length [m] End of the hole [m]

AI1540 7361.59 3498.8 224.86 555.3

AI1520 7168.99 4425.04 210.01 220
AI1529 7867.82 5400.09 85.96 269
AI1524 7860.49 5008.06 54.62 162.7
AI1495 7860.49 5008.06 54.62 353
Pit samples 7620 5696 0 -314
North pit samples 7370.253 4682.803 0 -286
Salmijärvi pit samples 8132.313 7082.839 0 -60

The sampling campaign was supervised by Aitik’s geologists in order to cover the most
frequently encountered lithologies. This was later verified by comparison to the literature and
discussion with Prof. Christina Wanhainen (personnal communication). The classification
into ore types was based on both the macroscopic texture and mineralogy.

16
Sampling in Leveäniemi

In the case of Leveäniemi, a similar distinction was made between metallurgical and ge-
ometallurgical samples. No previous geometallurgical ore types classification had been done
in Leveäniemi before, but a similar deposit (Malmberget) had been classified (Lund, 2013).
Geometallurgical ore types have been established within this work based on the samples col-
lected during the campaigns and laboratory tests, see Table 2.3. Process samples have been
collected as well for mass balancing and process simulation purposes.

Table 2.5: Drill cores sampled in Leveäniemi (the location is given in the local coordinate
system)

Hole ID North East Relative length [m] End of the hole [m]

LEV14022 7505512.549 182860.045 294.438 250.8

LEV14009 7505602.664 182977.017 305.799 350.1
LEV14072 7505406.09 182947.821 259.067 200.3
LEV14057B 7505306.901 182960.797 259.975 330.0
LEV14062B 7505406.629 182941.571 259.606 354.0
LEV14011 7505612.687 182875.954 295.980 251.5
SE Pit - - - -

Geologists from Leveäniemi ensured that the most frequently encountered lithologies
were collected. The samples were compared to documented cases available in literature or
in internal documents and classified according to their mineralogy and lithology.

17
Chapter 3

Methods, concepts and tools

3.1 Drill core scanning

In order to acquire textural information, the samples were first scanned with a high resolution
RGB camera. Images of the original data set are taken with a Canon Pro 815 camera and a
pixel size of 0.12 mm2, with a stable light source and images are stored as RAW files then
converted to TIFF files. For each ore type between 5 and 10m of drill cores from different
parts of Aitik’s ore body are scanned to obtain representative color images of each type. The
scale is kept constant. The data set is then randomly cropped and rotated to generate im-
ages of different sizes and 3 color channels for feature extraction and classification purposes.
The complete dataset includes 1509 images in 15 textural classes. The samples have been
collected from a selection of drill cores provided by the mine geologists. To ensure that the
dataset is representative of the textural variations in the ore body, the selected samples were
further studied and compared to previous works on Aitik deposit (Wanhainen et al., 2006,
2012). For image processing purposes, the Fiji image processing package is used (Schindelin
et al., 2012). WEKA (Hall et al., 2009), Python and MATLAB (Mathworks, 2016) are used for
the machine learning parts.

3.2 Geometallurgical tests

Tests used within this work were oriented towards the behavior of particles in grinding and
concentration. The geometallurgical comminution test (GCT) was developed at LTU and
aimed at providing an estimate of the Bond Work Index of the ore. It has been already suc-

18
cessfully applied to many deposits and was chosen, as it is one of the fastest ways to obtain
data for a large mass of samples (Mwanga et al., 2016). The parameters of the test procedure
are summarized in Table 3.1.

Table 3.1: Experimental procedure for GCT tests

Mass 220 g
Initial P80 3.35 mm
Equipment GCT mill
Conditions Dry
Step 1 1 min
Step 2 2 min
Step 3 5 min
Step 4 8 min
Step 5 11 min

The data collected consist of particle size distributions and the energy consumption after
each step. This allows the estimation of the Bond work index.
In addition to the GCT tests, Davis Tube (DT) tests (1921) and wet low-intensity magnetic
separation tests (WLIMS) were performed and taken into account for the construction of
geometallurgical classes. DT results indicate the amount of impurities in the magnetite con-
centrate and are a proxy for the maximum grade that can be obtained by magnetic separation
in the plant. WLIMS results show how the recovery of magnetite can vary when different ore
types are used as feed.

3.3 Methods for texture quantification and classification

The methods developed within this thesis are mainly based on two different axes.

Image processing

The first one are image processing methods from which a numerical model of a physical real-
ity is obtained. This step is crucial since it determines the choice of tools used afterwards. A
particular attention was given to the question of stationarity and ergodicity which originally

19
have been developed in geostatistics but are relevant in any spatial statistics problem (Chiles
and Delfiner, 2009, Lantuéjoul, 2013, Matheron, 1970). Additionally, signal processing pro-
vided tools to understand the separation of scales in mineral textures. The most basic method
is the Fourier transform which allows, for instance, to study the frequency of patterns in a rock
sample. The quantitative study of shape is also provided by decomposing the boundaries
in frequency components. Moreover, Fourier methods pave the way for practitioners to use
spectral methods in general like wavelets. In turn, wavelets lead to a general understanding
of scales and the need for an appropriate basis for texture analysis. While such a basis might
exist for each separate texture field, it is generally difficult to arbitrary choose one that is valid
for an entire dataset.

Machine learning

The previous paragraph leads to the second axis, machine learning methods. At the cross-
point between geostatistics, image processing, spectral methods and machine learning, the
convolution operator is placed. Instead of choosing an appropriate wavelet basis, one can
choose to learn it from the data itself. By minimizing the error between examples provided
by a mineralogist or a geologist and their predictions, these methods use the data to detect
regularities. Depending on the architecture of such models, one can choose to either perform
a segmentation per mineral phases, reduce textures to a vector of numbers or even combine
different data sources to produce an estimation of the ore performance in the processing
plant (in terms of grade and recovery for example). While the approach taken in geostatistics
or stereology is deductive by nature, as a set of hypotheses are held true and laws are derived
from them, the typical approach in machine learning is inductive or data-based. In the con-
clusion chapter some of these aspects and their implications are briefly discussed.
Among the most popular textural features extraction methods, two main approaches are se-
lected: the gray-level co-occurrence matrix (GLCM) statistics and local binary patterns (LBP).
GLCM descriptors have been introduced by Haralick et al. (1973) and recently applied to drill
core characterization by Becker et al. (2016). The frequency at which a given pixel value is in
a neighborhood of another (specified by a distance or offset and a direction) is represented
as a frequency table called co-occurrence matrix. It can be made rotation-invariant by aver-
aging the co-occurrence values in several directions (usually at an angle of 0, 45, 90, and 135
degrees).

20
Another class of descriptors are wavelets. As an extension of the frequency analysis per-
formed by Fourier filter banks, wavelet filters decompose the signal of an image into different
frequencies at different scales. The implementation used in this study uses a Haar wavelet
multi-scale decomposition (Antonini et al., 1992, Mallat, 1989a,b) and statistics derived from
the GLCM at each scale are combined. The advantage of the method is to separate different
scales present in the initial texture. The resulting feature vector contains the variance, correla-
tion, energy and homogeneity statistics of the GLCM at each scale of a wavelet decomposition
using Haar filters. Physically, it represents the association of minerals or phases with each
other at different scales. These features are chosen to consider scale properties of textures.
LBP have been introduced by Ojala et al. (2002) and have found multiple applications in
image analysis such as in face description (Ahonen et al., 2006) and specifically in texture
analysis either as such (Li et al., 2015), or in a modified version (Hegenbart and Uhl, 2015,
Li et al., 2015). The descriptors are based on gray-level images and a discrete circular neigh-
borhood. In the case of color images, the LBP can be computed on the image converted to
gray-scale values or on each color channel separately. Physically, it represents the local geo-
metric patterns observed on textures. These features are chosen to obtain local information
that is not sensitive to rotation.

Once the acquired images have been quantified, either by using wavelets, Haralick descrip-
tors, gray level or co-occurrence matrices, these textural indicators can be used to classify
ore textures into different types. Formally, a classification problem is to assign a label to
an image. In the case of geometallurgy, the problem becomes either to assign a lithological
class to an image, or a visual description to an image (for example the number and size of
veins) but also to combine a variety of input data about the ore (images of the ore at different
scales, mineralogy, results from grinding tests, recovery in a laboratory scale test) about the
ore to create groups of samples that behave similarly in the comminution and concentration
process. These classes are referred to as geometallurgical classes.

Two approaches exist for automated texture classification: (i) the supervised, which relies
on two subsets of the initial data: the training set and the test set and (ii) the unsupervised,
which only uses unlabeled data as an input (images not associated with a class). When a
supervised method is used, the training set consists of labeled data (images associated with a
class) from which the model learns a classification function. A test set consists of unlabeled

21
data to which the classifier is applied. Usually, the model learns to discriminate by minimiz-
ing an error. As such, supervised algorithms applied to images are equivalent to a dimension
reduction. Supervised classifiers are useful when knowledge about the classification already
exists. This is the case when classifying mineral textures from drill core images. Trained geol-
ogists are experts in classification as they are able to extract most of the relevant features in a
limited time. Additional information can be taken into account in the features, for example,
grinding tests results (BWi), modal mineralogy or chemical assays. In that case, the classi-
fication does not depend on visible information only, making it suitable for geometallurgy.
Classifiers have different performances and must be validated. A common way is to split the
initial data set between a training and a testing set. If the classifier is built on the training
set then applied to the testing set, one can compare the predicted class to the actual one. A
confusion matrix can be built, which has the predicted classes as columns and the actual
classes as rows, the values are the number of instances predicted to be of each class.

The main idea behind classification algorithms is to find a transformation that reduces
the dimension and defines boundaries between images of different classes. Classifiers are
functions that assign a texture class to each feature vector. Two main families are used,
unsupervised and supervised classifiers. Unsupervised methods rely on the data itself to
learn a classification whereas supervised methods learn a classification based on examples
provided by the user. Supervised classification can be thought of as teaching by examples: a
sufficient number of example pairs (feature vector, actual texture class) is provided and the
classifier builds a function that assigns a class number to the feature vector. In general, the
obtained models map an input to an output but cannot always be described by a linear or
analytic function. This approach leads to a classifier model that is trained to reproduce and
generalize the behavior of a human operator that provided the examples as a training set.

Unsupervised methods

In practice, unsupervised methods do not guarantee a classification that respects a preexist-

ing descriptive classification for example done by geologists. K-Nearest Neighbors methods
(Aggarwal, 2015, Bandyopadhyay and Pal, 1991) are unsupervised and rely on the calcula-
tion of a distance between feature vectors. Starting with the center of each class randomly
distributed in the feature space, a feature vector is assigned to the class that has the closest

22
center. The key of the algorithm is to measure distances between vectors in relatively low
dimensions compared to distances between images, which is a much more complex task.
After one step, the location of the center of each class is updated. By iteration, each feature
vector will be assigned to a class. This method can be thought of as grouping textural images
that have similar numerical features.

Supervised methods

Random Forests (Breiman, 1996), Support vector machine (Pérez-Ortiz et al., 2016), artificial
neural networks (Hebb, 1949, Rosenblatt, 1961) and their convolutional counterparts (LeCun
et al., 1995) are part of the supervised methods. These techniques are useful when examples
of correct classification already exist. In that case, the models can learn a classification func-
tion from the data. Random forests are based on the construction of decision trees. Each
decision tree is built by selecting, at each node, a criterion that maximizes the information
gain regarding the textural class. By doing so, each tree offers a sequence of binary decisions
that lead to a full classification. In random forests methods, the class of a sample is chosen
based on a vote of many decision trees.
A support vector machine (SVM) is based on finding a hyperplane that maximizes the dis-
tance between feature vectors belonging to different classes. Since the boundary between
classes can be complex in high dimensions, a function that maps the feature vector into a
space in which the boundary is linear is computed by optimization (Boser et al., 1992). Once
the features have been transformed, planes that separate each class are built. Each plane is
described by its normal vector (support vector) and projecting each feature on this vector
separates the data into classes.
Artificial neural networks (ANN) construct a function which is a series of weighted sums. The
weights are organized in layers that define the order in which the sums are computed. The
input layer has the same number of weights as the components of the feature vector and the
last layer has one weight for each class. Several layers can be placed in between and connect
the input to the output, they are called hidden layers. At first, the weights can be random
numbers. When an example data is provided, the difference between the actual class to which
the image belongs and the predicted class by the network provides an error. As the number
of examples increases, the network adapts the weights to minimize the classification error.

23
Conclusions

No single method can guarantee its performance either for extracting feature vectors from
an image, nor for classifying textures. As a result, different combinations of features and
classifiers are tested. Since the methods are applied to drill core scanning, two parameters are
taken into account. Accuracy is the main objective in order to provide a stable classification
that makes sense from a geological and mineralogical perspective. The second important
parameter is the time needed to extract features and classify the image. In a drill core scanning
context, the processing time should be reasonably short to maximize the amount of material
that can be scanned per day.

24
3.4 Sequential decisions

Geometallurgical models are often defined as a set of tools to improve knowledge of an ore
body by measuring relevant information for the concentration process (Lund and Lamberg,
2014). The objective is to gain information but also to use this information to improve the
grade or recovery of the mineral of interest, optimize the throughput or reduce the environ-
mental footprint of the process. In a geometallurgical program, a large amount of data from
different sources is acquired, structured and made available as a basis for decision-making in
mining and process planning. Sequential decisions have to be made in order to maximize an
objective function, usually aggregating economic objectives, production and environmental
targets, or to minimize a regret function. Besides the selection of different processing options,
these decisions often involve dealing with uncertainties with respect to the variations of the
ore feed and other input variables, model parameters and the model structure. A special
challenge arises from cases where certain information is not available yet, i.e. the decision
must be taken a priori to knowing the outcome of the selected processing option. Numerical
methods based on bandit models have been proposed for this type of problems but, to the
authors’ knowledge, have not been made available yet to geometallurgy.

In a situation where several ore bodies are mined and several processing options exist, as
shown on Figure 3.1, the steps between the mines and the beneficiation process (i.e. steps as
ore handling, transportation and blending) are neglected as a first approximation. A further
simplification is to consider only the grade of one mineral of interest and that every ore body
has the same grade distribution resulting in the same feed properties to all the processes.
Additionally, all the processes are considered stationary with respect to time. The actual
performance of each process is unknown but represented as a statistical distribution. The
task to be solved is then to maximize an economic objective function through assigning an
ore block at a time (sequentially) to the processing option that has the highest expected payoff
while not knowing the performance of the other processes. At each discrete time step, the
plant operator (or software agent) can choose one of the processes in which the ore block
will be processed. Once the block has been processed, a performance indicator relative to
the chosen process (e.g. grade and recovery of the concentrate, presence or absence of a
deleterious element) is observed by the plant manager (also called agent since it can act

25
on the system). With time, the available information about each process increases. The
model describing this procedure is commonly called “multi-armed bandits” (MAB) setting,
by analogy to a player in a casino who can choose between several slot machines equipped
with levers called arms.

Figure 3.1: General approach to sequential decisions in geometallurgy

The multi-armed bandit is a simple model which exemplifies the exploitation-exploration

trade-off in reinforcement learning. Solutions to this problem have numerous practical ap-
plications from sequential clinical trials to web-page ad placement (Bubeck et al., 2012). This
work uses the stochastic setting in which an agent is given access to a collection of unknown
reward distributions (arms). The agent sequentially selects a reward distribution to sample
from, to maximize its cumulative reward. One of the most widely used strategies for stochas-
tic multi-armed bandits is the Upper Confidence Bound (UCB) algorithm, which is based
on the principle of optimism in the face of uncertainty (Agrawal, 1995, Auer, 2002, Lai and
Robbins, 1985). Garivier and Cappé’s KL-UCB algorithm utilizes tighter upper confidence
bounds to provide a policy with sharper regret bounds and a superior empirical performance
(Garivier and Cappé, 2011). One can adapt the definitions to geometallurgy. Regarding ore
variability or the spatial side of geometallurgy, different cases occur :

• One ore type (perfect blending strategy)

26
 • One ore type with variability

• Several ore types

On the process side of the model, several cases exist. The main distinction between processes
is their time invariance or not. A time-invariant process will produce the same output given
the same input regardless of the time. In this work it is called time-stationary or stationary
(S).

• One stationary process (S-process)

• One non-stationary process (NS-process)

• Several S-processes

• Several NS-processes

If a bandit model is chosen, a natural aim is to develop a strategy (or policy) to maximize
the overall performance in the long run. An alternative aim would be to identify the best
process with high confidence. In both cases, the concepts of reinforcement learning are ap-
plied, i.e. the acquisition of new knowledge from the process (called exploration) and the
subsequent improvement of decisions based on the existing knowledge (called exploitation).
At the beginning, decisions must be taken without knowledge about any of the processes.
The software agent must explore the different processes and, once the number of observa-
tions increases, exploit the correct process. Objectives used within process optimization
and decision-making can be either of binary or continuous type, thus requiring different
mathematical formulations for the description of the observations, see Table 3.2.

Table 3.2: Different rewards for geometallurgy

Objec- Distribu- Values Examples

tive tions

Binary Bernoulli Discrete ∈ 0, 1 Quality requirements of the concentrate,

environmental criteria
Real Exponen- Continuous Performance score : grade, recovery, throughput
tial ∈ [0, 1]

27
 To describe the problem, the usual terms and notation from the bandits’ literature are
used:

• At each time step, a decision is taken about which block will be assigned to which
process

• The response of each process to each block is random to some extent. In particular, the
probability distributions are unknown

• The reward can only be observed after a decision was made

The best arm is noted a ∗ and its mean reward is µa ∗ . The overall performance of a policy
can be measured by the cumulative regret RT which is the sum of the differences between the
rewards obtained by the agent and the rewards that the agent could have obtained by always
selecting the arm yielding the highest reward. The challenge is to build “good” policies for
a specific class of rewards. This has been approached by Lai and Robbins (1985) for a single
parameter exponential family and extended in the works of Burnetas and Katehakis (1997).
In short, to achieve a cumulative regret of the order of t α with α ∈ [0, 1], any sub-optimal
l n(t )
arm that has a mean µa must be selected on average more often than K L(µa ,µa ∗
), where KL
is the Kullback-Leibler divergence that quantifies the difficulty of making a choice between
two arms. These policies are called asymptotically optimal which means that this property is
desirable in practice since it provides bounds for the cumulative regret when the number of
time steps is large.

One of the early examples of designing an optimal strategy (or policy) in a sequential
decision process can be found in the works of Thompson (1933). The problem addressed
is to determine whether a given medical treatment is better than another. The success or
failure of a treatment is modelled as a Bernoulli trial. The solution adopted is to model a
priori probabilities, draw a sample from both distributions, choose the treatment that has
the highest probability of success and update the a priori distributions accordingly. This
approach is Bayesian, i.e. unknowns are described by probability laws while the observations
are non-random – as opposed to frequentists methods. Thompson’s sampling is asymptoti-
cally optimal (Kaufmann et al., 2012).

28
 When different processing paths exist, two extreme strategies can be adopted: either ran-
domly explore all the possible choices (pure exploration) or choose one process and always
select it (pure exploitation). In the first case, one can see that the optimal process is chosen
on average as much as any sub-optimal process. In the second case, either the best process is
selected, and, in that case, it is a successful strategy, or another process is chosen and the re-
ward will always be sub-optimal. Additionally, the probability to choose the optimal process
from the beginning decreases with the number of available choices. This is referred to as the
exploration-exploitation trade-off. One can see that a reasonable strategy should combine
exploration and exploitation.
A simple yet effective approach to the exploration-exploitation trade-off is the epsilon-greedy
1−²
strategy in which the arm with the highest mean is selected with probability K and the oth-
²
ers with a probability K where K is the number of arms. Epsilon is often chosen around 0.1
which means that the probability of choosing the arm with the highest mean reward is 0.9
and 0.1 for the other arms. Epsilon-greedy strategies solve the trade-off by assigning a proba-
bility for exploration and exploitation. Another family of policies uses confidence bounds on
estimates of the mean of each arm. As an example, the UCB algorithm chooses the arm that
has the highest upper confidence bound (hence UCB). The construction of the interval itself
is detailed in Auer and Fischer (2002) but in essence, the upper bound consists of the sum of
q
2l og (t )
an exploitation term µa (t ) (empirical estimate of the mean) and an exploration term N a (t )

where Na is the number of times the arm “a” was selected. This approach is frequentist since
it is based on an estimation of the mean that improves when the arm is selected. A variant
based on UCB and the Kullback-Leibler divergence called KL-UCB offers better guarantees
in terms of regret (Garivier and Cappé, 2011).
Applications of the algorithms explained above to mineral processing and geometallurgy are
possible in several cases presented in Table 3.3. A distinction is made between time stationary
(S) and non-stationary (NS) processes and ore classes. The former represents a process that
has a steady performance in time. This assumption holds when a short time is considered but
when longer times are studied, the process is often non-stationary. With regard to ores, the
distinction can be made between a constant ore class (one class of constant physical prop-
erties), a steady or targeted ore class (where the feed to the process has the same physical
properties on average) and k S-classes where the ore is split into different geometallurgical
classes representing various process performances, compare the examples given in Table 3.2.

29
 Table 3.3: Different cases for decision-making

1 S-process 1 NS-Process k S-Processes k NS-Processes

1 ore steady dynamic steady+scenarios dy-

namic+scenarios
1 variable steady+scenarios dy- case 1 case 2
ore namic+scenarios
N ore types steady+scenarios dynamic + case 3 case 4
scenarios

As presented in Table 3.3, one can simplify these cases and evaluate whether commercially
available software does offer a solution. Conceptually, a difference can be made about sim-
plifications of the ore type and the process classification. This is supported by previous work
in geometallurgy (Lishchuk et al., 2015b). Classification into ore types has been done using
dimension reduction techniques such as principal component analysis (PCA) (Andrusiewicz
et al., 2014, Keeney, 2010). Linear dimension reduction offers the advantage that distances
have an intuitive meaning in terms of physical properties such as mineralogy. Non-linear
mappings can also be used and theoretical guarantees exist regarding the existence of such
a mapping and the preservation of distances (Kłopotek, 2017). A stationary process has the
same statistical distribution with time. In mineral processing plants, this represents stable
operating conditions, i.e. a plant that does not display a drastically different average behavior
in time. An example of a stationary signal is a constant measurement affected by statistical
noise. An example of a non-stationary process is a plant affected by seasonal factors like
temperature (Lin, 1989). Several of these cases can be solved using existing numerical tools,
namely mineral processes simulators. Current software allows the simulation of stationary or
non-stationary (also called dynamic) processes. Additionally, different scenarios can simu-
late changes in the feed composition and process parameters (Outotec, 2019). The question
of optimal allocation is not directly addressed in these simulations. Moreover, the construc-
tion and calibration of such models rely on sampling and mass balancing, and simulations
are based on unit models that often include empirical terms (King, 2001). While empirical
parameters allow calibration, they depend on measurements. Their values can be regarded

30
as estimates of a real probability distribution. This hypothesis points to the use of stochastic
settings to model mineral beneficiation processes.
The cases where multiple processes exist are more complex (see for example cases 1 to 4 in
Table 3.3). From now on, no difference is made between simulations of processes (used for
example during the design phase of a process) or an actual plant (in that case, the strategies
can either complement the design of experiments or act as a guide to support decisions).
If an interface exists between the plant data (or the simulation software) and the software
applying the strategies, the methods introduced are useful in both cases. When several pro-
cessing routes exist, the number of parameters in each plant increases, which leads to either
more unknowns (in the case of actual plants) or more calibration data required (in the case of
simulations). If the plant parameters are controlled by an automation system or if statistical
variability is introduced in the plant model, the problem becomes intractable. To overcome
this difficulty when decisions are made, bandit algorithms are useful.
Since each ore type has a direct impact on the process performance, the available information
should be used by the policies. Bandit methods where side information is used are referred
to as contextual bandits and an extensive literature is available on the subject (Bubeck et al.,
2012, Reeve et al., 2018). Geometallurgical ore types are supposed to behave similarly in
the process since they were classified based on metallurgical tests and texture. As such, two
samples with close geometallurgical indicators should have the same optimal process. A
logical choice would be to look at the k-nearest neighbors (kNN) and use a similar strategy.
A relevant method referred to as kNN-UCB, based on nearest neighbors and UCB, offers a
straightforward implementation (Reeve et al., 2018). Contextual data in geometallurgy can be
a normalized vector that contains mineral grades, grindability indicators, textural properties
or a scaled version of a combination of all these properties. Geometallurgical ore types are
usually a classification of high-dimensional data.
The kNN-UCB algorithm proceeds in several steps:

1. For each arm (process), the number of neighbors k is chosen to minimize a measure of
uncertainty

2. Given a context x at time t, it selects the k neighbors around x

3. A local upper confidence bound on the reward functions is computed

4. The arm with the highest bound is selected

31
Several hypotheses are needed to establish optimality bounds, i.e. hypotheses that suggest
smooth and bounded reward functions. For the problem stated here, a simplified set of
assumptions involves:

• Dimension assumption: The context data is supported by a d-dimensional manifold

where d is smaller than the dimension of the context. This assumption means that
the vectors that contain the mineralogy and grinding indicators can be projected on a
surface of lower dimension than the vectors themselves.

• Margin assumption: The arms have rewards which are different enough.

• Lipschitz assumption: For similar context vectors the rewards are similar. Another way
to interpret the assumption is that the reward functions should be varying smoothly.

• Rewards are bounded in [0, 1].

Conclusions

In order to test the different methods, reward functions that respect the respective hypotheses
of the algorithms are constructed. Once the reward functions exist, three algorithms are
tested : random decisions as the worse case scenario, an epsilon-greedy strategy that mimics
a reasonable human behavior and a UCB policy. Depending on the reward function, its values
and its variables, different adjustments are done to the simple UCB method. The first variant
is the sliding-window that let the algorithm forget after a specific period. As a result, in non-
stationary cases, the mean reward is only be influenced by a limited variability. The second
variant is to limit the calculation of the mean reward to a local neighborhood by a kNN step.
This represents the fact that when different ore types are present, similar ores should behave
similarly.

32
Chapter 4

Numerical results and discussion

4.1 Geometallurgical classification

Geometallurgical classes were defined using the lithology, the chemical assays, the modal
mineralogy and the grindability tests. The developed classification has been further tested
in several cases, for Aitik see Tiu (2017) and Bilal (2017), and for Leveäniemi see for example
Singh (2017) and Lishchuk et al. (2018b).

4.1.1 Mineralogical characterization

In both cases, epoxy resin mounts were prepared and analyzed by Auto-SEM with EDS for min-
eral classification. The complete mineralogy was obtained for both deposits, Leveäniemi and
Aitik. The data was in line with both previous works and literature data for similar deposits
(Lund, 2013, Wanhainen and Martinsson, 1999). Figure 4.1 shows that the mineralization
in Leveäniemi is dominated by iron oxides, the distribution between hematite, magnetite
and maghemite could not be obtained from the SEM, but by using chemical assays including
Satmagan values. Gangue minerals are mostly silicates with an important part of micas and
plagioclase. These silicates have been shown to undergo entrainment and entrapment during
the WLIMS concentration in the plant.

33
 100%

90%

80%

70%

60%

50%

40%

30%

20%

10%

0%
L1 L6 L9 L12 L13 L2 L3 L4 L5 L7 L8 L10 L11

Ap Cal Amph/Px Mica Na-Pl Fsp Qtz Fe-oxides Ti-Mgt Sulph Other Silic Others

Figure 4.1: Modal mineralogy of the Leveäniemi samples

In the case of Aitik, the chalcopyrite grade is low but differences in the gangue mineralogy
are observed. The association of biotite and plagioclase was already noticed in the drill cores.
The presence of magnetite in the samples as a result of the IOCG overprint is noted. Pyrite
and Pyrrhotite are important in this case since pyrite must be separated from chalcopyrite
during froth flotation.

34
 100%

90%

80%

70%

60%

50%

40%

30%

20%

10%

0%
C1 C2 C3A C3B C4 C5 C6 C7A C7B C8 C9A C9B C10A C10B C11

Ccp Py Po Mgt Hem Qtz K-Fsp Pl Bt Hbl Chl

Ms Epi Tit Gar Cal Tur Ap Zrn Ilm Scp Others

Figure 4.2: Modal mineralogy of the Aitik samples

Obtaining the modal mineralogy required a preliminary step using XRD. The main min-
erals were first identified to guide the SEM runs.

4.1.2 Grindability tests

Each sample was tested using the GCT procedure. The tests were repeated three times and
at each step, the complete particle size distribution (PSD) and energy consumption was ob-
tained. The BWi estimates for Leveäniemi and Aitik are displayed on Figures 4.3 and 4.4
respectively.

35
 Figure 4.3: GCT tests results for the Leveäniemi samples

In the case of Leveäniemi, the estimated BWi values are quite similar. Low values of this
order of magnitude is reported in literature (Mwanga, 2016).

36
 Figure 4.4: GCT tests results for the Aitik samples

Aitik samples display more variability with a clear outlier (C1) which is the pegmatite
class.

37
4.2 Texture quantification and classification

After geometallurgical classes are obtained from the tests and the lithology, textural infor-
mation can be added. The main object of interest in measuring textures is the mineral map
(Koch et al., 2019c). It is a model of the assemblage of minerals and pores in space. Since it is
very generic, numerical methods from other disciplines such as spatial statistics and image
processing are used.
The first step is to obtain a mineral map which is a spatial field in which each phase (includ-
ing the background and the voids) is coded by a number. One of the tools developed in this
work aims at acquiring such information in a time and cost effective way (Koch et al., 2019a).
Visible light imaging combined with training data proved to perform well when the model
is deposit-specific. The training set consists of groups of pixels identified by a geologist as a
same mineral phase. The group of pixels is then given the label corresponding to the mineral.
By repeating the process for all the targeted minerals a dataset is built. This gives the oppor-
tunity to supervised methods to learn from human experience. The results after training are
displayed as a color overlay on the initial image, see 4.5.

Figure 4.5: Example of mineral phases segmentation for a drill core from Aitik

Based only on visual inspection, any quantitative validation of the method is difficult. But
by comparing the results with the ones obtained by a SEM-EDS, an evaluation is possible. A
small region of the drill core is selected, drilled out of the core and polished for SEM scanning.
Since the resolution is much higher on the microscope, details are lost but the overall geome-
try is respected and most phases are correctly identified. On Figure 4.6 an example of result
is presented. Visual inspection is sufficient to identify a major flaw in the method but not to
quantify the approximation error.

38
Figure 4.6: Qualitative comparison between mineral identification methods. Top center: SEM
results, Top right: visible light results, Bottom: initial sample

Once both mineral maps have been obtained, scaled and aligned, statistics are computed.
The simplest way is to compare the proportion of each phase one by one. Because each phase
has a label, the calculation is limited to counting the number of pixels that have the same
values. Once all the phases have been counted, a bar chart is plotted as seen on Figure 4.7.

39
 Comparison between SEM-EDS and visible light image
estimated mineralogy
0.3
0.25
Wt %/100

0.2
0.15
0.1
SEM-EDS
0.05
Visible
0

Main minerals

Figure 4.7: Quantitative estimation of the error

Once a mineral map is obtained, three main variables are relevant

• the shape

• the size of the particles or phases

• the mineral associations

Shape can be effectively described by using Fourier descriptors, as an example see 4.8. This
is achieved by considering successive distances from the center of the particle as values of a
function. A Fourier transform of this function is performed and the coefficients are calculated.
By selecting coefficients from low frequencies to high, scales of the particle’s boundary can be
isolated. By taking the inverse transform of these coefficients the shape can be reconstructed.

40
 Original particle Fourier reconstruction as a Number of
function of the number of Fourier
descriptors descriptors
1

10

50

200

Figure 4.8: Fourier shape descriptors of one particle

One of the key aspects of the Fourier transform is that it can both measure a shape and,
from the coefficients, reconstruct the shape. From an information theory point of view this
represents a compression of shape information. Ideally, quantified textural information
should allow a perfect reconstruction of the initial texture. In fact, the Fourier transform
is linked to another operator called convolution (Katznelson, 2004). Convolution filters are
tensors that can detect patterns similar to themselves in images. Instead of using prede-
fined filters to analyze images, they can be learned in convolutional neural networks during
the training phase. Figure 4.9 displays some filters from a trained network. Several textu-
ral features likes the presence of veins and their orientation or disseminated phases can be
observed.

41
Figure 4.9: Convolution filters learned during texture classification applied to one of the
image patches

Classical methods to measure particles from image analysis rely on one or two axes of an
equivalent ellipsoid. Size information can also be obtained by morphological filtering. By
using a sequence of morphological closures and openings, one can construct a morphological
sieve that can produce particle (or grain) size distributions. This presents the advantage of
not relying on any assumption about the shapes present (Serra and Soille, 1994).
Association of minerals is usually done by computing the difference in pixel values at two
locations. Results obtained for the texture quantification in most cases are numerical, as
such, representing them is not trivial. Nevertheless, one can try to project them into a lower
dimension space by applying a principal component analysis (PCA) for example, a Fourier
transform or a non-linear stochastic method for example t-SNE (van der Maaten and Hinton,
2008). It should be noted that not every projection algorithm provides the adequate structure

42
to define a distance between projections. Another way to look at it is from the vectors formed
by the extracted values. If placed in a space equipped with a norm, the appropriate basis, and
a product, vectors may be projected to yield lower dimensional data. The idea of projection is
at the core of this work since it allows to take all the acquired information, from the drill core
characterization but also from grindability and metallurgical tests, and use them in either
spatial models or process models. As an example of a projection, three selected frequencies
of local binary patterns values were used and plotted. The result given in Figure 4.10 indicates
that data is clustered according to their textural types.

Figure 4.10: Selected frequencies of LBP for different textures. Each color labels one texture

Once the features have been extracted, a classification is done by using machine learning
algorithms. The main results show that these classifiers (supervised or not) are able to suc-
cessfully classify Aitik textures to their class. All the data is tested with each classifier using
15-fold cross validation. For this purpose, the original sample is split into sub-samples. One
of these is kept as validation data while the rest is used to train the model. The procedure is
repeated 15 times so that all the data has been used for validation. The results are averaged
over the 15 folds to produce a single estimate. This is the list of selected data sets :

43
 • LBP: Local binary patterns

• HaarLevel2: GLCM statistics from Haar wavelet decomposition at level two

• HaarLevel3: GLCM statistics from Haar wavelet decomposition at level three

• LBP+HaarLevel2: Combination of LBP and HaarLevel2 descriptors

and classifiers:

44
 • RF: Random Forest

• SVM: Support vector machine with sequential minimal optimization (SMO)

• kNN-p2: k-Nearest Neighbors classifier with a Euclidean distance

• kNN-p1: k-Nearest Neighbors classifier with a Manhattan distance

• kNN-p1.5: k-Nearest Neighbors classifier with a Minkowski distance (p=1.5)

• ANN: Artificial neural network with three hidden layers

Different metrics are used to compare classifiers. In this case, the Matthews correlation
coefficient (Matthews, 1975) is presented in Table 4.1. It is a statistic derived from the matrix
that has predicted classes as rows and actual classes as columns and contains the count
of each sample that is classified (confusion matrix). It indicates the correlation between
the predicted class and the observed class. When combined with the time taken for the
acquisition of the textural descriptors, an overall score can be computed. Since it represents
the trade-off between accuracy and time needed, it must be low for a fast and accurate drill
core scan. Each set of pairs (features,classifier) was ranked accordingly and presented on
Table 4.1.

Table 4.1: Overall performance of different combinations of features and classifiers

Descriptors RF SVM kNN-p2 kNN-p1 kNN-p1.5 ANN

LBP 0.87 0.72 0.85 0.85 0.82 0.82

HaarLevel2 0.64 0.64 0.54 0.56 0.28 0.63

HaarLevel3 0.52 0.52 0.41 0.43 0.00 0.51

LBP+HaarLevel2 0.57 0.57 0.48 0.48 0.25 0.57

4.3 Data fusion

In the Leveäniemi case, since metallurgical tests have been done, the fusion of image data
from the drill cores and the iron grade could be used to estimate the performance in terms
of grade, recovery and throughput in the plant. A convolutional neural network, designed

45
specifically for data fusion has been developed (Koch et al., 2019b). Figure 4.11 shows the
architecture of such a model.

INPUTS

Iron assays
Drill core images
[wt%]

PRE-PROCESSING

Cropping Normal fitting

Frame
[128x128x3,1]

MODEL

Convolutional Neural
network network

Last layer of neurons

DT WLIMS
BWi
Iron Iron
estimate
grade recovery
[kWh/t]
[wt%] [wt%]

OUTPUTS

Figure 4.11: Structure of the fusion model

It takes as inputs (i) the RGB images and (ii) the iron grade assays. While the image branch
of the network performs convolutions to extract meaningful features at different scales, the

46
other branch uses densely connected layers of neurons to learn the iron grade distribution.
These two branches are then concatenated and connected to one dense layer. That layer is
then connected to three neurons, one that predicts the GCT Bond Index, another one that
predicts the WLIMS recovery and the last one that predicts the DT iron grade. As a result, a
fusion model that can use both drill core information and grade information is obtained.
The model was trained and predictions were obtained for each variable as displayed on Fig-
ures 4.12 , 4.13 and 4.14. Training the model requires the minimization of the error. This is
usually done via a stochastic gradient descent algorithm (Bottou, 2010) or one of its many
variants like Adam (Kingma and Ba, 2014).

Figure 4.12: Scatter plot of the estimates of the GCT Bond Work Index by the fusion model

The predictions in the case of GCT values indicates that the model seems to predict values
only between 6 and 9.5 at least for the validation set. One reason could be the lack of training
data with values outside this range or the presence of unusual points in the testing set.

47
 Figure 4.13: Scatter plot of the estimates of the WLIMS recovery by the fusion model.

Figure 4.14: Scatter plot of the estimates of the DT iron grade by the fusion model.

The vertical groups of points that can be seen on Figure 4.13 are the result of the measure-
ments being done per geometallurgical class. As a result, each sample of the same class has
values close to the average in the case of low standard deviations. The extremely low recovery

48
values come from breccia samples.
The obtained model is stored as one file that contains both the structure of the network and
its weights. Combined with a client-server interface, the model can be used in an industrial
context where the server hosts the model and the client receives the predictions. To train such
a model, two approaches are available: a local one where the server that is trained is located
on site or a remote server. The advantages of the first setting are to offer access to hardware
in case of failure, but the main weakness is the presence of an additional computer inside an
industrial environment, potentially in presence of dust and vibrations. In the second setting,
one limitation is the available data bandwidth for sending and receiving the batches of data
over through a network. The model should be trained regularly on a need basis or when the
textural ore types change dramatically. Regular incremental and full backups can be done
with any modern deep learning platform but a strategy on data and model maintenance is
advisable.

4.4 Sequential decisions

Given the amount of data collected and analyzed, it is possible to use in conjunction with
Bandits algorithms to make optimal decision when the performances vary as a result of
a different process or a different ore type or both (Koch and Rosenkranz, 2019). For each
case presented on Table 3.3, numerical experiments have been conducted to compare the
different strategies or policies. The cumulative regret is used as the main metric since it is
representative of the optimality of a given policy. The process models were created in HSC
Chemistry and calibrated with plant surveys and mass balances. Simplified models were
derived from the full models to reduce the computation time. The policies and plots were
implemented in MATLAB 2018b. When working with models for flow sheet selection or design
purposes, simulations can run in parallel. Each experiment consisted of 50000 blocks and was
run 150 times due to the stochastic nature of the policies, and intervals of confidence were
calculated from the repeats. These are not to be confused with upper confidence bounds that
are calculated inside the UCB algorithm during each round. Two process models are used,
one involving a magnetic separation flow sheet for processing iron oxide ores and a froth
flotation circuit for copper sulphide beneficiation. Context is provided by the experimental
data.

49
 06KV001 06SE011 06KV002 06SE021 06KV003 06SE031
ROM_1

KV1_5 SE1_CONC_6 KV2_9 SE2_CONC_10 KV3_13 SE3_CONC_14

SE2_TAIL_11 SE3_TAIL_1_15
SE1_TAIL_7

SE2_TAIL_B_19
06DO001 06DO002

SE3_TAIL_B_22

SE3_TAIL_1_16

Figure 4.15: Flowsheet of the iron process at the time of sampling

The first case is an application of bandits to Leveäniemi ore and the concentration process
presented in Figure 4.15. The simplest situation occurs when the properties of the feed ore
are stationary. In practice, this is achieved by using the mine scheduling software to stabilize
the feed grade for example according to the iron and vanadium grade (Vos, 2018). From the
process point of view, one can construct a concentrate quality indicator that aggregates the
presence of desired minerals and the absence of detrimental elements (vanadium or silicate
for example). The product specification for silica can be modelled as a Bernoulli distribution
giving 1 if the concentrate is within the requirements, 0 else. The agent can choose between
4 different processes modelled as

• R1 = B (0.6) : a circuit equivalent to the flowsheet on Figure 4.15 with a silica reverse
flotation step to remove any silicate from the concentrate

• R2 = B (0.4): this is the circuit presented on Figure 4.15

• R3 = B (0.3): a circuit equivalent to the flowsheet on Figure 4.15 with only two grinding
stage

• R4 = B (0.1): this circuit is equivalent to the flowsheet on Figure 4.15 with only one
grinding stage

where B µ is a Bernoulli distribution of mean µ. Since these distributions are unknown

¡ ¢

to the agent, different policies can be employed:

• A random selection of the process

50
 • An epsilon-greedy strategy

• KL-UCB policy

In this case, the benefits of having a good policy are clear. Random choices lead to the highest
regret, followed by the epsilon-greedy strategy. Policies based on confidence bounds perform
best. Figure 4.16 depicts the cumulative regret as a function of time for the selected strategies.
For the KL-UCB strategy even the interval of confidence over the 150 repeats is shown, thus
illustrating that the optimality is given with a high statistical reliability.

Figure 4.16: Cumulative regret for the iron ore case with stationary Bernoulli rewards

In the second case about Leveäniemi, each process presents temporal variations while
the feed varies around a controlled average. When the performances of a process are not
known in advance, especially when the variation cannot be predicted, computing the aver-
age will fail since the period of the variations is not known. While the period can be estimated,
the presence of noise can make the estimation difficult if the amplitude is small. Moreover,
any other change to the plant (maintenance or automatic control) could influence the esti-
mates. Using sliding windows reduces the difficulty but still requires the window’s size to
be calibrated. Figure 4.17 shows the regret over time for non-stationary processes. While a
random choice is still the worst strategy, the epsilon-greedy strategies seem to fail at adapting
to changing rewards using the same time window as UCB strategies.

51
 Figure 4.17: Cumulative regret for the iron ore case with non-stationary Bernoulli rewards

The third case and fourth cases are applications of bandits algorithms to Aitik. The con-
centration process is a series of froth flotation cells with intermediate grinding stages, see
Figure 4.18. The results from the geometallurgical characterization are used to provide con-
text as vectors. Once contextual data is introduced, the optimal decision depends on the
block.
FD4000
Ing malm 2

BL4101 FA4110-13 FA4114-18 FA4119-22

Ing malm 1 IngCond2
IngFlot2 MpRå2 Stream 2 Mp2

NSTailings

Högsvavliga
KV4210
S kon L2
RåConc2
RepScavB
S kon L1
RåConc1
Bl4001 FA4014-18

MpRå1
IngCond1 IngFlot1 Mp1

ConcSvavel

FA4010-13 FA4019-22
MpRep+S
RåCycUF ScavCon1

MpRepScavA
ScavCon2
RåCycOF
MpRepScavB MpRevScav
RåRepB RepScavA 2
PS3 MpRåRepB
Cl1Feed
RåKvUt

FA4510-14 FA4530-33 FA4550-52

RåRepKonB KoncRep2 KoncRep3 KoncRep4

RåRepKonA

RepScavA

MpRåRepA
MpRep4

RåRepA
MpRep3

MpRep2
MpRep

FD4560

RepScavKoncA

RepScavKoncB

Figure 4.18: Flowsheet of the Aitik process at the time of sampling

52
 As seen before, a purely random strategy fails here as well as presented in Figure 4.19,
but the difference between epsilon-greedy strategies and UCB is much larger. This result
could be explained by the fact that kNN-UCB takes the context into account by taking similar
decisions for close ore types. Epsilon-greedy strategies rely on the empirical estimate of
the mean only, which makes it difficult to distinguish between the impact of a different ore
type and variations in the process. In practice blending can offer a solution and reduce the
problem to previous cases in which the grade of the feed is targeted. However, it is not always
possible due to mining or equipment constraints.

Figure 4.19: Cumulative regret for the copper ore case with context and stationary rewards

This case is the most realistic one in which both the ore type and the process change in
time. The proposed solution is to combine kNN-UCB with a sliding time window. The choice
of the window depends on previous experience of the process. This case is not relevant when
simulations are used since calibrating the simulations would mean knowing variation period.
Figure 4.20 indicates that, despite the non-stationary situation, SW-KNN-UCB outperforms
the other policies as well.

53
Figure 4.20: Cumulative regret for the copper ore case with context and non-stationary re-
wards

To select the appropriate method in each case, a decision tree is constructed, compare
4.21. Based on the distinction between ore variability and process variability, it is a practical
guide to the deployment of the algorithms in an industrial context.The computational time
grows from UCB to kNN-UCB methods since neighborhoods must be calculated. If a very
large number of contextual points are used, it might be useful to reduce the maximum num-
ber of neighbors used by the kNN step. Since only one block is used per time step, the history
of the system will grow as well. This will result in a large memory allocation to compute the
mean rewards for each step in a given neighborhood. For this reason, a sliding-window can
be used when the process is stationary as a way to limit the time horizon of the algorithm and
reduce the memory footprint.

54
 Sequential
Problem allocation
problem

Ore class Targeted Variable

Process Steady-state Dynamic Steady-state Dynamic

Algorithm UCB SW-UCB kNN-UCB SW-kNN-UCB

Figure 4.21: Selection of the right method depending on the conditions

55
Chapter 5

Conclusion

The different steps presented in this work lead from quantitative characterization to decision-
making in a geometallurgical model. The combination of these steps has proved to be an
effective way to use available data from different sources for enhanced characterization,
whereas decision models have provided algorithms for finding optimal production strate-
gies even if the data is noisy or if parts of the information are missing.
Textural information can be acquired at any scale, however, it does not guarantee that tex-
tural features are similar at different scales. From a geological and mineralogical point of
view, large scale textures have their origins in the assemblage of minerals at micro-scale. But
if the association index or grain size is measured from a mineral map obtained from instru-
ments investigating different scales (e.g. by SEM and drill core scanner for example), spatial
resolution of the sensors and wavelength of the probing rays will produce different results.
This challenge can be solved by carefully reviewing the assumptions that are made when
one builds a mineral map from digital signals. In that case, for instance, the integral range
(Lantuéjoul, 2013) for example can produce an empirical law of change of scale. The second
computational method is to use wavelets since they explicitly take scale into account. The
main challenge is to capture both linear and non-linear scale interactions. CNNs that com-
bine convolution filters, pooling operations and non-linear activation functions, regularize
complex data by learning these inter- actions and linearize the groups of transformation. In
mineralogy, groups of symmetry play a key role in the classification of minerals. What CNNs
do follows the same idea, i.e. from an image of a drill core or a mineral map, it learns the sym-
metries. This involves the group of all the transforms (translation, rotation, scale and small
deformations) that do not change the geometallurgical class of the image. To summarize,

56
CNNs are not only a model because they produce a quantitative representation of textures
and capture their variability. They are a predictive model since during training they can be
asked to produce a classification into geometallurgical ore classes or a regression on geomet-
allurgical variables. Moreover, they can use images and numerical as inputs. These models
are a natural link between the geologists, the mineralogist and the process metallurgist by
using mineralogy, chemistry and images to make predictions. Of course, the limitations in-
herent to supervised models are present in the fusion model presented. This refers to the risk
of overfitting, the large amount of training data needed but also the lack of theoretical guide-
lines regarding the architecture. One solution to the problem could be to use Auto-encoders
that are unsupervised versions of CNNs. They do not require human expertise to be trained
since their output is the same as their training input. The trick is to use a low dimensional
layer at the center of the net (bottleneck) and train the net to reproduce its input from this
lower dimensional data (Hinton and Salakhutdinov, 2006). Once information from differ-
ent sources has been combined and used in a process model either by making predictions
about the performance using for example grade, recovery of the mineral of interest in the
concentrate, throughput of the feed material and environmental indicators such as the mass
of CO 2 -equivalent produced per ton of metal or concentrate , the difficulty is to take the best
decisions based on the available data.
The action of CNNs is beneficial by reducing the dimension to 2D or 3D, which makes pos-
sible and an interpretation in terms of mineralogy and lithology. Additionally, it addresses
the issue of variability by smoothing it, thus ensuring that the assumptions for kNN-UCB are
verified but also providing a metric for textural data. Some of these networks can be designed
specifically for interpolation based on that smooth lower-dimensional space. For example
it is possible not only to compute a distance between geometallurgical classes but also to
interpolate between them, offering new perspectives for the population of a block model by
interpolating between classes directly.
Different concepts and the appropriate tools have been developed and tested to reach the
objectives of this thesis. As a summary of the conducted work, the following answers to
the research questions can be given. With respect to textures and their classification and
quantification, it was observed that

• A texture in the context of geometallurgy is the association of different minerals in space.

While this definition is broad, it eliminates the need for an explicit scale and connects

57
 different disciplines such as spatial statistics, image analysis, mineralogy, geology and
mineral processing. Once assumptions are made, the computational tools become
distinct but the fundamental properties of a texture remain the same

• Textures can be characterized by their shape, size and association information. Efficient
computational tools have been developed for each property. A convenient object to
manipulate mineral textures is the mineral map. It is a model based on measurements
that can be used to estimate textural characteristics but also to simulate the production
of particles. These particles can then be used for process simulation.

• High-dimensional data can be represented by using textural features or convolutional

neural networks. The advantage of convolutional networks is their ability to learn ap-
propriate features with or without supervision. They are used for classification or the
prediction of numerical variables and can accept different types of data as inputs.

• Once textures have been quantified and projected, it is possible to define a metric space.
A distance between two different samples or images can be computed and a topology
emerges. The space of projected features is an intermediate representation between
the spatial model and the process model. Finding the appropriate projections and
understanding the topology of textures are the key to future industrial applications in
geometallurgy.

Regarding the process modelling for geometallurgy, the following answers can be given about
a recommendation system using the available information to choose from different process-
ing routes in a sequential context.

• Particles are the link between the mine and the beneficiation process and tools devel-
oped to study their shape, size and mineral association are the same for intact rocks.
Depending on the level of details, models based on mineral particles should be pre-
ferred for geometallurgy applications

• Process models can use information from different sources involving the mineralogy,
grindability, particle size and liberation. As a result, once the models have been cali-
brated, the population of mineral particles can be estimated for any stream of a pro-
cessing circuit.

58
 • Variability description can be separated between the ore and the process. If the prop-
erties of the ore and the operating variables of the process are considered as random
variables, statistical tools are relevant. Texture classification reduces the variability
of the ore by projecting it to a smoother space of lower dimension, while stochastic
bandits handle the variability of the process by estimating its average performance.

• To support enhanced decisions in changing contexts, several policies and numerical

methods exist. Among these, UCB strategies offer strong theoretical guarantees. The
algorithms have been implemented and successfully applied to geometallurgical data.

• The hypotheses required by UCB methods are not very restrictive. In general, a reward
function with values in [0, 1] can be used. When the context is taken into account, then
the reward functions must be smooth.

By the successive integration of data from textures to the market performance, efficient com-
putational methods have been introduced and new perspectives have been opened for ge-
ometallurgy.

59
Bibliography

Aasly, K. and Ellefmo, S. (2014). Geometallurgy applied to industrial minerals operations.

Mineralproduksjon, 5:21–34.

Aggarwal, C. C. (2015). Data classification : algorithms and applications. Chapman &

Hall/CRC.

Agrawal, R. (1995). Sample mean based index policies by O(log n) regret for the multi-armed
bandit problem. Advances in Applied Probability, 27(4):1054–1078.

Ahonen, T., Hadid, A., and Pietikäinen, M. (2006). Face description with local binary pat-
terns: Application to face recognition. IEEE Transactions on Pattern Analysis and Machine
Intelligence, 28(12):2037–2041.

Alruiz, O., Morrell, S., Suazo, C., and Naranjo, A. (2009). A novel approach to the geometallur-
gical modelling of the Collahuasi grinding circuit. Minerals Engineering, 22(12):1060–1067.

Andrusiewicz, M., Edraki, M., Tungpalan, K., Manlapig, E., Wightman, E., Keeney, L., An-
drusiewicz, M., Keeney, L., Wightman, E., and Edraki, M. (2014). An integrated approach
of predicting metallurgical performance relating to variability in deposit characteristics.
Minerals Engineering, 71:49–54.

Antonini, M., Barlaud, M., Mathieu, P., and Daubechies, I. (1992). Image Coding Using
Wavelet Transform. IEEE transactions on image processing : a publication of the IEEE
Signal Processing Society, I(2):194–205.

Auer, P. (2002). Using confidence bounds for exploitation-exploration trade-offs. Journal of

Machine Learning Research, 3(Nov):397–422.

60
Bandyopadhyay, S. and Pal, S. K. (1991). Classification and Learning Using Genetic Al-
gorithms. In Proceedings of the International Joint Conference on Artificial Intelligence,
volume 2, pages 651–656.

Barbery, G. (1991). Mineral liberation: measurement, simulation and practical use in mineral
processing. Québec: Éditions GB.

Barbery, G. and Leroux, D. (1988). Prediction of particle composition distribution after

fragmentation of heterogeneous materials. International Journal of Mineral Processing,
22(1-4):9–24.

Bilal, D. (2017). Master thesis: Geometallurgical estimation of comminution indices for

porphyry copper deposit applying mineralogical approach. Master thesis, Luleå University
of Technology.

Boliden (2018). Boliden Annual and Sustainability Report 2018. Technical report, New
Boliden AB.

Bonnici, N. K. (2012). Doctoral thesis: The mineralogical and textural characteristics of

copper-gold deposits related to mineral processing attributes. Doctoral thesis, University of
Tasmania.

Boser, B. E., Guyon, I. M., and Vapnik, V. N. (1992). A Training Algorithm for Optimal Mar-
gin Classifiers. In Proceedings of the Fifth Annual Workshop on Computational Learning
Theory, COLT ’92, pages 144–152, New York, NY, USA. ACM.

Bottou, L. (2010). Large-scale machine learning with stochastic gradient descent. In Proceed-
ings of COMPSTAT’2010, pages 177–186. Springer.

Breiman, L. (1996). Out-of-bag estimation. Technical report, University of California, Berke-

ley, CA.

Bubeck, S., Cesa-Bianchi, N., Bubeck, S., and Cesa-Bianchi, N. O. (2012). Regret Analysis of
Stochastic and Nonstochastic Multi-armed Bandit Problems. Foundations and Trends in
Machine Learning, 5(1):1–122.

Burnetas, A. N. and Katehakis, M. N. (1997). Optimal adaptive policies for Markov decision
processes. Mathematics of Operations Research, 22(1):222–255.

61
Chiles, J.-P. and Delfiner, P. (2009). Geostatistics: modeling spatial uncertainty, volume 497.
John Wiley & Sons.

Davis, E. W. (1921). Magnetic concentration of iron ore. Mines Experiment Station Bulletin,
XXIV(9):56–59.

Donskoi, E., Suthers, S. P., Campbell, J. J., and Raynlyn, T. (2008). Modelling and optimization
of hydrocyclone for iron ore fines beneficiation - using optical image analysis and iron ore
texture classification. International Journal of Mineral Processing, 87:106–119.

Eriksson, B. and Hallgren, U. (1975). Description of the geological maps Vittangi NV, NO, SV,
SO. Sver. Geol. Unders., page 203.

Frietsch, R. (1966). Berggrund och malmer I Svappavaarafältet, Norra Sverige. Technical

Report 604, SGU, Stockholm.

Garivier, A. and Cappé, O. (2011). The KL-UCB Algorithm for Bounded Stochastic Bandits
and Beyond. In Sham Kakade, U. v. L., editor, 24th Annual Conference on Learning Theory.
JMLR : Workshop and Conference Proceedings, volume 19, pages 359–376.

Gay, S. L. (1999). Numerical verification of a non-preferential-breakage liberation model.

International Journal of Mineral Processing, 57(2):125–134.

Gay, S. L. (2004). Simple texture-based liberation modelling of ores. Minerals Engineering,

17(11):1209–1216.

Guntoro, P. I., Ghorbani, Y., Koch, P.-H., and Rosenkranz, J. (2019). X-ray Microcomputed
Tomography (µCT) for Mineral Characterization: A Review of Data Analysis Methods. Min-
erals, 9(3):183.

Hall, M., Frank, E., Holmes, G., Pfahringer, B., Reutemann, P., and Witten, I. H. (2009). The
WEKA data mining software. SIGKDD Explorations Newsletter, 11(1):10.

Haralick, R. M., Shanmugam, K., and Dinstein, I. (1973). Textural features for image classifi-
cation. IEEE Transactions on Systems, Man and Cybernetics, smc 3:610–621.

Hebb, D. O. (1949). Organization of Behavior. Wiley, New York.

62
Hegenbart, S. and Uhl, A. (2015). A scale- and orientation-adaptive extension of Local Binary
Patterns for texture classification. Pattern Recognition, 48(8):2633–2644.

Hinton, G. E. and Salakhutdinov, R. R. (2006). Reducing the dimensionality of data with

neural networks. Science, 313(5786):504–507.

Katznelson, Y. (2004). An introduction to harmonic analysis. Cambridge University Press.

Kaufmann, E., Korda, N., and Munos, R. (2012). Thompson Sampling: An Asymptotically
Optimal Finite-Time Analysis BT - Algorithmic Learning Theory. In Bshouty, N. H., Stoltz,
G., Vayatis, N., and Zeugmann, T., editors, Proceedings of the International Conference on
Algorithmic Learning Theory, pages 199–213, Berlin, Heidelberg. Springer Berlin Heidel-
berg.

Keeney, L. (2010). Doctoral thesis: The Development of a Novel Method for Integrating Geomet-
allurgical Mapping and Orebody Modelling. Doctoral thesis, University of Queensland.

King, R. (1979). A model for the quantitative estimation of mineral liberation by grinding.
International Journal of Mineral Processing, 6(3):207–220.

King, R. (2001). Modeling and simulation of mineral processing systems. Reed Educational
and Professional Publishing Ltd, Woburn, 1 edition.

Kingma, D. P. and Ba, J. (2014). Adam: A method for stochastic optimization. arXiv preprint
arXiv:1412.6980.

Kłopotek, M. A. (2017). Machine Learning Friendly Set Version of Johnson-Lindenstrauss

Lemma. Technical report, Institute of Computer Science of the Polish Academy of Sciences.

Koch, P.-H. (2017). Licentiate thesis: Particle generation for geometallurgical process modeling.
Monograph, Luleå University of Technology.

Koch, P.-H. (2018). A numerical study of the effects of microwave pre-treatment on value
liberation from a zinc ore. In Proceedings of the 29th International Mineral Processing
Congress, Moscow, Russian Federation.

Koch, P.-H., Lamberg, P., and Rosenkranz, J. (2015). How to build a process model in a
geometallurgical program? In Andre-Mayer, A., Cathelineau, M., Muchez, P., Pirard, E., and
Sindern, S., editors, Proceedings of the 13th SGA Biennial meeting, page 85, Nancy, France.

63
Koch, P.-H., Lund, C., Lishchuk, V., Lamberg, P., and Rosenkranz, J. (2017). Automated
classification of drill cores textures for geometallurgy. In Proceedings of Process Mineralogy
2017, Cape Town, South Africa.

Koch, P.-H., Lund, C., and Rosenkranz, J. (2019a). Automated drill core mineralogical charac-
terization method for texture classification and modal mineralogy estimation for geomet-
allurgy. Minerals Engineering, 136:99–109.

Koch, P.-H., Lund, C., and Rosenkranz, J. (2019b). Data fusion for enhanced prediction of
geometallurgical variables.

Koch, P.-H. and Rosenkranz, J. (2017). Texture-based liberation models for comminution. In
Proceedings of the Conference in Minerals Engineering, pages 83–96, Luleå, Sweden.

Koch, P.-H. and Rosenkranz, J. (2019). Sequential decision-making in geometallurgy. pages

1–22.

Koch, P.-H., Tiu, G., Semsari, P., Parian, M., and Lishchuk, V. (2019c). Textures in geometal-
lurgy : from characterization to applications. pages 1–28.

Lai, T. and Robbins, H. (1985). Asymptotically Efficient Adaptive Allocation Rules. Advances
in Applied Mathematics, 6:4–22.

Lamberg, P. (2011). Particles – the bridge between geology and metallurgy. In Proceedings of
the Conference in Minerals Engineering (CME), pages 1–16, Luleå.

Lantuéjoul, C. (2013). Geostatistical simulation: models and algorithms. Springer Science &
Business Media.

Lawson, K. (2017). Industry Questions&Answers : Data mining. resourceful, 1(12):8 – 9.

LeCun, Y., Bengio, Y., and Others (1995). Convolutional networks for images, speech, and
time series. In Arbib, M. A., editor, The handbook of brain theory and neural networks, page
1995. MIT Press.

Li, W., Chen, C., Su, H., and Du, Q. (2015). Local Binary Patterns and Extreme Learning
Machine for Hyperspectral Imagery Classification. IEEE Transactions on Geoscience and
Remote Sensing, 53(7):3681–3693.

64
Lin, I. J. (1989). The effect of seasonal variations in temperature on the performance of
mineral processing plants. Minerals Engineering, 2(1):47–54.

Lishchuk, V., Koch, P.-H., Lund, C., and Lamberg, P. (2015a). The geometallurgical framework.
Malmberget and Mikheevskoye case studies. In Proceedings of the XVth Conference of PhD
students and Young Scientists, pages 57–66, Szlarska Poreba (Poland).

Lishchuk, V., Lamberg, P., and Lund, C. (2015b). Classification of geometallurgical programs
based on approach and purpose. 13th Biennial SGA Meeting. Mineral resources in a sus-
tainable world. Volume 4., pages 1431–1434.

Lishchuk, V., Lund, C., Koch, P.-H., Gustafsson, M., and Pålsson, B. I. I. (2018a). Geomet-
allurgical characterization of Leveäniemi iron ore - Unlocking the patterns. Minerals
Engineering, 131(November 2018):325–335.

Lishchuk, V., Lund, C., Koch, P.-H., Pålsson, B. I., and Gustafsson, M. (2018b). Geometallurgi-
cal characterisation of Leveäniemi iron ore – Unlocking the patterns. Minerals Engineering,
131:325–335.

LKAB (2018). LKAB’s annual and sustainability report 2018. Technical report, LKAB, Luleå.

Lund, C. (2013). Doctoral thesis: Mineralogical, chemical and textural characterisation of the
Malmberget iron ore deposit for a geometallurgical model. composite, Luleå University of
Technology.

Lund, C. and Lamberg, P. (2014). Geometallurgy–A tool for better resource efficiency. Euro-
pean Geologist, 37(May):39–43.

Lund, C., Lamberg, P., and Lindberg, T. (2013). Practical way to quantify minerals from chem-
ical assays at Malmberget iron ore operations–An important tool for the geometallurgical
program. Minerals Engineering, 49(49):7–16.

Mallat, S. (1989a). A theory for multiresolution signal decomposition: the wavelet represen-
tation. IEEE transactions on pattern analysis and machine.

Mallat, S. (1989b). Multifrequency Channel Decomposition of Images and Wavelet Model.

IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37(no. 12):2091.

65
Matheron, G. (1970). La théorie des variables regionalisées et ses applications. Technical
report, ENSP, Paris.

Matthews, B. (1975). Comparison of the predicted and observed secondary structure of T4

phage lysozyme. Biochimica et Biophysica Acta (BBA) - Protein Structure, 405(2):442–451.

McQuiston, F. and Bechaud, L. (1968). Metallurgical Sampling and Testing. Surface Mining.
New York: The American Institute of Mining, Metallurgical, and Petroleum Engineers, pages
103–121.

Mwanga, A., Rosenkranz, J., and Lamberg, P. (2016). Development and experimental
validation of the Geometallurgical Comminution Test (GCT). Minerals Engineering,
108(November):109–114.

Mwanga, A.-R. (2016). Doctoral thesis: Development of a geometallurgical testing framework

for ore grinding and liberation properties. Doctoral thesis, Luleå University of Technology,
Luleå.

Ojala, T., Pietikainen, M., and Maenpaa, T. (2002). Multiresolution gray-scale and rotation
invariant texture classification with local binary patterns. IEEE Transactions on Pattern
Analysis and Machine Intelligence, 24(7):971–987.

Outotec (2019). HSC Chemistry 9 - Chemical Reaction and Equilibrium Software with Exten-
sive Thermochemical Database and Flowsheet Simulation.

Parian, M., Lamberg, P., and Rosenkranz, J. (2016). Developing a particle-based process
model for unit operations of mineral processing–WLIMS. International Journal of Mineral
Processing, 154:53–65.

Pérez-Barnuevo, L., Lévesque, S., and Bazin, C. (2018). Drill core texture as geometallurgical
indicator for the Mont-Wright iron ore deposit (Quebec, Canada). Minerals Engineering,
122(February):130–141.

Pérez-Barnuevo, L., Pirard, E., and Castroviejo, R. (2013). Automated characterisation of

intergrowth textures in mineral particles. A case study. Minerals Engineering, 52:136–142.

Pérez-Ortiz, M., Jiménez-Fernández, S., Gutiérrez, P., Alexandre, E., Hervás-Martínez, C., and
Salcedo-Sanz, S. (2016). A Review of Classification Problems and Algorithms in Renewable
Energy Applications. Energies, 9(8):607.

66
Reeve, H. W. J., Mellor, J., and Brown, G. (2018). The k -Nearest Neighbour UCB Algorithm
for Multi-Armed Bandits with Covariates. In JMLR: Workshop and Conference Proceedings,
pages 1–29.

Rosenblatt, F. (1961). Principles of Neurodynamics: Perceptrons and the Theory of Brain

Mechanisms. Technical report, Cornell Aeronautical Lab Inc Buffalo NY, Buffalo, NY.

Runge, I. C. (1998). Mining economics and strategy. SME.

Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., Preibisch, S.,
Rueden, C., Saalfeld, S., Schmid, B., Tinevez, J.-Y., White, D. J., Hartenstein, V., Eliceiri, K.,
Tomancak, P., and Cardona, A. (2012). Fiji: an open-source platform for biological-image
analysis. Nature Methods, 9(7):676–682.

Schouwstra, R., De Vaux, D., Muzondo, T., and Prins, C. (2013). A geometallurgical approach
at Anglo American Platinum’s Mogalakwena operation. In MEI, editor, Proceedings of the
second AusIMM international geometallurgy conference, pages 85–92, Brisbane, Queens-
land, Autralia. The Australasian Institute of Mining and Metallurgy (AusIMM).

Serra, J. P. and Soille, P. (1994). Mathematical Morphology and Its Applications to Image
Processing, volume 2.

Singh, K. (2017). Master thesis: A Geometallurgical Forecast Model For Predicting Concen-
trate Quality in WLIMS Process for Leveäniemi Ore. Master thesis, Luleå University of
Technology.

Thompson, W. R. (1933). On the likelihood that one unknown probability exceeds another in
view of the evidence of two samples. Biometrika, 25(3/4):285–294.

Tiu, G. (2017). Master thesis : Classification of Drill Core Textures for Process Simulation in
Geometallurgy: Aitik Mine, New Boliden. Master thesis, Luleå University of Technology,
Luleå.

Tungpalan, K., Wightman, E., and Manlapig, E. (2015). Relating mineralogical and textural
characteristics to flotation behaviour. Minerals Engineering, 82:136–140.

van den Boogaart, K. and Tolosana-Delgado, R. (2018). Predictive Geometallurgy: An Interdis-

ciplinary Key Challenge for Mathematical Geosciences. In B.S. Daya, S., Qiuming, C., and

67
 Frits, A., editors, Handbook of Mathematical Geosciences: Fifty Years of IAMG, chapter 33,
pages 673–686. Springer International Publishing, Cham.

van der Maaten, L. and Hinton, G. (2008). Visualizing data using t-SNE. Journal of machine
learning research, 9(Nov):2579–2605.

Vann, J., Jackson, J., Coward, S., and Dunham, S. (2011). The Geomet Curve – A Model
for Implementation of Geometallurgy. In Proceedings of the 1st AusIMM International
Geometallurgy Conference, pages 35–43, Brisbane, Queensland, Autralia.

Vos, F. (2017). Doctoral thesis: The effect of mineral grain textures at particle surfaces on
flotation response. Doctoral thesis, University of Queensland.

Vos, K. (2018). Master thesis: Medium term planning with Evolution Comparison with current
method at LKAB Svappavaara. Master thesis, Luleå University of Technology.

Wanhainen, C. (2005). On the origin and evolution of the Palaeoproterozoic Aitik Cu-Au-Ag
deposit, Northern Sweden. PhD thesis, Luleå University of Technology.

Wanhainen, C., Billström, K., and Martinsson, O. (2006). Age, petrology and geochemistry of
the porphyritic Aitik intrusion, and its relation to the disseminated Aitik Cu-Au-Ag deposit,
northern Sweden. GFF, 128(4):273–286.

Wanhainen, C., Broman, C., and Martinsson, O. (2003). The Aitik Cu-Au-Ag deposit in
northern Sweden: A product of high salinity fluids. Mineralium Deposita, 38:715–726.

Wanhainen, C., Broman, C., Martinsson, O., and Magnor, B. (2012). Modification of a Palaeo-
proterozoic porphyry-like system: Integration of structural, geochemical, petrographic,
and fluid inclusion data from the Aitik Cu–Au–Ag deposit, northern Sweden. Ore Geology
Reviews, 48:306–331.

Wanhainen, C. and Martinsson, O. (1999). Geochemical characteristics of host rocks to the

Aitik Cu-Au deposit. In C. J. Stanley, editor, Proceedings of the 5th biennial SGA meet-
ing and the 10th quadrennial IAGOD symposium, pages 1443–1446, Rotterdam. Balkema
Publishers, A.A. / Taylor & Francis The Netherlands.

68
Appended papers
Paper I
 Automated drill core mineralogical characterization method for
texture classification and modal mineralogy estimation for
Geometallurgy
Pierre-Henri Koch*, Cecilia Lund, Jan Rosenkranz

* Corresponding author: [email protected]

Luleå University of Technology

Division of Minerals and Metallurgical Engineering (MiMeR)

Luleå, Sweden

Abstract
In geometallurgy, a process model operating at the mineral liberation level needs quantitative textural
information about the ore. The utilization of this information within process modeling and simulation will
increase the quality of the predictions.

In this study, descriptors derived from color images and machine learning algorithms are used to group
drill core intervals into textural classes and estimate mineral maps by automatic pixel classification.
Different descriptors and classifiers are compared, based on their accuracy and capacity to be automated.
Integration of the classifier approach with mineral processing simulation is also demonstrated. The
quantification of textural information for mineral processing simulation introduced new tools towards an
integrated information flow from the drill cores to a geometallurgical model.

The approach has been verified by comparing traditional geological texture classification against the one
obtained from automatic methods. The tested drill cores are sampled from a porphyry copper deposit
located in Northern Sweden.

1. Introduction
Drill core logging by geologists commonly includes a qualitative description of the lithology, mineralogy
and ore textures. It is not only a cost and time consuming method, it also relies on each geologist
knowledge and experience (Bonnici et al., 2008). With recent technological advances in drill core scanning
and industrial devices (e.g. HyLogger and SisuROCK) fast identification of some of the mineral phases to
generate mineral maps at macro scale has been facilitated. This information gives the opportunity to
quantify this textural information for comparison and classification purposes (Schodlok et al., 2016;
Specim, 2013).

Textural information is critical for detailed modeling of the downstream process. The impact of particles
of different textures on hydrocyclones classification, sintering and ultrasonic treatment of iron ore has
been studied by Donskoi et al. (2016). Other studies show the impact of texture on froth flotation (Bushell,
2012; Tungpalan et al., 2015).

Quantified mineralogical information (modal mineralogy and mineral texture) can be obtained at different
scales (on particles or on intact ore). On particles, methods to quantify the mineralogy by using the routine
chemical assay to convert the elemental components into mineral grade have been developed and tested
on iron ore operations (Lund et. al 2013). Common methods to study textures at the particle level include

1
scanning electron microscopy (SEM) and optical microscopy. Textures have been analysed on optical
images by Donskoi et al. (2007). For drill cores, several measurement methods exist for elemental analyses
for example, X-Ray fluorescence (XRF) (Estrada Calderon and Zagerholm, 2012). On intact ore on the other
hand, mineral (micro-) textures are difficult to measure and quantify. Moreover, high-resolution
hyperspectral scans still do not provide adequate additional mineralogical information for the technique
to be implemented at a mine production scale due to its reliance on human interpretation. Therefore it is
associated with an interpretation bias (Arqué Armengol, 2015). A new logging methodology is therefore
necessary to facilitate the daily tasks of geologists which also allow cost-effective, non-destructive and
fast logging of large intervals of drill cores.

The development of technology in image acquisition and image processing combined with the increasing
computing power enables the integration of a wide range of sensors in drill core scanning devices.
However, the acquired signals should be relevant to the problem at hand. A geometallurgical model
consists of two main parts: a spatial model and a process model. Drill core samples are adequate for
geometallurgy since they are collected in any case for exploration and resource estimation purposes.
Different small-scale tests can be conducted on such samples, for example in comminution (Mwanga,
2014) but also in concentration (Vos et al. 2014).

Texture is among the parameters that influence the processing behavior of an ore and has been carefully
studied in different deposits (Bonnici, 2012; Lund et al., 2015; Pérez-Barnuevo et al., 2013; Tungpalan et
al., 2015).

A general approach based on particles and their mineralogy has been described by Lamberg (2011). The
idea is to use the information from the geological model (spatial model) to generate particles and simulate
the comminution and concentration process to forecast the products. The results can be converted into
performance indicators relevant to a specific deposit and plant. These indicators can then be stored in the
spatial model as a tool for mine planning and process design. The entire chain forms a geometallurgical
model as shown in Figure 1.

Performance
Particles
indicators
•Spatial model •Process
•Texture model •Production
breakage forecast
model
Mineralogy
Particles behavior
Textures

Figure 1 : Geometallurgical model based on particles (Lamberg, 2011), modified

In this context, the present work aims at using drill cores to acquire information about texture. Recently,
interest in automated drill core scanning has increased, with extensions towards 3D textural analysis
(Becker et al., 2016), and correlations between the textural class and their behavior in the process (Nguyen
et al. 2016, Pérez-Barnuevo et al., 2016).

2
A systematic approach to texture classification is developed and used to study the potential of several
textural descriptors (or features) and classifiers in an automated logging process.

2. Material and methods

2.1 Data set and software
Images of the original dataset are taken with a Canon Pro 815 camera and a pixel size of 0.12 mm2, with
a stable light source and images are stored as RAW files then converted to TIFF files. For each ore type
between 5 and 10m of drill cores from different parts of the ore body are scanned to obtain representative
color images of each type. The scale is kept constant.

The data set is then randomly cropped and rotated to generate images of different sizes and 3 color
channels for feature extraction and classification purposes. The complete dataset includes 1509 images
in 15 textural classes. In general, the definition of classes depends on the deposit and its variability. The
samples have been collected from a selection of drill cores provided by the mine geologists. To ensure
that the dataset is representative of the textural variations in the ore body, the selected samples were
further studied and compared to previous work on Aitik deposit (Wanhainen et al., 2012, 2006).

For image processing purposes, Fiji (Schindelin et al., 2012) is used. WEKA (Hall et al., 2009), Python and
MATLAB (Mathworks, 2016) are used for the machine learning part of the work. HSC Sim 9 (Outotec, 2016)
is used for process simulation.

2.2 Description of the textural classes

The selected case study is the Aitik Cu-Au-Ag mine, operated by New Boliden, located south of Gällivare,
northern Sweden, which is the largest copper mine in Europe. Aitik operates with a head grade of 0.23%
Cu, 1.8 g/t Ag and 0.13 g/t Au with chalcopyrite as the main copper mineral (Wanhainen et al., 2012).
Previous studies have described and classified in detail all the rock units of the Cu-Au-Ag deposit
(Wanhainen et al., 2012, 2006) to unravel the genetic evolution. Each lithology is described in terms of
mineralogy and qualitative textural description such as grain size, shape and mineral association.

Spatially distributed drill core samples of all rock units in the mine were sampled and mapped based on
their geological features and their estimated metallurgical performance. Fifteen textural classes
representing the hanging wall, ore zone and footwall are, on a macroscopic scale, classified into different
textural and metallurgical classes giving broad ranges of textures for each rock unit (Table 1).
Table 1 : Rock units from Wanhainen et al. (2006, 2012) divided into geometallurgical textural classes

Mine units Rock units/lithology Geometallurgical textural classes

Footwall Feldspar-biotite-amphibole gneiss 3,5

Quartz monzodiorite 10,11,12,13,14,15

Felsic intrusive rocks Micro-quartz monzodiorite 10, 12

Pegmatite veins 1

Ore zone Garnet-bearing biotite schist 7,8,9

Quartz-muscovite-sericite schist 4,6

3
 Hanging wall Feldspar-biotite-amphibole gneiss 2

Mafic intrusion rock Meta-gabbro -

Geometallurgical classes have been defined by grouping samples with a similar mineralogy and Bond work
index (Figure 2). More details and alternative classifications can be found in Tiu (2017).

Estimate of BWi ± 2σ
15.00 14.70

14.00
Bond Ball mill work Index [kWh/t]

13.00

12.00 11.34 11.45

10.70 10.76
11.00 10.31
10.05
10.00 9.41 9.17
9.00 8.31 8.34 8.48
7.99
8.00

7.00

6.00

5.00
1 2 3 4 6 8 9 10 11 12 13 14 15
Geometallurgical class

Figure 2 : Bond Work Index values for each geometallurgical class

In this case study, two textures, one from the mineralized footwall and one from the main ore zone are
described in detail and used to demonstrate pixel classification (Table 2).
Table 2 : Example of texture classes used for pixel classification.

Textural class 11 Quartz monzodiorite – Footwall

A medium gray-white, with grain size about 0.3–0.6 mm, and a

plagioclase-porphyritic matrix of plagioclase, biotite, K-feldspar, and
quartz. The granoblastic texture presents partial sericitic alteration
and microcline grains are usually intergrown within the matrix.

Chalcopyrite, pyrite or magnetite, occur as veinlets and as

disseminated grains.

1 cm

Textural class 9 Garnet-bearing biotite schist/gneiss – Ore zone

4
 The microcrystalline matrix of feldspar biotite, amphibole and
quartz with a grain size of average < 0.2 mm contains
approximately 40% of plagioclase.

Finely disseminated chalcopyrite and/or larger patches of

chalcopyrite and pyrite are the most common mineralization styles
in the biotite schist/gneiss.

1 cm

2.3 Automated drill core mineralogical characterization process

The first step is to acquire images of halved drill cores with sufficiently high resolution. Features are then
extracted from the image and stored in a data set. Depending on the values of the descriptors, a
classification algorithm assigns a texture class to each image. Once the class is known, the adequate pixel
classifier is chosen, and the mineral phases of the drill core are identified. If the textural class is known
and the modal mineralogy is estimated, the knowledge is sufficient to build a stream file for process
simulation. The simulation software (in this case HSC Sim, v9.0 Outotec) uses the constructed feed to
calculate every stream including the concentrate and tailings stream (in terms of recovery per mineral,
flow rate, particle size distribution). Figure 2 illustrates the procedure.

Image Feature Texture Pixel Feed

acquisition extraction classification classification construction
•Color •Batch process •Fixed size •Segmentation •Modal
calibration •Write datasets patches algorithm mineralogy
•Filtering •Deposit •Mineral •Comminution
•Cropping specific identification properties

Figure 3 : Texture classification process

2.3.1 Feature extraction

Among the most popular textural features extraction methods, two main approaches are selected: grey-
level co-occurrence matrix (GLCM) statistics and local binary patterns (LBP). GLCM descriptors have been
introduced by Haralick et al. (1973) and recently applied to drill core characterization by Becker et al.
(2016). The frequency at which a given pixel value is in a neighborhood (specified by a distance or offset
and a direction) of another is represented as a frequency matrix called co-occurrence matrix. It can be
made rotation invariant by averaging the co-occurrence values in several directions (usually at an angle
of 0, 45, 90, and 135 degrees).

Another class of descriptors is wavelets. As an extension of the frequency analysis performed by Fourier
filter banks, wavelet filters decompose the signal of an image into different frequencies at different scales.
In this implementation, a Haar wavelet multi-scale decomposition (Antonini et al., 1992; Mallat, 1989)
and statistics derived from the GLCM at each scale are combined. The advantage of the method is to
separate different scales present in the initial texture. The resulting feature vector contains the variance,
correlation, energy and homogeneity statistics of the GLCM at each scale of a wavelet decomposition

5
using Haar filters. Physically it represents the association of minerals or phases with each other at different
scales. These features are chosen to consider scale properties of textures.

LBP have been introduced by Ojala and Pietikainen (2002) and have found multiple applications in image
analysis such as in face description (Ahonen and Hadid, 2006) and specifically in texture analysis either as
such (Li et al., 2015) or in a modified version (Hegenbart and Uhl, 2015; Liao et al., 2009). The descriptors
are based on gray-level images and a discrete circular neighborhood. In the case of color images, the LBP
can be computed on the image converted to grey-scale values or on each color channel separately.
Physically it represents the local geometric patterns observed on textures. These features are chosen to
obtain local information that is not sensitive to rotation.

Figure 4 : Representation of the feature space for 3 bins of the local binary patterns’ histograms. The colors represent a class from
2 to 16 (1 being a non-classified sample).

The result of feature extraction is a vector that contains the texture descriptors associated with a given
input image. When several descriptors are combined, the resulting feature vector is the concatenation of
each individual descriptor feature vector. Figure 4 presents an example of LBP features limited to their 3
first coordinates with respect to their textural class.

2.3.2 Classification algorithms

The main purpose of this work is to classify images into textural classes. The dimensionality of images is
very high due to their dimensions and color depth information. The main idea behind classification
algorithms is to find a transformation that reduces the dimension and defines boundaries between
images of different classes.

6
Classifiers are functions that assigns to each feature vector a texture class. Two main families are used:
unsupervised and supervised classifiers. Unsupervised methods rely on the data itself to learn a
classification whereas supervised methods learn a classification based on examples provided by the user.
Supervised classification can be thought as teaching by examples: a sufficient number of example pairs
(feature vector, actual texture class) is provided and the classifier builds a function that assigns a class
number, ranging between 1 and 15, to the feature vector. In general, the obtained models map an input
to an output but cannot always be described by a linear or analytic function. This approach leads to a
classifier model that is trained to reproduce and generalize the behavior of the human operator that
provided the examples (training set).

In practice, unsupervised methods do not guarantee a classification that respects a pre-existing

descriptive classification for example done by geologists.

K-Nearest Neighbors methods (Aggarwal, 2015; Bandyopadhyay and Pal, 1991) are unsupervised and are
based on the calculation of a distance between the feature vectors. Starting with the center of each class
randomly distributed in the feature space, a feature vector is assigned to the class that has the closest
center. The key of the algorithm is to measure distances between vectors in relatively low dimensions
compared to distances between images which is a much more complex task. After one step, the location
of the center of each class is updated. By iteration, each feature vector will be assigned to a class. This
method can be thought as grouping textural images that have similar numerical features.

Random Forests (Breiman, 2001), Support vector machine (Pérez-Ortiz et al., 2016) and artificial neural
networks (Hebb, 1949; White and Rosenblatt, 2006) are part of the supervised methods. These techniques
are useful when examples of correct classification exist already. In that case, the models can learn a
classification function from the data. Random forests are based on the construction of decision trees. Each
decision tree is built by selecting at each node a criterion that maximizes the information gain regarding
the textural class. By doing so, each tree offers a sequence of binary decisions that lead to a full
classification. In random forests methods, the class of a sample is chosen based on a vote of many decision
trees.

Support vector machine (SVM) is based on finding a hyperplane that maximizes the distance between
feature vectors belonging to different classes. Since the boundary between classes can be complex in high
dimensions, a function that maps the feature vector into a space in which the boundary is linear is
computed by optimization (Boser et al., 1992). Once the features have been transformed, planes that
separate each class is built. Each plane is described by its normal vector (support vector) and projecting
each feature on this vector separates the data into classes.

Artificial neural networks (ANN) is a method based on the construction of a function represented as a
series of weighted sums. The weights are organized in layers that defines the order in which the sums are
computed. The input layer has the same number of weights as the components of the feature vector and
the last layer has one weight for each class (15). Several layers can be placed in between and connect the
input to the output, they are called hidden layers. At first, the weights can be random numbers. When an
example data is provided, the difference between the actual class to which the image belongs and the
predicted class by the network will provide and error. As the number of examples increases, the network
will adapt the weights to minimize the classification error.

7
 2.3.3 Pixel classification
Once the acquired image is classified, the set of probable minerals can be reduced. It is thus possible to
classify pixels from the original image into minerals. A machine learning approach is chosen, based on a
random forest classifier and image features. This procedure allows the pixel classifier to learn from the
geologist and improve every time new training data is provided.

Using Fiji/ImageJ image processing software and its trainable Weka segmentation plugin (Arganda-
Carreras et al., 2017) for pixel classification the following steps are taken (Table 3).
Table 3 : Steps to train and validate WEKA trainable pixel classifier on drill core images.

Pixel classifier training and validation procedure

1 Open Fiji/ImageJ
2 Open a digital image in tiff format
3 Start the trainable Weka segmentation plug-in
4 Select the mineral phases used for training by using point, polygon or freehand selection tool.
Use the original drill core to ensure valid training data set
5 Train the model with the defined training areas from step 4
6 Improve the classification by adding new training areas in the image
7 Repeat stage 5 until satisfactory results are achieved
8 Save the model and the training data set
9 Validate the model with a new image of the same texture and load the existing model/classifier
10 Apply the model. If necessary, add or remove a mineral class and re-train the model

3. Experimental results
This section presents the results of the numeric experiments for each of the steps presented in Figure 3.
The focus besides the accuracy of the classification is to evaluate the potential of each method in terms
of automation. The main factors considered are the speed of feature extraction and the speed of
classification once the features have been extracted. The results presented for texture classification are
performed on all classes and all images.

3.1 Feature extraction performances

The capacity of texture indicators to be automated depends on the resources needed to extract features.
Table 4 presents the total time for the whole image data set but also the time taken to extract a feature
vector for each image.
Table 4 : Time performances for the selected feature sets

Method Number of features Total CPU time [s] CPU time [ms/image]
LBP 32 65.22 43.22

HaarLevel2 128 217.98 144.45

HaarLevel3 192 299.23 198.29

LBP+HaarLevel2 160 268.94 178.22

The GLCM features include the wavelet decomposition which increases their CPU time.

8
 3.2 Texture classification
Each data set is tested with each classifier using 15-fold cross validation. The original sample is split into
15 subsamples (of size 100 in our case, one of these is kept as validation data while the rest is used to
train the model. The procedure is repeated 15 times so that all the data has been used for validation. The
results are averaged over the 15 folds to produce a single estimate. This is the list of selected data sets

• LBP: Local binary patterns

• HaarLevel2: GLCM statistics from Haar wavelet decomposition at level 2
• HaarLevel3: GLCM statistics from Haar wavelet decomposition at level 3
• LBP+HaarLevel2: Combination of LBP and HaarLevel2 descriptors

and classifiers:

• RF: Random Forest

• SVM: Support vector machine with sequential minimal optimization (SMO)
• kNN-p2: k-Nearest Neighbors classifier with a Euclidean distance
• kNN-p1: k-Nearest Neighbors classifier with a Manhattan distance
• kNN-p1.5: k-Nearest Neighbors classifier with a Minkowski distance (p=1.5)
• ANN: Artificial neural network with 3 hidden layers

Different metrics can be used to compare classifiers, in this case Matthews correlation coefficient is
presented in Table 5. It is a statistic derived from the matrix that has predicted classes as rows and actual
classes as columns and contains the count of each sample that is classified (confusion matrix). It indicates
the correlation between the predicted class and the observed class (Matthews, 1975).
Table 5 : Matthews correlation coefficient. Values in bold are significantly different at p=0.05

RF SVM kNN-p2 kNN-p1 kNN-p1.5 ANN

LBP 0.97 ± 0.06 0.80 ± 0.14 0.95 ± 0.07 0.95 ± 0.07 0.96 ± 0.07 0.91 ± 0.11

HaarLevel2 0.97 ± 0.06 0.97 ± 0.05 0.90 ± 0.10 0.90 ± 0.10 0.92 ± 0.09 0.95 ± 0.08

HaarLevel3 0.97 ± 0.06 0.97 ± 0.06 0.91 ± 0.09 0.91 ± 0.09 0.92 ± 0.09 0.96 ± 0.07

LBP+HaarLevel2 0.98 ± 0.05 0.98 ± 0.05 0.91 ± 0.09 0.91 ± 0.09 0.94 ± 0.08 0.98 ± 0.05

By adding the time needed to extract features from one image to the time needed to classify that image,
the total time to process one image can be estimated. If these values are transformed into a time score
with value 0 for the longest time and 1 for 0 second of processing, a global score can be built as the
product of Matthews correlation coefficient and the time score. This helps ranking the performance of a
given combination of features and classifiers, the objective being to be as close as possible to 1. The equal
weights on classification and time scores is justified by the balance needed between accuracy and time
needed in an in-line industrial measurement system. For different applications, different scores can be
built by constructing a score function with variable weights or even other statistics.

9
Table 6 : Global score for classifiers and data sets. Values in bold indicate the best scores. See the text for
abbreviations.

Descriptors RF SVM kNN-p2 kNN-p1 kNN-p1.5 ANN

LBP 0.87 0.72 0.85 0.85 0.82 0.82

HaarLevel2 0.64 0.64 0.54 0.56 0.28 0.63

HaarLevel3 0.52 0.52 0.41 0.43 0.00 0.51

LBP+HaarLevel2 0.57 0.57 0.48 0.48 0.25 0.57

The results shown on Table 6 indicate that the Random Forest classifier and local binary patterns is a
suitable combination for classifying Aitik textures with a reasonable time and performance balance.

3.3 Pixel classification

As an example of pixel classification, the two classes introduced in section 2.2 are presented.

3.3.1 Training of pixel classifier

The first pixel classifier is trained using the coarser grained plagioclase-porphyritic quartz monzodiorite
rock type. In total four different minerals are identified. The major minerals are plagioclase, quartz and
biotite occurring as granular matrix cross-cutting by smaller veins of quartz. Minor amount of chalcopyrite
exists as smaller patches. To train the classifier, three different images of the same texture class are used
(Figure 5) until the out-of-bag error is small.

1 cm
Figure 5 : A trained pixel classifier using three different images of quartz monzodiorite textural type.

3.3.2 Verification of single pixel classifier

The stepwise procedure described above is applied to the garnet-bearing biotite schist/gneiss rock type.
Five different images of the same texture are used to train the classifier and the image presented in Figure
6 is used only to verify the trained pixel classifier. The procedure of verification is done visually by a
geologist. While this tells nothing about the absolute accuracy of the method, it makes it reproducible
and less prone to variability. Proper validation is currently investigated by using scanning electron
microscope (SEM) and X-Ray tomography (XRT). In this case, the drill core belongs to the same textural
class as the one used for training the pixel classifier but was not used during training.

10
A matrix of fine-grained feldspar and biotite is identified with larger veins and patches of chalcopyrite,
pyrite, as well as irregularly distributed magnetite (Figure 6). The mineral map produced by the pixel
classifier is in accordance with a naturally formed geological texture for both smaller features, such as
disseminated grains of chalcopyrite and larger veinlets of magnetite.

1 cm

Figure 6 Verification of a pixel classifier trained on five different images of the same texture.

11
Figure 7 : Verification of visible light pixel classification. Top left: SEM-EDS mineralogy. Top right: Visible light estimated
mineralogy. Bottom: Visible light image

The validation of pixel classifier is done by comparing the results from SEM and visible light pixel
classification. On Figure 7 , one can observe that, despite the different pixel resolution and the lower
accuracy of the visible light classification.

The mineralogy estimates from visible light segmentation can be compared to the SEM-EDS
measurements. One can observe on Figure 8 that some minerals are underestimated by segmentation
and reported to the quartz-felspar fraction. The smoothing effect from the lower resolution has an impact
on the accuracy of the estimates. Differences between segmentation results can be compared by
computing the histograms of each image using the number of observed phases as bins and fit a linear
function f(x) = x +b. The result of the fit by least-squares is presented on Figure 9 and shows a low term b
of -0.0188 which supports a good fit except for the quartz-felspar fraction.

12
 Comparison between SEM-EDS and visible light image
estimated mineralogy
0.3
0.25
0.2
Wt %

0.15
0.1
SEM-EDS
0.05
Visible
0

Main minerals

Figure 8 : Comparison between mineralogy estimates from SEM and visible light images from Figure 7
Comparison between mineralogy estimates

0.25

0.2

0.15
Visible [wt%]

0.1

0.05
Visible vs. SEM-EDS
f(x) = x - 0.01138
99 % Confidence bounds
0

0 0.05 0.1 0.15 0.2

SEM-EDS [wt%]

Residuals

0.1

0.08

0.06

0.04
Visible [wt%]

0.02

-0.02

-0.04

0 0.05 0.1 0.15 0.2

SEM-EDS [wt%]

Figure 9 : Quantitative measurement of the difference between the estimated mineralogy by visible light image and SEM-EDS
measurements

13
The results suggest that visible light image can be used for phase identification after textural classification
constrained the mineralogy. Training the segmentation model requires mineralogists’ and geologists’
expertise but once trained, it can be automated. Further validation has been continued in a master thesis
at Luleå University of Technology by Tiu (2017) and the study concluded that drill core photographs offers
a cheap and fast way to produce a mineral map if the training set is large enough and accurate.

4 Discussion
The case study has revealed several critical steps. The first one is the definition of appropriate
geometallurgical classes based on rock types but also process performance (either estimated from
laboratory tests or plant survey). The underlying assumption is that a correlation exists between the meso
texture observed on the drill core surface, its micro texture as seen under a microscope and its behavior
in comminution and concentration. There is a need for a measure of similarity between different scales
of textures.

The second challenge is to define the boundaries of the approach: in this case study, a good knowledge
of the geology is needed as well as variability in textures and enough drill core images to use machine
learning tools. This challenge is directly linked to the process model: each geometallurgical class must be
tested at a level sufficient to build a stream file and run a process simulation based on a feed derived from
texture classification. Validation of the pixel classification algorithm is important since the verification
presented in this work is based on observations presenting the same risk of bias as drill core logging by
visual inspection.

Benchmarking different combinations of textural features and classifier is a tool for decision. The first
requirement is a suitable metric that includes both texture discrimination accuracy and speed
requirements. This study proposes a global score in this study, but other options should be investigated.

5 Towards a geometallurgical model with textures

Once drill core intervals can be classified into textural classes, this information can be used to build a
geometallurgical model. Assuming a reliable model exists for the process and blending strategies are
known, the appropriate stream can be constructed for the feed. The integration of textural information
into the spatial model (block model) is still a challenge but different approaches have been studied in
literature.

Textural classes can be used for process simulation. If chemical assays of the current drill core are
available, the appropriate element-to-mineral conversion (EMC) recipe can be selected (with correct
minerals and the exact mineral compositions if available) (Lund et al. 2013). In that case, the quality of the
modal mineralogy estimate is improved. In absence of chemical assays for every drill core, its textural
class can be used to choose the appropriate pixel classifier and estimate the mineralogy at the surface of
the drill core. Grinding tests results and modal mineralogy provide enough information to simulate the
concentration processes at the 2D level (i.e. minerals by size, Lamberg and Lund, 2012).

Integration of the mineral processing parameters into a spatial model has multiple approaches due to
their variety: recovery, throughput, energy and reagent consumption, degree of liberation, texture etc.
Having different natures, mineral processing parameters or geometallurgical parameters will have
different properties when propagated in the spatial model. These variables can be classified as continuous
(e.g., elemental and mineral grades, Davis tube recovery, A*b, solubility ratio, liberation) or categorical

14
(e.g., lithology, hydrothermal alteration, ore type, textural class). For sampling and modelling purposes it
can be important to distinguish intrinsic properties called primary parameters (e.g., mass, color, specific
gravity) from process-dependent ones, called response (e.g., throughput, recovery, grindability)
parameters. This classification was introduced by Coward et al. (2009) and is known as the “Primary-
response framework”.

Limitations of spatial modelling of the geometallurgical parameters caused by nonlinearity or non-

additivity have been highlighted by multiple scholars (Coward et al., 2009; Deutsch et al., 2016; Keeney
and Walters, 2011; Newton and Graham, 2011; Richmond and Shaw, 2009; Walters, 2011). Therefore,
many (such as Richmond and Shaw (2009), Keeney and Walters (2011), Dunham and Vann (2007), van den
Boogaart et al., 2013) agree that traditional geostatistical methods which are based on weighted averages
have limited application for the geometallurgical block model. Both Newton and Graham (2011) and
Walters (2009) distinguish two possible origins of the additivity problem: mathematical and machine
influence. Boogaart et al. (2013) also mention how non-linearity originates in non-linear scales at which
those parameters are measured.

One possible solution for overcoming non-additivity issue is spatial simulation of the ore body (Carrasco
et al., 2008; Deutsch, 2013; van den Boogaart et al., 2013). The advantage of the simulation approach over
the conventional one is that value calculated from several simulations is an average and therefore is a
more reliable estimate of the nonlinear parameter. In addition, van den Boogaart et al. (2013) suggest to
use compositional and geometric geostatistics.
Table 7 Geometallurgical indices deployment into the 3D spatial model – current state

Parameters Measured/claimed* Method Comments Case Reference

impact study
Throughput (via Open pit optimization, Multivariate geo- Metallurgical South Deutsch et
low A*b and high 300 days of the statistical simulation properties are American al., 2015
BWi) operation forecast using an intrinsic super frequently copper–
(probabilistic). Monte secondary approach, unequally molybdenum
Carlo simulation was Gaussian sequential sampled and porphyry
used. simulation, global sampled at deposit.
kriging, scales much
larger than
typical assay
measurements.
Handling
multivariate data
was a
challenge.

Grade of Gain per ton was Conditional simulations Correct Non-spatial Boogaart et
unspecified optimized modelling of the simulated al., 2013
element/mineral, complete data
liberation grain random field of
sizes geometallurgical
properties is
required.
Material Gain of selling at Co-kriging and - A high-grade Tolosana-
composition selected choice for the geostatistical iron ore Delgado et
(deleterious, minable blocks (sell as simulation deposit al., 2015
hematite, quartz, lump, sell as high
and shale) quality, sell as low
quality)

15
 Parameters Measured/claimed* Method Comments Case Reference
impact study
Multivariate Mine planning, Nearest neighbor, - Cadia East, Newton and
mineralogical scheduling activities indicator kriging, Red Dog Graham,
groups (MMG) stochastic, ordinary 2011
defined via PCA, kriging trend analysis
A*b, BMWi (with and without zone
control for A*b and
BMWi)
Cu recovery, Mine planning, mine Sequential Gaussian Reliance on the Olympic Boisvert et
U3O8 recovery, optimization, and plant simulation with a multivariate Dam al., 2013
acid performance subsequent regression Gaussian
consumption, net optimization* fit to predict plant distribution after
recovery, drop performance. univariate
weight index, and Linear regression transformation
bond mill work model was based on was a limitation
index four super variables factor.
derived from 112
original variables.
Throughput, Net present values, Process simulation - Synthetic Lishchuk et
recovery internal rate of return applied to every block. deposit al., 2016a,
Nearest neighbor, model 2016b
linear prediction for
data propagating.

Most common techniques for propagating geometallurgical parameters into a block model include:
geostatistical interpolation (kriging, inverse distance), geostatistical simulation, regression, machine
learning, domaining, classification or their combinations. Some of those are listed in the Table 7 and
mention relevant case studies. Until recently, not much work had been conducted on deployment of the
textural data in a spatial model. Texture as it is understood in this study, is a multivariate descriptor of the
ore and multivariate predictor for the process simulation. Texture carries information about both primary
and response parameters and that is why it cannot be regarded as another dimension reducing method
(such as principle component analysis method for example). This poses three important requirements to
textural classes. The first one is that texture should include a distinctive description of the intrinsic ore
properties thus being easily separated from other textures either by a trained geologist or by the
automation system. Next, it should represent a significant difference in processing properties. The third
one is to have a representative volume of the defined texture in the ore body. Without fulfilling those
requirements spatial modelling will not yield into a reliable forecast. The easiest way to approach solution
is testing hypotheses on a synthetic deposit models such as proposed by Lishchuk (2016). An advantage
of using synthetic data is the complete knowledge of the simulated ground truth.

6 Conclusion
Comparison of combination of different textural descriptors (features) with different classifiers in terms
of accuracy but also speed (with automation and integration in a geometallurgical program in mind)
suggests that the combination of local binary patterns (features) and random forests (classifier) is suitable
for classifying the selected textures of Aitik. However, further validation on new samples or ore types is
required. Using a new dataset would allow a proper validation of the current models. If needed, the
existing models can be updated until the requirements for industrial application are met.

The proper validation of pixel classifiers would require not only a comparison of the frequency of each
phase but also textural features like association index and grain shape.

16
The approach with a limited number of features and a simple but efficient classifier allows reaching
performances close to a real-time drill core texture classification system. Depending on the objectives of
the industrial application, a score based on a weighted average of the time and classification
performances can be built. The value for geometallurgy is to maximize the information extracted from
samples available at an early stage. Building a stream file for an integrated simulation is possible, however
the link between meso- and micro-texture should be studied in more detail to simulate the process at the
liberation level. The combination of automated drill core scanning with a process and a spatial model
inside a geometallurgical model should be the scope of future studies.

7 Acknowledgments
The authors wish to thank Boliden for their support in this study, Prof. C. Wanhainen and Glacialle Tiu for
interesting discussions about the Aitik deposit and the two anonymous reviewers for their valuable
comments. This work is part of PREP project funded by VINNOVA.

8 References
Aggarwal, C.C., 2015. Data classification : algorithms and applications, Series: Chapman & Hall/CRC data mining and
knowledge discovery series ; 35. Chapman & Hall/CRC.

Ahonen, T., Hadid, A., Pietikäinen, M., 2006. Face description with local binary patterns: Application to face
recognition. IEEE Transactions on Pattern Analysis and Machine Intelligence 28, 2037–2041.
doi:10.1109/TPAMI.2006.244

Antonini, M., Barlaud, M., Mathieu, P., Daubechies, I., 1992. Image Coding Using Wavelet Transform. IEEE
transactions on image processing : a publication of the IEEE Signal Processing Society I, 194–205.
doi:10.1109/83.136597

Arganda-Carreras, I., Kaynig, V., Rueden, C., Eliceiri, K.W., Schindelin, J., Cardona, A., Seung, H.S., 2017. Trainable
Weka Segmentation: a machine learning tool for microscopy pixel classification. Bioinformatics (Oxford
University Press) btx180.

Arqué Armengol, A., 2015. Towards Automated Logging of Ore: Positive Identification of Sulphides in the Ores of
Agnico Eagle Kittilä and New Boliden Mines. Luleå University of Technology.

Bandyopadhyay, S., Pal, S.K., 1991. Classification and Learning Using Genetic Algorithms, in: Proceedings of the
International Joint Conference on Artificial Intelligence. pp. 651–656.

Becker, M., Jardine, M.A., Miller, J.A., Harris, M., 2016. X-ray Computed Tomography: A geometallurgical tool for 3D
textural analysis of drill core?, in: Proceedings of the 3rd AusIMM International Geometallurgy Conference.
pp. 15–16.

Boisvert, J.B., Rossi, M.E., Ehrig, K., Deutsch, C. V., 2013. Geometallurgical Modeling at Olympic Dam Mine, South
Australia. Mathematical Geosciences 45, 901–925. doi:10.1007/s11004-013-9462-5

Bonnici, N., Hunt, J.A., Walters, S.G., Berry, R., Collett, D., 2008. Relating textural attributes to mineral processing:
Developing a more effective approach for the Cadia East Cu-Au porphyry deposit., in: Proceedings of the Ninth
International Congress for Applied Mineralogy Conference (ICAM). pp. 4–5.

Bonnici, N.K., 2012. The mineralogical and textural characteristics of copper-gold deposits related to mineral
processing attributes. University of Tasmania.

Boser, B.E., Guyon, I.M., Vapnik, V.N., 1992. A Training Algorithm for Optimal Margin Classifiers, in: Proceedings of
the Fifth Annual Workshop on Computational Learning Theory, COLT ’92. ACM, New York, NY, USA, pp. 144–
152. doi:10.1145/130385.130401

17
Breiman, L., 2001. Random forests. Machine Learning 45, 5–32. doi:10.1023/A:1010933404324

Bushell, C., 2012. The PGM flotation predictor: Predicting PGM ore flotation performance using results from
automated mineralogy systems. Minerals Engineering 36–38, 75–80. doi:10.1016/j.mineng.2012.02.016

Carrasco, P., Chiles, J.-P., Seguret, S., 2008. Additivity, Metallurgical Recovery, and Grade, in: Proceedings of the 8th
International Geostatistics Congress. Santiago, Chile.

Coward, S., Vann, J., Dunham, S., Stewart, M., 2009. The Primary-Response Framework for Geometallurgical
Variables, in: Proceedings of the Seventh International Mining Geology Conference. Perth, WA, Australia, pp.
109–113.

Deutsch, J.L., Palmer, K., Deutsch, C. V., Szymanski, J., Etsell, T.H., 2016. Spatial Modeling of Geometallurgical
Properties: Techniques and a Case Study. Natural Resources Research 25, 161–181. doi:10.1007/s11053-015-
9276-x

Deutsch, C. V, 2013. Geostatistical Modelling of Geometallurgical Variables – Problems and Solutions, in: Proceedings
of the Second AusIMM International Geometallurgy Conference. Brisbane, pp. 7–15.

Donskoi, E., Poliakov, A., Holmes, R., Suthers, S., Ware, N., Manuel, J., Clout, J., 2016. Iron ore textural information
is the key for prediction of downstream process performance. Minerals Engineering 86, 10–23.
doi:10.1016/j.mineng.2015.11.009

Donskoi, E., Suthers, S.P., Fradd, S.B., Young, J.M., Campbell, J.J., Raynlyn, T.D., Clout, J.M.F., 2007. Utilization of
optical image analysis and automatic texture classification for iron ore particle characterisation. Minerals
Engineering 20, 461–471. doi:10.1016/j.mineng.2006.12.005

Dunham, S., Vann, J., 2007. Geometallurgy , Geostatistics and Project Value — Does Your Block Model Tell You What
You Need to Know?, in: Project Evaluation Conference. Melbourne, Vic, 19-20 June 2007. pp. 19–20.

Estrada Calderon, E., Zagerholm, I., 2012. Development of a system for X-ray analysis within the mining industry.

Hall, M., Frank, E., Holmes, G., Pfahringer, B., Reutemann, P., Witten, I.H., 2009. The WEKA data mining software.
SIGKDD Explorations Newsletter 11, 10. doi:10.1145/1656274.1656278

Haralick, R.M., Shanmugam, K., Dinstein, I., 1973. Textural features for image classification. IEEE Transactions on
Systems, Man and Cybernetics smc 3, 610–621. doi:10.1109/TSMC.1973.4309314

Hebb, D.O., 1949. Organization of Behavior. Wiley, New York.

Hegenbart, S., Uhl, A., 2015. A scale- and orientation-adaptive extension of Local Binary Patterns for texture
classification. Pattern Recognition 48, 2633–2644. doi:10.1016/j.patcog.2015.02.024

Keeney, L., Walters, S.G., 2011. A methodology for geometallurgical mapping and orebody modelling, in: Proceedings
of the 1st AusIMM International Geometallurgy Conference. Brisbane, pp. 217–225.

Lamberg, P., 2011. Particles–the bridge between geology and metallurgy, in: Proceedings of Conference in Minerals
Engineering (CME). Citeseer, Luleå, pp. 1–16.

Lamberg, P., Lund, C., 2012. Taking Liberation Information into a Geometallurgical Model-Case Study Malmberget,
Northern Sweden, in: Proceedings of Process Mineralogy’12. Cape Town.

Li, W., Chen, C., Su, H., Du, Q., 2015. Local Binary Patterns and Extreme Learning Machine for Hyperspectral Imagery
Classification. IEEE Transactions on Geoscience and Remote Sensing 53, 3681–3693.
doi:10.1109/TGRS.2014.2381602

Liao, S., Law, M.W.K., Chung, A.C.S., 2009. Dominant local binary patterns for texture classification. IEEE Transactions
on image processing 18, 1107–1118. doi:10.1109/TIP.2009.2015682

18
Lishchuk, V., 2016. Geometallurgical Programs – Critical Evaluation of Applied Methods and Techniques. Luleå
University of Technology.

Lishchuk, V., Lamberg, P., Lund, C., 2016a. Evaluation of sampling in geometallurgical programs through synthetic
deposit model, in: Proceedings of the XXVI International Mineral Processing Congress (IMPC).

Lishchuk, V., Lund, C., Lamberg, P., 2016b. Development of a synthetic ore deposit model for geometallurgy, in:
Proceedings of the 3rd AusIMM International Geometallurgy Conference. pp. 1–30.

Lund, C., Lamberg, P., Lindberg, T., 2015. Development of a geometallurgical framework to quantify mineral textures
for process prediction. Minerals Engineering 82, 61–77. doi:10.1016/j.mineng.2015.04.004

Mallat, S., 1989. Multifrequency Channel Decomposition of Images and Wavelet Model. IEEE Transactions on
Acoustics, Speech, and Signal Processing 37, 2091.

Mathworks, 2016. MATLAB.

Matthews, B.W., 1975. Comparison of the predicted and observed secondary structure of T4 phage lysozyme.
Biochimica et Biophysica Acta (BBA) - Protein Structure 405, 442–451. doi:10.1016/0005-2795(75)90109-9

Mwanga, A., 2014. Test Methods for Characterising Ore Comminution Behavior in Geometallurgy. Luleå University
of Technology, Luleå.

Newton, M.J., Graham, J.M., 2011. Spatial Modelling and Optimisation of Geometallurgical Indices, in: Proceedings
of the 1st AusIMM International Geometallurgy Conference. pp. 5–7.

Nguyen, A., Jackson, J., Nguyen, K., Manlapig, E., 2016. A new semi-automated method to rapidly evaluate the
processing variability of the ore body, in: Proceedings of the 3rd International Geometallurgy Conference.
Perth, WA, pp. 145–151.

Ojala, T., Pietikainen, M., Maenpaa, T., 2002. Multiresolution gray-scale and rotation invariant texture classification
with local binary patterns. IEEE Transactions on Pattern Analysis and Machine Intelligence 24, 971–987.
doi:10.1109/TPAMI.2002.1017623

Outotec, 2016. HSC Sim.

Pérez-Barnuevo, L., Lévesque, S., Michaud, D., Bazin, C., Longuépée, H., 2016. Geometallurgical characterization of
drill core textural patterns: a case study from the Mont-Wright iron ore deposit, in: IMPC 2016: XXVIII
International Mineral Processing Congress Proceedings. Québec.

Pérez-Barnuevo, L., Pirard, E., Castroviejo, R., 2013. Automated characterisation of intergrowth textures in mineral
particles. A case study. Minerals Engineering 52, 136–142. doi:10.1016/j.mineng.2013.05.001

Pérez-Ortiz, M., Jiménez-Fernández, S., Gutiérrez, P., Alexandre, E., Hervás-Martínez, C., Salcedo-Sanz, S., 2016. A
Review of Classification Problems and Algorithms in Renewable Energy Applications. Energies 9, 607.
doi:10.3390/en9080607

Richmond, A., Shaw, W.J., 2009. Geometallurgical Modelling- Quo Vadis?, in: Proceedings of the Seventh
International Mining Geology Conference. pp. 115–118.

Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., Preibisch, S., Rueden, C., Saalfeld,
S., Schmid, B., Tinevez, J.-Y., White, D.J., Hartenstein, V., Eliceiri, K., Tomancak, P., Cardona, A., 2012. Fiji: an
open-source platform for biological-image analysis 9, 676–682. doi:10.1038/nmeth.2019

Schodlok, M.C., Whitbourn, L., Huntington, J., Mason, P., Green, A., Berman, M., Coward, D., Connor, P., Wright, W.,
Jolivet, M., Martinez, R., 2016. HyLogger-3, a visible to shortwave and thermal infrared reflectance
spectrometer system for drill core logging: functional description. Australian Journal of Earth Sciences 63, 929–
940. doi:10.1080/08120099.2016.1231133

19
Specim, 2013. SisuROCK [WWW Document]. URL http://www.specim.fi/products/sisurock-hyperspectral-core-
imaging-station/ (accessed 8.10.17).

Tiu, G., 2017. Master thesis : Classification of Drill Core Textures for Process Simulation in Geometallurgy: Aitik Mine,
New Boliden. Luleå University of Technology, Luleå, Sweden.

Tolosana-Delgado, R., Mueller, U., van den Boogaart, K., Ward, C., Gutzmer, J., 2015. Improving processing by
adaption to conditional geostatistical simulation of block compositions. Journal of the Southern African
Institute of Mining and Metallurgy 115, 13–26.

Tungpalan, K., Wightman, E., Manlapig, E., 2015. Relating mineralogical and textural characteristics to flotation
behaviour. Minerals Engineering 82, 136–140. doi:10.1016/j.mineng.2015.02.005

van den Boogaart, K., Konsulke, S., Delgado, R.T., 2013. Non-Linear Geostatistics for Geometallurgical Optimisation.
Proceedings of the Second AusIMM International Geometallurgy Conference 253–257.

Vos, C.F., Stange, W., Bradshaw, D.J., 2014. A new small-scale test to determine flotation performance – Part 1:
Overall performance. Minerals Engineering 66–68, 62–67. doi:10.1016/j.mineng.2014.04.015

Walters, S.G., 2011. Initiatives to Support Geometallurgical Mapping and Modelling, in: Proceedings of the 1st
AusIMM International Geometallurgy Conference. Brisbane, pp. 273–278.

Walters, S.G., 2009. New Research Initiatives in Geometallurgical Integration – Moving Towards a Common
Operating Language, in: 7th International Mining Geology Conference. pp. 19–22.

Wanhainen, C., Billström, K., Martinsson, O., 2006. Age, petrology and geochemistry of the porphyritic Aitik
intrusion, and its relation to the disseminated Aitik Cu-Au-Ag deposit, northern Sweden. GFF 128, 273–286.
doi:10.1080/11035890601284273

Wanhainen, C., Broman, C., Martinsson, O., Magnor, B., 2012. Modification of a Palaeoproterozoic porphyry-like
system: Integration of structural, geochemical, petrographic, and fluid inclusion data from the Aitik Cu–Au–Ag
deposit, northern Sweden. Ore Geology Reviews 48, 306–331. doi:10.1016/j.oregeorev.2012.05.002

White, B.W., Rosenblatt, F., 2006. Principles of Neurodynamics: Perceptrons and the Theory of Brain Mechanisms.
The American Journal of Psychology 76, 705. doi:10.2307/1419730

20
Paper II
 Textures in geometallurgy: from characterization to applications

Pierre-Henri Koch*,a, Glacialle Tiub, Parisa Semsaria, Mehdi Pariana, Viktor Lishchuka,
Yousef Ghorbania, Jan Rosenkranza

* Corresponding author: [email protected]

a Luleå University of Technology, Division of Minerals and Metallurgical Engineering

b Luleå University of Technology, Division of Geosciences and Environmental Engineering

SE-97187 Luleå, Sweden

Abstract
Recent studies in the field of geometallurgy have identified challenges compared to the quantification of
texture and its applications. For practitioners, it is relevant to know whether textures play a role in the
problem at hand, which parameters can be measured and their interpretation in terms of minerals and
particles. This study offers new perspectives, relevant tools and recent industrial case studies regarding
textures in geometallurgy. Clarifications regarding the hypotheses and models one takes when acquiring
a mineral map both in two and three dimensions are also presented. Contributions from different fields
such as spatial statistics, image processing but also machine learning are used to rationalize the workflow
from ore sampling to predictive modelling. Textural information is available in many cases and contains
relevant data about the deposit itself but also the behavior of particles in comminution and concentration
circuits. As such, extracting its value is a step towards a more efficient production in the mining and
mineral processing industry.

1 Introduction and motivation

According to Brodie, Fettes, Harte, & Schmid, (2007) and references within, structure, fabric and texture
are the terms to describe geometrical properties of rock elements. According to Wills and Finch (2015),
the texture refers to the limited number of rock properties such as size, dissemination, association, and
shape of the minerals within the ore. The current state-of-the-art techniques for identifying, classifying,
and processing textural information are done by processing digital images with image processing software
(Berry et al., 2008; Dutta et al., 2014). The images are often available as classified mineral maps produced
with quantitative/automated scanning electron microscopy (Gottlieb et al., 2000; Gy, 2003), light imaging
(Pirard et al., 2007) and hyperspectral imaging (Burley et al., 2017; Fox et al., 2017; Ramanaidou and Wells,
2011; Schodlok et al., 2016). As an addition to images, non-geometric information (see Table 1) can be
included which suggests an extended definition of textures. Depending on the additional information it
has been referred to as geochemistry-mineralogy-texture (Parbhakar-Fox et al., 2013), geometallurgical
texture (Koch, 2017), or GEM-type (Lund, 2013). Textures aggregated with process information and
expressed as a simple numerical parameter could then be modelled by geostatistics to estimate the
metallurgical response on a block-by-block basis in an ore block model (Williams, 2013). In the context of
a geometallurgical model, an approach based on the simulation of orebodies is provided by Lishchuk,
Lamberg, & Lund (2016)

1
 Table 1 Some properties which can be aggregated with textural information

Properties Sub-category Data collected

Geological Chemical List of chemical elements

Mineralogical List of minerals
Physical e.g., hardness, density
Processing Breakages Energy, PSD, liberation
Concentrations Recovery, mass pull
Qualities Grades, penalties
Environmental Acid rock drainage
Metallurgical Qualities Recoveries, mass pull
Environmental CO2 emission
Reflective Visual e.g., color, luster, fluorescence
Multispectral Mineral phases
Spatial Location of origin Coordinates x, y, z
Mechanical Rock quality designation (RQD) RQD
Geophysical Conductivity Conductivity/Resistivity
Susceptibility Susceptibility

The motivation of this paper is to provide a structured introduction to the quantitative study of textures
for the purpose of applications in geometallurgy and mineral processing. Tools to measure, describe,
characterize, simulate and link textures to the processes will be introduced. While it aims to be as general
as possible, it represents an entry point rather than an exhaustive compilation of all work ever done in
this field. Recent references have been preferred over only historical ones, as they offer practical utilities
such as current software names and case studies. The qualitative approach to textures, while useful in
geology and mineralogy (Taylor, 2010), is outside the scope of this study.
The main objective is to present a new concept of texture that has emerged in recent years from the
geometallurgy community. Geometallurgy bridges the gap between different mining disciplines (geology,
mining engineering, mineral processing, environmental and finance). Though there are several
geometallurgical approaches, recent studies gear towards a data-centered approach wherein data
collected are directly transformed to machine-understandable information that allows direct analysis and
modelling (Koch et al., 2015; Koch, 2017; Lishchuk et al., 2018). Thus, it is important to have a quantifiable

2
data, such as textures and mineralogy, which can easily be translated and correlated to parameters in the
mine and mill operation.
In this study, the first key is to clarify how the concept of texture was defined by different authors and
propose an operative definition along with innovative tools to characterize ore textures using quantifiable
parameters that can be used to build texture models.

2 Background
The etymology of the word texture appears to be the Latin verb texĕre, which means to compose, to
interweave, and to construct. It is interesting to note that the progressive structuration of a medium can
relate to the genesis of an ore by physical and chemical processes. The idea of a texture as the result of
an iterative process is a basis for texture simulation. As a historical reference, early descriptions of rock
textures (specifically color, homogeneity and hardness) can be found in Theophrastus’s Περὶ λίθων (Eyles,
1958).
In metallurgy and mineralogy, the term “texture” is defined as the preferred crystallographic orientation
of grains or minerals while the general structure of a rock in a wider sense is referred to as “fabric”.
However, given the recent advances in image processing and its contributions to the study of both intact
rock masses and particles, the restriction of the term texture to a preferred orientation of crystals seems
too restrictive. Moreover, the use of the term texture is nowadays widespread in mineral processing and
process metallurgy at various scales. Therefore, when used in a geometallurgical context, the ambiguity
is eliminated.
Scale has been described as one of the main properties of a texture (Koch, 2017; Tiu, 2017). In this study,
texture is defined as the distribution of minerals and pores in a defined volume. This definition relies on
the scale of observation. It is directly linked to the smallest homogeneous volume considered. This volume
depends on the spatial physical measurements of the sample properties. Different scales have been
proposed in Bonnici (2012): macro-scale, where the minimum volume is of the order of 1 cubic meter ,
meso-scale, at a scale of 1 cubic centimeter (Bonnici et al., 2008; Lund et al., 2015; Pérez-Barnuevo et al.,
2016; Wightman et al., 2014) and micro-scale, with a minimum scale of 1 cubic micrometer (Pérez-
Barnuevo et al., 2013, 2012; Pirard et al., 2007).
The definition adopted in this work is a combination of geometrical features (size, shape and patterns)
and physical features (mineral phases being present or absent). As a result, attention will be paid to the
notion of a texture model that considers as much information as possible. Acquisition and representation
of texture data defines the first model often referred to as the mineral map.
The first step in the study of a texture is to acquire the necessary information to construct a mineral map.
Many options exist to acquire mineralogical information depending on the spatial resolution needed, the
accuracy, size of the sample and costs.

2.1 Mineral map acquisition methods

In acquisition methods, two main factors contribute to the obtained map: the spatial resolution which
defines the scale at which we study the map and the energy range of the observed signal. The resolution

3
will be defined as the length of the pixel or volume of the voxel in the case of 3D data. The energy will
provide information about the physics of the interaction between the imaging equipment and the sample.
In general, the wavelength measured by the device is related to the scale of interest: too short
wavelengths are associated with small objects, as for example in the case of copper that has an atomic
radius of the order of 140 pm which can produce a Kα signal at 8 kEV that can be measured in X-ray
fluorescence. X-ray tomography as a non-destructive analysis technique is highly applicable for the 3D
imaging and analysis of geological materials in the matter of quantification and qualification. This scanning
method has been used for identification, measurement and processing of texture and mineralogy (Becker
et al., 2016; Denison et al., 1997; Godel et al., 2013; Tsafnat et al., 2009). X-ray tomography is based on
density and linear attenuation values and, therefore, it is applicable to identifying the minerals that have
a large density contrast. A summary of the different instruments for texture characterization is given in
Table 2.
Table 2 : A list of common instruments for acquisition of textures

Method Instrument Application for References Dimension Notes

texture
characterization
XRF X-ray Chemical assays (Barnes et al., 2D/3D 1D spectrum
Fluorescence: (elements) 2018; Fisher
mapping or gun et al., 2015;
Togami et al.,
2000)
XRD X-ray Crystalline (Wicks et al., 2D 1D spectrum
diffractometer: structure 1995) contains
crystals or powder information
about the 3D
structure
Reflectance Optical Minerals, grain (Berry et al., 2D 1D reflectance
microscope size and 2008) spectrum at
equipped with orientation different
wavelengths
Cameras (visible
light or
multispectral)
Electron Scanning electron Chemical/Mineral, (Ginibre et al., 2D/1D
beam microscope (SEM), grain size and 2002)
transmission orientation (Krekeler and
electron Guggenheim,
microscope (TEM) 2008)

4
 UV Hyper/multi- Mineral (Cloutis et al., 2D Indicated for
spectral cameras 2008) oxides,
hydroxides
and clay
minerals due
to the metal-
O charge
transfer
absorption
IR Hyper/multi- Minerals, (Lyon and 2D 1D spectrum
spectral cameras Molecules, water Burns, 1963)
content
Terahertz Terahertz sensors Mineral mapping, (Hu and Nuss, 2D 1D spectrum
porosity, water 1995; Yu et
content al., 2010)
X-ray X-ray computed Density (Jardine et al., 3D Scalar field of
tomography micro- 2018) attenuation
tomography
(µXCT)
X-rays Mineral and (Artioli et al., Combination
Tomography crystalline 2010) of XRT and
3D
and structure XRD
Diffraction

2.2 Mineral map as a model

Once a characteristic signal has been measured and processed (i.e. noise removal and filters), each signal
can be assigned to a location on the sample. If the measurement yields direct information about minerals,
a mineral map is readily obtained. In other cases, the comparison of the data with databases offers a
correlation between the signal and the mineral phases so that after matching the data to the
corresponding minerals present at the surface, a mineral map can also be obtained. In the case of optical
microscopy or SEM measurements, greyscale images are used, and a mineral map can be produced by
thresholding the pixel values of the image (Partio et al., 2002; Pirard et al., 2007). Recently, the method
of grey level co-occurrence matrix (GLCM) has been applied to 3D X-ray computed tomography of drill
cores (Becker et al., 2016; Jardine et al., 2018).
Limitations regarding mineral maps exist. Ideally, the signal should have the same dimension as the
mineral map. This is not the case when a 2D mineral map is to be matched with 3D measurements for
example. Studies done by Reyes (2017) and Tiu (2017) compared the data extracted from mineral maps
generated from 2D and 3D data (Figure 1). Both studies showed significant stereological effect.

5
Corrections have been developed in stereology but if the object to be studied is in 3D, additional methods
should be applied to decrease uncertainty (e.g. X-ray tomography combined with SEM). The quality of the
data heavily depends on the method. In general, an error associated with the measurement can be
obtained, both in terms of mineral identification as well as spatial resolution.

Figure 1 : Comparison of 2D and 3D-generated mineral maps from Tiu (2017). (A) SEM-BSE image and its (B) equivalent µXCT
image slice overlain by generated mineral map with the extracted modal mineralogy shown in C. Results show a comparative
result between the SEM and 2D µXCT equivalent slice but a large underestimation of 2D analysis compared to 3D, indicating
stereological effect.

Computers can store mineral maps in various ways, typically by coding each mineral phase by an integer
and representing the data on a Euclidean 2D or 3D grid in which the smallest element is a pixel or a voxel.
Alternatively, a graph with nodes and edges can be used. Mineral maps can represent an intact rock
sample or a particulate sample. In the case of particles, attention should be paid to touching particles and
to the representativeness of the field that is measured. Sampling particulate systems has been studied by
(Gy, 1979) and a useful toolbox for the practitioner can be found in Petersen, Minkkinen, & Esbensen
(2005).

6
3 Characterization of textures
The main object to characterize textures is a 2D or 3D field for which a measurement is available at each
pixel or voxel location. Both particles and intact rock sample mineral maps can be studied using the same
tools. A preliminary step in the case of particles would be to properly code the background as a separate
phase and keep in mind that the area that includes both background and mineral phases, is part of the
surface properties of the particles.
Two major points of view can be taken when quantifying the textural properties of a mineral map. The
first one is to only carry out measurements to discriminate between textures. In this case, no model is
required, and no assumption is made regarding the statistical distribution of the measured features. Since
this approach is simplistic, it presents several limitations. The first risk is to make underlying assumptions
without explicitly mentioning them, e.g. when measuring the mean grain size of a given texture based on
a limited number of epoxy-mounted samples. In this case, one assumption is that the mean particle size
does not change when the size of the field under study is extended. A second assumption is that for the
same texture class, the value of the mean does not change excessively by analyzing other samples of the
same class. In other words, statistics measured over a mineral map need to have certain properties to
have a physical meaning. An additional risk is to disconnect the extraction of features from a mineral map
and their implications in terms of genesis or possible behavior in a process. In that case, the risk of
selecting meaningless features for minerals engineering applications increases. These two concepts have
been discussed, both in the image processing community (Serra, 1982) and in the geostatistics community
(Chiles and Delfiner, 2009). It is interesting to note that mathematical morphology and geostatistics had
similar early concerns regarding variability and scales (Matheron and Serra, 2002).
The second point of view is to consider any spatial measurement as a random field and choose which
hypotheses are relevant for the case at hand. As such, the field is already a model of a texture.

3.1 Stationarity
In the case of phase elements in space, strict stationarity means that the same physical process is repeated
in space at a given distance (or lag). A less restrictive class of models integrates second-order stationarity,
meaning that the increments of the random function are stationary. Figure 2 shows an example of a
stationary and non-stationary synthetic mineral map in which each mineral phase is displayed as a color.

7
Figure 2 : left: an example of a texture that can be modeled with a stationary model. right: texture that should not be modeled
using a stationary model.

3.2 Ergodicity
Ergodicity is needed whenever one wishes to estimate statistics of a random function based on only one
example. For example, in geostatistics, the observed mineral phases are a unique case and cannot be
replicated. Assuming ergodicity means that one can infer properties of the underlying random function
from which only one realization is observed. In the study of statistics inside the same ore type, assuming
ergodicity means that one assumes that all the samples from the same ore type derive from the same
random function. It is the responsibility of the user to choose an appropriate model depending on the
application. If the estimation of a mean is needed, then a stationary model may be required. Ergodicity
can be used in the case where a mean is used, for example to distinguish different textures. In that case,
ergodicity ensures that different physical phenomena are not pooled in a same texture class.

3.3 Global properties

To evaluate meaningful global properties of a mineral map, one must make assumptions regarding the
underlying model.
A practical tool to evaluate the validity of the model assumptions was provided by Lantuéjoul (1991). The
idea is to relate the variance measured in spatial subdomains to ergodicity and stationarity. In practice,
the mineral map is split into smaller subdomains (patches) and the dispersion variance is evaluated as a
function of the size of the subdomains. The slope in logarithmic scale provides an empirical law of change
of scales. Figure 3 shows four examples of simulated fields derived from models with different scale
properties. One can clearly see that in the case of the red points, the dispersion variance varies linearly
and has a slope lower than 1. The blue points suggest a similar case with a slightly worse situation when
the size of the patches is small. In these two cases, one can safely estimate the mean of a phase under the
assumption that the mineral map is one realization of a stationary, ergodic random function. One can
further choose the correct patch size for estimation purposes at a given variance. The two other cases in

8
green and yellow represent the case where either the maximum size of the patches is not large enough
or that the underlying stationary model either has an infinite integral range or is not stationary. The case
of an infinite integral range will result in a higher variance, but the other case could mean that a stationary
model is not appropriate (in that case, a spatial mean for the entire field is not meaningful). All the tools
and discussion introduced in this work will focus on stationary cases.

Figure 3 : Dispersion variance as a function of the size of subdomains. Solid lines are experimental data points and dashed
lines present the empirical law of change of scale: c/K where c is the origin and K the size of the patches. Axes are
logarithmic.

3.4 Geometry
Usually, the dimension and shape measurements apply to particles and not to intact textures. But in
certain cases, when an ore texture presents a grain-like structure, the study of dimension and shape can
be done on mineral grains rather than on separate particles (for instance by thresholding a mineral map
with N phases to a binary mineral map containing the phase of interest and the rest). Figure 4 shows a
drill core from the Aitik mine, a porphyry copper deposit in Northern Sweden. The first image is the high-
resolution visible light image, the second one is a mineral map obtained by segmentation and the bottom
image is the mineral map after binary masking chalcopyrite only.

9
1 cm

Figure 4 : Example of thresholding and binary mask selective to chalcopyrite

Among the different approaches to characterize size and shape of particles, the development of
mathematical morphology by Matheron and Serra in the 1970 was important. It offered a range of tools
describing geometrical features on binary sets by applying morphological filters (like opening and closure).

3.4.1 Size
Direct measurements of size based on mineral maps are based on a sphere of equivalent diameter. It
relies on an appropriate distance function defined on the mineral map. The simplest descriptors are the
minimum (dmin), maximum (dmax) and mean diameters of a particle or grain. Additionally, other size
descriptors have been developed for example

• Feret’s diameter: measured as the distance between two parallels planes tangent to the particle
at opposite points. A global value can be calculated by averaging the distance over different pairs
of opposite points
• Martin’s diameter: measured as the length of a segment that separates the particle in two equal
areas
Feret and Martin diameters depend on the orientation of the measurement. To obtain a representative
measurement one should measure many particles (see 3.3).
An indirect measurement of shape relies on the phase specific surface area (PSSA). Using the mean
intercept length L, it is defined as PSSA = 4/L. Assuming the particles are cubes or spheres, the mean size
is given by DPSSA=6/PSSA.

10
3.4.2 Shape
A structured review of general shape descriptors can be found in Zhang and Lu (2004). For mineral maps,
an extensive study of shape is provided by Pirard (1993) in which the author describes desirable properties
for shape descriptors
• Independence with respect to other descriptors: avoid redundant measurements
• Sensitivity: descriptors should discriminate between particles of different shapes
• Specificity: each shape descriptor should capture one precise physical property
• Additivity: the mean of descriptors should have a physical interpretation
• Invariance: descriptors should be invariant with respect to rotation and translation
The problem of finding descriptors that verify all these properties is not obvious. Moreover, regarding
invariance, one can argue that scale invariance is missing. As an example of these properties, the shape
information recovered based on the number of Fourier descriptors is presented in Table 3 . The impact of
a slight deformation on the Fourier shape descriptors is illustrated in Table 4. The figures have been
produced in MATLAB based on an existing code (Peyré, 2011).
Table 3 : Recovery of shape information as a function of the number of Fourier descriptors

Original particle Fourier reconstruction as a Number of

function of the number of Fourier
descriptors descriptors
1

10

50

200

11
 Table 4 : Impact of slight deformations on the amplitude of the Fourier descriptors

Contour of the particle Amplitude of the first 20 Fourier Description

shape descriptors
Original
particle

Slightly
deformed
particle

Lower
resolution
of the
contour

In 3D, a natural extension of Fourier shape descriptors are spherical harmonics (Garboczi and Bullard,
2017; Mollon and Zhao, 2014; Su and Yan, 2018).

3.4.3 Association
In the context of mineral processing and comminution, the association between different mineral phases
or mineral phases and voids (in the case of cracks or the estimation of exposed surface) is a topic of
interest. By nature, association measurements rely on statistics measured over a mineral map. As such,
they are related to the field of spatial statistics in which the distribution of a random variable in space is
studied. The basis for association measurements is to compare two or more different points in space. An
example of two points comparison is the experimental semi-variogram γ(h) or V(h). To illustrate the
similarity between these methods, their formulas are given as Equation 1. They all have in common that
the squared difference is computed between two values at different locations.

12
 𝑁𝑁(ℎ)
𝟐𝟐
1 ′
𝑉𝑉(ℎ) = � �𝐺𝐺(𝒙𝒙) − 𝐺𝐺�𝒙𝒙 (ℎ)��
2𝑁𝑁(ℎ)
𝑖𝑖=0
𝑃𝑃−1
(1)
𝐿𝐿𝐿𝐿𝑃𝑃𝑃𝑃,𝑅𝑅 = � 𝑠𝑠�𝑔𝑔𝑝𝑝 − 𝑔𝑔0 �2𝑝𝑝
𝑝𝑝=0

𝑛𝑛

𝑈𝑈(𝒙𝒙) = 𝛽𝛽 � (𝑥𝑥𝑠𝑠 − 𝑥𝑥𝑡𝑡 )2 + 𝛼𝛼 �(𝑥𝑥𝑠𝑠 − 𝜇𝜇𝑠𝑠 )2

𝑐𝑐=𝑠𝑠,𝑡𝑡 𝑘𝑘=0

where
V(h) is the semi-variogram where h is the lag and x the location (Chiles and Delfiner, 2009)
LBP is the value computed in local binary patterns where g is the greyscale value of pixels (Ojala et al.,
2002)
and
U(x) is the potential function in Potts model used for segmentation or simulation of textures (Koch and
Rosenkranz, 2017).
The most commonly used two points association indicators are grey-level co-occurrence matrices (GLCM,
Becker et al., 2016; Jardine et al., 2018). It has been shown that GLCM and semi-variograms are related
(van der Sanden and Hoekman, 2005). Furthermore, under the right hypotheses, the Wiener-Khintchine
theorem (Khintchine, 1934; Wiener, 1930) links the covariance with the Fourier transform of the random
variable.
Variants of GLCM have been introduced, some of them using extended neighborhoods (Koch, 2017) or
focusing on the association at the boundary of phases (Parian et al., 2018).

4 Classification of textures
Once features have been derived from the texture representation, classification offers a natural way to
discriminate between different ore types. Fundamentally, a reduction of dimension is applied to the set
of features. In simple cases, it is a projection of a vector into a lower dimension space. The extracted
features can be used for automated textural classification using machine-learning methods.
In many cases, the classification is correlated with the behavior of the resulting particles in the
concentration process. It implies that some physical properties of the ore have been captured and helps
in discriminating different populations of particles produced by comminution.
A study done by Tiu (2017) employed an automated textural classification and segmentation of drill core
images through feature extraction and machine-learning method. The textural classes were also subject
to geometallurgical comminution tests (Bilal, 2017), where they showed significant differences in their
comminution behavior (Figure 5).

13
Figure 5 : Bond Work Index (BWI) and copper content of different textural classes from Tiu (2017). Examples of different textural
classes (drill core photos and mineral maps) are shown on the right.

5 Applications of mineral maps

5.1 Simulation of textures
Once a mineral map has been obtained and a representation has been chosen, different options exist to
simulate a texture that has similar properties. This topic is a research field on its own at the intersection
of image processing (Galerne et al., 2011), spatial statistics (Mariethoz and Lefebvre, 2014), compression
and geometry (Mollon and Zhao, 2014).The underlying idea is derived from spectral analysis, i.e. the
features measured should contain enough information to simulate a map that can reproduce them. These
methods are useful in process simulation when an initial estimate of liberation is needed to simulate an
ore beneficiation process at liberation level. Once the particles are obtained, advanced unit models and
particle tracking are required in order to fully use this information (Cárdenas, 2017; Hannula et al., 2018;
Lamberg and Vianna, 2007).

5.2 Comminution
The nature of the particles produced after blasting, crushing, and grinding have long been a topic of
interest in mineral liberation studies. The size, shape and mineral association in particles links
mineralogy, geology and geometallurgy since these properties are inherited from the intact ore to the
concentrate. As such, particles are a bridge between the mine and the process that close the loop of
geometallurgy (Lamberg, 2011).

14
Therefore, to complete the chain of geometallurgy, a study on the liberation distribution after
comminution from a texture point of view was conducted (Koch and Rosenkranz, 2017; Mariano, 2016;
Parian et al., 2018). Predicting particle populations after comminution started by constructing a mineral
map of ore texture. The particle size distribution after breakage was used for generating various sizes of
rectangles that were superimposed over the mineral map. For positioning the rectangles over the mineral
map, the breakage probability was applied. The assumption was that the preferential breakage in the
phase is the cause of the liberation of minerals and random breakage can yield both liberated and non-
liberated particles (Figure 6). Comparison of results from liberation prediction and experimental works
shows positive correlation (Figure7) confirming that texture is a significant factor on particle liberation
and consequently downstream processing (Parian et al., 2018).

100 100

90 90

80 80

70 70

60 60
Cumulative mass (wt%)
Cumulative mass (wt%)

50 50

40 40

30 30

1680-3350 (µm) 20 1680-3350 (µm)

20
850-1680 (µm) 850-1680 (µm)
425-850 (µm) 425-850 (µm)
10 10
300-425 (µm) 300-425 (µm)

0 0
0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100

Mgt (wt%) in particle Mgt (wt%) in particle

Figure 6 : Modeled liberation distribution of magnetite in four size fractions for random breakage (left) and 100% preferential
breakage (right).

100 100

90 90

80 80

70 70

60 60
Cumulative yield (%)
Cumulative yield (%)

50 50

40 40

30 30

1680-3350 (µm) 20 1680-3350 (µm)

20
850-1680 (µm) 850-1680 (µm)
425-850 (µm) 425-850 (µm)
10 10
300-425 (µm) 300-425 (µm)

0 0
0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100

Magnetite (wt%) liberation Magnetite (wt%) liberation

Figure 7 : Liberation distribution of magnetite particles from experimental (left) and forecasting particles from mineral map
(right) in four size fractions.

15
The particle heterogeneity (i.e. texture at the micro-scale) affects and controls the breakage properties of
a given material. Comparison between the breakage behavior of various textures for iron ores revealed
that grain size and mineral composition can empirically approximate the specific rate of breakage
(Mwanga et al., 2015). Experiments showed that materials with similar modal composition (magnetite
grade) but different texture have different specific rates of breakage (Figure 8). Observations further
showed that the breakage rate of particles decreases as the particle size reaches the grain size of a mineral
(Figure 9). The grade by particle size pattern gives an estimate on liberation size and can be used to choose
the right size class that enables liberation of minerals and efficient comminution (Mwanga et al., 2015).
0.4

0.35 Experimental FAR

Specific rate of breakage (per minute)

Modeled FAR
0.3 Experimental HAR

Modeled HAR
0.25 Experimental Mgt-ore

Modeled Mgt-ore
0.2 Experimental Hmt-ore

Modeled Hmt-ore
0.15 Experimental KA3

Modeled KA3
0.1 Experimental F

Modeled F

0.05

0
100 1000
Particle size (µm)

Figure 8 : Comparison of the experimental and modeled specific rate of breakage by using Austin model for texturally different
iron ore samples (Mwanga et al., 2015)

16
 B1

Size at which breakage of FeO grade start to decrease (µm)

1200

FAR

HAR

y = 205.85e0.002x
R² = 0.9821
KA3
D3

120
0 200 400 600 800 1000 1200

Grain size of the dominating i.e. FeO (µm)

Figure 9 : Correlation between grain sizes of the dominating mineral vs. size at which the breakage efficiency starts to decrease
for five texturally different iron ore samples (Mwanga et al., 2015)

5.3 Separation
One of the case studies where textures have been directly related to the behavior in ore concentration
includes the Leveäniemi mine, a large apatite-iron ore deposit in the Norrbotten region, Sweden.
Leveäniemi is operated as an open pit mine after it was re-opened in 2015 (Gustafsson, 2016; Niiranen,
2015). Several ore types can be observed in this area that differ massive magnetite, for example
calcareous magnetite ore, hematite-altered ore, and brecciated ore (Frietsch, 1966). Lishchuk, Lund, Koch,
Gustafsson, & Pålsson (2019) conducted a geometallurgical study for the Leveäniemi deposit, which
resulted in predictive process models and an ore classification based on mineral and textural variations
(Table 5).

17
 Table 5 A classification of Leveäniemi ore based on mineral and textural variations (Lishchuk, 2019).

Iron wt%
Fe- Main associated
ID Ore classes in drill Textural type
mineral minerals
core
13s <20 Disseminated/Veiny Mgt+Bt+Ab+Qtz
1s Semi-massive Veiny/Granular Mgt+Bt+Kfsp+Qtz
Mgt
6s ore 20-40 Patchy Mgt+Bt+Amph+Qtz+Ap
9s Banded Mgt+Bt+Kfsp+Qtz
Rich semi-
12s Mgt 40-50 Coarse-grained amphibole rich Mgt+Ap+Ab+Bt+Qtz
massive ore
Fine-grained amphibole/Calcite
5m Mgt+Amph+Ap+Cal
rich
Coarse-grained
8m Mgt+Amph+Bt+Qtz+Cal
amphibole/Calcite rich
10m Coarse-grained amphibole rich Mgt+Amph+Ap
Massive ore Mgt
11m >50 Fine-grained amphibole rich Mgt+Amph+Ap
2m Mgt+Amph+Ap
4m Coarse-grained apatite rich Mgt+Amph+Ap
7m Mgt+(Amph)+Ap
3m Hem Fine-grained hematite Mgt+Ap
Abbreviations: s -semi-massive; m-massive; Mgt-Magnetite; Hem-Hematite; Amph-Amphibole; Cal-
Calcite; Bt-Biotite; Ab-Albite Qtz-Quartz; K-Fsp-K-Feldspar.

18
The connection between geological and processing properties of the ore was done by developing a new
ore estimator for classifying magnetite ore denoted by 𝑿𝑿𝑳𝑳𝑳𝑳𝑳𝑳. The 𝑋𝑋𝐿𝐿𝐿𝐿𝐿𝐿 classifies the ore by considering
both its composition and the behavior in magnetic separation tests done with a Davis tube (DT) (Davis,
1923; Lishchuk et al., 2018; Niiranen, 2015). 𝑋𝑋𝐿𝐿𝐿𝐿𝐿𝐿 equals zero when the available magnetite in the ore
reports to the Davis Tube concentrate as pure magnetite (Lishchuk et al., 2018), it can be calculated
from the iron grade of the feed and the mass pull :
𝑋𝑋𝐿𝐿𝐿𝐿𝐿𝐿 = 𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹 − 𝑀𝑀𝑀𝑀 ∙ 0.724 (2)
Positive values of 𝑋𝑋𝐿𝐿𝐿𝐿𝐿𝐿 indicates a lack of iron in the Davis Tube concentrate. This means that a part of
the iron is lost to the tailings. Negative values of 𝑋𝑋𝐿𝐿𝐿𝐿𝐿𝐿 occur when the mass pull (MP) is larger than
expected. In this case, it means that the concentrate is not as clean as expected (Lishchuk et al., 2018).
The textural difference between the ore types can be depicted by plotting 𝑿𝑿𝑳𝑳𝑳𝑳𝑳𝑳 against liberation of Fe-
oxides, see Figure 10, which shows the massive and semi-massive ore samples from two separate clusters.
In the Leveäniemi case, geological textural types clearly overlap with geometallurgical classes defined for
the ore response in the Davis tube test.

Figure 10.Ore classification with 𝑋𝑋𝐿𝐿𝐿𝐿𝐿𝐿 (Lishchuk, 2019).

As a result, such a classification could be combined with the classification developed by Lund (2013),
providing geological textural types and geometallurgical (GEM) classes as shown on Table 6.

19
 Table 6. Description of the Leveäniemi ore samples according to the geometallurgical ore types (GEM-type) (Lund, 2013).

Ore type Ore mineral GEM type GEM sub-type GEM code Leveäniemi
samples
Massive Mgt Amph+Ap+Bt+Py Amph 2b 2, 8
Hem Amph 2b 3
Mgt Ap 2c 4, 5, 7, 10, 11
Rich semi-massive Mgt Qtz - 3 12
Semi-massive Mgt Amph+Ap+Bt+Py Bt 2d 1, 6
Mgt Qtz - 3 9, 13

Similar results have been obtained for the characterization of two Li-pegmatite deposits being part of
Keliber’s lithium project, located in Kaustinen, Ostrobothnia, Finland (Guiral Vega, 2018). Based on this
work, a methodology in four main stages has been proposed: (i) Identification of mineral features at
different scales: This stage includes the mineral characterization both in macro scale, through drill core
login, and in micro scale through the application of micro analytical techniques such as SEM-EDS (see
Table 7), (ii) Mineral quantification from drill cores: a model, based on image analysis, was created for the
automated identification and subsequent quantification of the spodumene contained in the drill cores,
(iii) Spatial modelling: a spatial model was proposed to indicate the idealised distribution of each textural
class within the ore body and (iv) Texture modelling: the relationship between textural classes and
processing behaviour has been proposed through the implementation of a machine learning model which
uses chemical composition of drill cores to forecast the Li2O flotation recovery classified by textural class
(Figure 11).

20
 Table 7. Textural classification scheme for Li-pegmatite deposits at Keliber Lithium Project (Guiral Vega, 2018).

Spd
Texture Typical Association Spd crystal grain Texture
Zone Ore type Class
type lithology minerals distribution size Sub-type
[cm]
Spd – Ab – Coarse
>4 C1
Qtz – Kfs – grain
Spodumene Ms ± Bt ± Ap Preferentially Medium
Foliated 2–4 C2
pegmatite ± Py ± Apy ± oriented grain
Sp ± Tur ± Fine
<2 C3
Spd Grt grain
pegmatite Spd – Ab – Coarse
Randomly >4 C4
Qtz – Kfs – grain
distributed
Spodumene Ms ± Bt ± Ap Medium
Isotropic (Not 2-4 C5
pegmatite ± Py ± Apy ± grain
preferentially
Sp ± Tur ± Fine
oriented) <2 C6
Grt grain
Muscovite
Qtz – Ms – Disseminated
granite / 0.1 –
Pegmatite Ab ± Spd ± altered Banded C7
Muscovite 0.7
Py crystals
pegmatite
Foliated Ab – Qtz –
Quartz- Ms ± Spd ± Disseminated
0.1 –
feldspar Bt ± Hbl ± altered Veiny C8
0.7
Barren veins Tur ± Grt ± crystals
pegmatite Ap
Muscovite Ms – Ab –
Disseminated
granite / Qtz ± Spd Bt 0.1 –
altered Granular C9
muscovite ± Py ± Apy ± 0.7
crystals
Isotropic pegmatite Sp
Quartz- Qtz – Ab ± Disseminated
0.1 –
feldspar Spd ± Kfs ± altered Massive C10
0.7
pegmatite Ms crystals
Pl – Hbl –
Andesite Qtz ± Bt ±
n.a. n.a. n.a. C11
porphyry Chl ± Py ±
Host rock n.a. Foliated Apy ± Sp
Qtz – Pl – Ms
Mica schist ± Bt ± Apy ± n.a. n.a. n.a. C12
Sp

21
 Figure 11. Forecasting of Li2O flotation recovery per textural class - Machine learning model. Error bars denote standard
deviation (i.e. ± std) (Guiral Vega, 2018).

This methodology considers texture as the main parameter for linking the ore features with the final
flotation recovery. Thus, it is capable to take information from drill cores and transform it into
downstream processing performance data. Further works using the same methodology will be applied in
the future on the mineral and textural characterization of different iron ores from Canada, to understand
the influence of mineral textures on the efficiency of magnetic separation processes.

6 Conclusions and recommendations

Through a brief review and several practical examples where textural information was effectively used to
enhance knowledge and prediction capability in geometallurgy contexts, the value that lies in textural
information has been highlighted. This refers to both texture information on intact rock material or after
crushing and grinding. Indeed, today, most of the information is available (for example using existing
instruments such as optical microscopes or SEM) and with new technologies the cost is constantly
decreasing. The appropriate numerical tools are widely available, both as commercial and open-source
software. Not measuring textures could result in a missed opportunity for better understanding of the ore
deposits.
An expertise in image processing is not required for most practical cases since available algorithms exist
to verify whether, for instance, stationarity is a correct hypothesis. Proper representative sampling can
be prepared, based either on Gy’s theory of sampling for particulate materials or based on a geological
model. The interpretation of the results and their validity with respect to the hypotheses, however, do
require expertise.

22
The main textural features are size, shape and association, for which relevant methods have been
presented. The data produced by each method can be used for visualization and classification purposes
but can also be directly used in process simulation. In several cases, textural features can control
comminution or concentration stages in a mineral processing plant. One of the key aspects is the
liberation spectrum of the particles, which influences both, the grade and recovery of the minerals of
interest.
Recent methods that combine the extraction of shape, size and association at different scales are for
example convolutional neural networks. While the interpretation of which parameter is learned and
where it is stored in such models is not trivial, they seem to offer the advantage of taking all the
parameters into account at different scales. The action of these networks is to linearize small
deformations by applying a sequence of multiscale convolutions similar to wavelet transforms and non-
linear activation functions (Mallat, 2016). Again, these algorithms are available as open-source packages
that allow both to classify textures but also combine textural information with operating parameters in a
geometallurgical program.

7 Acknowledgements
The authors wish to thank LKAB, Boliden Mineral AB and Outotec (Sweden) AB for their support in this
study. This work is part of PREP project funded by VINNOVA within the Strategic innovation program for
the Swedish mining and metal producing industry (SIP-STRIM).

8 References
Artioli, G., Cerulli, T., Cruciani, G., Dalconi, M.C., Ferrari, G., Parisatto, M., Rack, A., Tucoulou, R., 2010. X-
ray diffraction microtomography (XRD-CT), a novel tool for non-invasive mapping of phase
development in cement materials. Analytical and Bioanalytical Chemistry 397, 2131–2136.
doi:10.1007/s00216-010-3649-0
Barnes, S.J., Piña, R., Le Vaillant, M., 2018. Textural development in sulfide-matrix ore breccias in the
Aguablanca Ni-Cu deposit, Spain, revealed by X-ray fluorescence microscopy. Ore Geology Reviews
95, 849–862. doi:10.1016/j.oregeorev.2018.03.004
Becker, M., Jardine, M.A., Miller, J.A., Harris, M., 2016. X-ray Computed Tomography: A geometallurgical
tool for 3D textural analysis of drill core?, in: Proceedings of the 3rd AusIMM International
Geometallurgy Conference. pp. 15–16.
Berry, R., Walters, S.G., McMahon, C., 2008. Automated mineral identification by optical microscopy, in:
Ninth International Congress for Applied Mineralogy (ICAM). The Australasian Institute of Mining
and Metallurgy: Melbourne, Melbourne, Australia, pp. 91–94.
Bilal, D., 2017. Master thesis: Geometallurgical estimation of comminution indices for porphyry copper
deposit applying mineralogical approach. Luleå University of Technology.
Bonnici, N., Hunt, J.A., Walters, S.G., Berry, R., Collett, D., 2008. Relating textural attributes to mineral
processing: Developing a more effective approach for the Cadia East Cu-Au porphyry deposit., in:
Proceedings of the Ninth International Congress for Applied Mineralogy Conference (ICAM). pp. 4–
5.

23
Bonnici, N.K., 2012. The mineralogical and textural characteristics of copper-gold deposits related to
mineral processing attributes. University of Tasmania.
Brodie, K., Fettes, D., Harte, B., Schmid, R., 2007. Structural terms including fault rock terms. Khb 1–14.
doi:10.1007/978-1-4939-7216-6
Burley, L.L., Barnes, S.J., Laukamp, C., Mole, D.R., Le Vaillant, M., Fiorentini, M.L., 2017. Rapid
mineralogical and geochemical characterisation of the Fisher East nickel sulphide prospects,
Western Australia, using hyperspectral and pXRF data. Ore Geology Reviews 90, 371–387.
doi:10.1016/j.oregeorev.2017.04.032
Cárdenas, E., 2017. Master thesis: Particle tracking in geometallurgical testing for Leveäniemi Iron ore ,
Sweden. Luleå University of Technology.
Chiles, J.-P., Delfiner, P., 2009. Geostatistics: modeling spatial uncertainty. John Wiley & Sons.
Cloutis, E.A., McCormack, K.A., Bell, J.F., Hendrix, A.R., Bailey, D.T., Craig, M.A., Mertzman, S.A., Robinson,
M.S., Riner, M.A., 2008. Ultraviolet spectral reflectance properties of common planetary minerals.
Icarus 197, 321–347. doi:https://doi.org/10.1016/j.icarus.2008.04.018
Davis, E.W., 1923. November 20. Magnetic concentrator. U.S. Patent 1,474,624.
Denison, C., Carlson, W.D., Ketcham, R.A., 1997. Three-dimensional quantitative textural analysis of
metamorphic rocks using high-resolution computed X-ray tomography: Part I. Methods and
techniques. Journal of Metamorphic Geology 15, 29–44. doi:10.1111/j.1525-1314.1997.00006.x
Dutta, S., Barat, K., Das, A., Das, S.K., Shukla, A.K., Roy, H., 2014. Characterization of micrographs and
fractographs of Custrengthened HSLA steel using image texture analysis. Measurement 47, 130–144.
Eyles, V.A., 1958. Theophrastus on Stones. Nature 181, 1686.
Fisher, L.A., Fougerouse, D., Cleverley, J.S., Ryan, C.G., Micklethwaite, S., Halfpenny, A., Hough, R.M., Gee,
M., Paterson, D., Howard, D.L., Spiers, K., 2015. Quantified, multi-scale X-ray fluorescence element
mapping using the Maia detector array: application to mineral deposit studies. Mineralium Deposita
50, 665–674. doi:10.1007/s00126-014-0562-z
Fox, N., Parbhakar-Fox, A., Moltzen, J., Feig, S., Goemann, K., Huntington, J., 2017. Applications of
hyperspectral mineralogy for geoenvironmental characterisation. Minerals Engineering 107, 63–77.
doi:10.1016/j.mineng.2016.11.008
Frietsch, R., 1966. Berggrund och malmer I Svappavaarafältet, Norra Sverige, SGU. Sveriges geologiska
undersökning (SGU). Serie C, Avhandlingar och uppsatser, 0082-0024; 604, Stockholm, Sweden.
Galerne, B., Gousseau, Y., Morel, J.M., 2011. Random phase textures: Theory and synthesis. IEEE
Transactions on Image Processing 20, 257–267. doi:10.1109/TIP.2010.2052822
Garboczi, E.J., Bullard, J.W., 2017. 3D analytical mathematical models of random star-shape particles via
a combination of X-ray computed microtomography and spherical harmonic analysis. Advanced
Powder Technology 28, 325–339. doi:https://doi.org/10.1016/j.apt.2016.10.014
Ginibre, C., Kronz, A., Wörner, G., 2002. High-resolution quantitative imaging of plagioclase composition
using accumulated backscattered electron images: New constraints on oscillatory zoning.
Contributions to Mineralogy and Petrology 142, 436–448. doi:10.1007/s004100100298
Godel, B., Barnes, S.J., Gürer, D., Austin, P., Fiorentini, M.L., 2013. Chromite in komatiites: 3D
morphologies with implications for crystallization mechanisms. Contributions to Mineralogy and
Petrology 165, 173–189. doi:10.1007/s00410-012-0804-y

24
Gottlieb, P., Wilkie, G., Sutherland, D., Ho-Tun, E., Suthers, S., Perera, K., Jenkins, B., Spencer, S., Butcher,
A., Rayner, J., 2000. Microtextural Mineralogy Testing & Analysis Using Quantitative Electron
Microscopy for Process Mineralogy Applications. The Journal of The Minerals, Metals & Materials
Society 52, 24–25.
Guiral Vega, J.S., 2018. Master thesis : Textural and Mineralogical Characterization of Li-pegmatite
Deposit, using Microanalytical and Image Analysis to Link Micro and Macro Properties of Spodumene
in Drill Cores. : Keliber Lithium Project, Finland. Luleå University of Technology.
Gustafsson, M., 2016. Master thesis: Element to Mineral Conversion of LKAB Deposits Element to Mineral
Conversion of LKAB deposits: Leveäniemi. Gruvberget Hematite and Sparre. Luleå University of
Technology.
Gy, P., 1979. Sampling of particulate materials : theory and practice, Developmen. ed. Elsevier Scientific
Publishing Company.
Gy, Y., 2003. Automated scanning electron microscope based mineral liberation analysis – An introduction
to the JKMRC/FEI mineral liberation analyser. Minerals Materials Characterization and Engineering
33–41.
Hannula, J., Kern, M., Luukkanen, S., Roine, A., van den Boogaart, K.G., Reuter, M.A., 2018. Property-based
modelling and simulation of mechanical separation processes using dynamic binning and neural
networks. Minerals Engineering 126, 52–63.
Hu, B.B., Nuss, M.C., 1995. Imaging with terahertz waves. Optics Letters 20, 1716.
doi:10.1364/OL.20.001716
Jardine, M.A., Miller, J.A., Becker, M., 2018. Coupled X-ray computed tomography and grey level co-
occurrence matrices as a method for quantification of mineralogy and texture in 3D. Computers and
Geosciences 111, 105–117. doi:10.1016/j.cageo.2017.11.005
Khintchine, A., 1934. Korrelationstheorie der stationären stochastischen Prozesse. Mathematische
Annalen 109, 604–615. doi:10.1007/BF01449156
Koch, P.-H., 2017. Licentiate thesis: Particle generation for geometallurgical process modeling. Luleå
University of Technology. Luleå University of Technology.
Koch, P.-H., Lamberg, P., Rosenkranz, J., 2015. How to build a process model in a geometallurgical
program?, in: Andre-Mayer, A., Cathelineau, M., Muchez, P., Pirard, E., Sindern, S. (Eds.), Proceedings
of the 13th SGA Biennial Meeting. Nancy, p. 85.
Koch, P.-H., Rosenkranz, J., 2017. Texture-based liberation models for comminution, in: Proceedings of
the Conference in Minerals Engineering. Luleå, Sweden, pp. 83–96.
Krekeler, M.P.S., Guggenheim, S., 2008. Defects in microstructure in palygorskite–sepiolite minerals: A
transmission electron microscopy (TEM) study. Applied Clay Science 39, 98–105.
doi:10.1016/J.CLAY.2007.05.001
Lamberg, P., 2011. Particles – the bridge between geology and metallurgy, in: Proceedings of the
Conference in Minerals Engineering (CME). Luleå, pp. 1–16.
Lamberg, P., Vianna, S., 2007. A technique for tracking multiphase mineral particles in flotation circuits,
in: Proceedings of the Meeting of the Southern Hemisphere on Mineral Technology (MSHMT).
Universidade Federal de Ouro Preto, Ouro Preto, Minas Gerais, Brazil, pp. 195–202.
Lantuéjoul, C., 1991. Ergodicity and integral range. Journal of Microscopy 161, 387–403.

25
 doi:10.1111/j.1365-2818.1991.tb03099.x
Lishchuk, V., 2019. Doctoral Thesis: Bringing predictability into a geometallurgical program: An iron ore
case study. Luleå University of Technology, Luleå, Sweden.
Lishchuk, V., Lamberg, P., Lund, C., 2016. Evaluation of sampling in geometallurgical programs through
synthetic deposit model, in: Proceedings of the XXVI International Mineral Processing Congress
(IMPC).
Lishchuk, V., Lund, C., Koch, P.-H., Gustafsson, M., Pålsson, B.I.I., 2018. Geometallurgical characterization
of Leveäniemi iron ore - Unlocking the patterns. Minerals Engineering 131, 325–335.
doi:10.1016/j.mineng.2018.11.034
Lund, C., 2013. Doctoral thesis: Mineralogical, chemical and textural characterisation of the Malmberget
iron ore deposit for a geometallurgical model. Luleå University of Technology Department.
Lund, C., Lamberg, P., Lindberg, T., 2015. Development of a geometallurgical framework to quantify
mineral textures for process prediction. Minerals Engineering 82, 61–77.
doi:10.1016/j.mineng.2015.04.004
Lyon, R.J.P., Burns, E.A., 1963. Analysis of rocks and minerals by reflected infrared radiation. Economic
Geology 58, 274–284. doi:10.2113/gsecongeo.58.2.274
Mallat, S., 2016. Understanding deep convolutional networks. Philosophical Transactions of the Royal
Society of London Series A 374, 20150203. doi:10.1098/rsta.2015.0203
Mariano, R., 2016. Measurement and modelling of the liberation and distribution of minerals in
comminuted ores. University of Queensland. doi:10.14264/uql.2016.1081
Mariethoz, G., Lefebvre, S., 2014. Bridges between multiple-point geostatistics and texture synthesis:
Review and guidelines for future research. Computers and Geosciences.
doi:10.1016/j.cageo.2014.01.001
Matheron, G., Serra, J., 2002. The birth of mathematical morphology, in: Proceedings of the 6th Intl. Symp.
Mathematical Morphology. Sydney, Australia, pp. 1–16.
Mollon, G., Zhao, J., 2014. 3D generation of realistic granular samples based on random fields theory and
Fourier shape descriptors. Computer Methods in Applied Mechanics and Engineering 279, 46–65.
doi:https://doi.org/10.1016/j.cma.2014.06.022
Mwanga, A., Rosenkranz, J., Lamberg, P., 2015. Testing of Ore Comminution Behavior in the
Geometallurgical Context—A Review. Minerals 5, 276–297. doi:10.3390/min5020276
Niiranen, K., 2015. Characterisation of the Kiirunavaara iron ore deposit for mineral processing with a
focus on the high silica ore type B2. Doctoral Thesis 133.
Ojala, T., Pietikainen, M., Maenpaa, T., 2002. Multiresolution gray-scale and rotation invariant texture
classification with local binary patterns. IEEE Transactions on Pattern Analysis and Machine
Intelligence 24, 971–987. doi:10.1109/TPAMI.2002.1017623
Parbhakar-Fox, A., Lottermoser, B., Bradshaw, D., 2013. Cost-effective means for identifying acid rock
drainage risks – integration of the geochemistry-mineralogy-texture approach and geometallurgical
techniques, in: The Second AusIMM International Geometallurgy Conference (GeoMet). Brisbane,
pp. 143–154.
Parian, M., Mwanga, A., Lamberg, P., Rosenkranz, J., 2018. Ore texture breakage characterization and
fragmentation into multiphase particles. Powder Technology 327, 57–69.

26
 doi:10.1016/j.powtec.2017.12.043
Partio, M., Cramariuc, B., Gabbouj, M., Visa, A., 2002. Rock texture retrieval using gray level co-occurrence
matrix. The Proceedings of the 5th Nordic Signal Conference 1–5.
Pérez-Barnuevo, L., Lévesque, S., Michaud, D., Bazin, C., Longuépée, H., 2016. Geometallurgical
characterization of drill core textural patterns: a case study from the Mont-Wright iron ore deposit,
in: IMPC 2016: XXVIII International Mineral Processing Congress Proceedings. Québec.
Pérez-Barnuevo, L., Pirard, E., Castroviejo, R., 2013. Automated characterisation of intergrowth textures
in mineral particles. A case study. Minerals Engineering 52, 136–142.
doi:10.1016/j.mineng.2013.05.001
Pérez-Barnuevo, L., Pirard, E., Castroviejo, R., 2012. Textural descriptors for multiphasic ore particles.
Image Analysis and Stereology 31, 175–184. doi:10.5566/ias.v31.p175-184
Petersen, L., Minkkinen, P., Esbensen, K.H., 2005. Representative sampling for reliable data analysis:
Theory of Sampling. Chemometrics and Intelligent Laboratory Systems 77, 261–277.
doi:10.1016/j.chemolab.2004.09.013
Peyré, G., 2011. The Numerical Tours of Signal Processing. Computing in Science & Engineering 13, 94–97.
doi:10.1109/MCSE.2011.71
Pirard, E., 1993. Doctoral thesis : Morphologie euclidienne des figures planes, application à l’analyse des
matériaux granulaires. Université de Liège.
Pirard, E., Lebichot, S., Krier, W., 2007. Particle texture analysis using polarized light imaging and grey level
intercepts. International Journal of Mineral Processing 84, 299–309.
doi:10.1016/j.minpro.2007.03.004
Ramanaidou, E.R., Wells, M.A., 2011. Hyperspectral Imaging of Iron Ores, in: Proceedings of the 10th
International Congress for Applied Mineralogy. pp. 1–5. doi:10.1007/978-3-642-27682-8
Schodlok, M.C., Whitbourn, L., Huntington, J., Mason, P., Green, A., Berman, M., Coward, D., Connor, P.,
Wright, W., Jolivet, M., Martinez, R., 2016. HyLogger-3, a visible to shortwave and thermal infrared
reflectance spectrometer system for drill core logging: functional description. Australian Journal of
Earth Sciences 63, 929–940. doi:10.1080/08120099.2016.1231133
Serra, J., 1982. Image Analysis and Mathematical Morphology. doi:10.2307/2531038
Su, D., Yan, W.M., 2018. 3D characterization of general-shape sand particles using microfocus X-ray
computed tomography and spherical harmonic functions, and particle regeneration using
multivariate random vector. Powder Technology 323, 8–23.
doi:https://doi.org/10.1016/j.powtec.2017.09.030
Taylor, R., 2010. Ore textures: recognition and interpretation. Springer Science & Business Media.
Tiu, G., 2017. Master thesis : Classification of Drill Core Textures for Process Simulation in Geometallurgy:
Aitik Mine, New Boliden. Luleå University of Technology, Luleå, Sweden.
Togami, S., Takano, M., Kumazawa, M., Michibayashi, K., 2000. An algorithm for the transformation of XRF
images into mineral-distribution maps. Canadian Mineralogist 38, 1283–1294.
doi:10.2113/gscanmin.38.5.1283
Tsafnat, N., Tsafnat, G., Jones, A.S., 2009. Automated mineralogy using finite element analysis and X-ray
microtomography. Minerals Engineering 22, 149–155. doi:10.1016/j.mineng.2008.06.003

27
van der Sanden, J.J., Hoekman, D.H., 2005. Review of relationships between grey-tone co-occurrence,
semivariance, and autocorrelation based image texture analysis approaches. Canadian Journal of
Remote Sensing 31, 207–213. doi:10.5589/m05-008
Wicks, F. j, Corbeil, M.-C., Back, M.E., Ramik, R.A., 1995. Microbeam X-ray diffraction in the analysis of
minerals and materials. The Canadian Mineralogist 33, 313–322.
Wiener, N., 1930. Generalized harmonic analysis 117–258. doi:10.1007/BF02546511
Wightman, E., Evans, C., Tungpalan, K., 2014. Measurement and interpretation of textural features at
meso-scale, in: Proceedings of the 3rd International Symposium on Process Mineralogy (Process
Mineralogy’14).
Williams, S., 2013. A Historical Perspective of the Application and Success of Geometallurgical
Methodologies. The second AUSIMM international geometallurgy conference / Brisbane, QLD, 30
September - 2 October 2013 37 37–47.
Wills, B., Napier-Munn, T., 2005. Mineral Processing Technology, An Introduction to the Practical Aspects
of Ore Treatment and Mineral Recovery, 7th ed, Wills’ Mineral Processing Technology (Seventh
Edition).
Yu, Z., Cui, Y., Wang, X., Zhang, Y., 2010. Identification of ores and gems using THz polarization information,
in: Proceedings of SPIE 7854, Infrared, Millimeter Wave, and Terahertz Technologies. pp. 78541U-
78541U–13. doi:10.1117/12.870541
Zhang, D., Lu, G., 2004. Review of shape representation and description techniques. Pattern Recognition
37, 1–19. doi:10.1016/j.patcog.2003.07.008

28
Paper III
 Data fusion for enhanced prediction of geometallurgical variables

Pierre-Henri Koch*,a, Cecilia Lunda, Jan Rosenkranza

* Corresponding author: [email protected]

a Luleå University of Technology, Division of Minerals and Metallurgical Engineering

Abstract
The integration of geological and metallurgical information has led to the definition of geometallurgical
ore types which can be characterized by different mineralogical and mineral processing properties.
However, the classification into geometallurgical ore types is often done manually and can be subject
to misclassification. To integrate geometallurgical variables acquired at different scales, a data fusion
concept has been developed which allows a more reliable and objective ore type classification. The
approach is based on convolutional neural networks for image-based information and dense neural
networks for vector information. The data set used for training and validation is based on an apatite
iron ore deposit and includes metallurgical test results from grindability tests as well as from wet low
intensity magnetic separation at laboratory scale. Various levels of noise and missing data were tested
as different data fusion scenarios. The proposed algorithm allows decreasing the subjectivity of ore
type classification and can be used to incorporate geometallurgical ore classes into an ore block model.

1 Introduction
Geometallurgy aims at building a predictive, spatial model that describes both, the geological and the
processing related variables of an ore body. Input parameters to these models usually comprise data
sets coming from different types of analyses, involving characterization methods as drill core logging,
automated mineralogy, metallurgical testing etc. Several challenges arise from the different types of
information as well as from the large quantity of data. Among these challenges, the aggregation of
data from different sources and measured at different scales has been identified as an important and
emerging field in geometallurgy where little work has been done so far (van den Boogaart and
Tolosana-Delgado, 2018). In the following, the problem of data fusion for geometallurgical variables is
addressed, one case study is provided, and the relevant tools and concepts are introduced.
In the context of an apatite iron ore deposit, data has been compiled to build a predictive model for
both, iron grade and recovery. Drill cores have been collected, cut in half, assayed and high-resolution
images were taken of the cut plane. A previous study has shown how the chemical assays could be
used to build such a model (Koch et al., 2019). But images could not be incorporated using classical
statistics. The objectives of this study are (i) to describe the data fusion concept, i.e. architecture of
models that accept both images and numerical data as inputs, (ii) to build and train such models for
predictions and (iii) to evaluate the model performance. The nature and strengths of these models are
briefly discussed as well as their potential applications in geometallurgy.

1
2 Materials and method
2.1 Deposit and samples
13 different ore types from the Leveäniemi apatite iron ore mine, located in Northern Sweden, were
sampled from drill cores and grab samples (>60 kg per ore type), representing the entire deposit. The
mine, operated by LKAB, is an open pit mine that was in operation between 1964 and 1983. Due to an
increasing market, production resumed in 2015. The mine is expected to operate until 2035. The
reserve (proven and probable) is estimated at 113 Mt with an average iron content of 44 % (LKAB,
2018)
The deposit is dominated by massive magnetite, but large volumes of semi-massive ore occur in the
surrounding rocks in up to 100 m wide zones. The notation m- or s- refers to the massive and semi-
massive classification as used by the mine geologists. In the middle part of the deposit, magnetite is
altered to hematite. The massive ore is high in iron (>60%) and the typical low-grade magnetite ore
(<40%) shows a higher grade of silica, both of them hosted in felsic-mafic volcanic rocks (Lishchuk et
al., 2018; Martinsson et al., 2016). The main iron bearing mineral is magnetite, but hematite is locally
present. The main gangue minerals are amphiboles, biotite, potassium-rich feldspar and plagioclase as
the major gangue minerals, as well as apatite and calcite. Traditionally, the ore has been described as
either massive or semi-massive (breccia). But recent studies have shown more subtle differences in
terms of chemistry, mineralogy and behavior in metallurgical tests (Gustafsson, 2016; Lishchuk et al.,
2018).
In Figure 1 and Table 1, the different samples are presented and assigned to an ore class from 1 to 13.
The samples are described according to their rock type but also by the massive/semi-massive
classification. Several sub-samples of each type were assayed, yielding an average iron grade and a
standard deviation. These grades will be considered as part of the input to the model.

Figure 1: Modal mineralogy for the samples (m: massive, s: semi-massive). The composition is
expressed in wt%.

2
 Table 1 : Description of the samples
Standard
Name Rock type Description Average iron grade [wt.%] deviation
[wt.%]
Biotite schist, Semi-massive,
1 30.72 2.21
dyke foliated
Massive,
2 64.84 1.80
noticeable weathering
3 Massive 62.13 3.00

4 Massive, foliated 65.70 1.97

Trachyandesite
Massive,
apatite rich,
5 60.39 3.09
biotite,
foliated
Semi-massive,
6 26.11 0.41
foliated
7 Massive 66.90 0.92
8 Massive 62.20 2.96
Trachyandesite,
Semi-massive,
9 near biotite schist 37.90 1.92
highly foliated
Massive,
calcite veining,
10 62.20 1.79
brittle,
amphibole
Trachyandesite
Massive,
11 64.16 1.84
biotite
12 Rich semi-massive 46.28 1.74
13 Semi-massive 16.16 3.43

3
In addition to the iron grade, approximately 5 meters of each ore type were imaged in the visible light
spectrum, using two light sources and a high-resolution camera. The motivation for using macro-scale
textures has been that these are cheap to acquire and may contain additional information about the
mineralogy or the behavior of the ore in the process (Koch et al., 2017; Lund et al., 2015). Figure 2
shows an example of two typical ore types, one semi-massive.

Figure 2 : Ore types 7 (left: massive) and 6 (right: semi-massive)

1 cm

2.2 Metallurgical tests

Three types of metallurgical tests were conducted on the samples. A small-scale test, called the
Geometallurgical Comminution Test (GCT), described in Mwanga et al. (2015), was used to estimate
the Bond work index (BWi) of the ore. The test is based on dry grinding with a standard charge in a
small CAPCO ball mill while measuring the energy consumption and change in particle size at given
time intervals. The main advantage of such a test is the small amount of material and the shorter time
needed to perform it compared to a full-scale bond work test. The second metallurgical test aimed at
determining the maximal grade of the concentrate using a Davis Tube (DT). A three-stage process with
0.1, 0.2 and 0.5 A and a duration of 130 seconds per run produced a magnetic concentrate that was
assayed. The third metallurgical test was a laboratory-scale wet low-intensity magnetic separation test
(WLIMS) with a flow of 0.5 kg/min for solids and 3l/min for water. The recovery of iron in the
concentrate was calculated from the products. Figure 3 illustrates the procedure.

METALLURGICAL TESTS

Crushing -3.35 mm

Grinding ~80µm P80 Grinding in GCT for 11

min

WLIMS Davis tube test

Figure 3 : Summary of the metallurgical tests. Orange boxes show the input of the model.

4
Each test covers distinct aspects of the variability expected in the orebody, i.e. variations of the
comminution properties, the grade of the concentrate and the recovery of iron. These variables have
been determined in 3 to 6 repetitions with sub-samples to estimate their mean and their standard
deviation, compare Tables 2 to 4.
Table 2: Bond work index from Geometallurgical Comminution Test (GCT, kWh/t)
Sample Mean [kWh/t] Standard deviation [kWh/t]
1 9.162 0.456
2 7.020 0.361
3 6.882 0.353
4 7.651 0.375
5 6.078 0.308
6 9.648 0.499
7 8.347 0.421
8 5.336 0.265
9 9.086 0.468
10 6.602 0.328
11 8.602 0.447
12 7.141 0.351
13 7.804 0.393

Table 3: Grade for different field strengths (Davis tube, wt.%)

Sample Mean [Fe wt.%] Standard deviation [Fe wt.%]
1 63.030 0.009
2 69.749 0.011
3 66.937 0.010
4 70.343 0.010
5 69.143 0.011
6 62.846 0.010
7 71.319 0.010
8 69.656 0.011
9 65.773 0.010
10 69.963 0.011
11 70.255 0.011
12 68.008 0.010
13 59.452 0.009

5
 Table 4: Recovery from WLIMS tests (wt.%)
Sample Mean [Fe Rec wt.%] Standard deviation [Fe Rec wt.%]
1 89.995 0.013
2 82.538 0.012
3 54.881 0.008
4 89.916 0.013
5 89.148 0.014
6 92.439 0.014
7 84.655 0.013
8 86.356 0.013
9 91.158 0.014
10 90.478 0.014
11 79.761 0.012
12 93.787 0.013
13 95.331 0.014

2.3 Modelling
The purpose of the developed model is to provide an estimate of the results from the laboratory tests
(i.e. GCT BWi [kWh/t], DT iron grade [wt.%] and WLIMS recovery [wt.%]) based on the iron grade of
the sample and image patches from drill cores of the same ore type. The general structure of the model
calculation is shown in Figure 4. The role of each step is briefly introduced in the following sub-sections.

2.3.1 Inputs to the model

Inputs to the model include chemical assays of the drill cores and metrics derived from the images of
half drill cores. Each image is represented as a matrix of double precision values with 3 channels (red,
blue, and green), stored as a TIFF file without any compression to avoid artifacts and noise.

2.3.2 Pre-processing of data and images

The pre-processing involves cropping for removing unwanted areas from the images and fitting of the
assay data to a normal law. The image size was limited to patches of 128 pixels by 128 pixels with 3
color channels. Then, a corresponding assay value from the drill core is assigned to each patch. The
results are presented as frames containing one image patch and a single assay value. It is important to
note that there is no feature extraction before using images to the model as this step will be done by
the model itself. Given the size and number of images, a total of 150 000 frames were generated,
evenly distributed between the 13 ore classes. At this stage, the dataset was already split into a training
set with 80 % of the total number of frames and a testing set with 10% of the data while a validation
set not used to train the network but only to test the ability of the model to make accurate predictions
based on unseen data.

6
Figure 4: General structure of the model calculations

7
2.3.3 Modelling approach
The model itself is a large artificial neural network referred to as ANN or NN (Hebb, 1949; White and
Rosenblatt, 2006). A neural network belongs to the group of supervised machine learning models that
can be used to approximate an arbitrary function. Supervised means that when the network is trained
on the combination of correct inputs and outputs (the training set) then the network can adapt its
parameters to produce an output from the testing data set, which is not too different from the output
of the training set.
In analogy to biological neural networks, artificial neural networks consist of layers of neurons
connected to each other. Each neuron has a weight wi and the weighted sum of the inputs to each
neuron with these weights produces the input to the next layer of neurons. The last layer of the
network produces the value of the function to estimate. In this case, the values to be estimated are
the DT iron grade, the WLIMS recovery and the BWi estimate.
The learning or training phase in neural network is based on the computation of the error (or loss)
function that measures the difference between the true value of a point of the training set with the
predicted value. The most commonly used method to minimize the loss is the stochastic gradient
descent in which the weights wi are modified according to the gradient of the error until the values of
the loss function become reasonably small (Hecht-Nielsen, 1989). A risk associated with training a
neural network is overfitting. This situation can happen when the number of parameters of the model
is large enough to memorize the training data and lose any capacity to generalize.
A variant of neural networks, in which every neuron is not connected to each other but seen as a series
of convolutional layers, is called a convolutional neural network (CNN or convnet). In a CNN, the
weights to learn are organized as a set of convolution operators (or filters). The operation of
convolution on images offers many advantages over traditional (or dense) neural networks:
• It reduces the number of weights to learn
• It preserves the geometry of the input: in the case of an image, the pixels keep their location,
• It is also invariant to a spatial shift of the input image.

Instead of carefully designed filters, CNNs offer a way to learn which filters are interesting depending
on the task at hand (LeCun and Bengio, 1998). This feature of the model developed here differs from
previous work on drill core texture classification in the context of geometallurgy (Pérez-Barnuevo et
al., 2018). The novelty of this approach is to greatly simplify the quantification of textures since no
selection must be done regarding feature selection. Additionally, the images can be normalized by
subtracting the mean of the image from pixel values and divide by the standard deviation of pixel
values. Convolutional neural networks are also robust with respect to noise and variations in
brightness.
The convolution filters learned can be used as an input to a dense layer of neurons, making the
connection between a dense neural network and a convolutional neural network possible. The
structure of the designed net has two branches, one convolutional branch that takes the image patches
as inputs and has a dense layer or 32 neurons as an output, and the dense branch that takes the iron
grade as an input and has a dense layer of 32 neurons as an output. The two final layers of each branch
are concatenated to a dense layer of 64 neurons.

8
2.3.4 Model output
The concatenated layer of 64 neurons is linked to a dense layer of 3 neurons (one for each variable
that shall be predicted). The output is a vector of the three 3 values. Each of the last neurons provides
an estimate for the bond work index, DT grade and WLIMS recovery, respectively.

3 Numerical results
Evaluating the performance of the model approach was done by first training the network to 256
frames during the training phase, and second by using the produced model for geometallurgical ore
class predictions based on the 128 frames of the testing data set. The testing data set was not part of
the training set.

3.1 Training the model

Training a convolutional neural network requires a loss function which measures how far the
predictions are away from the training data during the training phase. Based on this function, the
gradient of the error is calculated, and the weights of the convolutional filters are adapted according
to a stochastic gradient descent (or one of its variants). The idea is to compute the weights that allow
the best prediction possible based on the training data. The quality of the model was evaluated with
the mean squared error. The final model had a mean squared error of and was selected as the final
model.
Training a model is a computationally expensive step and requires a recent graphical processing unit
(GPU) to perform well. In this case, the network was trained in 12 hours using a laptop and a recent
GPU. Depending on the data available and the complexity of the models, one can use cloud-based or
distributed computations on powerful servers to minimize the training time.

3.2 Model predictions for the validation set

The goal was to predict the results of the laboratory tests (GCT BWi [kWh/t], DT iron grade [wt%] and
iron recovery in WLIMS [wt%]) representative of the concentration process. A summary of the quality
of the predictions is presented in Table 5. The main measure of interest is the mean squared error
since it was used as the loss function for training the network. Its penalization of extreme values is
visualized in Figure 5, Figure 6 and Figure 7. The overall quality of the predictions is good but weaker
for the estimates of the Bond Work index from GCT tests as showed by a lower coefficient of
determination.
Table 5 : Summary of the error measures of the model
Measure GCT estimates DT grade estimates WLIMS recovery estimates
Mean absolute error 0.554 kWh/t 0.673 wt.% 1.791 wt.%
Mean squared error 0.586 kWh/t 0.8157 wt.% 6.843 wt.%
Explained variance 0.682 0.959 0.96
Coefficient of 0.674 0.932 0.96
determination

9
 Figure 5 : Scatter plot of GCT estimates

The predictions in the case of GCT values indicates that the model seems to predict values only
between 6 and 9.5 at least for the validation set. One reason could be the lack of training data with
values outside this range or the presence of unusual points in the testing set.

Figure 6 : Scatter plot of DT estimates

10
 Figure 7 : Scatter plot of WLIMS estimates

The vertical groups of points that can be seen on Figure 7 are the result of the measurements being
done per geometallurgical class. As a result, each sample of the same class has values close to the
average in the case of low standard deviations. The extremely low recovery values come from semi-
massive samples (breccia).

3.3 Predictions in presence of noise in the input

To evaluate how much the image information is needed for the predictions, different cases were
evaluated where normally distributed noise was added to the inputs (the image part of the frame or
the iron grade measurement). Table 6 shows that the model is not overly sensitive to noise since at
the 5 % level only minor changes are visible in the coefficient of determination. This is the case
irrespective of the image part.
Table 6 : Impact of the presence of noise in the input variables
Measure Value for Value for Value for WLIMS
GCT DT estimates
estimates estimates
Coefficient of determination without noise 0.674 0.932 0.960
Coefficient of determination with random 0.645 0.915 0.845
geometric transformation of images
Coefficient of determination with 5 % noise in 0.550 0.917 0.854
the images
Coefficient of determination with 5 % noise in 0.553 0.919 0.799
the iron assay

11
4 Discussion
The combination of geometric features relative to drill core textures and chemical information by the
developed data fusion concept improves the predictive capacity of the model. While the previous study
used textural features that were manually selected, the present study offers an alternative that not
only skip the manual selection step but also incorporates chemical information in the same model. The
chemical information about the drill cores represents the smallest level of detail. On the other hand,
images of drill cores have information from a very coarse scale. One interpretation of the model is that
data fusion happening inside the network integrates information from different scales. This shows a
direction for future developments for this class of models. One could include spatial information from
a coarser scale such as a block model and still manage to produce estimates of metallurgical responses
at a lab or full scale. The notion of geometallurgical ore types (or classes) is not explicit anymore but
still present in the weights of the global model. A careful study of the weights of the network could
reveal which textural features are learned from the data but also where the dimension has been
reduced the most, i.e. how many ore types are used by the network for prediction.
A major weakness inherent to supervised learning methods is the amount of data needed to train the
model and its capacity to generalize from the examples to unseen data (the generalization problem).
In this case, chemical assays and visible light images were readily available in sufficient quantity but
this may not be the case for other types of data. The validation of the presented model was done on a
128 randomly sampled frames and sample size could have a significant impact on the error measures.
The sum of the mean-squared losses for all predicted variable in the training set was lower than 10.
This value is the same as the one measured during testing as presented on Figure 8.
General guidelines about the structure of convolutional neural networks (CNN) can be found in
literature (Szegedy et al., 2016). However, specifications of the CNN rely on trial-and-error, i.e., the
architecture of the different neuron layers, the dimension of the convolution filters, the connections
between layers and the choice of the optimization technique. In this case, a simple grid search
algorithm was performed using a python package called Keras for training the model and another one
named Talos to search for the optimal parameters (Chollet and others, 2015; Kotila, 2019).
Training the model requires the minimization of the error. This is usually done via a stochastic gradient
descent algorithm (Bottou, 2010) or one of its many variants like Adam (Kingma and Ba, 2014).
The obtained model is stored as one file that contains both the structure of the network and its
weights. Combined with a client-server interface, the model can be used in an industrial context where
the server will host the model and the client will receive the predictions. To train such a model, two
approaches are available: a local one where the server that is trained is located on site or a remote
server. The advantages of the first setting are to offer access to hardware in case of failure, but the
main weakness is the presence of the computer inside an industrial environment potentially in
presence of dust and vibrations. In the second setting, one limitation is the available data bandwidth
for sending and receiving the batches of data over Internet. The training of the model should be done
regularly on a need basis or when the textural ore types change dramatically. Regular incremental and
full backups can be done with any modern deep learning platform but a strategy on data and model
maintenance is advisable.

12
Figure 8 : Evolution of the error as a function of the number of training step. Each step provides a batch
of 32 frames to the model. On x axis is the value of the mean squared error and on the y axis the number
of steps.
In addition, the features that the model learns can be analyzed and can serve as a basis for genetic and
lithology studies. Since convolution filters can be seen as patterns that produce a high value at a pixel
that presents a neighborhood similar to the filter, the inspection of the filters learned by the model
will display the patterns the model has identified as relevant for prediction. Figure 9 shows one image
patch from class 3 after successive convolutions in the model. One can see different patterns for
example veins, vertical structures or disseminated phases that are analyzed by the filters.

13
Figure 9 : Examples of the result of the convolution of one sample image with some of the convolution
filters learned during training. Low values are indicated by a blue color and high values tend to red.

14
5 Conclusion
The study has shown that it is possible to design, train and evaluate the predictive capacity of a
geometallurgical model that accepts both image and vector data within a data fusion concept. The
performance of the model depends on both the quality of the input data and the architecture of the
network. In its core, the capacity of neural networks to concatenate or add the weights of two different
branches is used for data fusion in the context of geometallurgy. While the design of such networks is
mostly empirical, the general structure of a multi-input deep network demonstrated a good capacity
of prediction.
The proposed concept is contributing to filling the existing gaps in efficient data processing and
numerical modelling for geometallurgy. It can be easily extended to other data types and additional
data sources. The selected models even allow uncertainty in the data. With the increasing
computational power available to both researchers and the mining industry, it can be expected that
such models become more common in geometallurgical applications. In an extension, spatial
coordinates could be used as part of the input data and the model could be used for spatial
interpolation purposes based on drill core images but also tomography and chemical assays.

6 Acknowledgments
The authors wish to thank LKAB, Boliden Mineral AB and Outotec (Sweden) AB for their support in this
study, Viktor Lishchuk for the interesting discussions and suggestions about machine learning. This
work is part of PREP project funded by VINNOVA within the Strategic innovation program for the
Swedish mining and metal producing industry (SIP-STRIM).

7 References
Bottou, L., 2010. Large-scale machine learning with stochastic gradient descent, in: Proceedings of
COMPSTAT’2010. Springer, pp. 177–186.
Chollet, F., others, 2015. Python package : Keras.
Gustafsson, M., 2016. Master thesis: Element to Mineral Conversion of LKAB Deposits Element to
Mineral Conversion of LKAB deposits: Leveäniemi. Gruvberget Hematite and Sparre. Master
Thesis. Luleå University of Technology.
Hebb, D.O., 1949. Organization of Behavior. Wiley, New York.
Hecht-Nielsen, R., 1989. Theory of the backpropagation neural network, in: Proceeding of the
International Joint Conference on Neural Networks. pp. 593–605 vol.1.
doi:10.1109/IJCNN.1989.118638
Kingma, D.P., Ba, J., 2014. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
Koch, P.-H., Lund, C., Lishchuk, V., Rosenkranz, J., Lamberg, P., 2017. Automated Classification of Drill
Core Textures for Geometallurgy, in: Proceedings of the Process Mineralogy Conference, Session
5. MEI, Cape Town.
Koch, P.-H., Lund, C., Rosenkranz, J., 2019. Automated drill core mineralogical characterization method
for texture classification and modal mineralogy estimation for geometallurgy. Minerals
Engineering 136, 99–109. doi:10.1016/J.MINENG.2019.03.008
Kotila, M., 2019. Python package: Talos.
LeCun, Y., Bengio, Y., 1998. The Handbook of Brain Theory and Neural Networks, in: Arbib, M.A. (Ed.),
. MIT Press, Cambridge, MA, USA, pp. 255–258.
Lishchuk, V., Lund, C., Koch, P.-H., Gustafsson, M., Pålsson, B.I.I., 2018. Geometallurgical

15
 characterization of Leveäniemi iron ore - Unlocking the patterns. Minerals Engineering 131, 325–
335. doi:10.1016/j.mineng.2018.11.034
LKAB, 2018. LKAB’s annual and sustainability report 2018. Luleå, Sweden.
Lund, C., Lamberg, P., Lindberg, T., 2015. Development of a geometallurgical framework to quantify
mineral textures for process prediction. Minerals Engineering 82, 61–77.
doi:10.1016/j.mineng.2015.04.004
Martinsson, O., Billström, K., Broman, C., Weihed, P., Wanhainen, C., 2016. Metallogeny of the
Northern Norrbotten Ore Province, Northern Fennoscandian Shield with emphasis on IOCG and
apatite-iron ore deposits. Ore Geology Reviews. doi:10.1016/j.oregeorev.2016.02.011
Mwanga, A., Rosenkranz, J., Lamberg, P., 2015. Testing of Ore Comminution Behavior in the
Geometallurgical Context—A Review. Minerals 5, 276–297. doi:10.3390/min5020276
Pérez-Barnuevo, L., Lévesque, S., Bazin, C., 2018. Automated recognition of drill core textures: A
geometallurgical tool for mineral processing prediction. Minerals Engineering 118, 87–96.
doi:https://doi.org/10.1016/j.mineng.2017.12.015
Szegedy, C., Vanhoucke, V., Ioffe, S., Shlens, J., Wojna, Z., 2016. Rethinking the inception architecture
for computer vision, in: Proceedings of the IEEE Conference on Computer Vision and Pattern
Recognition. pp. 2818–2826.
van den Boogaart, K.G., Tolosana-Delgado, R., 2018. Predictive Geometallurgy: An Interdisciplinary Key
Challenge for Mathematical Geosciences, in: B.S. Daya, S., Qiuming, C., Frits, A. (Eds.), Handbook
of Mathematical Geosciences: Fifty Years of IAMG. Springer International Publishing, Cham,
Switzerland, pp. 673–686. doi:10.1007/978-3-319-78999-6
White, B.W., Rosenblatt, F., 2006. Principles of Neurodynamics: Perceptrons and the Theory of Brain
Mechanisms. The American Journal of Psychology 76, 705. doi:10.2307/1419730

16
Paper IV
 Sequential decision-making in geometallurgy
Pierre-Henri Koch*,a, Jan Rosenkranza
* Corresponding author: [email protected]

a Luleå University of Technology, Division of Minerals and Metallurgical Engineering

SE-97187 Luleå, Sweden

Abstract

Geometallurgy as a multi-disciplinary field has been applied at various levels in different operations. By
linking the ore performance in mineral beneficiation processes to the ore block model, it supports
estimating the value of a block before it is mined. Efforts in the classification of the ore into
geometallurgical classes have led to a better understanding of the entire value chain. While classification
provides a convenient tool for forecasting and visualization purposes, it simplifies the actual complexity
of an ore body. In mining and process planning, sequential decisions are made to maximize an objective
function or equivalently minimize a regret function. Using available information from geology or
metallurgical test work, an optimal strategy can be found using tools from the machine learning
community.

In this study, a framework based on machine learning to maximize the use of such classifications for
sequential decision-making is proposed. The concepts of reinforcement learning and bandit algorithms
offer powerful tools to explore and exploit different optimization strategies. In certain cases, theoretical
guarantees about the performance of given methods can be obtained by regret bounds.

Based on existing models of a porphyry copper deposit and an iron ore deposit, this study presents a
methodology and different available algorithms to maximize an objective function that depends on a high
number of variables and in the presence of noise or uncertainty in the models. Different numerical
experiments provide a basis for discussion and comparison to human decisions. The hypotheses relative
to each algorithm are discussed in relation to the mineral processing models.

1 Introduction
Geometallurgical models are often defined as a set of tools to improve knowledge of an ore body by
measuring relevant information for the concentration process (Lund and Lamberg, 2014). The objective is
to gain information but also to use this information to improve the grade or recovery of the mineral of
interest, optimize the throughput or reduce the environmental footprint of the process. In a
geometallurgical program, a large amount of data from different sources is acquired, structured and made
available as a basis for decision-making for mining and process planning. Sequential decisions have to be
made in order to maximize an objective function, usually aggregating economic objectives, production
and environmental objectives, or to minimize a regret function, respectively. Besides the selection of
different processing options, these decisions often involve dealing with uncertainties with respect to the
variations of the ore feed and other input variables, model parameters and the model structure. A special
challenge arises from cases where certain information is not available yet, i.e. the decision must be taken

1
a priori to knowing the outcome of the selected processing option. Numerical methods based on bandit
models have been proposed for this type of problems but, to the authors’ knowledge, have not been
made available to the field of geometallurgy yet.

The first objective of this paper is to study the dimensions and regularity of the available data, a second
objective is to identify and present strategies as algorithms that offer a good performance in terms of
recommendations for mineral processing. Finally, the implications of the investigated methods as well as
their limitations for geometallurgy are presented.

1.1 Stochastic bandit model

In a situation where several ore bodies are mined, and several processing options exist, as shown on Figure
1. The steps between the mines and the beneficiation process (i.e. steps as ore handling, transportation
and blending) are neglected as a first approximation. A further simplification is to consider only the grade
of one mineral of interest and that every ore body has the same grade distribution resulting in the same
feed properties to all the processes. Additionally, all the processes are considered stationary with respect
to time. The actual performance of each process is unknown but represented as a statistical distribution.

The task to be solved is then to maximize an economic objective function through assigning an ore block
at a time (sequentially) to the processing option that has the highest expected payoff while not knowing
the performance of the other processes.

At each discrete time step, the virtual plant operator (or software agent) can choose one of the processes
in which the ore block will be processed. Once the block has been processed, a performance indicator
relative to the chosen process (e.g. grade and recovery of the concentrate, presence or absence of a
deleterious element) is observed by the plant manager (also called agent since it can act on the system).
With time, the information about each process increases. The model describing this procedure is
commonly called “multi-armed bandits” (MAB) setting, by analogy to a player in a casino who can choose
between several slot machines equipped with levers called arms.

2
 Figure 1 : Bandit model of the resource allocation in mineral processing

1.2 Stochastic bandit problem and solution strategies

If a bandit model is chosen, a natural aim is to develop a strategy (or policy) to maximize the overall
performance in the long run. An alternative aim would be to identify the best process with high
confidence. In both cases, the concepts of reinforcement learning are applied, i.e. the acquisition of new
knowledge from the process (so-called exploration) and the subsequent improvement of decisions based
on the existing knowledge (so-called exploitation). At the beginning, decisions must be taken without
knowledge about any of the processes. The software agent must explore the different processes and,
once the number of observations increases, exploit the correct process. The objectives used within
process optimization and decision-making can be either of binary or continuous type, thus requiring
different mathematical formulations for the description of the observations, see Table 1.
Table 1: Objective types

Objective Parametric Non-parametric Examples of geometallurgical

distributions variables
Binary: 1 if the Bernoulli Discrete function Presence or absence of a penalty
product is within in {0,1} element in the concentrate,
specifications, meet or not environmental
otherwise 0 requirements, revenue higher
than a specified threshold or a
logical combination of these
Continuous: a score in Exponential family Bounded Grade or recovery of the
[0,1] that represents continuous concentrate, throughput,
the performance of function in [0,1] amount of penalty elements, or a
the process combination (sum) of these
(normalized between 0 and 1).

To describe the problem, the usual terms and notation from the bandits’ literature are used:

3
 - At each time step t = 1, …, n the agent (plant manager) choses an arm (sends a mass of ore to one
process or another) At ϵ {1,…,K}
- The agent observes a reward Xt (performance indicator for the process) from the chosen arm At.
The rewards are independent and identically distributed. The reward measures in this case the
performance of the chosen process.
- The policy ∏ of the agent is the decision rule that selects the choice At based on the previous
observations (A1,X1,…,At-1,Xt-1)
- The best arm is noted a* and its mean reward is µa*
- The overall performance of a policy can be measured by the cumulative regret RT which is the sum
of the differences between the rewards obtained by the agent and the rewards that the agent
could have obtained by always selecting the arm yielding the highest reward

The challenge is to build “good” policies for a specific class of rewards. This has been approached by Lai
and Robbins (1985) for a single parameter exponential family and extended in the works of Burnetas and
Katehakis (1997). In short, to achieve a cumulative regret of the order of tα with α ϵ [0,1], any suboptimal
arm that has a mean µa must be selected on average more often than ln (t)/KL(µa, µa*), where KL is the
Kullback-Leibler divergence that quantifies the difficulty of making a choice between two arms. These
policies are called asymptotically optimal which means that this property is desirable in practice since it
provides bounds for the cumulative regret when the number of time steps is large.

One of the early examples of designing an optimal strategy (or policy) in a sequential decision process can
be found in the works of Thompson (1933). The problem addressed is to determine whether a given
medical treatment is better than another. The main hypothesis is to model the success or failure of a
treatment as a Bernoulli trial. The solution adopted is to model a-priori probabilities, draw a sample from
both distributions, choose the treatment that has the highest probability of success and update the a-
priori distributions accordingly. This approach is Bayesian, i.e. unknowns are described by probability
distributions while the observations are non-random – as opposed to frequentists methods. Thompson’s
sampling is asymptotically optimal (Kaufmann et al., 2012).

When different processing paths exist, two extreme strategies can be adopted: either randomly explore
all the possible choices (pure exploration) or choose one process and always choose it (pure exploitation).
In the first case, one can see that the optimal process will be chosen on average as much as any sub-
optimal process. In the second case, either the best process is selected, and, in that case, it is a successful
strategy, or another process is chosen, and the reward will always be suboptimal. Additionally, the
probability to choose the optimal process from the beginning decreases with the number of available
choices. This is referred to as the exploration-exploitation tradeoff. One can see that a reasonable strategy
should combine exploration and exploitation.

A naïve yet effective approach to the exploration-exploitation tradeoff is the epsilon-greedy strategy in
which the arm with the highest mean is selected with probability (1-ε)/K and the others with a probability
ε/K where K is the number of arms. Epsilon is often chosen around 0.1 which means that the probability
of choosing the arm with the highest mean reward is 0.9 and 0.1 for the other arms. Epsilon-greedy
strategies solve the tradeoff by assigning a probability for exploration and exploitation.

Another family of policies are based on the construction of intervals of confidence on estimates of the
mean of each arm. As an example, the algorithm UCB chooses the arm that has the highest upper
confidence bound (hence UCB). The construction of the interval itself is detailed in Auer and Fischer (2002)

4
but in essence the upper bound consists of the sum of an exploitation term µ𝑎 (𝑡) (empirical estimate of
2log(t)
the mean) and an exploration term √ where Na is the number of times the arm “a” was selected.
N a (t)

This approach is frequentist since it is based on an estimation of the mean that improves when the arm is
selected. A variant based on UCB and the Kullback-Leibler divergence called KL-UCB offers better
guarantees in terms of regret (Garivier and Cappé, 2011).

1.3 Applications in mineral processing

Applications of the algorithms explained above to mineral processing and geometallurgy are possible in
several cases presented in Table 2. A distinction is made between time stationary (S) and non-stationary
(NS) processes and ore classes. The former represents a process that has a steady performance in time.
This assumption holds when a short time is considered but when longer times are studied, the process is
often non-stationary. With regard to ores, the distinction can be made between a constant ore class (one
class of constant physical properties), a steady or targeted ore class (where the feed to the process has
the same physical properties on average) and k S-classes where the ore is split into different
geometallurgical classes representing different process performances, compare the examples given in
Table 1.
Table 2 : Different cases for resource allocation in geometallurgy

Ore type 1 S-process 1 NS-process k S-processes k NS-processes

1 class 1 S-simulation 1 NS-simulations k S-simulations k NS-simulations
1 S-class Simulation Simulation scenarios Case 1 Case 2
scenarios
k S-classes Simulation Simulation scenarios Case 3 Case 4
scenarios

As presented in Table 2, one can simplify these cases and evaluate whether commercially available
software offers a solution. Conceptually, a difference can be made about simplifications of the ore type
and the process classification. This is supported by previous work in geometallurgy (Lishchuk et al., 2015).
Classification into ore types has been done using dimension reduction techniques such as principal
component analysis (PCA) (Andrusiewicz et al., 2014; Keeney, 2010). Linear dimension reduction offers
the advantage that distances have an intuitive meaning in terms of physical properties such as mineralogy.
Non-linear mappings can also be used and theoretical guarantees exist regarding the existence of such a
mapping and the preservation of distances (Kłopotek, 2017).

A stationary process has the same statistical distribution with time. In mineral processing plants, this
represents stable operating conditions, i.e. a plant that does not display drastically different average
behavior in time. An example of a stationary signal is a constant measurement affected by statistical noise.
An example of a non-stationary process is a plant affected by seasonal factors like temperature (Lin, 1989).

Several of these cases can be solved using existing numerical tools, namely mineral processes simulators.
Current software allows the simulation of stationary or non-stationary (also called dynamic) processes.
Additionally, different scenarios can simulate changes in the feed composition and process parameters
(Outotec, 2016). The question of optimal allocation is not directly addressed in these simulations,
moreover, the construction and calibration of such models relies on sampling and mass balancing, and

5
simulations are based on unit models that often include an empirical term (King, 2012). While empirical
terms and parameters allow calibration, they depend on measurements. Measurements can be regarded
as estimates of a real probability distribution. This hypothesis points to the use of stochastic settings to
model mineral beneficiation processes.

The cases where multiple processes exist are more complex (cases 1 to 4 in Table 2). From now on, no
difference will be made between simulations of processes (used for example during the design phase of
a process) or an actual plant (in that case, the strategies can either complement the design of experiments
or act as a guide to support decisions). If an interface exists between the plant data (or the simulation
software) and the implementation of algorithm for the strategies, the methods introduced can be used in
both cases. When several processing routes exist, the number of parameters in each plant increases,
which leads to either more unknowns (in the case of actual plants) or more calibration data required (in
the case of simulations). If the plant parameters are controlled by an automation system or if statistical
variability is introduced in the plant model, the problem becomes intractable. To overcome this difficulty
when decisions are made, bandit algorithms are useful.

2 Processes and algorithms

2.1 Industrial case studies and data acquisition
The data sets are describing two different industrial case studies:

a) Iron oxide ore: an iron ore deposit for which the modal mineralogy and chemistry of the feed is
available and where the ore has been classified into 8 geometallurgical classes or ore types as
presented in Table 3. Details regarding the geometallurgical characterization and classification can be
found in (Lishchuk et al., 2018). The classification is based on the Bond Work Index, Wet low-intensity
magnetic separation (WLIMS) recovery and grade of the Davis Tube concentrate (DT).

6
 Table 3 : Iron ore samples from Leveäniemi (Lishchuk et al., 2018)

Name Type Rock type Ore types Description of the selected sample

Biotite schist,
1 3, 5 Semi-massive, foliated
dyke

2 5 Massive, noticeable weathering

3 7 Massive

Massive, foliated, cleaned from: biotite schist

4
Drill Trachyandesite 5 (dyke), oxidation

5 core Massive, apatite rich, biotite, foliated

6 3, 6, 8 Semi-massive, foliated

7 6 Massive

8 Trachyandesite, 1 Massive
near biotite
9 schist 2 Semi-massive, highly foliated

10 5 Massive; calcite veining, brittle, amphiboles

11a 5 Massive, biotite

Grab Trachyandesite
11b 4 Rich semi-massive
samples
11c 3 Semi-massive

Names of the ore types adopted at Leveäniemi mine: (1) Western Massive Magnetite Zone >50%;
(2) Western Breccia Magnetite Zone <50%; (3) Eastern Massive Magnetite Zone <50%; (4) Eastern
Massive Magnetite Zone 50-60% (high grade breccia); (5) Eastern Massive Magnetite Zone >60%; (6)
Shallow Zone (Massive Magnetite); (7) Hematite Zone; (8) Low Grade Disseminated Mineralisation

7
 The process consists of wet low intensity magnetic separation with grinding stages, see Figure 2. A
process model was built and calibrated to experimental data in HSC Chemistry software and a
simplified model was derived. The simplification consists in constructing an analytical function that
allows a faster computation of the reward function compared to simulating the detailed flowsheet.
This is not a requirement for the method to work since the computation of a reward function can be
done in HSC Sim for example. The simplified models are detailed in section 2.2.

Figure 2 : Process flowsheet of the iron ore deposit at the time of sampling

b) Copper sulphide ore: a porphyry copper for which the mineralogy and chemistry of the feed is
available and where the ore has been classified into 15 geometallurgical classes or ore types as
presented on Table 4. Details regarding the geometallurgical characterization and classification can
be found in (Tiu, 2017). The classification is based on mineralogy, Bond Work Index and ore texture.
Table 4 : Porphyry copper samples

Name Type Rock type Ore types Mining unit

1 Feldspar-biotite-amphibole gneiss 11 Footwall
2 Quartz monzodiorite 10A,8,9
3 Micro-quartz monzodiorite 10B,12
4 Pegmatite veins 1 Dyke
5 Drill core Biotite/amphibole plagioclase 3A, 3B, 4,5,7 Main ore zone
gneiss
6 Quartz-muscovite-sericite schist 6
7 Feldspar-biotite-amphibole gneiss 2 Hanging wall
8 Meta-gabbro - Mafic intrusive

The process consists of froth flotation (rougher, scavenger and cleaner) with regrinding and
classification stages, compare Figure 3. A process model was built and calibrated to experimental data

8
 in HSC Chemistry software and a simplified model was derived. The simplified model is detailed in
section 2.2.

Figure 3 : Process flowsheet of the porphyry copper deposit at the time of sampling

2.2 Algorithms and software

MATLAB (Mathworks, 2016) has been used to implement the different decision-making policies and HSC
Chemistry 9 software (Outotec, 2016) has been used for running the unit-based process simulation. The
evaluation of the policies and statistical parts were conducted in MATLAB.

2.2.1 Case 1: Iron ore - Several steady state processing options, targeted feed grade
This case in an application of bandits to the iron ore example. The simplest situation is when the properties
of the feed ore are stationary; in practice, this is achieved by using the mine scheduling software to
stabilize the feed grade for example according to the iron and vanadium grade (Vos, 2018). From the
process point of view, one can construct a concentrate quality indicator that aggregates the presence of
desired minerals and the absence of detrimental elements (vanadium or silicate for example). The product
specification for silica can be modelled as a Bernoulli distribution: 1 if the concentrate is within the
requirements, 0 else. The agent can choose between 4 different processes modelled as

• R1=B(0.6): a circuit equivalent to the flowsheet on Figure 2 with a silica reverse flotation step to
remove any silicate from the concentrate
• R2=B(0.4): this is the circuit presented on Figure 2
• R3=B(0.3): a circuit equivalent to the flowsheet on Figure 2 with only two grinding stage
• R4=B(0.1): this circuit is equivalent to the flowsheet on Figure 2 with only one grinding stage

where 𝐵(𝜇) is a Bernoulli distribution of mean 𝜇

9
Since these distributions are unknown to the agent, different policies can be employed:

• A random selection of the process

• An epsilon-greedy strategy
• KL-UCB policy

The random selection represents an exploration-only strategy, obviously this should not be the optimal
method. The epsilon-greedy strategy is equivalent to a reasonable human strategy where the agent
empirically estimates which process to exploit but is still exploring every now and then. UCB is a strategy
based on the construction of intervals of confidence and selects the option that has the higher upper
confidence bounds (Agrawal, 1995).

2.2.2 Case 2: Iron ore - Several non-steady state processing options, targeted feed grade
When processes are not stationary, classical UCB policies are not appropriate anymore since the rewards
change with time (Hartland et al., 2006). A natural solution is to consider only a limited time during which
the confidence bounds are estimated. In that case, the policies forget the past after some time. As a result,
a sliding-window variant of UCB (SW-UCB) can be designed by choosing a time window. These algorithms
are described and analyzed in terms of optimality in Garivier and Moulines (2008). Works in non-
stationary bandits make a distinction between abrupt variations and smooth ones in the reward functions.
In mineral processing, this depends on the stability of the process. An example of smooth variations are
seasonal temperature differences affecting flotation and maintenance interruptions where the reward
could suddenly go to zero. As such, variations of the models used in Case 1 are introduced, with different
non-stationary behaviors. In this case, the iron ore concentration processes are not stable in time, due to
different operating conditions or instabilities.

• R1=Alternate with a period of 1000 time steps between B(0.5) and B(0.7)
• R2=Alternate with a period of 100 time steps between B(0.6) and B(0.8)
• R3=Alternate with a period of 150 time steps between B(0.5) and B(0.6)
• R4=Alternate with a period of 150 time steps between B(0.7) and B(0.9)

where 𝐵(𝜇) is a Bernoulli distribution of mean 𝜇

The adapted policies in this case are

• A random selection of the process

• A sliding-window epsilon-greedy strategy (SW-Epsilon-greedy)
• A sliding-window UCB policy (SW-UCB)

The SW-Epsilon-greedy strategy stands for the case of a periodical estimate of the best process on
average.

10
2.2.3 Case 3: Copper ore – Several steady state processing options, non-targeted feed grade
In this case, different ore types have been identified that behave differently during processing. Each ore
type has a direct impact on the process performance and therefore the available information should be
used by the policies. Bandit methods where side information is used are referred to as contextual bandits
and an extensive literature is available on the subject (Bubeck et al., 2012). Geometallurgical ore types
are supposed to behave in a similar way in the process since they were classified based on metallurgical
tests and texture. As such, two samples with close geometallurgical indicators should have the same
optimal process. A logical choice would be to look at the k-nearest neighbors (kNN) and use a similar
strategy. A relevant method referred to as kNN-UCB, based on nearest neighbors and UCB, offers a
straightforward implementation (Reeve et al., 2018).

Contextual data in geometallurgy can be a normalized vector that contains mineral grades, grindability
indicators, textural properties or a scaled version of a combination of all these properties.
Geometallurgical ore types are usually a classification of high-dimensional data. One method to reduce
dimension for classification purposes is k-nearest neighbors. As an example, Figure 4 displays contextual
data for the copper case (silicates grade, chalcopyrite grade, pyrite grade, grindability index and sulfur
grade) projected to three dimensions, with X, Y and Z being linear combinations of all five parameters. A
surface was fitted to the projected data to show that the data is supported on a lower dimension manifold.

Figure 4 : Representation of contextual data projected into a 3D space. The black circle indicates a neighborhood. All the points
inside the neighborhood are considered as k-nearest neighbors. The axis x, y and z are non-linear combinations of the initial
contextual data.

11
The kNN-UCB proceeds in several steps:

• For each arm (process), the number of neighbors k is chosen to minimize a measure of
uncertainty
• Given a context x at time t, it selects the k neighbors around x
• A local upper confidence bound on the reward functions is computed
• The arm with the highest bound is selected

Several hypotheses are needed to establish optimality bounds, i.e. hypotheses that suggest smooth and
bounded reward functions. For the problem stated here, a simplified set of assumptions involves:

• Dimension assumption: the context data is supported by a d-dimensional manifold where d is smaller
than the dimension of the context. This assumption means that the vectors that contain the
mineralogy and grinding indicators can be projected on a surface of lower dimension than the vectors
themselves as displayed on Figure 4.
• Margin assumption: the arms have rewards which are different enough
• Lipschitz assumption: for similar context vectors the rewards are similar. Another way to interpret the
assumption is that the reward functions should be varying smoothly (in a continuous way).
• Rewards are bounded in [0,1]

The application to the porphyry ore deposit is straightforward: a context is built by normalizing a vector
that contains the pyrite grade p (wt%), the chalcopyrite grade c (wt%), the sulfur grade s (wt.%) and a
grindability index g (kWh/t), giving the entire context vector denoted by x = (p, c, s, g). For each of the 15
ore types, a multivariate normal distribution P has been fitted based on the intra-class variability of the
samples. At each time step t one of the 50 000 blocks belonging to one of the 13 classes is observed as
context by the algorithms. In this case, the different options are different settings in the flotation circuit.
The reference case displayed on Figure 3 is modelled by R3 whereas increased grinding alternatives (i.e.
finer grinding resulting in better liberation) are represented as R1 and R2 and suboptimal operation as R4 .

The agent can choose between the following rewards

• R1=P(75,x)
• R2=P(70,x)
• R3=P(60,x)
• R4=P(50,x)

where

P is the recovery performance

𝑃 (b, p, c, s, g) = 𝑃(𝑏, 𝒙) =  0.01((𝑏 + 0.7𝑔) + 𝑘1 (0.35𝑐) − 𝑘2 (0.35𝑐 + 0.53𝑝) − 𝑘3 (0.35𝑐)2 (1)

and

Q a measure for grindability performance

k4 − 𝑔
Q=  (2)
k4

12
The reward is then given by

𝑅 = P  + 𝑄 + 0.01 ∗ 𝑁(0,1) (3)

where N(0,1) is the normal distribution with a mean of 0 and a standard deviation of 1 and ki are constants.

The selected policies in this case are:

• A random selection of the process

• A contextual epsilon-greedy strategy (C-Epsilon-greedy)
• A kNN-UCB strategy (kNN-UCB)

2.2.4 Case 4 Copper ore – Several non-steady state processing options, non-targeted feed grade
This case is treated as a contextual variant of case 2, by adapting algorithms from the previous subsection
and using a sliding time window. This represents seasonal variations in the performance in flotation. The
models are modified accordingly

• R1=Alternate with a period of 1000 time steps between P(75,x) and P(70,x)
• R2=Alternate with a period of 100 time steps between P(60,x) and P(80,x)
• R3=Alternate with a period of 150 time steps between P(55,x) and P(65,x)
• R4=Alternate with a period of 150 time steps between P(55,x) and P(70,x)

The selected policies in this case are

• A random selection of the process

• A sliding-window contextual epsilon-greedy strategy (SW-C-Epsilon-greedy)
• A sliding-window KNN-UCB strategy (SW-KNN-UCB)

3 Numerical experimental results from case study evaluations

For each case, numerical experiments have been conducted to compare the different strategies or
policies. The cumulative regret is used as the main metric since it is representative of the optimality of a
given policy. The process models were created in HSC Chemistry and calibrated with plant surveys. The
policies and plots were implemented in MATLAB 2018b. Each experiment consisted of 50000 blocks and
was run 150 times due to the stochastic nature of the policies and intervals of confidence were calculated
from the repeats. These are not to be confused with upper confidence bounds that are calculated inside
the UCB algorithm during each round.

The cumulative regret is defined as

𝑛𝑏 𝑜𝑓 𝑏𝑙𝑜𝑐𝑘𝑠
𝑅(𝑡) = ∑𝑡=1 (µ∗ − 𝜇𝑎 ) (4)

Where 𝜇∗ is the optimal process at time t and 𝜇𝑎 is the process selected by the policy. This value
measures how close one policy is to always make the best choice. For a perfect strategy, the cumulative
regret will always be zero. For any other policy the goal is to keep it small. From the plant point of view
this represents how close to the optimal operating point the decisions are for each ore type, at each time
step. The value of the regret can also be used to benchmark decisions at each point in time as it allows to
quantify optimality.

13
These algorithms can run on a recent laptop. As an example, on a laptop, 50 000 blocks were run 150
times in about 25 minutes to produce Figure 7. In practice, a decision is taken in less than a second.

3.1 Case 1
In this case, the benefits of having a policy are clear. Random choices lead to the highest regret, followed
by the epsilon-greedy strategy. Policies based on confidence bounds perform best. Figure 5 depicts the
cumulative regret as a function of time for the selected strategies. For the KL-UCB strategy even the
interval of confidence over the 150 repeats is shown, thus illustrating that the optimality is given with a
high statistical reliability.

Figure 5 : Regrets for case 1

3.2 Case 2
In this case, each process presents temporal variations while the feed varies around a controlled average.
When the performances of a process are not known in advance, especially when the variation cannot be
predicted, computing the average will fail since the period of the variations is not known. While the period
can be estimated, the presence of noise can make the estimation difficult if the amplitude is small.
Moreover, any other change to the plant (maintenance or automatic control) could influence the
estimates. Using sliding windows reduces the difficulty but still requires the window’s size to be calibrated.

14
Figure 6 shows the regret over time for non-stationary processes. While a random choice is still the worst
strategy, the epsilon-greedy strategies seem to fail at adapting to changing rewards using the same time
window as UCB strategies.

Figure 6 : Regrets for Case 2

3.3 Case 3
Once contextual data is introduced, the optimal decision depends on the block. As seen before, a purely
random strategy will fail here as well, but the difference between epsilon-greedy strategies and UCB is
much larger. This result could be explained by the fact that kNN-UCB takes the context into account by
taking similar decisions for close ore types. Epsilon-greedy strategies rely on the empirical estimate of the
mean only, which makes it difficult to distinguish between the impact of a different ore type and variations
in the process. In practice blending can offer a solution and reduce the problem to cases 1 or 2 in which
the grade of the feed is targeted. However, it is not always possible due to mining or equipment
constraints.

15
 Figure 7 : Regrets for case 3

To have a better view at the difference in behavior rather than the global one, it is useful to switch to a
logarithmic y axis to see the order of magnitude of the improvement by using UCB, see Figure 8.

Figure 8 : Regrets for case 3 with a logarithmic y axis

16
3.4 Case 4
This case is the most realistic one in which both the ore type and the process change in time. The proposed
solution is to combine kNN-UCB with a sliding time window. The choice of the window will depend on
previous experience of the process. This case is not relevant when simulations are used since calibrating
the simulations would mean knowing the period of the variations. Despite the non-stationary situation,
SW-KNN-UCB outperforms the other policies as well.

Figure 9 : Regrets for case 4

The application of bandit algorithms seems to provide a solution for complex decision-making in
geometallurgy. These methods do not only offer guarantees in terms of regret, but they also provide
suggestions for action. Regret is useful to quantify optimality but for each case, at each step a decision
with the probable highest reward is given. One can visualize the decisions in two different ways, either as
a table that gives the best decision at each time step as displayed on Table 5 or as a plot showing the
distribution of the choices for each strategy as seen on Figure 10.

17
 Table 5 : Decisions suggested by KNN-UCB for Case 3

Time Choose P1 Choose P2 Choose P3 Choose P4

1 0 1 0 0
2 0 0 1 0
3 0 0 0 1
4 1 0 0 0
5 1 0 0 0
6 1 0 0 0
7 1 0 0 0

Figure 10 : Choices distribution for case 3 by KNN-UCB, Random and Epsilon-greedy strategies

A major advantage of UCB methods is that they provide confidence bounds that are adjusted at every
time step. With time, the bounds become narrower which results in higher certainty about the choices
proposed by the method. On Figure 11 are shown box plots of the upper confidence bounds calculated
by kNN-UCB in case 3.

18
Figure 11 : Upper confidence bounds of the reward for kNN-UCB for 500 timesteps. When time increases the values shift from 0
to tighter bounds.

Despite an overall good performance, the hypotheses on the models are strong. In the non-contextual
case, theoretical results show that UCB strategies offer an optimal behavior but only in the case of
exponential family distributions. As such, other parametric families or even non-parametric models offer
no strong theoretical guarantees in terms of cumulative regret. For relatively simple cases, a good model
of the process can be obtained but the general case is missing. The reduction from a full flowsheet model
to a Bernoulli distribution is made simple by the existence of scenarios in most simulation software.

The contextual case is delicate due to the hypotheses, especially the Lipschitz continuity assumption. This
means that similar geometallurgical ore types will behave similarly in the process. This is justified if these
ore types produce the same particle size distribution, the same liberation spectrum and the same particles
geometry. If a sufficient number of tests have been done and if the mineralogy is well known, it can be
assumed but it remains a hypothesis.

A weaker hypothesis could be to use Hölder continuity instead which allows not so smooth function for
rewards, but no guarantee has been presented in literature in that case. This would allow slightly different
populations of particles, some noise in the context data or irregular process performance. Empirical
results suggest that it could be enough or at least that the performance is like the Lipschitz case. In
practice, a tool to characterize the regularity of a signal (for example the recovery or throughput of a
concentrate in the plant) is based on wavelet leaders (Serrano and Figliola, 2009) and provides an estimate
of both the global and local Hölder exponent. Incidentally, it offers a tool to quantify the stability of a
process performance with time by estimating its local regularity.

19
The existence of a lower dimensional manifold that supports the contextual data is another challenge. In
general, the problem of finding the right representation for high-dimensional data is hard but in the case
of ore physical properties, if the assumption of some form of geological continuity is accepted, then the
existence of such a manifold is implied. In any case, one can plot the representation and check whether
the assumption holds (see for instance Figure 4).

4 Conclusion
Bandit problems formulation for decision support in geometallurgy is a promising approach due to several
advantages. The main one is that a perfect knowledge of the process is not needed since in the stochastic
setting, bandit algorithms perform well in face of uncertainty. Therefore, it can be used either on
calibrated models with uncertainties, or directly on the process if a reward function can be obtain solely
based on the measured performances. An additional advantage for contextual methods is the fact that
each decision can be explained and motivated based on observations and available context. Experts could
then determine whether their experience gives similar policies.

The reward function also offers flexibility since one can construct a weighted sum of several parameters
including environmental scores and energy consumption for example (this was already present in case 3
and 4 in which a grindability index was included in the calculation of reward functions). If several ore types
are available and the problem is to blend the ores to obtain the best performances, one can use bandit
algorithms to optimally choose between blends. In that case the different processes (or arms) will be the
same flowsheet with the same operating parameters but different performances as a function of the ore
type. This case can be further adapted to time-varying blends that depend on mine production. From an
environmental point of view, if a life-cycle assessment score is calculated (for example the CO2 emitted
for each ton of concentrate) it can be included in the reward functions and provide guidance when
selecting a process. However, these methods do not calculate the weights assigned to each component
of a reward function. The construction of the reward functions is critical since it will determine which
choice is optimal and it is a human-based decision.

In practice, be it for full scale tests or simulated tests, the selection of the method is the same and is
illustrated in Figure 12.

20
 Sequential
Problem allocation
problem

Ore class Targeted Variable

Process Steady-state Dynamic Steady-state Dynamic

Algorithm UCB SW-UCB kNN-UCB SW-kNN-UCB

Figure 12 : Selection of the right method depending on the conditions

The integration of such methods on a mining or processing site is straightforward if element grades and
throughput are measured in the plant since it can be packed as an independent software package. As for
now, these methods are not included in mineral circuit simulation software to the best knowledge of the
authors but again in the case of HSC Chemistry, the data structure is appropriate and could be
implemented. One limitation lies in the connection with the mining scheduling software since it is not
always possible to have a choice in terms of blending or which material to feed at a certain point but even
in the case of a single process it is a valuable tool to evaluate decisions.

Limitations on theoretical guarantees on optimality can be an obstacle to the adoption of bandit

algorithms for geometallurgy. However, the results presented in this study suggest that it always provides
a better policy than pure random decisions (guided by other factors) or epsilon-greedy strategies that are
comparable to human-based decisions (estimates of the mean for each arm and a low probability of
exploration). One solution could be to use these policies and compare them against current policies over
a limited time horizon. Cumulative regret offers a quantitative assessment of policies.

As a conclusion, the proposed numerical methods provide good strategies for sequential resource
allocation in geometallurgy and tools to evaluate these policies over time.

5 Acknowledgements
The authors wish to thank LKAB, Boliden Mineral AB and Outotec (Sweden) AB for their support in this
study, Erdogan Kol for the interesting discussions about process simulation and flexibility. This work is
part of PREP project funded by VINNOVA within the Strategic innovation program for the Swedish mining
and metal producing industry (SIP-STRIM).

21
6 References
Agrawal, R., 1995. Sample mean based index policies by O(log n) regret for the multi-armed bandit problem.
Advances in Applied Probability 27, 1054–1078.

Andrusiewicz, M., Edraki, M., Tungpalan, K., Manlapig, E., Wightman, E., Keeney, L., Andrusiewicz, M., Keeney, L.,
Wightman, E., Edraki, M., 2014. An integrated approach of predicting metallurgical performance relating to
variability in deposit characteristics. Minerals Engineering 71, 49–54. doi:10.1016/j.mineng.2014.10.004

Auer, P., Cesa-Bianchi, N., Fischer, P., 2002. Finite-time Analysis of the Multiarmed Bandit Problem. Machine
Learning 47, 235–256. doi:10.1109/CEC.2016.7744106

Bubeck, S, Cesa-Bianchi, N., Bubeck, Sébastien, Cesa-Bianchi, N.O., 2012. Regret Analysis of Stochastic and
Nonstochastic Multi-armed Bandit Problems. Foundations and Trends in Machine Learning 5, 1–122.
doi:10.1561/2200000024

Burnetas, A.N., Katehakis, M.N., 1997. Optimal adaptive policies for Markov decision processes. Mathematics of
Operations Research 22, 222–255. doi:10.1287/moor.22.1.222

Garivier, A., Cappé, O., 2011. The KL-UCB Algorithm for Bounded Stochastic Bandits and Beyond, in: Sham Kakade,
U. von L. (Ed.), 24th Annual Conference on Learning Theory. JMLR : Workshop and Conference Proceedings.
pp. 359–376.

Garivier, A., Moulines, E., 2008. On Upper-Confidence Bound Policies for Non-Stationary Bandit Problems.

Hartland, C., Gelly, S., Baskiotis, N., Teytaud, O., Sebag, M., 2006. Multi-armed bandit, dynamic environments and
meta-bandits. Environments 1–14.

Kaufmann, E., Korda, N., Munos, R., 2012. Thompson Sampling: An Asymptotically Optimal Finite-Time Analysis BT
- Algorithmic Learning Theory, in: Bshouty, N.H., Stoltz, G., Vayatis, N., Zeugmann, T. (Eds.), . Springer Berlin
Heidelberg, Berlin, Heidelberg, pp. 199–213.

Keeney, L., 2010. The Development of a Novel Method for Integrating Geometallurgical Mapping and Orebody
Modelling.

King, R.P., 2012. Modeling and simulation of mineral processing systems. Elsevier.

Kłopotek, M.A., 2017. Machine Learning Friendly Set Version of Johnson-Lindenstrauss Lemma.

Lai, T.L., Robbins, H., 1985. Asymptotically Efficient Adaptive Allocation Rules. Advances in Applied Mathematics 6,
4–22.

Lin, I.J., 1989. The effect of seasonal variations in temperature on the performance of mineral processing plants.
Minerals Engineering 2, 47–54. doi:10.1016/0892-6875(89)90064-2

Lishchuk, V., Lamberg, P., Lund, C., 2015. Classification of geometallurgical programs based on approach and
purpose. 13th Biennial SGA Meeting. Mineral resources in a sustainable world. Volume 4. 1431–1434.

Lishchuk, V., Lund, C., Koch, P.-H., Gustafsson, M., Pålsson, B.I.I., 2018. Geometallurgical characterization of
Leveäniemi iron ore - Unlocking the patterns. Minerals Engineering 131, 325–335.
doi:10.1016/j.mineng.2018.11.034

Lund, C., Lamberg, P., 2014. Geometallurgy–A tool for better resource efficiency. European Geologist 37, 39–43.

Mathworks, 2016. MATLAB.

Outotec, 2016. HSC Sim.

22
Reeve, H.W.J., Mellor, J., Brown, G., 2018. The k -Nearest Neighbour UCB Algorithm for Multi-Armed Bandits with
Covariates, in: JMLR: Workshop and Conference Proceedings. pp. 1–29.

Serrano, E., Figliola, A., 2009. Wavelet Leaders: A new method to estimate the multifractal singularity spectra.
Physica A 388, 2793–2805. doi:10.1016/j.physa.2009.03.043

Thompson, W.R., 1933. On the likelihood that one unknown probability exceeds another in view of the evidence of
two samples. Biometrika 25, 285–294.

Tiu, G., 2017. Master thesis : Classification of Drill Core Textures for Process Simulation in Geometallurgy: Aitik Mine,
New Boliden. Luleå University of Technology, Luleå, Sweden.

Vos, K., 2018. Master thesis: Medium term planning with Evolution. Luleå Technical University.

23
Department of Civil, Environmental & Natural Resources Engineering
Division of Minerals and Metallurgical Engineering

ISSN 1402-1544
ISBN 978-91-7790-346-8 (print)
ISBN 978-91-7790-347-5 (pdf)

Luleå University of Technology 2019