notebook.community

Edit and run



In [3]:

    
using Seismic, PyPlot

# Explore application of the parabolic Radon transform on a real marine CMP gather



In [5]:

    
# 1- Download data set from web site (1 NMO-corrected CMP gather from the GOM)

download("http://seismic.physics.ualberta.ca/data/gom_cdp_nmo.su",
         "gom_cdp_nmo.su")

# 2- Convert SU to internal SeismicJulia format

SegyToSeis("gom_cdp_nmo.su", "gom_cdp_nmo.seis", format="su", input_type="ieee",
           swap_bytes=true)

# 3- Read data d and trace headers h

d, h = SeisRead("gom_cdp_nmo.seis")
i1 = 801; i2 = 1201
d = d[i1:i2, :]       # Use subset of the data
imute = findin(d, 0)  # Extract indexes of muted data
offset = Seismic.ExtractHeader(h, "h")
offset = offset*0.3048  # I like SI units 
href = maximum(abs.(offset))
nt = size(d, 1)
dt = h[1].d1









    



  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0





    



number of traces: 92
number of samples per trace: 1751






    



100  650k  100  650k    0     0   9.9M      0 --:--:-- --:--:-- --:--:-- 10.0M






    Out[5]:





0.004f0



In [7]:

    
# 4- Model data in the Radom domain via inersion with parabolic Radon transform
p  = collect(linspace(-0.3, 0.9, 200))
param = Dict(:order=>"parab", :dt=>dt, :p=>p, :h=>offset, :href=>href,
             :flow=>0.0, :fhigh=>125.0)
m = SeisRadonFreqInv(d; mu=0.01, param...)

# 5- Filter primaries and keep multiples in Radon gather
mf = copy(m)
pcut = 0.1 
icut = find(p .< pcut)
mf[:, icut] = 0

# 6- Forward model multiples from Radon panel 
d_mult = SeisRadonFreqFor(mf, nt; param...)
   
# 7- Compute primaries = input data - modelled multiples
d_prim = d - d_mult
d_prim[imute] = 0;



In [8]:

    
# 8- Plot input data, Radon panel, and estimated primaries

figure(1, figsize=(9,6));

t1 = (i1-1)*dt

subplot(131)
SeisPlot(d, title="Input CMP", xlabel="Offset [m]", ylabel="Time [s]",
         cmap="PuOr", ox=offset[1], dx=offset[2]-offset[1], oy=t1, dy=dt,
         fignum=1)

subplot(132)
SeisPlot(m[1:nt, :], title="Radon panel", xlabel="Residual moveout [s]",
         cmap="PuOr", ox=p[1], dx=p[2]-p[1], oy=t1, dy=dt, fignum=1)

subplot(133)
SeisPlot(d_prim, title="Estimated primaries", xlabel="Offset [m]", cmap="PuOr", 
         ox=offset[1], dx=offset[2]-offset[1], oy=t1, dy=dt, fignum=1)

tight_layout()



In [27]:

    
# 9- Compare near offsets wiggles for QC

figure(2, figsize=(9,6));

t1 = (i1-1)*dt


SeisPlot(hcat(d[:,1:15], zeros(nt,2), d_prim[:,1:15]), style="wiggles",
         title="Original                Primaries", xlabel="Trance number", ylabel="Time [s]",
         cmap="PuOr", ox=1, dx=1, oy=t1, dy=dt, fignum=1)









    












    Out[27]:





1-element Array{PyCall.PyObject,1}:
 PyObject <matplotlib.lines.Line2D object at 0x138d85ad0>