STFINV library: seek source wavelet correction filter
Information for end-users
Contents of this page

The library libstfinv for end-users

Introduction

The purpose of this library is to provide methods for the derivation of source wavelet correction filters in approaches to full waveform inversion. Given a set of recorded data and a set of synthetic data (typically, but not necessarilly the impulse response of the subsurface) a source wavelet correction filter is obtained due to some optimization citerion. The synthetic waveforms are convolved with this filter wavelet and the convolved synthetics as well as the wavelet itself are returned to the user.

The source wavelet correction filter in this context not necessarily is the actual force time history of the source used in the experiment or a similar quantity of physical meaning. The source wavelet correction filter simply is the wavelet which minimizes the misfit between synthetic and recorded waveforms due to some misfit condition, if the synthetics are concolved with this wavelet. In particular this implies that the synthetics not necessarily need be the impulse response (Greens function) of the subsurface, they may simply be synthetic waveforms computed for some generic source wavelet (like a Ricker wavelet). The derived source wavelet correction filter then has to be understood with respect to this generic wavelet.

How to construct parameter strings

A specific engine is selected by passing a parameter string to the library interface like function initstfinvengine() or the constructor of the interface class stfinv::STFEngine::STFEngine(). This parameter string may further contain parameters to control the execution mode of the engine.

The parameter string starts with an ID-sequence identifying the desired engine. This ID-sequences are stored in stfinv::STFEngineFDLeastSquares::ID for the Fourier domain least squares engine. In the parameter string the ID-sequence is terminated by a colon (:).

After selecting the desired engine, the interface function strips of the ID-sequence as well as the colon from the parameter string and initializes the engine, passing the references to user workspace as well as the rest of the parameter string. The rest of the parameter string may consist of several control parameters being separated by colons (:). Each control parameter may just be a flag (switch to turn an option on) or may come along with a parameter value. The value of the parameter is separated by an equal sign (=).

Examples
  • To select frequency domain least squares (stfinv::STFEngineFDLeastSquares) and shift the returned source correction wavelet by 0.4s and switch on verbose mode, pass the following parameter string:
    fdlsq:tshift=0.4:verbose 
  • To select stfinv::STFEngineIdentity and to switch on debug level 4:
    ident:DEBUG=4 
  • To select stfinv::STFEngineFDLeastSquares, apply offset dependent weights and use a power of two to speed up the FFT:
    fdlsq:pow2:exp=1.4 

End-user usage information text files

Summary and detailed description for overall library

# this is <usage/stfinv_summary_usage.txt>
# ----------------------------------------------------------------------------
# 
STFINV -- create and apply a source wavelet correction filter
=============================================================

Purpose
-------
Create and apply a source wavelet correction filter for full-waveform
inversion of field data.
# ----- END OF usage/stfinv_summary_usage.txt ----- 
# this is <usage/stfinv_description_usage.txt>
# ----------------------------------------------------------------------------
#
STFINV -- create and apply a source wavelet correction filter
=============================================================

Different methods are provided for the derivation of source wavelet correction
filters in approaches to full waveform inversion. Given a set of recorded data
and a set of synthetic data (typically, but not necessarilly the expected
impulse response of the subsurface) a source wavelet correction filter is
obtained by application of a user-selectable optimization citerion. The
synthetic waveforms are convolved with this filter wavelet and the convolved
synthetics as well as the wavelet itself are returned to the user.

The effective time history of the seismic source used in field recordings is
not well known in most cases. This applies in particular to transient sources
(like explosives or hammer blows). The so-called 'source-time-function' might
even vary from shot to shot. For this reason it is not possible to use an
appropriate source-time-function in the initial simulation of synthetic data
in an approach of full-waveform inversion. However, after synthetic data have
been calculated using a generic source-time-function, a correction filter can
be constructed such that an improved source-time-function will reduce the
misfit to the recorded data.
# ----- END OF usage/stfinv_description_usage.txt ----- 

Summary and detailed description in handle class

# this is <stfinvany_summary_usage.txt>
# ----------------------------------------------------------------------------
#
How to control the procedure
----------------------------
Procedures can be selected by passing the procedures ID together with options
an parameters controlling the algorithm used.

Syntax of parameter string:
  ID:option:option=value:option=value
# ----- END OF stfinvany_summary_usage.txt ----- 
# this is <stfinvany_description_usage.txt>
# ----------------------------------------------------------------------------
# This file should address the selection of the different procedures. The
# overall concept of the library is explained in the base class.
How to control the procedure
----------------------------
A specific procedure is selected by passing a parameter string. This string
may further contain parameters to control the way the procedure will be
applied.

The parameter string starts with an ID-sequence identifying the desired
procedure. See the list available procedures for ID strings. In the parameter
string the ID-sequence is terminated by a colon (:).

