Wave-equation migration velocity analysis. II. Subsalt imaging examples

Subsalt imaging is strongly dependent on the quality of the velocity model. However, rugose salt bodies complicate wavefield propagation and lead to subsalt multipathing, illumination gaps and shadow zones, which cannot be handled correctly by conventional traveltime‐based migration velocity analysis (MVA). We overcome these limitations by the wave‐equation MVA technique, introduced in a companion paper, and demonstrate the methodology on a realistic synthetic data set simulating a salt‐dome environment and a Gulf of Mexico data set. We model subsalt propagation using wave paths created by one‐way wavefield extrapolation. Those wave paths are much more accurate and robust than broadband rays, since they inherit the frequency dependence and multipathing of the underlying wavefield. We formulate an objective function for optimization in the image space by relating an image perturbation to a perturbation of the velocity model. The image perturbations are defined using linearized prestack residual migration, thus ensuring stability, relative to the first‐order Born approximation assumptions. Synthetic and real data examples demonstrate that wave‐equation MVA is an effective tool for subsalt velocity analysis, even when shadows and illumination gaps are present.     

Extraction of amplitude-preserving angle gathers based on vector wavefield reverse-time migration

Angle-domain common-image gathers (ADCIGs) transformed from the shotdomain common-offset gathers are input to migration velocity analysis (MVA) and prestack inversion. ADCIGs are non-illusion prestack inversion gathers, and thus, accurate. We studied the extraction of elastic-wave ADCIGs based on amplitude-preserving elastic-wave reversetime migration for calculating the incidence angle of P-and S-waves at each image point and for different source locations. The P-and S-waves share the same incident angle, namely the incident angle of the source P-waves. The angle of incidence of the source P-wavefield was the difference between the source P-wave propagation angle and the reflector dips. The propagation angle of the source P-waves was obtained from the polarization vector of the decomposed P-waves. The reflectors’ normal direction angle was obtained using the complex wavenumber of the stacked reverse-time migration (RTM) images. The ADCIGs of P-and S-waves were obtained by rearranging the common-shot migration gathers based on the incident angle. We used a horizontally layered model, the graben medium model, and part of the Marmousi-II elastic model and field data to test the proposed algorithm. The results suggested that the proposed method can efficiently extract the P-and S-wave ADCIGs of the elastic-wave reverse-time migration, the P-and S-wave incident angle, and the angle-gather amplitude fidelity, and improve the MVA and prestack inversion.     

双平方根波动方程偏移速度分析

传统的剩余校正（RMO）偏移速度分析方法基于走时原理,在陡倾角和欠照明地区,因为不能得到充分的角度域信息而失效.本文将展示一种基于波场延拓理论的偏移速度分析方法,即波动方程偏移速度分析（WEMVA）.这种方法先利用成像优化方法获得剩余成像,再利用剩余成像反演剩余速度.此类方法继承了波动方程偏移方法的优点和缺点.波动方程偏移速度分析是一种线性反演方法,它要求对Born近似的展开序列作一阶截断.高阶部分的丢失必然带来巨大的截断误差,因此剩余成像必须也进行线性化,以适应大速度扰动和大延拓步长.因此,在此类算法中,剩余成像的获取和线性化是偏移速度分析的关键.在叠前偏移算子中,因为双平方根算子的数学表达式更为简洁,所以本文基于对波动方程偏移速度分析初步讨论,并通过模型验证其原理.     

Migration velocity analysis using pre‐stack wave fields

Using both image and data domains to perform velocity inversion can help us resolve the long and short wavelength components of the velocity model, usually in that order. This translates to integrating migration velocity analysis into full waveform inversion. The migration velocity analysis part of the inversion often requires computing extended images, which is expensive when using conventional methods. As a result, we use pre‐stack wavefield (the double‐square‐root formulation) extrapolation, which includes the extended information (subsurface offsets) naturally, to make the process far more efficient and stable. The combination of the forward and adjoint pre‐stack wavefields provides us with update options that can be easily conditioned to improve convergence. We specifically use a modified differential semblance operator to split the extended image into a residual part for classic differential semblance operator updates and the image (Born) modelling part, which provides reflections for higher resolution information. In our implementation, we invert for the velocity and the image simultaneously through a dual objective function. Applications to synthetic examples demonstrate the features of the approach.     

弹性矢量波成像域联合层析反演

在地震弹性矢量波场框架下,推导了多波联合层析速度反演方程以及走时残差与角道集剩余曲率的转换关系式,提出了一种利用成像域角道集更新P波、S波速度的走时层析反演方法.其实现过程可以概括为:将弹性波多分量数据作为输入,基于高斯束实现矢量波场成像并提取角道集,利用层析反演方程求解慢度更新量,最终获得多波联合反演结果.模型试算和实际资料处理验证了该方法的反演效果,能够为弹性矢量波联合深度偏移提供高质量的叠前速度场.     

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.     

Tomographic velocity inversion for ADCIGs in areas with a rugged surface

Pre-stack depth migration velocity analysis is one of the key techniques influencing image quality. As for areas with a rugged surface and complex subsurface, conventional prestack depth migration velocity analysis corrects the rugged surface to a known datum or designed surface velocity model on which to perform migration and update the velocity. We propose a rugged surface tomographic velocity inversion method based on angle-domain common image gathers by which the velocity field can be updated directly from the rugged surface without static correction for pre-stack data and improve inversion precision and efficiency. First, we introduce a method to acquire angle-domain common image gathers (ADCIGs) in rugged surface areas and then perform rugged surface tomographic velocity inversion. Tests with model and field data prove the method to be correct and effective.     

波动方程叠前深度偏移直接产生角道集

基于单程波深度延拓算法,得到震源波场在成像点的入射角度,结合对地层倾角的估计,获得入射地震波与界面法线的夹角.通过运用"保幅"的反褶积成像条件和考虑累加的炮数,解决了炮点覆盖不均匀导致的成像幅值误差问题,进而建议了炮域波动方程叠前深度偏移直接产生角道集的方法和流程.与基于空间移动或时移成像条件的波动方程叠前深度偏移提取角道集的方法相比,本文建议的方法只需少量额外的存储空间,又可补偿观测系统非均匀覆盖对成像幅值的影响;其增加的计算量与炮域偏移算法相比几乎可以忽略.文中算例表明,本文方法提取的角道集可为叠前反演提供较精确的AVO振幅特性.此外,就改善地震成像效果本身而言,提取角道集使得可在波动方程叠前深度偏移中应用剩余动校和拉伸切除技术,从而可更好地保持高频成分并提高成像的信噪比.     

Preserved‐amplitude angle domain migration by shot‐receiver wavefield continuation

We present preserved‐amplitude downward continuation migration formulas in the aperture angle domain. Our approach is based on shot‐receiver wavefield continuation. Since source and receiver points are close to the image point, a local homogeneous reference velocity can be approximated after redatuming. We analyse this approach in the framework of linearized inversion of Kirchhoff and Born approximations. From our analysis, preserved‐amplitude Kirchhoff and Born inverse formulas can be derived for the 2D case. They involve slant stacks of filtered subsurface offset domain common image gathers followed by the application of the appropriate weighting factors. For the numerical implementation of these formulas, we develop an algorithm based on the true amplitude version of the one‐way paraxial approximation. Finally, we demonstrate the relevance of our approach with a set of applications on synthetic datasets and compare our results with those obtained on the Marmousi model by multi‐arrival ray‐based preserved‐amplitude migration. While results are similar, we observe that our results are less affected by artefacts.     

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

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

3D angle-domain common-image gathers for migration velocity analysis

