Changeset 8905 in ntrip for trunk/BNC/src/PPP/pppFilter.cpp

Timestamp:

Mar 18, 2020, 11:13:50 AM (4 years ago)

Author:

stuerze

Message:

some developments regarding PPP, not completed!

File:

• trunk/BNC/src/PPP/pppFilter.cpp (modified) (18 diffs)

trunk/BNC/src/PPP/pppFilter.cpp

@@ -34 +34 @@
 // Constructor
 ////////////////////////////////////////////////////////////////////////////
-t_pppFilter::t_pppFilter() {
+t_pppFilter::t_pppFilter(t_pppObsPool* obsPool) {
   _parlist = 0;
+  _numSat = 0;
+  _obsPool = obsPool;
 }
 
@@ -46 +48 @@
 // Process Single Epoch
 ////////////////////////////////////////////////////////////////////////////
-t_irc t_pppFilter::processEpoch(t_pppObsPool* obsPool) {
-
+t_irc t_pppFilter::processEpoch(int num) {
   _numSat     = 0;
   const double maxSolGap = 60.0;
@@ -57 +58 @@
   // Vector of all Observations
   // --------------------------
-  t_pppObsPool::t_epoch* epoch = obsPool->lastEpoch();
+  t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch();
   if (!epoch) {
     return failure;
@@ -75 +76 @@
   string epoTimeStr = string(_epoTime);
 
+  //--
   // Set Parameters
   // --------------
@@ -91 +93 @@
   for (unsigned ii = 0; ii < params.size(); ii++) {
     const t_pppParam* par1 = params[ii];
-
     _x0[ii] = par1->x0();
-
     int iOld = par1->indexOld();
     if (iOld < 0) {
@@ -113 +113 @@
   predictCovCrdPart(QFltOld);
 
+
+  // Pre-Process Satellite Systems separately
+  // ----------------------------------------
+  bool preProcessing = false;
+  if (OPT->_obsModelType == OPT->DCMcodeBias ||
+      OPT->_obsModelType == OPT->DCMphaseBias) {
+    preProcessing = true;
+    for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
+      char system = OPT->systems()[iSys];
+      t_prn refPrn = (_obsPool->getRefSatMapElement(system))->prn();
+      vector<t_pppSatObs*> obsVector;
+      for (unsigned jj = 0; jj < allObs.size(); jj++) {
+        if (allObs[jj]->prn().system() == system) {
+          obsVector.push_back(allObs[jj]);
+        }
+      }
+      if (processSystem(OPT->LCs(system), obsVector, refPrn,
+                        epoch->pseudoObsIono(), preProcessing) != success) {
+        return failure;
+      }
+    }
+  }
+
   // Process Satellite Systems separately
   // ------------------------------------
   for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     char system = OPT->systems()[iSys];
+    t_prn refPrn = (_obsPool->getRefSatMapElement(system))->prn();
     unsigned int num = 0;
     vector<t_pppSatObs*> obsVector;
@@ -126 +150 @@
     }
     LOG << epoTimeStr << " SATNUM " << system << ' ' << right << setw(2) << num << endl;
-    if ( processSystem(OPT->LCs(system), obsVector) != success ) {
+    if (processSystem(OPT->LCs(system), obsVector, refPrn,
+                      epoch->pseudoObsIono(), preProcessing) != success) {
       return failure;
     }
   }
-
-  cmpDOP(allObs);
-
-  _parlist->printResult(_epoTime, _QFlt, _xFlt);
-  _lastEpoTimeOK = _epoTime;  // remember time of last successful epoch processing
+  if (_obsPool->epoReProcessing()) {
+    // set A1 und A2 abhängig von num
+    // if num == 1 => A1
+    // if num >  1 => A2
+  }
+  else {
+    cmpDOP(allObs);
+    _parlist->printResult(_epoTime, _QFlt, _xFlt);
+    _lastEpoTimeOK = _epoTime;  // remember time of last successful epoch processing
+  }
+
   return success;
 }
@@ -141 +172 @@
 ////////////////////////////////////////////////////////////////////////////
 t_irc t_pppFilter::processSystem(const vector<t_lc::type>& LCs,
-                                 const vector<t_pppSatObs*>& obsVector) {
-
+                                 const vector<t_pppSatObs*>& obsVector,
+                                 const t_prn& refPrn,
+                                 bool pseudoObsIonoAvailable,
+                                 bool preProcessing) {
+  qDebug() << "======t_pppFilter::processSystem=======";
   LOG.setf(ios::fixed);
 
   // Detect Cycle Slips
   // ------------------
-  if (detectCycleSlips(LCs, obsVector) != success) {
+  if (detectCycleSlips(LCs, obsVector, refPrn, preProcessing) != success) {
     return failure;
   }
 
-  ColumnVector            xSav       = _xFlt;
-  SymmetricMatrix         QSav       = _QFlt;
-  string                  epoTimeStr = string(_epoTime);
+  unsigned usedLCs = LCs.size(); //qDebug() << "usedLCs: " << usedLCs;
+  if (OPT->_pseudoObsIono && !pseudoObsIonoAvailable) {
+    usedLCs -= 1;  // GIM not used
+  }
+  ColumnVector               xSav       = _xFlt;
+  SymmetricMatrix            QSav       = _QFlt;
+  string                     epoTimeStr = string(_epoTime);
   const vector<t_pppParam*>& params     = _parlist->params();
-  unsigned                maxObs     = obsVector.size() * LCs.size();
+  unsigned                   maxObs     = obsVector.size() * usedLCs;
+  //unsigned                   maxObs     = 2 * usedLCs;
+
+  if (OPT->_pseudoObsIono && pseudoObsIonoAvailable) { // stecDiff w.r.t refSat
+    maxObs -= 1;
+  }
+  qDebug() << "par.size()      : " << _parlist->nPar() << " LCs.size(): " << usedLCs;
+  qDebug() << "obsVector.size(): " << obsVector.size() << " maxObs    : " << maxObs;
 
   // Outlier Detection Loop
   // ----------------------
   for (unsigned iOutlier = 0; iOutlier < maxObs; iOutlier++) {
+    qDebug() << "iOutlier: " << iOutlier;
 
     if (iOutlier > 0) {
@@ -168 +214 @@
     // First-Design Matrix, Terms Observed-Computed, Weight Matrix
     // -----------------------------------------------------------
+    qDebug() << "A(" << maxObs << "," << _parlist->nPar() << ")";
     Matrix                AA(maxObs, _parlist->nPar());
     ColumnVector          ll(maxObs);
     DiagonalMatrix        PP(maxObs); PP = 0.0;
+    //TETSPLOT
+    for (unsigned iPar = 0; iPar < params.size(); iPar++) {
+      const t_pppParam* par = params[iPar];
+      cout << "  " << par->toString() <<endl;
+    }
+    cout << endl;
 
     int iObs = -1;
@@ -176 +229 @@
     vector<t_lc::type>   usedTypes;
     for (unsigned ii = 0; ii < obsVector.size(); ii++) {
-      t_pppSatObs* obs = obsVector[ii];
+      t_pppSatObs* obs = obsVector[ii]; qDebug() << "SATELLITE: " << obs->prn().toString().c_str() << "isRef: " << obs->isReference();
       if (!obs->outlier()) {
-        for (unsigned jj = 0; jj < LCs.size(); jj++) {
+        for (unsigned jj = 0; jj < usedLCs; jj++) {
           const t_lc::type tLC = LCs[jj];
+          if (tLC == t_lc::GIM &&  obs->isReference()) {continue;}
           ++iObs;
           usedObs.push_back(obs);
@@ -186 +240 @@
             const t_pppParam* par = params[iPar];
             AA[iObs][iPar] = par->partial(_epoTime, obs, tLC);
+            cout << setw(10) << par->partial(_epoTime, obs, tLC);
           }
+          cout << endl;
           ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1));
           PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC));
+          //qDebug() << "ll[iObs]: " << ll[iObs];
+          //qDebug() << "PP[iObs]: " << PP[iObs];
         }
       }
@@ -229 +287 @@
       t_pppSatObs* obs = usedObs[maxOutlierIndex];
       t_pppParam* par = 0;
-      LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' '
-          << obs->prn().toString()                        << ' '
-          << setw(8) << setprecision(4) << maxOutlier << endl;
+      if (!preProcessing) {
+        LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' '
+            << obs->prn().toString()                        << ' '
+            << setw(8) << setprecision(4) << maxOutlier << endl;
+      }
       for (unsigned iPar = 0; iPar < params.size(); iPar++) {
         t_pppParam* hlp = params[iPar];
-        if (hlp->type() == t_pppParam::amb && hlp->prn()  == obs->prn() &&
+        if (hlp->type() == t_pppParam::amb &&
+            hlp->prn()  == obs->prn() &&
             hlp->tLC()  == usedTypes[maxOutlierIndex]) {
           par = hlp;
         }
       }
-      if (par) {
+      if      (par && preProcessing) {
+        if (par->prn() == refPrn) {
+          _obsPool->setEpoReProcessing(true);
+        }
+      }
+      else if (par && !preProcessing) {
         if (par->ambResetCandidate()) {
           resetAmb(par->prn(), obsVector, &QSav, &xSav);
@@ -248 +314 @@
         }
       }
-      else {
+      else if (!par && !preProcessing){
         obs->setOutlier();
       }
@@ -256 +322 @@
     // ---------------
     else {
-      for (unsigned jj = 0; jj < LCs.size(); jj++) {
-        for (unsigned ii = 0; ii < usedObs.size(); ii++) {
-          const t_lc::type tLC = usedTypes[ii];
-          t_pppSatObs*        obs = usedObs[ii];
-          if (tLC == LCs[jj]) {
-            obs->setRes(tLC, vv[ii]);
-            LOG << epoTimeStr << " RES "
-                << left << setw(3) << t_lc::toString(tLC) << right << ' '
-                << obs->prn().toString().substr(0,3) << ' '
-                << setw(8) << setprecision(4) << vv[ii] << endl;
+      if (!preProcessing) {
+        for (unsigned jj = 0; jj < LCs.size(); jj++) {
+          for (unsigned ii = 0; ii < usedObs.size(); ii++) {
+            const t_lc::type tLC = usedTypes[ii];
+            t_pppSatObs*        obs = usedObs[ii];
+            if (tLC == LCs[jj]) {
+              obs->setRes(tLC, vv[ii]);
+              LOG << epoTimeStr << " RES "
+                  << left << setw(3) << t_lc::toString(tLC) << right << ' '
+                  << obs->prn().toString().substr(0,3) << ' '
+                  << setw(8) << setprecision(4) << vv[ii] << endl;
+            }
           }
         }
@@ -279 +347 @@
 ////////////////////////////////////////////////////////////////////////////
 t_irc t_pppFilter::detectCycleSlips(const vector<t_lc::type>& LCs,
-                                    const vector<t_pppSatObs*>& obsVector) {
+                                    const vector<t_pppSatObs*>& obsVector,
+                                    const t_prn& refPrn,
+                                    bool preProcessing) {
 
   const double            SLIP       = 20.0;  // slip threshold
   string                  epoTimeStr = string(_epoTime);
-  const vector<t_pppParam*>& params     = _parlist->params();
+  const vector<t_pppParam*>& params  = _parlist->params();
 
   for (unsigned ii = 0; ii < LCs.size(); ii++) {
@@ -294 +364 @@
         // ---------------------------------
         bool slip = false;
-
         if (obs->slip()) {
           LOG << epoTimeStr  << "cycle slip set (obs) " << obs->prn().toString()  << endl;
@@ -312 +381 @@
           slip = true;
         }
-        _slips[obs->prn()]._biasJumpCounter = obs->biasJumpCounter();
 
         // Slip Set
         // --------
         if (slip) {
-          resetAmb(obs->prn(), obsVector);
-        }
-
+          if (preProcessing) {
+            if (obs->prn() == refPrn) {
+              _obsPool->setEpoReProcessing(true);
+            }
+          }
+          else {
+            resetAmb(obs->prn(), obsVector);
+          }
+        }
         // Check Pre-Fit Residuals
         // -----------------------
@@ -328 +402 @@
             AA[iPar] = par->partial(_epoTime, obs, tLC);
           }
-
           double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA);
           double vv = DotProduct(AA, _xFlt) - ll;
-
           if (fabs(vv) > SLIP) {
-            LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
-                << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
-            resetAmb(obs->prn(), obsVector);
+            if (preProcessing) {
+              if (obs->prn() == refPrn) {
+                _obsPool->setEpoReProcessing(true);
+              }
+            }
+            else {
+              LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
+                  << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
+              resetAmb(obs->prn(), obsVector);
+            }
           }
         }