2019 | 72 / Special Issue | 111–126 | 19 Figs. | 8 Tabs. |  www.geologia-croatica.hr
 Journal of the Croatian Geological Survey 
and the Croatian Geological Society 

1 INTRODUCTION
Ore grade estimation is an important aspect of mine value evalu-
ation. The inverse distance weighted (IDW) and ordinary Krig-
ing (OK) are the most commonly used ore grade interpolation 
methods (NOVAK ZELENIKA et al., 2017). The IDW method is 
a widely used simple and convenient ore grade estimation 
method. The smoothness of IDW interpolation results is directly 
affected by the number of samples included in the estimation. The 
OK method is unbiased and has evolved into a series of estima-
tion methods including linear and nonlinear processes. The OK 
method is one of the most important Kriging methods, giving re-
sults with high smoothness (WANG et al., 2016; KUMARI & 
SINGH, 2018). In addition, OK uses a variogram to estimate the 
weight of the samples, but the variogram can only use the infor-
mation between the two points. Multi-point geostatistics uses a 
training image instead of a variogram, which considers informa-
tion from multiple samples, which is advantageous compared to 
the OK method. The single normal equation simulation (SNE-
SIM) is the most classic multi-point geostatistics (MPG) method. 
Scholars developed a series of improvements based on the SNE-
SIM method (ZHANG et al., 2015; HANSEN & LE, 2016). The 
simulation of patterns (SIMPAT) method compares the similarity 
between data events and data patterns in the training image and 
replaces data events with the most similar pattern (BURC AR-
PAT & CAERS, 2007; YANG, 2014). The SIMPAT method used 
the similarity calculation to compare the similarity between the 
data events and the data patterns in the training image and re-
places the data events with the most similar pattern. The SIMPAT 
method weakened the data stability requirement compared to the 
SNESIM method, which had some advantages; however, some 
errors in the data pattern selection and target continuity still ex-
ist. For example, sampling results from edge probabilities in-
crease the estimation uncertainty by not considering data fully. 
A multi-point geostatistical simulation based on the reservoir 
skeleton (SMPS) used the shape, type, and scale of the target body 

 
Multi-point geostatistics for ore grade estimation

Yu-Chen Song1, Zhan-Ning Liu1, Hai-Dong Meng1 and Xiao-Yan Yu1

1  Inner Mongolia University of Science and Technology, Institute of Mining Engineering, Baotou 014010, P.R. China;  
(songyuchen@imust.edu.cn; 494570031@qq.com; haidongm@imust.edu.cn; 421351047@qq.com) 

doi: 10.4154/gc.2019.23  

Abstract
A multi-point geostatistical method for ore grade estimation is introduced in order to fully utilize 
existing sampling information. A block model is used to construct a new three-dimensional train-
ing image instead of a variogram. Data events and pattern matching is improved, and the direc-
tionality of the data template is considered in the matching. The inverse distance weighted 
method is used to make up for the lack of multi-point geostatistics. The research improves the 
reliability of multi-point geostatistical estimation. Optimal estimation results for Li2O and Ta2O5 
come from the inverse distance weighted, ordinary Kriging, and multi-point geostatistical me-
thods. Multi-point geostatistical estimation results are compared with those of the inverse distance 
weighted and ordinary Kriging methods. Deviation, trend, and variogram analyses are used to 
assess the effect of multipoint geostatistical estimation. This study shows that reducing the sam-
ples participating in the estimation can reduce the maximum and minimum deviation of the es-
timated grade to a certain extent. The grade distribution pattern is the primary factor affecting 
minimum and maximum deviation. This study proves the reliability and accuracy of the multi-
point geostatistical method for ore grade estimation.

contained in the skeleton model to constrain and guide the multi-
point prediction results of the target simulation (YIN et al., 
2008a, 2008b) and was mainly used for river phase mode simu-
lation. Since the skeletal model for the target body was predicted 
before the simulation, simulation uncertainty can be reduced, and 
the continuity of the target body can be improved. The FILTER-
SIM method was used to remedy the long processing time of the 
SNESIM method (ZHANG et al., 2006). The FILTERSIM 
method determined the appropriate classification and dimensiona-
lity reduction of data events in the image while maintaining es-
timation accuracy. The dimensionality reduction process was im-
plemented using a mode filter combined with the FILTERSIM 
method (YU et al., 2016). Mode filters selected for a specific line ar 
combination of mode cell values include directional average, 
 gradient, and curvature. The improved method was primarily 
used for reservoir phase pattern simulation (WU & FU, 2016; 
STRAUBHAAR et al., 2016). Simulations were mostly pattern 
simulations (VRIES et al., 2009; HUYSMANS et al., 2014), 
comprehensive simulations (LOU et al., 2015; GARDET et al., 
2016), parameter optimization (GAO et al. 2016), and computa-
tional efficiency research (ZUO et al., 2016; ABDOLLAHIFARD 
& NASIR, 2017). There are few studies on the patterns of con-
tinuous variables, and even fewer have been conducted on the 
calculation of resource reserves using the MPG method.

Here, the SNESIM method is adapted and used for ore grade 
estimation. Multi-point geostatistics are compared with statistics 
from the OK and the IDW methods to confirm the validity of the 
results. Deviation analysis, trend analysis, and variogram analy-
sis are used to examine the estimation characteristics.

2 RESEARCH METHODS 
This section introduces the IDW, OK, and MPG methods as well 
as an improved MPG method, which overcomes the uncertainty 
in matching data templates to data events by incorporating the 
IDW method.

Article history:

Manuscript received March 06, 2019 
Revised manuscript accepted July 19, 2019 
Available online December 20, 2019

 
 
 
 
Keywords: Training image, Grade, Multi-point 
geostatistics, Variogram



G
eo

lo
gi

a 
C

ro
at

ic
a

Geologia Croatica 72 / Special Issue112

2.1. Inverse distance weighted (IDW) estimation
The inverse distance weighted (IDW) method is a geometric 
space interpolation method (spatial statistics) and is one of the 
most commonly used spatial interpolation methods. It uses the 
inverse relationship of the distance index powers of known neigh-
bour values to estimate grid point values and can be expressed 
by the equation (LI et al., 2008):

 

1

1

1

1
=

=

=
∑

∑

n

iP
ii

n

P
ii

Z
d

Z

d

     (1)

where Z is the estimation value for the point to be estimated. Zi 
is the i-th (i = 1, 2, ..., n) sample participating in the estimation. 
di is the distance between the point to be estimated and the sam-
ple. P is the power of the distance. The value of P significantly 
affects the interpolation result, and its selection criterion is the 
minimum average absolute error.

2.2. Kriging estimation
It is important to understand the theoretical explanation of the 
ordinary Kriging method, its experimental variogram, and theo-
retical variogram.

2.2.1. Variogram
The variogram is an important tool in the OK method and can be 
used alone or as the basis for the Kriging method to study the 
spatial variability of ore grades. The experimental variogram can 
be written as:

 
[ ]21( , ) ( ) ( )

2
γ = − +x h E Z x Z x h

 
    (2)

where γ(x, h) is the variogram value. h is the lag distance or 
distance between samples. Z(x) and Z(x+h) are the observed 
values of the regionalized variable Z at positions x and x+h, 
respectively.

The experimental variogram provides the spatial variation 
of the sample grade. However, an experimental variogram can-
not be directly used in Kriging estimation (GEIGER et al., 2016). 
It is necessary to fit the experimental variogram and obtain a 
theo retical curve after fitting, which is the theoretical variogram. 
The theoretical variogram can be used in Kriging estimation, and 
theoretical variograms used in this study are based on spherical, 
power function, and fractal models.

1. Spherical model

 
3

0 3

0

0 0

3 1( ) 0
2 2

C

γ

=


 = + ⋅ − ⋅ < ≤     
 + >

h

h hh C C a h
a a

C h a

     (3)

