THE NON-LINEAR INVERSION OF SEISMIC WAVEFORMS CAN BE PERFORMED EITHER BY TIME EXTRAPOLATION OR BY DEPTH EXTRAPOLATION1

The classical aim of non-linear inversion of seismograms is to obtain the earth model which, for null initial conditions and given sources, best predicts the observed seismograms. This problem is currently solved by an iterative method: each iteration involves the resolution of the wave equation with the actual sources in the current medium, the resolution of the wave equation, backwards in time, with the current residuals as sources; and the correlation, at each point of space, of the two wavefields thus obtained. Our view of inversion is more general: we want to obtain a whole set of earth model, initial conditions, source functions, and predicted seismograms, which are the closest to some a priori values, and which are related through the wave equation. It allows us to justify the previous method, but it also allows us to set the same inverse problem in a different way: what is now searched for is the best fit between calculated and a priori initial conditions, for given sources and observed surface displacements. This leads to a completely different iterative method, in which each iteration involves the downward extrapolation of given surface displacements and tractions, down to a given depth (the‘bottom’), the upward extrapolation of null displacements and tractions at the bottom, using as sources the initial time conditions of the previous field, and a correlation, at each point of the space, of the two wavefields thus obtained. Besides the theoretical interest of the result, it opens the way to alternative numerical methods of resolution of the inverse problem. If the non-linear inversion using forward-backward time propagations now works, this non-linear inversion using downward-upward extrapolations will give the same results but more economically, because of some tricks which may be used in depth extrapolation (calculation frequency by frequency, inversion of the top layers before the bottom layers, etc.).     

Multi‐block finite‐difference method for 3D elastodynamic simulations in anisotropic subhorizontally layered media

Prediction of elastic full wavefields is required for reverse time migration, full waveform inversion, borehole seismology, seismic modelling, etc. We propose a novel algorithm to solve the Navier wave equation, which is based on multi‐block methodology for high‐order finite‐difference schemes on curvilinear grids. In the current implementation, the blocks are subhorizontal layers. Smooth anisotropic heterogeneous media in each layer can have strong discontinuities at the interfaces. A curvilinear adaptive hexahedral grid in blocks is generated by mapping the original 3D physical domain onto a parametric cube with horizontal layers and interfaces. These interfaces correspond to the main curvilinear physical contrast interfaces of a subhorizontally layered formation. The top boundary of the parametric cube handles the land surface with smooth topography. Free‐surface and solid–solid transmission boundary conditions at interfaces are approximated with the second‐order accuracy. Smooth media in the layers are approximated up to sixth‐order spatial schemes. All expected properties of the developed algorithm are demonstrated in numerical tests using corresponding parallel message passing interface code.     

NON-LINEAR AVO INVERSION FOR A STACK OF ANELASTIC LAYERS1

Parameters in a stack of homogeneous anelastic layers are estimated from seismic data, using the amplitude versus offset (AVO) variations and the travel-times. The unknown parameters in each layer are the layer thickness, the P-wave velocity, the S-wave velocity, the density and the quality factor. Dynamic ray tracing is used to solve the forward problem. Multiple reflections are included, but wave-mode conversions are not considered. The S-wave velocities are estimated from the PP reflection and transmission coefficients. The inverse problem is solved using a stabilized least-squares procedure. The Gauss-Newton approximation to the Hessian matrix is used, and the derivatives of the dynamic ray-tracing equation are calculated analytically for each iteration. A conventional velocity analysis, the common mid-point (CMP) stack and a set of CMP gathers are used to identify the number of layers and to establish initial estimates for the P-wave velocities and the layer thicknesses. The inversion is carried out globally for all parameters simultaneously or by a stepwise approach where a smaller number of parameters is considered in each step. We discuss several practical problems related to inversion of real data. The performance of the algorithm is tested on one synthetic and two real data sets. For the real data inversion, we explained up to 90% of the energy in the data. However, the reliability of the parameter estimates must at this stage be considered as uncertain.     

Separation of P- and S-wavefields from wide-angle multicomponent OBC data for a basalt model

Wide-angle multicomponent ocean-bottom cable (OBC) data should further enhance sub-basalt imaging by using both compressional and converted shear wavefields. The first step in analysing multicomponent OBC data is to decompose the recorded wavefields into pure P- and pure S-wavefields, and extract the upgoing P- and S-waves. This paper presents a new scheme to separate P- and S-wavefields from wide-angle multicomponent OBC data in the τ–p domain. By considering plane-wave components with a known horizontal slowness, the P- and S-wavefields are separated into the directions of observed P- and S-wave oscillations using the horizontal and vertical components of the data. The upgoing P- and S-waves are then extracted from the separated P- and S-wavefields. The parameters used in the separation are the seismic wave velocities and the density at the receiver location, which can be estimated from the first reflection phase observed on the horizontal and vertical components. Numerical tests on synthetic data for a plane-layered model show good performance and demonstrate the accuracy of the scheme. Separation of wavefields from a basalt model is performed using synthetic wide-angle multicomponent OBC data. The results show that both near-offset and wide-angle reflections and conversions from within and below basalt layers are enhanced and clearly identified on the separated wavefields.     

RESISTIVITY SOUNDING ON A HORIZONTALLY STRATIFIED MULTI-LAYERED EARTH*

The electrical potential due to a point source of current placed on the ground surface is studied for a multi-layered earth consisting of homogeneous overburden of constant conductivity over a stack of transition layers where conductivity varies with depth according to power and exponential laws in even and odd layers, respectively. The general recursion relations are derived and are used to obtain expressions for the apparent resistivities for Schlumberger and Wenner electrode arrays. Their asymptotic behavior has been studied. The solutions for some particular cases are given: (i) odd layers with conductivities exponentially varying with depth while all even layers (and the first) have constant conductivities; (ii) even layer conductivities varying as a power law while odd number layers are of constant conductivity; and (iii) any two successive layers as transition layers and all others having constant conductivities. Further it is shown that Patella's theory is a particular case of the present study. It is concluded, therefore, that the present treatment is more general as all earlier models consisting of trasition layers can be derived from this study.     

Reflected waves in finely layered tilted orthorhombic media

Upscaling in seismics is a homogenization of finely layered media in the zero-frequency limit. An upscaling technique for arbitrary anisotropic layers has been developed by Schoenberg and Muir. Applying this technique to a stack of layers of orthorhombic (ORT) symmetry whose vertical symmetry planes are aligned, results in an effective homogeneous layer with orthorhombic symmetry. If the symmetry planes in a horizontal orthorhombic layer are rotated with respect to vertical, the medium is referred to as tilted orthorhombic (TOR) medium, and the stack composed of TOR layers in zero-frequency limit will produce an effective medium of a lower symmetry than orthorhombic. We consider a P-wave that propagates through a stack of thin TOR layers, then it is reflected (preserving the mode) at some interface below the stack, and then propagates back through the same stack. We propose to use a special modified medium for the upscaling in case of this sequential down- and up-propagation: each TOR layer in the stack is replaced by two identical TOR layers whose tilt angles have the opposite algebraic sign. In this modified medium, one-way propagation of a seismic wave (any wave mode) is equivalent to propagation of a pure-mode reflection in the original medium. We apply this idea to study the contribution from an individual layer from the stack and show how the approach can be applied to a stack of TOR layers. To demonstrate the applicability of the model, we use well log data for the upscaling. The model we propose for the upscaling can be used in well-seismic ties to correct the effective parameters obtained from well log data for the presence of tilt, if latter is confirmed by additional measurements (for example, borehole imaging).     

Time evolution of the first‐order linear acoustic/elastic wave equation using Lie product formula and Taylor expansion

