1
0
mirror of https://github.com/gnss-sdr/gnss-sdr synced 2025-09-07 21:36:03 +00:00

refactor: consolidate per-satellite observations into a new sat_obs_t struct

This commit is contained in:
pedromiguelcp
2025-07-03 19:05:55 +01:00
committed by Carles Fernandez
parent bbe313b0f6
commit e767417c64
8 changed files with 99 additions and 103 deletions

View File

@@ -238,8 +238,8 @@ PVT.kf_system_ecef_pos_sd_m=0.001 ; static scenario
PVT.kf_system_ecef_pos_sd_m=0.001
;# VECTOR TRACKING CONFIG
PVT.enable_pvt_vtl=true ; enable/disable VTL
PVT.enable_pvt_output_vtl=false ; enable/disable receiver to use VTL solution
PVT.enable_pvt_closure_vtl=false ; enable/disable VTL feedback to tracking loops
PVT.enable_pvt_output_vtl=true ; enable/disable receiver to use VTL solution
PVT.enable_pvt_closure_vtl=true ; enable/disable VTL feedback to tracking loops
PVT.vtl_kinematic=false ; enable/disable kinematic model
PVT.vtl_init_pos_ecef_sd_m=10.0
PVT.vtl_init_vel_ecef_sd_ms=5.0

View File

