Full-waveform inversion using seislet regularization

Because of inaccurate, incomplete, and inconsistent waveform records, full-waveform inversion (FWI) in the framework of a local optimization approach may not have a unique solution, and thus it remains an ill-posed inverse problem. To improve the robustness of FWI, we have developed a new model regularization approach that enforced the sparsity of solutions in the seislet domain. The construction of seislet basis functions requires structural information that can be estimated iteratively frommigration images.We implement FWI with seislet regularization using nonlinear shaping regularization and impose sparseness by applying soft thresholding on the updated model in the seislet domain at each iteration of the data-fitting process. The main extra computational cost of the method relative to standard FWI is the cost of applying forward and inverse seislet transforms at each iteration. This cost is almost negligible compared with the cost of solving wave equations. Numerical tests using the synthetic Marmousi model demonstrate that seislet regularization can greatly improve the robustness of FWI by recovering high-resolution velocity models, particularly in the presence of strong crosstalk artifacts from simultaneous sources or strong random noise in the data.


INTRODUCTION
Full-waveform inversion (FWI) is a data-fitting procedure used to construct high-resolution quantitative subsurface models by exploiting full information in the observed data (Lailly, 1983;Tarantola, 1984;Virieux and Operto, 2009). However, the problem is inherently ill-posed and suffers from artifacts that could be falsely interpreted as "geological structures." One way to mitigate this problem is by adding regularization or preconditioning to guide the inversion toward a model consistent with a priori geologic or geophysical constraints. The choice among different regularization techniques depends on the specific problem, and it can be governed by the need to preserve or emphasize particular desired features of the model (Loris et al., 2007). For example, if one is mostly interested in large-scale features, it is natural to introduce a regularizing constraint or penalty term involving spatial derivatives. On the other hand, if one seeks to find solutions that are sparse with respect to a given basis, this can be achieved by imposing a sparsity constraint involving appropriate transforms, such as Fourier, wavelet, or curvelet transforms (Loris et al., 2007). This regularization procedure finds solutions that are faithfully represented by a relatively small number of nonzero coefficients in the transformed domain.
We propose to impose sparsity regularization on the model in the seislet domain (Fomel and Liu, 2010) to improve the robustness of FWI. We refer to this regularization as seislet regularization. A model in the seislet domain is expressed using basis functions aligned along locally planar structures. Other researchers have previously used the sparsity of velocity model in other transform domains for velocity model building. For instance, Loris et al. (2007) apply l 1 -norm regularization in a wavelet basis to solve global seismic tomography problems, allowing for the possibility of sharp discontinuities superimposed on a smoothly varying background model. Li et al. (2012) compute the model updates from random subsets of data and use curvelet-domain sparsity promotion to suppress crosstalk between different sources. Curvelets are appropriate for seismic data because they provide a provably optimal decomposition of wave-propagation operators (Candès and Demanet, 2005). In this paper, sparsity regularization for a velocity model is implemented by imposing a soft thresholding on the updated model in the seislet domain. Compared with other transform domains, seislets exhibit superior data compression and sparsity for events with dominant local slopes . The classic digital wavelet transform is equivalent to the seislet transform with an erroneous zero slope (Fomel and Liu, 2010).
Seislet regularization allows us to build a model that fits the data and has a strong tendency to be sparse in the seislet basis. Because seislet basis functions are aligned along locally planar structures, this helps to attenuate random noise and build geologically consistent models. Different approaches in generating geologically sensible models for seismic inversion have been proposed before. Guitton et al. (2012) use a directional Laplacian filter as a model preconditioning operator in FWI to smooth gradients along geologic dips. Ma et al. (2012) propose to invert for a sparse velocity model in FWI, and connect the sparse and dense models through image-guided interpolation. Xue et al. (2016) incorporate linear shaping regularization (Fomel, 2007) into least-squares reverse time migration (RTM) and use structure-enhancing filtering to mitigate artifacts caused by simultaneous-source or incomplete data. Compared with these methods, the proposed method offers an alternative formulation with a highly efficient implementation. It is formulated as the sparsity constraint on the model in the seislet domain, and it is implemented by imposing a soft thresholding on the updated model at each iteration.
In this paper, we first introduce FWI with and without seislet regularization, and illustrate their implementation differences. Then, we use two numerical examples to verify the effectiveness of the proposed method in improving the robustness of FWI by suppressing artifacts caused by encoded data and random noise. Both examples are based on the 2D Marmousi synthetic model (Bourgeois et al., 1991).

