preprocess for RTM

toolbox/base.py

@@ -318,7 +318,10 @@ def __screenprint(self):
 
         # pick and mute
         if self.mute_late_arrival:
-            print('Data processing  : time window, %.2f s after the first break' % self.mute_late_window)
+            if self.mute_late_window > 0:
+                print('Data processing  : time window, %.2f s after the first break' % self.mute_late_window)
+            else:
+                print('Data processing  : time window, %.2f s mute first arrivals for RTM' % abs(self.mute_late_window))
         else:
             print('Data processing  : no time window')
 

toolbox/fielddata.py

@@ -0,0 +1,221 @@
+###############################################################################
+#
+# SWIT: Seismic Waveform Inversion Toolbox
+#
+# by Haipeng Li at USTC, [email protected]
+# 
+# June, 2021  
+#
+# Field data module
+#
+###############################################################################
+
+
+from multiprocessing import Pool, cpu_count
+import numpy as np
+import obspy
+
+from plot import plot_wavelet
+from tools import array2su, get_offset, get_su_parameter, loadsu, savesu, smooth1d, su2array
+
+
+def load_field_data(simu, data_list, x_beg, x_end):
+    ''' Load field data (pre-processed), vertical geophone
+    '''
+
+    nproc    = np.min((simu.system.mpiproc, cpu_count()//2))
+    homepath = simu.system.homepath
+
+    # get prepared
+    srcn = len(data_list)
+    dt = simu.model.dt
+    nt = simu.model.nt
+
+    print('Field data: loading data and adjusting to current acquisition')
+
+    pool = Pool(nproc)
+    receiver_range = [pool.apply_async(load_field_data_serial, (homepath, data_list[isrc], isrc, x_beg, x_end, dt, nt, )) for isrc in range(srcn)]
+    pool.close()
+    receiver_range = np.array([p.get() for p in receiver_range])
+
+    pool.join()
+
+    simu.receiver.rec_beg = receiver_range[:,0]
+    simu.receiver.rec_end = receiver_range[:,1]
+
+    print('Field data: %d shotgathers in total' %(srcn))
+    print('Field data: %d time steps, %.2f ms samping interval\n'%(nt, dt*1000))
+    print('*****************************************************\n')
+
+
+
+def load_field_data_serial(homepath, loadfile, isrc, x_beg, x_end, dt, nt):
+
+        # set path
+        savefile = homepath + 'data/obs/src%d_sg.su' %(isrc + 1)
+
+        # load data
+        trace = loadsu(loadfile)
+        _, _, dt_data = get_su_parameter(trace)
+
+        # select traces in calculate domain
+        trace = obspy.Stream(traces = [tr for tr in trace
+                if  tr.stats.su.trace_header.group_coordinate_x - x_beg >= 0.0
+                and tr.stats.su.trace_header.group_coordinate_x - x_end <= 0.0])
+
+        # modify header to adjust to current acquisition geometry and for obspy plotting
+        for tr in trace:
+            tr.stats.su.trace_header.group_coordinate_x  = int(tr.stats.su.trace_header.group_coordinate_x  - x_beg)
+            tr.stats.su.trace_header.source_coordinate_x = int(tr.stats.su.trace_header.source_coordinate_x - x_beg)
+            tr.stats.distance = tr.stats.su.trace_header.group_coordinate_x - tr.stats.su.trace_header.source_coordinate_x
+
+        # interpolate if necessary
+        if dt_data != dt:
+            for tr in trace:
+                tr.resample(1.0/dt)
+
+        # ensure the same time steps
+        for tr in trace:
+            t0 = tr.stats.starttime 
+            tr.trim(starttime=t0, endtime=t0 + dt*(nt-1), pad=True, nearest_sample=True, fill_value=0.0)
+
+        # save file to the working folder
+        savesu(savefile, trace)
+
+        # record the receiver range
+        x0 = trace[ 0].stats.su.trace_header.group_coordinate_x
+        x1 = trace[-1].stats.su.trace_header.group_coordinate_x
+
+        return np.array([x0, x1])
+
+
+def source_wavelet_process(stf, stf_dt=0.001, use_dt = 0.001, use_nt=4001, shift = 0.0, lowpass = 0, highpass = 50, taper_beg=0.005):
+    ''' source wavelet process
+    '''
+    # convert to su
+    stf = array2su(1, stf_dt, stf)
+    t0 = stf[0].stats.starttime 
+    stf.resample(1.0/use_dt)
+    stf.detrend('demean')
+    stf.detrend('linear')
+    # filter
+    if lowpass == 0:
+        stf.filter('lowpass', freq=highpass)#, corners=4, zerophase=True)
+    elif 0 < lowpass < highpass:
+        stf.filter('bandpass', freqmin=lowpass, freqmax=highpass)#, corners=4, zerophase=True)
+    else:
+        pass
+    stf.trim(starttime=t0, endtime=t0 + use_dt*(use_nt-1), 
+             pad=True, nearest_sample=True, fill_value=0.)
+    stf[0].data[:] = smooth1d(stf[0].data[:], span = 5)
+    # apply shift
+    if shift > 0:
+        shift = int(shift//use_dt)
+        stf[0].data[0:-1-shift] =  stf[0].data[shift:-1]
+
+    stf.taper(taper_beg, type='hann', side='left')
+
+    stf = stf[0].data[:]
+    if stf.size != use_nt:
+        ValueError('wrong size of the source wavelet')
+
+    return stf
+
+
+def source_inversion(simu, inv_offset=10000):
+    ''' source signature inversion by (Parrat, 1999)
+    '''
+    # get prepared
+    homepath = simu.system.homepath
+    stf_now = simu.source.wavelet
+    srcx = simu.source.xyz[:,0]
+    recx = simu.receiver.xyz[:,:,0]
+    srcn = simu.source.n
+    dt = simu.model.dt
+
+    stf_inv = np.zeros_like(stf_now)
+
+    for isrc in range(srcn):
+
+        # load data and set offset
+        obs = loadsu(homepath + 'data/obs/src%d_sg_proc.su'%(isrc+1))
+        syn = loadsu(homepath + 'data/syn/src%d_sg_proc.su'%(isrc+1))
+        if not np.array_equal(get_offset(obs), get_offset(syn)):
+            raise ValueError("offset not consistant.")
+        offset = get_offset(syn)
+        obs = su2array(obs)
+        syn = su2array(syn)
+
+        # select data for source inversion
+        rec_used = np.argwhere(abs(offset) < inv_offset)
+        n1 = int(rec_used[0])
+        n2 = int(rec_used[-1])
+        data_obs = np.squeeze(obs[n1:n2, :]).T
+        data_syn = np.squeeze(syn[n1:n2, :]).T
+        Do = np.fft.fft(data_obs, axis=0)  # frequency domain
+        Dm = np.fft.fft(data_syn, axis=0)  # frequency domain
+
+        # current source wavelet
+        src = np.squeeze(stf_now[isrc, :])
+        S = np.fft.fft(np.squeeze(src), axis=0) # frequency domain
+
+        # check
+        if abs(np.sum(Dm * np.conj(Dm))) == 0:
+            raise ValueError("No trace for source inversion, check for the reason.")
+
+        # source inversion
+        A = np.sum(np.conj(Dm)*Do, axis=1) / np.sum(Dm * np.conj(Dm), axis=1)
+        temp = np.real(np.fft.ifft(A*S[:]))
+        temp = temp / np.max(abs(temp))
+        stf_inv[isrc, :] = temp
+
+    # process the source wavelet
+    stf_now = convert_wavelet_su(dt, stf_now, srcx)
+    stf_inv = convert_wavelet_su(dt, stf_inv, srcx)
+
+    # save source wavelet plots
+    plot_wavelet(simu, stf_now, 'stf_now', scale=1.0, color='k', plot_dx=5000, t_end = 1.0)
+    plot_wavelet(simu, stf_inv, 'stf_inv', scale=1.0, color='r', plot_dx=5000, t_end = 1.0)
+
+    # save source wavelet data
+    np.savetxt(homepath+'outputs/stf_now.dat', su2array(stf_now))
+    np.savetxt(homepath+'outputs/stf_inv.dat', su2array(stf_inv))
+
+    print('Source inversion finished\n')
+
+
+def set_frequency_band(zmax, hmax, f_begin, band_num):
+    ''' Multi-scale implemation based on Boonyasiriwat et al, 2009.
+        f(n+1)    = f(n) / alpha_min,       where 
+        alpha_min = z/sqrt( h**2 + z**2),   and 
+        z is maximum depth to be imaged,    and
+        h is maximum half-offset.
+    '''
+
+    alpha_min = zmax / np.sqrt(zmax**2 + hmax**2)
+    fre_band = np.zeros(band_num)
+    fre_band[0] = f_begin
+    for iband in range(band_num):
+        if iband > 0:
+            fre_band[iband] = fre_band[iband-1] / alpha_min
+
+    # set two decimal
+    fre_band = np.around(fre_band, decimals=2)
+
+    return fre_band
+
+
+def convert_wavelet_su(dt, wavelet, srcx):
+    ''' Convert wavelet array to the SU streame
+    '''
+    srcn  = np.size(wavelet, 0)
+
+    wavelet_su = array2su(srcn, dt, wavelet)
+
+    ishot = 0
+    for iwvlt in wavelet_su:
+        iwvlt.stats.distance = srcx[ishot]
+        ishot+=1
+
+
+    return wavelet_su

toolbox/inversion.py

@@ -47,8 +47,8 @@ def inversion(simu, optim, inv_model):
         # synthetic data from the current model
         forward(simu, simu_type='syn', savesnap=1)
 
-        plot_trace(simu, 'syn', simu_type='syn', suffix='', src_space=1, trace_space=5, scale = 0.8, color='b')
-        plot_trace(simu, 'syn-proc', simu_type='syn', suffix='_proc', src_space=1, trace_space=5, scale = 0.8, color='b')
+        # plot_trace(simu, 'syn', simu_type='syn', suffix='', src_space=1, trace_space=5, scale = 0.8, color='b')
+        # plot_trace(simu, 'syn-proc', simu_type='syn', suffix='_proc', src_space=1, trace_space=5, scale = 0.8, color='b')
 
         # process the synthetic data
         process_workflow(simu, optim, simu_type='syn')
@@ -97,15 +97,10 @@ def rtm(simu, optim, inv_model):
     # synthetic data from the current model
     forward(simu, simu_type='syn', savesnap=1)
 
-    plot_trace(simu, 'syn', simu_type='syn', suffix='', src_space=1, trace_space=5, scale = 0.8, color='b')
-    plot_trace(simu, 'syn-proc', simu_type='syn', suffix='_proc', src_space=1, trace_space=5, scale = 0.8, color='b')
-
-    # process the synthetic data
-    process_workflow(simu, optim, simu_type='syn')
-
     # compute the RTM image
     inv_scheme['g_now'] = adjoint(simu, optim) 
 
+    # plot and save
     plot_rtm(simu, optim, inv_scheme)
 
     print('\n-----------  RTM end  -----------\n')
@@ -285,7 +280,7 @@ def backtrack(simu, optim, inv_scheme):
     ftol = 1e-4        # control the accuracy of the line search routine
     wolfe = 0.9        # coefficient for the Wolfe condition
     search_max = 6     # line search max iteration
-    vmax_thresh = 200  # when direction is too large
+    vmax_thresh = 300  # when direction is too large
     vmin_thresh = 20   # when direction is too small
 
     # get parameters

toolbox/misfit.py

@@ -94,10 +94,13 @@ def adjoint_source_serial(homepath, isrc, misfit_type):
     '''
 
     obs = loadsu(homepath + 'data/obs/src%d_sg_proc.su'%(isrc+1))
-    syn = loadsu(homepath + 'data/syn/src%d_sg_proc.su'%(isrc+1))
 
-    if get_su_parameter(obs) != get_su_parameter(syn):
-        raise ValueError('obs and syn are not consistent.')
+    # we do not use syn data in RTM
+    if misfit_type.lower() not in ['rtm']:
+        syn = loadsu(homepath + 'data/syn/src%d_sg_proc.su'%(isrc+1))
+
+        if get_su_parameter(obs) != get_su_parameter(syn):
+            raise ValueError('obs and syn are not consistent.')
 
     # Waveform difference L2-norm (Tarantola, 1984)
     if misfit_type.lower() in ['waveform']:
@@ -117,7 +120,7 @@ def adjoint_source_serial(homepath, isrc, misfit_type):
 
     # Reverse Time Migration
     elif misfit_type.lower() in ['rtm']:
-        adj = adjoint_source_rtm(obs, syn)
+        adj = adjoint_source_rtm(obs)
 
     return adj
 
@@ -240,7 +243,7 @@ def adjoint_source_global_correlation(obs, syn):
 
 
 
-def adjoint_source_rtm(obs, syn):
+def adjoint_source_rtm(obs):
     ''' Reverse Time Migration
     '''
     # parameters
@@ -250,6 +253,10 @@ def adjoint_source_rtm(obs, syn):
     for irec in range(recn):
         adj[irec,:] = obs[irec].data
 
+        # # differentiate the seismic trace twice
+        # adj[irec,:] = np.diff(adj[irec,:], prepend = adj[irec,0])
+        # adj[irec,:] = np.diff(adj[irec,:], prepend = adj[irec,0])
+
     return adj
 
 

toolbox/plot.py

@@ -18,7 +18,7 @@
 import matplotlib.pyplot as plt
 from multiprocessing import Pool
 
-from tools import loadsu, add_su_header, convert_wavelet_su
+from tools import loadsu, add_su_header, convert_wavelet_su, savebinfloat32
 
 ### Plot acquisition geometry
 def plot_geometry(simu):
@@ -251,12 +251,12 @@ def plot_inv_scheme(simu, optim, inv_scheme):
     plot_model2D(simu, dire.reshape(nx, nz).T, -dirc_caxis, dirc_caxis, 'dire-%03d' % it, colormap = 'seismic')
 
     if optim.iter == 1 :
-        plot_trace(simu, 'syn-proc-initial-model', simu_type = 'syn', suffix='_proc', src_space=1, trace_space=5, scale=0.8, color='r')
+        plot_trace(simu, 'syn-proc-initial', simu_type = 'syn', suffix='_proc', src_space=1, trace_space=5, scale=0.8, color='r')
     elif optim.iter == optim.maxiter:
         data_misfit = np.loadtxt('./outputs/misfit_data.dat')
         data_misfit = data_misfit / data_misfit[0]
         plot_misfit(simu, data_misfit, 'data')
-        plot_trace(simu, 'syn-proc-final-model', simu_type = 'syn', suffix='_proc', src_space=1, trace_space=5, scale=0.8, color='b')
+        plot_trace(simu, 'syn-proc-final', simu_type = 'syn', suffix='_proc', src_space=1, trace_space=5, scale=0.8, color='b')
     else:
         pass
 
@@ -277,6 +277,9 @@ def plot_rtm(simu, optim, inv_scheme):
 
     plot_model2D(simu, grad.reshape(nx, nz).T, -grad_caxis, grad_caxis, 'RTM-image', colormap = 'gray')
 
+    # save RTM image
+    savebinfloat32(simu.system.homepath + 'outputs/gradient/RTM-image.bin', inv_scheme['g_now'])
+
 
 
 def my_seismic_cmap():