# Smoothing parameters and time-lapse inversion¶

In this notebook, we will see what it the impact of the different smoothing parameters and how we can use them to perform time-lapse inversion.

import matplotlib.pyplot as plt
import numpy as np
import sys
sys.path.append('../src')
from emagpy import Problem


## Smoothing parameters¶

There are three different smoothing parameters: - alpha: linked to the vertical smoothing - beta: linked to lateral smoothing - gamma: linked to smoothing with the first inverted survey (usefull for time-lapse inversion)

k = Problem()
k.surveys[0].df = k.surveys[0].df[:30] # truncate dataset
k.show()

Removing 1 NaN from survey


### Vertical smoothing (alpha)¶

The alpha parameter regulates the vertical smoothing of the profile. It is needed to regularize the inversion. It plays an important role in smooth inversion with many layers such as used with Problem.invertGN(). To determine the optimal value of alpha a common method is to plot the so-called ‘L-curve’ that show the size of the model and data misfit as a function of alpha. It can be plotted using the Problem.lcurve() method.

k.lcurve()


From the L-curve, we look for the alpha value that offers the best trade-off between the model and data misfit. This value is usually taken at the maximum curvature point of the curve. In this case it would corresponds to 0.07. This value can then be used in the call to Problem.invert(alpha=alpha).

The graphs below shows the impact of different values of alpha. Larger alpha values will give a smoother model, while smaller alpha values can show very high/low EC values.

fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))
for i, alpha in enumerate([0, 0.07, 0.5]):
ax = axs[i]
k.invert(alpha=alpha)
k.showResults(ax=ax)
ax.set_title('alpha = {:.2f}'.format(alpha))

30/30 inverted


### Lateral smoothing (beta)¶

We can see that the data present some similarity between continguous measurements. To force lateral smoothing we can use the beta parameters. It will penalize the objective function with the difference between two continuous profiles. Note that this can be dangerous in some conditions, for instance if the first profile is badly inverted.

The figure below show the effect of lateral smoothing.

fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))
for i, beta in enumerate([0, 0.1, 0.5]):
ax = axs[i]
k.invert(beta=beta)
k.showResults(ax=ax)
ax.set_title('beta = {:.2f}'.format(beta))

30/30 inverted


### Reference survey smoothing (gamma) or time-lapse inversion¶

gamma controls the smoothing between the first survey and subsequent survey. This way a quasi time-lapse inversion can be performed.

k = Problem()
for s in k.surveys:
s.df = s.df[:20]

# show apparent values collected
fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))
for i in range(len(k.surveys)):
ax = axs[i]
k.show(index=i, ax=ax)
ax.set_ylim([12, 42])
k.setInit(depths0=[0.7], fixedDepths=[False])
bounds = [(0.1, 1),(2, 60),(2, 60)]

---------------------------------------------------------------------------

FileNotFoundError                         Traceback (most recent call last)

5 for s in k.surveys:
6     s.df = s.df[:20]

/builds/hkex/emagpy/src/emagpy/Problem.py in createTimeLapseSurvey(self, fnames, targetProjection, unit)
162             targetProjection = self.projection
163         for fname in fnames:
--> 164             self.createSurvey(fname, targetProjection=targetProjection, unit=unit)
165
166

/builds/hkex/emagpy/src/emagpy/Problem.py in createSurvey(self, fname, freq, hx, targetProjection, unit)
107             targetProjection = self.projection
108
--> 109         survey = Survey(fname, freq=freq, hx=hx, targetProjection=targetProjection, unit=unit)
110
111         # remove NaN from survey

/builds/hkex/emagpy/src/emagpy/Survey.py in __init__(self, fname, freq, hx, targetProjection, unit)
170         self.projection = None # store the project
171         if fname is not None:
173             if freq is not None:
174                 self.freqs = np.ones(len(self.coils))*freq

/builds/hkex/emagpy/src/emagpy/Survey.py in readFile(self, fname, sensor, targetProjection, unit)
201         if fname.find('.DAT')!=-1:
202             delimiter = '\t'
--> 203         df = pd.read_csv(fname, delimiter=delimiter)
204         self.readDF(df, name, sensor, targetProjection, unit)
205

/usr/local/lib/python3.7/site-packages/pandas/io/parsers.py in read_csv(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, squeeze, prefix, mangle_dupe_cols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, dialect, error_bad_lines, warn_bad_lines, delim_whitespace, low_memory, memory_map, float_precision, storage_options)
603     kwds.update(kwds_defaults)
604
606
607

