空间统计ch9variogram精品课件

1、Chapter 9 (Semi)Variogram modelsGiven a geostatistical model, Z(s), its variogram g(h) is formally defined aswhere f(s, u) is the joint probability density function of Z(s) and Z(u).For an intrinsic random field, the variogram can be estimated using the method of moments estimator, as follows:where

2、h is the distance separating sample locations si and si+h, N(h) is the number of distinct data pairs. In some circumstances, it may be desirable to consider direction in addition to distance. In isotropic case, h should be written as a scalar h, representing magnitude.Note: In literature the terms v

3、ariogram and semivariogram are often used interchangeably. By definition g(h) is semivariogram and the variogram is 2g(h).1Robust variogram estimatorVariogram provides an important tool for describing how the spatial data are related with distance. As we have seen it is defined in terms of dissimila

4、rity in data values between two locations separated by a distance h. It is noted that the moment estimator given in the previous page is sensitive to outliers in the data. Thus, sometimes robust estimators are used. The widely used robust estimator is given by Cressie and Hawkins (1980):The motivati

5、on behind this estimator is that for a Gaussian process, we haveBased on the Box-Cox transformation, it is found that the fourth-root of 12 is more normally distributed.* Cressie, N. and Hawkins, D. M. 1980. Robust estimation of the variogram, I. Journal of the International Association for Mathemat

6、ical Geology 12:115-125.2Variogram parametersThe main goal of a variogram analysis is to construct a variogram that best estimates the autocorrelation structure of the underlying stochastic process. A typical variogram can be described using three parameters:Nugget effect represents micro-scale vari

7、ation or measurement error. It is estimated from the empirical variogram at h = 0.Range is the distance at which the variogramreaches the plateau, i.e., the distance (if any)at which data are no longer correlated.Sill is the variance of the random field V(Z),disregarding the spatial structure. It is

8、 theplateau where the variogram reaches at therange, g(range).hg(h)02468100.0range h = 5nugget = 0.2sill = 1.035 m20200.5 m0.5 mSetting variogram parametersConstruction of a variogram requires consideration of a few things:An appropriate lag increment for h It defines the distance at which

9、the variogram is calculated.A tolerance for the lag increment It establishes distance bins for the lag increments to accommodate unevenly spaced observations.The number of lags over which the variogram will be calculated The number of lags in conjunction with the size of the lag increment will defin

10、e the total distance over which a variogram is calculated.A tolerance for angle It determines how wide the bins will span.Two practical rules:It is recommended that h is chosen as suchthat the number of pairs is greater than 30.2.The distance of reliability for anexperimental variogram is h D/2, whe

11、reD is the maximum distance over the field of data.4Computing variogramsAn experimental variogram is calculated using the R function (in package gstat): variogrm(pH1,locgx+gy, soil87.dat)# gx: list or vector of x-coordinates# gy: list or vector of y-coordinates# pH: list or vector of a response vari

12、able 5Covariogram and CorrelogramCovariogram (analogous to covariance) and correlogram (analogous to correlation coefficient) are another two useful methods for measuring spatial correlation. They describe similarity in values between two locations.Covariogram:Its estimator:where is the sample mean.

13、At h = 0, (0) is simply the finite variance of the random field. It is straightforward to establish the relationship: The correlogram is defined as6Properties of the moment estimator for variogramIt is unbiased:If Z(s) is ergodic, as n . This means that the moment estimator approaches the true value

14、 for the variogram as the size of the region increases. The estimator is consistent. The moment estimator converges in distribution to a normal distribution as n , i.e., it is approximately normally distributed for large samples.For Gaussian processes, the approximate variance-covariance matrix of i

