onefreqcsdetector.hpp

gps源代码

/**
 * @file OneFreqCSDetector.hpp
 * This is a class to detect cycle slips using observables in just one frequency.
 */

#ifndef ONEFREQCSDETECTOR_GPSTK
#define ONEFREQCSDETECTOR_GPSTK

//============================================================================
//
//  This file is part of GPSTk, the GPS Toolkit.
//
//  The GPSTk is free software; you can redistribute it and/or modify
//  it under the terms of the GNU Lesser General Public License as published
//  by the Free Software Foundation; either version 2.1 of the License, or
//  any later version.
//
//  The GPSTk 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 Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with GPSTk; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//  
//  Dagoberto Salazar - gAGE ( http://www.gage.es ). 2007
//
//============================================================================



#include "DataStructures.hpp"



namespace gpstk
{

    /** @addtogroup GPSsolutions */
    //@{


    /** This is a class to detect cycle slips using observables in just one frequency.
     * This class is meant to be used with the GNSS data structures objects
     * found in "DataStructures" class.
     *
     * A typical way to use this class follows:
     *
     * @code
     *   RinexObsStream rin("ebre0300.02o");
     *
     *   gnssRinex gRin;
     *   OneFreqCSDetector markCSC1;
     *
     *   while(rin >> gRin) {
     *      gRin >> markCSC1;
     *   }
     * @endcode
     *
     * The "OneFreqCSDetector" object will visit every satellite in the GNSS data
     * structure that is "gRin" and will decide if a cycle slip has happened in the
     * given observable.
     *
     * By default, the algorithm will use C1 and L1 observables, and the LLI1 index.
     * The result (a 0 if a cycle slip is found, 1 otherwise) will be stored in the
     * data structure as the CSL1 index. Note that these data types may be changed
     * using the appropriate methods. For example:
     *
     * @code
     *    markCSC1.setCodeType(TypeID::P2);
     *    markCSC1.setPhaseType(TypeID::L2);
     *    markCSC1.setLLIType(TypeID::LLI2);
     *    markCSC1.setResultType(TypeID::CSI2);
     * @endcode
     *
     * This algorithm will compute the bias between code and phase, and will compare
     * it with a mean bias that is computed on the fly. If the current bias exceeds a
     * given threshold, then a cycle slip is declared. The algorithm will also use
     * the corresponding LLI index (and the RINEX epoch flag, if present) to guide
     * its decision.
     *
     * The threshold, as well as the filter window size and the maximum allowed 
     * time interval between two successive measures, may be tuned with their 
     * corresponding methods. For instance:
     *
     * @code
     *    markCSC1.setMaxNumSigmas(3.5);
     *    markCSC1.setMaxWindowSize(20);
     * @endcode
     *
     * Please be aware that the window size should not be too big, because 
     * other factors (such as the ionospheric drift) may show up in the bias, affecting
     * the algorithm. When using 1 Hz data sampling, a window size between 60 and 100
     * samples will be fine. 
     *
     * When used with the ">>" operator, this class returns the same incoming
     * data structure with the cycle slip index inserted along their corresponding
     * satellites. Be warned that if a given satellite does not have the 
     * observations required, it will be summarily deleted from the data
     * structure.
     *
     * @sa LICSDetector.hpp and MWCSDetector.hpp for more information.
     *
     * \warning Cycle slip detectors are objets that store their internal state,
     * so you MUST NOT use the SAME object to process DIFFERENT data streams.
     *
     */    
    class OneFreqCSDetector
    {
    public:

        /// Default constructor, setting default parameters and C1 and L1 observables.
        OneFreqCSDetector() : codeType(TypeID::C1), phaseType(TypeID::L1), lliType(TypeID::LLI1), resultType(TypeID::CSL1), deltaTMax(31.0), maxWindowSize(60), maxNumSigmas(4.5), defaultBiasSigma(4.0) {};


