GS-GVINS: A Tightly-integrated GNSS-Visual-Inertial
Navigation System Augmented by 3D Gaussian Splatting

To achieve a consistently accurate and reliable navigation solution in complex real-world scenarios, multi-sensor fusion is widely recognized as the most effective approach. The tightly coupled integration of GNSS, visual, and inertial data enables the joint optimization of estimated parameters by fusing raw measurements and leveraging all available error sources from each sensor.

In recent years, multiple GVINS frameworks have been developed and published. [29] represents the first attempt to integrate raw GNSS measurements of pseudorange and Doppler shift into an optimization-based visual-inertial SLAM system. It demonstrates superior performance compared to existing visual-inertial SLAM and GNSS Single-Point-Positioning (SPP), especially in GNSS-denied environments. GVINS [3] introduces a coarse-to-fine initialization technique that accurately establishes the real-time transformation between global measurements and local states, significantly reducing the GNSS state initialization time. Extensive evaluations validate its accuracy and robustness compared to other SOTA visual-inertial-odometry (VIO) systems. However, GVINS supports only the SPP algorithm for GNSS factors, limiting its ability to fully exploit the potential of GNSS positioning. To better utilize the advantages of INS, IC-GVINS [4] augments visual feature tracking and landmark triangulation by incorporating precise INS information. This integration significantly improves GVINS performance in high-dynamic conditions and complex environments. However, IC-GVINS only performs GNSS solution-level integration with raw inertial and visual measurements in a batch optimization framework. To fully harness the capabilities of GNSS, GICI-LIB [2] incorporates nearly all GNSS measurement error sources in its sensor fusion. Its RTK-RRR processing mode delivers superior pose estimation performance compared to other SOTA GNSS-Visual-Inertial navigation systems, particularly in challenging environments.

Despite these innovations, the primary focus of advancements has been on GNSS and INS. However, further improvement of GNSS or INS performance is increasingly challenging and presents a bottleneck. Meanwhile, visual sensors remain low-cost alternatives, but visual-based navigation still heavily relies on feature tracking and the minimization of reprojection errors. These areas require further attention to enhance overall navigation capabilities in GVINS. With the rapid evolution of graphical computing resources, more sophisticated visual representations and computer graphics techniques are being developed and can now be efficiently applied to address a variety of engineering challenges. These advancements open new opportunities to push the boundaries of visual navigation.

Traditional SLAM systems rely on points [30], surfels [31] or voxel grids [32] as scene representations, enabling direct and fast computations. However, these methods often struggle to achieve high-fidelity mapping due to their fixed spatial resolution and the lack of correlation among the 3D primitives. Neural-based SLAM [33][34][35] offers improved mapping quality but suffers from computationally expensive training processes, making it unsuitable for many real-time applications. In contrast, 3DGS has demonstrated superiority in efficient, real-time, and high-resolution image rendering. This is achieved by using expressive anisotropic 3D Gaussians as scene representations and leveraging tile-based rasterization techniques. As a result, 3DGS-based SLAM systems outperform traditional SLAM methods by enabling rapid photo-realistic rendering, and achieve superior performance and efficiency in both mapping and tracking. GS-SLAM [22] is the first to incorporate real-time 3DGS into the SLAM pipeline using RGB-D rendering. It integrates an adaptive expansion strategy to add or remove noisy 3D Gaussians, which enhances the mapping performance. Additionally, GS-SLAM introduces an effective coarse-to-fine method for selecting reliable 3D Gaussians, enhancing camera pose optimization. MonoGS [24] represents the first application of 3DGS in real-time monocular SLAM. It replaces the offline Structure-from-Motion (SfM) process used in the original 3DGS algorithm with direct optimization against 3D Gaussians, taking advantage of their wide convergence basin for robust camera tracking. MonoGS also introduces a geometric verification and regularization technique to address ambiguities in incremental 3D dense reconstruction. SplaTAM [25] proposes a single RGB-D camera SLAM framework that utilizes a silhouette mask to capture the presence of scene density. Both online tracking and mapping evaluations show that SplaTAM outperforms existing methods. Despite these advancements, 3DGS-based SLAM faces challenges, including error accumulation during pose tracking, scale ambiguity in monocular 3DGS SLAM, and degradation of pose tracking caused by inconsistent or erroneous 3D Gaussian maps. Our navigation system takes advantage of the tight integration of 3DGS with GNSS and INS to address these challenges.

In GS-GVINS, a factor graph representation is employed to describe the non-linear Least-Squares (LSQ) problem, where the estimated parameters and measurements are represented as nodes, and the residuals are represented as edges. Figure 2 illustrates the factor graph of the LSQ problem in GS-GVINS, where the Gaussian states are represented by yellow rectangles as they are optimized independently from the FGO. The solution of the LSQ problem can be defined by (1).

The state vector is defined as $\bm{\chi}=[\bm{\chi}_{r},\bm{\chi}_{I},\bm{\chi}_{c}]$, where the subscripts $r$, $I$, $c$ and $GS$ correspond to GNSS receiver, INS, camera, and 3DGS, respectively. Their full elements are shown in (2), (3) and (4). ${}^{B}t^{T}_{r}$ represents the GNSS antenna extrinsic parameters, ${}^{W_{G}}p^{T}$ and ${}^{W_{G}}v^{T}$ denote the position and velocity of the vehicle body in Earth-Centered Earth-Fixed (ECEF) frame, $df_{r}$ represents the frequencies for each satellite system, $N^{s}_{rr_{b},i}$ refers to single-difference carrier phase ambiguities. For the INS states, ${}^{W}p^{T}$, ${q^{W}_{B}}^{T}$ and ${}^{W}v^{T}$ represent the position, orientation and velocity of the vehicle body in East-North-Up (ENU) world frame. Its origin is defined as the first estimated GNSS position during system initialization. $b_{a}$ and $b_{g}$ are the biases of accelerometer and gyroscope. For the feature-based camera states, ${}^{B}t_{c}$ and $q^{B}_{C}$ represent the camera’s extrinsic parameters, and ${}^{W}p_{l}$ refers to the 3D positions of landmarks in ENU frame. $\bm{z}$ denotes the measurements and $h$ are the non-linear measurement models. $\bm{z}_{p}$ and $H_{p}$ represent the pseudo-measurements and linearized Jacobian of the prior information, which are computed during marginalization. Note that $\bm{z}_{c}$, $h_{c}$ and $\bm{\chi}_{c}$ only refer to the information related to the image-feature-based reprojection error factor. In GS-GVINS, 3D Gaussian parameters are optimized independently from the FGO, considering only the INS states in the 3DGS measurement model.

3DGS factor applies the principle of differentiable 3DGS to provide additional information for pose estimation. In our monocular case, the $\mathcal{L}_{1}$ photometric error can be modeled based on the current Gaussian map and the estimated camera pose, illustrated in (5).

where $I(\mathcal{G},\bm{T}_{CW})$ represents the rendered image from the 3D Gaussian map $\mathcal{G}$ based on $\bm{T}_{CW}$. And $\bar{I}$ is the ground truth image from observation. In the loss computation, optimized affine brightness parameters are utilized to account for varying exposure conditions. An RGB boundary mask is applied to focus the penalty on edge pixels while disregarding non-edge regions. Additionally, the opacity parameter derived from rendering is incorporated to penalize pixels with low-opacity values, ensuring a more robust loss evaluation. We optimize the estimated camera pose by minimizing the photometric error in (5).

In the forward process of 3DGS, the RGB-D point clouds extracted from the image are initialized as 3D anisotropic Gaussians $\mathcal{G}$. In our outdoor monocular setting, depth estimates are initialized with a large integer value, added by noise with high variance. Each Gaussian $\mathcal{G}^{i}$ is characterized by its optical properties: color $\bm{c}^{i}$ and opacity $\bm{\alpha}^{i}$, and its geometric properties: mean (3D position of center) $\bm{\mu}^{i}_{W}$ and covariance (ellipsoidal shape) $\bm{\Sigma}^{i}$, defined in world space.

The view-dependent radiance described by spherical harmonics (SHs) are omitted in this work for performance efficiency and simplicity.

After 3D Gaussian initialization, the rasterization iterates over each $\mathcal{G}^{i}$ and project it from 3D world space $\mathcal{G}(\bm{\mu}_{W},\bm{\Sigma}_{W})$ into 2D image space $\mathcal{G}^{\prime}(\bm{\mu}_{I},\bm{\Sigma}_{I})$ through a projective transformation, illustrated by (7).

where $\pi$ is the projection matrix and $\bm{T}_{CW}\in SE(3)$ represents the camera pose that transform world frame into camera frame. $\bm{J}$ is the Jacobian of the affine approximate of the projective transformation and $\bm{W}$ is the rotation part of $\bm{T}_{CW}$. Given the 3D position of the center of each $\mathcal{G}^{i}$, we formulate a list $\mathcal{N}$ that sorts all the overlapping Gaussians in 3D camera-view space with respect to the ray of each pixel, based on their depths. By traversing $\mathcal{N}$ from front to back, the final pixel color in the 2D image space is synthesized through the volumetric rendering along a ray, as shown in (8)

where $c_{n}$ is the optimized color of $\mathcal{G}^{n}\in\mathcal{N}$, $\alpha_{n}^{\prime}$ is the opacity value in projected Gaussian space, explained by (9), and $\alpha_{j}^{\prime}$ represents the opacity of each previous Gaussian in the list. The per-pixel contribution from the listed Gaussians are decayed in the order of sequence, by the transmittance based on all previous blended Gaussians. In (9), $x^{\prime}$ and $\mu_{n}^{\prime}$ are the coordinates of image pixel and the projected Gaussian ’s center. To facilitate the rasterization efficiency, tile-based parallel rendering is applied in CUDA programming architecture.

In the backward process of 3DGS, the calculated loss value is used to compute the gradients through backward-propagation.

For mapping, both optical parameters and geometric parameters of 3D Gaussians are optimized through iterative rendering against the training view and first-order Stochastic Gradient Descent (SGD) techniques. Since the Gaussians’ attributes do not have correlation with parameters of GNSS or INS, the mapping is conducted independently from the FGO. Keyframes are selected based on the co-visibility of Gaussian map and the relative translation of vehicle. If the intersection of constructed Gaussians between the current frame and the last keyframe drops below a threshold, or if the relative translation since the last keyframe surpass the median depth of the current Gaussian map, the current frame is registered as a new keyframe. A small window of keyframes is maintained to capture the scene in the surrounding environment. A keyframe will be removed from the window if its overlap with the latest keyframe drops below a threshold, along with its associated 3D Gaussians, which will also be removed from the map. Furthermore, we prune redundant Gaussians based on their visibility. In the outdoor monocular case, accurately estimating the 3D positions of many Gaussians remains challenging, and those associated with dynamic objects should be excluded from the map. These Gaussians usually violate multi-view consistency and will quickly vanish during optimization. To refine the 3D Gaussians map during optimization, we follow the same procedure as the original 3DGS work for adaptive density control. We further adopt the isotropic regularization techniques from [24] to penalize the Gaussians which are highly elongated along the viewing direction and cause artefacts during rendering, further improve the tracking accuracy.

During tracking, both camera pose and the affine brightness parameters are optimized. To reduce the computational overhead of automatic differentiation, the explicit analytical Jacobian of $\bm{SE}(3)$ camera pose with respect to 3D Gaussians, derived from [24], is integrated into our FGO pipeline. In (7), both $\bm{\mu}_{I}$ and $\bm{\Sigma}_{I}$ are differentiable with respect to $\bm{T}_{CW}$. The minimal Jacobians of 6-D world-camera transformation can be derived using Lie algebra as:

