- Email: [email protected]
Improved amplitude multi-one-way modeling method
Wave Motion 43 (2005) 99–115
Improved amplitude multi-one-way modeling method D. Kiyashchenkoa,∗ , R.-E. Plessixb , B. Kashtana , V. Troyana a
Saint-Petersburg State University, Laboratory of Elastic Media Dynamics, Ulianovskaya Street 1, 198504 St. Petersburg, Russian Federation b Shell International E&P, P.O. Box 60, 2280 AB Rijswijk, The Netherlands Received 29 September 2004; received in revised form 3 May 2005; accepted 9 June 2005 Available online 2 September 2005
Abstract Classical paraxial wave equations, also called one-way wave equations, lead to efficient numerical schemes in threedimensional spaces for computing the wave field. The kinematics is well approximated by these algorithms, however the dynamics is generally not correctly modeled when the medium is inhomogeneous. This prevents from using one-way wave equation scheme for solution of inverse problems. In order to improve estimation of the amplitudes, a new multi-one-way scheme is introduced. This approach is based on an iterative solution of the factorized two-way wave equation with a right-hand side incorporating the information about medium heterogeneities. The numerical scheme is a succession of several classic oneway schemes. The proposed approach takes into account vertical and horizontal velocity variations of the medium and allows modeling reflected waves. Numerical results for several two-dimensional examples with smooth and rough velocity contrasts are compared with finite-difference solutions of the two-way wave equation. It is shown that the multi-one-way scheme significantly reduces the error in the amplitude estimates of the wave field. © 2005 Elsevier B.V. All rights reserved. PACS: 43.20.Bi Keywords: One-way equation; True amplitude
1. Introduction Wide range of inverse problems require accurate solutions of the wave equation, for example in geophysics or underwater acoustics. The finite-difference solution of the full wave equation still remains expensive in threedimensional spaces and is not yet used for practical applications [1]. The high-frequency (ray-based) solutions ∗
Corresponding author. Tel.: +7 8124284317; fax: +7 8124284317. E-mail addresses: [email protected] (D. Kiyashchenko), [email protected] (R.-E. Plessix).
0165-2125/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.wavemoti.2005.06.003
100
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
suffer from some limitations over complex media. The paraxial approximation of the wave equation then gives an alternative approach [2]. This approach consists of particularizing one direction and extrapolating the wave field step by step along this direction. This approach is also called classic one-way because the wave field contains only the up-going (or down-going) waves, the reflected waves propagating in the direction opposite to preferred one are not modeled. In contrast, the full wave equation is called two-way wave equation in this paper. One of the drawbacks of the classic one-way wave equation approximation is the large errors in the amplitude estimates of the wave field. For example, in exploration geophysics these errors generally prevent us from using the classic one-way scheme in preserved amplitude migration algorithms. Migration is the technique which is commonly used in seismic exploration for mapping the time-domain data (seismograms) to the depth domain. This gives a structural image of the subsurface. The preserved amplitude migration allows to relate the amplitudes of the imaged reflectors with media parameters. Generally, a migration algorithm requires the computation of the incident down-going wave field propagating from the source point and the back-propagated down-going wave field of the seismic data [3]. A preserved amplitude migration requires a good amplitude estimates of those wave fields. The idea of the one-way approach can be found in pioneer works of Leontovich and Fock [4] and Tappert [5], who considered the solution of the wave equation in a preferred propagation direction. Claerbout [6,3] introduced this approach in exploration geophysics. One-way wave equations are based on an approximation of the two-way dispersion curve. This dispersion curve is obtained from the time–space domain wave equation by temporal and spatial Fourier transforms and represents the wave equation in the frequency–wavenumber domain. This derivation assumes a constant velocity field. The wavenumber in the preferred direction is expressed in terms of wavenumbers in the other directions. The one-way wave equation in the spatial coordinates is then formally derived by replacing the wavenumbers by the spatial derivatives in the dispersion curve. This one-way equation involves the non-local square root operator. In order to obtain an efficient numerical scheme for spatially varying background velocity and to compute the wave field step by step along the preferred direction with a marching approach, a first possibility is to approximate the square root operator with local operators, for instance with Pad´e fractions leading to the classic 15◦ , 45◦ or 60◦ approximation scheme (see [3,7–9]). A second possibility is to first extrapolate the wave field in the frequency–wavenumber domain assuming a constant velocity model and then to apply a correction using Born approximation. This corresponds to the phase-shift methods [10,11] and the screen propagator methods [12–16]. Up to a given angle from the preferred direction, the one-way wave equation methods correctly model the kinematics of waves. Unfortunately, the dynamic properties (amplitudes) are not accurately estimated with most of the one-way schemes. Zhang et al. [17] have shown that the classic one-way eikonal equation give the same traveltimes as the two-way eikonal equation, but it is not the case for the transport equations, which govern the amplitudes in high-frequency approximation. In the literature, few papers have addressed the problem of the amplitude behavior of the one-way scheme with spatially varying velocity. Bamberger et al. [18] suggested a method to take into account the transmission losses of the normal incident waves through the artificial wavefield discontinuities introduced during the marching approach. Zhang et al. [17] proposed to add a term to the one-way wave equation to make the eikonal and transport equations of the latter matching those obtained from the two-way wave equations. Interesting formulations of one-way and two-way wave propagation problems have also been suggested by Wapenaar and Grimbergen [19] and De Hoop and Gautesen [20]. The approach is based on a decomposition of the acoustic wave equation into two coupled one-way equations for upward and downward propagating waves. When the medium is homogeneous, these equations are decoupled. In the case of not too large velocity gradients, the usual one-way approximation is used. This approach allows an iterative solution to take into account wave transmission and reflection phenomena. However, for three-dimensional applications the expressions for the square root operator proposed in these papers may lead to very expensive computational procedure. In this paper, a different approach is derived to correct for the amplitude errors of the classic one-way wave equation due to velocity variations. The two-way wave equation is rewritten as the product of two one-way wave equations plus an error (correction) term. An iterative solution is then conducted using the perturbation theory. The scheme consists of a number of classic one-way wave equations with right-hand sides. The right-hand sides take into account the error term. This approach corrects for the errors due to the velocity variations and models the primary
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
101
reflections with rough velocity models. This method is called multi-one-way method. The use of the perturbation theory assumes small velocity contrasts. In practical applications, three one-way equations are solved and the error term is computed. This means that this method is about five times more expensive than the classic one-way method. The outline of the paper is following. First the derivation of the classic one-way wave equation is recalled and the multi-one-way method is explained. Then three two-dimensional examples are described to demonstrate the relevance of the approach and its limitations with smooth and rough velocity models. In Appendix A, the modeling of the transmission and reflection coefficients with the multi-one-way method is analyzed for the case of a plane interface.
2. Two-way, one-way and multi-one-way modeling methods 2.1. The one-way approach One-way wave equation is a paraxial approximation of the full two-way wave equation after the choice of a preferred direction of the wave propagation. In this work, the preferred direction is the vertical axis (z-direction). The constant density two-way acoustic wave equation in the space-frequency domain reads 2 ∂ ω2 (1) + 2 + x u = f (x, z, ω), ∂z2 c ∂ ∂ ∂ with x = ∂x 2 in two dimensions and x = ∂x2 + ∂y2 in three dimensions. c(x, z) is the velocity, u(x, z) the pressure field and f (x, z, ω) the source term. The vector x has two components, x ≡ (x, y), in three dimensions and one component, x ≡ x, in two dimensions. The derivation of the one-way wave equation is based on the factorization of the two-way operator 2 ∂ 2 (2) + Λ A= ∂z2 2
2
2
into the product of up-going and down-going wave one-way operators: ∂ ∂ L= + iΛ − iΛ , (3) ∂z ∂z ω . where Λ(x, z) = k2 + x with k(x, z) = c(x,z) L equals A only when the velocity does not depend on the depth coordinate z. Otherwise the square root operator, ∂ Λ, depends on z, and does not commute with the operator ∂z . In fact L satisfies ∂ ∂Λ ∂ + iΛ − iΛ = A − i . (4) L= ∂z ∂z ∂z In practice, a local approximation, Λa , of the square root operator, Λ, is used. The product of the up-going and down-going one-way wave operators, with the local approximation of the square root, gives ∂ ∂ ∂Λa La = + iΛa − iΛa = A + (Λ2a − Λ2 ) − i . (5) ∂z ∂z ∂z With the error operator, E, defined by E = −i
∂Λa + (Λ2a − Λ2 ), ∂z
(6)
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
102
this gives La − E = A.
(7)
The classic one-way approach neglects the error term, E. Assuming no source term for the time being (f (x, z, ω) = 0), this leads to ∂ ∂ (8) + iΛa − iΛa = 0 ∂z ∂z or if the derivation is repeated by interchanging the up- and down-going one-way wave operators, then ∂ ∂ − iΛa + iΛa = 0. ∂z ∂z
(9)
From Eqs. (8) and (9) the classic one-way equations are obtained: ∂u = iΛa u, ∂z
(10)
with = 1 for the classic down-going one-way wave equation and = −1 for the up-going one. These one-way approximations do not correctly estimate the amplitudes as shown by Zhang et al. [17]. In their study, they demonstrate that the travel times obtained from the one-way wave equation satisfy the two-way eikonal equation, but the amplitudes of the one-way solution do not satisfy the two-way transport equation. This proves that even in the high-frequency sense, the amplitudes are not correctly modeled with the classic one-way method in inhomogeneous media. The error term, E, explains the differences between the one-way and the two-way solutions. The expression a (6) has two parts. The first one, −i ∂Λ ∂z , is non-zero if vertical velocity variations are present. The second one, (Λ2a − Λ2 ), is non-zero for large incidence angles, i.e. large angles between the vertical (the z-direction) and the wave propagation directions, due to the square root approximation, and in the presence of lateral velocity variations. 2.2. Multi-one-way modeling method In order to obtain a better amplitude estimation of the wave field with an approach based on the one-way wave equation, the (two-way) solution, u, of La u = Eu
(11)
is approximated using the perturbation method. The solution u0 of the down-going one-way wave equation ∂u0 = iΛa u0 ∂z can be considered as a “zero-order” approximation of the solution of Eq. (11). A “first-order” approximation, u0 + u1 , is obtained when the correction term u1 satisfies ∂ ∂ + iΛa − iΛa u1 = Eu0 . ∂z ∂z
(12)
(13)
Eq. (13) can be solved in two steps. First the solution, u1/2 , of the up-going one-way wave equation with the right-hand side (source term) equal to the error term evaluated at u0 is computed: ∂ (14) + iΛa u1/2 = Eu0 . ∂z
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
103
Secondly the solution u1 of the down-going one-way wave equation with u1/2 as the right-hand side is computed: ∂ − iΛa u1 = u1/2 . (15) ∂z The set of the one-way wave equations, Eqs. (12), (14) and (15), should be solved sequentially in order to obtain a “first-order” approximation of the solution of the two-way wave equation. Therefore, this approach is called multione-way modeling method. Appropriate initial boundary conditions for the above mentioned equations should be defined. 2.3. Boundary conditions for the one-way equations in multi-one-way scheme In the previous subsections the study of the one-way wave equations was performed in the spatial domain with f (x, z, ω) = 0 in Eq. (1). To take into account the source term it is sufficient to consider the case f (x, z, ω) = δ(z)fs (x), with δ(z) a dirac delta-function in z-coordinate. The general solution can be obtained with the superposition principle. With the source term, Eq. (11) becomes La u = Eu + δ(z)fs (x). Neglecting the error term leads to the equation for the “zero-order” approximation, u0 : ∂ ∂ + iΛa − iΛa u0 = δ(z)fs (x). ∂z ∂z
(16)
(17)
Theoretically Eq. (17) is solved in two steps. The first step consists of solving the up-going one-way wave equation: ∂ (18) + iΛa g = δ(z)fs (x). ∂z The second step consists of solving the down-going one-way wave equation: ∂ − iΛa u0 = g. ∂z
(19)
The radiation free condition at z = zmax , the bottom of the computational domain, gives g(zmax ) = 0, which is the initial condition of the up-going Eq. (18). Therefore g(z) = 0 for all z > 0 for a source term f which is non-zero only at z = 0, the top of computational domain. In this way only Eq. (19) needs to be solved with the appropriate initial boundary condition. So, the first step of the multi-one-way scheme consists only in solving a down-going one-way wave equation. To derive the initial condition of Eq. (12) (or Eq. (19)), the first equation of the multi-one-way scheme, we follow [7]. As the right-hand side of Eq. (17) is proportional to δ(z)fs (x) and the left-hand side contains a second derivative in z, the first order derivative of u0 at z = 0 is discontinuous and satisfies: 1 ∂u0 = ± fs (x). (20) ∂z z=±0 2 With Eq. (12) this gives iΛa u0 |z=+0 = 21 fs (x),
(21)
or u0 |z=+0 =
1 Λa Λ−2 a fs (x). 2i
(22)
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
104
−2 Λ−2 a fs (x) can be approximated by Λ fs (x) which is the solution of a two-dimensional two-way wave equation in three dimensions (one-dimensional two-way wave equation in two dimensions). The initial boundary condition for Eq. (12) is then
u0 |z=+0 =
1 Λa h(x), 2i
(23)
2 where h(x) is the solution of x + c2 ω(x,0) h(x) = fs (x) with absorbing boundary conditions. The second equation of the multi-one-way approach, Eq. (14), is an up-going one-way wave equation. The initial boundary condition is deduced from the radiation condition at the bottom of computational domain (z = zmax ). This condition means that the wave field u1 should not contain up-going waves at the bottom of computational domain. Then u1/2 |z=zmax satisfies ∂ u1/2 |z=zmax = = 0. (24) − iΛa u1 ∂z z=zmax The third equation of the multi-one-way approach, Eq. (15), is a down-going one-way wave equation. The initial boundary condition is obtained assuming that the wave field u1 does not contain reflections from the top of the computational domain at z = 0. This means that all plane waves contained in the wave field u1 should propagate upward at z = 0. This gives ∂ =0 (25) + iΛa u1 ∂z z=+0 and with Eq. (15) − 2iΛa u1 |z=+0 = u1/2 |z=+0 .
(26)
The initial boundary condition u1 |z=0 is then computed as follows: 1 (27) Λ a h1 , 2i 2 where h1 is the solution of x + c2 ω(x,0) h1 = u1/2 |z=+0 . The initial boundary conditions for the multi-one-way set of Eqs. (12), (14), and (15) are then given by Eqs. (23), (24), and (27). u1 |z=+0 = −
2.4. Solution of one-way equations with non-zero source term The multi-one-way method consists of successive solutions of three one-way wave equations. The last two are inhomogeneous, as their right-hand side (source term) is non-zero. It is interesting to discuss in more details the solution of these equations. The second one-way wave equation of the multi-one-way approach is an up-going one-way wave equation with Eu0 as source term: ∂ (28) + iΛa u1/2 = Eu0 , ∂z with the initial boundary condition u1/2 (zmax ) = 0. We define for all ζ in [0, zmax ] the function e0 by e0 [ζ](x, z) = (Eu0 )(x, ζ)δ(z − ζ)
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
and v1/2 [ζ] by ∂ + iΛa v1/2 [ζ] = e0 [ζ], ∂z
105
(29)
with v1/2 [ζ](z) = 0 for z > ζ. At z = ζ an initial boundary condition is needed. v1/2 [ζ](z) at z = ζ is discontinuous due to the presence of the delta-function in the right-hand side of Eq. (29). The integration of Eq. (29) over a small interval (ζ − ε, ζ + ε) gives z=ζ+ε dz iΛa v1/2 [ζ](z) = Eu0 (x, ζ). (30) v1/2 [ζ](ζ + ε) − v1/2 [ζ](ζ − ε) + z=ζ−ε
The initial condition is obtained when ε goes to 0: v1/2 [ζ]|z=ζ−0 = −(Eu0 )(x, ζ)
(31)
and for z < ζ the function v1/2 [ζ] satisfies ∂v1/2 [ζ] = −iΛa v1/2 [ζ]. ∂z Using the superposition principle, u1/2 (z) is given by zmax u1/2 (z) = v1/2 [ζ] dζ. z
(32)
(33)
The third one-way wave equation of the multi-one-way approach is a down-going one-way wave equation with u1/2 as the source term: ∂ (34) − iΛa u1 = u1/2 . ∂z As previously, we define for all ζ in [0, zmax ] the function e1/2 [ζ] by e1/2 [ζ](x, z) = u1/2 (x, ζ)δ(z − ζ) and the function v1 by ∂ − iΛa v1 [ζ] = e1/2 [ζ]. ∂z
(35)
We assume that v1 [ζ](z) = 0 for z < ζ. Following the derivation similar as for Eq. (29), the initial boundary condition for v1 [ζ] is v1 [ζ]|z=ζ+0 = u1/2 [x, ζ]
(36)
and for z > ζ, v1 [ζ] satisfies ∂v1 [ζ] = iΛa v1 [ζ]. ∂z The superposition principle gives z u1 = u10 + v1 [ζ] dζ. 0
(37)
(38)
106
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
u10 is the solution of the down-going one-way equation with the initial boundary condition (27). It is necessary to add this term because the initial condition for u1 is not zero. In practice it is not necessary to compute the function v1/2 [ζ] (or v1 [ζ]) independently for each ζ, because only u1/2 and u1 are of interest. It is sufficient to add the initial boundary condition term (31) (or (36)) to the wave field at each step of the upward (or downward) one-way wave propagation. In such a way, the integration over ζ is performed during the wave field extrapolation procedure. Three steps of the multi-one-way modeling method are summarized as follows: 1. Solve the down-going one-way wave equation (12) with the initial boundary condition (23). As a result, the “zero-order” approximation of the wave field, u0 , is computed from the top to the bottom of the model, and the error term, Eu0 , is evaluated. 2. Solve the up-going one-way wave equation (14) with the right-hand side, Eu0 , and the initial boundary condition (24) at the bottom of the model (z = zmax ). As a result, u1/2 is computed from the bottom to the top of the model. 3. Solve the down-going one-way wave equation (15) with the right-hand side u1/2 and the initial boundary condition (27) at the top of the model. As the result, u1 is obtained.
Fig. 1. Low velocity lens model: the wave fields computed with the classic one-way scheme (a), with the multi-one-way scheme (b) and with the two-way scheme (c). The source is located at the top-left corner of the model.
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
107
The “first-order” approximation of the total field is obtained by adding u0 and u1 . It may be possible to compute higher-order approximations but in this paper we restrict the consideration only to the approximation of the “first order”. The multi-one-way modeling method is derived without any assumption on the local approximation, Λa , of the square root operator Λ. Therefore, this algorithm can be implemented with any arbitrary one-way scheme. For the numerical applications of this paper, the down-going and up-going one-way wave equations are solved with a finite-difference scheme. The square root operator is approximated by Pad´e fractions and the splitting method is used [7].
3. Numerical examples Three two-dimensional numerical examples are now described, first a small smoothed lens model to evaluate the efficiency of the new method in the case of a laterally and vertically varying velocity model and to show the importance of the two parts of the error term, secondly a small blocky wedge model to show how the method behaves with rough models, and finally the SEG/EAGE salt model to demonstrate the interest of the method on a more realistic model and to study its limitations. The results obtained with the multi-one-way approach are compared with those computed with the classic one-way method and the two-way method. The latter serves as reference. The comparisons are mainly carried out on the pressure wave field in the frequency domain. 3.1. Smoothed lens model In this first example the velocity model is the lens of 1 km of diameter centered in the middle of the computational domain with a velocity of 1.2 km/s embedded in a background velocity of 1.5 km/s. The boundary between the lens and the background has been smoothed. The pressure wave field at 15 Hz obtained with the classic one-way, the
Fig. 2. Low velocity lens model: the relative error of the amplitudes of the wave field at 5 Hz for the classic one-way result (a), for the multione-way result without the correction relative to the approximation of the square root operator (b), and for the complete multi-one-way result (c). The source is located at the middle-top of the model.
108
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
multi-one-way, and the two-way methods are displayed in Fig. 1. The shadow zone due to the low velocity lens is better modeled with the multi-one-way method than with the classic one-way method. To quantify the differences the relative errors between the two-way result and the classic one-way or the multi-one-way results are computed. The relative error of the classic one-way method reaches 40% (Fig. 2a). With the multi-one-way method the relative error is divided by a factor 5 and is smaller than 8% (Fig. 2c). If in the error term the part Λa Λa − ΛΛ is ignored, the relative error reaches 18% (Fig. 2b). This results show that the two parts of the error term are equally important when the velocity varies in lateral direction. Similar results have been obtained with other smoothed models. 3.2. Blocky wedge model The second model is a wedge model with two interfaces (Fig. 3a). The source is located at the top-left corner. The relative errors of the amplitudes of the wave field with respect to the two-way solution are displayed in Fig. 3b and c for a 5 Hz source frequency. Away from the interfaces and up to an incident angle of about 60◦ , the relative error is reduced by a factor of about 5 to 2–4% with the multi-one-way approach. Close to the interfaces and for large horizontal distances from the source (also referred as offset further in the text), the relative error is also reduced, but is still relatively large, about 40%. In fact at large offsets, not only reflected waves but also refracted wave exist. Both are modeled with the two-way scheme, but the multi-one-way model only the reflected waves. This explains why the relative errors remains large in the vicinity of the interfaces. The multi-one-way scheme still brings some improvements because the classic one-way scheme models none of these waves. To study the reflected waves, the pressure wave fields are computed at the top of the computational model for a range of frequencies and Fourier
Fig. 3. Wedge layered model (a). The relative errors of the amplitude of wave field with respect to the two-way solution of the classic one-way solution (b) and the multi-one-way solution (c).
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
109
Fig. 4. Wedge model: the time response caused by the reflection on the first interface with the multi-one-way scheme (a) and the two-way scheme (b), and the time response caused by the reflection on the second interface with the multi-one-way scheme (c) and the two-way scheme (d).
transform is applied to obtain the time series. The seismograms computed with the multi-one-way and the two-way methods are plotted in Fig. 4. For comparison, the two reflected events have been separated and the direct wave between the source and the receivers is not displayed. The traveltimes of these two events are identical with both modeling methods. To compare the dynamics, the maximum amplitudes are picked along each reflected event and displayed in Fig. 5. For both reflected events, the amplitudes obtained with the multi-one-way scheme match relatively well the ones obtained with the two-way scheme. Some discrepancies appear at large offsets. The errors on the second event are slightly larger because of a cumulative effect. In Appendix A we demonstrate that for simple layered model the multi-one-way method correctly models the reflected and transmitted waves, thanks to the three-step approach. 3.3. SEG/EAGE salt model The last example is based on the SEG/EAGE salt body example [21]. A two-dimensional slide, Fig. 6a, has been extracted from the three-dimensional model and defines the rough velocity model. A smoothed version is created
110
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
Fig. 5. Wedge model: the amplitudes of the first reflected event (a) and of the second one (b) picked on the time data generated by the multi-one-way (solid line) and the two-way (dashed line) modeling.
by applying a Gaussian filter (Fig. 6b). The multi-one-way scheme is expected to be less efficient with strong rough contrasts, such as in this case, because it uses the perturbation approach. The wave fields computed over the smoothed velocity model at 14 Hz with the classic one-way, the multi-oneway and the two-way schemes are displayed in Fig. 7. The explosive source is located at the middle of the top of the model (x = 6 km). Qualitatively, these plots show that the multi-one-way result contains more features than the one-way result and is closer to the two-way result. To quantify this remark, the transmission through the salt body is studied. The relative differences between the classic one-way or the multi-one-way and the two-way pressure fields computed at the depth z = 3.8 km are plotted in Fig. 8 for a range of frequencies. Whereas the relative errors of the classic one-way scheme range around 30%, the relative errors of the multi-one-way scheme vary around 5%. However below the edges of the salt body, around x = 4 and 8 km, where the lateral velocity variations are large, the relative error of the multi-one-way result remains large. This may be explained by the presence in the two-way results of waves traveling along the salt flanks which cannot be modeled by the one-way schemes. Using now the rough velocity model, Fig. 6a, time domain seismograms have been generated with the multi-oneway and the two-way schemes and plotted in Fig. 9a and b. Qualitatively the multi-one-way scheme models most of the reflected events existing in the two-way results. However it is noticed that some of the amplitudes are different. To quantify this remark, the maximum amplitudes of the most energetic event, corresponding to the reflection on the top of the salt, have been picked on both seismograms and compared in Fig. 9c. The large discrepancies are probably a consequence of the large variations between the sediment velocity at about 2000 m/s and the salt velocity at 4500 m/s. This is a limitation of the proposed approach because of the use of the perturbation theory.
Fig. 6. SEG/EAGE salt models.
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
111
Fig. 7. Smoothed SEG/EAGE salt model: the wave fields computed with the classic one-way scheme (a), with the multi-one-way scheme (b) and with the two-way scheme (c).
Fig. 8. Smoothed SEG/EAGE salt model: the relative error vs. frequency and lateral position of the amplitude spectrum recorded at the depth z = 3.8 km with respect to the two-way solution for the classic one-way (a) and the multi-one-way (b) modeling.
112
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
Fig. 9. Seismograms obtained on the rough SEG/EAGE salt model with the two-way (a) and multi-one-way (b) modeling. The amplitudes of the most energetic reflected events obtained with the two-way (dashed) scheme and the multi-one-way (solid) scheme (c).
4. Conclusions and discussion The multi-one-way method has been developed to improve the estimation of the amplitudes of the wave fields based on one-way equations. It is an iterative method built using the perturbation approach. It consists of the set of one-way equations with right-hand sides. The right-hand sides are constructed from the error operator, the difference between the two-way wave operator and the product of the up-going and down-going one-way wave operators. The proposed approach takes into account vertical and horizontal velocity variations of the medium. In practice the new scheme is five times more expensive than the classic one-way scheme, because three one-way equations are solved and the error term is computed. The errors in the amplitudes of the wave field are almost one order of magnitude smaller with the new approach than with the classic one-way one. This is shown for smooth and rough velocity models. With the multi-one-way approach, reflection seismic data can be generated. As long as the velocity contrasts around the interfaces are small, the reflected and the transmitted amplitudes are well estimated even if some cumulative errors may appear.
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
113
However, because of the use of the perturbation method, large velocity contrasts still induce significant errors in the estimations of the reflected amplitudes in the seismograms. The multi-one-way modeling method, detailed in this paper, can be applied to inverse problems such as seismic migration in the same way as classic one-way methods.
Acknowledgements The authors are thankful to reviewer of this paper, Dr. Norman Bleistein, for his comments which allowed us to improve the manuscript quality. This research was supported by the grants CRDF RG0-1318(3), RFBR 02-05-65081 and A03-2.13-366.
Appendix A The transmission and reflection coefficients modeled by the multi-one-way scheme are analyzed in this appendix. The model consists of two homogeneous half spaces separated by a flat interface at depth z = z0 . The source is located at depth z = 0. The wavenumber for the upper layer (z < z0 ) is denoted by k1 and for the lower layer (z > z0 ) by k2 . With the Fourier transform U(kx , z) =
∞
−∞
u(x, z) exp(−ikx x) dx,
the multi-one-way equations become in the horizontal wavenumber domain: • Step 1: 2iΛU0 |z=0 =
1 , 2
∂U0 (z) = iΛU0 (z). ∂z
(39)
∂U1/2 (z) = −iΛU1/2 (z) − iU0 (z0 )(Λ2 − Λ1 )δ(z − z0 ) ∂z
(40)
• Step 2: U1/2 (zmax ) = 0,
(the second term of the right-hand side corresponds to the error term, EU0 ). • Step 3: − 2iΛU1 (z)|z=0 = U1/2 (0),
∂U1 (z) = iΛU1 (z) + U1/2 (z), ∂z
(41)
with Λ = k(z)2 − kx2 . Λ1 = k12 − kx2 is the restriction of Λ in the upper half space and Λ2 = k22 − kx2 is the restriction of Λ in the lower half space. The dependence on the horizontal wavenumber, kx , and the angular frequency, ω, have not been written to simplify the notation.
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
114
The analytical solutions are z exp (i 0 Λ(z ) dz ) U0 (z) = , U1/2 (z) = iU0 (z0 )(Λ2 − Λ1 )θ(z0 − z) exp (iΛ1 (z0 − z)), 2iΛ1 z min(z,z0 ) z U1/2 (0) exp i Λ(z ) dz + dz U1/2 (z ) exp i U1 (z) = − Λ(z ) dz . 2iΛ 0 0 z
(42)
In the upper half space, z < z0 , this gives U0 (z) =
exp (iΛ1 z) , 2iΛ1
U1 (z) =
Λ1 − Λ2 exp (iΛ1 (2z0 − z)) . 2Λ1 2iΛ1
U1/2 (z) =
Λ2 − Λ1 exp (iΛ1 (2z0 − z)), 2Λ1 (43)
U1 is the scattered (reflected) wave field. If we define Ra =
Λ1 −Λ2 2Λ1 ,
we have in the upper layer
U1 (z) = Ra U0 (z0 ) exp (iΛ1 (z0 − z)).
(44)
Ra is the reflection coefficient modeled by the multi-one-way scheme and is an approximation of the plane wave 1 −Λ2 reflection coefficient, R = Λ Λ1 +Λ2 , for small velocity contrasts and small horizontal wavenumbers (i.e. small inci-
1 dence angles). It is interesting to notice that the difference between U1/2 and U1 is just the factor 2iΛ . However the phases of U1/2 and U1 are identical. The multi-one-way scheme, as described in this paper, mainly corrects the amplitudes. In the lower half space, z > z0 , the analytical solutions of Eqs. (39)–(41) are
U0 (z) =
exp (iΛ1 z0 + iΛ2 (z − z0 )) , 2iΛ1
U1/2 (z) = 0,
U1 = Ra
exp (iΛ1 z0 + iΛ2 (z − z0 )) . 2iΛ1
(45)
The total field, namely the transmitted field, after the interface at z = z0 is equal to U0 + U1 and we have U0 + U1 = (1 + Ra )U0 .
(46)
1 + Ra is an approximation of the transmission coefficient for small incident angles and small contrasts, since the plane wave transmission coefficient is 1 + R [22]. In this study the square root operator, Λ, and not its approximation has been used, which restricts the validity of this demonstration to small incidence angles. However this appendix demonstrates that the multi-one-way scheme allows estimating of the amplitudes of reflected and transmitted waves in a simple layered earth for small incidence angles and small velocity contrasts and explains why reflection data can be generated with this scheme.
References [1] [2] [3] [4]
W.A. Mulder, R.-E. Plessix, A comparison between one-way and two-way wave-equation migration, Geophysics 69 (6) (2004) 1491–1504. D. Lee, A.D. Pierce, Parabolic equation development in recent decade, J. Comput. Acoust. 3 (2) (1995) 95–173. J.F. Claerbout, Imaging the Earth’s Interior, Blackwell Scientific Publications, 1983. M.A. Leontovich, V.A. Fock, Solution of the problem of propagation of electromagnetic waves along the earth’s surface by the method of parabolic equation, J. Exp. Theoret. Phys. 16 (1946) 557–573. [5] F. Tappert, The parabolic approximation method, in: J.B. Keller, J.S. Papadakis (Eds.), Wave Propagation and Underwater Acoustics, American Institute of Physics, vol. 70, Springer, New York, 1977, pp. 224–280. [6] J.F. Claerbout, Coarse grid calculation of waves in inhomogeneous media with application to delineation of complicated seismic structure, Geophysics 35 (1970) 407–418.
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
115
[7] F. Collino, P. Joly, Splitting of operators, alternate directions, and paraxial approximations for three-dimensional wave equation, SIAM J. Sci. Comput. 16 (5) (1999) 1019–1048. [8] M. Kern, A non-split 3D migration algorithm, in: S. Hassanzadeh (Ed.), Mathematical Methods in Geophysical Imaging, SPIE, 1995. [9] D. Ristow, T. Ruhl, 3-D implicit finite-difference migration by multi-way splitting, Geophysics 62 (2) (1997) 554–567. [10] J. Gazdag, Wave equation migration with the phase shift method, Geophysics 43 (1978) 1342–1351. [11] J. Gazdag, P. Squazerro, Migration of seismic data by phase shift plus interpolation, Geophysics 48 (1984) 124–131. [12] R.-S. Wu, Wide-angle elastic wave one-way propagation in heterogeneous media and an elastic wave complex-screen method, J. Geophys. Res. 99 (B1) (1994) 751–766. [13] R.-S. Wu, L.-J. Huang, Reflected wave modeling in heterogeneous acoustic media using De Wolf approximation, in: S. Hassenzadeh (Eds.), Mathematical Methods in Geophysical Imaging, SPIE, 1995, pp. 176–186. [14] X.-B. Xie, R.-S. Wu, 3D elastic wave modeling using the complex screen method, in: Proceedings of the SEG 66th Annual Meeting, 1996. [15] S. Jin, R.-S. Wu, Prestack depth migration using a hybrid pseudo-screen propagator, in: Proceedings of the SEG 69th Annual Meeting, 1999. [16] S. Jin, Experimenting with the hybrid pseudo-screen migration, in: Proceedings of the SEG 70th Annual Meeting, 2000. [17] Yu. Zhang, G. Zhang, N. Bleistein, Wave equation migration arising from true amplitude one-way equations, Prob. Inverses 19 (2003) 1113–1138. [18] A. Bamberger, B. Engquist, L. Halpern, P. Joly, Parabolic wave equation approximations in heterogeneous media, SIAM J. Appl. Math. 48 (1) (1991) 99–128. [19] C. Wapenaar, J. Grimbergen, Reciprocity theorems for one-way wavefields, Geophys. J. Int. 127 (1996) 169–177. [20] M. De Hoop, A. Gautesen, Uniform asymptotic expansion of the square root Helmholtz operator and the one-way wave propagator, SIAM J. Appl. Math. 63 (3) (2003) 77–800. [21] F. Aminzadeh, J. Brac, T. Kunz, 3-D Salt and Overthrust Model, S.E.G., Tulsa, Oklahoma, 1997. [22] K. Aki, P.G. Richard, Quantitative Seismology, Theory and Methods, Freeman abd Company, San Francisco, 1980.
Improved amplitude multi-one-way modeling method D. Kiyashchenkoa,∗ , R.-E. Plessixb , B. Kashtana , V. Troyana a
Saint-Petersburg State University, Laboratory of Elastic Media Dynamics, Ulianovskaya Street 1, 198504 St. Petersburg, Russian Federation b Shell International E&P, P.O. Box 60, 2280 AB Rijswijk, The Netherlands Received 29 September 2004; received in revised form 3 May 2005; accepted 9 June 2005 Available online 2 September 2005
Abstract Classical paraxial wave equations, also called one-way wave equations, lead to efficient numerical schemes in threedimensional spaces for computing the wave field. The kinematics is well approximated by these algorithms, however the dynamics is generally not correctly modeled when the medium is inhomogeneous. This prevents from using one-way wave equation scheme for solution of inverse problems. In order to improve estimation of the amplitudes, a new multi-one-way scheme is introduced. This approach is based on an iterative solution of the factorized two-way wave equation with a right-hand side incorporating the information about medium heterogeneities. The numerical scheme is a succession of several classic oneway schemes. The proposed approach takes into account vertical and horizontal velocity variations of the medium and allows modeling reflected waves. Numerical results for several two-dimensional examples with smooth and rough velocity contrasts are compared with finite-difference solutions of the two-way wave equation. It is shown that the multi-one-way scheme significantly reduces the error in the amplitude estimates of the wave field. © 2005 Elsevier B.V. All rights reserved. PACS: 43.20.Bi Keywords: One-way equation; True amplitude
1. Introduction Wide range of inverse problems require accurate solutions of the wave equation, for example in geophysics or underwater acoustics. The finite-difference solution of the full wave equation still remains expensive in threedimensional spaces and is not yet used for practical applications [1]. The high-frequency (ray-based) solutions ∗
Corresponding author. Tel.: +7 8124284317; fax: +7 8124284317. E-mail addresses: [email protected] (D. Kiyashchenko), [email protected] (R.-E. Plessix).
0165-2125/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.wavemoti.2005.06.003
100
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
suffer from some limitations over complex media. The paraxial approximation of the wave equation then gives an alternative approach [2]. This approach consists of particularizing one direction and extrapolating the wave field step by step along this direction. This approach is also called classic one-way because the wave field contains only the up-going (or down-going) waves, the reflected waves propagating in the direction opposite to preferred one are not modeled. In contrast, the full wave equation is called two-way wave equation in this paper. One of the drawbacks of the classic one-way wave equation approximation is the large errors in the amplitude estimates of the wave field. For example, in exploration geophysics these errors generally prevent us from using the classic one-way scheme in preserved amplitude migration algorithms. Migration is the technique which is commonly used in seismic exploration for mapping the time-domain data (seismograms) to the depth domain. This gives a structural image of the subsurface. The preserved amplitude migration allows to relate the amplitudes of the imaged reflectors with media parameters. Generally, a migration algorithm requires the computation of the incident down-going wave field propagating from the source point and the back-propagated down-going wave field of the seismic data [3]. A preserved amplitude migration requires a good amplitude estimates of those wave fields. The idea of the one-way approach can be found in pioneer works of Leontovich and Fock [4] and Tappert [5], who considered the solution of the wave equation in a preferred propagation direction. Claerbout [6,3] introduced this approach in exploration geophysics. One-way wave equations are based on an approximation of the two-way dispersion curve. This dispersion curve is obtained from the time–space domain wave equation by temporal and spatial Fourier transforms and represents the wave equation in the frequency–wavenumber domain. This derivation assumes a constant velocity field. The wavenumber in the preferred direction is expressed in terms of wavenumbers in the other directions. The one-way wave equation in the spatial coordinates is then formally derived by replacing the wavenumbers by the spatial derivatives in the dispersion curve. This one-way equation involves the non-local square root operator. In order to obtain an efficient numerical scheme for spatially varying background velocity and to compute the wave field step by step along the preferred direction with a marching approach, a first possibility is to approximate the square root operator with local operators, for instance with Pad´e fractions leading to the classic 15◦ , 45◦ or 60◦ approximation scheme (see [3,7–9]). A second possibility is to first extrapolate the wave field in the frequency–wavenumber domain assuming a constant velocity model and then to apply a correction using Born approximation. This corresponds to the phase-shift methods [10,11] and the screen propagator methods [12–16]. Up to a given angle from the preferred direction, the one-way wave equation methods correctly model the kinematics of waves. Unfortunately, the dynamic properties (amplitudes) are not accurately estimated with most of the one-way schemes. Zhang et al. [17] have shown that the classic one-way eikonal equation give the same traveltimes as the two-way eikonal equation, but it is not the case for the transport equations, which govern the amplitudes in high-frequency approximation. In the literature, few papers have addressed the problem of the amplitude behavior of the one-way scheme with spatially varying velocity. Bamberger et al. [18] suggested a method to take into account the transmission losses of the normal incident waves through the artificial wavefield discontinuities introduced during the marching approach. Zhang et al. [17] proposed to add a term to the one-way wave equation to make the eikonal and transport equations of the latter matching those obtained from the two-way wave equations. Interesting formulations of one-way and two-way wave propagation problems have also been suggested by Wapenaar and Grimbergen [19] and De Hoop and Gautesen [20]. The approach is based on a decomposition of the acoustic wave equation into two coupled one-way equations for upward and downward propagating waves. When the medium is homogeneous, these equations are decoupled. In the case of not too large velocity gradients, the usual one-way approximation is used. This approach allows an iterative solution to take into account wave transmission and reflection phenomena. However, for three-dimensional applications the expressions for the square root operator proposed in these papers may lead to very expensive computational procedure. In this paper, a different approach is derived to correct for the amplitude errors of the classic one-way wave equation due to velocity variations. The two-way wave equation is rewritten as the product of two one-way wave equations plus an error (correction) term. An iterative solution is then conducted using the perturbation theory. The scheme consists of a number of classic one-way wave equations with right-hand sides. The right-hand sides take into account the error term. This approach corrects for the errors due to the velocity variations and models the primary
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
101
reflections with rough velocity models. This method is called multi-one-way method. The use of the perturbation theory assumes small velocity contrasts. In practical applications, three one-way equations are solved and the error term is computed. This means that this method is about five times more expensive than the classic one-way method. The outline of the paper is following. First the derivation of the classic one-way wave equation is recalled and the multi-one-way method is explained. Then three two-dimensional examples are described to demonstrate the relevance of the approach and its limitations with smooth and rough velocity models. In Appendix A, the modeling of the transmission and reflection coefficients with the multi-one-way method is analyzed for the case of a plane interface.
2. Two-way, one-way and multi-one-way modeling methods 2.1. The one-way approach One-way wave equation is a paraxial approximation of the full two-way wave equation after the choice of a preferred direction of the wave propagation. In this work, the preferred direction is the vertical axis (z-direction). The constant density two-way acoustic wave equation in the space-frequency domain reads 2 ∂ ω2 (1) + 2 + x u = f (x, z, ω), ∂z2 c ∂ ∂ ∂ with x = ∂x 2 in two dimensions and x = ∂x2 + ∂y2 in three dimensions. c(x, z) is the velocity, u(x, z) the pressure field and f (x, z, ω) the source term. The vector x has two components, x ≡ (x, y), in three dimensions and one component, x ≡ x, in two dimensions. The derivation of the one-way wave equation is based on the factorization of the two-way operator 2 ∂ 2 (2) + Λ A= ∂z2 2
2
2
into the product of up-going and down-going wave one-way operators: ∂ ∂ L= + iΛ − iΛ , (3) ∂z ∂z ω . where Λ(x, z) = k2 + x with k(x, z) = c(x,z) L equals A only when the velocity does not depend on the depth coordinate z. Otherwise the square root operator, ∂ Λ, depends on z, and does not commute with the operator ∂z . In fact L satisfies ∂ ∂Λ ∂ + iΛ − iΛ = A − i . (4) L= ∂z ∂z ∂z In practice, a local approximation, Λa , of the square root operator, Λ, is used. The product of the up-going and down-going one-way wave operators, with the local approximation of the square root, gives ∂ ∂ ∂Λa La = + iΛa − iΛa = A + (Λ2a − Λ2 ) − i . (5) ∂z ∂z ∂z With the error operator, E, defined by E = −i
∂Λa + (Λ2a − Λ2 ), ∂z
(6)
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
102
this gives La − E = A.
(7)
The classic one-way approach neglects the error term, E. Assuming no source term for the time being (f (x, z, ω) = 0), this leads to ∂ ∂ (8) + iΛa − iΛa = 0 ∂z ∂z or if the derivation is repeated by interchanging the up- and down-going one-way wave operators, then ∂ ∂ − iΛa + iΛa = 0. ∂z ∂z
(9)
From Eqs. (8) and (9) the classic one-way equations are obtained: ∂u = iΛa u, ∂z
(10)
with = 1 for the classic down-going one-way wave equation and = −1 for the up-going one. These one-way approximations do not correctly estimate the amplitudes as shown by Zhang et al. [17]. In their study, they demonstrate that the travel times obtained from the one-way wave equation satisfy the two-way eikonal equation, but the amplitudes of the one-way solution do not satisfy the two-way transport equation. This proves that even in the high-frequency sense, the amplitudes are not correctly modeled with the classic one-way method in inhomogeneous media. The error term, E, explains the differences between the one-way and the two-way solutions. The expression a (6) has two parts. The first one, −i ∂Λ ∂z , is non-zero if vertical velocity variations are present. The second one, (Λ2a − Λ2 ), is non-zero for large incidence angles, i.e. large angles between the vertical (the z-direction) and the wave propagation directions, due to the square root approximation, and in the presence of lateral velocity variations. 2.2. Multi-one-way modeling method In order to obtain a better amplitude estimation of the wave field with an approach based on the one-way wave equation, the (two-way) solution, u, of La u = Eu
(11)
is approximated using the perturbation method. The solution u0 of the down-going one-way wave equation ∂u0 = iΛa u0 ∂z can be considered as a “zero-order” approximation of the solution of Eq. (11). A “first-order” approximation, u0 + u1 , is obtained when the correction term u1 satisfies ∂ ∂ + iΛa − iΛa u1 = Eu0 . ∂z ∂z
(12)
(13)
Eq. (13) can be solved in two steps. First the solution, u1/2 , of the up-going one-way wave equation with the right-hand side (source term) equal to the error term evaluated at u0 is computed: ∂ (14) + iΛa u1/2 = Eu0 . ∂z
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
103
Secondly the solution u1 of the down-going one-way wave equation with u1/2 as the right-hand side is computed: ∂ − iΛa u1 = u1/2 . (15) ∂z The set of the one-way wave equations, Eqs. (12), (14) and (15), should be solved sequentially in order to obtain a “first-order” approximation of the solution of the two-way wave equation. Therefore, this approach is called multione-way modeling method. Appropriate initial boundary conditions for the above mentioned equations should be defined. 2.3. Boundary conditions for the one-way equations in multi-one-way scheme In the previous subsections the study of the one-way wave equations was performed in the spatial domain with f (x, z, ω) = 0 in Eq. (1). To take into account the source term it is sufficient to consider the case f (x, z, ω) = δ(z)fs (x), with δ(z) a dirac delta-function in z-coordinate. The general solution can be obtained with the superposition principle. With the source term, Eq. (11) becomes La u = Eu + δ(z)fs (x). Neglecting the error term leads to the equation for the “zero-order” approximation, u0 : ∂ ∂ + iΛa − iΛa u0 = δ(z)fs (x). ∂z ∂z
(16)
(17)
Theoretically Eq. (17) is solved in two steps. The first step consists of solving the up-going one-way wave equation: ∂ (18) + iΛa g = δ(z)fs (x). ∂z The second step consists of solving the down-going one-way wave equation: ∂ − iΛa u0 = g. ∂z
(19)
The radiation free condition at z = zmax , the bottom of the computational domain, gives g(zmax ) = 0, which is the initial condition of the up-going Eq. (18). Therefore g(z) = 0 for all z > 0 for a source term f which is non-zero only at z = 0, the top of computational domain. In this way only Eq. (19) needs to be solved with the appropriate initial boundary condition. So, the first step of the multi-one-way scheme consists only in solving a down-going one-way wave equation. To derive the initial condition of Eq. (12) (or Eq. (19)), the first equation of the multi-one-way scheme, we follow [7]. As the right-hand side of Eq. (17) is proportional to δ(z)fs (x) and the left-hand side contains a second derivative in z, the first order derivative of u0 at z = 0 is discontinuous and satisfies: 1 ∂u0 = ± fs (x). (20) ∂z z=±0 2 With Eq. (12) this gives iΛa u0 |z=+0 = 21 fs (x),
(21)
or u0 |z=+0 =
1 Λa Λ−2 a fs (x). 2i
(22)
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
104
−2 Λ−2 a fs (x) can be approximated by Λ fs (x) which is the solution of a two-dimensional two-way wave equation in three dimensions (one-dimensional two-way wave equation in two dimensions). The initial boundary condition for Eq. (12) is then
u0 |z=+0 =
1 Λa h(x), 2i
(23)
2 where h(x) is the solution of x + c2 ω(x,0) h(x) = fs (x) with absorbing boundary conditions. The second equation of the multi-one-way approach, Eq. (14), is an up-going one-way wave equation. The initial boundary condition is deduced from the radiation condition at the bottom of computational domain (z = zmax ). This condition means that the wave field u1 should not contain up-going waves at the bottom of computational domain. Then u1/2 |z=zmax satisfies ∂ u1/2 |z=zmax = = 0. (24) − iΛa u1 ∂z z=zmax The third equation of the multi-one-way approach, Eq. (15), is a down-going one-way wave equation. The initial boundary condition is obtained assuming that the wave field u1 does not contain reflections from the top of the computational domain at z = 0. This means that all plane waves contained in the wave field u1 should propagate upward at z = 0. This gives ∂ =0 (25) + iΛa u1 ∂z z=+0 and with Eq. (15) − 2iΛa u1 |z=+0 = u1/2 |z=+0 .
(26)
The initial boundary condition u1 |z=0 is then computed as follows: 1 (27) Λ a h1 , 2i 2 where h1 is the solution of x + c2 ω(x,0) h1 = u1/2 |z=+0 . The initial boundary conditions for the multi-one-way set of Eqs. (12), (14), and (15) are then given by Eqs. (23), (24), and (27). u1 |z=+0 = −
2.4. Solution of one-way equations with non-zero source term The multi-one-way method consists of successive solutions of three one-way wave equations. The last two are inhomogeneous, as their right-hand side (source term) is non-zero. It is interesting to discuss in more details the solution of these equations. The second one-way wave equation of the multi-one-way approach is an up-going one-way wave equation with Eu0 as source term: ∂ (28) + iΛa u1/2 = Eu0 , ∂z with the initial boundary condition u1/2 (zmax ) = 0. We define for all ζ in [0, zmax ] the function e0 by e0 [ζ](x, z) = (Eu0 )(x, ζ)δ(z − ζ)
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
and v1/2 [ζ] by ∂ + iΛa v1/2 [ζ] = e0 [ζ], ∂z
105
(29)
with v1/2 [ζ](z) = 0 for z > ζ. At z = ζ an initial boundary condition is needed. v1/2 [ζ](z) at z = ζ is discontinuous due to the presence of the delta-function in the right-hand side of Eq. (29). The integration of Eq. (29) over a small interval (ζ − ε, ζ + ε) gives z=ζ+ε dz iΛa v1/2 [ζ](z) = Eu0 (x, ζ). (30) v1/2 [ζ](ζ + ε) − v1/2 [ζ](ζ − ε) + z=ζ−ε
The initial condition is obtained when ε goes to 0: v1/2 [ζ]|z=ζ−0 = −(Eu0 )(x, ζ)
(31)
and for z < ζ the function v1/2 [ζ] satisfies ∂v1/2 [ζ] = −iΛa v1/2 [ζ]. ∂z Using the superposition principle, u1/2 (z) is given by zmax u1/2 (z) = v1/2 [ζ] dζ. z
(32)
(33)
The third one-way wave equation of the multi-one-way approach is a down-going one-way wave equation with u1/2 as the source term: ∂ (34) − iΛa u1 = u1/2 . ∂z As previously, we define for all ζ in [0, zmax ] the function e1/2 [ζ] by e1/2 [ζ](x, z) = u1/2 (x, ζ)δ(z − ζ) and the function v1 by ∂ − iΛa v1 [ζ] = e1/2 [ζ]. ∂z
(35)
We assume that v1 [ζ](z) = 0 for z < ζ. Following the derivation similar as for Eq. (29), the initial boundary condition for v1 [ζ] is v1 [ζ]|z=ζ+0 = u1/2 [x, ζ]
(36)
and for z > ζ, v1 [ζ] satisfies ∂v1 [ζ] = iΛa v1 [ζ]. ∂z The superposition principle gives z u1 = u10 + v1 [ζ] dζ. 0
(37)
(38)
106
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
u10 is the solution of the down-going one-way equation with the initial boundary condition (27). It is necessary to add this term because the initial condition for u1 is not zero. In practice it is not necessary to compute the function v1/2 [ζ] (or v1 [ζ]) independently for each ζ, because only u1/2 and u1 are of interest. It is sufficient to add the initial boundary condition term (31) (or (36)) to the wave field at each step of the upward (or downward) one-way wave propagation. In such a way, the integration over ζ is performed during the wave field extrapolation procedure. Three steps of the multi-one-way modeling method are summarized as follows: 1. Solve the down-going one-way wave equation (12) with the initial boundary condition (23). As a result, the “zero-order” approximation of the wave field, u0 , is computed from the top to the bottom of the model, and the error term, Eu0 , is evaluated. 2. Solve the up-going one-way wave equation (14) with the right-hand side, Eu0 , and the initial boundary condition (24) at the bottom of the model (z = zmax ). As a result, u1/2 is computed from the bottom to the top of the model. 3. Solve the down-going one-way wave equation (15) with the right-hand side u1/2 and the initial boundary condition (27) at the top of the model. As the result, u1 is obtained.
Fig. 1. Low velocity lens model: the wave fields computed with the classic one-way scheme (a), with the multi-one-way scheme (b) and with the two-way scheme (c). The source is located at the top-left corner of the model.
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
107
The “first-order” approximation of the total field is obtained by adding u0 and u1 . It may be possible to compute higher-order approximations but in this paper we restrict the consideration only to the approximation of the “first order”. The multi-one-way modeling method is derived without any assumption on the local approximation, Λa , of the square root operator Λ. Therefore, this algorithm can be implemented with any arbitrary one-way scheme. For the numerical applications of this paper, the down-going and up-going one-way wave equations are solved with a finite-difference scheme. The square root operator is approximated by Pad´e fractions and the splitting method is used [7].
3. Numerical examples Three two-dimensional numerical examples are now described, first a small smoothed lens model to evaluate the efficiency of the new method in the case of a laterally and vertically varying velocity model and to show the importance of the two parts of the error term, secondly a small blocky wedge model to show how the method behaves with rough models, and finally the SEG/EAGE salt model to demonstrate the interest of the method on a more realistic model and to study its limitations. The results obtained with the multi-one-way approach are compared with those computed with the classic one-way method and the two-way method. The latter serves as reference. The comparisons are mainly carried out on the pressure wave field in the frequency domain. 3.1. Smoothed lens model In this first example the velocity model is the lens of 1 km of diameter centered in the middle of the computational domain with a velocity of 1.2 km/s embedded in a background velocity of 1.5 km/s. The boundary between the lens and the background has been smoothed. The pressure wave field at 15 Hz obtained with the classic one-way, the
Fig. 2. Low velocity lens model: the relative error of the amplitudes of the wave field at 5 Hz for the classic one-way result (a), for the multione-way result without the correction relative to the approximation of the square root operator (b), and for the complete multi-one-way result (c). The source is located at the middle-top of the model.
108
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
multi-one-way, and the two-way methods are displayed in Fig. 1. The shadow zone due to the low velocity lens is better modeled with the multi-one-way method than with the classic one-way method. To quantify the differences the relative errors between the two-way result and the classic one-way or the multi-one-way results are computed. The relative error of the classic one-way method reaches 40% (Fig. 2a). With the multi-one-way method the relative error is divided by a factor 5 and is smaller than 8% (Fig. 2c). If in the error term the part Λa Λa − ΛΛ is ignored, the relative error reaches 18% (Fig. 2b). This results show that the two parts of the error term are equally important when the velocity varies in lateral direction. Similar results have been obtained with other smoothed models. 3.2. Blocky wedge model The second model is a wedge model with two interfaces (Fig. 3a). The source is located at the top-left corner. The relative errors of the amplitudes of the wave field with respect to the two-way solution are displayed in Fig. 3b and c for a 5 Hz source frequency. Away from the interfaces and up to an incident angle of about 60◦ , the relative error is reduced by a factor of about 5 to 2–4% with the multi-one-way approach. Close to the interfaces and for large horizontal distances from the source (also referred as offset further in the text), the relative error is also reduced, but is still relatively large, about 40%. In fact at large offsets, not only reflected waves but also refracted wave exist. Both are modeled with the two-way scheme, but the multi-one-way model only the reflected waves. This explains why the relative errors remains large in the vicinity of the interfaces. The multi-one-way scheme still brings some improvements because the classic one-way scheme models none of these waves. To study the reflected waves, the pressure wave fields are computed at the top of the computational model for a range of frequencies and Fourier
Fig. 3. Wedge layered model (a). The relative errors of the amplitude of wave field with respect to the two-way solution of the classic one-way solution (b) and the multi-one-way solution (c).
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
109
Fig. 4. Wedge model: the time response caused by the reflection on the first interface with the multi-one-way scheme (a) and the two-way scheme (b), and the time response caused by the reflection on the second interface with the multi-one-way scheme (c) and the two-way scheme (d).
transform is applied to obtain the time series. The seismograms computed with the multi-one-way and the two-way methods are plotted in Fig. 4. For comparison, the two reflected events have been separated and the direct wave between the source and the receivers is not displayed. The traveltimes of these two events are identical with both modeling methods. To compare the dynamics, the maximum amplitudes are picked along each reflected event and displayed in Fig. 5. For both reflected events, the amplitudes obtained with the multi-one-way scheme match relatively well the ones obtained with the two-way scheme. Some discrepancies appear at large offsets. The errors on the second event are slightly larger because of a cumulative effect. In Appendix A we demonstrate that for simple layered model the multi-one-way method correctly models the reflected and transmitted waves, thanks to the three-step approach. 3.3. SEG/EAGE salt model The last example is based on the SEG/EAGE salt body example [21]. A two-dimensional slide, Fig. 6a, has been extracted from the three-dimensional model and defines the rough velocity model. A smoothed version is created
110
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
Fig. 5. Wedge model: the amplitudes of the first reflected event (a) and of the second one (b) picked on the time data generated by the multi-one-way (solid line) and the two-way (dashed line) modeling.
by applying a Gaussian filter (Fig. 6b). The multi-one-way scheme is expected to be less efficient with strong rough contrasts, such as in this case, because it uses the perturbation approach. The wave fields computed over the smoothed velocity model at 14 Hz with the classic one-way, the multi-oneway and the two-way schemes are displayed in Fig. 7. The explosive source is located at the middle of the top of the model (x = 6 km). Qualitatively, these plots show that the multi-one-way result contains more features than the one-way result and is closer to the two-way result. To quantify this remark, the transmission through the salt body is studied. The relative differences between the classic one-way or the multi-one-way and the two-way pressure fields computed at the depth z = 3.8 km are plotted in Fig. 8 for a range of frequencies. Whereas the relative errors of the classic one-way scheme range around 30%, the relative errors of the multi-one-way scheme vary around 5%. However below the edges of the salt body, around x = 4 and 8 km, where the lateral velocity variations are large, the relative error of the multi-one-way result remains large. This may be explained by the presence in the two-way results of waves traveling along the salt flanks which cannot be modeled by the one-way schemes. Using now the rough velocity model, Fig. 6a, time domain seismograms have been generated with the multi-oneway and the two-way schemes and plotted in Fig. 9a and b. Qualitatively the multi-one-way scheme models most of the reflected events existing in the two-way results. However it is noticed that some of the amplitudes are different. To quantify this remark, the maximum amplitudes of the most energetic event, corresponding to the reflection on the top of the salt, have been picked on both seismograms and compared in Fig. 9c. The large discrepancies are probably a consequence of the large variations between the sediment velocity at about 2000 m/s and the salt velocity at 4500 m/s. This is a limitation of the proposed approach because of the use of the perturbation theory.
Fig. 6. SEG/EAGE salt models.
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
111
Fig. 7. Smoothed SEG/EAGE salt model: the wave fields computed with the classic one-way scheme (a), with the multi-one-way scheme (b) and with the two-way scheme (c).
Fig. 8. Smoothed SEG/EAGE salt model: the relative error vs. frequency and lateral position of the amplitude spectrum recorded at the depth z = 3.8 km with respect to the two-way solution for the classic one-way (a) and the multi-one-way (b) modeling.
112
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
Fig. 9. Seismograms obtained on the rough SEG/EAGE salt model with the two-way (a) and multi-one-way (b) modeling. The amplitudes of the most energetic reflected events obtained with the two-way (dashed) scheme and the multi-one-way (solid) scheme (c).
4. Conclusions and discussion The multi-one-way method has been developed to improve the estimation of the amplitudes of the wave fields based on one-way equations. It is an iterative method built using the perturbation approach. It consists of the set of one-way equations with right-hand sides. The right-hand sides are constructed from the error operator, the difference between the two-way wave operator and the product of the up-going and down-going one-way wave operators. The proposed approach takes into account vertical and horizontal velocity variations of the medium. In practice the new scheme is five times more expensive than the classic one-way scheme, because three one-way equations are solved and the error term is computed. The errors in the amplitudes of the wave field are almost one order of magnitude smaller with the new approach than with the classic one-way one. This is shown for smooth and rough velocity models. With the multi-one-way approach, reflection seismic data can be generated. As long as the velocity contrasts around the interfaces are small, the reflected and the transmitted amplitudes are well estimated even if some cumulative errors may appear.
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
113
However, because of the use of the perturbation method, large velocity contrasts still induce significant errors in the estimations of the reflected amplitudes in the seismograms. The multi-one-way modeling method, detailed in this paper, can be applied to inverse problems such as seismic migration in the same way as classic one-way methods.
Acknowledgements The authors are thankful to reviewer of this paper, Dr. Norman Bleistein, for his comments which allowed us to improve the manuscript quality. This research was supported by the grants CRDF RG0-1318(3), RFBR 02-05-65081 and A03-2.13-366.
Appendix A The transmission and reflection coefficients modeled by the multi-one-way scheme are analyzed in this appendix. The model consists of two homogeneous half spaces separated by a flat interface at depth z = z0 . The source is located at depth z = 0. The wavenumber for the upper layer (z < z0 ) is denoted by k1 and for the lower layer (z > z0 ) by k2 . With the Fourier transform U(kx , z) =
∞
−∞
u(x, z) exp(−ikx x) dx,
the multi-one-way equations become in the horizontal wavenumber domain: • Step 1: 2iΛU0 |z=0 =
1 , 2
∂U0 (z) = iΛU0 (z). ∂z
(39)
∂U1/2 (z) = −iΛU1/2 (z) − iU0 (z0 )(Λ2 − Λ1 )δ(z − z0 ) ∂z
(40)
• Step 2: U1/2 (zmax ) = 0,
(the second term of the right-hand side corresponds to the error term, EU0 ). • Step 3: − 2iΛU1 (z)|z=0 = U1/2 (0),
∂U1 (z) = iΛU1 (z) + U1/2 (z), ∂z
(41)
with Λ = k(z)2 − kx2 . Λ1 = k12 − kx2 is the restriction of Λ in the upper half space and Λ2 = k22 − kx2 is the restriction of Λ in the lower half space. The dependence on the horizontal wavenumber, kx , and the angular frequency, ω, have not been written to simplify the notation.
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
114
The analytical solutions are z exp (i 0 Λ(z ) dz ) U0 (z) = , U1/2 (z) = iU0 (z0 )(Λ2 − Λ1 )θ(z0 − z) exp (iΛ1 (z0 − z)), 2iΛ1 z min(z,z0 ) z U1/2 (0) exp i Λ(z ) dz + dz U1/2 (z ) exp i U1 (z) = − Λ(z ) dz . 2iΛ 0 0 z
(42)
In the upper half space, z < z0 , this gives U0 (z) =
exp (iΛ1 z) , 2iΛ1
U1 (z) =
Λ1 − Λ2 exp (iΛ1 (2z0 − z)) . 2Λ1 2iΛ1
U1/2 (z) =
Λ2 − Λ1 exp (iΛ1 (2z0 − z)), 2Λ1 (43)
U1 is the scattered (reflected) wave field. If we define Ra =
Λ1 −Λ2 2Λ1 ,
we have in the upper layer
U1 (z) = Ra U0 (z0 ) exp (iΛ1 (z0 − z)).
(44)
Ra is the reflection coefficient modeled by the multi-one-way scheme and is an approximation of the plane wave 1 −Λ2 reflection coefficient, R = Λ Λ1 +Λ2 , for small velocity contrasts and small horizontal wavenumbers (i.e. small inci-
1 dence angles). It is interesting to notice that the difference between U1/2 and U1 is just the factor 2iΛ . However the phases of U1/2 and U1 are identical. The multi-one-way scheme, as described in this paper, mainly corrects the amplitudes. In the lower half space, z > z0 , the analytical solutions of Eqs. (39)–(41) are
U0 (z) =
exp (iΛ1 z0 + iΛ2 (z − z0 )) , 2iΛ1
U1/2 (z) = 0,
U1 = Ra
exp (iΛ1 z0 + iΛ2 (z − z0 )) . 2iΛ1
(45)
The total field, namely the transmitted field, after the interface at z = z0 is equal to U0 + U1 and we have U0 + U1 = (1 + Ra )U0 .
(46)
1 + Ra is an approximation of the transmission coefficient for small incident angles and small contrasts, since the plane wave transmission coefficient is 1 + R [22]. In this study the square root operator, Λ, and not its approximation has been used, which restricts the validity of this demonstration to small incidence angles. However this appendix demonstrates that the multi-one-way scheme allows estimating of the amplitudes of reflected and transmitted waves in a simple layered earth for small incidence angles and small velocity contrasts and explains why reflection data can be generated with this scheme.
References [1] [2] [3] [4]
W.A. Mulder, R.-E. Plessix, A comparison between one-way and two-way wave-equation migration, Geophysics 69 (6) (2004) 1491–1504. D. Lee, A.D. Pierce, Parabolic equation development in recent decade, J. Comput. Acoust. 3 (2) (1995) 95–173. J.F. Claerbout, Imaging the Earth’s Interior, Blackwell Scientific Publications, 1983. M.A. Leontovich, V.A. Fock, Solution of the problem of propagation of electromagnetic waves along the earth’s surface by the method of parabolic equation, J. Exp. Theoret. Phys. 16 (1946) 557–573. [5] F. Tappert, The parabolic approximation method, in: J.B. Keller, J.S. Papadakis (Eds.), Wave Propagation and Underwater Acoustics, American Institute of Physics, vol. 70, Springer, New York, 1977, pp. 224–280. [6] J.F. Claerbout, Coarse grid calculation of waves in inhomogeneous media with application to delineation of complicated seismic structure, Geophysics 35 (1970) 407–418.
D. Kiyashchenko et al. / Wave Motion 43 (2005) 99–115
115
[7] F. Collino, P. Joly, Splitting of operators, alternate directions, and paraxial approximations for three-dimensional wave equation, SIAM J. Sci. Comput. 16 (5) (1999) 1019–1048. [8] M. Kern, A non-split 3D migration algorithm, in: S. Hassanzadeh (Ed.), Mathematical Methods in Geophysical Imaging, SPIE, 1995. [9] D. Ristow, T. Ruhl, 3-D implicit finite-difference migration by multi-way splitting, Geophysics 62 (2) (1997) 554–567. [10] J. Gazdag, Wave equation migration with the phase shift method, Geophysics 43 (1978) 1342–1351. [11] J. Gazdag, P. Squazerro, Migration of seismic data by phase shift plus interpolation, Geophysics 48 (1984) 124–131. [12] R.-S. Wu, Wide-angle elastic wave one-way propagation in heterogeneous media and an elastic wave complex-screen method, J. Geophys. Res. 99 (B1) (1994) 751–766. [13] R.-S. Wu, L.-J. Huang, Reflected wave modeling in heterogeneous acoustic media using De Wolf approximation, in: S. Hassenzadeh (Eds.), Mathematical Methods in Geophysical Imaging, SPIE, 1995, pp. 176–186. [14] X.-B. Xie, R.-S. Wu, 3D elastic wave modeling using the complex screen method, in: Proceedings of the SEG 66th Annual Meeting, 1996. [15] S. Jin, R.-S. Wu, Prestack depth migration using a hybrid pseudo-screen propagator, in: Proceedings of the SEG 69th Annual Meeting, 1999. [16] S. Jin, Experimenting with the hybrid pseudo-screen migration, in: Proceedings of the SEG 70th Annual Meeting, 2000. [17] Yu. Zhang, G. Zhang, N. Bleistein, Wave equation migration arising from true amplitude one-way equations, Prob. Inverses 19 (2003) 1113–1138. [18] A. Bamberger, B. Engquist, L. Halpern, P. Joly, Parabolic wave equation approximations in heterogeneous media, SIAM J. Appl. Math. 48 (1) (1991) 99–128. [19] C. Wapenaar, J. Grimbergen, Reciprocity theorems for one-way wavefields, Geophys. J. Int. 127 (1996) 169–177. [20] M. De Hoop, A. Gautesen, Uniform asymptotic expansion of the square root Helmholtz operator and the one-way wave propagator, SIAM J. Appl. Math. 63 (3) (2003) 77–800. [21] F. Aminzadeh, J. Brac, T. Kunz, 3-D Salt and Overthrust Model, S.E.G., Tulsa, Oklahoma, 1997. [22] K. Aki, P.G. Richard, Quantitative Seismology, Theory and Methods, Freeman abd Company, San Francisco, 1980.