Enhanced EKF-Based Time Calibration for GNSS/UWB Tight Integration

Tight integration of low-cost ultrawideband (UWB) ranging sensors with mass-market Global Navigation Satellite System (GNSS) receivers is gaining attention as a high-accuracy positioning strategy for consumer applications dealing with challenging environments. However, due to independent clocks embedded in commercial-off-the-shelf (COTS) chipsets, the time-scales associated with sensor measurements are misaligned, leading to inconsistent data fusion. Centralized, recursive filtering architectures can compensate for this offset and achieve accurate state estimation. In line with this, a GNSS/UWB tight integration scheme based on an extended Kalman filter (EKF) is developed that performs online time calibration of the sensors’ measurements by recursively modeling the GNSS/UWB time-offset as an additional unknown in the system state-space model. Furthermore, a double-update filtering model is proposed that embeds optimizations for the adaptive weighting of UWB measurements. Simulation results show that the double-update EKF algorithm can achieve a horizontal positioning accuracy gain of 41.60% over a plain EKF integration with uncalibrated time-offset and of 15.43% over the EKF with naive time-offset calibration. Moreover, a real-world experimental assessment demonstrates improved root-mean-square error (RMSE) performance of 57.58% and 31.03%, respectively.


I. INTRODUCTION
N OWADAYS, hundreds of applications in the massmarket segment are pushing the demand for continuous and dependable positioning, navigation and timing (PNT) services, and geo-positioning sensors such as the global navigation satellite system (GNSS) are key enablers to satisfy these growing needs [1], [2]. While different applications often pinpoint diverse requirements, enhanced accuracy and availability of the navigation solution are crucial for reliable positioning and guidance of agents. As a matter of fact, the robustness of the standalone GNSS technology is severely compromised in harsh environments, such as dense urban areas or wooded zones; the interplay between manifold impairing phenomena contributes to degrading the quality of GNSS standalone PNT solutions [3], [4], [5]. Hence, the design of GNSS-based navigation units in the massmarket segment has been increasingly addressing customizable, embedded architectures integrating low-cost ranging sensors [6], [7]. The latest proliferation in consumer electronics of ultrawide band (UWB) impulse radio (IR) transceivers-featuring small size and low power consumption-makes it an appealing candidate for hybridization with GNSS signals [8], [9]. As a carrier-free, spread-spectrum communication technology transmitting nonsinusoidal pulses with nanosecond life-cycle, UWB can ensure centimeter-level accurate ranging in dense, cluttered environments, thanks to supreme time resolution and remarkable obstacle penetration capabilities [10]. For outdoor applications, on-site deployment of georeferenced UWB beacons-acting as local environment landmarks-can help relative localization and dynamic tracking of mobile receivers in areas with reduced GNSS availability [11].
This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Traditional signal-processing methods for centralized sensor fusion, such as recursive maximum a posteriori (MAP) filters based on Bayesian estimation algorithms, assume accurately timestamped and synchronized measurements from different sensors to deliver trustworthy state-estimation performance [12], [13]. However, in the framework of low-cost GNSS/UWB tight integration, commercial UWB transceivers and mass-market GNSS chipsets operate as self-contained subsystems carrying independent, asynchronous clocks. Hence, GNSS and UWB measurements are timestamped with respect to different time-scales, and timing disturbances such as latency jitters or spurious clock drifts can affect the timestamping precision with respect to the true sampling instant. Moreover, at each estimation epoch, a lag exists between the measurements' timestamp and the time instant at which the integration filter (also known as hybridization, or navigation filter) processes the measurements. In particular, this lag is different for GNSS and UWB measurements because sampling rates are different and clocks are independent. Consequently, an unknown time-offset exists between the GNSS observables' timestamps and the UWB measurements' timestamps which models the relative misalignment between the GNSS and UWB time-scales as well as the shift of these scales with respect to the integration filter time-scale [14]. If neglected, this time-offset injects an inconsistency bias into the hybridization algorithm which may jeopardize state-estimation accuracy [15], [16]. This study aims at proposing a novel time calibration method to guarantee a consistent GNSS/UWB tight integration [12], thus enhancing the benefits introduced by UWB superior ranging accuracy over short distances.
In multisensor systems' literature, different methodologies have been explored to handle synchronization among selfcontained sensor units. On the one hand, hardware-level synchronization-especially popular in the framework of GNSS-aided strapdown inertial navigation system (INS) systems [17], [18]-typically exploits the pulse-per-second (PPS) signal from the GNSS receiver as triggering reference, and it cross-references the timing signals from the coupled sensors in order to establish a shared event base in the integration engine [19], [20]. As a preferred alternative when integrating commercial off-the-shelf (COTS) sensors which feature limited access to the hardware, software-based strategies are put in place; they operate a time calibration process involving accurate time-offset estimation in the integration architecture by leveraging sensors' measurement models [21], [22]. Consequently, software-based solutions can vary significantly depending on the sensors which are taken into account and the associated sensor integration models. In the existing literature, most of the proposed techniques address time calibration as a registration task, where timeoffsets are estimated through an offline, preprocessing step. For example, [23], [24], [25] suggested a continuous-time batch formulation of the time calibration problem, fitting with the framework of maximum likelihood estimation. Moreover, Kelly et al. [26] and Voges and Wagner [27] performed visual-inertial time calibration by temporally aligning orientation curves sensed by the independent sensors. Yet, joint multisensor optimization-based calibration strategies are explored in [28], [29]. As opposed to offline techniques, online temporal calibration via filtering-based methods (i.e., filtering-based calibration) is a promising alternative [30], [31]; this strategy models and recursively estimates the time-offset as an additional state-vector unknown under the hybridization filter state-space formulation [31], [32]. In this article, filtering-based time calibration is implemented in the context of GNSS/UWB tight integration and, at the time of writing, applications of online GNSS/UWB calibration are still missing in the literature.
After exploring the impact of uncalibrated time-offset on state-estimation performance via plain extended Kalman filter (EKF) hybridization, an EKF-based filtering model supporting time-offset calibration is established for GNSS/UWB tight integration. Nonetheless, the naive EKF-based time calibration model suffers pitfalls that can compromise integration performance under peculiar kinematic conditions [30]. In fact, the local identifiability [33] of the time-offset, that is, the accurate and unique estimation of such an unknown timing parameter based on the available observables-can be undermined by the relative geometry between the mobile receiver and the UWB nodes, thus leading to an ineffective time calibration. Therefore, a novel, double-update EKF architecture with an optimized weighting of UWB covariance statistics is put forward in order to mitigate the impact of the relative UWB anchors' geometry on accurate time-offset estimation. The proposed double-update EKF model is experimentally assessed to demonstrate the improved positioning accuracy against both plain EKF integration (i.e., without time-offset calibration) and EKF integration with naive time-offset calibration. To this end, after identifying critical kinematic conditions leading to performance deterioration induced by uncalibrated time-offset, experiments with simulated data are carried out over sample vehicular scenarios characterized by considerable variability of the observable processes. In addition, the effectiveness and practicability of the proposed time calibration methodology for GNSS/UWB tight integration are further validated with a real-world assessment of a vehicular trajectory in a suburban area.
The article outline is organized as follows: Section II lays the mathematical foundations for a discrete-time filtering model fitting with GNSS/UWB tight integration. After that, Section III introduces a convenient GNSS/UWB time-offset formulation and carries out a mathematical analysis to highlight the time-offset propagation on EKF state-estimation error. Then, Section IV, after enhancing the tight integration model for EKF-based time-offset calibration, presents the novel, double-update EKF with an adaptive weighting of UWB measurements. Eventually, Section V experimentally quantifies the navigation accuracy degradation caused by plain EKF tight integration and carries out a statistical performance assessment of the proposed calibration methodology both with simulated and real-world data.

II. EKF-BASED GNSS/UWB TIGHT INTEGRATION
Addressing the problem of recursive estimation of the timevarying state of a dynamic system, along with the continuous flow of noisy observation information, the hybridization filter builds upon a discrete-time system state space (DSS) formulation. The latter models the dynamic state evolution of the tracked mobile agent and the measurements, jointly with the associated noise statistics [34].

A. GNSS Code-Based Ranging Model
In the observation model for GNSS/UWB tight integration, GNSS ranging data include raw pseudorange and Doppler-shift observables. As a matter of fact, pseudorange measurements carry GNSS satellites-to-receiver ranging information corrupted by the contributions of receiver clock bias, atmospheric delay, and other impairments affecting the GNSS signal-inspace (SIS) propagation [35]. After compensation of modeled bias components, the pseudorange measurement equation for the i th GNSS satellite at generic epoch k can be written in metric units as [36] where r (i) k true, geometric range from satellite i to receiver at epoch k; r is the kth epoch position vector of satellite i expressed in Earth-centered Earth-fixed (ECEF) coordinates; r u,k = (r x,k , r y,k , r z,k ) is the kth epoch receiver position vector expressed in ECEF coordinates; δt u,k is the receiver clock bias term at epoch k; (i) G,k is the residual error in the pseudorange measurement from i th satellite.
Doppler-shift measurements, instead, bring information about the relative dynamics between each GNSS satellite and the receiver; hence, they are relevant to the estimation of both the receiver velocity and the receiver clock drift. By differentiating (1), the Doppler-shift model for the i th GNSS satellite at generic epoch k can be obtained in terms of satellite-to-receiver range variation rate [36] wherė k is the rate of change of the true, geometric range from satellite i to the receiver at epoch k; δ f u,k is the receiver clock drift; (i) G,k is the residual error in the pseudorange-rate measurement from the i th satellite.
As a matter of fact, the range-rateṙ (i) k can be thought of as the projection of satellite-to-receiver velocity vector on the transmitter-receiver line-of-sight (LOS) [36] where v is the kth epoch receiver velocity vector expressed in ECEF coordinates; I (i) k is the kth epoch LOS unit vector from the receiver position to the i th satellite position expressed in ECEF coordinates.
Moreover,ρ (i) G,k can also be expanded as [36] ρ (i)  (3) and (4) into (2), the Doppler measurement model can be expressed as a function of the unknown receiver velocity and clock-drift states [36]: B. UWB Measurement Model UWB leverages time-based protocols to pursue peer-topeer ranging. Adopting a time-of-arrival (TOA)-based ranging model [37], the baseline UWB ranging equation to the j th UWB anchor at generic epoch k can be written as where ρ ( j ) U,k is the measured two-way range from the receiver to the jth UWB anchor node; r U,k ) is the kth epoch position vector of UWB anchor j expressed in ECEF coordinates; ( j ) U,k is the residual error due to additive noise, non-LOS (NLOS) propagation, and further unmodeled effects [38].

