Public documentattion

SeisAddNoise(d, snr; <keyword arguments>)

Add noise to an N-dimensional input data d. Should specify the signal-to-noise ratio level snr. Noise can be band limited using kewyord L.

Arguments

• d::Array{Real, N}: N-dimensional data.
• snr::Real: signal-to-noise ratio.
• db::Bool=false: Flag is false if snr is given by amplitude. Flag is true if

snr is given in dB.

• pdf::AbstractString="gaussian": random noise probability distribution:

"gaussian" or "uniform".

• L::Int=1: averaging operator length to band-limit the random noise.

Examples

julia> using PyPlot
julia> w = Ricker(); wn = SeisAddNoise(w, 2); plot(w); plot(wn);
MeasureSNR(w, wn)
julia> d = SeisHypEvents(); dn = SeisAddNoise(d, 1.0, db=true, L=9);
SeisPlotTX([d dn]); MeasureSNR(d, dn, db=true)

Credits: Juan I. Sabbione, 2016

SeisLinearEvents(; <keyword arguments>)

Generate up to five dimensional data consisting of linear events.

Arguments

• ot=0.0: first sample for the time axis in secs.
• dt=0.004: time sampling interval in secs.
• nt=500: number of time samples.
• ox1=0.0: first sample for the first spatial dimension in meters.
• dx1=10.0: sample interval for the first spatial dimension in meters.
• nx1=100: number of samples for the first spatial dimension.
• ox2=0.0: first sample for the second spatial dimension in meters.
• dx2=10.0: sample interval for the second spatial dimension in meters.
• nx2=1: number of samples for the second spatial dimension.
• ox3=0.0: second sample for the third spatial dimension in meters.
• dx3=10.0: sample interval for the third spatial dimension in meters.
• nx3=1: number of samples for the third spatial dimension.
• ox4=0.0: third sample for the fourth spatial dimension in meters.
• dx4=10.0: sample interval for the fourth spatial dimension in meters.
• nx4=1:number of samples for the fourth spatial dimension.
• tau=[1.0, 1.5]: intercept traveltimes for each event.
• p1=[0.0001,-0.0003]: Dip of events on the first dimension
• p2,p3,p4=[0, 0]: Dip of events on the following dimensions
• amp=[1.0,-1.0]: amplitudes for each linear event.
• f0=20.0: central frequency of wavelet for each linear event.

Example

julia> using SeisPlot
julia> d = SeisLinearEvents(); SeisPlotTX(d);

Credits: Aaron Stanton, 2015

SeisParabEvents(; <keyword arguments>)

Generate up to five dimensional data consisting of parabolic events.

Arguments

• ot=0.0: first sample for the time axis in secs.
• dt=0.004: time sampling interval in secs.
• nt=500: number of time samples.
• ox1=0.0: first sample for the first spatial dimension in meters.
• dx1=10.0: sample interval for the first spatial dimension in meters.
• nx1=100: number of samples for the first spatial dimension.
• ox2=0.0: first sample for the second spatial dimension in meters.
• dx2=10.0: sample interval for the second spatial dimension in meters.
• nx2=1: number of samples for the second spatial dimension.
• ox3=0.0: second sample for the third spatial dimension in meters.
• dx3=10.0: sample interval for the third spatial dimension in meters.
• nx3=1: number of samples for the third spatial dimension.
• ox4=0.0: third sample for the fourth spatial dimension in meters.
• dx4=10.0: sample interval for the fourth spatial dimension in meters.
• nx4=1:number of samples for the fourth spatial dimension.
• tau=[1.0, 1.6]: intercept traveltimes for each event.
• `p1=[0.0000,-0.0001]
• `p2=[0.0003, 0.0002]
• `p3=[-0.0001,-0.0001]
• `p4=[0.0001,-0.0000]
• amp=[1.0,-1.0]: amplitudes for each parabolic event.
• wavelet="ricker": wavelet used to model the parabolicr events.
• f0=[20.0]: central frequency of wavelet for each parabolic event.

Example

julia> using SeisPlot
julia> d = SeisParabEvents(); SeisPlotTX(d);

