Digital Analysis of Geophysical Signals and Waves

5.2 Filtering of Electromagnetic Fields

Because of the sensitivity of the Earth’s electromagnetic (EM) properties to temperature, fluids, ore minerals, and lithologic variations there is a remarkably broad application of EM techniques in geophysics. They range from monitoring environmental problems, exploring for energy and mineral occurrences, studying the deep interior structure of the Earth (and other planets), and investigating near-surface regions of potential natural and man-made hazards (e.g., earthquakes, volcanic eruptions, unexploded ordnance).

The governing equations in EM geophysics are the classical Maxwell’s equations. They, together with the so-called constitutive relations, and appropriate boundary conditions, determine the behavior of electric and magnetic fields in matter. All EM geophysical measurements are recorded in the time (t) domain and space domain, (x, y, z). In some EM methods the processing and analysis is performed in the time-space domain whereas in others these steps are done in the temporal frequency domain. Thus, EM methods in geophysics are referred to as time domain EM (TDEM) or temporal frequency domain EM (FDEM) even though as we emphasize, the measurements are always made in the time domain. Therefore, it is of utmost importance to appreciate and understand EM geophysical signals in time, space, and frequency domains. Our discussion of the filtering of potential fields in Section 5.1 focused on the space-spatial frequency domains. This section on EM fields will address the time-temporal frequency domains only so we’re going to simply use frequency to mean temporal frequency in our discussion.

EM geophysical signals are generated artificially by radar and inductively or directly coupled transmitters. Natural EM sources are of two major types: 1) complex interactions between the time-varying solar wind and the Earth’s geomagnetic field and 2) world-wide lightning discharges. The frequencies used in EM geophysics extend over 20 decades from below 10^-4 Hz to beyond the visible spectrum. Lower frequencies penetrate deepest into the Earth and allow penetration to many 10s of km.

Maxwell’s Equations

We do not intend to present the principles of EM theory as embodied in Maxwell’s equations. Two good treatments aimed specifically at geophysics are found in Ward and Hohmann (1988) and West and Macnae (1991). What we want to do here is apply the principles presented in previous sections of the DAGSAW web site to relationships in EM that are not often couched in terms of the filtering operations that they really are. We believe that this will provide insight beyond the conventional presentations. You will be the judge.

First we list in the table below the EM quantities that we will use and their units. We list them in the order they appear in our discussions below. You’ve seen a table like this before if you’ve taken any physics class with EM. Capital letters are almost universally used for most symbols, irrespective of whether or not they are in the time or frequency domain. Bold quantities denote vectors. As we said in Section 3.1, in spectral analysis it is customary to use lower case symbols for the time (or space) domain and upper case for the corresponding functions of frequency. The distinction can be made by using capital letters only and explicitly including t or f as arguments, e.g., E(t) and E(f) would be the electric field Fourier transform pair E(t)E(f). We have elected to use lower and upper case symbols, for time and frequency domains, respectively, as done by Ward and Hohmann (1988). We’ll also explicitly include the independent variables for additional emphasis so the electric field Fourier transform pair is e(t)E(f).

For illustration we will only consider the two Maxwell’s equations that deal with time-varying fields (the other two equations deal with steady or static fields). The appropriate equations are written in differential form in the usual way as:

and

Of course, the nonzero curl operations, e. g., tells us that a field, in this case E, circulates about a source that generates it. The source of the curl isin equation 5.2.1. and is in equation 5.2.2. Visualizations of the circulations described by Maxwell’s curl equations are animated in Figure 5.2 following sketches presented by West and Macnae (1991).

Following our adopted notation procedure, since Maxwell’s curl equations 5.2.1 and 5.2.2 are clearly in the time domain, i.e., they contain time derivatives, they are written as

and

We have added the position vector, r to denote (x, y, z) space domain dependency too. Those of you who remember your EM theory know that equations 5.2.1 and 5.2.6 are expressions of Faraday’s law. Equations 5.2.2 and 5.2.7 are Ampere’s law with Maxwell’s displacement current density term,added to the original conduction current density, J. Using our knowledge in Section 4.5 of the expression for the derivative operation in the frequency domain, we can take the Fourier transform of Maxwell’s equations 5.2.6 and 5.2.7 yielding

and

Now, let’s look at the constitutive relations which are frequently substituted into Maxwell’s equations. You’ve undoubtedly seen them as:  

