APPLIED GEOPHYSICS, Vol.14, No.4 (December 2017), P. 523–528, 9 Figures. DOI:10.1007/s11770-017-0629-6
Time-domain wavefield reconstruction inversion* Li Zhen-Chun1,2, Lin Yu-Zhao♦1,2, Zhang Kai1,2, Li Yuan-Yuan1,2, and Yu Zhen-Nan1 Abstract: Wavefield reconstruction inversion (WRI) is an improved full waveform inversion theory that has been proposed in recent years. WRI method expands the searching space by introducing the wave equation into the objective function and reconstructing the wavefield to update model parameters, thereby improving the computing efficiency and mitigating the influence of the local minimum. However, frequency-domain WRI is difficult to apply to real seismic data because of the high computational memory demand and requirement of time-frequency transformation with additional computational costs. In this paper, wavefield reconstruction inversion theory is extended into the time domain, the augmented wave equation of WRI is derived in the time domain, and the model gradient is modified according to the numerical test with anomalies. The examples of synthetic data illustrate the accuracy of time-domain WRI and the low dependency of WRI on low-frequency information. Keywords: Wavefield reconstruction, waveform inversion, augmented wave equation, timedomain inversion
Introduction The conventional method used for full waveform inversion (FWI) (Tarantola, 1984; Pratt, 1999) obtains the distribution of underground media parameters by conducting an inversion between synthetic and observed data. However, the objective function of FWI is a strong nonlinear function that is easily affected by local minima. The absence and inaccuracy of low-frequency information can lead to a mismatch between the travel time of the simulated record and that of real seismic data causing the inversion to fall into local minima. To solve this problem, Shin and Cha (2009) extended FWI to the Laplace domain to invert low-frequency data, and Moghaddam and Mulder (2012) proposed the use of
improved objective function inversion to reduce lowfrequency information requirements. To mitigate the influence of the local minimum on FWI, van Leeuwen and Herrmann (2013) proposed a new inversion method, known as wavefield reconstruction inversion (WRI), by expanding the target space. The method extends the solution space to the data and model spaces to increase the accuracy of the solution and to improve computational efficiency. In addition, van Leeuwen and Herrmann added the wave equation to the objective function in the frequency domain and calculated the model gradient by reconstructing the wavefield. This method did not require storing or calculating the forward and adjoint wavefields; therefore, the calculation efficiency is greatly improved. The WRI is stable and only weakly influenced by the
Manuscript received by the Editor May 19, 2016; revised manuscript received June 16, 2017. *This work was supported by the National Natural Science Foundation of China (Nos. 41374122 and 41504100). 1. School of Geosciences, China University of Petroleum, Qingdao 266580, China. 2. Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China. ♦Corresponding author: Lin Yu-Zhao (Email:
[email protected]) © 2017 The Editorial Department of APPLIED GEOPHYSICS. All rights reserved.
523
Time-domain wavefield reconstruction inversion local minimum when the solution space is expanded. van Leeuwen and Herrmann (2013) optimized the WRI theory and provided a physical meaning for the reconstructed wavefield along with a value for the balance factor in the objective function; subsequently Fang et al. (2015) proposed a seismic wavelet evaluation method for WRI. The WRI is usually conducted in the frequency domain; however, inversion in the time domain requires less computational memory. In addition, time-domain inversion theory is more instructive for practical applications. In this paper, the WRI is extended to the time domain, and the augmented wave equation and model gradient are deduced. An anomaly model is used to illustrate the stability of WRI, and Marmousi tests are used to demonstrate the feasibility of WRI. Finally, the advantages and limitations of this method are analyzed by the inversion of seismic parameters.
Frequency-domain WRI According to waveform-matching criteria, the fullwave inversion method can be used to obtain parameters of underground media using amplitude, travel time, and phase information. The traditional frequency domain FWI objective function can be expressed as (Pratt et al., 1998),
min I (m) m
M
1
¦2
2
Pi Ai (m)1 qi di , 2
i 1
(1)
where m denotes model parameters, d denotes the observed data, i is the frequency group, M is the total frequency, P is the shot operator (to ensure the shot location corresponds to the observed data), and A is the forward modeling operator A(m)=Ȧ2 diag (m)+L (where ω is frequency and L is the Laplace operator). The objective function of FWI is the functional under the L2 norm, and the function can be solved using the adjoint state method (Plessix, 2006) as follows:
mI
M
¦ Z diag (u ) v . 2 i
*
i
i
(2)
i 1
However, solving the above equation necessitates computing and storing the forward wavefields for all shot points at each frequency, Ai (m)ui qi and the adjoint wavefield v: Ai (m)* vi Pi * ( Pu i i d i ) . The inversion method therefore has high computational requirements and a low computational efficiency. Furthermore, in 524
the absence of an accurate initial model, the inversion results are susceptible to local minima. To solve this problem, van Leeuwen and Herrmann (2013) proposed an inversion method based on wavefield reconstruction by modifying the objective function, where the state equation is added to the objective function, and the wavefield and model are used as parameters to process the inversion. The general mathematical equation for the objective function in the frequency domain WRI is (van Leeuwen and Herrmann, 2013),
min IO (m, u ) m ,u
M
1 Pu ¦ i i di i 1 2
2 2
O2 2
2
Ai (m)ui qi 2 . (3)
On the right side of above equation, the first term is the misfit term between synthesized data and observed data in the classical FWI and the second term is the penalty term. The state equation is incorporated into the objective function to limit the model and wavefield, where λ is the balance factor. In deriving the wavefield by the objective function, the expression of wavefield reconstruction can be obtained by (van Leeuwen and Herrmann, 2013) [O 2 A(m)* A(m) Pi* Pi ]U i
O 2 A(m)* qi Pi* d ,
(4)
where * denotes the matrix complex conjugate and U is the reconstructed wavefield. Equation (4) is an augmented wave equation; the wavefield reconstructed by the equation satisfies the requirements of the state equation and contains information of the observed data. However, high computational memory is required to solve this complex equation. After reconstruction, the wavefield is no longer coupled with the model parameters, and the updated model gradient can be obtained directly from (van Leeuwen and Herrmann, 2013)
wJ wm
Zi2 diag (U i )* ( A(m)U i qi ).
(5)
When the gradient is solved, constructing the adjoint wavefield and storing the forward wavefield seems unnecessary. Only the residual of the wavefield is to reconstrut wavefield, and the calculation efficiency is greatly improved.
Time-domain WRI There is a relatively low demand for calculating
Li et al. memory in forward modeling or for inversion in the time domain, and the inversion method used in the time domain is more practical than in the frequency domain. The WRI is extended to the time domain, and the augmented wave equation and model gradient are stated in this section.
Time-domain augmented equation The objective function of WRI can be written as (van Leeuwen and Herrmann, 2013)
min IO (m, u ) m ,u
M
1 ui di ¦ i 1 2
2 2
O2 2
2
Ai (m)ui qi 2 . (6)
The wavefield reconstruction problem in the time domain can be approximated as a norm L2 minimization problem, and the wavefield update gradient can be obtained by deriving the wavefield from the objective function as follows:
gu
wJ wui
wA(m)ui (ui di ) O ( A(m)ui qi ). wui 2
(7a)
The second term of the wavefield gradient is 0, and therefore the reconstructed wavefield can be expressed as
Ui
ui D gu
ui D (ui d i ).
(7b)
Equation (7) is a time-domain augmented wave equation. When the same wave equation is used in forward and inversion methods during in numerical tests, the most accurate wavefield is obtained by setting the step length, α, to 1. The solution of augmented equation is limited by the state equation when dealing with real seismic data, and we can thus only use a part of the observed data that satisfies the state equation, which weakens the nonlinearity of inversion. Besides, the wavefield reconstruction method does not requires as much computational memory compared with that of the frequency domain. Moreover, the WRI dose not require the same storage memory for the forward modeling wavefield of the initial model compared with the classical FWI.
M
gm
¦ O 2 i 1
2 w 2U i [ A(m)U i qi ]. v3 wt 2
(8)
The gradient equation is found to exhibit a linear relation with the model perturbation, and the relation is proven using the following procedure. Replace the source term in the gradient expression by
qi
A(m)ui ,
A(m)U i qi
A(m)G u.
(9) (10)
The above equation decomposes the wave field residuals, δu, by forward modeling the operator. The wavefield residuals are caused by the difference between the real and initial models. Therefore, m G m u G u is taken into the wave equation, A(m)G u
2 w 2u G m. m3 wt 2
(11)
The above formula expresses the model perturbation and wavefield residual. In equation (11), the right side of the model, m, and the initial wavefield, u, can be treated as constant; thus A(m)G u v G m.
(12)
As shown, the model gradient expression has a linear relation with the velocity perturbation. The derivative of wavefield with time can be treated as a high-pass filter in the same way as those in FWI. The derivative term of wavefield versus time can be seen as a high-pass filter, which can be migrated when replaced as a normal wavefield term (Figure 4), and the part in brackets is linear to model perturbation. Furthermore, WRI updates the model based on the reconstructed wavefield and later determines the optimal solution from the data and model spaces, which broadens the search space of the solution and greatly reduces the influence of the local minimum.
Numerical examples
Time-domain WRI gradient After solving the augmented equation, the model parameters and wavefield exhibit a more linear behavior than classical FWI,, and the model gradient can be obtained by deriving the wavefield from the gradient expression as follows:
To demonstrate the accuracy of time-domain WRI, the proposed algorithm is tested using an anomaly model and a complex model extracted from the Marmousi model. 525
Time-domain wavefield reconstruction inversion
Anomaly model
the initial model (Figure 2). The grid interval for both horizontal and vertical directions is 10 m, and a Ricker wavelet with a main frequency of 10 Hz is applied to the model.
Figure 1 shows the true velocity model, which has a reflection layer and three anomalies. A homogeneous background model with one reflection layer is used as Distance (m) 1000
2000
3000
5000
1000
2700 2600
500
2500
1000
2400
Depth (m)
2300
1500
2200
2000
2100 2000
2500
1900
3000
1800
Distance (m) 3000
5000
1000
1000
150 100
1500
50
2200
1500
2100
2000
2000
2500
1900
3000
1800
high- and low-frequency formation (Figure 4a), and the energy of high-frequency data is much stronger than the low-frequency data, which causes an unstable inversion. To solve this problem, the filter term is replaced with a reconstructed wavefield, and the optimal gradient is shown in Figure 4b. The gradient contains a larger amount of low-wave number information, which makes the inversion more convergent. As illustrated in Figure 5, the WRI is capable of retrieving the position and energy of the anomalies, which proves the accuracy of this method.
2000
3000
Distance (m) 4000
5000
1000 250
1000
150 100
1500
50 0
2000
2500
-50
2500
-50
2500
-100
Distance (m)
Depth (m)
1000 1500 2000 2500 3000
2000
3000
5000 #104 7 6 5 4 3 2 1 0 -1 -2 -3
1000
100 50 0 -50 -100
Fig.3c Wavefield residual.
2000
3000
4000
5000 0.4
500
0.2 0
1000
Fig.4a Gradient calculated from equation (8).
526
150
Distance (m) 4000
Depth (m)
1000
5000
3000
Fig. 3b A reconstructed snapshot captured at 1200 ms.
Fig.3a A true snapshot captured at 1200 ms.
4000
1500
2000
3000
3000
1000
0
-100
2000
500
200
2000
3000
2400
1000
500
200
Depth (m)
Depth (m)
500
Velocity (m/s)
5000
2300
Distance (m) 4000
250
500
4000
Fig.2 An initial model.
Snapshots of the true model were captured in forward modeling at 1200 ms, and the wavefield was relatively well reconstructed. Figure 3a shows a wavefield snapshot of the true model, and Figure 3b shows the snapshot of reconstructed wavefield. To verify the accuracy of the augmented equation, the process of wavefield residual was performed and is shown in Figure 3c. Compared with the true snapshot, the anomalies in the reconstructed wavefield exhibit only a few differences, and the diferences can be neglected. As previously discussed, the model gradient includes both a 2000
3000
500
Fig.1 A true model.
1000
2000
Depth (m)
Depth (m)
Distance (m)
Velocity (m/s)
4000
-0.2
1500
-0.4
2000
-0.6
2500
-0.8
3000
-1
Fig.4b An optimal gradient.
Li et al. Distance (m) 1000
2000
3000
4000
Marmousi model
Velocity (m/s)
5000
In the following tests, we adopted the system same as that used in the anomalies model. An acoustic seismic modeling is performed with a high-order finitedifference method in the time domain. A section (Figure 6) measuring 5.0 × 3.0 km is adopted from the original Marmousi model. The smoothed model is taken as the initial model (Figure 7).
Depth (m)
2600
500
2500
1000
2400 2300
1500
2200 2100
2000
2000
2500
1900
3000
1800
Fig.5 The final result.
1000
Distance (m) 2000 3000
4000
5000
Velocity (m/s)
1000
3600
2800 2600 2400
2000
1000 Depth (m)
Depth (m)
2400
3000
1500
2200
1500 2000
2000
2200
1800
2000
2500
2500
1800
1600
1600
3000
3000
Fig.6 A Marmousi model.
Fig.7 The initial model.
propagated to the boundary, and the energy is extremely weak at the receiving end. Therefore, it is not possible to precisely reconstruct the wavefield, and it is difficult to describe the middle part of the model. However, WRI is less nonlinear, which enables a reliable velocity estimate to be obtained within 10 iterations (Figure 9).
When reconstructing the wavefield, strong direct waves or transmitted waves propagate as second sources, which cause extra energy clusters at the gradient boundary (upper left corner and upper right corner of Figure 8). However, the intermediate model is too complex, and therefore some waves cannot be
1000
Distance (m) 3000
2000
5000 Velocity (m/s)
500
3200
1000
4000
2600
3400
500
Distance (m) 2000 3000
Distance (m) 4000
5000 Velocity (m/s)
1000
2000
3000
4000
5000 Velocity (m/s)
0.2
500
0
1000
-0.2
1500
-0.4
2000
-0.6 -0.8
2500
3500
1000 Depth (m)
Depth (m)
500
3000
1500 2500
2000 2000
2500
-1
3000
3000
Fig.8 The first gradient.
Conclusions In this paper, WRI is extended to the time domain,
1500
Fig.9 The final result.
and a time-domain augmented wave equation and model gradient are deduced. The solution for the augmented wave equation in the time domain also contains limitations of observed data and the state equation, and 527
Time-domain wavefield reconstruction inversion it is thus possible to make complete use of real seismic data, which can be accurately predicted. After wavefield reconstruction, the wavefield and model parameters exhibit linear behavior, and the velocity field can be well retrieved. Reconstruction of the time-domain wavefield requires less computation memory and can be easily applied to actual production. Although the relevant theory for WRI is not yet mature at this stage, the theory discussed herein shows advantages with respect to improving the computing efficiency, calculating memory requirements, and providing high-resolution inversions. The results of this paper will be useful for studies of 3D inversion and multi-parameter inversion.
Acknowledgements We are grateful to SWPI laboratory staff for their support and providing use of Matlab software. We also wish to thank Zhang Hua and Yuan San-Yi for their constructive comments.
References De Hoop, A. T., 1960, A modification of Cagniard’s method for solving seismic pulse problems: Applied Scientific Research, Section B, 8(1), 349–356. Fang, Z., and Herrmann, F., 2015, Source estimation for wavefield reconstruction Inversion: 77th Annual International Conference and Exhibition, EAGE, Extended Abstracts, 1854–1857. Moghaddam, P. P., and Mulder, W. A., 2012, The diagonalator: inverse data space full waveform inversion: SEG Technical Program, Expanded Abstracts, 1–6. Mosegaard, K., and Tarantola, A., 1995, Monte Carlo sampling of solutions to inverse problems: Journal of Geophysical Research: Solid Earth, 100(B7), 12431– 12447. Mulder, W. A., and Hak, B., 2009, Simultaneous imaging of velocity and attenuation perturbations from seismic data is nearly impossible: 71th Conference & Technical Exhibition, EAGE, Extended Abstracts, S043. Peters, B., Herrmann, F. J., and Van, L. T., 2014, Waveequation Based Inversion with the Penalty MethodAdjoint-state Versus Wavefield-reconstruction Inversion: 76th Annual International Conference and Exhibition,
528
EAGE, Extended Abstracts, 3002–3005. Plessix, R. E., 2006, A review of the adjoint-state method for computing the gradient of a functional with geophysical applications: Geophysical Journal International, 167(2), 495–503. Pratt, R. G., 1999, Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model: Geophysics, 64(3), 888–901. Pratt, R. G., Shin C., and Hick, G. J., 1998, Gauss–Newton and full Newton methods in frequency–space seismic waveform inversion: Geophysical Journal International, 133(2), 341–362. Shen, P., and Symes, W. W., 2008, Automatic velocity analysis via shot profile migration: Geophysics, 73(5), VE49–VE59. Shin, C., and Cha, Y. H., 2009, Waveform inversion in the Laplace—Fourier domain: Geophysical Journal International, 177(3), 1067–1079. Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation: Geophysics, 49(8), 1259– 1266. Tarantola, A., and Valette, B., 1982, Generalized nonlinear inverse problems solved using the least squares criterion: Reviews of Geophysics, 20(2), 219–232. van Leeuwen, T., and Herrmann, F. J., 2013, Mitigating local minima in full-waveform inversion by expanding the search space: Geophysical Journal International, 195(1), 661–667. van Leeuwen, T., Herrmann, F. J., and Peters, B., 2010, A new take on FWI: Wavefield reconstruction inversion: 76th Annual International Conference and Exhibition, EAGE, Extended Abstracts, 2651–2654. van Leeuwen, T., and Mulder, W. A., 2010, A correlationbased misfit criterion for wave-equation traveltime tomography: Geophysical Journal International, 182(3), 1383–1394. Virieux, J., and Operto, S., 2009, An overview of fullwaveform inversion in exploration geophysics: Geophysics, 74(6), WCC1–WCC26.
Professor Li Zhen-Chun obtained an M.S. in Applied Geophysics in 1991 from the China University of Petroleum (East China) and a Ph.D. in Geophysics in 2000 from Tongji University. He currently works at the School of Geosciences of China University of Petroleum (East China). His main research interests are seismic modeling and imaging.
