Joint inversion of Rayleigh‐wave dispersion data and vertical electric sounding data: synthetic tests on characteristic sub‐surface models

In the traditional inversion of the Rayleigh dispersion curve, layer thickness, which is the second most sensitive parameter of modelling the Rayleigh dispersion curve, is usually assumed as correct and is used as fixed a priori information. Because the knowledge of the layer thickness is typically not precise, the use of such a priori information may result in the traditional Rayleigh dispersion curve inversions getting trapped in some local minima and may show results that are far from the real solution. In this study, we try to avoid this issue by using a joint inversion of the Rayleigh dispersion curve data with vertical electric sounding data, where we use the common‐layer thickness to couple the two methods. The key idea of the proposed joint inversion scheme is to combine methods in one joint Jacobian matrix and to invert for layer S‐wave velocity, resistivity, and layer thickness as an additional parameter, in contrast with a traditional Rayleigh dispersion curve inversion. The proposed joint inversion approach is tested with noise‐free and Gaussian noise data on six characteristic, synthetic sub‐surface models: a model with a typical dispersion; a low‐velocity, half‐space model; a model with particularly stiff and soft layers, respectively; and a model reproduced from the stiff and soft layers for different layer‐resistivity propagation. In the joint inversion process, the non‐linear damped least squares method is used together with the singular value decomposition approach to find a proper damping value for each iteration. The proposed joint inversion scheme tests many damping values, and it chooses the one that best approximates the observed data in the current iteration. The quality of the joint inversion is checked with the relative distance measure. In addition, a sensitivity analysis is performed for the typical dispersive sub‐surface model to illustrate the benefits of the proposed joint scheme. The results of synthetic models revealed that the combination of the Rayleigh dispersion curve and vertical electric sounding methods in a joint scheme allows to provide reliable sub‐surface models even in complex and challenging situations and without using any a priori information.     

Iteratively re‐weighted and refined least squares algorithm for robust inversion of geophysical data

A robust metric of data misfit such as the ?₁‐norm is required for geophysical parameter estimation when the data are contaminated by erratic noise. Recently, the iteratively re‐weighted and refined least‐squares algorithm was introduced for efficient solution of geophysical inverse problems in the presence of additive Gaussian noise in the data. We extend the algorithm in two practically important directions to make it applicable to data with non‐Gaussian noise and to make its regularisation parameter tuning more efficient and automatic. The regularisation parameter in iteratively reweighted and refined least‐squares algorithm varies with iteration, allowing the efficient solution of constrained problems. A technique is proposed based on the secant method for root finding to concentrate on finding a solution that satisfies the constraint, either fitting to a target misfit (if a bound on the noise is available) or having a target size (if a bound on the solution is available). This technique leads to an automatic update of the regularisation parameter at each and every iteration. We further propose a simple and efficient scheme that tunes the regularisation parameter without requiring target bounds. This is of great importance for the field data inversion where there is no information about the size of the noise and the solution. Numerical examples from non‐stationary seismic deconvolution and velocity‐stack inversion show that the proposed algorithm is efficient, stable, and robust and outperforms the conventional and state‐of‐the‐art methods.     

Misfit functionals in Laplace‐Fourier domain waveform inversion,with application to wide‐angle ocean bottom seismograph data

In seismic waveform inversion, non‐linearity and non‐uniqueness require appropriate strategies. We formulate four types of L2 normed misfit functionals for Laplace‐Fourier domain waveform inversion: i) subtraction of complex‐valued observed data from complex‐valued predicted data (the ‘conventional phase‐amplitude’ residual), ii) a ‘conventional phase‐only’ residual in which amplitude variations are normalized, iii) a ‘logarithmic phase‐amplitude’ residual and finally iv) a ‘logarithmic phase‐only’ residual in which the only imaginary part of the logarithmic residual is used. We evaluate these misfit functionals by using a wide‐angle field Ocean Bottom Seismograph (OBS) data set with a maximum offset of 55 km. The conventional phase‐amplitude approach is restricted in illumination and delineates only shallow velocity structures. In contrast, the other three misfit functionals retrieve detailed velocity structures with clear lithological boundaries down to the deeper part of the model. We also test the performance of additional phase‐amplitude inversions starting from the logarithmic phase‐only inversion result. The resulting velocity updates are prominent only in the high‐wavenumber components, sharpening the lithological boundaries. We argue that the discrepancies in the behaviours of the misfit functionals are primarily caused by the sensitivities of the model gradient to strong amplitude variations in the data. As the observed data amplitudes are dominated by the near‐offset traces, the conventional phase‐amplitude inversion primarily updates the shallow structures as a result. In contrast, the other three misfit functionals eliminate the strong dependence on amplitude variation naturally and enhance the depth of illumination. We further suggest that the phase‐only inversions are sufficient to obtain robust and reliable velocity structures and the amplitude information is of secondary importance in constraining subsurface velocity models.     

P‐wave stacking‐velocity tomography for VTI media

A major complication caused by anisotropy in velocity analysis and imaging is the uncertainty in estimating the vertical velocity and depth scale of the model from surface data. For laterally homogeneous VTI (transversely isotropic with a vertical symmetry axis) media above the target reflector, P‐wave moveout has to be combined with other information (e.g. borehole data or converted waves) to build velocity models for depth imaging. The presence of lateral heterogeneity in the overburden creates the dependence of P‐wave reflection data on all three relevant parameters (the vertical velocity V_P0 and the Thomsen coefficients ε and δ) and, therefore, may help to determine the depth scale of the velocity field. Here, we propose a tomographic algorithm designed to invert NMO ellipses (obtained from azimuthally varying stacking velocities) and zero‐offset traveltimes of P‐waves for the parameters of homogeneous VTI layers separated by either plane dipping or curved interfaces. For plane non‐intersecting layer boundaries, the interval parameters cannot be recovered from P‐wave moveout in a unique way. Nonetheless, if the reflectors have sufficiently different azimuths, a priori knowledge of any single interval parameter makes it possible to reconstruct the whole model in depth. For example, the parameter estimation becomes unique if the subsurface layer is known to be isotropic. In the case of 2D inversion on the dip line of co‐orientated reflectors, it is necessary to specify one parameter (e.g. the vertical velocity) per layer. Despite the higher complexity of models with curved interfaces, the increased angle coverage of reflected rays helps to resolve the trade‐offs between the medium parameters. Singular value decomposition (SVD) shows that in the presence of sufficient interface curvature all parameters needed for anisotropic depth processing can be obtained solely from conventional‐spread P‐wave moveout. By performing tests on noise‐contaminated data we demonstrate that the tomographic inversion procedure reconstructs both the interfaces and the VTI parameters with high accuracy. Both SVD analysis and moveout inversion are implemented using an efficient modelling technique based on the theory of NMO‐velocity surfaces generalized for wave propagation through curved interfaces.     

Assessment of the reliability of 2D inversion of apparent resistivity data

The reliability of inversion of apparent resistivity pseudosection data to determine accurately the true resistivity distribution over 2D structures has been investigated, using a common inversion scheme based on a smoothness‐constrained non‐linear least‐squares optimization, for the Wenner array. This involved calculation of synthetic apparent resistivity pseudosection data, which were then inverted and the model estimated from the inversion was compared with the original 2D model. The models examined include (i) horizontal layering, (ii) a vertical fault, (iii) a low‐resistivity fill within a high‐resistivity basement, and (iv) an upfaulted basement block beneath a conductive overburden. Over vertical structures, the resistivity models obtained from inversion are usually much sharper than the measured data. However, the inverted resistivities can be smaller than the lowest, or greater than the highest, true model resistivity. The substantial reduction generally recorded in the data misfit during the least‐squares inversion of 2D apparent resistivity data is not always accompanied by any noticeable reduction in the model misfit. Conversely, the model misfit may, for all practical purposes, remain invariant for successive iterations. It can also increase with the iteration number, especially where the resistivity contrast at the bedrock interface exceeds a factor of about 10; in such instances, the optimum model estimated from inversion is attained at a very low iteration number. The largest model misfit is encountered in the zone adjacent to a contact where there is a large change in the resistivity contrast. It is concluded that smooth inversion can provide only an approximate guide to the true geometry and true formation resistivity.     

The use of low frequencies in a full‐waveform inversion and impedance inversion land seismic case study

Velocity model building and impedance inversion generally suffer from a lack of intermediate wavenumber content in seismic data. Intermediate wavenumbers may be retrieved directly from seismic data sets if enough low frequencies are recorded. Over the past years, improvements in acquisition have allowed us to obtain seismic data with a broader frequency spectrum. To illustrate the benefits of broadband acquisition, notably the recording of low frequencies, we discuss the inversion of land seismic data acquired in Inner Mongolia, China. This data set contains frequencies from 1.5–80 Hz. We show that the velocity estimate based on an acoustic full‐waveform inversion approach is superior to one obtained from reflection traveltime inversion because after full‐waveform inversion the background velocity conforms to geology. We also illustrate the added value of low frequencies in an impedance estimate.     

Linearized seismic waveform inversion using a multiple re-weighted least-squares method with QR preconditioning

Linearized inversion methods such as Gauss‐Newton and multiple re‐weighted least‐squares are iterative processes in which an update in the current model is computed as a function of data misfit and the gradient of data with respect to model parameters. The main advantage of those methods is their ability to refine the model parameters although they have a high computational cost for seismic inversion. In the Gauss‐Newton method a system of equations, corresponding to the sensitivity matrix, is solved in the least‐squares sense at each iteration, while in the multiple re‐weighted least‐squares method many systems are solved using the same sensitivity matrix. The sensitivity matrix arising from these methods is usually not sparse, thus limiting the use of standard preconditioners in the solution of the linearized systems. For reduction of the computational cost of the linearized inversion methods, we propose the use of preconditioners based on a partial orthogonalization of the columns of the sensitivity matrix. The new approach collapses a band of co‐diagonals of the normal equations matrix into the main diagonal, being equivalent to computing the least‐squares solution starting from a partial solution of the linear system. The preconditioning is driven by a bandwidth L which can be interpreted as the distance for which the correlation between model parameters is relevant. To illustrate the benefit of the proposed approach to the reduction of the computational cost of the inversion we apply the multiple re‐weighted least‐squares method to the 2D acoustic seismic waveform inversion problem. We verify the reduction in the number of iterations in the conjugate'gradient algorithm as the bandwidth of the preconditioners increases. This effect reduces the total computational cost of inversion as well.     

Enforcing smoothness and assessing uncertainty in non-linear one-dimensional prestack seismic inversion

Estimation of elastic properties of rock formations from surface seismic amplitude measurements remains a subject of interest for the exploration and development of hydrocarbon reservoirs. This paper develops a global inversion technique to estimate and appraise 1D distributions of compressional‐wave velocity, shear‐wave velocity and bulk density, from normal‐moveout‐corrected PP prestack surface seismic amplitude measurements. Specific objectives are: (a) to evaluate the efficiency of the minimization algorithm (b) to appraise the impact of various data misfit functions, and (c) to assess the effect of the degree and type of smoothness criterion enforced by the inversion. Numerical experiments show that very fast simulated annealing is the most efficient minimization technique among alternative approaches considered for global inversion. It is also found that an adequate choice of data misfit function is necessary for a reliable and efficient match of noisy and sparse seismic amplitude measurements. Several procedures are considered to enforce smoothness of the estimated 1D distributions of elastic parameters, including predefined quadratic measures of length, flatness and roughness. Based on the general analysis of global inversion techniques, we introduce a new stochastic inversion algorithm that initializes the search for the minimum with constrained random distributions of elastic parameters and enforces predefined autocorrelation functions (semivariograms). This strategy readily lends itself to the assessment of model uncertainty. The new global inversion algorithm is successfully tested on noisy synthetic amplitude data. Moreover, we present a feasibility analysis of the resolution and uncertainty of prestack seismic amplitude data to infer 1D distributions of elastic parameters measured with wireline logs in the deepwater Gulf of Mexico. The new global inversion algorithm is computationally more efficient than the alternative global inversion procedures considered here.     

