**Part 3**

## **GNSS Errors Mitigation and Modelling**

356 Global Navigation Satellite Systems – Signal, Theory and Applications

Parra, S. & Angel, J. (2005). *Low cost navigation system for UAV's*, Aerospace Science and

Cramera, M. (2006). *The ADS40 Vaihingen/Enz geometric performance test*, ISPRS Journal of Photogrammetry and Remote Sensing, Volume 60, Issue 6, pp. 363-374 Kunii, Y. & Chikatsu, H. (2001). *Application of 3-Million Pixel Armature Camera for the 3D* 

Shapiro, R. (1978). *Direct linear transformation method for three-dimensional cinematography*, Res.

Kumagai, H., Kubo, Y., Kihara, M. & Sugimoto, S. (2002). *DGPS/INS/VMS Integration for* 

Kuamgai, H., Kindo, T., Kubo, Y. & Sugimoto, S. (2000). *DGPS/INS/VMS Integration for High* 

Takagi, M. & Shimoda, H. (2004). *Handbook of image analysis*, University of Tokyo Press Chen, T., Shibasaki, R. & Murai, S. (2003). *Development and Calibration of the Airborne Three-*

Photogrammetry and Remote Sensing PE&RS, vol.69, No.1,pp. 71-78

Photogrammetry and Remote Sensing, vol.41, no.4 pp. 77-84

*Modeling of Historical Structures*, Asian Journal of GEOINFORMATICS Vol. 2, No. 1,

*High Accuracy Land-Vehicle Positioning*, Journal of the Japan Society of

*Accuracy Land-Vehicle Positioning*, Proceedings of the Institute of Navigation, GPS-

*Line Scanner (TLS) Imaging System*, Journal of the American Society for

Technology, Issue 6, Volume 9, pp.504-516

pp. 39-48

Quart. 49, pp. 197-205

2000, Salt Lake

**15** 

*China* 

Shuanggen Jin

**– Methods and Results** 

**GNSS Atmospheric and Ionospheric Sounding** 

The GPS atmospheric and ionospheric delays have been considered as an error source for a long time. In 1992 when the GPS became fully operational, Ware (1992) suggested limb sounding the Earth atmosphere using GPS atmospheric delay signals. In April 1995, the small research satellite of Microlab-1 was successfully put into a Low Earth Orbit (LEO) to validate the GPS radio occultation method (Feng and Herman, 1999). Since then, the GPS/Meteorology Mission (GPS/MET) has been widely used to produce accurate, all weather pressure, temperature, density profiles in the troposphere and the ionospheric total electron content (TEC) as well as electron density profiles (Rocken, 1997; Hajj and Romans, 1998; Syndergaard, 2000), to improve weather analysis and forecasting, monitor climate change, and monitor ionospheric events. While traditional observing instruments, e.g. water vapour radiometer (WVR), incoherent scatter radars (ISR), ionosonde, topside sounders onboard satellites, in situ rocket and satellite observations, are expensive and also partly restricted to either the bottomside ionosphere or the lower part of the topside ionosphere (usually lower than 800 km), such as ground based radar ionospheric measurements. While GPS satellites in high altitude orbits (~20,200 km) are capable of providing details on the structure of the entire atmosphere, even the plasma-sphere. Moreover, GPS is a low-cost, allweather, near real time, and high-temporal resolution (1~30s) technique. Therefore, GPS is a powerful tool to sound the atmosphere and ionosphere as well as their application in

The tropospheric delay of GPS signal through the neutral atmosphere was one of major error sources in navigation and positioning, which contributes a bias in height of several centimetres (Tregoning et al. 1998). Nowadays, GPS has been used to determine the zenith tropospheric delay (ZTD) (Jin and Park, 2005) through mapping functions (Niell, 1996). The ZTD is the integrated refractivity along a vertical path through the neutral

> 6 0 *ZTD c N s ds* τ

10 ( ) <sup>∞</sup> <sup>−</sup> = = ∫ (1)

**1. Introduction** 

meteorology, climate and space weather.

**2. Tropospheric sounding** 

atmosphere:

*Shanghai Astronomical Observatory, Chinese Academy of Sciences* 

## **GNSS Atmospheric and Ionospheric Sounding – Methods and Results**

Shuanggen Jin

*Shanghai Astronomical Observatory, Chinese Academy of Sciences China* 

## **1. Introduction**

The GPS atmospheric and ionospheric delays have been considered as an error source for a long time. In 1992 when the GPS became fully operational, Ware (1992) suggested limb sounding the Earth atmosphere using GPS atmospheric delay signals. In April 1995, the small research satellite of Microlab-1 was successfully put into a Low Earth Orbit (LEO) to validate the GPS radio occultation method (Feng and Herman, 1999). Since then, the GPS/Meteorology Mission (GPS/MET) has been widely used to produce accurate, all weather pressure, temperature, density profiles in the troposphere and the ionospheric total electron content (TEC) as well as electron density profiles (Rocken, 1997; Hajj and Romans, 1998; Syndergaard, 2000), to improve weather analysis and forecasting, monitor climate change, and monitor ionospheric events. While traditional observing instruments, e.g. water vapour radiometer (WVR), incoherent scatter radars (ISR), ionosonde, topside sounders onboard satellites, in situ rocket and satellite observations, are expensive and also partly restricted to either the bottomside ionosphere or the lower part of the topside ionosphere (usually lower than 800 km), such as ground based radar ionospheric measurements. While GPS satellites in high altitude orbits (~20,200 km) are capable of providing details on the structure of the entire atmosphere, even the plasma-sphere. Moreover, GPS is a low-cost, allweather, near real time, and high-temporal resolution (1~30s) technique. Therefore, GPS is a powerful tool to sound the atmosphere and ionosphere as well as their application in meteorology, climate and space weather.

## **2. Tropospheric sounding**

The tropospheric delay of GPS signal through the neutral atmosphere was one of major error sources in navigation and positioning, which contributes a bias in height of several centimetres (Tregoning et al. 1998). Nowadays, GPS has been used to determine the zenith tropospheric delay (ZTD) (Jin and Park, 2005) through mapping functions (Niell, 1996). The ZTD is the integrated refractivity along a vertical path through the neutral atmosphere:

$$ZTD = c\tau = 10^{-6} \int\_0^{\circ} \text{N}(s)ds\tag{1}$$

GNSS Atmospheric and Ionospheric Sounding – Methods and Results 361

0 60 120 180 240 300 360

y(t)=1.071\*104-4.094\*t+65.090\*sin(2\*π\*t+1.013)+0.224\*sin(π\*t+0.385)

Longitude (deg)

1998 1999 2000 2001 2002 2003 2004 2005 2006

Time (year)

The processing software must resolve or model the orbital parameters of the satellites, solve for the transmitter and receiver positions, account for ionospheric delays, solve for phase

Fig. 2. ZTD time series at TOW2 station, Australia. The solid line is the fitting results,

consisting of a linear decrease and seasonal components

**2.2 Zenith tropospheric delay retrieval** 


2300

2350

2400

2450

2500

2550

ZTD (mm)

2600

2650

2700

2750

2800

Fig. 1. The distribution of global IGS GPS sites



0

Latitude (deg)

30

60

90

where c is the speed of light in a vacuum, *τ* is the delay measured in units of time and N is the neutral atmospheric refractivity. The N is empirically related to standard meteorological variables as (Davis et al. 1985)

$$N = k\_1 \rho + k\_2 \frac{P\_w}{Z\_w T} + k\_3 \frac{P\_w}{Z\_w T^2} \tag{2}$$

where ( 1,2,3) *<sup>i</sup> k i* = is constant, ρ is the total mass density of the atmosphere, *Pw* is the partial pressure of water vapor, *Zw* is a compressibility factor near unity accounting for the small departures of moist air from an ideal gas, and T is the temperature in degrees Kelvin. The integral of the first term of equation (2) is the hydrostatic component ( ) *Nh* and the integral of the remaining two terms is the wet component (*Nw* ) . Thus, ZTD is the sum of the hydrostatic or dry delay (ZHD) and non-hydrostatic or wet delay (ZWD), respectively. The dry component ZHD is related to the atmospheric pressure at the surface, while the wet component ZWD can be transformed into the precipitable water vapor (PWV) and plays an important role in energy transfer and in the formation of clouds via latent heat, thereby directly or indirectly influencing numerical weather prediction (NWP) model variables (Bevis et al, 1994; Tregoning et al. 1998; Manuel et al., 2001). Therefore, the Zenith Tropopospheric Delay (ZTD) is an important parameter of the atmosphere, which reflects the weather and climate processes, variations, and atmospheric vertical motions, etc.

In the last decade, ground-based GPS receivers have been developed as all-weather, high spatial-temporal resolution and low-cost remote sensing systems of the atmosphere (Bevis et al., 1994; Manuel et al., 2001), as compared to conventional techniques such as satellite radiometer sounding, ground-based microwave radiometer, and radiosondes (Westwater, 1993). With independent data from other instruments, in particular water vapor radiometers, it has been demonstrated that the total zenith tropospheric delay or integrated water vapor can be retrieved using ground based GPS observations at the same level of accuracy as radiosondes and microwave radiometers (Elgered et al. 1997; Tregoning et al. 1998). Currently, the International GPS Service (IGS) has operated a global dual frequency GPS receiver observation network with more than 350 permanent GPS sites since 1994 (Beutler et al., 1999). It provides a high and wide range of scales in space and time to study seasonal and secular variations of ZTD as well as its possible climate processes.

#### **2.1 GPS data and analysis**

The IGS (International GPS Service) was formally established in 1993 by the International Association of Geodesy (IAG), and began routine operations on January 1, 1994 (Beutler et al. 1999). The IGS has developed a worldwide network of permanent tracking stations with more than 350 GPS sites, and each equipped with a GPS receiver, providing raw GPS orbit and tracking data as a data format called Receiver Independent Exchange (RINEX). All available near-real-time global IGS observation data archived in the Global IGS Data Center, which contributes to geodesy and atmosphere research activities in a global scale. Here the globally distributed 150 IGS sites with better continuous observations are selected with spanning at least four years of measurements (Figure 1), and most sites observations are from 1994 to 2006.

360 Global Navigation Satellite Systems – Signal, Theory and Applications

where c is the speed of light in a vacuum, *τ* is the delay measured in units of time and N is the neutral atmospheric refractivity. The N is empirically related to standard meteorological

*P P Nk k k*

partial pressure of water vapor, *Zw* is a compressibility factor near unity accounting for the small departures of moist air from an ideal gas, and T is the temperature in degrees Kelvin. The integral of the first term of equation (2) is the hydrostatic component ( ) *Nh* and the integral of the remaining two terms is the wet component (*Nw* ) . Thus, ZTD is the sum of the hydrostatic or dry delay (ZHD) and non-hydrostatic or wet delay (ZWD), respectively. The dry component ZHD is related to the atmospheric pressure at the surface, while the wet component ZWD can be transformed into the precipitable water vapor (PWV) and plays an important role in energy transfer and in the formation of clouds via latent heat, thereby directly or indirectly influencing numerical weather prediction (NWP) model variables (Bevis et al, 1994; Tregoning et al. 1998; Manuel et al., 2001). Therefore, the Zenith Tropopospheric Delay (ZTD) is an important parameter of the atmosphere, which reflects the weather and climate processes, variations, and

In the last decade, ground-based GPS receivers have been developed as all-weather, high spatial-temporal resolution and low-cost remote sensing systems of the atmosphere (Bevis et al., 1994; Manuel et al., 2001), as compared to conventional techniques such as satellite radiometer sounding, ground-based microwave radiometer, and radiosondes (Westwater, 1993). With independent data from other instruments, in particular water vapor radiometers, it has been demonstrated that the total zenith tropospheric delay or integrated water vapor can be retrieved using ground based GPS observations at the same level of accuracy as radiosondes and microwave radiometers (Elgered et al. 1997; Tregoning et al. 1998). Currently, the International GPS Service (IGS) has operated a global dual frequency GPS receiver observation network with more than 350 permanent GPS sites since 1994 (Beutler et al., 1999). It provides a high and wide range of scales in space and time to study seasonal and secular variations of ZTD as well as its possible

The IGS (International GPS Service) was formally established in 1993 by the International Association of Geodesy (IAG), and began routine operations on January 1, 1994 (Beutler et al. 1999). The IGS has developed a worldwide network of permanent tracking stations with more than 350 GPS sites, and each equipped with a GPS receiver, providing raw GPS orbit and tracking data as a data format called Receiver Independent Exchange (RINEX). All available near-real-time global IGS observation data archived in the Global IGS Data Center, which contributes to geodesy and atmosphere research activities in a global scale. Here the globally distributed 150 IGS sites with better continuous observations are selected with spanning at least four years of measurements (Figure 1), and most sites observations are

=+ + ρ

ρ

12 3 2 *w w w w*

*Z T Z T*

is the total mass density of the atmosphere, *Pw* is the

(2)

variables as (Davis et al. 1985)

where ( 1,2,3) *<sup>i</sup> k i* = is constant,

atmospheric vertical motions, etc.

climate processes.

from 1994 to 2006.

**2.1 GPS data and analysis** 

Fig. 1. The distribution of global IGS GPS sites

Fig. 2. ZTD time series at TOW2 station, Australia. The solid line is the fitting results, consisting of a linear decrease and seasonal components

#### **2.2 Zenith tropospheric delay retrieval**

The processing software must resolve or model the orbital parameters of the satellites, solve for the transmitter and receiver positions, account for ionospheric delays, solve for phase

GNSS Atmospheric and Ionospheric Sounding – Methods and Results 363

<sup>10</sup> log 5- 15.5

where *P* is the pressure in Pascals and *h* is the height in millimeters. Based on the Eq. (4), ZHD can be expressed as (5 /15.5) 2.28 \* 10 <sup>−</sup>*<sup>h</sup>* . Therefore, ZTD at all GPS sites can be

where the units of *ZTD* and *h* are in millimeters, respectively. Comparing GPS-derived ZTD

0 60 120 180 240 300 360

The GPS ZTD time series have been analyzed for 4-12 years at globally distributed 150 GPS sites. Using the least square the fitting parameters of all GPS site are obtained, including trend and seasonal variation terms (Jin et al., 2007). The mean secular ZTD variation trend is about 1.5±0.001 mm/yr. Figure 4 shows the distribution of the secular ZTD variation trends at all GPS sites as the yearly increase or decrease. It can be seen that the trends are positive in most parts of the Northern Hemisphere and negative in most parts of the Southern Hemisphere (excluding positive in Antarctic), corresponding to a systematic increase or decrease of ZTD. It is interesting to note that the downtrend in Australia is larger than other regions. This downtrend of ZTD is probably due to the highly deserted in Australia. In addition, the ZTD variation trend decreases with increasing altitude, and furthermore, the ZTD trends are almost symmetrical with altitude. This indicates that the secular ZTD variations are larger at the lower altitude and at the higher altitude the secular ZTD variations hardly increase or decrease. In addition, the sum of downward and upward

Longitude (deg)

with the empirical formula estimations, it has shown a good consistency.

approximately expressed as:


**2.4 Trend analysis** 

Fig. 3. Distribution of mean ZTD at global IGS sites



0

Latitude (deg)

30

60

90

*<sup>P</sup>* <sup>≈</sup> *<sup>h</sup>* (5)

1500

1720

1940

2160

2380

2600 mm

(5 /15.5) 2.28 \* 10 *<sup>h</sup> ZTD* <sup>−</sup> = (6)

cycle ambiguities and the clock drifts in addition to solving for the tropospheric delay parameters of interest. We use the GAMIT software (King and Bock 1999), which solves for the ZTD and other parameters using a constrained batch least squares inversion procedure. In addition, this study uses the newly recommended strategies (Byun et al. 2005) to calculate ZTD time series with temporal resolution of 2 hours from 1994 to 2006. The GAMIT software parameterizes ZTD as a stochastic variation from the Saastamoinen model (Saastamoinen 1972), with piecewise linear interpolation in between solution epochs. GAMIT is very flexible in that it allows a priori constraints of varying degrees of uncertainty. The variation from the hydrostatic delay is constrained to be a Gauss-Markov process with a specified power density of 2 cm / *hour* , referred to below as the "zenith tropospheric parameter constraint". We designed a 12-hour sliding window strategy in order to process the shortest data segment possible without degrading the accuracy of ZTD estimates. The ZTD estimates are extracted from the middle 4 hours of the window and then move the window forward by 4 hours. Finally, the ZTD time series from 1994 to 2006 are obtained at globally distributed 150 IGS sites with temporal resolution of 2 hours. For example, Figure 2 shows the times series of zenith total delay (ZTD) (upper) at TOW2 station, Australia.

#### **2.3 Global mean zenith tropospheric delay**

The ZTD consists of the hydrostatic delay (ZHD) and wet delay (ZWD). The ZHD can be well calculated from surface meteorological data, ranging 1.5-2.6 meters, which accounts for 90% of ZTD. It derives from the relationship with hydrostatic equilibrium approximation for the atmosphere. Under hydrostatic equilibrium, the change in pressure with height is related to total density at the height h above the mean sea level by

$$d\mathfrak{p} = -\rho(\mathfrak{h})\mathfrak{g}(\mathfrak{h})d\mathfrak{h} \tag{3}$$

where ρ( ) and *h* g(h) are the density and gravity at the height h. It can be further deduced as

$$\text{ZHD=kpc0}\tag{4}$$

where k is constant (2.28 mm/hPa) and p0 is the pressure at height h0 (Davis et al. 1985). It shows that the ZHD is proportional to the atmospheric pressure at the site. The ZWD is highly variable due possibly to varying climate, relating to the temperature and water vapour. The mean ZTD values at all GPS sites are shown in Figure 3 as a color map. It has noted that lower ZTD values are found at the areas of Tibet (Asia), Andes Mountain (South America), Northeast Pacific and higher latitudes (Antarctica and Arctic), and the higher ZTD values are concentrated at the areas of middle-low latitudes. In addition, the ZTD values decrease with increasing altitude, which is due to the atmospheric pressure variations with the height increase. Atmospheric pressure is the pressure above any area in the Earth's atmosphere caused by the weight of air. Air masses are affected by the general atmospheric pressure within the mass, creating areas of high pressure (anti-cyclones) and low pressure (depressions). Low pressure areas have less atmospheric mass above their locations, whereas high pressure areas have more atmospheric mass above their locations. As elevation increases, there are exponentially, fewer and fewer air molecules. Therefore, atmospheric pressure decreases with increasing altitude at a decreasing rate. The following relationship is a first-order approximation to the height (http://www.chemistrydaily.com/ chemistry/Atmospheric\_pressure):

$$
\log\_{10} P \approx 5 \text{-} \frac{h}{15.5} \tag{5}
$$

where *P* is the pressure in Pascals and *h* is the height in millimeters. Based on the Eq. (4), ZHD can be expressed as (5 /15.5) 2.28 \* 10 <sup>−</sup>*<sup>h</sup>* . Therefore, ZTD at all GPS sites can be approximately expressed as:

$$\text{ZTD} = 2.28 \, \text{\*} \, 10^{(5-h/15.5)} \tag{6}$$

where the units of *ZTD* and *h* are in millimeters, respectively. Comparing GPS-derived ZTD with the empirical formula estimations, it has shown a good consistency.

Fig. 3. Distribution of mean ZTD at global IGS sites

#### **2.4 Trend analysis**

362 Global Navigation Satellite Systems – Signal, Theory and Applications

cycle ambiguities and the clock drifts in addition to solving for the tropospheric delay parameters of interest. We use the GAMIT software (King and Bock 1999), which solves for the ZTD and other parameters using a constrained batch least squares inversion procedure. In addition, this study uses the newly recommended strategies (Byun et al. 2005) to calculate ZTD time series with temporal resolution of 2 hours from 1994 to 2006. The GAMIT software parameterizes ZTD as a stochastic variation from the Saastamoinen model (Saastamoinen 1972), with piecewise linear interpolation in between solution epochs. GAMIT is very flexible in that it allows a priori constraints of varying degrees of uncertainty. The variation from the hydrostatic delay is constrained to be a Gauss-Markov process with a specified power density of 2 cm / *hour* , referred to below as the "zenith tropospheric parameter constraint". We designed a 12-hour sliding window strategy in order to process the shortest data segment possible without degrading the accuracy of ZTD estimates. The ZTD estimates are extracted from the middle 4 hours of the window and then move the window forward by 4 hours. Finally, the ZTD time series from 1994 to 2006 are obtained at globally distributed 150 IGS sites with temporal resolution of 2 hours. For example, Figure 2 shows

the times series of zenith total delay (ZTD) (upper) at TOW2 station, Australia.

related to total density at the height h above the mean sea level by

The ZTD consists of the hydrostatic delay (ZHD) and wet delay (ZWD). The ZHD can be well calculated from surface meteorological data, ranging 1.5-2.6 meters, which accounts for 90% of ZTD. It derives from the relationship with hydrostatic equilibrium approximation for the atmosphere. Under hydrostatic equilibrium, the change in pressure with height is

> *dp h g h dh* = −ρ

 ZHD=kp0 (4) where k is constant (2.28 mm/hPa) and p0 is the pressure at height h0 (Davis et al. 1985). It shows that the ZHD is proportional to the atmospheric pressure at the site. The ZWD is highly variable due possibly to varying climate, relating to the temperature and water vapour. The mean ZTD values at all GPS sites are shown in Figure 3 as a color map. It has noted that lower ZTD values are found at the areas of Tibet (Asia), Andes Mountain (South America), Northeast Pacific and higher latitudes (Antarctica and Arctic), and the higher ZTD values are concentrated at the areas of middle-low latitudes. In addition, the ZTD values decrease with increasing altitude, which is due to the atmospheric pressure variations with the height increase. Atmospheric pressure is the pressure above any area in the Earth's atmosphere caused by the weight of air. Air masses are affected by the general atmospheric pressure within the mass, creating areas of high pressure (anti-cyclones) and low pressure (depressions). Low pressure areas have less atmospheric mass above their locations, whereas high pressure areas have more atmospheric mass above their locations. As elevation increases, there are exponentially, fewer and fewer air molecules. Therefore, atmospheric pressure decreases with increasing altitude at a decreasing rate. The following relationship is a first-order approximation to the height (http://www.chemistrydaily.com/

( ) and *h* g(h) are the density and gravity at the height h. It can be further deduced as

()() (3)

**2.3 Global mean zenith tropospheric delay** 

chemistry/Atmospheric\_pressure):

where ρ

> The GPS ZTD time series have been analyzed for 4-12 years at globally distributed 150 GPS sites. Using the least square the fitting parameters of all GPS site are obtained, including trend and seasonal variation terms (Jin et al., 2007). The mean secular ZTD variation trend is about 1.5±0.001 mm/yr. Figure 4 shows the distribution of the secular ZTD variation trends at all GPS sites as the yearly increase or decrease. It can be seen that the trends are positive in most parts of the Northern Hemisphere and negative in most parts of the Southern Hemisphere (excluding positive in Antarctic), corresponding to a systematic increase or decrease of ZTD. It is interesting to note that the downtrend in Australia is larger than other regions. This downtrend of ZTD is probably due to the highly deserted in Australia. In addition, the ZTD variation trend decreases with increasing altitude, and furthermore, the ZTD trends are almost symmetrical with altitude. This indicates that the secular ZTD variations are larger at the lower altitude and at the higher altitude the secular ZTD variations hardly increase or decrease. In addition, the sum of downward and upward

GNSS Atmospheric and Ionospheric Sounding – Methods and Results 365

0 60 120 180 240 300 360


Latitude (deg)

The mean amplitude of semiannual ZTD variations is much smaller than annual variations with about 10 mm. The amplitudes of the semiannual oscillations are much smaller on the Southern Hemisphere than on the Northern Hemisphere. The distribution of the semiannual variation phase with the latitude has no clear symmetry. For example, at the latitudes of 40°N-50°N in the Northern Hemisphere, the semiannual phase is about 30° (about January), while at the latitudes of 40°S-50°S in the Southern Hemisphere, the semiannual phase is

Fig. 6. Distribution of annual ZTD variation amplitude with the latitude

5

25

45

65

85

105

125 mm

Longitude (deg)

Fig. 5. Annual variation amplitude of ZTD at globally distributed 150 GPS sites


0

20

40

60

Amplitude (mm)

about 330° (November).

80

100

120



0

Latitude (deg)

30

60

90

trends at globally distributed GPS sites is almost zero, which possibly indicates that the secular variation is in balance at a global scale, but subjecting to unevenly distributed GPS stations, etc. It need further be confirmed with much denser GPS network in the future. These secular ZTD variation characteristics reflect the total variations of surface atmospheric pressure, temperature and relative humidity, atmospheric vertical motions, etc.

Fig. 4. Secular trend of ZTD variations at global IGS sites. The red upwards arrows represent the increase of secular ZTD variations and the green downwards arrows stand for the decrease of secular ZTD variations

#### **2.5 Seasonal cycles**

Meanwhile, the seasonal components are also obtained using least square at annual and semiannual scales, which can be used to study the seasonal cycle, including amplitude and phase shift. The fitted phase shift is used to determine in which month the seasonal maximum takes place. The annual variation of ZTD ranges from 25 to 75 mm depending on the site, and the average amplitude is about 50 mm at most sites (Fig. 5). The annual variation amplitudes of ZTD at the IGS sites near Oceanic coasts are generally larger than in the continental inland. In addition, larger amplitudes of annual ZTD variation are mostly found at middle-low latitudes (near 20S° and 40N°), and the amplitudes of annual ZTD variation are especially smaller at higher latitudes (e.g. Antarctic and Arctic) and the equator areas (see Fig. 6). Sites on the eastern Atlantic and northeast Pacific coasts have lower annual variations, probably because of the moderating effect of the ocean on climate. Sites on the lee side of the Alps have higher annual variation, possibly due to the combined effects of a rain shadow in the winter and high moisture from the Mediterranean in the summer (Haase et al., 2003; Deblonde et al., 2005). Figure 7 shows the annual phase distribution with the latitude, where phase values are counted as clockwise from the north. It can be seen that the phase of annual ZTD variation is almost found at about 60° (about February, summer) in the Southern Hemisphere and at about 240° (about August, summer) in the Northern Hemisphere, which is just a half-year difference.

364 Global Navigation Satellite Systems – Signal, Theory and Applications

trends at globally distributed GPS sites is almost zero, which possibly indicates that the secular variation is in balance at a global scale, but subjecting to unevenly distributed GPS stations, etc. It need further be confirmed with much denser GPS network in the future. These secular ZTD variation characteristics reflect the total variations of surface atmospheric

0 60 120 180 240 300 360

2 mm/yr Scale

Longitude (deg)

Fig. 4. Secular trend of ZTD variations at global IGS sites. The red upwards arrows represent the increase of secular ZTD variations and the green downwards arrows stand for the decrease

Meanwhile, the seasonal components are also obtained using least square at annual and semiannual scales, which can be used to study the seasonal cycle, including amplitude and phase shift. The fitted phase shift is used to determine in which month the seasonal maximum takes place. The annual variation of ZTD ranges from 25 to 75 mm depending on the site, and the average amplitude is about 50 mm at most sites (Fig. 5). The annual variation amplitudes of ZTD at the IGS sites near Oceanic coasts are generally larger than in the continental inland. In addition, larger amplitudes of annual ZTD variation are mostly found at middle-low latitudes (near 20S° and 40N°), and the amplitudes of annual ZTD variation are especially smaller at higher latitudes (e.g. Antarctic and Arctic) and the equator areas (see Fig. 6). Sites on the eastern Atlantic and northeast Pacific coasts have lower annual variations, probably because of the moderating effect of the ocean on climate. Sites on the lee side of the Alps have higher annual variation, possibly due to the combined effects of a rain shadow in the winter and high moisture from the Mediterranean in the summer (Haase et al., 2003; Deblonde et al., 2005). Figure 7 shows the annual phase distribution with the latitude, where phase values are counted as clockwise from the north. It can be seen that the phase of annual ZTD variation is almost found at about 60° (about February, summer) in the Southern Hemisphere and at about 240° (about August, summer) in the Northern Hemisphere, which is just a half-year difference.


of secular ZTD variations

**2.5 Seasonal cycles** 



0

Latitude (deg)

30

60

90

pressure, temperature and relative humidity, atmospheric vertical motions, etc.

Fig. 5. Annual variation amplitude of ZTD at globally distributed 150 GPS sites

Fig. 6. Distribution of annual ZTD variation amplitude with the latitude

The mean amplitude of semiannual ZTD variations is much smaller than annual variations with about 10 mm. The amplitudes of the semiannual oscillations are much smaller on the Southern Hemisphere than on the Northern Hemisphere. The distribution of the semiannual variation phase with the latitude has no clear symmetry. For example, at the latitudes of 40°N-50°N in the Northern Hemisphere, the semiannual phase is about 30° (about January), while at the latitudes of 40°S-50°S in the Southern Hemisphere, the semiannual phase is about 330° (November).

GNSS Atmospheric and Ionospheric Sounding – Methods and Results 367

(a) mm

(b) mm

0

0

2

4

6

8

2

4

6

8


Longitude (o)


Longitude (o)

Fig. 8. Diurnal variation amplitudes (mm). (a) from GPS-derived ZTD and (b) from COADS

component ZHD accounts for approximately 90% of ZTD, ZTD is strongly correlated with surface pressure p0 at the site. It can be seen that if subdiurnal surface pressure varies by 1 hPa, the scale factor predicts a subdiurnal ZTD variation with amplitude of 2.28 mm. The

surface pressure data adjusted by a scale factor (2.28 mm/hPa)





0

Latitude (o)

30

60

90



0

Latitude (o)

30

60

90

Although the ZWD is small accounting for about 10 percent of the total zenith tropospheric delay (ZTD), the seasonal cycle in ZTD is due primarily to the wet component (ZWD). Furthermore, the seasonal variation phases of ZTD are almost consistent with surface temperature variations with the correlation coefficient of about 0.8. This reflects that annual and semiannual variations of ZTD are due mainly to the ZWD variations, about 80% in the surface temperature and 20% mainly in water vapour variations.

Fig. 7. Distribution of annual variation phase with the latitude. The phases are counted as clockwise from the north

#### **2.6 Diurnal and semidiurnal ZTD cycles**

The 4-7 years of GPS-derived ZTD time series has been used to analyze the diurnal and semidiurnal ZTD cycles and their features. Fig. 8a shows a colour coded map of diurnal ZTD amplitudes derived from the GPS data. The diurnal cycle (S1) has amplitudes between 0.2 and 10.9 mm with an uncertainty of about 0.5 mm. The diurnal ZTD amplitudes reduce with increasing latitude with the largest amplitudes appearing in the low-latitude equatorial areas, in particular in tropical Asia and the Gulf of Mexico. At these low latitudes, amplitudes of up to 10.9 mm are observed, while the high latitude areas reveal generally lower diurnal ZTD amplitudes. The peak values of the diurnal cycles occur spreading over the whole day (Fig. 9a). For the European stations there appears to be a preference for the second half of the day. At the semidiurnal cycle (S2), amplitudes between 0.1 and 4.3 mm with an uncertainty of about 0.2 mm are observed. Similar to the diurnal results mentioned above, the largest semidiurnal amplitudes are also found in low-latitude equatorial areas. The first peak of the semidiurnal cycle occurs typically around local noon. These diurnal and semidiurnal cycles of ZTD may be due to certain short time scale physical processes such as diurnal convection, atmospheric tides, general circulation and the coupling between the lower and the middle and upper atmosphere.

From Eq. (4), the zenith hydrostatic delay (ZHD) can be written as ZHD=2.28 p0. The scale factor k varies less than 1% even under severe weather conditions. As the hydrostatic 366 Global Navigation Satellite Systems – Signal, Theory and Applications

Although the ZWD is small accounting for about 10 percent of the total zenith tropospheric delay (ZTD), the seasonal cycle in ZTD is due primarily to the wet component (ZWD). Furthermore, the seasonal variation phases of ZTD are almost consistent with surface temperature variations with the correlation coefficient of about 0.8. This reflects that annual and semiannual variations of ZTD are due mainly to the ZWD variations, about 80% in the


Latitude (deg)

The 4-7 years of GPS-derived ZTD time series has been used to analyze the diurnal and semidiurnal ZTD cycles and their features. Fig. 8a shows a colour coded map of diurnal ZTD amplitudes derived from the GPS data. The diurnal cycle (S1) has amplitudes between 0.2 and 10.9 mm with an uncertainty of about 0.5 mm. The diurnal ZTD amplitudes reduce with increasing latitude with the largest amplitudes appearing in the low-latitude equatorial areas, in particular in tropical Asia and the Gulf of Mexico. At these low latitudes, amplitudes of up to 10.9 mm are observed, while the high latitude areas reveal generally lower diurnal ZTD amplitudes. The peak values of the diurnal cycles occur spreading over the whole day (Fig. 9a). For the European stations there appears to be a preference for the second half of the day. At the semidiurnal cycle (S2), amplitudes between 0.1 and 4.3 mm with an uncertainty of about 0.2 mm are observed. Similar to the diurnal results mentioned above, the largest semidiurnal amplitudes are also found in low-latitude equatorial areas. The first peak of the semidiurnal cycle occurs typically around local noon. These diurnal and semidiurnal cycles of ZTD may be due to certain short time scale physical processes such as diurnal convection, atmospheric tides, general circulation and the coupling between

From Eq. (4), the zenith hydrostatic delay (ZHD) can be written as ZHD=2.28 p0. The scale factor k varies less than 1% even under severe weather conditions. As the hydrostatic

Fig. 7. Distribution of annual variation phase with the latitude. The phases are counted as

surface temperature and 20% mainly in water vapour variations.

0

**2.6 Diurnal and semidiurnal ZTD cycles** 

the lower and the middle and upper atmosphere.

60

120

180

Phase (deg)

clockwise from the north

240

300

360

Fig. 8. Diurnal variation amplitudes (mm). (a) from GPS-derived ZTD and (b) from COADS surface pressure data adjusted by a scale factor (2.28 mm/hPa)

component ZHD accounts for approximately 90% of ZTD, ZTD is strongly correlated with surface pressure p0 at the site. It can be seen that if subdiurnal surface pressure varies by 1 hPa, the scale factor predicts a subdiurnal ZTD variation with amplitude of 2.28 mm. The

GNSS Atmospheric and Ionospheric Sounding – Methods and Results 369

Atmospheric Research (NCAR) (DS464.0; (http://dss.ucar.edu/datasets/ds464.0). Notably, these results are adjusted by a scale factor 2.28 mm/hPa, which are accounted for in the comparison. These pressure data origin from more than 8000 land and ocean weather stations including the Global Telecomuncation System (GTS) and marine reports from the Comprehensive Ocean-Atmosphere Data Set (COADS) (Dai and Wang 1999). The plots of diurnal and semidiurnal ZTD cycles show general similarities, indicating that the diurnal and semidiurnal atmospheric tides are probably the main driver of the diurnal and

The GPS consists of a constellation of 24 operating satellites in six circular orbits 20,200 km above the Earth at an inclination angle of 55º with a 12-h period. The satellite transmits two frequencies of signals (f1 = 1575.42 MHz and f2 = 1227.60 MHz). The equations of carrier phase (L) and code observations (pseudorange P) of double frequency GPS can be expressed

, , ,, , , , ( )( ) *i i i i i ii L d dc b <sup>k</sup> <sup>j</sup>*

, , , , , , , ( ) *i ii i i i P d d c dd <sup>k</sup> <sup>j</sup>*

where superscript i and subscript j represent the satellite and ground-based GPS receiver,

is the satellite or receiver clock offset, *b* is the phase delay of satellite and receiver instrument

ambiguity of carrier phase, and *ε* is the other residuals. From Eq.(7) and (8), the ionospheric delay can be determined, which is useful for ionospheric delay correction and space weather.

The ionospheric delay can be determined from the double frequency GPS phase and code

41 2 2 2 4 1 2 1 1 40.3 () ( ,) *i i L F j j <sup>z</sup> VTEC <sup>s</sup> <sup>B</sup> f f*

41 2 2 2 4 1 2 1 1 40.3 () ( ,) *i i PP P j j F z VTEC s b f f*

*j j dq* − +− *dq dq dq* . The Differential Code Biases (*b*4) can be obtained through GPS

where N is the epoch of GPS observation (Jin et al., 2008). For the TEC representation, a

carrier phase observations, and *B*4 can be obtained through the formula, 4 44

λ

where *F z*( ) is the mapping function, *B*4 is ( 4 11 1 22 2 ( )( )) *ii ii B bN bN* =− + + +

⎛ ⎞ <sup>=</sup> −= − ⎜ ⎟ <sup>+</sup> ⎝ ⎠

⎛ ⎞ <sup>=</sup> − =− − ⎜ ⎟ <sup>+</sup> ⎝ ⎠

ionospheric and tropospheric delays, respectively, *c* is the speed of light in vacuum space,

bias, *<sup>q</sup> d* is the code delay of satellite and receiver instrumental bias,

*k k <sup>j</sup> ion k <sup>j</sup> trop j <sup>j</sup> k k <sup>j</sup> Nk <sup>j</sup>* = =− + + − − +

 *ion k <sup>j</sup> trop j j q <sup>k</sup> <sup>q</sup> <sup>k</sup> j j* = + + + −+ + + τ τ

ττ

is the distance between satellite i and GPS receiver j, and *ion trop d d* are the

is the total carrier phase between the satellite and receiver, *N* is the

 λ

β

β

 λ ε

(7)

is the carrier

(9)

(10)

( )/

∑ + − ,

*p Lb N*

*jj jj* , and 4 *b* is

1

*i*

=

*N*

τ

(8)

λ

semidiurnal ZTD variations derived from GPS (Jin et al., 2009).

 ρ

λφ

ρ

φ φ

**3. Ionospheric sounding** 

ρ

φ

**3.1 2-D ionospheric imaging** 

(pseudorange) observations as

1 2 12 ( )( ) *i i*

as follows:

respectively,

wavelength,

Fig. 9. Time of diurnal peak values at local time (LT: hour) where at each GPS sites longitude the Sun is at its highest elevation at 12:00 LT. (a) from global IGS GPS observations and (b) from COADS surface pressure data

GPS-derived S1 and S2 signals in ZTD are compared with three-hour surface synoptic pressure observations from 1997 to 2007 that are archived at the National Center for Atmospheric Research (NCAR) (DS464.0; (http://dss.ucar.edu/datasets/ds464.0). Notably, these results are adjusted by a scale factor 2.28 mm/hPa, which are accounted for in the comparison. These pressure data origin from more than 8000 land and ocean weather stations including the Global Telecomuncation System (GTS) and marine reports from the Comprehensive Ocean-Atmosphere Data Set (COADS) (Dai and Wang 1999). The plots of diurnal and semidiurnal ZTD cycles show general similarities, indicating that the diurnal and semidiurnal atmospheric tides are probably the main driver of the diurnal and semidiurnal ZTD variations derived from GPS (Jin et al., 2009).

## **3. Ionospheric sounding**

368 Global Navigation Satellite Systems – Signal, Theory and Applications

(a) Hour

(b) Hour

0

0

6

12

18

24

6

12

18

24


Longitude (o)


Longitude (o)

Fig. 9. Time of diurnal peak values at local time (LT: hour) where at each GPS sites longitude the Sun is at its highest elevation at 12:00 LT. (a) from global IGS GPS observations and

GPS-derived S1 and S2 signals in ZTD are compared with three-hour surface synoptic pressure observations from 1997 to 2007 that are archived at the National Center for



(b) from COADS surface pressure data



0

Latitude 9o

)

30

60

90



0

Latitude (o)

30

60

90

The GPS consists of a constellation of 24 operating satellites in six circular orbits 20,200 km above the Earth at an inclination angle of 55º with a 12-h period. The satellite transmits two frequencies of signals (f1 = 1575.42 MHz and f2 = 1227.60 MHz). The equations of carrier phase (L) and code observations (pseudorange P) of double frequency GPS can be expressed as follows:

$$L\_{k,j}^i = \lambda\_k \phi\_{k,j}^i = \rho - d\_{ion,k,j}^i + d\_{trop,j}^i + \mathfrak{c}(\boldsymbol{\pi}^i - \boldsymbol{\tau}\_j) - \lambda\_k (b\_{k,j}^i + N\_{k,j}^i) \tag{7}$$

$$P\_{k,j}^{i} = \rho + d\_{ion,k,j}^{i} + d\_{trop,j}^{i} + c(\boldsymbol{\pi}^{i} - \boldsymbol{\pi}\_{j}) + d\_{q,k}^{i} + d\_{q,k,j} + \boldsymbol{\pi}\_{j}^{i} \tag{8}$$

where superscript i and subscript j represent the satellite and ground-based GPS receiver, respectively, ρ is the distance between satellite i and GPS receiver j, and *ion trop d d* are the ionospheric and tropospheric delays, respectively, *c* is the speed of light in vacuum space, τ is the satellite or receiver clock offset, *b* is the phase delay of satellite and receiver instrument bias, *<sup>q</sup> d* is the code delay of satellite and receiver instrumental bias, λ is the carrier wavelength, φ is the total carrier phase between the satellite and receiver, *N* is the ambiguity of carrier phase, and *ε* is the other residuals. From Eq.(7) and (8), the ionospheric delay can be determined, which is useful for ionospheric delay correction and space weather.

#### **3.1 2-D ionospheric imaging**

The ionospheric delay can be determined from the double frequency GPS phase and code (pseudorange) observations as

$$L\_4 = \phi\_{1j}^i - \phi\_{2j}^i = -40.3 \left( \frac{1}{f\_1^2} - \frac{1}{f\_2^2} \right) F(z) V TEC(\beta, s) + B\_4 \tag{9}$$

$$P\_4 = P\_{1j}^i - P\_{2j}^i = 40.3 \left( \frac{1}{f\_1^2} - \frac{1}{f\_2^2} \right) \text{F}(\text{z}) V \text{TEC}(\beta, \text{s}) + b\_4 \tag{10}$$

where *F z*( ) is the mapping function, *B*4 is ( 4 11 1 22 2 ( )( )) *ii ii B bN bN* =− + + + λ λ *jj jj* , and 4 *b* is 1 2 12 ( )( ) *i i j j dq* − +− *dq dq dq* . The Differential Code Biases (*b*4) can be obtained through GPS carrier phase observations, and *B*4 can be obtained through the formula, 4 44 1 ( )/ *N i p Lb N* = ∑ + − , where N is the epoch of GPS observation (Jin et al., 2008). For the TEC representation, a

GNSS Atmospheric and Ionospheric Sounding – Methods and Results 371

the MART performs a multiplicative modification in each iteration, and thus the inversion results are always positive. Therefore, MART has the advantage over ART in determining the electron densities that avoid unreasonable negative values and is the one used in this

*k ij a*

λ

*<sup>k</sup>* until the residual does not change (see Figure 10). This

*<sup>k</sup>* =0. 01 has been chosen for all iterations. In addition, it is

*x*

ε<− <sup>+</sup><sup>1</sup> . *<sup>k</sup> xAy*

Reconstruction solution

Yes

(15)

*<sup>k</sup>* < , and the inner product of the

λ

λ*<sup>k</sup>* is

*<sup>k</sup>* value is identified

λ

<sup>1</sup> . ,

*a x*

⎛ ⎞ ⎜ ⎟ <sup>=</sup> ⎜ ⎟ ⎝ ⎠

where yi is the ith observed STEC in a column of m measurements, xj is the jth resulted cell electron density in a column of n unknowns, aij is the length of link i that lies in cell j,

vectors x and ai is the simulated STEC for the ith path. The electron density matrix x is therefore corrected iteratively by the ratio of the measured STEC and the simulated STEC

where the solution converges quickly with a reasonable number of iterations and the

noted that any iterative algorithm requires an initial condition before the iteration begins. Due to the poor STEC geometry, the initialization could be extremely important for the tomographic reconstruction. In practice, the closer the initial condition is to the true electron density distribution, the more accurate the reconstruction will be. Here the latest IRI-2007 model (http://nssdcftp.gsfc.nasa.gov/models/ionospheric/iri/iri2007) is used as an initial

*<sup>k</sup> x <sup>k</sup>*+<sup>1</sup>

*k i*

*i xa y*

. 〉〈

⎜ ⎜ ⎝

⎛

*ij ak*

λ

⎟ ⎟ ⎠

⎞

*k k i j j k i*

*<sup>y</sup> x x*

+

study. Basically, the MART algorithm is iterated cyclically:

the relaxation parameter at the kth iteration with 0 1 <

λ

λ

relaxation parameter value is chosen from experience in which the best

No (iteration)

with a relaxation parameter of

residuals are a minimum. Here

guess for the reconstruction iteration.

0 *x*

Fig. 10. Flow chart of MART

single layer model (SLM) ionosphere approximation was used. SLM assumes that all the free electrons are contained in a shell of infinitesimal thickness at altitude H (generally 350 km above the Earth). A mapping function is used to convert the slant TEC into the vertical TEC (VTEC) as shown:

$$F(z) = \sqrt{(1 - R\cos(90 - z)/(R + H))}\tag{11}$$

where R is Earth radius, H is SLM height, and z is satellite zenith angle. When using the above mapping function, F(z), one can obtain VTEC values at the ionosphere pierce points (IPPs). The GPS-derived TEC can correct ionospheric delay for microwave techniques and monitor space weather events.

#### **3.2 3-D ionospheric tomography reconstruction**

The STEC is defined as the line integral of the electron density as expressed by:

$$STEC = \int\_{R\_{\text{matrix}}}^{R\_{\text{stable}}} N\_e(\mathcal{X}, \rho, h) ds \tag{12}$$

where (,,) *N h <sup>e</sup>* λ ϕ is the ionospheric electron density, λ , ϕ and *h* are the longitude, latitude and height, respectively. To obtain *Ne* , the ionosphere is divided into grid pixels with a small cell where the electron density is assumed to be constant, so that the STEC in Eq.(4) along the ray path *i* can be approximately written as a finite sum over the pixels *j* as follows:

$$STEC\_i = \sum\_{j=1}^{M} a\_{ij} n\_j \tag{13}$$

where *ij a* is a matrix whose elements denote the length of the path-pixel intersections in the pixel *j* along the ray path *i* , and *nj* is the electron density for the pixel *j* . Each set of STEC measurements along the ray paths from all observable satellites at consecutive epochs are combined with the ray path geometry into a linear expression:

$$Y = A\mathfrak{x} + \mathfrak{z} \tag{14}$$

where Y is a column of m measurements of STEC, x is a column of n electron density unknowns for cells in the targeted ionosphere region, and A is an m × n normal matrix with elements *ij a* . The unknown electron densities *x* can be estimated by the ionospheric tomographic reconstruction technique. Many tomography algorithms are used in different ways, e.g. algebraic reconstruction technique (Gordon et al., 1970). One of the most common approaches is the algebraic reconstruction technique (ART), which was first introduced in Computerized Ionospheric Tomography (CIT) by Austen et al. (1986). This is an iterative procedure for solving a linear equation. A modified version of ART is the so-called multiplicative ART (MART), where the correction in each iteration is obtained by making a multiplicative modification to *x* (Raymund et al., 1990; Tsai et al., 2002). The ART generally produces estimates of the unknown parameters by minimization of the L2 norm, while the MART follows maximum entropy criteria and thus underlies different statistics. In addition, 370 Global Navigation Satellite Systems – Signal, Theory and Applications

single layer model (SLM) ionosphere approximation was used. SLM assumes that all the free electrons are contained in a shell of infinitesimal thickness at altitude H (generally 350 km above the Earth). A mapping function is used to convert the slant TEC into the vertical

where R is Earth radius, H is SLM height, and z is satellite zenith angle. When using the above mapping function, F(z), one can obtain VTEC values at the ionosphere pierce points (IPPs). The GPS-derived TEC can correct ionospheric delay for microwave techniques and

(,,) *satellite*

λ ϕ

> λ , ϕ

*receiver R <sup>e</sup> <sup>R</sup> STEC N h ds* =

latitude and height, respectively. To obtain *Ne* , the ionosphere is divided into grid pixels with a small cell where the electron density is assumed to be constant, so that the STEC in Eq.(4) along the ray path *i* can be approximately written as a finite sum over the pixels *j*

1

ε

where *ij a* is a matrix whose elements denote the length of the path-pixel intersections in the pixel *j* along the ray path *i* , and *nj* is the electron density for the pixel *j* . Each set of STEC measurements along the ray paths from all observable satellites at consecutive epochs are

*Y Ax* = +

where Y is a column of m measurements of STEC, x is a column of n electron density unknowns for cells in the targeted ionosphere region, and A is an m × n normal matrix with elements *ij a* . The unknown electron densities *x* can be estimated by the ionospheric tomographic reconstruction technique. Many tomography algorithms are used in different ways, e.g. algebraic reconstruction technique (Gordon et al., 1970). One of the most common approaches is the algebraic reconstruction technique (ART), which was first introduced in Computerized Ionospheric Tomography (CIT) by Austen et al. (1986). This is an iterative procedure for solving a linear equation. A modified version of ART is the so-called multiplicative ART (MART), where the correction in each iteration is obtained by making a multiplicative modification to *x* (Raymund et al., 1990; Tsai et al., 2002). The ART generally produces estimates of the unknown parameters by minimization of the L2 norm, while the MART follows maximum entropy criteria and thus underlies different statistics. In addition,

*M i ij j j STEC a n* =

The STEC is defined as the line integral of the electron density as expressed by:

is the ionospheric electron density,

combined with the ray path geometry into a linear expression:

*Fz R z R H* ( ) (1 cos(90 ) ( )) =− − + (11)

∫ (12)

<sup>=</sup> ∑ (13)

(14)

and *h* are the longitude,

TEC (VTEC) as shown:

monitor space weather events.

where (,,) *N h <sup>e</sup>* λ ϕ

as follows:

**3.2 3-D ionospheric tomography reconstruction** 

the MART performs a multiplicative modification in each iteration, and thus the inversion results are always positive. Therefore, MART has the advantage over ART in determining the electron densities that avoid unreasonable negative values and is the one used in this study. Basically, the MART algorithm is iterated cyclically:

$$\mathbf{x}\_{j}^{k+1} = \mathbf{x}\_{j}^{k} \cdot \left(\frac{\mathbf{y}\_{i}}{\left\{a\_{i}, \mathbf{x}^{k}\right\}}\right)^{\mathbf{x}\_{k} a\_{ij}} \tag{15}$$

where yi is the ith observed STEC in a column of m measurements, xj is the jth resulted cell electron density in a column of n unknowns, aij is the length of link i that lies in cell j, λ*<sup>k</sup>* is the relaxation parameter at the kth iteration with 0 1 < λ*<sup>k</sup>* < , and the inner product of the vectors x and ai is the simulated STEC for the ith path. The electron density matrix x is therefore corrected iteratively by the ratio of the measured STEC and the simulated STEC with a relaxation parameter of λ*<sup>k</sup>* until the residual does not change (see Figure 10). This relaxation parameter value is chosen from experience in which the best λ*<sup>k</sup>* value is identified where the solution converges quickly with a reasonable number of iterations and the residuals are a minimum. Here λ*<sup>k</sup>* =0. 01 has been chosen for all iterations. In addition, it is noted that any iterative algorithm requires an initial condition before the iteration begins. Due to the poor STEC geometry, the initialization could be extremely important for the tomographic reconstruction. In practice, the closer the initial condition is to the true electron density distribution, the more accurate the reconstruction will be. Here the latest IRI-2007 model (http://nssdcftp.gsfc.nasa.gov/models/ionospheric/iri/iri2007) is used as an initial guess for the reconstruction iteration.

Fig. 10. Flow chart of MART

GNSS Atmospheric and Ionospheric Sounding – Methods and Results 373

GPS IRI-2007 Ionosonde

0 2 4 6 8 10

Electron density (1011 el/m3)

As the F2-layer peak electron density value (denoted as NmF2, proportional to the square of the F2 layer critical frequency foF2) is of great influence on the shape of the ionospheric electron density profile Ne(h), and also probably related to various physical processes of ionospheric activities as well as the F2-layer peak height (hmF2), the NmF2 and hmF2 are essential parameters for monitoring ionospheric activities and understanding the nature of the ionosphere. Although many studies of geomagnetic storms have been carried out in the US and Europe (e.g. Foster et al. 2004; Yin et al. 2004; Goncharenko et al. 2007), investigations on ionospheric behaviour to storms in Asia are relatively few due to a lack of dense sets of GPS observations, etc. In this paper, the responses of the GPS-derived NmF2 and hmF2 to the super geomagnetic storm (20 November 2003) over South Korea using data from the dense Korean GPS Network are described. First, ground-based dual-frequency GPS observations are used to produce electron density profiles using the ionospheric tomography technique, which are verified by independent ionosonde data. Then the responses of the key ionospheric F2-layer parameters NmF2 and hmF2 to the 20 November 2003 geomagnetic storm are investigated over South Korea to gain insights into the effects of

In the following, the 3-D ionospheric disturbances during the large November 20th 2003 geomagnetic storm are investigated and analyzed using GPS data in South Korea. The geomagnetic storm Dst, Kp and AE indices on 20-21 November 2003 obtained from the World Data Center in Kyoto (http://swdcdb.kugi.kyoto-u.ac.jp/) showed a strong

Fig. 11. Comparison of the electron density profiles derived from the ground-based GPS tomography reconstruction (solid line), ionosonde observation at Anyang stations (37.39ºN, 126.95ºE) (dot line) and IRI-2007 estimation (dashed line) at 9:00 UT on 1 October 2003

100

different physical conditions and processes.

geomagnetic storm on November 20th, 2003 (Figure 2).

200

300

400

Height (km)

500

600

700

800

The ionospheric reconstruction algorithm MART can integrate the STEC from all available GPS receivers to all GPS satellites visible from each site of the KGN network above a userspecified elevation cut-off angle (usually 15°). The unknown electron density profile is expressed in 4-D (longitude-latitude-height and time) voxel basis functions over the following grid: longitude 124E°-130E° in 1° increments, latitude 33N°-39N° in 0.5° increments, altitude 100 km - 1000 km in 25 km increments and time: 1 h increments of linear change in the electron density per voxel. As there are a sparse number of ions in the ionosphere above 1000 km, the effect on the inversion for ionospheric electron density profiles is very small and the ionosphere is only considered up to an altitude of 1000 km. Furthermore, it is faster to invert the unknown ionospheric density parameters because of the reduced number of unknown variables. In addition, the fewer leaving rays from the ionospheric space of above defined latitude and longitude range are not used, but are useful to further obtain ionospheric profiles outside the latitudinal/longitudinal boundaries space. Using the STEC of all ray paths passing the ionospheric grid cells from the Korean GPS network, the 4-D ionospheric electron density profiles can be derived through the tomography reconstruction algorithm. To verify the reliability of GPS ionospheric tomography reconstruction results, the available ionosonde station (Anyang) in South Korea provides an independent comparison with the GPS tomographically reconstructed electron density profiles. The electron density profiles at 25 km height steps from GPS reconstruction and ionosonde data match well in October and November with a root-mean-square (RMS) of 11 3 0.3 10 / × *el m* . For example, Figure 11 shows a comparison of the GPS reconstructed electron density profile at 9:00 UT on 1 October 2003 with the available ionosonde data at Anyang station (37.39°N, 126.95°E) and the profiles from the IRI-2007 model. It can be seen that the GPS tomographically reconstructed density profile is in a good agreement with ionosonde data and the IRI-2007 model, but is closer to the ionosonde, which confirms the validity of our GPS ionospheric reconstruction approach (Jin et al., 2006 and 2008).

#### **3.3 F-2 layer ionospheric response to storm**

It is well known that geomagnetic storms may profoundly affect the global ionosphere and upper atmosphere, inducing great variations in such parameters as the Total Electron Content (TEC), the F2-layer peak density (NmF2) and its height (hmF2). These influences vary with location, season, local time and solar activity. The responses of the ionosphere to geomagnetic storms have been studied for several decades using moderately priced to expensive instrumentation, such as ionosondes and Incoherent Scatter Radar (ISR) (Lei et al. 2004). However, it is well known that ionospheric storms have a global impact on ionization, and under very disturbed conditions the ionospheric response to severe storms often presents significant changes in the distribution of ionization with latitude and altitude. Furthermore, ionosondes cannot measure the topside ionosphere and sometimes suffer from absorption during storms, whereas Incoherent Scatter Radars have geographical limitations. Nowadays the GPS satellites, being in high altitude orbits (~20,200 km), are very useful for studying the structure of the entire ionosphere, even the plasmasphere. Moreover, GPS is a low-cost, all-weather, near real time, and high-resolution atmospheric sounding technique. Therefore, GPS has been widely used to monitor the ionosphere (e.g. Jin et al. 2004-2007; Yin et al., 2004).

372 Global Navigation Satellite Systems – Signal, Theory and Applications

The ionospheric reconstruction algorithm MART can integrate the STEC from all available GPS receivers to all GPS satellites visible from each site of the KGN network above a userspecified elevation cut-off angle (usually 15°). The unknown electron density profile is expressed in 4-D (longitude-latitude-height and time) voxel basis functions over the following grid: longitude 124E°-130E° in 1° increments, latitude 33N°-39N° in 0.5° increments, altitude 100 km - 1000 km in 25 km increments and time: 1 h increments of linear change in the electron density per voxel. As there are a sparse number of ions in the ionosphere above 1000 km, the effect on the inversion for ionospheric electron density profiles is very small and the ionosphere is only considered up to an altitude of 1000 km. Furthermore, it is faster to invert the unknown ionospheric density parameters because of the reduced number of unknown variables. In addition, the fewer leaving rays from the ionospheric space of above defined latitude and longitude range are not used, but are useful to further obtain ionospheric profiles outside the latitudinal/longitudinal boundaries space. Using the STEC of all ray paths passing the ionospheric grid cells from the Korean GPS network, the 4-D ionospheric electron density profiles can be derived through the tomography reconstruction algorithm. To verify the reliability of GPS ionospheric tomography reconstruction results, the available ionosonde station (Anyang) in South Korea provides an independent comparison with the GPS tomographically reconstructed electron density profiles. The electron density profiles at 25 km height steps from GPS reconstruction and ionosonde data match well in October and November with a root-mean-square (RMS) of 11 3 0.3 10 / × *el m* . For example, Figure 11 shows a comparison of the GPS reconstructed electron density profile at 9:00 UT on 1 October 2003 with the available ionosonde data at Anyang station (37.39°N, 126.95°E) and the profiles from the IRI-2007 model. It can be seen that the GPS tomographically reconstructed density profile is in a good agreement with ionosonde data and the IRI-2007 model, but is closer to the ionosonde, which confirms the validity of our GPS ionospheric reconstruction approach

It is well known that geomagnetic storms may profoundly affect the global ionosphere and upper atmosphere, inducing great variations in such parameters as the Total Electron Content (TEC), the F2-layer peak density (NmF2) and its height (hmF2). These influences vary with location, season, local time and solar activity. The responses of the ionosphere to geomagnetic storms have been studied for several decades using moderately priced to expensive instrumentation, such as ionosondes and Incoherent Scatter Radar (ISR) (Lei et al. 2004). However, it is well known that ionospheric storms have a global impact on ionization, and under very disturbed conditions the ionospheric response to severe storms often presents significant changes in the distribution of ionization with latitude and altitude. Furthermore, ionosondes cannot measure the topside ionosphere and sometimes suffer from absorption during storms, whereas Incoherent Scatter Radars have geographical limitations. Nowadays the GPS satellites, being in high altitude orbits (~20,200 km), are very useful for studying the structure of the entire ionosphere, even the plasmasphere. Moreover, GPS is a low-cost, all-weather, near real time, and high-resolution atmospheric sounding technique. Therefore, GPS has been widely used to monitor the ionosphere (e.g. Jin et al. 2004-2007; Yin

(Jin et al., 2006 and 2008).

et al., 2004).

**3.3 F-2 layer ionospheric response to storm** 

Fig. 11. Comparison of the electron density profiles derived from the ground-based GPS tomography reconstruction (solid line), ionosonde observation at Anyang stations (37.39ºN, 126.95ºE) (dot line) and IRI-2007 estimation (dashed line) at 9:00 UT on 1 October 2003

As the F2-layer peak electron density value (denoted as NmF2, proportional to the square of the F2 layer critical frequency foF2) is of great influence on the shape of the ionospheric electron density profile Ne(h), and also probably related to various physical processes of ionospheric activities as well as the F2-layer peak height (hmF2), the NmF2 and hmF2 are essential parameters for monitoring ionospheric activities and understanding the nature of the ionosphere. Although many studies of geomagnetic storms have been carried out in the US and Europe (e.g. Foster et al. 2004; Yin et al. 2004; Goncharenko et al. 2007), investigations on ionospheric behaviour to storms in Asia are relatively few due to a lack of dense sets of GPS observations, etc. In this paper, the responses of the GPS-derived NmF2 and hmF2 to the super geomagnetic storm (20 November 2003) over South Korea using data from the dense Korean GPS Network are described. First, ground-based dual-frequency GPS observations are used to produce electron density profiles using the ionospheric tomography technique, which are verified by independent ionosonde data. Then the responses of the key ionospheric F2-layer parameters NmF2 and hmF2 to the 20 November 2003 geomagnetic storm are investigated over South Korea to gain insights into the effects of different physical conditions and processes.

In the following, the 3-D ionospheric disturbances during the large November 20th 2003 geomagnetic storm are investigated and analyzed using GPS data in South Korea. The geomagnetic storm Dst, Kp and AE indices on 20-21 November 2003 obtained from the World Data Center in Kyoto (http://swdcdb.kugi.kyoto-u.ac.jp/) showed a strong geomagnetic storm on November 20th, 2003 (Figure 2).

GNSS Atmospheric and Ionospheric Sounding – Methods and Results 375

Monthly median

GPS by this study

Ionosonde

0 5 10 15 20 25

Time (UT: hour)

Fig. 13. The hmF2 variations from monthly median, ionosonde and GPS observations on 20

Normally the increase/loss rate of F region electron density depends mainly on the molecular nitrogen concentration [N2] and atomic oxygen concentration [O] [8]. However, the O/N2 ratio obtained by the GUVI instrument on board the TIMED satellite doesn't show significant changes in South Korea where the increased NmF2 was observed, indicating that the increased NmF2 in South Korea is not caused by changes in neutral composition, and other possible non-chemical effects, such as dynamical changes of vertical ion motions induced by winds and E×B drifts, tides and waves in the mesosphere/lower thermosphere (MLT) region, which can be dynamically coupled upward to generate

The GPS signals are delayed due to the effects of dry gas and water vapor when propagating through the neutral atmosphere. The hydrostatic delay is proportional to the surface pressure and the wet delay is a key parameter in atmospheric radiation, hydrological cycle, energy transfer and the formation of clouds via latent heat. Thereby, the total tropospheric delay (ZTD) is an important parameter of the atmosphere, which directly or indirectly reflects the weather and climate processes, variations, and atmospheric vertical motions, etc. Comparing the time series of the zenith tropospheric delay (ZTD), zenith hydrostatic delay (ZHD), zenith wet delay (ZWD), surface temperature, pressure and relative humidity, it has been noted that the ZHD is highly proportional to the atmospheric pressure at the site and relatively stable and the ZWD is positively correlated with the

250

ionospheric perturbations and oscillations.

300

350

400

hmF2 (Km)

November 2003

**4. Conclusion** 

450

500

550

Fig. 12. The geomagnetic storm index (Dst) (upper), Kp (middle) and AE indices (bottom) on 20-21 November 2003

As the F2-layer peak electron density (NmF2) and its height (hmF2) are main parameters of the ionospheric electron density profile Ne (h), the behaviour of the ionospheric F2-layer to the storm was investigated over South Korea in terms of the NmF2 and hmF2. Here the peak density (NmF2) and its corresponding height (hmF2) are obtained from the groundbased GPS observations using the MART reconstructed tomography technique. The monthly median value of GPS reconstructed electron density profiles during the quiet days is regarded as the reference and the deviation of ionospheric NmF2 and hmF2 can reflect the ionospheric behaviors during the geomagnetic storm. It has shown that the GPS-derived NmF2 has a disturbance at 9:00 UT and then increases from 10:00 UT until 19:00 UT. The corresponding hmF2 also suddenly rises from 8:00 UT when the storm just started, and reaches the maximum height at about 16:00 UT with a maximum Kp value of 9, and then gradually descends until 21:00 UT (Fig.13), which are also supported by anther independent ionosonde measurement at Anyang station.

Fig. 13. The hmF2 variations from monthly median, ionosonde and GPS observations on 20 November 2003

Normally the increase/loss rate of F region electron density depends mainly on the molecular nitrogen concentration [N2] and atomic oxygen concentration [O] [8]. However, the O/N2 ratio obtained by the GUVI instrument on board the TIMED satellite doesn't show significant changes in South Korea where the increased NmF2 was observed, indicating that the increased NmF2 in South Korea is not caused by changes in neutral composition, and other possible non-chemical effects, such as dynamical changes of vertical ion motions induced by winds and E×B drifts, tides and waves in the mesosphere/lower thermosphere (MLT) region, which can be dynamically coupled upward to generate ionospheric perturbations and oscillations.

#### **4. Conclusion**

374 Global Navigation Satellite Systems – Signal, Theory and Applications

0:00 8:00 16:00 24:00 8:00 16:00 24:00

0:00 8:00 16:00 24:00 8:00 16:00 24:00

0:00 8:00 16:00 24:00 8:00 16:00 24:00

20 Nov. 2003 21 Nov. 2003

Fig. 12. The geomagnetic storm index (Dst) (upper), Kp (middle) and AE indices (bottom) on

As the F2-layer peak electron density (NmF2) and its height (hmF2) are main parameters of the ionospheric electron density profile Ne (h), the behaviour of the ionospheric F2-layer to the storm was investigated over South Korea in terms of the NmF2 and hmF2. Here the peak density (NmF2) and its corresponding height (hmF2) are obtained from the groundbased GPS observations using the MART reconstructed tomography technique. The monthly median value of GPS reconstructed electron density profiles during the quiet days is regarded as the reference and the deviation of ionospheric NmF2 and hmF2 can reflect the ionospheric behaviors during the geomagnetic storm. It has shown that the GPS-derived NmF2 has a disturbance at 9:00 UT and then increases from 10:00 UT until 19:00 UT. The corresponding hmF2 also suddenly rises from 8:00 UT when the storm just started, and reaches the maximum height at about 16:00 UT with a maximum Kp value of 9, and then gradually descends until 21:00 UT (Fig.13), which are also supported by anther independent

�400

0

0

ionosonde measurement at Anyang station.

1000

20-21 November 2003

2000

AE (nT)

3000

3

6

9

�200

Dst (nT)

Kp

0

The GPS signals are delayed due to the effects of dry gas and water vapor when propagating through the neutral atmosphere. The hydrostatic delay is proportional to the surface pressure and the wet delay is a key parameter in atmospheric radiation, hydrological cycle, energy transfer and the formation of clouds via latent heat. Thereby, the total tropospheric delay (ZTD) is an important parameter of the atmosphere, which directly or indirectly reflects the weather and climate processes, variations, and atmospheric vertical motions, etc. Comparing the time series of the zenith tropospheric delay (ZTD), zenith hydrostatic delay (ZHD), zenith wet delay (ZWD), surface temperature, pressure and relative humidity, it has been noted that the ZHD is highly proportional to the atmospheric pressure at the site and relatively stable and the ZWD is positively correlated with the

GNSS Atmospheric and Ionospheric Sounding – Methods and Results 377

possible non-chemical effects, such as dynamical changes of vertical ion motions induced by winds and E×B drifts, tides and waves in the mesosphere/lower thermosphere (MLT) region, which can be dynamically coupled upward to generate ionospheric perturbations

Ground-based and space borne GPS observations have been widely used in atmospheric sounding, including sensing tropopsheric precipitable water vapor (PWV), ionospheric total electron content (TEC) and atmospheric profile information (e.g., pressure, temperature, humidity, tropopause and ionospheric electron density). These observations have facilitated greater advancements in meteorology, climatology, numerical weather model, atmospheric science and space weather (e.g., Jin et al., 2007; Jin et al., 2009; Schmidt, et al. 2010). For example, the dual frequency ground GPS array could detect ionospheric response and its processes during large geomagnetic storms (Jin et al., 2008). Meanwhile, ground GPS also observed the plasma bubbles and retrieved reliable propagation characteristics of the depletions without assumptions about the mapping of the depletion along magnetic field

lines to large latitudinal distances, comparable with airglow data (Haase et al., 2011).

profile information and evolution processes of the atmosphere and ionosphere.

This work was supported by the Shanghai Pujiang Talent Program, National Natural Science Foundation of China (Grant No.11043008) and Key Direction Project of Chinese

Austen, J. R.; S. J. Franke, and C. H. Liu, Application of computerized tomography

techniques to ionospheric research, in Radio Beacon contributions to the study of ionization and dynamics of the ionosphere and to corrections to geodesy and

Additionally, space-borne GPS receivers have proven very successful in making high vertical resolution and global atmospheric measurements using the radio occultation (RO) technique The existing GPS radio occultation (RO) missions have been widely used to estimate the detailed vertical profile information, including pressure, temperature, gravity waves and sporadic E-layers as well as their variation characteristics, particularly six satellites of the Taiwan/US FORMOSAT-3/COSMIC (FORMOsa SATellite mission-3/Constellation Observing System for Meteorology, Ionosphere and Climate) mission with more than 2000 radio occultation profiles per day. Schmidt et al. (2010) was the first to observe upper tropospheric warming and lower stratospheric cooling using GPS RO data (2001-2009). Although a number of progresses in atmospheric and ionospheric sensing have been made using GPS RO missions in the past few years, they still do not satisfy actual requirements for short-time scales and higher temporal-spatial resolution monitoring together with ground GNS observations. For instance, the tropospheric or ionospheric profile information cannot be directly estimated from GPS tomography due to the lack of enough line-of-sight GPS signals passing each grid cell (Jin et al., 2006 and 2008; Nesterov and Kunitsyn, 2011). Moreover, most current RO satellite missions are approaching their end of operations. With the increase of future GNSS satellite constellations and more GNSS RO missions, the goal of improved temporal-spatial resolution will enable more detailed

and oscillations.

**5. Acknowledgement** 

**6. References** 

Academy of Sciences (Grant No. KJCX2- EW-T03).

temperature and correlated with the relative humidity. The mean correlation coefficient between ZWD and surface temperature is about 0.80 and the correlation coefficient between ZTD and ZWD is about 0.95, reflecting a good agreement between ZTD and ZWD variations. Therefore, the seasonal cycles of the ZTD are due primarily to the wet component (ZWD), especially in the surface temperature, even though the wet delay is only 10% of the total delay (ZTD).

The lower mean ZTD values are located at the areas of higher altitude (e.g. Tibet, Asia and Andes Mountain, South America) and higher latitude areas (Antarctica and Arctic), and the higher mean ZTD values are concentrated at the areas of middle-low latitudes. The mean ZTD decreases with increasing altitude at an exponentially decreasing rate due to the atmospheric pressure decreasing with the height increasing. The mean secular ZTD variation trend is about 1.5±0.001 mm/yr at all GPS sites. The secular variations are systematically increasing in most parts of the Northern Hemisphere and decreasing in most parts of the Southern Hemisphere (excluding increasing in Antarctic). The ZTD trends are almost symmetrically decreasing with increasing altitude, while the sum of trends at globally distributed GPS sites is almost zero, possibly reflecting that the secular ZTD variation is in balance at a global scale. The annual variation of ZTD ranges from 25 to 75 mm depending on the site, and the mean amplitude is about 50 mm at most sites. The annual variation amplitudes of ZTD at the IGS sites near Oceanic coasts are generally larger than in the continental inland. Larger amplitudes of annual ZTD variation are mostly found at middle-low latitudes (near 20S° and 40N°), and the smaller amplitudes of annual ZTD variation are located at higher latitudes (e.g. Antarctic and Arctic) and the equator areas. The phase of annual ZTD variation is almost about 60° (about February, summer) in the Southern Hemisphere and at about 240° (about August, summer) in the Northern Hemisphere. The mean amplitude of semiannual ZTD variations is about 10 mm, much smaller than annual variations. The significant semiannual variations with a consistent phase of about 30° (about January) are at above 30°N in the Northern Hemisphere and the amplitudes of semiannual variations in other parts are not significant. In addition, the higher frequency variability (RMS of the ZTD residuals) ranges from 22 to 40 mm of delay, once again primarily due to the wet component. The variability depends on altitude of the station. Inland stations tend to have lower variability and stations at ocean and coasts have higher variability. This is because these stations in particular are located in a region well known for large abrupt changes in the weather, such as Indian, West Pacific and West Atlantic oceans.

In addition, the responses of the key ionospheric F2-layer parameters (NmF2 and hmF2) to the 20 November 2003 super storm have been studied using the GPS ionospheric tomography technique over South Korea. A strong increase of NmF2 during this storm has been found, accompanied by a significant hmF2 uplift, which is also confirmed by independent ionosonde observations. The uplift of the F2 layer is mainly associated with a strong eastward electric field. The increase of electron density in the F2-layer peak depends mainly on the molecular nitrogen concentration [N2] with some contribution from molecular oxygen concentration [O2], while the production rate depends on the atomic oxygen concentration [O]. However, the O/N2 ratio from the GUVI instrument on board the TIMED satellite shows no significant change during this geomagnetic storm. It suggests that the increase in NmF2 is not caused by changes in neutral composition, but is related to other 376 Global Navigation Satellite Systems – Signal, Theory and Applications

temperature and correlated with the relative humidity. The mean correlation coefficient between ZWD and surface temperature is about 0.80 and the correlation coefficient between ZTD and ZWD is about 0.95, reflecting a good agreement between ZTD and ZWD variations. Therefore, the seasonal cycles of the ZTD are due primarily to the wet component (ZWD), especially in the surface temperature, even though the wet delay is only

The lower mean ZTD values are located at the areas of higher altitude (e.g. Tibet, Asia and Andes Mountain, South America) and higher latitude areas (Antarctica and Arctic), and the higher mean ZTD values are concentrated at the areas of middle-low latitudes. The mean ZTD decreases with increasing altitude at an exponentially decreasing rate due to the atmospheric pressure decreasing with the height increasing. The mean secular ZTD variation trend is about 1.5±0.001 mm/yr at all GPS sites. The secular variations are systematically increasing in most parts of the Northern Hemisphere and decreasing in most parts of the Southern Hemisphere (excluding increasing in Antarctic). The ZTD trends are almost symmetrically decreasing with increasing altitude, while the sum of trends at globally distributed GPS sites is almost zero, possibly reflecting that the secular ZTD variation is in balance at a global scale. The annual variation of ZTD ranges from 25 to 75 mm depending on the site, and the mean amplitude is about 50 mm at most sites. The annual variation amplitudes of ZTD at the IGS sites near Oceanic coasts are generally larger than in the continental inland. Larger amplitudes of annual ZTD variation are mostly found at middle-low latitudes (near 20S° and 40N°), and the smaller amplitudes of annual ZTD variation are located at higher latitudes (e.g. Antarctic and Arctic) and the equator areas. The phase of annual ZTD variation is almost about 60° (about February, summer) in the Southern Hemisphere and at about 240° (about August, summer) in the Northern Hemisphere. The mean amplitude of semiannual ZTD variations is about 10 mm, much smaller than annual variations. The significant semiannual variations with a consistent phase of about 30° (about January) are at above 30°N in the Northern Hemisphere and the amplitudes of semiannual variations in other parts are not significant. In addition, the higher frequency variability (RMS of the ZTD residuals) ranges from 22 to 40 mm of delay, once again primarily due to the wet component. The variability depends on altitude of the station. Inland stations tend to have lower variability and stations at ocean and coasts have higher variability. This is because these stations in particular are located in a region well known for large abrupt changes in the weather, such as Indian, West Pacific and West

In addition, the responses of the key ionospheric F2-layer parameters (NmF2 and hmF2) to the 20 November 2003 super storm have been studied using the GPS ionospheric tomography technique over South Korea. A strong increase of NmF2 during this storm has been found, accompanied by a significant hmF2 uplift, which is also confirmed by independent ionosonde observations. The uplift of the F2 layer is mainly associated with a strong eastward electric field. The increase of electron density in the F2-layer peak depends mainly on the molecular nitrogen concentration [N2] with some contribution from molecular oxygen concentration [O2], while the production rate depends on the atomic oxygen concentration [O]. However, the O/N2 ratio from the GUVI instrument on board the TIMED satellite shows no significant change during this geomagnetic storm. It suggests that the increase in NmF2 is not caused by changes in neutral composition, but is related to other

10% of the total delay (ZTD).

Atlantic oceans.

possible non-chemical effects, such as dynamical changes of vertical ion motions induced by winds and E×B drifts, tides and waves in the mesosphere/lower thermosphere (MLT) region, which can be dynamically coupled upward to generate ionospheric perturbations and oscillations.

Ground-based and space borne GPS observations have been widely used in atmospheric sounding, including sensing tropopsheric precipitable water vapor (PWV), ionospheric total electron content (TEC) and atmospheric profile information (e.g., pressure, temperature, humidity, tropopause and ionospheric electron density). These observations have facilitated greater advancements in meteorology, climatology, numerical weather model, atmospheric science and space weather (e.g., Jin et al., 2007; Jin et al., 2009; Schmidt, et al. 2010). For example, the dual frequency ground GPS array could detect ionospheric response and its processes during large geomagnetic storms (Jin et al., 2008). Meanwhile, ground GPS also observed the plasma bubbles and retrieved reliable propagation characteristics of the depletions without assumptions about the mapping of the depletion along magnetic field lines to large latitudinal distances, comparable with airglow data (Haase et al., 2011).

Additionally, space-borne GPS receivers have proven very successful in making high vertical resolution and global atmospheric measurements using the radio occultation (RO) technique The existing GPS radio occultation (RO) missions have been widely used to estimate the detailed vertical profile information, including pressure, temperature, gravity waves and sporadic E-layers as well as their variation characteristics, particularly six satellites of the Taiwan/US FORMOSAT-3/COSMIC (FORMOsa SATellite mission-3/Constellation Observing System for Meteorology, Ionosphere and Climate) mission with more than 2000 radio occultation profiles per day. Schmidt et al. (2010) was the first to observe upper tropospheric warming and lower stratospheric cooling using GPS RO data (2001-2009). Although a number of progresses in atmospheric and ionospheric sensing have been made using GPS RO missions in the past few years, they still do not satisfy actual requirements for short-time scales and higher temporal-spatial resolution monitoring together with ground GNS observations. For instance, the tropospheric or ionospheric profile information cannot be directly estimated from GPS tomography due to the lack of enough line-of-sight GPS signals passing each grid cell (Jin et al., 2006 and 2008; Nesterov and Kunitsyn, 2011). Moreover, most current RO satellite missions are approaching their end of operations. With the increase of future GNSS satellite constellations and more GNSS RO missions, the goal of improved temporal-spatial resolution will enable more detailed profile information and evolution processes of the atmosphere and ionosphere.

## **5. Acknowledgement**

This work was supported by the Shanghai Pujiang Talent Program, National Natural Science Foundation of China (Grant No.11043008) and Key Direction Project of Chinese Academy of Sciences (Grant No. KJCX2- EW-T03).

## **6. References**

Austen, J. R.; S. J. Franke, and C. H. Liu, Application of computerized tomography techniques to ionospheric research, in Radio Beacon contributions to the study of ionization and dynamics of the ionosphere and to corrections to geodesy and

GNSS Atmospheric and Ionospheric Sounding – Methods and Results 379

Jin SG, Wang JL, Zhang HP, Zhu WY (2004), Real-time monitoring and prediction of the

Jin S.G., and P.H. Park (2005), A new precision improvement of zenith tropospheric delay

Jin SG, Park J, Wang J, Choi B, (2006), Electron density profiles derived from ground-based

Jin SG, Park J, Cho J, Park P (2007), Seasonal variability of GPS-derived Zenith Tropospheric Delay (1994-2006) and climate implications. J Geophys Res 112: D09110,

Jin, S.G., O. Luo, and P. Park (2008), GPS observations of the ionospheric F2-layer behavior

Jin, S.G., O. Luo, and S. Gleason (2009), Characterization of diurnal cycles in ZTD from a

King R.W. and Y. Bock (1999), Documentation for the GAMIT GPS Analysis Software, Mass.

Lei J, Liu L, Wan W, Zhang S (2004) Modeling the behavior of ionosphere above Millstone

Manuel H., J. Juan, J. Sanz, et al. (2001), A New Strategy for Real-Time Integrated Water

Nesterov, I.A., V.E., Kunitsyn, GNSS radio tomography of the ionosphere: The problem

Niell A.E. (1996), Global Mapping Functions for the Atmospheric Delay at Radio

Raymund TD, Austen JR, Franke SJ (1990) Application of computerized tomography to the

Rocken, C. Analysis and validation of GPS/MET data in the neutral atmosphere. Journal of

Saastamoinent J. (1973), Contribution to the theory of atmospheric refraction, Bulletin

Schmidt, T., Wickert, J., and Haser, A. (2010), Variability of the upper troposphere and lower

temperatures, Adv. Space Res., 46(2),150-161,doi:10.1016/j.asr.2010.01.021 Syndergaard, S. On the Ionosphere Calibration in GPS Radio Occultation Measurements.

Tsai LC, Liu CH, Tsai WH, Liu CT (2002) Tomographic imaging of the ionosphere using the GPS/MET and NNSS data. J. Atmos. Sol. Terr. Phys. 64, 2003-2011. Tregoning P., R. Boers and D. O'Brien (1998) Accuracy of Absolute Precipitable Water Vapor

stratosphere observed with GPS radio occultation bending angles and

Estimates from GPS Observations, Journal of Geophysical Research, 103(28), 701-

with essentially incomplete data, Advances in Space Research,

Wavelengths, Journal of Geophysical Research, 101(B2), 3227-3246.

investigation of ionospheric structures. Radio Science 25: 771-789.

Geophysical. Research, 102: 29849-29866, 1997.

during the 20th November 2003 geomagnetic storm over South Korea, J. Geod.,

decade of global GPS observations, J. Geod., 83(6), 537-545, doi: 10.1007/s00190-

Hill during the September 21-27, 1998 storm. Journal of Atmospheric and Solar-

Vapor Determination in WADGOPS Networks, Geophy. Res. Lett., 28(17), 3267-

estimates by GPS, Current Science, 89 (6): 997-1000.

82(12), 883-892, doi: 10.1007/s00190-008-0217-x.

Inst. of Technol., Cambridge Mass.

Terrestrial Physics 66: 1093-1102.

doi: 10.1016/j.asr.2010.11.034, 2011.

Geodesique, 107, 13-34.

Radio Science, 35(3): 865-883, 2000.

GPS observations. Journal of navigation 59(3): 395-401.

Astrophysics 28(3): 331-337.

doi: 10.1029/2006JD007772.

008-0264-3.

3270.

710.

total ionospheric electron content (TEC) by means of GPS. Chinese Astronomy and

technical workshop, A. Tauriainen, Eds. pp. 25-35, Ouluensis Universitas, Oulu, Finland, 1986.


378 Global Navigation Satellite Systems – Signal, Theory and Applications

Beutler G., M. Rothacher, S. Schaer, T.A. Springer, J. Kouba, R.E. Neilan (1999), The

Bevis M., S. Businger, S. Chiswell, et al. (1994), GPS Meteorology: Mapping Zenith Wet

Byun, S.H., Y. Bar-Sever, and G.Gendt (2005), The new tropospheric product of the

Dai A, Wang J (1999) Diurnal and Semidiurnal Tides in Global Surface Pressure Fields.

Dai A, Wang J, Ware RH, Van Hove T (2002) Diurnal variation in water vapor over North

Davis, J.L., T.A. Herring, I. Shapiro, A. Rogers, and G.Elgered (1985), Geodesy by radio

Deblonde G., S. Macpherson, Y. Mireault , et al.(2005), Evaluation of GPS precipitable water

Elgered, G., et al. (1997), Measuring Regional Atmospheric Water Vapor Using the Swedish Permanent GPS Network, Geophysical Research Letters, 24, 2663-2666. Feng, D., Herman B.M. Remotely sensing the earth's atmosphere using the Global

Foster J, Coster A, Erickson P, Rich F, Sandel B (2004) Stormtime observations of the flux of plasmaspheric ions to the dayside cusp/magnetopause, Geophys Res Lett 31:

Goncharenko LP, Foster J, Coster A, Huang C, Aponte N, Paxton L(2007) Observations of

Gordon R, Bender R, Therman G (1970) Algebraic Reconstruction Techniques (ART) for

Haase, J., M. Ge, H.Vedel and E. Calais (2003),Accuracy and Variability of GPS Tropospheric

Hagemann S., L. Bengtsson and G. Gendt (2003), On the determination of atmospheric water

Hajj, G.A., Romans, L.L. Ionospheric Electron Density Profiles Obtained with the Global

Haase, J.S., T. Dautermann,, M.J. Taylor, N. Chapagain, E. Calais,, D. Pautet, Propagation of

vapor from GPS measurements, J. Geophys. Res., 108 (D21): 4678,

Delays Onto Precipitable Water, J. App. Meteor., 33, 379-386.

J Geophys Res 107(D10): 4090, doi:10.1029/ 2001JD000642.

Sciences, Adv. Space Res. 23(4) 631-635, 1999.

of Navig., Long Beach, Calif.

length, Radio Sci., 20(6), 1593-1607.

and Oceanic Technology, 16, 989-1002, 1999.

Applied Meteorology,42(11): 1547-1568.

doi: 10.1016/j.asr.2010.09.025, 2011.

J Atmos Sci 56, 3874-3891.

doi:10.1029/2004GL020082.

doi:10.1029/2002JD003235.

175-190, 1998.

1272.

471-481

Finland, 1986.

technical workshop, A. Tauriainen, Eds. pp. 25-35, Ouluensis Universitas, Oulu,

International GPS Service (IGS): An Interdisciplinary Service in Support of Earth

International GNSS service, paper presentaed at 2005 ION GNSS conference, Inst.

America and its implications for sampling errors in radiosonde humidity.

interferometryL Effects of atmospheric modeling errors on esimatetes of baseline

over Canada and the IGS network, Journal of Applied Meteorology, 44 (1): 153-166.

Positioning System (GPS)- The GPS/MET data analysis. Journal of Atmospheric

a positive storm phase on September 10, 2005, J Atmos Solar Terr Phys 69: 1253-

three Dimensional Electron Micoscopy and X-ray Photography. J Theor.Biol 29:

Delay Measurements of Water Vapor in the Western Mediterranean, Journal of

Positioning System: Results from the GPS/MET Experiment. Radio Science, 33(1):

plasma bubbles observed in Brazil from GPS and airglow data. Adv. Space Res.


**16** 

*Neustrelitz, Germany* 

**Ionospheric Propagation Effects on GNSS** 

*German Aerospace Center (DLR), Institute of Communications and Navigation* 

The ionosphere is the ionized part of the earth's atmosphere lying between about 50 km and several earth radii (Davies, 1990) whereas the upper part above about 1000 km height up to the plasmapause is usually called the plasmasphere. Solar extreme ultraviolet (EUV) radiation at wave lengths < 130 nm significantly ionizes the earth's neutral gas. In addition to photoionisation by electromagnetic radiation also energetic particles from the solar wind and cosmic rays contribute to the ionization. The ionized plasma can affect radio wave propagation in various ways modifying characteristic wave parameters such as amplitude, phase or polarization (Budden, 1985; Davies, 1990). The interaction of the radio wave with the ionospheric plasma is one of the main reasons for the limited accuracy and vulnerability

A trans-ionospheric radio wave propagating through the plasma experiences a propagation delay / phase advance of the signal causing a travel distance or time larger / smaller than the real one. The reason of the propagation delay can be realized considering the nature of the refractive index which depends on the density of the ionospheric plasma. The refractive index (n ≠ 1) of the ionosphere is not equal to that of free space (n = 1). This causes the propagation speed of radio signals to differ from that in free space. Additionally, spatial gradients in the refractive index cause a curvature of the propagation path. Both effects lead in sum to a delay / phase advance of satellite navigation signals in comparison to a free

The variability of the ionospheric impact is much larger compared to that of the troposphere. The ionospheric range error varies from a few meters to many tens of meters at the zenith, whereas the tropospheric range error varies between two to three meters at the zenith (Klobuchar, 1996). The daily variation of the ionospheric range error can be up to one

After removal of the Selective Availability (SA, i.e., dithering of the satellite clock to deny full system accuracy) in 2000, ionosphere becomes the single largest error source for Global Navigation Satellite Systems (GNSS) users, especially for high-accuracy (centimeter millimeter) applications like the Precise Point Positioning (PPP) and Real Time Kinematic (RTK) positioning. Fortunately, the ionosphere is a dispersive medium with respect to the

**1. Introduction** 

space propagation.

order of magnitude (Klobuchar, 1996).

in satellite based positioning or time estimation.

**Signals and New Correction Approaches** 

M. Mainul Hoque and Norbert Jakowski


## **Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches**

M. Mainul Hoque and Norbert Jakowski *German Aerospace Center (DLR), Institute of Communications and Navigation Neustrelitz,* 

*Germany* 

## **1. Introduction**

380 Global Navigation Satellite Systems – Signal, Theory and Applications

Wagner, C.. Kloko, J. The value of ocean reflections of GPS signals to enhance satellite

Westwater, E. R. (1993), Ground-based microwave remote sensing of meteorological

Yin P, Mitchell CN, Spencer PS, Foster JC (2004) Ionospheric electron concentration imaging

2003.

Ed., J. Wiley and Sons, Inc., 145-213.

L12806, doi: 10.1029/2004GL019899.

altimetry: data distribution and error analysis. Journal of Geodesy, 74, 128-138,

variables. Atmospheric Remote Sensing by Microwave Radiometry, M. A. Janssen,

using GPS over the USA during the storm of July 2000. Geophys Res Lett 31:

The ionosphere is the ionized part of the earth's atmosphere lying between about 50 km and several earth radii (Davies, 1990) whereas the upper part above about 1000 km height up to the plasmapause is usually called the plasmasphere. Solar extreme ultraviolet (EUV) radiation at wave lengths < 130 nm significantly ionizes the earth's neutral gas. In addition to photoionisation by electromagnetic radiation also energetic particles from the solar wind and cosmic rays contribute to the ionization. The ionized plasma can affect radio wave propagation in various ways modifying characteristic wave parameters such as amplitude, phase or polarization (Budden, 1985; Davies, 1990). The interaction of the radio wave with the ionospheric plasma is one of the main reasons for the limited accuracy and vulnerability in satellite based positioning or time estimation.

A trans-ionospheric radio wave propagating through the plasma experiences a propagation delay / phase advance of the signal causing a travel distance or time larger / smaller than the real one. The reason of the propagation delay can be realized considering the nature of the refractive index which depends on the density of the ionospheric plasma. The refractive index (n ≠ 1) of the ionosphere is not equal to that of free space (n = 1). This causes the propagation speed of radio signals to differ from that in free space. Additionally, spatial gradients in the refractive index cause a curvature of the propagation path. Both effects lead in sum to a delay / phase advance of satellite navigation signals in comparison to a free space propagation.

The variability of the ionospheric impact is much larger compared to that of the troposphere. The ionospheric range error varies from a few meters to many tens of meters at the zenith, whereas the tropospheric range error varies between two to three meters at the zenith (Klobuchar, 1996). The daily variation of the ionospheric range error can be up to one order of magnitude (Klobuchar, 1996).

After removal of the Selective Availability (SA, i.e., dithering of the satellite clock to deny full system accuracy) in 2000, ionosphere becomes the single largest error source for Global Navigation Satellite Systems (GNSS) users, especially for high-accuracy (centimeter millimeter) applications like the Precise Point Positioning (PPP) and Real Time Kinematic (RTK) positioning. Fortunately, the ionosphere is a dispersive medium with respect to the

Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches 383

frequencies (> 100 MHz), the refractive index mainly depends on the electron density, the strength and direction of the geomagnetic field in relation to the ray path. Thus, the spatial distribution of the electron density along the ray path and corresponding geomagnetic field

For high frequency (HF) radio waves with frequencies *f* > 100 MHz the phase refractive index *n* can be derived from the Appleton – Hartree formula as (Appleton, 1932; Bassiri &

1 1 cos

=− ± − + + <sup>⎢</sup> <sup>⎥</sup> <sup>⎢</sup> <sup>⎥</sup> <sup>⎣</sup> <sup>⎦</sup>

where *fp* is the plasma frequency, *fg* is the gyro frequency, *ε*0 is the free space permittivity, *B* is the geomagnetic induction, Θ is the angle between the wave propagation direction and the geomagnetic field vector *B*, and *e*, *ne*, *m* are the electron charge, density and mass, respectively. The wave with the upper (+) sign in Eq. (1) is called the ordinary wave and is left-hand circularly polarized, whereas the wave with the lower (-) sign is called the extraordinary wave and is right-hand circularly polarized (Hartmann & Leitinger, 1984). The GPS signals are transmitted in right-hand circular polarization (Parkinson & Gilbert,

Equation (1) indicates that the phase refractive index is less than the unity resulting in a phase velocity that is greater than the speed of light in vacuum (i.e., phase advance). Therefore, the integration of *n* along a signal path gives a measure of the range / travel time between a receiver and a satellite that is smaller than the geometric distance / travel time in

To compute group delay measurements, the group refractive index *ngr* should be considered. The expression for *ngr* can be determined by the relationship *ngr* = *n*+*f*(*dn*/*df*).

= + + ++ <sup>⎢</sup> <sup>⎥</sup> <sup>⎢</sup> <sup>⎥</sup> <sup>⎣</sup> <sup>⎦</sup>

Equation (2) indicates that the group refractive index is greater than the unity resulting in a group velocity that is less than the speed of light. Therefore, the integration of *ngr* along a signal path gives a measure of the range / travel time that is greater than the geometric distance / travel time in the vacuum. Therefore, when GNSS signals propagate through the ionosphere, the carrier-phase experiences an advance and the code experiences a group delay. The carrier-phase pseudoranges are measured too short and the code pseudoranges are measured too long compared to the geometric range between a satellite and a receiver.

2 2 2 2

Θ

23 4 cos 3 1 1 cos

2 4 2 *p pg p p gr g f ff f f n f ff f*

2 2 2 2

234 cos

*f ff f f n f fff* Θ

224 2 *p pg p p*

( )

( )

Θ

2 2

∓ (2)

⎡ ⎤

Θ

(1)

2 2

⎡ ⎤

*g*

relationships determine the ionospheric impact on the electromagnetic wave.

**2.1 Ionospheric refractive index** 

Hajj, 1993)

in which

1983).

the vacuum.

( ) 22 2 *p e* / 4 <sup>0</sup> *f ne m* = π ε

/ 2( ) *gf* = *eB m*π

radio wave; therefore, the magnitude of the ionospheric delay depends on the signal frequency. The advantage is that an elimination of the major part of the ionospheric refraction through a linear combination of dual-frequency observables is possible. However, inhomogeneous plasma distribution and anisotropy cause higher order nonlinear effects which are not removed in this linear approach. Mainly the second and third order ionospheric terms (in the expansion of the refractive index) and errors due to bending of the signal remain uncorrected. They can be several tens of centimeters of range error at low elevation angles and during high solar activity conditions.

Brunner & Gu (1991) were pioneers to compute higher order ionospheric effects and developing correction for them. Since then higher order ionospheric effects have been studied by different authors during last decades, e.g., Bassiri & Hajj (1993), Jakowski et al. (1994), Strangeways & Ioannides (2002), Kedar et al. (2003), Fritsche et al. (2005), Hawarey et al. (2005), Hoque & Jakowski (2006, 2007, 2008, 2010b), Hernández-Pajares et al. (2007), Kim & Tinin (2007, 2011), Datta-Barua et al. (2008), Morton et al. (2009), Moore & Morton (2011). The above literature review shows that higher order ionospheric terms are less than 1% of the first order term at GNSS frequencies. Hernández-Pajares et al. (2007) found submillimeter level shifting in receiver positions along southward direction for low latitude receivers and northward direction for high latitude receivers due to the second order term correction. Fritsche et al. (2005) found centimeter level correction in GPS satellite positions considering higher order ionospheric terms. Elizabeth et al. (2010) investigated the impacts of the bending terms described by Hoque & Jakowski (2008) on a Global Positioning System (GPS) network of ground receivers. They found the bending correction for the dualfrequency linear GPS L1-L2 combination to exceed the 3 mm level in the equatorial region. Kim & Tinin (2011) found that the systematic residual ionospheric errors can be significantly reduced (under certain ionospheric conditions) through triple frequency combinations. All these studies were conducted to compute higher order ionospheric effects on GNSS signals for ground-based reception. Recently Hoque & Jakowski (2010b, 2011) investigated the ionospheric impact on GPS occultation signals received onboard Low Earth Orbiting (LEO) CHAMP (CHAllenging Minisatellite Payload) satellite.

In this chapter, the first and higher order ionospheric propagation effects on GNSS signals are described and their estimates are given at different level of ionospheric ionization. Multi-frequency ionosphere-free and geometry-free solutions are studied and residual terms in the ionosphere-free solutions are computed. Different correction approaches are discussed for the second and third order terms, and ray path bending correction. Additionally, we have proposed new approaches for correcting straight line of sight (LoS) propagation assumption error, i.e., ray path bending error for ground based GNSS positioning. We have modelled the excess path length of the signal in addition to the LoS path length and the total electron content (TEC) difference between a curved and LoS paths as functions of signal frequency, ionospheric parameters such as TEC and TEC derivative with respect to the elevation angle. We have found that using the TEC derivative in addition to the TEC information we can improve the existing correction results.

## **2. Ionospheric propagation effects**

Quantitatively, the propagation of a radio wave through the ionospheric plasma is described by the refractive index of the ionosphere (Appleton-Hartree formula). At high frequencies (> 100 MHz), the refractive index mainly depends on the electron density, the strength and direction of the geomagnetic field in relation to the ray path. Thus, the spatial distribution of the electron density along the ray path and corresponding geomagnetic field relationships determine the ionospheric impact on the electromagnetic wave.

#### **2.1 Ionospheric refractive index**

For high frequency (HF) radio waves with frequencies *f* > 100 MHz the phase refractive index *n* can be derived from the Appleton – Hartree formula as (Appleton, 1932; Bassiri & Hajj, 1993)

$$n = 1 - \frac{f\_p^2}{2f^2} \pm \frac{f\_p^2 f\_g \cos\Theta}{2f^3} - \frac{f\_p^2}{4f^4} \left[\frac{f\_p^2}{2} + f\_g^2 \left(1 + \cos^2\Theta\right)\right] \tag{1}$$

in which

382 Global Navigation Satellite Systems – Signal, Theory and Applications

radio wave; therefore, the magnitude of the ionospheric delay depends on the signal frequency. The advantage is that an elimination of the major part of the ionospheric refraction through a linear combination of dual-frequency observables is possible. However, inhomogeneous plasma distribution and anisotropy cause higher order nonlinear effects which are not removed in this linear approach. Mainly the second and third order ionospheric terms (in the expansion of the refractive index) and errors due to bending of the signal remain uncorrected. They can be several tens of centimeters of range error at low

Brunner & Gu (1991) were pioneers to compute higher order ionospheric effects and developing correction for them. Since then higher order ionospheric effects have been studied by different authors during last decades, e.g., Bassiri & Hajj (1993), Jakowski et al. (1994), Strangeways & Ioannides (2002), Kedar et al. (2003), Fritsche et al. (2005), Hawarey et al. (2005), Hoque & Jakowski (2006, 2007, 2008, 2010b), Hernández-Pajares et al. (2007), Kim & Tinin (2007, 2011), Datta-Barua et al. (2008), Morton et al. (2009), Moore & Morton (2011). The above literature review shows that higher order ionospheric terms are less than 1% of the first order term at GNSS frequencies. Hernández-Pajares et al. (2007) found submillimeter level shifting in receiver positions along southward direction for low latitude receivers and northward direction for high latitude receivers due to the second order term correction. Fritsche et al. (2005) found centimeter level correction in GPS satellite positions considering higher order ionospheric terms. Elizabeth et al. (2010) investigated the impacts of the bending terms described by Hoque & Jakowski (2008) on a Global Positioning System (GPS) network of ground receivers. They found the bending correction for the dualfrequency linear GPS L1-L2 combination to exceed the 3 mm level in the equatorial region. Kim & Tinin (2011) found that the systematic residual ionospheric errors can be significantly reduced (under certain ionospheric conditions) through triple frequency combinations. All these studies were conducted to compute higher order ionospheric effects on GNSS signals for ground-based reception. Recently Hoque & Jakowski (2010b, 2011) investigated the ionospheric impact on GPS occultation signals received onboard Low Earth Orbiting (LEO)

In this chapter, the first and higher order ionospheric propagation effects on GNSS signals are described and their estimates are given at different level of ionospheric ionization. Multi-frequency ionosphere-free and geometry-free solutions are studied and residual terms in the ionosphere-free solutions are computed. Different correction approaches are discussed for the second and third order terms, and ray path bending correction. Additionally, we have proposed new approaches for correcting straight line of sight (LoS) propagation assumption error, i.e., ray path bending error for ground based GNSS positioning. We have modelled the excess path length of the signal in addition to the LoS path length and the total electron content (TEC) difference between a curved and LoS paths as functions of signal frequency, ionospheric parameters such as TEC and TEC derivative with respect to the elevation angle. We have found that using the TEC derivative in addition

Quantitatively, the propagation of a radio wave through the ionospheric plasma is described by the refractive index of the ionosphere (Appleton-Hartree formula). At high

elevation angles and during high solar activity conditions.

CHAMP (CHAllenging Minisatellite Payload) satellite.

**2. Ionospheric propagation effects** 

to the TEC information we can improve the existing correction results.

$$\begin{aligned} f\_p^2 &= n\_e e^2 \left/ \left(4\pi^2 \varepsilon\_0 m\right)\right. \\ f\_\S &= eB \left/ \left(2\pi m\right)\right. \end{aligned}$$

where *fp* is the plasma frequency, *fg* is the gyro frequency, *ε*0 is the free space permittivity, *B* is the geomagnetic induction, Θ is the angle between the wave propagation direction and the geomagnetic field vector *B*, and *e*, *ne*, *m* are the electron charge, density and mass, respectively. The wave with the upper (+) sign in Eq. (1) is called the ordinary wave and is left-hand circularly polarized, whereas the wave with the lower (-) sign is called the extraordinary wave and is right-hand circularly polarized (Hartmann & Leitinger, 1984). The GPS signals are transmitted in right-hand circular polarization (Parkinson & Gilbert, 1983).

Equation (1) indicates that the phase refractive index is less than the unity resulting in a phase velocity that is greater than the speed of light in vacuum (i.e., phase advance). Therefore, the integration of *n* along a signal path gives a measure of the range / travel time between a receiver and a satellite that is smaller than the geometric distance / travel time in the vacuum.

To compute group delay measurements, the group refractive index *ngr* should be considered. The expression for *ngr* can be determined by the relationship *ngr* = *n*+*f*(*dn*/*df*).

$$m\_{gr} = 1 + \frac{f\_p^2}{2f^2} \mp \frac{f\_p^2 f\_g \cos\Theta}{f^3} + \frac{3f\_p^2}{4f^4} \left[\frac{f\_p^2}{2} + f\_g^2 \left(1 + \cos^2\Theta\right)\right] \tag{2}$$

Equation (2) indicates that the group refractive index is greater than the unity resulting in a group velocity that is less than the speed of light. Therefore, the integration of *ngr* along a signal path gives a measure of the range / travel time that is greater than the geometric distance / travel time in the vacuum. Therefore, when GNSS signals propagate through the ionosphere, the carrier-phase experiences an advance and the code experiences a group delay. The carrier-phase pseudoranges are measured too short and the code pseudoranges are measured too long compared to the geometric range between a satellite and a receiver.

Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches 385

ray paths through the ionosphere on their way to a receiver and thus the TEC along a *f*1 path will be different from that along a *f*2 path and also from that along the LoS path. Considering this, TEC in Eq. (7) is separated into *TECLoS* and Δ*TECbend* where *TECLoS* is the TEC along the LoS and Δ*TECbend* is the difference between TECs along a curved path and the LoS path. The term Δ*TECbend* represents TEC contribution due to ray path bending only, i.e., the second

The observables are travel time or ranges which are deduced from measured time or phase differences based on a comparison between received signals and receiver generated signals. Thus, the ranges are biased by satellite and receiver clock errors, instrumental biases and atmospheric effects, and therefore, called pseudoranges. The code pseudorange (Ψ) and carrier-phase pseudorange (Φ) at a selected frequency can be described by observation

=+ − + + + + + +

=+ − − + + + + + +

*c dt dT d d d dq dQ N* ( ) *I A MP* ( )

where *ρ* is the geometric distance between a satellite and a receiver, *c* is the velocity of light, *dt* and *dT* are the satellite and receiver clock errors, respectively, *dI* and *dIgr* are the ionospheric effects on carrier-phase and code pseudoranges, respectively, *dA* is the atmospheric (tropospheric delay) effect, (*dMP*)*Ψ* and (*dMP*)*Φ* are multipath effects on code and carrier-phase pseudoranges, respectively, *dq* and *dQ* are the instrumental biases of the satellite and the receiver, respectively, *λ* is the carrier wavelength, *N* is the integer carrierphase ambiguity, and *εΨ* and *εΦ* are the rest errors. The carrier-phase pseudorange is expressed in units of length (meters) instead of cycles. However, it can be expressed in

For simplicity we confine our interest to only ionospheric effects. Thus, the code and carrier-

 ρ

where *f* is the signal frequency. In case of GPS L1, L2 and L5 signals *f* = 1575.42, 1227.6 and 1176.45 MHz, respectively. To take into account the ray path bending on observables, the

As already mentioned, ionosphere is a dispersive medium, i.e., the ionospheric propagation delay is frequency dependent. Therefore, one very popular way to get rid of ionospheric

 ρ

*c dt dT d d d dq dQ* ( ) *Igr A MP* ( )

Ψ

Φ

234 *len len Igr I I <sup>p</sup> <sup>q</sup> <sup>u</sup> d d <sup>d</sup>*

=+ + =+ + + + (12)

=− + =− − − + (13)

*fff*

<sup>234</sup> 2 3 *len len I I I <sup>p</sup> <sup>q</sup> <sup>u</sup> d d <sup>d</sup>*

*fff*

 ε

 λ Ψ

εΦ

(10)

(11)

and third order terms are not considered in TEC estimation by Eq. (7).

cycles dividing simply by the signal's wave length (*λ* meter/cycle).

Ψρ

Φρ

*<sup>I</sup> d* is introduced in Eqs. (12) and (13).

**2.4.1 First order ionosphere-free combination** 

**2.4 Multi-frequency combinations** 

**2.3 Ionospheric effects on GNSS observables** 

equations in units of length as

Ψ

Φ

phase pseudoranges can be simplified as

term *len*

ρ

ρ

Considering ionospheric refraction the geometric distance (Euclidean line) or true range *ρ* between a transmitting satellite *S* and a receiver *R* can be written in units of length as

$$\rho = L + \int\_{S}^{R} (1 - n) ds - d\_{l}^{len} \tag{3}$$

where the optical distance *R S L nds* <sup>=</sup> ∫ is the line integral of the refractive index between the

satellite and the receiver along the ray path, ( ) 1 *R S* <sup>−</sup> *n ds* ∫ is the ionospheric group delay and *len*

*<sup>I</sup> d* is the excess path length of the signal in addition to the geometric path length caused by the ray path bending and defined by

$$d\_I^{len} = \bigcap\_{S}^{R} \text{ds} - \rho \tag{4}$$

where *R S ds* ∫ is the curved path length in the vacuum. The travel time of the signal can be computed dividing the expression of *ρ* (Eq. 3) simply by the speed of light.

#### **2.2 Group delay and phase advance**

Assuming a right hand circularly polarized signal, the ionospheric group delay *dIgr* and carrier phase advance *dI* can be written in units of length as (using Eqs. (1) and (2))

$$d\_{lgr} = d\_{lgr}^{(1)} + d\_{lgr}^{(2)} + d\_{lgr}^{(3)} = \underset{\mathcal{S}}{\int} (n\_{\mathcal{S}^r} - 1) \text{ds} = \frac{p}{f^2} + \frac{q}{f^3} + \frac{u}{f^4} \tag{5}$$

$$d\_I = d\_l^{(1)} + d\_l^{(2)} + d\_l^{(3)} = \underset{S}{\int}{(1 - n)} ds = \frac{p}{f^2} + \frac{q}{2f^3} + \frac{u}{3f^4} \tag{6}$$

$$p = K \int n\_c ds = K \cdot TEC = K \left( TEC\_{LoS} + \mathcal{A}TEC\_{bend} \right) \tag{7}$$

$$q = 2.2566 \times 10^{12} \int n\_e B \cos \Theta \cdot ds \tag{8}$$

$$
\mu = 2437 \int n\_e^2 ds + 4.74 \times 10^{22} \int n\_e B^2 (1 + \cos^2 \Theta) ds\tag{9}
$$

where *K* = *e*2/(8π<sup>2</sup>*ε*0*m*) = 40.3 m3s-2, the integration of *ne* along signal paths ∫*n ds <sup>e</sup>* is called the total electron content TEC and measured in TEC units (1 TECU = 1016 electrons/m2). The terms ( ) <sup>1</sup> *Igr <sup>d</sup>* / ( ) <sup>1</sup> *<sup>I</sup> <sup>d</sup>* , ( ) <sup>2</sup> *Igr <sup>d</sup>* / ( ) <sup>2</sup> *<sup>I</sup> <sup>d</sup>* and ( ) <sup>3</sup> *Igr <sup>d</sup>* / ( ) <sup>3</sup> *<sup>I</sup> d* in Eq. (5) / (6) are the first, second and third order ionospheric group delays / phase advances, respectively. Due to the dispersive nature of the ionosphere, satellite signals transmitted on different frequencies travel along different ray paths through the ionosphere on their way to a receiver and thus the TEC along a *f*1 path will be different from that along a *f*2 path and also from that along the LoS path. Considering this, TEC in Eq. (7) is separated into *TECLoS* and Δ*TECbend* where *TECLoS* is the TEC along the LoS and Δ*TECbend* is the difference between TECs along a curved path and the LoS path. The term Δ*TECbend* represents TEC contribution due to ray path bending only, i.e., the second and third order terms are not considered in TEC estimation by Eq. (7).

#### **2.3 Ionospheric effects on GNSS observables**

384 Global Navigation Satellite Systems – Signal, Theory and Applications

Considering ionospheric refraction the geometric distance (Euclidean line) or true range *ρ* between a transmitting satellite *S* and a receiver *R* can be written in units of length as

> (1 ) *R*

> > *R*

*S*

*<sup>I</sup> d* is the excess path length of the signal in addition to the geometric path length caused

*R len I*

*S d ds* = −

Assuming a right hand circularly polarized signal, the ionospheric group delay *dIgr* and

*R*

*S p q <sup>u</sup> d d d d n ds*

*R*

*S p q <sup>u</sup> d d d d n ds*

*p K n ds K TEC K TEC TEC* = =⋅ = + *e L* ( *oS*

<sup>12</sup> 2.2566 10 cos *<sup>e</sup> <sup>q</sup>* <sup>=</sup> × ⋅ ∫*n B ds*

2 222 2 2437 4.74 10 (1 cos ) *u n ds* = +× + *e e n B*

the total electron content TEC and measured in TEC units (1 TECU = 1016 electrons/m2). The

order ionospheric group delays / phase advances, respectively. Due to the dispersive nature of the ionosphere, satellite signals transmitted on different frequencies travel along different

*Igr <sup>d</sup>* / ( ) <sup>3</sup>

carrier phase advance *dI* can be written in units of length as (using Eqs. (1) and (2))

() () () ( ) <sup>123</sup>

() () () ( ) <sup>123</sup>

*Igr Igr Igr Igr gr*

*I III*

computed dividing the expression of *ρ* (Eq. 3) simply by the speed of light.

ρ

*ds* ∫ is the curved path length in the vacuum. The travel time of the signal can be

*S*

ρ

*R*

*S*

satellite and the receiver along the ray path, ( ) 1

by the ray path bending and defined by

**2.2 Group delay and phase advance** 

where the optical distance

*len*

where

*R*

*S*

where *K* = *e*2/(8

*Igr <sup>d</sup>* / ( ) <sup>1</sup>

terms ( ) <sup>1</sup>

π

*<sup>I</sup> <sup>d</sup>* , ( ) <sup>2</sup>

*Igr <sup>d</sup>* / ( ) <sup>2</sup>

*<sup>I</sup> <sup>d</sup>* and ( ) <sup>3</sup>

*len I*

=+ − − *L n ds d* ∫ (3)

<sup>−</sup> *n ds* ∫ is the ionospheric group delay and

∫ (4)

*L nds* <sup>=</sup> ∫ is the line integral of the refractive index between the

<sup>234</sup> 1

<sup>234</sup> 1

Δ*bend* ) ∫ (7)

Θ

<sup>2</sup>*ε*0*m*) = 40.3 m3s-2, the integration of *ne* along signal paths ∫*n ds <sup>e</sup>* is called

*fff* = + + = − =++ ∫ (5)

2 3

Θ*ds* ∫ ∫ (9)

*<sup>I</sup> d* in Eq. (5) / (6) are the first, second and third

(8)

*fff* = + + = − =+ + ∫ (6) The observables are travel time or ranges which are deduced from measured time or phase differences based on a comparison between received signals and receiver generated signals. Thus, the ranges are biased by satellite and receiver clock errors, instrumental biases and atmospheric effects, and therefore, called pseudoranges. The code pseudorange (Ψ) and carrier-phase pseudorange (Φ) at a selected frequency can be described by observation equations in units of length as

$$\Psi' = \rho + c\left(dt - dT\right) + d\_{\rm lgr} + d\_A + \left(d\_{\rm MP}\right)\_{\rm pr} + d\eta + dQ + \varepsilon\_{\Psi} \tag{10}$$

$$
\mathcal{Q} = \mathcal{\rho} + \mathcal{c}\left(dt - dT\right) - d\_I + d\_A + \left(d\_{MP}\right)\_{\mathcal{\Phi}} + dq + d\mathcal{Q} + N\mathcal{X} + \varepsilon\_{\Phi} \tag{11}
$$

where *ρ* is the geometric distance between a satellite and a receiver, *c* is the velocity of light, *dt* and *dT* are the satellite and receiver clock errors, respectively, *dI* and *dIgr* are the ionospheric effects on carrier-phase and code pseudoranges, respectively, *dA* is the atmospheric (tropospheric delay) effect, (*dMP*)*Ψ* and (*dMP*)*Φ* are multipath effects on code and carrier-phase pseudoranges, respectively, *dq* and *dQ* are the instrumental biases of the satellite and the receiver, respectively, *λ* is the carrier wavelength, *N* is the integer carrierphase ambiguity, and *εΨ* and *εΦ* are the rest errors. The carrier-phase pseudorange is expressed in units of length (meters) instead of cycles. However, it can be expressed in cycles dividing simply by the signal's wave length (*λ* meter/cycle).

For simplicity we confine our interest to only ionospheric effects. Thus, the code and carrierphase pseudoranges can be simplified as

$$d\Psi = \rho + d\_{lgr} + d\_l^{len} = \rho + \frac{p}{f^2} + \frac{q}{f^3} + \frac{u}{f^4} + d\_l^{len} \tag{12}$$

$$\rho \Phi = \rho - d\_I + d\_I^{len} = \rho - \frac{p}{f^2} - \frac{q}{2f^3} - \frac{u}{3f^4} + d\_I^{len} \tag{13}$$

where *f* is the signal frequency. In case of GPS L1, L2 and L5 signals *f* = 1575.42, 1227.6 and 1176.45 MHz, respectively. To take into account the ray path bending on observables, the term *len <sup>I</sup> d* is introduced in Eqs. (12) and (13).

#### **2.4 Multi-frequency combinations**

#### **2.4.1 First order ionosphere-free combination**

As already mentioned, ionosphere is a dispersive medium, i.e., the ionospheric propagation delay is frequency dependent. Therefore, one very popular way to get rid of ionospheric

Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches 387

Assuming the same measurement noise on each signal, it can be shown that the carrierphase or code noise will be amplified by a factor of 2.98 for the GPS L1-L2 combination, whereas for the L1-L5 and L2-L3 combinations amplification factors are 2.59 and 16.64, respectively (see Hoque & Jakowski, 2010a). The amplification factor is inversely proportional to the separation of combination frequencies. Since the frequency separation is

Since the first order ionospheric effect on carrier-phase and code pseudoranges (see Eq. 12 and 13) is the same in magnitude but opposite in sign, computing the sum of carrier-phase and code pseudoranges would theoretically eliminate the first order ionospheric term in single frequency measurements. However, the resulting observable would inherit the high

Receiving signals on three coherent frequencies will allow triple-frequency combinations to eliminate the first and the second order ionospheric terms. The third order ionospheric term and errors due to ray path bending are not fully removed in this approach. Such a triplefrequency combination can be written as (combining code / carrier-phase pseudoranges Eq.

( ) ( ) ( ) ()( )

( ) ( ) ( )()( )

− − − =− − + ⎣ ⎦ 

− − − =+ + + ⎣ ⎦ 

 ρΔ

 ρΔ

= −− − ⎡ ⎤ ⎣ ⎦ (23)

= −− − <sup>⎡</sup> <sup>⎤</sup> <sup>⎣</sup> <sup>⎦</sup> (25)

⎪ ⎭

11 22 11 33 3

11 22 11 33 3

 Φ

*Af f Bf f s s s <sup>C</sup>*

( ) *TEC tr* ( ) *bend*31 21 *bend* ( ) *bend bend <sup>K</sup> s B TEC TEC A TEC TEC*

( ) ( 2 3 )

( ) ( )( ) 22 22

*len tr s Bfd fd Afd fd <sup>C</sup>*

*f f <sup>A</sup> f f*

*f f <sup>B</sup> f f Cff f f f f*

2 3 3 *tr <sup>u</sup> <sup>f</sup> <sup>f</sup> <sup>s</sup> C ff*

( )( )

= − ++ ⎪

1 2 1 2 1 3 1 3 12 3 1 2 3

33 11 22 11 1 *len len len len*

<sup>⎫</sup> <sup>=</sup> <sup>⎪</sup> <sup>−</sup> <sup>⎪</sup>

<sup>⎪</sup> <sup>=</sup> <sup>⎬</sup> <sup>−</sup> <sup>⎪</sup>

*Af f Bf f s s s <sup>C</sup>*

ρ

( )

( )

ΔΔ

<sup>−</sup> <sup>=</sup> (24)

 Δ

 Δ

*gr tr TEC tr tr len tr RRE*

*tr TEC tr tr len tr RRE*

 Δ

 Δ

(26)

by Eq. (7), for details

(21)

(22)

relatively large for the L1-L5 combination, the amplification factor is the smallest.

code noise and the carrier-phase ambiguity and is therefore, practically not suitable.

**2.4.2 Second order ionosphere-free combination** 

see Hoque & jakowski, 2010a).

Ψ

Φ

1

Δ

In which

 Ψ

> Φ

*C*

Δ

(12) / Eq. (13) measured on *f*1, *f*2 and *f*3 frequencies and substituting

22 22

 Ψ

22 22

 Φ

ΔΔ

3

Δ

⎡ ⎤

⎡ ⎤

<sup>1</sup> <sup>3</sup>

 Ψ

effects is to compute the so called first order ionosphere-free combination of carrier-phase or code pseudoranges measured on two frequencies. However, the second and third order ionospheric terms and errors due to bending of the signal remain uncorrected in this approach. Such a dual-frequency combination can be written in units of length as (combining code / carrier-phase pseudoranges Eq. (12) / Eq. (13) measured on *f*1 and *f*<sup>2</sup> frequencies and substituting *p* by Eq. (7), for details see Hoque & Jakowski, 2008)

$$\frac{f\_1^2}{f\_1^2 - f\_2^2} \Psi\_1 - \frac{f\_2^2}{f\_1^2 - f\_2^2} \Psi\_2 = \rho \underbrace{-\mathbf{A} \mathbf{s}\_{\text{TEC}} - 2\mathbf{A} \mathbf{s}\_2 - 3\mathbf{A} \mathbf{s}\_3 - \mathbf{A} \mathbf{s}\_{\text{len}}}\_{\text{RRE}\_{\text{yr}}} \tag{14}$$

$$\frac{f\_1^2}{f\_1^2 - f\_2^2} \Phi\_1 - \frac{f\_2^2}{f\_1^2 - f\_2^2} \Phi\_2 = \rho \underbrace{+ \mathbf{A} \mathbf{s}\_{\text{TEC}} + \mathbf{A} \mathbf{s}\_2 + \mathbf{A} \mathbf{s}\_3 - \mathbf{A} \mathbf{s}\_{\text{len}}}\_{\text{RRE}} \tag{15}$$

$$\text{As}\_{\text{TEC}} = \frac{\text{K}(\text{TEC}\_2 - \text{TEC}\_1)}{\left(f\_1^2 - f\_2^2\right)} = \frac{\text{K}(\text{A}\text{TEC}\_{bend2} - \text{A}\text{TEC}\_{bend1})}{\left(f\_1^2 - f\_2^2\right)}\tag{16}$$

$$\text{TEC}\_{1,2} = \int n\_e ds = \left( \text{TEC}\_{LoS} + \mathcal{A} \text{TEC}\_{bend1,2} \right) \tag{17}$$

$$\Delta \mathbf{s}\_2 = \frac{q}{2f\_1 f\_2 \left(f\_1 + f\_2\right)}\tag{18}$$

$$\Delta \mathbf{s}\_3 = \frac{\boldsymbol{\mu}}{3f\_1^2 f\_2^2} \tag{19}$$

$$\Delta \mathbf{s}\_{\text{let}} = \frac{d\_2^{\text{len}} f\_2^2 - d\_1^{\text{len}} f\_1^2}{\left(f\_1^2 - f\_2^2\right)} \tag{20}$$

where Ψ1, Ψ2 and Φ1, Φ1 are the measured code and carrier-phase pseudoranges on *f*1 and *f*<sup>2</sup> frequencies, respectively, Δ*s*2 and Δ*s*3 are the dual-frequency second and third order residual terms, respectively. The TEC along a *f*1 path will be different from that along a *f*2 path due to ray path bending. Due to the same reason the excess path length will not be the same for both signals. Therefore, they will not be cancelled out in the ionosphere-free solution. Thus, the terms Δ*sTEC* and Δ*slen* in Eq. (14) and (15) refer to the dual-frequency residual errors due to TEC difference and excess path length, respectively. Their expressions are given by Eqs. (16) and (20). The quantities Δ*TECbend*1 and Δ*TECbend*2 are the differences between TECs along curved and LoS paths and the quantities 1 *len d* and 2 *len d* are the differences between curved and LoS path lengths for *f*1 and *f*2 signals, respectively. The RRE and RREgr are the total residual range errors in the carrier-phase and code combinations, respectively.

The disadvantages of such combinations (Eqs. 14, 15) are that i) the observation noise is increased by a factor depending on frequencies involved in the combination, ii) the ambiguity term of the carrier-phase combination is no more an integer value and iii) only the first order term is eliminated, i.e., higher order terms remain uncorrected. Moreover, this method cannot be applied to single-frequency receivers.

Assuming the same measurement noise on each signal, it can be shown that the carrierphase or code noise will be amplified by a factor of 2.98 for the GPS L1-L2 combination, whereas for the L1-L5 and L2-L3 combinations amplification factors are 2.59 and 16.64, respectively (see Hoque & Jakowski, 2010a). The amplification factor is inversely proportional to the separation of combination frequencies. Since the frequency separation is relatively large for the L1-L5 combination, the amplification factor is the smallest.

Since the first order ionospheric effect on carrier-phase and code pseudoranges (see Eq. 12 and 13) is the same in magnitude but opposite in sign, computing the sum of carrier-phase and code pseudoranges would theoretically eliminate the first order ionospheric term in single frequency measurements. However, the resulting observable would inherit the high code noise and the carrier-phase ambiguity and is therefore, practically not suitable.

#### **2.4.2 Second order ionosphere-free combination**

Receiving signals on three coherent frequencies will allow triple-frequency combinations to eliminate the first and the second order ionospheric terms. The third order ionospheric term and errors due to ray path bending are not fully removed in this approach. Such a triplefrequency combination can be written as (combining code / carrier-phase pseudoranges Eq. (12) / Eq. (13) measured on *f*1, *f*2 and *f*3 frequencies and substituting ρ by Eq. (7), for details see Hoque & jakowski, 2010a).

$$\frac{1}{C} \left[ A \left( \Psi\_1 f\_1^2 - \Psi\_2 f\_2^2 \right) - B \left( \Psi\_1 f\_1^2 - \Psi\_3 f\_3^2 \right) \right] = \rho \underbrace{+ \left( \mathbf{A}\_{\mathrm{TEC}} \right)\_{\mathrm{tr}} + \mathbf{3} \left( \mathbf{A}\_{\mathrm{3}} \right)\_{\mathrm{tr}} + \left( \mathbf{A}\_{\mathrm{lcn}} \right)\_{\mathrm{tr}}}\_{\mathrm{\{RRE}\_{\mathrm{gr}}\}\_{\mathrm{tr}}} \tag{21}$$

$$\frac{1}{C} \left[ A \left( \Phi\_1 f\_1^2 - \Phi\_2 f\_2^2 \right) - B \left( \Phi\_1 f\_1^2 - \Phi\_3 f\_3^2 \right) \right] = \rho \underbrace{- \left( \mathcal{A} \mathfrak{s}\_{\text{TEC}} \right)\_{\text{tr}} - \left( \mathcal{A} \mathfrak{s}\_3 \right)\_{\text{tr}} + \left( \mathcal{A} \mathfrak{s}\_{\text{len}} \right)\_{\text{tr}}}\_{\text{(RRE)}\_{\text{tr}}} \tag{22}$$

In which

386 Global Navigation Satellite Systems – Signal, Theory and Applications

effects is to compute the so called first order ionosphere-free combination of carrier-phase or code pseudoranges measured on two frequencies. However, the second and third order ionospheric terms and errors due to bending of the signal remain uncorrected in this approach. Such a dual-frequency combination can be written in units of length as (combining code / carrier-phase pseudoranges Eq. (12) / Eq. (13) measured on *f*1 and *f*<sup>2</sup>

frequencies and substituting *p* by Eq. (7), for details see Hoque & Jakowski, 2008)

 Ψ ρΔ

22 22 1 2 2 3

22 22 1 2 2 3

 Φ ρΔ

*f f s s ss*

*f f s sss*

− =+ + + − − − 

( ) ( )

*f f f f*

1,2 1,2 ( ) *TEC n ds TEC TEC* = = + ∫ *e LoS b*

( ) <sup>2</sup> 12 1 2 2

> 3 2 2 1 2 3 *u*

*df df <sup>s</sup> f f*

where Ψ1, Ψ2 and Φ1, Φ1 are the measured code and carrier-phase pseudoranges on *f*1 and *f*<sup>2</sup> frequencies, respectively, Δ*s*2 and Δ*s*3 are the dual-frequency second and third order residual terms, respectively. The TEC along a *f*1 path will be different from that along a *f*2 path due to ray path bending. Due to the same reason the excess path length will not be the same for both signals. Therefore, they will not be cancelled out in the ionosphere-free solution. Thus, the terms Δ*sTEC* and Δ*slen* in Eq. (14) and (15) refer to the dual-frequency residual errors due to TEC difference and excess path length, respectively. Their expressions are given by Eqs. (16) and (20). The quantities Δ*TECbend*1 and Δ*TECbend*2 are the differences between TECs along

*ff f f*

*f f*

*len len*

( )

*len d* and 2

and LoS path lengths for *f*1 and *f*2 signals, respectively. The RRE and RREgr are the total

The disadvantages of such combinations (Eqs. 14, 15) are that i) the observation noise is increased by a factor depending on frequencies involved in the combination, ii) the ambiguity term of the carrier-phase combination is no more an integer value and iii) only the first order term is eliminated, i.e., higher order terms remain uncorrected. Moreover, this

2 2 22 11 2 2 1 2

*<sup>q</sup> <sup>s</sup>*

*s*

<sup>−</sup> <sup>=</sup> <sup>−</sup>

Δ

*len*

residual range errors in the carrier-phase and code combinations, respectively.

Δ

<sup>=</sup> <sup>+</sup>

Δ

*K TEC TEC K TEC TEC*

2 2 2 2 1 2 1 2

− =− − − − − − 

2 3 *gr TEC len RRE*

 Δ Δ

*TEC len RRE*

> Δ

(14)

(15)

*end* (17)

(18)

(20)

*len d* are the differences between curved

= (19)

 Δ

 Δ Δ Δ

<sup>−</sup> <sup>−</sup> <sup>=</sup> <sup>=</sup> − − (16)

2 1 2 1

Δ

( ) ( ) *bend bend*

Δ

2 2 1 2

12 12

2 2 1 2

12 12

*ff ff* Φ

*TEC*

*s*

curved and LoS paths and the quantities 1

method cannot be applied to single-frequency receivers.

Δ

*ff ff* Ψ

$$\left(\left(\text{As}\_{\text{TEC}}\right)\_{\text{fr}} = \frac{K}{\mathbb{C}} \Bigg[B \left(\text{ATEC}\_{\text{bend3}} - \text{ATEC}\_{\text{bend1}}\right) - A \left(\text{ATEC}\_{\text{bend2}} - \text{ATEC}\_{\text{bend1}}\right)\right] \tag{23}$$

$$\left(\text{As}\_3\right)\_{tr} = \frac{\mu}{3\text{C}} \frac{\left(f\_2 - f\_3\right)}{f\_2 f\_3} \tag{24}$$

$$\left(\left(\Delta s\_{len}\right)\_{tr} = \frac{1}{C} \left[ B\left(f\_3^2 d\_3^{len} - f\_1^2 d\_1^{len}\right) - A\left(f\_2^2 d\_2^{len} - f\_1^2 d\_1^{len}\right) \right] \tag{25}$$

$$\begin{aligned} A &= \frac{f\_1 f\_2}{f\_1 - f\_2} \\ B &= \frac{f\_1 f\_3}{f\_1 - f\_3} \\ C &= f\_1 (f\_2 - f\_3)(f\_1 + f\_2 + f\_3) \end{aligned} \tag{26}$$

Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches 389

Fig. 1. Examples of global TEC maps and corresponding ionospheric range errors at GPS L1 during night time (0 UT) and day time (14 UT) (http://swaciweb.dlr.de). The dots represent

IPP locations

where *K* = 40.3 m3s-2, Ψ and Φ are the code and carrier-phase pseudoranges, and their subscripts correspond to measured signals on *f*1, *f*2 and *f*3 frequencies, (Δ*s*3)*tr* is the third order residual term and (Δ*sTEC*)*tr* and (Δ*slen*)*tr* are the residual terms due to TEC difference and excess path length, respectively. The quantities Δ*TECbend* and *dI len* are the TEC and path length differences between curved and LoS paths and and their subscripts correspond to received signals on frequencies *f*1, *f*2 and *f*3. The (*RRE*)*tr* and (*RREgr*)tr are the total residual range errors in the triple-frequency carrier-phase and code pseudorange combinations, respectively.

However, as already mentioned, such a multiple frequency combination amplifies all uncorrelated errors or noises (multipath and noise). Assuming the same measurement noise on each signal, it can be shown that the noise will be amplified by a factor of 33.7 in the GPS L1-L2-L5 combination (see Hoque & Jakowski, 2010a).

The Galileo system will transmit signals on four frequencies E2-L1-E1, E5a, E5b and E6 (1575.42, 1176.45, 1207.14 and 1278.75 MHz, respectively). Simultaneous reception of four signals will allow quadruple-frequency combinations to eliminate the first, second and third order ionospheric terms. Such a combination would theoretically eliminate higher order ionospheric terms successfully from the range equation. However, the noise will be amplified by a factor of about 626.13 in the E2L1E1-E5a-E5b-E6 combination (assuming the same measurement noise on each signal) which is about two orders larger than a dualfrequency factor. Therefore, a quadruple-frequency combination is barely pragmatic. However, if the frequency separation is large (e.g., combinations between 4-8 GHz C band and 1-2 GHz L band frequencies), the amplification factor will be small. In such cases, measurements on four frequencies may be useful.

#### **2.4.3 Geometry-free combination**

When microwave signals are transmitted on two frequencies, all the nondispersive effects, e.g., tropospheric delay, satellite and receiver clock offsets, antenna phase centre offsets and variations etc., manipulate the signals on both frequencies in the same way – apart from the ionosphere. Therefore, by differencing code / carrier-phase pseudoranges measured on two frequencies, all non-dispersive terms including *ρ* will be cancelled out giving the estimate of TEC along ray paths as (combining code / carrier-phase pseudoranges Eq. (12) / Eq. (13) measured on *f*1 and *f*2 frequencies and substituting ρ by Eq. (7) and neglecting the second and higher order terms)

$$\text{TECC} = \frac{f\_1^2 f\_2^2}{\mathcal{K} \left(f\_1^2 - f\_2^2\right)} \left[ \left(\Psi\_2 - \Psi\_1\right) + \text{noise}\_{\Psi\_2 - \Psi\_1} \right] \tag{27}$$

$$TEC = \frac{f\_1^2 f\_2^2}{K \left(f\_1^2 - f\_2^2\right)} \left[ \left(\Phi\_1 - \Phi\_2\right) + B\_{\text{ambiquity}} + \text{noise}\_{\Phi\_1 - \Phi\_2} \right] \tag{28}$$

in which *Bambiguity* = *λ*2*N*2 - *λ*1*N*1 is the carrier-phase ambiguity constant where *λ*1, *λ*2 are wave lengths and *N*1, *N*2 are integer ambiguities measured on *f*1 and *f*2 frequencies, 2 1 *noise*Ψ Ψ<sup>−</sup> and 1 2 *noise*Φ Φ<sup>−</sup> are noises (e.g., thermal noise etc.) in code and carrier-phase combinations, respectively. For simplicity different terms such as inter-frequency satellite and receiver biases and multipath effects are not considered.

388 Global Navigation Satellite Systems – Signal, Theory and Applications

where *K* = 40.3 m3s-2, Ψ and Φ are the code and carrier-phase pseudoranges, and their subscripts correspond to measured signals on *f*1, *f*2 and *f*3 frequencies, (Δ*s*3)*tr* is the third order residual term and (Δ*sTEC*)*tr* and (Δ*slen*)*tr* are the residual terms due to TEC difference

length differences between curved and LoS paths and and their subscripts correspond to received signals on frequencies *f*1, *f*2 and *f*3. The (*RRE*)*tr* and (*RREgr*)tr are the total residual range errors in the triple-frequency carrier-phase and code pseudorange combinations,

However, as already mentioned, such a multiple frequency combination amplifies all uncorrelated errors or noises (multipath and noise). Assuming the same measurement noise on each signal, it can be shown that the noise will be amplified by a factor of 33.7 in the GPS

The Galileo system will transmit signals on four frequencies E2-L1-E1, E5a, E5b and E6 (1575.42, 1176.45, 1207.14 and 1278.75 MHz, respectively). Simultaneous reception of four signals will allow quadruple-frequency combinations to eliminate the first, second and third order ionospheric terms. Such a combination would theoretically eliminate higher order ionospheric terms successfully from the range equation. However, the noise will be amplified by a factor of about 626.13 in the E2L1E1-E5a-E5b-E6 combination (assuming the same measurement noise on each signal) which is about two orders larger than a dualfrequency factor. Therefore, a quadruple-frequency combination is barely pragmatic. However, if the frequency separation is large (e.g., combinations between 4-8 GHz C band and 1-2 GHz L band frequencies), the amplification factor will be small. In such cases,

When microwave signals are transmitted on two frequencies, all the nondispersive effects, e.g., tropospheric delay, satellite and receiver clock offsets, antenna phase centre offsets and variations etc., manipulate the signals on both frequencies in the same way – apart from the ionosphere. Therefore, by differencing code / carrier-phase pseudoranges measured on two frequencies, all non-dispersive terms including *ρ* will be cancelled out giving the estimate of TEC along ray paths as (combining code / carrier-phase pseudoranges Eq. (12) / Eq. (13)

ρ

Ψ Ψ

> Φ Φ

( ) ( ) 2 1

( ) ( ) 1 2

<sup>−</sup> are noises (e.g., thermal noise etc.) in code and carrier-phase combinations, respectively. For simplicity different terms such as inter-frequency satellite and receiver

<sup>−</sup> <sup>=</sup> <sup>⎡</sup> −+ + <sup>⎤</sup> <sup>⎣</sup> <sup>⎦</sup> <sup>−</sup>

in which *Bambiguity* = *λ*2*N*2 - *λ*1*N*1 is the carrier-phase ambiguity constant where *λ*1, *λ*2 are wave lengths and *N*1, *N*2 are integer ambiguities measured on *f*1 and *f*2 frequencies, 2 1 *noise*

*ambiguity*

<sup>−</sup> <sup>=</sup> <sup>⎡</sup> <sup>−</sup> <sup>+</sup> <sup>⎤</sup> <sup>⎣</sup> <sup>⎦</sup> <sup>−</sup>

2 1 2 2

Ψ Ψ

*f f TEC noise*

*f f TEC B noise*

Φ Φ

2 2 1 2

1 2

1 2 2 2

*Kf f*

2 2 1 2

1 2

*Kf f*

by Eq. (7) and neglecting the second

(27)

(28)

Ψ Ψ<sup>−</sup> and

*len* are the TEC and path

and excess path length, respectively. The quantities Δ*TECbend* and *dI*

L1-L2-L5 combination (see Hoque & Jakowski, 2010a).

measurements on four frequencies may be useful.

measured on *f*1 and *f*2 frequencies and substituting

biases and multipath effects are not considered.

**2.4.3 Geometry-free combination** 

and higher order terms)

1 2 *noise*Φ Φ

respectively.

Fig. 1. Examples of global TEC maps and corresponding ionospheric range errors at GPS L1 during night time (0 UT) and day time (14 UT) (http://swaciweb.dlr.de). The dots represent IPP locations

Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches 391

The vertical TEC may vary between 1 TECU and 300 TECU depending on a number of factors such as local time, geographic location, season, solar activity level etc. The frequency dependence of the first order ionospheric term has been plotted for elevations 5° and 30° in Fig. 2 at different levels of ionospheric ionization characterized by vertical TECs such as i) 250 and 150 TECU correspond to TEC during extreme space weather conditions, ii) 50 TECU corresponds to mid latitude day time and iii) 5 TECU corresponds to mid latitude night time

Fig. 2. Frequency dependence of the first order term at different levels of ionospheric

L-band frequencies (1 – 2 GHz) during times of high TEC at low elevation angles.

As Fig. 2 demonstrates, the first order ionospheric term can be more than 100 m at GNSS

Higher order ionospheric terms include the second and third order ionospheric terms, and the excess path length. Equations (5, 6, 8) indicate that the second order term depends on the electron density *ne* and as well as on the geomagnetic induction *B*. The electron gyro frequency *fg* = *eB*/(2*πm*) is usually less than 1.4 MHz. The value of *B* can be derived as ~5x10-5 Tesla for *fg* = 1.4 MHz and considered constant throughout the propagation. Thus, for the worst case condition with *fg* = 1.4 MHz and Θ = 0, the second order term can be simplified as

( ) <sup>7</sup>

3 11.28 10 *Igr d TEC f*

approximation, the frequency dependence of the second order term at different levels of

Figure 3 shows that during the worst case conditions the second order ionospheric term on code observables can be as big as about 500 millimeters at GNSS L-band frequencies. Due to

*Igr d* is measured in meters, TEC in electrons/m2 and *f* in Hz. Using the above

<sup>×</sup> <sup>=</sup> (30)

2

ionospheric ionization and elevation angles has been plotted in Fig. 3.

TEC.

ionization and elevation angles.

**3.1.1 Second order term** 

(using Eqs. 5 and 8)

where ( ) <sup>2</sup>

In case of cycle slips (a jump in carrier-phase ambiguity constant due to loss of signal tracking results in discontinuous arcs of phase data) in the phase data, the wide-lane combination method of Blewitt (1987) can be applied for the correction. While the TEC estimated by the carrier-phase difference Eq. (28) is precise and smooth but biased by an unknown phase ambiguity constant, the TEC estimated by the code pseudorange difference Eq. (27) is noisy and less precise but not ambiguous. In order to obtain an absolute and precise estimate of TEC, the accurate phase measurements needs to be levelled to the calibrated absolute code measurements by a least square method.

To derive an elevation independent vertical TEC from a slant TEC measurement, the ionosphere is assumed to be composed of a single thin layer at a representative height of about 350, 400 or 450 km from the earth's surface. The intersection point between a slant ray path and the thin layer is called an ionospheric piercing point (IPP). A mapping function is used to convert the slant STEC to vertical VTEC at the IPP or vice versa (details of the derivation is given in Hoque & Jakowski, 2008).

$$\text{STEC} \mid \text{VTEC} \approx \frac{\left(h\_m + R\_E\right)}{\sqrt{\left(h\_m + R\_E\right)^2 - \left(R\_h + R\_E\right)^2 \cos^2 \beta}}\tag{29}$$

where *hm* is the height of maximum electron density (varies between 250 -450 km), *RE* is the earth's mean radius (~ 6371 km), *Rh* is the receiver height above the earth's surface and *β* is the elevation angle.

Based on similar techniques using observation from more than hundred worldwide GNSS ground stations, German Aerospace Center (DLR) Neustrelitz computes vertical TEC estimates at numerous IPPs worldwide. Thus, TEC maps are produced by assigning IPP measurements to homogeneous latitude and longitude grid points as shown in Fig. 1. European and global TEC maps and 1-hour-ahead forecasts are distributed via the operational space-weather and ionosphere data service SWACI (Space Weather Application Center Ionosphere, http://swaciweb.dlr.de, see also Jakowski et al., 2011) to the international community with an update rate of 5 minutes. The advantage of such services is that single frequency GNSS users can correct the ionospheric propagation effect in near real time.

## **3. Estimation of ionospheric effects**

#### **3.1 First- and higher-order ionospheric terms**

Equations (12) and (13) indicate that the signal delay caused by the first order term is equal in magnitude but opposite in sign on GNSS carrier-phase and code pseudoranges, i.e., the carrier-phase pseudorange is advanced while code pseudorange is retarded. The first order term is directly proportional to the TEC encountered by the satellite signal during its travel through the ionosphere and inversely proportional to the square of the signal frequency. The first order term includes about 99% of the total ionospheric effect. Therefore, if the frequency and link related slant TEC are known, the first order propagation effect can easily be computed and corrected. If the TEC map is available, the slant TEC can be computed simply multiplying the vertical TEC at the IPP by the mapping function (e.g., Eq. 29).

The vertical TEC may vary between 1 TECU and 300 TECU depending on a number of factors such as local time, geographic location, season, solar activity level etc. The frequency dependence of the first order ionospheric term has been plotted for elevations 5° and 30° in Fig. 2 at different levels of ionospheric ionization characterized by vertical TECs such as i) 250 and 150 TECU correspond to TEC during extreme space weather conditions, ii) 50 TECU corresponds to mid latitude day time and iii) 5 TECU corresponds to mid latitude night time TEC.

Fig. 2. Frequency dependence of the first order term at different levels of ionospheric ionization and elevation angles.

As Fig. 2 demonstrates, the first order ionospheric term can be more than 100 m at GNSS L-band frequencies (1 – 2 GHz) during times of high TEC at low elevation angles.

#### **3.1.1 Second order term**

390 Global Navigation Satellite Systems – Signal, Theory and Applications

In case of cycle slips (a jump in carrier-phase ambiguity constant due to loss of signal tracking results in discontinuous arcs of phase data) in the phase data, the wide-lane combination method of Blewitt (1987) can be applied for the correction. While the TEC estimated by the carrier-phase difference Eq. (28) is precise and smooth but biased by an unknown phase ambiguity constant, the TEC estimated by the code pseudorange difference Eq. (27) is noisy and less precise but not ambiguous. In order to obtain an absolute and precise estimate of TEC, the accurate phase measurements needs to be levelled to the

To derive an elevation independent vertical TEC from a slant TEC measurement, the ionosphere is assumed to be composed of a single thin layer at a representative height of about 350, 400 or 450 km from the earth's surface. The intersection point between a slant ray path and the thin layer is called an ionospheric piercing point (IPP). A mapping function is used to convert the slant STEC to vertical VTEC at the IPP or vice versa (details of the

<sup>+</sup> <sup>≈</sup>

where *hm* is the height of maximum electron density (varies between 250 -450 km), *RE* is the earth's mean radius (~ 6371 km), *Rh* is the receiver height above the earth's surface and *β* is

Based on similar techniques using observation from more than hundred worldwide GNSS ground stations, German Aerospace Center (DLR) Neustrelitz computes vertical TEC estimates at numerous IPPs worldwide. Thus, TEC maps are produced by assigning IPP measurements to homogeneous latitude and longitude grid points as shown in Fig. 1. European and global TEC maps and 1-hour-ahead forecasts are distributed via the operational space-weather and ionosphere data service SWACI (Space Weather Application Center Ionosphere, http://swaciweb.dlr.de, see also Jakowski et al., 2011) to the international community with an update rate of 5 minutes. The advantage of such services is that single frequency GNSS users can correct the ionospheric propagation effect in near real

Equations (12) and (13) indicate that the signal delay caused by the first order term is equal in magnitude but opposite in sign on GNSS carrier-phase and code pseudoranges, i.e., the carrier-phase pseudorange is advanced while code pseudorange is retarded. The first order term is directly proportional to the TEC encountered by the satellite signal during its travel through the ionosphere and inversely proportional to the square of the signal frequency. The first order term includes about 99% of the total ionospheric effect. Therefore, if the frequency and link related slant TEC are known, the first order propagation effect can easily be computed and corrected. If the TEC map is available, the slant TEC can be computed simply multiplying the vertical TEC at the IPP by the mapping function (e.g., Eq.

( ) ( )( ) 2 2 <sup>2</sup>

*h R*

*mE hE*

*h R RR*

+ −+

*m E*

cos

β

(29)

calibrated absolute code measurements by a least square method.

derivation is given in Hoque & Jakowski, 2008).

**3. Estimation of ionospheric effects** 

**3.1 First- and higher-order ionospheric terms** 

the elevation angle.

time.

29).

/

*STEC VTEC*

Higher order ionospheric terms include the second and third order ionospheric terms, and the excess path length. Equations (5, 6, 8) indicate that the second order term depends on the electron density *ne* and as well as on the geomagnetic induction *B*. The electron gyro frequency *fg* = *eB*/(2*πm*) is usually less than 1.4 MHz. The value of *B* can be derived as ~5x10-5 Tesla for *fg* = 1.4 MHz and considered constant throughout the propagation. Thus, for the worst case condition with *fg* = 1.4 MHz and Θ = 0, the second order term can be simplified as (using Eqs. 5 and 8)

$$d\_{lg}^{(2)} = \frac{11.28 \times 10^7}{f^3} \text{TEC} \tag{30}$$

where ( ) <sup>2</sup> *Igr d* is measured in meters, TEC in electrons/m2 and *f* in Hz. Using the above approximation, the frequency dependence of the second order term at different levels of ionospheric ionization and elevation angles has been plotted in Fig. 3.

Figure 3 shows that during the worst case conditions the second order ionospheric term on code observables can be as big as about 500 millimeters at GNSS L-band frequencies. Due to

Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches 393

For a GNSS user in the northern hemisphere the magnitude of the second order term is the largest when the signal is received from a satellite in southward direction. However, for a user in the southern hemisphere the scenario is reversed, i.e., the largest effect is observed when the signal is received from a satellite in northward direction. Figure 4 shows that the magnitude of the second order term and its sign differ depending on the user location on the earth and direction of the signal reception. Therefore, such non systematic effects cannot be cancelled out by averaging GNSS measurements over long period at a certain user location.

0.6577*NmTEC* (obtained by analytical integration of the Chapman layer, Hoque & Jakowski, 2008, see also Brunner & Gu, 1991; Hartmann & Leitinger, 1984; Leitinger & Putz, 1988) where *Nm* is the maximum ionospheric ionization. Therefore, assuming the worst case condition with *fg* = 1.4 MHz and Θ = 0, the third order term can be simplified as (using Eqs.

> ( ) ( ) <sup>3</sup> <sup>14</sup> <sup>4</sup> *Igr* 1602.81 2.37 10 *<sup>m</sup> TEC d N*

where VTEC is the TEC in vertical direction and *H* is the atmospheric scale height. The parameter *H* can be assumed as 70 km for a rough estimation of the third order ionospheric term. Using the above approximations (Eqs. 31 and 32) the frequency dependence of the third order term at different levels of ionospheric ionization and elevation angles has been

*Igr d* is measured in meters, TEC is the slant TEC and measured in electrons/m2, *f* in Hz and *Nm* is measured in m-3. If the vertical TEC is known, *Nm* can be computed assuming a Chapman profile for the ionosphere by the following expression (Hoque & Jakowski, 2007).

∫*n ds <sup>e</sup>* (see Eq. 9) which can be simplified as

*f* = +× (31)

*VTEC HN* = 4.13 *<sup>m</sup>* (32)

*Igr d* at different levels of ionospheric

**3.1.2 Third order term** 

5 and 9)

where ( ) <sup>3</sup>

plotted in Fig. 5.

The third order term depends on the integral <sup>2</sup>

Fig. 5. Frequency dependence of the third order term ( ) <sup>3</sup>

ionization and elevation angles

the dependency on *B* field, the second order term depends on the receiver's geographic / geomagnetic position and direction of the signal reception. Such dependencies are given in Fig. 4. The simulation has been made using a two dimensional ray tracing program (Hoque & Jakowski, 2008) which includes International Geomagnetic Reference Field (IGRF) (Mandea & Macmillan, 2000) model for magnetic field computation along ray paths. A single layered Chapman profile (Eq. 40, see Rishbeth & Garriott, 1969) with a maximum ionization of 4.96×1012 m-3 at 350 km altitude and atmospheric scale height of 70 km has been used, and the corresponding vertical TEC is found 143 TECU.

Fig. 3. Frequency dependence of the second order term ( ) <sup>2</sup> *Igr d* at different levels of ionospheric ionization and elevation angles

Fig. 4. Azimuth dependency of the second order term at GPS L1 frequency for elevations 5°, 30° and 60° and VTEC = 143 TECU. The receiver position is considered at geographic 25° N, 10° E and at its geomagnetic conjugate position 3.4° S, 11° E. The symbols N, E, S and W correspond to the geographic north, east, south and west directions, respectively

For a GNSS user in the northern hemisphere the magnitude of the second order term is the largest when the signal is received from a satellite in southward direction. However, for a user in the southern hemisphere the scenario is reversed, i.e., the largest effect is observed when the signal is received from a satellite in northward direction. Figure 4 shows that the magnitude of the second order term and its sign differ depending on the user location on the earth and direction of the signal reception. Therefore, such non systematic effects cannot be cancelled out by averaging GNSS measurements over long period at a certain user location.

#### **3.1.2 Third order term**

392 Global Navigation Satellite Systems – Signal, Theory and Applications

the dependency on *B* field, the second order term depends on the receiver's geographic / geomagnetic position and direction of the signal reception. Such dependencies are given in Fig. 4. The simulation has been made using a two dimensional ray tracing program (Hoque & Jakowski, 2008) which includes International Geomagnetic Reference Field (IGRF) (Mandea & Macmillan, 2000) model for magnetic field computation along ray paths. A single layered Chapman profile (Eq. 40, see Rishbeth & Garriott, 1969) with a maximum ionization of 4.96×1012 m-3 at 350 km altitude and atmospheric scale height of 70 km has

Fig. 4. Azimuth dependency of the second order term at GPS L1 frequency for elevations 5°, 30° and 60° and VTEC = 143 TECU. The receiver position is considered at geographic 25° N, 10° E and at its geomagnetic conjugate position 3.4° S, 11° E. The symbols N, E, S and W correspond to the geographic north, east, south and west directions, respectively

*Igr d* at different levels of ionospheric

been used, and the corresponding vertical TEC is found 143 TECU.

Fig. 3. Frequency dependence of the second order term ( ) <sup>2</sup>

ionization and elevation angles

The third order term depends on the integral <sup>2</sup> ∫*n ds <sup>e</sup>* (see Eq. 9) which can be simplified as 0.6577*NmTEC* (obtained by analytical integration of the Chapman layer, Hoque & Jakowski, 2008, see also Brunner & Gu, 1991; Hartmann & Leitinger, 1984; Leitinger & Putz, 1988) where *Nm* is the maximum ionospheric ionization. Therefore, assuming the worst case condition with *fg* = 1.4 MHz and Θ = 0, the third order term can be simplified as (using Eqs. 5 and 9)

$$d\_{gr}^{(3)} = \left(1602.81N\_m + 2.37 \times 10^{14}\right) \frac{TEC}{f^4} \tag{31}$$

where ( ) <sup>3</sup> *Igr d* is measured in meters, TEC is the slant TEC and measured in electrons/m2, *f* in Hz and *Nm* is measured in m-3. If the vertical TEC is known, *Nm* can be computed assuming a Chapman profile for the ionosphere by the following expression (Hoque & Jakowski, 2007).

$$VTEC = 4.13HN\_m \tag{32}$$

where VTEC is the TEC in vertical direction and *H* is the atmospheric scale height. The parameter *H* can be assumed as 70 km for a rough estimation of the third order ionospheric term. Using the above approximations (Eqs. 31 and 32) the frequency dependence of the third order term at different levels of ionospheric ionization and elevation angles has been plotted in Fig. 5.

Fig. 5. Frequency dependence of the third order term ( ) <sup>3</sup> *Igr d* at different levels of ionospheric ionization and elevation angles

Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches 395

We see that at low elevation angle (< 15°), Δ*sTEC* is the largest and it decreases very rapidly with increasing the elevation angle. The second order term Δ*s*2 is determined for an azimuth angle 180° at a receiver position at geographic 50° N and 15° E. Although Δ*s*2 is less than the Δ*sTEC* at low elevation angles, it exceeds Δ*sTEC* at higher elevation angles (> 25°). The Δ*s*<sup>2</sup> does not reduce significantly with increasing the elevation angle and therefore, it cannot be ignored even at zenith. The excess path length Δ*slen* decreases with increasing the elevation angle very rapidly and vanishes at zenith. The third order term Δ*s*3 is small (< 5 mm) but it can be bigger than Δ*sTEC* and Δ*slen* at very high (> 60°) elevation angles. We find that the magnitude of the RREgr is much higher than the RRE. This is mainly due to different signs of Δ*sTEC* and Δ*slen* in the RRE and RREgr expressions. In case of the code combination Eq. (14), Δ*slen* is additive to other terms whereas it is subtractive in the carrier-phase combination Eq. (15). Additionally, the Δ*s*2 and Δ*s*3 are two and three times higher in the code combination

Fig. 7. Residual terms in the dual-frequency GPS L1-L2 combination for an ionospheric

For the same ionospheric ionization, the residual terms in the triple-frequency GPS L1-L2-L5 combination (Eqs. 21-22) are plotted in Fig. 8. It shows that the magnitude of (RREgr)tr is much higher than the magnitude of (RRE)tr. The reason is already discussed for the dual-

We find that the GPS L1-L2 residual terms Δ*s*3, Δ*sTEC* and Δ*slen* are about 2.4 times higher than the GPS L1-L2-L5 residual terms (Δ*s*3)*tr*, (Δ*sTEC*)*tr* and Δ*slen*. The sum of all residual terms, i.e., *RRE* and *RREgr* are found to be more than three times higher for the L1-L2

Comparing the dual- and triple-frequency carrier-phase combinations Eq. (15) and Eq. (22), we see that the signs of Δ*sTEC* and Δ*s*3 are positive in the dual-frequency combination whereas their signs are negative in the triple-frequency combination. However, the sign of

compared to the phase combination.

ionization of VTEC = 250 TEC units

combination than the L1-L2-L5 combination.

frequency case.

Figure 5 shows that the third order term on code observables can be as big as 50 mm at low elevation angles during times of high TEC.

#### **3.2 Estimate of LoS propagation assumption error**

Due to the ray path bending satellite signals propagate through curvature paths instead of straight line of sight paths. However, a curvature path length and the corresponding LoS path length are not equal rather the curvature path is slightly longer than the LoS one. The difference between them is defined as the excess path length and it can be computed by the following formula given by Jakowski et al. (1994).

$$d\_{I}^{len} = \frac{b\_1}{f^4} \left( \frac{1}{\left(1 - b\_2 \cos^2 \beta\right)^{1/2}} - 1 \right) \text{TEC}^2 \tag{33}$$

where *b*1 = 2.495×108, *b*2 = 0.8592 and *β* is the elevation angle. The excess path length *dI len* will be estimated in millimeters when *β* is measured in radians, *f* is in MHz and TEC is in TEC units. The frequency dependence of the excess path length has been plotted in Fig. 6.

Fig. 6. Frequency dependence of the excess path length at different levels of ionospheric ionization and elevation angles

Figure 6 shows that at the L2 frequency, the excess path length can be as big as 100 mm at low elevation angles during times of high TEC such as VTEC = 250 TECU.

#### **3.3 Estimates of residual terms in the ionosphere-free solution**

Although residual terms in ionosphere-free solutions are less than 1% of the first order ionospheric effect, they cannot be ignored if centimeter / millimeter level accuracy is required in GNSS positioning and timing applications. A plot showing comparison of dualfrequency GPS L1-L2 residual terms is given in Fig. 7 for better understanding of their relative influences on precise range estimation. For this, the ray tracing tool has been used in which the ionosphere is modelled by a Chapman layer with a peak density of 7.75×1012 m-3 at 350 km altitude and scale height of 78 km, and corresponding VTEC = 250 TECU.

394 Global Navigation Satellite Systems – Signal, Theory and Applications

Figure 5 shows that the third order term on code observables can be as big as 50 mm at low

Due to the ray path bending satellite signals propagate through curvature paths instead of straight line of sight paths. However, a curvature path length and the corresponding LoS path length are not equal rather the curvature path is slightly longer than the LoS one. The difference between them is defined as the excess path length and it can be computed by the

( )

1 cos

where *b*1 = 2.495×108, *b*2 = 0.8592 and *β* is the elevation angle. The excess path length *dI*

units. The frequency dependence of the excess path length has been plotted in Fig. 6.

Fig. 6. Frequency dependence of the excess path length at different levels of ionospheric

low elevation angles during times of high TEC such as VTEC = 250 TECU.

**3.3 Estimates of residual terms in the ionosphere-free solution** 

Figure 6 shows that at the L2 frequency, the excess path length can be as big as 100 mm at

Although residual terms in ionosphere-free solutions are less than 1% of the first order ionospheric effect, they cannot be ignored if centimeter / millimeter level accuracy is required in GNSS positioning and timing applications. A plot showing comparison of dualfrequency GPS L1-L2 residual terms is given in Fig. 7 for better understanding of their relative influences on precise range estimation. For this, the ray tracing tool has been used in which the ionosphere is modelled by a Chapman layer with a peak density of 7.75×1012 m-3

at 350 km altitude and scale height of 78 km, and corresponding VTEC = 250 TECU.

*f b*

4 1/2 <sup>2</sup> 2

*<sup>b</sup> <sup>d</sup> TEC*

1 2

⎛ ⎞ ⎜ ⎟ <sup>=</sup> <sup>−</sup>

<sup>−</sup> ⎝ ⎠

be estimated in millimeters when *β* is measured in radians, *f* is in MHz and TEC is in TEC

β

<sup>1</sup> <sup>1</sup>

(33)

*len* will

elevation angles during times of high TEC.

**3.2 Estimate of LoS propagation assumption error** 

following formula given by Jakowski et al. (1994).

ionization and elevation angles

*len I*

We see that at low elevation angle (< 15°), Δ*sTEC* is the largest and it decreases very rapidly with increasing the elevation angle. The second order term Δ*s*2 is determined for an azimuth angle 180° at a receiver position at geographic 50° N and 15° E. Although Δ*s*2 is less than the Δ*sTEC* at low elevation angles, it exceeds Δ*sTEC* at higher elevation angles (> 25°). The Δ*s*<sup>2</sup> does not reduce significantly with increasing the elevation angle and therefore, it cannot be ignored even at zenith. The excess path length Δ*slen* decreases with increasing the elevation angle very rapidly and vanishes at zenith. The third order term Δ*s*3 is small (< 5 mm) but it can be bigger than Δ*sTEC* and Δ*slen* at very high (> 60°) elevation angles. We find that the magnitude of the RREgr is much higher than the RRE. This is mainly due to different signs of Δ*sTEC* and Δ*slen* in the RRE and RREgr expressions. In case of the code combination Eq. (14), Δ*slen* is additive to other terms whereas it is subtractive in the carrier-phase combination Eq. (15). Additionally, the Δ*s*2 and Δ*s*3 are two and three times higher in the code combination compared to the phase combination.

Fig. 7. Residual terms in the dual-frequency GPS L1-L2 combination for an ionospheric ionization of VTEC = 250 TEC units

For the same ionospheric ionization, the residual terms in the triple-frequency GPS L1-L2-L5 combination (Eqs. 21-22) are plotted in Fig. 8. It shows that the magnitude of (RREgr)tr is much higher than the magnitude of (RRE)tr. The reason is already discussed for the dualfrequency case.

We find that the GPS L1-L2 residual terms Δ*s*3, Δ*sTEC* and Δ*slen* are about 2.4 times higher than the GPS L1-L2-L5 residual terms (Δ*s*3)*tr*, (Δ*sTEC*)*tr* and Δ*slen*. The sum of all residual terms, i.e., *RRE* and *RREgr* are found to be more than three times higher for the L1-L2 combination than the L1-L2-L5 combination.

Comparing the dual- and triple-frequency carrier-phase combinations Eq. (15) and Eq. (22), we see that the signs of Δ*sTEC* and Δ*s*3 are positive in the dual-frequency combination whereas their signs are negative in the triple-frequency combination. However, the sign of

Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches 397

average peak height. Later Kedar et al. (2003) assumed the ionosphere as a single layer at a 400 km altitude for estimating the effect of the second order GPS ionospheric correction on receiver positions. Hernandez-Pajares et al. (2007) considered the ionosphere as a single layer at a 450 km altitude to estimate the impact of the second order ionospheric term on

profile shape which is not available to the GNSS users; they only have TEC information

at the IPP are very suitable for practical use. However, such assumptions lead up to 2 mm

for the magnetic field component and consider it constant throughout the propagation.

along any receiver-to-satellite link geometries inside European geographic latitude 30 – 65°

12

cos

222

*i i i*

( ) ( ) ( )

,,, , , , ,

 Θ

αα

As an alternative approach, we (Hoque & Jakowski, 2007) assume an average value *B*cos

errors in the second order ionospheric term computation (Hoque & Jakowski, 2008).

12 1 2 1.1283 10

1 11 2 *B y ry r* cos cos

ς βφλ

= = =

ς βφ

ς βφ

1 2 1

Inside the ray tracing program, the IGRF model has been used to compute cos *B*

 =− + − − α

> *r a r b y c*

The parameters *r*1, *r*2 and *y*1 are the functions of the receiver-to-satellite elevation angle *β*, geographic latitude *φ* and longitude *λ* at the receiver position. The quantity α is the receiverto-satellite azimuth angle and α' is the modified azimuth angle. The quantities *ai*, *bi* and *ci* are the polynomial coefficients. Thirty polynomial coefficients have been derived for the European region (30° - 65° N, 15° W- 45° E) by least squares fitting of ray tracing results.

paths. For details and values of the polynomial coefficients we refer to Hoque & Jakowski (2007). Using such a correction formula and knowing the TEC value, the second order term can be corrected to the 2-3 millimeter accuracy level for a vertical TEC level of 100 TEC units. The formula can be adapted for other geographic regions too after deriving new set of

It has been found that the second term of Eq. (9) is less than the first term by about 1-2

orders of magnitude. As already discussed in the section 3.1.2, the integral <sup>2</sup>

*ff f f*

( ) *s BTEC*

<sup>×</sup> = ⋅⋅ <sup>+</sup>

along ray paths. Therefore, assumptions of a thin ionospheric layer and cos *B*

Based on simulation studies we derived a correction formula for the *B*cos

2

Δ

Θ

along ray paths requires the knowledge of the ionospheric

Θ

Θ

(35)

sin 2 cos ′ (36)

(37)

Θ

along ray

∫*n ds <sup>e</sup>* can be

computation

computation

Θ

geodetic estimates.

The computation of cos *B*

N and longitude 15° W – 45° E.

polynomial coefficients.

**4.2 Third order term correction** 

In which

Θ

Δ*slen* is negative in the dual-frequency combination and positive in the triple-frequency combination. Again, the magnitude of Δ*sTEC* is higher than the magnitude of Δ*slen* for both combinations. As a result, the triple-frequency (RRE)tr is found to be negative, i.e., the corrected *ρ* is longer than the uncorrected one whereas the dual-frequency RRE is found to be positive, i.e., the corrected *ρ* is shorter than the uncorrected *ρ*. Similarly, it can be shown that (RREgr)tr is positive in the triple-frequency combination and RREgr is negative in the dual-frequency combination. These relations are true for combination frequencies *f*1 > *f*2 > *f*3.

Fig. 8. Residual terms in the triple-frequency GPS L1-L2-L5 combination for an ionospheric ionization of VTEC = 250 TEC units

### **4. Correction of higher order ionospheric terms**

#### **4.1 Second order term correction**

The estimation of the second order term requires computation of the geomagnetic induction and its direction with respect to the propagation direction along ray paths. Since this computation is very cumbersome, a common practice is to assume the ionosphere as a single thin layer at a certain altitude and compute *B*cosΘ at the IPP and consider it constant throughout the propagation. Thus, the second order term coefficient *q* (Eq. 8) can be written as

$$q \approx 2.2566 \times 10^{12} \,\text{ ${}^{\circ}$ } \,\text{cos}\,\Theta^{\circ} \text{TEC} \tag{34}$$

where TEC is the total electron content along ray paths, *B\** is the magnitude of *B*, and Θ\* is the angle between the magnetic field vector and the wave direction at the IPP (the symbol \* denotes the values at the single layer).

Bassiri and Hajj (1993) were the first to propose such a single thin layered ionosphere for the second order ionospheric correction; they choose the 300 km as a representative global average peak height. Later Kedar et al. (2003) assumed the ionosphere as a single layer at a 400 km altitude for estimating the effect of the second order GPS ionospheric correction on receiver positions. Hernandez-Pajares et al. (2007) considered the ionosphere as a single

layer at a 450 km altitude to estimate the impact of the second order ionospheric term on geodetic estimates.

The computation of cos *B* Θ along ray paths requires the knowledge of the ionospheric profile shape which is not available to the GNSS users; they only have TEC information along ray paths. Therefore, assumptions of a thin ionospheric layer and cos *B* Θ computation at the IPP are very suitable for practical use. However, such assumptions lead up to 2 mm errors in the second order ionospheric term computation (Hoque & Jakowski, 2008).

As an alternative approach, we (Hoque & Jakowski, 2007) assume an average value *B*cosΘ for the magnetic field component and consider it constant throughout the propagation. Based on simulation studies we derived a correction formula for the *B*cosΘ computation along any receiver-to-satellite link geometries inside European geographic latitude 30 – 65° N and longitude 15° W – 45° E.

$$\text{As}\_2 = \frac{1.1283 \times 10^{12}}{f\_1 f\_2 (f\_1 + f\_2)} \cdot \overline{B \cos \Theta} \cdot \text{TEC} \tag{35}$$

In which

396 Global Navigation Satellite Systems – Signal, Theory and Applications

Δ*slen* is negative in the dual-frequency combination and positive in the triple-frequency combination. Again, the magnitude of Δ*sTEC* is higher than the magnitude of Δ*slen* for both combinations. As a result, the triple-frequency (RRE)tr is found to be negative, i.e., the corrected *ρ* is longer than the uncorrected one whereas the dual-frequency RRE is found to be positive, i.e., the corrected *ρ* is shorter than the uncorrected *ρ*. Similarly, it can be shown that (RREgr)tr is positive in the triple-frequency combination and RREgr is negative in the dual-frequency combination. These relations are true for combination frequencies

Fig. 8. Residual terms in the triple-frequency GPS L1-L2-L5 combination for an ionospheric

The estimation of the second order term requires computation of the geomagnetic induction and its direction with respect to the propagation direction along ray paths. Since this computation is very cumbersome, a common practice is to assume the ionosphere as a single thin layer at a certain altitude and compute *B*cosΘ at the IPP and consider it constant throughout the propagation. Thus, the second order term coefficient *q* (Eq. 8) can be written

12 \* \* *q* ≈ × 2.2566 10 cos *B TEC*

where TEC is the total electron content along ray paths, *B\** is the magnitude of *B*, and Θ\* is the angle between the magnetic field vector and the wave direction at the IPP (the symbol \*

Bassiri and Hajj (1993) were the first to propose such a single thin layered ionosphere for the second order ionospheric correction; they choose the 300 km as a representative global

Θ

(34)

*f*1 > *f*2 > *f*3.

ionization of VTEC = 250 TEC units

**4.1 Second order term correction** 

denotes the values at the single layer).

as

**4. Correction of higher order ionospheric terms** 

$$\overline{B\cos\Theta} = -y\_1 \cos a + \left| \sqrt{r\_1^2 - y\_1^2 \sin^2 a} \right| - 2r\_2 \cos a' \tag{36}$$

$$\begin{aligned} r\_1 &= \mathfrak{z}\left(\mathcal{J}, \phi, \mathcal{A}, a\_i\right) \\ r\_2 &= \mathfrak{z}\left(\mathcal{J}, \phi, b\_i\right) \\ y\_1 &= \mathfrak{z}\left(\mathcal{J}, \phi, c\_i\right) \end{aligned} \tag{37}$$

The parameters *r*1, *r*2 and *y*1 are the functions of the receiver-to-satellite elevation angle *β*, geographic latitude *φ* and longitude *λ* at the receiver position. The quantity α is the receiverto-satellite azimuth angle and α' is the modified azimuth angle. The quantities *ai*, *bi* and *ci* are the polynomial coefficients. Thirty polynomial coefficients have been derived for the European region (30° - 65° N, 15° W- 45° E) by least squares fitting of ray tracing results. Inside the ray tracing program, the IGRF model has been used to compute cos *B* Θ along ray paths. For details and values of the polynomial coefficients we refer to Hoque & Jakowski (2007). Using such a correction formula and knowing the TEC value, the second order term can be corrected to the 2-3 millimeter accuracy level for a vertical TEC level of 100 TEC units. The formula can be adapted for other geographic regions too after deriving new set of polynomial coefficients.

#### **4.2 Third order term correction**

It has been found that the second term of Eq. (9) is less than the first term by about 1-2 orders of magnitude. As already discussed in the section 3.1.2, the integral <sup>2</sup> ∫*n ds <sup>e</sup>* can be

Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches 399

*d2TEC*/*dβ2* have been plotted as functions of elevation angle in Fig. 9. The *dTEC*/*dβ* has been calculated dividing the TEC difference between two measurement epochs by the corresponding elevation angle difference. Then, *d2TEC*/*dβ2* has been calculated dividing the *dTEC*/*dβ* difference between two measurement epochs by the corresponding elevation angle

Comparing plots in Fig. 9, we see that although the dependency of *dlen* on the *dTEC*/*dβ* is not straight forward, its dependency on the *d2TEC*/*dβ2* is obvious at low elevation angles (< 20°). Thus, the magnitude of the *dlen* depends on the magnitude of TEC as well as on the magnitude of *d2TEC*/*dβ2*. Considering this, functional dependencies have been studied separately for different parameters to develop correction formulas. For this, ray tracing calculation has been carried out for different geometrical and ionospheric conditions varying elevation and Chapman layer parameters *H*, *Nm* and *hm*. Thus, the following

*len* correction.

ray tracing Eq. (33)

H = 60 km

4

3

2

Excess Path /cm

1

0

1.5

 2

1

0.5

 *TEC/dβ* /(TECU/deg )

*2*

0

<sup>1</sup> 1 cos

500

H = 60 km H = 80 km

H = 80 km

ray tracing Eq. (33)

Elevation /deg 0 30 60 90

*len*, TEC (see right scale), *dTEC*/*dβ* and *d2TEC*/*dβ2* for

<sup>2</sup> <sup>2</sup>

β β (41)

H = 60 km H = 80 km

( )

2

1 cos

1 2

β

*<sup>a</sup> d TEC d TEC <sup>a</sup> f d a*

2

4 1/2 2 3

<sup>−</sup> ⎝ ⎠ ⎝ ⎠

⎛ ⎞ ⎛ ⎞ ⎜ ⎟ <sup>=</sup> <sup>−</sup> <sup>+</sup> ⎜ ⎟

375

250

125

0

Elevation /deg Elevation /deg 0 30 60 90 0 30 60 90

*d*-0.5

*2*

TEC/TECU

difference.

3

2

1

Excess Path /cm

0

20

15

10

*dTEC/dβ*/(TECU/deg)

5

0

formula has been obtained for the *dI*

Elevation /deg 0 30 60 90

Fig. 9. Elevation angle dependence of *dI*

*len I*

the Chapman layer parameter *H* = 60 and 80 km

simplified as 0.6577*NmTEC*. Thus, the third order residual term can be approximated by the first term only as (using Eqs. 9 and 19)

$$\Delta \mathbf{s}\_3 = \frac{534.27}{f\_1^2 f\_2^2} N\_m T \text{EC} \tag{38}$$

The third order term Δ*s*3 will be measured in meters when *f* is measured in Hz and the maximum ionization *Nm* and TEC in m-3 and electrons/m2, respectively.

#### **4.3 New approaches for correcting LoS propagation assumption errors**

#### **4.3.1 Excess path length correction**

As we have seen in the section 3.2, the excess path length *dI len* can be computed by Eq. (33). There is another formula published by Hoque & Jakowski (2008) for the excess path length computation.

$$d\_{l}^{len} = \frac{7.5 \times 10^{-5} \exp\left(-2.13 \beta\right) \text{TEC}^2}{f^4 H \left(l m\right)^{1/8}}\tag{39}$$

where *dI len* is measured in meters, TEC is in TEC units, frequency *f* in GHz, atmospheric scale height *H* and maximum ionization height *hm* in kilometers and elevation *β* in radians. Comparing both formulas we see that Eq. (33) requires only TEC and elevation information as inputs whereas Eq. (39) additionally requires ionospheric parameters *H* and *hm*. However, these parameters are not easy to estimate in practical cases.

Both the correction formulas are derived based on simulation studies using Chapman profiles for the ionosphere. The Chapman profile (Rishbeth & Garriott, 1969) has been proved very useful for modeling ionospheric correction. It describes the electron density distribution *ne* as a function of height *h* in the ionosphere as

$$m\_{\varepsilon}(h) = N m \exp(0.5(1 - z - \exp(-z)))\tag{40}$$

where *Nm* is the maximum ionization and *z* = (*h* - *hm*)/*H* in which *hm* is the height of maximum ionization and *H* is the atmospheric scale height.

We have found that the correction by Eq. (33) shows the best performance for the atmospheric scale height *H* = 70 km. However, when the scale height is too low (e.g., *H* = 60 km) or too large (e.g., *H* = 80 km), its performance degrades especially at low elevation angles (see Fig. 9). Our present investigation shows that its performance can be improved by taking into account the *dI len* dependency on the rate of change of TEC with respect to the elevation angle. In order to find their dependencies, the excess path length has been computed by the ray tracing program considering Chapman profiles with different *H* = 60 and 80 km. The signal frequency *f* = 1227.6 MHz, parameters *hm* = 350 km and *Nm* = 4.96×1012 m-3 are kept constant in each case. The total electron content in the vertical direction is found 123 and 164 TEC units, respectively. The obtained *dlen*, TEC, and the first and second order TEC derivatives with respect to the elevation angle *dTEC*/*dβ* and 398 Global Navigation Satellite Systems – Signal, Theory and Applications

simplified as 0.6577*NmTEC*. Thus, the third order residual term can be approximated by the

The third order term Δ*s*3 will be measured in meters when *f* is measured in Hz and the

There is another formula published by Hoque & Jakowski (2008) for the excess path length

7.5 10 exp 2.13 *len*

*TEC <sup>d</sup> f H hm*

scale height *H* and maximum ionization height *hm* in kilometers and elevation *β* in radians. Comparing both formulas we see that Eq. (33) requires only TEC and elevation information as inputs whereas Eq. (39) additionally requires ionospheric parameters *H* and *hm*.

Both the correction formulas are derived based on simulation studies using Chapman profiles for the ionosphere. The Chapman profile (Rishbeth & Garriott, 1969) has been proved very useful for modeling ionospheric correction. It describes the electron density

where *Nm* is the maximum ionization and *z* = (*h* - *hm*)/*H* in which *hm* is the height of

We have found that the correction by Eq. (33) shows the best performance for the atmospheric scale height *H* = 70 km. However, when the scale height is too low (e.g., *H* = 60 km) or too large (e.g., *H* = 80 km), its performance degrades especially at low elevation angles (see Fig. 9). Our present investigation shows that its performance can be improved by

elevation angle. In order to find their dependencies, the excess path length has been computed by the ray tracing program considering Chapman profiles with different *H* = 60 and 80 km. The signal frequency *f* = 1227.6 MHz, parameters *hm* = 350 km and *Nm* = 4.96×1012 m-3 are kept constant in each case. The total electron content in the vertical direction is found 123 and 164 TEC units, respectively. The obtained *dlen*, TEC, and the first and second order TEC derivatives with respect to the elevation angle *dTEC*/*dβ* and

( )

β

<sup>−</sup> × − <sup>=</sup> (39)

( ) exp(0.5(1 exp( ))) *n h Nm <sup>e</sup>* = − *z z* − − (40)

*len* dependency on the rate of change of TEC with respect to the

5 2

( )

*len* is measured in meters, TEC is in TEC units, frequency *f* in GHz, atmospheric

4 1/8

= (38)

*len* can be computed by Eq. (33).

3 2 2 1 2 534.27 *s N TEC <sup>m</sup> f f*

Δ

maximum ionization *Nm* and TEC in m-3 and electrons/m2, respectively.

As we have seen in the section 3.2, the excess path length *dI*

*I*

distribution *ne* as a function of height *h* in the ionosphere as

maximum ionization and *H* is the atmospheric scale height.

However, these parameters are not easy to estimate in practical cases.

**4.3 New approaches for correcting LoS propagation assumption errors** 

first term only as (using Eqs. 9 and 19)

**4.3.1 Excess path length correction**

computation.

where *dI*

taking into account the *dI*

*d2TEC*/*dβ2* have been plotted as functions of elevation angle in Fig. 9. The *dTEC*/*dβ* has been calculated dividing the TEC difference between two measurement epochs by the corresponding elevation angle difference. Then, *d2TEC*/*dβ2* has been calculated dividing the *dTEC*/*dβ* difference between two measurement epochs by the corresponding elevation angle difference.

Comparing plots in Fig. 9, we see that although the dependency of *dlen* on the *dTEC*/*dβ* is not straight forward, its dependency on the *d2TEC*/*dβ2* is obvious at low elevation angles (< 20°). Thus, the magnitude of the *dlen* depends on the magnitude of TEC as well as on the magnitude of *d2TEC*/*dβ2*. Considering this, functional dependencies have been studied separately for different parameters to develop correction formulas. For this, ray tracing calculation has been carried out for different geometrical and ionospheric conditions varying elevation and Chapman layer parameters *H*, *Nm* and *hm*. Thus, the following formula has been obtained for the *dI len* correction.

Fig. 9. Elevation angle dependence of *dI len*, TEC (see right scale), *dTEC*/*dβ* and *d2TEC*/*dβ2* for the Chapman layer parameter *H* = 60 and 80 km

$$d\_{l}^{len} = \frac{a\_1}{f^4} \left( \frac{1}{\left(1 - a\_2 \cos^2 \beta\right)^{1/2}} - 1 \right) \text{TEC}^2 + a\_3 \left(\frac{d^2 \text{TEC}}{d\beta^2}\right)^2 \cos \beta \tag{41}$$

Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches 401

*TEC TEC*

where Δ*TECbend* is measured in TECU, atmospheric scale height *H* is in km, the maximum ionization height *hm* is in km, signal frequency *f* is in GHz, TEC is in TECU and elevation angle *β* is in radians. Again, it requires the knowledge of the ionospheric parameters *H* and *hm*. If actual parameters are not known, the formula may not be useful in practical purposes. Therefore, in the present work, we have looked for a correction formula depending only on the TEC, elevation angle and second order derivative of TEC with respect to the elevation. We have found that the formula Eq. (41) can be used for such purposes multiplying simply

*bend*

( )

2

1 cos

ray tracing Eq. (42) Eq. (43) Eq. (43) approx.

0 15 30 45

H = 60 km

Fig. 11. Comparison of Δ*TECbend* correction formulas with ray tracing results

1 2

*<sup>c</sup> d TEC TEC TEC c f d c*

β

2 1/2 2 <sup>3</sup> <sup>2</sup>

where *c*1 = 1.2963, *c*2 = 0.8260, *c*3 = 0.0496. The Δ*TECbend* will be computed in TEC units when *β* is measured in radians, *f* is in MHz and TEC is in TEC units and *d*2*TEC*/*dβ*2 in TECU/deg2. The polynomial coefficients are derived based on a nonlinear fit with ray tracing results in

<sup>−</sup> ⎝ ⎠ ⎝ ⎠

Elevation /deg Elevation /deg

The elevation angle dependence of Δ*TECbend* has been plotted in Fig. 11 for the proposed correction formula Eq. (43) as well as for the Eq. (42). Also ray tracing results are plotted for comparisons. Comparing Δ*TECbend* computed by the Eq. (43) and Eq. (42) with ray tracing results, we see that at higher *H* values (e.g., *H* = 80 km) both correction results are comparable. However, the performance of the new approach significantly degrades at lower

As already mentioned, the derivative *d2TEC*/*dβ2* is very sensitive to TEC gradients. Considering this, another set of coefficients have been determined excluding the derivative term in Eq. (43). In this case, we have found that *c*1 = 1.4563, *c*2 = 0.8260 and *c*3 is set to zero.

ΔTEC /TECU

bend

0.3

0.2

0.1

0

⎛ ⎞ ⎛ ⎞ ⎜ ⎟ <sup>=</sup> <sup>−</sup> <sup>+</sup> ⎜ ⎟

Δ

by *f*2 and determining new coefficients.

Δ

least square senses as before.

0.2

0.15

0.1

bend ΔTEC /TECU

0.05

0

*H* values (e.g., *H* = 60 km).

*bend*

3 2

<sup>1</sup> 1 cos

β

<sup>−</sup> × − <sup>=</sup> (42)

<sup>2</sup> <sup>2</sup>

β  β

0 15 30 45

(43)

H = 80 km

2 0.3 1.108 10 exp( 2.1844 )

*f Hhm*

where *a*1 = 2.6123×108, *a*2 = 0.8260, *a*3 = 6.64. The *dI len* will be computed in millimeters when *β* is measured in radians, *f* is in MHz, TEC is in TEC units and *d*2*TEC*/*dβ*2 in TECU/deg2. The polynomial coefficients are derived based on a nonlinear fit with ray tracing results in least square senses.

The elevation angle dependence of *dI len* has been plotted in Fig. 10 using the proposed correction formula Eq. (41) as well as by Eqs. (33) and (39). In addition, ray tracing results are plotted for comparisons. Comparing *dI len* computed by the Eq. (33) and Eq. (41) with ray tracing results, we see that at higher *H* values (e.g., *H* = 80 km) the correction given by the Eq. (41) performs better. However, its performance degrades at lower *H* values (e.g., *H* = 60 km), especially around 7 – 21° elevation angle.

We find that the Eq. (39) gives the best performance. However, it requires ionospheric parameters *H* and *hm* as inputs which are not known to the GNSS users. Inaccurate assumption of ionospheric parameters may give erroneous estimation of *dI len*. We see that at *H* = 80 km the correction given by the new approach Eq. (41) is even comparable to the correction given by Eq. (39).

Fig. 10. Comparison of excess path length correction formulas with ray tracing results

While using the proposed correction Eq. (41), it should be remembered that due to its dependency on the *d2TEC*/*dβ2* term, it is very sensitive to TEC gradients or irregularities. In such cases the correction given by Eq. (33) is recommended for use.

#### **4.3.2 ΔTECbend correction**

Due to the ray path bending satellite signals propagate in curvature paths instead of straight LoS paths. However, TECs along a curvature path and the corresponding LoS path are not the same rather the TEC along the curvature path is slightly larger than the LoS one. The difference between them is defined as the Δ*TECbend* (see Eq. 7). The Δ*TECbend* can be computed by the following formula given by Hoque & Jakowski (2008).

400 Global Navigation Satellite Systems – Signal, Theory and Applications

is measured in radians, *f* is in MHz, TEC is in TEC units and *d*2*TEC*/*dβ*2 in TECU/deg2. The polynomial coefficients are derived based on a nonlinear fit with ray tracing results in least

correction formula Eq. (41) as well as by Eqs. (33) and (39). In addition, ray tracing results

tracing results, we see that at higher *H* values (e.g., *H* = 80 km) the correction given by the Eq. (41) performs better. However, its performance degrades at lower *H* values (e.g., *H* = 60

We find that the Eq. (39) gives the best performance. However, it requires ionospheric parameters *H* and *hm* as inputs which are not known to the GNSS users. Inaccurate

*H* = 80 km the correction given by the new approach Eq. (41) is even comparable to the

4

3

2

Excess Path /cm

1

0

0 15 30 45 0 15 30 45 Elevation /deg Elevation /deg

Fig. 10. Comparison of excess path length correction formulas with ray tracing results

While using the proposed correction Eq. (41), it should be remembered that due to its dependency on the *d2TEC*/*dβ2* term, it is very sensitive to TEC gradients or irregularities. In

Due to the ray path bending satellite signals propagate in curvature paths instead of straight LoS paths. However, TECs along a curvature path and the corresponding LoS path are not the same rather the TEC along the curvature path is slightly larger than the LoS one. The difference between them is defined as the Δ*TECbend* (see Eq. 7). The Δ*TECbend* can be

assumption of ionospheric parameters may give erroneous estimation of *dI*

ray tracing Eq. (33) Eq. (39) Eq. (41)

H = 60 km

such cases the correction given by Eq. (33) is recommended for use.

computed by the following formula given by Hoque & Jakowski (2008).

*len* will be computed in millimeters when *β*

*len* has been plotted in Fig. 10 using the proposed

*len* computed by the Eq. (33) and Eq. (41) with ray

*len*. We see that at

H = 80 km

where *a*1 = 2.6123×108, *a*2 = 0.8260, *a*3 = 6.64. The *dI*

The elevation angle dependence of *dI*

correction given by Eq. (39).

**4.3.2 ΔTECbend correction** 

Excess Path /cm

3

2

1

0

are plotted for comparisons. Comparing *dI*

km), especially around 7 – 21° elevation angle.

square senses.

$$\text{ATEC}\_{bend} = \frac{1.108 \times 10^{-3} \exp(-2.1844 \beta) \text{TEC}^2}{f^2 \text{Hlim}^{0.3}} \tag{42}$$

where Δ*TECbend* is measured in TECU, atmospheric scale height *H* is in km, the maximum ionization height *hm* is in km, signal frequency *f* is in GHz, TEC is in TECU and elevation angle *β* is in radians. Again, it requires the knowledge of the ionospheric parameters *H* and *hm*. If actual parameters are not known, the formula may not be useful in practical purposes. Therefore, in the present work, we have looked for a correction formula depending only on the TEC, elevation angle and second order derivative of TEC with respect to the elevation. We have found that the formula Eq. (41) can be used for such purposes multiplying simply by *f*2 and determining new coefficients.

$$\text{ATEC}\_{bend} = \frac{c\_1}{f^2} \left( \frac{1}{\left(1 - c\_2 \cos^2 \beta \right)^{1/2}} - 1 \right) \text{TEC}^2 + c\_3 \left( \frac{d^2 \text{TEC}}{d\beta^2} \right)^2 \cos \beta \tag{43}$$

where *c*1 = 1.2963, *c*2 = 0.8260, *c*3 = 0.0496. The Δ*TECbend* will be computed in TEC units when *β* is measured in radians, *f* is in MHz and TEC is in TEC units and *d*2*TEC*/*dβ*2 in TECU/deg2. The polynomial coefficients are derived based on a nonlinear fit with ray tracing results in least square senses as before.

Fig. 11. Comparison of Δ*TECbend* correction formulas with ray tracing results

The elevation angle dependence of Δ*TECbend* has been plotted in Fig. 11 for the proposed correction formula Eq. (43) as well as for the Eq. (42). Also ray tracing results are plotted for comparisons. Comparing Δ*TECbend* computed by the Eq. (43) and Eq. (42) with ray tracing results, we see that at higher *H* values (e.g., *H* = 80 km) both correction results are comparable. However, the performance of the new approach significantly degrades at lower *H* values (e.g., *H* = 60 km).

As already mentioned, the derivative *d2TEC*/*dβ2* is very sensitive to TEC gradients. Considering this, another set of coefficients have been determined excluding the derivative term in Eq. (43). In this case, we have found that *c*1 = 1.4563, *c*2 = 0.8260 and *c*3 is set to zero.

Ionospheric Propagation Effects on GNSS Signals and New Correction Approaches 403

Hoque, M. M. & Jakowski, N. (2006). Higher-order ionospheric effects in precise GNSS

Hoque, M. M. & Jakowski, N. (2007). Mitigation of higher order ionospheric effects on GNSS users in Europe. GPS Solut., Vol. 12, No. 2, doi: 10.1007/s10291-007-0069-5 Hoque, M. M. & Jakowski, N. (2008). Estimate of higher order ionospheric errors in GNSS positioning, Radio Sci., Vol. 43, No. RS5008, doi: 10.1029/2007RS003817 Hoque, M. M. & Jakowski, N. (2010a). Higher Order Ionospheric Errors in Modernized

Hoque, M. M. & Jakowski, N. (2010b). Higher order ionospheric propagation effects on GPS

Hoque, M. M. & Jakowski N. (2011). Ionospheric bending correction for GNSS radio

Jakowski, N., Porsch, F. & Mayer, G. (1994). Ionosphere-Induced-Ray-Path Bending Effects

Positionierung, Navigation und Kommunikation, Vol. SPN 1/94, pp. (6-13) Jakowski, N., Mayer, C., Hoque, M. M. & Wilken, V. (2011). TEC Models and Their Use in

Kedar, S., Hajj, G., Wilson, B. & Heflin, M. (2003). The effect of the second order GPS

Kim, B. C. & Tinin M. V. (2007). Contribution of ionospheric irregularities to the error of

Kim, B.C. & Tinin, M.V. (2011). Potentialities of multifrequency ionospheric correction in

Klobuchar, J. A. (1996). Ionospheric Effects on GPS, In: *Global Positioning System: Theory and* 

Leitinger R., Putz, E. (1988). Ionospheric refraction errors and observables, In: *Atmospheric* 

Mandea M. & Macmillan, S. (2000). International Geomagnetic Reference Field—the eighth

Moore, R. C. & Morton, Y.T. (2011). Magneto-ionic polarization and GPS signal propagation

Morton, Y. T., Zhou, Q. & van Graas, F. (2009). Assessment of second-order ionosphere error

occultation signals. Radio Sci, Vol. 46, No. RS0D06, pp. (9),

Ionosphere Monitoring. Radio Sci., doi:10.1029/2010RS004620

Institute of Aeronautics & Astronautics, ISBN 156347106X

generation. *Earth Planets Space*, Vol. 52, No. 12, pp. (1119–1124)

Radio Sci., Vol. 44, No. RS1002, doi:10.1029/2008RS003888

Monograph 12, School of Surveying, UNSW, Sydney

006-0106-0

012-6

doi:10.1016/j.asr.2010.02.013

(1829), DOI 10.1029/2003 GL017639

10.1007/s00190-006-0099-8

10.1007/s00190-010-0425.

doi:10.1029/2010RS004583

positioning. Journal of Geodesy, Vol. 81, No. 4, pp. (259-268), DOI 10.1007/s00190-

GPS and Future Galileo Systems. In: Global Positioning Systems, Asphaug V. & Sorensen E. (Eds.), pp. (1-28), Nova Science Publishers, Inc., ISBN 978-1-60741-

radio occultation signals. J. Adv. Space Res., Vol. 46, No. 2, pp. (162-173),

in Precise Satellite Positioning Systems. Zeitschrift für Satellitengestützte

ionospheric correction on receiver positions. Geophys Res Lett, Vol. 30, No. 16, pp.

dual-frequency GNSS positioning. J Geod, Vol. 81, pp. (189–199), DOI

Global Navigation Satellite Systems. J Geodesy, Vol. 85, No. 3, pp. (159-169), doi:

*Applications*, Vol I, Parkinson, B. W. & Spilker, J. J. (Eds.), pp. (485-515), American

*effects on the geodetic space measurements*, Brunner, F. K. (Ed.), pp. (81-102),

through the ionosphere, Radio Sci, Vol. 46, No. RS1008, doi: 10.1029/2010RS004380

in GPS range observables using Arecibo incoherent scatter radar measurements,

The elevation angle dependency for such an approximation of Eq. (43) is also plotted in Fig. 11.

As already mentioned, after removal of the Selective Availability, the ionosphere becomes the single largest error source for GNSS error budgets. Fortunately, a dual-frequency ionosphere-free combination can remove about 99% of the ionospheric effects; thanks to the dispersive nature of the ionosphere. Although higher order residual terms are less than 1% of the first order term, they can be many centimeters during times of high TEC and represent large errors in geodetic measurements especially in precise point positioning. This chapter gives estimation of higher order ionospheric terms at different level of ionospheric ionization and discusses different correction approaches for them.

#### **5. References**


402 Global Navigation Satellite Systems – Signal, Theory and Applications

The elevation angle dependency for such an approximation of Eq. (43) is also plotted in

As already mentioned, after removal of the Selective Availability, the ionosphere becomes the single largest error source for GNSS error budgets. Fortunately, a dual-frequency ionosphere-free combination can remove about 99% of the ionospheric effects; thanks to the dispersive nature of the ionosphere. Although higher order residual terms are less than 1% of the first order term, they can be many centimeters during times of high TEC and represent large errors in geodetic measurements especially in precise point positioning. This chapter gives estimation of higher order ionospheric terms at different level of ionospheric

Appleton, E. V. (1932). Wireless studies of the ionosphere, *Proceeding of Instn. Elect. Engrs*.,

Bassiri, S. & Hajj, G. A. (1993). Higher-order ionospheric effects on the global positioning

Blewitt, G. (1987). New approaches to GPS carrier-phase ambiguity resolution. Proceeding of XIX General Assembly of the IUGG, Vancouver, Canada, August 10-22 Brunner, F. K. & Gu, M. (1991). An improved model for the dual frequency ionospheric

Budden, K. G. (Ed.). (1985). *The Propagation of Radio Waves: the theory of radio waves of low power in the ionosphere and magnetosphere*, Cambridge University Press, Cambridge. Datta-Barua, S., Walter, T., Blanch, J. & Enge, P. (2008). Bounding higher-order ionosphere

Elizabeth, J. P., Matt, A. K., Philip, M. & David, A. L. (2010). A first look at the effects of

Fritsche, M., Dietrich, R., Knöfel, C., Rülke, A., Vey, S., Rothacher, M. & Steigenberger, P.

Hartmann, G. K. & Leitinger, R. (1984). Range errors due to ionospheric and tropospheric

Hawarey, M., Hobiger, T. & Schuh, H. (2005). Effects of the 2nd order ionospheric terms on

Hernandez-Pajares, M., Jaun, J. M., Sanz, J. & Orus, R. (2007). Second order ionospheric term

system observables and means of modeling them. manuscripta geodaetica, Vol. 18,

correction of GPS observations. manuscripta geodaetica, Vol. 16, No. 3, pp. (205-

errors for the dual-frequency GPS user. Radio Sci., Vol. 43, No. RS5010, pp. (15),

ionospheric signal bending on a globally processed GPS network. J Geod, Vol. 84,

(2005). Impact of higher-order ionospheric terms on GPS estimates. Geophys Res

effects for signal frequencies above 100 MHz. Bull. Geod, Vol 58, No. 2, pp. (109-

VLBI measurements. Geophys Res Lett, Vol. 32, No. 11, L11304, DOI

in GPS: implementation and impact on geodetic estimates. Journal of Geophysical

ionization and discusses different correction approaches for them.

Vol. 7, No. 21, pp. (257-265), 10.1049/pws.1932.0027

Davies, K. (Ed.). (1990). *Ionospheric Radio*, Peter Peregrinus Ltd, London.

Lett, Vol. 32, No. 23, L23311, DOI 10.1029/2005GL024342

Research, Vol. 112, No. B08417, doi:10.1029/2006JB004707

pp. (491–499), DOI 10.1007/s00190-010-0386-2

Fig. 11.

**5. References** 

214)

136)

No. 6, pp. (280-289)

doi:10.1029/2007RS003772

10.1029/2005GL022729


**0**

**17**

*Finland*

**Multipath Mitigation Techniques for**

**Satellite-Based Positioning Applications**

*Department of Communications Engineering, Tampere University of Technology*

The ever-growing public interest on location and positioning services has originated a demand for a high performance Global Navigation Satellite System (GNSS), such as the Global Positioning System (GPS) or the future European satellite navigation system, Galileo. The performance of GNSS is subject to several errors, such as ionosphere delay, troposphere delay, receiver noise and multipath. Among all these errors, multipath is the main limiting factor in precision-oriented GNSS applications. The reception of multipath creates a bias into the time delay estimate of the Delay Lock Loop (DLL) of a conventional navigation receiver, which eventually leads to an error in the receiver's position estimate. In order to mitigate the multipath influence on navigation receivers, the multipath problem has been approached from several directions. Among them, the use of special multipath limiting antennas (i.e., choke ring or multi-beam antennas), the post-processing techniques to reduce carrier multipath, the carrier smoothing to reduce code multipath, and the code tracking techniques based on receiver internal correlation function are the most prominent approaches. In this chapter, the discussion is mainly focused on the correlation-based multipath mitigation techniques at the receiver side; since the correlation-based multipath mitigation approach is by far the most convenient and popular way to deal with multipath error for a stand-alone GNSS receiver. The classical correlation-based code tracking structure used in GNSS is based on a feedback delay estimator and is implemented via a feedback loop. The most known feedback delay estimator is the Early-Minus-Late (EML) DLL technique, where two correlators spaced at one chip from each other are used in the receiver in order to form a discriminator function, whose zero crossings determine the path delays of the received signal Baltersee et al. (2001), Bischoff et al. (2002), Chen & Davisson (1994), Fine & Wilson (1999), Fock et al. (2001), Laxton (1996). The classical EML fails to cope with multipath propagation Dierendonck et al. (1992), Simon et al. (1994). Therefore, several enhanced EML-based techniques have been introduced in the literature during the last two decades in order to mitigate the impact of multipath, especially in closely spaced path scenarios. One class of these enhanced EML techniques is based on the idea of narrowing the spacing between the early and late correlators, i.e., narrow EML (nEML) or narrow correlator Dierendonck et al. (1992), Irsigler & Eissfeller (2003), McGraw & Braasch (1999). The choice of correlator spacing depends on the receiver's available front-end bandwidth along with the associated sampling frequency Betz & Kolodziejski (2000). Correlator spacings in the range of 0.05 to 0.2 chips are

commercially available for nEML based GPS receivers Braasch (2001).

**1. Introduction**

Mohammad Zahidul H. Bhuiyan and Elena Simona Lohan


## **Multipath Mitigation Techniques for Satellite-Based Positioning Applications**

Mohammad Zahidul H. Bhuiyan and Elena Simona Lohan *Department of Communications Engineering, Tampere University of Technology Finland*

#### **1. Introduction**

404 Global Navigation Satellite Systems – Signal, Theory and Applications

Parkinson, B. W. & Gilbert, S. W. (1983). NAVSTAR: Global Positioning System - Ten Years Later. Proc IEEE, Vol. 71, No. 10, pp (1177-1186), 10.1109/PROC.1983.12745 Rishbeth, H. & Garriott O. K. (Eds.). (1969*). Introduction to ionospheric physics*, Academic,

Strangeways, H. J. & Ioannides, R. T. (2002). Rigorous calculation of ionospheric effects on

Geoph Hung, Vol. 37, No. 2-3, pp. (281-292)

GPS Earth-Satellite paths using a precise path determination methods. Acta Geod

New York

The ever-growing public interest on location and positioning services has originated a demand for a high performance Global Navigation Satellite System (GNSS), such as the Global Positioning System (GPS) or the future European satellite navigation system, Galileo. The performance of GNSS is subject to several errors, such as ionosphere delay, troposphere delay, receiver noise and multipath. Among all these errors, multipath is the main limiting factor in precision-oriented GNSS applications. The reception of multipath creates a bias into the time delay estimate of the Delay Lock Loop (DLL) of a conventional navigation receiver, which eventually leads to an error in the receiver's position estimate. In order to mitigate the multipath influence on navigation receivers, the multipath problem has been approached from several directions. Among them, the use of special multipath limiting antennas (i.e., choke ring or multi-beam antennas), the post-processing techniques to reduce carrier multipath, the carrier smoothing to reduce code multipath, and the code tracking techniques based on receiver internal correlation function are the most prominent approaches.

In this chapter, the discussion is mainly focused on the correlation-based multipath mitigation techniques at the receiver side; since the correlation-based multipath mitigation approach is by far the most convenient and popular way to deal with multipath error for a stand-alone GNSS receiver. The classical correlation-based code tracking structure used in GNSS is based on a feedback delay estimator and is implemented via a feedback loop. The most known feedback delay estimator is the Early-Minus-Late (EML) DLL technique, where two correlators spaced at one chip from each other are used in the receiver in order to form a discriminator function, whose zero crossings determine the path delays of the received signal Baltersee et al. (2001), Bischoff et al. (2002), Chen & Davisson (1994), Fine & Wilson (1999), Fock et al. (2001), Laxton (1996). The classical EML fails to cope with multipath propagation Dierendonck et al. (1992), Simon et al. (1994). Therefore, several enhanced EML-based techniques have been introduced in the literature during the last two decades in order to mitigate the impact of multipath, especially in closely spaced path scenarios. One class of these enhanced EML techniques is based on the idea of narrowing the spacing between the early and late correlators, i.e., narrow EML (nEML) or narrow correlator Dierendonck et al. (1992), Irsigler & Eissfeller (2003), McGraw & Braasch (1999). The choice of correlator spacing depends on the receiver's available front-end bandwidth along with the associated sampling frequency Betz & Kolodziejski (2000). Correlator spacings in the range of 0.05 to 0.2 chips are commercially available for nEML based GPS receivers Braasch (2001).

in order to determine accurately the shape of the multipath corrupted correlation function. According to Townsend et al. (1995), MEDLL provides superior long delay multipath mitigation performance compared to nEML at the cost of multi-correlator based tracking

Multipath Mitigation Techniques for Satellite-Based Positioning Applications 407

A new technique to mitigate multipath by means of correlator reference waveform was proposed in Weill (1997). This technique, referred to as Second Derivative correlator, generates a signal correlation function which has a much narrower width than a standard correlation function, and is therefore capable of mitigating multipath errors over a much wider range of secondary path delays. The narrowing of the correlation function is accomplished by using a specially designed code reference waveform (i.e. the negative of the second order derivative of correlation function) instead of the ideal code waveform used in almost all existing receivers. However, this new technique reduces the multipath errors at the expense of a moderate decrease in the effective Signal-to-Noise Ratio (SNR) due to the effect of narrowing the correlation function. A similar strategy, named as Slope Differential (SD), is based on the second order derivative of the correlation function Lee et al. (2006). It is shown in Lee et al. (2006) that this technique has better multipath performance than nEML and Strobe Correlator. However, the performance measure was solely based on the theoretical MEE curves, thus its

A completely different approach to mitigate multipath error is used in NovAtel's recently developed Vision Correlator Fenton & Jones (2005). The Vision Correlator (VC) is based on the concept of Multipath Mitigation Technique (MMT) developed in Weill (2002). It can provide a significant improvement in detecting and removing multipath signals as compared to other standard multipath resistant code tracking algorithms (for example, PAC of NovAtel). However, VC has the shortcoming that it requires a reference function shape to be used to fit the incoming data with the direct path and the secondary path reference signals. The reference function generation has to be accomplished a-priori, and it must incorporate the issues related

Several advanced multipath mitigation techniques were also proposed in Bhuiyan (2011), Granados et al. (2005), Granados & Rubio (2000), Lohan et al. (2006). These techniques, in general, offer better tracking performance than the traditional DLL at a cost of increased complexity. However, the performance of these techniques have not yet been evaluated in

The rest of this chapter is organized as follows. Multipath propagation phenomena is described first, followed by a description on the influence of signal and receiver parameters on multipath error. The following section provides an elaborate description on different multipath mitigation techniques starting from the conventional state-of-the-art techniques to the relatively complex advanced multipath mitigation techniques. An extensive literature review on the contemporary research on multipath mitigation techniques are also provided. The performance evaluations of some of the multipath mitigation techniques are shown in Section 6 via multipath error envelopes and also via simulations in multipath fading channel model. Finally, some general conclusions are drawn based on the discussions provided in

potential benefit in more realistic multipath environment is still an open issue.

to Radio Frequency (RF) distortions introduced by the front-end.

more realistic multipath channel model with real GNSS signals.

structure.

earlier sections.

Another family of discriminator-based DLL variants proposed for GNSS is the so-called Double-Delta (ΔΔ) technique, which uses more than 3 correlators in the tracking loop (typically, 5 correlators: two early, one in-prompt and two late) Irsigler & Eissfeller (2003). The ΔΔ technique offers better multipath rejection in medium-to-long delay multipath Hurskainen et al. (2008), McGraw & Braasch (1999) in good Carrier-to-Noise density ratio (C/N0). Couple of well-known particular cases of ΔΔ technique are the High Resolution Correlator (HRC) McGraw & Braasch (1999), the Strobe Correlator (SC) Garin & Rousseau (1997), Irsigler & Eissfeller (2003), the Pulse Aperture Correlator (PAC) Jones et al. (2004) and the modified correlator reference waveform Irsigler & Eissfeller (2003), Weill (2003). One other similar tracking structure is the Multiple Gate Delay (MGD) correlator Bello & Fante (2005), Bhuiyan (2006), Fante (2003), Fante (2004), where the number of early and late gates and the weighting factors used to combine them in the discriminator are the parameters of the model, and can be optimized according to the multipath profile as illustrated in Hurskainen et al. (2008). While coping better with the ambiguities of Binary Offset Carrier (BOC) correlation function, the MGD provides slightly better performance than the nEML at the expense of higher complexity and is sensitive to the parameters chosen in the discriminator function (i.e., weights, number of correlators and correlator spacing) Bhuiyan (2006), Hurskainen et al. (2008).

Another tracking structure closely related to ΔΔ technique is the Early1 / Early2 (E1/E2) tracker, initially proposed in Dierendonck & Braasch (1997), and later described in Irsigler & Eissfeller (2003). In E1/E2 tracker, the main purpose is to find a tracking point on the correlation function that is not distorted by multipath. As reported in Irsigler & Eissfeller (2003), E1/E2 tracker shows some performance improvement over ΔΔ technique only for very short delay multipath for GPS L1 Coarse / Acquisition (C/A) signal.

Another feedback tracking structure is the Early-Late-Slope (ELS) Irsigler & Eissfeller (2003), which is also known as Multipath Elimination Technique (MET) Townsend & Fenton (1994). The simulation results performed in Irsigler & Eissfeller (2003) showed that ELS is outperformed by HRC with respect to Multipath Error Envelopes (MEEs), for both Binary Phase Shift Keying (BPSK) and Sine BOC(1,1) (SinBOC(1,1)) modulated signals.

A new multipath estimation technique, named as A-Posteriori Multipath Estimation (APME), is proposed in Sleewaegen & Boon (2001), which relies on a-posteriori estimation of the multipath error tracking. Multipath error is estimated independently in a multipath estimator module on the basis of the correlation values from the prompt and very late correlators. According to Sleewaegen & Boon (2001), the multipath performance of GPS L1 C/A signal is comparable with that of the Strobe Correlator: slight improvement for very short delays (i.e., delays less than 20 meters), but rather significant deterioration for medium delays.

In Phelts & Enge (2000a), a fundamentally different approach is adopted to solve the problem of multipath in the context of GNSS. The proposed technique, named as Tracking Error Compensator (TrEC), utilizes the multipath invariant properties of the received correlation function in order to provide significant performance benefits over nEML for narrow-band GPS receivers Phelts & Enge (2000a), Phelts & Enge (2000b).

One of the most promising advanced multipath mitigation techniques is the Multipath Estimating Delay Lock Loop (MEDLL) Nee (1992), Nee et al. (1994), Townsend et al. (1995) implemented by NovAtel for GPS receivers. MEDLL is considered as a significant evolutionary step in the receiver-based attempt to mitigate multipath. It uses many correlators 2 Will-be-set-by-IN-TECH

Another family of discriminator-based DLL variants proposed for GNSS is the so-called Double-Delta (ΔΔ) technique, which uses more than 3 correlators in the tracking loop (typically, 5 correlators: two early, one in-prompt and two late) Irsigler & Eissfeller (2003). The ΔΔ technique offers better multipath rejection in medium-to-long delay multipath Hurskainen et al. (2008), McGraw & Braasch (1999) in good Carrier-to-Noise density ratio (C/N0). Couple of well-known particular cases of ΔΔ technique are the High Resolution Correlator (HRC) McGraw & Braasch (1999), the Strobe Correlator (SC) Garin & Rousseau (1997), Irsigler & Eissfeller (2003), the Pulse Aperture Correlator (PAC) Jones et al. (2004) and the modified correlator reference waveform Irsigler & Eissfeller (2003), Weill (2003). One other similar tracking structure is the Multiple Gate Delay (MGD) correlator Bello & Fante (2005), Bhuiyan (2006), Fante (2003), Fante (2004), where the number of early and late gates and the weighting factors used to combine them in the discriminator are the parameters of the model, and can be optimized according to the multipath profile as illustrated in Hurskainen et al. (2008). While coping better with the ambiguities of Binary Offset Carrier (BOC) correlation function, the MGD provides slightly better performance than the nEML at the expense of higher complexity and is sensitive to the parameters chosen in the discriminator function (i.e., weights, number

of correlators and correlator spacing) Bhuiyan (2006), Hurskainen et al. (2008).

Phase Shift Keying (BPSK) and Sine BOC(1,1) (SinBOC(1,1)) modulated signals.

short delay multipath for GPS L1 Coarse / Acquisition (C/A) signal.

GPS receivers Phelts & Enge (2000a), Phelts & Enge (2000b).

Another tracking structure closely related to ΔΔ technique is the Early1 / Early2 (E1/E2) tracker, initially proposed in Dierendonck & Braasch (1997), and later described in Irsigler & Eissfeller (2003). In E1/E2 tracker, the main purpose is to find a tracking point on the correlation function that is not distorted by multipath. As reported in Irsigler & Eissfeller (2003), E1/E2 tracker shows some performance improvement over ΔΔ technique only for very

Another feedback tracking structure is the Early-Late-Slope (ELS) Irsigler & Eissfeller (2003), which is also known as Multipath Elimination Technique (MET) Townsend & Fenton (1994). The simulation results performed in Irsigler & Eissfeller (2003) showed that ELS is outperformed by HRC with respect to Multipath Error Envelopes (MEEs), for both Binary

A new multipath estimation technique, named as A-Posteriori Multipath Estimation (APME), is proposed in Sleewaegen & Boon (2001), which relies on a-posteriori estimation of the multipath error tracking. Multipath error is estimated independently in a multipath estimator module on the basis of the correlation values from the prompt and very late correlators. According to Sleewaegen & Boon (2001), the multipath performance of GPS L1 C/A signal is comparable with that of the Strobe Correlator: slight improvement for very short delays (i.e., delays less than 20 meters), but rather significant deterioration for medium delays.

In Phelts & Enge (2000a), a fundamentally different approach is adopted to solve the problem of multipath in the context of GNSS. The proposed technique, named as Tracking Error Compensator (TrEC), utilizes the multipath invariant properties of the received correlation function in order to provide significant performance benefits over nEML for narrow-band

One of the most promising advanced multipath mitigation techniques is the Multipath Estimating Delay Lock Loop (MEDLL) Nee (1992), Nee et al. (1994), Townsend et al. (1995) implemented by NovAtel for GPS receivers. MEDLL is considered as a significant evolutionary step in the receiver-based attempt to mitigate multipath. It uses many correlators in order to determine accurately the shape of the multipath corrupted correlation function. According to Townsend et al. (1995), MEDLL provides superior long delay multipath mitigation performance compared to nEML at the cost of multi-correlator based tracking structure.

A new technique to mitigate multipath by means of correlator reference waveform was proposed in Weill (1997). This technique, referred to as Second Derivative correlator, generates a signal correlation function which has a much narrower width than a standard correlation function, and is therefore capable of mitigating multipath errors over a much wider range of secondary path delays. The narrowing of the correlation function is accomplished by using a specially designed code reference waveform (i.e. the negative of the second order derivative of correlation function) instead of the ideal code waveform used in almost all existing receivers. However, this new technique reduces the multipath errors at the expense of a moderate decrease in the effective Signal-to-Noise Ratio (SNR) due to the effect of narrowing the correlation function. A similar strategy, named as Slope Differential (SD), is based on the second order derivative of the correlation function Lee et al. (2006). It is shown in Lee et al. (2006) that this technique has better multipath performance than nEML and Strobe Correlator. However, the performance measure was solely based on the theoretical MEE curves, thus its potential benefit in more realistic multipath environment is still an open issue.

A completely different approach to mitigate multipath error is used in NovAtel's recently developed Vision Correlator Fenton & Jones (2005). The Vision Correlator (VC) is based on the concept of Multipath Mitigation Technique (MMT) developed in Weill (2002). It can provide a significant improvement in detecting and removing multipath signals as compared to other standard multipath resistant code tracking algorithms (for example, PAC of NovAtel). However, VC has the shortcoming that it requires a reference function shape to be used to fit the incoming data with the direct path and the secondary path reference signals. The reference function generation has to be accomplished a-priori, and it must incorporate the issues related to Radio Frequency (RF) distortions introduced by the front-end.

Several advanced multipath mitigation techniques were also proposed in Bhuiyan (2011), Granados et al. (2005), Granados & Rubio (2000), Lohan et al. (2006). These techniques, in general, offer better tracking performance than the traditional DLL at a cost of increased complexity. However, the performance of these techniques have not yet been evaluated in more realistic multipath channel model with real GNSS signals.

The rest of this chapter is organized as follows. Multipath propagation phenomena is described first, followed by a description on the influence of signal and receiver parameters on multipath error. The following section provides an elaborate description on different multipath mitigation techniques starting from the conventional state-of-the-art techniques to the relatively complex advanced multipath mitigation techniques. An extensive literature review on the contemporary research on multipath mitigation techniques are also provided. The performance evaluations of some of the multipath mitigation techniques are shown in Section 6 via multipath error envelopes and also via simulations in multipath fading channel model. Finally, some general conclusions are drawn based on the discussions provided in earlier sections.

• Front-end filter bandwidth (i.e., precorrelation bandwidth),

• Type of discriminator used to run the DLL (i.e., nEML, HRC, etc.),

• Amplitudes, delays and phases of multipath signals with respect to the LOS signal, etc. The type of signal modulation basically determines the shape of the correlation function. For example, BPSK is used to modulate GPS L1 C/A signal, which has a single significant tracking peak within ±1 chip delay from the correct code delay, whereas CBOC modulations (i.e. CBOC(+) for data channel and CBOC(-) for pilot channel) are used to modulate Galileo E1 signal, each of which has more than one significant tracking peak within ±1 chip delay from the correct code delay. Non-coherent (i.e., absolute value of the correlation function) correlation functions for the above modulations are shown in Fig. 2, where the extra peaks

Multipath Mitigation Techniques for Satellite-Based Positioning Applications 409

Correlation shapes for different signal modulations

BPSK CBOC(+) CBOC(−)

−1 −0.5 0 0.5 1

can be clearly observed in case of CBOC modulations. This is the situation in the most ideal single path scenario. The situation would get far worse in the presence of multipath signals, for example, in a typical fading channel model with a two to four path assumption. Fig. 3 shows the distorted correlation shapes of different signal modulations in a two path static channel with path delays [0 0.1] chips and with path powers [0 -6] dB. As seen in Fig 3, the presence of an additional peak (in case of CBOC(+) and CBOC(-)) due to multipath imposes a challenge for the signal acquisition and tracking techniques to lock to the correct peak. If the receiver fails to lock to the correct peak, a multipath error in the order of few tens of meters is

The front-end filter bandwidth used for band-limiting the received signal also has some impact on the correlation shape. The bandwidth, if not chosen sufficiently high, may round off the correlation peak as well as flatten the width of the correlation function, as shown in Fig. 4. For this particular reason, the choice of correlator spacing depends on the receiver's available front-end bandwidth (and off course, on the sampling frequency), that follows the

Fig. 2. Non-coherent correlation functions for different signal modulations.

Code Delay [chips]

• Correlator spacing used in the code tracking,

0

0.2

0.4

0.6

Non−coherent Correlation Function

of no surprise.

0.8

1

• Code chipping rate,

• Number of multipath signals,

### **2. Multipath propagation**

Multipath propagation occurs mostly due to reflected GNSS signals from surfaces (such as buildings, metal surfaces etc.) near the receiver, resulting in one or more secondary propagation paths. These secondary path signals, which are superimposed on the desired direct path signal, always have a longer propagation time and can significantly distort the amplitude and phase of the direct path signal. This eventually leads to a deformation in the correlation function as shown in Fig. 1, where a direct LOS signal is added constructively

Fig. 1. Received correlation function in two path static channel model, path delays: [0 0.5] chips, path amplitudes: [0 -3] dB, in-phase combination.

with an in-phase (i.e., 0◦ phase difference), delayed (0.5 chips delayed) and attenuated (-3 dB attenuated) version of it to form a compound signal. The deformed correlation shape introduces an error bias in the pseudorange measurement that resulted in a degraded positioning performance.

In severe multipath environments like those in dense urban areas, it may be possible that the LOS signal is obstructed completely and only the reflected signals are present. These multipath effects on the code phase measurements are most crucial, and the multipath error can reach up to a few tens of meters, or a couple of hundred at most Gleason & Egziabher (2009). Moreover, unlike other error sources, multipath cannot be reduced through differential processing, since it decorrelates spatially very rapidly. All these issues are the main driving factors for the research conducted in the context of this thesis striving for an optimum correlation-based multipath mitigation technique in terms of mitigation performance as well as implementation complexity.

#### **3. Influence of signal and receiver parameters on multipath error**

The way multipath affects the tracking and navigation performance of a receiver depends on a number of signal and receiver parameters. Among them, the most influential parameters are:

• Type of signal modulation,


4 Will-be-set-by-IN-TECH

Multipath propagation occurs mostly due to reflected GNSS signals from surfaces (such as buildings, metal surfaces etc.) near the receiver, resulting in one or more secondary propagation paths. These secondary path signals, which are superimposed on the desired direct path signal, always have a longer propagation time and can significantly distort the amplitude and phase of the direct path signal. This eventually leads to a deformation in the correlation function as shown in Fig. 1, where a direct LOS signal is added constructively

Impact of multipath on correlation shape

LOS Signal NLOS Signal Combined Signal

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2

Fig. 1. Received correlation function in two path static channel model, path delays: [0 0.5]

with an in-phase (i.e., 0◦ phase difference), delayed (0.5 chips delayed) and attenuated (-3 dB attenuated) version of it to form a compound signal. The deformed correlation shape introduces an error bias in the pseudorange measurement that resulted in a degraded

In severe multipath environments like those in dense urban areas, it may be possible that the LOS signal is obstructed completely and only the reflected signals are present. These multipath effects on the code phase measurements are most crucial, and the multipath error can reach up to a few tens of meters, or a couple of hundred at most Gleason & Egziabher (2009). Moreover, unlike other error sources, multipath cannot be reduced through differential processing, since it decorrelates spatially very rapidly. All these issues are the main driving factors for the research conducted in the context of this thesis striving for an optimum correlation-based multipath mitigation technique in terms of mitigation performance as well

The way multipath affects the tracking and navigation performance of a receiver depends on a number of signal and receiver parameters. Among them, the most influential parameters

**3. Influence of signal and receiver parameters on multipath error**

Code Delay [chips]

**2. Multipath propagation**

0

chips, path amplitudes: [0 -3] dB, in-phase combination.

0.2

0.4

0.6

0.8

Correlation Function

positioning performance.

as implementation complexity.

• Type of signal modulation,

are:

1

1.2

1.4


The type of signal modulation basically determines the shape of the correlation function. For example, BPSK is used to modulate GPS L1 C/A signal, which has a single significant tracking peak within ±1 chip delay from the correct code delay, whereas CBOC modulations (i.e. CBOC(+) for data channel and CBOC(-) for pilot channel) are used to modulate Galileo E1 signal, each of which has more than one significant tracking peak within ±1 chip delay from the correct code delay. Non-coherent (i.e., absolute value of the correlation function) correlation functions for the above modulations are shown in Fig. 2, where the extra peaks

Fig. 2. Non-coherent correlation functions for different signal modulations.

can be clearly observed in case of CBOC modulations. This is the situation in the most ideal single path scenario. The situation would get far worse in the presence of multipath signals, for example, in a typical fading channel model with a two to four path assumption. Fig. 3 shows the distorted correlation shapes of different signal modulations in a two path static channel with path delays [0 0.1] chips and with path powers [0 -6] dB. As seen in Fig 3, the presence of an additional peak (in case of CBOC(+) and CBOC(-)) due to multipath imposes a challenge for the signal acquisition and tracking techniques to lock to the correct peak. If the receiver fails to lock to the correct peak, a multipath error in the order of few tens of meters is of no surprise.

The front-end filter bandwidth used for band-limiting the received signal also has some impact on the correlation shape. The bandwidth, if not chosen sufficiently high, may round off the correlation peak as well as flatten the width of the correlation function, as shown in Fig. 4. For this particular reason, the choice of correlator spacing depends on the receiver's available front-end bandwidth (and off course, on the sampling frequency), that follows the

strongly influences the resulting multipath performance. Generally speaking, a narrower correlator spacing leads to a reduced multipath error and a tracking jitter error, as long as

Multipath Mitigation Techniques for Satellite-Based Positioning Applications 411

The code chipping rate determines the chip length (*Tc*), which ultimately decides the resulting ranging error caused by the multipath. This means that a signal with a larger chip length results in a smaller multipath error contribution. That is why, the modernized GPS L5 signal can offer ten times smaller multipath error contribution than the legacy GPS L1 C/A signal,

The remaining multipath related parameters (i.e., amplitudes, delays, phases and number of multipath signals) depend on the multipath environment, and have direct influence on the tracking performance of the receiver. These parameters are generally used to define a simulation model (for example, multipath fading channel model) in order to analyze the

A unified multi-correlator based delay tracking structure is developed to fit all the multipath mitigation techniques in one common tracking structure. In a multi-correlator based delay tracking structure, a bank of correlators is generated, unlike the conventional DLL-based tracking structure, where only few correlators (i.e., in the range of three to seven complex correlators depending on the type of techniques) are used. This large number of correlators are required by the advanced multipath mitigation techniques in order to estimate the channel properties and to take a decision on the correct LOS code delay. As shown in Fig. 5, after

the necessary front-end processing, and after the carrier has been wiped-off, the received post-processed signal is passed through a bank of correlators. The NCO (Numerically Controlled Oscillator) and PRN generator block produces a bank of early and late versions of replica codes based on the delay of the LOS signal *<sup>τ</sup>*, the correlator spacing <sup>Δ</sup>, and the number of correlators *M*. In case of an EML tracking loop, the corresponding early-late spacing is equal to 2Δ. The received signal is correlated with each replica in the correlator bank, and the output of the correlator bank is a vector of samples in the correlation envelope. Therefore, we obtain the correlation values for the range of ±*M*Δ chips from the prompt correlator, where *M* is the number of correlators and Δ is the correlator spacing between successive correlators. The various code tracking techniques (named as Discriminator in Fig. 5) utilize the correlation

sufficient front-end bandwidth is ensured Dierendonck et al. (1992).

as it has ten times higher chipping rate than that of L1.

performance of different multipath mitigation techniques.

Fig. 5. Block diagram for multi-correlator based DLL implementation.

**3.1 Multi-correlator based delay tracking structure**

Fig. 3. Non-coherent correlation functions for different signal modulations in two path static channel.

Fig. 4. Non-coherent correlation functions for BPSK modulated GPS L1 C/A signal in different front-end bandwidths.

relation: the more the bandwidth, the smaller the correlator spacing. As mentioned in Betz & Kolodziejski (2000), the early-late spacing Δ*EL* (i.e., twice the correlator spacing) is related to the front-end bandwidth (double-sided) *BW* and the code chip rate *fchip* according to the following equation:

$$
\Delta\_{EL} \ge \frac{f\_{clip}}{\mathcal{BW}} \tag{1}
$$

The type of the discriminator and the correlator spacing used to form the discriminator function (i.e., nEML, HRC, SC, etc.) determines the behavior of the code discriminator that 6 Will-be-set-by-IN-TECH

Correlation shapes in two path static channel

BPSK CBOC(+) CBOC(−)

> 16 MHz 8 MHz 4 MHz 2 MHz

*BW* (1)

−1 −0.5 0 0.5 1

Fig. 3. Non-coherent correlation functions for different signal modulations in two path static

−1 −0.5 0 0.5 1

relation: the more the bandwidth, the smaller the correlator spacing. As mentioned in Betz & Kolodziejski (2000), the early-late spacing Δ*EL* (i.e., twice the correlator spacing) is related to the front-end bandwidth (double-sided) *BW* and the code chip rate *fchip* according to the

<sup>Δ</sup>*EL* <sup>≥</sup> *fchip*

The type of the discriminator and the correlator spacing used to form the discriminator function (i.e., nEML, HRC, SC, etc.) determines the behavior of the code discriminator that

Fig. 4. Non-coherent correlation functions for BPSK modulated GPS L1 C/A signal in

Code Delay [chips]

Code Delay [chips]

Correlation shapes for different front−end bandwidths

0

0

0.2

0.4

0.6

Non−coherent Correlation Function

different front-end bandwidths.

following equation:

0.8

1

0.2

0.4

0.6

Non−coherent Correlation Function

channel.

0.8

1

strongly influences the resulting multipath performance. Generally speaking, a narrower correlator spacing leads to a reduced multipath error and a tracking jitter error, as long as sufficient front-end bandwidth is ensured Dierendonck et al. (1992).

The code chipping rate determines the chip length (*Tc*), which ultimately decides the resulting ranging error caused by the multipath. This means that a signal with a larger chip length results in a smaller multipath error contribution. That is why, the modernized GPS L5 signal can offer ten times smaller multipath error contribution than the legacy GPS L1 C/A signal, as it has ten times higher chipping rate than that of L1.

The remaining multipath related parameters (i.e., amplitudes, delays, phases and number of multipath signals) depend on the multipath environment, and have direct influence on the tracking performance of the receiver. These parameters are generally used to define a simulation model (for example, multipath fading channel model) in order to analyze the performance of different multipath mitigation techniques.

#### **3.1 Multi-correlator based delay tracking structure**

A unified multi-correlator based delay tracking structure is developed to fit all the multipath mitigation techniques in one common tracking structure. In a multi-correlator based delay tracking structure, a bank of correlators is generated, unlike the conventional DLL-based tracking structure, where only few correlators (i.e., in the range of three to seven complex correlators depending on the type of techniques) are used. This large number of correlators are required by the advanced multipath mitigation techniques in order to estimate the channel properties and to take a decision on the correct LOS code delay. As shown in Fig. 5, after

Fig. 5. Block diagram for multi-correlator based DLL implementation.

the necessary front-end processing, and after the carrier has been wiped-off, the received post-processed signal is passed through a bank of correlators. The NCO (Numerically Controlled Oscillator) and PRN generator block produces a bank of early and late versions of replica codes based on the delay of the LOS signal *<sup>τ</sup>*, the correlator spacing <sup>Δ</sup>, and the number of correlators *M*. In case of an EML tracking loop, the corresponding early-late spacing is equal to 2Δ. The received signal is correlated with each replica in the correlator bank, and the output of the correlator bank is a vector of samples in the correlation envelope. Therefore, we obtain the correlation values for the range of ±*M*Δ chips from the prompt correlator, where *M* is the number of correlators and Δ is the correlator spacing between successive correlators. The various code tracking techniques (named as Discriminator in Fig. 5) utilize the correlation

(2003). ΔΔ technique offers better multipath rejection in medium-to-long delay multipath in good C/N0 Hurskainen et al. (2008); McGraw & Braasch (1999). Couple of well-known particular cases of ΔΔ technique are the High Resolution Correlator (HRC) McGraw & Braasch (1999), the Strobe Correlator (SC) Garin & Rousseau (1997); Irsigler & Eissfeller (2003), the Pulse Aperture Correlator (PAC) Jones et al. (2004) and the modified correlator reference waveform Irsigler & Eissfeller (2003); Weill (2003). One other similar tracking structure is the Multiple Gate Delay (MGD) correlator Bello & Fante (2005); Bhuiyan (2006); Fante (2003; 2004); Jie (2010), where the number of early and late gates and the weighting factors used to combine them in the discriminator are the parameters of the model, and can be optimized according to the multipath profile as illustrated in Hurskainen et al. (2008); Jie (2010). While coping better with the ambiguities of BOC correlation function, the MGD provides slightly better performance than the nEML at the expense of higher complexity and is sensitive to the parameters chosen in the discriminator function (i.e., weights, number of correlators and correlator spacing) Bhuiyan (2006); Hurskainen et al. (2008); Jie (2010). In Hurskainen et al. (2008), it is also shown that ΔΔ technique is a particular case of MGD implementation.

Multipath Mitigation Techniques for Satellite-Based Positioning Applications 413

Another feedback tracking structure is the Early-Late-Slope (ELS) Irsigler & Eissfeller (2003), which is also known as Multipath Elimination Technique (MET) Townsend & Fenton (1994). The ELS is based on two correlator pairs at both sides of the correlation function's central peak with parameterized spacing. Once both slopes are known, they can be used to compute a pseudorange correction that can be applied to the pseudorange measurement. However, simulation results performed in Irsigler & Eissfeller (2003) showed that ELS is outperformed by HRC with respect to Multipath Error Envelopes (MEEs), for both BPSK and SinBOC(1,1) modulated signals. An Improved ELS (IELS) technique was proposed by the Author in Bhuiyan et al. (2008), which introduced two enhancements to the basic ELS approach. The first enhancement was the adaptation of random spacing between the early and the late correlator pairs, while the later one was the utilization of feedforward information in order to determine the most appropriate peak on which the IELS technique should be applied. It was shown in Bhuiyan et al. (2008) that IELS performed better than nEML only in good C/N0 for BPSK and SinBOC(1,1) modulated signals in case of short-delay multipath, but still had

A new multipath estimation technique, named as A-Posteriori Multipath Estimation (APME), is proposed in Sleewaegen & Boon (2001), which relies on a-posteriori estimation of multipath error. Multipath error is estimated independently in a multipath estimator module on the basis of the correlation values from the prompt and very late correlators. The performance of APME in multipath environment is comparable with that of the Strobe Correlator: a slight improvement for very short delays (i.e., delays less than 20 meters), but rather significant

One of the most promising state-of-art multipath mitigation techniques is the Multipath Estimating Delay Lock Loop (MEDLL) Nee (1992); Nee et al. (1994); Townsend et al. (1995)

**4.3 Early-late-slope**

poorer performance than HRC.

**4.4 A-Posteriori multipath estimation**

**4.5 Multipath estimating delay lock loop**

deterioration for medium delays Sleewaegen & Boon (2001).

values as input, and generate the estimated LOS delay as output, which is then smoothed by a loop filter. A loop filter is generally used to improve the code delay estimate, reducing the noise present at the output of the discriminator, and to follow the signal dynamics. The order of a loop filter determines the ability of the filter to respond to different types of dynamics, whereas the filter bandwidth ensures that a low bandwidth leads to a good filtering with a high amount of noise filtered, but also requires that the dynamics of the signals are not too high. The loop bandwidth is usually determined by the coefficients of the filter, and can thus be considered as a design parameter for the filter. In accordance with Kaplan & Hegarty (2006), the code loop filter is a 1st order filter, which can be modeled as:

$$
\hat{\pi}(k+1) = \hat{\pi}(k) + \omega\_0 d(k), \tag{2}
$$

where *ω*<sup>0</sup> is calculated based on the loop filter bandwidth, *Bn*. As shown in Bhuiyan & Lohan (2010), the multi-correlator based tracking structure being used by the advanced multipath mitigation techniques, offers a superior tracking performance to the traditional nEML DLL at the cost of higher number of correlators.

#### **4. State-of-the-art multipath mitigation techniques**

The GNSS community has started the correlation-based multipath mitigation studies since early 1990s with the advent of Narrow Correlator (NC) or narrow Early-Minus-Late (nEML) DLL Dierendonck et al. (1992). This section highlights some of the most prominent state-of-the-art techniques, which have gained a lot of interest in the research community by now.

#### **4.1 Early-minus-late delay lock loop**

The classical correlation-based code tracking structure used in a GNSS receiver is based on a feedback delay estimator and is implemented via a feedback loop. The most known feedback delay estimator is the Early-Minus-Late (EML) DLL, where two correlators spaced at one chip from each other, are used in the receiver in order to form a discriminator function, whose zero crossings determine the path delays of the received signal Baltersee et al. (2001); Bischoff et al. (2002); Chen & Davisson (1994); Fine & Wilson (1999); Fock et al. (2001); Laxton (1996); Lohan (2003). The classical EML usually fails to cope with multipath propagation Dierendonck et al. (1992). Therefore, several enhanced EML-based techniques have been introduced in the literature for last two decades in order to mitigate the impact of multipath, especially in closely spaced path scenarios. A first approach to reduce the influences of code multipath is based on the idea of narrowing the spacing between the early and late correlators, i.e., nEML or narrow correlator Dierendonck et al. (1992); Fenton (1995); Fenton & Dierendonck (1996). The choice of correlator spacing depends on the receiver's available front-end bandwidth along with the associated sampling frequency Betz & Kolodziejski (2000). Correlator spacings in the range of 0.05 to 0.2 chips are commercially available for nEML based GPS receivers Braasch (2001).

#### **4.2 Double delta (**ΔΔ**) technique**

Another family of discriminator-based DLL variants proposed for GNSS receivers is the so-called Double Delta (ΔΔ) technique, which uses more than three correlators in the tracking loop (typically, five correlators: two early, one in-prompt and two late) Irsigler & Eissfeller (2003). ΔΔ technique offers better multipath rejection in medium-to-long delay multipath in good C/N0 Hurskainen et al. (2008); McGraw & Braasch (1999). Couple of well-known particular cases of ΔΔ technique are the High Resolution Correlator (HRC) McGraw & Braasch (1999), the Strobe Correlator (SC) Garin & Rousseau (1997); Irsigler & Eissfeller (2003), the Pulse Aperture Correlator (PAC) Jones et al. (2004) and the modified correlator reference waveform Irsigler & Eissfeller (2003); Weill (2003). One other similar tracking structure is the Multiple Gate Delay (MGD) correlator Bello & Fante (2005); Bhuiyan (2006); Fante (2003; 2004); Jie (2010), where the number of early and late gates and the weighting factors used to combine them in the discriminator are the parameters of the model, and can be optimized according to the multipath profile as illustrated in Hurskainen et al. (2008); Jie (2010). While coping better with the ambiguities of BOC correlation function, the MGD provides slightly better performance than the nEML at the expense of higher complexity and is sensitive to the parameters chosen in the discriminator function (i.e., weights, number of correlators and correlator spacing) Bhuiyan (2006); Hurskainen et al. (2008); Jie (2010). In Hurskainen et al. (2008), it is also shown that ΔΔ technique is a particular case of MGD implementation.

### **4.3 Early-late-slope**

8 Will-be-set-by-IN-TECH

values as input, and generate the estimated LOS delay as output, which is then smoothed by a loop filter. A loop filter is generally used to improve the code delay estimate, reducing the noise present at the output of the discriminator, and to follow the signal dynamics. The order of a loop filter determines the ability of the filter to respond to different types of dynamics, whereas the filter bandwidth ensures that a low bandwidth leads to a good filtering with a high amount of noise filtered, but also requires that the dynamics of the signals are not too high. The loop bandwidth is usually determined by the coefficients of the filter, and can thus be considered as a design parameter for the filter. In accordance with Kaplan & Hegarty

where *ω*<sup>0</sup> is calculated based on the loop filter bandwidth, *Bn*. As shown in Bhuiyan & Lohan (2010), the multi-correlator based tracking structure being used by the advanced multipath mitigation techniques, offers a superior tracking performance to the traditional nEML DLL at

The GNSS community has started the correlation-based multipath mitigation studies since early 1990s with the advent of Narrow Correlator (NC) or narrow Early-Minus-Late (nEML) DLL Dierendonck et al. (1992). This section highlights some of the most prominent state-of-the-art techniques, which have gained a lot of interest in the research community by

The classical correlation-based code tracking structure used in a GNSS receiver is based on a feedback delay estimator and is implemented via a feedback loop. The most known feedback delay estimator is the Early-Minus-Late (EML) DLL, where two correlators spaced at one chip from each other, are used in the receiver in order to form a discriminator function, whose zero crossings determine the path delays of the received signal Baltersee et al. (2001); Bischoff et al. (2002); Chen & Davisson (1994); Fine & Wilson (1999); Fock et al. (2001); Laxton (1996); Lohan (2003). The classical EML usually fails to cope with multipath propagation Dierendonck et al. (1992). Therefore, several enhanced EML-based techniques have been introduced in the literature for last two decades in order to mitigate the impact of multipath, especially in closely spaced path scenarios. A first approach to reduce the influences of code multipath is based on the idea of narrowing the spacing between the early and late correlators, i.e., nEML or narrow correlator Dierendonck et al. (1992); Fenton (1995); Fenton & Dierendonck (1996). The choice of correlator spacing depends on the receiver's available front-end bandwidth along with the associated sampling frequency Betz & Kolodziejski (2000). Correlator spacings in the range of 0.05 to 0.2 chips are commercially available for nEML based GPS receivers Braasch (2001).

Another family of discriminator-based DLL variants proposed for GNSS receivers is the so-called Double Delta (ΔΔ) technique, which uses more than three correlators in the tracking loop (typically, five correlators: two early, one in-prompt and two late) Irsigler & Eissfeller

*<sup>τ</sup>*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) = *<sup>τ</sup>*(*k*) + *<sup>ω</sup>*0*d*(*k*), (2)

(2006), the code loop filter is a 1st order filter, which can be modeled as:

the cost of higher number of correlators.

**4.1 Early-minus-late delay lock loop**

**4.2 Double delta (**ΔΔ**) technique**

now.

**4. State-of-the-art multipath mitigation techniques**

Another feedback tracking structure is the Early-Late-Slope (ELS) Irsigler & Eissfeller (2003), which is also known as Multipath Elimination Technique (MET) Townsend & Fenton (1994). The ELS is based on two correlator pairs at both sides of the correlation function's central peak with parameterized spacing. Once both slopes are known, they can be used to compute a pseudorange correction that can be applied to the pseudorange measurement. However, simulation results performed in Irsigler & Eissfeller (2003) showed that ELS is outperformed by HRC with respect to Multipath Error Envelopes (MEEs), for both BPSK and SinBOC(1,1) modulated signals. An Improved ELS (IELS) technique was proposed by the Author in Bhuiyan et al. (2008), which introduced two enhancements to the basic ELS approach. The first enhancement was the adaptation of random spacing between the early and the late correlator pairs, while the later one was the utilization of feedforward information in order to determine the most appropriate peak on which the IELS technique should be applied. It was shown in Bhuiyan et al. (2008) that IELS performed better than nEML only in good C/N0 for BPSK and SinBOC(1,1) modulated signals in case of short-delay multipath, but still had poorer performance than HRC.

#### **4.4 A-Posteriori multipath estimation**

A new multipath estimation technique, named as A-Posteriori Multipath Estimation (APME), is proposed in Sleewaegen & Boon (2001), which relies on a-posteriori estimation of multipath error. Multipath error is estimated independently in a multipath estimator module on the basis of the correlation values from the prompt and very late correlators. The performance of APME in multipath environment is comparable with that of the Strobe Correlator: a slight improvement for very short delays (i.e., delays less than 20 meters), but rather significant deterioration for medium delays Sleewaegen & Boon (2001).

#### **4.5 Multipath estimating delay lock loop**

One of the most promising state-of-art multipath mitigation techniques is the Multipath Estimating Delay Lock Loop (MEDLL) Nee (1992); Nee et al. (1994); Townsend et al. (1995)

**5.2 Second derivative correlator**

**5.3 Peak tracking**

fading channel.

**5.4 Teager Kaiser operator**

defined as Hamila et al. (2003):

A new technique to mitigate multipath by means of correlator reference waveform was proposed in Weill (1997). This technique, referred to as Second Derivative correlator, generates a signal correlation function which has a much narrower width than a standard correlation function, and is therefore capable of mitigating multipath errors over a much wider range of secondary path delays. The narrowing of the correlation function is accomplished by using a specially designed code reference waveform (i.e. the negative of the second order derivative of correlation function) instead of the ideal code waveform used in almost all existing receivers. However, this new technique reduces the multipath errors at the expense of a moderate decrease in the effective Signal-to-Noise Ratio (SNR) due to the effect of narrowing the correlation function. A similar strategy, named as Slope Differential (SD), is based on the second order derivative of the correlation function Lee et al. (2006). It is shown in Lee et al. (2006) that this technique has better multipath performance than nEML and Strobe Correlator. However, the performance measure was solely based on the theoretical MEE curves, thus its

Multipath Mitigation Techniques for Satellite-Based Positioning Applications 415

Peak Tracking (PT) based techniques, namely PT based on 2nd order Differentiation (PT(Diff2)) and PT based on Teager Kaiser (PT(TK)), were proposed in Bhuiyan et al. (2008). Both the techniques utilize the adaptive thresholds computed from the estimated noise variance of the channel in order to decide on the correct code delay. The adaptive thresholds are computed according to the equations given in Bhuiyan et al. (2008). After that, the advanced techniques generate the competitive peaks which are above the computed adaptive thresholds. The generation of competitive peaks for PT(Diff2) technique is shown in Fig. 6 in two path Nakagami-m fading channel. For each of the competitive peak, a decision variable is formed based on the peak power, the peak position and the delay difference of the peak from the previous delay estimate. Finally, the PT techniques select the peak which has the maximum weight as being the best LOS candidate. It was shown in Bhuiyan et al. (2008) that PT(Diff2) has superior multipath mitigation performance over PT(TK) in two to five path Nakagami-*m*

The Teager Kaiser based delay estimation technique is based on the principle of extracting the signal energy of various channel paths via a nonlinear TK operator Hamila (2002), Hamila et al. (2003). The output Ψ*TK*(*x*(*n*)) of TK operator applied to a discrete signal *x*(*n*), can be

If a non-coherent correlation function is used as an input to the nonlinear TK operator, it can then signal the presence of a multipath component more clearly than looking directly at the correlation function. At least three correlation values (in-prompt, early and very early) are required to compute TK operation. But usually, TK based delay estimation utilizes a higher number of correlators (for example, 193 correlators were used in the simulations reported in

− 1 2

Ψ*TK*(*x*(*n*)) = *x*(*n* − 1)*x*∗(*n* − 1) (3)

[*x*(*n* − 2)*x*∗(*n*) + *x*(*n*)*x*∗(*n* − 2)]

potential benefit in more realistic multipath environment is still an open issue.

implemented by NovAtel for GPS receivers. MEDLL uses several correlators per channel in order to determine accurately the shape of the multipath-corrupted correlation function. Then, a reference correlation function is used in a software module in order to determine the best combination of LOS and NLOS components (i.e., amplitudes, delays, phases and number of multipath). An important aspect of MEDLL is an accurate reference correlation function that can be constructed by averaging the measured correlation functions over a significant amount of total averaging time Nee et al. (1994). However, MEDLL provides superior multipath mitigation performance than nEML at a cost of expensive multi-correlator based delay tracking structure.

#### **4.6 Vision correlator**

A completely different approach to mitigate multipath error is used in NovAtel's recently developed Vision Correlator Fenton & Jones (2005). The Vision Correlator (VC) is based on the concept of Multipath Mitigation Technique (MMT) developed in Weill (2002). It can provide a significant improvement in detecting and removing multipath signals as compared to other standard multipath resistant code tracking algorithms (for example, PAC of NovAtel). However, VC has the shortcoming that it requires a reference function shape to be used to fit the incoming data with the direct path and the secondary path reference signals. The reference function generation has to be accomplished a-priori, and it must incorporate the issues related to Radio Frequency (RF) distortions introduced by the front-end.

## **5. Advanced multipath mitigation techniques**

The advanced multipath mitigation techniques generally require a significant number of correlators (i.e., tens of correlators) in order to estimate the channel characteristics, which are then used to mitigate the multipath effect. Some of the most promising advanced multipath mitigation techniques are presented in the following sub-sections. In most occasions, references are made to the original publications in order to avoid too many technical details.

#### **5.1 Non-coherent multipath estimating delay lock loop**

MEDLL is considered as a significant evolutionary step in the correlation-based multipath mitigation approach. Moreover, MEDLL has stimulated the design of different maximum likelihood based implementations for multipath mitigation. One such variant is the non-coherent MEDLL, developed by the authors, as described in Bhuiyan et al. (2008). The classical MEDLL is based on a maximum likelihood search, which is computationally extensive. The authors implemented a non-coherent version of MEDLL that reduces the search space by incorporating a phase search unit, based on the statistical distribution of multipath phases. It was shown in Bhuiyan et al. (2008) that the performance of this suggested approach depends on the number of random phases considered; meaning that the larger the number is, the better the performance will be. But this will also increase the processing burden significantly. The results reported in Bhuiyan et al. (2008), show that the non-coherent MEDLL provides very good performance in terms of RMSE, but has a rather poor Mean-Time-to-Lose-Lock (MTLL) as compared to the conventional DLL techniques.

#### **5.2 Second derivative correlator**

10 Will-be-set-by-IN-TECH

implemented by NovAtel for GPS receivers. MEDLL uses several correlators per channel in order to determine accurately the shape of the multipath-corrupted correlation function. Then, a reference correlation function is used in a software module in order to determine the best combination of LOS and NLOS components (i.e., amplitudes, delays, phases and number of multipath). An important aspect of MEDLL is an accurate reference correlation function that can be constructed by averaging the measured correlation functions over a significant amount of total averaging time Nee et al. (1994). However, MEDLL provides superior multipath mitigation performance than nEML at a cost of expensive multi-correlator

A completely different approach to mitigate multipath error is used in NovAtel's recently developed Vision Correlator Fenton & Jones (2005). The Vision Correlator (VC) is based on the concept of Multipath Mitigation Technique (MMT) developed in Weill (2002). It can provide a significant improvement in detecting and removing multipath signals as compared to other standard multipath resistant code tracking algorithms (for example, PAC of NovAtel). However, VC has the shortcoming that it requires a reference function shape to be used to fit the incoming data with the direct path and the secondary path reference signals. The reference function generation has to be accomplished a-priori, and it must incorporate the issues related

The advanced multipath mitigation techniques generally require a significant number of correlators (i.e., tens of correlators) in order to estimate the channel characteristics, which are then used to mitigate the multipath effect. Some of the most promising advanced multipath mitigation techniques are presented in the following sub-sections. In most occasions, references are made to the original publications in order to avoid too many technical details.

MEDLL is considered as a significant evolutionary step in the correlation-based multipath mitigation approach. Moreover, MEDLL has stimulated the design of different maximum likelihood based implementations for multipath mitigation. One such variant is the non-coherent MEDLL, developed by the authors, as described in Bhuiyan et al. (2008). The classical MEDLL is based on a maximum likelihood search, which is computationally extensive. The authors implemented a non-coherent version of MEDLL that reduces the search space by incorporating a phase search unit, based on the statistical distribution of multipath phases. It was shown in Bhuiyan et al. (2008) that the performance of this suggested approach depends on the number of random phases considered; meaning that the larger the number is, the better the performance will be. But this will also increase the processing burden significantly. The results reported in Bhuiyan et al. (2008), show that the non-coherent MEDLL provides very good performance in terms of RMSE, but has a rather poor Mean-Time-to-Lose-Lock (MTLL) as compared to the conventional DLL techniques.

to Radio Frequency (RF) distortions introduced by the front-end.

**5. Advanced multipath mitigation techniques**

**5.1 Non-coherent multipath estimating delay lock loop**

based delay tracking structure.

**4.6 Vision correlator**

A new technique to mitigate multipath by means of correlator reference waveform was proposed in Weill (1997). This technique, referred to as Second Derivative correlator, generates a signal correlation function which has a much narrower width than a standard correlation function, and is therefore capable of mitigating multipath errors over a much wider range of secondary path delays. The narrowing of the correlation function is accomplished by using a specially designed code reference waveform (i.e. the negative of the second order derivative of correlation function) instead of the ideal code waveform used in almost all existing receivers. However, this new technique reduces the multipath errors at the expense of a moderate decrease in the effective Signal-to-Noise Ratio (SNR) due to the effect of narrowing the correlation function. A similar strategy, named as Slope Differential (SD), is based on the second order derivative of the correlation function Lee et al. (2006). It is shown in Lee et al. (2006) that this technique has better multipath performance than nEML and Strobe Correlator. However, the performance measure was solely based on the theoretical MEE curves, thus its potential benefit in more realistic multipath environment is still an open issue.

#### **5.3 Peak tracking**

Peak Tracking (PT) based techniques, namely PT based on 2nd order Differentiation (PT(Diff2)) and PT based on Teager Kaiser (PT(TK)), were proposed in Bhuiyan et al. (2008). Both the techniques utilize the adaptive thresholds computed from the estimated noise variance of the channel in order to decide on the correct code delay. The adaptive thresholds are computed according to the equations given in Bhuiyan et al. (2008). After that, the advanced techniques generate the competitive peaks which are above the computed adaptive thresholds. The generation of competitive peaks for PT(Diff2) technique is shown in Fig. 6 in two path Nakagami-m fading channel. For each of the competitive peak, a decision variable is formed based on the peak power, the peak position and the delay difference of the peak from the previous delay estimate. Finally, the PT techniques select the peak which has the maximum weight as being the best LOS candidate. It was shown in Bhuiyan et al. (2008) that PT(Diff2) has superior multipath mitigation performance over PT(TK) in two to five path Nakagami-*m* fading channel.

#### **5.4 Teager Kaiser operator**

The Teager Kaiser based delay estimation technique is based on the principle of extracting the signal energy of various channel paths via a nonlinear TK operator Hamila (2002), Hamila et al. (2003). The output Ψ*TK*(*x*(*n*)) of TK operator applied to a discrete signal *x*(*n*), can be defined as Hamila et al. (2003):

$$\begin{aligned} \Psi\_{TK}(\mathbf{x}(n)) &= \mathbf{x}(n-1)\mathbf{x}^\*(n-1) \\ &- \frac{1}{2} [\mathbf{x}(n-2)\mathbf{x}^\*(n) + \mathbf{x}(n)\mathbf{x}^\*(n-2)] \end{aligned} \tag{3}$$

If a non-coherent correlation function is used as an input to the nonlinear TK operator, it can then signal the presence of a multipath component more clearly than looking directly at the correlation function. At least three correlation values (in-prompt, early and very early) are required to compute TK operation. But usually, TK based delay estimation utilizes a higher number of correlators (for example, 193 correlators were used in the simulations reported in

Δ) according to the following equation:

at a C/N0 of 50 dB-Hz.

0

0.2

0.4

0.6

**Peak Threshold**

**Normalized Correlation Function**

0.8

1

1.2

paths only for the middle delay index (i.e., ( *<sup>M</sup>*+<sup>1</sup>

*<sup>τ</sup><sup>W</sup>* <sup>=</sup> <sup>±</sup>(*<sup>M</sup>* <sup>−</sup> <sup>1</sup>)

Multipath Mitigation Techniques for Satellite-Based Positioning Applications 417

For example, in Bhuiyan & Lohan (2010), 193 correlators were used with a correlator spacing of 0.0208 chips, resulting in a code delay window range of ±2 chips with respect to prompt correlator. Therefore, the LOS delay estimate can be anywhere within this ±2 chips window range. The ideal non-coherent reference correlation functions are generated for up to *L*max

delay index is 97). These ideal correlation functions for the middle delay index are generated off-line and saved in a look-up table in memory. In real-time, RSSML reads the correlation values from the look-up table, translates the ideal reference correlation functions at the middle delay index to the corresponding candidate delay index within the code delay window, and then computes the Minimum Mean Square Error (MMSE) for that specific delay candidate. Instead of considering all possible LOS delays within a predefined code delay window as delay candidates, the search space is first reduced to some competitive peaks which are generated based on the computed noise thresholds as explained in Bhuiyan & Lohan (2010). This will eventually reduce the processing time required to compute the MMSE (i.e., MMSE needs to be computed only for the reduced search space). An example is shown in Fig. 7, where RSSML estimates a best-fitted non-coherent correlation function at a cost of 3.6 · <sup>10</sup>−<sup>4</sup> MMSE in a two-path Rayleigh channel with path delays [0 0.35] chips, path powers [0 -2] dB

**2 Path Rayleigh Channel, Path Delays:[ 0 0.35 ] Chips, CNR: 50 dBHz**

−2 −1.5 −1 −0.5 <sup>0</sup> 0.5 <sup>1</sup> 1.5 <sup>2</sup> −0.2

Fig. 7. Estimated and received non-coherent correlation functions in two path Rayleigh

channel, path delay: [0 0.35] chips, path power: [0 -2] dB, C/N0: 50 dB-Hz.

**Code Delays [Chips]**

<sup>2</sup> <sup>Δ</sup> (4)

<sup>2</sup> )-th delay index; for *M* = 193, the middle

**Received Correlation Function Estimated Correlation Function**

**Residual Error**

Two Path Nakagami−m fading channel, path delays: [0 0.75122] chips, C/No :100dB−Hz

Fig. 6. Generation of competitive peaks for PT(Diff2) technique.

Bhuiyan & Lohan (2010)) and is sensitive to the noise dependent threshold choice. Firstly, it computes the noise variance, which is then used to compute an adaptive threshold. The peaks which are above the adaptive threshold are considered as competitive peaks. Among all the competitive peaks, TK selects the delay associated to that competitive peak which has the closest delay difference from the previous delay estimate.

#### **5.5 Reduced Search Space Maximum Likelihood delay estimator**

A Reduced Search Space Maximum Likelihood (RSSML) delay estimator is another good example of maximum likelihood based approach, which is capable of mitigating the multipath effects reasonably well at the expense of increased complexity. The RSSML, proposed by the Author in Bhuiyan et al. (2009) and then further enhanced in Bhuiyan & Lohan (2010), attempts to compensate the multipath error contribution by performing a nonlinear curve fit on the input correlation function which finds a perfect match from a set of ideal reference correlation functions with certain amplitude(s), phase(s) and delay(s) of the multipath signal. Conceptually, a conventional spread spectrum receiver does the same thing, but for only one signal (i.e., the LOS signal). With the presence of multipath signal, RSSML tries to separate the LOS component from the combined signal by estimating all the signal parameters in a maximum likelihood sense, which consequently achieves the best curve fit on the received input correlation function. As mentioned in Bhuiyan & Lohan (2010), it also incorporates a threshold-based peak detection method, which eventually reduces the code delay search space significantly. However, the downfall of RSSML is the memory requirement which it uses to store the reference correlation functions.

In a multi-correlator based structure, the estimated LOS delay, theoretically, can be anywhere within the code delay window range of ±*τ<sup>W</sup>* chips, though in practice, it is quite likely to have a delay error around the previous delay estimate. The code delay window range essentially depends on the number of correlators (i.e., *M*) and the spacing between the correlators (i.e., Δ) according to the following equation:

12 Will-be-set-by-IN-TECH

:100dB−Hz

Correlation Function

Diff2 Diff2 Threshold Noise Threshold

Two Path Nakagami−m fading channel, path delays: [0 0.75122] chips, C/No

−3 −2 −1 0 1 2 3

Bhuiyan & Lohan (2010)) and is sensitive to the noise dependent threshold choice. Firstly, it computes the noise variance, which is then used to compute an adaptive threshold. The peaks which are above the adaptive threshold are considered as competitive peaks. Among all the competitive peaks, TK selects the delay associated to that competitive peak which has

A Reduced Search Space Maximum Likelihood (RSSML) delay estimator is another good example of maximum likelihood based approach, which is capable of mitigating the multipath effects reasonably well at the expense of increased complexity. The RSSML, proposed by the Author in Bhuiyan et al. (2009) and then further enhanced in Bhuiyan & Lohan (2010), attempts to compensate the multipath error contribution by performing a nonlinear curve fit on the input correlation function which finds a perfect match from a set of ideal reference correlation functions with certain amplitude(s), phase(s) and delay(s) of the multipath signal. Conceptually, a conventional spread spectrum receiver does the same thing, but for only one signal (i.e., the LOS signal). With the presence of multipath signal, RSSML tries to separate the LOS component from the combined signal by estimating all the signal parameters in a maximum likelihood sense, which consequently achieves the best curve fit on the received input correlation function. As mentioned in Bhuiyan & Lohan (2010), it also incorporates a threshold-based peak detection method, which eventually reduces the code delay search space significantly. However, the downfall of RSSML is the memory requirement which it

In a multi-correlator based structure, the estimated LOS delay, theoretically, can be anywhere within the code delay window range of ±*τ<sup>W</sup>* chips, though in practice, it is quite likely to have a delay error around the previous delay estimate. The code delay window range essentially depends on the number of correlators (i.e., *M*) and the spacing between the correlators (i.e.,

Delay Error [Chips]

−0.5

0

0.5

Competitive Peaks

Fig. 6. Generation of competitive peaks for PT(Diff2) technique.

the closest delay difference from the previous delay estimate.

uses to store the reference correlation functions.

**5.5 Reduced Search Space Maximum Likelihood delay estimator**

Normalized Non−coherent Correlation Function

Normalized Diff2

1

$$
\tau\_W = \pm \frac{(M-1)}{2} \Delta \tag{4}
$$

For example, in Bhuiyan & Lohan (2010), 193 correlators were used with a correlator spacing of 0.0208 chips, resulting in a code delay window range of ±2 chips with respect to prompt correlator. Therefore, the LOS delay estimate can be anywhere within this ±2 chips window range. The ideal non-coherent reference correlation functions are generated for up to *L*max paths only for the middle delay index (i.e., ( *<sup>M</sup>*+<sup>1</sup> <sup>2</sup> )-th delay index; for *M* = 193, the middle delay index is 97). These ideal correlation functions for the middle delay index are generated off-line and saved in a look-up table in memory. In real-time, RSSML reads the correlation values from the look-up table, translates the ideal reference correlation functions at the middle delay index to the corresponding candidate delay index within the code delay window, and then computes the Minimum Mean Square Error (MMSE) for that specific delay candidate. Instead of considering all possible LOS delays within a predefined code delay window as delay candidates, the search space is first reduced to some competitive peaks which are generated based on the computed noise thresholds as explained in Bhuiyan & Lohan (2010). This will eventually reduce the processing time required to compute the MMSE (i.e., MMSE needs to be computed only for the reduced search space). An example is shown in Fig. 7, where RSSML estimates a best-fitted non-coherent correlation function at a cost of 3.6 · <sup>10</sup>−<sup>4</sup> MMSE in a two-path Rayleigh channel with path delays [0 0.35] chips, path powers [0 -2] dB at a C/N0 of 50 dB-Hz.

Fig. 7. Estimated and received non-coherent correlation functions in two path Rayleigh channel, path delay: [0 0.35] chips, path power: [0 -2] dB, C/N0: 50 dB-Hz.

**−2 −1.5 −1 −0.5 0 0.5 1 1.5 2**

Fig. 8. A non-coherent S-curve for CBOC(-) modulated single path static channel. Bhuiyan,

discrimination is applied to the TK output. The motivation for this combined approach comes from the fact that, when we apply TK operation to the non-coherent correlation function, it usually makes the main lobe of the non-coherent correlation function (after TK operation) much steeper than before. This eventually reduces the effect of multipath in case of TK based

The multipath performance of some of the discussed techniques are evaluated in different simulation models, i.e., the semi-analytical simulation in a static channel model and the Matlab-based simulation in a multipath fading channel model. Each simulation model is

The most typical way to evaluate the performance of a multipath mitigation technique is via Multipath Error Envelopes (MEE). Typically, two paths, either in-phase or out-of-phase, are assumed to be present, and the multipath errors are computed for multipath delays up to 1.2 chips at maximum, since the multipath errors become less significant after that. The upper multipath error envelope can be obtained when the paths are in-phase and the lower multipath error envelope when the paths are out-of-phase (i.e., 180◦ phase difference). In MEE analysis, several simplifying assumptions are usually made in order to distinguish the performance degradation caused by the multipath errors only. Such assumptions include zero Additive-White-Gaussian-Noise (AWGN), ideal infinite-length PRN codes, and zero residual

**Code Delay [chips]**

**Correct zero crossing**

**Nearest ambiguous**

**CBOC(−) signal, Single Path Static Channel,** Δ**EL=0.1 chips**

Multipath Mitigation Techniques for Satellite-Based Positioning Applications 419

**zero crossings for HRC Nearest ambiguous**

**zero crossings for nEML**

**nEML HRC**

**−0.8**

nEML (TK+nEML) as compared to nEML.

briefly described first before presenting the results.

Zhang & Lohan (2010)

**6. Performance analysis**

**6.1 Semi-analytical simulation**

**−0.6**

**−0.4**

**−0.2**

**Non−coherent S−curve**

**0**

**0.2**

**0.4**

**0.6**

**0.8**

#### **5.6 Slope-based multipath estimator**

A simple slope-based multipath mitigation technique, named as Slope-Based Multipath Estimator (SBME), in Bhuiyan, Lohan & Renfors (2010). Unlike the multipath mitigation techniques discussed above, SBME does not require a huge number of correlators, rather it only requires an additional correlator (as compared to a conventional DLL) at the late side of the correlation function. In fact, SBME is used in conjunction with a nEML tracking loop Bhuiyan, Lohan & Renfors (2010). It first derives a multipath estimation equation by utilizing the correlation shape of the ideal normalized correlation function, which is then used to compensate for the multipath bias of a nEML tracking loop. The derivation of the multipath estimation equation for BPSK modulated GPS L1 C/A signal can be found in Bhuiyan, Lohan & Renfors (2010). It is reported in Bhuiyan, Lohan & Renfors (2010) that SBME has superior multipath mitigation performance than nEML in closely spaced two path channel model.

## **5.7 C/N**0 **-based two-stage delay tracker**

A C/N0 based two-stage delay tracker is a combination of two individual tracking techniques, namely nEML and HRC (or MGD). The tracking has been divided into two stages based on the tracking duration and the received signal strength (i.e., C/N0). At the first stage of tracking (for about 0.1 seconds or so), the two-stage delay tracker always starts with a nEML tracking loop, since it begins to track the signal with a coarsely estimated code delay as obtained from the acquisition stage. And, at the second or final stage of tracking (i.e., when the DLL tracking error is around zero), the two-stage delay tracker switches its DLL discriminator from nEML to HRC (or MGD), since HRC (or MGD) has better multipath mitigation capability as compared to nEML. While doing so, it has to be ensured that the estimated C/N0 level meets a certain threshold set by the two-stage tracker. This is because of the fact that HRC (or MGD) involves one (or two in case of MGD) more discrimination than NEML, which makes its discriminator output much noisier than nEML. It has been empirically found that a C/N0 threshold of 35 dB-Hz can be a good choice, as mentioned in Bhuiyan, Zhang & Lohan (2010). Therefore, at this fine tracking stage, the two-stage delay tracker switches from nEML to HRC (or MGD) only when the estimated C/N0 meets the above criteria (i.e., C/N0 threshold is greater than 35 dB-Hz).

An example non-coherent S-curve is shown in Fig. 8 for CBOC(-) modulated signal in single path static channel Bhuiyan, Zhang & Lohan (2010). The nearest ambiguous zero crossings for HRC (around ±0.16 chips) is much closer as compared to that of nEML (around ±0.54 chips) in this particular case. This indicates the fact that the probability of locking to any of the side peaks is much higher for HRC than that of nEML, especially in the initial stage of tracking when the code delay may not necessarily be near the main peak of the correlation function. This is the main motivation to choose a nEML tracking at the initial stage for a specific time duration (for example, 0.1 seconds or so). This will eventually pull the delay tracking error around zero after the initial stage.

#### **5.8 TK operator combined with a nEML DLL**

A combined simplified approach with TK operator and a nEML DLL was implemented in Bhuiyan & Lohan (2010), in order to justify the feasibility of having a nEML discrimination after the TK operation on the non-coherent correlation function. In this combined approach, TK operator is first applied to the non-coherent correlation function, and then nEML

Fig. 8. A non-coherent S-curve for CBOC(-) modulated single path static channel. Bhuiyan, Zhang & Lohan (2010)

discrimination is applied to the TK output. The motivation for this combined approach comes from the fact that, when we apply TK operation to the non-coherent correlation function, it usually makes the main lobe of the non-coherent correlation function (after TK operation) much steeper than before. This eventually reduces the effect of multipath in case of TK based nEML (TK+nEML) as compared to nEML.

### **6. Performance analysis**

14 Will-be-set-by-IN-TECH

A simple slope-based multipath mitigation technique, named as Slope-Based Multipath Estimator (SBME), in Bhuiyan, Lohan & Renfors (2010). Unlike the multipath mitigation techniques discussed above, SBME does not require a huge number of correlators, rather it only requires an additional correlator (as compared to a conventional DLL) at the late side of the correlation function. In fact, SBME is used in conjunction with a nEML tracking loop Bhuiyan, Lohan & Renfors (2010). It first derives a multipath estimation equation by utilizing the correlation shape of the ideal normalized correlation function, which is then used to compensate for the multipath bias of a nEML tracking loop. The derivation of the multipath estimation equation for BPSK modulated GPS L1 C/A signal can be found in Bhuiyan, Lohan & Renfors (2010). It is reported in Bhuiyan, Lohan & Renfors (2010) that SBME has superior multipath mitigation performance than nEML in closely spaced two path channel model.

A C/N0 based two-stage delay tracker is a combination of two individual tracking techniques, namely nEML and HRC (or MGD). The tracking has been divided into two stages based on the tracking duration and the received signal strength (i.e., C/N0). At the first stage of tracking (for about 0.1 seconds or so), the two-stage delay tracker always starts with a nEML tracking loop, since it begins to track the signal with a coarsely estimated code delay as obtained from the acquisition stage. And, at the second or final stage of tracking (i.e., when the DLL tracking error is around zero), the two-stage delay tracker switches its DLL discriminator from nEML to HRC (or MGD), since HRC (or MGD) has better multipath mitigation capability as compared to nEML. While doing so, it has to be ensured that the estimated C/N0 level meets a certain threshold set by the two-stage tracker. This is because of the fact that HRC (or MGD) involves one (or two in case of MGD) more discrimination than NEML, which makes its discriminator output much noisier than nEML. It has been empirically found that a C/N0 threshold of 35 dB-Hz can be a good choice, as mentioned in Bhuiyan, Zhang & Lohan (2010). Therefore, at this fine tracking stage, the two-stage delay tracker switches from nEML to HRC (or MGD) only when the estimated C/N0 meets the above criteria (i.e., C/N0 threshold is

An example non-coherent S-curve is shown in Fig. 8 for CBOC(-) modulated signal in single path static channel Bhuiyan, Zhang & Lohan (2010). The nearest ambiguous zero crossings for HRC (around ±0.16 chips) is much closer as compared to that of nEML (around ±0.54 chips) in this particular case. This indicates the fact that the probability of locking to any of the side peaks is much higher for HRC than that of nEML, especially in the initial stage of tracking when the code delay may not necessarily be near the main peak of the correlation function. This is the main motivation to choose a nEML tracking at the initial stage for a specific time duration (for example, 0.1 seconds or so). This will eventually pull the delay tracking error

A combined simplified approach with TK operator and a nEML DLL was implemented in Bhuiyan & Lohan (2010), in order to justify the feasibility of having a nEML discrimination after the TK operation on the non-coherent correlation function. In this combined approach, TK operator is first applied to the non-coherent correlation function, and then nEML

**5.6 Slope-based multipath estimator**

**5.7 C/N**0 **-based two-stage delay tracker**

greater than 35 dB-Hz).

around zero after the initial stage.

**5.8 TK operator combined with a nEML DLL**

The multipath performance of some of the discussed techniques are evaluated in different simulation models, i.e., the semi-analytical simulation in a static channel model and the Matlab-based simulation in a multipath fading channel model. Each simulation model is briefly described first before presenting the results.

#### **6.1 Semi-analytical simulation**

The most typical way to evaluate the performance of a multipath mitigation technique is via Multipath Error Envelopes (MEE). Typically, two paths, either in-phase or out-of-phase, are assumed to be present, and the multipath errors are computed for multipath delays up to 1.2 chips at maximum, since the multipath errors become less significant after that. The upper multipath error envelope can be obtained when the paths are in-phase and the lower multipath error envelope when the paths are out-of-phase (i.e., 180◦ phase difference). In MEE analysis, several simplifying assumptions are usually made in order to distinguish the performance degradation caused by the multipath errors only. Such assumptions include zero Additive-White-Gaussian-Noise (AWGN), ideal infinite-length PRN codes, and zero residual

where AME(*xct*) is the mean of the absolute multipath error for the successive path delay *xct*. Now, the running average error for each particular delay in the range [0;1.167] chips can be

Multipath Mitigation Techniques for Satellite-Based Positioning Applications 421

*i* ∑ *i*=1

where *i* is the successive path delay index, and RAE(*xct*) is the RAE for the successive path

The RAE curves for CBOC(-) modulated Galileo E1C signal (i.e., pilot channel) is shown in Fig. 9. It is obvious from Fig. 9 that the RSSML and PT(Diff2) show the best performance in terms of RAE as compared to other techniques in this noise-free two to five paths static channel model. Among other techniques, TK+nEML showed very good performance followed by SBME, HRC and nEML. The SBME coefficient and the late slope at very late spacing of 0.0833 chips were determined according to Bhuiyan, Lohan & Renfors (2010) for a 24.552 MHz front-end bandwidth (double-sided). For the above configuration, the SBME coefficient is

It is worth to mention here that the RAE analysis is quite theoretical from two perspectives: firstly, the delay estimation is a one-shot estimate, and does not really include any tracking loop in the process; and secondly, the analysis is usually carried out with an ideal noise free assumption. These facts probably explain the reason why an algorithm which performs

CBOC(−) signal, 2 to 5 paths, BW 24.552 MHz

nEML HRC TK+nEML PT(Diff2) SBME RSSML

0 0.2 0.4 0.6 0.8 1 1.2 1.4

Multipath delay [chips]

very well with respect to RAE may not necessarily provide the same performance in a more realistic closed loop fading channel model, especially in the presence of more than two channel paths. However, MEE or RAE analysis has been widely used by the research community as an important tool for analyzing the multipath performance due to simpler implementation, and also due to the fact that it is hard to isolate multipath from other GNSS error sources in

Fig. 9. Running average error curves for CBOC(-) modulated Galileo E1C signal.

AME(*xct*)

*<sup>i</sup>* ,

(8)

RAE(*xct*) =

computed as follows:

0.007 and the late slope is −4.5.

0 0.5 1 1.5 2 2.5 3 3.5 4

Code tracking error [meters]

delay *xct*.

real life.

Doppler. Under these assumptions, the correlation R*rx*(*τ*) between the reference code of the modulation type MOD (for example, BPSK or CBOC(-)) and the received MOD-modulated signal via an *L*-path channel can be written as:

$$\mathcal{R}\_{\rm rx}(\tau) = \sum\_{l=1}^{L} \mathfrak{a}\_{l} e^{j\theta\_{l}} \mathcal{R}\_{\rm MCD}(\tau - \tau\_{l}) \tag{5}$$

where *αl*, *θl*, *τ<sup>l</sup>* are the amplitude, phase, and delay, respectively, of the *l*-th path; and RMOD(*τ*) is the auto-correlation function of a signal with the modulation type MOD. The analytical expressions for MEEs become complicated in the presence of more than two paths due to the complexity of channel interactions. Therefore, an alternative Monte-Carlo simulation-based approach is presented here, in accordance with Bhuiyan & Lohan (2010), for multipath error analysis in more than one path scenarios (i.e, for *L* ≥ 2). First, a sufficient number of random realizations, *N*random are generated (i.e., in the simulation, we choose *N*random equals to 1000), and then we look at the absolute mean error for each path delay over *N*random points. The objective is to analyze the multipath performance of some of the proposed advanced techniques along with some conventional DLLs in the presence of more than two channel paths, which may occur in urban or indoor scenarios.

The following assumptions are made while running the simulation for generating the RAE curves Hein et al. (2006). According to Hein et al. (2006), RAE is computed from the area enclosed within the multipath error and averaged over the range of the multipath delays from zero to the plotted delay values. In the simulation, the channel follows a decaying Power Delay Profile (PDP), which can be expressed by the equation:

$$\mathfrak{a}\_{l} = \mathfrak{a}\_{l} \exp^{-\mu(\tau\_{l} - \tau\_{1})},\tag{6}$$

where (*τ<sup>l</sup>* − *τ*1) �= 0 for *l* > 1, *μ* is the PDP coefficient (assumed to be uniformly distributed in the interval [0.05; 0.2], when the path delays are expressed in samples). The channel path phases *θ<sup>l</sup>* are uniformly distributed in the interval [0; 2*π*] and the number of channel paths *L* is uniformly distributed between 2 and *Lmax*, where *Lmax* is set to 5 in the simulation. A constant successive path spacing *xct* is chosen in the range [0; 1.167] chips with a step of 0.0417 chips (which defines the multipath delay axis in the RAE curves). It is worth to mention here that the number of paths is reduced to only one LOS path when *xct* = 0. The successive path delays can be found using the formula *τ<sup>l</sup>* = *lxct* in chips. Therefore, for each channel realization (which is a combination of amplitudes �*α* = *α*1,..., *αL*, phases �*θ* = *θ*1,..., *θL*, fixed path spacings, and the number of channel paths L), a certain LOS delay is estimated *<sup>τ</sup>*<sup>1</sup>(�*α*,�*θ*, *<sup>L</sup>*) from the zero crossing of the discriminator function (i.e., *<sup>D</sup>*(*τ*) = 0), when searched in the linear range of *D*(*τ*) in case of conventional DLLs, or directly from the auto-correlation function in case of advanced multi-correlator based techniques. The estimation error due to multipath is *<sup>τ</sup>*ˆ1(�*α*,�*θ*, *<sup>L</sup>*) <sup>−</sup> *<sup>τ</sup>*1, where *<sup>τ</sup>*<sup>1</sup> is the true LOS path delay. The RAE curves are generated in accordance with Hein et al. (2006). RAE is actually computed from the area enclosed within the multipath error and averaged over the range of the multipath delays from zero to the plotted delay values. Therefore, in order to generate the RAE curves, the Absolute Mean Error (AME) is computed for all *N*random random points via eqn. 7:

$$\text{AME}(\mathbf{x}\_{ct}) = \text{mean}(\left| \hat{\pi}\_1(\vec{a}, \vec{\theta}, L) - \pi\_1 \right|), \tag{7}$$

16 Will-be-set-by-IN-TECH

Doppler. Under these assumptions, the correlation R*rx*(*τ*) between the reference code of the modulation type MOD (for example, BPSK or CBOC(-)) and the received MOD-modulated

where *αl*, *θl*, *τ<sup>l</sup>* are the amplitude, phase, and delay, respectively, of the *l*-th path; and RMOD(*τ*) is the auto-correlation function of a signal with the modulation type MOD. The analytical expressions for MEEs become complicated in the presence of more than two paths due to the complexity of channel interactions. Therefore, an alternative Monte-Carlo simulation-based approach is presented here, in accordance with Bhuiyan & Lohan (2010), for multipath error analysis in more than one path scenarios (i.e, for *L* ≥ 2). First, a sufficient number of random realizations, *N*random are generated (i.e., in the simulation, we choose *N*random equals to 1000), and then we look at the absolute mean error for each path delay over *N*random points. The objective is to analyze the multipath performance of some of the proposed advanced techniques along with some conventional DLLs in the presence of more than two channel

The following assumptions are made while running the simulation for generating the RAE curves Hein et al. (2006). According to Hein et al. (2006), RAE is computed from the area enclosed within the multipath error and averaged over the range of the multipath delays from zero to the plotted delay values. In the simulation, the channel follows a decaying Power

*<sup>α</sup><sup>l</sup>* = *<sup>α</sup><sup>l</sup>* exp−*μ*(*τl*−*τ*1)

where (*τ<sup>l</sup>* − *τ*1) �= 0 for *l* > 1, *μ* is the PDP coefficient (assumed to be uniformly distributed in the interval [0.05; 0.2], when the path delays are expressed in samples). The channel path phases *θ<sup>l</sup>* are uniformly distributed in the interval [0; 2*π*] and the number of channel paths *L* is uniformly distributed between 2 and *Lmax*, where *Lmax* is set to 5 in the simulation. A constant successive path spacing *xct* is chosen in the range [0; 1.167] chips with a step of 0.0417 chips (which defines the multipath delay axis in the RAE curves). It is worth to mention here that the number of paths is reduced to only one LOS path when *xct* = 0. The successive path delays can be found using the formula *τ<sup>l</sup>* = *lxct* in chips. Therefore, for each channel realization (which is a combination of amplitudes �*α* = *α*1,..., *αL*, phases �*θ* = *θ*1,..., *θL*, fixed path spacings, and the number of channel paths L), a certain LOS delay is estimated *<sup>τ</sup>*<sup>1</sup>(�*α*,�*θ*, *<sup>L</sup>*) from the zero crossing of the discriminator function (i.e., *<sup>D</sup>*(*τ*) = 0), when searched in the linear range of *D*(*τ*) in case of conventional DLLs, or directly from the auto-correlation function in case of advanced multi-correlator based techniques. The estimation error due to multipath is *<sup>τ</sup>*ˆ1(�*α*,�*θ*, *<sup>L</sup>*) <sup>−</sup> *<sup>τ</sup>*1, where *<sup>τ</sup>*<sup>1</sup> is the true LOS path delay. The RAE curves are generated in accordance with Hein et al. (2006). RAE is actually computed from the area enclosed within the multipath error and averaged over the range of the multipath delays from zero to the plotted delay values. Therefore, in order to generate the RAE curves, the Absolute Mean

*<sup>α</sup>le<sup>j</sup>θ<sup>l</sup>* <sup>R</sup>MOD(*<sup>τ</sup>* <sup>−</sup> *<sup>τ</sup>l*) (5)

, (6)

), (7)

*L* ∑ *l*=1

signal via an *L*-path channel can be written as:

paths, which may occur in urban or indoor scenarios.

Delay Profile (PDP), which can be expressed by the equation:

Error (AME) is computed for all *N*random random points via eqn. 7:

AME(*xct*) = mean(

 

*<sup>τ</sup>*<sup>1</sup>(�*α*,�*θ*, *<sup>L</sup>*) <sup>−</sup> *<sup>τ</sup>*<sup>1</sup>

 

R*rx*(*τ*) =

where AME(*xct*) is the mean of the absolute multipath error for the successive path delay *xct*. Now, the running average error for each particular delay in the range [0;1.167] chips can be computed as follows:

$$\text{RAE}(\mathbf{x}\_{ct}) = \frac{\sum\_{i=1}^{l} \text{AME}(\mathbf{x}\_{ct})}{i},\tag{8}$$

where *i* is the successive path delay index, and RAE(*xct*) is the RAE for the successive path delay *xct*.

The RAE curves for CBOC(-) modulated Galileo E1C signal (i.e., pilot channel) is shown in Fig. 9. It is obvious from Fig. 9 that the RSSML and PT(Diff2) show the best performance in terms of RAE as compared to other techniques in this noise-free two to five paths static channel model. Among other techniques, TK+nEML showed very good performance followed by SBME, HRC and nEML. The SBME coefficient and the late slope at very late spacing of 0.0833 chips were determined according to Bhuiyan, Lohan & Renfors (2010) for a 24.552 MHz front-end bandwidth (double-sided). For the above configuration, the SBME coefficient is 0.007 and the late slope is −4.5.

It is worth to mention here that the RAE analysis is quite theoretical from two perspectives: firstly, the delay estimation is a one-shot estimate, and does not really include any tracking loop in the process; and secondly, the analysis is usually carried out with an ideal noise free assumption. These facts probably explain the reason why an algorithm which performs

Fig. 9. Running average error curves for CBOC(-) modulated Galileo E1C signal.

very well with respect to RAE may not necessarily provide the same performance in a more realistic closed loop fading channel model, especially in the presence of more than two channel paths. However, MEE or RAE analysis has been widely used by the research community as an important tool for analyzing the multipath performance due to simpler implementation, and also due to the fact that it is hard to isolate multipath from other GNSS error sources in real life.

performance in this two to five paths closely spaced multipath channel. Among other techniques, PT(Diff2) and HRC have better performance only in good C/N0 (around 40 dB-Hz and onwards). It can also be observed that the proposed SBME and TK+nEML do not bring any advantage in the tracking performance as compared to nEML in this multipath fading channel model. Here also, the SBME coefficient and the late slope were set to 0.007 and −4.5,

Multipath Mitigation Techniques for Satellite-Based Positioning Applications 423

CBOC(−) signal, 2 to 5 paths, BW 24.552 MHz

nEML HRC TK+nEML PT(Diff2) SBME RSSML

30 35 40 45

[dB−Hz]

C/N0

Fig. 10. RMSE vs. C/N0 plot for CBOC(-) modulated Galileo E1C signal in two to five path

This chapter addressed the challenges encountered by a GNSS signal due to multipath propagation. A wide range of correlation-based multipath mitigation techniques were discussed and the performance of some of these techniques were evaluated in terms of running average error and root-mean-square error. Among the analyzed multipath mitigation techniques, RSSML, in general, achieved the best multipath mitigation performance in moderate-to-high C/N0 scenarios (for example, 30 dB-Hz and onwards). The other techniques, such as PT(Diff2) and HRC showed good multipath mitigation performance only in high C/N0 scenarios (for example, 40 dB-Hz and onwards). The other new technique SBME offered slightly better multipath mitigation performance to the well-known nEML DLL at the cost of an additional correlator. However, as the GNSS research area is fast evolving with many potential applications, it remains a challenging topic for future research to investigate the feasibility of these multipath mitigation techniques with the multitude of signal modulations, spreading codes, and spectrum placements that are (or are to be)

Baltersee, J., Fock, G. & Schulz-Rittich, P. (2001). Adaptive code-tracking receiver for

direct-sequence Code Division Multiple Access (CDMA) communications over

respectively.

0

Rayleigh fading channel.

**7. Conclusion**

proposed.

**8. References**

10

20

30

RMSE [meters]

40

50

60

70

#### **6.2 Matlab-based simulation**

Simulation has been carried out in closely spaced multipath environments for CBOC(-) modulated Galileo E1C (i.e., pilot channel) signal for a 24.552 MHz front-end bandwidth. The simulation profile is summarized in Table 1. Rayleigh fading channel model is used in the simulation, where the number of channel paths follows a uniform distribution between two and five. The successive path separation is random between 0.02 and 0.35 chips. The channel paths are assumed to obey a decaying PDP according to Eqn. 6, where (*τ<sup>l</sup>* − *τ*1) �= 0 for *l* > 1, and the PDP coefficient *μ* = 0.1. The received signal is sampled such that there are 48 samples per chip. The received signal duration is 800 milliseconds (ms) or 0.8 seconds for each particular C/N0 level. The tracking errors are computed after each *Nc* ∗ *Nnc* ms (in this case, *Nc* ∗ *Nnc* = 20 ms) interval. In the final statistics, the first 600 ms are ignored in order to remove the initial error bias that may come from the delay difference between the received signal and the locally generated reference code. Therefore, for the above configuration (i.e., code loop filter parameters and the first path delay of 0.2 chips), the left-over tracking errors after 600 ms are mostly due to the effect of multipath. The simulation has been carried out for 100 random realizations, which give a total of 10\*100=1000 statistical points, for each C/N0 level. The RMSE of the delay estimates are plotted in meters, by using the relationship:

$$\text{RMSE}\_{\text{m}} = \text{RMSE}\_{\text{chips}} c T\_{\text{c}} \tag{9}$$


where *c* is the speed of light, *Tc* is the chip duration, and RMSEchips is the RMSE in chips.

Table 1. Simulation profile description

RMSE vs. C/N0 plot for the given multipath channel profile is shown in Fig. 10. It can be seen from Fig. 10 that the proposed RSSML clearly achieves the best multipath mitigation

<sup>1</sup> Not all the correlators are used in all the tracking algorithms. For example, nEML only requires 3 correlators.

performance in this two to five paths closely spaced multipath channel. Among other techniques, PT(Diff2) and HRC have better performance only in good C/N0 (around 40 dB-Hz and onwards). It can also be observed that the proposed SBME and TK+nEML do not bring any advantage in the tracking performance as compared to nEML in this multipath fading channel model. Here also, the SBME coefficient and the late slope were set to 0.007 and −4.5, respectively.

Fig. 10. RMSE vs. C/N0 plot for CBOC(-) modulated Galileo E1C signal in two to five path Rayleigh fading channel.

## **7. Conclusion**

18 Will-be-set-by-IN-TECH

Simulation has been carried out in closely spaced multipath environments for CBOC(-) modulated Galileo E1C (i.e., pilot channel) signal for a 24.552 MHz front-end bandwidth. The simulation profile is summarized in Table 1. Rayleigh fading channel model is used in the simulation, where the number of channel paths follows a uniform distribution between two and five. The successive path separation is random between 0.02 and 0.35 chips. The channel paths are assumed to obey a decaying PDP according to Eqn. 6, where (*τ<sup>l</sup>* − *τ*1) �= 0 for *l* > 1, and the PDP coefficient *μ* = 0.1. The received signal is sampled such that there are 48 samples per chip. The received signal duration is 800 milliseconds (ms) or 0.8 seconds for each particular C/N0 level. The tracking errors are computed after each *Nc* ∗ *Nnc* ms (in this case, *Nc* ∗ *Nnc* = 20 ms) interval. In the final statistics, the first 600 ms are ignored in order to remove the initial error bias that may come from the delay difference between the received signal and the locally generated reference code. Therefore, for the above configuration (i.e., code loop filter parameters and the first path delay of 0.2 chips), the left-over tracking errors after 600 ms are mostly due to the effect of multipath. The simulation has been carried out for 100 random realizations, which give a total of 10\*100=1000 statistical points, for each C/N0 level. The RMSE of the delay estimates are plotted in meters, by using the relationship:

where *c* is the speed of light, *Tc* is the chip duration, and RMSEchips is the RMSE in chips.

Channel model Rayleigh fading channel

Path Power Decaying PDP with *μ* = 0.1

Path Phase Random between 0 and 2*π*

Filter Type Finite Impulse Response (FIR)

RMSE vs. C/N0 plot for the given multipath channel profile is shown in Fig. 10. It can be seen from Fig. 10 that the proposed RSSML clearly achieves the best multipath mitigation

<sup>1</sup> Not all the correlators are used in all the tracking algorithms. For example, nEML only requires 3

Path Spacing Random between 0.02 and 0.35 chips

Parameter Value

Samples per Chip, *Ns* 48

Filter Order 6 Coherent Integration, *Nc* 20 ms Non-coherent Integration, *Nnc* 1 block Initial Delay Error ± 0.1 chips First Path Delay 0.2 chips Code Tracking Loop Bandwidth 2 Hz Code Tracking Loop Order 1st order

Table 1. Simulation profile description

correlators.

Number of paths between 2 to 5

E-L Spacing, Δ*EL* 0.0833 chips Number of Correlators, *M* 193<sup>1</sup>

Double-sided Bandwidth, BW 24.552 MHz

RMSE*<sup>m</sup>* = RMSEchips*cTc* (9)

**6.2 Matlab-based simulation**

This chapter addressed the challenges encountered by a GNSS signal due to multipath propagation. A wide range of correlation-based multipath mitigation techniques were discussed and the performance of some of these techniques were evaluated in terms of running average error and root-mean-square error. Among the analyzed multipath mitigation techniques, RSSML, in general, achieved the best multipath mitigation performance in moderate-to-high C/N0 scenarios (for example, 30 dB-Hz and onwards). The other techniques, such as PT(Diff2) and HRC showed good multipath mitigation performance only in high C/N0 scenarios (for example, 40 dB-Hz and onwards). The other new technique SBME offered slightly better multipath mitigation performance to the well-known nEML DLL at the cost of an additional correlator. However, as the GNSS research area is fast evolving with many potential applications, it remains a challenging topic for future research to investigate the feasibility of these multipath mitigation techniques with the multitude of signal modulations, spreading codes, and spectrum placements that are (or are to be) proposed.

#### **8. References**

Baltersee, J., Fock, G. & Schulz-Rittich, P. (2001). Adaptive code-tracking receiver for direct-sequence Code Division Multiple Access (CDMA) communications over

Fenton, P. C. & Dierendonck, A. J. V. (1996). Pseudorandom noise ranging receiver which

Multipath Mitigation Techniques for Satellite-Based Positioning Applications 425

spacing between early and late correlators, NovAtel Patent, US 5 495 499. Fenton, P. C. & Jones, J. (2005). The theory and performance of NovAtel Inc.'s Vision

Fine, P. & Wilson, W. (1999). Tracking algorithms for GPS offset carrier signals, *Proc. of ION*

Fock, G., Baltersee, J., Schulz-Rittich, P. & Meyr, H. (2001). Channel tracking for RAKE

Garin, L. & Rousseau, J. M. (1997). Enhanced Strobe Correlator multipath rejection for code

Granados, G. S. & Rubio, J. A. F. (2000). Time-delay estimation of the line-of-sight signal

Granados, G. S., Rubio, J. A. F. & Prades, C. F. (2005). ML estimator and hybrid beamformer for

Hamila, R. (2002). *Synchronization and multipath delay estimation algorithms for digital receivers*,

Hamila, R., Lohan, E. & Renfors, M. (2003). Subchip multipath delay estimation for downlink

Hein, G. W., Avila-Rodriguez, J.-A., Wallner, S., Pratt, A. R., Owen, J., Issler, J. L., Betz, J. W.,

Irsigler, M. & Eissfeller, B. (2003). Comparison of multipath mitigation techniques with

Jie, Z. (2010). *Delay tracker for Galileo CBOC modulated signals and their Simulink-based*

Jones, J., Fenton, P. & Smith, B. (2004). Theory and Performance of the Pulse Aperture

Kaplan, E. D. & Hegarty, C. J. (2006). *Understanding GPS: Principles and Applications*, second

Laxton, M. (1996). *Analysis and simulation of a new code tracking loop for GPS multipath mitigation*,

Lee, C., Yoo, S., Yoon, S. & Kim, S. Y. (2006). A novel multipath mitigation scheme based

Lohan, E. S. (2003). *Multipath delay estimators for fading channels with applications in CDMA receivers and mobile positioning*, PhD thesis, Tampere University of Technology.

on slope differential of correlator output for Galileo systems, *Proceedings of 8th*

*implementations*, Master's thesis, Tampere University of Technology.

Correlator. Tech. Rep., NovAtel, Calgary, Alberta, Canada.

*International Conference on Advanced Communication Technology*.

Master's thesis, Air Force Institute of Technology.

Correlator, *Proc. of ION GNSS*, Long Beach, USA, pp. 2178–2186.

and carrier, *Proc. of ION GPS*, Kansas City, USA, pp. 559–568. Gleason, S. & Egziabher, D. G. (2009). *GNSS Applications and Methods*, Artech House.

*NTM*, San Diego, USA.

*Comm.* 19(12): 2420–2431.

*(EUSIPCO)*, Tampere, Finland.

*Processing* 53(3): 1194–1208.

*DOI:10.1155/2008/785695* .

edn, Artech House.

7(1): 125–128.

PhD thesis, Tampere University of Technology.

compensates for multipath distortion by dynamically adjusting the time delay

receivers in closely spaced multipath environments, *IEEE Journal on Sel. Areas in*

in a multipath environment, *Proc. of the 10th European Signal Processing Conference*

multipath and interference mitigation in GNSS receivers, *IEEE Transactions on Signal*

WCDMA systems based on Teager-Kaiser operator, *IEEE Communications Letters*

Hegarty, C. J., Lenahan, L. S., Rushanan, J. J., Kraay, A. L. & Stansell, T. A. (2006). MBOC: The new optimized spreading modulation recommended for GALILEO L1 OS and GPS L1C, *Proc. of Position, Location and Navigation Symposium*, pp. 883–892. Hurskainen, H., Lohan, E. S., Hu, X., Raasakka, J. & Nurmi, J. (2008). Multiple

gate delay tracking structures for GNSS signals and their evaluation with Simulink, SystemC, and VHDL, *International Journal of Navigation and Observation,*

consideration of future signal structures, *Proc. of ION GNSS*, OR, USA, pp. 2584–2592.

multipath fading channels and method for signal processing in a RAKE receiver, US Patent Application Publication, US 2001/0014114 A1 (Lucent Technologies).


20 Will-be-set-by-IN-TECH

Patent Application Publication, US 2001/0014114 A1 (Lucent Technologies). Bello, P. A. & Fante, R. L. (2005). Code tracking performance for novel unambiguous M-code

Betz, J. W. & Kolodziejski, K. R. (2000). Extended theory of early-late code tracking

Bhuiyan, M. Z. H. (2006). *Analyzing code tracking algorithms for Galileo Open Service signal*,

Bhuiyan, M. Z. H. (2011). *Analysis of Multipath Mitigation Techniques for Satellite-based Positioning Applications*, PhD thesis, Tampere University of Technology. Bhuiyan, M. Z. H. & Lohan, E. S. (2010). Advanced multipath mitigation techniques

Bhuiyan, M. Z. H., Lohan, E. S. & Renfors, M. (2008). Code tracking algorithms for mitigating

Bhuiyan, M. Z. H., Lohan, E. S. & Renfors, M. (2009). A reduced search space maximum

positioning, *Proc. of the 13th International Association of Institute of Navigation*. Bhuiyan, M. Z. H., Lohan, E. S. & Renfors, M. (2010). A slope-based multipath estimation

*Special Session on Circuits, Systems and Algorithms for Next Generation GNSS*. Bhuiyan, M. Z. H., Zhang, J. & Lohan, E. S. (2010). Enhanced delay tracking performance of

Bischoff, R., Häb-Umbach, R., Schulz, W. & Heinrichs, G. (2002). Employment of a multipath

Braasch, M. S. (2001). Performance comparison of multipath mitigating receiver architectures, *Proc. of IEEE Aerospace Conference*, Big Sky, Mont, USA, pp. 1309–1315. Chen, K. C. & Davisson, L. D. (1994). Analysis of a SCCL as a PN code tracking loop, *IEEE*

Dierendonck, A. J. V. & Braasch, M. S. (1997). Evaluation of GNSS receiver correlation

Dierendonck, A. J. V., Fenton, P. C. & T., F. (1992). Theory and performance of narrow correlator spacing in a GPS receiver, *Journal of The Institute of Navigation* 39(3). Fante, R. (2003). Unambiguous tracker for GPS Binary-Offset Carrier signals, *Proc. of ION*

Fante, R. (2004). Unambiguous First-Order Tracking Loop M-Code, MITRE Technical Report.

Fenton, P. C. (1995). Pseudorandom noise ranging receiver which compensates for multipath

distortion by making use of multiple correlator time delay spacing, NovAtel Patent,

*Conference on Global Navigation Satellite Systems, ENC GNSS 2010*.

*Transactions on Communications* 42(11): 2942–2946.

time discriminators, *Proc. of ION NTM*, San Diego, USA.

Master's thesis, Tampere University of Technology.

*on Advances in Signal Processing, DOI: 10.1155/2008/863629* .

*Observation, DOI:10.1155/2010/412393* .

47(3): 211–226.

Vol. 4, pp. 1844–1848.

USA, pp. 207–215.

MTR 94B0000040.

US 5 414 729.

*NTM*, Albuquerque, USA.

multipath fading channels and method for signal processing in a RAKE receiver, US

for a bandlimited GPS receiver, *Navigation: Journal of the Institute of Navigation*

for satellite-based positioning applications, *International Journal of Navigation and*

multipath effects in fading channels for satellite-based positioning, *EURASIP Journal*

likelihood delay estimator for mitigating multipath effects in satellite-based

technique for mitigating short-delay multipath in GNSS receivers, *Proc. of ISCAS in*

a C/N0-based two-stage tracker for GNSS receivers, *Proc. of the European Navigation*

receiver structure in a combined Galileo/UMTS receiver, *Proc. of IEEE VTC Spring*,

processing techniques for multipath and noise mitigation, *Proc. of ION NTM*, CA,


22 Will-be-set-by-IN-TECH

426 Global Navigation Satellite Systems – Signal, Theory and Applications

Lohan, E. S., Lakhzouri, A. & Renfors, M. (2006). Feedforward delay estimators in adverse

McGraw, G. A. & Braasch, M. S. (1999). GNSS multipath mitigation using Gated and High

Nee, R. D. J. V. (1992). The multipath estimating delay lock loop, *Proc. of IEEE Second*

Nee, R. D. J. V., Sierveld, J., Fenton, P. C. & Townsend, B. R. (1994). The multipath estimating

Phelts, R. E. & Enge, P. (2000a). The multipath invariance approach for code multipath

Phelts, R. E. & Enge, P. (2000b). Multipath mitigation for narrowband receivers, *Proc. of the*

Simon, M. K., Omura, J. K., Scholtz, R. A. & Levitt, B. K. (1994). *Spread Spectrum Communication*

Sleewaegen, J. M. & Boon, F. (2001). Mitigating short-delay multipath: a promising new

Townsend, B. R. & Fenton, P. C. (1994). A practical approach to the reduction of pseudorange multipath errors in a L1 GPS receiver, *Proc. of ION GPS*, UT, USA, pp. 143–148. Townsend, B. R., Nee, R. D. J. V., Fenton, P. C. & Dierendonck, A. J. V. (1995). Performance

Weill, L. R. (1997). A GPS multipath mitigation by means of correlator reference waveform

Weill, L. R. (2002). Multipath mitigation using modernized GPS signals: How good can it get?,

Weill, L. R. (2003). Multipath mitigation - how good can it get with new signals?, *GPS World*

evaluation of the multipath estimating delay lock loop, *Proc. of ION NTM*, Anaheim,

*Advances in Signal Processing, Article ID 50971* .

Japan, pp. 39–42.

USA.

16(6): 106–113.

*Division of the Insitute of Navigation*, San Diego, USA.

*Location and Navigation Symposium*, Vol. 1, pp. 246–251.

mitigation, *Proc. of ION GPS*, UT, USA, pp. 2376–2384.

technique, *Proc. of ION GPS*, UT, USA, pp. 204–213.

design, *Proc. of the of ION NTM*, CA, USA, pp. 197–206.

*Proc. of the of ION GPS*, CA, USA, pp. 493–505.

*Handbook*, revised edition edn, McGraw-Hill Inc, New York.

*IEEE PLANS 2000*, CA, USA, pp. 30–36.

multipath propagation for Galileo and modernized GPS signals, *EURASIP Journal on*

Resolution Correlator concepts, *Proc. of the National Technical Meeting of the Satellite*

*International Symposium on Spread Spectrum Techniques and Applications*, Yokohama,

delay lock loop: Approaching theoretical accuracy limits, *Proc. of IEEE Position*

## *Edited by Shuanggen Jin*

Global Navigation Satellite System (GNSS) plays a key role in high precision navigation, positioning, timing, and scientific questions related to precise positioning. This is a highly precise, continuous, all-weather, and real-time technique. The book is devoted to presenting recent results and developments in GNSS theory, system, signal, receiver, method, and errors sources, such as multipath effects and atmospheric delays. Furthermore, varied GNSS applications are demonstrated and evaluated in hybrid positioning, multi-sensor integration, height system, Network Real Time Kinematic (NRTK), wheeled robots, and status and engineering surveying. This book provides a good reference for GNSS designers, engineers, and scientists, as well as the user market.

Photo by liulolo / iStock

Global Navigation Satellite Systems - Signal, Theory and Applications

Global Navigation

Satellite Systems

Signal, Theory and Applications

*Edited by Shuanggen Jin*