Changeset 11129 in ntrip

Timestamp:

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

Author:

stuerze

Message:

attitude check added to the SSR QC

Location:

trunk/BNC/src

Files:

• bncsettings.cpp (modified) (1 diff)
• bncwindow.cpp (modified) (6 diffs)
• bncwindow.h (modified) (1 diff)
• ssrQc.cpp (modified) (9 diffs)
• ssrQc.h (modified) (2 diffs)

trunk/BNC/src/bncsettings.cpp

@@ -144 +144 @@
         setValue_p("ssrQcIntr", "1 hour");
         setValue_p("ssrQcBiasFile", "");
+        setValue_p("ssrQcAntexFile", "");
         setValue_p("ssrRefSignals", "");
         // Feed Engine

trunk/BNC/src/bncwindow.cpp

@@ -303 +303 @@
   _ssrQcBiasFileChooser = new qtFileChooser(0, qtFileChooser::File);
   _ssrQcBiasFileChooser->setFileName(settings.value("ssrQcBiasFile").toString());
+  _ssrQcAntexFileChooser = new qtFileChooser(0, qtFileChooser::File);
+  _ssrQcAntexFileChooser->setFileName(settings.value("ssrQcAntexFile").toString());
   _ssrQcIntrComboBox = new QComboBox();
   _ssrQcIntrComboBox->setEditable(false);
@@ -933 +935 @@
   cLayout->addWidget(new QLabel("QC Bias File"), 8, 0);
   cLayout->addWidget(_ssrQcBiasFileChooser, 8, 1, 1, 30);
-  cLayout->addWidget(new QLabel(""), 9, 1);
-  cLayout->setRowStretch(8, 999);
+  cLayout->addWidget(new QLabel("QC ANTEX File"), 9, 0);
+  cLayout->addWidget(_ssrQcAntexFileChooser, 9, 1, 1, 30);
+  cLayout->addWidget(new QLabel(""), 10, 1);
+  cLayout->setRowStretch(9, 999);
   cgroup->setLayout(cLayout);
 
@@ -1606 +1610 @@
   _ssrQcRefLineEdit->setWhatsThis(tr("<p>Specify the reference for the comparison of the orbits and clocks of the checked streams: 'CONSENSUS' for the median of all other streams, or the name of a stream, e.g. 'INTERNAL' for the combination.</p><p>Default is 'CONSENSUS'. <i>[key: ssrQcRef]</i></p>"));
   _ssrQcStreamsLineEdit->setWhatsThis(tr("<p>Specify the Broadcast Correction streams to be checked, as a comma separated list of stream names (e.g. SSRA00CNE1,SSRA00GMV0), or 'ALL' for all decoded streams. A stream with a leading '!' is excluded from the QC and from the consensus, e.g. 'ALL,!SIRGAS200001' for streams in another reference frame.</p><p>Default is 'ALL'. <i>[key: ssrQcStreams]</i></p>"));
+  _ssrQcAntexFileChooser->setWhatsThis(tr("<p>Specify an ANTEX file (e.g. the current igs20.atx) to enable the attitude check of the SSR QC: for streams sending the satellite yaw angle, the clock change of the stream relative to the other ACs during noon and midnight turns is compared with the change predicted from the difference of its yaw angle to BNC's attitude model ('Computed', see section 'Combination'). The ANTEX file provides the satellite block types the model needs.</p><p>Default is an empty option field, meaning no attitude check. <i>[key: ssrQcAntexFile]</i></p>"));
   _ssrQcBiasFileChooser->setWhatsThis(tr("<p>Specify a Bias SINEX file (OSB or DSB, e.g. a CAS DCB or a GFZ OSB product) as reference for the code bias check of the SSR QC: the DCBs of each stream are compared with the DCBs of this file, in addition to the comparison with the other ACs. A sign opposite to the reference for most DCBs indicates a stream using the SINEX instead of the SSR sign convention.</p><p>Default is an empty option field, meaning the code biases are compared with the other ACs only. <i>[key: ssrQcBiasFile]</i></p>"));
   _ssrQcIntrComboBox->setWhatsThis(tr("<p>Select the interval of the SSR QC reports. <i>[key: ssrQcIntr]</i></p>"));
@@ -1887 +1892 @@
   delete _ssrQcIntrComboBox;
   delete _ssrQcBiasFileChooser;
+  delete _ssrQcAntexFileChooser;
   delete _corrPortLineEdit;
   delete _outPortLineEdit;
@@ -2419 +2425 @@
   settings.setValue("ssrQcIntr", _ssrQcIntrComboBox->currentText());
   settings.setValue("ssrQcBiasFile", _ssrQcBiasFileChooser->fileName());
+  settings.setValue("ssrQcAntexFile", _ssrQcAntexFileChooser->fileName());
   settings.setValue("corrPort", _corrPortLineEdit->text());
   // Feed Engine
@@ -2959 +2966 @@
     enableWidget(!_ssrQcPathLineEdit->text().isEmpty(), _ssrQcRefLineEdit);
     enableWidget(!_ssrQcPathLineEdit->text().isEmpty(), _ssrQcBiasFileChooser);
+    enableWidget(!_ssrQcPathLineEdit->text().isEmpty(), _ssrQcAntexFileChooser);
   }
 

trunk/BNC/src/bncwindow.h

@@ -155 +155 @@
     QLineEdit* _ssrQcRefLineEdit;
     qtFileChooser* _ssrQcBiasFileChooser;
+    qtFileChooser* _ssrQcAntexFileChooser;
     QLineEdit* _miscMountLineEdit;
     QLineEdit* _miscPortLineEdit;

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;

trunk/BNC/src/ssrQc.h

@@ -16 +16 @@
 ////////////////////////////////////////////////////////////////////////////
 class bncBiasSnx;
+class bncAntex;
 
 class t_ssrQc : public QObject {
@@ -171 +172 @@
   QMap<QString, t_vtecStat>            _vtecStat;  // stream
   bncTime                              _nextVtecCmp;
+
+  // Attitude check: for streams sending the yaw angle, the clock change of
+  // the stream relative to the consensus during a turn (yaw difference to
+  // BNC's attitude model changing by more than ATT_ONSET) is compared with
+  // the change predicted by the phase wind-up, lambda_IF * dpsi / 2pi
+  class t_yawLast {          // latest SSR yaw angle of a stream and satellite
+   public:
+    bncTime time;
+    double  yaw  = 0.0;      // [rad]
+    double  rate = 0.0;      // [rad/s]
+  };
+  class t_attSample {
+   public:
+    double t;                // [s] (tKey)
+    double dc;               // clock difference to the consensus [m]
+    double p;                // predicted clock difference, unwrapped [m]
+    double dpsi;             // yaw difference model - stream, unwrapped [rad]
+  };
+  class t_attTrack {         // per stream and satellite
+   public:
+    bool                 active   = false;
+    bool                 hasLast  = false;
+    double               lastDPsi = 0.0;
+    double               lastT    = 0.0;
+    QVector<t_attSample> hist;   // recent samples before a turn
+    QVector<t_attSample> event;  // samples of the current turn
+    double               dc0 = 0.0, p0 = 0.0, dpsi0 = 0.0;  // before the turn
+  };
+  void   attitudeSample(const QString& staID, const QString& prn, const bncTime& tt,
+                        const ColumnVector& xx, const ColumnVector& vv, double dc);
+  void   attitudeEvent(const QString& staID, const QString& prn, t_attTrack& tr, const QString& end);
+  bncAntex*                                       _antex;      // option ssrQcAntexFile
+  QMap<QString, QMap<QString, t_yawLast> >        _yawLast;    // stream, PRN
+  QMap<QString, QMap<QString, t_attTrack> >       _attTrack;   // stream, PRN
+  QMap<QString, double>                           _psiModel;   // PRN: model yaw at _psiModelT
+  double                                          _psiModelT = 0.0;
+  QMap<QString, QMap<char, QStringList> >         _attEvents;  // stream, system: report lines
 };