## Abstract

In this work, a new time-sequential positioning and fault detection method is developed for dual-frequency, multi-constellation Advanced Receiver Autonomous Integrity Monitoring (ARAIM). Unlike conventional “snapshot” ARAIM, sequential ARAIM exploits changes in satellite geometry at the cost of slightly higher computation and memory loads. From the perspective of users on Earth, the motion of any given GNSS satellite is small over short time intervals. But the accumulated geometry variations of redundant satellites from multiple GNSS can be substantial. This paper quantifies performance benefits brought by satellite motion to ARAIM. It specifically addresses the following challenges: (a) defining raw GNSS code and carrier error models over time, (b) designing estimators and fault detectors exploiting geometric diversity for positioning, cycle ambiguity estimation, and integrity evaluation, and (c) formulating these algorithms in a computationally efficient implementation. Performance improvements provided by sequential ARAIM over snapshot ARAIM are evaluated by worldwide availability analysis for aircraft approach navigation.

## 1 INTRODUCTION

This paper describes the design, analysis, and evaluation of a new time-sequential positioning and fault detection method for Advanced Receiver Autonomous Integrity Monitoring (ARAIM) using the dual-frequency, multiconstellation Global Navigation Satellite System (GNSS). The new approach differs from prior work on “snapshot” (or instantaneous) ARAIM algorithms^{1-3} in that it also exploits satellite motion. This provides observability on constant measurement bias errors^{4} and is used in this work to estimate floating (real valued) carrier phase cycle ambiguities, thereby improving navigation accuracy and integrity. The new approach is formulated as a straightforward, computationally efficient augmentation of snapshot ARAIM.

With the modernization of GPS, the full deployment of GLONASS, and the emergence of Galileo and Beidou, a greatly increased number of redundant ranging signals become available, which has recently drawn a renewed interest in RAIM. RAIM exploits redundant GNSS measurements to achieve self-contained fault detection at the user receiver.^{5,6} In particular, RAIM can help relax requirements on ground-based integrity monitors. For example, researchers in the European Union and in the United States are investigating ARAIM for worldwide vertical guidance of aircraft.^{1-3}

One of the primary tasks in RAIM is to evaluate integrity risk or, equivalently, the protection levels (which are probabilistic bounds on positioning errors). Integrity risk is the probability of undetected faults causing unacceptably large positioning errors. Multiple research efforts have recently been conducted to design optimal estimators and detectors that minimize the integrity risk in ARAIM, while meeting specified continuity and accuracy criteria.^{7,8} These methods have been employed in the ARAIM Milestone Reports 2^{2} (MS2), and Milestone Report 3^{3} (MS3) to identify the circumstances under which dual-frequency GPS/Galileo could satisfy LPV-200 requirements globally. The localizer precision vertical (LPV) requirements are set to support vertical navigation during approach operations down to 200-ft altitude above the ground. The ARAIM MS2-3^{2,3} shows that worldwide coverage of LPV-200 is achievable using optimal snapshot ARAIM algorithms for a wide range of nominal measurement error and fault parameters.

However, the MS2-3^{2,3} also points out cases where LPV200 is not achievable using dual-frequency GPS and Galileo, for example, when the nominal constellations are even slightly depleted—eg, by a single satellite in each. In addition, the alert limit (ie, the limit on acceptable positioning errors that defines hazardous situations) is 35 m for LPV-200, which is much larger than, for example, the category II precision approach alert limit requirement of 10 m.^{9} Thus, given that the methods used in Blanch et al.^{7} and Joerger et al.^{8} reach the best achievable performance and unless additional satellites from GLONASS or BeiDou are incorporated, snapshot ARAIM algorithms cannot provide global service better than LPV-200.

In response, this work explores the potential benefit of new integrity monitoring methods that exploit satellite motion in dual-frequency, multi-constellation ARAIM.

Unlike snapshot ARAIM algorithms, where carrier-smoothed code (CSC) measurements from multiple satellites are combined at one instant in time, position estimates obtained, for example, from a batch estimator (or finite-interval estimator^{10}) are directly derived from time sequences of raw measurements. The batch is more computation and memory expensive, but it gives the means to exploit satellite motion over short time intervals. The time interval is limited by the minimum mission duration for an ARAIM-equipped aircraft. The shortest mission duration would be for an aircraft compelled to land just after takeoff, which would include taxiing, takeoff, go-around, and landing. A nominal batch duration of 10 minutes is therefore assumed, and periods of up to 20 minutes are considered in sensitivity analyses.

In the next section of this paper, previous research efforts are outlined, which specifically use geometry change for carrier phase–based positioning. These prior references have in common that they rely on large changes in geometry from a small number of ranging sources. Conversely, in this paper, we exploit the combination of small angular variations from many space vehicles (SV).

The third section of this paper addresses the key challenge of measurement error modeling over time. The ARAIM error models given in Working Group C^{1} provide values of the instantaneous CSC measurement error standard deviations due to satellite clock and orbit ephemeris, troposphere, multipath, and receiver noise. For each of these error sources, this paper builds upon the work of Joerger et al.^{11} to model the time correlation of errors affecting raw code and carrier measurements, which are the inputs to the CSC filter used in snapshot ARAIM.^{1-3} An analysis of 9 months of experimental data is carried out to validate the new temporal GPS satellite clock and orbit ephemeris error models.

These raw measurement error models are then processed in a batch estimator and detector. In this implementation, a batch estimator is preferred over a Kalman filter (KF) because KF-based integrity risk evaluation would require banks of KFs as in Brenner.^{12} The batch estimator runs sequentially using a sliding-window mechanism. A compact, computationally efficient formulation is introduced using raw carrier and CSC measurements taken at a few, infrequent sample times within the fixed batch interval. This algorithm exploits two fundamental estimation principles: (a) Raw code errors due to multipath and receiver noise are averaged out using CSC; (b) redundant satellite motion is captured using raw carrier at a few, temporally separated sample times. This second principle is not leveraged in snapshot ARAIM.

Fault detection is performed using a batch multiple hypothesis solution separation (MHSS) algorithm. The batch MHSS approach is similar to snapshot ARAIM, except that a time sequence of measurements is processed rather than a single set of CSC data. The proposed formulation is a straightforward augmentation of snapshot ARAIM that only requires a few past-time measurements.

The fifth section of this paper presents a step-by-step analysis of the new method. The section first shows that the raw measurement error models match the ARAIM assumptions on CSC errors. It then evaluates, in comparison with snapshot ARAIM, the integrity risk reduction provided by satellite motion using only three samples over a nominal 10-minute batch period, which is possible using the computationally efficient formulation.

