Sea surface shape derivation above the seismic streamer

The rough sea surface causes perturbations in the seismic data that can be significant for time‐lapse studies. The perturbations arise because the reflection response of the non‐flat sea perturbs the seismic wavelet. In order to remove these perturbations from the received seismic data, special deconvolution methods can be used, but these methods require, as input, the time varying wave elevation above each hydrophone in the streamer. In addition, the vertical displacement of the streamer itself must also be known at the position of each hydrophone and at all times. This information is not available in conventional seismic acquisition. However, it can be obtained from the hydrophone measurements provided that the hydrophones are recorded individually (not grouped), that the recording bandwidth is extended down to 0.05 Hz and that data are recorded without gaps between the shot records. The sea surface elevation, and also the wave‐induced vertical displacement of the streamer, can be determined from the time‐varying pressure that the sea waves cause in the hydrophone measurements. When this was done experimentally, using a single sensor seismic streamer without a conventional low cut filter, the wave induced pressure variations were easily detected. The inversion of these experimental data gives results for the sea surface elevation that are consistent with the weather and sea state at the time of acquisition. A high tension approximation allows a simplified solution of the equations that does not demand a knowledge of the streamer tension. However, best results at the tail end of the streamer are obtained using the general equation.     

Operator-oriented CRS interpolation

In common‐reflection‐surface imaging the reflection arrival time field is parameterized by operators that are of higher dimension or order than in conventional methods. Using the common‐reflection‐surface approach locally in the unmigrated prestack data domain opens a potential for trace regularization and interpolation. In most data interpolation methods based on local coherency estimation, a single operator is designed for a target sample and the output amplitude is defined as a weighted average along the operator. This approach may fail in presence of interfering events or strong amplitude and phase variations. In this paper we introduce an alternative scheme in which there is no need for an operator to be defined at the target sample itself. Instead, the amplitude at a target sample is constructed from multiple operators estimated at different positions. In this case one operator may contribute to the construction of several target samples. Vice versa, a target sample might receive contributions from different operators. Operators are determined on a grid which can be sparser than the output grid. This allows to dramatically decrease the computational costs. In addition, the use of multiple operators for a single target sample stabilizes the interpolation results and implicitly allows several contributions in case of interfering events. Due to the considerable computational expense, common‐reflection‐surface interpolation is limited to work in subsets of the prestack data. We present the general workflow of a common‐reflection‐surface‐based regularization/interpolation for 3D data volumes. This workflow has been applied to an OBC common‐receiver volume and binned common‐offset subsets of a 3D marine data set. The impact of a common‐reflection‐surface regularization is demonstrated by means of a subsequent time migration. In comparison to the time migrations of the original and DMO‐interpolated data, the results show particular improvements in view of the continuity of reflections events. This gain is confirmed by an automatic picking of a horizon in the stacked time migrations.     

Fast modelling of marine surface multiples

In the development and testing of water‐surface multiple‐removal algorithms, it is valuable to have accurate synthetic seismograms which exhibit multiples, for which the multiple‐free solution is known. A method is presented for constructing 2D and 3D solutions of the acoustic wave equation in water, by combining the solution from a primary source with other scaled solutions of secondary sources, which simulate diffractors. The computation involves function evaluation rather than numerical solution of differential equations and is consequently accurate and comparatively fast. The analytic formulae on which the method is based give insights into methods for multiple removal. Generalized reflection coefficients, defined on a horizontal plane above the diffractors, are derived and used to construct the integral equations which are the basis for many multiple‐removal schemes.     

基于匹配滤波的多次波压制方法研究

常规滤波方法对浅海多次波压制不能取得理想效果.通过改变模型表面边界条件,求解波动方程实现多情景地震数值模拟.利用自由边界条件模拟时获得含海平面强反射多次波地震记录,而利用吸收边界条件时得到理想条件下无多次波地震记录.对比分析多次波对有效地震信号的影响得到海平面强反射多次波改变了有效地震信号的振幅、频率和相位.利用归一化方法计算得到压制多次波滤波器算子,并作用于含海平面强反射多次波地震记录.叠加地震剖面对比分析表明,归一化滤波器算子较好的压制了地震记录中的多次波,提高地震资料分辨率.最后,根据多次波产生条件,给出了利用实际区域简化地质模型的地震数值模拟结果计算归一化滤波器算子的方法,从而利用该方法压制实际地震数据中多次波.     

Rough seas and time-lapse seismic

Time‐lapse seismic surveying has become an accepted tool for reservoir monitoring applications, thus placing a high premium on data repeatability. One factor affecting data repeatability is the influence of the rough sea‐surface on the ghost reflection and the resulting seismic wavelets of the sources and receivers. During data analysis, the sea‐surface is normally assumed to be stationary and, indeed, to be flat. The non‐flatness of the sea‐surface introduces amplitude and phase perturbations to the source and receiver responses and these can affect the time‐lapse image. We simulated the influence of rough sea‐surfaces on seismic data acquisition. For a typical seismic line with a 48‐fold stack, a 2‐m significant‐wave‐height sea introduces RMS errors of about 5–10% into the stacked data. This level of error is probably not important for structural imaging but could be significant for time‐lapse surveying when the expected difference anomaly is small. The errors are distributed differently for sources and receivers because of the different ways they are towed. Furthermore, the source wavelet is determined by the sea shape at the moment the shot is fired, whereas the receiver wavelet is time‐varying because the sea moves significantly during the seismic record.     

Point‐spread functions for interferometric imaging

Interferometric redatuming is a data‐driven method to transform seismic responses with sources at one level and receivers at a deeper level into virtual reflection data with both sources and receivers at the deeper level. Although this method has traditionally been applied by cross‐correlation, accurate redatuming through a heterogeneous overburden requires solving a multidimensional deconvolution problem. Input data can be obtained either by direct observation (for instance in a horizontal borehole), by modelling or by a novel iterative scheme that is currently being developed. The output of interferometric redatuming can be used for imaging below the redatuming level, resulting in a so‐called interferometric image. Internal multiples from above the redatuming level are eliminated during this process. In the past, we introduced point‐spread functions for interferometric redatuming by cross‐correlation. These point‐spread functions quantify distortions in the redatumed data, caused by internal multiple reflections in the overburden. In this paper, we define point‐spread functions for interferometric imaging to quantify these distortions in the image domain. These point‐spread functions are similar to conventional resolution functions for seismic migration but they contain additional information on the internal multiples in the overburden and they are partly data‐driven. We show how these point‐spread functions can be visualized to diagnose image defocusing and artefacts. Finally, we illustrate how point‐spread functions can also be defined for interferometric imaging with passive noise sources in the subsurface or with simultaneous‐source acquisition at the surface.     

Water‐bottom multiple attenuation by Kirchhoff extrapolation

Despite being less general than 3D surface‐related multiple elimination (3D‐SRME), multiple prediction based on wavefield extrapolation can still be of interest, because it is less CPU and I/O demanding than 3D‐SRME and moreover it does not require any prior data regularization. Here we propose a fast implementation of water‐bottom multiple prediction that uses the Kirchhoff formulation of wavefield extrapolation. With wavefield extrapolation multiple prediction is usually obtained through the cascade of two extrapolation steps. Actually by applying the Fermat’s principle (i.e., minimum reflection traveltime) we show that the cascade of two operators can be replaced by a single approximated extrapolation step. The approximation holds as long as the water bottom is not too complex. Indeed the proposed approach has proved to work well on synthetic and field data when the water bottom is such that wavefront triplications are negligible, as happens in many practical situations.     

Scattering diagrams in seismic imaging: More insights into the construction of virtual events and internal multiples

