Medium to Long Range Kinematic GPS Positioning with Position-Velocity-Acceleration Model Using Multiple Reference Stations

In order to obtain precise kinematic global positioning systems (GPS) in medium to large scale networks, the atmospheric effects from tropospheric and ionospheric delays need to be properly modeled and estimated. It is also preferable to use multiple reference stations to improve the reliability of the solutions. In this study, GPS kinematic positioning algorithms are developed for the medium to large-scale network based on the position-velocity-acceleration model. Hence, the algorithm can perform even in cases where the near-constant velocity assumption does not hold. In addition, the estimated kinematic accelerations can be used for the airborne gravimetry. The proposed algorithms are implemented using Kalman filter and are applied to the in situ airborne GPS data. The performance of the proposed algorithms is validated by analyzing and comparing the results with those from reference values. The results show that reliable and comparable solutions in both position and kinematic acceleration levels can be obtained using the proposed algorithms.


Introduction
For the past two decades, global positioning systems (GPS) has been widely used for precise positioning in various engineering and scientific fields. While its application areas in the early stage were limited to the determination of a position at a static environment, more and more applications which require kinematic positioning have been developed. For the kinematic GPS positioning, either single or multi-reference stations can be used for the data processing. The multi-reference station approach generally offers better positioning accuracy and reliability [1,2]. Therefore, the multi-reference station approach is recommended when precise kinematic positioning is required. However, many distance-dependent errors, such as atmospheric and ionospheric delays, are not fully removed through the double-differencing technique when the baseline length increases. Thus, these errors should be properly modeled for the medium to long range kinematic positioning [3][4][5]. GPS kinematic positioning is usually performed by using the Kalman filter, and so called Position-Velocity (PV) model can be adopted when the object is assumed to move with a nearly constant velocity. In addition, Position-Velocity-Acceleration (PVA) model can be introduced for the cases where the near-constant velocity assumption does not hold [6]. This means that the PVA model can provide not only the flexibility in estimation procedure but also the position, velocity, and kinematic acceleration of the moving object simultaneously.
It is notable that one of the popular applications utilizing the GPS derived acceleration is the airborne gravimetry. The airborne gravity survey is performed by measuring the kinematic accelerations of the aircraft using GPS. Then, the gravity information is extracted using the well-known equation as shown in Equation (1). The Equation (1) expressed in inertial frame shows that the kinematic accelerations consist of the gravitational acceleration and the specific force sensed by an accelerometer [7,8]; x g a = +  (1) where x  is kinematic acceleration; g is gravitational acceleration; a is specific force.
In general, the GPS measurements are processed in relative positioning mode to obtain accurate and reliable positioning results in 3-D space. Then, the kinematic accelerations are computed by taking second-order time-derivatives of positions [9,10]. This means that two steps of data processing are required to get the kinematic accelerations from the positions, which decrease the efficiency in a computational point of view. To overcome the drawback of this method, one can also use the PVA model approach so that the position, velocity, and acceleration of the aircraft can be determined at the same time. An efficient PVA model for Kalman filter is proposed and demonstrated using the measurements from the electronic tacheometer, i.e., geodetic kinematic measurements [11]. However, it should be noted that the experiments are conducted in an indoor environment, and the measurement used in the experiments are not from the GPS.
In this paper, the PVA model to directly determine not only the positions but also the kinematic accelerations is proposed, and the validation of the algorithm is demonstrated with in situ GPS measurements. The proposed algorithm is designed for the medium to long baseline scenario so that it can perform properly even in a survey for a large area under unfavorable surveying condition. The atmospheric effects such as tropospheric and ionospheric delays are modeled in the Kalman filter for the medium to long baseline scenario. In addition, multi-baseline approach is adopted for the improvement of the accuracy and reliability of the solutions. The algorithms are developed and numerical test are performed to validate the proposed algorithms by analyzing the positioning results and comparing the kinematic accelerations with those from the second-order time-derivatives of positions.

