sconstruct

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

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


from rsfproj import *


#----------------
# Spatial 1D grid
#----------------
#
# Master parameters
#
# nx = number of points in X
# dx = 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


pgrid = {'nx':512, 'ox':0., 'dx':1.}

Flow('spacegrid',None,'spike nsp=1 mag=1 n1=%(nx)d d1=%(dx)g o1=%(ox)g | put label1=x unit1=m' % pgrid)


# Synthetic 1D Gaussian white noise
#
# b   = correlation length scale in x (m)
# nu  = Hurst exponent -0.5 < and < 0.8
# gmu = Gaussian white noise mean
# gvr = Gaussian white noise variance = (standard deviation)^2
# grn = noise range
# gsd = seed for random generator


psyn = {'b':3., 'nu':0.3, 'gmu':0., 'gvr':1., 'grn':1., 'gsd':1213}

Flow('wgauss','spacegrid','noise mean=%(gmu)g range=%(grn)g rep=y seed=%(gsd)g type=y var=%(gvr)g' % psyn)


# 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) = S*(u)
# Periodicity S(ix+NX) = S(ix)
# 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] is nx/2+1
# [fx] = dfx*(-nx/2+1:1:nx/2)


Flow('fwgauss','wgauss','fft1 sym=y')

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

Flow('vkfilt','fwgauss','math output="(1+(%(b)g*x1)^2)^(-0.25-0.5*%(nu)g)"' % psyn)
Flow('rvkfilt','vkfilt','add abs=y | real | put label1=kx unit1=1/m')

# Filtering and inverse Fourier transform

Flow('fcgauss',['vkfilt','fwgauss'],'math r=${SOURCES[0]} p=${SOURCES[1]} type=complex output="r*p"')
Flow('rfcgauss','fcgauss','add abs=y | real | put label1=kx unit1=1/m')
Flow('cgauss','fcgauss','fft1 sym=y inv=y')

# Plots

Plot('wgauss','graph title="Gaussian white noise"')
Plot('fwgauss','add abs=y | real | put label1=kx unit1=1/m | graph title="Gaussian white noise spectrum"')
Plot('rfcgauss','graph title="Correlated Gaussian spectrum"')
Plot('rvkfilt','graph title="Filter VK spectrum"')
Plot('cgauss','put label1=x unit1=m | graph title="Correlated Gaussian noise"')
Plot('vkcfilt','rvkfilt rfcgauss',
     '''
     cat ${SOURCES[0:2]} axis=2 | graph title="Correlated Gaussian spectrum"
     ''',stdin=0)

Result('panel1','wgauss fwgauss vkcfilt cgauss','TwoRows',vppen='xsize=10 ysize=10')


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


Flow('llbfilt','rvkfilt','math output="log(1+((%(b)g)*x1)^2)"' % psyn)
Flow('lvkfilt','rvkfilt','math output="log(input)"')
Flow('lfcgauss','rfcgauss','math output="log(input)"')
Plot('llfilt','lvkfilt llbfilt',
     '''
     cmplx ${SOURCES[0:2]} |
     put label1='Ln(1+b2kx2)' unit1= label2='Ln[F(kx)]' | graph title="Ln[F(kx)] = p*Ln(1+b2kx2)"
     ''',stdin=0)
Plot('llgauss','lfcgauss llbfilt',
     '''
     cmplx ${SOURCES[0:2]} |
     put label1='Ln(1+b2kx2)' unit1= label2='Ln[F(kx)]' | graph title="Ln[F(kx)] = p*Ln(1+b2kx2)"
     ''',stdin=0)

Result('panel2','rvkfilt rfcgauss llfilt llgauss','TwoRows',vppen='xsize=10 ysize=10')


# -----------------------------------
# Inversion for stochastic parameters
# -----------------------------------
#
# Estimation of von Karman filter logarithm in spectral domain
# Separate least square for exponent and amplitude
# Newton algorithm on nonlinear parameter b*b

# Scaling and initial guess

# Iterative inversion

Flow('irfcgauss','rfcgauss','karman verb=y niter=100')
Plot('fitfilt','irfcgauss rfcgauss',
     '''
     cat ${SOURCES[0:2]} axis=2 |
     put label1=kx unit1=1/m | graph title="Filter estimated spectrum F(kx)"
     ''',stdin=0)

Result('panel4','fitfilt llgauss rvkfilt llfilt','TwoRows',vppen='xsize=10 ysize=10')


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


End()