In [1]:

    

using Seismic, PyPlot # required deps


    

In [2]:

    

# First let's view an inline out of the "real" seismic data:

seismic, seismic_h = SeisRead("dat/image"); # load data
SeisPlot(seismic[:,:,1], cmap = "seismic")  # plot data


    

    













    
Out[2]:





PyObject <matplotlib.image.AxesImage object at 0x7f0c89b0c710>

In [3]:

    

# Then let's assume we've picked a velocity model
# and converted it to interval velocities.
# We'll assume there are interpretive errors in
# the model.

# An inline from our picked vel model looks like this:

vel_initial, vel_h = SeisRead("dat/vel_incorrect"); # read data
SeisPlot(vel_initial[:,:,1], cmap = "magma")        # plot data


    

    













    
Out[3]:





PyObject <matplotlib.image.AxesImage object at 0x7f0c892efd30>

In [4]:

    

# Obviously, this model is incorrect, or there would be three
# seismic events in the section. Let's forward model this 
# incorrect vel model and view the output modeled seismic:

include("mod/modeler.jl");                                   # model seismic from our bad vels
seis_incorrect, seismic_h = SeisRead("dat/image_incorrect"); # load data
SeisPlot(seis_incorrect[:,:,1], cmap = "seismic")            # plot data


    

    




grep: ../dat/vel_correct: No such file or directory






    




LoadError: LoadError: failed process: Process(`grep in= ../dat/vel_correct`, ProcessExited(2)) [2]
while loading /home/gram/tle_fwi/mod/modeler.jl, in expression starting on line 11
while loading In[4], in expression starting on line 5

 in pipeline_error at process.jl:568
 in readbytes at process.jl:515
 in readall at process.jl:520
 in SeisRead at /home/gram/.julia/v0.4/Seismic/src/Utils/SeisRead.jl:33
 in include at ./boot.jl:261
 in include_from_node1 at ./loading.jl:320

In [5]:

    

# Therefore, if we tried to migrate with our bad velocity 
# model, we'd get a kinematically inaccurate result. Hence
# we need to improve our velocity model somehow. So let's
# use full waveform inversion.

# Let's input our poorly picked velocity model and the 
# seismic data into an FWI algorithm and take a look
# at the updated velocity model to see if there was any
# improvement:

include("fwi/fwi.jl")                            # prep deps
fwi("dat/vel_incorrect", "dat/image", 
"dat/wav", "dat/updated_vel", .03, 5, 0);     # FWI algorithm

vel_updated, vel_h = SeisRead("dat/updated_vel") # load data
SeisPlot(vel_updated[:,:,1], cmap = "magma")        # plot data


    

    




3D Poststack full waveform inversion algorithm


In this program FWI updates a velocity model one sample at a time,
one trace at a time. The velocity sample is perturbed positively,
negatively, and not at all. Then the algorithm models new seismic
from all three of the perturbed velocity models, and compares them to
the input (real) seismic data. Whichever perturbation generates the
least amount of error is chosen to be correct, and input into the 
updated velocity model.
This algorithm accepts six inputs in the following order:

	1. Initial velocity model
	2. Poststack 3D seismic volume in Seis format
	3. One dimensional wavelet in Seis format
	4. Output file name (updated vel model in Seis format)
	5. Velocity update increment percentage in decimal
	6. Maximum number of velocity update iterations
	7. Verbosity (see note)

If verbose is 0 operation is silent. If verbose is 1 updates will
print to stdout. If verbose is 2 debugging info will print to stdout.
Here is an example to get you started on the syntax:


fwi("../dat/vel_incorrect", "../dat/image", "../dat/wav", "../dat/updated_vel", .1, 5, 0)


Please note this algorithm is written for clarity not speed, so
use forgivingly small velocity models.







    




grep: ../dat/vel_incorrect: No such file or directory






    




LoadError: failed process: Process(`grep in= ../dat/vel_incorrect`, ProcessExited(2)) [2]
while loading In[5], in expression starting on line 12

 in pipeline_error at process.jl:568
 in readbytes at process.jl:515
 in readall at process.jl:520
 in SeisRead at /home/gram/.julia/v0.4/Seismic/src/Utils/SeisRead.jl:33
 in fwi at /home/gram/tle_fwi/fwi/fwi.jl:37

In [6]:

    

# We can see here that the vel model has been improved considerably
# by full waveform inversion. For comparison, let's view
# the insitu velocities.
# (This is the model the author used to create the "real" seismic)


vel_true, vel_h = SeisRead("dat/vel_correct"); # read data
SeisPlot(vel_true[:,:,1], cmap = "magma")         # plot data


    

    













    
Out[6]:





PyObject <matplotlib.image.AxesImage object at 0x7f0c84b77710>

In [7]:

    

# As you can see, the FWI updated vel model is very close
# to the true velocities. In this case it's actually within 3%.


    

In [ ]: