Changeset 11083 in ntrip


Ignore:
Timestamp:
Oct 12, 2026, 3:32:43 PM (9 hours ago)
Author:
stuerze
Message:

further fixes and improvements

Location:
trunk/BNC
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/BNC/CHANGELOG.md

    r11082 r11083  
    55- 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
    66- ADDED: SIS (radial - clock, satellite-specific clock offset included) per epoch and SISRMS/SISSig per satellite and system in the SP3 comparison
     7- 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
    78- 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
    89- ADDED: Combination: clocks of AC satellites without code biases of both reference signals are combined, but excluded from the clock datum (hourly message)
    … …  
    2223- CHANGED: Replays of raw files ('--file') are reproducible, including PPP with ambiguity resolution; corrections pending at the end of a raw file are written
    2324- CHANGED: Example configurations 17 and 18 (sampling interval 5 s, ACs CNE, CHC, BKG, CAS bias file), ANTEX file igs20.atx updated
     25- 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
    2426- CHANGED: Legacy PPP_SSR_I module removed
     27- 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
    2528- 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
    2629- 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)

  • trunk/BNC/Example_Configs/21_CompareSp3.bnc

    r10968 r11083  
    107107sslClientCertPath=
    108108sslIgnoreErrors=0
     109ssrQcIntr=1 hour
     110ssrQcPath=
     111ssrQcRef=CONSENSUS
     112ssrQcStreams=ALL
    109113startTab=5
    110114statusTab=0

  • trunk/BNC/scripts/compareSP3Plot.pl

    r10126 r11083  
    1515# - at first one offset 'per epoch for all satellites' and
    1616# - secondly one offset 'per satellite for all epochs'.
     17# The signal-in-space difference SIS (radial minus clock difference, only the
     18# epoch offset removed, i.e. including the satellite offsets) is given in the
     19# last column of BNC's epoch lines.
    1720# ========================================================================
    1821# Uses
    … …  
    7073    my (%COMPARISON);
    7174    my (@array);
    72     my ($sys, $sat, $epo, $radial, $along, $out, $drange, $clk, $clkRed);
     75    my ($sys, $sat, $epo, $radial, $along, $out, $drange, $clk, $clkRed, $sis);
    7376    # -----------------------------------------------------------------------------
    7477    # create pdf for plot results
    … …  
    9598            $clk    = $array[5];#print "$clk\n";
    9699            $clkRed = $array[6];#print "$clkRed\n";
     100            $sis    = $array[8];# not present if the clock comparison is missing
    97101            if ($out !~ /^\.$/) {
    98102              push @{$COMPARISON{$sys}{$sat}{EPOCH }}, $epo;
    … …  
    102106              push @{$COMPARISON{$sys}{$sat}{CLK   }}, $clk;
    103107              push @{$COMPARISON{$sys}{$sat}{CLKRED}}, $clkRed;
     108            }
     109            if (defined $sis && $sis =~ /^-?\d/) {
     110              push @{$COMPARISON{$sys}{$sat}{EPOCH_SIS}}, $epo;
     111              push @{$COMPARISON{$sys}{$sat}{SIS      }}, $sis;
    104112            }
    105113        }
    … …  
    280288        }
    281289    }
     290    ######## SIS #####################
     291    if (("ALL" ~~ @plotTypes ) || ("SIS" ~~ @plotTypes)) {
     292        print "Plot SIGNAL-IN-SPACE DIFFERENCES (RADIAL MINUS CLOCK) ...\n";
     293        $page = $pdf->page(); $page->mediabox('A4');
     294        $headline = sprintf("Signal-in-space differences, radial minus clock, incl. satellite offsets (%s)", $inputFilename);
     295        $headline_text = $page->text;
     296        $headline_text->font($font1, 11/pt);
     297        $headline_text->translate(15/mm, 280/mm);
     298        $headline_text->text($headline);
     299        $y = $y0+$dy;
     300        #SYSTEM
     301        foreach my $key_sys (sort keys %COMPARISON) {
     302            my @sats = grep { exists $COMPARISON{$key_sys}{$_}{SIS} } sort keys %{$COMPARISON{$key_sys}};
     303            next unless @sats;
     304            print "    for $key_sys\n";
     305            my @datasets_s; # init datasets
     306            my $pngName_s  = sprintf ("%s_SIS_%s.png", $inputFilename, $key_sys);
     307            # SIS contains the satellite offsets: y range from the data
     308            my $maxAbs = 0.1;
     309            foreach my $key_sat (@sats) {
     310                my ( $min, $max ) = minmax @{$COMPARISON{$key_sys}{$key_sat}{SIS}};
     311                $maxAbs = abs($min) if abs($min) > $maxAbs;
     312                $maxAbs = abs($max) if abs($max) > $maxAbs;
     313            }
     314            my $chart_s = Chart::Gnuplot->new(
     315                          output => $pngName_s, terminal => 'png',title => "SIS (radial minus clock)",
     316                          ylabel => "SIS Difference [m]", yrange => ["-$maxAbs", "$maxAbs"],
     317                          xlabel => "Time [h]", timeaxis => 'x', xtics => {labelfmt => '%H:%M', rotate => '-270',},
     318                          legend => { position => "outside right",},
     319                          grid => 'on',);
     320            #SATELLITE
     321            foreach my $key_sat (@sats) {
     322                my $colour = sprintf("%d", (substr $key_sat, -2)+1);
     323                my $dataset_s = Chart::Gnuplot::DataSet->new(
     324                                xdata => \@{$COMPARISON{$key_sys}{$key_sat}{EPOCH_SIS}},
     325                                ydata => \@{$COMPARISON{$key_sys}{$key_sat}{SIS}},
     326                                title => "$key_sat",
     327                                timefmt => '%Y-%m-%d_%H:%M:%S',
     328                                linetype => $colour,
     329                                style => " lines ",);
     330                push (@datasets_s, $dataset_s);
     331            }
     332            $chart_s->plot2d(@datasets_s);
     333            $y = $y - $dy;
     334            if($y < 30/mm) {
     335                $page = $pdf->page(); $page->mediabox('A4');
     336                $y = $y0;
     337            }
     338            $png = $page->gfx(); die("Unable to find image file: $!") unless -e $pngName_s;
     339            $png->image($pdf->image_png($pngName_s), $x, $y, $width, $height);
     340        }
     341    }
    282342    $pdf->save();
    283343    $pdf->end();
    … …  
    295355PARAMETER:
    296356
    297       --plotTypes     possible: ALL => RAO and CLOCK; CLK => CLOCK and CLOCKRED; RAO => RAO only
    298       --logFiles      comma separated list of BNC's combination logfiles
     357      --plotTypes     possible: ALL => RAO, CLOCK and SIS; CLK => CLOCK and CLOCKRED; RAO => RAO only;
     358                      SIS => signal-in-space differences (radial minus clock) only
     359      --logFiles      comma separated list of BNC's SP3 comparison logfiles
    299360      --help          show help contents
    300361

  • trunk/BNC/src/RTCM3/RTCM3coDecoder.cpp

    r11074 r11083  
    12751275    if (bdt) {
    12761276      epoSecBds += 14.0;
    1277       if (epoSecBds > 604800.0) {
     1277      if (epoSecBds >= 604800.0) {
    12781278        epoSecBds -= 7.0*24.0*60.0*60.0;
    12791279      }
    … …  
    12891289  }
    12901290
     1291  // The epoch is given as seconds of week (or of day): move it by whole days
     1292  // to the one closest to the current time. Folding in 12 h steps (as
     1293  // before) stopped at the edge of the +-12 h window, e.g. at the week
     1294  // change a correction of Saturday 23:59:55 received after midnight ended
     1295  // up at Sunday 11:59:55.
    12911296  if (_lastTime.valid()) {
    1292     double maxDiff = 12 * 3600.0;
    1293     while (_lastTime < currentTime - maxDiff) {
    1294       _lastTime = _lastTime + maxDiff;
    1295     }
    1296     while (_lastTime > currentTime + maxDiff) {
    1297       _lastTime = _lastTime - maxDiff;
     1297    const double halfDay = 12 * 3600.0;
     1298    const double day     = 24 * 3600.0;
     1299    while (_lastTime < currentTime - halfDay) {
     1300      _lastTime = _lastTime + day;
     1301    }
     1302    while (_lastTime > currentTime + halfDay) {
     1303      _lastTime = _lastTime - day;
    12981304    }
    12991305  }

  • trunk/BNC/src/RTCM3/clock_and_orbit/clock_orbit_rtcm.cpp

    r11042 r11083  
    9797            break;
    9898          case CLOCKORBIT_SATGLONASS:
    99             E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
     99            E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
    100100            E_GLONASS_IOD(co->Sat[i].IOD)
    101101            break;
    … …  
    161161            break;
    162162          case CLOCKORBIT_SATGLONASS:
    163             E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
     163            E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
    164164            break;
    165165        }
    … …  
    221221            break;
    222222          case CLOCKORBIT_SATGLONASS:
    223             E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
     223            E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
    224224            E_GLONASS_IOD(co->Sat[i].IOD)
    225225            break;
    … …  
    294294            break;
    295295          case CLOCKORBIT_SATGLONASS:
    296             E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
     296            E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
    297297            break;
    298298        }
    … …  
    334334            break;
    335335          case CLOCKORBIT_SATGLONASS:
    336             E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
     336            E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
    337337            break;
    338338        }
    … …  
    385385            break;
    386386          case CLOCKORBIT_SATGLONASS:
    387             E_GLONASS_SATELLITE_ID(cb->Sat[i].ID)
     387            E_GLONASS_SATELLITE_ID_OLD(cb->Sat[i].ID)
    388388            break;
    389389        }
    … …  
    442442            break;
    443443          case CLOCKORBIT_SATGLONASS:
    444             E_GLONASS_SATELLITE_ID(pb->Sat[i].ID)
     444            E_GLONASS_SATELLITE_ID_OLD(pb->Sat[i].ID)
    445445            break;
    446446        }

  • trunk/BNC/src/bnchelp.html

    r11081 r11083  
    43374337  </p>
    43384338  <p>
    4339     To compare satellite clocks provided by the two files, BNC first converts coordinate differences dX,dY,dZ
    4340     into along track, out-of-plane, and radial components. It then corrects the clock differences for the radial
    4341     components
    4342     of coordinate differences. RMS values of clock differences are finally calculated after introducing at first one
    4343     offset
    4344     'per epoch for all satellites' and secondly one offset 'per satellite for all epochs'.
     4339    To compare the orbits, BNC converts the coordinate differences dX,dY,dZ into radial, along-track and out-of-plane
     4340    components, using the satellite velocity derived from neighbouring epochs.
     4341  </p>
     4342  <p>
     4343    The clocks of two products usually differ by an offset common to all satellites of a system, which changes from
     4344    epoch to epoch (clock datum), and by a constant offset per satellite (e.g. different realizations of the satellite
     4345    code biases). For each system, BNC therefore estimates from all clock differences in a least-squares adjustment
     4346    one offset 'per epoch for all satellites' and one offset 'per satellite for all epochs' (column 'Offset'), with
     4347    the sum of the satellite offsets constrained to zero. The clock residuals of this adjustment ('clk') are then
     4348    corrected for the radial component of the orbit differences ('clkRed'). The signal-in-space difference ('SIS',
     4349    radial orbit difference minus clock difference) removes only the epoch offset, i.e. it contains the satellite
     4350    offsets.
     4351  </p>
     4352  <p>
     4353    The epoch offsets are eliminated from the normal equations before the satellite offsets are solved for; the
     4354    result is identical to that of the complete adjustment. This keeps the comparison fast also for long or densely
     4355    sampled files: comparing two daily files with 30-second sampling (2880 epochs) takes a few seconds.
    43454356  </p>
    43464357  <p><img src="IMG/Figure15.png" width=800 /></p>

  • trunk/BNC/src/combination/bnccomb.cpp

    r11079 r11083  
    15251525    // 30 min and are valid for ±900 s around the reference time (TOC).
    15261526    // t_ephGlo::position() uses Runge-Kutta integration from TOC; beyond
    1527     // ±900 s the integration error grows unchecked and propagates into both
    1528     // the combined clock estimate and the SP3 position.
     1527    // that the integration error grows and propagates into both the
     1528    // combined clock estimate and the SP3 position. A margin of 60 s is
     1529    // allowed, as ACs switch to the next ephemeris a few seconds before its
     1530    // window begins (otherwise satellites drop out at every switch).
    15291531    if (corr->_eph->prn().system() == 'R') {
    15301532      double age = epoTime - corr->_eph->TOC();
    1531       if (fabs(age) > 900.0) {
     1533      if (fabs(age) > 960.0) {
    15321534        emit newMessage(QString().asprintf(
    15331535          "bncComb: skip %s, GLONASS eph age %.0fs",

  • trunk/BNC/src/orbComp/sp3Comp.cpp

    r11076 r11083  
    178178  }
    179179
    180   int nPar = epochs.size() + clkSats.size();
    181   SymmetricMatrix NN(nPar); NN = 0.0;
    182   ColumnVector    bb(nPar); bb = 0.0;
    183 
    184   // Create Matrix A'A and vector b
    185   // ------------------------------
    186   for (unsigned ie = 0; ie < epochs.size(); ie++) {
     180  // Least-squares adjustment of an epoch-specific offset (common to all
     181  // satellites) and a satellite-specific offset (common to all epochs), the
     182  // sum of the satellite offsets constrained to zero (weight 1). The normal
     183  // equations of the epoch offsets are diagonal, so these offsets are
     184  // eliminated (Schur complement) and only the satellite offsets are solved
     185  // for; the epoch offsets follow from them. The solution is identical to
     186  // that of the full system, but the effort grows only linearly with the
     187  // number of epochs (before: cubic).
     188  // -----------------------------------------------------------------------
     189  int nEpo = epochs.size();
     190  int nSat = clkSats.size();
     191  int nPar = nEpo + nSat;
     192
     193  vector<vector<int> >    satIdx(nEpo);   // satellites (index) of each epoch
     194  ColumnVector    nEp(nEpo); nEp = 0.0;   // number of satellites per epoch
     195  ColumnVector    bEp(nEpo); bEp = 0.0;   // sum of the differences per epoch
     196  SymmetricMatrix NS(nSat);  NS  = 1.0;   // regularization: sum of satellite offsets
     197  ColumnVector    bS(nSat);  bS  = 0.0;
     198  for (int ie = 0; ie < nEpo; ie++) {
    187199    const map<t_prn, double>& dc = epochs[ie]->_dc;
    188     Matrix       AA(dc.size(), nPar); AA = 0.0;
    189     ColumnVector ll(dc.size());       ll = 0.0;
    190     map<t_prn, double>::const_iterator it;
    191     int ii = -1;
    192     for (it = dc.begin(); it != dc.end(); it++) {
    193       const t_prn& prn = it->first;
    194       if (satIndex(clkSats, prn) != -1) {
    195         ++ii;
    196         int index = epochs.size() + satIndex(clkSats, prn);
    197         AA[ii][ie]    = 1.0; // epoch-specfic offset (common for all satellites)
    198         AA[ii][index] = 1.0; // satellite-specific offset (common for all epochs)
    199         ll[ii]        = it->second;
    200       }
    201     }
    202     SymmetricMatrix dN; dN << AA.t() * AA;
    203     NN += dN;
    204     bb += AA.t() * ll;
    205   }
    206 
    207   // Regularize NN
    208   // -------------
    209   RowVector HH(nPar);
    210   HH.columns(1, epochs.size())      = 0.0;
    211   HH.columns(epochs.size()+1, nPar) = 1.0;
    212   SymmetricMatrix dN; dN << HH.t() * HH;
    213   NN += dN;
    214 
    215   // Estimate Parameters
    216   // -------------------
    217   ColumnVector xx = NN.i() * bb;
     200    for (map<t_prn, double>::const_iterator it = dc.begin(); it != dc.end(); it++) {
     201      int is = satIndex(clkSats, it->first);
     202      if (is != -1) {
     203        satIdx[ie].push_back(is);
     204        nEp[ie]     += 1.0;
     205        bEp[ie]     += it->second;
     206        NS(is+1, is+1) += 1.0;
     207        bS[is]      += it->second;
     208      }
     209    }
     210  }
     211
     212  // Eliminate the epoch offsets
     213  // ---------------------------
     214  for (int ie = 0; ie < nEpo; ie++) {
     215    const vector<int>& idx = satIdx[ie];
     216    for (unsigned ia = 0; ia < idx.size(); ia++) {
     217      bS[idx[ia]] -= bEp[ie] / nEp[ie];
     218      for (unsigned ib = 0; ib <= ia; ib++) {
     219        NS(idx[ia]+1, idx[ib]+1) -= 1.0 / nEp[ie];
     220      }
     221    }
     222  }
     223
     224  // Estimate Parameters (epoch offsets first, then satellite offsets)
     225  // ---------------------------------------------------------------
     226  ColumnVector yy = NS.i() * bS;
     227  ColumnVector xx(nPar);
     228  for (int ie = 0; ie < nEpo; ie++) {
     229    double sum = bEp[ie];
     230    for (unsigned ia = 0; ia < satIdx[ie].size(); ia++) {
     231      sum -= yy[satIdx[ie][ia]];
     232    }
     233    xx[ie] = sum / nEp[ie];
     234  }
     235  for (int is = 0; is < nSat; is++) {
     236    xx[nEpo + is] = yy[is];
     237  }
    218238
    219239  // Compute clock residuals

  • trunk/BNC/tests/README.md

    r11075 r11083  
    7171  received stream again: `upload_rtnet` (IGS-SSR encoding, SP3/CLK/BIA files,
    7272  decoded `.ssr`), `upload_eph` (RTCM3 ephemerides to RINEX nav) and
    73   `upload_raw` (raw forwarding to RINEX nav).
     73  `upload_raw` (raw forwarding to RINEX nav). `upload_rtnet_rtcm` encodes the same input to RTCM-SSR
     74  (message numbers 1057-1068, 1240-1270) and decodes it again.
    7475* SSR decoding (`ssr_*`): 3 min slices of real streams replayed through
    7576  `03_BrdcCorr.bnc`. They compare the `.ssr` files and the `BNC_DEBUG_SSR`
    … …  
    7980  * `ssr_rtcm_SSRA00CNE0`, `ssr_rtcm_WHU1_OSBC`: RTCM 1059-1270; these test
    8081    whichever class `RTCM3Decoder.cpp` selects for non-4076 SSR.
     82* `ssr_weekchange`: 2 min across the GPS week change (SSRA00CNE1, SSRA00CHC1);
     83  corrections of Saturday received after midnight must keep their epochs.
    8184* `ssr_qc`, `ssr_qc_select`: SSR QC (31_SsrQC, one file per stream, 5 min interval;
    8285  all streams resp. `ssrQcStreams SSRA00CNE1,SSRA00GMV0`) of the 7 AC streams of
Note: See TracChangeset for help on using the changeset viewer.