Changeset 11067 in ntrip

Timestamp:

Oct 5, 2026, 12:33:42 PM (23 hours ago)

Author:

stuerze

Message:

some more fixes

Location:

trunk/BNC/src

Files:

• PPP/pppClient.cpp (modified) (1 diff)
• PPP/pppSatObs.cpp (modified) (1 diff)
• combination/bnccomb.cpp (modified) (15 diffs)
• combination/bnccomb.h (modified) (3 diffs)

trunk/BNC/src/PPP/pppClient.cpp

@@ -61 +61 @@
   if (!_opt->_antexFileName.empty()) {
     _antex = new bncAntex(_opt->_antexFileName.c_str());
+    *_log << _antex->info().toStdString() << endl;
   }
   if (!_opt->_blqFileName.empty()) {

trunk/BNC/src/PPP/pppSatObs.cpp

@@ -672 +672 @@
       }
       else {
-        _model._antPCO[ii] += PPP_CLIENT->antex()->satCorr(prn, frqType, _model._elTx, _model._azTx, found);
+        _model._antPCO[ii] += PPP_CLIENT->antex()->satCorr(prn, frqType, _model._elTx, _model._azTx, found, _time);
+        if (!found && !PPP_CLIENT->antex()->hasSatEntry(prn, _time) &&
+            PPP_CLIENT->antex()->warnNoSatEntry(prn, _time)) {
+          LOG << "ANTEX: no satellite antenna entry valid at " << _time.datestr() << ' ' << _time.timestr()
+              << " for " << _prn.toString() << " - ANTEX file outdated?" << endl;
+        }
         if (OPT->_isAPC && found) {
           // the PCOs as given in the satellite antenna correction for all frequencies
           // have to be reduced by the PCO of the respective reference frequency
           if      (_prn.system() == 'G') {
-            _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::G1, _model._elTx, _model._azTx, found);
+            _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::G1, _model._elTx, _model._azTx, found, _time);
           }
           else if (_prn.system() == 'R') {
-            _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::R1, _model._elTx, _model._azTx, found);
+            _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::R1, _model._elTx, _model._azTx, found, _time);
           }
           else if (_prn.system() == 'E') {
-            _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::E1, _model._elTx, _model._azTx, found);
+            _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::E1, _model._elTx, _model._azTx, found, _time);
           }
           else if (_prn.system() == 'C') {
-            _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::C2, _model._elTx, _model._azTx, found);
+            _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::C2, _model._elTx, _model._azTx, found, _time);
           }
         }

trunk/BNC/src/combination/bnccomb.cpp

@@ -19 +19 @@
 #include <sstream>
 #include <map>
+#include <algorithm>
 
 #include "bnccomb.h"
@@ -290 +291 @@
       _antex = 0;
     }
+    else {
+      emit newMessage("bncComb: " + _antex->info().toLatin1(), true);
+    }
   }
 
@@ -853 +857 @@
     // Check satellite code biases
     // ----------------------------
+    QMap<t_frequency::type, double> codeBiasesRefSig;
     if (_satCodeBiases.contains(acName)) {
       QMap<t_prn, t_satCodeBias>& storage = _satCodeBiases[acName];
       if (storage.contains(clkCorr._prn)) {
         _newCorr->_satCodeBias = storage[clkCorr._prn];
-        QMap<t_frequency::type, double> codeBiasesRefSig;
         for (unsigned ii = 1; ii < cmbRefSig::cIF; ii++) {
           t_frequency::type frqType = cmbRefSig::toFreq(sys, static_cast<cmbRefSig::type>(ii));
@@ -870 +874 @@
           }
         }
-        // Code biases present, but not for both reference signals (e.g.
-        // Galileo HAS): the AC clocks cannot be referred to the reference
-        // signals and are used as they are - warn (at most once per hour
-        // per AC and system)
-        if (codeBiasesRefSig.size() < 2) {
-          QString key = acName + sys;
-          bncTime& warned = _refBiasWarned[key];
-          if (!warned.valid() || _newCorr->_time - warned >= 3600.0) {
-            warned = _newCorr->_time;
-            QString refSig;
-            for (unsigned ii = 1; ii < cmbRefSig::cIF; ii++) {
-              t_frequency::type frqType = cmbRefSig::toFreq(sys, static_cast<cmbRefSig::type>(ii));
-              refSig += QString(" C%1%2").arg(t_frequency::toString(frqType)[1])
-                                         .arg(cmbRefSig::toAttrib(sys, static_cast<cmbRefSig::type>(ii)));
-            }
-            emit newMessage("bncComb: " + acName.toLatin1() + " provides no code biases for the " +
-                            QByteArray(1, sys) + " reference signals" + refSig.toLatin1() +
-                            " - its clocks are combined without referring them to these signals", false);
+        // only with the biases of both reference signals (one of them alone
+        // does not refer the clocks to the reference signals)
+        if (codeBiasesRefSig.size() == 2) {
+          map<t_frequency::type, double> codeCoeff;
+          double channel = double(_newCorr->_eph->slotNum());
+          cmbRefSig::coeff(sys, cmbRefSig::cIF, channel, codeCoeff);
+          map<t_frequency::type, double>::const_iterator it;
+          for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
+            t_frequency::type frqType = it->first;
+            double codeCoeff          = it->second;
+            _newCorr->_satCodeBiasIF += codeCoeff * codeBiasesRefSig[frqType];
           }
-        }
-        map<t_frequency::type, double> codeCoeff;
-        double channel = double(_newCorr->_eph->slotNum());
-        cmbRefSig::coeff(sys, cmbRefSig::cIF, channel, codeCoeff);
-        map<t_frequency::type, double>::const_iterator it;
-        for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
-          t_frequency::type frqType = it->first;
-          double codeCoeff          = it->second;
-          _newCorr->_satCodeBiasIF += codeCoeff * codeBiasesRefSig[frqType];
-        }
-        if (codeBiasesRefSig.size() == 2) {
           _acBiasIF[acName][prnStr.mid(0,3)] = _newCorr->_satCodeBiasIF;
         }
@@ -906 +892 @@
         }
         _newCorr->_satCodeBias._bias.clear();
+      }
+    }
+    // Without the code biases of both reference signals (none at all, or
+    // only one of them, e.g. Galileo HAS for GPS), the AC clock cannot be
+    // referred to the reference signals: it is combined as it is, but
+    // excluded from the clock datum (see createAmat) - warn at most once
+    // per hour per AC and system
+    _newCorr->_refBiasesOK = (codeBiasesRefSig.size() == 2);
+    if (!_newCorr->_refBiasesOK) {
+      _acBiasIF[acName].remove(prnStr.mid(0,3));
+      QString key = acName + sys;
+      bncTime& warned = _refBiasWarned[key];
+      if (!warned.valid() || _newCorr->_time - warned >= 3600.0) {
+        warned = _newCorr->_time;
+        QString missing;
+        for (unsigned ii = 1; ii < cmbRefSig::cIF; ii++) {
+          t_frequency::type frqType = cmbRefSig::toFreq(sys, static_cast<cmbRefSig::type>(ii));
+          if (!codeBiasesRefSig.contains(frqType)) {
+            missing += QString(" C%1%2").arg(t_frequency::toString(frqType)[1])
+                                        .arg(cmbRefSig::toAttrib(sys, static_cast<cmbRefSig::type>(ii)));
+          }
+        }
+        emit newMessage("bncComb: " + acName.toLatin1() + " provides no code biases for the " +
+                        QByteArray(1, sys) + " reference signal(s)" + missing.toLatin1() +
+                        " (e.g. " + prnStr.mid(0,3).toLatin1() + ") - these clocks are combined"
+                        " without referring them to the reference signals and, as long as other ACs"
+                        " provide referred clocks, excluded from the clock datum", false);
       }
     }
@@ -1087 +1100 @@
     if (_method == filter) {
       checkBiasConsistency(epoTime, sys, out);
+      checkPersistentOutliers(epoTime, sys, out);
     }
     printResults(epoTime, out, masterCorr);