The key processes in marine seismic imaging include (i) removing from seismic data all seismic events (free-surface multiples and ghosts) which contain at least one reflection at the sea surface in their wave-propagation path, and leaving those with no reflection at the free surface (internal multiples and primaries), (ii) removing events with at least two reflections in the subsurface (internal multiples), and leaving events with only one reflection in the subsurface (primaries), and then (iii) locating the scattering points and reflectors inside the subsurface which are the sources of primaries and internal multiple events. All these processes are here explained, derived, and optimized via scattering diagrams (diagrammatica) in a way similar to the way the quantum field theory is often explained via Feynman diagrams. Our discussion of the removal of events with free-surface reflections from the data will be brief, as the diagrammatica of these events are now well understood.The main focus of this paper is the diagrammatica of internal multiples and primaries. Although these events do not contain any reflection at the sea surface, it is important to reconstruct them with scattering points near the sea surface, where seismic data are recorded. So our diagrammatica of primaries and internal multiples include events which are not directly recorded in seismic data but which can be constructed from seismic data. These events have allowed us to construct scattering diagrams of primaries and internal multiples with scattering points near the sea surface. Furthermore, these new diagrammatica of internal multiples and primaries can be used to remove internal multiples from the data.     

Diffraction imaging by uniform asymptotic theory and double exponential fitting

Seismic diffracted waves carry valuable information for identifying geological discontinuities. Unfortunately, the diffraction energy is generally too weak, and standard seismic processing is biased to imaging reflection. In this paper, we present a dynamic diffraction imaging method with the aim of enhancing diffraction and increasing the signal‐to‐noise ratio. The correlation between diffraction amplitudes and their traveltimes generally exists in two forms, with one form based on the Kirchhoff integral formulation, and the other on the uniform asymptotic theory. However, the former will encounter singularities at geometrical shadow boundaries, and the latter requires the computation of a Fresnel integral. Therefore, neither of these methods is appropriate for practical applications. Noting the special form of the Fresnel integral, we propose a least‐squares fitting method based on double exponential functions to study the amplitude function of diffracted waves. The simple form of the fitting function has no singularities and can accelerate the calculation of diffraction amplitude weakening coefficients. By considering both the fitting weakening function and the polarity reversal property of the diffracted waves, we modify the conventional Kirchhoff imaging conditions and formulate a diffraction imaging formula. The mechanism of the proposed diffraction imaging procedure is based on the edge diffractor, instead of the idealized point diffractor. The polarity reversal property can eliminate the background of strong reflection and enhance the diffraction by same‐phase summation. Moreover,the fitting weakening function of diffraction amplitudes behaves like an inherent window to optimize the diffraction imaging aperture by its decaying trend. Synthetic and field data examples reveal that the proposed diffraction imaging method can meet the requirement of high‐resolution imaging, with the edge diffraction fully reinforced and the strong reflection mostly eliminated.     

复杂地表边界元-体积元波动方程数值模拟

复杂近地表引起来自深部构造的地震反射信号振幅和相位的异常变化,是影响复杂近地表地区地震资料品质的主要原因.本文采用边界元-体积元方法,通过求解含复杂地表的波动积分方程,来模拟地震波在复杂近地表构造中的传播.其中,边界元法模拟地形起伏和表层地质结构对地震波传播的影响;体积元法模拟起伏地表下非均质低降速层的影响.与其他数值...     

Internal multiple attenuation by Kirchhoff extrapolation

Surface‐related multiple elimination is the leading methodology for surface multiple removal. This data‐driven approach can be extended to interbed multiple prediction at the expense of a huge increase of the computational burden. This cost makes model‐driven methods still attractive, especially for the three dimensional case. In this paper we present a methodology that extends Kirchhoff wavefield extrapolation to interbed multiple prediction. In Kirchhoff wavefield extrapolation for surface multiple prediction a single round trip to an interpreted reflector is added to the recorded data. Here we show that interbed multiples generated between two interpreted reflectors can be predicted by applying the Kirchhoff wavefield extrapolation operator twice. In the first extrapolation step Kirchhoff wavefield extrapolation propagates the data backward in time to simulate a round trip to the shallower reflector. In the second extrapolation step Kirchhoff wavefield extrapolation propagates the data forward in time to simulate a round trip to the deeper reflector. In the Kirchhoff extrapolation kernel we use asymptotic Green's functions. The prediction of multiples via Kirchhoff wavefield extrapolation is possibly sped up by computing the required traveltimes via a shifted hyperbola approximation. The effectiveness of the method is demonstrated by results on both synthetic and field data sets.     

Correlation‐based reflection waveform inversion by one‐way wave equations