We propose a new numerical solution to the first‐order linear acoustic/elastic wave equation. This numerical solution is based on the analytic solution of the linear acoustic/elastic wave equation and uses the Lie product formula, where the time evolution operator of the analytic solution is written as a product of exponential matrices where each exponential matrix term is then approximated by Taylor series expansion. Initially, we check the proposed approach numerically and then demonstrate that it is more accurate to apply a Taylor expansion for the exponential function identity rather than the exponential function itself. The numerical solution formulated employs a recursive procedure and also incorporates the split perfectly matched layer boundary condition. Thus, our scheme can be used to extrapolate wavefields in a stable manner with even larger time‐steps than traditional finite‐difference schemes. This new numerical solution is examined through the comparison of the solution of full acoustic wave equation using the Chebyshev expansion approach for the matrix exponential term. Moreover, to demonstrate the efficiency and applicability of our proposed solution, seismic modelling results of three geological models are presented and the processing time for each model is compared with the computing time taking by the Chebyshev expansion method. We also present the result of seismic modelling using the scheme based in Lie product formula and Taylor series expansion for the first‐order linear elastic wave equation in vertical transversely isotropic and tilted transversely isotropic media as well. Finally, a post‐stack migration results are also shown using the proposed method.     

ON GETTING REFLECTION COEFFICIENTS FROM WAVES *

A horizontally layered non-absorptive system of homogeneous layers may be specified by giving the reflection coefficients at each interface. Provided the layers have equal vertical travel time and a perfect reflector as a free surface, the reflection coefficients are generally reconstructed from the reflected pulses by way of solving simultaneous equations of the Toeplitz matrix form with the Levinson recursion method. There exists an alternative approach to solving this problem which by simple reasoning immediately turns out the (Levinson) recursion scheme. The method is based on formulas that relate to solving the forward problem. It resembles Kunetz's (1962) original inverse solution in as much as the computation of the reflection coefficients is based on the idea of separating the contribution of a primary from the sum of all multiples.     

A direct method of calculating interval velocities and layer thicknesses from wide-angle seismic reflection times

Wide-angle reflections are now routinely recorded in high resolution explosion seismics to study the crustal structure. Use of Dix's hyperbolic approximation to the nonhyperbolic wide-angle reflection travel times causes major errors in the determination of interval velocities and layer thicknesses of a stack of horizontal velocity layers. Here we propose a layer stripping method to directly calculate the interval velocities and layer thicknesses in a vertically heterogeneous earth from the strong and reliable wide-angle reflected events. Synthetic reflection travel times, at wide-angle range, for a given velocity model, contaminated by some random errors, have been used to demonstrate the reliability of the algorithms to determine the interval velocities and thicknesses of various layers. The method has also been tested on two field examples along two deep seismic sounding (DSS) profiles with well identified wide-angle reflection travel times, which illustrates the practical feasibility of the proposed method.     

Perfectly matched layer boundary conditions for the second-order acoustic wave equation solved by the rapid expansion method

We derive a governing second-order acoustic wave equation in the time domain with a perfectly matched layer absorbing boundary condition for general inhomogeneous media. Besides, a new scheme to solve the perfectly matched layer equation for absorbing reflections from the model boundaries based on the rapid expansion method is proposed. The suggested scheme can be easily applied to a wide class of wave equations and numerical methods for seismic modelling. The absorbing boundary condition method is formulated based on the split perfectly matched layer method and we employ the rapid expansion method to solve the derived new perfectly matched layer equation. The use of the rapid expansion method allows us to extrapolate wavefields with a time step larger than the ones commonly used by traditional finite-difference schemes in a stable way and free of dispersion noise. Furthermore, in order to demonstrate the efficiency and applicability of the proposed perfectly matched layer scheme, numerical modelling examples are also presented. The numerical results obtained with the put forward perfectly matched layer scheme are compared with results from traditional attenuation absorbing boundary conditions and enlarged models as well. The analysis of the numerical results indicates that the proposed perfectly matched layer scheme is significantly effective and more efficient in absorbing spurious reflections from the model boundaries.     

INDUCED POLARIZATION RESPONSE OF A HORIZONTALLY MULTILAYERED EARTH WITH NO RESISTIVITY CONTRAST*

