Chapter 5: Fourier-domain Modeling & Transformations

Jiajia Sun

5 Chapter 5: Fourier-domain Modeling & Transformations

5.1. background

5.1.1. Motivation

Up until this point, we have dealt with 2D gravity and magnetic data maps in the spatial domain. In this chapter, we will be looking at different way in dealing with these data: in the frequency domain. A lot of data processing techniques require the data to be in frequency domain due to its simpler and more efficient mathematical operators. The motivation to studying Fourier domain modeling for modeling any potential fields data are the following:

The expression of potential field in frequency domain gives us a different perspective to look at our data
Much of the processing and analysis of potential field data are done in frequency domain
Potential field is related to the source distribution by a convolution with the Green’s function (will be explained later on)

In this chapter, we will examine Fourier expression of simple potential functions and Fourier expression of field due to simple sources.

5.1.2. frequency and phase

Before going over Fourier Transforms, we first need to know the basic theory behind it. There are several basic concepts that we will first cover before going through Fourier transforms.

5.1.2.1. Waveforms

A spring-mass system, as we see in the figure below, oscillates throughout time. The displacement of the mass from the oscillatory behavior can be graphed with respect to time as we see in the figure. It turns out this movement can be expressed as a sinusoidal function. We call this a waveform.

The amplitude of the waveform depends on the constant before the sin function. For example, the above figure shows an waveform with an amplitude of 2 units. The blue value in the below indicates the amplitude.

f(t)= 2.0 *sin (2* π*2*t)

If there is a function that has the form of f(t) = 1.0*sin(2*π*2*t), then the waveform will be exactly the same as the aforementioned example, except that the amplitude would be halved, as we can see in the blue line in the following figure:

The period of the waveform depends on the value preceding t inside the sin function . One period represents the amount of time it takes for the waveform to undergo one cycle (1 peak and 1 trough). For example, the red text in the following waveform determines the period:

f(t) = 2.0 * sin(2* π *2*t)

If there is a function that has the form of f(t) = 2.0*sin(2*π*1*t), then the waveform will be exactly the same as the aforementioned example, except that the period is twice of that, as we can see in the blue line in the following figure:

5.1.2.2. Frequency

The frequency of the waveform can be determined by the period. Frequency represents the number of cycles in one second, and the unit is (/s) or Hertz(Hz). Frequency can also be mathematically defined as:

f = 1/T

Where T is the period.

If we have waveforms that have their peaks and troughs closer together, we say that the waveform has a high frequency. If they are far apart, then the waveform has low frequency.

Test your understanding: frequency

What is the frequency of this waveform? What about the period?

5.1.2.3. Phase

Phase is determined by the argument inside the sin function. For example, the following figure shows two waveforms of different phases. The first function below is for the red line in the graph, and the second one is for the blue line in the graph.

f(t) = 2.0 * sin(2 * π * 2 * t)

Notice that the two waveforms have the same amplitude and frequency, but one seems to be delayed compared to the other. This is because these two waveforms have different phases.

Phase is the argument inside the sine(or cosine) function for a waveform. Given a fixed amplitude and frequency, phase determines when the peaks (or troughs) occur. Simply, it specifies where in the cycle is it oscillating at t=0. The unit for phase is in degrees or radians.

There are two terms that are related to phase: in-phase and out-of-phase.

If the two peaks of two signals with the same frequency are in the exact same alignment at the same time, then they are called in-phase. If the two peaks of two signals with the same frequency are not in an exact alignment at the same time, then they are out-of-phase. The following figure gives out a visualization of the two difference.

The important concept that we want to take home from this is the fact that a sine wave and a cosine wave are 90 degrees out-of-phase with each other. Thus, we can express a waveform in two ways: using a sine wave or a cosine wave.

5.1.2.4 General notation

The more general notation of a waveform can be expressed as:

$f(t) = A sin(\omega t + \phi)$

or

$f(t) = A cos(\omega t + \phi)$

Where A is the amplitude, and $\omega$ is the angual frequency in radians per seconds. The angular frequency can be found from the equation:

$\omega = \frac{2\pi}{T} = 2\pi f$

Main takeaways from this frequency and phase

Any periodically oscillating sinusoidal wave can be characterized by three fundamental properties

Amplitude(A)
Frequency( $\omega = 2\pi f$ )
Phase( $\phi$ )

A wave can be expressed using the following notation:

f(t) = A sin(ωt + φ)

or

f(t) = A cos(ωt + φ)

5.1.3. complex variable

A complex variable can be expressed as such:

$z= a+ib$

where a and b are real numbers, and i is the imaginary unit that is equal to $\sqrt(-1)$ , thus $i^2 = -1$

Thus, a complex number can be separated into two:

a: real part
b: imaginary part

A complex number can be viewed as a point in the complex plane, as the following figure:

5.1.3.1 Absolute value and argument

We know that a complex number can be expressed as such:

$z = a + ib$

The absolute value(modulus) of this complex number is thus:

$|z| = r = \sqrt(a^2 + b^2)$

The argument(phase) is the angle between the vector and the positive real axis. This makes more sense from the complex plane representation in 6.1.3.

The argument is thus expressed as:

$arg(z) = \theta = atan2(y,x)$

where y is the projection of the complex number in the imaginary axis, and x is the projection of the complex number in the real axis. Therefore:

$a = r cos(\theta), b = rsin(\theta)$

and thus:

$z = rcos(\theta) + i rsin(\theta)$

Test your understanding: Complex Number

Given a complex number:

$z = 1 + i\sqrt(3)$

Calculate its modulus and argument!

5.1.3.2. Euler’s formula

Euler’s formula is the fundamental relationship between the trigonometric functions and the complex exponential function. This can be mathematically expressed as:

$e^{i\theta} = cos(\theta) + i sin(\theta)$

From Euler’s formula, we can express the complex number as:

$z = rcos(\theta) + i rsin(\theta) = re^{i\theta}$

Therefore, z can also be expressed as:

$z = re^{i\theta}$

In other words, a complex number can be expressed as either with a trigonometric function or an exponential function.

Test your understanding: Euler’s Formula

Let $\phi = \frac{\pi}{2}$

Calculate $e^{i\theta}$

Hint: Just plug in $\theta = \frac{\pi}{2}$ to the function $e^{i\theta} = cos(\theta) + i rsin(\theta)$

Ans: i

Consider that a mathematical notation for a cosine wave is:

$f(t) = A.cos(\theta)$

where A is the amplitude and $\theta = \omega t + \phi$ is phase.

A different way to to express that wave is to use the following notation:

$f(t) = Ae^{i\theta}$

The representation of this in the complex plane is as follows;

(Note that if you were to see something that looks like $Ae^{i\theta}$ , please always keep in mind that this represents sinusoidal waves with amplitude A and phase $\theta$ )

5.1.3.3. Deriving a complex number

Deriving an exponential function is much simpler than deriving a trigonometric function. Consider the following problem;

$\frac{\partial e^{i\omega t}}{\partial t} = i \omega e^{i\omega t}$

Now, remember that:

$i = e^{i\frac{\pi}{2}}$

Therefore:

$\frac{\partial e^{i\omega t}}{\partial t} = i \omega e^{i\omega t} = e^{i\frac{\pi}{2}}\omega e^{i\omega t}$

$\frac{\partial e^{i\omega t}}{\partial t} = \omega e^{i(\omega t + \frac{\pi}{2})}$

Notice that the argument in the exponential has changed from $i\omega t$ to $i(\omega t + \frac{\pi}{2})$

The conclusion to this is that when you take a time derivative, the phase will change by $\frac{\pi}{2}$ .

5.1.3.4. Why bother with complex numbers?

The first question that you might ask is: why do we need to know about complex numbers? The reason is: It turns out, the whole Fourier theory was built upon complex variables.

Other than that, complex numbers offers a lot of mathematical conveniences, such as:

No need to deal with trigonometric functions
Exponentials are easier to manipulate
Some problems can be solved easier using complex number/complex analysis (examples: quantum physics, conformal transformations, AC circuits, etc.)

If you are interested to learn more about complex numbers, here are some resources that might be of interest:
• Complex Analysis Made Simple on youtube.
• MIT OpenCourseWare: Development of the
complex numbers
• Functions of complex variables in StackExchange
•Intuitive Arithmetic with Complex Numbers (Betterexplained.com)

5.2. Fourier transform

5.2.1. waveforms and amplitude spectrum

Virtually, any real world waveforms can be represented as a sum of sins, no matter how complicated they may look.

For example, take a look on this waveform:

The blue waveform may look complicated, but it is actually the sum of four other sinusoidal functions in the above graph, with the equation for each line being;

Red line: $f(t) = 0.3 sin(\frac{2\pi}{T}t + \frac{\pi}{2})$

Green line: $f(t) = 0.6sin(\frac{4\pi}{T}t)$

Magenta line: $f(t) = 0.2 sin(\frac{6\pi}{T}t)$

Cyan line: $f(t) = 0.3sin(\frac{8\pi}{T}t - \frac{\pi}{2})$

This proves that no wonder how complex they are, waveforms are just a sum of sinusoidal waves.

If we were to plot a waveform that is the sum of two sinusoidal waves, and plot the magnitude of each wave in terms of their frequency, we get something like this:

http://www.continuummechanics.org/fourierxforms.html

The signal results from the sum of two sine waves. These two waves have different amplitudes, frequencies, and phases, and this difference is made more clear in the amplitude spectrum on the right side of the above figure.

Another way to look at the difference between frequency domain and time domain is with the following figure:

The figure above shows two different ways of viewing the combination of these 3 waveforms. In time domain, it is the summation of these three waveforms (with noise), and thus the perceived signal in the time domain looks more complicated. When viewing the three waveforms in the frequency domain, it is in the form of three separate peaks that tells you the frequency and amplitude of said waveforms.

5.2.2. definition

Mathematically speaking, any waveform function can be expressed as:

$f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i\omega t} dw$

The term $F(\omega)e^{i\omega t}$ is a sinusoidal wave, with the amplitude being $F(\omega)$ . This essentially says that any waveform is a sum of many sinusoidal waves with different amplitudes, phases, and frequencies.

The Fourier Transform allows us to find the constituent frequencies given a signal in time domain. This is mathematically expressed as:

$F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-i\omega t} dt$

It decomposes a signal into a series of sinusoidal waves of different amplitudes, frequencies, and phases. This helps identify what the sine and cosine components that make up the signal are.

The inverse Fourier transform is thus:

$f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i\omega t} dw$

In terms of notation, we refer the Fourier transform of a function $f(t)$ as the capital $F$ , and the inverse Fourier transform as $F^{-1}$ . Putting this in mathematical perspective:

Fourier Transform: $F[f(t)] = F(\omega)$

Inverse Fourier Transform: $F^{-1}[F(\omega)] = f(t)$

If you want to know more about Fourier transform(FT), here’s some optional reading material on FT:

5.2.3. from time domain to spatial domain

For now, you might notice that everything that we have talked about so far is in time domain. However, potential field data is in the spatial domain. This does not seem to be a problem because Fourier transform also applies in the spatial domain. The Fourier transform in time and spatial domain can be seen in the following figure:

Notice that they are both the same, except for the fact that the Fourier transform of a in time(t) domain is a function of angular frequency ( $\omega$ ), and in spatial(x) domain, the Fourier transform is a function of wavenumber(k). Thus, everything in spatial domain is the same as in time domain, except we’re dealing with space (x) instead of time (t), and the Fourier transform involves wavenumber (k) instead of angular frequency ( $\omega$ ).

As we mentioned before, frequency is a measure of how many cycles per unit of time. Wavenumber is the measure of how many cycles per unit of distance. The following equation applies for wavenumber:

$k= \frac{2\pi }{\lambda}$

$f =\frac{1}{\lambda}$

Where $\lambda$ is the wavelength (one cycle in unit of distance).

Another concern is that the data maps that we have dealt with are 2-dimensional. Therefore, if we were to take the Fourier transform of such spatial data maps, then we need to do a 2D Fourier transform in spatial domain. The following is such equation:

Note that $k_x$ is the wavenumber in the x axis and $k_y$ is the wavenumber in the y axis respectively.

5.2.4. DC-component

There is one interesting result that we can see from the Fourier transform:

$F(k) = \int_{-\infty}^{\infty} f(x)e^{-ikx}$

If we have k=0, that means the Fourier transform would be:

$F(0) = \int_{-\infty}^{\infty} f(x)e^{-i.0.x} = \int_{-\infty}^{\infty} f(x) dx$

Fourier transform of a function f(x) evaluated at k=0 is simply the integral of this function over the entire x axis. And it is always real if f(x) is real. We call this the DC Component. Thus, the DC component is the Fourier transform of a function f(x) evaluated at k=0.

5.2.5. Amplitude and phase

In general, Fourier transform of $F(k)$ is a complex function with a real and an imaginary part:

$F(k) = Re(F(k)) + i Im(F(k))$

Which also can be written as:

$F(k) = |F(k)|e^{i\theta(k)}$

where:

$|F(k)| = \sqrt(Re(F(k))^2 + (Im(F(k)))^2)$

and

$\theta(k) = arctan\frac{Im(F(k))}{Re(F(k))}$

5.2.6. properties

There are several properties of Fourier transform. We will cover each one in this section.

5.2.6.1. Linearity

If we have a two functions in the spatial domain, with each having a Fourier transform:

$f_1(x) \leftrightarrow F_1(k)$

$f_2(x) \leftrightarrow F_2(k)$

Then the summation of the two functions in the spatial domain is also the summation in the Fourier domain as such:

$a_1f_1(x) + a_2f_2(x) \leftrightarrow a_1F_1(k)+a_2F_2(k)$

5.2.6.2. Scaling

If we have a function $f(x)$ with its Fourier transform $F(k)$ : $f(x) \leftrightarrow F(k)$ , then if you stretch your $f(x)$ by 0.5, we would have an $F(k)$ with narrower band and higher amplitude. This is mathematically expressed as:

$f(ax) \leftrightarrow \frac{1}{|a|} F(\frac{k}{a})$

If a>1, we shrink the signal (the signal contains more high frequency content), then the Fourier transform of this new signal will contain higher frequencies (hence the $\frac{k}{a}$ )
If a<1, we stretch the signal (the signal contains more low-frequency contents) then the Fourier transform of this new signal will contain lower frequencies

An example of this property can be seen in the following graph.

Thus, a broad anomaly will have a narrower amplitude spectrum than a narrow anomaly. Because the width of an anomaly is directly related to the depth of its source, we can expect that the narrowness of a Fourier-transformed anomaly will also be related to the depth of the source.

Scaling Property

If you stretch in spatial domain $\leftrightarrow$ you are shrinking in Fourier domain

If you shrink in spatial domain $\leftrightarrow$ you are stretching in Fourier domain

5.2.6.3. Shifting

If you have a function such that $f(x) \leftrightarrow F(k)$ , then:

$f(x-x_0) \leftrightarrow F(k)e^{-ikx_0}$

Shifting a function along the x axis in the spatial domain is equivalent to adding a linear phase factor to the function’s Fourier transform. The amplitude spectrum is unaffected.

5.2.6.4. Differentiation

If we have a function such that: $f(x) \leftrightarrow F(x)$

then:

$\frac{d}{dx} f(x) \leftrightarrow ikF(k)$

As mentioned before in the previous section, the phase will change by $\frac{\pi}{2}$ . Thus, the amplitudes for low-frequency components will be suppressed, and the amplitudes for high-frequency components will be amplified.

The concept of differentiation in the Fourier domain is applied in potential field data. Specifically, for derivative-based anomaly enhancement. For example, the following figure shows a bouguer anomaly and Gzz data map. We can see that by taking the derivative of the gravity data, the anomalies are more enhanced and thus we see the boundaries of the anomaly clearer.

Another way to enhance the anomaly is to take the horizontal gradient of the magnitude as seen in the following picture. We can see that the boundaries are more accentuated in the data map on the right.

Bouguer anomaly

Tensor data

With how well derivatives are in enhancing the anomalies, it is no surprise that we would want to do Fourier transform to these data maps. This is because, as mentioned before, the math involved in derivation in the fourier domain is much more simpler than in spatial domain.

Summary of Fourier Transform Properties

The properties of Fourier transform are:

Linearity: Suppose that $f_1(x) \leftrightarrow F_1(k)$ and $f_2(x) \leftrightarrow F_2(k)$ , then:

$a_1f_1(x) + a_2f_2(x) \leftrightarrow a_1F_1(k)+a_2F_2(k)$

Scaling: Suppose that $f(x) \leftrightarrow F(k)$ , then:

$f(ax) \leftrightarrow \frac{1}{|a|} F(\frac{k}{a})$

Shifting: Suppose that $f(x) \leftrightarrow F(k)$ , then:

$f(x-x_0) \leftrightarrow F(k)e^{-ikx_0}$

Differentiation: Suppose that $f(x) \leftrightarrow F(x)$ , then:

$\frac{d}{dx} f(x) \leftrightarrow ikF(k)$

5.2.6.5. Fourier transform of 1/r

If you still remember the gravity potential, the equation is shown:

$U(p)=\gamma\frac{dm}r$

where $\gamma$ is the gravitational constant, $dm$ is a small mass over the $r$ .

Similar to the magnetic scalar potential:

$V(p)=\frac{\mu_0}{4\mathrm\pi}\frac{\mathbf{m}\cdot\widehat{ \mathbf{r}}}{r^2}=-\frac{\mu_0}{4\mathrm\pi}\mathbf{m}\cdot\nabla_p\frac1r$

We observe that both gravity potential and magnetic scalar potential have the $\frac{1}r$ . Before we apply Fourier transform to gravity potential and magnetic scalar potential, we need to figure out the Fourier transform of $\frac{1}r$ . Specifically, we need to make use of Bessel function and Hankel transform. But the derivation is mathematically involved, so we will skip the detailed derivations. For those who are interested, please refer to Blakely (1996, p271-273).

After a series complicated derivation, we obtain:

$\mathcal F(\frac1r)=2\mathrm\pi\frac{\mathrm e^{\left|\mathrm k\right| ({\mathrm z}_0-\mathrm z^{'})}} {\left|\mathrm k\right|}$

where, $z^'$ > $z_0$ , $\left|\mathrm k\right| \neq 0$ , $k$ is wave number.

The equation above is the results for Fourier transform of $\frac{1}r$ .

5.3. Terminology

5.3.1. Amplitude and phase, amplitude sepctrum and phase spectrum

We can use an equation to do Fourier transform (shown in below).

$\mathcal F(k)=\int_{-\infty}^{+\infty}{f(x)e^{-ikx}}\operatorname dx$

where $f(x)$ is a signal. IN most general cases, Fourier transform is a complex function with a real and an imaginary part.

$\mathcal F(k)=Re\mathcal F(k)+iIm\mathcal F(k)$

Also, if we apply Euler’s formulation, the Fourier transform can be expresses as:

$\mathcal F(k)=\left|\mathcal F(k)\right|e^{i\Theta(k)}$

where,

$\left|\mathcal F(k)\right|=\sqrt{(Re\mathcal F{(k))}^2+(Im\mathcal F{(k))}^2}$

$\Theta(k)=arc\tan\frac{Im\mathcal F{(k)}}{Re\mathcal F(k)}$

The $\left|\mathcal F(k)\right|$ is called amplitude and $\Theta(k)$ is named phase. If we want to transform a signal from spatial domain to the Fourier domain, we need to know both amplitude and phase.

$\left|\mathcal F(k)\right|$ is simply called amplitude spectrum, which is a function of wave number $k$ . In similar, we also called $\Theta(k)$ phase spectrum, which means how phase change with the wave number $k$ .

5.3.2. Energy and power spectrum

The total energy of a real function f(x). It is an integral of entire real axis of the square of amplitude, which measures how much energy in function $\mathcal F(k)$ .

$E=\int_{-\infty}^{+\infty}\left|\mathcal F(k)\right|^2\operatorname dk$

where $\left|\mathcal F(k)\right|$ is Fourier transform which is a complex number. $\left|\mathcal F(k)\right|^2$ is also called energy density function over frequency, which tells us, how energy of this function distribute with different wave numbers over frequency band. Sometimes, $\left|\mathcal F(k)\right|^2$ is named energy spectral density, power spectral density (PSD) or power spectrum (PS).

Power spectrum: the voice waveform over time (left) has a broad audio power spectrum (right). https://en.wikipedia.org/wiki/File:Voice_waveform_and_spectrum.png

We can see some high amplitude and low amplitude in the figure above (left), if we do the integration of square of amplitude, we will obtain the power spectrum (right).

5.4. Gravity foR a point mass in Fourier domain

Recall Fourier transform of 1/r

Based on the Fourier transform of 1/r, we can calculate the 3D gravity modeling in Fourier domain with point mass.

$U=\gamma\frac\mu r$

where $\mu$ is point mass, U is gravitational potential. Then we apply Fourier transform, and will obtain:

$\mathcal F(U)=\mathcal F(\gamma\frac\mu r)=\gamma\mu\mathcal F(\frac1r)$

After applying Fourier transform of 1/r, it is easily to get:

$\mathcal F(\frac1r)=2\mathrm\pi\frac{\mathrm e^{\left|\mathrm k\right| ({\mathrm z}_0-\mathrm z^{'})}} {\left|\mathrm k\right|}$

$\mathcal F(U)=2\mathrm\pi\gamma\mu\frac{\mathrm e^{\left|\mathrm k\right| ({\mathrm z}_0-\mathrm z^{'})}} {\left|\mathrm k\right|}$

where, $z^'$ > $z_0$ , $\left|\mathrm k\right| \neq 0$ .

We currently measures the vertical component of gravity field Gz, which is derivative of gravity potential in z direction.

$g_z=\frac{\partial U}{\partial z}=\gamma\mu\frac{\partial{}}{\partial z}\frac1r$

Then, we apply Fourier transform,

$\mathcal F(g_z)=\mathcal F(\gamma\mu\frac\partial{\partial z}\frac1r)=\gamma\mu\mathcal F(\frac\partial{\partial z}\frac1r)$

The derivative in z direction has nothing to do in x, y direction, so we can directly take the derivative outside, and obtain the following equation:

$\mathcal F(g_z)=\gamma\mu\frac\partial{\partial z}\mathcal F(\frac1r)=2\mathrm \pi \gamma\mu e^{\left|\mathrm k\right|({\mathrm z}_0-\mathrm z^{'})}$

Thus, we obtain the gravity measure (z component) after Fourier transform:

$\mathcal F(g_z)=2\mathrm \pi \gamma\mu e^{\left| \mathrm k\right|({z}_0- z^{'})}$

Next, we are going to make a few observations about the equation, which will help us better understand the gravity.

Power spectrum of gravity anomalies caused by a point mass of 1 kg at depth 1 km (assuming observation height 0 km).

The x axis is wave number, y axis is normalized power spectrum of gravity anomalies with a point mass. This figure thus tell us how much energy there is at each wavenumber. The most obvious thing here is exponential decay. Another observation is maximum energy occurs at k=0 (zero wavenumber). When k = 0, we can obtain:

$\mathcal F(g_z)=2\mathrm \pi \gamma\mu$

where, $\gamma$ is a constant value, $\mu$ is mass.

We can easily observe that the power spectrum value at k=0 is proportional to the total mass.

Observation from the power spectrum in log scale, we noticed that it is a straight line with negative slope. Besides, the energy decreases exponentially with increasing wavenumber. The black lines show that the x axis is 1, y axis is -1, so the slope is -1. We knew that the depth is negative slope, so the depth for the anomaly in the figure above is 1 m.

The rate of decreases on log scale is equal to the elevation difference between the source body and observation plane.

$\mathcal F(g_z)=2\mathrm \pi \gamma\mu e^{\left| \mathrm k\right|({z}_0- z^{'})}$

$\log(\mathcal F(g_z))=\log(2\mathrm{πγμ})+\left|k\right|(z_0-z^{'})$

If we set $h=z_0-z^{'}$ ,

$\log(\mathcal F(g_z))=\log(2\mathrm{πγμ})+\left|k\right|h$

The equation above mathematically illustrates that the rate of decrease depends on the depth of the source body.

The figures above have same depth but different masses, we noticed the trend is completely same, which means the mass does not change the decay rate and the decay rate has only to do with h.

Power spectrum of gravity anomalies caused by a point mass of 1 kg at depths 1 km, 2km, 6km, and 10km (observation height 0 km).

Here we created a synthetic model with same mass but located in different depth. The red line is shallowest one. We noticed from the figure above is if the source body is deeper, the energy decay faster and faster. Specifically, look at the source body located at 10 km, the energy decay to 0 at 0.2 radians/km. In comparison, if source body located at 2 km (green) line, the energy decays to 0 when wavenumber is about 1.2 radians/km.

Then, we plot 4 cases in log scale, shown in below:

In similar trend, the source body is deeper, the energy decay is faster and faster, where energy decay is represented by the slope. The slope equal to -h, so the slope is smaller, the depth is deeper and the energy decay is faster.

The summaries of two figures above

The decay rate depends upon depth.
The deeper the source, the faster the energy decays.
As source become deeper, energy at higher – wave numbers become smaller.
The gravity anomaly is approximately band limited.

In previous experiment, we fixed observation at 0 km, and varied the depth of source bodies: 1km, 2km, 6 km and 10 km. If we fix the source body at a depth of 1 km and change the observation height: 0km, -1 km, -5 km, and -9km, we will observe the exact same thing (shown in the figure below).

(top left) changing observation heights. (top right) changing source body depth. (bottom left) changing observation heights in log scale. (bottom right) changing source body depth in log scale.

From the figure above (top left and bottom left), we noticed that when source depth fixed, the decay rate depends upon observation height. The higher the observation height, the faster the energy decays. As observation become higher, energy at higher, wave numbers become smaller.

This observation plays a critical role in data acquisition that is: if you want to resolve fine details (i.e., high wave number information) of your targets, you need to be closer to your target. If using the airborne platform, we have to fly as lower as we can.

Here are several examples.

This is one field gravity data map simulated from five source bodies with fixed depth. The data map shown above is 1 meter above the surface. We can easily see the features at data maps.

In the top right figure, we keep source body shape, depth same and the only thing changed is observation height. The observation height is changed from 1 m to 100 m. Compared with previous map, we noticed that anomalies become much smooth. But we can still see a few features of source body. We can still tell there something subsurface. Keeping increasing observation height to 200 m, this map is smooth and it is always explained as one source body rather five.

The observation from the above figures is the height is higher, the gravity anomalies become smoother that means less high frequency information. This phenomenon can be explained as the higher the observation height, the faster the energy decays and wave numbers become smaller.

Another example is:

From the figure above, we noticed that the height of observation is increased, the power spectrum becomes more smooth. The more details about the figures above will be discussed in the following sections.

A quick summary

We have made quite a few observations based on this simple equation.
These insights help understand gravity well.
This is one of the reasons why we look at gravity data in Fourier doamin.

5.5. FOURIER DOMAIN MODELING OF MAGNETIC DATA DUE TO A DIPOLE

Let’s move on to the magnetic data for the simplest case of considering a magnetic dipole.

5.5.1. Fourier transform of magnetic potential

For the magnetic scalar potential $V(P)$ , which is the dot product of two vectors expressed in the following equation:

$V(P)=-\frac{{{\mu }_{0}}}{4\pi }\mathbf{m}\cdot {{\nabla }_{P}}\frac{1}{r}$

Where the dipole moment is expressed as its magnitude and directions in 3D by $\mathbf{m}=m{{[{{\hat{m}}_{x}},{{\hat{m}}_{y}},{{\hat{m}}_{z}}]}^{T}}$ .

According to the definition of the dot product, the above magnetic scalar potential expression can be rearranged into the following equation:

$V(P)=-\frac{{{\mu }_{0}}}{4\pi }m\left( {{{\hat{m}}}_{x}}\frac{\partial }{\partial x}\frac{1}{r}+{{{\hat{m}}}_{y}}\frac{\partial }{\partial y}\frac{1}{r}+{{{\hat{m}}}_{z}}\frac{\partial }{\partial z}\frac{1}{r} \right)$

If we apply the Fourier transform to the magnetic potential in the above equation, through using properties such as the linearity and differentiation, we will then have the derivations as shown below. Since we are taking the Fourier transform in 2D x-y domain, the z-direction differentiation is not relevant here.

$F\left[ V \right]=-\frac{{{\mu }_{0}}}{4\pi }m\cdot \left( {{{\hat{m}}}_{x}}i{{k}_{x}}F\left[ \frac{1}{r} \right]+{{{\hat{m}}}_{y}}i{{k}_{y}}F\left[ \frac{1}{r} \right]+{{{\hat{m}}}_{z}}\frac{\partial }{\partial z}F\left[ \frac{1}{r} \right] \right)$

After above derivation, we are left with familiar expression of $F\left[ \frac{1}{r} \right]$ , which equals $F\left( \frac{1}{r} \right)=2\pi \frac{{{e}^{\left| k \right|\left( {{z}_{0}}-{{z}^{'}} \right)}}}{\left| k \right|}$ , with z’ > z0.

Therefore, after substituting it into the Fourier transform equation, the simplified expression for the Fourier transform of magnetic potential is as following, consisting of four parts, the constant number $-\frac{{{\mu }_{0}}}{2}$ , the magnitude of dipole moment $m$ , the complex number ${{\Theta }_{m}}$ , and the exponential expression ${{e}^{\left| k \right|\left( {{z}_{0}}-{{z}^{'}} \right)}}$ which only depends on the source body depth.

$F\left[ V \right]=-\frac{{{\mu }_{0}}}{2}m{{\Theta }_{m}}{{e}^{\left| k \right|\left( {{z}_{0}}-{{z}^{'}} \right)}}$

Where the complex number

${{\Theta }_{m}}={{\hat{m}}_{z}}+i\frac{{{{\hat{m}}}_{x}}{{k}_{x}}+{{{\hat{m}}}_{y}}{{k}_{y}}}{\left| k \right|}$

depends only on the orientation of the dipole.

5.5.2. Fourier transform of total-field anomaly

Practically in the field survey, most likely we will collect total-field anomaly $\Delta T$ as the magnetic data. Therefore, we can derive the vector B field by taking the negative gradient of the potential V, which can be mathematically expressed as:

$\mathbf{B}=-{{\nabla }_{p}}V$

That is, any component (such as ${{B}_{x}},{{B}_{y}},{{B}_{z}},\Delta T$ ) can be derived from the directional derivative of the magnetic potential.

As we have talked before, the total field anomaly $\Delta T$ is simply the projection of the Bfield onto the inducing field direction, which is $\mathbf{\hat{f}}={{\left[ {{{\hat{f}}}_{x}},{{{\hat{f}}}_{y}},{{{\hat{f}}}_{z}} \right]}^{T}}$ . Therefore, mathematically, the total field anomaly is equal to the dot product of the inducing field direction with the Bfield, as shown below:

$\Delta T=-\mathbf{\hat{f}}\cdot {{\nabla }_{p}}V=-{{\hat{f}}_{x}}\frac{\partial }{\partial x}V-{{\hat{f}}_{y}}\frac{\partial }{\partial y}V-{{\hat{f}}_{z}}\frac{\partial }{\partial z}V$

Similarly, after applying the Fourier transform to the total field anomaly, we will get the expression as follows:

$\mathbf{F}\left[ \Delta T \right]=-{{{\hat{f}}}_{x}}i{{k}_{x}}\mathbf{F}\left[ V \right]-{{{\hat{f}}}_{y}}i{{k}_{y}}\mathbf{F}\left[ V \right]-{{{\hat{f}}}_{z}}\frac{\partial }{\partial z}\mathbf{F}\left[ V \right]$

After simplifying this expression, we will have the Fourier transform as shown below, which consists of four parts: the constant number, the magnitude of dipole moment $m$ (i.e., it depends on the source strength), the two directional dependent complex numbers ${{\Theta }_{m}}$ and ${{\Theta }_{f}}$ , and the depth-relevant exponential expression ${{e}^{\left| k \right|\left( {{z}_{0}}-{{z}^{'}} \right)}}$ .

$\mathbf{F}\left[ \Delta T \right]=\frac{{{\mu }_{0}}}{2}m{{\Theta }_{m}}{{\Theta }_{f}}\left| k \right|{{e}^{\left| k \right|\left( {{z}_{0}}-{{z}^{'}} \right)}}$

Where another complex number is ${{\Theta }_{f}}={{\hat{f}}_{z}}+i\frac{{{{\hat{f}}}_{x}}{{k}_{x}}+{{{\hat{f}}}_{y}}{{k}_{y}}}{\left| k \right|}$ , and it depends only on inducing field direction;

And the complex number ${{\Theta }_{m}}={{\hat{m}}_{z}}+i\frac{{{{\hat{m}}}_{x}}{{k}_{x}}+{{{\hat{m}}}_{y}}{{k}_{y}}}{\left| k \right|}$ depends only on the orientation of the dipole.

5.5.3. Fourier transform of total-field anomaly in polar coordinates

In order to better understand the Fourier transform of the total field anomaly expression derived above, let’s introduce its representation in polar system.

In the formulas of the complex number ${{\Theta }_{m}}$ ,

${{\Theta }_{m}}={{\hat{m}}_{z}}+i\frac{{{{\hat{m}}}_{x}}{{k}_{x}}+{{{\hat{m}}}_{y}}{{k}_{y}}}{\left| k \right|}$

it consists of ${{k}_{x}}$ and ${{k}_{y}}$ , which are wavenumbers in x and y directions. If we express ${{k}_{x}}$ and ${{k}_{y}}$ in polar coordinates, graphically explained in the figure below, we will have:

${{k}_{x}}=\left| k \right|\cos \phi$

${{k}_{y}}=\left| k \right|\sin \phi$

where $\left| k \right|$ is the radial wavenumber.

After replacing the ${{k}_{x}}$ and ${{k}_{y}}$ in polar coordinates, the ${{\Theta }_{m}}$ can then be written as:

${{\Theta }_{m}}={{\hat{m}}_{z}}+i\left( {{{\hat{m}}}_{x}}\cos \phi +{{{\hat{m}}}_{y}}\sin \phi \right)$

Another way of representing ${{\Theta }_{m}}$ is by using the Euler’s formula through amplitude and phase as a function of the angle $\phi$ , which can be re-written as

${{\Theta }_{m}}={{A}_{m}}\left( \phi \right){{e}^{i{{P}_{m}}\left( \phi \right)}}$

Similarly, we can follow the same procedure to represent ${{\Theta }_{f}}$ , which we will get:

${{\Theta }_{f}}={{A}_{f}}\left( \phi \right){{e}^{i{{P}_{f}}\left( \phi \right)}}$

Now, let us substitute the polar representations of ${{\Theta }_{m}}$ and ${{\Theta }_{f}}$ into the Fourier transform expression of the total field anomaly, we will get the following expression:

$\mathbf{F}\left[ \Delta T \right]=\frac{{{\mu }_{0}}}{2}m{{A}_{m}}\left( \phi \right){{e}^{i{{P}_{m}}\left( \phi \right)}} {{A}_{f}}\left( \phi \right){{e}^{i{{P}_{f}}\left( \phi \right)}}\left| k \right|{{e}^{\left| k \right|\left( {{z}_{0}}-{{z}^{'}} \right)}}$

5.5.3.1. Interpretations of the first part of the amplitude spectrum

If we only focus on the amplitude spectrum term in the above expression, which is a function of wavenumber $k$ and angle $\phi$ as shown below:

$\mathbf{F}\left[ \Delta T \right]=\frac{{{\mu }_{0}}}{2}m{{A}_{m}}\left( \phi \right){{A}_{f}}\left( \phi \right)\left| k \right|{{e}^{\left| k \right|\left( {{z}_{0}}-{{z}^{'}} \right)}}$

The first part of this amplitude spectrum $\frac{{{\mu }_{0}}}{2}m{{A}_{m}}\left( \phi \right){{A}_{f}}\left( \phi \right)$ only varies with the angle $\phi$ , since once we fix the strength of the source body the magnitude $m$ will become a constant. Therefore, given a fixed angle $\phi$ , along any ray originating from the origin, this part of the amplitude spectrum value is constant, that is, there is no decay on the ray since all values on this ray are equal.

Thus, we can express the first part as $C(\phi)$ :

$C\left( \phi \right)=\frac{{{\mu }_{0}}}{2}m{{A}_{m}}\left( \phi \right){{A}_{f}}\left( \phi \right)$

In brief summary so far, the amplitude spectrumterm can be written as:

$A\left( k,\phi \right)=C\left( \phi \right)\left| k \right|{{e}^{\left| k \right|\left( {{z}_{0}}-{{z}^{'}} \right)}}$

5.5.3.2. Interpretations of the second part of the amplitude spectrum

Now let us focus on the second part of the amplitude spectrumterm $\left| k \right|{{e}^{\left| k \right|\left( {{z}_{0}}-{{z}^{'}} \right)}}$ , which has nothing to do with angle $\phi$ , but it is related with the observation height ${{z}_{0}}$ and the depth ${{z}^{'}}$ , and the radial wavenumber ${\left| k \right|$ . Then how does it look like?

The graph plotted below is what it looks like as a function of the radial wavenumber ${\left| k \right|$ . This figure implies that the energy increases to a peak value along with the increase of the wavenumber before starting decreasing.

Therefore, this second part of the amplitude spectrum determines how the amplitude spectrum changes along a ray and determines the shape of the energy change in that ray.

Moreover, since this part has nothing to do with angle $\phi$ , the shape of the amplitude spectrum along any ray is the same; and this shape is determined by the depth of the source body ${{z}^{'}}$ .

Therefore, by looking at the shape of the amplitude spectrum along any ray, we can infer depth of the source body, no matter which ray we are looking at, since they all have the exact same shapes!

In addition, the maximum energy does not occur at wavenumber k=0, but at a larger radial wavenumber that is dependent on the depth of the source dipole; once the wavenumbers exceeding this peak energy wavenumber, energy starts decaying monotonically and approaching exponential decay at higher wavenumbers.

5.6. Upward continuation

5.6.1. Concepts of upward continuation

Upward continuation is the process of using originally observed data to calculate the gravity or magnetic anomaly response that would be observed at locations above the original observation surface. Shifting the observation location to a higher elevation removes high-frequency signals from the near surface due to signal attenuation with distance. This allows for the focus of lower-frequency signals from deeper parts of a survey.

5.6.2. Upward continuation in spatial domain

Recall Green’s third identity:

$V\left( P \right)=\frac{1}{4\pi }\iint\limits_{Surf}{\left( \frac{1}{r}\frac{\partial V}{\partial \mathbf{\hat{n}}}-V\frac{\partial }{\partial \mathbf{\hat{n}}}\frac{1}{r} \right)}dS$

We previously stated that this is the theoretical basis for the upward continuation and the equivalent source technique. If we consider $V$ to be harmonic, Green’s third identity is a representation formula which a potential field can be calculated at any point simply from the behavior of the field at its boundaries. This also means that no knowledge of the sources is required, except that none are located within the region. Unfortunately, the equation above requires not only the values of $V$ on the surface, but also the values of its vertical derivative. These values are unlikely to be available in most practical applications. Fortunately, through some mathematical manipulation we can eliminate the derivative from the equation, but it is fairly complex so we won’t show that here.

Eventually, we will obtain the following equation:

$V(x,y,z_0 - \Delta z) = \frac{\Delta z}{2\pi}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\frac{V(x',y',z_0)}{[(x-x')^2 + (y-y')^2 + \Delta z^2]^{3/2}}dx'dy'$

$\Delta z > 0$

This equation is the upward-continuation integral. It shows how to calculate the value of a potential field at any point above a level, horizontal surface at $z_0$ from a complete knowledge of the field on the surface. This means upward continuation can be done in the spatial domain using this integral.

At this point, we will need to cover convolution. The convolution of two functions $f(x)$ and $g(x)$ is:

$h(x) = \int\limits_{-\infty}^{\infty} f(x')g(x-x')dx'$

We can also convert this equation to the Fourier domain if we assume that $h(x) \leftrightarrow H(k), f(x) \leftrightarrow F(k)$ and $g(x) \leftrightarrow G(k)$ :

$H(k) = F(k)G(k)$

We will also need to examine convolution in 2D for upward continuation. In 2D, the convolution of two functions f(x) and g(x) in the spatial domain becomes:

$h(x,y) = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} f(x',y')g(x-x',y-y')dx'dy'$

If we assume that $h(x,y) \leftrightarrow H(k_x,k_y), f(x,y) \leftrightarrow F(k_x,k_y)$ and $g(x,y) \leftrightarrow G(k_x,k_y)$ , the Fourier domain equation is:

$H(k_x,k_y) = F(k_x,k_y)G(k_x,k_y)$

This is important as we go back to our previously established upward continuation equation:

$V(x,y,z_0 - \Delta z) = \frac{\Delta z}{2\pi}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\frac{V(x',y',z_0)}{[(x-x')^2 + (y-y')^2 + \Delta z^2]^{3/2}}dx'dy'$

$\Delta z > 0$

We can define

$\psi(x,y,\Delta z) = \frac{\Delta z}{2\pi}\frac{1}{[x^2 + y^2 + \Delta z^2]^{3/2}}$

to obtain the following equation:

$V(x,y,z_0 - \Delta z) = \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}V(x',y',z_0)\psi(x-x',y-y',\Delta z)dx'dy'$

This equation is 2D convolution!

5.6.3. upward continuation in fourier domain

We can find the Fourier representation of upward continuation using the following equation:

$F[V_U] = F[V]F[\psi]$

All we need now is an analytical expression of F[\psi], and we will skip the derivations here to get the following equation:

$F[\psi] = e^{-\Delta z |k|}$

$\Delta z > 0$

At this point, we can find the Fourier transform of the upward-continued field in three steps:

Transform the observed data via FFT (Fast Fourier Transform)
Multiply by the exponential operator (e^{-h\sqrt{w_x^2 + w_y^2}})
Apply the inverse transform to obtain the upward-continued field

5.6.4. Example

In the example below, we will be showing the influence that upward continuation has on clean and noisy data. The figure below shows the difference between the shape of the input data in 2D and 3D between its original elevation and its upward continuation.

The following figure shows the same figures after noise has been introduced to the data. It is clear that the noise has a significant influence on the input data at its original data, but has been largely reduced by the upward continuation.

The figure below shows the magnetic response of the input data at various elevations. Pay attention to the smoothing of the overall feature along with the reduction in high-frequency data as the elevation increases.

The same figure was reproduced using noisy data. Notice the significant difference between the accurate and noisy data at the lower elevations, and the reduction in noise as the responses are recorded at higher elevations of upward continuation.

5.7. Reduction to pole

5.7.1. Fundamentals of reduction to pole

Reduction to Pole (RTP) is the process of adjusting recorded magnetic data so that it appears as if it was recorded if the Earth’s inducing magnetic field was vertical. The figure below shows this adjustment, and note the difference in terms of symmetry between the unadjusted and adjusted data. It is important to note that RTP is difficult to do at low magnetic inclinations (next to the equator).

The equation for RTP in the Fourier domain can be listed as the following:

$F[\Delta T_r] = F[\Delta T]F[\psi_r]$

$F[\psi] = \frac{1}{\Theta_m\Theta_f} = \frac{|k|^2}{(\vec{R}\cdot\hat{f})(\vec{R}\cdot\hat{m})}$

$\vec{R} = (ik_x,ik_y,|k|)$

where $\hat{f}$ is the inducing field direction and $\hat{m}$ is the magnetization direction

We can find the Fourier transform of RTP using the following three steps:

Transform the observed data via FFT (Fast Fourier Transform)
Multiply by $\frac{|k|^2}{(\vec{R}\cdot\hat{f})(\vec{R}\cdot\hat{m})}$
Apply the inverse transform to obtain RTP

5.7.2. Examples of reduction to pole

https://gpg.geosci.xyz/content/magnetics/magnetics_processing.html

https://em.geosci.xyz/content/case_histories/balboa/processing.html

5.8. Issues with fourier modelling

Throughout this chapter, we have gone through different methods in Fourier domain modelling using FFT (Fast Fourier Transform). However, FFT-based processing methods have two major prerequisites on data:

The observation surface needs to be planar.
Data needs to be interpolated to uniform grid.

These two prerequisites/assumptions is not made in practice because:

Data are often acquired on undulating surface. For example, in an airborne survey with a helicopter, the flight height will be made constant to 50 m above the underlying terrain. However, the underlying terrain has a lot of hills and troughs. As such, the data measured by the helicopter is not from a planar surface because it follows the topography.
Majority of data are acquired along flight lines or scattered stations. What this means is that the survey grid is not uniform in practice. For example, suppose an airborne survey using a helicopter. The helicopter will fly in one direction and obtain some measurements. After going through one survey line, it will move to the next survey line, which is around 100m away from the first survey line, and then it measures the data along that survey line. Now, the problem with this is that the data measured along the survey line (let’s say y direction, measurement every 5m) is much finer than data measured adjacent to the survey lines (let’s say x direction, measurement every 100m). Thus, you can have a grid of 5m x 100m. Clearly this grid is not uniform, and thus FFT will not work with this survey grid. One way to solve this problem is to use interpolation method but this has other problems as well
The process of interpolation is a heavy processing step. Interpolation can be a dangerous process because it can create artifacts that does not exist in the data. Processed data maps would have many artifacts that would lead to misinterpretation.

In the next chapter, we will be going through a different method that does not produce artifacts in the data. This method is called the equivalent source method.

License

Icon for the Creative Commons Attribution-NonCommercial 4.0 International License

5.1. background

5.1.1. Motivation

5.1.2. frequency and phase

5.1.2.1. Waveforms

5.1.2.2. Frequency

5.1.2.3. Phase

5.1.2.4 General notation

5.1.3. complex variable

5.1.3.1 Absolute value and argument

5.1.3.2. Euler’s formula

5.1.3.3. Deriving a complex number

5.1.3.4. Why bother with complex numbers?

5.2. Fourier transform

5.2.1. waveforms and amplitude spectrum

5.2.2. definition

5.2.3. from time domain to spatial domain

5.2.4. DC-component

5.2.5. Amplitude and phase

5.2.6. properties

5.2.6.1. Linearity

5.2.6.2. Scaling

5.2.6.3. Shifting

5.2.6.4. Differentiation

5.2.6.5. Fourier transform of 1/r

5.3. Terminology

5.3.1. Amplitude and phase, amplitude sepctrum and phase spectrum

5.3.2. Energy and power spectrum

5.4. Gravity foR a point mass in Fourier domain

5.5. FOURIER DOMAIN MODELING OF MAGNETIC DATA DUE TO A DIPOLE

5.5.1. Fourier transform of magnetic potential

5.5.2. Fourier transform of total-field anomaly

5.5.3. Fourier transform of total-field anomaly in polar coordinates

5.5.3.1. Interpretations of the first part of the amplitude spectrum

5.5.3.2. Interpretations of the second part of the amplitude spectrum

5.6. Upward continuation

5.6.1. Concepts of upward continuation

5.6.2. Upward continuation in spatial domain

5.6.3. upward continuation in fourier domain

5.6.4. Example

5.7. Reduction to pole

5.7.1. Fundamentals of reduction to pole

5.7.2. Examples of reduction to pole

5.8. Issues with fourier modelling

License

Share This Book