Methodology
The Kalman filter is an optimal recursive estimator which uses a system's dynamics model and sequential measurements to estimate the system's unknown states in a minimum variance sense. The Equation (2) shows the discrete Kalman filter equations categorized into two components, i.e., the prediction of state vector through the state-transition matrix and the update of the state vector using the measurements [12,13].
where subscript k indicates the epoch; x k is state vector; Φ k is state-transition matrix; w k is process noise vector; Q k is covariance matrix for the processes noise; z k is measurement vector; H k is design matrix; v k is measurement noise vector; R k is covariance matrix for the measurement noise.

State Vector and System Equation
The objective of the proposed algorithm is to determine the kinematic accelerations of the aircraft in unfavorable conditions. Therefore, the algorithm is designed for the medium to long baseline scenario implementing multi-baseline data processing. This means that the proper modeling for the atmospheric effects such as tropospheric and ionospheric delays is required and, as a consequence, the zenith wet delay (ZWD) residuals and the double-differenced (DD) ionospheric delays are included in the state vector as shown in Equation , 1 where subscript i and j denote reference station and the moving aircraft, respectively; superscript k and  are the satellite indices; x , y , z are the position components; x  , y  , z  and The coordinates of reference stations are included in the state vector and estimated with an assumption of random constant stochastic process, so the transition matrix for those are given as in Equation (5).
The state-transition matrix for position-velocity-acceleration components of the aircraft can be derived with the assumption that the kinematic acceleration is Gaussian white noise process [6,11,14]: For the modeling of tropospheric effect, either random walk or first-order Gauss-Markov process is used because either process is known to be sufficient to describe the stochastic property of tropospheric effect [4]. Here, the zenith wet delay residual, i.e., increment with respect to a priori value, is assumed to be random walk stochastic process and the transition matrix is: The remaining atmospheric effect, i.e., DD ionospheric delays, is modeled with first-order Gauss-Markov process. The first-order Gauss-Markov process can describe a large number of physical processes with reasonable accuracy [15]. However, it should be noted that the DD ionospheric delays are correlated with not only time but also baseline length. Hence, the DD ionospheric delays are modeled with the first-order Gauss-Markov process which incorporates the effect of both the time and baseline length changes as shown by [3].
where T is first-order correlation time for the DD ionosphere; δ is baseline length change; 1500 km D = is the first-order correlation distance.
The DD ambiguities for L1 and wide-lane measurements can be modeled as a random constant process as the case of reference stations' coordinates. It should be noted that the ambiguity resolution is not attempted because it is an issue outside the scope of this study in long range static or kinematic positioning.
For the Kalman filtering, the process noise covariance matrix should be properly modeled as it confines the variability of the state estimates. The discrete form of covariance matrices for the reference stations' coordinates and the position-velocity-acceleration of aircraft can be written as shown in Equations (9) and (10), respectively [6,11].
The initial spectral density, acc q for the acceleration states can be determined using the initial kinematic solutions by computing the standard deviation of the time-differenced accelerations [14]. In this study, acc q is set to 1 m 2 /s 4 /Hz for the data processing through the analysis of the spectral density from the estimated acceleration data. As mentioned already, the zenith wet delay residual is modeled with a random walk process so that the corresponding process noise covariance matrix can be expressed as following: It should be noted that the spectral density, ZWD q Δ may be assigned with different values for reference stations and aircraft because the behavior of the zenith wet delay residuals are expected to be different from each other. The spectral density for zenith wet delay, ZWD q Δ are determined by using the empirical auto-covariance functions method proposed by [4]. Then, 3 × 10 −10 m 2 /Hz and 5 × 10 −8 m 2 /Hz are chosen for the spectral densities for reference stations and aircraft, respectively. The process noise covariance matrix for DD ambiguities for L1 and wide-lane is modeled with zero matrix. The covariance matrix for DD ionospheric delay is provided by [3] as shown in Equation (12).

GPS Measurement Equations
The DD GPS measurements are widely used for the relative positioning due to the fact that various systematic errors can be eliminated by differencing operation [16]. In other words, many error sources in the GPS measurements, for example, receiver and satellite clock errors, and atmospheric effects can be removed or significantly reduced in the DD measurements. However, the spatially correlated atmospheric errors should be modeled as the spatial separation between the two receivers increases.
The Equation (13) shows the DD measurement equations for both L1 and L2 frequencies [17,18] where, subscript i and j indicate receivers while superscript k and  denote satellites; ( ) measurements on L1 and L2, respectively; 1 P , 2 P are code pseudo-range measurements on L1 and L2, respectively; ρ is the geometric distance between the receiver and satellite; T is tropospheric delay, 1 f and 2 f are carrier frequencies of L1 and L2, respectively; The DD GPS measurement equations shown in Equation (13) are nonlinear forms; thus, the linearization step is required. After linearization of Equation (13), the design matrix, H has the matrix form as shown in Equation (14).
As shown in Equation (14), the components of the design matrix can be categorized into partial derivatives with respect to the coordinates for i (reference station) and j (aircraft), ZWD residual, DD ambiguities, and DD ionospheric delays.
Hence, the components corresponding to the coordinates of reference station and rover in the design matrix can be expressed as following: , , The tropospheric delay consists of dry and wet components. Approximately 90% of the total zenith delay (TZD) is from the dry component, and about 10% is from the wet component [16]. It is well known that the dry component can be modeled with high accuracy while the wet component is less predictable. Therefore, the ZWD residual with respect to a priori value is estimated in this approach. Equation (16) shows the tropospheric delays between the receiver and satellite.
where ZDD is modeled zenith dry delay; ZWD is modeled zenith wet delay; ZWD Δ is ZWD residual; m is mapping function; subscript D and W indicate "dry" and "wet", respectively.
The modeled ZDD , ZWD and their mapping functions are computed by using the Saastamoinen model [19]. The tropospheric delays in Equation (16) can be used to make the DD forms and the components of the design matrix correspond to the ZWD Δ is as following: To obtain the covariance matrix for the DD measurements, R k , the error propagation law is applied to the covariance matrix for the un-difference (UD) measurements. The magnitudes of the noises for code and phase pseudo-range observations are set to decimetre-level and milimetre-level, respectively [14].

