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
datadir = '../src/examples/'
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.createSurvey(datadir + 'cover-crop/coverCrop.csv')
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()
k.createTimeLapseSurvey([datadir + 'timelapse-wheat/170316.csv',
datadir + 'timelapse-wheat/170427.csv',
datadir + 'timelapse-wheat/170516.csv'])
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)
<ipython-input-1-5c138adf3fe3> in <module>
2 k.createTimeLapseSurvey([datadir + 'timelapse-wheat/170316.csv',
3 datadir + 'timelapse-wheat/170427.csv',
----> 4 datadir + 'timelapse-wheat/170516.csv'])
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:
--> 172 self.readFile(fname, targetProjection=targetProjection, unit=unit)
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
--> 605 return _read(filepath_or_buffer, kwds)
606
607
/usr/local/lib/python3.7/site-packages/pandas/io/parsers.py in _read(filepath_or_buffer, kwds)
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()
k.createTimeLapseSurvey([datadir + 'timelapse-wheat/170316.csv',
datadir + 'timelapse-wheat/170427.csv',
datadir + 'timelapse-wheat/170516.csv'])
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: 