From 1c94e180d806c6b4f8957a577023cc05292d16b2 Mon Sep 17 00:00:00 2001 From: MDelescluse <143093360+MDelescluse@users.noreply.github.com> Date: Thu, 24 Aug 2023 16:31:32 +0200 Subject: [PATCH] ajout de l'estimation de Q (#2) * ajout de l'estimation de Q * Camel cases modifs * with some documentation --- Project.toml | 4 + docs/src/index.md | 3 + examples/FausseTrace.jl | 89 ++++++++++++ examples/GetQTest.jl | 86 +++++++++++ examples/{Source_Day1.jl => SourceDay1.jl} | 0 src/GenFakeData.jl | 64 ++++++++ src/SeismicQ.jl | 6 + src/Spec.jl | 161 +++++++++++++++++++++ 8 files changed, 413 insertions(+) create mode 100644 examples/FausseTrace.jl create mode 100644 examples/GetQTest.jl rename examples/{Source_Day1.jl => SourceDay1.jl} (100%) create mode 100644 src/GenFakeData.jl create mode 100644 src/Spec.jl diff --git a/Project.toml b/Project.toml index 5a514f8..6b7af22 100644 --- a/Project.toml +++ b/Project.toml @@ -5,8 +5,12 @@ version = "0.1.0" [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" SymbolicNumericIntegration = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/docs/src/index.md b/docs/src/index.md index e123948..6103a40 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -5,4 +5,7 @@ ## Function Documentation ```@docs Ricker +Spec +ComputeQgraph +Getfreq ``` \ No newline at end of file diff --git a/examples/FausseTrace.jl b/examples/FausseTrace.jl new file mode 100644 index 0000000..f26fc4f --- /dev/null +++ b/examples/FausseTrace.jl @@ -0,0 +1,89 @@ +# Function that generates a seismic trace (Ricker function), attenuating along time, +# and plot a receiver gather (received waves at different geophone positions) + +using SeismicQ, Plots, SpecialFunctions, LinearAlgebra, Printf + +function main() + + # Functtion that create a vector of time and + function FausseTrace(x,Δt,Nt,t,Vp,Vs,αp,αs,𝑓₀,t₀) + + tₓp = x/Vp + tₓs = x/Vs + + # Storage + time = zeros(Nt) + acc = zeros(Nt) + + # Time loop + for it=1:Nt + + # Compute Ricker function + t += Δt + a = 0.5 * exp(-αp*x)*Ricker(t, tₓp+t₀, 𝑓₀)+ 0.5 * exp(-αs*x)*Ricker(t, tₓs+t₀, 𝑓₀) + + # For visualisation purpose + time[it] = t + acc[it] = a + end + + return(time',acc') + + end + + # Geophone position [m] + listₓ = 0:100:5000; + + # Time domain + Δt = 1e-3 + Nt = 2000 + # Central frequency of the source [Hz] + 𝑓₀ = 10. + t₀ = 1.0/𝑓₀ + t = -t₀ + + # Velocities + Vp = 7000 # m/s + Vs = 4000 + αp = 2e-4 + αs = 4e-4 # (π * f)/ (Q * V) + + time_axis = t:(Δt*Nt-t)/(Nt-1):Δt*Nt ; + time_vec = zeros(size(listₓ,1),Nt); + acc_vec = zeros(size(listₓ,1),Nt); + time_vec[1,:],acc_vec[1,:]= FausseTrace(listₓ[1,1],Δt,Nt,t,Vp,Vs,αp,αs,𝑓₀,t₀); + dist_axis = 0:1: size(listₓ,1)-1; # distances des geophones pour l instant 1,2,3,... + + for i=1:size(listₓ,1)-1 + time_vec[i+1,:],acc_vec[i+1,:] = FausseTrace(listₓ[i+1,1],Δt,Nt,t,Vp,Vs,αp,αs,𝑓₀,t₀) + end + + seismic_matrix = hcat(listₓ,acc_vec); + @show size(seismic_matrix) + #= + # Visualisation si peu de positions + fig1 = plot(layout = (size(listₓ,1),1)) + p1 = plot!(fig1[1],time_vec[1,:], acc_vec[1,:], xlabel="t", ylabel="a") + + for i=1:size(listₓ,1)-1 + p1 = plot!(fig1[i+1],time_vec[i+1,:], acc_vec[i+1,:], xlabel="t", ylabel="a") + end + display(fig1) + =# + + # Visualisation of receiver gather + #p2 = heatmap(dist_axis,time_axis,acc_vec', + p2 = heatmap(listₓ,time_axis,acc_vec', + color=palette(:RdBu, 100, rev=true), + clim=(-0.5, 0.5), + cbar=true, + label=" ", + yflip=true, + title = "Receiver gather", + xlabel = "Geophone position [m]", + ylabel = "Time [s]") + display(plot(p2)) # equivalent du drawnow + + +end +main() diff --git a/examples/GetQTest.jl b/examples/GetQTest.jl new file mode 100644 index 0000000..b3872cd --- /dev/null +++ b/examples/GetQTest.jl @@ -0,0 +1,86 @@ +using SeismicQ,Plots + +function main() + +#P and S velocities +Vp=7000 +Vs=4000 + +SeismicMatrix=GenMatrix(Vp,Vs) + +#distances from source +d1=3000 +d2=5000 + +#Number of samples and time increment of traces +Nt=2000 +dt=1e-3 + +#min and max frequencies to compute Q +fqmin=5.0 +fqmax=20.0 + +Tmax=dt*Nt +time=[dt*i for i=1:Nt] + +println("Trace length is $Tmax s using $Nt samples") + +ind1=Int(trunc(d1/100)) +ind2=Int(trunc(d2/100)) + +rec1=SeismicMatrix[ind1,2:Nt+1] +rec2=SeismicMatrix[ind2,2:Nt+1] + +println("Trace 1 at offset $d1 m is found at index $ind1") +println("Trace 2 at offset $d2 m is found at index $ind2") + +ampmin=-0.15 +ampmax=0.3 + +p1=plot(time,rec1,xlabel="time (s)",ylabel="acc. (m.s⁻²)",ylim=(ampmin,ampmax),c=:"blue") +p2=plot(time,rec2,xlabel="time (s)",ylabel="acc. (m.s⁻²)",ylim=(ampmin,ampmax),c=:"red") + +# second part rec1: +println("Trace 1:") + +ampP,ampS,frequP,frequS=Spec(rec1,0.1,dt,Nt,0.4,0.75,0.35) + +ymax=maximum(ampP) + +p3=plot(frequP,ampP,title="P",xlabel="f(Hz)",ylabel="Amp",ylim=(0,ymax),c=:"blue") +p4=plot(frequS,ampS,title="S",xlabel="f(Hz)",ylabel="Amp",ylim=(0,ymax),c=:"blue") + +#plot(p1,p3,p4,layout=(2,2)) + +# second part rec2: +println("Trace 2:") + +ampP2,ampS2,frequP2,frequS2=Spec(rec2,0.1,dt,Nt,0.7,1.25,0.35) + +p5=plot(frequP2,ampP2,title="P2",xlabel="f(Hz)",ylabel="Amp",ylim=(0,ymax),c=:"red") +p6=plot(frequS2,ampS2,title="S2",xlabel="f(Hz)",ylabel="Amp",ylim=(0,ymax),c=:"red") + + +@show frequP +indf=findall(x->x>fqmin && x 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 +``` +""" +function Getfreq(dt,Nt) + f=[1/(Nt*dt)*i for i=1:Nt] + return f +end + +@doc raw""" + Spec(trace,t0,dt,Nt,tp,ts,\Deltaph) + +Returns the amplitude spectrums of the P ans S phases where amplitudes are significant + +```Input + 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) + \Deltaph: width of the phases (s) +``` + +```Output + Vec2dP: modulus of the fft coeffs of P phase where larger than max(Vec2dp)/10 + Vec2dS: modulus of the fft coeffs of S phase where 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-repl +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 +``` +""" +function Spec(trace,t0,dt,Nt,tp,ts,Δph) + + # define all indexes for P and S phases + iP=trunc(Int,(tp+t0)/dt) + iS=trunc(Int,(ts+t0)/dt) + idelta=trunc(Int,Δph/dt)+1 + + # mute necessary segments for both phases + traceP=copy(trace) + traceS=copy(trace) + traceP[1:iP-1].=0. + traceS[1:iS-1].=0. + traceP[iP+idelta:Nt].=0. + traceS[iS+idelta:Nt].=0. + + println("phase P is picked at $tp s") + println("phase S is picked at $ts s") + println("phase width is set to $Δph") + # Compute Fourier coeffs + ampP=zeros(Nt) + ampcP=fft(traceP) + @. ampP=abs(ampcP) + + ampS=zeros(Nt) + ampcS=fft(traceS) + @. ampS=abs(ampcS) + + # zoom in non zero coeffs.: P + ind=findall(x->x>(maximum(ampP)/10),ampP) + + ind_halfP=zeros(Int(trunc(length(ind)/2))) + ind_halfP=ind[1:Int(trunc(length(ind)/2))] + + Vec2dP=ampP[ind_halfP] + + # zoom in non zero coeffs.: S + indS=findall(x->x>(maximum(ampS)/10),ampS) + + ind_halfS=zeros(Int(trunc(length(ind)/2))) + ind_halfS=ind[1:Int(trunc(length(ind)/2))] + + Vec2dS=ampS[ind_halfS] + + fs=Getfreq(dt,Nt) + + frequP=zeros(length(ind_halfP)) + frequP=fs[ind_halfP] + + frequS=zeros(length(ind_halfS)) + frequS=fs[ind_halfS] + + return Vec2dP,Vec2dS,frequP,frequS +end + +@doc raw""" + 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 + +```Equation + Can't write it for now (no LateX in online Editor?) +``` + +# Examples +```julia-repl +julia> TBD +``` +""" +function ComputeQgraph(amp1,amp2,freq,index,d1,d2,V) + + x=zeros(length(index)) + y=zeros(length(index)) + + + cpt=0 + for i in index + cpt += 1 + x[cpt]=freq[i] + y[cpt]=V*log(amp2[i]/amp1[i])/(π*(d2-d1)) + # Qp[cpt]=-frequP[i]*π*(d2-d1)/(Vp*log(ampP2[i]/ampP[i])) + end + + return x,y +end