## Abstract

Our objective in this paper is to assess all errors in NavIC’s pseudorange so that they can be mitigated or accounted for appropriately in the dynamics model to achieve high-accuracy navigation. Multipath errors in the L5 and S1 pseudoranges are different; hence, oscillatory multipath is removed using sidereal repeatability for these frequencies separately for accurate iono delay estimation. Iono delay estimates using uncompensated pseudoranges are 7-m greater than the normal night time estimates due to multipath. When oscillations are removed, the difference reduces to ≤3 m. Dual-frequency iono delay estimates using code phase are within two-sigma of the grid-based estimates. A comparison with the NTCM-GL and Klobuchar iono delay estimates is presented for the low-latitude region in India. The ephemeris line-of-sight errors blended with the unknown multipath bias errors are estimated using a low-pass filter and are found to lie at the edges of the two-sigma signal-in-space error.

## 1 INTRODUCTION

NavIC, short for Navigation with Indian Constellation, is a new regional navigation system for India.^{1} The satellites mean longitudes, maximum eccentricities, inclinations, and PRNs are presented in Table 1. There are publications concerned with specific error sources of NavIC signals,^{2-4} but this paper presents the first detailed investigation of its errors in pseudorange and carrier phase observables in its two frequencies. NavIC satellite signals of PRN1 in a geosynchronous orbit (GSO) are analyzed and its short-term clock stability is studied in Thoelert et al.^{2} Assessment of broadcast ephemerides of NavIC satellites by comparing with orbits estimated with laser ranging is provided in Montenbruck et al.^{3} The comparison includes assessment of the clock bias and range errors for two NavIC satellites. Dual-frequency (DF) iono delay estimate using NavIC signals is studied in Desai and Shah,^{4} and it is concluded that the DF estimate is close to the iono delay interpolated from the vertical iono grid estimates transmitted in the NavIC message. The authors did not encounter significant multipath^{5} but there is a ∼3-m difference in the DF iono delay estimate and grid-based estimate. In our research, the errors in the observables caused by NavIC satellites’ clock and instrument delays, multipath at the receiver site, atmospheric delays, and ephemeris errors are all estimated. The carrier phase wind-up error is not compensated for in this work. Accord’s NavIC/GPS antenna and the receiver used in this study are shown in the left and right plots of Figure 1. The receiver provides raw navigation data, ie, subframe (SF) level data from which the primary and secondary navigation data transmitted by each satellite can be read. The error corrections, PVT algorithms, and internal hardware details of the receiver are not known to the user.

This paper is organized as follows. In Section 2, time-varying satellite clock corrections, instrument delays for dual-frequency and single-frequency users, and relativistic corrections, parsed from the NavIC messages, are illustrated. Code phase multipath errors are analyzed in Section 3. It is known that the multipath error have a non-zero mean oscillations. Time-varying frequencies of multipath errors caused by horizontal reflectors surrounding the receiver along with the corresponding periods versus elevation angles of a GEO and a GSO NavIC satellite are illustrated. Oscillatory multipath errors are eliminated with well-known sidereal repeatability.^{6,7} For real-time removal of multipath error, its spectrum similarity with signal-to-noise ratio (C/N_{o}) is explored, following Sleewaegen.^{7}

Section 4 deals with single-frequency estimation of iono delay in NavIC signal and its variance. For single-frequency NavIC users, iono delay at the reference nodes provided by NavIC signals at an interval of 288 seconds^{4} is compared with the delay provided by GAGAN (GPS and GEO Augmented Navigation) system for aviation in India. The iono delay at an ionospheric pierce point (IPP) is interpolated from the surrounding grid iono delays using a bilinear relationship. However, this method ignores variances of the estimates at the nodes. So following Lapucha and Huff,^{8} a technique is introduced based on variance-dependent weights to estimate the IPP iono delay and compared with the standard interpolation scheme. Klobuchar and NTCM-GL (Global Neustrelitz TEC Model)^{9} based slant iono delays are also compared with the grid-based estimates.

While iono delay estimates of the grid nodes provided by NavIC signals are noise-free and multipath-free, they are not current since they are updated every 288 seconds,^{4} and they may have considerable variance. These values can be compared with the instantaneous iono delays calculated from dual-frequency observables. This is shown in Section 5 where iono delay is first estimated using pseudorange observables at two frequencies. Julien et al^{10} have shown that this, however, introduces a blended DF multipath error in the estimate. Also, to eliminate the amplified DF code noise, DF carrier phase observables are often used, which introduces ambiguities.^{10} A complementary filter therefore is used in Section 5 to estimate this ambiguity. But the multipath bias error is shown to still corrupt the iono delay estimates. We also show how the multipath oscillations in NavIC signals, if not removed, affect the iono delay estimate and illustrate smoother estimates with their removal. Finally, we compare these estimates with single-frequency grid-based, Klobuchar and NTCM-GL estimates considered in Section 4. Tropospheric delay is considered in Section 6 using the Collins model.

Section 7 presents methods to estimate (a) discrete ephemeris errors of NavIC satellites at orbit update epochs and their line-of-sight components, (b) unknown slowly-varying multipath bias summed with the line-of-sight component of the ephemeris error and range and receiver clock bias estimates errors, and (c) variance of the Accord receiver electronic noise. A low-pass filter of bandwidth 20 times the orbit update frequency is used for the bias estimation, and a complementary high-pass filter is used to estimate the noise variance. The LOS-components of the ephemeris errors in (a) and the biases in (b) are compared with the user range accuracy (or signal-in-space) standard deviations provided in the NavIC signals, and found to be, respectively, within +/− *σ* and 2*σ* of the signal-in-space, implying that the residual unknown bias in pseudorange observables corrected for all errors as shown in the preceding sections is small, within statistical uncertainty of the signal-in-space error. Section 8 concludes the paper.

### 1.1 Measurement equations

