ajout de l'estimation de Q

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"

docs/src/index.md

@@ -5,4 +5,7 @@
 ## Function Documentation
 ```@docs
 Ricker
+Spec
+ComputeQgraph
+Getfreq
 ```

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()

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<fqmax,frequP)
+@show indf
+
+#Qp=zeros(length(indf))
+
+xp,yp=ComputeQgraph(ampP,ampP2,frequP,indf,d1,d2,Vp)
+
+p7=plot(xp,yp,xlabel="f(Hz)",ylim=(-1.0,0.),title="Qgraph for P phase",c=:"green")
+
+xs,ys=ComputeQgraph(ampS,ampS2,frequS,indf,d1,d2,Vs)
+
+@show xs,ys
+
+p8=plot(xs,ys,xlabel="f(Hz)",ylim=(-1.0,0.),title="Qgraph for S phase",c=:"green")
+
+
+plot(p1,p2,p3,p5,p4,p6,p7,p8,layout=(4,2),legend=false)
+#plot(p7,legend=false)
+
+end
+
+main()

src/GenFakeData.jl

@@ -0,0 +1,64 @@
+# 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 GenMatrix(Vp,Vs)
+
+    # 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);
+
+    return seismic_matrix
+end

src/SeismicQ.jl

@@ -3,4 +3,10 @@ module SeismicQ
 include("Sources.jl")
 export Ricker
 
+include("Spec.jl")
+export Getfreq,Spec,ComputeQgraph
+
+include("GenFakeData.jl")
+export GenMatrix
+
 end # module SeismicQ