Jointly filtering and regularizing seismic data using space-varying FIR filters

Array forming in seismic data acquisition can be likened to FIR filtering. Misplacement of the receivers used to record seismic waves can lead to degraded performance with respect to the filtering characteristics of the array. We propose two methods for generating linear space-varying filters that take receiver misplacements into account and demonstrate their performance on synthetic data.


I. INTRODUCTION
Variations in the sampling interval when sampling a signal are often difficult or impossible to prevent.This might make processing the sampled signal problematic, if the related tools can only handle uniformly sampled signals.The motivation for our work comes from the field of seismic data acquisition.A very common practice in this field, is to sum together the (possibly weighted) output signal of multiple seismic receivers.This process is known as array forming and is used to improve the signal-to-noise ratio and to reduce the amount of recorded data.
Being a weighted summation of the output of a finite number of receivers, array forming can be likened to FIR filtering.More specifically, we can view array forming as the application of a linear space-invariant (LSI) filter, since usually the same set of weights is used for the array elements of all arrays in a field.
A problem arises, however, when the receivers are misplaced due to e.g., terrain difficulties.Usually the array weights are designed for a specific geometrical layout of the array elements.When this geometrical layout is violated, it can prove detrimental for the filtering capabilities of the array, as shown in [1].Fortunately, advances in acquisition hardware have enabled us a) to record the output of each individual receiver and b) to know with high (but limited) accuracy the actual location of each receiver.This makes more sophisticated techniques viable for array forming/filtering that can compensate for irregularities in sampling.
A number of solutions have been suggested for the problem of array forming/filtering nonuniformly sampled data, which can be roughly divided in three categories.Methods of the first category interpolate an FIR filter that is designed to filter uniformly sampled data.We shall call this FIR filter the prototype filter.An example is given in [2], where the prototype filter is interpolated to the actual locations of the receivers.These interpolated filter coefficients (or equivalently, the array weights) are then reweighted based on the sampling density at the area around each receiver.We shall refer to this method as geometry-compensating filtering (GCF).The second category involves methods that approximate the outputs of the prototype filter applied to the uniformly sampled data.An example is given in [3], which uses the projections onto convex sets (POCS) framework applied to nonuniformly sampled data.The third category involves reconstruction of the data at the regular sampling locations.The prototype filter can then be applied to the reconstructed data.The methods given in [4], [5] are examples of data reconstruction.
The goal of this work is to propose two methods that generate a linear space-varying (LSV) FIR filter suitable for filtering the nonuniformly sampled data.We will refer to these methods as Method A and Method B. The LSV FIR filter designed by any of these two methods generates the filter output at equi-spaced intervals.In this respect, both methods combine filtering and regularization in one operator.The difference is that • Method A approximates the prototype filter in the spatial domain, which has been already designed, • Method B approximates the ideal response of the prototype filter in the wavenumber domain 1 .In other words, Method B also skips the intermediate step of designing the prototype filter.
Compared to existing works, Method A has a similar flavor as the approach of [2] since they both interpolate a prototype filter to the actual locations of the receivers.However, the interpolation in [2] is driven by the geometry of the receiver arrays while in this paper, the interpolation is based on the band-limited assumption on the received signal, which is also utilized in [4], [5].In reality, such a band-limited assumption is often valid, and the corresponding interpolation process can yield a better approximation to the prototype filter response in the wavenumber domain as will be demonstrated in the paper.

II. PROPOSED ALGORITHMS
For simplicity, we only consider filtering along one spatial dimension x.The continuous data is represented by d(x); the data samples gathered at the N locations of the receivers x j , j = 0, . . ., N − 1 are denoted by d(x j ).The nominal locations of the receivers are defined on the grid xn = n∆x while the actual locations of the receivers are assumed to lie on the denser grid xm = mδx = m(∆x/M ) with M being a positive integer.This is not overly restrictive, since M can be large and the precise receiver locations are known with high, but limited, accuracy.We also define the indicator function s( xm ) to take the value s( xm ) = 1 when a receiver is present at xm and the value zero otherwise.Suppose a prototype LSI FIR filter has already been defined on the nominal grid, whose ith tap is denoted as h(x i ).We assume that h(x i ) = 0 if xi < 0 or xi ≥ L f ∆x, where L f is referred to as the spatial support of the FIR filter.

A. Method A
The filter to be designed in Method A is represented as a set of LSV FIR filters, whose filter taps l ( xm ) for m = 0, 1, . . ., N M − 1 are defined on the dense grid.The spatial support of g l ( xm ) is the same as that of h(x i ), therefore [4] imposes a different FIR constraint on the filter).We desire that the output of one such filter at output location xl should be identical as if the prototype filter were applied on uniformly sampled data.In other words, The output locations xl lie on the nominal grid, i.e., xl = l∆x = lM δx.Utilizing the band-limited assumption on d( xm ) it can be shown that where N = 2P + 1 (a similar expression can be derived when N is even).We can exchange the order of summation and arrive at Substituting (3) in (1) yields We rewrite (4) in matrix-vector form as where h l is an N ×1 vector with h(x l − xn ) as its n-th element; g l is an N M × 1 vector with g l (x l − xm ) as its m-th element; S is an N M × N M diagonal matrix with s( xm ) as its m-th diagonal element; Q is an N M × N matrix with sincd(N ; xo , xm ) as its (m, o)th element and d is an N × 1 vector with d(x n ) as its n-th element.A sufficient condition for (5) to hold is , for which a suitable g l can be found by solving the following least-squares problem where gH l is formed by removing its elements corresponding to the zeros of S. Similarly, Q is constructed after removing the rows of Q corresponding to the zero columns in S.This eliminates S from (6).In order to limit the spatial support of g l , we remove the elements of gl that correspond to elements m of g l for which m < lM or m ≥ (l + L f )M holds, thus forming gl .The corresponding rows of Q are removed as well, to form Q. The problem has a closed-form solution given by A different FIR filter g l has to be calculated for each output location xl .The solution can be seen as a composition of two operations: the first operation is h H l QH = ( Qh l ) H , which can be interpreted as an interpolation of the filter to the actual locations of the receivers.The second operation ( Q QH ) −1 deconvolves the effects of nonuniform sampling.

B. Method B
In comparison to Method A, Method B does not rely on the knowledge of the prototype filter h l , which needs to be pre-designed.To this end, we will try to approximate the behavior of the prototype filter in the wavenumber domain.As a first step, let us use a circular convolution operator to describe the target LSV FIR filter as well as the prototype filter.Accordingly, (1) should be adapted to the form where g l ( xm ) stands for the LSV FIR filter from Method B. In deriving the above, we have used the assumption that the output locations xl lie on the nominal grid, i.e., xl = l∆x = lM δx.Just like Method A, we require that g l ( xm ) have a spatial support in the interval [0, L f ∆x).In other words, To get rid of the dependence on h(x i ) in (8), we resort to the wavenumber domain.Suppose the wavenumber response of the prototype filter is known as h w (k) with −π ≤ k ≤ π.If we use H w to denote an N × N diagonal matrix whose nth diagonal element is given by the sample h w (k) at k = −π + n 2π N for n = 0, 1, . . ., N − 1, then the wavenumberdomain counterpart of (8) can be expressed as where we have coined an N ×N M matrix G to represent the circular convolution of the filter with the (l, m)th element of G given by g l ( x(lM−m)mod(NM) ); F stands for the N -point DFT matrix, and S and Q are defined in (5).A sufficient condition for (9) to hold is The right-hand side of (10) can be broken down to three parts: • G is the LSV filter.Its output is defined on the nominal grid and its input on the dense grid.• F acts on the columns of G to produce the DFT of the LSV filter g l (x l − xm ) with respect to xl .• SQF H acts of the rows of G to produce the nonuniform DFT of g l (x l − xm ) with respect to xm .The product of FG SQF H is a wavenumber connection matrix 2 [6], which can be viewed as the wavenumber response of the LSV FIR filter.Note that the off-diagonal elements in FG SQF H are in most cases non-zero, which means that the relation between the spectra of the data before and after filtering in the wavenumber domain is not simply an elementwise multiplication for nonuniform sampling.
From (10), we formulate the following least-squares problem min where W is a weighting matrix that can apply an individual weight to each element of H w − FG SQF H with the symbol denoting the Hadamard (element-wise) product.Including such a weighting matrix is beneficial, for instance, to enable a better trade-off between different approximation errors in the passband, stopband as well as the "do-not-care" (transition) zones.
As in Method A, the matrix S can be eliminated by removing all the columns of G and rows of Q with the same index as the elements of the diagonal of S that have the value zero.Let G and Q respectively be the reduced form of G and Q resulting from this column-and row-removal.Then (11) can be rewritten as 2 In [6] its continuous counterpart is called "frequency connection function".We use the term "wavenumber connection matrix" for consistency with the rest of the terminology in this paper.This problem can be restated in its vectorized form as Here vec(A) returns a column-vector that stacks the columns of the matrix A. The function diag(v) returns a diagonal matrix with the vector v on its main diagonal.
Using the identity vec(ABC) = (C T ⊗ A) vec(B), where ⊗ denotes the Kronecker product, (13) can be rewritten as where * denotes the complex conjugate of a matrix.Let The elements of g and rows of Ũ that should be removed due to the limited spatial support of the filter g l ( xn ) are given by the indexes of those elements of g that correspond to the zero elements of g l ( xm ).If we call g and Ũ the results after the corresponding row and element removal, the solution is given by G can be constructed from the elements of g and can be applied to the nonuniformly sampled data.

