GPS Filtering to Minimize Ionosphere Divergence Error for Safe Aircraft Landing

By

Shiladitya Sen B.Tech., Manufacturing Sc. & Engineering Indian Institute of Technology, Kharagpur, 2006

SUBMITTED TO THE DEPARTMENT OF MECHANICAL ENGINEERING IN PARTIAL FULLFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE IN MECHANICAL ENGINEERING AT TUFTS UNIVERSITY

May 2008 © Tufts University. All rights reserved.

Signature of Author ………………………………………………………………………... Department of Mechanical Engineering May 19, 2007

Certified By………………………………………………………………………………… Jason H. Rife Assistant Professor of Mechanical Engineering Thesis Supervisor 1

2

GPS Filtering to Minimize Ionosphere Divergence Error for Safe Aircraft Landing by Shiladitya Sen Submitted to the Department of Mechanical Engineering on May 1, 2008, in partial fulfillment of the requirements for the Degree of Master of Science in Mechanical Engineering

Abstract Aircraft use navigation guides to help pilots land especially in low visibility conditions. The commonly used Instrument Landing System is very effectively yet expensive and inflexible. The Global Positioning System (GPS) can be used to design potentially cheaper and more versatile navigation systems. The Local Area Augmented System (LAAS) is one of the proposed systems for navigation. However it is yet to be certified by the Federal Aviation Administration (FAA) owing to its susceptibility to certain rare events which can potentially cause a system failure in LAAS. LAAS implements the idea of using differential corrections from a nearby reference station to increase the accuracy of the position estimate of the user. One of the biggest threats to LAAS, which undermines the concept of differential corrections, is an ionospheric storm caused by high solar activity. The ionosphere is a spherical shell of ionized gases over earth and affects radiowave communication passing through it. Although ionospheric storms are rare events, under certain worst cases they can be undetectable and might cause a missed landing. Therefore to ensure the safety of the usage of LAAS, this problem is addressed in this research. Errors due to the ionospheric storm are magnified by the standard real-time smoothing filter in LAAS, the Hatch Filter. The Hatch Filter leverages the two components of GPS signals, namely the code and carrier measurement, and provides excellent attenuation of noise under nominal conditions. However under the influence of the ionospheric storm, the Hatch Filter suffers a delay in its position estimate and hence produces a systematic bias. This thesis designs two alternative real-time smoothing filters with the motivation of increasing insensitivity towards ionospheric storms without sacrificing the “nominal day” noise attenuation. The first design uses mathematical rigor to optimize a general linear digital filter for the least maximum transient error under the influence of an ionospheric storm. The optimized linear filter is found to reduce the maximum transient error by a maximum of 20% during the offline simulation test conducted on GPS data. The second filter is designed with the motivation of reducing steady state errors caused by the ionospheric storm. This filter is nonlinear and is able to overcome the limitations of linear filters (which suffer degraded transient error while reducing steady state errors). In the offline experiments the nonlinear filter attained zero steady state error in each case and also reduced the maximum transient error as compared to the Hatch Filter. Thesis Supervisor: Jason H. Rife Title: Assistant Professor of Mechanical Engineering, Tufts University 3

4

Acknowledgements Delving into the intricacies of GPS navigation has been one of the most wonderfully rewarding experiences of my life. It has been challenging at all times to deal with a project of such far-reaching consequences but the intellectual stimulation it has provided me has been worth the journey. And the ride would not be what it has been if not for all the support that I have received from some truly remarkable people. First and foremost, I would like to acknowledge the unwavering support and insightful guidance of Professor Jason Rife. Without his resourcefulness and creative inputs, the project would have faltered many a time. His enthusiasm for the project never flagged and he always encouraged me to come up with new ideas. It has been a pleasure to work with him. Through my interactions with him, I have grown as a researcher and as a person. I would also like to thank Professor Anil Saigal who helped me find projects matching my interests. His suggestions regarding various department matters as well as day to day matters helped me settle very well at Tufts University. I owe a great deal of gratitude to my undergraduate studies advisor Professor A. K. Chattopadhyay for teaching me the ABCs of research. Under him I gained a lot of insight into thinking scientifically. He also helped me in finding my field of interest within mechanical engineering. My friends, from Kolkata and Kharagpur, deserve a special mention. My conversations and emails with each and every one of them have let me lead a normal existence. I cannot thank them enough for that. Two of them deserve a special mention in this regard. Titli, especially provided me with a constant sounding board. Her enthusiasm in listening to the maths behind GPS must be commended. A special token of thanks needs to go out to Ishan during the period when he was at MIT. He has given me a helping hand, right from groceries to MEMS. Our conversations must have made Cingular rethink their phone plans. Finally, I do not know where and how to start thanking my parents. There are no words to express my feelings of gratitude, and really anything I will put down does not do justice to their efforts. They have always – with a capital A – believed in me and my pursuits. Never have they had anything but words of encouragement and support. Baba and Ma, I would like you to know that one cannot have better mentors, friends and guides than you have been.

5

6

Table of Contents ABSTRACT...............................................................................................................3 ACKNOWLEDGEMENTS.............................................................................................5 TABLE OF CONTENTS................................................................................................7 LIST OF FIGURES........................................................................................................9 LIST OF TABLES.......................................................................................................11 CHAPTER 1: INTRODUCTION...................................................................................12 CHAPTER 2: GPS ERRORS........................................................................................21 CHAPTER 3: DATA SETS...........................................................................................27 CHAPTER 4: OPTIMUM LINEAR FILTER DESIGN......................................................33 CHAPTER 5: NONLINEAR DIVERGENCE ELIMINATION (NLDE) FILTER DESIGN........55 CHAPTER 6: CONCLUSION.......................................................................................72 REFERENCES............................................................................................................74 APPENDIX................................................................................................................76

7

8

List of figures Fig. 1.1:

Trilateration

13

Fig. 1.2:

GPS satellite constellation

15

Fig. 1.3:

Divergence bias error in Hatch filter

18

Fig. 2.1:

(a) The Klobuchar ionospheric model for nominal ionospheric delay (b)Nominal Ionospheric delay on 17 January 2008 at CORS station ZOA1

23 23

Fig. 2.2:

(a) Nominal ionospheric delays (b) Nominal ionospheric delay with an anomalous event due to ionospheric Storm

24

Fig. 2.3:

Ionosphere threat model

25

Fig. 3.1:

Sample Rinex file

28

Fig. 3.2:

Code measurement noise on simultaneously all visible satellites

30

Fig. 3.3:

GPS 18 USB

31

Fig. 3.4:

Pseudorange noise from experimental data of Satellite PRN 32

32

Fig. 4.1:

(a) Nominal GPS errors (sensor noise + nominal ionospheric delay) (b) Anomalous GPS errors (sensor noise + anomalous ionospheric delay + divergence bias error)

34

Fig. 4.2:

Coefficient of Hatch Filter impulse function

36

Fig. 4.3:

Code bias and CUC bias due to anomalous ionospheric storm

37

Fig. 4.4:

CUC divergence error Ecuc,i|k

38

Fig. 4.5:

Tracking error for OL filters for different lengths

42

Fig. 4.6:

Coefficients of the MTOL filter and the Hatch Filter

46

Fig. 4.7:

Maximum Divergence bias Error as a Function of White-Gaussian Noise Rejection for Hatch and MTOL Filter

47

Maximum Ionosphere Anomaly Tracking Error as a Function of Actual (Correlated) Noise Rejection for Hatch and MTOL Filters

48

Power Spectrum of the transfer functions of the two filters

50

Fig. 4.8:

Fig. 4.9:

9

24

34

Fig. 4.10:

Power Spectra for GPS Pseudorange Input Data, Normalized to a Peak Value of One

51

Fig. 4.11:

Ratio of Anomalous Ionosphere Bias for MTOL and Hatch Filters. Results are plotted as a Function of Output Noise Attenuation (for Filters of Different Lengths) Applied to Eight Satellites in View 52

Fig. 4.12:

Ratio of Output Noise Standard Deviations for MTOL and Hatch Filters (with Same Maximum Ionosphere Bias) as a Function of Hatch Time Constant. Data Plotted for All 8 Satellites in View

52

Fig. 4.13:

Comparison of transients of Hatch Filter and ZSEL filter of similar noise output

54

Fig. 5.1:

Comparison of Hatch Filter, MTOL filter, ZSEL filter

56

Fig. 5.2:

Example of the curve fit on ionospheric delay data

57

Fig. 5.3:

(a) Instantaneous ionopheric gradient change modelled by 1 slope and 2 slope Method 58 (b) Gradual ionopheric gradient change modelled by 1 slope and 2 slope Method

58

Fig. 5.4:

Block diagram for NLDE filter

59

Fig. 5.5:

Block diagram of the nonlinear divergence estimator in the NLDE filter

59

Fig. 5.6:

Response comparison of NLDE filter and Hatch Filter of same length to an ionospheric event

62

Fig. 5.7:

Noise characteristic comparison of NLDE filter and Hatch Filter of same length

62

Fig. 5.8:

(a) Power spectrum of noise from code measurement input (b) Power spectrum of noise from Hatch Filter (c) Power spectrum of noise from NLDE filter

64 64 64

Fig. 5.9:

Response of NLDE filter to pseudorange with ionospheric event

65

Fig. 5.10:

Dependence of Noise output ratio and Maximum transient error on Hatch Filter length M and minimum length L

67

Fig. 5.11:

Comparison of NLDE and Hatch Filter

69

Fig. 5.12:

Comparison of Hatch Filter, optimized linear filter and NLDE filter done on PRN17 data

70

10

List of tables Table 2.1:

GPS errors

22

Table 3.1:

USB packet structure

31

Table 5.1:

Analyzed parameter space

68

11

Chapter 1 Introduction Navigation has been important to man all throughout history. Being able to navigate with accuracy and precision is important for commercial vehicles to move smartly and safely from point A to point B. While for land transportation navigation is helpful to find the shortest and the quickest routes, for aircraft safety is the bigger issue. Navigation becomes extremely important for aircraft while landing in low visibility conditions. To safely land planes in those conditions the pilots needs to know his/her exact position, velocity and acceleration relative to the exact location of the runway at the airport. In other words the aircraft position with respect to the airport needs to be calculated in real time. Thus we require a real-time position sensor. The Instrument Landing System is presently used as a navigational aid to pilots for landing. The ILS is a ground based instrument approach system which uses a combination of radio signals and high intensity light arrays to guide aircraft to the runway [1]. The Local Area Augmentation System (LAAS) provides an alternative to ILS and has a potential to be more flexible and much less expensive as it leverages Global Positioning System (GPS). The challenge of certifying LAAS is to prove that the navigation sensor will not cause a missed landing or an accident. Although GPS, by itself, does not provide such precision in position estimation, using differentially corrected LAAS measurements was shown to provide reliability in most cases. However, as the confidence bounds required for LAAS is a probability of approximately one in a billion, it is necessary to consider even very rare fault modes to assess system safety. Unless LAAS proves reliability under all conditions, it is unlikely to replace ILS. The most severe system anomalies are caused by ionospheric storms [2]. Although these events may occur only once in an eleven year solar cycle, they can produce errors of tens of meters in magnitude. LAAS relies on a spatially uniform ionosphere assumption which works fine under nominal conditions but fails during these events. The challenge is to either monitor and detect these anomalies reliably or design the system to be insensitive to this occurrence. Hence this thesis attempt to design signal processing based solutions to the “ionospheric storm” problem faced in LAAS. The next five sections discuss, in more detail, the principles involved in the working of these GPS based systems, the signal structure, the ionospheric threat, previous work and the specific contribution of this thesis.

1.1 Satellite navigation systems: GPS, LAAS and WAAS GPS works on the principle that position can be determined by measuring the distances to three points at known locations. This principle is sometimes called “trilateration” [3, 4]. Figure 1.1 pictorially

12

describes the concept of trilateration. The distance between the unknown point P and the reference station S1 is R1. This distance is also called range of P from S1. Hence P is a point on the circle with center at S1 and radius R1. Similar circles can also be drawn for S2 and S3 with radius R2 and R3 respectively. The intersection of these three circles is P.

Figure 1.1: Trilateration

To estimate the distances R1, R2, R3 in figure 1.1, the system uses the time taken by electromagnetic waves to travel between the sources and the receiver locate at P. This requires the clock at both the ends to be synchronized. However since this is impossible to ensure on a global scale, it is safe to assume that the clock at the source’s end is lagging/leading the clock at the receiver end. The use of synced atomic clocks, which are very precise and stable, for the reference stations (satellites in the case of GPS) ensure that all the satellite have the same time reference. Hence looking back at trilateration, we have another unknown which is the time difference between the receiver clock at P and all the satellites. Thus the number of unknowns increase to four, which logically leads to the need of another reference station or satellite. This means the GPS really needs a minimum of four visible satellites for position calculation. Due to economic constraints, users generally use less expensive clocks with quartz oscillators which are much less stable and do not stay synced. Thus the distance measurements aren’t really true ranges but pseudoranges, subject to a timing error that is consistent on all GPS channels. GPS is widely used in aviation for non-precision flight [5]. For precision applications like aircraft landing, however, the accuracy of stand-alone GPS is unsuitable. Atmospheric effects (due to the ionosphere and troposphere) are the dominant error source for stand-alone GPS and may introduce positioning errors of 4 m or more (one-sigma). To enable sufficient accuracy for precision applications, differential GPS (or DGPS) is commonly used. These systems employ reference stations, whose position is already known, to estimate common-mode GPS errors. Because atmospheric errors are, under nominal conditions, highly correlated over distances of hundreds of kilometers, DGPS almost entirely eliminates these errors. In the United States, the 13

Federal Aviation Authority (FAA) has commissioned two types of DGPS for aviation applications: the Local Area Augmentation System (LAAS) and the Wide Area Augmentation System (WAAS). These systems not only improve accuracy, but also provide monitoring to detect and exclude satellites experiencing anomalous conditions [6]. LAAS is a navigation system intended to support safe aircraft landing. In LAAS the reference station is actually located at the airport, and it serves approximately a 50 mile radius. As the differential corrections are highly accurate for the relatively small area of operation, the position accuracy is normally within 50cm. WAAS uses approximately 30 ground stations stationed in the US, Canada and Mexico. Corrections are interpolated on to a grid and broadcast to users in North America via a pair of communications satellites in geosynchronous orbit. These geosynchronous satellites broadcast differential corrections via GPS like signals. Due to the more widespread nature of WAAS, the accuracy in position obtained is only about 3m (one sigma). LAAS and WAAS perform different functions and are complementary in nature. While WAAS can be used over the entire North American continent and provide en route navigation as well as precision approach (but not precision landing), LAAS is useful for precision landing in zero-visibility conditions within a few miles of the base station.

1.2 GPS Signals Each GPS satellite transmits continuously using two radio frequencies in the L-band referred to as Link 1 (L1) and Link 2 (L2) [3, 7]. The L-band covers frequencies from 1 GHz to 2 GHz. The frequency of L1 is 1575.42 MHz and that of L2 is 1227.60 MHz. L1 is for civil users while L2 is for military users. Each signal consists of three parts. The first part is called the carrier signal. It is a radio frequency (RF) sinusoidal signal. The second part is the code signal from which the pseudoranges can be calculated. It is a binary sequence of zeros and ones which let the user determine the signal transit time. The third part is the navigation data which is also a binary coded message consisting of data about the satellite’s position (orbital position), the satellite’s health, and clock bias parameters. The baseline GPS satellite constellation consists of 24 satellites. At present there are 31 satellites in orbit. Generally 6-8 satellites are visible to a receiver at any point of time. Figure 1.2 shows the orbits of the GPS satellites. All the satellites are approximately at an altitude of 24,000 km above the sea level. Let us look into how GPS users calculate position. The GPS receiver receives the signal from all the satellites in view above the horizon. By identifying the satellite specific code, the GPS receiver determines which signal was transmitted from each satellite. The receiver continuously estimates signal transmission times from each satellite subject to an unknown timing bias associated with its own internal clock. Assuming that the speed of the signal is the speed of light, the range of the receiver is determined as the product of speed and time of travel.

14

Figure 1.2: GPS satellite constellation

To find the three-dimensional position and the time bias for a GPS receiver, a minimum of four equations are required. As more than 4 satellites are almost always visible (each satellite contributing one equation), a least-squares solution can be used to reduce measurement error. An iterative (Newton-Rhapson) solution for the position is used since the equations are nonlinear. The nonlinear equation associated with each satellite has the following form.         … 1.1

Equation 1.1 is written for the kth satellite. ρc is the range calculated from the code sequence for a particular satellite. The variable b is the error in range due to the time offset between the satellite clock and the GPS receiver clock while ε is the noise associated with the range. The actual range of the satellite is denoted by r and is given in equation 1.2, which is a nonlinear function of the user position (xr(i), yr(i), zr(i)) and the satellite position (x(k), y(k), z(k)).                        … 1.2 Newton-Raphson iteration generally converges within 10-12 iterations when the receiver is first activated and in a single iteration when a prior estimate of position is available. In concept, the carrier signal can also be used to determine user position. The GPS receiver gets the carrier sinusoidal waves continuously. The Doppler shift of the carrier wave can be used to determine the relative velocity of the satellite and the user. This velocity can be integrated to determine a change in range over time, as the wavelength of the carrier phase is known (about 19cm for L1). The constant of integration cannot be found, however, as the initial range is not known. Due to the structure of the GPS signal, this ambiguity is always an integer number of wavelengths and can thus be written as λ*N where λ is the wavelength of the signal and N is known as the integer ambiguity.

