# Optimizing covariance estimation model for collaborative integrity monitoring in heterogeneous receiver satellite environments

Victor Vince, Alexandre Vervisch-Picois, Jose Rubio-Hernan, Dominique Heurguier

► To cite this version:

Victor Vince, Alexandre Vervisch-Picois, Jose Rubio-Hernan, Dominique Heurguier. Optimizing covariance estimation model for collaborative integrity monitoring in heterogeneous receiver satellite environments. Proceedings of the 37th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2024), Institut of navigation, Sep 2024, BALTIMORE, MD, United States. 10.33012/2024.19910. hal-04769640

HAL Id: hal-04769640

https://hal.science/hal-04769640v1

Submitted on 18 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimizing Covariance Estimation Model for Collaborative Integrity Monitoring in Heterogeneous Receiver Satellite Environments

Victor VINCE, *Thales SIX, SAMOVAR, Institut Polytechnique de Paris, Telecom SudParis*

Dominique HEURGUIER, *Thales SIX*

Alexandre VERVISCH-PICOIS, *SAMOVAR, Institut Polytechnique de Paris, Telecom SudParis*

Jose Manuel RUBIO HERNAN, *SAMOVAR, Institut Polytechnique de Paris, Telecom SudParis*

## BIOGRAPHY

**Victor Vince** : Victor Vince is a first-year PhD student at Institut Polytechnique de Paris, affiliated with Telecom SudParis and Thales SIX. His research focuses on developing innovative collaborative detection methods for situations where the GNSS system is compromised. Additionally, he is working on collaborative solutions for localization without using GNSS. His goal is to create robust and reliable approaches to ensure accurate navigation and localization even in the absence of GNSS signals. With a background in geophysics, he completed an internship at Thales, working on a localization system for his final year project. Contact: victor.vince@thalesgroup.com

**Dominique HEURGUIER** : Dominique Heurguier received his Engineering Degree from the "Ecole Nationale Superieure des Telecommunications" (Telecom Paris) in 1986. It also received, in 1987, a "Diplome d'Etudes Approfondies" in operational research. After been engineer for AERO company in robotic, he joined in 1998 Thales Communications and Security (Thales SIX) where he was involved in Spectrum Monitoring, then in Communication Electronic Warfare and now in Navigation warfare. He is also teacher at EUROSAE in charge of the GNSS alternative PNT solutions module « AED136 ».

**Alexandre VERVISCH PICOIS** : Alexandre Vervisch Picois received the diploma of engineer in Telecommunications in 2004 from the "Institut National des Télécommunications," France. He received the Ph.D. degree in electrical engineering, computer science and telecommunications from the University of Paris and "Institut National des Télécommunications" in 2010. He worked for 6 months in Canada on the design of a combined GPS and GALILEO software receiver. He is a member of the Navigation Group of the Department of Electronics and Physics at Telecom Sud-Paris Institut Télécom, France. He has been working since 2004 in the GPS field, on projects such as indoor location techniques and multipath mitigation methods. Since 2012, he has integrated the UMR Samovar. His current interests are related to GNSS repeater and pseudolite based indoor positioning and more specifically to mitigation of the interference limiting the accuracy and the availability of the system.

**Jose Manuel RUBIO HERNAN** : Jose Manuel RUBIO HERNAN received his engineering degree in telecommunications, specializing in electronics, from the Polytechnic University of Madrid (UPM), Spain, in 2012, and his Ph.D. in telecommunications from Paris VI, Sorbonne Universities, France, in 2017. Before his Ph.D., he worked on fiber optic deployment for a French company in France. Currently, he is an associate professor in the Department of Physics and Electronics at Télécom SudParis, IP Paris, France. His research interests include ICT safety and security in critical infrastructures and theoretical control approaches for detecting attacks against cyber-physical systems; GNSS receiver clock drift for spoofing attack detection; Indoor localization; And cognitive radio paradigm for spectrum management in the new generation of radio networks for connected objects, focused on physical and link layers.

## ABSTRACT

For several decades, GNSS systems have been vulnerable to spoofing attacks, which have become increasingly sophisticated and difficult to detect. This article builds on collaborative detection methods, particularly the Collaboration-Enhanced Receiver Integrity Monitoring (CERIM), which generalizes the concept of Receiver Autonomous Integrity Monitoring (RAIM) for multiple receivers. CERIM detects measurement faults using residuals from Position, Velocity, and Time estimations. This paper extends CERIM enabling the identification of satellite failures or spoofing attempts through anomalies in the global monitor
statistic derived from multiple receivers’ residuals.

Covariance matrices, essential for monitoring statistics, enhance the sensitivity of detecting system errors and malicious attacks. The literature on estimating these matrices is limited, often assuming unrealistic conditions such as receivers seeing the same number of satellites. This research addresses the challenge of estimating covariance matrices in scenarios where each receiver observes a variable number of satellites, leading to fluctuating dimensions over time. The study focuses on algorithms using both parity space and the complete satellite visibility space, known as Baseline and Common Residual algorithms.

