Detection and extraction of orientation-and-scale-dependent information from two-dimensional GPR data with tuneable directional wavelet filters

The Ground Probing Radar (GPR) is a valuable tool for near surface geological, geotechnical, engineering, environmental, archaeological and other work. GPR images of the subsurface frequently contain geometric information (constant or variable-dip reflections) from various structures such as bedding, cracks, fractures, etc. Such features are frequently the target of the survey; however, they are usually not good reflectors and they are highly localized in time and in space. Their scale is therefore a factor significantly affecting their detectability. At the same time, the GPR method is very sensitive to broadband noise from buried small objects, electromagnetic anthropogenic activity and systemic factors, which frequently blurs the reflections from such targets.This paper introduces a method to de-noise GPR data and extract geometric information from scale-and-dip dependent structural features, based on one-dimensional B-Spline Wavelets, two-dimensional directional B-Spline Wavelet (BSW) Filters and two-dimensional Gabor Filters. A directional BSW Filter is built by sidewise arranging s identical one-dimensional wavelets of length L, tapering the s-parallel direction (span) with a suitable window function and rotating the resulting matrix to the desired orientation. The length L of the wavelet defines the temporal and spatial scale to be isolated and the span determines the length over which to smooth (spatial resolution). The Gabor Filter is generated by multiplying an elliptical Gaussian by a complex plane wave; at any orientation the temporal or spatial scale(s) to be isolated are determined by the wavelength. λ of the plane wave and the spatial resolution by the spatial aspect ratio γ, which specifies the ellipticity of the support of the Gabor function. At any orientation, both types of filter may be tuned at any frequency or spatial wavenumber by varying the length or the wavelength respectively. The filters can be applied directly to two-dimensional radargrams, in which case they abstract information about given scales at given orientations. Alternatively, they can be rotated to different orientations under adaptive control, so that they remain tuned at a given frequency or wavenumber and the resulting images can be stacked in the LS sense, so as to obtain a complete representation of the input data at a given temporal or spatial scale.In addition to isolating geometrical information for further scrutiny, the proposed filtering methods can be used to enhance the S/N ratio in a manner particularly suitable for GPR data, because the frequency response of the filters mimics the frequency characteristics of the source wavelet. Finally, signal attenuation and temporal localization are closely associated: low attenuation interfaces tend to produce reflections rich in high frequencies and fine-scale localization as a function of time. Conversely, high attenuation interfaces will produce reflections rich in low frequencies and broad localization. Accordingly, the temporal localization characteristics of the filters may be exploited to investigate the characteristics of signal propagation (hence material properties). The method is shown to be very effective in extracting fine to coarse scale information from noisy data and is demonstrated with applications to noisy GPR data from archaeometric and geotechnical surveys.     

A FREQUENCY DOMAIN MAPPING APPROXIMATION OF MOVEOUT FILTERS WITH APPLICATIONS*

Wavenumber aliasing is the main limitation of conventional optimum least-squares linear moveout filters: it prevents adequate reject domain weighting for efficient coherent noise rejection. A general frequency domain multichannel filter design technique based on a one-to-one mapping method between two-dimensional (2D) space and one-dimensional (1D) space is presented. The 2D desired response is mapped to the 1D frequency axis after a suitable sorting of the coefficients. A min-max or Tchebycheff approximation to the desired response is obtained in the 1D frequency domain and mapped back to the 2D frequency domain. The algorithm is suitable for multiband 2D filter design. No aliasing damage is inherent in the linear moveout filters designed using this technique because the approximation is done in the frequency-wavenumber (f, k)-domain. Linear moveout filters designed by using the present coefficient mapping technique achieve better pass domain approximations than the corresponding conventional least-squares filters. Compatible reject domain approximations can be obtained from suitable mappings of the origin coefficient of the desired (f k)-response to the 1D frequency axis. The (fk)-responses of linear moveout filters designed by using the new technique show equi-ripple behavior. Synthetic and real data applications show that the present technique is superior to the optimum least-squares filters and straight stacking in recovering and enhancing the signal events with relatively high residual statics. Their outputs also show higher resolution than those of the optimum least-squares filters.     

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.     

LEAST-SQUARES FILTERS WITHOUT TRANSIENT ERRORS: AN EXAMINATION OF THE ERRORS IN LEAST-SQUARES FILTER DESIGN*

A new approach has been developed for the design of cross-equalization filters by the least-squares method. The filters estimated by this new exact method are subject to only two types of error: bias and random error. Cross-equalization filters estimated by a more conventional least-squares method are further subject to “transient error”. This type of error becomes important when designing filters from a data gate of a length comparable with the length of the filter, i.e., less than four times the length of the filter. The effect of altering various design parameters has been investigated for the new method. It has been found that the proportion of bias in the filter decreases as the effective filter length increases, whereas the random error in the filter decreases with increase in either the signal-to-noise ratio of the data or the ratio of the data duration to the filter length. The level of whitening applied to the auto-correlation matrix before inversion was not found to be a critical design parameter. Also, two techniques have been tested for reducing any anomalous d.c. component in the calculated filter.     

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.     

The dynamic response of a fluid‐filled borehole under a normal point load on the free surface

The dynamic response of a semi‐infinite fluid‐filled borehole embedded in an elastic half‐space under a concentrated normal surface load is analysed in the long‐wavelength limit. The solution of the problem is obtained with integral transforms in the form of a double integral with respect to the slowness and frequency. The partial P‐ and SV‐wave responses are further transformed to path integrals along Cagniard paths in the complex slowness plane. Unlike the traditional Cagniard‐de Hoop technique based on the Laplace transform of time dependence, this paper is based on the Fourier transform. The tube‐wave response is presented as a causal integral over a slowness range. The resultant representation in the time‐domain is suitable for the numerical evaluation of the complete response in the fluid‐filled borehole, especially at large distances. Asymptotic analysis of seismic phases arising in the borehole is performed on the basis of the obtained solution. The complete asymptotic wavefield consists in P‐ and SV‐waves, the Rayleigh wave and the low‐frequency Stoneley (tube) wave. Pressure synthetics obtained by the use of the asymptotic formulas are shown to be in good agreement with straightforward calculations.     

自适应空间分区裂步傅立叶偏移方法

本文提出了一种新的偏移方法——自适应空间分区裂步傅立叶（ASDSSF）偏移方法。该方法将剥层相位移方法的思想推广到裂步傅立叶偏移方法，使之当速度场出现强间断时也能精确而有效地成像。原理上ASDSSF偏移属于多参考慢度（MRS）偏移方法，本文的重点是，在不损失精度的同时选取比同类MRS偏移方法更少的参考慢度。我们根据全局速度函数的变化和误差控制参数来选择参考慢度，同一个参考慢度所对应的速度构成一个分区，每一个空间分区可以由几个空间上不连续的子分区组成，从而有效地减小了参考慢度的个数。每一延拓步的参考慢度以及参考慢度的个数和如何构建分区都是根据速度函数自动生成，因此更为合理。为了消除速度场强间断产生的人为噪音，设计了简单有效的f-k域的光滑滤波。我们对一个生成的二维叠前模型和SEG／EAEG盐丘模型进行了试算。     

A Scheme for Initial Beam Deployment for the International Monitoring System Arrays

v--vThe International Monitoring System (IMS) includes a diverse set of seismic arrays with different configurations. These configurations have apertures ranging from less than 1 to more than 25 km and minimum interelement spacings varying from 0.1 to 3.6 km. This paper presents a scheme for initial beam deployment for this variety of seismic arrays. Beamforming is equivalent to a spatiotemporal bandpass filter of which passband is defined by the minimum and maximum wavenumbers, which are functions of the geometry configuration of the array. Deployment of steered-beams for signal detection is based on the wavenumber resolution of the array, slowness and frequency distributions of seismic phases, and coherence properties of seismic signals and noises among sensors. Within the wavenumber passband, all possible slowness values are determined by the resolution for each frequency band, and those that are outside the range of seismological interest are excluded. The appropriate azimuthal distribution for each selected slowness is determined from the azimuthal resolution. Using this approach, detection beams for each array are rationally deployed in the slowness-azimuth and frequency domain.     

ω-k FILTER DESIGN*

