SeismicQ.jl

Ricker(t, t₀, 𝑓₀)

Compute the Ricker function for t, t₀ and 𝑓₀.

\[f = (1.0 - 2.0 (\pi 𝑓₀ (t - t₀))^2) \exp(-π^2 𝑓₀^2 (t - t₀)^2)\]

Examples

julia>  Ricker(0.0, 0.0, 0.0)
1.0

Ricker(x, x₀, t, t₀, 𝑓₀)

Compute the Ricker function with a dirac.

\[f = \delta(x₀) \exp{(-((x-x₀)^2)/2.0/σ₀^2)} (1.0 - 2.0 (\pi 𝑓₀ (t - t₀))^2) \exp(-π^2.0𝑓₀^2.0(t - t₀)^2)\]

Examples

julia> Ricker(2., 2., 0.0, 0.0, 0.0)
0.1353352832366127

Ricker(x, x₀, t, t₀, 𝑓₀, σ₀)

Compute the Ricker function with 1D spatial support.

\[f = \exp{(-((x-x₀)^2)/2.0/σ₀^2)} (1.0 - 2.0 (\pi 𝑓₀ (t - t₀))^2) \exp(-π^2.0𝑓₀^2.0(t - t₀)^2)\]

Examples

julia> Ricker(2., 0., 0.0, 0.0, 0.0, 1)
0.1353352832366127

Ricker(x, x₀, y, y₀, t, t₀, 𝑓₀, σ₀)

Compute the Ricker function with 2D spatial support.

\[f = \exp{(-((x-x₀)^2)/2.0/σ₀^2)} (1.0 - 2.0 (\pi 𝑓₀ (t - t₀))^2) \exp(-π^2.0𝑓₀^2.0(t - t₀)^2)\]

Examples

julia> Ricker(2., 0., 2., 0., 0.0, 0.0, 0.0, 1)
0.01831563888873418

θs(G,Δt)

Deviatoric part of the elastic constitutive update: provides the shear modulus G for use in

\[Δτ = 2G Δε\]

θs(G,η,Δt)

Deviatoric part of the Maxwell visco-elastic constitutive update: used to compute

\[Δτ = 2 G (DeN / (1+DeN)) Δε + Δτ_{RELAX}\]

with η the shear viscosity, and DeN a numerical Deborah number, equal to

\[DeN = η / (G Δt)\]

χs(G,Δt)

Relaxation part of the deviatoric elastic constitutive update: by definition equal to 1 (no relaxation)

χs(G,η,Δt)

Relaxation part of the deviatoric Maxwell visco-elastic constitutive update: for use in the calculation of

\[Δτ = (DeN / (1+DeN)) τ_\mathrm{OLD} + Δτ_\mathrm{CONSTITUTIVE}\]

with DeN a numerical Deborah number, equal to

\[DeN = η / (G Δt)\]

and η the shear viscosity.

θb(K,Δt)

Volumetric part of the elastic constitutive update: provides the bulk modulus K for use in

\[Δp = -K ∇⋅v Δt\]

θb(K,ηb,Δt)

compute effective bulk viscosity for a Kelvin visco-elastic constitutive update

χb(K,Δt)

compute effective bulk viscosity linear elastic update: 0.

χb(K,ηb,Δt)

compute effective bulk viscosity for a Kelvin visco-elastic constitutive update: here obviously it is for a linear update

Spec(trace,t0,dt,Nt,tp,ts,\Deltaph)

Returns the amplitude spectrums of the P ans S phases where amplitudes are significant

    trace: the trace recorded at a given virtual station, starting at -t0, with time sampling dt, Nt samples
    tp: picked P-wave phase (s)
    ts: picked S-wave phase (s)
    Δph: width of the phases (s) 

    Vec2dP: modulus of the fft coeffs of P phase where they are larger than max(Vec2dp)/10 
    Vec2dS: modulus of the fft coeffs of S phase where they are larger than max(Vec2dp)/10
    FrequP: corresponding frequency vector to plot P-wave amp. spectrum
    FrequS: corresponding frequency vector to plot S-wave amp. spectrum

Examples

julia>  (AmpP,AmpS,fP,fS)=Spec(trace1,0.1,1e-3,2000,0.7,1.25,0.35))
Trace length is 2.0 s using 2000 samples
Trace 1 at offset 3000 m is found at index 30
Trace 2 at offset 5000 m is found at index 50
Trace 1:
phase P is picked at 0.4 s
phase S is picked at 0.75 s
phase width is set to 0.35
Trace 2:
phase P is picked at 0.7 s
phase S is picked at 1.25 s
phase width is set to 0.35
+ Figure with 8 subplots

(x,y)=ComputeQgraph(amp1,amp2,freq,index,d1,d2,V)

Returns x and y to plot Q-graph for a given phase using the results of Spec() for two different traces from two stations

    y=V/(π(d2-d1))*ln(amp2(f)/amp1(f)
    x=f

    Considering two attenuated waves at two different distances from the source (d1<d2), we obtain the following relation:
    V/(π(d2-d1))*ln(amp2(f)/amp1(f)=-f/Q+V*n/(\pi(d2-d1))*ln(d1/d2)

    After plotting (plot(x,y)), a flat curve indicates infinite Q (no attenuation). This is what happens with the "fake" data used
    for the developpment of these routines. In an homogeneous Q medium the graph should be linear. In more complex attenuation media,
    the graph could be non-linear?

    Velocity is considered constant: the current code does not account for possible dispersion...

Examples

julia>  

Getfreq(dt,Nt)

Returns the frequency vector corresponding to the fft results

    dt is the time sample in s
    Nt is the number of samples in the trace given to fft()
    The frequency increment is 1/(Nt*dt)
    Vector goes beyond Nyquist Frequency as fft also does. Only the first half of the vectors are useful for us

Examples

julia>  Getfreq(1e-3,2000)
2000-element Vector{Float64}:
    0.5
    1.0
    1.5
    2.0
    2.5
    3.0
    3.5
    4.0
    ⋮
  996.5
  997.0
  997.5
  998.0
  998.5
  999.0
  999.5
 1000.0

time_vec,acc_vec = GenAttenuatedRicker(listₓ,Δt,Nt,Vp,Vs,αp,αs,𝑓₀)  ;

Function that generates a seismic trace (Ricker wavelet) of P and S waves. This wave is attenuated along time and the function creates a matrix of wave amplitude as a function of time and distance to the source. The function takes as inputs: listₓ : vector of geophone distances to the source [m] Δt : time step of the wave signal [s] Nt : number of time steps of the wave signal Vp : P-wave velocity [m/s] Vs : S-wave velocity [m/s] αp : attenuation factor for P-wave αs : attenuation factor for S-wave 𝑓₀ : central frequency of the source [Hz]

and return: timevec : vector containing the time steps of the received wave signal accvec : matrix of wave acceleration at each geophone position

Examples

julia>  time_vec,acc_vec = GenAttenuatedRicker(0:1000:5000,1e-3,2000,7000,4000,2e-4,4e-4,10.0) 

PlotReceiverGather(listₓ,time_vec,acc_vec)

Function that plots a receiver gather The function takes as inputs: listₓ : vector of geophone distances to the source [m] timevec : vector containing the time steps of the wave signal [s] accvec : matrix of wave acceleration at each geophone position

Examples

julia>  PlotReceiverGather(0:1000:5000,time_vec,acc_vec)