Inference of locally varying anisotropy fields from diverse data sources

Computers & Geosciences 82 (2015) 170–182

Contents lists available at ScienceDirect

Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo

Inference of locally varying anisotropy fields from diverse data sources Maksuda Lillah a,b, Jeff B. Boisvert a,b,n a b

University of Alberta, Civil and Environmental Engineering, Centre for Computational Geostatistics (CCG), Canada 3-133 Markin/CNRL Natural Resources Engineering Facility (NREF), University of Alberta, Edmonton, Alberta, Canada T6G 2G7

art ic l e i nf o

a b s t r a c t

Article history: Received 23 August 2014 Received in revised form 16 April 2015 Accepted 27 May 2015 Available online 10 June 2015

A locally varying anisotropy (LVA) field characterizes the magnitude and direction of anisotropy in a modelling domain and can provide valuable information to geostatistical modeling methodologies that account for non-linear spatial features, in particular, second order non-stationarities that are observed in complex geological formations. In this framework, the primary variable to model (i.e. grade, porosity, concentration, etc) is measured at discrete spatial locations and there may or may not be exhaustively sampled gridded secondary information. Existing geostatistical techniques require knowledge of the magnitude of anisotropy (correlation lengths) and local orientations (strike, dip, and plunge) exhaustively in the domain. Inference of these parameters from samples of the primary variable is difficult when anisotropy varies locally; the correlation lengths and orientations of the discrete primary variable must be estimated at all locations in a domain and are not globally constant. Estimation of the required LVA field parameters depends on the available data and is unique to the type of information available. This work proposes three different methodologies for LVA field inference, including (1) use of a correlated exhaustively sampled secondary variable (2) use of direct point measurements of the orientation of the primary variable, and (3) use of a geological model that displays the desired LVA characteristics. Often a gridded exhaustively sampled secondary variable is available from geophysical surveys and displays similar spatial features as the primary variable. An adaptive moving window technique is proposed to infer the required LVA field from correlated exhaustive secondary information. The local anisotropy orientation is obtained from the gradient of the secondary information and the appropriate window size is calculated automatically, which is important when the scale of features vary in the domain. Secondly, direct angle measurements of orientation may be available from dip meter or borehole images, however, these data are axial in nature (when plunge is assumed to be 0°) or they represent a 3D coordinate system (when considering strike, dip and plunge); preprocessing is required before the LVA field can be estimated from point measurements of orientation because of the wrapping effect at 0°. The proposed technique treats point measurements of orientation as quaternions and assigns weights to nearby samples to estimate the exhaustive LVA field. In addition to the proposed methodologies for exhaustive and point data, a skeletonization technique is implemented to determine orientations from geological interpretations. The geological object is reduced to a spline curve and the direction of the curve is used as the local LVA field orientation. The suite of proposed LVA inference techniques are implemented in both 2D and 3D and illustrated on several data sets. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Locally varying anisotropy Stationarity Geostatistics Modeling Orientation Resource estimation

1. Introduction Complex nonlinear features such as fluvial channel systems, vein type deposits, uranium role fronts, etc, are abundant in nature and are characteristic of many earth science modeling problems; often a second order stationarity assumption (the correlation function is constant everywhere in the modeling domain) is

n Corresponding author at: University of Alberta, Civil and Environmental Engineering, Centre for Computational Geostatistics (CCG), Canada. E-mail address: [email protected] (J.B. Boisvert).

http://dx.doi.org/10.1016/j.cageo.2015.05.015 0098-3004/& 2015 Elsevier Ltd. All rights reserved.

inappropriate because of complex local features. Moreover, the complexity of natural deposits may preclude the use of redefining the domains (subdomaining) such that the resulting domains can be considered stationary; the domains would be so small there would be inadequate data for statistical inference. Correctly accounting for stationarity in these situations is critical for accurate geostatistical modeling. Anisotropy refers to properties in geological formations that show preferential directions of continuity, i.e. properties are more continuous in one orientation than in another. Consider an undulating structure where the deposit continuity spreads erratically in various directions. In such formations no single direction of

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182

anisotropy can be imposed, instead the dispersion is complex and changes locally. If information on this non-stationary anisotropy is available, it can be used in estimation or simulation to generate more geologically realistic predictions. Information on anisotropy can be incorporated by generating a locally varying anisotropy (LVA) field that exhaustively defines the magnitude and direction of anisotropy. In this context, it is assumed that there are some point measurements of the primary variable of interest, lacking this information, it is inappropriate to implement the majority of available geostatistical techniques. Moreover, it is assumed that the primary variable of interest is sparsely sampled; if the variable is very densely sampled nearly any estimation or simulation technique could be applied and the local features would be carried through to the models; in this case, the models are entirely controlled by the existing data rather than the modeling technique. There is no standard density to define ‘densely sampled’ and this depends heavily on the spatial structure of the variable. Typically this is assessed visually by applying kriging or inverse distance with a modeled covariance function and observing the spatial features of the models. If the desired LVA features are not present, advanced techniques that explicitly consider LVA should be considered. It is also assumed that there is some justification for believing that second order stationarity is not appropriate for the variable of interest. In geostatistical modeling, second order stationarity is more of a modeling decision than an assumption, McLennan (2007) provides clear guidance on assessing the decision of stationarity. Often, this decision is made based on geological knowledge of the deposit as well as assessing all available data in the domain and looking for subsets that violate the assumption of stationarity. Interested readers are referred to McLennan (2007) for more specific tests. There is significant interest in incorporating LVA into geological

171

modeling. The characterization of anisotropy is important in assessing variable continuity, establishing spatial interdependences, making decisions of stationarity, boundary modeling and identifying geological features, all of which lead to better strategic decisions such as well/drill hole placements and ultimately building better geostatistical models for resource estimation. Previous studies have demonstrated an increase in model accuracy when LVA has been correctly incorporated into modeling (Boisvert and Deutsch, 2011). There are numerous geostatistical techniques that can be implemented to incorporate LVA; however, inference of the input LVA field is still difficult and is the main focus of this work. Several variants of kriging and simulation have been proposed that utilize anisotropic modeling. Stroet and Snepvangers (2005) proposed local anisotropy kriging (LAK) which iteratively combines a gradient algorithm with kriging. The information of locally varying or curvilinear features is determined automatically from the point samples of the primary variable of interest using a local gradient function, but is only implemented in 2D and cannot use an external LVA field. The benefit if using an external LVA field is that additional information from secondary data, geological interpretations, expert knowledge and point measurements of orientation can be integrated into the modeling of the primary variable of interest, this is not possible with LAK. Pardo-Iquiza and Chica-Olmo (2007) also propose a universal cokriging methodology that can incorporate local orientation measurements and primary data in estimation mode. The standard/traditional approach for incorporating an LVA field into kriging or sequential Gaussian simulation, is to consider a local search when estimating or simulating (Deutsch and Lewis, 1992; Xu, 1996 to name a few). Practical implementations of kriging and SGS use local search neighbourhoods to restrict the number of data used. To incorporate LVA the gridded correlation function parameters (from the exhaustively modeled LVA field) are assumed to be stationary within this window when estimating or simulating at a given

Fig. 1. When estimating the gray locations with the surrounding data (red points) the anisotropy at the estimation location is applied everywhere and kriging is performed (Boisvert and Deutsch, 2011). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

172

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182

location (Fig. 1). This works well for incorporating LVA fields that are smoothly changing across large distances in the domain; however, if there are local features of interest that causes the LVA field to vary at the scale of the sampling, the local search window cannot account for these features. Boisvert and Deutsch (2011) use the shortest non-linear path to incorporate LVA prior to kriging to better account for the local features inherent in many complex deposits. Many of the techniques for incorporating LVA into geostatistical modeling require an exhaustive model of the local anisotropy and there is a need to develop a suit of techniques that can be used to infer this LVA field if these non-stationary geostatistical techniques are to be implemented effectively. In this work the LVA field is parameterized by five scalar values (strike, dip, plunge, major anisotropy ratio and minor anisotropy ratio). Anisotropy directions are completely defined in 3D by the three angle parameters: strike (α), dip (β) and plunge (φ) (Deutsch and Journel, 1998). The magnitude of anisotropy is parameterized with respect to the principal or major direction of continuity, by two ratios of the range: r1 ¼minor correlation length/major correlation length. r2 ¼vertical correlation length/major correlation length. These five scalar parameters must be known at every location in the modeling domain and constitute the LVA field described in this work. For simplification, these parameters are estimated on a Cartesian grid as the majority of geostatistical techniques use a Cartesian grid and within each cell the parameters are assumed constant; however the methodologies presented in this work could be easily adapted to estimate the LVA field at any arbitrary point in space if a grid free modeling methodology was desired. The constraints on this field differ slightly depending on the LVA methodology implemented. If the standard local search reorientation technique is used (Deutsch and Lewis, 1992) the 5 parameters must be used in conjunction with a covariance function (often a spherical variogram model) that is positive definite. If the technique of Boisvert and Deutsch (2011) is used, there are no theoretical constraints on these parameters, they need not be differentiable, nor correlated, nor positive definite; however, experience has shown that poor models are generated if the parameters are not smoothly varying, at least at the scale of the primary data sampling density. Consequently, in this work, the LVA field is parameterized by the exhaustive estimation of these five parameters on a Cartesian grid. This does not constitute a local covariance function; rather the field is simply the five parameters described above. In 2D, only the strike and a single ratio are required to delineate the LVA field. This work focuses on three different sources of information that can be used to infer the LVA field (1) exhaustive data such as outcrop images and geophysical surveys, (2) point data such as dip meter measurements, and (3) geological models that have been created prior to building geostatistical property models. First, a structure tensor is used to characterize the LVA for exhaustive data. Next, an adaptive window scheme is presented to consider features at different scales as this is a significant limitation when using a fixed moving window. Following this, 3D point data are estimated exhaustively with consideration given to the axial nature of the data (ignoring plunge) and the coordinate nature of the data (strike, dip and plunge). Finally, techniques to extract orientations from geological models are explored.

required variable. A correlated secondary variable from the exhaustive survey data can be used when the spatial features of the secondary data are deemed representative of the spatial features of the primary data of interest. Analyzing the intensity changes from exhaustive data (or calculating gradients of gray values in outcrop images) is an important step for LVA field generation. 2.1. Estimation of local orientation using gradients The first step in detecting orientations from exhaustive data is to search for local symmetries in the intensity distribution, which is the basis for finding a local orientation (HauBecker and Jahne, 1996). The term local orientation means that the change in intensity values show a distinct pattern or orientation and this orientation can be shown by a vector in 2D or orientation angle in higher dimensions. Before the formulation of tensors in 2D and 3D, there is a need to discuss how clear dominant directions can be identified and how a single direction is obtained. First, the gradient field, g(u), is calculated from the exhaustive data. The gradient transformation of the image in 2D has two vector components which contain the information of change in signal intensities. The gradient direction indicates the maximal change in intensity and is normal to the oriented structure. If gradient vectors cannot identify a dominant direction, the feature is assumed to be isotropic. Consider a local neighbourhood A around an arbitrary point u0, and a gradient field ∇g(u) for every point u in the domain. The objective is to find a mean vector a(u0) that unifies all the directional vectors of g(u) to give a single local orientation and is invariant to a rotation of π. It is important that the mean a(u0) maximize the average angles with respect to all direction vectors, g(u), in A (Feng and Milanfar, 2002). A schematic representation of the problem is shown in Fig. 2 A projection of g(u) onto a(u0) or the inner product of the vectors gives a measure of the proximity of the two. Moreover, since antipodal directions are the same, the absolute value or square of the inner products can be taken. In order to find a(u0) that lies along the feature and in the direction of the least change in intensity field, the problem is formulated as the following minimization problem (Eq. 2):

p(u) = g(u)⋅a(u0) = gT(u)⋅a(u0) Q(u) =

∫A(u ) p(u)2du

(1)

o

⎛ min(Q ) = min⎜ ⎝

⎞

∫A(u ) p(u)2du⎟⎠ o

Q can be expressed of in terms of a symmetric tensor T and the vector a(u0).

2. Data type 1: LVA field estimation from exhaustive secondary data Geophysical remote sensing data from seismic and other surveys provide low-resolution, bulk property measurements of the domain which are used in determining regional continuities for a

(2)

Fig. 2. Determining the orientation of local structure within a window.

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182

173

Q(u0) = aT(u0)T(u0)a(u0) where, T(u0) =

Fig. 3. Greyscale image of a folded outcrop (source: http://gsfc.nasa.gov/Sect2/ Sect2_1a.html) with the gradient based implementation (arrows) overlain obtained from local windows (bottom left).

(3) T

∫A(u ) g(u)⋅g (u)du o

Tensors have a wide range of application in many fields as they are generalizations of vectors in higher dimensions. They are an important tool in image processing in perceiving textures, edge detection from local changes in contours, intensity signals and directions of adjacent patterns. The motivation for such analysis comes from the fact that physical quantities of complex natural phenomenon – stress at a point on rock mass due to internal forces, electrical or thermal conductivity in anisotropic medium, diffusion, etc – cannot be understood by the use of scalars or vectors, they must be formulated as tensors (Cammoun et al., 2009). In a geostatistical setting, a tensor can be used to characterize spatial anisotropy and can be extracted from remote sensing data. This data can come from seismic, magnetic or gravitational surveys, outcrop images, geological interpretations, etc. It is assumed that the measured bulk variable has the correct spatial characteristics and can be used as a proxy for the underlying LVA field of the variable of interest. T (Eq. 3) expresses the outer product of the gradient vectors with itself, which is the result of rotation invariance (Cyganek and Siebert, 2009). From Eq. 3 it can be deduced that the minimum is found when a(u0) is the eigenvector corresponding to the smallest eigenvalue of T. The task of determining a(u0) is equivalent to eigen-decomposition of the symmetric tensor T, and the dimensionality of T follows that of the vector g(u) (Cyganek and Siebert, 2009) allowing for similar analysis in higher dimensions. First, intensity gradients ∇g(u) are computed from the exhaustive data. Then for a local neighbourhood A, about a node u0 the tensor components are:

Tij(u0) =

∫A(u −u) ∇g(u)∇gT(u)du

(4)

0

Eq. (4) is used when the intensity changes are continuous. An equivalent expression for discrete data takes the form:

^ Tij(u0) = F Di , Dj

(

)

(5) th

Dk is a discrete differentiating operator in the k axis, and F() is usually a Gaussian filter used for smoothing. Once the tensors are constructed in two and three dimensions, the next step is to extract the dominant direction through eigenvalue decomposition. The eigenvectors obtained preserve the invariance to a rotation of π. So, for a given symmetric tensor T the desired orientation of a can be obtained numerically as:

Fig. 4. Above: a large window size is selected for the domain. While orientations for larger scale features are detected correctly, continuities at a finer scale are not determined correctly. Below: continuity at a larger scale cannot be determined and is considered isotropic (R¼ 0).

⎡T11 T12 ⎤ ⎡T − T11⎤ ⎥ a= ⎢ 22 ⎥ T=⎢ T T ⎣ 21 22 ⎦ ⎣ 2T12 ⎦

(6)

There are several ways to obtain a, for instance by rotating T until it lies at the coordinates of its principal axis or by singular

Fig. 5. Impact of noise on orientation estimation and reliability. Plot dimensions are 500  300 units.

174

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182

Fig. 6. LVA field generated using large (above) and small (below) windows.

Fig. 7. LVA field generated using adaptive windows. The shape and dimensions of the search window at three locations is indicated.

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182

175

Fig. 8. Left: two data points are considered such that they point in opposite direction in the same axis i.e. antipodal pairs. Without preprocessing the input the axial nature of the point data is not considered in kriging. Right: estimation is performed by considering the nearby combination of point data that has the smallest spherical angle variance.

Fig. 9. Estimating at an unknown location (?) with three conditioning data located at u1, u2 and u3.

Fig. 10. The shaded area represents the shape of the underlying feature. Vectors represent point measurements of the features orientation at discrete locations. In this case there are 27 measurements of orientation.

Table 1 Angle variance for all data (Fig. 9) combinations. 1 Indicates the orientation has been flipped. u1

u2

u3

Variance

0 1 0 0

0 0 1 0

0 0 0 1

0.60 0.60 0.58 0.01

value decomposition (SVD) (Cyganek and Siebert, 2009). Feng and Milanfar (2002) implemented a principal component analysis (PCA) based feature orientation technique where local orientation is established by performing SVD onto a collection of gradient vectors. A central feature of the technique is to subdivide the gradient field into N blocks (Eq. 7) and perform SVD decomposition on the gradient vectors of each block. In 2D the following matrix is constructed for the individual blocks: Fig. 11. LVA field generated by recombining the kriged X, Y and Z components.

⎡ ∇f (1)T ⎤ ⎥ ⎢ ⎢ ∇f (2)T ⎥ G=⎢ ⎥ ⎢ ⋮ ⎥ ⎢ ∇f (N )T ⎥ ⎦ ⎣ where, ∇f (k ) = ⎡⎣∂f xk , yk /∂x T

(

)

(7) ⎤ ∂f xk , yk /∂y⎦ are the gradients at

(

)

each point k, and G is an N  2 matrix. The SVD decomposition of the form G ¼USVT, where U is an N  2 matrix with orthogonal columns; S is a diagonal matrix containing the singular values in sorted order; and V is 2  2 orthogonal matrix in which the columns contain information about

176

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182

Fig. 12. Orientations found using quaternions and inverse distance weighting of two samples: sample A at x ¼1, y¼ 0 and z ¼0 has a strike ¼45, dip ¼ 20, plunge¼ 10 and sample B at x ¼6, y¼0 and z ¼0 has a strike dip and plunge of 0.

Fig. 13. A synthetic 3D model used to demonstrate the quaternion averaging result, point measurements of orientation are shown as red vectors. (Left) Input samples. (Right) Model populated with estimated LVA orientations, estimated orientations are show as red vectors. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

dominant image orientations. The following steps are used to determine the LVA field orientation in the spatial domain: 1. Calculate the gradients at each point in the image of the model; Dk is the gradient in the kth coordinate direction. 2. Find necessary products of Di and Dj for every location, and if necessary apply a low pass filter to obtain the components of the tensor matrix T. 3. Similarly, a PCA based analysis can be performed, and the gradient vectors are arranged in an N  2 matrix G (Eq. 7) In 2D the tensor matrix for an image location (i,j) is: ⎡F D, D F D, D ⎤ i j ⎥ ⎢ ( i i) l Tij(u0) = ⎢ ⎥ ⎢⎣ F Dj , Di F Dj , Dj ⎥⎦ 4. Local orientations are obtained from the tensor matrix through eigenvalue decomposition. In the PCA implementation, similar calculations are performed by considering SVD of G.

(

( ) (

) )

The gradient method is applied to an outcrop image (Fig. 3) and the associated LVA field using the gradient based method shows that local orientations are largely consistent with the geologic structure. The gradient (change in the variable in a particular direction) is directly related to concept of the orientation of the

major and minor directions of continuity. The direction of maximum change is captured by the direction with the largest change in values (highest gradient) and thus represents the minimum direction of continuity. Conversely, the direction of minimum change is captured by the direction with the smallest change in values (lowest gradient) and thus represents the minimum direction of continuity. Note that this tensor formulation is only presented here as one methodology for calculating a gradient of a local window. There are other techniques that could be used to obtain this orientation, such as: the moment of inertia (Hassanpour, 2007), but we have found the gradient to be more stable, or; Hristopulos (2002) but this is only implemented in 2D and requires a simplex search optimization at each window that would be CPU intensive for large 3D grids, or; Spiliopoulos et al. (2011) however they require differentiability of a local covariance function, no covariance function is assumed in the above workflow. 2.2. Adaptive window for exhaustive data With the gradient feature extraction methodology discussed previously, the nature of the LVA field is largely dependent on the size of the local neighbourhood. The selection of window size becomes even more relevant when geological deposits display

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182

177

Fig. 14. Cross sections taken from the 3D model (Fig. 13). Slices taken at y¼ 21, 26 and 31.

different continuity at different scales. To highlight the problem of choosing an appropriate window size, consider a synthetic example where the feature scales vary (Fig. 4). When a large window is selected, features at that scale can be effectively captured (Fig. 4 above); however, the appropriateness of this search area is questionable for smaller features. When scales of geological continuities differ, orientations at every scale cannot be determined correctly by using a single preset window dimension (Fig. 4 below). Instead, if the search areas are allowed to vary such that they adapt to the underlying feature scale, a continuous LVA field can be generated that accounts for all apparent features, regardless of scale. A similar idea is implemented in Katkovnik et al. (2002) where an intersection of confidence intervals (ICI) rule is used to select an optimal window size or signal bandwidth for filtering image signals. Here, an ellipse is used to define the window or local neighbourhood (an ellipsoid in 3D). By locally redefining search windows, local neighbourhoods can be scaled and oriented in a way such that the scale dependent features are accounted for. The parameters of the adaptive window are the rotation angle (θ) the semi major and minor axes lengths. The user selects a range of parameters and the local reliability (R) for every parameter combination is calculated; the window with the largest reliability is selected. Reliability is a measure of confidence in the estimated orientation (Eq. 8). Reliability could also be considered a relativistic measure of the magnitude of the local anisotropy.

R=

λ1 − λ2 , forλ1 > λ2 λ1 + λ2

(8)

where λ1 and λ2 are the two non-zero eigenvalues of the symmetric

tensor Tij. R is limited between 0 and 1, and is a measure of how dominant one direction of intensity change is compared to another. The range of correlation in a particular direction relates to the magnitude of the eigenvalues; thus this ratio of eigenvalues relates to the strength of anisotropy in the local window. When λ1 is similar to λ2 R is approximately 0 and indicates an isotropic area; when λ1 ≫ λ2, R approaches 1 indicating a strong degree of anisotropy. Fig. 5 shows estimation of the dominant feature direction is robust to low levels of noise while the reliability drops; increasing noise results in a window that is approaching an isotropic state. The optimal window selection process searches over all user defined window parameter combinations to find the ellipsoid with highest reliability measure. The user selects a number of window sizes, defined by window length in x, y and z. This process is done for every local window defined by the modeling grid used. Calculation of all parameter combinations at each location may be computationally excessive for large grids or many window parameters. Instead a partial optimization is implemented where a smaller range of parameters are searched (Fig. 6). The steps involved in the adaptive window selection with partial optimization are: 1. Consider the exhaustive set parameters in choosing the local window for the first grid node (u0) and store the local orientation that corresponds to the highest reliability. 2. When considering an adjacent node (ui þ 1), the optimized values of the previous node (ui) are used as the initial window

178

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182

Fig. 15. Iterative thinning by K3M; different stages of skeletonization are shown.

parameters. Consider window size parameters smaller and larger than this initial size until the optimum reliability is found. This assume convexity of R. When considering the small window size only (Fig. 6 above) the orientations inside the larger features are modeled incorrectly. When considering the large window size only (Fig. 6 below) the small scale features are not reproduced well, Fig. 7 shows how the search windows vary throughout the field. Reliability is used as an index of confidence or accuracy of the orientation derived from a structure tensor. In this application, a reliability of 1 indicates the presence of very strong local anisotropy. As the gradient vectors are aligned along a clear dominant direction, the eigenvalue associated to the smallest gradient change predictably, λ2 approaches zero and R-1. When λ1 E λ2 there is no preferred direction of local intensity change and R-0, the region is assumed isotropic.

3. Data type 2: LVA field estimation from discrete point measurements of orientation Direct angle measurements at discrete locations may be available from outcrops, dipmeter measurements, borehole cameras and full-bore formation micro imagers (FMI). Inference of LVA from such point data is not straightforward and requires preprocessing since orientations can be axial in nature. In the first proposed methodology the common case of plunge¼ 0° is

considered, simplifying the inference problem as only strike and dip are required to fully define local orientation. Standard geostatistical algorithms consider axial (or undirected) orientations; strike α is identical to α þ π. In dealing with directional data in two dimensions, a standard procedure is to decompose the angles into circular coordinate vectors (C ¼cosα and S ¼sinα) pertaining to angles that wrap on the unit circle. This convention cannot be applied directly as axial data are 0r α r180°. Mardia and Jupp (2000) suggest doubling the data such that observations are 0 r α r360° and through standardization, data wraps at 2π. Prior to doubling, data are restricted to [0, π]. Standard circular transformations can then be applied to the doubled data. Once the observations are in circular components, estimation (or simulation) is done using kriging (or simulation). It is necessary to reproduce the correct correlation between vector components C and S, and a strong relationship may warrant cokriging or cosimulation (Boisvert, 2010). Three dimensional data are more complex as angle ‘doubling’ is not sufficient to consider the axial nature of data defined by a strike and dip. In this case, both the strike and dip angles wrap at 0°, there is no preprocessing of the angles that would allow for an appropriate decomposition into coordinates that can be estimated with standard techniques. Each input is defined by strike (α) and dip (β) and is decomposed into spherical polar coordinates. The direction X is decomposed as:

X(u) = (sinαcosβ, cosαcosβ, sinβ ) Using kriging to interpolate an exhaustive LVA field with this

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182

179

Fig. 16. Thinning in 3D using the 8-subiteration technique: Upper left: the model after smoothing by MAPS. Upper right: associated centerline. Below: resulting LVA field generated from the centerlines.

decomposition would give erroneous results as points with a strike 45° and 225° are recognized opposite whereas they lie in the same axis and are considered the same orientation in the case of geostatistical modeling, assuming symmetrical covariance functions (Fig. 8 left). One solution is to orient all of the directional vectors in the correct half-sphere so that the above decomposition into components would result in appropriate estimates of orientation. Thus, there is a need to identify which orientation (flipped or not flipped) results in the correct estimation of orientation. The variance of a set of angles (Mardia and Jupp, 2000), provides a good measure in choosing axial orientations. For any number of sampled points the aim is to align data on the same axis or along axes with the smallest variance. It is important to distinguish this angle variance from the typical kriging variance used in estimation. This angle variance is only used to determine which orientation results in the most consistent estimate for a given location. The angle variance is

computed for all combinations of orientations of the surrounding conditioning data found within the estimation search window during kriging. Consider three sample points being used to estimate at an unknown location (Fig. 9). The proposed technique is to determine which axial orientation would be consistent for all three locations. There are a total of 3 different combinations and the angle variance can be used to determine the appropriate orientation. Consider flipping one location at a time and recording the associated spherical angle variance. The combination where location u3 is ‘flipped’ (consider the opposite axis) has the minimum angle variance (Table 1). More combinations are required for considering a larger number of conditioning data, but flipping one orientation at a time and passing through the conditioning data once is sufficient to find the optimal number of flips; this process is not computationally demanding. The combination with the minimum variance is selected and traditional estimation of the vector components is implemented. The methodology is

180

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182

Fig. 17. Parameters for program gen_lva.

Fig. 18. Parameters for program quat_avg.

Fig. 19. Parameters for program k3m.

implemented on a synthetic 3D example with 27 discrete measurements of orientation (Fig. 10) resulting in the LVA field in Fig. 11. The benefit of considering flipping angles is that any geostatistical estimation or simulation methodology could be applied to generate the exhaustive LVA field. The drawback is that a plunge angle cannot be considered. In the cases where a plunge angle is required, the LVA field orientation can be defined using quaternions. Quaternions are 4-tuple entities and are routinely used in determining attitudes of airplanes and spacecrafts by combining data from several attitude sensors. Quaternion representation is based on Euler’s rotational theorem whereby the relative position of any two arbitrary, rigid coordinate systems can be represented

by a single rotation about a fixed axis (Grobekatthofer and Yoon, 2012). The fixed axis or rotation axis is e, and the single rotation angle is denoted θ (Eq. 9).

⎡ qs ⎤ ⎡ q0 ⎤ ⎡ ⎛ ⎞⎤ ⎢ ⎥ ⎢ ⎥ ⎢ cos⎜ θ ⎟ ⎥ ⎝2⎠ ⎥ ⎡ s ⎤ ⎢ qx ⎥ ⎢ q1 ⎥ ⎢ q=⎣⎢ ⎥⎦ = ⎢ ⎥ = ⎢ ⎥ = qy v q2 ⎢ ⎛ θ ⎞⎥ ⎢ ⎥ ⎢ ⎥ ⎢e sin⎜ ⎟⎥ ⎝ 2 ⎠⎦ ⎢⎣ qz ⎥⎦ ⎣ q3 ⎦ ⎣

(9)

The three rotation angles (strike, dip and plunge) can be easily converted into a quaternion and a simple interpolation technique can be used to estimate at unsampled locations (Markley et al., 2007). The problem reduces to finding a 4  4 matrix M:

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182 n

M≜ ∑ wiqi qiT

(10)

i=1

where, wi are local weights determined using the inverse distance method. The desired estimated quaternion is the eigenvector of M (Eq. 11) corresponding to the maximum eigenvalue (Markley et al., 2007). Strike, dip and plunge are found locally by back transforming the estimated quaternion:

uαβϕ

⎡ 2 2 2 2⎤ ⎡ α ⎤ ⎢atan2( − 2q1q3 + 2q0q2, q3 − q2 − q1 − q0 ⎥ ⎥ = ⎢β⎥ = ⎢ atan2 2q2q3 + 2q1q0 ⎢ ⎥ ⎢ ⎥ ⎣ϕ⎦ ⎢ 2 2 2 2⎥ ⎣atan2( − 2q1q2 + 2q0q3, q2 − q3 − q0 − q1 ⎦

(

181

in 3D; the order of removing the border points is symmetric and therefore the skeleton is geometrically at the ‘middle’ of the original object. The implementation of Min (2011) is used. It is recommended that before applying thinning techniques the data are smoothed; MAPS (Deutsch, 1998) is used in this case to generate a more consistent and smooth centerline. In generating the LVA field, the centerlines/skeletons are fit with splines and the local orientations are found from the spline gradients (Fig. 16).

5. Conclusion

)

(11)

This method can generate a smooth interpolation of the three anisotropy parameters (Fig. 12) with a 3D example also provided (Figs. 13 and 14).

4. Data type 3: LVA field estimation from a geological solid model It is common to generate a geological solid/categorical model when there are different rock types, facies or depositional environments in a modeling domain. In many cases, the primary variable of interest follows along the features of this geological model. A morphological thinning technique is used to generate skeletons of the original geological model which are used to indicate local orientations of properties within the geological model. This technique is only appropriate when the variable of interest (grade, concentration, etc) is known to parallel the features of the geological model. Thinning algorithms are well known image processing techniques that have numerous applications in pattern recognition, data compression and data storage (Gonzalez and Woods, 2007). It is a morphological process that codes specified object pixels to the background in binary images and gives a compact shape representation (skeleton). Once boundary points have been identified, an object is eroded layer by layer until only a skeleton is left. The skeleton is the medial axes of an object. In Euclidean space the median axis is a set or loci of centers of maximal inscribed spheres (circles in 2D) (O’Rour, 1998). Maximal spheres are such that they are contained within the object but not entirely covered by another sphere. The exact computation of a medial axes is non-trivial (Wang and Basu, 2007, 2005). The methodology presented here is a template based iterative thinning technique that results in the approximate medial axis. The first criteria of thinning is that the object be binary mapped with the object having a value of 1 and the remainder of the model 0. The pixels containing zeros form the background and are in a stable condition, i.e. a point with zero value is never transformed to 1. Most thinning algorithms are parallel (a set of object voxels are coded to zero at the same time), but differ in the manner of organizing an iteration step. The 2D thinning algorithm K3M (Saeed et al., 2010), is a sequential iterative method that produces a thinned one voxel width skeleton and is used here. Geometry conservation is an important aim in all thinning algorithms; K3M (Saeed et al., 2010) removes the image pixels while preserving the image structure. It consists of a series of six phases which subsequently decide which pixels to keep and which to remove. The advantages of the algorithm are that it preserves the right angles at the line interconnections, which ultimately results in better correspondence between the original and the modified images (Fig. 15). The methodology proposed by Palagyi and Kuba (1999) is used

Second order stationarity is a critical decision in any geostatistical modeling methodology; the methodologies proposed herein can be used to generate parameters for many geostatistical methodologies that explicitly account for non-stationary features using LVA fields. The main conclusion of this work is that the generation of the LVA field is very data specific. Every deposit is unique and selecting an appropriate LVA generation methodology depends on the form of the information regarding orientation and magnitude of the LVA. Some caution should be exercised if these methods are applied to data that do not average linearly (such as permeability, grindability, etc.), such variables should be considered carefully. Although the generation of the LVA field is deposit specific and dependent on the type of available data, multiple techniques were presented that cover a wide range of potential sources for LVA field generation. If an exhaustive secondary variable is available and is deemed to have the same spatial features as the variable of interest, the gradient method with a moving window is recommended. When the scale of the features vary (even slightly) within a domain, it is important to consider the adaptive window technique. To clarify, this optimum moving window size does not relate to local kriging search windows. The kriging search window should relate to the range of correlation locally, rather than an optimum reliability which is used for determining local anisotropies; or the kriging neighbourhood should be restricted based on first order stationarity concerns. The difficulty in incorporating point axial data was examined and the implementation in both 2D and 3D was discussed for various cases. When the plunge angle can be ignored, the estimation of the orientations can be carried out by preprocessing the data and using any estimation/simulation technique. Quaternions are preferred when considering plunge but the problem formulation restricts the methodologies available for estimation of the quaternion components. If the only information available is in the form of a geologic interpretation, skeletonization of complex formations yields a compact shape representation that is consistent with the geological structure and can be used to infer the LVA field. This assumes that the LVA of the variable of interest follows along the modeled geological structures.

Acknowledgement The authors would like to thank member companies with the Center for Computational Geostatistics (CCG) for their financial support of this research.

Appendix Gen_lva generates an LVA field from 2D and 3D exhaustive data using gradient, PCA and moment of inertia (Hassanpour, 2007)

182

M. Lillah, J.B. Boisvert / Computers & Geosciences 82 (2015) 170–182

techniques. Additionally it has the option to use adaptive windows in determining LVA. If the adaptive window option is selected the user inputs the major range, ratios and angles of the ellipsoid window (Lines 16–21). Lines 9–11 define a standard GSLIB grid (Deutsch and Journel, 1998). The output file (Line 8) contains three angles and two ratios (Line 13–14) needed to describe the direction and magnitude of anisotropy. This is the same format that the suit of LVA programs presented in Boisvert and Deutsch (2011) requires. (Fig. 17) A standard GSLIB grid is used (Lines 7–9). Similar to gen_lva, this program allows the user to define the anisotropy magnitude for the domain (Lines 11–12) which are the necessary parameters to completely describe the LVA field. Search neighbourhoods are the traditional moving window whose dimension must be specified in Line 10 and must always be an even number to maintain a square search. Interpolation of direction angle measurements through quaternions is implemented using inverse distance using the desired power (ω) (Line 13). (Fig. 18) A GSLIB style program k3m (Fig. 19) performs skeletonization of 2D compact geological bodies. The boundaries of these structures must be well defined and the program expects that the bodies are generally well connected. Prior to thinning the geobody must be binary coded where all object nodes are assigned a value of 1. If there are irregularities such as missing data or jarring at the boundaries the object must be preprocessed with a smoothing program. Line 5 specifies the binary input file; line 6 is the output file which contains a 2D LVA field. In addition to thinning, k3m also cleans the centerline of any artefacts and fits smooth splines for computing the LVA orientation. The grid dimensions in line 7–8 and search window size in line 10 are required. An additional parameter (line 9) refers to the number of discrete nodes in the thinned skeleton that should be skipped during spline fitting and a user can select the ideal value by trial and error.

References Boisvert, J., 2010. Geostatistics with Locally Varying Anisotropy. University of Alberta, Edmonton, p. 175 (Ph.D. Dissertation). Boisvert, J., Deutsch, C., 2011. Programs for kriging and sequential Gaussian simulation with locally varying anisotropy using non-Euclidean distances. Comput. Geosci. 37 (4), 495–510. Cammoun, L., Castaño-Moraga, C.A., Muñoz-Moreno, E., Sosa-Cabrera, D., Acar, B., Rodriguez-Florido, M.A., Brun, A., Knutsson, H., Thiran, J.P., 2009. A review of tensors and tensor signal processing. Tensors in Image Processing and Computer Vision Editors: Assoc.Prof. Santiago Aja-Fernández, Assoc.Prof. Rodrigo de Luis García, Dr. Dacheng Tao, Dr. Xuelong Li ISBN: 978-1-84882-298-6 (Print)

978-1-84882-299-3 (Online). Cyganek, B., Siebert, J.P., 2009. An Introduction to 3D Computer Vision Technique and Algorithms. John Wiley & Sons Ltd., West Sussex, UK, p. 504. Deutsch, C., 1998. Cleaning categorical variable (lithofacies) realizations with maximum a-posteriori selection. Comput. Geosci. 24 (6), 551–562. Deutsch, C., Journel, A., 1998. Gslib: Geostatistical Software Library and User’s Guide. Oxford University Press, New York, p. 384. Deutsch, C., Lewis, R., 1992. Advances in the practical implementation of indicator geostatistics. In: Proceedings of the 23rd International APCOM Symposium. AZ. Publisher, Society of Mining Engineers, Tucson, pp. 169–179. Feng, X., Milanfar, P., 2002. Multiscale principal components analysis for image local orientation estimation. In: Proceedings of the 36th Asilomar Conference on Signals, Systems and Computers. Pacific Grove, CA, vol. 1, pp. 478–482. Grobekatthofer, K., Yoon, Z., 2012. Introduction into Quaternions for Spacecraft Attitude Representation. Technical University of Berlin, Berlin, Germany 〈http:// www.tu-berlin.de/fileadmin/fg169/miscellaneous/Quaternions.pdf〉. Hassanpour, M., 2007. Tools for Multivariate Modeling of Permeability Tensors and Geometric Parameters for Unstructured Grids (M.Sc. thesis). University of Alberta, Edmonton, p. 61. HauBecker, H., Jahne, B., 1996. A tensor approach for local structure analysis in multi-dimensional images. In: Proceedings of 3D Image Analysis and Synthesis ’96 (Proceedings in Artificial Intelligence 3), pp. 171–178. Hristopulos, D.T., 2002. New anisotropic covariance models and estimation of anisotropic parameters based on the covariance tensor identity. Stoch. Environ. Res. Risk Assess. 16 (1), 43–62. Katkovnik, V., Egiazarian, K., Astola, J., 2002. Adaptive window size image denoising based on Intersection of Confidence Intervals (ICI) rule. J. Math. Imaging Vis. 16, 223–235. Mardia, K., Jupp, P., 2000. Directional Statistics. John Wiley and Sons, West Sussex, England, p. 456. Markley, L.F., Cheng, Y., Crassidis, J., Oshman, Y., 2007. Quaternion veraging. J. Guid. Control Dyn. 30 (4), 1193–1196. McLennan, J.A., 2007. The Decision of Stationarity (Ph.D. thesis). University of Alberta, Edmonton, p. 191. Min, P., 2011. Thinvox version 0.35. URL: 〈http://www.cs.princeton.edu/  min/thin vox/〉. O’Rour, J., 1998. Computational Geometry in C, 2nd ed. Cambridge University Press, Cambridge, New York, p. 392. Palagyi, K., Kuba, A., 1999. Directional 3D Thinning Using 8 Subiterations. Discret. Geom. Comput. Imag. LNCS 1568, 325–336. Pardo-Igzquiza, E., Chica-Olmo, M., 2007. KRIGRADI: a cokriging program for estimating the gradient of spatial variables from sparse data. Comput. Geosci. 33 (4), 497–512 493. Saeed, K., Marek, T., Rybnik, M., Adamski, M., 2010. K3M: a Universal algorithm for image skeletonization and a review of thinning techniques. Int. J. Appl. Math. Comput. Sci. 20 (3), 317–335. Stroet, C., Snepvangers, J., 2005. Mapping curvilinear structures with local anisotropy kriging. Math. Geol. 37 (6), 635–649. Spiliopoulos, I., Hristopulos, D.T., Petrakis, M.P., Chorti, A., 2011. A multigrid method for the estimation of geometric anisotropy in environmental data from sensor networks. Comput. Geosci. 37 (3), 320–330. Wang, T., Basu, A., 2005. An Improved Fully Parallel 3D Thinning Algorithm. Department of Computing Science, University of Alberta, Canada. Wang, T., Basu, A., 2007. A note on ‘A fully parallel 3D thinning algorithm and its applications’. Pattern Recognit. Lett. 28, 501–506. Xu, W., 1996. Conditional curvilinear stochastic simulation using pixel-based algorithms. Math. Geol. 28 (7), 937–949.