# High-resolution seismic data regularization and wavefield separation

High-resolution seismic data regularization and wavefield separation Summary We present a new algorithm, non-equispaced fast antileakage Fourier transform (NFALFT), for irregularly sampled seismic data regularization. Synthetic tests from 1-D to 5-D show that the algorithm may efficiently remove leaked energy in the frequency wavenumber domain, and its corresponding regularization process is accurate and fast. Taking advantage of the NFALFT algorithm, we suggest a new method (wavefield separation) for the detection of the Earth's inner core shear wave with irregularly distributed seismic arrays or networks. All interfering seismic phases that propagate along the minor arc are removed from the time window around the PKJKP arrival. The NFALFT algorithm is developed for seismic data, but may also be used for other irregularly sampled temporal or spatial data processing. Fourier analysis, Computational seismology, Wave propagation 1 INTRODUCTION The sampling theorem was well developed nearly 70 yr ago (Shannon 1949), and has played a crucial role in signal processing and communication. If a function is frequency bandlimited or pre-filtered with the maximum angular frequency ωmax, it is completely determined by regularly sampled points spaced $$\Delta t = \frac{\pi }{{{\omega _{{\rm{max}}}}}}$$ seconds apart. For regularly sampled signals, Cooley & Tukey (1965) developed the revolutionary fast Fourier transform (FFT) algorithm to accurately and quickly conduct data analysis in the temporal frequency and spatial wavenumber domains. The study of sampling reached a very mature state and became mathematically oriented (e.g. Butzer 1983). In the modern Hilbert-space formulation, Shannon's sampling has been re-interpreted as an orthogonal projection onto the subspace of bandlimited functions with the sampling basis functions (e.g. sinc-function) orthogonal with each other (Unser 2000). However, for irregular non-equispaced sampling, Shannon's sampling theorem is not appropriate and the basis function orthogonality disappears. The fundamental requirement for the FFT is no longer satisfied, and thus the frequency spectrum has to be calculated directly from the original direct Fourier transform (DFT). The slowness of the DFT is not the fatal issue. It is energy leakage from one frequency into others that makes the calculated spectrum inaccurate without further manipulation (e.g. Wehlau & Leung 1964; Press et al.1993). The technique of data regularization is designed to rebuild the correct frequency spectrum and then transform the irregularly sampled data into regularly sampled data. Temporal spectral energy leakage due to the irregular time interval was documented early in astrophysical observations (Wehlau & Leung 1964; Gray & Desikachary 1973). In order to extract very weak frequency components, the technique of pre-whitening was adopted to rebuild the original spectrum in the time or frequency domain. It selects the largest amplitude frequency in the spectrum and subtracts the corresponding sinusoid from the data, then a new spectrum computation is followed free of this frequency and its leaked energy. The new spectrum is searched for the largest amplitude frequency again, and the process is repeated until the noise level is reached. In exploration seismology, due to natural, cultural and economic reasons, the effects of irregular and inadequate spatial sampling of seismic data are very important problem in seismic data processing. Both marine and land data, especially wide azimuth in five dimensions (e.g. shot xs, shot ys, reciever xr, receiver yr; time t), are usually quite irregular and sparse in at least two of the four spatial dimensions due to acquisition and processing costs (Trad 2009). However, seismic migration imaging theories require regular sampling in order to make the correct constructive and destructive interference from data (Abma et al.2007). Fourier finite difference-based wave equation migration explicitly requires regular sampling along each spatial dimension (Ristow & Ruhl 1997). In earthquake seismology, in order to improve the detection of nuclear tests worldwide, seismic arrays have been built, since 1950s to suppress noise and enhance the signal-to-noise ratio. Array apertures range from less than 3 km (e.g. NORES) to more than 100 km (e.g. NORSAR). Most of the geometries of seismic arrays are in certain patterns (e.g. circle, cross, rectangular, etc.), often with irregular spatial sampling. In the case of regional seismic networks, apertures can be larger than hundreds of kilometres (e.g. GRSN in Germany and F-net in Japan), and station distributions are intrinsically irregular. The energy leakage caused by spatial irregularity may affect the measurement of signal slowness and backazimuth in the frequency wavenumber domain with the seismic array/network data (Capon 1969; Kvarna & Doornbos 1986). In particular, if the signal is weaker than the dominant interfering phases, the leaked energy by the stronger phases may bury the weaker signal in the frequency–wavenumber domain (Wehlau & Leung 1964), even if there is no random background noise. Irregular and sparse spatial sampling causes not only energy leakage but also aliasing in the frequency domain. A variety of methods have been suggested for interpolation of spatial aliased data. For instance, Spitz (1991) utilizes the predictability of linear events to obtain an interpolation operator based on the spatial prediction filter in the f–x domain, Abma & Claebout (1995) attenuate random noise with a prediction filter in t–x and f–x domains, and Gulunay (2003) uses non-aliased lower temporal frequencies to de-alias higher frequencies in the f–k domain. However, these methods cannot handle irregular sampling. Regularization of irregularly and sparsely sampled seismic data is primarily conducted in the Fourier domain. The principle is to obtain the correct Fourier coefficients from the irregular data, and then the data regularization may be readily done by an inverse Fourier transform. Fourier reconstruction (e.g. Hindriks & Duijndam 2000; Zwartjes & Sachi 2007) and antileakage Fourier Transform (ALFT, Xu et al.2005, 2010) are two different categories of the method recently developed. Fourier reconstruction is an inverse process for each temporal frequency slice. Because the number of irregularly and sparsely sampled data is usually less than that of inverted Fourier coefficients globally or locally, this kind of inversion is basically ill-conditioned. Least-squares Fourier reconstruction (Hindriks & Duijndam 2000) treats the inversion as a mix-determined issue. It is difficult to find the optimal inverse parameter to balance the misfitting term and the model norm term, and the inverted model cannot exactly recover the original input data. Fourier reconstruction with sparse inversion (Zwartjes & Sacchi 2007) adopts the hybrid norm inversion to solve the aliasing issue using the non-aliased lower frequencies as guide, but the introduction of two inverse parameters makes the inversion even harder in practice. Damped least-norm Fourier inversion (Jin 2010) treats the inversion as a pure underdetermined inverse problem. No inverse parameter is required, and it may completely recover the original input data in theory, but a damping matrix is required to be as close as possible to the final Fourier coefficient model. The ALFT algorithm was first proposed for seismic data regularization by Xu & Pham (2004). In principle, it is similar to the technique of pre-whitening first suggested in astrophysics (Wehlau & Leung 1964; Gray & Desikachary 1973), but ALFT has been implemented for multidimensional seismic data regularizations. It may be used to process 2-D/3-D common offset and common azimuth data (Xu & Pham 2004), and 5-D seismic data by the cascading two passes of 3-D common shot and common receiver regularizations (Xu et al.2005). The accuracy of ALFT regularization is closely related to the size of spatial subwindows. The larger the subwindow, the denser the sampling in the spatial spectrum, and thus the energy leakage can be better removed. However, the large window not only significantly increases computational cost but also makes de-aliasing very difficult. So far, for both Fourier reconstruction and ALFT methods, the most efficient way to handle aliasing is to adopt small spatial subwindows and thus the events inside are basically regarded as linear (Schonewille et al.2009; Jin 2010). In order to solve the conflict between accuracy and computational efficiency for ALFT regularization, Xu et al. (2010) developed a new version of the ALFT algorithm. They introduced a single model parameter inversion method to directly increase the sampling number in the wavenumber domain without changing the spatial input data subwindow. High-resolution accurate regularization for irregularly sampled data is very important not only in exploration seismology and astrophysics, but also in earthquake seismology. Seismic arrays and regional seismic networks are designed to record wavefields propagating along different directions. If multiple wavefields arrive at nearly the same time, they interfere with each other, and thus complicate their separation in the time–space domain, making it difficult to conduct the slowness and backazimuth analysis for individual wavefields based on seismic records. Even if a single wavefield is recorded without noise, the irregular distribution of stations with respect to epicentral distance is still a problem. The resulting energy leakage in the frequency–wavenumber spectrum may hinder us from correctly constraining the slowness and backazimuth of the signals, especially when arrays are small. In order to detect weak signals from background noise, a number of stacking methods (e.g. Array beamforming, Nth root process, Weighted stack methods and Double beam technique) have been developed for seismic array and network data analysis (Muirhead & Datt 1976; Kruger et al.1993; Schimmel & Paulssen 1997; Kennett 2000). The fundamental assumption for these methods is that in the individual traces the noise around the signal is incoherent and no other coherent interfering energy arrives nearby. The stacking technique may significantly improve the signal-to-noise ratio, and plays an important role in nuclear test monitoring and the analysis of subtle anomalies at the core–mantle boundary (Kruger et al.1993) and the constraint of the density contrast at the inner core boundary (Koper & Dombrovskaya 2005). Among all of the weak signals, detection of the Earth's inner core shear wave phase PKJKP, the direct evidence of the solidity of the inner core, is the most challenging (e.g. Bullen 1951; Shearer et al.2011). Based on the IASP91 velocity model (Kennett & Engdahl 1991) and the high-resolution direct solution method (DSM, Takeuchi et al.1996; Kawei et al.2006), synthetic simulation shows that the PKJKP energy is so weak that it is well below the background noise level. The quasi-sphericity of the Earth makes the observation of PKJKP even harder. When the weak PKJKP travels along the major arc of the Earth's big circle to a seismic array or network, a number of other phases arrive at the same time along the opposite minor arc. These interfering phases are much weaker than the first arrival PKIKP and usually unidentified, but they are coherent and stronger than the inner core shear wave phase PKJKP. Current seismic stacking methods cannot handle this kind of situation. This may be one of the main reasons why the observation of PKJKP is so difficult even with high-quality broad-band seismic arrays or networks. In this paper, we present a new method for high-resolution and high-speed seismic data regularization. We first introduce the non-equispaced FFT (NFFT), ALFT algorithms, and then illustrate 1-D/2-D spatial and 5-D temporal and spatial synthetic tests. In the end, we suggest a possible wavefield separation method for the seismic array observation of the inner core shear wave phase PKJKP. 2 NON-EQUISPACED FAST FOURIER TRANSFORM For irregularly sampled seismic data, we may use the non-equispaced DFT (NDFT) to conduct the regularization, but the computational cost is heavy (if not infeasible), especially in the case of multidimensional data processing. Here, we briefly introduce the NDFT and then focus on the flexible and efficient NFFT (Potts et al.2001). Similar to the FFT, the NFFT is another way to fast calculate the frequency spectrum. Their difference is that the FFT works only for the regularly sampled data, and the NFFT can accurately and efficiently handle irregularly sampled data. The data sampling range of dimension $$d \in \mathbb{N}\$$is given by the torus   \begin{eqnarray} &&{{\mathbb{T}^d}:\ = \Big\{ {\boldsymbol{x}} = {{({x_t})}_{t = 0, \ldots ,d - 1}} \in {\mathbb{R}^d}:}\nonumber\\ &&{\quad\qquad-\, \frac{1}{2} \le {x_t} < \frac{1}{2},\ t = 0, \ldots ,d - 1 \Big\}} \end{eqnarray} (1)the sampling set is given by   $$\mathcal{X}: = \left\{ {{{\boldsymbol{x}}_j} \in {\mathbb{T}^d}:j\ = \ 0, \ldots ,M - 1} \right\}$$ (2)and the frequencies $${\boldsymbol{k}} \in {\mathbb{Z}^d}$$ are in the multi-index set   \begin{eqnarray} &&{{I_N}: = \Big\{ {\boldsymbol{k\ }} = {{\left( {{k_t}} \right)}_{t = 0, \ldots ,d - 1}}\ \in {\mathbb{Z}^d}:}\nonumber\\ &&{\qquad\quad- \frac{{{N_t}}}{2} \le {k_t} < \frac{{{N_t}}}{2},\ t\ = \ 0, \ldots ,d - 1 \Big\}} \end{eqnarray} (3)where $$N\ = {({N_t})_{t = 0, \ldots ,d - 1\! }},\ {N_t} \in 2\mathbb{N}$$ is the even multiband limit. For a finite number of given Fourier coefficients, $${f_{\boldsymbol{k}}} \in \mathbb{C},{\boldsymbol{\ k}} \in {I_N}$$, the forward NDFT may be expressed as   $${f_{\boldsymbol{k}}} = \mathop \sum \nolimits_{j = 0}^{M - 1} w\left( {{{\boldsymbol{x}}_j}} \right){f_j}{e^{ - 2\pi i{\boldsymbol{k}}{{\boldsymbol{x}}_{\boldsymbol{j}}}}}$$ (4)where w(xj) is the integral weight and $$\mathop \sum \nolimits_{j\ = \ 0}^{M - 1} w\ ( {{{\boldsymbol{x}}_j}} ) = 1.0$$. It takes $$\mathcal{O}( {M\,{\rm{N}}} )$$ arithmetical operations. The core of the ALFT technique is to clean leaked energy in the frequency domain, thus the forward NFFT, which transforms the irregularly sampled data from spatial domain to wavenumber domain, plays the key role in accuracy and speed of the ALFT regularization. Here, we summarize the forward NFFT algorithm as follows: $$Input{:}\ d,\ M \in \mathbb{N},\ N \in 2{N^d}$$  \begin{eqnarray} &&{{{\boldsymbol{x}}_j} \in {\left[ { -\frac{1}{2},\frac{1}{2}} \right]^d},\ j\ = \ 0, \ldots ,M - 1,\quad and\quad {f_j} \in \mathbb{C},}\nonumber\\ &&{\quad j\ = \ 0, \ldots ,M - 1} \end{eqnarray} (5) 1: For l ∈ In compute  $${g_{\boldsymbol{l}}}: = \sum\limits_{j \in {I_{n,m}}\left( l \right)} {f_j}\psi \left( {{{\boldsymbol{x}}_j} - {{\boldsymbol{n}}^{ - 1}} \odot {\boldsymbol{l}}} \right)$$ (6) 2: For k ∈ In compute d dimensional forward FFT  $${g_{\boldsymbol{k}}}: = \sum\limits_{{\boldsymbol{l}} \in {I_n}} {g_l}{e^{ - 2\pi i{\boldsymbol{k}}\left( {{{\boldsymbol{n}}^{ - 1}} \odot {\boldsymbol{l}}} \right)}}$$ (7) 3: For k ∈ IN compute  $${f_{\boldsymbol{k}}}: = \frac{{{g_{\boldsymbol{k}}}}}{{\left| {{I_n}} \right|{c_k}\left( \varphi \right)}}$$ (8) output:  fk, k ∈ IN $$\mathcal{O}( {| N |log| N | + M} )\ arithmetic\ operations$$ where xj is the normalized input sampling range, ψ is a one periodic window function, ck is the Fourier coefficients of the periodic window function and  In, m(l): ={j = 0, …, M − 1: l − m ≤ n⊙xj ≤ l + m}. 3 ANTILEAKAGE FOURIER TRANSFORM ALGORITHM Two versions of ALFT algorithms were developed for the multidimensional regularization of irregularly sampled seismic data. The first ALFT (Xu et al.2005) is similar to the pre-whitening technique (Wehlau & Leung 1964; Gray & Desikachary 1973) (hereafter it is referred as ALFT05), and the second ALFT is a one-parameter inversion related technique (Xu et al.2010) (hereafter it is referred as ALFT10). Here, we introduce both algorithms for theoretical integrity and convenience in the following discussion in this paper. ALFT05 Set all initial Fourier components to zero. Calculate all Fourier coefficients of original/residual input data with the NDFT (eq. 4). Find the wavenumber kmaxand its Fourier coefficient$${f_{{{\boldsymbol{k}}_{\rm {max}}}}}$$ with the maximum energy, and keep adding the Fourier coefficient$${f_{{{\boldsymbol{k}}_{\rm {max}}}}}$$to the initial Fourier components. Calculate the wavenumber kmaxcomponent in the original/residual input data   \begin{equation*}{f^{{{\boldsymbol{k}}_{\rm {max}}}}}( {{x_j}} ) = {f_{{{\boldsymbol{k}}_{\rm {max}}}}}{e^{2\pi i{\boldsymbol{k}}{{\boldsymbol{x}}_j}}}\end{equation*} Update the input by subtracting the wavenumber kmax component   \begin{equation*}{f^u_j} = {f_j}\ - {f^{{{\boldsymbol{k}}_{\rm {max}}}}}( {{{\boldsymbol{x}}_j}} ),\quad j\ = \ 0, \ldots ,M - 1\end{equation*} Iterate steps (b) to (e) until the leaked energy disappear or is under a threshold. ALFT10 Similar to ALFT05, Fourier coefficients are instead calculated with the inversion by minimizing the objective function   $$E\ = \sum\limits_{j = 1}^{M - 1} {\left\| {{f_j} - {f_{{\boldsymbol{k'}}}}{e^{2\pi i{\boldsymbol{k'}}{{\boldsymbol{x}}_j}}}} \right\|}$$ (9)where k' is the oversampled wavenumber in spatial frequency domain. ALFT05 is an iterative technique. The computation can be extremely heavy for the entire process of regularization with the NDFT. The adoption of the NFFT may significantly reduce the computation time, but the accuracy of the resulting spectrum is still a concern. In other words, ALFT05 cannot effectively remove the leaked energy in the spatial spectrum, especially when small spatial subwindows are selected in order to handle the sampling alias. ALFT10 is designed to improve the spectrum accuracy by significantly oversampling the wavenumber k' N times (e.g. N = 32) in eq. (9). It means that the new sampling rate Δk'is N times denser than in an equivalent DFT case. The adoption of oversampled wavenumber k' may help remove the leaked energy in the spectrum domain. However, what this inversion-oriented ALFT10 algorithm obtains in the gain of accuracy, it loses in speed, because the inversion takes N2 more time than the NDFT in ALFT05 to calculate the individual frequency spectrum. 4 NON-EQUISPACED FAST ANTILEAKAGE FOURIER TRANSFORM ALGORITHM Here, we suggest a new algorithm for the high accuracy and high-speed regularization of irregularly sampled data based on NFFT and ALFT05. Instead of directly oversampling the wavenumber in the frequency domain as ALFT10 does, we introduce the d dimensional input sampling parameter vector $${\boldsymbol{\sigma }} \in \mathbb{N}$$ and change the NFFT input sampling set (eq. 5) as   $${\boldsymbol{x}}_j^{\prime} = \ {{\boldsymbol{x}}_j} \odot \frac{1}{{\boldsymbol{\sigma }}},\quad {\boldsymbol{x}}_j^{\boldsymbol{=}} \in {\left[ { - \frac{1}{2},\frac{1}{2}} \right]^d}\quad j\ = \ 0, \ldots ,M - 1$$ (11) After being divided by σ, the original sampling set xj and the new sampling set $${\boldsymbol{x}}_j^{\prime}$$ are still in the normalized ranges with the same number of sampling input. However, the difference is that $$x_j^{\prime}$$ is now equivalent to a subset of zero-padded new set of size σ · M. The introduction of the input sampling parameter σ increases the frequency sampling rate by (σi − 1) times in ith dimension, and thus the resulting new frequency resolution is $$\mathop \prod \nolimits_{i\ = \ 0}^{d - 1} ({\sigma _i} - 1)$$ times higher. This means that the frequency of the maximum energy in each iteration now has a much greater chance to locate at or very close to a frequency sampling point in the spectrum. Therefore, we may not only remove the leaked energy more efficiently in smaller number of iterations, but also obtain a more accurate spectrum in the end. NFFT is a fast and good approximation to NDFT with a relative error usually less than 0.001 (Potts et al.2001). Therefore, in practice, we always use NFFT rather than NDFT. For the high-dimensional (especially 5-D) regularization, in order to significantly reduce the computational memory and time cost, we suggest a modification of the original NFFT algorithm. Given the fact that the input sampling locations are always the same for different temporal frequency slices and different iterations in a subwindow, we just need to calculate those one-periodic window functions in eq. (6) one time and then save them inside a shared memory for the parallel computation. For convenience in the following tests and discussion, hereafter we refer to our new algorithm as the non-equispaced fast ALFT (NFALFT). 5 SYNTHETIC TESTS In order to quantitatively examine the performance of the NFALFT algorithm, three sets of analytical function are adopted to calculate the regularly and irregularly sampled synthetic data, respectively. Both ALFT05 and NFALFT regularizations are conducted for the irregularly sampled synthetic data. Then, we compare the regularized data with the regularly calculated data from these analytical functions. We start with a 1-D sinusoidal function y = 5.0*cos (0.5*x) + 2.5*cos(x). If the range of the x window is infinite, there are supposed to be two δ-functions in the spatial frequency (i.e. wavenumber) spectrum. One is at frequency 0.5 rad m−1 with amplitude 5.0, and the other is at frequency 1.0 rad m−1 with amplitude 2.5. If the range of the x window is limited and the sampling is regular, there are supposed to be two limit-sized pulses in the frequency spectrum with smaller amplitudes. If the range of the x window is limited and the sampling is irregular (randomly sampled 64 points in the x range from 0 to 31) (Fig. 1), the frequency spectrum is much more complicated and energy is leaked across the whole frequency range (the green line in Fig. 2a). After applying the antileakage processing with ALFT05 algorithm, most of the leaked energy is removed (the red line in Fig. 2a), but the resulting regularized results (the red line in Fig. 2b) do not match the analytical function (the blue line in Fig. 2b) well. This means that the ALFT05 algorithm cannot handle the leaked energy accurately in the frequency spectrum. Figure 1. View largeDownload slide The analytical function y = 5.0*cos (0.5*x) + 2.5*cos(x) (in red line) and the random sampling (in black cross). Figure 1. View largeDownload slide The analytical function y = 5.0*cos (0.5*x) + 2.5*cos(x) (in red line) and the random sampling (in black cross). Figure 2. View largeDownload slide (a) Frequency amplitude spectrum with NDFT (in green) and ALFT05 processing (in red). (b) Regularization with ALFT05 processing (in red) and the original analytical function (in blue). Figure 2. View largeDownload slide (a) Frequency amplitude spectrum with NDFT (in green) and ALFT05 processing (in red). (b) Regularization with ALFT05 processing (in red) and the original analytical function (in blue). Using the NFALFT algorithm, we take the spatial input sampling parameter σ = 32 to oversample the wavenumber k up to 31 times in the frequency domain. After applying the antileakage processing with NFALFT, just a number of discrete frequency with non-zero energy are left (the red line in Fig. 3a). The largest four amplitudes are 4.78 at k = 1.00 rad m−1, 2.44 at k = 0.50 rad m−1, 0.37 at k = 0.48 rad m−1 and 0.04 at 0.37 rad m−1, respectively. Energy at other frequencies are zeros or very close to zeros. The resulting regularized results by NFALFT are identical to the calculated ones from the analytical function (Fig. 3b). The thresholds for both ALFT05 and NFALFT antileakage processing are set the same (0.001). The NFALFT regularization reaches the threshold in 10 iterations, on the other hand, the ALFT05 algorithm took 1093 iterations to reach the threshold. The NFALFT is more than 2 times faster than the ALFT05, even though the frequency is 32 times oversampled relative to the ALFT05. Figure 3. View largeDownload slide (a) Frequency amplitude spectrum with NDFT (in green) and NFALFT processing (in red). (b) Regularization after NFALFT processing (in red) and the original analytical function (in blue). The two curves overlap with each other and thus only the red one is visible. Figure 3. View largeDownload slide (a) Frequency amplitude spectrum with NDFT (in green) and NFALFT processing (in red). (b) Regularization after NFALFT processing (in red) and the original analytical function (in blue). The two curves overlap with each other and thus only the red one is visible. The energy leakage in the frequency domain caused by the spatial sampling irregularity and the principal of antileakage algorithms has been demonstrated in the above 1-D sinusoidal synthetic test. In the next two tests, we will focus on the comparison of regularization with ALFT05 and AFALFT algorithms. In the second synthetic test, we use the 2-D sinusoidal analytical function z = sin (x) *cos(y). Its contour map and corresponding 102 × 102 regular sampling grids are presented in Fig. 4. We randomly remove 50 per cent of the grid to create a sparser and irregular sampled grid (Fig. 5a). Figure 4. View largeDownload slide (a) Contour map of the 2-D analytical function z = sin (x) *cos (y). (b) Regular sampling with 102 × 102 gridpoints for the analytical function. Figure 4. View largeDownload slide (a) Contour map of the 2-D analytical function z = sin (x) *cos (y). (b) Regular sampling with 102 × 102 gridpoints for the analytical function. Figure 5. View largeDownload slide (a) Irregular sampling by randomly removing 50 per cent of the original sampling points. (b) Contour map of the 2-D regularization result with the ALFT05 algorithm. Figure 5. View largeDownload slide (a) Irregular sampling by randomly removing 50 per cent of the original sampling points. (b) Contour map of the 2-D regularization result with the ALFT05 algorithm. For the NFALFT algorithm, we take the spatial input sampling parameter σ  = (4, 4) to oversample the wavenumber k up to 15 times in the frequency domain. After applying the ALFT05 algorithm, the contour map plotted with the regularized data (Fig. 5b) does not match the original one from the analytical data (Fig. 4a). The distortion at the boundaries is significant. The regularization results look better in the middle, but small kinks are distributed along each of the contour lines. After applying the NFALFT algorithm, we recover nearly the same contour map as with the regularized data (Fig. 6). There are slight disconnections at contour lines of value 0. The thresholds for both ALFT05 and NFALFT antileakage processing are set as the same (0.001). The NFALFT regularization reaches the threshold in 77 iterations, on the other hand, the ALFT05 algorithm took 2733 iterations to reach the threshold. NFALFT is more than 14 times faster than ALFT05, even though the frequency is 16 times oversampled relative to ALFT05. Figure 6. View largeDownload slide Contour map of the 2-D regularization result with the NFALFT algorithm. Figure 6. View largeDownload slide Contour map of the 2-D regularization result with the NFALFT algorithm. Our third synthetic test is conducted with three 5-D hyperplane analytical functions.   $$A\!:3{x_1} + 2{x_2} - 2{x_3} - 2{x_4} - 25000t + 11600\ = \ 0$$ (12)  $$B\!:3{x_1} + 2{x_2} + 8{x_3} + 8{x_4} - 25000t + 2100\ = \ 0\$$ (13)  $$C\!:3{x_1} + 2{x_2} - 4{x_3} - 4{x_4} - 25000t + 14\ 000\ = \ 0$$ (14)where x1, x2, x3, and x4 are four spatial coordinates, which may be shot x, shot y, receiver x and receiver y, or offset x, offset y, cdp x and cdp y in exploration seismology and t is the temporal coordinate. For a regular 16 × 16 × 16 × 16  sampling in all four spatial dimensions, we show the sampling grid in the x3x4 plane in Fig. 7(a) and the corresponding profile of the 5-D analytical hyperplane in Fig. 8(a). In order to get irregular 16 × 16 × 16 × 16 sampling (Fig. 7b), we shift the original sampling points randomly in each spatial dimension. For the NFALFT algorithm, we take the spatial input sampling parameter σ = (2, 2, 2, 2) to oversample the wavenumber k up to 15 times in the frequency domain. After applying the NFALFT algorithm, the regularized hyperplanes (Fig. 8b) are basically the same as the ones calculated from the analytical functions with the regular sampling grids (Fig. 8c). The computation of the ALFT05 algorithm does not converge to threshold (0.001), and so here we do not provide a quantitative comparison of the computational cost between ALFT05 and NFALFT for this 5-D test. Figure 7. View largeDownload slide (a) The regular sampling with 16 × 16 gridpoints in x3x4 plane for the 5-D hypeplane function in eqs (12)–(14). (b) The random sampling in x3x4 plane in the 5-D hypeplane function. Figure 7. View largeDownload slide (a) The regular sampling with 16 × 16 gridpoints in x3x4 plane for the 5-D hypeplane function in eqs (12)–(14). (b) The random sampling in x3x4 plane in the 5-D hypeplane function. Figure 8. View largeDownload slide Comparison of original data and regularized data in x4 − t plane. (a) Original data calculated with eqs (12)–(14). (b) Regularized data by means of the NFALFT algorithm. (c) Amplitude difference between the two panels of (a) and (b). Figure 8. View largeDownload slide Comparison of original data and regularized data in x4 − t plane. (a) Original data calculated with eqs (12)–(14). (b) Regularized data by means of the NFALFT algorithm. (c) Amplitude difference between the two panels of (a) and (b). 6 SEISMIC WAVEFIELD SEPARATION Strong earthquakes (e.g. Mw > 6.5) usually trigger seismic wave propagation around the Earth along the minor and the major arcs, respectively. This is good for the global/regional waveform tomography (e.g. Li & Romanowicz 1996) with denser ray path coverage, but it is a challenge for seismic array observation. The inner core shear wave PKJKP is a unique phase to verify the solidity of the Earth's centre. However, due to the high attenuation in the inner core and the poor P-to-S and S-to-P conversions at the inner core boundary, PKJKP is so weak that it only barely stands out of the background noise (e.g. Shearer et al.2011). PKJKP travels along the major arc. Within its theoretical arrival time window, there are always a number of unidentified weak phases that arrive along the minor arc at nearly the same time. These minor arc interference phases are much weaker than the first arrival PKIKP, but they are still strong in comparison with the potential PKJKP, making the observation of PKJKP very difficult with seismic array/network stacking methods in the temporal–spatial domain. Given the fact that PKJKP and the interference phases travel along opposite directions in a great circle of the Earth, one should be able to separate seismic wavefields along the major arc from those along the minor arc in the temporal frequency and spatial wavenumber f–k domain. PKJKP appears in the negative wavenumber spectrum and those minor arc interference phases are in the positive wavenumber spectrum. However, in practice, due to the spatial irregularity of station distribution in arrays/networks, the energy leakage of strong minor arc wavefields into weak major arc wavefields in the f–k spectrum makes the wavefield separation difficult without further manipulation. In the following PKJKP synthetic analyses, we first show the energy leakage in the f–k domain caused by the irregular spatial distribution of a seismic array, and then demonstrate the efficiency of the NFALFT algorithm in the removal of leaked energy and wavefield separation. In order to obtain an irregularly distributed seismic array, in the epicentral distance range from 138° to 140°, we first set up a linear array of 40 stations with a regular station interval 5.55 km along the backazimuthal direction of 90°, and then randomly shift each station from the second to the 39th within the range of ± 2.6 km along the array direction. Synthetic seismograms (Fig. 9a) are computed with the high-resolution DSM (Kawei et al.2006) for an event Mw = 7.0 with respect to the IASP91 velocity model (Kennett & Engdahl 1991). The high corner frequency is 0.5 Hz, and the inner core shear quality factor is taken as 350, which is much larger than 85 first suggested in the PREM model (Dziewonski & Anderson 1981). Around the theoretical arrival time of PKJKP, we cannot identify any coherent seismic phase dipping from right to left, even after the data regularization (Fig. 9b). PKJKP is buried in the wavefields from the minor arc. Figure 9. View largeDownload slide Synthetic seismograms computed with the direct solution method (DSM) for the inner core shear wave phase PKJKP in the epicentral distance range of 138° − 140°. The theoretical arrival times of PKJKP are ∼1710 s (red dashed lines). (a) Seismograms from irregularly distributed 40 stations. (b) Regularized seismograms of 40 stations with regular interval of 5.55 km by means of the NFALFT algorithm. Figure 9. View largeDownload slide Synthetic seismograms computed with the direct solution method (DSM) for the inner core shear wave phase PKJKP in the epicentral distance range of 138° − 140°. The theoretical arrival times of PKJKP are ∼1710 s (red dashed lines). (a) Seismograms from irregularly distributed 40 stations. (b) Regularized seismograms of 40 stations with regular interval of 5.55 km by means of the NFALFT algorithm. In the f–k domain, the energy propagating along the minor arc is dominant in the positive wavenumber region, meanwhile, the energy leakage from the positive wavenumber into the negative wavenumber is also large due to the irregular spatial sampling (Fig. 10a). The leaked energy is so strong that it not only severely distorts the PKJKP waveforms but also make the PKJKP amplitude much weaker (Fig. 11a). After applying the NFALFT algorithm, we remove the leaked energy in the f–k domain (Fig. 10b). The wavefields along the major arc are well separated from the wavefields along the minor arc in both f–k domain (Fig. 10b) and x–t domain (Figs 11b and 12b). After this wavefield separation, we clearly see the PKJKP phase in individual traces in the major arc wavefield without stacking (Fig. 11b). There are also a number of other major arc phases. Two of them around relative time 1730 s, we identify as pPKJKP + PKPPKPdf. There are more minor arc phases in the time window (Fig. 12b). It is difficult to identify them unambiguously, but one of them arrives at nearly the same time as PKJKP does but with larger amplitudes. Figure 10. View largeDownload slide The temporal frequency versus spatial wavenumber f–k spectra. (a) The spectrum calculated for the irregularly distributed seismograms (Fig. 9a) with the NDFT method. The energy leakage is remarkable at each temporal frequency slice. (b) Removal of the leaked energy by the NFALFT algorithm. All seismic energy propagating along the minor arc is on the right, and all energy propagating along the major arc is on the left. Figure 10. View largeDownload slide The temporal frequency versus spatial wavenumber f–k spectra. (a) The spectrum calculated for the irregularly distributed seismograms (Fig. 9a) with the NDFT method. The energy leakage is remarkable at each temporal frequency slice. (b) Removal of the leaked energy by the NFALFT algorithm. All seismic energy propagating along the minor arc is on the right, and all energy propagating along the major arc is on the left. Figure 11. View largeDownload slide Seismic phases propagating along the major arc after the wavefield separation. Dashed red lines are the theoretical PKJKP arrival time. (a) Waveform distortion caused by the leaked energy from the minor arc wave energy. (b) More coherent and stronger PKJKP waveforms after the removal of the leaked energy by the NFALFT algorithm. Figure 11. View largeDownload slide Seismic phases propagating along the major arc after the wavefield separation. Dashed red lines are the theoretical PKJKP arrival time. (a) Waveform distortion caused by the leaked energy from the minor arc wave energy. (b) More coherent and stronger PKJKP waveforms after the removal of the leaked energy by the NFALFT algorithm. Figure 12. View largeDownload slide Seismic phases propagating along the minor arc after the wavefield separation. (a) The leaked energy also affects the strong minor arc wavefields. (b) The minor arc wavefields after the NFALFT regularization. One phase arrives at nearly the same time as PKJKP does (red dashed lines). Figure 12. View largeDownload slide Seismic phases propagating along the minor arc after the wavefield separation. (a) The leaked energy also affects the strong minor arc wavefields. (b) The minor arc wavefields after the NFALFT regularization. One phase arrives at nearly the same time as PKJKP does (red dashed lines). 7 DISCUSSION NFALFT and ALFT05 are iteration-based antileakage algorithms. Their computational costs mainly depend on the number of iterations for them to reach the threshold (e.g. 0.001). The efficiency of the leaked energy removal in each iteration is closely related to the frequency sampling rate in each dimension. The denser the sampling rate, the more likely the leaked energy of the dominant frequency may be removed in a single iteration. Otherwise, the iteration process can take a long time before reaching the threshold, and meanwhile the resulting accuracy of the frequency spectrum is poor. ALFT10 is an inversion-based antileakage algorithm to remove the leaked energy in the frequency spectrum. It directly increases the sampling rate in each dimension of the frequency spectrum and conducts a one-parameter least-squares inversion for each frequency in the spectrum. Thus, the computational cost of ALFT10 is $${{\boldsymbol \sigma }^{\boldsymbol 2}}$$ greater than that of NDFT. On the other hand, NFALFT takes the NFFT to compute the frequency spectrum, and thus the arithmetical operation contrast between ALFT10 and NFALFT is approximately $$\mathcal{O}( {{\sigma ^2}M \cdot {\rm{N}}} )$$ to $$\mathcal{O}( {| {\sigma \cdot N} |{\rm{log}}| {\sigma \cdot N} | + M} )$$. If the oversampling rates are the same for NFALFT and ALFT10, their accuracies are the same in theory. In practice, however, NFALFT may be more accurate because it is based on the NFFT, which provides a more accurate frequency spectrum than the inversion method does. The iteration number is usually less than 10 per cent of the total number of frequencies in the spectrum for NFALFT. Especially, for the linear events in high-dimensional data, the corresponding frequency spectrum is very sparse, and so the iteration number of NFALFT is even smaller. The iteration threshold is the only parameter we need to specify for the NFALFT algorithm, and usually we take it as a constant (e.g. 0.001) for all regularization. In comparison with the Fourier reconstruction inversion methods (e.g. Zwartjes & Sacchi 2007; Jin 2010), there are not trial-and-error inverse parameters or pre-conditional starting models required. For the high-dimensional irregularly sampled data, the irregularity and sparseness usually co-exist in local regions when we use smaller spatial subwindows. We may handle the alias by increasing the range of wavenumber spectrum and use the non-aliased low temporal frequency to guide the de-aliasing in the high temporal frequency (Schonewille et al.2009). PKJKP is a unique seismic phase. It travels along the major arc, and thus it stays in the negative wavenumber region of the f–k spectrum. Its slowness is around 1.4 s per degree in the epicentral distance range of our interest. That is quite small in comparison with those of interfering signals. Therefore, even if there are many interfering signals appearing around the PKJKP arrival time window in the t–x domain, PKJKP will be well separated from them in the f–k domain. The arrival along the red dashed line in Fig. 9(b) is not PKJKP itself. As we show in Figs 11(b) and 12(b), besides of the PKJKP phase from the major arc, there is a stronger phase arriving at the same time from the minor arc. That is also the reason why we cannot see consistent signals along the red dashed line. Given that there is always background noise around the PKJKP arrival time, stacking is usually required for the observation of PKJKP in field data (Cao et al.2005; Cao & Romanowicz 2009). 8 CONCLUSIONS We introduce a new data regularization method, NFALFT, based on the non-equispaced NFFT and the ALFT. Synthetic tests from 1-D to 5-D indicate that the NFALFT algorithm is efficient for irregularly sampled seismic data regularization. Both NFALFT and ALFT05 are iteration-based methods and they remove the leaked energy directly from the frequency spectrum. Because the frequency sampling rate for NFALFT is much denser, it is more accurate and faster than ALFT05. The NFALFT algorithm is also faster than the inversion-based ALFT10 algorithm. For the seismic array observation of the Earth's inner core phase PKJKP, we suggest the wavefield separation method to remove all the interfering phase from the minor arc. In practice, all seismic array or regional networks are irregularly distributed with respect to the epicentral distance. Therefore, the NFALFT algorithm-based wavefield separation method provides a new potential way to detect the elusive inner core shear wave phase. ACKNOWLEDGEMENTS The author Aimin Cao deeply appreciates the support of Hamilton Endowment at Southern Methodist University. We thank the editor and two anonymous reviewers for their comments, which help us improve this paper. The source codes are freely available for academic studies by email contact (acao@smu.edu or aimin.f.cao@gmail.com). REFERENCES Abma R., Claerbout J.F., 1995. Lateral prediction for noise attenuation by t-x and f-x techniques, Geophysics  60( 6), 1887– 1896. https://doi.org/10.1190/1.1443920 Google Scholar CrossRef Search ADS   Abma R., Kelley C., Kaldy J., 2007. Sources and treatments of migration-introduced artifacts and noise, in 77th Annual International Meeting, SEG, Expanded Abstracts , pp. 2349– 2353. Bullen K.E., 1951. Theoretical amplitudes of the seismic phase PKJKP, Astron. Soc., Geophys. Suppl.  6 163– 167. https://doi.org/10.1111/j.1365-246X.1951.tb06275.x Google Scholar CrossRef Search ADS   Butzer P.L., 1983. A survey of the Whittaker-Shannon sampling theorem and some of its extensions, J. Math. Res. Exposit.  3 185– 212. Cao A., Romanowicz B., Takeucki N., 2005. An observation of PKJKP: inferences on inner core shear properties, Science  308( 5727), 1453– 1455. https://doi.org/10.1126/science.1109134 Google Scholar CrossRef Search ADS PubMed  Cao A., Romanowicz B., 2009. Constraints on shear wave attenuation in the Earth's inner core from an observation of PKJKP, Geophys. Res. Lett.  36( 9), L09301, doi:10.1029/2009GL038342. https://doi.org/10.1029/2009GL038342 Google Scholar CrossRef Search ADS   Capon J., 1969. High-resolution frequency-wavenumber spectrum analysis, Proc. IEEE  57( 8), 1408– 1418. https://doi.org/10.1109/PROC.1969.7278 Google Scholar CrossRef Search ADS   Cooley J.W., Tukey J.W., 1965. An algorithm for the machine calculation of complex Fourier series, Math. Comp.  19( 90), 297– 297. https://doi.org/10.1090/S0025-5718-1965-0178586-1 Google Scholar CrossRef Search ADS   Dziewonski A., Anderson D., 1981. Preliminary reference earth model, Phys. Earth planet. Inter.  25( 4), 297– 356. https://doi.org/10.1016/0031-9201(81)90046-7 Google Scholar CrossRef Search ADS   Gray D.F., Desikachary K., 1973. A new approach to periodogram analyses, Astrophys. J.  181 523– 530. https://doi.org/10.1086/152068 Google Scholar CrossRef Search ADS   Gulunay N., 2003. Seismic trace interpolation in the Fourier transform domain, Geophysics  68( 1), 355– 369. https://doi.org/10.1190/1.1543221 Google Scholar CrossRef Search ADS   Hindriks K., Duijndam A., 2000. Reconstruction of 3-D seismic signals irregularly sampled along two spatial coordinates, Geophysics  65( 1), 253– 263. https://doi.org/10.1190/1.1444716 Google Scholar CrossRef Search ADS   Jin S., 2010. 5D seismic data regularization by a damped least-norm Fourier inversion, Geophysics  75( 6), WB103– WB111. https://doi.org/10.1190/1.3505002 Google Scholar CrossRef Search ADS   Kawai K., Takeuchi N., Geller R.J., 2006. Complete synthetic seismograms up to 2 Hz for transversely isotropic spherically symmetric media, Geophys. J. Int.  164( 2), 411– 424. https://doi.org/10.1111/j.1365-246X.2005.02829.x Google Scholar CrossRef Search ADS   Kennett B.L.N., Engdahl E.R., 1991. Traveltimes for global earthquake location and phase identification, Geophys. J. Int.  105( 2), 429– 465. https://doi.org/10.1111/j.1365-246X.1991.tb06724.x Google Scholar CrossRef Search ADS   Kennett B.L.N., 2000. Research note Stacking three-component seismograms, Geophys. J. Int.  141( 1), 263– 269. https://doi.org/10.1046/j.1365-246X.2000.00048.x Google Scholar CrossRef Search ADS   Koper K.D., Dombrovskaya M., 2005. Seismic properties of the inner core boundary from PKiKP/P amplitude ratios, Earth planet. Sci. Lett.  237( 3–4), 680– 694. https://doi.org/10.1016/j.epsl.2005.07.013 Google Scholar CrossRef Search ADS   Krüger F., Weber M., Scherbaum F., Schlittenhardt J., 1993. Double beam analysis of anomalies in the core-mantle boundary region, Geophys. Res. Lett.  20( 14), 1475– 1478. https://doi.org/10.1029/93GL01311 Google Scholar CrossRef Search ADS   Kværna T., Doornbos D.J., 1986. An integrated approach to slowness analysis with arrays and three-component stations, NORSAR Semiannu. Techn. Summ.  1 2– 85. Li X.D., Romanowicz B., 1996. Global mantle shear velocity model developed using nonlinear asymptotic coupling theory, J. geophys. Res.  101( B10), 22245– 22272. https://doi.org/10.1029/96JB01306 Google Scholar CrossRef Search ADS   Muirhead K.J., Datt R., 1976. The N-th root process applied to seismic array data, Geophys. J. Int.  47( 1), 197– 210. https://doi.org/10.1111/j.1365-246X.1976.tb01269.x Google Scholar CrossRef Search ADS   Potts D., Steidl G., Tasche M., 2001. Fast Fourier transforms for nonequispaced data: a tutorial, in Modern Sampling Theory: Mathematics and Applications , pp. 247– 270, eds, Benedetto J.J., Ferreira P.J., Birkhauser, Boston. Google Scholar CrossRef Search ADS   Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P., 1993. Numerical Recipes: The art of Scientific Computing  Cambridge University Press. Ristow D., Ruhl T., 1994. Fourier finite–difference migration, Geophysics  59( 12), 1882– 1893. https://doi.org/10.1190/1.1443575 Google Scholar CrossRef Search ADS   Schimmel M., Paulssen H., 1997. Noise reduction and detection of weak, coherent signals through phase-weighted stacks, Geophys. J. Int.  130( 2), 497– 505. https://doi.org/10.1111/j.1365-246X.1997.tb05664.x Google Scholar CrossRef Search ADS   Schonewille M., Kladtke A., Vigner A., 2009. Anti-alias anti-leakage Fourier transform, in 79th Annual International Meeting, SEG, Expanded Abstract , pp. 3249– 3253. Shannon C.E., 1949. Communication in the presence of noise, Proc. IRE  37( 1), 10– 21. https://doi.org/10.1109/JRPROC.1949.232969 Google Scholar CrossRef Search ADS   Shearer P.M., Rychert C.A., Liu Q., 2011. On the visibility of the inner-core shear wave phase PKJKP at long periods, Geophys. J. Int.  185( 3), 1379– 1383. https://doi.org/10.1111/j.1365-246X.2011.05011.x Google Scholar CrossRef Search ADS   Spitz S., 1991. Seismic trace interpolation in the F-X domain, Geophysics  56( 6), 785– 794. https://doi.org/10.1190/1.1443096 Google Scholar CrossRef Search ADS   Takeuchi N., Geller R.J., Cummis P.R., 1996. Highly accurate P-SV complete synthetic seismograms using modified DSM operators, Geophys. Res. Lett.  23( 10), 1175– 1178. https://doi.org/10.1029/96GL00973 Google Scholar CrossRef Search ADS   Trad D., 2009. Five-dimensional interpolation: recovering from acquisition constraints, Geophysics  74( 6), V123– V132. https://doi.org/10.1190/1.3245216 Google Scholar CrossRef Search ADS   Unser M., 2000. Sampling-50 years after Shannon, Proc. IEEE  88( 4), 569– 587. https://doi.org/10.1109/5.843002 Google Scholar CrossRef Search ADS   Wehlau W., Leung K.C., 1964. The multiple periodicity of Delta Delphini, Astrophys. J.  139 843– 863. https://doi.org/10.1086/147820 Google Scholar CrossRef Search ADS   Xu S., Pham D., 2004. Seismic data regularization with anti-leakage Fourier transform, in 66th Annual International Meeting, EAGE, Extended Abstracts , D032. Xu S., Zhang Y., Pham D., Lambare G., 2005. Antileakage Fourier transform for seismic data regularization, Geophysics  70( 4), V87– V95. https://doi.org/10.1190/1.1993713 Google Scholar CrossRef Search ADS   Xu S., Zhang Y., Lambare G., 2010. Antileakage Fourier transform for seismic data regularization in higher dimensions, Geophysics  75( 6), WB113– WB120. https://doi.org/10.1190/1.3507248 Google Scholar CrossRef Search ADS   Zwartjes P.M., Sacchi M.D., 2007. Fourier reconstruction of nonuniformly sampled, aliased seismic data, Geophysics  72( 1), V21– V32. https://doi.org/10.1190/1.2399442 Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# High-resolution seismic data regularization and wavefield separation

, Volume 213 (1) – Apr 1, 2018
11 pages

/lp/ou_press/high-resolution-seismic-data-regularization-and-wavefield-separation-46fbDrNUW7
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy009
Publisher site
See Article on Publisher Site

### Abstract

Summary We present a new algorithm, non-equispaced fast antileakage Fourier transform (NFALFT), for irregularly sampled seismic data regularization. Synthetic tests from 1-D to 5-D show that the algorithm may efficiently remove leaked energy in the frequency wavenumber domain, and its corresponding regularization process is accurate and fast. Taking advantage of the NFALFT algorithm, we suggest a new method (wavefield separation) for the detection of the Earth's inner core shear wave with irregularly distributed seismic arrays or networks. All interfering seismic phases that propagate along the minor arc are removed from the time window around the PKJKP arrival. The NFALFT algorithm is developed for seismic data, but may also be used for other irregularly sampled temporal or spatial data processing. Fourier analysis, Computational seismology, Wave propagation 1 INTRODUCTION The sampling theorem was well developed nearly 70 yr ago (Shannon 1949), and has played a crucial role in signal processing and communication. If a function is frequency bandlimited or pre-filtered with the maximum angular frequency ωmax, it is completely determined by regularly sampled points spaced $$\Delta t = \frac{\pi }{{{\omega _{{\rm{max}}}}}}$$ seconds apart. For regularly sampled signals, Cooley & Tukey (1965) developed the revolutionary fast Fourier transform (FFT) algorithm to accurately and quickly conduct data analysis in the temporal frequency and spatial wavenumber domains. The study of sampling reached a very mature state and became mathematically oriented (e.g. Butzer 1983). In the modern Hilbert-space formulation, Shannon's sampling has been re-interpreted as an orthogonal projection onto the subspace of bandlimited functions with the sampling basis functions (e.g. sinc-function) orthogonal with each other (Unser 2000). However, for irregular non-equispaced sampling, Shannon's sampling theorem is not appropriate and the basis function orthogonality disappears. The fundamental requirement for the FFT is no longer satisfied, and thus the frequency spectrum has to be calculated directly from the original direct Fourier transform (DFT). The slowness of the DFT is not the fatal issue. It is energy leakage from one frequency into others that makes the calculated spectrum inaccurate without further manipulation (e.g. Wehlau & Leung 1964; Press et al.1993). The technique of data regularization is designed to rebuild the correct frequency spectrum and then transform the irregularly sampled data into regularly sampled data. Temporal spectral energy leakage due to the irregular time interval was documented early in astrophysical observations (Wehlau & Leung 1964; Gray & Desikachary 1973). In order to extract very weak frequency components, the technique of pre-whitening was adopted to rebuild the original spectrum in the time or frequency domain. It selects the largest amplitude frequency in the spectrum and subtracts the corresponding sinusoid from the data, then a new spectrum computation is followed free of this frequency and its leaked energy. The new spectrum is searched for the largest amplitude frequency again, and the process is repeated until the noise level is reached. In exploration seismology, due to natural, cultural and economic reasons, the effects of irregular and inadequate spatial sampling of seismic data are very important problem in seismic data processing. Both marine and land data, especially wide azimuth in five dimensions (e.g. shot xs, shot ys, reciever xr, receiver yr; time t), are usually quite irregular and sparse in at least two of the four spatial dimensions due to acquisition and processing costs (Trad 2009). However, seismic migration imaging theories require regular sampling in order to make the correct constructive and destructive interference from data (Abma et al.2007). Fourier finite difference-based wave equation migration explicitly requires regular sampling along each spatial dimension (Ristow & Ruhl 1997). In earthquake seismology, in order to improve the detection of nuclear tests worldwide, seismic arrays have been built, since 1950s to suppress noise and enhance the signal-to-noise ratio. Array apertures range from less than 3 km (e.g. NORES) to more than 100 km (e.g. NORSAR). Most of the geometries of seismic arrays are in certain patterns (e.g. circle, cross, rectangular, etc.), often with irregular spatial sampling. In the case of regional seismic networks, apertures can be larger than hundreds of kilometres (e.g. GRSN in Germany and F-net in Japan), and station distributions are intrinsically irregular. The energy leakage caused by spatial irregularity may affect the measurement of signal slowness and backazimuth in the frequency wavenumber domain with the seismic array/network data (Capon 1969; Kvarna & Doornbos 1986). In particular, if the signal is weaker than the dominant interfering phases, the leaked energy by the stronger phases may bury the weaker signal in the frequency–wavenumber domain (Wehlau & Leung 1964), even if there is no random background noise. Irregular and sparse spatial sampling causes not only energy leakage but also aliasing in the frequency domain. A variety of methods have been suggested for interpolation of spatial aliased data. For instance, Spitz (1991) utilizes the predictability of linear events to obtain an interpolation operator based on the spatial prediction filter in the f–x domain, Abma & Claebout (1995) attenuate random noise with a prediction filter in t–x and f–x domains, and Gulunay (2003) uses non-aliased lower temporal frequencies to de-alias higher frequencies in the f–k domain. However, these methods cannot handle irregular sampling. Regularization of irregularly and sparsely sampled seismic data is primarily conducted in the Fourier domain. The principle is to obtain the correct Fourier coefficients from the irregular data, and then the data regularization may be readily done by an inverse Fourier transform. Fourier reconstruction (e.g. Hindriks & Duijndam 2000; Zwartjes & Sachi 2007) and antileakage Fourier Transform (ALFT, Xu et al.2005, 2010) are two different categories of the method recently developed. Fourier reconstruction is an inverse process for each temporal frequency slice. Because the number of irregularly and sparsely sampled data is usually less than that of inverted Fourier coefficients globally or locally, this kind of inversion is basically ill-conditioned. Least-squares Fourier reconstruction (Hindriks & Duijndam 2000) treats the inversion as a mix-determined issue. It is difficult to find the optimal inverse parameter to balance the misfitting term and the model norm term, and the inverted model cannot exactly recover the original input data. Fourier reconstruction with sparse inversion (Zwartjes & Sacchi 2007) adopts the hybrid norm inversion to solve the aliasing issue using the non-aliased lower frequencies as guide, but the introduction of two inverse parameters makes the inversion even harder in practice. Damped least-norm Fourier inversion (Jin 2010) treats the inversion as a pure underdetermined inverse problem. No inverse parameter is required, and it may completely recover the original input data in theory, but a damping matrix is required to be as close as possible to the final Fourier coefficient model. The ALFT algorithm was first proposed for seismic data regularization by Xu & Pham (2004). In principle, it is similar to the technique of pre-whitening first suggested in astrophysics (Wehlau & Leung 1964; Gray & Desikachary 1973), but ALFT has been implemented for multidimensional seismic data regularizations. It may be used to process 2-D/3-D common offset and common azimuth data (Xu & Pham 2004), and 5-D seismic data by the cascading two passes of 3-D common shot and common receiver regularizations (Xu et al.2005). The accuracy of ALFT regularization is closely related to the size of spatial subwindows. The larger the subwindow, the denser the sampling in the spatial spectrum, and thus the energy leakage can be better removed. However, the large window not only significantly increases computational cost but also makes de-aliasing very difficult. So far, for both Fourier reconstruction and ALFT methods, the most efficient way to handle aliasing is to adopt small spatial subwindows and thus the events inside are basically regarded as linear (Schonewille et al.2009; Jin 2010). In order to solve the conflict between accuracy and computational efficiency for ALFT regularization, Xu et al. (2010) developed a new version of the ALFT algorithm. They introduced a single model parameter inversion method to directly increase the sampling number in the wavenumber domain without changing the spatial input data subwindow. High-resolution accurate regularization for irregularly sampled data is very important not only in exploration seismology and astrophysics, but also in earthquake seismology. Seismic arrays and regional seismic networks are designed to record wavefields propagating along different directions. If multiple wavefields arrive at nearly the same time, they interfere with each other, and thus complicate their separation in the time–space domain, making it difficult to conduct the slowness and backazimuth analysis for individual wavefields based on seismic records. Even if a single wavefield is recorded without noise, the irregular distribution of stations with respect to epicentral distance is still a problem. The resulting energy leakage in the frequency–wavenumber spectrum may hinder us from correctly constraining the slowness and backazimuth of the signals, especially when arrays are small. In order to detect weak signals from background noise, a number of stacking methods (e.g. Array beamforming, Nth root process, Weighted stack methods and Double beam technique) have been developed for seismic array and network data analysis (Muirhead & Datt 1976; Kruger et al.1993; Schimmel & Paulssen 1997; Kennett 2000). The fundamental assumption for these methods is that in the individual traces the noise around the signal is incoherent and no other coherent interfering energy arrives nearby. The stacking technique may significantly improve the signal-to-noise ratio, and plays an important role in nuclear test monitoring and the analysis of subtle anomalies at the core–mantle boundary (Kruger et al.1993) and the constraint of the density contrast at the inner core boundary (Koper & Dombrovskaya 2005). Among all of the weak signals, detection of the Earth's inner core shear wave phase PKJKP, the direct evidence of the solidity of the inner core, is the most challenging (e.g. Bullen 1951; Shearer et al.2011). Based on the IASP91 velocity model (Kennett & Engdahl 1991) and the high-resolution direct solution method (DSM, Takeuchi et al.1996; Kawei et al.2006), synthetic simulation shows that the PKJKP energy is so weak that it is well below the background noise level. The quasi-sphericity of the Earth makes the observation of PKJKP even harder. When the weak PKJKP travels along the major arc of the Earth's big circle to a seismic array or network, a number of other phases arrive at the same time along the opposite minor arc. These interfering phases are much weaker than the first arrival PKIKP and usually unidentified, but they are coherent and stronger than the inner core shear wave phase PKJKP. Current seismic stacking methods cannot handle this kind of situation. This may be one of the main reasons why the observation of PKJKP is so difficult even with high-quality broad-band seismic arrays or networks. In this paper, we present a new method for high-resolution and high-speed seismic data regularization. We first introduce the non-equispaced FFT (NFFT), ALFT algorithms, and then illustrate 1-D/2-D spatial and 5-D temporal and spatial synthetic tests. In the end, we suggest a possible wavefield separation method for the seismic array observation of the inner core shear wave phase PKJKP. 2 NON-EQUISPACED FAST FOURIER TRANSFORM For irregularly sampled seismic data, we may use the non-equispaced DFT (NDFT) to conduct the regularization, but the computational cost is heavy (if not infeasible), especially in the case of multidimensional data processing. Here, we briefly introduce the NDFT and then focus on the flexible and efficient NFFT (Potts et al.2001). Similar to the FFT, the NFFT is another way to fast calculate the frequency spectrum. Their difference is that the FFT works only for the regularly sampled data, and the NFFT can accurately and efficiently handle irregularly sampled data. The data sampling range of dimension $$d \in \mathbb{N}\$$is given by the torus   \begin{eqnarray} &&{{\mathbb{T}^d}:\ = \Big\{ {\boldsymbol{x}} = {{({x_t})}_{t = 0, \ldots ,d - 1}} \in {\mathbb{R}^d}:}\nonumber\\ &&{\quad\qquad-\, \frac{1}{2} \le {x_t} < \frac{1}{2},\ t = 0, \ldots ,d - 1 \Big\}} \end{eqnarray} (1)the sampling set is given by   $$\mathcal{X}: = \left\{ {{{\boldsymbol{x}}_j} \in {\mathbb{T}^d}:j\ = \ 0, \ldots ,M - 1} \right\}$$ (2)and the frequencies $${\boldsymbol{k}} \in {\mathbb{Z}^d}$$ are in the multi-index set   \begin{eqnarray} &&{{I_N}: = \Big\{ {\boldsymbol{k\ }} = {{\left( {{k_t}} \right)}_{t = 0, \ldots ,d - 1}}\ \in {\mathbb{Z}^d}:}\nonumber\\ &&{\qquad\quad- \frac{{{N_t}}}{2} \le {k_t} < \frac{{{N_t}}}{2},\ t\ = \ 0, \ldots ,d - 1 \Big\}} \end{eqnarray} (3)where $$N\ = {({N_t})_{t = 0, \ldots ,d - 1\! }},\ {N_t} \in 2\mathbb{N}$$ is the even multiband limit. For a finite number of given Fourier coefficients, $${f_{\boldsymbol{k}}} \in \mathbb{C},{\boldsymbol{\ k}} \in {I_N}$$, the forward NDFT may be expressed as   $${f_{\boldsymbol{k}}} = \mathop \sum \nolimits_{j = 0}^{M - 1} w\left( {{{\boldsymbol{x}}_j}} \right){f_j}{e^{ - 2\pi i{\boldsymbol{k}}{{\boldsymbol{x}}_{\boldsymbol{j}}}}}$$ (4)where w(xj) is the integral weight and $$\mathop \sum \nolimits_{j\ = \ 0}^{M - 1} w\ ( {{{\boldsymbol{x}}_j}} ) = 1.0$$. It takes $$\mathcal{O}( {M\,{\rm{N}}} )$$ arithmetical operations. The core of the ALFT technique is to clean leaked energy in the frequency domain, thus the forward NFFT, which transforms the irregularly sampled data from spatial domain to wavenumber domain, plays the key role in accuracy and speed of the ALFT regularization. Here, we summarize the forward NFFT algorithm as follows: $$Input{:}\ d,\ M \in \mathbb{N},\ N \in 2{N^d}$$  \begin{eqnarray} &&{{{\boldsymbol{x}}_j} \in {\left[ { -\frac{1}{2},\frac{1}{2}} \right]^d},\ j\ = \ 0, \ldots ,M - 1,\quad and\quad {f_j} \in \mathbb{C},}\nonumber\\ &&{\quad j\ = \ 0, \ldots ,M - 1} \end{eqnarray} (5) 1: For l ∈ In compute  $${g_{\boldsymbol{l}}}: = \sum\limits_{j \in {I_{n,m}}\left( l \right)} {f_j}\psi \left( {{{\boldsymbol{x}}_j} - {{\boldsymbol{n}}^{ - 1}} \odot {\boldsymbol{l}}} \right)$$ (6) 2: For k ∈ In compute d dimensional forward FFT  $${g_{\boldsymbol{k}}}: = \sum\limits_{{\boldsymbol{l}} \in {I_n}} {g_l}{e^{ - 2\pi i{\boldsymbol{k}}\left( {{{\boldsymbol{n}}^{ - 1}} \odot {\boldsymbol{l}}} \right)}}$$ (7) 3: For k ∈ IN compute  $${f_{\boldsymbol{k}}}: = \frac{{{g_{\boldsymbol{k}}}}}{{\left| {{I_n}} \right|{c_k}\left( \varphi \right)}}$$ (8) output:  fk, k ∈ IN $$\mathcal{O}( {| N |log| N | + M} )\ arithmetic\ operations$$ where xj is the normalized input sampling range, ψ is a one periodic window function, ck is the Fourier coefficients of the periodic window function and  In, m(l): ={j = 0, …, M − 1: l − m ≤ n⊙xj ≤ l + m}. 3 ANTILEAKAGE FOURIER TRANSFORM ALGORITHM Two versions of ALFT algorithms were developed for the multidimensional regularization of irregularly sampled seismic data. The first ALFT (Xu et al.2005) is similar to the pre-whitening technique (Wehlau & Leung 1964; Gray & Desikachary 1973) (hereafter it is referred as ALFT05), and the second ALFT is a one-parameter inversion related technique (Xu et al.2010) (hereafter it is referred as ALFT10). Here, we introduce both algorithms for theoretical integrity and convenience in the following discussion in this paper. ALFT05 Set all initial Fourier components to zero. Calculate all Fourier coefficients of original/residual input data with the NDFT (eq. 4). Find the wavenumber kmaxand its Fourier coefficient$${f_{{{\boldsymbol{k}}_{\rm {max}}}}}$$ with the maximum energy, and keep adding the Fourier coefficient$${f_{{{\boldsymbol{k}}_{\rm {max}}}}}$$to the initial Fourier components. Calculate the wavenumber kmaxcomponent in the original/residual input data   \begin{equation*}{f^{{{\boldsymbol{k}}_{\rm {max}}}}}( {{x_j}} ) = {f_{{{\boldsymbol{k}}_{\rm {max}}}}}{e^{2\pi i{\boldsymbol{k}}{{\boldsymbol{x}}_j}}}\end{equation*} Update the input by subtracting the wavenumber kmax component   \begin{equation*}{f^u_j} = {f_j}\ - {f^{{{\boldsymbol{k}}_{\rm {max}}}}}( {{{\boldsymbol{x}}_j}} ),\quad j\ = \ 0, \ldots ,M - 1\end{equation*} Iterate steps (b) to (e) until the leaked energy disappear or is under a threshold. ALFT10 Similar to ALFT05, Fourier coefficients are instead calculated with the inversion by minimizing the objective function   $$E\ = \sum\limits_{j = 1}^{M - 1} {\left\| {{f_j} - {f_{{\boldsymbol{k'}}}}{e^{2\pi i{\boldsymbol{k'}}{{\boldsymbol{x}}_j}}}} \right\|}$$ (9)where k' is the oversampled wavenumber in spatial frequency domain. ALFT05 is an iterative technique. The computation can be extremely heavy for the entire process of regularization with the NDFT. The adoption of the NFFT may significantly reduce the computation time, but the accuracy of the resulting spectrum is still a concern. In other words, ALFT05 cannot effectively remove the leaked energy in the spatial spectrum, especially when small spatial subwindows are selected in order to handle the sampling alias. ALFT10 is designed to improve the spectrum accuracy by significantly oversampling the wavenumber k' N times (e.g. N = 32) in eq. (9). It means that the new sampling rate Δk'is N times denser than in an equivalent DFT case. The adoption of oversampled wavenumber k' may help remove the leaked energy in the spectrum domain. However, what this inversion-oriented ALFT10 algorithm obtains in the gain of accuracy, it loses in speed, because the inversion takes N2 more time than the NDFT in ALFT05 to calculate the individual frequency spectrum. 4 NON-EQUISPACED FAST ANTILEAKAGE FOURIER TRANSFORM ALGORITHM Here, we suggest a new algorithm for the high accuracy and high-speed regularization of irregularly sampled data based on NFFT and ALFT05. Instead of directly oversampling the wavenumber in the frequency domain as ALFT10 does, we introduce the d dimensional input sampling parameter vector $${\boldsymbol{\sigma }} \in \mathbb{N}$$ and change the NFFT input sampling set (eq. 5) as   $${\boldsymbol{x}}_j^{\prime} = \ {{\boldsymbol{x}}_j} \odot \frac{1}{{\boldsymbol{\sigma }}},\quad {\boldsymbol{x}}_j^{\boldsymbol{=}} \in {\left[ { - \frac{1}{2},\frac{1}{2}} \right]^d}\quad j\ = \ 0, \ldots ,M - 1$$ (11) After being divided by σ, the original sampling set xj and the new sampling set $${\boldsymbol{x}}_j^{\prime}$$ are still in the normalized ranges with the same number of sampling input. However, the difference is that $$x_j^{\prime}$$ is now equivalent to a subset of zero-padded new set of size σ · M. The introduction of the input sampling parameter σ increases the frequency sampling rate by (σi − 1) times in ith dimension, and thus the resulting new frequency resolution is $$\mathop \prod \nolimits_{i\ = \ 0}^{d - 1} ({\sigma _i} - 1)$$ times higher. This means that the frequency of the maximum energy in each iteration now has a much greater chance to locate at or very close to a frequency sampling point in the spectrum. Therefore, we may not only remove the leaked energy more efficiently in smaller number of iterations, but also obtain a more accurate spectrum in the end. NFFT is a fast and good approximation to NDFT with a relative error usually less than 0.001 (Potts et al.2001). Therefore, in practice, we always use NFFT rather than NDFT. For the high-dimensional (especially 5-D) regularization, in order to significantly reduce the computational memory and time cost, we suggest a modification of the original NFFT algorithm. Given the fact that the input sampling locations are always the same for different temporal frequency slices and different iterations in a subwindow, we just need to calculate those one-periodic window functions in eq. (6) one time and then save them inside a shared memory for the parallel computation. For convenience in the following tests and discussion, hereafter we refer to our new algorithm as the non-equispaced fast ALFT (NFALFT). 5 SYNTHETIC TESTS In order to quantitatively examine the performance of the NFALFT algorithm, three sets of analytical function are adopted to calculate the regularly and irregularly sampled synthetic data, respectively. Both ALFT05 and NFALFT regularizations are conducted for the irregularly sampled synthetic data. Then, we compare the regularized data with the regularly calculated data from these analytical functions. We start with a 1-D sinusoidal function y = 5.0*cos (0.5*x) + 2.5*cos(x). If the range of the x window is infinite, there are supposed to be two δ-functions in the spatial frequency (i.e. wavenumber) spectrum. One is at frequency 0.5 rad m−1 with amplitude 5.0, and the other is at frequency 1.0 rad m−1 with amplitude 2.5. If the range of the x window is limited and the sampling is regular, there are supposed to be two limit-sized pulses in the frequency spectrum with smaller amplitudes. If the range of the x window is limited and the sampling is irregular (randomly sampled 64 points in the x range from 0 to 31) (Fig. 1), the frequency spectrum is much more complicated and energy is leaked across the whole frequency range (the green line in Fig. 2a). After applying the antileakage processing with ALFT05 algorithm, most of the leaked energy is removed (the red line in Fig. 2a), but the resulting regularized results (the red line in Fig. 2b) do not match the analytical function (the blue line in Fig. 2b) well. This means that the ALFT05 algorithm cannot handle the leaked energy accurately in the frequency spectrum. Figure 1. View largeDownload slide The analytical function y = 5.0*cos (0.5*x) + 2.5*cos(x) (in red line) and the random sampling (in black cross). Figure 1. View largeDownload slide The analytical function y = 5.0*cos (0.5*x) + 2.5*cos(x) (in red line) and the random sampling (in black cross). Figure 2. View largeDownload slide (a) Frequency amplitude spectrum with NDFT (in green) and ALFT05 processing (in red). (b) Regularization with ALFT05 processing (in red) and the original analytical function (in blue). Figure 2. View largeDownload slide (a) Frequency amplitude spectrum with NDFT (in green) and ALFT05 processing (in red). (b) Regularization with ALFT05 processing (in red) and the original analytical function (in blue). Using the NFALFT algorithm, we take the spatial input sampling parameter σ = 32 to oversample the wavenumber k up to 31 times in the frequency domain. After applying the antileakage processing with NFALFT, just a number of discrete frequency with non-zero energy are left (the red line in Fig. 3a). The largest four amplitudes are 4.78 at k = 1.00 rad m−1, 2.44 at k = 0.50 rad m−1, 0.37 at k = 0.48 rad m−1 and 0.04 at 0.37 rad m−1, respectively. Energy at other frequencies are zeros or very close to zeros. The resulting regularized results by NFALFT are identical to the calculated ones from the analytical function (Fig. 3b). The thresholds for both ALFT05 and NFALFT antileakage processing are set the same (0.001). The NFALFT regularization reaches the threshold in 10 iterations, on the other hand, the ALFT05 algorithm took 1093 iterations to reach the threshold. The NFALFT is more than 2 times faster than the ALFT05, even though the frequency is 32 times oversampled relative to the ALFT05. Figure 3. View largeDownload slide (a) Frequency amplitude spectrum with NDFT (in green) and NFALFT processing (in red). (b) Regularization after NFALFT processing (in red) and the original analytical function (in blue). The two curves overlap with each other and thus only the red one is visible. Figure 3. View largeDownload slide (a) Frequency amplitude spectrum with NDFT (in green) and NFALFT processing (in red). (b) Regularization after NFALFT processing (in red) and the original analytical function (in blue). The two curves overlap with each other and thus only the red one is visible. The energy leakage in the frequency domain caused by the spatial sampling irregularity and the principal of antileakage algorithms has been demonstrated in the above 1-D sinusoidal synthetic test. In the next two tests, we will focus on the comparison of regularization with ALFT05 and AFALFT algorithms. In the second synthetic test, we use the 2-D sinusoidal analytical function z = sin (x) *cos(y). Its contour map and corresponding 102 × 102 regular sampling grids are presented in Fig. 4. We randomly remove 50 per cent of the grid to create a sparser and irregular sampled grid (Fig. 5a). Figure 4. View largeDownload slide (a) Contour map of the 2-D analytical function z = sin (x) *cos (y). (b) Regular sampling with 102 × 102 gridpoints for the analytical function. Figure 4. View largeDownload slide (a) Contour map of the 2-D analytical function z = sin (x) *cos (y). (b) Regular sampling with 102 × 102 gridpoints for the analytical function. Figure 5. View largeDownload slide (a) Irregular sampling by randomly removing 50 per cent of the original sampling points. (b) Contour map of the 2-D regularization result with the ALFT05 algorithm. Figure 5. View largeDownload slide (a) Irregular sampling by randomly removing 50 per cent of the original sampling points. (b) Contour map of the 2-D regularization result with the ALFT05 algorithm. For the NFALFT algorithm, we take the spatial input sampling parameter σ  = (4, 4) to oversample the wavenumber k up to 15 times in the frequency domain. After applying the ALFT05 algorithm, the contour map plotted with the regularized data (Fig. 5b) does not match the original one from the analytical data (Fig. 4a). The distortion at the boundaries is significant. The regularization results look better in the middle, but small kinks are distributed along each of the contour lines. After applying the NFALFT algorithm, we recover nearly the same contour map as with the regularized data (Fig. 6). There are slight disconnections at contour lines of value 0. The thresholds for both ALFT05 and NFALFT antileakage processing are set as the same (0.001). The NFALFT regularization reaches the threshold in 77 iterations, on the other hand, the ALFT05 algorithm took 2733 iterations to reach the threshold. NFALFT is more than 14 times faster than ALFT05, even though the frequency is 16 times oversampled relative to ALFT05. Figure 6. View largeDownload slide Contour map of the 2-D regularization result with the NFALFT algorithm. Figure 6. View largeDownload slide Contour map of the 2-D regularization result with the NFALFT algorithm. Our third synthetic test is conducted with three 5-D hyperplane analytical functions.   $$A\!:3{x_1} + 2{x_2} - 2{x_3} - 2{x_4} - 25000t + 11600\ = \ 0$$ (12)  $$B\!:3{x_1} + 2{x_2} + 8{x_3} + 8{x_4} - 25000t + 2100\ = \ 0\$$ (13)  $$C\!:3{x_1} + 2{x_2} - 4{x_3} - 4{x_4} - 25000t + 14\ 000\ = \ 0$$ (14)where x1, x2, x3, and x4 are four spatial coordinates, which may be shot x, shot y, receiver x and receiver y, or offset x, offset y, cdp x and cdp y in exploration seismology and t is the temporal coordinate. For a regular 16 × 16 × 16 × 16  sampling in all four spatial dimensions, we show the sampling grid in the x3x4 plane in Fig. 7(a) and the corresponding profile of the 5-D analytical hyperplane in Fig. 8(a). In order to get irregular 16 × 16 × 16 × 16 sampling (Fig. 7b), we shift the original sampling points randomly in each spatial dimension. For the NFALFT algorithm, we take the spatial input sampling parameter σ = (2, 2, 2, 2) to oversample the wavenumber k up to 15 times in the frequency domain. After applying the NFALFT algorithm, the regularized hyperplanes (Fig. 8b) are basically the same as the ones calculated from the analytical functions with the regular sampling grids (Fig. 8c). The computation of the ALFT05 algorithm does not converge to threshold (0.001), and so here we do not provide a quantitative comparison of the computational cost between ALFT05 and NFALFT for this 5-D test. Figure 7. View largeDownload slide (a) The regular sampling with 16 × 16 gridpoints in x3x4 plane for the 5-D hypeplane function in eqs (12)–(14). (b) The random sampling in x3x4 plane in the 5-D hypeplane function. Figure 7. View largeDownload slide (a) The regular sampling with 16 × 16 gridpoints in x3x4 plane for the 5-D hypeplane function in eqs (12)–(14). (b) The random sampling in x3x4 plane in the 5-D hypeplane function. Figure 8. View largeDownload slide Comparison of original data and regularized data in x4 − t plane. (a) Original data calculated with eqs (12)–(14). (b) Regularized data by means of the NFALFT algorithm. (c) Amplitude difference between the two panels of (a) and (b). Figure 8. View largeDownload slide Comparison of original data and regularized data in x4 − t plane. (a) Original data calculated with eqs (12)–(14). (b) Regularized data by means of the NFALFT algorithm. (c) Amplitude difference between the two panels of (a) and (b). 6 SEISMIC WAVEFIELD SEPARATION Strong earthquakes (e.g. Mw > 6.5) usually trigger seismic wave propagation around the Earth along the minor and the major arcs, respectively. This is good for the global/regional waveform tomography (e.g. Li & Romanowicz 1996) with denser ray path coverage, but it is a challenge for seismic array observation. The inner core shear wave PKJKP is a unique phase to verify the solidity of the Earth's centre. However, due to the high attenuation in the inner core and the poor P-to-S and S-to-P conversions at the inner core boundary, PKJKP is so weak that it only barely stands out of the background noise (e.g. Shearer et al.2011). PKJKP travels along the major arc. Within its theoretical arrival time window, there are always a number of unidentified weak phases that arrive along the minor arc at nearly the same time. These minor arc interference phases are much weaker than the first arrival PKIKP, but they are still strong in comparison with the potential PKJKP, making the observation of PKJKP very difficult with seismic array/network stacking methods in the temporal–spatial domain. Given the fact that PKJKP and the interference phases travel along opposite directions in a great circle of the Earth, one should be able to separate seismic wavefields along the major arc from those along the minor arc in the temporal frequency and spatial wavenumber f–k domain. PKJKP appears in the negative wavenumber spectrum and those minor arc interference phases are in the positive wavenumber spectrum. However, in practice, due to the spatial irregularity of station distribution in arrays/networks, the energy leakage of strong minor arc wavefields into weak major arc wavefields in the f–k spectrum makes the wavefield separation difficult without further manipulation. In the following PKJKP synthetic analyses, we first show the energy leakage in the f–k domain caused by the irregular spatial distribution of a seismic array, and then demonstrate the efficiency of the NFALFT algorithm in the removal of leaked energy and wavefield separation. In order to obtain an irregularly distributed seismic array, in the epicentral distance range from 138° to 140°, we first set up a linear array of 40 stations with a regular station interval 5.55 km along the backazimuthal direction of 90°, and then randomly shift each station from the second to the 39th within the range of ± 2.6 km along the array direction. Synthetic seismograms (Fig. 9a) are computed with the high-resolution DSM (Kawei et al.2006) for an event Mw = 7.0 with respect to the IASP91 velocity model (Kennett & Engdahl 1991). The high corner frequency is 0.5 Hz, and the inner core shear quality factor is taken as 350, which is much larger than 85 first suggested in the PREM model (Dziewonski & Anderson 1981). Around the theoretical arrival time of PKJKP, we cannot identify any coherent seismic phase dipping from right to left, even after the data regularization (Fig. 9b). PKJKP is buried in the wavefields from the minor arc. Figure 9. View largeDownload slide Synthetic seismograms computed with the direct solution method (DSM) for the inner core shear wave phase PKJKP in the epicentral distance range of 138° − 140°. The theoretical arrival times of PKJKP are ∼1710 s (red dashed lines). (a) Seismograms from irregularly distributed 40 stations. (b) Regularized seismograms of 40 stations with regular interval of 5.55 km by means of the NFALFT algorithm. Figure 9. View largeDownload slide Synthetic seismograms computed with the direct solution method (DSM) for the inner core shear wave phase PKJKP in the epicentral distance range of 138° − 140°. The theoretical arrival times of PKJKP are ∼1710 s (red dashed lines). (a) Seismograms from irregularly distributed 40 stations. (b) Regularized seismograms of 40 stations with regular interval of 5.55 km by means of the NFALFT algorithm. In the f–k domain, the energy propagating along the minor arc is dominant in the positive wavenumber region, meanwhile, the energy leakage from the positive wavenumber into the negative wavenumber is also large due to the irregular spatial sampling (Fig. 10a). The leaked energy is so strong that it not only severely distorts the PKJKP waveforms but also make the PKJKP amplitude much weaker (Fig. 11a). After applying the NFALFT algorithm, we remove the leaked energy in the f–k domain (Fig. 10b). The wavefields along the major arc are well separated from the wavefields along the minor arc in both f–k domain (Fig. 10b) and x–t domain (Figs 11b and 12b). After this wavefield separation, we clearly see the PKJKP phase in individual traces in the major arc wavefield without stacking (Fig. 11b). There are also a number of other major arc phases. Two of them around relative time 1730 s, we identify as pPKJKP + PKPPKPdf. There are more minor arc phases in the time window (Fig. 12b). It is difficult to identify them unambiguously, but one of them arrives at nearly the same time as PKJKP does but with larger amplitudes. Figure 10. View largeDownload slide The temporal frequency versus spatial wavenumber f–k spectra. (a) The spectrum calculated for the irregularly distributed seismograms (Fig. 9a) with the NDFT method. The energy leakage is remarkable at each temporal frequency slice. (b) Removal of the leaked energy by the NFALFT algorithm. All seismic energy propagating along the minor arc is on the right, and all energy propagating along the major arc is on the left. Figure 10. View largeDownload slide The temporal frequency versus spatial wavenumber f–k spectra. (a) The spectrum calculated for the irregularly distributed seismograms (Fig. 9a) with the NDFT method. The energy leakage is remarkable at each temporal frequency slice. (b) Removal of the leaked energy by the NFALFT algorithm. All seismic energy propagating along the minor arc is on the right, and all energy propagating along the major arc is on the left. Figure 11. View largeDownload slide Seismic phases propagating along the major arc after the wavefield separation. Dashed red lines are the theoretical PKJKP arrival time. (a) Waveform distortion caused by the leaked energy from the minor arc wave energy. (b) More coherent and stronger PKJKP waveforms after the removal of the leaked energy by the NFALFT algorithm. Figure 11. View largeDownload slide Seismic phases propagating along the major arc after the wavefield separation. Dashed red lines are the theoretical PKJKP arrival time. (a) Waveform distortion caused by the leaked energy from the minor arc wave energy. (b) More coherent and stronger PKJKP waveforms after the removal of the leaked energy by the NFALFT algorithm. Figure 12. View largeDownload slide Seismic phases propagating along the minor arc after the wavefield separation. (a) The leaked energy also affects the strong minor arc wavefields. (b) The minor arc wavefields after the NFALFT regularization. One phase arrives at nearly the same time as PKJKP does (red dashed lines). Figure 12. View largeDownload slide Seismic phases propagating along the minor arc after the wavefield separation. (a) The leaked energy also affects the strong minor arc wavefields. (b) The minor arc wavefields after the NFALFT regularization. One phase arrives at nearly the same time as PKJKP does (red dashed lines). 7 DISCUSSION NFALFT and ALFT05 are iteration-based antileakage algorithms. Their computational costs mainly depend on the number of iterations for them to reach the threshold (e.g. 0.001). The efficiency of the leaked energy removal in each iteration is closely related to the frequency sampling rate in each dimension. The denser the sampling rate, the more likely the leaked energy of the dominant frequency may be removed in a single iteration. Otherwise, the iteration process can take a long time before reaching the threshold, and meanwhile the resulting accuracy of the frequency spectrum is poor. ALFT10 is an inversion-based antileakage algorithm to remove the leaked energy in the frequency spectrum. It directly increases the sampling rate in each dimension of the frequency spectrum and conducts a one-parameter least-squares inversion for each frequency in the spectrum. Thus, the computational cost of ALFT10 is $${{\boldsymbol \sigma }^{\boldsymbol 2}}$$ greater than that of NDFT. On the other hand, NFALFT takes the NFFT to compute the frequency spectrum, and thus the arithmetical operation contrast between ALFT10 and NFALFT is approximately $$\mathcal{O}( {{\sigma ^2}M \cdot {\rm{N}}} )$$ to $$\mathcal{O}( {| {\sigma \cdot N} |{\rm{log}}| {\sigma \cdot N} | + M} )$$. If the oversampling rates are the same for NFALFT and ALFT10, their accuracies are the same in theory. In practice, however, NFALFT may be more accurate because it is based on the NFFT, which provides a more accurate frequency spectrum than the inversion method does. The iteration number is usually less than 10 per cent of the total number of frequencies in the spectrum for NFALFT. Especially, for the linear events in high-dimensional data, the corresponding frequency spectrum is very sparse, and so the iteration number of NFALFT is even smaller. The iteration threshold is the only parameter we need to specify for the NFALFT algorithm, and usually we take it as a constant (e.g. 0.001) for all regularization. In comparison with the Fourier reconstruction inversion methods (e.g. Zwartjes & Sacchi 2007; Jin 2010), there are not trial-and-error inverse parameters or pre-conditional starting models required. For the high-dimensional irregularly sampled data, the irregularity and sparseness usually co-exist in local regions when we use smaller spatial subwindows. We may handle the alias by increasing the range of wavenumber spectrum and use the non-aliased low temporal frequency to guide the de-aliasing in the high temporal frequency (Schonewille et al.2009). PKJKP is a unique seismic phase. It travels along the major arc, and thus it stays in the negative wavenumber region of the f–k spectrum. Its slowness is around 1.4 s per degree in the epicentral distance range of our interest. That is quite small in comparison with those of interfering signals. Therefore, even if there are many interfering signals appearing around the PKJKP arrival time window in the t–x domain, PKJKP will be well separated from them in the f–k domain. The arrival along the red dashed line in Fig. 9(b) is not PKJKP itself. As we show in Figs 11(b) and 12(b), besides of the PKJKP phase from the major arc, there is a stronger phase arriving at the same time from the minor arc. That is also the reason why we cannot see consistent signals along the red dashed line. Given that there is always background noise around the PKJKP arrival time, stacking is usually required for the observation of PKJKP in field data (Cao et al.2005; Cao & Romanowicz 2009). 8 CONCLUSIONS We introduce a new data regularization method, NFALFT, based on the non-equispaced NFFT and the ALFT. Synthetic tests from 1-D to 5-D indicate that the NFALFT algorithm is efficient for irregularly sampled seismic data regularization. Both NFALFT and ALFT05 are iteration-based methods and they remove the leaked energy directly from the frequency spectrum. Because the frequency sampling rate for NFALFT is much denser, it is more accurate and faster than ALFT05. The NFALFT algorithm is also faster than the inversion-based ALFT10 algorithm. For the seismic array observation of the Earth's inner core phase PKJKP, we suggest the wavefield separation method to remove all the interfering phase from the minor arc. In practice, all seismic array or regional networks are irregularly distributed with respect to the epicentral distance. Therefore, the NFALFT algorithm-based wavefield separation method provides a new potential way to detect the elusive inner core shear wave phase. ACKNOWLEDGEMENTS The author Aimin Cao deeply appreciates the support of Hamilton Endowment at Southern Methodist University. We thank the editor and two anonymous reviewers for their comments, which help us improve this paper. The source codes are freely available for academic studies by email contact (acao@smu.edu or aimin.f.cao@gmail.com). REFERENCES Abma R., Claerbout J.F., 1995. Lateral prediction for noise attenuation by t-x and f-x techniques, Geophysics  60( 6), 1887– 1896. https://doi.org/10.1190/1.1443920 Google Scholar CrossRef Search ADS   Abma R., Kelley C., Kaldy J., 2007. Sources and treatments of migration-introduced artifacts and noise, in 77th Annual International Meeting, SEG, Expanded Abstracts , pp. 2349– 2353. Bullen K.E., 1951. Theoretical amplitudes of the seismic phase PKJKP, Astron. Soc., Geophys. Suppl.  6 163– 167. https://doi.org/10.1111/j.1365-246X.1951.tb06275.x Google Scholar CrossRef Search ADS   Butzer P.L., 1983. A survey of the Whittaker-Shannon sampling theorem and some of its extensions, J. Math. Res. Exposit.  3 185– 212. Cao A., Romanowicz B., Takeucki N., 2005. An observation of PKJKP: inferences on inner core shear properties, Science  308( 5727), 1453– 1455. https://doi.org/10.1126/science.1109134 Google Scholar CrossRef Search ADS PubMed  Cao A., Romanowicz B., 2009. Constraints on shear wave attenuation in the Earth's inner core from an observation of PKJKP, Geophys. Res. Lett.  36( 9), L09301, doi:10.1029/2009GL038342. https://doi.org/10.1029/2009GL038342 Google Scholar CrossRef Search ADS   Capon J., 1969. High-resolution frequency-wavenumber spectrum analysis, Proc. IEEE  57( 8), 1408– 1418. https://doi.org/10.1109/PROC.1969.7278 Google Scholar CrossRef Search ADS   Cooley J.W., Tukey J.W., 1965. An algorithm for the machine calculation of complex Fourier series, Math. Comp.  19( 90), 297– 297. https://doi.org/10.1090/S0025-5718-1965-0178586-1 Google Scholar CrossRef Search ADS   Dziewonski A., Anderson D., 1981. Preliminary reference earth model, Phys. Earth planet. Inter.  25( 4), 297– 356. https://doi.org/10.1016/0031-9201(81)90046-7 Google Scholar CrossRef Search ADS   Gray D.F., Desikachary K., 1973. A new approach to periodogram analyses, Astrophys. J.  181 523– 530. https://doi.org/10.1086/152068 Google Scholar CrossRef Search ADS   Gulunay N., 2003. Seismic trace interpolation in the Fourier transform domain, Geophysics  68( 1), 355– 369. https://doi.org/10.1190/1.1543221 Google Scholar CrossRef Search ADS   Hindriks K., Duijndam A., 2000. Reconstruction of 3-D seismic signals irregularly sampled along two spatial coordinates, Geophysics  65( 1), 253– 263. https://doi.org/10.1190/1.1444716 Google Scholar CrossRef Search ADS   Jin S., 2010. 5D seismic data regularization by a damped least-norm Fourier inversion, Geophysics  75( 6), WB103– WB111. https://doi.org/10.1190/1.3505002 Google Scholar CrossRef Search ADS   Kawai K., Takeuchi N., Geller R.J., 2006. Complete synthetic seismograms up to 2 Hz for transversely isotropic spherically symmetric media, Geophys. J. Int.  164( 2), 411– 424. https://doi.org/10.1111/j.1365-246X.2005.02829.x Google Scholar CrossRef Search ADS   Kennett B.L.N., Engdahl E.R., 1991. Traveltimes for global earthquake location and phase identification, Geophys. J. Int.  105( 2), 429– 465. https://doi.org/10.1111/j.1365-246X.1991.tb06724.x Google Scholar CrossRef Search ADS   Kennett B.L.N., 2000. Research note Stacking three-component seismograms, Geophys. J. Int.  141( 1), 263– 269. https://doi.org/10.1046/j.1365-246X.2000.00048.x Google Scholar CrossRef Search ADS   Koper K.D., Dombrovskaya M., 2005. Seismic properties of the inner core boundary from PKiKP/P amplitude ratios, Earth planet. Sci. Lett.  237( 3–4), 680– 694. https://doi.org/10.1016/j.epsl.2005.07.013 Google Scholar CrossRef Search ADS   Krüger F., Weber M., Scherbaum F., Schlittenhardt J., 1993. Double beam analysis of anomalies in the core-mantle boundary region, Geophys. Res. Lett.  20( 14), 1475– 1478. https://doi.org/10.1029/93GL01311 Google Scholar CrossRef Search ADS   Kværna T., Doornbos D.J., 1986. An integrated approach to slowness analysis with arrays and three-component stations, NORSAR Semiannu. Techn. Summ.  1 2– 85. Li X.D., Romanowicz B., 1996. Global mantle shear velocity model developed using nonlinear asymptotic coupling theory, J. geophys. Res.  101( B10), 22245– 22272. https://doi.org/10.1029/96JB01306 Google Scholar CrossRef Search ADS   Muirhead K.J., Datt R., 1976. The N-th root process applied to seismic array data, Geophys. J. Int.  47( 1), 197– 210. https://doi.org/10.1111/j.1365-246X.1976.tb01269.x Google Scholar CrossRef Search ADS   Potts D., Steidl G., Tasche M., 2001. Fast Fourier transforms for nonequispaced data: a tutorial, in Modern Sampling Theory: Mathematics and Applications , pp. 247– 270, eds, Benedetto J.J., Ferreira P.J., Birkhauser, Boston. Google Scholar CrossRef Search ADS   Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P., 1993. Numerical Recipes: The art of Scientific Computing  Cambridge University Press. Ristow D., Ruhl T., 1994. Fourier finite–difference migration, Geophysics  59( 12), 1882– 1893. https://doi.org/10.1190/1.1443575 Google Scholar CrossRef Search ADS   Schimmel M., Paulssen H., 1997. Noise reduction and detection of weak, coherent signals through phase-weighted stacks, Geophys. J. Int.  130( 2), 497– 505. https://doi.org/10.1111/j.1365-246X.1997.tb05664.x Google Scholar CrossRef Search ADS   Schonewille M., Kladtke A., Vigner A., 2009. Anti-alias anti-leakage Fourier transform, in 79th Annual International Meeting, SEG, Expanded Abstract , pp. 3249– 3253. Shannon C.E., 1949. Communication in the presence of noise, Proc. IRE  37( 1), 10– 21. https://doi.org/10.1109/JRPROC.1949.232969 Google Scholar CrossRef Search ADS   Shearer P.M., Rychert C.A., Liu Q., 2011. On the visibility of the inner-core shear wave phase PKJKP at long periods, Geophys. J. Int.  185( 3), 1379– 1383. https://doi.org/10.1111/j.1365-246X.2011.05011.x Google Scholar CrossRef Search ADS   Spitz S., 1991. Seismic trace interpolation in the F-X domain, Geophysics  56( 6), 785– 794. https://doi.org/10.1190/1.1443096 Google Scholar CrossRef Search ADS   Takeuchi N., Geller R.J., Cummis P.R., 1996. Highly accurate P-SV complete synthetic seismograms using modified DSM operators, Geophys. Res. Lett.  23( 10), 1175– 1178. https://doi.org/10.1029/96GL00973 Google Scholar CrossRef Search ADS   Trad D., 2009. Five-dimensional interpolation: recovering from acquisition constraints, Geophysics  74( 6), V123– V132. https://doi.org/10.1190/1.3245216 Google Scholar CrossRef Search ADS   Unser M., 2000. Sampling-50 years after Shannon, Proc. IEEE  88( 4), 569– 587. https://doi.org/10.1109/5.843002 Google Scholar CrossRef Search ADS   Wehlau W., Leung K.C., 1964. The multiple periodicity of Delta Delphini, Astrophys. J.  139 843– 863. https://doi.org/10.1086/147820 Google Scholar CrossRef Search ADS   Xu S., Pham D., 2004. Seismic data regularization with anti-leakage Fourier transform, in 66th Annual International Meeting, EAGE, Extended Abstracts , D032. Xu S., Zhang Y., Pham D., Lambare G., 2005. Antileakage Fourier transform for seismic data regularization, Geophysics  70( 4), V87– V95. https://doi.org/10.1190/1.1993713 Google Scholar CrossRef Search ADS   Xu S., Zhang Y., Lambare G., 2010. Antileakage Fourier transform for seismic data regularization in higher dimensions, Geophysics  75( 6), WB113– WB120. https://doi.org/10.1190/1.3507248 Google Scholar CrossRef Search ADS   Zwartjes P.M., Sacchi M.D., 2007. Fourier reconstruction of nonuniformly sampled, aliased seismic data, Geophysics  72( 1), V21– V32. https://doi.org/10.1190/1.2399442 Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations