IDENTIFICATION OF SEISMIC REFLECTIONS USING SINGULAR VALUE DECOMPOSITION*

Singular value decomposition (SVD) is applied to the identification of seismic reflections by using two different models: the impulse response model, where a seismic trace is assumed to consist of a known signal pulse convolved with a reflection coefficient series plus noise, and the delayed pulse model, where the seismic signal is assumed to consist of a small number of delayed pulses of known shape and with unknown amplitudes and arrival times. SVD clearly shows how least-squares estimation of the reflection coefficients may become unstable, since a division by the singular values is required. Two methods for stabilizing this procedure are investigated. The inverse of the singular values may be replaced by zeros when they are less than a given threshold. This is called the SVD cut-off method. Alternatively, we may use ridge regression which in filter design corresponds to assuming white noise. Statistical methods are used to compute an optimal SVD cut-off level and also to compute an optimal weighting parameter in ridge regression. Numerical studies indicate that the use of SVD cut-off or ridge regression stabilizes the least-squares procedure, but that the results are inferior to maximum-likelihood estimation where the noise is assumed to be filtered white noise. For the delayed pulse model, we use a linearization procedure to iteratively update the estimates of both the reflection amplitudes and the arrival times. In each step, the optimal SVD cut-off method is used. Confidence regions for the estimated reflection amplitudes and arrival times are also computed. Synthetic data examples demonstrate the effectiveness of this method. In a real data example, the maximum-likelihood method assuming an impulse response model is first used to obtain initial estimates of the number of reflections and their amplitudes and traveltimes. Then the iterative procedure is used to obtain improved estimates of the reflection amplitudes and traveltimes.     

Random noise attenuation using a structure-oriented weighted singular value decomposition

Singular value decomposition (SVD) is a useful method for random noise suppression in seismic data processing. A structure-oriented SVD (SOSVD) approach which incorporates structure prediction to the SVD filter is effcient in attenuating noise except distorting seismic events at faults and crossing points. A modified SOSVD approach using a weighted stack, called structure-oriented weighted SVD (SOWSVD), is proposed. In this approach, the SVD filter is used to attenuate noise for prediction traces of a primitive trace which are produced via the plane-wave prediction. A weighting function related to local similarity and distance between each prediction trace and the primitive trace is applied to the denoised prediction traces stacking. Both synthetic and field data examples suggest the SOWSVD performs better than the SOSVD in both suppressing random noise and preserving the information of the discontinuities for seismic data with crossing events and faults.

     

A SEMBLANCE-GUIDED MEDIAN FILTER1

Reiter , E.C., Toksoz , M.N. and Purdy , G.M. 1992. A semblance-guided median filter. Geophysical Prospecting 41 , 15–41. A slowness selective median filter based on information from a local set of traces is described and implemented. The filter is constructed in two steps, the first being an estimation of a preferred slowness and the second, the selection of a median or trimmed mean value to replace the original data point. A symmetric window of traces defining the filter aperture is selected about each trace to be filtered and the filter applied repeatedly to each time point. The preferred slowness is determined by scanning a range of linear moveouts within the user-specified slowness passband. Semblance is computed for each trial slowness and the preferred slowness selected from the peak semblance value. Data points collected along this preferred slowness are then sorted from lowest to highest and in the case of a pure median filter, the middle point(s) selected to replace the original data point. The output of the filter is therefore quite insensitive to large amplitude noise bursts, retaining the well-known beneficial properties of a traditional 1D median filter. Energy which is either incoherent over the filter aperture or lies outside the slowness passband, may be additionally suppressed by weighting the filter output by the measured peak semblance. This approach may be used as a velocity filter to estimate coherent signal within a specified slowness passband and reject coherent energy outside this range. For applications of this type, other velocity estimators may be used in place of our semblance measure to provide improved velocity estimation and better filter performance. The filter aperture may also be extended to provide increased velocity estimation, but will result in additional lateral smearing of signal. We show that, in addition to a velocity filter, our approach may be used to improve signal-to-noise ratios in noisy data. The median filter tends to suppress the amplitude of random background noise and semblance weighting may be used to reduce the amplitude of background noise further while enhancing coherent signal. We apply our method to vertical seismic profile data to separate upgoing and downgoing wavefields, and also to large-offset ocean bottom hydrophone data to enhance weak refracted and post-critically reflected energy.     

