An improved Hatch filter and its application in kinematic positioning with single-frequency GPS

Measurement 146 (2019) 868–878

Contents lists available at ScienceDirect

Measurement journal homepage: www.elsevier.com/locate/measurement

An improved Hatch filter and its application in kinematic positioning with single-frequency GPS Qinghua Zhang a,b,c,⇑, Zhengsheng Chen c, Fengjuan Rong a, Yang Cui d a

Army Engineering University of PLA, Nanjing, Jiangsu 210007, China Nanjing University of Science and Technology, Nanjing, Jiangsu 210094, China c State Key Laboratory of Geo-information Engineering, Xian 710054, China d Army Logistics University of PLA, Chongqing 401311, China b

a r t i c l e

i n f o

Article history: Received 19 February 2019 Received in revised form 5 July 2019 Accepted 6 July 2019 Available online 10 July 2019 Keywords: GNSS Kinematic positioning Hatch filtering Ionosphere Single frequency

a b s t r a c t The Hatch filter has been widely used in GNSS kinematic positioning and navigation. However, the classic carrier smoothing of code pseudorange observations often results in filter divergence and accuracy degradation in the positioning solution due to the influence of ionospheric delay. One solution is to reduce the impact of ionospheric delay by using the Klobuchar model, which has low correction accuracy. Another solution is to use differential ionospheric delay correction data provided by third parties, which have a high correction accuracy. However, this solution requires the support of obtaining precise ionospheric data through a customized data network. In this study, the authors model the ionospheric delay from the observations themselves and build a mathematical model of the ionospheric delay and its variations. A new and improved Hatch filtering algorithm is proposed. The carrier smoothing of code pseudorange observations in this algorithm has high accuracy and stability because it eliminates the influence of ionospheric delay and variation simultaneously. A sliding average window with an overlap segment is designed for this algorithm, which eliminates the discontinuity in the smoothed data. The analyses in this paper have been performed by means of a set of experimental tests using Trimble, some NETR9 receivers, and a GNSS data processing platform called GNSSer. The executed field tests cover two low kinematic experiments: a ship experiment and a trolley experiment. The results show that the improved Hatch filtering algorithm eliminates most of the influence of the ionosphere and obtains carrier-smoothing pseudorange observations with smoothness and high accuracy. In the kinematic experiment of the ship position, single-frequency pseudorange observations are used to establish an observation equation. The horizontal accuracy of the ship position is approximately 50 cm, and the vertical accuracy is less than 1 m, which is an order of magnitude higher than the accuracy of the classic Hatch filter. However, the kinematic positioning of the trolley is better, the accuracy is within 10 cm in the horizontal and 10–15 cm in the vertical. Both tests show consistent results that the improved Hatch filter has excellent single-frequency, single-epoch positioning performance under low kinematic conditions. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Global Navigation Satellite System (GNSS) pseudorange observations are the pseudo-distance between a navigation satellite receiver and a satellite, and the accuracy of these observations is low. The accuracy of the carrier phase observation can reach the order of millimeters, as it needs to repair cycle slips and fix ambiguities [1,2]. The carrier smoothing of code pseudorange observations combines the advantages of both and can greatly improve ⇑ Corresponding author at: Army Engineering University of PLA, Nanjing, Jiangsu 210007, China. E-mail address: [email protected] (Q. Zhang). https://doi.org/10.1016/j.measurement.2019.07.030 0263-2241/Ó 2019 Elsevier Ltd. All rights reserved.

the accuracy of pseudorange observations without introducing ambiguity. The classic carrier smoothing of the code pseudorange algorithm was proposed by Ron Hatch [3], which is also known as the Hatch filter algorithm. This algorithm relies on the fact that the ionospheric delay and phase delay are substantially equal in a one order term, and the variation in ionospheric delay between adjacent epochs is small. The ambiguity and ionospheric delay are obtained by the superposition averaging of multiple epochs to further improve the accuracy of pseudorange observations. The classic Hatch filter can improve the accuracy of pseudorange observations significantly, but it is more severely affected by the ionosphere and leads to divergence in the smoothed pseudorange. In view of the shortcomings of the classic carrier smoothing of

869

Q. Zhang et al. / Measurement 146 (2019) 868–878

the code pseudorange algorithm, some scholars have improved this algorithm. Several improved strategies for the original Hatch filter have been proposed. Hatch proposed using individual epoch weights [4]; Lachapelle proposed a weight reduction algorithm based on epochs [5]; and Meyerhof [6], Hwang [7] and HofmannWellenhof [8] proposed similar algorithms. Since the ionospheric delay is usually unknown to single-frequency GNSS users, accurate ionospheric delay is typically used for high-integrity applications from the Local Area Augmentation System (LAAS) or the Wide Area Augmentation System (WAAS). Park [9] uses the ionospheric physical information provided by the differential base station to establish the ionospheric model and propose an optimal smooth pseudorange considering the ionospheric changes in the ionosphere in different regions, days, seasons, and satellite elevation angles. The differentiation GPS (DGPS) single-frequency receiver algorithm has higher positioning accuracy and robustness than the classic Hatch filter. Kim [10], Lee [11,12,13] and Byungwoon [14] use linear regression to model the rate of change in ionospheric delay, comprehensively considering the characteristics of the satellite height cutoff angle and multipath noise, and propose an adaptive filtering algorithm that can improve GPS static positioning. Walter and Luo et al. from Stanford University [15,16] noted that single-frequency WAAS and LAAS receivers mounted on kinematic aircraft are often affected by changes in the ionospheric layer, especially at the equator, the Antarctic and the Arctic. An algorithm was proposed for this region to mitigate the effects of the total electron content (TEC). Helena designed an algorithm to reduce the accuracy attenuation caused by an ionospheric gradient change in the Hatch filter in kinematic mode by using a Kalman filter [17]. To solve the error problem of carrier phase smoothing in the foundation enhancement system, Huang [18] calculated the width of the smoothing window by using the satellite elevation angle, ionospheric variation and the distance change in the user to the reference station and reduced the error of the classic Hatch filter with time constants of 100 s by 68%. To cope with the divergence problem of the classic Hatch filter, Zhang [19] proposed using the ionospheric delay variation and the satellite elevation angle to adaptively change the width of the smoothing window. Almost all Hatch filters and their variants are linear, but when considering the final pseudo-distance estimation of measurement noise statistics, some scholars have proposed using nonlinear methods for the calculation, but this method has limited improvement in precision and increased computational complexity [20]. Cheng [21], Halimi [22] and Zhou [23] used Doppler measurements to smooth the phase observations to reduce the effect of ionospheric delay change on kinematic positioning. The current research results include the use of part of a linear fitting ionospheric delay for subtraction processing [9–14], the use of external differential information to correct the ionospheric delay [15,16], and external physical and geometric parameters to estimate the width of the ionospheric correction window [9,18,19]. These algorithms require external data support, where the correction algorithm is cumbersome, and the accuracy improvement is not obvious. In this paper, the authors try to solve the problem that the ionosphere is not fully modelled in the classic single-frequency Hatch filter algorithm. In the second chapter, an error analysis of the carrier-smoothing data model for the code pseudorange is carried out. The third chapter analyses the influence of the ionospheric delay on smooth pseudorange observations. The fourth chapter introduces the improved Hatch filtering model and the sliding window model with overlapping segments and describes the calculation process of the algorithm in detail. Chapter 5 contains two GNSS kinematic positioning tests (sea and ground) to verify the accuracy of the proposed algorithm. The sixth chapter concludes the advantages of the algorithm and discusses the accuracy of the algorithm.

2. Hatch filter and its accuracy If the influences of tropospheric delay, clock biases in the satellite/receiver and the multipath are ignored, the GNSS carrier phase observation and pseudorange observation can be expressed as:




Lk ¼ qk  Ik  N þ eL

ð1Þ

Pk ¼ qk þ Ik þ eP

where Lk and Pk represent the carrier phase distance observation value and the pseudo-distance observation, respectively. qk is the geometric distance between the receiver and the satellite. N is the ambiguity in the carrier phase observation (represented by distance), which absorbs some hardware delay. I is the ionospheric delay (represented by distance). eL and ep are observation noise. 2.1. Classical mathematical model of the Hatch filter By subtracting the two observation equations in formula (1) and absorbing the observed noise into the observation, we can obtain the following formula:

Pk  Lk ¼ 2Ik þ N

ð2Þ

Ambiguity and ionospheric delay are considered constant over a period of time (k epochs). The sum of two times the ionospheric delay and ambiguity is represented by A. 

A ¼ 2 I þN

ð3Þ

Ak is the average of A in k epochs, and i is the ordinal number of epochs.

Ak ¼

k 1X ðPi  Li Þ k i¼1

ð4Þ

The smooth pseudorange on epoch k can be obtained by formulas (2) and (4) 

Pk ¼ Lk þ Ak ¼ Lk þ

k 1X ðP i  L i Þ k i¼1

ð5Þ

The recursive formula of carrier phase smoothing pseudodistance is obtained by subtracting the adjacent two epoch elements. 



Pk ¼ P k1 þ ðLk  Lk1 Þ þ ðAk  Ak1 Þ

ð6Þ

Where

Ak  Ak1 ¼

1 ðPk  Lk  Ak1 Þ k

ð7Þ

Placing (7) into the formula (6) and we can obtain the following formula 

Pk ¼

1 Pk þ k

   1  Pk1 þ Lk  Lk1 1 k

ð8Þ

Formula (8) is a weighted average formula, the essence of which is the recursion of the adjacent two epochs. At the same time, we find that the first two items on the right-hand side of the formula (6) are the same as the last item on the right-hand side of formula (8), which includes only the smoothing value Ak-1 of the former k-1 ^



epochs. P k is the estimated smooth pseudorange Pk , as follows. ^



P k ¼ Pk1 þ ðLk  Lk1 Þ 





ð9Þ ^

where I k1 ¼ I k and P1 ¼ P k ¼ P1 . Formulas (8) can be expressed in the following form.

870

Q. Zhang et al. / Measurement 146 (2019) 868–878



^

Pk ¼ wk Pk þ ð1  wk ÞP k

ð10Þ

where wk is the weight corresponding to the measured pseudorange Pk of the kth epoch, which is the reciprocal of the number of epochs wk = 1/k. The number k of fixed epochs needs to be m; that is, when k > m, wk = 1/m. Further, if wk = 0, the formulas (9) and (10) are equal. 2.2. Theoretical accuracy analysis of HATCH filtering Let the variance in the carrier phase observation be r2L and the variance in the pseudorange observation be r2P . According to the law of error propagation and ignoring the relationship between the pseudorange and phase observations, the smooth pseudorange variance can be obtained by formulas (4) and (5).

r2P k ¼ r2L þ Since

1 2 r þ k P

r2L



ð11Þ

r2P  r2L , then

r rP k  pPffiffiffi k

ð12Þ

Formula (12) shows that the smoothed pseudorange accuracy is improved by k times, that is, the higher the number of epochs participating in smoothing, the higher the precision is. Let the accuracy of the C/A code of the GPS L1 frequency be 3 m and the accuracy of the P code be 0.3 m [24]. The convergence of the theoretical accuracy is shown in Fig. 1. Fig. 1 shows that if the number of epochs increases, the accuracy of the smooth pseudorange observation is continuously improved, but the speed of the improvement is continuously decreased. The windows of different lengths are specified for both the C/A code and the P code, and the smoothing precision is calculated. The results are shown in Table 1. Table 1 shows that the accuracy of the carrier smoothing of the code pseudorange can reach the decimeter and sub-decimeter levels quickly, regardless of the C/A code or the P code. This is enough for most kinematic navigation applications. In addition,

since the effective period of the satellite observed by the mid-latitude station is usually only 3–6 h, data at the 1 s sampling rate are up to approximately 20,000 epochs. Table 1 shows that the theoretical accuracy limit of the carrier-smoothing pseudorange observation (C/A code), with a sampling rate of 1 s, is 0.02 m. In addition, the accuracy of the carrier smoothing of the code pseudorange observation value is proportional to the reciprocal of the square root of the number of epochs, that is, the higher the number of epochs is, the higher the smoothing precision. This shows that the observation data with a high sampling rate are more conducive to the improvement in precision. If the sampling rate is 0.1 s with a C/A code, a smoothing accuracy of 0.1 decimeters can be achieved in only 1.5 min. 3. Influence of the ionosphere on Hatch filtering 3.1. Hatch filtering error caused by ionospheric delay variation From formula (2), the current epoch is k, the starting epoch is 1, the ionosphere of the starting epoch is I1, and DIi;1 is the difference in the ionospheric delay between epoch i and epoch 1. Taking into account ionospheric changes, the following formulas exist:

Ak ¼

k k 1X 2X ðP i  Li Þ ¼ A þ DIi;1 k i¼1 k i¼1

ð13Þ

Let the ionospheric delay residual of Ak be IAk ; then,

A ¼ Ak  I A k

ð14Þ

According to formula (13),

I Ak ¼

k 2X DIi;1 k i¼1

ð15Þ

Formula (14) shows that there is a fixed difference between the estimated value Ak of the carrier smoothing of code pseudorange and the true value A, that is, the ionospheric delay variation. This difference is related to the sum of the ionospheric delay variations in each epoch participating in smoothing. The size of this

Fig. 1. Theoretical precision of carrier-smoothing of code pseudorange.

871

Q. Zhang et al. / Measurement 146 (2019) 868–878

3.2. Analysis of the evolution in the pattern of the ionosphere

Table 1 Smoothing accuracy of the pseudorange in windows of different lengths. Number of epochs

Length of the window

C/A code (3 m)

P code (0.3 m)

1 30 60 120 300 600 1200 3600 7200 14,400 21,600

1s 30 s 1 min 2 min 5 min 10 min 20 min 1h 2h 4h 6h

3.000 0.548 0.387 0.274 0.173 0.122 0.087 0.050 0.035 0.025 0.020

0.300 0.055 0.039 0.027 0.017 0.012 0.009 0.005 0.004 0.003 0.002

L1  L2 ¼

Let



ð16Þ

Let

IP



k

¼ 2DIk;1  IAk

k 1X ¼ 2 DIk;1  DIi;1 k i¼1



   bL1r  bL2r  bL1s  bL2s c þ ðk1 N1  k2 N2 Þ   2 2  A 1=f 1  1=f 2 þ eL1 L2

DL ¼ L1  L2 ,

! ð17Þ

From formula (16), there is the following formula:

Dbr ¼ bL1r  bL2r ,

DL ¼ ðDbr  Dbs Þc þ DN  DIA þ eL1 L2

ð18Þ

I is the correction of the original smooth pseudorange with respect to the smoothed pseudo-range of the ionospheric delay variation. The formula derivation in this section shows that for the smooth ^ k after correcting the ionosphere of the specified pseudorange P 

epoch, if IP is positive, the P k value will be too small, and vice versa.

ð20Þ

Also let DB ¼ ðDbr  Dbs Þc þ DN, then

DL ¼ DB  DIA þ eL1 L2

ð21Þ

In addition, DIA ¼ DB  DL þ eL1 L2 . That is, the inter-frequency carrier phase difference can obtain the hardware delay difference and the ionospheric delay difference values. In formula (21), DB can be regarded as a constant in the short term (i.e., within the visible range of the satellite). By multiplying both sides of (21) by the same coefficient, the sum of the distances for each ionospheric frequency and hardware delay can be obtained, which is represented by the variable IB.

