Siesmic Reverse Time Migration

Saleh Al Nasser, Geophysics EAPS

Supermarket Acquistion System

juicy-watermelon.jpg

Land Seismic Acquisition System

LAND-Seismic-Acquisition-Sercel.jpg

shot_gather.jpg

In [1]:
using Distributed
using Makie
using QuantEcon
using LinearAlgebra
using Interact
Unable to load WebIO. Please make sure WebIO works for your Jupyter client.
In [2]:
addprocs(4)
Out[2]:
4-element Array{Int64,1}:
 2
 3
 4
 5
In [3]:
Threads.nthreads()
Out[3]:
4

Forward Modelling

In [4]:
@everywhere function F2d(v,model,dx,nt,dt)
    (nz,nx)=size(model)
    #data=Array(Float64, nx,nt)
   data=zeros(nx,nt)
    fdm=zeros(nz,nx,3)
    
    #Boundary Absorbing Model
    iz=1:20 
    boundary = (exp.(-( (0.015 .* (20 .- iz)).^2 ) )).^10
    boundary = boundary'
    
    #Forward-T Modeling 
    fdm[:,:,2] = model;
    data[:,1]  = model[1,:];
    
    a = (v.*dt/dx).^2;    #wave equation coefficient
    b = 2 .-4 .*a;
    
    iz=2:(nz-1)
    ix=2:(nx-1)
    izb=1:nz-20
    
    snapshot=zeros(nz,nx,nt)
 for it=2:nt

    
  fdm[iz,ix,3]=b[iz,ix].*fdm[iz,ix,2]-fdm[iz,ix,1]+a[iz,ix].*(fdm[iz,ix.+1,2]+fdm[iz,ix.-1,2]+
            fdm[iz.+1,ix,2]+fdm[iz.-1,ix,2])
        
        fdm[iz,1,3] = b[iz,1].*fdm[iz,1,2] - fdm[iz,1,1] +
        a[iz,1].*(fdm[iz,2,2] + fdm[iz.+1,1,2] + fdm[iz.-1,1,2]);
        fdm[iz,nx,3] = b[iz,nx].*fdm[iz,nx,2] - fdm[iz,nx,1] +
        a[iz,nx].*(fdm[iz,nx.-1,2] + fdm[iz.+1,nx,2] +
        fdm[iz.-1,nx,2]);
  #      
        fdm[1,ix,3] = b[1,ix].*fdm[1,ix,2] -  fdm[1,ix,1] + 
        a[1,ix].*(fdm[2,ix,2] + fdm[1,ix.+1,2] + fdm[1,ix.-1,2]);
        fdm[nz,ix,3]= b[nz,ix].*fdm[nz,ix,2]- fdm[nz,ix,1] + 
        a[nz,ix].*(fdm[nz.-1,ix,2] + fdm[nz,ix.+1,2] + fdm[nz,ix.-1,2]);
        
        
        fdm[1 ,1 ,3] = b[1 , 1].*fdm[1 ,1 ,2] -fdm[1 ,1 ,1] + 
        a[1 , 1]*(fdm[2,1,2] + fdm[1,2,2]);
        fdm[nz,1 ,3] = b[nz, 1].*fdm[nz,1 ,2] -fdm[nz,1 ,1] + 
        a[nz, 1]*(fdm[nz,2,2] +fdm[nz.-1,1,2]);
        fdm[1 ,nx,3] = b[1 ,nx].*fdm[1 ,nx,2] -fdm[1 ,nx,1] + 
        a[1 ,nx]*(fdm[1,nx.-1,2] +fdm[2,nx,2]);
        fdm[nz,nx,3] = b[nz,nx].*fdm[nz,nx,2] -fdm[nz,nx,1] + 
        a[nz,nx]*(fdm[nz.-1,nx,2] +fdm[nz,nx.-1,2]);
                        
        fdm[:,:,1] = fdm[:,:,2];
        fdm[:,:,2] = fdm[:,:,3];
