```## # Angle-gathers after wave-equation migration (ziggy model) ## from rsf.proj import * import pplot # ------------------------------------------------------------ # MIGRATION parameters # ------------------------------------------------------------ par = { 'ox':10925,'dx':150, 'nx':500, 'oz':6100, 'dz':25, 'oh':-3000,'dh':150 } # ------------------------------------------------------------ # from ft to km ft2km = 0.3024/1000 par['ox']=par['ox']*ft2km par['dx']=par['dx']*ft2km par['oz']=par['oz']*ft2km par['dz']=par['dz']*ft2km par['oh']=par['oh']*ft2km par['dh']=par['dh']*ft2km # ------------------------------------------------------------ par['zmin']=2 par['zmax']=9 par['xmin']=4 par['xmax']=21 IRATIO = ' labelsz=6 screenratio=0.5 screenht=7 min2=%g max2=%g ' % (par['xmin'],par['xmax']) CRATIO = ' labelsz=5 screenratio=3.5 screenht=10.5 ' LABEL = " wantaxis=n " # ------------------------------------------------------------ def igrey(custom,par): return ''' grey labelrot=n wantaxis=y wanttitle=n grid=y gridcol=6 title="" pclip=97 label1="z" unit1=km label2="x" unit2=km min1=%g max1=%g %s ''' % (par['zmin'],par['zmax'],custom) def p1x3(plot,p0,p1,p2,ys,xs,x0,x1,x2): j0 = '_' + plot+ p0 j1 = '_' + plot+ p1 j2 = '_' + plot+ p2 Plot(j0,p0,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x0)) Plot(j1,p1,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x1)) Plot(j2,p2,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x2)) Plot (plot,[j0,j1,j2],'Overlay') Result(plot,[j0,j1,j2],'Overlay') def p1x4(plot,p0,p1,p2,p3,ys,xs,x0,x1,x2,x3): j0 = '_' + plot+ p0 j1 = '_' + plot+ p1 j2 = '_' + plot+ p2 j3 = '_' + plot+ p3 Plot(j0,p0,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x0)) Plot(j1,p1,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x1)) Plot(j2,p2,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x2)) Plot(j3,p3,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x3)) Plot (plot,[j0,j1,j2,j3],'Overlay') Result(plot,[j0,j1,j2,j3],'Overlay') def p1x6(plot,p0,p1,p2,p3,p4,p5,ys,xs,xc): j0 = '_' + plot+ p0 j1 = '_' + plot+ p1 j2 = '_' + plot+ p2 j3 = '_' + plot+ p3 j4 = '_' + plot+ p4 j5 = '_' + plot+ p5 Plot(j0,p0,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs, 0-1)) Plot(j1,p1,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs, xc-1)) Plot(j2,p2,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,2*xc-1)) Plot(j3,p3,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,3*xc-1)) Plot(j4,p4,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,4*xc-1)) Plot(j5,p5,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,5*xc-1)) Plot (plot,[j0,j1,j2,j3,j4,j5],'Overlay') Result(plot,[j0,j1,j2,j3,j4,j5],'Overlay') # ------------------------------------------------------------ # ------------------------------------------------------------ # from SGY to RSF # ------------------------------------------------------------ velo = 'sigsbee2a_migvel.sgy' vstr = 'sigsbee2a_stratigraphy.sgy' Fetch(velo,'sigsbee') Fetch(vstr,'sigsbee') Flow('zvelo tzvelo ./vhead ./bvhead',velo, ''' segyread tape=\$SOURCE tfile=\${TARGETS[1]} hfile=\${TARGETS[2]} bfile=\${TARGETS[3]} ''',stdin=0) Flow('zvstr tzvstr ./shead ./bshead',vstr, ''' segyread tape=\$SOURCE tfile=\${TARGETS[1]} hfile=\${TARGETS[2]} bfile=\${TARGETS[3]} ''',stdin=0) # ------------------------------------------------------------ # VELOCITY # ------------------------------------------------------------ # stratigraphic slowness Flow('vstr','zvstr', ''' scale rscale=%g | put d1=%g o2=%g d2=%g label1=z label2=x '''% (ft2km,25*ft2km, 10025*ft2km, 25*ft2km)) Flow('vstrpad','vstr', ''' window n1=1 f1=1200 | spray axis=1 n=100 o=0 d=%g ''' % (25*ft2km)) Flow('s0','vstr vstrpad', ''' cat axis=1 space=n \${SOURCES[1]} | window | math "output=1/input" | transp | spray axis=2 n=1 o=0 d=1 | put label2=y | window squeeze=n f3=244 ''' % par) # migration slowness Flow('velo','zvelo', ''' scale rscale=%g | put d1=%g o2=%g d2=%g label1=z label2=x '''% (ft2km,25*ft2km, 10025*ft2km, 37.5*ft2km)) Flow('velopad','velo', ''' window n1=1 f1=1200 | spray axis=1 n=100 o=0 d=%g ''' % (25*ft2km)) Flow('s1','velo velopad', ''' cat axis=1 space=n \${SOURCES[1]} | window | math "output=1/input" | transp | spray axis=2 n=1 o=0 d=1 | put label2=y | window squeeze=n f3=244 ''') # wrong migration velocity Flow('sm','s1', ''' window squeeze=n f3=350 | math output=0.1 | pad beg3=350 | math output=input+1 | smooth rect3=31 ''') Flow('s2','s1 sm', 'math s=\${SOURCES[0]} m=\${SOURCES[1]} output=s*m') # ------------------------------------------------------------ # MIGRATION # ------------------------------------------------------------ SLOlist = '0','2' CIGlist = '7','9','11','13','15','17' IMGlist = 'x','t' cigpar = {'npt':500,'opt':4000,'dpt':80, 'na':400, 'oa':-2.0, 'da':0.01, 'mina':0, 'maxa':1.3 } cigpar['opt']=cigpar['opt']*ft2km cigpar['dpt']=cigpar['dpt']*ft2km # ------------------------------------------------------------ # CIG to slant-stack def tcig2ssk(): return ''' slant adj=y p0=%g np=%d dp=%g ''' % (cigpar['opt'],cigpar['npt'],cigpar['dpt']) def hcig2ssk(): return ''' slant adj=y p0=%g np=%d dp=%g ''' % (cigpar['oa'],cigpar['na'],cigpar['da']) # slant-stack to angle def tssk2ang(): return ''' tshift cos=y a0=0 na=150 da=0.45 velocity=\${SOURCES[1]} dip=\${SOURCES[2]} ''' def hssk2ang(): return ''' tan2ang a0=0 na=150 da=0.45 ''' # ------------------------------------------------------------ for i in SLOlist: for k in IMGlist: SRcig = 'SRC' + i + k Fetch( SRcig+'.hh','sigsbee') if(k=='t'): Flow(SRcig,SRcig+'.hh', ''' window squeeze=n | put o1=%(ox)g d1=%(dx)g o3=%(oz)g d3=%(dz)g ''' % par) if(k=='x'): Flow(SRcig,SRcig+'.hh', ''' window squeeze=n | put o1=%(ox)g d1=%(dx)g o3=%(oz)g d3=%(dz)g o4=%(oh)g d4=%(dh)g ''' % par) if(0): cig = 'cig' + i + k Flow(cig,SRcig,'transp plane=15 memsize=500 | window') env = 'env' + i + k Flow(env,cig,'envelope') pik = 'pik' + i + k Flow(pik,env,'pick rect1=10 rect2=100') slc = 'slc' + i + k Flow(slc,[cig,pik],'slice pick=\${SOURCES[1]} | window') Result(slc,igrey('pclip=95'+IRATIO+LABEL,par)) # ------------------------------------------------------------ # RESULTS # ------------------------------------------------------------ for i in SLOlist: slo = 's' + i Plot (slo,'window | transp |' +igrey('pclip=99 allpos=y bias=0.22'+IRATIO+LABEL,par)) Result(slo,'window | transp |' +igrey('pclip=99 allpos=y bias=0.22'+IRATIO+LABEL,par)) vel = 'v' + i Flow(vel,slo,'window | transp | smooth rect1=51 | math output=1/input') for k in IMGlist: # image / CIGs SRimg = 'SRI' + i + k SRcig = 'SRC' + i + k Flow(SRimg,SRcig,'window n4=1 min4=0 | transp') Plot (SRimg,igrey('pclip=95'+IRATIO+LABEL,par)) Result(SRimg,igrey('pclip=95'+IRATIO+LABEL,par)) IMGSLO = 'IMGSLO' + i + k pplot.p2x1(IMGSLO,SRimg,slo,0.75,0.75,-6.5) # dip field SRdip = 'DIP' + i + k Flow(SRdip,SRimg,'dip rect1=3 rect2=3 order=3 liter=100') Result(SRdip,igrey('pclip=99 color=j'+IRATIO+LABEL,par)) SRmap = 'MAP' + i + k Flow(SRmap,[SRdip,slo],'transp | remap1 pattern=\${SOURCES[1]} | transp') # image after plane-wave distruction SRpwd = 'PWD' + i + k Flow(SRpwd,[SRimg,SRdip],'pwd dip=\${SOURCES[1]}') Result(SRpwd,igrey('pclip=99'+IRATIO+LABEL,par)) # dip-corrected velocity vco = 'v' + i + k Flow(vco,[SRmap,vel], ''' math v=\${SOURCES[1]} output="v/cos(atan(input*%(dz)g/%(dx)g))" ''' % par) for l in CIGlist: ikl = i + k + l # image z = 'z' + ikl lmin = float(l) - 20 * float('%(dx)g' % par) lmax = float(l) + 20 * float('%(dx)g' % par) Flow(z,SRcig, ''' window n4=1 min4=0 min1=%g max1=%g | transp | bandpass fhi=10 ''' % (lmin,lmax) ) Plot(z,z,igrey('labelsz=5 screenratio=6 screenht=10.5'+LABEL,par)) # velocity v = 'v' + ikl Flow(v,vco,'window n2=1 min2=%g' % float(l) ) # structural dip d = 'd' + ikl Flow(d,SRdip,'window n2=1 min2=%g' % float(l) ) # offset gathers h = 'h' + ikl Flow(h,SRcig, ''' window n1=1 min1=%g | bandpass fhi=10 ''' % float(l)) # angle gathers a = 'a' + ikl # slant-stack c = 'c' + ikl # angle-gather cdip = 'cdip' + ikl # ccln = 'ccln' + ikl # dip noise separation if(k=='x'): Flow(a,h,'bandpass fhi=10 | '+hcig2ssk()) Flow(c,a, hssk2ang()) Plot(h,igrey('label2="h" unit2=km' +CRATIO+LABEL,par)) Plot(a,igrey('label2="tan" unit2="\F10 q\F3 " min2=0'+CRATIO+LABEL,par)) if(k=='t'): Flow(a,h,'bandpass fhi=10 | '+tcig2ssk()) Flow(c,[a,v,d], tssk2ang()) Plot(h,igrey('label2="\F10 t\F3 " unit2=s ' +CRATIO+LABEL,par)) Plot(a,igrey('label2="\F10 n\F3 " unit2="km/s" max2=5'+CRATIO+LABEL,par)) # clean-up noise Flow(cdip,c,'twodip2 eps=10 lam=1 p0=-1 q0=0 both=n') Flow(ccln,[c,cdip],'pwdsigk dips=\${SOURCES[1]} verb=y eps=1 | window n3=1 f3=1') Plot (c ,igrey('label2="\F10 q\F3 " unit2="\^o\_"' +CRATIO+LABEL,par)) Plot (ccln,igrey('label2="\F10 q\F3 " unit2="\^o\_"' +CRATIO+LABEL,par)) # < clean0t11.rsf sfwindow f3=1 > clean.rsf # < clean.rsf sfenvelope | sfthreshold pclip=100 | sfmask min=0.1 | sfdd type=float | sfadd mode=p clean.rsf | sfstack | sfgraph | tube # ------------------------------------------------------------ for i in SLOlist: for l in CIGlist: SRoff = 'SRoff' + i + '-' + l SRang = 'SRang' + i + '-' + l SRxic = 'SRxic' + i + '-' + l SRtic = 'SRtic' + i + '-' + l zx = 'z' + i + 'x' + l zt = 'z' + i + 't' + l hx = 'h' + i + 'x' + l ht = 'h' + i + 't' + l ax = 'a' + i + 'x' + l at = 'a' + i + 't' + l cx = 'ccln' + i + 'x' + l ct = 'ccln' + i + 't' + l p1x3(SRxic,hx,ax,cx,1,1,-1,-4,-7) p1x3(SRtic,ht,at,ct,1,1,-1,-4,-7) p1x3(SRoff,zx,hx,ht,1,1,-1,-3,-6) p1x3(SRang,zx,cx,ct,1,1,-1,-3,-6) SRt = 'SRt' + i + '-' + l p1x4(SRt,zt,ht,at,ct,1,1,-1,-3,-6,-9) SRx = 'SRx' + i + '-' + l p1x4(SRx,zx,hx,ax,cx,1,1,-1,-3,-6,-9) # SRic = 'SRic' + i + '-' + l # pplot.p2x1(SRic,SRxic,SRtic,0.5,0.5,-9) # SRcig = 'SRcig' + i + '-' + l # pplot.p2x1(SRcig,SRoff,SRang,0.5,0.5,-9) p1x6('allxoff', 'h0x7', 'h0x9', 'h0x11', 'h0x13', 'h0x15', 'h0x17',0.75,0.75,-3) p1x6('allxang','ccln0x7','ccln0x9','ccln0x11','ccln0x13','ccln0x15','ccln0x17',0.75,0.75,-3) p1x6('alltoff', 'h0t7', 'h0t9', 'h0t11', 'h0t13', 'h0t15', 'h0t17',0.75,0.75,-3) p1x6('alltang','ccln0t7','ccln0t9','ccln0t11','ccln0t13','ccln0t15','ccln0t17',0.75,0.75,-3) pplot.p2x1('alloff','allxoff','alltoff',0.65,0.65,-7) pplot.p2x1('allang','allxang','alltang',0.65,0.65,-7) # ------------------------------------------------------------ End()```

 sfsegyread sfscale sfput sfwindow sfspray sfcat sfmath sftransp sfpad sfsmooth sfgrey sfdip sfremap1 sfpwd sfbandpass sfslant sftan2ang sftwodip2 sfpwdsigk sftshift