Least‐squares reverse‐time migration in a matrix‐based formulation

This paper describes least‐squares reverse‐time migration. The method provides the exact adjoint operator pair for solving the linear inverse problem, thereby enhancing the convergence of gradient‐based iterative linear inversion methods. In this formulation, modified source wavelets are used to correct the source signature imprint in the predicted data. Moreover, a roughness constraint is applied to stabilise the inversion and reduce high‐wavenumber artefacts. It is also shown that least‐squares migration implicitly applies a deconvolution imaging condition. Three numerical experiments illustrate that this method is able to produce seismic reflectivity images with higher resolution, more accurate amplitudes, and fewer artefacts than conventional reverse‐time migration. The methodology is currently feasible in 2‐D and can naturally be extended to 3‐D when computational resources become more powerful.     

Correlation‐based reflection waveform inversion by one‐way wave equations

Reflection full waveform inversion can update subsurface velocity structure of the deeper part, but tends to get stuck in the local minima associated with the waveform misfit function. These local minima cause cycle skipping if the initial background velocity model is far from the true model. Since conventional reflection full waveform inversion using two‐way wave equation in time domain is computationally expensive and consumes a large amount of memory, we implement a correlation‐based reflection waveform inversion using one‐way wave equations to retrieve the background velocity. In this method, one‐way wave equations are used for the seismic wave forward modelling, migration/de‐migration and the gradient computation of objective function in frequency domain. Compared with the method using two‐way wave equation, the proposed method benefits from the lower computational cost of one‐way wave equations without significant accuracy reduction in the cases without steep dips. It also largely reduces the memory requirement by an order of magnitude than implementation using two‐way wave equation both for two‐ and three‐dimensional situations. Through numerical analysis, we also find that one‐way wave equations can better construct the low wavenumber reflection wavepath without producing high‐amplitude short‐wavelength components near the image points in the reflection full waveform inversion gradient. Synthetic test and real data application show that the proposed method efficiently updates the background velocity model.     

频率域粘弹性声波透射波形速度反演

在用稀疏矩阵的LU分解技术对频率域粘弹性声波方程进行直接求解的基础上,根据失配函数二范数最小准则, 用预条件梯度类方法对粘弹性声波介质的速度结构进行了逐频反演. 局部非均匀介质模型和层状介质模型速度结构反演的实验结果表明,不同频率能够反映地下介质的多尺度物性结构（低频数据对应与介质物性的大尺度结构）,用低频反演结果作为高频反演的初值逼近这一顺序模式,能大大改善反演过程中解的非唯一性. 而且,在反演过程中用Hess矩阵的对角线元素来做梯度类方法的预条件算子,能够吸收了高斯牛顿法的二次收敛优势, 使得本文算法具有较快的收敛速度.      

Visco-acoustic transmission waveform inversion of velocity structure in space-frequency domain

According to the least square criterion of minimizing the misfit between modeled and observed data,this paper provides a preconditioned gradient method to invert the visco-acoustic velocity structure on the basis of using sparse matrix LU factorization technique to directly solve the visco-acoustic wave forward problem in space-frequency domain.Numerical results obtained in an inclusion model inversion and a layered homogeneous model inversion demonstrate that different scale media have their own frequency ...     

最小二乘反演迭代算法在压制海上变深度缆采集数据虚反射中的应用

海水面的虚反射(鬼波)引起海上拖缆采集数据陷波,导致地震记录频带变窄,而近年发展的变深度缆采集技术,具有多样的陷波特征,通过专门的去虚反射处理方法可获得宽频数据.本文基于已有研究成果,将最小二乘反演迭代压制虚反射算法应用于某海上变深度缆宽频处理.基于频率波数域镜像记录生成方法获得镜像炮集记录,并采用最小二乘解从变深度缆原始和镜像炮集记录中提取上行波.针对镜像炮集记录生成受初始速度模型精度的影响,使得某深度缆接收的上行波和下行波之间的实际延迟时间存在误差,采用最小二乘反演迭代算法最优化计算下行波与上行波之间的平均延迟时间和上行波记录,并采用时空数据窗口滑动克服延迟时间随炮检距和目的层深度变化问题.合成数据及某海上实际变深度缆数据处理测试结果表明,该方法能较好地压制变深度缆由海水面产生的虚反射,能达到拓宽地震记录频带目的.     

Wave equation least square imaging using the local angular Hessian for amplitude correction

Local angular Hessian can be used to improve wave equation least square migration images. By decomposing the original Hessian operator into the local wavenumber domain or the local angle domain, the least square migration image is obtained as the solution of a linearized least‐squares inversion in the frequency and local angle domains. The local angular Hessian contains information about the acquisition geometry and the propagation effects based on the given velocity model. The inversion scheme based on the local angular Hessian avoids huge computation on the exact inverse Hessian matrix. To reduce the instability in the inversion, damping factors are introduced into the deconvolution filter in the local wavenumber domain and the local angle domain. The algorithms are tested using the SEG/EAGE salt2D model and the Sigsbee2A model. Results show improved image quality and amplitudes.     

Stability of source mechanisms inverted from P‐wave amplitude microseismic monitoring data acquired at the surface

We study the stability of source mechanisms inverted from data acquired at surface and near‐surface monitoring arrays. The study is focused on P‐wave data acquired on vertical components, as this is the most common type of acquisition. We apply ray modelling on three models: a fully homogeneous isotropic model, a laterally homogeneous isotropic model and a laterally homogeneous anisotropic model to simulate three commonly used models in inversion. We use geometries of real arrays, one consisting in surface receivers and one consisting in ‘buried’ geophones at the near‐surface. Stability was tested for two of the frequently observed source mechanisms: strike‐slip and dip‐slip and was evaluated by comparing the parameters of correct and inverted mechanisms. We assume these double‐couple source mechanisms and use quantitatively the inversion allowing non‐double‐couple components to measure stability of the inversion. To test the robustness we inverted synthetic amplitudes computed for a laterally homogeneous isotropic model and contaminated with noise using a fully homogeneous model in the inversion. Analogously amplitudes computed in a laterally homogeneous anisotropic model were inverted in all three models. We show that a star‐like surface acquisition array provides very stable inversion up to a very high level of noise in data. Furthermore, we reveal that strike‐slip inversion is more stable than dip‐slip inversion for the receiver geometries considered here. We show that noise and an incorrect velocity model may result in narrow bands of source mechanisms in Hudson's plots.     

Iterative resolution estimation in least‐squares Kirchhoff migration

We apply iterative resolution estimation to least‐squares Kirchhoff migration. Reviewing the theory of iterative optimization uncovers the common origin of different optimization methods. This allows us to reformulate the pseudo‐inverse, model resolution and data resolution operators in terms of effective iterative estimates. When applied to Kirchhoff migration, plots of the diagonal of the model resolution matrix reveal low illumination areas on seismic images and provide information about image uncertainties. Synthetic and real data examples illustrate the proposed technique and confirm the theoretical expectations.     

Improved Monte Carlo inversion of surface wave data

Inversion of surface wave data suffers from solution non‐uniqueness and is hence strongly biased by the initial model. The Monte Carlo approach can handle this non‐uniqueness by evidencing the local minima but it is inefficient for high dimensionality problems and makes use of subjective criteria, such as misfit thresholds, to interpret the results. If a smart sampling of the model parameter space, which exploits scale properties of the modal curves, is introduced the method becomes more efficient and with respect to traditional global search methods it avoids the subjective use of control parameters that are barely related to the physical problem. The results are interpreted drawing inference by means of a statistical test that selects an ensemble of feasible shear wave velocity models according to data quality and model parameterization. Tests on synthetic data demonstrate that the application of scale properties concentrates the sampling of model parameter space in high probability density zones and makes it poorly sensitive to the initial boundary of the model parameters. Tests on synthetic and field data, where boreholes are available, prove that the statistical test selects final results that are consistent with the true model and which are sensitive to data quality. The implemented strategies make the Monte Carlo inversion efficient for practical applications and able to effectively retrieve subsoil models even in complex and challenging situations such as velocity inversions.     

