3.2.4. The Construction of the Test Statistic

The targets detected by the single frame point cloud are matched with the targets detected by the global map, and the test statistic is constructed based on the mean position deviation of the matched targets, where the number of matched targets is *f*.

The real position of the *j*th matched target *X<sup>j</sup>* can be obtained from the 3D LiDAR map at the *k*th epoch without faults. The observed position *X<sup>j</sup> <sup>k</sup>* can be obtained through the vehicle position *P<sup>k</sup> car*, which is estimated by GNSS/INS integrated positioning, and the relative positioning *L<sup>k</sup> target*, which is estimated by single frame point cloud target detection. The real position *X<sup>j</sup>* and the observed position *X<sup>j</sup> <sup>k</sup>* are approximately equal, as shown below:

$$X^{\dot{j}} \approx X\_k^{\dot{j}} = P\_{\text{carr}}^k + L\_{\text{tarçet}}^k \tag{10}$$

The observed position *X<sup>j</sup> <sup>i</sup>* of the *j*th matched target at the *i*th epoch with GNSS faults can also be obtained, and a large error is present between the 3D LiDAR map position *X<sup>j</sup>* and the matched position *X<sup>j</sup> i* .

$$X^j \neq X\_i^j = P\_{car}^i + L\_{target}^i \tag{11}$$

The test statistic is then constructed by calculating the mean horizontal difference Δ*Xi* between the real target position *X<sup>j</sup>* and the observed target position *X<sup>j</sup> <sup>i</sup>* at each epoch for all *f* matched targets, i.e.,

$$
\Delta X\_i = \frac{1}{f} \sum\_{j=1}^{f} \left( X^j - X\_i^j \right)\_{horizontal} \tag{12}
$$

With the equations described in this section, the test statistic is constructed. In the next section, a threshold calculation approach is provided.

#### 3.2.5. The Threshold Constructed by an Adaptive Sliding Window

An adaptive sliding window approach is proposed to construct a dynamic threshold for fault detection. An adaptive weight sequence is set up to determine the weight of each value in the sliding window.

The threshold of the *i*th epoch is constructed by selecting the mean position deviation Δ*X* over the previous m epochs to construct a sliding window, which is defined as {Δ*Xi*−*<sup>m</sup>* ··· Δ*Xi*−2, Δ*Xi*−1}. The adaptive weight sequence is selected as {*αn*} = {*α*1, *<sup>α</sup>*2, ··· *<sup>α</sup>m*}, which satisfies the condition *<sup>m</sup>* ∑ *i*=1 *α<sup>n</sup>* = 1. Each value in the adaptive weight sequence {*αn*} is determined by a geometric sequence:

$$b^1 + b^2 + b^3 \cdots + b^m = \left(1 - b^{m+1}\right) / (1 - b) \tag{13}$$

As 0 < *b* < 1, Equation (13) can be transformed into

$$\left(b^1 + b^2 + b^3 \cdot \dots + b^m\right)(1 - b) / \left(1 - b^{m+1}\right) = 1\tag{14}$$

Assuming that (1 − *b*)/ <sup>1</sup> <sup>−</sup> *<sup>b</sup>m*+<sup>1</sup> = *δ*, Equation (14) is redefined as

$$\left(b^1 + b^2 + b^3 \cdot \dots + b^m\right)\delta = \delta b + \delta b^2 \cdot \dots + \delta b^m = 1\tag{15}$$

Thus, the weight sequence can be determined by

$$\{\mathfrak{a}\_{\mathsf{il}}\} = \{\mathfrak{a}\_1, \mathfrak{a}\_2 \cdot \cdots \cdot, \mathfrak{a}\_{\mathsf{m}}\} = \left\{\delta b, \delta b^2 \cdot \cdots \cdot, \delta b^{\mathsf{m}}\right\} \tag{16}$$

The threshold can be constructed as

$$\mathbf{T\_D} = \sum\_{n=1}^{m} a\_n \Delta X\_{i-n} \tag{17}$$

To avoid false alarms when the threshold is too small while considering the measurement errors of LiDAR and the lateral and longitudinal positioning accuracy standards stated in [44], we set a threshold of 0.5 m when there are no NLOS and multipath effects.

Therefore, the final threshold is given by

$$\text{TD}\_{\text{D}} = \max \left\{ 0.5, \sum\_{n=1}^{m} \alpha\_{n} \Delta X\_{i-n} \right\} \tag{18}$$

Finally, with the test statistic and threshold, the hypothesis test is constructed as

