Seeking Advice on FWI using DENISE with OBN Field Data: Challenges and Simulation Instability

Hello everyone @daniel-koehn @ovcharenkoo @ADharaUTEXAS123007 @pplotn, I hope you’re doing well. I would like to seek your advice and assistance regarding Full Waveform Inversion (FWI) using DENISE with OBN field data. I have learned that you all are familiar with and have expertise in DENISE and FWI, which is why I am reaching out for your guidance. I have encountered some challenges and would greatly appreciate any assistance you can provide. Thank you in advance for your help!

Here is a brief overview of what I have accomplished in the past month:

First, I selected a source line that closely aligns with the receiver line, ensuring that the receiver interval remains constant. Then, I projected the receiver line onto the source line.
I performed a 3-component pre-stack time migration using the OBN data. During this process, I obtained Vp and Vs values in the time domain, which I subsequently converted into the depth domain. Here is the overview of my acquisition (to save time and for testing purposes, I used only 16 shots).
I utilized the raw data, meaning that I did not apply any conventional seismic preprocessing to my data (though I’m uncertain if this is the correct approach).
To generate synthetic data, I used the velocity field and density (rho) field. I’m attaching both my synthetic data and the actual field data for your reference. I observed that the Marmousi example also includes a water layer. However, when I set the water layer in the Vs field, my water depth is approximately 80 meters (same as the receiver positions), while my source position is at around 6 meters. Unfortunately, I couldn’t obtain the desired synthetic data until I moved the source position to 12.5 meters.

During my attempts to achieve FWI using DENISE, the first epoch performed well. However, upon entering the second epoch, I encountered an error stating “Simulation is unstable.” I understand that the reason behind this issue is likely due to my lack of experience with FWI, specifically in obtaining the gradient field. As a newcomer to this field, I am struggling to analyze the exact cause of this occurrence.


*******************************************************************************
This is program DENISE Black-Edition                                           
Parallel 2-D elastic Finite Difference FWI code                                

Forward/FWI/RTM codes written by D. Koehn and D. De Nil                        
2D isotropic PSV forward code partly based on FDVEPS written by T. Bohlen      
Institute of Geosciences, Kiel University, Germany                           

See README.md file and LICENSE.md for redistribution conditions.               
*******************************************************************************

**Message from check_model_phys (printed by PE 0):

-----------------------  DENISE operation mode  ----------------------
MODE=1: Time-domain FWI is applied.

-----------------------  DENISE Physics  ----------------------
PHYSICS=1: Solve 2D isotropic elastic PSV problem.

This is the log-file generated by PE 0 

**Message from initprocs (printed by PE 0):
Size of subarrays in gridpoints:
IENDX= 628
IENDY (vertical) = 250

**Message from initprocs (written by PEEEP 0):
Processor locations in the 2D logical processor array
MYID    POS(1):left,right   POS(2): top, bottom
0       0: 2,1          0: 9,3 

**Message from write_par (printed by PE 0):

------------------------- Processors ------------------------ Number of PEs in x-direction (NPROCX): 3 Number of PEs in vertical direction (NPROCY): 4 Total number of PEs in use: 12

----------------------- Discretization --------------------- Number of gridpoints in x-direction (NX): 1884 Number of gridpoints in y-direction (NY): 1000 Grid-spacing (DH): 6.250000e+00 meter Time of wave propagation (T): 4.000000e+00 seconds Timestep (DT): 1.000000e-03 seconds Number of timesteps: 4000

------------------------- FD ORDER ----------------------------- FDORDER = 8 MAXRELERROR = 1

------------------------- SOURCE ----------------------------- reading source positions, time delay, centre frequency and initial amplitude from ASCII-file ./source/OBNtest.dat

wavelet of source: 1st derivative of Gaussian

------------------------- RECEIVER -------------------------- reading receiver positions from single file ./receiver/OBN_ReceiverFile.dat

reference_point_for_receiver_coordinate_system: x=0.000000 y=0.000000 z=0.000000 ------------------------- FREE SURFACE ------------------------ free surface at the top of the model !

------------------------- CPML --------------------- width of absorbing frame is 10 gridpoints. CPML damping applied. Damping velocity in the PML frame in m/s: 1500.000000 . Frequency within the PML frame in Hz: 20.000000 npower: 4.000000 k_max: 1.000000 No periodic boundary condition. ------------------------- MODEL-FILES ------------------------- names of model-files: shear wave velocities: start/OBN_model.vs tau for shear waves: start/OBN_model.ts density: start/OBN_model.rho compressional wave velocities: start/OBN_model.vp tau for P-waves: start/OBN_model.tp

------------------------- Q-APROXIMATION -------------------- Number of relaxation mechanisms (L): 0 The L relaxation frequencies are at:
Hz Value for tau is : 0.000010

----------------------- SEISMOGRAMS ---------------------- seismograms of x- and y-component of particle velocity. output-files: su/OBN_RealData_x.su su/OBN_RealData_y.su The data is written in IEEE SU-format . samplingrate of seismic data: 0.001000 s Number of samples per trace: 4000

----------------------- DENISE elastic specific parameters ----------------------

Maximum number of iterations: 100 location of the measured seismograms : su/OBNTest/OBN_FinalRealdata

INVMAT1=1: Inversion parameters are vp, vs and rho. QUELLTYPB=1: Inversion of x and y component.

Shots used for step length estimation: TESTSHOT_START = 1 TESTSHOT_END = 31 TESTSHOT_INCR = 10

Cosine Taper used : 0 Log file for misfit in each iteration step: LOG_TEST.dat

Output of inverted models to: model/DFB/modelTest

--------------- Gradient tapering ------------------- SWS_TAPER_GRAD_VERT=1: Vertical taper applied. (GRADT1=21, GRADT2=25, GRADT3=490, GRADT4=500)

SWS_TAPER_GRAD_HOR=1: Horizontal taper applied. (GRADT1=21, GRADT2=25, GRADT3=490, GRADT4=500, EXP_TAPER_GRAD_HOR=2.000000)

SWS_TAPER_GRAD_SOURCES=0: No taper around the sources applied.

SWS_TAPER_CIRCULAR_PER_SHOT=0: No taper around the sources applied.

--------------- Gradient smoothing with 2D-Gaussian filter ------------------- GRAD_FILTER=0: Gradients are not filtered.

--------------- Limits of model parameters ------------------- VPLOWERLIM = 1400.000000 VPUPPERLIM = 5000.000000 VSLOWERLIM = 0.000000 VSUPPERLIM = 4000.000000 RHOLOWERLIM = 0.000000 RHOUPPERLIM = 3000.000000

--------------- Optimization method ------------------- GRAD_METHOD=2: LBFGS NLBFGS=20

--------------- Model smoothing ------------------- MODEL_FILTER=0: vp and vs models are not filtered after each iteration step.

--------------- Trace kill ------------------- TRKILL=0: No trace kill is applied

--------------- Trace normalization ------------------- NORMALIZE=0: No normalization of measured and synthetic seismograms.

--------------- Reduce size of inversion grid ------------------- Every 3 time sample is used for the calculation of the gradients.

--------------- Step length estimation ------------------- EPS_SCALE = 0.020000 STEPMAX = 6 SCALEFAC = 2.000000

--------------- Reverse Time Modelling ------------------- RTMOD=0: No Reverse Time Modelling applied.

--------------- Gravity Modelling and Inversion ------------------- No Gravity Modelling and Inversion applied. GRAVITY=0 Boundary in x-direction [gridpoints] NGRAVB = 500 Boundary in z-direction [m] NZGRAV = 200000 Modelling and Inversion of Gravity Data. GRAV_TYPE=1 Self-defined Model is used as Background Density. BACK_DENSITY=2 background density file: gravity/background_density.rho

Reading receiver positions from single file: ./receiver/OBN_ReceiverFile.dat Message from function receiver (written by PE 0): Number of receiver positions found: 225

**Message from main (printed by PE 0): Size of local grids: NX=628 NY=250 Each process is now trying to allocate memory for: Dynamic variables: 3.13 MB Static variables: 3.76 MB Seismograms: 2.26 MB Buffer arrays for grid exchange: 0.07 MB Storage of forward modelled wavefields: 3994.71 MB Gradients, updated material parameters: 13.77 MB Residual seismograms: 3.39 MB Full seismograms: 10.30 MB Network Buffer for MPI_Bsend: 0.13 MB

Total memory required: 4021.22 MB.

... memory allocation for PE 0 was successfull.

Reading source positions, time-shift, centre frequency and amplitude from file: ./source/OBNtest.dat Number of source positions specified in ./source/OBNtest.dat : 16 Maximum frequency defined in ./source/OBNtest.dat: 2.00e+01 Hz All sources will be modelled individually because of RUN_MULTIPLE_SHOTS=1! Message from function sources (written by PE 0):

...reading model information from modell-files... Vp: start/OBN_model.vp

 Vs:
 start/OBN_model.vs

 Density:
 start/OBN_model.rho

PE 0 is writing model to start/OBN_model.denise.pi.0.0

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to start/OBN_model.denise.pi Opening model files: start/OBN_model.denise.pi.??? ... finished. Copying... ... finished. Use ximage n1=1000 < start/OBN_model.denise.pi label1=Y label2=X title=start/OBN_model.denise.pi to visualize model.

PE 0 is writing model to start/OBN_model.denise.mu.0.0

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to start/OBN_model.denise.mu Opening model files: start/OBN_model.denise.mu.??? ... finished. Copying... ... finished. Use ximage n1=1000 < start/OBN_model.denise.mu label1=Y label2=X title=start/OBN_model.denise.mu to visualize model.

PE 0 is writing model to start/OBN_model.denise.rho.0.0

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to start/OBN_model.denise.rho Opening model files: start/OBN_model.denise.rho.??? ... finished. Copying... ... finished. Use ximage n1=1000 < start/OBN_model.denise.rho label1=Y label2=X title=start/OBN_model.denise.rho to visualize model.

**Message from checkfd (printed by PE 0): Minimum and maximum P-wave and S-wave velocities within subvolumes: MYID Vp_min(f=fc) Vp_max(f=inf) Vs_min(f=fc) Vsmax(f=inf) 0 1.459000e+03 2.090969e+03 0.000000e+00 1.344821e+03 Global values for entire model: Vp_max= 3.251180e+03 m/s Vs_min=0.000000e+00 m/s

------------------ CHECK FOR GRID DISPERSION -------------------- To satisfactorily limit grid dispersion the number of gridpoints per minimum wavelength (of S-waves) should be 6 (better more). Here the minimum wavelength is assumed to be minimum model phase velocity (of S-waves) at maximum frequency of the source devided by maximum frequency of the source. Maximum frequency of the source is approximately 40.00 Hz The minimum wavelength (of S-waves) in the following simulation will be 0.000000e+00 meter. Thus, the recommended value for DH is 0.000000e+00 meter. You have specified DH= 6.250000e+00 meter.

W A R N I N G M E S S A G E: Grid dispersion will influence wave propagation, choose smaller grid spacing (DH).

----------------------- CHECK FOR STABILITY --------------------- The following simulation is stable provided that

 p=cmax*DT/DH < 1/(sqrt(2)*gamma),

where cmax is the maximum phase velocity at infinite frequency and gamma = sum(|FD coeff.|) In the current simulation cmax is 3251.18 m/s .

DT is the timestep and DH is the grid size.

In this simulation the stability limit for timestep DT is 1.009956e-03 seconds . You have specified DT= 1.000000e-03 s. The simulation will be stable.

----------------------- ABSORBING BOUNDARY ------------------------ Width (FW) of absorbing frame should be at least 10 gridpoints. You have specified a width of 10 gridpoints. =========================================== FWI-stage 1 of 4 =========================================== Density is inverted from iteration step 0.

Vp is inverted from iteration step 0.

Vs is inverted from iteration step 0.

Qs is inverted from iteration step 0.

Smoothing (spatial filtering) of the gradients: SPATFILTER=0: Gradients are not smoothed.

--------------- Frequency filtering ------------------- TIME_FILT=1: Time domain filtering is applied Applying FWI with lowpass filter for corner frequency 2.000000 Hz Order of lowpass filter is: 6

--------------- Time windowing and damping ------------------- TIMEWIN=0: No time windowing and damping is applied

--------------- termination of the program ------------------- Misfit change during the last two iterations is smaller than 1.000000 percent.

--------------- Energy preconditioning ------------------- EPRECOND = 3 - Hessian approximation by Plessix & Mulder (2004)

EPRECOND = 3 - energy preconditioning deactivated

--------------- Gradient calculation ------------------- Misfit function: LNORM==1 corresponds to L1 Norm LNORM==2 corresponds to L2 Norm LNORM==3 corresponds to Cauchy LNORM==4 corresponds to SECH LNORM==5 corresponds to global correlation LNORM==6 corresponds to envelope objective function

Used LNORM=2 N_ORDER=0

--------------- No STF inversion -------------------

--------------- TRACE normalization -------------------

NORMALIZE=1

--------------- No Offset mute -------------------

--------------- Scale density update ------------------- SCALERHO = 0.500000

--------------- Scale Qs update ------------------- SCALEQS = 1.000000

---------------- No Joint Inversion ----------------

               TDFWI ITERATION 1     of 100

**Message from matcopy (written by PE 0): Copy material properties at inner boundaries ... finished (real time: 0.00 s). Vp_avg = 2486 Vs_avg = 1473 rho_avg = 2098 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/OBN_model_source_signal.0.su.shot1

* Starting simulation (forward model) for shot 1 of 16 **

Calculate residuals

PE 0 is writing 225 seismograms (vx) to su/OBN_RealData_x.su.shot1.it1 PE 0 is writing 225 seismograms (vy) to su/OBN_RealData_y.su.shot1.it1

* Starting simulation (adjoint wavefield) ** 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/OBN_model_source_signal.0.su.shot2

* Starting simulation (forward model) for shot 2 of 16 **

Calculate residuals

. . . . PE 2 outputs source time function in SU format to start/OBN_model_source_signal.2.su.shot16

* Starting simulation (forward model) for shot 16 of 16 **

Calculate residuals

* Starting simulation (adjoint wavefield) **

**Message from taper_grid (printed by PE 0): Coefficients for gradient taper are now calculated.

PE 0 is writing model to taper/taper_coeff_vert.bin.0.0

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to taper/taper_coeff_vert.bin Opening model files: taper/taper_coeff_vert.bin.??? ... finished. Copying... ... finished. Use ximage n1=1000 < taper/taper_coeff_vert.bin label1=Y label2=X title=taper/taper_coeff_vert.bin to visualize model.

**Message from taper_grid (printed by PE 0): Coefficients for gradient taper are now calculated. 21 25 490 500 4 10 5

PE 0 is writing model to taper/taper_coeff_hor.bin.0.0

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to taper/taper_coeff_hor.bin Opening model files: taper/taper_coeff_hor.bin.??? ... finished. Copying... ... finished. Use ximage n1=1000 < taper/taper_coeff_hor.bin label1=Y label2=X title=taper/taper_coeff_hor.bin to visualize model.

**Message from taper_grid (printed by PE 0): Coefficients for gradient taper are now calculated.

PE 0 is writing model to taper/taper_coeff_vert.bin.0.0

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

**Message from taper_grid (printed by PE 0): Coefficients for gradient taper are now calculated. 21 25 490 500 4 10 5

PE 0 is writing model to taper/taper_coeff_hor.bin.0.0

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

**Message from taper_grid (printed by PE 0): Coefficients for gradient taper are now calculated.

PE 0 is writing model to taper/taper_coeff_vert.bin.0.0

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

**Message from taper_grid (printed by PE 0): Coefficients for gradient taper are now calculated. 21 25 490 500 4 10 5

