Investigating a time‐shift extended imaging condition in a Kirchhoff pre‐stack depth migration algorithm

Extracting true amplitude versus angle common image gathers is one of the key objectives in seismic processing and imaging. This is achievable to different degrees using different migration techniques (e.g., Kirchhoff, wavefield extrapolation, and reverse time migration techniques) and is a common tool in exploration, but the costs can vary depending on the selected migration algorithm and the desired accuracy. Here, we investigate the possibility of combining the local‐shift imaging condition, specifically the time‐shift extended imaging condition, for angle gathers with a Kirchhoff migration. The aims are not to replace the more accurate full‐wavefield migration but to offer a cheaper alternative where ray‐based methods are applicable and to use Kirchhoff time‐lag common image gathers to help bridge the gap between the traditional offset common image gathers and reverse time migration angle gathers; finally, given the higher level of summation inside the extended imaging migration, we wish to understand the impact on the amplitude versus angle response. The implementation of the time‐shift imaging condition along with the computational cost is discussed, and results of four different datasets are presented. The four example datasets, two synthetic, one land acquisition, and a marine dataset, have been migrated using a Kirchhoff offset method, a Kirchhoff time‐shift method, and, for comparison, a reverse time migration algorithm. The results show that the time‐shift imaging condition at zero time lag is equivalent to the full offset stack as expected. The output gathers are cleaner and more consistent in the time‐lag‐derived angle gathers, but the conversion from time lag to angle can be considered a post‐processing step. The main difference arises in the amplitude versus offset/angle distribution where the responses are different and dramatically so for the land data. The results from the synthetics and real data show that a Kirchhoff migration with an extended imaging condition is capable of generating subsurface angle gathers. The same disadvantages with a ray‐based approach will apply using the extended imaging condition relative to a wave equation angle gather solution. Nevertheless, using this approach allows one to explore the relationship between the velocity model and focusing of the reflected energy, to use the Radon transformation to remove noise and multiples, and to generate consistent products from a ray‐based migration and a full‐wave equation migration, which can then be interchanged depending on the process under study.     

Wavefield interpolation in 3D large‐step Fourier wavefield extrapolation

Extrapolating wavefields and imaging at each depth during three‐dimensional recursive wave‐equation migration is a time‐consuming endeavor. For efficiency, most commercial techniques extrapolate wavefields through thick slabs followed by wavefield interpolation within each thick slab. In this article, we develop this strategy by associating more efficient interpolators with a Fourier‐transform‐related wavefield extrapolation method. First, we formulate a three‐dimensional first‐order separation‐of‐variables screen propagator for large‐step wavefield extrapolation, which allows for wide‐angle propagations in highly contrasting media. This propagator significantly improves the performance of the split‐step Fourier method in dealing with significant lateral heterogeneities at the cost of only one more fast Fourier transform in each thick slab. We then extend the two‐dimensional Kirchhoff and Born–Kirchhoff local wavefield interpolators to three‐dimensional cases for each slab. The three‐dimensional Kirchhoff interpolator is based on the traditional Kirchhoff formula and applies to moderate lateral velocity variations, whereas the three‐dimensional Born–Kirchhoff interpolator is derived from the Lippmann–Schwinger integral equation under the Born approximation and is adapted to highly laterally varying media. Numerical examples on the three‐dimensional salt model of the Society of Exploration Geophysicists/European Association of Geoscientists demonstrate that three‐dimensional first‐order separation‐of‐variables screen propagator Born–Kirchhoff depth migration using thick‐slab wavefield extrapolation plus thin‐slab interpolation tolerates a considerable depth‐step size of up to 72 ms, eventually resulting in an efficiency improvement of nearly 80% without obvious loss of imaging accuracy. Although the proposed three‐dimensional interpolators are presented with one‐way Fourier extrapolation methods, they can be extended for applications to general migration methods.     

Poststack diffraction imaging using reverse‐time migration

We propose a method for imaging small‐scale diffraction objects in complex environments in which Kirchhoff‐based approaches may fail. The proposed method is based on a separation between the specular reflection and diffraction components of the total wavefield in the migrated surface angle domain. Reverse‐time migration was utilized to produce the common image gathers. This approach provides stable and robust results in cases of complex velocity models. The separation is based on the fact that, in surface angle common image gathers, reflection events are focused at positions that correspond to the apparent dip angle of the reflectors, whereas diffracted events are distributed over a wide range of angles. The high‐resolution radon‐based procedure is used to efficiently separate the reflection and diffraction wavefields. In this study, we consider poststack diffraction imaging. The advantages of working in the poststack domain are its numerical efficiency and the reduced computational time. The numerical results show that the proposed method is able to image diffraction objects in complex environments. The application of the method to a real seismic dataset illustrates the capability of the approach to extract diffractions.     

Accelerating the T‐matrix approach to seismic full‐waveform inversion by domain decomposition

A seismic variant of the distorted Born iterative inversion method, which is commonly used in electromagnetic and acoustic (medical) imaging, has been recently developed on the basis of the T‐matrix approach of multiple scattering theory. The distorted Born iterative method is consistent with the Gauss–Newton method, but its implementation is different, and there are potentially significant computational advantages of using the T‐matrix approach in this context. It has been shown that the computational cost associated with the updating of the background medium Green functions after each iteration can be reduced via the use of various linearisation or quasi‐linearisation techniques. However, these techniques for reducing the computational cost may not work well in the presence of strong contrasts. To deal with this, we have now developed a domain decomposition method, which allows one to decompose the seismic velocity model into an arbitrary number of heterogeneous domains that can be treated separately and in parallel. The new domain decomposition method is based on the concept of a scattering‐path matrix, which is well known in solid‐state physics. If the seismic model consists of different domains that are well separated (e.g., different reservoirs within a sedimentary basin), then the scattering‐path matrix formulation can be used to derive approximations that are sufficiently accurate but far more speedy and much less memory demanding because they ignore the interaction between different domains. However, we show here that one can also use the scattering‐path matrix formulation to calculate the overall T‐matrix for a large model exactly without any approximations at a computational cost that is significantly smaller than the cost associated with an exact formal matrix inversion solution. This is because we have derived exact analytical results for the special case of two interacting domains and combined them with Strassen's formulas for fast recursive matrix inversion. To illustrate the fact that we have accelerated the T‐matrix approach to full‐waveform inversion by domain decomposition, we perform a series of numerical experiments based on synthetic data associated with a complex salt model and a simpler two‐dimensional model that can be naturally decomposed into separate upper and lower domains. If the domain decomposition method is combined with an additional layer of multi‐scale regularisation (based on spatial smoothing of the sensitivity matrix and the data residual vector along the receiver line) beyond standard sequential frequency inversion, then one apparently can also obtain stable inversion results in the absence of ultra‐low frequencies and reduced computation times.     

Image‐warping waveform tomography

Imaging the change in physical parameters in the subsurface requires an estimate of the long wavelength components of the same parameters in order to reconstruct the kinematics of the waves propagating in the subsurface. One can reconstruct the model by matching the recorded data with modeled waveforms extrapolated in a trial model of the medium. Alternatively, assuming a trial model, one can obtain a set of images of the reflectors from a number of seismic experiments and match the locations of the imaged interfaces. Apparent displacements between migrated images contain information about the velocity model and can be used for velocity analysis. A number of methods are available to characterize the displacement between images; in this paper, we compare shot‐domain differential semblance (image difference), penalized local correlations, and image‐warping. We show that the image‐warping vector field is a more reliable tool for estimating displacements between migrated images and leads to a more robust velocity analysis procedure. By using image‐warping, we can redefine the differential semblance optimization problem with an objective function that is more robust against cycle‐skipping than the direct image difference. We propose an approach that has straightforward implementation and reduced computational cost compared with the conventional adjoint‐state method calculations. We also discuss the weakness of migration velocity analysis in the migrated‐shot domain in the case of highly refractive media, when the Born modelling operator is far from being unitary and thus its adjoint (migration) operator poorly approximates the inverse.     

