diff --git a/user/psava/Msmopick.py b/user/psava/Msmopick.py index 51381347a..1376d6135 100755 --- a/user/psava/Msmopick.py +++ b/user/psava/Msmopick.py @@ -50,8 +50,6 @@ def lapop2D(nx,nz): n1 = Fin.int ("n1") o1 = Fin.float("o1") d1 = Fin.float("d1") -l1 = 0.5 * n1*d1 # xcorr time window -z1 = int( abs(o1)/d1) nd = Fin.size(1); # number of traces nm = nd @@ -62,67 +60,52 @@ def lapop2D(nx,nz): Fou.put("o1",0.0) Fou.put('d1',1.0) -Fou.put("n2",3) # temporary - +Fou.put("n2",3) Fou.put("n3",1) Fou.put("n4",1) # ------------------------------------------------------------ # pick max value and weight # ------------------------------------------------------------ -din = np.zeros(n1,'f') +din = np.zeros(n1,'f') # xcorrelation pck = np.zeros(nd,'f') # picks wgh = np.zeros(nd,'f') # weights +# ------------------------------------------------------------ +l1 = 0.5 * n1*d1 # xcorr time window + for i in range(nd): Fin.read(din) - # raw picks - #pck[i] = o1 + np.argmax(din) * d1 - #wgh[i] = np.max(din) - - din += 1.0 # avoid zero weights - pck[i] = o1 + z1 * d1 - wgh[i] = din[z1] - for j in range(n1): - if wgh[i] < din[j]: - pck[i] = o1 + j * d1 - wgh[i] = din[j] + pck[i] = o1 + np.argmax(din) * d1 # raw picks + wgh[i] = np.max(din) # data weight -# bias weights +# pick bias pckmed = np.median(pck) -for i in range(nd): - #wgh[i]*= np.exp(- abs(pck[i]) ) - wgh[i] *= np.cos( 0.5*np.pi * abs((pck[i]-pckmed)/l1) ) +# data weight +for i in range(nd): + wgh[i] *= np.cos( 0.25*np.pi * abs((pck[i]-pckmed)/l1) ) wgh /= np.max(wgh) +wgh = 1 * np.power(wgh-np.min(wgh),wpo) # ------------------------------------------------------------ -# smooth picks +# setup inverse problem # ------------------------------------------------------------ - -xx = np.arange(nd) -xS = xx[0] -xE = xx[nd-1] - -nb = np.max([1,nd//10]) -#mbar = pck*0 # reference model -mbar = pck*0 + pckmed +mbar = pck*0 + pckmed # reference model dbar = pck # rough picks -wgh = np.power(wgh-np.min(wgh),wpo) - WDop = spdiags(wgh,[0],nd,nd) # data weight operator Gop = idnop1D(nd) # mapping operator -Rop = lapop1D(nm) # Laplacian operator +Rop = lapop1D(nm) # regularization operator # ------------------------------------------------------------ # L-curve -ne = par.int('ne',1) +# ------------------------------------------------------------ +ne = par.int ('ne',1) oe = par.float('oe',0.0) de = par.float('de',+0.1) ee = np.arange(oe, oe+ne*de, de) -#print(ee,file=sys.stderr) se = np.power(10,ee) ME = np.zeros( (nm,ne), dtype='float') @@ -130,35 +113,36 @@ def lapop2D(nx,nz): rm = np.zeros( ne, dtype='float') for ie in range(ne): - WMop = Rop / se[ie] # setup regularization - # solve IP + # scale the regularization operator + WMop = Rop / se[ie] + + # solve the IP modE = spsolve( (WDop*Gop).T * (WDop*Gop) + WMop.T * WMop , \ (WDop*Gop).T * WDop * dbar + WMop.T * WMop * mbar) - ME[:,ie] = modE # store model + + # store the model + ME[:,ie] = modE # compute residual norms rd[ie] = np.linalg.norm( WDop * (Gop * modE - dbar)) rm[ie] = np.linalg.norm( Rop * ( modE - mbar)) -if rd > 0: - rdn = rd / np.max(rd) # normalized residual norm -else: - rdn = 1.0 -if rm > 0: - rmn = rm / np.max(rm) -else: - rmn = 1.0 -rrn = rdn**2 + rmn**2 # distance from origin +rdn = rd / np.max(rd) # normalized the data residual norm +rmn = rm / np.max(rm) # normalized the model residual norm +rrn = rdn**2 + rmn**2 # compute the distance from origin + je = np.argmin(rrn) # index of the optimal model -spk = ME[:,je] # smooth picks +spk = ME[:,je] # optimal model + +# print(je,ee[je], file=sys.stderr) # ------------------------------------------------------------ -# write picks +# write picks and weights # ------------------------------------------------------------ Fou.write(spk) # smooth picks Fou.write(pck) # rough picks -Fou.write(wgh) +Fou.write(wgh) # xcor pick weight Fin.close() Fou.close()