stfinvfdleastsquares.h 9.95 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
/*! \file stfinvfdleastsquares.h
 * \brief least squares in the frequency domain (prototypes)
 * 
 * ----------------------------------------------------------------------------
 * 
 * \author Thomas Forbriger
 * \date 06/05/2011
 * 
 * least squares in the frequency domain (prototypes)
 * 
 * Copyright (c) 2011 by Thomas Forbriger (BFO Schiltach) 
 *
 * ----
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version. 
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 * ----
 *
 * 
 * REVISIONS and CHANGES 
 *  - 06/05/2011   V1.0   Thomas Forbriger
 *  - 30/09/2011   V1.1   implemented handling of additional time series pairs
 *  - 04/10/2011   V1.2   renamed engine
34
 *  - 14/10/2015   V1.3   new end-user usage functions
35
36
37
38
39
40
41
42
 * 
 * ============================================================================
 */

// include guard
#ifndef STFINV_STFINVFDLEASTSQUARES_H_VERSION

#define STFINV_STFINVFDLEASTSQUARES_H_VERSION \
43
  "STFINV_STFINVFDLEASTSQUARES_H   V1.3"
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249

#include<stfinv/stfinvfourier.h>

namespace stfinv {

  /*! \brief Fourier domain least squares engine
   * \ingroup engines
   *
   * \par Concept behind this engine
   * If
   * - \f$d_{lk}\f$ is the Fourier coefficient of recorded data at Frequency
   *   \f$f_l\f$ and receiver \f$k\f$ at offset \f$r_k\f$,
   * - \f$s_{lk}\f$ is the Fourier coefficient of the corresponding
   *   synthetics and
   * - \f$q_l\f$ is that of the sought source tim function,
   * .
   * then this engine will minimize the objective function
   * \f[
   *   E=\sum\limits_{l,k}\left|w_{lk}\,
   *      \left(d_{lk}-s_{lk}q_l\right)
   *   \right|^2+\sum\limits_{l}\lambda^2\left|q_l\right|^2
   *   =\chi^2+\psi^2
   * \f]
   * with respect to the real part \f$q_l^\prime\f$ and the
   * imaginary part \f$q_l^{\prime\prime}\f$ of 
   * \f[
   *   q_l=q_l^\prime+i\,q_l^{\prime\prime}.
   * \f]
   * In the above expression
   * \f[
   *   \chi^2=\sum\limits_{l,k}\left|w_{lk}\,
   *           \left(d_{lk}-s_{lk}q_l\right)
   *            \right|^2
   * \f]
   * is the data misfit with weights \f$w_{lk}\f$ and
   * \f[
   *   \psi^2=\sum\limits_{l}\lambda^2\left|q_l\right|^2
   * \f]
   * is used for regularization and will introduce a water-level in the
   * deconvolution.
   * \f$\lambda\f$ will balance both contributions.
   * The conditions
   * \f[
   *   \frac{\partial E}{\partial q_l^\prime}\stackrel{!}{=}0
   *   \quad\wedge\quad
   *   \frac{\partial E}{\partial q_l^{\prime\prime}}\stackrel{!}{=}0
   * \f]
   * result in (Forbriger, 2001, appendix A.3)
   * \f[
   *   q_l=\frac{
   *     \eta^2\sum\limits_{k}f_k^2\,s_{kl}^\ast\,d_{kl}
   *   }{
   *     \lambda^2+\eta^2\sum\limits_{k}f_k^2\,s_{kl}^\ast\,s_{kl}
   *   }
   *   \quad\forall\, l
   * \f]
   * where
   * \f[
   *   w_{lk}=\eta\,f_k
   * \f]
   * and \f$f_k\f$ is a receiver specific weighting factor.
   * Now \f$\eta\f$ and \f$\lambda\f$ have to be used to balance the
   * regularization.
   * We aim to specify a waterlevel as a fraction of synthetic data energy.
   *  
   * \par Setting up the waterlevel
   * The misfit equals one if the scaled energy of the residual
   * \f$d_{lk}-s_{lk}q_l\f$ equals the scaled energy of the synthetics
   * \f$s_{lk}\f$ and
   * \f[
   *   \eta^2=\frac{1}{\sum\limits_k f_k^2\sum\limits_l \left|s_{lk}\right|^2}
   * \f]
   * is the reciprocal of the scaled energy of the synthetics.
   * If we then choose
   * \f[
   *   \frac{\lambda^2}{\eta^2}=\frac{\epsilon^2}{N\eta^2}=
   *     \frac{\epsilon^2}{N}\sum\limits_k f_k^2\sum\limits_{l=0}^{N-1}
  *     \left|s_{lk}\right|^2
   * \f]
   * where \f$N\f$ is the number of frequencies, then \f$\epsilon^2\f$
   * will specify a waterlevel as a fraction of the scaled energy of the
   * synthetics.
   *
   * \par Using Parceval's Theorem to calculate signal energy
   * Parceval's Theorem for a signal \f$a(t)\f$ and its Fourier transform 
   * \f$\tilde{a}(\omega)\f$ is
   * \f[
   *   \int\limits_{-\infty}^{+\infty}\bigl|a(t)\bigr|^2\,\textrm{d} t=
   *   \int\limits_{-\infty}^{+\infty}\bigl|\tilde{a}(\omega)\bigr|^2\,
   *     \frac{\textrm{d} \omega}{2\pi}.
   * \f]
   * If \f$S_{jk}\f$ are the time series samples corresponding to the Fourier
   * coefficients \f$\tilde{s}_{lk}\f$ and \f$\Delta t\f$ is the sampling
   * interval then
   * \f[
   *   \sum\limits_{k=0}^{M-1}\left|S_{jk}\right|^2\,\Delta t=
   *   \sum\limits_{l=0}^{M-1}\left|\tilde{s}_{lk}\right|^2\,\frac{1}{M\,\Delta t},
   * \f]
   * where \f$M=2N\f$ is the number of samples in the time series.
   * In the above calculation the energy sum only uses the positive
   * frequencies and
   * \f[
   *   \sum\limits_k f_k^2\sum\limits_{l=0}^{N-1}\left|\tilde{s}_{lk}\right|^2
   *   =
   *     N\,(\Delta t)^2\,
   *     \sum\limits_k f_k^2
   *     \sum\limits_{j=0}^{2N-1}\left|S_{jk}\right|^2.
   * \f]
   * Fourier coefficients \f$s_{lk}\f$ calculated by the
   * stfinv::STFFourierDomainEngine are not scaled (see documentation of
   * libfourierxx and libfftw3), such that
   * \f[
   *   \Delta t\,s_{lk}=\tilde{s}_{lk}
   * \f]
   * (both, \f$s_{lk}\f$ and \f$\tilde{s}_{lk}\f$ are Fourier coefficients).
   * Consequently
   * \f[
   *   \sum\limits_k f_k^2\sum\limits_{l=0}^{N-1}\left|s_{lk}\right|^2
   *   =
   *     N\,
   *     \sum\limits_k f_k^2
   *     \sum\limits_{j=0}^{2N-1}\left|S_{jk}\right|^2.
   * \f]
   *
   * \par Final calculation recipe
   * The solution to our problem is
   * \f[
   *   q_l=\frac{
   *     \sum\limits_{k}f_k^2\,s_{lk}^\ast\,d_{lk}
   *   }{
   *     \epsilon^2\,\sum\limits_k f_k^2
   *                \sum\limits_{j=0}^{2N-1}\left|S_{jk}\right|^2
   *           +\sum\limits_{k}f_k^2\,s_{lk}^\ast\,s_{lk}
   *   }
   *   \quad\forall\, l,
   * \f]
   * where
   * \f[
   *   \sum\limits_{j=0}^{2N-1}\left|S_{jk}\right|^2
   * \f]
   * is the sum of the squared sample values \f$S_{jk}\f$ of the synthetic
   * time series for receiver \f$k\f$, \f$f_k\f$ are the scaling factors
   * provided by stfinv::STFBaseEngine::weight(), and \f$\epsilon^2\f$
   * is the water level parameter passed to STFEngineFDLeastSquares.
   *
   * \par References
   * Forbriger, T., 2001. Inversion flachseismischer Wellenfeldspektern.
   *   PhD thesis, University of Stuttgart. 
   *   http://elib.uni-stuttgart.de/opus/volltexte/2001/861/
   *
   * \par Why was it renamed?
   * This engine was called 'blind deconvolution' engine.
   * It was renamed to fourier domain least squares engine in October 2011.
   *
   * This engine understands the recorded data as being the synthetic data
   * being convoled with the STF (source correction filter)
   * and having some noise
   * added.
   * The true impulse response of the subsurface could be obtained by
   * deconvolution of the recorded data with the STF.
   * Neither the impulse response of the subsurface nor the STF are known. 
   * For this reason the STF has to be found by minimizing an objective
   * function.
   * For this reason the it is called 'blind' deconvolution.
   *
   * However, the approach of libstfinv is to convolved the synthetics with
   * the STF to reducde the misfit to the recorded data, not to deconvolve the
   * recorded data. 
   * For this reason, maybe, we should better call the approach 'blind
   * convolution'.
   * With respect to image processing techniques, all engines in libstfinv do
   * some kind of 'blind convolution', not only the blind deconvolution
   * engine.
   * Consequently this engine should better be called 'least-squares engine'.
   */
  class STFEngineFDLeastSquares: public stfinv::STFFourierDomainEngine {
    public:
      //! \brief typedef to refer to base class
      typedef stfinv::STFFourierDomainEngine Tbase;
      //! \brief ID used to select this engine
      static const char* const ID;
      //! \brief short description of this engine
      static const char* const description;
      /*! \brief Constructor.
       */
      STFEngineFDLeastSquares(const stfinv::Tvectoroftriples& triples,
                                  const stfinv::Waveform& stf,
                                  const std::string& parameters)
        : Tbase(triples, stf, parameters), Mwaterlevel(1.e-3)
      { this->initialize(); }
      /*! \brief Constructor.
       */
      STFEngineFDLeastSquares(const stfinv::Tvectoroftriples& triples,
                                  const stfinv::Waveform& stf,
                                  const stfinv::Tvectorofpairs& pairs,
                                  const std::string& parameters)
        : Tbase(triples, stf, pairs, parameters), Mwaterlevel(1.e-3)
      { this->initialize(); }
      //! \brief abstract base requires virtual destructor
      virtual ~STFEngineFDLeastSquares() { }
      //! \brief Start engine 
      virtual void exec();
      //! \brief print online help
      virtual void help(std::ostream& os=std::cout) const;
      //! \brief print online help
      static void classhelp(std::ostream& os=std::cout);
250
251
252
253
      //! \brief print detailed description
      virtual void usage(std::ostream& os=std::cout) const;
      //! \brief print detailed description
      static void classusage(std::ostream& os=std::cout);
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
      //! \brief return name of engine
      virtual const char* name() const;
    private:
      //! \brief initialize work space
      void initialize();

      // member data
    private:
      //! \brief waterlevel
      double Mwaterlevel;
  }; // class STFEngineFDLeastSquares

} // namespace stfinv

#endif // STFINV_STFINVFDLEASTSQUARES_H_VERSION (includeguard)

/* ----- END OF stfinvfdleastsquares.h ----- */