Changeset 11117 in ntrip

Timestamp:

Oct 22, 2026, 5:22:15 PM (24 hours ago)

Author:

stuerze

Message:

VTEC check added in ssr qc

Location:

trunk/BNC

Files:

• CHANGELOG.md (modified) (2 diffs)
• src/bnchelp.html (modified) (1 diff)
• src/ssrQc.cpp (modified) (10 diffs)
• src/ssrQc.h (modified) (4 diffs)

trunk/BNC/CHANGELOG.md

@@ -4 +4 @@
 - ADDED: **SSR QC** of decoded Broadcast Correction streams (options 'QC Directory', 'QC Streams', 'QC Reference', 'QC Interval'; example configuration '31_SsrQC.bnc'). One daily report file per stream with epochs, declared and actual update intervals, gaps, latency, completeness w.r.t. healthy broadcast ephemerides, IOD matching, code and phase bias signals, yaw angle availability, missing code biases of the reference signals, clock jumps and largest orbit corrections
 - ADDED: Stream comparison in the SSR QC: orbits (radial, along, cross), clocks referred to the reference signals and signal-in-space differences (SIS = radial - clock) of each stream compared with the consensus of the other ACs (one vote per AC, own AC left out, combination products not voting) or with a selected stream (e.g. 'INTERNAL'), as RMS and standard deviation per satellite and system
+- ADDED: VTEC check in the SSR QC (pseudo system 'V'): update interval, gaps and latency of the VTEC messages, model (layers, heights, degree/order), VTEC range on a global grid and comparison with the consensus of the other AC solutions (or the 'QC Reference' stream) in TECU; ionosphere streams (IONO...) are assigned to their AC solution
 - ADDED: Code bias check in the SSR QC: the DCBs of each stream are compared with each other AC and optionally with a Bias SINEX file (option 'QC Bias File'); a stream with the opposite sign for most DCBs (SINEX instead of SSR sign convention) is reported
 - ADDED: SIS (radial - clock, satellite-specific clock offset included) per epoch and SISRMS/SISSig per satellite and system in the SP3 comparison
@@ -39 +40 @@
 - FIXED: GPS attitude model ('Computed'): at the end of a noon/midnight turn the yaw jumped to the nominal yaw (e.g. by 95 deg); the satellite now keeps yawing at its maximum rate until it has caught up with the nominal yaw
 - FIXED: GLONASS attitude model ('Computed'): the yaw was frozen while the nominal yaw rate exceeded the maximum rate and then jumped to the nominal yaw; now following Dilssner et al. (2011), the satellite yaws at the maximum rate (0.25 deg/s) during noon turns until it has caught up, and in the Earth's shadow until it reaches the nominal yaw at shadow exit
+- FIXED: Source table window: streams with an incomplete source table entry (e.g. without the fields authentication, fee, bitrate, misc) were not listed; they are now shown as long as the entry contains the fields up to 'nmea'
 - FIXED: Upload of raw streams lost data chunks
 - FIXED: RTCM3 decoder: a 0xD3 byte followed by a zero length field at the end of the received data caused an endless loop until std::bad_alloc, and the stream was deleted

trunk/BNC/src/bnchelp.html

@@ -5235 +5235 @@
     clock 0.689 (STD 0.008), radial 0.020, along 0.036, cross 0.022, SIS (radial - clock) 0.687 (STD 0.008)
   G satellites deviating from the reference, RMS SIS [m]: G02(0.38) G16(0.33) G20(0.91) G26(3.36)
+</pre>
+  <p>
+    <b>Ionospheric models (VTEC):</b> Streams with VTEC messages (e.g. the ionosphere streams 'IONO...', or SSR
+    streams also carrying VTEC) get the pseudo system 'V' (VTEC) in the report: epochs, declared and actual update interval,
+    gaps and latency of the VTEC messages, the model (number of layers, layer heights, degree and order of the
+    spherical harmonics), and the range of the VTEC values on a global grid (10 degrees in latitude from -80 to 80
+    degrees, 15 degrees in longitude), including the percentage of negative values, which are set to zero as in PPP.
+    Every 300 seconds (with a delay of 120 seconds), the latest model of each stream not older than twice its update
+    interval (at least 120 seconds, at most one hour) is evaluated on this grid and compared with the median of the
+    models of the other AC solutions (one vote per solution, own AC left out, combinations like 'IONO00IGS1' not
+    voting; at least one other solution), or with the selected 'QC Reference' stream. The mean difference (a bias
+    of the model) and the RMS are given in TEC units, e.g.
+  </p>
+  <pre>
+  V   VTEC            60        60 /  60.0 /   60.0     0             0       1.0 /    5.0      1
+  V VTEC model: 1 layer(s): 450 km (degree/order 15/15)
+  V VTEC on a global 10 x 15 deg grid: 2.1 ... 65.7 TECU, negative (set to zero) 0.0 %
+  V VTEC comparison with consensus of the other ACs (1.0 solutions on average), 12 epochs: mean 3.47 TECU, RMS 7.55 TECU, max |diff| 24.0 TECU
 </pre>
   <p>

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

trunk/BNC/src/ssrQc.h

@@ -10 +10 @@
 // Quality control of the decoded SSR streams: epochs and update intervals,
 // latency, completeness with respect to the broadcast ephemerides, IOD
-// matching, code bias signals, clock jumps and orbit correction size. A
+// matching, code bias signals, clock jumps, orbit correction size and the
+// ionospheric (VTEC) models. A
 // report per stream and system is written to the QC logfile at the end of
 // each interval and at the end of the program.
@@ -29 +30 @@
   void slotNewCodeBiases(QList<t_satCodeBias> codeBiases);
   void slotNewPhaseBiases(QList<t_satPhaseBias> phaseBiases);
+  void slotNewTec(t_vTec vTec);
 
  signals:
@@ -34 +36 @@
 
  private:
-  enum e_type {ORB = 0, CLK = 1, CBIAS = 2, PBIAS = 3};
+  enum e_type {ORB = 0, CLK = 1, CBIAS = 2, PBIAS = 3, VTEC = 4};
 
   class t_epoStat {          // per stream, system and message type
@@ -147 +149 @@
   bncBiasSnx*               _bsxRef;    // reference code biases (option ssrQcBiasFile)
   QString                   _bsxName;
+
+  // VTEC check: model content, values on a global grid and comparison with
+  // the consensus of the other ACs (or the reference stream)
+  class t_vtecStat {         // per stream
+   public:
+    unsigned numCmp  = 0;     // compared epochs
+    unsigned numSol  = 0;     // sum of the solutions in the consensus
+    unsigned numDiff = 0;     // compared grid values
+    double   sDiff   = 0.0;   // sum of the differences [TECU]
+    double   sDiff2  = 0.0;   // sum of squares [TECU^2]
+    double   maxAbs  = 0.0;   // [TECU]
+    unsigned numGrid = 0;     // evaluated grid values
+    unsigned numNeg  = 0;     // negative grid values
+    double   minVal  = 0.0;   // [TECU]
+    double   maxVal  = 0.0;
+    QString  model;           // layers, heights, degree/order of the latest model
+  };
+  void          compareVtec(const bncTime& tt);
+  static double vtecAt(const t_vTec& vTec, double lat, double lon, const bncTime& tt, bool& negative);
+  QMap<QString, QMap<double, t_vTec> > _vtecHist;  // stream, time
+  QMap<QString, t_vtecStat>            _vtecStat;  // stream
+  bncTime                              _nextVtecCmp;
 };