Index: /trunk/BNC/src/PPP/pppClient.cpp =================================================================== --- /trunk/BNC/src/PPP/pppClient.cpp (revision 11066) +++ /trunk/BNC/src/PPP/pppClient.cpp (revision 11067) @@ -61,4 +61,5 @@ if (!_opt->_antexFileName.empty()) { _antex = new bncAntex(_opt->_antexFileName.c_str()); + *_log << _antex->info().toStdString() << endl; } if (!_opt->_blqFileName.empty()) { Index: /trunk/BNC/src/PPP/pppSatObs.cpp =================================================================== --- /trunk/BNC/src/PPP/pppSatObs.cpp (revision 11066) +++ /trunk/BNC/src/PPP/pppSatObs.cpp (revision 11067) @@ -672,19 +672,24 @@ } else { - _model._antPCO[ii] += PPP_CLIENT->antex()->satCorr(prn, frqType, _model._elTx, _model._azTx, found); + _model._antPCO[ii] += PPP_CLIENT->antex()->satCorr(prn, frqType, _model._elTx, _model._azTx, found, _time); + if (!found && !PPP_CLIENT->antex()->hasSatEntry(prn, _time) && + PPP_CLIENT->antex()->warnNoSatEntry(prn, _time)) { + LOG << "ANTEX: no satellite antenna entry valid at " << _time.datestr() << ' ' << _time.timestr() + << " for " << _prn.toString() << " - ANTEX file outdated?" << endl; + } if (OPT->_isAPC && found) { // the PCOs as given in the satellite antenna correction for all frequencies // have to be reduced by the PCO of the respective reference frequency if (_prn.system() == 'G') { - _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::G1, _model._elTx, _model._azTx, found); + _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::G1, _model._elTx, _model._azTx, found, _time); } else if (_prn.system() == 'R') { - _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::R1, _model._elTx, _model._azTx, found); + _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::R1, _model._elTx, _model._azTx, found, _time); } else if (_prn.system() == 'E') { - _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::E1, _model._elTx, _model._azTx, found); + _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::E1, _model._elTx, _model._azTx, found, _time); } else if (_prn.system() == 'C') { - _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::C2, _model._elTx, _model._azTx, found); + _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::C2, _model._elTx, _model._azTx, found, _time); } } Index: /trunk/BNC/src/combination/bnccomb.cpp =================================================================== --- /trunk/BNC/src/combination/bnccomb.cpp (revision 11066) +++ /trunk/BNC/src/combination/bnccomb.cpp (revision 11067) @@ -19,4 +19,5 @@ #include #include +#include #include "bnccomb.h" @@ -290,4 +291,7 @@ _antex = 0; } + else { + emit newMessage("bncComb: " + _antex->info().toLatin1(), true); + } } @@ -853,9 +857,9 @@ // Check satellite code biases // ---------------------------- + QMap codeBiasesRefSig; if (_satCodeBiases.contains(acName)) { QMap& storage = _satCodeBiases[acName]; if (storage.contains(clkCorr._prn)) { _newCorr->_satCodeBias = storage[clkCorr._prn]; - QMap codeBiasesRefSig; for (unsigned ii = 1; ii < cmbRefSig::cIF; ii++) { t_frequency::type frqType = cmbRefSig::toFreq(sys, static_cast(ii)); @@ -870,34 +874,16 @@ } } - // Code biases present, but not for both reference signals (e.g. - // Galileo HAS): the AC clocks cannot be referred to the reference - // signals and are used as they are - warn (at most once per hour - // per AC and system) - if (codeBiasesRefSig.size() < 2) { - QString key = acName + sys; - bncTime& warned = _refBiasWarned[key]; - if (!warned.valid() || _newCorr->_time - warned >= 3600.0) { - warned = _newCorr->_time; - QString refSig; - for (unsigned ii = 1; ii < cmbRefSig::cIF; ii++) { - t_frequency::type frqType = cmbRefSig::toFreq(sys, static_cast(ii)); - refSig += QString(" C%1%2").arg(t_frequency::toString(frqType)[1]) - .arg(cmbRefSig::toAttrib(sys, static_cast(ii))); - } - emit newMessage("bncComb: " + acName.toLatin1() + " provides no code biases for the " + - QByteArray(1, sys) + " reference signals" + refSig.toLatin1() + - " - its clocks are combined without referring them to these signals", false); + // only with the biases of both reference signals (one of them alone + // does not refer the clocks to the reference signals) + if (codeBiasesRefSig.size() == 2) { + map codeCoeff; + double channel = double(_newCorr->_eph->slotNum()); + cmbRefSig::coeff(sys, cmbRefSig::cIF, channel, codeCoeff); + map::const_iterator it; + for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) { + t_frequency::type frqType = it->first; + double codeCoeff = it->second; + _newCorr->_satCodeBiasIF += codeCoeff * codeBiasesRefSig[frqType]; } - } - map codeCoeff; - double channel = double(_newCorr->_eph->slotNum()); - cmbRefSig::coeff(sys, cmbRefSig::cIF, channel, codeCoeff); - map::const_iterator it; - for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) { - t_frequency::type frqType = it->first; - double codeCoeff = it->second; - _newCorr->_satCodeBiasIF += codeCoeff * codeBiasesRefSig[frqType]; - } - if (codeBiasesRefSig.size() == 2) { _acBiasIF[acName][prnStr.mid(0,3)] = _newCorr->_satCodeBiasIF; } @@ -906,4 +892,31 @@ } _newCorr->_satCodeBias._bias.clear(); + } + } + // Without the code biases of both reference signals (none at all, or + // only one of them, e.g. Galileo HAS for GPS), the AC clock cannot be + // referred to the reference signals: it is combined as it is, but + // excluded from the clock datum (see createAmat) - warn at most once + // per hour per AC and system + _newCorr->_refBiasesOK = (codeBiasesRefSig.size() == 2); + if (!_newCorr->_refBiasesOK) { + _acBiasIF[acName].remove(prnStr.mid(0,3)); + QString key = acName + sys; + bncTime& warned = _refBiasWarned[key]; + if (!warned.valid() || _newCorr->_time - warned >= 3600.0) { + warned = _newCorr->_time; + QString missing; + for (unsigned ii = 1; ii < cmbRefSig::cIF; ii++) { + t_frequency::type frqType = cmbRefSig::toFreq(sys, static_cast(ii)); + if (!codeBiasesRefSig.contains(frqType)) { + missing += QString(" C%1%2").arg(t_frequency::toString(frqType)[1]) + .arg(cmbRefSig::toAttrib(sys, static_cast(ii))); + } + } + emit newMessage("bncComb: " + acName.toLatin1() + " provides no code biases for the " + + QByteArray(1, sys) + " reference signal(s)" + missing.toLatin1() + + " (e.g. " + prnStr.mid(0,3).toLatin1() + ") - these clocks are combined" + " without referring them to the reference signals and, as long as other ACs" + " provide referred clocks, excluded from the clock datum", false); } } @@ -1087,4 +1100,5 @@ if (_method == filter) { checkBiasConsistency(epoTime, sys, out); + checkPersistentOutliers(epoTime, sys, out); } printResults(epoTime, out, masterCorr); @@ -1121,4 +1135,11 @@ if (checkOrbits(epoTime, sys, out, masterCorr) != success) { return failure; + } + + // Statistics of persistent outliers + // --------------------------------- + for (int ii = 0; ii < corrs(sys).size(); ii++) { + const cmbCorr* corr = corrs(sys)[ii]; + _outlierStat[sys][corr->_acName + ' ' + corr->_prn.mid(0,3)].numEpo += 1; } @@ -1161,4 +1182,6 @@ out << " Outlier" << "\n"; + _outlierStat[sys][corrs(sys)[maxResIndex-1]->_acName + ' ' + + corrs(sys)[maxResIndex-1]->_prn.mid(0,3)].residuals.push_back(maxRes); _QQ[sys] = QQ_sav; delete corrs(sys)[maxResIndex-1]; @@ -1254,11 +1277,17 @@ sxy += (xx[ii] - mx) * (yy[ii] - my); } + QString key = acName + sys; double spread = sqrt(sxx / nn); - if (spread < MIN_SPREAD || syy <= 0.0) { - continue; - } - double slope = sxy / sxx; - double corrCoeff = sxy / sqrt(sxx * syy); - if (slope < MIN_SLOPE || corrCoeff < MIN_CORR) { + double slope = (sxx > 0.0) ? sxy / sxx : 0.0; + double corrCoeff = (sxx > 0.0 && syy > 0.0) ? sxy / sqrt(sxx * syy) : 0.0; + if (spread < MIN_SPREAD || slope < MIN_SLOPE || corrCoeff < MIN_CORR) { + if (_biasInconsistent.remove(key)) { + QString msg = QString("%1 %2 code biases consistent with its clocks again, " + "AC used for the clock datum again").arg(acName).arg(sys); + out << epoTime.datestr().c_str() << " " << epoTime.timestr().c_str() + << " " << msg << "\n"; + emit newMessage("bncComb: " + msg.toLatin1(), true); + resetSatOffsets(sys); + } continue; } @@ -1275,11 +1304,73 @@ .arg(acName).arg(sys).arg(refSig).arg(nn) .arg(spread, 0, 'f', 3).arg(slope, 0, 'f', 2).arg(corrCoeff, 0, 'f', 2); + if (!_biasInconsistent.contains(key)) { + _biasInconsistent.insert(key); + msg += ", AC excluded from the clock datum"; + _biasCheckWarned.remove(key); + resetSatOffsets(sys); + } out << epoTime.datestr().c_str() << " " << epoTime.timestr().c_str() << " " << msg << "\n"; - QString key = acName + sys; bncTime& warned = _biasCheckWarned[key]; if (!warned.valid() || epoTime - warned >= 3600.0) { warned = epoTime; emit newMessage("bncComb: " + msg.toLatin1(), true); + } + } +} + +// Report AC satellites rejected as outliers in most epochs +//////////////////////////////////////////////////////////////////////////// +// A persistent clock difference of one AC for a satellite cannot be absorbed +// by its satellite offset, as the offsets of all ACs of the satellite sum up +// to zero: the correction is rejected again in every epoch. Once per hour +// and system, AC satellites rejected in at least half of their epochs are +// reported. +//////////////////////////////////////////////////////////////////////////// +void bncComb::checkPersistentOutliers(bncTime epoTime, char sys, QTextStream& out) { + + const double INTERVAL = 3600.0; + const unsigned MIN_EPO = 10; + + if (!_outlierStatStart[sys].valid()) { + _outlierStatStart[sys] = epoTime; + return; + } + if (epoTime - _outlierStatStart[sys] < INTERVAL) { + return; + } + QMapIterator it(_outlierStat[sys]); + while (it.hasNext()) { + it.next(); + const t_outlierStat& stat = it.value(); + unsigned numOut = stat.residuals.size(); + if (stat.numEpo < MIN_EPO || 2 * numOut < stat.numEpo) { + continue; + } + QVector res = stat.residuals; + std::sort(res.begin(), res.end()); + double median = (numOut % 2) ? res[numOut/2] : 0.5 * (res[numOut/2-1] + res[numOut/2]); + QString msg = QString("%1 rejected as outlier in %2% of %3 epochs in the last hour " + "(median residual %4 m)") + .arg(it.key()).arg(100.0 * numOut / stat.numEpo, 0, 'f', 0) + .arg(stat.numEpo).arg(median, 0, 'f', 2); + out << epoTime.datestr().c_str() << " " << epoTime.timestr().c_str() << " " << msg << "\n"; + emit newMessage("bncComb: " + msg.toLatin1(), true); + } + _outlierStat[sys].clear(); + _outlierStatStart[sys] = epoTime; +} + +// Re-initialize the satellite offsets of all ACs of a system (after a change +// of the clock datum) +//////////////////////////////////////////////////////////////////////////// +void bncComb::resetSatOffsets(char sys) { + for (int iPar = 1; iPar <= _params[sys].size(); iPar++) { + cmbParam* pp = _params[sys][iPar-1]; + if (pp->type == cmbParam::offACSat) { + pp->xx = 0.0; + _QQ[sys].Row(iPar) = 0.0; + _QQ[sys].Column(iPar) = 0.0; + _QQ[sys](iPar,iPar) = pp->sig0 * pp->sig0; } } @@ -1411,6 +1502,17 @@ if (_antex->satCoMcorrection(corr->_prn, Mjd, xc.Rows(1,3), vv, dx, attMode, extYaw) != success) { - dx = 0; - emit newMessage("bncComb: antenna not found " + corr->_prn.mid(0,3).toLatin1(), false); + // without the satellite antenna offset, the orbit cannot be + // converted between antenna phase center and center of mass + if (_antex->warnNoSatEntry(corr->_prn, epoTime)) { + emit newMessage("bncComb: no satellite antenna entry valid at " + + QByteArray(epoTime.datestr().c_str()) + " " + QByteArray(epoTime.timestr().c_str()) + + " for " + corr->_prn.mid(0,3).toLatin1() + + " in ANTEX file - ANTEX file outdated? Satellite not combined", true); + } + orbCorrections.pop_back(); + clkCorrections.pop_back(); + delete corr; + it.remove(); + continue; } QString glonassYawLog = _antex->takeGlonassYawLog(); @@ -1537,5 +1639,34 @@ int maxSat = _cmbSysPrn[sys]; - const int nCon = (_method == filter) ? 1 + maxSat : 0; + // AC clocks not referred to the reference signals are excluded from the + // zero-sum conditions of the satellite offsets (the clock datum): those of + // ACs whose code biases are inconsistent with their clocks (see + // checkBiasConsistency) and those without the code biases of both + // reference signals. This applies as long as another AC provides referred + // corrections; the sum of the excluded offsets of an AC is constrained to + // zero instead. + QSet excludedOff; // "AC PRN" (full PRN) + QStringList excludedACs; // ACs with excluded offsets + if (_method == filter) { + bool referredAC = false; + for (int ii = 0; ii < corrs(sys).size(); ii++) { + const cmbCorr* corr = corrs(sys)[ii]; + if (_biasInconsistent.contains(corr->_acName + sys) || !corr->_refBiasesOK) { + excludedOff.insert(corr->_acName + ' ' + corr->_prn); + if (!excludedACs.contains(corr->_acName)) { + excludedACs << corr->_acName; + } + } + else { + referredAC = true; + } + } + if (!referredAC) { + excludedOff.clear(); + excludedACs.clear(); + } + } + + const int nCon = (_method == filter) ? 1 + maxSat + excludedACs.size() : 0; AA.ReSize(nObs+nCon, nPar); AA = 0.0; @@ -1589,4 +1720,6 @@ // } int iCond = 1; + QSet offExcluded; // offsets excluded from the zero-sum condition of their satellite + QSet acsInDatum; // ACs with at least one offset in the zero-sum conditions // GNSS for (unsigned iGnss = 1; iGnss <= _cmbSysPrn[sys]; iGnss++) { @@ -1595,4 +1728,5 @@ ++iCond; PP(nObs+iCond) = Ph; + QList offAll, offConsistent; for (int iPar = 1; iPar <= _params[sys].size(); iPar++) { cmbParam* pp = _params[sys][iPar-1]; @@ -1601,4 +1735,33 @@ pp->type == cmbParam::offACSat && pp->prn == prn) { + offAll << iPar; + if (!excludedOff.contains(pp->AC + ' ' + pp->prn)) { + offConsistent << iPar; + } + } + } + // satellite observed by inconsistent ACs only: all its offsets are used + const QList& offUsed = offConsistent.isEmpty() ? offAll : offConsistent; + for (int iPar : offAll) { + if (offUsed.contains(iPar)) { + AA(nObs+iCond, iPar) = 1.0; + acsInDatum.insert(_params[sys][iPar-1]->AC); + } + else { + offExcluded << iPar; + } + } + } + // ACs with all their offsets excluded from the zero-sum conditions above: + // the sum of their offsets is constrained to zero (otherwise not + // separable from the AC's epoch offset) + for (const QString& acName : excludedACs) { + ++iCond; + if (acsInDatum.contains(acName)) { + continue; + } + PP(nObs+iCond) = Ph; + for (int iPar : offExcluded) { + if (_params[sys][iPar-1]->AC == acName) { AA(nObs+iCond, iPar) = 1.0; } Index: /trunk/BNC/src/combination/bnccomb.h =================================================================== --- /trunk/BNC/src/combination/bnccomb.h (revision 11066) +++ /trunk/BNC/src/combination/bnccomb.h (revision 11067) @@ -117,4 +117,5 @@ _dClkResult = 0.0; _satCodeBiasIF = 0.0; + _refBiasesOK = false; _lambdaIF = 0.0; _satYawAngle = 0.0; @@ -136,4 +137,5 @@ QString _acName; double _satCodeBiasIF; + bool _refBiasesOK; // code biases of both reference signals available double _lambdaIF; double _satYawAngle; @@ -274,5 +276,15 @@ QMap _biasCheckLast; // last bias/clock consistency check per system QMap _biasCheckWarned; // last inconsistency warning per AC and system + QSet _biasInconsistent; // AC+system excluded from the clock datum void checkBiasConsistency(bncTime epoTime, char sys, QTextStream& out); + void checkPersistentOutliers(bncTime epoTime, char sys, QTextStream& out); + struct t_outlierStat { + t_outlierStat() : numEpo(0) {} + unsigned numEpo; // epochs with a correction of the AC for the satellite + QVector residuals; // residuals of the epochs it was rejected + }; + QMap > _outlierStat; // per system, key "AC PRN" + QMap _outlierStatStart; // start of the statistics interval + void resetSatOffsets(char sys); bool allACsBeyond(const bncTime& epoTime) const; bncTime _resTime;