@@ -1121 +1135 @@
   if (checkOrbits(epoTime, sys, out, masterCorr) != success) {
     return failure;
+  }
+
+  // Statistics of persistent outliers
+  // ---------------------------------
+  for (int ii = 0; ii < corrs(sys).size(); ii++) {
+    const cmbCorr* corr = corrs(sys)[ii];
+    _outlierStat[sys][corr->_acName + ' ' + corr->_prn.mid(0,3)].numEpo += 1;
   }
 
@@ -1161 +1182 @@
 
       out << "  Outlier" << "\n";
+      _outlierStat[sys][corrs(sys)[maxResIndex-1]->_acName + ' ' +
+                        corrs(sys)[maxResIndex-1]->_prn.mid(0,3)].residuals.push_back(maxRes);
       _QQ[sys] = QQ_sav;
       delete corrs(sys)[maxResIndex-1];
@@ -1254 +1277 @@
       sxy += (xx[ii] - mx) * (yy[ii] - my);
     }
+    QString key = acName + sys;
     double spread = sqrt(sxx / nn);
-    if (spread < MIN_SPREAD || syy <= 0.0) {
-      continue;
-    }
-    double slope = sxy / sxx;
-    double corrCoeff = sxy / sqrt(sxx * syy);
-    if (slope < MIN_SLOPE || corrCoeff < MIN_CORR) {
+    double slope     = (sxx > 0.0) ? sxy / sxx : 0.0;
+    double corrCoeff = (sxx > 0.0 && syy > 0.0) ? sxy / sqrt(sxx * syy) : 0.0;
+    if (spread < MIN_SPREAD || slope < MIN_SLOPE || corrCoeff < MIN_CORR) {
+      if (_biasInconsistent.remove(key)) {
+        QString msg = QString("%1 %2 code biases consistent with its clocks again, "
+                              "AC used for the clock datum again").arg(acName).arg(sys);
+        out << epoTime.datestr().c_str() << " " << epoTime.timestr().c_str()
+            << " " << msg << "\n";
+        emit newMessage("bncComb: " + msg.toLatin1(), true);
+        resetSatOffsets(sys);
+      }
       continue;
     }
@@ -1275 +1304 @@
                   .arg(acName).arg(sys).arg(refSig).arg(nn)
                   .arg(spread, 0, 'f', 3).arg(slope, 0, 'f', 2).arg(corrCoeff, 0, 'f', 2);
+    if (!_biasInconsistent.contains(key)) {
+      _biasInconsistent.insert(key);
+      msg += ", AC excluded from the clock datum";
+      _biasCheckWarned.remove(key);
+      resetSatOffsets(sys);
+    }
     out << epoTime.datestr().c_str() << " " << epoTime.timestr().c_str()
         << " " << msg << "\n";
-    QString key = acName + sys;
     bncTime& warned = _biasCheckWarned[key];
     if (!warned.valid() || epoTime - warned >= 3600.0) {
       warned = epoTime;
       emit newMessage("bncComb: " + msg.toLatin1(), true);
+    }
+  }
+}
+
+// Report AC satellites rejected as outliers in most epochs
+////////////////////////////////////////////////////////////////////////////
+// A persistent clock difference of one AC for a satellite cannot be absorbed
+// by its satellite offset, as the offsets of all ACs of the satellite sum up
+// to zero: the correction is rejected again in every epoch. Once per hour
+// and system, AC satellites rejected in at least half of their epochs are
+// reported.
+////////////////////////////////////////////////////////////////////////////
+void bncComb::checkPersistentOutliers(bncTime epoTime, char sys, QTextStream& out) {
+
+  const double   INTERVAL = 3600.0;
+  const unsigned MIN_EPO  = 10;
+
+  if (!_outlierStatStart[sys].valid()) {
+    _outlierStatStart[sys] = epoTime;
+    return;
+  }
+  if (epoTime - _outlierStatStart[sys] < INTERVAL) {
+    return;
+  }
+  QMapIterator<QString, t_outlierStat> it(_outlierStat[sys]);
+  while (it.hasNext()) {
+    it.next();
+    const t_outlierStat& stat = it.value();
+    unsigned numOut = stat.residuals.size();
+    if (stat.numEpo < MIN_EPO || 2 * numOut < stat.numEpo) {
+      continue;
+    }
+    QVector<double> res = stat.residuals;
+    std::sort(res.begin(), res.end());
+    double median = (numOut % 2) ? res[numOut/2] : 0.5 * (res[numOut/2-1] + res[numOut/2]);
+    QString msg = QString("%1 rejected as outlier in %2% of %3 epochs in the last hour "
+                          "(median residual %4 m)")
+                  .arg(it.key()).arg(100.0 * numOut / stat.numEpo, 0, 'f', 0)
+                  .arg(stat.numEpo).arg(median, 0, 'f', 2);
+    out << epoTime.datestr().c_str() << " " << epoTime.timestr().c_str() << " " << msg << "\n";
+    emit newMessage("bncComb: " + msg.toLatin1(), true);
+  }
+  _outlierStat[sys].clear();
+  _outlierStatStart[sys] = epoTime;
+}
+
+// Re-initialize the satellite offsets of all ACs of a system (after a change
+// of the clock datum)
+////////////////////////////////////////////////////////////////////////////
+void bncComb::resetSatOffsets(char sys) {
+  for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
+    cmbParam* pp = _params[sys][iPar-1];
+    if (pp->type == cmbParam::offACSat) {
+      pp->xx = 0.0;
+      _QQ[sys].Row(iPar)    = 0.0;
+      _QQ[sys].Column(iPar) = 0.0;
+      _QQ[sys](iPar,iPar)   = pp->sig0 * pp->sig0;
     }
   }
