00001 #ifndef C1DRECOBJECTIVE_H 00002 #define C1DRECOBJECTIVE_H 00003 #include "PlottableObjective.h" 00004 #include "RecCalc.h" 00005 #include "SeismicDataComp.h" 00006 #include <boost/shared_ptr.hpp> 00007 00008 namespace gplib 00009 { 00010 //! Calculate the misfit between observed receiver function for a given 1D model by calculating a synthetic receiver function from that model 00011 /*! The constructor takes a few essential parameters that are expected not to change during the 00012 * program (in most cases inversion), e.g. the data to fit. See also #C1DRecObjective for the necessary 00013 * constructor parameters. 00014 * 00015 * There are two ways to calculate the Misfit for a given model, either calling GeneralObjective#CalcPerformance, 00016 * or, for parallel programs, calling the three function PreParallel, SafeParallel and PostParallel in that order. 00017 * In both cases you have to give the parameter #member, which is a ::ttranscribed of size 2n, where 00018 * n is the number of layers of the model. The first n entries are the layer thicknesses in km and the last n entries 00019 * are the S-wave velocities in km/s. Densities and P-Wave velocities are calculated from S-Velocities by hard-coded relationships. 00020 */ 00021 class C1DRecObjective: public PlottableObjective 00022 { 00023 private: 00024 //! The object holding the calculated synthetic data 00025 SeismicDataComp RecSynthData; 00026 //! The slowness used in the calculation 00027 double slowness; 00028 //! The object holding the 1D model for forward calculation 00029 ResPkModel Model; 00030 //! The object used to calculate the synthetics 00031 RecCalc RecCalculator; 00032 //! The relative errorfloor 00033 double errorlevel; 00034 //! The absolute errorfloor is calculated from the relative errorfloor 00035 double errorvalue; 00036 //! Poisson's ratio vp/vs, this might become part of the model 00037 double poisson; 00038 //! We might not want to use the whole receiver function for misfit calculation, the startpoint is calculated from the starttime in SetTimeWindow 00039 int startpoint; 00040 //! We might not want to use the whole receiver function for misfit calculation, the endpoint is calculated from the endtime in SetTimeWindow 00041 int endpoint; 00042 //! The measured data is stored in a shared pointer. 00043 /*! For parallel implementations we often work with several instances of this 00044 * class that all work with the same ObservedData. To save memory we use a shared_pointer 00045 * and to avoid collisions SeismicDataComp is const */ 00046 boost::shared_ptr<const SeismicDataComp> ObservedData; 00047 protected: 00048 const SeismicDataComp &GetObservedData() 00049 { 00050 return *ObservedData; 00051 } 00052 public: 00053 //! return a pointer to a copy of the current object 00054 virtual C1DRecObjective *clone() const 00055 { 00056 return new C1DRecObjective(*this); 00057 } 00058 //! Set the time window used for misfit calculations, start and end are in seconds 00059 void SetTimeWindow(const double start, const double end); 00060 //! Set the errorlevel for fit, this is relative to the maximum amplitude, not for each individual data point 00061 void SetErrorLevel(const double level) 00062 { 00063 errorlevel = level; 00064 const double recmaxamp = *max_element( 00065 ObservedData->GetData().begin(), ObservedData->GetData().end()); 00066 // we specify the error relative to the maximum amplitude 00067 //this is usually the initial correlation peak 00068 // here we calculate its value 00069 errorvalue = errorlevel * recmaxamp; 00070 } 00071 //! Set poisson's ratio, at the moment the same for all layers, used for calculating P-velocity 00072 void SetPoisson(const double ratio) 00073 { 00074 poisson = ratio; 00075 } 00076 //! Write the synthetic data to a sac file with name filename, makes only sense after calculating the misfit 00077 virtual void WriteData(const std::string &filename) 00078 { 00079 RecSynthData.WriteAsSac(filename); 00080 } 00081 //! Write the current model to ascii file for calculations 00082 void WriteModel(const std::string &filename) 00083 { 00084 Model.WriteModel(filename); 00085 } 00086 //! Write the current model to ascii file for plotting 00087 void WritePlot(const std::string &filename) 00088 { 00089 Model.WritePlot(filename); 00090 } 00091 //! We have to write runfiles before parallel execution 00092 virtual void PreParallel(const ttranscribed &member); 00093 //! We also clean up files serially 00094 virtual double PostParallel(const ttranscribed &member); 00095 //! Calculate the misfit between the data calculated from model vector member and measured data given in the constructor. 00096 virtual void SafeParallel(const ttranscribed &member); 00097 C1DRecObjective(const C1DRecObjective &Old); 00098 C1DRecObjective& operator=(const C1DRecObjective& source); 00099 //! The constructor needs a few essential parameters 00100 /*! @param TheRecData Shared pointer to the measured receiver function 00101 * @param myshift the time shift used for calculating the measured receiver function 00102 * @param mysigma sigma used for calculating the measured receiver function 00103 * @param myc water level used for calculating the measured receiver function 00104 * @param myslowness slowness used for calculating the measured receiver function 00105 * @param method The method used to calculate the observed data. Can be specdiv or iterdecon 00106 * @param normalized Is the measured data normalized to an initial correlation peak of 1 00107 */ 00108 C1DRecObjective(boost::shared_ptr<const SeismicDataComp> TheRecData, 00109 const int myshift, const double mysigma, const double myc, 00110 const double myslowness, const RecCalc::trfmethod method = 00111 RecCalc::specdiv, const bool normalized = true); 00112 00113 virtual ~C1DRecObjective(); 00114 friend class AbsVelRecObjective; 00115 }; 00116 } 00117 #endif // C1DRECOBJECTIVE_H
1.5.8