1. INTRODUCTION

The seismic migration process is one of the most computational expensive steps in seismic data processing. The main purpose of seismic migration is to relocate the collected seismic events in the correct spatial position and collapse the diffraction energy. The Kirchhoff migration is a widely used imaging method, being a flexible and simple method not requiring a high-resolution velocity model. However, when the vertical and horizontal velocity variations are strong, Kirchhoff migration can cause significant imaging issues beyond these variations. More complex imaging methods, such as Reverse Time Migration (RTM) or Wavefield Extrapolation Migration (WEM), have been proposed for imaging complex structures. However, those methods are based on the two-way or one-way wave extrapolation and, therefore, they demand more computational resources than the Kirchhoff migration.

The size of seismic data to be migrated can be hundreds of Terabytes [¹]. This amount of data imposes an l/O bottleneck between the main memory (disk storage system) and the node memory [²],[³].

Some recent computational efforts have been developed to mitigate this I/O bottleneck in seismic applications [³],[⁴].

On the other hand, some works have focused on performing the Kirchhoff migration in a compressed domain to reduce the amount of data to be processed. In [⁵], the samples that represent local extrema in a length window equivalent to the seismic wavelet, are selected and added to a vector that will be migrated. A sorting procedure is used for selecting the samples. In [⁶], the seismic traces are compressed into coefficients by a wavelet decomposition and then, a 2D Kirchhoff migration is applied according to their time location. In [⁷], a 3D Kirchhoff migration on compressed real data is developed using a similar strategy as in [⁶]. In [⁸] a Kirchhoff 3D prestack migration is proposed by compressing the traveltime table, which allows faster access to desirable parts of the traveltime table. These works usually require a post-migration process to improve the quality of the migrated image.

A different approach to perform the Kirchhoff migration on compressed data is through the Matching Pursuit algorithm (MP). Given a dictionary Φ∈R^{N×M}of redundant functions called atoms, it is possible to approximate a signal f∈R^N by a linear combination of k elements from a dictionary, with k << M [⁹]. In [¹⁰], the seismic traces are decomposed by MP using the Ricker wavelets as atoms and then, the migration process is performed atom by atom. In [¹¹], an estimated wavelet is used as atom to avoid noise sensitivity. The time and amplitude values of the selected atoms are migrated as samples of the original data. In [¹²], a multiwavelet-based approach is used to solve the problems regarding stretching and aliasing in prestack Kirchhoff migration. In [¹³], a weighted matching pursuit (WMP) is proposed to perform a high dimensional interpolation, which enables to make a distinction between linear events and spatial aliasing in the FK domain.

In this paper, we compress the seismic data into coefficients obtained by three different Matching Pursuit algorithms: Matching Pursuit (MP), Orthogonal Matching Pursuit (OMP) and Orthogonal Least Squares (OLS). These algorithms have been applied in different seismic applications such as decomposition [¹⁴], [¹⁵] and denoising process [¹⁶]. For all methods, a 10Hz - Ricker wavelet was selected as atom to decompose the seismic data, as it leads to a better approximation of the seismic data than other functional family [¹⁷]. Next, the coefficients are used through a linear combination process to compute the amplitude contributions required by the travel time tables directly in the migrated image. The Ricker waveform equation and some parameters from the compression stage were introduced into the Kirchhoff operator to transform from the Ricker domain to the travel time tables domain during the migration process.

The main advantage of this strategy lies in changing most of the memory access operations by mathematical operations (sums, subtractions, multiplications, etc.). These mathematical operations are related to the compression and decompression algorithms. With this strategy, we seek to reduce the number of transfer operations between the main memory and the local memory. These transfer operations are thousands-fold slower than the processing speed of the parallel co-processors. Therefore, the proposed strategy can reduce the impact of the I/O bottleneck and can improve the performance of seismic applications.

The proposed method was tested using 2D synthetic noiseless data, obtaining considerable compression ratio (up to 20:1) without materially affecting the attributes of the image such as frequency spectrum and amplitude values. For noisy seismic traces, the compression ratio is limited by the Signal-to-Noise Ratio (SNR) of the data and the type of noise. In fact, for real noisy seismic traces, compression can be considered a denoising process [¹⁸].

This paper is organized as follows: Section 2 includes a brief description of the three Matching Pursuit methods applied to the seismic data compression problem. In Section 3, a modified Kirchhoff migration approach is proposed. . Next, in Section 4, we present the results of the proposed method using synthetic data. Lastly, the conclusions are provided in Section 5.

2. FUNDAMENTALS

DECOMPOSING PROCESS BY MATCHING PURSUIT ALGORITHM

MATCHING PURSUIT (MP)

MP is a greedy algorithm that enables the decomposition of a signal into a linear expansion of waveforms, chosen from a redundant dictionary of functions Φ [⁹]. The waveforms are discrete-time signals called atoms, denoted by g_y (t), such that ∥gγ(t)∥2=1,, where Y represents operations of scaling, frequency modulation, and translations. The algorithm selects in each iteration the waveforms that are best correlated to the signal by a succession of orthogonal projections of f(t) on Φ. The decomposition process can be described as follows:

where R ⁱ f is the residual after approximating f on Φ in each step, Φ _Γk represents the columns selected from the dictionary, Γ^k is the set of indices of the selected columns, and α^k is the solution vector. Assuming Rf=f, the algorithm iteratively decomposes the calculated residuals until a stopping criterion. In each step, ∥Rf ∥₂ is minimized by choosing a g ) such that l(Rf,g_Yi >l is maximum.

ORTHOGONAL MATCHING PURSUIT (OMP)

OMP was developed as an improvement of the Matching Pursuit algorithm [¹⁹],[²⁰]. The OMP method is based on the same criterion as MP to select the atoms from the dictionary but, the process to compute the solution vector is different. In each iteration, the selected atoms are used to find a set of coefficients that minimize the distance to f. This process can be described as follows:

where Φ^† _Γk t is the pseudoinverse of Φ_Γk.

ORTHOGONAL LEAST SQUARES (OLS)

The OLS method computes the solution vector as OMP, althoughthe selection of the atom is different than OMP and MP. In each iteration OLS search the atom g_yi that results in the minimum least square error of the residue, taking into account the previously atoms selected. Thus, the selection process can be described as:

where Γ ^k is the set of indices of the selected columns.

As mentioned, the Ricker waveform was selected as atom to decompose the seismic traces, defined as:

where f is the frequency of the signal and t _p is the delay of maximum peak.

KIRCHHOFF MIGRATION ON MP COEFFICIENTS

STANDARD KIRCHHOFF MIGRATION

The Kirchhoff migration is derived from the solution of the wave equation. The solution is performed by applying Green's function [¹²] and it can be described in discrete form by [²¹]:

where M(x,z) is the migrated image in depth and I(r _j, S _k ,tt) represents the samples of the input data required by the travel time tables tt. It should be mentioned that, to perform the Kirchhoff migration it is necessary to compute the traveltime tables first on a given velocity model. The numerical solution of the Eikonal equation by finite differences is one of the methods used to perform travel time [²²]. The Eikonal equation is defined as:

where s(x,z) is the slowness or reciprocal velocity. The value of points of the travel time tables represent the time spent by the wave to get there from a specific source. According to those points, the samples of the input data are mapped in depth to obtain the subsurface image.

3 PROPOSED KIRCHHOFF MIGRATION

In this paper, the relationship between the size of seismic data and the size of the compressed data is defined as compression ratio (CR):

where nt is the number of samples per seismic trace, nx is the number of receivers or seismic traces, and ns is the number of shots.

The proposed method is based on the fact that the coefficients y^k and α^k contain coded Ricker waveforms. The example shown in Figure 1 illustrates this idea.

Figure 1 OMP process on a seismic trace using 1 atom. 

The OMP process was performed on a seismic trace where the stopping criterion was 1 atom. As already mentioned, the output would be a y₀ and its corresponding % These coefficients represent a Ricker waveform located at some point with an amplitude (the red signal). Thus, the proposed method performs mathematical operations on the coefficients during the migration process in order to map the red signal. The proposed method can be divided as follows:

Transform from Ricker to travel time domain: It is necessary to transform the y^k into time domain because of the Kirchhoff operator. However, it is also possible to directly transform into travel time domain instead of time, i.e, just compute the samples required by the travel time tables and not all the reconstructed signal. This can be done by including y^k, α^k into the Ricker equation Equation (4) as follows:

where tt are the samples required by the travel time tables and t _d is the time where the maximum peak occurs in the first atom from the dictionary.

Migrating the atoms: The traditional Kirchhoff method searches and maps samples from time to depth, which are obtained by applying the Equation (8). During the mapping process, it is necessary to bear in mind that the computed samples are part of a combination of Ricker waveforms associated with a respective source and receptor. Accordingly, the Equation (5) is modified as follows

Figure 2 shows an example of the proposed method applied on two coefficients. The red line represents the samples of the signal required by the travel time tables that, in most of the cases, are less than the total number of samples.

Figure 2 Illustration of the proposed method on two coefficients. 

However, most of the computed samples have a value of zero, which did not contribute energy to the image. Therefore, a more fitting strategy is to compute only the portion of the Ricker waveform that provides energy to the image, i.e, the samples located in the lobules of the Ricker. Hence, we reduce the amount of mathematical operations to be performed and save computational efforts during the migration process. Thus, the equation (9) can be rewritten as follows:

where ttY _ip represents the samples required by the travel time tables within the period of each atom.

It is worth noting that, although the development of this strategy requires the reconstruction of some parts of the atoms, the samples of (ttY _ip) are directly computed in the image. Thus, we achieve one of the main objectives of this work, i.e. performing the migration process over compressed data. In this way, the reconstruction (decompress) stage before/after performing the seismic migration is eliminated.

4. RESULTS

In this section, the results of the proposed Kirchoff migration with compressed data for both, synthetic and real velocity models, are discussed.

We used a synthetic model composed by six reflectors (subsurface layers) with a double seismic fault, each one with a different inclination angle (see Figure 3a). The model size is 3500x3000 [m], and the parameters used to model the seismic traces were Δx = Δz = 12.5 [m], Δt = 2 [ms] with a source frequency of 10 [Hz]. In this experiment, the seismic data were compressed with a CR ≈ 76. Figure 3 presents the migrated images obtained by the standard and the proposed Kirchhoff migration using OMP to compress the seismic data.

Figure 3 a) Six-layer Velocity Model. Migrated images obtained by b) the standard migration, and c) the proposed Kirchh of migration. 

Note that the main purpose of the migration process was succesfully achieved, i.e, with a CR of 76, our method correctly locates the reflectors, whichpreserve continuity. Moreover, figure 4 shows the results obtained by the three Matching Pursuit algorithms in terms of SNR and perceptual error in the amplitude of the reflectors for different values of CR. The amplitude error was computed based on the 12-norm.

Figure 4 Relationship between a) SNR vs CR, and b) the percentual error in the amplitude of reflectors vs CR obtained from the 6-layer model. 

In terms of SNR, OLS and OMP showed an exponential decrease behavior, while MP was approximately linear (figure 4a). The results show that OLS obtained the best results in terms of SNR, with a maximum of 41 [dB] for a compression ratio of 10. In terms of % error, OLS obtained the best results as compared with other methods, with a minimum error close to 0.08 % (figure 4b). The three MP algorithms presented a nearly linear behavior.

Figure 5 presents the Fourier spectrum of the migrated images obtained by the standard Kirchhoff migration and our method. Figure 5a depicts the Fourier spectrum of Figure 3b, and figure 5b shows the Fourier spectrum of the migrated image obtained with OMP for a CR of 20. Notice that the shape as well as magnitude are well preserved. The magnitude error is less than 0.1%.

Figure 5 Frequency spectrum of the migrated images obtained by the a) standard Kirchhoff migration and b) proposed Kirchhoff migration using OMP with a CR of 20. 

Second, we also tested the proposed method on a well-known real velocity model such as Hess/SEG 2D, which presents a high complex structure. This model size is 7.2x14.4 [km], with eight layers and a saline dome (figure 6a). A total amount of 81 shot gathers were taken to process and Δx = Δz = 20 [m], Δt = 2 [ms] and 10 [Hz] for the frequency of the source.

Figure 6 a) Hess/SEG velocity model. Migrated images obtained by a) the standard Kirchhoff migration and b) the proposed Kirchhoff migration using OMP 

Second, we also tested the proposed method on a well-known real velocity model such as Hess/SEG 2D, which presents a high complex structure. This model size is 7.2x14.4 [km], with eight layers and a saline dome (figure 6a). A total amount of 81 shot gathers were taken to process and Ax = Az = 20 [m], At = 2 [ms] and 10 [Hz] for the frequency of the source.

with a CR = 58.

Figure 6b presents the migrated image obtained by the standard Kirchhoff migration and figure 6c shows the image obtained by the proposed migration using OMP. It should be noted that, for the complex velocity model, the continuity, number and position of reflectors were preserved, in comparison with the traditional mode. The amplitude and frequency content are also preserved, with magnitude errors below 3 % and, the SNR between the migrated mages was about 14 [dB],

CONCLUSIONS

For the purpose of this work, we introduced a strategy to perform the Kirchhoff migration over compressed data by using three different Matching Pursuit algorithms. The proposed Kirchhoff operator requires fewer memory access operations and more mathematical operations, which could save memory requirements in a real scenario. The decompression stage is eliminated by directly computing the coded information of the coefficients into the migrated image. Synthetic experimental results suggest compressing the seismic data up to 20:1 to preserve the seismic attributes. Nonetheless,, if the interest is to preserve the position and continuity of the subsurface structure, the compression ratio can increase above 50. Future work will focus on the use of our method on real data. We also want to extrapolate our method to a 3D Kirchhoff migration by using a GPU-based cluster. It will also focus on the research of the feasibility to develop a different type of migration (eg. RTM) over compressed seismic data.