sconstruct

国外免费地震资料处理软件包

from rsfproj import *
import math, sys

# Velocity model
########################################
v0 = 0.5        # velocity at the origin
gz = 0.8        # vertical gradient
gx = 0.6        # horizontal gradient
nz = 601        # depth samples
nx = 1201       # lateral samples
dz = 0.5/(nz-1) # depth sampling
dx = 1.0/(nx-1) # lateral sampling

# v(z,x) = v0 + gz*z + gx*x
###########################
Flow('vel',None,
     '''
     math n1=%d d1=%g n2=%d d2=%g
     output="%g+%g*x1+%g*x2"
     ''' % (nz,dz,nx,dx,v0,gz,gx))

# v(z,x) = sqrt(v0^3/(v0 + 2*gz*z + 2*gx*x))
########################################

# !!!!!!!!!! UNCOMMENT BELOW !!!!!!!!!!

## Flow('vel',None,
##        '''
##        math n1=%d d1=%g n2=%d d2=%g
##        output="%g*sqrt(%g/(%g-2*%g*x1-2*%g*x2))"
##        ''' % (nz,dz,nx,dx,v0,v0,v0,0.125*gz,0.125*gx))

Plot('vel',
     '''
     grey color=j allpos=y screenht=6 screenratio=%g
     bias=%g wanttitle=n barreverse=y
     label1=Depth unit1=km label2=Lateral unit2=km
     scalebar=y barlabel=Velocity barunit="km/s" 
     ''' % ((nz-1.)/(nx-1.),v0))

# Analytical traveltime
########################

# !!!!!!!!!! CHANGE BELOW !!!!!!!!!!!!!

g = math.hypot(gx,gz) # sqrt(gx*gx + gz*gz)
Flow('time','vel',
     '''
     math output="acosh(1 + %g*(x1*x1+x2*x2)/input)/%g"
     ''' % (0.5*g*g/v0,g))
Plot('time',
     '''
     contour scalebar=y plotcol=7 c0=0 dc=0.1 nc=30
     screenht=6 screenratio=%g plotfat=5 wanttitle=n
     ''' % ((nz-1.)/(nx-1.)))

# Plot wavefronts overlayed on the velocity model
#################################################
Result('exact','vel time','Overlay')

errs = []
cpus = []
sizes = []
spaces = []

samples = range(1,nz/30)

# !!!!!!!!! CHANGE BELOW !!!!!!!!!!

eikonal = '%s yshot=0 order=1 br1=%g br2=%g' % \
          (WhereIs('sfeikonal'),20*dz,20*dx)

if sys.platform[:5] == 'sunos':
    time = '''
    (/usr/bin/time %s < $SOURCE > ${TARGETS[1]} ) 2> cpu.out &&
    (awk 'NR < 5 && NR > 2 {print $2}' cpu.out;
    ''' % eikonal
else:
    time = '''
    (/usr/bin/time -f "%%S %%U" %s < $SOURCE > ${TARGETS[1]}
    ) >& cpu.out && (tail -1 cpu.out;
    ''' % eikonal

for sample in samples:

    # Subsample velocity and analytical traveltime
    ##############################################
    vel = 'vel%d' % sample
    tim = 'tim%d' % sample
    window = 'window j1=%d j2=%d' % (sample,sample)
    
    Flow(vel,'vel',window)
    Flow(tim,'time',window)

    size = 'siz%d' % sample
    n = float(nx*nz)/(sample*sample)
    Flow(size,None,'spike n1=1 mag=%g' % n)
    sizes.append(size)

    space = 'spc%d' % sample
    h = sample*math.hypot(dx,dz)
    Flow(space,None,'spike n1=1 mag=%g' % h)
    spaces.append(space)


    # Solve the eikonal equation, record CPU time
    #############################################
    eik = 'eik%d' % sample
    cpu = 'cpu%d' % sample
    
    Flow([cpu,eik],vel,time + \
         '''
         echo in=cpu.rsf n1=2 data_format=ascii_float)
         > cpu.rsf && dd form=native < cpu.rsf |
         stack axis=1 > ${TARGETS[0]} &&
         /bin/rm cpu.out cpu.rsf
         ''',
         stdin=0,stdout=-1)
    cpus.append(cpu)

    # Compute error
    ###############
    err = 'err%d' % sample
    Flow(err,[eik,tim],
         '''
         math ref=${SOURCES[1]} output="abs(input-ref)" |
         stack axis=2 norm=y | stack axis=1 norm=y 
         ''')
    errs.append(err)

# Concatenate results
#####################
cat = 'cat axis=1 ${SOURCES[1:%d]}' % len(samples)
   
Flow('space',spaces,cat)
Flow('logspace','space','math output="log(input)" ')

Flow('err',errs,cat)
Flow('logerr','err','math output="log(input)" ')

Flow('size',sizes,cat)
Flow('cpu',cpus,cat)

graph = 'cmplx ${SOURCES[0:2]} | graph wanttitle=n '
for case in (0,1):
    Plot('err'+str(case),'space err',graph + '''
    label1="Grid Spacing" unit1=km label2=Error unit2=s
    ''',stdin=0)
    Plot('logerr'+str(case),'logspace logerr',graph + '''
    label1="Log(Grid Spacing)" unit1= label2="Log(Error)" unit2=
    ''',stdin=0)
    Plot('cpu'+str(case),'size cpu',graph + '''
    label1="Grid Size" unit1= label2="CPU Time" unit2=s
    ''',stdin=0)
    graph = graph + 'symbol=x plotcol=2 '

# Plot error versus spacing
###########################
Plot('err','err0 err1','Overlay')

# Plot error versus spacing on a logarithmic scale
##################################################
Plot('logerr','logerr0 logerr1','Overlay')

Result('err','err logerr','SideBySideAniso')

# Plot CPU time
###############
Result('cpu','cpu0 cpu1','Overlay')

End()