455
456     # Create the parser.
--> 457     parser = TextFileReader(filepath_or_buffer, **kwds)
458
459     if chunksize or iterator:

/usr/local/lib/python3.7/site-packages/pandas/io/parsers.py in __init__(self, f, engine, **kwds)
812             self.options["has_index_names"] = kwds["has_index_names"]
813
--> 814         self._engine = self._make_engine(self.engine)
815
816     def close(self):

/usr/local/lib/python3.7/site-packages/pandas/io/parsers.py in _make_engine(self, engine)
1043             )
1044         # error: Too many arguments for "ParserBase"
-> 1045         return mapping[engine](self.f, **self.options)  # type: ignore[call-arg]
1046
1047     def _failover_to_python(self):

/usr/local/lib/python3.7/site-packages/pandas/io/parsers.py in __init__(self, src, **kwds)
1860
1861         # open handles
-> 1862         self._open_handles(src, kwds)
1863         assert self.handles is not None
1864         for key in ("storage_options", "encoding", "memory_map", "compression"):

/usr/local/lib/python3.7/site-packages/pandas/io/parsers.py in _open_handles(self, src, kwds)
1361             compression=kwds.get("compression", None),
1362             memory_map=kwds.get("memory_map", False),
-> 1363             storage_options=kwds.get("storage_options", None),
1364         )
1365

/usr/local/lib/python3.7/site-packages/pandas/io/common.py in get_handle(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)
645                 encoding=ioargs.encoding,
646                 errors=errors,
--> 647                 newline="",
648             )
649         else:

FileNotFoundError: [Errno 2] No such file or directory: '../src/examples/timelapse-wheat/170316.csv'

# inversion with gamma = 0 (no smoothing between surveys)
k.invert(gamma=0, bnds=bounds)
fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))
for i in range(3):
ax = axs[i]
k.showResults(index=i, ax=ax, maxDepth=1.5)

bounds = [(0.1, 1), (2, 60), (2, 60)]
Survey 1/3
20/20 inverted
Survey 2/3
20/20 inverted
Survey 3/3
20/20 inverted

# inversion with constrain to the first survey
k.invert(gamma=0.1, bnds=bounds)
fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))
for i in range(3):
ax = axs[i]
k.showResults(index=i, ax=ax, maxDepth=1.5)

bounds = [(0.1, 1), (2, 60), (2, 60)]
Survey 1/3
20/20 inverted
Survey 2/3
20/20 inverted
Survey 3/3
20/20 inverted


### Observations¶

In the second inversion with between survey constrain we can clearly see that the depths obtained during the first inversion constrain the depth of the other survey. We can also notice that the last two profiles are badly inverted in the first survey and so this propagate for the other subsequent surveys.

## Inverting change in ECa¶

An other way of producing time-lapse inversion is to actually invert for the change in ECa. This can be done using the CS function and the Problem.invertGN() function. See appendix of Whalley et al. (2017) for more details.

k = Problem()
for s in k.surveys:
s.df = s.df[:20]

k.computeApparentChange()
# show change in apparent values collected
fig, axs = plt.subplots(1, 3, sharex=True, figsize=(12,3))
for i in range(len(k.surveys)):
ax = axs[i]
k.show(index=i, ax=ax)
if i > 0:
ax.set_ylabel('Change in EC [mS/m]')
fig.tight_layout()

Trimming surveys and only keep common positions
Matching positions between surveys for time-lapse inversion...20 in common...done in 0.011s
Computing relative ECa compared to background (1st survey).

k.setInit(depths0=np.linspace(0.2, 2, 10)) # use a smooth model
k.invertGN(alpha=0.1)
fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))
for i in range(3):
ax = axs[i]
if i == 0: # reference survey
k.showResults(index=i, ax=ax)
else:
k.showResults(index=i, ax=ax, cmap='bwr_r', vmin=-8, vmax=2)
fig.axes[-1].set_ylabel('EC change [mS/m]')

Survey 1/3
20/20 inverted
Survey 2/3
20/20 inverted
Survey 3/3
20/20 inverted


Download python script: nb_smoothing.py

Download Jupyter notebook: nb_smoothing.ipynb

View the notebook in the Jupyter nbviewer

Run this example interactively: