Non-monotone spectral projected gradient method applied to full waveform inversion

The seismic inversion problem is a highly non‐linear problem that can be reduced to the minimization of the least‐squares criterion between the observed and the modelled data. It has been solved using different classical optimization strategies that require a monotone descent of the objective function. We propose solving the full‐waveform inversion problem using the non‐monotone spectral projected gradient method: a low‐cost and low‐storage optimization technique that maintains the velocity values in a feasible convex region by frequently projecting them on this convex set. The new methodology uses the gradient direction with a particular spectral step length that allows the objective function to increase at some iterations, guarantees convergence to a stationary point starting from any initial iterate, and greatly speeds up the convergence of gradient methods. We combine the new optimization scheme as a solver of the full‐waveform inversion with a multiscale approach and apply it to a modified version of the Marmousi data set. The results of this application show that the proposed method performs better than the classical gradient method by reducing the number of function evaluations and the residual values.     

Source-independent waveform inversion for attenuation estimation in anisotropic media

In previous publications, we presented a waveform-inversion algorithm for attenuation analysis in heterogeneous anisotropic media. However, waveform inversion requires an accurate estimate of the source wavelet, which is often difficult to obtain from field data. To address this problem, here we adopt a source-independent waveform-inversion algorithm that obviates the need for joint estimation of the source signal and attenuation coefficients. The key operations in that algorithm are the convolutions (1) of the observed wavefield with a reference trace from the modelled data and (2) of the modelled wavefield with a reference trace from the observed data. The influence of the source signature on attenuation estimation is mitigated by defining the objective function as the ℓ₂-norm of the difference between the two convolved data sets. The inversion gradients for the medium parameters are similar to those for conventional waveform-inversion techniques, with the exception of the adjoint sources computed by convolution and cross-correlation operations. To make the source-independent inversion methodology more stable in the presence of velocity errors, we combine it with the local-similarity technique. The proposed algorithm is validated using transmission tests for a homogeneous transversely isotropic model with a vertical symmetry axis that contains a Gaussian anomaly in the shear-wave vertical attenuation coefficient. Then the method is applied to the inversion of reflection data for a modified transversely isotropic model from Hess. It should be noted that due to the increased nonlinearity of the inverse problem, the source-independent algorithm requires a more accurate initial model to obtain inversion results comparable to those produced by conventional waveform inversion with the actual wavelet.     

Offset‐variable density improves acoustic full‐waveform inversion: a shallow marine case study

We have previously applied three‐dimensional acoustic, anisotropic, full‐waveform inversion to a shallow‐water, wide‐angle, ocean‐bottom‐cable dataset to obtain a high‐resolution velocity model. This velocity model produced an improved match between synthetic and field data, better flattening of common‐image gathers, a closer fit to well logs, and an improvement in the pre‐stack depth‐migrated image. Nevertheless, close examination reveals that there is a systematic mismatch between the observed and predicted data from this full‐waveform inversion model, with the predicted data being consistently delayed in time. We demonstrate that this mismatch cannot be produced by systematic errors in the starting model, by errors in the assumed source wavelet, by incomplete convergence, or by the use of an insufficiently fine finite‐difference mesh. Throughout these tests, the mismatch is remarkably robust with the significant exception that we do not see an analogous mismatch when inverting synthetic acoustic data. We suspect therefore that the mismatch arises because of inadequacies in the physics that are used during inversion. For ocean‐bottom‐cable data in shallow water at low frequency, apparent observed arrival times, in wide‐angle turning‐ray data, result from the characteristics of the detailed interference pattern between primary refractions, surface ghosts, and a large suite of wide‐angle multiple reflected and/or multiple refracted arrivals. In these circumstances, the dynamics of individual arrivals can strongly influence the apparent arrival times of the resultant compound waveforms. In acoustic full‐waveform inversion, we do not normally know the density of the seabed, and we do not properly account for finite shear velocity, finite attenuation, and fine‐scale anisotropy variation, all of which can influence the relative amplitudes of different interfering arrivals, which in their turn influence the apparent kinematics. Here, we demonstrate that the introduction of a non‐physical offset‐variable water density during acoustic full‐waveform inversion of this ocean‐bottom‐cable field dataset can compensate efficiently and heuristically for these inaccuracies. This approach improves the travel‐time match and consequently increases both the accuracy and resolution of the final velocity model that is obtained using purely acoustic full‐waveform inversion at minimal additional cost.     

Comparison of waveform inversion, part 1: conventional wavefield vs logarithmic wavefield

This is the first in a series of three papers focused on using variants of a logarithmic objective function approach to full waveform inversion. In this article, we investigate waveform inversion using full logarithmic principles and compare the results with the conventional least squares approach. We demonstrate theoretically that logarithmic inversion is computational similar to the conventional method in the sense that it uses exactly the same back‐propagation technology as used in least‐squares inversion. In the sense that it produces better results for each of three numerical examples, we conclude that logarithmic inversion is also more robust. We argue that a major reason for the inherent robustness is the fact that the logarithmic approach produces a natural scaling of the amplitude of the residual wavefield by the amplitude of the modelled wavefield that tends to stabilize the computations and consequently improve the final result. We claim that any superiority of the logarithmic inversion is based on the fact that it tends to be tomographic in the early stage of the inversion and more dependent on amplitude differences in the latter stages.     

不依赖源子波的跨孔雷达时间域波形反演

波形反演是近年来较热门的反演方法,其分辨率可以达到亚波长级别.在波形反演的实际应用中,源子波的估计十分重要.传统方法使用反褶积来估计源子波并随着反演过程更新,该方法在合成数据波形反演中效果较好,但在实际数据反演过程中存在一系列的问题.由于实际数据信噪比较低,在源子波估计过程中需要大量的人为干涉,且结果并不一定可靠.本文使用一种基于褶积波场的新型目标函数,令反演过程不再依赖源子波.详细推导了针对跨孔雷达波形反演的梯度及步长公式,实现介电常数和电导率的同步反演.针对一个合成数据模型同时反演介电常数和电导率,结果表明该方法能够反演出亚波长尺寸异常体的形状和位置.接着,将该方法应用到两组实际数据中,并与基于估计源子波的时间域波形反演结果进行比较.结果表明不依赖源子波的时间域波形反演结果分辨率更高,也更准确.     

基于弹性波全波形反演的主被动源多分量混采地震数据速度建模

在常规地震采集中,被动源地震波场往往被视为噪声而去除,这就造成了部分有用信息的丢失.在目标区进行主动源和被动源弹性波地震数据的多分量混合采集,并对两种数据进行联合应用,使其在照明和频带上优势互补,能显著提高成像和反演的质量.本文针对两种不同类型的主被动源混采地震数据,分别提出了相应的联合全波形反演方法.首先,针对主动源与瞬态被动源弹性波混采地震数据,为充分利用被动源对深部照明的优势,同时有效压制被动震源点附近的成像异常值,提出了基于动态随机组合的弹性波被动源照明补偿反演策略.然后,针对低频缺失主动源与背景噪声型被动源弹性波混采地震数据,为充分利用被动源波场携带的低频信息,并避免对被动源的定位和子波估计,提出了基于地震干涉与不依赖子波算法的弹性波主被动源串联反演策略.最后,分别将两种方法在Marmousi模型上进行反演测试.结果说明,综合利用主动源和被动源弹性波混采地震数据,不仅能增强深部弹性参数反演效果,还能更好地构建弹性参数模型的宏观结构,并有助于缓解常规弹性波全波形反演的跳周问题.     

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.     

Migration velocity analysis and waveform inversion

Least‐squares inversion of seismic reflection waveform data can reconstruct remarkably detailed models of subsurface structure and take into account essentially any physics of seismic wave propagation that can be modelled. However, the waveform inversion objective has many spurious local minima, hence convergence of descent methods (mandatory because of problem size) to useful Earth models requires accurate initial estimates of long‐scale velocity structure. Migration velocity analysis, on the other hand, is capable of correcting substantially erroneous initial estimates of velocity at long scales. Migration velocity analysis is based on prestack depth migration, which is in turn based on linearized acoustic modelling (Born or single‐scattering approximation). Two major variants of prestack depth migration, using binning of surface data and Claerbout's survey‐sinking concept respectively, are in widespread use. Each type of prestack migration produces an image volume depending on redundant parameters and supplies a condition on the image volume, which expresses consistency between data and velocity model and is hence a basis for velocity analysis. The survey‐sinking (depth‐oriented) approach to prestack migration is less subject to kinematic artefacts than is the binning‐based (surface‐oriented) approach. Because kinematic artefacts strongly violate the consistency or semblance conditions, this observation suggests that velocity analysis based on depth‐oriented prestack migration may be more appropriate in kinematically complex areas. Appropriate choice of objective (differential semblance) turns either form of migration velocity analysis into an optimization problem, for which Newton‐like methods exhibit little tendency to stagnate at nonglobal minima. The extended modelling concept links migration velocity analysis to the apparently unrelated waveform inversion approach to estimation of Earth structure: from this point of view, migration velocity analysis is a solution method for the linearized waveform inversion problem. Extended modelling also provides a basis for a nonlinear generalization of migration velocity analysis. Preliminary numerical evidence suggests a new approach to nonlinear waveform inversion, which may combine the global convergence of velocity analysis with the physical fidelity of model‐based data fitting.     

Full waveform inversion of SH‐ and Love‐wave data in near‐surface prospecting

We develop a two‐dimensional full waveform inversion approach for the simultaneous determination of S‐wave velocity and density models from SH ‐ and Love‐wave data. We illustrate the advantages of the SH/Love full waveform inversion with a simple synthetic example and demonstrate the method's applicability to a near‐surface dataset, recorded in the village ?achtice in Northwestern Slovakia. Goal of the survey was to map remains of historical building foundations in a highly heterogeneous subsurface. The seismic survey comprises two parallel SH‐profiles with maximum offsets of 24 m and covers a frequency range from 5 Hz to 80 Hz with high signal‐to‐noise ratio well suited for full waveform inversion. Using the Wiechert–Herglotz method, we determined a one‐dimensional gradient velocity model as a starting model for full waveform inversion. The two‐dimensional waveform inversion approach uses the global correlation norm as objective function in combination with a sequential inversion of low‐pass filtered field data. This mitigates the non‐linearity of the multi‐parameter inverse problem. Test computations show that the influence of visco‐elastic effects on the waveform inversion result is rather small. Further tests using a mono‐parameter shear modulus inversion reveal that the inversion of the density model has no significant impact on the final data fit. The final full waveform inversion S‐wave velocity and density models show a prominent low‐velocity weathering layer. Below this layer, the subsurface is highly heterogeneous. Minimum anomaly sizes correspond to approximately half of the dominant Love‐wavelength. The results demonstrate the ability of two‐dimensional SH waveform inversion to image shallow small‐scale soil structure. However, they do not show any evidence of foundation walls.     

基于归一化能量谱目标函数的全波形反演方法

基于L2范数的常规全波形反演目标函数是一个强非线性泛函,在反演过程中容易陷入局部极小值.本文提出归一化能量谱目标函数来缓解全波形反演过程中的强非线性问题,同时能够有效地缓解噪声和震源子波不准等因素的影响.能量谱目标函数是通过匹配观测数据与模拟数据随频率分布的能量信息来实现最小二乘反演的,其忽略了地震数据波形与相位变化的细节特征,这在反演的过程中能够有效缓解波形匹配错位等问题.数值测试结果表明,基于归一化能量谱目标函数在构建初始速度模型、抗噪性和缓解震源子波依赖等方面都优于归一化全波形反演目标函数.金属矿模型测试结果表明,即使地震数据缺失低频分量,基于归一化能量谱目标函数的全波形反演方法在像金属矿这样的强散射介质反演问题上同样具有一定的优势.     

Exploring some issues in acoustic full waveform inversion

