陈景波---中国科学院地质与地球物理研究所.pdf

GEOPHYSICAL RESEARCH LETTERS, VOL. 31, L06613, doi:10.1029/2004GL019429, 2004 Optimization approximation with separable variables for the one-way wave operator Jing-Bo Chen and Hong Liu Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, China Received 7 January 2004; revised 19 February 2004; accepted 24 February 2004; published 23 March 2004. [1] An optimization approximation with separable variables for the one-way wave operator is presented in this Letter. This new method approximates the oneway wave operator by products of functions in space variables and functions in wave number variables by means of optimization approximation with separable variables. This approximation enables us to use FFT algorithm which is independent of space variables while suffering no problem of branch points present INDEX TERMS: 0902 in the generalized-screen method. Exploration Geophysics: Computational methods, seismic; 0999 Exploration Geophysics: General or miscellaneous; 3230 Mathematical Geophysics: Numerical solutions. Citation: Chen, J.-B., and H. Liu (2004), Optimization approximation with separable variables for the one-way wave operator, Geophys. Res. Lett., 31, L06613, doi:10.1029/2004GL019429. 1. Introduction [2] Consider the 3-D acoustic one-way wave equation for upcoming wave in the frequency-space domain @U ð z; x; y; wÞ ¼i @z sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ w2 @ 2 @2 U ð z; x; y; wÞ; þ þ c2 @x2 @y2 ð1Þ where U is the wavefield, w is the angular frequency, x is the lateral coordinate along the in-line direction, y is the lateral coordinate along the cross-line direction and c = c(x, y, z) is the velocity function. [3] Now we consider the one-way propagator associated with equation (1), i.e., the solution of equation (1) with initial condition U ð z ¼ z0 ; x; y; wÞ ¼ dð x x0 ; y y0 Þ: For a sufficiently small vertical step Dz = z z0 (thin slab) and using the high-frequency approximation, the one-way thin slab propagator is given by de Hoop et al. [2000]: 1 gð z; x; y; z ; x ; y Þ ’ 4p2 0 0 0 Z " sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ # w2 exp i kx2 þ ky2 Dz 2 cðz; x; yÞ exp i kx ð x x0 Þ þ ky ð y y0 Þ dkx dky ; ð2Þ Copyright 2004 by the American Geophysical Union. 0094-8276/04/2004GL019429$05.00 where kx, ky are wave number and 1 z ¼ z0 þ Dz: 2 For laterally homogeneous thin slab, i.e., c(z, x, y) is independent of x, y, the propagator reduces to Gazdag’s phase-shift operator [Gazdag, 1978]. In this case, the computation of equation (2) requires only one two-dimensional FFT. For inhomogeneous thin slab, however, the computation of equation (2) requires one two-dimensional FFT for each different velocity c(z, x, y). This means a considerable computational effort. The split-step Fourier method introduced by Stoffa et al. [1990] requires much less computational cost by using a simple correction term applied in the w, x domain to deal with lateral velocity variations. But this approach only works well for smooth velocity variations and near vertical propagation angles. [4] Le Rousseau and de Hoop [2001] developed a scalar generalized-screen method which generalizes the phasescreen and the split-step Fourier methods to increase their accuracies with large and rapid lateral variations. Using two Taylor approximation and a perturbation hypothesis, this approach approximates the one-way wave operator by products of functions in space variables and functions in wave number variables. This approximation enables the dependency of equation (2) on x, y to be taken out of the integral thus resulting in a simplification of the computation. In spite of its great success, this method suffers a problem of branch points, and an integral contour deformation in the complex plane is needed. [5] Recently, Song [2001] (an English translation of this paper can be obtained from chenjb@mail.igcas.ac.cn) suggested a theoretical method of expressing a multi-variable real function by products of single-variable functions. In fact, Song’s method also works for complex functions. In this Letter, we will present a numerical implementation approach of Song’s method for complex functions. Based on the numerical approach, we will obtain an approximation of the one-way wave operator which attains the goal of the generalized-screen method but suffers no problem of branch points. [6] The new technique developed in this Letter can be used to construct fast 3-D wave-equation prestack depth migration algorithms. Now seismic imaging has been commonly applied to regions where geologic complexities are present. 3-D wave-equation prestack depth migration algorithms play a very important role in imaging regions with geologic complexities. However, because of the huge prestack data, the computational efficiency of the algorithms is in great demand. Therefore, the fast prestack depth migration algorithms based on this new technique L06613 1 of 4 CHEN AND LIU: OPTIMIZATION FOR ONE-WAY WAVE OPERATOR L06613 will be of great significance in imaging complex geologic regions. 2. Optimization Approximation With Separable Variables for the One-way Wave Operator We use repeated rectangle formula to integrate the left-hand side of equation (8), i.e, use rectangle formula on each rectangle [ui, ui+1] [kj, kj+1], and obtain Z bZ d Aðul ; k ÞAð~u; k Þ*fð~uÞ d~u dk a [7] We approximate the one-wave operator in the frequency-wave number domain c ¼ m X n Z uiþ1Z kjþ1 X ui i¼1 j¼1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Aðu; k Þ ¼ exp i u2 k 2 dz ; L06613 ð3Þ m n X X Aðul ; k ÞAð~u; k Þ*fð~uÞ d~u dk kj ! DuDkA ul ; kj A ui ; kj * fi : ð9Þ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ w where u = cð x;y;z , k = kx2 þ ky2 . Þ [8] The optimization approximation with separable variables for equation (3) is to find functions f(u), y(k) and a complex number l such that Set F = (f1, f2,. . ., fm)T. Using equations (8) and (9), we obtain the equation satisfied by F: ~f ~ ðuÞy ~ ðk Þ* kL2 ; k Aðu; k Þ lfðuÞyðk Þ* kL2 ¼ min k Aðu; k Þ l FF F ¼ jlj2 F : i¼1 j¼1 ~ y; ~ l ~ f; ð4Þ ð10Þ Here F is a matrix with entries: ~ C, and where * denotes the complex conjugate, l2 fli ¼ ~ ðuÞ kL2 ¼ 1 ; ~2 f ~ ðuÞ : f ~ ðuÞ 2 L2 ½a; b; k f f n X DuDkA ul ; kj A ui ; kj *; l; i ¼ 1; 2; . . . ; m: j¼1 ~2 y ~ ðk Þ : y ~ ðk Þ 2 L2 ½c; d ; k y ~ ðk Þ kL2 ¼ 1 : y Further, let A = (ai,j) be a matrix with entries: Using Lagrange multiplier, it can be easily proved that the solution to equation (4) is the eigenfunction corresponding to the eigenvalue with maximum modula of the following dual integral equation system ai;j ¼ A ui ; kj ; i ¼ 1; 2; . . . ; m; j ¼ 1; 2; . . . ; n: Then we have Z d F ¼ AAH ; ð11Þ Aðu; k Þyðk Þdk ¼ lfðuÞ; c ð5Þ Z b Aðu; kÞ*fðuÞdu ¼ l*yðk Þ: a In general, the analytical solution of equation (5) is not available, and the system (5) can only be solved numerically. To obtain the numerical solution, we transform equation (5) into the following two independent self-adjoint integral equations Z bZ d a ð6Þ c Z bZ d a Aðu; k ÞAð~ u; k Þ*fð~ uÞ d~ u dk ¼ jlj2 fðuÞ; Aðu; k Þ*A u; ~ k y ~ k du d ~ k ¼ jlj2 yðk Þ: ð7Þ c Now we use two-dimensional numerical integration to solve equations (6) and (7). First consider the integral equation (6). Consider partitions of intervals [a, b] and [c, d] with nodes: ui ¼ a þ ði 1ÞDu; i ¼ 1; 2; . . . ; m þ 1; kj ¼ c þ ð j 1ÞDk; j ¼ 1; 2; . . . ; n þ 1; ba ; m dc Dk ¼ : n Du ¼ where AH denotes the conjugate transposition of A. For simplicity, we have incorporated the area factor of the rectangular into the eigenvalue and keep the notation for the eigenvalue unchanged. [9] Now consider the integral equation (7). Set Y = (y1, y2,. . ., yn)T. In the same way as before, we obtain the equation satisfied by Y: a Aðul ; k ÞAð~ u; k Þ*fð~ uÞ d~ u dk ¼ jlj2 fl ; l ¼ 1; 2; . . . ; m: ð12Þ G ¼ AH A: ð13Þ where From equations (10) – (13), we can draw a conclusion that F and Y are the left and right singular vector of A corresponding to the maximum singular value l1 respectively. Using the power method (also called vector iteration method [Stoer and Bulirsch, 1993]), we can easily obtain F from equation (10), and then we have Y = AHF. Let f(1)(u) and y(1)(k) denote the interpolation function of F and Y respectively. [10] Now we obtain the optimization approximation with separable variables for A(u, k): Set f(ul) = fl, l = 1, 2,. . ., m + 1. From equation (6), we have Z bZ d Gy ¼ jlj2 Y; Aðu; k Þ ’ l1 fð1Þ ðuÞyð1Þ ðk Þ*: To increase accuracy, set c ð8Þ 2 of 4 A1 ðu; k Þ ¼ Aðu; k Þ l1 fð1Þ ðuÞyð1Þ ðk Þ*: CHEN AND LIU: OPTIMIZATION FOR ONE-WAY WAVE OPERATOR L06613 Figure 1. The first approximation of the one-way wave operator. The amplitude which varies with the wave number is shown. Figure 3. The convergence of (14). The amplitude which varies with the wave number is shown. [11] The approximation (14) satisfies We can obtain the optimization approximation with separable variables for A1(u, k) by using the same method as used for A(u, k): ð2Þ r X A ui ; kj ¼ ll fðlÞ ðui ÞyðlÞ kj *: ð2Þ Thus, we have the second order approximation Further, as r ! 1, we have Aðu; k Þ ’ l1 fð1Þ ðuÞyð1Þ ðk Þ* þ l2 fð2Þ ðuÞyð2Þ ðk Þ*: Aðu; k Þ ¼ Repeating this process, we finally obtain s X ll fðlÞ ðuÞyðlÞ ðk Þ*; ð15Þ l¼1 A1 ðu; k Þ ’ l2 f ðuÞy ðkÞ*: Aðu; k Þ ’ L06613 1 X ll fðlÞ ðuÞyðlÞ ðk Þ* ð16Þ l¼1 ð14Þ where s \ r and r is the rank of A. which converges at exponential rate [Song, 2001]. This expansion is closely related to the separable approximation of integral kernel [Pipkin, 1991]. Figure 2. The first approximation of the one-way wave operator. The phase which varies with the wave number is shown. Figure 4. The convergence of (14). The phase which varies with the wave number is shown. l¼1 3 of 4 CHEN AND LIU: OPTIMIZATION FOR ONE-WAY WAVE OPERATOR L06613 L06613 Table 1. Errors for Different s s=1 Error s=2 2 1.58439 10 s=3 3 1.53274 10 [12] Substituting equation (14) into equation (2), we obtain gð z; x; y; z0 ; x0 ; y0 Þ ’ 3.05616 10 3. Numerical Experiments [14] Now we perform some numerical experiments. We 40p 40p p , b = 1500 and c = 0, d = 25 . Set m = 40, n = 100. take a = 2500 In Figures 1 and 2, we show the comparison between l1f(1)(u)y(1)(k)* and A(u, k). We compare their amplitude (Figure 1) and phase (Figure 2) which vary with k for different u. We see that the first approximation l1f(1)(u)y(1)(k)* preserves the basic shape of A(u, k). The accuracy varies with u. This just demonstrates that l1f(1)(u)y(1)(k)* is the global approximation of A(u, k). [15] In Figure 3, we show the approximation (14) for different s at u38 for the amplitude. Figure 4 shows the 9.55636 105 corresponding results for the phase. At accuracy of 106, the iteration number of the power method is 3. Table 1 shows the following relative error for different s: Z s 1 X ðl Þ l f ð u Þ yðlÞ ðk Þ* l 4p2 l¼1 exp i kx ð x x0 Þ þ ky ð y y0 Þ dkx dky : ð17Þ The computation of equation (17) requires s two-dimensional FFTs which are independent of u and therefor of x, y. Usually for a small number s, the approximation (14) can achieve very good accuracy. This means a great simplification of equation (2) by using equation (17). Furthermore, there is no problem of branch points in equation (17). [13] Remark. In numerically solving equation (6), other high order numerical integration formulas can be chosen. The key point in choosing integration formulas is to insure that the resulting matrix F in equation (10) is an Hermitian matrix. For example, by direct calculations, it can be readily seen that applying the mid-point numerical integration formula results in an Hermitian matrix while applying the trapezoidal numerical integration formula does not lead to an Hermitian matrix. s=4 4 2 100 s X 1 X ðl Þ ðl Þ Error ¼ ll f ðu38 Þy ðki Þ* ; Aðu38 ; ki Þ T i¼1 l¼1 ð18Þ P 2 where T = 100 i¼1 jA(u38, ki)j . [16] For s = 4, the approximation (14) is already good enough and basically agrees with A(u, k). For different u, we can draw the same conclusion. Therefore, in the present example, the computation of equation (2) requires 40 FFTs while we only use 4 FFTs by using equation (17). [17] Acknowledgments. This work has been supported partially from Chinese Academy of Science with Key Project of Knowledge innovation KZCX1-SW-18 and Chinese National Scientific Foundation with Key Project 49894190. References de Hoop, M. V., J. H. Le Rousseau, and R.-S. Wu (2000), Generalization of the phase-screen approximation for the scattering of acoustic waves, Wave Motion, 31, 43 – 70. Gazdag, J. (1978), Wave equation migration with the phase-shift method, Geophysics, 43, 1342 – 1351. Le Rousseau, J. H., and M. V. de Hoop (2001), Modelling and imaging with the scalar generalized-screen algorithms in isotropic media, Geophysics, 66, 1551 – 1568. Pipkin, A. C. (1991), A Course on Integral Equations, Springer-Verlag, New York. Song, J. (2001), The optimization expression of functions and manifolds in high dimensions by ones in low dimensions (in Chinese), Chin. Sci. Bull., 46, 977 – 984. An English translation of this paper can be obtained at request to chenjb@mail.com.igcas.ac.cn. Stoer, J., and R. Bulirsch (1993), Introduction to Numerical Analysis, 2nd ed., Springer-Verlag, New York. Stoffa, P. L., J. T. Fokkema, R. M. de Luna Freir, and W. P. Kessinger (1990), Split-step Fourier migration, Geophysics, 55, 410 – 421. J.-B. Chen and H. Liu, Institute of Geology and Geophysics, Chinese Academy of Sciences, P.O. Box 9825, Beijing 100029, China. (chenjb@ mail.igcas.ac.cn; liuhong@mail.igcas.ac.cn) 4 of 4