ELASTIC EXTRAPOLATION OF PRIMARY SEISMIC P- AND S-WAVES1

The elastic Kirchhoff-Helmholtz integral expresses the components of the monochromatic displacement vector at any point A in terms of the displacement field and the stress field at any closed surface surrounding A. By introducing Green's functions for P- and S-waves, the elastic Kirchhoff-Helmholtz integral is modified such that it expresses either the P-wave or the S-wave at A in terms of the elastic wavefield at the closed surface. This modified elastic Kirchhoff-Helmholtz integral is transformed into one-way elastic Rayleigh-type integrals for forward extrapolation of downgoing and upgoing P- and S-waves. We also derive one-way elastic Rayleigh-type integrals for inverse extrapolation of downgoing and upgoing P- and S-waves. The one-way elastic extrapolation operators derived in this paper are the basis for a new prestack migration scheme for elastic data.     

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.     

THE GENERALIZED SPECTRAL FUNCTION OF A PLANE LAYERED MEDIUM*

The spectral function of a plane layered medium, which represents the net downgoing energy in the first layer due to a normally incident impulsive plane wave, plays an important role in the solution of the one-dimensional inverse problem in reflection seismology. Hitherto the extension to non-normal incidence was known only for a medium with free surface. By giving the extension for arbitrary surface reflection coefficients, this paper fills a gap.     

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.     

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.     

基于改进线性Radon变换的井孔地震上下行波场分离方法

反射波场分离是井孔地震资料处理中极其重要的一个环节,波场分离的质量直接影响成像结果的精度.不管是VSP还是井间地震资料,其反射波时距曲线都近似直线型,根据这一特征,本文提出一种改进的线性Radon变换方法来进行井孔资料的反射波上下行波场分离.该方法基于频率域线性Radon变换,通过引入一个新的变量λ来消除变换算子对频率的依赖性,避免了求取每一频率分量对应的不同变换算子,显著降低了计算成本;文中在求解该方法对应的最小二乘问题时,引入了发展较为成熟的高分辨率Radon变换技术来进一步提高波场分离的精度.采用本文方法进行井孔地震资料的上下行波场分离可以在保证分离精度的前提下有效地提高计算效率.根据上下行波在λ-f域内分布的特殊性,设计简单的滤波算子就可实现上下行波场的分离.最后通过合成数据试算以及实际资料处理（VSP数据和井间地震数据）验证了该方法的可行性和有效性.     

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.     

Spectral characteristics of natural and artificial seismic events in the Lop Nor test site, China

Seismic discriminants based on the spectral seismogram and spectral magnitude techniques have been tested to discriminate between three events; a nuclear explosion which took place in Lop Nor, China with m _b 6.1 and two earthquakes from the closest area with m _b 5.5 and 5.3, respectively. The spectral seismogram of the three events shows that the frequency content of the nuclear explosion differs from that of the earthquakes where the P wave is richier in high frequency content in the nuclear explosion than the corresponding earthquakes. It is also observed that the energy decays more rapidly for the nuclear explosion than for the earthquakes. Furthermore, the spectral magnitudes reveal significant differences in the spectra between the nuclear explosion and the two earthquakes. The observed differences appear to be quite enough to provide a reliable discriminant. The estimated stress drop from the magnitude spectra indicates a higher stress drop of the nuclear explosion relative to the earthquakes of the same tectonic region.     

Comment on ‘Aspects of 1D seismic modelling using the Goupillaud principle’ by Evert Slob and Anton Ziolkowski1

The paper by Slob and Ziolkowski (1993) is apparently a comment on my paper (Szaraniec 1984) on odd-depth structure. In fact the basic understanding of a seismogram is in question. The fundamental equation for an odd-depth model and its subsequent deconvolution is correct with no additional geological constraints. This is the essence of my reply which is contained in the following points.

