Changeset 11083 in ntrip for trunk/BNC/src

Timestamp:

Oct 12, 2026, 3:32:43 PM (2 weeks ago)

Author:

stuerze

Message:

further fixes and improvements

Location:

trunk/BNC/src

Files:

• RTCM3/RTCM3coDecoder.cpp (modified) (2 diffs)
• RTCM3/clock_and_orbit/clock_orbit_rtcm.cpp (modified) (7 diffs)
• bnchelp.html (modified) (1 diff)
• combination/bnccomb.cpp (modified) (1 diff)
• orbComp/sp3Comp.cpp (modified) (1 diff)

trunk/BNC/src/RTCM3/RTCM3coDecoder.cpp

@@ -1275 +1275 @@
     if (bdt) {
       epoSecBds += 14.0;
-      if (epoSecBds > 604800.0) {
+      if (epoSecBds >= 604800.0) {
         epoSecBds -= 7.0*24.0*60.0*60.0;
       }
@@ -1289 +1289 @@
   }
 
+  // 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;
     }
   }

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

@@ -97 +97 @@
             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 +161 @@
             break;
           case CLOCKORBIT_SATGLONASS:
-            E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+            E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
             break;
         }
@@ -221 +221 @@
             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 +294 @@
             break;
           case CLOCKORBIT_SATGLONASS:
-            E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+            E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
             break;
         }
@@ -334 +334 @@
             break;
           case CLOCKORBIT_SATGLONASS:
-            E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+            E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
             break;
         }
@@ -385 +385 @@
             break;
           case CLOCKORBIT_SATGLONASS:
-            E_GLONASS_SATELLITE_ID(cb->Sat[i].ID)
+            E_GLONASS_SATELLITE_ID_OLD(cb->Sat[i].ID)
             break;
         }
@@ -442 +442 @@
             break;
           case CLOCKORBIT_SATGLONASS:
-            E_GLONASS_SATELLITE_ID(pb->Sat[i].ID)
+            E_GLONASS_SATELLITE_ID_OLD(pb->Sat[i].ID)
             break;
         }

trunk/BNC/src/bnchelp.html

@@ -4337 +4337 @@
   </p>
   <p>
-    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.
+  </p>
+  <p>
+    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.
+  </p>
+  <p>
+    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.
   </p>
   <p><img src="IMG/Figure15.png" width=800 /></p>

trunk/BNC/src/combination/bnccomb.cpp

@@ -1525 +1525 @@
     // 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",

trunk/BNC/src/orbComp/sp3Comp.cpp

@@ -178 +178 @@
   }
 
-  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<vector<int> >    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<t_prn, double>& dc = epochs[ie]->_dc;
-    Matrix       AA(dc.size(), nPar); AA = 0.0;
-    ColumnVector ll(dc.size());       ll = 0.0;
-    map<t_prn, double>::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<t_prn, double>::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<int>& 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