15

The accuracy of the navigation solution is dependent on the magnitude of ερ(k) in equation 1.1. The smaller its value the more accurate will be the position solution. The error terms in the code and carrier pseudorange measurements are summarized by the following equations.             !  … 1.3

#          $  !$  %& $ … 1.4

ρ = pseudorange from code measurements ϕ = pseudorange from carrier measurements r = true range λ = wavelength of carrier signal N = integer ambiguity δtu = receiver clock bias δts = satellite clock bias Iρ, Iϕ = Delay due to the Ionosphere Tρ, Tϕ = Delay due to the Troposphere ερ, εϕ= errors due to receiver noise and multipath

These equations show that the carrier and code measurements are the algebraic sum of the true range of the satellite from the GPS receiver plus errors. The term c[δtu – δts] represents the error in range due to the error in time estimated by the satellite and receiver clocks and is same as the variable ‘b’ in equation 1.1. Together, the ionosphere and the troposphere delays generate the majority of the GPS measurement error. These systematic atmospheric effects will be discussed in significant detail in the next chapter. Random errors (ερ, εϕ ) make up a much smaller fraction of the total GPS measurement error. It turns out that although ϕ has an ambiguity λN, its noise εϕ is about 100 times smaller than ερ, the noise in code measurement. Hence carrier measurements provide precision (tightly clustered errors) but not accuracy (due to the unknown integer ambiguity). To take advantage of the precision of the carrier signal, aviation applications frequently combine the code and carrier measurements using the Hatch filter, which is defined by (1.5).    

1 91    : ;  <=   #   #<=  … 1.5 9 9

ρs = smoothed range M = filter length (M is generally set as 100 for 1 Hz signal)

The Hatch filter is a recursive filter which is basically a linear combination of the raw pseudorange measurement at an epoch ti and the smoothed pseudorange from the previous epoch, updated (to account for user motion) to the current epoch using carrier measurements. As the carrier measurements are less noisy, this filter enables moving platforms to significantly attenuate noise without introducing a lag in the range measurement. If the ionosphere and troposphere delays are subtracted from both the carrier and the code measurements using differential corrections, the smoothed ranges can achieve an accuracy of tens of 16

centimeters (one-sigma). This nominal level of accuracy is more than sufficient to support precision aircraft landing in zero visibility.

1.3 Ionosphere Threat Under nominal conditions, the ionosphere delay can be fully removed using DGPS. The ionosphere is the earth’s uppermost layer of atmosphere from about 50km to 1000km above sea level. The ionosphere consists of charged particles ionized by solar radiation. The ionosphere delay changes with time of day and level of solar activity. Ionosphere delays are highly correlated for users spread over a large region (e.g. over the same time zone). Under rare ionosphere storm conditions, DGPS fails to remove ionosphere errors. Storm conditions correspond to high levels of ionization caused by intense solar radiation. These conditions are more likely to occur during a period of elevated solar activity which repeats with each eleven year solar cycle. As a result of intense solar activity, steep spatial gradients may occur in the ionosphere. As a consequence of these gradients, a DGPS reference station may measure significantly different ionosphere delays than a mobile GPS user. Under these conditions, differential corrections broadcast by the DGPS reference station do not cancel out the user ionosphere error.

In addition to this spatial gradient error, ionosphere storms also affect the Hatch Filter. Since the ionospheric delay in the code measurement is exactly equal and opposite of that in the carrier measurement, the two signals start diverging from each other. This divergence introduces a bias error into the Hatch Filter at each time step. These biases are accumulated into the filter’s smoothed pseudorange estimate. Figure 1.3 shows how, for a given ionospheric gradient, the

smoothed range output by the Hatch Filter lags behind the true range.

17

50 'Iono' profile Hatch filter output

Ionospheric delay (m) ---->

40

30

B

20

10

0

-10

0

500

1000

1500 2000 Time (s) ---->

2500

3000

3500

Figure 1.3: Divergence bias error in Hatch filter

In steady state, the accumulated bias, labeled B in the figure, can be computed using

equation 1.6 [8].

M  29  1

N ∆ … 1.6 N

To ensure safe landing, new methods are required to mitigate the divergence bias, B. In concept, B could be decreased by reducing the filter time M. However, reducing M also increases the filter’s output noise. Because landing safety requirements do not tolerate an increased level of filter output noise, this approach is poorly suited for precision landing applications. An alternative approach to removing the ionosphere error is to use multiple frequencies [9]. The rate of change of ionospheric delay can be shown to be a function of the frequency of the signal. Hence if the same pseudorange is measured on two frequencies, the ionospheric delay Iρ can be calculated and subtracted out. Unfortunately, dual-frequency signal processing methods are not suited for aviation applications, because the L2 signal (primarily a military signal) cannot be robustly tracked by civilian receivers and because the L2 signal (unlike the L1 signal) does not lie in a portion of the frequency band protected for aviation applications by international law. For these reasons, a single frequency method of mitigating ionosphere errors is desired for civilian aviation applications. It should also be noted that reference-station monitoring is not a sufficient solution for ionosphere threat mitigation. The weakness of ground-based monitoring is that it cannot assure timely alerts against all storms. In particular, the storm gradient might affect only the mobile user and not the

18

ground station. In this case, the gradient would introduce a severe user error that would remain completely unobservable to the ground-based reference station.

1.4 Previous work The errors due to ionospheric storms are single largest threat to LAAS and have stopped it from being commercially operational for over three years now. It is the only threat that DGPS cannot mitigate. There are broadly three approaches to protect LAAS from ionospheric storms. The first approach is to bound the errors by using larger σ values to model errors [10, 11, 12]. However this is a conservative approach in which availability is sacrificed to protect users from anomalous events. The second approach is to monitor for these anomalous events. Dwarakanath V. Simili and Boris Pervan designed a Code-Carrier Divergence (CCD) monitor to detect these extreme gradients [13]. With the CCD monitor an affected satellite can be detected and ‘turned off’ from position calculation of the aircraft. Todd Walter, Sam Pullen, Jason Rife, Jiwon Seo, and Per Enge suggested using a local area monitor (LAM) independently at the airport to augment WAAS to detect ionospheric storms without adding any noise to the nominal GPS signal [14]. Todd Walter also proposed a monitor which used two filters of different length and on the basis of their difference in the reading of the Ionospheric delay, detected ionospheric ramp [8]. The third approach is to change the signal processing block in such a way that is becomes insensitive to ionospheric storms. Hiroyuki Konno, Sam Pullen, Jason Rife, and Per Enge designed two filters to eliminate ionospheric errors and divergence errors [15]. Ionospheric delay is dependent on electron content in the ionosphere and on signal frequency. Given there are two frequencies at which the GPS information is broadcast it is possible to solve for the ionospheric delay. However there are two potential problems to using two frequencies. Firstly the proposed L5 frequency band as the second frequency band is a long way away from certification by FAA. Secondly by using the two frequencies, the resulting noise in the signal is increased. A Kalman filtering estimator structure has also been proposed to estimate ionospheric delay, which could be used as a monitored variable or a correction term of the Hatch Filter [16]. Much less research has been done to improve the signal processing block at the receiver end using only one frequency. While working with one frequency the absolute ionospheric delay cannot be calculated and eliminated without introducing the high noise of unfiltered code measurements. The Hatch Filter has been found to reduce errors excellently under all conditions except under anomalous ionospheric condition during which it adds errors. Hence the real challenge is to find a filter which will preserve the highly desirable characteristics of the Hatch Filter while improving its less desirable properties.

1.5 Contribution The primary aim of this thesis is to formulate new filter algorithms that do not elevate nominal random noise but that reduce the divergence error, B, using only signals broadcast on the GPS L1 frequency. To this end, the thesis derives two new filtering approaches to serve as alternatives to the conventional Hatch Filter for precision landing applications.

19

The first approach taken is to optimize a general linear filter to give the least divergence error B. This approach is discussed in detail and experimentally validated in chapter 4. In this method the ease of implementation of the Hatch Filter is sacrificed, but advancement in computation power and data storage capacity justify this approach. In the second approach, discussed in chapter 5, the ionosphere profile is estimated using a nonlinear, two-slope approach. After estimating the ionosphere profile the divergence error B can be estimated using equation 1.6 and subtracted from the smoothed range output of the Hatch Filter. The first approach provides a 10-20% reduction in the size of the ionosphere divergence error relative to that of the Hatch Filter (for matched levels of output noise). The second approach provides a similar reduction in the size of the ionosphere divergence error, but also drives divergence errors to zero in steady state to reduce the time that landing aircraft are exposed to ionosphere storm threats.

20

Chapter 2 GPS errors 2.1 Sources of errors The major sources of errors experienced by GPS can be classified into three categories. The first type of error arises from the inaccuracies of data broadcasted by GPS satellites. The second type of error is related to the precision and accuracy of the instrument at the user end and is called user error. The third type of error arises from the uncertainties of the medium through which GPS signal travels. This error is called the propagation error and is the largest of the three sources of error. The errors of GPS can be due to error in the estimation of the speed of electromagnetic waves, c, transit time at the satellite, receiver time or the estimation of the position of the satellite itself. As described in chapter 1 the message from the GPS satellite contains the transmission time, t0, when the message is broadcasted. Let t be the time when this message is received at the GPS receiver. Hence the pseudorange ρ of the satellite from the receiver is calculated as given in equation 2.1.     Q  … 2.1 Where c = speed of light in vacuum = 299792458 m/s The control segment of the GPS architecture broadcasts two types of parameters, satellite clock data and ephemeris data. The ephemeris data describe the orbit and position of the GPS satellites. The satellite clock and orbital positioning errors associated with these parameters are in general accurate within a meter. The clock suffers from phase bias, frequency bias and frequency drift. A kalman filter is used to estimate and predict these values [3]. The ephemeris data or the position co-ordinates of the satellite are also vulnerable to errors. GPS satellites follow perturbed keplerian orbits which are specified by thirteen parameters. With these thirteen parameters the co-ordinates of the GPS satellites can be predicted as a function of time. However the error in this calculation grows as the parameters become old. Hence it is necessary to update the parameter at a particular frequency. In general, the control segment of GPS updates the parameters once every 3 hours. The control segment error is eliminated to a large extent in Differential GPS (DGPS). The user error can be divided in two parts, receiver noise and multipath. Receiver noise is introduced by the antenna, amplifiers, cables, multi-access noise (two signals broadcasted in the same frequency bandwidth interfering) and signal quantization noise. Multipath error occurs due the fact that the GPS signal can reach a receiver directly as well as after reflecting from some object, causing interference at the receiver. On land the GPS signal can reflect from the ground, trees etc, while in air when the receiver is on an aircraft, the GPS signal can reflect from the aircraft’s body. To minimize the multipath error, the surfaces from which the GPS signal can reflect are reduced. Also, signal processing algorithms are implemented at the receiver to reduce multipath noise. Typical multipath noise in pseudorange 21

measurements are 0.1-5m depending on the surroundings and the quality of the receiver. Aircraft GPS receivers experience multipath error less than 1 meter. In comparison, carrier multipath errors are in the range of 1-5cm. By using the carrier and the pseudorange measurements in union, the Hatch filter does a very good job of mitigating multipath and receiver noise for users in motion. The third type of error is due to the medium through which the GPS error propagates. The speed of light in the atmosphere differs slightly from the speed of light in a vacuum. These deviations introduce small errors, in that the GPS processing algorithms assume a fixed wave speed, and somewhat larger errors, in that continuous variations in the speed of light in the atmosphere cause the ray path to bend and increase its length. Some path deviation occurs in the troposphere, the lowest portion of the earth’s atmosphere extending from sea level to about 10Km. In this region the speed of light depends on the moisture content of the air, which is a function both of altitude and weather. Even larger deviations occur in the ionosphere, a region of ionized particles, extending from about 50 km to 100 km. Due to the ionosphere’s high electron content (electrons in a column along the ray path), the radio waves ‘slow’ down and bend. The resulting measurement errors, called ionospheric delays, are related to electron content as follows. $    

40.3 R !ST … 2.2 U

40.3 R !ST … 2.3 U

Where TEC is the total electron content: the number of electrons in a tube of 1 m2 cross-section extending from satellite to receiver, and f is the frequency of the radio wave which is broadcasted. Under nominal conditions, troposphere and ionosphere errors can almost entirely be eliminated through use of differential GPS (as discussed in Chapter 1). Table 2.1 summarizes the standard deviations of errors from different sources.

Control segment errors

User errors

Propagation medium errors

Satellite clock and Ephemeris error

3m

Receiver Noise and Multipath (pseudorange measurements) Receiver Noise and Multipath (carrier measurements)

~1 m

Ionospheric delay (modeling error) Tropospheric delay (modeling error)

5-10 m or more 5-10 cm

Table 2.1: GPS errors

22

1-5 cm

2.2 Ionosphere Divergence A key issue for GPS processing is that the ionosphere introduces equal and opposite delays on the code and carrier signals, as indicated by equations 2.2 and 2.3. As a result the carrier and the code pseudorange measurements diverge from each other when the ionosphere delay changes. Due to this phenomenon, the Hatch filter can develop a lag in its prediction of range which is called the divergence bias error. GPS users experience significant changes in ionosphere delay, even over the course of a nominal day. The total electron content depends on the ionospheric ‘weather’ on any particular day and the elevation angle seen by the user. The elevation dependence results from the fact that the signal from a satellite on the horizon will travel a longer path through the ionosphere than that from a satellite vertically overhead. Thus, satellites at low elevations experience greater delays than satellites at high elevations.

12

10 9

10

7 6

Ionospheric Delay (m) --->

Ionospheric Zenith Delay (m)

8

A2

5 4 3

8

6

4

2

2

½A4

1 0

0

5

10 15 Local time (h)

20

0

25

Figure 2.1a:The Klobuchar ionospheric model for nominal ionospheric delay

23

0

0.5

1

1.5

2 2.5 Time (s) --->

3

3.5

4 4

x 10

Figure 2.1b: Nominal Ionospheric delay on 17 January 2008 at CORS station ZOA1

10

10

9

9

8

8 Ionospheric delay (m) --->

Ionospheric delay (m) --->

Differential GPS eliminates divergence errors introduced by ionosphere divergence under nominal conditions. Ionosphere variations for an ‘average’ day are represented by the Klobuchar model [3, 19, 20]. Figure 2.1a shows a pictorial representation of the Klobuchar model. In the illustration, the parameter A2 indicates the maximum amplitude of the delay (which occurs in early afternoon) while A4 represents the length of the daylight hours. During daytime the sun ionizes the gas in the ionosphere and hence the TEC increases above its stable night time value. Figure 2.1b shows the ionospheric delay computed from actual data, recorded by CORS station ZOA1 on 17th January 2008 for 10 hours starting at GMT 0000. Since the ionosphere changes experienced by users in the same time zone are nearly identical, divergence biases can be eliminated under nominal conditions by applying differential corrections from a nearby ground station.

7 6 5

7 6 5

4

4

3

3

2

0

1000

2000

3000 4000 Time (s) --->

5000

6000

2

7000

0

1000

2000

3000 4000 Time (s) --->

5000

6000

7000

Figure 2.2b: Nominal ionospheric delay with an anomalous event due to ionospheric storm

Figure 2.2a: Nominal ionospheric delay

Differential GPS does not necessarily eliminate the divergence bias under anomalous ionosphere conditions. On some days the ionospheric weather deviates from the ‘average day’ by a significant amount. Due to high solar activity the ionosphere may experience strongly varying TEC. These phenomena are called ionospheric ‘storms’. Storms give rise to very high spatial and temporal gradients of ionospheric delay. Figure 2.2 compares nominal changes in ionosphere delay (2.2a) to changes which can occur in the event of an ionospheric storm (2.2b). The nominal ionosphere gradient is drawn from the Klobuchar model for the period of 2 hours shown by the red box in figure 2.1 a. For short periods of time, like 1-2 hours, the Klobuchar model can be approximated by a straight line. However, the straight line approximation is not valid in the event of a storm. Storm-induced divergence bias errors may exceed 10m and are not generally eliminated by differential GPS since the user may observe a storm that the ground station does not (or vice versa). The magnitude of ionosphere-induced errors makes them the single biggest threat to GPS integrity.

24

2.3 Ionospheric threat model

Gradient = 400mm/Km 100 m/s

LGF 50 km Figure 2.3: Ionosphere threat model

Ionosphere storm models have been developed by several researchers [21, 22] to support safety analyses of precision landing using LAAS. These models are based on observations of large ionosphere walls detected with arrays of GPS receivers positioned across North America [17, 23]. The standard ionosphere storm model consists of a steep traveling wavefront between regions of low and high ionosphere delay. The wavefront is modeled as a piecewise linear curve with a specified ionosphere delay gradient, wavefront speed and wavefront width. This standard model is illustrated for representative parameters in Figure 2.3. For the purposes of safety analysis, storm parameters (delay gradient, wavefront speed and width) are bounded by a threat space derived from GPS data. Due to an ionospheric storm the ionospheric delay, which is directly proportional to the total electron content (TEC), can have a gradient of 2 - 400 mm/km. This ionospheric ‘wall’ can have a speed of 0 – 1000 m/s. The threat model limits establish a 25 km lower bound on front width and no upper bound. In addition to describing the storm, the threat model also includes parameters that describe aircraft motion. Aircraft generally approach the airport with an average speed of 60-100 m/s. The compass angle between the ionosphere ‘wall’ velocity and the aircraft velocity can be 0-3600. Safety analyses combine these parameters to bound worst case scenarios. In practice, some logical simplifications are implemented to make the analysis tractable. In the worst case scenario the ionospheric ‘wall’ only affects the aircraft. In this scenario, the ionospheric ‘wall’ comes from behind the aircraft and overtakes it. Accordingly, ground monitors at the reference station [13] are unable to detect the storm. Both the ionospheric ‘wall’ and the aircraft have a 25

