Mathematical Geology, Vol. 32, No. 6, 2000
Understanding Anisotropy Computations1 2 and . Siska2
Most descriptions of anisotropy make reference to reduced distances and conversion of anisotropic models to isotropic counterparts and equations are presented for a certain class of range-anisotropic models. Many descriptions state that sill anisotropy is modelled using a range-anisotropic structure having a very elongated ellipse. The presentations typically have few or no intervening steps. Stu- dents and applied researchers rarely follow these presentations and subsequently regard the programs that compute anisotropic variograms as black-boxes, the contents of which are too complex to try to understand. We provide the geometry necessary to clarify those computations. In so doing, we provide a general way to model any type of anisotropy (range, sill, power, slope, nugget) on an ellipse. We note cases in the literature in which the printed descriptions of anisotropy on an ellipse do not match the stated or coded models. An example is provided in which both range- and sill-anisotropic structures are fitted to the experimental variogram values from an elevation data set using the provided equations and weighted nonlinear regression. The original variogram values are plotted with the fitted surfaces to view the fit and anisotropic structure in many directions at once.
KEY WORDS: ellipse, nested structures, range, sill, slope, variogram.
Copyright By PowCoder代写 加微信 powcoder
INTRODUCTION
A number of books (Isaaks and Srivastava, 1989; Pannatier, 1996; Goovaerts, 1997; Myers, 1997) present the concept of geometric anisotropy using the re- duced distance technique and “conversion to an isotropic model.” While not al- gebraically incorrect, this conversion process is confusing to most students and many practitioners, who rarely follow these presentations and subsequently regard the programs that compute anisotropic variograms as black-boxes, the contents of which are too complex to try to understand. The purpose of this paper is to help clarify the geometry behind and computations for anisotropy. For those needing a strong foundation in all aspects of variogram estimation, such an understand- ing is requisite. Those interested only in fitting variogram models, and who are neither interested in the interpretation of model parameters nor will have cause to code their own programs, can continue to rely on available programs. We will
1Received 9 December 1998; accepted 29 June 1999.
2Department of Forest Science, Texas A&M University, College Station, Texas, 77843-2135. e-mail:
0882-8121/00/0800-0683$18.00/1 ⃝C 2000 International Association for Mathematical Geology
684 Eriksson and 1. Examples of variograms with anisotropy; A, sill anisotropy in direction 30, B, range anisotropy in direction 30, C, sill in direction 30 and range anisotropy in direction 90, and D, nested range anisotropic structures in directions 90 and 150.
restrict ourselves to the two-dimensional case; extension to three dimensions is straightforward.
There are three concepts that are muddled together in the presentation by Isaaks and Srivastava (1989) of anisotropy—namely, geometric anisotropy, zonal anisotropy, and nested structures. Zimmerman (1993) further noted the incon- sistent use, in the literature, of the terms geometric and zonal anisotropy. He suggested, instead, adopting the terms nugget anisotropy, range anisotropy, and sill anisotropy. Examples of range and sill anisotropies are shown for a spherical model in Figures 1A and 1B. In each case, the direction of the maximum, sill or range, is 30◦ north of east. Figure 1C combines range and sill anisotropies for a spherical model. The direction of maximum sill is 30◦ and the direction of maxi- mum range is 90◦. Complex variogram models are often constructed as the sum of simpler models. The simpler models are called nested structures. Figure 1D shows a variogram constructed as the sum of a spherical model with maximum range 100 in direction 90 and another spherical with maximum range 50 in direction 150. The axes in Figures 1A–1D are lag-distance h, angle of separation φ, and variogram value 2γ (h, φ). We assume that angles are measured counterclockwise from east. Some variogram and Kriging implementations use angles measured in degrees azimuth (clockwise from north).
Understanding Anisotropy
RANGE ANISOTROPY
Consider first the case of range anisotropy (Fig. 1B). It is assumed that the sill c of all directional variograms is the same. The basic idea behind most anisotropic modeling algorithms is that there is some separation direction θ in which the directional variogram has maximum range amax. Perpendicular to that direction, θ ± π/2, the range is assumed to be minimum, amin. The maximum and minimum ranges define an ellipse (Fig. 2) that lies in the (x, y) coordinate system, where x and y denote the rectangular locational coordinates, usually eastings and northings. That is, x is the lag distance in the x direction and y is the lag distance in the y direction. The range aφ of the directional variogram for any separation angle φ is assumed to lie on the ellipse.
Let u and v represent coordinates in a system oriented along the major and minor axes of the ellipse. Then the equation of the ellipse is
(u/a )2+(v/a )2=1 (1) max min
We can convert the rectangular uv coordinates to and from polar coordinates via the relationships
aξ= u2+v2 u=aξcosξ ξ = tan−1(v/u) v = aξ sin ξ
where ξ is the angle of separation relative to the uv coordinate system and aξ is the range in that direction (Fig. 2).
Figure 2. Geometry of a range-anisotropic ellipse as it relates to points pi and p j separatedbydistancehij inthedirectionφij.Therangeindirectionφij isaφij.
686 Eriksson and than directly specifying amax and amin, most implementations of vari- ogram modeling require that the user specify the magnitude a of either maximum or minimum anisotropy, and a factor η that relates a to amax and to amin. If η > 1, then amin =a and amax =ηa. Similarly, if η < 1, then amax =a and amin =ηa. For concreteness assume, for now, that η < 1. Upon substituting into (1) we see that the ellipse can be represented as aξ2cos2ξ +aξ2sin2ξ =1 (2) a2 (aη)2 and that the range in direction ξ, relative to the uv coordinate system, is aξ =aη η2cos2ξ+sin2ξ (3) Now consider two points pi and p j having locational coordinates (xi , yi ), (x , y ), and translate the two points to the origin p′ of the ellipse (Fig. 2). The = y2 + x2 and the separation angle, relative to the separation distance is h xy coordinate system, is φi j = tan−1(yi j /xi j ). Relative to the uv coordi- nate system the separation angle is ξi j = φi j − θ . From (3), the range in direction φij is a =a a a2 cos2(φ −θ)+a2 sin2(φ −θ) φij max min or using the anisotropy ratio, min ij max ij a =aη η2cos2(φ −θ)+sin2(φ −θ) (4) φij ij ij The i j subscripting is used to emphasize the fact that these distances and angles are defined by, and calculated for, the pair of points pi and p j . The range parameter, usually symbolized as a in an isotropic variogram equa- tion, can be replaced by Equation (4) to obtain an anisotropic variogram model. Figure 1B, for example, was created simply by plotting √ √ 5 1.5 h 0.32 cos2(φ−30)+sin2(φ−30) −0.5 h 0.32 cos2(φ−30)+sin2(φ−30) 3 100·0.3 100·0.3 = forh≤100·0.3/ 0.32 cos2(φ−30)+ sin2(φ−30) 5 otherwise. The sill is c = 5, the range, in direction θ = 30, is amax = 100, the anisotropy ratio is η = 0.3, so the minimum range, in direction 120, is amin = 30. Notice that the anisotropic variogram model is formally a function of both lag distance h and separation angle φ. Understanding Anisotropy 687 Making Sense of the Reduced Distance Approach As in the previous example of creating Figure 1B, Equation (4) can be used directly in the calculation of variogram values. However, most authors and algo- rithms use the reduced distance and standardized model approach. One could, of course, begin with Equation (4) to obtain the computational algorithm as presented, for example, by Isaaks and Srivastava (1989) or Goovaerts (1997), but it is more instructive to take another look at the geometry. The uv coordinate system is related to the xy coordinate system via the rotation matrix Again translate the two points pi and pj to the center of the anisotropy ellipse such that pi′ coincides with the center, and apply the rotation matrix to p′j : R= cosθ sinθ −sinθ cosθ uij = cosθ sinθ xij (5) vij −sinθ cosθ yij The intersection of line v = (vi j /ui j )u, defined by points pi′ and p′j , with the ellipse defines the point pφi j having distance aφi j from the center of the ellipse. So substituting v = (vi j /ui j )u into the Equation (1) of an ellipse, we find that u2 v2 ij+ij =1 ij max min Solving for u and substituting into the equation of the line gives u2 v2 u=u ij+ij φij ij a2 a2 max min u2 v2 v=v ij+ij φij ij a2 a2 max min 688 Eriksson and have been added to u and v to emphasize the fact that these are the uv coordinates of the point pφi j on the ellipse at separation angle φi j . Noting that u2 + v2 = x2 + y2 = h and using the definitions of ij uφij and vφij , the range in direction φij is found to be a = u2 +v2 =h a a a2 u2 +a2 v2 (6) φij φij φij ij max min min ij max ij This is a reexpression of Equation (4) based on the rotated coordinates. Another useful expression for aφi j in terms of the original coordinates is a =h bx2−bxy+by2 (7) φij ij 1 ij 2 ij ij 3 ij b =(cosθ/a )2 +(sinθ/a )2, b =2sinθ cosθ(1/a2 −1/a2 1maxmin2 minmax Zimmerman (1993) in his Equation (8). = (sin θ/a )2 + (cos θ/a )2. This is, essentially, the form used by min A number of the common variogram models incorporate the range as the divisor in the ratio h′ = hi j /aφi j . This term is referred to as the reduced distance. From (6) it is clear then that h′= ij+ij (8) ij a2 a2 The reduced distances can be used directly in the Gaussian, spherical, and expo- nential variogram models. For example, Figure 1B can be generated by using the rotation matrix R to convert (x,y) coordinates to (u,v) coordinates and then plotting 5(1.5 2γ(h,φ) = (u/100)2 + (v/30)2 − 0.5( for h ≤ u2 +v2 (30u)2 + (100v)2 otherwise. Combining the Steps (u/100)2 + (v/30)2)3) The computations leading to h′ can be combined as follows: ij ′ uφij = 1/amax 0 cosθ sinθ xij v′ 0 1/a −sinθ cosθ y φij min ij Understanding Anisotropy 689 and h′ = (u′ )2 + (v′ )2. This is the computational form presented by Isaaks ij φij φij and Srivastava (1989) and repeated by Myers (1997). It differs from Goovaerts (1997) only in that he assumed that angles were measured clockwise from the north and used the relationships between a, η, amax, and amin, rather than amax and amin, themselves. Journel and Huijbregts (1978) actually rotate back to the xy orientation before computing the “reduced distances,” but that reorientation is unnecessary. Many authors (e.g., Isaaks and Srivastava 1989; Pannatier, 1996; Goovaerts, 1997; Myers, 1997; Deutsch and Journel, 1998) state, effectively, that the above computations are equivalent to a rotation of axes, followed by the reduction of theseparationdistancetoh′ ,whichallowstheuseofanisotropicmodel.Thisis ij misleading. It is true that, given φij, the use of h′ =h′ /1 in place of the hij/aφij ij ij term(s) in a variogram model is equivalent to using a transformed variogram model having range 1. However, h′ is, itself, a function of φi j . In no way are all direc- tional models combined into a single isotropic model that has a range of 1. Im- plicitly,ananisotropicmodelisafunctionofbothhij andφij andcannotbemade isotropic. In practice, the use of Equation (6) or (7) is to be preferred over that of Equation (4) because they are computationally more efficient. Sine and cosine evaluations are slow operations and the use of (4) requires at least one sine and one cosine evaluation (depending on how it is coded) for each pair of points requiring a variogram estimate. If the R matrix is permanently stored, then calculation of aφi j in Equation (6) requires only one sine and one cosine evaluation. This fact may be what has lead to the use of the reduced distance concept and the resulting confusion concerning “equivalent isotropic” models. If one is simply fitting a vari- ogram model and is not concerned with the interpretation of the model parameters, then Equation (7) is preferred, as it requires neither the translation nor the rotation of axes. SILL ANISOTROPY With sill anisotropy, the sill varies with direction. Clark (1980), McBratney and Webster (1986), Isaaks and Srivastava (1989), among others, indicate that sill anisotropy is not nearly as common as range anisotropy. Sill anisotropy has been dealt with in a rather ad hoc manner. While range anisotropy has been clearly dealt with by modeling the range as an ellipse, rather than modeling the sill directly, it is usually assumed that directionally varying sills can be accommodated by using two or more range-anisotropic structures, one with a very large anisotropy ratio. 690 Eriksson and a Range-Anisotropic Structure with a Large Anisotropy Ratio The approach to fitting a sill-anisotropic variogram surface using a range- anisotropic structure with a large anisotropy ratio η can be thought of as a sequential process. The motivation for the approach is explained with the help of Figure 3. Assume, for the moment, that the variogram model displays no range anisotropy. An isotropic model that well fits the experimental variogram in the direction of minimum sill is computed. In Figure 3, this is represented in each panel by the lower of the two surfaces. The spherical model was used for both structures. The direction of minimum sill is 120◦ and the common range is 40. The sill in this direction is five. Now a second structure is added to the first. The sill of the second structure is taken to be the difference between the maximum and minimum sills of the directional variograms. In the case of Figure 3, the maximum and minimum sills are eight and five, respectively. The minimum range in the second structure is 40 and the maximum is set to a very large value. In Figure 3 the anisotropy ratio was set to 1000 so the maximum range was 40,000. If we look at the form of the spherical model the reasoning for using the elongated range ellipse approach to sill modelling becomes clear. Dropping the i j Figure 3. Variogram surfaces displaying sill anisotropy modeled using the ellongated ellipse approace. Parameters for the first structure, all four panels, are c1 = 5, a1 = 40, and η1 = 1. Parameters for the second structure are c2 = 3, θ2 = 120◦, η2 = 1000, and A, a2 = 40, B, a2 = 10, C, a2 =80, D, a2 =80. Understanding Anisotropy 691 subscripting used in the previous section, the range-anisotropic spherical model takes the form 2γ(h,φ)=c(1.5h/aφ −0.5(h/aφ)3) where, as before, aφ is the range in direction aφ . There will be a maximum lag distance that any application will use. It will be determined by either the boundaries of the geographical region under consideration or by the size of the window used for kriging, or both. In panels A–C, it is assumed to be 100. As φ approaches the direction of minimum sill for the second structure, which is also the direction of maximum range, the denominator in the h/aφ terms becomes very large and h/aφ → 0. Thus the contribution of the second structure, in the direction of minimum sill, is negligible. It is largest in the direction of maximum sill. The overall model for Figure 3A is 2γ(h,φ)=5(1.5(h/40)−0.5(h/40)3)+3(1.5(h/aφ)−0.5(h/aφ)3) where aφ is computed as in the previous section, with amin = 40 and η = 1000. We began the description of this example with the assumption that the under- lying variogram model is not range-isotropic. From Figure 3A, it is apparent that for angles of about 0–80, and at 120, the range seems to be close to 40. From about 0 to 80 the sill is constant at 8. The range in directions 80–170 is significantly different from 40 and the sills appear to change with direction. In fact, by the very act of using a range-anisotropic second structure, the range changes continuously from 40 to 40,000 and the sill is constant at 8. This is not an argument against the use of range-anisotropic structures as models of sill-anisotropic processes provid- ing the fitted variogram fits the experimental variogram well in all directions over the range of experimental data. Figures 3B and C differ from Figure 3A only in the value used for the minimum range of the second structure. For Figure 3B it was set to 10 and in Figure 3C it was 80. The variogram surface in Figure 3C is generally smooth throughout except for the bend at the range of the first structure and for the crease in the direction of minimum sill. Indeed, models of this form always produce a crease in the direction of minimum sill. Figure 3D is identical to Figure 3C except that the length of the lag distance axis was increased to 400 to better indicate the constant sill and the crease in the direction of the minimum sill. Modeling the Sill as an Ellipse The elongated ellipse procedure, just described, has become the standard for fitting variogram models to experimental variograms exhibiting directionally varying sills. Sill-anisotropic variograms, alternatively, can be fitted using the same procedure used to model range anisotropy. Let cmax be the maximum sill, cmin be the minimum sill, ν = cmin/cmax be the sill anisotropy ratio, and θc be the direction of maximum sill. Then the sill in 692 Eriksson and Siska direction φ can be modeled as an ellipse using c =c ν/ ν2cos2(φ−θ)+sin2(φ−θ) φmax c c b1x2 − b2xy + b3y2 u2φ + vφ2 = h · cmaxcmin/ (cminu)2 + (cmaxv)2 u = cosθc sinθc xij v −sinθc cosθc yij b =(cosθ/c )2+(sinθ/c )2 1 c max c min b =2sinθ cosθ 1 c2 −1 c2 2 c c min max b =(sinθ/c )2+(cosθ/c )2 3 c max c min As an example, Figure 1A was generated by plotting 2γ (h, φ) = (3/ (3/8)2 cos2(φ − 30) + sin2(φ − 30))(1.5(h/100) − 0.5(h/100)3) for different values of h and φ. Indeed, any directionally varying parameter, whether it be the range, sill, slope, nugget, or whatever, can be modeled using this approach provided that the assumption that the directionally changing param- eters lie approximately on an ellipse is reasonable. Zimmerman (1993) expressed concern regarding the use of sill-anisotropic models because they are not second-order stationary. The question of sill-anisotropy and stationarity will be addressed in another paper. POWER MODELS As usually written, the power model has the form 2γ (h) = cha , where 0 < a < 2. As just noted, either or both of the parameters in this model can be modeled as an ellipse yielding, for example, 2γ(h,φ) = cφhaφ (9) Pannatier (1996) and Deutsch and Journel (1998) both state that their software handles anisotropy in the correct manner by calculating an anisotropic distance Understanding Anisotropy 693 and leaving the parameter a unchanged. Examination of the GSLIB code, however, and empirical experience with VARIOWIN-fitted models reveals that, in the case of η < 1, both programs use the equation 2γ(h,φ) = c′ (h′)a = c′ (c′ /c′ )aha max max max φ c′ =c′ c′ c′2 cos2(φ−θ′)+c′2 sin2(φ−θ′) (10) φ max min min max is an elliptically varying parameter. This is because these programs compute h′ as u′ 2 + (v′/η′)2. The result follows from Equation (6). The meaning of the primes will be made clear in a moment. No other programs were checked, but we suspect that this usage is common. There is nothing inherently wrong with this model since c =c′ (c′ /c′ )a can be considered 程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com