After selecting the desired procedure, the interface function strips off the
ID-sequence as well as the first colon from the parameter string. The
remainder may consist of several control parameters being separated by colons
(:). Each control parameter may just be a flag (switch to turn on an option)
or may come along with a parameter value. The value is separated from the
parameter by an equal sign (=).

Where several values in an argument to a parameter must be separated (like in
the 'irtap' option of the Fourier domain procedures) white space ( ), commas
(,), and semicolons (;) are allowed as field delimiters, at your convenience.

Examples:
- To select Fourier domain least squares and shift the returned source
  correction filter wavelet by 0.4s and switch on verbose mode, pass the
  following parameter string:
    fdlsq:tshift=0.4:verbose
- To select the identity procedure and switch on debug level 4:
    ident:DEBUG=4
- To select Fourier domain least squares, apply offset dependent weights and
  use a power of two to speed up the FFT:
    fdlsq:pow2:exp=1.4
#
# ----- END OF stfinvany_description_usage.txt ----- 

Summary and detailed description in base class

# this is <stfinvbase_summary_usage.txt>
# ----------------------------------------------------------------------------
#
Common control parameters
-------------------------
Options and parameters in common for all procedures:
  verbose       produce verbose output (if implemented)
  DEBUG=l       produce debug output with level l
  exp=k         apply offset dependent weights to signals
# ----- END OF stfinvbase_summary_usage.txt ----- 
# this is <stfinvbase_description_usage.txt>
# ----------------------------------------------------------------------------
# Common control parameters
# -------------------------
Options and parameters in common for all procedures:
  verbose       produce verbose output (if implemented)
  DEBUG=l       produce debug output with level l
  exp=k         apply offset dependent weights to signals

Due to the amplitude decay with offset to the source signals from receivers at
larger offset would contribute less to the optimization criterion for which
the source wavelet correction filter is constructed. The option exp provides
means to add a weight factor ((r/1m)**k) to each signal, where r is the
receiver to source offset. This is used to compensate the decrease in signal
amplitude. If the energy in the original signal decays with ((-r/1m)**(2*k))
all traces will contribute at a similar level to the derived correction filter
after application of the weight factors.
# ----- END OF stfinvbase_description_usage.txt ----- 

Summary and detailed description in Fourier class

# this is <stfinvfourier_summary_usage.txt>
# ----------------------------------------------------------------------------
#
Procedures in the Fourier domain
--------------------------------
Options and parameters in common for procedures in the Fourier domain:
  fpow2              use power of two for number of coefficients
  fdiv=d             use integer multiple of d for number of coefficients
  fpad=f             padding factor
  tshift=d           delay source correction filter wavelet by d (in seconds)
                     this is applied only to the waveform of the impulse
                     response of the correction filter after application to
                     synthetic data (has no effect on convolved data)
  irtap=t1,t2,t3,t4  taper impulse response of correction filter
                     this is applied to the correction filter before
                     application to synthetic data (affects convolved data)
# ----- END OF stfinvfourier_summary_usage.txt ----- 
# this is <stfinvfourier_description_usage.txt>
# ----------------------------------------------------------------------------
#
# Procedures in the Fourier domain
# --------------------------------
Options and parameters in common for procedures in the Fourier domain:
  fpow2              use power of two for number of coefficients
  fdiv=d             use integer multiple of d for number of coefficients
  fpad=f             padding factor
  tshift=d           delay source correction filter wavelet by d (in seconds)
                     this is applied only to the waveform of the impulse
                     response of the correction filter after application to
                     synthetic data (has no effect on convolved data)
  irtap=t1,t2,t3,t4  taper impulse response of correction filter
                     this is applied to the correction filter before
                     application to synthetic data (affects convolved data)

These options define the number of samples N used for the FFT (Fast Fourier
Transform). This number N should be larger than the number of samples M in the
original input time series to avoid wrap-around. If fpow2 is set, N will be
the next power of 2 larger than M*f. Else if fdiv is set, N will be the next
integer multiple of d larger than M*f. If fdiv is not set explicitely, a
default value for d (commonly 100) is used. If option fpad ist used, N will be
f times larger than without padding. Without explicitely setting fpad, a
default value for f is used (which commonly equals 1.5). When defining the
number of samples N, first padding is considered (fpad), then the either
selection of a power of two (pow2) or the divisor criterion (fdiv) is applied.
The latter is only applied, if pow2 ist not selected.

Input time series with M samples will be padded with (N-M) zeros to create the
time series which actually will be transformed to the Fourier domain. Upon
inverse FFT the additional (N-M) samples of the resulting time series will be
discarded before returning the M remaining samples to the caller. Note, that
this is a form of implicite taper. In particular the caller will not obtain
exactly the filter response, which was used for convolution internally.