whereγ(h) is the value of experimental variogram. C0 is the nug-
get value. h is the lag distance. C is the arch height. C0+C is the 
sill value. a is the range.

2. Power function model

 ( ) , 0 2θγ θ= < <h Ah      (4)

where θ is a power exponent. A is a constant. Other parameters 
are the same as Eq. (3).

3. Fractal model
The variogram of the fractal Brownian motion is expressed 

as:

 
2( )γ = H

Hh V h      (5)

where: VH is a constant. h is the lag distance. H is the discontinu-
ity index (Hurst index).

2.2.2. Ordinary Kriging (OK) estimation
If x be any point to be evaluated in the study area, then the ore 
grade at point x can be expressed as:

 

*

1

( ) ( )λ
=

=∑
n

V i i
i

Z x Z x
 

    (6)

where Z*
V(x) is an estimation value of the point to be estimated. 

xi is the i-th most recent sample around the point to be estimated. 
Z(xi) is the grade of the i-th most recent sample within the range 
to be estimated. λi is the weight.

In Eq. (6), the estimation depends on the weight coefficient 
λi, while solving for the optimal λi value becomes essential to es-
timation. The best estimation is considered to be unbiased and 
has the least variance (optimal condition) in the OK. If the 
E[Z(x)]=m of the ore grade Z(x) is known, the calculation is called 
a simple Kriging estimation. If E[Z(x)] of Z(xi) is unknown, the 
calculation is called ordinary Kriging (OK).

2.3. Multi-point geostatistical estimation
2.3.1. Principle of the SNESIM
The classic SNESIM method takes multi-point geostatistics from 
theory to practicality, which is of great significance (STREBELLE, 
2017). The SNESIM uses data events to estimate the local proba-
bility of the point to be estimated. For two-point statistics, the 
probability that the variable type is sk state is written as:

 [ ]( , ) ( ; ) ( ; )Φ = × +h k E I u k I u h k      (7)

where Φ(h, k) represents the probability that variable types are sk 
states at locations u and u+h. sk represents the possible k states of 
the attribute S. k=1,…, K. h is the distance between two points. 
I(u; k) represents the probability when the state at position u is sk. 
I(u+h; k) represents the probability when the state at position u+h 
is sk.

The sample tn is determined by n vectors (hα, α=1, ..., n). The 
probability that n data points (u+h1, ..., u+hn) are simultaneously 
sk can be expressed by the following equation:

 1

( , ) ( ; )α
α=

 
Φ = + 

  
∏

n

h k E I u h k
 

    (8)

where Φ(h, k) represents the probability that multiple point attri-
butes are sk.

(u+h1, ..., u+hn) take the values of the probability of sk1, ...,skn 
respectively, which is written as:

 
1 1

1

( , , , , , ) ( ; )α α
α=

 
Φ = + 

  
∏ 

n

n nh h k k E I u h k
 

    (9)

where Φ(h1,…,hn, k1,…,kn) is the probability that multiple points 
have different attribute states.

By analyzing the definition of a data event, the probability 
of the data event dn is the probability of the above probability 
event, which is written as:

 
1 1

( )
( , , , , , )Φ = 

n
n n

n

c d
h h k k

N  
    

 
(10)



G
eologia C

roatica
Song et al.: Multi-point geostatistics for ore grade estimation 113

where c(dn) is the number of times the data event dn appears in 
the training image. Nn is the number of data events contained in 
the training image.

Assuming that the data event constitutes the indicator varia-
ble, D is written as:

 

1 ( ) , 1, ,
0

α α α= ∀ =
= 


kS u s n
D

Other
      (11)

where S(uα)=skα, which indicates a case where the attribute S is 
the state k at location uα.

The value of the estimation center point Ak can be written as:

 

1 ( )
0

=
= 


k
k

S u s
A

Other
     (12)

where Ak is an indicator variable associated with the occurrence 
of state sk at location u.

Using the simple Kriging method, the probability that any 
point will have a given attribute type is written as:

 [ ]( )1| 1 1λ= = = + −      k kProb A D E A E D       (13)

where Ak=1 is an indicator value corresponding to datum value 
sk at location u. D=1 is the observed data event. E[Ak] is the 
probability of the (unknown) state at location u to be sk. E[D] is 

the probability of the conditioning data event occurring. λ is the 
single weight.

The weight λ can be determined by: 

[ ]Var ,λ =   kD Cov A D , and 

[ ], = −          k k kCov A D E A D E A E D , 

[ ] [ ] [ ]( )Var 1= −D E D E D .

Thus there is:

 

[ ]
[ ] [ ]( )1

λ
−      =
−

k kE A D E A E D

E D E D
     (14)

Then there is:

[ ]( )
[ ]

[ ]

[ ] [ ]

1| 1 1

1, 1
1

= = = + −      

−      = +  

= =      = =
=

k k

k k
k

k k

Prob A D E A E D

E A D E A E D
E A

E D

E A D Prob A D
E D Prob D  

(15)

With the single normal equation described in Eq. (15), the 
conditional probability of the point attribute state sk can be cal-
culated. The attribute state value is obtained by the quotient of 

Figure 1. The SNESIM Execution Process.



G
eo

lo
gi

a 
C

ro
at

ic
a

Geologia Croatica 72 / Special Issue114

two multi-point probabilities, which are crucial to obtaining the 
conditional probability. Therefore, the probability of Prob[D=1] 
can be obtained by calculating the number of repetitions of D in 
the training image.

[ ] ( )1 11 , , , , , ( )= = Φ = n n n nProb D h h k k c d N is simi-
lar to the derivation of D probability. Prob[Ak=1, D=1] can also 
be obtained by calculating the repetitions of Ak=1 and D=1 in the 
training image: 

 

( )
1, 1= = ≈  

k n
k

n

c d
Prob A D

N  
    (16)

where ck(dn) represents the number of times the data event is at 
the centre position, whose attribute is sk, and is repeated in the 
training image. From the calculation of c(dn) and ck(dn), the proba-
bility can be obtained from the search tree. Therefore, the MPG 
method can be used to estimate the point value.

The conditional probability distribution of the data event at 
the point to be estimated is extracted from the search tree before 
the estimation is performed at location u, and the local condi-
tional probability is calculated for this location. The conditional 
probability distribution function at the point to be estimated is 
obtained. If the conditional probability at the point to be esti-
mated is 0, the global probability is extracted, instead of the local 
conditional probability. If the data event at the point to be esti-
mated rarely appears in the training image, the data point (sample 
grade point) furthest from the point to be estimated is removed, 
the number of data events repetitions at the point to be evaluated 
is increased, and the reliability of the local probability estimation 
is improved.

2.3.2. SNESIM improvement
In order to achieve accurate ore grade estimation, sample grade 
values need to be discretized in to multiple intervals. The number 
of discrete intervals influences estimation accuracy. The training 
image replaces the variogram in the multi-point geostatistical 
analysis. Training images used in the precious stratified phase 
model simulation are not suitable for ore grade estimation. There-
fore, it is important to build a new training image for ore grade 
estimation.

1. Sample processing
Ore grade value is a continuous variable and can be discretized 
using the following expression:

 

min

max min

 −
= × − 

k
p pp rp N

p p  
    (17)

where p is an arbitrary grade value. pmax is the maximum value 
of the samples. pmin is the minimum value of the samples. rp in-
dicates an integer up. N is the number of discrete intervals. pk is 
the discretized grade value.

In the estimation, the attribute (S) is divided into k states, 
where k is equivalent to the number (N) of discrete intervals. The 
estimation result of continuous variables is still a continuous vari-
able. The probability that the attribute to be evaluated is Sk is 
ck(dn)/c(dn) in MPG. The authors use the combined weighting 
valu es of S1, S2, ..., Sk, ..., Sn, and perform inverse discretized pro-
cessing on the value of the estimated attribute with probability 
values used as the weights. This process can be expressed as:

 
 

( )

( ) ( )
k

1 max min
max min min

c

2
=

×
−

= × − + −
× ×

∑
N

n k
k

E
n

d p
p p

p p p p
c d N N

(18)

where pE is an estimation grade. N is the number of possible 
 values of the estimation point, which is the same as the number 
of discrete intervals. pk is the grade value that may be obtained.

In theory, more discrete intervals lead to more accurate es-
timates. In fact, too many discrete intervals will reduce the 
matching rate between data events and data patterns, which re-
duces estimation accuracy. Multiple discrete interval schemes are 
used to evaluate the estimated effects and select the optimal es-
timates in subsequent studies. 

2. Directional problems and improvements in data event 
matching

The data template shown in Fig. 2. is used in multi-point geosta-
tistical estimation. According to the SNESIM method, the esti-
mation works by matching the six samples closest to the point to 
be estimated with six ore grade data (data pattern) obtained from 
scanning the training image; however, u1, u2, ..., u6 are directional 
in the data template (Fig. 2), and it is not appropriate to only con-
sider the 6 closest points to match the data pattern. Therefore, the 
authors consider the directional characteristics of the data event 
and samples. The estimation searches for more surrounding data, 
matches the data event direction with the template direction, and 
discards the samples with poor matching in some directions. The 
data template uses six data points corresponding to the six direc-
tions, which are the positive and negative directions of X, Y, and 
Z, respectively. (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), and 
(0, 0, -1) indicate the 6 directions in the template. The vector hα 
(α=1, ..., n) is used to indicate the directionality of the data. The 
six data points from the n dataset that best match the six direc-
tions in the template are selected for the estimation. Matching a 
single data with one direction of the template is written as:

 

[ ]
( )( )

( ) ( )

3

1

3 32 2

1 1

, =

= =

− −

=

− −

∑

∑ ∑

j j
j

j j
j j

h h l l

corr h l

h h h l

 

    (19)

where corr[h,l] represents the relationship coefficient between 
data events and the data template in a certain direction. hj and lj 
indicate the three-dimensional coordinates of the data points in 
the data event and data template, respectively. h and l  are the 
average values of hj and lj.

3. Data template and training image
Previous studies concerning the influence of the data template 
size on estimation results have shown that large data templates 
do not significantly contribute to improving the estimation accu-
racy (HU & CHUGUNOVA, 2008). Therefore, a small data tem-
plate is used in this study to save computing time (Fig. 2). The 
data template is used to scan the training image of porcelain clay 
to obtain the corresponding data pattern and store it in the search 
tree. According to the SNESIM method, the estimation principle 
is shown in Fig. 1. Local conditional probabilities are obtained 
by matching data events and data patterns. Monte Carlo sampling 
is then used to determine the point value to be estimated. Points 
for which the conditional probability distribution cannot be ob-
tained by matching, are estimated from the global probability us-
ing the Monte Carlo sampling method. In practice, the number 
of data events that do not get a conditional probability by match-
ing with the model is relatively small. In this case, using the 
global probability instead of the local conditional probability is 



G
eologia C

roatica
Song et al.: Multi-point geostatistics for ore grade estimation 115

likely to introduce errors into the estimation; therefore, the study 
uses the IDW method to estimate points not matched by the 
model and compensate for the uncertainty caused by global 
proba bility sampling. This method can also be regarded as a hy-
brid estimation method (M-I) of the MPG and the IDW.

The authors use the M-I, IDW, and OK methods to estimate 
the Li2O and Ta2O5 grades in porcelain clay and compare the 
estimation results with the sample grade by comparing deviation, 
trend, and variogram analyses. In the estimation process, the es-
tablished 3-dimensional ore body model is used, and the block 
unit (block model) is divided based on the ore body model to store 
the estimated ore grade information. The block model and drill 
sample information are exported and saved in .txt format. MAT-
LAB is used to read block information, implement the above es-
timates, and output the results.

3. STUDY AREA
An introduction to the regional geology, deposit geology, and ore 
grade characteristics of the porcelain ore mining areas is pre-
sented. 

3.1. Regional Geology of the Study Area
The study area is located in the transition zone between the Yang-
tze and South China fold systems. The Jiuling uplift comprises 
the northern part of the study area, the central part is the Pingle 
depression and the south is the Wugong uplift. The uplift and de-
pression contain two regional fault zones. The Yifeng-Jingdezhen 
and Pingxiang-Guangfeng faults control the distribution charac-
teristics of the stratigraphic facies, magmatic rocks, and minerals 
in the area.

3.1.1. Stratum
The Jiu-Ling uplift, located within the Jiu-Gang uplift, is a set of 
large-scale oceanic turbidite deposits more than 2,000 metres 
thick composed of volcanic sediments. The Sinian system in the 
Wugong Uplift is a series of siliceous iron-bearing rocks with 
low-grade metamorphosed sandstone and siltstone with mud-
stone and conglomerate. The thickness of the Sinian system is 
2000-10000 metres. Located on both sides of the Pingle depres-
sion and the uplift area, the Devonian-Triassic system is com-
posed of a set of marine carbonates and terrigenous clastic sedi-
mentary deposits, with a thickness of 3000-4000 metres. The 

fault basin next to the Shenda fault contains red clastic rocks and 
basalt with a thickness of 5000-7000 metres.

3.1.2. Construction
The study area is located in the southwestern Jiuling-Mekong-
shan thrust belt and the southern wing of the Jiuling compound 
anticline. It is located between the Yifeng-Jingdezhen deep fault 
zone a buried one and the Tongguxixiang-Fengxinshangfu fault.

The Jiuling uplift in the study area is composed of tight line ar 
and double inverted folds. Fold axes trend northeast accompanied 
by a large-scale ductile shear zone. The fault structure is developed, 
and there are three large-scale faults in the north-north-east, 
north-east, and north-west directions. Most of the faults are 
 pre-ore-formation structures, and the north-east and north-north-
east faults have associated mineralization. In the Yanshanian 
 period, rare metal ore-bearing hydrothermal fluids altered the 
alderite-type lithium-bearing porcelain ore body.

3.1.3. Magmatic rocks
From the Proterozoic to the Mesozoic, there has been abundant 
sub-magmatic activity, especially the Yanshanian magmatism 
associated with minerals containing W, Sn, and Rare Earth 
Elements (REEs). Due to the multi-stage nature of the magmatic 
activity, ore bodies commonly formed in complex rock masses 
or complex rock belts in the Jurassic or Early Cretaceous. The 
east channel of the Yifeng County-Gaoling-Baishuidong area is 
a complex rock belt dominated by Middle to Late Jurassic rocks. 
The lithology of the lithium-bearing porcelain ore body in this 
area is medium-fine-muddy-sodium feldspar-granite, and the 
genetic type is an altered alderite-type alkaline granite. In ad-
dition, a variety of associated rock drips and veins are developed 
in the area.

3.2. Porcelain clay area geology
The mining area, located in the Jiu-Gang uplift, has a complex 
lithology. Magmatic rocks are widely exposed in the mining area. 
Two lithium-bearing porcelain ore bodies exist in the mining 
area, named Li-1 and Li-2. Li-1 is the main ore body and is located 
in the eastern part of the mining area. Li-2 is a sub-mineral body 
located in the western part of the mining area. The ore bodies 
occur in the shape of irregular veins or long strips, and the overall 
trend is approximately north-south. The ore bodies trend to the 
south, and the dip angle is steep and irregular. The surface of the 
F3 fault in the mining area is about 150m at the surface of the 16-
24 line, and the drilling control tends to be about 180m long. West 
of the F3 fault, the surface of the 12-0 line shows a width of 
approximately 150m, and the western end extends for 1000m and 
becomes thinner in the vicinity of the 0 line. The 12-line control 
tends to be 100m long, and the 0-line control tends to be 30m 
long. The Li-2 ore body runs from the 13th to the 5th line in the 
southeast-northwest direction. The mineralization is uneven, 
showing irregular veins trending north, with an inclination of 
45~60° and an along- strike length of 340m. The northwest 
mineralization is wide, and it is cut off from the southeast by the 
5th line.

