Geostatistical traveltime tomography in elliptically anisotropic media

In geological materials, anisotropy may arise due to different mechanisms and can be found at different scales. Neglecting anisotropy in traveltime tomographic reconstruction leads to artefacts that can obscure important subsurface features. In this paper, a geostatistical tomography algorithm to invert cross‐hole traveltime data in elliptically anisotropic media is presented. The advantages of geostatistical tomography are that the solution is regularized by the covariance of the model parameters, that known model parameters can be used as constraints and fitted exactly or within a prescribed variance and that stochastic simulations can be performed to appraise the variability of the solution space. The benefits of the algorithm to image anisotropic media are illustrated by two examples using synthetic georadar data and real seismic data. The first example confirms suspected electromagnetic anisotropy in the vadose zone caused by relatively rapid water content variations with respect to wavelength at georadar frequencies. The second presents how sonic log data can be used to constrain the inversion of cross‐well seismic data and how geostatistical simulations can be used to infer parameter uncertainty. Results of both examples show that considering anisotropy yields a better fit to the data at high ray angles and reduces reconstruction artefacts.     

Target-level waveform inversion: a prospective application of the convolution-type representation for the acoustic wavefield

Nowadays, full-waveform inversion, based on fitting the measured surface data with modelled data, has become the preferred approach to recover detailed physical parameters from the subsurface. However, its application is computationally expensive for large inversion domains. Furthermore, when the subsurface has a complex geological setting, the inversion process requires an appropriate pre-conditioning scheme to retrieve the medium parameters for the desired target area in a reliable manner. One way of dealing with both aspects is by waveform inversion schemes in a target-oriented fashion. Therefore, we propose a prospective application of the convolution-type representation for the acoustic wavefield in the frequency–space domain formulated as a target-oriented waveform inversion method. Our approach aims at matching the observed and modelled upgoing wavefields at a target depth level in the subsurface, where the seismic wavefields, generated by sources distributed above this level, are available. The forward modelling is performed by combining the convolution-type representation for the acoustic wavefield with solving the two-way acoustic wave-equation in the frequency–space domain for the target area. We evaluate the effectiveness of our inversion method by comparing it with the full-domain full-waveform inversion process through some numerical examples using synthetic data from a horizontal well acquisition geometry, where the sources are located at the surface and the receivers are located along a horizontal well at the target level. Our proposed inversion method requires less computational effort and, for this particular acquisition, it has proven to provide more accurate estimates of the target zone below a complex overburden compared to both full-domain full-waveform inversion process and local full-waveform inversion after applying interferometry by multidimensional deconvolution to get local-impulse responses.     

Towards large-scale stochastic refraction tomography: a comparison of three evolutionary algorithms

The main goal of this study is to assess the potential of evolutionary algorithms to solve highly non-linear and multi-modal tomography problems (such as first arrival traveltime tomography) and their abilities to estimate reliable uncertainties. Classical tomography methods apply derivative-based optimization algorithms that require the user to determine the value of several parameters (such as regularization level and initial model) prior to the inversion as they strongly affect the final inverted model. In addition, derivative-based methods only perform a local search dependent on the chosen starting model. Global optimization methods based on Markov chain Monte Carlo that thoroughly sample the model parameter space are theoretically insensitive to the initial model but turn out to be computationally expensive. Evolutionary algorithms are population-based global optimization methods and are thus intrinsically parallel, allowing these algorithms to fully handle available computer resources. We apply three evolutionary algorithms to solve a refraction traveltime tomography problem, namely the differential evolution, the competitive particle swarm optimization and the covariance matrix adaptation–evolution strategy. We apply these methodologies on a smoothed version of the Marmousi velocity model and compare their performances in terms of optimization and estimates of uncertainty. By performing scalability and statistical analysis over the results obtained with several runs, we assess the benefits and shortcomings of each algorithm.     

Joint inversion of surface wave dispersion curves and reflection travel times via multi-objective evolutionary algorithms

Due to the character of seismic energy generation and propagation, shallow high-resolution seismic-reflection surveys often fail in the identification of the shallowest horizons and, due to the limited offsets, accuracy of velocity analyses is often not very high.In recent years, Rayleigh wave dispersion analysis have proved to have good potential also for near-surface applications but dispersion curve inversion and related uncertainty evaluation pose serious problems to a completely stand-alone application.In order to overcome these problems a joint inversion scheme is proposed, which is based on the identification of the Pareto front, performed in the framework of a Multi-Objective Evolutionary Algorithm (MOEA). Seismic data considered to design the two objectives are the Rayleigh wave dispersion curve and reflection travel times.We initially analyse a set of synthetic cases and evaluate the obtained results. A significant improvement of the retrieved models is observed as long as reflection travel times are added to the dispersion curve alone.Furthermore, the proposed methodology also provides relevant indications about the consistency of the overall inversion process. In fact, the distribution of the models in the objective space, the trend of the objectives over the passing generations and the evolution of the Pareto front can provide useful information to evaluate the provisional tentative interpretation (number of strata and reflector identification) inherently adopted for the data inversion.On the basis of the results obtained from the tests on the synthetic datasets, the analyses of a field dataset are interpreted as possible evidence of lateral heterogeneities.     

A genetic algorithm for filter design to enhance features in seismic images

We present a novel method to enhance seismic data for manual and automatic interpretation. We use a genetic algorithm to optimize a kernel that, when convolved with the seismic image, appears to enhance the internal characteristics of salt bodies and the sub‐salt stratigraphy. The performance of the genetic algorithm was validated by the use of test images prior to its application on the seismic data. We present the evolution of the resulting kernel and its convolved image. This image was analysed by a seismic interpreter, highlighting possible advantages over the original one. The effects of the kernel were also subject to an automatic interpretation technique based on principal component analysis. Statistical comparison of these results with those from the original image, by means of the Mann‐Whitney U‐test, proved the convolved image to be more appropriate for automatic interpretation.     

Estimation of a non-stationary prior covariance from seismic data

Non-stationarity in statistical properties of the subsurface is often ignored. In a classical linear Bayesian inversion setting of seismic data, the prior distribution of physical parameters is often assumed to be stationary. Here we propose a new method of handling non-stationarity in the variance of physical parameters in seismic data. We propose to infer the model variance prior to inversion using maximum likelihood estimators in a sliding window approach. A traditional, and a localized shrinkage estimator is defined for inferring the prior model variance. The estimators are assessed in a synthetic base case with heterogeneous variance of the acoustic impedance in a zero-offset seismic cross section. Subsequently, this data is inverted for acoustic impedance using a non-stationary model set up with the inferred variances. Results indicate that prediction as well as posterior resolution is greatly improved using the non-stationary model compared with a common prior model with stationary variance. The localized shrinkage predictor is shown to be slightly more robust than the traditional estimator in terms of amplitude differences in the variance of acoustic impedance and size of local neighbourhood. Finally, we apply the methodology to a real data set from the North Sea basin. Inversion results show a more realistic posterior model than using a conventional approach with stationary variance.     

Constraining the anisotropy structure of the crust by joint inversion of seismic reflection travel times and wave polarizations

In the context of wide-angle seismic profiling, the determination of the physical properties of the Earth crust, such as the elastic layer depth and seismic velocity, is often performed by inversion of P- and/or S-phases propagation data supplying the geometry of the medium (reflector depths) or any other structural parameter (P- or S-wave velocity, density...). Moreover, the inversion for velocity structure and interfaces is commonly performed using only seismic reflection travel times and/or crustal phase amplitudes in isotropic media. But it is very important to utilize more available information to constrain the non-uniqueness of the solution. In this paper, we present a simultaneous inversion method of seismic reflection travel times and polarizations data of transient elastic waves in stratified media to reconstruct not only layer depth and vertical P-wave velocity but also the anisotropy feature of the crust based on the estimation of the Thomsen’s parameters. We carry out a checking with synthetic data, comparing the inversion results obtained by anisotropic travel-time inversion to the results derived by joint inversion of seismic reflection travel times and polarizations data. The comparison proves that the first procedure leads to biased anisotropic models, while the second one fits nearly the real model. This makes the joint inversion method feasible. Finally, we investigate the geometry, P-wave velocity structure and anisotropy of the crust beneath Southeastern China by applying the proposed inversion method to previously acquired wide-angle seismic data. In this case, the anisotropy signature provides clear evidence that the Jiangshan-Shaoxing fault is the natural boundary between the Yangtze and Cathaysia blocks.     

A hybrid approach to compute seismic travel times in three-dimensional tetrahedral meshes

We propose an optimized method to compute travel times for seismic inversion problems. It is a hybrid method combining several approaches to deal with travel time computation accuracy in unstructured meshes based on tetrahedral elementary cells. As in the linear travel time interpolation method, the proposed approach computes travel times using seismic ray paths. The method operates in two sequential steps: At a first stage, travel times are computed for all nodes of the mesh using a modified version of the shortest path method. The difference with the standard version is that additional secondary nodes (called tertiary nodes) are added temporarily around seismic sources in order to improve accuracy with a reasonable increase in computational cost. During the second step, the steepest travel time gradient method is used to trace back ray paths for each source–receiver pair. Travel times at each receiver are then recomputed using slowness values at the intersection points between the ray path and the traversed cells. A number of numerical tests with an array of different velocity models, mesh resolutions and mesh topologies have been carried out. These tests showed that an average relative error in the order of 0.1% can be achieved at a computational cost that is suitable for travel time inversion.     

Non‐uniqueness with refraction inversion – the Mt Bulga shear zone

The tau‐p inversion algorithm is widely employed to generate starting models with many computer programs that implement refraction tomography. However, this algorithm can frequently fail to detect even major lateral variations in seismic velocities, such as a 50 m wide shear zone, which is the subject of this study. By contrast, the shear zone is successfully defined with the inversion algorithms of the generalized reciprocal method. The shear zone is confirmed with a 2D analysis of the head wave amplitudes, a spectral analysis of the refraction convolution section and with numerous closely spaced orthogonal seismic profiles recorded for a later 3D refraction investigation. Further improvements in resolution, which facilitate the recognition of additional zones with moderate reductions in seismic velocity, are achieved with a novel application of the Hilbert transform to the refractor velocity analysis algorithm. However, the improved resolution also requires the use of a lower average vertical seismic velocity, which accommodates a velocity reversal in the weathering. The lower seismic velocity is derived with the generalized reciprocal method, whereas most refraction tomography programs assume vertical velocity gradients as the default. Although all of the tomograms are consistent with the traveltime data, the resolution of each tomogram is comparable only with that of the starting model. Therefore, it is essential to employ inversion algorithms that can generate detailed starting models, where detailed lateral resolution is the objective. Non‐uniqueness can often be readily resolved with head wave amplitudes, attribute processing of the refraction convolution section and additional seismic traverses, prior to the acquisition of any borehole data. It is concluded that, unless specific measures are taken to address non‐uniqueness, the production of a single refraction tomogram that fits the traveltime data to sufficient accuracy does not necessarily demonstrate that the result is either correct, or even the most probable.     

Some remarks on Q-compensated sparse deconvolution without knowing the quality factor Q

The subsurface media are not perfectly elastic, thus anelastic absorption, attenuation and dispersion (aka Q filtering) effects occur during wave propagation, diminishing seismic resolution. Compensating for anelastic effects is imperative for resolution enhancement. Q values are required for most of conventional Q-compensation methods, and the source wavelet is additionally required for some of them. Based on the previous work of non-stationary sparse reflectivity inversion, we evaluate a series of methods for Q-compensation with/without knowing Q and with/without knowing wavelet. We demonstrate that if Q-compensation takes the wavelet into account, it generates better results for the severely attenuated components, benefiting from the sparsity promotion. We then evaluate a two-phase Q-compensation method in the frequency domain to eliminate Q requirement. In phase 1, the observed seismogram is disintegrated into the least number of Q-filtered wavelets chosen from a dictionary by optimizing a basis pursuit denoising problem, where the dictionary is composed of the known wavelet with different propagation times, each filtered with a range of possible values. The elements of the dictionary are weighted by the infinity norm of the corresponding column and further preconditioned to provide wavelets of different values and different propagation times equal probability to entry into the solution space. In phase 2, we derive analytic solutions for estimates of reflectivity and Q and solve an over-determined equation to obtain the final reflectivity series and Q values, where both the amplitude and phase information are utilized to estimate the Q values. The evaluated inversion-based Q estimation method handles the wave-interference effects better than conventional spectral-ratio-based methods. For Q-compensation, we investigate why sparsity promoting does matter. Numerical and field data experiments indicate the feasibility of the evaluated method of Q-compensation without knowing Q but with wavelet given.     

Estimating free parameters in 2D inversion: example of gravity inversion in a rifted basement

We advance a principle directed to assigning numerical values to free parameters usually present in inversion methods. It may be formulated as: ‘Optimum estimates of free parameters in an inversion procedure must lead, in tests using synthetic data, to solutions whose geometrical expression reflects the main qualitative or semiquantitative geological characteristic of the study area.’ To this end, the interpreter should (i) specify a typical anomalous source geometry which incorporates the most relevant geological information for the study area, (ii) compute the corresponding gravity anomaly and (iii) invert the anomaly for the source geometry finding the numerical values of the free parameters that lead to a solution closest to the typical source. Application of the above methodology to synthetic and real data from the basement relief of a rift basin has asserted its efficacy.     

A Bayesian model for lithology/fluid class prediction using a Markov mesh prior fitted from a training image

We consider a Bayesian model for inversion of observed amplitude variation with offset data into lithology/fluid classes, and study in particular how the choice of prior distribution for the lithology/fluid classes influences the inversion results. Two distinct prior distributions are considered, a simple manually specified Markov random field prior with a first-order neighbourhood and a Markov mesh model with a much larger neighbourhood estimated from a training image. They are chosen to model both horizontal connectivity and vertical thickness distribution of the lithology/fluid classes, and are compared on an offshore clastic oil reservoir in the North Sea. We combine both priors with the same linearized Gaussian likelihood function based on a convolved linearized Zoeppritz relation and estimate properties of the resulting two posterior distributions by simulating from these distributions with the Metropolis–Hastings algorithm. The influence of the prior on the marginal posterior probabilities for the lithology/fluid classes is clearly observable, but modest. The importance of the prior on the connectivity properties in the posterior realizations, however, is much stronger. The larger neighbourhood of the Markov mesh prior enables it to identify and model connectivity and curvature much better than what can be done by the first-order neighbourhood Markov random field prior. As a result, we conclude that the posterior realizations based on the Markov mesh prior appear with much higher lateral connectivity, which is geologically plausible.     

First‐arrival traveltime tomography with modified total‐variation regularization

First‐arrival traveltime tomography is a robust tool for near‐surface velocity estimation. A common approach to stabilizing the ill‐posed inverse problem is to apply Tikhonov regularization to the inversion. However, the Tikhonov regularization method recovers smooth local structures while blurring the sharp features in the model solution. We present a first‐arrival traveltime tomography method with modified total‐variation regularization to preserve sharp velocity contrasts and improve the accuracy of velocity inversion. To solve the minimization problem of the new traveltime tomography method, we decouple the original optimization problem into the two following subproblems: a standard traveltime tomography problem with the traditional Tikhonov regularization and a L₂ total‐variation problem. We apply the conjugate gradient method and split‐Bregman iterative method to solve these two subproblems, respectively. Our synthetic examples show that the new method produces higher resolution models than the conventional traveltime tomography with Tikhonov regularization, and creates less artefacts than the total variation regularization method for the models with sharp interfaces. For the field data, pre‐stack time migration sections show that the modified total‐variation traveltime tomography produces a near‐surface velocity model, which makes statics corrections more accurate.     

The utilization of the double focal transformation for sparse data representation and data reconstruction

In many cases, seismic measurements are coarsely sampled in at least one dimension. This leads to aliasing artefacts and therefore to problems in the subsequent processing steps. To avoid this, seismic data reconstruction can be applied in advance. The success and reliability of reconstruction methods are dependent on the assumptions they make on the data. In many cases, wavefields are assumed to (locally) have a linear space–time behaviour. However, field data are usually complex, with strongly curved events. Therefore, in this paper, we propose the double focal transformation as an efficient way for complex data reconstruction. Hereby, wavefield propagation is formulated as a transformation, where one‐way propagation operators are used as its basis functions. These wavefield operators can be based on a macro velocity model, which allows our method to use prior information in order to make the data decomposition more effective. The basic principle of the double focal transformation is to focus seismic energy along source and receiver coordinates simultaneously. The seismic data are represented by a number of localized events in the focal domain, whereas aliasing noise spreads out. By imposing a sparse solution in the focal domain, aliasing noise is suppressed, and data reconstruction beyond aliasing is achieved. To facilitate the process, only a few effective depth levels need to be included, preferably along the major boundaries in the data, from which the propagation operators can be calculated. Results on 2D and 3D synthetic data illustrate the method's virtues. Furthermore, seismic data reconstruction on a 2D field dataset with gaps and aliased source spacing demonstrates the strength of the double focal transformation, particularly for near‐offset reflections with strong curvature and for diffractions.     

三维复杂层状介质中多震相走时联合反演成像

采用新近提出的改进型不规则最短路径多次波射线追踪正演技术,结合共轭梯度法求解带约束的阻尼最小二乘最优化反演问题,讨论了三维复杂层状模型中利用多震相走时资料进行联合反演成像的技术方法.考虑到不同震相种类走时的拾取误差不同,反演算法中引入了不同震相种类数据的权系数;另外,考虑到同时反演速度模型和反射界面起伏中不同参数变化对走时影响程度的不同, Jacobi偏导矩阵元素中引入了不同种类参数的归一化因子.几种数值模拟实例表明: 多震相走时的联合或同时反演成像是一种提高走时成像空间分辨率,进而降低重建速度模型失真度行之有效的方法技术.     

Linear inversion via the discrete wavelet transform pseudoinverse

Classical least‐squares techniques (Moore–Penrose pseudoinverse) are covariance based and are therefore unsuitable for the solution of very large‐scale linear systems in geophysical inversion due to the need of diagonalisation. In this paper, we present a methodology to perform the geophysical inversion of large‐scale linear systems via the discrete wavelet transform. The methodology consists of compressing the linear system matrix using the interesting properties of covariance‐free orthogonal transformations, to design an approximation of the Moore–Penrose pseudoinverse. We show the application of the discrete wavelet transform pseudoinverse to well‐conditioned and ill‐conditioned linear systems. We applied the methodology to a general‐purpose linear problem where the system matrix has been generated using geostatistical simulation techniques and also to a synthetic 2D gravimetric problem with two different geological set‐ups, in the noise‐free and noisy cases. In both cases, the discrete wavelet transform pseudoinverse can be applied to the original linear system and also to the linear systems of normal equations and minimum norm. The results are compared with those obtained via the Moore–Penrose and the discrete cosine transform pseudoinverses. The discrete wavelet transform and the discrete cosine transform pseudoinverses provide similar results and outperform the Moore–Penrose pseudoinverse, mainly in the presence of noise. In the case of well‐conditioned linear systems, this methodology is more efficient when applied to the least‐squares system and minimum norm system due to their higher condition number that allows for a more efficient compression of the system matrix. Also, in the case of ill‐conditioned systems with very high underdetermined character, the application of the discrete cosine transform to the minimum norm solution provides very good results. Both solutions might differ on their regularity, depending on the wavelet family that is adopted. These methods have a general character and can be applied to solve any linear inverse problem arising in technology, particularly in geophysics, and also to non‐linear inversion by linearisation of the forward operator.     

Geomechanical property estimation of unconventional reservoirs using seismic data and rock physics

An extension of a previously developed rock physics model is made that quantifies the relationship between the ductile fraction of a brittle/ductile binary mixture and the isotropic seismic reflection response. By making a weak scattering (Born) approximation and plane wave (eikonal) approximation, with a subsequent ordering according to the angles of incidence, singular value decomposition analyses are performed to understand the stack weightings, number of stacks, and the type of stacks that will optimally estimate two fundamental rock physics parameters – the ductile fraction and the compaction and/or diagenesis. It is concluded that the full PP stack, i.e., sum of all PP offset traces, and the “full” PS stack, i.e., linear weighted sum of PS offset traces, are the two optimal stacks needed to estimate the two rock physics parameters. They dominate over both the second‐order amplitude variation offset “gradient” stack, which is a quadratically weighted sum of PP offset traces that is effectively the far offset traces minus the near offset traces, and the higher order fourth order PP stack (even at large angles of incidence). Using this result and model‐based Bayesian inversion, the seismic detectability of the ductile fraction (shown by others to be the important rock property for the geomechanical response of unconventional reservoir fracking) is demonstrated on a model characteristic of the Marcellus shale play.     

基于Born波路径的高斯束初至波波形反演

为了提高表层速度反演精度,本文提出了一种新的波形反演方法.该方法只利用初至波波形信息以减少波形反演对初始模型的依赖性,降低反演多解性与稳定性.由于只利用初至波波形信息,所以该方法利用高斯束计算格林函数和正演波场,以减少正演计算量.为了避免庞大核函数的存储,该方法基于Born波路径,利用矩阵分解算法实现方向与步长的累加计算.将此基于Born波路径的初至波波形反演方法应用于理论模型实验,并与声波方程全波形反演和初至波射线走时层析方法相对比,发现该方法的反演效果略低于全波形反演方法,但明显优于传统初至波射线走时层析方法,而计算效率却与射线走时层析相当.同时,相对于全波形反演,本文方法对初始模型的依赖性也有所降低.