Two-dimensional band-pass filters can be constructed by a simple extension of the theory of one-dimensional band-pass filters. Similarly to the one-dimensional analogue the shape of the two-dimensional filter is important in determining its effectiveness. The band-pass filter formulation can be further refined so that the filter will concentrate its rejection energies in certain areas of the ω, k plane. Such band-pass, band-reject filters are found by solving a set of simultaneous equations.     

Filter formulation and wavefield separation of cross-well seismic data

Multichannel filtering to obtain wavefield separation has been used in seismic processing for decades and has become an essential component in VSP and cross-well reflection imaging. The need for good multichannel wavefield separation filters is acute in borehole seismic imaging techniques such as VSP and cross-well reflection imaging, where strong interfering arrivals such as tube waves, shear conversions, multiples, direct arrivals and guided waves can overlap temporally with desired arrivals. We investigate the effects of preprocessing (alignment and equalization) on the quality of cross-well reflection imaging wavefield separation and we show that the choice of the multichannel filter and filter parameters is critical to the wavefield separation of cross-well data (median filters, f–k pie-slice filters, eigenvector filters). We show that spatial aliasing creates situations where the application of purely spatial filters (median filters) will create notches in the frequency spectrum of the desired reflection arrival. Eigenvector filters allow us to work past the limits of aliasing, but these kinds of filter are strongly dependent on the ratio of undesired to desired signal amplitude. On the basis of these observations, we developed a new type of multichannel filter that combined the best characteristics of spatial filters and eigenvector filters. We call this filter a ‘constrained eigenvector filter’. We use two real data sets of cross-well seismic experiments with small and large well spacing to evaluate the effects of these factors on the quality of cross-well wavefield separation. We apply median filters, f–k pie-slice filters and constrained eigenvector filters in multiple domains available for these data sets (common-source, common-receiver, common-offset and common-midpoint gathers). We show that the results of applying the constrained eigenvector filter to the entire cross-well data set are superior to both the spatial and standard eigenvector filter results.     

ω-k FILTER DESIGN*

Two-dimensional band-pass filters can be constructed by a simple extension of the theory of one-dimensional band-pass filters. Similarly to the one-dimensional analogue the shape of the two-dimensional filter is important in determining its effectiveness. The band-pass filter formulation can be further refined so that the filter will concentrate its rejection energies in certain areas of the ω, k plane. Such band-pass, band-reject filters are found by solving a set of simultaneous equations.     

地震台阵上信号方位角和慢度的时、频域估计方法比较

地震台阵常用频率-波数分析法来估计地震信号的后方位角和慢度。尽管有关的算法是众所周知的,但不同的实现在某些情况下可能导致不同的结果。如应用在宽频带的澳大利亚WRA台阵的记录时,标准的f-k分析方法往往给出不正确的结果。本文发现错误是由FFT的频谱泄漏效应引起的,如在进行FFT之前对原始数据进行高通滤波则可以有效地避免这样的错误。进一步对在时域中直接估算信号方位角和慢度的一种算法进行了分析,比较发现时域方法的计算速度和精度都不低于频域方法,且在某些低信噪比的情况下,前者可以给出更可靠的结果。     

Ground roll attenuation using a time-frequency dependent polarization filter based on the S transform

The ground roll and body wave usually show significant differences in arrival time, frequency content, and polarization characteristics, and conventional polarization filters that operate in either the time or frequency domain cannot consider all these elements. Therefore, we have developed a time-frequency dependent polarization filter based on the S transform to attenuate the ground roll in seismic records. Our approach adopts the complex coefficients of the S transform of the multi-component seismic data to estimate the local polarization attributes and utilizes the estimated attributes to construct the filter function. In this study, we select the S transform to design this polarization filter because its scalable window length can ensure the same number of cycles of a Fourier sinusoid, thereby rendering more precise estimation of local polarization attributes. The results of applying our approach in synthetic and real data examples demonstrate that the proposed polarization filter can effectively attenuate the ground roll and successfully preserve the body wave.     

Detection of regional phases of seismic body waves using an array of three-component sensors