#        
   
     for ixb = 1:20
         fdm[izb,ixb,1] = boundary[ixb].*fdm[izb,ixb,1];
         fdm[izb,ixb,2] = boundary[ixb].*fdm[izb,ixb,2];
         ixb2 = nx.-20 .+ixb;
         fdm[izb,ixb2,1] = boundary[nx.-ixb2.+1].*fdm[izb,ixb2,1];
         fdm[izb,ixb2,2] = boundary[nx.-ixb2.+1].*fdm[izb,ixb2,2];
         izb2 = nz.-20 .+ixb;
         fdm[izb2,:,1] = boundary[nz.-izb2.+1].*fdm[izb2,:,1];
         fdm[izb2,:,2] = boundary[nz.-izb2.+1].*fdm[izb2,:,2];
    end
        data[:,it] = fdm[1,:,2];
        snapshot[:,:,it] = fdm[:,:,2];
    #data = data[21:end-20,:]';
    end 
    return snapshot, data
    end
In [5]:
function ricker(f,n,dt,t0,t1)



T = dt*(n-1);
t = 0:dt:T;
tau = t.-t0;

    (t1,t2) = meshgrid(tau,t.-t1);
 s = (1 .-(t1.^2+t2 .^2).*f^2 .*pi^2).*exp.(-(t1.^2+t2 .^2).*pi^2 .*f^2);
        rw = s;
    return rw,t
end
Out[5]:
ricker (generic function with 1 method)
In [6]:
nz = 200; nx = 200;
dz = 5 ; dx = 5 ;
x = (0:nx-1)*dx;
z = (0:nz-1)*dz;
In [7]:
velo1 = 2000*ones(nz,nx);
velo1[Int(round(nz/2)):Int(round(nz/2)+4),Int(round(nx/2)):Int(round(nx/2)+4)] = ones(5,5).*3000;
(g,h)=size(velo1[Int(round(nz/2)):end,1:end] )
#velo[Int(round(nz/2)):end,1:end] = ones(g,h).*3000;
In [8]:
scene = Scene()
 heatmap!(scene, velo1)
scene
 axis = scene[Axis] # get axis
 axis[:names][:axisnames] = ("Distance", "Depth")
