Changeset 6281 in ntrip


Ignore:
Timestamp:
Nov 1, 2014, 8:55:34 AM (9 years ago)
Author:
mervart
Message:
 
Location:
trunk/BNC/src/rinex
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/BNC/src/rinex/reqcanalyze.cpp

    r6280 r6281  
    171171        t_satObs satObs;
    172172        t_rnxObsFile::setObsFromRnx(obsFile, _currEpo, rnxSat, satObs);
    173         t_qcObs qcObs;
    174         if (setQcObs(satObs, qcObs) == success) {
    175           qcEpo._qcObs[satObs._prn] = qcObs;
    176           updateQcSat(qcObs, _qcFile._qcSat[satObs._prn]);
    177         }
     173        t_qcObs& qcObs = qcEpo._qcObs[satObs._prn];
     174        setQcObs(qcEpo._epoTime, xyzSta, satObs, qcObs);
     175        updateQcSat(qcObs, _qcFile._qcSat[satObs._prn]);
    178176      }
    179177      _qcFile._qcEpo.push_back(qcEpo);
    180178    }
     179
     180    analyzeMultipath();
    181181
    182182    preparePlotData(obsFile);
    … …  
    264264//
    265265////////////////////////////////////////////////////////////////////////////
    266 t_irc t_reqcAnalyze::setQcObs(const t_satObs& satObs, t_qcObs& qcObs) {
    267 
    268   if (satObs._prn.system() == 'R') {
    269     t_ephGlo* ephGlo  = 0;
    270     for (int ie = 0; ie < _ephs.size(); ie++) {
    271       if (_ephs[ie]->prn() == satObs._prn) {
    272         ephGlo = dynamic_cast<t_ephGlo*>(_ephs[ie]);
    273         break;
    274       }
    275     }
    276     if (ephGlo) {
     266void t_reqcAnalyze::setQcObs(const bncTime& epoTime, const ColumnVector& xyzSta,
     267                             const t_satObs& satObs, t_qcObs& qcObs) {
     268
     269  t_eph* eph = 0;
     270  for (int ie = 0; ie < _ephs.size(); ie++) {
     271    if (_ephs[ie]->prn() == satObs._prn) {
     272      eph = _ephs[ie];
     273      break;
     274    }
     275  }
     276  if (eph) {
     277    ColumnVector xc(4);
     278    ColumnVector vv(3);
     279    if (xyzSta.size() && eph->getCrd(epoTime, xc, vv, false) == success) {
     280      double rho, eleSat, azSat;
     281      topos(xyzSta(1), xyzSta(2), xyzSta(3), xc(1), xc(2), xc(3), rho, eleSat, azSat);
     282      qcObs._eleSet = true;
     283      qcObs._azDeg  = azSat  * 180.0/M_PI;
     284      qcObs._eleDeg = eleSat * 180.0/M_PI;
     285    }
     286    if (satObs._prn.system() == 'R') {
    277287      qcObs._slotSet = true;
    278       qcObs._slotNum = ephGlo->slotNum();
    279     }
    280   }
    281 
    282   bool okFlag = false;
     288      qcObs._slotNum = eph->slotNum();
     289    }
     290  }
     291
     292  // Check Gaps
     293  // ----------
     294  if (qcObs._lastObsTime.valid()) {
     295    double dt = epoTime - qcObs._lastObsTime;
     296    if (dt > 1.5 * _qcFile._interval) {
     297      qcObs._gapL1 = true;
     298      qcObs._gapL2 = true;
     299    }
     300  }
     301  qcObs._lastObsTime = epoTime;
    283302
    284303  // Availability and Slip Flags
    … …  
    288307  double P1 = 0.0;
    289308  double P2 = 0.0;
    290 
    291309  for (unsigned iFrq = 0; iFrq < satObs._obs.size(); iFrq++) {
    292310    const t_frqObs* frqObs = satObs._obs[iFrq];
    293311    if      (frqObs->_rnxType2ch[0] == '1') {
    294312      if (frqObs->_phaseValid) {
    295         L1              = frqObs->_phase;
     313        L1            = frqObs->_phase;
    296314        qcObs._hasL1  = true;
    297315        qcObs._slipL1 = frqObs->_slip;
    … …  
    307325              (satObs._prn.system() == 'E' && frqObs->_rnxType2ch[0] == '5') ) {
    308326      if (frqObs->_phaseValid) {
    309         L2             = frqObs->_phase;
    310         qcObs._hasL2 = true;
     327        L2            = frqObs->_phase;
     328        qcObs._hasL2  = true;
    311329        qcObs._slipL2 = frqObs->_slip;
    312330      }
    … …  
    320338  }
    321339
    322   // Compute the Multipath
    323   // ----------------------
    324   if (L1 != 0.0 && L2 != 0.0) {
     340  // Compute the Multipath Linear Combination
     341  // ----------------------------------------
     342  if (L1 != 0.0 && L2 != 0.0 && P1 != 0.0 && P2 != 0.0) {
    325343    double f1 = 0.0;
    326344    double f2 = 0.0;
    … …  
    341359    L2 = L2 * t_CST::c / f2;
    342360
    343     if (P1 != 0.0) {
    344       qcObs._MP1 = P1 - L1 - 2.0*f2*f2/(f1*f1-f2*f2) * (L1 - L2);
    345       okFlag = true;
    346     }
    347     if (P2 != 0.0) {
    348       qcObs._MP2 = P2 - L2 - 2.0*f1*f1/(f1*f1-f2*f2) * (L1 - L2);
    349       okFlag = true;
    350     }
    351   }
    352 
    353   if (okFlag) {
    354     return success;
    355   }
    356   else {
    357     return failure;
    358   }
    359 }
    360 
    361 //
    362 ////////////////////////////////////////////////////////////////////////////
    363 void t_reqcAnalyze::preparePlotData(const t_rnxObsFile* obsFile) {
    364 
    365   ColumnVector xyzSta = obsFile->xyz();
    366 
    367   QVector<t_polarPoint*>* dataMP1  = new QVector<t_polarPoint*>;
    368   QVector<t_polarPoint*>* dataMP2  = new QVector<t_polarPoint*>;
    369   QVector<t_polarPoint*>* dataSNR1 = new QVector<t_polarPoint*>;
    370   QVector<t_polarPoint*>* dataSNR2 = new QVector<t_polarPoint*>;
     361    qcObs._rawMP1 = P1 - L1 - 2.0*f2*f2/(f1*f1-f2*f2) * (L1 - L2);
     362    qcObs._rawMP2 = P2 - L2 - 2.0*f1*f1/(f1*f1-f2*f2) * (L1 - L2);
     363    qcObs._mpSet  = true;
     364  }
     365}
     366
     367//
     368////////////////////////////////////////////////////////////////////////////
     369void t_reqcAnalyze::analyzeMultipath() {
    371370
    372371  const double SLIPTRESH = 10.0;  // cycle-slip threshold (meters)
    373372  const double chunkStep = 600.0; // 10 minutes
    374 
    375   bncSettings settings;
    376   QString reqSkyPlotSystems = settings.value("reqcSkyPlotSystems").toString();
    377   bool plotGPS = false;
    378   bool plotGlo = false;
    379   bool plotGal = false;
    380   if      (reqSkyPlotSystems == "GPS") {
    381     plotGPS = true;
    382   }
    383   else if (reqSkyPlotSystems == "GLONASS") {
    384     plotGlo = true;
    385   }
    386   else if (reqSkyPlotSystems == "Galileo") {
    387     plotGal = true;
    388   }
    389   else {
    390     plotGPS = true;
    391     plotGlo = true;
    392     plotGal = true;
    393   }
    394373
    395374  // Loop over all satellites available
    … …  
    406385         chunkStart < _qcFile._endTime; chunkStart += chunkStep) {
    407386
    408       // Chunk (sampled) Epoch
    409       // ---------------------
    410       _qcFile._qcEpoSampled.push_back(t_qcEpo());
    411       t_qcEpo& qcEpoSampled = _qcFile._qcEpoSampled.back();
    412       t_qcObs& qcObsSampled = qcEpoSampled._qcObs[prn];
    413 
    414       QVector<double> MP1;
    415       QVector<double> MP2;
     387      bncTime chunkEnd = chunkStart + chunkStep;
     388
     389      QVector<t_qcObs*> obsVec;
     390      QVector<double>   MP1;
     391      QVector<double>   MP2;
    416392   
    417393      // Loop over all Epochs within one Chunk of Data
    418394      // ---------------------------------------------
    419       bncTime prevTime;
    420395      for (int iEpo = 0; iEpo < _qcFile._qcEpo.size(); iEpo++) {
    421396        t_qcEpo& qcEpo = _qcFile._qcEpo[iEpo];
    422         if (qcEpo._epoTime < chunkStart) {
    423           continue;
    424         }
    425         if (!qcEpo._qcObs.contains(prn)) {
    426           continue;
    427         }
    428      
    429         t_qcObs& qcObs = qcEpo._qcObs[prn];
    430      
    431         // Compute the Azimuth and Zenith Distance
    432         // ---------------------------------------
    433         if (chunkStart == qcEpo._epoTime) {
    434           qcEpoSampled._epoTime = qcEpo._epoTime;
    435           qcEpoSampled._PDOP    = qcEpo._PDOP;
    436 
    437           if (qcObs._slotSet) {
    438             qcObsSampled._slotSet = true;
    439             qcObsSampled._slotNum = qcObs._slotNum;
    440           }
    441           if (xyzSta.size()) {
    442             t_eph* eph = 0;
    443             for (int ie = 0; ie < _ephs.size(); ie++) {
    444               if (_ephs[ie]->prn() == prn) {
    445                 eph = _ephs[ie];
    446                 break;
    447               }
    448             }
    449             if (eph) {
    450               ColumnVector xc(4);
    451               ColumnVector vv(3);
    452               if (eph->getCrd(qcEpo._epoTime, xc, vv, false) == success) {
    453                 double rho, eleSat, azSat;
    454                 topos(xyzSta(1), xyzSta(2), xyzSta(3), xc(1), xc(2), xc(3), rho, eleSat, azSat);
    455                 qcObsSampled._azDeg  = azSat  * 180.0/M_PI;
    456                 qcObsSampled._eleDeg = eleSat * 180.0/M_PI;
    457                 qcObs._azDeg         = azSat  * 180.0/M_PI;
    458                 qcObs._eleDeg        = eleSat * 180.0/M_PI;
    459               }
     397        if (chunkStart <= qcEpo._epoTime && qcEpo._epoTime < chunkEnd) {
     398          if (qcEpo._qcObs.contains(prn)) {
     399            t_qcObs& qcObs = qcEpo._qcObs[prn];
     400            obsVec << &qcObs;
     401            if (qcObs._mpSet) {
     402              MP1 << qcObs._rawMP1;
     403              MP2 << qcObs._rawMP2;
    460404            }
    461405          }
    462406        }
    463      
    464         // Check Interval
    465         // --------------
    466         if (prevTime.valid()) {
    467           double dt = qcEpo._epoTime - prevTime;
    468           if (dt > 1.5 * _qcFile._interval) {
    469             qcObsSampled._gapL1 = true;
    470             qcObsSampled._gapL2 = true;
    471           }
    472         }
    473         prevTime = qcEpo._epoTime;
    474      
    475         // Check L1 and L2 availability
    476         // ----------------------------
    477         if (qcObs._hasL1) {
    478           qcObsSampled._hasL1 = true;
    479         }
    480         if (qcObs._hasL2) {
    481           qcObsSampled._hasL2 = true;
    482         }
    483      
    484         // Check Minimal Signal-to-Noise Ratio
    485         // -----------------------------------
    486         if ( qcObs._SNR1 > 0 && (qcObsSampled._SNR1 == 0 || qcObsSampled._SNR1 > qcObs._SNR1) ) {
    487           qcObsSampled._SNR1 = qcObs._SNR1;
    488         }
    489         if ( qcObs._SNR2 > 0 && (qcObsSampled._SNR2 == 0 || qcObsSampled._SNR2 > qcObs._SNR2) ) {
    490           qcObsSampled._SNR2 = qcObs._SNR2;
    491         }
    492      
    493         // Check Slip Flags
    494         // ----------------
    495         if (qcObs._slipL1) {
    496           qcObsSampled._slipL1 = true;
    497         }
    498         if (qcObs._slipL2) {
    499           qcObsSampled._slipL2 = true;
    500         }
    501      
    502         MP1 << qcObs._MP1;
    503         MP2 << qcObs._MP2;
    504       }
    505 
    506       // Compute the Multipath
    507       // ---------------------
    508       if ( MP1.size() > 0 && MP2.size() > 0 &&
    509            ( (prn.system() == 'G' && plotGPS                         ) ||
    510              (prn.system() == 'R' && plotGlo && qcObsSampled._slotSet) ||
    511              (prn.system() == 'E' && plotGal                         ) ) ) {
    512 
     407      }
     408
     409      // Compute the multipath mean and standard deviation
     410      // -------------------------------------------------
     411      if (MP1.size() > 1) {
    513412        double meanMP1 = 0.0;
    514413        double meanMP2 = 0.0;
    … …  
    521420
    522421        bool slipMP = false;
     422
    523423        double stdMP1 = 0.0;
    524424        double stdMP2 = 0.0;
    … …  
    535435
    536436        if (slipMP) {
    537           qcObsSampled._slipL1 = true;
    538           qcObsSampled._slipL2 = true;
     437          stdMP1 = 0.0;
     438          stdMP2 = 0.0;
    539439          qcSat._numSlipsFound += 1;
    540440        }
    541441        else {
    542           stdMP1 = sqrt(stdMP1 / MP1.size());
    543           stdMP2 = sqrt(stdMP2 / MP2.size());
    544           (*dataMP1)  << (new t_polarPoint(qcObsSampled._azDeg, 90.0 - qcObsSampled._eleDeg, stdMP1));
    545           (*dataMP2)  << (new t_polarPoint(qcObsSampled._azDeg, 90.0 - qcObsSampled._eleDeg, stdMP2));
    546         }
    547       }
    548    
    549       // Signal-to-Noise Ratio Plot Data
    550       // -------------------------------
     442          stdMP1 = sqrt(stdMP1 / (MP1.size()-1));
     443          stdMP2 = sqrt(stdMP2 / (MP2.size()-1));
     444        }
     445
     446        for (int ii = 0; ii < obsVec.size(); ii++) {
     447          t_qcObs* qcObs = obsVec[ii];
     448          if (slipMP) {
     449            qcObs->_slipL1 = true;
     450            qcObs->_slipL2 = true;
     451          }
     452          else {
     453            qcObs->_stdMP1 = stdMP1;
     454            qcObs->_stdMP2 = stdMP2;
     455          }
     456        }
     457      }
     458    }
     459  }
     460}
     461
     462//
     463////////////////////////////////////////////////////////////////////////////
     464void t_reqcAnalyze::preparePlotData(const t_rnxObsFile* obsFile) {
     465
     466  QVector<t_polarPoint*>* dataMP1  = new QVector<t_polarPoint*>;
     467  QVector<t_polarPoint*>* dataMP2  = new QVector<t_polarPoint*>;
     468  QVector<t_polarPoint*>* dataSNR1 = new QVector<t_polarPoint*>;
     469  QVector<t_polarPoint*>* dataSNR2 = new QVector<t_polarPoint*>;
     470
     471  bncSettings settings;
     472  QString reqSkyPlotSystems = settings.value("reqcSkyPlotSystems").toString();
     473  bool plotGPS = false;
     474  bool plotGlo = false;
     475  bool plotGal = false;
     476  if      (reqSkyPlotSystems == "GPS") {
     477    plotGPS = true;
     478  }
     479  else if (reqSkyPlotSystems == "GLONASS") {
     480    plotGlo = true;
     481  }
     482  else if (reqSkyPlotSystems == "Galileo") {
     483    plotGal = true;
     484  }
     485  else {
     486    plotGPS = true;
     487    plotGlo = true;
     488    plotGal = true;
     489  }
     490
     491  // Loop over all observations
     492  // --------------------------
     493  for (int iEpo = 0; iEpo < _qcFile._qcEpo.size(); iEpo++) {
     494    t_qcEpo& qcEpo = _qcFile._qcEpo[iEpo];
     495    QMapIterator<t_prn, t_qcObs> it(qcEpo._qcObs);
     496    while (it.hasNext()) {
     497      it.next();
     498      const t_prn&   prn   = it.key();
     499      const t_qcObs& qcObs = it.value();
    551500      if ( (prn.system() == 'G' && plotGPS) ||
    552501           (prn.system() == 'R' && plotGlo) ||
    553502           (prn.system() == 'E' && plotGal) ) {
    554         (*dataSNR1) << (new t_polarPoint(qcObsSampled._azDeg, 90.0 - qcObsSampled._eleDeg, qcObsSampled._SNR1));
    555         (*dataSNR2) << (new t_polarPoint(qcObsSampled._azDeg, 90.0 - qcObsSampled._eleDeg, qcObsSampled._SNR2));
     503
     504        (*dataSNR1) << (new t_polarPoint(qcObs._azDeg, 90.0 - qcObs._eleDeg, qcObs._SNR1));
     505        (*dataSNR2) << (new t_polarPoint(qcObs._azDeg, 90.0 - qcObs._eleDeg, qcObs._SNR2));
     506
     507        (*dataMP1)  << (new t_polarPoint(qcObs._azDeg, 90.0 - qcObs._eleDeg, qcObs._stdMP1));
     508        (*dataMP2)  << (new t_polarPoint(qcObs._azDeg, 90.0 - qcObs._eleDeg, qcObs._stdMP2));
    556509      }
    557510    }

  • trunk/BNC/src/rinex/reqcanalyze.h

    r6277 r6281  
    7676      _gapL2    = false;
    7777      _slotSet  = false;
     78      _eleSet   = false;
     79      _mpSet    = false;
    7880      _slotNum  = 0;
    79       _MP1      = 0.0;
    80       _MP2      = 0.0;
     81      _rawMP1   = 0.0;
     82      _rawMP2   = 0.0;
     83      _stdMP1   = 0.0;
     84      _stdMP2   = 0.0;
    8185      _SNR1     = 0.0;
    8286      _SNR2     = 0.0;
    … …  
    8488      _azDeg    = 0.0;
    8589    }
    86     t_irc set(const t_satObs& obs);
    87     bool   _hasL1;
    88     bool   _hasL2;
    89     bool   _slipL1;
    90     bool   _slipL2;
    91     bool   _gapL1;
    92     bool   _gapL2;
    93     bool   _slotSet;
    94     int    _slotNum;
    95     double _MP1;
    96     double _MP2;
    97     double _SNR1;
    98     double _SNR2;
    99     double _eleDeg;
    100     double _azDeg;
     90    bncTime _lastObsTime;
     91    bool    _hasL1;
     92    bool    _hasL2;
     93    bool    _slipL1;
     94    bool    _slipL2;
     95    bool    _gapL1;
     96    bool    _gapL2;
     97    bool    _slotSet;
     98    int     _slotNum;
     99    bool    _mpSet;
     100    double  _rawMP1;
     101    double  _rawMP2;
     102    double  _stdMP1;
     103    double  _stdMP2;
     104    double  _SNR1;
     105    double  _SNR2;
     106    bool    _eleSet;
     107    double  _eleDeg;
     108    double  _azDeg;
    101109  };
    102110 
    … …  
    134142    QMap<t_prn, t_qcSat> _qcSat;
    135143    QVector<t_qcEpo>     _qcEpo;
    136     QVector<t_qcEpo>     _qcEpoSampled;
    137144  };
    138145
    … …  
    149156  void   updateQcSat(const t_qcObs& qcObs, t_qcSat& qcSat);
    150157
    151   t_irc  setQcObs(const t_satObs& satObs, t_qcObs& qcObs);
     158  void   setQcObs(const bncTime& epoTime, const ColumnVector& xyzSta,
     159                  const t_satObs& satObs, t_qcObs& qcObs);
     160
     161  void   analyzeMultipath();
    152162
    153163  void   preparePlotData(const t_rnxObsFile* obsFile);
Note: See TracChangeset for help on using the changeset viewer.