User:ChaseBilleaudeau/sandbox

User:>/Sandbox

F-k filtering is a 2-D data processing technique that is based on a relationship between two domains: a time-space (t-x) domain and its Fourier dual, an f-k domain. This technique filters coherent noise for wavefield separation to produce better quality data. In short, f-k filtering relies on signal and noise separation via the 2-D Fourier transform. The region in the f-k domain where noise is present is set to be zero and the signal region is set to be one. Then the filtered shot record is restored by the inverse Fourier transform.

Introduction
Geophysicists often conduct seismic reflection surveys in order to understand the subsurface geology of a given area. The method of acquisition may vary, but in general, an array of geophones or hydrophones are arranged linearly in order to collect the waves generated by a source. The shot point sends multiple waves through the lithology, and the waves become reflected, refracted, or dispersed. Typically, waves that may be encountered include: P-waves, S-waves, Rayleigh waves, Love waves, and Stoneley waves. Filtering techniques have been developed to remove coherent noise; thus, imaging is improved.

Inventions and Important Discoveries
The first time interest was shown in propagating waves beneath Earth’s surface dates back to as early as 132 BC in China, with the invention of the seismoscope. Centuries later, Filippo Cecchi invented the seismometer, which was the first time-recording seismograph. Seismic recording instruments continued to evolve, and as a result, imaging improved. In 1906, the Earth’s core was discovered through seismic imagery by Richard Oldham. Three years later, another important discovery took place – the Moho discontinuity.

Computer Age
In the 1960's, Seismology began its movement from needle-and-paper recordings to recording and processing through computers. Along with the inception of computers, came the introduction of filtering techniques. These filtering techniques were rapidly evolving, and many of them are still used today, such as f-k filtering. f-k filtering was first introduced to the seismic industry through papers published by Embree et al. (1963) and by Fail and Grau (1963). With the advancement of technology, in 1967, 3D seismic was brought to industry. ExxonMobil conducted the first 3D seismic survey over the Friendswood field near Houston, Texas. Processing such a large amount of input traces, roughly a half million, took over 2 years for the technology at that time. Most recently, 4D, or time-lapse, seismic has been developed and is widely used in the oil and gas industry. As survey methods continue to rapidly evolve, processing methods must keep pace.

Mathematical and Physical Explanations
We would introduce two design methods of an f–k filter to suppress an incident wave whose transmitted position and received position are known.

Method I

The basis for f-k filtering lies in understanding the relationship between phase velocity vp, frequency f, and wavelength &lambda;:


 * $$v = \lambda f $$.

Frequency is related to wavenumber by the following relation:


 * $$v = \frac{f}{k}$$

The picture to the right demonstrates a wavefront approaching an array of geophones spaced equal distance apart:

The wavefront approaches the array of geophones with the apparent velocity:


 * $$v_\alpha = \frac {v} {\sin_\alpha}$$ where $$v_\alpha \ge v $$

By the above equations, we can get the formula for apparent wavelength:


 * $$\lambda_\alpha = \frac {v_\alpha} {f}$$.

And by the relationship between frequency and wavenumber, we can get the formula for the apparent wavenumber:


 * $$k_\alpha = \frac {f} {v_\alpha} $$



The figure above shows four velocities (V1, V2, V3, V4). Between V2 and V3 lies a range of coherent noise that should be eliminated to improve data quality.

In order to eliminate this noise, we take the 2-D signal, represented as F(&omega;, kα), and apply an f-k filter represented by the function H(&omega;, kα) to it to get a filtered spectrum G(&omega;, kα):


 * $$G(\omega,k_\alpha) = F(\omega,k_\alpha) \cdot H(\omega,k_\alpha)$$

where ω is the angular frequency and the f-k filter is given by the function


 * $$H(\omega,k_\alpha) = \begin{cases} 0(v_{2} \le \frac {\nu} {k_\alpha} \le v_{3}) \\ 1(0 \le \frac {\nu} {k_\alpha} \le v_{1}) \end{cases}$$

When apply 2-D inverse Fourier transform to G(ω, kα), we get a time-spatial radar profile without the incident wave component as follows:


 * $$g(t,x_r) = \int \int G(\omega,k_x) \cdot e^{j(\omega t + k_x x_r)} d \omega d k_x$$.

Method II

In this method, we apply a time shift so that the incidence would arrive at different geophones at the same time. Thus, the apparent velocity of the incident wave grows infinitely large, and the f–k spectrum of the incident wave becomes a spectrum whose slope is infinitely large. The following is a mathematical explanation of the implementation.

