The Laguerre finite difference one-way equation solver

The Laguerre finite difference one-way equation solver

Andrew V. Terekhov Institute of Computational Mathematics and Mathematical Geophysics, 630090, Novosibirsk, Russia Novosibirsk State University, 630090, Novosibirsk, Russia

This paper presents a new finite difference algorithm for solving the 2D one-way wave equation with a preliminary approximation of a pseudo-differential operator by a system of partial differential equations. As opposed to the existing approaches, the integral Laguerre transform instead of Fourier transform is used. After carrying out the approximation of spatial variables it is possible to obtain systems of linear algebraic equations with better computing properties and to reduce computer costs for their solution. High accuracy of calculations is attained at the expense of employing finite difference approximations of higher accuracy order that are based on the dispersion-relationship-preserving method and the Richardson extrapolation in the downward continuation direction. The numerical experiments have verified that as compared to the spectral difference method based on Fourier transform, the new algorithm allows one to calculate wave fields with a higher degree of accuracy and a lower level of numerical noise and artifacts including those for non-smooth velocity models. In the context of solving the geophysical problem the post-stack migration for velocity models of the types Syncline and Sigsbee2A has been carried out. It is shown that the images obtained contain lesser noise and are considerably better focused as compared to those obtained by the known Fourier Finite Difference and Phase-Shift Plus Interpolation methods. There is an opinion that purely finite difference approaches do not allow carrying out the seismic migration procedure with sufficient accuracy, however the results obtained disprove this statement. For the supercomputer implementation it is proposed to use the parallel dichotomy algorithm when solving systems of linear algebraic equations with block-tridiagonal matrices.


rmkRemark \newproofpfProof \newproofpotProof of Theorem LABEL:thm2

1 Introduction

The wave equation describes waves propagating in all directions. A one-way wave equation, also known as paraxial or parabolic, wave equation is an equation describing only downward propagating waves or only upward propagating waves for a positive or a negative vertical component of propagation velocity, respectively. Mathematical models based on the wave equation are often considered in problems of seismic prospecting [1, 2, 3], ocean acoustics [4, 5, 6] as well as for setting non-reflecting boundary conditions [7, 8, 9], whereas initially the ideas were formulated for solving the problem of electromagnetic wave propagation along the Earth’s surface[10].

For solving the one-way equation it is possible to use finite difference methods [11, 12]. To this end, the pseudo-differential square-root operator was first approximated by an optimized series which has its origin in a continued fraction expansion [11, 13, 14, 15]. A finite difference method can handle strong lateral velocity variations, however it is dip-limited and has significant numerical dispersions and artifacts. A conventional approach to reducing numerical dispersions is in applying finite difference schemes of a high order of accuracy, but when solving the one-way equation there arise essential computational problems. First, there is a numerical instability of calculations in using schemes of high orders of accuracy for the approximation of spatial derivatives towards the wave field extrapolation. Another difficulty is in that for decreasing the time of calculation, a practical approach, independently proposed in [16] and in [17], to solving extrapolation equation (4) is implemented by the Marchuk-Strang splitting [18]. However for an inhomogeneous velocity model such an approach provides only the first accuracy order as related to a step in depth because the operators to be split are not mutually commutative. In order to attain second order of accuracy, the splitting procedure should be specially symmetrized [18] thus twice increasing computer costs. The abandonment of using the Marchuk-Strang splitting procedure as regards this algorithm demands higher computer costs, especially for a 3D case when it is necessary to use iterative procedures for solving systems of linear algebraic equations (SLAEs). Another approach to reducing numerical dispersion and computing artifacts when applying the finite difference approximation is in the use of the filtration procedures [19] as well as of the nonlinear procedures [20] similar to the TVD-schemes concepts for hyperbolic equations.

For reducing the numerical dispersion when solving geophysics problems the Fourier-based methods were developed [21, 22, 23]. These methods only approximately take into account the lateral velocity variations as being obtained on the basis of non-identical mathematical transformations thus unpredictedly affecting the quality of solution. Nevertheless, for real practical tasks it has been just these methods that make possible to carry out calculations with sufficiently large mesh sizes. In this case numerical artifacts generated by the Fourier-based methods are essentially smaller than those generated by finite difference methods of second order of accuracy. On the other hand, the advantage of the finite difference approximation is in that it possesses convergence and is spatially localized, thus allowing the evaluation of the quality of solutions that are obtained with a sequence of imbedded meshes. In [24], an alternative approach to obtaining solutions similar to that of the one-way equation is proposed. This approach is a combination of various ideas related to the one-way equations, half-space stiffness relation, special finite element discretization and complex coordinate stretching. However, the question of constructing a stable difference method of a high accuracy order for the classical one-way equation is still an open question and urgent both from the standpoint of applications and methods.

In this paper we propose a new finite difference algorithm of high order of accuracy for solving the one-way equation. As opposed to the generally accepted approach we use the integral Laguerre transform [25] instead of Fourier transform with respect to time. As for the approximation of spatial derivatives we employ only finite difference schemes. The approach in question was investigated in solving dynamic problems for the acoustics and elasticity equations [26, 27, 28, 29]. Applying the Laguerre transform as well as Fourier transform with respect to time enables us not to consider time-centering of the method of calculating the wave field that is required in analyzing finite difference approximations of the temporal derivative. Also, the Laguerre transform makes possible to obtain SLAEs with considerably better computational properties, than with using Fourier transform. Decreasing computer costs for solving SLAEs allows the abandonment of the Marchuk-Strang splitting and hence an increase both in stability and accuracy of the algorithm as a whole. In this case the approximation error relatively the mesh size towards the extrapolation of the wave field is of fourth accuracy order thus considerably reducing the numerical dispersion as compared to algorithms accurate to second order. As a result, a widespread opinion that a purely finite difference approach does not allow obtaining high quality solutions to the one-way equation will be invalidated.

2 Description of the method

2.1 Governing wave equations

The one-way equation can be obtained from the acoustics equation


where is the field variable, is the wave velocity, the vertical direction is the extrapolation direction, i.e., the direction of one-way propagation and the positive axis is directed downward, i.e., toward increasing depth. Applying Fourier transforms to equation (1) over and yields


where is a wave component at the radial frequency , is the horizontal wave number. Then, the 2D one-way acoustic wave equation for upcoming waves in the frequency domain can be expressed as:


where is the imaginary unit. In practice, the square-root operator is approximated by an optimized series which has its origin in a continued fraction expansion that can be represented by the ratios of polynomials [11].


and the polynomial coefficients for the propagation angle should be optimized[13, 30].

Let us write down the one-way equation for the spatial-temporal domain . To this end it is convenient to introduce the auxiliary functions[14, 15] of the form

so that equation (4) can be rewritten as

Applying the inverse Fourier transform to the two latter equations, we come to a system of equations determining the down-continuation process


2.2 Temporal approximation

Let us consider the direct and the inverse Laguerre transforms[31]


Here are the orthonormal Laguerre functions, is Laguerre polynomial degree, is the order of the Laguerre functions, and is the transformation parameter. Assuming , we can show [31, 32] that


where for we have


Then applying the Laguerre transform (6) to equations (5a),(5b) we obtain the following system of equations for the calculation of the -th coefficient of expansion:


where , and the index denotes number of a term of series (6).

The expansion parameters from (4),(5) for are chosen the following:  , , , , , , for which, as shown in [30, 2], such approximation is valid up to the angles of degrees. It is also possible to select large values , thus asymptotically approaching to , however in this case computer costs will essentially increase.

The method in question can be considered as analog of the spectral-difference one based on Fourier transform, but the part of ”frequency”  is played by the parameter , which determines the degree of the Laguerre polynomials. An important property is in that opposed to the Fourier method the operator of problem (9) is independent of the parameter separating harmonics . Hence it appears possible to solve in the domain the common system of equations for different right-hand sides whereas this property is not valid in the domain . After carrying out the spatial approximation of equations (9), the multiple solution to the SLAE with a common matrix allows one to construct an effective computing procedure for solving a difference problem.

A disadvantage of the Laguerre transform is the absence of the fast transformation algorithm. However, taking into account the fact that input data are set only along the upper surface (), and the inverse transformation is done for a fixed time instant, the total cost of the direct and inverse transformations appears to be minor as compared to that needed for calculation of coefficients of decomposition (6) from the solution to problem (9).

2.3 Spatial approximation

To approximate equation (9a) we will use unconditionally stable the Crank-Nicolson scheme[33] of second order of accuracy:


where the difference operator is of the form


Formally, the presence of a singular component in the solution [14, 15] does not allow attaining a high order of convergence for finite difference schemes, but for practical purposes the calculations have shown (see Section Numerical Experiments) that schemes of high orders of accuracy make possible to obtain solutions with a considerably lesser number of mesh nodes because such schemes more precisely reproduce the dispersion law. In deriving equations (5) the velocity model was assumed to be homogeneous, although it also yields satisfactory results for inhomogeneous media. In the latter case this model correctly keeps kinematics of waves, but not their amplitudes. For a wide range of problems such an approximate one-way model is admissible, as the correct account of amplitudes essentially increases computer costs [34]. Based on the above-said for approximating it is reasonable to use the dispersion-relationship-preserving method (DRP) by Tam and Webb [35], whose main idea is in the following. According to the Fourier derivative rule, values of optimized coefficients in (11) are defined as solution to the problem of minimizing the error functional in the space of wave numbers

where is some weight function, is a maximum accurate wave number. This approach and its various modifications [36, 37, 38] make possible to decrease the number of mesh nodes and to preserve high accuracy of calculations as compared to conventional difference schemes obtained with the Taylor expansion in series [39].

2.4 Solution of the SLAEs

Solving the difference equations as a stage of the algorithm proposed demands high computer costs, therefore its applicability depends on the speed of solving the corresponding SLAEs. As opposed to Fourier transform, the coefficients of the Laguerre expansion in series (6) are dependent in a recurrent manner (see (9)). Hence, for a fixed for different it is required to solve SLAEs many times with a common real matrix and different right-hand sides. This allows the use of the methods based on decomposition, where factorization of a matrix for all is done only once.

Let us write down the difference problem (10) in the form of a SLAE as


where is the unit matrix. Employing the Schur complement [40], the mesh functions can be calculated through the solution to the following reduced SLAE


where the functions are defined as


Making use of the matrix property [41]


we can make a simplification



Multiplying equation (16) by the matrix and taking into consideration the commutative property of , we obtain the governing equation for the calculation of mesh functions


Matrix (17) is banded and can be explicitly represented without calculation of the matrices thus allowing us to apply efficient algorithms for solving SLAEs based on -decomposition. Hence, we have to calculate different decompositions for (17), where is the number of mesh nodes in the direction . Note, also, that opposed to (16), for calculating the vector in the right-hand side (17) the inversion of the matrix is not required.

After the calculation of the mesh functions before turning to calculating the functions , the functions can be expressed as


Equation (18) is obtained applying the property of (15) to (14). The function should be calculated at the integer mesh nodes, because this simplifies the implementation of the Richardson extrapolation algorithm considered in the next Section.

2.5 Increasing accuracy of the spatial approximation via the Richardson extrapolation

According to the studies [27, 28, 29] it is needed to have a few tens of mesh nodes per a minimum wavelength in order to provide a reasonable accuracy of calculating schemes of second order for real spatial-temporal scales. In order to reduce computer costs, we increase the accuracy of the algorithm up to based on the Richardson extrapolation [42, 43].

Let the mesh functions , , defined on the meshes with the mesh steps and , be solutions to problem (10). Then the linear combination


approximates the solution to problem (9) on the mesh accurate to . The Richardson extrapolation occurs more rarely than difference schemes of high accuracy orders. However, preliminary calculations have shown that multi-step methods of the Adams-Moulton and the Adams-Bashfort types do not provide stability even for small steps . The implicit Runge-Kutta methods of high accuracy orders are not efficient because they demand multiple calculation of the right-hand side of the equation to be solved. The instability of methods of high orders is due to the presence of a singular component in the solution, but according to computational experiments the Richardson extrapolation makes possible to stabilize such instability.

Fourth order downward-continuation algorithm with the global correction To calculate the mesh functions accurate to , the following is necessary:

  1. Based on the cubic splines interpolate values of the functions , preset on the mesh , into nodes .

  2. Based on the cubic splines interpolate values of the functions ,, preset on the mesh , into semi-integer nodes .

  3. On the mesh , applying equation (17), calculate the solution .

  4. On the mesh , applying equation (17), calculate the solution .

  5. Based on the Richardson extrapolation, correct the mesh function with (19).

  6. After the calculation of the mesh functions for all depths , the functions can be defined from the solution to equations (18).

  7. Turn to the calculation of , etc. coefficients of the expansion of the Laguerre series.

Also, we can formulate the second algorithm where the extrapolated value of (19) is calculated immediately after calculating the function , , , then the extrapolated value is used at the next step for calculating , , . This algorithm will be called ”the downward-continuation algorithm with the local correction”. If for calculating the solutions and one uses a numerically stable algorithm (the Crank-Nicolson scheme in our case), the extrapolation procedure based on the global correction will also be numerically stable, whereas in the general case, the stability is not ensured with the use of the local correction [43].

2.6 Parallel consideration

The possibility of the efficient use of modern supercomputer systems is one of obligatory demands imposed upon the development of new numerical algorithms. Therefore let us consider specific features of the parallel implementation of the approach proposed. As opposed to Fourier transform, applying the Laguerre transform does not allow the calculation in parallel of the function for different , as th coefficient of series expansion (6) recurrently depends on the -th coefficient (9). Hence, it is needed to parallelize the algorithm at the stage of solving problem (10). In this case, the main difficulties are in solving SLAEs (17) and (18), whereas the parallel operation of multiplication by the matrices is not consuming from the communication standpoint. Let us consider the 2D data decomposition (Fig. 11a), where computing nodes contain processors having the shared memory access allowing implementation of in-node inter-process communications considerably faster than the communications among the nodes. For solving SLAEs (17),(18) with banded matrices as a parallel algorithm it is reasonable to use the parallel dichotomy algorithm, which was developed for tridiagonal matrices [44] and block-tridiagonal matrices [28]. With respect to the number of arithmetical operations, the dichotomy algorithm is comparable with other available algorithms; however the time needed for inter-process communications is considerably less in the dichotomy algorithm as compared to other algorithms. This is because the implementation of the dichotomy process on a supercomputer reduces to calculating the sum of series for distributed data. The commutative and associative properties of addition enable a considerable reduction in the total computation time with the use of inter-processor interaction optimization algorithms. First, the parallel dichotomy algorithm was developed for solving SLAEs with the same tridiagonal matrix but different right-hand sides. In [45], the dichotomy algorithm was applied to solving SLAEs with the Toeplitz tridiagonal matrices. It was shown that the Toeplitz tridiagonal matrices were able to effectively solve SLAEs both with one and several right-hand sides. In [27, 29], the dichotomy algorithm was applied to implement a spectral-difference method to calculate acoustic and elastic wave fields. The authors could effectively use from up to processors for one calculation and to obtain a highly accurate numerical solution of the dynamic problem of elasticity theory. Thus, the dichotomy algorithm is a powerful instrument for solving SLAEs with tridiagonal and block-tridiagonal (banded) matrices.

The parallel dichotomy algorithm makes possible to use thousands of processors, however its efficiency depends on the number of simultaneously solved equations and the size of a SLAE matrix. However, such dependence is considerably weaker as compared to other parallel algorithms for solving tridiagonal and block-tridiagonal SLAEs. Due to dependence of the solution on it is impossible to solve problems (17),(18) for all depths in parallel, which would allow reducing communication costs at the expense of striping the inter-processor exchanges. Therefore the parallelization for -direction should be done within one computational node profiting from the shared memory for the fast data exchange among processors thus allowing a considerable reduction of communication costs.