AUTOMATIC EDITING OF NOISY SEISMIC DATA1

Seismic data often contain traces that are dominated by noise; these traces should be removed (edited) before multichannel filtering or stacking. Noise bursts and spikes should be edited before single channel filtering. Spikes can be edited using a running median filter with a threshold; noise bursts can be edited by comparing the amplitudes of each trace to those of traces that are nearby in offset-common midpoint space. Relative amplitude decay rates of traces are diagnostic of their signal-to-noise (S/N) ratios and can be used to define trace editing criteria. The relative amplitude decay rate is calculated by comparing the time-gated trace amplitudes to a control function that is the median trace amplitude as a function of time, offset, and common midpoint. The editing threshold is set using a data-adaptive procedure that analyses a histogram of the amplitude decay rates. A performance evaluation shows that the algorithm makes slightly fewer incorrect trace editing decisions than human editors. The procedure for threshold setting achieves a good balance between preserving the fold of the data and removing the noisiest traces. Tests using a synthetic seismic line show that the relative amplitude decay rates are diagnostic of the traces’S/N ratios. However, the S/N ratios cannot be accurately usefully estimated at the start of processing, where noisy-trace editing is most needed; this is the fundamental limit to the accuracy of noisy trace editing. When trace equalization is omitted from the processing flow (as in amplitude-versus-offset analysis), precise noisy-trace editing is critical. The S/N ratio of the stack is more sensitive to type 2 errors (failing to reject noisy traces) than it is to type 1 errors (rejecting good traces). However, as the fold of the data decreases, the S/N ratio of the stack becomes increasingly sensitive to type 1 errors.     

Median filter behaviour with seismic data1

Median filters may be used with seismic data to attenuate coherent wavefields. An example is the attenuation of the downgoing wavefield in VSP data processing. The filter is applied across the traces in the ‘direction’ of the wavefield. The final result is given by subtracting the filtered version of the record from the original record. This method of median filtering may be called ‘median filtering operated in subtraction’. The method may be extended by automatically estimating the slowness of coherent wavefields on a record. The filter is then applied in a time- and-space varying manner across the record on the basis of the slowness values at each point on the record. Median filters are non-linear and hence their behaviour is more difficult to determine than linear filters. However, there are a number of methods that may be used to analyse median filter behaviour: (1) pseudo-transfer functions to specific time series; (2) the response of median filters to simple seismic models; and (3) the response of median filters to steps that simulate terminating wavefields, such as faults on stacked data. These simple methods provide an intuitive insight into the behaviour of these filters, as well as providing a semiquantitative measurement of performance. The performance degradation of median filters in the presence of trace-to-trace variations in amplitude is shown to be similar to that of linear filters. The performance of median filters (in terms of signal distortion) applied obliquely across a record may be improved by low-pass filtering (in the t-dimension). The response of median filters to steps is shown to be affected by background noise levels. The distortion of steps introduced by median filters approaches the distortion of steps introduced by the corresponding linear filter for high levels of noise.     

半径-斜率域加权反假频地震数据重建

为减小地震数据缺失给地震后续处理工作带来的影响,需要对地震数据进行插值重建.针对反假频插值重建这个难点问题,进行了相关研究,并由此提出了一种改进的R-P（半径-斜率）域加权反假频地震数据插值重建方法.该方法将F-K（频率-波数）谱变换到R-P域,在R-P域设计一个权函数并将其作用于每次的迭代插值过程.通过模型数据和实际数据的测试,证明了该方法具有较好的反假频插值重建能力.     

L'EGALISATION SPECTRALE: UN MOYEN D'AMELIORER LA QUALITE DES DONNEES SISMIQUES*