Angle‐domain common‐image gathers (ADCIGs) are an essential tool for migration velocity analysis (MVA). We present a method for computing ADCIGs in 3D from the results of wavefield‐continuation migration. The proposed methodology can be applied before or after the imaging step in a migration procedure. When computed before imaging, 3D ADCIGs are functions of the offset ray parameters (p, p) ; we derive the geometric relationship that links the offset ray parameters to the aperture angle γ and the reflection azimuth φ. When computed after imaging, 3D ADCIGs are directly produced as functions of γ and φ. The mapping of the offset ray parameters (p, p) into the angles (γ, φ) depends on both the local dips and the local interval velocity; therefore, the transformation of ADCIGs computed before imaging into ADCIGs that are functions of the actual angles is difficult in complex structure. By contrast, the computation of ADCIGs after imaging is efficient and accurate even in the presence of complex structure and a heterogeneous velocity function. On the other hand, the estimation of the offset ray parameters (p, p) is less sensitive to velocity errors than the estimation of the angles (γ, φ). When ADCIGs that are functions of the offset ray parameters (p, p) are adequate for the application of interest (e.g. ray‐based tomography), the computation of ADCIGs before imaging might be preferable. Errors in the migration velocity cause the image point in the angle domain to shift along the normal to the apparent geological dip. By assuming stationary rays (i.e. small velocity errors), we derive a quantitative relationship between this normal shift and the traveltime perturbation caused by velocity errors. This relationship can be directly used in an MVA procedure to invert depth errors measured from ADCIGs into migration velocity updates. In this paper, we use it to derive an approximate 3D residual moveout (RMO) function for measuring inconsistencies between the migrated images at different γ and φ. We tested the accuracy of our kinematic analysis on a 3D synthetic data set with steeply dipping reflectors and a vertically varying propagation velocity. The tests confirm the accuracy of our analysis and illustrate the limitations of the straight‐ray approximation underlying our derivation of the 3D RMO function.     

The effectiveness of a pseudo-inverse extended Born operator to handle lateral heterogeneity for imaging and velocity analysis applications

Wave equation–based migration velocity analysis techniques aim to construct a kinematically accurate velocity model for imaging or as an initial model for full waveform inversion applications. The most popular wave equation–based migration velocity analysis method is differential semblance optimization, where the velocity model is iteratively updated by minimizing the unfocused energy in an extended image volume. However, differential semblance optimization suffers from artefacts, courtesy of the adjoint operator used in imaging, leading to poor convergence. Recent findings show that true amplitude imaging plays a significant role in enhancing the differential semblance optimization's gradient and reducing the artefacts. Here, we focus on a pseudo-inverse operator to the horizontally extended Born as a true amplitude imaging operator. For laterally inhomogeneous models, the operator required a derivative with respect to a vertical shift. Extending the image vertically to evaluate such a derivative is costly and impractical. The inverse operator can be simplified in laterally homogeneous models. We derive an extension of the approach to apply the full inverse formula and evaluate the derivative efficiently. We simplified the implementation by applying the derivative to the imaging condition and utilize the relationship between the source and receiver wavefields and the vertical shift. Specifically, we verify the effectiveness of the approach using the Marmousi model and show that the term required for the lateral inhomogeneity treatment has a relatively small impact on the results for many cases. We then apply the operator in differential semblance optimization and invert for an accurate macro-velocity model, which can serve as an initial velocity model for full waveform inversion.     

单程波角度域内压制多次波偏移假象

多次波偏移中的假象主要来自于不同地震事件之间的互相关,由于这种互相关满足成像条件,很难直接在偏移过程中去除.但是对于准确的速度模型,真实的成像结果在角度域内应该是平直的.根据这个判断准则,可以在角度域内移除多次波偏移中的假象.本文以数据自相关偏移为例,提出了在单程波多次波偏移中移除假象的主要流程:首先在在单程波偏移过程中高效地提取角度域共成像点道集,然后对角度域共成像点道集应用高分辨率的抛物线型Radon变换,用合适的切除函数处理后,反变换回到角度域,最后叠加各个角度成分,得到偏移结果.Marmousi模型的合成数据测试表明,这种方法可以很好地压制多次波偏移过程中产生的假象,有效地提高成像结果的信噪比.     

Image‐ray Tomography

