fkkstrip: [Kx,Ky] Filtering
by Paul G. Garossino June, 1993
Introduction
Recent work on 3D Post Stack DeMult [psdm3d] made necessary a filtering algorithm to remove or enhance flat-lying events from a cube of data. The resulting routine, fkkstrip, has been found to also be useful as a general [Kx,Ky] filter.
The Algorithm
At the heart of the routine is a radial implementation of the Bessel Jo(r) function. The filter is designed in [x,y] and applied through a two dimensional convolution in the [x,y,t] domain. An optional radial cosine taper may be applied to the filter. Optional tapering and padding may also be applied to the input data independently in both the x and y directions.
Program output is either [input - input * filter] or optionally [input * filter].
Data Example
A high cut filter has been applied to a post stack 3D dataset from the Younis project in the Gulf of Suez. The data arrived in LI [Line Index] sorted order (Fig. 1) and was converted to time [TI] slice sorted order (Fig. 2) prior to filtering using:
ttds3d -N LIslice -ND txy -OD xyt -O timeSlice
The data was then filtered (Fig. 3) using:
fkkstrip -N timeSlice -O timeSliceFkkstrip -l 3 -pass
Subtraction of the output from the input results in a difference plot (Fig. 4) exposing the data removed during filtering. The effects of the filter are made obvious when the data are viewed in the Fourier domain. This was accomplished using:
The output records (Figs. 5 & 6) display the 2 dimensional spatial Fourier transforms of the input data both pre and post filter application. The axes are labeled in samples and extend from the positive to negative Nyquist frequency in both the Kx and Ky directions. The data of interest is located between the origin and 0.5 Nyquist in both orientations.
The post fkkstrip spectrum (Fig. 6) shows the level of attenuation achieved.
The filtering may have been a little too severe as what appears to be a coherent event [marked by the arrows in Fig. 5] has been attenuated. This energy corresponds to high frequency events on the input TI slice. Though probably not optimum, for this example no further parameterization of the filter was attempted. In production use it is expected that the processor would quality control the filter design so as not to remove any primary reflections of interest.
Should a time variant implementation be required one could use a combination of fftxy and polymute by picking a suite of polygons to apply to the transformed TI sorted data. In this example however a single fkkstrip time invariant filter was applied to the entire TI volume.
The DI Sort
To return the data to the input LI presentation it is often interesting to first construct a DI sort using:
timeSliceFkkstrip -ND xyt -OD tyx -OD IsliceFkkstrip
The dramatic reduction in noise seen on this presentation (Fig's 7 & 8) is quite remarkable when considering the minimal smoothing observed on the time slices. The difference plot (Fig. 9) reveals that mainly random noise, but also some coherent energy, had been removed. This coherent energy corresponds to events that dip steeply in time and thus land within the reject frequency range on the TI sorted data.
The LI Sort
The LI sorted volume was then generated using:
ttds3d -N DIsliceFkkstrip -ND xyt -OD txy -O LIslice Fkkstrip
Comparison of the pre and post fkkstrip LI slice data
(Fig's 1 & 10) reveals that the TI slice processing has been very effective at attenuating both random noise and steeply dipping events. Some attenuation of energy at shallower dips is also evident especially on the difference plot (Fig. 11). This was disconcerting, however, an analysis of the DI slice data (Fig's 7, 8 & 9) revealed that these events are associated with steeply dipping arrivals in the crossline direction which fall outside the pass band of the filter applied to the TI slice data.
The filter may have been more effective if convolved with the original dataset as the original crossline sample interval was 25 meters [the processed volume was sampled using a 50 meters crossline sample interval]. It would also have been instructive so see the fourier spectrum of the data prior to decimation in order to determine if an anti-alias filter was required in advance of resampling.
A comparison of the Fourier transforms of the LI records before and after filtering (Fig's 12 &13) shows that the random noise level has been lowered not just at high wavenumbers [as expected], but also over the pass band of the filter. This both explains the unexpected improvement in signal continuity and illustrates the power of time slice processing.
Filter Parameterization
The command line entries affecting the rejection characteristics of the filter are:
-l the filter length. This parameter specifies the number of points [N], on one side of the filter. The filter used will contain N2 number of points. This parameter greatly affects the spacial rejection characteristics of the filter. The greater the number of points in the filter the lower the spatial cutoff frequency.
-i the filter intercept. This parameter controls the number of lobes of the Bessel function used in the filter. A value less than 0.35 is recommended. Values greater that this may be used but could result in an aliased filter design.
-ftaper the percent of the filter edge to multiply by a cosine taper. The greater this value the lower the resulting sidelobes on the frequency response of the filter.
The command line entries affecting data conditioning prior to filtering are:
-xtaper,-ytaper the x and y edge tapers. These entries are the number of samples at the edge of the input TI slice over which to use a cosine taper.
-xpad,-ypad the x and y edge pads. These entries specify the number of samples to use as a pad around the input data.
Parameterization Examples
To illustrate the effect of filter length the following was run:
in -O out -l 7
fkkstrip -N in
-O out
-l 21
fkkstrip -N in
-O out
-l 81
The 2-d amplitude spectrum (Fig's 14, 15 & 16) of the output data illustrates that as the filter length is increased the spatial cutoff frequency is reduced.
To demonstrate the effect of the intercept parameter the following was run:
in -O out -i 0.005 -l 31
fkkstrip -N in
-O out
-i 0.20
-l 31
In this instance the resulting amplitude spectrums (Fig's 17 & 18) reveal that as the intercept is increased the depth of the first sidelobe in [Kx,Ky] is decreased.
To illustrate the effect of the cosine taper the following was run:
fkkstrip -N in -O out -ftaper 5
fkkstrip -N in
-O out
-ftaper 50
Without a taper there is no guarantee that the filter edges are anywhere close to zero amplitude. As the taper width is increased the filter becomes more stable [the amount of filter ringing in (Kx,Ky) is reduced] (Fig's 19 & 20).
Pitfalls
1. There is no provision in the program guarding against aliasing. It is possible to chose an intercept value that results in an aliased filter. For instance:
fkkstrip -N in -O out -i 10.0
A quick look at the 2-d amplitude spectrum of the output (Fig. 21) will verify that this has occurred.
2. Since fdslice has no way of remembering the original data indexing it is necessary to run hdrswap on the final output accessing the headers of the original input dataset. It is also necessary to restore the line header entry for sample interval SmpInt [using utop] as it will be set to unity by fdslice.
Runtime Considerations
The Younis dataset processed contained 252 Mbytes of data and required 35 minutes [sparc 10] for the application of the filter. The various sorting runs took an additional 90 minutes. The greater the filter length specified the longer the run time as the number of points in the filter is the square of this length. Also the larger the pad used the longer the run time though this has a second order effect when compared to the effects of the filter length.
Miscellaneous
The filter has been applied to conventional 2 dimensional seismic data by treating the prestack records as a post stack 3 dimensional volume. In the experiment the data were moveout corrected, sorted to common mid point and time sliced. The results of the filtering showed promise in this type of usage for dipping event rejection, however, more work is required to understand all the implications of this type of filtering.
The filter has also been applied to map data in order to provide regional / residual identification. This approach has met with great success.
This field of image processing of geophysical data is quite young although image processing itself has a long history in the medical, EE and space sciences fields. If you can think of some other off the wall use for this type of filter give us a call and we will endeavour to implement your algorithm in a timely fashion.
Return to the USP Home Page.
Frame Doc to html conversion by Tony
Garossino