In land seismic surveys spectrum equalization can increase the quality of seismic data in a selected frequency band. The power of lower frequencies in the spectrum of input traces is generally greater than that of higher frequencies, particularly in land seismic surveys because of ground roll. In order to improve the quality of seismic data it is necessary to raise the energy of higher frequencies to the same level as that of lower frequencies, without alteration of the phases. The first step of the method is to compute the amplitude spectrum of each input trace to determine a weighting function which is then applied to the amplitude spectrum in order to balance it. The function is the inverse of the short wavelength variation of the amplitude spectrum. The short wavelength variation can be obtained by interpolation between average values of the modulus of the amplitude spectrum computed in narrow bands within a selected band of frequencies. Another way of obtaining the short wavelength variation is to apply a low-pass filter to the amplitude spectrum. The calculations are readily performed in the frequency domain by the Fourier transform. Spectrum equalization is automatically adjusted to each trace and does not modify the average amplitude in the time domain. However, as the frequency band and energy of the ground roll both vary according to the distance from the shot, spectrum equalization tends to make the spectrum of output traces independent of the offset distance. The use of spectrum equalization before any two-dimensional filtering improves ground roll elimination. Continuity and resolution of horizons are also increased by spectrum equalization before CDP stack. Several examples of applications of spectrum equalization to seismic land and marine surveys are shown.     

THREE-DIMENSIONAL OPTIMUM MULTICHANNEL VELOCITY FILTERS*

Design procedures and characteristics of three stacking filters are discussed which may find application in various three-dimensional velocity filtering problems. These filters are derived in the time-domain as optimum multichannel Wiener filters. Random stationary functions are assumed as stochastic models for the seismic traces. All power and crosspower spectra which are the basic elements of the multichannel normal equations are statistically averaged according to specific three-dimensional considerations. Various properties of the input traces may be incorporated in the design of the optimum filters. With fairly general assumptions about the input these filters are deterministic in the sense that they are applicable to a broad class of input traces with similar statistics in amplitudes and arrival times of signals.     

基于非平稳相似性系数的构造导向滤波及断层检测方法

不连续地质体（如断层）的自动检测一直以来都是叠后地震数据解释中的关键问题之一,尤其在三维情况中尤为重要.然而,大多数边缘检测和相干算法都对随机噪声很敏感,随机噪声衰减是叠后地震数据解释的另一个主要问题.针对构造保护去噪和断层检测问题,本文基于非平稳相似性系数完善一种构造导向滤波方法并且提出一种自动断层检测方法,形成了一套匹配的处理技术.该构造导向滤波既能够有效地衰减随机噪声又可以很好地保护地震资料中的断层等信息不被破坏,增强地震剖面中弯曲、倾斜同相轴的连续性.根据地震数据局部倾角走向,利用相邻道构建当前地震道的预测,通过预测道的叠加得到参考道,计算预测道与参考道之间的非平稳相似性系数可以设计出数据驱动的加权中值滤波.另一方面,预测道与原始道之间的非平稳相似性系数能够用于带有断层指示性的相干分析.这两种方法都基于构造预测和非平稳相似性系数,但是使用不同的调节参数和处理方案.理论模型和实际数据的处理结果证明了本文提出构造导向滤波和断层检测方法的有效性.     

模糊嵌套多级中值滤波方法及其在地震数据处理中的应用

多级中值滤波器在地震数据处理应用中,其滤波长度越大,消噪效果越好,但同时也会破坏有效信息. 如果为了保护有效信息而减少滤波长度,又会造成大量的噪声不能消除. 本文提出一种新的模糊嵌套多级中值滤波器,设计一个阀值作为判断参数,使滤波器能够在消除随机噪声时采用长滤波器滤波,而在保留有效信息时采用短滤波器滤波,从而既能很好地消除随机噪声,又能最大限度地保护有效信息,保留有效信息的细节结构. 经过模型分析和实际资料处理都取得了很好的效果.     

A COMPARISON OF STATISTICAL AND DETERMINISTIC WIENER DECONVOLUTION OF DEEP-TOW SEISMIC DATA*

