Changeset 8905 in ntrip for trunk/BNC/src/PPP/pppFilter.cpp
- Timestamp:
- Mar 18, 2020, 11:13:50 AM (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/BNC/src/PPP/pppFilter.cpp
r8329 r8905 34 34 // Constructor 35 35 //////////////////////////////////////////////////////////////////////////// 36 t_pppFilter::t_pppFilter( ) {36 t_pppFilter::t_pppFilter(t_pppObsPool* obsPool) { 37 37 _parlist = 0; 38 _numSat = 0; 39 _obsPool = obsPool; 38 40 } 39 41 … … 46 48 // Process Single Epoch 47 49 //////////////////////////////////////////////////////////////////////////// 48 t_irc t_pppFilter::processEpoch(t_pppObsPool* obsPool) { 49 50 t_irc t_pppFilter::processEpoch(int num) { 50 51 _numSat = 0; 51 52 const double maxSolGap = 60.0; … … 57 58 // Vector of all Observations 58 59 // -------------------------- 59 t_pppObsPool::t_epoch* epoch = obsPool->lastEpoch();60 t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch(); 60 61 if (!epoch) { 61 62 return failure; … … 75 76 string epoTimeStr = string(_epoTime); 76 77 78 //-- 77 79 // Set Parameters 78 80 // -------------- … … 91 93 for (unsigned ii = 0; ii < params.size(); ii++) { 92 94 const t_pppParam* par1 = params[ii]; 93 94 95 _x0[ii] = par1->x0(); 95 96 96 int iOld = par1->indexOld(); 97 97 if (iOld < 0) { … … 113 113 predictCovCrdPart(QFltOld); 114 114 115 116 // Pre-Process Satellite Systems separately 117 // ---------------------------------------- 118 bool preProcessing = false; 119 if (OPT->_obsModelType == OPT->DCMcodeBias || 120 OPT->_obsModelType == OPT->DCMphaseBias) { 121 preProcessing = true; 122 for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) { 123 char system = OPT->systems()[iSys]; 124 t_prn refPrn = (_obsPool->getRefSatMapElement(system))->prn(); 125 vector<t_pppSatObs*> obsVector; 126 for (unsigned jj = 0; jj < allObs.size(); jj++) { 127 if (allObs[jj]->prn().system() == system) { 128 obsVector.push_back(allObs[jj]); 129 } 130 } 131 if (processSystem(OPT->LCs(system), obsVector, refPrn, 132 epoch->pseudoObsIono(), preProcessing) != success) { 133 return failure; 134 } 135 } 136 } 137 115 138 // Process Satellite Systems separately 116 139 // ------------------------------------ 117 140 for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) { 118 141 char system = OPT->systems()[iSys]; 142 t_prn refPrn = (_obsPool->getRefSatMapElement(system))->prn(); 119 143 unsigned int num = 0; 120 144 vector<t_pppSatObs*> obsVector; … … 126 150 } 127 151 LOG << epoTimeStr << " SATNUM " << system << ' ' << right << setw(2) << num << endl; 128 if ( processSystem(OPT->LCs(system), obsVector) != success ) { 152 if (processSystem(OPT->LCs(system), obsVector, refPrn, 153 epoch->pseudoObsIono(), preProcessing) != success) { 129 154 return failure; 130 155 } 131 156 } 132 133 cmpDOP(allObs); 134 135 _parlist->printResult(_epoTime, _QFlt, _xFlt); 136 _lastEpoTimeOK = _epoTime; // remember time of last successful epoch processing 157 if (_obsPool->epoReProcessing()) { 158 // set A1 und A2 abhängig von num 159 // if num == 1 => A1 160 // if num > 1 => A2 161 } 162 else { 163 cmpDOP(allObs); 164 _parlist->printResult(_epoTime, _QFlt, _xFlt); 165 _lastEpoTimeOK = _epoTime; // remember time of last successful epoch processing 166 } 167 137 168 return success; 138 169 } … … 141 172 //////////////////////////////////////////////////////////////////////////// 142 173 t_irc t_pppFilter::processSystem(const vector<t_lc::type>& LCs, 143 const vector<t_pppSatObs*>& obsVector) { 144 174 const vector<t_pppSatObs*>& obsVector, 175 const t_prn& refPrn, 176 bool pseudoObsIonoAvailable, 177 bool preProcessing) { 178 qDebug() << "======t_pppFilter::processSystem======="; 145 179 LOG.setf(ios::fixed); 146 180 147 181 // Detect Cycle Slips 148 182 // ------------------ 149 if (detectCycleSlips(LCs, obsVector ) != success) {183 if (detectCycleSlips(LCs, obsVector, refPrn, preProcessing) != success) { 150 184 return failure; 151 185 } 152 186 153 ColumnVector xSav = _xFlt; 154 SymmetricMatrix QSav = _QFlt; 155 string epoTimeStr = string(_epoTime); 187 unsigned usedLCs = LCs.size(); //qDebug() << "usedLCs: " << usedLCs; 188 if (OPT->_pseudoObsIono && !pseudoObsIonoAvailable) { 189 usedLCs -= 1; // GIM not used 190 } 191 ColumnVector xSav = _xFlt; 192 SymmetricMatrix QSav = _QFlt; 193 string epoTimeStr = string(_epoTime); 156 194 const vector<t_pppParam*>& params = _parlist->params(); 157 unsigned maxObs = obsVector.size() * LCs.size(); 195 unsigned maxObs = obsVector.size() * usedLCs; 196 //unsigned maxObs = 2 * usedLCs; 197 198 if (OPT->_pseudoObsIono && pseudoObsIonoAvailable) { // stecDiff w.r.t refSat 199 maxObs -= 1; 200 } 201 qDebug() << "par.size() : " << _parlist->nPar() << " LCs.size(): " << usedLCs; 202 qDebug() << "obsVector.size(): " << obsVector.size() << " maxObs : " << maxObs; 158 203 159 204 // Outlier Detection Loop 160 205 // ---------------------- 161 206 for (unsigned iOutlier = 0; iOutlier < maxObs; iOutlier++) { 207 qDebug() << "iOutlier: " << iOutlier; 162 208 163 209 if (iOutlier > 0) { … … 168 214 // First-Design Matrix, Terms Observed-Computed, Weight Matrix 169 215 // ----------------------------------------------------------- 216 qDebug() << "A(" << maxObs << "," << _parlist->nPar() << ")"; 170 217 Matrix AA(maxObs, _parlist->nPar()); 171 218 ColumnVector ll(maxObs); 172 219 DiagonalMatrix PP(maxObs); PP = 0.0; 220 //TETSPLOT 221 for (unsigned iPar = 0; iPar < params.size(); iPar++) { 222 const t_pppParam* par = params[iPar]; 223 cout << " " << par->toString() <<endl; 224 } 225 cout << endl; 173 226 174 227 int iObs = -1; … … 176 229 vector<t_lc::type> usedTypes; 177 230 for (unsigned ii = 0; ii < obsVector.size(); ii++) { 178 t_pppSatObs* obs = obsVector[ii]; 231 t_pppSatObs* obs = obsVector[ii]; qDebug() << "SATELLITE: " << obs->prn().toString().c_str() << "isRef: " << obs->isReference(); 179 232 if (!obs->outlier()) { 180 for (unsigned jj = 0; jj < LCs.size(); jj++) {233 for (unsigned jj = 0; jj < usedLCs; jj++) { 181 234 const t_lc::type tLC = LCs[jj]; 235 if (tLC == t_lc::GIM && obs->isReference()) {continue;} 182 236 ++iObs; 183 237 usedObs.push_back(obs); … … 186 240 const t_pppParam* par = params[iPar]; 187 241 AA[iObs][iPar] = par->partial(_epoTime, obs, tLC); 242 cout << setw(10) << par->partial(_epoTime, obs, tLC); 188 243 } 244 cout << endl; 189 245 ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1)); 190 246 PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC)); 247 //qDebug() << "ll[iObs]: " << ll[iObs]; 248 //qDebug() << "PP[iObs]: " << PP[iObs]; 191 249 } 192 250 } … … 229 287 t_pppSatObs* obs = usedObs[maxOutlierIndex]; 230 288 t_pppParam* par = 0; 231 LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' ' 232 << obs->prn().toString() << ' ' 233 << setw(8) << setprecision(4) << maxOutlier << endl; 289 if (!preProcessing) { 290 LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' ' 291 << obs->prn().toString() << ' ' 292 << setw(8) << setprecision(4) << maxOutlier << endl; 293 } 234 294 for (unsigned iPar = 0; iPar < params.size(); iPar++) { 235 295 t_pppParam* hlp = params[iPar]; 236 if (hlp->type() == t_pppParam::amb && hlp->prn() == obs->prn() && 296 if (hlp->type() == t_pppParam::amb && 297 hlp->prn() == obs->prn() && 237 298 hlp->tLC() == usedTypes[maxOutlierIndex]) { 238 299 par = hlp; 239 300 } 240 301 } 241 if (par) { 302 if (par && preProcessing) { 303 if (par->prn() == refPrn) { 304 _obsPool->setEpoReProcessing(true); 305 } 306 } 307 else if (par && !preProcessing) { 242 308 if (par->ambResetCandidate()) { 243 309 resetAmb(par->prn(), obsVector, &QSav, &xSav); … … 248 314 } 249 315 } 250 else {316 else if (!par && !preProcessing){ 251 317 obs->setOutlier(); 252 318 } … … 256 322 // --------------- 257 323 else { 258 for (unsigned jj = 0; jj < LCs.size(); jj++) { 259 for (unsigned ii = 0; ii < usedObs.size(); ii++) { 260 const t_lc::type tLC = usedTypes[ii]; 261 t_pppSatObs* obs = usedObs[ii]; 262 if (tLC == LCs[jj]) { 263 obs->setRes(tLC, vv[ii]); 264 LOG << epoTimeStr << " RES " 265 << left << setw(3) << t_lc::toString(tLC) << right << ' ' 266 << obs->prn().toString().substr(0,3) << ' ' 267 << setw(8) << setprecision(4) << vv[ii] << endl; 324 if (!preProcessing) { 325 for (unsigned jj = 0; jj < LCs.size(); jj++) { 326 for (unsigned ii = 0; ii < usedObs.size(); ii++) { 327 const t_lc::type tLC = usedTypes[ii]; 328 t_pppSatObs* obs = usedObs[ii]; 329 if (tLC == LCs[jj]) { 330 obs->setRes(tLC, vv[ii]); 331 LOG << epoTimeStr << " RES " 332 << left << setw(3) << t_lc::toString(tLC) << right << ' ' 333 << obs->prn().toString().substr(0,3) << ' ' 334 << setw(8) << setprecision(4) << vv[ii] << endl; 335 } 268 336 } 269 337 } … … 279 347 //////////////////////////////////////////////////////////////////////////// 280 348 t_irc t_pppFilter::detectCycleSlips(const vector<t_lc::type>& LCs, 281 const vector<t_pppSatObs*>& obsVector) { 349 const vector<t_pppSatObs*>& obsVector, 350 const t_prn& refPrn, 351 bool preProcessing) { 282 352 283 353 const double SLIP = 20.0; // slip threshold 284 354 string epoTimeStr = string(_epoTime); 285 const vector<t_pppParam*>& params 355 const vector<t_pppParam*>& params = _parlist->params(); 286 356 287 357 for (unsigned ii = 0; ii < LCs.size(); ii++) { … … 294 364 // --------------------------------- 295 365 bool slip = false; 296 297 366 if (obs->slip()) { 298 367 LOG << epoTimeStr << "cycle slip set (obs) " << obs->prn().toString() << endl; … … 312 381 slip = true; 313 382 } 314 _slips[obs->prn()]._biasJumpCounter = obs->biasJumpCounter();315 383 316 384 // Slip Set 317 385 // -------- 318 386 if (slip) { 319 resetAmb(obs->prn(), obsVector); 320 } 321 387 if (preProcessing) { 388 if (obs->prn() == refPrn) { 389 _obsPool->setEpoReProcessing(true); 390 } 391 } 392 else { 393 resetAmb(obs->prn(), obsVector); 394 } 395 } 322 396 // Check Pre-Fit Residuals 323 397 // ----------------------- … … 328 402 AA[iPar] = par->partial(_epoTime, obs, tLC); 329 403 } 330 331 404 double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA); 332 405 double vv = DotProduct(AA, _xFlt) - ll; 333 334 406 if (fabs(vv) > SLIP) { 335 LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' ' 336 << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl; 337 resetAmb(obs->prn(), obsVector); 407 if (preProcessing) { 408 if (obs->prn() == refPrn) { 409 _obsPool->setEpoReProcessing(true); 410 } 411 } 412 else { 413 LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' ' 414 << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl; 415 resetAmb(obs->prn(), obsVector); 416 } 338 417 } 339 418 }
Note:
See TracChangeset
for help on using the changeset viewer.