velocity but since the rate of change of ionosphere delay produces the error, it is the relative velocity of the aircraft with respect to the ionospheric ‘wall’ that is important. For mathematical convenience the aircraft speed was taken to be 100 m/s and the anomalous ionosphere spatial gradient was taken as 400 mm/Km. Because of the extended smoothing time associated with the Hatch Filter, the divergence bias error persists for several minutes even after the aircraft exits from the ionosphere front. The piecewise linear model for ionosphere storms is adopted in this thesis to enable design of algorithms for divergence bias mitigation. Equation (2.5) gives the mathematical representation of this ionospheric model. Without loss of generality, the time when the user transitions into or out of the wavefront can be identified as the epoch k = 0. The constant d is the divergence rate of the carrier and pseudorange measurements; the divergence rate is equal to twice the ionospheric delay gradient. The variable Δt is the discrete time interval between successive measurements. This thesis focuses on values of Δt at or near the value used for LAAS (0.5 second).  d   I 0 +  ∆t  k Ik =  2   I0 

k >0

... (2.5)

k ≤0

To get the worst case it is assumed here that the aircraft remains under the influence of the ionosphere gradient till it lands and the LGF is hit by this ionospheric ‘wall’ only at the last instant. For this worst case scenario the rate of change of ionospheric delay is 40 mm/s (ionospheric gradient times relative velocity = 400 mm/Km X 100 m/s) and the divergence bias error is 8 m for a 100 second Hatch filter. It is logical that although the ground station is initially oblivious to the anomalous event, as the aircraft approaches the airport, ultimately the ground station also sees the same ionospheric delay as the aircraft. Hence the ionospheric delay is corrected by differential corrections as the aircraft comes closer to the airport. However the divergence bias error, which persists for sometime even after the aircraft comes out of the ionospheric gradient, is capable of introducing some serious errors in the GPS estimated position of the aircraft. Small variations in aircraft speed and ionospheric ‘wall’ speed can vary this error from 5m to 10m. Clearly this error is unacceptable to LAAS. Hence the thesis aims to develop filters that produce smaller divergence bias errors than the Hatch filter without sacrifice the noise performance of the Hatch filter.

26

Chapter 3 Data sets As discussed in the previous chapter GPS signal errors consists of measurement errors, receiver noise, multipath, ionospheric delay and tropospheric delay. In this thesis ionospheric delay is the error being focused upon, the rest of the errors are treated as unmodeled errors. However it is important to analyze the effect of the unmodeled error on the designed digital filters. These errors cannot be modeled as white noise and hence need to be deal separately. Hence it is necessary to test the filters on real GPS data. Data was obtained from two sources. The first source was online data which is uploaded by reference stations which monitor GPS. These reference stations have standard GPS receivers of high accuracy and equipped multipath elimination algorithms. The second source was a self-made experimental setup to acquire GPS signals autonomously. Although due to the surrounding this experiment is likely to encounter more multipath errors and receiver noise, it provides a more severe test for the digital filters to perform. Hence the experimentally acquired data set augments the vigorous analysis of the proposed digital algorithms.

3.1 Web based GPS archives Reference stations are set up all over the US to monitor and record continuously code and carrier measurements. The National Geodetic Survey (NGS) operates two network of continuously operating reference stations (CORS) and their recorded data is available through the internet for researchers. To make the data readable, strict protocols and formats of data need to be maintained. One of the popular formats is called the “Receiver Independent Exchange Format” (rinex) format. The rinex data files begin with a header which consists of useful information like the date, frequency of GPS signals, data provided, xyz position co-ordinates of the reference station GPS receiver, sampling rate etc. Subsequently, for each measurement epoch, pseudorange and carrier measurements are listed for each visible satellite. The CORS data sets used for this work were sampled at 1 Hz. A sample rinex file is shown in figure 3.1. In the header it is mentioned that the data recorded are L1, L2 (pseudoranges in L1 and L2 frequency), P1, P2, C1 (carrier phase measurements in L1 frequency), S1 and S2. The date is given as “07 5 9” which means 9th May 2007. The time is 0 hours, 0 minutes and 0 seconds (midnight or beginning of the day) given by “0 0 0.0000000”. “10G” signifies that 10 satellites are visible at the current epoch and then the PRN of each visible satellite is given, each followed by “G”. After specifying the visible satellites the data (L1, L2, P1, P2, C1, S1, S2) are given separated by tabs. For this thesis, rinex data from reference station ZOA1 in Oakland, CA were used extensively. Matlab was used to process the rinex files and extract the code and the carrier measurements for each satellite.

27

Figure 3.1: Sample Rinex file

The noise of code measurements for all the visible satellites from ZOA1 on August 28, 2006 during the first hour of the day indicate that the levels of multipath and noise levels differ among satellites. PRN 8, 11, 17, 24 and 26 were the high elevation satellites and their noise levels were below 1m of 1σ level. Low elevations satellites, namely PRN 9, 27, 28 had relatively high noise of 1σ greater than 1.5m. As observed in figure 3.2 PRN 17 was severely affected by relatively high frequency multipath while PRN 24 and 28 represent relatively lower frequency multipath errors.

28

1. PRN 8

2. PRN 9

3

6 4

2 2 0 Cose Noise (m)

Code Noise (m)

1

0

-2 -4 -6

-1

-8

-2 -10

-3

-12

0

500

1000

1500 2000 Time (s)

2500

3000

3500

0

500

1000

1500 2000 Time (s)

2500

3000

3500

2500

3000

3500

2500

3000

3500

1. PRN 17

3. PRN 11 2.5

6

2

4

1.5 1 Code Noise (m)

Code Noise (m)

2

0

0.5 0 -0.5

-2

-1 -1.5

-4 -2

-6

-2.5

0

500

1000

1500 2000 Time (s)

2500

3000

3500

0

500

1000

1500 2000 Time (s)

6. PRN 26

5. PRN 24 15

3

10

2

5

Code Noise (m)

Code Noise (m)

1 0

-5

0

-1 -10

-2

-15

-20

0

500

1000

1500 2000 Time (s)

2500

3000

3500

29

-3

0

500

1000

1500 2000 Time (s)

8. PRN 28 2.5

3

2

2

1.5

1

1 Code Noise (m)

Code Noise (m)

7. PRN 27 4

0 -1 -2

0.5 0 -0.5

-3

-1

-4

-1.5

-5

-2

-6

0

500

1000

1500 2000 Time (s)

2500

3000

3500

-2.5

0

500

1000

1500 2000 Time (s)

2500

3000

Figure 3.2: Code measurement noise on simultaneously all visible satellites

3.2 Experimentally acquired data For additional validation, GPS data were also acquired on the Tufts campus using a low-cost commercial GPS receiver. The clock on board has a high decay or gain rate and needs to be reset frequently. This results in errors in the position solution of the GPS receiver. The receiver is also required to maintain a lock on the carrier signal to ‘count’ the number of cycles received. A typical Doppler tracking error, called a cycle slip can result in an error of one wavelength which is around 19 centimeter. Hence the noise of the GPS data acquired by the commercial GPS receiver is expected to be noisier than CORS data. 3.2.1 Device specification The Garmin GPS 18 USB was the GPS sensor used for this experiment. As the name suggests it interfaces to the computer with an available USB Port. It is a 12 channel GPS receiver which means that it has the capability to receive GPS data continuously from 12 satellites simultaneously. It operates using a FLASHbased program in non-volatile memory. It is pretty compact with a size of 61 mm in diameter and a weight of 100 grams. It takes an input of 4.4 - 5.5 Volts and has current characteristics of 5mA @ 5V. Its data acquisition rate is approximately 1 record per second. Its position accuracy is less than 15meters with 95% probability and its velocity accuracy is 0.1 knots RMS value at steady state. The GPS 18 USB is compatible with both USB 2.0 full-speed protocol as well as with USB 1.1 full-speed protocol. The device driver required for this hardware was downloaded from http://www8.garmin.com/support/commProtocol.html. The GPS 18 USB communicates with the computer using data packets. The packet has the structure as described in Table 3.1.

30

3500

Figure 3.2: GPS 18 USB

Byte number 0 1-3 4-5 6-7 8-11 12+

Byte Description Packet Type Reserved Packet ID Reserved Data Size Data

Notes USB protocol layer = 0, Application layer = 20 Must be set to 0 *set to 10 to specify it’s a command Must be set to 0 *set to 2 for commands *set to 49 or 110 to receiver pseudorange and carrier data Table 3.1: USB packet structure

For our experiments, as the application layer is always used, the 0th byte of the packet is set to 20. The asterisks denote that for the present experiment those values were used but for other applications, other values can be used for those bytes. Each record ends with the data record, which can be of any length. The GPS 18 USB requires a trigger to start receiver data collection. Two commands are needed to turn on the sensor. A packet {20, 0, 10, 0, 2, 49} turns on the position record, which gives the receiver’s location in latitude, longitude and altitude. This packet is used for verification that the GPS sensor is working. Packet {20, 0, 10, 0, 2, 110} turns on the receiver measurement record which contains the pseudoranges and carrier measurement records. For the current experiment a laptop with standard Pentium Dual Core processor of speed 1.6 GHz was used, and the program for communicating with the GPS sensor was written in C and can be found in the Appendix. 3.2.2 Data representation Using the procedure described above, pseudorange and carrier phase data were collected from an arbitrarily chosen satellite (in this case satellite PRN 32). The noise in pseudorange was estimated by fitting a four degree polynomial curve to 100 seconds of data. 36 continuous chunks of data were analyzed in this procedure to obtain an average noise performance over an hour. Figure 3.4 shows the noise from an experimentally based data and its standard deviation was 1.755m. The experimental data is significantly noisier than the web-based data due to the quality of the receiver. The experimental data also suffered a great deal from cycle slips which introduced errors in the carrier phase measurement.

31

6

Pseudorange noise (m) ---->

4

2

0

-2

-4

-6

-8 0

500

1000

1500 2000 Time (s) ---->

2500

3000

3500

Figure 3.3: Pseudorange noise from experimental data of Satellite PRN 32

32

Chapter 4 Optimum Linear Filter Design 4.1 Problem analysis For LAAS to be certified by the FAA, it has to be shown that the probability of a missed landing is less that once in a billion. This requires us to consider even the rare possibility of an ionospheric storm and remove the errors caused by it. These storms introduce large errors which differential GPS like LAAS cannot remove. Reducing the sensitivity of LAAS to these errors is the major challenge to prove reliability. This thesis proposes two signal processing methods that mitigate the ionosphere threat to LAAS. The first is an optimal linear approach, described in this chapter, which minimizes the maximum transient error (which directly impacts standard availability analysis [17]). The second is a nonlinear approach, described in the next chapter, which attempts to obtain a similar maximum transient error while minimizing error in steady-state. The ionosphere divergence error that affects LAAS users depends strongly on the filter used to smooth GPS measurements. The standard filter used in LAAS to attenuate receiver noise and multipath is the Hatch Filter. Under nominal conditions the Hatch Filter does an excellent job and the resulting output noise in generally within 1σ equal to 50 cm. The data sets discussed in Chapter 3 are representative of this, but that both the CORS data and the GPS18 data are higher noise than can be achieved with specialized hardware used in LAAS where nominal noise is 1 σ equal to 20 cm. However in the presence of high ionospheric gradients, the Hatch Filter suffers from significant errors. The general structure of the GPS signal is given in equation 1.3 and 1.4. Studying equations 2.3 and 2.4, it is observed that the ionospheric delay in pseudorange measurement is equal and opposite of that in the carrier signal measurement. Hence the GPS signal structure at epoch k can be restructured as in equation 4.1 and 4.2.       !    V

,

… 4.1

#      !    %&  V$, … 4.2 All the symbols have their usual meanings (see Chapter 1). Both observables are also subject to systematic delays caused by the ionosphere (Ik), by the troposphere (Tk) and by the clock difference between the satellite constellation and the user clock (δtk). The key thing to note here is the sign reversal of Ik for the two observables and the integer ambiguity N in the carrier signal measurement. Let us take a closer look at the errors caused by the ionospheric divergence. Ionosphere divergence occurs even on nominal days in the absence of severe ionosphere weather. Figure 4.1a shows GPS code measurement errors on a typical day. On a nominal day, the ionosphere delay changes relatively slowly (over a time scale of hours). To separate out the ionosphere error (subject to an integer bias), the smoothed pseudorange needs to be estimated. It is calculated by fitting a 5th order polynomial to the nuisance factor (2Ik-λN), obtained by subtracting the carrier signal from the pseudorange as shown in 33

equation 4.3, and then adding the polynomial to the much less noisy carrier signal to get a good estimate of the smoothed code measurement. This estimate can be used to compute the ionosphere divergence rate and the associated divergence error introduced by filtering. Even on a nominal day, the divergence error can be significant, as illustrated by Figure 4.1a. These errors have no impact on LAAS applications, however, as differential corrections mitigate errors due to nominal ionospheric delay.   #  2  %&  V

,

 V$, … 4.3

Divergence errors that occur during ionosphere storms are especially dangerous, because differential corrections do not necessarily remove the associated errors. Figure 4.1b illustrates an analysis of anomalous ionospheric activity during a simulated ionosphere storm. The simulated ionosphere storm is obtained by adding the threat model (equation 2.5) to the raw data illustrated in Figure 4.1a. Consider the case in which an aircraft observes an ionosphere anomaly (Figure 4.1b) while the ground reference station does not (Figure 4.1a). Although the Hatch Filter again smoothes out the sensor noise, in the process it introduces a tracking error which cannot be differentially corrected. Because the filter has “memory,” the divergence error remains for a time even if the aircraft leaves the region containing the ionosphere anomaly. Thus, even as the aircraft draws in closer to the airport, where the propagation medium will eventually become the same for both receivers, the divergence bias error remains in the system unchecked. Given the large size of these divergence errors, they are the biggest threat to LAAS integrity. Hence the motivation here is to optimize a linear filter to get the minimum transient divergence bias error. 3

14 Original error due to nominal Ionospheric delay Error after passign through the Hatch filter

2.5

Original error due to anomalous ionospheric delay Error after passing through the Hatch filter

12 10

2 Errors (m) ---->

Errors (m) ---->

8 1.5

1

6 4 2

0.5 0 0

-0.5 0

-2

500

1000

1500 2000 Time (s) ---->

2500

3000

3500

-4

0

500

1000

1500 2000 Time (s) ---->

2500

3000

3500

Figure 4.1b: Anomalous GPS errors (sensor noise + anomalous ionospheric delay + divergence bias error)

Figure 4.1a: Nominal GPS errors (sensor noise + nominal ionospheric delay)

4.2 Hatch Filter – A linear filter The optimization process minimizes the maximum transient error for a given noise attenuation level. To get the optimal linear filter, first a generalized linear filter needs to be designed. A linear filter by definition is a linear combination of its inputs. Hence it is logical to look more closely at what input the Hatch Filter takes. The Hatch Filter uses the smooth (noiseless) characteristics of the carrier phase 34

measurements and the absolute value of range from the code measurements. Restated in equation 4.4 the Hatch Filter attenuates the high noise level of the pseudorange measurement by multiplying it with a constant γ which is equivalent to 1/M in equation 1.5, and updates the change in range by the change in the carrier phase measurement. The carrier phase measurement difference adds little noise to the output and in the process also gets rid of the integer ambiguity N. W,  X  1  XYW,<=  #  #<= Z … 4.4

The divergence error V , and noise attenuation for the Hatch Filter can be characterized using its impulse function. Any linear operator can be expressed as the convolution of an impulse function with the input signal or signals. Expressing the Hatch Filter as a convolution of an impulse function with the input code and carrier signals gives the following. ]

W,  #  [ \W, <  #<  … 4.5 ^Q

To further condense this form of the Hatch Filter a new variable Carrier-Updated-Code (CUC) is introduced in equation 4.6. Then the Hatch Filter can be expressed as in equation 4.7. The CUC observable applies a carrier-phase difference that updates the range measurement from epoch i to match the value of the current range, rk, for epoch k [1]. T_T,  <  #  #< … 4.6 ]

W,  [ \W, T_T, … 4.7 ^Q

The CUC variable updates the lag in the aged pseudorange measurement by the highly accurate carrier phase difference. Equation 4.8 compares the current pseudorange measurement and the CUC variable. Although the range measurements match, the ionospheric delay characteristics and the noise characteristics are slightly different. Given that the carrier noise is an order of magnitude less than the code measurement noise, the CUC noise is approximately the same as the code measurement noise. The ionospheric delay bias under nominal conditions is taken care of by differential corrections. T_T.    2<     V

,

 V$,  V$,  … 4.8

To find out the value of αH,I , equation 4.3 needs to be expanded and the coefficients of the CUC variable needs to be studied. With a little bit of mathematical manipulation after expanding equation 4.4, one can derive equation 4.9. W,  X1  XQ 
35

Looking at the coefficients of equation 4.9 i.e. the coefficients of the CUC variable for i=0 and i=1, it can be derived using induction that the value of αH,I is given in equation 4.10 (as pictorially represented in figure 4.2). \W,  X1  X … 4.10 In the actual implementation of the Hatch Filter, the recursive form in equation 4.3 is used rather than that given in equation 4.7. This is because in the recursive format, only data at two different epochs need to be considered at each time while in the impulse function format, in theory, data at an infinite number of epochs need to be considered for each computation. Even an approximate implementation of the Hatch impulse function that truncates filter length must still process as many as 2/γ coefficients or more (in order to ensure minimal truncation errors). However in comparison to the contemporary computational and memory storage technology available to users, the requirement of extra storage memory and computation in the impulse function approach is insignificant. 0.01 0.009 0.008

Coefficient value --->

0.007 0.006 0.005 0.004 0.003 0.002 0.001 0

0

100

200

300 400 Coefficient index, i --->

500

600

Figure 4.2: Coefficient of Hatch Filter impulse function

36

700

After having generalized the Hatch Filter, it is important to study its response to the anomalous ionospheric storm model defined in chapter 2. Using equations 2.5, 4.1, 4.2 and 4.6 the impact of the ionospheric gradient on the CUC variable is depicted in figure 4.3. The resulting ionosphere bias in the CUC observable, Ii,k , reflects both the positive contribution of the code and the negative contribution of the carrier.

60

Code Error, I Iono Error Relative to Io, normalized by d∆t

~

40

CUC Error, Ii|k

20

0

i=0 k = 20

-20

k = 40 k = 60

-40

k = 80 k = 100

-60 -50

0 50 Absolute Time Index, k -i

100

Figure 4.3: Code bias and CUC bias due to anomalous ionospheric storm

In the figure 4.3, the horizontal axis is the absolute time index associated with the aged CUC measurements. Index (k) is the current time epoch and index (i) is the number of epochs that the CUC is aged by. As the absolute time index (k) increases, the ionosphere delay error in the CUC variable at a certain time (i.e. at constant (k-i) value) becomes increasingly negative. This is because with increasing time a CUC variable at a constant time (k-i) needs a longer update from the carrier measurements resulting in more negative ionospheric delay. The mathematical formulation of the CUC variable in equation 4.11 verifies the same fact.  k   I 0 + ( d ∆t )  2 − i  k ≥ k − i ≥ 0     d   I%i|k =  I 0 −  ∆t  k k ≥ 0 ≥ k − i ... (4.11) 2    I0 k <0  

The CUC variable divergence error can be found by subtracting the CUC variable from the code measurement. It should be kept in mind that the ionospheric anomaly was assumed to start at k=0. Using equation 4.1 and 4.6, it is observed that all the other biases drop out from the divergence error 37

except the ionospheric error. Equation 4.12 and figure 4.4 also show that at any epoch k, as the CUC variable is update by more aged carrier phase measurement data, the divergence error increases. However once the CUC variable is updated by a carrier phase data older than the epoch k=0, the divergence error becomes constant. This is logical because ionospheric model assumed that there is no ionospheric delay before epoch k=0. Ecuc,i|k = E ( ρ k − CUCi|k ) = I k − I%i|k  ( d ∆t ) i  = ( d ∆t ) k  0 

k ≥ k −i ≥ 0 k ≥ 0 ≥ k −i k <0

... (4.12)

The tracking error is computed as the difference between the smooth range, output of the Hatch Filter, and the code measurements. SW,  SY  W, Z … 4.13 Now writing the output of the Hatch Filter in term of the CUC variables and using the fact that all the coefficients of the impulse function of the Hatch Filter add up to one, the tracking error can be recalculated from equation 4.13 to get equation 4.14.

CUC Divergence Error, normalized by d∆ t

120

100

k = 100

80

k = 80

60

k = 60

40

k = 40

20

k = 20

0

k=0

-20 -50

0 50 Absolute Time Index, k -i

100

Figure 4.4: CUC divergence error Ecuc,i|k

] SW,  SY∑] ^Q \W, Y  T_T, ZZ  ∑^Q \W, S,| … 4.14

Using the values of αH,I from equation 4.10 and substituting it in equation 4.14, the steady-state tracking error is found by letting k tend to infinity.

38

lim SW, 

g]

1X N∆ … 4.15 X

4.3 Optimized Linear (OL) filter To find a general linear filter, the recursive nature of the Hatch Filter is relaxed. This gives the flexibility to alter the coefficients of a linear filter to be able to minimize the ionosphere induced tracking error, also called the divergence bias error. The structure of the generalized linear filter is given in equation 4.16. This filter takes the same input as the Hatch Filter. i

hij  [ \ T_T, … 4.16 ^Q

The form of the generalized convolution notation of equation 4.16 is analogous to that of equation 4.7, except that the H subscript (indicating the Hatch Filter) has been replaced by subscript GLF (indicating Generalized Linear Filter) and, more significantly, that the length of the generalized impulse function is specified as L+1. Hence to formulate the OL filter, the optimal values of α and L need to be determined. This optimization problem is subjected to two constraints. This first one comes from the constraint that the coefficients must add up to one to conserve the total value of the actual range of the GPS receiver from the satellite. In other words the DC gain should be unity. i

k=  [ \  1  0 … 4.17 ^Q

The second constraint comes from the fact that the noise output characteristics of the OL filter should be same as the Hatch Filter. Noise output characteristics are measured by Γ which is the square of the standard deviation of the output noise to the input noise (assuming input noise variance to be unity). i

k  [ \   l  0 … 4.18 ^Q

This constraint can be derived by noting that the resultant variance of independent Gaussian distributions is the sum of their individual variances. Since the noise at each epoch is scaled by the weighting coefficient, αi, the output variance due to each epoch (for a unity variance input) is αi2. The constraint G2 is derived assuming that each measurement is independent of all the previous measurements. In other words, the noise in the system is assumed to be white. An equivalent output noise constraint for correlated noise can also be derived, given that the noise correlation is known [18]. In this case, a correlated noise constraint, G2,corr, can be derived by convolving the filter impulse function with the impulse function for the correlation process, ηk .

39

i

i

^Q

^Q



k,m  [ n[ \ o< p  l  0 … 4.19 Although recorded data sets indicate correlation in the noise of the code and carrier measurements, the dominate frequencies of the correlated noise vary significantly from one data set to another. To make the optimized filter robust and not effective on only a particular modeled correlation process, ηk, the constraint G2 is used for the optimization process, rather than G2,corr. To compare the off-nominal performance of the OL filter to the Hatch Filter, it is desirable to match the value of the output variance, Γ, for both filters. The value of Γ for the Hatch Filter that satisfies G2, as described by (4.18), can be evaluated as the sum of the squares of the Hatch Filter coefficients, which is a geometric series. In general the sum of an infinitely decaying geometric series has the following form. ]

[ q <= 

r^Q

q … 4.20 1

Using this relation the Γ of the Hatch Filter is given below. The general form of (4.20) can be adapted to match the form of the Hatch coefficients of (4.10) by replacing the variable a with γ2 and the variable r with (1- γ)2. The result is that the constraint (4.18) can be written in the following manner for the Hatch Filter. l

X X  … 4.21 1  1  X 2  X

Using this relationship, the Hatch coefficients can be determined for any desired level of noise attenuation by substituting (4.21) back into (4.10). This provides a way to find a Hatch Filter of the same noise output characteristics as the optimized filter. Hence we next move to the optimization procedure. For a given level of noise attenuation, Γ, the LAAS application requires as small a divergence error as possible. An optimized filter can be obtained by framing a minimization problem that results in the smallest possible maximum transient error induced by the anomalous ionospheric model (equation 2.5), subjected to the constraints G1 and G2. The cost function, F, for the optimization problem consists of the maximum tracking error Ek of all the epochs during the transient phase i.e. from zero (when the ionospheric event starts) to L (when the full length of the filter is exposed to the threat). To determine the optimal filter coefficients, αk, the cost function must be minimized subject to the filter constraints for unity DC gain (G1) and for output filter noise (G2). mins|k= , k  ,

s  maxbSr \ c Uu v w 0, x … 4.22

Lagrange multipliers are used to solve this optimization problem. The method of Lagrange Multipliers notes that the minimum value of a cost function occurs either where the function’s gradient is zero or on the boundaries of its domain. Constraints that are equal to a constant may be scaled and added to

40

the cost function without changing the location of the minimum. Thus for the case with no domain boundaries, a cost function, F, and two constraints, G1 and G2, the minimum cost must occur where ys  %= yk=  % yk  0 … 4.23

Since this is digital system, only discrete epochs are being considered here. Hence the problem is nontrivial, and needs to be solved in two stages. Initial a particular epoch is chosen (epoch L which the transient error just reaches its steady state value) and the maximization function of equation 4.22 is taken out. Later for the full optimization problem, the maximization function is considered. At epoch L, when the divergence bias error reaches its steady state value, equation 4.14 can be used to state the steady state cost function, Fss, which is the divergence bias error in steady state. mins |k= , k  ,

i

s  [ \ S,|i … 4.24 ^Q

Using equation 4.12, the steady-state cost function, Fss, can be calculated as the following. i

s  N∆ [ z\ … 4.25 ^Q

Substituting the values in equation 4.23 and scaling the arbitrary constant λ1 and λ2 by (LdΔt)-1 we get the following. xN∆<=

{s {k= {k z  %=  %   %=  2% \  0 … 4.26 {\ {\ {\ x

Rearranging equation 4.26 the linear coefficients of the OL filter are found. \ 

z %=  … 4.27 2% x 2%

The relationship suggests that the coefficients are a linear function of index i, and they fall on a straight line with a constant slope. To fully determine the value of the coefficients the two constraints G1 and G2 can be used to solve out the values of λ1 and λ2. i

[ ^Q i

[ ^Q

z %=    1  0 … 4.28 2% x 2%

z %=     l  0 … 4.29 2% x 2%

Using the following identities the values of λ1 and λ2 can be solved to get the following.

41

1 i x |x  2} x  1 xx  1  [z  ,[z  … 4.30 2 3 i

^Q

%= 

^Q

1 2 x1 x2   % , %  ~ … 4.31 2 x1 2 12xlx  l  1

Equation 4.27 combined with equation 4.31 gives the value of the coefficients. Since the Lagrange multiplier technique obtains an extremum, it is necessary to determine whether that extremum is a minimum or a maximum. In this case, it is possible to demonstrate that the two solutions for the Lagrange multiplier values correspond to a minimum as well as a maximum. Specifically the positive value of λ2 results in maximization while the negative value of λ2 results in minimization.

-3

8

x 10

Uniform

Trapezoid

Zero Steady-State Error

L = L*

L = Lmax

Value for L* Value for L

6 Filter Coefficients

Triangle

4 2 0 -2

L = Lmin

Lmin < L < L*

Transient Error, normalized by d∆t

150

100

50

0

0

200 400 600 Current Time index, k

800

0

200 400 600 Current Time index, k

800

0

200 400 600 Current Time index, k

800

0

200 400 600 Current Time index, k

800

Figure 4.5: Tracking error for OL filters for different lengths

There is still one parameter to be determined, which is the length (or order) of the filter, L. It is observed that by increasing the length of the filter, the steady-state error decreases. The shortest filter that meets the constraints is uniform, with all its coefficients equal to Γ -1. The uniform filter also has the highest steady state error and figure 4.5 shows the steady state error slowly decreases as the length of the filter increases. At order of the filter L = L*, the filter is triangular in shape (third filter in figure 4.5) and its 42

steady-state error is minimum of all filters whose maximum transient error occur at steady state. As the length further increases some of the coefficients become negative and the maximum transient error start to occur sometime before the steady state. Although for filters with order L > L*, the steady-state error is smaller than the triangular filter (order L = L*), their maximum transient error was more than the steady-state error of the triangular filter. At L = Lmax (fourth filter in figure 4.5), the steady-state divergence bias error is found to be zero. In the next few paragraphs it will be shown analytically that the triangular filter indeed is the optimal filter with the best transient. Accordingly, it is relevant to derive an equation for L*. The value of L* can be found using the argument that the last coefficient i.e αL = 0. Hence from equation 4.31 it is clear that λ1 should be equal to 1. Using the relation between λ1 and λ2 from equation 4.31, L* is solved out as in equation 4.32. Since L* is an integer the solution is rounded down to the nearest integer. The steady state error is calculated using equation 4.25 and the identities in 4.30. xR  U€uu 

1 9 1 2  4  l  ‚  ‚ … 4.32 3l 4 2

s 

xR  1 N∆ … 4.33 3

For the case in which Γ results in an integer value of L* without rounding, the steady-state error of equation 4.33 may be expanded using equation 4.32 and has the following form. s 1 9 1  2  4  l  ‚  … 4.34 N∆ 9l 4 2 Moving on to the second part of the optimization problem, the transient error at all the epochs leading up to the steady state needs to be considered. This results in the optimization problem to be done on each epoch separately. The optimized filter will be characterized by the fact that the minimum tracking error Ecuc,j found by optimizing the filter at epoch j will lower than all other Ecuc,j, i ≠ j. Hence the minimization problem from equation 4.24 can be restated with the addition of a new inequality constraint in equation 4.35 below. For v w 0, x, minYsr |k= , k , S,r ƒ S,„r Z , sr  S,r \  … 4.35 This optimization problem is functionally equivalent to the original transient optimization problem of equation 4.22, with the maximization operation transformed into an inequality constraint. Using equation 4.35, it is possible to step through each epoch of the transient, identifying the best possible filter which results in a maximum transient error at the epoch of interest. The new constraint Ecuc,j ≥ Ecuc,i for i ≠ j, is treated in the following manner. If the solution of the optimization problem given in equation 4.35 using Lagrange multiplier automatically satisfies the new inequality constraint then this constraint is of no consequence. If, however, the unconstrained Lagrange 43

Multiplier solution does not naturally fall into the domain of the inequality constraint, then the optimal constrained solution must lie on the boundary where the inequality becomes equality. In that case the Lagrange multiplier approach much take into account another equality constraint, with Ecuc,j = Ecuc,i for i ≠ j. Hence let’s first solve the unconstraint case. Applying equation 4.23 to the current optimization problem at a particular epoch j, the following equation is obtained. xN∆<=

{sr {k= {k  %=  %  0 … 4.36 {\ {\ {\

Using equation 4.12, the above relation can be evaluated to get the form below.

 1  λ1 , i< j  i −   2 Lλ2  2λ2 αi =  ... (4.37) λ1  1  , j − i ≥ j  2 Lλ  2λ2 2   The coefficient are similar to the form obtained in equation 4.27 except the second relation where i ≥ j. This is the case when these coefficients apply to measurements taken before the onset of the ionospheric storm. As shown in figure 4.4, the divergence error is constant in that region and hence αi being constant in that region is not a surprise. The filter that results from equation 4.37 cannot satisfy the strict inequality constraint that Ecuc,j > Ecuc,i for i ≠ j. This can be proved by contradiction. Assuming that the maximum transient error occurs at epoch j, one can say that Ej should be greater than both its neighboring points. i

0 … S,r  S,r<=  [ N∆\ … 4.38 ^r

i

0 … S,r  S,r†=   [ N∆\ … 4.39 ^r†=

Since all of the coefficients for i ≥ j are identical, and since j is a general index in [0, L] the strict inequality constraints (4.37) and (4.38) may be rewritten as follows. 0 ‡ \i … 4.40

0 ‡ \i … 4.41

As a coefficient cannot be simultaneously positive and negative, the filter transient cannot be a strict maximum at epoch j. As the unconstrained problem did not yield a satisfactory answer, the only possible valid result lies on the boundary of the constraint, which is Ecuc,j = Ecuc,i for at least one i ≠ j. More 44

specifically, since equations (4.38) and (4.39) indicate that the error at epoch j cannot be a strict local maximum, the error at epoch j must be equal to at least one of its neighboring error values. This results in the equality of either equation 4.38 or 4.39. Both of these equalities give rise to the same constraint given below. i

kˆ  [ \  0 … 4.42 ^r

The filter that minimizes the error at epoch j, given that epoch j is a local (if not a global) maximum can be determined by extending the Lagrange multiplier approach using this third constraint. xN∆<=

ˆ

{sr {k  [ %  0 … 4.43 {\ {\ ^=

For coefficient αi where i < j, the G3 constraint does not impact the structure of the Lagrange multipliers, and the coefficient are same as in equation 4.27. All subsequent coefficients (i > j) can be shown to be identical. Since they must also sum up to zero (by constraint G3), all coefficients for i > j are zero. Hence the length of the filter can be shortened from L+1 to j as coefficients of value zero have no effect on the filter design. This implies that the optimal filter (if it exists!) will have its maximum transient error at the last instant of its transient phase and will be equal to the steady-state error. The constraint G3 only ensures that the error at epoch j is a local maximum. However to satisfy the optimization problem in equation 4.35, Ecuc,j is required to be the global maximum. It should also be observed that G3 constraints the local maximum to be at the end of the transient (i.e. at epoch j-1). For a filter with all positive coefficients, error increases monotonically as more coefficients are exposed to the ionosphere anomaly. Hence, the global maximum error occurs only at the last epoch in the transient. Thus, if j is less than L*+1, the largest length for which the filter coefficients are all positive, the optimization process is complete. The result of the optimization is a triangular shaped filter, with filter coefficients obeying a linear relationship over a filter length j given by equation 4.27 (with Lagrange multiplier constants given by equation 4.31 and length given by 4.32). This triangular shaped filter minimizes the maximum transient from the modeled ionosphere anomaly, and it is hence designated as the Maximum Transient Optimized Linear (MTOL) filter shown in the third part of figure 4.5. The comparison between the coefficients of the Hatch Filter and the MTOL filter is given in figure 4.6.

45

0.01 Hatch filter MTOL filter

0.009 0.008

Coefficient value

0.007 0.006 0.005 0.004 0.003 0.002 0.001 0

0

100

200

300 400 Coefficient index, i

500

600

700

Figure 4.6: Coefficients of the MTOL filter and the Hatch Filter

4.4 Maximum Transient Optimized Linear (MTOL) filter compared to Hatch Filter 4.4.1 Theoretical Comparison To evaluate how the MTOL filter compares to the Hatch Filter, it is instructive to compare the maximum transient error for each of the filters (for a fixed level of noise attenuation). As discussed in the last section the MTOL filter is given equations 4.27, 4.31 and 4.32. As the output noise characteristics Γ tends to be very small in practical operations, the maximum transient error for the MTOL can be evaluated by considering the small Γ limit. In this limit, the filter coefficients are the following. lim \ 

‰gQ

2 z :1  R ; … 4.44 R x x

Using equation 4.32 the filter order L* can be determined to be 4/3 Γ in the low Γ limit. For a digital filter, this length must be rounded down in practice. Using these values of L* and alpha, the divergence bias error can also be computed as below. lim S,iR 

‰gQ

4 N∆ … 4.45 9l

Applying a similar approximation to the Hatch Filter, using equations 4.15 and 4.21, the steady state divergence error of Hatch Filter can be shown to be as follows. lim SW 

‰gQ

1 N∆ … 4.46 2l 46

Hence taking the ratios of equation 4.45 and 4.46, a reduction of 11.11% is expected for the divergence error of the MTOL filter as compared to the Hatch Filter. Figure 4.7 pictorially represents this, showing equal error reduction for all noise attenuation levels. lim 

‰gQ

S,iR 8   … 4.47 SW 9

Although this result was derived for the limiting case of infinitesimally small output noise (Γ→0), the result is nonetheless accurate for finite output noise values in the range relevant for GPS applications. It is also important to remember the assumption of white noise at this point. The parameter Γ describes noise attenuation given that the input is a white Gaussian signal. Because actual GPS data features noise inputs that are correlated in time, the performance of the MTOL filter should be evaluated on actual GPS data sets.

Steady-State Error, normalized by d∆ t

50 Hatch filter MTOL filter

45 40 35 30 25 20 15 10 5 0.01

0.02

0.03

0.04

0.05

0.06

0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 Ratio of Output to Input Variance, Γ

0.07

Figure 4.7: Maximum Divergence bias Error as a Function of White-Gaussian Noise Rejection for Hatch and MTOL Filter

4.4.2 Experimental Comparison Analysis of GPS data indicates that the MTOL filter outperforms the Hatch Filter for a wide range of filter parameters. A theoretical comparison of the MTOL filter to the Hatch Filter predicted an 11% reduction of divergence bias error at equivalent noise level. Experimental results indicated improvements from 0% to 20% which was dependant on the dominant frequencies in noise. The MTOL filter was tested on GPS data. The GPS data was obtained from the Continuously Operated Reference Station (CORS) database observed at Oakland, CA reference station designated as ZOA1. A random day was chosen, which in this case was August 28, 2006 and the GPS data received from all the visible satellites during the first hour of the day was analyzed. Taking the advantages of the superposition principle for linear filters, the nominal noise performance and the ionospheric anomaly 47

response of the MTOL filter and the Hatch Filter were analyzed separately. For each data set both the filter lengths were experimentally set to achieve similar noise performance (i.e. the same value of Γ) and then the divergence bias error of the filters for the given ionospheric model was calculated analytically using equations 4.14 and 4.33. Analysis of this GPS data confirms that the optimized filter outperforms the conventional Hatch Filter for typical applications (i.e. for cases with a Hatch time constant of 100 s or more). For shorter smoothing times (with a Hatch time constant of 50 s), the optimized filter does not always provide an advantage. The strong performance of the optimized filter at moderate and long smoothing times, and its weaker performance at shorter smoothing times, relates directly to time correlation in GPS pseudorange data.

Satellite 8

Satellite 9 180 Hatch filter MTOL filter

160

Steady-State Error, normalized by d∆ t

Steady-State Error, normalized by d∆ t

180

140 120 100 80 60 40 Smoothing time=100sec for Hatch filter

20 0 0.04

0.05

0.06

0.07 0.08 0.09 0.1 0.11 0.12 Ratio of noise output to input variance

0.13

Hatch filter MTOL filter

160 140 120 100 80 60 40 Smoothing time=100sec for Hatch filter

20 0 0.01

0.14

0.02

Satellite 11

0.07

250

Hatch filter MTOL filter

160

Steady-State Error, normalized by d∆ t

Steady-State Error, normalized by d∆ t

0.06

Satellite 17

180

140 120 100 80 60 40 Smoothing time=100sec for Hatch filter

20 0 0.04

0.03 0.04 0.05 Ratio of noise output to input variance

0.05

0.06

0.07 0.08 0.09 0.1 0.11 0.12 Ratio of noise output to input variance

0.13

Hatch filter MTOL filter 200

150

100

50 Smoothing time=100sec for Hatch filter 0 0.06

0.14

48

0.07

0.08

0.09 0.1 0.11 0.12 0.13 0.14 Ratio of noise output to input variance

0.15

0.16

Satellite 24

Satellite 26 200 Hatch filter MTOL filter

140 120 100 80 60 40 Smoothing time=100sec for Hatch filter

20 0 0.06

0.07

0.08

0.09 0.1 0.11 0.12 0.13 0.14 Ratio of noise output to input variance

0.15

Steady-State Error, normalized by d∆ t

Steady-State Error, normalized by d∆ t

160

160 140 120 100 80 60 40 Smoothing time=100sec for Hatch filter

20 0 0.02

0.16

Hatch filter MTOL filter

180

0.03

0.04

Satellite 27

0.11

0.12

150 Hatch filter MTOL filter

180 160 140 120 100 80 60 40

Smoothing time=100sec for Hatch filter

20

0.03

0.04 0.05 0.06 Ratio of noise output to input variance

0.07

Steady-State Error, normalized by d∆t

Steady-State Error, normalized by d∆ t

0.1

Satellite 28

200

0 0.02

0.05 0.06 0.07 0.08 0.09 Ratio of noise output to input variance

Hatch filter MTOL filter

100

50

Smoothing time=100sec for Hatch filter 0 0.04

0.08

0.05

0.06

0.07 0.08 0.09 0.1 0.11 0.12 Ratio of noise output to input variance

0.13

0.14

Figure 4.8: Maximum Ionosphere Anomaly Tracking Error as a Function of Actual (Correlated) Noise Rejection for Hatch and MTOL Filters

GPS data from all the visible eight satellites were tested for the ionospheric storm anomaly. The tradeoff between output noise attenuation under nominal conditions and bias minimization under anomalous ionosphere conditions is similar whether analyzing correlated GPS noise or a white-noise model. Inspite of these similarities, noticed between theoretical performance (figure 4.7) and experimental performance (figure 4.8), the effect of time correlation of noise is evident for satellites 8, 11 and 17 in figure 4.8. For most of the satellites the MTOL filter showed a positive improvement. However the improvement is more pronounced at lower Γ values (which imply longer filter lengths). As the noise parameter Γ was increased, the improvement of the MTOL filter bias error over that of Hatch Filter decreased. For a typical Hatch Filter (γ = 0.01), the improvement ranged from 0% to 20%.

49

To understand the effect of correlated noise on the two filters, it is necessary to study their power spectra. Figure 4.9 shows the normalized power spectrum of a 100 second long (γ = 0.01) Hatch Filter and a 264 seconds long (L*=264) MTOL filter that have the same output noise performance Γ for white noise. The figure clearly shows that the MTOL filter attenuates higher frequency noise better while the Hatch Filter attenuates lower frequency noise better. For this particular case the “cross-over” frequency occurs at 3⋅10-3 Hz. More generally the cross-over frequency depends on the length of the two filters. For a 150 second long Hatch Filter and an equivalent 398 seconds long MTOL filter, the cross-over frequency is at 2⋅10-3 Hz. By contrast if the Hatch Filter is shortened to 50 seconds for which the equivalent MTOL filter is 132 seconds long, then the cross-over frequency changes to 6⋅10-3 Hz.

0

10

MTOL filter Hatch filter

Cross-over frequency

-1

Amplitude

10

-2

10

-3

10

-4

10

-3

10

-2

10 Frequency (Hz)

-1

10

0

10

Figure 4.9: Power Spectrum of the transfer functions of the two filters

Hence whether the MTOL filter is better than the Hatch Filter is dependent on the dominant frequency present in the time-correlated noise in the GPS signal. If the dominant frequencies are below the crossover frequency then Hatch Filter performs better while in the situation when the dominant frequency of the GPS signal is above the cross-over frequency the MTOL filter does better.

50

On analyzing the power spectrum of the code signal for a characteristic satellite, the dominant frequencies were found be in the range of 10-3Hz to 10-1Hz. Figure 4.10 shows the power spectrum of the GPS code signal from Satellite PRN 9 on 28th August 2006 during the first hour, as observed from ZOA1. Although this power spectrum is for a particular satellite it is representative of general GPS noise. 0

10

-1

10

-2

Normalized Power

10

-3

10

-4

10

-5

10

-6

10

-7

10 -4 10

-3

10

-2

10 Frequency (Hz)

-1

10

0

10

Figure 4.10: Power Spectra for GPS Pseudorange Input Data, Normalized to a Peak Value of One

The specific performance benefits of the optimized filter vary somewhat from one satellite to another. These variations depend on the noise spectra for individual satellites. To illustrate this effect, filters of a range of smoothing times were applied to all satellites in view during the span of the data set. For each case, output noise was assessed by comparing the filter output to the carrier signal plus the fifth-order curve fit to the nuisance parameters, 2I-λN (and their resulting divergence). Accordingly, the filter output noise was essentially zero mean for each satellite. Figure 4.11 plots the ionosphere bias for each satellite as a function of noise attenuation under nominal conditions. The vertical (ionosphere bias) axis is plotted as a ratio of the bias for Hatch and optimized filters with the same level of GPS noise attenuation. Thus values below unity (solid horizontal line) indicate that the optimized filter outperforms the Hatch Filter. The white-Gaussian noise (WGN) prediction, that the bias ratio is 8/9, is shown as a dashed line for reference. For moderate levels of noise attenuation (0.15 or better) the optimized filter nearly always outperforms the Hatch Filter. For any particular level of noise rejection, there is a significant spread in the bias ratio around the predicted level of 8/9. In some cases, the optimized filter delivers almost no performance advantage (ratio of 1). In other cases the optimized filter reduces ionosphere biases by as much as 20% or more.

51

1.3

PRN 8 PRN 11 PRN 17 PRN 24 PRN 26 PRN 27 PRN 28 PRN 9

1.2

1.3 1.2 1.1 Unity WGN PRN 8 PRN 11 PRN 17 PRN 24 PRN 26 PRN 27 PRN 28 PRN 9

1 0.9 0.8 0.7

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

Γactual

Figure 4.11: Ratio of Anomalous Ionosphere Bias for MTOL and Hatch Filters. Results are Plotted as a Function of Output Noise Attenuation (for Filters of Different Lengths) Applied to Eight Satellites in View.

σ opt/σ Hatch for Matched Max. Iono Error

Max. Iono Error for Optimum Filter Normalized by Hatch Filter

1.4

1.1 1 0.9 0.8 0.7 0.6 0.5 0.4

0

50

100

150 200 250 Hatch Time Constant, γ-1 (s)

300

350

400

Figure 4.12: Ratio of Output Noise Standard Deviations for MTOL and Hatch Filters (with Same Maximum Ionosphere Bias) as a Function of Hatch Time Constant. Data Plotted for All 8 Satellites in View.

Another way of looking at the analysis is represented in figure 4.12. In this case the maximum divergence errors of both the filters are matched and the ratio of their noise (standard deviation) is studied. Ratios of less than unity indicate that the MTOL filter outperforms the Hatch Filter in noise characteristics. It is evident that above a Hatch time constant of 150 seconds the MTOL filter is almost always better than the Hatch Filter, while for Hatch time constant below 50 seconds the Hatch Filter is better half the time.

4.5 Summary To mitigate filter divergence that results from carrier-smoothing of code measurements during an anomalous ionosphere storm, this chapter discussed an alternative linear filter. An optimization process determined that the minimum-error transient was obtained for a filter with an impulse function triangular in shape. This form was shown, using the method of Lagrange Multipliers, to have the lowestpossible maximum transient error possible across the entire space of linear filters with unity DC gain and a specified level of white-Gaussian noise attenuation. The optimized filter can be shown analytically to reduce the ionosphere-induced bias by 11% relative to a Hatch Filter for the same level of white-Gaussian noise attenuation. The optimized filter does, however, require more computational effort than the Hatch Filter, since the optimized filter cannot be implemented recursively. To assess the performance of the optimized filter on correlated data, the filter was applied to smooth eight GPS datasets. For short filter smoothing times (50 s Hatch time constant), the optimized filter did not outperform the conventional Hatch Filter. For moderate and long smoothing times, however, the optimized filter provides a clear advantage over the conventional Hatch Filter. For instance, for the standard Hatch Filter smoothing time (100 s time constant), the optimized filter can match the nominal

noise attenuation of the Hatch Filter while reducing the magnitude of an anomalous ionosphere bias by 0 to 20%, or more.

4.6 Zero Steady state Error Linear (ZSEL) filter – extension of MTOL filter The ZSEL filter is an extension of the MTOL filter to get a steady state error of zero. It is important to note here that ZSEL is not an optimized filter with the least instantaneous transient error and zero steady state error. Rather it has the same structure of the coefficients as the MTOL filter, but with the length changed appropriately for zero steady state error. Since MTOL filter has the least instantaneous transient tracking error, it is obvious to expect that the maximum transient error for the ZSEL filter is more than that of the MTOL filter. The reason behind this is that in order to get a zero steady state error; some of the coefficients need to be negative. However as all the coefficients have to sum up to zero, some coefficient will have to be increased to compensate for the negative coefficients. Because of this effect the transient divergence error increases. To mathematically compute the length of the ZSEL filter, the steady state error, defined in equation 4.24, was set to zero. This resulted in equation 4.47, given below, which is used along with equations 4.16, 4.17, 4.26 and 4.30 to solve the value of L for the ZSEL filter. Using mathematical manipulation the value of L is given below. x  U€uu Š

4  3l  √l   16l  16 Œ … 4.47 2l

The ZSEL filter has output noise characteristics, Γ and for that particular value of Γ, the length of the filter is defined as above. For instance, ZSEL filter of length 797 seconds has similar white noise characteristics as a Hatch Filter of 100 seconds length. The transient response of the ZSEL filter, given below in figure 4.12, to the ionospheric event shown in the inset (400mm/Km), illustrates the higher instantaneous divergence error of ZSEL filter compared to the MTOL filter. Although this optimized linear filter has the desired characteristics of zero steady state error, it increases the maximum divergence bias error by 19%. The ZSEL filter brings to light the point that in order to reduce the steady state error, high transient errors needs to be allowed. The understanding of the ZSEL filter is useful in understanding importance of the NLDE filter that will be introduced in the next chapter.

53

10 45

Hatch Filter ZSEL filter

40

8

35

Ionosphere delay --->

30

Divergence bias error (m)

6 4

25 20 15 10 5

2

0 -5

0

500

1000

0

1500 2000 Time --->

2500

3000

3500

-2 -4 -6 -8 -10

0

500

1000

1500 2000 Time (s)

2500

3000

3500

Figure 4.13: Comparison of transients of Hatch Filter and ZSEL filter of similar noise output

54

Chapter 5 Nonlinear Divergence Elimination (NLDE) Filter Design 5.1 Motivation The MTOL filter significantly decreases transient ionosphere error as compared to the Hatch Filter given requirements for nominal filter noise performance. Accordingly, the MTOL filter offers significant potential to improve LAAS performance as assessed through standard LAAS availability modeling [17, 21, 22] where the maximum transient error is the criteria of performance. Standard LAAS availability modeling tools, however, neglect an important detail which is the duration the divergence error persists. Because LAAS specifications [24] require users to apply a Hatch Filter to GPS measurements, these simulations fundamentally assume that users are exposed to extended duration errors during an ionospheric anomaly. In fact, it is possible to construct filters for which the steady-state error induced by an ‘iono’ ramp is zero. This type of filter would limit the time that approaching aircraft are exposed to large errors and reduce the overall risk of an unsafe aircraft landing. Although current availability modeling tools only analyze maximum transient error, there would be a strong impetus to consider exposure time were a suitable filter designed with zero steady-state error. The problem with filters that introduce zero steady-state error is that there is typically a tradeoff between the size of the max transient error and the size of the steady-state error. Reduction in steadystate error tends to increase max transient error. This was illustrated in figure 4.5. It was proved in the last chapter that the linear filter that minimizes max transient error (MTOL filter) reaches this max error condition in steady state. On the other hand a linear filter that decreases steady-state error to zero (ZSEL filter) has the detrimental effect of increasing the maximum transient. For safe aircraft landing applications, the challenge is to identify a class of filters that can reduce steadystate error without degrading nominal noise or the maximum transient. This is theoretically not possible with a linear filter, as proved last chapter and will be illustrated with an example in the next section. Therefore to answer this challenge we must look to formulate nonlinear filters. In particular, this chapter will look at a nonlinear filtering method for ionosphere bias estimation and subtraction based on identification of the parameters for a piecewise-linear ionospheric error model. Analysis shows that, for a certain range of allowable output noise, this class of filter can achieve the desired goals, reducing steady-state error to zero for a fixed level of maximum transient error and nominal noise.

5.2 An example of the limitations of linear filters The optimization procedure described in chapter 4 minimized the maximum instantaneous error. The resulting MTOL filter had this error as long as the LAAS receiver was exposed to the ionospheric delay gradient. If the time for which the aircraft were exposed to this maximum instantaneous error were also 55

minimized, then the risk to LAAS would be significantly reduced. Ideally the steady-state error should be zero with the maximum transient error and nominal noise minimized. Satellite 9

Maximum transient divergence bias error (m)

14 Hatchfilter (theorically calculated) MTOL filter (optimized for same gamma and best transient ZSEL filter (optimized for same gamma and zero steady state error)

12

10

8

6

4

2

0 0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

0.11

Noise Characteristic, Γ

Figure 5.1: Comparison of Hatch Filter, MTOL filter, ZSEL filter

Using this thought process; the MTOL filter was extended to drive the steady-state error to zero. The MTOL filter length was optimized for zero transient to get the ZSEL filter. On testing the ZSEL filter on GPS data with an ionospheric gradient of 400 mm/km superimposed, the maximum errors were found to increase at fixed noise levels. This observation is shown in figure 5.1 where the graph of Γ vs. divergence bias error shifts upward. This figure is based on archived GPS data from ZOA1 on 28th August 2006 during the first hour of the day of satellite 9. This means that the ZSEL filter partially achieved success in reducing the steady-state error to zero, but at the cost of increasing the maximum transient error. In the worst-case scenario the aircraft might reach its decision distance of 5 km from the airport exactly at the time of maximum transient. In this situation the ZSEL filter is worse than both the Hatch Filter and the MTOL filter. Therefore to improve the performance of the MTOL filter the time of exposure to the divergence bias error needs to be reduced without escalating the maximum instantaneous error as compared to the MTOL filter. The non-linear filter should be designed to achieve a zero steady-state divergence bias error subjected to three constraints. Firstly the proposed filter should not worsen availability during the transient phase as compared to the Hatch Filter. Secondly the filter should have comparable noise characteristics to the Hatch Filter. Thirdly the filter should have as short a settling time as possible to minimize the time of exposure to the divergence effects.

5.3 Non Linear Divergence Elimination (NLDE) algorithm One of the reasons of the poor performance of linear filters is that that they inherently treat a ramp as a straight line and ignore the starting time of the ionospheric ramp. This is not an accurate representation 56

of the ionospheric model in chapter 2 which is actually piecewise linear. An alternative approach is to estimate the starting time of the ionospheric anomaly along with the slope of the ramp, which is inherently a non-linear process. A novel idea of estimating the starting point of an ionospheric anomaly forms the core structure of the NLDE algorithm. The concept allows the filter to smooth the noise occurring before and after the ionospheric anomaly independently. This is achieved by varying the length of the lines fitted to the ionospheric data to exhaustively search for the event. The NLDE algorithm works on a number of ionosphere data points, P, buffered over a certain time period, and fits a piecewise linear curve to the data. This piecewise linear curve consists of two straight lines of variable lengths as shown in the figure 5.2. Although real ionospheric events do not necessarily have a sharply defined ‘start time’, the twoslope approach still provided a more accurate model of the ionosphere as shown in figure 5.3 than a one-slope approach. Both the straight lines are fitted using the least square method. The intersection point of these two lines is the start of the ionospheric event. The second of the two lines represents the ionospheric ramp and hence its slope should give an accurate estimate of the divergence rate. Knowing the divergence rate and the filter length, the bias error of the Hatch Filter can be calculated and hence it can be subtracted out. Best fit - Actual transition time at 200s 4 Ionospheric delay Fitted model

3

Ionospheric delay (m)

2

1

0 X: 206 Y: -0.1189

-1

-2

-3 0

50

100

150 Time (s)

200

250

Figure 5.2: Example of the curve fit on ionospheric delay data

57

300

Figure 5.3: a) Instantaneous ionopheric gradient change modelled by 1 slope and 2 slope method b) Gradual ionopheric gradient change modelled by 1 slope and 2 slope method

In operation, the NLDE generates a two-slope estimate of the ionosphere divergence error and subtracts this estimate from the user's smoothed pseudorange estimate, as illustrated in the block diagram of Figure 5.4. The total ionospheric delay error Iρ is calculated from the code and the carrier measurements. Iρ is the input to the nonlinear divergence estimator subroutine. The non-linear divergence estimator block estimates the divergence bias error from equation 1.6 and adds it to the Hatch Filter to compensate for the tracking error of the Hatch Filter. The nonlinear divergence estimator does an exhaustive search of the best fit to a buffered Iρ data points. Its block diagram is given in figure 5.5. A particular curve, comprising two straight lines, is the best fit to the ionospheric delay data when the sum of the distances of all the data points to this curve is least. Hence the first norm residual r is calculated as the sum of all the distances of the data points to the two slope model at the same epoch. Once the set of slopes (D1 and D2) and residuals (R) corresponding to set of transition times (S) is determined, the element with the least residual error r is the best fit. M  29  1

N ∆ … 5.1 N

58

This algorithm calculates the slope of the ionospheric delay which is half the divergence, d, defined in the ionospheric model. In other words it is a nonlinear w way of differentiating the Iρ term in equation (5.1) which is a restatement of equation (1.6) done here for convenience. Differentiation ifferentiation leads to the amplification of the higher frequency components present in the signal. To reduce the noise, the divergence estimate is passed through an Infinite Impulse Response (IIR) smoothing filter. However the advantages of the NLDE algorithm come at a price of computational time. While the processing time for the linear filters are in micro seconds, the processing time for the NLDE filter is in milliseconds, both being calculated on a 1.6 GHz dual core laptop processor. However since the GPS data used for analysis has a frequency of only 1 Hz, the filter can be realized in real time. The filter also needs memory to bufferr the GPS data. This is another aspect that is different from the Hatch Filter. The Hatch Filter did not need to store any data. However as memory on board is pretty common and cheap now a days, this new aspect pose no major problem in the implementation o off the NLDE filter.

Figure 5.4: Block diagram for NLDE filter

Figure 5.5: Block diagram of the nonlinear divergence estimator in the NLDE filter

S ≡ Set of transition times from 1 to P-L where P is the number of buffered epochs and L is smallest length of the lines s ≡ Member of S d1, d2 ≡ Slope on either side of the transition time D1 ≡ Set of d1 corresponding to each transition time s D2 ≡ Set of d2 corresponding to each transition time s r ≡ Residual error of the least square fit of the two slopes from Iρ R ≡ Set of all r corresponding to each transition time s

59

The ionosphere affects the code measurements in exactly the opposite way it affects the carrier measurements. Hence Iρ = - Iϕ. Solving for Iρ from equations (1.3) and (1.4) we get the expression for Iρ.   # !  !$ %&      … 5.2 2 2 2 We neglect the troposphere effect here as it is two orders smaller that the ionospheric effect. Using equation (5.1) and (5.2) the divergence bias error Bt is given below in equation (5.3). Note that λN is a constant. M 

29  1 N N# :  ; ∆ … 5.3 2 N N

GPS implementation uses discrete sampling whose frequency in this case is 1Hz. Thus the divergence bias error at any time t is calculated by equation (5.4) M 

29  1   <=   # #<=  … 5.4 2

The divergence bias error Bt calculated by equation (5.4) is at least twice as noisy as pseudorange data due to the summing of two noise code measurements. Hence instead of directly differentiating (ρ – ϕ/2, it is fitted by two straight lines as shown in figure 5.2. Each of the straight lines is determined by the least squares method as shown below. Let the event be at time ts.

Since GPS noise is not white, the slope of very small straight lines for the two slope model will be dominated by multipath rather than ionospheric delay. Therefore it is important to introduce a concept of minimum length, L, of each line to reduce the susceptibility of the slopes to multipath. As the slope of the first line is unimportant in the calculation of bias Bt as will be seen ahead, the minimum length constraint for the first line is relaxed. Hence the transition time s varies from 1 to P-L. Let the equation of the line fitting the (ρ – ϕ) for the first s seconds be y  dxe — 1 – N ’ ”–  “ – –[  • =

<=



[  š ™ = ™  ™ [   ™ ˜ =



— [ š  – ™ – = ™ … 5.5 – ™ –[   ™ • = ˜

The slope d and offset e are solved by equation (5.5) where xi is the time and yi is the value of (ρ – ϕ)/2 at that time. The second slope is calculated in a similar way using the (ρ – ϕ)/2 data from time = s+1 to P-L, which is the number of epochs the data buffered. After the two slopes are determined the fitted curve Ct is found. A fitness residual r is calculates as in equation (5.6). ž

  [ ›œ ^=

  #   T › … 5.6 2 60

The one-norm form residual was use to estimate the best fit. The transition time s is varied from 1 to P-L and r is calculated for all the curves. L is the minimum length assigned to each line to assure some amount of smoothing while estimating the rate of change of ionospheric delay. The curve with the least r is the best fit. The slope of the second line in the best fitted curve is the best estimate of

¡¢ ¡

from which the divergence bias error Bt can be calculated using equation 5.1. To reduce the noise being added to the Hatch Filter Bt is filtered by a smoothing filter F which in this case is a first order IIR filter. Since the divergence error is a lag, it is added to the Hatch Filter according to the equation 5.7. It should be noted here that ρhf,s is the smoothed range for the Hatch Filter while ρNLDE,s is the smoothed range from the NLDE filter with the divergence bias correction. ži§¨,  

1  1   :1  ; Y©ª, <=  #   # <= Z  M, … 5.7 9 9

5.4 Comparison of the NLDE and Hatch Filters

In comparison to a Hatch Filter of an equivalent length M, the NLDE filter produces an output with a substantially smaller (zero) steady-state divergence error in exchange for an increased random noise level under nominal conditions. Figure 5.6 gives a better understanding of the NLDE filter. It compares the divergence bias error of the Hatch Filter to that of the NLDE filter. In this particular case the length of the Hatch Filter is set to 70 seconds. The buffering time in the NLDE filter is set as 300 seconds. The minimum length L of each line with the curve fitting occurring in the NLDE filter is set as 60 seconds. The length of the Hatch Filter within the NLDE filter is also set as 70 seconds. The length of the IIR first order filter with the NLDE filter is set to be 200 seconds. The choice of parameters is explained in the next section. An ionosphere profile with an ionospheric gradient of 400 mm/km as shown in the inset of figure 5.6, which is rough estimate of the largest threat recorded as yet [17, 22, 23, 25], is the input to both these filters. While the Hatch Filter suffers from a divergence bias error of 5.56m, the NLDE filter Hatch Filter reduces the divergence bias error to zero in steady state while having the maximum error of only 2.55m.

61

3 NLDE filter Hatch Filter

2

0 -1 50

-2

40

Ionosphere error ---->

Divergence bias error (m)

1

-3 -4

30

20

10

0

-5

-10

-6

0

0

200

200

400

600

400

800

1000 1200 Time ---->

600

1400

1600

800

1800

2000

1000 1200 Time (s)

1400

1600

1800

2000

Figure 5.6: Response comparison of NLDE filter and Hatch Filter of same length to an ionospheric event 3 Input noise Noise output from NLDE filter Noise output from Hatch Filter

2

Noise (m)

1

0

-1

-2

-3

0

500

1000

1500 2000 Time (s)

2500

3000

3500

Figure 5.7: Noise characteristic comparison of NLDE filter and Hatch Filter of same length

62

The impact of the NLDE on filter accuracy can be assessed through the analysis of GPS data. Figure 5.6 shows the noise performance of the NLDE filter on nominal GPS data for CORS of a high elevation satellite and compares it with that of the Hatch Filter. To get the noise on each signal it is necessary to know the correct range. The exact range is not readily available, but one can estimate it subject to an unobservable bias, which may be set to zero without loss of generality for the purposes of error analysis. The unobservable bias is due to the integer ambiguity of the carrier measurements. By subtracting equation (1.4) from (1.3) one gets a (2I – λN) factor. The nominal ionospheric delay change very slowly during the day [3, 20] and hence can approximated by a polynomial curve. Also as N is a constant, a reasonably accurate smoothing of the (2I – λN) can be obtained by fitting a fifth order polynomial through it. Now taking advantage of the fact that the carrier measurement is much less noisy than the code, add the smoothed (2I – λN) to equation (2) to get a smoothed reference code. The noise is estimated by subtracting off this reference code from the smoothed code measurement from the Hatch Filter and the NLDE filter. In this particular case the Hatch Filter and the NLDE filter parameters were the same as in the previous section. The 1σ value of the input noise was 0.66m while that of the output from the Hatch Filter is 0.24m and 0.34m for the NLDE filter. That the NLDE noise is significantly higher than the Hatch Filter noise should not come as a surprise because the correction term that the NLDE filter introduces some noise which adds on to the noise of the Hatch Filter. It is noticeable that in this case the noise is dominated by multipath. The NLDE filter is affected by the trends of multipath more severely than the Hatch Filter. Multipath is generally negligibly small in carrier phase. Hence in the input of raw ionospheric delay Iρ, the multipath trends also creep in after being scaled by a factor of 0.5. As a result the NLDE also estimates the trends of multipath while calculating the rate of change of Iρ. Therefore the NLDE filter noise out is dominated by the multipath trends than the Hatch Filter, which is observable in figure 5.7.

63

a) The power spectrum of input noise

0

b) Power spectrum of noise from Hatch Filter

0

10

10

-1

10

-1

Normalised power

-2

10

-2

10

-3

10

-3

10

-4

10

-4

10 -4 10

-5

-3

10

-2

-1

10 Frequency (Hz)

10

0

10

10

-4

-3

10

10

-2

10 Frequency (Hz)

-1

10

c) Power spectrum of noise from NLDE filter

0

10

-1

10

Normalised Power

Normalised Power

10

-2

10

-3

10

-4

10

-5

10

-4

10

-3

10

-2

10 Frequency (Hz)

-1

10

0

10

Figure 5.8: Power spectrum of noise from (a) code measurement input (b) Hatch Filter (c) NLDE filter

The power spectrum of the output noise of the Hatch Filter and the NLDE filter show that they more or less attenuate similar frequency noise. Figure 5.8(a) shows that the input noise has more or less an uniform power of noise through frequency range of 10-3Hz to 10-2Hz. Figure 5.8(b) is the power spectrum noise output from the Hatch Filter which shows that the Hatch Filter attenuates all frequencies above 10-2Hz. Figure 5.8(c) is the power spectrum of the noise output from the NLDE filter which is very similar to that of the Hatch Filter. Although for nonlinear systems deduction from a particular power spectrum is valid only for that particular signal (as superposition does not hold) similar tests on other data also resulted in comparable Hatch and NLDE power spectrum. Visual inspection

64

0

10

suggested that the periods of oscillations of the errors were similar but the amplitude of the oscillations for the NLDE was bigger.

3 NLDE filter Hatch Filter of same length as NLDE filter Hatch filter of some noise output as NLDE filter

2

Divergence Bias Error (m)

1 0 -1 -2 -3 -4 -5 -6

0

200

400

600

800

1000 1200 Time (s)

1400

1600

1800

2000

Figure 5.9: Response of NLDE filter to pseudorange with ionospheric event

It is observed that for the given data set NLDE filter gives lesser maximum instantaneous divergence bias error than that of a Hatch Filter with same noise output. Since the NLDE filter is nonlinear in nature, superposition cannot be used. Hence it is important to check NLDE performance taking the noise and the ionospheric event together. Although at similar length the NLDE maximum divergence bias error is about 50% less than that of the Hatch Filter, the output noise of NLDE filter is more than that of the Hatch Filter. Therefore NLDE filter should be compared with two Hatch Filters, one of same length and the other of similar noise output characteristics. The 1σ noise level of the NLDE filter output is 0.34m which was match by a Hatch Filter of 36 seconds length. Figure 5.8 compares the output of the three filters. The parameters of the NLDE filter are kept the same as they were stated at the beginning of this section. Although for this data set figure (5.9) showed that NLDE filter outperforms the Hatch Filter, it was noted that the maximum transient error for NLDE filter for the ‘iono’ profile alone was 2.55m which increase to 2.7m when GPS noise was added. Although for the data with the GPS noise, the transient was longer. Hence this suggests that the performance of the NLDE filter is dependent on the GPS noise characteristics. Due to the nonlinear nature of the NLDE filter, the maximum transient of the NLDE from the noise free experiment cannot be assumed to be an average of maximum transients got from “noisy 65

data” experiments. Therefore it is essential to obtain a statistical average of the maximum transient error of NLDE filter.

5.5 Choice of parameters To optimize the performance of the NLDE filter, we evaluated its performance for a range of parameter values. The parameters that control the performance of the NLDE are the buffer time N, minimum line length L, Hatch Filter length M and the IIR filter length F. Although it was observed that increasing the buffer time N, the noise in the system reduces, as a side effect, the computation time linearly increases with the increase of N. Offline experiments showed that for N=500, the computation time is almost a second per data point. Since, in general, electronics on an aircrafts is slower that common laptop, it is expected that computation time will further increase for real applications. Also it should be noted that LAAS usually operates as 2 Hz, gave 0.5 seconds for each computation. Therefore buffering 500 data point is impractical for aerospace applications at this point. Computationally more efficient algorithms may be used in future to gain further noise reduction. For the current analysis the buffering time N was kept at 300 data points. The two parameters that are directly used in the nonlinear divergence estimator block diagram (figure 5.5) are the Hatch Filter length and the minimum length L. Due to nonlinearity in the NLDE filter we need to analyze data from many satellites to statistically formulate the trends.

Dependance of Γ on Hatch Filter length keeping L at 60s Noise output to into ratio, Γ

a)

0.3 PRN 8

0.25

PRN 9

0.2

PRN 11

0.15

PRN 17

0.1

PRN 24

0.05

PRN 26

0 0

20

40

60

80

Hatch Filter length, M (s)

66

100

PRN 27 PRN 28

Dependence of maximum transient error E on L keeping Hatch Filter length M at 60s Maximum Transient Error, E (m)

b) 5

PRN 9

3

PRN 11

2 1

PRN 17

0

PRN 24 0

20

40

60

80

100

PRN 26 PRN 27

Minimum Line Length L (s)

Dependence of maximum transient error E on Hatch Filter length keeping L at 60s

c) Maximum Transient Error, E (m)

PRN 8

4

5 PRN 8

4

PRN 9

3

PRN 11

2 1

PRN 17

0

PRN 24 0

20

40

60

80

Hatch Filter length, M (s)

100

PRN 26 PRN 27

Figure 5.10: Dependence of Noise output ratio and Maximum transient error on Hatch Filter length M and minimum length L

The noise output ratio Γ was observed to decrease for increasing values of L for all satellites and all Hatch Filters. Γ was also observed to decrease slowly with the increase of Hatch Filter length. Thus from noise performance point of view higher values of minimum line length L should be increased as should be Hatch Filter length, M. The maximum transient divergence error E was studied to find its dependence on L and M. It was found to increase with the increase in both L and M. As low errors are desirable, from the maximum transient error point of view, lower values of L and M are desirable. The parameters chosen for analysis in this paper are given in table 5.1. Values were selected empirically to provide reasonable performance. Further performance enhancement may be possible through a formal parameter optimization.

67

No. of buffered points (P epochs)

Minimum length (L epochs)

Hatch Filter length (M epochs)

Smoothing filter length (F epochs)

300 300 300 300 300 300 300 300 300

60 60 60 90 90 90 120 120 120

70 70 70 90 90 90 120 120 120

200 400 800 200 400 800 200 400 800

Table 5.2: Analyzed parameter space

5.6 Experimental comparison of NLDE filter to other filters The maximum transient divergence bias error of the NLDE filter was found to be less than that of the Hatch Filter for similar noise performance. The maximum transient error of the NLDE filter was statistically calculated since due to nonlinearity the transient error calculated on noiseless ionospheric storm model can be different from that observed with noisy data. However preliminary analysis in the last section showed only minor deviation from predictions using superposition. Hence to calculate the maximum transient error the ionospheric storm was placed at different random starting points on a GPS data set and the maximum divergence error in a window of ten epoch centered about the point of the theoretical “noiseless” occurrence of maximum transient was recorded. This window accommodated for small shifts, if any, in the occurrence of maximum transient divergence error due to the nonlinearity of the filter. The average of all the divergence errors, each corresponding to a different random starting point, was considered to be a good estimate of the maximum transient error of the NLDE filter. PRN 9

PRN 8 10

Hatch Filter NLDE Filter

9

Maximum transient divergence bias error (m)

Maximum transient divergence bias error (m)

10

8 7 6 5 4 3 2

Smoothing time = 100s for Hatch Filter

1 0 0.1

0.15

0.2

0.25

0.3

0.35

68

Hatch Filter NLDE Filter

9 8 7 6 5 4 3 2

Smoothing time=100s for Hatch Filter

1

0 0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

0.11

PRN 11

PRN 17 Hatch Filter NLDE Filter

9

Maximum transient divergence bias error (m)

Maximum transient divergence bias error (m)

10

8 7 6 5 4 3 2 Smoothing time=100s for Hatch Filter

1 0 0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

Hatch Filter NLDE Filter

9 8 7 6 5 4 3 Smoothing time=100s for Hatch Filter

2 1

0.22

0.07

0.08

Noise Characteristic, Γ

0.09

0.1

PRN 24 Maximum transient divergence bias error (m)

Maximum transient divergence bias error (m)

Hatch Filter NLDE Filter

9 8 7 6 5 4 3 Smoothing time=100s for Hatch Filter

1 0 0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

8 7 6 5 4 3

Smoothing time=100s for Hatch Filter

0.1

0.12

0.17

7 6 5 4 3 2

Smoothing time=100s for Hatch Filter

1

Maximum transient divergence bias error (m)

Maximum transient divergence bias error (m)

Hatch Filter NLDE Filter

0.08

0.16

Hatch Filter NLDE Filter

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

0.22

0.24

PRN 28

10

9

0.06

0.15

Noise Characteristic, Γ

PRN 27

0 0.04

0.14

8

0 0.04

0.22

10

1

0.13

9

Noise Characteristic, Γ

2

0.12

PRN 26

10

10

2

0.11

Noise Characteristic, Γ

0.14

0.16

0.18

0.2

Noise Characteristic, Γ

Hatch Filter NLDE Filter

9 8 7 6 5 4 3 2 1 0 0.05

Smoothing time=100s for Hatch Filter

0.1

0.15

0.2

0.25

Noise Characteristic, Γ

Figure 5.11: Comparison of NLDE and Hatch Filter

69

0.3

0.35

0.4

The NLDE filter was found to perform better than Hatch Filter on most of the GPS signals analyzed. To measure the performance of the two filters, all of them were tuned to give the same Γ value at nominal conditions. Γ is a measure of noise attenuation expressed as the ratio of the variance of output noise to that of input noise. The plots of maximum transient error versus Γ for the 8 satellites in figure 5.11 shows the relative performance of the Hatch Filter and the NLDE. However the NLDE filter does not beat the MTOL filter for the given set of parameters. This is well illustrated in figure 5.12 which shows the general trend. However, while the maximum transient errors of the Hatch Filter and the MTOL filter are also equal to their steady state errors, the NLDE filter transient errors decay to zero in steady state. It should be noted here that the errors for the Hatch Filter and MTOL filter were computed analytically by invoking superposition. However since the number of data points in the case of NLDE filter are finite, the statistically averaged maximum transient error is only an approximate. For better comparison the maximum transient error of the Hatch Filter was also confirmed statistically to check for validity of the analysis.

10 Hatch Filter MTOL Filter NLDE Filter

9 8 7 6 5 4 3 2 1 0 0.05

Smoothing time=100s for Hatch Filter

0.1

0.15

0.2

0.25

Figure 5.12: Comparison of Hatch Filter, optimized linear filter and NLDE filter done on PRN 17 data

Figure 5.11 and 5.12 clearly shows that NLDE filter is better than the Hatch Filter and comparable to the MTOL filter. Experimentation with various satellites at different times of the day has resulted in the maximum reduction of 29% of divergence bias error while in the worst case; NLDE gave a 14% increase in the divergence bias error. In fact in the best case scenario the NLDE filter was even better than the MTOL filter. The results are limited by the number of satellites examined and also by the errorbars on each data point. The reduction in maximum transient error is an added improvement given that the NLDE filter drives this maximum error to zero. The high performance of the NLDE filter came at computational cost. Due to the exhaustive search of an ionospheric anomaly at every instant, the computational time requirement is very high. The algorithm computation order was approximately O(P) i.e. the buffer time. For the offline tests, the time taken per 70

epoch per satellite was approximately 0.1 seconds for a 1.6GHz laptop. However several numerical strategies might be taken in future to ensure faster convergence in the NLDE block while finding the best curve fit. The settling time to attain zero steady-state divergence bias varied from 300 epochs to 900 epochs. The settling time was found to be dependent on the choice of parameters especially on F. The smoothing, filter being an IIR filter, produced exponential decay of the maximum transient error. Future works include experimenting with alternate smoothing filters with faster transient responses in place of the IIR filter. The NLDE filter design seemed to deviate only marginally from the predictions made using superposition. Although all the modeled ionospheric storms had a gradient of 400 mm/km, similar offline experiments were also conduction for ionospheric storm model with a gradient of 200 mm/km. The NLDE filter gave similar reduction of divergence bias for 200 mm/km.

71

Chapter 6 Conclusion GPS systems have been around for more than 30 years. They have proved to be a great position sensor because they use absolute positioning rather than dead reckoning (estimating current position based on previously estimated position, known speed and elapsed time). They are successfully used to navigate cars where road maps are available to aid the accuracy of the system. Commercial GPS systems have an accuracy of about 3 meters which is sufficient, along with routing algorithms, to provide reliable navigation. However for aircraft navigation, the safety concerns are much higher. In particular, the divergence bias error, which occurs during rare ionosphere anomalies, is the single greatest threat to the integrity of differentially aided GPS landing systems, such as the FAA’s LAAS. Two methods were proposed to minimize the impact of ionosphere divergence on differentially aided GPS aircraft landing. In the first approach a linear filter, named the MTOL filter, was designed to minimize the maximum transient error for a given level of noise attenuation. An optimization technique was used to derive this alternative to the standard Hatch filter. To keep the system robust the noise was assumed to be white and a Lagrange Multiplier optimization process was adopted to formulate a filter with the least instantaneous transient error. The resulting MTOL filter was triangular in shape, and the predicted improvement on the divergence bias error was calculated as 11%. The maximum transient error of the MTOL filter occurred at the final instant of its transient phase and was equal to the steady-state error. On conducting offline experiments on GPS data with the MTOL filter, the improvement in divergence bias error was found to vary from 0 to more than 20%. On further investigation of the properties of the MTOL filter, it was found that it outperformed the Hatch filter only when the frequency of the noise was above a certain “cut off” frequency. The cut-off frequency was found to vary with the amount of smoothing of receiver noise the filter was designed for. For the standard Hatch filter of length 100 seconds, the cut-off frequency was determined to be 3·10-3Hz. The power spectrum of typical GPS noise showed that the noise frequency was of the same order as the cut-off frequency and hence the MTOL performance was found to vary with different GPS signals. The MTOL filter produced significant improvement in most cases. Although the MTOL filter was more difficult to implement than the recursive Hatch filter, its computation requirement was found to be insignificant compared to modern computing technologies available. In the second approach a nonlinear filter, named the NLDE filter was designed, reduce the maximum transient error by the same amount as the MTOL filter (at similar noise performance) but also to reduce the duration of the maximum error by driving the steady-state error to zero. Although the MTOL filter was able to reduce the maximum transient divergence bias error, it had a nonzero steady state error. Current LAAS analysis techniques consider only the maximum transient error, not the duration over 72

which the error lasts. Reducing the time that users are exposed to the divergence error (by reducing the steady-state error) reduces the probability of a LAAS integrity failure. A non-linear approach was adopted, in which the ionospheric model was estimated with a piece-wise linear curve which had the capability of detecting the starting time of an ionosphere anomaly. The resulting NLDE filter was tested on GPS data and compared to the Hatch and MTOL filter. The NLDE filter was able to achieve a similar maximum transient error reduction as did the MTOL filter compared to Hatch filter, but the NLDE filter had an added advantage of having a zero steady-state error. In its working range it has settling time from 200 seconds to 1000 seconds, dependent on the amount of noise reduction required. Therefore in choosing its parameters, a trade-off between reducing the risk of LAAS integrity failure and allowing more noise in the system needs to be considered. The major limitation of the NLDE filter is that it cannot be tuned for arbitrary levels of code noise reduction under nominal conditions. The nonlinear filter achieves best performance when the noise level under nominal conditions is allowed to match that of a Hatch Filter with a 50 second time constant. The better performance of the NLDE filter came at the cost of computation time, although it was argued that given the operating frequency of LAAS is 2 Hz, the system would have enough time for the complex computations associated with the NLDE filter. One area of improvement required for the NLDE filter is its transient time. The NLDE filter is quite appropriate to medium ionospheric delay gradients like 100-200 mm/km which are difficult to detect reliably. To further reduce the risk posed by ionospheric storm, especially the higher anomalous spatial gradients (300-400 mm/km) where the NLDE filter is less effective, a monitor might be used in unison with the NLDE filter to detect and eliminate the threats posed by the high ionospheric delay gradients. A reasonable approach to build the monitor might be to remove the transients due to nominal divergence using a subset of the NLDE algorithm which would reduce the standard deviation of divergence statistics (code – carrier). With a good implementation of a monitor, the combined system of the NLDE filter and the monitor would be able to reduce the ionospheric threat to a large extent. There is a lot of scope for future work in this field. Fine tuning of the NLDE filter is possible. The first order recursive block in the NLDE filter can be replaced by a Butterworth filter or Chebyshev filter to improve its transient response. Also there is scope of analytical manipulation (using numerical methods) the NLDE filter parameters to reduce the computation time. One of the assumptions that went into this thesis was that the ionospheric anomaly had a sudden change. Removing this constraint will help build a more ‘real’ ionospheric model which will allow for NLDE filter to assume the ionosphere model as 3piecewise linear model which might improve the noise characteristics. Another approach to deal with the ionospheric threat is to remove the ionospheric delay completely from the system. Preliminary work was done to modify the Hatch Filter to completely remove the ionospheric delay from the range estimation. This resulted is increased noise and a transient error in the system. Although a longer smoothing time could be used to reduce the output noise, the major drawback of the system was that the initial ionospheric delay could not eliminated. One possibility of solving the initial ionospheric delay is to use measurements from more than one epoch and solve out the ionospheric delays assuming that the delays don’t change significantly in a couple of seconds.

73

References [1] “Navigation Integrity for Aircraft Precision Landing using the Global Positioning System”, Boris Steven Pervan – Ph.D. Dissertation, Stanford University, 1996 [2] “Ionospheric Correlation Analysis for WAAS: Quiet and Stormy”, Andrew Hansen, Juan Blanch, Todd Walter, Per Enge, Institute of Navigation's GPS conference, Salt Lake City, UT, September 2000 [3] “Global Positioning System: Signals Measurements and Performance”, Pratap Mishra, Per enge, Ganga-Jamuna Press, 2001 [4] “Land-Vehicle Navigation Using GPS”, Eric Abbott, David Powell, Proceedings of the IEEE, Vol. 87(1), Pg. 145 – 162. [5] “Interference Direction Finding for Aviation Applications of GPS”, Konstantin Gromov, Dennis Akos, Sam Pullen, Per Enge, Bradford Parkinson, Boris Pervan, Institute of Navigation's GPS conference, Nashville, TN, September 1999 [6] “System Overview, Recent Developments, and Future Outlook for WAAS and LAAS”, Sam Pullen, Todd Walter, and Per Enge, Tokyo University of Mercantile Marine GPS Symposium, Tokyo, Japan, November 11 - 13, 2002 [7] “High Integrity Carrier Phase Navigation for Future LAAS Using Multiple Civilian GPS Signals”, Jaewoo Jung, Institute of Navigation's GPS conference, Nashville, TN, September 1999 [8] “The Effects of Large Ionospheric Gradients on Single Frequency Airborne Smoothing Filters for WAAS and LAAS”, Todd Walter, Seebany Datta-Barua, Juan Blanch, Per Enge, Institute of Navigation's National Technical Meeting, San Diego CA, January 2004 [9] “Enhanced DifferentialGPS Carrier-Smoothed Code Processing Using Dual-Frequency Measurements”, Patrick Y. Hwang, Gary A. Mcgraw, and John R. Bader, Journal of The Institute of Navigation, Vol 46(2), Pg 127-138 [10] “Overbounding SBAS and GBAS Error Distributions with Excess-Mass Functions”, Jason Rife, Todd Walter, J. Blanch, GNSS 2004, Sydney, 6-8 December, 2004 [11] “Formulation of a Time-Varying Maximum Allowable Error for Ground-Based Augmentation Systems”, Jason Rife, R. Eric Phelts, Institute of Navigation's National Technical Meeting, Monterey CA, January 2006 [12] “Paired Overbounding and Application to GPS Augmentation”, Jason Rife, Sam Pullen, Boris Pervan and Per Enge, Position Location and Navigation Symposium, 2004. PLANS 2004, 26-29 April 2004, Pg 439-446.

74

[13] “Code-Carrier Divergence Monitoring for the GPS Local Area Augmentation System”, Dwarakanath V. Simili, Boris Pervan, Position, Location, And Navigation Symposium, 2006 IEEE/ION, 25-27 April 2006, pg 483-493. [14] “The Advantages of Local Monitoring and VHF Data Broadcast for SBAS”, Todd Walter, Sam Pullen, Jason Rife, Jiwon Seo, Per Enge, European Navigation Conference GNSS 2005, Munich, Germany, July 2005 [15] “Ionosphere Monitoring Methodology for Hybrid Dual-Frequency LAAS”, Hiroyuki Konno, Sam Pullen, Jason Rife, Per Enge, ION Institute of Navigation Global Navigation Satellite Systems Conference, Fort Worth, TX, September 2006 [16] United States Patent 5867411, “Kalman filter ionospheric delay estimator”, Rajendra Kumar [17] “LAAS Ionosphere Spatial Gradient Threat Model and Impact of LGF and Airborne Monitoring”, Ming Luo, Sam Pullen, Jed Dennis, Hiroyuki Konno, Gang Xie, Todd Walter, Per Enge, Seebany DattaBarua and Tom Dehel, ION GPS 2003 [18] “Symmetric Overbounding of Correlated Errors”, J. Rife and D. Gebre-Egziabher, NAVIGATION, Journal of the Institute of Navigation, in press. [19] “Using IRI for the computation of Ionospheric Corrections for Altimeter Data Analysis”, D. Bilitzn, C. Koblinsky, B. Beckley, S. Zia, R. Williamson, Advances in Space Research, Vol. 15(2), Feb 1995, Pg. 113119 [20] "Ionospheric Effects on GPS", J. Klobuchar, in Global Positioning System: Theory and Applications, vol. I, B. Parkinson and J. Spilker, Eds. Washington, D.C: AIAA, 1996, pp. 485-515. [21] “Ionosphere Threat to LAAS: Updated Model, User Impact, and Mitigations”, Ming Luo, Sam Pullen, Alexandru Ene, Di Qiu, Todd Walter, and Per Enge, Proceedings of ION GNSS 2004 [22] “LAAS Study of Slow-Moving Ionosphere Anomalies and Their Potential Impacts”, Ming Luo, Sam Pullen, Seebany Datta-Barua, Godwin Zhang, Todd Walter, and Per Enge, Proceedings of Institute of Navigation GNSS 2005 [23] “Using WAAS Ionospheric Data to Estimate LAAS Short Baseline Gradients”, Seebany Datta-Barua, Todd Walter, Sam Pullen, Ming Luo, Juan Blanch, and Per Enge, Proceedings of ION 2002 National Technical Meeting [24] “WAAS MOPS: Practical Examples”, Todd Walter, ION NTM '99, San Diego, CA, January 1999 [25] “Modeling the 20 November 2003 Ionosphere Storm with GRACEModeling the 20 November 2003 Ionosphere Storm with GRACE”, Seebany Datta-Barua, Todd Walter, Sam Pullen, Per Enge, ION Institute of Navigation Global Navigation Satellite Systems Conference, Fort Worth, TX, September 2007

75

Appendix Program for communication with Garmin GPS 18 USB The program was written in C using Visual Studio Package. The program first initializes the GPS device. For that it requires a handle to the USB port. For this it uses functions given in the header files which are windows.h, initguid.h, setupapi.h and winioctl.h . After getting the handle to the USB port the program sends the initialization packet which is {0, 0, 0, 5, 0, 0}. Once the GPS 18 USB is initialized, the program sends the command packets to turn on position record and receiver measurement record. As soon as this is done, the GPS sensor send two packets every second, one position record and the other receiver measurement record. These data packets come in the same format as the structure given in table 3.1. Packet ID 51 signifies position record while packet ID 52 signifies receiver position record. A pointer to the chunk of data bytes is stored in 12th byte of the USB packet. Since the data is really only a series of numbers, it is important to parse it properly to be able to interpret the data. The structure of the position record and the receiver measurement record was got from the specifications sheets available at Garmin’s website given in the previous section. After parsing, the data was stored in files opened in ‘append’ mode, noting the satellite PRN, time, pseudorange measurement and carrier phase measurement.

C Code #include #include #include #include #include #include #include

// You may need to explicitly link with setupapi.lib

DEFINE_GUID(GUID_DEVINTERFACE_GRMNUSB, 0x2c9c45c2L, 0x8e7d, 0x4c08, 0xa1, 0x2d, 0x81, 0x6b, 0xba, 0xe7, 0x22, 0xc0); #define IOCTL_ASYNC_IN CTL_CODE (FILE_DEVICE_UNKNOWN, 0x850, METHOD_BUFFERED, FILE_ANY_ACCESS) #define IOCTL_USB_PACKET_SIZE CTL_CODE (FILE_DEVICE_UNKNOWN, 0x851, METHOD_BUFFERED, FILE_ANY_ACCESS) #define MAX_BUFFER_SIZE 4096 #define ASYNC_DATA_SIZE 64 //------------------------------------------------------------------------HANDLE gHandle; DWORD gUSBPacketSize; //------------------------------------------------------------------------#pragma pack( push, 1) typedef struct { unsigned char mPacketType;

76

unsigned char mReserved1; unsigned short mReserved2; unsigned short mPacketId; unsigned short mReserved3; unsigned long mDataSize; unsigned short mData[1]; } Packet_t; #pragma pack( pop) #pragma pack( push, 1) typedef struct { unsigned long cycles; double pr; unsigned short phase; char slp_dtct; unsigned char snr_dbhz; unsigned char svid; char valid; } cpo_rcv_sv_data; #pragma pack( pop) #pragma pack( push, 1) typedef struct { double rcvr_tow; short rcvr_wn; cpo_rcv_sv_data sv[12]; } cpo_rcv_data; #pragma pack( pop) #pragma pack( push, 1) typedef struct { char svid; unsigned short snr; unsigned char elev; unsigned short azmth; unsigned char status; } cpo_sat_data; #pragma pack( pop) #pragma pack( push, 1) typedef struct { cpo_sat_data sat[12]; } cpo_all_sat_data; #pragma pack( pop) #pragma pack( push, 1) typedef struct { float alt; float epe; float eph; float epv; short fix; double gps_tow; double lat; double lon; float long_vel; float lat_vel;

77

float alt_vel; float msl_hght; short leap_sec; long grmn_days; } cpo_pvt_data; #pragma pack( pop) #pragma pack( push,1 ) typedef struct { short wn; BYTE junk1[2]; float toc; float toe; float af0; float af1; float af2; float ura; double e; double sqrta; double dn; double m0; double w; double omg0; double i0; float odot; float idot; float cus; float cuc; float cis; float cic; float crs; float crc; unsigned char iod; char svid; BYTE junk2[2]; } SDM_spc_eph_type; #pragma pack( pop) //------------------------------------------------------------------------void SendPacket(Packet_t aPacket) { DWORD theBytesToWrite = sizeof( aPacket ) - 1 + aPacket.mDataSize; DWORD theBytesReturned = 0; WriteFile(

gHandle, &aPacket, theBytesToWrite, &theBytesReturned, NULL );

// If the packet size was an exact multiple of the USB packet // size, we must make a final write call with no data if( theBytesToWrite % gUSBPacketSize == 0 ) { WriteFile( gHandle, 0, 0,

78

&theBytesReturned, NULL ); } } //------------------------------------------------------------------------// Gets a single packet. Since packets may come simultaneously through // asynchrous reads and normal (ReadFile) reads, a full implementation // may require a packet queue and multiple threads. Packet_t* GetPacket() { Packet_t* thePacket = 0; DWORD theBufferSize = 0; BYTE* theBuffer = 0; for( ; ; ) { // Read async data until the driver returns less than the // max async data size, which signifies the end of a packet BYTE theTempBuffer[ASYNC_DATA_SIZE]; BYTE* theNewBuffer = 0; DWORD theBytesReturned = 0; DeviceIoControl(

gHandle, IOCTL_ASYNC_IN, 0, 0, theTempBuffer, sizeof( theTempBuffer ), &theBytesReturned, NULL );

theBufferSize += ASYNC_DATA_SIZE; theNewBuffer = (BYTE*) malloc( theBufferSize ); memcpy( theNewBuffer, theBuffer, theBufferSize - ASYNC_DATA_SIZE ); memcpy( theNewBuffer + theBufferSize - ASYNC_DATA_SIZE, theTempBuffer, ASYNC_DATA_SIZE ); free( theBuffer ); theBuffer = theNewBuffer; if( theBytesReturned != ASYNC_DATA_SIZE ) { thePacket = (Packet_t*) theBuffer; break; } } // If this was a small "signal" packet, read a real // packet using ReadFile if( thePacket->mPacketType == 0 && thePacket->mPacketId == 2 ) { BYTE* theNewBuffer = (BYTE*) malloc( MAX_BUFFER_SIZE );

79

DWORD theBytesReturned = 0; free( thePacket ); printf("small packet\n\n"); // A full implementation would keep reading (and queueing) // packets until the driver returns a 0 size buffer. ReadFile( gHandle, theNewBuffer, MAX_BUFFER_SIZE, &theBytesReturned, NULL ); return (Packet_t*) theNewBuffer; } else { //printf("big packet\n\n"); return thePacket; } }

//------------------------------------------------------------------------void Initialize() { // Make all the necessary Windows calls to get a handle // to our USB device DWORD theBytesReturned = 0; PSP_INTERFACE_DEVICE_DETAIL_DATA theDevDetailData = 0; SP_DEVINFO_DATA theDevInfoData = { sizeof( SP_DEVINFO_DATA ) }; Packet_t theStartSessionPacket = { 0, 0, 0, 5, 0 , 0 }; Packet_t* thePacket = 0; HDEVINFO theDevInfo = SetupDiGetClassDevs( (GUID*) &GUID_DEVINTERFACE_GRMNUSB, NULL, NULL, DIGCF_PRESENT | DIGCF_INTERFACEDEVICE ); SP_DEVICE_INTERFACE_DATA theInterfaceData; theInterfaceData.cbSize = sizeof( theInterfaceData ); if( !SetupDiEnumDeviceInterfaces( theDevInfo, NULL, (GUID*) &GUID_DEVINTERFACE_GRMNUSB, 0, &theInterfaceData ) && GetLastError() == ERROR_NO_MORE_ITEMS ) { gHandle = 0;

80

return; } SetupDiGetDeviceInterfaceDetail( theDevInfo, &theInterfaceData, NULL, 0, &theBytesReturned, NULL ); theDevDetailData = (PSP_INTERFACE_DEVICE_DETAIL_DATA) malloc( theBytesReturned ); theDevDetailData->cbSize = sizeof( SP_INTERFACE_DEVICE_DETAIL_DATA ); SetupDiGetDeviceInterfaceDetail( theDevInfo, &theInterfaceData, theDevDetailData, theBytesReturned, NULL, &theDevInfoData ); gHandle = CreateFile( theDevDetailData->DevicePath, GENERIC_READ | GENERIC_WRITE, 0, NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL ); free( theDevDetailData ); // Get the USB packet size, which we need for sending packets DeviceIoControl( gHandle, IOCTL_USB_PACKET_SIZE, 0, 0, &gUSBPacketSize, sizeof( gUSBPacketSize ), &theBytesReturned, NULL ); // Tell the device that we are starting a session. SendPacket( theStartSessionPacket ); // Wait until the device is ready to the start the session for( ; ; ) { thePacket = GetPacket(); if( thePacket->mPacketType == 0 && thePacket->mPacketId == 6 ) { break; }

81

free( thePacket ); } free( thePacket ); } //-------------------------------------------------------------------------int _tmain(int argc, _TCHAR* argv[]) { Packet_t theProductDataPacket = { 20, 0, 0, 254, 0 , 0 }; Packet_t theCommandPacket1 = { 20, 0, 0, 10, 0 , 2, 49 }; Packet_t theCommandPacket2 = { 20, 0, 0, 10, 0 , 2, 110 }; Packet_t* thePacket = 0; cpo_pvt_data* pvt_data = 0; cpo_all_sat_data* sat_data = 0; cpo_rcv_data* rcv_data = 0; SDM_spc_eph_type* spc_data = 0; int i,j; short k=0; short* count; long seconds,sec,t1; char line_str[64] ={0}; FILE* file_handle = NULL; seconds = clock(); printf("%d\n",sizeof(long)); printf("%d\n",sizeof(Packet_t)); printf("%d\n",sizeof(cpo_pvt_data)); printf("%d\n",sizeof(cpo_all_sat_data)); printf("%d\n",sizeof(cpo_rcv_data)); printf("%d\n",sizeof(SDM_spc_eph_type)); Initialize(); if( gHandle == 0 ) { printf( "%s", "No device" ); return 0; } // Tell the device to send product data SendPacket( theProductDataPacket ); // Get the product data packet for( ; ; ) { thePacket = GetPacket(); if( thePacket->mPacketType == 20 && thePacket->mPacketId == 255 ) { break; } free( thePacket ); } // Print out the product description printf( "Version info : %s\n\n", (char*) &thePacket->mData[4] );

82

thePacket = GetPacket(); SendPacket( theCommandPacket1 ); SendPacket( theCommandPacket2 ); // Get the PVT data packet j=1; file_handle = fopen("Code.txt","w"); for(i=0 ;i<11000 ;i++ ) { thePacket = GetPacket(); sec = clock (); t1 = (sec-seconds); //printf( "Time :%ld\n", (sec-seconds) ); if (thePacket->mPacketId==114) { sat_data = &thePacket->mData; printf( "Sat pkt Sat3 SVID: %d\n\n", sat_data>sat[2].svid); } else if (thePacket->mPacketId==51) { pvt_data = &thePacket->mData; printf( "PVT latitute: %lf\n", pvt_data->lat); printf( "PVT longitude: %lf\n", pvt_data->lon); printf( "PVT altitude: %lf\n\n", pvt_data->alt); } else if (thePacket->mPacketId==52) { rcv_data = &thePacket->mData; printf( "Pseudorange : %lf\n", ((rcv_data->sv)+j)>pr); printf( "Carrier Cycles : %ld\n", ((rcv_data->sv)+j)>cycles); printf( "Phase

: %d\n", ((rcv_data-

>sv)+j)->phase); printf( "Validity

: %d\n", ((rcv_data->sv)+j)-

>valid); printf( "SVID

: %d\n", ((rcv_data-

printf( "Time

: %d\n\n", rcv_data-

>sv)+j)->svid); >rcvr_tow); sprintf(line_str,"%lf\t%ld\t%d\t%d\t%d\t%d\n",((rcv_data->sv)+j)>pr,((rcv_data->sv)+j)->cycles,((rcv_data->sv)+j)->phase,((rcv_data->sv)+j)>valid,((rcv_data->sv)+j)->svid,rcv_data->rcvr_tow); printf("%s\n",line_str); fwrite(line_str,1,sizeof(line_str),file_handle); } free( thePacket ); } fclose(file_handle); return 0; }

83

84

GPS Filtering to Minimize Ionosphere Divergence ...

System (GPS) can be used to design potentially cheaper and more versatile navigation .... Maximum Ionosphere Anomaly Tracking Error as a Function of Actual ..... does a very good job of mitigating multipath and receiver noise for users in motion. ...... [4] “Land-Vehicle Navigation Using GPS”, Eric Abbott, David Powell, ...

1MB Sizes 1 Downloads 148 Views

Recommend Documents

Simulating the Ionosphere - GitHub
Sep 30, 2009 - DEFINITION: Approximating measurements at intermediate scales/positions from scattered measurements. We have sparse measurements.

Application to Divergence Estimation Introduction ...
2. ( ) , where 1 and 2 are densities and is a smooth convex function. • Includes KL and Renyi- divergences. • Important in statistics, machine learning, and information theory. ▫ Convergence rates are unknown for most divergence

Cheap Tk102 Gps Track Locator Gps Mount Tracker Gps Bracket ...
Cheap Tk102 Gps Track Locator Gps Mount Tracker Gp ... ket For Dji Phantom 4 Quadcopter Free Shipping.pdf. Cheap Tk102 Gps Track Locator Gps Mount ...

Distributed Ionosphere Monitoring by Collaborating ...
This effect is a compelling advantage of the proposed method, but it will not ..... randomly distributed over a circular area of 50 km in radius is shown in figure Figure 4. In the specific ..... As an illustration of the reduction in complexity due

A Multi-Carrier Collaborative Solution to Minimize Connectivity-l.pdf ...
A MULTI-CARRIER COLLABORATIVE SOLUTION TO MINIMIZE. CONNECTIVITY-LOSS. A Thesis. presented to. the Faculty of California Polytechnic State University,. San Luis Obispo. In Partial Fulfillment. of the Requirements for the Degree. Master of Science in

Appendix to "Reconciling the divergence in ... - André Kurmann
Mar 23, 2017 - ... 2002 North American Industry Classification System (NAICS), the U.S. economy ... LPC covers the non-farm business sector of the U.S economy. ...... To construct a sample analog of S, we use the Newey-West estimate of S:.

Auto-Scaling to Minimize Cost and Meet Application ...
Motivation: real world application. 6. • Huge memory requirement: ... Powerful computing capability: hundreds or thousands or even more of iterations.

Scheduling to Minimize Gaps and Power Consumption
in execution. For both versions in a multiprocessor system, we develop a ...... 2qth place of extra-interval will be filled and the second unit in extra-interval will be ...

Google Message Filtering - PDFKUL.COM
ABOUT GOOGLE APPS. Google Apps is a suite of applications that includes Gmail, Google Calendar. (shared calendaring), Google Talk. (instant messaging and voice over IP),. Google Docs & Spreadsheets (online document hosting and collaboration),. Google

Download The Lean Farm: How to Minimize Waste ...
sound business practices apply whether you are producing cars or carrots. ... finding that incorporating the best new ideas from business into their farming.