up [pdf]
# Nonlinear estimation of von Karman autocorrelation in 1D Fourier domain for well data
#
# July 2008
#
# Thomas Jules Browaeys 
# Bureau of Economic Geology
# University of Texas at Austin
# mailto:jules.browaeys@beg.utexas.edu


from rsf.proj import *
from rsf.recipes.beg import server as private
from math import pi

# ----------------------------------------------------------------------
# WELL LOG DATA RHO,VP,VS
# ----------------------------------------------------------------------

# WELLS A & B
# DATA (07): z(col=0), rho(col=5), VP(col=6), VS(col=7)

# WELLS C & D
# DATA (04): z(col=0), rho(col=1), VP(col=2), VS(col=3)

# Data files
dfiles = {'A':'07-23-36-10.txt',
          'B':'07-26-36-10.txt',
          'C':'04-02-37-10.txt',
          'D':'04-13-36-10.txt'
          }

Fetch(dfiles.values(),'apache',private)

# Grids 1-D and data type
pgrids = {'A':{'nx':7033, 'ny':8, 'ox':155.875, 'dx':0.125},
          'B':{'nx':7044, 'ny':8, 'ox':157.125, 'dx':0.125},
          'C':{'nx':6921, 'ny':4, 'ox':142.5,   'dx':0.125},
          'D':{'nx':7113, 'ny':4, 'ox':161.75,  'dx':0.125}
          }

# Frequency restriction and plot parameters
fcbpms = {
    # Well A - VP
    'A1':{'bp':0.7847,   'fc1':0.6784,      'fc2':1.4850278,
          'bp1':6.698,   'bp2':0.349,       'bp3':1.,
          'col':6,
          'optsig':'title="Sonic log \F5 V\F2 \_P\^ \k100 \k100 Well 1" format2=%3.1f',
          'optrfs':'title="Sonic log V\_P\^ spectrum - Well 1" format2=%3.1f',
          'optest':'title="Estimation on sonic log \F5 V\F2 \_P\^ \k100 \k100 Well 1"',
          'optesr':'title="Restriction for sonic log \F5 V\F2 \_P\^ \k100 \k100 Well 1"',
          'min2sc':-0.9, 'max2sc':0.9},
    # Well A - VS
    'A2':{'bp':1.0466,   'fc1':0.5069,      'fc2':1.430250,
          'bp1':5.9241,  'bp2':0.2833,      'bp3':1.,
          'col':7,
          'optsig':'title="Sonic log V\_S\^ - Well 1" format2=%3.1f',
          'optrfs':'title="Sonic log V\_S\^ spectrum - Well 1" format2=%3.1f',
          'optest':'title="Estimation - V\_S\^ - Well 1"',
          'optesr':'title="Restriction - V\_S\^ - Well 1"',
          'min2sc':-1.1, 'max2sc':1.1},
    # Well B - VP
    'B1':{'bp':1.89926,  'fc1':0.6133,      'fc2':1.681,
          'bp1':5.34,    'bp2':1.0,         'bp3':1.,
          'col':6,
          'optsig':'title="Sonic log V\_P\^ - Well 2" format2=%3.1f',
          'optrfs':'title="Sonic log V\_P\^ spectrum - Well 2" format2=%3.1f',
          'optest':'title="Estimation - V\_P\^ - Well 2"',
          'optesr':'title="Restriction - V\_P\^ - Well 2"',
          'min2sc':-1.1, 'max2sc':1.1},
    # Well B - VS
    'B2':{'bp':2.84062,  'fc1':0.68047651,  'fc2':1.518522,
          'bp1':3.0759,  'bp2':4.5,         'bp3':1.,
          'col':7,
          'optsig':'title="Sonic log V\_S\^ - Well 2" format2=%3.1f',
          'optrfs':'title="Sonic log V\_S\^ spectrum - Well 2" format2=%3.1f',
          'optest':'title="Estimation - V\_S\^ - Well 2"',
          'optesr':'title="Restriction - V\_S\^ - Well 2"',
          'min2sc':-1.4, 'max2sc':1.4},
    # Well C - VP
    'C1':{'bp':1.3377,   'fc1':0.5198378,   'fc2':1.854343528,
          'bp1':7.21752, 'bp2':0.3,         'bp3':1.,
          'col':2,
          'optsig':'title="Sonic log V\_P\^ - Well 3" format2=%3.1f',
          'optrfs':'title="Sonic log V\_P\^ spectrum - Well 3" format2=%3.1f',
          'optest':'title="Estimation - V\_P\^ - Well 3"',
          'optesr':'title="Restriction - V\_P\^ - Well 3"',
          'min2sc':-0.9, 'max2sc':0.9},
    # Well C - VS
    'C2':{'bp':1.246,    'fc1':0.558081343, 'fc2':1.99418703,
          'bp1':5.01091, 'bp2':0.352044,    'bp3':1.,
          'col':3,
          'optsig':'title="Sonic log \F5 V\F2 \_S\^ \k100 \k100 Well 3" format2=%3.1f  n2tic=5 o2num=-1.0 d2num=0.5',
          'optrfs':'title="Sonic log V\_S\^ spectrum - Well 3" format2=%3.1f',
          'optest':'title="Estimation on sonic log \F5 V\F2 \_S\^ \k100 \k100 Well 3"',
          'optesr':'title="Restriction for sonic log \F5 V\F2 \_S\^ \k100 \k100 Well 3"',
          'min2sc':-1.4, 'max2sc':1.4},
    # Well D - VP
    'D1':{'bp':0.638323, 'fc1':0.706749585, 'fc2':1.82077179,
          'bp1':2.5841,  'bp2':0.19704,     'bp3':1.,
          'col':2,
          'optsig':'title="Sonic log V\_P\^ - Well 4" format2=%3.1f',
          'optrfs':'title="Sonic log V\_P\^ spectrum - Well 4" format2=%3.1f',
          'optest':'title="Estimation - V\_P\^ - Well 4"',
          'optesr':'title="Restriction - V\_P\^ - Well 4"',
          'min2sc':-0.9, 'max2sc':0.9},
    # Well D - VS
    'D2':{'bp':0.57290,  'fc1':0.706749585, 'fc2':1.82077179,
          'bp1':2.5841,  'bp2':0.19704,     'bp3':1.,
          'col':3,
          'optsig':'title="Sonic log V\_S\^ - Well 4" format2=%3.1f',
          'optrfs':'title="Sonic log V\_S\^ spectrum - Well 4" format2=%3.1f',
          'optest':'title="Estimation - V\_S\^ - Well 4"',
          'optesr':'title="Restriction - V\_S\^ - Well 4"',
          'min2sc':-0.9, 'max2sc':0.9}
          }


