Changeset 11117 in ntrip for trunk/BNC/src/ssrQc.cpp

Timestamp:

Oct 22, 2026, 5:22:15 PM (2 weeks ago)

Author:

stuerze

Message:

VTEC check added in ssr qc

File:

• trunk/BNC/src/ssrQc.cpp (modified) (10 diffs)

trunk/BNC/src/ssrQc.cpp

@@ -60 +60 @@
   const double DEV_MIN        = 0.05;  // [m] deviating satellites: RMS above max(DEV_MIN,
   const double DEV_FACTOR     = 3.0;   //     DEV_FACTOR x median RMS of the system)
+  const double VTEC_CMP_SAMPL = 300.0;  // [s] sampling of the VTEC comparison
+  const double VTEC_CMP_DELAY = 120.0;  // [s] VTEC epochs compared when data this much newer arrived
+  const double VTEC_HIST      = 3600.0; // [s] VTEC models kept for the comparison
+  const double VTEC_GRID_LAT  = 10.0;   // [deg] grid of the VTEC evaluation, latitudes -80 ... 80
+  const double VTEC_GRID_LON  = 15.0;   // [deg]
 
   double tKey(const bncTime& tt) {
     return tt.gpsw() * 604800.0 + tt.gpssec();
   }
-  const char*  TYPE_NAME[]    = {"ORBIT", "CLOCK", "CODE_BIAS", "PHASE_BIAS"};
+  const char*  TYPE_NAME[]    = {"ORBIT", "CLOCK", "CODE_BIAS", "PHASE_BIAS", "VTEC"};
 
   double median(QVector<double> vv) {
@@ -155 +160 @@
   connect(BNC_CORE, SIGNAL(newPhaseBiases(QList<t_satPhaseBias>)),
           this, SLOT(slotNewPhaseBiases(QList<t_satPhaseBias>)), Qt::DirectConnection);
+  connect(BNC_CORE, SIGNAL(newTec(t_vTec)),
+          this, SLOT(slotNewTec(t_vTec)), Qt::DirectConnection);
   // combined corrections (stream INTERNAL)
   if (bncComb::exists()) {
@@ -175 +182 @@
 ////////////////////////////////////////////////////////////////////////////
 QString t_ssrQc::acOf(const QString& staID) {
-  static const QRegularExpression re("^(?:SSR|OSB)[AC][0-9]{2}([A-Z0-9]{3})[0-9]$");
+  static const QRegularExpression re("^(?:SSR[AC]|OSB[AC]|IONO)[0-9]{2}([A-Z0-9]{3})[0-9]$");
   QRegularExpressionMatch mm = re.match(staID);
   return mm.hasMatch() ? mm.captured(1) : staID;
@@ -185 +192 @@
 ////////////////////////////////////////////////////////////////////////////
 QString t_ssrQc::solOf(const QString& staID) {
-  static const QRegularExpression re("^(?:SSR|OSB)[AC]([0-9]{2})([A-Z0-9]{3})[0-9]$");
+  static const QRegularExpression re("^(?:SSR[AC]|OSB[AC]|IONO)([0-9]{2})([A-Z0-9]{3})[0-9]$");
   QRegularExpressionMatch mm = re.match(staID);
   return mm.hasMatch() ? mm.captured(2) + mm.captured(1) : staID;
@@ -222 +229 @@
            << " completeness w.r.t. healthy broadcast ephemerides," << endl
            << "IOD matching, code and phase bias signals, yaw angle, code biases of the reference signals, clock jumps > "
-           << MAX_CLK_JUMP << " m, largest orbit corrections" << endl
+           << MAX_CLK_JUMP << " m, largest orbit corrections, VTEC models (system V)" << endl
            << "Reference signals of the clocks: " << t_ssrRefSig::currentString().toStdString() << endl;
     }
@@ -278 +285 @@
     compareEpoch(_nextCmpEpoch);
     _nextCmpEpoch = _nextCmpEpoch + CMP_SAMPL;
+  }
+  // VTEC comparison
+  if (!_vtecHist.isEmpty()) {
+    if (!_nextVtecCmp.valid()) {
+      double sec = ceil(time.gpssec() / VTEC_CMP_SAMPL) * VTEC_CMP_SAMPL;
+      _nextVtecCmp = bncTime(time.gpsw(), sec);
+    }
+    while (time - _nextVtecCmp >= VTEC_CMP_DELAY) {
+      compareVtec(_nextVtecCmp);
+      _nextVtecCmp = _nextVtecCmp + VTEC_CMP_SAMPL;
+    }
   }
   if (time - _intervalStart >= _interval) {
@@ -983 +1001 @@
 }
 
+// Ionospheric model (VTEC)
+////////////////////////////////////////////////////////////////////////////
+void t_ssrQc::slotNewTec(t_vTec vTec) {
+  QMutexLocker locker(&_mutex);
+  if (vTec._layers.empty()) {
+    return;
+  }
+  QString sta = QString::fromStdString(vTec._staID);
+  if (!_excluded.contains(sta) || sta == _reference) {
+    QMap<double, t_vTec>& hist = _vtecHist[sta];
+    hist[tKey(vTec._time)] = vTec;
+    while (!hist.isEmpty() && tKey(vTec._time) - hist.firstKey() > VTEC_HIST) {
+      hist.erase(hist.begin());
+    }
+  }
+  if (selected(vTec._staID)) {
+    newEpoch(vTec._staID, 'V', VTEC, vTec._time, vTec._updateInt, "VTEC");
+    QStringList layers;
+    for (const t_vTecLayer& layer : vTec._layers) {
+      layers << QString("%1 km (degree/order %2/%3)").arg(layer._height / 1000.0, 0, 'f', 0)
+                .arg(layer._C.Nrows() - 1).arg(layer._C.Ncols() - 1);
+    }
+    _vtecStat[sta].model = QString("%1 layer(s): %2").arg(vTec._layers.size()).arg(layers.join(", "));
+  }
+  checkInterval(vTec._time);
+}
+
+// VTEC [TECU] of a model at a geocentric latitude and longitude [rad] and
+// epoch: spherical harmonics in the sun-fixed frame as used in PPP, sum of
+// the layers (negative layer values set to zero as in PPP)
+////////////////////////////////////////////////////////////////////////////
+double t_ssrQc::vtecAt(const t_vTec& vTec, double lat, double lon, const bncTime& tt, bool& negative) {
+  double epoch = fmod(tt.gpssec(), 86400.0);
+  double lonS  = fmod(lon + (epoch - 50400.0) * M_PI / 43200.0, 2.0 * M_PI);
+  double vtec  = 0.0;
+  negative = false;
+  for (const t_vTecLayer& layer : vTec._layers) {
+    int    N  = layer._C.Nrows() - 1;
+    int    M  = layer._C.Ncols() - 1;
+    double vv = 0.0;
+    for (int n = 0; n <= N; n++) {
+      for (int m = 0; m <= min(n, M); m++) {
+        double pnm = associatedLegendreFunction(n, m, sin(lat));
+        double fac = (m == 0) ? sqrt(2.0 * n + 1.0)
+                              : sqrt(2.0 * (2.0 * n + 1.0) * factorial(n - m) / factorial(n + m));
+        vv += (layer._C[n][m] * cos(m * lonS) + layer._S[n][m] * sin(m * lonS)) * pnm * fac;
+      }
+    }
+    if (vv < 0.0) {
+      negative = true;
+      vv = 0.0;
+    }
+    vtec += vv;
+  }
+  return vtec;
+}
+
+// VTEC comparison: each model (latest one not older than twice its update
+// interval) evaluated on a global grid and compared with the median of the
+// other AC solutions (one vote per solution, own AC left out, combination
+// products not voting) or with the reference stream
+////////////////////////////////////////////////////////////////////////////
+void t_ssrQc::compareVtec(const bncTime& tt) {
+  double tk = tKey(tt);
+
+  QMap<QString, const t_vTec*> models;
+  QMapIterator<QString, QMap<double, t_vTec> > itH(_vtecHist);
+  while (itH.hasNext()) {
+    itH.next();
+    QMap<double, t_vTec>::const_iterator it = itH.value().upperBound(tk);
+    if (it == itH.value().constBegin()) {
+      continue;
+    }
+    --it;
+    const t_vTec& vTec = it.value();
+    double interval = (vTec._updateInt < unsigned(ssrUpdateInt.size()) && ssrUpdateInt[vTec._updateInt] > 0)
+                    ? ssrUpdateInt[vTec._updateInt] : VTEC_CMP_SAMPL;
+    double maxAge = min(3600.0, max(120.0, 2.0 * interval));
+    if (tk - it.key() <= maxAge) {
+      models[itH.key()] = &vTec;
+    }
+  }
+  if (models.isEmpty()) {
+    return;
+  }
+
+  // grid values of all models
+  QVector<QPair<double, double> > grid;
+  for (double lat = -80.0; lat <= 80.0 + 1e-6; lat += VTEC_GRID_LAT) {
+    for (double lon = -180.0; lon < 180.0 - 1e-6; lon += VTEC_GRID_LON) {
+      grid << QPair<double, double>(lat * M_PI / 180.0, lon * M_PI / 180.0);
+    }
+  }
+  QMap<QString, QVector<double> > val;
+  QMapIterator<QString, const t_vTec*> itM(models);
+  while (itM.hasNext()) {
+    itM.next();
+    QVector<double>& vv = val[itM.key()];
+    vv.resize(grid.size());
+    bool selectedSta = selected(itM.key().toStdString());
+    t_vtecStat* st = selectedSta ? &_vtecStat[itM.key()] : 0;
+    for (int ii = 0; ii < grid.size(); ii++) {
+      bool negative = false;
+      vv[ii] = vtecAt(*itM.value(), grid[ii].first, grid[ii].second, tt, negative);
+      if (st) {
+        if (st->numGrid == 0 || vv[ii] < st->minVal) st->minVal = vv[ii];
+        if (st->numGrid == 0 || vv[ii] > st->maxVal) st->maxVal = vv[ii];
+        st->numGrid += 1;
+        if (negative) st->numNeg += 1;
+      }
+    }
+  }
+
+  // comparison
+  for (const QString& staID : models.keys()) {
+    if (!selected(staID.toStdString())) {
+      continue;
+    }
+    QMap<QString, QStringList> perSol;   // solution: streams
+    for (const QString& oo : models.keys()) {
+      if (oo == staID) continue;
+      if (_reference == "CONSENSUS") {
+        if (!votes(oo) || acOf(oo) == acOf(staID)) continue;
+      }
+      else if (oo != _reference) {
+        continue;
+      }
+      perSol[solOf(oo)] << oo;
+    }
+    if (perSol.isEmpty()) {
+      continue;
+    }
+    t_vtecStat& st = _vtecStat[staID];
+    st.numCmp += 1;
+    st.numSol += perSol.size();
+    for (int ii = 0; ii < grid.size(); ii++) {
+      QVector<double> solVal;
+      for (const QStringList& streams : perSol) {
+        QVector<double> vv;
+        for (const QString& oo : streams) vv << val[oo][ii];
+        solVal << median(vv);
+      }
+      double diff = val[staID][ii] - median(solVal);
+      st.numDiff += 1;
+      st.sDiff   += diff;
+      st.sDiff2  += diff * diff;
+      st.maxAbs   = max(st.maxAbs, fabs(diff));
+    }
+  }
+}
+
+
 // Write the report of the current interval and reset the statistics
 ////////////////////////////////////////////////////////////////////////////
@@ -1029 +1199 @@
           << "  Latency med/max [s]  Sats" << endl;
     for (char sys : systems) {
-      for (int type = ORB; type <= PBIAS; type++) {
+      for (int type = ORB; type <= VTEC; type++) {
         QString kk = QString("%1 %2 %3").arg(staID).arg(sys).arg(type);
         if (!_epoStat.contains(kk)) {
@@ -1059 +1229 @@
     // Per system: completeness, ephemeris, biases, plausibility
     for (char sys : systems) {
+      // Ionosphere (VTEC): model, values on the grid, comparison
+      if (sys == 'V') {
+        QStringList lines;
+        QString kk = QString("%1 %2 %3").arg(staID).arg(sys).arg(int(VTEC));
+        if (_epoStat.contains(kk)) {
+          const t_epoStat& stat = _epoStat[kk];
+          double medInt = median(stat.intervals);
+          if (stat.nominal > 0 && stat.intervals.size() > 1 &&
+              (medInt > 1.5 * stat.nominal || medInt < stat.nominal / 1.5)) {
+            lines << QString("VTEC: declared update interval %1 s, actual %2 s")
+                     .arg(stat.nominal).arg(medInt, 0, 'f', 1);
+          }
+        }
+        if (_vtecStat.contains(staID)) {
+          const t_vtecStat& vs = _vtecStat[staID];
+          if (!vs.model.isEmpty()) {
+            lines << "VTEC model: " + vs.model;
+          }
+          if (vs.numGrid) {
+            lines << QString("VTEC on a global %1 x %2 deg grid: %3 ... %4 TECU, negative (set to zero) %5 %")
+                     .arg(VTEC_GRID_LAT, 0, 'f', 0).arg(VTEC_GRID_LON, 0, 'f', 0)
+                     .arg(vs.minVal, 0, 'f', 1).arg(vs.maxVal, 0, 'f', 1)
+                     .arg(100.0 * vs.numNeg / vs.numGrid, 0, 'f', 1);
+          }
+          if (vs.numDiff) {
+            QString refText = (_reference == "CONSENSUS")
+                            ? QString("consensus of the other ACs (%1 solutions on average)")
+                              .arg(double(vs.numSol) / vs.numCmp, 0, 'f', 1)
+                            : QString("stream %1").arg(_reference);
+            double mean = vs.sDiff / vs.numDiff;
+            double rms  = sqrt(vs.sDiff2 / vs.numDiff);
+            lines << QString("VTEC comparison with %1, %2 epochs: mean %3 TECU, RMS %4 TECU, max |diff| %5 TECU")
+                     .arg(refText).arg(vs.numCmp).arg(mean, 0, 'f', 2).arg(rms, 0, 'f', 2)
+                     .arg(vs.maxAbs, 0, 'f', 1);
+          }
+          else {
+            lines << "VTEC not compared: no other model available";
+          }
+        }
+        for (const QString& line : lines) {
+          out << "  " << sys << " " << line.toStdString() << endl;
+        }
+        continue;
+      }
       QMap<QString, unsigned> missOrb, missClk, missCB, missPB, iodUnknown, unhealthy, noEph,
                               jumps, missRef, noYaw, yawPartly;
@@ -1266 +1480 @@
   _cmpStat.clear();
   _cmpEpochs.clear();
+  for (QMap<QString, t_vtecStat>::iterator it = _vtecStat.begin(); it != _vtecStat.end(); ++it) {
+    t_vtecStat fresh;
+    fresh.model = it.value().model;
+    it.value() = fresh;
+  }
 
   // reset, keeping the state needed across intervals