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.     

ODD-DEPTH STRUCTURES FOR DETERMINISTIC DECONVOLUTION AND SEISMOGRAM TESTING*

In odd-depth structure the two-way traveltime to each boundary is constrained to be an odd integer. The odd-depth property of a model is exposed to possible refutation under a seismogram test. Test function is a simple transformation of a synthetic seismogram. For an odd-depth model the test function has identically the value 1. The testability of a synthetic seismogram over an odd-depth structure provides a method of deterministic deconvolution. There is no need of specialized assumptions, like the minimum-phase property, about the source wavelet. The deconvolution may be performed in the absence of the early segment of a seismogram.     

ON DIRECT RECOVERY OF THE IMPULSE RESPONSE*

The deconvolution equation is solved in the z-transform domain directly for an impluse response. The principal assumption is the odd-depth model: two-way traveltimes to the boundaries are constrained to be odd integers only. It is further assumed that the length of the wavelet sequence is known to be less than half the length of the data sequence. An inverse of the impulse response is constrained by the zero samples of the source function. The resulting underdetermined set of equations is supplemented with the equations provided by the odd-depth model. The impulse response is found from the inverse by polynomial division.     

A FORWARD GENERATED SYNTHETIC SEISMOGRAM FOR EQUAL-DELAY LAYERED MEDIA MODELS*

A forward solution for the reflection response of a parallel stratified lossless medium characterized by discrete reflection coefficients and unequal layer delays, for a normally incident pressure source signal, is presented. The notation, which details the reflection history of each wavelet in a response record, facilitates systematic enumeration of all terms in the reflection impulse response model, the determination of compact closed form expressions for amplitudes and delays of multiply reflected wavelets, and the aggregation of dynamic analog groups. An equal delay time constraint on layer thicknesses leads then to the reflection sequence or synthetic seismogram structure as an infinite sum of wavelets by their order of reflection.     

STABILIZATION OF NORMAL-INCIDENCE SEISMOGRAM INVERSION REMOVING THE NOISE-INDUCED BIAS*

A main problem in computing reflection coefficients from seismograms is the instability of the inversion procedure due to noise. This problem is attacked for two well-known inversion schemes for normal-incidence reflection seismograms. The crustal model consists of a stack of elastic, laterally homogeneous layers between two elastic half-spaces. The first method, which directly computes the reflection coefficients from the seismogram is called “Dynamic Deconvolution”. The second method, here called “Inversion Filtering”, is a two-stage procedure. The first stage is the construction of a causal filter by factorization of the spectral function via Levinson-recursion. Filtering the seismogram is the second stage. The filtered seismogram is a good approximation for the reflection coefficients sequence (unless the coefficients are too large). In the non-linear terms of dynamic deconvolution and Levinson-recursion the noise could play havoc with the computation. In order to stabilize the algorithms, the bias of these terms is estimated and removed. Additionally incorporated is a statistical test for the reflection coefficients in dynamic deconvolution and the partial correlation coefficients in Levinson-recursion, which are set to zero if they are not significantly different from noise. The result of stabilization is demonstrated on synthetic seismograms. For unit spike source pulse and white noise, dynamic deconvolution outperforms inversion filtering due to its exact nature and lesser computational burden. On the other hand, especially in the more realistic bandlimited case, inversion filtering has the great advantage that the second stage acts linearly on the seismogram, which allows the calculation of the effect of the inversion procedure on the wavelet shape and the noise spectrum.     

VELOCITY-STACK PROCESSING1

A conventional velocity-stack gather consists of constant-velocity CMP-stacked traces. It emphasizes the energy associated with the events that follow hyperbolic traveltime trajectories in the CMP gather. Amplitudes along a hyperbola on a CMP gather ideally map onto a point on a velocity-stack gather. Because a CMP gather only includes a cable-length portion of a hyperbolic traveltime trajectory, this mapping is not exact. The finite cable length, discrete sampling along the offset axis and the closeness of hyperbolic summation paths at near-offsets cause smearing of the stacked amplitudes along the velocity axis. Unless this smearing is removed, inverse mapping from velocity space (the plane of stacking velocity versus two-way zero-offset time) back to offset space (the plane of offset versus two-way traveltime) does not reproduce the amplitudes in the original CMP gather. The gather resulting from the inverse mapping can be considered as the model CMP gather that contains only the hyperbolic events from the actual CMP gather. A least-squares minimization of the energy contained in the difference between the actual CMP gather and the model CMP gather removes smearing of amplitudes on the velocity-stack gather and increases velocity resolution. A practical application of this procedure is in separation of multiples from primaries. A method is described to obtain proper velocity-stack gathers with reduced amplitude smearing. The method involves a t²-stretching in the offset space. This stretching maps reflection amplitudes along hyperbolic moveout curves to those along parabolic moveout curves. The CMP gather is Fourier transformed along the stretched axis. Each Fourier component is then used in the least-squares minimization to compute the corresponding Fourier component of the proper velocity-stack gather. Finally, inverse transforming and undoing the stretching yield the proper velocity-stack gather, which can then be inverse mapped back to the offset space. During this inverse mapping, multiples, primaries or all of the hyperbolic events can be modelled. An application of velocity-stack processing to multiple suppression is demonstrated with a field data example.     

