Frequency-Wavenumber Domain Elastic Full Waveform Inversion with a Multistage Phase Correction

: Elastic full waveform inversion (EFWI) is essential for obtaining high-resolution multi-parameter models. However, the conventional EFWI may suffer from severe cycle skipping without the low-frequency components in elastic seismic data. To solve this problem, we propose a multistage phase correction-based elastic full waveform inversion method in the frequency-wavenumber domain, which we call PC-EFWI for short. Speciﬁcally, the seismic data are ﬁrst split using 2-D sliding windows; for each window, the seismic data are then transformed into the frequency-wavenumber domain for PC-EFWI misﬁt. In addition, we introduced a phase correction factor in the PC-EFWI misﬁt. In this way, it is possible to reduce phase differences between measured and synthetic data to mitigate cycle skipping by adjusting the phase correction factor in different scales. Numerical examples with the 2-D Marmousi model demonstrate that the frequency-wavenumber domain PC-EFWI with multistage strategy is an excellent way to reduce the risk of EFWI cycle skipping and build satisfactory start models for the conventional EFWI.


Introduction
The characterization of the velocity models is significant for improving seismic migration and geology interpretation [1][2][3][4][5]. Full waveform inversion (FWI) uses the kinematic and dynamic information of seismic waves and is considered an effective means for velocity model building [6][7][8][9][10]. However, the characteristics of the FWI misfit show strong nonlinearity [11][12][13][14][15][16][17]. Suppose the seismic data lack low-frequency components and the background velocity is very different from the accurate model. In that case, the FWI misfit may be trapped in local minima, leading to severe cycle skipping [18][19][20][21][22][23]. In addition, the physical properties of the subsurface are very complex, and the acoustic wave equationbased FWI cannot correctly obtain the multi-parameter models [11,[24][25][26][27][28]. Therefore, given the quality of existing seismic data, the proposed new strategy and method comprise an important way to solve cycle skipping for FWI.
The cycle skipping destroys the advantage of high-resolution acoustic FWI and elastic FWI, which has attracted wide attention [29][30][31][32][33]. In recent years, many new methods and strategies have been proposed to mitigate cycle skipping of FWI by building a better initial model. The multiscale strategy is commonly used for FWI to obtain high-resolution velocity models [8,[34][35][36][37]. However, the low-frequency seismic data are usually costly and contaminated by noise, leading to a multiscale strategy that usually does not work.

Window-Based Frequency-Wavenumber Domain Phase Correction
In this section, we try to utilize a window-based 2-D Fourier transform to develop a frequency-wavenumber domain phase correction-based elastic full waveform inversion (PC-EFWI). In this case, we first need to define a window-based 2-D Fourier transform, where F 2D [ · ] and F * 2D [ · ] mean the window-based 2-D Fourier forward and inverse transform, ω = 2π f , f means the frequency, k x is the wavenumber in x direction; u(t, x) are the seismic data, g(t, x) is a 2-D sliding window,ũ(τ, h, ω, k x ) are the window-based frequency-wavenumber domain seismic data; and τ and h are the displacements of t and x, respectively.
The seismic data can easily be separated into phase and amplitude components in the frequency-wavenumber domain as follows, whereũ,d are the frequency-wavenumber domain synthetic and measured data, respectively. Therefore, the phase differences betweenũ andd can be denoted as follows, where ψ u , ψ d are the phase information in the frequency-wavenumber domain. In this way, the phase correction can denoted as, where ε is a phase correction factor that controls the phase correction amount. If ε = 0.5, there are no phase differences between the corrected measured and synthetic data (Im ln ). If ε = 0, it hasũ c =ũ,d c =d. Therefore, it is possible to adjust the phase differences between the corrected measured and synthetic data using the phase correction factor. The horizontal component of the measured elastic seismic data and the 2-D window are shown in Figure 1. In this case, we first use sliding windows to split seismic data and obtain Figure 2a,d. Equation (1) is then used to obtain the frequency-wavenumber domain seismic data. After that, we set ε = 0.5 with Equation (4), and we can obtain the phasecorrected seismic data. In addition, we also show the phase-corrected seismic waveforms in Figure 3. Comparing the original waveforms and phase-corrected waveforms demonstrates that the frequency-wavenumber domain phase correction strategy can effectively reduce the phase differences between the measured and synthetic data.
are the frequency-wavenumber domain synthetic and measured data, respectively. Therefore, the phase differences between ũ and d can be denoted as follows, where d u   , are the phase information in the frequency-wavenumber domain. In this way, the phase correction can denoted as, where  is a phase correction factor that controls the phase correction amount. If , there are no phase differences between the corrected measured and synthetic data ( Therefore, it is possible to adjust the phase differences between the corrected measured and synthetic data using the phase correction factor. The horizontal component of the measured elastic seismic data and the 2-D window are shown in Figure 1. In this case, we first use sliding windows to split seismic data and obtain Figure 2a (1) is then used to obtain the frequency-wavenumber domain seismic data. After that, we set 5 . 0 =  with Equation (4), and we can obtain the phase-corrected seismic data. In addition, we also show the phase-corrected seismic waveforms in Figure 3. Comparing the original waveforms and phase-corrected waveforms demonstrates that the frequency-wavenumber domain phase correction strategy can effectively reduce the phase differences between the measured and synthetic data.

Review of Elastic Full Waveform Inversion
The elastic wave equation is the basis of EFWI. In isotropic media, the second order partial differential-based 2D elastic wave equations can be written as, where λ and µ are Lamé coefficients; ρ is density, which is taken as a constant value (ρ = 2 g/cm 3 ) in the tests of EFWI. In Equation (5), the u = [u x , u z ] and σ = [σ xx , σ zz , σ xz ] are displacement components and stress components, respectively. The misfit of EFWI is to measure the elastic seismic data residuals between the measured and synthetic data as follows [8,11]: where ns, nr indicate the number of shots and receivers, respectively; u, d are the timedomain measured and synthetic elastic seismic data, which contain horizontal and vertical The EFWI misfit can be minimized by updating the elastic parameters with an optimization algorithm. The gradients of the elastic parameters can be expressed as follows, where ← u z are the displacements of the forward and back propagated stress components, respectively. In addition, Lamé coefficients are related to the P and S velocities (V p , V s ), which can be denoted as follows, Therefore, the elastic parameters can be denoted by the P and S velocities m = v p , v s . The gradient for V p and V s can be expressed as follows,

Frequency-Wavenumber Domain Phase Correction-Based EFWI
The nonlinearity of the conventional EFWI misfit makes the inverted results easily fall into local minima. Therefore, a frequency-wavenumber domain phase correctionbased EFWI (PC-EFWI) is proposed to solve this problem, which helps to reduce the phase differences between measured and synthetic data. The proposed PC-EFWI misfit is as follows, Phase correction factor ε ∈ [0, 0.5] is used to control the amount of phase correction. In Equation (10), if we set ε = 0.5, the misfit only matches the amplitude information, which becomes envelope inversion in the frequencywavenumber domain [38]. In this way, we can obtain a better initial model with ε = 0.5 for conventional EFWI. Then, we can gradually reduce the phase correction factor to match the phase and amplitude differences at the same time to obtain detailed structures. The partial derivative of the PC-EFWI misfit is as follows, where * means the complex conjugation. Taking a further derivative, we have Remote Sens. 2022, 14, 5916 The gradient of the PC-EFWI misfit is as follows, According to Equation (1), the final gradient of the PC-EFWI misfit is as follows, Additionally, the adjoint source is as follows, where β is a small positive number to avoid dividing over zero. Therefore, the PC-EFWI gradients can be obtained by correlating forward and backward wave-fields.

Numerical Test
In this section, we apply complex modified Marmousi models to test the frequencywavenumber domain phase correction-based elastic full waveform inversion (PC-EFWI). The Marmousi model consists of various complex geological structures, such as anticlines, faults, angular unconformities, and oil reservoirs. Figure 4 shows that the oil reservoir is deep in the Marmousi model and overlaid by complex geological structures. However, the conventional EFWI is challenging to invert the deep oil reservoir, especially when seismic data lack low-frequency components. In this case, the Marmousi velocity models are 1.2 km × 3.6 km, and their initial models are shown in Figure 5. The acquisition system consists of 40 seismic sources and 240 receivers on the top of the velocity model. This article uses a low−cut 5 Hz Ricker wavelet with a peak frequency of 12 Hz to test the low-frequency dependence of the PC-EFWI method. In addition, we first use low−frequency data (5-10 Hz) to recover macro structures. Then, we use it as the initial model for high−frequency data (5-15 Hz) conventional EFWI.
The conventional EFWI results are shown in Figure 6. Compared to the accurate Marmousi models, the conventional EFWI results produce a cluster of anomaly values and artifacts on the left side of the Marmousi velocity models. This is mainly because the conventional EFWI suffers from severe cycle skipping. Then, we use high-frequency seismic data to invert the detailed geological structures with Figure 6 as the initial model. According to Figures 6 and 7, if the low-frequency EFWI produces cycle skipping, the high-frequency data-based EFWI will also produce errors. Therefore, it is essential to use new methods to recover the low-wavenumber components of the velocity models. especially when seismic data lack low-frequency components. In this case, the Marmousi velocity models are km km 6 . 3 2 .
1  , and their initial models are shown in Figure 5. The acquisition system consists of 40 seismic sources and 240 receivers on the top of the velocity model. This article uses a low−cut 5 Hz Ricker wavelet with a peak frequency of 12 Hz to test the low-frequency dependence of the PC-EFWI method. In addition, we first use low−frequency data (5-10 Hz) to recover macro structures. Then, we use it as the initial model for high−frequency data (5-15 Hz) conventional EFWI.   especially when seismic data lack low-frequency components. In this case, the Marmousi velocity models are km km 6 . 3 2 .
1  , and their initial models are shown in Figure 5. The acquisition system consists of 40 seismic sources and 240 receivers on the top of the velocity model. This article uses a low−cut 5 Hz Ricker wavelet with a peak frequency of 12 Hz to test the low-frequency dependence of the PC-EFWI method. In addition, we first use low−frequency data (5-10 Hz) to recover macro structures. Then, we use it as the initial model for high−frequency data (5-15 Hz) conventional EFWI.  The conventional EFWI results are shown in Figure 6. Compared to the accurate Marmousi models, the conventional EFWI results produce a cluster of anomaly values and artifacts on the left side of the Marmousi velocity models. This is mainly because the conventional EFWI suffers from severe cycle skipping. Then, we use high-frequency seismic data to invert the detailed geological structures with Figure 6 as the initial model. According to Figures 6 and 7, if the low-frequency EFWI produces cycle skipping, the high-frequency data-based EFWI will also produce errors. Therefore, it is essential to use new methods to recover the low-wavenumber components of the velocity models.  Figure 8. Compared to Figure 6, the geological macro structures of the Marmousi model inverted by PC-EFWI (Figure 8) seem better than the conventional EFWI result, because when the phase correction factor 5 . 0 =  , the PC-EFWI misfit becomes a pure amplitude misfit, similar to the envelope misfit with a solid ability to recover macro structures. In addition, we can gradually reduce the phase correction factor to match the phase differences between the measured and synthetic data. For example, we then set 3 . 0 =  and use Figure 8 as the initial velocity model to obtain detailed geological structures, shown in Figure 9. After that, we can use the inverted background Marmousi model (Figure 9) as the initial model for the conventional EFWI to invert detailed geological structures (Figures 10 and 11). A comparison of Figure 11 and Figure 7 shows that the PC-EFWI can better recover the oil reservoir, and the deep geological structures can be seen clearly. The conventional EFWI results are shown in Figure 6. Compared to the accurate Marmousi models, the conventional EFWI results produce a cluster of anomaly values and artifacts on the left side of the Marmousi velocity models. This is mainly because the conventional EFWI suffers from severe cycle skipping. Then, we use high-frequency seismic data to invert the detailed geological structures with Figure 6 as the initial model. According to Figures 6 and 7, if the low-frequency EFWI produces cycle skipping, the high-frequency data-based EFWI will also produce errors. Therefore, it is essential to use new methods to recover the low-wavenumber components of the velocity models.  Figure 8. Compared to Figure 6, the geological macro structures of the Marmousi model inverted by PC-EFWI (Figure 8) seem better than the conventional EFWI result, because when the phase correction factor 5 . 0 =  , the PC-EFWI misfit becomes a pure amplitude misfit, similar to the envelope misfit with a solid ability to recover macro structures. In addition, we can gradually reduce the phase correction factor to match the phase differences between the measured and synthetic data. For example, we then set 3 . 0 =  and use Figure 8 as the initial velocity model to obtain detailed geological structures, shown in Figure 9. After that, we can use the inverted background Marmousi model (Figure 9) as the initial model for the conventional EFWI to invert detailed geological structures (Figures 10 and 11). A comparison of Figure 11 and Figure 7 shows that the PC-EFWI can better recover the oil reservoir, and the deep geo- Now, the low-frequency seismic data-based PC-EFWI results with phase correction factor ε = 0.5 are shown in Figure 8. Compared to Figure 6, the geological macro structures of the Marmousi model inverted by PC-EFWI (Figure 8) seem better than the conventional EFWI result, because when the phase correction factor ε = 0.5, the PC-EFWI misfit becomes a pure amplitude misfit, similar to the envelope misfit with a solid ability to recover macro structures. In addition, we can gradually reduce the phase correction factor to match the phase differences between the measured and synthetic data. For example, we then set ε = 0.3 and use Figure 8 as the initial velocity model to obtain detailed geological structures, shown in Figure 9. After that, we can use the inverted background Marmousi model (Figure 9) as the initial model for the conventional EFWI to invert detailed geological structures (Figures 10 and 11). A comparison of Figures 7 and 11 shows that the PC-EFWI can better recover the oil reservoir, and the deep geological structures can be seen clearly. (a) (b) Figure 11. High-frequency data (5-10 Hz) based EFWI results using Figure 10 as the initial velocity model; (a) Vp; (b) Vs.
The inverted velocity profiles with conventional EFWI and PC-EFWI are shown in Figure 12. Figure 12a,b show that the conventional EFWI has severe cycle skipping, and the final inverted models are deplorable. The comparison results in Figure 12a-d show that the Vp and Vs of PC-EFWI + EFWI results are more similar to the accurate elastic velocity of Vp and Vs. The velocity comparison in Figure 12 shows that the frequen- The inverted velocity profiles with conventional EFWI and PC-EFWI are shown in Figure 12. Figure 12a,b show that the conventional EFWI has severe cycle skipping, and the final inverted models are deplorable. The comparison results in Figure 12a-d show that the Vp and Vs of PC-EFWI + EFWI results are more similar to the accurate elastic velocity of Vp and Vs. The velocity comparison in Figure 12 shows that the frequen- The inverted velocity profiles with conventional EFWI and PC-EFWI are shown in Figure 12. Figure 12a,b show that the conventional EFWI has severe cycle skipping, and the final inverted models are deplorable. The comparison results in Figure 12a-d show that the Vp and Vs of PC-EFWI + EFWI results are more similar to the accurate elastic velocity of Vp and Vs. The velocity comparison in Figure 12 shows that the frequen- (a) (b) Figure 11. High-frequency data (5-10 Hz) based EFWI results using Figure 10 as the initial velocity model; (a) Vp; (b) Vs.
The inverted velocity profiles with conventional EFWI and PC-EFWI are shown in Figure 12. Figure 12a,b show that the conventional EFWI has severe cycle skipping, and the final inverted models are deplorable. The comparison results in Figure 12a-d show that the Vp and Vs of PC-EFWI + EFWI results are more similar to the accurate elastic velocity of Vp and Vs. The velocity comparison in Figure 12 shows that the frequen- The inverted velocity profiles with conventional EFWI and PC-EFWI are shown in Figure 12. Figure 12a,b show that the conventional EFWI has severe cycle skipping, and the final inverted models are deplorable. The comparison results in Figure 12a-d show that the V p and V s of PC-EFWI + EFWI results are more similar to the accurate elastic velocity of V p and V s . The velocity comparison in Figure 12 shows that the frequency-wavenumber domain PC-EFWI misfit can successfully build better initial models for EFWI and mitigate the cycle skipping. cy-wavenumber domain PC-EFWI misfit can successfully build better initial models for EFWI and mitigate the cycle skipping.

Discussion
In this article, the 2-D window size for PC-EFWI is an important parameter, which is similar to the window size for Short-time Fourier transform (STFT). In STFT, the time and frequency resolution is a trade-off problem. Although a narrow-width window results in a better resolution in the time domain, it generates a poor resolution in the frequency domain, and vice versa. Similarly, although a narrow-width 2-D window results in a better exhibition of local-scale characteristics in the time-space domain, it can not encompass a complete waveform of seismic events, and vice versa. Therefore, the choice

Discussion
In this article, the 2-D window size for PC-EFWI is an important parameter, which is similar to the window size for Short-time Fourier transform (STFT). In STFT, the time and frequency resolution is a trade-off problem. Although a narrow-width window results in a better resolution in the time domain, it generates a poor resolution in the frequency domain, and vice versa. Similarly, although a narrow-width 2-D window results in a better exhibition of local-scale characteristics in the time-space domain, it can not encompass a complete waveform of seismic events, and vice versa. Therefore, the choice of 2-D window size involves a certain amount of experience. In this article, the window length in the time direction is 0.36 s, and the window length in the distance direction is 0.6 km (Figure 2). In addition, the PC-EFWI misfit needs to calculate the forward and inverse window based 2-D Fourier transform both of synthetic and measured data, which requires a substantial amount of computing resources. Fortunately, a large sliding step length can greatly improve calculation efficiency. Additionally, the sliding step length of the 2D window is 1/10 of the window length in the time and distance directions. In this way, an appropriate sliding step length can not only ensure the accuracy of the window based 2-D Fourier transform, but also improve the calculation efficiency.

Conclusions
We proposed a frequency-wavenumber domain PC-EFWI to invert for the lowwavenumber components of subsurface structures. With a phase correction factor adopted in the PC-EFWI misfit, we can flexibly adjust the phase differences between the measured and synthetic data to reduce the nonlinearity of the EFWI misfit, which mitigates the cycle skipping problem to a large extent. For the multistage strategy, we first use a large phase correction factor to eliminate phase differences in the PC-EFWI misfit, which is suitable to invert the macro structures of the complex Marmousi models. Then, we can gradually reduce the phase correction factor to form a more robust multistage inversion strategy to obtain detailed geological structures. Numerical tests show that the frequency-wavenumber domain PC-EFWI method can successfully build accurate background velocity models and mitigate cycle skipping of the EFWI.