The derived correction filter in some cases can contain acausal components.
This means that the impulse response is non-zero for negative time values.
Since by definition, the impulse response is output for the time interval of
the input data, these acausal components can remain unnoticed. The option
tshift can be used to shift the impulse response as obtained by inverse FFT in
order to expose acausal components. The time shift is only applied to the
impulse response of the correction filter before returning the impulse
response to the caller. It does not affect the convolved synthetics. If the
output impulse response is discussed with respect to the input synthetics, the
option tshift effectively shifts the temporal origin away from the origin time
specified in the input data.

A time domain taper can be applied to the impulse response of the correction
filter by using option irtap. Four time values are given in units of seconds:
t1, t2, t3, and t4. They must be in increasing order and (t4-t1) must be
smaller than the total duration of the time series used to represent signals
internally. Times value are allowed to be negative. Time series are understood
to be periodic (due to discrete Fourier transformation). Prior to application
of the correction filter to the time series passed to the algorithm, the
correction filter is transformed to the time domain, tapered, and then
transformed to the Fourier domain again. The values of the taper are:
  0                              if       t <  t1
  0.5-0.5*cos(pi*(t-t1)/(t2-t1)) if t1 <= t <= t2
  1                              if t2 <  t <  t3
  0.5+0.5*cos(pi*(t-t3)/(t3-t4)) if t3 <= t <= t4
  0                              if       t >  t4

Time values are given in the same unit in which the sampling interval is given
in the input time series. I.e. if sampling interval is specified as a fraction
of seconds (which is standard) then all time values passed as parameters are
also given as fractions or multiples of seconds.
# ----- END OF stfinvfourier_description_usage.txt ----- 

Summary and detailed description for scaling procedure

# this is <stfinvidentity_summary_usage.txt>
# ----------------------------------------------------------------------------
#
Procedure: Scale with amplitude factor
--------------------------------------
Scale signals by a scalar factor (i.e. keep signal shape).

Options and parameters:
  scaleenergy   if flag is set: scale energy
# ----- END OF stfinvidentity_summary_usage.txt ----- 
# this is <stfinvidentity_description_usage.txt>
# ----------------------------------------------------------------------------
#
Procedure: Scale with amplitude factor
--------------------------------------
This procedure does not modify the waveform of the synthetic data. It
convolves the signales with a discrete delta pulse so to speak. Optionally the
synthetics are scaled with an amplitude factor such that the weighted average
signal energy of all traces in the scaled synthetics equals that of the
recordings. Offset dependent weights are applied in this case. The
appropriately scaled discrete delta pulse is returned as correction filter.

Options and parameters:
  scaleenergy   if flag is set: scale energy
# ----- END OF stfinvidentity_description_usage.txt ----- 

Summary and detailed description for Fourier least squares procedure

# this is <stfinvfdleastsquares_summary_usage.txt>
# ----------------------------------------------------------------------------
#
Procedure: Least squares in the frequency domain
------------------------------------------------
Find optimal correction filter by waterlevel-deconvolution in the Fourier
domain (i.e. least squares with regularization).

Options and parameters:
  waterlevel=l  waterlevel to be applied for regularization.
# ----- END OF stfinvfdleastsquares_summary_usage.txt ----- 
# this is <stfinvfdleastsquares_description_usage.txt>
# ----------------------------------------------------------------------------
#
Procedure: Least squares in the frequency domain
------------------------------------------------
This procedure calculates a least squares fit in the Fourier domain, which
equals a waterlevel-deconvolution. A waterlevel as a fraction of the signal
energy of the input synthetics is applied. If selected, offset dependent
weights will be applied to the signals from different receivers.

Options and parameters:
  waterlevel=l  waterlevel to be applied for regularization.

The procedure is implemented in the Fourier domain. Fourier coefficients for
the correction filter are constructed such that the residual between Fourier
coefficients of recorded signals and Fourier coefficients of synthetic signals
after application of the correction filter are minimized in a least-squares
sense. The least-squares solution is subject to a damping constraint. If the
weighted power spectral density of the synthetics at a given frequency drops
below the l-th fraction (waterlevel parameter) of the average (over all
frequencies) power spectral density, the filter coefficients are damped
(artificially made smaller) by the damping constraint (i.e. waterlevel as seen
from the perspective of deconvolution). 

Start with l=0.01 as a reasonable initial value for the waterlevel. The best
fitting waterlevel must be searched by trial and error and depends of signal
noise and on how well the synthetics describe the actually observed wave
propagtion.

The theory behind the Fourier domain least squares procedure is outlined by
Lisa Groos (2013, Appendix F, page 146). She also describes a way to find an
approrpiate waterlevel by application of the L-curve criterion (Groos, 2013,
Appendix G, page 148).

Lisa Groos. 2013. 2D full waveform inversion of shallow seismic Rayleigh waves.
  Dissertation, Karlsruhe Institute of Technology.
  http://nbn-resolving.org/urn:nbn:de:swb:90-373206
# ----- END OF stfinvfdleastsquares_description_usage.txt -----