and

In a seminal, but unappreciated, paper by Fuller and Ward (1970) the point is made that a frequency dependency of μ, σ, and ε means the equations 5.2.7, 5.2.8, and 5.2.9 are in the frequency domain, i. e.,

and

Here, we recognize from the convolution theorem (Section 4.7) that such multiplications in the frequency domain are filtering operations, therefore, they are convolutions in the time domain. So, in the time domain, the constitutive relations are

and

where the corresponding impulse response functions are the EM properties μ(r,t), σ(r,t), and ε(r,t) as as sketched in Figure 5.3 using the “black box” filtering description. The system response (or transfer) functions for these filters are μ(r,f), σ(r,f), and ε(r,f), respectively, in equations 5.2.10 to 5.2.12. Note that we have not used lower and upper case symbols for these EM properties. Lower case Greek symbols are used in both the time and frequency domains; they are italic in the time domain.

 

As pointed out by Fuller and Ward (1970), equations 5.2.13 - 5.2.15 have interesting ramifications for the time domain Maxwell’s equations 5.2.6 and 5.2.7 which, after substitution of the appropriate constitutive equations, take on the forms

and

Clearly, Maxwell’s curl equations in the time domain must be expressed with these convolutions if μ, σ, and ε are functions of frequency. We also confirm, from Section 3.3, that since real quantities in the time domain are generally complex in the frequency domain, the possibility of frequency dependent, complex μ, σ, and ε follows. The common case where μ, σ, and ε are assumed to be real constants, independent of frequency, is discussed by Fuller and Ward (1970). Since the inverse Fourier transform of  a real constant is a delta function (Figure 3.7), this means that μ(r,t), σ(r,t), and ε(r,t) are all delta functions and the convolutions in equations 5.2.16 and 5.2.17 yield the familiar form of equations 5.2.3 and 5.2.4.

Magnetotellurics

To emphasize how ubiquitous filtering is in EM geophysics we turn to the magnetotelluric (MT) method introduced briefly in the Preface, Section 1.1 (Figure 1.7). MT employs a frequency range that spans both high frequency lightning (1 f 10⁴ Hz) and low frequency (1 f 10^-4 Hz) solar-induced, magnetospheric sources (Figure 5.4). Therefore, MT has shallow to deep applications throughout the entire Earth’s crust and upper mantle (from ~10 m to 100 km). MT has been used in the SAGE program to study midcrustal conductive zones, geologic structure, and groundwater resources. Figure 5.5 shows what it looks like when SAGE students happily make EM measurements in north-central New Mexico.

With the MT method there is immediate Fourier transformation from time domain recordings into the frequency domain. This is accomplished using the FFT (discussed in Section 3.6) or an equivalent method. We will not dwell on the details but, rest assured that the issues inherent in the time and space sampling theorem (Section 2.3) and finite length data (Section 2.4) are properly addressed within the context of the discrete Fourier transform as discussed in Section 3.5. For the present we will look at important relationships inherent in MT that are easily recognized as filtering operations. As we were told in the previous section, the underlying interrelationships are given by Maxwell’s classical equations, the constitutive relations, and the appropriate boundary conditions. The three EM properties of a material listed in the table above, namely the magnetic permeability, μ; electric conductivity, σ; and dielectric permittivity, ε as functions of x, y, and z are what we ultimately seek in EM measurements. These properties are always assumed to be constants over the MT frequency band. The importance of these properties is emphasized when one recalls that the velocity, c of light in free-space is equal to 1/(ε₀μ₀)^1/2 where the subscript 0 indicates free-space values. One could argue that the intrinsic or characteristic impedance of free space, Z₀ which equals the ratio (μ₀/ε₀)^1/2 is of equal importance. Unlike the relations for free-space where the electric conductivity σ = 0, velocity and intrinsic impedance of Earth materials are functions of all three EM properties.

Simply stated, MT is the study of interrelations of the natural, time varying electric and magnetic fields at the surface of the Earth. Usually these fields are measured in orthogonal directions, x and y at selected locations over time windows lasting from seconds to several days. The time records of the electric, e(t) and magnetic, h(t) (or b(t)) fields at a recording site, as presented in Figure 5.4, are individual samples of all possible (past, present, and later) records of an infinite, unknown random process inherent in solar wind fluctuations and lightning.

Normally, a single sample time record would not be suitable for representing an entire random process. However, in some instances it is possible to derive meaningful statistical information about the random process from the analysis of a single, arbitrary sample record. Such information could be spectral density functions (Section 3.2), autocorrelation and cross-correlation functions (Section 4.6), or energy density functions (Section 4.6). But, it is important that there exists a measure of MT signals that is statistically stationary so measurements will yield EM Earth properties that do not vary depending on the time the measurements are recorded. Clearly the electric or magnetic fields alone do not meet this requirement although they probably are stationary in the time average, i. e., all records measured at different times have zero mean when averaged over long enough time windows.

Homogeneous Earth Impedance. As it turns out, the “trick” needed for stationary results in MT is to take the ratio of the surface electric to magnetic fields in the frequency domain. Moreover, it’s easy to show that for a homogeneous, isotropic Earth, the ratio of these orthogonal fields in the frequency domain is identical to the intrinsic impedance, Z(f) of the material:

This simple equation (and its extension below to multi-dimensional and anisotropic Earths) is the most important equation in MT! Inherent in equation 5.2.18 are the assumptions that the MT source fields are plane waves normally (vertically) incident on a planar Earth. These assumptions have been experimentally verified to be correct enough in all but obvious circumstances, e.g., near lightning strikes, in high magnetic latitudes during auroral activity, or when a spherical Earth model must be retained.

Let’s expose the true nature of equation 5.2.18 by simply writing it as

Since we are now keenly aware what a multiplication in the frequency domain means, we immediate recognize this as a filtering operation. It’s time domain convolution equivalent is

We remind you that we are continuing the custom of using lower case expressions for the time domain variables even though this is rarely observed in EM literature. Ward and Hohmann (1988) is one exception.

In filtering language, equation 5.2.20 defines a linear system where the magnetic field, h_y(t) is the input signal; the electric field, e_x(t) is the output; and z(t) is the impulse response function. Z(f) in equation 5.2.19 is the system response (or transfer) function in this case so the impedance is the “black box” whether it’s expressed in either time or frequency domain. This is schematically shown in Figure 5.6 in the time domain. Fascinating don’t you think?

We said in Section 4.5 that when we know the impulse response function (or the system response or transfer function) for a filtering operation we know everything about it, so it’s not surprising that the intrinsic impedance in EM has this role. The impedance in the frequency domain is the ultimate MT result from which all analysis proceeds. We, of course, know that Z(f) is complex so we can express it as an amplitude and phase or as real and imaginary parts (Section 3.2). Traditionally in MT, the impedance phase and a proxy for the impedance amplitude are used for plotting and analysis purposes. The function of the impedance amplitude that used is called the apparent resistivity, ρ_a given by

A representative plot of the apparent resistivity and impedance phase versus period (the reciprocal of frequency) from a single site in the SAGE field area in New Mexico is presented in Figure 5.7.

We need not concern ourselves here with the derivation of equation 5.2.21 other than to say that it results from an assumption that the Earth is homogeneous and isotropic. If the Earth was homogeneous and isotropic then the apparent resistivity would be the Earth’s true resistivity and the apparent resistivity curve in Figure 5.7 would be a constant value. Derivation of equation 5.2.21 also assumes two more things: 1) the magnetic permeability is that of free space, μ₀ since μ rarely differs much from this in Earth materials and 2) the electric permittivity, ε can be neglected because the conductivity property dominates. These assumptions result in the final objective of MT measurements to be the three-dimensional (3-D) distribution of conductivity usually expressed as its reciprocal, the resistivity, ρ(x, y, z). Going from the phase of Z(f) and ρ_a(f) to a 3-D Earth resistivity model requires forward and/or inverse modeling which will not be discussed here.

Homogeneous Earth Boundary Conditions. By applying the boundary conditions that the tangential electric and magnetic fields are continuous across a layer interface we determine the partitioning of incident, i; reflected, r; and transmitted, t fields at a boundary. For a plane EM wave, normally incident on a homogeneous, isotropic Earth with a planar air-Earth interface, the reflected and transmitted electric and magnetic fields in the frequency domain obey the relations:

and

The subscripts 0i, 0r, and 1t on the E and H fields first refer to the medium, either 0 (free space) or 1 (Earth), respectively. Then i, r, and t denote the incident, reflected, and transmitted field amplitudes, respectively. Z₀(f) and Z₁(f) are the characteristic impedances of free space, 0 and the Earth, 1. Now, since we recognize that all of these equations are multiplications in the frequency domain we know they’re filtering operations. The complex transfer functions in equations 5.2.22 to 5.2.25 above are obviously functions of the characteristic impedances of the two media. As such, they tell the whole story relating the incident, reflected, and transmitted fields as transfer functions must do. For example, interpreting equation 5.2.22 as a filtering operation means that in the frequency domain, the incident electric field, E_0i (f) is the input and reflected electric field, E_0r(f) is the output. The “black box” transfer function that performs this input-output relationship is called the reflection coefficient. Using the notation RC_E(f) for the electric field reflection coefficient we see that it is, as expected, equal to the ratio of the reflected to incident electric field amplitudes in frequency domain i. e., 

Of course, we can use the reflection coefficient notation to inverse Fourier transform equation 5.2.22 and write it in the time domain as a convolution,

Figure 5.8 presents both the frequency and time domain “Black box” views of equations 5.2.22 and 5.2.27.

We now rewrite equations 5.2.23 to 5.2.25 (as we did for equation 5.2.22 with equation 5.2.26) using reflection, RC and transmission coefficients, TC of the electric, E and magnetic, H fields to yield:

and

Along with equation 5.2.26, there is a wealth of very neat, practical information in these equations. For example, because Z₀(f) = 376.6 ohms which is >>|Z₁(f)| for the Earth at MT frequencies (often by orders of magnitude) you should be able to convince yourself that the following approximations are justified: RC_E ~ -1, TC_E ~0, RC_H~+1, and TC_H ~+2. If these relations were strictly true, rather than approximations, there would be no MT method because there would be total reflection of the incident electric and magnetic fields with no EM penetration into the Earth. In fact, the Earth is nearly a mirror to MT fields but a tiny fraction does penetrate and can be used to probe the Earth electromagnetically, often to many 10s of km depth.

Figure 5.9 presents such a result from New Zealand where non-homogeneous Earth electrical properties are imaged far below the occurrence of earthquakes. This reminds one of us when upon describing to an electrical engineer colleague what MT was, he said that what were doing was impossible! Here’s to the impossible.

Non-Homogeneous Earth Impedance. The Earth need not be entirely homogeneous and isotropic for an identical form of the reflection and transmission equations above to be completely valid. A case in point is that of a an isotropic, one-dimensional (1-D or layered) Earth where the complex impedance as defined by equation 5.2.18 theoretically contains all the necessary information for modeling to obtain the layered resistivity values as a function of depth, ρ(z). This impedance measured at the surface in the general case of homogeneity or 1-D is called the surface impedance.

In the case of multi-dimensional (2-D or 3-D) or anisotropic 1-D Earths the simple relation in equation 5.2.19 does not hold. Experimentally the EM fields are not plane waves but there is still linear system coupling between the horizontal EM fields expressed by

or using tensor notation as

Now, there are four complex, frequency domain impedance elements, Z_xx, Z_xy, Z_yx, and Z_yy that make the life of an MT researcher much more difficult. They describe the 2 by 2 complex, tensor impedance in equation 5.2.32; the first subscript refers to the electric field and the second to the magnetic field as related in equation 5.2.31. Each of the equations in 5.2.31 can be represented as a two-input/single-output linear system filter, i. e., in the time domain we have the convolutions:

The “black box” view for the single output, ex(t) with two inputs, hx(t) and hy(t) in equation 5.2.33  is given in Figure 5.10.

Since the Earth is always 3-D to some extent, the coupled equations in equation 5.2.31 must be solved for four impedance elements each yielding corresponding apparent resistivities and impedance phases. But, since there are only two impedance equations with four unknown impedance values in equations 5.2.31, we have an under-determined system that appears to be unsolvable. As it turns out, at least two more independent equations of the form of 5.2.31 will be available if multiple time-series signals are recorded and they have different EM source polarizations. In fact, it is statistically desirable to have many recordings from different EM source polarizations so averaging can be used to reduce measurement noise. This actually yields an over-determined system of equations at each recording site. In this case, least-squares estimates of the four impedance elements involve auto- and cross-powers (Section 4.6), averaging over both multi-records and frequency bands, and other noise-reducing schemes (e.g., Sims et al. 1971; Gamble et al., 1979). These procedures are excellent applications of digital signal processing but they’re beyond the scope of this introduction to the MT method.