A small-aperture seismic array consisting of seven three-component seismometers carried out continuous measurements of regional seismicity in a selected area of the Nizhni Novgorod nuclear power plant during four months of 2013. Automatic signal detection using beamforming was applied separately for each motion component. Two horizontal components were transformed into radial and transverse components for the given values of the velocity and azimuth of the plane wave front. We have investigated the dependence of the coherence of microseismic noise on frequency, azimuth, and slowness, as well as determining the level of cross-correlation between signals on separate channels in order to estimate expected improvement in the signal-to-noise ratio, which is crucial for signal detection. Most signals detected by the seismic array from regional sources are associated with quarry blasts. Using repetitive explosions at seven quarries, we have quantitatively estimated and compared the increase in detection efficiency of regional seismic phases using a three-component small aperture seismic array and a subarray of vertical sensors. Horizontal sensors showed a higher efficiency in the detection of transverse waves, while the subarray of vertical sensors missed S-waves from certain events. For one of the nearby quarries, the vertical subarray missed up to 25% of events (5 of 20). The results of the investigation point to the need for the use of three-component seismic arrays for the study of regional seismicity.     

REALIZATION OF SHARP CUT-OFF FREQUENCY CHARACTERISTICS ON DIGITAL COMPUTERS*

Sharp cut-off frequency filtering is carried out in the discrete time domain on digital computers. A convolution of the digital filter impulse response with the sampled input yields the output. For practical reasons, the length of the filter inpulse response, corresponding to the number of filter coefficients, is limited, and consequently the resulting frequency characteristic will no longer be identical to that originally specified. This is analogous to synthesising some specified frequency characteristic with a finite number of resistive, capacitative and inductive components. In Part I of this paper, we examine the effect of approximating the sharp cut-off frequency characteristic best in a mean square sense by an impulse response of finite length. The resulting frequency characteristic corresponds to the truncated impulse response of the specified frequency characteristic. It has a cut-off slope proportional to, and a mean square error inversely proportional to, the length of the impulse response, and is a biassed odd function about the cut-off frequency point. Because of the Gibbs phenomenon for discontinuous functions, the resulting frequency characteristic will always have a maximum overshoot with respect to the specified characteristic of ± 9%, regardless of the length of the corresponding impulse response. Equal length truncated impulse responses of specified filters with different cut-off frequencies yield frequency characteristics which are almost identical about their respective cut-off points. Now on a log frequency scale (as against a linear frequency scale implied previously) such characteristics may be made almost identical about the respective cut-off points by having the truncated impulse responses composed of an equal number of zero crossings. Results for the low-pass filter are applicable to the high-pass and band-pass characteristics. In the latter case, the mean square error is double that for a single slope characteristic (low-pass or high-pass) and the slopes at both edges of the passband are approximately equal in magnitude to the length of the impulse response (linear frequency scale). Part II of this paper is concerned with reducing the ± 9% overshoot that results from the discontinuous nature of the sharp cut-off frequency characteristic and which is not dependent on the length of the truncated impulse response. The reduction is achieved, at the expense of the steepness of cut-off for the resulting frequency characteristic, by the use of functions which weight the truncated impulse response of the specified frequency characteristic. These functions are called apodising functions. Among other variables, the length of the truncated weighted impulse response will determine the amount of maximum overshoot since the effective frequency characteristic being approximated is no longer a discontinuous function. The digital realization of the finite length impulse responses of Parts I and II is discussed in Part III, together with the optimum partially specified digital filter approximation to the desired frequency characteristic.     

Body-wave dispersion: Measurement and interpretation

Generalizing previous studies on short-period data, it is shown that body-wave dispersion can be measured from broad-band records of earthquakes of moderate magnitude. The method is based on the direct measurement of the arrival time of the frequency components of a seismic wave, and the arrival time is defined by its expectation value. The frequency components of the signal are obtained through a narrow band-pass filtering process. Previous to any interpretation, a correction of the arrival time for instrument response and group delay of the filter is needed. In the first step, body-wave dispersion is related to an absorption band to account for intrinsic attenuation, and thereafter we generalize this interpretation by considering a cascade of filters to account for medium parameters (attenuation and a layered crust) and source parameters (source time function and finiteness of fault). An inversion scheme to obtain the filter parameters can be devised by following, in a formal way, the same procedure as for the case of surface wave dispersion.     