Wave-equation migration velocity analysis. I. Theory

We present a migration velocity analysis (MVA) method based on wavefield extrapolation. Similarly to conventional MVA, our method aims at iteratively improving the quality of the migrated image, as measured by the flatness of angle‐domain common‐image gathers (ADCIGs) over the aperture‐angle axis. However, instead of inverting the depth errors measured in ADCIGs using ray‐based tomography, we invert ‘image perturbations’ using a linearized wave‐equation operator. This operator relates perturbations of the migrated image to perturbations of the migration velocity. We use prestack Stolt residual migration to define the image perturbations that maximize the focusing and flatness of ADCIGs. Our linearized operator relates slowness perturbations to image perturbations, based on a truncation of the Born scattering series to the first‐order term. To avoid divergence of the inversion procedure when the velocity perturbations are too large for Born linearization of the wave equation, we do not invert directly the image perturbations obtained by residual migration, but a linearized version of the image perturbations. The linearized image perturbations are computed by a linearized prestack residual migration operator applied to the background image. We use numerical examples to illustrate how the backprojection of the linearized image perturbations, i.e. the gradient of our objective function, is well behaved, even in cases when backprojection of the original image perturbations would mislead the inversion and take it in the wrong direction. We demonstrate with simple synthetic examples that our method converges even when the initial velocity model is far from correct. In a companion paper, we illustrate the full potential of our method for estimating velocity anomalies under complex salt bodies.     

Kinematic characteristics and the influence of reference velocities of phase‐shift‐plus‐interpolation and extended‐split‐step‐Fourier migration methods

The phase‐shift‐plus‐interpolation and extended‐split‐step‐Fourier methods are wavefield‐continuation algorithms for seismic migration imaging. These two methods can be applied to regions with complex geological structures. Based on their unified separable formulas, we show that these two methods have the same kinematic characteristics by using the theory of pseudodifferential operators. Numerical tests on a Marmousi model demonstrate this conclusion. Another important aspect of these two methods is the selection of reference velocities and we explore the influence of the selection of reference velocities by comparing the geometric progression method and the statistical method. We show that the geometric progression method is simple but does not take into account the velocity distribution while the statistical approach is relatively complex but reflects the velocity distribution.     

Stacking angle‐domain common‐image gathers for normalization of illumination

Unequal illumination of the subsurface highly impacts the quality of seismic imaging. Different image points receive different folds of reflection‐angle illumination, which can be caused by irregular acquisition or by wave propagation in complex media. Illumination problems can deteriorate amplitudes in migrated images. To address this problem, we present a method of stacking angle‐domain common‐image gathers, in which we use local similarity with soft thresholding to determine the folds of local illumination. Normalization by local similarity regularizes local illumination of reflection angles for each image point of the subsurface model. This approach compensates for irregular illumination by selective stacking in the image space, regardless of the cause of acquisition or propagation irregularities. Additional migration is not required because the methodology is implemented in the reflection angle domain after migration. We use two synthetic examples to demonstrate that our method can normalize migration amplitudes and effectively suppress migration artefacts.     

基于单程波方程的角度域保幅偏移(英文)

传统叠前深度偏移只能够提供地下的构造信息,但工业界在需要构造信息的同时还要与地下界面反射系数成比例的振幅信息。最近几年,基于单程波方程的保幅叠前深度偏移算法有了一定的发展,但是,基于炮域、单程波的保幅型叠前深度偏移必须应用反褶积型的成像条件,这种成像条件在构造复杂、速度变化剧烈的地区会出现不稳定现象。基于角度域的保幅深度偏移克服了这一不稳定性缺点的同时,还域的保幅深度偏移,模型和实际资料的试算分析验证该思路方法的正确性和有效性。     

Velocity model updating in prestack Kirchhoff time migration for PS converted waves: Part II – Application

We have developed a practical approach for updating the velocity of PS converted waves based on the inverse normal‐moveout common‐image‐point gather obtained from prestack Kirchhoff time migration. We have integrated all the steps involved in updating the migration velocity model into an interactive tool and have applied this approach to a real seismic data set from the Alba Field in the North Sea. Based on experience in handling the real data, we discuss various practical aspects of updating the velocity model, including: what kind of initial velocity model should be used; which parameters in the velocity model should be updated; and how to update them. Application of prestack Kirchhoff time migration to the data set using the updated velocity model produces an improved image of the Alba Field.     

Common‐shot Fresnel beam migration based on wave‐field approximation in effective vicinity under complex topographic conditions

Gaussian beam migration is a versatile imaging method for geologically complex land areas, which overcomes the limitation of Kirchhoff migration in imaging multiple arrivals and has no steep‐dip limits of one‐way wave‐equation migration. However, its imaging accuracy depends on the geometry of Gaussian beam that is determined by the initial parameter of dynamic ray tracing. As a result, its applications in exploration areas with strong variations in topography and near‐surface velocity are limited. Combined with the concept of Fresnel zone and the theory of wave‐field approximation in effective vicinity, we present a more robust common‐shot Fresnel beam imaging method for complex topographic land areas in this paper. Compared with the conventional Gaussian beam migration for irregular topography, our method improves the beam geometry by limiting its effective half‐width with Fresnel zone radius. Moreover, through a quadratic travel‐time correction and an amplitude correction that is based on the wave‐field approximation in effective vicinity, it gives an accurate method for plane‐wave decomposition at complex topography, which produces good imaging results in both shallow and deep zones. Trials of two typical models and its application in field data demonstrated the validity and robustness of our method.     

2D共炮时间域高斯波束偏移

针对传统射线方法在奇异区成像精度不高,而2D频率域高斯波束叠前深度偏移需要计算成像点处每个频率的格林函数,影响计算效率的问题,本文通过使用复走时代替实走时,改变频率域下成像公式的积分顺序,给出了在时间域下进行高斯波束偏移的方法和计算公式.本文使用复杂数值模型验证了2D时间域高斯波束叠前偏移方法的正确性,并同传统射线偏移成像结果做了对比.对比结果表明时间域高斯波束偏移在成像精度上优于传统射线偏移.     

Source–receiver two‐way wave extrapolation for prestack exploding‐reflector modelling and migration

Most modern seismic imaging methods separate input data into parts (shot gathers). We develop a formulation that is able to incorporate all available data at once while numerically propagating the recorded multidimensional wavefield forward or backward in time. This approach has the potential for generating accurate images free of artiefacts associated with conventional approaches. We derive novel high‐order partial differential equations in the source–receiver time domain. The fourth‐order nature of the extrapolation in time leads to four solutions, two of which correspond to the incoming and outgoing P‐waves and reduce to the zero‐offset exploding‐reflector solutions when the source coincides with the receiver. A challenge for implementing two‐way time extrapolation is an essential singularity for horizontally travelling waves. This singularity can be avoided by limiting the range of wavenumbers treated in a spectral‐based extrapolation. Using spectral methods based on the low‐rank approximation of the propagation symbol, we extrapolate only the desired solutions in an accurate and efficient manner with reduced dispersion artiefacts. Applications to synthetic data demonstrate the accuracy of the new prestack modelling and migration approach.     

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.     

Converted‐wave beam migration with sparse sources or receivers

Gaussian beam depth migration overcomes the single‐wavefront limitation of most implementations of Kirchhoff migration and provides a cost‐effective alternative to full‐wavefield imaging methods such as reverse‐time migration. Common‐offset beam migration was originally derived to exploit symmetries available in marine towed‐streamer acquisition. However, sparse acquisition geometries, such as cross‐spread and ocean bottom, do not easily accommodate requirements for common‐offset, common‐azimuth (or common‐offset‐vector) migration. Seismic data interpolation or regularization can be used to mitigate this problem by forming well‐populated common‐offset‐vector volumes. This procedure is computationally intensive and can, in the case of converted‐wave imaging with sparse receivers, compromise the final image resolution. As an alternative, we introduce a common‐shot (or common‐receiver) beam migration implementation, which allows migration of datasets rich in azimuth, without any regularization pre‐processing required. Using analytic, synthetic, and field data examples, we demonstrate that converted‐wave imaging of ocean‐bottom‐node data benefits from this formulation, particularly in the shallow subsurface where regularization for common‐offset‐vector migration is both necessary and difficult.     