Reflection full waveform inversion can update subsurface velocity structure of the deeper part, but tends to get stuck in the local minima associated with the waveform misfit function. These local minima cause cycle skipping if the initial background velocity model is far from the true model. Since conventional reflection full waveform inversion using two‐way wave equation in time domain is computationally expensive and consumes a large amount of memory, we implement a correlation‐based reflection waveform inversion using one‐way wave equations to retrieve the background velocity. In this method, one‐way wave equations are used for the seismic wave forward modelling, migration/de‐migration and the gradient computation of objective function in frequency domain. Compared with the method using two‐way wave equation, the proposed method benefits from the lower computational cost of one‐way wave equations without significant accuracy reduction in the cases without steep dips. It also largely reduces the memory requirement by an order of magnitude than implementation using two‐way wave equation both for two‐ and three‐dimensional situations. Through numerical analysis, we also find that one‐way wave equations can better construct the low wavenumber reflection wavepath without producing high‐amplitude short‐wavelength components near the image points in the reflection full waveform inversion gradient. Synthetic test and real data application show that the proposed method efficiently updates the background velocity model.     

How rough sea affects marine seismic data and deghosting procedures

Most seismic processing algorithms generally consider the sea surface as a flat reflector. However, acquisition of marine seismic data often takes place in weather conditions where this approximation is inaccurate. The distortion in the seismic wavelet introduced by the rough sea may influence (for example) deghosting results, as deghosting operators are typically recursive and sensitive to the changes in the seismic signal. In this paper, we study the effect of sea surface roughness on conventional (5–160 Hz) and ultra‐high‐resolution (200–3500 Hz) single‐component towed‐streamer data. To this end, we numerically simulate reflections from a rough sea surface using the Kirchhoff approximation. Our modelling demonstrates that for conventional seismic frequency band sea roughness can distort results of standard one‐dimensional and two‐dimensional deterministic deghosting. To mitigate this effect, we introduce regularisation and optimisation based on the minimum‐energy criterion and show that this improves the processing output significantly. Analysis of ultra‐high‐resolution field data in conjunction with modelling shows that even relatively calm sea state (i.e., 15 cm wave height) introduces significant changes in the seismic signal for ultra‐high‐frequency band. These changes in amplitude and arrival time may degrade the results of deghosting. Using the field dataset, we show how the minimum‐energy optimisation of deghosting parameters improves the processing result.     

Insight into the marine controlled‐source electromagnetic signal propagation

In marine controlled‐source electromagnetic (CSEM) surveys the subsurface is explored by emitting low‐frequency signals from an electric dipole source close to the sea‐bed. The main goal is often to detect and describe possible thin resistive layers beneath the sea‐bed. To gain insight into how CSEM signals propagate, it is informative to study a stratified model. The electromagnetic field is then given in terms of integrals over TE‐ and TM‐polarized plane‐wave constituents. An asymptotic evaluation of the field integrals for large propagation distances results in explicit spatial expressions for the field components and the derived expressions can be used to analyse how the CSEM signals propagate. There are two major signal pathways in a standard CSEM model. One of these pathways is via the thin resistive layer and the resulting response is accounted for by a pole in the reflection response for the TM mode. The signal is propagating nearly vertically down to the resistor from the source, then guided while attenuated along the reservoir, before propagating nearly vertically up to the receiver. The response is slightly altered by the sea‐bed interface and further modified in shallow water due to multiple reflections between the sea‐surface and sea‐bed at both the source and receiver sides. The other major signal pathway is via the resistive air half‐space, the so‐called airwave. The airwave is generated by the TE mode and interacts with the subsurface via vertically propagating signals reflected between the sea‐surface and subsurface at both the source and receiver sides.     

The common focal point and common reflection surface methodologies: subsets or complements?

The common focal point (CFP) method and the common reflection surface (CRS) stack method are compared. The CRS method is a fast, highly automated procedure that provides high S/N ratio simulation of zero‐offset (ZO) images by combining, per image point, the reflection energy of an arc segment that is tangential to the reflector. It uses smooth parametrized two‐way stacking operators, based on a data‐driven triplet of attributes in 2D (eight parameters in 3D). As a spin‐off, the attributes can be used for several applications, such as the determination of the geometrical spreading factor, multiple prediction, and tomographic inversion into a smooth background velocity model. The CFP method aims at decomposing two‐way seismic reflection data into two full‐aperture one‐way propagation operators. By applying an iterative updating procedure in a half‐migrated domain, it provides non‐smooth focusing operators for prestack imaging using only the energy from one focal point at the reflector. The data‐driven operators inhibit all propagation effects of the overburden. The CFP method provides several spin‐offs, amongst which is the CFP matrix related to one focal point, which displays the reflection amplitudes as measured at the surface for each source–receiver pair. The CFP matrix can be used to determine the specular reflection source–receiver pairs and the Fresnel zone at the surface for reflection in one single focal point. Other spin‐offs are the prediction of internal multiples, the determination of reflectivity effects, velocity‐independent redatuming and tomographic inversion to obtain a velocity–depth model. The CFP method is less fast and less automated than the CRS method. From a pointwise comparison of features it is concluded that one method is not a subset of the other, but that both methods can be regarded as being to some extent complementary.     

A Tailored Computation of the Mean Dynamic Topography for a Consistent Integration into Ocean Circulation Models

Geostrophic surface velocities can be derived from the gradients of the mean dynamic topography—the difference between the mean sea surface and the geoid. Therefore, independently observed mean dynamic topography data are valuable input parameters and constraints for ocean circulation models. For a successful fit to observational dynamic topography data, not only the mean dynamic topography on the particular ocean model grid is required, but also information about its inverse covariance matrix. The calculation of the mean dynamic topography from satellite-based gravity field models and altimetric sea surface height measurements, however, is not straightforward. For this purpose, we previously developed an integrated approach to combining these two different observation groups in a consistent way without using the common filter approaches (Becker et al. in J Geodyn 59(60):99–110, 2012; Becker in Konsistente Kombination von Schwerefeld, Altimetrie und hydrographischen Daten zur Modellierung der dynamischen Ozeantopographie 2012). Within this combination method, the full spectral range of the observations is considered. Further, it allows the direct determination of the normal equations (i.e., the inverse of the error covariance matrix) of the mean dynamic topography on arbitrary grids, which is one of the requirements for ocean data assimilation. In this paper, we report progress through selection and improved processing of altimetric data sets. We focus on the preprocessing steps of along-track altimetry data from Jason-1 and Envisat to obtain a mean sea surface profile. During this procedure, a rigorous variance propagation is accomplished, so that, for the first time, the full covariance matrix of the mean sea surface is available. The combination of the mean profile and a combined GRACE/GOCE gravity field model yields a mean dynamic topography model for the North Atlantic Ocean that is characterized by a defined set of assumptions. We show that including the geodetically derived mean dynamic topography with the full error structure in a 3D stationary inverse ocean model improves modeled oceanographic features over previous estimates.     

FDM Simulation of Seismic Waves, Ocean Acoustic Waves, and Tsunamis Based on Tsunami-Coupled Equations of Motion

We have developed a new, unified modeling technique for the total simulation of seismic waves, ocean acoustic waves, and tsunamis resulting from earthquakes, based on a finite difference method simulation of the 3D equations of motion. Using the equilibrium between the pressure gradient and gravity in these equations, tsunami propagation is naturally incorporated in the simulation based on the equations of motion. The performance of the parallel computation for the newly developed tsunami-coupled equations using a domain partitioning procedure shows a high efficiency coefficient with a large number of CPU cores. The simulation results show how the near-field term associated with seismic waves produced by shallow earthquakes leads to a permanent coseismic deformation of the ground surface, which gives rise to the initial tsunami on the sea surface. Propagation of the tsunami along the sea surface as a gravity wave, and ocean acoustic waves in seawater with high-frequency multiple P-wave reflections between the free surface and sea bottom, are also clearly demonstrated by the present simulations. We find a good agreement in the tsunami waveform between our results and those obtained by other simulations based on an analytical model and the Navier–Stokes equations, demonstrating the effectiveness of the tsunami-coupling simulation model. Based on this simulation, we show that the ratio of the amplitude of ocean acoustic waves to the height of the tsunami, both of which are produced by the earthquake, strongly depends on the rise time of the earthquake rupture. This ratio can be used to obtain a more detailed understanding of the source rupture processes of subduction zone earthquakes, and for implementing an improved tsunami alert system for slow tsunami earthquakes.     

