WAVE FIELD EXTRAPOLATION TECHNIQUES FOR INHOMOGENEOUS MEDIA WHICH INCLUDE CRITICAL ANGLE EVENTS. PART I: METHODS USING THE ONE-WAY WAVE EQUATIONS*

Part I of this series starts with a brief review of the fundamental principles underlying wave field extrapolation. Next, the total wave field is split into downgoing and upgoing waves, described by a set of coupled one-way wave equations. In cases of limited propagation angles and weak inhomogeneities these one-way wave equations can be decoupled, describing primary waves only. For large propagation angles (up to and including 90°) an alternative choice of sub-division into downgoing and upgoing waves is presented. It is shown that this approach is well suited for modeling as well as migration and inversion schemes for seismic data which include critical angle events.     

THE SPECTRAL FUNCTION OF A LAYERED SYSTEM AND THE DETERMINATION OF THE WAVEFORMS AT DEPTH*

Consider the mathematical model of a horizontally layered system subject to an initial downgoing source pulse in the upper layer and to the condition that no upgoing waveforms enter the layered system from below the deepest interface. The downgoing waveform (as measured from its first arrival) in each layer is necessarily minimum-phase. The net downgoing energy in any layer, defined as the difference of the energy spectrum of the downgoing wave minus the energy spectrum of the upgoing wave, is itself in the form of an energy spectrum, that is, it is non-negative for all frequencies. The z-transform of the autocorrelation function corresponding to the net downgoing energy spectrum is called the net downgoing spectral function for the layer in question. The net downgoing spectral functions of any two layers A and B are related as follows: the product of the net downgoing spectral function of layer A times the overall transmission coefficient from A to B equals the product of the net downgoing spectral function of layer B times the overall transmission coefficient from B to A. The net downgoing spectral function for the upper layer is called simply the spectral function of the system. In the case of a marine seismogram, the autocorrelation function corresponding to the spectral function can be used to recursively generate prediction error operators of successively increasing lengths, and at the same time the reflection coefficients at successively increasing depths. This recursive method is mathematically equivalent to that used in solving the normal equations in the case of Toeplitz forms. The upgoing wave-form in any given layer multiplied by the direct transmission coefficient from that layer to the surface is equal to the convolution of the corresponding prediction error operator with the surface seismogram. The downgoing waveform in this given layer multiplied by the direct transmission coefficient from that layer to the surface is equal to the convolution of the corresponding hindsight error operator (i.e., the time reverse of the prediction error operator) with the surface seismogram.     

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.     

DECOMPOSITION OF MULTICOMPONENT SEISMIC DATA INTO PRIMARY P- AND S-WAVE RESPONSES1

Inversion of multicomponent seismic data can be subdivided in three main processes: (1) Surface-related preprocessing (decomposition of the multicomponent data into ‘primary’ P-and S-wave responses). (2) Prestack migration of the primary P- and S-wave responses, yielding the (angle-dependent) P-P, P-S, S-P and S-S reflectivity of the subsurface. (3) Target-related post-processing (transformation of the reflectivity into the rock and pore parameters in the target). This paper deals with the theoretical aspects of surface-related preprocessing. In a multicomponent seismic data set the P- and S-wave responses of the subsurface are distorted by two main causes: (1) The seismic vibrators always radiate a mixture of P- and S-waves into the subsurface. Similarly, the geophones always measure a mixture of P- and S-waves. (2) The free surface reflects any upgoing wave fully back into the subsurface. This gives rise to strong multiple reflections, including conversions. Therefore, surface-related preprocessing consists of two steps: (1)Decomposition of the multicomponent data (pseudo P- and S-wave responses) into true P- and S-wave responses. In practice this procedure involves (a) decomposition per common shot record of the particle velocity vector into scalar upgoing P- and S-waves, followed by (b) decomposition per common receiver record of the traction vector into scalar downgoing P- and S-waves. (2) Elimination of the surface-related multiple reflections and conversions. In this procedure the free surface is replaced by a reflection-free surface. The effect is that we obtain ‘primary’ P-and S-wave responses, that contain internal multiples only. An interesting aspect of the procedure is that no knowledge of the subsurface is required. In fact, the subsurface may have any degree of complexity. Both the decomposition step and the multiple elimination step are fully determined by the medium parameters at the free surface only. After surface-related preprocessing, the scalar P- and S-wave responses can be further processed independently by existing scalar algorithms.     

Elastic reverse-time migration based on amplitude-preserving P- and S-wave separation