The derivatives of camera pose lay on the manifold for minimal parametrization, ensuring the dimension of the Jacobians match the exact degrees of freedom of $\bm{SE}(3)$. The partial derivatives with respect to $\bm{T}\in SE(3)$ can be described by (12)

where $\tau\in\mathit{se}(3)$ represents the 6-D twists on the tangent space of manifold, $\circ$ is a group composition, $\mathrm{Log}$ and $\mathrm{Exp}$ are the logarithmic and exponential mappings between Lie Group and Lie algebra. Therefore, the Jacobians of camera pose in (10) and (11) are derived as:

where ^× represents the skew symmetric matrix form, $\bm{\mathrm{W}}^{\times}_{:,i}$ indexes the $i$th column of the rotation matrix.

The 3DGS factor computes analytical Jacobians of camera pose using the $\mathcal{L}_{1}$ photometric loss by rendering the 3D Gaussian map at viewpoint and compare it with the ground truth image. The quality of rendering strongly depends on the accuracy of the optimized 3D Gaussian map. A Gaussian map optimized with a higher mapping loss introduces more rendering errors in camera tracking, resulting in less reliable pose estimation.

In the window-based mapping strategy, all keyframes within the window contribute to the map optimization at a given viewpoint. However, in complex large-scale and continuous driving scenarios, it is infeasible to optimize a 3D Gaussian from multiple viewing directions. The optimized 3D Gaussians may appear distorted when rendered from different viewing directions later. This is a typical issue that can degrade tracking performance due to the resulting erroneous rendering.

To account for the uncertainties from the 3DGS mapping optimization, we introduce a weighting scheme for 3DGS factor based on its mapping loss. The full weighting scheme can be described by (14).

where $\bm{\sigma}^{2}_{\mathrm{3DGS},i}$ is the variance of 3DGS factor at the $i_{th}$ frame, $\mathcal{W}$ is the current keyframe window, $k$ represents the $k_{th}$ keyframe, $E^{k}_{photometric}$ is the photometric loss computed using (5), $\lambda_{isotropic}$ is the pre-defined isotropic weighting factor, $s_{j}$ and $\tilde{s}_{j}$ are scaling parameters and its mean for the $j_{th}$ Gaussian in the map. This considers the contribution of error sources from all the keyframes in the current window and the isotropic regularization of over-stretched Gaussians.

When the vehicle approaches extreme dynamic states, such as near-static conditions or high-dynamic motion, camera tracking performance during non-keyframes can be impacted by noisy rendering. Erroneous renderings are more likely to occur due to positional errors in 3D Gaussians and the predominantly forward-oriented motion of the vehicle. These errors may cause 3D Gaussians to appear misplaced or excessively close to the viewpoint, leading to rendering artifacts such as blurriness, holes, or floating elements, as shown in Figure 8.

The optimization of the 3D Gaussian map relies on all keyframes within the current window. However, the probabilities of keyframe registration decrease significantly when the vehicle undergoes rapid deceleration or comes to a halt. As a result, the current Gaussian map is updated much less frequently, as no modifications occur within the keyframe window. Additionally, due to the inherent scale ambiguity in a monocular setting, the optimized map may appear distorted when viewed from different angles. This distortion becomes more pronounced during high-dynamic motion, where larger displacements between sequential viewpoints introduce relatively more significant rendering artifacts, particularly in forward-oriented motion. As a result, rendering from extreme dynamic viewpoints can lead to high photometric loss during tracking.

To address this issue, we propose a motion-aware 3D Gaussian pruning mechanism. This approach removes redundant 3D Gaussians from the map that could introduce noise in the rendering during extreme motion. The pruning process is guided by the computed relative translation between the consecutive frames, and the accumulated opacity $\alpha$ from the rays that capture the most recently inserted 3D Gaussians (from the latest keyframe). This mechanism is illustrated in (16).