Seismic tomography is a well‐established approach to invert smooth macro‐velocity models from kinematic parameters, such as traveltimes and their derivatives, which can be directly estimated from data. Tomographic methods differ more with respect to data domains than in the specifications of inverse‐problem solving schemes. Typical examples are stereotomography, which is applied to prestack data and Normal‐Incidence‐Point‐wave tomography, which is applied to common midpoint stacked data. One of the main challenges within the tomographic approach is the reliable estimation of the kinematic attributes from the data that are used in the inversion process. Estimations in the prestack domain (weak and noisy signals), as well as in the post‐stack domain (occurrence of triplications and diffractions leading to numerous conflicting dip situations) may lead to parameter inaccuracies that will adversely impact the resulting velocity models. To overcome the above limitations, a new tomographic procedure applied in the time‐migrated domain is proposed. We call this method Image‐Incident‐Point‐wave tomography. The new scheme can be seen as an alternative to Normal‐Incidence‐Point‐wave tomography. The latter method is based on traveltime attributes associated with normal rays, whereas the Image‐Incidence‐Point‐wave technique is based on the corresponding quantities for the image rays. Compared to Normal‐Incidence‐Point‐wave tomography the proposed method eases the selection of the tomography attributes, which is shown by synthetic and field data examples. Moreover, the method provides a direct way to convert time‐migration velocities into depth‐migration velocities without the need of any Dix‐style inversion.     

Wavefield tomography based on local image correlations

The estimation of a velocity model from seismic data is a crucial step for obtaining a high‐quality image of the subsurface. Velocity estimation is usually formulated as an optimization problem where an objective function measures the mismatch between synthetic and recorded wavefields and its gradient is used to update the model. The objective function can be defined in the data‐space (as in full‐waveform inversion) or in the image space (as in migration velocity analysis). In general, the latter leads to smooth objective functions, which are monomodal in a wider basin about the global minimum compared to the objective functions defined in the data‐space. Nonetheless, migration velocity analysis requires construction of common‐image gathers at fixed spatial locations and subsampling of the image in order to assess the consistency between the trial velocity model and the observed data. We present an objective function that extracts the velocity error information directly in the image domain without analysing the information in common‐image gathers. In order to include the full complexity of the wavefield in the velocity estimation algorithm, we consider a two‐way (as opposed to one‐way) wave operator, we do not linearize the imaging operator with respect to the model parameters (as in linearized wave‐equation migration velocity analysis) and compute the gradient of the objective function using the adjoint‐state method. We illustrate our methodology with a few synthetic examples and test it on a real 2D marine streamer data set.     

Residual curvature migration velocity analysis for angle domain common imaging gathers

Pre-stack depth migration velocity analysis is one of the keys to influencing the imaging quality of pre-stack migration. In this paper we cover a residual curvature velocity analysis method on angle-domain common image gathers (ADCIGs) which can depict the relationship between incident angle and migration depth at imaging points and update the migration velocity. Differing from offset-domain common image gathers (ODCIGs), ADCIGs are not disturbed by the multi-path problem which contributes to imaging artifacts, thus influencing the velocity analysis. On the basis of horizontal layers, we derive the residual depth equation and also propose a velocity analysis workflow for velocity scanning. The tests to synthetic and field data prove the velocity analysis methods adopted in this paper are robust and valid.     

Wave‐equation migration velocity analysis with time‐shift imaging

Wave‐equation migration velocity analysis is a technique designed to extract and update velocity information from migrated images. The velocity model is updated through the process of optimizing the coherence of images migrated with the known background velocity model. The capacity for handling multi‐pathing of the technique makes it appropriate in complex subsurface regions characterized by strong velocity variation. Wave‐equation migration velocity analysis operates by establishing a linear relation between a slowness perturbation and a corresponding image perturbation. The linear relationship and the corresponding linearized operator are derived from conventional extrapolation operators and the linearized operator inherits the main properties of frequency‐domain wavefield extrapolation. A key step in the implementation is to design an appropriate procedure for constructing an image perturbation relative to a reference image that represents the difference between the current image and a true, or more correct image of the subsurface geology. The target of the inversion is to minimize such an image perturbation by optimizing the velocity model. Using time‐shift common‐image gathers, one can characterize the imperfections of migrated images by defining the focusing error as the shift of the focus of reflections along the time‐shift axis. The focusing error is then transformed into an image perturbation by focusing analysis under the linear approximation. As the focusing error is caused by the incorrect velocity model, the resulting image perturbation can be considered as a mapping of the velocity model error in the image space. Such an approach for constructing the image perturbation is computationally efficient and simple to implement. The technique also provides a new alternative for using focusing information in wavefield‐based velocity model building. Synthetic examples demonstrate the successful application of our method to a layered model and a subsalt velocity update problem.     

