reproducible version

book/tccs/rtseis/SConstruct

@@ -0,0 +1,5 @@
+from rsf.tex import *
+import os, sys
+os.environ['PSTEXPENOPTS']='color=y fat=3 fatmult=1.5 serifs=n'
+
+End(lclass='geophysics',options='reproduce',color='ALL',hires='ALL',use='hyperref,listings,moreverb,amsmath,amsfonts')

book/tccs/rtseis/gath/SConstruct

@@ -0,0 +1,168 @@
+from rsf.proj import *
+from rsf.recipes.beg import server as private
+
+import random, math
+
+random.seed(2018)
+
+# Data size
+n2=128
+
+# Fetch data
+Fetch('elfgath.rsf','elf',private)
+
+# Display data
+Flow('gath','elfgath',
+     'dd form=native | put unit1=s unit2=m label2=Offset d2=25 o2=100')
+Result('gath','grey title=')
+
+# Prepare for predictive painting
+Flow('mask','gath','math output=1 | pad n1=1000')
+Flow('gath-pad','gath','pad n1=1000')
+
+Flow('dip','gath-pad mask',
+     'dip rect1=10 rect2=10 p0=0 pmin=0 mask=${SOURCES[1]} order=2')
+Flow('trace','dip','window n2=1 | math output=x1')
+Result('dip',
+       '''
+       window n1=800 |
+       grey  color=j title="Slope" 
+       ''')
+
+# RT with single reference trace
+Flow('rt','dip trace','pwpaint seed=${SOURCES[1]} eps=0.05 order=2')
+
+
+# RT-seislet transform
+Flow('seis', 'gath-pad rt',
+     'rtseislet rt=${SOURCES[1]} eps=0.1 adj=y inv=y unit=y')
+Result('seis',
+       '''
+       put o2=0 d2=1 |
+       window n1=800 |
+       grey  title="RT-seislet Transform" label2=Scale unit2= 
+       ''')
+
+
+# Inverse RT-seislet transform
+Flow('seisinv', 'seis rt',
+     'rtseislet rt=${SOURCES[1]} eps=0.1 inv=y unit=y')
+Result('seisinv', 
+       '''
+       window n1=800 |
+       grey  title="Inverse seislet Transform" 
+       ''')
+
+# Inverse RT-seislet transform using 1% significant coefficients
+Flow('seisrec1','seis rt',
+     '''
+     threshold pclip=1 |
+     rtseislet rt=${SOURCES[1]}
+     eps=0.1 inv=y unit=y 
+     ''')
+Result('seisrec1',
+       '''
+       mutter v0=2600 |
+       window n1=800 |
+       grey  title="Inverse RT-seislet Transform (1%)" 
+       ''')
+
+
+# Denoise using different scales
+max=int(math.log(n2)/math.log(2))
+for m in range(max):
+    scale = int(math.pow(2,m))
+    seis  = 'seis%d' % scale
+    Flow(seis,'seis rt',
+         '''
+         cut f2=%d | 
+         rtseislet rt=${SOURCES[1]} eps=0.1 inv=y unit=y
+         ''' % scale)
+
+Result('seis16',
+       '''
+       mutter v0=2600 |
+       window n1=800 |
+       grey unit1=s unit2=m  title="Denoising result" 
+       ''')
+
+# Denoise
+seis = 'seis%d' % 16 
+rt   = 'rt'
+Result('seisdeno',[seis,rt],
+       '''
+       rtseislet rt=${SOURCES[1]} eps=0.1 adj=y inv=y unit=y |
+       put o2=0 d2=1 |
+       window n1=800 |
+       grey unit1=s title="RT-seislet Transform" label2=Scale unit2=
+       ''')
+
+# Noise section
+Flow('diff',[seis,'gath'],
+     '''
+     mutter v0=2600 |
+     window n1=800 |
+     add scale=1,-1 ${SOURCES[1]} 
+     ''')
+Result('diff','grey unit1=s title="Noise" label2=Scale unit2=')
+
+# Interpolation
+# Resample by 2
+Flow('rand','seis',
+     '''
+     window n1=800 | 
+     noise rep=y var=20000 seed=2006 | 
+     mutter v0=2600 | pad n1=1000
+     ''')
+
+Flow('seisrand','seis rand','cat axis=2 ${SOURCES[1]} | put d2=12.5')
+
+Flow('rt2','rt',
+     'transp | remap1 d1=12.5 o1=100 n1=256 | transp')
+
+Flow('gath2','seisrand rt2',
+     'rtseislet rt=${SOURCES[1]} eps=0.1 inv=y unit=y')
+
+# Resample by 4
+Flow('rand2','seisrand',
+     '''
+     window n1=800 | 
+     noise rep=y var=20000 seed=2006 | 
+     mutter v0=2600 | pad n1=1000
+     ''')
+
+Flow('seisrand2','seisrand rand2','cat axis=2 ${SOURCES[1]} | put d2=6.25')
+Result('seisrand2',
+    '''
+    put o2=0 d2=1 | 
+    window n1=800 | 
+    grey unit1=s title="Seislet Transform" label2=Scale unit2=
+    ''')
+
+Flow('rt4','rt2',
+     'transp | remap1 d1=6.25 o1=100 n1=512 | transp')
+
+Flow('gath4','seisrand2 rt4',
+     'rtseislet rt=${SOURCES[1]} eps=0.01 inv=y unit=y')
+Result('gath4',
+       '''
+       mutter v0=2600 |
+       window n1=800 |
+       grey unit1=s unit2=m title="Interpolated shot gather"
+       ''')
+
+# FFT
+Flow('fft','gath','addtrace ratio=4 | fft1 | fft3')
+Result('fft',
+    '''
+    math output="abs(input)" | real | 
+    grey wanttitle=
+    ''')
+Flow('fft2','gath4','window n1=800 | fft1 | fft3')
+Result('fft2',
+    '''
+    math output="abs(input)" | real |
+    grey wanttitle=
+    ''')
+
+End()

book/tccs/rtseis/noise/SConstruct

@@ -0,0 +1,127 @@
+from rsf.proj import *
+
+# Size of the model
+n1=200
+n2=256
+
+# Generate and display the model
+Flow('sigmoid',None,
+        '''
+        sigmoid n1=%d n2=%d d2=.008 |
+        smooth rect1=3 diff1=1 adj=1 | smooth rect1=3 |
+        put label2=Distance
+        ''' % (n1,n2) )
+Result('sigmoid', 'grey title=Input')
+
+
+# Sensitivity to noise
+noise_var = [0.05, 0.1, 0.2, 0.4, 0.8, 1]
+snr_rt = []
+snr_pwd = []
+snr_noise = []
+
+# Normalize the model
+Flow('sigmoid-norm','sigmoid','scale axis=12 | scale dscale=6')
+
+for i in range(len(noise_var)):
+    var = noise_var[i]
+    sigmoid     = 'sigmoid-noise%d' % i
+    pad         = 'pad-noise%d' % i
+    sigmoidpad  = 'sigmoidpad-noise%d' % i
+    dip         = 'dip-noise%d' % i
+    dippad      = 'dippad-noise%d' % i
+    seed        = 'seed-noise%d' % i
+    rt          = 'rt-noise%d' % i
+    rtseis      = 'rtseis-noise%d' % i
+    rtseisrec   = 'rtseisrec-noise%d' % i
+    pwdseis     = 'pwdseis-noise%d' % i
+    pwdseisrec  = 'pwdseisrec-noise%d' % i
+    diffrt      = 'diff-rt-noise%d' % i
+    snrrt       = 'snr-rt%d' % i
+    diffpwd     = 'diff-pwd-noise%d' % i
+    snrpwd      = 'snr-pwd%d' %  i
+    n           = 'noise%d' % i
+    snrnoise    = 'snr-noise%d' % i
+
+    snr_rt.append(snrrt)
+    snr_pwd.append(snrpwd)
+    snr_noise.append(snrnoise)
+
+    # Add noise
+    Flow(sigmoid,'sigmoid-norm','noise seed=2019 var=%g' % var)
+    Flow(n,[sigmoid,'sigmoid'],'add ${SOURCES[1]} scale=-1,1')
+    Flow(snrnoise,None,"math n1=1  output='%g' " % var)
+
+    # Prepare for RT-seislet transform
+    Flow(pad,sigmoid,'math output=1 | pad beg1=50 end1=50')
+    Flow(sigmoidpad,sigmoid,'pad beg1=50 end1=50 | bandpass fhi=60')
+    Flow(dippad,[sigmoidpad,pad],
+         '''
+         dip order=2 p0=0 verb=y niter=10 rect1=3 rect2=3
+         mask=${SOURCES[1]} 
+         ''')
+    Flow(seed,dippad,'window n2=1 | math output=x1')
+    Flow(dippad,[sigmoidpad,pad],
+         '''
+         dip order=2 p0=0 verb=y niter=10 rect1=10 rect2=10
+         mask=${SOURCES[1]}
+         ''')
+    # Compute RT
+    picks=[]
+    for i0 in (5,10,80,128,156,175,250,251,255):
+        pick = 'pick-noise%d-%d' % (i,i0)
+        picks.append(pick)
+
+        # RT with single reference trace
+        Flow(pick,[dippad,seed],
+             'pwpaint order=2 seed=${SOURCES[1]} i0=%d eps=1' % i0)
+    np = len(picks)
+    Flow(rt,picks,
+         'add ${SOURCES[1:%d]} | scale dscale=%g' % (np,1.0/np))
+
+    # RT-seislet transform
+    Flow(rtseis,[sigmoidpad,rt],
+         'rtseislet rt=${SOURCES[1]} eps=0.1 adj=y inv=y unit=y type=b')
+    Flow(rtseisrec,[rtseis,rt],
+         '''
+         threshold pclip=10 |
+         rtseislet rt=${SOURCES[1]} eps=0.1 inv=y unit=y type=b | 
+         window n1=200 min1=0
+         ''')
+    # Compute SNR
+    Flow(diffrt,[rtseisrec,'sigmoid-norm'],'add ${SOURCES[1]} scale=-1,1')
+    Flow(snrrt,[rtseisrec,diffrt],'snr2 noise=${SOURCES[1]}')
+
+    # Estimate dips for PWD-seislet tranform
+    Flow(dip,sigmoid,'dip rect1=10 rect2=10 p0=0 pmin=-100')
+
+    # PWD-seislet transform
+    Flow(pwdseis,[sigmoid,dip],
+         'seislet dip=${SOURCES[1]} eps=0.1 adj=y inv=y unit=y type=b')
+    Flow(pwdseisrec,[pwdseis,dip],
+         '''
+         threshold pclip=10 |
+         seislet dip=${SOURCES[1]} eps=0.1 inv=y unit=y type=b
+         ''')
+    # Compute SNR
+    Flow(diffpwd,[pwdseisrec,'sigmoid-norm'],'add ${SOURCES[1]} scale=-1,1')
+    Flow(snrpwd,[pwdseisrec,diffpwd],'snr2 noise=${SOURCES[1]}')
+
+Flow('snr-rt',snr_rt,'cat axis=1 ${SOURCES[1:%d]}' % len(snr_rt))
+Flow('snr-pwd',snr_pwd,'cat axis=1 ${SOURCES[1:%d]}' % len(snr_pwd))
+Flow('snr-noise',snr_noise,'cat axis=1 ${SOURCES[1:%d]}' % len(snr_noise))
+Flow('snr-rt-noise','snr-noise snr-rt','cmplx ${SOURCES[1]}')
+Flow('snr-pwd-noise','snr-noise snr-pwd','cmplx ${SOURCES[1]}')
+Flow('snr','snr-rt-noise snr-pwd-noise','cat axis=2 ${SOURCES[1]}')
+
+Result('snr',
+       '''
+       graph wanttitle=n dash=0,1 screenratio=1. labelsz=9 plotfat=9
+       label2="S/N" unit2=dB label1="Noise variances" unit1=
+       ''')
+
+# Example with noise
+Result('sigmoid-noise0','grey title="Data with noise"')
+Result('rtseisrec-noise0','grey title="Reconstructed data"')
+
+End()