Wiener ‘spiking’ deconvolution of seismic traces in the absence of a known source wavelet relies upon the use of digital filters, which are optimum in a least-squares error sense only if the wavelet to be deconvolved is minimum phase. In the marine environment in particular this condition is frequently violated, since bubble pulse oscillations result in source signatures which deviate significantly from minimum phase. The degree to which the deconvolution is impaired by such violation is generally difficult to assess, since without a measured source signature there is no optimally deconvolved trace with which the spiked trace may be compared. A recently developed near-bottom seismic profiler used in conjunction with a surface air gun source produces traces which contain the far-field source signature as the first arrival. Knowledge of this characteristic wavelet permits the design of two-sided Wiener spiking and shaping filters which can be used to accurately deconvolve the remainder of the trace. In this paper the performance of such optimum-lag filters is compared with that of the zero-lag (one-sided) operators which can be evaluated from the reflected arrival sequence alone by assuming a minimum phase source wavelet. Results indicate that the use of zero-lag operators on traces containing non-minimum phase wavelets introduces significant quantities of noise energy into the seismic record. Signal to noise ratios may however be preserved or even increased during deconvolution by the use of optimum-lag spiking or shaping filters. A debubbling technique involving matched filtering of the trace with the source wavelet followed by optimum-lag Wiener deconvolution did not give a higher quality result than can be obtained simply by the application of a suitably chosen Wiener shaping filter. However, cross correlation of an optimum-lag spike filtered trace with the known ‘actual output’ of the filter when presented with the source signature is found to enhance signal-to-noise ratio whilst maintaining improved resolution.     

Spectral factorization technique for estimation of an ARMA operator for multichannel deconvolution of seismic data

A new spectral factorization method is presented for the estimation of a causal as well as a causally invertible ARMA operator from the correlation sequence of seismic traces. The method has been implemented for multichannel deconvolution of seismic traces with the aim of exploiting the trace-to-trace correlation that exists within seismograms. A layered earth model with a small reflectivity sequence has been considered, and the seismic traces have been considered as the output of a linear system driven by white noise reflection coefficient sequences. The present method is the concatenation of three algorithms, namely Kung's method for state variable ( F , G , H ) realization using a singular value decomposition (SVD) algorithm, Faurre's technique for computation of the strong spectral factor and Leverrier's algorithm for ARMA representation of the spectral factor. The inverted ARMA operator is used as a recursive filter for deconvolution of seismic traces. In the example shown, two traces with a covariance sequence of 160 ms length have been considered for multichannel deconvolution of stacked seismic traces. The results presented, when compared with those obtained from a conventional deconvolution algorithm, have shown encouraging prospects.     

A procedure to reduce side lobes of reflection wavelets: A contribution to low frequency information

Side lobes of the wavelets arise from the lack of low frequency content in a reflection wavelet. They tend to increase the time span of an individual reflection event and interfere with the other primary reflections or side lobes. Furthermore, their trace-by-trace consistency may produce pseudo-reflections and may cause misinterpretations of the side lobes as weak reflections.A procedure in order to improve the low frequency content of the seismic traces by suppressing the side lobe amplitudes based on the complex trace envelope is proposed. Using the average energies of the seismic trace and its envelope, the polarity table of the trace is obtained and used to correct the phase of the envelope. The resultant trace is termed “side lobe reduced (SLR) trace”. The method can be applied to the stack or migrated seismic data by a trace-by-trace basis. The only required parameter of the method is the moving average operator length which is used to calculate average energies of the input traces. In general, shorter operator lengths yield better results when the dominant frequency of the input increases.Results from synthetics and real seismic data sets show that the procedure improves the low frequency components of the input trace and side lobes in the output SLR trace are significantly suppressed. The method may be considered as a seismic amplitude attribute, which aids the interpreter to obtain the true seismic signature of the geological formations by removing the side lobes of the wavelet and restoring the low frequency components if the lower frequencies of deeper reflections are of primary concern.     

Enhancement of Deep Seismic Reflections in Pre-stack Data by Adaptive Filtering