Finally, a performance analysis is carried out for aircraft approach applications using ARAIM with dual-frequency GPS and Galileo satellite measurements. The results quantify the substantial reduction in integrity risk achieved using batch ARAIM as compared to the baseline snapshot ARAIM algorithm. Worldwide availability maps show that, unlike snapshot ARAIM, batch ARAIM meets LPV-200 requirements even with depleted constellations. Under conditions described in the paper, batch ARAIM can even satisfy much more stringent requirements, including a 10-m alert limit.

## 2 BACKGROUND ON GEOMETRY CHANGE

This work aims at exploiting changes in satellite geometry to obtain fast and accurate estimates of carrier phase cycle ambiguities.^{4} The Integrity Beacon Landing System (IBLS) devised in the early 1990s^{13} and experimentally demonstrated in Cohen et al.^{14} and Gutt et al.^{15} was an early implementation of this principle for aircraft precision approach and landing. GPS signal transmitters serving as pseudo-satellites (“pseudolites”) placed on the ground along the airplane’s trajectory provided additional ranging sources and a large geometry change as the receiver’s downward-looking antenna flew over the installation. By 2000, Rabinowitz et al.^{16} designed a receiver capable of tracking carrier phase measurements from GPS and from the low Earth orbiting (LEO) telecommunication constellation GlobalStar. Using GlobalStar satellites’ large range variations, precise cycle ambiguity resolution and positioning were achieved within 5 minutes. Another example of LEO satellite-augmented GPS was analyzed in 2010,^{11} where ranging signals from Iridium satellites were combined with GPS to provide worldwide carrier phase positioning and fault detection using RAIM.

In each of the above references, greatly improved positioning performance was achieved by exploiting the fast relative angular motion between user receiver and the few ranging sources (LEO SVs or pseudolites) in view during short (5-10 minutes long) mission durations. Conversely, in this work, we only consider measurements from medium earth orbiting (MEO) GNSS satellites, which are slowly moving from the perspective of a user near Earth’s surface. It is the multiplicity of ranging sources from several GNSS constellations that is exploited here to achieve a significant accumulated geometry change.

To reinforce this idea, Appendix A introduces a relative positioning problem inspired from the kinematic surveying problem given in Hwang.^{17} This qualitative exercise helps identify the driving mechanisms of carrier phase–based estimation. In Hwang,^{17} differential carrier phase measurements from a single SV were used to estimate a one-dimensional horizontal baseline. We slightly modify this problem by including measurements from multiple satellites to estimate a vertical displacement with respect to an initial position, as illustrated in Figure 1. Appendix A shows that the variance of the baseline estimate error is given by

1

where is the variance of the vertical position estimate error

is the variance of the carrier phase ranging measurement error (assumed constant and equal for all SVs)

^{i}*θ*_{0} is the zenith angle for SV *i* at some initial time 0

^{i}*θ*_{1} is the zenith angle for SV *i* at a later time 1

Equation (1) describes two fundamental mechanisms by which can be reduced. First, the larger the angular variation between ^{i}*θ*_{0} and ^{i}*θ*_{1} is, the smaller becomes.^^{17} Second, the larger the number of visible SVs *n* is, the smaller becomes. In multi-constellation ARAIM, even though ^{i}*θ* variations may be small over short time periods, the number *n* of contributing terms in the denominator of Equation (1) can be large enough to provide significant sigma_V^2 reduction.

Further results in Appendix A show that the motion of low-elevation satellites contributes to vertical positioning more than geometry changes of high-elevation SVs.

## 3 RAW MEASUREMENT ERROR MODEL DERIVATION

Two fundamental principles are exploited in batch ARAIM. The first one is code noise averaging using carrier-based smoothing. This is also achieved in snapshot ARAIM via Hatch filters, which are geometry-free filters that average code minus carrier over time. The second principle, which is not leveraged in snapshot ARAIM, is cycle ambiguity estimation using changes in geometry.

Both snapshot and batch estimation are represented in Figure 2. Let *n* be the number of satellites in view and *q* be the number of samples, collected from batch start time at epoch 1 to current time *q*. For snapshot estimation, raw code measurements ^{i}*ρ*_{Q} for SV *i* (as indicated by the left superscript) are smoothed using carrier data ^{i}*φ*_{Q} sampled at all times 1 to *q* (right subscript “*Q*,” where *Q* ≡ 1,…,*q*). The resulting carrier-smoothed code measurements at time *q* for *i* = 1,…,*n*, are then used, for example, in a weighted least-squares estimator^{1} or in a modified non-least-squares estimator,^{2,3,8} to determine the user’s position at one instant in time. In contrast, the batch estimator directly uses raw code and carrier signals from all SVs at all times to simultaneously estimate user position at time *q* and floating cycle ambiguities.

The first key challenge in this work is that the ARAIM measurement error models in Working Group C^{1-3} were only established for CSC code , using large data collection campaigns in Murphy et al.^{18,19} However, evaluating batch ARAIM performance requires raw measurement error models for ^{i}*ρ*_{Q} and ^{i}*φ*_{Q}. This work builds upon prior research on ranging error time correlation in Joerger et al^{11} to derive statistical models over time for raw data used in batch ARAIM.

The linearized ionosphere error–free code and carrier phase measurement equations for satellite *i* at time *k*, respectively, This comments is in relation to the Equation (2) throughout the paper: can we have the second line of equations to be to the right of the equal sign?... (these equations should really be a single line, but the columns are not large enough) are as follows:

2

where ^{i}**e**_{k} is the 3 × 1 line-of-sight vector in a local reference frame (eg, north-east-down or NED) for satellite *i* at time *k*

**x**_{k} is the 3 × 1 user position vector in NED

*τ*_{k} is the receiver clock offset

^{i}*ε*_{T,k} is the tropospheric error

^{i}*ε*_{E,k} is the SV clock and orbit ephemeris error

^{i}*b*_{ρ,k} is a nominal code bias mainly due to code signal deformation

^{i}*b*_{φ,k} is a nominal carrier bias, which appears in the paired overbounding process^{20,21}

^{i}*η* is the carrier phase cycle ambiguity (constant over time)

^{i}*ε*_{MP,ρ,k} and ^{i}*ε*_{MP,φ,k} respectively, are code and carrier errors due to multipath

^{i}*ε*_{RN,ρ,k} and respectively, are code and carrier receiver noise terms

The following subsections describe models for each individual source of error.

### 3.1 Tropospheric delay

We model the tropospheric delay for SV *i* at time *k* as follows:

3

where ^{i}*v*_{ZTD,k} is the zenith tropospheric delay

^{i}*θ*_{k} is the elevation angle for SV *i* at time *k*

^{i}*c*_{T,k} is the tropospheric zenith-to-slant mapping coefficient:

Equation (3) is the instantaneous error contribution assumed in snapshot ARAIM^{1-3}: ^{i}*v*_{ZTD,k} is zero-mean normally distributed with variance . We use the notation where *σ*_{ZTD} = 0.12 m.^{1} In this paper, ^{i}*v*_{ZTD,k} is a random time sequence (white noise) independent across SVs *i*, for *i* = 1,…,*n*. Unlike in previous work,^{22,23} we do not model the tropospheric error correlation over time or across satellites. This is a conservative assumption that will be refined in future work using tropospheric error data.

### 3.2 Satellite clock and orbit ephemeris error

The satellite clock and orbit ephemeris error equation for satellite *i* at time *k* is as follows:

4

where ^{i}*b*_{E} is an unknown, constant bias

^{i}*g*_{E} is an unknown, constant gradient (ramp)

*t*_{1} is the first time epoch of the batch interval

*t*_{k} represents any time *k* during the batch interval, such that 1 ≤ *k* ≤ *q*, where *q* is the last epoch of the batch interval

^{i}*ε*_{RES,k} residual errors not captured by the “bias plus ramp” model

The first term in Equation (4) captures the instantaneous uncertainty in ^{i}*ε*_{E,k}. It is therefore set according to ARAIM assumptions as , nominally, *σ*_{URA} = 1 m.^{2,3} It can easily be shown that a time-invariant bias ^{i}*b*_{E} affecting a sequence of raw code measurements input to a Hatch filter causes the same error ^{i}*b*_{E} in the resulting output CSC .

However, the ARAIM error model does not specify whether nominal SV clock and orbit errors can be assumed constant over any particular time interval. To account for temporal error variation, a second term in Equation (4) is added, which is a ramp over time with an unknown but constant gradient ^{i}*g*_{E}, accounting for linear variations from the initial value over (*t*_{k} − *t*_{1}). Joerger et al.^{11} give evidence that periodic variations of orbit errors are on the order of the MEO GNSS orbital period, which supports the use of a simple linear model over 10 minutes. It also cites Parkinson and Spilker,^{24} Warren and Raquet,^{25} and Pervan and Gratton^{26} to establish a distribution on the rate of change of orbit/clock errors: . In addition, Section 3.5.2 of the GPS Standard Positioning Service Performance Standard (SPS PS)^{27} describes the signal-in-space (SIS) user range rate error (URRE), which supports the model of linear orbit and clock variations. The GPS SPS PS^{27} does not specify an integrity performance standard yet.

To validate Equation (4), a preliminary experimental data analysis is carried out. Similar to prior work,^{2,25,28} precise GPS satellite orbit and clock estimates from the National Geospatial-Intelligence Agency (NGA) are considered truth reference data and are compared to broadcast GPS ephemerides archived by the Crustal Dynamics Data Information System (CDDIS). The focus in this work is on the validation of error model characteristics over time.

In this analysis, a set of data from January 4, 2015, to September 19, 2015, is processed, for Block IIF satellites including PRNs 1, 3, 6, 9, 25, 26, 27, and 30 (the total number of orbit data points is 572 520). The bias-plus-ramp model in Equation (4) is fit to truth-minus-broadcast data. Because NGA orbit data are only provided every 5 minutes, the fit interval is selected to be 30 minutes long. This is larger than the example 10- and 20-minute batch periods assumed in the next sections of this paper. It will provide conservative results because residual fitting errors are larger over 30 minutes than they would be over 20 minutes.

Based on this data, empirical parameter distributions for ^{i}*b*_{E} and ^{i}*g*_{E} are established. These are then bounded in the cumulative distribution function (CDF) sense^{20,21} using Gaussian distributions. The standard deviations of the overbounding Gaussian functions are given in Table 1 for the clock error contribution and for the three-dimensional orbit errors. Orbit errors are expressed in a local-level, satellite-fixed reference frame, in terms of the in-track, cross-track, and radial components.^{24} Because GPS satellites are at altitudes of about 20 000 km, user receivers near Earth’s surface are affected by ranging errors that are mostly due to the radial orbit and clock components. Fortunately, orbit radial and clock errors are significantly smaller than orbit in-track and cross-track components, as shown in Table 1.

To get a conservative estimate of the error parameter variance for ^{i}*b*_{E} while taking into account the worst-case geometry between an SV and a user receiver near Earth’s surface, we use the following equation:

5

where 0.24 is a multiplier accounting for the worst-case projection of non-radial orbit error components; , , and , respectively, are the variances of the clock, orbit radial, in-track, and cross-track error components contributing to ^{i}*b*_{E}. The same can be done for the variance of ^{i}*g*_{E}.

Equation (5) is a formulation similar to the ones found in Warren and Raquet,^{25} Department of Defense GPS NAVSTAR,^{27} and Heng^{29} but modified to evaluate variances instead of errors. The worst-case error projections are often referred to signal-in-space range error (SISRE) for ^{i}*b*_{E} and SIS range rate error (SISRRE) for ^{i}*g*_{E}. Our 9-month data analysis provides values of *σ*_{BE} = 0.61m and *σ*_{GE} = 1.8 · 10^{−4}m/s, which are both smaller than the assumptions mentioned above. It is therefore conservative, in the upcoming availability analysis, to use the values assumed in prior work^{1-3,11}: *σ*_{BE} = *σ*_{URA} = 1m and (in future work, we will process more data to check whether it is safe to assume lower *σ*_{BE} and *σ*_{BE} values).

In addition, the distribution of residual errors ^{i}*ε*_{RES,k} is plotted for the orbit radial and in-track components in Figures 3 and 4, respectively. ^{i}*ε*_{RES,k} is obtained by removing the best fit bias-plus-ramp model from the truth-minus-broadcast data. The figures show that ^{i}*ε*_{RES,k} is not negligible as compared to other error sources affecting carrier measurements. Since we do not know whether these residual errors are due to errors in the NGA truth data, we will account for them in the availability analysis, assuming

6

where the standard deviation of 0.056 m was obtained using the same formula as in Equation (5) but applied to ^{i}*ε*_{RES,k}, using the bounding standard deviations in the rightmost column of Table 1. Future work will include a detailed analysis of the time characteristics of ^{i}*ε*_{RES,k}, which, for now, is modeled as white noise for samples taken several minutes apart (we will use a 5-minute sample time interval in Section 6 later in this paper).

### 3.3 Multipath and antenna group delay error

Time-correlated raw code and carrier measurement errors due to multipath reflections are modeled as first-order Gauss-Markov processes (GMP) with time constant *T*_{MP}. Using this model, Chan et al.^{30} provide an expression for the steady-state CSC multipath error and dual-frequency antenna group delay variance in terms of *T*_{MP}, the Hatch filter (HF) smoothing time constant *T*_{HF}, the HF sampling interval *T*_{S}, and the raw code and carrier error variances and , respectively. is given by

7

where *α* = *T*_{HF}/*T*_{S} and.

To date, ARAIM analyses have assumed *T*_{HF} = 100 seconds and an associated elevation-dependent model for . For GPS, the following equation is given in Working Group C^{1}:

8

where

is expressed in meters for SV *i* at time *k*

^{i}*θ*_{k} is the elevation angle in degrees

*c*_{IF} is the ionosphere-free measurement combination multiplier:

*f*_{L1} and *f*_{L5}, respectively, are L1 and L5 frequencies.

Equation (8) accounts for the combination of multipath and antenna group delay.^{31,32} Dual-frequency antennas assumed in this work have much lower group delay than single-frequency antennas.^{33}

A similar model is provided in Working Group C^{1} for Galileo. Equation (8) was established using experimental data collected and analyzed in Murphy et al.^{18} using *T*_{HF} = 100 seconds and *T*_{S} = 0.5 second. We assume, based on experimental data analysis in Joerger et al.^{11} and Parkinson and Spilker,^{24} that is typically 100 times smaller than . Equation (7) still contains two unknowns: and *T*_{MP}.

Fortunately, a second set of experimental data was analyzed in Murphy et al.,^{19} which provides a different value of the CSC standard deviation using *T*_{HF} = 30 seconds. This second relationship provides the means to solve for both unknowns in Equation (7). It follows that

9

and

10

This model assumes a fixed *T*_{MP} even though *T*_{MP} values actually vary.^{34} In future work, we will evaluate the sensitivity of our proposed batch ARAIM integrity risk bound to *T*_{MP} values (preliminary analyses suggest that the bound is not very sensititive to *T*_{MP} because ^{i}*ε*_{T,k} is the dominating term in the carrier phase error budget). In parallel, we will use the method in Langel^{35} to derive an upper bound on the positioning error based on unknown, bounded values of *T*_{MP}.

### 3.4 Receiver noise

The CSC receiver noise standard deviation assumed in ARAIM^{1} is given by

11

where is expressed in meters for SV *i* at time *k*

^{i}*θ*_{k} is the elevation angle in degrees

Receiver noise is time uncorrelated. We follow the same approach as for multipath errors in Equations (7) to (9), which is simplified because *T*_{MP} no longer needs to be determined and *β* = 0. Raw code and carrier standard deviations are thus expressed as follows:

12

Throughout this section, we made the assumption that the error models are robust and accurate over the entire 10- to 20-minute–long batch interval. References are cited in Joerger et al.^{11} to support this assumption, but further experimental validation will have to be carried out in future work. Also, undetected and unrepaired cycle slips that impact both snapshot and batch ARAIM are not accounted for in this paper, the emphasis being on comparing snapshot to batch ARAIM. The only sources of error left unaddressed in Equation (2) are ^{i}*b*_{ρ,k} and ^{i}*b*_{φ,k}, which are dealt with in the next section.

## 4 BATCH ESTIMATOR AND DETECTOR DESIGN

In this section, the raw measurement error models derived above are incorporated in a batch estimator and in a batch MHSS RAIM fault detection algorithm that enables integrity risk evaluation. This raw-data–based implementation, or “raw batch,” is then re-formulated in a much more computationally efficient process, referred to as “batch ARAIM,” which only requires about three samples over 10 minutes to provide equivalent performance to the full raw batch.

### 4.1 Batch measurement equation

For each SV *i*, for *i* = 1,…,*n*, raw code and carrier measurements from times 1 to *q* are, respectively, stacked in *q* × 1 vectors ^{i}**ρ** and ^{i}**φ**. These vectors are then arranged in a batch measurement vector:

13

This vector can be expressed in terms of state variables and measurement error vectors as follows:

14

where, for an example dual-constellation GPS/Galileo implementation, **u** is a vector of positions and GPS and Galileo receiver clock offsets at all times

**η** is an *n* × 1 vector of cycle ambiguities

**s**_{ERR} is a 2*n* × 1 vector of constant error states

**0**_{a × b} is an *a* × *b* matrix of zeros

**v**_{T,E,RNM,φ} and **v**_{T,E,RNM,ρ} respectively, are raw code and carrier errors due to residual troposphere delay, residual clock and orbit ephemeris errors, receiver noise, and multipath.

Vector **u** is constructed using the following equations:

15

where *τ*_{GPS,k} and *τ*_{GAL,k} are the receiver clock offsets for GPS and Galileo, respectively, assuming that the time offset between the two constellations is unknown.

The vector of constant cycle ambiguities is defined as follows:

16

The vector of error states is given by

where **b**_{E} and **g**_{E} are the *n* × 1 vectors of constant clock and orbit ephemeris biases and gradients for all *n* satellites constructed following the exact same pattern as **η** in Equation (16).

Error states **s**_{ERR} are included in the state vector, not because their estimated values are of particular interest but because state augmentation is a practical way to incorporate measurement error dynamics. Prior knowledge on these error states is captured in a diagonal state information matrix (inverse of covariance matrix) with diagonal elements:

where**1**_{a × b} is an *a* × *b* matrix of ones.

This *a priori* knowledge of state estimate errors can directly be incorporated in the estimator by adding up information matrices as described in Section 2.1.2 of Joerger et al.,^{11} or it can be included by measurement vector augmentation, ie, assuming pseudo-measurements on **s**_{ERR} as in Joerger et al.^{11} and Joerger.^{36}

For clarity of presentation, the batch state coefficient matrices **G**, **H**_{N}, and **H**_{ERR} are constructed for example GPS/Galileo geometries where satellites are visible throughout the batch duration. Procedures to handle satellites coming in and out of sight are straightforward in batch implementations and are described in Joerger.^{36}

Position and receiver clock state coefficients for SV *i* at time *k* are arranged in the following vectors:

if SV *i* is a GPS satellite,

if SV *i* is a Galileo satellite, where line-of-sight vectors were defined in Equation (2). We then build geometry matrices over time, one satellite at a time:

17

and we stack these matrices for all satellites:

Next, the cycle ambiguity and error state coefficient matrices take the following forms:

. The measurement error covariance matrix **V** of vector requires additional steps to account for code/carrier measurement error correlation; both are described in Appendix B.

The 2*qn* × 1 vector represents nominal carrier and code biases for all SVs at all times for vector [**φ**^{T} **ρ**^{T}]^{T}. The elements of **b**_{φ} are significantly smaller than those of **b**_{ρ} because we assume that signal deformation does not affect the carriers. **b**_{φ} is included to account for non-Gaussian errors following the paired overbounding method.^{21} **b**_{ρ} and **b**_{φ} are unknown, but their elements are bounded by known values.

Next, we consider the batch measurement Equation (14), with an additive 2*qn* × 1 batch fault vector **f**, which will have to be detected. The resulting batch observation equation can be rewritten in a standard form as follows:

18

where **z** is the batch measurement vector

**H** is the batch observation matrix

**x** is the batch state vector

**v** is the batch measurement error vector:

**v**~*N*(**0,V**)**b** is the batch bias-bounded vector accounting for non-Gaussian components of errors:

### 4.2 Batch estimator design

The batch estimate for the state of interest (eg, for the vertical position coordinate, which is of primary interest in aircraft approach navigation) at the current time *q* is defined as follows:

19

where **s**_{0} is the 2*qn* × 1 vector of batch estimator coefficients. Subscript “0” designates the all-in-view or “full-set” solution (the same notations, with additional details, are used in Joerger et al.^{37}) The estimator is a weighted least-squares estimator modified using Equations (35) to (37) in Working Group C^{3} to reduce integrity risk. The full-set estimate error is noted , where *x* is the true value of the state of interest. *ε*_{0} is such that

20

where *b*_{0} is the impact on state estimation of nominal bias-bounded errors. The ARAIM error model in Working Group C^{1,2} assumes that the elements of **b** for CSC measurements can be bounded by a maximum value noted *b*_{nom}: *b*_{nom} = 0.75 m in MS2.^{2} In this paper, we conservatively assume that nominal ranging biases can be different for all satellites at all times, that they mostly affect code due to signal deformation, and that they also account for small biases in the common clock and orbit ephemeris paired bounding process. It follows that *b*_{0} can be bounded by

21

where || denotes the element-wise absolute value operator. Matrix **A** is defined as follows:

where **I**_{a} is the *a* × *a* identity matrix.

### 4.3 Batch detector design

A multiple-hypothesis solution separation (MHSS) batch RAIM method^{1-3,12,37} is adopted for detection of **f**. Let *h* be the number of fault hypotheses that need to be monitored against (refer to Working Group C^{3} and Joerger et al.^{37} for details on how to determine *h*). A set of mutually exclusive, exhaustive hypotheses *H*_{i}, for *i* = 0,…,*h*, is considered. Under *H*_{i}, a number *n*_{i} of measurements are simultaneously impacted by the fault. The fault-free subset solution, which excludes these *n*_{i} measurements, is written as follows: , where **s**_{i} is the 2*qn* × 1 vector of the subset solution’s batch estimator coefficients with zeros for elements corresponding to the *n*_{i} faulted measurements.^{37} Under *H*_{i}, the estimation error *ε*_{i} of is such that

22

where

23

The batch MHSS test statistics are then defined as follows:

24

Δ_{i} is normally distributed with variance .^{3} In snapshot ARAIM, hypotheses *H*_{i}, for *i* = 1,…,*h*, are defined by satellites being faulted or not. Batch ARAIM additionally opens the possibility for a sub-subset of data being either faulted or fault free over time. It is then worth asking whether a larger number of measurements being faulted over time necessarily cause a higher integrity risk. On the one hand, more faulted measurements may increase the estimation error, but, on the other hand, it may also provide more opportunities for detection.

### 4.4 Computationally efficient batch implementation

In this section, two additional parameters are introduced. Let *T*_{B} be the batch period, ie, the finite time interval over which measurements are processed: *T*_{B} determines the amount of change in SV geometry (eg, we will use *T*_{B} = 10 minutes and *T*_{B} = 20 minutes in the next sections of the paper). Also, let *T*_{S} be the sampling interval, ie, the time between raw samples within the batch (eg, *T*_{S} = 0.5 second).

Batch Equations (13) to (25) assume that multi-GNSS measurements are processed at regular *T*_{S} intervals over *T*_{B}. The dimensions of the resulting raw batch matrices and vectors are therefore very large, and matrix inversions, needed, for example, to calculate **s**_{0} in Equation (19), are extremely time-consuming. Alternatives to batches were considered, including Kalman filters (KF). But KF innovation-based RAIM requires computation of the *P*_{HMI} maximizing failure mode slope,^{38,39} and KF-based MHSS requires banks of KFs, which is cumbersome.^{12}

In response, we derive a computationally efficient, reduced batch implementation, illustrated in Figure 5, based on raw carrier phase and CSC measurements taken at few, infrequent sample times during *T*_{B}. This algorithm exploits the same two fundamental estimation principles as the raw batch: (a) Code noise is averaged out using CSC obtained from Hatch filters (HF); (b) redundant satellite motion is captured using raw carrier phase measurements at few, temporally separated sample times. This implementation is a straightforward augmentation of snapshot ARAIM that incorporates a few samples over time instead of a set of data at one instant in time.

Let *T*_{RBS} be the reduced batch sampling interval within the batch period *T*_{B}. Because in this reduced batch process, noise averaging is performed in separate HFs, *T*_{RBS} can be selected much larger than the raw measurement sampling period *T*_{S} (*T*_{RBS} >> *T*_{S}).

For example, in Section 5, we use *T*_{B} = 10 minutes and *T*_{RBS} = 5 minutes. The resulting number of sample times drops from *q* = 1200 for the raw batch (for *T*_{B} = 10 minutes, *T*_{S} = 0.5 second) down to *q* = 3 for the reduced batch.

The reduced batch ARAIM MHSS method is the same as for raw batch in Equations (13) to (25), except for two differences. First, Equation (14) becomes

25

where *q* is the number of sample times, much smaller for the reduced batch than for the raw batch

is a *qn* × 1 vector of CSC measurements

is a *qn* × 1 vector of residual ephemeris and tropospheric error and of CSC receiver noise and multipath error

Then, off-diagonal elements of the covariance matrix of vector account for the correlation of raw carrier measurements with CSC. Complete analytical expressions of the correlation between CSC and raw carrier measurements are given in Appendix C. It is worth noting that the 5-minute *T*_{RBS} value used in the next sections of this paper is very large as compared to the multipath correlation time constant (*T*_{MP} = 80 seconds) and to the HF smoothing time constant (*T*_{HF} = 100 seconds). It follows, as explained in Appendix C, that the expression of can be simply and accurately approximated using blocks of diagonal matrices.

In summary, the reduced batch ARAIM MHSS method presents the following characteristics.

a. It is much more computation and memory efficient than the raw batch.

b. Since it is also based on CSC, it only requires a minor augmentation of the existing snapshot ARAIM algorithms.

The reduced batch implementation in Figure 5 is used in Section 6. From now on, it is referred to as “reduced-batch ARAIM” or simply “batch ARAIM,” as opposed to the “raw batch” in Equation (14), which will be used shortly to analyze raw measurement error models.

### 4.5 Integrity and continuity risk evaluation

The integrity risk, or probability of hazardous misleading information *P*_{HMI}, is upper bounded as follows^{37}:

26

where *ℓ* is the alert limit (AL) that defines hazardous situations: in MS2,^{2} the vertical AL is *ℓ* =35m

*P*_{Hi} is the prior probability of *H*_{i} occurrence

*H*_{0} is the fault-free hypothesis

*H*_{i} for *i* = 1,…,*h* are the fault hypotheses corresponding to faults on subset measurement *i* (including single-satellite and multi-satellite faults)

Under fault-free hypothesis *H*_{0}, the detection threshold *τ*_{i} is set based on an allocated continuity risk requirement *C*_{REQ} (specified in Working Group C^{2,3}) to limit the probability of false alarms. *τ*_{i} can be defined as follows:

27

where the function *Q*^{−1}{} is the inverse tail probability distribution of the standard normal distribution.

The following section analyzes the probability distributions of *ε*_{0}, *ε*_{i}, and Δ_{i} for raw batch versus snapshot ARAIM and for varying values of *T*_{B}.

## 5 ANALYSIS OF MEASUREMENT ERROR MODELS AND OF SNAPSHOT VS BATCH APPROACHES

This section describes a step-by-step analysis to show that, given *T*_{B} of 10 minutes or larger, estimation performance improves substantially for batch ARAIM as compared to snapshot ARAIM. The raw batch described in Equation (14), based on raw measurements, is used throughout this section. Also, the computationally efficient reduced batch in Equation (25) is shown to be equivalent in integrity performance to the raw batch.

### 5.1 “Frozen geometry” analysis: Raw versus CSC measurement error model comparison

The comparison between snapshot ARAIM and raw batch ARAIM in Figure 2 illustrates that the main difference between the two approaches is the use of SV motion. Assuming a constant, “frozen geometry” should therefore provide identical performance for both implementations. This observation can be used to confirm that the raw measurement error models in Equations (1) to (12) are consistent with ARAIM CSC models in the ARAIM MS1-2.^{1,2}

To do so, the transient response of the Hatch filter (HF) must also be accounted for, because CSC error models in^{1,2}MS2[2], Murphy et al.^{19} and De Cleene^{20} are specified assuming that the HF has reached steady state. The steady-state CSC multipath error standard deviation is expressed in Equation (7) in terms of *α* = *T*_{HF}/*T*_{S}. In parallel, in appendix A of Joerger,^{36} the variance of the raw batch position estimate assuming frozen geometry is expressed in terms of the number of samples *n*_{S} (*n*_{S}=*T*_{B}/*T*_{S}). These two expressions are used to find that equivalent performance is expected when *T*_{B} = 2*T*_{HF}, ie, *T*_{B} = 200 seconds since *T*_{HF} = 100 seconds.^{19}

In the following Figures 6–8, the fault-free means and standard deviations of the random variables in Equation (26) are evaluated over 24 hours at an example Chicago location (25.5^{°}N, −80.1^{°}E), assuming dual-frequency measurements from GPS and Galileo. Nominal simulation parameters are further described in the next section where they become relevant. The means and standard deviations for raw batch ARAIM are *σ*_{0}, *b*_{0}, *σ*_{i}, *b*_{i}, and *σ*_{Δi}, which are defined in Equations (20) to (25). The same quantities for snapshot ARAIM are, respectively, noted , and . The legend is the same in all three figures. Figure 6 shows the results for *T*_{B} = 200 seconds, for frozen geometries. Snapshot versus raw batch results match closely, which confirms that the raw measurement error models derived earlier are consistent with ARAIM assumptions in MS2.^{2}

### 5.2 Impact of satellite motion

The fault-free estimation error means and standard deviations are plotted in Figure 7 assuming the same simulation parameters as in Figure 6, but the geometry is no longer frozen. Figure 7 shows the substantial impact of satellite motion over *T*_{B} = 200 seconds on estimation performance. It confirms that bias observability using SV motion is a performance-driving principle. This effect is further accentuated in Figure 8 (especially for *b*_{0} and *b*_{i}), where the batch interval is increased from *T*_{B} = 200 seconds to *T*_{B} = 600 seconds.

Figure 9 also shows the significant *P*_{HMI} reduction obtained using raw batch ARAIM as compared to snapshot ARAIM. In Figure 9, the fraction of time where the *P*_{HMI} curves are below the dashed line is the availability. In this case, availability is 99.7% for raw batch ARAIM and is much lower for snapshot ARAIM (65%). In addition, Figure 9 shows that raw batch ARAIM and reduced batch ARAIM have overlapping *P*_{HMI} curves, even though reduced batch ARAIM is much more computation and memory efficient.

The integrity risk *P*_{HMI} is plotted in Figure 9 for raw batch ARAIM (*T*_{B} = 600 seconds, *T*_{S} = 0.5 seconds) and for reduced batch ARAIM (*T*_{B} = 600 seconds, *T*_{RBS} = 300 seconds), over 24 hours, at the example Chicago location. The dashed horizontal line in Figure 9 is the integrity risk requirement *I*_{REQ}, specified in Working Group C^{2} to be *I*_{REQ} = 0.98 · 10^{−7} for the vertical position coordinate. Figure 9 also assumes that the vertical alert limit in Equation (24) is *ℓ* = 10m, that prior probabilities of satellite and constellation faults, noted *P*_{Hi} in Equation (24), respectively, are *P*_{sat} = 10^{−5} and *P*^{const} = 10^{−8}, and that the continuity risk requirement in Equation (27) is *C*_{REQ} = 3.9 · 10^{−6}.

## 6 BATCH ARAIM AVAILABILITY ANALYSIS

This section evaluates the global availability of LPV-200 navigation requirements to support localizer precision vertical aircraft approach operations down to 200 ft above the ground, using ARAIM with dual-frequency measurements from GPS and Galileo. Nominal simulation parameters are defined in Working Group C^{2} and include the following:

a5

^{°}satellite elevation maskvalues for all ARAIM measurement error model parameters in Equations (1) to (10) and (21)

navigation requirements: I

_{REQ}= 0.98.10^{-7},*C*_{REQ}= 3.9 · 10^{−6},*ℓ*= 35mprior probabilities of faults:

*P*_{sat}=10^{−5},*P*_{const}= 10^{−4}reduced batch period and sampling interval:

*T*_{B}= 10 minutes,*T*_{RBS}= 5 minutesnominal constellations, comprising 24 GPS satellites and 24 Galileo SVs

^{40}

These parameters are modified below to evaluate performance sensitivity. For example, to account for potential satellite outages, “depleted” constellations of 24-1 GPS satellites and 24-1 Galileo SVs are simulated.^{40} “Optimistic” constellations of 27 GPS and 27 Galileo SVs are considered as well.^{40} Additional requirements, including effective monitor threshold (EMT) and fault-free accuracy requirements,^{2} are included in the simulation but not discussed in this paper as they only have a minor impact on overall availability.

### 6.1 Batch vs snapshot ARAIM availability using depleted constellations

Figures 10 and 11 display availability maps for a 10^{°} × 10^{°} latitude-longitude grid of locations, for “depleted” 24-1 GPS and 24-1 Galileo constellations, for satellite geometries simulated at regular 10-minute intervals over a 24-hour period. Availability is computed at each location as the fraction of time where the *P*_{HMI} bound is lower than *I*_{REQ}.

In the figures, availability is color coded: White corresponds to a value of 100%, and black represents 80%. Constant availability contours are also displayed. The gray areas in Figure 10 indicate that snapshot ARAIM is clearly outperformed by batch ARAIM in Figure 11. Even more substantial differences will be observed in the next subsection.

The worldwide availability metric given in the figure captions is the weighted coverage of 99.5% availability: Coverage is defined as the percentage of grid point locations exceeding 99.5% availability. The coverage computation is weighted at each location by the cosine of the location’s latitude because grid point locations near the equator represent larger areas than near the poles. Figures 10 and 11 show that the coverage of 99.5% availability increases from 87.2% for snapshot ARAIM to 93.3% for batch ARAIM, assuming *T*_{B} = 600 seconds for *ℓ* = 35 m. For context, based on the assumptions described in Section 4 of MS3,^{3} coverage of 99% availability using Satellite Based Augmentation Systems (SBAS) is 11% as of 2015 and may reach 45% in 2026 using dual-frequency SBAS.

### 6.2 Batch vs snapshot ARAIM availability for a 10-m alert limit

Figures 12 and 13 evaluate the potential of batch ARAIM to meet requirements that are more stringent than LPV200, including a tight alert limit: *ℓ* = 10 m. One key assumption to achieve high availability is that the prior probability of constellation-wide faults must be reduced to *P*_{const} = 10^{−8}. This assumption is also made for GPS in en route operation using horizontal ARAIM (see appendix B in Working Group C^{2}) and may be justifiable using the guidelines in Walter and Blanch.^{41} These guidelines exploit additional information from the Air Navigation Service Provider (ANSP) ground segment. In Figures 12 and 13, availability maps assume nominal 24-satellite GPS and Galileo constellations. The same color code as in Figures 10 and 11 is employed. Again, batch ARAIM with *T*_{B} = 1200 seconds in Figure 13 provides a dramatic improvement as compared to snapshot ARAIM in Figure 12.

Table 2 lists worldwide coverage of 99.5% availability and of 95% availability (given in parentheses), for the above configurations, for *T*_{B} = 600 seconds and assuming *T*_{B} = 1200 seconds, and for depleted, nominal, and optimistic constellations. The table quantifies the global performance improvement brought by batch ARAIM as compared to snapshot ARAIM (snapshot ARAIM results match the coverage reported in tables 2 and 3 in Working Group C^{2} for *ℓ* = 35 m). This improvement comes at the cost of slightly higher computation and memory loads.

Further improvement can be obtained by lengthening the batch duration. For example, for *T*_{B} = 1800 seconds, *ℓ* = 10 m, *P*_{const} = 10^{−8}, and assuming a nominal constellation, the coverage of 99.5% availability reaches 95% as compared to 75.5% using *T*_{B} = 1200 seconds and to 0% for snapshot ARAIM. But the batch duration is limited by mission duration and by the period of validity of the error models. Another source of improvement is the potential future addition of other GNSS constellations (eg, GLONASS and BeiDou), which will increase the accumulated impact of redundant satellite motion. In this case, new measurement error and fault models would have to be derived and validated because current assumptions on GPS and Galileo performance may not hold for other constellations (eg, the definition of GLONASS satellite faults is not rigorously consistent with that of GPS^{42}).

It must also be noted that, for aircraft navigation standards that are more stringent than LPV-200 and that include *ℓ* = 10 m, additional requirements are typically involved, eg, on the communication link between ANSP ground segment and aircraft. Such considerations are beyond the scope of this paper.

## 7 CONCLUSION

In this paper, a new ARAIM integrity monitoring method was devised, which exploits the motion of satellites from multiple GNSS constellations. Raw measurement error models over time were established and were shown to be consistent with conventional snapshot ARAIM assumptions on carrier-smoothed code measurements in MS1-3.^{1-3} These raw measurements were then incorporated in batch estimation and solution separation fault detection processes that constitute minor augmentations to snapshot ARAIM.

The impact of satellite motion on batch ARAIM was analyzed as a function of batch period and then quantified globally in comparison with conventional snapshot ARAIM. The proposed batch ARAIM implementation is slightly more computation and memory expensive than snapshot ARAIM. However, it can provide dramatic performance improvements both when aiming to achieve LPV-200 requirements using depleted constellations and when trying to meet a much more stringent 10-m alert limit.

The next step of this research is to pursue the validation and refinement of raw measurement error models over time using large amounts of experimental data. Further measurement error and performance sensitivity analyses will be carried out considering additional GNSS (GLONASS and BeiDou).

## HOW TO CITE THIS ARTICLE

Joerger M, Pervan B. Multi-constellation ARAIM exploiting satellite motion. *NAVIGATION-US*. 2020;67:235–253. https://doi.org/10.1002/navi.334

## ACKNOWLEDGEMENTS

The authors would like to thank the Federal Aviation Administration (FAA) for their support of this research. However, the opinions in this paper are our own and do not necessarily represent those of any other person or organization.

## APPENDIX A: ILLUSTRATIVE “STATIC SURVEYING” EXAMPLE USING MULTIPLE SATELLITES

The analysis described in this appendix and illustrated in Figure 1 is a direct extension of the “static surveying” problem given in Hwang.^{17} In this case, the measurement vector is augmented to include differential carrier phase signals ^{i}Δ*φ*_{0} and ^{i}Δ*φ*_{1} (at times 0 and 1 as indicated by subscripts) not only from a single SV but also from satellites *i*, for *i* =1,…,*n*. The measurement equation becomes

A.1

where

and Δ*x*_{V} is the vertical baseline distance to be estimated (displayed in Figure 1)

**I**_{n} is an *n* × *n* identity matrix

^{i}*η* carrier phase cycle ambiguity for SV *i*

The carrier phase measurement noise vector is assumed zero-mean normally distributed with covariance matrix . The least-squares estimate error variance for Δ*x*_{V} is given by

A.2

where **0**_{a} × ^{b} is an *a* × *b* matrix of zeros

Using a popular matrix inversion formula, becomes

A.3

which, using the above definitions of **G**_{0} and **G**_{1}, can be written as follows:

A.4

where ^{i}*θ*_{0} is the zenith angle for SV *i* at some initial time 0 ^{i}*θ*_{1} is the zenith angle for SV *i* at a later time 1

Equation (A.4) is discussed in the body of the paper (under Equation (1)).

Equation (A.4) can be further analyzed considering small angular variations over an infinitesimally small time interval *δt*. The reasoning is pursued in the information form because information (inverse of variance) can directly be added up, and it avoids discussing the special case of ^{i}*θ*_{0} = ^{i}*θ*_{1} when no information is provided by carrier only. The change in square root information can be expressed as follows:

A.5

Using a first-order Taylor series approximation for small angular variations about *θ* of the second term in the sum and dividing both sides by *δt*, Equation (A.5) becomes

A.6

Low-elevation satellites have a zenith angle ^{i}*θ* approaching 90^{°}. Equation (A.6) shows that the information contribution is larger for low-elevation satellites than for high-elevation SVs, both because of the sin^{i}*θ* term and because is larger at low elevations.

## APPENDIX B: BATCH MEASUREMENT ERROR COVARIANCE MATRIX DESCRIPTION

In this appendix, consistent with Section 4.1, the measurement noise covariance matrix is constructed, for example, GPS/Galileo geometries where satellites are visible throughout the batch duration. The measurement error covariance matrix is the sum of two contributions:

**V**_{T,E} is constructed using blocks of four identical diagonal matrices accounting for the fact that code and carrier measurements are impacted by the same troposphere and ephemeris errors. **V**_{T,E} can be expressed as follows:

where in Equation (6) and **V**_{T} is a diagonal matrix (assuming that errors from different satellites are independent from each other as in Working Group C^{2}) with the ((*i* − 1)· *q + k*)th diagonal term being as described in Equation (3) and in the text that follows.

For the reduced batch, **V**_{RNM} is described in Appendix C. For the raw batch, **V**_{RNM} is block diagonal, each block corresponding to observations from the same SV over time (to facilitate this process, we grouped measurements satellite per satellite in Equation (17)). Within each block, off-diagonal components capture the time correlation due to multipath error, which is modeled as a first-order Gauss-Markov process. Thus, within the block corresponding to satellite *i*, the off-diagonal element accounting for the error correlation of sample time *t*_{k} and is expressed as , where Δ*t*_{kl} =| *t*_{k} − *t*_{l}| and where is the (code/carrier) multipath error variance defined in Equation (9) and *T*_{MP} is the correlation time constant in Equation (10). The quantities and in Equation (12) are also added to the diagonal elements of **V**_{RNM} to account for code phase and carrier phase uncorrelated receiver noise, respectively.

## APPENDIX C: MEASUREMENT ERROR CORRELATION IN REDUCED BATCH IMPLEMENTATION

This appendix provides analytical expressions of the correlation between CSC and raw carrier measurements. These expressions are needed to express the measurement error covariance matrix of vector in Equation (25). is constructed similar to **V** = **V**_{T,E}+**V**_{RNM} in Appendix B, but the dimensions of **V**_{T,E} and the term **V**_{RNM} differ. We use the notation .

Equations are provided assuming time-correlated raw code and carrier phase noise modeled as first-order Gauss-Markov processes with a time constant *T*_{MP}:

C.1 C.2

where *v*_{ρ,k} and *v*_{φ,k} are white random sequences, and

C.3

These models are used to account for multipath (assuming *T*_{MP} = 80 seconds) and receiver noise (assuming *β* = 0).

We can express the raw carrier phase measurement time correlation over one time step as follows:

C.4

where *E*{} is the expected value operator and where we defined . For the reduced batch implementation, measurements are taken at infrequent time intervals, much larger than the raw data sampling period. Let us define the number *K* of raw sample times between two successive reduced batch samples. *K* is given by

C.5

For example, for values used in performance analysis, *K* = 600. Over *T*_{RBS}, Equation (C.4) becomes

C.6

For *K*=600, *T*_{S} = 0.5 second, and *T*_{MP} = 80 seconds, we can use the following approximation: .

In addition, the HF CSC measurement equation is expressed as follows:

C.7

where

C.8

Equation (C.7) is used to evaluate the correlation between and *ε*_{φ,k}, which can be written as follows:

At steady state, , we obtain the following expression:

C.9

The carrier phase measurement variance changes slowly over time: It is elevation dependent and varies with satellite motion. However, it is constant over *T*_{S} = 0.5 second, so that Equation (C.9) becomes

C.10

where

C.11

The same assumptions and the same type of derivations are used to establish the following results, which are needed to populate :

C.12 C.13

One of the most challenging terms to derive is the time correlation of CSC measurements over time. Time correlation is caused by both multipath and the HF, which keeps a finite memory of past measurements. A full derivation is not provided here to limit the length of the paper, but the following result was obtained:

C.14

where we used the notations

C.15 C.16

It is worth noting that when *T*_{RBS} >> *T*_{MP} and *T*_{RBS} >> *T*_{HF}, which is the case in our performance analysis, the correlation terms in Equations (C.6) and (C.12) to (C.14) approach 0. The only non-negligible terms in are the known quantities , and in Equation (C.10).

## Footnotes

**Funding information**Federal Aviation Administration

- Received March 14, 2019.
- Revision received June 14, 2019.
- Accepted August 5, 2019.

- © 2020 Institute of Navigation

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.