Results and Discussion
An airborne gravity survey was conducted in South Korea to develop a new geoid model in 2009. A Cessna Grand Caravan was used for the survey and it was flown with the help of autopilot at a constant altitude, i.e., 10,000 feet, in order to get as smooth a flight as possible. The GPS data were collected from both the GPS receiver on board and six Continuously Operating Reference Station (CORS) with 1 s data interval. Figure 1 shows one of the trajectories of aircraft (blue line) and the locations of reference CORS used to demonstrate the proposed algorithm. The trajectories of the aircraft started and ended at the Gimpo (GIMP) airport and the total length of the aircraft trajectory is about 1200 km. The coordinates of the reference stations in ITRF2005 frame are determined using the BERNESE GPS data processing software [20], and the results are shown in Table 1. The coordinates of reference stations are tightly constrained by assigning very small variance values (i.e., 1 × 10 −10 m 2 ) to the initial covariance matrix.
For the medium to long baseline data processing scenario, the SHAO reference station is included in the data processing. Hence, the diameter of the CORS GPS network reaches up to approximately 900 km. The GPS data was collected on 11 January 2009 and the data span was 4 h and 18 min.  As described in Section 2, the kinematic accelerations together with positions and velocities are estimated through the Kalman filter approach. The multi-baseline data processing is performed with radial type of network configuration. To evaluate the filter efficiency, the convergence of the covariance matrix Pk is examined. Figure 2a presents the variations of the trace of Pk with respect to time. As shown in Figure 2a, the trace of Pk converges rapidly with respect to time, which indicates the reliability of the filter. It is also notable that two peaks in Figure 2a correspond to the epochs at which new satellites are observed. Then, new states for DD ionospheric delays and ambiguities are included in state vector, and predefined variance values, i.e., (0.5 m) 2 for DD ionospheric delay and 1 × 10 10 for DD ambiguity, are assigned to the new states.     Figure 3b, the constant velocities in most of the data span for airborne gravity survey are retained. The peaks observed in Figure 3c correspond to the aircraft maneuvering for takeoff, turnaround, and landing. Table 2 shows the statistics of the estimated accelerations which ranges from about −4.5 m/s 2 to 4.5 m/s 2 . The atmospheric effects are estimated in forms of ZWD residuals and DD ionospheric delays, respectively. As explained in Section 2, the ZWD residuals are estimated at each station including the aircraft and then final TZDs are computed by adding the modeled values to the estimated ZWD residuals. Figure 4a shows the final TZDs obtained from the data processing. As shown in Figure 4a    presents the differences between the estimated coordinates using the proposed method and the reference values from the Waypoint software. The differences between the two solutions range from about −0.15 m to 0.1 m as shown in Table 3. From the results, one can conclude that the filter based on the proposed algorithm is working properly.  Then, the second-order time-derivatives of determined positions are applied to obtain the kinematic accelerations. It should be mentioned that taking the time-derivatives is performed after the data fitting with B-spline function. In the final step, the differences between the two kinematic accelerations are computed as shown in Figure 6. This procedure is applied to only 1 h data span, i.e., 1-2 h (UTC) because no significant dynamics of the aircraft is observed. This also indicates that the actual gravity survey is performed during that time span. Table 4 shows the statistical characteristics of the differences between the estimated kinematic accelerations and the reference values. The differences between the two solutions range from −0.15 m/s 2 to 0.16 m/s 2 and the standard deviation of the differences is about 0.03 m/s 2 . The mean value of the differences is 10 −5 m/s 2 level for each component.

As shown in
From the results, one can state that the comparable kinematic accelerations for airborne gravity survey can be obtained by using the proposed algorithmIt should be mentioned that the kinematic acceleration for airborne gravimetry needs the accuracy at least at the level of 10 −4 m/s 2 to detect the gravity signal at the level of 10 mGal. In airborne gravimetry, however, a smoother is applied to the acceleration to extract meaningful gravity signal. Although the extraction of the gravity signal is out of the scope in this study, a B-spline smoother with window size of 120 s are applied to both accelerations and verified that the differences reside at the level of 10 −4 m/s 2 . Therefore, one can conclude that the proposed algorithm is reliable for the application of kinematic accelerations. The validation for the proposed algorithm in this study, however, is performed with the dataset of one trajectory collected under the relatively low dynamic conditions. Therefore, it might be necessary to process more dataset under the high dynamic environments in the future.

Conclusions
In this study, we proposed the kinematic GPS positioning algorithm using multiple reference stations for medium to large scale networks. The Kalman filter with PVA dynamic model is used for the kinematic positioning so that the kinematic acceleration information can be obtained simultaneously. For the long range kinematic applications, the tropospheric and ionospheric delays are modeled with random walk and first-order Gauss-Markov process models, respectively. The algorithm is implemented and tested with in situ airborne GPS data collected at 10,000 ft altitude with 1 s data interval. For the validation of the proposed algorithm, the positioning results are compared with the reference values first. The computed differences between the estimated and reference values range from −0.17 m to 0.11 m, which indicates the validity of the proposed algorithm. The estimated kinematic accelerations from the proposed algorithms are also compared with those from the second-order time-derivatives of the reference positions. The computed differences between the two solutions range from −0.15 m/s 2 to 0.16 m/s 2 and the standard deviation of the differences is about 0.03 m/s 2 . From these results, we can conclude that comparable kinematic accelerations, which can be used for airborne gravimetry, are obtainable using the proposed algorithm.