Anisotropic Models

The KRIGE2D Procedure

Anisotropic Models

In all the theoretical models considered previously, the lag distance h entered as a scalar value. This implies that the correlation between the spatial process at two point pairs P₁,P₂ is dependent only on the separation distance h = | P₁P₂ |, not on the orientation of the two points. A spatial process { $Z(r), r \in D$ } with this property is called isotropic, as is the associated covariance or semivariogram.

However, real spatial phenomena often show directional effects. Particularly in geologic applications, measurements along a particular direction may be highly correlated, while the perpendicular direction shows little or no correlation. Such processes are called anisotropic. Refer to Journel and Huijbregts (1978, section III.B.4) for more details.

There are two types of anisotropy. The simplest type occurs when the same covariance form and scale parameter c₀ is present in all directions but the range a₀ changes with direction. In this case, there is a single sill, but the semivariogram reaches the sill in a shorter lag distance along a certain direction.

This type of anisotropy is called "geometric" and is discussed in the following section.

Geometric Anisotropy

Geometric anisotropy is illustrated in Figure 34.10, where an anisotropic Gaussian semivariogram is plotted. The two curves displayed in this figure are generated using a₀=1 in the NE -SW direction and a₀=3 in the SE -NW direction.

Figure 34.10: Geometric Anisotropy with Major Axis along E -W Direction

As you can see from the figure, the SE -NW curve gets "close" to the sill at approximately h=2, while the NE -SW curve does so at h=6. The ratio of the shorter to longer distance is (2/6) = (1/3). This is the value to use in the RATIO= parameter in the MODEL statement in PROC KRIGE2D. Since the longer, or major, axis is in the NE -SW direction, the ANGLE= parameter in the MODEL statement in PROC KRIGE2D should be 45^o (all angles are measured clockwise from north).

The terminology associated with geometric anisotropy is that of ellipses. To see how this comes about, consider the following hypothetical set of calculations.

Let { $Z(r), r \in D$ } be a geometrically anisotropic process, and assume that there are sufficient data points to calculate an experimental semivariogram at a large number of angle classes $\theta \in \{0,\delta \theta, 2\delta \theta, ... ,180^o$ }. At each of these angles $\theta$ , the experimental semivariogram is plotted and the effective range $r_{\epsilon}$ is recorded. A diagram, in polar coordinates, of $(r_{\epsilon},\theta)$ yields an ellipse, with the major axis in the direction of the largest $r_{\epsilon}$ and the minor axis perpendicular. Denote the largest $r_{\epsilon}$ by $r_{\epsilon}^{max}$ , the smallest by $r_{\epsilon}^{min}$ , and their ratio by

$R = \frac{r_{\epsilon}^{min}}{r_{\epsilon}^{max}}$

By a rotation, a new set of axes are aligned along the major and minor axis. Then, a rescaling elongates the minor axis so its length equals that of the major axis of the ellipse.

First, the angle $\theta$ of the major axis of the ellipse (measured clockwise from north) is transformed to standard Cartesian orientation or counter-clockwise from the x-axis (east). Let $\varphi=90^o-\theta$ denote the transformed angle. The matrix to transform the distance h is in terms of $\varphi$ and the ratio R and it is given by

$H= ( \cos(\varphi) & \sin(\varphi) \-\sin(\varphi)/R & \cos(\varphi)/R)$

For a given point pair P₁P₂, with coordinates (x₁,y₁), (x₂,y₂), the transformed interpair distance is computed by first transforming the components $\delta x=x_1-x_2$ and $\delta y=y_1-y_2$ by

$( \delta x' \\delta y') = H ( \delta x \\delta y)$

The transformed interpair distance is then

$h' = \sqrt{(\delta x')^2+(\delta y')^2}$

The original semivariogram, a function of both h and $\theta$ , is then transformed to a function only of h':

$\hat{\gamma}(h') = \gamma(h,\theta)$

This single semivariogram is then used for kriging purposes.

The interpretation of the major/minor axis in the case of geometric anisotropy is that the direction of the major axis is the direction in which the spatial process ${\{Z(r), r \in D\}}$ is most highly correlated; the process is least correlated in the perpendicular direction.

In some cases, these directions are known a priori. This can occur in mining applications where the geology of a region is known in advance. In most cases, however, nothing is known about possible anisotropy. Depending on the amount of data available, using four to six directions is usually sufficient to determine the presence of anisotropy and find the approximate major/minor axis directions.

The most convenient way of performing this is to use the NDIR= option in the COMPUTE statement in PROC VARIOGRAM to obtain a separate experimental semivariogram for each direction. After determining the direction of the major axis, use a DIRECTIONS statement on a subsequent run of PROC VARIOGRAM with this direction and its perpendicular direction. For example, if the initial run of PROC VARIOGRAM with NDIR=6 in the COMPUTE statement indicates that $\theta=45^o$ is the major axis (has the largest $r_{\epsilon}$ ), then rerun PROC VARIOGRAM with

   DIRECTIONS 45,135;

Then, determine the ratio of $r_{\epsilon}$ for the minor and major axis for the RATIO= parameter in the COMPUTE statement of PROC KRIGE2D. This ratio is $\le 1$ for modeling geometric anisotropy. In the other type of anisotropy, zonal anisotropy, the RATIO= parameter is set to a large number for reasons explained in the following section.

Zonal Anisotropy

In zonal anisotropy, either the form covariance structure or the parameter c₀ (or both) is different in different directions. In particular, the sill is different for different directions. In geologic applications, this is the more common type of anisotropy. It is not possible to transform such a structure into an isotropic semivariogram.

Instead, nesting and geometric anisotropy are used together to approximate zonal anisotropy. For example, suppose the spatial process has a correlation structure in the N -S direction described by $\gamma_{z,1}$ , a spherical model with sill at c₀=6 and range a₀=2, while in the E -W direction the correlation structure, described by $\gamma_{z,2}$ , is again a spherical model but with sill at c₀=3 and range a₀=1.

You can approximate this structure in PROC KRIGE2D by specifying two nested models with large RATIO= values. In particular, the appropriate MODEL statement is

   MODEL FORM=(S,S) ANGLE=(0,90) SCALE=(6,3) 
         RANGE=(2,1) RATIO=(1E8,1E8);

The large values of the RATIO= parameter for each nested structure have the effect of an "infinite" range parameter in the direction of the minor axis. Hence, there is no variation in $\gamma_{z,1}$ in the E -W direction and no variation in $\gamma_{z,2}$ in the N -S direction.

Anisotropic Nugget Effect

Note that an isotropic nugget effect can be approximated by using nested models, with one of the nested models having a small range. Applying a geometric anisotropy specification to this nested structure results in an anisotropic nugget effect.

Chapter Contents
Previous
Next
Top