DETERMINATION OF RESISTIVITY SOUNDING FILTERS BY THE WIENER-HOPF LEAST-SQUARES METHOD*

It has been found that the Wiener-Hopf least-squares method is a very successful tool for the determination of resistivity sounding filters. The values of the individual filter coefficients differ quite appreciably from those obtained by the Ghosh procedure. These differences in the filter coefficients, however, have only a negligible effect on the output of the filter. It seems that these differences in the coefficients correspond to a filter function of a rather narrow frequency band around the Nyquist frequency, which is only very weakly present in the input and output functions.     

基于波干涉原理的套管井阵列声波处理方法及其应用

在套管井测井中,沿着钢管产生的拉伸波(即管波)在声波数据中产生很强的干扰,使得地层声速的提取十分困难.基于管波与地层声波的干涉原理,本文提出了一种声波处理方法.管波的传播特征是已知的,且管波和地层波都由同一声源激发,因此,可以用管波来改造测井声波数据.对改造后数据采用希尔伯特变换作进一步处理,得到与地层信息相关的解析信号,从该信号中即可提取出地层声速.通过对理论模拟数据的处理和现场测井数据的应用,验证了新方法的正确性和有效性.     

OPTIMUM DIGITAL FILTERING AND INVERSE FILTERING IN THE FREQUENCY DOMAIN*

Two distinct filters are developed in the frequency domain which represent an attempt to increase the resolution of fine structure contained in the signal whilst keeping the expected filtered noise energy within reasonable bounds. A parameter termed the White Noise Amplification is defined and used together with a measure of the deconvolved pulse width in order to provide a more complete characterisation of the filters. Each of the two main types of frequency domain filters discussed varies in properties with respect to a single adjustable parameter. This may be contrasted with a time domain Wiener filter which in general has three variables: length, delay and an adjustable noise parameter or weight. The direct frequency domain analogue of the Wiener filter is termed a gamma-Fourier filter, and is shown to have properties which span the range from those of a spiking filter with zero least square error at one extreme, to those of a matched filter at the other extreme of its variable parameter's range. The second type of filter considered—termed the modulated Gaussian filter—is similarly shown to be a perfect spiking filter at one extreme of its parameter range, but adopts the properties of an output energy filter at the other extreme.     

Coherence methods in mapping AVO anomalies and predicting P-wave and S-wave impedances

Filters for migrated offset substacks are designed by partial coherence analysis to predict ‘normal’ amplitude variation with offset (AVO) in an anomaly free area. The same prediction filters generate localized prediction errors when applied in an AVO‐anomalous interval. These prediction errors are quantitatively related to the AVO gradient anomalies in a background that is related to the minimum AVO anomaly detectable from the data. The prediction‐error section is thus used to define a reliability threshold for the identification of AVO anomalies. Coherence analysis also enables quality control of AVO analysis and inversion. For example, predictions that are non‐localized and/or do not show structural conformity may indicate spatial variations in amplitude–offset scaling, seismic wavelet or signal‐to‐noise (S/N) ratio content. Scaling and waveform variations can be identified from inspection of the prediction filters and their frequency responses. S/N ratios can be estimated via multiple coherence analysis. AVO inversion of seismic data is unstable if not constrained. However, the use of a constraint on the estimated parameters has the undesirable effect of introducing biases into the inverted results: an additional bias‐correction step is then needed to retrieve unbiased results. An alternative form of AVO inversion that avoids additional corrections is proposed. This inversion is also fast as it inverts only AVO anomalies. A spectral coherence matching technique is employed to transform a zero‐offset extrapolation or near‐offset substack into P‐wave impedance. The same technique is applied to the prediction‐error section obtained by means of partial coherence, in order to estimate S‐wave velocity to P‐wave velocity (V_S/V_P) ratios. Both techniques assume that accurate well ties, reliable density measurements and P‐wave and S‐wave velocity logs are available, and that impedance contrasts are not too strong. A full Zoeppritz inversion is required when impedance contrasts that are too high are encountered. An added assumption is made for the inversion to the V_S/V_P ratio, i.e. the Gassmann fluid‐substitution theory is valid within the reservoir area. One synthetic example and one real North Sea in‐line survey illustrate the application of the two coherence methods.