To increase the number of processors participating in one calculation it seems reasonable to carry out parallelization in -direction. In this case calculations can be effectively implemented by the conveyor principle, i.e. as soon as at the node number the functions are calculated, their values and those of their second derivatives at the lower boundary of a sub-domain are communicated to the computational node number . At the node with number , the calculation of functions with numbers starts, while at the node –functions with numbers . The value of the second derivative at the boundary of two sub-domains is required for setting boundary conditions when constructing the cubic interpolation spline within the Richardson extrapolation procedure. Taking into account the fact that the number of terms of series (6) essentially exceeds the number of computational nodes, the load of such computational conveyor will be optimal. The conveyor communication costs will be minor because only one inter-node exchange for calculation of the functions for one value is required.

3 Numerical experiments

Let us discuss a few tests that would allow the evaluation of the quality of the solution to be obtained as compared with other known algorithms. Numerical procedures were implemented in Fortran-90 using the MPI library. The calculation was performed on the ”Lomonosov”  supercomputer of the Moscow State University. The supercomputer comprises Intel Xeon X5570 fore-core processors operating at GHz in the Infiniband QDR communication environment. Each computational node contains two processors and 12 GB of RAM.

The coefficients of difference scheme (11) were chosen as follows[38]: , , , , , , , to provide the twelfth approximation order. The calculations have shown that as compared to classical difference schemes of high accuracy orders based on the Taylor expansion, the DRP approach makes possible to decrease the number of nodes in the difference scheme approximately by the factor of two with comparable accuracy. Taking into consideration the fact that computer costs for decomposition of matrices (17),(18) are proportional to the third degree of the matrix band width and forward and backward substitution of decomposition have costs proportional to the second degree, the effect of applying the DRP schemes is significant. Also, this allows one to decrease the demands for a minimum required random access memory for storing the matrices and .

3.1 Impulse Response

In the first test we illustrate analyzing the accuracy by impulse responses. For the calculation we used a homogeneous medium model with the speed and the size , the mesh steps being and . A point source was located at the center of the upper surface; the time dependence was as follows:


where . As compared to the Fourier transform, where the basis functions are uniquely defined, the two parameters, and , should be set for using the Laguerre transform (6). These parameters were experimentally chosen based on the analysis of the convergence rate of the Fourier–Laguerre series for the shifted function with , where is the upper boundary of the time interval for which the wave field is calculated. The parameters and are chosen such that the function with in the mean-quadratic norm is approximated accurate to . It should be regarded that as the value of grows, the values of the spatial derivatives of the functions will also increase. For this reason, the size of the mesh has to be simultaneously decreased in order to prevent growth of the error of the spatial approximation as the parameter grows. The number of addends in series (6) was for ; the expansion parameters were and .

Figures 1a,b present snapshots of the wave field from the point source for equation (4), which were obtained employing the finite difference algorithm (FD) of second accuracy order using the Marchuk-Strang splitting. This software procedure has been implemented in the available program package Seismic Unix. It is evident (Figs. 1a,b) that such an approach does not provide any accuracy of calculations due to multiple numerical noise and artifacts. In this case the transfer from the mesh with the step (Fig. 1a) to the mesh with the step (Fig. 1b) does not allow an increase in calculation accuracy up to the acceptable level.

A conceptually other situation arises (Figs. 1c,d,f,e) when the proposed Laguerre finite difference method (LFD) is used. It is clear that the new approach excludes the appearance of any numerical noise and the approximation error manifests as numerical dispersion which, as expected, is more distinct for the Crank-Nicolson scheme of second accuracy order (Figs. 1c,d). In Figs. 1e,f it is clear that applying the Richardson extrapolation makes possible to essentially improve the quality of solution. In addition, as for the second order scheme and for the fourth order of accuracy the transfer to a finer mesh allows diminishing the numerical dispersion. Thus, the key point in this case is the presence of convergence with decreasing the mesh step, which is not observed for calculations (Figs. 1a,b).

Figure 1: Snapshots for the wave field at for the homogeneous velocity model and impulse source . FD method, 2nd order (a)  and (b) , LFD method, 2nd order (c)  and (d) , LFD method, 4th order (e)  and (f) 

In spite of the absence of smoothness in the solution due to singularities [14, 15], the use of the method of fourth order of accuracy allows a decrease in numerical dispersion as compared to the scheme of second order of accuracy. However one should take into account the fact that with an equal number of mesh points the algorithm based on the Richardson extrapolation requires three times as much arithmetical operations as the Crank-Nicolson scheme. This is explained by the necessity of solving equations (10) on two meshes with the steps and . The question arises whether it is possible just to increase the number of mesh nodes by the factor of three and to make use of the second order scheme to obtain a result compatible in accuracy with the fourth order scheme. Fig. 2 shows the dependence of a wave field on the coordinate along the straight line ”Slice” (Fig. 1f). As is obvious, when the mesh step for the method of second accuracy order is three times smaller than that for the fourth order method (computer costs being approximately the same) the second order method yields less accurate results. Thus, it is reasonable to apply more complicated procedures based on the Richardson extrapolation, requiring large computer costs. It is also evident that the algorithm which is based on the local extrapolation inserts additional dispersion errors but better preserves amplitudes. Nevertheless, we are not going to use it because it is unstable for inhomogeneous media and inserts distinct artifacts associated with numerical dispersion into the image of a wave field.

Figure 2: Dependence of the wave field on the coordinate along the straight line ”Slice” (Fig. 1) for different meshes for the LFD method.