C. GNSS/UWB State-Space Model
For the GNSS/UWB tight integration architecture under study, a total-state implementation is considered which estimates absolute properties of the system [39], [40]. Under the established framework of a discrete-time, EKF-based MAP tracking filter model [41], [42], [43], the state-vector at generic epoch k can be defined where a u,k = (a x,k , a y,k , a z,k ) is the receiver acceleration vector expressed in ECEF coordinates. As proven in Section III, GNSS/UWB time calibration is especially required in high-dynamic scenarios. Therefore, a constant acceleration model [44] is necessary to ensure enough accuracy of the system model to finely track the state-vector evolution.

1) GNSS/UWB System Model:
The discretization method of the state-transition function from a continuous-time linear time-invariant (LTI) system is given in [40]. Given a constant acceleration model for state dynamics [44], the state-transition matrix F k -obtained as the first-order truncation of the powerseries expansion of the linearized system matrix about the state-vector estimate-can be written as where t is the time between consecutive epochs (i.e., the time discretization step) and I 3×3 is the third-order identity matrix.
The system (or process) noise covariance matrix Q kwhich can also be obtained from the inference method when discretizing the state-transition function [40]-gathers any disturbance in the state characterization [42], and it can be split into two terms. By fixing a Cartesian direction (e.g., x-direction) in the ECEF reference frame, the process noise covariance for the corresponding navigation components (7) can be expressed as where S j x is the power spectral density (PSD) of acceleration noise in the ECEF-frame for the x-direction. Specifically, formulation (9) is valid under the assumption that process noise realizations come from a band-limited white noise process as long as the sampling rate is much less than the double-sided process noise bandwidth [40]. Equivalent modeling applies to the other Cartesian directions as well. Similarly, the process noise covariance for the timing parameters δt u,k δ f u,k T in (7) can be written as where S t and S f are the PSD of clock-bias and clock-drift noise, respectively.
2) GNSS/UWB Observation Model: Moving ahead in the DSS characterization, the measurement vector for GNSS/UWB tight integration embeds raw GNSS pseudorange and Doppler-shift observables (see Section II-A) as well as UWB-based auxiliary ranging information (see Section II-B); at epoch k, it can be formulated as where n is the number of tracked GNSS satellites and m is the number of auxiliary UWB anchors. Then, the linearized observation matrix-which models how measurements can affect the dynamic system state and which is determined from the known properties of the system [40]-can be written as In (12), the term H G n×3 identifies the Jacobian matrix resulting from first-order linearization of (1) and (5). It can be written as where r (i) is the true geometric range from the receiver to satellite i . In particular, the i th row of H G n×3 collects the Cartesian components of the unit steering vector pointing from the receiver position to the i th satellite position. Similarly, the term H U m×3 is the Jacobian matrix for UWB ranging model (6) and, in row j , it collects the Cartesian components of the unit steering vector from the receiver position to the j th UWB anchor position. Overall, H U m×3 has a similar structure to H G n×3 .

III. IMPACT OF TIME-OFFSET ON GNSS/UWB TIGHT INTEGRATION ACCURACY
In this section, after establishing a mathematical formulation for GNSS/UWB time-offset, a theoretical analysis is carried out to investigate its impact on EKF-based state estimation.