8  f2 f2  > < IB1 ¼  f 2 f2 2 DB þ IL1 ¼  f 2 f2 2 DL  eL1 L2 1 2 1 2 2 2   > : IB2 ¼  2f 1 2 DB þ IL ¼  2f 1 2 DL  eL L 2 1 2 f f f f 1



^ k ¼ P k þ I P P k

  2 2 DIA ¼ A 1=f 1  1=f 2 ,

ð19Þ

Dbs ¼ bL1s  bL2s , DN ¼ k1 N 1  k2 N 2 Then (19) can be simplified to

difference is twice the average of the ionospheric delay variation in each epoch relative to the first epoch. Since A is a constant, if IAk is positive, Ak becomes smaller, which further results in a smaller smoothing pseudo-margin, and vice versa. From formulas (2) and (5), a smoothed pseudorange with an ionospheric delay variation correction can be obtained.

^ k ¼ Pk þ 2DIk;1  IA P k

To study the evolution in the pattern of the ionosphere, this paper uses the results of a dual-frequency carrier data calculation for analysis. Two carrier observations of different frequencies are subtracted to obtain formula (19)

2

1

ð22Þ

2

The pattern of error propagation shows that the sum of the obtained ionospheric and hardware delay IB is on the same order of magnitude as the carrier-phase observation L, which can reach the millimeter level. For the GPS L1 signal, there is

  IL1 þ 1:5457DB ¼ 1:5457 DL  eL1 L2

ð23Þ

k

Fig. 2. Ionospheric delay variation (4–8 h) (unit: m).

Table 2 Statistical results of the ionospheric delay variation in the day (absolute value). period (hours)

Maximum value (m)

Minimum value (m)

Average value (m)

RMS

Data amount

0–4 4–8 8–12 12–16 16–20 20–24

0.056 0.045 0.067 0.047 0.0426 0.045

0.0004 0.0004 0.0001 0.0003 0.0002 0.0004

0.0016 0.0013 0.0022 0.0009 0.0011 0.0013

0.003645 0.003002 0.003670 0.003613 0.003064 0.003241

111,692 123,694 126,463 133,118 134,061 130,727

872

Q. Zhang et al. / Measurement 146 (2019) 868–878

By determining the difference between the epochs, the rate of change in the ionosphere can be accurately determined as follows:

8  f2  > < DIL1 ;k ¼ IL1 ;k  IL1 ;k1 ¼  f 2 f2 2 DLk  DLk1  eDLk DLk1 1 2 2   > : DIL ;k ¼ IL ;k  IL ;k1 ¼  2f 1 2 DLk  DLk1  eDL DL 2 2 2 k k1 f f 1

ð24Þ

2

For the GPS L1 signal, there is

  DIk ¼ Ik  Ik1 ¼ 2:5457 DLk  DLk1  eDLk DLk1

Fig. 4. Segmented sliding window.

ð25Þ

Differences in the carrier-phase observations of the two frequencies can eliminate most of the error, leaving only the ionospheric delay, hardware delay, and ambiguity equivalent distance. 3.3. Analysis of the ionospheric delay variation with IGS data Since kinematic positioning requires GNSS observation data with a high sampling rate, this paper uses the observation data of a 1 s sampling rate provided by the International GNSS Service (IGS) [25] to analyze the ionospheric delay variation. The L1 band data of the KIRI station (172.92E-1.35N) are analyzed below using a period of 4 h for a total of 6 time periods. The height cutoff angle of the satellite is 15 degrees. Fig. 2 shows the ionospheric delay variation (can be calculated by Eq. (25)) measured by different satellites during the second time period. Table 2 shows the statistical results for all time periods in the day. By analyzing Fig. 2 and Table 2, we can find that (1) the ionospheric delay variation in each satellite is on the same magnitude; (2) the average variation in all satellites is between 2 and 5 cm, and the maximum variation in the change does not exceed 10 cm; and (3) for high-precision kinematic positioning, when using Hatch filtering, the ionospheric delay variation in each observation must be considered and cannot be estimated as a constant. 4. Modelling ionospheric delay variation and improved Hatch filter 4.1. Self-modelling of the ionospheric delay variation The ionospheric delay variation is the main cause of divergence in the classic Hatch filter. Therefore, to reduce the impact of ionospheric delay variation, the impact of the ionosphere must be analyzed and deducted. The ionospheric delay variation can be estimated using the Klobuchar model in the navigation message but with lower precision. Without the use of an external precision ionospheric model, this paper uses a single-frequency carrier and pseudorange combination to model the ionospheric delay variation and ultimately eliminate its effects. Since the ionospheric delay variation established in this paper does not depend on any external data, only the measured carrier and pseudorange observations are used; therefore, it is called self-modelling. The following formula can be obtained from formula (2):

Pk  Lk ¼ 2Ik þ N þ ePk

ð26Þ

Fig. 5. Segmented sliding window with overlap.

Fig. 6. Overlapping relationship between windows.

In the above formula, the pseudorange subtracts the sum of the ionosphere and the ambiguity equal to 2 times the carrier distance and includes the hardware delay error, which is constant in the short term. The change in this value is directly reflected in the change in the ionosphere. Therefore, directly subtracting the difference between the pseudorange and a carrier distance, we can find the pattern of the ionospheric delay variation, as seen in formula (27):

IN ¼ I k þ

N Pk  Lk ePk ¼  2 2 2

ð27Þ

In formula (27), the ionospheric delay variation can be calculated using the observations of multiple epochs. The general form of a unitary n-degree polynomial with epoch t as an independent variable is as follows:

f ðt Þ ¼

n X

ai tii

ð28Þ

i¼0

where f(t) is the ionospheric fitting value at time t, ai is the i-th fitting term coefficient of the polynomial, and a0 is a constant equivalent to the initial ionospheric delay I1. By dividing both sides of formula (26) by 2, the following formula is obtained.

yk ¼

P k  Lk N eP ¼ I1 þ þ 2 2 2

ð29Þ

when the observed epoch is tk, an error equation can be established according to formula (28).

vk ¼ Fig. 3. Epoch sliding window.

n X i¼0

ai t ii  yk

ð30Þ

Q. Zhang et al. / Measurement 146 (2019) 868–878

873

Fig. 7. Flow chart of the improved Hatch filtering algorithm.

where vk is the observed residual of the epoch tk. If there are multiple epoch observations, the matrix of the error equation is expressed as

2 3 ^0 1 t 10 a 6 7 6a 6 1 t 11 6 ^1 7 7; A ¼ 6 ^¼6 a 6 . .. 6 . 7 6 ... 4 . 5 . 4 ^i a 1 t 1n 2

3

2 3 y0 7 2 n7 7 6 t1 t1 6 y1 7 7 6 7 ; L ¼ 7 . .. . 6 . 7 4 . 5 . .. 7 5 yn t2n t nn t20

t n0

ð31Þ

The coefficient vector is estimated by the least squares method as follows:

 1 b ¼ AT A AT L a

ð32Þ

Using formulas (28) and (32), the fitting values of the ionospheric delay for any epoch can be obtained. On this basis, the differential calculation can be performed to obtain a carrier phase smoothing measurement with higher accuracy. Considering the ionospheric delay and its rate of change, this paper uses a second-order polynomial fit. The specific window of the observed epoch will be analyzed in Section 4.2.

4.2. Optimization of the ionospheric modelling window To estimate the two parameters of the ionospheric delay and its rate of change, the quadratic polynomial regression algorithm is used in this paper. For the observation selection of the estimation parameter n, the data in a window are primarily used for solving, which is also called the sliding window fitting method. However, if the amount of observation data is large and the window is large, the calculation amount of the conventional sliding window fitting window becomes very large and time consuming. Therefore, the methods of ‘‘epoch sliding window” and ‘‘segmented sliding window” are utilized in practice. The sliding window for each epoch uses a data block that processes the number of epochs n, expands into a data block of m epoch sizes centered on an epoch and is referred to as window data of size m. As shown in Fig. 3, one of the small cells represents the data of one epoch, and the green filled cell is the current epoch. In the data processing, the epoch moves forward or backward until all data are processed. This method is called the ‘‘epoch sliding window”. However, as the window increases, the amount of computational effort for this method will become increasingly larger and even difficult to calculate.

874

Q. Zhang et al. / Measurement 146 (2019) 868–878

10

De

Dn

Du

Posioning error(m)

5 0 -5 -10 -15 -20 5:40

6:10

6:40

7:10

7:40

8:10

Epoch(GMT, hour) Fig. 11. Positioning error of the classic Hatch filter.

3

De

Dn

Du

Fig. 8. GNSS reference station and route of the GNSS kinematic station.

17

Ionospheric delay(m)

1 0 -1 -2

19 G02 G07 G13 G20

15 13 11

G05 G09 G17 G25

-3

G06 G12 G19 G29

5:40

6:10

6:40

7:10

7:40

8:10

Epoch(GMT, hour) Fig. 12. Positioning error of the improved Hatch filter (5 min window).

9 7 5

3 1 5:40

6:10

6:40

7:10

7:40

8:10

Epoch(GMT,hour) Fig. 9. Ionospheric delay (from the CODE spherical harmonic model).

3 Du

Dn

2

Posioning error(m)

Posioning error(m)

2

De

To solve the above problem, a segmented sliding window is used. Divide all data into N windows of size m epochs, as shown in Fig. 4. Among them, the windows do not overlap and are connected end to end and stitched together. This method can process large amounts of data quickly. However, the independent window causes an artificially independent cut, resulting in a fracture of the fitted data. To prevent data from framing, the authors propose a segmented sliding window with overlap. This method has the characteristics of both the epoch sliding window and the segmented sliding window, which reduces the amount of data calculation and ensures the continuity between the windows, as shown in Fig. 5. The overlapping relationship of each window is shown in Fig. 6, and the green part is an overlap (Fig. 7).

1

4.3. Calculation process of the improved Hatch filter 0 -1 -2 -3

5:40

6:10

6:40

7:10

7:40

8:10

Epoch(GMT, hour) Fig. 10. Positioning error of the original pseudorange.

The algorithm flow chart of the improved Hatch filter proposed by the authors is described below. In the first step, the singlefrequency carrier and original pseudorange observation data are read and composed of carrier-smoothing pseudo-side measurements. In the second step, two cache windows are created: one is the ionosphere fitting window (this study uses a fixed 20-min window fit), and the other is a smooth pseudorange window. In the third step, the self-modelled ionospheric delay variation is used to improve the accuracy of the carrier observation. In the

875

Q. Zhang et al. / Measurement 146 (2019) 868–878 Table 3 Positioning errors in ENU (±m). Method

E

N

U (vertical)

Horizontal

Original pseudorange Hatch filter The improved Hatch filter

0.523 3.578 2.868 2.073 1.608 1.098 0.563 0.314 0.218 0.227

1.567 2.350 2.824 2.143 1.702 1.180 0.631 0.383 0.286 0.294

2.942 6.691 2.804 2.151 1.746 1.359 0.940 0.813 0.840 0.877

1.652 4.281 4.025 2.982 2.341 1.612 0.845 0.495 0.359 0.371

1h 40 min 30 min 20 min 10 min 5 min 2 min 1 min

5. Kinematic positioning analysis based on improved Hatch filtering

Fig. 13. The RMS of the positioning errors in the ENU.

fourth step, an improved calculation of the smooth pseudorange observation is established. In the fifth step, the improved smooth pseudo-measurement is used to perform the positioning solution. The above steps can be cycled for calculation.

To test the improved Hatch filtering algorithm proposed in this paper, the authors conducted two kinematic positioning experiments. In the first experiment, a GNSS receiver was fixed on a ship, while the second experiment used a trolley to conduct kinematic experiments on a flat playground. In both experiments, the difference technique was used to obtain the higher precision position value, which was used as the true value to evaluate the accuracy of the improved Hatch filtering algorithm proposed in this paper. The detailed results and comparative analysis of the two experiments are shown below. In both experiments, the Trimble Net-R9 R9 GPS receiver was used [26]. The improved Hatch filter algorithm is implemented on the GNSSer (GNSS Data Parallel Processer, V1.4) software platform [27,28]. GNSSer is a research-oriented GNSS data processing software developed by our team that aims to provide GNSS computing services with high precision, parallelism and cloud mode.

Fig. 14. Measurement trolley used in the kinematic experiment.

876

Q. Zhang et al. / Measurement 146 (2019) 868–878

1

Posioning error(m)

De

Dn

Du

0.5

0

-0.5

-1 Fig. 15. Motion trajectory.

4:25

4:35

4:45

4:55

Epoch(GMT, hour) Fig. 18. Positioning error of the classic Hatch filter.

G05 G15

Ionospheric delay (m)

7

G10 G18

G12 G20

G13 G21

5

3

1 4:25

4:35

4:45

4:55

Epoch(GMT, hour) Fig. 16. Ionospheric delay (from the CODE spherical harmonic model).

Fig. 17. Positioning error of the original pseudorange.

the kinematic station and the base station, and the positioning result is used as the true value of the running track of the kinematic station to evaluate the accuracy of the improved Hatch filtering algorithm proposed in this paper. Using the original pseudorange, the classic Hatch filter and the improved Hatch filter proposed in this paper, the pseudorange of the L1 frequency of the kinematic station is pseudo-ranged by epoch. The difference between the calculated position results of the three methods and the double-difference result (true value) is used to calculate the coordinate difference in the north, east and vertical directions. GPT2/VMF1 model is used in tropospheric correction [29]. The ionospheric correction uses the results of the non-difference non-combination calculations, and the troposphere uses model corrections. Fig. 9 shows the ionospheric delay of the satellite provided by the CODE spherical harmonic model [28]. The ionospheric delay of each satellite is at the meter level, and there are irregular changes. Figs. 10–12 show the positioning results of the three algorithm methods. In Fig. 10, the positioning results of the original pseudorange are at the meter level. As shown in Fig. 11, due to the influence of ionospheric delay, the classic Hatch filter has a divergence phenomenon with the accumulation of ionospheric delay. As shown in Fig. 12, the improved Hatch filter proposed in this paper can suppress divergence. At the same time, the positioning accuracy is improved to the sub-meter level. To evaluate the accuracy of the improved Hatch filter for different size windows, 8 windows were selected in this experiment for the positioning solution. The statistical results are shown in Table 3. For comparison, the original pseudorange and the classic Hatch filtering positioning results are also counted. The comparison results are shown in Fig. 13. In the kinematic positioning, the window precision of 2 min to 5 min is the best, the highest level of precision is 35.9 cm, and the height is 81.3 cm.

5.1. Trajectory determination of a ship

5.2. Kinematic experiment of a trolley

The location of the experiment was offshore of Jiangsu, China, and was carried out at 5:38 to 8:44 UTC on September 29, 2017. Maritime kinematic experiments were conducted for a total of 3 h, with a total journey of approximately 32 km. Among them, a GNSS receiver was installed as a kinematic station on the ship, and another GNSS receiver was installed as a reference station near the departure terminal (as shown in Fig. 8). The carrier doubledifference observation is established by the observation data of

In the land kinematic positioning test, a small flatbed (Fig. 14) was used to make a circular motion around a football field, with a total of 4.5 laps and a distance of approximately 1.8 km (Fig. 15). The experimental time is 4:25–5:00 UTC for a total of 35 min, and the data sampling rate is 1 s. In addition to the GNSS receiver on the trolley being used as a kinematic station, a static reference station is also set up in the southeast direction 800 m away from the centre of the football field to obtain the high-

877

Q. Zhang et al. / Measurement 146 (2019) 868–878

precision trajectory of the trolley as a true value through differential observation (Fig. 16). Similar to the shipboard data, the data of the land vehicle in Fig. 9 show that the ionospheric delay of each satellite is at the meter level, and there are irregular changes. Using the pseudorange of L1, the position solution is performed by the epoch parameter adjustment method. The ionospheric correction uses the results of the single-frequency ionospheric delay calculation of the static station and uses three smooth pseudorange algorithms to perform the positioning solution. The double-difference solution is used as the true value of the kinematic positioning data to evaluate the accuracy of the three methods. Figs. 17–19 show the positioning results of the three methods. 1

Posioning error (m)

De

Dn

Du

0.5

0

6. Conclusion and discussion -0.5

-1 4:25

4:35

4:45

4:55

Epoch(GMT, hour) Fig. 19. Positioning error of the improved Hatch filter (5 min window).

0.8 E

N

U

0.6

RMS (m)

As shown in Fig. 17, the positioning result using the original pseudorange is much smaller in the horizontal direction than it is in the vertical direction. As shown in Fig. 18, due to the influence of ionospheric delay, the conventional Hatch filter has obvious divergence with the accumulation of ionospheric delay. As shown in Fig. 19, the improved Hatch filter proposed in this paper can be used very quickly. The convergence of the positioning results can be obtained, and the positioning accuracy can be improved to the centimeter level (Fig. 20). To evaluate the accuracy of the improved Hatch filter when using different sizes of windows, five windows are selected for the positioning calculation. The statistical results are shown in Table 4. For comparison, the original pseudorange and the classic Hatch filtering positioning results are also counted. The positioning accuracy with the improved Hatch filtering algorithm is an order of magnitude higher than that of the original pseudorange and classic Hatch filter. The positioning result can be 6–7 cm in the horizontal direction and 10 cm in the vertical direction, showing excellent positioning performance. The five kinds of windows can achieve extremely high positioning accuracy and are not sensitive to the window size.

0.4

0.2

0.0

The improved HATCH filter

Fig. 20. The RMS of the positioning errors in the ENU.

Pseudorange positioning can only achieve positioning results at the metre level. The classic Hatch filter uses the carrier observations to smooth the pseudorange but is susceptible to divergence caused by ionospheric variations. In the two kinematic experiments in this paper, the classic Hatch filter has divergence and reduces the accuracy of positioning. The improved Hatch filtering algorithm proposed in this paper can eliminate the influence of most ionospheric changes by using carrier observation data for regression modelling without external ionospheric correction data or differential methods. At the same time, for the method proposed in this paper, a sliding window with overlapping segments is designed to solve the problem of discontinuity in the carrier phase smoothing data. In the kinematic experiment of shipborne data, the singlefrequency GPS data and the proposed algorithm can achieve better than 50 cm results in the horizontal direction, and the data in the vertical direction can be better than 1 m, which is an order of magnitude better than the original pseudorange and classic Hatch filter results. While the results of the land kinematic experiment are higher, the horizontal displacement can be better than 10 cm, and the vertical direction can obtain a precision of 10 cm–15 cm. For the single-frequency non-differential pseudorange solution, the algorithm in this paper achieves high accuracy. The algorithm in this paper is expected to provide an efficient and highly accurate solution for single-frequency low kinematic GNSS applications. The classic Hatch filter is actually designed for static or low kinematic receivers, and the proposed improved Hatch filter has excellent performance. However, its performance in environments that are highly kinematic or suffering from cycle jumps needs further validation.

Table 4 Positioning errors in the ENU (±m). Method

E

N

U (vertical)

Horizontal

Original pseudorange Hatch filter The improved Hatch filter

0.195 0.134 0.0445 0.0445 0.0450 0.0488 0.0523

0.278 0.289 0.0445 0.0446 0.0444 0.0497 0.0549

0.730 0.468 0.1036 0.1045 0.1034 0.1285 0.1444

0.340 0.319 0.0629 0.0630 0.0633 0.0697 0.0758

20 min 10 min 5 min 2 min 1 min

878

Q. Zhang et al. / Measurement 146 (2019) 868–878

Declaration of Competing Interest None declared.

[12] [13]

Acknowledgment This work was carried under funding from The National Natural Science Foundation of China (Grant nos. 41604024, 41804006).

[14] [15]

[16]