3.3. Three-dimensional geological model
3D geological models are an important way to visualize ore 
bodi es as well as an important process for future mine resource 
eva luation. Based on existing exploration data from the mining 
area, a 3-dimensional geological model and a geological database 
of the porcelain ore body are established (Fig. 3). The geological 

Figure 2. Data template.



G
eo

lo
gi

a 
C

ro
at

ic
a

Geologia Croatica 72 / Special Issue116

database is built with the help of Surpac software. Establishment 
of the geological database includes three main processes: data 
structure construction, exploration engineering data preprocess-
ing, and data importation. After the database is built, a 3D display 
of the borehole can be constructed. The ore body profile is deline-
ated on the profile of the exploration line formed by the borehole. 
A 3D ore body model is built through the developed three-dimen-
sional framework of the ore body. The geological database is 
mainly used to store the acquired exploration engineering data, 

and the sample grade combination can be performed on the stored 
exploration engineering data.

3.4. Porcelain clay samples
Firstly, the original Li2O and Ta2O5 sample grades were combined 
using the Surpac software sample combination tool. The 
composite sample is finally produced, and the Li2O and Ta2O5 
grades of porcelain clay samples are determined based on the 
geological database (Table 1). Fig. 4 and 5 show the Li2O and 
Ta2O5 grade histograms, respectively.

4. VARIOGRAM AND TRAINING IMAGES
Descriptions of the ore body test variogram calculation using the 
OK method follow. The fitting parameters of the spherical vario-
gram and fractal variogram are determined. The block model is 
used as the training image for grade estimation using the MPG 
method, and block models of different sizes are established to test 
the effect of block size on the estimation results.

4.1. Variogram fitting for porcelain clay
The spherical, power function, and fractal models were used to 
fit the experimental variogram for the sample grade of porcelain 
clay, and fitting parameters for the Li2O and Ta2O5 grades in the 
X, Y, and Z directions were obtained. In the spherical model fit-

Figure 3. Three-dimensional model of the porcelain clay ore body.

Figure 4. Li2O sample grade distribution.

Table 1. Li2O and Ta2O5 grade statistics results of porcelain clay.

Name Min Max Average Variance
Standard 
Deviation

Median Kurtosis Skewness
Coefficient of 

variation

Li2O (%) 0.2343 0.8679 0.51165 0.01193 0.10923 0.5183 -0.1208 -0.0729 0.21349

Ta2O5 (10-6) 15.436 99.40 38.935 65.5316 8.09516 38.6467 2.08147 15.5978 0.20791

Figure 5. Ta2O5 sample grade distribution.

Table 2. Fitting parametersforLi2O and Ta2O5 grade in the spherical model.

Name Direction C0 C C0+C a

Li2O

X 0.0045 0.0043 0.0088 37.15

Y 0.0045 0.0043 0.0088 37.15

Z 0.0029 0.0058 0.00087 67.61

Ta2O5

X 18.05 48.31 66.36 27.69

Y 18.05 48.31 66.36 27.69

Z 6.17 51.94 58.11 32.76



G
eologia C

roatica
Song et al.: Multi-point geostatistics for ore grade estimation 117

ting, the existing spherical model fitting tool in MATLAB was 
used. Table 2. shows the spherical model fitting parameters for 
the Li2O and Ta2O5 grades. Table 3. shows fitting parameters for 
the power function and fractal models.

The fractal and power function models were used for vario-
gram fitting. First, the fractal and power function models were 
converted into logarithmic form. Then, the SPSS tool was used 
to fit the experimental variogram. Table 3. gives the fitting parame-
ters of the fractal and power function models. The fractal model 
belongs to the variogram model without sill. According to the fit-
ting for Li2O grade in the fractal and power function models, the 
R-square is above 0.89 and the significance (Sig.) is below 0.05 
in the fractal model; hence, its fit is better than that the power ex-

ponential function model. Based on the fitting of Ta2O5 grade, 
the fit of the fractal model is better than the spherical model. The 
R-squares in the X direction and Y direction are above 0.78, and 
the significance is below 0.05. So the fitting effect is acceptable. 
The R-square of the fractal model is 0.91, and the significance is 
below 0.05; hence, the fitting effect is better than the spherical 
model.

4.2. Training image construction
The training image is a representation of the a priori geological 
model. Previous studies have indicated that a training image can 
be constructed by unconditional simulations, analog layer simu-
lations, or using digitized images sketched by geologists. It is dif-

Table 3. Fitting parameters forLi2O and Ta2O5grade in the power function model and the fractal model.

Name Direction Model
Model evaluation Model parameter

R-square F df1 df2 Sig. Constant b1 b2

Li2O

X 
power 0.777 62.888 1 18 0 -5.584 0.217 --

fractal 0.894 71.652 2 17 0 -6.049 0.623 -0.074

Y 
power 0.777 62.888 1 18 0 -5.584 0.217 --

fractal 0.894 71.652 2 17 0 -6.049 0.623 -0.074

Z 
power 0.899 160.768 1 18 0 -6.487 0.425 --

fractal 0.954 174.703 2 17 0 -7.065 0.929 -0.091

Ta2O5

X 
power 0.599 26.925 1 18 0 3.103 0.279 --

fractal 0.787 31.429 2 17 0 2.24 1.031 -0.137

Y 
power 0.599 26.925 1 18 0 3.103 0.279 --

fractal 0.787 31.429 2 17 0 2.24 1.031 -0.137

Z 
power 0.805 74.218 1 18 0 2.372 0.432 --

fractal 0.91 86.075 2 17 0 1.508 1.186 -0.137

Figure 6. Variogram fitting results of Li2O grade. A, B, and C are the fit using the spherical model in the X, Y, and Z directions, respectively. D, E, and F are the fit us-
ing the fractal model in the X, Y, and Z directions respectively.



G
eo

lo
gi

a 
C

ro
at

ic
a

Geologia Croatica 72 / Special Issue118

ficult to directly determine the grade distribution pattern of the 
ore body in grade estimation. It is important to establish a simple 
and effective training image for ore body grade estimation; there-
fore, the authors assign the sample point grade to the nearest 
block in the block model and build a 3-dimensional training 
ima ge based on the block model. The establishment of the block 

model is primarily accomplished by selecting the block range, 
determining the block size, adding block attributes, adding block 
constraints, and block assignment. The block range in this study 
should be slightly larger than the ore body range. A block model 
is constructed using a variety of block sizes (Table 4). Block as-
signments are primarily determined using the nearest distance 
method. Table 4. shows the block size parameters in five training 
images. This study uses multiple training images to obtain seve-
ral estimated grades and determines the optimal Li2O and Ta2O5 
grade estimation parameters.

The block size in Figs. 7. and 8. is 150m × 75m × 4m.
The length, width, and height of the training images of por-

celain clay are 1050m, 675m, and 224m, respectively. Figs. 8. and 
9. are training images of Li2O and Ta2O5, respectively.

Figure 7. Variogram fitting results of Ta2O5 grad. A, B, and C are the fit of the spherical model in the X, Y, and Z directions, respectively. D, E, and F are the fit of the 
fractal model in the X, Y, and Z directions, respectively.

Table 4. Block sizes of training images (pixel number).

Image name Image 1 Image 2 Image 3 Image 4 Image 5

X-direction (m) 75 75 75 38 38

Y-direction (m) 150 75 75 75 75

Z-direction (m) 4 4 2 4 2

Figure 8. Li2O training image.



G
eologia C

roatica
Song et al.: Multi-point geostatistics for ore grade estimation 119

5. ESTIMATION RESULTS AND ANALYSIS
The estimated features from the M-I method were analyzed, and 
optimal estimation results for Li2O and Ta2O5 determined. The 
estimation result of the M-I method was compared with the IDW 

and OK methods. The optimal estimation results of Li2O and 
Ta2O5 grades of porcelain clay using the IDW, OK, and M-I me-
thods were compared. The deviations, trends, and variogram ana-
lysis were examined in order to estimate ore grade quality.

Figure 9. Ta2O5 training image

Table 5. Estimated grade statistics for Li2O.

Training 
images

Discrete 
intervals

Min (%) Max (%) Average (%) variance
Standard 
deviation

median Kurtosis Skewness
Coefficient 
of variation

Image 1

4 0.277 0.824 0.480 0.012 0.112 0.472 -0.030 -0.909 0.233

6 0.277 0.824 0.483 0.011 0.104 0.498 -0.220 -0.451 0.215

8 0.274 0.824 0.485 0.012 0.107 0.512 -0.423 -0.464 0.221

10 0.266 0.824 0.487 0.010 0.101 0.484 0.037 -0.610 0.207

12 0.261 0.824 0.488 0.010 0.101 0.495 -0.106 -0.545 0.208

14 0.277 0.824 0.485 0.011 0.105 0.483 -0.137 -0.660 0.216

Image 2

4 0.277 0.798 0.477 0.012 0.108 0.472 -0.035 -0.766 0.228

6 0.277 0.824 0.486 0.012 0.110 0.498 -0.153 -0.437 0.226

8 0.274 0.824 0.488 0.012 0.110 0.512 -0.344 -0.390 0.225

10 0.266 0.824 0.497 0.011 0.103 0.519 -0.100 -0.462 0.206

12 0.274 0.824 0.489 0.010 0.102 0.517 -0.171 -0.595 0.208

14 0.277 0.824 0.490 0.011 0.106 0.498 -0.192 -0.625 0.216

Image 3

4 0.314 0.798 0.481 0.011 0.107 0.472 -0.084 -0.743 0.223

6 0.277 0.824 0.487 0.012 0.108 0.498 -0.244 -0.438 0.222

8 0.274 0.824 0.484 0.012 0.111 0.512 -0.241 -0.432 0.230

10 0.266 0.818 0.496 0.010 0.101 0.519 -0.145 -0.650 0.203

12 0.277 0.824 0.490 0.010 0.102 0.525 -0.279 -0.590 0.207

14 0.277 0.824 0.489 0.011 0.106 0.520 -0.273 -0.679 0.217

Image 4

4 0.277 0.798 0.483 0.012 0.112 0.472 -0.051 -0.845 0.231

6 0.277 0.824 0.482 0.013 0.113 0.498 -0.116 -0.614 0.234

8 0.274 0.828 0.490 0.013 0.114 0.512 -0.147 -0.082 0.233

10 0.266 0.836 0.497 0.010 0.101 0.519 -0.151 -0.440 0.204

12 0.277 0.842 0.489 0.010 0.101 0.518 -0.169 -0.493 0.206

14 0.277 0.824 0.489 0.011 0.105 0.503 -0.201 -0.591 0.215

Image 5

4 0.314 0.798 0.484 0.012 0.109 0.472 -0.091 -0.834 0.225

6 0.277 0.824 0.486 0.012 0.111 0.498 -0.240 -0.584 0.228

8 0.274 0.824 0.484 0.013 0.113 0.512 -0.247 -0.479 0.233

10 0.266 0.836 0.496 0.010 0.101 0.519 -0.123 -0.662 0.205

12 0.277 0.842 0.491 0.010 0.102 0.525 -0.286 -0.544 0.208

14 0.277 0.845 0.491 0.012 0.109 0.524 -0.143 -0.433 0.221



G
eo

lo
gi

a 
C

ro
at

ic
a

Geologia Croatica 72 / Special Issue120

5.1. Estimation results
Firstly, the estimation results of the M-I method using different 
training images were calculated, and the optimal estimation re-
sults determined. Similarly, the optimal results of Li2O and Ta2O5 
estimated using the IDW and OK methods are given.

5.1.1. M-I estimation results
Li2O and Ta2O5 grades in porcelain clay ore were evaluated using 
the M-I method. The estimation characteristics of the M-I method 
were analyzed using the deviation between valuation grade and 
sample grade.

1. Results and analysis of Li2O grade valuation
Table 5. shows the results of the Li2O grade estimation using the 
M-I method.

From the data event and data pattern matching rate in Fig. 
10, the number of discrete intervals increases, which reduces the 
successful matching probability. At this time, more data points 
use the overall probability instead of the local probability in the 
valuation. In the M-I method, the number of points to be esti-
mated using the IDW method increases as the number of discrete 
intervals increases, which has a gradual upward trend. In terms 
of influence of the number of pixels in the training image on 
matching rate, the high-pixel training image has a somewhat im-
proved matching.

Fig. 11A. shows the minimum deviation in Li2O valuation 
grade. In the M-I method, increasing number of pixels in the 
training image has no noticeable effect on the estimated mini-
mum deviation. For training image 1 with 12 discrete intervals, 
the smallest minimum deviation is 11.27%.

Fig. 11B. shows the maximum deviation in Li2O valuation 
grade. As the number of discrete intervals increases, the maxi-Figure 10. Data event and data pattern matching rates of Li2O.

Table 6. Estimated grade statistics for Ta2O5.

Training 
images

Discrete 
intervals

Min (10-6) Max (10-6) Average (10-6) variance
Standard 
deviation

median (10-6) Kurtosis Skewness
Coefficient of 

variation

Image 1

4 21.886 88.441 38.460 63.911 8.110 34.785 0.376 0.613 0.211

6 18.287 88.441 37.841 44.640 6.753 36.934 -0.136 3.058 0.178

8 18.145 88.441 38.277 34.898 8.305 38.009 -0.100 5.948 0.217

10 18.016 88.441 37.857 42.792 6.696 38.654 -0.094 3.501 0.177

12 17.586 88.441 38.490 50.829 7.199 39.084 -0.298 2.255 0.187

14 17.279 88.441 38.765 46.878 7.176 39.391 0.169 3.495 0.185

Image 2

4 21.886 88.441 37.787 54.518 7.505 34.785 0.401 1.184 0.199

6 18.287 88.441 37.803 48.845 7.027 36.934 -0.272 2.563 0.186

8 18.287 88.441 38.473 43.656 8.562 38.009 -0.023 3.132 0.223

10 18.016 88.441 37.854 50.432 7.139 38.654 -0.516 2.957 0.189

12 17.586 88.441 38.147 56.644 7.562 39.084 -0.605 2.416 0.198

14 17.279 88.441 39.070 46.273 7.139 39.391 -0.107 2.867 0.183

Image 3

4 19.167 88.441 37.543 54.591 7.396 34.785 0.402 1.270 0.197

6 18.287 88.441 37.734 44.287 6.695 36.934 -0.098 3.096 0.177

8 18.287 88.441 38.107 40.285 8.115 38.009 0.027 3.705 0.213

10 18.016 88.441 37.966 41.358 6.589 38.654 -0.118 3.410 0.174

12 17.586 88.441 38.150 44.446 6.667 39.084 -0.298 3.251 0.175

14 17.279 88.441 38.483 41.369 6.432 39.391 -0.084 3.854 0.167

Image 4

4 24.291 88.441 39.378 91.267 9.559 46.923 -0.448 -1.116 0.243

6 22.433 88.441 38.305 44.749 6.828 36.427 0.677 2.703 0.178

8 18.350 88.441 38.083 46.580 6.928 41.675 -0.093 2.193 0.182

10 18.287 88.441 38.140 45.505 6.818 36.427 -0.349 2.486 0.179

12 18.263 88.441 38.457 42.186 9.487 39.926 -0.072 3.140 0.247

14 18.287 88.441 38.510 44.251 6.818 36.427 0.026 2.582 0.177

Image 5

4 24.291 88.441 39.713 90.120 9.499 46.923 -0.497 -1.049 0.239

6 19.167 88.441 37.478 51.473 7.174 36.427 0.362 2.101 0.191

8 18.350 88.441 38.155 50.497 7.205 41.675 0.086 1.928 0.189

10 18.287 88.441 38.048 46.319 6.877 36.427 -0.259 2.543 0.181

12 18.263 88.441 38.096 42.926 8.272 39.926 0.029 2.997 0.217

14 18.287 88.441 38.278 42.450 6.682 36.496 -0.046 2.989 0.175



G
eologia C

roatica
Song et al.: Multi-point geostatistics for ore grade estimation 121

mum deviation decreases. When the resolution of the training 
image is high, the decrease in maximum deviation becomes more 
obvious, indicating that increasing the training image pixels can 
reduce the maximum deviation. When training image 5 is used 
and the number of discrete intervals is 14, the smallest maximum 
deviation is -2.61%.

Fig. 11C. shows the Li2O valuation grade average deviation, 
indicating that a larger number of training image pixels can re-
duce maximum deviation. As the number of discrete intervals 
increases, average deviation decreases first and then increases. 
An inflection point occurs when the number of discrete intervals 
is 10. As the matching ratio decreases (Fig.10), the number of dis-
crete intervals and estimated deviation will increase. When train-
ing image 4 is used and the number of discrete intervals is 10, the 
smallest average deviation is -2.88%.

The deviation of estimated Li2O grade standard deviation is 
estimated using the M-I method (Fig. 11D). On the one hand, the 
deviation of the standard deviation variation trend is affected by 
the number of discrete intervals; on the other hand, it is affected 
by the matching rate of data events to data patterns. When train-
ing image 5 is used and the number of discrete intervals is 14, the 
Li2O grade minimum deviation using the M-I method is -0.54%.

The four graphs A, B, C, and D are the estimated deviations 
of the minimum, maximum, average, and deviation of standard 
deviation of the estimated Li2O grade, respectively.

Comprehensive analysis of Fig. 11. shows that the M-I method 
can effectively improve Li2O grade estimation, and the number 
of discrete intervals can improve estimation accuracy. In addi-
tion, the excessive number of calculation intervals (10 or more) 
reduces the valuation effect. When the number of training image 
pixels is large, the success rate between the data event and data 
template can be improved, which helps to improve and stabilize 
the evaluation. The M-I method with training image 5 and 10 dis-
crete intervals has the best evaluation. It is estimated that the 
minimum deviation of the Li2O grade is 13.52%, the maximum 
deviation is -3.65%, the average deviation is -3.14%, and the 
deviation of the standard deviation is -7.10%.

2. Results and analysis of Ta2O5 grade valuation
Table 6. shows the Ta2O5 grade estimated using the M-I method.

Figure 12. shows the data event and data pattern matching 
rate for the M-I method. Increasing the number of training image 
pixels can improve the matching rate. The highest number of 
training image pixels has the highest matching rate, and the train-
ing image with the lowest number of pixels has a lower matching 
rate. Upgrading the training image in the Ta2O5 grade evaluation 
can help reduce the use of global probability sampling in the val-
uation process and thus improve valuation accuracy. Li2O grade 
estimation shows a similar relationship between matching ratio 
and number of training image pixels. The effect of the training 
image on matching efficiency is related to the sample distribution 
form or sampling method. The M-I method is affected by the 

Figure 11. Deviation in Li2O grade estimation based on integrated approach.



G
eo

lo
gi

a 
C

ro
at

ic
a

Geologia Croatica 72 / Special Issue122

sample grade distribution or sampling method when estimating 
the ore grade.

Figure 13. presents an analysis of the estimated Ta2O5 de-
viation when using the M-I method. The effect of training image 
pixels on the minimum deviation is not significant (Fig. 13A). 

After the number of discrete intervals increases, the minimum 
deviation decreases. When using training image 1, 2, or 3 with 
14 discrete intervals, the smallest minimum deviation is 11.94%.

Figure 13B shows the maximum deviation of the Ta2O5 va-
luation grade. The number of discrete intervals and the pixel size 
of the training image have no affect on the maximum deviation. 
The maximum deviation is consistently 11.03%. Increasing the 
number of training image pixels fails to reduce the average de-
viation (Fig. 13C), indicating that the valuation results are af-
fected by the IDW method when using the M-I method. The ave-
rage deviation of the M-I method is smaller than the MPG 
method, indicating that the M-I method is still superior to the 
simple MPG method. When using training image 2 with 14 dis-
crete intervals, the smallest minimum deviation is 0.35%.

Figure 13D. shows the deviation of the standard deviation of 
the Ta2O5 valuation grade. Under the M-I method, the deviation 
of standard deviation fluctuates greatly with the number of dis-
crete intervals. However, compared with the MPG method, the 
M-I method shows a significant reduction in deviation of the stan-
dard deviation.

This study concludes that the use of the M-I method, train-
ing image 4, and 14 discrete intervals provides the best results. It 
is estimated that the minimum deviation of the Ta2O5 grade is 
18.31%, the maximum deviation is -11.03%, the average deviation 
is -2.15%, and the deviation of the standard deviation is 2.19%. 
The Ta2O5 valuation shows that the M-I method is closer to the 

Figure 12. Data event and data pattern matching rates of Ta2O5.

Figure 13. Deviation in Ta2O5 grade estimation based on integrated approach.



G
eologia C

roatica
Song et al.: Multi-point geostatistics for ore grade estimation 123

sample Ta2O5 grade distribution than the MPG method, indicat-
ing the M-I method is superior to the MPG method. The M-I 
method is not scientific and rigorous enough to use global sam-
pling instead of overall sampling for data events that cannot be 
matched. Increasing the number of discrete intervals reduces the 
matching rate between data events and data templates, which in-
creases the probability of using global sampling to generate local 
samples. This increases the inaccuracy of the valuation. Increas-
ing the number of training image pixels can improve the match-
ing rate to a certain extent but cannot fundamentally solve the 
problems outlined above. The IDW is used instead of the global 
sample estimate for points that cannot be estimated using local 
sampling. Although the IDW method has a certain smoothing ef-
fect, its valuation results are better than estimates determined us-
ing global sampling. Because data events around the point to be 
estimated cannot be used, it is commonly the case that there is a 
small probability of data events; therefore, it is unscientific to use 
global sampling instead of local sampling.

The four graphs A, B, C, and D show the estimated devi-
ations of the minimum, maximum, average, and deviation of 
standard deviation of the estimated Ta2O5 grade, respectively.

5.1.2. Optimal estimation results
The optimal estimated grades of Li2O and Ta2O5 are shown in 
Tables 7 and 8, respectively. IDW1 and IDW2 respectively rep-
resent the optimal estimation results when the three samples and 
five samples of Li2O and Ta2O5 were included in the estimation 
of IDW. Optimal estimation results for Li2O and Ta2O5 grade oc-
curred when using the OK method and when the fractal model is 
used. When image 5 is used and the discrete interval is divided 
into 10 pieces, the Li2O grade has the best M-I estimation result. 
When image 4 is used and the discrete interval is divided into 14 
pieces, the Ta2O5 grade has the best estimation result.

5.2. Estimation results analysis
A comparison is presented of the estimation method characteris-
tics of the IDW1, IDW2, OK, and M-I methods. Deviation analy-
sis, trend analysis, and variogram analysis were performed on 
Li2O and Ta2O5 grades. Finally, the advantages of the M-I method 
were compared with IDW1, IDW2 and the OK methods.

5.2.1. Deviation analysis
The estimated deviation of the optimal estimated grade of porce-
lain clay was analyzed. The minimum deviation of Li2O esti-

