Satellite Autonomous Clock Fault Monitoring with InterSatellite Ranges Using Euclidean Distance Matrices

Iiyama Keidai; Neamati Daniel,; Gao Grace

doi:10.33012/navi.764

Satellite Autonomous Clock Fault Monitoring with InterSatellite Ranges Using Euclidean Distance Matrices

Iiyama Keidai
Iiyama Keidai
¹Department of Aeronautics and Astronautics, Stanford University, California, United States
- Find this author on Google Scholar
- Find this author on PubMed
- Search for this author on this site
- For correspondence: [email protected]
Neamati Daniel,
Neamati Daniel,
¹Department of Aeronautics and Astronautics, Stanford University, California, United States
- Find this author on Google Scholar
- Find this author on PubMed
- Search for this author on this site
- For correspondence: [email protected]
and Gao Grace
Gao Grace
¹Department of Aeronautics and Astronautics, Stanford University, California, United States
- Find this author on Google Scholar
- Find this author on PubMed
- Search for this author on this site
- For correspondence: [email protected]

NAVIGATION: Journal of the Institute of Navigation
May 2026,
73
navi.764;
DOI: https://doi.org/10.33012/navi.764

Abstract

To support reliable positioning, navigation, and timing services for lunar satellite constellations, this paper proposes a framework for detecting satellite clock phase jumps using dual one-way inter-satellite range measurements. The method models the constellation as a weighted graph and exploits rigidity properties to detect faults. Specifically, satellite clock faults are identified by monitoring the realizability of vertex-redundantly rigid subgraphs through the singular values of geometric-centered Euclidean distance matrices (GCEDMs). We establish the graph-topological conditions required for fault detection and derive the statistical distribution of the GCEDM test statistic under both nominal and faulty conditions. To address sparse link connectivity, we also introduce an ephemeris-augmented formulation that combines measured and estimated ranges and performs fault detection on fault-detectable subgraphs. Simulation results using a terrestrial Global Positioning System constellation and a notional lunar constellation demonstrate that the proposed approach provides robust fault detection under both dense and sparse measurement conditions.

Keywords

1 INTRODUCTION

To meet the growing need for robust positioning, navigation, and timing (PNT) services at the lunar surface and lunar orbits, the National Aeronautics and Space Administration and its international partners are collaborating to develop LunaNet (Israel et al., 2020), a network of networks providing data relay, PNT, detection, and science services. In LunaNet, these services are provided through cooperation among interoperable systems that evolve over time to meet the growing needs for these services, efficiently establishing a reliable, sustainable, and scalable network.

It is crucial to monitor LunaNet navigation satellites for the reliable operation of safety-critical missions, such as lunar landing and human missions on the lunar surface. The quality of lunar navigation signals can be compromised by various system faults, such as clock runoffs (e.g., phase and frequency jumps (Weiss et al., 2010)), unflagged maneuvers, failures in satellite payload signal generation components, and code-carrier incoherence (Walter et al., 2018). While the LunaNet Relay Service Documents state that robustness of the navigation signal should be a key consideration for LunaNet (National Aeronautics and Space Administration, 2022), the specific methodology for monitoring faults on LunaNet satellites has not yet been solidified.

Currently implemented integrity monitoring methods for terrestrial global navigation satellite systems (GNSSs) can be categorized into three types: (1) ground-based monitoring, (2) receiver autonomous integrity monitoring (RAIM), and (3) satellite-based integrity monitoring (SAIM). The first method relies on a network of ground-based monitoring stations at known positions to track and collect GNSS signals, which are then collected and processed at a central computing center to compute corrections and integrity information (Misra et al., 1993). The system fault alerts are provided as integrity information from GNSS augmentation systems, such as a satellite-based augmentation system (van Diggelen, 2009). The second method, RAIM, uses redundant measurements from GNSS satellites to detect and exclude faulty satellites on the receiver side (Misra et al., 1993). The third method, SAIM, computes the integrity information onboard the navigation satellite to detect satellite-related faults such as clock anomalies and ephemeris errors (Xu et al., 2011).

While the first two methods are well established and widely adopted for integrity monitoring of terrestrial GNSSs, their adoption may be more limited for lunar constellations. The first method requires a network of ground-based monitoring stations, which would be expensive and impractical to deploy widely on the lunar surface, especially for near-term constellations. The second method, RAIM, requires at least five satellites in view to detect (six satellites in view to exclude) a faulty satellite, which limits its availability in lunar constellations with a smaller number of navigation satellites than terrestrial GNSSs.

The third method, SAIM, has several advantages compared with the other two methods, including (1) a faster time to alert enabled by eliminating the latency associated with data collection, processing, and uplink on the ground, (2) the use of precise inter-satellite range (ISR) measurements (BeiDou satellites have demonstrated 10-cm accuracy (Zhang et al., 2024)) that are unaffected by the ionosphere and troposphere, and (3) monitoring that can be operated autonomously within the constellation and independently of ground monitoring. The autonomous nature of SAIM is particularly valuable for early lunar deployments, when lunar surface stations will be scarce or absent.

SAIM has been implemented in the BeiDou GNSS to detect clock phase and frequency jumps by comparing the outputs of multiple onboard clocks (Cao et al., 2019; Chen et al., 2022), at the cost of adding extra hardware. Several other works have also proposed SAIM methods utilizing ISR measurements, which detect ephemeris anomalies by directly comparing the ISR measurement and the predicted range from the ephemeris (Rodriguez-Perez et al., 2011; Wolf, 2000; Xu et al., 2011). However, these methods process ISR measurements on a per-satellite basis; as a result, their minimal detectable bias (MDB) is dependent on ephemeris accuracy (Xu et al., 2011). This dependence may prove challenging for lunar constellations, where orbit determination and time synchronization are more difficult than for terrestrial GNSSs. Another potential SAIM method for solving this problem is to centrally aggregate all ISR measurements, form a linear system of residual equations, and apply data-snooping techniques (Baarda, 1968; Imparato et al., 2018) to isolate faulty satellites.

In this paper, we propose a new algorithm to detect phase jumps in the satellite clocks using dual one-way ISR measurements. The proposed method detects faults by identifying biases injected into ISR measurements whenever a satellite’s onboard clock experiences a phase jump between two dual one-way links. We prove that we are able to detect faulty satellites (satellites with clocks that experienced phase jumps) by assessing the realizability of 5-clique (a fully connected graph of 5 nodes) subgraphs, as these subgraphs remain rigid (graphs that cannot be susceptible to continuous flexing are called rigid graphs) after the removal of any vertex from the graph (Motevallian et al., 2015). The concept of the algorithm is illustrated in Figure 1.

FIGURE 1

Concept of the proposed fault detection algorithm

Our proposed method monitors singular values of the geometric-centered Euclidean distance matrix (GCEDM) (Dokmanic et al., 2015) constructed from the ISR measurements, to check if the 5-clique subgraphs are realizable in three-dimensional (3D) space. In particular, we provide analytical expressions of the distribution of the squared fourth singular value under the presence of noise and biases in the range measurement and design hypothesis tests to detect and identify a satellite clock with phase jumps.

The contributions of this paper are summarized below. This work is based on our prior conference paper presented at the 2024 Institute of Navigation GNSS+ Conference (Iiyama, Neamati, & Gao, 2024), but includes key theoretical and simulation updates beyond our prior publication.

We establish the necessary graph topology, specifically, vertex-redundant rigidity, for unambiguous fault detection from ISR measurements.
We prove key rank properties of EDMs and GCEDMs and derive the distribution of the squared fourth singular value under noise and bias.
We develop a hypothesis test based on the fourth singular value of the GCEDM to detect and identify satellite clock phase jumps.
We validate our approach via simulations of a Global Positioning System (GPS) constellation and a notional lunar constellation, analyzing the impact of hyperparameters and fault magnitudes and comparing fault detection and false alarm rates against residual-based baseline methods.

Our method is not intended to replace decades of research in existing fault detection methods. Rather, our method provides new theoretical and practical insights to view the fault detection problem from a different perspective.

This paper is arranged as follows. Section 2 introduces the formulation of the dual one-way ISR measurements. Section 3 describes two baseline clock jump detection methods that use ISR measurement residuals computed from the ephemeris. Section 4 describes the proposed fault detection algorithm based on singular values of the GCEDM. Section 5 provides simulation results for the GPS constellation and a notional lunar constellation. The paper concludes with Section 6.

2 MEASUREMENTS FROM INTER-SATELLITE LINKS

2.1 Dual One-Way ISR Measurement

We consider two satellites establishing a dual one-way inter-satellite link (ISL), as shown in Figure 2, to obtain ISR measurements. In this section, we formulate the ISR measurement model for the dual one-way link, based on work by Zhang et al. (2024).

FIGURE 2

Dual one-way ISL measurement between two satellites. (a) Dual one-way link between the two satellites when there is no phase jump at either satellite. Two pseudorange measurements ρ_ij(t₁) and ρ_ij(t₂) are obtained. (b) Dual one-way link between the two satellites when there is a phase jump at satellite i. When the phase jump occurs at $- Δ t_{1} < t < t_{2}$ , a bias $c δ t_{i}$ is added to the pseudorange measurement ${\tilde{ρ}}_{j i} (t_{2})$ .

Satellite i transmits the signal to satellite j, received at t₁, and satellite j transmits the signal to satellite i, received at t₂.

Let i, j be the indices of the two satellites establishing the link and x_i, x_j ∊ℝ³ be the position vectors of the two satellites. The two observation equations for the dual one-way link can be expressed as follows:

$ρ_{i j} (t_{1}) = ‖ x_{j} (t_{1}) - x_{i} (t_{1} - Δ t_{1}) ‖ + c \cdot [τ_{j} (t_{1}) - τ_{i} (t_{1} - Δ t_{1})] + c \cdot [Δ t_{j}^{r x} + Δ t_{i}^{t x}] + ω_{1}$ 1

$ρ_{j i} (t_{2}) = ‖ x_{i} (t_{2}) - x_{j} (t_{2} - Δ t_{2}) ‖ + c \cdot [τ_{i} (t_{2}) - τ_{j} (t_{2} - Δ t_{2})] + c \cdot [Δ t_{i}^{r x} + Δ t_{j}^{t x}] + ω_{2}$ 2

Here, ρ_ij, ρ_ji represent the measured pseudoranges for the two links, where satellites i and j receive the signal at t₁ and t₂, respectively. Δt₁ and Δt₂ are the signal travel times, τ_i and τ_j are the clock offsets, t^rx and t^tx represent hardware delays of the receiving channel and the transmitting channel, and ω₁ and ω₂ are the measurement noises. We assume that the phase center offset and relativistic effects are corrected; thus, these terms are not included in the equations.

To eliminate common errors between the two links, the two measurements are transformed into the common time t₀ using the following equations:

$ρ_{i j} (t_{0}) = ρ_{i j} (t_{1}) + Δ ρ_{i j} = ‖ x_{j} (t_{0}) - x_{i} (t_{0}) ‖ + c \cdot [τ_{j} (t_{0}) - τ_{i} (t_{0})] + c \cdot [Δ t_{j}^{r x} + Δ t_{i}^{t x}] + ω_{1, 0}$ 3

$ρ_{j i} (t_{0}) = ρ_{j i} (t_{2}) + Δ ρ_{j i} = ‖ x_{i} (t_{0}) - x_{j} (t_{0}) ‖ + c \cdot [τ_{i} (t_{0}) - τ_{j} (t_{0})] + c \cdot [Δ t_{i}^{r x} + Δ t_{j}^{t x}] + ω_{2, 0}$ 4

where $Δ ρ_{i j}$ and $Δ ρ_{j i}$ are the corrections of satellite position and clocks, expressed as follows:

$Δ ρ_{i j} = ‖ x_{j} (t_{0}) - x_{i} (t_{0}) ‖ - ‖ x_{j} (t_{1}) - x_{j} (t_{1} - Δ t_{1}) ‖ + c \cdot [τ_{j} (t_{0}) - τ_{i} (t_{0}) - τ_{j} (t_{1}) + τ_{i} (t_{1} - Δ t_{1})]$ 5

$Δ ρ_{j i} = ‖ x_{i} (t_{0}) - x_{j} (t_{0}) ‖ - ‖ x_{i} (t_{2}) - x_{j} (t_{2} - Δ t_{2}) ‖ + c \cdot [τ_{i} (t_{0}) - τ_{j} (t_{0}) - τ_{i} (t_{1}) + τ_{j} (t_{2} - Δ t_{2})]$ 6

The corrections $Δ ρ_{i j}$ and $Δ ρ_{j i}$ are obtained by using the predicted satellite ephemeris. Therefore, the accuracy of the corrections depends on the accuracy of the satellite velocity and the satellite clock drift of the ephemeris. For BeiDou satellites, the accuracy of the predicted satellite velocity and clock drift is better than 0.1mm/s and 1 x 10^-13 s/s, respectively, resulting in $Δ ρ_{i j}$ , $Δ ρ_{j i}$ <1.0 cm when t₂ – t₁ <3 s (Xie et al., 2022). For LunaNet satellites, the specification for the sum of velocity and clock drift error (converted to velocity via multiplication by the speed of light) in the ephemeris is 1.2 mm/s (for 3 σ); thus, the error of the correction terms can be on the order of centimeters (National Aeronautics and Space Administration, 2022). Whenever further accuracy is required or velocity ephemeris information is not available, we can also utilize the ISR rate measurements between the satellites for correction (Iiyama et al., 2024).

The sum and differences of Equations (5) and (6) can be expressed as follows:

$r_{i j} (t_{0}) = \frac{ρ_{i j} (t_{0}) + ρ_{j i} (t_{0})}{2} = ‖ x_{j} (t_{0}) - x_{i} (t_{0}) ‖ + {\bar{b}}_{i j}^{range} + \frac{ω_{1, 0} + ω_{2, 0}}{2}$ 7

$τ_{i j} (t_{0}) = \frac{ρ_{i j} (t_{0}) - ρ_{j i} (t_{0})}{2} = c \cdot [τ_{j} (t_{0}) - τ_{i} (t_{0})] + {\bar{b}}_{i j}^{clock} + \frac{ω_{1, 0} - ω_{2, 0}}{2}$ 8

where r_ij (t₀) and τ_ij (t₀) are the clock-free and geometry-free observations, respectively. b_ij and b_ji are the bias terms due to hardware delays, given as follows:

${\bar{b}}_{i j}^{range} = c \cdot \frac{Δ t_{j}^{r x} + Δ t_{i}^{t x} + Δ t_{i}^{r x} + Δ t_{j}^{t x}}{2}$ 9

${\bar{b}}_{i j}^{clock} = c \cdot \frac{Δ t_{j}^{r x} + Δ t_{i}^{t x} - Δ t_{i}^{r x} - Δ t_{j}^{t x}}{2}$ 10

The bias terms $b_{i j}^{r a n g e}$ and $b_{i j}^{c l o c k}$ are usually constant over a short term and can therefore be estimated using filtering techniques. In this paper, we assume that these biases can be estimated within certain errors, and the error terms are included in the measurement noise, resulting in the following two equations:

$\begin{array}{l} r_{i j} (t_{0}) = ‖ x_{i} - x_{j} ‖ + ω_{range} & clock-free combination \end{array}$ 11

$\begin{array}{l} τ_{i j} (t_{0}) = c [τ_{j} (t_{0}) - τ_{i} (t_{0})] + ω_{clock} & geometry-free combination \end{array}$ 12

For ISLs in BeiDou, the observed ω_range and ω_clock can be approximated as Gaussian distributions with standard deviations on the order of centimeters.

2.2 Satellite Clock Phase Jumps

Atomic clocks on navigation satellites occasionally experience a sudden jump in the clock phase. When the clock of satellite i experiences a phase jump of δt_i (delay) between t₁ – Δt₁ and t₂, as illustrated in Figure 2(b), the corrupted pseudorange measurement ${\tilde{ρ}}_{j i} (t_{2})$ can be rewritten as follows:

$\begin{array}{l} {\tilde{ρ}}_{j i} (t_{2}) & = & ‖ x_{i} (t_{2}) - x_{j} (t_{2} - Δ t_{2}) ‖ + c \cdot [τ_{i} (t_{2}) - τ_{j} (t_{2} - Δ t_{2})] + c \cdot [Δ t_{i}^{r x} + Δ t_{j}^{t x}] + c \cdot δ t_{i} + ω_{2} \\ = & ρ_{j i} (t_{2}) + c \cdot δ t_{i} \end{array}$ 13

The calculated correction term is as follows:

$\begin{array}{l} Δ {\tilde{ρ}}_{j i} & = & ‖ x_{i} (t_{0}) - x_{j} (t_{0}) ‖ - ‖ x_{i} (t_{2} + δ t_{i}) - x_{j} (t_{2} + δ t_{i} - Δ t_{2}) ‖ \\ + c \cdot [τ_{i} (t_{0}) - τ_{j} (t_{0}) - τ_{i} (t_{1}) + τ_{j} (t_{2} + δ t_{i} - Δ t_{2})] \\ = & Δ ρ_{i j} + δ Δ ρ_{j i} \end{array}$ 14

Here, δΔρ_ji is the error in the correction term due to the error in the time stamping of the reception time t₂ to t₂ + δt_i, which can be represented as follows:

$\begin{array}{l} δ Δ ρ_{j i} & = & (‖ x_{i} (t_{2}) - x_{j} (t_{2} - Δ t_{2}) ‖ - ‖ x_{i} (t_{2} + δ t_{i}) - x_{j} (t_{2} + δ t_{i} - Δ t_{2}) ‖) \\ + c \cdot [τ_{j} (t_{2} + δ t_{i} - Δ t_{2}) - τ_{j} (t_{2} - Δ t_{2})] \\ \approx & [- {\dot{ρ}}_{j i} (t_{2}) + c \cdot {\dot{τ}}_{j} (t_{2} - Δ t_{2})] δ t_{i} \end{array}$ 15

${\dot{ρ}}_{j i}$ represents the range rate of the link between j and i, and ${\dot{τ}}_{j}$ represents the clock drift of the onboard clock at satellite j. Therefore, the corrected pseudorange at t₀ is as follows:

$\begin{array}{l} {\tilde{ρ}}_{j i} (t_{0}) & = & {\tilde{ρ}}_{j i} (t_{2}) + Δ {\tilde{ρ}}_{j i} = (ρ_{j i} (t_{2}) + c δ t_{i}) + (Δ ρ_{i j} + δ Δ ρ_{j i}) \\ \approx & (ρ_{i j} (t_{2}) + Δ ρ_{i j}) + δ t_{i} \cdot [c - {\dot{ρ}}_{j i} (t_{2}) + c {\dot{τ}}_{j} (t_{2} - Δ t_{2})] \\ \approx & ρ_{i j} (t_{0}) + c \cdot δ t_{i} (∵ {\dot{ρ}}_{j i} ≪ c, {\dot{τ}}_{j} ≪ 1) \end{array}$ 16

Similarly, when the clock of satellite j experiences a phase jump of δt_j (delay) between t₁ and t₂ – Δt₂, ρ_ji can be rewritten as follows:

${\tilde{ρ}}_{j i} (t_{0}) \approx ρ_{i j} (t_{0}) - c \cdot δ t_{j}$ 17

Therefore, the clock-free and geometry-free combination considering satellite clock jumps is as follows:

$\begin{array}{l} r_{i j} (t_{0}) = ‖ x_{i} (t_{0}) - x_{j} (t_{0}) ‖ + f_{i} - f_{j} + ω_{range} & clock-free combination \end{array}$ 18

$\begin{array}{l} τ_{i j} (t_{0}) = c [τ_{j} (t_{0}) - τ_{i} (t_{0})] + f_{i} + f_{j} + ω_{clock} & geometry-free combination \end{array}$ 19

$f_{i} = {\begin{cases} \frac{c}{2} \cdot δ t_{i} & clock jump δ t_{i} occurs at satellite i for t within t_{1} - Δ t < t < t_{2} \\ 0 & satellite i's clock is normal \end{cases}$ 20

$f_{j} = {\begin{cases} \frac{c}{2} \cdot δ t_{j} & clock jump δ t_{j} occurs at satellite j for t within t_{1} < t < t_{2} - Δ t \\ 0 & satellite j's clock is normal \end{cases}$ 21

Note that if the phase jump δt_i occurs at $t < t_{1} - Δ t_{1}$ or t > t₂ and δt_j occurs at t < t₁ or $t > t_{2} - Δ t_{2}$ , the errors cancel when ${\tilde{ρ}}_{i j} (t_{0})$ and ${\tilde{ρ}}_{j i} (t_{0})$ are added or subtracted. Consequently, these errors do not appear in r_ij (t₀) or τ_ij (t₀). Therefore, to improve the likelihood of detecting clock jumps, we should either increase the frequency of the ISL measurements or extend the interval between t₁ and t₂. However, note that increasing the interval t₂ – t₁ will lead to larger correction errors (Δρ_ij and Δρ_ji), which, in turn, result in greater ω_range and ω_clock values.

3 BASELINE METHODS: CLOCK FAULT DETECTION USING ISR RESIDUALS

In this section, we introduce two baseline methods for detecting clock jump faults based on existing literature. The two methods use the range residuals, which are computed from the ISL range measurements and the ephemeris information.

3.1 Fault Detection Using the Expected Range from Ephemeris

The first method is based on the ephemeris fault detection method proposed by Rodríguez-Pérez et al. (2011) and Xu et al. (2011). While the original method was proposed for ephemeris fault detection, it can also be utilized for clock phase jump detection.

Let g_ij be the range residual of the i-th ISL at time t₀, which is defined as follows:

$\begin{array}{l} g_{i j} & = & r_{i j} (t_{0}) - ‖ {\hat{x}}_{i} (t_{0}) - {\hat{x}}_{j} (t_{0}) ‖ \\ = & ‖ x_{i} (t_{0}) - x_{j} (t_{0}) ‖ - ‖ {\hat{x}}_{i} (t_{0}) - {\hat{x}}_{j} (t_{0}) ‖ + f_{i} - f_{j} + ω_{range} \end{array}$ 22

where ${\hat{x}}_{i} (t_{0})$ is the estimated position (from ephemeris) of the i-th satellite at time t₀.

The test statistics for the i-th satellite are defined as follows:

$T_{i} = \sum_{j = 1}^{l_{i}} {(\frac{g_{i j}^{2}}{\sqrt{2 σ_{r}^{2} + σ_{m}^{2}}})}^{2}$ 23

where l_i is the number of ISLs that involve the i-th satellite and $σ_{r}^{2}$ and $σ_{m}^{2}$ are the variances of the position estimates in the ephemeris and the variance of ISL range measurement noise (ω_range), respectively.

We define the null hypothesis H_0,i as the absence of faults in the i-th satellite and the alternative hypothesis H_a,i as the presence of a fault in the i-th satellite. The test statistic T_i is distributed as a chi-squared distribution with l_i degrees of freedom under the null hypothesis H₀ . The critical value for the test $X_{α}^{2} (l_{i}, 0)$ for the false alarm rate α is given as follows:

$χ_{α}^{2} (l_{i}, 0) = P [χ^{2} (l_{i}, 0) \leq χ_{α}^{2} (l_{i}, 0)] = α$ 24

where X²(r, λ) is a non-central chi-squared distribution with r degrees of freedom and the non-centrality parameter λ. Based on work by Xu et al. (2011), the testing procedure is as follows:

$Accept H_{0} if {max}_{j \in {1, \dots, n}} T_{i} < \sqrt{χ_{α}^{2} (l_{i}, 0)}$ 25

otherwise:

$Accept H_{k} where k = {argmax}_{j \in {1, \dots, n}} \frac{T_{j}}{\sqrt{l_{j}}}$ 26

The minimum detectable bias, MDB(k), for H_k is computed as follows:

$M D B (k) = \sqrt{λ (α, γ, 1) \cdot (2 σ_{r}^{2} + σ_{m}^{2})}$ 27

The non-centrality parameter $λ (α, γ, l_{i})$ is computed as a function of the false alarm rate α and detection power γ (equal to the missed-detection rate of 1 – γ), by solving for the following equation:

$P_{M D} = P [χ^{2} (1, λ) \leq χ_{α}^{2} (1, 0)] = 1 - γ$ 28

Because of the correlation between g_i,j. and g_i,k (i ≠ k) via the common ephemeris error, the distribution of the test statistic under the null and alternate hypotheses will not exactly follow a chi-squared and non-central chi-squared distribution, respectively. Based on this observation, Xu et al. (2011) proposed to compute a conservative threshold for the test statistic via numerical simulation by using the correlation between two links g_i,j. and g_i,k as follows:

$ρ (g_{i j}, g_{i k}) \leq \frac{σ_{r}^{2}}{2 σ_{r}^{2} + σ_{m}^{2}}$ 29

In this paper, we used Imhof’s method (Imhof, 1961) to compute the critical value for the test statistic. The detailed procedure is explained in Appendix A.

The advantage of this “ephemeris comparison” approach is that the fault detection is fully distributed among the satellites. In addition, the method does not require a particular geometry of links to perform fault detection: if the ephemeris accuracy is sufficiently high, the fault can be detected by directly comparing the expected and measured range. The test is also able to detect faults in the ephemeris information if they are present. However, the disadvantage of this fault detection method is that the MDB is dependent on the ephemeris error σ_r, as shown in Equation (27).

3.2 Data Snooping with Baarda’s w-Test Statistic

The test in Section 3.1 does not use all of the ISL measurements available in the constellation. By testing the null hypothesis on a complete set of measurements from all ISLs, we can reduce the MDB when the measurement errors are smaller than the ephemeris errors.

Let n be the number of satellites in the constellation, $X \in ℝ^{3 n}$ be the flattened vector of spacecraft positions at time t₀, and $Y \in ℝ^{m}$ be the vector of m ISR measurements at time t₀. We also define b as the bias in the range measurement (equal to cδt/2 in Equations (20) and (21)) from one faulty satellite f that experienced a clock phase jump δt.

The mapping between X and Y can be represented as the following nonlinear function h:

$Y = h (X) + c_{f} b + ω, ω \sim N (0_{m}, Σ_{y y})$ 30

where $\sum_{y y} = σ_{m}^{2} I_{m}$ is the covariance of the range measurements. c_f is a known vector that maps the faulty satellites to the corrupted range measurements, as follows:

$c_{f}^{(i)} = {\begin{cases} 1 & (satellite f is the transmitter of link i) & \\ (clock phase jump at satellite f occurred during link i) \\ - 1 & (satellite f is the receiver of link i) & \\ (clock phase jump at satellite f occurred during link i) \\ 0 & otherwise \end{cases}$ 31

where $c_{f}^{(i)}$ is the i-th element of cf. By linearizing the above equation around the ephemeris position $\bar{X}$ , we obtain the following:

$y = Y - h (\bar{X}) \approx H (X - \overset{―}{X}) + c_{f} b + ω = H x + c_{f} b + ω$ 32

where $y = Y - h (\bar{X})$ , $x = X - \bar{X}$ , and $H \in ℝ^{m \times 3 n}$ is the Jacobian matrix.

For the test, we consider n + 1 hypotheses: the null hypothesis H₀ where model errors are absent, and the alternative hypotheses H₁, …, H_n, which correspond to the faulty satellite being f = 1, f = 2, …, f = n, respectively. Baarda’s w-test statistic for testing H₀ against a single alternative H_k(k = 1, …, n) is given as follows (de Jong & Teunissen, 2014):

$w_{k} = \frac{c_{k}^{⊤} Σ_{y y}^{- 1} P_{A}^{⊥} y}{\sqrt{c_{k}^{⊤} Σ_{y y}^{- 1} P_{A}^{⊥} c_{k}}}$ 33

Here, the projector $P_{A}^{⊥}$ is defined as follows:

$P_{A}^{⊥} = I_{m} - H {(H^{⊤} Σ_{yy}^{- 1} H)}^{- 1} H^{⊤} Σ_{yy}^{- 1}$ 34

and:

$c_{k}^{(i)} = {\begin{cases} 1 & satellite k is the transmitter of link i \\ - 1 & satellite k is the receiver of link i \\ 0 & otherwise \end{cases}$ 35

Here, we assume that the clock phase jump at satellite k corrupted all of the ISLs that involve satellite k. The testing procedure is summarized as follows (Imparato et al., 2018):

$Accept H_{0} if {max}_{j \in {1, \dots, n}} | w_{j} | < \sqrt{χ_{α}^{2} (1, 0)}$ 36

otherwise:

$Accept H_{k} if | w_{k} | = {max}_{j \in {1, \dots, n}} | w_{j} | \geq \sqrt{χ_{α}^{2} (1, 0)}$ 37

The MDB for H_k is computed as follows (de Jong & Teunissen, 2014):

$MDB (k) = \sqrt{\frac{λ (α, γ, 1)}{c_{k}^{⊤} Σ_{y y}^{- 1} P_{A}^{⊥} c_{k}}}$ 38

Iterating over n alternate hypotheses using Baarda’s w-test statistic (which is called data snooping) enables one to detect the existence of biases by considering the effect of biases on the full set of observations. Each test is optimal for any bias magnitude b, whenever the assumptions of the test are true (Lehmann & VoB-Böhme, 2017). At the same time, this test also has some limitations:

In c_k, we assume that the biases are injected into all ISR measurements that contain a satellite with a clock phase jump. However, this assumption is not always true because clock jumps may have occurred outside the transmission and reception time for some links (the second condition in Equation (31) is not satisfied), based on the satellite link schedule (Iiyama et al., 2024). This situation would result in having smaller fault links than the assumed model (c_k) and having reduced detection power. If we want to rigorously handle this uncertainty of c_k, we would need to consider all possible combinations of the links as alternate hypothesis, but this approach would be computationally intractable.
As further discussed in Section 4.1, this method fundamentally detects faults by checking the consistency among redundant measurements. Therefore, a sufficient level of measurement redundancy is required. When the measurement graph is sparse and redundancy cannot be guaranteed, certain fault modes become unobservable, and the method may fail to detect them. In these cases, it may be necessary to rely on direct comparison with ephemeris-based ranges, which do not depend on internal measurement consistency.

Despite these limitations, there remains significant room for methodological improvement to mitigate them. For example, instead of formulating n hypotheses at the satellite (vertex) level, one could construct hypotheses at the edge (range) level, explicitly modeling which individual ISRs are corrupted. A subsequent inference layer could then estimate the faulty satellites from the set of detected faulty edges. Such a hierarchical formulation would better accommodate partial-link corruption. We leave the development and evaluation of these extensions to future work.

4 PROPOSED CLOCK FAULT DETECTION ALGORITHM USING EUCLIDEAN DISTANCE MATRICES

In this section, we propose a satellite clock fault detection algorithm based on rigid graph theory. The method models the set of ISR measurements as a weighted graph and determines whether the set of range measurements is realizable (embeddable) in 3D Euclidean space. Instead of analyzing the realizability of the entire graph, the proposed method analyzes a set of subgraphs to detect faults. The method directly processes the range measurements and therefore does not require any prior information about the satellite positions or clocks.

This section is structured as follows. First, in Section 4.1, we discuss the topology required to detect clock faults based on graph embeddability. In particular, we prove that the graph must be vertex-redundantly rigid (the minimum graph with such a property in three dimensions is a fully connected graph of five nodes) to detect biases in range measurements from the graph’s embeddability. Next, in Section 4.2, we show that the embeddability of a graph can be determined by the properties of the GCEDM, without solving for the position of the nodes. Finally, in Section 4.3.2, we propose a fault detection test based on the GCEDM distributions of fully connected subgraphs in the graph.

4.1 Required Graph Topology for Fault Detection

4.1.1 Definitions

Consider a graph of n nodes, which correspond to satellites. Let $x_{i} \in ℝ^{d}, (i \in {1, \dots, n})$ be the position of the i-th node and $X \in ℝ^{d \times n}$ be a matrix collecting these points. In this paper, we are interested in the 3D case, d = 3. In this section and the latter sections, we use the following notation for range measurements between the two satellites i and j, instead of Equation (18):

$r_{i j} = ‖ x_{i} - x_{j} ‖ + f_{i j} + w_{i j}$ 39

where f_ij = f_i – f_j is the bias of the range measurements added by the clock phase jumps and w_ij is the measurement noise.

Based on the graph and observed ranges r_ij, we define a weighted graph $G = 〈 V, E, W 〉$ . A realization (embedding) of G is d-dimensional Euclidean distance space, ℝ^d, defined as a function Φ that maps V into ℝ^d, such that for each edge $e = (v_{i}, v_{j}) \in E$ , $W (e) = r_{i j}$ (Jian et al., 2010). A weighted graph $G = 〈 V, E, W 〉$ is called d-embeddable when there is a realization (embedding) $Φ : V \to ℝ^{d}$ that maps the edges to points (x₁, …, x_n) in ℝ^d space.

We define a “faulty satellite” as a satellite that experiences a clock phase jump and corrupts at least one range measurement containing itself. A weighted graph is fault-disprovable if and only if the embeddability of G implies that it contains no faulty satellites and the non-embeddability of G implies that it contains at least one faulty satellite.

A graph is called rigid if it has no continuous deformation other than rotation, translation, and reflection while preserving the relative distance constraints between the vertices (Laman, 1970). Otherwise, the graph is called flexible. A graph is called k-vertex-redundantly rigid if it remains rigid after the removal of any (k – 1) vertices.

Flexible, rigid, and two-vertex redundantly rigid graphs in ℝ2. (a) Flexible graph. (b) Rigid graph. (c) Two-vertex-redundantly rigid graph.

FIGURE 3

Flexible, rigid, and two-vertex redundantly rigid graphs in ℝ². (a) Flexible graph. (b) Rigid graph. (c) Two-vertex-redundantly rigid graph.

4.1.2 Embeddability of Graphs with Faulty Satellites

In this section, we prove the sufficient condition for a graph to be fault-disprovable in ℝ^d. For the proofs of this subsection only, we ignore the noise added to the measurements.

Proposition 4.1. Suppose that there is at most k – 1 faulty satellites on a graph $G = 〈 V, E, W 〉$ that is k-vertex (k ≥ 2) redundantly rigid in ℝ^d. If G is d-embeddable, then G contains no faulty vertices (satellites) with probability 1.

Proof. See Appendix B.

Proposition 4.2. Given a weighted graph G, if G is k-vertex-redundantly rigid (k ≥ 2), then G is fault-disprovable at probability 1, if there are at most k – 1 faulty satellites.

Proof. For an arbitrary graph G, if G is not embeddable, then G contains at least one faulty satellite. Therefore, from Proposition 4.1, if G is k-vertex-redundantly rigid (k ≥ 2), then G is fault-disprovable at probability 1, if there are at most k – 1 faulty satellites.

4.2 Properties of the Euclidean Distance Matrix

In the previous section, we proved that by using vertex-redundantly rigid graphs, we can detect faults by analyzing their embeddability. Analyzing the embeddability of an arbitrary graph is known to be NP-hard (Hendrickson, 1992), but if the graph is fully connected, we can evaluate it by analyzing the ranks of the GCEDM. Note that fully connected graphs of more than five nodes are 2-vertex-redundantly rigid in ℝ³.

In this section, we introduce and prove several properties of a Euclidean distance matrix (EDM) and GCEDM constructed from the observed ranges between satellites. We also show how we can use the singular values of the GCEDM to evaluate whether a (fully connected) (sub)graph is embeddable in ℝ³.

4.2.1 Geometric-Centered Euclidean Distance Matrix

Using the set of observed ranges r_ij, we construct an EDM in which the elements are equivalent to the square of the observed ranges $(r_{j i}^{2})$ .

Let $D^{n, d, m} \in ℝ^{n \times n}$ be an EDM without measurement noise (W_ij = 0) and with m nodes out of n nodes being faulty. For example, a noiseless EDM with six satellites (in 3D space) with two faulty satellites is denoted as D^6,3,2. Similarly, we define ${\tilde{D}}^{n, d, m} \in ℝ^{n \times n}$ as an EDM with measurement noise and with m nodes out of n nodes having a fault. A GCEDM G^n,d,m (or ${\tilde{G}}^{n, d, m}$ for a noisy EDM ${\tilde{D}}^{n, d, m}$ ) is constructed from the EDM via the following operation:

$G^{n, d, m} = - \frac{1}{2} J^{n} D^{n, d, m} J^{n}$ 40

Here, Jⁿ is the geometric-centering matrix, as follows:

$J^{n} = I^{n} - \frac{1}{n} 11^{⊤}$ 41

where 1 ∊ ℝⁿ is a vector of ones. When no fault or noise is present (m = 0), G^n,d,0 will be positive semi-definite, and its rank will satisfy rank(G^n,d,0) = rank(X^⊤X) ≤ d (Dokmanic et al., 2015).

4.2.2 Rank of EDMs with Faulty Satellites

Knowles and Gao (2023) observed that there are more than three non-zero singular values when the GCEDM is constructed from an EDM with faulty satellites. An example is provided in Figure 4. However, mathematical proofs were not provided in their paper. In this section, we provide and prove several properties of EDMs and GCEDMs that are corrupted by faults and measurement noise.

FIGURE 4

Log plot of the singular values of an example GCEDM for n = 12, d = 3. (a) Singular values of a GCEDM with one faulty satellite. The fourth and fifth singular values increase compared with the non-fault case when the fault magnitude is sufficiently larger than the noise magnitude. (b) Singular values of a GCEDM with two faulty satellites. The fourth, fifth, sixth, and seventh singular values increase compared with the non-fault case when the fault magnitude is sufficiently larger than the noise magnitude.

The singular values are ordered in descending order of their magnitude. The blue circles show the singular values when there is no fault or noise (G^12,3,0), the orange squares show the singular values when there is noise (of scale σ_w = 10^–6) added to the range measurements ( ${\tilde{G}}^{12, 3, 0}$ ), and the green crosses show the singular values when both noise and one or two faults (of scale $\bar{b} = 10^{- 3}$ ) are added ( ${\tilde{G}}^{12, 3, 2}$ ). The 12 points are identical for all cases, which are sampled randomly inside a unit cube. When noise is present, the singular values increase except for the last singular value, which is close to zero (theoretically zero). When m faults with magnitudes that are larger than the noise magnitude are added, the singular values up to the d + 2m index increase compared with the non-fault case.

Proposition 4.3. The rank of the EDM D^n,d,m satisfies the following:

$rank (D^{n, d, m}) \leq min (d + 2 + 2 m, n)$ 42

Proof. See Appendix C.

Proposition 4.4. The rank of the GCEDM $G^{n, d, m} = - \frac{1}{2} J^{n} D^{n, d, m} J^{n}$ satisfies the following:

$rank (G^{n, d, m}) \leq min (d + 2 m, n - 1)$ 43

Proof. See Appendix D.

Corollary 4.4.1. The GCEDM ${\tilde{G}}^{n, d, m} = - \frac{1}{2} J_{n} {\tilde{D}}^{n, d, m} J_{n}$ constructed from an EDM with the edges corrupted by Gaussian noise almost surely has the following rank:

$rank ({\tilde{G}}^{n, d, m}) = n - 1$ 44

Proof. See Appendix E.

4.2.3 Distribution of the Fourth Singular Value of the GCEDM

The propositions proved in the previous sections indicate that we can detect faulty satellites by observing the increase in the fourth singular value (for the 3D case) of the GCEDMs constructed from ISR measurements. We use the following test statistic γ_test to monitor whether a faulty satellite exists in a given subgraph:

$γ_{test} = λ_{4}^{2}$ 45

where $λ_{\bar{i}}$ is the i-th singular value of the GCEDM. If the ISR measurements are completely noiseless, we can detect whether a faulty satellite exists within the graph by checking whether γ_test of the GCEDMs is not 0, because GCEDMs have rank 3 when no fault exists. However, when the ISR measurements are noisy, the fourth singular value increases regardless of faults, as shown in Figure 4 and Corollary 4.4.1. Therefore, to detect faults in the presence of noise, we must set a threshold based on the singular value distributions of the noisy (but non-fault) GCEDMs to separate the fault from the noise.

For n = 5, we can compute the distribution of the test statistic analytically. The distribution of the test statistic γ_test/s² follows a centralized chi-squared distribution with one degree of freedom, as follows:

$\frac{γ_{test}}{s^{2}} \sim χ^{2} (1, 0) for n = 5$ 46

Here, the squared scale of the test statistic s² is given by the following equation:

$s^{2} = \sum_{a = 1}^{2} \sum_{b = 1}^{2} \sum_{i = 1}^{5} \sum_{j = 1}^{5} {(σ_{m} {\tilde{D}}_{i j})}^{2} \cdot {({\hat{U}}_{i, a} {\hat{V}}_{j, b} + {\hat{U}}_{j, a} {\hat{V}}_{i, b})}^{2}$ 47

$\hat{U} = J_{5} {\tilde{U}}_{:, 4 : 5} \in R^{5 \times 2} {\tilde{U}}_{:, 4 : 5} : 4 th and 5th column of \tilde{U} \in R^{5 \times 5}$ 48

$\hat{V} = J_{5} {\tilde{V}}_{:, 4 : 5} \in R^{5 \times 2} {\tilde{V}}_{:, 4 : 5} : 4 th and 5th column of \tilde{V} \in R^{5 \times 5}$ 49

$\tilde{U} \tilde{Σ} {\tilde{V}}^{⊤} = SVD (\tilde{G}) SVD : singular value decomposition, γ_{test} = {\tilde{Σ}}_{4, 4}^{2} = {(λ_{4})}^{2}$ 50

Above, $\tilde{G}$ and $\tilde{D} \in ℝ^{5 t i m e s 5}$ are the (noisy) GCEDM and the distance matrix constructed from the measured ranges, respectively. The derivation is given in Appendix F. The squared scale of the distribution is proportional to the noise variance $σ_{m}^{2}$ (Figure 5). When the variance of noise is different for each measurement, Equation (47) can be modified to compute s² by replacing σ_m with (σ_m)_ij (the noise for each range measurement).

Distribution of the test statistic γtest and its scale s2 (squared scale of λ4) for different noise levels σm. (a) 3D points used for singular value simulation. (b) Distribution of the squared fourth singular value for different noise levels σm. (c) Scale of the chi-squared distribution s, which is proportional to the noise magnitude σm2. The distribution of the test statistic γtest follows a chi-squared distribution with one degree of freedom. The scale s2 is proportional to the noise magnitude σm2.

FIGURE 5

Distribution of the test statistic γ_test and its scale s² (squared scale of λ₄) for different noise levels σ_m. (a) 3D points used for singular value simulation. (b) Distribution of the squared fourth singular value for different noise levels σ_m. (c) Scale of the chi-squared distribution s, which is proportional to the noise magnitude $σ_{m}^{2}$ .

The distribution of the test statistic γ_test follows a chi-squared distribution with one degree of freedom. The scale s² is proportional to the noise magnitude $σ_{m}^{2}$ .

When there is a clock jump at satellite S_f, the distribution of the test statistic ${(γ_{t e s t})}_{s_{f}}$ follows a non-central chi-squared distribution with one degree of freedom, with its squared scale s² given by the same equation as above and non-centrality parameter $λ_{s_{f}}$ :

$\frac{{(γ_{test})}_{s_{f}}}{s^{2}} \sim χ^{2} (1, λ_{s_{f}}) for n = 5$ 51

The non-centrality parameter $λ_{s_{f}}$ is given by the following equation:

$λ_{s_{f}} = \frac{1}{s^{2}} {‖ {\hat{U}}^{⊤} (F \circ D) \hat{V} ‖}_{F}^{2}$ 52

$F = {\begin{cases} f_{i j} & i = s_{f} or j = s_{f} \\ 0 & otherwise \end{cases}$ 53

where $‖ \cdot ‖_{F}$ represents the Frobenius norm. This result indicates that when all biases added by faults have the same magnitude, the non-centrality parameter λ is proportional to the square of the fault magnitude, as shown in Figure 6 and Figure 7. The derivation is given in Appendix G.

FIGURE 6

Distribution of the test statistic γ_test for different fault magnitudes and faulty satellites

The top row shows the points in 3D space, where the red points and edges indicate the faulty satellite and its connected edges, respectively. The distribution of the test statistic γ_test follows a non-central chi-squared distribution with one degree of freedom. The non-centrality parameter λ increases as the fault magnitude increases.

The non-centrality parameter λ of the distribution of the test statistic is proportional to the square of the fault magnitude b¯2.

FIGURE 7

The non-centrality parameter λ of the distribution of the test statistic is proportional to the square of the fault magnitude ${\bar{b}}^{2}$ .

Using Equations (46) and (51), we can compute the MDB for satellite i as follows:

$MDB (i) = \sqrt{\frac{\bar{λ}}{λ_{i}} \cdot s^{2}} where P [χ^{2} (1, \bar{λ}) \leq χ_{α}^{2} (1, 0)] = 1 - γ$ 54

where α is the false alarm rate, γ is the detection power, s² is the squared scale provided in Equation (47), and λ_i is the non-centrality parameter for the satellite i fault case, as provided in Equation (52).

Figure 8 shows the computed MDBs for the two baseline methods (ephemeris comparison and data snooping) and the EDM-based fault detection, using the same geometry as in Figures 5 and 6. The MDB of the data snooping has the lowest value owing to its optimality, whereas the MDB of the ephemeris comparison method largely depends on the ephemeris accuracy. The MDB of the EDM-based detection depends on the relative geometry of each point with respect to other points, which affects the non-centrality parameter $λ_{s_{f}}$ .

FIGURE 8

Comparison of MDBs for different fault detection methods

The topology of the points and links is the same as in Figure 5. The detection power γ is fixed to 0.8, and the noise level is fixed to 1e-8. λ_s in the figure represents λ_sf in Equation (53).

If desired, the test statistic can be extended for n ≥ 6 as follows:

$λ_{test} = \sum_{i = 4}^{n - 1} λ_{i}^{2}$ 55

The readers are referred to Appendix H for the distribution of this test statistic under nominal conditions and the existence of faults.

4.3 Fault Detection Algorithm Using the EDM

4.3.1 Clique Listing

To construct an EDM and GCEDM, we need a set of range measurements between all pairs of n satellites. However, owing to occultation by the planetary body and attitude constraints, the satellites in the constellation are not fully connected with ISRs in most cases. Therefore, the first step of the fault detection algorithm is to find a set of k-clique subgraphs, which are fully connected subgraphs of k satellites. For the graph to be 2-vertex-redundantly rigid and to use the derived distribution in Equation (46), we are interested in cliques with k = 5. An example of 5-clique subgraphs is shown in Figure 9.

FIGURE 9

Examples of 5-cliques within a graph of 13 nodes that are not fully connected. (a) Topology of the entire graph. (b) 5-clique subgraph 1. (c) 5-clique subgraph 2. (d) 5-clique subgraph 3.

Various exact algorithms for finding all k-cliques have been proposed (Li et al., 2020). In this paper, we use the Chiba–Nishizeki algorithm (Arbo) owing to its simplicity of implementation (Chiba & Nishizeki, 1985). For each node v, Arbo expands the list of k-cliques by recursively creating a subgraph induced by v’s neighbors. The readers are referred to the work by Li et al. (2020) for the pseudocode and summary of the k-clique listing algorithm. Note that the clique finding algorithm does not have to be executed online. Because the availability of the ISRs can be predicted beforehand when scheduling the ISLs using a coarse predicted orbit, we can estimate the topology of the ISLs for future time epochs. We can compute the list of k-cliques of the predicted topologies for future time epochs (e.g., one orbit) and uplink them to satellites intermittently. If some of these links were actually not available in orbit for some reason, we can remove the k-cliques that contain the missing links from the list.

4.3.2 Online Fault Detection

We define a “faulty satellite” as a satellite that has a clock jump and a “fault-free satellite” as a satellite that does not have a clock jump. We design a test where we consider the following hypotheses:

$\begin{array}{l} H_{0} : All satellites are fault-free (null hypothesis) \\ H_{i} : Satellite i is a faulty satellite \end{array}$ 56

Let $L_{t}^{G}$ be the list of all 5-cliques of the satellite network as a graph at time step t, and let N_s be the total number of satellites. We propose an online fault detection algorithm that uses $L_{t}^{G}$ . The procedure of the detection algorithm is as follows.

Algorithmml:

1. First, we compute the scaled test statistic $γ_{t e s t, G'} λ_{4}^{2} / s^{2}$ for each 5-clique subgraph G' in $L_{t}^{G}$ , using Equation (47).

2. For each satellite i in the graph, repeat the following steps (a) to (c):

(a) Compute the sum of all of the scaled test statistics for subgraphs that do not include satellite i:

$γ_{test}^{i} = \sum_{G^{'} \in L_{t}^{G}, i \notin G^{'}} γ_{test, G^{'}}$ 57

(b) Compute the threshold ${\bar{γ}}_{test}^{i}$ for the sum of the scaled test statistics $γ_{t e s t}^{i}$ for a given false alarm rate α and (user-defined) threshold margin η_α:

${\bar{γ}}_{test}^{i} = η_{α} \cdot χ_{α}^{2} (N_{G, i}, 0)$ 58

where N_G,i is the number of subgraphs that do not include satellite i,η_α is the threshold margin, and $X_{α}^{2} (N_{G, i}, 0)$ is the α-quantile of the central chi-squared distribution with N_G,i degrees of freedom.

$γ_{test}^{i} = \frac{γ_{test}^{i}}{{\bar{γ}}_{test}^{i}}$ 59

3. If $γ_{test}^{i} < 1$ for each satellite i, we accept the null hypothesis; otherwise, we reject the null hypothesis. We identify the faulty satellite S_f as the one with the smallest $γ_{test}^{i}$ :

$s_{f} = {argmin}_{i} γ_{test}^{i}$ 60

The proposed fault detection algorithm has two hyperparameters that affect its performance.

False Alarm Rate α: This is the probability of a false alarm, which is the probability of detecting a satellite as faulty when it is actually fault-free.

Threshold Margin η_α: This margin inflates the threshold ${\bar{γ}}_{test}$ to $η_{α} \cdot {\bar{γ}}_{test}$ , to account for the increase in the variance of the sum of the test statistics, due to the correlation between the test statistics of the subgraphs.

We must set the threshold margin η_α to a value larger than 1, to account for the increase in the variance of the distribution of $γ_{t e s t}^{i}$ , owing to the correlation between the test statistics of the subgraphs. The correlation between the test statistics of the subgraphs is caused by the fact that the same set of range measurements is used to compute the test statistics between different subgraphs. We compute the test statistics for subgraphs that do not contain satellite i, instead of subgraphs that contain i, to avoid having strong correlations between the test statistics of the subgraphs, which would require a larger threshold margin η_α.

The above test only considers a maximum of one faulty satellite. To consider multiple faulty satellites, the test could potentially be expanded by running the algorithm iteratively and removing the detected faulty satellite from the graph after each iteration until no fault is detected (Knowles & Gao, 2024). Alternatively, we can consider $N_{s} C_{N_{f}}$ combinations of faulty satellites, where N_f is the number of faulty satellites, and apply the test for each combination. The detection performance in the case of multiple faults will be assessed in future work.

One practical consideration of the proposed approach is that, for satellite constellations, the fourth singular value of the GCEDM can become extremely small (~10^–7) because of the large inter-satellite distances (on the order of 10⁶ m) compared with range measurement accuracy (on the order of sub-meters to meters). This scenario can lead to values that approach machine precision when operating in single precision, necessitating the use of double-precision arithmetic.

4.3.3 Fault-Detectable Graphs with Ephemeris Augmentation

For scenarios in which the links are sparse or the number of satellites is small, the graph may lack a sufficient number of 5-cliques to perform robust fault detection using only ISR measurements. To mitigate this issue, we can augment the graph by “filling in” missing edges using range estimates derived from satellite ephemerides. This hybrid approach utilizes a combination of measured edges (ISR) and computed edges (ephemeris) to construct redundantly rigid subgraphs for fault detection.

However, a fault (specifically, a clock phase jump) appears only in the physical ISR measurements and does not affect the ranges computed purely from ephemeris data. Consequently, for an augmented graph to reliably detect faults, the topology of the measured edges must satisfy a specific connectivity condition.

Proposition 4.5. Consider a redundantly rigid weighted graph $G = 〈 V, E, W 〉$ , where the edge set is partitioned into measured edges E_m (ISR measurements) and computed edges E_e (derived from ephemeris), such that E = E_m u E_e. Let a fault at satellite v ∊ V be defined as a bias that corrupts the weights of all incident edges in E_m. The graph G is fault-disprovable (i.e., the graph becomes unrealizable under a fault) if and only if every vertex v ∈ V is incident to at least one measured edge in E_m.

Proof. See Appendix I.

We denote the k-node subgraphs that satisfy the connectivity condition in Proposition 4.2 (i.e., every vertex is incident to at least one measured edge) as fault-detectable k-node subgraphs.

For scenarios in which the graph is sparse or the number of satellites is small, the standard 5-clique list $L_{t}^{G}$ described in Section 4.3.2 may be empty or insufficient for robust detection. In such cases, we replace $L_{t}^{G}$ with the list of fault-detectable 5-node subgraphs, denoted as $L_{t}^{FD}$ . These subgraphs are constructed by identifying sets of 5 nodes and “filling in” missing measured edges with computed edges derived from the ephemeris to form fully connected 5-cliques.

To account for the difference in error characteristics between measured and computed edges, we modify the computation of the test statistic variance s². In Equation (50) (and the generalized form in Equation (117)), the uniform measurement noise variance ${(σ_{m})}_{i j}^{2}$ is replaced by an edge-dependent variance ${(σ_{link})}_{i j}^{2}$ :

${(σ_{link})}_{i j}^{2} = {\begin{cases} σ_{m}^{2} & if (i, j) \in E_{m} & (measured edge) \\ 2 σ_{r}^{2} & if (i, j) \in E_{e} & (computed edge) \end{cases}$ 61

where $σ_{m}^{2}$ is the variance of the ISR measurement noise and $σ_{r}^{2}$ is the variance of the satellite position error (assuming independent position errors, the variance of the computed range error is approximately 2 $σ_{r}^{2}$ ).

This substitution drastically increases the number of subgraphs available for fault detection, particularly in sparse constellations or during orbital phases with limited visibility. It will also enable the detection of ephemeris faults. However, this improved availability comes at a cost: the MDB for these subgraphs becomes sensitive to the magnitude of ephemeris errors (σ_r), rather than being driven solely by the typically smaller measurement noise (σ_m).

5 SATELLITE FAULT DETECTION SIMULATION

5.1 Simulation Configuration

5.1.1 Orbits and Links

To validate the proposed fault detection algorithm, we consider the following two constellations.

GPS constellation of 31 satellites: The satellite orbits are initialized using two-line element (TLE) data from 2025/4/3 from the CelesTrak website (Kelso, 2024). The 3D plots of the constellation and the available ISLs at t = 0 are shown in Figure 11(a).
A notational lunar constellation of 9 satellites: We consider a constellation of 9 satellites in three orbital planes of elliptical lunar frozen orbit (ELFO) with a 30-h orbital period, as a potential extension of the Lunar Communications Relay and Navigation Systems (Ryden & Volle, 2025). The orbital elements of the constellation are shown in Table 1, and the 3D plots of the constellation and the available ISLs at t = 0 are shown in Figure 12(a).

For both constellations, the orbits are propagated in the two-body propagator without a consideration of perturbations. The availability of the links is calculated by assuming that the links are available when the link between the satellites is not occulted by the central body and the angle between the line of sight and the satellite-body center vector is less than ϕ_max, as illustrated in Figure 10. The cutoff angle ϕ_max is set to 60° for the Earth case (Zhang et al., 2024) and 90° for the Moon case. To compute body occultation, we set an additional mask of 1000 km for the Earth case and 100 km for the Moon case to avoid the link being corrupted by ionospheric delays or terrain blockage, respectively. We assume that we can obtain range measurements from all available links at each time step. The scheduling of the links will be investigated in future work.

FIGURE 10

The visibility of the link between two satellites is calculated by assuming that the link is visible when it is not occulted by the central body (including the atmospheric mask) and the angle between the line of sight and the satellite-body center vector is less than ϕ_max. w.r.t.: with respect to

A scatterplot of the number of 5-cliques in the constellation containing each satellite is shown in Figure 11(b) and Figure 12(b) for the Earth and Moon cases, respectively. For the lunar constellation, we observe that the number of available 5-cliques is significantly lower (below 30), limiting the use of EDM-based fault detection from ISR alone. Therefore, for the lunar case, we used the ephemeris augmentation described in Section 4.3.3 to increase the number of subgraphs that can be used for fault detection. The number of fault-detectable 5-node subgraphs is shown in Figure 12(c).

FIGURE 11

A GPS constellation around Earth with 31 satellites. (a) Orbital planes with satellite position (black points) and links (orange lines) at the initial time epoch. The constellation consists of 31 satellites in six orbital planes. (b) Scatterplot of the number of 5-cliques for each satellite in the constellation over a full orbit. The number of 5-cliques that contain pseudorandom noise (PRN) 1, 10, and 15 are shown in blue, orange, and green, respectively. Each satellite has more than 30 self-containing subgraphs throughout the orbit.

The orbit is initialized using TLE data from CelesTrak.

View this table:

TABLE 1 Orbital Elements of the Waker/ELFO Hybrid Lunar Constellation (Case 2)

FIGURE 12

A lunar constellation around the Moon with 9 satellites. (a) Orbital planes with satellite position (black points) and links (orange lines) at the initial time epoch. The constellation consists of 9 satellites in the Walker constellation with three planes. (b) Scatterplot of the number of self-containing 5-cliques for each satellite in the constellation. The numbers of 5-cliques that contain PRN 1, 2, and 3 are shown in blue, orange, and green, respectively. (c) Scatterplot of the number of self-containing fault-detectable 5-satellite subgraphs for each satellite in the constellation. All satellites have the same number of self-containing fault-detectable subgraphs.

5.1.2 Parameters

We tested different combinations of the following detection and simulation parameters:

Test methods: “ephemeris comparison,” “data snooping,” and “EDM”
Fault magnitudes: $\bar{b}$ = 2,4,..., 20 m
Ratio of fault links (the probability that the link connected to a faulty satellite is biased): r_f = 0.2,1.0
Number of faulty satellites: N_f = 0, 1
Ephemeris accuracy: σ_m = {1 m, 2 m} for Earth, σ_m = {2 m, 3 m, 4 m} for Moon
Input probability of false alarm for the fault detection algorithms:
α = 0.001, 0.002, 0.003,0.005,0.008,0.013,0.022,0.036, 0.06,0.1

For each of these cases, we ran 5000 Monte Carlo analyses with randomly sampled initial times (within one orbital period) and a set of faulty satellites. For all cases, we set the standard deviation of the two-way noise measurement error to σ_m = 0.5 m, based on the results of the BeiDou constellation (Zhang et al., 2024). An inflation rate of η_α = 1.5 is used for the Earth scenario, and η_α = 5.0 is used for

the lunar scenario.

5.1.3 Evaluation Metrics

We compute the following metrics to evaluate the performance of the fault detection algorithm:

$Probability of Missed Detection P_{m d} = \frac{F N}{T P + F N}$ 62

$Probability of False Alarm P_{f a} = \frac{F P}{F P + T N}$ 63

where TP, FN, FP, and TN are the total number of true positives (faulty satellite classified as faulty), false negatives (faulty satellite classified as non-faulty), false positives (non-faulty satellite classified as faulty), and true negatives (non-faulty satellite classified as non-faulty). A lower P_md is desired to avoid missing faulty satellites, and a lower P_fa is desired to prevent false alarms, detecting normal satellites as a fault.

5.2 Fault Detection Results

5.2.1 Case 1: GPS Constellation

Figure 13 summarizes the fault detection performance for the GPS constellation under different parameter combinations.

FIGURE 13

Fault detection results for the Earth case with zero or one faulty satellites. (a) [Earth] Probability of false alarm when there are no faults. (b) [Earth] 1 fault, α: 0.001, fault ratio: 1. (c) [Earth] 1 fault, α: 0.001, fault ratio: 0.2. (d) [Earth] 1 fault, a: 0.013, fault ratio: 0.2. (e) [Earth] P_fa vs P_md for different fault magnitudes (from left to right: 2 m, 4 m, 8 m, 16 m) and input P_fa (varied inside each subplot).

When there is no fault, the P_fa of all methods is below the user-defined P_fa (α = 0.001-0.01). When the fault ratio is 1.0, the data-snooping method shows the lowest P_md and P_fa, followed by the EDM approach and ephemeris approach.

a) No Fault

When no satellites are faulty, the false alarm probability P_fa of all methods remains below the user-specified threshold α = 0.001 – 0.1. Among the three approaches, the data-snooping method and the EDM-based method achieve the lowest P_fa, indicating strong robustness in fault-free conditions.

b) Single Faulty Satellite with Fault Ratio = 1

When the fault ratio is 1, meaning that all ISRs involving the faulty satellite are corrupted, the data-snooping method yields the best performance. Owing to the optimality of its statistical test under this assumption, it achieves nearly perfect fault detection (P_fa = P_md = 0) for fault magnitudes as small as 2 m, as shown in Figure 13(e). The P_fa – P_md trade-off curve of the EDM test lies between those of data snooping and the ephemeris comparison method, as illustrated in the same figure.

The performance of the ephemeris comparison method improves as the ephem-eris accuracy increases, as shown in Figure 13(b). For α = 0.001 and ephemeris accuracy σ_m = 1m, perfect fault detection (P_fa = P_md = 0) is achieved when the fault magnitude exceeds 5 m. In contrast, for a lower ephemeris accuracy (α_m = 2 m), a larger fault magnitude (approximately 8 m) is required to drive both P_fa and P_md to zero.

c) Single Faulty Satellite with Fault Ratio = 0.2

When the fault ratio is reduced to 0.2, corresponding to a scenario in which only a subset of links is biased (e.g., owing to link scheduling), the assumptions underlying the data-snooping test are violated. Consequently, data snooping is no longer optimal. Nevertheless, for small fault magnitudes, the data-snooping method still provides the most favorable P_fa – P_md trade-off among the three methods.

For larger fault magnitudes (e.g., ≥ 15 m), the ephemeris comparison method outperforms the other two approaches. This behavior arises because, when corrupted edges are sparse, both the EDM and data-snooping methods exhibit a higher probability of misattributing the fault to a non-faulty satellite. As in the previous cases, the EDM method consistently demonstrates intermediate performance between ephemeris comparison and data snooping.

5.2.2 Case 2: Lunar Constellation

Figure 14 summarizes the fault detection performance for the lunar constellation under different parameter combinations. Because the ephemeris augmentation strategy described in Section 4.3.3 is employed for the EDM approach, its detection performance also becomes explicitly dependent on the assumed ephemeris accuracy.

FIGURE 14

Fault detection results for the lunar case with zero or one faulty satellites. (a) [Moon] Probability of false alarm when there are no faults. (b) [Moon] 1 fault, α : 0.001, fault ratio: 1. (c) [Moon] 1 fault, α : 0.001, fault ratio: 0.2. (d) [Moon] 1 fault, α : 0.013, fault ratio: 0.2. (e) [Moon] P_fa vs P_md for different fault magnitudes (from left to right: 2 m, 4 m, 8 m, 16 m) and input Pfa (varied inside each subplot).

When there is no fault, the P_fa of all methods is below the user-defined Pf_a (αrid="R06"= 0.001-0.1). For a fault ratio of 1.0, the data-snooping method shows the lowest P_md and Pf_a, followed by the EDM approach and then the ephemeris approach. When the fault ratio is 0.2, the EDM and ephemeris approaches provide better P_fa and P_md values for large fault magnitudes.

a) No Fault

Consistent with the GPS case, when no satellites are faulty, the false alarm probability P_fa of all methods remains below the user-specified threshold α = 0.001 – 0.1. Among the three approaches, the data-snooping method and the EDM-based method achieve the lowest Pf_a, indicating robust behavior under fault-free lunar operating conditions.

b) Single Faulty Satellite with Fault Ratio = 1

When the fault ratio is 1.0, the data-snooping method again yields the best overall performance, owing to the optimality of its statistical test under the assumption that all links associated with the faulty satellite are corrupted. However, because of the reduced number and sparsity of ISLs in the lunar constellation compared with the GPS case, a larger fault magnitude (approximately ≥ 10 m) is required for data snooping to achieve nearly perfect fault detection (P_fa ≈ P_md ≈ 0), as shown in Figures 14(b) and 14(e).

As in the terrestrial case, the P_fa – P_md trade-off curve of the EDM approach lies between those of data snooping and the ephemeris comparison method. The performance of the ephemeris comparison method improves with increasing ephemeris accuracy, although its sensitivity to ephemeris errors is more pronounced in the lunar case owing to weaker geometric redundancy.

c) Single Faulty Satellite with Fault Ratio = 0.2

When the fault ratio is reduced to 0.2, corresponding to scenarios in which only a subset of ISLs is corrupted, the assumptions underlying the data-snooping test are violated. As a result, the performance of data snooping degrades significantly, with an elevated P_md observed even for moderate fault magnitudes.

In contrast, both the EDM and ephemeris comparison methods exhibit improved robustness in this sparse-fault regime. For larger fault magnitudes (e.g., ≥ 8 – 10 m), these two approaches achieve lower P_md while maintaining P_fa near or below the prescribed threshold. The EDM method consistently demonstrates intermediate performance, outperforming data snooping for sparse faults while remaining less sensitive to ephemeris inaccuracies than the ephemeris comparison method when high-accuracy ephemerides are not available.

6 CONCLUSION

This paper presented a novel framework for detecting satellite clock phase jumps using ISR measurements, grounded in rigid graph theory and EDM analysis. The proposed method models the constellation as a weighted graph and evaluates the realizability of vertex-redundantly rigid subgraphs to identify clock jumps that manifest as structured biases in the ISR.

We established the necessary graph-topological conditions for fault detectability and derived analytical expressions for the distribution of singular values of the GCEDM under both nominal and faulty conditions. These results enabled the development of statistically rigorous hypothesis tests with controllable false alarm probabilities and analytically characterized detection performance. To address sparse measurement geometries, we introduced an ephemeris-augmented formulation that increases fault detectability when inter-satellite connectivity is limited, at the cost of increased sensitivity to ephemeris uncertainty. Simulation studies using both a terrestrial GPS constellation and a notional lunar constellation demonstrated that the proposed method remains effective across a wide range of fault magnitudes, link topologies, and operational scenarios.

Overall, the proposed approach offers a geometry-driven perspective on rangebased fault detection through the lens of graph rigidity, complementing conventional residual-based methods. Future work will investigate extensions that incorporate modifications to data-snooping methods, multi-fault detection, and integration with link scheduling to maximize fault observability and detection performance.

HOW TO CITE THIS ARTICLE:

Iiyama, K., Neamati, D., & Gao, G. (2026). Satellite autonomous clock fault monitoring with inter-satellite ranges using euclidean distance matrices. NAVIGATION, 73. https://doi.org/10.33012/navi.764

A | COMPUTING THE TEST THRESHOLD USING IMHOF’S METHOD

In this section, we describe the procedure for obtaining the critical value for the hypothesis test in Section 3.1. The correlation between two residuals that share the same satellites can be described as follows (Xu et al., 2011):

$ρ (g_{i, j}, g_{i, k}) = \frac{cov (g_{i, j}, g_{i, k})}{\sqrt{cov (g_{i, j}, g_{i, j}) \cdot cov (g_{i, k}, g_{i, k})}} = \frac{β σ_{r}^{2}}{2 σ_{r}^{2} + σ_{m}^{2}}$ 64

Here, β is the cosine of the angle between two line-of-sight vectors e_t,j, e_i,k, that satisfies the following:

$\cos (ϕ_{max}) < \cos (β) < \cos (0) = 1$ 65

The covariance matrix Σ of the l_i normalized residuals ${\tilde{g}}_{i, j} = g_{i, j} / \sqrt{2 σ_{r}^{2} + σ_{m}^{2}}$ obtained at satellite i can be described as follows:

$Σ_{i, j} = {\begin{cases} 1 & i = j \\ \frac{β σ_{r}^{2}}{2 σ_{r}^{2} + σ_{m}^{2}} & i \neq j \end{cases}$ 66

Based on Imhof’s method, the cumulative distribution function value of the test statistic $T_{i} = \sum_{j = 1}^{l_{i}} {\tilde{g}}_{i j}^{2}$ is given as follows (Imhof, 1961):

$F_{Σ} (q) = P (T_{i} \leq q) = \frac{1}{2} - \frac{1}{π} \int_{0}^{\infty} ℑ (e^{- i t q} \cdot ϕ (t)) \frac{d t}{t}$ 67

Here, ℑ(⋅) denotes the imaginary part, and the characteristic function ϕ (t) is defined as follows:

$ϕ (t) = - \prod_{k = 1}^{l_{i}} \sqrt{1 - 2 i t λ_{k}}$ 68

where 𝜆_k (k = 1,..., l_i) are the eigenvalues of the covariance matrix Σ. However, it is computationally expensive to evaluate Equation (67) for every satellite at every time step. Therefore, Xu et al. (2011) suggested using a conservative bound for the cross-terms of Σ_i,j, as follows:

${\tilde{Σ}}_{i, j} = {\begin{cases} 1 & i = j \\ \frac{σ_{r}^{2}}{2 σ_{r}^{2} + σ_{m}^{2}} & i \neq j \end{cases}$ 69

Then, we can pre-compute $F_{\tilde{Σ}} (q)$ for l_t = 1,...,N_s and a pre-determined false alarm rate α = 1 – q and store them in a table. The test procedure using this conservative bound is as follows:

$Accept H_{0} if {max}_{j \in {1, \dots, n}} T_{i} < F_{\tilde{Σ}} (1 - α)$ 70

otherwise:

$Accept H_{k} where k = {argmax}_{j \in {1, \dots, n}} \frac{T_{j}}{\sqrt{l_{j}}}$ 71

B | PROOF OF PROPOSITION 4.1

Proof. Let r_uv represent the observed inter-note distances. Consider a k-vertex-redundantly rigid graph G with weight W (e_uv ) = r_uv that is d-embeddable. Let q: V ∈ ℝ³ be a generic realization (i.e., it does not lie on a plane or line, unless forced by the ranges) of the points with true inter-node distances ${\bar{r}}_{u v}$ given by the following:

${\bar{r}}_{u v} = ‖ q_{u} - q_{v} ‖$ 72

and p: V ∈ ℝ³ be the generic realization of the points with the observed inter-note distances r_uv:

$r_{u v} = ‖ p_{u} - p_{v} ‖ = ‖ q_{u} - q_{v} ‖ + f_{u v}$ 73

where f_uv ≠ 0 if u or v is the faulty satellite and f_uv = 0 otherwise.

We assume the following statement:

$G contains m (1 \leq m \leq k - 1) faulty nodes$ 74

Without a loss of generality, we can assume that the indices of the faulty nodes are 1,..., m for m ≤ k-1. Let G'= G\{1,..., m} be a graph without the m ≤ k-1 faulty nodes. Because G is k-vertex-redundantly rigid, G' is rigid in ℝ^d. In addition, because G' contains no fault nodes, we have the following:

${\bar{r}}_{u v} = r_{u v}, u, v \neq {1, \dots, m}$ 75

Therefore, we may choose coordinates so that we have the following:

$p_{j} = q_{j} j \neq {1, \dots, m}$ 76

For each fault node $i \in {1, \dots, m}$ , its neighbor $j \in N (i), N (i) = {k ∣ (i, k) \in G^{'}}$ satisfies the following:

$‖ p_{u} - p_{v} ‖ = ‖ q_{i} - q_{j} ‖ + f_{i j} (f_{i j} \neq 0)$ 77

Hence, p_i must lie in the following intersection:

$⋂_{j \in N (i)} {x \in R^{3} : ‖ x - q_{j} ‖ = ‖ q_{i} - q_{j} ‖ + f_{i j}}$ 78

Geometrically, this is the common intersection of ||𝒩(i)|| spheres in ℝ³. It is known that every vertex in a k-vertex-redundantly rigid graph must have a degree of at least d + k –1 (Jordán et al., 2022). Because m ≤ k -1, each fault node has at maximum k – 2 edges with other fault nodes (which are not included in G'); therefore, we have ||𝒩(i)|| ≥ (d + k – l)–(k – 2 ) = 3 – 1 + 2 = 4.

Let {j₁,..., j₄} ∈ 𝒩(i), with the following:

$r = {[‖ q_{i} - q_{j_{1}} ‖ + f_{i j_{1}} \dots ‖ q_{i} - q_{j_{4}} ‖ + f_{i j_{4}}]}^{⊤} \in R^{4}$ 79

We define a smooth map:

$F : R^{3} \to R^{4} F (x) = [\begin{array}{l} ‖ x - q_{i_{1}} ‖ & \dots & ‖ x - q_{i_{4}} ‖ \end{array}]$ 80

The spheres intersect if and only if the following holds:

$r \in F (R^{3})$ 81

Because the Jacobian of F can have a rank of at most 3, its image F(ℝ³) ⊂ ℝ⁴ is at most a 3D immersed submanifold of ℝ⁴. Any smooth embedding of a k-dimensional manifold into ℝⁿ has an n-dimensional Lebesgue measure of 0 whenever k <n (Guillemin & Pollack, 1974). Since the probability law of the random vector r is assumed to be absolutely continuous with respect to the 4-dimensional Lebesgue measure, the probability that Equation (81) holds is 0. This contradicts the statement that the graph G is d-embeddable with the observed ranges, at probability 1. Because we derived a falsehood, the assumption of Equation (74) that the graph contains 1 ≤ m ≤ k – 1 fault nodes is false, at probability 1. Therefore, if the graph is k-vertex rigid, is d-embeddable, and contains at most k – 1 faulty nodes, the graph does not contain any fault nodes, at probability 1.

C | PROOF OF PROPOSITION 4.3

Proof. As discussed by Dokmanić et al. (2015), the noiseless and faultess EDM, $D_{i j}^{n, d, 0} = {‖ x_{i} - x_{j} ‖}^{2}$ , based on a collection of points $X \in R^{d \times n}$ , can be constructed from Equation (82):

$D^{n, d, 0} = 1 diag {(X^{⊤} X)}^{⊤} - 2 X^{⊤} X + diag (X^{⊤} X) 1^{⊤}$ 82

where $diag (X^{⊤} X) \in R^{n}$ is the vector formed from the diagonal entries of $X^{⊤} X$ . The first and last matrices of Equation (82) are rank-1 matrices by construction. These rank-1 matrices each contribute to the matrix rank, except for degenerate cases, where the points are equally spread apart ( $X 1 = 0$ ). The middle term includes the Gram matrix $X^{⊤} X$ , which is rank d except for degenerate cases where the points lie on a lower-dimensional hyperplane (i.e., $rank (X^{⊤} X) = 2$ if the points lie on a plane). Thus, in line with the work of Dokmanić et al. (2015), $rank (D^{n, d, 0}) \leq d + 2$ , where the condition holds with equality outside of the aforementioned degenerate cases.

During a fault, we encounter cross-terms with the bias following Equation (83):

$\begin{array}{l} D_{i j}^{n, d, m} & = & {(‖ X_{i} - X_{j} ‖ + f_{i j})}^{2} = {‖ X_{i} - X_{j} ‖}^{2} + 2 ‖ X_{i} - X_{j} ‖ f_{i j} + f_{i j}^{2} \\ = & {‖ X_{i} - X_{j} ‖}^{2} + f_{i j} (2 ‖ X_{i} - X_{j} ‖ + f_{i j}) \\ = & {‖ X_{i} - X_{j} ‖}^{2} + f_{i j} (s_{i j} + f_{i j}) \end{array}$ 83

where s_i,j is two times the element-wise square root of the corresponding entry in the EDM. As a matrix, we can write $S_{i j}^{n, d} = s_{i j}$ We can write the bias term $F_{i j}^{n, m} = b_{i j}$ as a matrix as well. Both S^n,d and F^n,m are real, symmetric matrices. Then, we arrive at Equation (84):

$D^{n, d, m} = D^{n, d, 0} + F^{n, m} \circ (S^{n, d} + F^{n, m})$ 84

where ∘ is the element-wise or Hadamard product. For the sake of illustration, consider a case with three faults of equal magnitude and opposite signs realized at satellites k = {2,3,4} for n >7, without a loss of generality:

$F^{n, {2, 3, 4}} = \bar{f} [\begin{array}{c} 0 & 1 & 1 & 1 & 0 & \dots & 0 \\ 1 & 0 & 2 & 2 & 1 & \dots & 1 \\ 1 & 2 & 0 & 2 & 1 & \dots & 1 \\ 1 & 2 & 2 & 0 & 1 & \dots & 1 \\ 0 & 1 & 1 & 1 & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 1 & 1 & 1 & 0 & \dots & 0 \end{array}] = \bar{f} [\begin{array}{c} 0 & 1 & 1 & 1 \\ 1 & 0 & 2 & 2 \\ 1 & 2 & 0 & 2 \\ 1 & 2 & 2 & 0 \\ 0 & 1 & 1 & 1 \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 1 & 1 & 1 \end{array}] [\begin{array}{c} 1 & 0 & 0 & 0 & 1 & \dots & 1 \\ 0 & 1 & 0 & 0 & 0 & \dots & 0 \\ 0 & 0 & 1 & 0 & 0 & \dots & 0 \\ 0 & 0 & 0 & 1 & 0 & \dots & 0 \end{array}]$ 85

In the rank-decomposition of $F^{n, m}$ , the first column in the matrix has a one in each fault entry and a zero otherwise. This column spans all of the $n - m$ fault-free columns. The remaining m columns correspond to the columns in $F^{n, m}$ with faults and have a zero at the fault entry. In the above example, the entry 2 is illustrative. In general, the entry will depend on the fault magnitudes, which may be different. With probability zero, the faults across satellites exactly cancel out and lead to degeneracy in the fault detection, as discussed below. $F^{n, m}$ achieves full rank when the number of faults reaches $m = n - 1$ . In general, $rank (F^{n, m}) = min (1 + m, n)$ if $m > 0$ and $rank (F^{n, 0}) = 0$ . By a similar argument, $rank (F^{n, m} \circ F^{n, m}) = rank (F^{n, m})$ because we will have one column indicating the fault entries and m columns associated with the columns of $F^{n, m} \circ F^{n, m}$ with faults, which are the same columns as $F^{n, m}$ . This overall rank relationship often holds even if the faults are of different non-zero magnitudes. However, if the faults cancel out, the rank will drop. For example, in the above case, if $b_{2} = 2 \bar{f}$ and $b_{3} = b_{4} = - \bar{f}$ , then $rank (F^{n, m}) = 3$ . Therefore, $rank (F^{n, m}) \leq min (1 + m, n)$ for arbitrary fault sizes. Nevertheless, the term $rank (F^{n, m} \circ F^{n, m})$ can still be rank $m + 1$ even if $rank (F^{n, m}) < m + 1$ , as the squaring can remove the linear cancellation.

However, the span of the columns of $F^{n, m} \circ S^{n, d}$ will generally span the columns of $F^{n, m} \circ F^{n, m}$ . Again, for the sake of illustration, we extend the above example with faults at satellites $k = {2, 3, 4}$ for $n > 7$ , without a loss of generality:

$\begin{aligned} F^{n, {2, 3, 4} \circ} \circ S^{n, d} & = & \bar{f} [\begin{array}{c} 0 & s_{12} & s_{13} & s_{14} & 0 & \dots & 0 \\ s_{12} & 0 & 2 s_{23} & 2 s_{24} & s_{25} & \dots & s_{2 n} \\ s_{13} & 2 s_{23} & 0 & 2 s_{34} & s_{35} & \dots & s_{3 n} \\ s_{14} & 2 s_{24} & 2 s_{34} & 0 & s_{45} & \dots & s_{4 n} \\ 0 & s_{25} & s_{35} & s_{45} & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & s_{2 n} & s_{3 n} & s_{4 n} & 0 & \dots & 0 \end{array}] \\ = & \bar{f} [\begin{array}{cccccc} 0 & s_{12} & s_{13} & s_{14} & 0 & 0 \\ s_{12} & 0 & 2 s_{23} & 2 s_{24} & s_{25} & s_{26} \\ s_{13} & 2 s_{23} & 0 & 2 s_{34} & s_{35} & s_{36} \\ s_{14} & 2 s_{24} & 2 s_{34} & 0 & s_{45} & s_{46} \\ 0 & s_{25} & s_{35} & s_{45} & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & s_{2 n} & s_{3 n} & s_{4 n} & 0 & 0 \end{array}] [\begin{array}{ccccccccc} 1 & 0 & 0 & 0 & 0 & 0 & Δ_{17} & \dots & Δ_{1 n} \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & \dots & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & \dots & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & \dots & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & Δ_{57} & \dots & Δ_{5 n} \\ 0 & 0 & 0 & 0 & 0 & 1 & Δ_{67} & \dots & Δ_{6 n} \end{array}] \end{aligned}$ 86

Here, $Δ_{i j}$ are leftover terms from the row-reduction. In the rank-decomposition of $F^{n, m} \circ S^{n, d}, m$ columns correspond to the columns in $F^{n, m}$ with faults and have a zero at the fault entry, just as before. However, now, one column is generally not enough to span the fault-free columns because the entries in $S^{n, d}$ are not linearly related. For example, ( $s_{12}, s_{13}, s_{14}$ ) is generally not a linear scaling of ( $s_{25}, s_{35}, s_{45}$ ). These fault-free columns constitute an m-dimensional subspace, which will require m columns to fully span. $F^{n, m} \circ S^{n, d}$ achieves full rank when the number of faults reaches $m = ⌈ n / 2 ⌉$ , where $⌈ \cdot ⌉$ is the ceiling function. At that point, the subspace of fault-free columns is large enough that the fault columns will contribute to the span. Therefore, $rank (F^{n, m} \circ S^{n, d}) \leq min (2 m, n)$ , where we do not have equality when $S^{n, d}$ is degenerate (i.e., if the faulty satellites are above or below a plane of fault-free satellites).

The columns needed to span $F^{n, m} \circ S^{n, d}$ are the same as those needed for $F^{n, m} \circ F^{n, m}$ , with the same sparsity pattern. Therefore, we will have the following:

$rank (F^{n, m} \circ (S^{n, d} + F^{n, m})) = rank (F^{n, m} \circ S^{n, d}) \leq min (2 m, n)$ 87

Lastly, the subspace spanned with the fault contributions is distinct from the subspace spanned by the EDM. Therefore, we have the following:

$rank (D^{n, d, m}) \leq min (d + 2 + 2 m, n)$ 88

where the equality holds in non-degenerate cases.

D | PROOF OF PROPOSITION 4.4

Proof. First, by rank inequality, rank(AB) ≤ min(rank(A), rank(B)) for arbitrary matrices A and B. The rank of the geometric-centering matrix is rank(Jⁿ) = n – 1. Thus, rank(G^n,d,k) ≤ n – 1. However, this bound is too loose when there are few faults. Expanding the expression for geometric centering yields Equation (89):

$\begin{aligned} G^{n, d, m} = & - \frac{1}{2} J^{n} D^{n, d, m} J^{n} \\ = & - \frac{1}{2} J^{n} (D^{n, d, 0} + F^{n, m} \circ (S^{n, d} + F^{n, m})) J^{n} \\ = & - \frac{1}{2} J^{n} (1 diag {(X^{⊤} X)}^{⊤} - 2 X^{⊤} X + diag (X^{⊤} X) 1^{⊤}) J^{n} \\ - \frac{1}{2} J^{n} (F^{n, m} \circ (S^{n, d} + F^{n, m})) J^{n} \end{aligned}$ 89

First, geometric centering removes the two rank-1 matrices, as shown in Equations (90) and (91) (Dokmanić et al., 2015):

$- \frac{1}{2} J^{n} (1 diag {(X^{⊤} X)}^{⊤}) J^{n} = - \frac{1}{2} (1 diag {(X^{⊤} X)}^{⊤} - 1 diag {(X^{⊤} X)}^{⊤}) (I^{n} - \frac{1}{n} 11^{⊤}) = 0$ 90

$- \frac{1}{2} J^{n} (diag (X^{⊤} X) 1^{⊤}) J^{n} = - \frac{1}{2} (I^{n} - \frac{1}{n} 11^{⊤}) (diag (X^{⊤} X) 1^{⊤} - diag (X^{⊤} X) 1^{⊤}) = 0$ 91

For the $- 2 X^{⊤} X$ term, notice that the mean of the points is $μ = \frac{1}{n} X 1 \in R^{d}$ . Using this property yields Equation (92), which is a Gram matrix of the form $X_{c}^{⊤} X_{c}$ for the point matrix $X_{c}$ centered about the origin (Dokmanić et al., 2015). This shift will not change the rank, meaning that $rank (X^{⊤} X) = rank (X_{c}^{⊤} X_{c})$ :

$- \frac{1}{2} J^{n} (- 2 X^{⊤} X) J^{n} = {(X - μ 1^{⊤})}^{⊤} (X - μ 1^{⊤}) = X_{c}^{⊤} X_{c}$ 92

The remaining term is $\frac{1}{2} J^{n} (F^{n, m} \circ (S^{n, d} + F^{n, m})) J^{n}$ . By rank inequality, we have the following:

$\begin{array}{l} rank (J^{n} (F^{n, m} \circ (S^{n, d} + F^{n, m})) J^{n}) & \leq & min (rank (F^{n, m} \circ (S^{n, d} + F^{n, m})), rank (J^{n})) \\ = & min (2 m, n - 1) \end{array}$ 93

Hence, we are left with the following:

$rank (G^{n, d, k}) \leq min (d + 2 m, n - 1)$ 94

E | PROOF OF COROLLARY 4.4.1

Proof. In terms of rank, the noise acts as many small faults on each satellite, with the same sparsity structure as F^n,m in Proposition 4.3, almost surely with m = n, as there is zero probability that the randomly sampled noise is exactly zero or exactly cancels out. Using Proposition 4.4, we obtain the following:

$rank ({\tilde{G}}^{n, d, m}) = rank (G^{n, d, n}) \leq min (d + 2 n, n - 1) = n - 1$ 95

Therefore, succinctly, $rank ({\tilde{G}}^{n, d, m}) = n - 1$ , almost surely.

F | DERIVATION OF THE SCALE OF THE CHI-SQUARED DISTRIBUTION

Let D ∈ ℝ^5x5 be the distance matrix between the nodes where D_i,j = ||x_i – xj ||, and let E ∈ ℝ^5x5 be the perturbation matrix (noise matrix) where $E_{i j} = ω_{i j} \sim N (0, σ_{m}^{2})$ . The observed noisy distance matrix is $\tilde{D} = D + E$ The EDM is D ∘ D, where ∘ is the Hadamard (element-wise) product. The observed (noisy) GCEDM is as follows:

$\begin{array}{l} \tilde{G} & = & - \frac{1}{2} J^{5} (\tilde{D} \circ \tilde{D}) J^{5} \\ = & - \frac{1}{2} J^{5} [(D \circ D) + 2 (D \circ E) + (E \circ E)] J^{5} \\ \approx & G + Δ G (∵ E ≪ D) \end{array}$ 96

where ΔG = J⁵ (D o E)J⁵. The second-order error term (E o E) is approximated as zero in the third line. Here, G has rank 3, and $\tilde{G}$ has rank 4, as proved in Corollary 4.4.1. Therefore, rank (ΔG) ≥ 4 – 3 = 1, by rank subadditivity. We denote the SVD of G as follows:

$G = U Σ V^{⊤} = [U_{1} ∣ U_{0}] [\begin{array}{cc} Σ_{3} & 0_{3 \times 2} \\ 0_{2 \times 3} & 0_{2 \times 2} \end{array}] {[V_{1} ∣ V_{0}]}^{⊤} (U_{1}, V_{1} \in R^{5 \times 3} and U_{0}, V_{0} \in R^{5 \times 2})$ 97

From the SVD of the perturbed matrix $\tilde{G} = \tilde{U} \tilde{Σ} {\tilde{V}}^{⊤}$ , we obtain the following:

$\begin{array}{l} \tilde{Σ} & = & {\tilde{U}}^{⊤} \tilde{G} \tilde{V} \\ = & (U + Δ U)^{⊤} (G + Δ G) (V + Δ V) \\ = & Σ + U^{⊤} (Δ G) V + U^{⊤} G (Δ V) + Δ U^{⊤} GV + [(Δ U)^{⊤} Δ GV + U^{⊤} Δ G Δ V + (Δ U)^{⊤} Δ G Δ V] \\ \approx Σ + U^{⊤} (Δ G) V + Σ V^{⊤} (Δ V) + Δ U^{⊤} U Σ (∵ second- and third-order perturbation \approx 0) \\ = & [\begin{array}{l} Σ_{3} + U_{1}^{⊤} Δ G V_{1} + Δ_{11} & U_{1}^{⊤} Δ {GV}_{0} + Δ_{12} \\ U_{0}^{⊤} Δ G V_{1} + Δ_{21} & U_{0}^{⊤} Δ {GV}_{0} \end{array}] \\ (Δ_{11}, Δ_{12}, Δ_{21} : perturbation from Σ V^{⊤} (Δ V) + Δ U^{⊤} U Σ) \\ = & [\begin{array}{l} Σ_{3} + U_{1}^{⊤} Δ G V_{1} + Δ_{11} & 0_{3 \times 2} \\ 0_{2 \times 3} & U_{0}^{⊤} Δ G V_{0} \end{array}] \end{array}$ 98

$U_{0}^{⊤} Δ G V_{0}$ has rank 1 because $\tilde{G}$ has rank 4. Thus, the fourth singular value of $\tilde{G}$ satisfies the following:

$λ_{4}^{2} = {\tilde{Σ}}_{44}^{2} = {‖ U_{0}^{⊤} Δ G V_{0} ‖}_{F}^{2}$ 99

We observe that the expectation of the matrix elements inside the norm is zero:

$E [{(U_{0}^{⊤} Δ G V_{0})}_{i j}] = {(U_{0}^{⊤} J^{5} (D \circ E [E]) J^{5} V_{0})}_{i j} = 0 (∵ E [E_{i j}] = 0)$ 100

Moreover, the term ${(U_{0}^{⊤} Δ G V_{0})}_{i j}$ is a linear combination of the elements of E, where $E_{i j}$ follows a Gaussian distribution. Because a linear combination of Gaussian variables is also Gaussian, the term itself is normally distributed. Combined with the zero-mean property, the squared fourth singular value $λ_{4}^{2}$ follows a scaled central chi-squared distribution with one degree of freedom $(k)$ , where $λ_{4} / s^{2}$ follows a chi-squared distribution. Because the mean of a chi-squared distribution can be written as $k s^{2} (k = 1), s^{2}$ is given by the following:

$\begin{array}{l} s^{2} & = & E [λ_{4}^{2}] \\ = & E [{‖ U_{0}^{⊤} J^{5} (D \circ E) J^{5} V_{0} ‖}_{F}^{2}] \\ = & E [{‖ {\hat{U}}^{⊤} (D \circ E) \hat{V} ‖}_{F}^{2}] (\hat{U} = J^{5} U_{0}, \hat{V} = J^{5} V_{0} \in R^{5 \times 2}) \\ = & \sum_{a = 1}^{2} \sum_{b = 1}^{2} (E [{(\sum_{i = 1}^{5} \sum_{j = 1}^{5} {\hat{U}}_{i, a} D_{i j} E_{i j} {\hat{V}}_{j, b})}^{2}]) \end{array}$ 101

The distribution of the product of the two elements in the noise matrices can be represented as follows:

$E [E_{i j} E_{k l}] = {\begin{cases} σ_{m}^{2} & if (i, j) = (k, l) or (i, j) = (l, k) \\ 0 & otherwise \end{cases}$ 102

Therefore, the squared scale of 𝜆₄ is as follows:

$s^{2} = \sum_{a = 1}^{2} \sum_{b = 1}^{2} \sum_{i = 1}^{5} \sum_{j = 1}^{5} σ_{m}^{2} D_{i j}^{2} {({\hat{U}}_{i, a} {\hat{V}}_{j, b} + {\hat{U}}_{j, a} {\hat{V}}_{i, b})}^{2}$ 103

In addition, by comparison of the right bottom element of Equation (98), we obtain the following:

${\tilde{U}}_{0}^{⊤} (G + Δ G) {\tilde{V}}_{0} = U_{0}^{⊤} Δ G V_{0}$ 104

where ${\tilde{U}}_{0}, {\tilde{V}}_{0}$ are the fourth and fifth columns of $\tilde{U}, \tilde{V}$ . Therefore, by plugging Equation (99) into Equation (104), we obtain the same scale as Equation (103) by replacing D and $U_{0}, V_{0}$ with the perturbed distance matrix $\tilde{D}$ and singular vectors ${\tilde{U}}_{0}, {\tilde{V}}_{0}$ . Therefore, we obtain the following:

$s^{2} = \sum_{a = 1}^{2} \sum_{b = 1}^{2} \sum_{i = 1}^{5} \sum_{j = 1}^{5} {(σ_{m} {\tilde{D}}_{i j})}^{2} \cdot {({\hat{U}}_{i, a} {\hat{V}}_{j, b} + {\hat{U}}_{j, a} {\hat{V}}_{i, b})}^{2} (\hat{U} = J^{5} {\tilde{U}}_{0}, \hat{V} = J^{5} {\tilde{V}}_{0} \in R^{5 \times 2})$ 105

G | DERIVATION OF THE NON-CENTRALITY PARAMETER WHEN A FAULT EXISTS

Let the faulty satellite be denoted as S_f, with all of the ranges connected to the faulty satellite biased by $\bar{f}$ . Using the same notation as in Appendix F, the perturbed distance matrix ${\tilde{D}}_{F}$ is as follows:

${\tilde{D}}_{f} = D + F + E$ 106

where:

$F_{i j} = {\begin{cases} f_{i j} i = s_{f} or j = s_{f} \\ 0 otherwise \end{cases}$ 107

Therefore, the perturbed GCEDM is as follows:

${\tilde{G}}^{f} = - \frac{1}{2} J^{5} ({\tilde{D}}_{f} \circ {\tilde{D}}_{f}) J^{5} \approx G + Δ G^{F} + Δ G$ 108

where:

$Δ G^{F} = J^{5} (F \circ D) J^{5}$ 109

From Collloary 4.4.1, ${\tilde{G}}^{f}$ has rank 4. Similar to the discussion in Appendix F, the fourth singular value of ${\tilde{G}}^{f}$ satisfies the following:

$λ_{4}^{2} = \sum_{44}^{2} = ‖ U_{0}^{T} Δ G V_{0} + U_{0}^{T} Δ G^{f} V_{0} ‖_{F}^{2}$ 110

We analyze the expectation of the matrix term inside the norm:

$\begin{array}{l} E [U_{0}^{⊤} (Δ G + Δ G^{f}) V_{0}] & = & E [U_{0}^{⊤} Δ G V_{0}] + U_{0}^{⊤} Δ G^{f} V_{0} \\ = & 0 + U_{0}^{⊤} J^{5} (D \circ F) J^{5} V_{0} \\ \neq & 0 (∵ D \circ F \neq 0) \end{array}$ 111

The underlying variable is the sum of a zero-mean Gaussian term ( $U_{0}^{⊤} Δ G V_{0}$ ) and a non-zero constant bias term ( $U_{0}^{⊤} Δ G^{f} V_{0}$ ). Because the underlying variable is Gaussian with a non-zero mean, its squared Frobenius norm $λ_{4}^{2}$ follows a non-central chi-squared distribution.

The mean of this non-central chi-squared distribution is as follows:

$\begin{array}{l} s^{2} (k + λ) & = & E [λ_{4}^{2}] = E [{‖ U_{0}^{T} Δ G V_{0} ‖}_{F}^{2}] + E [\sum_{i = 1}^{2} \sum_{j = 1}^{2} {(U_{0}^{T} Δ G V_{0})}_{i j} {(U_{0}^{T} Δ G^{F} V_{0})}_{i j}] + {‖ U_{0}^{T} Δ G^{F} V_{0} ‖}_{F}^{2} \\ = & s^{2} + {‖ U_{0}^{T} Δ G^{F} V_{0} ‖}_{F}^{2} (∵ Eq. (101), E [{(U_{0}^{T} Δ G V_{0})}_{i j}] = 0) \\ = & s^{2} (1 + \frac{1}{s^{2}} {‖ U_{0}^{⊤} J (D \circ F) J V_{0} ‖}_{F}^{2}) \end{array}$ 112

where s² is the squared scale given in Equation (105).

Therefore, the non-centrality parameter is given by the following:

$λ_{s_{f}} = \frac{1}{s^{2}} {‖ U_{0}^{⊤} J^{5} (D \circ F) J^{5} V_{0} ‖}_{F}^{2}$ 113

H | GENERALIZATION OF THE TEST STATISTIC DISTRIBUTION FOR N ≥ 5

We generalize the derivation from Appendix F to an arbitrary number of satellites n ≥ 5. We define the generalized test statistic γ_test as the total energy of the noise subspace singular values:

$γ_{test} = \sum_{k = 4}^{n - 1} {\tilde{λ}}_{k}^{2}$ 114

where ${\tilde{λ}}_{k}$ represents the singular values of the observed (noisy) GCEDM $\tilde{G}$ . This statistic is equivalent to the squared Frobenius norm of the matrix projected onto the active noise subspace.

Because the true noiseless matrix G is unknown, we approximate the active noise subspace using the singular vectors of the measured matrix $\tilde{G}$ . Let the SVD of the observed matrix be $\tilde{G} = \tilde{U} \tilde{Σ} {\tilde{U}}^{⊤}$ . We partition the singular vectors $\tilde{U}$ as follows:

$\tilde{U} = [\begin{matrix} \underset{{\tilde{U}}_{s}}{\underset{︸}{\begin{matrix} {\tilde{u}}_{1} & {\tilde{u}}_{2} & {\tilde{u}}_{3} \end{matrix}}} & \underset{{\tilde{U}}_{noise}}{\underset{︸}{\begin{matrix} {\tilde{u}}_{4} & \dots & {\tilde{u}}_{n - 1} \end{matrix}}} & \underset{{\tilde{u}}_{1}}{\underset{︸}{{\tilde{u}}_{n}}} \end{matrix}]$ 115

Here, ${\tilde{U}}_{S}$ spans the signal subspace. The vectors ${\tilde{u}}_{4}, \dots, {\tilde{u}}_{n - 1}$ span the active noise subspace. The final vector ${\tilde{u}}_{n}$ corresponds to the structural null direction parallel to 1. Owing to geometric centering ( $\tilde{G} 1 = 0$ ), the singular value associated with ${\tilde{u}}_{n}$ remains strictly zero.

We approximate the projection of the noise perturbation $Δ G$ using this observed basis ${\tilde{U}}_{noise}$ . Let $M_{noise} \approx {\tilde{U}}_{noise}^{⊤} Δ G {\tilde{U}}_{noise}$ . Because $Δ G$ is linear in the Gaussian measurement noise $E, M_{noise}$ is a symmetric matrix of zero-mean Gaussian variables. The sum of squares of its independent elements follows a chi-squared distribution with k degrees of freedom:

$k = \frac{(n - 4) (n - 3)}{2}$ 116

(e.g., k = 1 for n = 5). The scaled test statistic follows $\frac{γ_{test}}{s^{2}} \sim χ^{2} (k, 0)$ , where the squared scale parameter s² is computed using the observed singular vectors ${\tilde{u}}_{a}$ (columns of ${\tilde{U}}_{noise}$ ) and the observed distances ${\tilde{D}}_{i j}$ :

$s^{2} \approx \frac{1}{k} \sum_{a = 4}^{n - 1} \sum_{b = 4}^{n - 1} (\sum_{i = 1}^{n} \sum_{j = 1}^{n} {(σ_{m} {\tilde{D}}_{i j})}^{2} {({\tilde{u}}_{i, a} {\tilde{u}}_{j, b} + {\tilde{u}}_{j, a} {\tilde{u}}_{i, b})}^{2})$ 117

This approximation holds because the perturbation of the singular vectors corresponding to the noise subspace is small (Equation (98)), allowing ${\tilde{U}}_{noise}$ (obtained from $\tilde{G}$ ) to serve as a valid approximation for the true noise subspace ${\hat{U}}_{N}$ (of G ) in the variance computation. Note that for $n = 5$ , Equation (117) reduces to the result derived in Equation (105).

When a fault occurs, the GCEDM includes a deterministic bias term $Δ G^{f} = J^{n} (D \circ F) J^{n}$ . The projection onto the noise subspace becomes $M_{total} = M_{noise} + M_{bias}$ , where $M_{bias} = {\hat{U}}_{N}^{⊤} Δ G^{f} {\hat{U}}_{N}$ . The test statistic $γ_{test} \approx {‖ M_{total} ‖}_{F}^{2}$ then follows a non-central chi-squared distribution:

$\frac{γ_{test}}{s^{2}} \sim χ^{2} (k, λ_{n c})$ 118

The non-centrality parameter 𝜆_nc can be computed as follows:

$λ_{n c} = \frac{1}{s^{2}} {‖ M_{b i a s} ‖}_{F}^{2} = \frac{1}{s^{2}} \sum_{a = 1}^{n - 4} \sum_{b = 1}^{n - 4} {({\hat{u}}_{a}^{⊤} Δ G^{f} {\hat{u}}_{b})}^{2}$ 119

For n = 5, this reduces to the result $λ_{s_{f}} = \frac{1}{s^{2}} {(u_{0}^{⊤} Δ G^{f} u_{0})}^{2}$ derived in Equation (113).

I PROOF OF PROPOSITION 4.5

Proof. We assume that graph G is redundantly rigid. A graph is fault-disprovable if the injection of a fault bias makes the set of edge constraints inconsistent (i.e., the graph cannot be embedded in Euclidean space).

Let v_f ∈ V be a satellite experiencing a clock fault. This fault introduces a bias b ≠ 0 to the measured ranges. The weights of the edges incident to v_f are given by the following:

$W (u, v_{f}) = {\begin{cases} ‖ X_{u} - X_{v_{f}} ‖ + b & if (u, v_{f}) \in E_{m} \\ ‖ X_{u} - X_{v_{f}} ‖ & if (u, v_{f}) \in E_{e} \end{cases}$ 120

where the measurement noise and ephemeris error are not considered in this proof.

Case 1: Sufficient Condition. Assume that v_f has at least one incident edge e_meas = (u, v_f) ∈ E_m. Because G is redundantly rigid, the distance between u and v_f is structurally constrained by the remaining edges in E \{e_meas}. The faulty weight W(e_meas) includes the bias b, creating a geometric contradiction with the constraints imposed by the rest of the graph. Thus, the graph becomes unrealizable, and the fault is detectable.

Case 2: Necessary Condition. Assume that v_f has no incident edges in E_m (i.e., deg_E (v_f) = 0). In this case, all edges connected to v_f are from the set E_e. The weights of these edges are derived solely from the ephemeris, which is independent of the instantaneous clock jump fault. Consequently, the geometry of v_f in the graph remains consistent with the ephemeris model, regardless of the actual clock fault. The graph remains realizable, and the fault cannot be detected by analyzing the rigidity or embeddability of G.

Therefore, for the graph to be fault-disprovable for any single satellite fault, it is necessary that every node v ∈ V have at least one incident edge in E_m.

Acknowledgments

We acknowledge Derek Knowles for reviewing this paper. This material is based upon work supported by the Nakajima Foundation and the National Science Foundation under grant no. DGE-1656518.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

REFERENCES

↵
Baarda, W. (1968). A testing procedure for use in geodetic networks. Netherlands Geodetic Commission, 2(5).
Google Scholar
↵
Cao, Y., Hu, X., Chen, J., Bian, L., Wang, W., Li, R., Wang, X., Rao, Y., Meng, X., & Wu, B. (2019). Initial analysis of the BDS satellite autonomous integrity monitoring capability. GPS Solutions, 23, 35. https://doi.org/10.1007/s10291-019-0829-z
Google Scholar
↵
Chen, L., Dai, Y., Gao, W., Cao, Y., Hu, Z., Ren, Q., Nie, X., Zheng, J., Shao, R., Pei, L., & Wang, L. (2022). Performance analysis of BDS-3 SAIM and enhancement research on autonomous satellite ephemeris monitoring. Remote Sensing, 14(15), 3543. https://doi.org/10.3390/rs14153543
Google Scholar
↵
Chiba, N., & Nishizeki, T. (1985). Arboricity and subgraph listing algorithms. SIAM Journal on Computing, 14(1), 210–223. https://doi.org/10.1137/0214017
Google Scholar
↵
de Jong, K., & Teunissen, P. J. G. (2014). Minimal detectable biases of GPS observations for a weighted ionosphere. Earth, Planets and Space, 52(10), 857–862. https://doi.org/10.1186/BF03352295
Google Scholar
↵
Dokmanić, I., Parhizkar, R., Ranieri, J., & Vetterli, M. (2015). Euclidean distance matrices: Essential theory, algorithms and applications. arXiv. https://doi.org/10.48550/arXiv.1502.07541
Google Scholar
↵
Guillemin, V., & Pollack, A. (1974). Differential topology (1st ed.). Prentice-Hall.
Google Scholar
↵
Hendrickson, B. (1992). Conditions for unique graph realizations. SIAM Journal on Computing, 21(1), 65–84. https://doi.org/10.1137/0221008
Google Scholar
↵
Iiyama, K., Neamati, D., & Gao, G. (2024, September). Autonomous constellation fault monitoring with inter-satellite links: A rigidity-based approach. In Proceedings of the 37th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2024), (pp. 2807–2824). https://doi.org/10.33012/2024.19935
Google Scholar
↵
Iiyama, K., Vila, G. C., & Gao, G. (2024, March). Contact plan optimization and distributed state estimation for delay tolerant satellite networks. In 2024 IEEE Aerospace Conference (pp. 1–13). https://doi.org/10.1109/AERO58975.2024.10521114
Google Scholar
↵
Imhof, J. (1961). Computing the distribution of quadratic forms in normal variables. Biometrika, 48(3–4), 419–426. https://doi.org/10.1093/biomet/48.3-4.419
Google Scholar CrossRef Web of Science
↵
Imparato, D., Teunissen, P. J., & Tiberius, C. C. J. M. (2018). Minimal detectable and identifiable biases for quality control. Survey Review, 51(367), 289–299. https://doi.org/10.1080/00396265.2018.1437947
Google Scholar
↵
Israel, D. J., Mauldin, K. D., Roberts, C. J., Mitchell, J. W., Pulkkinen, A. A., Cooper, L. V. D., Johnson, M. A., Christe, S. D., & Gramling, C. J. (2020, March). LunaNet: A flexible and extensible lunar exploration communications and navigation infrastructure. In 2020 IEEE Aerospace Conference (pp. 1–14). https://doi.org/10.1109/AERO47225.2020.9172509
Google Scholar
↵
Jian, L., Yang, Z., & Liu, Y. (2010, March). Beyond triangle inequality: Sifting noisy and outlier distance measurements for localization. In 2010 Proceedings IEEE INFOCOM (pp. 1–9). https://doi.org/10.1109/INFCOM.2010.5462019
Google Scholar
↵
Jordán, T., Poston, C., & Roach, R. (2022). Extremal families of redundantly rigid graphs in three dimensions. Discrete Applied Mathematics, 322, 448–464. https://doi.org/10.1016/j.dam.2022.03.006
Google Scholar
↵
Kelso, T. (2024). A new way to obtain GP data (aka TLEs). Retrieved April 28, 2025, from https://celestrak.org/NORAD/documentation/gp-data-formats.php
Google Scholar
↵
Knowles, D., & Gao, G. (2023). Euclidean distance matrix-based rapid fault detection and exclusion. NAVIGATION, 70(1). https://doi.org/10.33012/navi.555
Google Scholar
↵
Knowles, D., & Gao, G. (2024). Greedy detection and exclusion of multiple faults using euclidean distance matrices. arXiv. https://doi.org/10.48550/arXiv.2404.12617
Google Scholar
↵
Laman, G. (1970). On graphs and rigidity of plane skeletal structures. Journal of Engineering Mathematics, 4(4), 331–340. https://doi.org/10.1007/BF01534980
Google Scholar CrossRef Web of Science
↵
Lehmann, R., & Voß-Böhme, A. (2017). On the statistical power of Baarda’s outlier test and some alternative. Journal of Geodetic Science, 7(1), 68–78. https://doi.org/10.1515/jogs-2017-0008
Google Scholar
↵
Li, R., Gao, S., Qin, L., Wang, G., Yang, W., & Yu, J. X. (2020). Ordering heuristics for k-clique listing. Proceedings of the VLDB Endowment, 13, 2536–2548. https://doi.org/10.14778/3407790.3407843
Google Scholar
↵
Misra, P., Bayliss, E., Lafrey, R., Pratt, M., & Muchnik, R. (1993). Receiver autonomous integrity monitoring (RAIM) of GPS and GLONASS. NAVIGATION, 40(1), 87–104. https://doi.org/10.1002/j.2161-4296.1993.tb02296.x
Google Scholar
↵
Motevallian, S. A., Yu, C., & Anderson, B. D. (2015). On the robustness to multiple agent losses in 2D and 3D formations: Multiple agent losses in 2d and 3d formations. International Journal of Robust and Nonlinear Control, 25(11), 1654–1687. https://doi.org/10.1002/rnc.3167
Google Scholar
↵
National Aeronautics and Space Administration. (2022). Lunar communications relay and navigation systems (LCRNS) preliminary lunar relay services requirements document (SRD). https://esc.gsfc.nasa.gov/static-files/ESC-LCRNS-REQ-0090%20Rev_B%2012-05-2022%20DCN001.pdf
Google Scholar
↵
Rodríguez-Pérez, I., García-Serrano, C., Catalán Catalán, C., García, A. M., Tavella, P., Galleani, L., & Amarillo, F. (2011). Inter-satellite links for satellite autonomous integrity monitoring. Advances in Space Research, 47(2), 197–212. https://doi.org/10.1016/j.asr.2010.07.019
Google Scholar
↵
Ryden, G., & Volle, M. (2025). NASA lunar communications relay and navigation systems (LCRNS) reference constellation 3.1 (tech. rep.). NASA Goddard Space Flight Center. https://esc.gsfc.nasa.gov/static-files/LCRNS_Reference_Constellation_White_Paper_03_2025.pdf
Google Scholar
↵
van Diggelen, F. (2009). A-GPS: Assisted GPS, GNSS, and SBAS. Artech House.
Google Scholar
↵
Walter, T., Gunning, K., Phelts, R. E., & Blanch, J. (2018). Validation of the unfaulted error bounds for ARAIM. NAVIGATION, 65(1), 117–133. https://doi.org/10.1002/navi.214
Google Scholar
↵
Weiss, M., Shome, P., & Beard, R. (2010, November). On-board GPS clock monitoring for signal integrity. In Proceedings of the 42nd Annual Precise Time and Time Interval Systems and Applications Meeting (pp. 465–480). https://www.ion.org/publications/abstract.cfm?articleID=10752
Google Scholar
↵
Wolf, R. (2000, September). Onboard autonomous integrity monitoring using intersatellite links. In Proceedings of the 13th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS 2000) (pp. 1572–1581) https://www.ion.org/publications/abstract.cfm?articleID=1564
Google Scholar
↵
Xie, X., Geng, T., Ma, Z., Chen, L., & Liu, J. (2022). Estimation and analysis of BDS-3 satellite yaw attitude using inter-satellite link observations. GPS Solutions, 26, 106. https://doi.org/10.1007/s10291-022-01290-8
Google Scholar
↵
Xu, H., Wang, J., & Zhan, X. (2011). GNSS satellite autonomous integrity monitoring (SAIM) using inter-satellite measurements. Advances in Space Research, 47(7), 1116–1126. https://doi.org/10.1016/j.asr.2010.11.026
Google Scholar
↵
Zhang, C., Geng, T., Xie, X., Zhao, Q., Li, T., Li, Z., & Meng, Y. (2024). Analysis of BDS intersatellite link ranging performance. Advances in Space Research, 73(10), 4955–4966. https://doi.org/10.1016/j.asr.2024.02.023
Google Scholar

PDF

Previous Next

In this Issue

Print
Download PDF
Article Alerts

Alerts for this Article

User Name *

Password *

Sign In to Email Alerts with your Email Address
Email *
Email Article

Email This Article

Thank you for your interest in spreading the word on NAVIGATION: Journal of the Institute of Navigation.
NOTE: We only request your email address so that the person you are recommending the page to knows that you wanted them to see it, and that it is not junk mail. We do not capture any email address.
Your Email *

Your Name *

Send To *

Enter multiple addresses on separate lines or separate them with commas.

You are going to email the following
Satellite Autonomous Clock Fault Monitoring with InterSatellite Ranges Using Euclidean Distance Matrices
Message Subject
(Your Name) has sent you a message from NAVIGATION: Journal of the Institute of Navigation

Message Body
(Your Name) thought you would like to see the NAVIGATION: Journal of the Institute of Navigation web site.

Your Personal Message
Citation Tools
Citation Tools
Satellite Autonomous Clock Fault Monitoring with InterSatellite Ranges Using Euclidean Distance Matrices

Iiyama Keidai

Iiyama Keidai

Neamati Daniel,

Neamati Daniel,

and Gao Grace

Gao Grace

Institute of Navigation

May 2026,

73

navi.764;

DOI: 10.33012/navi.764

Citation Manager Formats

BibTeX
Bookends
EasyBib
EndNote (tagged)
EndNote 8 (xml)
Medlars
Mendeley
Papers
RefWorks Tagged
Ref Manager
RIS
Zotero

Share
Share This Article
Satellite Autonomous Clock Fault Monitoring with InterSatellite Ranges Using Euclidean Distance Matrices

Iiyama Keidai

Iiyama Keidai

Neamati Daniel,

Neamati Daniel,

and Gao Grace

Gao Grace

Institute of Navigation

May 2026,

73

navi.764;

DOI: 10.33012/navi.764

Share This Article:
Post
Like 0
Bookmark this Article

SAVE ITEM IN SELECTED FOLDER.

User Name *

Password *

Remember my user name & password.

Jump to section

No citing articles found.

Google Scholar

Show more Original Article

PDF Help

Keywords

fault detection, inter-satellite links, lunar positioning, navigation, timing (PNT), satellite constellation

Satellite Autonomous Clock Fault Monitoring with InterSatellite Ranges Using Euclidean Distance Matrices

Abstract

1 INTRODUCTION

2 MEASUREMENTS FROM INTER-SATELLITE LINKS

2.1 Dual One-Way ISR Measurement

2.2 Satellite Clock Phase Jumps

3 BASELINE METHODS: CLOCK FAULT DETECTION USING ISR RESIDUALS

3.1 Fault Detection Using the Expected Range from Ephemeris

3.2 Data Snooping with Baarda’s w-Test Statistic

4 PROPOSED CLOCK FAULT DETECTION ALGORITHM USING EUCLIDEAN DISTANCE MATRICES

4.1 Required Graph Topology for Fault Detection

4.1.1 Definitions

4.1.2 Embeddability of Graphs with Faulty Satellites

4.2 Properties of the Euclidean Distance Matrix

4.2.1 Geometric-Centered Euclidean Distance Matrix

4.2.2 Rank of EDMs with Faulty Satellites

4.2.3 Distribution of the Fourth Singular Value of the GCEDM

4.3 Fault Detection Algorithm Using the EDM

4.3.1 Clique Listing

4.3.2 Online Fault Detection

4.3.3 Fault-Detectable Graphs with Ephemeris Augmentation

5 SATELLITE FAULT DETECTION SIMULATION

5.1 Simulation Configuration

5.1.1 Orbits and Links

5.1.2 Parameters

5.1.3 Evaluation Metrics

5.2 Fault Detection Results

5.2.1 Case 1: GPS Constellation

5.2.2 Case 2: Lunar Constellation

6 CONCLUSION

HOW TO CITE THIS ARTICLE:

A | COMPUTING THE TEST THRESHOLD USING IMHOF’S METHOD

B | PROOF OF PROPOSITION 4.1

C | PROOF OF PROPOSITION 4.3

D | PROOF OF PROPOSITION 4.4

E | PROOF OF COROLLARY 4.4.1

F | DERIVATION OF THE SCALE OF THE CHI-SQUARED DISTRIBUTION

G | DERIVATION OF THE NON-CENTRALITY PARAMETER WHEN A FAULT EXISTS

H | GENERALIZATION OF THE TEST STATISTIC DISTRIBUTION FOR N ≥ 5

I PROOF OF PROPOSITION 4.5

Acknowledgments

REFERENCES

In this Issue

NAVIGATION: Journal of the Institute of Navigation

Alerts for this Article

Email This Article

Citation Tools

Satellite Autonomous Clock Fault Monitoring with InterSatellite Ranges Using Euclidean Distance Matrices

Citation Manager Formats

Share This Article

Satellite Autonomous Clock Fault Monitoring with InterSatellite Ranges Using Euclidean Distance Matrices

SAVE ITEM IN SELECTED FOLDER.

Jump to section

Related Articles

Cited By...

More in this TOC Section

Similar Articles

Keywords