Index: /trunk/BNC/CHANGELOG.md
===================================================================
--- /trunk/BNC/CHANGELOG.md (revision 11082)
+++ /trunk/BNC/CHANGELOG.md (revision 11083)
@@ -5,4 +5,5 @@
- ADDED: Stream comparison in the SSR QC: orbits (radial, along, cross), clocks referred to the reference signals and signal-in-space differences (SIS = radial - clock) of each stream compared with the consensus of the other ACs (one vote per AC, own AC left out, combination products not voting) or with a selected stream (e.g. 'INTERNAL'), as RMS and standard deviation per satellite and system
- ADDED: SIS (radial - clock, satellite-specific clock offset included) per epoch and SISRMS/SISSig per satellite and system in the SP3 comparison
+- CHANGED: SP3 comparison: the clock offsets are estimated with the epoch offsets eliminated from the normal equations (identical result); a daily file with 30-second sampling is now compared in seconds instead of hours
- ADDED: Combination: check of the consistency of each AC's code biases with its clocks; an AC whose satellite offsets follow the ionosphere-free combination of its reference signal biases is reported and excluded from the clock datum
- ADDED: Combination: clocks of AC satellites without code biases of both reference signals are combined, but excluded from the clock datum (hourly message)
@@ -22,5 +23,7 @@
- CHANGED: Replays of raw files ('--file') are reproducible, including PPP with ambiguity resolution; corrections pending at the end of a raw file are written
- CHANGED: Example configurations 17 and 18 (sampling interval 5 s, ACs CNE, CHC, BKG, CAS bias file), ANTEX file igs20.atx updated
+- CHANGED: Combination: GLONASS corrections are accepted up to 960 s (before 900 s) from the reference time of their broadcast ephemeris, as ACs switch to the next ephemeris shortly before its validity begins
- CHANGED: Legacy PPP_SSR_I module removed
+- FIXED: SSR corrections of the old GPS week received after the week change (e.g. Saturday 23:59:55 received at 00:00:03) were assigned an epoch 12 h later (Sunday 11:59:55); affected all streams once per week
- FIXED: IGS-SSR epoch time (IDF003) of GLONASS and BDS is GPS time; decoded corrections were placed 18 s (GLONASS) and 14 s (BDS) late, and the encoder wrote the wrong convention
- FIXED: Combination: sign of the ionosphere-free combination of the SSR code biases used to refer the AC clocks to the reference signals (SSR and SINEX sign conventions differ)
Index: /trunk/BNC/Example_Configs/21_CompareSp3.bnc
===================================================================
--- /trunk/BNC/Example_Configs/21_CompareSp3.bnc (revision 11082)
+++ /trunk/BNC/Example_Configs/21_CompareSp3.bnc (revision 11083)
@@ -107,4 +107,8 @@
sslClientCertPath=
sslIgnoreErrors=0
+ssrQcIntr=1 hour
+ssrQcPath=
+ssrQcRef=CONSENSUS
+ssrQcStreams=ALL
startTab=5
statusTab=0
Index: /trunk/BNC/scripts/compareSP3Plot.pl
===================================================================
--- /trunk/BNC/scripts/compareSP3Plot.pl (revision 11082)
+++ /trunk/BNC/scripts/compareSP3Plot.pl (revision 11083)
@@ -15,4 +15,7 @@
# - at first one offset 'per epoch for all satellites' and
# - secondly one offset 'per satellite for all epochs'.
+# The signal-in-space difference SIS (radial minus clock difference, only the
+# epoch offset removed, i.e. including the satellite offsets) is given in the
+# last column of BNC's epoch lines.
# ========================================================================
# Uses
@@ -70,5 +73,5 @@
my (%COMPARISON);
my (@array);
- my ($sys, $sat, $epo, $radial, $along, $out, $drange, $clk, $clkRed);
+ my ($sys, $sat, $epo, $radial, $along, $out, $drange, $clk, $clkRed, $sis);
# -----------------------------------------------------------------------------
# create pdf for plot results
@@ -95,4 +98,5 @@
$clk = $array[5];#print "$clk\n";
$clkRed = $array[6];#print "$clkRed\n";
+ $sis = $array[8];# not present if the clock comparison is missing
if ($out !~ /^\.$/) {
push @{$COMPARISON{$sys}{$sat}{EPOCH }}, $epo;
@@ -102,4 +106,8 @@
push @{$COMPARISON{$sys}{$sat}{CLK }}, $clk;
push @{$COMPARISON{$sys}{$sat}{CLKRED}}, $clkRed;
+ }
+ if (defined $sis && $sis =~ /^-?\d/) {
+ push @{$COMPARISON{$sys}{$sat}{EPOCH_SIS}}, $epo;
+ push @{$COMPARISON{$sys}{$sat}{SIS }}, $sis;
}
}
@@ -280,4 +288,56 @@
}
}
+ ######## SIS #####################
+ if (("ALL" ~~ @plotTypes ) || ("SIS" ~~ @plotTypes)) {
+ print "Plot SIGNAL-IN-SPACE DIFFERENCES (RADIAL MINUS CLOCK) ...\n";
+ $page = $pdf->page(); $page->mediabox('A4');
+ $headline = sprintf("Signal-in-space differences, radial minus clock, incl. satellite offsets (%s)", $inputFilename);
+ $headline_text = $page->text;
+ $headline_text->font($font1, 11/pt);
+ $headline_text->translate(15/mm, 280/mm);
+ $headline_text->text($headline);
+ $y = $y0+$dy;
+ #SYSTEM
+ foreach my $key_sys (sort keys %COMPARISON) {
+ my @sats = grep { exists $COMPARISON{$key_sys}{$_}{SIS} } sort keys %{$COMPARISON{$key_sys}};
+ next unless @sats;
+ print " for $key_sys\n";
+ my @datasets_s; # init datasets
+ my $pngName_s = sprintf ("%s_SIS_%s.png", $inputFilename, $key_sys);
+ # SIS contains the satellite offsets: y range from the data
+ my $maxAbs = 0.1;
+ foreach my $key_sat (@sats) {
+ my ( $min, $max ) = minmax @{$COMPARISON{$key_sys}{$key_sat}{SIS}};
+ $maxAbs = abs($min) if abs($min) > $maxAbs;
+ $maxAbs = abs($max) if abs($max) > $maxAbs;
+ }
+ my $chart_s = Chart::Gnuplot->new(
+ output => $pngName_s, terminal => 'png',title => "SIS (radial minus clock)",
+ ylabel => "SIS Difference [m]", yrange => ["-$maxAbs", "$maxAbs"],
+ xlabel => "Time [h]", timeaxis => 'x', xtics => {labelfmt => '%H:%M', rotate => '-270',},
+ legend => { position => "outside right",},
+ grid => 'on',);
+ #SATELLITE
+ foreach my $key_sat (@sats) {
+ my $colour = sprintf("%d", (substr $key_sat, -2)+1);
+ my $dataset_s = Chart::Gnuplot::DataSet->new(
+ xdata => \@{$COMPARISON{$key_sys}{$key_sat}{EPOCH_SIS}},
+ ydata => \@{$COMPARISON{$key_sys}{$key_sat}{SIS}},
+ title => "$key_sat",
+ timefmt => '%Y-%m-%d_%H:%M:%S',
+ linetype => $colour,
+ style => " lines ",);
+ push (@datasets_s, $dataset_s);
+ }
+ $chart_s->plot2d(@datasets_s);
+ $y = $y - $dy;
+ if($y < 30/mm) {
+ $page = $pdf->page(); $page->mediabox('A4');
+ $y = $y0;
+ }
+ $png = $page->gfx(); die("Unable to find image file: $!") unless -e $pngName_s;
+ $png->image($pdf->image_png($pngName_s), $x, $y, $width, $height);
+ }
+ }
$pdf->save();
$pdf->end();
@@ -295,6 +355,7 @@
PARAMETER:
- --plotTypes possible: ALL => RAO and CLOCK; CLK => CLOCK and CLOCKRED; RAO => RAO only
- --logFiles comma separated list of BNC's combination logfiles
+ --plotTypes possible: ALL => RAO, CLOCK and SIS; CLK => CLOCK and CLOCKRED; RAO => RAO only;
+ SIS => signal-in-space differences (radial minus clock) only
+ --logFiles comma separated list of BNC's SP3 comparison logfiles
--help show help contents
Index: /trunk/BNC/src/RTCM3/RTCM3coDecoder.cpp
===================================================================
--- /trunk/BNC/src/RTCM3/RTCM3coDecoder.cpp (revision 11082)
+++ /trunk/BNC/src/RTCM3/RTCM3coDecoder.cpp (revision 11083)
@@ -1275,5 +1275,5 @@
if (bdt) {
epoSecBds += 14.0;
- if (epoSecBds > 604800.0) {
+ if (epoSecBds >= 604800.0) {
epoSecBds -= 7.0*24.0*60.0*60.0;
}
@@ -1289,11 +1289,17 @@
}
+ // The epoch is given as seconds of week (or of day): move it by whole days
+ // to the one closest to the current time. Folding in 12 h steps (as
+ // before) stopped at the edge of the +-12 h window, e.g. at the week
+ // change a correction of Saturday 23:59:55 received after midnight ended
+ // up at Sunday 11:59:55.
if (_lastTime.valid()) {
- double maxDiff = 12 * 3600.0;
- while (_lastTime < currentTime - maxDiff) {
- _lastTime = _lastTime + maxDiff;
- }
- while (_lastTime > currentTime + maxDiff) {
- _lastTime = _lastTime - maxDiff;
+ const double halfDay = 12 * 3600.0;
+ const double day = 24 * 3600.0;
+ while (_lastTime < currentTime - halfDay) {
+ _lastTime = _lastTime + day;
+ }
+ while (_lastTime > currentTime + halfDay) {
+ _lastTime = _lastTime - day;
}
}
Index: /trunk/BNC/src/RTCM3/clock_and_orbit/clock_orbit_rtcm.cpp
===================================================================
--- /trunk/BNC/src/RTCM3/clock_and_orbit/clock_orbit_rtcm.cpp (revision 11082)
+++ /trunk/BNC/src/RTCM3/clock_and_orbit/clock_orbit_rtcm.cpp (revision 11083)
@@ -97,5 +97,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
E_GLONASS_IOD(co->Sat[i].IOD)
break;
@@ -161,5 +161,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
break;
}
@@ -221,5 +221,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
E_GLONASS_IOD(co->Sat[i].IOD)
break;
@@ -294,5 +294,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
break;
}
@@ -334,5 +334,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
break;
}
@@ -385,5 +385,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(cb->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(cb->Sat[i].ID)
break;
}
@@ -442,5 +442,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(pb->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(pb->Sat[i].ID)
break;
}
Index: /trunk/BNC/src/bnchelp.html
===================================================================
--- /trunk/BNC/src/bnchelp.html (revision 11082)
+++ /trunk/BNC/src/bnchelp.html (revision 11083)
@@ -4337,10 +4337,21 @@
- To compare satellite clocks provided by the two files, BNC first converts coordinate differences dX,dY,dZ
- into along track, out-of-plane, and radial components. It then corrects the clock differences for the radial
- components
- of coordinate differences. RMS values of clock differences are finally calculated after introducing at first one
- offset
- 'per epoch for all satellites' and secondly one offset 'per satellite for all epochs'.
+ To compare the orbits, BNC converts the coordinate differences dX,dY,dZ into radial, along-track and out-of-plane
+ components, using the satellite velocity derived from neighbouring epochs.
+
+
+ The clocks of two products usually differ by an offset common to all satellites of a system, which changes from
+ epoch to epoch (clock datum), and by a constant offset per satellite (e.g. different realizations of the satellite
+ code biases). For each system, BNC therefore estimates from all clock differences in a least-squares adjustment
+ one offset 'per epoch for all satellites' and one offset 'per satellite for all epochs' (column 'Offset'), with
+ the sum of the satellite offsets constrained to zero. The clock residuals of this adjustment ('clk') are then
+ corrected for the radial component of the orbit differences ('clkRed'). The signal-in-space difference ('SIS',
+ radial orbit difference minus clock difference) removes only the epoch offset, i.e. it contains the satellite
+ offsets.
+
+
+ The epoch offsets are eliminated from the normal equations before the satellite offsets are solved for; the
+ result is identical to that of the complete adjustment. This keeps the comparison fast also for long or densely
+ sampled files: comparing two daily files with 30-second sampling (2880 epochs) takes a few seconds.
Index: /trunk/BNC/src/combination/bnccomb.cpp
===================================================================
--- /trunk/BNC/src/combination/bnccomb.cpp (revision 11082)
+++ /trunk/BNC/src/combination/bnccomb.cpp (revision 11083)
@@ -1525,9 +1525,11 @@
// 30 min and are valid for ±900 s around the reference time (TOC).
// t_ephGlo::position() uses Runge-Kutta integration from TOC; beyond
- // ±900 s the integration error grows unchecked and propagates into both
- // the combined clock estimate and the SP3 position.
+ // that the integration error grows and propagates into both the
+ // combined clock estimate and the SP3 position. A margin of 60 s is
+ // allowed, as ACs switch to the next ephemeris a few seconds before its
+ // window begins (otherwise satellites drop out at every switch).
if (corr->_eph->prn().system() == 'R') {
double age = epoTime - corr->_eph->TOC();
- if (fabs(age) > 900.0) {
+ if (fabs(age) > 960.0) {
emit newMessage(QString().asprintf(
"bncComb: skip %s, GLONASS eph age %.0fs",
Index: /trunk/BNC/src/orbComp/sp3Comp.cpp
===================================================================
--- /trunk/BNC/src/orbComp/sp3Comp.cpp (revision 11082)
+++ /trunk/BNC/src/orbComp/sp3Comp.cpp (revision 11083)
@@ -178,42 +178,62 @@
}
- int nPar = epochs.size() + clkSats.size();
- SymmetricMatrix NN(nPar); NN = 0.0;
- ColumnVector bb(nPar); bb = 0.0;
-
- // Create Matrix A'A and vector b
- // ------------------------------
- for (unsigned ie = 0; ie < epochs.size(); ie++) {
+ // Least-squares adjustment of an epoch-specific offset (common to all
+ // satellites) and a satellite-specific offset (common to all epochs), the
+ // sum of the satellite offsets constrained to zero (weight 1). The normal
+ // equations of the epoch offsets are diagonal, so these offsets are
+ // eliminated (Schur complement) and only the satellite offsets are solved
+ // for; the epoch offsets follow from them. The solution is identical to
+ // that of the full system, but the effort grows only linearly with the
+ // number of epochs (before: cubic).
+ // -----------------------------------------------------------------------
+ int nEpo = epochs.size();
+ int nSat = clkSats.size();
+ int nPar = nEpo + nSat;
+
+ vector > satIdx(nEpo); // satellites (index) of each epoch
+ ColumnVector nEp(nEpo); nEp = 0.0; // number of satellites per epoch
+ ColumnVector bEp(nEpo); bEp = 0.0; // sum of the differences per epoch
+ SymmetricMatrix NS(nSat); NS = 1.0; // regularization: sum of satellite offsets
+ ColumnVector bS(nSat); bS = 0.0;
+ for (int ie = 0; ie < nEpo; ie++) {
const map& dc = epochs[ie]->_dc;
- Matrix AA(dc.size(), nPar); AA = 0.0;
- ColumnVector ll(dc.size()); ll = 0.0;
- map::const_iterator it;
- int ii = -1;
- for (it = dc.begin(); it != dc.end(); it++) {
- const t_prn& prn = it->first;
- if (satIndex(clkSats, prn) != -1) {
- ++ii;
- int index = epochs.size() + satIndex(clkSats, prn);
- AA[ii][ie] = 1.0; // epoch-specfic offset (common for all satellites)
- AA[ii][index] = 1.0; // satellite-specific offset (common for all epochs)
- ll[ii] = it->second;
- }
- }
- SymmetricMatrix dN; dN << AA.t() * AA;
- NN += dN;
- bb += AA.t() * ll;
- }
-
- // Regularize NN
- // -------------
- RowVector HH(nPar);
- HH.columns(1, epochs.size()) = 0.0;
- HH.columns(epochs.size()+1, nPar) = 1.0;
- SymmetricMatrix dN; dN << HH.t() * HH;
- NN += dN;
-
- // Estimate Parameters
- // -------------------
- ColumnVector xx = NN.i() * bb;
+ for (map::const_iterator it = dc.begin(); it != dc.end(); it++) {
+ int is = satIndex(clkSats, it->first);
+ if (is != -1) {
+ satIdx[ie].push_back(is);
+ nEp[ie] += 1.0;
+ bEp[ie] += it->second;
+ NS(is+1, is+1) += 1.0;
+ bS[is] += it->second;
+ }
+ }
+ }
+
+ // Eliminate the epoch offsets
+ // ---------------------------
+ for (int ie = 0; ie < nEpo; ie++) {
+ const vector& idx = satIdx[ie];
+ for (unsigned ia = 0; ia < idx.size(); ia++) {
+ bS[idx[ia]] -= bEp[ie] / nEp[ie];
+ for (unsigned ib = 0; ib <= ia; ib++) {
+ NS(idx[ia]+1, idx[ib]+1) -= 1.0 / nEp[ie];
+ }
+ }
+ }
+
+ // Estimate Parameters (epoch offsets first, then satellite offsets)
+ // ---------------------------------------------------------------
+ ColumnVector yy = NS.i() * bS;
+ ColumnVector xx(nPar);
+ for (int ie = 0; ie < nEpo; ie++) {
+ double sum = bEp[ie];
+ for (unsigned ia = 0; ia < satIdx[ie].size(); ia++) {
+ sum -= yy[satIdx[ie][ia]];
+ }
+ xx[ie] = sum / nEp[ie];
+ }
+ for (int is = 0; is < nSat; is++) {
+ xx[nEpo + is] = yy[is];
+ }
// Compute clock residuals
Index: /trunk/BNC/tests/README.md
===================================================================
--- /trunk/BNC/tests/README.md (revision 11082)
+++ /trunk/BNC/tests/README.md (revision 11083)
@@ -71,5 +71,6 @@
received stream again: `upload_rtnet` (IGS-SSR encoding, SP3/CLK/BIA files,
decoded `.ssr`), `upload_eph` (RTCM3 ephemerides to RINEX nav) and
- `upload_raw` (raw forwarding to RINEX nav).
+ `upload_raw` (raw forwarding to RINEX nav). `upload_rtnet_rtcm` encodes the same input to RTCM-SSR
+ (message numbers 1057-1068, 1240-1270) and decodes it again.
* SSR decoding (`ssr_*`): 3 min slices of real streams replayed through
`03_BrdcCorr.bnc`. They compare the `.ssr` files and the `BNC_DEBUG_SSR`
@@ -79,4 +80,6 @@
* `ssr_rtcm_SSRA00CNE0`, `ssr_rtcm_WHU1_OSBC`: RTCM 1059-1270; these test
whichever class `RTCM3Decoder.cpp` selects for non-4076 SSR.
+* `ssr_weekchange`: 2 min across the GPS week change (SSRA00CNE1, SSRA00CHC1);
+ corrections of Saturday received after midnight must keep their epochs.
* `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