Estimation and application of near‐surface full waveform redatuming operators

Land seismic data quality can be severely affected by near‐surface anomalies. The imprint of a complex near‐surface can be removed by redatuming the data to a level below the surface, from where the subsurface structures are assumed to be relatively smooth. However, to derive a velocity‐depth model that explains the propagation effects of the near‐surface is a non‐trivial task. Therefore, an alternative approach has been proposed, where the redatuming operators are obtained in a data‐driven manner from the reflection event related to the datum. In the current implementation, the estimation of these redatuming operators is done in terms of traveltimes only, based on a high‐frequency approximation. The accompanying amplitudes are usually derived from a local homogeneous medium, which is obviously a simplification of reality. Such parametrization has produced encouraging results in the past but cannot completely remove the near‐surface complexities, leaving artefacts in the redatumed results. In this paper we propose a method that estimates the redatuming operators directly from the data, i.e., without using a velocity model, in a full waveform manner, such that detailed amplitude and phase variations are included. The method directly outputs the inverse propagation operators that are needed for true‐amplitude redatuming. Based on 2D synthetic data it is demonstrated that the resulting redatuming quality is improved and artefacts are reduced.     

Retrieving virtual reflection responses at drill‐bit positions using seismic interferometry with drill‐bit noise

In the field of seismic interferometry, researchers have retrieved surface waves and body waves by cross‐correlating recordings of uncorrelated noise sources to extract useful subsurface information. The retrieved wavefields in most applications are between receivers. When the positions of the noise sources are known, inter‐source interferometry can be applied to retrieve the wavefields between sources, thus turning sources into virtual receivers. Previous applications of this form of interferometry assume impulsive point sources or transient sources with similar signatures. We investigate the requirements of applying inter‐source seismic interferometry using non‐transient noise sources with known positions to retrieve reflection responses at those positions and show the results using synthetic drilling noise as source. We show that, if pilot signals (estimates of the drill‐bit signals) are not available, it is required that the drill‐bit signals are the same and that the phases of the virtual reflections at drill‐bit positions can be retrieved by deconvolution interferometry or by cross‐coherence interferometry. Further, for this case, classic interferometry by cross‐correlation can be used if the source power spectrum can be estimated. If pilot signals are available, virtual reflection responses can be obtained by first using standard seismic‐while‐drilling processing techniques such as pilot cross‐correlation and pilot deconvolution to remove the drill‐bit signatures in the data and then applying cross‐correlation interferometry. Therefore, provided that pilot signals are reliable, drill‐bit data can be redatumed from surface to borehole depths using this inter‐source interferometry approach without any velocity information of the medium, and we show that a well‐positioned image below the borehole can be obtained using interferometrically redatumed reflection responses with just a simple velocity model. We discuss some of the practical hurdles that restrict the application of the proposed method offshore.     

Rough seas and statistical deconvolution

The rough‐sea reflection‐response varies (1) along the streamer (2) from shot to shot and (3) with time along the seismic trace. The resulting error in seismic data can be important for time‐lapse imaging. One potential way of reducing the rough‐sea receiver error is to use conventional statistical deconvolution, but special care is needed in the choice of the design and application windows. The well‐known deconvolution problem associated with the non‐whiteness of the reflection series is exacerbated by the requirement of an unusually short design window – a requirement that is imposed by the non‐stationary nature of the rough‐sea receiver wavelet. For a synthetic rough‐sea data set, with a white 1D reflection series, the design window needs to be about 1000 ms long, with an application window about 400 ms long, centred within the design window. Although such a short design window allows the deconvolution operator to follow the time‐variation of the rough‐sea wavelet, it is likely to be too short to prevent the non‐whiteness of the geology from corrupting the operator when it is used on real data. If finely spatial‐sampled traces are available from the streamer, the design window can be extended to neighbouring traces, making use of the spatial correlations of the rough‐sea wavelet. For this ‘wave‐following’ approach to be fruitful, the wind (and hence the dominant wave direction) needs to be roughly along the line of the streamer.