$$\begin{cases} \mathbb{H}\_0: no\ fault, & \Delta X\_i \le \mathcal{T}\_\mathcal{D} \\ \mathbb{H}\_1:fault, & \Delta X\_i > \mathcal{T}\_\mathcal{D} \end{cases} \tag{19}$$

With the hypothesis test given by Equation (19), GNSS positioning faults are detected in this study.

#### **4. LiDAR Aided Measurement Noise Estimation Adaptive Filter Algorithm**

*4.1. Existing Adaptive Filter Algorithms*

4.1.1. Single Fading Factor Adaptive Filter

Li et al. proposed an adaptive filter against outliers in [32] and designed a sliding window with length *n* to save the innovation of the filter γ*k*. If the GNSS fault is detected at the *k*th epoch, the innovation γ*<sup>k</sup>* is replaced by the average value of the last *n* innovation sequences and calculated as follows:

$$\mathbf{y}\_{k} = (\mathbf{y}\_{k-1} + \mathbf{y}\_{k-2} \cdot \dots + \mathbf{y}\_{k-n}) / n \tag{20}$$

γ*<sup>k</sup>* is brought into the KF calculation, and the modified innovation of the filter is saved in the sliding data window.

#### 4.1.2. Optimal Fading Factor Adaptive Filter

When the GNSS of an integrated navigation system has abnormal errors, the observation noise increases, which leads to an increase in the state covariance matrix **P**. Therefore, Geng et al. proposed an optimal fading adaptive filter algorithm, adding a fading factor *Sk* to the state prediction step of the state covariance matrix **P** [45]:

$$\mathbf{P}\_{k/k-1} = \mathbf{S}\_k \boldsymbol{\Phi}\_{k/k-1} \mathbf{P}\_{k-1} \boldsymbol{\Phi}\_{k/k-1}^T + \mathbf{Q}\_{k-1} \tag{21}$$

The convergence criterion of the KF is:

$$\mathbf{y}\_k^T \mathbf{y}\_k \le \kappa \text{tr}\left(\mathbb{E}\left(\mathbf{y}\_k \mathbf{y}\_k^T\right)\right) \tag{22}$$

where γ*<sup>k</sup>* is the KF innovation, *κ* is the regulating coefficient and *κ* ≥ 1. The strictest convergence condition is *κ* = 1, and γ*<sup>T</sup> <sup>k</sup>* γ*<sup>k</sup>* takes the minimum value satisfying Equation (22).

To realize the optimal fading factor, it must meet the strictest convergence condition:

$$\kappa = 1,\ \mathbf{y}\_k^T \mathbf{y}\_k = \text{tr}\left(\mathbb{E}\left(\mathbf{y}\_k \mathbf{y}\_k^T\right)\right) \tag{23}$$

This factor can be obtained from the properties of the innovation sequence:

$$\mathbb{E}\left(\mathbf{y}\_k \mathbf{y}\_k^T\right) = \mathbf{H}\_k \mathbf{P}\_{k/k-1} \mathbf{H}\_k^T + \mathbf{R}\_k \tag{24}$$

We substitute Equation (21) into Equations (23) and (24) to obtain:

$$\mathbf{y}\_k^T \mathbf{y}\_k = \mathbf{S}\_k \mathbf{H}\_k \boldsymbol{\Phi}\_{k/k-1} \mathbf{P}\_{k-1} \boldsymbol{\Phi}\_{k/k-1}^T \mathbf{H}\_k^T + \mathbf{H}\_k \mathbf{Q}\_{k-1} \mathbf{H}\_k^T + \mathbf{R}\_k \tag{25}$$

Equation (25) is traced and simplified, and the formula of *Sk* is obtained:

$$S\_k = \frac{\mathbf{y}\_k^T \mathbf{y}\_k - \text{tr}\left(\mathbf{H}\_k \mathbf{Q}\_{k-1} \mathbf{H}\_k^T + \mathbf{R}\_k\right)}{\text{tr}\left(\mathbf{H}\_k \boldsymbol{\Phi}\_{k/k-1} \mathbf{P}\_{k-1} \boldsymbol{\Phi}\_{k/k-1}^T \mathbf{H}\_k^T\right)}\tag{26}$$

It can be seen from Equation (26) that when the observed data are abnormal, the sum squares of the innovation γ*<sup>T</sup> <sup>k</sup>* γ*<sup>k</sup>* increase. Since the other parameters of the system remain unchanged in Equation (26), the current *Sk* is increased correspondingly, which is equivalent to increasing the state error covariance matrix **P** in Equation (21).
