A Precise Heave Determination System Using Time-Differenced GNSS Carrier Phase Measurements

Abstract

In this study, a system that precisely determines the heave of ship hull was designed using time-differenced GNSS carrier phase measurement, and the performance was examined. First, a technique that calculates precise position relative to the original position based on TDCP measurement for point positioning using only one receiver was implemented. Second, to eliminate the long-cycle drift error occurring due to the measurement error that has not been completely removed by time-differencing, an easily implementable high-pass filter was designed, and the optimum coefficient was determined through an experiment. In a static experiment based on the precise heave measurement system implemented using low-cost commercial GNSS receiver and PC, the heave could be measured with a precision of 2 cm standard deviation. In addition, in a dynamic experiment where it moved up and down with an amplitude of 48 cm and a cycle of 20 seconds, precise heave without drift error could be determined. The system proposed in this study can be easily used for many applications, such as the altitude correction of fish detection radar.

1. INTRODUCTION

The heave of a ship refers to motion in the altitude direction, and the heave measurement is used to correct the ship hull motion due to waves in the case of offshore drilling, fish detection using radar, and marine exploration (Kuchler et al. 2011). Methods for measuring heave include a method that measures altitude change using a barometer, a method using an acceleration sensor, and a method using the Global Navigation Satellite System (GNSS) (Lee et al. 2015). In the method that measures altitude using a barometer, the altitude is measured based on the fact that atmospheric pressure decreases as the altitude from the mean sea level altitude increases. However, the atmospheric pressure at the reference sea level changes depending on the day, and it is significantly affected by the temperature and weather (Jan et al. 2008). The method using an acceleration sensor is a traditional method, and the heave is measured by obtaining displacement through the double integral of acceleration. However, to avoid the error accumulation occurring in the integration process, a high-cost precise sensor needs to be used (Kielland & Hagglund 1995, BEng 2007). In recent years, the method that measures heave using GNSS is widely used. Depending on the measurement type and utilization method, the technique using GNSS can be classified into the absolute positioning using code measurement, the relative positioning using the Real Time Kinematic (RTK) method (Kielland & Hagglund 1995, Godhavn 2000), and the measurement of altitude change based on the Time-Differenced Carrier Phase (TDCP) method (BEng 2007, Traugott 2011, Lee et al. 2015, Cho et al 2017). The absolute positioning using code measurement has a precision of several m, and thus cannot be used for precise heave determination. The RTK method can obtain precise relative altitude using carrier phase measurement. However, an additional cost is incurred since a reference station needs to be installed; and in the case of a ship that performs independent work on the sea, the application is difficult since the installation of a reference station is difficult.

In this study, a technique that determines heave by obtaining the altitude change based on the TDCP method was used. Using TDCP method, the heave can be determined using only one GNSS receiver. BEng (2007) reported that the heave value determined by the TDCP method was accurate for a short period of time, but the long-cycle drift error occurring due to the residual ionospheric, tropospheric, and multipath errors that had not been eliminated by the time differencing could be effectively eliminated using 276th-order high pass filter (HPF). However, considering the complexity of implementation and the delay of output, it is difficult to apply the 276th-order filter to a real-time system. In this study, a heave determination system based on the TDCP method using simple first-order HPF was designed, and the performance was examined. The heave determination system was implemented using ublox EVK-M8T receiver and Matlab, and the performance of the designed system was examined through an experiment.

The contents of this study are as follows. In Chapter 2, a technique that calculates heave using the TDCP model was derived. In Chapter 3, the heave measurement system was implemented using Matlab, and the problem of drift error occurrence was identified through an experiment. In Chapter 4, first-order HPF was additionally applied to resolve the problem mentioned in Chapter 3, and it was found that the drift error could be eliminated. Lastly, the conclusions were summarized in Chapter 5.

2. HEAVE DETERMINATION TECHNIQUE USING TIME-DIFFERENCED GNSS CARRIER PHASE MEASUREMENT

2.1 Carrier Phase Measurement Model

The pseudorange obtained using carrier phase can be expressed by Eq. (1).

\(\Phi(t)=r(t)+cB(t)+cb(t)+\lambda N-I(t)+T(t)+w(t)\)                                                                           (1)

where \(\Phi(t)\) denotes the carrier phase measurement measured by the receiver, \(r(t)\) the actual distance between the satellite and the receiver, \(B(t)\) the receiver clock error, \(b(t)\) the satellite clock error, \(\lambda\) the wavelength of the carrier frequency, \(N\) the integer ambiguity, \(I(t)\) the ionospheric delay error, and \(T(t)\) the tropospheric delay error. The ionospheric delay error is represented as a negative number since the size is identical to that of code measurement but the sign is different. \(w(t)~N(0,\sigma_\Phi^2)\) is the Additive White Gaussian Noise with a mean of 0 and a variance of \(\sigma_\Phi^2\).

