Home

DENISE - Black Edition User Manual - Christian

1. 62 7 1 1 Acquisition geometry FD 65 7 1 2 Elastic wave propagation in the complex Marmousi model 65 7 1 3 FWI of the complex Marmousi 68 7 1 4 Marmousi II Benchmarks Chapter 1 Introduction The aim of Full Waveform Tomography FWT is to estimate the elastic material parameters in the underground This can be achieved by minimizing the misfit energy between the modelled and field data using a gradient optimization approach Because the FWT uses the full information content of each seismogram structures below the seismic wave length can be resolved This is a tremendous improvement in resolution compared to traveltime tomography Pratt et al 2002 The concept of full waveform tomography was originally developed by Albert Tarantola in the 19805 for the acoustic isotropic elastic and viscoelastic case Tarantola 1984b a 1986 1988 First numerical implementations were real ized at the end of the 1980s Gauthier et al 1986 Mora 1987 Pica et al 1990 but due to limited computational resources the application was restricted to simple 2D synthetic test problems and small near offset datasets At the begining of the 1990s the original time domain formulation was transfered to a robust frequency domain approach Pratt and Worthington 1990 Pratt 1990 With the increasing performance of supercomputers mod
2. djm OK uj Xo t 3 sources The terms only depending on time t and the positions x can be moved infront of the sum over the receivers M a e x t gt dt Gi Xa t x t duj xa t sources y rs x t ws av cia x 6 t 3 28 sources M t dm Ef av Ci x t u xa t sources Defining the wavefield pT dt t x 0 t 3 29 CHAPTER 3 THE ADJOINT PROBLEM 25 Eqs 3 28 can be written as gt x t Vi sources aw M Js dte t 61m a 3 30 sources QV y ra dt im sources The wavefield V is generated by propagating the residual data du from the receiver positions backwards in time through the elastic medium To obtain a more symmetric expression for the density gradient let us integrate the density gradient in 3 30 by parts M le x sources gt er wen Seien sources 3 31 According to Eqs 2 12 the field uj x t satisfies initial conditions of rest 0 0 and 0 0 0 The field V x t satisfies final conditions of rest V x 0 Therefore M de x 0 x an x t 20 RATS t 3 32 sources sources Writing out the implicit sums in the gradients of the Lam parameters and dy in Eqs 3 30 _
3. 32 How to find an optimum 1 3 3 Calculation of the gradient direction E a Be SG EA ei ta 3 4 Gradients for different model 2 2 2 3 5 Gradients for the stress velocity 3 6 Estimation of an optimum step length un 3 7 Nonlinear Conjugate Gradient 3 8 JT he elas c 5 28 wa Pe Bet RU Tea RN Bou k 4 Source Wavelet Inversion 5 Getting Started 2 1 JREQUIFEMENIS a Ano ode Are ae RISE Pendet EAM eh a ra ak Boe deep rx bou epp 5 1 2 How to run DENISE on the NEC Linuxcluster at RZKiel 25 2 Installation Furl ua SIE d A S RH eps uie 5 3 Compilation of DENISE 4 222 22 ep eder 5 4 Running the program o scs e RES SUE den e 5 5 Postprocessing ke mobs CR Re Br eue a Cu RUE 6 Definition of parameters for the modeling and inversion code 6 1 Input file with fixed parameters 6 2 Workflow file with variable inversion parameters FWI workflowinp cess 18 18 19 20 26 26 27 29 32 33 3 7 Example 1 the elastic Marmousi2 model 62 7 1 complex Marmousi2 model
4. 3 23 leads to 7 u u x t t Git tox mda Ot d ETA xS OA 3 24 V t gg 56x Rt rd 3 24 0 T Gij n dt Ox 2 x t x t en A t CHAPTER 3 THE ADJOINT PROBLEM 24 Utilization of Eq 3 24 to solve the forward problem is known as the Born approximation In waveform tomography the Born approximation is not used to solve the forward problem Instead the full elastic wave equation is solved Equation 3 24 has the same form as the desired expression for the forward problem Eqs 3 13 u dV dm 3 25 Therefore the Frech t kernels O for the individual material parameters can be identified as Ou T u dt Gi x tx 0 0 Ou i2 W RE NR DA 18 jx RE t J 3 26 Ou 36 a i dt 3 x t x yt Vem A t 0 By definition the adjoint of the operator 3 25 can be written as 2 5 are sources a t 3 27 Because a linear operator and its transpose have the same kernels the only difference arise in the variables of sum integration which are complementary Inserting the integral kernels 3 26 in Eq 3 27 yields M dt Gi t Uy sources y UN is En a Ci Kt t 1m X 8 1 01m 9 uj xo t sources gt ri av Gi x e 6
5. starting model 2000 3000 4000 500 1000 1500 2000 2500 1500 2000 2500 P wave velocity m s S wave velocity m s Density kg m 72 Figure 7 6 Depth profiles at xp 3 5 km top and xp2 6 4 km bottom of the starting model FWT result are compared with the true model for the Marmousi II model P wave velocity left S wave velocity center and density right CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL 73 Start Model FWT Result 200 300 channel channel Initial Residuals Final Residuals 200 channel channel Evolution of the Residual Energy True Model o o gt a Normalized Residual Energy e 100 200 300 Iteration Step No channel Figure 7 7 Seismic sections shot 50 y component for the Marmousi II model a starting model b FWT result c initial residuals d final residuals f true model and e evolution of the residual energy CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL 74 PCG depth V m s 2500 12000 2 1500 a 11000 500 EPRECOND 3 m s 2500 E 2000 2 1500 a 1000 500 LBFGS EPRECOND 3 V m s 2500 E 2000 2 2 1500 1000 500 Distance km Figure 7 8 Influence of different Hessian approximations on the Vs inversion results A PCG optimization with a simple linear scaling of the gradient with d
6. are the Lam parameters and v is the particle velocity vector The locations of the receivers may either be specified in a separate file REC_FILE or in this parameter file If READREC 1 receiver locations are read from the ASCII file REC_FILE Each line contains the coordinates of one receiver the first two number in each line indicate the horizontal x and the vertical y coordinate of each receiver position To give an example of a receiver file the following 3 lines specify 3 receivers located at constant depth 2106 0 m However the receiver coordinates change in x direction starting at 540 m and therefore lining up along the x axis 540 0 2106 0 1080 0 2106 0 1620 0 210640 These receiver coordinates in REC_FILE are shifted by REFREC 1 REFREC 2 into the and y direction respectively This allows for completely moving the receiver spread without modifying REC_FILE This may be useful for the simulation of moving profiles in reflection seismics If READREC 0 the receiver locations must be specified in the parameter file In this case it is assumed that the receivers are located along a straight line The first receiver position is defined by XREC1 YREC1 and the last receiver position by YREC1 see Figure 6 1 The spacing between receivers is NGEOPH grid points Receivers are always located on full grid indices i e a receiver that is located between two grid points will be shifted by the FD program to the cl
7. NORMALIZE 2 Normalize for each trace the maximum of the field data relative to the maximum of the syn thetic data Normalize the estimated source wavelet to its maximum amplitude CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 60 Offset Windowing In some cases the application of an offset window be useful to achieve a layer stripping update of the model parameters from top to bottom e OFFSET MUTE 1 mutes all traces with an offset larger than OFFSETC far offset mute OFFSET MUTE 2 mutes all traces with an offset smaller than OFFSETC near offset mute Density model update restriction Because changes of the density model are in most cases smaller than velocity changes the step length for the density update can be systematically reduced by a factor SCALERHO see Eq 3 49 24 E E BIZ Z idi NIS E cms gt lt a 8 5 2 2 0 01 1 0 0 1 71612 1e2 0 0 1 75 0 600 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1110 0 2 9 6 2 le2 0 0 1 75 0 600 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 00 9 01612 le2 0 0 1 75 0 600 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 0 0 1 7 6 2 1e2 0 0 2 55 0 600 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0
8. INCR X or y direction by REC INCR Y In this example the re ceivers are moved in x direction by 80 m for each shot increment More complex acquisition geometries e g streamer and OBS can be implemented by assuming that only the first N_STREAMER receivers contribute to the streamer and are subsequently moved All other receiver positions in the receiver file REC_FILE are fixed Seismograms Seismograms samplingrate in timesteps NDT 1 data format_ SU 1 ASCII 2 BINARY 3 output files for seismograms 1 particle velocities if SEISMO 1 or SEISMO 4 filename_for_Vx_ SEIS_FILE_VX su DENISE_US_x su filename_for_Vy_ SEIS_FILE_VY su DENISE_US_y su curl and div of wavefield if SEISMO 3 or SEISMO 4 filename_for_curl_ SEIS_FILE_CURL su 2layer_rot su filename_for_div_ SEIS_FILE_DIV su 2layer_div su pressure field hydrophones if SEISMO 2 or SEISMO 4 filename_for_pressure_ SEIS_FILE_P su full_wave_forward_p su p If SEISMO gt O seismograms recorded at the receiver positions are written to the corresponding output files The sampling rate of the seismograms is NDT DT seconds In case of a small time step interval and a high number of time steps it might be useful to choose a high NDT in order to avoid a unnecessary detailed sampling of the seismograms and consequently large files of seismogram data
9. V fena fi 1 2 Balf raya fi 3 2 8 1 8 27 ore dh 8 3 728 31 x The weights y can be calculated by the following approach The factor in front of the partial derivative on the LHS of Eq 2 10 should equal 1 therefore 61 382 1 on the RHS of Eq 2 10 should vanish 1 3 The coefficients in front of As B1 2783 0 The weights can be estimated by solving the matrix equation 1 3 41 1 27 bo 0 The resulting coefficients 81 9 8 and G2 1 24 Therefore the 4th order backward and forward operators are Of 1 a 21 fina fi Balfiva 1 forward operator X li 1 2 2 11 f 1 e ai A fii Balfiz fi 2 backward operator X 111 2 The coefficients 6 in the FD operator are called Taylor coefficients The accuracy of higher order FD operators can be improved by seeking for FD coefficients that approximate the first derivative in a certain frequency range Holberg 1987 These numerically optimized coefficients are called Holberg coefficients 2 2 3 Initial and Boundary Conditions To find a unique solution of the problem initial and boundary conditions have to be defined The initial conditions for the elastic forward problem are uj x t 0 x t 2 12 ee O for all x V at 0 For the geophysical application two types of boundary conditions are very important 1 Horizonta
10. 2 2 2 m 9 Karlsruher Institut f r Technologie DENISE User Manual Christian Albrechts Universit t Kiel Germany and Karlsruher Institut f r Technologie KIT Germany Version 1 1 April 24 2015 Authors The DENISE code was first developed by Daniel K hn Denise De Nil and Kurzmann at the Christian Albrechts Universit t Kiel and TU Bergakademie Freiberg Germany from 2005 to 2009 The forward code is based on the viscoelastic FD code fdveps now SOFI2D by Bohlen 2002 Different external libraries for timedomain filtering are used The copyright of the source codes are held by different persons cseife c cseife h lib_stfinv lib_aff lib_ fourier Copyright c 2005 by Thomas Forbriger BFO Schiltach cseife_deriv c cseife_gauss c cseife_rekfl c cseife_rfk c and cseife_tides c Copyright c 1984 by Erhard Wielandt This algorithm was part of seife f A current version of seife f can be obtained from http www software for seismometry de The Matlab implementation of a few SU routines mainly used to read and write SU files in data pre processing are Copyright C 2008 Signal Analysis and Imaging Group For more information http www geo phys ualberta ca saig SeismicLab Author M D Sacchi Since then it has been further developed maintained and improved by contributions from in alphabetical order Lisa Groos Sven Heider Martin Sch fer Linbin Zhang add other devel
11. Sources 5 DA du sources 3 40 2Vp Xl HET sources 2pV_pdA The gradients for Vs and are calculated in a similar way so the gradients in terms of seismic velocities can be written as OVp 2pV p l 4pV50 20V s u 3 41 dpvear Vi 2V2 0A 3 5 Gradients for the stress velocity code The elastic forward problem and back propagation of the data residuals is solved by using a time domain stress velocity finite difference FD code Therefore the displacements in 3 36 have to be replaced by stresses and particle velocities v e g Shipp and Singh 2002 ee sources Se fa 7 an 4 2 3 42 Jas ow were and gt are the stresses of the forward and backpropagated wavefield respectively The displacements Y in the density gradient are calculated from the particle velocities by numerical integration CHAPTER 3 THE ADJOINT PROBLEM 27 3 6 Estimation of an optimum step length un The choice of the step length in Eq 3 10 is crucial for the convergence of the steepest gradient optimization method I demonstrate this using a very familiar test problem for optimization routines the Rosenbrock function with two unkown parameters Rosenbrock 1960 Fig 3 4 amp x y 1 x 100 y 2 3 43 The aim is to find the minimum of this function located at the point 1 1 which
12. Optimization method ses Optimization Method gradient method PCG 1 LBFGS 2 GRADIENT 3 _ GRAD_METHOD 2 save NLBFGS updates during LBFGS optimization 20 During FWI the misfit function can be minimized by different optimization methods Currently a preconditioned conjugate gradient PCG the quasi Newton method limited memory Broyden Fletcher Goldfarb Shanno 1 BFGS method see e g Nocedal and Wright 2006 and a simple gradient method can be applied When using the 1 BFGS optimization the last NLBFGS updates are stored Reduce inversion grid Reduce inversion grid use_only_every_DTINV_time_sample_for_gradient_calculation_ DTINV 3 To reduce the memory requirements during an inversion one can define that only every DTINV time sample is used for the calculation of the gradients To set this parameter appropriately one has to keep in mind the Nyquist criterion to avoid aliasing effects CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 56 Step length estimation Step length estimation maximum_model_change_of_maximum_model_value_ EPS_SCALE maximum number of attemps to find a step length STEPMAX SCALEFAC 4 0 testshots for step length est TESTSHOT_START TESTSHOT_END TESTSHOT_INCR 1 o ds O For the step length estimation a parabolic line search method proposed by Sourbier et al 20
13. lcseif 3333333 The program snapmerge that is used to merge the snapshots see below can be compiled with make snapmerge Since this is not a MPI program no MPI functions are called the MPI libraries are not required and any standard compiler like gec and cc can be used to compile this program The executables denise and snapmerge are located the directory bin You can also compile DENISE and snapmerge from the par directory using the shell script compile DENISE sh 5 4 Running the program Each DENISE run reads the required parameters from the parameter files par DENISE inp and par FWI_workflow inp A detailed description of the parameters can be found in chapter 6 The command to start a simulation on 8 processor with the lowest priority of 19 in order to allow working on the a workstation while running a simulation is as follows Please note that we assume you have navigated to the folder DENISE par and all parameter files are located in this directory mpirun 8 nice 19 bin denise DENISE inp FWI_workflow inp If you use LAMMPI the command lamboot v lamhost must be executed on node 0 which is the PE where par lamhosts is the file containing IP addresses of all computing nodes It if often useful to save the standard output of the program for later reference The screen output may be saved to DENISE out using mpirun 8 nice 19 bin denise DENISE inp FWI_work
14. Keep in mind that the application of FWT requires NDT 1 Possible output formats of the seismograms are SU ASCII and BINARY It is recommended to use SU format for saving the seismograms The main advantage of this format is that the time step interval NDT DT and the acquisition geometry shot and receiver locations are stored in the corresponding SU header words Also additional header words like offset are set by DENISE This format thus facilitates a further visualization and processing of the synthetic seismograms Note however that SU cannot handle sampling rates smaller than 1 0e 6 seconds and the number of samples is limited to about 32 000 In such cases you should increase the sampling rate by increasing NDT If this is impossible for example because the Nyquist criterion is violated you must choose a different output format ASCII or binary CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 53 Monitoring the simulation each PE is printing log information to LOG_FILE MYID log file_for_information_about_progress_of_program_ LOG_FILE log 2layer log info of processing element zero to stdout yes 1 no 0 LOG 1 DENISE can output a lot of useful information about the modeling parameters and the status of the modeling process etc The major part of this information is output by PE 0 If LOG 1 PE 0 writes this info to stdout i e on the screen of your shell This is generally recommended to m
15. _SOURCES 1 only the final gradients that means the gradients obtained by the summation of the gradients of each shots are multiplied with a taper that decreases the gradients at all shot positions Therefore one looses the update information at the source positions To avoid this one can use SWS_TAPER_CIRCULAR_PER_SHOT 1 In this case the gradients of the single shots are preconditioned with a window that only damps the gradient at the current shot position This is done before the summation of all gradients to keep model update information at the shot positions The actual tapers are generated by the function taper grad c and taper grad shot c respectively The circular taper around the source positions decrease from a value of one at the edge of the taper to a value of zero at the source position The damping shape can be defined by an error function SRTSHAPE 1 or a log function SRTSHAPE 2 The radius of the taper is defined in meter by SRTRADIUS Note that this radius must be at least 6 gridpoints With the parameter FILTSIZE one can extend the region where the taper is zero around the source The taper is set to zero in a square region of 2 FILTSIZE 1 times 2 FILTSIZE 1 gridpoints All preconditioning matrices applied to the gradients are saved in the par directory with the file names taper vert bin taper coeff horz bin and taper coeff sources bin To apply an externally defined taper on the gradients in DENISE the parameter SWS TAPER
16. 0 0 01 1110 0 2 91612 1e2 0 0 2 55 0 600 600 0 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 00 9 0 6 2 1e2 0 0 2 55 0 600 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 0 0 9 0 6 2 1e2 0 0 3 55 0 600 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1110 0 9 0 6 2 1e2 0 0 5 15 0 600 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 0 0 1 7 6 2 1 2 0 0 1 75 600 0 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 00 2 9 6 2 1e2 0 0 1 75 600 0 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 10 0 9 0 6 2 1e2 0 0 1 75 600 0 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 0 0 1 71612 1e2 0 0 2 55 600 0 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 10 0 2 9 6 2 1e2 0 0 2 55 600 0 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 10 0 9 0 6 2 1e2 0 0 2 55 600 0 600 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 0 0 9 0 6 2 1e2 0 0 3 55 600 0 600 O 160 160 3 2 O 2000 0 le 1 0 O 10 0 0 0 0 01 1 10 0 9 0 6 2 1e2 0 0 5 15 600 0 600 0 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01
17. 1 0 0 1 7 6 2 le2 0 0 1 75 600 600 0 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 10 0 2 9 6 2 1e2 0 0 1 75 600 600 0 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 10 0 9 0 6 2 1e2 0 0 1 75 600 600 0 0 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 10 0 1 7 6 2 1e2 0 0 2 55 600 600 0 0 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 10 0 2 9 6 2 1e2 0 0 2 55 600 600 0 O 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 10 0 9 0 6 2 1e2 0 0 2 55 600 600 0 0 160 160 3 2 O 2000 0 le 1 O O 10 0 0 0 0 01 1 10 0 9 0 6 2 1e2 0 0 3 55 600 600 0 160 160 3 2 O 2000 0 le 1 0 O 10 0 0 0 0 01 1 10 0 9 0 6 2 1e2 0 0 5 15 600 600 0 0 160 160 3 2 O 2000 0 le 1 0 O 10 0 0 0 Table 6 1 Example of a complex hierachical multiparameter FWI work flow definition modified after Kurzmann 2012 The colors represent different parameter groups HGOO NOISAHANI ONITHCOW FHL 5 HO NOLLINIHHG 9 AHLAVHO Chapter 7 Example 1 the elastic Marmousi2 model Developed in the 1990s by the French Petroleum Institute IFP Versteeg 1994 the Marmousi model is a widel
18. 1 seismic impedances Zp Zs rho INVMAT 1 2 or Lam parameters A u INVMATI 3 If models are read from binary files appropriate file extensions are required for the different models see section 6 1 Depending on the data different components of the seismic sections can be backpropagated For two il 2 CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 54 component data x and y component set QUELLTYPB 1 y component only QUELLTYPB 2 x component only QUELLTYPB 3 and pressure component QUELLTYPB 4 To estimate the step length u for the model update requires the calculation of the misfit function and therefore the solution of the forward problem for all shots involved However to save computation time the calculation of the misfit function can be restricted to certain shots using only shots between TESTSHOT_START and TESTSHOT_END with a shot increment TESTSHOT_INCR The parameters TAPER TAPERLENGTH and INVTYPE are debug parameters and should not be changed The parameters JACOBIAN and GRADTI to GRADT4 are explained in section 6 1 and section 6 1 respectively Definition of the gradient taper geometry SAR Definition of gradient taper geometry Vertical taper apply vertical taper yes 1 SWS TAPER GRAD VERT 0 Horizontal taper apply horizontal taper yes 1 SWS TAPER GRAD HOR 0 exponent
19. CHAPTER 5 GETTING STARTED 37 N name of batch job o file name for standard output 4 the requested batch class An example for a job file can be found in the DENISE jobs directory The job can be submitted with sungwXXX rzcluster qsub DENISE job With sungwXXX rzcluster qstat you can check the status of your Jobs and with sungwXXX rzcluster 15 adel job id you can cancel a submitted or running job where lt job d gt denotes the number at the first column of the status information f e sungwXXX nesh fe2 jobs qstat RequestID ReqNam UserName Queu Pri SIT E Memory CEU Elapse M Job 459654 ace ssio DENISE sungwXXX clmedium 0 RUN 13 89G 1077746 75 69303 Y Y X 470371 ace ssio SAVA sungwXXX clmedium 0 QUE 0 00B o 00 D X Y X sungwXXX8nesh fe2 jobs qdel 459654 would kill the first job in the queue For further information I again refer to the homepage of the RZ Kiel https www rz uni kiel de hpc nec cluster html 5 2 Installation DENISE consists of four different packages The source code A collection of benchmark models and pre postprocessing tools Matlab scripts for data pre processing The manuals for the different code versions Start with unpacking the source code package e g by tar zxvf DENISE tgz and changing to the directory DENISE cd DENISE you will find different subdirectories svn The software is updated using the vers
20. E LED Ad sources 5 3 33 Ef ELLE and neglecting all wavefield components and derivatives in z direction leads to OW OW UE er 3 34 Y fa Once OUy Be OV OW Sources id Ox Using the definition of the strain tensor e we get 5 i OV OW sources 3 35 3 dt i 0 OW 0 QV sources Finally the gradients for the Lam parameters A and the density can be written as 2 Ouy OV es OW sources Ou Ouy OV OV du OV du AV t x y x ES y 2 y y p fe E Ox Ox 3 36 sources Ou OV Qu OV 2 x a sources CHAPTER 3 THE ADJOINT PROBLEM 26 3 4 Gradients for different model parametrizations The gradients in terms of other material parameters Mpew can be calculated by applying the chain rule on the Frech t kernel in the adjoint problem Eq 3 27 Om X dt gt gt u 3 37 SOUTCES Using the relationships between P wave velocity Vp S wave velocity Vs the Lam parameters and density p x42 V voy 3 38 p p pV2 29V2 pV 3 39 or The gradient for V can be written as 9u gt PE V OnOV 7 Op x 25
21. R Pratt Full waveform tomography for lithospheric imaging results from a blind test in a realistic crustal model Geophys J Int 168 133 151 2007 Brossier Imagerie sismique deux dimensions des visco lastiques par inversion des formes d ondes developpements methodologiques et applications PhD thesis Universite de Nice Sophia Antipolis 2009 R Brossier S Operto and J Virieux Two dimensional seismic imaging of the Valhall model from synthetic OBC data by frequency domain elastic full waveform inversion chapter 461 pages 2293 2297 2009 Y Choi and T Alkhalifah Application of multi source waveform inversion to marine streamer data using the global correlation norm Geophysical Prospecting 60 748 758 2012 Y Choi D Min and C Shin Frequency domain elastic full waveform inversion using the new pseudo hessian matrix Experience of elastic marmousi 2 synthetic data Bull Seis Soc Am 98 2402 2415 2008a Y Choi D Min and C Shin Two dimensional waveform inversion of multi component data in acoustic elastic coupled media Geophysical Prospecting 56 863 881 2008b R Courant K Friedrichs and H Lewy ber die partiellen Differenzengleichungen der mathematischen Physik Mathematische Annalen 100 32 74 1928 Courant Friedrichs and Lewy On the partial difference equations of mathematical physics IBM Journal pages 215 234 March 1967 H Denli and L Huang Double di
22. Spike 6 Klauder 7 QUELLART 6 point source explosive 1 force in x 2 force in y 23 rotated force 4 QUELLTYP 3 angle of rotated directed source relative to y direction in degree ANGLE 0 0 depth of plane wave excitation no 0 in meter PLANE WAVE DEPTH 0 0 dip of plane wave from vertical in degrees PHI 0 0 SIGNAL FILE wavelet wavelet US opt dat duration of Klauder wavelet in seconds TS 8 0 read source positions from SOURCE FILE yes 1 SRCREC 1 SOURCE FILE source source dat run multiple shots defined in SOURCE FILE yes 1 RUN MULTIPLE SHOTS 1 corner frequency of highpass filtered spike FC SPIKE 1 5 0 corner frequency of lowpass filtered spike FC SPIKE 2 15 0 order of Butterworth filter ORDER SPIKE 5 6 built in wavelets of the seismic source are available The corresponding time functions are defined in src wavelet c You may modify the time functions in this file and recompile to include your own analytical wavelet or to modify the shape of the built in wavelets Ricker wavelet t 1 5 f 1 5 f t 1 2 with ait LU 4 6 2 Fuchs M ller wavelet fm t sin 27 t ta fe 0 5sin 4r t ta fe if t ta ta 1 fc else fm t 0 6 3 sin wavelet s3 t 0 75mf sin m t ta f if t tata 4 1 fc else s3 t 0 6 4 First Gaussian derivative wavelet gd t 2z 2 t ta exp 7 f2 t ta 6 5 Ban
23. elastic viscoelastic wave propagation pages 1277 1280 Houston Texas 1995 H Rosenbrock An automatic method for finding the greatest or least value of a function The Computer Journal 3 175 184 1960 T Sears S Singh and P Barton Elastic full waveform inversion of multi component OBC seismic data Geophysical Prospecting 56 6 843 862 2008 D Sheen K Tuncay C Baag and P Ortoleva Time domain Gauss Newton seismic waveform inversion in elastic media Geophys J Int 167 1373 1384 2006 C Shin and W Ha A comparison between the behavior of objective functions for waveform inversion in the frequency and Laplace domains Geophysics 73 5 VE119 VE133 2008 C Shin K Yoon K Marfurt K Park D Yang H Lim S Chung and S Shin Efficient calculation of a partial derivative wavefield using reciprocity for seismic imaging and inversion Geophysics 6 1856 1863 2001 R Shipp and S Singh Two dimensional full wavefield inversion of wide aperture marine seismis streamer data Geophys J Int 151 325 344 2002 F Sourbier S Operto J Virieux P Amestoy and J L Excellent FWT2D a massively parallel program for frequency domain full waveform tomography of wide aperture seismic data part 1 algorithm Computer amp Geosciences 35 487 496 2009a F Sourbier S Operto J Virieux P Amestoy and J L Excellent FWT2D a massively parallel program for frequency domain full waveform tomography of
24. from the data to the model space Eq 3 14 is equal to the gradient of the residual energy Eq 3 12 m Y fa Y FE Y fa D xim 3 15 Om if the perturbation of the data space is interpretated as data residuals So the approach to estimate the gradient direction E m can be split into 3 parts 1 Find a solution to the forward problem L m 2 Identify the Frech t kernels Ou Om 3 Use the property that a linear operator L and it s adjoint L have the same kernels and calculate the gradient direction by using OE dm u This is a very general approach Now we apply this approach to the equations of motion for an elastic medium The following derivation is much easier when assuming a general elastic medium first and introduce the isotropy later on Therefore the elastic forward problem Eqs 2 3 can be written as u P Ot OX 7ij Cijki ki Tij 1 Qu Ou 3 16 2 Ox initial and boundary conditions 0ij fi ij where p denotes the density oj the stress tensor ei the strain tensor c the stiffness tensor fi Tij source terms for volume and surface forces respectively In the next step every parameter and variable in the elastic wave equation is perturbated by a first order perturbation as shown in Fig 3 3 dui p pt op Tij 05 3 17 Cii Cijkl ij ij CHAP
25. inp CHAPTER 5 GETTING STARTED 42 Depending on the model size the merge process may take a few seconds or hours The output should read like this pressure files snap test bin p writing merged snapshot file to snap test bin p Opening snapshot files snap test bin p finished Copying finished Use xmovie n1 100 n2 100 snap test bin p loop 1 labell Y label2 X title g to play movie Chapter 6 Definition of parameters for the modeling and inversion code The geometry of the FD grid and all parameters for the wavefield simulation and inversion have to be defined in a parameter file which we name in this case DENISE par DENISE inp Parameters changing during the waveform inversion are defined in a separate file which we name in this case DENISE par FWI_workflow inp This allows the flexible combination of different inversion parameters and therefore an implementation of complex FWI workflows In the following we will explain every input parameter section in detail 6 1 Input file with fixed parameters DENISE inp Most lines in the parameter file are organized as follows description_of_parameter_ VARNAME _ switches parameter value where VARNAME denotes the name of the global variable in which the value is saved in all functions of the pro gram The possible values are described in switches A comment line is indicated by a on the very first position of a
26. of depth scaling for preconditioning EXP TAPER GRAD HOR 2 0 Circular taper around all sources not at receiver positions apply cylindrical taper yes 1 SWS TAPER GRAD SOURCES 0 apply cylindrical taper per shot yes 1 SWS TAPER CIRCULAR PER SHOT 1 1 function 2 10g function SRTSHAPE 1 radius in m SRTRADIUS 5e 3 gt minimum for SRTRADIUS is 5x5 gridpoints filtsize in gridpoints FILTSIZE 25 read taper from file yes 1 SWS TAPER FILE 0 Different preconditioning operators can be created and applied to the gradients using the function taper_grad c to improve convergence speed depth resolution and define which parts ofthe model should be updated To apply a ver tical or a horizontal taper one has to set the switches SWS_TAPER_GRAD_VERT and SWS_TAPER_GRAD_HOR to 1 respectively The parameters for the vertical and the horizontal window are defined by the input file paramters GRADTI GRADT2 GRADT3 and GRADT4 Please have a look at the function taper_grad c directly to obtain more information about the actual definition of the tapers In case of SWS_TAPER_GRAD_HOR 1 the gradient can also be scaled with the following depth dependent preconditioning operator P depth XP TAPER_GRAD_HOR 6 8 It is also possible to apply cylindrical tapers around the source positions This can be done by setting the switch SWS_TAPER_GRAD_SOURCES or SWS_TAPER_CIRCULAR_PER_SHOT to 1 If one uses SWS_TAPER_GRAD
27. recorded by the streamer in the water column As an example Fig 7 7 f shows the seismic section of the y component for shot 50 Beside the direct wave and a strong reflection from the seafloor numerous small reflection events from the thrust fault system are dominating the seismic section CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL 67 Pressure Pa Pressure Pa 25 3 3 i 2 5 x km x km Pressure Pa Pressure Pa 25 0 5 2 25 x km Pressure Pa x km Pressure Pa 25 3 x km Figure 7 3 Pressure wavefield excited by shot 50 for the elastic Marmousi2 model at 6 different time steps CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL 68 7 1 3 FWI of the complex Marmousi model Due to the far offset acquisition geometry we use seismic velocities as model parameters for the inversion K hn et al 2012 To generate a starting model which accurately describes the long wavelength part of the material parameters the true models m2 V Vs and m1 u are filtered using a spatial 2D Gaussian filter 1 is E 1 x x T Y Msmooth X y 2x2 J dx dy m x no yY jew 222 with a correlation length A 800 0 m As a result all the small scale structures vanish and only the large scale structures are present Fig 7 4 This starting model is comparable with the resolution of a Laplace domain waveform inversion Shin and Ha 2008 Because the Marmousi II model is
28. the wavefield variables at the correct positions the variables are not placed on the same grid points but staggered by half of the spatial grid point distance Virieux 1986 and Levander 1988 figure 2 1 shows the distribution of the material parameters and wavefield variables on the spatial grid X xx yy uy Py 9 uo i 1 4 6 gt u dh y a Y 1 1 G 1 j D lt gt dh Figure 2 1 Grid geometry for a standard staggered grid SSG in Cartesian coordinates as suggested by Virieux 1986 and Levander 1988 To guarantee the stability of the standard staggered grid SSG code Lam parameter and density have to be averaged harmonically and arithmetically Moczo et al 2004 Bohlen and Saenger 2006 respectively n il o un eu oe eee een exlilli Ici 1 2 6 el 1 20 CHAPTER 2 THEORETICAL BACKGROUND 9 The discretization of the linear stress strain relationship in 2 3 at time step leads to the following system of equations for simplicity I skip the time index n gt ee Uxx jlli e dh P E dh tli Elio uy Elf Jl 2 er dh sll 1 i na ux j li 5 2 7 osi lli slo AGILE vedi us I 2 lt gt Gl Well edil 2 lt gt s The discretization of the momentum equation in
29. 00 0 01 s Computing timestep 3000 of 3462 Message from update v printed by PE 0 Updating particle velocities finished real timet 000 3 particle velocity exchange between PEs finished real time 0 00 s Message from update s printed by PE 0 Updating stress components finished real time 0 00 s stress exchange between PEs finished real time 0 00 s total real time for timestep 3000 0 01 s PE 0 is writing 11 seismograms vx to Su DENISE US x sn shotl iti1 PE 0 is writing 11 seismograms vy to Su DENISE US y su shotl itl j Info from main written by PE 0 CPU time of program per PE 17 seconds Total real time of program 18 08 seconds Average times for velocity update 0 003 seconds stress update 0 002 seconds velocity exchange 0 000 seconds stress exchange 0 000 seconds timestep 0 005 seconds 5 5 Postprocessing The wavefield snapshots can be merged using the program snapmerge The program snapmerge is not a MPI program Therefore it can be executed without MPI and the mpirun command You can run snapmerge on any PC since a MPI environment e g LAM is not required You may therefore copy the snapshot outputs of the different nodes to another non MPI computer to merge the files together snapmerge reads the required information from the DENISE parameter file Simply type bash 2 05b DENISE par gt bin snapmerge DENISE
30. 09a b Brossier 2009 and Nocedal and Wright 2006 is implemented For this step length estimation only two further test forward modelings are needed The vector L2t contains the misfit values and the vector epst contains the corresponding step length During the forward modeling of an iteration step the misfit norm of the data residuals is calculated for the shots defined by TESTSHOT_START TESTSHOT_END and TESTSHOT_INC The value L2t 1 then contains the misfit from the forward modeling and the corresponding epst 1 value is 0 0 The step lengths for the different parameters are defined as EPSILON EPS_SCALE m_max grad_max EPSILON epst i m_max grad_max where m_max is the maximum value of the corresponding model parameter in the whole model and grad_max is the maximum absolute value of the gradient For a better definition of the parabola the improved line search is now trying to estimate a steplength epst 2 with L2t 2 lt L2t 1 If the code is not able to find an appropiate steplength using the user defined value EPS_SCALE f e EPS_SCALE 0 01 1 change in terms of m_max grad_max the code divides this steplength by the variable SCALEFAC and calculates the misfit norm again If this search fails after STEPMAX attempts DENISE exits with an error message If the algorithm has found an appropriate value for epst 2 it is trying to estimate a steplength epst 3 with L2t 3 gt L2t 2 by increasing the steplength EPS_SCALE EPS_SCALE
31. 100 05 0 01 1 00 200 6 0 12 00 oo 160 16013 2 0 20000 le 1 O 100 05 Table 7 2 FWI workflow definition for the inversion of the Marmousi model data The colors represent different parameter groups IHHOW CISQOWYAVW OLLSV IH TIdNAVXA 7 69 CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL 70 m s 4000 E 3000 X 2 2000 1000 1 2 3 4 5 6 7 8 9 10 Distance km m s 2500 2000 4 1500 1000 1 2 3 4 5 6 7 8 9 Distance km Density kg m 2500 2000 Depth km 1500 Distance km Figure 7 4 Initial models for the Marmousi II model CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL 71 m s 4000 3000 Depth km 2000 1000 Depth km Distance km 2500 2000 Depth km 1500 Distance km Figure 7 5 Results of the elastic FWT for the Marmousi II model CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL Profile 1 Depth km in N true model FWT result starting model 1500 2000 2500 3000 3500 4000 4500 500 1000 1500 2000 2500 1500 2000 P wave velocity m s S wave velocity m s Density kg m Profile 2 Depth km true model FWT result
32. 2 3 leads to the following system of equations 0 318 056 5 ett alti owbl 57 owbl 51 08 1 1 1 dt 29 w o 2 2 2 Accuracy of FD operators The derivation of the FD operators in the last section was a simple replacement of the partial derivatives by finite differences In the following more systematic approach the first derivative of a variable f at a grid point i is estimated by a Taylor series expansion Jastram 1992 Of 1 2k Ux E qi tit k 1 2 amp 1 2 1 Y k 4 dh 1 2N TH ao aan TO i 1 2 For an operator with length 2N N equations are added with a weight N 1 a 5 6 1 2 fi i2 1 k 1 N N 1 21 1 1 1 1 k z dh l 1 se 1 1 xCD N DA k 1 2 9 O dh N 1 The case N 1 leads to the FD operator derived in the last section which has a length of 2N 2 The Taylor series is truncated after the first term O dh Therefore this operator is called 2nd order FD operator which refers to the truncation error of the Taylor series and not to the order of the approximated derivative To understand equation 2 9 2 THEORETICAL BACKGROUND 10 better we estimate a 4th order FD operator This operator has the length 2N 4 or N 2 The sums Eq 2 9 lead to of 1 61 30
33. 24 _ 1 Ov OV Ot u 2 OX Ox This expression is called Stress Velocity formulation For simple cases 2 3 and 2 4 can be solved analytically More complex problems require numerical solutions One possible approach for a numerical solution is described in the next section CHAPTER 2 THEORETICAL BACKGROUND 8 2 2 Solution of the elastic wave equation by finite differences 2 2 1 Discretization of the wave equation For the numerical solution of the elastic equations of motion Eqs 2 3 to be discretized in time and space on a grid The particle displacement u the stresses the Lam parameters and are calculated and defined at discrete Cartesian coordinates x i dh j dh and discrete times t n dt dh denotes the spatial distance between two adjacent grid points and dt the difference between two successive time steps Therefore every grid point is located in the interval i N 1 NX j N 1 NY and n N 1 NT where NX NY and NT are the number of discrete spatial grid points and time steps respectively Finally the partial derivatives are replaced by finite difference FD operators Two types of operators can be distinguished forward and backward operators D D The derivative of a function y with respect to a variable x can be approximated by the following operators Dty forward operator dh B ti 1 Fen D y r backward operator To calculate the spatial derivatives of
34. 40 m beneath the free surface The source signature is a butterworth low pass filtered spike with a corner frequency of 15 Hz and order 5 The model has the dimensions 10 km x 3 48 km Using an 8th order spatial FD operator the model can be discretized with 500 x 174 gridpoints in x and z direction with a spatial gridpoint distance of 20 0 m The time is discretized using DT 2 7 ms thus for a recording time of T 6 0 s 2222 time steps are needed 7 1 2 Elastic wave propagation in the complex Marmousi model To generate the measured field data for the Marmousi model copy the files listed in table 7 1 from their origin in the benchmark directory to their destination to define the input and elastic model parameters as well as the acquisition geometry filename origin directory destination directory DENISE_marm_OBC inp benchmark Marmousi classical_FWT PCG DENISE trunk par FWI workflow marmousi inp benchmark Marmousi classical FWT PCG DENISE trunk par receiver OBC dat benchmark Marmousi classical_FWT receiver DENISE trunk par receiver source OBC VSPdat benchmark Marmousi classical_FWT source DENISE trunk par source marmousi II marine vp benchmark Marmousi classical_FWT start Vel DENISE trunk par start marmousi Il marine vs benchmark Marmousi classical_FWT start Vel DENISE trunk par start marmousi II marine rho benchmark Marmousi classical_FWT start Vel DENISE trunk par start marmousi II smooth2 vp benchmark Marmousi c
35. FILE has to be set to 1 Each model parameter requires a taper file which should be located in the par directory and named as taper bin for the Vp model taper u bin for the Vs model and taper rho bin for the density model CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 55 Output of inversion results fp seses Output of inverted models output_of_models_ INV_MODELFILE model model_Test first_model_to_be_saved_ nfstart 2000 increment_between_saved_models_ nf 2000 eec Output of gradients first_gradient_to_be_saved nfstart_jac 2000 increment_between_saved_gradients_ nf_jac 1 Options can produce conflicts in the current version of DENISE It is highly recommended to use the default values Limits for the model parameters Upper and lower limits for model parameters upper_limit_for_vp lambda_ VPUPPERLIM 2600 ower_limit_for_vp lambda_ VPLOWERLIM 500 upper_limit_for_vs mu_ VSUPPERLIM 1600 ower limit for vs mu VSLOWERLIM 1000 upper limit for rho RHOUPPERLIM 5000 ower limit for rho RHOLOWERLIM 0 The six limits for the model parameters specify the minimum and maximum values which may be achieved by the elastic inversion Here known a priori information can be used Depending on the choice of the parameter INVMATI either vp and vs or lambda and mu is limited
36. ICAL BACKGROUND 500 500 N 1000 1000 1500 1500 L 2000 2000 2500 2500 a a 8 3000 8 3000 3500 3500 4000 4000 4500 4500 5000 5000 1000 2000 3000 4000 5000 1000 2000 3000 4000 5000 Distance m Distance m 500 500 1000 1000 1500 1500 2000 2000 2500 2500 8 3000 8 3000 3500 3500 4000 4000 4500 4500 5000 5000 XS 1000 2000 3000 4000 5000 1000 2000 3000 4000 5000 Distance m Distance m 500 500 1000 1000 1500 1500 2000 2000 2500 E 2500 5 a 8 3000 8 3000 3500 3500 4000 4000 4500 4500 5000 5000 1000 2000 3000 4000 5000 Distance m 1000 2000 3000 4000 5000 Distance m 15 Figure 2 3 The influence of grid dispersion in FD modeling Spatial sampling of the wavefield using n 16 top n 4 center and n 2 gridpoints bottom per minimum wavelength Amin CHAPTER 2 THEORETICAL BACKGROUND 16 2 3 2 The Courant Instability Beside the spatial the temporal discretization has to satisfy a sampling criterion to ensure the stability of the FD code If a wave is propagating on a discrete grid then the timestep dt has to be less than the time for the wave to travel between two adjacent grid points with grid spacing dh For an elastic 2D grid this means mathematically dt lt di hV2V max where Vmax is the maximum velocity in the model The factor h depends on the order of the FD operator and can ea
37. PARAMETERS FOR THE MODELING AND INVERSION CODE 59 e EPRECOND 1 approximates the inverse of the Hessian by the absolute value of the forward wavefield Shin et al 2001 1 e EPRECOND 3 approximates the inverse of the Hessian by a zero lag correlation of the absolute value of the forward wavefield with an approximation of the receiver Greens function contribution Plessix and Mulder 1 2004 max __ min 1 xs x t asian TE asi 2 2 where max x denote the minimum maximum receiver and source positions Misfit definition Different objective functions can have a significant impact on the nonlinearity of the inverse problem Changing the misfit function between the modelled data u and field data d does only change the backpropagated residuals in the FWI algorithm LNORM 2 sets the misfit function to the classical L2 norm of the data residuals Eq 3 3 ns nr Er2 a 2 uij dij In this case the misfit is scaled with the energy of the measured seismograms e LNORMSS sets the misfit function to the global correlation norm Choi and Alkhalifah 2012 ns nr gt gt lis Source wavelet inversion As discussed in chapter 4 the estimation of the source wavelet is vital for a successful FWI With STF_INV 1 a source wavelet inversion by a stabilized Wiener deconvolution is activated This requires one additional forward model run per shot to estimate the Gr
38. RUN MULTIPLE SHOTS 1 will start individual model runs from 1 1 to i NSRC with source locations and properties defined at line 1 of the SOURCE FILE The options PLANE WAVE DEPTH TS are obsolete and not supported in the current version of DENISE 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 48 a Source Signals with fc 50 Hz Ricker solid FM dashed sin dotted 15r AN if 0 5 g 2 H 1 0 54 1 i 1 1 bi Voy Va 15 1 1 1 1 1 j 0 10 20 30 40 50 60 Time ms b Amplitude Spectrum ds 0 9 0 8 y 07 5 0 6 o 3 05 E lt 0 4 0 3 0 2 1 014r 1 1 0 1 1 1 0 50 100 150 Frequency Hz Phase Spectrum unwrapped 0 10 te ati ne 20 5 30 8 2 40h 50L 60 270 1 1 50 100 150 Frequency Hz Figure 6 2 Plot of built in source wavelets equations 6 2 6 3 6 4 for a center frequency of f 50 Hz TS 1 fe 0 025 Ricker signal solid Fuchs M ller signal dashed sin signal dotted a Time function b amplitude spectrum c phase spectrum CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 49 Model input Model read_model_from_MFILE 1 READMOD 0 MFILE model test If READMOD 1 the P wave S wave and density model grids are read from external binary files MFILE defines the basic file na
39. SCALEFAC If a corresponding value epst 3 can be found after STEPMAX forward modellings DENISE can fit a parabola through the 3 points L2t i epst i and estimates an optimum step length at the minimum of the parabola If the L2 value L2t 3 after STEPMAX forward models is still smaller than L2t 2 the optimum steplength estimated by parabolic fitting will be not larger than epst 3 Trace killing deeem Trace killing apply trace killing 1 TRKILL 0 TRKILL FILE trace kill trace kill dat To mute noisy or unwanted traces during FWI the parameter TRKILL is introduced If TRKILL is set to 1 all traces defined in the parameter file TRKILL FILE are muted The file should include a mute table where the rows have the same lengths as the number of receivers and the columns reflects the number of sources A 1 ONE indicates a mute of the trace while a 0 ZERO means that this trace should NOT be muted Time damping posos Time windowing and damping files_with_picked_times_ PICKS_FILE picked_times picks_ If time damping of the seismograms is activated in the workflow file by setting TIMEWIN 0 picked times of spe cific seismic phases like first arrivals for each source and receiver must be specified in a seperate file The folder and CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 57 file name can be set with the parameter PICKS_FILE The files must be named like thi
40. TER 3 THE ADJOINT PROBLEM 23 These substitutions yield new equations of motion describing the displacement perturbations and stress perturba tions a function of new source terms Af and 07 du Pa Afi 3 18 Se 00u p 00u m OX Ox initial and boundary conditions The new source terms are Hu Af 6 3 19 P 9g and AT 3 20 Two points are important to notice 1 Eq 3 18 states that every change of a material parameter acts as a source Eq 3 19 and Eq 3 20 but the perturbated wavefield is propagating in the unperturbated medium 2 The new wave equation 3 18 has mathematically the same form as the unperturbated elastic wave equation and hence its solution can be obtained in terms of Green s functions G of the elastic wave equation The solution of the perturbated elastic equations of motion 3 18 in terms of the elastic Green s function G x t x t can be written as f av t x t Af x t p 3 21 T av de 0 Ox Substituting the force and traction terms given by Eqs 3 19 and 3 20 into Eq 3 21 yields after some rearranging u x t pui t x v tap 3 22 Fal ar x t x yt t ejam Introducing isotropy via Aki and Richards 1980 p 23
41. ameter GAMMA Definition of multiparameter inversion The parameters INV VP ITER INV VS ITER INV RHO ITER define from which iteration step the correspond ing parameters Vp Vs and p are updated Setting a parameter to O updates the model during all iterations while setting a parameter larger than ITERMAX no model update is applied Combinations of these parameters allows the implementation of a simultaneous or hierachical inversion workflow Spatial filtering of gradients To suppress short wavelength artefacts below the source and receiver positions the gradients can be smoothed SPATFILTER 1 applies a wavenumber domain damping with a Gaussian function B kx ky g kx ky exp WD_DAMP ky EN to the gradients g kx ky The amount of damping can be controlled by the parameter WD_DAMP SPATFILTER 2 applies a damped least squares technique to the gradients The size of the filter is defined by WD_DAMP in y direction and WD DAMPI in x direction Preconditioning To accelerate the convergence speed of the optimization method and avoid the convergence in a local minimum am plitude loss with depth due to geometrical spreading and reflections in the upper model parts have to be compensated In case of Quasi Newton or Full Newton methods these effects are corrected by the inverse Hessian For conjugate gradient methods different approximations of the inverse Hessian can be used as preconditioning operator CHAPTER 6 DEFINITION OF
42. annel channel Figure 3 1 Definition of data residuals du CHAPTER 3 THE ADJOINT PROBLEM 20 no general rule for an optimum preconditioning operator but two very simple operators are described in more detail chapter for a cross well acquisition geometry and in chapter 7 1 3 for a reflection geometry OE m z P 3 8 Om By replacing m in Eq 3 4 with Eq 3 8 we get OE t P E mo mj M1 54 3 9 The optimum model parameters can be found along the negative gradient direction of the residual energy The starting point m is not a particular point so the update function can be applied to every point in the parameter space my OE Mn 1 Mn EE 3 10 3 3 Calculation of the gradient direction ER To estimate the gradient direction OE 0m the residual energy is rewritten as 1 1 2 ju y fa Su Get 3 11 sources receiver Starting model my Final model mp Figure 3 2 Schematic sketch of the residual energy at one point in space as a function of two model parameters m and ma The blue dot denotes the starting point in the parameter space while the red cross marks a minimum of the objective function CHAPTER 3 THE ADJOINT PROBLEM 21 After derivation with respect to a model parameter m we get sources receiver obs mod _ E 2 _ su 3 12 sources receiver dumed m M fa 5 n sources receiver 3 12 can be rela
43. article velocities SEISMO 2 pressure hydrophones SEISMO 3 curl and div SEISMO 4 everything read receiver positions from file yes 1 READREC 1 REC FILE receiver receiver us dat reference point for receiver coordinate system REFREC 0 0 0 0 if READREC 1 the following three lines are ignored position of first receiver in m XRECIl YREC1 690 0 2300 0 position of last receiver in m XREC2 YREC2 790 0 300 0 distance between two adjacent receivers in gridpoints NGEOPH 6 If SEISMO gt 0 seismograms are saved on hard disk If SEISMO equals 1 x and y component of particle velocity will be written according to parameters specified in Chapter 6 1 If SEISMO 2 pressure sum of the diagonal com ponents of the stress tensor recorded at the receiver locations receivers are hydrophones is written if SEISMO 3 the curl and divergence are saved The curl and divergence of the particle velocities are useful to separate between P and S waves in the snapshots of the wavefield DENISE calculates the divergence and the magnitude of the curl of the particle velocity field according to Dougherty and Stephen 1988 The motivation for this is as follows According to Morse and Feshbach Morse and Feshbach 1953 the energy of P and S wave particle velocities is respectively Ep A 2 div and p rot o _ 6 7
44. ax_relative_error 1 The order ofthe used FD operator is defined by the option FDORDER With the option max_relative_error the user can switch between Taylor max_relative_error 0 and Holberg max_relative_error 1 4 FD coefficients of different accuracy The chosen FD operator and FD coefficients have an influence on the numerical stability and grid dispersion see chapter 2 3 1 Discretization 2 D Grid number of gridpoints in x direction NX 800 number of gridpoints in y direction 200 distance between gridpoints in m DH 2 0e 4 These lines specify the size of the total numerical grid Figure 6 1 NX and NY give the number of grid points in the x and y direction respectively and DH specify the grid spacing in x and y direction The size of the total internal grid in meters in x direction is NX DH and in y direction NY DH To allow for a consistent domain decomposition NX NPROCX and NY NPROCY must be integer values To avoid numerical dispersion the wavefield must be discretized with a certain number of gridpoints per wave length The number of gridpoints per wavelength required depends on the order of the spatial FD operators used in the simulation see section 2 3 1 In the current FD software 2nd 4th 6th 8th 10th and 12th order operators are implemented The criterion to avoid numerical dispersion is defined as Us min DH lt Eh 6 1 with 23 ARA the smallest wavelength propagating t
45. ceivers should certainly be placed outside A convolutional perfectly matched layer CPML boundary condition is used The PML implementation is based on the following papers Komatitsch and Martin 2007 width of the absorbing frame of FW 10 20 grid points should be sufficient For the optimal realization of the PML boundary condition you have to specify the dominant signal frequency FPML occurring during the wave simulation This is usually the center source frequency FC spec ified in the source file DAMPING specifies the attenuation velocity in m s within the PML DAMPING should be approximately the propagation velocity of the dominant wave near the model boundaries In some cases it is usefull to apply periodic boundary conditions see section 2 2 3 IF BOUNDARY 1 no ab sorbing boundaries are installed at the left right sides of the grid Instead wavefield information is copied from left to right and vice versa Therefore a wave which leaves the model at the left side enters the model again at the right side and vice versa Wavetield snapshots Snapshots output of snapshots SNAP gt 0 0 output of particle velocities SNAP 1 output of pressure field SNAP 2 output of curl and divergence energy SNAP 3 output of both particle velocities and energy SNAP 4 first snapshot in sec TSNAP1 3e 6 last snapshot in sec TSNAP2 1 5e 4 increment in TSNAPINC 3e 6 increment x dire
46. cements of the field data are vd 1by 1 Update the real and imaginary parts of the source wavelet with the Newton method OE n 1 n Hn a H 1 far f n Of 4 1 where 28 and 28 are the Hessian matrix and gradient vector for e and f respectively With the objective f function 1 E 5 Sor vA vP v 4 2 the gradients and Hessian can be explictly calculated OE die Mec ae di d rder PS by rcv c OPE 32 re 4 3 PE_ fde OE _ Ot ie d 33 CHAPTER 4 SOURCE WAVELET INVERSION 34 Inserting Eq 4 3 in Eq 4 1 leads to av rev r by dy 8n41 mer T f 2 Xp av rdy r m by To stabilize the inversion a Marquardt Levenberg regularization is required _1 6541 en Ha 1 7 OE fayi f Ha Ad 4 ap where Af are damping factors and I the unity matrix Therefore Eq 4 4 can be written as a by 1 8n41 EG de 4 5 f av rdy r NE by TES PR m d2 x The values of the damping factors can be expressed as fractions of the maximum value in the denominators Ae Ar stf Max 6 2 d 4 6 r This approach is a stabilized Wiener deconvolution Chapter 5 Getting Started In the follo
47. center frequency and covers a broader frequency range You may also use your own time function as the source wavelet for instance the signal of the first arrival recorded by a geophone at near offsets Specify QUELLART 3 and save the samples of your source wavelet in ASCII format in SIGNAL FILE SIGNAL FILE should contain one sample per line It should thus look like The time interval between the samples must equal the time step interval DT of the FD simulation see above Therefore it might be necessary to resample interpolate a given source time function with a smaller sample rate You may use the matlab script mfiles resamp m to resample your external source signal to the required sampling interval The following source types are availabe explosive sources that excite compressional waves only QUELLTYP 1 and point forces in the x and y direction QUELLTYP 2 3 The force sources excite both P and S waves The explosive source is located at the same position as the diagonal elements of the stress tensor i e at i j Figure 2 1 The forces are located at the same position as the corresponding components of particle velocity Figure 2 1 If x y denotes the position at which the source location is defined in source dat then the actual force in x direction is located at x DX 2 y and the actual force in y direction is located at x y DY 2 With QUELLTYP 4 a custom directive force can be defined by a force angle between y and x The ang
48. compared with the true model in Fig 7 6 The results contain lot of small details fine layers which are completely absent in the initial model are resolved The thrust faults and reef structures in the deeper part of the model are imaged also very well AII hydrocarbon reservoirs can be identifed and even structures at the scale of the FD grid are resolved The shear wave velocity model could also be reconstructed even though only streamer data and therefore only P wave information is used This can be explained by the replacement of the soft seabed in the original Marmousi II model by the hard seabed used in this study which leads to significant P S conversions and S P conversions at the seafloor Therefore the hard seabed produces a significant footprint in the shear wave velocity model Even the density a parameter which can be hardly estimated from seismic data could be recovered from the seismic wavefield Keep in mind though that the density image is based not only on the density information but contains also structural and information due to the ambiguity investigated by the CTS test problem K hn et al 2012 The quality of the inversion results is also evident in the seismic sections of shot 50 vertical component plotted in Fig 7 7 Notice the direct wave the reflection from the ocean bottom a few multiples and the dominant reflection from the interface between low velocity and high velocity sediments but the lack of ot
49. ction IDX 1 increment y direction 1 data format SNAP FORMAT ASCII 2 BINARY 3 3 basic filename SNAP FILE snap waveform forward If SNAP gt 0 wavefield information particle velocities pressure or curl and divergence of particle velocities for the entire model is saved on the hard disk assure that enough free space is on disk Each PE is writing his sub volume to disk The filenames the basic SNAP_FILE plus an extension that indicates the number in the logical processor array see Figure 6 1 i e the with number PEno writes his wavefield SNAP FILE PEno The first snapshot is written at TSNAPI seconds of seismic wave traveltime to the output files the second CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 51 at TSNAP1 TSNAPINC seconds etc The last snapshots contains wavefield at TSNAP2 seconds Note that the file sizes increase during the simulation The snapshot files might become quite LARGE It may therefore be necessary to reduce the amount of snapshot data by increasing IDX IDY and or TSNAPINC In order to merge the separate snapshot of each PE after the comletion of the wave modeling you can use the program snapmerge see Chapter 5 2 section src The bash command line to merge the snapshot files can look like this bin snapmerge DENISE inp Receivers Receiver output_of_seismograms_ SEISMO 1 SEISMO 0 seismograms SEISMO 1 p
50. d by a red cross the starting point by a blue point The number of iterations is 30 CHAPTER 3 THE ADJOINT PROBLEM 32 3 8 elastic FWT algorithm In summary the FWT algorithm consists of the following steps 1 Define a starting model m in the parameter space This model should represent the long wavelength part of the underground very well because the FWT code is only capable to reconstruct structures at or below the dominant seismic wavelength due to its slow convergence speed the nonlinearity of the problem and the inherent use of the Born approximation to calculate the gradient direction 2 At iteration step do a For each shot solve the forward problem stated in Eq 3 16 for the actual model m to generate a synthetic dataset uM d and the wavefield u x t b Calculate the residual seismograms du 1172994 u bs for the x and y components of the seismic data c Generate the wavefield W x t by backpropagating the residuals from the receiver postions d Calculate the gradients 6m of each material parameter according to Eqs 3 36 e To increase the convergence speed an appropriate preconditioning operator P is applied to the gradient dm P m 3 57 Examples of simple preconditioning operator are given in chapter for a cross well acquisition geometry and in chapter 7 1 3 for a reflection geometry f For a further increase of the convergence speed calculate the conjugate gradient direction f
51. d maximum velocity values occure in the viscoelastic model The frequency dependence of attenuation 1 and phase velocity a function of frequency may be calculated using the Matlab functions in the directory mfiles Free surface FE Surfac surface_ yes 1 FREE SURF 1 CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 50 A plane stress free surface is applied at the top of the global grid if FREE_SURF 0 using the imaging method proposed by Levander 1988 Note that the free surface is always located at y 0 Boundary conditions PML Boundary quadratic damping applied width_of_absorbing_frame_ in_gridpoints _ No lt 0 _ FW 10 Damping velocity in CPML in m s DAMPING 1277 0 Frequency within the PML Hz FPML 100 0 npower 4 0 k max PML 1 0 apply periodic boundary condition at edges BOUNDARY no20 left and right 1 0 The boundary conditions are applied on each side face and the bottom face of the model grid If FREE SURF 0 the boundary conditions are also applied at the top face of the model grid Note that the absorbing frames are always located INSIDE the model space i e parts of the model structure are covered by the absorbing frame in which no physically meaningful wavefield propagates You should therefore consider the frame width when you design the model structure and the acquisition geometry shot and re
52. dlimited spike wavelet A spike bandlimited by a low pass or band pass butterworth filter to avoid grid dispersion If FC SPIKE 1 lt 0 0 a low pass filtered spike with upper corner frequency SPIKE 2 and order OR DER SPIKE is calculated If FC SPIKE 1 gt 0 0 a band pass filtered spike with lower corner frequency SPIKE 1 and upper corner frequency SPIKE 2 with order ORDER SPIKE is calculated Klauder wavelet A Klauder wavelet represents the autocorrelation of a linearly swept frequency modulated sinusoidal signal used in Vibroseis TkT TS klau t real sin ES with 7 t 1 5 FC SPIKE 1 ta 6 6 CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 47 with SPIKE 2 FC SPIKE 1 TS rate of change of frequency with time fo SPIKE 2 FC SPIKE 1 2 midfrequency of bandwidth i v l In these equations t denotes time and fe is the center frequency 4 is a time delay which can be defined for each source position in SOURCE FILE Note that the symmetric zero phase Ricker signal is always delayed by 1 0 fe which means that after one period the maximum amplitude is excited at the source location These 5 source wavelets and the corresponding amplitude spectra for a center frequency of fe 50 Hz and a delay of tg 0 are plotted in Figure 6 2 Note the delay of the Ricker signal described above The Fuchs Miiller wavelet has a slightly higher
53. e parabolic fit Starting from 0 0 the test step lengths and ps are calculated to satisfy the following criteria L25 mtest2 Mn lt 121 Mtest Mn 3 50 23 Mtest3 Mn gt L22 Mtest2 Mn This approach leads to smoother decrease of the objective function but also increases the number of required forward models 3 7 Nonlinear Conjugate Gradient Method To increase the convergence speed in narrow valleys it would be better to update the model at iteration step n not exactly along the gradient direction m but along the conjugate direction c c Om B10C 1 3 51 The first iteration step n 1 consists of a model update along the steepest descent direction ma My 10011 3 52 For all subsequent iteration steps gt 1 the model is updated along the conjugate direction My41 Mn 3 53 where c The weighting factor can be calculated in different ways Nocedal and Wright 2006 CHAPTER 3 THE ADJOINT PROBLEM 30 Residual Energy 250 200 150 y gt 100 x gt Figure 3 6 Results of the convergence test for the Rosenbrock function using the pure gradient method The minimum is marked by a red point the starting position by a blue point The number of iterations is 200 The optimum step length is calculated at each iteration by the parabola fitting algorithm Note the criss cross pattern of the updates i
54. ed 2 order finite difference operator sin 25429 2 2 FE 2 20 du 2 1 dx m GY x 0 dux 3Ax u x dx x 0 x 0 Using the Nyquist Shannon sampling theorem it should be sufficient to sample the wavefield with 2 2 In table 2 1 the numerical solutions of 2 20 and the analytical solution 2 19 are compared for different sample intervals A n where n is the number of gridpoints per wavelength For the case n 2 which corresponds to the Nyquist Shannon theorem the numerical solution is dus x 0 4 0 which is not equal with the analytical solution 27 A refinement of the spatial sampling of the wavefield results in an improvement of the finite difference solution For n 16 the numerical solution is accurate to the second decimal place The effect of a sparsly sampled pressure field is illustrated in figure 2 3 for a homogeneous block model with stress free surfaces The dimensions of the FD grid are fixed and the central frequency of the source signal is increased systematically When using a spatial sampling of 16 grid points per minimum wavelength figure 2 3 top the wavefronts are sharply defined For n 4 grid points a slight numerical dispersion of the wave occurs figure 2 3 center This effect is obvious when using the Nyquist criterion n 2 figure 2 3 bottom Since the numerical calculated wavefield seem to be dispersive this numer
55. eens function solution for the actual model The parameter STF defines how many traces in the vicinity of the shot point are used In case of dispersive wavefields it is recommended to limit the source wavelet inversion only to the near offset traces to avoid the interpretation of model parameter changes as source wavelet Traces with maximum offsets OFFSETC STF are used for the wavelet inversion EPS_STF denotes the regularization parameter defined in Eq 4 6 So far the source wavelet will be only estimated from the vertical component data In case of STF_INV 2 the first arrival of the synthetic data will be automatically picked by an STA LTA picker and the amplitudes of the field data before the first arrival and after the first arrival plus a time window defined by TW STF are exponentially damped The amount of damping is controlled GAMMA STF TW STF and GAMMA STF currently have to be set in time window stf c Normalization Normalize seismic data and source wavelet during source wavelet inversion This can be required if the amplitudes of source wavelets for each shot show strong variations The normalization leads to an equalization of the shot contributions to the gradient NORMALIZE 0 No normalization of seismic data and source wavelet during source wavelet inversion NORMALIZE 1 Trace normalization of synthetic and field data to the maximum amplitude of each trace Normalize the estimated source wavelet to its maximum amplitude
56. el of the force must be specified in the SOURCE FILE after AMP This force is not aligned along the main directions The parameter ANGLE is without any use The locations of multiple sources must be defined in an external ASCII file SOURCE FILE that has the following format NSRC ZSRC TD FC SOURCE AZIMUTH SOURCE TYPE NSRC lines In the following lines you can define certain parameters for each source point the first line must be the overall number of sources NSRC XSRC is the x coordinate of a source point in meter YSRC is the y coordinate of a source point in meter ZSRC is the z coordinate should always be set to 0 0 because DENISE is a 2D code TD is the excitation time time delay for the source point in seconds FC is the center frequency of the source signal in Hz and AMP is the maximum amplitude of the source signal Optional parameter The SOURCE AZIMUTH if SOURCE TYPE is 4 The SOURCE AZIMUTH is the angle between the y and x direction and with SOURCE TYPE if SOURCE is set here the value of SOURCE TYPE in the input file is ignored The SOURCE FILE sources source dat that defines an explosive source at 2 2592 0 m and y 2106 0 m with a center frequency of 5 Hz no time delay is 2592 0 0 0 2106 0 0 0 5 0 140 If the option RUN_MULTIPLE_SHOTS 0 in the parameter file all shot points defined in the SOURCE FILE are excitated simultaneously in one simulation Setting
57. ep in mind that the file system WORK will not be automatically backuped so do a manual backup from time to time After compiling the code see section 5 3 you can define and start a batch job with a shell script like this bin ksh 5 1 elapstim req 48 00 00 walltime PBS 1 cputim 768 00 00 akkumulated CPU time per nod PBS 1 memsz job 100gb RAM requirement PBS p 1 number of nodes PRS T intmpi job type intmpi for Intel MPI PBS 1 cpunum job 16 number of cores node PBS N DENISE name of the batch job PBS DENISE out file for standard output PES Y standard and error output PBS q clmedium batch class opt intel composer xe 2013 spl bin compilervars sh intel64 opt intel impi 4 1 1 036 intel64 bin mpivars sh opt modules Modules 3 2 6 init ksh cd SWORK DENISE_PSV par mpirun SNOSITI PIOPIS enp 15 bin denise DENISE marm OBC inp FWI workflow marmousi inp The individual parameters and possible batch job classes are described in more detail on the homepage of the RZ Kiel https www rz uni kiel de hpc nec cluster html The most important parameters are elapstim req which defines how long the job will actually run cputim_job the accumulated CPU time per node memsz job the required memory per node b the number of nodes mpich the job type in case of Intel MPI you have to choose mpich cpunum job number of CPUs per node
58. epth top Hessian approximation by Plessix and Mulder 2004 with PCG centre and LBFGS bottom CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL 75 7 1 4 Marmousi II Benchmarks Due to the massive amount of required forward models the Marmousi II inversion is a perfect benchmark problem A few representative benchmark results are shown in table 7 3 Due to continous bugfixing and improvements in the source code the computation times are not comparable PGI 11 4 MPICH2 1 8 8 20 3 Intel 14 0 0 Intel MPI 4 1 Table 7 3 Marmousi II benchmark results Bibliography K Aki and P Richards Quantitative seismology W H Freeman and Company 1980 S al Hagrey D K hn and W Rabbel Geophysical assessments of renewable gas energy compressed in geologic pore storage reservoirs SpringerPlus 3 1 1 16 2014 J Blanch J Robertsson and W Symes Modeling of a constant Q Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique Geophysics 60 1 176 184 1995 T Bohlen Interpretation of Measured Seismograms by Means of Viscoelastic Finite Difference Modelling PhD thesis Kiel University 1998 T Bohlen Parallel 3 D viscoelastic finite difference seismic modelling Computers amp Geosciences 28 8 887 899 2002 T Bohlen and E Saenger Accuracy of heterogeneous staggered grid finite difference modeling of Rayleigh waves Geophysics 71 4 T109 T115 2006 A Brenders and
59. equency domain Choi et al 2008a b Brossier 2009 In order to extract information about the structure and composition of the crust from seismic observations it is necessary to be able to predict how seismic wavefields are affected by complex structures Since exact analytical solutions to the wave equations do not exist for most subsurface configurations the solutions can be obtained only by numerical methods For iterative calculations of synthetic seismograms with limited computer resources fast and accurate modeling methods are needed The FD modeling inversion program DENISE subwavelength DEtail resolving Nonlinear Iterative SEismic inversion is based on the FD approach described by Virieux 1986 and Levander 1988 The present program DENISE has the following extensions is efficently parallelized using domain decomposition with MPI Bohlen 2002 considers viscoelastic wave propagation effects like attenuation and dispersion Robertsson et al 1994 Blanch et al 1995 Bohlen 2002 employs higher order FD operators CHAPTER 1 INTRODUCTION 5 applies Convolutional Perfectly Matched Layer boundary conditions at the edges of the numerical mesh Ko matitsch and Martin 2007 In the following sections we give an extensive description of the theoretical background the different input pa rameters and show a few benchmark modeling and inversion applications 1 1 Citation If you use this code for your own researc
60. erately sized problems could be inverted with frequency domain approaches A spectacular result to prove the application of acoustic FWT on laboratory scale was presented by Pratt 1999 for ultrasonic tomography measurements on a simple block model In a numerical blind test Brenders and Pratt 2007 achieved a very good agreement between their inversion result and the unkown true P wave velocity model The paral lelization and performance optimizations of the frequency domain approach see e g Sourbier et al 2009a Sourbier et al 2009b lead to a wide range of acoustic FWT applications for problems on different scales from the global scale crustal scale over engineering and near surface scale down to laboratory scale Pratt 2004 Beside the application to geophysical problems the acoustic FWT is also used to improve the resolution in medical cancer diagnostics Pratt et al 2007 However all these examples are restricted to the inversion of the acoustic material parameters P wave velocity density and additionally the viscoacoustic damping Q for the P waves Even today the independent 2D FWT of all three isotropic elastic material parameters is still a challenge Most elastic ap proaches invert for P wave velocity only and use empirical relationships to deduce the distribution of S wave velocity and density Shipp and Singh 2002 Sheen et al 2006 Recently some authors also investigated the independent multiparameter FWT in the fr
61. face only horizontal particle displacements should be used because vertical derivatives over the free surface lead to instabilities Levander 1988 The vertical derivative of the y displacement u can be replaced by using the boundary condition at the free surface 2 Ali 0 Uyy u 7 Therefore the stress o can be written as c p A 2u Uxx 2 15 2 16 2 Absorbing Boundary Conditions Due to limited computational resources the FD grid has to be as small as possible To model problems with an infinite extension in different directions e g a full or half space problem an artificial absorbing boundary condition has to be applied A very effective way to damp the waves near the boundaries are Perfectly Matched Layers PMLs This can be achieved by a coordinate stretch of the wave equations in the frequency domain Komatitsch and Martin 2007 The coordinate stretch creates exponentially decaying plane wave solutions in the absorbing boundary frame The PML s are only reflectionless if the exact wave equation is solved As soon as the problem is discretized for example using finite differences you are solving an approximate wave equation and the analytical perfection of the PML is no longer valid To overcome this shortcoming the wavefield is damped by the damping function log a L Vpmi 2 17 where denotes the typical P wave veloc
62. fference elastic waveform tomography in the time domain In SEG Technical Program Expanded Abstracts pages 2302 2306 2009 M Dougherty and R Stephen Seismic energy partitioning and scattering in laterally heterogeneous ocean crust Pure Appl Geophys 128 1 2 195 239 1988 R Fletcher and C M Reeves Function minimization by conjugate gradients The Computer Journal 7 2 149 154 1964 76 BIBLIOGRAPHY 77 Gauthier J Virieux and A Tarantola Two dimensional nonlinear inversion of seismic waveforms numerical results Geophysics 51 7 1387 1403 1986 M Hestenes and E Stiefel Methods of conjugate gradients for solving linear systems Journal of Reasearch of the National Bureau of Standards 49 6 409 436 1952 O Holberg Computational aspects of the choice of operator and sampling interval for numerical differentiation in lage scale simulation of wave phenomena Geophysical Prospecting 35 629 655 1987 C Jastram Seismische Modellierung mit Finiten Differenzen h herer Ordnung auf einem Gitter mit vertikal variieren dem Gitterabstand PhD thesis Universit t Hamburg 1992 D K hn Time Domain 2D Elastic Full Waveform Tomography PhD thesis Kiel University 2011 available at http nbn resolving de urn nbn de gbv 8 diss 67866 D Kohn D De Nil A Kurzmann A Przebindowska and T Bohlen On the influence of model parametrization in elastic full waveform tomography Geophysical Journal Internatio
63. fficients of this matrix equation are formally defined by lb 3 47 In the FWT code the solution vector x is calculated by Gaussian elimination In the following the step length at the extremum of the parabola defines the optimum step length opi denoted as green square Fig 3 5 This optimum step length is Hopt 3 48 The application of the variable step length calculation to the Rosenbrock test problem is shown in Fig 3 6 The number of required iteration steps to reach the minimum is reduced by a factor 80 when compared with the constant step length gradient method The only problem remaining is the slow convergence speed in the small valley of the Rosenbrock function due to the fact that the update occurs along the gradient direction of the objective function resulting in a criss cross pattern This behaviour can be avoided by applying a nonlinear conjugate gradient method chapter 3 7 In case of the FWT algorithm the three test step lengths for the individual material parameter classes are calculated by scaling the maximum of the gradient to the maximum of the actual models max An pu Pmax 6A _ Pix 3 49 Mp P Srho max dpn CHAPTER 3 THE ADJOINT PROBLEM Residual Energy E xo Residual Energy E 250 200 150 100 150 250 200 150 100 150 28 Figure 3 4 Results of the convergence test for the Rosenb
64. flow inp gt DENISE out After the output of geometry and model parameters the code starts the time stepping and displaying timing infor mation MYID 0 Starting simulation forward model for shot 1 of 1 xxxxxxx x xx Number of samples nts in source file 3462 Number of samples nts in source file 3462 Message from function wavelet written by PE 0 1 source positions located in subdomain of PE 0 have been assigned with a source signal PE 0 outputs source time function in SU format to start source signal 0 su shotl Computing timestep 1000 of 3462 Message from update v printed by PE 0 Updating particle velocities finished real time 0 00 s CHAPTER 5 GETTING STARTED 41 particle velocity exchange between PEs finished real time 0 00 Message from update s printed by PE 0 Updating stress components finished real time 0 00 s stress exchange between PEs finished real time 0 00 s total real time for timestep 1000 0 01 s Computing timestep 2000 of 3462 Message from update v printed by PE 0 Updating particle velocities 2 finished real times 0 00 particle velocity exchange between PEs finished real time 0 00 s x Message from update s printed by PE 0 Updating stress components finished real time 0 00 s stress exchange between PEs finished real time 0 00 s total real time for timestep 20
65. g instructions you must first install the LAM MPI software On SUSE LINUX systems the package can be installed with yast The software can also be downloaded from http www lam mpi org A good documentation of LAM MPI is available at this site Before you can run the DENISE software on more than one node you must boot LAM This requires that you can connect to all of your nodes with ssh without a password To enable ssh connection without a password you must generate authentication keys on each node using the command ssh keygen t rsa These are saved in the file ssh id_rsa pub Copy all authentication keys into one file named authorized_keys and copy this file on all nodes into the directory ssh Before you can boot LAM you must specify the IP addresses of the different processing elements in an ASCII file An example file is par lamhosts Each line must contain one IP address The first IP number corresponds to node 0 the second line to node 1 and so forth Note that in older LAM MPI implementations the mpirun command to run the FD programs must always be executed on node 0 i e you must log on node 0 first line in the file par lamhosts and run the software on this node In the last stable release of LAM MPI the node 0 just has to be listed in the lamhosts file the order does not matter To boot lam just do lamboot v par lamhosts To run the software in serial on a single PC just type lamboot without specifying a lamhosts file You can still
66. h please cite at least one article written by the developers of the package for instance Kohn 2011 Kohn et al 2012 or XX add more references here and or other articles from http www geophysik uni kiel de dkoehn publications htm The corresponding BibT X entries may be found in file doc USER_MANUAL thesis bib 1 2 Support The development of the code was suppported by the Christian Albrechts Universit t Kiel TU Bergakademie Freiberg Deutsche Forschungsgemeinschaft DFG Bundesministerium fiir Bildung und Forschung BMBF Bundesminis terium fiir Umwelt Naturschutz Bau und Reaktorsicherheit BMU the Wave Inversion Technology WIT Con sortium and the Verbundnetz Gas AG VNG The code was tested and optimized at the computing centres of Kiel University TU Bergakademie Freiberg TU Chemnitz TU Dresden the Karlsruhe Institute of Technology KIT and the Hochleistungsrechenzentrum Nord HLRN 1 2 Acknowledgments and contact We thank for constructive discussions and further code improvements Anna Przebindowska Karlsruhe Institute of Technology Olaf Hellwig TU Bergakademie Freiberg Dennis Wilken and Wolfgang Rabbel Christian Albrechts Universit t Kiel Please e mail your feedback questions comments and suggestions to Daniel K hn dkoehn AT geophysik uni kiel de Chapter 2 Theoretical Background 21 Equations of motion for an elastic medium The propagation of waves in a general ela
67. her events beyond the first arrivals in the seismic section of the starting model The initial data residuals show no direct wave anymore so it is fitted perfectly The residuals only contain events from small scale model features The seismic sections of the FWT result and the true model are nearly identical so the final data residuals are very small reflection events are fitted perfectly The normalized misfit function for the different frequency bands decreases very fast The influence of different Hessian approximations on the inversion results is shown for the V model in Fig 7 8 PCG optimization with a simple linear scaling of the gradient with depth top shows a bad model resolution in the deeper parts of the model The Hessian approximation by Plessix and Mulder 2004 with PCG centre improves the depth resolution The combination with the LBFGS method bottom adds an additional depth resolution and overall improvement of the model sharpness ale E 7 E i gt x B LX 8812 2 8 gt 2 le lei 212 1 un m Y E gt J gt es n lt E E E RR IB EISlElEIZTZ ZIE EIE 5 E 5 8 00171100120 6 0 1 2 00 00 0 0 0 0 160 160 3 2 0 20000 1 1 0 0 100 105 001 1 00 50 6 01142100 00 0 160 16013 2 0 2000 le 1 O 100 05 0 01 1 00 100 6 0 16 00 o 160 16013 2 0 20000 le 1
68. hrough the model vs min denotes the minimum shear wave velocity in the model and f 1 T S is the center frequency of the source wavelet The program assumes that the maximum frequency of the source signal is approximately two times the center frequency The center frequency is approximately one over the duration time TS The value of n for different FD operators is tabulated in table 2 2 The criterion 6 1 is checked by the FD software If the criterion is violated a warning message will be displayed in the DENISE output section CHECK FOR GRID DISPERSION Please note that the FD code will NOT terminate due to grid dispersion only a warning is given in the output file Time stepping Time Stepping time_of_wave_propagation_ in_sec _ TIME 1 8e 4 timestep_ in_seconds _ DT 5 2e 8 CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 46 The propagation time of seismic waves in the entire model is TIME time stepping interval DT has to fulfill the stability criterion 2 22 in section 2 3 2 The program checks these criteria for the entire model outputs a warning message if these are violated stops the program and will output the time step interval for a stable model run Sources Soure Shape_of_source signal ricker 1 fumue 2 SOURCE FILE 3 SIN 3 4 Gauss deriv 5
69. ical artefact is called grid dispersion To avoid the occurence of grid dispersion the following criteria for the spatial grid spacing dh has to be satisfied dh lt Amin V min 2 21 n nf Here Amin denotes the minimum wavelength Vmin the minimum velocity in the model and fmax is the maximum frequency of the source signal Depending on the accuracy of the used FD operator the parameter n is different In table 2 2 n is listed for different FD operator lengths and types Taylor and Holberg operators The Holberg coefficients are calculated for a minimum dispersion error of 0 196 at 3fmax For short operators n should be choosen relatively large so the spatial grid spacing is small while for longer FD operators n is smaller and the grid spacing can be larger CHAPTER 2 THEORETICAL BACKGROUND 14 n Ax fm Gelee analytical 2 6 283 2 2 2 4 0 4 2 4 5 657 8 A 8 6 123 16 2 16 6 2429 32 A 32 6 2731 Table 2 1 Comparison of the analytical solution Eq 2 19 with the numerical solution Eq 2 20 for different grid spacings Ax A n FDORDER n Taylor n Holberg 2nd 12 12 4th 8 8 32 6th 6 4 77 8th 5 3 69 10th 5 3 19 12th 4 2 91 Table 2 2 The number of grid points per minimum wavelength n for different orders 2nd 12th and types Taylor and Holberg of FD operators For the Holberg coefficients n is calculated for a minimum dispersion error of 0 1 at 3fmax CHAPTER 2 THEORET
70. if TIMELAPSE 1 DATA DIR should be the directory containing the data differences between time 0 and tl Seismograms of synthetic data at t0 DATA DIR TO su CAES spike time 0 DENISE CAES If TIMELAPSE 1 the spatial FWI is replaced by a double difference time lapse FWI Denli and Huang 2009 al Hagrey et al 2014 In this case DATA_DIR defines the data differences between the baseline data at time tO and the time lapse data at tl For existing SU files the data differences can be calculated with the shell script time_lapse_data_diff sh in the par directory The location of the baseline data can be defined with DATA_DIR_TO Elastic Reverse Time Migration f Elastic Reverse Time Migration apply_reverse_time_migration_ yes 1 _ RTM 0 If RTM 1 and INVMAT 0 an elastic Reverse Time Migration RTM is applied for the field data defined in the directory DATA_DIR If time lapse mode is activated TIMELAPSE 1 the time lapse data will be migrated The workflow file section 6 2 should only contain one stage The resulting migrated seismic sections are written to the directory JACOBIAN Currently Reverse Time Migration is only defined for the L2 Norm L2NORM 2 CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 58 6 2 Workflow file with variable inversion parameters FWI_workflow inp Complex FWI workflows can be designed with the input file shown in table 7 2 Each line represents a FWI stage with a spec
71. ific combination of different inversion parameters defined in the columns Abort criterion Beside the parameter ITERMAX a second abort criterion is implemented in DENISE which is using the relative misfit change within the last two iterations The relative misfit of the current iteration step and the misfit of the second to last iteration step is calculated with regard to the misfit of the second to last iteration step If this relative change is smaller than PRO the inversion aborts or proceeds to the next inversion stage Frequency filtering To tame the nonlinearity of the inversion problem Butterworth frequency filters can be applied to the source wavelet and field data e TIME_FILT 1 applies a lowpass frequency filter with an upper corner frequency FC high TIME FILT 2 applies a bandpass frequency filter with a lower corner frequency FC low and upper corner frequency FC high The order of the Butterworth filter is defined by the parameter ORDER Time damping Multiple or complex reflections can significantly increase the nonlinearity of the inverse problem Different time damping strategies are implemented in DENISE to TIMEWIN 1 reads traveltime picks of the first arrival from the PICKS FILES defined in the parameter file section 6 A constant time delay TWIN can be applied to each pick TIMEWIN 2 applies a time damping from a constant time TWIN for all receivers and shots The amount of damping can be defined by the par
72. ion control system subversion SVN This directory contains internal informa tion on recent software updates etc You can ignore this directory a clean release should not have any svn directories at all Each subdirectory described below also contains a subdirectory SVN which can be ignored as well bin This directory contains all executable programs generally DENISE and snapmerge These executables are generated using the command make lt program gt see below jobs This directory contains Batch scripts to submit SAVA modelling inversion runs on HPCs with PBS batch system libcseife This directory contains external contributions to DENISE for the implementation of a Butterworth frequency filter mfiles Here some Matlab routines m files are stored These Matlab programs can be used to find optimal relaxation frequen cies to approximate a constant qapprox m or to plot Q as a function of frequency for certain relaxation frequencies CHAPTER 5 GETTING STARTED 38 and value of tau qplot m For further details on the theory behind these algorithms we refer to the Phd thesis of T Bohlen Bohlen 1998 and to the paper in which the so called tau method is described Blanch et al 1995 The parameter tau is used in this software to define the level of attenuation par Parameter files for DENISE modelling and inversion src This directory contains the complete source codes 5 3 Compilation of DENISE Before compiling DENISE you have to c
73. is surrounded by a very narrow valley We start the search for the minimum at 0 5 0 5 An obvious first choice would be a constant step length Fig 3 4 top shows the convergence after 16000 iteration steps of the steepest descent method when choosing a step length 2e 3 Note the large model update during the first iteration step when the gradient of the Rosenbrock function is large After reaching the narrow valley the gradient is much smaller and as a result the model updates are also decreasing This leads to a very slow convergence speed Especially near the minimum the model updates become very small When choosing a larger step length 2e 3 Fig 3 4 bottom the model update is larger even when the gradient is small but the code fails to converge at all Instead it is trapped in a narrow part of the valley To solve this problem a variable step length is introduced For three test step lengths 444 and three test models are calculated Mtest1 Mn jdm Mtest2 Mn 20m 3 44 Meest3 Mn ua m and the corresponding L2 norms 1 21 125 and 1 25 are estimated Fig 3 5 The true misfit function yellow line can be approximated by fitting a parabola through the three points L2 au bui c 3 45 where i 1 2 3 anda b are the unkown coefficients This system of equations can be written as matrix equation uU 1 L2 ps 1 c L23 or Ax b 3 46 The unknown coe
74. ity of the medium in the absorbing boundary frame a 1 x 1074 and L is the thickness of the absorbing boundary layer A comparison between the exponential damping and the PML boundary is shown in Fig 2 2 The PMLs are damping the seismic waves by a factor 5 10 more effective than the absorbing boundary frame CHAPTER 2 THEORETICAL BACKGROUND Time 1206 6 ms Time 1206 6 ms 12 Figure 2 2 Comparison between exponential damping left column and PML right column absorbing boundary conditions for a homogeneous full space model CHAPTER 2 THEORETICAL BACKGROUND 13 2 3 Numerical Artefacts and Instabilities To avoid numerical artefacts and instabilities during a FD modelling run spatial and temporal sampling conditions for the wavefield have to be satisfied These will be discussed in the following two sections 2 3 1 Grid Dispersion The first question when building a FD model is What is the maximum spatial grid point distance dh for a correct sampling of the wavefield To answer this question we take a look at this simple example The particle displacement in x direction is defined by a sine function x sin 272 2 18 u sin r where denotes the wavelength When calculating the derivation of this function analytically at x 0 and setting 1 m we get 2m o X cos 27 A In the next step the derivation is approximated numerically by stagger
75. kefile for DENISE edit here source code for model generation MODEL_EL half_space c bin Compiler LAM CC hcc CRAY CC cc ON Linux cluster running LAM options for DENISE CU hcc LFLAGS 1m lmpi lcseif CFLAGS 03 SFLAGS L libcseife IFLAGS I libcseife On CRAY T3E CHAPTER 5 GETTING STARTED On NEC cluster with Intel MPI 39 LFLAGS 1m lcseif fir eblstudeat CFLAGS O3 xAVX ipo fno fnalias restrict SFLAGS L libcseife L sfs fs4 work sh2 sungw331 fftw 3 3 4 1lib IFLAGS I libcseife I sfs fs4 work sh2 sungw331 fftw 3 3 4 include On MARWIN CC mpiicc LFLAGS 1m lcseif lftftw3 Letdct CFLAGS 03 SFLAGS L libcseife IFLAGS I libcseife On Desktop computer with LinuxMint 17 OpenMPI and gcc 4 8 2 CC mpicc LFLAGS m lcseif lfftw3 lstdos4 CFLAGS O3 fno stack protector SFLAGS L libcseife IFLAGS I libcseife On HLRN system CC mpcc FLFLAGS 1m ALTIX CFLAGS mp 03 LFLAGS 1mpi i static after this line no further editing should be necessary To compile the program DENISE you must change to the src directory and execute bash 2 05b DENISE src gt make denise The following or a similar output sh
76. l Free Surface The interface between the elastic medium and air at the surface is very important when trying to model surface waves or multiple reflections in a marine environment Since all stresses in the normal direction at this interface vanish Oxy 0 0 2 13 this boundary condition is called stress free surface Two types of implementations are common In the im plicit defintion of the free surface a small layer with the acoustic parameters of air Vp 300 m s Vs 0 0 m s 1 25 kg m is placed on top of the model One advantage of the implicit definition of the free surface is the easy implementation of topography on the FD grid however to get accurate results for surface waves or multiples this approach requires a fine spatial sampling of the FD grid near the free surface An explicit free surface can be implemented by using the mirroring technique by Levander which leads to stable and accurate CHAPTER 2 THEORETICAL BACKGROUND 11 solutions for plain interfaces Levander 1988 Robertsson et al 1995 If the planar free surface is located at grid point 7 h the stress at this point is set to zero and the stresses below the free surface are mirrored with an inverse sign Oyy h i 0 Oyy h 1 1 1 1 1 1 1 1 Oxy h 51 5 oxy h 5 1 5 2 14 fll 3 1 Oxy 9 az 21 5 When updating the stress component ox A 24 Uxx Auyy at the free sur
77. lassical_FWT start Vel DENISE trunk par start Il smooth2 vs benchmark Marmousi classical_FWT start Vel DENISE trunk par start marmousi II smooth2 rho benchmark Marmousi classical_FWT start Vel DENISE trunk par start Table 7 1 Origin and destination directories of the files required for the modeling and inversion of the Marmousi model The optimum FD modeling section 7 1 1 and inversion parameters are already defined in the input file DENISE_marm_spike inp However before running the forward code check if the following parameters are set correctly READMOD 1 MFILE start marmousi_II_marine 10 Go to par directory and run the forward code 15 cores of the cluster by typing mpirun np 15 bin denise DENISE_marm_OBC inp FWI_workflow_marmousi inp or submit a job file on the NEC cluster The resulting seismograms for the v and v component are written to DENISE trunk par su You can check the results f e for the 1st shot with the SU command sugain qbal 1 lt su DENISE_MARMOUSI_y su shotl itl suximage Generate a sub directory for the data in DENISE trunk par su mkdir MARMOUSI spike CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL 66 To move the seismograms from DENISE trunk par su to DENISE runk par su MARMOUSL spike modify the shell script move_x_y sh in the par directory bin csh set 1 while x lt 101 mv su DENISE_MARMOUSI_x su sho
78. line The meaning of the different parameters is described in the following Domain decomposition Domain Decomposition number_of_processors_in_x direction_ NPROCX 2 number_of_processors_in_y direction_ NPROCY 2 Parallelization is based on domain decomposition see Figure 6 1 i e each processing element PE updates the wavefield within his portion of the grid The model is decomposed by the program into sub grids After de composition each processing elements PE saves only his sub volume of the grid NPROCX and NPROCY spec ify the number of processors in x y direction respectively Figure 6 1 The total number of processors thus is NP NPROCX NPROCY This value must be specified when starting the program with the mpirun command mpirun np lt NP gt bin DENISE DENISE inp see section 5 4 If the total number of processors in DENISE inp and the command line differ the program will terminate immediately with a corresponding error message Obviously the total number of PEs NPROCX NPROCY used to decompose the model should be less equal than the total number of CPUs which are available on your parallel machine If you use LAM and decompose your model in more domains than CPUs are available two or more domains will be updated on the same CPU the program will not terminate and will produce the correct results However this is only efficient if more than one processor is available on each node In order to reduce the amount
79. me that is expanded by the following extensions P wave model vp S wave model vs density model rho In the example above the model files thus are model test vp P wave velocity model model test vs S wave velocity model and model test rho density model In these files each material parameter value must be saved as 32 bit 4 byte native float Velocities must be in meter second density values kg m The fast dimension is the y direction See src readmod c The number of samples for the entire model in the x direction is NX the number of values in the y direction is NY The file size of each model file thus must be NX NY 4 bytes You may check the model structure using the SU command ximage ximage nl lt NY gt lt model test vp Itis also possible to read Qp and Qs grid files to allow for spatial variable attenuation For this you must uncom ment a few lines in src readmod c and generate the corresponding binary files If READMOD O the model is generated on the fly by DENISE i e it is generated internally before the time loop starts If READMOD 0 this function is called in DENISE c and therefore must be specified in src Makefile at the top of src Makefile see section 5 3 If you change this file for example to change the model structure you need to re compile DENISE by changing to the src directory and make denise Q approximation E Q approximation Number of relaxation mecha
80. n the narrow valley and the slow convergence speed near the minimum 1 Fletcher Reeves Fletcher and Reeves 1964 dm 6m FR _ n n Sn mr 3 54 2 Polak Ribiere Polak and Ribi re 1969 m dm 6m _1 PR 3 55 9 mI Sm en 3 Hestenes Stiefel Hestenes and Stiefel 1952 dm m gie se 3 56 m 1 I use the very popular choice 0 BPR which provides an automatic direction reset This is important because subsequent search directions lose conjugacy requiring the search direction to be reset to the steepest descent direction Note that the conjugate gradient method doesn t require any additional computational time because only the gradient m at two subsequent iterations has to be known The application of the nonlinear conjugate gradient method combined with the variable step length calculation to the Rosenbrock function is shown in Fig 3 7 The criss cross pattern of the steepest descent method has vanished and the conjugate gradient method already converges after 30 iterations compared with 200 iteration steps of the pure gradient method CHAPTER 3 THE ADJOINT PROBLEM 31 Residual Energy E 250 200 150 100 Figure 3 7 Results of the convergence test for the Rosenbrock function using the nonlinear conjugate gradient method where the optimum step length is calculated with the parabolic fitting algorithm The minimum is marke
81. nal 191 1 325 345 2012 D Komatitsch and R Martin An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation Geophysics 72 5 155 167 2007 A Kurzmann Applications of 2D and 3D full waveform tomography in acoustic and viscoacoustic complex media PhD thesis Karlsruhe Institute of Technology KIT 2012 available at http nbn resolving de urn nbn de swb 90 344211 A Levander Fourth order finite difference P SV seismograms Geophysics 53 11 1425 1436 1988 Martin Wiley and Marfurt Marmousi2 An elastic upgrade for Marmousi The Leading Edge 25 156 166 2006 P Moczo J Kristek and L Halada The Finite Difference Method for Seismologists An Introduction Comenius University Bratislava 2004 P Mora Nonlinear two dimensional elastic inversion of multioffset seismic data Geophysics 52 1211 1228 1987 P Morse and H Feshbach Methods of theoretical physics McGraw Hill Book Company New York 1953 J Nocedal and S Wright Numerical Optimization Springer New York 2006 A Pica J Diet and A Tarantola Nonlinear inversion of seismic reflection data in a laterally invariant medium 55 3 284 282 1990 R E Plessix and W A Mulder Frequency domain finite difference amplitude preserving migration Geophysical Journal International 157 3 975 987 2004 E Polak and G Ribi re sur la convergence de m thodes de di
82. nisms L 0 L Relaxation frequencies FL 5 0 Tau value for entire model TAUO 0 00001 These lines may be used to define an overall level of intrinsic viscoelastic attenuation of seismic waves In case of L 0 a purely elastic simulation is performed no absorption The frequency dependence of the intrinsic Quality factor Q w is defined by the L relaxation frequencies FL fj 27 1 and one value see equation 5 in Bohlen 2002 For a single relaxation mechanism L 1 Q 2 7 Bohlen 1998 Blanch et al 1995 Bohlen 2002 If the model is generated on the fly the value of TAU can be assigned to all gridpoints for both P and S waves Thus intrinsic attenuation is homogeneous and equal for P and S waves Q w Qs w However it is possible to simulate any spatial distribution of absorption by assigning the gridpoints with different Q values by reading external grid files for Q P waves and Qs S waves see src readmod c or by generating these files on the fly see section 22 Small Q values lt 50 may lead to significant amplitude decay and velocity dispersion Please note that due to dispersive media properties the viscoelastic velocity model is defined for the reference frequency only In denise this reference frequency is specified as the center source frequency At the exact reference frequency elastic and viscoelastic models are equivalent As a consequence slightly smaller and larger minimum an
83. o computational restrictions the original Marmousi II model could not be used because the very low S wave velocities in the sediments would require a too small spatial sampling of the model Therefore new S wave velocities are calculated so that the Poisson ratio is not larger than 0 25 so the soft seabed is replaced by a hard seabed Sears et al 2008 and Brossier et al 2009 have shown the difficulties associated with soft seabed environments for elastic FWT Additionally the size of the Marmousi II model is reduced from 17 km x 3 5 km to 10 km x 3 48 km figure 7 2 62 CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL 63 Marmousi2 Geology 2 4 6 8 10 12 14 16 x km Figure 7 1 Marmousi2 model geology CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL 64 Depth km EL TRE 6 7 Distance km Depth km Distance km Density kg m 2500 2000 Depth km 1500 Distance km Figure 7 2 The reduced and modified complex Marmousi2 model used for the elastic FWT CHAPTER 7 EXAMPLE 1 THE ELASTIC MARMOUSI2 MODEL 65 7 1 1 Acquisition geometry and FD model The acquisition geometry consists of a Ocean Bottom Cable OBC located on the seafloor at depth of 460 m below the free surface The OBC contains 400 two component geophones with a spatial spacing of 20 m recording the particle velocities v For the synthetic dataset 100 airgun shots are excited The sources are towed
84. of data that needs to be exchanged between PEs you should decompose the model into more or less cubic sub grids In our example we use 2 PEs in each direction NPROCX NPROCY 2 The total number of PEs used by the program is NPROC NPROCX NPROCY 4 43 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 44 X free surface or absorbing boundary NX Figure 6 1 Geometry of the numerical FD grid using 2 processors in x direction NPROCX 2 and 2 processors in y direction NPROCY 2 Each processing element PE is updating the wavefield in its domain At the top of the numerical mesh the PEs apply a free surface boundary condition if FREE SURF I otherwise an absorbing boundary condition PML The width of the absorbing frame is FW grid points The size of the total grid is NX grid points in x direction and NY gridpoints in y direction The size of each sub grid thus is NX NPROCX x NY NPROCY gridpoints The origin of the Cartesian coordinate system x y is at the top left corner of the grid CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 45 Order of the FD operator FD order Order of FD coefficients values 2 4 12 FD_ORDER 2 Maximum relative group velocity error E minimum number of grid points per shortest wavelength is defined by FD_ORDER values 0 Taylor coefficients 1 Holberg coeff t E D 2 0 5 3 1 0 4 3 0 m
85. ompile the additional library for timedomain filtering In the DENISE libcseife directory simply use the shell script bash 2 05b DENISE libcseife gt make It is necessary to preinstall FFTW Fastest Fourier Transform in the West http www fftw org The source code of DENISE is located in the directory DENISE src To compile DENISE the name of the model function has to be entered in the MAKEFILE Depending on your MPI environment MPI distribution you may need to modify the compiler options in src Makefile For a few typical platforms the compiler options are available in src Makefile It is often useful to enable a moderate level of optimization typically 03 The highest level of optimization O4 can lead to a strong performance improvement For example the optimization option 04 of the hcc LAM compiler leads to a speedup of DENISE of approximately 30 per cent Even though keep in mind that O4 can also lead to crashes and compilation errors when used in combination with certain compilers Linbin Zhang suggested to use Advanced Vector Extensions AV X on supported CPUs by using the following compiler flags LFLAGS 1m lcseif fir isto CFLAGS O3 xAVX ipo fno fnalias restrict Depending e g on the size of the inversion problem used spatial FD operators a reduction of computation time of up to 20 could be measured No other changes in the Makefile should be necessary Ma
86. on problem 18 CHAPTER 3 THE ADJOINT PROBLEM 19 3 2 How to find an optimum model Figure 3 2 shows a schematic sketch of the residual energy at one point in space as a function of two model parameters The colors represent different values of the residual energy Red areas represent models with high residual energy which do not fit the data while the blue parts are good fitting models with low residual energies The aim is to find the minimum of the residual energy marked by the red cross Starting at a point m3 A1 x 1 in the parameter space we want to find the minimum by updating the material parameters in an iterative way m m 3 4 along the search direction m with the step length To find the optimum search direction m we expand the residual energy E m near the starting point in a Taylor series 1 1 1 and set the derivative of Eq 3 5 with respect to m zero Q E m m OE EN _ 3m om E PS Which finally leads to EN 1 9 im da 22 in where OE Om denotes the steepest descent direction of the objective function and H the inverse Hessian matrix The inverse Hessian matrix for the elastic problem is often singular and can only be calculated with high computational costs Therefore the inverse Hessian matrix is approximated by a preconditioning operator P There is channel ch
87. onitor the modeling process You may want to save this screen info to an output file by adding 7 DENISE out or tee DENISE out to your starting command If LOG 1 all other processes with PE number PEno greater than zero will write their information to LOG FILE PEno If you specify LOG 2 PE 0 will also output information to LOG FILE 0 As a consequence only little information is written directly to the screen of your shell On supercomputers where you submit modeling jobs to a queuing system as batch jobs LOG 2 may be advantageous In case of LOG 2 you may still watch the simulation by checking the content of LOG FILE 0 for example by using the Unix commands more or tail After finishing the program the timing infor mation is written to the ASCII file log test log 0 timings This feature is useful to benchmark your local PC cluster or supercomputer If LOG 0 no output from node 0 will be written neither to stdout nor to an LOG file There will be also no output of timing information to the ASCII file log test log 0 timings Checkpointing Checkpoints read_wavefield_from_checkpoint_file_ yes 1 no 0 _ CHECKPTREAD 0 save wavefield to checkpoint file yes 1 no 0 CHECKPTWRITE 0 checkpoint file CHECKPTFILE tmp checkpoint DENISE These options are obsolete and are not be supported in the current version of DENISE General inversion parameters General DENISE inversion parame
88. opers here in the future License DENISE is free software you can redistribute it and or modify it under the terms of the GNU General Public License as published by the Free Software Foundation version 2 0 of the License only DENISE is distributed in the hope that it will be useful but WITHOUT ANY WARRANTY without even the implied warranty of MERCHANTABIL ITY or FITNESS FOR A PARTICULAR PURPOSE See the GNU General Public License for more details You should have received a copy of the GNU General Public License along with DENISE See file COPYING and or http www gnu org licenses gpl 2 0 html The authors of DENISE are listed in the AUTHORS section Contents 1 Introduction Di Citation uio eee So xps ERE Ita Tu REN SS ES 1 2 Support s S S muse Dede Pega eth OOS ped RES 2 Theoretical Background 2 1 Equations of motion for an elastic medium 0 2 2 Solution of the elastic wave equation by finite differences 2 2 1 Discretization of the Wave equation ee 2 2 2 Accuracy of FD operators s ku nn rn ee 2 2 3 Initial and Boundary Conditions 2 3 Numerical Artefacts and 2 3 1 Gnd Dispersion ee he he a 23 2 The Courant Instability 54 mos A Kuren ara 3 The adjoint problem 3 1 Whatisan op imum 12
89. or iteration steps n gt 2 c mP B cy 1 with dc dm 3 58 where the weighting factor cu 3 59 by Polak Ribi re is used The convergence of the Polak Ribi re method is guaranteed by choosing max GP 0 g Estimate the step length un by the line search algorithm described in chapter 3 6 h Update the material parameters using the gradient method Mn 1 my Un CH 3 60 If the material parameters are not coupled by empirical relationships it is important to update all three elastic material parameters at the same time otherwise strong artifacts may dominate the inversion result especially in the case of very complex media 3 Iftheresidual energy E is smaller than a given value stop the iteration Otherwise continue with the next iteration step Chapter 4 Source Wavelet Inversion So far we assumed that only the material parameters are unkown while the characteristics of the sources are perfectly known For the application of FWI to a field dataset the source wavelet has to be estimated In frequency domain the source wavelet s has the complex value s e if where 1 y 1 e and f are the real and imaginary parts respectively The seismograms of the vertical displacements of the modelled data can be described by vr Cy r idy r e if where id y denotes the spike response of the vertical displacement and the receiver location Similar the seismograms of the vertical displa
90. osest next grid point It is not yet possible to output seismograms for arbitrary receiver locations since this would require a certain wavefield interpolation CHAPTER 6 DEFINITION OF PARAMETERS FOR THE MODELING AND INVERSION CODE 52 It is important to note that the actual receiver positions defined in REC_FILE or in DENISE inp may vary by DX 2 and or DY 2 and or DZ 2 due to the staggered positions of the particle velocities and stress tensor components In our example we specify 100 receiver location Due to the small size of the model most of the specified receiver positions will be located inside this absorbing boundary if ABS 2 see Chapter 6 1 A corresponding warning mes sage will appear If you choose to read the receiver location from REC_FILE receiver receiver dat READREC 1 only 10 receivers locations are considered The list of receivers specified in file receiver receiver dat is equivalent to the parameters in the input file with only a larger distance between adjacent receivers NGEOPH 10 Towed streamer Towed streamer parameters for towed streamer acquisition The first N STREAMER receivers in REC FILE belong to streamer 31 Cable increment per shot REC INCR X 80 0 Cable increment per shot REC INCR Y 0 A streamer geometry can be defined by moving the receiver positions defined in the receiver file REC_FILE in accordance with the sources in x direction by
91. ould occur 03 sources c I libesei te mpice 03 e solvelin o I liboseife 03 apat_filt 1 libeseife pios 07 e BpliESTC E Fa flihesseire pie 03 splitsro back I liboseife npice 03 c splitrec c I libeoseife mias eco stalta c I liboseife mpicc 03 c step length est c I libcseife mpicc 03 e step length estl c I libcseife mplec 0S g These Te mpice 03 taper ce I liboseife mpice 03 c taper grad c I liboseife ples 03 6 tapes grad shot o I liboseife mpice 03 timedomain filt c I libeseife mpicc 03 c time window c I libcseife mpios 03 util c 1 libeoselfe CHAPTER 5 GETTING STARTED 40 03 c wavelet c I libcseife 03 wavelet_stf c 1 libcseife picc 03 c writemod c I libcseife pice 03 c write par c I libeseife pios 03 c writedsk o I libeseife picc 03 c zero fdveps c I libcseife picc L libcseife denise o calc mat change o calc mat change test o calc res o calc opt step o calc opt step test o calc energy o catseis o checkfd ssg elastic o conv FD o psource o holbergcoeff o comm ini o exchange v o exchange s o exchange L2 0 fft o fft filt o forward mod o snap ssg o seismo ssg o surface elastic 2nd o writemod o write par o writedsk o zero fdveps o o bin denis lm
92. quite complicated an additional constraint is applied during the inversion To stabilize the inversion possible density values are restricted between 1000 kg m and 3000 kg m using hard constraints and the density step length calculated during the model update systematically reduced by a factor SCALERHO 0 5 Otherwise geophysically unrealistic density values might occur in the model To increase the convergence of FWI the inverse Hessian is approximated by a simple linear scaling of the gradient with depth For the Marmousi II model the inversion is separated into four parts which cover different frequency ranges with maximum frequencies of 2 5 10 and 20 Hz respectively This inversion strategy is assembled in the workflow file FWI workflow marmousi inp Before running the FWI change the following parameters in the file DENISE marm OBC inp 1 1 MFILE start marmousi_II_smooth2 0 which defines the initial model and switches the DENISE code from modeling to inversion mode Run the inverse code on 15 cores of the cluster by typing is mpirun np 15 bin denise DENISE marm OBC inp FWI workflow marmousi inp or submit a job file on the NEC cluster The whole inversion requires approximately 1 day The final inversion results after 285 iterations for the velocity parametrization are shown in Fig 7 5 Additionally depth profiles at 1 3 5 km and xy 6 4 km of the starting model and inversion result are
93. rections conjugu es Revue Francaise d Informatique et de Recherche Op rationnelle 16 35 43 1969 R Pratt Inverse theory applied to multi source cross hole tomography Part II Elastic wave equation method Geo physical Prospecting 38 311 329 1990 R Pratt Seismic waveform inversion in the frequency domain Part 1 Theory and verification in a physical scale model Geophysics 64 888 901 1999 R Pratt Velocity models from frequency domain waveform tomography Past present and future In 66th EAGE conference and exhibition Expanded Abstracts pages 181 182 Paris France 2004 R Pratt and M Worthington Inverse theory applied to multi source cross hole tomography Part I Acoustic wave equation method Geophysical Prospecting 38 287 310 1990 BIBLIOGRAPHY 78 R Pratt Gao Zelt and A Levander The limits and complementary nature of and waveform tomogra phy In International Conference of Sub basalt imaging Expanded Abstracts pages 181 182 Cambridge England 2002 R Pratt L Huang N Duric and P Littrup Sound speed and attenuation imaging of breast tissue using waveform tomography of transmission ultrasound data 2007 J Robertsson J Blanch and W Symes Viscoelastic finite difference modeling Geophysics 59 9 1444 1456 1994 J Robertsson A Levander W Symes and K Holliger A comparative study of free surface boundary conditions for finite difference simulation of
94. rock function The minimum is marked with ared point the starting position with a blue point The maximum number of iterations is 16000 The step length jj varies between 2e 3 top and 6 1 3 bottom CHAPTER 3 THE ADJOINT PROBLEM 29 Normalized L2 Norm inimum of the parabo Aur step lengtl Step length u Figure 3 5 Line search algorithm to find the optimum step length opt The true misfit function yellow line is approximated by a parabola fitting values of the objective function for 3 different step length Because changes of the density model are in most cases smaller than velocity changes the step length for the density update be systematically reduced by a factor s no All material parameters can be updated simultaneously or ac cording to a hierachical strategy To save computational time the corresponding L2 norms are calculated for a few representative shots The number of shots depends on the complexity of the problem and the signal noise ratio of the data For the acoustic case the step length estimation by parabolic fitting works very well and leads to a smooth de crease of the misfit function during the FWT Kurzmann 2007 personal communication For the multiparameter elastic FWT the misfit function consists of more local minima and therefore the decrease of the objective function is not as smooth as in the acoustic case Brossier 2009 proposed a more intensive bracketing stage before applying th
95. s PICKS_FILE _sourcenumber dat So the number of sources in SRCREC must be equal to the number of files Each file must contain the picked times for every receiver Other important parameters are set in the workflow file see section 6 2 Name of the misfit log file deeem MISFIT LOG FILE log_file_for_misfit_evolution_ MISFIT_LOG_FILE LOG_TEST dat The name of the misfit log file can be changed with the parameter MISFIT_LOG_FILE The columns of the misfit log file contain information about the step length and misfit function values acquired during the step length estimation and the stage number nstage opteps epst 1 epst 2 epst 3 L2t 1 L2t 2 L2t 3 L2t 1 nstage When a frequency filter is applied information about the corner frequencies are also written to the misfit log file opteps epst 1 epst 2 epst 3 L2t 1 L2t 2 L2t 3 L2t 1 FC low FC high nstage 2D Gaussian gradient smoothing filter Definition of smoothing the Jacobians with 2D Gaussian apply spatial filtering yes 1 GRAD FILTER 0 filter length in gridpoints FILT SIZE GRAD 10 If GRAD_FILTER 1 smooth the gradients for the different material parameters using a 2D Gaussian with a filter length defined by FILT_SIZE_GRAD in gridpoints Time lapse FWI mode FWT double difference time lapse mod activate time lapse mode yes 1 TIMELAPSE 0
96. s but also the amplitudes which contain information on the distribution of the elastic material parameters in the underground To achieve this goal three problems have to be solved 1 What is an optimum model 2 How can this model be found 3 Is this model unique or are other models existing which could explain the data equally well 3 1 What is an optimum model In reflection seismics the i component of the elastic displacement field ui Xs Xr t excited by sources located at x will be recorded by receivers at x at time t For a given distribution of the material parameters the forward problem Eq 2 3 can be solved by finite differences section 2 2 The result is a model data set u d This modelled data be compared with the field data 998 If the misfit or data residuals du 17194 1958 figure 3 1 between the modelled and the field data is small the model can explain the data very well If the residuals are large the model cannot explain the data The misfit can be measured by a vector norm 1 which is defined for p 1 2 as 1 Llp 3 1 The special case L is defined as max du P 3 2 The L2 norm 1 E du 3 3 has special physical meaning It represents the residual elastic energy contained in the data residuals An optimum model can be found in a minimum of the residual energy Therefore the optimum model is the solution of a nonlinear optimizati
97. sily calculated by summing over the weighting coefficients 5 2 22 h 5 B 2 23 In table 2 3 h is listed for different FD operator lengths and types Taylor and Holberg operators Criterion 2 22 is called Courant Friedrichs Lewy criterion Courant et al 1928 Courant et al March 1967 figure 2 4 shows the evolution of the pressure field when the Courant criterion is violated After a few time steps the amplitudes are growing to infinity and the calculation becomes unstable FDORDER h Taylor h Holberg 2nd 1 0 1 0 4th 7 6 1 184614 6th 149 120 1 283482 8th 2161 1680 1 345927 10th 53089 40320 1 387660 12th 1187803 887040 1 417065 Table 2 3 The factor in the Courant criterion for different orders 2nd 12th and types Taylor and Holberg of FD operators CHAPTER 2 THEORETICAL BACKGROUND 17 T 0 8ms T 1 5ms 3 0 17 0 27 0 77 4 0 8 1 0 8 J 0 0 2 0 4 0 6 0 8 1 0 0 2 0 4 0 6 0 8 1 x m x m 2 3ms T 3 0ms 0 0 2 0 4 0 6 0 8 1 0 0 2 0 4 0 6 0 8 1 x m x m Figure 2 4 Temporal evolution of the Courant instability In the colored areas the wave amplitudes are extremly large Chapter 3 The adjoint problem The aim of full waveform tomography is to find an optimum model which can explain the data very well It should not only explain the first arrivals of specific phases of the seismic wavefield like refractions or reflection
98. specify more PEs in the FD software but all processes will be executed on the same CPU Thus this only makes sense if you run the software on a multicore system In this case you should boot it without a lamhosts file and specify a total number of 2 processing elements PEs To shut down LAM again before logout use the command amclean v However it is not necessary to shut down restart LAM after each logout login 35 CHAPTER 5 GETTING STARTED 36 5 1 2 How to run DENISE on NEC Linuxcluster at RZ Kiel Before you can run DENISE on the Linux cluster at the computing centre in Kiel you have initialize Intel MPI and Intel compilers and assure that the different nodes can communicate password free This has to be done only once 1 Add the following lines to your bashrc in your HOME directory to intialize Intel MPI and the Intel compiler opt intel composer xe 2013 spl bin compilervars sh intel64 opt intel impi 4 1 1 036 intel64 bin mpivars sh 2 To setup a password free communication between the different nodes generate a pair of authentication keys for ssh with sungwXXX nesh fe2 ssh keygen t dsa You can accept the default values by hitting return 3 Copy the file SHOME ssh id_dsa pub to SHOME ssh authorized_keys Because DENISE can produce up to a few GB of data output don t run the code from the home directory To submit a batch job it is required that DENISE is located in the WORK directory Ke
99. stic medium can be described by a system of coupled linear partial differential equations They consist of the equations of motion Ov ij P 0 OX f 2 1 which simply state that the momentum of the medium the product of density p and the displacement velocity v can be changed by surface forces described by the stress tensor oj or body forces fi These equations describe a general medium like gas fluid solid or plasma The material specific properties are introduced by additional equations which describe how the medium reacts when a certain force is applied In the isotropic elastic case this can be described by a linear stress strain relationship 2uei u 1 Ou 2 2 22 2 OX Ox ij where and are the Lam parameters ei the strain tensor and u the displacement Using v om 2 1 and 2 2 can be transformed into a system of second order partial differential equations Pu EE 3 Ot B OX Oij A00 276 2 3 1 Ou d Qu Ej E 2 OX This expression is called Stress Displacement formulation Another common form of the elastic equations of motion can be deduced by taking the time derivative of the stress strain relationship and the strain tensor in Eq 2 3 Since the Lam parameters and yu do not depend on time Eq 2 3 can be written as _ p Ot 3 OX 00 2 4 a ap t
100. t x itl su MARMOUSI spike DENISE MARMOUSI x su shot x mv su DENISE MARMOUSI y su shot x itl su MARMOUSI spike DENISE MARMOUSI y su shot x 4 mv su DENISE MARMOUSI p su shot x itl su MARMOUSI spike DENISE MARMOUSI p su shot x set x expr x 1 end Fig 7 3 shows the development of the pressure wavefield excited by shot 50 for the central part of the complex elastic Marmousi2 model at 6 different time steps The P wave is traveling from the source through the water column 100 0 ms and is reflected at the seafloor T 400 0 ms In the elastic subseafloor medium the wavefield becomes very complex The layers in the steep thrust fault system produce numerous reflections and internal multiples T 600 0 ms Additionally strong diffracted waves are generated at the sharp corners of the thrust faults between the disturbed high velocity sediment blocks within the thrust faults and the surrounding low velocity sediments At the free surface strong multiple reflections occur 800 0 ms The wavefront of the direct wave is quite deformed due to strong velocity contrasts within the thrust fault system After 1500 ms nearly all kinds of waves which can be found in the literature are present Reflections refractions diffractions internal multiples or interface waves The trapped gas sand reservoirs C2 and produce strong reflections and mode conversions This complexity is also visible in the seismic section
101. ted to the mapping of small changes from the data to the model space and vice versa figure 3 3 A small change in the model space m e g one model parameter at one point in space will result in a small Model Space x coordinate gt x coordinate gt in Y fa Y s si sources receiver Figure 3 3 Mapping between model and data space and vice versa perturbation of the data space e g one wiggle in the seismic section If the Frech t derivative 2 om is known all the small perturbations in model space can be integrated over the model volume V to calculate the total change in data space Tarantola 20051 Xy t f 2 5m 3 13 or introducing the linear operator ie L m dV sm V om In a similar way small changes in the data space can be integrated to calculate the total change in the model space m Tarantola 2005 fer E3 3 14 om sources receiver or as operator equation L dw CHAPTER 3 THE ADJOINT PROBLEM 22 In this case the Frech t derivative gu is replaced by it s adjoint counterpart Zu Note that du 4 and dm 6m so there is no unique way to map perturbations from the model to the data space or vice versa Because the operator L is linear the kernel of L and it s adjoint counterpart L are identical see chapter 5 4 2 in Tarantola 2005 u m Therefore the mapping
102. ters number_of_TDFWI_iterations_ ITERMAX 100 output_of_jacobian_ JACOBIAN jacobian jacobian_Test seismograms_of_measured_data_ DATA_DIR su plexiglas DENISE_plexiglas cosine_taper_ yes 1 no 0 _ TAPER 0 taper_length_ in_rec_numbers _ TAPERLENGTH 5 Inverse_Type_ gradient 1 complete 2 _ INVTYPE 2 gradient_taper_geometry_ GRADT1 GRADT2 GRADT3 GRADT4 5 10 510 520 type_of_material_parameters_ Vp Vs rho 1 Zp 2s rho 2 lam mu rho 3 INVMAT1 forward_modelling_only_ yes 10 _FWT_ yes 0 INVMAT 0 point source backpropagation x and y comp 1 y comp 2 x comp 3 p comp 4 QUELLTYPB testshots for step length est TESTSHOT START TESTSHOT END TESTSHOT INCR 1 21 10 This section covers some general inversion parameters The maximum number of iterations are defined by ITER MAX The switch INVMAT controls if only the forward modeling code should be used INVMAT 10 e g to calculate synthetic seismograms or a complete FWT run INVMAT 0 In case of INVMAT 10 the parameters in the workflow file section 6 2 are ignored but a workflow file still has to be defined The seismic sections of the real field data are located in the DATA DIR As noted in section 3 4 the gradients can be expressed for different model parameterizations The switch INVMATI defines which parameterization should be used seismic velocities and density Vp Vs rho INVMATIz
103. wide aperture seismic data part 2 numerical examples and scalability analysis Computer amp Geosciences 35 496 5 14 2009b A Tarantola Inversion of seismic reflection data in the acoustic approximation Geophysics 49 1259 1266 1984a A Tarantola Linearized inversion of seismic reflection data Geophysical Prospecting 32 998 1015 1984b A Tarantola A strategy for nonlinear elastic inversion of seismic reflection data Geophysics 51 1893 1903 1986 A Tarantola Theoretical background for the inversion of seismic waveforms including elasticity and attenuation 128 365 399 1988 A Tarantola Inverse Problem Theory SIAM 2005 R Versteeg The marmousi experience Velocity model determination on a complex data set The Leading Edge 13 927 936 1994 J P SV wave propagation in heterogeneous media velocity stress finite difference method Geophysics 51 4 889 901 1986
104. wing sections we give a short description of the different modeling parameters options and how the program is used in a parallel MPI environment 5 1 Requirements The parallelization employs functions of the Message Passing Interface MPI MPI has to be installed when compiling and running the DENISE software At least two implementations exist for Unix based networks OpenMPI MPICH2 and Intel MPI The LAM MPI implementation is no longer supported by the developers However currently all four implementations work with DENISE OpenMPI MPICH2 and Intel MPI are MPI programming environments and development systems for heterogeneous computers on a network As of the time of writing we get the best performance out of DENISE by using Intel MPI together with the latest Intel Compiler on a NEC Linux Cluster With MPI a dedicated cluster or an existing network computing infrastructure can act as a parallel computer Fast network infiniband connections and RAM access are the most important issuses for a good scaling of the DENISE code The latest version of OpenMPI can be obtained from http www open mpi org MPICH2 is available at http www unix mcs anl gov mpi mpich LAM MPI can be downloaded here http www lam mpi org the commerical Intel MPI from here https software intel com en us intel mpi library 5 1 1 LAM Even though outdated LAM MPI will be used to illustrate the setting up of the MPI implementation In order to reproduce the followin
105. y used test problem for seismic imaging techniques Beside the original acoustic version of the model an elastic version was developed by Martin et al 2006 This model contains both simple approximately 1D and complex geological structures In the following the performance of the FWT code will be tested for the complex part of a modified Marmousi II model using a parametrization with seismic velocities 7 1 complex Marmousi2 model The Marmousi2 model Fig 7 1 consists of a 500 m thick water layer above an elastic subseafloor model The sediment model is very simple near the left and right boundaries but rather complex in the centre At both sides the subseafloor is approximately horizontally layered while steep thrust faults are disturbing the layers in the centre of the model Embedded in the thrust fault system and layers are small hydrocarbon reservoirs figure 7 1 Martin et al 2006 One shallow gas sand in a simple structural area One relatively shallow oil sand in a structural simple area B Four faulted trap gas sands at varying depths C1 C2 C3 C4 Two faulted trap oil sands at medium to deep depths D1 D2 One deep oil and gas sand anticlinal trap 1 2 Water wet sand The deeper parts of the model consist of salt and reef structures The thrust fault system and the reef structures are not easy to resolve by conventional first arrival tomography so it is an ideal test model for the FWT Due t

DENISE - Black Edition User Manual - Christian

Contents

Download Pdf Manuals

Related Search

Related Contents