from rsf.proj import * from findradius1 import findradius1 from rsf.recipes.beg import server ############################################################# # Load Data and First Look ############################################################# # Fetch Data from local directory #Fetch(['i455_legacy3d_pstm_4ms.sgy','i455_sc2dpoststkmig_1ms.sgy'], # 'data',top='',server='local') # Fetch Data data = {'legacy':'i455_legacy3d_pstm_4ms.sgy', 'hires':'i455_sc2dpoststkmig_1ms.sgy'} for key in data.keys(): Fetch(data[key],'apache',server) # Convert from SEGY format to Madgascar RSF format Flow(['legacy','tlegacy','legacy.asc','legacy.bin'], 'i455_legacy3d_pstm_4ms.sgy', ''' segyread tfile=${TARGETS[1]} hfile=${TARGETS[2]} bfile=${TARGETS[3]} | put label2=Cross-line o2=6638 d2=1 ''') Flow(['hires','thires','hires.asc','hires.bin'], 'i455_sc2dpoststkmig_1ms.sgy', ''' segyread tfile=${TARGETS[1]} hfile=${TARGETS[2]} bfile=${TARGETS[3]} | put label2=Cross-line o2=6638 d2=1 ''') # Display Data Result('hires','grey color=seismic title="Hires Image" ') Result('legacy','grey color=seismic title="Legacy Image" ') ############################################################# # Subsampling ############################################################# # Initial # Hires: dt=1ms , ns=4001 # Legacy: dt=4ms , ns=1001 # Sample both at dt=2ms, ns=2001 Flow('hires2','hires','bandpass fhi=250 | window j1=2') # Subsampling Flow('legacy2','legacy','spline n1=2001 d1=0.002 o1=0') # interpolation # interpolate legacy: dt=1ms , ns=4001 Flow('legacy4','legacy','spline n1=4001 d1=0.001 o1=0') Result('legacy4','grey color=seismic title="Legacy Image" ') ############################################################# # Spectra ############################################################# Flow('legacy2-spec','legacy2','spectra all=y') Flow('hires2-spec','hires2','spectra all=y') Flow('legacy-spec','legacy','spectra all=y') # display spectra Result('spectra','hires2-spec legacy2-spec', 'cat axis=2 ${SOURCES[1]} | graph title=Spectra') # display normalized spectra Result('nspectra','hires2-spec legacy2-spec', ''' cat axis=2 ${SOURCES[1]} | scale axis=1 | window max1=130 | graph title="Normalized Spectra" label2="Amplitude" ''') ############################################################# # BANDPASS FILTERING ############################################################# # Filter frequencies (Filter low at 18 and high at 38) lof = 18 # other value: 21 hif = 38 # other value: 35 Flow('hires-filt', 'hires2', 'bandpass fhi='+str(hif)+' flo='+str(lof)) Flow('legacy-filt', 'legacy2', 'bandpass fhi='+str(hif)+' flo='+str(lof)) Result('hires-filt', 'grey color=seismic title="Inline hires filtered"') Result('legacy-filt', 'grey color=seismic title="Inline legacy filtered"') # Spectra Flow('legacy-filt-spec', 'legacy-filt', 'spectra all=y') Flow('hires-filt-spec', 'hires-filt', 'spectra all=y') # Display spectra Result('filtspectra','hires-filt-spec legacy-filt-spec', ''' cat axis=2 ${SOURCES[1]} | graph title="Filtered Spectra" ''') # Display normalized spectra Result('filtnspectra','hires-filt-spec legacy-filt-spec', ''' cat axis=2 ${SOURCES[1]} | scale axis=1 | graph title="Filtered Normalized Spectra" ''') ################################################################# # NON-STATIONARY SMOOTHING APPLIED TO HIRES - theoretical radius ################################################################# # measure local frequency for key in ['legacy','hires']: freq = key+'-locfreq' Flow(freq,key,'iphase order=10 rect1=%d rect2=80 hertz=y complex=y | put label=Frequency unit=Hz' % (80,80)[key=='hires']) Result(freq,'grey mean=y color=j scalebar=y title="%s Local Frequency" ' % key.capitalize()) #clip=35 bias=35 # Absolute Difference in local frequencies Flow('freqdif','legacy-locfreq hires-locfreq','remap1 n1=4001 d1=0.001 o1=0 | add scale=-1,1 ${SOURCES[1]} | math output="abs(input)"') Result('freqdif','grey color=j minval=0 maxval=50 bias=25 scalebar=y title="initial difference in local frequency" ') # non-stationary radius calculated theoretically from math import pi scale=6.75 Flow('rect','legacy-locfreq hires-locfreq', ''' remap1 n1=4001 d1=0.001 o1=0 | math f1=${SOURCES[1]} output="sqrt(%g*(1/(input*input)-1/(f1*f1)))/%g" ''' % (scale,2*pi*0.001)) Result('rect','grey color=viridis mean=y title="Smoothing Radius" scalebar=y barlabel=Radius barunit=samples') # smooth hires with non-starionary radius Flow('hires-smooth','hires rect','nsmooth1 rect=${SOURCES[1]} | window j1=4') Flow('hires-smooth-2','hires rect','nsmooth1 rect=${SOURCES[1]}') Result('hires-smooth', 'grey color=seismic title="Hires Smooth"') # get spectra Flow('hires-smooth-spec','hires-smooth','spectra all=y') Result('hires-smooth-spec','hires-smooth-spec legacy-spec', ''' cat axis=2 ${SOURCES[1]} | scale axis=1 | window max1=120 | graph title="Normalized Spectra after Smoothing" label2="Amplitude" ''') # local frequency of smoothed hires Flow('hires-smooth-locfreq','hires-smooth','iphase order=10 rect1=40 rect2=80 hertz=y complex=y | put label=Frequency unit=Hz') Result('hires-smooth-locfreq','grey mean=y color=j scalebar=y title="Hires Local Frequency Smoothed" ') # Difference in local frequencies Flow('freqdif-filt','legacy-locfreq hires-smooth-locfreq','add scale=-1,1 ${SOURCES[1]} | math output="abs(input)"') Result('freqdif-filt','grey color=j minval=0 maxval=50 bias=25 scalebar=y title="Difference after Smoothing" ') ############################################################# # non-stationary radius estimation ############################################################# D1 = 'hires' D2 = 'legacy4' niter = 5 recmin = 60 rec1 = recmin*(2**(niter-3)) rec2 = recmin*(2**(niter-3)) findradius1(D1, D2, # original data (high) and filtered data (low) niter=niter, # number of corrections rec1=rec1, rec2=rec2, # smooth division parameters recmin=recmin, # smooth division minimum radius (for decreasing radius) recmax=1000, # smooth division radius 1st iteration (choose large value to estimate the stationary radius) rec1f=80, rec2f=80, # smooth local frequency parameters maxvalf=50, # maximum value for plotting local frequency difference biasf=25, # bias for plotting local frequency difference theor=True, # use smoothed version of theoretical smoothing radius as starting model rectheor=700, # smoothing radius for theoretical radius scale=scale, # scale for theoretical smoothing radius d1=0.001, # interval spacing in 1st dimension (time) initial=25, # initial value for constant smoothing radius titlehigh="High-Resolution", # original data (high) title titlelow="Legacy", # filtered data (low) title localfreq=False, # estimate based on local frequency attribute it=0) # correction number End() |
sfsegyread sfput sfgrey sfbandpass sfwindow | sfspline sfspectra sfcat sfgraph sfscale | sfiphase sfremap1 sfadd sfmath sfnsmooth1 | sfsmooth sfnsmooth sfndsmooth sfdivn |