The Inverse-Scattering Imaging Condition (ISIC) is an imaging condition for Reverse-Time Migration (RTM) that attempts to recover the medium reflectivity. It is theoretically based on the asymptotic inverse to the Born approximation and can be represented in several theoretically approximately equivalent forms. Its application leads to more reliable reflectivity estimates and strongly reduces backscattering artifacts. In this work, we demonstrate that an ISIC formulation involving a Laplacian filter can be used as an effective preconditioning for Least-Squares RTM (LSRTM). The Laplacian-filter ISIC does not increase the computational cost over conventional imaging conditions. Our numerical experiments using synthetic seismic data from the Marmousi II model demonstrate that this preconditioning leads to faster convergence and superior final images of LSRTM, both in the image domain (ID-LSRTM) and data domain (DD-LSRTM), in this way actually reducing computational cost and turnaround time.

Full waveform inversion (FWI), a powerful geophysical technique for subsurface imaging through seismic velocity-model construction, relies on numerical optimization, thus requiring the computation of derivatives for an objective function. This paper proposes a discrete development for accurate computation of the gradient and Hessian-vector product, providing second-order optimization benefits like higher convergence rates and improved resolution. The approach is a promising alternative for computing the gradient and Hessian action in time-domain FWI, applicable to various geophysical problems. Computational costs and memory requirements are comparable to the Adjoint-State Method and more avorable than Automatic Differentiation. While efficient automatic differentiation algorithms have transformed gradient computation in applications like FWI, challenges may arise in 3D due to unforeseen memory allocations. Our approach addresses this by exploring the reverse mode differentiation algorithm, mapping temporary memory allocations and computational complexity. By means of introducing auxiliary fields all involved wavefield evolutions can be carried out with the very same evolution scheme, in this way simplifying the implementation and focusing the performance improvement effort in a single routine thus reducing the maintenance cost of these algorithms, especially when using GPU implementations.

Seismic inversion methods are crucial for understanding subsurface structures; however, the presence of noise in the data can negatively impact the results of these methods. To address this issue, state-of-art applications employ a variety of techniques to enhance the signal over the noise, with researchers developing increasingly sophisticated denoising methods and combining them into new procedures to further improve the signal-to-noise ratio (SNR). While some methodologies operate on a single scale, the curvelet transform is a multi-scale transform useful for decomposing seismic signals into multi-resolution and multi-direction elements. Here, we evaluate the effectiveness of denoising by means of curvelet thresholding as a preconditioning method for poststack seismic data in a 2D acoustic inversion processing using a Bayesian framework. Our application of the curvelet thresholding method on the Marmousi model and a real dataset from a Brazilian offshore basin demonstrates that it can successfully eliminate random noise. Even the use of a hard global threshold, dependent on curvelet scale and orientation, led to improvements in the deepest parts of the models. However, we observed a decrease in the SNR in the presence of soft rocks with pronounced absorption, which are typical in the shallowest regions. Future work will have to explore alternative methods for selecting coefficients that can robustly incorporate changes in the seismic wavelet with depth.

The inverse-scattering imaging condition (ISIC) for reverse time migration (RTM) aims at recovering amplitudes proportional to seismic reflectivity. It has been derived as the high-frequency asymptotic inverse of Born modeling, which justifies its being called a true-amplitude imaging condition. It involves the temporal and spatial derivatives of the up- and downgoing wavefields, in this way generalizing the conventional crosscorrelation imaging condition. The temporal derivations can be redistributed between different wavefield contributions, in this way deriving a set of different implementational forms of the ISIC. By making use of the wave equation for the up- and downgoing wavefields, one can substitute the time derivatives by the Laplacian operator. This provides a theoretical foundation for a popular filter for reducing the backscattering artifacts in RTM. Using Born data from a simple three-layer model and the Marmousi II model as well as the Sigsbee2b data, we have determined that the theoretical equivalence of the equations leads to similar but not identical images. Our numerical tests indicate that the ISIC versions using spatial derivatives are the most economical approach, and that the images obtained with the second time derivative of the source wavefield indicate slightly improved resolution over the other implementations, making the combination of these two characteristics the best choice.

The numerical simulation of wave propagation can be represented by a propagator matrix applied to previous instances of the wavefield. Using the sparsity of the propagator matrix to approximate it by a low-rank representation, one can increment the wavefield’s phase from one time instance to another. Though this procedure does not pay attention to the amplitudes of the seismic waves, it is important to understand its dynamic properties. Here, we evaluate the amplitudes obtained by the low-rank method in the simulation of 2D acoustic wave propagation. In homogeneous media, where theoretical expressions for the wavefield are available, the method provides not only an excellent kinematic approximation, but also reliable amplitudes. For a single horizontal reflector below a homogeneous overburden, the reflection coefficients approximated by the low-rank method are of the same quality or slightly superior to those obtained by a second-order finite-difference (FD) method (implementation from SU). However, in more generally inhomogeneous media, our tests showed larger discrepancies between FD and low-rank modeling results. While comparing unfavorably with FD regarding computation time for small models, its quasi-linear scaling with model size makes the low-rank method superior for large models. Moreover, a generalization to more complex, e.g., anisotropic, media is straightforward.

The inverse scattering imaging condition (ISIC) for reverse-time migration (RTM) aims at recovering amplitudes proportional to seismic reflectivity. It has been derived as the high-frequency asymptotic inverse of Born modeling, which justifies it being called a true-amplitude imaging condition. It involves the temporal and spatial derivatives of the up- and downgoing wavefields, in this way generalizing the conventional crosscorrelation imaging condition. The temporal derivations can be redistributed between the different wavefield contributions, in this way coming up with a set of different implementational forms of the ISIC. By making use of the wave equation for the up- and downgoing wavefields, one can substitute the time derivatives by the Laplacian operator. This provides a theoretical foundation for a popular filter for reducing backscattering artifacts in RTM. Using data from a simple three-layer model as well as the Marmousi II and Sigsbee2A data, we demonstrate that the theoretical equivalence of the equations leads to similar, but not identical, images. Our numerical tests indicate that the ISIC using spatial derivatives are the most economic ones, and that the images obtained with the second time derivative of the source wavefield show slightly improved resolution over the other implementations, making the combination of these two characteristics the best choice.

We have developed a procedure to derive low-rank evolution operators in the mixed space-wavenumber domain for modeling the qP Born-scattered wavefield at perturbations of an anisotropic medium under the pseudoacoustic approximation. To approximate the full wavefield, this scattered field is then added to the reference wavefield obtained with the corresponding low-rank evolution operator in the background medium. Being built upon a Hamiltonian formulation using the dispersion relation for qP-waves, this procedure avoids pseudo-S-wave artifacts and provides a unified approach for linearizing anisotropic pseudoacoustic evolution operators. Therefore, it is immediately applicable to any arbitrary class of anisotropy. As an additional asset, the scattering operators explicitly contain the sensitivity kernels of the Born-scattered wavefield with respect to the anisotropic medium parameters. This enables direct access to important information such as its offset dependence or directional characteristics as a function of the individual parameter perturbations. For our numerical tests, we specify the operators for a mildly anisotropic tilted transversely isotropic (TTI) medium. We validate our implementation in a simple model with weak contrasts and simulate reflection data in the BP TTI model to indicate that the procedure works in a more realistic scenario. The Born-scattering results indicate that our procedure is applicable to strongly heterogeneous anisotropic media. Moreover, we use the analytical capabilities of the kernels by means of sensitivity tests to demonstrate that using two different medium parameterizations leads to different results. The mathematical formulation of the method is such that it allows for an immediate application to least-squares migration in pseudoacoustic anisotropic media.

Redatuming is used to relocate sources and receivers to a new depth level.This aim can be achieved by different techniques, for example the interferometric method or the wavefield continuation.In this work, we propose to use redatuming by backpropagating the seismic records.To apply this methodology, we need to know the medium above the redatuming level.To validate the theory discussed, we used 2D synthetic examples with different complexities in the overburden, and we compare the results to those of correlation-based redatuming.It was possible to note that the backpropagation-based redatuming presented similarities to the correlation-based one, regarding the events positioning and in the addition of artifacts.We also note that both techniques presented almost the same sensibility to an inexact velocity model.Despite the similarities, the method based on the wavefield continuation required significantly less computation time when we compare it to the interferometric technique.

Joint migration inversion (JMI) is a method based on one-way wave equations that aims at fitting seismic reflection data to estimate an image and a background velocity. The depth-migrated image describes the high spatial-frequency content of the subsurface and, in principle, is true amplitude. The background velocity model accounts mainly for the large spatial-scale kinematic effects of the wave propagation. Looking for a deeper understanding of the method, we briefly review the continuous equations that compose the forward-modeling engine of JMI for acoustic media and angle-independent scattering. Then, we use these equations together with the first-order adjoint-state method to arrive at a new formulation of the model gradients. To estimate the image, we combine the second-order adjoint-state method with the truncated-Newton method to obtain the image updates. For the model (velocity) estimation, in comparison to the image update, we reduce the computational cost by adopting a diagonal preconditioner for the corresponding gradient in combination with an image-based regularizing function. Based on this formulation, we build our implementation of the JMI algorithm. Our image-based regularization of the model estimate allows us to carry over structural information from the estimated image to the jointly estimated background model. As demonstrated by our numerical experiments, this procedure can help to improve the resolution of the estimated model and make it more consistent with the image.

Eikonal solvers have important applications in seismic data processing and inversion, the so-called image-guided methods. To this day, in image-guided applications, the solution of the eikonal equation is implemented using partial-differential-equation solvers, such as fast-marching or fast-sweeping methods. We have found that alternatively, one can numerically integrate the dynamic Hamiltonian system defined by the image-guided eikonal equation and reconstruct the solution with image-guided rays. We evaluate interesting applications of image-guided ray tracing to seismic data processing, demonstrating the use of the resulting rays in image-guided interpolation and smoothing, well-log interpolation, image flattening, and residual-moveout picking. Some of these applications make use of properties of the ray-tracing system that are not directly obtained by eikonal solvers, such as ray position, ray density, wavefront curvature, and ray curvature. These ray properties open space for a different set of applications of the image-guided eikonal equation, beyond the original motivation of accelerating the construction of minimum distance tables. We stress that image-guided ray tracing is an embarrassingly parallel problem that makes its implementation highly efficient on massively parallel platforms. Image-guided ray tracing is advantageous for most applications involving the tracking of seismic events and imaging-guided interpolation. Our numerical experiments using synthetic and real data sets indicate the efficiency and robustness of image-guided rays for the selected applications.

The true-amplitude (TA) imaging condition for reverse-time migration (RTM) is based on a combination of temporal and spatial derivatives of the upand downgoing wavefields. By means of partial integrations (or redistribution of the frequency factors in the frequency domain, we derive several alternative expressions for this imaging condition. Interestingly, the temporal derivatives can be completely replaced by spatial derivatives and temporal integrations. In this way, one version of the TA imaging condition makes use of the Laplacian operator, in this way relating to a common way of removing backscattering artifacts in RTM. We demonstrate by means of numerical examples using the Marmousi and Sigsbee2A data that the quality of the migrated image strongly depends on the version chosen for implementation. The best quality is achieved with a version that combines second derivatives of the source wavefield with the Laplacian operator. INTRODUCTION Reverse-time migration (RTM) is a seismic imaging method based on the full (two-way) wave equation (Schultz and Sherwood, 1980; McMechan, 1983; Baysal et al., 1983; Sun and McMechan, 2001; Yan and Sava, 2008). In the same way as other wave-equation based migration techniques, it makes use of an image condition, the most basic form of which is simple cross-correlation of the upand downgoing wavefields (Claerbout, 1971). In the early days of seismic imaging, RTM was not of much practical use because of its high computational cost and the presence of strong low-frequency artifacts from backscattering if the velocity model contains sharp velocity contrasts. Its use has gained much popularity in the first decade of this century, after the advent of more powerful computers and new technologies to remove the backscattering artifacts (Yoon et al., 2004; Fletcher et al., 2006; Guitton et al., 2006). Most of these techniques rely on modified imaging conditions (see, e.g., Yoon and Marfurt, 2006; Costa et al., 2009; Luo et al., 2009, 2010; Zhu et al., 2009). In the same context, Kiyashchenko et al. (2007) and Op’t Root et al. (2012) derive the true-amplitude (TA) imaging condition for reverse-time migration (RTM). The purpose of the TA imaging condition is to remove the backscattering artifacts and to provide image amplitudes that are proportional to reflection coefficients. Unfortunately, in the form presented by these authors, the TA imaging condition does not allow for an efficient implementation in the time domain. For a time-domain implementation, it must be recast into a different form. In this paper, we derive a number of theoretically equivalent forms of the approximate TA imaging condition that can be efficiently implemented in the time domain. In a similar way to Douma et al. (2010), we show that the true-amplitude imaging condition for RTM can be reformulated into a version containing the Laplacian operator. This operator is frequently used in seismic imaging without a profound theoretical basis to remove the low-frequency backscattering artifacts from RTM images. We numerically evaluate the derived time-domain versions of the TA imaging condition by comparing the resulting migrated images of the Marmousi and Sigsbee2A data. 16 Annual WIT report 2019 TRUE-AMPLITUDE IMAGING CONDITION According to Op’t Root et al. (2012), the TA imaging condition is given in the frequency domain by Ir(x) = 1 2π ∑ s ∫ ω dω 1 (−iω)PsPs [ PsPr − c(x) ω2 ∇Ps · ∇Pr ] , (1) where Ps = Ps(ω,x;xs) and Pr = Pr(ω,x;xs) are the (downgoing) source and (upgoing) receiver wavefields for a source at xs, downward continued to the imaging point x. Moreover, the bar over a symbol denotes the complex conjugate operation. Equation (1) is slightly different from the one of Op’t Root et al. (2012). For simplicity, we have assumed that the source wavelet is a (possibly band-limited) delta-function, the effects of which are acceptable in the final migrated image. Therefore, we have combined in equation (1) the Green’s function and source wavelet in the formula of Op’t Root et al. (2012) into the source wavefield Ps. Moreover, we have made the sum over all sources explicit. Finally, the different sign of the factor (−iω) in the denominator of the above equation results from our use of the following definition of the Fourier transform pair, f(ω) = ∫ ∞ −∞ dt f(t)e e f(t) = 1 2π ∫ ∞ −∞ dω f(ω)e−iωt . (2) When trying to implement the TA imaging condition in the time domain, one recognizes that its basic form, equation (1) is not very favorable. For an efficient time-domain implementation, it must be recast into a more adequate form. Below, we derive a number of theoretically equivalent forms. We then compare them numerically by looking at the corresponding images. Implementational forms The advantage of the frequency domain representation of the TA imaging condition in equation (1) is that the time derivatives are represented by factors (−iω). Thus, it immediately allows us to recognize that these factors can be rather freely redistributed among the wavefield terms. Making use of this freedom, our first rewrite moves the ω in the denominator of the spatial-derivatives term, where it would represent a two-fold time integration, to a position in front of the parentheses. This results in Ir(x) = ∑ s ∫ ω dω 1 (−iω)(−ω2)PsPs [ −ωPsPr + c(x)∇Ps · ∇Pr ] ,

As it becomes increasingly important for migration methods to provide compensation for illumination and amplitude recovery, particularly, if the goal is to successfully recover a reflectivity of the medium, an imaging condition like simple crosscorrelation, which destroys the image amplitude, is unacceptable. For this reason, several alternative forms of imaging conditions have emerged in the recent past that are evaluated by the quality of the output amplitudes and artifacts produced. In this work we study a set of imaging conditions with illumination compensation. We also present new stabilized least-squares image conditions and compare them to previously proposed forms. Our numerical experiments on a simple horizontal interface model using a vertically inhomogeneous velocity model and on the Marmousi data set show that they produce satisfactory results. The general observation from the overall comparison is that the stabilized total least-squares imaging condition produced the best image with the least migration artifacts and the least affected migration amplitudes. Its image quality comes very close to the one of the simple crosscorrelation imaging condition, however with correctly recovered relative amplitudes.

Kirchhoff-type, isochrone-stack demigration is the natural asymptotic inverse to classical Kirchhoff or diffraction-stack migration. Both stacking operations can be performed in true amplitude by an appropriate selection of weight functions. As Kirchhoff migration is usually understood as the inverse process to Kirchhoff modeling, the natural question arises whether Kirchhoff demigration is identical to seismic forward modeling. The answer is that it is not, but these processes are closely enough related to enable the use of demigration for modeling purposes. All that has to be done is to implicitly construct a depth section as if obtained from a previous true-amplitude Kirchhoff migration.

Path-integral migration is a method for creating a migrated image without previous knowledge of the true velocity model by summing the migrated images from a representative set of velocity models. This concept can be expanded to automatically extract a velocity model, a technic called Migration Velocity Analysis by Double Path-Integral Migration (MVA by DPIM). In MVA by DPIM, a second, weighted image is created, with the weight containing the velocities used in the individual migrations. Division of the two images provides the velocity model. Here we discuss several practical aspects of implementing MVA by DPIM, ranging from the parametrization of the weight function to the stabilization of the division and the selection of only meaningful velocities. By means of tests on a wide range of velocity models, we find a robust implementation of the method. Presentation Date: Wednesday, September 18, 2019 Session Start Time: 1:50 PM Presentation Start Time: 2:40 PM Location: Poster Station 11 Presentation Type: Poster

We introduce a data-driven stacking technique that transforms 2D/2.5D prestack multicoverage data into a common-offset (CO) section. We refer to this new process, which is based on the offsetcontinuation operation (OCO), as offset-continuation stacking or briefly OCO stack. Similarly to the CMP and CRS stacks, the OCO stack does not rely on an a-priori velocity model but provides velocity information itself. The original OCO method is a seismic configuration transform designed to simulate a seismic section as if obtained with a certain source-receiver offset using the data measured with another offset. Since OCO is dependent on the velocity model used in the process, it can be combined with stacking techniques for a set of models, thus allowing for the extraction of velocity information. The algorithm is based on so-called OCO trajectories, which are related to the concepts of image waves and velocity rays. We theoretically relate the OCO trajectories to the kinematic properties of OCO image waves that describe the continuous transformation of the common-offset reflection event from one offset to another. Based on OCO trajectories, we then formulate a horizon-based velocity analysis method, where root mean square (RMS) velocities and local event slopes are determined by stacking along event horizons. INTRODUCTION By definition, the Offset-Continuation Operation (OCO) is an operator that transforms common offset (CO) seismic gathers from one constant offset to another (Deregowski and Rocca, 1981). It is an important tool for imaging in a complex medium. Possible applications of OCO include velocity analysis, commonreflection point (CRP) stacking, dip moveout (DMO), migration to zero offset (MZO), interpolation of missing data, amplitude variation with offset (AVO) studies, and geometrical-spreading correction (see, e.g., Salvador and Savelli, 1982; Bolondi et al., 1982, 1984; Fomel, 1994, 2003; Santos et al., 1997). Since OCO is a configuration transform, its objective is to simulate a seismic section using as input the data measured with another configuration. As discussed by Hubral et al. (1996a) and mathematically demonstranted by Tygel et al. (1996), any configuration transform can be thought of as being composed of a migration and a subsequent demigration after changing a configuration parameter. Configurations transforms have already been used for several purposes in seismic processing such as MZO (Tygel et al., 1998; Bleistein et al., 1999), source continuation operation (SCO) (Bagaini and Spagnolini, 1993, 1996), azimuth moveout (AMO) (Biondi et al., 1998), DMO (Hale, 1984; Canning and Gardner, 1996; Collins, 1997; Black et al., 1993), common-source (CS)-DMO (Schleicher and Bagaini, 2004), data reconstruction (Bagaini et al., 1994; Stolt, 2002; Chemingui and Biondi, 2002), and velocity analysis (Silva, 2005; Coimbra et al., 2012). For data of very low signal-to-noise ratio (S/N) or acquisitions with very low fold, conventional common-midpoint (CMP) processing might not provide stacked sections of sufficient quality. In such situations, alternative processing sequences are necessary to improve the data quality. The OCO stack represents such an alternative path for the processing of reflection-seismic data. Its key element is the construction of common-offset stacked sections together with coherency sections and sections of kinematic 60 Annual WIT report 2012 and dynamic wavefield attributes. The OCO stacking surface is composed of so-called OCO trajectories (Coimbra et al., 2012). Such a trajectory requires only two parameters (local event slope and stacking velocity) to describe the seismic reflection event in the multi-coverage data. Neighbouring trajectories can be located by event tracking in the stacked section or described by a third curvature-related parameter. Using these parameters, the method stacks the data along a predicted traveltime curve that approximates the CRP event. Since the parameters, and thus the predicted traveltime curve, are updated from the data at each offset, the approximation is better than by conventional methods that adjust the approximate traveltime expression at some initial point. The purpose of this paper is to establish a consistent processing chain that is based entirely on the OCO stack, relying on identical assumptions at all steps. METHOD The OCO stack is a multiparameter stacking procedure similar to its relatives, the CMP and CRS stacks and multifocusing. It automatically determines stacking attributes based on a coherence measure applied at every common-offset sample of the data. Since these attributes vary with time for the same event, the OCO stacked section is free of normal moveout (NMO)-stretch (Perroud and Tygel, 2004). The main advantages of the OCO stack are twofold. Firstly, it is not limited to a zero-offset stacked section like the CMP stack. Secondly, for the 2D/2.5D case as discussed here, the OCO stack needs at most two parameter in addition to stacking velocity, even for the construction of stacked common-offset sections. There are other multiparameter stacking methods (Gelchinsky et al., 1999; Jäger et al., 2001; Zhang et al., 2001; Hertweck et al., 2007; Fomel and Kazinnik, 2012), which are based on stacks data from multiple CMP locations. As a result, they considerably improve the signal-to-noise ratio. However, these methods require the estimation of more data parameters than conventional CMP processing, in addition to the conventional stacking velocity. For instance, the zero-offset CRS method requires two additional parameters and common-offset CRS requires four of them. Moreover, due to multi-coverage some events such as diffractions, far-offset faults and strong dips can disappear. In this section, we derive the theoretical basis for the OCO stack. It is based on the kinematic behaviour of the OCO transformation as described by the OCO image-wave equation (Hubral et al., 1996b). Image-wave for OCO The OCO image-wave equation was derived through image-wave theory from the kinematic behaviour of the OCO transformation (Hubral et al., 1996b). It is a second order linear partial differential equation, which can be written as ht ( ∂U ∂h2 + 4 V 2 ∂U ∂t2 ) + ( t + 4h V 2 ) ∂U ∂h∂t − ht U ∂ξ2 = Θ ( ξ, t, h, U, ∂U ∂ξ , ∂U ∂t , ∂U ∂h ) . (1) Equation (1) describes the behavior of an artificial (non-physical) process of transforming reflection seismic data U(ξ, t, h) in the offset-midpoint-time domain as a certain kind of “wave propagation”. In this case it is the record of the seismic reflection that “propagates”as a function of half-offset h. In equation (1), ξ and t are the midpoint and time coordinate of the reflection event under consideration. The velocity V is assumed to be a constant average velocity that is known a priori. We will refer to V as the OCO velocity. Its relationship to the RMS velocity is discussed below. Equation (1) belongs to the class of linear hyperbolic equations, when t > 0 and h 6= 0, with the halfoffset h acting as propagation variable (i.e., equivalent to time in conventional wave propagation). Equation (1) describes a wave-like propagation in the offset direction that Hubral et al. (1996b) termed image-wave propagation. We use the OCO image-wave equation (1) to obtain the trajectory of single point under variation of the half-offset. Formally, we can think of the solution to equation (1) as being approximated by an expression that is analogous to the one used in ray theory, i.e., the leading term of a high-frequency asymptotic (WKBJtype) approximation for a reflected wave recorded on a seismogram of the form U(ξ, t, h) = A(ξ, t)F (h−H(ξ, t)), (2) Annual WIT report 2012 61

Amplitude anomalies imply strong amplitude variations over relatively short distances. Thus, a question fundamental to their quantitative interpretation asks for the influence of lower amplitudes on nearby higher amplitudes and vice versa, particularly in the context of post-migration AVO analysis. This question is directly related to the resolving power of seismic migration as a function of source-receiver offset. Horizontal resolution can be quantified in the time domain by means of the region around the migrated reflection point that is influenced by the migrated elementary wave. To obtain a numerical estimate for the mentioned zone of horizontal influence after migration, we investigate the migration output at a chosen depth point in the vicinity of the specular reflection point for a simple model of a horizontal interface with a vertical fault. We find that the region of influence before migration is well approximated by the projected Fresnel zone, where the half-period is replaced by an effective wavelet length. Spatial resolution after migration depends on the reflection angle rather than source-receiver offset. Thus, in principle, achievable resolution does not depend on reflector depth. As expected, migration improves the resolution for the usual seismic range of offsets. The achievable resolution remains almost the same for reflection angles up to about 30 degrees, but then strongly decreases. In consequence, a large-offset AVO/AVA analysis may lead to wrong results. INTRODUCTION Amplitude anomalies along a seismic reflector are a principal hydrocarbon indicator. An amplitudevariations-with-offset (AVO) analysis of the bright or dim spot can often increase the usefulness of these indicators. By their very nature, amplitude anomalies are spatially localized. Therefore, strong amplitude variations along the reflector occur at their boundaries, often over rather short distances. Due to the limited frequency content of the seismic waves, this means that the amplitudes within the anomaly are influenced by the different adjacent amplitudes. The quesAnnual WIT report 2001 73 migrated reflector image M x S G r

The application of an imaging condition in wave equation shot profile migration is important to provide illumination compensation and amplitude recovery. Particularly for true-amplitude waveequation migration algorithms, a stable imaging condition is essential to successfully recover the medium reflectivity. We continue our study of a set of image conditions with illumination compensation by application to the Marmousi data. The imaging conditions are evaluated by the quality of the stacked migrated sections. The most stable of the tested imaging condition with illumination compensation divides the crosscorrelation of the upand downgoing wavefields by the autocorrelation of the downgoing wavefield. Smoothing imaging conditions, which work perfectly in vertically inhomogeneous media, tend to fail for laterally varying velocities.

Anisotropic cracked media have been widely investigated in many theoretical and experimental studies. In this work, we have performed ultrasonic surveys to investigate the influence of source frequency on elastic parameters (the Thomsen parameter γ and shear-wave attenuation) of fractured anisotropic media. Under controlled conditions, we prepared anisotropic models containing penny-shaped rubber inclusions in a solid epoxy resin matrix with crack density that ranges from 0 to 6.2 %. Two of the three cracked models have 10 layers and the last one has 17 layers. The number of uniform rubber inclusions per layer was from 0 up to 100. S-wave splitting measurements have shown that scattering effects are more prominent in models where the crack aperture to seismic wavelength ratio ranges from 1.6 to 13.3 than in other models where the ratio varied from 2.3 to 23. The model with large cracks gave a magnitude of attenuation 3 times higher compared with another model that had small inclusions. These results indicate that elastic scattering, intrinsic and scattering attenuation (Q−1 in and Q−1 s respectively), and velocity dispersion directly interfere in shear wave splitting, which in turn is a function of crack size and source frequency. INTRODUCTION Wave propagation in anisotropic cracked and fractured media has motived many studies in seismic exploration of hydrocarbons reservoirs. Because of the geologic complexities exhibited by anisotropic media, reliable conclusions about elastic properties are usually difficult to achieve with accuracy from field data. On the other hand, laboratory measurements have been shown to be a useful tool for modeling conditions present in the field, helping to reduce uncertainty about elastic parameters in numerical methods. It is well known that numerical simulation of cracked media is computationally and mathematically expensive and intense (Hudson, 1981; Crampin, 1981; Hudson et al., 2001). Furthermore, when scattering effects are taken into account, these costs become even more significant (Willis, 1964; Mal, 1970; Yang and Turner, 2003, 2005). Nonetheless, some difficulties of anisotropic modeling can be overcome using experimental scaled physical modeling. Assad et al. (1992, 1996), Wei (2004) and Wei et al. (2007) established an experimental relationship between crack density and shear velocity based on theoretical predictions by Hudson (1981). Melia and Carison (1984) carried out a series of experiments in anisotropic samples to investigate P-wave dispersion in anisotropic layered media as a function of the concentration of different layered materials as well as the thickness of the layers. Based on the same approach, Marion et al. (1994) and Rio et al. (1996) showed the influence of short and long wavelengths in stratified media as well as wave velocity dispersion and multiple scattering. Other sets of experimental observations performed by Rathore et al. (1995) and Peacock et al. (1994) demonstrated the feasibility of the ultrasonic approach to investigate artificially cracked porous media. Annual WIT report 2011 205 Figure 1: (a) From right to left: Reference model M1 (uncracked) and cracked models M2, M3; (b) model M4. Also shown are the orientations of the coordinate systems. All wave measurements were made in the Y direction. Using experimental data obtained by Rathore et al. (1995), the theoretical predictions of Thomsen (1995) for aligned cracks in porous rock received a strong support. More recently, experiments by Tillotson et al. (2011) have suggested the possible use of shear wave data to discriminate fluids on the basis of viscosity variations. In anisotropic cracked media, the frequency response is influenced by the size of the heterogeneities. However, quantification of this influence is still desirable. To better understand the influence of frequency on cracked materials, we conducted a series of experiments aimed at extending previous approaches by using a shear-wave source with different frequencies: low frequency (LF = 90 kHz), intermediate frequency (IF = 431 kHz) and high frequency (HF = 840 kHz). We carried out experiments on a reference model without inclusions and three other models with different inclusion sizes, thereby simulating different crack densities. In this arrangement, shear-wave splitting was observed with different magnitudes as a function of frequency. Our results show that effects associated with intrinsic (Q−1) and scattering (Q−1 s ) attenuation (Gorich and Muller, 1987; Tselentis, 1998) interfere directly with shear wave splitting, which in turn is related to crack density. Furthermore, we observed that the anisotropic parameter γ (Thomsen, 1986) varies with frequency and crack size. For these purposes, we quantified attenuation using the frequency shift method (Quan and Harris, 1997). EXPERIMENTAL PROCEDURE The construction of the cracked samples as well as the ultrasonic measurements were carried out at the Allied Geophysical Laboratories (AGL) at the University of Houston. Model preparation Under controlled conditions, we constructed three cracked models (M2, M3, and M4) with different crack densities and one uncracked model (M1) for reference. Pictures of all models are shown in Figure 1. Model M4 has five different points that can be analyzed. The same distance between layers (0.5 cm for M2 and M4 and 0.25 cm for M3) was ensured by using the same volume of epoxy resin poured for each layer. 206 Annual WIT report 2011 Figure 2: (a) Device developed for S-wave polarization rotation. (b) Sketch of experiment used for seismogram records. Each layer with inclusions was added to the model and air was extracted using a vacuum pump to avoid inhomogeneities in the epoxy resin. The crack density ε in the cracked models was determined by ε = Nπrh V , (1) where N is total number of inclusions, r is their radius, h is inclusions’ thickness (aperture of cracks), and, finally, V is the total model volume. Equation (1) is a modification of the relation of Hudson (1981) for crack density estimation. The ratio of compressional wave velocity between solid epoxy and neoprene was around 1.5 and for solid epoxy and silicone rubber was about 2.25. The S-wave velocity in rubber was difficult to determine because of the low shear modulus of this material. The parameters of the included rubber cracks in each model are displayed in Table 1. Model Crack Measuring Number Diameter Aperture Cracks Aspect density (%) length (cm) of layers (cm) (cm) per layer ratio M1 Isotropic 7.31 ± 0.02 0 0 0 M2 4.5 7.29 ± 0.02 10 0.7 0.091 36 0.13 M3 3.8 7.32 ± 0.02 17 0.4 0.051 90 0.12 M4-1 6.0 7.64 ± 0.02 10 0.7 0.091 30 0.13 M4-3 5.2 7.74 ± 0.02 10 0.44 0.091 80 0.20 M4-5 4.2 7.74 ± 0.02 10 0.32 0.091 100 0.28 Table 1: Physical parameters of models M1, M2, M3 and M4 Ultrasonic measurements Over these models, we carried out ultrasonic measurements using the Ultrasonic Research System at AGL with the pulse transmission technique. The sampling rate per channel for all experiments was 10 MHz. Figure 2a shows a device developed for recording S-wave seismograms. The source and receiver transducers were arranged on opposing sides of the model, separated by measuring length (see Table 1). The initial shear-wave polarization was parallel to the cracks. Changes in polarization were achieved by rotating both transducers by 10 degrees at a time until polarization was again parallel (i.e., 0 to 180 degrees) to the XZ plane (see Figure 2b). In total, 19 traces were recorded in each seismic section with 20fold stack to eliminate ambient noise. The polarizations of 0 and 180 degrees correspond to the fast S-wave (S1) and 90 degrees corresponds to the slow S-wave (S2). Figure 3a and b shows the S-wave signature sources and Fourier amplitude spectra of the three sources used to obtain the data shown in this paper. We performed a Gaussian non-linear fit in each frequency Annual WIT report 2011 207 0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.3 0.6 0.9 1.2 0.0 0.3 0.6 0.9 1.2 -1.0 -0.5 0.0 0.5 1.00 5 10 15 20 25 30 S wave 89 kHz S wave 386 kHz S wave 805 kHz

Three-dimensional wave-equation migration techniques are still quite expensive because of the huge matrices that need to be inverted. Several techniques have been proposed to reduce this cost by splitting the full 3D problem into a sequence of 2D problems. We compare the performance of splitting techniques for stable 3D Fourier Finite-Difference (FFD) migration techniques in terms of image quality and computational cost. The FFD methods are complex Padé FFD and FFD plus interpolation, and the compared splitting techniques are two and four-way splitting as well as alternating four-way splitting, i.e., splitting into the coordinate directions at one depth and the diagonal directions at the next depth level. From numerical examples in homogeneous and inhomogeneous media, we conclude that alternate four-way splitting yields results of the same quality as full four-way splitting at the cost of two-way splitting. INTRODUCTION Because of its superiority in areas of complex geology, wave-equation migration is substituting Kirchhoff migration in practice. However, while Kirchhoff migration counts on more than 30 years of technological development, wave-equation migration methods still need to be improved in various aspects. One of these aspects is the efficient implementation of three-dimensional wave-equation migration. The application of a three-dimensional wave-equation migration technique adds the problem of computational cost to those of stability and precision of the chosen migration algorithm. To speed up migration techniques like finite-difference (FD) (Claerbout, 1971) or Fourier finite-difference (FFD) migration (Ristow and Rühl, 1994), a technique known as splitting is frequently used. In this context, splitting means the separation of a single-step 3D migration into two 2D passes within planes parallel to the horizontal coordinate axes, usually the inline and crossline directions (Brown, 1983). When the splitting is applied to the implicit FD migration operator in such a way that the resulting equations are solved alternatingly in the inline and crossline directions, the resulting FD scheme is known as an Alternating-Direction-Implicit (ADI) scheme. This procedure has the drawback of being incorrect for strongly dipping reflectors, resulting in large positioning errors for this type of reflectors when the dip direction is away from the coordinate directions and thus outside the migration planes. This imprecision leads to numerical anisotropy, i.e., a migration operator that acts quite differently in different directions. To improve this behaviour while retaining the advantages of a rather low computation cost, different procedures have been proposed over the years. Ristow (1980, see also Ristow and Rühl, 1997) proposed to perform, in addition to the 2D migration in the coordinate planes, also 2D migrations in the diagonal directions between the coordinate axes. Kitchenside (1988) used phase-shift migration plus an additional FD propagation step of the residual field to reduce the splitting error. Graves and Clayton (1990) proposed the implementation of a phase-correction operator using FD and incorporating a damping function to guarantee the stability of the 3D FD migration scheme. 28 Annual WIT report 2010 Inverting the idea of Kitchenside (1988), who propagated the field using phase shift and the residual using FD, Li (1991) proposed to use conventional FD migration plus a residual field correction by phase shift to improve the migrated image quality. Without any need to modify the conventional 3D FD migration, the Li correction adds a phase-shift filter at certain steps of the downward extrapolation. This technique corrects not only for the splitting error, but also for the positioning error of steeply dipping reflectors. Collino and Joly (1995) solved a family of new 3D one-way wave equations by the ADI method. These equations significantly reduce the numerical anisotropy, but are approximately four times as expensive as conventional two-way splitting. Wang (2001) developed an alternative method to improve the precision of the FD solution of the one-way wave equation. To guarantee stability and efficiency, he keeps the implicit FD scheme and the alternation of directions, but interpolates between the ADI solution and the wavefield before each step of extrapolation. He calls the resulting method ADI plus interpolation (ADIPI). As a drawback, this ADIPI can produce instabilities in the presence of strong lateral velocity variations. Zhou and McMechan (1997) proposed a 45◦ one-way wave equation that can be expressed as a system of differential equations of first and second order (Zhang et al., 1988), and factored into a product of two one-dimensional terms corresponding to the lateral directions (Mitchell and Griffiths, 1980; Graves and Clayton, 1990). A big asset of the method is that the conventional FD extrapolation can be used with very little modification. In this way, the efficiency of conventional splitting is preserved, without adding the necessity of any error compensation. However, the method is also instable for strong lateral velocity contrasts and needs rather heavy model smoothing. Biondi (2002) showed that FFD migration is more precise than other methods that use implicit finite differences like pseudoscreen propagators (Jin et al., 1999) and high-angle screen propagators (Xie and Wu, 1998). Given that the computational complexity of all three methods is approximately the same, FFD migration is more attractive than the others. Unfortunately, when conventional FFD migration is applied in the presence of strong velocity contrasts, it can generate numerical instabilities, too. To overcome the problem of instabilities in models with strong lateral velocity contrasts, Biondi (2002) presented a correction to the FFD method that avoids stability problems. To derive it, he adapted a theory of Godfrey et al. (1979) and Brown (1979), which improves the stability of the 45◦ equation. The corrected FFD method is unconditionally stable for arbitrary velocity variations, as much in the velocity model as in the reference velocity. Particularly, and differently from conventional FFD migration, it is unconditionally stable even if the reference velocity is smaller than the model velocity. This new property allows for the application of the interpolation technique, conventionally used to improve phase-shift and split-step migration (Gazdag and Sguazzero, 1984) but impossible in FFD migration, because it needs propagation with a larger and a smaller reference velocity. The resulting migration technique is called FFD plus interpolation, or shortly FFDPI. Another, computationally less expensive method to stabilize FFD migration in the presence of strong lateral velocity contrasts was proposed by Amazonas et al. (2007). It substitutes the real Padé approximation (Bamberger et al., 1988) used in the derivation of FFD migration (Ristow and Rühl, 1994) by its complex version (Millinazzo et al., 1997). In this way, the incorrect treatment of near horizontal and slightly evanescent waves of the real Páde approximation is improved, leading to a more stable FFD algorithm, shortly referred to as complex Padé FFD (CPFFD) migration. In this work, we study possibilities of efficiently implementing these stable FFD migration techniques in 3D. We implemented and compared splitting techniques for FFDPI (Biondi, 2002) and CPFFD (Amazonas et al., 2007) migration. Our numerical tests indicate that a very robust, highly efficient, and satisfactorily accurate method is alternate four-way splitting, i.e., splitting into the coordinate directions at one extrapolation step and into the diagonal directions at the next step. THEORETICAL BACKGROUND The one-way wave equation The one-way wave equation (Leontovich and Fock, 1946) can be derived starting from the scalar wave equation, which for a homogeneous medium is given by ∂p(x, t) ∂z2 + ∂p(x, t) ∂x2 + ∂p(x, t) ∂y2 − 1 c2 ∂p(x, t) ∂t2 = 0 , (1) Annual WIT report 2010 29 where p(x, t) is the scalar wavefield and c = c(x) is the spatially varying wave velocity. For moderately varying media, where velocity derivatives can be neglected, Fourier transform in time and horizontal coordinates x and y allows to represent equation (1) as ∂P (kx, ky, z, ω) ∂z2 − (−iω) 2 c2 ( 1− c 2 ω2 (k x + k 2 y) ) P (kx, ky, z, ω) = 0 . (2) Equation (2) can be factorized into [ ∂P (kx, ky, z, ω) ∂z − (−iω) c √ 1− c 2 ω2 (k2 x + k2 y) ] × [ ∂P (y, kx, ω) ∂z + (−iω) c √ 1− c 2 ω2 (k2 x + k2 y) ] P (kx, ky, z, ω) = 0 . (3) The two differential operators in equation (3) represent, when taken alone, one-way wave equations that describe up and downgoing waves. For migration, the one-way wave equation of interest is the one describing downgoing waves, i.e., ∂P (kx, ky, z, ω) ∂z = (−iω) c √ 1− c 2 ω2 (k2 x + k2 y)P (kx, ky, z, , ω) . (4) Inverse Fourier transform in the horizontal wavenumbers kx and ky yields then formally ∂P (x, ω) ∂z = (−iω) c(x) √ 1 + c2(x) ω2 ( ∂2 ∂x2 + ∂2 ∂y2 ) P (x, ω) . (5) The actual restrictions that apply to equation (5) in inhomogeneous media are much less severe than the above derivation indicates. Of course, for the formal representation (5) to make practical sense, the square root of the differential operator needs to be approximated in terms of numerically executable operations. Expansion of the square root A well-used possibility for the approximation of the square root in the one-way wave equation (5) in terms of numerically executable operations is an expansion into a Padé series (Bamberger et al., 1988) √ 1 + Z ≈ 1 + N ∑ n=1 anZ 1 + bnZ (6) where the Padé coefficients are an = 2 2N + 1 sin ( nπ 2N + 1 ) , and bn = cos ( nπ 2N + 1 ) . (7) This approximation is used in most practical FD migration schemes. Depending on the number N of terms used in the expansion, it gives rise to the so-called 15◦, 45◦, or 60◦ migrations. However, when the interest is on accurate imaging up to very high

Seismic Migration by downward continuation using the one-way wave equation approximations has two shortcomings: imaging steep dip reflectors and handling evanescent waves. Complex Padé approximations allow a better treatment of evanescent modes, stabilizing finite-difference migration without requiring special treatment for the migration domain boundaries. Imaging steep dip reflectors can be improved using several terms in the Padé expansion. We discuss the implementation and evaluation of wide-angle complex Padé approximations for finite-difference and Fourier finite-difference migration methods. The dispersion relation and the impulsive response of the migration operator provide criteria to select the number of terms and coefficients in the Padé expansion. This assures stability for a prescribed maximum propagation direction. The implementations are validated on the Marmousi model dataset and SEG/EAGE salt model data. INTRODUCTION Wave equation migration algorithms have a better performance than ray-based migration when the velocity model has strong lateral velocity variations. Among several existing algorithms for wave-equation migration, finite-difference (FD) and Fourier finite-difference (FFD) migrations (Ristow and Rühl, 1994) can provide wide-angle approximations for the one-way continuation operators, thus improving the imaging of steep dip reflectors. However, standard (real-valued) FD and FFD migrations cannot handle evanescent waves correctly (Millinazzo et al., 1997). As a consequence, FFD algorithms tend to become numerically unstable in the presence of high velocity variations (Biondi, 2002). To overcome this limitation, Biondi (2002) proposed an unconditionally stable extension for the FFD algorithm. Earlier, Millinazzo et al. (1997) proposed a different approach to treating these evanescent modes in ocean acoustic applications. They introduced an extension of the Padé approximation which they called complex Padé. The complex Padé expansion was used previously in applied geophysics. Zhang et al. (2003) used the method in finite-difference migration. However, their implementation were not suited for wide angles. Later, Zhang et al. (2004) proposed a split-step migration based on complex Padé. In this paper, we study the use of the complex Padé expansion for wide-angle FD and FFD pre-stack depth migration algorithms. The expansion is evaluated numerically, comparing its approximation to the exact one-way operator. Based on studying the impulse response of the migration operator, we propose a prescription to choose parameters for wide-angle complex FD and FFD algorithms. The algorithms are validated on synthetic datasets from the Marmousi and SEG/EAGE salt models. Annual WIT report 2006 173 METHODOLOGY Complex Padé approximation The one-way wave equation for downward continuation of the acoustic wavefield reads (Ristow and Rühl, 1994) ∂P (x, ω) ∂x3 = (−iω) c(x) √ 1 + c2(x) ω2 ∂2 ∂x1 P (x, ω) . (1) where P (x, ω) is the pressure wavefield, c(x) is the medium propagation velocity. For vertically inhomogeneous media, the operator above has an exact representation in the Fourier domain (Gazdag (1978)). For laterally inhomogeneous media, a formal representation for this operator is based on the Padé expansion (Bamberger et al., 1988) √ 1 + Z = 1 + N ∑ n=1 anZ 1 + bnZ , (2) where Z ≡ c 2(x) ω2 ∂2 ∂x1 . The coefficients an and bn are (Bamberger et al., 1988) an = 2 2N + 1 sin nπ 2N + 1 and bn = cos nπ 2N + 1 . (3) If Z < −1 in equation (2), the left side is a pure imaginary number while the right side remains a realvalued quantity. In other words, the approximation breaks down. Physically, this means that representation (2) cannot properly handle evanescent modes. This causes numerical instabilities and is responsible for the unstable behavior of the FFD algorithm in the presence of high velocity variations (Biondi, 2002). To overcome these limitation, Millinazzo et al. (1997) proposed a complex representation of the Padé expansion in equation (2). They achieve this goal by rotating the branch cut of the square root in the complex plane. Their final expression is √ 1 + Z ≈ Rα,N (Z) = C0 + N ∑ n=1 AnZ 1 + BnZ , (4) where An ≡ ane −iα/2 [1 + bn(e−iα − 1)] , Bn ≡ bne −iα 1 + bn(e−iα − 1) , and C0 = e [ 1 + N ∑ n=1 an(e − 1) [1 + bn(e−iα − 1)] ] ; , with an and bn as defined in equation (3). An and Bn are the complex Padé coefficients and α is the rotation angle of the branch cut of the square root in the complex plane. Complex FD and FFD migration We use the complex Padé approximation (4) to represent the one-way continuation operator. Using this approximation the downward continuation operator for finite difference migration is ∂P (x, ω) ∂x3 = (−iω) c(x) C0 + N ∑ n=1 An c2(x) ω2 ∂2 ∂x1 1 + Bn c2(x) ω2 ∂2 ∂x1  P (x, ω) . (5) The FFD approximation for the downward continuation operator is deduced following the derivation proposed by Ristow and Rühl (1994). The final result is p √ 1 + X2 ≈ √ 1 + p2X2 + C0(p− 1) + N ∑ n=1 Anp(1− p)X 1 + σBnX , (6) 174 Annual WIT report 2006 −3 −2 −1 0 1 2 3 −1 −0.5 0 0.5 1 1.5 k1/k k3 /k Real part FD complex Pade Exact operator −3 −2 −1 0 1 2 3 0 0.5 1 1.5 2 k1/k k3 /k Imaginary part Figure 1: Complex Padé FD approximation for the dispersion relation of the one-way wave equation, computed with three terms and α = 90. where p ≡ cr c(x) is the ratio between the actual propagation velocity, c(x), and the propagation velocity in an homogeneous background medium, cr. Moreover, X ≡ ( c ω )2 ∂2 ∂x1 and σ = 1 + p + p. Based on numerical experiments comparing the exact operator and the FFD approximation, we propose to use σ = 1 + p for wide-angle approximations, instead of σ = 1 + p + p. NUMERICAL EVALUATIONS To better understand the involved approximations, we numerically evaluate the dispersion relations of the wide-angle FD and FFD approximations and compare them with the exact dispersion relation. We complement our numerical evaluation computing the impulse response of our approximations in a homogeneous medium. Finally, we compare the impulse response of the proposed wide-angle complex FD and FFD algorithms with the unconditionally FFD algorithm proposed by Biondi (2002). Dispersion relation ure 1 shows the comparison of the exact dispersion relation with its FD approximation 2. The FD approximation was calculated using the first three terms of the series with a rotation angle α = 90. The complex Padé approximation fits the real and imaginary part of the dispersion relation almost perfectly. In other words, it correctly represents the evanescent modes. ure 2 shows the comparison of three FFD approximations for the one-way wave equation dispersion relation, with the exact dispersion curve: FFD using real coefficients, complex FFD with σ = 1 + p + p and complex FFD using σ = 1 + p. The FFD approximations were determined using p = 0.5 and three terms in the Padé expansion. For the complex Padé approximation, we used a rotation angle of α = 45. The real part of the FFD operator using the real-valued Padé expansion (2) is clearly affected by spurious oscillations in the evanescent region. The complex FFD attenuates evanescent modes and, moreover, Annual WIT report 2006 175 −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1