Two-dimensional inversion of magnetotelluric data with consecutive use of conjugate gradient and least-squares solution with singular value decomposition algorithms

I investigated the two‐dimensional magnetotelluric data inversion algorithms in studying two significant aspects within a linearized inversion approach. The first one is the method of minimization and second one is the type of stabilizing functional used in parametric functionals. The results of two well‐known inversion algorithms, namely conjugate gradient and the least‐squares solution with singular value decomposition, were compared in terms of accuracy and CPU time. In addition, magnetotelluric data inversion with various stabilizers, such as L2‐norm, smoothing, minimum support, minimum gradient support and first‐order minimum entropy, were examined. A new inversion algorithm named least‐squares solution with singular value decomposition and conjugate gradient is suggested in seeing the outcomes of the comparisons carried out on least‐squares solutions with singular value decomposition and conjugate gradient algorithms subject to a variety of stabilizers. Inversion results of synthetic data showed that the newly suggested algorithm yields better results than those of the individual implementations of conjugate gradient and least‐squares solution with singular value decomposition algorithms. The suggested algorithm and the above‐mentioned algorithms inversion results for the field data collected along a line crossing the North Anatolian Fault zone were also compared each other and results are discussed.     

Non‐uniqueness with refraction inversion – a syncline model study

Non‐uniqueness occurs with the 1D parametrization of refraction traveltime graphs in the vertical dimension and with the 2D lateral resolution of individual layers in the horizontal dimension. The most common source of non‐uniqueness is the inversion algorithm used to generate the starting model. This study applies 1D, 1.5D and 2D inversion algorithms to traveltime data for a syncline (2D) model, in order to generate starting models for wave path eikonal traveltime tomography. The 1D tau‐p algorithm produced a tomogram with an anticline rather than a syncline and an artefact with a high seismic velocity. The 2D generalized reciprocal method generated tomograms that accurately reproduced the syncline, together with narrow regions at the thalweg with seismic velocities that are less than and greater than the true seismic velocities as well as the true values. It is concluded that 2D inversion algorithms, which explicitly identify forward and reverse traveltime data, are required to generate useful starting models in the near‐surface where irregular refractors are common. The most likely tomogram can be selected as either the simplest model or with a priori information, such as head wave amplitudes. The determination of vertical velocity functions within individual layers is also subject to non‐uniqueness. Depths computed with vertical velocity gradients, which are the default with many tomography programs, are generally 50% greater than those computed with constant velocities for the same traveltime data. The average vertical velocity provides a more accurate measure of depth estimates, where it can be derived. Non‐uniqueness is a fundamental reality with the inversion of all near‐surface seismic refraction data. Unless specific measures are taken to explicitly address non‐uniqueness, then the production of a single refraction tomogram, which fits the traveltime data to sufficient accuracy, does not necessarily demonstrate that the result is either ‘correct’ or the most probable.     

A robust approach to time‐to‐depth conversion and interval velocity estimation from time migration in the presence of lateral velocity variations

The problem of conversion from time‐migration velocity to an interval velocity in depth in the presence of lateral velocity variations can be reduced to solving a system of partial differential equations. In this paper, we formulate the problem as a non‐linear least‐squares optimization for seismic interval velocity and seek its solution iteratively. The input for the inversion is the Dix velocity, which also serves as an initial guess. The inversion gradually updates the interval velocity in order to account for lateral velocity variations that are neglected in the Dix inversion. The algorithm has a moderate cost thanks to regularization that speeds up convergence while ensuring a smooth output. The proposed method should be numerically robust compared to the previous approaches, which amount to extrapolation in depth monotonically. For a successful time‐to‐depth conversion, image‐ray caustics should be either nonexistent or excluded from the computational domain. The resulting velocity can be used in subsequent depth‐imaging model building. Both synthetic and field data examples demonstrate the applicability of the proposed approach.