The measurement noise of carrier phase measurement (mm level) is significantly smaller than that of code measurement (cm~m level), and thus precise positioning is possible using carrier phase measurement. However, for the use of carrier phase measurement, the integer ambiguity included in the carrier phase measurement needs to be calculated. When time differencing is used, the integer ambiguity is removed, and thus the problem of integer ambiguity can be avoided.

2.2 Receiver Heave Determination Using Time-Differenced Carrier Phase Measurement

Time differencing is a method that performs the differencing of measurements that are continuous in terms of time. The time-differenced GNSS carrier phase measurement for sufficiently close times, \(t_1\) and \(t_2\), can be expressed by Eq. (2).

\(\Phi\left(t_{12}\right)=r\left(t_{12}\right)+cB\left(t_{12}\right)+w\left(t_{12}\right)\)                                                                                                (2)

where \(•t12=•t2-•t1\). For a short time interval, the changes in the satellite clock error, ionospheric delay, and tropospheric delay can be ignored, and thus they are removed by time differencing. When the receiver maintains tracking, the integer ambiguity is removed by differencing as it is a constant. In the above measurement, when the satellite position is expressed as \(S(t)=[X(t),Y(t),Z(t)]\) and the receiver position as \(u(t)=[x(t),y(t),z(t)]\), the distance between the satellite and the receiver at time \(t\) can be expressed as \(r(t)=S(t)-u(t)\). In this regard, \(•\) is the Euclidean norm, and represents the distance.

To obtain the user position using the measurement of Eq. (2), \(r(t_2)\) is linearized at the receiver position of the previous time \(u(t_1)=[x(t_1),y(t_1),z(t_1)]\), it can be expressed as shown in Eq. (3).

\(r(t_2)=S(t2)-u(t1)+g(t2)⋅δu(t2)\)                                                                                              (3)

where \(g(t_2)=\left[\begin{matrix}\frac{X(t_2)-x(t_1)}{r(t_1)}&\frac{Y(t_2)-y(t_1)}{r(t_1)}&\frac{Z(t_2)-z(t_1)}{r(t_1)}\\\end{matrix}\right]\) and \(\delta u(t_2)=\left[\begin{matrix}\delta x(t_2)&\delta y(t_2)&\delta z(t_2)\\\end{matrix}\right]\). Based on Eqs. (2) and (3), Eq. (4) is obtained.

\(\Phi\left(t_{12}\right)=S(t2)-u(t1)-S(t1)-u(t1)+g(t2)⋅δu(t2)+cBt12+wt12\)                                                                    (4)

when \(\Delta S(t_2)=S(t2)-u(t1)-S(t1)-u(t1)\), which can be calculated, is transposed to the left-hand side, it can be expressed as shown in Eq. (5).

\(\Phi\left(t_{12}\right)-\Delta S\left(t_2\right)=g\left(t_2\right)\cdot\delta u\left(t_2\right)+cB\left(t_{12}\right)+w^i\left(t_{12}\right)\)                                                                                (5)

The unknowns to be calculated are the change in the receiver position, \(\delta u(t_2)=\left[\begin{matrix}\delta x(t_2)&\delta y(t_2)&\delta z(t_2)\\\end{matrix}\right]\), and the change in the receiver clock error, \(cB\left(t_{12}\right)\). For the calculation of the four unknowns, a linearized time differencing equation such as Eq. (5) needs to be made using the carrier phase measurements of more than or equal to four satellites, in order to obtain the solution. The change in the receiver clock error is identical to the clock drift; and in this study, they were interchangeably used \((cB\left(t_{12}\right)=c\dot{B}\left(t_2\right))\).

The linearized time differencing equation for \(n(\geq4)\) carrier phase measurements can be expressed by Eq. (6). In this regard, the superscript represents the satellite number. For the ease of expression, Eq. (6) can be simply expressed by a vector-matrix equation, as shown in Eq. (7).

