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

Timestamp:

Oct 28, 2026, 9:31:18 AM (37 hours ago)

Author:

stuerze

Message:

attitude check added to the SSR QC

File:

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

trunk/BNC/src/ssrQc.cpp

@@ -49 +49 @@
 #include "bncsettings.h"
 #include "combination/bnccomb.h"
+#include "bncantex.h"
 
 using namespace std;
@@ -60 +61 @@
   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 ATT_ONSET      = 20.0 * M_PI / 180.0; // [rad] change of the yaw difference starting a turn
+  const double ATT_STABLE     =  2.0 * M_PI / 180.0; // [rad] yaw difference stable again ...
+  const double ATT_STABLE_T   = 600.0;   // [s] ... for this time: end of the turn
+  const double ATT_REF_FROM   = 900.0;   // [s] reference before the turn: samples from 15 min
+  const double ATT_REF_TO     = 300.0;   //     to 5 min before the current epoch
+  const double ATT_MAX_GAP    = 600.0;   // [s] longer gaps interrupt a turn
+  const double ATT_MAX_TURN   = 10800.0; // [s] a turn is closed after 3 h at the latest
+  const double ATT_MIN_PRED   = 0.03;    // [m] smaller predicted clock changes are not assessed
   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
@@ -150 +159 @@
       delete _bsxRef;
       _bsxRef = 0;
+    }
+  }
+  _antex = 0;
+  QString antexFile = settings.value("ssrQcAntexFile").toString();
+  expandEnvVar(antexFile);
+  if (!antexFile.isEmpty()) {
+    _antex = new bncAntex();
+    if (_antex->readFile(antexFile) != success || _antex->numSatPrns() == 0) {
+      emit newMessage(("SSR QC: wrong ANTEX file " + antexFile + " - no attitude check").toLatin1(), true);
+      delete _antex;
+      _antex = 0;
     }
   }
@@ -247 +267 @@
   }
   delete _bsxRef;
+  delete _antex;
 }
 
@@ -488 +509 @@
     if (yawProvided) {
       stat.numYaw++;
+      t_yawLast& yl = _yawLast[QString::fromStdString(bias._staID)][QString::fromStdString(bias._prn.toString())];
+      double dt = 0.0;
+      if (bias._updateInt < unsigned(ssrUpdateInt.size())) {
+        dt = 0.5 * ssrUpdateInt[bias._updateInt];   // reference time in the middle of the interval
+      }
+      yl.time = bias._time + dt;
+      yl.yaw  = bias._yaw;
+      yl.rate = bias._yawRate;
     }
   }
@@ -986 +1015 @@
         if (dClk.contains(prn)) {
           double dc = dClk[prn] - off;
+          if (_antex && _yawLast.value(staID).contains(prn)) {
+            attitudeSample(staID, prn, tt, crd[staID][prn].xx, crd[staID][prn].vv, dc);
+          }
           stat.numClk += 1;
           stat.sClk   += dc * dc;
@@ -999 +1031 @@
     }
   }
+}
+
+// Attitude check: one sample per comparison epoch for a stream sending the
+// yaw angle. The yaw difference of BNC's attitude model ('Computed') minus
+// the stream is followed continuously over time; the predicted clock
+// difference of the stream to an AC following the model is
+// lambda_IF * dpsi / 2pi (see combination, dC_att). A turn starts when the
+// yaw difference changes by more than ATT_ONSET against the samples 5-15
+// minutes before, and ends when it is stable again for ATT_STABLE_T.
+////////////////////////////////////////////////////////////////////////////
+void t_ssrQc::attitudeSample(const QString& staID, const QString& prn, const bncTime& tt,
+                             const ColumnVector& xx, const ColumnVector& vv, double dc) {
+  const t_yawLast& yl = _yawLast[staID][prn];
+  double age = tt - yl.time;
+  if (!yl.time.valid() || fabs(age) > 120.0) {
+    return;
+  }
+  double psiAC = yl.yaw + yl.rate * age;
+
+  // model yaw, once per satellite and epoch (model state in time order)
+  double tk = tKey(tt);
+  if (tk != _psiModelT) {
+    _psiModel.clear();
+    _psiModelT = tk;
+  }
+  if (!_psiModel.contains(prn)) {
+    ColumnVector dx(3);
+    double psiM = 0.0;
+    double Mjd = tt.mjd() + tt.daysec() / 86400.0;
+    if (_antex->satCoMcorrection(prn, Mjd, xx, vv, dx, bncAntex::ATT_COMPUTED, 0.0,
+                                 "ssrQc", &psiM) != success) {
+      return;
+    }
+    _psiModel[prn] = psiM;
+  }
+
+  // ionosphere-free wavelength of the reference signals, c / (f1 + f2)
+  QString sig1, sig2;
+  t_frequency::type frq1, frq2;
+  if (!refSignals(prn[0].toLatin1(), sig1, sig2, frq1, frq2)) {
+    return;
+  }
+  double lambda = t_CST::c / (t_CST::freq(frq1, 0) + t_CST::freq(frq2, 0));
+
+  t_attTrack& tr = _attTrack[staID][prn];
+  if (tr.hasLast && tk - tr.lastT > ATT_MAX_GAP) {
+    if (tr.active) {
+      attitudeEvent(staID, prn, tr, "interrupted by a gap");
+    }
+    tr = t_attTrack();
+  }
+  double dpsi = remainder(_psiModel[prn] - psiAC, 2.0 * M_PI);
+  if (tr.hasLast) {
+    dpsi = tr.lastDPsi + remainder(dpsi - tr.lastDPsi, 2.0 * M_PI);
+  }
+  tr.hasLast  = true;
+  tr.lastDPsi = dpsi;
+  tr.lastT    = tk;
+  t_attSample smp;
+  smp.t    = tk;
+  smp.dc   = dc;
+  smp.p    = lambda * dpsi / (2.0 * M_PI);
+  smp.dpsi = dpsi;
+
+  if (!tr.active) {
+    tr.hist << smp;
+    while (!tr.hist.isEmpty() && tk - tr.hist.first().t > ATT_REF_FROM) {
+      tr.hist.removeFirst();
+    }
+    // reference: samples 5 to 15 minutes before
+    QVector<double> rDc, rP, rPsi;
+    for (const t_attSample& hh : tr.hist) {
+      if (tk - hh.t >= ATT_REF_TO) {
+        rDc << hh.dc; rP << hh.p; rPsi << hh.dpsi;
+      }
+    }
+    if (rPsi.size() >= 3 && fabs(dpsi - median(rPsi)) > ATT_ONSET) {
+      tr.active = true;
+      tr.dc0    = median(rDc);
+      tr.p0     = median(rP);
+      tr.dpsi0  = median(rPsi);
+      for (const t_attSample& hh : tr.hist) {
+        if (tk - hh.t < ATT_REF_TO) tr.event << hh;  // the beginning of the turn
+      }
+      tr.hist.clear();
+    }
+    return;
+  }
+
+  tr.event << smp;
+  // end: yaw difference stable for ATT_STABLE_T, or the turn too long
+  double tStart = tr.event.first().t;
+  bool   stable = (tk - tStart >= ATT_STABLE_T);
+  if (stable) {
+    for (int ii = tr.event.size() - 1; ii >= 0 && tk - tr.event[ii].t <= ATT_STABLE_T; ii--) {
+      if (fabs(tr.event[ii].dpsi - dpsi) > ATT_STABLE) {
+        stable = false;
+        break;
+      }
+    }
+  }
+  if (stable || tk - tStart > ATT_MAX_TURN) {
+    attitudeEvent(staID, prn, tr, stable ? "" : "closed after 3 h");
+    double keepPsi = tr.lastDPsi, keepT = tr.lastT;
+    tr = t_attTrack();
+    tr.hasLast  = true;
+    tr.lastDPsi = keepPsi;
+    tr.lastT    = keepT;
+    tr.hist << smp;
+  }
+}
+
+// Assessment of a turn: observed clock change of the stream relative to the
+// consensus against the predicted one (regression slope, RMS reduction)
+////////////////////////////////////////////////////////////////////////////
+void t_ssrQc::attitudeEvent(const QString& staID, const QString& prn, t_attTrack& tr,
+                            const QString& end) {
+  if (tr.event.size() < 3) {
+    return;
+  }
+  double sxx = 0.0, sxy = 0.0, syy = 0.0, srr = 0.0, maxX = 0.0, maxPsi = 0.0;
+  for (const t_attSample& ss : tr.event) {
+    double x = ss.p  - tr.p0;
+    double y = ss.dc - tr.dc0;
+    sxx += x * x; sxy += x * y; syy += y * y; srr += (y - x) * (y - x);
+    if (fabs(x) > fabs(maxX)) maxX = x;
+    if (fabs(ss.dpsi - tr.dpsi0) > fabs(maxPsi)) maxPsi = ss.dpsi - tr.dpsi0;
+  }
+  int    nn  = tr.event.size();
+  // changes at the end of the turn: mean of the last 5 minutes
+  double tEnd = tr.event.last().t, xEnd = 0.0, yEnd = 0.0;
+  int    nEnd = 0;
+  for (const t_attSample& ss : tr.event) {
+    if (tEnd - ss.t <= 300.0) {
+      xEnd += ss.p - tr.p0; yEnd += ss.dc - tr.dc0; ++nEnd;
+    }
+  }
+  xEnd /= nEnd; yEnd /= nEnd;
+  double rmsObs = sqrt(syy / nn);
+  double rmsRes = sqrt(srr / nn);
+  double slope  = (sxx > 0.0) ? sxy / sxx : 0.0;
+
+  bncTime t1, t2;
+  t1.set(int(tr.event.first().t / 604800.0), fmod(tr.event.first().t, 604800.0));
+  t2.set(int(tEnd / 604800.0), fmod(tEnd, 604800.0));
+  QString verdict;
+  if (fabs(maxX) < ATT_MIN_PRED) {
+    // the consensus is assumed to follow BNC's attitude model; a larger clock
+    // change shows that the other ACs model this turn differently
+    verdict = (rmsObs >= ATT_MIN_PRED) ? "clock change not explained by the yaw difference "
+                                         "(the other ACs model this turn differently?)"
+                                       : "predicted change too small to assess";
+  }
+  else if (slope > 0.7 && slope < 1.3 && rmsRes < 0.5 * rmsObs) {
+    verdict = "SSR yaw consistent with the clocks";
+  }
+  else {
+    verdict = "SSR yaw NOT consistent with the clocks";
+  }
+  QString line = QString("attitude check %1 %2-%3: yaw difference to the model %4 deg, clock change vs. "
+                         "consensus %5 cm, predicted %6 cm (max %7 cm), slope %8, clock RMS %9 -> %10 cm: %11")
+                 .arg(prn).arg(QString::fromStdString(t1.timestr(0))).arg(QString::fromStdString(t2.timestr(0)))
+                 .arg(maxPsi * 180.0 / M_PI, 0, 'f', 0)
+                 .arg(yEnd * 100.0, 0, 'f', 1).arg(xEnd * 100.0, 0, 'f', 1).arg(maxX * 100.0, 0, 'f', 1)
+                 .arg(slope, 0, 'f', 2).arg(rmsObs * 100.0, 0, 'f', 1).arg(rmsRes * 100.0, 0, 'f', 1)
+                 .arg(verdict);
+  if (!end.isEmpty()) {
+    line += " (" + end + ")";
+  }
+  _attEvents[staID][prn[0].toLatin1()] << line;
 }
 
@@ -1406 +1608 @@
         }
       }
+      for (const QString& line : _attEvents.value(staID).value(sys)) {
+        lines << line;
+      }
       addSats(QString("clock jumps > %1 m").arg(MAX_CLK_JUMP), jumps);
       if (!maxRadialPrn.isEmpty()) {
@@ -1489 +1694 @@
   _cmpStat.clear();
   _cmpEpochs.clear();
+  _attEvents.clear();
   for (QMap<QString, t_vtecStat>::iterator it = _vtecStat.begin(); it != _vtecStat.end(); ++it) {
     t_vtecStat fresh;