The pseudorange and carrier phase measurement equations at L5, composed of true range *r*, user clock bias Δ*t _{r}*, satellite clock bias Δ

*t*, ephemeris error component

_{s}*E*along the satellite-receiver line-of-sight, ionospheric delay

*I*, tropospheric delay

*T*, and satellite and user receiver instrument delays

*k*and

^{s}*k*, are given by

^{r}1 2

Carrier phase multipath *m*_{l5} magnitude is less than a decimeter whereas *M*_{l5} has meters of delay.^{11,12} Similarly, the carrier phase random noise *β*_{l5} is less than the code phase random noise *ν*_{l5} by two orders of magnitude. *N*_{l5} represents an unknown ambiguity. Before these measurements can be used to determine the receiver position and its clock bias Δ*t _{r}* with some desired accuracy, the errors, shown in the right sides of Equations (1) and (2), need to be removed or mitigated.

## 2 NAVIC SATELLITES’ CLOCK CORRECTIONS AND INSTRUMENT DELAYS

To begin with, the satellite clock bias Δ*t _{s}* is removed using the quadratic coefficients (

*a*

_{ƒ0},

*a*

_{ƒ1},

*a*

_{ƒ2}) in the navigation message using the formula (Appendix A of ISRO

^{1}and Appendix C.1 of Farrell

^{13}):

3

where *t _{oc}* is the time of clock, Δ

*t*is the relativistic correction, and

_{rel}*t*is the time at which the clock corrections are calculated by the user. Methods to compensate for the satellite instrument delays ( and ) for single- and dual-frequency users are presented.

_{sv}### 2.1 Dual-frequency users

The *a*_{ƒ0} parameter in Equation (3) inherently contains additional terms to compensate for the satellite instrument delays in iono-free (IF) pseudorange combination, defined by Equation (4). [Sec. 6.2.1.5 in Ref. 1]

4

where . Thus, DF users compensate for the satellite instrument delays through satellite clock corrections by forming iono-free pseudoranges.

### 2.2 Single-frequency users

Single-frequency users need to use the total group delay (*T _{GD}*) parameter in the navigation message to remove the satellite instrument delays from the L5 and S1 pseudoranges using Equations (5) and (6), respectively [Sec. 6.2.1.5 in Ref. 1].

5 6

Δ*t _{s}* is calculated using Equation (3), and

*T*is read from the navigation message. If the instrument delays are ignored, the range errors for S1 and L5 users for a typical value of

_{GD}*T*= 1.86 nanoseconds will be 0.56 m and 2.5 m, respectively.

_{GD}### 2.3 Illustration

The plots of satellite clock corrections for PRN3 (GEO) and PRN9 (GSO) are shown in the left and right columns of Figure 2.

The clock coefficients *a*_{ƒ0} and *a*_{ƒ1} are shown in the first two rows for about 24-h period with orbit updates every 2 hours. *a*_{ƒ2} is zero for the duration we studied. The clock offsets *a*_{ƒ0} for PRN3 and PRN9 are both negative as seen from the first row. The clock bias is within −0.2 to −0.1 milliseconds. The relativistic correction is sinusoidal [Section C.1 in Ref. 13], as shown in the third row, with eccentric anomaly of the satellite varying from −180^{◦} to 180^{◦} in 24 hours, with an amplitude of about 6 nanoseconds proportional to the eccentricity of the orbits. The total group delay parameter, for single-frequency users, is −1.397 nanoseconds and −1.8625 nanoseconds for PRN3 and PRN9, respectively, as shown in the fourth row of Figure 2. The usage of this parameter in correcting the satellite clock bias shows a difference with the Accord receiver’s estimate; this is investigated further in the next paragraph. The resulting satellite clock bias (Δ*t _{s}*) is in the range of 0.1 to 0.15 milliseconds, consistent with the assessment of clock bias for PRN1 and PRN2 in Montenbruck et al.

^{3}The maximum correction for the study duration at the time of clock updates is around 10 nanoseconds (3 m) as seen from the last plot. The discrete jumps seen in this plot are due to inherent changes in

*a*

_{ƒ0}at orbit updates and not due to the effect of

*a*

_{ƒ1}for the 2-h duration prior to the orbit updates. If the later was the case, then

*a*

_{ƒ1}= 0.03 ns/s would cause a clock offset of 216 nanoseconds (much greater than the 10 ns) for the 2-h duration.

Differences between the satellite clock correction for L5 users (using Equation (5)) and the Accord estimate are shown in Figure 3 for PRN3 (left) and PRN9 (right). The top row shows the difference with Accord when *γ* = 1 (which is incorrect) in Equation (5), and the bottom row shows the difference when *γ* = 4.48 (which is correct for NavIC frequencies). As seen from the top row, the difference is negligible (in millimeters) for both the satellites for *γ* = 1, whereas we get a significant difference of 1.46 m and 1.94 m for PRN3 and PRN9, respectively, for *γ* = 4.48. This may happen if L5 and S1 corrections applied within the receiver or reported in the file are inadvertently interchanged. A discrepancy in broadcast TGD is reported in Montenbruck et al^{3} for PRN1 and PRN2 where better accuracy is found if TGD is not applied for the L5 observable.

## 3 CODE PHASE MULTIPATH ERROR

In this section, multipath errors affecting the pseudorange observables are investigated. In the first part, we explore multipath error’s frequency and time period due to a single dominant ground reflection and show its correlation with C/N_{0}. In the second part, we illustrate repetition of the multipath error every sidereal day for a stationary receiver. We use this temporal correlation to remove the oscillatory part of the multipath error and illustrate the results.

### 3.1 Multipath error frequencies and periods due to horizontal reflectors

Multipath error frequency caused by ground reflection at L5 frequency is [Equation 5.164, Ref. 14]:

7

where *h* is the height of the antenna from the ground or horizontal reflector, *E* is the elevation angle of the satellite, and *λ*_{l5} is the wavelength of the carrier. The height of the antenna is 7.7 m.

The frequencies and periods for two satellites, PRN5 (GSO) and PRN3 (GEO), for 24 hours are shown in Figure 4. PRN5’s frequency and period are shown in the first row of Figure 4. It has a maximum frequency of 1.5 mHz, which is also the case for other GSO satellites. The minimum period is 12 minutes (shown in the middle column). The plot of frequency versus elevation angle, shown in the last column, is a closed curve because the elevation angle rises to a peak value and then decreases as the satellite sets. Similar plots for PRN3 are shown in the last row of Figure 4. As GEO satellite elevation angle changes are less than those for GSO satellites, the maximum frequency for PRN3 is limited to 0.10 mHz, seen in the first column of the last row. The GEO satellites have a minimum period of about 100 to 180 minutes, with PRN3 having a minimum period of 180 minutes, as seen in the middle column. The closed curve of PRN3’s frequency versus elevation angle, shown in the last column last row, is symmetric about the x-axis as compared to that of PRN5. The multipath error time periods calculated analytically in this section are validated in Section 3.4 using the estimate of multipath obtained from the pseudorange and carrier phase observables (calculated in Section 3.2).

### 3.2 Sidereal multipath error oscillations

Code phase multipath error is estimated using a linear combination^{11,15} of pseudorange (Equation (1)) and carrier phase (Equation (2) in terms of length) observables of two frequencies. However, this estimate is corrupted with receiver noise and unknown, real ambiguity terms. This is evident by observing the measurement equation of the code phase multipath error at frequency 1 (frequencies are denoted by subscripts 1 and 2):

8

where denotes the code phase multipath blended with the receiver noise at frequency 1 and *K*_{1} is a blend of dual frequency ambiguities expressed in lengths. *K*_{1} is a constant as long as there are no cycle slips.

#### 3.2.1 Illustration

A plot of the sum for all PRNs on 22 April 2018 is shown in the top left of Figure 5. Sinusoidal oscillations at diverse frequencies caused by the changes in elevation angles of the NavIC satellites are visible in the plot, especially for PRN5 (red) with a peak-to-peak amplitude of around 3 m for an elevation angle change of ∼14^{◦} seen in the bottom left plot of Figure 5. PRN4’s elevation angle (blue) changes even more (over a range of 30^{◦}), so the frequency of oscillation of the biased multipath is even higher than that of the oscillations in the PRN5 signals. Consistent with these observations, the multipath (plus a constant) of PRN7, a GEO satellite, the elevation angles of which change by 2^{◦} to 3^{◦}, oscillates at a low-frequency with min-max amplitude difference of 5 m. Another GEO satellite, PRN6, also has minimal oscillations due to its near constant elevation angle. In addition to elevation angles, the local surroundings (nearby buildings) affect multipath errors significantly. The sky plot, shown in the top right plot of Figure 5, reveals that PRN4, 5, and 7 covering azimuth angles from 45^{◦} (northeast) to 135^{◦} (southeast) have higher amplitude multipath oscillations than the PRN2, 3, and 6 satellites have. The bottom right plot in Figure 5 is discussed below.

### 3.3 Oscillatory multipath error removal

A multipath error (*M*) is likely to consist of some bias (*M _{bias}*) and an oscillatory (

*M*) component

_{osc}^{11,16}(subscript

*l*5 is dropped henceforth for clarity). In the case of stationary receivers, multipath error almost repeats the next sidereal day (Sec. 15.8.3, Ref. 17) and this repetition is illustrated in the bottom right plot of Figure 5 over three consecutive days (April 22, 23, and 24). The plots look alike except for a discrete jump on 23 April, between 00:00 and 01:00 hours, possibly due to cycle slips, but the rest of the plot follows a similar pattern as the other 2 days. Because the errors are repeatable, the previous day’s time-differenced (

*pd*superscript for previous day) observable using Equation (8) can be subtracted from the current day’s L5 pseudorange observable to mitigate the oscillatory error. Mathematically, the corrected pseudorange observable (with a superscript

*c*) will be

9

where *j* is the current day’s time index, related to previous day’s time index *i* by *j* = *i* + 86164. *M _{j,osc}* and , which are a sidereal day apart, are nearly the same and hence nearly cancel each other. They are not identical because the orbit period is not exactly 86164 seconds, which is for a circular orbit, whereas the real GEO and GSO orbits have time-varying eccentricity countering solar radiation pressures continuously. Hence, the residual oscillatory error (

*δM*) in Equation (9) will be small but non-zero. The multipath bias term

_{j,osc}*M*in the next day pseudorange observable also remains. Further, the previous day’s oscillatory multipath error at the first time instant, whereat the time-differencing operation started, remains as an error. Observing that there are two random code phase noises in Equation (9), the standard deviation of the random sum noise in the compensated pseudorange observable will be times the

_{j,bias}*σ*of the raw code phase noise.

#### 3.3.1 Illustration

Figure 6 illustrates elimination of the multipath oscillations in the PRN5 pseudorange. The top plot of Figure 6 shows the multipath error plus a constant bias in the raw pseudorange observable calculated using Equation (8). The middle plot shows mostly successful removal of the oscillations using the above algorithm, though some residual oscillations remain, as anticipated above. The 13^{◦} change in elevation angle of PRN5 is shown in the bottom plot. Thus, the algorithm for removal of the oscillatory part of the multipath error is validated.

Multipath error and the associated ambiguities in the raw S1 pseudorange are shown in the top plot of Figure 7. The oscillations visible in the top plot are removed using sidereal repeatability, and the results are shown in the bottom plot of Figure 7. Two points are evident from this figure: (1) Comparing the top plots of Figure 6 and Figure 7, we observe that the multipath errors of L5 and S1 frequencies are different; hence, they will not cancel each other in DF code phase iono delay estimation, affecting the accuracy (more on this in Section 4); (2) a higher code (random) noise *σ* = 29.3 cm is observed in the raw S1 pseudorange as compared to raw L5 pseudorange *σ* = 18.3 cm. The *σ* of the raw pseudorange measurement is calculated by dividing the standard deviation of the observable in Equation (9), shown in the Figure 6 and Figure 7, in the interval where they are minimal oscillations, by .

### 3.4 Multipath error and C/N_{0}

Breivik et al^{18} and Sleewaegen^{7} showed a strong spectrum correlation between multipath error and signal-to-noise ratio and developed an algorithm to estimate the former with the latter using a single-frequency receiver and mitigated it by half. We will explore this spectrum similarity for NavIC signals now. Multipath error calculated using Equation (8) for PRN5 (GSO) observable at the Hub building on 22 August 2018 is shown in the top plot of Figure 8. There are four low-frequency cycles in the estimate in the 3-h interval from 12:00 to 15:00 hours local time, yielding an average period of 180 minutes∕4 = 45 minutes. This matches with about the average of the periods shown in the last plot of Figure 8 from 12:00 to 15:00 hours of the multipath error oscillations calculated with Equation (7) for *h* = 7.7 m for elevation angle shown in the third plot of Figure 8. The C/N_{0} (averaged every 120 seconds) of NavIC signal from this satellite is shown in the second plot. A low-frequency spectrum similarity with the multipath error is evident on comparison. Its online estimation and mitigation following Breivik et al^{18} will be pursued in the future.

## 4 IONOSPHERIC DELAY ESTIMATION AND ITS VARIANCE FOR SINGLE-FREQUENCY USERS

This sections deals with mitigation of iono delay for single-frequency (SF) users. NavIC provides grid-based vertical iono corrections in its navigation message for SF users. We investigate the vertical delays and their confidence limits at the grid points as well as at the ionospheric pierce point (IPP). GAGAN, a geostationary augmented GPS navigation system similar to WAAS for SF users, transmits vertical delays at the same grid points as the delays provided by NavIC. A comparison of these two vertical grid corrections is presented in this section. Klobuchar coefficients are also transmitted in NavIC in case the grid corrections are unusable, though its efficiency in a low-latitude region like ours (Indore ∼22^{◦}) needs to be researched. Hence, we also use the NTCM-GL^{9} model to estimate slant iono delay. A comparison of all the methods is presented.

### 4.1 Nodal vertical delay estimates, their variances, and comparison with GAGAN estimates

NavIC transmits vertical delay estimates and their accuracy (expressed using GIVEI [Grid Ionospheric Vertical Error Indicator]) in message type 5 for L5 users through a set of 90 grid points (at 5^{◦} × 5^{◦} spacing) in the Indian region. A plot of vertical iono delay at grid point 1 and its 2*σ _{GIVE}* (95.5% error vertical bar) are shown in the top plot of Figure 9 for 24 hours. The maximum iono delay estimate occurs between 2:00 and 3:00 PM local time. As PRN1 (GSO) moves, its IPP viewed from Ahmedabad site moves also, resulting in a change of latitude-longitude of the grid point, portrayed in the bottom plot of Figure 9.

GAGAN’s vertical iono delay estimate at the same node 1 considered above, scaled by the factor , is shown in the top plot (in green) of Figure 9. We observe that GAGAN underpredicts the iono delay by about 3 m through all 24 hours, but this difference is within 2*σ* band of the NavIC grid-based iono delay estimates.

### 4.2 Grid-based ionospheric pierce point delay estimate and its variance

Iono delay at the IPP and its variance are calculated from the vertical grid delays (shown above) at the four surrounding nodes. Two algorithms are presented below.

#### 4.2.1 Using interpolation of grid iono delays

The vertical delay estimate at PRN1’s IPP using the well-known bilinear interpolation [Appendix D, Ref. 1] is shown in the left plot of Figure 10. With regard to the variance of the IPP delay estimate, the NavIC ICD specifies calculating the standard deviation from the four surrounding *σ _{GIVE,i}* using the same bilinear interpolation as done for the IPP vertical delay estimate:

10

This, however, appears oversimplistic because one should instead determine the variance of the IPP iono delay from the variances of the surrounding nodes using

11

assuming that these variances are uncorrelated, as 5^{◦} (grid points separation) latitude distance (555 km or 300 nautical miles) may be long enough for the iono delays to be uncorrelated. A comparison of the standard deviations, Equations (10) and (11), are shown in the left plot of Figure 10. The linearly interpolated sigma (blue circle bar) is found to be conservative compared to the sigma (magenta triangle bar) calculated from the variance of the nodes.

#### 4.2.2 Using standard deviation based weights

While determining the vertical iono delay at an IPP with bilinear interpolation, the variances of the nodal iono delay estimates are not considered. Intuitively, though, this consideration appears warranted since the nodal iono delay estimate with a larger variance should contribute commensurately lesser to the IPP delay. Lapucha and Huff^{8} provide one such interpolation scheme, in Equation (12), for multisite DGPS applied to ships. In Equation (12), the weights *W _{i}*, defined in Equation (13), are assigned to a grid point i;

*D*in Equation (13) denotes the distance between IPP and the grid point.

_{i}12 13

Both the distance *D _{i}* and the

*it*h node

*σ*, defined earlier, are in the denominator of the weight

_{GIVE,i}*W*so that the larger this product for a surrounding node, the smaller its iono delay contribution. The IPP vertical delay thus calculated is illustrated in Figure 10 by the plot on the right. In the present case, the IPP delay based on Equations (12) and (13) differs little from the delay calculated using bilinear interpolation and illustrated in the left plot of Figure 10 because the variances of the four grid points also differ little. The sigma calculated is also smaller (right plot of Figure 10) than the linearly interpolated sigma suggested in the ICD.

_{i}### 4.3 Slant iono delays using NavIC grid-based corrections, Klobuchar coefficients and NTCM-GL model

The slant iono delay is computed from the vertical IPP delay using the mapping function described in Appendix D of the ISRO^{1} and from Klobuchar coefficients (Message Type 11) using Appendix H (also ISRO).^{1} A comparison of these delays with an estimate using a model called NTCM-GL^{9} is shown in Figure 11. These results are commented on in the next section after the true iono delay calculated from the DF L5 and S1 pseudoranges is presented. The NTCM-GL F10.7 parameter value for our data set is 70, and the magnetic north pole latitude and longitude are +86.45^{◦} and +175.35^{◦}, respectively.

## 5 IONOSPHERIC DELAY ESTIMATES USING DUAL-FREQUENCY SIGNALS

While the iono delay presented in the previous section is supplied by the NavIC system for single-frequency users to remove the delay from the signal, the message could be old because the delay is updated every 288 seconds^{4} and the estimates have some specified variances. Dual-frequency users, on the other hand, can calculate the instantaneous iono delay in the code phase observables on the fly for greater navigation accuracy. So it is worthwhile to compare the two and that is done below. To reduce the noise in the DF code estimates, smoothing is performed with carrier phase measurements. The effect of ignoring multipath in DF pseudorange observables on the iono delay estimate is illustrated also.

### 5.1 Estimation with code phase

The measurement equation of L5 iono delay calculated using DF pseudoranges,^{12,13} with multipath included, is

14

which shows that the multipath error *M _{dƒ}* and the blended code phase noises at the two frequencies corrupt the iono delay estimates. Further, the noise in the estimate, given by , is enhanced. Assuming

*σ*to be the uncorrelated standard deviation of

*ν*, the standard deviation of the noise in the iono delay estimate for NavIC (L5-S1) and GPS (L1-L2) are 1.82

*σ*and 2.1

*σ*, respectively.

#### 5.1.1 Illustrations

##### Comparison with grid-based single-frequency estimates

The DF code phase estimate of iono delay is noisy and biased by ∼6 m due to the multipath error. To aid in comparison with the iono delay estimates with the single-frequency methods, the biased iono delay estimate is shifted down by 6 m, shown in yellow in Figure 11. This shifted down estimate lies within the 2*σ* limits of the single-frequency NavIC grid-based iono delay estimate (red) in the figure. This suggests that it may be preferable to use the single-frequency iono delay estimate, which is inherently free from multipath error, but further corroboration is warranted. The estimate using the Klobuchar coefficient, shown in green, follows the DF estimate well near noon hours (12-15 hours), but the error increases during quiet times. Also, during 8:00 to 10:00 hours and 17:00 to 20:00 hours, the Klobuchar estimate is outside the 2*σ* band from below. The NTCM-GL estimate, on the other hand, is within the 2*σ* limits except touching the lower limit during 16:00 to 17:00 hours. The peak value of NTCM-GL delay occurs at 13:00 hours and, for the Klobuchar model, at 14:00 hours, nearly the same as for the DF delay. The Klobuchar delay, however, falls off the peak more rapidly than the DF delay.

##### Bias multipath appearing as iono delay

The DF iono delay code estimate shown in Figure 11 is biased due to the multipath term in the right side of Equation (14). This is illustrated by comparing concurrent iono delays for two PRNs at the two locations in our campus and at two other sites in ISI, Kolkata (∼1300 km away) and SVNIT, Surat (∼300 km away) on 20 August 2018 in Figure 12. The spectrum similarity in these four iono delays is apparent, differing only by biases. Evidently the Workshop, IITI, has the highest multipath bias, and the other three have far lesser throughout the interval. We observe in Figure 12 that the iono delay at the Hub site is less than at the Workshop terrace by nearly a constant value, and the iono delays at ISI, Kolkata, differ from the iono delay at the Hub site by some time-varying bias that becomes negative around midnight and the morning hours. The iono delay at SVNIT, Surat, also shows negative iono delay similar to ISI, Kolkata, during midnight and morning hours. Code iono delays are, by definition, positive, and their negative estimate is caused by local unknown constant multipath errors with 180^{◦} phase error relative to the direct signal. This clearly indicates that a constant multipath error affects iono delay estimates for stationary receivers significantly. Estimation of this constant multipath error is beyond the scope of this paper, but will be pursued in the future.

### 5.2 Smoothing of the iono delay estimates with carrier phase observables and ambiguity resolution

The iono delay *I*_{l5}, calculated using DF carrier phase (CP) measurements,^{12,13} is precise but ambiguous as seen from its measurement equation:

15

Since the CP multipath errors are at centimeter level, these are ignored in Equation (15). The bias *B* on the right side of Equation (15) is an unknown real constant consisting of the difference between the L5 and S1 ambiguities as long as there are no cycle slips in either frequencies. Also, the noise in the synthesized carrier phase observable (*n*2), Equation (15), is two orders of magnitude less than *n*1 in the synthesized code observable Equation (14). So these two synthetic observables can be used to arrive at low-noise estimates of the two unknowns, iono delay *I*_{l5} and the real constant *B*, with a complementary filter shown in Figure 13 [Example 5.4, Ref. 13]. The two estimates, however, will be affected by the multipath error *M _{dƒ}* in Equation (14), modulated by the bandwidth of the low-pass filter because the filtering process does not incorporate the multipath error model explored in Section 3.

^{19}

The recursive equations of the complementary KF in Figure 13 to estimate *B* and *I*_{l5} are

16 17

Because the recursive index *k*, indicating the discrete measurement instants and also the inverse of the filter gain under white noise assumption, continues to rise, Montenbruck et al^{20} suggest limiting *k* to 50 for 1-second measurement intervals. This limit is indicative of the steady-state Kalman gain equal to the ratio carrier/code phase noise standard deviations, about 1/50.

#### 5.2.1 Illustration

Figure 14 shows the estimates of iono delay using the method presented above for two cases: (1) without removing the oscillatory multipath error from the pseudorange measurements, shown in the left column of Figure 14; (2) with oscillatory multipath removed from the pseudorange measurement shown in the right column of Figure 14. This side-by-side comparison clearly shows how the iono delay estimate is affected by the multipath error.

#### 5.2.2 Iono delay estimate with multipath error present

The iono delays estimated from the raw pseudoranges using DF code observables, Equation (14), and the smoothed DF observables, Equation (17), are shown in violet and blue colors, respectively, in the first plot on the left column of Figure 14. The iono delay calculated using smoothed DF observables are less noisy but follow the trend of the iono delay estimated from the raw DF code observables. The Accord receiver also estimates the iono delay through proprietary methods and conveys it through a set of text files. This estimate, shown in yellow in the same plot, is not visible as it is buried under the smoothed DF results.

The *B _{k}* estimate, Equation (16), shown in the second plot of the left column, is supposed to be a constant in the absence of cycle slips. But presently, it varies between −11.5 m and −5.8 m, neglecting the initial transients. This variation is due to the multipath errors in pseudorange observables as shown in Equation (14). We verify this by using the estimates of and from Equation (8) and similar equation for S1 and multiplying these with the frequency term, . This is the magnitude by which the multipath errors in raw pseudorange observables are enhanced by the DF code combinations as observed in the third plot in the left column of Figure 14. As expected, the variation pattern of the multipath error with time in the third plot and those of the estimates of the DF iono delay and the constant

*B*in the top two resembles one another. Thus, if oscillatory multipath errors are ignored, the estimate of DF code-based iono delay can vary from 15 m to 10 m, a variation of 5 m in 30 minutes as seen in the top left plot of Figure 14 during 22:00 to 22:30 local time.

The residual of *B _{k}*, that is, , is shown in the last plot of the left column along with the 3

*σ*bounds; the sigma is calculated by substituting

*σ*

_{νl5}= 18.3 cm and

*σ*

_{νs1}= 29.3 cm in Equation (14). The statistical values of the noise are chosen from the conclusions drawn from Section 3.

#### 5.2.3 Iono delay estimates with multipath oscillations removed

The iono delay estimates calculated using the DF code, Equation (14), and the smoothed DF observables, Equation (17), using oscillatory multipath error-removed pseudoranges (Section 3.3), are shown in violet and blue colors in the first plot on the right column of Figure 14. The smoothed iono delay estimate has less variation, ∼2 m in the 30 minutes, compared to 5 m of iono delay seen in the top left plot using the raw pseudoranges. During the 8-hour duration, the difference in the iono delay peak-to-peak value is less than 3 m compared to ∼7 m in the top left plot. This is due to the fact that *B _{k}* from the oscillatory multipath compensated pseudoranges has less oscillations than those of the raw pseudoranges; this is verified by observing

*B*in the second row on the right plot. Also, the left column results have both constant and oscillatory multipath errors, whereas the right column results have just the unknown constant multipath error. The

_{k}*M*term is suspected to be the reason for the high iono delay of ∼12 m during nighttime in the top plot (right column) in Figure 14. The NavIC antenna is mounted on the terrace of the Workshop building with other dish antennas nearby, so the neighboring metallic structures are suspected to have enhanced

_{bias}*M*. When the antenna is moved to the Hub building (∼200 m away) with no metallic surfaces nearby, the nighttime iono delay value is around 1 to 2 m as shown in Figure 12.

_{bias}Even after compensating for the oscillatory multipath errors, the resultant iono delay estimate has some residual multipath error. This is evident by observing the third plot in the right column of Figure 14 which shows the enhancement of the remaining multipath error and noise on the DF iono estimation as explained earlier for the raw pseudoranges. This plot matches well with the DF code based iono estimate (violet) shown in the first plot on the right column, thus clearly showing the inclusion of residual errors in the iono estimate. The bottom plot shows the measurement residual of *B _{k}* along with its 3

*σ*bound calculated by substituting

*σ*

_{νl5,comp}= 25.9 cm and

*σ*

_{νs1,comp}= 41.4 cm in Equation (14). The standard deviation of the oscillatory multipath compensated pseudorange noise is times that of the raw pseudorange, as can be inferred from Equation (9).

Iliopoulos et al^{21} present a method to jointly estimate the multipath error and the iono delay. However, this can only be implemented by advanced tracking in the receiver. This work will be pursued in the future.

## 6 TROPOSPHERIC DELAY ESTIMATION

Signal delay caused by the tropospheric region is estimated using the Collins model.^{22} In this method, the slant tropo delay *T*_{slant} needed to compensate for range error is calculated from vertical tropo delay *T*_{vertical} using the elevation *El* based mapping function (Ref. 22) of the receiver-satellite vector using Equation (18). Slant tropo-delay for PRN2 is composed of dry and wet components contributing 2.3 to 4.7 m and 0.2 m, respectively. A plot of the slant tropo-delay errors is shown in Figure 15 for all NavIC satellites.

18

## 7 LINE-OF-SIGHT COMPONENT OF EPHEMERIS ERROR AND ITS VARIANCES

In this section, we present two techniques to estimate line-of-sight (LoS) components of ephemeris errors of the navigation satellites. In Section 7.1, we present a method to calculate the ephemeris error along the line-of-sight vector and the satellite clock error from the orbital parameters at update epochs in the navigation message. In Section 7.2, we consider the smoothed iono-free pseudorange measurement, remove the estimated receiver-satellite range, remove the satellite and the receiver clock biases, and tropo delay, and pass the remainder separately through both a low-pass filter and a complementary high-pass filter. The low-pass filter outputs the LoS component of the ephemeris error plus any unknown bias, and the high-pass filter outputs the remaining zero-mean high-frequency noise the variance of which is then determined.

### 7.1 Ephemeris error LOS components at orbit update epochs

The radial (R), transverse (T), and normal (N) components of ephemeris error, using the satellite’s new and old radial locations, longitudes, and orbit inclinations, in the orbital plane at instant *k* (when new orbital parameters arrive) are

19

where *r _{k}* is the corrected radius of the satellite calculated as shown in Table C.2 of Farrell.

^{13}The superscripts

*new*implies

*r*is calculated using the new orbital elements available at

_{k}*k*and

*old*implies that

*r*is calculated using the old elements (

_{k}*k*− 7200, if orbits are updated every 2 hours) propagated to the

*k*th instant.

The transverse component of the ephemeris error is calculated using Ω* _{i}*, the inertial longitude of the satellite.

20

where Ω_{0} is the Right Ascension of the Ascending Node (RAAN) at the reference time when orbital parameters are available, is the rate of change of RAAN, *t _{k}* is the difference between

*k*and reference time, and

*u*is the argument of latitude of the satellite. The normal component of the ephemeris error is calculated using the inclination of the orbital plane. is the new inclination angle available at

_{k}*k*, and is the inclination at

*k*propagated using the orbital elements at

*k*− 7200. The ephemeris error vector in the RTN frame is transformed to the ECEF frame using the PQW matrix defined below [Sec. 2.2.3, Ref. 23]:

21

where *Ω _{k}* is the longitude of the ascending node from the Greenwich meridian,

*u*is the argument of latitude, and

_{k}*i*is the inclination of the orbit. This vector when projected onto the line-of-sight unit vector through its azimuth and elevation angles to the satellite gives us the estimate of ephemeris error in the pseudorange observable.

_{k}The error in the satellite clock bias is calculated from the just updated clock correction and the previous clock correction at the same instant; thus,

22

where Δ*t _{s}* is calculated using Equation (3) and Equation (5).

#### 7.1.1 Illustration

The radial, transverse, and normal components of the ephemeris error calculated are shown in the top three plots of Figure 16 for PRN7 for four discrete events of orbit updates over the duration of some 8 hours.

The radial error is less than about 8 cm while the tangential error is mostly less than 2 m. The normal error component is larger, about 28 m on one occasion. The orbital update period is 2 hours. The maximum radial, transverse, and normal errors reported in Ganeshan^{24} are 4 m, 26 m, and 16 m, respectively. The line-of-sight projection of these vectors *E* shown in the fourth plot varies between −1.3 and 1.5 m. Montenbruck et al^{3} reports a monthly root-mean-square line-of-sight error of 7.6 m and 9.1 m for PRN1 and PRN2, respectively. It is to be noted that these discrete estimates of *E* in the fourth plot are without any multipath bias. The satellite clock bias errors at the arrival of new clock parameters are shown in the last plot of Figure 16. It varies between −0.5 m and 0.4 m for PRN7 which is smaller than the ephemeris error contribution.

### 7.2 Estimation of line-of-sight component of ephemeris error plus bias and its variance

Now that we have a reference value of ephemeris error, we move to estimating the ephemeris error in the pseudorange observable from which the oscillatory multipath and iono delay has been removed. The pseudorange is compensated for the instantaneous estimated range , receiver and satellite clock bias, and Collins-modeled tropospheric delay *T*, to arrive at a signal composed of the ephemeris error component *E* along the line-of-sight, bias due to the multipath and range and clock bias estimation errors and and random noise *ζ* shown below.

23

Observing the above equation, we see that the slowly varying multipath errors and the time-varying range and receiver clock biases will influence the ephemeris error estimate. Because the satellites orbits are updated no sooner than every 15 minutes (about 1-mHz frequency), the ephemeris errors of the satellites grow between updates due to the low-frequency unaccounted satellite accelerations and drop to lower values at each update. This leads one to infer that the ephemeris errors would have some low-frequency spectrum. The NavIC pseudorange observable, on the other hand, is available at 1-second intervals, so the random noise will be of wide spectrum up to the frequency 0.5 Hz. This spectrum dissimilarity suggests that the low-frequency ephemeris error *E* summed with multipath bias and any other unidentified biases could be extracted using a low-pass filter, whereas the random noise could be filtered out with a complementary high-pass filter.

Hence, we design a low-pass filter whose input is the left side of Equation (23) to extract the low-frequency error *E* summed with biases and a complementary high-pass filter to extract the noise. Angular frequency of the ephemeris update at 15-minute intervals is . The low-pass filter with its bandwidth *ω _{lpƒ}* set to 20 ×

*ω*(time period 45 s) is chosen to minimize the phase delay through the filter. Thus, we select the low-pass filter transfer function:

_{ephm}24

The corresponding digital filter is arrived at by substituting (bilinear transform) in *H _{lpƒ}* (

*s*), where

*T*is the sample period.

_{s}The inverse z-transform provides the following difference equation in the time domain, necessary for implementation in a computer:

25

where *V _{i}, V_{o}* are the input and output signal and

*k*represents a time sample. The input signal is the left side of Equation (23) as

*T*= 1 second.

_{s}#### 7.2.1 Illustrations

The estimated ephemeris error components along the line-of-sight plus biases for geosynchronous (GSO) and geostationary (GEO) PRNs are shown in the first and second plots of Figure 17, respectively. The ephemeris LOS error estimate in both plots has a negative bias of 3 to 6 m. The User Range Accuracy (URA, Section 6.2.1.4 in Ref. 1) of each satellite, which provides 1*σ* error bound for the Signal-in-Space range error, is shown in the top two plots. The ephemeris error estimates for all PRNs lie at the edge of this error bound and sometimes exceed it (for example, PRN4 in the top plot between 0 and 2 hours).

This higher error estimate is due to the bias and residual errors in range and clock bias estimate seen in Equation (23). The pseudorange measurements used in this section were collected on the terrace of the Workshop building which had significant multipath bias error as shown in Figure 12. We expect the algorithm to provide a reasonable estimate of *E* if multipath bias errors are negligible.

The orbital/clock parameters are nominally updated once every 2 hours. But there are provisions for faster updates during station keeping operations.^{25} The issue of new parameters is conveyed through Issue of Data Ephemeris and Clock (IODEC) which ranges between 0 and 255. IODEC variation between 160 and 255 signifies the frequent update mode (Table 22, Ref. 1). PRN6 parameters are updated regularly (every 900 s) as seen from the IODEC values (225-255) in the bottom plot of Figure 17. This is likely due to the drift in PRN6’s onboard clock. All other satellites are updated once every 2 hours during this study.

The ephemeris LOS-component error plus some inevitable and unknown multipath biases for PRN7, shown in the middle plot of Figure 17, is compared with the unbiased ephemeris LOS error estimate shown in the fourth plot of Figure 16. This comparison in Figure 18 clearly shows that the bias in the ephemeris estimate (blue) of 6 m exceeds the 1*σ* = 2.5 m limits (green bar) of the unbiased ephemeris estimate (red marker) at the update epochs. The contribution of the satellite clock bias error and the ephemeris LOS error (green marker) also lies within 1*σ* for the 8-hour duration that we studied.

### 7.3 Signal noise variance estimation

The difference equation of the high-pass filter to extract the random noise in the pseudorange is

26

The mean and standard deviation of the residual noise at the filter’s output are 0.05 cm and 0.48 cm, respectively. These values indicate the high precision of the smoothed observable .

## 8 CONCLUSIONS

All known errors affecting NavIC pseudorange observables are investigated in the paper, and our findings are summarized below. Two-hour updates of the NavIC satellites clock error coefficients are illustrated. It is shown that at update epochs, the clock coefficients correction can cause a step change in range by as much as 2 m. Corrections for satellite instrument delays, expressed in terms of total group delay for single frequency users, are illustrated. There appears to be an inadvertent error in the output of the Accord’s receiver in the calculation of this correction for single frequency users. This error is illustrated. Since NavIC satellites have small intentional eccentricity, amplitude of sinusoidal relativistic corrections in 24 hours can be as high as 6 nanoseconds.

The multipath oscillation frequency due to a ground reflector is 0 to 0.1 mHz for NavIC signals from a GEO satellite and 0 to 1.6 mHz (16-fold) for a GSO satellite. The amplitude of multipath oscillations is found to be as high as 3 m peak-to-peak at times with a period of ∼50 minutes but are mostly removed by utilizing their sidereal repeatability. Also, the multipath errors in L5 and S1 pseudoranges are different; hence, it is important to correct the error in each signal separately before using it for iono-delay or position estimation. For real-time removal of multipath errors, spectrum similarity of multipath error with C/No of the code phase observable is demonstrated.

For the days and sites considered in this study, the maximum iono delay of the NavIC signals at the surrounding nodes during the day is found to be 8 m with *σ* = 1.4 to 1.8 m. The iono delay of the GPS signal through the GAGAN system is found to be 2 m less than the NavIC grid delay but within its 1 to 2 *σ*. The slant iono delay calculated using Klobuchar coefficients and NTCM-GL mostly lie within 2*σ* of the grid-based slant delay with the Klobuchar-based estimates deviating more often. The multipath oscillations in the iono delay smoothed by dual-frequency phase observables are seen to be as high as 5 to 6 m peak-to-peak. When these oscillations are removed using sidereal repeatability, the blended, dual-frequency multipath bias, ∼12-13 m, masquerades as an iono delay estimate during night hours. While the Klobuchar model predicts iono delay at night ∼2 m, the estimated multipath oscillations-free iono delay at four different sites are found to vary from −3 m to +10 m in the same hours. Clearly, the effect of multipath oscillations and bias on iono delay mitigation and subsequently on position estimation will be significant and cannot be ignored. The tropo delay for all PRNs lies within 2 to 6 m with the dry component contributing the most.

The ephemeris error vector projection on the line-of-sight vector for PRN7 (a GEO at 129.5^{◦} E) using the step corrections in orbital parameters is found to be within −1.3 m and +1.5 m for the 8-hour duration. The satellite clock bias errors are less than 0.6 m. The sum of the line-of-sight components of the orbit and the satellite clock errors is well within the 1*σ* = 2.5 m broadcast ephemeris accuracy conveyed by URA. The same ephemeris LOS error calculated using low-pass filter is −6 m, exceeding the 2*σ* bound. This is caused by the multipath bias in the smoothed pseudorange. A real-time estimation of multipath error invoking its relationship with signal-to-noise ratio following Sleewaegen^{7} and Breivik et al^{18} will be pursued in the future.

## HOW TO CITE THIS ARTICLE

Azeez A, Hablani H. Assessment of errors in NavIC observables for stationary receivers. *NAVIGATION-US*. 2020;67:347–364. https://doi.org/10.1002/navi.355

## ACKNOWLEDGEMENTS

We are grateful to SAC/ISRO for providing the Accord receivers and also for funding this study. Especially, we would like to thank Shri K.S. Parikh, Shri N.M. Desai, Shri. Atul Shukla, Shir, Rajat Acharya, Shri. Ashish Shukla, Shri. V.K. Tank, and Shri. Nishkam Jain of SAC, ISRO, for their support and guidance. We are also grateful to Mr. Mehul Desai, SVNIT, for sharing Surat’s iono delay measurements. We would like to thank Dr Surendra for providing GAGAN vertical ionospheric data. Our gratitude also goes to Dr Saurabh Das, IITI for sharing the iono delays from Kolkata. We would like to thank Shri. Subramanya Ganesh and Shri. Suresh Dakkumalla of ISTRAC, ISRO, for clarifying our questions on NavIC orbit updates.

## Footnotes

**Funding information**Indian Space Research Organisation, Grant/Award Number: NGP-17

- Received June 20, 2019.
- Revision received November 5, 2019.
- Accepted December 16, 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.