```from rsf.proj import * import math from rsf.recipes.tpx import FPX ### plotting functions def section(title,label1='Time',unit1='s',min1=5.8,max1=8.0,extra=" "): return ''' window min1=%g max1=%g | grey title="%s" label1="%s" unit1="%s" label2=Distance unit2=m %s ''' % (min1,max1,title,label1,unit1,extra) def sectionw(title,label1='Time',unit1='s',min1=5.5,max1=7.0): return ''' window min1=%g max1=%g max2=5500 | grey title="%s" label1="%s" unit1="%s" label2=Distance unit2=m ''' % (min1,max1,title,label1,unit1) def plotvel(title,min1=5.5,max1=8.0): return ''' grey min1=%g max1=%g minval=1480 maxval=2000 bias=1700 clip=200 allpos=n color=j scalebar=y barreverse=y barunit=m/s barlabel="Velocity" title="%s" label2="Distance" unit2=m ''' %(min1,max1,title) def plotdip(title,min1=5.5,max1=8.0): return ''' grey min1=%g max1=%g color=j scalebar=y barlabel="Dip" title="%s" label2="Distance" unit2=m ''' %(min1,max1,title) def plotdipw(title,label1='Time',unit1='s',min1=5.5,max1=7.0): return ''' window min1=%g max1=%g max2=5500 | grey title="%s" color=j scalebar=y barlabel="Dip" label1="%s" unit1="%s" label2=Distance unit2=m ''' % (min1,max1,title,label1,unit1) frame1=(int)( (7.0-5.5)/0.002 ) frame2=(int)( (2000)/16.667 ) frame3=(int)( (1500-1400)/20 ) print frame1 print frame2 print frame3 def plotdip3d(title,min1=5.5,max1=8.0,fr1=frame1,fr2=frame2,fr3=frame3): return ''' window min1=%g max1=%g max2=5500 | byte | grey3 flat=n color=j scalebar=n barlabel="Dip" title="%s" label2="Distance" label3="Velocity" unit3="m/s" unit2=m frame1=%d frame2=%d frame3=%d point1=.6 point2=.7 ''' %(min1,max1,title,fr1,fr2,fr3) def plotpick(min2=1400,max2=2600): return ''' graph yreverse=y transp=y plotcol=7 plotfat=7 pad=n min2=%g max2=%g wantaxis=n wanttitle=n '''%(min2,max2) def plot3d(title,min1=5.5,max1=8.0,fr1=frame1,fr2=frame2,fr3=frame3): return ''' window min1=%g max1=%g max2=5500 | byte | grey3 flat=n scalebar=n barlabel="Dip" title="%s" label2="Distance" label3="Velocity" unit3="m/s" unit2=m frame1=%d frame2=%d frame3=%d point1=.6 point2=.7 ''' %(min1,max1,title,fr1,fr2,fr3) ####### Download data Fetch('Nshots.su','nankai') ####### Convert to RSF Flow('shots tshots','Nshots.su', 'suread suxdr=y tfile=\${TARGETS[1]}') ####### DC Removal Flow('mean','shots','stack axis=1 | spray axis=1 n=5500 o=0.0 d1=0.002') Flow('shotsdc','shots mean','add scale=1,-1 \${SOURCES[1]}') ####### Bandpass Filtering Flow('shotsf','shotsdc','bandpass flo=10') ####### Windowing Data Flow('shotsw','shotsf','window min1=5.5') # windowing 5.5 and bandpass filtering significantly # improves surface-consistent analysis ####### Mask zero traces Flow('mask0','shotsw','mul \$SOURCE | stack axis=1 | mask min=1e-20') Flow('shots0','shotsw mask0','headerwindow mask=\${SOURCES[1]}') # update a database Flow('tshots0','tshots mask0','headerwindow mask=\${SOURCES[1]}') ####### Surface consistent # Average trace amplitude Flow('arms','shots0', 'mul \$SOURCE | stack axis=1 | math output="log(input)"') # shot/offset indeces: fldr and tracf Flow('indexshot','tshots0','window n1=1 f1=2') Flow('offsets4index','tshots0',' headermath output=offset | dd type=float | window') Flow('offsetindex','offsets4index','math output="abs(input) - 170" | dd type=int') # receiver/midpoint Flow('midpoint','tshots0','window n1=1 f1=5') Flow('cmps4index','tshots0',' headermath output=cdp | dd type=float | math output="input*16.667" | window') Flow('recv','cmps4index offsets4index','add scale=1,0.5 \${SOURCES[1]} | math output="input - 13799" | dd type=int') Flow('index','indexshot offsetindex', ''' cat axis=2 \${SOURCES[1]} ''') Flow('extindex','index midpoint', ''' cat axis=2 \${SOURCES[1]} ''') Flow('extindrecv','extindex recv', ''' cat axis=2 \${SOURCES[1]} ''') def plot(title): return ''' spray axis=1 n=1 | intbin head=\${SOURCES[1]} yk=fldr xk=tracf | window | grey title="%s" label2="Shot Number" unit2= label1="Offset Number" unit1= scalebar=y ''' % (title) def plotb(title,bias=-5): return ''' spray axis=1 n=1 | intbin head=\${SOURCES[1]} yk=fldr xk=tracf | window | grey title="%s" label2="Shot Number" unit2= label1="Offset Number" unit1= scalebar=y clip=3 bias=%g ''' % (title,bias) # Display in shot/offset coordinates Flow('varms','arms tshots0','spray axis=1 n=1 | intbin head=\${SOURCES[1]} yk=fldr xk=tracf | window') Plot('varms','arms tshots0',plotb('Log-Amplitude')) prog = Program('surface-consistent.c') sc = str(prog[0]) # recv index # get model dimensions Flow('recvmodel',['arms','extindrecv',sc], './\${SOURCES[2]} index=\${SOURCES[1]} verb=y') # find a term Flow('recvsc',['arms','extindrecv',sc,'recvmodel'], ''' conjgrad ./\${SOURCES[2]} index=\${SOURCES[1]} mod=\${SOURCES[3]} niter=150 ''') # project to a data space Flow('recvscarms',['recvsc','extindrecv',sc], './\${SOURCES[2]} index=\${SOURCES[1]} adj=n') Plot('recvvscarms','recvscarms tshots0', plotb('Source, Offset, CDP, Recv S-C Log(A)')) # compute difference Flow('recvadiff','arms recvscarms','add scale=1,-1 \${SOURCES[1]}') Plot('recvadiff','recvadiff tshots0',plot('s,h,cdp,r difference')) ### apply to traces to all times - no windowing is considered Flow('ampl','recvscarms', 'math output="exp(-input/2)" | spray axis=1 n=5500 d=0.002 o=0') Flow('shotsf0','shotsf mask0','headerwindow mask=\${SOURCES[1]}') Flow('shots-preproc','shotsf0 ampl','mul \${SOURCES[1]}') Plot('shots-preproc','shots-preproc','window n2=100 | grey min1=6.0 max1=8.0 title="Shots Preproc"') Plot('shots-raw','shots0','window n2=100 | grey min1=6.0 max1=8.0 title="Shots Raw"') ### make stack figures to compare # Resample to 4 ms Flow('subsampled','shots-preproc', ''' bandpass fhi=125 | window j1=2 ''') #Result('spectra','subsampled', # 'spectra all=y | graph title="Subsampled Spectra"') #Result('spectra-check','shots-preproc', # 'spectra all=y | graph max1=160 title="Spectra Check"') Flow('cmps mask','subsampled tshots0', ''' intbin xk=tracf yk=cdp head=\${SOURCES[1]} mask=\${TARGETS[1]} | put o3=900 d3=16.667 label3=Distance unit3=m | pow pow1=2 ''') Flow('cmps-raw mask-raw','shots0 tshots0', ''' intbin xk=tracf yk=cdp head=\${SOURCES[1]} mask=\${TARGETS[1]} | put o3=900 d3=16.667 label3=Distance unit3=m ''') Flow('offset-file','tshots0', ''' sfheadermath output=offset | dd type=float | intbin xk=tracf yk=cdp head=\$SOURCE ''') Flow('watervel','cmps','window n2=1 | math output="1500"') Flow('nmo-wat-vel-preproc','cmps offset-file mask watervel', ''' nmo half=n offset=\${SOURCES[1]} mask=\${SOURCES[2]} velocity=\${SOURCES[3]} ''') Flow('nmo-wat-vel-raw','cmps-raw offset-file mask watervel', ''' nmo half=n offset=\${SOURCES[1]} mask=\${SOURCES[2]} velocity=\${SOURCES[3]} ''') Flow('stack-wat-vel-preproc','nmo-wat-vel-preproc','stack') Flow('stack-wat-vel-raw','nmo-wat-vel-raw','stack') #Result('stack-preproc','stack-wat-vel-preproc',section('Stack 1.5 km/s Preproc')) #Result('stack-raw','stack-wat-vel-raw',section('Stack 1.5 km/s Raw')) ### DMO stacking nv=60 Flow('stacks','cmps offset-file mask', ''' window min1=4 | stacks half=n v0=1400 nv=%g dv=20 offset=\${SOURCES[1]} mask=\${SOURCES[2]} '''%nv,split=[3,'omp']) #>>> # Consider stacking without time windowing to check # if slope decomposition is affected #Flow('nw-stacks','cmps offset-file mask', # ''' # stacks half=n v0=1400 nv=%g dv=20 # offset=\${SOURCES[1]} mask=\${SOURCES[2]} # '''%nv,split=[3,'omp']) ### interface p = 100 min1=5.5 max1=8.0 mute = ''' mutter v0=130 x0=1300 t0=4.0 half=n inner=n | mutter x0=1400 v0=20 t0=5.0 half=n inner=y | mutter v0=2500 x0=1400 t0=5.8 half=n inner=n | mutter v0=500 x0=1400 t0=7.0 half=n inner=y ''' pick = ''' pick rect1=50 rect2=20 vel0=1480 ''' mutec = ''' mutter v0=130 x0=1300 t0=4.0 half=n inner=n | mutter x0=1400 v0=20 t0=5.0 half=n inner=y | mutter v0=2500 x0=1400 t0=5.8 half=n inner=n | mutter v0=500 x0=1400 t0=7.0 half=n inner=y ''' pickc = ''' pick rect1=80 rect2=20 vel0=1400 ''' # Compute envelope for picking (stack without DMO) Flow('envelope-stacks','stacks','envelope | scale axis=2',split=[3,'omp']) Flow('vpick-stacks','envelope-stacks',mutec + ' | ' + pickc) # Take a slice Flow('slice-stacks','stacks vpick-stacks','slice pick=\${SOURCES[1]}') #Result('vpick-stacks',plotvel('NMO velocity')) #Result('stacks-check-w','stacks','window n2=1 min2=1499 | window max2=5500 | grey title="Stack 1500 m/s (extracted)" min1=5.5 max1=7.0') #Result('stacks-check','stacks','window n2=1 min2=1499 | grey title="Stack 1500 m/s (extracted)" min1=4.0') Flow('stackst','stacks','costaper nw3=20') # Apply double Fourier transform (cosine transform) # pad n3=601 | Flow('cosft','stackst','cosft sign1=1 sign3=1') # Transpose f-v-k to v-f-k Flow('transp','cosft','transp',split=[3,'omp']) # Fowler DMO: mapping velocities Flow('map','transp', ''' math output="x1/sqrt(1+0.25*x3*x3*x1*x1/(x2*x2))" | cut n2=1 ''') Flow('fowler','transp map','iwarp warp=\${SOURCES[1]} | transp', split=[3,'omp']) # Inverse Fourier transform # | window n3=401 # does not improve the result severely Flow('dmo','fowler','cosft sign1=-1 sign3=-1') #>>> #Flow('nw-stackst','nw-stacks','costaper nw3=20') #Flow('nw-cosft','nw-stackst','cosft sign1=1 sign3=1') #Flow('nw-transp','nw-cosft','transp',split=[3,'omp']) #Flow('nw-map','nw-transp', # ''' # math output="x1/sqrt(1+0.25*x3*x3*x1*x1/(x2*x2))" | # cut n2=1 # ''') #Flow('nw-fowler','nw-transp nw-map','iwarp warp=\${SOURCES[1]} | transp', # split=[3,'omp']) #Flow('nw-dmo','nw-fowler','cosft sign1=-1 sign3=-1') # Looks like DMO helps to preserve diffractions # dipping flanks from right to left right part of a section #Result('dmo-wat-vel','dmo','window n2=1 min2=1500 | ' + section('DMO 1500 m/s (extracted)')) #Result('dmo-wat-vel-w','dmo','window n2=1 min2=1500 | ' + sectionw('DMO 1500 m/s (extracted)')) # Compute envelope for picking Flow('envelope','dmo','envelope | scale axis=2',split=[3,'omp']) Flow('vpick','envelope',mutec + ' | ' + pickc) #Result('vpick',plotvel('DMO velocity')) # Take a slice Flow('slice','dmo vpick','slice pick=\${SOURCES[1]}') # 'grey title="Nankai DMO Stack" ' Result('slice',section('DMO stack')) # 'window max2=5500 | grey title="DMO stack (extracted)" min1=5.5 max1=7.0' #Result('slice-w','slice',sectionw('DMO Stacked Data')) #>>> #Flow('nw-envelope','nw-dmo','envelope | scale axis=2',split=[3,'omp']) #Flow('nw-vpick','nw-envelope',mutec + ' | ' + pickc) #Result('nw-vpick',plotvel('nw DMO velocity')) #Flow('nw-slice','nw-dmo nw-vpick','slice pick=\${SOURCES[1]}') #Result('nw-slice',section('nw DMO stack')) ### PWD on 1.5 km/s image Flow('vc15','slice', ''' cosft sign2=1 | vczo v0=0.0 nv=1 dv=1500 | window | cosft sign2=-1 ''') #>>> #Flow('nw-vc15','nw-slice', # ''' # cosft sign2=1 | # vczo v0=0.0 nv=1 dv=1500 | # window | # cosft sign2=-1 # ''') #Result('vc15',section("VC 1.5 km/s")) #Result('vc15-w','vc15',sectionw("VC 1.5 km/s")) #Result('vc15-w','vc15','window max2=5500 | grey title="VC 1500 m/s" min1=5.5 max1=7.0') rect1=30 rect2=20 rect3=30 Flow('dip-vc','vc15','dip rect1=%d rect2=%d' % (rect1, rect2)) #Result('dip-vc',plotdip("VC 1.5 km/s Dip")) #Result('dip-vc-w','dip-vc',plotdipw("VC 1.5 km/s Dip")) Flow('pwd-vc','vc15 dip-vc','pwd dip=\${SOURCES[1]}') #Result('pwd-vc-w','pwd-vc',sectionw("VC 1.5 km/s PWD")) #Result('pwd-vc',section("VC 1.5 km/s PWD")) Flow('dif','pwd-vc', ''' cosft sign2=1 | vczo v0=1500.0 nv=1 dv=-1500 | window | cosft sign2=-1 ''') #Result('inv-vc15-w','dif',sectionw("inv VC 1.5 km/s PWD")) #Result('inv-vc15','dif',section("inv VC 1.5 km/s PWD")) Result('dif',section('Separated Diffractions')) #>>> #Flow('nw-dip-vc','nw-vc15','dip rect1=%d rect2=%d' % (rect1, rect2)) #Result('nw-dip-vc',plotdip("nw VC 1.5 km/s Dip")) #Result('nw-dip-vc-w','nw-dip-vc',plotdipw("nw VC 1.5 km/s Dip")) #Flow('nw-pwd-vc','nw-vc15 nw-dip-vc','pwd dip=\${SOURCES[1]}') #Result('nw-pwd-vc-w','nw-pwd-vc',sectionw("nw VC 1.5 km/s PWD")) #Result('nw-pwd-vc',section("nw VC 1.5 km/s PWD")) #Flow('nw-dif','nw-pwd-vc', # ''' # cosft sign2=1 | # vczo v0=1500.0 nv=1 dv=-1500 | # window | # cosft sign2=-1 # ''') #Result('nw-inv-vc15-w','nw-dif',sectionw("nw inv VC 1.5 km/s PWD")) #Result('nw-inv-vc15','nw-dif',section("nw inv VC 1.5 km/s PWD")) ### Lets perform slope decomposition # it is deep not sure if padding is needed pad=1800#3000 Flow('warp','dif','t2warp pad=%i'%(pad)) Plot('dif-spectrum','dif','sffft1 | sfmath output="abs(input)" | sfreal | stack | sfgraph title="Spectrum before warping"') Plot('spectrum','warp','sffft1 | sfmath output="abs(input)" | sfreal | stack | sfgraph title="Spectrum after warping"') #Result('warp','grey title="Warp"') #>>> #Flow('nw-warp','nw-dif','t2warp') #Plot('nw-dif-spectrum','nw-dif','sffft1 | sfmath output="abs(input)" | sfreal | stack | sfgraph title="nw Spectrum before warping"') #Plot('nw-spectrum','nw-warp','sffft1 | sfmath output="abs(input)" | sfreal | stack | sfgraph title="nw Spectrum after warping"') np=151 p0=-0.014 nw=901 #?# -0.014 is it reasonable? #?# do not understand why it is so small #?# dt/dx = 0.5:0.002/1000:16.667 = tgalpha = 5 #?# is it a problem that we cut a dataset? #?# Luke actually windows Barrolka #?# I pad to 3000 and get 3000 instead of 3500 #?# may be should pad more than 3500 - 5000 and window frequencies #?# try the above tests #?# start with 0.005 may be - looks this is sufficient #?# might make sence to perform v0 migration before slope decomposition #?# it should reduce p0 FPX('fpx','warp',np=np,p0=p0,nw=nw,v0=0.0) Flow('txp','fpx','fft1 inv=y | t2warp inv=y | transp plane=23 memsize=1000') #Result('txp',section("TXP")) #Result('tx','txp','stack axis=3 | ' + section("TX")) #>>> looks like slope decomposition does not care if we cut data or not #FPX('nw-fpx','nw-warp',np=np,p0=p0,nw=1441,v0=0.0) #Flow('nw-txp','nw-fpx','fft1 inv=y | t2warp inv=y | transp plane=23 memsize=1000') #Result('nw-txp',section("nw TXP")) #Result('nw-tx','nw-txp','stack axis=3 | ' + section("nw TX")) ### Perform OVC with 1.5 km/s and see if it makes sence on CIGs Flow('fkp','fpx','transp plane=23 memsize=1000 | fft3 axis=2') #>>> #Flow('nw-fkp','nw-fpx','transp plane=23 memsize=1000 | fft3 axis=2') ### trying to use psovc to avert additional transp ### #Flow('f_hall_kp','fpx','fft3 axis=3 | transp plane=24 memsize=30000') #Flow('f_hs_pk','fpx','sfwindow n4=1 f4=%d | fft3 axis=3 | transp plane=24 memsize=1000 | transp plane=34 memsize=1000'%(offset_num)) ### trying to use psovc to avert additional transp ### # creating offset dimension #Flow('fhpk','fpx','spray axis=2 n=1 o=0.0 d=0.1 | fft3 axis=4') #?# had some troubles to parallelize it Flow('ovc-fkp-15','fkp','ovczo nv=%d dv=%g v0=%g' % (1,1500,0.0)) Flow('ovc-15','ovc-fkp-15','window | fft3 axis=2 inv=y | fft1 inv=y | t2warp inv=y | transp plane=23 memsize=5000') #?# looks reasonable and even p0=-0.014 looks to be too wide #Result('ovc-15-cig','ovc-15','window n3=1 min3=2516.7 | grey min1=6.0 max1=8.0 title="OVC 1.5 km/s CIG"') Flow('ovc-15-image','ovc-15','stack') #Result('ovc-15-image',section("OVC 1.5 km/s")) #>>> we still have sharp slope change on the left #Flow('nw-ovc-fkp-15','nw-fkp','ovczo nv=%d dv=%g v0=%g' % (1,1500,0.0)) #Flow('nw-ovc-15','nw-ovc-fkp-15','window | fft3 axis=2 inv=y | fft1 inv=y | t2warp inv=y | transp plane=23 memsize=5000') #Result('nw-ovc-15-cig','nw-ovc-15','window n3=1 min3=2516.7 | grey min1=6.0 max1=8.0 title="nw OVC 1.5 km/s CIG"') ### Performing actual OVC nv = 60 dv = 20 v0 = 1400 Flow('ovc-1st-step','fkp','ovczo nv=%d dv=%g v0=%g | window' % (1,v0,0.0)) #np #Flow('ovc-fkp','ovc-1st-step','ovczo nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[3,'omp'],reduce='cat axis=4') Flow('ovc-fkp','ovc-1st-step','ovczo nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[3,np],reduce='cat axis=4') ### trying to use psovc to avert additional transp - not sure if it is worth it ### #Flow('vc_f_hs_pk_psovc','f_hs_pk','psovc nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[4,1024]) #Flow('vc_f_hall_kp_psovc','f_hall_kp','psovcp nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[4,np]) #?# f-v-k-p -> t-v-x-p, t-x-v-p - because fft3 seems to have a problem when applied to a different axis Flow('txvp','ovc-fkp', ''' transp plane=23 memsize=1000 | fft3 axis=2 inv=y | fft1 memsize=100 inv=y | t2warp inv=y ''',split=[4,np]) # this is where fft3 fails #Flow('tvxp','ovc-fkp','fft3 axis=3 inv=y | fft1 memsize=100 inv=y | t2warp inv=y',split=[4,np]) Flow('ovc','txvp','stack axis=4') #Flow('ovc','ovc-fkp','stack axis=4 | fft3 axis=3 inv=y | fft1 inv=y | t2warp inv=y') #Result('ovc',section("OVC")) #?# sharp boundary on CIG - zero on the left (from around +/-0.005) and non-zero on the right #Result('ovc-cig','txvp','window n2=1 min2=2516.7 | transp plane=23 | grey min1=6.0 max1=8.0 title="OVC CIG"') ### Semblance calculation Flow('txv1','txvp','stack norm=n axis=4') Flow('txv2','txvp','mul \$SOURCE | stack norm=n axis=4') #Result('txvp', # ''' # window n3=1 f3=1 | # pow pow1=1 | # byte | # grey3 flat=n frame1=250 frame2=250 frame3=100 # ''') #Plot('txv1','grey gainpanel=all title="Velocity Continuation" ',view=1) #Plot('txv2','grey gainpanel=all title="What is txv2 check" ',view=1) # optimal rect1=5 rect2=5 # rect1=10 Flow('sem','txv1 txv2', ''' mul \${SOURCES[0]} | divn den=\${SOURCES[1]} rect1=5 rect2=5 niter=25 | clip2 lower=0 ''') # I dont know how to automate it #Flow('max.txt','sem','attr want=max | awk "{print \$2}" --out=\$TARGET',stdout=0) #Flow('max',None,'spike mag=`attr