The induced polarization response of a horizontally multilayered earth with no resistivity contrast can rapidly be calculated on a desk calculator or minicomputer for any electrode array. The formulation is a simple series summation of the products of weighting coefficients and the true induced polarization responses for each of the layers. The coefficients are directly derivable from the corresponding resistivity model. This series approach to IP formulation was originally described by Seigel but has not been treated extensively in the present-day geophysical literature. This method can be applied to either time or frequency domain induced polarization measurements. Once the coefficients are known, apparent induced polarization response can readily be obtained by judicious substitution of known, suspected, or assumed values of the true induced polarization of each layer. Basic formulation is presented for the IP potential coefficients (pole-pole or two array) with no resistivity contrast between the layers. From these coefficients, response of any number of layers for any electrode array can be obtained by suitable differentiation. Some examples of Wenner array for a three-layered earth and dipole-dipole array for a four-layered earth are used to illustrate the application. The results of this technique are valid for many natural situations of modest resistivity contrast. However, they definitely cannot be used if there are highly contrasting resistivity layers present. Such an approach is conceptually simple and is useful for survey planning, checking or setting the “depth-of-penetration”of a given array. For field induced polarization data that fits reasonably well to the no-resistivity-contrast model, this simple approach facilitates quantitative interpretation.     

Mapping Coastal Aquifers by Joint Inversion of DC and TEM Soundings-Three Case Histories

Electrical and electromagnetic methods are well suited for coastal aquifer studies because of the large contrast in resistivity between fresh water-bearing and salt water-bearing formations. Interpretation models for these aquifers typically contain four layers: a highly resistive unsaturated zone; a surficial fresh water aquifer of intermediate resistivity; an underlying conductive, salt water saturated aquifer; and resistive substratum. Additional layers may be added to allow for variations in lithology within the fresh water and salt water layers. Two methods are evaluated: direct current resistivity and time domain electromagnetic soundings. Use of each method alone produces nonunique solutions for resistivities and/or thicknesses of the different layers. We show that joint inversion of vertical electric and time domain electromagnetic soundings produces a more tightly constrained interpretation model at three test sites than is produced by inversion methods applied to each data set independently.     

Seismic Wavefields in the Deep Seafloor Area from a Submarine Landslide Source

We use the finite difference method to simulate seismic wavefields at broadband land and seafloor stations for a given terrestrial landslide source, where the seafloor stations are located at water depths of 1,900–4,300 m. Our simulation results for the landslide source explain observations well at the seafloor stations for a frequency range of 0.05–0.1 Hz. Assuming the epicenter to be located in the vicinity of a large submarine slump, we also model wavefields at the stations for a submarine landslide source. We detect propagation of the Airy phase with an apparent velocity of 0.7 km/s in association with the seawater layer and an accretionary prism for the vertical component of waveforms at the seafloor stations. This later phase is not detected when the structural model does not consider seawater. For the model incorporating the seawater, the amplitude of the vertical component at seafloor stations can be up to four times that for the model that excludes seawater; we attribute this to the effects of the seawater layer on the wavefields. We also find that the amplification of the waveform depends not only on the presence of the seawater layer but also on the thickness of the accretionary prism, indicating low amplitudes at the land stations and at seafloor stations located near the trough but high amplitudes at other stations, particularly those located above the thick prism off the trough. Ignoring these characteristic structures in the oceanic area and simply calculating the wavefields using the same structural model used for land areas would result in erroneous estimates of the size of the submarine landslide and the mechanisms underlying its generation. Our results highlight the importance of adopting a structural model that incorporates the 3D accretionary prism and seawater layer into the simulation in order to precisely evaluate seismic wavefields in seafloor areas.     

Vertical Discretization Impact in Numerical Modeling of Matrix Diffusion in Contaminated Groundwater

