Structural similarity regularization scheme for multiparameter seismic full waveform inversion
Maokun Li, Schlumberger, Lin Liang, Schlumberger, Aria Abubakar, Schlumberger, and Peter M. van den Berg, Delft University of Technology
SUMMARY
We introduce a new regularization scheme for multiparameter seismic full-waveform inversion (FWI). Using this scheme, we can constrain spatial variations of parameters which are having a weak sensitivity with the one that having a good sensitivity to the measurement, assuming that these parameters have similar-ities in their structures. In seismic FWI, we apply this scheme when inverting the P-wave velocity and mass density simulta-neously. Results from numerical tests show that this scheme may significantly improve the reconstruction of the mass den-sity. Since we obtain an improved mass-density distribution, the inverted P-wave velocity is also enhanced. Hence, we also obtain a better data fit. As numerical examples we show in-versions of both vertical seismic profiling (VSP) and surface seismic measurements.
INTRODUCTION
Seismic full waveform inversion (FWI) has attracted much at-tention recently in seismic data processing. It has been suc-cessfully used in the process of seismic velocity modeling to improve the resolution of velocity models (for examples, see Virieux et al. (2009); Vigh and Starr (2008); Plessix (2009); Routh et al. (2011)). Compared with the travel-time tomog-raphy approach, FWI optimizes velocity models in the sur-vey domain by matching the entire recorded seismic wave-form with the simulated wavewave-form, including both amplitude and time delay. Therefore, it is commonly agreed that we can extract more information from the waveforms by using FWI (for example, see Virieux and Operto (2009)). FWI algorithms usually include two parts: forward modeling and inversion al-gorithms. The forward modeling algorithms are usually based on acoustic or elastic wave equations. The solvers are imple-mented on computers using numerical methods such as finite difference, finite element, or integral equation. In FWI, we usually use nonlinear inversion algorithms to update the veloc-ity and densveloc-ity models iteratively, such as the nonlinear con-jugate gradient (NLCG) method, the quasi-Newton method, and the Gauss-Newton method. Because the sizes of seis-mic datasets and the sizes of velocity models can both be very large, computation of forward modeling and inversion in FWI require intensive computing resources.
As it is well-known, seismic FWI is an ill-posed problem, i.e. the solution can be non-existent, non-unique, and non-stable. One way to mitigate this problem is by adding a priori in-formation (constraint). This can be done through precondi-tioning (Operto et al. (2006)) or regularization in the inversion process. The application of regularization scheme to FWI has been studied in works by Hu et al. (2009); Ram´ırez and Lewis (2010); Guitton (2011); Asnaashari et al. (2013). Usually in acoustic and/or elastic FWI, the sensitivity of the mass density
to the measurement is weak, compared to the sensitivity of the P-wave velocity to the measurement data. Therefore, the mass density is often assumed to be constant, a known distribution in the inversion, or correlated to the P-wave velocity through an empirical formula such as Gardner.
In this proceeding paper, we focus on the acoustic FWI and we inverted the P-wave velocity and mass density simultaneously. The extension to the elastic FWI has been done and will be discussed in the presentation. Inverting these two parameters together can help to improve the reconstruction because the model is closer to the true physical model. However, their sen-sitivities to the measurement are not the same. Therefore, the inverted image of the mass density is usually noisy and lacks details. To improve the reconstruction of the mass density, a structural similarity constraint between P-wave velocity and mass density is a reasonable assumption. Sharp changes in the P-wave velocity also correlate to changes in the mass density. Therefore, we can improve the reconstruction of mass density by regularizing the mass density with a structural information derived from the P-wave velocity during the inversion process. In the framework of the multiplicative regularization (as intro-duced in van den Berg et al. (1999) and enhanced in van den Berg and Abubakar (2001)), this can be achieved by
redesign-ing the weight of the weighted L2-norm regularization function
for the mass density. The advantage of this type of regular-ization is that we avoid a bias in the inversion due to the as-sumption of the mass density distribution. Moreover it is more flexible because it allows a dynamic structural constraint from the P-wave velocity. This constraint is updated automatically when the P-wave velocity is updated in the inversion process. We tested this regularization scheme by inverting data from both vertical seismic profile (VSP) and surface seismic sur-veys. The numerical results show this scheme can further im-prove the structure of the reconstructed mass density. FORMULATION
In this work, we model the seismic wave propagation using the acoustic wave equation in the frequency domain as follows:
1 ρ ∂ ∂x ∂ p ∂x +1 ρ ∂ ∂z ∂ p ∂z +ω 2 λ p = − jωq, (1)
where p is pressure,ρ is mass density, ω is angular frequency,
q is source term, andλ is Lam´e modulus. We solve
Equa-tion 1 using the frequency-domain finite-difference (FDFD) approach with a fourth-order accuracy. This equation is dis-cretized using rectangular grids and the discrete pressure val-ues are defined at the center of each cell. We used the per-fectly matched layer (PML) of B´erenger (1994) as the bound-ary condition. Here, we consider solving the problem in 2D. After discretization, for each frequency, the linear system of
FWI using a structural similarity regularization scheme equations is solved using the multifrontal LU decomposition
method of Davis and Duff (1997). The details of our forward modeling algorithm can be found in Pan et al. (2012). In the inversion process, we iteratively minimize a multiplica-tive cost function as introduced by van den Berg et al. (1999) as follows
Φn(m) = Φd(m) Φmn(m), (2)
whereΦdis the normalized data misfit cost function measuring
the difference between the measured datad and the simulated
datas(m) generated from the model parameters m. The
func-tionΦm
n(m) is the regularization cost function at the n-th
itera-tion. In the inversion algorithm,m represents both normalized
P-wave velocity (VP) and normalized mass density (ρ) as the
following: m(r) =mVP(r),mρ(r) , (3) where mVP(r) = VP(r) VP,0(r), mρ(r) = ρ(r) ρ0(r), (4)
in which VP,0(r) and ρ0(r) are P-wave velocity and mass
den-sity of the initial model used in the inversion. The data misfit cost function can be written as
Φd(m) = XNF k=1 PNS j=1PNi=1R|Wi, j,k[di, j,k− si, j,k(m)]|2 2NFPNj=1S PNi=1R |Wi, j,kdi, j,k|2 = 1 2 kWd[d − s(m)]k2. (5)
Here, we assume that there are NSsources and NR receivers
in the survey. After converting the data into the frequency
do-main, we select NF frequencies for the inversion. di, j,kis the
measurement data at the i-th receiver, radiated from the j-th
source at the k-th frequency.s(m) represents the vector of the
simulated data, with its elements si, j,k. Wdis a diagonal data
weighting matrix with its elements Wi, j,krepresenting the data
weighting coefficient. For the regularization functionΦm
n(m),
we can use either the L2-norm or weighted L2-norm
regular-ization as described in van den Berg and Abubakar (2001). At the n-th iteration, it is given by
Φm n(m) = Z Db 2 n,VP(r) h |∇mVP(r)|2+δn2 i dr + Z Db 2 n,ρ(r)h ∇mρ(r) 2+δn2idr, (6)
where D is the inversion domain. The weighting factor bn,VP(r)
is given by b2n,VP(r) =Z 1 D h ∇mn,VP(r) 2+δn2 i dr (7)
for the L2-norm regularizer and
b2n,VP(r) = 1 V 1 ∇mn,VP(r) 2+δn2, V = Z Ddr (8)
for the weighted L2-norm regularizer. We can write similar
formulas for bn,ρ(r). The L2-norm regularizer tends to
fa-vor smooth profiles, while the weighted L2-norm regularizer is
known for its ability to preserve edges. The area of the
inver-sion domain is denoted by V . The parameterδ2
n is a non-zero
constant that is chosen to be equal to:
δ2
n =Φ d(mn)
∆x∆z , (9)
where∆xand∆zare the widths of the discretization cells in the
x and z directions. More details on this regularization function can also be found in Abubakar et al. (2008, 2009).
The structure similarity regularization scheme can be
imple-mented by employing the weighted L2-norm regularization
func-tion by weighting the variafunc-tion of the mass density with the variation of P-wave velocity as follows:
Φm n(m) = Z Db 2 n,VP(r) h |∇mVP(r)|2+δn2 i dr + Z Db 2 n,ρ,VP(r) h ∇mρ(r) 2+δn2idr, (10) where b2
n,VP(r) is defined as in equation 7 for the L2-norm
reg-ularization and as in equation 8 for the weighted L2-norm
reg-ularization. The weighting function b2
n,VP,ρ(r) is defined as b2n,ρ,VP(r) = 1 V 1 ∇mn,VP(r) 2+δn2, (11)
for both L2-norm and weighted L2-norm regularization.
In this formulation,Φm
n(m) is still dimensionless because both
mVP(r) and mρ(r) are normalized by the initial model and they
are dimensionless. In this regularization function the variation of the mass density in the n-th inversion iteration is constrained by the variation of the P-wave velocity of the (n − 1)-th inver-sion iteration. This enforces the mass density spatial variation to be similar as the one of the P-wave velocity. One thing to notice here is that this regularization only constrains the struc-ture profile in the mass density, but not the values of the mass density. The detailed derivation of the gradient and Hessian of this regularization function is similar as the one described in Abubakar et al. (2008).
In the FWI algorithm, we can employ the Gauss-Newton or preconditioned NLCG method to minimize the cost function in equation 2. For more details about our inversion algorithm implementation, please see Hu et al. (2009, 2011).
NUMERICAL EXAMPLES VSP data inversion
As the first example, we consider inverting a dataset from a synthetic VSP survey. The true model is shown in Figure 1. In this configuration we have a target with size of 200 × 100 m located 200 m away from the well and its depth is 450 m. In this VSP survey, we have 32 receivers, which are evenly distributed along the vertical well located at x = 0 m with depth ranging from 300 m to 600 m. There are 79 sources deployed on the surface (z = 0 m). Their horizontal distances from the well are within 100 m. Both sources and receivers are modeled
using monopoles. In the inversion, we invert data at 2, 3, and 5 Hz simultaneously. We use a perturbed background model as the initial model for both P-wave velocity and mass density as shown in Figure 2. This provides additional errors in the simulated data. We use the Gauss-Newton inversion algorithm
and the weighted L2-norm regularization in the inversion.
Figure 3 shows the inverted model using the weighted L2-norm
regularization for both P-wave velocity and mass density. We observe a good reconstruction of the P-wave velocity. Both location and size of the reconstructed target agree well with the true model. However, the mass density is not well recon-structed because its sensitivity to this dataset is weak. In order to improve the reconstruction of the mass density, we rerun the inversion again using the structural similarity regulariza-tion funcregulariza-tion. The inverted model is shown in Figure 4. We observe that the reconstruction of the mass density is signifi-cantly improved. Meanwhile, the reconstruction of the P-wave velocity is also slightly improved. The size and shape of the target are now nearly identical to the true model. This example clearly shows the value of the structural similarity regulariza-tion constraint. It can improve the reconstrucregulariza-tion by constrain-ing the weakly sensitive parameters usconstrain-ing the parameters with a reasonably good sensitivity to the measurement setup pro-vided that their structures have similarity.
(a) P-wave velocity (b) Mass density
Figure 1: True model
(a) P-wave velocity (b) Mass density
Figure 2: Initial model Surface seismic data inversion
In this example, we test the structural similarity regularization scheme using data generated from the well-known Marmousi model as shown in Figure 5. In the survey, we evenly deploy 58 monopole sources and 115 receivers on the surface of the inversion domain. The distances between neighboring sources
(a) P-wave velocity (b) Mass density
Figure 3: Inverted model using the weighted L2-norm
regular-ization
(a) P-wave velocity (b) Mass density
Figure 4: Inverted model using the weighted L2-norm
regular-ization and the structural similarity regularregular-ization
are 160 m and the distances between neighboring receivers are 80 m. The inversion domain is gridded into 460 × 150 cells with a cell size of 20 × 20 m. We apply multi-frequency se-quential inversion of the data at 1.5, 3, 6, 10, and 16 Hz. We add 5% random white noise to the data. We use the initial model with linearly increasing P-wave velocity and density values as shown in Figure 6. We employ the preconditioned nonlinear conjugate gradient method to invert the data. The maximum number of iterations is set to 100.
We first invert this data using the L2-norm regularization for
both the P-wave velocity and the mass density. The inverted model is shown in Figure 7. In the second step, we invert these data by applying the structural similarity regularization scheme to the mass density. Figure 8 shows the inverted model. Comparing these two results, we observe that the structural similarity constraint regularization may help to improve the structures in both P-wave velocity and mass density.
CONCLUSION
We developed a regularization scheme for enforcing the struc-tural similarities of P-wave velocity and mass density in the frequency-domain acoustic full-waveform inversion. This ap-proach assumes that the mass density usually has the same spa-tial variation as the P-wave velocity. By using this scheme, we can improve the reconstruction of mass density since it has a weaker sensitivity to the measurement than P-wave veloc-ity. This scheme can also be used in other multiparameter in-version algorithms, such as the elastic FWI and multi-physics
FWI using a structural similarity regularization scheme
(a) P-wave velocity
(b) Mass density
Figure 5: True model
(a) P-wave velocity
(b) Mass density
Figure 6: Initial model
(a) P-wave velocity
(b) Mass density
Figure 7: Inverted model using the L2-norm regularization
(a) P-wave velocity
(b) Mass density
Figure 8: Inverted model using the L2-normal regularization
and the structural similarity regularization
inversion. In the presentation we aim also to show some field data inversion results as well as the elastic multiparameter FWI results.
ACKNOWLEDGMENTS
The authors thank Dr. Guangdong Pan from Schlumberger-Doll Research for providing the 2D acoustic forward modeling and inversion code.
EDITED REFERENCES
Note: This reference list is a copy-edited version of the reference list submitted by the author. Reference lists for the 2013 SEG Technical Program Expanded Abstracts have been copy edited so that references provided with the online metadata for each paper will achieve a high degree of linking to cited sources that appear on the Web.
REFERENCES