We define the nearest receiving point, which is represented as (xrn, zr), to the source to be our reference point. Then, the propagation time between the source (xt, zt) and our reference point (xrn, zr) is given by the function
 * $$t_{ref} = \frac {\sqrt {(x_t - x_{rn})^2 + (z_t - z_r)^2}} {v}$$.

The following function presents an arrival time of the incidence using the position of the receiver:
 * $$t_(x_r) = \frac {\sqrt {(x_t - x_r)^2 + (z_t - z_r)^2}} {v}$$.

Then, differences of the arrival time which should be negated are calculated by
 * $$T(x_r) = t(x_r) - t_{ref}$$.

Using the above value. we can describe the time shift as
 * $$F'(\omega,x_r) = F(\omega,k_x) \cdot e^{+j \omega t + k_x x_r} d \omega d k_x$$.

The frequency-spatial signal applied to the time shift is substituted to the F(ω, xr) of the function
 * $$G(\omega,k_\alpha) = F(\omega,k_\alpha) \cdot H(\omega,k_\alpha)$$.

Then we can get an f–k spectrum. In the f–k spectrum, the spectrum component which should be suppressed is distributed on the frequency axis.

Fig. 1. shows the data obtained through a bistatic radar system by hitting a metallic sphere in a free space. These are the data before the time shift is applied. (a)Radar profile. (b)f–k spectrum. Fig. 2. shows an example, which is obtained by applying the time shift under the concept of this method to the data set shown in Fig. 1. In Fig. 2(a), the direct wave arrive at the receivers at the same time. And its apparent velocity becomes infinitely large. As a result, a spectrum component of the direct wave in the f–k domain is distributed along the frequency axis, spatial frequency k = 0, as shown in Fig. 2(b). Then, we apply a filter to reject a dc component in a spatial frequency direction. This filter in the f-k domain is given by
 * $$H(\omega,k_\alpha) = \begin{cases} 0(k_\alpha=0) \\ 1(k_\alpha \ne 0) \end{cases}$$.

Then a reverse time shift after the filtering is needed.

Compared to the filter in Method I, the 2-D filter in Method II is more advanced. It does not change with the incident angle. Therefore, a suppression of a spectrum other than the direct wave should not change even when the incident angle is very small.

Vertical Seismic Profiling (VSP)
In VSP, f-k filtering is used to separate upgoing and downgoing waves in the time-depth domain. This is possible because the propagation direction of the waves will affect their assigned value, being either positive or negative (N. Hayashi 2010). Essentially, when the upgoing and downgoing waves are separated, one wave spectrum is suppressed while the other is preserved.


 * $$V_{dir}*V_{ref}<0$$

Vertical Radar Profiling (VRP)
Unlike VSP, VRP implements f-k filtering in the f-k domain. Here, f-k filtering is applied to the separation of direct and reflected waves. In the f-k domain, direct and reflected waves have opposite values.

Separation and Enhancement of Split S-waves
S-waves can be separated from multicomponent shot records and enhanced through f-k filtering. The S-wave events are separated by using their differential moveout, or the time waves take to travel between each trace. Typically, the delay between split wavelets is between 20 to 40 ms for 10 HZ wavelets (M. Boulfoul and D.R. Watts 1994). The lag can be maximized by using the rotation of the receivers to find the angle that gives a global maximum crosscorrelation coefficient between two S-waves (M. Boulfoul and D.R. Watts 1994).

Ground-roll attenuation
In land-based seismic surveys, ground roll, a type of surface wave, is very problematical. Ground roll is generated by low-velocity, low-frequency, high-amplitude Rayleigh waves. This surface wave can be suppressed by applying an f-k filter in the frequency-wavenumber domain. The ground roll is mapped as lines in the f-k domain and can be removed by using a simple 2D band-pass filter (Porsani et al. 2010).

Limitations
f-k filtering can sometimes produce undesirable results when choosing a reject band for the filtering of seismic data. For example, significant signal distortion and/or the introduction of filtering artifacts may occur (March and Bailey, 1983).

Signal Distortion
Signal distortion is mainly caused by the overlapping of the reject band of the f-k filter with the signal region of the f-k domain. This results in images with “worm-like” appearances. To understand how this occurs, consider the effects applying the Fourier transform to the spike of a signal (Duncan 1994). When the Fourier transform is applied to a spike of a signal, a flat frequency spectrum is produced. This spike can be distorted if some of the frequencies are set to zero and the inverse Fourier transform is applied (Duncan 1994). Signal distortion is proportional to the size of the reject zone, and its effects can be reduced by choosing a velocity filter with a narrow reject band in the f-k domain.

Artifacts
Artifacts are attributed to the Fourier-based method of f-k filtering. Essentially, the occurrence of artifacts is due to the band-limited nature of the filtering process. In order to limit the interference of these filter-created artifacts, a wide-pass band filter approach can be taken.

Summary of Limitations
The f-k reject band can be widened or narrowed, but the aforementioned effects may occur. By widening the reject band, artifacts in the dataset can be reduced, but this negatively impacts the stratigraphic and structural detail due to overlapping of the signal and noise in the f-k domain. If the reject band is narrowed, real data may be attenuated. For best results, a balance must be found, and the width of the reject band will largely depend on the scope of work for the seismologist.

Introduction
O. C. Ogiesoba et al. (2011) conducted a seismic study in the Pampo Field, Campos Basin, offshore Brazil. The area of study covered roughly 36 sq. kilometers. Numerous carbonate systems run throughout the Campos Basin, and these systems pose quite a challenge for seismic interpretation. The geologic boundaries can be difficult to find due to the varying nature of the diagenetic processes that control the porosity and permeability of these rock units (Ogiesoba et. al. 2011). Not only did the carbonate systems mixed with the siliciclastic sediments pose a challenge, but so did the linear and random noise created by the various features within the basin.

Applications and Processes
In order to remove some of these effects, f-k and τ-ρ filtering techniques were chosen because of their effectiveness and economic viability. f-k filtering was used to attenuate high-angle coherent noise, and τ-ρ filtering was used to improve event continuity.

When performing f-k filtering, all seismic data undergoes a transformation — horizontal events in the x-t (offset-time) domain transform to vertical events, and vertical shifts to horizontal. Similar apparent velocities are placed into a single linear event in the f-k domain (Ogiesoba et. al. 2011). From the f-k domain, any of these events can be muted or removed.

O. C. Ogiesoba et. al. (2011) went through the seismic data in a step-by-step process in order to apply the f-k filter.

Step 1: The entire 3D dataset was checked at set intervals for linear noise, and the associated dips were computed in milliseconds per trace.

Step 2: The Nyquist wavenumber was computed in cycles per kilometer according to:


 * $$k_\text{N} = \frac{1}{2\Delta x}$$

Step 3: The Nyquist frequency was computed in cycles per second according to:


 * $$f_\text{N} = \frac{1}{2\Delta t}$$

Step 4: The filter was applied based on the parameters from steps 1 - 3.

Step 5: The output was examined, and if linear noise was still present, the dip was recomputed and the filter reapplied.

High-angled noise was steeply dipping in the crossline direction, and consequently, there was approximately zero dip in the inline direction; therefore, the inline direction was disregarded, and the crossline was used to process the Pampo field data, with an observed noise dip range at 18 ms/trace. With trace spacing intervals at 25m and temporal sampling at 4ms, the Nyquist wavenumber and Nyquist frequency were calculated at 20 cycles/km and 125 Hz, respectively (Ogiesoba et. al. 2011). The methodology in tackling the field data processing involved the following procedures:

Step 1: A portion of the data was selected to be tested in a set of parameters.

Step 2: The data underwent a series of testing, applying different minimum and maximum dips.

Step 3: In LandmarkTM, the “Reject” option was selected to reject dips outside of the following parameters: ±18, ±10, ±8, ±6, ±6/-16, and ±3/-8 ms/trace.

Step 4: The original data was subtracted by the rejected noise to obtain the newly filtered data.

Conclusion
The best result came from the ±6 ms/trace dip filtering. Much of the coherent noise was removed without sacrificing too much real data. τ-ρ filtering was then applied to further improve the seismic imagery by improving the event continuity.

Introduction
The main objective of the Bayou Corne reflection seismic surveys is to locate aquifer sands which may be charged with free gas. Two sites, sites 1 and 2, have been selected to undergo a series of seismic testing (Fig 1). The gas content at these sites is largely unknown. Shear wave data, acquired from previous surveys, indicate there are two prominent reflector bodies. The first reflector body is fairly shallow, lying at 7-10m below ground surface (bgs). The second reflector body is at a much greater depth. It is located 40-50m bgs and possibly corresponds to the top of the aquifer sands.

Data Collection
The Bayou Corne dataset was acquired on September 22, 2012 from land seismic reflection surveys in Assumption Parish, Louisiana. The seismic reflection surveys were conducted mechanically, with the use of a hammer and steel plate. An array of geophones was placed linearly from west to east spaced one meter apart. The source was initially placed one meter from the last geophone, then subsequently moved twenty-four meters eastward after every succession of recording. The data can be found in the directory ‘/home/refseis12/BayouCorne/seismics/data/092212/H/1/su’ located in the Louisiana State University seismic reflection laboratory on machine LGC10 with network address ‘lgc10.geol.lsu.edu.’

Program
dipfilter.sh is a shell script that attenuates coherent noise from data. It contains six individual programs, each performing separate tasks.

Justification
First, a working directory with the dataset must be placed into the code. To perform this action, '/home/cbille1/BayouCorne/seismics/data/092212/H/1/su' is set to the directory, and data is extracted from the file SH_diff, containing files 1009 to 1024. The horizontal and vertical units are then normalized. This is achieved by setting the vertical (dt) and horizontal (dx) units to 1. tracf in ‘susort <$DATA_DIR/$file.su tracf’ is the header value that is used to reorder the traces. The reordered file now has traces ordered in an increasing value of tracf. Only part of the dataset needs to be used, so ‘suwind tmin=0.0 tmax=2.0’ is applied to take all samples between 0 to 2 seconds and window those samples out. ‘sugain agc=1 wagc=0.2’ adjusts for changes in signal strength in time along individual traces. This adjustment leads to better image quality.

To further increase image quality, coherent noise is removed by finding the central slope of the Rayleigh wave. The reflector of a given slope is killed by creating a band pass filter. The central slope is rejected by setting the amp values to ‘amp=1,0,0,1.’ Two slopes on either side of the central slope are chosen to be killed. Also chosen, are two additional slopes where the effects of dip filtering are stopped. In this case, ‘slopes=3,14,21,42’ is applied. After sorting the traces and applying ‘sugain agc=1 wagc=0.2,’ the dip filter is applied by filtering out the rejected range of slopes according to the code ‘sudipfilt dt=$dt dx=$dx slopes=$slopes amps=$amps bias=$bias.’ suximage displays the newly filtered data, with the clip nulling out any traces above the amplitude number assigned.

'''The windows of the 6 programs are adjusted to specified heights and widths and are placed in a non-overlapping 3x2 format for easily accessible viewing. The directory is set to use files from SH_diff.'''

This program normalizes the horizontal and vertical units, and the data is plotted without any dip filtering applied.

'''In this program, the aim is to remove coherent noise by finding the central slope of the Rayleigh wave. Once a central slope is found, in this case it is approximately 18 ms/trace, two additional slopes expanding in either direction are selected. Another set of slopes are then chosen to stop the effects of dip filtering. In order to reject the dips surrounding the central slope, the amplitudes are set to ‘amps=1,0,0,1’.'''

This program, similar to the previous example, attenuates noise, but the noise being attenuated is negative.

An image of the F-K spectrum is produced before dip filtering is applied.

An image of the F-K spectrum is produced after dip filtering is applied.

An image of the F-K spectrum is produced after negative dip filtering is applied.

Results




Conclusion
The application of the dipfilter.sh program substantially improves image quality. The Rayleigh wave appears to have a slope at approximately 18m/s. What is found, is that when the slope is set to  ‘slope=3, 14, 21, 42’, there is minimal loss to real data, and much of the noise is attenuated.

f-k-k
f-k-k can be applied to both the conventional geometry and the crosshole geometry with source depth and receiver depth as the independent spatial variables. This filtering contains less coherent noise, and consequently, shows improved continuity of reflectors. It is a well-designed and practical technique for wavefield separation in crosshole seismic reflection processing.

t-f-k
The t-f-k filtering method is used to localize the apparent velocity panel of a seismic record in a time domain. This method allows for better filtering by creating new reject zones in varying time domains. Creating an exact reject zone, by use of f-k filtering, for highly scattered waves produced by ground roll presents a difficult issue. By adding a time domain (t-f-k), a more refined reject zone can be constructed to remove excess noise (Askari and Siahkoohi 2008).

x-f-k
The x-f-k filtering method, like the t-f-k method, enhances seismic imaging. Instead of using a panel of time, a panel of space is added to filter apparent velocities in a varying space domain (Askari and Siahkoohi 2008).