A. Time-Offset Definition
As mentioned in Section I, given independent rates and time-scales between the GNSS receiver clock and the UWB transceiver clock, an unknown time offset exists between the timestamps of raw GNSS observables and the timestamps of auxiliary UWB ranging measurements [31].
At a generic epoch k, the tight integration filter combines GNSS and UWB measurements in the observation model at a time instant t k tagged to the integration time-scale. Then, two quantities are identified as follows.
1) δt G,k expresses the unknown lag between the available set of GNSS measurements (timestamped in the GNSS time-scale) and t k . 2) δt U,k expresses the unknown lag between the available set of auxiliary UWB measurements (timestamped in the UWB time-scale) and t k . Then, the GNSS/UWB time-offset is defined as which expresses the misalignment of the timestamp associated with the set of auxiliary UWB measurements (tagged to the UWB time-scale) with respect to the GNSS time-scale. The time-offset t d,k includes the shift-under the GNSS timescale-between the timestamping times of GNSS and UWB measurements and it can be either a positive or a negative quantity at any epoch.
As a matter of fact, t d,k accounts for two time-varying effects as follows.
1) The relative misalignment between the GNSS and UWB measurements because of the different sampling rates and clocks of the independent sensors. 2) The lag of GNSS and UWB time-scales with respect to the time-scale of the centralized processing unit running the GNSS/UWB tight integration algorithm. While the former effect embeds nonidealities of individual sensor clocks (e.g., clock drift and latency jitters), the latter effect is contributed by manifold sources such as data-transfer latencies, hardware-level processing, and software overhead in the centralized processing unit. Furthermore, under the assumption of high UWB sampling rate and small drift of sensor clocks, t d,k can be modeled as a constant between consecutive estimation epochs. Hence, the GNSS/UWB timeoffset prediction at epoch k matches with the a posteriori timeoffset update at the previous epoch.
A graphical interpretation of the described framework is provided in Fig. 1. The three subplots are referred to as a common time-scale, that is, the integration time-scale. The top and middle subplots show the instants at which GNSS and UWB measurements are dumped, respectively. The bottom subplot shows which measurements the integration filter is processing. In case the measurements were provided at a high rate from at least one of the sensors, and the integration took place at a time t k at which the low-rate measurements are available, this time-offset would be negligible. However, in this work, the aim is to directly estimate t d,k as part of the system state, hence proposing a low-complexity strategy that relaxes constraints on measurement rates or related assumptions.