15、s available (Cressie 1985).* Cressie, N. 1985. Fitting variogram models by weighted least squares. Mathematical Geology 17:563-586.7Properties of the moment estimator for covarianceThe covariance: C(h) = cov(Z(si), Z(si+h)The moment estimator:Properties:The moment estimator for the covariance is bia

16、sed. The bias arises because the covariance function for the residuals, is not the same as the covariance function for the errors,For a second-order stationary random field, the moment estimator for the covariance is consistent: (h)C(h) almost surely as n . However, the convergence is slower than th

17、e varigogram.For a second-order stationary random field, the moment estimator is approximately normally distributed.Properties 1 and 2 are the reasons why the variogram is preferred over the covariance function (and correlogram) in modeling geostatistical data. 8pHFsrfxpyp01002003004000.00.

18、4ypxp01002003004000.00.000.25Variogram before detrending Variogram after detrending01002003004005000200400600800010020030040050002004006008009Variogram modelsThere are two reasons we need to fit a model to the empirical variogram:Spatial prediction (kriging) requires estimates of the var

19、iogram g(h) for those hs which are not available in the data.The empirical variogram cannot guarantee the variance of predicted values to be positive. A variogram model can ensure a positive variance.Various parametric variogram models have been used in the literature. The follows are some of the mo

20、st popular ones.Linear model where c0 is the nugget effect. The linear variogramhas no sill, and so the variance of the process is infinite. The existence of a linear variogram suggests a trend inthe data, so you should consider fitting a trend to thedata, modeling the data as a function of the coor

21、dinates (trend surface analysis).hg(h)10Power model - where c0 is the nugget effect. The power variogram has no sill, so the variance of the process is infinite. The linear variogram is a special case of the power model. Similarly, the existence of a linear variogram suggests a trend in the data, so

22、 you should consider fitting a trend to the data, modeling the data as a function of the coordinates (trend surface analysis).hg(h)a 111Exponential model - where c0 is the nugget effect. The sill is c0+c1. The range for the exponential model is defined to be 3a at which the variogram is of 95% of th

23、e sill.Gaussian model - where c0 is the nugget effect. c0+c1 is the sill. The range is 3a. This model describes a random field that is considered to be too smooth and possesses the peculiar property that Z(s) can be predicted without error for any s on the plane.hg(h)hg(h)12Logistic model (rational

24、quadratic model) - where c0 is the nugget effect. The sill is c0+a/b. The range for the exponential model isSpherical model - where c0 is the nugget effect. The sill is c0+c1. The range for the spherical model can be computed by setting g(h) = 0.95(c0+c1).for 0 h afor h ahg(h)hg(h)13Parameter estima

25、tionThere are commonly two ways to fit a variogram model to an empirical variogram. Assume the variogram model g(h; q), where q is an unknown parameter vector. For example, for the exponential variogram model q = (c0, c1, a).Ordinary least squares method The OLS estimator for q is obtained by findin

26、g that minimizesThe OLS estimation can be easily implemented in R using function optim or nls. Initial values for q are required, these values can be obtained from the empirical variogram.Notes: OLS estimation assumes that- does not depend on the lag distance hi- for all pairs of lag distances hi hi

27、.Both assumptions are violated. The variance and the covariance depend on the number of pairs of sites used to compute the empirical variogram (see Cressie 1985).These violations do not contribute significantly to the bias of the parameter estimation.14Weighted least squares estimatorThe WLS estimat

28、or for q is obtained by finding that minimizeswhereSo that,To note that the WLS estimator is more precise (has a smaller variance) than the OLS estimator.Model selection criteria: Select a model with the smallest residual sum of squares or AIC or log-likelihood ratio, but pay a particularly attentio

29、n to the goodness-of-fit at short distance lags (important for efficient spatial prediction).15Splus implementationFitvar.s(dt, c, sill, range, model, wt=F)# dt: list of distances and sample variograms obtained from the function variogrm# c: initial estimate of the nugget effect# sill: initial estim

30、ate of the sill# range: initial estimate of the range# model:= exp, fits an exponential model to the sample variogram #= gau, fits a Gaussian model to the sample variogram#= sph, fits a spherical model to the sample variogram#= lin, fits a linear model to the sample variogrampH.variog_variogrm(soil8

31、7.dat$gx,soil87.dat$gy,soil87.dat,5,nint=30,dmax=400)pH.exp_fitvar.s(pH.variog,0.12,0.22,300,model=exp)x_seq(0,380,1)lines(expvar(x,pH.exp),col=5)Note: In the function fitvar, wt=F (i.e., OLS) each sample equally contributes to the objective function Q(q), while wt=T (i.e., WLS) Q(q) is weighted in

32、proportion to the number of obs. used in computing the sample variance. Thus, locations based on a few obs. will not carry as much weight compared to the one based on a large number of obs.16Fractals The concept of dimension Geometric objects are traditionally viewed and measured in the Euclidean sp

33、ace, e.g., line, rectangle and cube, with dimension D = 1, 2, and 3, respectively.However, many phenomena in nature (e.g., clouds, snow flakes, tree architecture) cannot be satisfactorily described using Euclidean dimensions. To describe the irregularity of such geometric objects (irregular geometri

34、c objectsare called fractals), we need to generalize theconcept of Euclidean dimension.The Hausdorff Dimension If we take an objectresiding in Euclidean dimension D and reduce itslinear size by 1/r in each spatial direction, thenumber of replicas of the original object wouldincrease to N = rD times.

35、 D = log(N)/log(r), is theHausdorff dimension, named after the Germanmathematician, Felix Hausdorff. The importantpoint is that in fractal dimension D need not be aninteger, it could be a fraction. It has proved usefulfor describing natural objects. D = 1D = 2D = 3r = 1r = 2r = 3N = 1N = 1N = 1N = 4

36、N = 8N = 2N = 3N = 9N = 2717Examples of geometric objects with non-integer dimensions1. Cantor set (dust) Begin with a line of length 1, called initiator. Then remove the middle third of the line, this step is called the generator, because it specifies a rule that is used to generate a new form. The

37、 generator could iterativelyinfinitely be applied to the remaining segmentsso that to generate a set of “dust”. The dustsare obviously neither points nor lines, but laysomewhere between them, thus has a dimensionbetween 0 and 1:D = log(N)/log(r) = log(2)/log(3) = 0.6309.2. Koch curve D = log(4)/log(

38、3) = 1.2618.InitiatorGenerator3. Sierpinski triangle D = log(3)/log(2) = 1.5850 18Self-similarity and smoothnessAn important property of a fractal is self-similarity, which refers to an infinite nesting of structure on all scales. It means that a substructure resembles the form of its superstructure

39、, e.g., leaf shape resembles branch shape, whereas branch resembles tree shape.Another important way to understand fractal dimension is that D is a smoothness measure of a spatial process/object (e.g., surface smoothness/roughness). When D = 1 (a line), or = 2 (a plane), the objects are smooth. For those objects whose Ds are between 1 or 2 (e.g., Koch curve or Sierpinski triangle), their smoothness varies between a line and a plane.Study surface growth and smoothness is increasingly becoming an important physic and biological su