Changeset 6285 in ntrip


Ignore:
Timestamp:
Nov 1, 2014, 2:15:28 PM (10 years ago)
Author:
mervart
Message:
 
Location:
trunk/BNC/src/rinex
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r6284 r6285  
    165165        t_qcObs& qcObs = qcEpo._qcObs[satObs._prn];
    166166        setQcObs(qcEpo._epoTime, xyzSta, satObs, qcObs);
    167         updateQcSat(qcObs, _qcFile._qcSat[satObs._prn]);
     167        updateQcSat(qcObs, _qcFile._qcSatSum[satObs._prn]);
    168168      }
    169169      _qcFile._qcEpo.push_back(qcEpo);
    … …  
    245245//
    246246////////////////////////////////////////////////////////////////////////////
    247 void t_reqcAnalyze::updateQcSat(const t_qcObs& qcObs, t_qcSat& qcSat) {
    248   if (qcObs._hasL1 && qcObs._hasL2) {
    249     qcSat._numObs += 1;
    250   }
    251   if (qcObs._slipL1 && qcObs._slipL2) {
    252     qcSat._numSlipsFlagged += 1;
     247void t_reqcAnalyze::updateQcSat(const t_qcObs& qcObs, t_qcSatSum& qcSatSum) {
     248
     249  for (int ii = 0; ii < qcObs._qcFrq.size(); ii++) {
     250    const t_qcFrq& qcFrq    = qcObs._qcFrq[ii];
     251    t_qcFrqSum&    qcFrqSum = qcSatSum._qcFrqSum[qcFrq._rnxType2ch];
     252    qcFrqSum._numObs += 1;
     253    if (qcFrq._slip) {
     254      qcFrqSum._numSlipsFlagged += 1;
     255    }
     256    if (qcFrq._gap) {
     257      qcFrqSum._numGaps += 1;
     258    }
    253259  }
    254260}
    … …  
    282288  }
    283289
    284   // Check Gaps
    285   // ----------
    286   if (qcObs._lastObsTime.valid()) {
    287     double dt = epoTime - qcObs._lastObsTime;
    288     if (dt > 1.5 * _qcFile._interval) {
    289       qcObs._gapL1 = true;
    290       qcObs._gapL2 = true;
    291     }
    292   }
    293   qcObs._lastObsTime = epoTime;
    294 
    295290  // Availability and Slip Flags
    296291  // ---------------------------
    297   double L1 = 0.0;
    298   double L2 = 0.0;
    299   double P1 = 0.0;
    300   double P2 = 0.0;
    301   for (unsigned iFrq = 0; iFrq < satObs._obs.size(); iFrq++) {
    302     const t_frqObs* frqObs = satObs._obs[iFrq];
    303     if      (frqObs->_rnxType2ch[0] == '1') {
    304       if (frqObs->_phaseValid) {
    305         L1            = frqObs->_phase;
    306         qcObs._hasL1  = true;
    307         qcObs._slipL1 = frqObs->_slip;
    308       }
    309       if (frqObs->_codeValid) {
    310         P1 = frqObs->_code;
    311       }
    312       if (frqObs->_snrValid) {
    313         qcObs._SNR1 = frqObs->_snr;
    314       }
    315     }
    316     else if ( (satObs._prn.system() != 'E' && frqObs->_rnxType2ch[0] == '2') ||
    317               (satObs._prn.system() == 'E' && frqObs->_rnxType2ch[0] == '5') ) {
    318       if (frqObs->_phaseValid) {
    319         L2            = frqObs->_phase;
    320         qcObs._hasL2  = true;
    321         qcObs._slipL2 = frqObs->_slip;
    322       }
    323       if (frqObs->_codeValid) {
    324         P2 = frqObs->_code;
    325       }
    326       if (frqObs->_snrValid) {
    327         qcObs._SNR2 = frqObs->_snr;
    328       }
    329     }
    330   }
    331 
    332   // Compute the Multipath Linear Combination
    333   // ----------------------------------------
    334   if (L1 != 0.0 && L2 != 0.0 && P1 != 0.0 && P2 != 0.0) {
    335     double f1 = 0.0;
    336     double f2 = 0.0;
    337     if      (satObs._prn.system() == 'G') {
    338       f1 = t_CST::freq(t_frequency::G1, 0);
    339       f2 = t_CST::freq(t_frequency::G2, 0);
    340     }
    341     else if (satObs._prn.system() == 'R') {
    342       f1 = t_CST::freq(t_frequency::R1, qcObs._slotNum);
    343       f2 = t_CST::freq(t_frequency::R2, qcObs._slotNum);
    344     }
    345     else if (satObs._prn.system() == 'E') {
    346       f1 = t_CST::freq(t_frequency::E1, 0);
    347       f2 = t_CST::freq(t_frequency::E5, 0);
    348     }
    349 
    350     L1 = L1 * t_CST::c / f1;
    351     L2 = L2 * t_CST::c / f2;
    352 
    353     qcObs._rawMP1 = P1 - L1 - 2.0*f2*f2/(f1*f1-f2*f2) * (L1 - L2);
    354     qcObs._rawMP2 = P2 - L2 - 2.0*f1*f1/(f1*f1-f2*f2) * (L1 - L2);
    355     qcObs._mpSet  = true;
     292  for (unsigned ii = 0; ii < satObs._obs.size(); ii++) {
     293
     294    const t_frqObs* frqObs = satObs._obs[ii];
     295
     296    qcObs._qcFrq.push_back(t_qcFrq());
     297    t_qcFrq& qcFrq = qcObs._qcFrq.back();
     298
     299    qcFrq._rnxType2ch = QString(frqObs->_rnxType2ch.c_str());
     300    qcFrq._SNR        = frqObs->_snr;
     301    qcFrq._slip       = frqObs->_slip;
     302
     303    // Check Gaps
     304    // ----------
     305    if (qcFrq._lastObsTime.valid()) {
     306      double dt = epoTime - qcFrq._lastObsTime;
     307      if (dt > 1.5 * _qcFile._interval) {
     308        qcFrq._gap = true;
     309      }
     310    }
     311    qcFrq._lastObsTime = epoTime;
     312
     313    // Compute the Multipath Linear Combination
     314    // ----------------------------------------
     315    if (frqObs->_phaseValid && frqObs->_codeValid) {
     316      t_frequency::type fA;
     317      t_frequency::type fB;
     318      if      (satObs._prn.system() == 'G') {
     319        if      (frqObs->_rnxType2ch[0] == '1') {
     320          fA = t_frequency::G1;
     321          fB = t_frequency::G2;
     322        }
     323        else if (frqObs->_rnxType2ch[0] == '2') {
     324          fA = t_frequency::G2;
     325          fB = t_frequency::G1;
     326        }
     327      }
     328      else if (satObs._prn.system() == 'R') {
     329        if      (frqObs->_rnxType2ch[0] == '1') {
     330          fA = t_frequency::R1;
     331          fB = t_frequency::R2;
     332        }
     333        else if (frqObs->_rnxType2ch[0] == '2') {
     334          fA = t_frequency::R2;
     335          fB = t_frequency::R1;
     336        }
     337      }
     338      else if (satObs._prn.system() == 'E') {
     339        if      (frqObs->_rnxType2ch[0] == '1') {
     340          fA = t_frequency::E1;
     341          fB = t_frequency::E5;
     342        }
     343        else if (frqObs->_rnxType2ch[0] == '5') {
     344          fA = t_frequency::E5;
     345          fB = t_frequency::E1;
     346        }
     347      }
     348
     349      if (fA != t_frequency::dummy && fB != t_frequency::dummy) {
     350        for (unsigned jj = 0; jj < satObs._obs.size(); jj++) {
     351          const t_frqObs* frqObsB = satObs._obs[jj];
     352          if (frqObsB->_rnxType2ch[0] == t_frequency::toString(fB)[1] &&
     353              frqObsB->_phaseValid && frqObsB->_codeValid) {
     354
     355            double f_a = t_CST::freq(fA, qcObs._slotNum);
     356            double f_b = t_CST::freq(fB, qcObs._slotNum);
     357
     358            double L_a = frqObs->_phase  * t_CST::c / f_a;
     359            double C_a = frqObs->_code;
     360            double L_b = frqObsB->_phase * t_CST::c / f_b;
     361
     362            qcFrq._setMP = true;
     363            qcFrq._rawMP = C_a - L_a - 2.0*f_b*f_b/(f_a*f_a-f_b*f_b) * (L_a - L_b);
     364            break;           
     365          }
     366        }
     367      }
     368
     369    }
    356370  }
    357371}
    … …  
    366380  // Loop over all satellites available
    367381  // ----------------------------------
    368   QMutableMapIterator<t_prn, t_qcSat> it(_qcFile._qcSat);
    369   while (it.hasNext()) {
    370     it.next();
    371     const t_prn& prn   = it.key();
    372     t_qcSat&     qcSat = it.value();
    373 
    374     // Loop over all Chunks of Data
    375     // ----------------------------
    376     for (bncTime chunkStart = _qcFile._startTime;
    377          chunkStart < _qcFile._endTime; chunkStart += chunkStep) {
    378 
    379       bncTime chunkEnd = chunkStart + chunkStep;
    380 
    381       QVector<t_qcObs*> obsVec;
    382       QVector<double>   MP1;
    383       QVector<double>   MP2;
     382  QMutableMapIterator<t_prn, t_qcSatSum> itSat(_qcFile._qcSatSum);
     383  while (itSat.hasNext()) {
     384    itSat.next();
     385    const t_prn& prn      = itSat.key();
     386    t_qcSatSum&  qcSatSum = itSat.value();
     387
     388    // Loop over all frequencies available
     389    // -----------------------------------
     390    QMutableMapIterator<QString, t_qcFrqSum> itFrq(qcSatSum._qcFrqSum);
     391    while (itFrq.hasNext()) {
     392      itFrq.next();
     393      const QString& frqType  = itFrq.key();
     394      t_qcFrqSum&    qcFrqSum = itFrq.value();
     395
     396
     397      // Loop over all Chunks of Data
     398      // ----------------------------
     399      for (bncTime chunkStart = _qcFile._startTime;
     400           chunkStart < _qcFile._endTime; chunkStart += chunkStep) {
     401
     402        bncTime chunkEnd = chunkStart + chunkStep;
     403
     404        QVector<t_qcFrq*> frqVec;
     405        QVector<double>   MP;
    384406   
    385       // Loop over all Epochs within one Chunk of Data
    386       // ---------------------------------------------
    387       for (int iEpo = 0; iEpo < _qcFile._qcEpo.size(); iEpo++) {
    388         t_qcEpo& qcEpo = _qcFile._qcEpo[iEpo];
    389         if (chunkStart <= qcEpo._epoTime && qcEpo._epoTime < chunkEnd) {
    390           if (qcEpo._qcObs.contains(prn)) {
    391             t_qcObs& qcObs = qcEpo._qcObs[prn];
    392             obsVec << &qcObs;
    393             if (qcObs._mpSet) {
    394               MP1 << qcObs._rawMP1;
    395               MP2 << qcObs._rawMP2;
     407        // Loop over all Epochs within one Chunk of Data
     408        // ---------------------------------------------
     409        for (int iEpo = 0; iEpo < _qcFile._qcEpo.size(); iEpo++) {
     410          t_qcEpo& qcEpo = _qcFile._qcEpo[iEpo];
     411          if (chunkStart <= qcEpo._epoTime && qcEpo._epoTime < chunkEnd) {
     412            if (qcEpo._qcObs.contains(prn)) {
     413              t_qcObs& qcObs = qcEpo._qcObs[prn];
     414              for (int iFrq = 0; iFrq < qcObs._qcFrq.size(); iFrq++) {
     415                t_qcFrq& qcFrq = qcObs._qcFrq[iFrq];
     416                if (qcFrq._rnxType2ch == frqType) {
     417                  frqVec << &qcFrq;
     418                  if (qcFrq._setMP) {
     419                    MP << qcFrq._rawMP;
     420                  }
     421                }
     422              }
    396423            }
    397424          }
    398425        }
    399       }
    400 
    401       // Compute the multipath mean and standard deviation
    402       // -------------------------------------------------
    403       if (MP1.size() > 1) {
    404         double meanMP1 = 0.0;
    405         double meanMP2 = 0.0;
    406         for (int ii = 0; ii < MP1.size(); ii++) {
    407           meanMP1 += MP1[ii];
    408           meanMP2 += MP2[ii];
    409         }
    410         meanMP1 /= MP1.size();
    411         meanMP2 /= MP2.size();
    412 
    413         bool slipMP = false;
    414 
    415         double stdMP1 = 0.0;
    416         double stdMP2 = 0.0;
    417         for (int ii = 0; ii < MP1.size(); ii++) {
    418           double diff1 = MP1[ii] - meanMP1;
    419           double diff2 = MP2[ii] - meanMP2;
    420           if (fabs(diff1) > SLIPTRESH || fabs(diff2) > SLIPTRESH) {
    421             slipMP = true;
    422             break;
    423           }
    424           stdMP1 += diff1 * diff1;
    425           stdMP2 += diff2 * diff2;
    426         }
    427 
    428         if (slipMP) {
    429           stdMP1 = 0.0;
    430           stdMP2 = 0.0;
    431           qcSat._numSlipsFound += 1;
    432         }
    433         else {
    434           stdMP1 = sqrt(stdMP1 / (MP1.size()-1));
    435           stdMP2 = sqrt(stdMP2 / (MP2.size()-1));
    436         }
    437 
    438         for (int ii = 0; ii < obsVec.size(); ii++) {
    439           t_qcObs* qcObs = obsVec[ii];
     426
     427        // Compute the multipath mean and standard deviation
     428        // -------------------------------------------------
     429        if (MP.size() > 1) {
     430          double meanMP = 0.0;
     431          for (int ii = 0; ii < MP.size(); ii++) {
     432            meanMP += MP[ii];
     433          }
     434          meanMP /= MP.size();
     435       
     436          bool slipMP = false;
     437       
     438          double stdMP = 0.0;
     439          for (int ii = 0; ii < MP.size(); ii++) {
     440            double diff = MP[ii] - meanMP;
     441            if (fabs(diff) > SLIPTRESH) {
     442              slipMP = true;
     443              break;
     444            }
     445            stdMP += diff * diff;
     446          }
     447       
    440448          if (slipMP) {
    441             qcObs->_slipL1 = true;
    442             qcObs->_slipL2 = true;
     449            stdMP = 0.0;
     450            stdMP = 0.0;
     451            qcFrqSum._numSlipsFound += 1;
    443452          }
    444453          else {
    445             qcObs->_stdMP1 = stdMP1;
    446             qcObs->_stdMP2 = stdMP2;
    447           }
    448         }
    449       }
    450     }
    451   }
     454            stdMP = sqrt(stdMP / (MP.size()-1));
     455          }
     456       
     457          for (int ii = 0; ii < frqVec.size(); ii++) {
     458            t_qcFrq* qcFrq = frqVec[ii];
     459            if (slipMP) {
     460              qcFrq->_slip = true;
     461            }
     462            else {
     463              qcFrq->_stdMP = stdMP;
     464            }
     465          }
     466        }
     467      } // chunk loop
     468    } // frq loop
     469  } // sat loop
    452470}
    453471

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

    r6284 r6285  
    6666 private:
    6767
     68  class t_qcFrq {
     69   public:
     70    t_qcFrq() {
     71      _slip  = false;
     72      _gap   = false;
     73      _setMP = false;
     74      _rawMP = 0.0;
     75      _stdMP = 0.0;
     76      _SNR   = 0.0;
     77    }
     78    QString _rnxType2ch;
     79    bncTime _lastObsTime;
     80    bool    _slip;
     81    bool    _gap;
     82    bool    _setMP;
     83    double  _rawMP;
     84    double  _stdMP;
     85    double  _SNR;
     86  };
     87
    6888  class t_qcObs {
    6989   public:
    7090    t_qcObs() {
    71       _hasL1    = false;
    72       _hasL2    = false;
    73       _slipL1   = false;
    74       _slipL2   = false;
    75       _gapL1    = false;
    76       _gapL2    = false;
    77       _slotSet  = false;
    78       _eleSet   = false;
    79       _mpSet    = false;
    80       _slotNum  = 0;
    81       _rawMP1   = 0.0;
    82       _rawMP2   = 0.0;
    83       _stdMP1   = 0.0;
    84       _stdMP2   = 0.0;
    85       _SNR1     = 0.0;
    86       _SNR2     = 0.0;
    87       _eleDeg   = 0.0;
    88       _azDeg    = 0.0;
     91      _slotSet = false;
     92      _eleSet  = false;
     93      _slotNum = 0;
     94      _eleDeg  = 0.0;
     95      _azDeg   = 0.0;
    8996    }
    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;
     97    bool             _slotSet;
     98    bool             _eleSet;
     99    int              _slotNum;
     100    double           _eleDeg;
     101    double           _azDeg;
     102    QVector<t_qcFrq> _qcFrq;
    109103  };
    110104 
    … …  
    116110  };
    117111 
    118   class t_qcSat {
     112  class t_qcFrqSum {
    119113   public:
    120     t_qcSat() {
     114    t_qcFrqSum() {
    121115      _numObs          = 0;
    122116      _numSlipsFlagged = 0;
    … …  
    129123    int _numGaps;
    130124  };
     125
     126  class t_qcSatSum {
     127   public:
     128    QMap<QString, t_qcFrqSum> _qcFrqSum;
     129  };
    131130 
    132131  class t_qcFile {
    133132   public:
    134133    t_qcFile() {clear();}
    135     void clear() {_qcSat.clear(); _qcEpo.clear();}
    136     bncTime              _startTime;
    137     bncTime              _endTime;
    138     QString              _antennaName;
    139     QString              _markerName;
    140     QString              _receiverType;
    141     double               _interval;
    142     QMap<t_prn, t_qcSat> _qcSat;
    143     QVector<t_qcEpo>     _qcEpo;
     134    void clear() {_qcSatSum.clear(); _qcEpo.clear();}
     135    bncTime                 _startTime;
     136    bncTime                 _endTime;
     137    QString                 _antennaName;
     138    QString                 _markerName;
     139    QString                 _receiverType;
     140    double                  _interval;
     141    QMap<t_prn, t_qcSatSum> _qcSatSum;
     142    QVector<t_qcEpo>        _qcEpo;
    144143  };
    145144
    … …  
    154153  void   analyzeFile(t_rnxObsFile* obsFile);
    155154
    156   void   updateQcSat(const t_qcObs& qcObs, t_qcSat& qcSat);
     155  void   updateQcSat(const t_qcObs& qcObs, t_qcSatSum& qcSatSum);
    157156
    158157  void   setQcObs(const bncTime& epoTime, const ColumnVector& xyzSta,
Note: See TracChangeset for help on using the changeset viewer.