where $\mathrm{t}_{i}$ is the relative translation of the $i_{th}$ frame with respect to its previous frame, and $\lambda_{min}$ and $\lambda_{max}$ are the predefined thresholds of relative translation for low-dynamic and high-dynamic case, respectively; $\mathcal{G}$ denotes the set of 3D Gaussians in the current map, while $\mathcal{G}_{redundant}$ represent the Gaussians to be pruned. $\mathcal{R}$ is the set of rays that intersect the newly inserted 3D Gaussians, with $k$ denoting an individual ray from this set. $R_{k}$ represents the set of Gaussians along the $k^{th}$ ray that are positioned in front of the newly inserted Gaussians. Finally, $\mathcal{G}_{j}$ refers to each individual 3D Gaussian along the $k^{th}$ ray.

If the displacement between consecutive viewpoints falls below or exceeds predefined thresholds, we infer that the vehicle is undergoing extreme dynamic motion. As described in (16), 3D Gaussians in the current map are pruned by removing those originating from older keyframes and preventing the integration of new 3D Gaussians from recent keyframes. This approach ensures that the latest keyframe, which generally provides a clearer view and a more accurate camera-to-3D Gaussian orientation, is prioritized in extreme motion scenarios.

After sorting the 3D Gaussians along the camera ray, we identify and remove those that cause the accumulated opacity $\alpha$ to reach 0.5 before blending with new Gaussians. This process reduces the density of closer 3D Gaussians in the vehicle’s heading direction, preventing unstable Gaussians from accumulating in regions with near-stationary motion or rapid motion variations.

For our self-collected data, we utilized a multi-sensor integrated platform equipped with hardware-based time synchronization (milliseconds level accuracy) to collect GNSS rover, inertial, and visual raw measurements. The platform consists of a u-blox F9P GNSS receiver, an ICM-20948 MEMS IMU, and two FLIR Blackfly-S RGB cameras. An STM32 Micro-Controller-Unit (MCU) manages system time synchronization using the GPS PPS signal, as well as data streaming and storage. A u-blox patch antenna is connected to the u-blox F9P receiver. The GNSS, inertial and visual data are collected in 1 Hz, 100 Hz and 10 Hz respectively. For groundtruth generation, the system is equipped with a NovAtel SPAN system comprising an OEM7 receiver and an Epson EG320N IMU. The time synchronization algorithm and the full details can be found in [37]. GNSS reference station data were collected from a Trimble R9 receiver located on the rooftop of the CCIT building at the University of Calgary. The baseline distance between the reference station and the rover consistently remained below 10 km. All data are stored on a laptop using ROS 2 [38].

To evaluate our proposed system, we replay the ROS bag data in a pseudo-real time configuration. The data are published via ROS topics and streamed into the GS-GVINS in the order of their message timestamps. The entire system operates within the environment of Advanced Computing Cluster (ARC) at the University of Calgary. During runtime, the system utilizes 4 out of 7 units of an NVIDIA A100-MIG GPU with 40 GB VRAM, 4 cores of an Intel(R) Xeon(R) Silver 4316 CPU running at 2.30 GHz, and 70 GB of RAM.

Figure 4 provides views of scenes from each validation dataset. For consistency, we categorize the environments based on urbanization levels into three types: open-sky, sub-urban, and urban. Two self-collected datasets are used to evaluate system performance in open-sky scenario; the first part of the UrbanNav-Deep dataset and the UrbanNav-Medium dataset are used to assess performance in sub-urban environments; while the second part of the UrbanNav-Deep dataset and the UrbanNav-Harsh dataset are used for evaluation in urban settings.