\(\left[\begin{matrix}\Phi^1\left(t_{12}\right)-\Delta S^1(t_2)\\\Phi^2\left(t_{12}\right)-\Delta S^2(t_2)\\\vdots\\\Phi^n\left(t_{12}\right)-\Delta S^n(t_2)\\\end{matrix}\right]=\left[\begin{matrix}g^1(t_2)&1\\g^2(t_2)&1\\\vdots&\vdots\\g^i(t_2)&1\\\end{matrix}\right]\left[\begin{matrix}\delta u(t_2)\\cB\left(t_{12}\right)\\\end{matrix}\right]+\left[\begin{matrix}w^1\left(t_{12}\right)\\w^2\left(t_{12}\right)\\\vdots\\w^n\left(t_{12}\right)\\\end{matrix}\right]\)                                                                          (6)

\(\Phi_{TD}=G\delta x_{TD}+w_{TD}\)                                                                                                                (7)

The TDCP measurement can be calculated by Eq. (8) using \(T_D\) operator.

\(\left[\begin{matrix}\Phi^1\left(t_{12}\right)\\\Phi^2\left(t_{12}\right)\\\vdots\\\Phi^n\left(t_{12}\right)\\\end{matrix}\right]=\left[\begin{matrix}-1&1&0&0&\cdots&0&0\\0&0&-1&1&\cdots&0&0\\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&0&\cdots&-1&1\\\end{matrix}\right]\left[\begin{matrix}\Phi^1\left(t_1\right)\\\Phi^1\left(t_2\right)\\\Phi^2\left(t_1\right)\\\Phi^2\left(t_2\right)\\\vdots\\\Phi^n\left(t_1\right)\\\Phi^n\left(t_2\right)\\\end{matrix}\right]=T_D\left[\begin{matrix}\Phi^1\left(t_1\right)\\\Phi^1\left(t_2\right)\\\Phi^2\left(t_1\right)\\\Phi^2\left(t_2\right)\\\vdots\\\Phi^n\left(t_1\right)\\\Phi^n\left(t_2\right)\\\end{matrix}\right]\)                                                                (8)

Based on Eq. (8), the covariance of the TDCP measurement is expressed as shown in Eq. (9), and this indicates that the variance of the error doubles due to time differencing.

\(\textrm{cov}{(}w_{TD})=\sigma_\Phi^2(T_DT_D^T)=2\sigma_\Phi^2I=Q_{TD}\)                                                                                        (9)

when the weighted least squares method is applied to the measurement of Eq. (7), the solution of Eq. (10) can be obtained.

\(\delta x_{TD}=(G^TQ_{TD}^{-1}G)^{-1}G^TQ_{TD}^{-1}\Phi_{TD}\)                                                                                          (10)

By adding the change in the receiver position, \(\delta u(t_2)\), obtained in Eq. (10) to the receiver position at the previous time, \(u(t_1)\), which is the linearization reference point, the current receiver position, \(u(t_2)\), is calculated as shown in Eq. (11).

\(u(t_2)=u(t_1)+\delta u(t_2)\)                                                                                                  (11)

The heave is motion in the altitude direction, and the altitude value is obtained by converting the value obtained in the Earth Centered Earth Fixed coordinate system into a value in the latitude/longitude/altitude coordinate system. As shown in Eq. (12), the heave is determined as the relative altitude relative to the initial altitude by subtracting the initial altitude, \(alt(0)\), from the current altitude, \(alt(t_2)\).

\(h(t_2)=alt(t_2)-alt(0)\)                                                                                                  (12)

3. IMPLEMENTATION AND TEST OF THE HEAVE DETERMINATION SYSTEM

3.1 Implementation of the Heave Determination System

The heave determination system consists of hardware and Matlab program, as shown in Fig. 1. For the hardware, ublox EVK-M8T receiver and PC were used. The ublox EVK-M8T receiver outputs the raw pseudorange measurement message (UBX-RXM-RAWX) and raw satellite ephemeris message (UBX-RXM-SFRBX) at a maximum of 4 Hz (ublox 2017). This was received in the PC using serial communication.

Fig. 1. Block diagram of heave determination system.

The Matlab program can be divided into the data processing unit and the heave determination unit. The data processing unit stores necessary information from the raw measurement and navigation message obtained from the receiver, and calculates the elevation angle of the satellite based on the satellite position and the position calculated using the code measurement of the receiver. When the calculated elevation angle is lower than the reference value, it is eliminated and not used. When the CN0 and receiver tracking status (trkstat) provided by the receiver are below the reference values, the corresponding satellite is also excluded from the navigation. The receiver tracking status (trkstat) detects the validity of code measurement, the validity of carrier phase measurement, half cycle slip, and subhalf cycle slip, so that a user will not use wrong measurement (ublox 2017). In many experiments using the ublox EVK-M8T receiver, cycle slip did not occur when wrong measurements were eliminated using the receiver tracking status. Therefore, in this study, an additional cycle slip detection technique was not considered.

The heave determination unit calculates the navigation solution using the TDCP method derived in Chapter 2. The initial position is obtained based on the absolute positioning technique; and for the subsequent measurements, the heave is calculated using the TDCP measurement.

3.2 Performance and Problem of the Implemented Heave Measurement System

Using the ublox EVK-M8T receiver, measurements in a static condition were collected for 54 minutes from 13:00, July 09, 2017 at the building of the School of Electronics Engineering, Chungbuk National University. Fig. 2 shows the satellite arrangement during the experiment period, the number of observed satellites, and the Position Dilution of Precision (PDOP).

Fig. 2. Skyplot and number of satellite, PDOP.

As shown in Fig. 3, to evaluate the performance of the proposed method, the heave values were calculated based on: 1) Absolute positioning technique using code measurement, 2) Absolute positioning technique using carrier smoothed code measurement, and 3) Positioning technique using TDCP measurement. For all the three techniques, the initially obtained altitude was used as the reference value, and the altitude difference was expressed as the heave. Methods 1) and 2) showed similar results because they calculated the solution based on the absolute positioning technique, though different measurements were used. However, for the result of the carrier smoothed code measurement, the obtained position was smooth and the amplitude of the error decreased due to the effect of the carrier measurement. On the other hand, the result of the time differencing technique Method 3) was different from those of Methods 1) and 2). When the TDCP measurement was used, the altitude drifted as time passed. This drift error is due to the accumulation of the error included in the estimated displacement during the process of obtaining the current position by accumulating the displacement estimated using the TDCP method onto the previous position. The cause for the error of the estimated displacement is thought to be the residual ionospheric/tropospheric error had not been removed by the time differencing and multipath error that had not been modeled in Eq. (2) Except for the drift error, it is expected that more precise result could be obtained compared to the other two methods.

Fig. 3. Comparison of heave determination method (Code absolute, Carrier smoothed code, TDCP).

4. PERFORMANCE IMPROVEMENT USING HPF

4.1 Design of HPF

As shown in the experiment described in Chapter 3, the altitude value that had been calculated using the time differencing drifted as the time passed. It was reported that this drift error could be effectively resolved using HPF (BEng 2007). However, BEng (2007) designed 276th-order FIR HPF using the Matlab filter design toolbox, and there are problems of complicated calculation and delayed output in practical implementation.

In this study, considering the implementation of a heave determination system using a low-cost microprocessor, first-order IIR HPF was designed as shown in Eq. (13). In this regard, \(h(k)\) is the heave obtained in Eq. (12), \(a\) is the weighting factor, \(x(k)\) is the internal state variable, and \(h_F(k)\) is the output of the HPF. The zero point of the designed HPF is located at 1, and the pole point is located at \(a\). As the \(a\) value decreases, the cutoff frequency of the HPF increases.

\(x(k+1)=ax(k)+(1-a)h(k)\)                                                                                             
\( h_F(k)=h(k)-x(k)\)                                                                                                     (13)

The first-order HPF designed in this study can be more easily implemented and has less output delay compared to the complicated HPF designed in BEng (2007). However, the performance is significantly affected by the accurate setting of the cutoff frequency and the gain of the passband. BEng (2007) used a cutoff frequency of 0.03 Hz for the design of the filter. Based on this value, the cutoff frequency was finally determined through an experiment in this study.

4.2 Experiment and Performance Evaluation

To evaluate the performance of the heave measurement system based on the time-differenced carrier measurement technique using HPF and to examine the elimination of drift error as a result of HPF application, two types of experiments were conducted: a dynamic experiment where the antenna was attached to an experiment apparatus moving up and down at a cycle of 20 seconds and an amplitude of 48 cm, and a static experiment where the antenna was fixed.

4.2.1 Dynamic experiment

To examine the heave measurement performance of the heave determination system using HPF in a dynamic environment, an experiment was conducted by attaching the antenna to an experiment apparatus moving up and down at a cycle of 20 seconds and an amplitude of 48 cm. The experiment was carried out on the rooftop of the laboratory of Samyung ENC Co., Ltd. GPS raw measurements and navigation message were collected using the ublox EVK-M8T receiver at 16:00, February 6, 2017, and the heave was measured through post-processing. Fig. 4 shows the satellite position at the time of the collection. The number of observed GPS satellites was 7~8. The PDOP was around 2 when the number of observed satellites was 8, but it was close to 4 when the number of observed satellites was 7.

Fig. 4. Sky plot and number of satellite, PDOP.