3D post‐stack interval velocity analysis with effective use of datuming

Interval velocity analysis using post‐stack data has always been a desire, mainly for 3D data sets. In this study we present a method that uses the unique characteristics of migrated diffractions to enable interval velocity analysis from three‐dimensional zero‐offset time data. The idea is to perform a standard three‐dimensional prestack depth migration on stack cubes and generate three‐dimensional common image gathers that show great sensitivity to velocity errors. An efficient ‘top‐down’ scheme for updating the velocity is used to build the model. The effectiveness of the method is related to the incorporation of wave equation based post‐stack datuming in the model building process. The proposed method relies on the ability to identify diffractions along redatumed zero‐offset data and to analyse their flatness in the migrated local angle domain. The method can be considered as an additional tool for a complete, prestack depth migration based interval velocity analysis.     

2.5D modelling, inversion and angle migration in anisotropic elastic media

2.5D modelling approximates 3D wave propagation in the dip‐direction of a 2D geological model. Attention is restricted to raypaths for waves propagating in a plane. In this way, fast inversion or migration can be performed. For velocity analysis, this reduction of the problem is particularly useful. We review 2.5D modelling for Born volume scattering and Born–Helmholtz surface scattering. The amplitudes are corrected for 3D wave propagation, taking into account both in‐plane and out‐of‐plane geometrical spreading. We also derive some new inversion/migration results. An AVA‐compensated migration routine is presented that is simplified compared with earlier results. This formula can be used to create common‐image gathers for use in velocity analysis by studying the residual moveout. We also give a migration formula for the energy‐flux‐normalized plane‐wave reflection coefficient that models large contrast in the medium parameters not treated by the Born and the Born–Helmholtz equation results. All results are derived using the generalized Radon transform (GRT) directly in the natural coordinate system characterized by scattering angle and migration dip. Consequently, no Jacobians are needed in their calculation. Inversion and migration in an orthorhombic medium or a transversely isotropic (TI) medium with tilted symmetry axis are the lowest symmetries for practical purposes (symmetry axis is in the plane). We give an analysis, using derived methods, of the parameters for these two types of media used in velocity analysis, inversion and migration. The kinematics of the two media involve the same parameters, hence there is no distinction when carrying out velocity analysis. The in‐plane scattering coefficient, used in the inversion and migration, also depends on the same parameters for both media. The out‐of‐plane geometrical spreading, necessary for amplitude‐preserving computations, for the TI medium is dependent on the same parameters that govern in‐plane kinematics. For orthorhombic media, information on additional parameters is required that is not needed for in‐plane kinematics and the scattering coefficients. Resolution analysis of the scattering coefficient suggests that direct inversion by GRT yields unreliable parameter estimates. A more practical approach to inversion is amplitude‐preserving migration followed by AVA analysis. SYMBOLS AND NOTATION A list of symbols and notation is given in Appendix D .     

Migration for velocity and attenuation perturbations

Migration maps seismic data to reflectors in the Earth. Reflections are not only caused by small‐scale variations of the velocity and density but also of the quality factor that describes attenuation. We investigated scattering due to velocity and attenuation perturbations by computing the resolution function or point‐spread function in a homogeneous background model. The resolution function is the migration image of seismic reflection data generated by a point scatterer. We found that the resolution function mixes velocity and attenuation parameter perturbations to the extent that they cannot be reconstructed independently. This is true for a typical seismic setting with sources and receivers at the surface and a buried scatterer. As a result, it will be impossible to simultaneously invert for velocity and attenuation perturbations in the scattering approach, also known as the Born approximation. We proceeded to investigate other acquisition geometries that may resolve the ambiguity between velocity and attenuation perturbations. With sources and receivers on a circle around the scatterer, in 2D, the ambiguity disappears. It still shows up in a cross‐well setting, although the mixing of velocity and attenuation parameters is less severe than in the surface‐to‐surface case. We also consider illumination of the target by diving waves in a background model that has velocity increasing linearly with depth. The improvement in illumination is, however, still insufficient to remove the ambiguity.