c
TU Delft, Delft The Netherlands, 2006
AN HELMHOLTZ ITERATIVE SOLVER FOR THE
THREE-DIMENSIONAL SEISMIC IMAGING PROBLEMS ?
Ren´e-Edouard Plessix∗
∗Shell International E&P,
Kesslerpark 1, 2288 GD, Rijswijk, The Netherlands e-mail:
web page: http://www.shell.com Key words: iterative solver, wave equation, Helmholtz, seismic
Abstract. A preconditioned iterative solver for the three-dimensional frequency-domain wave equation applied to seismic problem is evaluated. The preconditioner corresponds to the approximate inverse of a heavily damped wave equation deduced from the (undamped) wave equation. The approximate inverse is computed with one multigrid cycle. Numerical results show that the method is robust and that the number of iterations roughly linearly increases with frequency, when the grid spacing is adapted to keep a constant number of discretization points per wavelength. To evaluate the interest of this new iterative solver, the number of floating point operations required for the frequency-domain imaging algorithms and the time-domain imaging algorithms is roughly approximated. This rough estimate indicates that the time-domain migration approach is more than one order of magnitude faster. The velocity model building problem based on a least-squares formulation and a scale separation approach is about three time faster in the frequency domain.
1 INTRODUCTION
Imaging the Earth from seismic data require many (thousands of) solutions of the wave equation because multi-shot and multi-receiver data are used to constrain the algorithm. To process large three-dimensional problems, high-frequency (ray) approximations are commonly made. With the increase of the computer power, finite-frequency solutions of the wave equation is now used when high-frequency approximations are not suitable. The finite-frequency solutions are mainly obtained with a finite-difference approach. In order to deal with three-dimensional problems on todays computers, paraxial (or one-way) approximations of the frequency-domain wave equation are used with the depth as the preferred propagation direction. This leads to efficient three-dimensional imaging algorithms. However the paraxial approximation is often inaccurate with large lateral velocity contrasts and with large distances between the source and the receiver.
without high-frequency or paraxial approximations are needed. The wave equation can be solved either in the frequency domain or in the time domain. The choice is mainly dictated by the computation time and memory requirements5. The time-domain wave equation
can be solved with a marching approach in time and an explicit scheme. Currently three-dimensional imaging algorithms, so-called reverse time migrations, are developed based on time-domain solutions and start to be used on large computer clusters. The frequency-domain wave equation, namely the Helmholtz equation in the acoustic case with constant density, leads to a large sparse linear matrix. With two-dimensional problems a direct solver, based on an LU decomposition, can be used. This is very attractive because only one LU decomposition is required to compute the solutions at multiple source locations. This leads to attractive and efficient imaging methods that are faster than the time-domain imaging methods5,7. In three-dimensional spaces, due to the size of the matrix,
the situation is different and an iterative solver is required.
In practice, three-dimensional frequency-domain solutions of the wave equation are not used because most of the iterative solvers are not robust at seismic frequency or are too expensive. The difficulty arises because the matrix of the linear system is not definite, namely it has positive and negative eigenvalues and because the model size corresponds to several hundreds of wavelength in the three space directions, which means that several hundreds of discretization points in each direction are needed. The efficiency of the iterative method depends on the preconditioner. Recently a new preconditioner based on a heavily damped wave equation and a multi-grid solver has been proposed by Erlangga et al.4. The iterative solver is a preconditioned BI-CGSTAB. The preconditioner corresponds
to the approximate inverse of a heavily damped wave-equation after one multigrid cycle. Preliminary results obtained with this solver on two-dimensional problems have been presented by Riyanti et al.11.
Here a three-dimension implementation of this new iterative solver of the wave equation is discussed together with its relevance for the imaging problem in exploration geophysics. Contrary to Erlangga et al.4,11 the multigrid method is implemented with a standard
full-weighting coarsening, a tri-linear interpolation, a matrix-free algorithm, and a 8th
2 Method
2.1 Iterative solver
The acoustic constant density wave equation corresponds to the Helmholtz equation:
− ω
2
v2(x)u(x, ω) − ∆u(x, ω) = s(x, ω) , (1)
with u the pressure field, s the source function, v the velocity, ω the angular frequency and x the spatial coordinates. k = ω
v is the wave number.
A finite-difference discretization with a J order scheme in space gives:
PJ j=1 aJ j(ul−j,m,p+ ul+j,m,p) + bJj(ul,m−j,p+ ul,m+j,p)+ cJ j(ul,m,p−j+ ul,m,p+j) + (dJ − kl,m,p2 )ul,m,p = sl,m,p . (2) aJ
j, bJj, cJj and dJ are the coefficients of the scheme. ul,m,p, sl,m,p, and kl,m,pare the discrete
values of the pressure field, the source function, and the wavenumber. Equation 2 is valid for i, j, k ∈ [J, nx − J] × [J, ny − J] × [J, nz − J]. Outside, the field values are set to
0. In practice, absorbing boundary layers are added to the domain to minick an infinite domain. In the boundary layers, a complex part to the wavenumber is added; this means that the wave field is damped. Therefore the field can be assumed to be 0 close to the boundaries.
In a matrix from, equation 2 reads
A(k)u = s , (3)
where u and s are now the pressure and source vectors. A is generally non-definite, because the left-hand side of equation 1 is the sum of a positive operator (−∆) and a negative operator (−k2I, with I the identity). With n the number of discretization points in one
spatial direction, A is an n3× n3 sparse matrix. With realistic geophysical examples, n is
between 200 and 1000.
Due to the size of the matrix, a direct solver cannot be used. The logical approach to solve the linear system 3 is an iterative method based on a Krylov method, for instance the BI-CGSTAB method. The efficiency of the solver relies on the preconditioner. The preconditioner system, with the preconditioner M, is
A(k)M−1v = f , with v = Mu . (4)
With the iterative solver the vector A(k)M−1v needs to be evaluated. This is carried out
in two steps; first the system v = Mu is solved, then the matrix-vector product A(k)u is computed. In order to efficiently solve the system v = Mu, Erlangga et al.4 proposed to
use one multigrid cycle and to consider M = A(k(βr+ iβi)), with i =
√
βi two (positive) real numbers. The choice of βr and βi should be such that M is close
enough to A and that the multigrid method is applicable to M. The first condition means that βr should be close to 1 and βi as small as possible. The second condition means that
βi should be large enough to avoid the difficulties associated with the non-definiteness of
A(k). It appears that βr close to 1 and βi between 0.25 and 1 generally lead to a correct
choice.
Physically M corresponds to the discretization of the wave equation with a wave number equal to k(βr+iβi), namely the discretization of a heavily damped wave equation. Indeed,
assuming k constant, the solution of the damped wave equation is a ei k(βr+iβi) ||x|| with a
the amplitude (geometrical divergence) term and ||.|| the L2-norm . βi induces an extra
damping equal to e−kβi||x||. Over a wavelength, i.e. for ||x|| = 2π
k , the extra damping
is 0.04, with βi = 0.5. Over a wavelength, 96% of the energy is lost due to this extra
damping. The heavily damped wave equation almost behaves like a diffusive equation since the energy hardly propagates after a wavelength. This explains why the multigrid method can be used on M.
2.2 Implementation
In order to handle large problems, the implementation of the method differs from the one of Erlangga et al.4. A matrix-free approach is used. The boundary conditions are
Dirichlet boundary conditions and absorbing layers are added to the model to mimic an infinite space.
The idea of multigrid is to estimate the oscillatory part of the solution on the fine grid and to estimate the smooth part of a coarser grid. With h the grid spacing on the fine grid, the discretization of the damped wave equation, equation 2 after replacing k by k(βr+ iβi), leads to the system:
vh = Mhuh . (5)
The superscript h means that the quantities are related to the grid with the spacing h. An approximation, wh, of uh, with the smoother operator, Sh(Mh, vh, uh
0), is computed:
wh = Sh(Mh, vh, uh0) . (6)
uh
0 is an initial guess for uh, generally uh0 = 0. In the implementation the smoother is one
iteration of the symmetric Gauss-Seidel method. wh contains an approximation of the
oscillatory part of the solution, but the smooth part of the solution is difficult to obtain from the Gauss-Seidel method. The residual, rh = vh− Mhwh, satisfies:
rh= Mheh, with eh = uh
− wh . (7)
When whapproximates the oscillatory part of uh, rhand ehhave mainly smooth variations.
rh and eh should then be correctly approximated on a coarser grid with a spacing 2h. With
operator from the coarse grid to the fine grid, the system rh = Mheh becomes: r2h= Ch2hM h P2hhe2h with r2h = Ch2hr h and eh = P2hhe2h . (8) In the implementation, the coarsening operator is the full-weigthing operator and the prolongation operator is the tri-linear interpolation.
The matrix-free implementation consists of approximating C2h
h MhP2hh by M2h, the
dis-cretization of the damped wave equation on the grid with a 2h spacing. Therefore the system on the coarse grid is
r2h= M2he2h . (9)
This equation can be solved with the same approach that the one for the system 5, namely by applying recursively the two-grid approach. The interpolation of e2h can add
suspicious oscillations in eh. To remove those oscillations, the smoother can be applied
on the approximate solution wh+ eh.
A multigrid cycle defines the way to go from the fine grid to the coarsest grid and to return to the fine grid. This gives an approximation of uh. To completely solve the
system 5 the cycle can be applied iteratively. In the implementation, the system u = Mv is approximated with one cycle. There are many possibilities to go back and forth from the fine grid to the coarsest grid. The full multigrid V-cycle is chosen2.
Notice that since the implementation is matrix-free, it is easy to implement high order schemes.
3 Three-dimensional numerical results
3.1 A salt dome example
To illustrate the method, the wave propagation through a salt dome model of 8920 m by 4440 m by 5100 m is computed. The velocity model is displayed in Figure 1. The real part of the pressure wave field at 10.625 Hz and 21.25 Hz is shown in Figures 2 and 3. The source is an explosive source located in the middle top of the model, at the point (4000 m,2000 m,10 m). The grid spacing is 40 m at 10.625 Hz and 20 m at 21.25 Hz. The minimum wavelength being 160 m at 10.625 hz and 80 m at 21.25 Hz, this corresponds to a discretization with 4 points per minimum wavelength, which is normally sufficient with a 8th order scheme. On each side, an absorbing layer with 35 points is added. Therefore
the model has 294×182×199 points at 10.625 Hz and 517×293×326 points at 21.25 Hz. For the preconditioner the wavenumber has been multiplied by βr+ iβi = 1 + 0.5i. It is
Figure 1: Salt velocity model. The velocity varies from 1700 m/s (blue) to 4900 m/s (brown).
Figure 2: Real part of the pressure field at 10.625 Hz. (The source is located at (4000 m,2000 m,10 m).)
3.2 Convergence results
Figure 3: Real part of the pressure field at 21.25 Hz. (The source is located at (4000 m,2000 m,10 m).) 0 50 100 150 200 250 300 350 10−6 10−5 10−4 10−3 10−2 10−1 100 101
Number of bicgstab iterations
Normalised error
Figure 4: Convergence of the BI-CGSTAB for the 10.625 Hz example, blue line, and the 21.25 Hz example, red line.
per wavelength has been kept constant. This result seems to indicate that the number of iterations is proportional to the average number of discretization points in one direction. This average number, called n, is simply obtained by taking the cubic root of the total number of discretization points.
100 150 200 250 300 350 400 450 50 100 150 200 250 300 350
Average number of points in one direction
Number of bicgstab iterations
Figure 5: Convergence results on three three-dimensional examples. The blue lines correspond to the model Figure 1, the black lines correspond to the SEG/EAGE salt model, and the red lines correspond to the Overthrust model. The thick lines correspond to the number of iterations versus the average number of discretization in one direction when the stopping criterion is 10−4and the dashed lines when
the stopping criterion is 10−5.
models, the SEG/EAGE salt model and the Overthrust model1 (not shown here). During
all these computations the number of points per minimum wavelength is kept constant. For the SEG/EAGE salt model and the Overthrust model, a discretization of 5 points per wavelength was used. In Figure 5 the number of BI-CGSTAB iterations versus n, the average number of discretization in one direction, is plotted for the three different models and for two stopping criteria values, 10−4 and 10−5. In practice, it appears that 10−4
is sufficient. Figure 5 numerically shows that the number of iterations roughly increases with the frequency because n is proportional to the frequency when the number of points per wavelength is kept constant. Since the multigrid solver is in O(n3), one BI-CGSTAB
iteration requires O(n3) iterations. This result then shows that the new iterative solver
for the Helmholtz equation is in O(n4).
Notice that the solver has converged on all the examples, which indicates that it is robust. For the SEG/EAGE model and the Overthrust model, the largest model that fits on 16 GB of RAM corresponds to the 12 Hz case and a model of about 80 million unknowns. 4 The application
The numerical results seem to indicate that the new iterative solver is more efficient and robust than the ones based on the ILU decomposition or separation of variables4,8.
However is it good enough for geophysical applications?
cross-correlation of the incident field with the propagated field, obtained by back-propagating the shot data into the Earth. Mathematically, this means to compute the gradient of the least-squares misfit function between the computed synthetics and the data. In practice, the problem is split into two parts. The first problem is the “mapping” in depth of the events recorded in time, this is the so-called migration. The smooth part of the velocity, corresponding to the low spatial frequencies, is assumed known. The second problem, called velocity model building problem, is the determination of the smooth part of the velocity.
The migration problem requires to process a sufficiently large frequency band to obtain an image with a good resolution. With offset data, it was nevertheless shown that the Nyquist criterion could be relaxed, meaning that the number of frequencies could be reduced6. Generally the migration problem does not require an iterative approach because
the least-squares functional is almost quadratic in the perturbations of the model around the known smooth part.
The velocity model building problem does require an iterative approach. However, the minimization of the least-squares misfit functional with respect to the smooth part of the velocity may be done per frequency, or per small overlapping frequency band7,10. The
scale separation is more difficult to carry out when working in the time domain3. It is
then useful to approximately count the number of floating point operations to determine in which domain the migration and the model building problem should be done.
4.1 Rough estimation of the number of operations of the time-domain ap-proach
It is assumed that the time-domain wave equation is solved with an explicit scheme, with a 8th order scheme in space and a 2nd order scheme in time. The evaluation of
a time step then requires about 70n3 operations (a multiplication and a addition are
counted for one operation and a division for 5.). The total number of operations is roughly 70 ntn3. For stability reason, the time sampling, ∆t, and the spacing, ∆x, satisfy
∆t = α1v∆xmax, with vmax the maximum velocity. α1 is generally smaller than 1 depending
on the particular implementation, here we assume α1 = 0.2. By definition the length of
the model is L = n∆x, and the recording time is T = nt∆t. In practice, the recording
time is also linked to the length of the model, for instance T = α2vminL , with vmin the
minimum velocity and α2 around 2. This gives nt= αα21
vmax
vminn. With
vmax
vmin ≃ 3, nt= 30n.
This crude calculation gives the total number of operations for the time-domain solver, Ntop:
4.2 Rough estimation of the number of operations of the frequency-domain approach
Whereas in the time-domain approach the number of grid points is fixed for a given frequency band (and related to the highest frequency modeled during the computation), the number of grid points in the frequency domain can be adapted to the current frequency to save computation time. Let us call n the (average) number of discretization points in one direction associated with a given frequency.
The cost of the frequency-domain approach depends on the cost of one BI-CGSTAB iteration and the number of iterations. One BI-CGSTAB iteration mainly consists of 12 vector-vector operations and 2 matrix-vector operations. Since the calculus is carried out in the complex space, the 12 vector-vector operations represent about 85 n3 (scalar)
operations. A matrix-vector operation is constituted by one multigrid cycle, to solve u = Mv, and one wave-operator evaluation Au. The cost of a multigrid cycle is mainly due to the computation done on the fine grid, since the cost is divided by 8 at each coarsening. More precisely, the cost of the multigrid cycle is the cost on the fine grid multiplied by P
j 81j = 87 ≃ 1.15. On the fine grid and with the current implementation,
the main operations are two Gauss-Seidel evaluations (corresponding to one post-smoother with a symmetric Gauss-Seidel), one wave-operator evaluation (related to the computation of the residual), 2 restrictions (one for the pressure field and one for the velocity), one prolongation (for the pressure field) and one vector-vector operation. A Gauss-Seidel evaluation is more or less similar to the wave-operator evaluation. With a 8thorder scheme,
this means about 140n3. A restriction or a prolongation take about 50n3operations. After
adding all these operations the number of floating point operations of one BI-CGSTAB iteration is roughly γ1n3, with γ1 ≃ 1500.
From the numerical results and the convergence plot, Figure 5, we have nit = γ2n, with
γ2 ≃ 0.5 and nit the number of BI-CGSTAB iterations.
This crude calculation gives the total number of operations for the iterative solver, Nfop:
Nfop = γn4, with γ ≃ 800 . (11)
4.3 The imaging problem
From a computational point of view the synthetics depends on n. The number of discretization points per direction depends on the frequency in the frequency domain and on the maximum frequency in the time domain. We can then write with c(n(f ), v), the synthetics:
B(n(f ), v)c(n(f ), v) = s , (12)
where B is the wave operator and f the current frequency in the frequency domain or the maximum frequency in the time domain. The least-squares misfit function is simply, with d the data:
J(v) = 1
2||c(n(f), v) − d||
and the gradient ∂J(v) ∂v =< λ(n(f ), v), ∂B(n(f ), v) ∂v c(n(f ), v) > (14) with B(n(f ), v)∗λ(n(f ), v) = (c − d) . (15)
<, > is the scalar product; in the time domain, this would corresponds to the cross-correlation. B(n(f ), v)∗ is the reverse time operator in the time domain and the adjoint
operator in the frequency domain. The system 15 is solved as the system 12. Therefore if the memory or I/O requirements are not taken into account, the computational cost of the gradient is roughly equal to twice the computational cost of the wave operator.
4.3.1 The migration problem
As mentioned before the migration problem requires to process a large enough fre-quency band. Let us call fmax the maximum frequency.
In the time domain, the least-squares misfit functional is 1
2||c(n(fmax), v) − d||
2 . (16)
Therefore the number of operations for the time-domain migration is roughly
Nt,migop ≃ 2Ntop= 2αn4(fmax) . (17)
.
In the frequency-domain, the least-squares misfit functional is
nf X j=1 1 2||c(n(fj), v) − d|| 2 , (18)
with nf the number of frequencies and fj ∈ [0, fmax]. When the number of grid points
per wavelength is kept constant, we have n(f ) = fj
fmaxn(fmax). For the migration, the full
frequency band needs to be used to obtain a good resolution. fj = j∆f , with ∆f the
frequency spacing and nf = fmax∆
f . Using equation 11 and the fact that
Pnf
j=1j4 ≃ j
5
5, the
number of operations for the frequency approach is γ5nfn4(fmax). The Nyquist theorem
tells us that ∆f = T1, with T the recording time. As mentioned previously T = α2n(fmax)∆x
vmin .
A discretization of 5 points per minimum wavelength gives 5∆x = vmin
nf∆f. Therefore
nf = α52n(fmax). The total number of operations for the migration problem with the
Helmholtz solver is then roughly
Nf,migop = 2βn(fmax)5 with β =
γα2
This rapid estimation of the number of operations required for the migration problem shows that the number of operations for the frequency-domain approach divided by the number of operations for the time-domain approach is about 0.04n. In practice, it is possible to reduce the number of frequencies (by a factor smaller than 2) in the frequency-domain approach if the data contains large offsets6. This estimation indicates that the
time-domain approach is roughly one order of magnitude faster for medium size problems and will be even more advantageous for larger problems.
4.4 The velocity model building problem
The velocity model building problem requires the minimization of the least-squares misfit functional. Unfortunately the least-squares functional has many local minima. Since only a gradient type of optimization is generally affordable, the initial guess should lie in the domain of attraction of the global minimum. The size of the domain of attraction is roughly proportional to the inverse of the maximum frequency present in the data3. To
take advantage of this behaviour, an idea is to first work with the low frequencies of the data and then to progressively increase the frequency content. In the frequency domain this means that the minimization is carried out per frequency (or per narrow frequency band)7,10. In the time domain, this means that the minimization is carried out per
low-pass filtered data. Given a frequency set (fj), the velocity model building problem then
consists of a sequence of minimization problems: ˜
v(fj) = argminv
1
2||c(n(fj), v) − d||
2 (20)
with ˜v(fj−1) as starting model.
The number of operations to compute c(n(fj), v) is given by equation 10 for the time
domain and equation 11 for the frequency domain. The two approaches have the same complexity. However, the number of operations of the frequency-domain approach is about a thrid of the one of the time-domain approach.
5 CONCLUSIONS
A three-dimensional iterative solver for the wave equation with a preconditioned BI-CGSTAB method has been studied in a geophysical context. The preconditioner consists of a multigrid cycle applied to a heavily damped wave equation. This leads to a ro-bust method even with a relatively high seismic frequency. A numerical study based on three three-dimensional examples shows that the number of iterations of the BI-CGSTAB methods roughly linearly increases with frequency. This confirms the results obtained on two-dimensional examples.
Helmholtz solver and the time-domain approach on an explicity scheme. Indeed, after a rough count of the number of floating point operations required for the computation of the wave-equation solution, it was notice that the time-domain approach will be at least one order of magnitude faster for realistic size problems. The situation is different for the velocity model building based on a least-squares formulation and a scale separation. For this problem, the two approaches have the same complexity, however the frequency-domain approach requires about a thrid of the number of floating point operations of the time-domain approach.
The imaging problem requires the processing of multi-shot data. It may then still be possible to speed-up the frequency-domain approach since the wave equation needs to be solved for multiple right hand sides. This remains an open question. In this study, the memory and I/O requirements have not been taken into account. In a practical implementation of the imaging algorithm, those aspects are also crucial.
REFERENCES
[1] F. Aminzadeh, J. Brac, and T. Kunz. 3-D Salt and Overthrust Models. SEG/EAGE 3-D Modeling Series No.1, SEG, (1997).
[2] W.L. Briggs, Van Emden Henson and S.F. McCormick. A multigrid tutorial - Second Edition, SIAM, (2000).
[3] C. Bunks, F.M. Saleck, S. Zaleski and G. Chavent. Multiscale seismic waveform inversion. Geophysics, 60, 1457-1473, (1995).
[4] Y.A. Erlangga, C.W. Oosterlee and C. Vuik. A novel multigrid based preconditioner for heterogeneous Helmholtz equation. SIAM J. Sci. Comp., Vol. 27, 1471–1492, (2006).
[5] K. Marfurt. Accuracy of finite-difference and finite-element modelling of the scalar and elastic wave equation. Geophysics, 49, 533–549, (1984).
[6] W.A. Mulder and R.-E. Plessix. How to choose a subset of frequencies in frequency-domain finite-difference migration. Geophysical Journal International, 158, 801-812 (2004).
[7] S. Operto, J., Virieux, J.X., Dessa and G. Pascal. High-resolution crustal imaging from multifold ocean bottom seismometers data by frequency-domain full-waveform tomography: application to the eastern Nankai trough. submitted to Journal of Geo-physical Research, (2005).
[9] R.-E. Plessix and W.A. Mulder. Frequency-domain finite-difference amplitude-preserving migration. Geophysical Journal International, 157, 957-987, (2004). [10] R.G. Pratt, Seismic waveform inversion in frequency domain. Part I: theory and
verification in a physical scale domain. Geophysics, 64, 888-901, (1999).