• 1 The discussion by Slob and Ziolkowski suffers from incoherence. On page 142 the Goupillaud (1961) paper is quoted: “… we must use a sampling rate at least double that… minimum interval…”. In the following analysis of such a postulated model Slob and Ziolkowski say that “… two constants are used in the model: Δt as sampling rate and 2Δt as two-way traveltime”. By reversing the Goupillaud postulation all the subsequent criticism becomes unreliable for the real Goupillaud postulation as well as the odd-depth model.
• 2 Slob and Ziolkowski take into consideration what they call the total impulse response. This is over and above the demands of the fundamental property of an odd-depth model. Following a similar approach I take truncated data in the form of a source function, S(z), convolved with a synthetic seismogram (earth impulse response), R?(z), the free surface being included. The problem of data modelling is a crucial one and will be discussed in more detail below. By my reasoning, however, the function may be considered as a mathematical construction introduced purely to work out the fundamental property. In this connection there is no question of this construction having a physical meaning. It is implicit that in terms of system theory, K(z) stands for what is known as input impedance.
• 3 Our understandings of data are divergent but Slob and Ziolkowski state erroneously that: “Szaraniec (1984) gives (21) as the total impulse response…”. This point was not made. This inappropriate statement is repeated and echoed throughout the paper making the discussion by Slob and Ziolkowski, as well as the corrections proposed in their Appendix A, ineffective. Thus, my equation (2) is quoted in the form which is in terms of the reflection response G^sc and holds true at least in mathematical terms. No wonder that “this identity is not valid for the total impulse response” (sic), which is denoted as G(z). None the less a substitution of G for G^sc is made in Appendix A, equation (A3). The equation numbers in my paper and in Appendix A are irrelevant, but (A3) is substituted for (32) (both numbers of equations from the authors’ paper). Afterwards, the mathematical incorrectness of the resulting equation is proved (which was already evident) and the final result (A16) is quite obviously different from my equation (2). However, the substitution in question is not my invention.
• 4 With regard to the problem of data modelling, I consider a bi-directional ID seismic source located just below the earth's surface. The downgoing unit impulse response is accompanied by a reflected upgoing unit impulse and the earth response is now doubled. The total impulse response for this model is thus given by where (—r₀) =— 1 stands for the surface reflection coefficient in an upward direction. Thus that is to say, the total response to a unit excitation is identical with the input impedance as it must be in system theory. The one-directional 1D seismic source model is in question. There must be a reaction to every action. When only the downgoing unit impulse of energy is considered, what about the compensation?
• 5 In more realistic modelling, an early part of a total seismogram is unknown (absent) and the seismogram is seen in segments or through the windows. That is why in the usual approach, especially in dynamic deconvolution problems, synthetic data in the presence of the free surface are considered as an equivalent of the global reflection coefficient. It is implicit that model arises from a truncated total seismogram represented as a source function convolved with a truncated global reflection coefficient.

Validation or invalidation of the truncation procedure for a numerically specified model may be attempted in the frame of the odd-depth assumption. My equations (22) and (23) have been designed for investigating the absence or presence of truncated energy. The odd-depth formalism allows the possibility of reconstructing an earlier part of a seismogram (Szaraniec 1984), that is to say, a numerical recovery of unknown moments which are unlikely designed by Slob and Ziolkowski for the data.     

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.     

DYNAMIC PREDICTIVE DECONVOLUTION*

Dynamic predictive deconvolution makes use of an entire seismic trace including all primary and multiple reflections to yield an approximation to the subsurface structure. We consider plane-wave motion at normal incidence in an horizontally layered system sandwiched between the air and the basement rock. Energy degradation effects are neglected so that the layered system represents a lossless system in which energy is lost only by net transmission downward into the basement or net reflection upward into the air; there is no internal loss of energy by absorption within the layers. The layered system is frequency selective in that the energy from a surface input is divided between that energy which is accepted over time by net transmission downward into the basement and the remaining energy that is rejected over time by net reflection upward into the air. Thus the energy from a downgoing unit spike at the surface as input is divided between the wave transmitted by the layered system into the basement and the wave reflected by the layered system into the air. This reflected wave is the observed seismic trace resulting from the unit spike input. From surface measurements we can compute both the input energy spectrum, which by assumption is unity, and the reflection energy spectrum, which is the energy spectrum of the trace. But, by the conservation of energy, the input energy spectrum is equal to the sum of the reflection energy spectrum and the transmission energy spectrum. Thus we can compute the transmission energy spectrum as the difference of the input energy spectrum and the reflection energy spectrum. Furthermore, we know that the layered system acts as a pure feedback system in producing the transmitted wave, from which it follows that the transmitted wave is minimum-delay. Hence from the computed energy spectrum of the transmitted wave we can compute the prediction-error operator that contracts the transmitted wave to a spike. We also know that the layered system acts as a system with both a feedback component and a feed-forward component in producing the reflected wave, that is, the observed seismic trace. Moreover, this feedback component is identical to the pure feedback system that produces the transmitted wave. Thus, we can deconvolve the observed seismic trace by the prediction-error operator computed above; the result of the deconvolution is the wave-form due to the feedforward component alone. Now the feedforward component represents the wanted dynamic structure of the layered system whereas the feedback component represents the unwanted reverberatory effects of the layered system. Because this deconvolution process yields the wanted dynamic structure and destroys the unwanted reverberatory effects, we call the process dynamic predictive deconvolution. The resulting feedforward waveform in itself represents an approximation to the subsurface structure; a further decomposition yields the reflection coefficients of the interfaces separating the layers. In this work we do not make the assumption as is commonly done that the surface as a perfect reflector; that is, we do not assume that the surface reflection coefficient has magnitude unity.     

TUTORIAL THE PROGRESSIVE ATTENUATION OF HIGH-FREQUENCY ENERGY IN SEISMIC REFLECTION DATA*

Seismic reflection data always exhibit a progressive loss of high-frequency energy with time. This effect is partly attributable to irreversible processes such as the conversion of elastic energy into heat (commonly known as absorption), and partly to reversible processes associated with interference between reflected waves arriving at different times. This paper looks only at reversible linear elastic effects at normal incidence and asks the following question: if there were no such absorption, would there still be a progressive loss of high-frequency energy? Using normal incidence and a layered elastic earth model we prove the following results. 1. The normal incidence response of a sequence of plane parallel elastic layers is non-white. 2. The pressure wave reflected by a layer that is thin compared with a wavelength is differentiated with respect to the incident wave. 3. The transmission response of a thin layer is consequently low-pass and the transmission response of a sequence containing many thin layers is very low-pass. 4. The well-known effect of the transport of acoustic energy by peg-leg multiples within thin layers is identical with this low-pass transmission response. 5. It follows that the high frequency energy is reflected back early in the seismogram. 6. By comparison, very low-frequencies are transmitted through the layered sequence easily and are reflected with difficulty. There is probably a lack of low-frequency energy in the reflection seismogram, by comparison with the spectrum of the incident plane wave. It follows that any meaningful evaluation of frequency-dependent absorption in seismic data cannot take place unless the frequency-dependent linear elastic effects are taken into account first.     

MEASUREMENTS OF EARTH ATTENUATION FROM DOWNHOLE AND SURFACE SEISMIC RECORDINGS*

Frequency-dependent attenuation of compressional waves within the earth has been estimated in the vicinity of wells from

• 1 spectral power ratios of the coherent events in separate time gates on the seismic section
• 2 matching a broadband synthetic trace with seismic data at the well, and
• 3 determining the operator that transforms one down(up) going pulse recorded in the well into another recorded at a deeper (shallower) level.

