mirrors/gnss-sdr
1
0
Fork 0
mirror of https://github.com/gnss-sdr/gnss-sdr synced 2026-10-17 04:52:36 +00:00
Code Issues Releases Wiki Activity

Let PVT.iono_model decide band combination in the single point solution

Browse Source
This commit is contained in:
Carles Fernandez
2026-09-10 12:20:38 +02:00
committed by Carles Fernandez
parent 2622d2666e
commit d84069121a
6 changed files with 298 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/CHANGELOG.md
+10
View File
@@ -79,6 +79,16 @@ All notable changes to GNSS-SDR will be documented in this file.
E5b and E6 observations always keep the clock reference of their configured
service. The reverse fallback (E1 using F/NAV when I/NAV is stale) requires
`use_unhealthy_sats=true`, since F/NAV carries no E1B health information.
- The PVT iono model now decides the pseudorange model for every satellite of a
system: only `PVT.iono_model=Iono-Free-LC` combines two bands, and any other
model uses the first band alone with its TGD/BGD (second band alone only when
the first is missing). Previously a satellite with two bands in the record was
silently switched to the iono-free combination while single-band satellites of
the same system kept the single-frequency model, mixing two clock references
within one solve and folding the receiver's uncalibrated inter-band delay
(e.g. differing input-filter group delays) into the solution, which could make
the chi-square test reject every epoch. The ISC-aware GPS L1/L5 combination is
now applied in `Iono-Free-LC` mode.
- Fixed Galileo single-frequency broadcast group-delay corrections in SPP and
PPP, selecting the E1/E5a or E1/E5b BGD from the active navigation service and
applying the ICD frequency scaling to E5a and E5b observations.
src/algorithms/libs/rtklib/rtklib_pntpos.cc
+27 -45
View File
@@ -411,7 +411,10 @@ double prange(const obsd_t *obs, const nav_t *nav, const double *azel,
int tgd_index = 0;
if (sys == SYS_GAL)
{
const unsigned char observation_code = obs->code[j] != CODE_NONE ? obs->code[j] : obs->code[i];
/* the BGD belongs to the band the returned range is
referenced to: the first band whenever it carries a
measurement, the second band only when it is alone */
const unsigned char observation_code = obs->code[i] != CODE_NONE ? obs->code[i] : obs->code[j];
tgd_index = galileo_bgd_index(observation_code, obs->sat, nav);
uses_galileo_bgd = tgd_index >= 0;
}
@@ -450,19 +453,37 @@ double prange(const obsd_t *obs, const nav_t *nav, const double *azel,
/* the DNAV clock is referenced to B3I: remove TGD1 from
the B1I pseudorange before combining (BDS SIS ICD) */
P1 -= P1_P2;
PC = (gamma_ * P1 - P2) / (gamma_ - 1.0);
}
else if ((sys == SYS_GPS || sys == SYS_QZS) && obs->code[j] == CODE_L5X)
{
/* L1 C/A + L5 dual-frequency correction, IS-GPS-705 20.3.3.3.1.2.2:
the gamma-weighted L1 term carries ISC_L1CA, not ISC_L5I5 */
PC = (P2 + ISCl5i - gamma_ * (P1 + ISCl1)) / (1.0 - gamma_) - P1_P2;
}
else
{
PC = (gamma_ * P1 - P2) / (gamma_ - 1.0);
}
PC = (gamma_ * P1 - P2) / (gamma_ - 1.0);
*iono_scale = 0.0;
}
////////////////////////////////////////////
else
{ /* single-frequency */
{
/* single-frequency: the iono model selected by the user applies to
every satellite of a system, so a second band present in the
record is never combined here (that is IONOOPT_IFLC's job). The
first band is used alone whenever it carries a measurement; the
second band is used alone only when the first one is missing.
Combining per satellite would mix two clock references within
one system as soon as any satellite lacks the second band, and
would also fold the receiver's uncalibrated inter-band delay
into the solution. */
if (obs->code[i] == CODE_NONE && (i == j || obs->code[j] == CODE_NONE))
{
return 0.0;
}
if (obs->code[i] != CODE_NONE && (i == j || obs->code[j] == CODE_NONE))
if (obs->code[i] != CODE_NONE)
{
if (sys == SYS_BDS && (obs->code[i] == CODE_L6I || obs->code[i] == CODE_L6Q))
{
@@ -516,40 +537,6 @@ double prange(const obsd_t *obs, const nav_t *nav, const double *azel,
PC = P2 - gamma_ * P1_P2;
}
}
/* dual-frequency */
else if (sys == SYS_GPS || sys == SYS_QZS) /* L1 + L2 */
{
if (obs->code[j] == CODE_L2S) /* L1 + L2 */
{
P1 += P1_C1; /* C1->P1 */
/* L1 pseudorange with LNAV clock: (dtSV)_L1 = dtSV - TGD,
see IS-GPS-200 20.3.3.3.3.2 */
PC = P1 - P1_P2;
// CNAV dual-frequency alternative (IS-GPS-200 30.3.3.3.1.1.2):
// PC = (P2 + ISCl2 - gamma_ * (P1 + ISCl1)) / (1.0 - gamma_) - P1_P2;
}
else if (obs->code[j] == CODE_L5X) /* L1 + L5 */
{
P1 += P1_C1; /* C1->P1 */
/* L1 C/A + L5 dual-frequency correction, IS-GPS-705 20.3.3.3.1.2.2:
the gamma-weighted L1 term carries ISC_L1CA, not ISC_L5I5 */
PC = (P2 + ISCl5i - gamma_ * (P1 + ISCl1)) / (1.0 - gamma_) - P1_P2;
*iono_scale = 0.0;
}
}
else if (sys == SYS_GAL || sys == SYS_GLO || sys == SYS_BDS) /* E1 + E5a / B1I + B3I */
{
P1 += P1_C1;
P2 += P2_C2;
if (sys == SYS_BDS)
{
/* the DNAV clock is referenced to B3I: remove TGD1
from the B1I pseudorange before combining */
P1 -= P1_P2;
}
PC = (gamma_ * P1 - P2) / (gamma_ - 1.0);
*iono_scale = 0.0;
}
}
if (opt->sateph == EPHOPT_SBAS)
{
@@ -852,12 +839,7 @@ int rescode(int iter, const obsd_t *obs, int n, const double *rs,
(*ns)++;
/* error variance */
double vmeasure = varerr(opt, &obs[i], azel[1 + i * 2], sys);
if (iono_scale == 0.0 && opt->ionoopt != IONOOPT_IFLC)
{
/* same noise amplification varerr applies for IONOOPT_IFLC */
vmeasure *= std::pow(3.0, 2.0);
}
const double vmeasure = varerr(opt, &obs[i], azel[1 + i * 2], sys);
var[nv++] = vmeasure + vare[i] + vmeas + vion + vtrp;
trace(4, "sat=%2d azel=%5.1f %4.1f res=%7.3f sig=%5.3f\n", obs[i].sat,
tests/test_main.cc
+1
View File
@@ -170,6 +170,7 @@ private:
#include "unit-tests/signal-processing-blocks/osnma/osnma_msg_receiver_test.cc"
#include "unit-tests/signal-processing-blocks/pvt/bds_tgd_iono_test.cc"
#include "unit-tests/signal-processing-blocks/pvt/galileo_nav_fallback_test.cc"
#include "unit-tests/signal-processing-blocks/pvt/galileo_e1_bgd_test.cc"
#include "unit-tests/signal-processing-blocks/pvt/geohash_test.cc"
#include "unit-tests/signal-processing-blocks/pvt/gps_tgd_isc_test.cc"
#include "unit-tests/signal-processing-blocks/pvt/nmea_printer_test.cc"
tests/unit-tests/signal-processing-blocks/pvt/bds_tgd_iono_test.cc
+7 -8
View File
@@ -128,7 +128,7 @@ TEST(BdsTgdIonoTest, SingleFrequencyB3IHasNoTgdAndDefersIonoScalingToRescode)
}
TEST(BdsTgdIonoTest, DualFrequencyB1IB3IRemovesTgd1BeforeCombining)
TEST(BdsTgdIonoTest, DualFrequencyB1IB3IBroadcastModeUsesB1IAlone)
{
nav_t nav{};
eph_t eph{};
@@ -152,14 +152,13 @@ TEST(BdsTgdIonoTest, DualFrequencyB1IB3IRemovesTgd1BeforeCombining)
double iono_scale = -1.0;
const double corrected_pseudorange = prange(&obs, &nav, azel, 0, &options, &variance, &iono_scale);
// The DNAV clock is referenced to B3I, so the iono-free combination must
// use the TGD1-corrected B1I pseudorange:
// PC = (gamma13*(P1 - c*TGD1) - P2) / (gamma13 - 1)
const double gamma13 = std::pow(B3I_WAVELENGTH_M / B1I_WAVELENGTH_M, 2.0);
const double p1_corr = obs.P[0] - SPEED_OF_LIGHT_M_S * TGD1_S;
const double expected = (gamma13 * p1_corr - obs.P[2]) / (gamma13 - 1.0);
// With a broadcast iono model the bands are never combined: the B1I
// single-frequency user algorithm applies, PC = P1 - c*TGD1, and the
// modeled iono delay is applied unscaled. The combination is reserved to
// IONOOPT_IFLC (next test).
const double expected = obs.P[0] - SPEED_OF_LIGHT_M_S * TGD1_S;
EXPECT_NEAR(expected, corrected_pseudorange, 1.0e-6);
EXPECT_DOUBLE_EQ(0.0, iono_scale);
EXPECT_DOUBLE_EQ(1.0, iono_scale);
}
tests/unit-tests/signal-processing-blocks/pvt/galileo_e1_bgd_test.cc
+248
View File
@@ -0,0 +1,248 @@
/*!
* \file galileo_e1_bgd_test.cc
* \brief Tests for Galileo Broadcast Group Delay selection in PVT
* \author Carles Fernandez-Prades, 2026. cfernandez(at)cttc.es
*
* -----------------------------------------------------------------------------
*
* GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
* This file is part of GNSS-SDR.
*
* Copyright (C) 2010-2026 (see AUTHORS file for a list of contributors)
* SPDX-License-Identifier: GPL-3.0-or-later
*
* -----------------------------------------------------------------------------
*/
#include "rtklib_pntpos.h"
#include "rtklib_ppp.h"
#include <gtest/gtest.h>
#include <cmath>
namespace
{
constexpr double BGD_E1_E5A_S = 2.0e-8;
constexpr double BGD_E1_E5B_S = 4.0e-8;
constexpr double E1_WAVELENGTH_M = SPEED_OF_LIGHT_M_S / FREQ1;
constexpr double E5A_WAVELENGTH_M = SPEED_OF_LIGHT_M_S / FREQ5;
constexpr double E5B_WAVELENGTH_M = SPEED_OF_LIGHT_M_S / FREQ7;
void set_galileo_test_ephemeris(nav_t& nav, eph_t& eph, int sat)
{
eph.sat = sat;
eph.code = 1; // I/NAV E1B data source
eph.tgd[0] = BGD_E1_E5A_S;
eph.tgd[1] = BGD_E1_E5B_S;
nav.n = 1;
nav.eph = &eph;
}
} // namespace
TEST(GalileoE1BgdTest, SingleFrequencyPositioningUsesE1E5bBgd)
{
nav_t nav{};
eph_t eph{};
obsd_t obs{};
prcopt_t options{};
const int sat = satno(SYS_GAL, 1);
ASSERT_GT(sat, 0);
set_galileo_test_ephemeris(nav, eph, sat);
nav.lam[sat - 1][0] = E1_WAVELENGTH_M;
nav.lam[sat - 1][2] = E5A_WAVELENGTH_M;
obs.sat = sat;
obs.P[0] = 24000000.0;
obs.code[0] = CODE_L1C;
options.ionoopt = IONOOPT_BRDC;
const double azel[2] = {0.0, 1.0};
double variance = 0.0;
double iono_scale = -1.0;
const double corrected_pseudorange = prange(&obs, &nav, azel, 0, &options, &variance, &iono_scale);
EXPECT_NEAR(obs.P[0] - SPEED_OF_LIGHT_M_S * BGD_E1_E5B_S, corrected_pseudorange, 1.0e-6);
}
TEST(GalileoE1BgdTest, SingleFrequencyE1FromFnavUsesE1E5aBgd)
{
nav_t nav{};
eph_t eph{};
obsd_t obs{};
prcopt_t options{};
const int sat = satno(SYS_GAL, 1);
ASSERT_GT(sat, 0);
set_galileo_test_ephemeris(nav, eph, sat);
eph.code = 2; // F/NAV E5a data source
nav.lam[sat - 1][0] = E1_WAVELENGTH_M;
nav.lam[sat - 1][2] = E5A_WAVELENGTH_M;
obs.sat = sat;
obs.P[0] = 24000000.0;
obs.code[0] = CODE_L1C;
options.ionoopt = IONOOPT_BRDC;
const double azel[2] = {0.0, 1.0};
double variance = 0.0;
double iono_scale = -1.0;
const double corrected_pseudorange = prange(&obs, &nav, azel, 0, &options, &variance, &iono_scale);
EXPECT_NEAR(obs.P[0] - SPEED_OF_LIGHT_M_S * BGD_E1_E5A_S, corrected_pseudorange, 1.0e-6);
}
TEST(GalileoE1BgdTest, SingleFrequencyE5aKeepsE1E5aBgd)
{
nav_t nav{};
eph_t eph{};
obsd_t obs{};
prcopt_t options{};
const int sat = satno(SYS_GAL, 1);
ASSERT_GT(sat, 0);
set_galileo_test_ephemeris(nav, eph, sat);
nav.lam[sat - 1][0] = E1_WAVELENGTH_M;
nav.lam[sat - 1][2] = E5A_WAVELENGTH_M;
obs.sat = sat;
obs.P[2] = 24000000.0;
obs.code[2] = CODE_L5X;
options.ionoopt = IONOOPT_BRDC;
const double azel[2] = {0.0, 1.0};
double variance = 0.0;
double iono_scale = -1.0;
const double corrected_pseudorange = prange(&obs, &nav, azel, 0, &options, &variance, &iono_scale);
const double gamma = std::pow(E5A_WAVELENGTH_M / E1_WAVELENGTH_M, 2.0);
const double expected = obs.P[2] - gamma * SPEED_OF_LIGHT_M_S * BGD_E1_E5A_S;
EXPECT_NEAR(expected, corrected_pseudorange, 1.0e-6);
// E5a-only measurement: modeled E1 iono must be scaled by (f_E1/f_E5a)^2
EXPECT_DOUBLE_EQ(gamma, iono_scale);
}
TEST(GalileoE1BgdTest, SingleFrequencyE5bUsesE1E5bBgd)
{
nav_t nav{};
eph_t eph{};
obsd_t obs{};
prcopt_t options{};
const int sat = satno(SYS_GAL, 1);
ASSERT_GT(sat, 0);
set_galileo_test_ephemeris(nav, eph, sat);
nav.lam[sat - 1][0] = E1_WAVELENGTH_M;
nav.lam[sat - 1][2] = E5B_WAVELENGTH_M;
obs.sat = sat;
obs.P[2] = 24000000.0;
obs.code[2] = CODE_L7X;
options.ionoopt = IONOOPT_BRDC;
const double azel[2] = {0.0, 1.0};
double variance = 0.0;
double iono_scale = -1.0;
const double corrected_pseudorange = prange(&obs, &nav, azel, 0, &options, &variance, &iono_scale);
const double gamma = std::pow(E5B_WAVELENGTH_M / E1_WAVELENGTH_M, 2.0);
const double expected = obs.P[2] - gamma * SPEED_OF_LIGHT_M_S * BGD_E1_E5B_S;
EXPECT_NEAR(expected, corrected_pseudorange, 1.0e-6);
// E5b-only measurement: modeled E1 iono must be scaled by (f_E1/f_E5b)^2
EXPECT_DOUBLE_EQ(gamma, iono_scale);
}
TEST(GalileoE1BgdTest, DualFrequencyE1E5aBroadcastModeUsesE1WithItsBgd)
{
nav_t nav{};
eph_t eph{};
obsd_t obs{};
prcopt_t options{};
const int sat = satno(SYS_GAL, 1);
ASSERT_GT(sat, 0);
set_galileo_test_ephemeris(nav, eph, sat);
nav.lam[sat - 1][0] = E1_WAVELENGTH_M;
nav.lam[sat - 1][2] = E5A_WAVELENGTH_M;
obs.sat = sat;
obs.P[0] = 24000000.0;
obs.P[2] = 24000008.0;
obs.code[0] = CODE_L1C;
obs.code[2] = CODE_L5X;
options.ionoopt = IONOOPT_BRDC;
const double azel[2] = {0.0, 1.0};
double variance = 0.0;
double iono_scale = -1.0;
const double corrected_pseudorange = prange(&obs, &nav, azel, 0, &options, &variance, &iono_scale);
// With a broadcast iono model the E5a range present in the record is not
// combined with E1: the E1 range is used alone with the BGD matching the
// ephemeris clock model (I/NAV here), and the modeled iono delay applies
// unscaled. Combining per satellite would mix clock references within the
// constellation whenever a satellite lacks E5a.
const double expected = obs.P[0] - SPEED_OF_LIGHT_M_S * BGD_E1_E5B_S;
EXPECT_NEAR(expected, corrected_pseudorange, 1.0e-6);
EXPECT_DOUBLE_EQ(1.0, iono_scale);
}
TEST(GalileoE1BgdTest, PppUsesE1E5bBgd)
{
nav_t nav{};
eph_t eph{};
obsd_t obs{};
prcopt_t options{};
const int sat = satno(SYS_GAL, 1);
ASSERT_GT(sat, 0);
set_galileo_test_ephemeris(nav, eph, sat);
nav.lam[sat - 1][0] = E1_WAVELENGTH_M;
nav.lam[sat - 1][1] = E5B_WAVELENGTH_M;
obs.sat = sat;
obs.P[0] = 24000000.0;
obs.L[0] = 100.0;
obs.code[0] = CODE_L1C;
options.ionoopt = IONOOPT_OFF;
const double position[3] = {0.0, 0.0, 0.0};
const double azel[2] = {0.0, 1.0};
double measurements[2] = {};
double variances[2] = {};
int cycle_slip = 0;
ASSERT_EQ(1, corrmeas(&obs, &nav, position, azel, &options, nullptr, nullptr, 0.0,
measurements, variances, &cycle_slip));
EXPECT_NEAR(obs.P[0] - SPEED_OF_LIGHT_M_S * BGD_E1_E5B_S, measurements[1], 1.0e-6);
}
TEST(GalileoE1BgdTest, PppUsesE1E5aBgdFromFnav)
{
nav_t nav{};
eph_t eph{};
obsd_t obs{};
prcopt_t options{};
const int sat = satno(SYS_GAL, 1);
ASSERT_GT(sat, 0);
set_galileo_test_ephemeris(nav, eph, sat);
eph.code = 2; // F/NAV E5a data source
nav.lam[sat - 1][0] = E1_WAVELENGTH_M;
nav.lam[sat - 1][1] = E5A_WAVELENGTH_M;
obs.sat = sat;
obs.P[0] = 24000000.0;
obs.L[0] = 100.0;
obs.code[0] = CODE_L1C;
options.ionoopt = IONOOPT_OFF;
const double position[3] = {0.0, 0.0, 0.0};
const double azel[2] = {0.0, 1.0};
double measurements[2] = {};
double variances[2] = {};
int cycle_slip = 0;
ASSERT_EQ(1, corrmeas(&obs, &nav, position, azel, &options, nullptr, nullptr, 0.0,
measurements, variances, &cycle_slip));
EXPECT_NEAR(obs.P[0] - SPEED_OF_LIGHT_M_S * BGD_E1_E5A_S, measurements[1], 1.0e-6);
}
tests/unit-tests/signal-processing-blocks/pvt/gps_tgd_isc_test.cc
+5 -4
View File
@@ -166,7 +166,7 @@ TEST(GpsTgdIscTest, SingleFrequencyL5ScalesModeledIonoByGamma15)
}
TEST(GpsTgdIscTest, DualFrequencyL1L5UsesIscL1caOnGammaWeightedTerm)
TEST(GpsTgdIscTest, IonoFreeL1L5UsesIscL1caOnGammaWeightedTerm)
{
nav_t nav{};
eph_t eph{};
@@ -183,7 +183,9 @@ TEST(GpsTgdIscTest, DualFrequencyL1L5UsesIscL1caOnGammaWeightedTerm)
obs.P[2] = 24000008.0;
obs.code[0] = CODE_L1C;
obs.code[2] = CODE_L5X;
options.ionoopt = IONOOPT_BRDC;
// Only the iono-free model combines bands; any other iono model uses L1
// alone (see DualFrequencyL1L2AppliesMinusTgdToL1Pseudorange)
options.ionoopt = IONOOPT_IFLC;
const double azel[2] = {0.0, 1.0};
double variance = 0.0;
@@ -238,8 +240,7 @@ TEST(SppVarerrTest, IonoFreeCombinationAmplifiesVarianceByUpstreamFactorNine)
{
// the ionosphere-free noise amplification is (f1^2/(f1^2-f2^2))^2 +
// (f2^2/(f1^2-f2^2))^2, approximately 8.9 for GPS L1/L2, which upstream
// RTKLIB rounds to 3^2 = 9; rescode's dual-frequency fallback applies the
// same factor, so this pins the canonical value
// RTKLIB rounds to 3^2 = 9; this pins the canonical value
prcopt_t opt{};
opt.err[1] = 0.003;
opt.err[2] = 0.003;