Imaging the PP- and PS-wave for the elastic vector wave reverse-time migration requires separating the P- and S-waves during the wave field extrapolation. The amplitude and phase of the P- and S-waves are distorted when divergence and curl operators are used to separate the P- and S-waves. We present a P- and S-wave amplitude-preserving separation algorithm for the elastic wavefield extrapolation. First, we add the P-wave pressure and P-wave vibration velocity equation to the conventional elastic wave equation to decompose the P- and S-wave vectors. Then, we synthesize the scalar P- and S-wave from the vector Pand S-wave to obtain the scalar P- and S-wave. The amplitude-preserved separated P- and S-waves are imaged based on the vector wave reverse-time migration (RTM). This method ensures that the amplitude and phase of the separated P- and S-wave remain unchanged compared with the divergence and curl operators. In addition, after decomposition, the P-wave pressure and vibration velocity can be used to suppress the interlayer reflection noise and to correct the S-wave polarity. This improves the image quality of P- and S-wave in multicomponent seismic data and the true-amplitude elastic reverse time migration used in prestack inversion.     

Joint viscoacoustic waveform inversion with the separated upgoing and downgoing wavefields in zero-offset vertical seismic profile data

The attenuation of seismic waves propagating in reservoirs can be obtained accurately from the data analysis of vertical seismic profile in terms of the quality-factor Q. The common methods usually use the downgoing wavefields in vertical seismic profile data. However, the downgoing wavefields consist of more than 90% energy of the spectrum of the vertical seismic profile data, making it difficult to estimate the viscoacoustic parameters accurately. Thus, a joint viscoacoustic waveform inversion of velocity and quality-factor is proposed based on the multi-objective functions and analysis of the difference between the results inverted from the separated upgoing and downgoing wavefields. A simple separating step is accomplished by the reflectivity method to obtain the individual wavefields in vertical seismic profile data, and then a joint inversion is carried out to make full use of the information of the individual wavefields and improve the convergence of viscoacoustic full-waveform inversion. The sensitivity analysis of the different wavefields to the velocity and quality-factor shows that the upgoing and downgoing wavefields contribute differently to the viscoacoustic parameters. A numerical example validates our method can improve the accuracy of viscoacoustic parameters compared with the direct inversion using full wavefield and the separate inversion using upgoing or downgoing wavefield. The application on real field data indicates our method can recover a reliable viscoacoustic model, which helps reservoir appraisal.     

Q-LOG DETERMINATION ON DOWNGOING WAVELETS AND TUBE WAVE ANALYSIS IN VERTICAL SEISMIC PROFILES1

Seismic attenuation introduces modifications in the wavelet shape in vertical seismic profiles. These modifications can be quantified by measuring particular signal attributes such as rise-time, period and shape index. Use of signal attributes leads to estimations of a seismic-attenuation log (Q-log). To obtain accurate signal attributes it is important to minimize noise influence and eliminate local interference between upgoing and downgoing waves at each probe location. When tube waves are present it is necessary to eliminate them before performing separation of upgoing and downgoing events. We used a trace-by-trace Wiener filter to minimize the influence of tube waves. The separation of upgoing and downgoing waves was then performed in the frequency domain using a trace-pair filter. We used three possible methods based on signal attribute measurements to obtain g-log from the extracted downgoing wavefield. The first one uses a minimum phasing filter and the arrival time of the first extremum. The two other methods determine the Q-factor from simple relations between the amplitudes of the first extrema and the pseudo-periods of the down-going wavelet. The relations determined between a signal attribute and traveltime over quality factor were then calibrated using field source signature and constant-Q models computed by Ganley's method. Q-logs thus obtained from real data are discussed and compared with geological information, specifically at reservoir level. Analysis of the tube wave arrivals at the level of the reservoir showed a tube wave attenuation that could not be explained by simple transmission effects. There was also a loss of signal coherence. This could be interpreted as tube wave diffusion in the porous reservoir, followed by dispersion. If this interpretation can be verified, tube wave analysis could lead to further characterization of porous permeable zones.     

陆上纵波源VSP资料中的纯横波