The accuracy of estimation of all three methods was insufficient to estimate attenuation over small depth intervals, and it was not possible to distinguish between the contribution due to internal multiples and that of genuine absorption with much confidence. Spectral ratios from (1) showed a smoother variation with frequency—and one more consistent with other estimates—when they were compensated for the spectra of the reflectivities over the time gates employed, but they did not provide more than a broad indication of attenuation over a substantial depth interval. Approach (2) was hampered by the restricted durations over which synthetic trace and seismic data can be reliably matched; approach (3) gave the best results. Here matching is a much more powerful tool than the spectral-ratio techniques that are commonly applied since it can yield the form of the attenuation operator, i.e., both its amplitude and phase response, together with properly defined measures of its accuracy, while at the same time it minimizes the influence of noise and local interference effects at each recording level. For seismic target depths where internal multiple activity was low the logarithms of the amplitude responses of the estimated attenuation operators decreased approximately linearly with frequency and the phase responses showed no significant dispersion. Application of approach (3) to downgoing and upgoing waves estimated from a vertical seismic profile revealed the importance of changes in frequency-dependent geophone coupling and their effect on values of Q determined from downgoing pulses only.     

Noise transfer in variable‐depth streamer deghosting

Single‐component towed‐streamer marine data acquisition records the pressure variations of the upgoing compressional waves followed by the polarity‐reversed pressure variations of downgoing waves, creating sea‐surface ghost events in the data. The sea‐surface ghost for constant‐depth towed‐streamer marine data acquisition is usually characterised by a ghost operator acting on the upgoing waves, which can be formulated as a filtering process in the frequency–wavenumber domain. The deghosting operation, usually via the application of the inverse Wiener filter related to the ghost operator, acts on the signal as well as the noise. The noise power transfer into the deghosted data is proportional to the power spectrum of the inverse Wiener filter and is amplifying the noise strongly at the notch wavenumbers and frequencies of the ghost operator. For variable‐depth streamer acquisition, the sea‐surface ghost cannot be described any longer as a wavenumber–frequency operator but as a linear relationship between the wavenumber–frequency representation of the upgoing waves at the sea surface and the data in the space–frequency domain. In this article, we investigate how the application of the inverse process acts on noise. It turns out that the noise magnification is less severe with variable‐depth streamer data, as opposed to constant depth, and is inversely proportional to the local slant of the streamer. We support this statement via application of the deghosting process to real and numerical random noise. We also propose a more general concept of a wavenumber–frequency ghost power transfer function, applicable for variable‐depth streamer acquisition, and demonstrate that the inverse of the proposed variable‐depth ghost power transfer function can be used to approximately quantify the action of the variable‐depth streamer deghosting process on noise.     

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.     

基于格林函数理论的波场预测和鬼波压制方法

作为一种特殊的噪声,鬼波对一次波的波形及频带宽度产生极大的影响,鬼波压制是提高海上地震资料分辨率及保真度的重要因素.以格林公式为基础,详细论述了基于格林函数理论的鬼波压制方法,在不需要地下介质信息的条件下,进行地震数据驱动鬼波压制,并根据"Double Dirichlet"(双狄利克雷)边界条件,预测压力波场和垂直速度波场.建立了基于格林函数理论鬼波压制的处理流程,数值模拟和实际资料处理结果表明,基于格林函数理论鬼波压制方法在很好地去除鬼波的同时极大地拓宽了地震资料的频带,尤其提升了低频端能量,有利于后续资料的处理解释.     

QUASI-ANALYTIC CONVOLUTION SOLUTION OF THE ELECTROMAGNETIC FIELD*

The objective of this study is to generate the separation-distance-domain (r-domain) transformation of the theoretically calculated wave number domain (m-domain) electromagnetic induction field component B_z(m, ω) of a stratified medium and to search for interpretive information which has been absent in the previously achieved numerical solutions of the problem. The r-domain kernel R?(r, ω) function defining the induction field appears to adequately reflect the layering and electrical properties of the medium if it is expressed as a function of the frequency if the source-receiver separation r is small with respect to the thickness of the first layer. However, exact values of the conductivity cannot be distinguished from those of the neighboring values unless a resistive basement layer is present. This feature is the result of the truncation in series representation of the kernel function R?(m, ω). However, this truncation is regarded as significant in the case of a conductive first layer. In m-domain static-zone studies, a conductive first layer slightly influences its r-domain correspondent. Although the computational cost of obtaining the kernel B(r, ω) by evaluation of the convolution in a cylindrical coordinate system is high, this semi-analytic solution is still superior to those based on the asymptotic assumptions.