DETECTION AND RESOLUTION OF THIN LAYERS: A MODEL SEISMIC STUDY1

The detection and resolution of a thin layer closely situated above a high-impedance basement are predominantly determined by both the frequency content of the incident seismic wavelet and the existence of the nearby high-impedance bedrock. The separation of the thin layer and the basement arrivals is investigated depending on the low-frequency content of the wavelet. The high-frequency content of the wavelet is kept constant. The initial wavelet spectrum with low frequencies has a rectangular shape. All wavelets used have zero-phase characteristics. Numerical and analogue seismic modelling techniques are used. The study is based on the geology of the Pachangchi Sandstone in West Taiwan. Firstly the resolution of a thin layer between two half-spaces is examined by applying the Ricker and De Voogd-Den Rooijen criteria. The lack of low-frequency components of the incident seismic wavelet reduces the shortest true two-way traveltime by about 20%. In addition, low-frequency components of the wavelet diminish the deviation between true and apparent two-way traveltime by about 65% for layer thicknesses in the transition from a thick to a thin layer. The second step deals with the influence of a high-impedance basement just below a thin layer on the detection and resolution of that thin layer. Reflected signal energies and apparent two-way traveltimes are considered. The reflected signal energy depends on the low-frequency content of the incident wavelet, the layer's thickness and the distance between the basement and the layer. This applies only to layers with thicknesses less than or equal to one-third of the mean wavelength in the layer, and a distance to basement in the range of one to one-half of the mean wavelength in the rock material between layer and basement. The minimum thin-layer thickness resolvable decreases with increasing distance to the basement; i.e. for a layer thickness of one-third of the mean wavelength in the layer the relative error of the two-way traveltime increases from 5% to 30%, if the distance is reduced from one to one-half of the mean wavelength in the material between the basement and the thin layer. Finally, a combination of vertical seismic profiling and downward-continuation techniques is presented as a preprocessing procedure to prepare realistic data for the detection and resolution investigation.     

TIME-DOMAIN COMPUTATION OF NON-NORMAL INCIDENCE WAVEFIELDS IN PLANE LAYERED MEDIA1

A new time-domain method is introduced for the calculation of theoretical seismograms which include frequency dependent effects like absorption. To incorporate these effects the reflection and transmission coefficients become convolutionary operators. The method is based on the communication theory approach and is applicable to non-normal incidence plane waves in flat layered elastic media. Wave propagation is simulated by tracking the wave amplitudes through a storage vector inside the computer memory representing a Goupillaud earth model discretized by equal vertical transit times. Arbitrary numbers of sources and receivers can be placed at arbitrary depth positions, while the computational effort is independent of that number. Therefore, the computation of a whole plane-wave vertical seismic profile is possible with no extra effort compared to the computation of the surface seismogram. The new method can be used as an aid to the interpretation of plane-wave decomposed reflection data where the whole synthetic vertical seismic profile readily gives the interpreter the correct depth position of reflection 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.     

2.5D full-wavefield viscoacoustic invercion1

Full-wavefield inversion for distributions of acoustic velocity, density and Q on a vertical slice through a25D model is implemented for common-source gathers in a cross-hole geometry. The wavefield extrapolation used is 3D, so all geometrical spreading, scattering, reflection, and transmission effects are correctly and automatically compensated for. In order to keep the number of unknowns tractable, application was limited to 2.5D models of known geometry; the latter assurnes a prior step, such as tomography, to fix the layer geometries. With the model geometry fixed, reliable solutions are obtained using synthetic data from only two independent source locations. Solutions from data with noisy and missing traces are comparable to those from noise-free data, but with higher residuals. When the source locations are spatially widely separated, conunon-source gathers may be summed and treated as a single wavefield to yield the same model estimates as when the individual source wavefields are treated separately, at substantially reduced cost. Inversions for full 3D parameter distributions can be handled with the same software, requiring only solution for more unknowns.     

速度分布统计特征与地震拟测井问题

根据声测井得出的速度分布统计特征,给出了多层构造反射系数序列的自相关函数形式。在反射系数强度不太大的条件下,导出存在于速度分布统计特征与多层构造反射传递函数之间的一种简单关系。据此,对目前常用的线性地震拟测井处理过程进行了具体分析,指出畸变恒不可免;在此基础上,提出了实际处理中力求减小畸变的可能途径。     

ON ERRORS OF FIT AND ACCURACY IN MATCHING SYNTHETIC SEISMOGRAMS AND SEISMIC TRACES*

A synthetic seismogram that closely resembles a seismic trace recorded at a well may not be at all reliable for, say, stratigraphic interpretation around the well. The most accurate synthetic seismogram is, in general, not the one that displays the smallest errors of fit to the trace but the one that best estimates the noise on the trace. If the match is confined to a short interval of interest or if the seismic reflection wavelet is allowed to be unduly long, there is considerable danger of forcing a spurious fit that treats the noise on the trace as part of the seismic reflection signal instead of making a genuine match with the signal itself. This paper outlines tests that allow an objective and quantitative evaluation of the accuracy of any match and illustrates their application with practical examples. The accuracy of estimation is summarized by the normalized mean square error (NMSE) in the estimated reflection signal, which is shown to be (/n)(P_N/P_S) where P_S/P_N is the signal-to-noise power ratio and n is the spectral smoothing factor. That is, the accuracy varies directly with the ratio of the power in the signal (taken to be the synthetic) to that in the noise on the seismic trace, and the smoothing acts to improve the accuracy of the predicted signal. The construction of confidence intervals for the NMSE is discussed. Guidelines for the choice of the spectral smoothing factor n are given. The variation of wavelet shape due to different realizations of the noise component is illustrated, and the use of confidence intervals on wavelet phase is recommended. Tests are described for examining the normality and stationarity of the errors of fit and their independence of the estimated reflection signal.     

LEVINSON INVERSION OF EARTHQUAKE GEOMETRY SH-TRANSMISSION SEISMOGRAMS IN THE PRESENCE OF NOISE*

The Kunetz-Claerbout equation for the acoustic transmission problem in a layered medium in its original form establishes the relation between the transmission and the reflec tion response for P-waves in an horizontally layered medium and with vertical incidence. It states that the reflection seismogram due to an impulsive source at the surface, is one side of the autocorrelation of the seismogram due to an impulsive source at depth and a surface receiver. By adapting Claerbout's formulation to the transmission of SH-waves, the Kunetz-Claerbout equation also holds for reflection and transmission coefficients dependent on the incident angle. Thus, earthquake geometry SH-transmission seismograms can be used to caculate corresponding pseudoreflection seismograms which, in turn, can be inverted for the impedance structure using the Levinson algorithm. If the average incidence angle is known, a geometrical correction on the resulting impedance model can improve the resolution of layer thicknesses. In contrast to the inversion of reflection seismograms, the Levinson algorithm is shown to yield stable results for the inversion of transmission seismograms even in the presence of additive noise. This noise stabilization is inherent to the Kunetz-Claerbout equation. Results of inverted SH-wave microearthquake seismograms from the Swabian Jura, SW Germany, seismic zone obtained at recording site Hausen im Tal have been compared with sonic-log data from nearby exploration drilling at Trochtelfingen. The agreement of the main structural elements is fair to a depth of several hundred metres.     

地震波在垂向不均匀介质中传播问题的某些探讨

本文计算了含有高速夹层介质中首波的理论地震图。通过分析得到,当高速夹层薄到一定程度时,就会产生干涉型首波,从而从一个侧面证明了射线理论的局限性。通过对地震波反射—折射系数能量守恒关系的分析,探讨了反射—折射系数大于1的可能性。最后,介绍了一种计算垂向不均匀介质中拉梅问题理论地震图的数值方法——有限差分法。     

Source/receiver array directivity effects on marine seismic attenuation measurements