转换波Kirchhoff叠前时间偏移的成像优化方案

叠前时间偏移是转换波数据处理流程中的重要技术环节.由于下行P波和上行S波对应着不同类型的传播速度和传播机理,使得转换波偏移的诸多实施环节和实现技巧均有别于常规的纵波偏移.本文就VTI介质转换波Kirchhoff叠前时间偏移实现过程中的几个关键技术环节进行研究分析,并给出了相应的实施方案.首先,依据常速假设推导了简化的转换波时间偏移振幅加权函数,在保证偏移的高效计算的同时,兼顾了成像振幅的准确性;其次,依据直射线近似假设,推导了转换波偏移的最大无假频频率计算公式,确保了算子假频的有效滤除;接下来,设计使用了依据共转换点对称的偏移孔径,保证了偏移孔径的优化选取.此外,本文还提出利用三次卷积插值进行成像网格上参数值的插值运算,以避免传统线性插值算法的非平滑性对偏移剖面中高频信息的不良影响.多个算例的应用效果证明,本文所给出的关键技术环节的实施方案可以有效的提升转换波叠前时间偏移的成像效果.     

Suppressing non‐Gaussian noises with scaled receiver wavefield for reverse‐time migration: comparison of different approaches

Numerical implementation of the gradient of the cost function in a gradient‐based full‐ waveform inversion (FWI) is essentially a migration operator used in wave equation migration. In FWI, minimizing different data residual norms results in different weighting strategies of data residuals at receiver locations prior to back‐propagation into the medium. In this paper, we propose different scaling methods to the receiver wavefield and compare their performances. Using time‐domain reverse‐time migration (RTM), we show that compared to conventional algorithms, this type of scaling is able to significantly suppress non‐Gaussian noise, i.e., outliers. Our tests also show that scaling by its absolute norm produces better results than other approaches.     

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 .     

Multipathing in three‐parameter common‐image gathers from reverse‐time migration

Reverse‐time migration gives high‐quality, complete images by using full‐wave extrapolations. It is thus not subject to important limitations of other migrations that are based on high‐frequency or one‐way approximations. The cross‐correlation imaging condition in two‐dimensional pre‐stack reverse‐time migration of common‐source data explicitly sums the product of the (forward‐propagating) source and (backward‐propagating) receiver wavefields over all image times. The primary contribution at any image point travels a minimum‐time path that has only one (specular) reflection, and it usually corresponds to a local maximum amplitude. All other contributions at the same image point are various types of multipaths, including prismatic multi‐arrivals, free‐surface and internal multiples, converted waves, and all crosstalk noise, which are imaged at later times, and potentially create migration artefacts. A solution that facilitates inclusion of correctly imaged, non‐primary arrivals and removal of the related artefacts, is to save the depth versus incident angle slice at each image time (rather than automatically summing them). This results in a three‐parameter (incident angle, depth, and image time) common‐image volume that integrates, into a single unified representation, attributes that were previously computed by separate processes. The volume can be post‐processed by selecting any desired combination of primary and/or multipath data before stacking over image time. Separate images (with or without artifacts) and various projections can then be produced without having to remigrate the data, providing an efficient tool for optimization of migration images. A numerical example for a simple model shows how primary and prismatic multipath contributions merge into a single incident angle versus image time trajectory. A second example, using synthetic data from the Sigsbee2 model, shows that the contributions to subsalt images of primary and multipath (in this case, turning wave) reflections are different. The primary reflections contain most of the information in regions away from the salt, but both primary and multipath data contribute in the subsalt region.