Credits: Mauricio D Sacchi, 2015

SeisHypEvents(; <keyword arguments>)

Generate two dimensional data consisting of hyperbolic events.

Arguments

• ot::Real=0.0: first sample for the time axis in secs.
• dt::Real=0.004: time sampling interval in secs.
• nt::Int=301: number of time samples.
• ox::Real=-1000.0: first sample for spatial dimension in meters.
• dx::Real=20.0: sample interval for the spatial dimension in meters.
• nx::Int=101: number of samples for the spatial dimension.
• tau::Vector{Real}=[0.2, 0.6, 0.9]: intercept traveltimes for each event.
• vel::Vector{Real}=[1500.0, 2000.0, 3000.0]: rms velocities in m/s
• apex::Vector{Real}=[0.0, 0.0, 0.0]: apex-shifts in meters.
• amp::Vector{Real}=[1.0, -1.0, 1.0]: amplitudes for each event.
• wavelet::AbstractString="ricker": wavelet used to model the events.
• f0::Vector{Real}=[20.0]: central frequency of wavelet for each event.

Output

• d::Array{Real, 2}: two dimensional data consisting of hyperbolic events.

Examples

julia> using SeisPlot
julia> d = SeisHypEvents(); SeisPlotTX(d);
julia> d = SeisHypEvents(apex=[100, 200, -300], f0=[30, 20, 15]);
SeisPlotTX(d);

SeisKolmogoroff(in)

Transform a wavelet into its minimum phase equivalent.

Arguments

• in::Array{Real,1}: input wavelet.

Example

julia> using PyPlot
julia> w = Ricker()
julia> wmin = SeisKolmogoroff(w)
julia> plot(w); plot(wmin)

Reference

• Claerbout, Jon F., 1976, Fundamentals of geophysical data processing.

McGraw-Hill Inc.

MeasureSNR(signal, noisy; db=false)

Measure the signal-to-noise ratio between the clean input signal and the contaminated input noisy.

Arguments

• signal::Array{Real, N}: N-dimensional clean signal. N must be <= 5.
• noisy::Array{Real, N}: N-dimensional noisy signal of same size as signal.
• db::Bool=false: Flag is false if the signal-to-noise ratio is measured by amplitude. Flag is true if snr is in dB.

Example

julia> d = SeisHypEvents(); dnoisy = SeisAddNoise(d, 2);
MeasureSNR(d, dnoisy)

SeisBandPass(in; <keyword arguments>)

Apply a bandpass filter to seismic data. Input and output data is 2D in tx domain. Filter is applied in fx domain.

Arguments

• in: Input 2D data array in tx domain. Time is first dimension.

Keyword arguments

• dt=0.004: time sampling interval in secs.
• fa=0: lower frequency cut [Hz]
• fb=0: lower frequency banpass [Hz]
• fc=60: upper frequency bandpass [Hz]
• fd=80: upper frequency cut [Hz]

Example

julia> d = SeisLinearEvents(); SeisPlotAmplitude(d,125,0.004);
julia> d_filter = SeisBandPass(d;dt=0.004,fa=2,fb=8,fc=12,fd=20); SeisPlotAmplitude(d_filter,125,0.004)
julia> SeisPlot([d d_filter],title="Data and Filtered data")

SeisBandPass(in,out,parameters; <keyword arguments>)

Arguments

• in::String: Input file - Seis format
• out::String: Output file - Seis format.
• parameters : list of the keyword arguments for the function SeisBandPass.

Keyword arguments

• group="gather" : Options are all, some or gather
• key=["imx","imy"] : Defines type of gather
• itrace=1 : Initial trace number
• ntrace=10000 : Total number of traces to process at once

Example

julia> param = Dict(:fa>=0,:fb=>0,:fc=>60,:fd=>80)
julia> SeisBandPass(filein,fileout, param,group="all")

SeisDecimate(in; <keyword arguments>)

Apply random or regular decimation to a multidimensional array of data. Input and output have the same dimensions.

Arguments