Attenuation of seismic waves, quantified by the seismic quality factor Q, holds important information for seismic interpretation, due to its sensitivity to rock and fluid properties. A recently published study of Q, based on surface seismic reflection data, used a modified spectral ratio approach (QVO), but both source and receiver responses were treated as isotropic, based on simple raypath arguments. Here, this assumption has been tested by computing apparent attenuation generated by frequency-dependent directivity of typical marine source and receiver arrays and acquisition geometries. Synthetic wavelet spectra were computed for reflected rays, summed over the first Fresnel zone, from the base of a single interval, 50–3000 m thick and velocity 2000 m/s, overlying a 2200 m/s half-space, and for offsets of 71–2071 m. The source and receiver geometry were those of an actual survey. The modelled spectra are clearly affected by directivity, most strongly because of surface ghosts. In general, the strong high-frequency component, produced by the array design, leads to apparently negative attenuation in individual reflection events, though this is dependent on offset and target depth. For shallow targets (less than 400–500 ms two-way traveltime (TWT) depth), apparent Q-values as extreme as ?50 to ?100 were obtained. For deeper target depths, the directivity effect is far smaller. The implications of the model study were tested on real data. QVO was applied to 20 true-spectrum-processed CMPs, in a shallow (405–730 ms TWT) and a deeper (1000–1300 ms TWT) interval, firstly using a measured far-field source signature (effectively isotropic), and secondly using computed directivity effects instead. Mean interval Q^?1-values for the deeper interval, 0.029 ± 0.011 and 0.027 ± 0.018 for conventional and directional processing, respectively, suggested no directivity influence on attenuation estimation. For the shallow interval (despite poor spectral signal-to-noise ratios and hence scattered attenuation estimates), directional processing removed directivity-generated irregularities from the spectral ratios, resulting in an improvement from Q^?1_int = ?0.036 ± 0.130 to a realistic Q^?1_int = 0.012 ± 0.030: different at 94% confidence level. Equivalent Q-values are: for the deeper interval, 35 and 37 for conventional and directional processing, respectively, and ?28 and 86 for the shallow interval. These results support the conclusions of the model studies, i.e. that source/receiver directivity has a negligible effect except for shallow targets (e.g. TWT depth ≤ 500 ms) imaged with conventional acquisition geometry. In such cases directivity corrections to spectra are strongly recommended.     

PROCESSING OF PS-REFLECTION DATA APPLYING A COMMON CONVERSION-POINT STACKING TECHNIQUE1

Converted waves require special data processing as the wave paths are asymmetrical. The CMP concept is not applicable for converted PS waves, instead a sorting algorithm for a common conversion point (CCP) has to be applied. The coordinates of the conversion points in a single homogeneous layer can be calculated as a function of the offset, the reflector depth and the velocity ratio v_P/ v_s. For multilayered media, an approximation for the coordinates of the conversion points can be made. Numerical tests show that the traveltime of PS reflections can be approximated with sufficient accuracy for a certain offset range by a two-term series truncation. Therefore NMO corrections can be calculated by standard routines which use the hyperbolic approximation of the reflection traveltime curves. The CCP-stacking technique has been applied to field data which have been generated by three vertical vibrators. The in-line horizontal components have been recorded. The static corrections have been estimated from additional P- and SH-wave measurements for the source and geophone locations, respectively. The data quality has been improved by processes such as spectral balancing. A comparison with the stacked results of the corresponding P- and SH-wavefield surveys shows a good coherency of structural features in P-, SH- and PS-time sections.     

起伏地形下的高精度反射波走时层析成像方法

全球造山带及中国大陆中西部普遍具有强烈起伏的地形条件.复杂地形条件下的地壳结构成像问题像一面旗帜引领了当前矿产资源勘探和地球动力学研究的一个重要方向.深地震测深记录中反射波的有效探测深度可达全地壳乃至上地幔顶部,而初至波通常仅能探测上地壳浅部.为克服和弥补初至波探测深度的不足,本文基于前人对复杂地形条件下初至波成像的已有研究成果,采用数学变换手段将笛卡尔坐标系的不规则模型映射到曲线坐标系的规则模型,并将快速扫描方法与分区多步技术相结合,发展了反射波走时计算和射线追踪的方法.进而利用反射波走时反演,实现起伏地形下高精度的速度结构成像,从而为起伏地形下利用反射波数据高精度重建全地壳速度结构提供了一种全新方案.数值算例从正演计算精度、反演中初始模型依赖性、反演精度、纵横向分辨率以及抗噪性等方面验证了算法的正确性和可靠性.     

Transmission + reflection anisotropic wave-equation traveltime and waveform inversion

A transmission + reflection wave-equation traveltime and waveform inversion method is presented that inverts the seismic data for the anisotropic parameters in a vertical transverse isotropic medium. The simultaneous inversion of anisotropic parameters and ε is initially performed using transmission wave-equation traveltime inversion method. Transmission wave-equation traveltime only provides the low-intermediate wavenumbers for the shallow part of the anisotropic model; in contrast, reflection wave-equation traveltime estimates the anisotropic parameters in the deeper section of the model. By incorporating a layer-stripping method with reflection wave-equation traveltime, the ambiguity between the background-velocity model and the depths of reflectors can be greatly mitigated. In the final step, multi-scale full-waveform inversion is performed to recover the high-wavenumber component of the model.  We use a synthetic model to illustrate the local minima problem of full-waveform inversion and how transmission and reflection wave-equation traveltime can mitigate this problem. We demonstrate the efficacy of our new method using field data from the Gulf of Mexico.