Les 3 less 3 py from numpy

Les 3

less 3. py from numpy import * import rsf. api as rsf import matplotlib. pyplot as plt import os from pylab import * outzr=rsf. Output('zr. rsf') выходной rsf файл для приемников outzs=rsf. Output('zs. rsf') выходной rsf файл для источников outdens=rsf. Output('dens. rsf') плотностная модель outvel=rsf. Output('velmod 2. rsf') скоростная модель outricker=rsf. Output('ricker 1. rsf') форма импульса outdata=rsf. Output('data. rsf') сейсмограмма dx=0. 01 шаг по скоростной модели dz=0. 01 nt=2001 кол-во отсчетов по времени ot=0 начальное время dt=0. 001 шаг по времени kt=50 положение пика импульса рикера – номер отсчета

########## velocity model X = arange(-1. 5, dx); nx=len(X); ox=X[0] Z = arange(0, 2, dz); nz=len(Z); oz=Z[0] v=zeros((nz, nx), 'f') for i in range(nz): for j in range(nx): v[i][j]=2+exp((-((X[j]-1)**2+(Z[i]-1. 25)**2)/(0. 1)))-exp((-((X[j]+0. 4)**2+(Z[i]-1. 1)**2)/(0. 15))) outvel. put('n 1', nz) outvel. put('n 2', nx) outvel. put('o 1', oz) outvel. put('o 2', ox) outvel. put('d 1', dz) outvel. put('d 2', dx) outvel. write(v[: ]) outvel. close() plt. imshow(v, extent=[-1. 5, 2, 0]) hold(True) plt. savefig("aa. png")

aa. png

######### receivers r 1=-0. 8; r 2=0. 8; dr=0. 01 параметры приемной линии r = arange(r 1, r 2, dr); nr=len(r) zr=zeros((nr, 2), 'f') for i in range(nr): zr[i][0]=r[i] zr[i][1]=0. 2 zrt=zr. T транспонируем для отрисовки plt. plot(zrt[0, : ], zrt[1, : ]) рисовать графику можно только по первой оси outzr. put('n 1', 2) outzr. put('n 2', nr) outzr. put('d 1', 1) outzr. put('d 2', dr) outzr. put('o 1', r 1) outzr. put('o 2', r 1) outzr. write(zr[: ]) outzr. close()

######### interfaces kz=0. 2; hz=1. 4 s 1=-1. 2; s 2=0. 2; ds=0. 1 параметры положения источников k 1=-0. 2; k 2=1. 2; dk=0. 1 s = arange(s 1, s 2, ds); ns=len(s) k = arange(k 1, k 2, ds); nk=len(k) zs=zeros((ns, 2), 'f') zk=zeros((nk, 2), 'f') for i in range(ns): zs[i][0]=s[i] zs[i][1]=kz*s[i]+hz zst=zs. T plt. plot(zst[0, : ], zst[1, : ]) for i in range(nk): zk[i][0]=k[i] zk[i][1]=0. 1*sin(k[i]*10)+1. 8 zkt=zk. T plt. plot(zkt[0, : ], zkt[1, : ])

sss=concatenate((zs, zk), axis=1) outzs. put('n 1', 2) outzs. put('n 2', ns+nk) outzs. put('d 1', 1) outzs. put('d 2', ds) outzs. put('o 1', s 1) outzs. put('o 2', s 1) outzs. write(sss[: ]) outzs. close() hold(False) plt. savefig('model. png') command='display model. png' os. system(command)

model. png Синим – приемники, красным и зеленым - источники

########## density model density=ones((nz, nx), 'f') outdens. put('n 1', nz) outdens. put('n 2', nx) outdens. put('d 1', dz) outdens. put('d 2', dx) outdens. put('o 1', oz) outdens. put('o 2', ox) outdens. write(density[: ]) outdens. close()

########## wavelet sp=rsf. sfspike(nsp='1', mag='1', n 1=nt, d 1=dt, o 1=ot, k 1=kt, output='x 1')[0 ] rick=rsf. ricker 1(frequency='20')[sp] transp=rsf. sftransp(plane='12')[rick] d=rsf. sfwiggle(title=‘Ricker')[rick] d. show() outricker. put('n 1', nt) outricker. put('d 1', dt) outricker. put('o 1', 0) outricker. write(transp[: ]) outricker. close()

Форма импульса

########### awefd 1 data=rsf. awefd 2 d(verb='y', vel='velmod 2. rsf', sou='zs. rsf', rec='zr. rsf', den ='dens. rsf', free='n', dabc='y', nb='100', wfl='wfld. rsf', snap='y', expl='y', jsnap='50')[transp] решаем волновое уравнение (прямая задача) outdata. put('n 2', nt) outdata. put('n 1', nr) outdata. put('d 2', dt) outdata. put('d 1', dr) outdata. put('o 2', ot) outdata. put('o 1', r 1) outdata. write(data[: ]) outdata. close() dataseism=rsf. sfgrey(title='Seisogramm')['data. rsf'] dataseism. show()

sfawefd 2 d vel data dens rec sou Impulse form sfawefd 2 d fronts

data. rsf

input = rsf. Input("wfld. rsf") n 1=input. int("n 1"); n 2=input. int("n 2"); n 3=input. int("n 3") snnap = zeros((n 3, n 2, n 1), 'f') input. read(snnap) sna 1=snnap[4, : ] plt. imshow(sna 1, extent=[0, 2, -1. 5, 1. 5]) plt. savefig('snap. png') command 1='display snap. png'; os. system(command 1) рисуем первый фронт

4 -й фронт (200 отсчетов по времени)

######### exploding-reflector reverse-time migration out=rsf. Output('d. rsf') dataw=rsf. sfwindow(f 2=kt, n 2=nt-kt)[data] отрезаем первые 50 элементов datapad=rsf. sfpad(end 2=kt)[dataw] дописываем в конец 50 нулевых элементов tdat 1=rsf. sfreverse(which='2', opt='i', verb='y')[datapad] разворачиваем во времени out. put('n 1', nr) out. put('n 2', nt) out. put('o 1', r 1) out. put('o 2', ot) out. put('d 1', dr) out. put('d 2', dt) out. write(tdat 1[: ]) out. close() tdatgr=rsf. sfgrey(bias='0')[tdat 1] tdatgr. show() junk=rsf. awefd 2 d(verb='y', sou='zr. rsf', rec='zr. rsf', vel='velmod 2. rsf', den='dens. rsf', free=' n', dabc='y', nb='100', jsnap=(nt-1), wfl='mwfld. rsf', snap='y')[tdat 1] миграция!!!

Развернутая сейсмограмма

sfawefd 2 d vel junk dens sfawefd 2 d Sou=rec Impulse form=seism Fronts=migr result

######## migr graph inputm = rsf. Input("mwfld. rsf") n 1=inputm. int("n 1"); n 2=inputm. int("n 2"); n 3=inputm. int("n 3") snnap 2 = zeros((n 3, n 2, n 1), 'f') inputm. read(snnap 2) plt. imshow(snnap 2[1, : ], extent=[0, 2, -1. 5, 1. 5]) plt. savefig('migr. png') command 2='display migr. png'; os. system(command 2)

Результат миграции