—Adaptive filters offer advantages over Wiener filters for time-varying processes. They are used for deconvolution of seismic data which exhibit non-stationary behavior, and seldom for noise reduction. Different algorithms for adaptive filtering exist. The least-mean-squares (LMS) algorithm, because of its simplicity, has been widely applied to data from different fields that fall outside geophysics. The application of the LMS algorithm to improve the signal-to-noise ratio in deep reflection seismic pre-stack data is studied in this paper. Synthetic data models and field data from the DEKORP project are used to this end.¶Three adaptive filter techniques, one-trace technique, two-trace technique and time-slice technique, are examined closely to establish the merits and demerits of each technique. The one-trace technique does not improve the signal-to-noise ratio in deep reflection seismic data where signal and noise cover the same frequency range. With the two-trace technique, the strongest noise reduction is achieved for small noise on the data. The filter efficiency decreases rapidly with increasing noise. Furthermore, the filter performance is poor upon application to common-midpoint (CMP) gathers with no normal-moveout (NMO) corrections. Application of the two-trace method to seismic traces before dynamic correction results in gaps in the signal along the reflection hyperbolas. The time-slice technique, introduced in this paper, offers the best answer. In this case, the one-trace technique is applied to the NMO-corrected gathers across all traces in each gather at each time to separate the low-wavenumber component of the signal in offset direction from the high-wavenumber noise component. The stacking velocities used for the dynamic correction do not need to be known very accurately because in deep reflection seismics, residual moveouts are small and have only a minor influence on the results of the adaptive time-slice technique. Noise reduction is more significant with the time-slice technique than with the two-trace technique. The superiority of the adaptive time-slice technique is demonstrated with the DEKORP data.     

VERTICAL SEISMIC PROFILING—SEPARATION OF P- AND S-WAVES*

In vertical seismic profile's (VSP's) shot with a large source offset, rays from shot to receiver can have large angles of incidence. Shear waves generated by the source and by conversions at interfaces are likely to be recorded by both the vertical and the horizontal geophones. Varying angles of incidence may give strong variations in the recorded amplitudes. Separation of P- and SV-waves and recovery of their full amplitudes are important for proper processing and interpretation of the data. A P-S separation filter for three-component offset VSP data is presented which performs this operation. The separation filter is applied in the k-f domain and needs an estimate of the P- and S-velocities along the borehole as input. Implementation and stability aspects of the filter are considered. The filter was tested on an 1800 m offset VSP and appeared to be robust. Large velocity variations along the borehole could be handled and results were superior to those obtained by velocity filtering.     

Multiple suppression by 2D filtering in the parabolic τ–p domain: a wave-equation-based method1

The filter for wave-equation-based water-layer multiple suppression, developed by the authors in the x-t, the linear τ-p, and the f-k domains, is extended to the parabolic τ-2 domain. The multiple reject areas are determined automatically by comparing the energy on traces of the multiple model (which are generated by a wave-extrapolation method from the original data) and the original input data (multiples + primaries) in τ-p space. The advantage of applying the data-adaptive 2D demultiple filter in the parabolic τ-p domain is that the waves are well separated in this domain. The numerical examples demonstrate the effectiveness of such a dereverberation procedure. Filtering of multiples in the parabolic τ-p domain works on both the far-offset and the near-offset traces, while the filtering of multiples in the f-k domain is effective only for the far-offset traces. Tests on a synthetic common-shot-point (CSP) gather show that the demultiple filter is relatively immune to slight errors in the water velocity and water depth which cause arrival time errors of the multiples in the multiple model traces of less than the time dimension (about one quarter of the wavelet length) of the energy summation window of the filter. The multiples in the predicted multiple model traces do not have to be exact replicas of the multiples in the input data, in both a wavelet-shape and traveltime sense. The demultiple filter also works reasonably well for input data contaminated by up to 25% of random noise. A shallow water CSP seismic gather, acquired on the North West Shelf of Australia, demonstrates the effectiveness of the technique on real data.     

AVOA techniques for fracture characterization