B. Mathematical Analysis
In light of the framework discussed in Section III-A, the GNSS/UWB measurement vector (11) collects observables that are not temporally consistent. In fact, at t k , the set of available GNSS ranging observables {ρ at a time instant which lags t k by δt U,k . Therefore, according to (14), {ρ Focusing on the auxiliary range to the j th UWB anchor, (6) can be reframed as which clearly highlights the mapping onto r U u,k . For ease of analysis, the residual error term, U,k , is neglected hereafter. Leveraging a continuous-time motion model for state dynamics [45], r G u,k can be related to r U u,k as follows: , v z (t)) expresses the instantaneous speed of the receiver and r,k identifies the displacement vector along the mobile user trajectory induced by t d,k . Given (16), the UWB range to the j th UWB anchor time-aligned to the can be expressed according to (17), as shown at the bottom of next page. Hence, a UWB ranging error term Substituting (17), and (15) at the numerator of (18), the following expression is obtained: which shows that, except for a positive scaling factor at the denominator, ρ,U,k for each of the m UWB anchors can be seen as an additive term to the set of auxiliary UWB range measurements. Consequently, it is possible to rewrite (11) as where y k collects time-aligned GNSS and UWB ranging observables. From (20), it is clear that t d,k introduces an additive error factor ρ,U,k in the GNSS/UWB observation model. This error then propagates on the a posteriori stateestimation error delivered by the tight integration filter.
The EKF innovation vector [41] can be written as wherex − u,k is the EKF a priori state estimate (i.e., the predicted state) at epoch k [41]. By expanding (21) according to (20), the EKF a posteriori state estimate (i.e., the updated state) at epoch k [41] is obtained aŝ with K k being the Kalman gain [41], [42] and ρ,U,k being the kth epoch UWB ranging error term. Eventually, the a posteriori state-estimation error due to t d,k can be obtained from (22) as wherex + u,k andx + u,k are the a posteriori state estimates obtained using y k and y k , respectively. Analyzing (23), both K k and ρ,U,k contribute to introduce errors in the integrated navigation solution. In particular, K k amplifies the UWB ranging error propagation on the a posteriori state estimation. In fact, when K k converges to an all-zeros matrix-a condition signaling that the EKF is trusting more the state prediction than sensors' observation information-the impact of ρ,U,k is largely mitigated. On the contrary, the state-estimation performance degradation due to t d,k exacerbates when the integration filter puts very high confidence on the observables' set. Moreover, according to (19), the components of ρ,U,k depend upon r,k which, in turn, is a function of v u (t) according to (16). Therefore, when the components of v u (t) take small values, the a posteriori state-estimation error contributed by t d,k reduces accordingly. Oppositely, ρ,x,k is expected to grow in high dynamics.
IV. ENHANCED EKF-BASED TIME-OFFSET CALIBRATION WITH DOUBLE-UPDATE FILTERING From the analysis presented in Section III-B, t d,k can affect the state-estimation performance as the tight integration filter processes observables which are not time-consistent. Hence, the baseline GNSS/UWB tight integration model of Section II is first improved in order to enable time calibration via EKF architecture. Then, a novel, double-update EKF filtering framework is put forward which enhances time calibration accuracy by adaptively accounting for the local identifiability of t d,k over consecutive epochs [33].

A. Improved GNSS/UWB Model for EKF-Based Time-Offset Calibration
An extended state vector is defined at epoch k where, compared to (7), the time-offset t d,k is introduced. Based on (24), and by reapplying linearization of the process function about the state-vector estimate and discretization [42], the state-transition matrix modifies as Using (25) in order to approximate the integral involved in the computation of r,k based on (16), (15) can be rewritten as which defines an improved UWB measurement model embedding t d,k . Examining (26), the receiver position information brought about by auxiliary UWB ranges (i.e., r U u,k ) gets compensated for the displacement locally induced by t d,k in order to geometrically match with the receiver position mapped by raw GNSS observables (i.e., r G u,k ). Moreover, the proposed modeling takes r G u,k as the unknown receiver position in the GNSS/UWB state-space formulation. Therefore, r G u,k matches with r u,k in (24) at every epoch. Furthermore, based on (26), modifications are required in the Jacobian matrix for UWB ranging H U m×3 , assuming m UWB anchors [8].

B. Double-Update Filtering With Adaptive Optimization
Comparing (26) against (15), the receiver position mapped by ρ ( j ) U,k is moved along the dynamic model using a specific, yet unknown, value of t d,k . Hence, the integration filter can exploit the DSS characterization to detect the unknown t d,k at t k based on the difference between ρ ( j ) For the j th UWB anchor, a function of the unknown state t d,k can be introduced which, by its definition, is proportional to the difference between ρ U,k . Given that, after a rough calibration of UWB transceiver clock [46], t d,k takes values of few tens of milliseconds, (27) can be expanded by substitution of (26) and (17) and a first-order approximation applies Apparently, f ( j ) (t d,k ) is a function that can drive the quality of t d,k estimation at t k . In fact, the sharper the envelope of f ( j ) (t d,k ) is, the more t d,k can be accurately and uniquely inferred from the set of available observables. Conversely, in case the first-order derivative of f ( j ) (t d,k ) approached a null value for some interval in the support, t d,k estimation would be jeopardized. By differentiating (28), it is obtained which highlights that the smaller is the inner product between the receiver velocity vector v u,k and the steering vector to the j th UWB anchor location h ( j ) U,k , the weaker is the local identifiability characterizing t d,k . This is equivalent to saying that the auxiliary UWB ranging observable associated with the j th UWB anchor brings little, if any, information to the integration filter to support accurate t d,k estimation. Fig. 2 displays the addressed geometrical framework. Moreover, when the inner product h ( j ) U,k · v u,k is small, t d,k would not even cause a significant mismatch between ρ ( j ) U,k and ρ ( j ) As a result, f ( j ) (t d,k ) would have little bearing on t d,k estimation. To cope with the aforementioned phenomenon, the covariance statistics of auxiliary UWB ranging measurements are adaptively weighted in the observation model for GNSS/UWB tight integration. Formally, accounting for the EKF measurement noise covariance matrix R k [41], the estimated variance of the auxiliary UWB range associated with the j th UWB anchor is amplified through the following epochdependent coefficient: where C k is an empirically set scaling parameter to rectify the instantaneous amplification effect of geometry on R k statistics. Focusing on the ranging contribution from the j th UWB anchor, (30) aims at weighting the degree of trust the integration filter should put on ρ ( j ) U,k to correct the a priori prediction of t d,k according to the identifiability conditions driven by the geometry-dependent behavior of (29). In particular, the smaller the sine of the angle between h k . In such a case, ρ ( j ) U,k is poorly informative to perform the update of t d,k , and its variance should be enhanced.
Nevertheless, according to Section IV-A, the adaptive weighting of UWB covariance statistics is expected to affect the a posteriori correction of all state variables in (24), not just t d,k . On this matter, it is essential to remark that, by its technological properties, UWB is able to deliver high-accuracy auxiliary ranging information [38], [47], [48] which can significantly contribute to the update of the navigation states in (24). It follows that the discussed weighting strategy must not affect the a posteriori estimation of these states given that (30) always operates an amplification of UWB ranging observables' variance. In light of the foregoing, a double-update EKF architecture is developed, which marginalizes the adaptive weighting on UWB covariance statistics to the a posteriori update of t d,k . For the remaining states in (24), instead, R k statistics are estimated by leveraging functional relationships with the available measurements [49]. The simplified block scheme in Fig. 3 illustrates the main stages of the proposed double-update model with embedded optimizations.

V. RESULTS ANALYSIS
In this section, experiments are presented to assess and validate both the theoretical analysis and the proposed time calibration methodology. In particular, Section V-B simulates the impact of uncalibrated t d under varying receiver kinematics when leveraging plain EKF tight integration. Then, Sections IV-B and V-C compare the simulated and real-world accuracy performance of the double-update EKF architecture against both the EKF model with naive t d calibration and the plain EKF hybridization. For ease of notation, the discretetime index k is omitted henceforth.

A. Experimental Setup and Methodology
To pursue the forthcoming assessment, sample vehicular scenarios are first generated through an radio frequency (RF) GNSS simulator-IFEN 1 network constellation simulator (NCS) Titan-in view of extracting the reference trajectory for the simulated vehicular target with a high position update rate of 250 Hz (i.e., 4 ms). The simulated trajectory-the Bernoullian Lemniscate-is centered in a location specified by its ellipsoidal World Geodetic System 1984 (WGS84) coordinates (latitude 45.063981 • , longitude 7.659017 • ) with a horizontal extension of 100 m. From the trajectory center, a network of three static UWB anchors is displaced at a fixed distance of 20 m and to a height of 5 m (measured along the vertical direction after taking local tangent plane approximation in correspondence with the trajectory center). A snapshot of the considered scenario is provided in Fig. 4(a). For the sake of rigorousness, it is remarked that the anchors' placement has been chosen with the sole purpose of establishing a uniform geometric distribution around the Lemniscate center. Nevertheless, despite being out of the scope of this study and left to future investigations, the geometry of the UWB network must be carefully taken into account since the relative localization estimates' accuracy at each time instant strongly depends on the local position of the mobile 1 Registered trademark. target with respect to the auxiliary UWB nodes [50]. Finally, the simulated GPS constellation with the motion track of the satellites is given in Fig. 4(b); an elevation mask of 15 • is used in order to prevent GNSS-based positioning performance deterioration from low-elevation satellites [36].
Based on the addressed scenario, simulations are run by configuring multiple receiver average speeds: 1, 2, 5, 10, 15, and 20 m/s. Over the different simulations, a shared time-span with a total duration of 3102 epochs is preserved to guarantee consistency in the comparison of the output results at different target speeds. Raw GNSS observables are logged via the NCS unit for GPS L1 C/A signals at 10 Hz rate (i.e., 100 ms) in the RINEX format. To prevent idealities, ionospheric effects are modeled inside the GNSS simulator. Logging of ephemeris data is also allowed by the NCS. On the UWB side, for each 250-Hz position fixed in the ground-truth, high-accuracy ranging data to the three UWB anchors are synthetically constructed by leveraging an empirical model proposed in [51], which has been preliminarily tested and validated in order to pursue the scopes of this research. Following this approach, UWB datasets with intentional t d have been constructed. The empirical parameter C k in (30) is set to 1 after some tests. Finally, by accounting for the GNSS/UWB tight integration framework treated in Section II, a C-language software implementing the EKF-based filtering architectures presented in this article is run in postprocessing on the retrieved datasets. For sake of completeness, Table I specifies relevant settings for  TABLE I  PARAMETER CONFIGURATION FOR THE PROCESS NOISE COVARIANCE  MATRIX Q k USED IN THE GNSS/UWB TIGHT INTEGRATION  POSTPROCESSING SOFTWARE matrix Q k used in the postprocessing software and common to all integration architectures. Moreover, the measurement covariance R k is constructed based on [49] and leveraging the synthetic UWB model [51]. Eventually, state-space domain initialization for the involved filters is obtained from a (WLSs) positioning solution [35].

B. Impact of Uncalibrated t d on Positioning Accuracy Under Varying Receiver Kinematics
The purpose of such analysis is to build a positioning error map which, depending on the mobile receiver kinematics and on the application-dependent accuracy performance requirements, can suggest whether it is worth or not implementing t d calibration for GNSS/UWB tight integration. To this end, given the sample vehicular scenarios at different average speeds (see Section V-A), synthetic UWB ranging datasets affected by diverse nominal t d 's are processed. Fig. 5 shows a heatmap chart of the positioning root-meansquare error (RMSE) for horizontal and vertical directions, respectively. Moreover, the line charts in Fig. 6 highlight both the horizontal and the vertical RMSE patterns over varying t d and are parameterized by different average receiver speeds. The positioning RMSE positively correlates with both receiver speed and t d , especially when accounting for the horizontal component (the upper plot of Fig. 6). More precisely, for t d values smaller than 10 ms, the measured RMSE is almost unaffected by the receiver kinematics. However, above 10 ms, it exhibits a nearly exponential increase with the average receiver speed. Looking at the overall horizontal RMSE pattern for small receiver speeds, the increase is nearly linear with t d , and this agrees with the mathematical analysis carried out in Section III-B. At a given estimation epoch, K k is mostly the same for different t d values. Besides, under low kinematics, the receiver speed can be considered constant over a time span equal to t d (which is of the order of tens of milliseconds). As a result, recalling (23), ρ,U,k and ρ,x,k increase proportionally with varying t d following a linear relationship. On the contrary, under higher kinematics, the receiver speed cannot be considered constant anymore due to the larger instantaneous acceleration. Therefore, the RMSE would not prompt a fixed slope anymore.
Concerning the vertical direction, instead, the RMSE pattern is apparently less regular. Overall, it grows for increasing t d , but it flattens off as low receiver kinematics are considered.  Furthermore, the vertical RMSE statistics are further penalized by the UWB anchors' geometry, and they fluctuate with the speed in low kinematic conditions as well.
Table II specifies the percentage proportion of horizontal RMSE contributed by uncalibrated t d to errors contributed by other sources. For simulations in low kinematics, t d marginally contributes to the measured RMSE. For instance, for a large t d of 100 ms, the RMSE increase amounts to only 12.44% at 1 m/s average speed. However, when the average receiver

C. Enhanced Double-Update EKF Error Statistics
In this section, with the aim of pursuing a state-estimation performance analysis in terms of both positioning error statistics and time-calibration accuracy, the simulated GNSS dataset for 20 m/s average receiver speed is postprocessed with a UWB ranging dataset affected by nominal t d = 40 ms. In fact, under these operating conditions, a remarkable positioning accuracy degradation has been measured for uncalibrated tight integration according to Table II.

1) Positioning Statistics and Time-Calibration Performance:
In Fig. 7, the time series of t d estimates are outlined for both EKF with naive t d calibration and double-update EKF. In addition, RMSE of t d estimates are summarized in Table III   for the considered time-calibration architectures. Apparently, the double-update EKF architecture delivers a smoother t d estimate (i.e., with lower variance) compared to the results obtained from the EKF with naive t d calibration. In particular, according to Table III for nominal t d = 40 ms, the doubleupdate EKF pursues a time-calibration accuracy-in RMSE terms-of 3.83 ms, which is higher than that achieved by the EKF with naive t d calibration (i.e., 4.90 ms) and which translates into an RMSE improvement of 21.77%. Fig. 8 highlights the time series of the positioning error in both the horizontal and the vertical directions. From a general standpoint, both filtering architectures implementing t d calibration substantially improve horizontal RMSE statistics compared to a plain EKF tight integration, although the vertical gain is moderate due to the poor UWB geometry. More in depth, accounting for the horizontal positioning accuracy in RMSE terms, the EKF with naive t d calibration reduces the error from 0.46 to 0.26 m compared to a plain EKF integration, and a further higher error reduction down to 0.18 m is obtained through the double-update EKF scheme.
Furthermore, Fig. 9 illustrates the empirical cumulative density function (ECDF) of the horizontal and the vertical positioning errors in East North Up (ENU) coordinates. What is more,    Eventually, horizontal summary statistics can be observed in Fig. 10 in terms of error mean, spread, and skewness characterizing the positioning estimates delivered by the analyzed filters. The superior performance of the proposed double-update EKF algorithm is hence assessed with respect to the other solutions.

