Seismic Imaging: An Overview and Parallel Implementation of Poststack Depth Migration

Seismic Imaging: An Ov ervie w and P arallel Implementation of Poststack Depth Migration Ahmad Shaw ahna ∗ , Syed Abdul Salam † , and Mayez Al-Mouhamed ∗ ∗ Department of Computer Engineering † Department of Electrical Engineering King Fahd Uni versity of Petroleum and Minerals, Dhahran-31261, KSA { g201206920, g201302650, mayez } @kfupm.edu.sa Abstract —Seismic migration is the core step of seismic data processing which is important for oil exploration. Poststack depth migration in frequency-space ( f − x ) domain is one of commonly used algorithms. The wav e-equation solution can be approximated as FIR filtering process to extrapolate the raw data and extract the subsurface image. Because of its computational complexity , its parallel implementation is encouraged. For calcu- lating the next depth lev el, pre vious depth level is r equired. So, this part cannot be parallelized because of data dependence. But at each depth level there is plenty of roam for parallelism and can be parallelized. In case of CUDA programming, each thread calculate a single pixel on the next depth plan. After calculating the next depth plan, we can calculate the depth row by summing over all the frequencies and calculating all the depth ro ws results in the final migrated image. The poststack depth migration is implemented in CUD A and its perf ormance is ev aluated with the sequential code with different problem sizes. Index T erms —Seismic Imaging, Migration, Parallel Comput- ing. f − x extrapolation and Poststack. I . I N T RO D U C T I O N Natural resources like oil, gas, coal and other minerals are important for life. These resources are buried inside earth (both land and marine). T o explore these resources we need a clear image of the earth sub-surface [1]. This image can be obtained using method of reflection seismology . In this method, artificial seismic wav es are generated and these wa ves are reflected from dif ferent geological layers which are recorded by receiv ers (e.g. geophones and hydrophones). The recorded data is in raw form and it needs further processing to obtain the image of earth sub-surface [1], [2]. Digital signal processing has played important role in many applications. Some examples are sonar , radar, medical, communications and seismology [3]. The actual application of signal processing in seismology began with the work of Geophysical Analysis Group at the Massachusetts Institute of T echnology (MIT) in the era 1960-65, it was historical milestones in seismic data processing [4]. T o obtain the image of earth subsurface some preprocessing is required to attenuate noise accompanying the data. By using an imaging technique, the time traces of the preprocessed shot records are transformed into depth traces. For this process a earth model is required and the imaging technique is called migration . After the imaging, a geologist can interpret the migrated section and can identify the layers and structures in the earth subsurface model. For instance, this interpretation can be used to make a decision about the position of a future borehole. If an error is made in the imaging technique and interpretation based on it is wrong, the borehole will miss its target. Therefore, a good quality of imaging technique (or seismic migration) is important. Seimic Migration can be performed in dif ferent domains and in each domain there are number of algorithms. There exists many migration (extrapolation) algorithms. One migration method is frequency-wa venumber ( f − k x ) algorithm [1], [5]. In this method the data in time-space ( t − x ) domain is first transformed into frequency-wa venumber ( f − k x ) domain. At each frequency sample, a complex-v alued FIR filters are ap- plied in the wavenumber -response domain. Another method is frequency-space ( f − x ), in which the data is transformed into frequency-space domain before migration, only the time axis is transformed to the frequency . In this method FIR extrapolation filters (kno wn as extrapolators) are applied in space domain via con volution [5], [6], [7]. By using conv olution, each designed filter output can be calculated independently (or in parallel). In addition, this method can easily be extended to 3-D depth migration. This method can accurately migrate the one-way wa vefields through strong laterally v arying media. In [8], 3D re verse time migration is implemented us- ing GPGPU. Reverse time migration is well-know seismic migration techniche but computationally expensi ve. Seismic migration for both time and depth is also implemented. There is huge performance gain for implementing Kirchhoff prestack time and depth migration on a cluster of 64 GPUs and 256 CPU cores [9]. No wadays, clusters are a v ailable with the GPU; so, its better to utilize all the resources. Some other pre vious work done for seismic signal processing (in different domains) using GPU (CUDA programming) can be found in [10], [11], [12]. In the next section, a brief literature revie w is presented. After that, the implementation of poststack depth migration is discussed and followed by the performance e valuation with the results from sequential code and CUD A code. Finally , the paper is concluded. I I . L I T E R AT UR E R E V I E W The earliest migration technique were graphical and was based on geometrical ideas de veloped. In 1954 Hagedoorn describes the process of seismic migration in terms of prop- agating wa vefronts and tries to av oid the use of non-physical ray paths [13]. Later on Huygens-Fresnel argued that the beam between source and receiver is at least a half wav elength wide, therefore, rather than rays it is better to work with propagating wa vefronts. Huygens’ principle is basis of migration [1], [14]. In seismic migration, wav e propagation effects, can cor - rectly determine the reflection points of the subsurface struc- tures [15]. Migration can be defined as the process of recon- structing a seismic section so that the reflection events are repositioned under their correct surface location at their correct vertical reflection (time or depth) location [1], [16]. The migration process removes the distorting effects of dipping reflectors from the seismic sections. In addition it also remo ves the diffracted arri vals which are resulted from sharp lateral discontinuities [17], [18]. Seismic imaging in the f − x domain is among the attractiv e methods for imaging the earth subsurface structures [1], [5], [6]. This method is implemented via spatial con v olution to obtain seismic images [1], [19], so each output sample can be computed independently , whenever parallel implementation is possible. Most importantly , one-way wa ve extrapolation in the f − x domain can accurately image the subsurface with strong lateral varying medium. For strong varying medium, short length filters are required. Also, the short length filters will reduce the computational cost of con volution. At the same time, to accommodate high propagation angles of wav efields, large length filters become desirable. Hence, the design prob- lem can be treated as an optimization problem between the two trade-of fs [20], [21]. The explicit f − x migration can be easily extended to the 3- D depth migration, which requires 2-D filters as extrapolators. The f − x extrapolation, can accurately image the laterally varying materials by using N length FIR digital filters. In the f − x extrapolation, the seismic wa vefields are spatially sampled, i.e. u ( x i , e j w r , z k ) from depth z k to z k +1 = z k + 4 z and the process can be performed independently for each frequency w r by spatial con volution with a pre-designed filter [22]. u ( x i , e j w r , z k +1 ) = ( N − 1) / 2 X n =( − N +1) / 2 h [ n ] u ( x i − n , e j w r , z k ) (1) In Equation 1, it is sho wn how to apply the designed h [ n ] filter response to extrapolate for each depth level. This process is recursive in nature as shown by the equation 1. Due to this recursiv e nature of algorithm the extrapolation process can be unstable, so accurate designs are required to av oid instability . In the equation 1, x i = i ∆ x and z k = k ∆ z for all i and k ∈ Z (set of integers).. This method can be easily extended to 3-D Seismic Migration. When explicit depth migration is treated as filtering process it resembles con volution and we can compute each sample independently . This method is less expensi ve and can handle lateral v ariations in velocity . I I I . I M P L E M E N TA T I O N In this section we discuss the implementation of Poststack depth migration in f − x domain. In the Algortihm 1, the shot gather is first transformed from t − x to f − x by taking the fourier transform along the time-axis. In order to calculate the image of the subsurface, we need to extrapolate the ra w data based on the wav e-equation. For each depth lev el, we repeat the process the process of extrapolation in f − x domain. It simply a filtering process. This filtering process corresponds to the solution of wave-equation in frequenc y-wav enumber ( f − k x ) domain. From previous depth le vel, ne xt depth level Algorithm 1 Explicit f − x poststack depth migration 1: procedure P O S T S T A C K Input shot gather i in ( t − x ) Fourier transform of shot gather along ( t − axis ) 2: for each depth level do 3: for each freqeuency do 4: for each space do Look-up T able (Predesigned 1-D FIR Filters) Downgoing wav e extrapolation 5: end for 6: end for Calculate depth row by summing over frequencies 7: end for Get the final depth migrated ( z − x ) section 8: end procedure Figure 1. Explicit depth wav efield extrapolation using 1-D FIR filters (courtesy of [23]). is calculated. The depth le vel ( z − x ) is obtained by summing ov er the all frequencies. The cutof f of the filter depends on wa venumber (which itself depends on frequency and velocity) which can vary from point to point. So, the filter cutoffs may vary for each point. Designing filters itself is optimization problem, it takes sometime. In order to speed-up the algorithm, we design the filters in advance and store it on the disk. Based on the filters’ cutoff, the required filter bank is accessed. In Figure 1, accessing the filter and applying it to the previous plan to get a pixel on the next plan is shown. W e repeat the process for all depth lev els/rows and for each frequency to get the final migrated image of the subsurface. In order to implement it on CUD A, two kernels are defined. First kernel is to calculate the next plan f − x from the pre vious using the extrapolation filters. Second kernel is to calculate the depth row from this new plan by summing ov er the frequency axis. Each thread is responsible for calculating the next plan pixel. W e mapped the space-axis to the blocks and frequency- axis to the threads. When the next plan is calculated, we need synchronization of threads for all the pixel. So, when the ne xt plan is calculated by a kernel for next plan, the kernel exits in order to ensure synchronization. Now the second kernel accumulates along the frequency axis and assign it to the depth row . Each depth row represents a ro w in the final migrated z − x section. This z − x section is the final image of the earth subsurface. I V . P E R F O R M A N C E E V A L UAT I ON A N D R E S U LT S In this section, the performance of CUD A code is ev aluated and the results from the sequential and CUD A codes are discussed. T able I sho ws a comparison of the run-time between the CPU and the GPU. W e obtained a speedup of 146X for the CUD A version compared to the CPU code for the largest grid. In Figure 2, the speedup is calculated with the problem size. As the problem size increases the CUDA code performs better . One main reason of lower speedup with smaller problem size is the overhead. So, for smaller problem size the speedup is not that high but as we increase the problem size ( N ), the parallel programs perform better and hence utilizes the GPU in better fashion. A comparison of the run-time for different number of threads and problem sizes are shown in T able II. In Figure 3, we varied the number of threads per block and plotted the speedup against the problem size ( N ). Again, the speedup increases with the problem size and we see that the highest speedup is for 256 threads per block. In this paper , we used benchmark SEG/EAGE salt model and its corresponding zero-of fset ra w data. In Figure 4, the T able I E X EC U T I ON T I ME ( M I NU T E S ) O F D I FF E RE N T M A T R IX S I Z ES . Problem Size CPU GPU 128 2 0.07 0.66 256 2 0.32 0.94 512 2 1.34 0.95 1024 2 12.25 1.14 2048 2 55.87 1.23 4096 2 215.03 1.47 T able II E X EC U T I ON T I ME ( M I NU T E S ) O F D I FF E RE N T M A T R IX S I Z ES A N D N U MB E R O F T H RE A D S . Number of Threads Problem Size 32 64 128 256 256 2 0.83 0.71 0.78 0.99 512 2 2.36 1.56 1.22 1.32 1024 2 3.11 1.99 1.55 1.52 2048 2 4.79 3.10 2.14 1.93 128 256 512 1024 2048 4096 0 30 60 90 120 150 Problem Size (N) Speedup Figure 2. Speed up over sequential program vs Problem size. 256 512 1024 2048 0 5 10 15 20 25 30 Problem Size (N) Speedup 32 64 128 256 Figure 3. Speed up with different number of threads vs Problem size. zero-offset data is sho wn. The zero-of fset data is obtained from the field data after some processing. In Figure 5, SEG/EAGE salt model is shown. In seismic signal processing, salt models are considered complex because of high velocity contrast and imaging beneath the salt is extremely difficult task. In Figures 6 and 7 both the migrated images from sequential Figure 4. SEG/EA GE salt zero-offset section. Figure 5. SEG/EA GE salt velocity model. code and CUD A code are compared. It seems that there is slight error in the final migrated images. W e didn’t achieve 100% correctness. But the depth migrated section is enough close to the original image. Figure 6. Migrated section using sequential code. Figure 7. Migrated section using CUD A code. V . C O N C L U S I O N In conclusion, the performance of the GPU is increasing with the problem size for this particular problem. As there is parallel in the algorithm at the same plan level; so, it is helpful to parallel the algorithm. At depth le vels, the migration process is totally dependent on the pre vious depth plan. Thus, it is important to wait for all the thread to finish and synchronize the parallel program at that le vel. In the future, this work can be easily extended to the prestack depth migration, which is relativ ely more computationally expensi ve. Another e xtension is the block le vel synchronization. W e used exit from the kernel for synchronization, it can also be achiev ed with some other alternativ e. V I . A C K N OW L E D G M E N T The authors would like to thank King Fahd Univ ersity of Petroleum and Minerals (KFUPM) for supporting this research and providing the computing facilities. R E F E R E N C E S [1] O. Y ilmaz, Ed., Seismic Data Analysis: Processing, Inversion, and Interpr etation of Seismic Data , 2nd ed. Society of Exploration Geophysicists, 2001. [2] R. P . Bording and L. R. Lines, Seismic modeling and imaging with the complete wave equation . Society of Exploration Geophysicists Tulsa, 1997. [3] D. E. Dudgeon and R. M. Mersereau, “Multidimensional digital signal processing, ” 1990. [4] B. Buttkus, Spectral Analysis and Filter Theory in Applied Geophysics: W ith 23 T ables . Springer Science & Business Media, 2000. [5] D. Hale, “Stable explicit depth extrapolation of seismic wavefields, ” Geophysics , vol. 56, pp. 1770 – 1777, 1991. [6] O. Holberg, “T owards optimum one-way wav e propagation, ” Geophys- ical Pr ospecting , vol. 36, pp. 99–114, 1988. [7] J. Thorbecke, “Common focus point technology , ” Ph.D. dissertation, Delft Univ ersity of T echnology , 1997. [8] G. Liu, Y . Liu, L. Ren, and X. Meng, “3d seismic reverse time migration on gpgpu, ” Computers & Geosciences , vol. 59, pp. 17–23, 2013. [9] J. Panetta, T . T eixeira, P . R. de Souza Filho, C. A. da Cunha Filho, D. Sotelo, F . M. Roxo da Motta, S. S. Pinheiro, A. L. Romanelli Rosa, L. R. Monnerat, L. T . Carneiro et al. , “ Accelerating time and depth seismic migration by cpu and gpu cooperation, ” International J ournal of P arallel Progr amming , vol. 40, no. 3, pp. 290–312, 2012. [10] H. Knibbe, Reduction of computing time for seismic applications based on the Helmholtz equation by Graphics Processing Units . TU Delft, Delft Univ ersity of T echnology , 2015. [11] S.-Q. W ang, X. Gao, and Z.-X. Y ao, “ Accelerating pocs interpolation of 3d irregular seismic data with graphics processing units, ” Computers & Geosciences , vol. 36, no. 10, pp. 1292–1300, 2010. [12] N. Moussa, “Seismic imaging using gpgpu accelerated re verse time migration, ” CS 315A P arallel Computer Architectur e and Pro gramming , 2009. [13] J. G. Hagedoorn, “ A process of seismic reflection interpretation, ” Geo- physical Pr ospecting , vol. 2, no. 2, pp. 85–127, 1954. [14] J. F . Clœrbout, “Imaging the earth’s interior , ” 1982. [15] E. A. Robinson and S. Treitel, Geophysical signal analysis . Prentice- Hall New Jersey , 1980, vol. 263. [16] P . Keare y , M. Brooks, and I. Hill, An introduction to geophysical explor ation . John Wiley & Sons, 2013. [17] D. E. Dudgeon and R. M. Mersereau, “Multidimensional digital signal processing, ” Prentice-Hall Signal Processing Series, Englewood Cliffs: Pr entice-Hall, 1984 , vol. 1, 1984. [18] L. J. Karam and J. H. McClellan, “Efficient design of digital filters for 2-d and 3-d depth migration, ” Signal Pr ocessing, IEEE T ransactions on , vol. 45, no. 4, pp. 1036–1044, 1997. [19] D. Bhardwaj, S. Y erneni, and S. Phadke, “Parallel computing in seismic data processing, ” in 3rd International P etr oleum Conf . and Exhibition (PETR O TECH-99) . Citeseer , 1999, pp. 279–285. [20] J. W . Thorbecke, K. W apenaar, and G. Swinnen, “Design of one-way wav efield extrapolation operators, using smooth functions in WLSQ optimization, ” Geophysics , vol. 69, no. 4, pp. 1037–1045, 2004. [21] W . A. Mousa, M. V . D. Baan, S. Boussakta, and D. McLernon, “Designing stable operators for explicit depth extrapolation of 2-D & 3-D wavefields using projections onto conv ex sets, ” Geophysics , vol. 74, no. 2, pp. P .S33–S45, March-April 2009. [22] D. Hale, “Stable explicit depth extrapolation of seismic wavefields, ” Geophysics , vol. 56, no. 11, pp. 1770–1777, 1991. [23] W . Mousa, “Imaging of the SEG/EAGE salt model seismic data us- ing sparse–finite-impulse-response wavefield extrapolation filters, ” Geo- science and Remote Sensing, IEEE T ransactions on , vol. 52, no. 5, pp. 2700–2714, 2014.