scene
Out[8]:
In [9]:
V=zeros(nz+20,nx+40)
V = [repeat(velo1[:,1],1,20) velo1 repeat(velo1[:,end],1,20)];
V = [V;(repeat(V[end,:],1,20))'];
In [10]:
# setup source function
dt = 0.9*minimum(minimum(dz./velo1/sqrt(2)));
vmin = minimum(velo1[:,:]);
nt = Int(round(sqrt((dx*nx)^2+(dz*nz)^2)/vmin/dt*1.2+1));
t = (0:nt-1)*dt;
f = 50;
In [11]:
# initial wavefield
ixs=100   

(rw,t)=ricker(f,nz+40,dt,dt*ixs,0);
rw = rw[1:nz+20,:];
    #generate shot record
       snapshot1,data1= @time F2d(V,rw,dx,nt,dt);
  7.461243 seconds (5.53 M allocations: 5.941 GiB, 19.81% gc time)

using Interact @manipulate for i=10:10:800 heatmap(snapshot[:,:,i]) end

@manipulate for ixs=1:20:200

(rw,t)=ricker(f,nz+40,dt,dt*ixs,0);

rw = rw[1:nz+20,:];

#generate shot record
   snapshot,data= F2d(V,rw,dx,nt,dt);

@manipulate for i=10:10:800 heatmap(snapshot[:,:,i]) end end

p1=heatmap(data[:,end:-1:1]) p2=heatmap(data[:,end:-1:1]) pscene = AbstractPlotting.vbox( AbstractPlotting.hbox(p1, p2), sizes = [1, 1, 1])

scene = Scene();

heat = heatmap!(scene, x,z,snapshot[end:-1:1,:,1]')[end]; N = 800 scene record(scene, "./Users/Slo0oH/Documents/Classes/animated_surface_and_wireframe.mp4",1:1:N) do i heat[3]=snapshot[end:-1:1,:,i]' end

In [49]:
scene = Scene();
heat1 = surface!(scene, x,z,snapshot1[end:-1:1,:,1]')[end];
N = 500 
scene
record(scene, "./Documents/Classes/animated_surface_and_wireframe.mp4",1:1:N) do i
    heat1[3]=snapshot1[end:-1:1,:,i]'  
end
failed process: Process(`ffmpeg -loglevel quiet -i '/var/folders/f1/rd9jgjwn3jgbfnyv0psw07480000gn/T/tmp39y39C/##video#366.mkv' -c:v libx264 -preset slow -crf 24 -pix_fmt yuv420p -c:a libvo_aacenc -b:a 128k -y ./Documents/Classes/animated_surface_and_wireframe.mp4`, ProcessExited(1)) [1]

Stacktrace:
 [1] error(::String, ::Base.Process, ::String, ::Int64, ::String) at ./error.jl:42
 [2] pipeline_error at ./process.jl:695 [inlined]
 [3] #run#505(::Bool, ::Function, ::Cmd) at ./process.jl:653
 [4] run at ./process.jl:651 [inlined]
 [5] save(::String, ::VideoStream) at /Users/Slo0oH/.julia/packages/AbstractPlotting/tmFCk/src/display.jl:273
 [6] record(::getfield(Main, Symbol("##15#16")), ::Scene, ::String, ::StepRange{Int64,Int64}) at /Users/Slo0oH/.julia/packages/AbstractPlotting/tmFCk/src/display.jl:333
 [7] top-level scope at In[49]:5
In [13]:
for k = 1:nt
    aux =  data1[:,k];
    amax = maximum((aux));
    data1[:,k] = data1[:,k]/amax;
   end
In [14]:
heatmap(data1[:,end:-1:1])
Out[14]:

scene2 = Scene(); datas=zeros(nx+40,nt)

heat = heatmap!(scene2, x,z,datas[:,:])[end]; N = 800 scene2

record(scene2, "./docs/media/animated_surface_and_wireframe.mp4",2:1:N) do i

 datas=zeros(nx+40,nt)
 datas[:,1:i]=data[:,1:i];
 heat[3]=datas[:,:].*100

end

Reverse-time Modelling

reverse.jpg

In [15]:
@everywhere function b2d(v,data,dx,nt,dt)
(nz,nx) = size(v);
(~,nt) = size(data);
fdm  = zeros(nz,nx,3);
    ss=zeros(nz,nx,3);
  
iz = 1:20;
boundary = (exp.(-( (0.015 .*(20 .-iz)).^2 ) )).^10;
    
fdm[1,:,1] = data[:,nt];
fdm[1,:,2] = data[:,nt-1];
fdm[1,:,3] = data[:,nt-2];
    
a = (v .* dt/dx) .^2;    
b = 2 .- 4 .* a;

    
ix   = 2:nx-1;       
ixb  = 1:20;         
ixb2 = nx-19:nx;
    
cz = 3;
   
snapshot = zeros(nz,nx,nt);

       for it = (nt-1):-1:1
    cz = cz .+1;
    bz = min(cz,nz);
        
        for iz = 1:bz
        fdm[iz,ixb,1] = boundary.*fdm[iz,ixb,1];
        fdm[iz,ixb,2] = boundary.*fdm[iz,ixb,2];
        fdm[iz,ixb2,1] = boundary[end:-1:1].*fdm[iz,ixb2,1];
        fdm[iz,ixb2,2] = boundary[end:-1:1].*fdm[iz,ixb2,2];
           
    end
    
            if bz >= (nz-19)
        for iz = nz-19:bz
            fdm[iz,:,1] = boundary[nz.-iz.+1].*fdm[iz,:,1];
            fdm[iz,:,2] = boundary[nz.-iz.+1].*fdm[iz,:,2];
        end
    end
        
         if bz == nz
        ez = nz .-1;
    else
        ez = bz;
    end
        
        
            iz = 1:bz;
    fdm[iz,ix,3] = fdm[iz,ix,3] - fdm[iz,ix,1];
   
    iz = 2:ez;
    fdm[iz,ix,2] = b[iz,ix].*fdm[iz,ix,1] + fdm[iz,ix,2] + a[iz,ix.+1].*fdm[iz,ix.+1,1] + 
        a[iz,ix.-1].*fdm[iz,ix.-1,1]+ a[iz.+1,ix].*fdm[iz.+1,ix,1] + a[iz.-1,ix].*fdm[iz.-1,ix,1];
    
  
    fdm[1,ix,2] = b[1,ix].*fdm[1,ix,1] + fdm[1,ix,2]+ a[1,ix.+1].*fdm[1,ix.+1,1] +
        a[1,ix.-1].*fdm[1,ix.-1,1]+ a[2,ix].*fdm[2,ix,1];
    
 
    if bz == nz
    
            
        fdm[nz,ix,2] = b[nz,ix,1].*fdm[nz,ix,1] + fdm[nz,ix,2] + a[nz,ix.+1].*fdm[nz,ix.+1,1] +
            a[nz,ix.-1].*fdm[nz,ix.-1,1] + a[nz.-1,ix].*fdm[nz.-1,ix,1];
        

        fdm[nz,1,2] = b[nz,1,1].*fdm[nz,1,1] + fdm[nz,1,2] +
            a[nz,2,1].*fdm[nz,2,1] + a[nz.-1,1,1].*fdm[nz.-1,1,1];
    end
    
    fdm[iz,1,2] = b[iz,1,1].*fdm[iz,1,1] + fdm[iz,1,2]+
        a[iz,2].*fdm[iz,2,1]+ a[iz.+1,1].*fdm[iz.+1,1,1] +
        a[iz.-1,1].*fdm[iz.-1,1,1];
    
    fdm[iz,nx,2] = b[iz,nx,1].*fdm[iz,nx,1] + fdm[iz,nx,2] + a[iz,nx.-1].*fdm[iz,nx.-1,1]+
        a[iz.+1,nx].*fdm[iz.+1,nx,1] + a[iz.-1,nx].*fdm[iz.-1,nx,1];
    

    fdm[1,1,2] = b[1,1,1].*fdm[1,1,1] + fdm[1,1,2]+
        a[1,2,1].*fdm[1,2,1] + a[2,1,1].*fdm[2,1,1];
    
  
    fdm[1,nx,2] = b[1,nx,1].*fdm[1,nx,1] + fdm[1,nx,2] +
        a[1,nx.-1,1].*fdm[1,nx.-1,1] + a[2,nx,1].*fdm[2,nx,1];
    

    fdm[:,:,1] = fdm[:,:,2];
    fdm[:,:,2] = fdm[:,:,3];
    
  
    if it > 2
        fdm[2:nz,:,3] = zeros(nz-1,nx);
        fdm[1,:,3] = data[:,it-2];
    end  
    
    snapshot[:,:,it] = fdm[:,:,1];
        end 
    model = fdm[:,:,1];
    return snapshot,model
    end
In [16]:
nz = Int(200); nx = Int(200);
dz = 10 ; dx = 10 ;
x = (0:nx-1)*dx;
z = (0:nz-1)*dz;
In [17]:
velo=zeros(nz,nx)
velo = 2000 .*ones(nz,nx);
velo[Int(round(nz/2)):Int(round(nz/2)+4),Int(round(nx/2)):Int(round(nx/2)+4)] = ones(5,5).*3000;
velo_const = 2000 .*ones(nz,nx);        
#velo[51:end,1:end] = 3000*ones(50,nx);
#velo[76:end,1:end] = 4000*ones(25,nx);
In [18]:
V = [repeat(velo[:,1],1,20) velo repeat(velo[:,end],1,20)];
V = [V;(repeat(V[end,:],1,20))'];
In [19]:
heatmap(V[end:-1:1,:]')
Out[19]:
In [20]:
Vc = [repeat(velo_const[:,1],1,20) velo_const repeat(velo_const[:,end],1,20)];
Vc = [Vc;(repeat(Vc[end,:],1,20))'];
In [21]:
# setup source function
dt = 0.9*minimum(minimum(dz./velo/sqrt(2)));
vmin = minimum(velo[:,:]);
nt = Int(round(sqrt((dx*nx)^2+(dz*nz)^2)/vmin/dt+1));
t = (0:nt-1)*dt;
f = 1000;
In [22]:
nt
Out[22]:
668
In [23]:
data = zeros(size(nt,nx));
snapshot_i = zeros(nz+20,nx+40,nt);
#snapshot_d = zeros(nz+20,nx+40,nt);
#snapshot_c = zeros(nz+20,nx+40,nt);
images = zeros(nz+20,nx+40);

#for ixs=1:10:100;
 ixs=100  
    (rw,t) = ricker(f,nz+40,dt,dt*ixs,0);
    rw = rw[1:nz+20,:];
    
 

  
    (snapshot_s,data) = F2d(V,rw,dx,nt,dt);
      ( snapshot_c,data_const) = F2d(Vc,rw,dx,nt,dt);
  
    
 data_refl = data .- data_const;


   @time    (snapshot_d,fdm) = b2d(Vc,data_refl,dx,nt,dt);
    
   
      for i = 1:nt
        snapshot_i[:,:,i] = snapshot_c[:,:,i] .* snapshot_d[:,:,i];
    end
   i=2
for i = 2:nt
        snapshot_i[:,:,i] = snapshot_i[:,:,i] .+ snapshot_i[:,:,i.-1];
    end
    images = images .+ snapshot_i[:,:,end];

#end
  8.482021 seconds (4.67 M allocations: 6.214 GiB, 14.72% gc time)

using SharedArrays

using DistributedArrays a = SharedArray{Float64}(10) @distributed for i = 1:10 a[i] = i end

In [50]:
scene = Scene();

heat = heatmap!(scene, x,z,snapshot_c[end:-1:1,:,1]')[end];
N = 200
scene
record(scene, "./docs/media/animated_surface_and_wireframe.mp4",1:1:N) do i
    heat[3]=snapshot_c[end:-1:1,:,i]'
end
failed process: Process(`ffmpeg -loglevel quiet -i '/var/folders/f1/rd9jgjwn3jgbfnyv0psw07480000gn/T/tmp2Tu1jq/##video#368.mkv' -c:v libx264 -preset slow -crf 24 -pix_fmt yuv420p -c:a libvo_aacenc -b:a 128k -y ./docs/media/animated_surface_and_wireframe.mp4`, ProcessExited(1)) [1]

Stacktrace:
 [1] error(::String, ::Base.Process, ::String, ::Int64, ::String) at ./error.jl:42
 [2] pipeline_error at ./process.jl:695 [inlined]
 [3] #run#505(::Bool, ::Function, ::Cmd) at ./process.jl:653
 [4] run at ./process.jl:651 [inlined]
 [5] save(::String, ::VideoStream) at /Users/Slo0oH/.julia/packages/AbstractPlotting/tmFCk/src/display.jl:273
 [6] record(::getfield(Main, Symbol("##17#18")), ::Scene, ::String, ::StepRange{Int64,Int64}) at /Users/Slo0oH/.julia/packages/AbstractPlotting/tmFCk/src/display.jl:333
 [7] top-level scope at In[50]:6
In [51]:
scene = Scene();

heat = heatmap!(scene, x,z,snapshot_d[end:-1:1,:,1]')[end];
N = 200
scene
record(scene, "./docs/media/animated_surface_and_wireframe.mp4",1:1:N) do i
    heat[3]=snapshot_d[end:-1:1,:,i]'
end
failed process: Process(`ffmpeg -loglevel quiet -i '/var/folders/f1/rd9jgjwn3jgbfnyv0psw07480000gn/T/tmpU8AEgF/##video#369.mkv' -c:v libx264 -preset slow -crf 24 -pix_fmt yuv420p -c:a libvo_aacenc -b:a 128k -y ./docs/media/animated_surface_and_wireframe.mp4`, ProcessExited(1)) [1]

Stacktrace:
 [1] error(::String, ::Base.Process, ::String, ::Int64, ::String) at ./error.jl:42
 [2] pipeline_error at ./process.jl:695 [inlined]
 [3] #run#505(::Bool, ::Function, ::Cmd) at ./process.jl:653
 [4] run at ./process.jl:651 [inlined]
 [5] save(::String, ::VideoStream) at /Users/Slo0oH/.julia/packages/AbstractPlotting/tmFCk/src/display.jl:273
 [6] record(::getfield(Main, Symbol("##19#20")), ::Scene, ::String, ::StepRange{Int64,Int64}) at /Users/Slo0oH/.julia/packages/AbstractPlotting/tmFCk/src/display.jl:333
 [7] top-level scope at In[51]:6
In [52]:
scene = Scene();

heat = heatmap!(scene, x,z,snapshot_i[end:-1:1,:,1]')[end];
N = 400
scene
record(scene, "./docs/media/animated_surface_and_wireframe.mp4",1:1:N) do i
    heat[3]=snapshot_i[end:-1:1,:,i]'
end
failed process: Process(`ffmpeg -loglevel quiet -i '/var/folders/f1/rd9jgjwn3jgbfnyv0psw07480000gn/T/tmplz9UFs/##video#370.mkv' -c:v libx264 -preset slow -crf 24 -pix_fmt yuv420p -c:a libvo_aacenc -b:a 128k -y ./docs/media/animated_surface_and_wireframe.mp4`, ProcessExited(1)) [1]

Stacktrace:
 [1] error(::String, ::Base.Process, ::String, ::Int64, ::String) at ./error.jl:42
 [2] pipeline_error at ./process.jl:695 [inlined]
 [3] #run#505(::Bool, ::Function, ::Cmd) at ./process.jl:653
 [4] run at ./process.jl:651 [inlined]
 [5] save(::String, ::VideoStream) at /Users/Slo0oH/.julia/packages/AbstractPlotting/tmFCk/src/display.jl:273
 [6] record(::getfield(Main, Symbol("##21#22")), ::Scene, ::String, ::StepRange{Int64,Int64}) at /Users/Slo0oH/.julia/packages/AbstractPlotting/tmFCk/src/display.jl:333
 [7] top-level scope at In[52]:6
In [27]:
p1 = heatmap(snapshot_c[end:-1:1,:,100]')
p2 = heatmap(snapshot_d[end:-1:1,:,100]')
p3 = heatmap(snapshot_i[end:-1:1,:,100]')
p4 = heatmap(snapshot_c[end:-1:1,:,200]')
p5 = heatmap(snapshot_d[end:-1:1,:,200]')
p6 = heatmap(snapshot_i[end:-1:1,:,200]')
p7 = heatmap(snapshot_c[end:-1:1,:,400]')
p8 = heatmap(snapshot_d[end:-1:1,:,400]')
p9 = heatmap(snapshot_i[end:-1:1,:,400]')
 t = Theme(align = (:left, :bottom), raw = true, camera = campixel!)
 title1 = text(t, "Forward modelling")
 title2 = text( t,"Reverse-time")
 title3 = text( t,"Crosscorrelation")
 pscene = AbstractPlotting.vbox(AbstractPlotting.hbox(p7,p4,p1,title1), AbstractPlotting.hbox(p8,p5,p2,title2),
    AbstractPlotting.hbox(p9,p6,p3,title3) 
 )
Out[27]: