Equivalent Sources

Equivalent Sources

Most potential field surveys gather data along scattered and uneven flight lines or ground measurements. For a great number of applications we may need to interpolate these data points onto a regular grid at a constant altitude. Upward-continuation is also a routine task for smoothing, noise attenuation, source separation, etc.

Both tasks can be done simultaneously through an equivalent sources [Dampney1969] (a.k.a equivalent layer). The equivalent sources technique consists in defining a finite set of geometric bodies (like point sources) beneath the observation points and adjust their coefficients so they generate the same measured field on the observation points. These fitted sources can be then used to predict the values of this field on any unobserved location, like a regular grid, points at different heights, any set of scattered points or even along a profile.

The equivalent sources have two major advantages over any general purpose interpolator:

  • it takes into account the 3D nature of the potential fields being measured by considering the observation heights, and

  • its predictions are always harmonic.

Its main drawback is the increased computational load it takes to fit the sources’ coefficients, both in terms of memory and computation time).

Harmonica has a few different classes for applying the equivalent sources techniques. Here we will explore how we can use the harmonica.EquivalentSources to interpolate some gravity disturbance scattered points on a regular grid with a small upward continuation.

We can start by downloading some sample gravity data over the Bushveld Igneous Complex in Southern Africa:

import ensaio
import pandas as pd

fname = ensaio.fetch_bushveld_gravity(version=1)
data = pd.read_csv(fname)
data
longitude latitude height_sea_level_m height_geometric_m gravity_mgal gravity_disturbance_mgal gravity_bouguer_mgal
0 25.01500 -26.26334 1230.2 1257.474535 978681.38 25.081592 -113.259165
1 25.01932 -26.38713 1297.0 1324.574150 978669.02 24.538158 -122.662101
2 25.02499 -26.39667 1304.8 1332.401322 978669.28 26.526960 -121.339321
3 25.04500 -26.07668 1165.2 1192.107148 978681.08 17.954814 -113.817543
4 25.07668 -26.35001 1262.5 1289.971792 978665.19 12.700307 -130.460126
... ... ... ... ... ... ... ...
3872 31.51500 -23.86333 300.5 312.710241 978776.85 -4.783965 -39.543608
3873 31.52499 -23.30000 280.7 292.686630 978798.55 48.012766 16.602026
3874 31.54832 -23.19333 245.7 257.592670 978803.55 49.161771 22.456674
3875 31.57333 -23.84833 226.8 239.199065 978808.44 5.116904 -20.419870
3876 31.37500 -23.00000 285.6 297.165672 978734.77 5.186926 -25.922627

3877 rows × 7 columns

The harmonica.EquivalentSources class works exclusively in Cartesian coordinates, so we need to project these gravity observations:

import pyproj

projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.values.mean())
easting, northing = projection(data.longitude.values, data.latitude.values)

Now we can initialize the harmonica.EquivalentSources class.

import harmonica as hm

equivalent_sources = hm.EquivalentSources(depth=10e3, damping=10)
equivalent_sources
EquivalentSources(damping=10, depth=10000.0)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

By default, it places the sources one beneath each data point at a relative depth from the elevation of the data point following [Cooper2000]. This relative depth can be set through the depth argument. Deepest sources generate smoother predictions (underfitting), while shallow ones tend to overfit the data.

Note

If instead we want to place every source at a constant depth, we can change it by passing depth_type="constant". In that case, the depth argument will be the exact depth at which the sources will be located.

The damping parameter is used to smooth the coefficients of the sources and stabilize the least square problem. A higher damping will create smoother predictions, while a lower one could overfit the data and create artifacts.

Now we can estimate the source coefficients through the harmonica.EquivalentSources.fit method against the observed gravity disturbance.

coordinates = (easting, northing, data.height_geometric_m)
equivalent_sources.fit(coordinates, data.gravity_disturbance_mgal)
EquivalentSources(damping=10, depth=10000.0)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Once the fitting process finishes, we can predict the values of the field on any set of points using the harmonica.EquivalentSources.predict method. For example, lets predict on the same observation points to check if the sources are able to reproduce the observed field.

disturbance = equivalent_sources.predict(coordinates)

And plot it:

import verde as vd
import matplotlib.pyplot as plt

# Get max absolute value for the observed gravity disturbance
maxabs = vd.maxabs(data.gravity_disturbance_mgal)

# Create a figure with two axes
fig, (ax1, ax2) = plt.subplots(figsize=(10, 12), nrows=2, ncols=1, sharex=True)

# Plot observed and predicted fields
tmp = ax1.scatter(
    easting,
    northing,
    c=data.gravity_disturbance_mgal,
    s=2,
    vmin=-maxabs,
    vmax=maxabs,
    cmap="seismic",
)
ax1.set_title("Observed gravity disturbance")
ax2.scatter(
    easting,
    northing,
    c=disturbance,
    s=2,
    vmin=-maxabs,
    vmax=maxabs,
    cmap="seismic",
)
ax2.set_title("Predicted gravity disturbance")

for ax in (ax1, ax2):
    ax.set_aspect("equal")
    plt.colorbar(tmp, ax=ax, label="mGal", pad=0.05, aspect=20, shrink=0.9)
plt.show()
../../_images/index_5_0.png

We can also grid and upper continue the field by predicting its values on a regular grid at a constant height higher than the observations. To do so we can use the verde.grid_coordinates function to create the coordinates of the grid and then use the harmonica.EquivalentSources.grid method.

First, lets get the maximum height of the observations:

data.height_geometric_m.max()
2167.498378036602

Then create the grid coordinates at a constant height of and a spacing of 2km; and use the equivalent sources to generate a gravity disturbance grid.

# Get region of the observations
region = vd.get_region(coordinates)

# Build the grid coordinates
grid_coords = vd.grid_coordinates(region=region, spacing=2e3, extra_coords=2.2e3)

# Grid the gravity disturbances
grid = equivalent_sources.grid(grid_coords, data_names=["gravity_disturbance"])
grid
<xarray.Dataset>
Dimensions:              (northing: 223, easting: 354)
Coordinates:
  * easting              (easting) float64 2.525e+06 2.527e+06 ... 3.231e+06
  * northing             (northing) float64 -2.816e+06 -2.814e+06 ... -2.372e+06
    upward               (northing, easting) float64 2.2e+03 2.2e+03 ... 2.2e+03
Data variables:
    gravity_disturbance  (northing, easting) float64 20.71 21.06 ... 12.31 12.08
Attributes:
    metadata:  Generated by EquivalentSources(damping=10, depth=10000.0)

And plot it

fig = plt.figure(figsize=(12, 8))

grid.gravity_disturbance.plot()
plt.gca().set_aspect("equal")
plt.show()
../../_images/index_8_0.png