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

▶
using gpgpu
SMART_READER_LITE
LIVE PREVIEW

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

Jul 11, 2023 160 likes •464 views

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

slide-1
SLIDE 1

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

slide-2
SLIDE 2

Outline

  • Introduction

○ Why is noise attenuation important? ○ Objective ○ Related Works

  • Structure-preserving Smoothing Method
  • Parallel Approach
  • Results
  • Conclusion

2

slide-3
SLIDE 3

Introduction

Why is noise attenuation important?

3

Original seismic data – crossline 905 of F3 block1

1Volume of the Netherlands offshore F3 block downloaded

from the Opendtect website

Zoom area from the input data and the corresponding filtered image.

slide-4
SLIDE 4

Introduction

Why is noise attenuation important?

4

Gaussian filter with radius 1.0 Original image Proposed filter

Preserve relevant structural features is fundamental

slide-5
SLIDE 5

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.

5

Input seismic data – inline 640 of the Kerry data1 Zoom area of the input data and the corresponding semblance attribute Filtered data and the corresponding semblance attribute

1 Kerry data from New Zealand provided for use by New

Zealand Petroleum and Minerals

slide-6
SLIDE 6

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.

6 Original Filtered 3x Original Filtered 3x Zoomed areas of time slice 396ms of F3 Block data

slide-7
SLIDE 7

Introduction

Objective

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

7

slide-8
SLIDE 8

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

slide-9
SLIDE 9

Structure-preserving Smoothing Method

Overview

– Pampanelli P., Seismic amplitude smoothing by anisotropic diffusion preserving structural features, 2015 – Use the instantaneous phase as horizon indicator attribute – Two steps

  • Compute the constraints
  • Assemble and solve the system of

equations

– Iterable

9

Input Volume STEP 1: Computing Seismic Attributes STEP 2: Applying The Anisotropic Diffusion Filter Filtered Volume Do you wish iterate again?

slide-10
SLIDE 10

Structure-preserving Smoothing Method

  • Compute Hilbert Transform

𝑍 𝑢

  • Compute Seismic Complex Trace

𝑎 𝑢 = 𝑌 𝑢 + 𝑗𝑍 𝑢

10

Input Volume STEP 1: Computing Seismic Attributes STEP 2: Applying The Anisotropic Diffusion Filter Filtered Volume Do you wish iterate again?

slide-11
SLIDE 11

Structure-preserving Smoothing Method

  • Compute Horizon Indicator Attribute¹

𝛼Φ = 𝜀Φ 𝜀𝑦 , 𝜀Φ 𝜀𝑧 𝜀Φ 𝜀𝑦 = 1 𝑌2 + 𝑍2 𝑌 𝜀𝑍 𝜀𝑦 − 𝑍 𝜀𝑌 𝜀𝑦 𝜀Φ 𝜀𝑧 = 1 𝑌2 + 𝑍2 𝑌 𝜀𝑍 𝜀𝑧 − 𝑍 𝜀𝑌 𝜀𝑧

where Φ is the instantaneous phase attribute

11

Input Volume STEP 1: Computing Seismic Attributes STEP 2: Applying The Anisotropic Diffusion Filter Filtered Volume

1 – Silva, P.M. et al., Horizon indicator attributes and applications. SEG Technical Program Expanded Abstracts 2012.

Do you wish iterate again?

slide-12
SLIDE 12

Structure-preserving Smoothing Method

  • Compute Horizon Indicator Attribute¹

12

Input Volume STEP 1: Computing Seismic Attributes STEP 2: Applying The Anisotropic Diffusion Filter Filtered Volume Do you wish iterate again? 𝛼Φ = 𝜀Φ 𝜀𝑦 , 𝜀Φ 𝜀𝑧

horizon instantaneous phase gradient 𝛼Φ

𝛼Φ⊥

slide-13
SLIDE 13

Structure-preserving Smoothing Method

  • Compute Fault Attribute¹

𝐺

𝑛𝑏𝑦 = 𝑁𝐵𝑌 𝛼𝑌. 𝛼Φ⊥ , 𝛼𝑍. 𝛼Φ⊥

where 𝛼𝑌 represents the normalized vector of the amplitude gradient and 𝐺

𝑛𝑏𝑦 ∈ −1, 1

13

Input Volume STEP 1: Computing Seismic Attributes STEP 2: Applying The Anisotropic Diffusion Filter Filtered Volume

1 – Pampaneli et al., A new fault-enhancement attribute based on first order directional derivatives of complex trace, EAGE, 2014.

Do you wish iterate again?

slide-14
SLIDE 14

Structure-preserving Smoothing Method

  • Diffusion Tensor

𝐸 = 𝛼Φ⊥ × 𝛼Φ⊥ 𝑈

  • Fault Preserving Factor

𝜁 = 𝐺

𝑛𝑏𝑦 2

14

Input Volume STEP 1: Computing Seismic Attributes STEP 2: Applying The Anisotropic Diffusion Filter Filtered Volume Do you wish iterate again?

slide-15
SLIDE 15

Structure-preserving Smoothing Method

  • The Anisotropic Diffusion Filter

– Solving the linear system Ax = b of size m x m, where m = w x h – Stencil 3x3 𝜀𝑣 𝜀𝑢 = 𝛼 ∙ 𝜁𝐸𝛼𝑣 𝑣𝑦,𝑧

𝑜+1 = 𝑣𝑦,𝑧 𝑜

+ ∆𝑢𝛼 ∙ 𝜁𝑦,𝑧

𝑜 𝐸𝑦,𝑧 𝑜 𝛼𝑣𝑦,𝑧 𝑜+1

15

Input Volume STEP 1: Computing Seismic Attributes STEP 2: Applying The Anisotropic Diffusion Filter Filtered Volume Do you wish iterate again?

slide-16
SLIDE 16

Parallel Approach

Compute Seismic Attributes (STEP 1)

16

CPU

  • FFTW3
  • OpenMP
  • 1 loop for each trace (Hilbert

Transform)

  • 1 loop with 1 inner loop

GPU

  • CUFFT
  • 2 kernels (Hilbert Transform)
  • 1 kernel to compute attributes
  • 1 thread per sample

Memory used in this step

Input seismic section

Hilbert Transform

Y(t) Imaginary part of complex trace Instantaneous phase gradient dxx dxy dyy

Horizon Atribute

Fault Attribute (𝜁) Diffusion tensor

slide-17
SLIDE 17

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)

17

A x = b

CSR NUMERICAL METHOD

Filtered image Input image Input section size 19x Input section size 1x Input section size 1x Input section size 10x + + +

=

Input section size 31x

A b x

̴31x input image size Memory used in this step

slide-18
SLIDE 18

Parallel Approach

Numerical Methods

  • Conjugate gradient (CG)

○ Symmetric linear systems ○ Symmetric and positive-definite matrix Ax = b => AT Ax = AT b

  • Biconjugate gradient stabilized (BiCGStab)
  • Generalized minimal residual (GMRES)

18

Stop criteria Error tolerance: 10-4 Limit of 100 iterations No preconditioner used

slide-19
SLIDE 19

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

slide-20
SLIDE 20

Results

HP Z820 Workstation

  • 64 GB RAM
  • Intel(R) Xeon(R) E5-2620
  • 6 cores (12 threads)

Tesla k40

  • 12 GB
  • 2880 processors

20

Inline 240 of the volume of the Netherlands offshore F3 block¹

  • 951 x 462

Linear system size: 439362 x 439362

1Volume of the Netherlands offshore F3 block

downloaded from the Opendtect website

slide-21
SLIDE 21

Results

21

Original slice Filtered Slice with 3 Iterations

Inline 240 of the volume of the Netherlands offshore F3 block

slide-22
SLIDE 22

Results

22

Original slice 1 iteration 3 iterations 5 iterations 10 iterations

Inline 240 of the volume of the Netherlands offshore F3 block

slide-23
SLIDE 23

Results

23

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

2x 4x 6x 8x 10x 12x 14x 16x 18x 1 2 3 4 5 6 7 8 9 10

Iterations

GPU Speedup using CG

MKL CUSP ViennaCL

slide-24
SLIDE 24

Results

24

Number of iterations BiCGSTAB MKL (ms) CUSP (ms) ViennaCL (ms) 1 358 60 66 2 686 88 116 3 1007 115 166 4

1333

140 215 5

1653

165 265 6

1977

189 313 7

2301

214 363 8

2618

240 413 9

2923

264 461 10

3258

289 511

2x 4x 6x 8x 10x 12x 14x 16x 18x 1 2 3 4 5 6 7 8 9 10

Iterations

GPU Speedup using BiCGStab

MKL CUSP ViennaCL

slide-25
SLIDE 25

Results

25

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

2x 4x 6x 8x 10x 12x 14x 16x 18x 1 2 3 4 5 6 7 8 9 10

Iterations

GPU Speedup using GMRES

MKL CUSP ViennaCL

slide-26
SLIDE 26

Results

26 4891 3258 3467 500 511 525 308 289 220 1000 2000 3000 4000 5000 6000 CG BiCGStab GMRES

Time (ms)

Filtering time with 10 iterations

MKL ViennaCL CUSP 524 358 360 64 66 66 62 60 52 100 200 300 400 500 600 CG BiCGStab GMRES

Time (ms)

Filtering time with 1 iteration

MKL ViennaCL CUSP

slide-27
SLIDE 27

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

slide-28
SLIDE 28

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

slide-29
SLIDE 29

Acknowledgements

Questions???

29