# Two migration methods based on paraxial equations in a 3D heterogeneous medium

Eliane Bécache, Francis Collino, Michel Kern, Patrick Joly

► To cite this version:

Eliane Bécache, Francis Collino, Michel Kern, Patrick Joly. Two migration methods based on paraxial equations in a 3D heterogeneous medium. Mathematical Methods in Geophysical Imaging III, 1995, San Diego, United States. pp.12, 10.1117/12.218503. hal-01316194

HAL Id: hal-01316194

https://inria.hal.science/hal-01316194v1

Submitted on 15 May 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# On two Migration Methods based on Paraxial Equations in a 3D Heterogeneous Medium

E. Bécache, F. Collino, M. Kern, P. Joly

INRIA, Domaine de Voluceau-Rocquencourt, BP 105, F-78153 Le Chesnay Cédex

## ABSTRACT

We review recent work on paraxial equation based migration methods for 3D heterogeneous media. Two different methods are presented: one deals directly with the classical paraxial equations, by solving a linear system at each step in depth. The other method derives new paraxial equations that lend themselves to splitting in the lateral variables, without losing either accuracy or isotropy. We also show how to incorporate Berenger's perfectly matched layers in this framework. We detail the discretization schemes, both for the full paraxial equations, and for the newly derived equations.

**Keywords:** Seismic migration, paraxial equations, splitting, finite difference, finite element, numerical dispersion

## 1 INTRODUCTION

This paper summarizes recent work done at INRIA on finite difference migration methods in 3D heterogeneous media. We present two methods, which differ chiefly in how they handle discretization of the lateral variables. Indeed, this is the main computational difficulty. Several authors$^{12,18,11}$ have proposed using splitting in the lateral variables. If done naively, this has the undesirable consequence of introducing difficult to control anisotropy. Indeed, the above authors find it necessary to add specially designed correction factors to regain some degree of isotropy. These corrections appear difficult to justify theoretically.

We deal with this issue in two different ways:

* We recognize that we have to solve a (large, sparse, non-hermitian) linear system at each step in depth, and we strive to do so efficiently. Thus we use modern iterative methods.

* Alternatively, we may try to keep the low costs of the splitting method, but do so in a *consistent* manner. As shown below, this is only possible if we introduce more than 2 directions of splitting. We are thus led to design a new family of paraxial approximations, specifically suited to splitting.

We also show how the ideas of “perfectly matched layers” of Berenger$^4$ can be incorporated into a migration program to provide efficient boundary conditions.

An outline of this paper is as follows. We begin by recalling the derivation of paraxial wave equations, and present the new paraxial equations. We then detail the procedure we use for back propagation : it is based on
splitting, then on a (conservative) Runge-Kutta method. The last section is the most important : it presents our discretization procedures in the lateral variables, first for the full wave equation (this is also presented in more detail in the companion paper¹⁶), then for the modified equations. In each case, we show numerical results to illustrate the method.

## 2 THE CONTINUOUS PROBLEM

### 2.1 The classical paraxial equations

We begin with a solution of the wave equation in the whole space,

$$ \frac{1}{c^2} \frac{\partial^2 v}{\partial t^2} - \frac{\partial^2 v}{\partial x_1^2} - \frac{\partial^2 v}{\partial x_2^2} - \frac{\partial^2 v}{\partial z^2} = 0 \qquad (1) $$

with appropriate boundary and initial conditions. This solution can be split up into two waves, an up-going wave and a down-going wave. The terminology up and down refers to the privileged z direction. In many applications one is interested in computing only one of these waves. In this paper, we are interested in the up-going wave. It is governed by the one way wave equation

$$ \frac{dv^+}{dz} + i \frac{\omega}{c} \left( 1 - \frac{c^2 |k|^2}{\omega^2} \right)^{\frac{1}{2}} v^+ = 0 \qquad (2) $$

with

$$ \hat{f}(k_1, k_2, z, \omega) = \iiint f(x_1, x_2, z, t) e^{i(k_1 x_1 + k_2 x_2 - \omega/c t)} dx_1 dx_2 dt. \qquad (3) $$

A major difficulty with this equation is that it corresponds to a non local pseudodifferential equation and therefore is not very tractable from a computational point of view.

Paraxial approximations consist in approximating the square root $(1 - \frac{c |k|^2}{\omega^2})^{1/2}$ with rational fractions so that the new equation is a local partial differential equation :

$$ \left\{ \begin{array}{l} (1 - |\kappa|^2)^{1/2} \rightarrow (1 - |\kappa|^2)^{1/2}_{ap} \\ \kappa = (\kappa_1, \kappa_2), \quad \kappa_1 = \frac{ck_1}{\omega}, \quad \kappa_2 = \frac{ck_2}{\omega}, \end{array} \right. \qquad (4) $$

this approximation being valid as long as $c|k|/\omega$ remains small enough, i.e. as long as the wave propagates close enough to the z-direction.

A general class of high order approximations has been proposed by Bamberger and al.¹ They are a rational fraction expansion of the form:

$$ (1 - |\kappa|^2)^{\frac{1}{2}} \approx 1 - \sum_{\ell=1}^{L} b^{\ell} \frac{|\kappa|^2}{1 - a^{\ell} |\kappa|^2}; \quad a^{\ell} \ge 0, b^{\ell} \ge 0, 1 \le \ell \le L. \qquad (5) $$

The number L specifies the degree of the approximation and the numbers $a^\ell$ and $b^\ell$ can be chosen so as to achieve various approximation criteria. One simple yet useful choice is obtained from a Padé expansion of the square root (see Engquist and Majda⁹).
Paraxial equations are usually handled in the frequency domain, so in writing the PDE obtained from (5) we will only do inverse Fourier transform in the space variables. This is specially convenient for migration purposes since we will simply have to add all the contributions to get the image, which is just zero time section of $v(x_1, x_2, z, t)$. Also, so as not to deal with a high-order PDE, we follow Bamberger and al.¹ and introduce auxiliary unknown functions, thus obtaining the following system of (L+1) coupled PDEs:

$$ \left\{ \begin{aligned} \frac{\partial v^+}{\partial z} + \frac{i\omega}{c} v^+ - \frac{i\omega}{c} \sum_{\ell=1}^{L} b^\ell \varphi_\ell = 0 \\ \frac{\omega^2}{c^2} \varphi_\ell + a^\ell \Delta \varphi_\ell = -\Delta v^+ && \ell = 1, \dots, L. \end{aligned} \right. \qquad (6) $$

System (6) describes paraxial wave equations of higher order in homogeneous media. The corresponding time-domain problem is well posed problem for strictly positive parameters $a^\ell$ and $b^\ell$ (see Trefethen and Halpern¹⁹). It should be noted that high order Taylor expansions are known to give rise to ill-posed problems.

The two most commonly used approximations are the so-called 15° and 45° approximations. The fifteen degree, or parabolic, paraxial equation is simply the first order Taylor expansion of the square root (it corresponds to $L = 1, b_1 = 1/2, a_1 = 0$). A more accurate approximation is obtained via the first Padé approximation, that is $L = 1, \beta_1 = 1/2, \alpha_1 = 1/4$. This gives the 45° equation.

In order to handle the transport term exactly, we make the classical Claerbout change of unknown functions

$$ u = v^+ e^{i\omega/cz}, \qquad (7) $$

and end up with the following system, which we find convenient to write in operator form (this formally eliminates the $\varphi$ unknown):

$$ \left\{ \begin{aligned} \frac{\partial u}{\partial z} - i \frac{\omega}{c} \sum_{\ell=1}^{L} A^\ell u = 0 \\ A^\ell = -b^\ell (a^\ell \Delta + \frac{\omega^2}{c^2})^{-1} \Delta. \end{aligned} \right. \qquad (8) $$

Equation (8) looks roughly like a Schrödinger equation, with a complicated evolution operator (for instance one could show, at least formally, that the energy of the solution is conserved). Since the evolution operator is written as a sum of simpler operators, it is a natural idea to use a splitting method to solve (8). Section 3.1 describes the algorithm in more details. Let us just mention here that after discretization in depth, we have to solve, at each depth step, a linear systems, with a large, sparse, complex valued, non-hermitian matrix. This can be done, as detailed in the companion paper,¹⁶ where we show how this system can be efficiently solved using modern iterative methods.

Alternatively, we can avoid the difficulty altogether, by introducing a new class of paraxial equations suitable for use with splitting in the lateral variables, and requiring only the solution of 1D PDEs. This is the topic of next section.

## 2.2 New paraxial approximations

Using splitting in the horizontal variables is not a new idea. In order to avoid the inversion of the large matrix when using the full paraxial approximation, several authors⁵,¹²,¹⁸,¹¹ have advocated approximating the square root in (4) with:

$$ (1 - |\kappa|^2)^{\frac{1}{2}} \approx 1 - \frac{1}{2} \frac{\kappa_1^2}{1 - \frac{1}{4} \kappa_1^2} - \frac{1}{2} \frac{\kappa_2^2}{1 - \frac{1}{4} \kappa_2^2}. \qquad (9) $$
which is only consistent with the forty-five degree equation in the $x_1 = 0$ and $x_2 = 0$ directions. The problem is that there is no exact decomposition of the form

$$ \hat{R}(\kappa_1, \kappa_2) = \hat{R}_1(\kappa_1) + \hat{R}_2(\kappa_2) \qquad (10) $$

and the conclusion of this brief analysis is that no alternating directions method by splitting exists for the forty-five degree paraxial equation.

The challenge is then to find a method which gives results of comparable accuracy to the forty-five degree equation and whose cost is comparable to that of the fifteen degree equation with splitting.

The basic idea is to introduce more than two directions of splitting. The paraxial equations constructed in Collino and Joly⁸ are derived from an approximation of the square root with rational fractions of the following form:

$$ (1 - |\kappa|^2)^{1/2} = (1 - (\kappa_1^2 + \kappa_2^2))^{1/2} \approx 1 - R(\kappa_1, \kappa_2) $$

with

$$ R(\kappa) = \sum_{j=1}^{N_D} \sum_{\ell=1}^{L} \frac{b_j^{\ell}(\kappa \cdot n_j)^2}{1 - a_j^{\ell}(\kappa \cdot n_j)^2} \qquad (11) $$

where $N_D$ corresponds to the number of directions, $L$ the number of fractions per direction, and $n_j$ the unit vector associated to the $j^{th}$ direction ($n_j = (\cos \alpha_j, \sin \alpha_j)$). It has been shown by Collino and Joly⁸ that the conditions on the coefficients $a_j^{\ell} > 0$; $b_j^{\ell} \ge 0$ ensures the well-posedness of these paraxial equations. Also, in order to keep the approximation error in the square root from blowing-up in certain directions, one has to impose $0 < a_j^{\ell} \le 1$.

Collino and Joly⁸ have constructed several families of approximations of the above type so as to achieve comparable accuracy to the classical forty-five or sixty degree approximations. In particular, they obtained a family of forty-five degree approximations depending on one parameter using four directions of splitting ($N_D = 4$) and one fraction per direction ($L = 1$). For instance, the choice $a_j = 1/3, b_j = 1/4, i = 1, \dots, 4$ gives them the “maxi-isotropic” forty-five degree approximation. More examples are given in the above paper.

The paraxial equation corresponding to (11) can now be written in a form analogous to (8):

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u}{\partial z} - i \frac{\omega}{c} \sum_{j=1}^{N_D} \sum_{\ell=1}^{L} A_j^{\ell} u = 0 \\[1.5em] \displaystyle A_j^{\ell} u = -b_j^{\ell} \left( \frac{\omega^2}{c^2} + a_j^{\ell} D_j^2 \right)^{-1} D_j^2 u \end{array} \right. \qquad (12) $$

where we have defined $D_j = n_j \cdot \vec{\nabla}$.

In practice, a mesh being given in the $x_1, x_2$ plane, there will be two main choices for $N_D$: $N_D = 4$ if the mesh is built from squares, and then the directions are given by the two coordinate axes and the two main diagonals, or $N_D = 3$ if the mesh is built from equilateral triangles, and the 3 directions are 60° apart.

Notice that each of the operators $A_j^{\ell}$ in (12) only involves a one-dimensional differential operator. Thus this new family equation lends itself to a splitting method in the horizontal variables, but as mentioned before, the splitting has the advantage of being consistent with the forty-five degree equation (for a proper choice of coefficients).

Last, notice that formally (8) and (12) have the same form:

$$ \frac{\partial u}{\partial z} - i \frac{\omega}{c} \sum_{m=1}^{N} B_m u = 0 \qquad (13) $$

with $B_m$ and $N$ being either $A^{\ell}$ and $L$ in (8) or $A_j^{\ell}$ and $N_D \times L$ in (12). This enables us, in the next section, to give a formal treatment of the splitting algorithm that embodies both types of equations.
## 2.3 Extension to heterogeneous media and design of absorbing boundary conditions

The extension to heterogeneous media as well as the treatment of the absorbing boundaries are simply obtained by using the same kind of modifications to the operators.

Paraxial equations in heterogeneous media have been proposed and analyzed by Bamberger and al.² Their approach was to define several criteria (both mathematical and geophysical), and to select among a general class of possible candidates the one that satisfied those criteria.

Their result gives a recipe which allows one to extend any paraxial equation to heterogeneous media. Thus equations (8) and (12) keep the same form, with the operators $A^l$ and $A_j^l$ becoming $-b(\frac{\omega^2}{c^2} + a\delta^c)^{-1}\delta^c$ and we have defined

$$ \delta^c = \begin{cases} \frac{1}{c}\mathrm{div}(c\nabla) & \text{for the full paraxial equation} \\ \frac{1}{c}D_j(cD_j) & \text{for the new paraxial equation in the } j^{\text{th}} \text{ direction} \end{cases} \qquad (14) $$

A problem of practical importance is the treatment of the lateral boundaries. It must be designed in such a way that the waves are absorbed when they reach the boundaries. Recently, Collino,⁷ proposed to adapt a novel technique, introduced for electromagnetism by J.P. Bérenger,⁴ This technique consists in designing an absorbing layer model denamed as perfectly matched layer (PML). It possesses the astonishing property to generate no reflexion at the interface between the free media and the artificial lossy medium, and the reflected waves are only due to the discretization of the model. This property allows one to use a very high damping parameter inside the layer, and consequently a short layer length, while still achieving a quasi-perfect absorption of the waves. Practically, the PML is very easy to implement : we change equation (14) to

$$ \delta^c = \begin{cases} \sum_{i=1}^{2} \frac{1}{cd(x_i)} \partial_{x_i} (cd(x_i)\partial_{x_i}) & \text{for the full paraxial equation} \\ \frac{1}{cd(x_j)} D_j (cd(x_j)D_j) & \text{for the new paraxial equation in the } j^{\text{th}} \text{ direction} \end{cases} \qquad (15) $$

with

$$ d(x) = \frac{i\omega}{i\omega + \sigma(x)}. \qquad (16) $$

$\sigma(x)$ is a positive function with support in the damped area. We close the system with a Dirichlet boundary condition at the end of the layer. The efficiency of the method has been assessed by numerical experiments : a quasi-perfect absorption has been obtained with only 4 extra-nodes at the boundary and an appropriate choice for $\sigma$.

# 3 EVOLUTION IN DEPTH

The construction of a discrete scheme for the evolution in depth proceeds in two steps. We first perform a splitting on equation (13), to obtain a family of problems with a single $B_m$ operator, then apply standard ODE techniques to this last problem.
## 3.1 Splitting

We briefly recall some classical points about splitting methods. The exact solution of (13) satisfies:

$$u(z + \Delta z) = \exp\left(\frac{i\omega}{c} \sum_{m=1}^{N} B_m \Delta z\right) u(z). \qquad (17)$$

The splitting methods consist in approximating the exponential of the sum of operators with the product of exponentials:

$$u_{ap}(z + \Delta z) \approx e^{\frac{i\omega}{c} B_N \Delta z} \times \dots \times e^{\frac{i\omega}{c} B_1 \Delta z} u(z) \qquad (18)$$

and are second order accurate with respect to the discretization in the z variable, in the case of an heterogeneous medium. Approximation (18) then leads naturally to define N intermediate values $w^m, m = 1, \dots, N$ satisfying:

$$\left\{ \begin{array}{l} \frac{dw_m}{dz} - \frac{i\omega}{c} B_m w_m = 0 \\ w_m(0) = w_{m-1}(\Delta z) \equiv w^{m-1}. \end{array} \right. \qquad (19)$$

Finally we have $u_{ap}(z + \Delta z) = w^N$.

The problem to be solved at each step is still an evolution problem in depth, but with a single operator. In the next subsection, the discretization in the evolution variable is performed in order to get a high order semi-discretized scheme.

## 3.2 Semi discretization in depth

The discretization in the transverse variable used here has been initially proposed by Joly and Kern.¹⁴ It is based on the fact that the exact solution of (19) satisfies $w^m = e^{iB_m \Delta z} w^{m-1}$, and on the relationship between Runge-Kutta methods and Padé approximations to the exponential (see for instance¹³). In order to get conservative schemes of order $2K$, the exponential is replaced with a Padé approximant on the following form:

$$\exp(ix) \approx \prod_{k=1}^{K} \frac{1 + r_k x}{1 + \bar{r}_k x}$$

The integration from 0 to $\Delta z$ is then done formally as follows:

$$w^m = \prod_{k=1}^{K} (I + \bar{r}_k \Delta z B_m)^{-1} (I + r_k \Delta z B_m) w^{m-1}$$

This procedure leads us to define K intermediate unknowns $w_k^{m-1}$ associated to each fraction, solutions of:

$$(I + \bar{r}_k \Delta z B_m) w_{k+1}^{m-1} = (I + r_k \Delta z B_m) w_k^{m-1}$$

We then set $w^m = w_K^{m-1}$.

The classical Crank-Nicolson second-order scheme is obtained for $K = 1$ and $r_1 = i/2$. Kern¹⁵ has also shown that this process can also be used to get non conservative schemes, which generalize the $\theta$-scheme to higher order schemes. Following Lambert,¹⁷ one can write family of third order schemes, which provide some damping with respect to the pure conservative fourth order scheme.
There remains the task of discretizing the operator $B_m$ with respect to the lateral variable. Let us recall that for the full paraxial equation, the $B_m$s involve 2D differential operators, whereas for the new equations, they involve 1D operators. Thus, after discretization in the lateral variables, the linear system to be solved will be much larger for the full equation than for the new equation. Moreover, we will see below that the matrices for the new equations are banded.

Finally, one elementary step for the full paraxial equation requires the solution of a system of the form

$$(I + d(\omega)\Delta_h)U = R$$

while the new equation requires the solution of several simpler systems of the form

$$(I + d_j(\omega)D_{j;h}^2)U = R,$$

where $\Delta_h$ and $D_{j;h}^2$ are approximations to the (heterogeneous) Laplacian, and the operator $D_j^2$ respectively. We devote next section to specifying how this approximation is achieved.

# 4 DISCRETIZATION IN THE LATERAL VARIABLES

In view of section 3.1, we now take for simplicity, and without loss of generality $L = 1$ and (for the new paraxial equation) one particular direction, which we denote by $x$. Thus we have to solve

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial w}{\partial z} - \frac{i\omega}{c}\varphi = 0 \\[1.5ex] \displaystyle \frac{\omega^2}{c}\varphi + \operatorname{div}(c\nabla(a\varphi + bw)) = 0 \end{array} \right. \qquad \text{for the full paraxial equation} \qquad (20) $$

or

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial w}{\partial z} - \frac{i\omega}{c}\varphi = 0 \\[1.5ex] \displaystyle \frac{\omega^2}{c}\varphi + \frac{\partial}{\partial x}\left(c\frac{\partial}{\partial x}(a\varphi + bw)\right) = 0. \end{array} \right. \qquad \text{for the new paraxial equation} \qquad (21) $$

We base the discretization of (20) and (21) on a variational formulation, which we recall below. This provides us with a systematic treatment of heterogeneous media, in a way that insures good numerical properties (ie. convergence and stability).

We let $\Omega$ be a rectangle in $\mathbb{R}^2$ for (20), or an interval in $\mathbb{R}$ for (21). We also denote by $\nabla$ the spatial derivative, irrespective of the space dimension. We set:

$$(u, v) = \int_{\Omega} u \bar{v} \, d\Omega \quad ; \quad m(u, v) = \int_{\Omega} \frac{1}{c} u \bar{v} \, d\Omega \quad ; \quad k(u, v) = \int_{\Omega} c \nabla u \nabla \bar{v} \, d\Omega$$

With these definitions, both problems (20) and (21) admit the variational formulation:

$$ \text{Find } (w, \varphi) \text{ such that } \left\{ \begin{array}{l} \displaystyle \frac{d}{dz}(w, \chi) - i\omega m(\varphi, \chi) = 0 \quad \forall \chi \\[1.5ex] \displaystyle \omega^2 m(\varphi, \psi) - k(a\varphi + bw, \psi) = 0 \quad \forall \psi. \end{array} \right. \qquad (22) $$

Choosing a discretization scheme amounts to replacing the unknown functions $w$ and $\varphi$ by vectors $w_h$ and $\varphi_h$
(h is the mesh size, and indicates discrete quantities), and the bilinear forms m and k by matrices

$$ \left\{ \begin{array}{l} \frac{dw_h}{dz} - i\omega M^* \varphi_h = 0 \\ (\omega^2 M_\alpha - aK) \varphi_h = bK w_h, \end{array} \right. \qquad (23) $$

with $M_\alpha = \alpha M^* + (1 - \alpha)U$, and where $M^*$ and $U$ are mass matrices, and $K$ is a stiffness matrix.

It is important to introduce two different mass matrices $M^*$ and $U$ for the following reasons:

* We need to keep $M^*$ diagonal because we will have to invert it when we eliminate $\varphi$

* On the other hand, one can show that a proper choice of $\alpha$ allows us to gain two orders of accuracy compared with the classical discretization (we come back to this point in sections 4.1 and 4.2).

After introducing discretization in depth, as in section 3.2, we have to solve a set of systems of the form:

$$ S_k w_{k+1} = \bar{S}_k w_k \qquad (24) $$

with $S_k = (\frac{\omega^2}{a} M_\alpha - K)(M^*)^{-1} + \bar{\tau}_k \frac{\omega \Delta z b}{a} K$.

We now give discretization methods for the two cases (20) and (21). We will illustrate the method with a simple example: the migration of a point source in a homogeneous medium. The computational domain is 1250 m in each of the horizontal directions, and 625 m in the vertical direction. The grid sizes are $h = \Delta z = 12.5$ m. The velocity is 1000 m/s. The source is located on the surface at the center of the computational domain, the time of the explosion is 0.5125 s (82% of the travel time to the bottom), and the cutoff frequency of the source is 50 Hz. This is a relatively severe test as the number of points per wavelength at the cutoff frequency is 1.5 (at the central frequency of 20 Hz, it is still around 4). This example is for a constant velocity. Of course the codes do not take any advantage of this fact. There still remains the question of how well our methods can handle varying velocities.

## 4.1 For the full paraxial equations

We base our discretization procedure on a rather special choice of finite elements. This is more fully analyzed in,¹⁶ where we show this choice can be interpreted as a finite difference scheme in the homogeneous medium case. Let us just mention here that the particular combination we select has good isotropy and accuracy properties.

We consider a regular grid in the $x-y$ plane. We can build three different finite element meshes: the right (resp. left) looking triangular mesh, leading to a P1 approximation, and the bilinear mesh leading to a Q1 approximation. Then, the stiffness matrix $K$ and the consistent mass matrix $U$ are computed as:

$$ K = 1/2(K^Q + K^P) \quad \text{and} \quad M = 1/2(M^Q + M^P) $$

where $K^P$ and $M^P$ are themselves the average of the matrices of the right- and left-looking P1 meshes. Also, the matrix $M^*$ is just the lumped version of $U$. As mentioned above, the parameter $\alpha$ allows us (at least in homogeneous media) to increase accuracy at no cost. It is similar to Claerbout's⁶ “$\gamma = 1/12$” trick (see also next section).

We now have to solve (24). This is a large linear system, whose matrix is sparse and non-hermitian. The size of the system (from $10^4$ to $10^6$ unknowns) prevents using direct methods, so we turn to iterative methods.
For symmetric positive linear systems, the (preconditioned) conjugate gradient method is currently considered as the “best” linear iterative solver. For non-symmetric (or non-hermitian) systems, the situation is less clear-cut. No method seems to work in all situations, although recent research has focused on Krylov subspace methods (see Freund and al.¹⁰ We have implemented a simple biconjugate gradient method. For details on the algorithm, see.¹⁶,¹⁵ It is interesting to note that, with a simple diagonal preconditioning, convergence is very rapid (less than 10 iterations) for all but the very few first frequencies. Since, for simplicity, we keep a fixed mesh for all frequencies, we have more than 20 points per wavelength for those troublesome frequencies. A simple explanation seems to be that, at low frequency, we are really solving an almost singular problem.

We show the results of this procedure on the example described above in figure 1. It compares a vertical slice of the migrated image along the *x* coordinate axis, for the second and fourth orders schemes in (x, y) (both are second order in depth). Given the stringency of the test, both schemes give qualitatively good results, though the fourth order scheme does not look any better than the second order one, even showing somewhat more dispersion. The fourth order method regains its superiority when looking at isotropy, as is shown in the companion paper.¹⁶

<figure><img src="image_2.png" /><figcaption>Figure 1: Vertical slices along the *x* axis</figcaption></figure>

## 4.2 For the new paraxial equations

In the homogeneous case, it is possible to get high order finite difference schemes with a reduced number of points by approximating the second-order derivative $d^2/dx^2$ in a proper way. The extension to heterogeneous media is then not clear. Through the variational approach, there is a systematic way to get schemes in a heterogeneous medium, with good properties (convergence, stability). The stiffness operator in the form $k$ to be approximated is then $(d/dx)^t(cd/dx)$. We thus have to find high order approximations of $d/dx$. The inconvenience is that the approximation of the operator $(d/dx)^t(cd/dx)$ requires then more discretization points than this of $d^2/dx^2$.

The discretization is built upon a variational finite differences method. It consists in approximating the operators with operators of finite difference type. Let us introduce some notations. The domain is discretized by a regular grid ($x_i = i\Delta x$) and we define the shifted grid by the nodes ($x_{i+1/2} = (i+1/2)\Delta x$). The approximation is then decomposed on the functions $\chi_i$ where $\chi_i$ is “almost” the characteristic function on $[x_{i-1/2}, x_{i+1/2}]$ or more precisely $\chi_i(x) = 1/\sqrt{h}$ for $x \in [x_{i-1/2}, x_{i+1/2}]$ and 0 elsewhere. At this step, we have two terms to approximate : the mass term associated to the operator $1/cI$ and the stiffness term associated to the operator $(d/dx)^t(cd/dx)$. We then define an approximation of order 2n of $d/dx$, denoted by $\partial_h^{[2n]}$ which leads to the approximate stiffness matrix $K_h^{[2n]} = k(\chi_i, \chi_j) = (c\partial_h^{[2n]}\chi_i, \partial_h^{[2n]}\chi_j)$.
As announced in the beginning of this section, we propose a “modified” scheme obtained from a linear combination of two different approximations of the mass term.

The classical scheme consists in keeping the exact operator $1/cI$ which gives the diagonal mass matrix $M_{ij}^* = m(\chi_i, \chi_j)$. It thus corresponds to the choice $\alpha = 1$ in (23).

For the modified scheme, we define approximations $\delta_h^{[2n]}$ of the operator Identity, of order $2n$ and the corresponding mass matrix $U_{ij}^{[2n]} = (\frac{1}{c}\delta_h^{[2n]}\chi_i, \delta_h^{[2n]}\chi_j)$. The resulting combination $M_\alpha^{[2n]} = \alpha M^* + (1-\alpha)U^{[2n]}$ is of course still a $2n$-order approximation of the mass form. System (23) can be rewritten after elimination of the auxiliary function as

$$ \frac{dw_h}{dz} = i\omega M^*(\omega^2 M_\alpha^{[2n]} - aK^{[2n]})^{-1}bK^{[2n]}w_h \quad (25) $$

After discretization in depth, we solve systems like (24), which have banded matrices, with a small bandwidth. This is the main advantage of the new method, compared to the one described in the previous section.

We have now one degree of freedom $\alpha$ and we can show that for any $n \ge 1$, there exists a choice $\alpha^{[2n]}$ such that the scheme becomes of order $2n+2$ in the homogeneous case. As a matter of fact, we can get an explicit expression for this value, $\alpha^{[2n]} = \frac{2n}{2n+1}$. For more details on the modified schemes, the reader can refer to Bécache et al.³

<figure><img src="image_3.png" /><figcaption>Figure 2: Comparison between several schemes; $x = x_s$</figcaption></figure>

With $n = 1$, the modified scheme corresponds to the classical Claerbout⁶'s scheme, with the relation $\alpha = (1-\gamma)/4$. The modification presented here plays on the mass matrix, but it can be interpreted as a modification on the stiffness matrix as in the Claerbout's scheme, by rewriting the modified mass matrix as

$$ M_\alpha^{[2n]} = (I + (1 - \alpha)(U^{[2n]} - M^*)(M^*)^{-1}) M^* $$
and multiplying the second equation in the approximate system by

$$ (I + (1 - \alpha)(U^{[2n]} - M^*)(M^*)^{-1})^{-1}, $$

in order to reobtain $M^*$ for the mass term. The stiffness matrix is then modified and becomes

$$ (I + (1 - \alpha)(U^{[2n]} - M^*)(M^*)^{-1})^{-1} K. $$

The corresponding value $\alpha^{[2]} = 2/3$ leads then to the well known $\gamma = 1/12$ choice and gives a fourth-order scheme with a matrix of bandwidth equal to 3 instead of 7 for the classical fourth-order scheme. For $n = 2$ and $\alpha = 4/5$, we get a sixth-order one with a matrix of bandwidth equal to 7 instead of 11 for the classical sixth-order scheme. The drawback is that the analysis in heterogeneous media through a-priory energy estimates that can be made for the classical schemes to show their well-posedness is not valid anymore for these modified schemes. However, on a numerical point of view, they seem to behave even better than the same order classical ones and this is also confirmed with a dispersion analysis.

We present a comparison between several schemes applied to the example described above in figure 2. As in the previous section, the figures represent vertical slices of the migrated image along the $x$ coordinate axis. Figures 2(a), (b) and (c) compare the results obtained with a second, fourth and sixth order discretization in the lateral variables and show a significant improvement. On the contrary, the use of higher order discretization in the depth variable does not seem to give better results (see figures (c) and (d)). But one should remember that this is a test with a very small number of points per wavelength, so the results are on the whole quite good.

# 5 CONCLUSION

We have presented two alternative ways to solve efficiently and accurately 3D paraxial equations.

The first way discretizes the usual full 3D equation with higher order schemes. The main difficulty lies in the solution of a large non-hermitian linear system. We have shown how this could be done efficiently by using modern iterative methods. As soon as the frequency is not too low, convergence is rapid enough to make this method competitive with splitting methods. There remains the question of efficiently solving for the low frequencies.

The second way avoids this issue altogether by introducing new paraxial equations that lend themselves to splitting with respect to the lateral variables in a consistent manner. This gives rise to the solution of a sequence of 1D linear systems for each splitting direction.

The cost of the first method is higher than that of the second one, but it gives rise to less numerical dispersion. An explanation is that the splitting method needs to work with the mesh along diagonal directions, for which the number of points per wavelength is somewhat less than that along the coordinate axes. A possible remedy is provided by the use of the modified higher order schemes, whose dispersion is less than that of the classical ones for a smaller computational cost, as shown in the last section.

We want to emphasize that our methods are based on sound computational procedures, and that this does not preclude an efficient solution.

# 6 ACKNOWLEDGMENTS

This research was carried out as part of the Prestack Structural Interpretation consortium project (PSI). The authors hereby acknowledge the support provided by the sponsors of this project.
# 7 REFERENCES

[1] A. Bamberger, B. Engquist, L. Halpern, and P. Joly. Higher order paraxial approximations for the wave equation. *SIAM J. on Appl. Math.*, 48:129-154, 1988.

[2] A. Bamberger, B. Engquist, L. Halpern, and P. Joly. Paraxial approximations in heterogeneous media. *SIAM Journal on Applied Mathematics*, 48:99-128, 1988.

[3] E. Bécache, F. Collino, and P. Joly. Higher order numerical schemes for the paraxial wave equation. Yearly report, PSI Consortium, 1995.

[4] J. Bérenger. A perfectly matched layer for the absorption of electromagnetic waves. *Journal of Comp. Physics.*, 114:185-200, 1994.

[5] D. L. Brown. Applications of operator separation in reflection seismology. *Geophysics*, 1983.

[6] J. F. Claerbout. *Imaging the Earth's Interior*. Blackwell Scientific Publication Co., Oxford, 1985.

[7] F. Collino. Perfectly matched absorbing layers for paraxial equations. Yearly report, PSI Consortium, 1995.

[8] F. Collino and P. Joly. Splitting of operators, alternate directions and paraxial approximations for the 3-d wave equation. submitted to *SIAM J. Sci. Comp.*, 1994.

[9] B. Engquist and A. Majda. Absorbing boundary conditions for the numerical simulation of waves. *Math. Comp.*, 31:629-651, 1977.

[10] R. W. Freund, G. H. Golub, and N. M. Nachtigal. Iterative solution of linear systems. *Acta Numerica*, 1:57-100, 1992.

[11] P. Froidevaux. First results of a 3D prestack migration program. In *60th Ann. Meeting. Soc. Expl. Geophys.*, 1990.

[12] R. W. Graves and R. W. Clayton. Modeling acoustic waves with paraxial extrapolators. *Geophysics*, 55:306-319, 1990.

[13] E. Hairer and G. Wanner. *Solving ordinary differential equations: 2: stiff and differential-algebraic problems*. Springer Verlag, 1991.

[14] P. Joly and M. Kern. Numerical methods for 3D migration. First yearly report, PSI Consortium, 1990.

[15] M. Kern. Numerical methods for the 3D paraxial wave equation in heterogeneous media, part iii. Yearly report, PSI Consortium, 1992.

[16] M. Kern. A nonsplit 3d migration algorithm. In S. Hassanzadeh, editor, *Mathematical Methods in Geophysical Imaging III*. SPIE, 1995.

[17] J. D. Lambert. *Numerical Methods for Ordinary Differential Systems*. John Wiley and Sons, 1991.

[18] Z. Li. Compensating finite-difference errors in 3D migration and modelling. *Geophysics*, 56:1650-1660, 1991.

[19] L. N. Trefethen and L. Halpern. Well posedness of one way equations and absorbing boundary conditions. *Math. Comp.*, 47:437-459, 1986.