The least‐squares error measures the difference between observed and modelled seismic data. Because it suffers from local minima, a good initial velocity model is required to avoid convergence to the wrong model when using a gradient‐based minimization method. If a data set mainly contains reflection events, it is difficult to update the velocity model with the least‐squares error because the minimization method easily ends up in the nearest local minimum without ever reaching the global minimum. Several authors observed that the model could be updated by diving waves, requiring a wide‐angle or large‐offset data set. This full waveform tomography is limited to a maximum depth. Here, we use a linear velocity model to obtain estimates for the maximum depth. In addition, we investigate how frequencies should be selected if the seismic data are modelled in the frequency domain. In the presence of noise, the condition to avoid local minima requires more frequencies than needed for sufficient spectral coverage. We also considered acoustic inversion of a synthetic marine data set created by an elastic time‐domain finite‐difference code. This allowed us to validate the estimates made for the linear velocity model. The acoustic approximation leads to a number of problems when using long‐offset data. Nevertheless, we obtained reasonable results. The use of a variable density in the acoustic inversion helped to match the data at the expense of accuracy in the inversion result for the density.     

Velocity model building from seismic reflection data by full‐waveform inversion

Full‐waveform inversion is re‐emerging as a powerful data‐fitting procedure for quantitative seismic imaging of the subsurface from wide‐azimuth seismic data. This method is suitable to build high‐resolution velocity models provided that the targeted area is sampled by both diving waves and reflected waves. However, the conventional formulation of full‐waveform inversion prevents the reconstruction of the small wavenumber components of the velocity model when the subsurface is sampled by reflected waves only. This typically occurs as the depth becomes significant with respect to the length of the receiver array. This study first aims to highlight the limits of the conventional form of full‐waveform inversion when applied to seismic reflection data, through a simple canonical example of seismic imaging and to propose a new inversion workflow that overcomes these limitations. The governing idea is to decompose the subsurface model as a background part, which we seek to update and a singular part that corresponds to some prior knowledge of the reflectivity. Forcing this scale uncoupling in the full‐waveform inversion formalism brings out the transmitted wavepaths that connect the sources and receivers to the reflectors in the sensitivity kernel of the full‐waveform inversion, which is otherwise dominated by the migration impulse responses formed by the correlation of the downgoing direct wavefields coming from the shot and receiver positions. This transmission regime makes full‐waveform inversion amenable to the update of the long‐to‐intermediate wavelengths of the background model from the wide scattering‐angle information. However, we show that this prior knowledge of the reflectivity does not prevent the use of a suitable misfit measurement based on cross‐correlation, to avoid cycle‐skipping issues as well as a suitable inversion domain as the pseudo‐depth domain that allows us to preserve the invariant property of the zero‐offset time. This latter feature is useful to avoid updating the reflectivity information at each non‐linear iteration of the full‐waveform inversion, hence considerably reducing the computational cost of the entire workflow. Prior information of the reflectivity in the full‐waveform inversion formalism, a robust misfit function that prevents cycle‐skipping issues and a suitable inversion domain that preserves the seismic invariant are the three key ingredients that should ensure well‐posedness and computational efficiency of full‐waveform inversion algorithms for seismic reflection data.     

时间二阶积分波场的全波形反演

通过对波场的时间二阶积分运算以增强地震数据中的低频成分,提出了一种可有效减小对初始速度模型依赖性的地震数据全波形反演方法—时间二阶积分波场的全波形反演方法.根据散射理论中的散射波场传播方程,推导出时间二阶积分散射波场的传播方程,再利用一阶Born近似对时间二阶积分散射波场传播方程进行线性化.在时间二阶积分散射波场传播方程的基础上,利用散射波场反演地下散射源分布,再利用波场模拟的方法构建地下入射波场,然后根据时间二阶积分散射波场线性传播方程中散射波场与入射波场、速度扰动间的线性关系,应用类似偏移成像的公式得到速度扰动的估计,以此建立时间二阶积分波场的全波形迭代反演方法.最后把时间二阶积分波场的全波形反演结果作为常规全波形反演的初始模型可有效地减小地震波场全波形反演对初始模型的依赖性.应用于Marmousi模型的全频带合成数据和缺失4Hz以下频谱成分的缺低频合成数据验证所提出的全波形反演方法的正确性和有效性,数值试验显示缺失4Hz以下频谱成分数据的反演结果与全频带数据的反演结果没有明显差异.     

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.     

Multistep inversion workflow for 3D long‐offset damped elastic waves in the Fourier domain

We present a new workflow for imaging damped three‐dimensional elastic wavefields in the Fourier domain. The workflow employs a multiscale imaging approach, in which offset lengths are laddered, where frequency content and damping of the data are changed cyclically. Thus, the inversion process is launched using short‐offset and low‐frequency data to recover the long spatial wavelength of the image at a shallow depth. Increasing frequency and offset length leads to the recovery of the fine‐scale features of the model at greater depths. For the fixed offset, we employ (in the imaging process) a few discrete frequencies with a set of Laplace damping parameters. The forward problem is solved with a finite‐difference frequency‐domain method based on a massively parallel iterative solver. The inversion code is based upon the solution of a least squares optimisation problem and is solved using a nonlinear gradient method. It is fully parallelised for distributed memory computational platforms. Our full‐waveform inversion workflow is applied to the 3D Marmousi‐2 and SEG/EAGE Salt models with long‐offset data. The maximum inverted frequencies are 6 Hz for the Marmousi model and 2 Hz for the SEG/EAGE Salt model. The detailed structures are imaged successfully up to the depth approximately equal to one‐third of the maximum offset length at a resolution consistent with the inverted frequencies.     

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.     

基于精确震源函数的解调包络多尺度全波形反演

本文提出解调包络方法来重构地震记录中缺失的低频信号,同时该方法能够降低全波形反演的非线性程度;提出伴随状态震源函数反演方法来得到精确的震源函数,并推导了梯度计算公式;解调包络方法结合低通滤波技术,实现了从低频到高频的多尺度反演策略,有效缓解了全波形反演的周波跳跃问题.数值算例证明了解调包络、伴随状态震源函数反演方法和低通滤波多尺度反演策略的可行性及优越性.震源函数反演精度测试结果表明:即使观测记录在缺失低频信息的情况下,也能反演得到精确的震源函数.缺失低频测试和抗噪能力测试结果表明:即使地震数据中缺失9Hz以下的低频信号或者信噪比极低的情况下,利用反演得到的精确震源函数进行解调包络多尺度全波形反演,同样可以得到高精度的全波形反演结果.与Hilbert包络全波形反演对比结果表明:解调包络在重构低频和降低伴随震源主频方面具有一定优势.     

地震数据的反射波动方程最小二乘偏移

基于反射波动方程,本文提出了一种估计地下反射率分布的地震数据最小二乘偏移方法.高频近似下,非齐次的一次反射波动方程的源项是由反射率与入射波场的时间一阶导数相互作用产生的.根据反射波动方程,利用线性最小二乘反演方法由地震反射数据重建出地下产生反射波的反射源,再结合波场正演计算出的地下入射波场,得到地下反射率分布的估计.在地下反射源的线性最小二乘反演重建中,我们采用迭代求解方法,并以地震波的检波器单向地下照明强度作为最小二乘优化问题中Hessian矩阵的近似.     

Seismic envelope inversion: reduction of local minima and noise resistance

Waveform inversion met severe challenge in retrieving long‐wavelength background structure. We have proposed to use envelope inversion to recover the large‐scale component of the model. Using the large‐scale background recovered by envelope inversion as new starting model, we can get much better result than the conventional full waveform inversion. By comparing the configurations of the misfit functional between the envelope inversion and the conventional waveform inversion, we show that envelope inversion can greatly reduce the local minimum problem. The combination of envelope inversion and waveform inversion can deliver more faithful and accurate final result with almost no extra computation cost compared to the conventional full waveform inversion. We also tested the noise resistance ability of envelope inversion to Gaussian noise and seismic interference noise. The results showed that envelope inversion is insensitive to Gaussian noise and, to a certain extent, insensitive to seismic interference noise. This indicates the robustness of this method and its potential use for noisy data.