A comparison of approaches to include outcrop information in overburden thickness estimation

The estimation of overburden sediment thickness is important in hydrogeology, geotechnics and geophysics. Usually, thickness is known precisely at a few sparse borehole data. To improve precision of estimation, one useful complementary information is the known position of outcrops. One intuitive approach is to code the outcrops as zero thickness data. A problem with this approach is that the outcrops are preferentially observed compared to other thickness information. This introduces a strong bias in the thickness estimation that kriging is not able to remove. We consider a new approach to incorporate point or surface outcrop information based on the use of a non-stationary covariance model in kriging. The non-stationary model is defined so as to restrict the distance of influence of the outcrops. Within this distance of influence, covariance parameters are assumed simple regular functions of the distance to the nearest outcrop. Outside the distance of influence of the outcrops, the thickness covariance is assumed stationary. The distance of influence is obtained thru a cross-validation. Compared to kriging based on a stationary model with or without zero thickness at outcrop locations, the non-stationary model provides more precise estimation, especially at points close to an outcrop. Moreover, the thickness map obtained with the non-stationary covariance model is more realistic since it forces the estimates to zero close to outcrops without the bias incurred when outcrops are simply treated as zero thickness in a stationary model.     

On a class of non-stationary, compactly supported spatial covariance functions

Globally supported covariance functions are generally associated with dense covariance matrices, meaning severe numerical problems in solution feasibility. These problems can be alleviated by considering methods yielding sparse covariance matrices. Indeed, having many zero entries in the covariance matrix can both greatly reduce computer storage requirements and the number of floating point operations needed in computation. Compactly supported covariance functions considerably reduce the computational burden of kriging, and allow the use of computationally efficient sparse matrix techniques, thus becoming a core aspect in spatial prediction when dealing with massive data sets. However, most of the work done in the context of compactly supported covariance functions has been carried out in the stationary context. This assumption is not generally met in practical and real problems, and there has been a growing recognition of the need for non-stationary spatial covariance functions in a variety of disciplines. In this paper we present a new class of non-stationary, compactly supported spatial covariance functions, which adapts a class of convolution-based flexible models to non-stationary situations. Some particular examples, computational issues, and connections with existing models are considered.     

基于阈值约束的协克里金法联合反演重力与重力梯度数据

常规协克里金方法反演重力或重力梯度数据具有抗噪性好、加入先验信息容易等优点,其反演的地下密度分布能够识别异常体中心位置,还原异常体基本形态,但反演图像光滑,分辨率低,这是由于常规方法估计的密度协方差矩阵全局发散、平稳.为了通过协克里金方法获得聚焦的密度分布需要改善密度协方差矩阵的性质.首先,本文推导了理论密度协方差公式,其性质表明,当理论模型聚焦分布时,其密度协方差矩阵是非平稳且聚焦分布的.为了打破常规协方差矩阵全局平稳、发散的特征,本文设置密度阈值处理协方差矩阵,通过不断更新协方差矩阵来迭代实现协克里金反演,最终得到相对聚焦的反演结果.用本文方法处理重力与重力梯度数据恢复两种密度模型,均得到了与正演模型匹配的反演结果;再将方法运用于文顿盐丘的实际测量重力与重力梯度数据,反演结果与已知的地质情况匹配较好.     

Plurigaussian modeling of geological domains based on the truncation of non-stationary Gaussian random fields

The plurigaussian model is used in mining engineering, oil reservoir characterization, hydrology and environmental sciences to simulate the layout of geological domains in the subsurface, while reproducing their spatial continuity and dependence relationships. However, this model is well-established only in the stationary case, when the spatial distribution of the domains is homogeneous in space, and suffers from theoretical and practical impediments in the non-stationary case. To overcome these limitations, this paper proposes extending the model to the truncation of intrinsic random fields of order k with Gaussian generalized increments, which allows reproducing spatial trends in the distribution of the geological domains. Methodological tools and algorithms are presented to infer the model parameters and to construct realizations of the geological domains conditioned to existing data. The proposal is illustrated with the simulation of rock type domains in an ore deposit in order to demonstrate its applicability. Despite the limited number of conditioning data, the results show a remarkable agreement between the simulated domains and the lithological model interpreted by geologists, while the conventional stationary plurigaussian model turns out to be unsuccessful.     

A non-stationary geostatistical approach to multigaussian kriging for local reserve estimation

Multigaussian kriging technique has many applications in mining, soil science, environmental science and other fields. Particularly, in the local reserve estimation of a mineral deposit, multigaussian kriging is employed to derive panel-wise tonnages by predicting conditional probability of block grades. Additionally, integration of a suitable change of support model is also required to estimate the functions of the variables with larger support than that of the samples. However, under the assumption of strict stationarity, the grade distributions and important recovery functions are estimated by multigaussian kriging using samples within a supposedly spatial homogeneous domain. Conventionally, the underlying random function model is required to be stationary in order to carry out the inference on ore grade distribution and relevant statistics. In reality, conventional stationary model often fails to represent complicated geological structure. Traditionally, the simple stationary model neither considers the obvious changes in local means and variances, nor is it able to replicate spatial continuity of the deposit and hence produces unreliable outcomes. This study deals with the theoretical design of a non-stationary multigaussian kriging model allowing change of support and its application in the mineral reserve estimation scenario. Local multivariate distributions are assumed here to be strictly stationary in the neighborhood of the panels. The local cumulative distribution function and related statistics with respect to the panels are estimated using a distance kernel approach. A rigorous investigation through simulation experiments is performed to analyze the relevance of the developed model followed by a case study on a copper deposit.     

Three-parameter Radon transform in layered transversely isotropic media

Radon transform is a powerful tool with many applications in different stages of seismic data processing, because of its capability to focus seismic events in the transform domain. Three-parameter Radon transform can optimally focus and separate different seismic events, if its basis functions accurately match the events. In anisotropic media, the conventional hyperbolic or shifted hyperbolic basis functions lose their accuracy and cannot preserve data fidelity, especially at large offsets. To address this issue, we propose an accurate traveltime approximation for transversely isotropic media with vertical symmetry axis, and derive two versions of Radon basis functions, time-variant and time-invariant. A time-variant basis function can be used in time domain Radon transform algorithms while a time-invariant version can be used in, generally more efficient, frequency domain algorithms. Comparing the time-variant and time-invariant Radon transform by the proposed basis functions, the time-invariant version can better focus different seismic events; it is also more accurate, especially in presence of vertical heterogeneity. However, the proposed time-invariant basis functions are suitable for a specific type of layered anisotropic media, known as factorized media. We test the proposed methods and illustrate successful applications of them for trace interpolation and coherent noise attenuation.     

A fully non-stationary linear coregionalization model for multivariate random fields

This paper introduces an extension of the traditional stationary linear coregionalization model to handle the lack of stationarity. Under the proposed model, coregionalization matrices are spatially dependent, and basic univariate spatial dependence structures are non-stationary. A parameter estimation procedure of the proposed non-stationary linear coregionalization model is developed under the local stationarity framework. The proposed estimation procedure is based on the method of moments and involves a matrix-valued local stationary variogram kernel estimator, a weighted local least squares method in combination with a kernel smoothing technique. Local parameter estimates are knitted together for prediction and simulation purposes. The proposed non-stationary multivariate spatial modeling approach is illustrated using two real bivariate data examples. Prediction performance comparison is carried out with the classical stationary multivariate spatial modeling approach. According to several criteria, the prediction performance of the proposed non-stationary multivariate spatial modeling approach appears to be significantly better.     

On a continuous spectral algorithm for simulating non-stationary Gaussian random fields

This paper presents an algorithm for simulating Gaussian random fields with zero mean and non-stationary covariance functions. The simulated field is obtained as a weighted sum of cosine waves with random frequencies and random phases, with weights that depend on the location-specific spectral density associated with the target non-stationary covariance. The applicability and accuracy of the algorithm are illustrated through synthetic examples, in which scalar and vector random fields with non-stationary Gaussian, exponential, Matérn or compactly-supported covariance models are simulated.     

New anisotropic covariance models and estimation of anisotropic parameters based on the covariance tensor identity

 Many heterogeneous media and environmental processes are statistically anisotropic. In this paper we focus on range anisotropy, that is, stochastic processes with variograms that have direction dependent correlation lengths and direction independent sill. We distinguish between two classes of anisotropic covariance models: Class (A) models are reducible to isotropic after rotation and rescaling operations. Class (B) models can be separated into a product of one-dimensional functions oriented along the principal axes. We propose a new Class (A) model with multiscale properties that has applications in subsurface hydrology. We also present a family of Class (B) models based on non-Euclidean distance metrics that are generated by superellipsoidal functions. Next, we propose a new method for determining the orientation of the principal axes and the degree of anisotropy, i.e., the ratio(s) of the correlation lengths. This information reduces the degrees of freedom of anisotropic variograms and thus simplifies the estimation procedure. In particular, Class (A) models are reduced to isotropic and Class (B) models to one-dimensional functions. Our method is based on an explicit relation between the second-rank slope tensor (SRST), which can be estimated from the data, and the covariance tensor. The procedure is conceptually simple and numerically efficient. It is more accurate for regular (on-grid) data distributions, but it can also be used for sparse (off-grid) spatial distributions. In the case of non-differentiable random fields the method can be extended using generalized derivatives. We illustrate its implementation with numerical simulations.     

A spectral approach to estimation and smoothing of continuous spatial processes

This paper is concerned with developing computational methods and approximations for maximum likelihood estimation and minimum mean square error smoothing of irregularly observed two-dimensional stationary spatial processes. The approximations are based on various Fourier expansions of the covariance function of the spatial process, expressed in terms of the inverse discrete Fourier transform of the spectral density function of the underlying spatial process. We assume that the underlying spatial process is governed by elliptic stochastic partial differential equations (SPDE's) driven by a Gaussian white noise process. SPDE's have often been used to model the underlying physical phenomenon and the elliptic SPDE's are generally associated with steady-state problems.A central problem in estimation of underlying model parameters is to identify the covariance function of the process. The cumbersome exact analytical calculation of the covariance function by inverting the spectral density function of the process, has commonly been used in the literature. The present work develops various Fourier approximations for the covariance function of the underlying process which are in easily computable form and allow easy application of Newton-type algorithms for maximum likelihood estimation of the model parameters. This work also develops an iterative search algorithm which combines the Gauss-Newton algorithm and a type of generalized expectation-maximization (EM) algorithm, namely expectation-conditional maximization (ECM) algorithm, for maximum likelihood estimation of the parameters.We analyze the accuracy of the covariance function approximations for the spatial autoregressive-moving average (ARMA) models analyzed in Vecchia (1988) and illustrate the performance of our iterative search algorithm in obtaining the maximum likelihood estimation of the model parameters on simulated and actual data.     

A latent class MDS model with spatial constraints for non-stationary spatial covariance estimation

Multidimensional scaling (MDS) has played an important role in non-stationary spatial covariance structure estimation and in analyzing the spatiotemporal processes underlying environmental studies. A combined cluster-MDS model, including geographical spatial constraints, has been previously proposed by the authors to address the estimation problem in oversampled domains in a least squares framework. In this paper is formulated a general latent class model with spatial constraints that, in a maximum likelihood framework, allows to partition the sample stations into classes and simultaneously to represent the cluster centers in a low-dimensional space, while the stations and clusters retain their spatial relationships. A model selection strategy is proposed to determine the number of latent classes and the dimensionality of the problem. Real and artificial data sets are analyzed to test the performance of the model.     

Scattering of elastic waves by a general anisotropic basin. Part 1: a 2D model

Scattering of incident plane harmonic pseudo P‐, SH‐, and SV‐waves by a two‐dimensional basin of arbitrary shape is investigated by using an indirect boundary integral equation approach. The basin and surrounding half‐space are assumed to be generally anisotropic, homogeneous, linearly elastic solids. No material symmetries are assumed. The unknown scattered waves are expressed as linear combinations of full‐space time‐harmonic two‐dimensional Green functions. Using the Radon transform, the Green functions are obtained in the form of finite integrals over a unit circle. An algorithm for the accurate and efficient numerical evaluation of the Green functions is discussed. A detailed convergence and parametric analysis of the problem is presented. Excellent agreement is obtained with isotropic results available in the literature. Steady‐state surface ground motion is presented for semi‐circular basins with generally anisotropic material properties. The results show that surface motion strongly depends upon the material properties of the basin as well as the angle of incidence and frequency of the incident wave. Significant mode conversion can be observed for general triclinic materials which are not present in isotropic models. Comparison with an isotropic basin response demonstrates that anisotropy is very important for assessing the nature of surface motion atop basins. Copyright © 2001 John Wiley & Sons, Ltd.     

轴对称各向异性介质波动方程深度偏移的对称非平稳相移方法

由所建立的三维qP波相速度表示式出发,导出并解析求解各向异性介质中的频散方程,得到三维各向异性介质中的相移算子,进而将以相移算子为基础的对称非平稳相移方法推广到各向异性介质,发展了一个三维各向异性介质的深度偏移方法. 文中使用的各向异性介质的速度模型与现行的各向异性构造的速度估计方法一致,将各向同性、弱各向异性及强各向异性统一在一个模型中. 所建立的各向异性介质对称非平稳相移波场延拓算子可以同时适应速度及各向异性参数横向变化;文中给出的算例虽然是针对二维VTI介质的,但所提出的算法同样适用于三维TI介质.     

Regularized Bayesian estimation for GEV-B-splines model

Large observed datasets are not stationary and/or depend on covariates, especially, in the case of extreme hydrometeorological variables. This causes the difficulty in estimation, using classical hydrological frequency analysis. A number of non-stationary models have been developed using linear or quadratic polynomial functions or B-splines functions to estimate the relationship between parameters and covariates. In this article, we propose regularised generalized extreme value model with B-splines (GEV-B-splines models) in a Bayesian framework to estimate quantiles. Regularisation is based on penalty and aims to favour parsimonious model especially in the case of large dimension space. Penalties are introduced in a Bayesian framework and the corresponding priors are detailed. Five penalties are considered and the corresponding priors are developed for comparison purpose as: Least absolute shrinkage and selection (Lasso and Ridge) and smoothing clipped absolute deviations (SCAD) methods (SCAD1, SCAD2 and SCAD3). Markov chain Monte Carlo (MCMC) algorithms have been developed for each model to estimate quantiles and their posterior distributions. Those approaches are tested and illustrated using simulated data with different sample sizes. A first simulation was made on polynomial B-splines functions in order to choose the most efficient model in terms of relative mean biais (RMB) and the relative mean-error (RME) criteria. A second simulation was performed with the SCAD1 penalty for sinusoidal dependence to illustrate the flexibility of the proposed approach. Results show clearly that the regularized approaches leads to a significant reduction of the bias and the mean square error, especially for small sample sizes (n < 100). A case study has been considered to model annual peak flows at Fort-Kent catchment with the total annual precipitations as covariates. The conditional quantile curves were given for the regularized and the maximum likelihood methods.     

基于保幅拉东变换的多次波衰减

为在去除多次波时有效保护地震一次反射波数据的AVO现象,给后续反演、解释提供准确的地震数据,本文提出了一种基于保幅拉东变换的多次波衰减方法,该方法是对常规抛物拉东变换的修改,把常规的稀疏拉东变换在拉东域分成两部分：一部分用于模拟零偏移距处的反射波能量,增加的另一部分用于模拟反射波振幅的AVO特性.该方法不仅考虑了反射波同相轴的形状,还考虑了反射波同相轴振幅幅度的变化,从而可把反射波信息进行有效转换,进而有利于多次波的消除,更好地恢复有效波的能量.在把地震数据由时间域转换到拉东域时,本文采用了IRLS算法实现保幅拉东算子的反演.模型数据和实际地震道集的试算分析表明,与常规拉东变换相比,保幅拉东变换在去除多次波的同时可有效保护一次反射波的AVO现象.     

Quantifying geological uncertainty for flow and transport modeling in multi-modal heterogeneous formations

In this work, we address the problem of characterizing the heterogeneity and uncertainty of hydraulic properties for complex geological settings. Hereby, we distinguish between two scales of heterogeneity, namely the hydrofacies structure and the intrafacies variability of the hydraulic properties. We employ multiple-point geostatistics to characterize the hydrofacies architecture. The multiple-point statistics are borrowed from a training image that is designed to reflect the prior geological conceptualization. The intrafacies variability of the hydraulic properties is represented using conventional two-point correlation methods, more precisely, spatial covariance models under a multi-Gaussian spatial law. We address the different levels and sources of uncertainty in characterizing the subsurface heterogeneity, and explore their effect on groundwater flow and transport predictions. Typically, uncertainty is assessed by way of many images, termed realizations, of a fixed statistical model. However, in many cases, sampling from a fixed stochastic model does not adequately represent the space of uncertainty. It neglects the uncertainty related to the selection of the stochastic model and the estimation of its input parameters. We acknowledge the uncertainty inherent in the definition of the prior conceptual model of aquifer architecture and in the estimation of global statistics, anisotropy, and correlation scales. Spatial bootstrap is used to assess the uncertainty of the unknown statistical parameters. As an illustrative example, we employ a synthetic field that represents a fluvial setting consisting of an interconnected network of channel sands embedded within finer-grained floodplain material. For this highly non-stationary setting we quantify the groundwater flow and transport model prediction uncertainty for various levels of hydrogeological uncertainty. Results indicate the importance of accurately describing the facies geometry, especially for transport predictions.     

时变地震谱估计的非参数方法评价及应用

针对工程中常见的非平稳地震动激励,进一步建立了时变功率谱估计的理论框架。首先,在多个样本情形下,采用信号处理中的时频分析方法对非平稳激励的时变谱估计理论进行了讨论。当激励样本数量有限时,对Priestley提出的估计方法进行了介绍。其次,以地震工程中常用的Kanai-Tajimi谱模型为目标,分别对均匀调制与一般调制谱的估计结果以及不同估计方法的精度与收敛性进行评价后,提出了能够方便工程应用的建议。最后,针对SMART-Ⅰ(Strong Motion Array in Taiwan,Phase Ⅰ)密集台阵第45次地震记录的多组三维样本,揭示了EW、NS、UD三个方向时变谱的典型特征以及各方向间的时变相干性。结果表明:与短时Fourier变换相比,复Morlet小波与广义谐和小波估计谱具有较好的精度与收敛性,Priestley方法估计单个样本具有优势;SMART-Ⅰ台阵地震动具有强度和频率的双重非平稳性,三个方向具有弱相干性。研究结论可为拓展谱估计理论的工程应用以及后续大跨结构多维多点非平稳地震响应的分析提供参考。     

Error covariance calculation for forecast bias estimation in hydrologic data assimilation

To date, an outstanding issue in hydrologic data assimilation is a proper way of dealing with forecast bias. A frequently used method to bypass this problem is to rescale the observations to the model climatology. While this approach improves the variability in the modeled soil wetness and discharge, it is not designed to correct the results for any bias. Alternatively, attempts have been made towards incorporating dynamic bias estimates into the assimilation algorithm. Persistent bias models are most often used to propagate the bias estimate, where the a priori forecast bias error covariance is calculated as a constant fraction of the unbiased a priori state error covariance. The latter approach is a simplification to the explicit propagation of the bias error covariance. The objective of this paper is to examine to which extent the choice for the propagation of the bias estimate and its error covariance influence the filter performance. An Observation System Simulation Experiment (OSSE) has been performed, in which ground water storage observations are assimilated into a biased conceptual hydrologic model. The magnitudes of the forecast bias and state error covariances are calibrated by optimizing the innovation statistics of groundwater storage. The obtained bias propagation models are found to be identical to persistent bias models. After calibration, both approaches for the estimation of the forecast bias error covariance lead to similar results, with a realistic attribution of error variances to the bias and state estimate, and significant reductions of the bias in both the estimates of groundwater storage and discharge. Overall, the results in this paper justify the use of the traditional approach for online bias estimation with a persistent bias model and a simplified forecast bias error covariance estimation.     

EMD Method Applied to Identification of Logging Sequence Strata

In this work, we compare Fourier transform, wavelet transform, and empirical mode decomposition (EMD), and point out that EMD method decomposes complex signal into a series of component functions through curves of local mean value. Each of Intrinsic Mode Functions (IMFs — component functions) contains all the information on the original signal. Therefore, it is more suitable for the interface identification of logging sequence strata.Well logging data reflect rich geological information and belong to non-linear and non-stationary signals and EMD method can deal with non-stationary and non-linear signals very well. By selecting sensitive parameters combination that reflects the regional geological structure and lithology, the combined parameter can be decomposed through EMD method to study the correlation and the physical meaning of each intrinsic mode function. Meanwhile, it identifies the stratigraphy and cycle sequence perfectly and provides an effective signal treatment method for sequence interface.     

Permissibility of fractal exponents and models of band-limited two-point functions for fGn and fBm random fields

The fractional Gaussian noise (fGn) and fractional Brownian motion (fBm) random field models have many applications in the environmental sciences. An issue of practical interest is the permissible range and the relations between different fractal exponents used to characterize these processes. Here we derive the bounds of the covariance exponent for fGn and the Hurst exponent for fBm based on the permissibility theorem by Bochner. We exploit the theoretical constraints on the spectral density to construct explicit two-point (covariance and structure) functions that are band-limited fractals with smooth cutoffs. Such functions are useful for modeling a gradual cutoff of power-law correlations. We also point out certain peculiarities of the relations between fractal exponents imposed by the mathematical bounds. Reliable estimation of the correlation and Hurst exponents typically requires measurements over a large range of scales (more than 3 orders of magnitude). For isotropic fractals and partially isotropic self-affine processes the dimensionality curse is partially lifted by estimating the exponent from measurements along fixed directions. We derive relations between the fractal exponents and the one-dimensional spectral density exponents, and we illustrate the relations using measurements of paper roughness.The author would like to acknowledge helpful comments from two anonymous referees.