Plot Overthrust RTM Results

Daniel Köhn

Kiel, 02/04/2016

Import necessary packages



In [ ]:

    
from __future__ import division
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from matplotlib.colors import LightSource, Normalize
from matplotlib.pyplot import gca
from pylab import rcParams
from matplotlib import rc
import scipy.ndimage.filters
from scipy.ndimage.filters import gaussian_filter
from scipy.signal import hilbert
from mpl_toolkits.axes_grid1 import make_axes_locatable
import pickle

Activate different post-processing options



In [ ]:

    
ENVELOPE=0;
GAUSSIAN=1;

FD grid dimensions



In [ ]:

    
DH = 25.0;
NX = 800;
NY = 186;

Define Axis



In [ ]:

    
x = np.arange(0.0, DH*NX, DH)
y = np.arange(0.0, DH*NY, DH)

Define fonts



In [ ]:

    
FSize = 20
font = {'color':  'black',
        'weight': 'bold',
        'size': FSize}
mpl.rc('xtick', labelsize=FSize) 
mpl.rc('ytick', labelsize=FSize) 
rcParams['figure.figsize'] = 12, 7

Read S-wave velocity model and RTM result



In [ ]:

    
#f = open ('../start/overthrust_start_smooth2.vp')
f = open ('../FWI_GERMAINE/29_09_2016_lbfgs_app_Hessian/overthrust_fwi_smooth.vp')
data_type = np.dtype ('float32').newbyteorder ('<')
vp = np.fromfile (f, dtype=data_type)
vp = vp.reshape(NX,NY)
vp = np.transpose(vp)
vp = np.flipud(vp)



In [ ]:

    
f1 = open ('29_09_2016_fwi/modelTest_p_image.bin')
#f1 = open ('21_08_2016_smooth2/modelTest_p_image.bin')
RTM = np.fromfile (f1, dtype=data_type)
RTM = RTM.reshape(NX,NY)
RTM = np.transpose(RTM)
RTM = np.flipud(RTM)
RTM = scipy.ndimage.filters.laplace(RTM) # suppress low-wavenumber artifacts in image

Scale RTM result with depth



In [ ]:

    
RTM_scale = np.zeros((NX,NY))
RTM_scale += np.flipud(y)**2
RTM*=RTM_scale.transpose()

Apply Gaussian filter



In [ ]:

    
if(GAUSSIAN==1):
    RTM = gaussian_filter(RTM, sigma=[1,3])

Calculate Envelope



In [ ]:

    
if(ENVELOPE==1):
    analytic_signal = hilbert(RTM,axis=-1)
    RTM = np.abs(analytic_signal)
    RTM = gaussian_filter(RTM, sigma=[3,1])

Plot $\alpha$-Blending of Vp FWI result (Jet) and Laplace-filtered RTM result (Gray)



In [ ]:

    
def do_plot(n, an, title):
    ax=plt.subplot(2, 1, n)
    extent = [0.0,NX*DH/1000.0,0.0,NY*DH/1000.0]
    #cmax = 2.5e11 # Overthrust (initial model)
    cmax = 3e-2 # Overthrust (FWI)
    cmin = -cmax

    plt.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
    plt.rc('text', usetex=True)

    if(n==1):
        im1 = plt.imshow(vp, cmap=plt.cm.jet, interpolation='nearest', extent=extent)
        plt.hold(True)

    im2 = plt.imshow(RTM, cmap=plt.cm.gray, alpha=.70, interpolation='bicubic',
                 extent=extent, vmin=cmin, vmax=cmax)

    a = gca()
    a.set_xticklabels(a.get_xticks(), font)
    a.set_yticklabels(a.get_yticks(), font)
    #plt.axis('scaled')
    plt.title(title, fontdict=font)
    plt.ylabel('Depth [km]', fontdict=font)
    if(n==2):
        plt.xlabel('Distance [km]', fontdict=font)
    plt.gca().invert_yaxis()
    plt.text(0.1, 0.6,an,fontdict=font,color='white')
    #cbar=plt.colorbar()
    #cbar.set_label('Vp[m/s]', fontdict=font, labelpad=1)

Plot SubPlots



In [ ]:

    
plt.close('all')
plt.figure()
do_plot(1, '(a)', r"Overthrust-RTM (smooth FWI result)")
do_plot(2, '(b)', r" ")
#plt.savefig('test.png', format='png', dpi=100)
plt.savefig('test.pdf', bbox_inches='tight', format='pdf')
plt.tight_layout()
plt.show()



In [ ]: