EMA2 – Engineering Modelling and Analysis II Assignment 2 – Linear systems
Assignment 2 – Linear systems in geostatistics – due Fri 1 Sept 2017, 5PM Designed by Exequiel Sepulveda, Michael Leonard and Dmitri Kavetski
1. BACKGROUND
Geostatistical modelling is a very general class of statistical techniques for analysing spatially correlated data. Spatially correlated data refers to data that exhibits spatial patterns. For example, rain does not just drizzle uniformly over the Earth, but occurs in “clumps” (storms) that, despite being highly variable in time and space, are nonetheless fairly localized.
In Earth Sciences, geostatistical methods are widely used to estimate the spatial distribution of resources of interest. For instance, in forestry, we can use a few measurements to estimate the amount of forest biomass over large areas – because the amount of biomass at nearby locations will generally be similar (though not exactly the same). Similarly, in mining, we might be interested in estimating the location and amount of mineral resources, again based on the assumption that the amounts of mineral resources at nearby locations are correlated (“similar”). For example, given a set of copper grade (concentration) samples at a few locations over an ore field, we can (approximately) estimate the distribution of copper grade throughout the entire field.
2. KRIGING
The most common geostatistical estimation method is called “kriging” (named after Danie Krige, a South African Mining Engineer who pioneered the field of geostatistics in the 1950s). Kriging assumes that the correlation between the properties of any two locations depends only on their distance from each other. Despite this assumption, kriging has a number of favourable statistical properties. In particular it is an “unbiased linear estimator with minimum variance” – meaning that it is generally accurate, does not suffer from systematic errors (biases), yet is computationally fast.
Consider a spatially distributed variable U(v), where U is the variable of interest (e.g., copper grade) and v is a location in space, e.g., v = {x, y, z} in 3D space. In kriging, the value U pred at
location v0 is expressed as a combination of observed values at other locations. More specifically: N
0U
iUi
i1
Upred(v) Uobs (1)
where U obs U obs ( v ) is the value of U measured at location i (with coordinates v ), N is the total iii
number of measurements, U is the average of U obs over all measured locations, and the λ’s are the “kriging” coefficients. The measured values are often referred to as “observations”, “samples”, etc.
The kriging coefficients λ { , ,…, }T are defined by the solution of the “kriging” system 12N
Aλ b (2) where A is the (symmetric) covariance matrix of U across the sampled locations, and b is the
covariance vector between U at the location of interest and U at the sampled locations. Page 1 of 7
EMA2 – Engineering Modelling and Analysis II
Assignment 2 – Linear systems
The terms in the linear system in eqn (2) are given as follows:
a1,N cov(U1,U1) cov(U1,U2 ) a cov(U,U) cov(U,U)
a1,1 a1,2 a a
cov(U1,UN ) cov(U,U )
A2,1 2,2
a a N,1 N,2
2,N 2 1 2 2 a cov(U ,U) cov(U ,U)
2 N (3) cov(U ,U )
N,N N 1 N 2 b cov(U,U)
N N
101 bb2 cov(U0,U2)
(4)
b cov(U ,U )
N0N
The covariance matrix summarizes the statistical relationship between the values of U at different
locations. Roughly speaking, the covariance between the values of U at locations i and j can be interpreted as follows: large positive values of cov(Ui ,U j ) indicate that Ui U j ; conversely, near-
zero values indicate no systematic dependence between Ui and U j .
The estimation of the covariance matrix is in itself a formidable task. In kriging, we make the
simplifying assumption that the covariance cov(Ui,Uj) depends solely on the distance between
locations i and j, and does not depend on the individual locations themselves. On average, this assumption is not unreasonable – as discussed earlier, we expect nearby locations to have similar characteristics, and we expect distant locations to have different characteristics.
The assumption that the covariance depends solely on the distance can be expressed as
cov(Ui ,U j ) L[vi , v j ] (5)
where L[vi , v j ] is the distance between locations i and j.
The function (L) gives the covariance as a function of distance L, with a common choice being
nugget sill when L 0
(L)sillexp 3L whenL0 (6)
range
where sill, nugget and range are spatial covariance parameters that must be given in advance. Note that the strict test for equality in eqn (6) is usually replaced by a suitably selected tolerance.
Page 2 of 7
EMA2 – Engineering Modelling and Analysis II Assignment 2 – Linear systems
You should carefully study equations (2)-(6), and convince yourself that they provide the necessary machinery to evaluate the kriging predictor in eqn (1).
Given the general kriging framework in eqns (1)-(6), different “flavours” of kriging are available. In this assignment, we will only consider “global” kriging, where all available data are used for the estimation of properties at any location of interest.
Algorithm 1 – Global kriging
INPUT: Dataset of measurements {Uobs,Uobs,…,Uobs} at corresponding locations {v ,v ,…,v } ;
12N 12N values of spatial covariance parameters (sill, nugget and range).
OUTPUT: Estimated value of U at a location of interest.
Step 1. Specify a location of interest, v0;
Step 2. Construct the linear system in eqn (2):
use eqns (5) and (6) to evaluate the LHS matrix A and RHS vector b;
Step 3. Solve the linear system in eqn (2) for λ,
e.g., using methods from the Week 3 lectures;
Step 4. Evaluate the kriging estimator U pred (v ) in eqn (1) using the ’s computed in Step 3. 0
In this case, we need to construct and solve the full N N linear system. Unfortunately, for large N, solving such large linear systems can be computationally prohibitive (why?).
Page 3 of 7
EMA2 – Engineering Modelling and Analysis II Assignment 2 – Linear systems
4. TASKS
You are provided with a dataset of copper grade values (expressed as a percentage, %) sampled at N = 90 irregularly spaced 3D locations (file “copper[samples].txt”), as shown in Figure 1.
140 120 100
80 60 40 20
0 350
2.5
2
1.5
400
0.5
1
600
500
300 250
400
300
Y
200 150
200
100 100 50
00X
Figure 1. Copper sample locations and values (colour-coded according to grade, %).
You are asked to estimate the copper grade on a 10 10 grid along a horizontal plane at 75 m
elevation (Fig 2). The Nreq = 100 required locations are listed in the file “copper[reqLoc].txt”.
Figure 2. 2D grid in the horizontal plane of interest.
The covariance parameters in eqn (6) have already been estimated by the previous contractor (now
rich and retired): sill = 0.35 (%2), nugget = 0.05 (%2) and range = 80 m. Page 4 of 7
Z
EMA2 – Engineering Modelling and Analysis II Assignment 2 – Linear systems
Like a good engineer, kicking and screaming, you agree to undertake the work in several stages. All work done should be succinctly described in a report (4-page max length).
Stage 1 – Preparatory work
1.1. Write the equation for calculating the average U of a set of U values
1.2. Write the equation for calculating the distance between a pair of points i and j in 3D space
1.3. Prove analytically that the covariance matrix A is symmetric (Hint: consider eqns (3) and (5), and compare the calculation of ai,j and aj,i)
1.4. Comment on the relative advantages and disadvantages of global kriging. Consider in particular the question of computational cost in relation to the size of the kriging system, the number of sampled locations, and the number of locations where kriging estimates are needed.
1.5. Consider the algorithm for global kriging (Algorithm 1), as listed above.
(a) Familiarize yourself with the algorithm and its sequence of computational steps. Make sure you understand the sizes of the matrices and vectors, how the various terms fit together, etc;
(b) Implement the algorithm in Excel and test it using the first N = 10 measurements in the file “copper[samples].txt” to predict U at the first 3 locations in the file “copper[reqLoc].txt”. This small test case will help you understand the algorithm without having to simultaneously design a general Fortran code, or having to manipulate gargantuan matrices;
1.6. Implement the algorithm in Fortran and ensure it reproduces the Excel calculations for the test case in Task 1.5b. Your code should meet the requirements detailed in Appendix A;
1.7. Comment on the advantages and disadvantages of implementing the kriging algorithm in Excel vs Fortran. Consider aspects such as clarity, convenience, efficiency, etc.
1.8. Which linear solver algorithm considered in EMA2 would be particularly well suited for solving the linear systems of equations arising in global kriging? Briefly explain why?
Stage 2 – Practical application
2.1. Estimate the copper grade at all requested locations using the complete dataset.
2.2. Produce a contour plot of the gridded values (e.g., using Excel or Matlab – see Appendix B). Briefly comment on the results – where should we start digging?
Page 5 of 7
EMA2 – Engineering Modelling and Analysis II Assignment 2 – Linear systems
Appendix A – Kriging code specifications
The code implementing kriging should conform to the following requirements:
(a) Read the input files listed above to load all required data. The kriging parameters should be specified as PARAMETERs in your code. No inputs should be read from screen/keyboard.
(b) Be written in Fortran in good programming style (indentation, comments, modularity) (c) Write a “results file” named “copper[results].txt” with the following format:
Line 1: single header line (containing a descriptive comment)
Line 2: number of locations contained in the file (i.e., number of requested locations)
Line 3: INPUTS_OK or INPUTS_INVALID depending on whether all inputs are valid
Line 4: KRIGING_OK or KRIGING_FAILED depending on whether kriging system(s) solved ok Line 5: Column headings for table to follow below
Lines 6 onward: indexed 3D location coordinates and kriging estimates (1 location per line)
A template for the results file is provided at MyUni (your numerical values will be different).
(d) Write intermediate calculations to a “diagnostics file” named “copper[diagnostics].txt”:
(1 line): single header line (containing a descriptive comment)
(1 line): number of samples used in the analysis (i.e., the value of N)
(1 line): value of U
(1 line): label “Matrix of distances”
(N lines): values of the distances of all points from each other, as a matrix (1 line): label “Covariance matrix A”
(N lines): values of the covariance matrix A, as a matrix
(1 line): number of requested locations (i.e., the value of Nreq)
(1 line): index of the location being reported (i.e., 1, 2, 3, … , up to Nreq)
(1 line): label “Covariance vector b”
(1 line): values of the covariance vector b, as a row vector
(1 line): KRIGING_OK or KRIGING_FAILED depending on success of kriging solution (1 line): label “Kriging coefficients Lambda”
(1 line): values of solution λ, as a row vector (1 line): label “RHS error”
(1 line): values of εb = b – Aλ, as a row vector
(this allows you to check how well your computed λ satisfies eqn (2)) Note: write these location-specific intermediate results for all requested locations
A template for the diagnostics file is provided at MyUni (your numerical values will be different).
(e) Able to detect and report obvious input errors (e.g., invalid number of locations) and potential runtime errors (e.g., inability to solve the kriging system), as described above.
(f) Use your own judgement for any code/procedure aspects that have not been specified above.
IMPORTANT: The ability to reproduce a requested file layout is a basic programming skill. Your output files will be regenerated on the robomarker’s supercomputers and loaded according to the format listed above, so follow it carefully. The right number in the wrong place is a wrong result!
Page 6 of 7
EMA2 – Engineering Modelling and Analysis II Assignment 2 – Linear systems
Appendix B – Producing Contour Plots
Contour plots are a very useful way of visualizing 3D data. For example, contour plots of mineral ore amounts across a region can summarize areas of higher extraction value.
Microsoft Excel
Microsoft Excel provides some rather basic facilities for producing contour and surface plots. These are certainly not too flashy compared with plots produced using more specialised graphics software, but are nevertheless quite useful and easy to produce.
The following steps describe how to generate a contour plot in Excel.
1. Highlight the 2D cell range you would like to produce the contour plot for;
2. From menu, select Insert > Chart > Chart type “Surface” > select “Contour plot”.
For further customization, and for other plots available in Excel, consult the Excel Help. Note that different versions of Excel might have different menu names / locations for these options.
Matlab
Matlab can produce nice-looking contour plots with a wide range of customization options.
The basic syntax is contour(U), which draws a contour plot of matrix U, where U is interpreted as a
matrix of heights with respect to the x-y plane.
You can also use contourf(U), which will produce a filled-in contour plot.
The matrix of values to be plotted can be loaded from a file (e.g., named ‘data.txt’), which can be accomplished using the command U=importdata(‘C:/path/data.txt’)– note that you may need to specify the path to the file. The data file should only contain comma-delimited numerical values in matrix format. You can also import data from files with a more complex layout, but would need to learn more complex commands and settings (see Note 3 below).
Note 1: You are not expected to learn or use Matlab in this course, but it is a useful piece of software, especially if you enjoy programming and would like to broaden your skills.
Note 2: The machines in the ECMS computer labs have Matlab installed.
Note 3: Online help on Matlab is available at http://www.mathworks.com.au/help/matlab/.
Page 7 of 7