for soniclog in fcbpms.keys():

    well = soniclog[0]
    tsonic = soniclog[1]

    dfile = dfiles[well]
    pgrid = pgrids[well]

    fcbpm = fcbpms[soniclog]

    pgrid['col'] = fcbpm['col']

    # Files
    sonic    = 'sonic'    + soniclog
    mscale   = 'mscale'   + soniclog
    msonic   = 'msonic'   + soniclog
    signal   = 'signal'   + soniclog
    fsignal  = 'fsignal'  + soniclog
    rfsignal = 'rfsignal' + soniclog

    # Sonic log signal
    Flow(sonic,dfile,
         '''
         echo n1=%(ny)d n2=%(nx)d o2=%(ox)g d2=%(dx)g in=$SOURCE data_format=ascii_float |
         dd form=native | transp |
         window n1=%(nx)d n2=1 f2=%(col)d |
         put label1=Depth unit1=m label2=Velocity unit2=m/s
         ''' % pgrid)

    # Plot
    Result(sonic,'graph title="Sonic well log"')

    # Shift and scale by mean signal
    Flow(mscale,sonic,'stack axis=1 norm=n | math output="input*%g"' % (1./pgrid['nx']))
    Flow(msonic,mscale,'spray axis=1 n=%d o=%g' % (pgrid['nx'],pgrid['ox']))
    Flow(signal,[sonic,msonic],'math r=${SOURCES[0]} m=${SOURCES[1]} output="r/m-1."')

    # Plot scaled signal
    Result(signal,
           '''
           graph min1=%g max1=%g min2=%g max2=%g label2="V/V\_0\^\k60 \F15 8\F2 1" unit2= %s parallel2=n font=2
           ''' % (pgrid['ox'],pgrid['ox']+(pgrid['nx']-1.)*pgrid['dx'],fcbpm['min2sc'],fcbpm['max2sc'],fcbpm['optsig'] ))

    # Discrete Fourier Transform (fft) (k=0,N-1)
    Flow(fsignal,signal,'fft1 sym=y')
    Flow(rfsignal,fsignal,'add abs=y | real | put label1=f unit1=1/m label2= unit2=')

    # Files
    irfsignal = 'irfsignal' + soniclog
    fitfilt   = 'fitfilt'   + soniclog
    illbfilt  = 'illbfilt'  + soniclog
    lifsignal = 'lifsignal' + soniclog
    lfsignal  = 'lfsignal'  + soniclog
    illsignal = 'illsignal' + soniclog
    llsignal  = 'llsignal'  + soniclog

    llisignalf = 'llisignalf' + soniclog


    # ------------------------------------------
    # INVERSION FOR VON KARMAN FILTER PARAMETERS 
    # ------------------------------------------
    # Estimation of von Karman filter logarithm in spectral domain
    # Separable least squares for exponent and amplitude
    # Newton algorithm on nonlinear parameter b*b

    Flow(irfsignal,rfsignal,'karman verb=y niter=100 x0=.2')

    # Plot spectrum and estimation
    Result(fitfilt,[rfsignal, irfsignal],
           '''
           cat ${SOURCES[0:2]} axis=2 |
           put label2= unit2= label1=f unit1=1/m |
           graph min1=0. min2=0. max1=1. max2=2.5 %s parallel2=n plotfat=1,5
           ''' % fcbpm['optrfs'],stdin=0)

    # Logarithmic plots
    Flow(illbfilt,irfsignal,'math output="log(1.+((%g)*x1)^2)"' % (2.*pi*fcbpm['bp']))
    Flow(lifsignal,irfsignal,'math output="log(input)"')
    Flow(lfsignal,rfsignal,'math output="log(input)"')
    Flow(illsignal,[illbfilt, lifsignal],'cmplx ${SOURCES[0:2]}',stdin=0)
    Flow(llsignal,[illbfilt, lfsignal],'cmplx ${SOURCES[0:2]}',stdin=0)

    # Plot nonlinear spectral estimation 
    Result(llisignalf,[llsignal, illsignal],
           '''
           cat ${SOURCES[0:2]} axis=2 |
           graph label1="Ln[1+(k\F5 b\F2 )\^2\_ ]" unit1= label2="Ln[F(k)]" unit2= 
           min1=0. max1=6.0 min2=-9. max2=1. %s parallel2=n plotfat=1,5 font=2
           ''' % fcbpm['optest'],stdin=0)


    # ---------------------------------------------
    # DOMAIN RESTRICTIONS FOR PARAMETERS ESTIMATION
    # ---------------------------------------------

    # Clipping frequencies are f = sqrt(exp( )-1)/(2.*pi*b)
    # Close initial guess -> converge

    for restricd in ('L', 'M', 'H'):

        rrfsignal  = 'rrfsignal'  + soniclog + restricd
        irrfsignal = 'irrfsignal' + soniclog + restricd
        rfitfilt   = 'rfitfilt'   + soniclog + restricd

        irllbfilt  = 'irllbfilt'  + soniclog + restricd
        lirfsignal = 'lirfsignal' + soniclog + restricd
        lrfsignal  = 'lrfsignal'  + soniclog + restricd
        irllsignal = 'irllsignal' + soniclog + restricd
        rllsignal  = 'rllsignal'  + soniclog + restricd

        rllisignalf = 'rllisignalf' + soniclog + restricd

        if restricd == 'L':

            # Restriction - Low frequency
            Flow(rrfsignal,rfsignal,'window max1=%g' % (fcbpm['fc1']))

        elif restricd == 'M':

            # Restriction - Medium frequency
            Flow(rrfsignal,rfsignal,'window min1=%g max1=%g' % (fcbpm['fc1'],fcbpm['fc2']))

        elif restricd == 'H':

            # Restriction - High frequency
            Flow(rrfsignal,rfsignal,'window min1=%g' % (fcbpm['fc2']))


        # Parameters spectral estimation 
        Flow(irrfsignal,rrfsignal,'karman verb=y niter=100 x0=.1')

        # Plot spectrum and estimation 
        Result(rfitfilt,[rrfsignal, irrfsignal],
               '''
               cat ${SOURCES[0:2]} axis=2 |
               put label2= unit2= label1=f unit1=1/m |
               graph min1=0. min2=0. max1=1. max2=2.5
               title="Spectrum F(k) restricted" parallel2=n plotfat=1,5           
               ''',stdin=0)

        # Logarithmic plots
        if restricd == 'L':
            Flow(irllbfilt,irrfsignal,'math output="log(1.+((%g)*x1)^2)"' % (2.*pi*fcbpm['bp1']))
        elif restricd == 'M':
            Flow(irllbfilt,irrfsignal,'math output="log(1.+((%g)*x1)^2)"' % (2.*pi*fcbpm['bp2']))
        elif restricd == 'H':
            Flow(irllbfilt,irrfsignal,'math output="log(1.+((%g)*x1)^2)"' % (2.*pi*fcbpm['bp3']))

        Flow(lirfsignal,irrfsignal,'math output="log(input)"')
        Flow(lrfsignal,rrfsignal,'math output="log(input)"')
        Flow(irllsignal,[irllbfilt, lirfsignal],'cmplx ${SOURCES[0:2]}',stdin=0)
        Flow(rllsignal,[irllbfilt, lrfsignal],'cmplx ${SOURCES[0:2]}',stdin=0)

        # Plot nonlinear spectral estimation
        Result(rllisignalf,[rllsignal, irllsignal],
               '''
               cat ${SOURCES[0:2]} axis=2 |
               graph label1="Ln[1+(k\F5 b\F2 )\^2\_ ]" unit1= label2="Ln[F(k)]" unit2= 
               min1=0. max1=5.5 min2=-9. max2=1. %s parallel2=n plotfat=1,5 font=2
               ''' % fcbpm['optesr'],stdin=0)

End()

sfdd
sftransp
sfwindow
sfput
sfgraph
sfstack
sfmath
sfspray
sffft1
sfadd
sfreal
sfkarman
sfcat
sfcmplx