sconstruct

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

# Inversion for von Karman spatial autocorrelation parameters in 2-D Fourier domain
#
# January 2007
#
# Thomas Jules Browaeys 
# Bureau of Economic Geology
# University of Texas at Austin
# mailto:[e-mail]


from rsfproj import *
from math import pi


# ---------------
# Spatial 2D grid
# ---------------

# Master parameters
#
# nx,ny = number of points in X,Y
# dx,dy = space data step sampling

# Slave parameters (spatial frequency content)
#
# fpx = 1/dx        = frequency space periodicity (Hz)
# lx  = nx*dx       = space period
# dfx = 1/lx        = frequency step
# fmx = 1/dx - 1/lx = maximum frequency
# fnx = 1/(2*dx)    = nx/(2*lx) = Nyquist frequency
#
# Signal X detecting content (Hz) =  dfx < fx < fnx
# Signal Y detecting content (Hz) =  dfy < fy < fny


pgrid = {'nx':256, 'ox':0., 'dx':1., 'ny':256, 'oy':0., 'dy':1.}


# Create spatial grid


Flow('spacegrid',None,
     '''
     spike nsp=1 mag=1 n1=%(nx)d d1=%(dx)g o1=%(ox)g n2=%(ny)d d2=%(dy)g o2=%(oy)g |
     put label1=x unit1=m label2=y unit2=m
     ''' % pgrid)


# ---------------------
# Synthetic 2-D surface
# ---------------------

# bx,by   = correlation length scale in x,y (m)
# nu      = Hurst exponent -0.5 < and < 0.5
# gmu     = Gaussian white noise mean
# gvr     = variance of Gaussian white noise = (standard deviation)^2
# grn     = noise range
# gsd     = seed for random generator


psyn = {'bx':10., 'by':100., 'nu':0.4, 'gmu':0., 'gvr':9., 'grn':1., 'gsd':1147}


# Create Gaussian white noise


Flow('wgauss','spacegrid','noise mean=%(gmu)g range=%(grn)g rep=y seed=%(gsd)g type=y var=%(gvr)g' % psyn)
Plot('wgauss','sfgrey color=j title="Gaussian white noise"')


# Discrete Fourier Transform (fft) (k=0,N-1)
#
# S(k*DF) = DL*sum(n=0,N-1) s(n*DL)*exp(-2*i*pi*k*n/N)
#
# s  = spatial signal
# S  = frequency signal
# N  = points number in space
# DF = frequency sampling step = 1/L
# L  = spatial length of signal = 1/DF
# DL = space sampling step = L/N = 1/(N*DF)
# FM = maximum frequency = 1/DL - 1/L
#
# Properties
#
# Symmetry = S(-u,-v) = S*(u,v)
# Periodicity = S(ix+NX,iy+NY) = S(ix,iy)
# Physical vectors (k=0,N-1)
# S(N-1) = S(-1) = S*(1)
# S(N-2) = S(-2) = S*(2)
# ...
# S(N/2+1) = S(-N/2-1) = S*(N/2-1)
# S(N/2)   = S(-N/2)   = S*(N/2)
#
# Dimension in Fourier space [0,nx/2][0,ny/2] = (nx/2+1)*(ny/2+1)
# [fx,fy] = dfx*(-nx/2+1:1:nx/2), dfy*(-ny/2+1:1:ny/2)


fft2 = 'sffft1 sym=y | sffft3 sym=y'
Flow('fwgauss','wgauss',fft2)
Plot('fwgauss',
     '''
     sfadd abs=y | sfreal | put label1=kx unit1=1/m label2=ky unit2=1/m |
     sfgrey color=j title="Gaussian white noise spectrum"
     ''')


# Amplitude of white spectrum


Flow('mfwgauss','fwgauss',
     '''
     add abs=y | real |
     stack axis=2 norm=n | math output="input*%g" |
     stack axis=1 norm=n | math output="input*%g" |
     spray axis=1 n=%d d=%g o=0 | spray axis=2 n=%d d=%g o=%g 
     ''' % (0.5/pgrid['ny'],1./(pgrid['nx']/2+1),pgrid['nx']/2+1,1./pgrid['nx'],2*pgrid['ny'],0.5/pgrid['ny'],-0.5/pgrid['dy']))


