program RWE2D !! . . Riemannian Wavefield Extrapolation program for migration !! . . in 2D generalized coordinates. !! !! The key concept is that wave-equation migration can be performed !! on almost any type of grid called Riemannian coordintes !! This program allows for user to input !! generalized source and receiver wavefields that are recorded on the !! first extrapolation step, and then propagate them into the user-specified !! computational mesh. !! !! . . General interal processing steps !! 1) Input wavefields, coordinate system, Velocity model !! 2) Calculate geometrical coefficients for coordinate system and velocity model !! 3) Propagate wavefields and apply imaging condition !! 4) Output Image as: a) interpolated to Cartesian; and b) In original coordinates !! !! !! . . Sample call in Makefile !! Image.H: rwf.H swf.H !! $(BINDIR)/RWE2D.x rwf="rwf.H" swf="swf.H" vel="vel.H"\ !! rays="rays.H" image="Image.H" Rimage="RayImage.H" \ !! nref=256 verbose=1 \ !! min_region_pct=1. del_dist=0.05 verb=1 \ !! niter_lloyd=25 ntap=10 nloop=4 migapp=1024 forward=1 \ !! norm=1 zmin=-2000 zmax=6000 xmin=0 xmax=20000 \ !! dxx=10 dzz=10 nsx=2 nsz=2 !! !! where requisite files: !! !! 1) "rwf" tag is the input receiver wavefield in nx-nw-ns format !! 2) "swf" tag is the input source wavefield in nx-nw-ns format !! 3) "vel" tag is the input velocity model in nz-nx format !! 4) "rays" tag is the input coordinate system in nz-nx-2 format !! where nz-nx-1 is the x-coordinate field and nz- !! nx-2 the z-coordinate field !! 5) "image" tag is the Cartesian output imagein cnz-cnx-ns format !! 6) "Riamge" tag is the Ray based ouput image in nz-nx-ns formate !! !! Required parameters: !! !! 1) nref=256 is the starting number of points for calculating reference velocities !! 2) verbose=1 is the level of verbosity for the RWE2D migration code !! 3) min_region=.1 controls how close a true coefficient is to a reference value !! before it is thrown away (in %) !! 4) del_dist=0.05 is the minimum distance between reference parameters (in %) !! 5) verb=1 is the verbage level of the Lloyds code !! 6) niter_lloyd=25 the number of iterations of Lloyds code !! 7) nloop=4 number of smoothing iterations for merging wavefields !! from different references !! 8) migapp=1024 migration aperature split evenly about shotpoint !! 9) forward=1 Forward/Backscattering switch (0/1) !! 10) norm=1 Normalize by Jacobian of matrix !! 11/12) zmin/zmax min/max output Cartesian image z-locations !! 13/14) xmin/xmax min/max output Cartesian image x-locations !! 15/16) dzz/dxx grid spacing in output Cartesian image !! !! Written by Jeff Shragge !! May 15/2006 !! contact: jshragge@gmail.com !! !!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Required modules use sep !! Basic SEPlib I/O use util_ !! handles axis I/O and aux file I/O use RWE2D_wemig !! Migration driver module use C2R2D !! Ray-to-Cartesian interpolation use geometry2D !! Calculate geometry coefficients for dispersion relation implicit none type(myaxis) :: ax,az,aw,ar,as !! axes types integer :: ix,iz,ii,iw,ishot(1),iis,mxl,mxr,oldn integer :: nref,nsx,nsz,norm,cnx,cnz,migapp real :: forward,xmin,zmin,xmax,zmax,dxx,dzz real, allocatable,dimension(:,:,:):: raysin !! Input rays real, allocatable,dimension(:,:) :: Vel,gnorm !! Velocity model, norm matrix complex, pointer ,dimension(:,: ):: rwf, swf !! receiver and source wavefiels real, allocatable,dimension(:,:) :: CCimg, Pimg!! Cartesian/ray coordinate images logical :: verbose,kin character(len=128) :: name !! Geometry calculations integer, allocatable :: mask(:,:),cutmask(:,:) real, allocatable :: refs(:,:,:), fields(:,:,:),cutfields(:,:,:) call sep_init(SOURCE) !! Routines to enable parallel processing with !! SEPlib "Parallel" command call getch_add_string("head=/dev/null") call getch_add_string("noheader=y") call set_no_putch() call sep_begin_prog() name="rays" !! Read in input dimensions from the "rwf" tag !! Determines the x grid parameters call from_aux("rwf","n1",ax%n) call from_aux("rwf","d1",ax%d) call from_aux("rwf","o1",ax%o) !! Determine the frequencies parameters call from_aux("rwf","n2",aw%n) call from_aux("rwf","o2",aw%o) call from_aux("rwf","d2",aw%d) !! Determine the number of shot parameters call from_aux("rwf","n3",as%n) call from_aux("rwf","d3",as%d) call from_aux("rwf","o3",as%o) !! Read in input dimensions from the "rays" tag !! Determines the grid for extrapolation direction call from_aux("rays","n1",az%n) call from_aux("rays","d1",az%d) call from_aux("rays","o1",az%o) !! Parameter input section. !! Call from command line, Makefile or a "par=pars.P" file call from_param("migapp",migapp,ax%n) call from_param("forward",forward,0.) !! Forward scattering option call from_param("nref",nref,256) !! Max # reference velocities call from_param("verbose",verbose,.false.) !! verbosity call from_param("kinematic", kin,.true.) !! Kinematic approximation call from_param("norm",norm,1) !! Whether (1) or not (0) to normalize by gnorm write(0,*) 'DEPTH',az%n,az%o,az%d write(0,*) 'X1 ',ax%n,ax%o,ax%d write(0,*) 'SHOTS',as%n,as%o,as%d write(0,*) 'FREQ ',aw%n,aw%o,aw%d !! Define output Cartesian grid call from_param("xmin",xmin,ax%o) call from_param("zmin",zmin,az%o) call from_param("xmax",xmax,(ax%n-1)*ax%d+ax%o) call from_param("zmax",zmax,(az%n-1)*az%d+az%o) call from_param("dxx",dxx,ax%d) call from_param("dzz",dzz,az%d) call from_param("nsx",nsx,3) call from_param("nsz",nsz,3) cnx = floor((xmax - xmin)/dxx)+1 cnz = floor((zmax - zmin)/dzz)+1 !! Create output header file associated with "image" tag !! . . Output Image in Cartesian coordinates call auxputch("n1","i",cnz ,"image") call auxputch("n2","i",cnx ,"image") call auxputch("n3","i",as%n,"image") call auxputch("d1","r",dzz ,"image") call auxputch("d2","r",dxx ,"image") call auxputch("d3","r",as%d,"image") call auxputch("o1","r",zmin,"image") call auxputch("o2","r",xmin,"image") call auxputch("o3","r",as%o,"image") call auxputch("esize","i",4,"image") !! Create output header file associated with "image" tag !! . . Output Image in RWE coordinates call auxputch("n1","i",az%n,"Rimage") call auxputch("n2","i",ax%n,"Rimage") call auxputch("n3","i",as%n,"Rimage") call auxputch("d1","r",az%d,"Rimage") call auxputch("d2","r",ax%d,"Rimage") call auxputch("d3","r",as%d,"Rimage") call auxputch("o1","r",az%o,"Rimage") call auxputch("o2","r",ax%o,"Rimage") call auxputch("o3","r",as%o,"Rimage") call auxputch("esize","i",4,"Rimage") !! Allocate required arrays allocate( rwf(migapp,aw%n), swf(migapp,aw%n) ) allocate( raysin(az%n,ax%n,2),Vel(az%n,ax%n),gnorm(az%n,ax%n) ) allocate( mask(ax%n ,az%n), refs(nref,4,az%n),fields(ax%n,4,az%n)) allocate( cutmask(migapp,az%n), cutfields(migapp,4,az%n) ) allocate( CCimg(cnz,cnx) , Pimg(az%n,migapp) ) !! Zero all arrays rwf=0.; swf=0.; gnorm = 0.; raysin=0 Vel=0.; mask=0; refs=0.; fields=0. cutmask=0.; cutfields=0.; CCimg=0.; Pimg=0. !! Read in rays file call sep_read(raysin,"rays") !! Read in Velocity file call sep_read(Vel ,"vel ") write(0,*) 'VELOCITY MIN/MAX', minval(Vel),maxval(Vel) !!---------------------------------------------------------------------------------- !! !! GEOMETRY SECTION !! !! . . Initialize module file call geometry_init(az%n,ax%n,az%d,ax%d,kin,nref,verbose) !! . . Calculate the appropriate coefficients !! . . INPUT: !! 1) raysin - the coordinate mesh used to image !! 2) Vel - velocity model !! !! . . OUTPUT: !! 1) fields - the coefficient fields for entire coordinate system (used in SSF) !! 2) refs - a sparse representation of 'fields"' used in the PSPI program !! 3) mask - an integer matrix listing the closest refs to each field value !! 4) gnorm - weight matrix for normalization call geometry_calc(raysin,Vel,fields,refs,mask,gnorm) call geometry_close() deallocate( Vel ) write(0,*) 'GNORM',minval(gnorm),maxval(gnorm) !!------------------------------------------------------------------------------- !! !! MIGRATION SECTION !! !! . . Initialize module and image space oldn=ax%n ax%n=migapp !! . . Initialize the wave-equation migration driver call wemig_init(ax,az,aw,ar,forward,kin,nref,verbose) !! Initialize the RWE to Cartesian interpolation routine call C2R_init(nsx,nsz,norm,dxx,dzz,cnx,cnz,migapp,az%n,xmin,zmin) do ii=1,as%n !! . . Loop over shot points swf = 0.; rwf = 0. call sep_read(rwf ,"rwf ") !! Read in only current receiver wavefield call sep_read(swf ,"swf ") !! Read in only currrent source wavefield Pimg = 0.; CCimg = 0. !! zero image !! Find shot location to determine migration aperature !! Not important when aperature is full width ishot = maxloc(cabs(swf(:,1))); iis=ishot(1) mxl = max( 1,iis-migapp/2) mxr = min(oldn,iis+migapp/2) if (mxl .eq. 1 ) mxr=migapp if (mxr .eq. oldn) mxl=oldn-migapp+1 write(0,*) 'ISHOT',iis,mxl,mxr write(0,*) 'RWF MIN/MAX', minval(abs(rwf)),maxval(abs(rwf)) write(0,*) 'SWF MIN/MAX', minval(abs(swf)),maxval(abs(swf)) !! Window fields and mask by aperture cutfields = fields(mxl:mxr,:,:) cutmask = mask (mxl:mxr,: ) !! Call to migration subroutine !! INPUT !! 1) rwf - receiver wavefield !! 2) swf - source wavefield !! 4) ii - the shot number !! 5) cutfields - the values of coefficients (for SSF) !! 6) refs - reference values for (PSPI/SSF) !! 7) cutmask - integer for closest reference value !! !! INPUT/OUTPUT !! 3) Pimg - the RWE-based image call wemig(rwf,swf,Pimg,ii,cutfields,refs,cutmask ) !! Write out RWE-based image field call sep_write(Pimg,"Rimage") !! RWE-to-Cartesian Interpolation Program write(0,*) "Pimg", maxval(Pimg),minval(Pimg) write(0,*) "Real Rays", maxval(raysin(:,:,1)),minval(raysin(:,:,1)) write(0,*) "Imag Rays", maxval(raysin(:,:,2)),minval(raysin(:,:,2)) write(0,*) "Gnorm", maxval(gnorm),minval(gnorm) call C2R_run(Pimg,raysin(:,mxl:mxr,:),CCimg,gnorm(:,mxl:mxr)) write(0,*) "CCimg",maxval(CCimg),minval(CCimg) !! Write out Cartesian-based image call sep_write(CCimg,"image") end do !! . . End Shot-point loop call sep_end_prog() call exit() end program RWE2D