PE 0 is writing model to taper/taper_coeff_hor.bin.0.0

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to jacobian/OBN/jacobian_Test_p_vp.old Opening model files: jacobian/OBN/jacobian_Test_p_vp.old.??? ... finished. Copying... ... finished. Use ximage n1=1000 < jacobian/OBN/jacobian_Test_p_vp.old label1=Y label2=X title=jacobian/OBN/jacobian_Test_p_vp.old to visualize model.

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to jacobian/OBN/jacobian_Test_p.old Opening model files: jacobian/OBN/jacobian_Test_p.old.??? ... finished. Copying... ... finished. Use ximage n1=1000 < jacobian/OBN/jacobian_Test_p.old label1=Y label2=X title=jacobian/OBN/jacobian_Test_p.old to visualize model.

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to jacobian/OBN/jacobian_Test_c.old Opening model files: jacobian/OBN/jacobian_Test_c.old.??? ... finished. Copying... ... finished. Use ximage n1=1000 < jacobian/OBN/jacobian_Test_c.old label1=Y label2=X title=jacobian/OBN/jacobian_Test_c.old to visualize model.

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to jacobian/OBN/jacobian_Test_p_vs.old Opening model files: jacobian/OBN/jacobian_Test_p_vs.old.??? ... finished. Copying... ... finished. Use ximage n1=1000 < jacobian/OBN/jacobian_Test_p_vs.old label1=Y label2=X title=jacobian/OBN/jacobian_Test_p_vs.old to visualize model.

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to jacobian/OBN/jacobian_Test_p_u.old Opening model files: jacobian/OBN/jacobian_Test_p_u.old.??? ... finished. Copying... ... finished. Use ximage n1=1000 < jacobian/OBN/jacobian_Test_p_u.old label1=Y label2=X title=jacobian/OBN/jacobian_Test_p_u.old to visualize model.

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to jacobian/OBN/jacobian_Test_c_u.old Opening model files: jacobian/OBN/jacobian_Test_c_u.old.??? ... finished. Copying... ... finished. Use ximage n1=1000 < jacobian/OBN/jacobian_Test_c_u.old label1=Y label2=X title=jacobian/OBN/jacobian_Test_c_u.old to visualize model.

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to jacobian/OBN/jacobian_Test_p_mrho.old Opening model files: jacobian/OBN/jacobian_Test_p_mrho.old.??? ... finished. Copying... ... finished. Use ximage n1=1000 < jacobian/OBN/jacobian_Test_p_mrho.old label1=Y label2=X title=jacobian/OBN/jacobian_Test_p_mrho.old to visualize model.

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to jacobian/OBN/jacobian_Test_p_rho.old Opening model files: jacobian/OBN/jacobian_Test_p_rho.old.??? ... finished. Copying... ... finished. Use ximage n1=1000 < jacobian/OBN/jacobian_Test_p_rho.old label1=Y label2=X title=jacobian/OBN/jacobian_Test_p_rho.old to visualize model.

**Message from mergemod (printed by PE 0): PE 0 starts merge of 12 model files

writing merged model file to jacobian/OBN/jacobian_Test_c_rho.old Opening model files: jacobian/OBN/jacobian_Test_c_rho.old.??? ... finished. Copying... ... finished. Use ximage n1=1000 < jacobian/OBN/jacobian_Test_c_rho.old label1=Y label2=X title=jacobian/OBN/jacobian_Test_c_rho.old to visualize model. MYID = 0 pimaxr = 3.251180e+03 gradmaxr = 0.000000e+00 MYID = 0 umaxr = 2.038207e+03 gradmaxr_u = 0.000000e+00 MYID = 0 rhomaxr = 2.265330e+03 gradmaxr_rho = 0.000000e+00

**Message from matcopy (written by PE 0): Copy material properties at inner boundaries ... finished (real time: 0.00 s).

***** Starting simulation (test-forward model) no. 2 for shot 1 of 16 (rel. step length 0.02000000) Message from function wavelet written by PE 0 1 source positions located in subdomain of PE 0 have been assigned with a source signal. Message from PE 6 R U N - T I M E E R R O R: Simulation is unstable ! ...now exiting to system. Message from PE 7 R U N - T I M E E R R O R: Simulation is unstable ! ...now exiting to system. Message from PE 8 Abort(1) on node 6 (rank 6 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 6 Abort(1) on node 7 (rank 7 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 7



I have updated my files to include detailed settings. Could you please provide me with some advice based on this information?Once again, I would like to express my deepest gratitude for your valuable time, expertise, and assistance. Your guidance will be immensely valuable to overcome the challenges I’m facing. Thank you all for your kind support.
[DENISE_OBN_fIeldData_Flow.txt](https://github.com/daniel-koehn/DENISE-Black-Edition/files/12163225/DENISE_OBN_fIeldData_Flow.txt)
[FWI_workflow_marmousi.txt](https://github.com/daniel-koehn/DENISE-Black-Edition/files/12163226/FWI_workflow_marmousi.txt)

Best regards,

Zhang

daniel-koehn / DENISE-Black-Edition

Seeking Advice on FWI using DENISE with OBN Field Data: Challenges and Simulation Instability #46

----------------------- SEISMOGRAMS ---------------------- seismograms of x- and y-component of particle velocity. output-files: su/OBN_RealData_x.su su/OBN_RealData_y.su The data is written in IEEE SU-format . samplingrate of seismic data: 0.001000 s Number of samples per trace: 4000