Digital Analysis of Geophysical Signals and Waves

5.1 Filtering of Potential Fields

In this section we apply what we have learned to specific geophysical signals and waves. Let’s begin with the signal that we have used repeatedly from the beginning of the web site, namely the magnetic profile (Figure 1.5) that was measured at the Los Alamos National Laboratory (LANL) during the SAGE program. We will perform several filtering operations on this data set; first let’s look at the upward continuation operation. As discussed in Section 4.5, this operation is typically applied to areal potential field data such as magnetic and gravity surveys measured on an x-y grid. Such 2-D data requires 2-D filtering as discussed in Section 4.5. However, a simplification occurs if the areal geophysical data are produced by strictly 2-D subsurface source distributions, for example, of density in the case of gravity anomalies. This can never truly occur in the actual Earth since the distributions would need to be invariant from –infinity to + infinity in a constant direction (the strike direction). However, some geological features such as ridges, linear fault strands, volcanic dikes, and rift structures can be approximately linear over survey-scale distances. In such cases the geophysical signals approximate 2-D and they are constant in some direction, let’s say the y direction. Then the 2-D signal, s(x, y) can be expressed as s(x) multiplied by one. The 2-D Fourier transform of s(x,y) in this case is δ(f_y)S(f_x).

Given this situation, i.e., where the anomaly patterns are quasi-linear in appearance (implying 2-D source distributions) data processing can precede using one-dimensional (1-D) Fourier analysis. We shall assume that this is true for the magnetic LANL SAGE data and use it to illustrate filtering of 1-D geophysical potential data.

Upward Continuation

Figure 5.1 interactively presents the filtering sequence in both space and time domains to perform upward continuation by 1-m of the LANL SAGE magnetic data. These input data, s(x) are assumed to be caused by an anomalous 2-D source field superimposed on a constant geomagnetic background field. For these data, N = 32, L_x = 17.6 m, and Δx = 0.55 m, therefore, Δf_x = 1/L_x ~ 0.057 cycles/m. The appropriate filter (the system response or transfer function) for upward continuation by h = 1-m from Section 4.5 is

To “walk through” the filtering operations using Figure 5.1 begin by clicking on the FFT arrow to perform the fast Fourier transform of the LANL SAGE input magnetic data, s(x). Now, the right-hand column under Frequency Domain presents the real and imaginary parts of the FFT, S(f_x) exactly in the order in which they are listed in FFT output. Instead of explicitly labeling the frequencies f_x = 0, + 0.057, + 2 x 0.057 , … on the frequency domain plots in Figure 5.1 we have simply noted the “frequency index” which runs from 0 through 31. For the LANL SAGE data the actual sequential values of frequency in the FFT output (Section 3.6) appropriate to these indices are listed in the table below:

The first value plotted in Figures 5.1a and b is the zero frequency (or DC) result. Then the positive frequency values increase from +1/L_x in increments of Δf_x = 1/L_x up to the Nyquist frequency, fN = + 1/2Δx. These values are followed by the negative frequency results in reverse order as discussed in Section 3.6 beginning with a frequency equal to Δf_x less than the Nyquist frequency, i.e., - (N-1)/2L_x. The negative frequency values are then incremented by Δf_x = 1/L_x up to - 1/L_x.

At this point we pause to point out several features of the FFT that were discussed earlier in the web site:

1. First, careful comparison of the real and imaginary FFT values of S(f_x) confirms the symmetry properties given in Section 3.3, That is, the real part of the Fourier transform of a real function has even symmetry and the imaginary part has odd symmetry.

2. We also understand why the real DC component of S(f_x) has a very large value in this case (917,000). It is because the zero frequency value is the integral (as explained in Section 3.1) of the entire input magnetic curve and the area under this curve, in units of nT-m, is large because all of the values exceed 50,000 nT (Figure 5.1a).

3. We now prove our observation in the preface of the web site (Section 1.1) that there is a relatively large spatial frequency component in the magnetic data at ~1 cycle/4 m = ~0.25 cycles/m. This is evident by the large spectral component in the real part of S(f_x) in Figure 5.1b at the frequency indices of 4 and 28. These have spatial frequencies of +4Δf_x ~ 0.23 cycles/m.

4. Finally, the fact that the real and imaginary values of the FFT are small in Figures 5.1b and c at high frequencies (at and near the Nyquist frequency) is evidence (but not conclusive) that the digitized magnetic signal is not seriously aliased. Higher frequency components in an unaliased signal are folded into this frequency range (and lower) as discussed in Section 2.3. Normally one does not know how serious an aliasing is but in our case we can see from Figure 2.4 that N = 32 sample points with a Nyquist frequency of 1/2Δx = 1/~0.91 cycles/m appears to adequately avoid serious aliasing.

 

At this juncture filtering of the input LANL SAGE magnetic data can precede as a convolution in the space domain or as a multiplication in the frequency domain as described in Section 4.4. Clicking on the filtering bar presents these options for 1-m upward continuation .

In the space domain the filtered output is obtained by a convolution () between the input data, s(x) (Figure 5.1a) and the 1-m upward continuation impulse response function, h_u(x) (Figure 5.1d). The filtered output, y(x) appears in the lower left-hand panel, Figure 5.2f.

To obtain the result in Figure 5.1f via spatial frequency domain filtering, the right-hand column first presents a single multiplication () between the complex Fourier transform, S(f_x) of the input data, s(x) (Figures 5.1b, c) and the 1-m upward continuation system response function, H_u(f_x) (equation 5.1.1 and Figure 5.1e). The real, exponential values of the 1-D upward continuation filter in the frequency domain (equation 5.1.1) are presented in Figure 5.1e in the order appropriate to the table above since this is necessary when using the FFT and IFFT. The complex result, Y(f_x) of the single multiplication appears in the lower two right-hand panels in Figure 5.1g and h. Y(f_x) is then inverse Fourier transformed, IFFT to yield the desired filtered, output result in the space domain, y(x) in the lower left-hand panel, Figure 5.1f. Remember, from our discussion in Section 3.5, that the consequence of digitizing in space and frequency domains is to produce periodic results in the other domain. Therefore, the plots in Figures 5.1f, g, and h must be regarded as one period of periodic functions in the space and frequency domains.

The filtered upward continued output (Figure 5.1f) meets our intuition of what a magnetic anomaly should look like farther away from the sources. That is, the output is smoother and the numerical values are less than the input (Figure 5.1a) because the filter in the frequency domain is a form (an exponential) of low pass filtering.  

Filtering in the frequency domain may appear to be a roundabout way of accomplishing a filtering task compared to the direct convolution operation in the space domain. The sequence of shifting then multiplying for all values of x (or t) in the convolution operation is computationally intensive especially when the number of points becomes large. This compares to a single multiplication in frequency domain filtering. Two key factors favor the frequency domain approach: 1) the FFT and IFFT necessary for  the frequency domain operation are computationally very fast operations and 2) design of filters is usually more direct and intuitive in the frequency domain.