A locally adaptive kernel regression method for facies delineation

1 A locally adaptive kernel regression method for facies delineation D. Fernàndez-Garcia, M. Barahona, C. Henri, X. Sanchez-Vila Hydrogeology Group (UPC-CSIC), Department of Geotechnical Engineering and Geosciences, Universitat Politècnica de Catalunya (UPC-BarcelonaTech), 08034 Barcelona, Spain Abstract Facies delineation is defined as the separation of geological units with distinct intrinsic characteristics (grain size, hydraulic conductivity, mineralogical composition). A major challenge in this area stems from the fact that only a few scattered pieces of hydrogeological information are available to delineate geological facies. Several methods to delineate facies are available in the literature, ranging from those based only on existing hard data, to those including secondary data or external knowledge about sedimentological patterns. This paper describes a methodology to use kernel regression methods as an effective tool for facies delineation. The method uses both the spatial and the actual sampled values to produce, for each individual hard data point, a locally adaptive steering kernel function, self-adjusting the principal directions of the local anisotropic kernels to the direction of highest local spatial correlation. The method is shown to outperform the nearest neighbor classification method in a number of synthetic aquifers whenever the available number of hard data is small and randomly distributed in space. In the case of exhaustive sampling, the steering kernel regression method converges to the true solution. Simulations ran in a suite of synthetic examples are used to explore the selection of kernel parameters in typical field settings. It is shown that, in practice, a rule of thumb can be used to obtain suboptimal results. The performance of the method is demonstrated to significantly improve when external information regarding facies proportions is incorporated. Remarkably, the method allows for a reasonable reconstruction of the facies connectivity patterns, shown in terms of breakthrough curves performance. 2 1. Introduction Image reconstruction has a long history in a number of disciplines such as satellite image mapping, shape recognition in robotics, face recognition, and license plate reading, among other uses [Bughin et al. 2008, Daoudi et al. 1999, Yang & Huang 1994, Lin & Chen 2007]. The topic can be loosely subdivided into two main groups: (a) The reconstruction of incomplete images where some of the pixels have no information; and (b) The reconstruction of noisy images, where some of the pixels display wrong information and the main problem is detecting and reclassifying the misclassified pixels. A good reconstruction work relies heavily on the presence of data and on an efficient reconstruction algorithm that can either complete information gaps, or else filter noisy signals. A particular case of reconstruction appears in subsurface hydrology, where the information relies on very few points (well logs), so that the initial available picture for reconstruction is mostly a black signal (meaning no information) with some sparse data scattered throughout the medium. Reconstruction is, thus, a really difficult and error prone task. Many methods for the interpolation of scattered data exist [Franke, 1982] and some of them have been used for geologic facies reconstruction [i.e., Ritzi et al., 1994, Guadagnini et al., 2004, Tartakovsky and Wohlberg, 2004, Wohlberg et al., 2006, Tartakovsky et al., 2007]. In particular, Tartakovsky et al. [2007] compared the fractional error obtained in two synthetic examples using three approaches: indicator kriging (IK) [Isaaks & Srivastava, 1990, Ritzi et al., 1994, Guadagnini et al., 2004], support vector machines (SVM) [Tartakovsky and Wohlberg 2004, Wohlberg et al., 2006] and nearest-neighbor classification (NNC) [Dixon, 2002]. Different sampling densities, ranging from 0.28% to 3.06%, and random sampling data generated following a 2D Poisson random process were used for comparison. Here sampling density refers to the proportion of pixels where hard data is available (pixels that are univocally classified). Their analysis indicated that NNC outperformed IK, in terms of proportion of correctly classified pixels, in both examples, and that SVM slightly outperformed NNC in one of the examples. There exist a number of reconstruction methods available in different disciplines that to our knowledge have never been used in geological facies reconstruction. A potential reason for this is that these methods were devised for the presence of massive data sets 3 that are never available in hydrogeology. One family of methods is based on kernel regression functions, widely used in signal theory for solving different problems such as image denoising, upscaling, interpolation, fusion, etc. Such methods have proved to be efficient for problems such as restoration and enhancement of noisy and/or incomplete sampled images. Even though regression methods have been used for reconstruction of images from extensive data sets, in principle, there is no reason not to use them when information is sparse. As an example, Takeda et al. [2007] tested a kernel regression method on an image reconstruction case in which only 15% of the pixels were informed, obtaining a very good reconstruction of a 2D image. Making an analogy between image reconstruction (from irregularly sampled data) and facies delineation (from scattered sampling points), we investigate the performance of a Steering Kernel Regression (SKR) method for the latter problem. The aim is to describe a methodology to use kernel regression as an effective tool for facies delineation, an application involving far less information available for image delineation from that for what it was originally developed (reconstruction). In doing this, we investigate the optimal tuning parameters to be used in the reconstruction of geological facies and their connectivity patterns. This paper is structured as follows; Section 2 briefly describes the fundamental concepts of facies reconstruction. Section 3 presents the details of the data-adapted kernel regression method. We test this method with respect to the NNC method in Section 4 by means of four synthetic images, here including the two figures profusely investigated by Tartakovsky et al. [2007] to allow for performance comparisons. 2. The concept of facies reconstruction The term facies is used in geology to differentiate among geological units on the basis of interpretive or descriptive characteristics, such as sedimentological conditions of formation, mineralogical composition, presence of fossils (biofacies), structures, grain size, etc. [Tarbuck et al., 2002]. In this work, we consider that each facies is a clear distinctive geology unit, understood in a descriptive sense. Keeping this in mind, facies reconstruction is defined as the process of assigning each unsampled point (eventually also the sampled ones if misclassification errors are admitted) to one facies. Formally, 4 for any given facies Fk, the reconstruction problem can be addressed using an indicator function defined as k k 1 (, ) 0 F IF otherwise ∈  =  x x (1) where the indicator variable I(x,Fk) is equal to 1 when a particular point in the domain, x, can be classified as belonging to facies Fk and zero otherwise. In this work we assume that the available data from the sampling points are clearly distinctive in order to be unmistakably classified as indicated in (1) without interpretation errors. From now on, we consider that only two facies are used for geological mapping. However, the method can be easily extended to any finite number of facies by direct superposition. Several methods have been proposed in the literature to estimate the spatial distribution of the indicator variable I(x,F1). Here we compile only three of such methods. The first one is indicator kriging (IK) [Journel, 1983], a method that provides a least-squares estimate of the probability that x belongs to F1 conditioned to nearby data. Once a threshold value is given, a distinction between categories (facies) can be done. The method relies on the theory of random functions to model the uncertainty of not having data at unknown locations. It accounts for the inherent spatial correlation of data but typically fails to properly estimate curvilinear geological bodies. Multiple point geostatistics [e.g., Strebelle, 2000] can overcome most of these problems by largely relying on an empirical multivariate distribution inferred from training images, i.e., under the assumption that significant information about the spatial distribution of facies is known from external sources (outcrops, modeling of sedimentological processes,…); these information is directly transferred to the final images. Alternatively, Support Vector Machine (SVM) methods are a set of popular tools for data mining tasks such as classification, regression, and novelty detection [Vapnik, 1963; Bennett and Campbell, 2000]. SVM takes a training data, i.e., a set of n data points Ji= J(xi,F1)∈{-1,1}, i=1,..,n, and separates them into two classes by delineating the hyperplane that has the largest distance to the nearest training data point of any class. Last, the nearest-neighbor classification (NNC) simply classifies each point in the domain by finding the nearest (not necessarily in the Euclidean sense) training point, assigning to the unsampled location the class corresponding to that training point. 5 A comparison of the three methods presented is provided in a series of papers by Tartakovsky and Wholberg [2004], Wholberg et al. [2006], and Tartakovsky et al. [2007]. Surprisingly, the NNC method outperformed the more sophisticated ones, i.e., SVM and IK, indicating the validity of the parsimony principle for this problem. Yet, the comparison between methods in such works was done only in terms of the number of misclassified points without considering other performance metrics, such as connectivity features inherent in geological facies that can strongly impact contaminant transport simulations (e.g., Fernàndez-Garcia et al., 2010). We consider this issue as non-ideal and in the next section we seek for a method that can actually represent the presence of connected geological bodies with elongated and curvilinear shapes. 3. Kernel regression approaches for facies classification Kernel regression methods have been developed in statistics to estimate the conditional expectation of a random variable without assumptions about its probability distribution function. These methods are well documented and summarized in the literature [e.g., Hardle, 1990; Simonoff, 1996; Li et al., 2007]. Suppose that we ignore the fact that the target classification output is a binary function I(x,F1). Instead, we consider that it is a continuous function that depends on the location x and a number of (yet unknown) parameters b=[b0,b1,…,bN]T. The regression model proposed here for facies classification assumes that the measured data Ii=I(xi,F1), i=1,…,n, can be expressed as iii mI ε + =); (bx , n i,..,1= , (2) where m(xi,b) is the regression function to be determined, and ε i are independent and identically distributed zero mean noise values. Kernel regression is a form of regression analysis in which the function m is exclusively dictated by the data, and not prespecified a priori (no model assumed). At each point x the conditional expected value of the dependent (indicator) variable can be estimated, i.e., m(x,b)=E[I(x,F1)]. The interest of kernel regression to facies reconstruction resides on the fact that the conditional expected value of the indicator variable is exactly the probability that the given facies F1 prevails at that location, since { } { } { } { } 1 1 11 ( , ) 1 Prob 0 Prob ProbEI F F F F=⋅∈+⋅∉=∈x x xx (3) 6 By definition, the probability of occurrence of a given facies is a continuous variable ranging between 0 and 1. In order to separate the data into classes or facies we must then establish a cut-off in the estimate of the indicator variable. This is similar to the facies reconstruction problem posed by the geostatistical indicator kriging approach. In this case, Ritzi et al. [1994] has suggested to define the boundary between facies by the isoline Prob{x∈Fk}=pk, where pk is estimated as either the global mean of the indicator values or the empirical relative volumetric fraction of the facies Fk. We propose here to use the same approach for classifying facies with regression methods. The benefits of such approach will be explored in section 4. Two kernel regression methods, namely the classical (CKR) and the adaptive steering (SKR) are presented next, and later their performance is compared in a number of synthetic cases. 3.1. Classical kernel regression (CKR) Let us consider a local Taylor expansion of the mean response m(x,b) of the indicator values around the estimation location x0, 22 0 01 2 3 4 5 6 7 ( ; ) ( ; , ) ' ' ' ' ' ' ' ' '...m m bbxbybzbx bxyby bxz≈ =++++ + + +xb xbx (4) where x’=x-x0 is the distance between any point and that being estimated, b0 is the mean response at x0, [b1,b2,b3]T is the gradient of the mean response at x0, and so on. The order of the polynomial is in principle arbitrary. Nonparametric regression generalizes the standard regression approach by locally estimating b at a given location x0 using only nearby data. This is done by weighting data located far away from the estimation location with a kernel function KH defined as ( ) xH H x 1 )det( 1 ) ( − =K K H (5) where H is a matrix that controls the degree of smoothing and is user dependent. The kernel associates a very low weight to points located far from the estimation point. Section 4 will explore the choice of kernel parameters for optimal facies reconstruction. The kernel function K is a continuous, bounded, and symmetric real function centered at zero that integrates to one and typically decays with distance. The choice of the kernel is known to not affect significantly the final solution and therefore a standard Gaussian 7 distribution is typically used for mathematical convenience. In n dimensions this is written as () /2 11 ( ) exp 2 2 T n K p  =  x xx (6) For any given estimation location x0, the principle of least squares expresses that one should choose as estimates of b those values that minimize the weighted sum of squared residuals, S(b), the residual being the difference between data values and model predictions, (7) Let us express equation (2) in matrix form, eXbI+= (8) where I=[I1,..,In]T, e=[ ε1, …,εn]T, and is a matrix composed of n rows and a number of columns that is associated with the degree of the polynomial chosen for b (i.e., in 3-D would be 4 for order 1, 10 for order 2,…) ' ' ' '2 ' ' '2 ' ' 111111111 ' ' ' '2 ' ' '2 ' ' 1 ... ... ... ... ... ... 1 ... nnnn nnn nn x y z x xy y xz x y z x xy y xz   =   X (9) Then, the optimization problem is written as (10) where W is a diagonal weight matrix given by { } )(K),...,(Kdiag 0nH01H xxxxW −−= (11) Setting ∂S(b)/∂bj=0 to each parameter bj we obtain the following solution ( ) WIX WXXb TT 1 ˆ − = (12) This solution is formally the same to that of standard regression but the matrices W and X depend now on the estimation location x0. Knowing the optimal estimate of b, the probability that x belongs to F1 can be estimated by 8 { } { } 1 1 010 0 ˆ Prob ( , ) ( , , )F EI F m F b∈= = =x xx xx (13) Let us define Weq by ( ) WXWXXeW TTT eq 1 1 − = (14) where e1 is a column vector with first element equal to one, and the rest equal to zero. Then, the Classical Kernel Regression (CKR) algorithm can be seen as a local weighted averaging of the data in which the probability that x belongs to F1 is determined by the following linear interpolation of indicator values 0 ˆT eq b= ⋅WI (15) Hence, Weq is a vector containing the equivalent weights of the indicator data values. The forms of these equivalent weights are exclusively dictated by the polynomial order chosen in (4). 3.2. Steering kernel regression (SKR) The SKR method comes as a direct extension of the CKR algorithm. Since the latter is nothing but a weighted average of indicator data values, the final regression estimate of Prob{x∈F1} only depends on the geometric configuration of the data, and therefore ignores the inherent correlations between data positions and their values. Takeda et al. [2007] developed a SKR algorithm to include key structural features into the estimated fields. The key idea behind the SKR algorithm is to modify the size and orientation of the regression kernel to assign more weight along the direction of highest local spatial correlation. The advantage of doing this to classify facies is the following: consider a point x0∈F1 located close to a facies boundary; the conventional CKR algorithm (symmetric spherical kernel) will estimate the probability that x0 belongs to F1 by equally considering both nearby samples of the same facies F1 and samples of other facies located beyond the boundaries. The SKR method is designed to adapt the regression kernel to the boundary isosurface so as to assign more weight to those samples belonging to the same facies. This way, the denoising is affected most strongly along the boundaries, rather than across them, resulting in a strong preservation of details in the final output. 9 The algorithm works by reorienting the smoothing matrix in the direction of the gradients of the mean response m(x,b) through a redefinition of the kernel matrix 2/1− = i steer i hCH (16) () ) ˆ ,( ) ˆ , (b x bx C j T j i mm ∇⋅ ∇ ≈ , ij w∈x (17) where the overbar stands for averaging over the mean response adjacent to xi, wi is the window search around xi, and h is a global smoothing parameter. In contrast to the CKR algorithm, the smoothing matrix Hsteer at each individual point xi depends now on the solution of the regression function m(x,F1). This makes the SKR method to be nonlinear in nature. Its application must be therefore iterative, starting with a first initial estimate of m(x,F1) computed, for instance, with the CKR method. This estimate is used to measure the dominant orientation of the local gradients, then used to sequentially steer the local kernel function through (17), resulting in elongated, ellipsoidal contours spread along the indicator isosurface (or isocurve in 2D). We must state that while the method is applicable to 3D reconstruction problems, here we present the details only for the 2D problems. The main reason is to be able to use the same synthetic examples available in the literature for geologic facies reconstruction using IK, SVM or NNC methods. Under these conditions, and from (16), the new form of the regression kernel is ( ) ( )        −− =− 2 0 0 0 2 exp 2 ) det( )( h h K ii T i i i Hsteer i xx Cxx C xx p (18) When estimating the covariance matrix Ci through (17), the resulting matrix can be rank deficient and unstable. To overcome this problem, a multiscale technique for estimating local gradients [Takeda et al., 2007] can be adopted. Let us consider the following matrix Gi formed by a collection of p estimated gradient vectors at the neighborhood of the sampled location xi           ∇ ∇ = ),( ... ),( 1 bx bx G p i m m , ,...,pjw ij 1 , =∈x (19) The singular value decomposition of Gi factorizes this matrix in the following form 16 correlation. As such, its reconstructed image (Figure 6b) fails to represent the central spatial continuity observed in facies F1, clearly extending from the northern to the southern boundaries. Instead, with only four iterations, the SKR(%) is able to correctly identify this spatial continuity of data values and properly represent the true connection north-south. Figure 3 illustrates the evolution of the local kernel functions associated to each data point in the same problem. In these images, the variable represented is the direct output data given by the SKR method without applying a classification strategy, and the progressive increase in the ratio of the two axes of the ellipse can be observed. In addition to the recognition of spatial continuity, the SKR(%) mehod is also capable of providing a measure of uncertainty in the delineation of the facies boundary. In principle this is not possible for any deterministic approach, such as that of the NNC algorithm. Figure 7 presents different maps to evaluate the uncertainty in the estimation corresponding to the same example already used previously. Interestingly, there is a very good correlation between low variance and high sampling density areas and viceversa. From this map, one can also delineate a safe zone for drawing the border between facies, plotted as gray areas in Figure 7e, those corresponding to values above 0.3 times the standard deviation. By visual inspection, a very good agreement between the results from the SKR(%) method (Figure 7e) and the original facies boundaries visible in Figure 7a, can be appreciated. 4.4. Impact on transport predictions In this paper, we contend that a key aspect to consider during the reconstruction of geological facies is the representation of connectivity. Even though the SKR method is shown to only slightly outperform the NNC in terms of volumetric fractional errors (see Figure 5), results demonstrated that the NNC is often not capable to properly describe the spatial continuity of the facies body. Solute transport simulations in a Monte Carlo framework were further performed to illustrate the impact that this effect can have on contaminant transport predictions. To do this, we considered the synthetic field presented in Figure 1a as a reference geological setting. The hydraulic conductivity is assumed to vary in space. A hydraulic conductivity of K=100 m/day and K=1 m/day was respectively assigned to the blue and red facies. Figure 8 shows the setup of the simulations. A non-reactive contaminant source was assumed to be originally located in a southern block region of size 5×5 m2 (it is assumed that each pixel has a length of 1 m). Groundwater is assumed at steady-state and moves from south to north along the 17 main facies direction driven by a hydraulic gradient of 0.001 in the x-direction and 0.002 in the y-direction. Prescribed heads are fixed at all boundaries according to this hydraulic gradient. Solute transport was simulated with a random walk code that solves the advectiondispersion equation [Fernàndez-Garcia et al., 2005; Henri and Fernàndez-Garcia, 2014]. Transport parameters were considered constant with a porosity of 0.3, a longitudinal dispersivity of 0.1 m, and a transverse dispersivity of 0.01 m. The effect of heterogeneity inside each facies was not considered to only focus on the reconstruction problem. Contaminant concentrations were observed at a control plane located at y=5 m. We then compare the transport simulations obtained with the reference hydraulic conductivity field with those resulting from the one hundred SKR(%) and NNC realizations generated using a sample density of 30 data points. The cumulative breakthrough curves are shown in Figure 9 normalized by the total mass injected. The ensemble of solutions provided by the SKR(%) and NNC methods is represented by the median and the 95% confidence interval (yellow region in this figure). Individual realizations are also depicted. Results clearly show that the SKR is more robust than the NNC method in terms of transport predictions. Even though the median solution provided by both methods is close to the true solution, the confidence interval of the SKR method is strikingly smaller than that obtained by the NNC, an effect that is more pronounced at late times. This indicates that the probability that reality is not properly represented by the SKR method is substantially smaller. Remarkably, this also reflects that, in many of the NNC realizations, the contaminant is forced to move through inexistent small permeability areas, resulting in artificial tailing and an artificial retardation. We also note that in some realizations, the poor south-north connection described by the NNC is such that the contaminant is partially exiting the system from the east and west boundaries without reaching at the control plane (note that some breakthrough curves do not contain all the mass injected). 5. Conclusions A non-parametric method, SKR, originally designed for image processing [Takeda et al. 2007], has been presented and tested for its application as a facies delineation algorithm. The performance of the method was compared with the nearest neighbor classification, 18 a method that has proven to be more efficient than others discussed in the literature [Tartakovsky et al., 2007]. Four synthetic scenarios were used for the comparison: two of them identical to the figures presented by Tartakovsky et al. [2007], and the other two figures are new for this work, one inspired on a cartographied river meander, and the other being a representation of a simple geometry. For each example different tests were studied ranging from very sparse to sparse number of data points available. Two variations of the SKR method were tested depending on whether additional information about the exact proportion of facies was introduced in the algorithm (SKR(%)) or not (SKR(0)). Our results indicate that the SKR(0) method had similar or lower fractional errors than those obtained with NNC, except for two cases (Figure 1(c) and (d), with a sampling density of 0.28%). The SKR(%) outperformed all methods, with improvements up to 5% in terms of reduction in misclassified points. The improvement is better in relative terms for the lowest sampling densities. This finding leads us to believe that the SKR(%) method would be an useful tool on real cases, when scattered and few sampling data points are expected. One of the major advantages of the SKR method is the quantification of the uncertainty in the delineation of the facies boundaries. In this context, we presented a method to stochastically generate variance maps that allows one to identify potential areas where a boundary between facies is more likely to exist. An example of application for one of the study cases is provided, leading to the delineation of an area over which there is most probably a boundary between facies. Acknowledgements This work has been supported by the Spanish Ministry of Science and Innovation through projects Consolider-Ingenio 2010 CSD2009-00065 and FEAR CGL2012-38120. XS acknowledges support of Program ICREA Acadèmia. Appendix: The nearest-neighbor classification (NNC) The nearest-neighbor classification (NNC) employed by Tartakovsky et al. [2007] is a k-nearest-neighbor classification [Hastie et al., 2001] in which the classification of a test point is determined by majority vote amongst the k nearest-neighbor points in the 19 training set, Tartakovsky et al. [2007] considered the case in which k=1, for which the classification of each point in the domain is determined by finding the nearest training point, and assigning the known class of that point. Given a set of training data points Ii=I(xi, Fk), i=1,…,n, the NNC classification for an arbitrary point x in the domain is computed as follows: (1) Define j as the index of the training data point, from the set { } 1 N ii x= , which is closest to query point x ; that is, 2 argmin ii j xx = − . Usually an Euclidean measure is prefer as distance metric, for simplicity, however, other metric can be used; (2) Assign the indicator function value of training data point j x (i.e., () j Ix ) as the indicator function value of query point x . This classification is simple to compute, and has no free parameters to estimate (no optimization of the method is possible). References Bennett, K.P., Campbell, C., 2000. Support vector machines: Hype or hallelujah? SIGKDD Explorations, 2(2). Bughin, E., Blanc-Feraud, L., Zerubia, J., 2008. Satellite image reconstruction from an irregular sampling. Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). 849-852. doi: 10.1109/ICASSP.2008.4517743 Daoudi, M., Ghorbel, F., Mokadem, A., Avaro, O., Sanson, H., 1999. Shape distances for contour tracking and motion estimation. Pattern Recognition, 1297-1306. Dixon, P.M., 2002. Nearest neighbour methods, in Encyclopedia of Environmetrics, vol. 3, edited by A. H. El-Shaarawi andW.W. Piegorsch, pp. 1370– 1383, John Wiley, New York. Guadagnini, L., Guadagnini, A., Tartakovsky, D.M., 2004. Probabilistic reconstruction of geologic facies. Journal of Hydrology, 294, 57-67. Feng, X., Milanfar, P., 2002. Multiscale principal components analysis for image local orientation estimation. 36th Asilomar Conf. Signals, Systems and Computers. 20 Fernàndez-Garcia, D., Illangasekare, T. H., Rajaram, H., 2005. Differences in the scaledependence of dispersivity estimated from temporal and spatial moments in physically and chemically heterogeneous porous media, Adv. Water Res., 28, 745-759, 2005. Fernàndez-Garcia, D., Trinchero, P., Sanchez-Vila, X., 2010. Conditional stochastic mapping of transport connectivity, Water Resour. Res., 46, Art. No. W10515. Franke, R., 1982. Scattered Data Interpolation: Tests of Some Methods. Mathematics of Computation, 38(157), 181-200. Hardle, W., 1990. Applied nonparametric regression, Econometric Society Monographs No. 19, Cambridge University Press, 333 p., Henri, C. V., Fernàndez-Garcia D., 2014. Toward efficiency in heterogeneous multispecies reactive transport modeling: A particle-tracking solution for firstorder network reactions, Water Resour. Res., 50, doi:10.1002/2013WR014956. Isaaks, E. H., Srivastava, R. M., 1990. An Introduction to Applied Geostatistics, Oxford Univ. Press, New York. Journel, A.G., 1983. Non-parametric estimation of spatial distribution. Mathematical Geology, 15(3). Li, Qi; Racine, Jeffrey S., 2007. Nonparametric Econometrics: Theory and Practice. Princeton University Press. Lin, S.C., Chen, C.T., 2008. Reconstructing vehicle license plate Image from low resolution images using nonuniform interpolation method. International Journal of Image Processing, 1(2). Ollero Ojeda, A., 1996. Dinámica de meandros y riesgos hidrogeomorfológicos en Alcalá de Ebro y Cabañas de Ebro (Zaragoza). In: IV Reunión de 21 Geomorfología, Grandal d’Anglade, A. and Pagés Valcarlos, J. Eds. Sociedad Española de Geomorfología. Ritzi, R.W. Jr., Jayne, D.F., Zahradnik, A.J.Jr., Field, A.A., Fogg, G.E., 1994. Geostatistical Modeling of Heterogeneity in Glaciofluvial, Buried-Valley Aquifers. Ground Water 32(4), 666-674. Simonoff, J.S., 1996. Smoothing Methods in Statistics. Springer. Strebelle, S., 2000. Sequential simulation drawing structures from training images: Unpublished doctoral dissertation, Stanford University, 200 p. Takeda, H., Farsiu, S., Milanfar, P., 2007. Kernel Regression for Image Processing and Reconstruction. IEEE Transactions on Image Processing, Vol. 16, No. 2, February 2007. Tarbuck, E.J., Lutgens, F.K., 2002. Earth: An Introduction to Physical Geology. Seventh Edition. Prentice Hall Tartakovsky, D.M., Wohlberg, B.E., 2004. Delineation of geologic facies with statistical learning theory. Geophysical Research Letters, Vol. 31, L18502, doi:10.1029/2004GL020864 Tartakovsky, D.M., Wohlberg, B., Guadagnini, A., 2007. Nearest-neighbor classification for facies delineation. Water Resour. Res., 43, W07201, doi: 10.1029/2007WR005968 Vapnik, V., Lerner, A., 1963. Pattern recognition using generalized portrait method. Automation and Remote Control, 24, 774–780. Wohlberg, B., Tartakovsky D.M., Guadagnini A., 2006. Subsurface characterization with support vector machines. IEEE Transactions on Geoscience and Remote Sensing, Vol. 44, No. 1, January 2006. 22 Yang, G., Huang, T.S., 1994. Human face detection in a complex background. Pattern Recognition. Volume 27, Issue 1, January 1994, Pages 53–63. 23 Figure 1. Synthetic fields used for facies delineation: a and b are the same figures presented by Tartakovsky et al. [2007]. We generated Figure 1 (c) and (d) considering a real case scenario (a meander from the Ebro river, Spain), and a simple geometric figure (circle). Blue and red colors indicate the two distinct facies. 24 Figure 2. Sensitivity analysis of the parameters needed to reconstruct geological facies with the SKR method: the local orientation analysis window ( w ), the regularization for the elongation parameter ( λ ), the structure sensitive parameter ( α ) and the global smoothing parameter (h). Blue dots indicate the different value choices for the calculation of the fractional errors and the red star indicates the value used for our calculations, coincidently with the lowest fractional error. 25 Figure 3. Iteration comparison: a) Original figure corresponding to Figure 1a. Random sampling points are shown as blue and red squares (example with a sample density of 30); b) Classical Kernel Regression results. The first, second and third iteration of the Steering Kernel is shown in c), d) and e), respectively.