Index: trunk/BNC/src/ssrQc.cpp =================================================================== --- trunk/BNC/src/ssrQc.cpp (revision 11073) +++ 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