Fig. 5 shows the result without the application of HPF and the result with the application of the 276th-order HPF in BEng (2007). As shown in Fig. 5, when only the TDCP measurement was used without the application of HPF, motion with an amplitude of 48 cm was observed, but it drifted as time passed. On the other hand, when the 276th-order HPF was used, the drift disappeared, and motion with an amplitude of 48 cm centered at 0 was clearly observed. In Fig. 5, the peak values were not constant because 1) the collection cycle of the measurement and the peak of the motion were not constant and 2) the apparatus that moves the antenna up and down using a motor was inaccurate. Considering this, it was shown that accurate heave could be measured using the 276th-order HPF. However, when the 276th-order HPF is used, the output is delayed by 276 data, and thus it cannot be applied to a system requiring real-time heave information. In the case of the first-order IIR HPF designed in this study, the output is delayed by one data, and thus it is advantageous for real-time applications. To determine the characteristics of the designed first-order HPF, the coefficient of the filter, \(a\), needs to be determined. Fig. 6 shows the results of the experiments using different \(a\) values. When \(a\)=0.99, the cutoff frequency was close to 0 and thus the effect of HPF was barely observed. On the other hand, when \(a\)=0.8, the drift disappeared, but the peak value decreased as the gain of the filter decreased. It was slightly improved when \(a\)=0.9, and the most desirable output was obtained when \(a\)=0.95. Similar results were obtained in other experiments; and based on this, it was fixed at \(a\)=0.95 in this study. When =0.95, the cutoff frequency was around 0.008 Hz, which was slightly lower than the value (0.03 Hz) used in BEng (2007).

Fig. 5. Result of dynamic experiment (HPF not applied vs. 276th order HPF).

Fig. 6. Heave result according to change of coefficient (Dynamic experiment).

Fig. 7a shows the result of the dynamic experiment when \(a\)=0.95. In Fig. 7, the blue color represents the heave measurement before the application of the HPF, and the red color represents the heave estimate after the application of the HPF. The result showed that the drift error before the application of the HPF was eliminated by the application of the HPF. In addition, a delay of 1 epoch occurred due to the application of the first-order HPF. Fig. 7b shows the result of the first-order HPF when \(a\)=0.95, and the result of the 276th-order HPF in BEng (2007). In Fig. 7, the blue color represents the heave estimate based on the 276th-order HPF of BEng (2007), and the red color represents the heave estimate based on the first-order HPF. For the 276th-order HPF, a delay of 276 epochs was observed; but in the case of the first-order HPF, output was provided after 1 epoch. In addition, as the outputs of the two HPFs were consistent within 4 cm, the first-order HPF was found to be advantageous for real-time applications.

Fig. 7. Comparison of heave result (Dynamic experiment).

4.2.2 Static experiment

To examine the precision and drift error elimination performance of the heave measurement system using the HPF, the measurement used in Chapter 3 was used. In this regard, \(a\)=0.95, and the results are shown in Fig. 8. In Fig. 8, the blue color represents the heave estimate before the application of the HPF, and the red color represents the heave estimate after the application of the HPF. In the static condition, the designed HPF also operated efficiently, and the obtained heave had a mean of -4.13 cm and a standard deviation of 2.06 cm, indicating that a very precise result could be obtained.

Fig.8. Heave result of static experiment (1st order HPF applied vs. not applied).

5. CONCLUSION

In this study, a system that calculates precise altitude using the TDCP method and determines heave based on the obtained altitude was implemented using a commercial receiver. In the case of the heave value determined by the TDCP method, drift errors occur due to the accumulation of the error. To resolve the problem of drift error occurrence, HPF was additionally applied. The HPF was designed as a simple first-order IIR filter considering the ease of real-time implementation; and to minimize the performance deterioration, the coefficient was selected through an experiment so that the cutoff frequency and passband gain could be similar to those of the 276th-order FIR filter. To evaluate the performance of the proposed system, a real-time heave determination system based on GNSS was implemented using ublox EVK-M8T receiver and Matlab, and the performance was examined through an actual experiment. In the case of the static experiment, a very precise result was obtained where the mean was -4.13 cm and the standard deviation was 2.06 cm. In the case of the up and down motion with a 48 cm amplitude, a precise result was also provided without a delay compared to the 276th-order filter.

Based on the technique proposed in this study, a precise heave measurement system can be implemented using only one GNSS receiver without the help of a reference station. Thus, it can be implemented through the addition of only software without additional hardware installation costs if an existing receiver installed for navigation purposes provides raw pseudorange measurement and navigation message. Therefore, the proposed technique can be easily used for various applications requiring heave measurement, such as the altitude correction of fish detection radar.