Different aspects of computational techniques for AVOA analysis (Amplitude Versus Offset and Azimuth) for fracture characterization are considered, in particular: using amplitudes instead of reflection coefficients, smoothing seismic data, and numerical methods for estimation of fracture directions. A new computational method and a new filter for smoothing are suggested. The different computational methods are compared in synthetic reflection surface data with noise, and without noise. Properties of the numerical methods in dependence on different sets of azimuths and offsets are obtained. It is shown a superiority of the new method.     

MAXIMUM LIKELIHOOD ESTIMATION OF SEISMIC REFLECTION COEFFICIENTS*

A seismic trace is modeled as a moving average (MA) process both in signal and noise: a signal wavelet convolved with a reflection coefficient series plus colored random noise. Seismic reflection coefficients can be estimated from seismic traces using suitable estimation algorithms if the input wavelet is known and vice versa. The maximum likelihood (ML) algorithm is used to estimate the system order and the reflection coefficients. The system order is related to the arrival time of the latest signal in a complex seismic reflection event. The least-squares (LS) method does not provide such information. The ML algorithm makes assumptions only about the Gaussian nature of the noise. It is better suited for seismic applications since the LS method inherits the white noise assumption. The Gauss-Newton (G-N) and Newton-Raphson (N-R) optimization algorithms are used to obtain the ML and the LS estimates. Reflection coefficient estimations are affected by the choice of sampling rate of seismic data. Theoretically, the optimum choice in system identification is the Nyquist rate. Experience with synthetic data confirms the theory. In practice, good estimates of reflection coefficients are possible only up to certain pulse separations (or, equivalently, orders). This is mostly due to numerical problems with the optimization algorithms used and partly due to the limited bandwidth of seismic signals. Good estimates from data simulated using three airgun array pulses recorded with 6–128 Hz filter setting are possible up to about 40.0 ms pulse separations. Successful estimations from pinchout and thin layer simulations and well controlled offshore “bright-spots” are given.     

INVERSION METHOD IN THE SLANT STACK DOMAIN USING AMPLITUDES OF REFLECTION ARRIVALS*

Almost all ray-tracing methods ignore the analysis of the amplitudes of seismic arrivals and therefore utilize only half of the available information. We propose a method which is a combination of ray-tracing imaging and transformation of the amplitudes of wide-aperture data. Seismic data in the conventional X-T domain are first transformed to the domain of intercept time τ and ray parameter p to recover the plane wave response. The next step is the derivation of a series of plane wave reflection coefficients, which are mapped as a function of τ and p. The reflection coefficients R(τ, p) for two arbitrarily chosen traces can then be used in our inversion method to derive a slowness-depth and a density-depth profile. It is shown that the inclusion of amplitudes of seismic arrivals (in this method, we consider the acoustic case) makes the inverse method highly stable and accurate. In a horizontally stratified medium one can recover separate profiles of velocity and density. Since this method utilizes large-offset data, it can be used for separate recovery of velocity and density to a greater depth.     

三维共偏移距叠前地震数据反假频射线束形成偏移成像

对稀疏/非规则采样或者低信噪比数据,射线束提取困难并伴随有假频产生,对叠加剖面和道集造成严重干扰.为了提升射线束偏移在稀疏和低信噪比地震数据采集中的成像效果,本文提出基于三角滤波的局部倾斜叠加波束形成偏移假频压制方法.射线束偏移首先将地震数据划分为超道集,经过部分NMO后转化为以射线束中心定义的共偏移距数据,倾斜叠加和反假频操作均在局部共中心点坐标上实现.时间域倾斜叠加是对地震数据的时移累加操作,三角低通滤波同样可以在时间域完成,在对地震数据进行因果和反因果积分后,亦为地震数据的时移累加.因此,三角低通滤波与倾斜叠加可在时间域结合同时完成,避免了频域滤波的正反傅里叶变换.本文在反假频公式中加入权重系数,用以对反假频的程度进行控制,达到分辨率和噪声压制的最佳折衷.以某海上三维实际数据为例,文中展示了反假频射线束形成对偏移叠加剖面和共成像点偏移距道集中的噪声进行了有效压制.