Two approaches for estimating covariance matrices are explored: the empirical approach and the theoretical approach. The empirical approach adjusts the matrices based on real data over time windows, accommodating the variable number of visible satellites. The theoretical approach uses classical variance estimations. Tests examine the impact of receiver proximity, comparing cases where receivers are far apart (over 50 km) and closely located (1-2 km). This comparison helps define geographical clusters for collaborative detection. Results indicate that the empirical approach is currently the most effective for anomaly detection, but it requires adaptation based on the number of visible satellites. A reflection on various spoofing attacks, from simple to sophisticated, aims to assess CERIM's performance, which is difficult to achieve with traditional individual methods.

The innovation of this article lies in identifying the factors that enable accurate estimation of covariance matrices under various constraints, such as different masking (changing dimensions), as well as the influence of the distance between receivers, which allows for a good distribution of false alarm probabilities in the monitoring statistics. Additionally, the study also aims to determine if CERIM is capable of detecting spoofing on one or more receivers in the network, as the residuals from the positioning calculations increase significantly in the presence of spoofing and should therefore normally be reflected in CERIM.

# I. INTRODUCTION

In a world where GNSS (Global Navigation Satellite System) plays an increasingly significant role in critical infrastructures, it is essential to enhance its resilience. This article aims to highlight the major contribution of exchanges between receivers that enable network resilience. CERIM (Rife, 2011) (Rife, 2012) appears to be a promising candidate as it allows the construction of a global monitor statistic that follows the principle of RAIM (Brown, 1996) but on a collaborative scale. This approach enables the issuance of an alert if the integrity of the GNSS is impacted across multiple receivers by any anomaly with increased sensitivity. These anomalies can be natural (e.g., satellite error, significant solar eruption) or intentional (e.g., GNSS signal spoofing or jamming).

To construct an effective monitoring statistic, meaning with a very low false alarm probability rate and a false alarm rate conform to the specified probability, it is essential to accurately estimate the covariance matrix. This covariance matrix can be estimated in several ways. This paper will focus on empirical and theoretical estimation, as well as the hybridization of this covariance matrix estimated both empirically and theoretically.

The organization of this article is as follows: the first section reviews previous works on the operation of different CERIM algorithms and the estimation of covariance matrices. The second section will present the various covariance matrix estimations, which will be tested on real data with receivers in two configurations: close proximity (distance between receivers less than 1 km) and distant (distance between receivers greater than 20 km). Finally, the third section will discuss their potential performance in detecting sophisticated spoofing on one or more receivers in the collaborative network. The article concludes with a brief summary of the work and a conclusion highlighting the key points raised by these covariance estimations, particularly regarding their spoofing detection performance.

# II. PREVIOUS WORK

## 1. CERIM: The Extension of RAIM

CERIM, which stands for Collaboration-Enhanced Receiver Integrity Monitoring, is a collaborative method aimed at establishing a monitoring statistic to highlight position drifts related to natural technical errors. The power of CERIM lies in its simplicity, requiring no cumbersome infrastructure other than a communication link between each receiver and a centralization of information. CERIM provides robust integrity control, particularly effective when receivers are close to each other and operate in challenging environments with multiple reflections, obstructions, etc. (Rife, 2011), (Rife, 2012). Essentially,
CERIM extends RAIM (Brown, 1996), or **Receiver Autonomous Integrity Monitoring**, which has proven useful for detecting satellite failures but remains vulnerable to certain attacks like spoofing (Psiaki and Humphreys, 2016). This paper aims to understand the limitations of CERIM, even though the collaborative aspect shows promise against sophisticated spoofing attacks.

RAIM and CERIM use residuals to verify the accuracy of positioning calculations. These residuals reflect an error redundancy, enabling the establishment of a monitoring statistic. When calculating a position using pseudo-ranges provided by the receiver, the employed method seeks to determine the position that best fits the observed data. Due to measurement uncertainties, the solution to an over-determined system of equations inevitably shows residuals, as it represents a "best compromise" considering the least square criterion solution for the entire set of equations. If any of the measurements deviate significantly from this compromise solution, it becomes obvious in the residual and thus in the RAIM calculations involved. For standard 3 dimensional Position and Time GNSS calculation, when the receiver receives signals from at least five satellites, the problem is over-determined, allowing for the calculation of residuals. These residuals pertain to both the measured pseudo-ranges and the clock bias.

Thus, during the last iteration of the least squares estimation (see (Kaplan and Hegarty, 2017)), the residual vector $\Delta\rho$ ($\rho$ symbolizes the pseudorange) is not zero simply because it is impossible to satisfy all the observables at time t, particularly due to uncertainties.

$$ \Delta\rho^j = \rho_{es}^j - \rho_{obs} = \mathbf{r} \qquad (1) $$

The residuals in the satellite space are interdependent because the subtraction of theoretical measurements, based on the common solution estimation, introduces an algebraic coupling. This residual vector $\mathbf{r}$ is as long as the number of satellites $K_l$ visible to receiver $l$. The method used in this paper is a least squares estimation (cf. (Kaplan and Hegarty, 2017)). $\mathbf{H}$ corresponds to the observation matrix of dimension $\mathbb{R}^{K_l \times 4}$.

$$ \Delta\mathbf{x} = (\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\Delta\rho \qquad (2) $$

This least squares method for calculating residuals will be used in all cases addressed in this paper.

## 2. Parity Algorithm

In this section we present the theoretical approaches developed in (Rife, 2011), (Rife, 2012), and (Yang and Rife, 2016).

Before introducing the next algorithms, it is necessary to define a new space: the parity space (Patton and Chen, 1991). Indeed, the residual vector has so far been represented in the satellite space $\mathbf{r} \in \mathbb{R}^{K_l}$. Since the residuals are constrained to be orthogonal to the four columns of $\mathbf{H}$, the residuals in the vector $\mathbf{r}$ are not independent random variables. This redundancy leads to a degenerate, non-invertible covariance matrix of the residuals. Thus, this representation of the residual vector has the drawback of introducing an algebraic dependency among the residuals by mixing the measurements necessary for calculating the solution with the redundant measurements. It is therefore advisable to transform the residuals into an alternative form that eliminates the redundancy by considering a space, called the parity space, where the parity relations are independent of the unknown solution. This alternative form of the residual vector is generally called the parity vector $\mathbf{p}_l$. This is achieved by considering the parity space in which the random parity residuals are independent.

$$ \mathbf{p}_l = \mathbf{N}_l^T \mathbf{r}_l \qquad (3) $$

where $\mathbf{N}_l^T$ is the null space matrix of the observation matrix $\mathbf{H}_l$ such that: $\mathbf{N}_l^T \mathbf{H}_l = 0$. It is of dimension $\mathbb{R}^{K_l \times K_l - 4}$.

The following monitor statistic is obtained:

$$ m_{base} = \sum_{l=1}^{L} \mathbf{p}_l^T \mathbf{Q}_l^{-1} \mathbf{p}_l \qquad (4) $$

According to (Yang and Rife, 2016), the matrix $\mathbf{Q}_l$ can be approximated as:

$$ \mathbf{Q}_l = \mathbf{N}_l^T \mathbf{R}_l \mathbf{N}_l \qquad (5) $$

where $\mathbf{R}_l$ corresponds to the covariance matrix of the residual vector $\mathbf{r}$ expressed in the satellite space for receiver $l$:

$$ \mathbf{R}_l = \mathbb{E}[\mathbf{r}_l \mathbf{r}_l^T] \qquad (6) $$
Another statistic can be calculated with the concatenated parity vector $\mathbf{p}_c$:

$$ \mathbf{p}_c^T = [\mathbf{p}_1^T \ \mathbf{p}_2^T \ \cdots \ \mathbf{p}_L^T] \qquad (7) $$

The statistic is then:

$$ m_{Cbase} = \mathbf{p}_c^T \mathbf{Q}_c^{-1} \mathbf{p}_c \qquad (8) $$

with the matrix $\mathbf{Q}_c$ corresponding to the covariance of the vector $\mathbf{p}_c$. The matrix $\mathbf{Q}_c$ can be constructed in blocks, where each block is the expected value of a pair of measurement sets, m and n. According to (Yang and Rife, 2016):

$$ \mathbf{Q}_c[i, j] = \mathbb{E} [\mathbf{p}_i \mathbf{p}_j^T] \qquad (9) $$

Following the explanation of (Yang and Rife, 2016), if the errors between the residuals are independent, then $\mathbf{Q}_c$ can be constructed as a block diagonal matrix from the covariance $\mathbf{Q}_l$ for each set of measurements n, giving:

$$ \mathbf{Q}_c[i, j] = \begin{cases} \mathbf{0} & i \neq j \\ \mathbf{Q}_n & i = j \end{cases} \qquad (10) $$

These monitoring statistics must have an anomaly detection threshold. An alert is triggered if:

$$ m_{base} > S_{base} \rightarrow alert \qquad (11) $$

If we assume that **nominal errors are independently and identically normally distributed**, the resulting monitoring statistic $m_{base}$ follows a $\chi^2$ distribution. With a fixed false alarm probability equal to $\alpha_c$, we can calculate the detection threshold:

$$ S_{base} = P_{\chi^2}^{-1} (1 - \alpha_c, d.d.l_{base}) \qquad (12) $$

With the total degrees of freedom of the Baseline monitoring statistic being:

$$ d.d.l_{base} = \sum_{l=1}^{L} (K_l - 4) \qquad (13) $$

The Baseline algorithm presents shortcomings regarding the masking effects experienced by the receivers. Indeed, the signal-to-noise ratio for a given satellite can deteriorate if a satellite not represented in a particular receiver's measurement set is included in $m_{base}$, thus increasing noise without a signal gain, according to (Yang and Rife, 2016).

### Common Residuals Algorithm:

This algorithm is designed to process data from receivers, each potentially tracking different sets of satellites. The algorithm must express residuals in a new comprehensive visible set, encompassing all distinct satellites tracked by at least one receiver. This set is referred to as the *all-view*, containing $K_{av}$ satellites.

Here, it is expected that a component of the residual vector is common across multiple receivers. Hence, it is crucial to separate the common residual from the specific residual to distinguish correlated and uncorrelated residual components among the receivers. The common residual is a projection of the global residual vector onto the set of satellites tracked by receiver $l$. The specific component of the residual vector is one that is not supposed to be correlated among the receivers (e.g., components due to multipath and thermal noise).

The notation $c_l$ is introduced to refer to the common residual in the parity space for receiver $l$, and $s_l$ to refer to the specific residual in the parity space for the same receiver. The sum of the common and specific components gives the complete residual vector in the parity space $\mathbf{p}_l$:

$$ \mathbf{p}_l = \mathbf{c}_l + \mathbf{s}_l \qquad (14) $$

The approximation made by (Yang and Rife, 2016) to obtain specific residuals is rather straightforward. It is necessary to have a good approximation of the influence of common residuals on the global residual. Recall that common residuals are generated by
error sources affecting all receivers, such as ionospheric and tropospheric effects. Specific errors directly involve receiver-related errors and the environment in which it operates, where it may experience varying levels of multipath and signal-to-noise ratio (SNR).

The monitor statistic of the common residual algorithm depends on the estimated common vector in the all-view space ($\hat{\mathbf{c}}_{av}$):

$$m_{rc} = \hat{\mathbf{c}}_{av}^{T} \mathbf{Q}_{\hat{\mathbf{c}}}^{-1} \hat{\mathbf{c}}_{av} \qquad (15)$$

According to (Rife, 2012), this monitor statistic is structured to improve performance as the number of receivers increases.

For each receiver at time t, we compute $\mathbf{A}_l$, which facilitates the transition from the parity space to the all-view space:

$$\mathbf{A}_l = \mathbf{N}_l^T \mathbf{P}_l \mathbf{N}_{av} \qquad (16)$$

The matrix $\mathbf{N}_{av}$ is computed using the observation matrix $\mathbf{H}$, which contains all satellites visible to the collaborative set of receivers. $\mathbf{N}_{av}$ has dimensions $\mathbb{R}^{K_{av} \times K_{av}-4}$.

The projection matrix, denoted $\mathbf{P}_n$, is used to match the column dimensions of the kernel matrix of each receiver to that of the entire visible set. This projection matrix is a permutation of the identity matrix, with rows removed for each satellite not tracked by the respective receiver $n$. Thus, it has dimensions $\mathbb{R}^{K_l \times K_{av}}$. The common residual can be linked to the concatenated parity vector through the matrix $\mathbf{A}_c$:

$$\mathbf{p}_c = \mathbf{A}_c \mathbf{c} + \mathbf{s}_c \qquad (17)$$

The matrices $\mathbf{A}_c^T$ and vector $\mathbf{s}_c^T$ concatenate terms from all receivers:

$$\mathbf{A}_c^T = [\mathbf{A}_1^T \mathbf{A}_2^T \cdots \mathbf{A}_L^T] \qquad (18) \\ \mathbf{s}_c^T = [\mathbf{s}_1^T \mathbf{s}_2^T \cdots \mathbf{s}_L^T]$$

According to (Rife, 2012), the all-view common residual $\hat{\mathbf{c}}_{av}$ is calculated:

$$\hat{\mathbf{c}}_{av} = \mathbf{A}^{\dagger} \mathbf{p} \qquad (19)$$

where $\mathbf{A}$ is the pseudo-inverse matrix weighted by the covariance of specific all-view residuals:

$$\mathbf{A}^{\dagger} = (\mathbf{A}_{c}^{T} \mathbf{Q}_{s} \mathbf{A}_{c})^{-1} \mathbf{A}_{c}^{T} \mathbf{Q}_{s}^{-1} \qquad (20)$$

with the covariance matrix of specific all-view residuals $\mathbf{Q}_s$:

$$\mathbf{Q}_s[i, j] = \mathbb{E}[\mathbf{s}_i \mathbf{s}_j^T] \qquad (21)$$

To obtain $\mathbf{Q}_s$, it is necessary to separate the common residual from the specific residual. Here, if we revisit equation 17, we obtain the following equation linking covariance matrices:

$$\mathbf{Q}_c = \mathbf{A}_l \mathbf{Q}_{av} \mathbf{A}_l^T + \mathbf{Q}_s \qquad (22)$$

Here, $\mathbf{Q}_{av}$ denotes the error associated with the influence of the troposphere and ionosphere on position estimation. This error is also influenced by parameters that change over time. It depends on the proximity between receivers; if it is greater than several kilometers, effects will influence specific residuals; conversely, if receivers are close, effects will influence common residuals. Thus, theoretically, distance has a significant effect on the estimation of the covariance matrix of specific residuals.

$$\mathbf{Q}_s[i, j] = \mathbb{E}[\mathbf{s}_i \mathbf{s}_j^T] = \begin{cases} \mathbf{0} & i \neq j \\ \mathbf{Q}_{s,i} & i = j \end{cases} \qquad (23)$$

In (Rife, 2012) (Yang and Rife, 2016), the matrix $\mathbf{Q}_{s,i}$ is approximated as $\mathbf{Q}_{s,i} \approx \mathbf{Q}_n$.

Therefore, it is necessary to remove common residuals using methods like simple differencing (The simple differencing method in GNSS involves using the difference in observations from the same satellite for two distinct receivers) (Kaplan and Hegarty,
2017), for example, to retain only specific residuals that are independent and specific to each receiver. However, this method requires that receivers observe the same satellites and be sufficiently close to each other.

To obtain the monitor statistic $m_{cr}$, the covariance matrix $Q_{\hat{c}}$ of $\hat{c}$ must be obtained, which corresponds to the covariance of the estimate:

$$Q_{\hat{c}} = \mathbb{E}[\hat{c} \hat{c}^T] \qquad (24)$$

It can be expressed differently by introducing the error covariance matrix $Q_{\epsilon\hat{c}}$ associated with the estimation process. Thus, we have:

$$Q_{\hat{c}} = Q_{av} + Q_{\epsilon\hat{c}} \qquad (25)$$

Using the previous methodology, we then obtain the expression for $Q_{\epsilon\hat{c}}$:

$$Q_{\epsilon\hat{c}} = (A_c^T Q_s A_c)^{-1} \qquad (26)$$

If the common error is removed from the positioning calculation, the approximation $Q_{\hat{c}} \approx Q_{\epsilon\hat{c}}$ is possible, provided all common influences on the error have been removed.

### III. ESTIMATIONS OF COVARIANCE MATRICES

This section discusses various estimations of covariance matrices, which are modified in various monitor statistics to establish probability distributions of false alarms (pfa). These pfa distributions are compared to a $\chi^2$ distribution to determine the best estimation. Here, two scenarios of real data are defined to evaluate the performances and influence of different covariance matrix estimations on the construction of monitor statistics Baseline and Common Residual with stable and unstable dimensions. These monitor statistics are estimated to obtain a distribution of correct pfa. The figures showing the distributions of pfa do not have a standardized x-axis, as it corresponds to the point from which the probability of false alarm is equal to 0.

In the context of signal detection theory, the probability of a false alarm (also known as the false positive rate) refers to the likelihood that a test incorrectly identifies a signal when there is none. This occurs when the test indicates a positive result for the presence of a signal (or condition) when it is actually absent.

#### 1. Theoretical Models

To construct theoretical covariance matrices, one can draw inspiration from existing solutions in RAIM. RAIM also uses the residual vector to establish its monitor statistic. In (Kaplan and Hegarty, 2017), the covariance matrix used is constructed as follows:

$$W^{-1} = \begin{bmatrix} \sigma_1^2 & 0 & \dots & 0 \\ 0 & \sigma_2^2 & \dots & 0 \\ \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & \sigma_K^2 \end{bmatrix} \qquad (27)$$

To achieve a pfa distribution similar to a $\chi^2$, it is imperative to find the $\sigma$ that minimizes the difference between the theoretical and actual curves. This method is quite crude and often does not perfectly match the $\chi^2$. Another way to calculate the theoretical matrix is to integrate into $\sigma$ the parameters affecting measurement uncertainties for each satellite such as elevation, ionosphere, etc., according to (Kaplan and Hegarty, 2017):

$$\sigma_i^2 = \sigma_{i,URA}^2 + \sigma_{i,URE}^2 + \sigma_{i,tropo}^2 + \sigma_{i,mp}^2 + \sigma_{i,revr}^2 + \sigma_{i,rgb}^2 \qquad (28)$$

These uncertainty parameters on GPS and GNSS in general come from (McGraw et al., 2000), (Blanch et al., 2013), (Blanch et al., 2012), (Joerger and Pervan, 2016), and are summarized in Table 1. This construction method is especially used to establish RAIM monitor statistics and limit false alarms (Joerger and Pervan, 2016).

Similarly, this construction adapts to the elevation of visible satellites, as mentioned in McGraw et al. (2000):

$$\sigma_i = a_0 + a_1 e^{-\theta_i / \theta_c} \qquad (29)$$

Parameters $a_0$, $a_1$, and $\theta_c$ are used, where $\theta$ corresponds to the satellite elevation for which the precision uncertainty is estimated. These parameters can be estimated in various ways. In the article by (Yang and Rife, 2016), they are optimally adjusted to the
<table><caption>**Table 1:** Table summarizing uncertainties on GPS and Galileo constellations</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Value</th></tr><tr><th>GPS</th><th>Galileo</th></tr></thead><tbody><tr><td>SV clock,orbit error (Blanch et al., 2012), σ<sub>URE</sub></td><td>0.5 m</td><td>0.67 m</td></tr><tr><td>SV clock,orbit error (Blanch et al., 2012), σ<sub>URA</sub></td><td>0.75 m</td><td>0.96 m</td></tr><tr><td>Fault free, maximun error</td><td>1 m</td><td>1 m</td></tr><tr><td>Residual tropo<sup>a</sup>, error, σ<sub>tropo</sub></td><td colspan="2">0.12 <sup>1.001</sup>⁄√<sub>0.002001+(sin(θ))<sup>2</sup></sub></td></tr><tr><td>Multipath error<sup>a</sup>, σ<sub>mp</sub></td><td>0.13 + 0.53 e<sup>-θ/10</sup></td><td>refer to table [A-1] from (Blanch et al., 2012)</td></tr><tr><td>Receiver noise<sup>a</sup>, error σ<sub>revr</sub></td><td>0.15 + 0.43 e<sup>-θ/6.9</sup></td><td>refer to table [A-1] from (Blanch et al., 2012)</td></tr></tbody></table>

data at time t. It is also possible to estimate these parameters using random sampling (Monte Carlo method). The order of magnitude of the parameters according to (McGraw et al., 2000) are: $a_0 = 0.35$ m, $a_1 = 1.20$ m, and $\theta_c = 15$ deg with $5 \le \theta_c \le 90$ deg.

Following, we will analyze the two configurations announced above, proximity and distant:

### a) Case n°1 close receivers

The chosen scenario consists of simultaneous acquisition of three receivers over a total duration of 15 minutes. Two receivers are mobile and move from point n°1 to point n°2 at an average speed of 5 km/h towards the fixed receiver SIRT (cf. Figure 1). These data were acquired on the Saclay plateau in Palaiseau in May 2024.

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Map showing the trajectories of the 2 mobile receivers and the fixed station SIRT (sources: *Geoportail*)</figcaption></figure>

The two mobile receivers are EVK-M8T models from Ublox. The fixed receiver is a NETR9 model from TRIMBLE.

<footnote>a $\theta$ is the satellite elevation angle in degrees.</footnote>
<figure><img src="image_3.png" /><figcaption>(a) pfa distribution from a monitor statistic constructed using a covariance matrix equal to $\sigma^2\mathbf{I}$ where $\sigma$ varies for each satellite and follows 29, $a_0 = 1.10$, $a_1 = 2.20$, $\theta_c = 15$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) pfa distribution from a monitor statistic constructed using a covariance matrix equal to $\sigma^2\mathbf{I}$ where $\sigma$ varies for each satellite and follows 29, $a_0 = 0.5$, $a_1 = 1.0$, $\theta_c = 15$.</figcaption></figure>

**Figure 2:** Figure showing distributions of false alarm probabilities depending on the monitor threshold. These distributions are established for the scenario in Figure 1, using CERIM Baseline monitor statistic constructed with the same covariance estimation method. These distributions are compared to the theoretical $\chi^2$ distribution.

Analyzing the results of Figure 2, we can conclude that, the residuals are not centred and reduced all the time, which deviates from the theoretical case since the $\chi^2$ test is valid only under these conditions. Moreover, achieving the simplistic vision of $\chi^2$ is almost never possible with real data, as the degrees of freedom, i.e., the number of visible satellites, constantly changes. Therefore, it will be important to first examine the behavior of monitor statistics with stable dimensions.

### b) Case n°2 distant receivers

The chosen scenario involves simultaneous acquisition of 3 receivers over a total duration of 7h hours. The distances between the stations are summarized in the Table 2. These 3 receivers are NETR9 models from TRIMBLE (cf. Figure 3). These receivers are base stations for RTK, providing high-quality data, which may also explain why the results are accurate.
<figure><img src="image_5.png" /><figcaption>**Figure 3:** Map showing the 3 fixed stations: AGDE, AGDS, SETE (sources: IGN site)</figcaption></figure>

<table><caption>**Table 2:** Table summarizing distances between the 3 receivers</caption><thead><tr><td>Stations</td><td>distance (km)</td></tr></thead><tbody><tr><td>AGDE - AGDS</td><td>1.918</td></tr><tr><td>AGDE - SETE</td><td>21.967</td></tr><tr><td>AGDS - SETE</td><td>20.525</td></tr></tbody></table>