# Stochastic process von Karman 2-D filter in spectral domain (Lord, 1954)


Flow('vkfilt','fwgauss','math output="(1+(%(bx)g*x1)^2+(%(by)g*x2)^2)^(-0.5-0.5*%(nu)g)"' % psyn)
Plot('vkfilt',
     '''
     sfreal | put label1=kx unit1=1/m label2=ky unit2=1/m |
     sfgrey color=j title="Filter spectrum F(kx,ky)"
     ''')


# Logarithmic plot for separable least square inversion method ("Conical Coolie Hat")


Flow('lvkfilt','vkfilt','sfreal | sfmath output="log(input)"')
Flow('llbfilt','lvkfilt','math output="log(1+(%(bx)g*x1)^2+(%(by)g*x2)^2)"' % psyn)
Plot('llfilt','lvkfilt llbfilt',
     '''
     cmplx ${SOURCES[0:2]} |
     put label1='Ln(1+bx2kx2+by2ky2)' unit1= label2='Ln[F(kx,ky)]' unit2= |
     sfgraph title="Ln[F(kx,ky)] = p*Ln(1+bx2kx2+by2ky2)"
     ''',stdin=0)


# Display panel 1


Result('panel1','wgauss fwgauss vkfilt llfilt','TwoRows',vppen='xsize=10 ysize=10')


# Spectral filtering of white noise and inverse 2D Fourier transform


Flow('fcgauss',['vkfilt','fwgauss'],'sfmath r=${SOURCES[0]} p=${SOURCES[1]} type=complex output="r*p"')
ifft2 = 'sffft3 sym=y inv=y | sffft1 sym=y inv=y'
Flow('cgauss','fcgauss',ifft2)

Flow('rfcgauss','fcgauss','sfadd abs=y | sfreal')
Plot('cgauss',
     '''
     put label1=x unit1=m label2=y unit2=m | sfgrey color=j title="Correlated Gaussian noise"
     ''')
Plot('rfcgauss',
     '''
     put label1=kx unit1=1/m label2=ky unit2=1/m | sfgrey color=j title="Correlated Gaussian noise spectrum"
     ''')


# Logarithmic plot for separable least square inversion method ("Conical Coolie Hat")


Flow('lrfcgauss',['rfcgauss','mfwgauss'],'math r=${SOURCES[0]} m=${SOURCES[1]} output="log(r)-log(m)"')
Plot('llgauss','lrfcgauss llbfilt',
     '''
     cmplx ${SOURCES[0:2]} |
     put label1='Ln(1+bx2kx2+by2ky2)' unit1= label2='Ln[F(kx,ky)]' unit2= |
     sfgraph title="Ln[F(kx,ky)] = p*Ln(1+bx2kx2+by2ky2)"
     ''',stdin=0)


# Display panels 2 and 3


Result('panel2','rfcgauss llgauss vkfilt llfilt','TwoRows',vppen='xsize=10 ysize=10')
Result('panel3','wgauss fwgauss cgauss rfcgauss','TwoRows',vppen='xsize=10 ysize=10')


# -----------------------------------
# Inversion for stochastic parameters
# -----------------------------------

# Estimation of von Karman filter logarithm in spectral domain
# Optimization of (bx,bxy,by,nu) parameters by separate least square and Gauss Newton algorithm


Flow('erfcgauss',['rfcgauss','mfwgauss'],
     '''
     math r=${SOURCES[0]} m=${SOURCES[1]} output="r/m" |
     sfkarman2 verb=y niter=100 a0=100. b0=5. c0=1000.
     ''')


# Display panel 4


Flow('prfcgauss',['erfcgauss','mfwgauss'],'math r=${SOURCES[0]} m=${SOURCES[1]} output="r*m"')

Plot('prfcgauss',
     '''
     put label1=kx unit1=1/m label2=ky unit2=1/m |
     sfgrey color=j title="Filter estimated spectrum F(kx,ky)"
     ''')

Result('panel4','prfcgauss llgauss vkfilt llfilt','TwoRows',vppen='xsize=10 ysize=10')

# -----------
# Termination
# -----------


End()