@@ -1411 +1502 @@
       if (_antex->satCoMcorrection(corr->_prn, Mjd, xc.Rows(1,3), vv, dx,
                                    attMode, extYaw) != success) {
-        dx = 0;
-        emit newMessage("bncComb: antenna not found " + corr->_prn.mid(0,3).toLatin1(), false);
+        // without the satellite antenna offset, the orbit cannot be
+        // converted between antenna phase center and center of mass
+        if (_antex->warnNoSatEntry(corr->_prn, epoTime)) {
+          emit newMessage("bncComb: no satellite antenna entry valid at " +
+                          QByteArray(epoTime.datestr().c_str()) + " " + QByteArray(epoTime.timestr().c_str()) +
+                          " for " + corr->_prn.mid(0,3).toLatin1() +
+                          " in ANTEX file - ANTEX file outdated? Satellite not combined", true);
+        }
+        orbCorrections.pop_back();
+        clkCorrections.pop_back();
+        delete corr;
+        it.remove();
+        continue;
       }
       QString glonassYawLog = _antex->takeGlonassYawLog();
@@ -1537 +1639 @@
   int maxSat = _cmbSysPrn[sys];
 
-  const int nCon = (_method == filter) ? 1 + maxSat : 0;
+  // AC clocks not referred to the reference signals are excluded from the
+  // zero-sum conditions of the satellite offsets (the clock datum): those of
+  // ACs whose code biases are inconsistent with their clocks (see
+  // checkBiasConsistency) and those without the code biases of both
+  // reference signals. This applies as long as another AC provides referred
+  // corrections; the sum of the excluded offsets of an AC is constrained to
+  // zero instead.
+  QSet<QString> excludedOff;  // "AC PRN" (full PRN)
+  QStringList   excludedACs;  // ACs with excluded offsets
+  if (_method == filter) {
+    bool referredAC = false;
+    for (int ii = 0; ii < corrs(sys).size(); ii++) {
+      const cmbCorr* corr = corrs(sys)[ii];
+      if (_biasInconsistent.contains(corr->_acName + sys) || !corr->_refBiasesOK) {
+        excludedOff.insert(corr->_acName + ' ' + corr->_prn);
+        if (!excludedACs.contains(corr->_acName)) {
+          excludedACs << corr->_acName;
+        }
+      }
+      else {
+        referredAC = true;
+      }
+    }
+    if (!referredAC) {
+      excludedOff.clear();
+      excludedACs.clear();
+    }
+  }
+
+  const int nCon = (_method == filter) ? 1 + maxSat + excludedACs.size() : 0;
 
   AA.ReSize(nObs+nCon, nPar);  AA = 0.0;
@@ -1589 +1720 @@
    // }
     int iCond = 1;
+    QSet<int>     offExcluded; // offsets excluded from the zero-sum condition of their satellite
+    QSet<QString> acsInDatum;  // ACs with at least one offset in the zero-sum conditions
     // GNSS
     for (unsigned iGnss = 1; iGnss <= _cmbSysPrn[sys]; iGnss++) {
@@ -1595 +1728 @@
       ++iCond;
       PP(nObs+iCond) = Ph;
+      QList<int> offAll, offConsistent;
       for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
         cmbParam* pp = _params[sys][iPar-1];
@@ -1601 +1735 @@
              pp->type == cmbParam::offACSat                 &&
              pp->prn == prn) {
+          offAll << iPar;
+          if (!excludedOff.contains(pp->AC + ' ' + pp->prn)) {
+            offConsistent << iPar;
+          }
+        }
+      }
+      // satellite observed by inconsistent ACs only: all its offsets are used
+      const QList<int>& offUsed = offConsistent.isEmpty() ? offAll : offConsistent;
+      for (int iPar : offAll) {
+        if (offUsed.contains(iPar)) {
+          AA(nObs+iCond, iPar) = 1.0;
+          acsInDatum.insert(_params[sys][iPar-1]->AC);
+        }
+        else {
+          offExcluded << iPar;
+        }
+      }
+    }
+    // ACs with all their offsets excluded from the zero-sum conditions above:
+    // the sum of their offsets is constrained to zero (otherwise not
+    // separable from the AC's epoch offset)
+    for (const QString& acName : excludedACs) {
+      ++iCond;
+      if (acsInDatum.contains(acName)) {
+        continue;
+      }
+      PP(nObs+iCond) = Ph;
+      for (int iPar : offExcluded) {
+        if (_params[sys][iPar-1]->AC == acName) {
           AA(nObs+iCond, iPar) = 1.0;
         }

trunk/BNC/src/combination/bnccomb.h

@@ -117 +117 @@
       _dClkResult                  = 0.0;
       _satCodeBiasIF               = 0.0;
+      _refBiasesOK                 = false;
       _lambdaIF                    = 0.0;
       _satYawAngle                 = 0.0;
@@ -136 +137 @@
     QString        _acName;
     double         _satCodeBiasIF;
+    bool           _refBiasesOK;   // code biases of both reference signals available
     double         _lambdaIF;
     double         _satYawAngle;
@@ -274 +276 @@
   QMap<char, bncTime>                        _biasCheckLast;  // last bias/clock consistency check per system
   QMap<QString, bncTime>                     _biasCheckWarned; // last inconsistency warning per AC and system
+  QSet<QString>                              _biasInconsistent; // AC+system excluded from the clock datum
   void checkBiasConsistency(bncTime epoTime, char sys, QTextStream& out);
+  void checkPersistentOutliers(bncTime epoTime, char sys, QTextStream& out);
+  struct t_outlierStat {
+    t_outlierStat() : numEpo(0) {}
+    unsigned        numEpo;   // epochs with a correction of the AC for the satellite
+    QVector<double> residuals; // residuals of the epochs it was rejected
+  };
+  QMap<char, QMap<QString, t_outlierStat> > _outlierStat;      // per system, key "AC PRN"
+  QMap<char, bncTime>                       _outlierStatStart; // start of the statistics interval
+  void resetSatOffsets(char sys);
   bool allACsBeyond(const bncTime& epoTime) const;
   bncTime                                    _resTime;