        /** Common constructor
         *
         * @param codeT         Type of code to be used.
         * @param dtMax         Maximum interval of time allowed between two successive epochs.
         * @param mwSize        Maximum  size of filter window, in samples.
         * @param mnSigmas      Maximum deviation allowed before declaring cycle slip (in number of sigmas).
         * @param dbSigma       Default value assigned to sigma when filter starts.
         */
        OneFreqCSDetector(const TypeID& codeT, const double& dtMax = 31.0, const int& mwSize = 60, const double& mnSigmas = 4.5, const double& dbSigma = 4.0) : codeType(codeT), deltaTMax(dtMax), maxNumSigmas(mnSigmas), defaultBiasSigma(dbSigma)
        {
            // Don't allow window sizes less than 1
            if (mwSize > 1) maxWindowSize = mwSize; else maxWindowSize = 60;

            switch ( codeT.type )
            {
                case TypeID::C1:
                    phaseType   = TypeID::L1;
                    lliType     = TypeID::LLI1;
                    resultType  = TypeID::CSL1;
                    break;
                case TypeID::C2:
                    phaseType   = TypeID::L2;
                    lliType     = TypeID::LLI2;
                    resultType  = TypeID::CSL2;
                    break;
                case TypeID::C5:
                    phaseType   = TypeID::L5;
                    lliType     = TypeID::LLI5;
                    resultType  = TypeID::CSL5;
                    break;
                case TypeID::C6:
                    phaseType   = TypeID::L6;
                    lliType     = TypeID::LLI6;
                    resultType  = TypeID::CSL6;
                    break;
                case TypeID::C7:
                    phaseType   = TypeID::L7;
                    lliType     = TypeID::LLI7;
                    resultType  = TypeID::CSL7;
                    break;
                case TypeID::C8:
                    phaseType   = TypeID::L8;
                    lliType     = TypeID::LLI8;
                    resultType  = TypeID::CSL8;
                    break;
                default:
                    phaseType   = TypeID::L1;
                    lliType     = TypeID::LLI1;
                    resultType  = TypeID::CSL1;
                };
        };


        /** Returns a satTypeValueMap object, adding the new data generated when calling this object.
         *
         * @param epoch     Time of observations.
         * @param gData     Data object holding the data.
         * @param epochflag Epoch flag.
         */
        virtual satTypeValueMap& Detect(const DayTime& epoch, satTypeValueMap& gData, const short& epochflag=0)
        {
            double value1(0.0);
            double value2(0.0);

            SatIDSet satRejectedSet;

            // Loop through all the satellites
            satTypeValueMap::iterator it;
            for (it = gData.begin(); it != gData.end(); ++it) 
            {
                try
                {
                    // Try to extract the values
                    value1 = (*it).second(codeType);
                    value2 = (*it).second(phaseType);
                }
                catch(...)
                {
                    // If some value is missing, then schedule this satellite for removal
                    satRejectedSet.insert( (*it).first );
                    continue;
                }
                // If everything is OK, then get the new value inside the structure
                // This way of doing it allows concatenation of several different cycle slip detectors
                (*it).second[resultType] += getDetection(epoch, (*it).first, (*it).second, epochflag, value1, value2);
                if ( (*it).second[resultType] > 1.0 ) (*it).second[resultType] = 1.0;
            }
            // Remove satellites with missing data
            gData.removeSatID(satRejectedSet);

            return gData;
        };


        /** Method to set the default code type to be used.
         * @param codeT     TypeID of code to be used
         */
        virtual void setCodeType(const TypeID& codeT)
        {
           codeType = codeT;
        };


        /// Method to get the default code type being used.
        virtual TypeID getCodeType() const
        {
           return codeType;
        };


        /** Method to set the default phase type to be used.
         * @param phaseT    TypeID of phase to be used
         */
        virtual void setPhaseType(const TypeID& phaseT)
        {
           phaseType = phaseT;
        };


        /// Method to get the default phase type being used.
        virtual TypeID getPhaseType() const
        {
           return phaseType;
        };


        /** Method to set the default LLI to be used.
         * @param lliT    LLI to be used
         */
        virtual void setLLIType(const TypeID& lliT)
        {
           lliType = lliT;
        };


        /// Method to get the default LLI being used.
        virtual TypeID getLLIType() const
        {
           return lliType;
        };