The design of a homotopy-based 1-D seismic FIR f − x wavefield extrapolation filters

Abstract This paper proposes a novel application of the homotopy method to design accurate non-causal complex-valued seismic FIR wavefield extrapolation digital filters for the purpose of seismic migration. The FIR filter design problem is solved using the scalar homotopy continuation method, since the filter design system of equations is overdetermined. In addition, the design approach was modified to satisfy the consistency conditions of the scalar homotopy method. Being a novel application to design FIR filters, the scalar homotopy method requires longer running design time compared to other existing methods such as the Weighted Least Squares (WLSQ) and the L1-norm methods. At the same time, the homotopy method is globally convergent and does not require the inversion of the Jacobian matrix, which is a significant advantage for practical purposes. Finally, the designed wavefield extrapolation digital filters were tested using the challenging Marmousi model data set, where prestack migration was performed successfully.


I. INTRODUCTION
Seismic migration is one of techniques for obtaining accurate images of inhomogeneous subsurface geological structures. One of the well established techniques of migration is the explicit depth wavefield extrapolation. This technique requires filtering seismic wavefields in the frequency-space (f − x) domain based on the following (normalized) one-dimensional (1-D) wavenumber response [1]: where k c is the normalized wavenumber cut-off, k x is the normalized horizontal wavenumber, and b is a constant that is equal to the ratio of the depth-sampling interval over the spatial sampling interval. Since, both the magnitude and the phase spectra are of even symmetry, therefore, (1) can be approximated by a non-causal complex-valued Finite Impulse Response (FIR) digital filter that is of even symmetry (see for example [1]- [8]). To obtain accurate seismic images, the accuracy of the FIR f − x digital filters requires a large number of coefficients. At the same time, short length filters are beneficial for strong lateral varying media.
To design these FIR filters, the homotopy continuation is used, which is a numerical method for calculating all solutions of a set of equations, particularly, nonlinear, [9], [10]. Such a method has found itself lending into many interesting applications within various disciplines (see for example [11]- [15]). The homotopy method involves, in principle, identifying a simple equation(s) with known solutions and slowly deforming it into a desired equation(s) with unknown solutions. A family of paths, during the deformation process, is defined from the known solutions of the simple equation to the desired solutions of the given equation(s). One way to track the solution paths numerically is through the so-called continuation method [10]. The homotopy method is known to be globally convergent in contrast to many other gradient based numerical methods [9], [10].
In this paper, it is proposed to design FIR f − x wavefield extrapolation digital filters that are accurate for approximating (1) using a homotopy method. This is considered a novel and new way of designing seismic migration FIR filters and can lend itself as a method to design other linear filters. Homotopy continuation methods requires solving system of equations that have equal number of unknowns to the number of equations. The problem of designing such filter requires solving an overdetermined problem, where the Jacobian matrix is not square, since the number of FIR filter coefficients are less than the number of wavenumber (frequency) samples. Therefore, it is suggested to use the scalar-based homotopy continuation approach suggested in [16], [17]. This method can be used for overdetermined systems of equations and does not require the inversion of the Jacobian matrix, while going through the algorithm iterations. Various designs of the seismic FIR wavefield extrapolation filters are shown and compared with designed filters via other algorithms such as the Weighted Least Squares (WLSQ) algorithm [6] and the L 1 -norm algorithm [7].
The main objective is to design an FIR f − x wavefield extrapolation filter that approximates (1). The transfer function for the symmetrical FIR f − x extrapolation filter is given by the noncausal transfer function: where M = N −1 2 and N is the length of the filter. Due to the even symmetry in both the magnitude and the phase responses in (1), the FIR coefficients h[n] belong to the set of complex numbers. Evaluating the FIR f −x extrapolation filter transfer function in (2) on the unit circle and performing few mathematical simplifications result in the following transfer function in the wavenumber k x domain: The design problem is to determine the set of the f − x FIR wavefield extrapolation filter coefficients h[n] such that H(e jkx ) given by (3) fits the ideal wavenumber response H(e jkx ) given by (1). This means, according to [7], that when sampling the wavenumber k x to obtain k xi (i = 1, 2, . . . , R), the response H(e kx i ) will exactly fit the desired H(e jkx ) at k xi (i = 1, 2, . . . , R). Hence, we can write the filter design problem as: or via the error function, where and, finally, Note that, in this case, Also, note that * denotes the Hermitian conjugate.

A. The Homotopy-based FIR f −x Extrapolation Filter Design Method
Let Γ(h, t) be a homotopy function, where t is a scalar called the homotopy parameter and t ∈ [0, 1] [9]. This function serves the objective of continuously deforming the function where h 0 is a simple solution. For example, it could be the least square solution of (4): Note that the homotopy function has to be continuous [9] so that the continuous deformation will start from:: and reach The homotopy method assumes that (5) at t = 0 is constructed in a way that the solutions are easy to obtain. Using numerical continuation, the parameter t is gradually varied until when t = 1 the homotopy function in (5) is equal to E(h).
There exists many forms of homotopy functions [9], [10]. Among the popular ones is the Newton homotopy function given as follows [9], [10]: Some computer codes have been developed for homotopy methods [13]. However, the seismic FIR extrapolation filter design problem in (5) represents an overdetermined system of equations, hence, cannot be solved using the homotopy function Γ(h, t). A possible way to overcome this challenge is to convert the vector based homotopy function in (5) into a scalar based homotopy continuation function, similar to the work shown in [16] and [17]. Note that (5) is equal to 0 and one can show that: Hence, we can transform the vector-based homotopy shown in (12) into a fictitious time dependent scalar function γ(h, t) given as [16], [17]: which holds for t ∈ [0, 1] since (13) is satisfied. Now, in order to guarantee that the homotopy path lies on the hyper-surface given by (14), the following consistency equation then can be shown to be: Assuming that the derivative of filter coefficients h with respect to t is in the direction of the gradient of (14), then a possible solution for such a derivative is given by [16]: in which the gradient of γ is assigned to be the driving force to adjust h. We can show that, from (14), the derivatives of the scaler homotopy function γ(h, t) with respect to h and t are given as follows, respectively: and Finally, employing an Euler forward scheme on (16) we can obtain [17]: where δt is homotopy scalar step size. This step size is suggested to be small to achieve robust designs. You can see from (19), that we do not require inverting the Jacobian matrix C, hence, this is really suitable for such an overdetermined problem.

B. The Proposed Homotopy-based FIR f − x Extrapolation Digital Filter Design Algorithm
The proposed Homotopy-based algorithm for designing the seismic migration f − x FIR digital filters is summarized as follows. The algorithm starts with an arbitrary complex-valued vector h 0 and choosing a filter length N , a wavenumber cutoff k c , and a suitable homotopy scalar step size δt. Then for the k th iteration, one will use the following steps: 1) Formulate the matrices C and H d based on (7) and (8), respectively. 2) Compute E(h 0 ) and E(h k ) using (5).
3) Obtain h k+1 using (19). The algorithm will stop once all the points within the path defined through the homotopy parameter are followed.

III. RESULTS
A design of the 1-D FIR wavefield extrapolation filter is presented here using the proposed homotopy method and is compared with those defined as a starting simple filter solution for the homotopy seismic migration filter design problem as well as filters designed using other methods. The results obtained in all the simulations shown in this paper are not restricted to the chosen parameters of the filters. For the wavenumber response figures, only half of the spectrum is shown, since the designed filters' spectra are of even symmetry with respect to the k x = 0 axis.
The FIR wavefield extrapolation filter was designed as shown in Figure 1 using the proposed algorithm. The filter was designed for a normalized wavenumber cut-off k c = 0.25, a homotopy step size δt = 10 −5 and with N = 35. The same figure shows the simple solution design which is basically the least square solution with an added 10% white noise. Figure 1 (a) shows the magnitude responses, which approximates the desired seismic wavefield extrapolation filter magnitude response for both the passband and the stopband (evanescent) parts (as given previously in (1)) with a very good design for such an application [1]. At the same time, Figure 1(b) shows the corresponding phase responses of the designed seismic FIR wavefield extrapolation filter that is very similar to the desired phase response based on (1) within its passband and stopband regions. Also, Figure 1 displays the paths of the homotopy function deforming the FIR wavefield extrapolation filter coefficients (c) real part and (d) imaginary part, respectively, of the initial seismic FIR wavefield extrapolation filter solution h 0 [n] into the final homotopy solution h[n] that were shown in Figure 1.
In addition, the designed seismic FIR wavefield extrapolation filter via homotopy was compared with two filters designed with the Weighted Least Squares (WLSQ) method reported in [6] and the L 1 -norm method reported in [7]. Both the homotopy and the L 1 -norm better approximate the desired wavenumber response, when compared to the WLSQ wavenumber response. In terms of design efficiency, both the L 1 -norm as well as the WLSQ methods are more efficient than the proposed homotopy method, as can be seen in Table  I. This is due to the small step size δt. If δt was selected to be larger, then the design run time is expected to be smaller but probably at the expense on the filter performance. The passband magnitude mean error is comparable for the homotopy and the L 1 -norm designs and are lower than that of the WLSQ. The L 1 provided the best passband phase mean error when compared with both the homotopy and the WLSQ, although the homotopy passband phase mean error is much lower than the WLSQ. Note that the running design time for the seismic FIR wavefield extrapolation filters is not the main concern for this application [2] because once the filters are designed, they will be stored in a look-up table, and the appropriate filter will be retrieved for usage according to the required wavenumber cut-off. It is worth investigation in future research to consider the effect of δt on the filter design performance.

IV. CONCLUSIONS
The paper presented a novel application of the homotopy continuation method for the design of the challenging 1-D complex-valued seismic FIR wavefield extrapolation digital filters. These filters can accurately be used to solve the wave equation that governs the propagation of seismic wavefields. The homotopy problem was converted from a vector-based into a scalar based problem to avoid the inversion of the Jacobian matrix, which is not square in the case of seismic FIR wavefield extrapolation filters. For the seismic extrapolation filtering problem, appropriate wavenumber responses were obtained at the expense of a longer running design time compared to the L 1 -norm and WLSQ methods. The additional design cost, however, resulted in better practical filters for this application when compared to that obtained using the WLSQ method and are comparable to the L 1 -norm method. This will help geophysicists deliver their migrated seismic data with a higher degree of accuracy using a new and novel method.

V. ACKNOWLEDGMENTS
The author thanks King Fahd University of Petroleum & Minerals (KFUPM) for sponsoring this work.