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.     

基于深度降噪自编码神经网络的中国大陆地壳厚度反演

本文采用基于数据驱动的深度降噪自编码网络构建了瑞雷面波群速度、相速度频散特性与地壳厚度的正反演函数关系,并利用最新频散模型反演了中国大陆的地壳厚度.对于神经网络架构体系的评价,除了考虑传统意义上的测试误差、训练误差之外,本文还用已知物理原理的正演结果与网络预测结果进行比较;在设计网络构架时,同时考虑地球模型和面波频散的...     

Multiscale linearized inversion of surface-wave dispersion curves

In recent years, surface-wave analysis method has been developed rapidly in many fields. Multichannel analysis of surface waves can provide near-surface one-dimensional shear-wave velocity profiles. Because linearized inversion of surface-wave dispersion curves relies heavily on the choice of the initial model, setting an inappropriate initial model can lead to poor inversion results, or even failure of inversion. However, it is difficult to establish a reasonable initial model without a priori information, which is unavailable in most cases. To cope with this problem, a multiscale linearized inversion method is proposed for surface-wave dispersion curves inversion. In contrast with the traditional single-scale linearized inversion, the key idea of the proposed multiscale surface-wave inversion method is the introduction of a merging and splitting process of layers. After every scale inversion, the merging and splitting operations automatically optimize the inversion model, making it gradually approach to a reasonable subsurface stratification. Multiscale surface-wave inversion method reduces the difficulty of establishing the initial model and has high computational efficiency. In addition, it has strong ability to identify high-velocity or low-velocity interlayers and thin layers, especially suited for the geological conditions with obvious stratification. In synthetic tests, the proposed method was compared with the single-scale surface-wave inversion and particle swarm optimization algorithm to demonstrate the effectiveness and practicability of multiscale surface-wave inversion method. We also applied the multiscale surface-wave inversion method to field seismic data acquired in Guizhou, China and Texas, USA. Borehole and crosshole test data were compared with the inversion results of field data to prove the reliability of the proposed method.     

Global surface wave inversion with model constraints‡

Practical applications of surface wave inversion demand reliable inverted shear‐wave profiles and a rigorous assessment of the uncertainty associated to the inverted parameters. As a matter of fact, the surface wave inverse problem is severely affected by solution non‐uniqueness: the degree of non‐uniqueness is closely related to the complexity of the observed dispersion pattern and to the experimental inaccuracies in dispersion measurements. Moreover, inversion pitfalls may be connected to specific problems such as inadequate model parametrization and incorrect identification of the surface wave modes. Consequently, it is essential to tune the inversion problem to the specific dataset under examination to avoid unnecessary computations and possible misinterpretations. In the heuristic inversion algorithm presented in this paper, different types of model constraints can be easily introduced to bias constructively the solution towards realistic estimates of the 1D shear‐wave profile. This approach merges the advantages of global inversion, like the extended exploration of the parameter space and a theoretically rigorous assessment of the uncertainties on the inverted parameters, with the practical approach of Lagrange multipliers, which is often used in deterministic inversion, which helps inversion to converge towards models with desired properties (e.g., ‘smooth’ or ‘minimum norm' models). In addition, two different forward kernels can be alternatively selected for direct‐problem computations: either the conventional modal inversion or, instead, the direct minimization of the secular function, which allows the interpreter to avoid mode identification. A rigorous uncertainty assessment of the model parameters is performed by posterior covariance analysis on the accepted solutions and the modal superposition associated to the inverted models is investigated by full‐waveform modelling. This way, the interpreter has several tools to address the more probable sources of inversion pitfalls within the framework of a rigorous and well‐tested global inversion algorithm. The effectiveness and the versatility of this approach, as well as the impact of the interpreter's choices on the final solution and on its posterior uncertainty, are illustrated using both synthetic and real data. In the latter case, the inverted shear velocity profiles are blind compared with borehole data.     

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.     

Estimation of shear wave velocity profiles by the inversion of spatial autocorrelation coefficients

The subsurface shear-wave velocity (Vs) is considered to be a key parameter for site characterization and assessment of earthquake hazard because of its great influence on local ground-motion amplification. Array microtremor measurements are widely used for the estimation of shear-wave velocities. Compared to other methods such as frequency-wavenumber (f-k) methods, the spatial autocorrelation (SPAC) method requires fewer sensors and thus is relatively easier to implement and gives robust estimations of shear-wave velocity profiles for depths down to a few hundred meters. The quantity derived from observed data is the SPAC coefficient, which is a function of correlation distance, frequency and phase velocity. Generally, estimation of Vs profiles is a two stage process: Estimation of the dispersion data from the SPAC coefficients and inversion of the dispersion data for shear-wave velocity structure. In this study, instead of inverting dispersion curves, a more practical approach is used; that is, observed SPAC coefficients are directly inverted for the S-wave velocities. A synthetic case and a field data application are presented to test the potential of the inversion algorithm. We obtain an iterative damped least-squares solution with differential smoothing. The differential smoothing approach constrains the change in shear-wave velocities of the adjacent layers and thus stabilizes the inversion.     

基于初至波走时层析成像的Tikhonov正则化与梯度优化算法

初至波走时层析成像是利用地震初至波走时和其传播的射线路径来反演地下介质速度的技术.该问题本质上是一个不适定问题,需要使用正则化方法并辅之以适当的最优化技巧.本文从数值优化的角度介绍了初至波走时层析成像的反演原理,建立了Tikhonov正则化层析成像反演模型并提出求解极小化问题的加权修正步长的梯度下降算法.该方法可以从速度模型的可行域中迭代找到一个最优解.数值试验表明,该方法是可行和有应用前景的.     

瑞雷波多阶模式频散曲线稀疏正则化反演方法研究

目前瑞雷波多阶模式频散曲线反演中仅考虑数据的拟合,缺乏对模型的约束,不能很好地刻画地层间断面的问题,针对此问题,研究了瑞雷波多阶模式频散曲线稀疏正则化反演方法.正演模拟基于广义反射-透射系数法,数值计算上采用一种快速求根方法,与二等分方法相比,能够在很短的时间内达到最优的收敛效果;反演建模时采用L1范数正则化方法对模型...     

Non-uniqueness in surface-wave inversion and consequences on seismic site response analyses

Surface-wave tests are based on the solution of an inverse problem for shear-wave velocity profile identification from the experimentally measured dispersion curve. The main criticisms for these testing methodologies are related to the inverse problem solution and arise from the possible equivalence of different shear-wave velocity profiles. In this paper, some implications of solution non-uniqueness for seismic response studies are investigated using both numerical simulations and experimental data. A Monte Carlo approach for the inversion problem has been used to obtain a set of equivalent shear-wave velocity models. This selection is based on a statistical test which takes into account both data uncertainty and model parameterization. This set of solutions (i.e., soil profiles) is then used to evaluate the seismic response with a conventional one-dimensional analysis. It is shown that equivalent profiles with respect to surface-wave testing are equivalent also with respect to site amplification, thus countering the criticism related to inversion uncertainty for the engineering use of surface-wave tests.     

Stepwise joint inversion of surface wave dispersion,Rayleigh wave ZH ratio,and receiver function data for 1D crustal shear wave velocity structure

Accurate determination of seismic velocity of the crust is important for understanding regional tectonics and crustal evolution of the Earth. We propose a stepwise joint linearized inversion method using surface wave dispersion, Rayleigh wave ZH ratio (i.e., ellipticity), and receiver function data to better resolve 1D crustal shear wave velocity (v _S) structure. Surface wave dispersion and Rayleigh wave ZH ratio data are more sensitive to absolute variations of shear wave speed at depths, but their sensitivity kernels to shear wave speeds are different and complimentary. However, receiver function data are more sensitive to sharp velocity contrast (e.g., due to the existence of crustal interfaces) and v _P/v _S ratios. The stepwise inversion method takes advantages of the complementary sensitivities of each dataset to better constrain the v _S model in the crust. We firstly invert surface wave dispersion and ZH ratio data to obtain a 1D smooth absolute v _S model and then incorporate receiver function data in the joint inversion to obtain a finer v _S model with better constraints on interface structures. Through synthetic tests, Monte Carlo error analyses, and application to real data, we demonstrate that the proposed joint inversion method can resolve robust crustal v _S structures and with little initial model dependency.     

扩展成像条件下的最小二乘逆时偏移

逆时偏移(RTM)是复杂介质条件下地震成像的重要手段.因受观测系统限制、上覆地层影响以及波场带宽有限等因素的影响,现行的常规RTM所采用的互相关成像条件通常对地下构造进行模糊成像.最小二乘逆时偏移(LSRTM)通过最小化线性Born近似正演数据和采集数据之间的波形差异,采用梯度类反演算法优化反射系数模型,获得的成像结果具有更高的分辨率和更可靠的振幅保真度.然而,基于波形拟合的LSRTM对背景速度模型的依赖性很强.误差太大的速度模型容易产生周波跳跃现象,导致LSRTM难以获得全局最优解.为了克服这一问题,本文基于扩展模型的思想,在线性Born近似下,推导得到RTM扩展成像条件.并基于最小二乘反演理论,提出扩展成像条件下的LSRTM方法.理论模型试算表明,本文方法不仅可以提供分辨率更高、振幅属性更为可靠的成像结果,而且能够在一定程度上消除速度误差对反演成像的影响.     

横向非均匀介质多道瑞雷波频散曲线正演

作为近地表横波速度结构成像的主要手段之一,面波多道分析法的正问题研究对现场观测系统设计及后续反演计算具有重要意义.目前面波频散曲线的正演主要分为两类:一是对水平层状介质中面波的本征值问题进行求解,该类方法计算效率高但较难考虑地下介质在横向上的不均匀性;二是基于波动方程的全波场模拟,该类方法在理论上可考虑任意复杂的地质模型但计算成本相对较高.本文基于振幅归一化加权的聚束分析,提出了一种适用于横向非均匀介质模型的多道瑞雷波频散曲线正演方法.首先,基于聚束分析的计算公式推导得到了经振幅归一化加权后输出功率谱中相速度与局部相速度之间的关系,然后通过黄金分割极值搜索算法计算得到了多道瑞雷波数据的理论频散曲线.数值分析结果表明,该算法能够快速地实现横向非均匀介质中多道瑞雷波频散曲线的正演计算,所求取的频散曲线与采用二维弹性波时间域有限差分模拟分析得到的结果误差较小,这在一定程度上说明了该计算方法的可靠性,从而可为面波多道分析法中的观测系统快速优化设计以及横向非均匀介质中频散曲线的反演解释提供理论支撑.     

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.     

Rayleigh wave dispersion curve inversion via genetic algorithms and Marginal Posterior Probability Density estimation

Surface wave dispersion curve inversion is a challenging problem for linear inversion procedures due to its highly non-linear nature and to the large numbers of local minima and maxima of the objective function (multi-modality). In order to improve the reliability of the inversion results, we implemented and tested a two-step inversion scheme based on Genetic Algorithms (GAs). The proposed scheme performs several preliminary “parallel” runs (first step) and a final global run using the previously-determined fittest models as starting population.In this work we focus on the inversion of shear-wave velocity and layer thickness while fixing compressional-wave velocity and density according to user-defined Poisson's ratios and velocity–density relationship respectively. The procedure can nonetheless perform the inversion under different degrees of regularization, depending on the a priori information and the desired degree of freedom of the system.Thanks to the large number of considered models, in addition to the fittest model, a mean model and its accuracy are evaluated by means of a statistical approach based on the estimation of the Marginal Posterior Probability Density (MPPD).We tested the proposed GA-based inversion scheme on three synthetic models reproducing a complex structure with low-to-moderate velocity cover (also including a low-velocity channel) lying over hard bedrock. For all the considered cases the bedrock velocity and depth were properly identified, and velocity inversion was reconstructed with minor uncertainties.The performed tests also investigate the influence of the first higher mode, the reduction of the frequency range of the considered dispersion curve as well as the use of different number of strata. While a limited frequency range of the dispersion curve (maximum frequency reduced from 80 to 40 Hz) does not seem to significantly limit the accuracy of the retrieved model, the adoption of the correct number of strata and the addition of the first higher mode help better focus the final solution.In conclusion, the proposed approach represents an improvement of a purely GA-based optimization scheme and the MPPD-based mean model typically offers a more significant and precise solution than the fittest one.Results of the inversion performed on a field data set were validated by borehole stratigraphy.     

Sparse Bayesian linearized amplitude-versus-angle inversion

The technique of seismic amplitude-versus-angle inversion has been widely used to estimate lithology and fluid properties in seismic exploration. The amplitude-versus-angle inversion problem is intrinsically ill-posed and generally stabilized by the use of L2-norm regularization methods but with drawback of smoothing important boundaries between adjacent layers. In this study, we propose a sparse Bayesian linearized solution for amplitude-versus-angle inversion problem to preserve the sharp geological interfaces. In this regard, a priori constraint term with two regularization functions is presented: the sparse constraint regularization and the low-frequency model information. In addition, to obtain high-resolution reflectivity estimation, the model parameters decorrelation technique combined with dipole decomposition method is employed. We validate the applicability of the presented method by both synthetic and real seismic data from the Gulf of Mexico. The accuracy improvement of the presented method is also confirmed by comparing the results with the commonly used Bayesian linearized amplitude-versus-angle inversion.     

基于面波频散的三维横波速度方位各向异性层析成像方法

地震各向异性是反映地球内部介质特性的重要指针之一。常用的横波分裂法和二维面波方位各向异性层析成像方法很难准确反映各向异性随深度的变化。将与周期相关的区域化面波方位各向异性转换成与深度相关的一维横波速度方位各向异性可以弥补深度信息不足的缺陷。现有三维横波速度各向异性研究多是通过两步方法来实现的,即逐个周期二维面波方位各向异性层析成像以及逐个格点一维横波速度方位各向异性反演。这种分步反演的方式既不利于三维先验约束的引入,也不利于利用原始观测拟合误差对三维模型进行直接评估。因此本文开发了基于面波频散曲线的三维横波速度方位各向异性层析成像方法,并编制了相关正演和反演程序。为了检测方法和程序的有效性,我们对规律分布的三维检测板模型进行了模拟测试。测试结果显示:该方法可以很好地恢复各向同性波速异常、各向异性相对强度和快波方向等三维结构信息;而且反演模型相对于参考模型明显改善了对观测数据的拟合,降低了对观测数据的均方根误差。但对各向同性理论模型进行各向异性反演时,在波速均匀区可产生小于0.5%的假各向异性幅值,在波速非均匀区该假的各向异性幅值会更大,浅部可达3.5%。因此在实际应用中需要谨慎解释（浅部）非均匀区的各向异性结果。      

多震相层析成像与在介质中的多波型射线跟踪

本文介绍了多震相的层析成像的思路和算法，通过穿透和反射走时可以同时作出２维和３维慢度（速度的倒数）重建。我们分析了在穿透和反射数据中确定速度和深度的不确定性，并认识到深度扰动对反射走时异常比慢度扰动更敏感。由不同波类型所提供的对速度和深度的约束，这个算法实际上减少了在一般反射层析成像在速度和深度之间的不确定性，并且也避免了在穿透层析成像中的不确定问题。线性化反演是通过从反射界面深度由分离速度参数迭代进行的。使用一个快速的２－Ｄ和３－Ｄ射线跟踪算法来计算穿透和反射走时和对幔度及反射界面深度的偏导数。深度和速度都用立方Ｂ样条函数来进行参数化。合成例子表明，当同时考虑穿透和反射时间，层析成像的结果得到改进。这个方法也应用到英国煤炭测量局（ＢｒｉｔｉｓｈＣｏａｌＭｅａｓｕｒｅｓ）沿跨线排列所记录的逆ＶＳＰ数据组。通过使用波形配合技术，用同时确定时间延迟和叠加权，可以自动拾取旅行时间。所观察到的逆ＶＳＰ层析成像可比周围介质具有较低速度的两个断层区域成像。断层的位置由附近的反射测线所确定。本文还讨论了在复杂２－Ｄ和３－Ｄ非均匀各向同性介质中地震射线跟踪方法。界面的几何形状和水平速度场都通过使用非均匀步长立方Ｂ－样条节点     

Full waveform inversion and the truncated Newton method: quantitative imaging of complex subsurface structures

Full waveform inversion is a powerful tool for quantitative seismic imaging from wide‐azimuth seismic data. The method is based on the minimization of the misfit between observed and simulated data. This amounts to the solution of a large‐scale nonlinear minimization problem. The inverse Hessian operator plays a crucial role in this reconstruction process. Accounting accurately for the effect of this operator within the minimization scheme should correct for illumination deficits, restore the amplitude of the subsurface parameters, and help to remove artefacts generated by energetic multiple reflections. Conventional minimization methods (nonlinear conjugate gradient, quasi‐Newton methods) only roughly approximate the effect of this operator. In this study, we are interested in the truncated Newton minimization method. These methods are based on the computation of the model update through a matrix‐free conjugate gradient solution of the Newton linear system. We present a feasible implementation of this method for the full waveform inversion problem, based on a second‐order adjoint state formulation for the computation of Hessian‐vector products. We compare this method with conventional methods within the context of 2D acoustic frequency full waveform inversion for the reconstruction of P‐wave velocity models. Two test cases are investigated. The first is the synthetic BP 2004 model, representative of the Gulf of Mexico geology with high velocity contrasts associated with the presence of salt structures. The second is a 2D real data‐set from the Valhall oil field in North sea. Although, from a computational cost point of view, the truncated Newton method appears to be more expensive than conventional optimization algorithms, the results emphasize its increased robustness. A better reconstruction of the P‐wave velocity model is provided when energetic multiple reflections make it difficult to interpret the seismic data. A better trade‐off between regularization and resolution is obtained when noise contamination of the data requires one to regularize the solution of the inverse problem.     

Middle and upper crust shear-wave velocity structure of the Chinese mainland

In order to give a more reliable shallow crust model for the Chinese mainland, the present study collected many short-period surface wave data which are better sensitive to shallow earth structures. Different from traditional two-step surface wave tomography, we developed a new linearized surface wave dispersion inversion method to directly get a 3D S-wave velocity model in the second step instead of inverting for 1D S-velocity profile cell by cell. We convert all the regionalized dispersions into linear constraints for a 3D S-velocity model. Checkerboard tests show that this method can give reasonable results. The distribution of the middle-and upper-crust shear-wave velocity of the Chinese mainland in our model is strongly heterogeneous and related to different geotectonic terrains. Low-velocity anomalies delineated very well most of the major sedimentary basins of China. And the variation of velocities at different depths gives an indication of basement depth of the basins. The western Tethyan tectonic domain (on the west of the 95°E longitude) is characterized by low velocity, while the eastern Tethyan domain does not show obvious low velocity. Since petroleum resources often distribute in sedimentary basins where low-velocity anomaly appears, the low velocity anomalies in the western Tethyan domain may indicate a better petroleum prospect than in its eastern counterpart. Besides, low velocity anomaly in the western Tethyan domain and around the Xing’an orogenic belt may be partly caused by high crustal temperature. The weak low-velocity belt along ~105°E longitude corresponds to the N-S strong seismic belt of central China.     

频率多尺度全波形速度反演

以二维声波方程为模型,在时间域深入研究了全波形速度反演.全波形反演要解一个非线性的最小二乘问题,是一个极小化模拟数据与已知数据之间残量的过程.针对全波形反演易陷入局部极值的困难,本文提出了基于不同尺度的频率数据的"逐级反演"策略,即先基于低频尺度的波场信息进行反演,得出一个合理的初始模型,然后再利用其他不同尺度频率的波场进行反演,并且用前一尺度的迭代反演结果作为下一尺度反演的初始模型,这样逐级进行反演.文中详细阐述和推导了理论方法及公式,包括有限差分正演模拟、速度模型修正、梯度计算和算法描述,并以Marmousi复杂构造模型为例,进行了MPI并行全波形反演数值计算,得到了较好的反演结果,验证了方法的有效性和稳健性.