References [1] E.D. Kaplan, Understanding GPS: principles and application, Artech House Mobile Commun 59 (5) (2005) 598–599. [2] B.W. Parkinson, P. Enge, P. Axelrad, et al., Global Positioning System: Theory and Applications, vol. II, American Institute of Aeronautics and Astronautics, 1996. [3] R. Hatch, The synergism of GPS code and carrier measurements, in: International Geodetic Symposium on Satellite Doppler Positioning. International Geodetic Symposium on Satellite Doppler Positioning, 1982, pp. 1213–1231. [4] R. Hatch. Dynamic Differential GPS at the centimeter level, in” Proceedings of 4th Int. Geodetic Symposium on Satellite Positioning, vol. vol 2, Austing, TX, April May 1986, pp. 1287–1298. [5] G. Lachapelle, J. Hagglund, W. Falkenberg, P. Bellemare, M. Casey, M. Eaton, GPS land kinematic positioning experiments, in: Proceedings of 4th Int. Geodetic Symposium Satellite Positioning, vol. 2, Austin, TX, April-May 1986, pp. 1327–1344. [6] S.L. Meyerhoff, A.G. Evans, Demonstration of the combined use of GPS pseudorange and Doppler measurements for improved dynamic positioning, in: Proceedings of 4th Int. Geodetic Symposium Satellite Positioning, vol. 2, Austin, TX, April-May 1986, pp. 1397–1409. [7] P.Y.C. Hwang, R.G. Brown, GPS navigation: combining pseudorange with continuous carrier phase using a Kalman filter, J Inst Navig 37 (2) (1990) 181– 196. [8] B. Hofmann-Wellenhof, H. Lichtenberger, J. Collins, Global Positioning System Theory and practice, Springer, Viena, 2001. [9] B. Park, K. Sohn, C. Kee, Optimal hatch filter with an adaptive smoothing window width, J. Navig. 61 (3) (2008) 435–454. [10] E. Kim, T. Walter, J.D. Powell, Adaptive carrier smoothing using code and carrier divergence, Proc. Natl. Tech. Meet. Inst. Navig. (2007). [11] H.K. Lee, J.G. Lee, G.I. Jee, An efficient GPS receiver algorithm for channelwise multipath detection and real-time positioning, in: Proceedings of the Institute

[17] [18]

[19]

[20] [21] [22]

[23] [24] [25]

[26] [27] [28]

[29]

of Navigation 2002 National Technical Meeting, San Diego, CA, pp. 265–276, 2002. H.K. Lee, C. Rizos, Position-domain Hatch filter for kinematic differential GPS/ GNSS, IEEE Trans. Aerospace Electr. Syst. 44 (1) (2004) 30–40. H.K. Lee, C. Rizos, Position-domain Hatch filter for kinematic differential GPS/ GNSS, IEEE Trans. Aerosp. Electron. Syst. 44 (1) (2008) 30–40. P. Byungwoon, S. Kyoungho, K. Changdon, Optimal Hatch filter with an adaptive smoothing window width, J. Navig. 61 (3) (2008) 20. T. Walter, S. Datta-Barua, J. Blanch, et al., The effects of large ionospheric gradients on single frequency airborne smoothing filters for WAAS and LAAS, Proc. Natl. Tech. Meet. Inst. Navig. (2004) 103–109. M. Luo, S. Pullen, J. Dennis, et al., LAAS Ionosphere spatial gradient threat model and impact of LGF and airborne monitoring, in: Proceedings of international technical meeting of the satellite division of the Institute of Navigation, 2003. L. Helena, S. Jari, T. Jarmo, Complementary Kalman filter for smoothing GPS position with GPS velocity, ION GPS/GNSS, 2003, Portland, OR. Zhenggang Huang, Huang, et al., A new optimal Hatch filter to minimize the effects of ionosphere gradients for GBAS, Chin. J. Aeronaut. 21 (6) (2008) 526– 532. X. Zhang, P. Huang, Optimal Hatch filter with an adaptive smoothing time based on SBAS, in: International Conference on Soft Computing in Information Communication Technology, 2014. K. Mazher, M. Tahir, K. Ali, GNSS pseudorange smoothing: linear vs non-linear filtering paradigm, in: Aerospace Conference, 2016. P. Cheng, Remarks on Doppler-aided smoothing of code ranges, J. Geod. 73 (1) (1999) 23–28. A. Halimi, C. Mailhes, J.Y. Tourneret, et al., Bayesian estimation of smooth altimetric parameters: application to conventional and delay/doppler altimetry, IEEE Trans. Geosci. Remote Sens. 54 (4) (2015) 2207–2219. Z. Zhou, B. Li, Optimal Doppler-aided smoothing strategy for GNSS navigation, GPS Sol. 21 (1) (2017) 197–210. P.K. Enge, The global positioning system: signals, measurements, and performance, Int. J. Wireless Inf. Networks 1 (2) (1994) 83–105. J.M. Dow, R.E. Neilan, C. Rizos, The international GNSS service in a changing landscape of global navigation satellite systems, J. Geod. 83 (7) (2009) 689-– 689. www.trimble.com/Infrastructure/Trimble-NetR9. www.gnsser.com. Linyang Li, Zhiping Lu, Zhengsheng Chen, et al., GNSSer: objected-oriented and design pattern-based software for GNSS data parallel processing, J. Spatial Sci. (2019), https://doi.org/10.1080/14498596.2019.1574245. P. Steigenberger, J. Boehm, V. Tesmer, Comparison of GMF/GPT with VMF1/ ECMWF and implications for atmospheric loading, J. Geod. 83 (10) (2009) 943– 951.