Figures 4 and 5 each group three graphs representing the probability of false alarms as a function of the fixed threshold. In the first case, the monitor statistic is weighted with $W_n = \sigma^2 I$, where $\sigma$ is chosen to closely match the corresponding $\chi^2$ curve (with the correct degrees of freedom). In the second case, the covariance matrix is constructed diagonally using $\sigma$ from 28, specific to each satellite visible to the receiver. In the third case, the covariance matrix diagonally incorporates 29, also adaptive for each receiver and satellite. **The fixed parameters are identical for all receivers.**

<figure><img src="image_6.png" /><figcaption>**Figure 4:** Figure grouping distributions of false alarm probabilities as a function of monitor threshold. These distributions are established using a CERIM Baseline monitor statistic, constructed with three different covariance estimations. These distributions are compared to the theoretical $\chi^2$ distribution.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Figure grouping distributions of false alarm probabilities as a function of monitor threshold. These distributions are established using a CERIM Residual Common monitor statistic, constructed with three different covariance estimations. These distributions are compared to the theoretical $\chi^2$ distribution.</figcaption></figure>

After study both scenarios, we verify that when the covariance matrices are estimated theoretically, the distributions of the *pfa* appear similar regardless of the algorithm and the distance between the receivers does not seem to have a significant impact, approximating the distribution of a $\chi^2$.

## 2. Empirical Models

Empirical models are built by memorizing previous iterations over a sliding window, which typically yields high-quality distributions of false alarm probabilities (*pfa*). The scenario used involves simultaneous acquisition of 3 receivers (cf Figure 3) over a total duration of 35 min. The distances between the stations are summarized in Table 2.
<figure><img src="image_8.png" /><figcaption>**Figure 6:** Figure grouping distributions of false alarm probabilities as a function of monitor threshold. These distributions are established using a CERIM Baseline monitor statistic, constructed with two covariance matrices estimated over a sliding window of 100 samples. These distributions are compared to the theoretical $\chi^2$ distribution.</figcaption></figure>

The use of the empirical model for estimating covariance matrices to establish the Baseline monitoring statistic allows us to achieve a very good distribution of *pfa* as shown in Figure 6. A better distribution is obtained for the covariance matrix estimated by blocks.

For the common residuals algorithm, the distributions obtained are similar to those of the Baseline algorithm and, in some cases, present a better distribution of *pfa*. However, it is necessary to remove the contributions of specific errors, which is relatively easy to do with a simple differencing method. This method requires that the receivers used to perform the differencing observe the same satellites, which is quite complex in an environment where receivers do not always see the same satellites or only have a few in common, making the simple differencing method nearly impossible. Therefore, we will use the Baseline method for covariance matrix estimations in unstable dimensions. To obtain the data shown in Figure 8, the scenario used involves the simultaneous use of 4 receivers (cf Figure 7) over a total duration of 7 hours. The distances between the stations are summarized in Table 3. These 4 receivers are NETR9 models from TRIMBLE.
<figure><img src="image_9.png" /><figcaption>**Figure 7:** Map showing the 4 fixed stations: AGDE, AGDS, NARB, SETE (sources: IGN site)</figcaption></figure>

<table><caption>**Table 3:** Table summarizing distances between the 4 receivers</caption><thead><tr><td>Stations</td><td>distance (km)</td></tr></thead><tbody><tr><td>AGDE - AGDS</td><td>1.918</td></tr><tr><td>AGDE - SETE</td><td>21.967</td></tr><tr><td>AGDE - NARB</td><td>41.561</td></tr><tr><td>AGDS - SETE</td><td>20.525</td></tr><tr><td>AGDS - NARB</td><td>42.678</td></tr><tr><td>SETE - NARB</td><td>62.979</td></tr></tbody></table>

It is easy to obtain a good estimation of the covariance matrix when the dimension is stable. However, the appearance of a new satellite complicates this task. Since we do not know the variations of the residuals associated with this satellite, we must theoretically estimate these variations. Here, we propose a hybrid estimation of the covariance matrix, combining theoretical and empirical estimates, to compensate for the peaks observed in our statistics when new satellites appear. Figure 8 shows a distribution of *pfa* for the Baseline algorithm, where the covariance matrix is estimated by blocks, combining empirical and theoretical estimates. On the Figure 8, the estimate are of good quality and show that this monitoring statistic is quite accurate.

<figure><img src="image_10.png" /><figcaption>**Figure 8:** Distribution of false alarm probabilities (*pfa*) over a period of 7 hours (25200 samples) from real data collected by four receivers (table 3), using a surveillance statistic constructed with an empirical block estimation of the covariance matrix. These distributions are derived using the CERIM Baseline, where the covariance matrix is estimated empirically over a sliding window of 100 samples and theoretically.</figcaption></figure>
## IV. SPOOFING PERFORMANCES ESTIMATION

Here, we will simulate a continuous spoofing attack on the receivers network presented in Table 3. The objective is to demonstrate the network's ability to detect an attack aiming to compromise the receivers positions. The spoofing scenario for this test involves introducing a simulated bias by modifying the pseudo-distances to gradually deviate the position calculated by the receiver. This corresponds to the attack of a coherent spoofer which would be synchronised on GNSS constellation. Four tests are conducted. The detection threshold is set to achieve a false alarm probability of $10^{-5}$.

<table><caption>**Table 4:** Summary of spoofing tests on a network of 4 receivers</caption><thead><tr><td>Number of spoofed receivers</td><td>Time of spoofing appearance on the monitoring statistic (seconds)</td></tr></thead><tbody><tr><td>1</td><td>1 sec</td></tr><tr><td>2</td><td>4 sec</td></tr><tr><td>3</td><td>18 sec</td></tr><tr><td>4</td><td>26 sec</td></tr></tbody></table>

The tests in Table 4 reveal that this type of spoofing without signal interruption is visible in this monitoring statistic. This test demonstrates the ability of the receiver network to detect spoofing attacks autonomously, without position monitoring. It seems that a successful spoofing of a receiver of this kind of network, would necessitate a coherent attack on each receiver of the network. The spoofer thus faces tougher task. This was not the case with a single receiver, where sophisticated spoofing with signal interruption was almost undetectable (Psiaki and Humphreys, 2016). Thus, in the case of sophisticated spoofing, CERIM, with a well-estimated covariance matrix as previously described, will potentially be able to detect if one or more receivers are being spoofed relative to the others. Additionally, the spatial separation of the receivers will enhance this detection, under the assumption that the influence of a spoofer is limited to a relatively confined area.

## V. SUMMARY

In a nutshell, this research enhances anomaly detection in GNSS data by proposing more precise and robust methods for estimating covariance matrices, while emphasizing the importance of individual receiver integrity monitoring in a collaborative environment, thereby ensuring system resilience. The study demonstrates some techniques for estimating covariance matrices, adapting them to fluctuating dimensions due to satellite masking, and increasing the sensitivity of surveillance statistics in challenging environments. Despite the study shown in the paper, further evaluation of simulated spoofing attacks is needed to confirm a high probability of detection and increased system resilience.

## VI. CONCLUSION AND FUTUR WORKS

Through this work, we have shown that there are numerous possibilities for estimating covariance matrices. The innovation in these results lies in the adaptation of covariance matrices to the fluctuating dimensions over time, due to the masking of certain satellites. In this article, we have adapted the covariance matrix and increased the sensitivity of the monitoring statistics in challenging environments. After these encouraging initial results on spoofing detection, we will study the impact of more sophisticated spoofing attacks. Using the approach proposed in this paper, we expect a good detection probability and increased system resilience against such attacks.

## REFERENCES

* Blanch, J., Walter, T., Enge, P., Lee, Y., Pervan, B., Rippl, M., and Spletter, A. (2012). Advanced raim user algorithm description: Integrity support message processing, fault detection, exclusion, and protection level calculation. In *Proceedings of the 25th International Technical Meeting of The Satellite Division of the Institute of Navigation (ION GNSS 2012)*, pages 2828–2849.

* Blanch, J., Walter, T., Enge, P., Wallner, S., Amarillo Fernandez, F., Dellago, R., Ioannides, R., Fernandez Hernandez, I.,
Belabbas, B., Spletter, A., et al. (2013). Critical elements for a multi-constellation advanced raim. *Navigation: Journal of The Institute of Navigation*, 60(1):53-69.

Brown, R. G. (1996). Receiver autonomous integrity monitoring. *Global Positioning System: Theory and applications.*, 2:143–165.

Joerger, M. and Pervan, B. (2016). Fault detection and exclusion using solution separation and chi-squared araim. *IEEE Transactions on Aerospace and electronic systems*, 52(2):726–742.

Kaplan, E. D. and Hegarty, C. (2017). *Understanding GPS/GNSS: principles and applications*. Artech house.

McGraw, G. A., Murphy, T., Brenner, M., Pullen, S., and Van Dierendonck, A. (2000). Development of the laas accuracy models. In *Proceedings of the 13th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GPS 2000)*, pages 1212–1223.

Patton, R. J. and Chen, J. (1991). A review of parity space approaches to fault diagnosis. *IFAC Proceedings Volumes*, 24(6):65-81.

Psiaki, M. L. and Humphreys, T. E. (2016). Gnss spoofing and detection. *Proceedings of the IEEE*, 104(6):1258-1270.

Rife, J. (2011). Collaboration-enhanced receiver integrity monitoring (cerim). In *2011 14th International IEEE Conference on Intelligent Transportation Systems (ITSC)*, pages 13–18. IEEE.

Rife, J. (2012). Collaboration-enhanced receiver integrity monitoring with common residual estimation. In *Proceedings of the 2012 IEEE/ION Position, Location and Navigation Symposium*, pages 1042–1053. IEEE.

Yang, L. and Rife, J. (2016). Estimating covariance models for collaborative integrity monitoring. In *Proceedings of the 29th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2016)*, pages 1103–1113.