III. RESULTS
The performance of Method A and Method B was examined using synthetically created seismic data.A portion of the spectral content of the synthesized data can be found at the lower part of Fig. 1.The reflections from the Earth's subsurface appear mostly on the lower part of the wavenumber spectrum in the range −0.02m −1 ≤ k ≤ 0.02m −1 .This region appears highlighted in all figures.All the plots have been smoothed by a 5-tap moving average filter.The peaks found at k = ±0.055m−1 are due to waves propagating along the Earth's surface.Their presence is often undesired and should be removed before further processing of the seismic data.Notice that when the data are not uniformly sampled, the wavenumber content appears smeared.This is a well-known side-effect introduced by nonuniform sampling [7].
We generated 100 different realizations of the receiver locations and filtered the nonuniformly sampled data using Methods A and B. Method A approximates the prototype filter whose wavenumber response is given in the upper part of Figure 1.In Method B, the prototype filter is not given but is supposed to have a similar passband and stopband region in the wavenumber domain as in Method A. The LSV FIR filters resulting form Method A and B have the same spatial support as that of the prototype filter.The average spectrum of the 100 realizations can be seen in Fig. 2. In addition to Methods A and B, an adapted version of the GCF was implemented 3 .
In addition, we also compared the proposed methods against 1) the ideal case where the prototype filter is applied on the uniformly sampled and 2) the case where the prototype filter is applied directly on nonuniform data.In the latter case, the spectral leakage induced by the irregularities is most obvious in the passband.In Fig. 2 we see that using the LSI filter on nonuniformly sampled data gives an output that differs almost 15dB at the edges of the passband.In contrast, Methods A and B give a filtered output that, on average, is closer to the ideal case in the low wavenumbers, exhibiting less than 5dB maximum difference from the uniformly sampled data case (Fig. 2) in the passband region.This is due to the fact that Methods A and B compensate for irregularities in sampling.The attenuation in the stopband, however, is less when using Methods A and B. This is because the FIR filters generated for each individual output will not, in general, have zeros at the same locations of their wavenumber responses.This leads to a more flat response in the stopband.GCF is somewhere in the middle, as it interpolates the filter to the locations of the receivers and compensates for sampling density, but does not deconvolve the effects of nonuniform sampling.
The standard deviation of the output spectrum can be seen in Fig. 3. Methods A and B exhibit in this case a significantly smaller standard deviation in the lower wavenumbers, for example, in Fig. 3, around 10dB lower than using the LSI filter.This means that Methods A and B may provide an output that is, generally, stable.

IV. CONCLUSION
We proposed two methods for generating LSV FIR filters suitable for filtering nonuniformly sampled data.The resulting filters yield a more accurate output in the passband than simply 3 The method in [2] works on data having two spatial dimensions and is adapted here to the single spatial dimension case.applying an LSI filter directly to nonuniformly sampled data.The output is also more stable, as it varies less for different realizations of the receiver locations.

Fig. 1 .
Fig. 1.Upper part: the wavenumber response of the prototype filter.Lower part: the wavenumber content for the data.