Understanding the effects of contaminants that can diffuse into low-permeability (“low-k”) zones is crucial for effective groundwater remedial decision-making. Because low-k zones can serve as low-level sources of contamination to more transmissive zones over time, an accurate evaluation of the impacts of matrix diffusion at contaminated sites is vital. This study compared numerical groundwater flow and transport simulations using MODFLOW/RT3D at a hypothetical site using three cases, each with increasing discretization of the vertical 10-m thick domain: (1) a coarse multilayer heterogeneous grid based on one layer for each of four different hydrogeological units, (2) a “low-resolution” discretization approach where the low-k units were divided into several sublayers giving the model 10 layers, and (3) a “high-resolution” numerical model with 199 layers that are a few centimeters thick. When comparing the results of each case, significant differences were observed between the discretizations used, even though all other model input data were identical. The conventional grid models (Cases 1 and 2) appeared to underestimate groundwater plume concentrations by a factor ranging from 1.1 to 36 when compared to the high-resolution grid model (Case 3), and underestimated predicted cleanup times by more than a factor of 10 for some of the hypothetical sampling points in the modeling domain. These results validate the implication of Chapman et al. (2012), that conventional vertical discretization of numerical groundwater flow and transport models at contaminated sites (with layers that are greater than 1 m thick) can lead to significant errors when compared to more accurate high-resolution vertical discretization schemes (layers that are centimeters thick).     

Application of a perfectly matched layer in seismic wavefield simulation with an irregular free surface

Recently, an effective and powerful approach for simulating seismic wave propagation in elastic media with an irregular free surface was proposed. However, in previous studies, researchers used the periodic condition and/or sponge boundary condition to attenuate artificial reflections at boundaries of a computational domain. As demonstrated in many literatures, either the periodic condition or sponge boundary condition is simple but much less effective than the well‐known perfectly matched layer boundary condition. In view of this, we intend to introduce a perfectly matched layer to simulate seismic wavefields in unbounded models with an irregular free surface. We first incorporate a perfectly matched layer into wave equations formulated in a frequency domain in Cartesian coordinates. We then transform them back into a time domain through inverse Fourier transformation. Afterwards, we use a boundary‐conforming grid and map a rectangular grid onto a curved one, which allows us to transform the equations and free surface boundary conditions from Cartesian coordinates to curvilinear coordinates. As numerical examples show, if free surface boundary conditions are imposed at the top border of a model, then it should also be incorporated into the perfectly matched layer imposed at the top‐left and top‐ right corners of a 2D model where the free surface boundary conditions and perfectly matched layer encounter; otherwise, reflections will occur at the intersections of the free surface and the perfectly matched layer, which is confirmed in this paper. So, by replacing normal second derivatives in wave equations in curvilinear coordinates with free surface boundary conditions, we successfully implement the free surface boundary conditions into the perfectly matched layer at the top‐left and top‐right corners of a 2D model at the surface. A number of numerical examples show that the perfectly matched layer constructed in this study is effective in simulating wave propagation in unbounded media and the algorithm for implementation of the perfectly matched layer and free surface boundary conditions is stable for long‐time wavefield simulation on models with an irregular free surface.     

Effect of anelastic patchy saturated sand layers on the reflection and transmission responses of a periodically layered medium

Based on analytic relations, we compute the reflection and transmission responses of a periodically layered medium with a stack of elastic shales and partially saturated sands. The sand layers are considered anelastic (using patchy saturation theory) or elastic (with effective velocity). Using the patchy saturation theory, we introduce a velocity dispersion due to mesoscale attenuation in the sand layer. This intrinsic anelasticity is creating frequency dependence, which is added to the one coming from the layering (macroscale). We choose several configurations of the periodically layered medium to enhance more or less the effect of anelasticity. The worst case to see the effect of intrinsic anelasticity is obtained with low dispersion in the sand layer, strong contrast between shales and sands, and a low value of the net‐to‐gross ratio (sand proportion divided by the sand + shale proportion), whereas the best case is constituted by high dispersion, weak contrast, and high net‐to‐gross ratio. We then compare the results to show which dispersion effect is dominating in reflection and transmission responses. In frequency domain, the influence of the intrinsic anelasticity is not negligible compared with the layering effect. Even if the main resonance patterns are the same, the resonance peaks for anelastic cases are shifted towards high frequencies and have a slightly lower amplitude than for elastic cases. These observations are more emphasized when we combine all effects and when the net‐to‐gross ratio increases, whereas the differences between anelastic and elastic results are less affected by the level of intrinsic dispersion and by the contrast between the layers. In the time domain, the amplitude of the responses is significantly lower when we consider intrinsic anelastic layers. Even if the phase response has the same features for elastic and anelastic cases, the anelastic model responses are clearly more attenuated than the elastic ones. We conclude that the frequency dependence due to the layering is not always dominating the responses. The frequency dependence coming from intrinsic visco‐elastic phenomena affects the amplitude of the responses in the frequency and time domains. Considering intrinsic attenuation and velocity dispersion of some layers should be analyzed while looking at seismic and log data in thin layered reservoirs.     

Improved amplitude preservation for prestack depth migration by inverse scattering theory

A prestack reverse time-migration image is not properly scaled with increasing depth. The main reason for the image being unscaled is the geometric spreading of the wavefield arising during the back-propagation of the measured data and the generation of the forward-modelled wavefields. This unscaled image can be enhanced by multiplying the inverse of the approximate Hessian appearing in the Gauss–Newton optimization technique. However, since the approximate Hessian is usually too expensive to compute for the general geological model, it can be used only for the simple background velocity model.We show that the pseudo-Hessian matrix can be used as a substitute for the approximate Hessian to enhance the faint images appearing at a later time in the 2D prestack reverse time-migration sections. We can construct the pseudo-Hessian matrix using the forward-modelled wavefields (which are used as virtual sources in the reverse time migration), by exploiting the uncorrelated structure of the forward-modelled wavefields and the impulse response function for the estimated diagonal of the approximate Hessian. Although it is also impossible to calculate directly the inverse of the pseudo-Hessian, when using the reciprocal of the pseudo-Hessian we can easily obtain the inverse of the pseudo-Hessian. As examples supporting our assertion, we present the results obtained by applying our method to 2D synthetic and real data collected on the Korean continental shelf.     

伪随机编码源激发下的时域电磁信号合成

将伪随机编码技术引入到人工源电磁法后,可以通过加大发射功率以及应用后续的相关处理技术来达到压制噪声、加大探测深度及提高分辨率的目的,因此引起了越来越多学者的关注及研究,但大多数研究集中在资料处理的相关技术上,对模拟电磁信号关注较少.然而,资料处理工作大多是从电磁信号出发的,模拟伪随机编码源激发下的电磁信号不但可以为资料处理环节提供理论数据,而且可以为检测资料处理的效果提供中间结果,因此,模拟电磁信号工作必不可少.本文根据获得接收信号的物理过程来实现伪随机编码源激发下的电磁信号合成.首先用解析公式获得特定地电结构的大地频率域响应,然后通过余弦变换得到时间域阶跃响应,接下来用阶跃响应的时间导数得到大地脉冲响应,通过将大地脉冲响应与伪随机编码源的褶积得到理想接收信号,最后,用低通滤波器来模拟发射设备和接收设备的频带限制,将之和噪声一起加到理想接收信号上,最终模拟出仿真的合成信号.通过和野外实际接收信号对比发现本文合成信号仿真度较高,可以服务于后续的数据处理环节.     

编码电磁测深

利用逆重复M序列伪随机信号良好的自相关特性,提出相关检测地电传输特性的编码电磁勘探法.分析了编码电磁测深原理及参数提取方法.同时记录发送电流信号和多收发距电磁场响应,对源信号和场信号进行相关运算,解卷积分离接收系统响应后可得到大地的频率特性或时间特性,以此实现地电断面的精细探测.在时间域,大地冲激响应和阶跃响应含有丰富的地电信息.通过冲激响应的峰值时间或阶跃响应的晚期渐近值可估计地电阻率分布.基于层状模型的大地冲激响应和阶跃响应正演计算结果表明,编码电磁测深法对大埋深薄层目标体有精细的分辨能力.可应用于能源、矿产资源、水资源、环境地质及工程地质勘察,有良好的应用前景.     

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.