# SEISGAMA: A Free C# Based Seismic Data Processing Software Platform

Authors: Theodosius Marwan Irnaka; Wahyudi Wahyudi; Eddy Hartantyo; Adien Akhmad Mufaqih; Ade Anggraini; Wiwit Suryanto
Journal: International Journal of Geophysics
 (2018)

---

## Abstract

Seismic reflection is one of the most popular methods in geophysical prospecting. Nevertheless, obtaining high resolution and accurate results requires a sophisticated processing stage. There are many open-source seismic reflection data processing software programs available; however, they often use a high-level programming language that decreases its overall performance, lacks intuitive user-interfaces, and is limited to a small set of tasks. These shortcomings reveal the need to develop new software using a programming language that is natively supported by Windows® operating systems, which uses a relatively medium-level programming language (such as C#) and can be enhanced by an intuitive user interface. SEISGAMA was designed to address this need and employs a modular concept, where each processing group is combined into one module to ensure continuous and easy development and documentation. SEISGAMA can perform basic seismic reflection processes. This ability is very useful, especially for educational purposes or during a quality control process (in the acquisition stage). Those processes can be easily carried out by users via specific menus on SEISGAMA’s main user interface. SEISGAMA has been tested, and its results have been verified using available theoretical frameworks and by comparison to similar commercial software.


---

## Body

## 1. Introduction

Geophysics is a branch of physics that studies physical properties of the earth, such as mass density, seismic wave speed, electrical resistance, and magnetic properties [1–3] from surface measurements. These physical properties can be further used to image subsurface geological conditions of the earth. One of the most popular methods in geophysics is the seismic reflection method. This method is very sensitive to a rock’s acoustic impedance variation, which is characterized by seismic waves reflected in each interface of the subsurface layer or structure. Acoustic impedance is denoted by the product of rock density and seismic wave speed within a certain medium. Seismic reflection is very popular, especially in the oil and gas industry, because this method yields high resolution subsurface images of geological structures [4]. The seismic reflection method is also useful in engineering, environmental, geohazard, and ground-water exploration, as discussed by Steeples and Miller [5].Recently, a large number of open-source software programs have emerged that can be used to process seismic reflection data; however, each program has its own disadvantages. For example, Seismic Un∗x [6] has a really good ability to process seismic reflection data, but it has no intuitive user interface and depends on a command-line interface. Templeton and Gough [7] developed a web-based seismic data processing tool based on the Seismic Un∗x platform. Through this web service, users can fill-out the form provided by the server, and the input data will be queued on a parallel processing scheduler application. In this sense, the online service is a good concept, but the performance is limited to user interaction via a form-based input with limited resources on the server side. DSISoft, which was developed by Beaty et al. [8], is a vertical seismic profiling data processor based on the MATLAB® programming language. MATLAB is a popular as programming language because it has a library of common mathematical and engineering functions. Unfortunately, it has a large computational cost as data size increases and is thus inefficient compared to a lower level programming language (e.g., C, Fortran, and C#). A comparison between MATLAB and C, made by Andrews [9], shows significant speed differences; He determines that a C program is up to 546 times faster than a MATLAB program. A more comprehensive seismic reflection data processing implementation in MATLAB was developed by Mousa and Al-Shuhail, [10]. Another example is SEISPRHO [11], software developed in Lazarus Pascal Language, which can be used to interpret seismic reflection and is enhanced with simple processing tools. This software is mainly developed for interpretation purposes.The recent popularity of Python, another high-level programming language that is free and reliable, has many features to handle numerical problems and has high quality graphics. It is increasingly being used to develop geophysical data processing software [12, 13]. One of the advantages of Python is its ability to be installed on a server so that the application can be accessed remotely (e.g., [14]). Schwehr [15] has developed an application calledseismic-py, which provides an infrastructure for creating and managing a Python library for each seismic data format. Rizvandi et al. [16] detailed the performance of an implementation of the Prestack Kirchhoff Time Migration in Python. There is currently no complete application that is built entirely with Python (i.e., ready to use without combining the Python library to another library), as this would require specialized Python programming skills, at least the basic ability to program in Python.To the best of our knowledge, there is no other free seismic reflection data processing software program built with the C# programming language that has a complete seismic processing toolbox, can process raw data from field recordings, and includes an intuitive and appealing user interface. The C# programming language was considered due to its native compatibility with the most popular operating system, that is, Windows. Netmarketshare states that approximately 81% of desktop users have a Windows XP or newer operating system installed on their computers [17]. The C# programming language is officially supported and periodically updated by Windows and included on a  .NET (dot NET) platform [18, 19]. Developing applications natively supported on Windows operating systems ensures long-term and continuous software development. Enhanced with an intuitive user interface, SEISGAMA can easily be used by beginner or advanced users, and it is particularly useful for undergraduate geophysics students as they are learning processes. The need for reliable and free seismic processing software has long been awaited, especially since the price of oil became variable. This software may ensure that oil exploration can be resumed with lower associated costs.


## 2. Methods

SEISGAMA’s development is divided into several development sections:(1) basic data processing, (2) intermediate data processing, and (3) advanced processing. This paper reports only the basic processing aspects of reflection seismic methods, and the advanced processing aspects will be discussed separately in another paper. The basic data processor that was developed in this research consists of amplitude correction, muting, F-x and F-K domain transform, velocity analysis, normal moveout (NMO) correction, and traces stacking. These processing steps combine into one trace with a high signal to noise ratio (SNR).A more detailed explanation of the processing steps follows.


### 2.1. Amplitude Correction

Amplitude correction was performed to compensate seismic signal weakening, which decreases exponentially as a function of time due to attenuation and geometrical spreading. Campillo et al. [20] explains that this phenomenon often holds in some regional seismic phases. The phenomenon is also described by geometrical spreading [21]. Amplitude attenuation may be exploited via an exponential or oscillatory equation; therefore, the amplitude correction can be formulated as an exponential equation (1). fAGC is seismic data after amplitude correction, f(t) is the data before amplitude correction, and β is an exponential constant. This equation is an example of amplitude correction for a data-independent method. After amplitude correction, the signal has amplitude with no attenuation as the recording time is increased.(1)fAGCt=fteβt.


### 2.2. Muting

The muting process removes unused parts of the seismic recording. This process is very important because in the seismic reflection method, only the reflection signal will be considered, and other signals, such as those from direct waves, ground-roll, or refraction, must be removed. The muting process is carried out by defining the area to be removed first, and then, SEISGAMA will automatically make the amplitude of that region equal to zero.


### 2.3.F-x and F-K Transformation

Seismic data is typically represented in time-domain versus (range or) distance (t-x). The F-x transformation is carried out by using a 1D Fourier transform and is applied to the transformation for each seismic trace in the data. This will transform the time series data into frequency domain data. The results are then displayed in terms of distance (X-axis) and frequency (Y-axis). The Fourier transformation is expressed by(2)Fω=∫-∞∞fte-iωtdt=FTft.The F-K transform utilizes the 2D Fourier transformation to transform a time-domain version of each seismic trace into the frequency domain in 1/s and transform the spatial domain of data into wave number (1/m). In this transformation, seismic data in time (t) versus distance (x) will be transformed into frequency (f) versus wave number (k). The 2D Fourier transform is an extension of the 1D Fourier transformation and expressed by(3)Fω=∫-∞∞∫-∞∞ft,xe-iωte-iktdtdx=FT2Dft,x.


### 2.4. Velocity Analysis

The reflection from a horizontal reflector of a seismic wave on the record for a particular seismic source will follow a parabolic equation(4)tx=t02+xv2.These characteristics are used to estimate the wave propagation speed by assuming the velocity of the medium. Parabolic response patterns can be detected using a shot-gather-based seismic semblance calculation (NE), as noted by Mousa and Al-Shuhail [10].A seismic semblance is formulated as the quotient between stacked and prestacked energy subsequently multiplied by the number of seismic traces (7) for each time window (w). Time window parameters, with width w, are required to perform an NE calculation. Here, N is the number of samples at the time window, w; M is the number of the traces on seismic trace collections; Aij is the amplitude at time i and trace j; Es is stacked energy; and Eu is prestacked energy.(5)Es=∑i=1N∑j=1MAij2,(6)Eu=∑i=1N∑j=1MAij2,(7)NE=EsM×Eu.


### 2.5. Normal Moveout Correction

A normal moveout (NMO) correction is carried out to restore the nonzero offset seismic reflectors to zero offset. The correction is done after we obtaint0 and vNMO for each reflector during velocity analysis process. The time delay after NMO correction is defined in (9). A NMO correction is suitable for a reflector having a slope of less than 15°. When a reflector has higher slope than that threshold, another correction (Dip Move Out correction) is required.(8)ΔtNMOx=tx-t0,(9)tNMOx≈x22t0vNMO2.NMO corrections can be categorized by three classifications: (1) if the NMO velocity (vNMO) is too high, then it will have a shape of a downward convex reflector, (2) if the vNMO is correct, then the reflector will be flat on t=t0, and (3) if vNMO is too low, then the reflector will be an upward convex reflector. During this process, the NMO correction of seismic signals often causes a distortion in the frequency domain. To overcome this problem, a stretching factor (10) is a provision threshold that is required. If a signal is found to have a stretching factor that exceeds the threshold, then a muting correction will be performed.(10)SNMOx=ΔtNMOxt0.


### 2.6. Stacking

Stacking is a process of merging many seismic traces that are already in a state of zero offset. The stacking process is done by quantifying the average amplitude of all signals (11). This stacking process is expected to improve SNR (signal to noise ratio) to factor √M [10].(11)tx=∑i=1MAixM.


## 2.1. Amplitude Correction

Amplitude correction was performed to compensate seismic signal weakening, which decreases exponentially as a function of time due to attenuation and geometrical spreading. Campillo et al. [20] explains that this phenomenon often holds in some regional seismic phases. The phenomenon is also described by geometrical spreading [21]. Amplitude attenuation may be exploited via an exponential or oscillatory equation; therefore, the amplitude correction can be formulated as an exponential equation (1). fAGC is seismic data after amplitude correction, f(t) is the data before amplitude correction, and β is an exponential constant. This equation is an example of amplitude correction for a data-independent method. After amplitude correction, the signal has amplitude with no attenuation as the recording time is increased.(1)fAGCt=fteβt.


## 2.2. Muting

The muting process removes unused parts of the seismic recording. This process is very important because in the seismic reflection method, only the reflection signal will be considered, and other signals, such as those from direct waves, ground-roll, or refraction, must be removed. The muting process is carried out by defining the area to be removed first, and then, SEISGAMA will automatically make the amplitude of that region equal to zero.


## 2.3.F-x and F-K Transformation

Seismic data is typically represented in time-domain versus (range or) distance (t-x). The F-x transformation is carried out by using a 1D Fourier transform and is applied to the transformation for each seismic trace in the data. This will transform the time series data into frequency domain data. The results are then displayed in terms of distance (X-axis) and frequency (Y-axis). The Fourier transformation is expressed by(2)Fω=∫-∞∞fte-iωtdt=FTft.The F-K transform utilizes the 2D Fourier transformation to transform a time-domain version of each seismic trace into the frequency domain in 1/s and transform the spatial domain of data into wave number (1/m). In this transformation, seismic data in time (t) versus distance (x) will be transformed into frequency (f) versus wave number (k). The 2D Fourier transform is an extension of the 1D Fourier transformation and expressed by(3)Fω=∫-∞∞∫-∞∞ft,xe-iωte-iktdtdx=FT2Dft,x.


## 2.4. Velocity Analysis

The reflection from a horizontal reflector of a seismic wave on the record for a particular seismic source will follow a parabolic equation(4)tx=t02+xv2.These characteristics are used to estimate the wave propagation speed by assuming the velocity of the medium. Parabolic response patterns can be detected using a shot-gather-based seismic semblance calculation (NE), as noted by Mousa and Al-Shuhail [10].A seismic semblance is formulated as the quotient between stacked and prestacked energy subsequently multiplied by the number of seismic traces (7) for each time window (w). Time window parameters, with width w, are required to perform an NE calculation. Here, N is the number of samples at the time window, w; M is the number of the traces on seismic trace collections; Aij is the amplitude at time i and trace j; Es is stacked energy; and Eu is prestacked energy.(5)Es=∑i=1N∑j=1MAij2,(6)Eu=∑i=1N∑j=1MAij2,(7)NE=EsM×Eu.


## 2.5. Normal Moveout Correction

A normal moveout (NMO) correction is carried out to restore the nonzero offset seismic reflectors to zero offset. The correction is done after we obtaint0 and vNMO for each reflector during velocity analysis process. The time delay after NMO correction is defined in (9). A NMO correction is suitable for a reflector having a slope of less than 15°. When a reflector has higher slope than that threshold, another correction (Dip Move Out correction) is required.(8)ΔtNMOx=tx-t0,(9)tNMOx≈x22t0vNMO2.NMO corrections can be categorized by three classifications: (1) if the NMO velocity (vNMO) is too high, then it will have a shape of a downward convex reflector, (2) if the vNMO is correct, then the reflector will be flat on t=t0, and (3) if vNMO is too low, then the reflector will be an upward convex reflector. During this process, the NMO correction of seismic signals often causes a distortion in the frequency domain. To overcome this problem, a stretching factor (10) is a provision threshold that is required. If a signal is found to have a stretching factor that exceeds the threshold, then a muting correction will be performed.(10)SNMOx=ΔtNMOxt0.


## 2.6. Stacking

Stacking is a process of merging many seismic traces that are already in a state of zero offset. The stacking process is done by quantifying the average amplitude of all signals (11). This stacking process is expected to improve SNR (signal to noise ratio) to factor √M [10].(11)tx=∑i=1MAixM.


## 3. Implementation

SEISGAMA is written using C# programming language provided by Microsoft Visual Studio 2015. Some scientific computing functions and the display are implemented using libraries from Ilnumerics Computing Engine and Ilnumerics Visualization Engine [22]. SEISGAMA has a modular design that aims to facilitate future development of the application. Some modules that are available in the application include the following:GamaseisProcessing andGamaseisView.GamaseisProcessing is on a SEISGAMA module that contains the reading and writing data processing and some physics and mathematical calculations of the seismic reflection method.GamaseisView is a module that controls the graphical user interface, including windows, graphs, seismic sections, labels, and buttons.GamaseisProcessing consists of scripts to read data in a SEG-Y format, performing amplitude correction, muting, velocity analysis, normal moveout correction, and stacking. The SEG-Y seismic data storage format is the most popular to date, and it is used by almost all commercial seismic software. The reading function of the SEG-Y format is already supported in any form; however, until now, this reading function has been limited to two-dimensional seismic reflection data. Other mathematical functions that have been included in this program are implemented in accordance with existing theory. Calculation of Fourier transforms on SEISGAMA uses libraries that have been prepared and optimized by Ilnumerics [22], which implements a special compilation of C language functions to implement the algorithm called the Fastest Fourier Transform from the West (FFTW) [23]. Utilizing a premade library, which has been tested and uses a lower level programming language, is expected to improve the computation time of the Fourier transform. Thus, the calculation of the transformation to F-x and F-K domain is mostly done using FFTW algorithm.GamaseisView consists of a set of views, views that govern the interface between user interaction with the software. Figure1 shows the main user interface of SEISGAMA where users can perform seismic data processing easily.Figure 1
SEISGAMA’s main graphical user interface.


## 4. Results and Discussion

### 4.1. Amplitude Correction

Amplitude correction was successfully implemented on SEISGAMA using a data-dependent exponential method. Users can change the exponential constantsβ, but the program provides a default value of β=16. Figure 2 shows the results of the amplitude correction using the implemented method. The amplitude of the reflectors, which are associated with t0=2400 ms, is not visible (Figure 2(a)), but it became clear after amplitude correction (Figure 2(b)). Amplitude correction not only amplifies the signal but also increases the amplitude of the noise. We can see in Figure 2(b) that refracted and direct waves became increasingly dominant compared to the reflected record after amplitude correction. Furthermore, upon real seismic recording, noncoherent noise will also be amplified after amplitude correction.Figure 2
Before (a) and after (b) amplitude correction using SEISGAMA.
(a)
(b)


### 4.2. Muting

In general, the response of the direct and refracted wave is assertive and intrusive to seismic signals. Therefore, the muting process is necessary to remove responses that degrade or disrupt the signal. In SEISGAMA, the muting process is accomplished by picking the appropriate unwanted points in the muting boundary (Figure3(a)). Then, the user is given a prompt to delete the top of the line which has been picked. The results of the muting process can be seen in Figure 3(b), where waves other than the reflection response have been erased.Figure 3
Muting process on SEISGAMA guided by white line (a) compared with muted shot-gather (b).
(a)
(b)


### 4.3.F-X Dan F-K Transformation

The transformation of seismic data from at-x domain (time and distance) into the F-x and F-K domains can be seen in Figure 4. In the F-x transformation, the vertical axis has been transformed into frequency units. From this display, we determine all of frequencies included in each seismic trace. In the example, the frequency content included on each seismic trace ranged between 100 and 400 Hz. The maximum of the frequency axis is the Nyquist frequency. In this case, the maximum frequency is 1250 Hz; thus the sampling frequency is 2500 Hz. In the F-K transformation, frequency components for each seismic trace are not known directly. Many methods can be performed in the F-K domain to determine specific frequency of the trace; one example is to determine the apparent velocity of the incoming signal [24].Figure 4
Shot-gather inF-x domain (a) and F-K domain (b).
(a)
(b)


### 4.4. Velocity Analysis

In this section, we describe velocity analysis performed on the set of NMO uncorrected traces or on the imperfect NMO corrections. Through analysis of the velocity, we can estimate the velocity of the medium which causes the parabolic phenomenon upon a reflective seismic response. Users are prompted to enter several parameters that will be tested: that is, the width of the gate time, the minimum source wave speed, and the maximum source wave speed. Figure5 shows the result of the velocity analysis applied to the figure on the left using the semblance seismic method. In this case, the user enters 25 ms as the time gate width, 1500 m/s as the minimum speed, and up to 3000 m/s as the maximum speed. The results of the speed analysis are picked based on the maximum values of the observed input for the NMO correction (Figure 6).Figure 5
Velocity analysis result (b) from shot-gather on (a).
(a)
(b)Figure 6
Velocity picking from velocity analysis.


### 4.5. Normal Moveout Correction

NMO correction is performed using the previously chosen parameters used in the velocity analysis. Velocity from the chosen result is associated with a linear relationship; therefore, the velocity is continuous each time. Figure7(b) shows the results of the NMO correction. A correct NMO correction, as shown, is found by flattening seismic reflectors. Figure 8(a) shows that the reflector has become flat at times over 500 ms, while at times below 500 ms the visible convex reflector is still down, which means the velocity is still too high.Figure 7
Uncorrected seismic traces (a) and NMO corrected seismic traces (b).
(a)
(b)Figure 8
NMO corrected seismic traces (a) and stacked seismic trace (b).
(a)
(b)


### 4.6. Stacking

In principle, the stacking process sums up all the seismic traces to improve the signal to noise ratio (high S/N). The stacking process can be carried out for any group of seismic traces, but in the refraction seismic method, this process must be done after the NMO correction, and it should be ensured that the reflector is already flat; therefore the stacking process really can improve the signal over noise. Figure8(b) showed the stacking result from SEISGAMA software.


## 4.1. Amplitude Correction

Amplitude correction was successfully implemented on SEISGAMA using a data-dependent exponential method. Users can change the exponential constantsβ, but the program provides a default value of β=16. Figure 2 shows the results of the amplitude correction using the implemented method. The amplitude of the reflectors, which are associated with t0=2400 ms, is not visible (Figure 2(a)), but it became clear after amplitude correction (Figure 2(b)). Amplitude correction not only amplifies the signal but also increases the amplitude of the noise. We can see in Figure 2(b) that refracted and direct waves became increasingly dominant compared to the reflected record after amplitude correction. Furthermore, upon real seismic recording, noncoherent noise will also be amplified after amplitude correction.Figure 2
Before (a) and after (b) amplitude correction using SEISGAMA.
(a)
(b)


## 4.2. Muting

In general, the response of the direct and refracted wave is assertive and intrusive to seismic signals. Therefore, the muting process is necessary to remove responses that degrade or disrupt the signal. In SEISGAMA, the muting process is accomplished by picking the appropriate unwanted points in the muting boundary (Figure3(a)). Then, the user is given a prompt to delete the top of the line which has been picked. The results of the muting process can be seen in Figure 3(b), where waves other than the reflection response have been erased.Figure 3
Muting process on SEISGAMA guided by white line (a) compared with muted shot-gather (b).
(a)
(b)


## 4.3.F-X Dan F-K Transformation

The transformation of seismic data from at-x domain (time and distance) into the F-x and F-K domains can be seen in Figure 4. In the F-x transformation, the vertical axis has been transformed into frequency units. From this display, we determine all of frequencies included in each seismic trace. In the example, the frequency content included on each seismic trace ranged between 100 and 400 Hz. The maximum of the frequency axis is the Nyquist frequency. In this case, the maximum frequency is 1250 Hz; thus the sampling frequency is 2500 Hz. In the F-K transformation, frequency components for each seismic trace are not known directly. Many methods can be performed in the F-K domain to determine specific frequency of the trace; one example is to determine the apparent velocity of the incoming signal [24].Figure 4
Shot-gather inF-x domain (a) and F-K domain (b).
(a)
(b)


## 4.4. Velocity Analysis

In this section, we describe velocity analysis performed on the set of NMO uncorrected traces or on the imperfect NMO corrections. Through analysis of the velocity, we can estimate the velocity of the medium which causes the parabolic phenomenon upon a reflective seismic response. Users are prompted to enter several parameters that will be tested: that is, the width of the gate time, the minimum source wave speed, and the maximum source wave speed. Figure5 shows the result of the velocity analysis applied to the figure on the left using the semblance seismic method. In this case, the user enters 25 ms as the time gate width, 1500 m/s as the minimum speed, and up to 3000 m/s as the maximum speed. The results of the speed analysis are picked based on the maximum values of the observed input for the NMO correction (Figure 6).Figure 5
Velocity analysis result (b) from shot-gather on (a).
(a)
(b)Figure 6
Velocity picking from velocity analysis.


## 4.5. Normal Moveout Correction

NMO correction is performed using the previously chosen parameters used in the velocity analysis. Velocity from the chosen result is associated with a linear relationship; therefore, the velocity is continuous each time. Figure7(b) shows the results of the NMO correction. A correct NMO correction, as shown, is found by flattening seismic reflectors. Figure 8(a) shows that the reflector has become flat at times over 500 ms, while at times below 500 ms the visible convex reflector is still down, which means the velocity is still too high.Figure 7
Uncorrected seismic traces (a) and NMO corrected seismic traces (b).
(a)
(b)Figure 8
NMO corrected seismic traces (a) and stacked seismic trace (b).
(a)
(b)


## 4.6. Stacking

In principle, the stacking process sums up all the seismic traces to improve the signal to noise ratio (high S/N). The stacking process can be carried out for any group of seismic traces, but in the refraction seismic method, this process must be done after the NMO correction, and it should be ensured that the reflector is already flat; therefore the stacking process really can improve the signal over noise. Figure8(b) showed the stacking result from SEISGAMA software.


## 5. Conclusions

SEISGAMA is a free seismic reflection data processing software program that has successfully processed 2D seismic data for amplitude correction, muting,F-x and F-K domain transform, velocity analysis, NMO correction, and trace stacking. SEISGAMA was built using C# programming language which is natively supported by a Windows operating system. A minimum  .NET 4.5 framework is required to run this software; therefore, a minimum Windows 7 Service Pack 1 is implicitly required. SEISGAMA can give an alternative software choice to users interested in geophysics that has the added benefit of Windows’s compatibility and an intuitive user interface. The basic procedure for the processing of seismic data of reflection has been implemented into SEISGAMA. For advanced seismic processing option, including a Prestack Depth Migration and Full-Waveform inversion, this is currently still under construction and will soon be provided in SEISGAMA. The software is freely available upon request by contacting the authors.



---

*Source: 2913591-2018-01-30.xml*