## 2 Materials and Methods
### Adaptive GNSS-RTK/LIO Fusion
The LiDAR factor constructed from scan-matching provides relative motion constraints. The motion increment derived from scan-matching is taken as the known measurements functioning as a reference, while that computed via sta  
Figure 3: The factor graph of the adaptive GNSS/LiDAR/IMU fusion procedure. The states to be estimated are the global position and the relative rotation of the vehicle. The LiDAR factor and the IMU factor are relative constraints of the 6D state. The GNSS factor is adaptively added depending on the GNSS-RTK availability assessment to provide absolute constrain of the position part. considered as observed measurements. The factor is designed to minimize the difference between them as follows:  
\[\mathbf{r}_{\mathcal{L},k}\Big{(}{}^{n}\mathbf{T}_{l_{k-1}},{}^{n}\mathbf{T}_{l_{ k}},\mathbf{X}\Big{)}=\bigg{(}\Big{(}{}^{n}\mathbf{T}_{l_{k-1}}\Big{)}^{-1}.{}^{n} \mathbf{T}_{l_{k}}\bigg{)}^{-1}\cdot\Big{(}\mathbf{x}_{k-1}{}^{-1}.\mathbf{x} _{k}\Big{)}, \tag{8}\]  
in which \({}^{n}\mathbf{T}_{l_{k}}\) are solved by scan-matching. The IMU factor derived from pre-integration is formulated in a similar way. The motion increment provided by IMU pre-integration and that derived from the states to be estimated are considered known measurements and observed measurements respectively. The difference between them namely the IMU factor is defined as below:  
\[\mathbf{r}_{B,k}\Big{(}{}^{n}\mathbf{Z}_{l_{k}}^{b_{k-1}},\mathbf{X}\Big{)}= \Big{(}{}^{n}\mathbf{Z}_{l_{k}}^{b_{k-1}}\Big{)}^{-1}\cdot\Big{(}\mathbf{x}_{k- 1}{}^{-1}.\mathbf{x}_{k}\Big{)}, \tag{9}\]  
in which \({}^{n}\mathbf{Z}_{l_{k}}^{b_{k-1}}\in\text{SE}(3)\) is the rotation and translation increment pre-integrated by IMU raw data with the bias corrected by local FGO. It is converted from the IMU frame to the ENU frame. The adaptively added GNSS factor once a GNSS-RTK solution is judged to be reliable by the GNSS-RTK availability assessment module directly restricts the position part of the state as follows:  
\[\mathbf{r}_{\mathcal{G},j}\big{(}{}^{n}\mathbf{g}_{j},\mathbf{X}\big{)}={}^{n }\mathbf{g}_{j}-{}^{n}\mathbf{p}_{j},j\in[0,w), \tag{10}\]  
where \({}^{n}\mathbf{g}_{j}\) is the known measurement of the position provided by GNSS-RTK. The difference between the position to be estimated and \({}^{n}\mathbf{g}_{j}\) is expected to be small. When the commercial-level GNSS receiver is exploited to generate location measurements under urban canyons, only 2D residuals respectively at the \"earth\" and the \"north\" direction are fused in the global optimization make more sense given the large error at the \"up\" direction. In short, the full objective function of the global optimization could be expressed as the sum of all the residuals in a window:  
\[ \begin{split}\mathbf{f}(\mathbf{X})=\sum\limits_{k=1,\cdots,w-1} \rho\bigg{(}\left\|\mathbf{r}_{\mathcal{L},k}\Big{(}{}^{n}\mathbf{T}_{l_{k-1} }{}^{n}\mathbf{T}_{l_{k}},\mathbf{X}\big{)}\right\|_{\mathbf{C},k}^{2}\bigg{)} \\ +\sum\limits_{k=1,\cdots,w-1}\rho\bigg{(}\left\|\mathbf{r}_{ \mathcal{B},k}\Big{(}{}^{n}\mathbf{Z}_{l_{k}}^{b_{k-1}},\mathbf{X}\big{)} \right\|_{\mathbf{C}_{\mathcal{B},k}}^{2}\bigg{)}\\ +\sum\limits_{j\in[0,w)}\rho\bigg{(}\left\|\mathbf{r}_{ \mathcal{G},j}\big{(}{}^{n}\mathbf{g}_{j},\mathbf{X}\big{)}\right\|_{\mathbf{ C}_{\mathcal{G},j}}^{2}\bigg{)}.\end{split}  \tag{11}\]  
In which \(\left\|(\cdot)\right\|_{\mathbf{\Sigma}}^{2}=(\cdot)^{T}\mathbf{\Sigma}^{-1}(\cdot)\) is the Mahalanobis distance. \(\rho(\cdot)\) is the Cauchy Kernel [54] for outlier suppression expressed as  
\[\rho(s)=c^{2}\log\Bigl{(}1+\frac{s}{c^{2}}\Bigr{)}, \tag{12}\]  
The control parameter \(c\) is set to 1 in our implementation [51]. \(\mathbf{C}_{\mathcal{L},k}\) is the covariance matrix weighting the confidence of the factor and is supposed to be constant in our implementation.