diff --git a/hypo2.py b/hypo2.py index 3d3f5c2..9d33839 100755 --- a/hypo2.py +++ b/hypo2.py @@ -15,12 +15,6 @@ from ttcrpy.rgrid import Grid3d -try: - import vtk -except ModuleNotFoundError: - print('VTK module not found, saving raypaths is disabled') - - # %% hypoloc def hypoloc(data, rcv, V, hinit, maxit, convh, verbose=False): """ @@ -148,8 +142,9 @@ def hypolocPS(data, rcv, V, hinit, maxit, convh, verbose=False): loc = hinit.copy() res = np.zeros((evID.size, maxit)) nev = 0 - + # set origin time to 0 and offset data accordingly + data = data.copy() hyp0_dt = {} for n in range(loc.shape[0]): hyp0_dt[loc[n, 0]] = loc[n, 1] @@ -307,6 +302,7 @@ def __init__(self, maxit, maxit_hypo, conv_hypo, Vlim, dmax, lagrangians, self._final_iteration = False # %% joint hypocenter - velocity + def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), Vpts=np.array([])): """ @@ -316,7 +312,7 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), ---------- par : instance of InvParams grid : instance of Grid3d - data : a numpy array with 3 columns + data : a numpy array with 5 columns 1st column is event ID number 2nd column is arrival time 3rd column is receiver index @@ -369,7 +365,7 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), sc = np.zeros(nsta) hyp0 = hinit.copy() - nnodes = grid.get_number_of_nodes() + nslowness = grid.nparams rcv_data = np.empty((data.shape[0], 3)) for ne in np.arange(nev): @@ -388,7 +384,7 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), hcal = np.column_stack((caldata[:, 0], np.zeros(caldata.shape[0]), caldata[:, 3:6])) tcal = caldata[:, 1] rcv_cal = np.empty((caldata.shape[0], 3)) - Msc_cal = [] + Lsc_cal = [] for nc in range(ncal): indr = np.nonzero(caldata[:, 0] == calID[nc])[0] nst = np.sum(indr.size) @@ -398,14 +394,14 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), tmp = np.zeros((nst, nsta)) for n in range(nst): tmp[n, int(1.e-6+caldata[indr[n], 2])] = 1.0 - Msc_cal.append(sp.csr_matrix(tmp)) + Lsc_cal.append(sp.csr_matrix(tmp)) else: ncal = 0 tcal = np.array([]) if np.isscalar(Vinit): - V = Vinit + np.zeros(nnodes) - s = np.ones(nnodes)/Vinit + V = Vinit + np.zeros(nslowness) + s = np.ones(nslowness)/Vinit else: V = Vinit s = 1./Vinit @@ -422,25 +418,28 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), resV = np.zeros(par.maxit+1) resAxb = np.zeros(par.maxit) - P = np.ones(nnodes) - dP = sp.csr_matrix((np.ones(nnodes), (np.arange(nnodes, dtype=np.int64), - np.arange(nnodes, dtype=np.int64))), shape=(nnodes, nnodes)) + P = np.ones(nslowness) + dP = sp.csr_matrix((np.ones(nslowness), (np.arange(nslowness, dtype=np.int64), + np.arange(nslowness, dtype=np.int64))), shape=(nslowness, nslowness)) + + Spmax = 1. / par.Vpmin + Spmin = 1. / par.Vpmax - deltam = np.zeros(nnodes+nsta) + deltam = np.zeros(nslowness+nsta) if par.constr_sc and par.use_sc: - u1 = np.ones(nnodes+nsta) - u1[:nnodes] = 0.0 + u1 = np.ones(nslowness+nsta) + u1[:nslowness] = 0.0 - i = np.arange(nnodes, nnodes+nsta) - j = np.arange(nnodes, nnodes+nsta) + i = np.arange(nslowness, nslowness+nsta) + j = np.arange(nslowness, nslowness+nsta) nn = i.size i = np.kron(i, np.ones((nn,))) j = np.kron(np.ones((nn,)), j) u1Tu1 = sp.csr_matrix((np.ones((i.size,)), (i, j)), - shape=(nnodes+nsta, nnodes+nsta)) + shape=(nslowness+nsta, nslowness+nsta)) else: - u1 = np.zeros(nnodes+nsta) - u1Tu1 = sp.csr_matrix((nnodes+nsta, nnodes+nsta)) + u1 = np.zeros(nslowness+nsta) + u1Tu1 = sp.csr_matrix((nslowness+nsta, nslowness+nsta)) if Vpts.size > 0: if par.verbose: @@ -448,6 +447,7 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), sys.stdout.flush() D = grid.compute_D(Vpts[:, 1:]) D1 = sp.hstack((D, sp.coo_matrix((Vpts.shape[0], nsta)))).tocsr() + Spts = 1. / Vpts[:, 0] else: D = 0.0 @@ -455,11 +455,11 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), print('Building regularization matrix K') sys.stdout.flush() Kx, Ky, Kz = grid.compute_K() - Kx1 = sp.hstack((Kx, sp.coo_matrix((nnodes, nsta)))).tocsr() + Kx1 = sp.hstack((Kx, sp.coo_matrix((nslowness, nsta)))).tocsr() KtK = Kx1.T * Kx1 - Ky1 = sp.hstack((Ky, sp.coo_matrix((nnodes, nsta)))).tocsr() + Ky1 = sp.hstack((Ky, sp.coo_matrix((nslowness, nsta)))).tocsr() KtK += Ky1.T * Ky1 - Kz1 = sp.hstack((Kz, sp.coo_matrix((nnodes, nsta)))).tocsr() + Kz1 = sp.hstack((Kz, sp.coo_matrix((nslowness, nsta)))).tocsr() KtK += par.wzK * Kz1.T * Kz1 nK = spl.norm(KtK) Kx = Kx.tocsr() @@ -483,22 +483,21 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), sys.stdout.flush() # compute vector C - cx = Kx * V - cy = Ky * V - cz = Kz * V + cx = Kx * s + cy = Ky * s + cz = Kz * s # compute dP/dV, matrix of penalties derivatives - for n in np.arange(nnodes): - if V[n] < par.Vpmin: - P[n] = par.PAp * (par.Vpmin-V[n]) + for n in np.arange(nslowness): + if s[n] < Spmin: + P[n] = par.PAp * (Spmin-s[n]) dP[n, n] = -par.PAp - elif V[n] > par.Vpmax: - P[n] = par.PAp * (V[n]-par.Vpmax) + elif s[n] > Spmax: + P[n] = par.PAp * (s[n]-Spmax) dP[n, n] = par.PAp else: P[n] = 0.0 dP[n, n] = 0.0 - if par.verbose: npel = np.sum(P != 0.0) if npel > 0: @@ -515,16 +514,16 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), indr = np.nonzero(data[:, 0] == evID[ne])[0] for i in indr: hyp[i, :] = hyp0[indh[0], :] - tcalc, rays, Mev = grid.raytrace(hyp, rcv_data, s, + tcalc, rays, Lev = grid.raytrace(hyp, rcv_data, s, return_rays=True, - compute_M=True) + compute_L=True) s0 = grid.get_s0(hyp) else: tcalc = np.array([]) if ncal > 0: - tcalc_cal, Mcal = grid.raytrace(hcal, rcv_cal, s, - compute_M=True) + tcalc_cal, Lcal = grid.raytrace(hcal, rcv_cal, s, + compute_L=True) else: tcalc_cal = np.array([]) @@ -545,10 +544,10 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), # initializing matrix M; matrix of partial derivatives of velocity dt/dV if par.verbose: - print(' Building matrix M') + print(' Building matrix L') sys.stdout.flush() - M1 = None + L1 = None ir1 = 0 for ne in range(nev): if par.verbose: @@ -573,45 +572,45 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), Q, _ = np.linalg.qr(H, mode='complete') T = sp.csr_matrix(Q[:, 4:]).T - M = sp.csr_matrix(Mev[ne], shape=(nst, nnodes)) + L = sp.csr_matrix(Lev[ne], shape=(nst, nslowness)) if par.use_sc: - Msc = np.zeros((nst, nsta)) + Lsc = np.zeros((nst, nsta)) for ns in range(nst): - Msc[ns, int(1.e-6+data[indr[ns], 2])] = 1. - M = sp.hstack((M, sp.csr_matrix(Msc))) + Lsc[ns, int(1.e-6+data[indr[ns], 2])] = 1. + L = sp.hstack((L, sp.csr_matrix(Lsc))) - M = T * M + L = T * L - if M1 is None: - M1 = M + if L1 is None: + L1 = L else: - M1 = sp.vstack((M1, M)) + L1 = sp.vstack((L1, L)) r1[ir1 + np.arange(nst2, dtype=np.int64)] = T.dot(r1a[indr]) ir1 += nst2 for nc in range(ncal): - M = Mcal[nc] + L = Lcal[nc] if par.use_sc: - M = sp.hstack((M, Msc_cal[nc])) - if M1 is None: - M1 = M + L = sp.hstack((L, Lsc_cal[nc])) + if L1 is None: + L1 = L else: - M1 = sp.vstack((M1, M)) + L1 = sp.vstack((L1, L)) if par.verbose: print(' Assembling matrices and solving system') sys.stdout.flush() - s = -np.sum(sc) # u1.T * deltam + ssc = -np.sum(sc) # u1.T * deltam - dP1 = sp.hstack((dP, sp.csr_matrix(np.zeros((nnodes, nsta))))).tocsr() # dP prime + dP1 = sp.hstack((dP, sp.csr_matrix(np.zeros((nslowness, nsta))))).tocsr() # dP prime # compute A & h for inversion - M1 = M1.tocsr() + L1 = L1.tocsr() - A = M1.T * M1 + A = L1.T * L1 nM = spl.norm(A) λ = par.λ * nM / nK @@ -628,12 +627,12 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), A += γ*tmp A += u1Tu1 - b = M1.T * r1 + b = L1.T * r1 tmp2x = Kx1.T * cx tmp2y = Ky1.T * cy tmp2z = Kz1.T * cz tmp3 = dP1.T * P - tmp = u1 * s + tmp = u1 * ssc b += - λ*tmp2x - λ*tmp2y - par.wzK*λ*tmp2z - γ*tmp3 - tmp if Vpts.shape[0] > 0: @@ -641,8 +640,7 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), nD = spl.norm(tmp) α = par.α * nM / nD A += α * tmp - b += α * D1.T * (Vpts[:, 0] - D*V) - + b += α * D1.T * (Spts - D*s) if par.verbose: print(' calling minres with system of size {0:d} x {1:d}'.format(A.shape[0], A.shape[1])) @@ -652,15 +650,15 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), deltam = x[0] resAxb[it] = np.linalg.norm(A*deltam - b) - dmean = np.mean(np.abs(deltam[:nnodes])) + dmean = np.mean(np.abs(deltam[:nslowness])) if dmean > par.dVp_max: if par.verbose: - print(' Scaling Vp perturbations by {0:e}'.format(par.dVp_max/dmean)) - deltam[:nnodes] = deltam[:nnodes] * par.dVp_max/dmean + print(' Scaling Slowness perturbations by {0:e}'.format(1./(par.dVp_max*dmean))) + deltam[:nslowness] = deltam[:nslowness] / (par.dVp_max*dmean) - V += deltam[:nnodes] - s = 1. / V - sc += deltam[nnodes:] + s += deltam[:nslowness] + V = 1. / s + sc += deltam[nslowness:] if par.save_V: if par.verbose: @@ -745,1263 +743,164 @@ def jointHypoVel(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), return hyp0, V, sc, (resV, resAxb) -def jointHypoVel_c(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), - Vpts=np.array([])): - """ - Joint hypocenter-velocity inversion on a regular grid +def _rl_worker(thread_no, istart, iend, par, grid, evID, hyp0, data, rcv, tobs, h_queue): + for ne in range(istart, iend): + h, indh = _reloc(ne, par, grid, evID, hyp0, data, rcv, tobs, thread_no) + h_queue.put((h, indh)) + h_queue.close() - Parameters - ---------- - par : instance of InvParams - grid : instance of Grid3d - data : a numpy array with 5 columns - 1st column is event ID number - 2nd column is arrival time - 3rd column is receiver index - *** important *** - data should sorted by event ID first, then by receiver index - rcv: : coordinates of receivers - 1st column is receiver easting - 2nd column is receiver northing - 3rd column is receiver elevation - Vinit : initial velocity model - hinit : initial hypocenter coordinate - 1st column is event ID number - 2nd column is origin time - 3rd column is event easting - 4th column is event northing - 5th column is event elevation - *** important *** - for efficiency reason when computing matrix M, initial hypocenters - should _not_ be equal for any two event, e.g. they shoud all be - different - caldata : calibration shot data, numpy array with 8 columns - 1st column is cal shot ID number - 2nd column is arrival time - 3rd column is receiver index - 4th column is source easting - 5th column is source northing - 6th column is source elevation - *** important *** - cal shot data should sorted by cal shot ID first, then by receiver index - Vpts : known velocity points, numpy array with 4 columns - 1st column is velocity - 2nd column is easting - 3rd column is northing - 4th column is elevation - Returns - ------- - loc : hypocenter coordinates - V : velocity model - sc : static corrections - res : residuals - """ +def _reloc(ne, par, grid, evID, hyp0, data, rcv, tobs, thread_no=None): - evID = np.unique(data[:, 0]) - nev = evID.size - if par.use_sc: - nsta = rcv.shape[0] - else: - nsta = 0 + if par.verbose: + print(' Updating event ID {0:d} ({1:d}/{2:d})'.format(int(1.e-6+evID[ne]), ne+1, evID.size)) + sys.stdout.flush() - sc = np.zeros(nsta) - hyp0 = hinit.copy() - nslowness = grid.nparams + indh = np.nonzero(hyp0[:, 0] == evID[ne])[0][0] + indr = np.nonzero(data[:, 0] == evID[ne])[0] - rcv_data = np.empty((data.shape[0], 3)) - for ne in np.arange(nev): - indr = np.nonzero(data[:, 0] == evID[ne])[0] - for i in indr: - rcv_data[i, :] = rcv[int(1.e-6+data[i, 2])] + hyp_save = hyp0[indh, :].copy() - if data.shape[0] > 0: - tobs = data[:, 1] - else: - tobs = np.array([]) + nst = np.sum(indr.size) - if caldata.shape[0] > 0: - calID = np.unique(caldata[:, 0]) - ncal = calID.size - hcal = np.column_stack((caldata[:, 0], np.zeros(caldata.shape[0]), caldata[:, 3:6])) - tcal = caldata[:, 1] - rcv_cal = np.empty((caldata.shape[0], 3)) - Lsc_cal = [] - for nc in range(ncal): - indr = np.nonzero(caldata[:, 0] == calID[nc])[0] - nst = np.sum(indr.size) - for i in indr: - rcv_cal[i, :] = rcv[int(1.e-6+caldata[i, 2])] - if par.use_sc: - tmp = np.zeros((nst, nsta)) - for n in range(nst): - tmp[n, int(1.e-6+caldata[indr[n], 2])] = 1.0 - Lsc_cal.append(sp.csr_matrix(tmp)) - else: - ncal = 0 - tcal = np.array([]) + hyp = np.empty((nst, 5)) + stn = np.empty((nst, 3)) + for i in range(nst): + hyp[i, :] = hyp0[indh, :] + stn[i, :] = rcv[int(1.e-6+data[indr[i], 2]), :] - if np.isscalar(Vinit): - V = Vinit + np.zeros(nslowness) - s = np.ones(nslowness)/Vinit - else: - V = Vinit - s = 1./Vinit + if par.hypo_2step: + if par.verbose: + print(' Updating latitude & longitude', end='') + sys.stdout.flush() + H = np.ones((nst, 2)) + for itt in range(par.maxit_hypo): + for i in range(nst): + hyp[i, :] = hyp0[indh, :] - if Vpts.size > 0: - if Vpts.shape[1] > 4: # check if we have Vs data in array - itmp = Vpts[:, 4] == 0 - Vpts = Vpts[itmp, :4] # keep only Vp data + tcalc, rays = grid.raytrace(hyp, stn, thread_no=thread_no, + return_rays=True) + s0 = grid.get_s0(hyp) + for ns in range(nst): + raysi = rays[ns] + S0 = s0[ns] - if par.verbose: - print('\n *** Joint hypocenter-velocity inversion ***\n') + d = (raysi[1, :]-hyp0[indh, 2:]).flatten() + ds = np.sqrt(np.sum(d*d)) + H[ns, 0] = -S0 * d[0]/ds + H[ns, 1] = -S0 * d[1]/ds - if par.invert_vel: - resV = np.zeros(par.maxit+1) - resAxb = np.zeros(par.maxit) + r = tobs[indr] - tcalc + x = lstsq(H, r) + deltah = x[0] - P = np.ones(nslowness) - dP = sp.csr_matrix((np.ones(nslowness), (np.arange(nslowness, dtype=np.int64), - np.arange(nslowness, dtype=np.int64))), shape=(nslowness, nslowness)) + if np.sum(np.isfinite(deltah)) != deltah.size: + try: + U,S,VVh = np.linalg.svd(H.T.dot(H)+1e-9*np.eye(2)) + VV = VVh.T + deltah = np.dot(VV, np.dot(U.T, H.T.dot(r))/S) + except np.linalg.linalg.LinAlgError: + print(' - Event could not be relocated, resetting and exiting') + hyp0[indh, :] = hyp_save + return hyp_save, indh - Spmax = 1. / par.Vpmin - Spmin = 1. / par.Vpmax + for n in range(2): + if np.abs(deltah[n]) > par.dx_max: + deltah[n] = par.dx_max * np.sign(deltah[n]) - deltam = np.zeros(nslowness+nsta) - if par.constr_sc and par.use_sc: - u1 = np.ones(nslowness+nsta) - u1[:nslowness] = 0.0 + new_hyp = hyp0[indh, :].copy() + new_hyp[2:4] += deltah + if grid.is_outside(new_hyp[2:].reshape((1, 3))): + print(' Event could not be relocated inside the grid ({0:f}, {1:f}, {2:f}), resetting and exiting'.format(new_hyp[2], new_hyp[3], new_hyp[4])) + hyp0[indh, :] = hyp_save + return hyp_save, indh - i = np.arange(nslowness, nslowness+nsta) - j = np.arange(nslowness, nslowness+nsta) - nn = i.size - i = np.kron(i, np.ones((nn,))) - j = np.kron(np.ones((nn,)), j) - u1Tu1 = sp.csr_matrix((np.ones((i.size,)), (i, j)), - shape=(nslowness+nsta, nslowness+nsta)) - else: - u1 = np.zeros(nslowness+nsta) - u1Tu1 = sp.csr_matrix((nslowness+nsta, nslowness+nsta)) + hyp0[indh, 2:4] += deltah - if Vpts.size > 0: - if par.verbose: - print('Building velocity data point matrix D') - sys.stdout.flush() - D = grid.compute_D(Vpts[:, 1:]) - D1 = sp.hstack((D, sp.coo_matrix((Vpts.shape[0], nsta)))).tocsr() - Spts = 1. / Vpts[:, 0] - else: - D = 0.0 - - if par.verbose: - print('Building regularization matrix K') - sys.stdout.flush() - Kx, Ky, Kz = grid.compute_K() - Kx1 = sp.hstack((Kx, sp.coo_matrix((nslowness, nsta)))).tocsr() - KtK = Kx1.T * Kx1 - Ky1 = sp.hstack((Ky, sp.coo_matrix((nslowness, nsta)))).tocsr() - KtK += Ky1.T * Ky1 - Kz1 = sp.hstack((Kz, sp.coo_matrix((nslowness, nsta)))).tocsr() - KtK += par.wzK * Kz1.T * Kz1 - nK = spl.norm(KtK) - Kx = Kx.tocsr() - Ky = Ky.tocsr() - Kz = Kz.tocsr() - else: - resV = None - resAxb = None - - if par.verbose: - print('\nStarting iterations') - - for it in np.arange(par.maxit): - - par._final_iteration = it == par.maxit-1 - - if par.invert_vel: - if par.verbose: - print('\nIteration {0:d} - Updating velocity model'.format(it+1)) - print(' Updating penalty vector') - sys.stdout.flush() - - # compute vector C - cx = Kx * s - cy = Ky * s - cz = Kz * s - - # compute dP/dV, matrix of penalties derivatives - for n in np.arange(nslowness): - if s[n] < Spmin: - P[n] = par.PAp * (Spmin-s[n]) - dP[n, n] = -par.PAp - elif s[n] > Spmax: - P[n] = par.PAp * (s[n]-Spmax) - dP[n, n] = par.PAp - else: - P[n] = 0.0 - dP[n, n] = 0.0 - if par.verbose: - npel = np.sum(P != 0.0) - if npel > 0: - print(' Penalties applied at {0:d} nodes'.format(npel)) - - if par.verbose: - print(' Raytracing') - sys.stdout.flush() - - if nev > 0: - hyp = np.empty((data.shape[0], 5)) - for ne in np.arange(nev): - indh = np.nonzero(hyp0[:, 0] == evID[ne])[0] - indr = np.nonzero(data[:, 0] == evID[ne])[0] - for i in indr: - hyp[i, :] = hyp0[indh[0], :] - tcalc, rays, Lev = grid.raytrace(hyp, rcv_data, s, - return_rays=True, - compute_L=True) - s0 = grid.get_s0(hyp) - else: - tcalc = np.array([]) - - if ncal > 0: - tcalc_cal, Lcal = grid.raytrace(hcal, rcv_cal, s, - compute_L=True) - else: - tcalc_cal = np.array([]) - - r1a = tobs - tcalc - r1 = tcal - tcalc_cal - if r1a.size > 0: - r1 = np.hstack((np.zeros(data.shape[0]-4*nev), r1)) - - if par.show_plots: - plt.figure(1) - plt.cla() - plt.plot(r1a, 'o') - plt.title('Residuals - Iteration {0:d}'.format(it+1)) - plt.show(block=False) - plt.pause(0.0001) - - resV[it] = np.linalg.norm(np.hstack((r1a, r1))) - - # initializing matrix M; matrix of partial derivatives of velocity dt/dV - if par.verbose: - print(' Building matrix L') - sys.stdout.flush() - - L1 = None - ir1 = 0 - for ne in range(nev): - if par.verbose: - print(' Event ID '+str(int(1.e-6+evID[ne]))) - sys.stdout.flush() - - indh = np.nonzero(hyp0[:, 0] == evID[ne])[0] - indr = np.nonzero(data[:, 0] == evID[ne])[0] - - nst = np.sum(indr.size) - nst2 = nst-4 - H = np.ones((nst, 4)) - for ns in range(nst): - raysi = rays[indr[ns]] - S0 = s0[indr[ns]] - - d = (raysi[1, :]-hyp0[indh, 2:]).flatten() - ds = np.sqrt(np.sum(d*d)) - H[ns, 1] = -S0 * d[0]/ds - H[ns, 2] = -S0 * d[1]/ds - H[ns, 3] = -S0 * d[2]/ds - - Q, _ = np.linalg.qr(H, mode='complete') - T = sp.csr_matrix(Q[:, 4:]).T - L = sp.csr_matrix(Lev[ne], shape=(nst, nslowness)) - if par.use_sc: - Lsc = np.zeros((nst, nsta)) - for ns in range(nst): - Lsc[ns, int(1.e-6+data[indr[ns], 2])] = 1. - L = sp.hstack((L, sp.csr_matrix(Lsc))) - - L = T * L - - if L1 is None: - L1 = L - else: - L1 = sp.vstack((L1, L)) - - r1[ir1 + np.arange(nst2, dtype=np.int64)] = T.dot(r1a[indr]) - ir1 += nst2 - - for nc in range(ncal): - L = Lcal[nc] - if par.use_sc: - L = sp.hstack((L, Lsc_cal[nc])) - if L1 is None: - L1 = L - else: - L1 = sp.vstack((L1, L)) - - if par.verbose: - print(' Assembling matrices and solving system') - sys.stdout.flush() - - ssc = -np.sum(sc) # u1.T * deltam - - dP1 = sp.hstack((dP, sp.csr_matrix(np.zeros((nslowness, nsta))))).tocsr() # dP prime - - # compute A & h for inversion - - L1 = L1.tocsr() - - A = L1.T * L1 - nM = spl.norm(A) - - λ = par.λ * nM / nK - - A += λ*KtK - - tmp = dP1.T * dP1 - nP = spl.norm(tmp) - if nP != 0.0: - γ = par.γ * nM / nP - else: - γ = par.γ - - A += γ*tmp - A += u1Tu1 - - b = L1.T * r1 - tmp2x = Kx1.T * cx - tmp2y = Ky1.T * cy - tmp2z = Kz1.T * cz - tmp3 = dP1.T * P - tmp = u1 * ssc - b += - λ*tmp2x - λ*tmp2y - par.wzK*λ*tmp2z - γ*tmp3 - tmp - - if Vpts.shape[0] > 0: - tmp = D1.T * D1 - nD = spl.norm(tmp) - α = par.α * nM / nD - A += α * tmp - b += α * D1.T * (Spts - D*s) - - if par.verbose: - print(' calling minres with system of size {0:d} x {1:d}'.format(A.shape[0], A.shape[1])) - sys.stdout.flush() - x = spl.minres(A, b) - - deltam = x[0] - resAxb[it] = np.linalg.norm(A*deltam - b) - - dmean = np.mean(np.abs(deltam[:nslowness])) - if dmean > par.dVp_max: - if par.verbose: - print(' Scaling Slowness perturbations by {0:e}'.format(1./(par.dVp_max*dmean))) - deltam[:nslowness] = deltam[:nslowness] / (par.dVp_max*dmean) - - s += deltam[:nslowness] - V = 1. / s - sc += deltam[nslowness:] - - if par.save_V: - if par.verbose: - print(' Saving Velocity model') - if 'vtk' in sys.modules: - grid.to_vtk({'Vp': V}, 'Vp{0:02d}'.format(it+1)) - - grid.set_slowness(s) - - if nev > 0: - if par.verbose: - print('Iteration {0:d} - Relocating events'.format(it+1)) - sys.stdout.flush() - - if grid.nthreads == 1 or nev < grid.nthreads: - for ne in range(nev): - _reloc(ne, par, grid, evID, hyp0, data, rcv, tobs) - else: - # run in parallel - blk_size = np.zeros((grid.nthreads,), dtype=np.int64) - nj = nev - while nj > 0: - for n in range(grid.nthreads): - blk_size[n] += 1 - nj -= 1 - if nj == 0: - break - processes = [] - blk_start = 0 - h_queue = Queue() - for n in range(grid.nthreads): - blk_end = blk_start + blk_size[n] - p = Process(target=_rl_worker, - args=(n, blk_start, blk_end, par, grid, evID, hyp0, data, rcv, tobs, h_queue), - daemon=True) - processes.append(p) - p.start() - blk_start += blk_size[n] - - for ne in range(nev): - h, indh = h_queue.get() - hyp0[indh, :] = h - - if par.invert_vel: - if nev > 0: - hyp = np.empty((data.shape[0], 5)) - for ne in np.arange(nev): - indh = np.nonzero(hyp0[:, 0] == evID[ne])[0] - indr = np.nonzero(data[:, 0] == evID[ne])[0] - for i in indr: - hyp[i, :] = hyp0[indh[0], :] - tcalc = grid.raytrace(hyp, rcv_data, s) - else: - tcalc = np.array([]) - - if ncal > 0: - tcalc_cal = grid.raytrace(hcal, rcv_cal, s) - else: - tcalc_cal = np.array([]) - - r1a = tobs - tcalc - r1 = tcal - tcalc_cal - if r1a.size > 0: - r1 = np.hstack((np.zeros(data.shape[0]-4*nev), r1)) - - if par.show_plots: - plt.figure(1) - plt.cla() - plt.plot(r1a, 'o') - plt.title('Residuals - Final step') - plt.show(block=False) - plt.pause(0.0001) - - resV[-1] = np.linalg.norm(np.hstack((r1a, r1))) - - r1 = r1.reshape(-1, 1) - r1a = r1a.reshape(-1, 1) - - if par.verbose: - print('\n ** Inversion complete **\n', flush=True) - - return hyp0, V, sc, (resV, resAxb) - - -def _rl_worker(thread_no, istart, iend, par, grid, evID, hyp0, data, rcv, tobs, h_queue): - for ne in range(istart, iend): - h, indh = _reloc(ne, par, grid, evID, hyp0, data, rcv, tobs, thread_no) - h_queue.put((h, indh)) - h_queue.close() - - -def _reloc(ne, par, grid, evID, hyp0, data, rcv, tobs, thread_no=None): - - if par.verbose: - print(' Updating event ID {0:d} ({1:d}/{2:d})'.format(int(1.e-6+evID[ne]), ne+1, evID.size)) - sys.stdout.flush() - - indh = np.nonzero(hyp0[:, 0] == evID[ne])[0][0] - indr = np.nonzero(data[:, 0] == evID[ne])[0] - - hyp_save = hyp0[indh, :].copy() - - nst = np.sum(indr.size) - - hyp = np.empty((nst, 5)) - stn = np.empty((nst, 3)) - for i in range(nst): - hyp[i, :] = hyp0[indh, :] - stn[i, :] = rcv[int(1.e-6+data[indr[i], 2]), :] - - if par.hypo_2step: - if par.verbose: - print(' Updating latitude & longitude', end='') - sys.stdout.flush() - H = np.ones((nst, 2)) - for itt in range(par.maxit_hypo): - for i in range(nst): - hyp[i, :] = hyp0[indh, :] - - tcalc, rays = grid.raytrace(hyp, stn, thread_no=thread_no, - return_rays=True) - s0 = grid.get_s0(hyp) - for ns in range(nst): - raysi = rays[ns] - S0 = s0[ns] - - d = (raysi[1, :]-hyp0[indh, 2:]).flatten() - ds = np.sqrt(np.sum(d*d)) - H[ns, 0] = -S0 * d[0]/ds - H[ns, 1] = -S0 * d[1]/ds - - r = tobs[indr] - tcalc - x = lstsq(H, r) - deltah = x[0] - - if np.sum(np.isfinite(deltah)) != deltah.size: - try: - U,S,VVh = np.linalg.svd(H.T.dot(H)+1e-9*np.eye(2)) - VV = VVh.T - deltah = np.dot(VV, np.dot(U.T, H.T.dot(r))/S) - except np.linalg.linalg.LinAlgError: - print(' - Event could not be relocated, resetting and exiting') - hyp0[indh, :] = hyp_save - return hyp_save, indh - - for n in range(2): - if np.abs(deltah[n]) > par.dx_max: - deltah[n] = par.dx_max * np.sign(deltah[n]) - - new_hyp = hyp0[indh, :].copy() - new_hyp[2:4] += deltah - if grid.is_outside(new_hyp[2:].reshape((1, 3))): - print(' Event could not be relocated inside the grid ({0:f}, {1:f}, {2:f}), resetting and exiting'.format(new_hyp[2], new_hyp[3], new_hyp[4])) - hyp0[indh, :] = hyp_save - return hyp_save, indh - - hyp0[indh, 2:4] += deltah - - if np.sum(np.abs(deltah) < par.conv_hypo) == 2: - if par.verbose: - print(' - converged at iteration '+str(itt+1)) - sys.stdout.flush() - break - - else: - if par.verbose: - print(' - reached max number of iterations') - sys.stdout.flush() - - if par.verbose: - print(' Updating all hypocenter params', end='') - sys.stdout.flush() - - H = np.ones((nst, 4)) - for itt in range(par.maxit_hypo): - for i in range(nst): - hyp[i, :] = hyp0[indh, :] - - tcalc, rays = grid.raytrace(hyp, stn, thread_no=thread_no, - return_rays=True) - s0 = grid.get_s0(hyp) - for ns in range(nst): - raysi = rays[ns] - S0 = s0[ns] - - d = (raysi[1, :]-hyp0[indh, 2:]).flatten() - ds = np.sqrt(np.sum(d*d)) - H[ns, 1] = -S0 * d[0]/ds - H[ns, 2] = -S0 * d[1]/ds - H[ns, 3] = -S0 * d[2]/ds - - r = tobs[indr] - tcalc - x = lstsq(H, r) - deltah = x[0] - - if np.sum(np.isfinite(deltah)) != deltah.size: - try: - U,S,VVh = np.linalg.svd(H.T.dot(H)+1e-9*np.eye(4)) - VV = VVh.T - deltah = np.dot(VV, np.dot(U.T, H.T.dot(r))/S) - except np.linalg.linalg.LinAlgError: - print(' Event could not be relocated, resetting and exiting') - hyp0[indh, :] = hyp_save - return hyp_save, indh - - if np.abs(deltah[0]) > par.dt_max: - deltah[0] = par.dt_max * np.sign(deltah[0]) - for n in range(1, 4): - if np.abs(deltah[n]) > par.dx_max: - deltah[n] = par.dx_max * np.sign(deltah[n]) - - new_hyp = hyp0[indh, 1:] + deltah - if grid.is_outside(new_hyp[1:].reshape((1, 3))): - print(' Event could not be relocated inside the grid ({0:f}, {1:f}, {2:f}), resetting and exiting'.format(new_hyp[1], new_hyp[2], new_hyp[3])) - hyp0[indh, :] = hyp_save - return hyp_save, indh - - hyp0[indh, 1:] += deltah - - if np.sum(np.abs(deltah[1:]) < par.conv_hypo) == 3: - if par.verbose: - print(' - converged at iteration '+str(itt+1)) - sys.stdout.flush() - break - - else: - if par.verbose: - print(' - reached max number of iterations') - sys.stdout.flush() - - if par.save_rp and 'vtk' in sys.modules and par._final_iteration: - if par.verbose: - print(' Saving raypaths') - filename = 'raypaths' - key = 'ev_{0:d}'.format(int(1.e-6+evID[ne])) - grid.to_vtk({key: rays}, filename) - - return hyp0[indh, :], indh - - -def jointHypoVelPS(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), - Vpts=np.array([])): - """ - Joint hypocenter-velocity inversion on a regular grid - - Parameters - ---------- - par : instance of InvParams - grid : instances of Grid3d - data : a numpy array with 5 columns - 1st column is event ID number - 2nd column is arrival time - 3rd column is receiver index - 4th column is code for wave phase: 0 for P-wave and 1 for S-wave - *** important *** - data should sorted by event ID first, then by receiver index - rcv: : coordinates of receivers - 1st column is receiver easting - 2nd column is receiver northing - 3rd column is receiver elevation - Vinit : tuple with initial velocity model (P-wave first, S-wave second) - hinit : initial hypocenter coordinate - 1st column is event ID number - 2nd column is origin time - 3rd column is event easting - 4th column is event northing - 5th column is event elevation - *** important *** - for efficiency reason when computing matrix M, initial hypocenters - should _not_ be equal for any two event, e.g. they shoud all be - different - caldata : calibration shot data, numpy array with 8 columns - 1st column is cal shot ID number - 2nd column is arrival time - 3rd column is receiver index - 4th column is source easting - 5th column is source northing - 6th column is source elevation - 7th column is code for wave phase: 0 for P-wave and 1 for S-wave - *** important *** - cal shot data should sorted by cal shot ID first, then by receiver index - Vpts : known velocity points, numpy array with 4 columns - 1st column is velocity - 2nd column is easting - 3rd column is northing - 4th column is elevation - 5th column is code for wave phase: 0 for P-wave and 1 for S-wave - - Returns - ------- - loc : hypocenter coordinates - V : velocity models (tuple holding Vp and Vs) - sc : static corrections - res : residuals - """ - - if grid.nthreads > 1: - # we need a second instance for parallel computations - grid_s = Grid3d(grid.x, grid.y, grid.z, grid.nthreads, cell_slowness=False) - else: - grid_s = grid - - evID = np.unique(data[:, 0]) - nev = evID.size - if par.use_sc: - nsta = rcv.shape[0] - else: - nsta = 0 - - sc_p = np.zeros(nsta) - sc_s = np.zeros(nsta) - hyp0 = hinit.copy() - nnodes = grid.get_number_of_nodes() - - # sort data by seismic phase (P-wave first S-wave second) - indp = data[:, 3] == 0.0 - inds = data[:, 3] == 1.0 - nttp = np.sum(indp) - ntts = np.sum(inds) - datap = data[indp, :] - datas = data[inds, :] - data = np.vstack((datap, datas)) - # update indices - indp = data[:, 3] == 0.0 - inds = data[:, 3] == 1.0 - - rcv_datap = np.empty((nttp, 3)) - rcv_datas = np.empty((ntts, 3)) - for ne in np.arange(nev): - indr = np.nonzero(np.logical_and(data[:, 0] == evID[ne], indp))[0] - for i in indr: - rcv_datap[i, :] = rcv[int(1.e-6+data[i, 2])] - indr = np.nonzero(np.logical_and(data[:, 0] == evID[ne], inds))[0] - for i in indr: - rcv_datas[i-nttp, :] = rcv[int(1.e-6+data[i, 2])] - - - if data.shape[0] > 0: - tobs = data[:, 1] - else: - tobs = np.array([]) - - if caldata.shape[0] > 0: - calID = np.unique(caldata[:, 0]) - ncal = calID.size - - # sort data by seismic phase (P-wave first S-wave second) - indcalp = caldata[:, 6] == 0.0 - indcals = caldata[:, 6] == 1.0 - nttcalp = np.sum(indcalp) - nttcals = np.sum(indcals) - caldatap = caldata[indcalp, :] - caldatas = caldata[indcals, :] - caldata = np.vstack((caldatap, caldatas)) - # update indices - indcalp = caldata[:, 6] == 0.0 - indcals = caldata[:, 6] == 1.0 - - hcalp = np.column_stack((caldata[indcalp, 0], np.zeros(nttcalp), caldata[indcalp, 3:6])) - hcals = np.column_stack((caldata[indcals, 0], np.zeros(nttcals), caldata[indcals, 3:6])) - - rcv_calp = np.empty((nttcalp, 3)) - rcv_cals = np.empty((nttcals, 3)) - - tcal = caldata[:, 1] - Msc_cal = [] - for nc in range(ncal): - indrp = np.nonzero(np.logical_and(caldata[:, 0] == calID[nc], indcalp))[0] - for i in indrp: - rcv_calp[i, :] = rcv[int(1.e-6+caldata[i, 2])] - indrs = np.nonzero(np.logical_and(caldata[:, 0] == calID[nc], indcals))[0] - for i in indrs: - rcv_cals[i-nttcalp, :] = rcv[int(1.e-6+caldata[i, 2])] - - if par.use_sc: - Mpsc = np.zeros((indrp.size, nsta)) - Mssc = np.zeros((indrs.size, nsta)) - for ns in range(indrp.size): - Mpsc[ns, int(1.e-6+caldata[indrp[ns], 2])] = 1. - for ns in range(indrs.size): - Mssc[ns, int(1.e-6+caldata[indrs[ns], 2])] = 1. - - Msc_cal.append(sp.block_diag((sp.csr_matrix(Mpsc), sp.csr_matrix(Mssc)))) - - else: - ncal = 0 - tcal = np.array([]) - - if np.isscalar(Vinit[0]): - Vp = Vinit[0] + np.zeros(nnodes) - s_p = np.ones(nnodes)/Vinit[0] - else: - Vp = Vinit[0] - s_p = 1./Vinit[0] - if np.isscalar(Vinit[1]): - Vs = Vinit[1] + np.zeros(nnodes) - s_s = np.ones(nnodes)/Vinit[1] - else: - Vs = Vinit[1] - s_s = 1./Vinit[1] - if par.invert_VsVp: - VsVp = Vs/Vp - V = np.hstack((Vp, VsVp)) - else: - V = np.hstack((Vp, Vs)) - - if par.verbose: - print('\n *** Joint hypocenter-velocity inversion -- P and S-wave data ***\n') - - if par.invert_vel: - resV = np.zeros(par.maxit+1) - resAxb = np.zeros(par.maxit) - - P = np.ones(2*nnodes) - dP = sp.csr_matrix((np.ones(2*nnodes), (np.arange(2*nnodes, dtype=np.int64), - np.arange(2*nnodes, dtype=np.int64))), - shape=(2*nnodes, 2*nnodes)) - - deltam = np.zeros(2*nnodes+2*nsta) - - if par.constr_sc and par.use_sc: - u1 = np.ones(2*nnodes+2*nsta) - u1[:2*nnodes] = 0.0 - u1[(2*nnodes+nsta):] = 0.0 - - i = np.arange(2*nnodes, 2*nnodes+nsta) - j = np.arange(2*nnodes, 2*nnodes+nsta) - nn = i.size - i = np.kron(i, np.ones((nn,))) - j = np.kron(np.ones((nn,)), j) - u1Tu1 = sp.csr_matrix((np.ones((i.size,)), (i, j)), - shape=(2*nnodes+2*nsta, 2*nnodes+2*nsta)) - else: - u1 = np.zeros(2*nnodes+2*nsta) - u1Tu1 = sp.csr_matrix((2*nnodes+2*nsta, 2*nnodes+2*nsta)) - - if Vpts.size > 0: - - if par.verbose: - print('Building velocity data point matrix D') - sys.stdout.flush() - - if par.invert_VsVp: - Vpts2 = Vpts.copy() - i_p = np.nonzero(Vpts[:, 4]==0.0)[0] - i_s = np.nonzero(Vpts[:, 4]==1.0)[0] - for i in i_s: - for ii in i_p: - d = np.sqrt(np.sum((Vpts2[i, 1:4]-Vpts[ii, 1:4])**2)) - if d < 0.00001: - Vpts2[i, 0] = Vpts[i, 0]/Vpts[ii, 0] - break - else: - raise ValueError('Missing Vp data point for Vs data at ({0:f}, {1:f}, {2:f})'.format(Vpts[i, 1], Vpts[i, 2], Vpts[i, 3])) - - else: - Vpts2 = Vpts - - if par.invert_VsVp: - D = grid.compute_D(Vpts2[:, 1:4]) - D = sp.hstack((D, sp.coo_matrix(D.shape))).tocsr() - else: - i_p = Vpts2[:, 4]==0.0 - i_s = Vpts2[:, 4]==1.0 - Dp = grid.compute_D(Vpts2[i_p, 1:4]) - Ds = grid.compute_D(Vpts2[i_s, 1:4]) - D = sp.block_diag((Dp, Ds)).tocsr() - - D1 = sp.hstack((D, sp.csr_matrix((Vpts2.shape[0], 2*nsta)))).tocsr() - else: - D = 0.0 - - if par.verbose: - print('Building regularization matrix K') - sys.stdout.flush() - Kx, Ky, Kz = grid.compute_K() - Kx = sp.block_diag((Kx, Kx)) - Ky = sp.block_diag((Ky, Ky)) - Kz = sp.block_diag((Kz, Kz)) - Kx1 = sp.hstack((Kx, sp.coo_matrix((2*nnodes, 2*nsta)))).tocsr() - KtK = Kx1.T * Kx1 - Ky1 = sp.hstack((Ky, sp.coo_matrix((2*nnodes, 2*nsta)))).tocsr() - KtK += Ky1.T * Ky1 - Kz1 = sp.hstack((Kz, sp.coo_matrix((2*nnodes, 2*nsta)))).tocsr() - KtK += par.wzK * Kz1.T * Kz1 - nK = spl.norm(KtK) - Kx = Kx.tocsr() - Ky = Ky.tocsr() - Kz = Kz.tocsr() - else: - resV = None - resAxb = None - - if par.invert_VsVp: - VsVpmin = par.Vsmin/par.Vpmax - VsVpmax = par.Vsmax/par.Vpmin - - if par.verbose: - print('\nStarting iterations') - - for it in np.arange(par.maxit): - - par._final_iteration = it == par.maxit-1 - - if par.invert_vel: - if par.verbose: - print('\nIteration {0:d} - Updating velocity model'.format(it+1)) - print(' Updating penalty vector') - sys.stdout.flush() - - # compute vector C - cx = Kx * V - cy = Ky * V - cz = Kz * V - - # compute dP/dV, matrix of penalties derivatives - for n in np.arange(nnodes): - if V[n] < par.Vpmin: - P[n] = par.PAp * (par.Vpmin-V[n]) - dP[n, n] = -par.PAp - elif V[n] > par.Vpmax: - P[n] = par.PAp * (V[n]-par.Vpmax) - dP[n, n] = par.PAp - else: - P[n] = 0.0 - dP[n, n] = 0.0 - for n in np.arange(nnodes, 2*nnodes): - if par.invert_VsVp: - if VsVp[n-nnodes] < VsVpmin: - P[n] = par.PAs * (VsVpmin-VsVp[n-nnodes]) - dP[n, n] = -par.PAs - elif VsVp[n-nnodes] > VsVpmax: - P[n] = par.PAs * (VsVp[n-nnodes]-VsVpmax) - dP[n, n] = par.PAs - else: - P[n] = 0.0 - dP[n, n] = 0.0 - else: - if Vs[n-nnodes] < par.Vsmin: - P[n] = par.PAs * (par.Vsmin-Vs[n-nnodes]) - dP[n, n] = -par.PAs - elif Vs[n-nnodes] > par.Vsmax: - P[n] = par.PAs * (Vs[n-nnodes]-par.Vsmax) - dP[n, n] = par.PAs - else: - P[n] = 0.0 - dP[n, n] = 0.0 - - if par.verbose: - npel = np.sum(P[:nnodes] != 0.0) - if npel > 0: - print(' P-wave penalties applied at {0:d} nodes'.format(npel)) - npel = np.sum(P[nnodes:2*nnodes] != 0.0) - if npel > 0: - print(' S-wave penalties applied at {0:d} nodes'.format(npel)) - - - if par.verbose: - print(' Raytracing') - sys.stdout.flush() - - if nev > 0: - hyp = np.empty((nttp, 5)) - for ne in np.arange(nev): - indh = np.nonzero(hyp0[:, 0] == evID[ne])[0] - indrp = np.nonzero(np.logical_and(data[:, 0] == evID[ne], indp))[0] - for i in indrp: - hyp[i, :] = hyp0[indh[0], :] - - tcalcp, raysp, Mevp = grid.raytrace(hyp, rcv_datap, s_p, - return_rays=True, - compute_M=True) - s0p = grid.get_s0(hyp) - - hyp = np.empty((ntts, 5)) - for ne in np.arange(nev): - indh = np.nonzero(hyp0[:, 0] == evID[ne])[0] - indrs = np.nonzero(np.logical_and(data[:, 0] == evID[ne], inds))[0] - for i in indrs: - hyp[i-nttp, :] = hyp0[indh[0], :] - tcalcs, rayss, Mevs = grid_s.raytrace(hyp, rcv_datas, s_s, - return_rays=True, - compute_M=True) - s0s = grid.get_s0(hyp) - - tcalc = np.hstack((tcalcp, tcalcs)) - s0 = np.hstack((s0p, s0s)) - rays = [] - for r in raysp: - rays.append(r) - for r in rayss: - rays.append(r) - - # Merge Mevp & Mevs - Mev = [None] * nev - for ne in np.arange(nev): - - Mp = Mevp[ne] - Ms = Mevs[ne] - - if par.invert_VsVp: - # Block 1991, p. 45 - tmp1 = Ms.multiply(np.tile(VsVp, (Ms.shape[0], 1))) - tmp2 = Ms.multiply(np.tile(Vp, (Ms.shape[0], 1))) - tmp2 = sp.hstack((tmp1, tmp2)) - tmp1 = sp.hstack((Mp, sp.csr_matrix(Mp.shape))) - Mev[ne] = sp.vstack((tmp1, tmp2)) - else: - Mev[ne] = sp.block_diag((Mp, Ms)) - - if par.use_sc: - indrp = np.nonzero(np.logical_and(data[:, 0] == evID[ne], indp))[0] - indrs = np.nonzero(np.logical_and(data[:, 0] == evID[ne], inds))[0] - - Mpsc = np.zeros((indrp.size, nsta)) - Mssc = np.zeros((indrs.size, nsta)) - for ns in range(indrp.size): - Mpsc[ns, int(1.e-6+data[indrp[ns], 2])] = 1. - for ns in range(indrs.size): - Mssc[ns, int(1.e-6+data[indrs[ns], 2])] = 1. - - Msc = sp.block_diag((sp.csr_matrix(Mpsc), sp.csr_matrix(Mssc))) - # add terms for station corrections after terms for velocity because - # solution vector contains [Vp Vs sc_p sc_s] in that order - Mev[ne] = sp.hstack((Mev[ne], Msc)) - - else: - tcalc = np.array([]) - - if ncal > 0: - tcalcp_cal, Mp_cal = grid.raytrace(hcalp, rcv_calp, s_p, - compute_M=True) - if nttcals > 0: - tcalcs_cal, Ms_cal = grid_s.raytrace(hcals, rcv_cals, s_s, - compute_M=True) - tcalc_cal = np.hstack((tcalcp_cal, tcalcs_cal)) - else: - tcalc_cal = tcalcp_cal - else: - tcalc_cal = np.array([]) - - r1a = tobs - tcalc - r1 = tcal - tcalc_cal - if r1a.size > 0: - r1 = np.hstack((np.zeros(data.shape[0]-4*nev), r1)) - - resV[it] = np.linalg.norm(np.hstack((tobs-tcalc, tcal-tcalc_cal))) - - if par.show_plots: - plt.figure(1) - plt.cla() - plt.plot(r1a, 'o') - plt.title('Residuals - Iteration {0:d}'.format(it+1)) - plt.show(block=False) - plt.pause(0.0001) - - # initializing matrix M; matrix of partial derivatives of velocity dt/dV - if par.verbose: - print(' Building matrix M') - sys.stdout.flush() - - M1 = None - ir1 = 0 - for ne in range(nev): + if np.sum(np.abs(deltah) < par.conv_hypo) == 2: if par.verbose: - print(' Event ID '+str(int(1.e-6+evID[ne]))) + print(' - converged at iteration '+str(itt+1)) sys.stdout.flush() + break - indh = np.nonzero(hyp0[:, 0] == evID[ne])[0] - indr = np.nonzero(data[:, 0] == evID[ne])[0] - - nst = np.sum(indr.size) - nst2 = nst-4 - H = np.ones((nst, 4)) - for ns in range(nst): - raysi = rays[indr[ns]] - S0 = s0[indr[ns]] - - d = (raysi[1, :]-hyp0[indh, 2:]).flatten() - ds = np.sqrt(np.sum(d*d)) - H[ns, 1] = -S0 * d[0]/ds - H[ns, 2] = -S0 * d[1]/ds - H[ns, 3] = -S0 * d[2]/ds - - Q, _ = np.linalg.qr(H, mode='complete') - T = sp.csr_matrix(Q[:, 4:]).T - M = T * Mev[ne] - - if M1 is None: - M1 = M - else: - M1 = sp.vstack((M1, M)) - - r1[ir1 + np.arange(nst2, dtype=np.int64)] = T.dot(r1a[indr]) - ir1 += nst2 - - for nc in range(ncal): - Mp = Mp_cal[nc] - if nttcals > 0: - Ms = Ms_cal[nc] - else: - Ms = sp.csr_matrix([]) - - if par.invert_VsVp: - if nttcals > 0: - # Block 1991, p. 45 - tmp1 = Ms.multiply(np.tile(VsVp, (Ms.shape[0], 1))) - tmp2 = Ms.multiply(np.tile(Vp, (Ms.shape[0], 1))) - tmp2 = sp.hstack((tmp1, tmp2)) - tmp1 = sp.hstack((Mp, sp.csr_matrix(Mp.shape))) - M = sp.vstack((tmp1, tmp2)) - else: - M = sp.hstack((Mp, sp.csr_matrix(Mp.shape))) - else: - M = sp.block_diag((Mp, Ms)) - - if par.use_sc: - M = sp.hstack((M, Msc_cal[nc])) - - if M1 is None: - M1 = M - else: - M1 = sp.vstack((M1, M)) - - if par.verbose: - print(' Assembling matrices and solving system') - sys.stdout.flush() - - s = -np.sum(sc_p) - - dP1 = sp.hstack((dP, sp.csr_matrix((2*nnodes, 2*nsta)))).tocsr() # dP prime - - # compute A & h for inversion - - M1 = M1.tocsr() - - A = M1.T * M1 - - nM = spl.norm(A) - λ = par.λ * nM / nK - - A += λ*KtK - - tmp = dP1.T * dP1 - nP = spl.norm(tmp) - if nP != 0.0: - γ = par.γ * nM / nP - else: - γ = par.γ - - A += γ*tmp - A += u1Tu1 - - b = M1.T * r1 - tmp2x = Kx1.T * cx - tmp2y = Ky1.T * cy - tmp2z = Kz1.T * cz - tmp3 = dP1.T * P - tmp = u1 * s - b += -λ*tmp2x - λ*tmp2y - par.wzK*λ*tmp2z - γ*tmp3 - tmp - - if Vpts2.shape[0] > 0: - tmp = D1.T * D1 - nD = spl.norm(tmp) - α = par.α * nM / nD - A += α * tmp - b += α * D1.T * (Vpts2[:, 0] - D*V) - - if par.verbose: - print(' calling minres with system of size {0:d} x {1:d}'.format(A.shape[0], A.shape[1])) - sys.stdout.flush() - x = spl.minres(A, b) - - deltam = x[0] - resAxb[it] = np.linalg.norm(A*deltam - b) - - dmean = np.mean(np.abs(deltam[:nnodes])) - if dmean > par.dVp_max: - if par.verbose: - print(' Scaling Vp perturbations by {0:e}'.format(par.dVp_max/dmean)) - deltam[:nnodes] = deltam[:nnodes] * par.dVp_max/dmean - dmean = np.mean(np.abs(deltam[nnodes:2*nnodes])) - if dmean > par.dVs_max: - if par.verbose: - print(' Scaling Vs perturbations by {0:e}'.format(par.dVs_max/dmean)) - deltam[nnodes:2*nnodes] = deltam[nnodes:2*nnodes] * par.dVs_max/dmean - - V += deltam[:2*nnodes] - Vp = V[:nnodes] - if par.invert_VsVp: - VsVp = V[nnodes:2*nnodes] - Vs = VsVp * Vp - else: - Vs = V[nnodes:2*nnodes] - s_p = 1./Vp - s_s = 1./Vs - sc_p += deltam[2*nnodes:2*nnodes+nsta] - sc_s += deltam[2*nnodes+nsta:] - - if par.save_V: - if par.verbose: - print(' Saving Velocity models') - if 'vtk' in sys.modules: - grid.to_vtk({'Vp': Vp}, 'Vp{0:02d}'.format(it+1)) - if par.invert_VsVp: - if 'vtk' in sys.modules: - grid.to_vtk({'VsVp': VsVp}, 'VsVp{0:02d}'.format(it+1)) - if 'vtk' in sys.modules: - grid.to_vtk({'Vs': Vs}, 'Vs{0:02d}'.format(it+1)) - - if nev > 0: + else: if par.verbose: - print('Iteration {0:d} - Relocating events'.format(it+1)) + print(' - reached max number of iterations') sys.stdout.flush() - if grid.nthreads == 1 or nev < grid.nthreads: - for ne in range(nev): - _relocPS(ne, par, (grid, grid_s), evID, hyp0, data, rcv, tobs, (s_p, s_s), (indp, inds)) - else: - # run in parallel - blk_size = np.zeros((grid.nthreads,), dtype=np.int64) - nj = nev - while nj > 0: - for n in range(grid.nthreads): - blk_size[n] += 1 - nj -= 1 - if nj == 0: - break - processes = [] - blk_start = 0 - h_queue = Queue() - for n in range(grid.nthreads): - blk_end = blk_start + blk_size[n] - p = Process(target=_rlPS_worker, - args=(n, blk_start, blk_end, par, (grid, grid_s), evID, hyp0, - data, rcv, tobs, (s_p, s_s), (indp, inds), h_queue), - daemon=True) - processes.append(p) - p.start() - blk_start += blk_size[n] - - for ne in range(nev): - h, indh = h_queue.get() - hyp0[indh, :] = h + if par.verbose: + print(' Updating all hypocenter params', end='') + sys.stdout.flush() - if par.invert_vel: - if nev > 0: - hyp = np.empty((nttp, 5)) - for ne in np.arange(nev): - indh = np.nonzero(hyp0[:, 0] == evID[ne])[0] - indrp = np.nonzero(np.logical_and(data[:, 0] == evID[ne], indp))[0] - for i in indrp: - hyp[i, :] = hyp0[indh[0], :] + H = np.ones((nst, 4)) + for itt in range(par.maxit_hypo): + for i in range(nst): + hyp[i, :] = hyp0[indh, :] - tcalcp = grid.raytrace(hyp, rcv_datap, s_p) + tcalc, rays = grid.raytrace(hyp, stn, thread_no=thread_no, + return_rays=True) + s0 = grid.get_s0(hyp) + for ns in range(nst): + raysi = rays[ns] + S0 = s0[ns] - hyp = np.empty((ntts, 5)) - for ne in np.arange(nev): - indh = np.nonzero(hyp0[:, 0] == evID[ne])[0] - indrs = np.nonzero(np.logical_and(data[:, 0] == evID[ne], inds))[0] - for i in indrs: - hyp[i-nttp, :] = hyp0[indh[0], :] - tcalcs = grid_s.raytrace(hyp, rcv_datas, s_s) + d = (raysi[1, :]-hyp0[indh, 2:]).flatten() + ds = np.sqrt(np.sum(d*d)) + H[ns, 1] = -S0 * d[0]/ds + H[ns, 2] = -S0 * d[1]/ds + H[ns, 3] = -S0 * d[2]/ds - tcalc = np.hstack((tcalcp, tcalcs)) - else: - tcalc = np.array([]) + r = tobs[indr] - tcalc + x = lstsq(H, r) + deltah = x[0] - if ncal > 0: - tcalcp_cal = grid.raytrace(hcalp, rcv_calp, s_p) - tcalcs_cal = grid_s.raytrace(hcals, rcv_cals, s_s) - tcalc_cal = np.hstack((tcalcp_cal, tcalcs_cal)) - else: - tcalc_cal = np.array([]) + if np.sum(np.isfinite(deltah)) != deltah.size: + try: + U,S,VVh = np.linalg.svd(H.T.dot(H)+1e-9*np.eye(4)) + VV = VVh.T + deltah = np.dot(VV, np.dot(U.T, H.T.dot(r))/S) + except np.linalg.linalg.LinAlgError: + print(' Event could not be relocated, resetting and exiting') + hyp0[indh, :] = hyp_save + return hyp_save, indh - r1a = tobs - tcalc - r1 = tcal - tcalc_cal - if r1a.size > 0: - r1 = np.hstack((np.zeros(data.shape[0]-4*nev), r1)) + if np.abs(deltah[0]) > par.dt_max: + deltah[0] = par.dt_max * np.sign(deltah[0]) + for n in range(1, 4): + if np.abs(deltah[n]) > par.dx_max: + deltah[n] = par.dx_max * np.sign(deltah[n]) - if par.show_plots: - plt.figure(1) - plt.cla() - plt.plot(r1a, 'o') - plt.title('Residuals - Final step') - plt.show(block=False) - plt.pause(0.0001) + new_hyp = hyp0[indh, 1:] + deltah + if grid.is_outside(new_hyp[1:].reshape((1, 3))): + print(' Event could not be relocated inside the grid ({0:f}, {1:f}, {2:f}), resetting and exiting'.format(new_hyp[1], new_hyp[2], new_hyp[3])) + hyp0[indh, :] = hyp_save + return hyp_save, indh - r1 = r1.reshape(-1, 1) - r1a = r1a.reshape(-1, 1) + hyp0[indh, 1:] += deltah - resV[-1] = np.linalg.norm(np.hstack((tobs-tcalc, tcal-tcalc_cal))) + if np.sum(np.abs(deltah[1:]) < par.conv_hypo) == 3: + if par.verbose: + print(' - converged at iteration '+str(itt+1)) + sys.stdout.flush() + break - if par.verbose: - print('\n ** Inversion complete **\n', flush=True) + else: + if par.verbose: + print(' - reached max number of iterations') + sys.stdout.flush() - return hyp0, (Vp, Vs), (sc_p, sc_s), (resV, resAxb) + if par.save_rp and 'vtk' in sys.modules and par._final_iteration: + if par.verbose: + print(' Saving raypaths') + filename = 'raypaths' + key = 'ev_{0:d}'.format(int(1.e-6+evID[ne])) + grid.to_vtk({key: rays}, filename) + + return hyp0[indh, :], indh -def jointHypoVelPS_c(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), - Vpts=np.array([])): +def jointHypoVelPS(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), + Vpts=np.array([])): """ Joint hypocenter-velocity inversion on a regular grid @@ -2075,6 +974,7 @@ def jointHypoVelPS_c(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), hyp0 = hinit.copy() # set origin time to 0 and offset data accordingly + data = data.copy() hyp0_dt = {} for n in range(hyp0.shape[0]): hyp0_dt[hyp0[n, 0]] = hyp0[n, 1] @@ -2605,13 +1505,13 @@ def jointHypoVelPS_c(par, grid, data, rcv, Vinit, hinit, caldata=np.array([]), dmax = np.max(np.abs(deltam[:nslowness])) if dmax > dVp_max: if par.verbose: - print(' Scaling P slowness perturbations by {0:e}'.format(1./(dVp_max*dmax))) - deltam[:nslowness] = deltam[:nslowness] / (dVp_max*dmax) + print(' Scaling P slowness perturbations by {0:e}'.format(dVp_max/dmax)) + deltam[:nslowness] = deltam[:nslowness] * dVp_max / dmax dmax = np.max(np.abs(deltam[nslowness:2*nslowness])) if dmax > dVs_max: if par.verbose: - print(' Scaling S slowness perturbations by {0:e}'.format(1./(dVs_max*dmax))) - deltam[nslowness:2*nslowness] = deltam[nslowness:2*nslowness] / (dVs_max*dmax) + print(' Scaling S slowness perturbations by {0:e}'.format(dVs_max/dmax)) + deltam[nslowness:2*nslowness] = deltam[nslowness:2*nslowness] * dVs_max / dmax s += deltam[:2*nslowness] s_p = s[:nslowness] @@ -2774,13 +1674,33 @@ def _relocPS(ne, par, grid, evID, hyp0, data, rcv, tobs, s, ind, thread_no=None) for itt in range(par.maxit_hypo): for i in range(nstp): hypp[i, :] = hyp0[indh, :] - tcalcp, raysp = grid_p.raytrace(hypp, stnp, s_p, thread_no, - return_rays=True) + + try: + tcalcp, raysp = grid_p.raytrace(hypp, stnp, s_p, thread_no, + return_rays=True) + except RuntimeError as rte: + if 'going outside grid' in str(rte): + print(' Problem while computing P-wave traveltimes, resetting and exiting') + hyp0[indh, :] = hyp_save + return hyp_save, indh + else: + raise rte + s0p = grid_p.get_s0(hypp) for i in range(nsts): hyps[i, :] = hyp0[indh, :] - tcalcs, rayss = grid_s.raytrace(hyps, stns, s_s, thread_no, - return_rays=True) + + try: + tcalcs, rayss = grid_s.raytrace(hyps, stns, s_s, thread_no, + return_rays=True) + except RuntimeError as rte: + if 'going outside grid' in str(rte): + print(' Problem while computing S-wave traveltimes, resetting and exiting') + hyp0[indh, :] = hyp_save + return hyp_save, indh + else: + raise rte + s0s = grid_s.get_s0(hyps) for ns in range(nstp): raysi = raysp[ns] @@ -2954,13 +1874,11 @@ def _relocPS(ne, par, grid, evID, hyp0, data, rcv, tobs, s, ind, thread_no=None) g2 = Grid3d(x, y, z, nthreads, cell_slowness=True) - testK = False - testP = False - testPS = False - testParallel = False + testK = True + testParallel = True addNoise = True testC = True - testCp = False + testCp = True testCps = True if testK: @@ -3091,7 +2009,7 @@ def _relocPS(ne, par, grid, evID, hyp0, data, rcv, tobs, s, ind, thread_no=None) plt.show() - if testP or testPS or testParallel: + if testParallel: g = g1 @@ -3206,8 +2124,6 @@ def Vz2(z): par = InvParams(maxit=3, maxit_hypo=10, conv_hypo=0.001, Vlim=Vlim, dmax=dmax, lagrangians=lagran, invert_vel=True, verbose=True) - if testParallel: - g = g1 tt = g.raytrace(src, rcv_data, slowness) @@ -3217,206 +2133,6 @@ def Vz2(z): print('done') - if testP: - - g = g1 - - tt += noise_variance*np.random.randn(tt.size) - - data = np.hstack((src[:, 0].reshape((-1, 1)), tt.reshape((-1, 1)), ircv_data)) - - hinit2, res = hypoloc(data, rcv, Vinit, hinit, 10, 0.001, True) - -# par.invert_vel = False -# par.maxit = 1 -# h, V, sc, res = jointHypoVel(par, g, data, rcv, Vp.flatten(), hinit2) - - h, V, sc, res = jointHypoVel(par, g, data, rcv, Vpinit, hinit2, caldata=caldata, Vpts=Vpts) - - plt.figure() - plt.plot(res[0]) - plt.show(block=False) - - err_xc = hinit2[:, 2:5] - h_true[:, 2:5] - err_xc = np.sqrt(np.sum(err_xc**2, axis=1)) - err_tc = hinit2[:, 1] - h_true[:, 1] - - err_x = h[:, 2:5] - h_true[:, 2:5] - err_x = np.sqrt(np.sum(err_x**2, axis=1)) - err_t = h[:, 1] - h_true[:, 1] - - plt.figure(figsize=(10, 4)) - plt.subplot(121) - plt.plot(err_x, 'o', label=r'$\|\|\Delta x\|\|$ = {0:6.5f}'.format(np.linalg.norm(err_x))) - plt.plot(err_xc, 'r*', label=r'$\|\|\Delta x\|\|$ = {0:6.5f}'.format(np.linalg.norm(err_xc))) - plt.ylabel(r'$\Delta x$') - plt.xlabel('Event ID') - plt.legend() - plt.subplot(122) - plt.plot(np.abs(err_t), 'o', label=r'$\|\|\Delta t\|\|$ = {0:6.5f}'.format(np.linalg.norm(err_t))) - plt.plot(np.abs(err_tc), 'r*', label=r'$\|\|\Delta t\|\|$ = {0:6.5f}'.format(np.linalg.norm(err_tc))) - plt.ylabel(r'$\Delta t$') - plt.xlabel('Event ID') - plt.legend() - - plt.show(block=False) - - V3d = V.reshape(g.shape) - - plt.figure(figsize=(10, 8)) - plt.subplot(221) - plt.pcolor(x, z, np.squeeze(V3d[:, 9, :].T), cmap='CMRmap',), plt.gca().invert_yaxis() - plt.xlabel('X') - plt.ylabel('Z') - plt.colorbar() - plt.subplot(222) - plt.pcolor(y, z, np.squeeze(V3d[8, :, :].T), cmap='CMRmap'), plt.gca().invert_yaxis() - plt.xlabel('Y') - plt.ylabel('Z') - plt.colorbar() - plt.subplot(223) - plt.pcolor(x, y, np.squeeze(V3d[:, :, 4].T), cmap='CMRmap') - plt.xlabel('X') - plt.ylabel('Y') - plt.colorbar() - - plt.show(block=False) - - plt.figure() - plt.plot(sc, 'o') - plt.xlabel('Station no') - plt.ylabel('Correction') - plt.show() - -# %% testPS - if testPS: - - g = g1 - - Vpts_s = np.array([[Vs, 0.100, 0.100, 0.001, 1.0], - [Vs, 0.100, 0.200, 0.001, 1.0], - [Vs, 0.200, 0.100, 0.001, 1.0], - [Vs, 0.200, 0.200, 0.001, 1.0], - [Vs, 0.112, 0.148, 0.011, 1.0], - [Vs, 0.152, 0.108, 0.005, 1.0], - [Vs, 0.152, 0.108, 0.075, 1.0], - [Vs, 0.192, 0.148, 0.011, 1.0]]) - - Vpts = np.vstack((Vpts, Vpts_s)) - - slowness_s = 1./Vs + np.zeros(g.get_number_of_nodes()) - - tt_s = g.raytrace(src, rcv_data, slowness_s) - - tt += noise_variance*np.random.randn(tt.size) - tt_s += noise_variance*np.random.randn(tt_s.size) - - # remove some values - ind_p = np.ones(tt.shape[0], dtype=bool) - ind_p[np.random.randint(ind_p.size, size=25)] = False - ind_s = np.ones(tt_s.shape[0], dtype=bool) - ind_s[np.random.randint(ind_s.size, size=25)] = False - - data_p = np.hstack((src[ind_p, 0].reshape((-1, 1)), tt[ind_p].reshape((-1, 1)), ircv_data[ind_p, :], np.zeros((np.sum(ind_p), 1)))) - data_s = np.hstack((src[ind_s, 0].reshape((-1, 1)), tt_s[ind_s].reshape((-1, 1)), ircv_data[ind_s, :], np.ones((np.sum(ind_s), 1)))) - - data = np.vstack((data_p, data_s)) - - tcal_s = g.raytrace(src_cal, rcv_cal, slowness_s) - caldata_s = np.column_stack((src_cal[:, 0], tcal_s, ircv_cal, src_cal[:, 2:], - np.ones(tcal_s.shape))) - caldata = np.vstack((caldata, caldata_s)) - - Vinit = (Vinit, 2.0) - - hinit2, res = hypolocPS(data, rcv, Vinit, hinit, 10, 0.001, True) - - par.save_V = True - par.save_rp = True - par.dVs_max = 0.01 - par.invert_VsVp = False - par.constr_sc = False - h, V, sc, res = jointHypoVelPS(par, g, data, rcv, (Vpinit, 2.0), hinit2, caldata=caldata, Vpts=Vpts) - - plt.figure() - plt.plot(res[0]) - plt.show(block=False) - - err_xc = hinit2[:, 2:5] - h_true[:, 2:5] - err_xc = np.sqrt(np.sum(err_xc**2, axis=1)) - err_tc = hinit2[:, 1] - h_true[:, 1] - - err_x = h[:, 2:5] - h_true[:, 2:5] - err_x = np.sqrt(np.sum(err_x**2, axis=1)) - err_t = h[:, 1] - h_true[:, 1] - - plt.figure(figsize=(10, 4)) - plt.subplot(121) - plt.plot(err_x, 'o', label=r'$\|\|\Delta x\|\|$ = {0:6.5f}'.format(np.linalg.norm(err_x))) - plt.plot(err_xc, 'r*', label=r'$\|\|\Delta x\|\|$ = {0:6.5f}'.format(np.linalg.norm(err_xc))) - plt.ylabel(r'$\Delta x$') - plt.xlabel('Event ID') - plt.legend() - plt.subplot(122) - plt.plot(np.abs(err_t), 'o', label=r'$\|\|\Delta t\|\|$ = {0:6.5f}'.format(np.linalg.norm(err_t))) - plt.plot(np.abs(err_tc), 'r*', label=r'$\|\|\Delta t\|\|$ = {0:6.5f}'.format(np.linalg.norm(err_tc))) - plt.ylabel(r'$\Delta t$') - plt.xlabel('Event ID') - plt.legend() - - plt.show(block=False) - - V3d = V[0].reshape(g.shape) - - plt.figure(figsize=(10, 8)) - plt.subplot(221) - plt.pcolor(x, z, np.squeeze(V3d[:, 9, :].T), cmap='CMRmap',), plt.gca().invert_yaxis() - plt.xlabel('X') - plt.ylabel('Z') - plt.colorbar() - plt.subplot(222) - plt.pcolor(y, z, np.squeeze(V3d[8, :, :].T), cmap='CMRmap'), plt.gca().invert_yaxis() - plt.xlabel('Y') - plt.ylabel('Z') - plt.colorbar() - plt.subplot(223) - plt.pcolor(x, y, np.squeeze(V3d[:, :, 4].T), cmap='CMRmap') - plt.xlabel('X') - plt.ylabel('Y') - plt.colorbar() - plt.suptitle('V_p') - - plt.show(block=False) - - V3d = V[1].reshape(g.shape) - - plt.figure(figsize=(10, 8)) - plt.subplot(221) - plt.pcolor(x, z, np.squeeze(V3d[:, 9, :].T), cmap='CMRmap',), plt.gca().invert_yaxis() - plt.xlabel('X') - plt.ylabel('Z') - plt.colorbar() - plt.subplot(222) - plt.pcolor(y, z, np.squeeze(V3d[8, :, :].T), cmap='CMRmap'), plt.gca().invert_yaxis() - plt.xlabel('Y') - plt.ylabel('Z') - plt.colorbar() - plt.subplot(223) - plt.pcolor(x, y, np.squeeze(V3d[:, :, 4].T), cmap='CMRmap') - plt.xlabel('X') - plt.ylabel('Y') - plt.colorbar() - plt.suptitle('V_s') - - plt.show(block=False) - - plt.figure() - plt.plot(sc[0], 'o', label='P-wave') - plt.plot(sc[1], 'r*', label='s-wave') - plt.xlabel('Station no') - plt.ylabel('Correction') - plt.legend() - plt.show() # %% testC if testC: @@ -3536,7 +2252,7 @@ def Vz2(z): par = InvParams(maxit=3, maxit_hypo=10, conv_hypo=0.001, Vlim=Vlim, dmax=dmax, lagrangians=lagran, invert_vel=True, - verbose=True, save_V=True) + verbose=True, save_V=True, show_plots=True) # %% testCp if testCp: @@ -3547,7 +2263,7 @@ def Vz2(z): hinit2, res = hypoloc(data, rcv, Vinit, hinit, 10, 0.001, True) - h, V, sc, res = jointHypoVel_c(par, g, data, rcv, Vpinit, hinit2, caldata=caldata, Vpts=Vpts) + h, V, sc, res = jointHypoVel(par, g, data, rcv, Vpinit, hinit2, caldata=caldata, Vpts=Vpts) plt.figure() plt.plot(res[0]) @@ -3652,10 +2368,11 @@ def Vz2(z): par.save_rp = True par.dVs_max = 0.01 par.invert_VsVp = False - par.constr_sc = False - h, V, sc, res = jointHypoVelPS_c(par, g, data, rcv, (Vpinit, 2.0), - hinit2, caldata=caldata, - Vpts=Vpts) + par.constr_sc = True + par.show_plots = True + h, V, sc, res = jointHypoVelPS(par, g, data, rcv, (Vpinit, 2.0), + hinit2, caldata=caldata, + Vpts=Vpts) plt.figure() plt.plot(res[0])