From d84069121a2850155eddba6a879da016b17876f3 Mon Sep 17 00:00:00 2001 From: Carles Fernandez Date: Wed, 9 Sep 2026 15:57:41 +0200 Subject: [PATCH] Let PVT.iono_model decide band combination in the single point solution --- docs/CHANGELOG.md | 10 + src/algorithms/libs/rtklib/rtklib_pntpos.cc | 72 ++--- tests/test_main.cc | 1 + .../pvt/bds_tgd_iono_test.cc | 15 +- .../pvt/galileo_e1_bgd_test.cc | 248 ++++++++++++++++++ .../pvt/gps_tgd_isc_test.cc | 9 +- 6 files changed, 298 insertions(+), 57 deletions(-) create mode 100644 tests/unit-tests/signal-processing-blocks/pvt/galileo_e1_bgd_test.cc diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 476a66341..e26bd11a0 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -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. diff --git a/src/algorithms/libs/rtklib/rtklib_pntpos.cc b/src/algorithms/libs/rtklib/rtklib_pntpos.cc index f47eafecb..31e86b0ba 100644 --- a/src/algorithms/libs/rtklib/rtklib_pntpos.cc +++ b/src/algorithms/libs/rtklib/rtklib_pntpos.cc @@ -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, diff --git a/tests/test_main.cc b/tests/test_main.cc index 0f634a1ff..9fb0b76b9 100644 --- a/tests/test_main.cc +++ b/tests/test_main.cc @@ -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" diff --git a/tests/unit-tests/signal-processing-blocks/pvt/bds_tgd_iono_test.cc b/tests/unit-tests/signal-processing-blocks/pvt/bds_tgd_iono_test.cc index 2946b6ae3..bad470bfa 100644 --- a/tests/unit-tests/signal-processing-blocks/pvt/bds_tgd_iono_test.cc +++ b/tests/unit-tests/signal-processing-blocks/pvt/bds_tgd_iono_test.cc @@ -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); } diff --git a/tests/unit-tests/signal-processing-blocks/pvt/galileo_e1_bgd_test.cc b/tests/unit-tests/signal-processing-blocks/pvt/galileo_e1_bgd_test.cc new file mode 100644 index 000000000..b934be5a5 --- /dev/null +++ b/tests/unit-tests/signal-processing-blocks/pvt/galileo_e1_bgd_test.cc @@ -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 +#include + +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); +} diff --git a/tests/unit-tests/signal-processing-blocks/pvt/gps_tgd_isc_test.cc b/tests/unit-tests/signal-processing-blocks/pvt/gps_tgd_isc_test.cc index d05024653..40c2239c8 100644 --- a/tests/unit-tests/signal-processing-blocks/pvt/gps_tgd_isc_test.cc +++ b/tests/unit-tests/signal-processing-blocks/pvt/gps_tgd_isc_test.cc @@ -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;