mated grade ranges from 13.52% to 20.24% (Fig. 14). The maxi-
mum deviation ranges from 2.47% to 6.30%. The Li2O grade 
estimation minimum deviation is significantly greater than the 
maximum deviation. Estimation results of the IDW1, IDW2, OK, 

Table 8. Optimal estimation results of Ta2O5 grade.

Method
Min

 (10-6) 
Max

 (10-6) 
Average

 (10-6) 
Variance

Standard 
deviation

Median
 (10-6) 

Kurtosis Skewness
Coefficient
 of variation

IDW1 16.6269 92.9562 38.4829 44.4826 6.6695 38.8171 -0.0167 5.6815 0.1733 

IDW2 18.3676 84.7443 38.6357 38.5344 6.2076 38.7802 -0.2967 3.4585 0.1607 

OK 18.4970 88.6346 38.6976 36.4862 6.0404 38.7437 -0.2018 4.0988 0.1561 

M-I 18.2870 88.4411 38.5101 44.2511 6.8176 36.4270 0.0257 2.5824 0.1770 

Table 7. Optimal estimation results of Li2O grade.

Method
Min
 (%) 

Max
 (%) 

Average
 (%) 

Variance
Standard 
deviation

Median
 (%) 

Kurtosis Skewness
Coefficient 
of variation

IDW1 0.2704 0.8465 0.5000 0.0127 0.1129 0.5101 -0.2713 -0.3560 0.2257 

IDW2 0.2761 0.8181 0.5001 0.0116 0.1077 0.5224 -0.3394 -0.3547 0.2153 

OK 0.2817 0.8133 0.4898 0.0105 0.1026 0.4972 -0.2813 -0.5665 0.2094 

M-I 0.2660 0.8362 0.4956 0.0103 0.1015 0.5194 -0.1233 -0.6617 0.2048 

Figure 15. Ta2O5 grade estimation deviation.

Fig. 14 Li2O grade estimation deviation.



G
eo

lo
gi

a 
C

ro
at

ic
a

Geologia Croatica 72 / Special Issue124

and M-I methods are similar, indicating that the large minimum 
deviation is not due to a change in estimation method and the re-
sults of the M-I method are accurate. The average deviations of 
the IDW1, IDW2, OK, and M-I methods are -2.28%, 2.26%, 
-4.28%, and -3.14%, respectively. The deviation of the standard 
deviation of the estimation results of the 4 methods is small. The 
standard deviation of the IDW2 method is the smallest (1.31%).

Figure 15. shows the estimated deviation of the Ta2O5 grade. 
The minimum deviation of the estimated grade of Ta2O5 ranges 
from 7.71% to 19.83%. The maximum deviation ranges from 
6.48% to 14.74%. The average deviation ranges from 0.61% to 
1.12%. The Ta2O5 grades estimated by the 4 methods are ex-
tremely similar. Although samples used in the estimation are 
from the same exploration project, the deviations of the Li2O and 
Ta2O5 grades are different for the same estimation methods, in-

dicating that sample grade distribution has an obvious influence 
on the estimation results. In addition, the estimated deviations of 
the 4 methods are very similar for the same grade estimation. 
This further illustrates that the estimation results of the improved 
multipoint geostatistics (M-I) are reliable and accurate.

5.2.2. Trend analysis
Trend analysis was performed on the Li2O and Ta2O5 estimated 
grades for porcelain clay. Figure 16. shows the trend of the Li2O 
grade change with elevation. Overall, the estimated ore grade is 
lower than the actual ore grade. The Li2O grade at the shallow 
(large) elevation of the ore body is higher than the deep (small 
elevation) grade. The Li2O grades estimated by the IDW1, IDW2, 
OK, and M-I methods have greater consistency with the eleva-
tion trends, and the estimated grades for the different methods 
are similar. Figure 17. shows the trend of the Ta2O5 grade with 

Figure 16. Trend analysis of Li2O grade. Figure 17. Trend analysis of Ta2O5 grade.

Figure 18. Variogram analysis of Li2O grade. A, B, and C are the variograms of the estimated and sample grade in the X, Y, and Z directions, respectively.



G
eologia C

roatica
Song et al.: Multi-point geostatistics for ore grade estimation 125

elevation. The Ta2O5 and Li2O estimation trends are similar. The 
Li2O and Ta2O5 grade is greater at higher elevation (shallow) than 
at lower elevation (deep).

5.2.3. Variogram analysis
The variograms of the estimated grade and sample grade of Li2O 
and Ta2O5 are shown in Figures 18. and 19. The IDW1 method 
used 3 samples in the estimation process. The IDW2 and OK 
methods used 5 samples. The number of samples included in the 
estimation using the IDW method in the M-I method is 4. The 
largest difference between the Li2O estimated grade and sample 
grade variograms is in the X and Y directions for the IDW2 
method. The OK method has the second largest gap. The vario-
gram based on the IDW1 method is closer to the sample grade 
than the M-I variogram. That is to say that the IDW1 variogram 
is closest to the sample grade in the X and Y directions. In the Z 
direction, the difference between the Li2O OK variogram and the 
sample grade variogram is the largest. The second largest gap is 
in the IDW2 variogram. Similarly, the IDW1 Li2O grade esti-
mated is closest to the sample grade, followed by M-I. The Li2O 
variograms indicate that the number of samples participating in 
the estimation has important effect on the estimated grade vario-
gram.

Variograms of the estimated Ta2O5 grade differ greatly from 
the sample variogram (Fig. 19). The difference among the esti-
mated grade variograms is small. The OK and IDW2 variograms 
are similar, and the IDW1 and M-I variograms are also similar. 
The IDW1 and M-I variograms are slightly more similar than the 
OK and IDW2 variograms. Overall, the variograms of the esti-
mated grade are significantly lower than that of the samples. The 
IDW2 and OK estimated grades have a stronger smoothing effect. 
Although the Li2O and Ta2O5 samples are from the same explo-
ration project, the gap between the estimated variograms and 
sample variograms is different, indicating that the sample grade 
distribution has an effect on the morphology of the variogram.

Comprehensive analysis of the Li2O and Ta2O5 grade vario-
grams indicates that the greater the number of samples partici-
pating in the estimation, the larger the difference between the 
sample grade variogram and the estimated grade variogram will 
be. If the same number of samples participate in the estimation, 
the estimation result of the OK method has a higher smoothness 
than that of the IDW2 method. If the number of samples partici-
pating in the estimation increases, the variogram of the IDW es-
timated grade becomes smoother. The sample distribution sig-
nificantly inf luences the estimated grade variogram. The 
estimated grade for the M-I method is similar to that of the IDW 
and OK methods, indicating that the ore grades estimated by the 
multi-point geostatistical method are reliable and accurate.

6. DISCUSSION
Valuation using the M-I method indicates that improving the 
training image resolution can improve the MPG valuation to 
some extent. However, increasing the training image resolution 
seriously affects the computational efficiency and does not com-
pletely solve the problem of matching between data events and 
data patterns. Matching efficiency between the data events and 
data patterns is an important factor affecting the accuracy of 
multi-point geostatistical estimation. In theory, increasing the 
number of discrete intervals can improve valuation accuracy, but 
it will reduce the matching rate between data events and data pat-
terns, which affects valuation accuracy. The M-I method can 
overcome the issue of inaccurate low-frequency grade informa-
tion in MPG estimation. The study produced the optimal results 
for different training images and discrete interval numbers when 
estimating using the M-I method.

This study evaluated the grades of Li2O and Ta2O5 in porce-
lain clay using an improved MPG method. Results were com-
pared with results from the OK and IDW method, demonstrating 
the feasibility, reliability, and accuracy of the M-I method for ore 
grade estimation. The stability and accuracy of multi-point geo-

Figure 19. Variogram analysis of Ta2O5 grade. A, B, and C are the variograms of the estimated and sample grade in the X, Y, and Z directions, respectively.



G
eo

lo
gi

a 
C

