sconstruct

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

##
 # Sigsbee 2A
 # Wave-equation migration: all shots
 ##

from rsfproj import *
import spmig

# ------------------------------------------------------------
# MIGRATION parameters
# ------------------------------------------------------------
par = {
    'nx':500,'ox':10925,'dx':150,
    'nt':1500,'ot':0,'dt':0.008,'kt':1,'ft':0,
#
    'nw':360,'ow':1, 'jw':1,
    'nw':1,  'ow':5, 'jw':1,
#
    'ns':100, 'js':5, 'fs':0,
    'ns':1,   'js':1, 'fs':250,
#
    'verb':'y','eps':0.01,'nrmax':5,'dtmax':0.00005,
    'tmx':16,'pmx':0,
    #
    'zmin':6100, 'zmax':30000
    }
par['xmin']=par['ox']
par['xmax']=par['ox'] + (par['nx']-1) * par['dx']
minmax = ' min2=%g max2=%g ' % (par['xmin'],par['xmax'])

def igrey(custom,par):
    return '''
    grey labelrot=n wantaxis=y wanttitle=n labelsz=6
    min2=%g max2=%g
    title="" pclip=97 label1="z(ft)" label2="x(ft)" %s
    ''' % (par['xmin'],par['xmax'],custom)
LABEL = " labelsz=12 titlesz=16 "
# ------------------------------------------------------------

# ------------------------------------------------------------
# read SEGy files
# convert to RSF format
# ------------------------------------------------------------
data = 'sigsbee2a_nfs.sgy'
vstr = 'sigsbee2a_stratigraphy.sgy'

Fetch(data,'sigsbee')
Fetch(vstr,'sigsbee')

Flow('zdata tzdata ./dhead ./bdhead',data,
     '''
     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)

# ------------------------------------------------------------
# SLOWNESS
# prepare slowness files (3-D format: x-y-z)
# ------------------------------------------------------------
Flow('vstr','zvstr',
     'put d1=25 o2=10025 d2=25 label1=z label2=x')

# migration slowness
Flow('slo','vstr',
     '''
     window |
     math "output=1/input" |
     transp |
     spray axis=2 n=1 o=0 d=1 |
     put label2=y |
     window squeeze=n      f3=244
     ''' % par)

Result('slo','window      | transp |'
       + igrey('title=s pclip=100 color=j allpos=y',par))


# datuming slowness
Flow('sdt','vstr',
     '''
     window |
     math "output=1/input" |
     transp |
     spray axis=2 n=1 o=0 d=1 |
     put label2=y |
     window squeeze=n n3=2 j3=244
     ''')

# ------------------------------------------------------------
# SHOTS
# prepare shots file (3-D format: t-o-s)
# ------------------------------------------------------------
Flow('ss','tzdata','dd type=float | headermath output="10925+fldr*150" | window')
Flow('oo','tzdata','dd type=float | headermath output="offset"         | window')
Flow('si','ss','math output=input/150')
Flow('oi','oo','math output=input/75')
Flow('os','oi si','cat axis=2 space=n ${SOURCES[1]} | transp | dd type=int')
Flow('sraw','zdata os',
     '''
     intbin head=${SOURCES[1]} xkey=0 ykey=1
     ''')
Flow('shot','sraw',
     '''
     put d2=75 d3=150 o3=10925 label1=t label2=o label3=s |
     mutter half=false t0=1.0 v0=6000
     ''')

# ------------------------------------------------------------
# WAVELET
# prepare source wavelet file (1-D format: t)
# ------------------------------------------------------------
Flow('wave',None,
     '''
     spike nsp=1 mag=1 k1=%(kt)d
     n1=%(nt)d d1=%(dt)g o1=%(ot)g
         ''' % par )    
Result('wave','wave','graph title="wavelet"')
# ------------------------------------------------------------

# ------------------------------------------------------------
# SHOT-RECORD migration
# ------------------------------------------------------------

# loop over migration configurations
for i in range(2):

    loc = par.copy()
    if  (i==0):
        loc['ow']=2
    elif(i==1):
        loc['nw']=1
        loc['ns']=3
        loc['fs']=50
        loc['js']=150

    # prepare source/receiver data
    dos = 'dos' + str(i)
    dor = 'dor' + str(i)
    spmig.wflds(dos,dor,'wave','shot',loc) # prepare S-R wavefields

    # datum source/receiver data in water
    dds = 'dds' + str(i)
    ddr = 'ddr' + str(i)
    spmig.datum(dds,ddr,'sdt',dos,dor,loc) #   datum S-R wavefields

    # migration
    img = 'img' + str(i)  # image (3-D format: x-y-z)
    cig = 'cig' + str(i)  # CIGs
    spmig.imagePW(img,cig,'slo',dds,ddr,loc)

    Result(img,'window | transp |'
           + igrey('title=i pclip=99',par))
        
# ------------------------------------------------------------


End()