Both Open-sky A and B datasets are collected near the University of Calgary, where A is captured near the CCIT building and B is primarily collected on a highway. The Sub-urban A dataset covers an area near Hung Hom in Hongkong, featuring side buildings, small tunnels, and loops. The Sub-urban B focus on the area near Tsim Sha Tsui Substation in the central area of Hongkong, encompassing medium-height buildings and numerous dynamic objects. The Urban A and B datasets are collected in deep urban environments where several GNSS signal blockages occur, resulting in frequent NLOS and multipath signals. Urban A is collected in the central area of Hung Hom, characterized by narrow roads surrounded by densely packed building and heavy traffic. Urban B is collected in an ultra-dense urban canyon with narrow streets flanked by high-rising buildings, bridges, pedestrians and dense vehicle traffic.

Figure 5 illustrates the GNSS status in both sub-urban and urban environments. The available satellite number, Geometric Dilution of Precision (GDOP), and GNSS RTK fix status are visualized. To exclude GNSS outlier measurements, thresholds are applied: a minimum GNSS signal-to-noise ratio (SNR) of 30.0 $dB$ and a minimum satellite elevation of 7.0 degrees. When the GNSS solution is unavailable, the satellite number is plotted as zero. A fix status of 3 indicates a float ambiguity resolution for the GNSS RTK solution.

To validate the performance of GS-GVINS, we computed the APE of the navigation solutions from GS-GVINS against the ground truth and benchmark them against the tightly coupled solutions from GICI-LIB and IC-GVINS. All datasets are evaluated with successful GNSS/IMU and visual initializations. Table 1 summarizes the APE results for each environment type using the mentioned datasets. The values in bold indicate the lowest APE in each category. For Urban-B data tested with IC-GVINS, the solution is not applicable as the estimator diverges.

It is evident that GS-GVINS outperforms both GICI-LIB and IC-GVINS in the majority of the datasets. In GNSS-friendly Open-sky A and B environments, GS-GVINS and GICI demonstrate competitive performance. In these scenarios, GNSS RTK provides high-accuracy absolute positioning, which predominantly influences the solution estimation. Since the GNSS RTK factor implementations in GS-GVINS is adopted from GICI-LIB, the improvements in absolute pose estimation are limited. Additionally, as inertial and visual navigation have minimal influences in these scenarios, the following discussion will focus on the sub-urban and urban environments.

In both the sub-urban and urban datasets, GS-GVINS demonstrates significant improvements of 3D translation estimation over GICI and IC-GVINS. In Sub-urban A, GS-GVINS reduces 3D translation RMSE by 4.05 % and 71.81 %, compared to GICI and IC-GVINS, respectively. In Sub-urban B, GS-GVINS further reduces translation errors by 17.26 % (GICI) and 71.93 % (IC-GVINS), along with a 14.62 % reduction in 3D rotation RMSE. In the more challenging environment, Urban A, GS-GVINS improves translation accuracy compared to GICI and IC-GVINS by 1.37 % and 39.44 %, respectively. Despite GS-GVINS achieves a slightly better rotation estimation than GICI, the lowest 3D rotation RMSE (0.913 degrees) is achieved by IC-GVINS. In Urban B, GS-GVINS achieves a 13.71 % reduction in 3D translation errors compared to GICI, further demonstrating its robustness in translation estimation under challenging conditions.

Despite of the improvements of 3D translation estimation by GS-GVINS, its improvements in 3D rotation estimation remain limited. This limitation in rotation accuracy is partly due to the depth initialization strategy for the monocular setting. To compromise the complex and dynamic driving environments, we uniformly initialize the depth of each new Gaussian based on the approximated median depth of the scene with added noise, then refine it during rasterization. However, this coarse depth initialization, combined with the lack of multi-view optimization for 3D Gaussians, can introduce rotation scale errors, which limit the accuracy of 3D rotation estimation in GS-GVINS.