@@ -1520,18 +1520,7 @@ bool Rtklib_Solver::get_PVT(const std::map<int, Gnss_Synchro> &gnss_observables_
}
}
int N_sv = valid_obs + glo_valid_obs;
std::vector<double> obs_pr_vec;
std::vector<double> tropo_m_vec;
std::vector<double> iono_m_vec;
std::vector<double> code_bias_m_vec;
std::vector<double> rho_vec;
double *rs;
double *dts;
rs = mat(6, N_sv);
dts = mat(2, N_sv);
result = rtkpos(&d_rtk, d_obs_data.data(), valid_obs + glo_valid_obs, &d_nav_data, obs_pr_vec, tropo_m_vec, iono_m_vec, code_bias_m_vec, rs, dts);
result = rtkpos(&d_rtk, d_obs_data.data(), valid_obs + glo_valid_obs, &d_nav_data);
if (result == 0)
{
@@ -1647,36 +1636,32 @@ bool Rtklib_Solver::get_PVT(const std::map<int, Gnss_Synchro> &gnss_observables_
int j = 0;
for (int n = 0; n < d_rtk.sol.ns; n++)
{
// one observation per satellite - exclude invalid observables
while (gnss_observables_iter->second.Pseudorange_m != obs_pr_vec[j])
// one observation per satellite (IFLC) - exclude invalid observables
while (gnss_observables_iter->second.Pseudorange_m != d_rtk.ssat[d_obs_data.at(j).sat - 1].sat_obs[0].prange)
{
if (obs_pr_vec[j] == 0)
++gnss_observables_iter;
if (!d_rtk.ssat[d_obs_data.at(j).sat - 1].vs)
{
j++;
}
else
{
++gnss_observables_iter;
}
}
const uint8_t idx6n = 6 * j;
const uint8_t idx2n = 2 * j;
// channel id
ch_id = gnss_observables_iter->second.Channel_ID;
vtl_data->rx_ch(ch_id) = ch_id;
// satellite info
vtl_data->sv_id(ch_id) = gnss_observables_iter->second.PRN;
vtl_data->sv_p.row(ch_id) = arma::Row<double>({rs[0 + idx6n], rs[1 + idx6n], rs[2 + idx6n]});
vtl_data->sv_v.row(ch_id) = arma::Row<double>({rs[3 + idx6n], rs[4 + idx6n], rs[5 + idx6n]});
vtl_data->sv_clk.row(ch_id) = arma::Row<double>({dts[0 + idx2n] * SPEED_OF_LIGHT_M_S, dts[1 + idx2n] * SPEED_OF_LIGHT_M_S});
double *svPosVel = d_rtk.ssat[d_obs_data.at(j).sat - 1].sat_obs[0].rs;
double *svClk = d_rtk.ssat[d_obs_data.at(j).sat - 1].sat_obs[0].dts;
vtl_data->sv_p.row(ch_id) = arma::Row<double>({svPosVel[0], svPosVel[1], svPosVel[2]});
vtl_data->sv_v.row(ch_id) = arma::Row<double>({svPosVel[3], svPosVel[4], svPosVel[5]});
vtl_data->sv_clk.row(ch_id) = arma::Row<double>({svClk[0] * SPEED_OF_LIGHT_M_S, svClk[1] * SPEED_OF_LIGHT_M_S});
vtl_data->sv_elev(ch_id) = d_rtk.ssat[d_obs_data.at(j).sat - 1].azel[1];
// atmosferic and code bias
vtl_data->tropo_bias(ch_id) = tropo_m_vec[j];
vtl_data->iono_bias(ch_id) = iono_m_vec[j];
vtl_data->code_bias(ch_id) = code_bias_m_vec[j];
vtl_data->tropo_bias(ch_id) = d_rtk.ssat[d_obs_data.at(j).sat - 1].sat_obs[0].dtrp;
vtl_data->iono_bias(ch_id) = d_rtk.ssat[d_obs_data.at(j).sat - 1].sat_obs[0].dion;
vtl_data->code_bias(ch_id) = d_rtk.ssat[d_obs_data.at(j).sat - 1].sat_obs[0].dcb;
// tracking info
vtl_data->CN0_dB_hz(ch_id) = gnss_observables_iter->second.CN0_dB_hz;

View File

@@ -1027,6 +1027,16 @@ typedef struct
} solopt_t;
typedef struct {
double rs[6]; /* satellite position|velocity (ecef) (m|m/s) */
double dts[2]; /* satellite clock bias/drift (s|s/s) */
double prange; /* pseudorange (m) */
double dtrp; /* tropospheric delay (m) */
double dion; /* ionospheric delay (m) */
double dcb; /* code bias delay (m) */
} sat_obs_t;
typedef struct
{ /* satellite status type */
unsigned char sys; /* navigation system */
@@ -1049,6 +1059,7 @@ typedef struct
double phw; /* phase windup (cycle) */
gtime_t pt[2][NFREQ]; /* previous carrier-phase time */
double ph[2][NFREQ]; /* previous carrier-phase observable (cycle) */
sat_obs_t sat_obs[NFREQ]; /* satellite observations */
} ssat_t;

View File

@@ -409,11 +409,10 @@ int rescode(int iter, const obsd_t *obs, int n, const double *rs,
const double *dts, const double *vare, const int *svh,
const nav_t *nav, const double *x, const prcopt_t *opt,
double *v, double *H, double *var, double *azel, int *vsat,
double *resp, int *ns, double *obs_pr, double *tropo_m, double *iono_m, double *code_bias_m)
double *resp, int *ns, double *pr, double *dion,
double *dtrp, double *dcb)
{
double r;
double dion;
double dtrp;
double vmeas;
double vion;
double vtrp;
@@ -443,7 +442,7 @@ int rescode(int iter, const obsd_t *obs, int n, const double *rs,
{
vsat[i] = 0;
azel[i * 2] = azel[1 + i * 2] = resp[i] = 0.0;
obs_pr[i] = 0;
pr[i] = 0;
if (!(sys = satsys(obs[i].sat, nullptr)))
{
@@ -486,7 +485,7 @@ int rescode(int iter, const obsd_t *obs, int n, const double *rs,
/* ionospheric corrections */
if (!ionocorr(obs[i].time, nav, obs[i].sat, pos, azel + i * 2,
iter > 0 ? opt->ionoopt : IONOOPT_BRDC, &dion, &vion))
iter > 0 ? opt->ionoopt : IONOOPT_BRDC, dion + i, &vion))
{
trace(4, "ionocorr error\n");
continue;
@@ -495,19 +494,17 @@ int rescode(int iter, const obsd_t *obs, int n, const double *rs,
/* GPS-L1 -> L1/B1 */
if ((lam_L1 = nav->lam[obs[i].sat - 1][0]) > 0.0)
{
dion *= std::pow(lam_L1 / LAM_CARR[0], 2.0);
dion[i] *= std::pow(lam_L1 / LAM_CARR[0], 2.0);
}
/* tropospheric corrections */
if (!tropcorr(obs[i].time, nav, pos, azel + i * 2,
iter > 0 ? opt->tropopt : TROPOPT_SAAS, &dtrp, &vtrp))
iter > 0 ? opt->tropopt : TROPOPT_SAAS, dtrp + i, &vtrp))
{
trace(4, "tropocorr error\n");
continue;
}
/* pseudorange residual */
v[nv] = P - (r + dtr - SPEED_OF_LIGHT_M_S * dts[i * 2] + dion + dtrp);
iono_m[i] = dion;
tropo_m[i] = dtrp;
v[nv] = P - (r + dtr - SPEED_OF_LIGHT_M_S * dts[i * 2] + dion[i] + dtrp[i]);
int j;
if (obs[i].code[0] != CODE_NONE)
{
@@ -521,8 +518,8 @@ int rescode(int iter, const obsd_t *obs, int n, const double *rs,
{
j = 2;
}
code_bias_m[i] = obs[i].P[j] - P;
obs_pr[i] = obs[i].P[j];
dcb[i] = obs[i].P[j] - P;
pr[i] = obs[i].P[j];
/* design matrix */
for (j = 0; j < NX; j++)
@@ -702,8 +699,8 @@ arma::vec rough_bancroft(const arma::mat &B_pass)
int estpos(const obsd_t *obs, int n, const double *rs, const double *dts,
const double *vare, const int *svh, const nav_t *nav,
const prcopt_t *opt, sol_t *sol, double *azel, int *vsat,
double *resp, char *msg, std::vector<double> &obs_pr_vec, std::vector<double> &tropo_m_vec,
std::vector<double> &iono_m_vec, std::vector<double> &code_bias_m_vec)
double *resp, char *msg, double *prange, double *dion, double *dtrp,
double *dcb)
{
double x[NX] = {0};
double dx[NX];
@@ -720,14 +717,6 @@ int estpos(const obsd_t *obs, int n, const double *rs, const double *dts,
int nv;
int ns;
char msg_aux[128];
double *obs_pr;
double *tropo_m;
double *iono_m;
double *code_bias_m;
obs_pr = mat(n + 4, 1);
tropo_m = mat(n + 4, 1);
iono_m = mat(n + 4, 1);
code_bias_m = mat(n + 4, 1);
trace(3, "estpos : n=%d\n", n);
@@ -771,7 +760,7 @@ int estpos(const obsd_t *obs, int n, const double *rs, const double *dts,
for (i = 0; i < MAXITR; i++)
{
/* pseudorange residuals */
nv = rescode(i, obs, n, rs, dts, vare, svh, nav, x, opt, v, H, var, azel, vsat, resp, &ns, obs_pr, tropo_m, iono_m, code_bias_m);
nv = rescode(i, obs, n, rs, dts, vare, svh, nav, x, opt, v, H, var, azel, vsat, resp, &ns, prange, dion, dtrp, dcb);
if (nv < NX)
{
@@ -826,13 +815,6 @@ int estpos(const obsd_t *obs, int n, const double *rs, const double *dts,
{
sol->stat = opt->sateph == EPHOPT_SBAS ? SOLQ_SBAS : SOLQ_SINGLE;
}
for (k = 0; k < n; k++)
{
obs_pr_vec.push_back(obs_pr[k]);
tropo_m_vec.push_back(tropo_m[k]);
iono_m_vec.push_back(iono_m[k]);
code_bias_m_vec.push_back(code_bias_m[k]);
}
free(v);
free(H);
free(var);
@@ -858,7 +840,9 @@ int estpos(const obsd_t *obs, int n, const double *rs, const double *dts,
int raim_fde(const obsd_t *obs, int n, const double *rs,
const double *dts, const double *vare, const int *svh,
const nav_t *nav, const prcopt_t *opt, sol_t *sol,
double *azel, int *vsat, double *resp, char *msg)
double *azel, int *vsat, double *resp, char *msg,
double *prange, double *dion, double *dtrp,
double *dcb)
{
obsd_t *obs_e;
sol_t sol_e = {{0, 0}, {}, {}, {}, '0', '0', '0', 0.0, 0.0, 0.0};
@@ -869,6 +853,10 @@ int raim_fde(const obsd_t *obs, int n, const double *rs,
double *vare_e;
double *azel_e;
double *resp_e;
double *prange_e;
double *dion_e;
double *dtrp_e;
double *dcb_e;
double rms_e;
double rms = 100.0;
int i;
@@ -893,6 +881,10 @@ int raim_fde(const obsd_t *obs, int n, const double *rs,
svh_e = imat(1, n);
vsat_e = imat(1, n);
resp_e = mat(1, n);
prange_e = mat(1, n);
dion_e = mat(1, n);
dtrp_e = mat(1, n);
dcb_e = mat(1, n);
for (i = 0; i < n; i++)
{
@@ -915,7 +907,7 @@ int raim_fde(const obsd_t *obs, int n, const double *rs,
std::vector<double> iono_m_vec;
std::vector<double> code_bias_m_vec;
if (!estpos(obs_e, n - 1, rs_e, dts_e, vare_e, svh_e, nav, opt, &sol_e, azel_e,
vsat_e, resp_e, msg_e, obs_pr_vec, tropo_m_vec, iono_m_vec, code_bias_m_vec))
vsat_e, resp_e, msg_e, prange_e, dion_e, dtrp_e, dcb_e))
{
trace(3, "raim_fde: exsat=%2d (%s)\n", obs[i].sat, msg);
continue;
@@ -953,7 +945,11 @@ int raim_fde(const obsd_t *obs, int n, const double *rs,
}
matcpy(azel + 2 * j, azel_e + 2 * k, 2, 1);
vsat[j] = vsat_e[k];
resp[j] = resp_e[k++];
resp[j] = resp_e[k];
prange[j] = prange_e[k];
dion[j] = dion_e[k];
dtrp[j] = dtrp_e[k];
dcb[j] = dcb_e[k++];
}
stat = 1;
*sol = sol_e;
@@ -1125,15 +1121,18 @@ void estvel(const obsd_t *obs, int n, const double *rs, const double *dts,
*-----------------------------------------------------------------------------*/
int pntpos(const obsd_t *obs, int n, const nav_t *nav,
const prcopt_t *opt, sol_t *sol, double *azel, ssat_t *ssat,
char *msg, std::vector<double> &obs_pr_vec, std::vector<double> &tropo_m_vec,
std::vector<double> &iono_m_vec, std::vector<double> &code_bias_m_vec, double *rs, double *dts)
char *msg)
{
prcopt_t opt_ = *opt;
// double *rs;
// double *dts;
double *rs;
double *dts;
double *var;
double *azel_;
double *resp;
double *prange;
double *dion;
double *dtrp;
double *dcb;
int i;
int stat;
std::vector<int> vsat(MAXOBS, 0);
@@ -1152,11 +1151,15 @@ int pntpos(const obsd_t *obs, int n, const nav_t *nav,
sol->time = obs[0].time;
msg[0] = '\0';
// rs = mat(6, n);
// dts = mat(2, n);
rs = mat(6, n);
dts = mat(2, n);
var = mat(1, n);
azel_ = zeros(2, n);
resp = mat(1, n);
prange = mat(1, n);
dion = mat(1, n);
dtrp = mat(1, n);
dcb = mat(1, n);
if (opt_.mode != PMODE_SINGLE)
{ /* for precise positioning */
@@ -1170,13 +1173,14 @@ int pntpos(const obsd_t *obs, int n, const nav_t *nav,
satposs(sol->time, obs, n, nav, opt_.sateph, rs, dts, var, svh.data());
/* estimate receiver position with pseudorange */
stat = estpos(obs, n, rs, dts, var, svh.data(), nav, &opt_, sol, azel_, vsat.data(), resp, msg, obs_pr_vec, tropo_m_vec,
iono_m_vec, code_bias_m_vec);
stat = estpos(obs, n, rs, dts, var, svh.data(), nav, &opt_, sol, azel_, vsat.data(), resp, msg, prange, dion,
dtrp, dcb);
/* raim fde */
if (!stat && n >= 6 && opt->posopt[4])
{
stat = raim_fde(obs, n, rs, dts, var, svh.data(), nav, &opt_, sol, azel_, vsat.data(), resp, msg);
stat = raim_fde(obs, n, rs, dts, var, svh.data(), nav, &opt_, sol, azel_, vsat.data(), resp, msg,
prange, dion, dtrp, dcb);
}
/* estimate receiver velocity with doppler */
if (stat)
@@ -1205,18 +1209,29 @@ int pntpos(const obsd_t *obs, int n, const nav_t *nav,
ssat[obs[i].sat - 1].azel[0] = azel_[i * 2];
ssat[obs[i].sat - 1].azel[1] = azel_[1 + i * 2];
ssat[obs[i].sat - 1].snr[0] = obs[i].SNR[0];
std::memcpy(ssat[obs[i].sat - 1].sat_obs[0].rs, &rs[i * 6], 6 * sizeof(double));
ssat[obs[i].sat - 1].sat_obs[0].prange = prange[i];
if (!vsat[i])
{
continue;
}
ssat[obs[i].sat - 1].vs = 1;
ssat[obs[i].sat - 1].resp[0] = resp[i];
ssat[obs[i].sat - 1].sat_obs[0].dts[0] = dts[i * 2 + 0];
ssat[obs[i].sat - 1].sat_obs[0].dts[1] = dts[i * 2 + 1];
ssat[obs[i].sat - 1].sat_obs[0].dion = dion[i];
ssat[obs[i].sat - 1].sat_obs[0].dtrp = dtrp[i];
ssat[obs[i].sat - 1].sat_obs[0].dcb = dcb[i];
}
}
// free(rs);
// free(dts);
free(rs);
free(dts);
free(var);
free(azel_);
free(resp);
free(prange);
free(dion);
free(dtrp);
free(dcb);
return stat;
}

View File

@@ -93,8 +93,8 @@ int rescode(int iter, const obsd_t *obs, int n, const double *rs,
const double *dts, const double *vare, const int *svh,
const nav_t *nav, const double *x, const prcopt_t *opt,
double *v, double *H, double *var, double *azel, int *vsat,
double *resp, int *ns, double *obs_pr, double *tropo_m, double *iono_m,
double *code_bias_m);
double *resp, int *ns, double *pr, double *dion,
double *dtrp, double *dcb);
/* validate solution ---------------------------------------------------------*/
int valsol(const double *azel, const int *vsat, int n,
@@ -105,15 +105,16 @@ int valsol(const double *azel, const int *vsat, int n,
int estpos(const obsd_t *obs, int n, const double *rs, const double *dts,
const double *vare, const int *svh, const nav_t *nav,
const prcopt_t *opt, sol_t *sol, double *azel, int *vsat,
double *resp, char *msg, std::vector<double> &obs_pr_vec,
std::vector<double> &tropo_m_vec, std::vector<double> &iono_m_vec,
std::vector<double> &code_bias_m_vec);
double *resp, char *msg, double *prange,double *dion, double *dtrp,
double *dcb);
/* raim fde (failure detection and exclution) -------------------------------*/
int raim_fde(const obsd_t *obs, int n, const double *rs,
const double *dts, const double *vare, const int *svh,
const nav_t *nav, const prcopt_t *opt, sol_t *sol,
double *azel, int *vsat, double *resp, char *msg);
double *azel, int *vsat, double *resp, char *msg,
double *prange, double *dion, double *dtrp,
double *dcb);
/* doppler residuals ---------------------------------------------------------*/
int resdop(const obsd_t *obs, int n, const double *rs, const double *dts,
@@ -144,8 +145,6 @@ void estvel(const obsd_t *obs, int n, const double *rs, const double *dts,
*/
int pntpos(const obsd_t *obs, int n, const nav_t *nav,
const prcopt_t *opt, sol_t *sol, double *azel, ssat_t *ssat,
char *msg, std::vector<double> &obs_pr_vec, std::vector<double> &tropo_m_vec,
std::vector<double> &iono_m_vec, std::vector<double> &code_bias_m_vec,
double *rs, double *dts);
char *msg);
#endif // GNSS_SDR_RTKLIB_PNTPOS_H

View File

@@ -2652,7 +2652,7 @@ void rtkinit(rtk_t *rtk, const prcopt_t *opt)
{
sol_t sol0 = {{0, 0}, {}, {}, {}, '0', '0', '0', 0.0, 0.0, 0.0};
ambc_t ambc0 = {{{0, 0}, {0, 0}, {0, 0}, {0, 0}}, {}, {}, {}, 0, {}};
ssat_t ssat0 = {0, 0, {0.0}, {0.0}, {0.0}, {'0'}, {'0'}, {'0'}, {'0'}, {'0'}, {}, {}, {}, {}, 0.0, 0.0, 0.0, 0.0, {{{0, 0}}, {{0, 0}}}, {{}, {}}};
ssat_t ssat0 = {0, 0, {0.0}, {0.0}, {0.0}, {'0'}, {'0'}, {'0'}, {'0'}, {'0'}, {}, {}, {}, {}, 0.0, 0.0, 0.0, 0.0, {{{0, 0}}, {{0, 0}}}, {{}, {}}, {}};
int i;
trace(3, "rtkinit :\n");
@@ -2762,8 +2762,7 @@ void rtkfree(rtk_t *rtk)
* notes : before calling function, base station position rtk->sol.rb[] should
* be properly set for relative mode except for moving-baseline
*-----------------------------------------------------------------------------*/
int rtkpos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav, std::vector<double> &obs_pr_vec, std::vector<double> &tropo_m_vec,
std::vector<double> &iono_m_vec, std::vector<double> &code_bias_m_vec, double *rs, double *dts)
int rtkpos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
{
prcopt_t *opt = &rtk->opt;
sol_t solb = {{0, 0}, {}, {}, {}, '0', '0', '0', 0.0, 0.0, 0.0};
@@ -2798,8 +2797,7 @@ int rtkpos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav, std::vector<d
time = rtk->sol.time; /* previous epoch */
/* rover position by single point positioning */
if (!pntpos(obs, nu, nav, &rtk->opt, &rtk->sol, nullptr, rtk->ssat, msg, obs_pr_vec, tropo_m_vec,
iono_m_vec, code_bias_m_vec, rs, dts))
if (!pntpos(obs, nu, nav, &rtk->opt, &rtk->sol, nullptr, rtk->ssat, msg))
{
errmsg(rtk, "point pos error (%s)\n", msg);
if (!rtk->opt.dynamics)
@@ -2841,8 +2839,7 @@ int rtkpos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav, std::vector<d
if (opt->mode == PMODE_MOVEB)
{ /* moving baseline */
/* estimate position/velocity of base station */
if (!pntpos(obs + nu, nr, nav, &rtk->opt, &solb, nullptr, nullptr, msg, obs_pr_vec, tropo_m_vec,
iono_m_vec, code_bias_m_vec, rs, dts))
if (!pntpos(obs + nu, nr, nav, &rtk->opt, &solb, nullptr, nullptr, msg))
{
errmsg(rtk, "base station position error (%s)\n", msg);
return 0;

View File

@@ -175,10 +175,7 @@ void rtkinit(rtk_t *rtk, const prcopt_t *opt);
void rtkfree(rtk_t *rtk);
int rtkpos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav,
std::vector<double> &obs_pr_vec, std::vector<double> &tropo_m_vec,
std::vector<double> &iono_m_vec, std::vector<double> &code_bias_m_vec,
double *rs, double *dts);
int rtkpos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav);
/** \} */

View File

@@ -602,15 +602,7 @@ void *rtksvrthread(void *arg)
}
/* rtk positioning */
rtksvrlock(svr);
std::vector<double> obs_pr_vec;
std::vector<double> tropo_m_vec;
std::vector<double> iono_m_vec;
std::vector<double> code_bias_m_vec;
double *rs;
double *dts;
rs = mat(6, 2);
dts = mat(2, 2);
rtkpos(&svr->rtk, obs.data, obs.n, &svr->nav, obs_pr_vec, tropo_m_vec, iono_m_vec, code_bias_m_vec, rs, dts);
rtkpos(&svr->rtk, obs.data, obs.n, &svr->nav);
rtksvrunlock(svr);
if (svr->rtk.sol.stat != SOLQ_NONE)