diff --git a/.DS_Store b/.DS_Store index cbc2f7c..8966778 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/toolbox/base.py b/toolbox/base.py index 4fb2880..cdcb170 100644 --- a/toolbox/base.py +++ b/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') diff --git a/toolbox/fielddata.py b/toolbox/fielddata.py new file mode 100644 index 0000000..9a6fb18 --- /dev/null +++ b/toolbox/fielddata.py @@ -0,0 +1,221 @@ +############################################################################### +# +# SWIT: Seismic Waveform Inversion Toolbox +# +# by Haipeng Li at USTC, haipengl@mail.ustc.edu.cn +# +# 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 diff --git a/toolbox/inversion.py b/toolbox/inversion.py index 73e0061..eeea62b 100644 --- a/toolbox/inversion.py +++ b/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 diff --git a/toolbox/misfit.py b/toolbox/misfit.py index 403a36f..dd5bc66 100644 --- a/toolbox/misfit.py +++ b/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 diff --git a/toolbox/plot.py b/toolbox/plot.py index b7e96c5..4fb25f6 100644 --- a/toolbox/plot.py +++ b/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(): diff --git a/toolbox/postprocess.py b/toolbox/postprocess.py index 7e28aed..b34ef4e 100644 --- a/toolbox/postprocess.py +++ b/toolbox/postprocess.py @@ -14,7 +14,7 @@ import scipy.signal from plot import plot_model2D -from tools import array2vector, smooth2d, vector2array +from tools import array2vector, smooth2d, vector2array, smooth1d def grad_precond(simu, optim, grad, forw, back): @@ -22,7 +22,6 @@ def grad_precond(simu, optim, grad, forw, back): ''' nx = simu.model.nx nz = simu.model.nz - vp = simu.model.vp vpmax = optim.vpmax marine_or_land = optim.marine_or_land grad_mute = optim.grad_mute @@ -41,15 +40,6 @@ def grad_precond(simu, optim, grad, forw, back): if grad_mute > 0: grad *= grad_taper(nx, nz, tapersize = grad_mute, thred = grad_thred, marine_or_land=marine_or_land) - if np.any(grad_mask == None): - pass - else: - if np.shape(grad_mask) != np.shape(grad): - raise('Wrong size of grad mask: the size of the mask should be identical to the size of vp model') - else: - grad *= grad_mask - - #apply the inverse Hessian if min(nx, nz) > 40: # set 40 grids in default span = 40 @@ -64,17 +54,30 @@ def grad_precond(simu, optim, grad, forw, back): precond = forw + back precond = precond / np.max(precond) precond[precond < epsilon] = epsilon - grad = grad / np.power(precond, 2) + grad = grad / np.power(precond, 1) # smooth the gradient if grad_smooth > 0: # exclude water-layer if marine_or_land in ['Marine', 'Offshore']: grad[:,grad_mute:] = smooth2d(grad[:,grad_mute:], span=grad_smooth) + # land gradient smooth else: grad = smooth2d(grad, span=grad_smooth) + + if np.any(grad_mask == None): + pass + else: + if np.shape(grad_mask) != np.shape(grad): + raise('Wrong size of grad mask: the size of the mask should be identical to the size of vp model') + else: + # apply mask + grad = grad * grad_mask + if optim.iter == 1: + plot_model2D(simu, grad_mask.T, 0., 1.0, 'grad_mask', 'jet') + # scale the gradient properly grad *= vpmax / abs(grad).max() diff --git a/toolbox/preprocess.py b/toolbox/preprocess.py index 3281262..5bdeba8 100644 --- a/toolbox/preprocess.py +++ b/toolbox/preprocess.py @@ -17,7 +17,7 @@ import obspy from scipy.signal import butter, filtfilt -from tools import add_su_header, get_su_parameter, loadsu, savesu, su2array +from tools import add_su_header, get_su_parameter, loadsu, savesu, su2array, array2su, smooth1d def process_workflow(simu, optim, simu_type = 'syn', use_first_break_in_su_header = False): @@ -99,7 +99,7 @@ def apply_mute(optim, trace, use_first_break_in_su_header): ''' apply time window and offset window mute. ''' # set parameters - length = 200 + length = 100 recn, nt, dt = get_su_parameter(trace) mute_late_arrival = optim.mute_late_arrival @@ -118,18 +118,28 @@ def apply_mute(optim, trace, use_first_break_in_su_header): if mute_offset_long: mute_offset(tr, mute_offset_long_dis, 'long') - # pick the first arrival and mute the late arrival + # pick the first arrival and mute the late or early arrival if mute_late_arrival: + + # decide mute late or early arrivals W.R.T first break + if mute_late_window > 0.0: + mutetype = 'late' + else: + mutetype = 'early' + + # mute window in second + mute_window = abs(mute_late_window) + # use the firsy arrival in the header if use_first_break_in_su_header: itrace = 0 pick = np.zeros(len(trace)) for tr in trace: - pick[itrace] = int((tr.stats.su.trace_header.mute_time_start_time_in_ms/1000 + mute_late_window)/dt) + pick[itrace] = int((tr.stats.su.trace_header.mute_time_start_time_in_ms/1000 + mute_window)/dt) itrace +=1 # pick the first arrival in a brutal manner else: - pick = brutal_picker(trace) + np.ceil(mute_late_window/dt) + pick = brutal_picker(trace) + np.ceil(mute_window/dt) itrace = 0 for tr in trace: @@ -137,7 +147,7 @@ def apply_mute(optim, trace, use_first_break_in_su_header): itmin = int(pick[itrace] - length/2) itmax = int(itmin + length) # apply mute - mute_arrival(tr, itmin, itmax, 'late', nt, length) + mute_arrival(tr, itmin, itmax, mutetype, nt, length) itrace +=1 @@ -295,3 +305,35 @@ def butter_highpass_filter(data, lowcut, dt, order=4): return data_hp + +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[:], window_len = 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 \ No newline at end of file diff --git a/toolbox/solver.py b/toolbox/solver.py index 564bcd0..3bf6eff 100644 --- a/toolbox/solver.py +++ b/toolbox/solver.py @@ -87,6 +87,7 @@ def adjoint(simu, optim): forw = loadbinfloat32(homepath+'data/syn/src0_illum_forw.bin') back = loadbinfloat32(homepath+'data/syn/src0_illum_back.bin') + # gradient precondtioning return grad_precond(simu, optim, grad, forw, back) diff --git a/toolbox/test.py b/toolbox/test.py new file mode 100644 index 0000000..999d76b --- /dev/null +++ b/toolbox/test.py @@ -0,0 +1,30 @@ +import sys +import matplotlib.pyplot as plt + +sys.path.append("/home/haipeng/PowerEdge4/data/SWIT-1.0-Open/toolbox") + +# import modules +import numpy as np +import base +from inversion import inversion, source_inversion +from preprocess import process_workflow +from solver import forward, source_wavelet +from tools import saveparjson, smooth2d, loadbinfloat32 + +vp_true = np.loadtxt('/home/haipeng/PowerEdge4/data/SWIT-1.0/thesis/Case13-Foothill/model/Foothill_801_291_25m_smooth25.dat') +grad = loadbinfloat32('/home/haipeng/PowerEdge4/data/SWIT-1.0/thesis/Case13-Foothill/outputs/gradient/grad-50.bin').reshape(801, 291) + +grad_mask = np.ones_like(vp_true) +grad_mask[vp_true==1500] = 0.0 + +grad_mask_new = np.zeros_like(grad_mask) +window_len = 9 +for i in range(801): + grad_mask_new[i, :] = smooth1d(grad_mask[i,:], window_len = window_len,window='hanning')[:-window_len+1] + +grad = grad * grad_mask * grad_mask_new + + +np.savetxt('/home/haipeng/Desktop/grad_new.dat', grad) +np.savetxt('/home/haipeng/Desktop/grad_mask.dat', grad_mask) +np.savetxt('/home/haipeng/Desktop/grad_mask_new.dat', grad_mask_new) \ No newline at end of file diff --git a/toolbox/tools.py b/toolbox/tools.py index b99f016..a64df58 100644 --- a/toolbox/tools.py +++ b/toolbox/tools.py @@ -276,6 +276,35 @@ def smooth2d(Z, span=10): return Z +def smooth1d(x,window_len=11,window='hanning'): + ''' smooth the data using a window with requested size. + ''' + + if x.ndim != 1: + raise ValueError("smooth only accepts 1 dimension arrays.") + + if x.size < window_len: + raise ValueError("Input vector needs to be bigger than window size.") + + if window_len<3: + return x + + if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']: + raise ValueError("Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'") + + s=np.r_[x[window_len-1:0:-1],x,x[-2:-window_len-1:-1]] + + #print(len(s)) + if window == 'flat': #moving average + w=np.ones(window_len,'d') + else: + w=eval('np.'+window+'(window_len)') + + y=np.convolve(w/w.sum(),s,mode='valid') + + return y + + def savetxt(filename, data, datatype): ''' save txt file '''