When implementing the Richardson extrapolation algorithm we considered selecting the interpolation procedure and its influence on the solution to be solved. When using the linear interpolation the LFD algorithm demonstrated a very strong anisotropic dissipation in -direction such that when passing a distance of several wavelengths the wave amplitude was tending to zero. The use of the parabolic or the cubic Lagrange interpolation somewhat improved the situation, but preservation of the amplitude was worse than for the cubic spline interpolation. Using the Lagrange interpolation of high orders as well as using splines of the orders exceeding the third caused instability of calculation. Thus, the choice of the cubic interpolation spline can be considered to be appropriate as other interpolation algorithms are either unstable or essentially dissipative.

Although the stability of the algorithm proposed has not been proved yet due to the difficulty in making the analytical analysis. In view of the following: asymmetry of matrices, the presence of singularities in a solution, applying the Richardson extrapolation procedure spline interpolation, etc. the algorithm in question in numerous experiments for different velocity models appeared to be more stable provided and a local contrast of the velocity medium model is not quite distinct. A higher level of stability as compared to other algorithms is explained by the fact that approximation errors in the direction of the wave field extrapolation are dominantly of dissipative character.

3.2 Migration procedure

To form the image of the Earth’s interior by the recorded seismic data one can use different algorithms of seismic migration [2, 11, 46, 47] including those of the wave field downward continuation algorithms. Let us consider the post-stack migration based on the model of explosive boundaries [1]. This model allows the evaluation of the algorithm proposed as compared to other known methods. As many studies have dealt with the problem in question, we will focus on its computational features. As there is no analytical solution, we will make a comparison by the FD method, with the numerical-analytical Fourier Finite Difference (FFD) method [23] and the Phase-Shift Plus Interpolation method (PSPI) [22] that are widespread in computational seismic prospecting. The two latter algorithms were developed instead of the FD method because it does not provide the required accuracy of calculation, and what is more important the computational stability when solving the one-way equation. However the FFD and PSPI methods approximately calculate the solution to the one-way equation for an inhomogeneous velocity model. The reason is that for the PSPI method the solution is obtained as a combination of several solutions for inhomogeneous models with different velocities. In the context of the FFD method, the input operator is represented as sum of two operators, one of them corresponding to a homogeneous medium with some averaged in horizontal velocity (usually, a minimum velocity is used as the reference velocity). The other operator represents a varying velocity correction. In the context of the Marchuk-Strang splitting procedure the first operator is inversed by analytical methods and the solution for correction is defined by finite difference methods of second order of accuracy. Although such a change is not similar to the original problem in terms of mathematics and is unstable for strongly contrast media, such an approach makes possible to considerably decrease computer costs. On the other hand, impossibility of evaluating the level of numerical errors on a sequence of imbedded meshes is an essential disadvantage of the FFD algorithm.

3.2.1 Syncline model

Theoretical seismograms (Fig. 3b) for the syncline model (Fig. 3a) were obtained with the help of the Gaussian beams algorithm [48, 49] implemented in the package Seismic Unix. For setting the boundary condition on the upper surface, the function for zero-offset section . Was expanded in series (6) with the parameters , and for . The calculations were carried out on meshes with the steps , and . According to the model of explosive boundaries the calculation velocities were set to be half the true velocity of the medium model.

As well as in the calculation of the wave field from the point source (Fig. 1) inadmissible level of numerical noise in the context of the migration problem for the FD method is observed (Fig. 4). In the FFD method, the noise level is much lower (Fig. 5), however, there are numerical artifacts and, in addition, the second horizon is insufficiently focused. When turning from the mesh with the step (Fig. 5a) to the mesh with the step (Fig. 5b) manifestation of artifacts remained as previously. Thus, the FFD algorithm can generate artifacts but does not allow their a posteriori evaluation on a sequence of imbedded meshes.

A conceptually different situation is observed with the LFD algorithm (Fig. 6), where the absence of high-frequency noise in the images obtained is distinct. The finite difference approximation errors for -direction manifest as effect of numerical dispersion, whereas for -direction the errors are of dissipative nature. The influence of these errors in turning to a finer mesh rapidly decreases, and the solution convergence of the mesh step makes possible to separate approximation errors from those of input data and inaccuracy of a model as it is. Hence, analyzing the results obtained when solving practical tasks can be essentially simplified. Additionally, calculations for a finer mesh with the step for all the three algorithms were carried out. The FD method for a small mesh step appeared to be unstable, whereas the quality of the image for the FFD method remained unchanged because the level of noise and artifacts was comparable with the calculation results with the mesh steps and . For the LFD method, the image in Fig. 6d does not visually differ from that obtained with the mesh step (Fig. 6c). Thus, with allowance for the convergence of the LFD algorithm, the result obtained with the mesh step can be considered to be sufficiently accurate.

Figure 3: (a) Syncline model and (b) zero-offset section.
Figure 4: Snapshots for the wave field at for the model and the zero-offset section in Fig. 3 for the FD method (a)  and (b) 
Figure 5: Snapshots for the wave field at for the model and the zero-offset section in Fig. 3 for the FFD method (a)  and (b) . Computational artifacts are indicated by arrows.
Figure 6: Snapshots for the wave field at for the model and the zero-offset section in Fig. 3 for the LFD method (a) , (b) , (c) , (d) 

3.2.2 Sigsbee Model

For Sigsbee2Amodel [50] (Fig. 7a), theoretical seismograms (Fig. 7b) were calculated using the algorithm of explosive boundaries implemented in the Madagascar package [51]. For setting the boundary condition, the function for zero-offset section was expanded in series (6) for , parameters and for . The calculations were done on the meshes with the steps , and .

According to the results of computational experiments, the FD method appeared to be unstable for the given velocity model and selected numerical parameters. Different degrees of smoothing the velocity medium model did not allow attaining the calculation stability for this method, therefore the FD method within this test is excluded from further consideration.

The PSPI method from the Seismic Unix package has shown unsatisfactory results because of bad energy focusing (Fig. 8a). This method was unable to correctly map boundaries of the complex salt shape (Fig. 8b). However note that the point type diffractors underlying the salt inclusion are correctly focused. The FFD method as compared to the PSPI has allowed obtaining a clearer image of the salt inclusion (Fig. 9b), but the image additionally contains numerous artifacts (Fig. 9a) due to approximate character of calculation formulas of the FFD method. The fact is, as well as in the previous test, an increase in resolving power of the mesh does not bring about an increase in calculation accuracy.

As opposed to the FFD and the PSPI methods, a conceptually different situation arises with the use of the LFD method: as is seen in Fig. 10, the LFD method does not contain additional artifacts but only those which are simultaneously present for the FFD and the PSPI algorithms, i.e. the artifacts being present in input data of the problem (zero-offset section). Using the DRP approach to determine coefficients of the finite difference approximation makes possible to attain high accuracy of calculation and to decrease numerical dispersion for -direction. For -direction, typical is numerical dissipation (Fig. 10a,b) that brings about smoothing an image. However, the numerical dissipation can be decreased at the cost of choosing a smaller step (Fig. 10c,d). Fourth approximation order due to the Richardson extrapolation allows the use of admissible steps for attaining the required calculation accuracy. Thus, errors caused by the LFD algorithm are controlled by choosing the step of a spatial mesh, whereas the choice of computational parameters for the PSPI and the FFD methods is a more ambiguous task, and decreasing a spatial mesh step for these methods does not really affect the imagery quality.

Figure 7: (a) The 2D Sigsbee2A salt model and (b) zero-offset section.
Figure 8: Snapshots for the wave field at for the model and the zero-offset section in Fig. 7 for the PSPI method for . Computational artifacts are indicated by arrows.
Figure 9: Snapshots for the wave field at for the model and the zero-offset section in Fig. 7 for the FFD method for . Computational artifacts are indicated by arrows.
Figure 10: Snapshots for the wave field at for the model and the zero-offset section in Fig. 7 for the LFD method (a),(b)  and (c),(d) . Computational artifacts are indicated by arrows.

Thus, in the context of the post-stack migration using the method in question allows obtaining a more accurate image of the Earth’s interior as compared to widespread mesh algorithms of solving the one-way equation. The method proposed can be also implemented for the pre-stack migration with a view to obtain high quality images. As opposed to numerical-analytical approaches requiring fulfillment of Fourier transform both with respect to time and with respect to space, the use of purely difference spatial approximation will enable us in the future to carry out calculations, first, for non-uniform meshes and, second, to implement a curvilinear boundary for the upper surface. The convergence of the new algorithm of the mesh step makes possible to evaluate the solution accuracy on a sequence of imbedded meshes thus separating approximation errors from errors of setting input data of the problem. For most of numerical migration methods it is necessary to preliminarily smooth a velocity medium model in order to provide the calculation stability. This brings about errors in the wave field kinematics and presents a severe problem in choosing a smoothing procedure. As the computational experiments show, the proposed LFD method for considered models and chosen computational parameters does not require preliminary smoothing of the medium model function, which counts in favor of a higher degree of stability as compared to existing mesh methods.

3.3 Performance consideration

As a matter of fact, a higher accuracy is not attained without supplementary computer costs. The single-processor implementation of the method proposed has shown that the LFD algorithm requires from to times more computational time than the PSPI, the FD, the FFD algorithms. The costs of the LFD algorithm for solving the SLAEs with banded matrices are essential, but the possibility of using optimized library functions for solving the SLAEs allows an increase in the method performance on the whole. Higher computer costs of the LFD algorithm can be justified because the FD method does not provide any accuracy and is weakly stable, and the FFD and the PSPI methods solve a certain approximate problem but not the one-way equation thus generating multiple noise and artifacts of the wave field.

Figure 11: (a) Data decomposition among processors. (b) Speedup versus the number of processors for different meshes and velocity models.

For testing the performance of the method proposed the calculation time for a different number of mesh nodes and the number of processors was assessed. The use of multi-processor systems makes absolute time costs insignificant because with a maximum number of processors the calculation time for all the test problems does not exceed several seconds (Table 1).

num proc.
12 247 1008 44 119 944
24 122 446 21 62 472
48 76 294 14.8 41 313
216 11 46 2.4 6.2 49.6
228 10 43 2.3 5.9 47.3
240 9.7 40 2.2 5.6 44.9
Table 1: Calculation time (seconds) versus the number of processors for meshes of different resolution for all the tests.

The dependence of the speedup value on the number of processors for all the tests is presented in Fig. 11b, indicating to the linear character of the dependence, for the first test attaining the super-linear speedup. A similar effect of super-linear speedup was observed when using the parallel dichotomy algorithm within the implementation of the alternating implicit direction method for solving the 2D Poisson equation[44]. Thus, the dichotomy algorithm implemented in the context of the approach proposed for solving the 2D one-way equation allows the effective use of a considerable number of processors. We can conclude that the chosen data decomposition of the problem (Fig. 11b) is fairly appropriate, as multiple exchanges for x-direction are instantly made at the expense of using the shared memory within one computational node, while slow inter-node communicative interactions for -direction for solving problem (17) are done only once.

4 Conclusion

One of the major results of this research is in managing to construct a stable, efficient, purely finite difference method of high accuracy order for solving the one-way equation. By now this has presented considerable difficulties in numerical analysis due to the necessity of providing high accuracy and at the same time the calculation stability. Replacing Fourier transform by the Laguerre transform after employing the finite difference approximation of spatial derivatives, using algebraic transformations made possible to obtain the banded SLAEs appropriate for solutions by direct methods, in particular, by the parallel dichotomy algorithm. The rejection from using the splitting of the Marchuk-Strang type for the problem operator has excluded appearance of numerous computing artifacts and has increased both the order of approximation and the stability of the method as a whole. The one-way equation does not allow a correct description of the wave field amplitude for inhomogeneous media, hence the optimization of values of finite difference scheme coefficients aimed at decreasing the phase errors is a reasonable approach. The numerical experiments have shown that such DRP schemes make possible to decrease the mesh step in horizontal approximately by the factor of two as compared to classical difference schemes based on the Taylor series. Implementation of the Richardson extrapolation procedure has not only increased the accuracy of the algorithm up to fourth order relative to the mesh step in depth, but also has provided the stability of the method as compared to using schemes of the Adams-Moulton type. Applying implicit Runge-Kutta methods of high orders of accuracy is unreasonable, as in this case it is required to solve SLAEs of much higher dimensions. In addition, as well as for the Richardson extrapolation, the interpolation of the wave field values for auxiliary mesh nodes will be needed.

The algorithm proposed can be generated to solving 3D problems as well. To this end, the efficient iterative procedure for solving SLAEs should be developed. This is because the parallel dichotomy algorithm, being a direct method, is not intended for banded matrices of a large width. Also, it is necessary to reject not only the splitting of the Marchuk-Strang type for the problem operator, but also the splitting along horizontal directions that is used for decreasing computer costs, but generates numerical anisotropy. The efficient algorithm for solving the 3D one-way equation will be presented in the next publication in the near future.

5 Acknowledgments.

The author is grateful to Gennady Zhernyak for fruitful discussions of the paper.


  • [1] J.F. Claerbout. Toward a unified theory of reflector mapping. Geophysics, 36(3):467–481, 1971.
  • [2] O. Yilmaz and S.M. Doherty. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data. Society of Exploration Geophysicist, 2001.
  • [3] D.A. Angus. The one-way wave equation: A full-waveform tool for modeling seismic body wave phenomena. Surveys in Geophysics, 35(2):359–393, 2013.
  • [4] F.D. Tappert. Wave Propagation and Underwater Acoustics, chapter The parabolic approximation method, pages 224–287. Springer Berlin Heidelberg, Berlin, Heidelberg, 1977.
  • [5] F.B. Jensen, W.A. Kuperman, M.B. Porter, and H. Schmidt. Computational Ocean Acoustics. Springer-Verlag New York, 2011.
  • [6] D. Lee, A.D. Pierce, and Er-C. Shang. Parabolic equation development in the twentieth century. Journal of Computational Acoustics, 08(04):527–637, 2000.
  • [7] E.L. Lindman. ”Free-space” boundary vonditions for the time dependent wave equation. Journal of Computational Physics, 18(1):66 – 78, 1975.
  • [8] B. Engquist and A. Majda. Absorbing boundary conditions for the numerical simulation of waves. Math. Comp., 31:629–651, 1977.
  • [9] L.N. Trefethen and L. Halpern. Well-posedness of one-way wave equations and absorbing boundary conditions. Mathematics of Computation, 47(176):421–435, 1986.
  • [10] M. Leontovich and V. Fock. Solution of the problem of propagation of electromagnetic waves along the earth’s surface by the method of parabolic equation. Acad. Sci. USSR J. Phys., 10:13–24, 1946. (Engl. transl).
  • [11] J.F. Claerbout. Imaging the Earth’s Interior. Blackwell Scientific Publications, Inc., Cambridge, MA, USA, 1985.
  • [12] D. Lee and G. Botseas. Ifd: An implicit finite-difference computer model for solving the parabolic equation. Technical Report 6659, Nav. Underwater Syst. Ctr, 1982.
  • [13] L. Halpern and L.N. Trefethen. Wide-angle one-way wave equations. J. Acoust. Soc. Am., 84(4):1397–404, 1988.
  • [14] A. Bamberger, B. Engquist, L. Halpern, and P. Joly. Higher order paraxial wave equation approximations in heterogeneous media. SIAM Journal on Applied Mathematics, 48(1):129–154, 1988.
  • [15] A. Bamberger, B. Engquist, L. Halpern, and P. Joly. Parabolic wave equation approximations in heterogenous media. SIAM Journal on Applied Mathematics, 48(1):99–128, 1988.
  • [16] G.I. Marchuk. Some application of splitting-up methods to the solution of mathematical physics problems. Applik Mat, 13(2):103–132, 1968.
  • [17] G. Strang. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis, 5(3):506–517, 1968.
  • [18] G.I. Marchuk. Handbook of numerical analysis. In P. G. Ciarlet and J. L. Lions, editors, Handbook of Numerical Analysis, volume 1, pages 197–460. Elsevior Science Publishers BV, 1990.
  • [19] C. Bunks. Effective filtering of artifacts for implicit finite-difference paraxial wave equation migration1. Geophysical Prospecting, 43(2):203–220, 1995.
  • [20] T. Fei and K. Larner. Elimination of numerical dispersion in finite-difference modeling and migration by flux-corrected transport. Geophysics, 60(6):1830–1842, 1995.
  • [21] J. Gazdag. Wave equation migration with the phase-shift method. Geophysics, 43(7):1342–1351, 1978.
  • [22] J. Gazdag and P. Sguazzero. Migration of seismic data by phase shift plus interpolation. Geophysics, 49(2):124–131, 1984.
  • [23] D. Ristow and T. Rühl. Fourier finite-difference migration. Geophysics, 59(12):1882–1893, 1994.
  • [24] M.N. Guddati. Arbitrarily wide-angle wave equations for complex media. Computer Methods in Applied Mechanics and Engineering, 195:65 – 93, 2006.
  • [25] M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover, New York, ninth dover printing, tenth gpo printing edition, 1964.
  • [26] B.G. Mikhailenko. Simulation of seismic wave propagation in heterogeneous media. Siberian J. of Numer. Mathematics, 6:415–429, 2003.
  • [27] A.G. Fatyanov and A.V. Terekhov. High-performance modeling acoustic and elastic waves using the parallel dichotomy algorithm. J. Comp. Phys., 230(5):1992–2003, 2011.
  • [28] A.V. Terekhov. A fast parallel algorithm for solving block-tridiagonal systems of linear equations including the domain decomposition method. Parallel Comput., 39(6-7):245–258, 2013.
  • [29] A.V. Terekhov. Spectral-difference parallel algorithm for the seismic forward modeling in the presence of complex topography. Journal of Applied Geophysics, 115(0):206–219, 2015.
  • [30] Myung W. Lee and Sang Y. Suh. Optimization of one-way wave equations. Geophysics, 50(10):1634–1637, 1985.
  • [31] L. Debnath and D. Bhatta. Integral Transforms and Their Applications, Second Edition. Taylor & Francis, 2006.
  • [32] B.G. Mikhailenko. Spectral laguerre method for the approximate solution of time dependent problems. Applied Mathematics Letters, 12:105–110, 1999.
  • [33] J. Crank and P. Nicolson. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Proc. Camb. Phil. Soc., 43(1):50–67, 1947.
  • [34] Y. Zhang, G. Zhang, and N. Bleistein. Theory of true-amplitude one-way wave equations and true-amplitude common-shot migration. Geophysics, 70(4):E1–E10, 2005.
  • [35] C.K.W. Tam and J.C. Webb. Dispersion-relation-preserving finite difference schemes for computational acoustics. Journal of Computational Physics, 107(2):262 – 281, 1993.
  • [36] W. Liang, Y. Wang, and C. Yang. Determining finite difference weights for the acoustic wave equation by a new dispersion-relationship-preserving method. Geophysical Prospecting, 63(1):11–22, 2015.
  • [37] C. Chu and P.L. Stoffa. Determination of finite-difference weights using scaled binomial windows. Geophysics, 77(3):W17–W26, 2012.
  • [38] J.-H. Zhang and Z.-X. Yao. Optimized explicit finite-difference schemes for spatial derivatives using maximum norm. Journal of Computational Physics, 250:511 – 526, 2013.
  • [39] A.A. Samarskii. The Theory of Difference Schemes. Marcel Dekker, 2001.
  • [40] R.W. Cottle. Manifestations of the schur complement. Linear Algebra and its Applications, 8(3):189 – 211, 1974.
  • [41] H.V. Henderson and S.R. Searle. On deriving the inverse of a sum of matrices. SIAM Review, 23(1):53–60, 1981.
  • [42] L.F. Richardson. On the approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 83(563):335–336, 1910.
  • [43] G. I. Marchuk and V. V. Shaidurov. Difference methods and their extrapolations. Appl. Math. Springer, New York, NY, 1983.
  • [44] A.V. Terekhov. Parallel dichotomy algorithm for solving tridiagonal system of linear equations with multiple right-hand sides. Parallel Comput., 36(8):423–438, 2010.
  • [45] A.V. Terekhov. A highly scalable parallel algorithm for solving toeplitz tridiagonal systems of linear equations. Journal of Parallel and Distributed Computing, 87:102–108, 2016.
  • [46] S.H. Gray, J. Etgen, J. Dellinger, and D. Whitmore. Seismic migration problems and solutions. Geophysics, 66(5):1622–1640, 2001.
  • [47] B. Biondi. 3D Seismic Imaging. Society of Exploration Geophysicists, 2006.
  • [48] M.M. Popov. A new method of computation of wave fields using gaussian beams. Wave Motion, 4, 1982.
  • [49] V. Cerveny. Gaussian beam synthetic seismograms. J. Geophys., 58:44–72, 1985.
  • [50] J. Paffenholz, B. McLain, J. Zaske, and P.J. Keliher. Subsalt multiple attenuation and imaging: Observations from the Sigsbee2B synthetic dataset, chapter 538, pages 2122–2125.
  • [51] S. Fomel, P. Sava, I. Vlad, Y. Liu, and V. Bashkardin. Madagascar: open-source software project for multidimensional data analysis and reproducible computational experiments. Journal of Open Research Software 1:e8, 2013.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description