THEORY
The objective function of standard FWI can be written as where m is the velocity model, d obs is the observed waveform data, and F stands for the nonlinear forward-modeling operator. A popular local optimization algorithm for minimizing this function is based on the nonlinear conjugate-gradient method (Mora, 1987;Tarantola, 1987), where the model is updated at the iteration n in the direction of s n , which is a linear combination of the gradient at iteration n, ∇Jðm n Þ, and the previous conjugate direction s n−1 . This optimization scheme can be formulated as follows: where scalar α n can be obtained by performing a line search and β n is designed to guarantee that s n and s n−1 are conjugate.
We propose to apply nonlinear shaping regularization (Fomel, 2008) to impose sparseness of the velocity model in the seislet domain and modify the iteration process in equation 2 as follows: (3) where S −1 and S stand for the inverse and forward seislet transforms, respectively. Here, T is the soft thresholding operator and has the following function expression: (4) where t is a positive threshold level. If the seislet transform were an unitary operator, the iteration process in equation 3 would be equivalent to using the linearized Bregman iteration (Yin et al., 2008;Cai et al., 2009) to minimize the following objective function: where λ is a regularization parameter. The difference of equation 3 to the iterative scheme described in equation 2 is that the new scheme involves a model shaping operator, which is the combination of the inverse seislet transform, soft thresholding operator, and forward seislet transform. The purpose of the shaping operator S −1 TS is to remove possible artifacts existing in the updated model and to gradually achieve a geologically meaningful model. In implementation of equation 3, step length α n is estimated by line search without including the shaping operator, and we apply the shaping operator to the updated model after the line search to get the new model m nþ1 . In this way, we just use the shaping operator once at each iteration.
The seislet transform requires the input of local slope because in the seislet domain, the model is represented by the basis functions that are aligned with local structures. Figure 1 shows a simple test to check the properties of the velocity model in the seislet domain. We compare seislet coefficients and basis functions with those of the digital wavelet transform, using the Marmousi model (Figure 1a). The dip required by the seislet transform is directly estimated from the Marmousi model in this test. Figure 1b and 1c shows the Marmousi model in the wavelet and seislet domains, respectively. We observe that the nonzero seislet coefficients focus at coarse scales, but the wavelet transform has small residual coefficients at fine scales. From this observation, we can conclude that the seislet transform has better data compression and sparsity than the wavelet transform, which is further verified in Figure 1d, in which a significantly faster decay of the seislet coefficients is evident. Figure 1e and 1f shows some randomly selected representative basis functions for the wavelet and seislet transforms, respectively. We can observe that the direction of the seislet basis functions follows structural direction, whereas the wavelet basis functions correspond to the zero slope.

Experiment setup
We use the classic 2D Marmousi model (Bourgeois et al., 1991) to test the effectiveness of the proposed method. The model was first subsampled from 2301 × 751 cells at 4 m cell size to 512 × 188 cells at 16 m cell size, as shown in Figure 1a. Figure 2a shows the initial model for FWI, obtained by smoothing the target model.

Xue et al.
A fixed-spread surface acquisition is used, involving 32 sources located at every 240 m. A Ricker wavelet with center frequency of 13 Hz is used to generate the synthetic data. Figure 2b shows a shot record with source location at 4 km. We perform frequency-domain FWI using the fast Helmholtz solver (Li et al., 2014). We used eight frequencies ranging from 4 to 11 Hz and carried out 10 iterations for each frequency. We perform two tests: one with encoded data and the other with noisy data.

Dip estimation
Before we conduct FWI with seislet regularization, we need to provide local-dip information for the seislet transform. We estimate the local dip from the RTM image using plane-wave destruction (Fomel, 2002). The dip estimation is computationally inexpensive compared with RTM. With the initial velocity, we can get an initial RTM image as shown in Figure 3a. Although some structures of the image appear at wrong locations, we can get an acceptable local dip map (Figure 3b). We use this initially estimated dip to perform the first 30 iterations (for the first three frequencies). After we finish the iterations at 6 Hz, a more accurate velocity is obtained, with which we can get a better RTM image. Then, we can estimate the second dip for the next 30 iterations (for the frequencies 7-9 Hz). Similarly, we can perform the third RTM (Figure 3c) after the inversions with 9 Hz frequency and get a much better dip (Figure 3d), which is used for the last 20 iterations (for the frequencies 10-11 Hz). Comparing Figure 3a and 3c, we can find that the previous position of structures in the central and deep parts has been corrected in the new RTM image. Note that we use unencoded data to perform RTM in the encoded data test to avoid crosstalk artifacts in the RTM images. In the noisy data test, the assumed noisy data are used to carry out RTM.

Encoded data test
To reduce the computational burden, one well-known technique is to build supershots by assembling several sources, which can re-duce the computational cost by a factor roughly equal to the number of sources assembled together (Morton and Ober, 1998;Romero et al., 2000). However, the source-combination technique may introduce crosstalk noise arising from the interference of waves generated at spatially adjacent sources. To remove the unwanted crosstalk noise, the blended data are typically paired with phaseencoding strategies; among them, dynamic phase encoding (Krebs et al., 2009;Baumstein et al., 2011;Ben-Hadj-Ali et al., 2011) is a particularly effective approach. Dynamic phase encoding in frequency-domain FWI involves changing the encoding code and building a new encoded superset at each iteration of each frequency inversion (Ben-Hadj-Ali et al., 2011).
In this test, we combine every eight equidistant shots in one supershot, creating four supershots in total. The first supershot contains shots with the following indices: 1, 5, 9, 13, 17, 21, 25, and 29. We perform three inversions. The first inversion is the standard FWI with the blended data without phase encoding, which means that data are directly blended together without any designed codes. The second inversion is the standard FWI with the dynamic phase encoding technique, and at each iteration, a new encoding code is generated to build a new super data set. The encoding code used in this study is the simple phase function e iγ , where γ is a random number in ½0; 2π. Figure 4 shows the real part of the four encoded supershots at the first iteration of the 6 Hz frequency data inversion. From each supershot, we can roughly observe eight large peaks and troughs, and their locations correspond to that of the original sources. The third inversion inverts the dynamic phase encoded data by using FWI with seislet regularization. We set the soft thresholda) b) c) d) Figure 3. Dip estimation from migration images: (a) RTM with initial velocity, (b) local dip map estimated from (a), which will be used for the first 30 iterations, (c) third RTM with updated velocity after 60 iterations (4-9 Hz data), and (d) dip map estimated from (c) and used for the last 20 iterations (10-11 Hz data).
ing parameter in the seislet regularization empirically to be 18%, meaning that 82% of the smaller seislet coefficients get removed at each iteration. Figure 5 shows the final results of the three inversions. All the results were obtained after 80 iterations. The result of standard FWI without phase encoding contains visible crosstalk artifacts, and the model is blurred by noise. Dynamic phase encoding can effectively suppress some of the artifacts, but there are some remaining artifacts. As shown in Figure 5c, seislet regularization leads to a noise-free and high-resolution model. Because this is a synthetic data inversion test, we can also display the evolution of the model misfits. As shown in Figure 5d, FWI with seislet regularization has a faster model convergence rate.

Noisy data test
To further test the robustness of the proposed method, we generate a noisy data set by adding strong random noise to the original time domain data. Then, we transform the noisy data to the frequency domain for inversion. We perform two inversions: standard FWI and FWI using seislet regularization, and we compare their results in Figure 6a and 6b. In this comparison, we find that seislet regularization suppresses the noise caused by the ambient noise in the data, and it helps to get a good inversion result. Figure 6c shows the model convergence curves, which also tell us that FWI using seislet regularization has a faster model convergence rate. Finally, we show the data convergence of FWI using seislet regularization for each frequency inversion in Figure 6d. We can observe that at each separate frequency, the proposed method exhibits a fast datafitting convergence.

CONCLUSION
With the observation that velocity models with planar structures appear sparse in the seislet domain, we introduce a sparsity regularization scheme with seislet transform to improve the robustness of FWI and to recover a noise-free, geologically consistent velocity model. We implement this seislet regularization by imposing soft thresholding on the updated model in the seislet domain at each iteration. We estimate the dip needed by the seislet transform from the migration image generated during the iterative inversion. Two numerical tests verify the effectiveness of FWI using seislet regularization.