• in: input data as 2D, or 3D,4D,5D tensors. The first dimension is time.

Keyword arguments

• mode="random": decimation mode. Random is default, else decimates uniformly.
• perc=50: percentage of traces equal to zero (decimated). Only for random mode
• incx1=1: consecutive traces zeroed in first dimension. Only for regular decimation
• incx2=1: consecutive traces zeroed in second dimension.
• incx3=1: consecutive traces zeroed in third dimension.
• incx4=1: consecutive traces zeroed in fourth dimension.

Example

julia> d = SeisLinearEvents(); deci = SeisDecimate(d);

Credits: Aaron Stanton,2017

SeisDelay(in; <keyword arguments>)

Apply a time delay to 2D data

Arguments

• in: input 2D data array in tx domain. First dimension is time.

Keyword arguments

• delay=0.1: desired time delay in secs.
• dt=0.004: time sampling interval in secs.

Example

julia> d = SeisLinearEvents(); deli = SeisDelay(d);

Credits: Aaron Stanton,2017

SeisDiff(d ; <keyword arguments>)

Apply differentiation to seismic traces in freq. It can also be used to apply a phase rotation

Arguments

• d: Input 2D data array in tx domain. First dimension is time.

Keyword arguments

• delay=0.1: desired time delay in seconds.
• pow=-2: order of derivative
• rot=0: constant phase shift or rotation

Example

julia> d = SeisLinearEvents(); SeisPlotTX(d);
julia> d2 = SeisDiff(d;dt=0.004,pow=1); SeisPlotTX(d);

SeisEnvelope(in; <keyword arguments>)

Calculate the envelope attribute of a group of input traces.

Arguments

• in: input data. A 2D Array where the first dimension is time.

Example

julia> using PyPlot
julia> dtsec = 0.002; w = Ricker(dt=dtsec);
julia> e = SeisEnvelope(w); t=dtsec*collect(0:1:length(w)-1);plot(t,w,t,e,"-r");xlabel("Time [s]")

Credits: Aaron Stanton,2017

SeisFKFilter(in; <keyword arguments>)

Removes energy from seismic data by applying filters to a gather in the Frequency-wavenumber domain

Arguments

• in: 2D input data in TX domain. First dimension is time.

Keyword arguments

• dt=0.004: time sampling interval in secs.
• dx=10: space sampling interval in meters.
• va=-2000,vb=-3000,vc=3000,vd=2000: apparent velocity corners to filter

Output

• out: Filtered data in TX domain

Example

julia> d = SeisLinearEvents(); df = SeisFKFilter(d,va=-4000.,vb=-6000.,vc=6000.,vd=4000.);

Credits: Aaron Stanton,2017

SeisGain(d ; <keyword arguments>)

Gain a group of traces. Input and output are 2D.

Arguments

• d::Array{Real,2}: two dimensional data.

Keyword arguments

• dt::Real=0.004: sampling interval in secs.
• kind::AbstractString="time": if kind="time", gain = t.^a . * exp(-bt); if kind="agc", automatic gain control is applied.
• coef::Vector{Real}=[2.0,0.0]: if kind="time", coef = [a,b]; if kind="agc", coef = [agc_gate]
• normal::Int=0: normal=0 no normalization; normal=1 normalize each trace by amplitude; normal=2 normalize each trace by rms value/

Example

julia> using PyPlot
julia> d = SeisHypEvents();
       dout = SeisGain(d, kind="agc", coef=[0.05]);
       SeisPlotTX([d dout]);

Credits: Aaron Staton Updates: Juan I. Sabbione, Mauricio D. Sacchi, Fernanda Carozzi

SeisGain(in,out,parameters; <keyword arguments>)

Gain a group of traces. Input and output are file names.

Arguments

• in::String: Input file - Seis format
• out::String: Output file - Seis format.
• parameters : list of the keyword arguments for the function SeisGain.

Keyword arguments

• group="gather" : Options are all, some or gather
• key=["imx","imy"] : Defines type of gather
• itrace=1 : Initial trace number
• ntrace=10000 : Total number of traces to process at once

Example

julia> param = Dict(:dt=>0.01,:kind=>"time",:coef=>[2.0,0.0],:normal=>0)
julia> SeisGain(filein,fileout, param,group="all")

SeisMute(d; <keyword arguments>) -> Array{Real,2}

Mute trace noise before the first arrival. The muting function depends on the offset between source and receiver and the velocity of the first layer as tmute = t0 + offset/v1.

Arguments

• d: Array of input traces.

Keyword arguments

• offset=[0.]: Vector of distances between source and receiver.
• tmute=0.: initial muting time
• vmute=1500.: Velocity of the first layer.
• taper=0.1: Taper of the muting function
• dt::Real=0.004: sampling interval in secs.

SeisMute(in,out,parameters; <keyword arguments>)

Mutes noise before first arrival of traces in a file. Saves the muted traces to an output file. The muting function depends on the offset between source and receiver and the velocity of the first layer as tmute = t0 + offset/v1.

Arguments

• in::String: Input file - Seis format.
• out::String: Output file - Seis format.
• parameters : list of the keyword arguments for the function SeisMute.

Keyword arguments

• group="gather" : Options are all, some or gather
• key=["imx","imy"] : Defines type of gather
• itrace=1 : Initial trace number
• ntrace=10000 : Total number of traces to process at once

Example

julia> param = Dict(:tmute=>0.0, :vmute=>10000, :taper=>0.05,:dt=>0.01)
julia> SeisMute(filein,fileout, param,group="some")

SeisNMO(in,out,parameters; <keyword arguments>)

Arguments

• in::String: Input file - Seis format.
• out::String: Output file - Seis format.
• parameters : list of the keyword arguments for the function SeisMute.

Keyword arguments

• group="gather" : Options are all, some or gather
• key=["imx","imy"] : Defines type of gather
• itrace=1 : Initial trace number
• ntrace=10000 : Total number of traces to process at once

Example

julia> param = Dict(:tmute=>0.0, :vmute=>10000, :taper=>0.05,:dt=>0.01)
julia> SeisMute(filein,fileout, param,group="some")

SeisPWD(in;<keyword arguments>)

Dip estimation by Plane Wave Destruction. see http://sepwww.stanford.edu/data/media/public/sep//prof/pvi.pdf Chapter 4.

Arguments

-w1=10 -w2=10 -dz_in=1 -dx_in=1 -format="angle"

  SeisProcessFile(in,out,operators,parameters;<keyword arguments>)

Run processing flows that read and write from disk

f is a function that has the following syntax: d2 = f(d1,param), where param is list of keyword arguments for the function. Note that f can be a vector of functions. They will be executed sequentially on the same group of traces.

Arguments

• in: input filename of type String or Array{String,1} - Seis format
• out: output filenames of type String or Array{String,1} - Seis format
• operators
• parameters

Keyword arguments

• group="gather" : Options are all, some or gather
• key=["imx","imy"] : Defines type of gather
• itrace=1 : Initial traces
• ntrace=10000 : Total number of traces to process at once

Example

Apply a bandpass filter to a seismic cube sequentially, by shot gather. Assume dt is equal to 0.002.

julia> operators = [SeisBandPass]
julia> param = [Dict(:dt=>0.002, :fa=>20,:fb=>30,:fc=>80,:fd=>90)]
julia> SeisProcessFile(filein,fileout,operators,param,key=["sx"])

  SeisProcessFile(in1,in2,out,operators,parameters;<keyword arguments>)

Run processing flows that read from 2 files and write to disk

Arguments

• in1::String: input filename - Seis format.
• in2::String: input filename - Seis format.
• out::String: output filename - Seis format.
• operators
• parameters

Keyword arguments

• group="gather" : Options are all, some or gather
• key=["imx","imy"] : Defines type of gather
• itrace=1 : Initial traces
• ntrace=10000 : Total number of traces to process at once

Example

Apply a bandpass filter to a seismic cube sequentially, by shot gather. Assume dt is equal to 0.002.

julia> operators = [SeisBandPass]
julia> param = [Dict(:dt=>0.002, :fa=>20,:fb=>30,:fc=>80,:fd=>90)]
julia> SeisProcessFile(filein,fileout,operators,param,key=["sx"])

SeisRadonFreqFor(in, nt; <keyword arguments>)

Transform a tau-p gather to a time-offset gather using a frequency domain forward parabolic or linear Radon operator.

Arguments

• in::Array{Float64,2}: 2D Radon panel, in[1:ntau,1:np], where ntau is the

number of intercept times and np the number of curvatures or ray parameters.

• nt::Int: number of time samples.

Keyword arguments

• order="parab": "parab" for parabolic transform, "linear"

for linear transform.

• dt=0.004: time sampling interval in seconds.
• h=collect(0.0:20.0:1000.0): offset vector; h[1:nh].
• href=0.0: reference offset for parabolic Radon Transform. If the

defautl value href=0.0 is given, href is set to max(abs(h)).

• p=collect(-0.05:0.01:2.2): p[1:np]. If order="parab", p

is a vector of residual moveout ("curvatures") at reference offset href in seconds; if order=linear, p is a vector of ray parameters in s/m.

• flow=0.0: minimum frequency in the data in Hz.
• fhigh=125.0: maximum frequency in the data in Hz.

Output

• d: data synthetized via forward Radon modeling, d[1:nt, 1:nh].

References

• Hampson, D., 1986, Inverse velocity stacking for multiple elimination:

Canadian Journal of Exploration Geophysics, 22, 44-55.

• Sacchi, M.D. and Ulrych, T.J., 1995, High-resolution velocity gathers and

offset space reconstruction: Geophysics, 60, 1169-1177.

SeisRadonFreqInv(d; <keyword arguments>)

Transform a Gather from time-offset gather to tau-p gather using a frequency domain inverse parabolic or linear Radon operator via least-squares inversion.

Arguments

• d::Array{Float64,2}: 2D data, d[1:nt,1:nh], where nt is number of

time samples and nh the number of receivers.

Keyword arguments

• order::"parab": "parab" for parabolic transform, "linear"

for linear transform.

• dt=0.004: sampling interval in seconds.
• h=collect(0.0:20.0:1000.0): offset vector h[1:nh].
• href=0.0: reference offset for parabolic Radon Transform. If the

defautl value href=0.0 is given, href is set to max(abs(h)).

• p=collect(-0.05:0.01:2.2): p[1:np]. If order="parab", p

is a vector of residual moveout ("curvatures") at reference offset href in seconds; if order=linear, p is a vector of ray parameters in s/m.

• flow=0.0: minimum frequency in the data in Hz.
• fhigh=125.0: maximum frequency in the data in Hz.
• mu=0.001: trade off parameter or damping for the L.S. inversion.

Output

• m: inverted Radon panel m[1:ntau, 1:np].

References

• Hampson, D., 1986, Inverse velocity stacking for multiple elimination:

Canadian Journal of Exploration Geophysics, 22, 44-55.

• Sacchi, M.D. and Ulrych, T.J., 1995, High-resolution velocity gathers and

offset space reconstruction: Geophysics, 60, 1169-1177.

SeisSincInterp1D(d,order)

Resample seismic traces via 1D sinc interpolation. The time series are sampled every dt secs, the output corresponds to series with time interval dt/order.

Arguments

• d::Array{Real,N}: N-dimensional data, first dimension is time
• order::Integer: order of interpolation 2,4

Examples

2-time upsampling of a wavelet

julia> order = 2; dt=0.004; w = Ricker(dt=dt, f0=20); t = collect(0:1:length(w)-1)*dt;
julia> wout = SeisSincInterp1D(w,order); tout = collect(0:1:length(wout)-1)*dt/order;
julia> plot(t,w); plot(tout,wout,"*")

4-time upsampling of a gather

julia> d = SeisLinearEvents(); di = SeisSincInterp1D(d,4);
julia>  SeisPlotTX(d,style="wiggles")
julia>  SeisPlotTX(di,style="wiggles")  

MAY 2018, MDS