ro
at

ic
a

Geologia Croatica 72 / Special Issue126

statistical estimation indicates that the new training image con-
struction method is feasible. Compared with the OK method, the 
use of training images instead of variograms improves the utili-
zation efficiency of sampling information. In other words, infor-
mational relationships between more samples are utilized. The 
improved method is theoretically superior to variograms that only 
utilize the pattern of information between two samples.

This study used deviation, trend, and variogram analyses to 
analyze the improved MPG method estimation features from dif-
ferent angles. It proves that the new training images improved 
data event matching methods, and the uncertainty problem pro-
cessing methods can meet the requirements of ore grade estima-
tion. An accurate estimate of the Li2O and Ta2O5 ore grades was 
achieved.

7. CONCLUSION
(1) This study introduced MPG into the field of ore grade estima-
tion, provides a new training construction method, and applies it 
to Li2O and Ta2O5 ore grade estimation in a porcelain clay de-
posit. A good estimation result is obtained using the M-I method.

(2) The matching of data events and data patterns was im-
proved. A new matching method was used to account for the di-
rectionality of the data template. The inverse distance weighted 
method was used to compensate for the lack of MPG. The im-
provements enhance the reliability of MPG estimates.

(3) This study produced the optimal estimation results for 
Li2O and Ta2O5 when the IDW, OK, and M-I methods were used. 
The grade distribution pattern is the main factor affecting the 
minimum and maximum deviation. If more samples are included 
in the estimation, the estimation results of the IDW and OK me-
thods will be smoother.

(4) Deviation, trend, and variogram analyses were used to 
study the estimated features of the inverse distance weighted, or-
dinary Kriging, and multi-point geostatistics. The reliability and 
accuracy of the M-I estimations is shown by comparing the 
estimations of the IDW and OK methods.

ACKNOWLEDGEMENT
The funding was provided by the China National Natural Science 
Foundation (NO. 71363040).

REFERENCES
ABDOLLAHIFARD, M.J. & NASIRI, B. (2017): Exploiting transformation-domain 

sparsity for fast query in multiple-point geostatistics.– Computational Geosciences, 
21/2, 1–11. doi: 10.1007/s10596-016-9612-1

BURC ARPAT, G. & CAERS, J. (2007): Conditional Simulation with Patterns.– Math. 
Geol., 39/02, 177–202.

GAO, S.C., TIAN, M. & SUN, Z. (2016): Multiple-point Geostatistics modeling factors 
analysis and optimization by using Markov Chain model.– Mathematics in Prac-
tice and Theory, 46/0, 202–211.

GARDET, C., RAVALEC, M.L. & GLOAGUEN, E. (2016): Pattern-based conditional 
simulation with a raster path: a few techniques to make it more efficient.– Stoch. 
Env. Res. Risk. A., 30/2, 429–446. doi: 10.1007/s00477-015-1207-1

GEIGER, J., JAKAB, N., CSÖKMEI, B., HORVÁTH, Z. & GELLÉRT, B. (2016): Sta-
tistical and geostatistical study of Rn and hydrocarbon components of a soil gas 
monitoring system: an application to surface hydrocarbon exploration.– Geologia 
Croatica, 69/2, 255–268. doi: 10.4154/gc.2016.21

HANSEN, T.M. & LE, T.V. (2016): Bach T. MPSLIB: A C++ class for sequential simu-
lation of multiple-point statistical models.– Soft. X., 5, 127–133.

HU, L. Y. & CHUGUNOVA, T. (2008): Multiple-point geostatistics for modeling sub-
surface heterogeneity: A comprehensive review.– Water. Resour. Res., 44, 1–14. 
doi: 10.1029/2008WR006993

HUYSMANS, M., ORBAN, P. & COCHET, E. (2014): Using Multiple-Point Geosta-
tistics for Tracer Test Modeling in a Clay-Drape Environment with Spatially Vari-
able Conductivity and Sorption Coefficient.– Math. Geosci., 46, 519–537. doi: 
10.1007/s11004-013-9502-1

KUMARI, M. & SINGH, C.K. (2018): GaRiRO: Gradient and residual integrated rank 
ordering of stations in rainfall monitoring network.– Earth. Sci. Inform., 11, 1–14. 
doi: 10.1007/s12145-018-0332-z

LI, Z.L., WANG, P. & ZHANG, X.L. (2008): Improvement and application of inverse 
distance weighted method.– Met. Min., 02, 88–92.

LUO, H.M., YANG, P.J. & WANG, C.J. (2015): Litho facies simulation based on mul-
tiple-point geostatistics multiple data joint constraints.– Oil Geophysical Prospect-
ing, 50/01, 162–169+19.

NOVAK ZELENIKA, K., VIDAČEK, R., ILIJAŠ, T. & PAVIĆ, P. (2017): Application 
of the deterministical and stochastical geostatistical methods in petrophysical 
mode ling, case study upper Pannonian reservoir in Sava depression.– Geologia 
Croatica, 70/2, 105–114. doi: 10.4154/gc.2017.10

STRAUBHAAR, S., RENARD, P. & MARIETHOZ, G. (2016): Conditioning Multiple-
point Statistics Simulations to Block Data.– Spatial Statistics, 16, 53–71. doi: 
10.1016/j.spasta.2016.02.005

STREBELLE, S. (2001): Sequential simulation drawing structures from training ima-
ges.– Palo Alto, PhD Thesis, Faculty of Science, Stanford University.

VRIES, L.M.D., CARRERA, J. & FALIVENE, O. (2009): Application of Multiple Point 
Geostatistics to Non-stationary Images.– Math. Geosci., 41/1, 29. doi: 10.1007/
s11004-008-9188-y

WANG, G.W., FENG, Y. & CARRANZA, M. (2016): Typomorphic characteristics of 
pyrite: Criteria for 3D exploration targeting in the Xishan gold deposit, China.– J. 
Goechem. Explor., 164, 136–163. doi: 10.1016/j.gexplo.2016.01.003

WU, T. & FU, B. (2016): Application of the multiple-point method to geological 
modeling in the Sulige gasfield.– Geol. Prospect., 52/05, 985–991.

YANG, P.J. (2014): Geostatistics inversion from two-point to multiple-point.– Chin. J. 
Geophys. Progr., 29/05, 2293–2300.

YIN, Y.S., WU, S.H. & ZHAI, R. (2008): The forecasting of fluvial facies reservoir dis-
tribution by using Simpat.– J. Southwest Petrol. Univ., 02, 19–22+184–185.

YIN, Y.S., WU, S.H. & ZHANG, C. (2008): M. Multi-point geostatistical method based 
on reservoir skeleton.– Sci. China Ser. D., S2, 157–164.

YU, S.Y., LI, S.H. & HE, Y.B. (2016): Multiple-point geostatistics algorithm bases on 
pattern scale-down cluster.– Acta. Petrol. Sinica, 37/11, 1403–1409.

ZHANG, T., BOMBARDE, S. & STREBELLE, S. B. (2006): 3D porosity modeling of 
a carbonate reservoir using continuous multiple point statistics simulation.– SPE 
Journal, 11/3, 375–379. doi: 10.2118/96308-PA

ZHANG, T., DU, Y. & HUSNG, T. (2015): A reconstruction method for spatial data us-
ing parallel SNESIM.– J. Comput. Dev., 52/06, 1431–1442. 

ZUO, C., TENG, Q.Z. & HE, X.H. (2016): A fast 3D reconstruction algorithm of mul-
tiple point statistics.– J. Sichuan Univ: Natl. Sci. Ed., 53/02, 337–346.

Author Contributions statement: Yu-Chen Song designed the experiment, Zhan-Ning Liu performed the experiment and wrote 
the paper, Xiao-Yan Yu processed the data and images, and Hai-Dong Meng revised the paper.

Competing interests statement: The authors declare that there are no competing interests between them, including financial 
and non-financial interests.