using GPGPU Joner Duarte jduartejr@tecgraf.puc-rio.br Outline - PowerPoint PPT Presentation

GPU Technology Conference 2016 – April, 4-7 – San Jose, CA, USA Structure-preserving Smoothing for Seismic Amplitude Data by Anisotropic Diffusion using GPGPU Joner Duarte jduartejr@tecgraf.puc-rio.br

Outline ● Introduction ○ Why is noise attenuation important? ○ Objective ○ Related Works ● Structure-preserving Smoothing Method ● Parallel Approach ● Results ● Conclusion 2

Introduction Why is noise attenuation important? Original seismic data – crossline Zoom area from the input data and 905 of F3 block 1 the corresponding filtered image. 1 Volume of the Netherlands offshore F3 block downloaded 3 from the Opendtect website

Introduction Why is noise attenuation important? Preserve relevant structural features is fundamental Original image Gaussian filter Proposed filter with radius 1.0 4

Introduction Why is noise attenuation important? Seismic attributes have been used to accelerate interpretations, and their results are directly related to the quality of the seismic data. Input seismic data – inline 640 of the Zoom area of the input data and Filtered data and the Kerry data 1 the corresponding semblance corresponding semblance attribute attribute 1 Kerry data from New Zealand provided for use by New 5 Zealand Petroleum and Minerals

Introduction Why is noise attenuation important? Curvature attributes are widely used to improve seismic interpretation process. The results of curvature attribute can be directly improved by our iterative noise attenuation process. Filtered 3x Original Original Filtered 3x Zoomed areas of time slice 396ms of F3 Block data 6

Introduction Objective Remove noise preserving relevant details such as structural and stratigraphic discontinuities in interactive time. 7

Related Works ● [1990] Perona and Malik, Scale-Space and Edge Detection Using Anisotropic Diffusion ● [2002] Hocker and Fehmers, Fast structural interpretation with structure-oriented filtering ● [2009] Dave Hale, Structure-oriented smoothing and semblance ● [2011] Baddari et al., Seismic Noise Attenuation by Means of an Anisotropic Non-linear Diffusion Filter 8

Structure-preserving Smoothing Method Overview Input Volume – Pampanelli P., Seismic amplitude STEP 1: Computing Seismic smoothing by anisotropic diffusion Attributes preserving structural features , 2015 – Use the instantaneous phase as STEP 2: Applying The horizon indicator attribute Anisotropic Diffusion Filter – Two steps • Compute the constraints Filtered • Assemble and solve the system of Volume Do you wish equations iterate again? – Iterable 9

Structure-preserving Smoothing Method • Compute Hilbert Transform Input Volume 𝑍 𝑢 STEP 1: Computing Seismic Attributes • Compute Seismic Complex Trace STEP 2: Applying The 𝑎 𝑢 = 𝑌 𝑢 + 𝑗𝑍 𝑢 Anisotropic Diffusion Filter Filtered Volume Do you wish iterate again? 10

Structure-preserving Smoothing Method • Compute Horizon Indicator Attribute¹ Input 𝛼Φ = 𝜀Φ 𝜀𝑦 , 𝜀Φ Volume 𝜀𝑧 STEP 1: Computing Seismic Attributes 𝜀Φ 𝑌 2 + 𝑍 2 𝑌 𝜀𝑍 1 𝜀𝑦 − 𝑍 𝜀𝑌 𝜀𝑦 = 𝜀𝑦 STEP 2: Applying The Anisotropic Diffusion Filter 𝜀Φ 𝑌 2 + 𝑍 2 𝑌 𝜀𝑍 1 𝜀𝑧 − 𝑍 𝜀𝑌 𝜀𝑧 = 𝜀𝑧 Filtered where Φ is the instantaneous phase Volume Do you wish attribute iterate again? 1 – Silva, P.M. et al., Horizon indicator attributes and applications. SEG Technical Program Expanded Abstracts 2012. 11

Structure-preserving Smoothing Method • Compute Horizon Indicator Attribute¹ Input 𝛼Φ = 𝜀Φ 𝜀𝑦 , 𝜀Φ Volume 𝜀𝑧 STEP 1: Computing Seismic Attributes instantaneous phase gradient 𝛼Φ STEP 2: Applying The 𝛼Φ ⊥ Anisotropic Diffusion Filter horizon Filtered Volume Do you wish iterate again? 12

Structure-preserving Smoothing Method • Compute Fault Attribute¹ Input Volume 𝑛𝑏𝑦 = 𝑁𝐵𝑌 𝛼𝑌. 𝛼Φ ⊥ , 𝛼𝑍. 𝛼Φ ⊥ 𝐺 STEP 1: Computing Seismic Attributes where 𝛼𝑌 represents the normalized STEP 2: Applying The vector of the amplitude gradient and Anisotropic Diffusion Filter 𝐺 𝑛𝑏𝑦 ∈ −1, 1 Filtered Volume Do you wish iterate again? 1 – Pampaneli et al., A new fault-enhancement attribute based on first order directional derivatives of complex trace, EAGE, 2014. 13

Structure-preserving Smoothing Method • Diffusion Tensor Input 𝐸 = 𝛼Φ ⊥ × 𝛼Φ ⊥ 𝑈 Volume STEP 1: Computing Seismic Attributes • Fault Preserving Factor 2 STEP 2: Applying The 𝜁 = 𝐺 𝑛𝑏𝑦 Anisotropic Diffusion Filter Filtered Volume Do you wish iterate again? 14

Structure-preserving Smoothing Method • The Anisotropic Diffusion Filter Input Volume – Solving the linear system Ax = b of size STEP 1: Computing Seismic m x m, where m = w x h Attributes – Stencil 3x3 STEP 2: Applying The 𝜀𝑣 Anisotropic Diffusion Filter 𝜀𝑢 = 𝛼 ∙ 𝜁𝐸𝛼𝑣 Filtered 𝑜+1 = 𝑣 𝑦,𝑧 𝑜 𝛼𝑣 𝑦,𝑧 𝑜 𝐸 𝑦,𝑧 𝑜+1 𝑜 𝑣 𝑦,𝑧 + ∆𝑢𝛼 ∙ 𝜁 𝑦,𝑧 Volume Do you wish iterate again? 15

Parallel Approach Compute Seismic Attributes ( STEP 1) CPU GPU • FFTW3 • CUFFT • OpenMP • 2 kernels (Hilbert Transform) • 1 loop for each trace (Hilbert • 1 kernel to compute attributes • 1 thread per sample Transform) • 1 loop with 1 inner loop Memory used in this step Fault Attribute ( 𝜁 ) Hilbert Horizon Y(t) Instantaneous Input seismic Transform Atribute Imaginary part of phase gradient section complex trace dxx dxy dyy Diffusion tensor 16

Parallel Approach Assemble and Solve the System of Equations (STEP 2) – System solver • 50~70% of total time of one iteration • Libraries for solving sparse linear systems • Sparse matrix format: CSR (Compressed Sparse Row) Memory used in this step ̴ 31x input image size A x b A x = b Input Input Input Input Input = 19x + 1x + 1x + 10x 31x section size section size section size section size section size Filtered image Input image CSR NUMERICAL METHOD 17

Parallel Approach Numerical Methods Stop criteria ● Conjugate gradient ( CG ) Error tolerance: 10 -4 Limit of 100 iterations ○ Symmetric linear systems ○ Symmetric and positive-definite matrix No preconditioner used Ax = b => A T Ax = A T b ● Biconjugate gradient stabilized ( BiCGStab ) ● Generalized minimal residual ( GMRES ) 18

Parallel Approach Libraries ● Intel MKL 11.3 ● CUSP 0.5.1 ● ViennaCL 1.7.0 ● There are other libraries that we may try in the future (MUMPS, cuSOLVER, PARALUTION, etc) 19

Results HP Z820 Workstation Inline 240 of the volume of the Netherlands offshore ● 64 GB RAM F3 block¹ ● Intel(R) Xeon(R) E5-2620 ● 951 x 462 ● 6 cores (12 threads) Tesla k40 Linear system size: 439362 x 439362 ● 12 GB 1 Volume of the Netherlands offshore F3 block ● 2880 processors downloaded from the Opendtect website 20

Results Inline 240 of the volume of the Netherlands offshore F3 block Filtered Slice with 3 Iterations Original slice 21

Results Inline 240 of the volume of the Netherlands offshore F3 block Original slice 1 iteration 3 iterations 5 iterations 10 iterations 22

Results CG Number of MKL CUSP ViennaCL GPU Speedup using CG iterations (ms) (ms) (ms) 18x 1 524 62 64 16x 14x 2 1009 92 113 12x 3 1549 119 162 10x 4 2028 147 210 8x 6x 5 2523 174 259 4x 6 3013 201 308 2x 7 3510 229 356 1 2 3 4 5 6 7 8 9 10 8 3976 255 405 Iterations 9 4429 282 451 MKL CUSP ViennaCL 10 4891 308 500 23

Results BiCGSTAB Number of MKL CUSP ViennaCL GPU Speedup using BiCGStab iterations (ms) (ms) (ms) 18x 1 358 60 66 16x 14x 2 686 88 116 12x 3 1007 115 166 10x 4 140 215 8x 1333 6x 5 165 265 1653 4x 6 189 313 1977 2x 7 214 363 2301 1 2 3 4 5 6 7 8 9 10 8 240 413 2618 Iterations 9 264 461 2923 MKL CUSP ViennaCL 10 289 511 3258 24

Results GMRES Number of MKL CUSP ViennaCL GPU Speedup using GMRES iterations (ms) (ms) (ms) 18x 1 360 52 66 16x 14x 2 706 71 117 12x 3 1051 90 168 10x 4 1395 109 220 8x 6x 5 1740 127 270 4x 6 2090 146 322 2x 7 2430 164 372 1 2 3 4 5 6 7 8 9 10 8 2781 183 424 Iterations 9 3125 201 474 MKL CUSP ViennaCL 10 3467 220 525 25

Results Filtering time with 10 Filtering time with 1 iterations iteration 6000 600 524 4891 5000 500 Time (ms) 4000 Time (ms) 3467 358 360 400 3258 3000 300 2000 200 66 66 64 62 60 100 52 1000 525 500 511 308 289 220 0 0 CG BiCGStab GMRES CG BiCGStab GMRES MKL ViennaCL CUSP MKL ViennaCL CUSP 26

Conclusions - Method that attenuates noise preserving structural features - Improves interpretation eficiency and also can enhance the results of other seismic attributes - Fast enough to be interactive 27

Conclusions - It can be used on demand during the interpretation process and also be integrated with tools like Cintiq tablet - CUSP library proved to be very fast for our problem with a high level interface that makes it very simple to use 28

Acknowledgements Questions??? 29