Index: /trunk/BNC/Example_Configs/31_SsrQC.bnc =================================================================== --- /trunk/BNC/Example_Configs/31_SsrQC.bnc (revision 11074) +++ /trunk/BNC/Example_Configs/31_SsrQC.bnc (revision 11075) @@ -23,5 +23,7 @@ corrPath= corrPort= -ssrQcFile=Output/SsrQC.txt +ssrQcPath=Output +ssrQcStreams=ALL +ssrQcRef=CONSENSUS ssrQcIntr=1 hour ephIntr=1 day Index: /trunk/BNC/src/bnccore.cpp =================================================================== --- /trunk/BNC/src/bnccore.cpp (revision 11074) +++ /trunk/BNC/src/bnccore.cpp (revision 11075) @@ -838,9 +838,13 @@ delete _ssrQc; _ssrQc = 0; bncSettings settings; - QString fileName = settings.value("ssrQcFile").toString(); - if (fileName.isEmpty()) { + QString dirName = settings.value("ssrQcPath").toString(); + if (dirName.isEmpty()) { return; } - expandEnvVar(fileName); + expandEnvVar(dirName); + QStringList streams = settings.value("ssrQcStreams").toString().split(QRegularExpression("[ ,;]+"), Qt::SkipEmptyParts); + if (streams.isEmpty()) { + streams << "ALL"; + } QStringList hlp = settings.value("ssrQcIntr").toString().split(' ', Qt::SkipEmptyParts); double interval = 3600.0; @@ -851,5 +855,5 @@ } } - _ssrQc = new t_ssrQc(fileName, interval); + _ssrQc = new t_ssrQc(dirName, interval, streams, settings.value("ssrQcRef").toString().trimmed()); } Index: /trunk/BNC/src/bnchelp.html =================================================================== --- /trunk/BNC/src/bnchelp.html (revision 11074) +++ /trunk/BNC/src/bnchelp.html (revision 11075) @@ -147,5 +147,5 @@             2.8.2 Interval
            2.8.3 Port
-             2.8.4 QC Logfile, QC Interval
+             2.8.4 QC Directory, QC Streams, QC Reference, QC Interval
            2.8.4 Wait for Full Corr Epoch
      2.9 Feed Engine
@@ -2245,6 +2245,6 @@ The configuration reads a Broadcast Ephemeris stream and three Broadcast Correction streams in IGS-SSR format (SSRA00CNE1, SSRA00CHC1) and RTCM-SSR format (Galileo HAS, SSRA00EUH0; needs an account - at the Galileo HAS caster). An hourly QC report per stream and system is saved in the QC logfile - 'SsrQC.txt'. The configuration can be applied to a raw file as well: + at the Galileo HAS caster). Hourly QC reports per system are saved in one file per stream in directory + 'Output'. The configuration can be applied to a raw file as well: bnc --conf 31_SsrQC.bnc --file <rawFile> --key ssrQcIntr "5 min" @@ -5066,11 +5066,15 @@

-

2.8.4 QC Logfile, QC Interval - optional

+

2.8.4 QC Directory, QC Streams, QC Reference, QC Interval - optional

BNC can check the quality of all Broadcast Correction streams it decodes, in real time or when reading a raw file - ('--file'). Specify the full path to a 'QC Logfile' to activate this function, and select the 'QC Interval' for - the reports (default: 1 hour). At the end of each interval and at the end of the program, a report is written per - stream and system. Broadcast ephemerides must be received as well, e.g. via stream 'BCEP00BKG0'. + ('--file'). Specify a 'QC Directory' to activate this function. 'QC Streams' selects the streams to be checked: + 'ALL' (default) for all decoded Broadcast Correction streams, or a comma separated list of stream names, e.g. + 'SSRA00CNE1,SSRA00GMV0'. 'QC Reference' selects the reference of the stream comparison (see below). Select + the 'QC Interval' for the reports (default: 1 hour). One daily file + '<stream>_S_<YYYYDDD>0000_01D_QC.txt' is written per stream; at the end of each interval and at the + end of the program, a report per system is appended. Broadcast ephemerides must be received as well, e.g. via + stream 'BCEP00BKG0'.

@@ -5115,5 +5119,5 @@ 

-Stream SSRA00GMV0
+=== SSR QC SSRA00GMV0 2026-10-02 11:40:00 - 2026-10-02 11:45:00 ===
   Sys Type        Epochs  Interval nom/med/max [s]  Gaps  Out-of-order  Latency med/max [s]  Sats
   G   ORBIT           47         1 /   5.0 /    5.0     0             0      11.0 /   12.0     32
@@ -5124,4 +5128,44 @@
   G ionosphere-free combination B_IF of the reference signal biases: std. dev. 0.586 m, max |B_IF| 2.058 m (32 satellites)
   G largest orbit corrections: radial 1.06 m (G11), along/cross 2.05 m (G11)
+
+

+ Stream comparison: In real time no ground truth is available. Therefore every 30 seconds (with a delay of 60 + seconds, so that the corrections of all streams have arrived), the orbits and clocks of each checked stream are + compared with a reference, selected with option 'QC Reference': +

+ +

+ For each stream the satellite positions and clocks are computed from the broadcast ephemeris and the latest + orbit and clock corrections at the comparison epoch. +

+ +

+ The report gives per system the RMS of the differences and, for clock and SIS, the standard deviation (STD) about + the mean of each satellite. Constant differences per satellite (a large RMS with a small STD) are typical between + ACs, as each AC realizes the satellite clocks and code biases with its own datum; a user applies the clocks and + biases of one AC consistently. Satellites whose RMS exceeds three times the median of the system (at least 5 cm) + are listed as deviating, e.g.: +

+ 
+  G comparison with consensus of the other streams (27 satellites), RMS (STD about the satellite means) [m]:
+    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)

Index: /trunk/BNC/src/bncmain.cpp =================================================================== --- /trunk/BNC/src/bncmain.cpp (revision 11074) +++ /trunk/BNC/src/bncmain.cpp (revision 11075) @@ -210,6 +210,8 @@ " corrIntr {File interval [character string: 1 min|2 min|5 min|10 min|15 min|30 min|1 hour|1 day]}\n" " corrPort {Output port [integer number]}\n" - " ssrQcFile {SSR QC logfile, full path [character string]}\n" - " ssrQcIntr {SSR QC report interval [character string: 5 min|15 min|30 min|1 hour|1 day]}\n" + " ssrQcPath {Directory for SSR QC files, one per stream [character string]}\n" + " ssrQcStreams {Streams to check [character string: ALL or comma separated list of stream names]}\n" + " ssrQcRef {Reference of the stream comparison [character string: CONSENSUS or a stream name, e.g. INTERNAL]}\n" + " ssrQcIntr {SSR QC report interval [character string: 5 min|15 min|30 min|1 hour|1 day]}\n" "\n" "Feed Engine Panel keys:\n" Index: /trunk/BNC/src/bncsettings.cpp =================================================================== --- /trunk/BNC/src/bncsettings.cpp (revision 11074) +++ /trunk/BNC/src/bncsettings.cpp (revision 11075) @@ -138,5 +138,7 @@ setValue_p("corrIntr", "1 day"); setValue_p("corrPort", ""); - setValue_p("ssrQcFile", ""); + setValue_p("ssrQcPath", ""); + setValue_p("ssrQcStreams", "ALL"); + setValue_p("ssrQcRef", "CONSENSUS"); setValue_p("ssrQcIntr", "1 hour"); // Feed Engine Index: /trunk/BNC/src/bncwindow.cpp =================================================================== --- /trunk/BNC/src/bncwindow.cpp (revision 11074) +++ /trunk/BNC/src/bncwindow.cpp (revision 11075) @@ -298,5 +298,7 @@ } _corrPortLineEdit = new QLineEdit(settings.value("corrPort").toString()); - _ssrQcFileLineEdit = new QLineEdit(settings.value("ssrQcFile").toString()); + _ssrQcPathLineEdit = new QLineEdit(settings.value("ssrQcPath").toString()); + _ssrQcRefLineEdit = new QLineEdit(settings.value("ssrQcRef").toString().isEmpty() ? QString("CONSENSUS") : settings.value("ssrQcRef").toString()); + _ssrQcStreamsLineEdit = new QLineEdit(settings.value("ssrQcStreams").toString().isEmpty() ? QString("ALL") : settings.value("ssrQcStreams").toString()); _ssrQcIntrComboBox = new QComboBox(); _ssrQcIntrComboBox->setEditable(false); @@ -309,5 +311,5 @@ _ssrQcIntrComboBox->setCurrentIndex(_ssrQcIntrComboBox->findText("1 hour")); } - connect(_ssrQcFileLineEdit, SIGNAL(textChanged(const QString&)), + connect(_ssrQcPathLineEdit, SIGNAL(textChanged(const QString&)), this, SLOT(slotBncTextChanged())); @@ -916,9 +918,13 @@ cLayout->addWidget(new QLabel("Port"), 3, 0); cLayout->addWidget(_corrPortLineEdit, 3, 1); - cLayout->addWidget(new QLabel("QC Logfile"), 4, 0); - cLayout->addWidget(_ssrQcFileLineEdit, 4, 1, 1, 30); - cLayout->addWidget(new QLabel("QC Interval"), 5, 0); - cLayout->addWidget(_ssrQcIntrComboBox, 5, 1); - cLayout->addWidget(new QLabel(""), 6, 1); + cLayout->addWidget(new QLabel("QC Directory"), 4, 0); + cLayout->addWidget(_ssrQcPathLineEdit, 4, 1, 1, 30); + cLayout->addWidget(new QLabel("QC Streams"), 5, 0); + cLayout->addWidget(_ssrQcStreamsLineEdit, 5, 1, 1, 30); + cLayout->addWidget(new QLabel("QC Reference"), 6, 0); + cLayout->addWidget(_ssrQcRefLineEdit, 6, 1, 1, 30); + cLayout->addWidget(new QLabel("QC Interval"), 7, 0); + cLayout->addWidget(_ssrQcIntrComboBox, 7, 1); + cLayout->addWidget(new QLabel(""), 8, 1); cLayout->setRowStretch(7, 999); cgroup->setLayout(cLayout); @@ -1588,5 +1594,7 @@ _corrPathLineEdit->setWhatsThis(tr("

Specify a directory for saving Broadcast Ephemeris Correction files.

If the specified directory does not exist, BNC will not create the files. [key: corrPath]

")); _corrIntrComboBox->setWhatsThis(tr("

Select the length of Broadcast Ephemeris Correction files. [key: corrIntr]

")); - _ssrQcFileLineEdit->setWhatsThis(tr("

Specify the full path to a logfile for the quality control of all decoded Broadcast Correction streams: epochs, update intervals, gaps, latency, completeness, IODs, code biases of the reference signals, clock jumps and orbit correction size.

Default is an empty option field, meaning that no SSR QC is done. [key: ssrQcFile]

")); + _ssrQcPathLineEdit->setWhatsThis(tr("

Specify a directory for the quality control of decoded Broadcast Correction streams: epochs, update intervals, gaps, latency, completeness, IODs, code and phase bias signals, code biases of the reference signals, clock jumps and orbit correction size. One daily file is written per stream.

Default is an empty option field, meaning that no SSR QC is done. [key: ssrQcPath]

")); + _ssrQcRefLineEdit->setWhatsThis(tr("

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.

Default is 'CONSENSUS'. [key: ssrQcRef]

")); + _ssrQcStreamsLineEdit->setWhatsThis(tr("

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.

Default is 'ALL'. [key: ssrQcStreams]

")); _ssrQcIntrComboBox->setWhatsThis(tr("

Select the interval of the SSR QC reports. [key: ssrQcIntr]

")); _corrPortLineEdit->setWhatsThis(tr("

BNC can produce Broadcast Ephemeris Corrections on your local host through an IP port.

Specify a port number here to activate this function. [key: corrPort]

")); @@ -1863,5 +1871,7 @@ delete _corrPathLineEdit; delete _corrIntrComboBox; - delete _ssrQcFileLineEdit; + delete _ssrQcPathLineEdit; + delete _ssrQcStreamsLineEdit; + delete _ssrQcRefLineEdit; delete _ssrQcIntrComboBox; delete _corrPortLineEdit; @@ -2391,5 +2401,7 @@ settings.setValue("corrPath", _corrPathLineEdit->text()); settings.setValue("corrIntr", _corrIntrComboBox->currentText()); - settings.setValue("ssrQcFile", _ssrQcFileLineEdit->text()); + settings.setValue("ssrQcPath", _ssrQcPathLineEdit->text()); + settings.setValue("ssrQcStreams", _ssrQcStreamsLineEdit->text()); + settings.setValue("ssrQcRef", _ssrQcRefLineEdit->text()); settings.setValue("ssrQcIntr", _ssrQcIntrComboBox->currentText()); settings.setValue("corrPort", _corrPortLineEdit->text()); @@ -2927,6 +2939,8 @@ enableWidget(enable, _corrIntrComboBox); } - if (sender() == 0 || sender() == _ssrQcFileLineEdit) { - enableWidget(!_ssrQcFileLineEdit->text().isEmpty(), _ssrQcIntrComboBox); + if (sender() == 0 || sender() == _ssrQcPathLineEdit) { + enableWidget(!_ssrQcPathLineEdit->text().isEmpty(), _ssrQcIntrComboBox); + enableWidget(!_ssrQcPathLineEdit->text().isEmpty(), _ssrQcStreamsLineEdit); + enableWidget(!_ssrQcPathLineEdit->text().isEmpty(), _ssrQcRefLineEdit); } Index: /trunk/BNC/src/bncwindow.h =================================================================== --- /trunk/BNC/src/bncwindow.h (revision 11074) +++ /trunk/BNC/src/bncwindow.h (revision 11075) @@ -151,5 +151,7 @@ QLineEdit* _ephPathLineEdit; QLineEdit* _corrPathLineEdit; - QLineEdit* _ssrQcFileLineEdit; + QLineEdit* _ssrQcPathLineEdit; + QLineEdit* _ssrQcStreamsLineEdit; + QLineEdit* _ssrQcRefLineEdit; QLineEdit* _miscMountLineEdit; QLineEdit* _miscPortLineEdit; Index: /trunk/BNC/src/combination/bnccomb.h =================================================================== --- /trunk/BNC/src/combination/bnccomb.h (revision 11074) +++ /trunk/BNC/src/combination/bnccomb.h (revision 11075) @@ -37,4 +37,5 @@ bncComb(const bncComb&) = delete; bncComb& operator=(const bncComb&) = delete; + static bool exists() {return instance != 0;} static void destruct() { delete instance; Index: /trunk/BNC/src/ssrQc.cpp =================================================================== --- /trunk/BNC/src/ssrQc.cpp (revision 11074) +++ /trunk/BNC/src/ssrQc.cpp (revision 11075) @@ -45,4 +45,6 @@ #include "bncutils.h" #include "bncconst.h" +#include "bncsettings.h" +#include "combination/bnccomb.h" using namespace std; @@ -51,4 +53,13 @@ const double MAX_CLK_JUMP = 0.5; // [m] clock change beyond its rate between consecutive epochs const double EPH_WARMUP = 120.0; // [s] no ephemeris checks at the start + const double CMP_SAMPL = 30.0; // [s] sampling of the stream comparison + const double CMP_DELAY = 60.0; // [s] epochs compared when data this much newer arrived + const double HIST_LENGTH = 300.0; // [s] corrections kept for the comparison + 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) + + double tKey(const bncTime& tt) { + return tt.gpsw() * 604800.0 + tt.gpssec(); + } const char* TYPE_NAME[] = {"ORBIT", "CLOCK", "CODE_BIAS", "PHASE_BIAS"}; @@ -91,8 +102,21 @@ // Constructor //////////////////////////////////////////////////////////////////////////// -t_ssrQc::t_ssrQc(const QString& fileName, double interval) : _ephUser(true) { - _fileName = fileName; +t_ssrQc::t_ssrQc(const QString& dirName, double interval, const QStringList& streams, + const QString& reference) : _ephUser(true) { + _dirName = dirName; + _reference = reference.isEmpty() ? QString("CONSENSUS") : reference; + bncSettings settings; + QStringList cmbStreams = settings.value("cmbStreams").toStringList(); + if (!cmbStreams.isEmpty()) { + _cmbMaster = cmbStreams[0].split(' ', Qt::SkipEmptyParts).value(0); + } _interval = interval > 0.0 ? interval : 3600.0; - _out = new ofstream(fileName.toLocal8Bit().data()); + for (const QString& sta : streams) { + if (sta.toUpper() == "ALL") { + _streams.clear(); + break; + } + _streams.insert(sta); + } connect(this, SIGNAL(newMessage(QByteArray,bool)), BNC_CORE, SLOT(slotMessage(const QByteArray,bool))); @@ -105,10 +129,42 @@ connect(BNC_CORE, SIGNAL(newPhaseBiases(QList)), this, SLOT(slotNewPhaseBiases(QList)), Qt::DirectConnection); - - *_out << "SSR QC of decoded correction streams, report interval " << _interval << " s" << endl - << "Checks: epochs and update intervals, gaps (interval > 2 x max(declared, median)), latency," - << " completeness w.r.t. healthy broadcast ephemerides," << endl - << "IOD matching, code biases of the reference signals, clock jumps > " << MAX_CLK_JUMP - << " m, largest orbit corrections" << endl; + // combined corrections (stream INTERNAL) + if (bncComb::exists()) { + connect(BNC_CMB, SIGNAL(newOrbCorrections(QList)), + this, SLOT(slotNewOrbCorrections(QList)), Qt::DirectConnection); + connect(BNC_CMB, SIGNAL(newClkCorrections(QList)), + this, SLOT(slotNewClkCorrections(QList)), Qt::DirectConnection); + } +} + +// Stream selected for the QC +//////////////////////////////////////////////////////////////////////////// +bool t_ssrQc::selected(const string& staID) const { + return _streams.isEmpty() || _streams.contains(QString::fromStdString(staID)); +} + +// Daily QC file of a stream: _S_0000_01D_QC.txt +//////////////////////////////////////////////////////////////////////////// +ofstream& t_ssrQc::streamOut(const QString& staID, const bncTime& time) { + unsigned year, month, day; + time.civil_date(year, month, day); + QDate date(year, month, day); + QString name = QDir(_dirName).filePath(QString("%1_S_%2%3%4_01D_QC.txt").arg(staID) + .arg(year).arg(date.dayOfYear(), 3, 10, QChar('0')).arg("0000")); + if (_outNames.value(staID) != name) { + delete _outs.value(staID, 0); + bool exists = QFile::exists(name); + ofstream* out = new ofstream(name.toLocal8Bit().data(), ios_base::out | ios_base::app); + if (!exists) { + *out << "SSR QC of stream " << staID.toStdString() << ", report interval " << _interval << " s" << endl + << "Checks: epochs and update intervals, gaps (interval > 2 x max(declared, median)), latency," + << " completeness w.r.t. healthy broadcast ephemerides," << endl + << "IOD matching, code and phase bias signals, code biases of the reference signals, clock jumps > " + << MAX_CLK_JUMP << " m, largest orbit corrections" << endl; + } + _outs[staID] = out; + _outNames[staID] = name; + } + return *_outs[staID]; } @@ -118,5 +174,7 @@ QMutexLocker locker(&_mutex); writeReport(true); - delete _out; + for (ofstream* out : _outs) { + delete out; + } } @@ -144,10 +202,16 @@ //////////////////////////////////////////////////////////////////////////// void t_ssrQc::checkInterval(const bncTime& time) { - if (!_startTime.valid()) { - _startTime = time; - } if (!_intervalStart.valid()) { double sec = floor(time.gpssec() / _interval) * _interval; _intervalStart = bncTime(time.gpsw(), sec); + } + // stream comparison of the epochs already complete for all streams + if (!_nextCmpEpoch.valid()) { + double sec = ceil(time.gpssec() / CMP_SAMPL) * CMP_SAMPL; + _nextCmpEpoch = bncTime(time.gpsw(), sec); + } + while (time - _nextCmpEpoch >= CMP_DELAY) { + compareEpoch(_nextCmpEpoch); + _nextCmpEpoch = _nextCmpEpoch + CMP_SAMPL; } if (time - _intervalStart >= _interval) { @@ -165,4 +229,11 @@ void t_ssrQc::newEpoch(const string& staID, char sys, e_type type, const bncTime& time, unsigned updateInt, const QString& prn) { + QString sta = QString::fromStdString(staID); + if (!_staStart.contains(sta)) { + _staStart[sta] = time; + } + if (!_staLast.contains(sta) || _staLast[sta] < time) { + _staLast[sta] = time; + } t_epoStat& stat = _epoStat[key(staID, sys, type)]; if (updateInt < unsigned(ssrUpdateInt.size())) { @@ -185,7 +256,8 @@ // Broadcast ephemeris the correction refers to //////////////////////////////////////////////////////////////////////////// -void t_ssrQc::checkEph(t_satStat& stat, const QString& prnInt, unsigned iod, const bncTime& time) { +void t_ssrQc::checkEph(t_satStat& stat, const QString& prnInt, unsigned iod, const bncTime& time, + const string& staID) { // not before all broadcast ephemerides may have been received - if (time - _startTime < EPH_WARMUP) { + if (time - _staStart[QString::fromStdString(staID)] < EPH_WARMUP) { return; } @@ -211,4 +283,8 @@ QMutexLocker locker(&_mutex); for (const t_orbCorr& corr : orbCorrections) { + if (selected(corr._staID) || corr._staID == _reference.toStdString()) { + storeHist(&corr, 0); + } + if (!selected(corr._staID)) continue; QString prn = QString::fromStdString(corr._prn.toString()); checkInterval(corr._time); @@ -216,5 +292,5 @@ t_satStat& stat = _satStat[QString::fromStdString(corr._staID) + ' ' + prn]; stat.numOrb += 1; - checkEph(stat, QString::fromStdString(corr._prn.toInternalString()), corr._iod, corr._time); + checkEph(stat, QString::fromStdString(corr._prn.toInternalString()), corr._iod, corr._time, corr._staID); double radial = fabs(corr._xr[0]); double alongCross = max(fabs(corr._xr[1]), fabs(corr._xr[2])); @@ -229,4 +305,8 @@ QMutexLocker locker(&_mutex); for (const t_clkCorr& corr : clkCorrections) { + if (selected(corr._staID) || corr._staID == _reference.toStdString()) { + storeHist(0, &corr); + } + if (!selected(corr._staID)) continue; QString prn = QString::fromStdString(corr._prn.toString()); checkInterval(corr._time); @@ -234,5 +314,5 @@ t_satStat& stat = _satStat[QString::fromStdString(corr._staID) + ' ' + prn]; stat.numClk += 1; - checkEph(stat, QString::fromStdString(corr._prn.toInternalString()), corr._iod, corr._time); + checkEph(stat, QString::fromStdString(corr._prn.toInternalString()), corr._iod, corr._time, corr._staID); // jump: change between consecutive epochs (same IOD) beyond the clock rate @@ -260,8 +340,12 @@ QMutexLocker locker(&_mutex); for (const t_satCodeBias& bias : codeBiases) { + bool sel = selected(bias._staID); + if (!sel && bias._staID != _reference.toStdString()) continue; char sys = bias._prn.system(); QString prn = QString::fromStdString(bias._prn.toString()); - checkInterval(bias._time); - newEpoch(bias._staID, sys, CBIAS, bias._time, bias._updateInt, prn); + if (sel) { + checkInterval(bias._time); + newEpoch(bias._staID, sys, CBIAS, bias._time, bias._updateInt, prn); + } QString staSys = QString("%1 %2").arg(bias._staID.c_str()).arg(sys); QMap values; @@ -302,4 +386,5 @@ QMutexLocker locker(&_mutex); for (const t_satPhaseBias& bias : phaseBiases) { + if (!selected(bias._staID)) continue; checkInterval(bias._time); newEpoch(bias._staID, bias._prn.system(), PBIAS, bias._time, bias._updateInt, @@ -312,14 +397,240 @@ } +// Reference point of the orbits of a stream: APC (SSRA...), CoM (SSRC...); +// the combined stream INTERNAL refers to that of its master AC (the first +// AC of the combination) +//////////////////////////////////////////////////////////////////////////// +char t_ssrQc::refPoint(const QString& staID) const { + QString sta = (staID == "INTERNAL") ? _cmbMaster : staID; + if (sta.startsWith("SSRA")) return 'A'; + if (sta.startsWith("SSRC")) return 'C'; + return '?'; +} + +// Keep the recent orbit and clock corrections of a stream +//////////////////////////////////////////////////////////////////////////// +void t_ssrQc::storeHist(const t_orbCorr* orb, const t_clkCorr* clk) { + if (orb) { + QMap& hist = _orbHist[QString::fromStdString(orb->_staID)] + [QString::fromStdString(orb->_prn.toInternalString())]; + hist[tKey(orb->_time)] = *orb; + while (!hist.isEmpty() && hist.firstKey() < tKey(orb->_time) - HIST_LENGTH) { + hist.erase(hist.begin()); + } + } + if (clk) { + QMap& hist = _clkHist[QString::fromStdString(clk->_staID)] + [QString::fromStdString(clk->_prn.toInternalString())]; + hist[tKey(clk->_time)] = *clk; + while (!hist.isEmpty() && hist.firstKey() < tKey(clk->_time) - HIST_LENGTH) { + hist.erase(hist.begin()); + } + } +} + +// Compare the orbits and clocks of all streams at epoch tt +//////////////////////////////////////////////////////////////////////////// +void t_ssrQc::compareEpoch(const bncTime& tt) { + double tk = tKey(tt); + + // Satellite positions and clocks of all streams (latest corrections <= tt) + QMap > crd; // stream, PRN (external) + QMapIterator > > itSta(_clkHist); + while (itSta.hasNext()) { + itSta.next(); + const QString& staID = itSta.key(); + QMapIterator > itPrn(itSta.value()); + while (itPrn.hasNext()) { + itPrn.next(); + const QString& prnInt = itPrn.key(); + QMap::const_iterator itC = itPrn.value().upperBound(tk); + if (itC == itPrn.value().constBegin()) { + continue; + } + const t_clkCorr& clk = (--itC).value(); + const QMap& orbHist = _orbHist[staID][prnInt]; + QMap::const_iterator itO = orbHist.upperBound(tk); + if (itO == orbHist.constBegin()) { + continue; + } + const t_orbCorr& orb = (--itO).value(); + t_eph* eph = 0; + t_eph* ephLast = _ephUser.ephLast(prnInt); + t_eph* ephPrev = _ephUser.ephPrev(prnInt); + if (ephLast && ephLast->IOD() == orb._iod) eph = ephLast; + else if (ephPrev && ephPrev->IOD() == orb._iod) eph = ephPrev; + if (!eph || clk._iod != orb._iod) { + continue; + } + eph->setOrbCorr(&orb); + eph->setClkCorr(&clk); + ColumnVector xc(6), vv(3); + if (eph->getCrd(tt, xc, vv, true) != success) { + continue; + } + t_satCrd& sc = crd[staID][QString::fromStdString(clk._prn.toString())]; + sc.xx = xc.Rows(1,3); + sc.vv = vv; + sc.clk = xc(4) * t_CST::c; + sc.refPoint = refPoint(staID); + // clock referred to the reference signals (as in the combination): + // + B_IF of the stream's code biases; the combination (INTERNAL) + // refers to them already + sc.clkOK = (staID == "INTERNAL"); + QString satKey = staID + ' ' + QString::fromStdString(clk._prn.toString()); + if (!sc.clkOK && _satStat.contains(satKey) && _satStat[satKey].refBiasOK) { + sc.clk += _satStat[satKey].biasIF; + sc.clkOK = true; + } + } + } + + // Satellites per system + QMap > satsSys; + QMapIterator > itCrd(crd); + while (itCrd.hasNext()) { + itCrd.next(); + for (const QString& prn : itCrd.value().keys()) { + satsSys[prn[0].toLatin1()].insert(prn); + } + } + + QMapIterator > itSys(satsSys); + while (itSys.hasNext()) { + itSys.next(); + const QSet& sats = itSys.value(); + QStringList staList; + for (const QString& staID : crd.keys()) { + for (const QString& prn : sats) { + if (crd[staID].contains(prn)) { + staList << staID; + break; + } + } + } + + // Clock datum: offset of each stream w.r.t. the per-satellite median of + // all streams (or w.r.t. the reference stream) + QMap offset; + for (const QString& staID : staList) { + QVector diffs; + for (const QString& prn : sats) { + if (!crd[staID].contains(prn) || !crd[staID][prn].clkOK) continue; + if (_reference == "CONSENSUS") { + QVector all; + for (const QString& oo : staList) { + if (oo != "INTERNAL" && crd[oo].contains(prn) && crd[oo][prn].clkOK) all << crd[oo][prn].clk; + } + if (all.size() >= 3) diffs << crd[staID][prn].clk - median(all); + } + else if (crd.value(_reference).contains(prn) && crd[_reference][prn].clkOK) { + diffs << crd[staID][prn].clk - crd[_reference][prn].clk; + } + } + if (!diffs.isEmpty()) { + offset[staID] = median(diffs); + } + } + + // Differences of each selected stream to its reference + for (const QString& staID : staList) { + if (!selected(staID.toStdString()) || staID == _reference) { + continue; + } + QMap dClk; + QMap dRac; + for (const QString& prn : sats) { + if (!crd[staID].contains(prn)) continue; + const t_satCrd& sc = crd[staID][prn]; + QStringList refs; + if (_reference == "CONSENSUS") { + for (const QString& oo : staList) { + // the combination (INTERNAL) is derived from the other streams + if (oo != staID && oo != "INTERNAL" && crd[oo].contains(prn)) refs << oo; + } + if (refs.size() < 2) continue; + refs.clear(); + for (const QString& oo : staList) { + if (oo != staID && oo != "INTERNAL" && crd[oo].contains(prn)) refs << oo; + } + } + else { + if (!crd.value(_reference).contains(prn)) continue; + refs << _reference; + } + QVector clks, xs[3]; + for (const QString& oo : refs) { + const t_satCrd& rc = crd[oo][prn]; + if (rc.clkOK && (oo == _reference || offset.contains(oo))) { + clks << rc.clk - (oo == _reference ? 0.0 : offset[oo]); + } + if (rc.refPoint == sc.refPoint && sc.refPoint != '?') { + for (int ii = 0; ii < 3; ii++) xs[ii] << rc.xx[ii]; + } + } + bool clkCmp = sc.clkOK && offset.contains(staID) && + (_reference == "CONSENSUS" ? clks.size() >= 2 : clks.size() == 1); + if (clkCmp) { + dClk[prn] = sc.clk - offset[staID] - median(clks); + } + else if (!sc.clkOK) { + _cmpStat[staID][prn].clkNotReferred += 1; + } + if (xs[0].size() == refs.size() || (xs[0].size() >= 2 && _reference == "CONSENSUS")) { + ColumnVector dx(3); + for (int ii = 0; ii < 3; ii++) dx[ii] = sc.xx[ii] - median(xs[ii]); + ColumnVector er = sc.xx / sc.xx.NormFrobenius(); + ColumnVector ec = crossproduct(sc.xx, sc.vv); + ec /= ec.NormFrobenius(); + ColumnVector ea = crossproduct(ec, er); + ColumnVector rac(3); + rac[0] = DotProduct(dx, er); + rac[1] = DotProduct(dx, ea); + rac[2] = DotProduct(dx, ec); + dRac[prn] = rac; + } + } + if (dClk.isEmpty() && dRac.isEmpty()) { + continue; + } + // remaining common clock offset of this stream and epoch + double off = dClk.isEmpty() ? 0.0 : median(dClk.values().toVector()); + _cmpEpochs[staID] += 1; + QSet prns; + for (const QString& prn : dClk.keys()) prns.insert(prn); + for (const QString& prn : dRac.keys()) prns.insert(prn); + for (const QString& prn : prns) { + t_cmpStat& stat = _cmpStat[staID][prn]; + if (dRac.contains(prn)) { + const ColumnVector& rac = dRac[prn]; + stat.numOrb += 1; + stat.sR += rac[0] * rac[0]; + stat.sA += rac[1] * rac[1]; + stat.sC += rac[2] * rac[2]; + } + if (dClk.contains(prn)) { + double dc = dClk[prn] - off; + stat.numClk += 1; + stat.sClk += dc * dc; + stat.mClk += dc; + if (dRac.contains(prn)) { + double sis = dRac[prn][0] - dc; + stat.numSis += 1; + stat.sSis += sis * sis; + stat.mSis += sis; + } + } + } + } + } +} + // Write the report of the current interval and reset the statistics //////////////////////////////////////////////////////////////////////////// void t_ssrQc::writeReport(bool final) { - if (!_out || _epoStat.isEmpty()) { + if (_epoStat.isEmpty()) { return; } bncTime endTime = final ? _lastTime : _intervalStart + _interval; - *_out << endl << "=== SSR QC " << _intervalStart.datestr() << ' ' << _intervalStart.timestr(0) - << " - " << endTime.datestr() << ' ' << endTime.timestr(0) - << (final ? " (end of data)" : "") << " ===" << endl; // Healthy and current broadcast ephemerides per system @@ -347,6 +658,13 @@ QList systems = itSta.value().values(); std::sort(systems.begin(), systems.end()); - *_out << endl << "Stream " << staID.toStdString() << endl - << " Sys Type Epochs Interval nom/med/max [s] Gaps Out-of-order" + ofstream& out = streamOut(staID, _intervalStart); + if (final) { + endTime = _staLast.value(staID, _lastTime); + } + out << endl << "=== SSR QC " << staID.toStdString() << ' ' + << _intervalStart.datestr() << ' ' << _intervalStart.timestr(0) + << " - " << endTime.datestr() << ' ' << endTime.timestr(0) + << (final ? " (end of data)" : "") << " ===" << endl + << " Sys Type Epochs Interval nom/med/max [s] Gaps Out-of-order" << " Latency med/max [s] Sats" << endl; for (char sys : systems) { @@ -369,5 +687,5 @@ double maxLat = stat.latencies.isEmpty() ? 0.0 : *std::max_element(stat.latencies.begin(), stat.latencies.end()); - *_out << " " << sys << " " << left << setw(10) << TYPE_NAME[type] << right + out << " " << sys << " " << left << setw(10) << TYPE_NAME[type] << right << setw(8) << stat.numEpo << setw(10) << stat.nominal << " /" << setw(6) << setprecision(1) << fixed << medInt @@ -484,15 +802,82 @@ .arg(maxRadial, 0, 'f', 2).arg(maxRadialPrn).arg(maxAlongCross, 0, 'f', 2).arg(maxAlongCrossPrn); } + // stream comparison + if (_cmpStat.contains(staID)) { + unsigned nSat = 0, nClk = 0, nOrb = 0, nSis = 0; + double sClk = 0.0, sR = 0.0, sA = 0.0, sC = 0.0, sSis = 0.0; + double vClk = 0.0, vSis = 0.0; // sums of squares about the satellite means + QMap rmsSis, rmsClk, rmsOrb; + QMap notReferred; + QMapIterator itC(_cmpStat[staID]); + while (itC.hasNext()) { + itC.next(); + if (itC.key()[0].toLatin1() != sys) continue; + const t_cmpStat& st = itC.value(); + if (st.numClk == 0 && st.numOrb == 0) { + if (st.clkNotReferred) notReferred[itC.key()] = 0; + continue; + } + ++nSat; + if (st.clkNotReferred) notReferred[itC.key()] = 0; + nClk += st.numClk; sClk += st.sClk; + nOrb += st.numOrb; sR += st.sR; sA += st.sA; sC += st.sC; + nSis += st.numSis; sSis += st.sSis; + if (st.numClk) vClk += st.sClk - st.mClk * st.mClk / st.numClk; + if (st.numSis) vSis += st.sSis - st.mSis * st.mSis / st.numSis; + if (st.numSis) rmsSis[itC.key()] = sqrt(st.sSis / st.numSis); + else if (st.numClk) rmsClk[itC.key()] = sqrt(st.sClk / st.numClk); + if (st.numOrb) rmsOrb[itC.key()] = sqrt((st.sR + st.sA + st.sC) / st.numOrb); + } + if (nSat) { + QString refText = (_reference == "CONSENSUS") ? QString("consensus of the other streams") + : QString("stream %1").arg(_reference); + QStringList parts; + if (nClk) parts << QString("clock %1 (STD %2)").arg(sqrt(sClk / nClk), 0, 'f', 3) + .arg(sqrt(max(0.0, vClk) / nClk), 0, 'f', 3); + if (nOrb) parts << QString("radial %1, along %2, cross %3") + .arg(sqrt(sR / nOrb), 0, 'f', 3).arg(sqrt(sA / nOrb), 0, 'f', 3) + .arg(sqrt(sC / nOrb), 0, 'f', 3); + if (nSis) parts << QString("SIS (radial - clock) %1 (STD %2)").arg(sqrt(sSis / nSis), 0, 'f', 3) + .arg(sqrt(max(0.0, vSis) / nSis), 0, 'f', 3); + lines << QString("comparison with %1 (%2 satellites), RMS (STD about the satellite means) [m]: %3") + .arg(refText).arg(nSat).arg(parts.join(", ")); + if (!nOrb) { + lines << "orbits not compared: no other stream with the same reference point (APC/CoM)"; + } + addSats("clocks not compared (code biases of the reference signals missing)", notReferred); + auto deviating = [&lines](const QString& what, const QMap& rms) { + if (rms.size() < 2) return; + double med = median(rms.values().toVector()); + QStringList dev; + QMapIterator itR(rms); + while (itR.hasNext()) { + itR.next(); + if (itR.value() > max(DEV_MIN, DEV_FACTOR * med)) { + dev << QString("%1(%2)").arg(itR.key()).arg(itR.value(), 0, 'f', 2); + } + } + if (!dev.isEmpty()) { + lines << QString("satellites deviating from the reference, RMS %1 [m]: %2").arg(what).arg(dev.join(" ")); + } + }; + deviating("SIS", rmsSis); + deviating("clock", rmsClk); + deviating("orbit (3D)", rmsOrb); + } + } for (const QString& line : lines) { - *_out << " " << sys << " " << line.toStdString() << endl; - } - } - } - _out->flush(); - - emit newMessage(QString("SSR QC: report %1 %2 written to %3") + out << " " << sys << " " << line.toStdString() << endl; + } + } + out.flush(); + } + + emit newMessage(QString("SSR QC: report %1 %2 of %3 stream(s) written to %4") .arg(QString::fromStdString(_intervalStart.datestr())) .arg(QString::fromStdString(_intervalStart.timestr(0))) - .arg(_fileName).toLatin1(), false); + .arg(streams.size()).arg(_dirName).toLatin1(), false); + + _cmpStat.clear(); + _cmpEpochs.clear(); // reset, keeping the state needed across intervals Index: /trunk/BNC/src/ssrQc.h =================================================================== --- /trunk/BNC/src/ssrQc.h (revision 11074) +++ /trunk/BNC/src/ssrQc.h (revision 11075) @@ -18,5 +18,6 @@ public: - t_ssrQc(const QString& fileName, double interval); + t_ssrQc(const QString& dirName, double interval, const QStringList& streams, + const QString& reference); ~t_ssrQc(); @@ -69,5 +70,6 @@ void newEpoch(const std::string& staID, char sys, e_type type, const bncTime& time, unsigned updateInt, const QString& prn); - void checkEph(t_satStat& stat, const QString& prnInt, unsigned iod, const bncTime& time); + void checkEph(t_satStat& stat, const QString& prnInt, unsigned iod, const bncTime& time, + const std::string& staID); void checkInterval(const bncTime& time); void writeReport(bool final); @@ -77,10 +79,49 @@ QMutex _mutex; - QString _fileName; - std::ofstream* _out; + bool selected(const std::string& staID) const; + + // Stream comparison: every CMP_SAMPL seconds, the orbits and clocks of each + // stream are compared with a reference: the consensus (median) of all + // other streams, or a selected stream (e.g. INTERNAL, the combination) + class t_cmpStat { // per stream and satellite + public: + t_cmpStat() : numClk(0), numOrb(0), numSis(0), clkNotReferred(0), + sClk(0.0), sR(0.0), sA(0.0), sC(0.0), sSis(0.0), mClk(0.0), mSis(0.0) {} + unsigned numClk; + unsigned numOrb; + unsigned numSis; + unsigned clkNotReferred; // epochs without reference signal biases + double sClk, sR, sA, sC, sSis; // sums of squares [m^2] + double mClk, mSis; // sums [m] + }; + class t_satCrd { // satellite position, velocity, clock of a stream + public: + ColumnVector xx; // [m] + ColumnVector vv; // [m/s] + double clk; // [m], referred to the reference signals (+ B_IF) if clkOK + bool clkOK; + char refPoint; // 'A' APC, 'C' CoM, '?' unknown + }; + void storeHist(const t_orbCorr* orb, const t_clkCorr* clk); + void compareEpoch(const bncTime& tt); + char refPoint(const QString& staID) const; + QString _reference; // "CONSENSUS" or stream name + QString _cmbMaster; // mountpoint of the combination's first AC + bncTime _nextCmpEpoch; + QMap > > _orbHist; // stream, PRN (internal), time + QMap > > _clkHist; + QMap > _cmpStat; // stream, PRN + QMap _cmpEpochs; // stream: compared epochs + std::ofstream& streamOut(const QString& staID, const bncTime& time); + + QString _dirName; + QSet _streams; // selected streams, empty: all + QMap _outs; // per stream + QMap _outNames; // per stream double _interval; bncTime _intervalStart; bncTime _lastTime; - bncTime _startTime; // first correction epoch + QMap _staStart; // first correction epoch per stream + QMap _staLast; // latest correction epoch per stream bncEphUser _ephUser; QMap _epoStat; // key "STA SYS TYPE" Index: /trunk/BNC/tests/README.md =================================================================== --- /trunk/BNC/tests/README.md (revision 11074) +++ /trunk/BNC/tests/README.md (revision 11075) @@ -79,5 +79,6 @@ * `ssr_rtcm_SSRA00CNE0`, `ssr_rtcm_WHU1_OSBC`: RTCM 1059-1270; these test whichever class `RTCM3Decoder.cpp` selects for non-4076 SSR. -* `ssr_qc`: SSR QC report (`ssrQcFile`, 5 min interval) of the 7 AC streams of +* `ssr_qc`, `ssr_qc_select`: SSR QC (31_SsrQC, one file per stream, 5 min interval; + all streams resp. `ssrQcStreams SSRA00CNE1,SSRA00GMV0`) of the 7 AC streams of `data/combi_multi_261002_5min.raw`. Index: /trunk/BNC/tests/cases/ssr_qc/case.json =================================================================== --- /trunk/BNC/tests/cases/ssr_qc/case.json (revision 11074) +++ /trunk/BNC/tests/cases/ssr_qc/case.json (revision 11075) @@ -1,4 +1,4 @@ { - "description": "31_SsrQC replay: SSR QC report of the 7 AC streams of data/combi_multi_261002_5min.raw, 5 min report interval", + "description": "31_SsrQC replay: SSR QC reports (one file per stream, ALL streams) of the 7 AC streams of data/combi_multi_261002_5min.raw, 5 min report interval", "config": "Example_Configs/31_SsrQC.bnc", "raw": "data/combi_multi_261002_5min.raw", @@ -11,8 +11,5 @@ ], "outputs": [ - { - "glob": "Output/SsrQC.txt", - "name": "ssrqc.txt" - } + "Output/*_QC.txt" ] } Index: /trunk/BNC/tests/cases/ssr_qc_select/case.json =================================================================== --- /trunk/BNC/tests/cases/ssr_qc_select/case.json (revision 11075) +++ /trunk/BNC/tests/cases/ssr_qc_select/case.json (revision 11075) @@ -0,0 +1,18 @@ +{ + "description": "31_SsrQC replay with a stream selection (ssrQcStreams SSRA00CNE1,SSRA00GMV0): files only for these streams", + "config": "Example_Configs/31_SsrQC.bnc", + "raw": "data/combi_multi_261002_5min.raw", + "exit_code": 5, + "timeout": 600, + "args": [ + "-key", + "ssrQcIntr", + "5 min", + "-key", + "ssrQcStreams", + "SSRA00CNE1,SSRA00GMV0" + ], + "outputs": [ + "Output/*_QC.txt" + ] +}