Figure 9 illustrates the 3D absolute translation error over GPS time for each dataset in the ENU directions. The red dashed vertical line represents the initialization time for GS-GVINS and GICI, as GS-GVINS adopts its initialization algorithm from GICI, resulting in identical completion times. The 3D Gaussian map is initialized in the next image frame after feature-based visual initialization. The blue vertical line marks the initialization time for IC-GVINS. Due to the different initialization strategies of GICI and IC-GVINS, and the dependency of IC-GVINS on the availability of GNSS RTK solutions, their initialization times vary. In Figure 9 (d), the absolute translation error of IC-GVINS is not included, as the estimator diverges.

In Figure 9 (a), a notable reduction in the vertical (up) translation error from GS-GVINS is observed indicated by the red dashed box. This improvement can also be observed near the right-side tail of the 3D trajectory in Figure 10. The corresponding driving scenario is depicted in the first, fourth and fifth rows of Figure 6, where the vehicle is going through an underground tunnel. Extreme camera exposure conditions are observed during this time period, where the tracked visual features can be lost or mismatched. However, 3DGS maintains high-fidelity rendering despite these conditions. Furthermore, the fourth row of Figure 6 highlights a featureless scene, which makes visual feature detection and tracking challenging. Despite this, 3DGS accurately renders the scene while preserving fine details. Figure 7 (a) displays the rendering from the opacity shader, where the opacity values extracted from the Gaussian rendering are visualized using a jet colormap (warmer colors indicate higher opacity values, and colder colors indicate lower opacity values). The comb-like patterns on the tunnel walls are precisely textured by the optimized opacity values. Although GNSS performance deteriorates during this period, as reflected in Figure 5, the 3DGS factor effectively bounds the vertical translation errors.

In Figure 9 (b), GS-GVINS significantly reduces translation errors in the ENU directions compared to GICI between 95610s and 95660s (black box). This period coincides with frequent fluctuations in GDOP values and a sharp decline in satellite numbers, indicating poor GNSS performance. The corresponding scenario is shown in the third row of Figure 6, where the vehicle is reaching to a stop. The motion-aware pruning mechanism aids in achieving accurate rendering in this case. Figure 7 (b) presents the opacity map, revealing a well-defined outline and structure of objects in the scene, which helps establish better zero-motion constraints from the 3DGS factor. Around 95680s, a sharp increase in GDOP value, combined with fewer than five available GNSS satellites, results in large error spikes in both east and north directions. However, GS-GVINS significantly reduces the east and vertical positioning error (red box) compared to both GICI and IC-GVINS, bringing the vertical positioning error close to zero. The second row of Figure 6 illustrates this environment, where dense tree canopies and advertising boards cast shadows over the left side of the street, making feature tracking difficult.

Figures 9 (c) and (d) display ENU error plots for more urbanized environments (Urban A and B). In Urban A, prolonged low satellite availability and extreme GDOP values (indicated by the red box) force the estimator to rely more on inertial and visual navigation. The last row of Figure 6 characterizes the street scene, featuring narrow roads, high-rise buildings, and parked cars, all of which are mapped by 3DGS. In this scenario, the 3DGS factor effectively constraints the growth of translation errors, as GNSS service is largely unavailable. As shown in Figure 9 (c), GS-GVINS improves both east and north translation errors to GICI and IC-GVINS as show in the red box. However, during the periods indicated by the black boxes, 3DGS factor seems to hinder translation error reduction. As the DOP values improve and the number of available satellites increases during this period, the influence of the 3DGS factor likely diminishes, weakening the impact of GNSS-based positioning during optimization. GS-GVINS achieves a significant reduction in north-direction error compared to GICI, as highlighted by the blue box, improving vehicle turning performance as shown in the third row of Figure 10 (near the end of the trajectory). Most of the cases discussed earlier also occur in Urban B, where GS-GVINS constrains the east and up error highlighted by the red box in Figure 9 (d). During these periods, with varying illumination caused by urban shadows and the presence of numerous pedestrians, the 3DGS factor effectively controls error drift in feature-based visual-inertial navigation.