Fitting a linear model of coregionalization for soil properties using simulated annealing

16 Pages·2003·0.16 MB·English

ARTICLE IN PRESS
Geoderma1962(2003)1–16
www.elsevier.com/locate/geoderma
Fitting a linear model of coregionalization for soil
properties using simulated annealing
R.M. Larka,*, A. Papritzb
aMathematicsandDecisionSystemsGroup,SilsoeResearchInstitute,WrestPark,
Silsoe,BedfordMK454HS,UK
bETHZu¨rich,SoilPhysics,InstituteofTerrestrialEcology,Grabenstrasse3,CH-8952Schlieren,Switzerland
Received5February2002;accepted8January2003
Abstract
Itisnotsimpletomodelcrossandauto-variogramstodescribethecovariationoftwoormoresoil
properties,sincethemodelsthatarefittedmustmeetcertainconstraints.Theseconstraintsaremost
readily met byfitting alinear model of coregionalization (LMCR).
This presents practical problems. Not all combinations of authorized variogram functions
constitute a LMCR. This paper presents a method for automated fitting of variogram functions to
auto and cross-variogram estimates, subject to the constraints of the LMCR. The method uses
simulated annealing to minimize a weighted sum of squares between the observed and modelled
variograms.
Themethodwasappliedtosomedataonsoil.Itwasfoundtoberobusttotheinitialchoiceof
variogram parameters. Practical methods for setting up a good cooling schedule for the simulated
annealing arediscussed.
D2003Elsevier ScienceB.V.All rights reserved.
Keywords:Geostatistics;Cross-variogram;Cokriging;Modelfitting
1. Introduction
Multivariategeostatisticalmethodshavebeenusedtoanalyzethespatialcovariationof
soil properties (Goovaerts, 1997; Webster et al., 1994). As in univariate geostatistics, the
central assumption is that the observed value of a variable u, z , at location x is a
u
* Correspondingauthor.Fax:[phone].
E-mailaddress:[email protected](R.M.Lark).
0016-7061/03/$-seefrontmatterD2003ElsevierScienceB.V.Allrightsreserved.
doi:10.1016/S0016-7061(03)00065-X
ARTICLE IN PRESS
2 R.M.Lark,A.Papritz/Geoderma1962(2003)1–16
realization of a random function, Z (x), which meets the intrinsic hypothesis (Matheron,
u
1962), so that we may define the auto-variogram, c (h), as a function of a lag vector h
u,u
only:
1 h i
c ðhÞ¼ E fZ ðxÞ(cid:5)Z ðxþhÞg2 : ð1Þ
u;u 2 u u
In multivariate geostatistics, we make a similar assumption that the covariance of the
differences{Z (x)(cid:5)Z (x+h)and{Z (x)(cid:5)Z (x+h)forvariablesuandvdependsonlyon
u u v v
the lag vector h, so that we may define the cross-variogram:
1
c ðhÞ¼ E½fZ ðxÞ(cid:5)Z ðxþhÞgfZ ðxÞ(cid:5)Z ðxþhÞg(cid:9): ð2Þ
v;u 2 v v u u
Whenautoandcross-variogramfunctionshavebeenobtaineditispossibletoestimate
the values of variables u and v jointly at unsampled sites by cokriging (e.g. Goovaerts,
1997). This is particularly useful when one of the variables (perhaps a cheaply measured
property) is more densely sampled than the other. In this case, cokriging may give more
precise estimates of the less-densely sampled variable than could be made by ordinary
(univariate) kriging from the same data (e.g. Atkinson et al., 1992).
Cross and auto-variograms are usually obtained in two steps. First, point estimates are
obtained using the standard estimator:
1 NXv;uðhÞ
cˆ ðhÞ¼ fz ðxÞ(cid:5)z ðx þhÞgfz ðxÞ(cid:5)z ðx þhÞg; ð3Þ
v;u 2N ðhÞ v i v i u i u i
v;u i¼1
wherethereareN (h)pairsofobservationsseparatedbylagh.Commonly,weignorethe
v,u
directional component of the lag vector and use a lag distance h. Note that the cross-
variogram can only beestimated in this way from data at sites where both variables have
been observed. If too few (or no) sites have observations of both variables, then the
coregionalization may be modelled non-parametrically (Ku¨nsch et al., 1997) or via the
pseudocross-variogram(Papritzetal.,1993;PapritzandFlu¨hler,1994;Lark,2002a).This
is not considered further in this paper.
Thesecond stepistofitcontinuous functionstothepointestimatesofthevariograms.
These are needed for cokriging, or planning a grid sample (McBratney and Webster,
1983). Not any continuous function will serve, since the functions must define a positive
definite covariance structure (i.e. where all linear combinations of the variables have a
positive variance). Thisiseasilyachievedintheunivariate case byfittingan‘authorized’
model (Journel and Huijbregts, 1978), such as linear combinations of some or all of the
nugget,exponentialandsphericalfunctions(WebsterandOliver,2001).Inthemultivariate
case, it is necessary but not sufficient to use authorized variogram functions. A model of
coregionalization comprising authorized functions fitted independently to the cross and
auto-variograms is not necessarily positive definite.
The usual strategy to ensure a positive definite model is to use a linear model of
coregionalization (LMCR) (Journel and Huijbregts, 1978; Goovaerts, 1997). Here we
assume that all the N constituent random functions Z (x) are linear combinations of a
w u
ARTICLE IN PRESS
R.M.Lark,A.Papritz/Geoderma1962(2003)1–16 3
commonset ofindependent randomfunctions ofmean zero and variance 1,Yl(x)(lisan
k
index not a power). Thus
XL Xnl
Z ðxÞ¼l þ al YlðxÞ: ð4Þ
u u u;k k
l¼0 k¼1
The coefficients al are specific to the random functions Z (x) and Yl(x). As stated
u,k u k
above, the random functions Yl(x) are mutually independent. However, if two such
k
functions have a common index, l, then they have the same spatial correlation structure.
There are L+1 such structures with variogram functions g(h), l=0,1,...,L, where:
l
1E(cid:4)fYlðxÞ(cid:5)YlðxþhÞgfYlVðxÞ(cid:5)YlVðxþhÞg(cid:5)¼ gðhÞ;l ¼lV; k: ¼kV ð5Þ
2 k k kV kV l
0;otherwise
There are therefore n random functions that share the variogram function g(h),
l l
although they are orthogonal.
It can be shown that the cross-semivariances of any two of the N intercorrelated
w
randomvariables,Z (x)andZ (x),couldbeexpressedasalinearcombinationoftheL+1
u v
basic variogram functions:
L
X
c ðxÞ¼ bl gðhÞ; ð6Þ
u;v u;v l
l¼0
where,
Xnl
bl ¼ al al : ð7Þ
u;v u;k v;k
k¼1
Thus the cross-variogram matrix,
2 c ðhÞ c ðhÞ ::: c ðhÞ 3
1;1 1;2 1;Nw
6 7
666 c2;1ðhÞ c2;2ðhÞ ::: c2;NwðhÞ 777
6 7
6 ::: 7
6 (cid:10) (cid:10) (cid:10) 7
6 7
GðhÞu6 7;
6 ::: 7
6 (cid:10) (cid:10) (cid:10) 7
6 7
6 7
6 ::: 7
6 (cid:10) (cid:10) (cid:10) 7
6 7
4 5
c ðhÞ c ðvÞ ::: c ðhÞ
Nw;1 Nw;2 Nw;Nw
can be written as
L
X
GðhÞ¼ gðhÞS; ð8Þ
l l
l¼0
ARTICLE IN PRESS
4 R.M.Lark,A.Papritz/Geoderma1962(2003)1–16
where
2 bl bl ::: bl 3
1;1 1;2 1;Nw
6 7
666 bl2;1 bl2;2 ::: bl2;Nw 777
6 7
6 ::: 7
6 (cid:10) (cid:10) (cid:10) 7
6 7
Slu6 7:
6 ::: 7
6 (cid:10) (cid:10) (cid:10) 7
6 7
6 7
6 ::: 7
6 (cid:10) (cid:10) (cid:10) 7
6 7
4 5
bl bl ::: bl
Nw;1 Nw;2 Nw;Nw
IftheL+1matricesS,coregionalizationmatrices,areallpositivesemi-definiteandthe
l
L+1 variogram functions g(h) are all authorized, then the LMCR has a positive definite
l
covariance structure.
To fit a LMCR, we must find a set of variogram functions and corresponding
coregionalizationmatricesthatoptimizeameasureofthefitofG(h)tothepointestimates
ofthevariograms,subject totheconstraintsonthecoregionalizationmatrices.Thisisnot
trivial, particularly as the number of variables is increased. Goovaerts (1997) describes
how a LMCR may be fitted by visual assessment and trial and error. A semi-automated
procedurewasdevisedbyGoulardandVoltz(1992).Ifthedistanceparametersofthebasic
variogramfunctionsaregiven,thenthisalgorithmwillfindestimatesofthecorresponding
coregionalization matrices optimizing a measure of the goodness of fit and meeting the
constraints on the coregionalization matrices.
ItwouldbeidealtoautomatethefittingofaLMCRtothesamedegreethatthefittingof
theunivariatelinearmodelofregionalizationisautomated(McBratneyandWebster,1986).
Thatistosay,havingspecifiedavariogramfunctionorfunctionsallparameterestimates,
including the distance parameters, are optimized by a statistical criterion. The problem in
thecaseoftheLMCRisensuringthatthefittedcoregionalizationmatricesareallpositive
semi-definite.Ku¨nschetal.(1997)describeawayoffittingtheLMCRusinganon-linear
regressionmethod(amodifiedLevenberg–Marquardtscheme,Pressetal.,1992)tofitthe
distance parameter and Goulard and Voltz’s (1992) algorithm to fit the coregionalization
matrices.Thisisdoneiterativelyinalternatestepsuntilconvergence,judgedbytherelative
changeinthesum-of-squarescriterioninthetwofittingstepsorbyaceilingonthenumber
ofiterations.Inthispaper,weproposeanalternativewayoffittingtheLMCR,andcompare
it with the algorithm of Ku¨nschet al.(1997) using some dataon thesoil.
ThenumericalmethodthatweproposeinthispaperforfittingtheLMCRissimulated
annealing. This is a method of optimization where constraints can be very simply
incorporated. Geostatisticians have used simulated annealing for simulating regionalized
variables (Deutsch and Journel, 1992), for maximum likelihood estimation of variogram
models(Pardo-Iguzquiza,1997)andforoptimizingthespatialdistributionofsamplepoints
(vanGroenigen,1999;Lark,2002b).Thispapershowshowsimulatedannealingcanbeused
tofittheLMCR.
ARTICLE IN PRESS
R.M.Lark,A.Papritz/Geoderma1962(2003)1–16 5
2. Theory and the algorithm
2.1. Simulated annealing
A detailed account of simulated annealing is given by Aarts and Korst (1989), but the
principles are now outlined. Consider the objective function S that we want to minimize
withrespecttoparametersp ,p ,...,p.Insimulatedannealing,thisisdoneanalogouslyto
1 2 i
thephysical process of amoltenmetalsolidifying asit is cooled.In this physical system,
the objective function is the lattice energy of the solid and the parameters are the relative
positionsoftheconstituentatoms.Theobjectivefunctionisataglobalminimumwhenthe
atoms are arranged in a regular crystalline structure. This will be achieved if the metal
cools sufficiently slowly (annealing). Locally optimal aggregations of atoms can be
disrupted during annealing because the heat energy of the system allows it to ‘jump’
over energy barriers. If, by contrast, the metal is cooled quickly then it will form an
amorphous or partly crystalline structure (quenching), that is it will ‘stick’ at a local
minimum some distance from the global minimum.
In simulated annealing the parameters of S are perturbed randomly. Over a time step,
t!t+1, a single change in a parameter occurs, resulting in a change, S !S . If
t t+1
S VS then the change is accepted (minimization problems). If S >S then the
t+1 t t+1 t
decision whether or not to accept the change is made randomly with probability p of
a
acceptance, where
pa ¼eSt(cid:5)Sctþ1: ð9Þ
This is the Metropolis criterion, central to simulated annealing. The probability of
acceptanceofaparticularincreaseinSisdeterminedbytheparameterc.Intheannealing
analogue, c is the temperature of the system. Increasing the temperature makes the
acceptance of a particular transition more likely.
Aseriesofrandomadjustmentstoasetofparameters,acceptedorrejectedaccordingto
Eq.(9),withafixedvalueofc,constitutesaMarkovchain.AnewMarkovchainmaybe
initiated by ‘cooling’ the system (reducing c). In the physical analogue of the molten
metal, after cooling to a fixed temperature, the system will tend, over many transitions
governedbytheMetropoliscriterion,toathermalequilibrium,suchthattheprobabilityof
a particular energy state is governed by the Boltzmann distribution (see Aarts and Korst,
1989).EachMarkovchainofthesimulatedannealingwilltendsimilarlytoanequilibrium,
the time (number of transitions) that this takes depending (inversely) on temperature.
The goal of simulated annealing is to take the system through a series of (approx-
imately)equilibratingMarkovchainsatsuccessivelylowertemperatures,untilthesystem
has been cooled sufficiently that no further changes are accepted and it is at an optimum
which is close to or equal to the global optimum. Global optimality is not guaranteed for
chains of finite length.
In practical terms, a good solution requires:
1. That c is large enough to allow the system to escape from local minima early in the
simulation, and is not reduced too rapidly.
ARTICLE IN PRESS
6 R.M.Lark,A.Papritz/Geoderma1962(2003)1–16
2. That each Markov chain is long enough to come close to equilibrium.
3. That sufficient Markov chains are run to bring the system to its optimum.
These conditions define the ‘cooling schedule’ for a simulated annealing. Various
approachestoselectionofacoolingschedulehavebeenproposed,thesimplestbeingdue
to Kirkpatrick et al. (1983), that is followed here:
1. AninitialvalueofcisfoundsothatnearlyallchangesareacceptedinthefirstMarkov
chain of the simulation.
2. The value of c is reduced from one chain to the next by a control parameter a , with
c
0<a <1, so that c =a c . Typically 0.8Va V0.99 (Aarts and Korst, 1989).
c m+1 c m c
3. A fixed length of the Markov chain is commonly predetermined. This requires an
element of trial and error.
4. ThealgorithmisterminatedwhensomenumberofMarkovchainshavebeencompleted
without any further change in the objective function.
More complex cooling schedules have been proposed (e.g. Aarts and Van Laarhoven,
1985), but these still require the more or less arbitrary definition of parameters of the
schedule.Trialanderrorwithparametersofthescheduleisnecessary,butthiscanbedone
easily in the case of fitting a model like the LMCR since the model fit can be assessed
rapidly by visual inspection, and other diagnostics to evaluate a particular annealing (see
Results and discussion).
In the following section, an objective function is defined that may be minimized to fit
the LMCR.
2.2. The objective function
The objective function is that proposed by Goulard and Voltz (1992). We denote by
Gˆ(h)(thematrixofestimatesofthesemivariancesandcross-semivariancesforlagh),and
i i
by G(h) we denote the corresponding matrix of semivariances and cross-semivariances
i
obtainedfromafittedLMCRwithparametersthatwewishtoestimate.GoulardandVoltz
(1992) proposed that the parameters be found that minimize WSS, the Euclidean-like
distance criterion over m lags:
WSS¼Xm wðhÞtraceh(cid:14)V(cid:12)GˆðhÞ(cid:5)GðhÞ(cid:13)(cid:15)2i; ð10Þ
i i i
i¼1
where w(h) is a positive weight and Vis some positive definite matrix. Here we follow
i
Goulard and Voltz (1992) by choosing for Va diagonal matrix with the inverse sample
variances,s(cid:5)2,asthediagonalelements.Thispreventsthosevariableswithlargevariances
u,u
fromhavinganundueinfluenceonthefittingofthemodel.Weconsideredtwooptionsfor
the weights, w(h). The first was to set them all to 1. The second was to set w(h) to the
i i
average number of pairs contributing to the estimated semivariances and cross-semi-
variances at lag h.
i
ARTICLE IN PRESS
R.M.Lark,A.Papritz/Geoderma1962(2003)1–16 7
2.3. Constrained optimization
2.3.1. Initialization
Havingspecifiedasetofauthorizedvariogramfunctions,wemustfindparametersthat
will minimize WSS subject to the constraints that the coregionalization matrices are all
positive semi-definite. To do this, an initial guess at the model parameters is supplied,
along with the cooling schedule: c —the initial system temperature, a , N —the number
1 c m
of perturbations of each model parameter in a single Markov chain, and N—a threshold
t
numberofMarkovchainssothat,ifthereisnochangeinWSSoverN successivechains,
t
the algorithm will stop.
It is also necessary to specify the random process whereby each parameter, h will be
perturbed during the simulated annealing. In this algorithm, a parameter is perturbed by
adding a random variable of mean zero. The random variable has a uniform density
functionovertheinterval[(cid:5)d ,d ],sod ,themaximumperturbationforeach
h,max h,max h,max
parameter, must be given.
AninitialguessoftheLMCRparametersisrequired.Thiscouldbegeneratedinseveral
ways.
(i) With relatively few (two or three) variable guesses of the values of the parameters
may be made by eye, checking the coregionalization matrices in a spreadsheet.
(ii) The distance parameter could be guessed by eye, then coregionalization matrices
obtained using Goulard and Voltz’s (1992) algorithm (this procedure was used by
Goovaerts, 1997).
(iii) If the simulated annealing method is reasonably robust, then a very poor initial
approximation to the LMCR, perhaps with a standard set of positive semi-definite
coregionalization matrices, could be used to initiate the program. The robustness of
the procedure is investigated empirically in this paper.
2.3.2. The first Markov chain
The value of WSS for the initial parameter values is calculated. Each parameter of the
LMCR is then perturbed randomlyin turn. After one perturbation, we first check that the
new set of parameters is legal, i.e. the distance parameters of the variogram are non-
negativeandthecoregionalizationmatricesarepositivesemi-definite.Ifthisisnotthecase
then the previous values of the parameter are restored and a new random perturbation of
the same parameters are drawn and assessed. A perturbation rejected because of the
constraintsoftheLMCRisnotcountedasastepintheMarkovchainofrandomchanges.
After a legal perturbation the change in WSS is evaluated. If WSS is reduced or
unchangedthenthenewvalueoftheparameterisaccepted.Otherwise,arandomdecision
on acceptance or rejection is made according to the Metropolis criterion (Eq. (9)). This
procedure is iterated until N legal random changes of each model parameter have been
t
accepted or rejected. The first Markov chain is then complete. The proportion of legal
changes accepted by the Metropolis criterion is calculated. If this is somewhere in the
range 0.90–0.99, then the algorithm continues (see next section), otherwise it is
terminated and a new initial temperature is specified (a larger value if a larger proportion
of accepted changes is required).
ARTICLE IN PRESS
8 R.M.Lark,A.Papritz/Geoderma1962(2003)1–16
2.3.3. Markov chains 2, 3,...
The system is cooled by changing the initial temperature from c to c =a c . The
1 2 c 1
second Markov chain then proceeds as above. This procedure is then iterated. At the
end of each Markov chain, the new value of WSS is compared to that from the
previous N chains. When there has been no change in the WSS over this number of
t
chains, the algorithm stops.
Intheremainderofthispaper,wedescribetheuseofthisalgorithmtofitanLMCRto
variograms estimated from data on trace metal concentrations in the topsoil.
3. Materials and methods
3.1. Soil data
The soil data that we used in this study were measurements of concentrations of trace
metals in the topsoil of the Swiss Jura, collected and described by Atteia et al. (1994). In
particular,weusedthedataoncadmium,nickelandzincinGoovaerts’spredictionsubset
ofthese data whicharelisted inan appendix toGoovaerts (1997).Goovaerts(1997) uses
these particular data to illustrate how the LMCR may be fitted using Goulard and Voltz’s
(1992) algorithm to obtain coregionalization matrices and fitting the distance parameters
of the variogram by eye.
3.2. Analysis
Autoandcross-variogramsfortheselectedvariableswereestimatedusingtheGAMV2
subroutine in the GSLIB library (Deutsch and Journel, 1992). By inspection, it was
decided to fit a LMCR with three components, a nugget, and two spherical variograms.
Fourinitialsetsofparameterswereusedtoinitiatetheprogram(Tables1aand1b).These
were as follows:
1. Aguessobtainedvisually,usingaspreadsheettoensurethattheinitialvaluesconstitute
a positive definite LMCR.
2. The parameters obtained by Goovaerts (1997) using Goulard and Voltz’s (1992)
algorithm and fitting the distance parameters by eye.
3. Twosetsofparametersdescribingmodelsthatfitthedataverypoorly.Thesewereused
to test the robustness of the fitting procedure.
Simulatedannealingwasused,withthefollowingcoolingschedule,afterKirkpatricket
al. (1983); a =0.975, N =25, N =15, and d =1.0 for all parameters. The initial
c m t h,max
temperature, c , was selected by trial and error so that 99% of the perturbations were
1
accepted in the initial Markov chain.
The algorithm of Ku¨nsch et al. (1997) was also used to fit the LMCR to the same
variogram estimates, using the same initial values and weighting criteria. As set out
above, the algorithm combines standard non-linear regression, using a gradient-based
optimization, with Goulard and Voltz’s algorithm to fit the coregionalization matrices, S,
l
ARTICLE IN PRESS
R.M.Lark,A.Papritz/Geoderma1962(2003)1–16 9
Table1a
ParametersoftheLMCRusedtoinitializethefitting
Visualguess
Guess Goovaerts
Cd Ni Zn Cd Ni Zn
S Cd 0.2 0.25
0
Ni 0.5 10 0.47 7.9
Zn 10 10 100 3.5 15.4 105
a 0.1 0.2
1
S Cd 0.4 0.43
1
Ni 1.0 5.0 0.49 2.0
Zn 10 20.0 400.0 8.3 24.5 398.0
a 1.5 1.3
2
S Cd 0.2 0.18
2
Ni 2.5 65.0 3.3 73.3
Zn 5.0 130.0 300.0 6.2 142.5 398.0
WSS 2.24e2 2.18e2
Guess—bestguessbyeyeusingaspreadsheettocheckthematrices.Goovaerts—resultsobtainedbyGoovaerts
(1997) S, S and S are, respectively, the coregionalization matrices for the nugget variogram function, the
0 1 2
sphericalvariogramfunctionofrangea,andthesphericalvariogramfunctionofrangea.Thequotedvalueof
1 2
WSSareforw(h)allequal.
i
subject to the condition of positive definiteness. The procedure starts with an initial non-
linear regression step where both the range parameters (including any further non-linear
parameters of the LMCR) and the entries b l of the coregionalization matrices are fitted
u,v
by (weighted) non-linear least-squares. With the non-linear parameters kept fixed, the
bl are then fitted by Goulard and Voltz’s algorithm to ensure that the coregionalization
u,v
matrices are positive definite. Given the updated estimates of the S, only the non-linear
l
parameters are then updated in the following non-linear regression steps. The alternation
Table1b
ParametersoftheLMCRusedtoinitializethefitting
Visualguess
Poor(1) Poor(2)
Cd Ni Zn Cd Ni Zn
S Cd 10.0 10.0
0
Ni 1.0 10.0 5.0 10.0
Zn 1.0 1.0 10.0 5.0 5.0 300.0
a 1.0 1.0
1
S Cd 10.0 10.0
1
Ni 1.0 10.0 5.0 50.0
Zn 1.0 1.0 10.0 10.0 20.0 500.0
a 2.0 2.0
2
S Cd 10.0 10.0
2
Ni 1.0 10.0 5.0 50.0
Zn 1.0 1.0 10.0 10.0 100.0 500.0
WSS 2.51e4 2.22e4
Thesearepoor-fittingparameters,chosentotesttherobustnessofthefittingmethod.
ARTICLE IN PRESS
10 R.M.Lark,A.Papritz/Geoderma1962(2003)1–16
Table2
EstimatedparametersoftheLMCR
Initialvalues:visualguess
SA NLR
Cd Ni Zn Cd Ni Zn
w(h ):equal
i
S Cd 0.135 0.132
0
Ni 0.382 6.42 0.342 6.56
Zn 0.457 1.18 1.79 0.726 2.48 4.05
a 0.099 0.104
1
S Cd 0.543 0.548
1
Ni 0.861 9.54 0.914 9.41
Zn 11.69 51.73 511.30 11.42 50.32 508.56
a 1.517 1.517
2
S Cd 0.180 0.177
2
Ni 3.09 67.75 3.082 67.74
Zn 6.00 130.51 399.59 6.01 130.70 399.86
WSS 1.916e2 1.916e2
w(h):meannumberofpairs
i
S Cd 0.226 0.138
0
Ni 0.441 7.01 0.343 6.15
Zn 2.51 3.55 28.44 0.922 1.73 6.24
a 0.108 0.099
1
S Cd 0.419 0.516
1
Ni 0.790 10.39 0.916 11.65
Zn 8.54 51.59 494.10 10.25 54.06 516.60
a 1.480 1.490
2
S Cd 0.208 0.200
2
Ni 3.07 65.29 3.04 64.95
Zn 6.99 126.05 380.32 6.87 125.53 381.11
WSS 2.248e5 2.246e5
Initialvalues:Goovaerts(1997)estimates
SA NLR
Cd Ni Zn Cd Ni Zn
w(h):equal
i
S Cd 0.102 0.115
0
Ni 0.256 6.41 0.321 6.44
Zn 0.856 2.74 8.59 0.355 1.20 1.10
a 0.104 0.098
1
S Cd 0.584 0.565
1
Ni 1.01 9.47 0.934 9.52
Zn 11.27 50.07 520.83 11.79 51.58 510.55
a 1.515 1.517
2
S Cd 0.171 0.177
2
Ni 3.07 67.87 3.08 67.75
Zn 6.03 130.69 400.97 6.01 130.72 400.97
WSS 1.916e2 1.916e2