Calculate TD shot gather from FD modelling results
Daniel Köhn
Kiel, 01.09.2017

Import necessary packages


In [ ]:
import numpy as np
import matplotlib as plt
import matplotlib.pyplot as plt
from matplotlib.pyplot import gca
import matplotlib as mpl
from pylab import rcParams
from matplotlib import rc
from scipy import interpolate
from scipy import signal
#%matplotlib notebook

Define parameters for data visualization


In [ ]:
clip = 5e-3   # data clip
nsrc = 1      # number of shots
ntr = 117     # number of receivers
shotno = 1    # extract FDFD data for shot shotno

# number of frequencies nf and frequency interval df
FC_low = 1e6
FC_high = 800e6

nf = 400
df = (FC_high - FC_low) / (nf-1)
nfnsrc = nf * nsrc

# TD parameters
TmaxTD = 150e-9       # maximum time TD
TmaxFD = 0.25 / df    # maximum time FD
dt = TmaxFD / nf      # time sampling FD

# maximum time sample of FD2TD corresponding to TmaxTD
nmaxFD = np.int(TmaxTD / dt)

# define time samples in time and frequency domain
tFD = np.arange(0.0, TmaxTD, dt)
offset = np.arange(0, ntr, 1)

Font properties


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

Time Damping


In [ ]:
S = 0.0

Read 2 cross TE-mode FD data


In [ ]:
name_FDFD = "../seis/2_cross_TE_p_stage_1.bin"
f = open (name_FDFD)
data_type = np.dtype ('float32').newbyteorder ('<')
data_FDFD = np.fromfile(f, dtype=data_type)
data_FDFD = data_FDFD.reshape(nf*nsrc,2*ntr)
data_FDFD = np.transpose(data_FDFD)
data_FDFD = np.flipud(data_FDFD)

Extract FD data for shot shotno


In [ ]:
FD_real = data_FDFD[0:2*ntr-1:2,shotno-1:nfnsrc:nsrc]
FD_imag = data_FDFD[1:2*ntr:2,shotno-1:nfnsrc:nsrc]
data_FDFD = None

# taper FD data
window = signal.tukey(nf,1)

for i in range (ntr):
    FD_real[i,:] *= window
    FD_imag[i,:] *= window

In [ ]:
# assemble complex data
FDFD = FD_real + 1j*FD_imag
FD_real = None
FD_imag = None
FDFD = np.transpose(FDFD)

tmp1 = np.concatenate((FDFD, np.zeros((nf, ntr)), np.zeros((nf, ntr)), np.flipud(FDFD)))
#tmp1 = np.concatenate((np.zeros((1, ntr)),FDFD,np.flipud(FDFD)))
FDFD = None
FDFD = tmp1

In [ ]:
#plt.imshow(np.real(FDFD))
#plt.show()

IFFT of FD data


In [ ]:
# transformation from FD to TD
FD2TD = np.fft.ifft(FDFD,axis=0)
FD2TD = np.real(FD2TD)
tmp1 = None

Extract FD2TD up to TmaxTD and resample FD2TD from dtFD to dt


In [ ]:
# extract FD2TD data up to TmaxTD
tmp1 = FD2TD[1:nmaxFD,:]
FD2DTD = None

# clean memory
tmp = None

In [ ]:
#plt.imshow(FD2TD, vmin=-clip, vmax=clip)
#plt.show()

In [ ]:
#print(np.shape(tFD), np.shape(offset), np.shape(tmp1))

Damp Time Domain Data


In [ ]:
#tmp2 = np.exp(-S * tTD)
#tdamp = np.tile(tmp2,(ntr,1))
#tdamp = np.flipud(tdamp.T)
#traces = traces * tdamp

#tmp2 = None
#tdamp = None

Normalize data


In [ ]:
#max_TD = np.max(traces)
max_FD = np.max(tmp1)
tmp1 = tmp1 / max_FD
FD2TD = np.flipud(tmp1)

tmp1 = None

Shuffle 2-cross TE model data


In [ ]:
tmp1 = FD2TD[:,0:29]
tmp2 = FD2TD[:,30:57]
tmp3 = FD2TD[:,58:87]
tmp4 = FD2TD[:,88:116]

FD2TD = None

tmp = np.concatenate((tmp1.T, tmp4.T, np.flipud(tmp2.T), np.flipud(tmp3.T)))
FD2TD = tmp.T

tmp = None
tmp1 = None
tmp2 = None
tmp3 = None
tmp4 = None

#plt.imshow(FD2TD)
#plt.show()

Plot shot gather


In [ ]:
def do_plot(n, model, cm, an, title, clip):
    
    extent = [np.min(offset), np.max(offset), 0.0, dt * nmaxFD * 1e9]

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

    im = plt.imshow(model, cmap=cm, interpolation='nearest', extent=extent, vmin=-clip, vmax=clip, aspect=0.5)

    a = gca()
    a.set_xticklabels(a.get_xticks(), font)
    a.set_yticklabels(a.get_yticks(), font)

    plt.title(title, fontdict=font)
    if n==1:
        plt.ylabel('Time [ns]', fontdict=font)
    plt.xlabel('Trace no.', fontdict=font)
    plt.gca().invert_yaxis()

In [ ]:
plt.close('all')
plt.figure()
do_plot(1, FD2TD, 'gray', '(b)', r"2-cross TE-mode model (GERMAINE)", clip)
plt.tight_layout()
output_file = "test.pdf"
plt.savefig(output_file, bbox_inches='tight', format='pdf')
plt.show()

In [ ]: