source: ntrip/trunk/BNC/src/PPP/pppParlist.cpp@ 9289

Last change on this file since 9289 was 9289, checked in by stuerze, 3 years ago

updates regarding inter system biases
  • Property svn:keywords set to Author Date Id Rev URL;svn:eol-style=native
  • Property svn:mime-type set to text/plain
File size: 16.4 KB
Line 
1/* -------------------------------------------------------------------------
2 * BKG NTRIP Client
3 * -------------------------------------------------------------------------
4 *
5 * Class: t_pppParlist
6 *
7 * Purpose: List of estimated parameters
8 *
9 * Author: L. Mervart
10 *
11 * Created: 29-Jul-2014
12 *
13 * Changes:
14 *
15 * -----------------------------------------------------------------------*/
16
17#include <cmath>
18#include <iostream>
19#include <sstream>
20#include <iomanip>
21#include <algorithm>
22#include <newmatio.h>
23
24#include "pppParlist.h"
25#include "pppSatObs.h"
26#include "pppStation.h"
27#include "bncutils.h"
28#include "bncconst.h"
29#include "pppClient.h"
30
31using namespace BNC_PPP;
32using namespace std;
33
34// Constructor
35////////////////////////////////////////////////////////////////////////////
36t_pppParam::t_pppParam(e_type type, const t_prn& prn, t_lc::type tLC,
37 const vector<t_pppSatObs*>* obsVector) {
38
39 _type = type;
40 _prn = prn;
41 _tLC = tLC;
42 _x0 = 0.0;
43 _indexOld = -1;
44 _indexNew = -1;
45 _noise = 0.0;
46 _ambInfo = 0;
47
48 switch (_type) {
49 case crdX:
50 _epoSpec = false;
51 _sigma0 = OPT->_aprSigCrd[0];
52 _noise = OPT->_noiseCrd[0];
53 break;
54 case crdY:
55 _epoSpec = false;
56 _sigma0 = OPT->_aprSigCrd[1];
57 _noise = OPT->_noiseCrd[1];
58 break;
59 case crdZ:
60 _epoSpec = false;
61 _sigma0 = OPT->_aprSigCrd[2];
62 _noise = OPT->_noiseCrd[2];
63 break;
64 case clkR:
65 _epoSpec = true;
66 _sigma0 = OPT->_noiseClk;
67 break;
68 case amb:
69 _epoSpec = false;
70 _sigma0 = OPT->_aprSigAmb;
71 _ambInfo = new t_ambInfo();
72 if (obsVector) {
73 for (unsigned ii = 0; ii < obsVector->size(); ii++) {
74 const t_pppSatObs* obs = obsVector->at(ii);
75 if (obs->prn() == _prn) {
76 double offGR = 0;
77 if (_prn.system() == 'R' && tLC != t_lc::MW) {
78 offGR = PPP_CLIENT->offGR();
79 }
80 double offGE = 0;
81 if (_prn.system() == 'E' && tLC != t_lc::MW) {
82 offGE = PPP_CLIENT->offGE();
83 }
84 double offGC = 0;
85 if (_prn.system() == 'C' && tLC != t_lc::MW) {
86 offGC = PPP_CLIENT->offGC();
87 }
88 _x0 = floor((obs->obsValue(tLC) - offGR -offGE - offGC - obs->cmpValue(tLC)) / obs->lambda(tLC) + 0.5);
89 break;
90 }
91 }
92 }
93 break;
94 case offGR:
95 _epoSpec = true;
96 _sigma0 = 1000.0;
97 _x0 = PPP_CLIENT->offGR();
98 break;
99 case offGE:
100 _epoSpec = true;
101 _sigma0 = 1000.0;
102 _x0 = PPP_CLIENT->offGE();
103 break;
104 case offGC:
105 _epoSpec = true;
106 _sigma0 = 1000.0;
107 _x0 = PPP_CLIENT->offGC();
108 break;
109 case trp:
110 _epoSpec = false;
111 _sigma0 = OPT->_aprSigTrp;
112 _noise = OPT->_noiseTrp;
113 break;
114 case ion:
115 _epoSpec = false;
116 _sigma0 = OPT->_aprSigIon;
117 _noise = OPT->_noiseIon;
118 break;
119 case cBias1:
120 case cBias2:
121 _epoSpec = false;
122 _sigma0 = OPT->_aprSigCodeBias;
123 _noise = OPT->_noiseCodeBias;
124 break;
125 case pBias1:
126 case pBias2:
127 _epoSpec = false;
128 _sigma0 = OPT->_aprSigPhaseBias;
129 _noise = OPT->_noisePhaseBias;
130 break;
131 }
132}
133
134// Destructor
135////////////////////////////////////////////////////////////////////////////
136t_pppParam::~t_pppParam() {
137 delete _ambInfo;
138}
139
140//
141////////////////////////////////////////////////////////////////////////////
142double t_pppParam::partial(const bncTime& /* epoTime */, const t_pppSatObs* obs,
143 const t_lc::type& tLC, const t_prn refPrn) const {
144
145 // Special Case - Melbourne-Wuebbena
146 // ---------------------------------
147 if (tLC == t_lc::MW && _type != amb) {
148 return 0.0;
149 }
150
151 const t_pppStation* sta = PPP_CLIENT->staRover();
152 ColumnVector rhoV = sta->xyzApr() - obs->xc().Rows(1,3);
153
154 map<t_frequency::type, double> codeCoeff;
155 map<t_frequency::type, double> phaseCoeff;
156 map<t_frequency::type, double> ionoCoeff;
157 obs->lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
158
159 switch (_type) {
160 case crdX:
161 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
162 return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.NormFrobenius();
163 case crdY:
164 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
165 return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.NormFrobenius();
166 case crdZ:
167 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
168 return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.NormFrobenius();
169 case clkR:
170 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
171 return 1.0;
172 case offGR:
173 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
174 return (obs->prn().system() == 'R') ? 1.0 : 0.0;
175 case offGE:
176 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
177 return (obs->prn().system() == 'E') ? 1.0 : 0.0;
178 case offGC:
179 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
180 return (obs->prn().system() == 'C') ? 1.0 : 0.0;
181 case amb:
182 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
183 else if ((OPT->_obsModelType == OPT->IF) ||
184 (OPT->_obsModelType == OPT->PPPRTK) ||
185 (OPT->_obsModelType == OPT->UncombPPP) ||
186 (OPT->_obsModelType == OPT->DCMcodeBias && !obs->isReference()) ||
187 (OPT->_obsModelType == OPT->DCMphaseBias && !obs->isReference()) ) {
188 if (obs->prn() == _prn) {
189 if (tLC == _tLC) {
190 return (obs->lambda(tLC));
191 }
192 else if (tLC == t_lc::lIF && _tLC == t_lc::MW) {
193 return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2);
194 }
195 else {
196 if (_tLC == t_lc::l1) {
197 return obs->lambda(t_lc::l1) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
198 }
199 else if (_tLC == t_lc::l2) {
200 return obs->lambda(t_lc::l2) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
201 }
202 }
203 }
204 }
205 break;
206 case trp:
207 if (tLC == t_lc::GIM) {
208 return 0.0;
209 }
210 else if (tLC == t_lc::Tz0) {
211 return 1.0;
212 }
213 else {
214 return 1.0 / sin(obs->eleSat());
215 }
216 case ion:
217 if (obs->prn() == _prn) {
218 if (tLC == t_lc::c1) {
219 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c1)];
220 }
221 else if (tLC == t_lc::c2) {
222 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c2)];
223 }
224 else if (tLC == t_lc::l1) {
225 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
226 }
227 else if (tLC == t_lc::l2) {
228 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
229 }
230 else if (tLC == t_lc::GIM) {
231 return -1.0;
232 }
233 }
234 if (tLC == t_lc::GIM && _prn == refPrn) {
235 return 1.0;
236 }
237 break;
238 case cBias1:
239 if (tLC == t_lc::c1) {
240 return 1.0;
241 }
242 else {
243 return 0.0;
244 }
245 break;
246 case cBias2:
247 if (tLC == t_lc::c2) {
248 return 1.0;
249 }
250 else {
251 return 0.0;
252 }
253 break;
254 case pBias1:
255 if (tLC == t_lc::l1) {
256 return 1.0;
257 }
258 else {
259 return 0.0;
260 }
261 break;
262 case pBias2:
263 if (tLC == t_lc::l2) {
264 return 1.0;
265 }
266 else {
267 return 0.0;
268 }
269 }
270 return 0.0;
271}
272
273//
274////////////////////////////////////////////////////////////////////////////
275string t_pppParam::toString() const {
276 stringstream ss;
277 switch (_type) {
278 case crdX:
279 ss << "CRD_X";
280 break;
281 case crdY:
282 ss << "CRD_Y";
283 break;
284 case crdZ:
285 ss << "CRD_Z";
286 break;
287 case clkR:
288 ss << "REC_CLK ";
289 break;
290 case offGR:
291 ss << "OGR ";
292 break;
293 case offGE:
294 ss << "OGE ";
295 break;
296 case offGC:
297 ss << "OGC ";
298 break;
299 case trp:
300 ss << "TRP ";
301 break;
302 case amb:
303 ss << "AMB " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
304 break;
305 case ion:
306 ss << "ION " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
307 break;
308 case cBias1:
309 case cBias2:
310 case pBias1:
311 case pBias2:
312 ss << "BIAS " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << "REC";
313 break;
314 }
315 return ss.str();
316}
317
318// Constructor
319////////////////////////////////////////////////////////////////////////////
320t_pppParlist::t_pppParlist() {
321}
322
323// Destructor
324////////////////////////////////////////////////////////////////////////////
325t_pppParlist::~t_pppParlist() {
326 for (unsigned ii = 0; ii < _params.size(); ii++) {
327 delete _params[ii];
328 }
329}
330
331//
332////////////////////////////////////////////////////////////////////////////
333t_irc t_pppParlist::set(const bncTime& epoTime,
334 const std::vector<t_pppSatObs*>& obsVector,
335 const QMap<char, t_pppRefSat*>& refSatMap) {
336
337 // Remove some Parameters
338 // ----------------------
339 vector<t_pppParam*>::iterator it = _params.begin();
340 while (it != _params.end()) {
341 t_pppParam* par = *it;
342
343 bool remove = false;
344
345 if (par->epoSpec()) {
346 remove = true;
347 }
348
349 else if (par->type() == t_pppParam::crdX ||
350 par->type() == t_pppParam::crdY ||
351 par->type() == t_pppParam::crdZ) {
352 if (par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 60.0)) {
353 remove = true;
354 }
355 }
356 else if (par->type() == t_pppParam::amb) {
357 char system = par->prn().system();
358 t_prn refPrn = t_prn();
359 if (OPT->_refSatRequired) {
360 refPrn = (refSatMap[system])->prn();
361 }
362 if ((par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 60.0)) ||
363 ( refPrn == par->prn())) {
364 remove = true;
365 }
366 }
367
368 if (remove) {
369 delete par;
370 it = _params.erase(it);
371 }
372 else {
373 ++it;
374 }
375 }
376
377 // Check whether parameters have observations
378 // ------------------------------------------
379 for (unsigned ii = 0; ii < _params.size(); ii++) {
380 t_pppParam* par = _params[ii];
381 if (par->prn() == 0) {
382 par->setLastObsTime(epoTime);
383 if (par->firstObsTime().undef()) {
384 par->setFirstObsTime(epoTime);
385 }
386 }
387 else {
388 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
389 const t_pppSatObs* satObs = obsVector[jj];
390 if (satObs->prn() == par->prn()) {
391 par->setLastObsTime(epoTime);
392 if (par->firstObsTime().undef()) {
393 par->setFirstObsTime(epoTime);
394 }
395 break;
396 }
397 }
398 }
399 }
400
401 // Required Set of Parameters
402 // --------------------------
403 vector<t_pppParam*> required;
404
405 // Coordinates
406 // -----------
407 required.push_back(new t_pppParam(t_pppParam::crdX, t_prn(), t_lc::dummy));
408 required.push_back(new t_pppParam(t_pppParam::crdY, t_prn(), t_lc::dummy));
409 required.push_back(new t_pppParam(t_pppParam::crdZ, t_prn(), t_lc::dummy));
410
411 // Receiver Clock
412 // --------------
413 required.push_back(new t_pppParam(t_pppParam::clkR, t_prn(), t_lc::dummy));
414
415 // GPS-GLONASS Clock Offset
416 // ------------------------
417 if (OPT->useSystem('R')) {
418 required.push_back(new t_pppParam(t_pppParam::offGR, t_prn(), t_lc::dummy));
419 }
420
421 // GPS-Galileo Clock Offset
422 // ------------------------
423 if (OPT->useSystem('E')) {
424 required.push_back(new t_pppParam(t_pppParam::offGE, t_prn(), t_lc::dummy));
425 }
426
427 // GPS-BDS Clock Offset
428 // ------------------------
429 if (OPT->useSystem('C')) {
430 required.push_back(new t_pppParam(t_pppParam::offGC, t_prn(), t_lc::dummy));
431 }
432
433 // Troposphere
434 // -----------
435 if (OPT->estTrp()) {
436 required.push_back(new t_pppParam(t_pppParam::trp, t_prn(), t_lc::dummy));
437 }
438
439 // Ionosphere
440 // ----------
441 if (OPT->_obsModelType == OPT->UncombPPP ||
442 OPT->_obsModelType == OPT->DCMcodeBias ||
443 OPT->_obsModelType == OPT->DCMphaseBias ) {
444 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
445 const t_pppSatObs* satObs = obsVector[jj];
446 required.push_back(new t_pppParam(t_pppParam::ion, satObs->prn(), t_lc::dummy));
447 }
448 }
449
450 // Ambiguities
451 // -----------
452 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
453 const t_pppSatObs* satObs = obsVector[jj];
454 if ((OPT->_obsModelType == OPT->IF) ||
455 (OPT->_obsModelType == OPT->PPPRTK) ||
456 (OPT->_obsModelType == OPT->UncombPPP) ||
457 (OPT->_obsModelType == OPT->DCMcodeBias && !satObs->isReference()) ||
458 (OPT->_obsModelType == OPT->DCMphaseBias && !satObs->isReference()) ) {
459 const vector<t_lc::type>& ambLCs = OPT->ambLCs(satObs->prn().system());
460 for (unsigned ii = 0; ii < ambLCs.size(); ii++) {
461 required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector));
462 }
463 }
464 }
465
466 // Receiver Code Biases
467 // --------------------
468 if (OPT->_obsModelType == OPT->DCMcodeBias) {
469 required.push_back(new t_pppParam(t_pppParam::cBias1, t_prn(), t_lc::c1));
470 required.push_back(new t_pppParam(t_pppParam::cBias2, t_prn(), t_lc::c2));
471 }
472
473 // Receiver Phase Biases
474 // ---------------------
475 if ((OPT->_obsModelType == OPT->DCMphaseBias) ||
476 (OPT->_obsModelType == OPT->PPPRTK) ) {
477 required.push_back(new t_pppParam(t_pppParam::pBias1, t_prn(), t_lc::l1));
478 required.push_back(new t_pppParam(t_pppParam::pBias2, t_prn(), t_lc::l2));
479 }
480
481 // Check if all required parameters are present
482 // --------------------------------------------
483 for (unsigned ii = 0; ii < required.size(); ii++) {
484 t_pppParam* parReq = required[ii];
485
486 bool found = false;
487 for (unsigned jj = 0; jj < _params.size(); jj++) {
488 t_pppParam* parOld = _params[jj];
489 if (parOld->isEqual(parReq)) {
490 found = true;
491 break;
492 }
493 }
494 if (found) {
495 delete parReq;
496 }
497 else {
498 _params.push_back(parReq);
499 }
500 }
501
502 // Set Parameter Indices
503 // ---------------------
504 sort(_params.begin(), _params.end(), t_pppParam::sortFunction);
505
506 for (unsigned ii = 0; ii < _params.size(); ii++) {
507 t_pppParam* par = _params[ii];
508 par->setIndex(ii);
509 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
510 const t_pppSatObs* satObs = obsVector[jj];
511 if (satObs->prn() == par->prn()) {
512 par->setAmbEleSat(satObs->eleSat());
513 par->stepAmbNumEpo();
514 }
515 }
516 }
517
518 return success;
519}
520
521//
522////////////////////////////////////////////////////////////////////////////
523void t_pppParlist::printResult(const bncTime& epoTime, const SymmetricMatrix& QQ,
524 const ColumnVector& xx) const {
525
526 string epoTimeStr = string(epoTime);
527 const t_pppStation* sta = PPP_CLIENT->staRover();
528
529 LOG << endl;
530
531 t_pppParam* parX = 0;
532 t_pppParam* parY = 0;
533 t_pppParam* parZ = 0;
534 for (unsigned ii = 0; ii < _params.size(); ii++) {
535 t_pppParam* par = _params[ii];
536 if (par->type() == t_pppParam::crdX) {
537 parX = par;
538 }
539 else if (par->type() == t_pppParam::crdY) {
540 parY = par;
541 }
542 else if (par->type() == t_pppParam::crdZ) {
543 parZ = par;
544 }
545 else {
546 int ind = par->indexNew();
547 double apr = (par->type() == t_pppParam::trp) ?
548 t_tropo::delay_saast(sta->xyzApr(), M_PI/2.0) : par->x0();
549 LOG << epoTimeStr << ' ' << par->toString() << ' '
550 << setw(10) << setprecision(4) << apr << ' '
551 << showpos << setw(10) << setprecision(4) << xx[ind] << noshowpos << " +- "
552 << setw(8) << setprecision(4) << sqrt(QQ[ind][ind]);
553 if (par->type() == t_pppParam::amb) {
554 LOG << " el = " << setw(6) << setprecision(2) << par->ambEleSat() * 180.0 / M_PI
555 << " epo = " << setw(4) << par->ambNumEpo();
556 }
557 LOG << endl;
558 }
559 }
560
561 if (parX && parY && parZ) {
562
563 ColumnVector xyz(3);
564 xyz[0] = xx[parX->indexNew()];
565 xyz[1] = xx[parY->indexNew()];
566 xyz[2] = xx[parZ->indexNew()];
567
568 ColumnVector neu(3);
569 xyz2neu(sta->ellApr().data(), xyz.data(), neu.data());
570
571 SymmetricMatrix QQxyz = QQ.SymSubMatrix(1,3);
572
573 SymmetricMatrix QQneu(3);
574 covariXYZ_NEU(QQxyz, sta->ellApr().data(), QQneu);
575
576 LOG << epoTimeStr << ' ' << sta->name()
577 << " X = " << setprecision(4) << sta->xyzApr()[0] + xyz[0] << " +- "
578 << setprecision(4) << sqrt(QQxyz[0][0])
579
580 << " Y = " << setprecision(4) << sta->xyzApr()[1] + xyz[1] << " +- "
581 << setprecision(4) << sqrt(QQxyz[1][1])
582
583 << " Z = " << setprecision(4) << sta->xyzApr()[2] + xyz[2] << " +- "
584 << setprecision(4) << sqrt(QQxyz[2][2])
585
586 << " dN = " << setprecision(4) << neu[0] << " +- "
587 << setprecision(4) << sqrt(QQneu[0][0])
588
589 << " dE = " << setprecision(4) << neu[1] << " +- "
590 << setprecision(4) << sqrt(QQneu[1][1])
591
592 << " dU = " << setprecision(4) << neu[2] << " +- "
593 << setprecision(4) << sqrt(QQneu[2][2])
594
595 << endl;
596 }
597}
598
Note: See TracBrowser for help on using the repository browser.