3) Horizontal Error Statistics for Varying UWB Ranging
Accuracy: From Section IV, the underlying principle of t d calibration involves the identification of the misalignment between the position information mapped by GNSS and UWB ranging observables by leveraging the system dynamic model. Therefore, the higher the UWB ranging accuracy is, the more accurate the identification of positions' misalignment should be. In turn, t d estimation accuracy is expected to improve.
To test the impact of UWB ranging accuracy degradation on tight integration performance, synthetic UWB datasets with amplified UWB ranging measurements' noise standard deviation are generated for integer scaling factors equal to 1 (i.e., no amplification), 2, and 4. In the analysis, the simulated kinematic scenario at 20 m/s average receiver speed is considered under fixed t d = 40 ms.

D. Real-World Test
To validate the performance of the proposed filtering-based calibration methodology for GNSS/UWB tight integration, a real-world experiment has been carried out about a car ride in a suburban area of the Metropolitan city of Turin, Italy. The maximum achieved average vehicle speed is about 10 m/s. Fig. 11 provides a snapshot of the considered scenario. Although the chosen experimental environment is not expected to severely degrade the quality of GNSS observables, the results of the following analysis are still valuable. In fact, the primary goal of this section is not just validating GNSS/UWB tight integration per se, but rather to emphasize the accuracy improvements GNSS/UWB tight integration can benefit from embedding the proposed time calibration strategy. Furthermore, good GNSS conditions would even be a worst case for GNSS/UWB tight integration because auxiliary UWB measurements might not bring remarkable accuracy gains.
A u-blox ZED-F9P high-precision module has been leveraged which integrates multiband GNSS and real-time kinematic (RTK) technology [52]. Such a positioning module is commonly used in the industrial navigation and robotics markets. In addition, a high-gain, multiband SinoGNSS AT340 geodetic antenna has been deployed [53]. GNSS noisy pseudorange and Doppler-shift measurements have been logged through the u-blox module for the GPS constellation at 10 Hz rate. In parallel, the multiband RTK solution from the high-precision GNSS module has been retrieved in order to grant a ground truth (i.e., reference trajectory) useful for the estimation of the error statistics in the tightly integrated   Fig. 12.
On the UWB side, consumer EVB1000 boards from Qorvo Inc., have been employed [54]. In particular, three UWB modules have been installed on tripods as static anchors and their positions have been estimated at subdecimeter level accuracy. In addition, a fourth UWB module has been operated as a tag and installed on the road vehicle's roof. The deployment of both the UWB tag and the GNSS antenna on the vehicle are further highlighted in Fig. 11. As a remark, efforts have been made to minimize the lever arm between the phase centers of the GNSS and the UWB antennas. UWB measurements to the network of static anchors have been logged through the EVB1000 tag module at 5-Hz rate. When the UWB tag forward a new measurement to the laptop via the serial port, the Universal Time Coordinated (UTC) time from the laptop is recorded as the timestamp of the corresponding UWB measurement. Besides, a time-offset of 100 ms is intentionally added to each measurement sample in the UWB dataset in order to enhance the magnitude of GNSS/UWB time-offset. Nonetheless, being added to all the UWB measurement samples by the same amount, it does not undermine the methodology assessment. Furthermore, due to the inaccuracies and nonidealities of the laptop's clock as well as the transmission delays affecting UWB measurements, it has not been possible to retrieve the ground truth for the GNSS/UWB time-offset. Anyhow, the proposed timecalibration techniques can still be validated by assessing positioning error statistics with respect to the ground truth. Fig. 13 shows the positioning error time-series both for the horizontal and vertical components in the local ENU frame. As a matter of fact, vertical statistics are penalized  by the geometry of the deployed network of UWB anchors. In addition, the vertical component of the vehicle speed approaches zero for most of the trajectory course. Focusing on the horizontal component, Table VIII summarizes the positioning RMSE in units of meter. The double-update EKF achieves horizontal RMSE improvements of 31.03% and 57.58% compared to the EKF with naive t d calibration and the uncalibrated tight integration via plain EKF, respectively. Eventually, horizontal summary statistics for the real-world dataset are highlighted in Fig. 14 in terms of error mean, spread, and skewness characterizing the positioning estimates delivered by the analyzed filters. At the 75th percentile, the double-update EKF achieves horizontal gains of 38.54% over the EKF with naive t d calibration and of 63.67% over the plain EKF. As such, a globally exhaustive evidence of the superior performance of the proposed double-update EKF algorithm for online GNSS/UWB time-calibration is given. It is worth noticing that the conditions of the experimental environment (i.e., few UWB anchors while generous open sky visibility of GNSS satellites) do not limit the validity of the test. The scenario is suboptimal in terms of maximum achievable accuracy of a tight GNSS/UWB architecture, but it shows the benefits introduced by the proposed time calibration technique.

VI. CONCLUSION
Time calibration is of great concern in GNSS/UWB tight integration leveraging centralized EKF hybridization, and this article has explored the impact of uncalibrated GNSS/UWB time-offset on state-estimation accuracy both theoretically and experimentally. Based on these premises, an EKFbased framework has been proposed to address time-offset calibration. First, an improved GNSS/UWB tight integration model has been presented to enable the modeling of the unknown time-offset as part of the EKF-based state-space formulation. Then, after pointing out criticalities in timeoffset estimation owing to local losses of identifiability, an enhanced, double-update EKF architecture has been put forward which adaptively weights UWB covariance statistics. Simulation results demonstrate that, under challenging kinematic conditions causing remarkable accuracy deterioration for uncalibrated GNSS/UWB tight integration, the EKF with naive time-offset calibration can reduce the horizontal positioning RMSE from 0.4625 to 0.2557 m, while the double-update EKF with statistical modeling optimizations can achieve higher RMSE reduction down to 0.1931m. Moreover, results obtained with a real-world dataset further assess the superior doubleupdate EKF performance by highlighting horizontal gains of 57.58% compared to the uncalibrated tight integration via plain EKF and of 31.03% over the EKF with naive time-offset calibration.