常规陆上VSP(Vertical Seismic Profiling)勘探普遍采用纵波震源激发,三分量检波器接收,主要利用的是纵波和转换横波信息。已有的研究表明,炸药震源在井下激发、可控震源在地面垂向振动,均会产生较强的纯纵波和一定强度的纯横波;泊松比差别较大的分界面有利于形成较强的透射转换横波。本文通过对激发形成的纯横波和下行转换形成的横波进行对比分析,认为纯横波的主频往往低于纯纵波的主频,而下行转换横波的主频通常接近纵波的主频。本文分别对两个陆上纵波源零偏和非零偏VSP资料进行分析,结果表明这些资料中普遍存在纯横波,只是横波的强弱存在不同程度的变化。利用纵波源零偏VSP资料,可以获得横波速度信。最后对VSP纵波和横波联合应用前景进行了分析,应该充分利用纵波源VSP资料中的横波信息。     

Topological invariants of field lines rooted to planes

Abstract

Two open curves with fixed endpoints on a boundary surface can be topologically linked. However, the Gauss linkage integral applies only to closed curves and cannot measure their linkage. Here we employ the concept of relative helicity in order to define a linkage for open curves. For a magnetic field consisting of closed field lines, the magnetic helicity integral can be expressed as the sum of Gauss linkage integrals over pairs of lines. Relative helicity extends the helicity integral to volumes where field lines may cross the boundary surface. By analogy, linkages can be defined for open lines by requiring that their sum equal the relative helicity.

With this definition, the linkage of two lines which extend between two parallel planes simply equals the number of turns the lines take about each other. We obtain this result by first defining a gauge-invariant, one-dimensional helicity density, i.e. the relative helicity of an infinitesimally thin plane slab. This quantity has a physical interpretation in terms of the rate at which field lines lines wind about each other in the direction normal to the plane. A different method is employed for lines with both endpoints on one plane; this method expresses linkages in terms of a certain Gauss linkage integral plus a correction term. In general, the linkage number of two curves can be put in the form L=r + n, |r|≦1J2, where r depends only on the positions of the endpoints, and n is an integer which reflects the order of braiding of the curves.

Given fixed endpoints, the linkage numbers of a magnetic field are ideal magneto-hydrodynamic invariants. These numbers may be useful in the analysis of magnetic structures not bounded by magnetic surfaces, for example solar coronal fields rooted in the photosphere. Unfortunately, the set of linkage numbers for a field does not uniquely determine the field line topology. We briefly discuss the problem of providing a complete and economical classification of field topologies, using concepts from the theory of braid equivalence classes.     

Seismic migration and absorbing boundaries with a one-way wave system for heterogeneous media1

A first-order one-way wave system has been created based on characteristic analysis of the acoustic wave system and optimization of the dispersion relation. We demonstrate that this system is equivalent to a third-order scalar partial-differential equation which, for a homogeneous medium, reduces to a form similar to the 45° paraxial wave equation. This system describes accurately waves propagating in a 2D heterogeneous medium at angles up to 75°. The one-way wave system representing downgoing waves is used for a modified reverse time migration method. As a wavefield extrapolator in migration, the downgoing wave system propagates the reflection events backwards to their reflectors without scattering at the discontinuities in the velocity model. Hence, images with amplitudes proportional to reflectivity can be obtained from this migration technique. We present examples of the application of the new migration method to synthetic seismic data where P-P reflections P-SV converted waves are present. Absorbing boundaries, useful in the generation of synthetic seismograms, have been constructed by using the one-way wave system. These boundaries absorb effectively waves impinging over a wide range of angles of incidence.     

Theory of convective saturation of Langmuir waves during ionospheric modification of a barium cloud

In recent experiments (Djuth, F. T., Sulzer, M. P., Elder, J. H. and Groves, K. M. (1995) Journal of Geophysical Research, 100, 17,347), a parametric decay instability was excited by an ordinarywave HF pump during an ionospheric chemical release from a rocket over Arecibo, PR, which created an artificial ‘barium ionosphere,’ with peak plasma frequency above the pump frequency, and a density gradient with a (short) 5 km scale length. Simultaneous incoherent scattering measurements revealed a strong initial asymmetry in the amplitudes of almost vertically upgoing versus downgoing measured plasma waves. We can account for this asymmetry in terms of linear convective saturation of parametrically unstable plasma waves propagating over a range of altitudes along geometric optics ray paths. Qualitative features of the frequency spectrum of the measured downgoing wave are in agreement with this model, although the theoretically predicted spectrum is narrower than observed. The observed altitude localization of the enhanced spectrum to a few range cells is consistent with the theory.     

Marine seismic wavefield measurement to remove sea-surface multiples

