From modeling to full waveform inversion: A hands-on tour using Madagascar

Reconstruct your source/incident wavefield backwards in time

As one can see from above, the computation of FWI gradient requires simultaneously accessing the source/incident wavefield and the adjoint wavefield at each time step. To achieve this goal, we may store the absorbing boundary at each time step when doing forward simulation, and then reconstruct the incident wavefield backwards in time using the stored boundary. According to equation (2), the reconstruction is easy

$$p^{k-1}=2p^{k}-p^{k+1}+\Delta t^2 v^2 \nabla^2 p^k$$ (9)

Therefore, we may employ the same subroutine by switching the role of wavefield $p^{k+1}$ and $p^{k-1}$ . The elements in the boundary does not follow the above equation but we can store it in forward simulation and re-inject them for backward reconstruction.

1. Redo forward modeling using the same model while storing the boundary at each time step

2. reverse propagate the incident wavefield backwards from final snapshot and stored boundary: at each time step, inject the corresponding boundary

3. check whether you perfectly reconstruct your incident wavefield at any time step kt, as shown in Figure 3.

check
Figure 3. The backward reconstructed wavefield is the same as the incident wavefield

The final SConstruct looks like

from rsf.proj import *

Flow('vel0',None,
     	'''
     	math output=1.6 n1=50 n2=200 d1=0.005 d2=0.005
     	label1=x1 unit1=km label2=x2 unit2=km 
     	''')

Flow('vel1',None,
     	'''
     	math output=1.8 n1=50 n2=200 d1=0.005 d2=0.005
     	label1=x1 unit1=km label2=x2 unit2=km 
     	''')
Flow('vel2',None,
     	'''
     	math output=2.0 n1=100 n2=200 d1=0.005 d2=0.005
     	label1=x1 unit1=km label2=x2 unit2=km 
     	''')
Flow('vel',['vel0','vel1','vel2'], 'cat axis=1 ${SOURCES[1:3]}')

Flow('shot check1 check2','vel',
	'''	
	sffbrec2d check1=${TARGETS[1]} check2=${TARGETS[2]} 
	csdgather=n fm=15 dt=0.001 ns=1 ng=200 nt=1100 ng=200 kt=550	
	sxbeg=100 szbeg=2 jsx=37 jsz=0 
	gxbeg=0 gzbeg=1 jgx=1 jgz=0
	''')

Result('shot','grey gainpanel=all title=shot')

Flow('diff','check1 check2','sfadd ${SOURCES[1]} scale=1,-1')

Plot('check1','grey title="snapshot forward" scalebar=y')
Plot('check2','grey title="snapshot backward" scalebar=y')
Plot('diff','grey title="difference" scalebar=y')

Result('check','check1 check2 diff','SideBySideIso')

End()

From modeling to full waveform inversion: A hands-on tour using Madagascar