We propose a new method for removing sea-surface multiples from marine seismic reflection data in which, in essence, the reflection response of the earth, referred to a plane just above the sea-floor, is computed as the ratio of the plane-wave components of the upgoing wave and the downgoing wave. Using source measurements of the wavefield made during data acquisition, three problems associated with earlier work are solved: (i) the method accommodates source arrays, rather than point sources; (ii) the incident field is removed without simultaneously removing part of the scattered field; and (iii) the minimum-energy criterion to find a wavelet is eliminated. Pressure measurements are made in a horizontal plane in the water. The source can be a conventional array of airguns, but must have both in-line and cross-line symmetry, and its wavefield must be measured and be repeatable from shot to shot. The problem is formulated for multiple shots in a two-dimensional configuration for each receiver, and for multiple receivers in a two-dimensional configuration for each shot. The scattered field is obtained from the measurements by subtracting the incident field, known from measurements at the source. The scattered field response to a single incident plane wave at a single receiver is obtained by transforming the common-receiver gather to the frequency–wavenumber domain, and a single component of this response is obtained by Fourier transforming over all receiver coordinates. Each scattered field component is separated into an upgoing wave and a downgoing wave using the zero-pressure condition at the water-surface. The upgoing wave may then be expressed as a reflection coefficient multiplied by the incident downgoing wave plus a sum of scattered downgoing plane waves, each multiplied by the corresponding reflection coefficient. Keeping the upgoing scattered wave fixed, and using all possible incident plane waves for a given frequency, yields a set of linear simultaneous equations for the reflection coefficients which are solved for each plane wave and for each frequency. To create the shot records that would have been measured if the sea-surface had been absent, each reflection coefficient is multiplied by complex amplitude and phase factors, for source and receiver terms, before the five-dimensional Fourier transformation back to the space–time domain.     

Characterization of seismic hazard and structural response by energy flux

Seismic safety of structures depends on the structure's ability to absorb the seismic energy that is transmitted from ground to structure. One parameter that can be used to characterize seismic energy is the energy flux. Energy flux is defined as the amount of energy transmitted per unit time through a cross-section of a medium, and is equal to kinetic energy multiplied by the propagation velocity of seismic waves. The peak or the integral of energy flux can be used to characterize ground motions. By definition, energy flux automatically accounts for site amplification. Energy flux in a structure can be studied by formulating the problem as a wave propagation problem. For buildings founded on layered soil media and subjected to vertically incident plane shear waves, energy flux equations are derived by modeling the building as an extension of the layered soil medium, and considering each story as another layer. The propagation of energy flux in the layers is described in terms of the upgoing and downgoing energy flux in each layer, and the energy reflection and transmission coefficients at each interface. The formulation results in a pair of simple finite-difference equations for each layer, which can be solved recursively starting from the bedrock. The upgoing and downgoing energy flux in the layers allows calculation of the energy demand and energy dissipation in each layer. The methodology is applicable to linear, as well as nonlinear structures.     

VERTICAL SEISMIC PROFILING IN COAL*

The intellection of seismic wave propagation in coal measures demands direct observation of the wavefield progression. Two vertical seismic profiles with high spatial and temporal sampling, were recently recorded in the Sydney Basin coalfields as part of an experimental coal seismic program. Static corrections and interval velocities were obtained by an automated system to determine first kicks and pulse rise times. Upgoing and downgoing waves were separated in the f—k-plane using a novel technique of contour slice filtering. The isolated upgoing waves clearly display reflections from the major coal seams within the stratigraphic sequence. The downgoing wave spectra were subjected to attenuation analysis. The deduced specific quality factor Q for Permian coal measure rocks lies in the range 20–70. Similar estimates were obtained in the time domain from measurements of pulse broadening. Synthetic VSP seismograms, computed using an exact recursive formulation, are an indispensable aid to interpretation. They illustrate the filtering effects of coal seams and sequences, and the effects of the contribution of internal and free-surface multiple reflections in the recorded wavetrains.     

频率-波数域递推计算海底多分量地震记录中的反射系数

在频率-波数域中采用解析法,解出多层条件下海底实测的多分量地震数据分解成上行和下行P波和S波的算法,导出海底各层地震反射系数随入射角变化(简称RVA)的递推计算公式,为海底多波多分量AVO弹性参数的反演及流体因子预测提供基础数据.合成数据的计算结果表明,本文给出的算法能较可靠地从海底多波多分量记录中提取RVA信息.     

DEPTH IMAGING OF OFFSET VERTICAL SEISMIC PROFILE DATA1

Depth migration consists of two different steps: wavefield extrapolation and imaging. The wave propagation is firmly founded on a mathematical frame-work, and is simulated by solving different types of wave equations, dependent on the physical model under investigation. In contrast, the imaging part of migration is usually based on ad hoc‘principles’, rather than on a physical model with an associated mathematical expression. The imaging is usually performed using the U/D concept of Claerbout (1971), which states that reflectors exist at points in the subsurface where the first arrival of the downgoing wave is time-coincident with the upgoing wave. Inversion can, as with migration, be divided into the two steps of wavefield extrapolation and imaging. In contrast to the imaging principle in migration, imaging in inversion follows from the mathematical formulation of the problem. The image with respect to the bulk modulus (or velocity) perturbations is proportional to the correlation between the time derivatives of a forward-propagated field and a backward-propagated residual field (Lailly 1984; Tarantola 1984). We assume a physical model in which the wave propagation is governed by the 2D acoustic wave equation. The wave equation is solved numerically using an efficient finite-difference scheme, making simulations in realistically sized models feasible. The two imaging concepts of migration and inversion are tested and compared in depth imaging from a synthetic offset vertical seismic profile section. In order to test the velocity sensitivity of the algorithms, two erroneous input velocity models are tested. We find that the algorithm founded on inverse theory is less sensitive to velocity errors than depth migration using the more ad hoc U/D imaging principle.     

On the boundary integral equations in the theory of elastic vibrations

Forward seismic problems are solved for elastic media by rigorous methods (i.e., methods with controllable accuracy). Analysis of the current state of research on this subject suggests that the most promising methods are based on integral and integro-differential equations, notwithstanding the rather modest results of their application to solving forward problems in the theory of elastic vibrations. The second Green integral theorem for seismic waves, formulated and proven in the paper, yields a system of two boundary (surface) integral equations for the displacement vector u(M ₀) and the normal (to the boundary surface) vector component of the stress tensor t_n(M ₀). The integrands of the surface integrals in terms of which the function t_n(M ₀) is expressed on both sides of the interface between the medium and the heterogeneity contain the second derivatives of the Green’s tensor functions ? ^e (M ₀, M) and ? ⁱ (M ₀, M), respectively, which are responsible for a cubic singularity (third-order singularity) if the integration point M coincides with the observation point M ₀. An original method of eliminating the cubic singularity proposed in the paper involves special tensor normalization of the integrals on the outer and inner sides of the interface and subsequent subtraction of one integral from another in order to construct the second integral equation.     

Quadratic normal moveouts of symmetric reflections in elastic media: A quick tutorial

The design of reflection traveltime approximations for optimal stacking and inversion has always been a subject of much interest in seismic processing. A most prominent role is played by quadratic normal moveouts, namely reflection traveltimes around zero-offset computed as second-order Taylor expansions in midpoint and offset coordinates. Quadratic normal moveouts are best employed to model symmetric reflections, for which the ray code in the downgoing direction coincides with the ray code in the upgoing direction in reverse order. Besides pure (non-converted) primaries, many multiply reflected and converted waves give rise to symmetric reflections. We show that the quadratic normal moveout of a symmetric reflection admits a natural decomposition into a midpoint term and an offset term. These, in turn, can be be formulated as the traveltimes of the one-way normal (N) and normal-incidence-point (NIP) waves, respectively. With the help of this decomposition, which is valid for propagation in isotropic and anisotropic elastic media, we are able to derive, in a simple and didactic way, a unified expression for the quadratic normal moveout of a symmetric reflection in its most general form in 3D. The obtained expression allows for a direct interpretation of its various terms and fully encompasses the effects of velocity gradients and Earth surface topography.     

Analysis of quality factors for Rayleigh channel waves

To facilitate investigation of the effect of imperfect elastic dissipation on thepropagation of Rayleigh-type channel waves and use of their quality factors in investigationsof the properties of coal seams, a simple method for calculating the quality factor QR isproposed in this paper. Introduction of complex velocities into the dispersion function allowscalculation of the dispersion function of Rayleigh-type channel waves in coal seams. By thecontrol variable method, we analyzed changes in QR with changes in coal seam thickness andP- and S-wave Q-factors within the coal seam and adjacent rock layers. The numerical resultsshow that the trend of the QR curve is consistent with the group velocity curve. The minimumQR value occurs at the Airy phase frequency; the Airy phase frequency decreases as coal seamthickness increases. The value of QR increases with increasing Qs2 （quality factor for S wavein coal seam）. We can compensate for the absorption of Rayleigh-type channel waves usingthe computed QR curve. Inversion of the QR curve can also be used to predict the thicknessesand litholoeies of coal seams.