A New Strategy for Improving the Accuracy of Aircraft Positioning Using DGPS Technique in Aerial Navigation

Abstract

In this paper a new mathematical algorithm is proposed to improve the accuracy of DGPS (Differential GPS) positioning using several GNSS (Global Navigation Satellites System) reference stations. The new mathematical algorithm is based on a weighting scheme for the following three criteria: weighting in function of baseline (vector) length, weighting in function of vector length error and weighting in function of the number of tracked GPS (Global Positioning System) satellites for a single baseline. The algorithm of the test method takes into account the linear combination of the weighting coefficients and relates the position errors determined for single baselines. The calculation uses a weighting scheme for three independent baselines denoted as $\left(1A,2A,3A\right)$. The proposed research method makes it possible to determine the resultant position errors for ellipsoidal BLh coordinates of the aircraft and significantly improves the accuracy of DGPS positioning. The analysis and evaluation of the new research methodology was checked for data from two flight experiments carried out in Mielec and Dęblin. Based on the calculations performed, it was found that in the flight experiment in Mielec, due to the application of the new research methodology, DGPS positioning accuracy improved from 55 to 94% for all the BLh components. In turn, in the flight experiment in Dęblin, the accuracy of DGPS positioning improved by 63–91%. The study shows that the highest DGPS positioning accuracy is seen when using weighting criterion II, the inverse of the square of the vector length error.

1. Introduction

The GNSS (Global Navigation Satellites System), as a satellite technology, is continuously developed and widely used in many areas of human life, including the entire aviation industry [1]. However, it is important to note that the benefits of the GNSS in aviation are enormous, ranging from improved flight safety to reduced fuel emissions for example for RNAV (Area Navigation) [2] or the economic factor itself through improved airport capacity and an increased number of flight operations [3]. This is especially important as the GNSS, as one of many on-board sensors, is now a fundamental navigation sensor on board every aircraft. This makes the implementation of this sensor in aircraft positioning and the selection of an appropriate positioning method or technique particularly important for aviation science. GNSS satellite technology allows the user to determine the position of the aircraft using absolute and differential methods [4]. Within the presented work, the Differential DGPS (Differential GPS) method [5] in the GPS (Global Positioning System) satellite system will be discussed. In the differential DGPS technique, the location of the base stations relative to the mobile on-board GPS receiver mounted in the aircraft plays a key role [6]. Moreover, the number of GNSS reference stations involved in DGPS positioning is also very important [7]. The calculation of the aircraft position in DGPS technique usually uses code observations at the L1 frequency of the GPS system, to ensure that it is possible to compare the pseudoranges between the GNSS reference station and the on-board GPS receiver [8]. This comparison of pseudoranges for the same instant of time has been called differential satellite observations. The differential process itself allows for the reduction in many systematic errors such as satellite clock error, hardware TGD (Timing Group Delay) or relativistic correction [9]. On the other hand, the influence of atmospheric delays is reduced thanks to pseudorange differentiation, which, consequently, also removes the ionospheric and tropospheric correction from the observation equation [10]. A very important parameter, i.e., receiver clock error is compensated in the observation equation in the DGPS method [11]. In fact, the compensation results in the determination of the receiver clock error difference from the observation equation in the DGPS method. As regards the position navigation parameters, the DGPS method determines three position components in the form of geocentric XYZ or ellipsoidal BLh coordinates [8].

2. Related Papers

The use of the differential DGPS technique in air navigation is crucial for at least several reasons. Many applications of the DGPS method in air navigation can be found in the research literature. In particular, the application of the DGPS technology was evident in the following research areas:

-

the use of the DGPS method as the primary method for aircraft positioning [6,12,13],

-

the use of the DGPS method to determine the reference position of the aircraft trajectory [14,15],

-

the use of the DGPS method as part of the GBAS (Ground Based Augmentation System) and the PA (Precision Approach) [16,17,18,19,20,21,22,23,24],

-

the application of the DGPS method in mapping applications, cartography and photogrammetry [25,26,27,28],

-

the application of the DGPS method in monitoring the ionosphere [29],

-

the use of the DGPS method in integration with various positioning sensors to determine flight altitude [30,31],

-

the determination of the accuracy and precision of DGPS positioning in air navigation [8,32,33,34],

-

the determination of the integrity of positioning using the DGPS technique [34,35],

-

the integration of DGPS and DGLONASS (Differential GLONASS) solutions [6,13,36].

As shown in the research literature, the use of the DGPS technique in air navigation has been thematically extensive. For navigation, the most relevant is the development and elaboration of new mathematical algorithms to improve the accuracy of aircraft positioning, also using the DGPS technique. Therefore, the development of a new position navigation solution can significantly contribute to the improvement of the performance of the determined coordinates of the aircraft derived from the DGPS solution. It may also result in the improved safety of flight and flight operations. Consequently, the main idea of the paper is to develop an optimal algorithm to improve the accuracy of DGPS positioning in aviation using GNSS reference stations. The publication shows a computational strategy for determining the accuracy of the resultant aircraft position using a DGPS solution with three GNSS reference stations. The concept of the algorithm is based on the use of a weighted average model for three weighting criteria. Namely, in the research, a weighting scheme was developed as a function of the baseline (vector) length, the vector length error and the number of tracked GPS satellites. Such a weighting strategy is important for several factors. First, it relates the geometry of the distribution of GNSS reference stations and the flight position of the aircraft by calculating the vector length. Secondly, the weighting scheme is also dependent on the number of GPS satellites taken for the DGPS solution for a single baseline. Finally, the weighting process uses the vector length error values, which are actually the determined average coordinate errors of the stochastic model from the DGPS solution for a single baseline. One can notice that the choice of weighting parameters is not random; therefore, indicating which weighting factor is optimal for improving the accuracy of DGPS positioning for several baselines will be the main task to be carried out in the submitted scientific paper. The obtained research results will only confirm which weighting variant is the most optimal and best for DGPS positioning using several GNSS reference stations. This will also allow further research into the development of modern and advanced algorithms optimising the mathematical model of DGPS positioning in air navigation.

In summary, our main research achievements in the paper relate to the following:

-

development of weighted average model algorithms to improve the accuracy of DGPS positioning using several GNSS reference stations,

-

implementation of a weighting model for the following three criteria: weighting as a function of baseline (vector) length, weighting as a function of vector length error, weighting as a function of the number of tracked GPS satellites,

-

definition of algorithms to calculate the resultant position errors,

-

conducting a formal evaluation of the research results obtained in percentage terms.

3. Research Method

Section 3 demonstrates the basic algorithm of the DGPS positioning method for a single baseline and a new mathematical model to improve the accuracy of DGPS positioning, including measurement weighting.

3.1. Basic DGPS Solution for a Single Baseline

The basic observation equation for the differential DGPS technique in the GPS navigation system for air navigation can be presented as follows [37,38,39]:

$$\Delta P_{1A}^{}={\rho }_{1A}^{}+\Delta I_{1A}^{}+\Delta T_{1A}^{}+\Delta M_{1A}$$

(1)

where $\Delta$—the operator of a single difference for GPS code measurements, it allows determining the difference in code measurements for two satellites tracked by one receiver; $\overrightarrow{1A}$—the vector in the space between the GNSS base station ($\overrightarrow{1}$) and the GPS rover receiver ($\overrightarrow{A}$) mounted onboard the aircraft; $P_{1A}^{}$—the value of single code difference (expressed in meters) on the vector $\overrightarrow{1A}$ at the L1 frequency in the GPS system; ${\rho }_{1A}^{}$—the geometric distance of the vector $\overrightarrow{1A}$ for the single code difference (expressed in geocentric coordinates XYZ) in the GPS system; ${\rho }_{1A}^{}={\rho }_1-{\rho }_A^{}$, ${\rho }_1^{}$—the geometric distance between the GNSS base station ($\overrightarrow{1}$) and GPS satellite; ${\rho }_1=\sqrt{{(X_1-X_{sat})}^2+{(Y_1-Y_{sat})}^2+{(Z_1-Z_{sat})}^2}$; $(X_1,Y_1,Z_1)$—the geocentric coordinates of the GNSS base station ($\overrightarrow{1}$); $(X_{sat},Y_{sat},Z_{sat})$—the geocentric coordinates of GPS satellite; ${\rho }_A^{}$—the geometric distance between the GPS rover receiver ($\overrightarrow{A}$) mounted onboard the aircraft and GPS satellite, ${\rho }_A=\sqrt{{(X_A-X_{sat})}^2+{(Y_A-Y_{sat})}^2+{(Z_A-Z_{sat})}^2}$; $(X_A,Y_A,Z_A)$—the geocentric coordinates of aircraft vehicle, e.g., GPS rover receiver ($\overrightarrow{A}$) mounted onboard the aircraft; $I_{1A}^{}$—the value of ionosphere delay at L1 frequency for single code difference in the GPS system; $T_{1A}$—the value of troposphere delay for single code difference in the GPS system; $M_{1A}$—the multipath effect and noise measurement at an L1 frequency for the GPS code measurements.

Equation (1) describes the mathematical model for DGPS positioning in air navigation for a single baseline. For this case, a single difference for GPS code observations on a vector ($\overrightarrow{1A}$) between a GNSS reference station and an on-board GPS receiver is considered. The position of the aircraft in geocentric XYZ coordinates $(X_A,Y_A,Z_A)$ is estimated using Kalman filtering, as written below [36,40]:

(a)

process of “prediction”:

$$\{\begin{array}{c}{x}_p^{}=A_0\times x_0^{}\\ P_p^{}=A_0\times P_0^{}\times A_0^T+Q_o^{}\end{array}$$

(2)

where

• $A_0$—the matrix of coefficients,
• $x_0^{}$—the estimated values of the designated parameters a priori from the previous step,
• $P_0^{}$—the estimated values of covariance a priori from the previous step,
• $x_p^{}$—a prediction of state value,
• $P_p^{}$—the predicted covariance values,
• $Q_0^{}$—the variance matrix of the noise of the measurement process.

(b)

process of “correction”:

$$\{\begin{array}{c}{K}_k^{}=P_p^{}\times H_{}^T\times {\left(H\times P_p^{}\times H_{}^T+R\right)}^{-1}\\ x_k^{}=x_p^{}+K_k^{}\times \left(z-H\times x_p^{}\right)\\ P_k^{}=\left(I-K_k^{}\times H\right)\times P_p^{}\end{array}$$

(3)

where

• $K_k^{}$—the Kalman gain matrix,
• $P_k^{}$—the covariance matrix of parameters determined a posteriori,
• $H$—the matrix of partial derivatives,
• $R$—the covariance matrix of measurements,
• $z$—the vector of measured values,
• $I$—the unit matrix,
• $x_k^{}$—the parameters determined a posteriori.

For the purposes of air navigation, the position of an aircraft is much better expressed in ellipsoidal BLh coordinates, as recorded below [41]:

$$\left[\begin{array}{c}{B}_A\\ L_A\\ h_A\end{array}\right]=\left[\begin{array}{c}arc\tan \left(\frac{Z_A}{p}+\frac{{\delta }_1\times tg{B}_{i-1}}{\sqrt{{\delta }_2\times t{g}^2{B}_{i-1}}}\right)\\ arc\tan \left(\frac{Y_A}{X_A}\right)\\ \frac{p}{\cos B}-R\end{array}\right]$$

(4)

where $R$—the radius of the curvature in the prime vertical, $R=\frac{a}{\sqrt{1-e^2\times {\sin }^{2}B}}$, $(a,b)$—the ellipsoid semi-major axis (a) and semi-minor axis (b), $e$—the first flattening, $e=\sqrt{\frac{a^2-b^2}{a^2}}$, $p=\sqrt{X_A^2+Y_A^2}$, ${\delta }_1=\frac{a\times e}{\rho \times \sqrt{1-e^2}}$, ${\delta }_2=\frac{1}{1-e^2}$, $i-1$—the previous iteration, $(B_A,L_A,h_A)$—the aircraft ellipsoidal coordinates, $B_A$—the geodetic latitude, $L_A$—the geodetic longitude, $h_A$—the ellipsoidal height.

Furthermore, for each aircraft position expressed by the BLh components, their mean error values are determined as shown in the scheme below [42]:

$$\left[\begin{array}{c}m{B}_A\\ m{L}_A\\ m{h}_A\end{array}\right]=\left[\begin{array}{c}\sqrt{m_{BLh}(1,1)}\\ \sqrt{m_{BLh}(2,2)}\\ \sqrt{m_{BLh}(3,3)}\end{array}\right]$$

(5)

where

• $m_{BLh}$—the variance–covariance matrix of parameters in the BLh ellipsoidal frame,
• $m_{BLh}=R\times P_k\times R^T$,
• $R$—the conversion matrix from the XYZ geocentric frame to the BLh geodetic frame,
• $m{B}_A$—the mean error for geodetic latitude B,
• $m{L}_A$—the mean error for geodetic longitude L,
• $m{h}_A$—the mean error for ellipsoidal height h.

3.2. New Solution to Improve Positioning Accuracy of DGPS Technique for Multiple Baselines

The assumption in Equation (1), however, has a geometric limitation, because the position of the aircraft is determined only from a single baseline, without controllability, which results in a lack of degrees of freedom and limits the verification of the positioning accuracy determination. A much better solution is to apply DGPS using several baselines. In the presented case, DGPS positioning accuracy was analysed on the basis of a solution using 3 baselines and a weight model. In the DGPS technique, 3 components of the aircraft position are determined; therefore, the minimum geometry of the system consisting of the GNSS reference station network and GPS receiver on-board should also have a minimum of 3. Additionally, such a configuration allows 2 degrees of freedom from 2 independent baselines, which means that the user always has control over the calculation of the determined coordinates of the aircraft position. This approach is very good, because it will always enable the elimination of coarse errors from navigation calculations.

The improvement of the accuracy of the determined position using the DGPS technique for several baselines should be based on a weighting model. In this paper, a weighting model for determining the positioning accuracy for the DGPS technique based on the solution of 3 independent baselines is presented. The weighting of measurements in the case of the differential DGPS technique was determined on the basis of three mathematical models taking into account the length of the baseline, the vector length error and the number of tracked GPS satellites. The scheme for determining the resultant accuracy of the aircraft position using the DGPS technique based on 3 independent baselines can be written as follows:

(a)

criterion I—weighting as a function of baseline (vector) length:

$$\{\begin{array}{c}dB=d{B}_{1A}\times w_{1A}+d{B}_{2A}\times w_{2A}+d{B}_{3A}\times w_{3A}\\ dL=d{L}_{1A}\times w_{1A}+d{L}_{2A}\times w_{2A}+d{L}_{3A}\times w_{3A}\\ dh=d{h}_{1A}\times w_{1A}+d{h}_{2A}\times w_{2A}+d{h}_{3A}\times w_{3A}\end{array}$$

(6)

$$\{\begin{array}{c}dB=d{B}_{1A}\times {\frac{1}{d}}_{1A}+d{B}_{2A}\times {\frac{1}{d}}_{2A}+d{B}_{3A}\times {\frac{1}{d}}_{3A}\\ dL=d{L}_{1A}\times {\frac{1}{d}}_{1A}+d{L}_{2A}\times {\frac{1}{d}}_{2A}+d{L}_{3A}\times {\frac{1}{d}}_{3A}\\ dh=d{h}_{1A}\times {\frac{1}{d}}_{1A}+d{h}_{2A}\times {\frac{1}{d}}_{2A}+d{h}_{3A}\times {\frac{1}{d}}_{3A}\end{array}$$

(7)

where $(dB,dL,dh)$—the resultant position accuracy of the aircraft in ellipsoidal BLh coordinates; $(d{B}_{1A},d{L}_{1A},d{h}_{1A})$—the errors of the aircraft DGPS position in ellipsoidal BLh coordinates determined for the baseline $1A$, relative to the RTK–OTF (Real Time Kinematic–On The Fly) solution [4,43]; $d{B}_{1A}=B_{1A}-Brt{k}_{1A}$; $d{L}_{1A}=L_{1A}-Lrt{k}_{1A}$; $d{h}_{1A}=h_{1A}-hrt{k}_{1A}$; $(B_{1A},L_{1A},h_{1A})$—the position of the aircraft in ellipsoidal BLh coordinates determined from the DGPS solution for the baseline $1A$ (see Equation (4)); $(Brt{k}_{1A},Lrt{k}_{1A},hrt{k}_{1A})$—the reference position of the aircraft in ellipsoidal BLh coordinates derived from the RTK–OTF differential technique for the baseline $1A$; $(d{B}_{2A},d{L}_{2A},d{h}_{2A})$—the errors of the aircraft DGPS position in ellipsoidal BLh coordinates determined for the baseline $2A$, relative to the RTK–OTF solution; $d{B}_{2A}=B_{2A}-Brt{k}_{2A}$; $d{L}_{2A}=L_{2A}-Lrt{k}_{2A}$; $d{h}_{2A}=h_{2A}-hrt{k}_{2A}$; $(B_{2A},L_{2A},h_{2A})$—the position of the aircraft in ellipsoidal BLh coordinates determined from the DGPS solution for the baseline $2A$ (see Equation (4)); $(Brt{k}_{2A},Lrt{k}_{2A},hrt{k}_{2A})$—the reference position of the aircraft in ellipsoidal BLh coordinates derived from the RTK–OTF differential technique for the baseline $2A$; $(d{B}_{3A},d{L}_{3A},d{h}_{3A})$—the errors of the aircraft DGPS position in ellipsoidal BLh coordinates determined for the baseline $3A$, relative to the RTK–OTF solution; $d{B}_{3A}=B_{3A}-Brt{k}_{3A}$; $d{L}_{3A}=L_{3A}-Lrt{k}_{3A}$; $d{h}_{3A}=h_{3A}-hrt{k}_{3A}$; $(B_{3A},L_{3A},h_{3A})$—the position of the aircraft in ellipsoidal BLh coordinates determined from the DGPS baseline solution for the baseline $3A$ (see Equation (4)); $(Brt{k}_{3A},Lrt{k}_{3A},hrt{k}_{3A})$—the reference position of the aircraft in ellipsoidal BLh coordinates derived from the RTK-OTF differential technique for the baseline $3A$; $(w_{1A},w_{2A},w_{3A})$—the weighting functions for individual baselines; $1A$—the baseline between the GNSS reference station ($\overrightarrow{1}$) and the on-board GPS receiver ($\overrightarrow{A}$); $2A$—the baseline between the GNSS reference station ($\overrightarrow{2}$) and the on-board GPS receiver ($\overrightarrow{A}$); $3A$—the baseline between the GNSS reference station ($\overrightarrow{3}$) and the on-board GPS receiver ($\overrightarrow{A}$); $d_{1A}$—the distance between the GNSS reference station ($\overrightarrow{1}$) and the on-board GPS receiver ($\overrightarrow{A}$) expressed in geocentric cartesian coordinates XYZ; $d_{1A}=\sqrt{{(X_1-X_A)}^2+{(Y_1-Y_A)}^2+{(Z_1-Z_A)}^2}$; $d_{2A}$—the distance between the GNSS reference station ($\overrightarrow{2}$) and the on-board GPS receiver ($\overrightarrow{A}$) expressed in geocentric cartesian coordinates XYZ; $d_{2A}=\sqrt{{(X_2-X_A)}^2+{(Y_2-Y_A)}^2+{(Z_2-Z_A)}^2}$; $d_{3A}$—the distance between the GNSS reference station ($\overrightarrow{3}$) and the on-board GPS receiver ($\overrightarrow{A}$) expressed in geocentric cartesian coordinates XYZ; $d_{3A}=\sqrt{{(X_3-X_A)}^2+{(Y_3-Y_A)}^2+{(Z_3-Z_A)}^2}$.

(b)

criterion II—weighting as a function of the error ellipse of the point position:

$$\{\begin{array}{c}dB=d{B}_{1A}\times w_{1A}+d{B}_{2A}\times w_{2A}+d{B}_{3A}\times w_{3A}\\ dL=d{L}_{1A}\times w_{1A}+d{L}_{2A}\times w_{2A}+d{L}_{3A}\times w_{3A}\\ dh=d{h}_{1A}\times w_{1A}+d{h}_{2A}\times w_{2A}+d{h}_{3A}\times w_{3A}\end{array}$$

(8)

$$\{\begin{array}{c}dB=d{B}_{1A}\times \frac{1}{m_{1A}^2}+d{B}_{2A}\times \frac{1}{m_{2A}^2}+d{B}_{3A}\times \frac{1}{m_{3A}^2}\\ dL=d{L}_{1A}\times \frac{1}{m_{1A}^2}+d{L}_{2A}\times \frac{1}{m_{2A}^2}+d{L}_{3A}\times \frac{1}{m_{3A}^2}\\ dh=d{h}_{1A}\times \frac{1}{m_{1A}^2}+d{h}_{2A}\times \frac{1}{m_{2A}^2}+d{h}_{3A}\times \frac{1}{m_{3A}^2}\end{array}$$

(9)

where $m_{1A}$—the vector length error for baseline $1A$, $m_{1A}=\sqrt{{(m_{B,1}^2+m_{B,A}^2)}^2+{(m_{L,1}^2+m_{L,A}^2)}^2+{(m_{h,1}^2+m_{h,A}^2)}^2}$, $(m_{B,1},m_{L,1},m_{h,1})$—the mean errors of the GNSS reference station coordinates ($\overrightarrow{1}$), $(m_{B,1}=0,m_{L,1}=0,m_{h,1}=0)$—the DGPS technique assumes error-free coordinates of the GNSS reference station ($\overrightarrow{1}$), $(m_{B,A},m_{L,A},m_{h,A})$—the mean errors of the determined aircraft coordinates ($\overrightarrow{A}$), $\left[\begin{array}{c}{m}_{B,A}=m{B}_A\\ m_{L,A}=m{L}_A\\ m_{h,A}=m{h}_A\end{array}\right]$ (see Equation (5)), $m_{2A}$—the vector length error for baseline $2A$, $m_{2A}=\sqrt{{(m_{B,2}^2+m_{B,A}^2)}^2+{(m_{L,2}^2+m_{L,A}^2)}^2+{(m_{h,2}^2+m_{h,A}^2)}^2}$, $(m_{B,2},m_{L,2},m_{h,2})$—the mean errors of the GNSS reference station coordinates ($\overrightarrow{2}$), $(m_{B,2}=0,m_{L,2}=0,m_{h,2}=0)$—the DGPS technique assumes error-free coordinates of the GNSS reference station ($\overrightarrow{2}$), $m_{3A}$—the vector length error for baseline $3A$, $m_{3A}=\sqrt{{(m_{B,3}^2+m_{B,A}^2)}^2+{(m_{L,3}^2+m_{L,A}^2)}^2+{(m_{h,3}^2+m_{h,A}^2)}^2}$, $(m_{B,3},m_{L,3},m_{h,3})$—the mean errors of the GNSS reference station coordinates ($\overrightarrow{3}$), $(m_{B,3}=0,m_{L,3}=0,m_{h,3}=0)$—the DGPS technique assumes error-free coordinates of the GNSS reference station ($\overrightarrow{3}$).

(c)

Criterion III—weighting as a function of the number of GPS satellites tracked:

$$\{\begin{array}{c}dB=d{B}_{1A}\times w_{1A}+d{B}_{2A}\times w_{2A}+d{B}_{3A}\times w_{3A}\\ dL=d{L}_{1A}\times w_{1A}+d{L}_{2A}\times w_{2A}+d{L}_{3A}\times w_{3A}\\ dh=d{h}_{1A}\times w_{1A}+d{h}_{2A}\times w_{2A}+d{h}_{3A}\times w_{3A}\end{array}$$

(10)

$$\{\begin{array}{c}dB=d{B}_{1A}\times {\frac{1}{ns}}_{1A}+d{B}_{2A}\times {\frac{1}{ns}}_{2A}+d{B}_{3A}\times {\frac{1}{ns}}_{3A}\\ dL=d{L}_{1A}\times {\frac{1}{ns}}_{1A}+d{L}_{2A}\times {\frac{1}{ns}}_{2A}+d{L}_{3A}\times {\frac{1}{ns}}_{3A}\\ dh=d{h}_{1A}\times {\frac{1}{ns}}_{1A}+d{h}_{2A}\times {\frac{1}{ns}}_{2A}+d{h}_{3A}\times {\frac{1}{ns}}_{3A}\end{array}$$

(11)

where

• $n{s}_{1A}$—the number of tracked GPS satellites for the baseline $1A$,
• $n{s}_{2A}$—the number of tracked GPS satellites for the baseline $2A$,
• $n{s}_{3A}$—the number of tracked GPS satellites for the baseline $3A$.

4. Research Test

The proposed calculation strategy for the mathematical Equations (6)–(11) was implemented and tested for the GPS kinematic data collected from an airborne experiment performed at the civil airport in Mielec, Poland. The test flight was performed over EPML (Northern Europe Poland Mielec) airport in Mielec, Poland using a Seneca Piper PA34-200T (manufacturer: Piper Aircraft, Vero Beach, FL, USA) aircraft. The test flight lasted from 09:55:00 a.m. to 11:10:00 a.m. according to GPST (GPS Time). Figure 1 shows a sketch of the flight trajectory of the Seneca Piper PA34-200T aircraft over the airport in Mielec and the location of the three GNSS base stations.

The aim of the test flight was to determine the basic parameters of the GNSS positioning quality for approach and landing procedures at Mielec airport. As part of the flight test, navigational calculations were also performed for the research method proposed in this paper. Navigational and observation data in RINEX format (Receiver Independent Exchange Format) from the on-board GNSS receiver mounted on Seneca Piper PA34-200T aircraft were used in the calculations. A Topcon HiperPro receiver was placed in the cockpit, which recorded GPS satellite data with an interval of 1 s. Virtual GPS reference stations generated in the positioning service POZGEO-D in ASG-EUPOS (Active Geodetic Network EUPOS) [44] were used as GNSS base stations.

Table 1. The coordinates of GNSS base stations.

Marker of GNSS Base Station	Latitude (B)	Longitude (L)	Ellipsoidal Height (h)
VirA ($\overrightarrow{1}$)	50°25′00″, 00000	20°35′00″, 00000	200,000 m
VirB ($\overrightarrow{2}$)	50°23′00″, 00000	21°43′00″, 00000	200,000 m
VirC ($\overrightarrow{3}$)	50°17′00″, 00000	21°09′00″, 00000	200,000 m

gives the ellipsoidal BLh coordinates of the GNSS base stations. The location of the GNSS reference stations was selected to ensure that the VirA, VirB and VirC stations were evenly distributed along the flight path of the Seneca Piper PA34-200T aircraft.

The distance between the GNSS base stations was as follows:

-

80.603 km between VirA and VirB (designation d1),

-

41.918 km between stations VirA and VirC (designation d2),

-

41.246 km between VirB and VirC (designation d3).

An illustrative sketch of the location of the VirA, VirB and VirC reference stations is shown in Figure 1.

The navigation calculations were divided into two stages. In the first stage, calculations for mathematical Equations (1)–(5) were performed. This stage was implemented in RTKLIB software [45]. RTKPOST application library was used in RTKLIB software, and the following calculation configuration was set up sequentially:

-

positioning mode: DGPS/DGNSS,

-

elevation mask: 5°,

-

type of filtration: forward Kalman filtration,

-

source of ionospheric correction: Klobuchar model,

-

source of tropospheric correction: Saastamoinen model,

-

ephemeris data source: GPS navigation message,

-

GNSS system: GPS system,

-

GPS observation type: C/A code on L1 frequency,

-

resulting coordinates of the aircraft position: ellipsoidal BLh coordinates,

-

base coordinates of the reference stations: catalogue coordinates of VirA, VirB and VirC reference stations according to Table 1,

-

calculation interval: 1 s.

On this basis, three reports were obtained from the RTKPOST calculation mode. The reports contain information on the ellipsoidal coordinates of the plane, the time of measurement, the number of GPS satellites tracked and the average errors of the determined coordinates. In the next stage, all three reports were imported to the Scilab software [46], in which a navigation application for calculations for the scheme of mathematical Equations (6)–(11) was developed. The entire source code with the positioning algorithm for the three weighting criteria was written in Scilab. It should be additionally emphasised that the flight reference position of the Seneca Piper PA34-200T aircraft was also determined in the RTKPOST mode, separately for the three baselines. The results of the calculations are presented in Section 5.

5. Results

The proposed weighted average model presented in Equations (6)–(11) gives the possibility to determine the resultant positioning accuracy of the aircraft for the DGPS technique. In the main, the model is based on the values of position errors $(d{B}_{1A},d{L}_{1A},d{h}_{1A},d{B}_{2A},d{L}_{2A},d{h}_{2A},d{B}_{3A},d{L}_{3A},d{h}_{3A})$ determined for individual baselines. Therefore, in the calculations in Scilab, the parameters $(d{B}_{1A},d{L}_{1A},d{h}_{1A},d{B}_{2A},d{L}_{2A},d{h}_{2A},d{B}_{3A},d{L}_{3A},d{h}_{3A})$ were determined first. The results of the calculations are shown in Figure 2, Figure 3 and Figure 4. The determined position errors of the B coordinate are as follows:

-

from −0.482 to 4.988 m for the baseline $1A$;

-

from −0.750 to 5.286 m for the baseline $2A$;

-

from −0.565 to 3.797 m for the baseline $3A$.

For the L component, the position errors for individual baselines were as follows:

-

from −0.993 to 0.677 m for the baseline $1A$;

-

from −1.160 to 0.950 m for the base line $2A$;

-

from −1.073 to 1.036 m for the baseline $3A$.

In turn, for the ellipsoidal height coordinate, the accuracy values were as follows:

-

from −2.910 to 10.857 m for the baseline $1A$;

-

from −3.296 to 10.269 m for the base line $2A$;

-

from −3.151 to 10.846 m for the base line $3A$.

Based on the position errors obtained, it can be seen that the ellipsoidal height h has the lowest accuracy, while the L component has the highest accuracy.

From the point of view of accuracy analyses, the key parameter is the determination of the geometric coefficient PDOP (Position Dilution of Precision) [47,48]. The determination of the PDOP is of fundamental importance for the determination of position errors, also in the DGPS technique. Figure 5 shows the results of the PDOP parameter for the flight experiment in Mielec. It is worth noting that for the initial measurement epochs, the PDOP values are high, similar to the position errors in Figure 2, Figure 3 and Figure 4. On the other hand, for the rest of measurement epochs, the PDOP values are equal to or less than two. Thus, it can be said that the higher the PDOP coefficient is, the lower the positioning accuracy is, which, in turn, translates into worse measurement conditions and a lower GNSS measurement quality [49]. This phenomenon is especially true for the h-component; it can be seen that the initial values of the position errors in Figure 2, Figure 3 and Figure 4 are up to about 11 m, while the PDOP values are up to about 4.

The algorithm of the DGPS method (see Equation (1)) is characterised by strong correlation between the coded measurements recorded by two GNSS receivers for the same moment of time. This is of course related to the differentiation of code observations and a comparison of the individual code observations from the same satellite for the same epoch. It can be said that differential corrections in the DGPS method stimulate this correlation, which is of course influenced by the measurement process itself and the geometric construction in the DGPS method. This, in turn, ultimately translates into the determination of the aircraft position coordinates and positioning accuracy for a single baseline. This correlation is of course characteristic for the baselines $1A$, $2A$ and $3A$. Thus, comparing with each other the obtained position errors, shown in Figure 2, Figure 3 and Figure 4, will actually show how strong this relationship is. Therefore, the relationship between the individual position errors for each baseline was determined, as recorded below:

$$\{\begin{array}{c}rdB=\{\begin{array}{c}\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{2A}}}\\ \frac{\overline{d{B}_{1A}}}{\overline{d{B}_{3A}}}\\ \frac{\overline{d{B}_{2A}}}{\overline{d{B}_{3A}}}\end{array}\\ rdL=\{\begin{array}{c}\frac{\overline{d{L}_{1A}}}{\overline{d{L}_{2A}}}\\ \frac{\overline{d{L}_{1A}}}{\overline{d{L}_{3A}}}\\ \frac{\overline{d{L}_{2A}}}{\overline{d{L}_{3A}}}\end{array}\\ rdh=\{\begin{array}{c}\frac{\overline{d{h}_{1A}}}{\overline{d{h}_{2A}}}\\ \frac{\overline{d{h}_{1A}}}{\overline{d{h}_{3A}}}\\ \frac{\overline{d{h}_{2A}}}{\overline{d{h}_{3A}}}\end{array}\end{array}$$

(12)

where

$(\overline{d{B}_{1A}},\overline{d{L}_{1A}},\overline{d{h}_{1A}},\overline{d{B}_{2A}},\overline{d{L}_{2A}},\overline{d{h}_{2A}},\overline{d{B}_{3A}},\overline{d{L}_{3A}},\overline{d{h}_{3A}})$—the average accuracy determined for individual baseline,

$(rdB,rdL,rdh)$—the parameters of the relationship between the average accuracy.

Table 2. The results of parameters $(rdB,rdL,rdh)$.

Parameter	Value	Relationship
$rdB$	0.643	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{2A}}}$
$rdB$	0.835	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{3A}}}$
$rdB$	1.326	$\frac{\overline{d{B}_{2A}}}{\overline{d{B}_{3A}}}$
$rdL$	0.856	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{2A}}}$
$rdL$	0.926	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{3A}}}$
$rdL$	1.081	$\frac{\overline{d{B}_{2A}}}{\overline{d{B}_{3A}}}$
$rdh$	0.883	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{2A}}}$
$rdh$	0.923	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{3A}}}$
$rdh$	1.046	$\frac{\overline{d{B}_{2A}}}{\overline{d{B}_{3A}}}$

shows the results of determining the parameters of the relationship $(rdB,rdL,rdh)$. On this basis, it can be seen that the relationship between the position errors between baseline $1A$ and $2A$, as well as $1A$ and $3A$ is smaller than one. On the other hand, the relationship between baseline $2A$ and $3A$ is larger than one. Going back to Figure 1 and the distribution of the GNSS reference stations, it can be seen that the distance between the VirB and VirC stations is the shortest. It can be concluded that the smaller the distance is between stations, the greater the relationship is between position errors.

For the developed calculation algorithm (Equations (6)–(11)), the key element is the determination of measurement weights for criteria I–III. The measurement weights were defined as a function of the baseline (vector) length, the vector length error and the number of GPS satellites tracked. The multiplicity of weighting functions is good, as it allows the user to determine the optimal measurement weighting strategy, which, in turn, translates into an improvement of the determined aircraft coordinates. Figure 6, Figure 7 and Figure 8 show the measurement weighting values for criteria I–III. In the case of weighting in scheme I, the values of measurement weights are from 0.013 to 0.239 for baseline $1A$, from 0.013 to 0.169 for baseline $2A$ and from 0.025 to 0.164 for baseline $3A$. It is worth noting a simple relationship that follows directly from Figure 6. Namely, the more the measurement weight increases, the shorter the geometric distance between the on-board receiver and the GNSS reference station will be. This relationship of course occurs in reverse and is alternating.

Figure 7 shows the results of the measurement weights for criterion II. For this criterion, the measurement weights range from 0.003 to 0.125 for baseline $1A$, from 0.003 to 0.100 for baseline $2A$ and from 0.010 to 0.096 for baseline $3A$. For weighting criterion II, the measurement weights vary inversely proportional to the squares of the mean errors, according to Equation (9). It can be seen that the highest value of measurement weights for criterion II reaches 0.125.

Figure 8 presents the results of the measurement weights for criterion III. In this criterion, the measurement weights range from 0.111 to 0.250 for all the baselines, $1A$, $2A$ and $3A$. For weighting criterion III, the measurement weights are, in most cases, equal and the same.

The following Figure 9, Figure 10 and Figure 11 show the results of determining the resultant values of the position errors determined according to Equations (7), (9) and (11). For the B coordinate, the values of the resultant position errors are as follows:

-

from −0.508 to 0.519 m for the weighting in criterion I,

-

from −0.259 to 0.265 m for the weighting in criterion II,

-

from −1.346 to 3.518 m for the weighting in criterion III.

Looking at Figure 9, it can be seen quite well how the measurement weights have reduced the position error values along the B axis. The smallest spread of position errors is $dB$ seen for weighting criterion II and the largest for weighting criterion III.

For the L component, the resultant position error values $dL$ are as follows:

-

from −0.488 to 0.091 m for the weighting in criterion I,

-

from −0.269 to 0.051 m for the weighting in criterion II,

-

from −1.122 to 0.223 m for the weighting in criterion III.

Comparing the results in Figure 9 and Figure 10, it can be observed that the position errors $dL$ show a smaller scatter than the parameters $dB$. This is quite an important message, saying that the positioning accuracy along the L-axis is the highest.

For the height coordinate h, the resultant position error values $dh$ are as follows:

-

from −0.078 to 1.334 m for the weighting in criterion I,

-

from −0.039 to 0.696 m for the weighting in criterion II,

-

from −0.332 to 6.764 m for the weighting in criterion III.

Analysing the results in Figure 9, Figure 10 and Figure 11, it can be said that the positioning accuracy is highest for the L component and lowest for the h coordinate. This only confirms that the dispersion of position errors is the highest for the h component, while the lowest for the L coordinate.

In the discussion of the results of the resultant position errors, a statement was made that the weighting scheme for criteria I–III significantly reduced the parameter values $(d{B}_{1A},d{L}_{1A},d{h}_{1A},d{B}_{2A},d{L}_{2A},d{h}_{2A},d{B}_{3A},d{L}_{3A},d{h}_{3A})$. It is, therefore, necessary to determine the percentage reductions for these position errors. For this purpose, the reductions in the resultant position errors $(dB,dL,dh)$ relative to the parameters $(d{B}_{1A},d{L}_{1A},d{h}_{1A},d{B}_{2A},d{L}_{2A},d{h}_{2A},d{B}_{3A},d{L}_{3A},d{h}_{3A})$, as shown in Equations (13)–(15), were calculated as follows:

$$\{\begin{array}{c}P{B}_I[\%]=\{\begin{array}{c}\frac{(\overline{d{B}_I}-\overline{d{B}_{1A}})}{\overline{d{B}_{1A}}}\times 100\%\\ \frac{(\overline{d{B}_I}-\overline{d{B}_{2A}})}{\overline{d{B}_{2A}}}\times 100\%\\ \frac{(\overline{d{B}_I}-\overline{d{B}_{3A}})}{\overline{d{B}_{3A}}}\times 100\%\end{array}\\ P{B}_{II}[\%]=\{\begin{array}{c}\frac{(\overline{d{B}_{II}}-\overline{d{B}_{1A}})}{\overline{d{B}_{1A}}}\times 100\%\\ \frac{(\overline{d{B}_{II}}-\overline{d{B}_{2A}})}{\overline{d{B}_{2A}}}\times 100\%\\ \frac{(\overline{d{B}_{II}}-\overline{d{B}_{3A}})}{\overline{d{B}_{3A}}}\times 100\%\end{array}\\ P{B}_{III}[\%]=\{\begin{array}{c}\frac{(\overline{d{B}_{III}}-\overline{d{B}_{1A}})}{\overline{d{B}_{1A}}}\times 100\%\\ \frac{(\overline{d{B}_{III}}-\overline{d{B}_{2A}})}{\overline{d{B}_{2A}}}\times 100\%\\ \frac{(\overline{d{B}_{III}}-\overline{d{B}_{3A}})}{\overline{d{B}_{3A}}}\times 100\%\end{array}\end{array}$$

(13)

$$\{\begin{array}{c}P{L}_I[\%]=\{\begin{array}{c}\frac{(\overline{d{L}_I}-\overline{d{L}_{1A}})}{\overline{d{L}_{1A}}}\times 100\%\\ \frac{(\overline{d{L}_I}-\overline{d{L}_{2A}})}{\overline{d{L}_{2A}}}\times 100\%\\ \frac{(\overline{d{L}_I}-\overline{d{L}_{3A}})}{\overline{d{L}_{3A}}}\times 100\%\end{array}\\ P{L}_{II}[\%]=\{\begin{array}{c}\frac{(\overline{d{L}_{II}}-\overline{d{L}_{1A}})}{\overline{d{L}_{1A}}}\times 100\%\\ \frac{(\overline{d{L}_{II}}-\overline{d{L}_{2A}})}{\overline{d{L}_{2A}}}\times 100\%\\ \frac{(\overline{d{L}_{II}}-\overline{d{L}_{3A}})}{\overline{d{L}_{3A}}}\times 100\%\end{array}\\ P{L}_{III}[\%]=\{\begin{array}{c}\frac{(\overline{d{L}_{III}}-\overline{d{L}_{1A}})}{\overline{d{L}_{1A}}}\times 100\%\\ \frac{(\overline{d{L}_{III}}-\overline{d{L}_{2A}})}{\overline{d{L}_{2A}}}\times 100\%\\ \frac{(\overline{d{L}_{III}}-\overline{d{L}_{3A}})}{\overline{d{L}_{3A}}}\times 100\%\end{array}\end{array}$$

(14)

$$\{\begin{array}{c}P{h}_I[\%]=\{\begin{array}{c}\frac{(\overline{d{h}_I}-\overline{d{h}_{1A}})}{\overline{d{h}_{1A}}}\times 100\%\\ \frac{(\overline{d{h}_I}-\overline{d{h}_{2A}})}{\overline{d{h}_{2A}}}\times 100\%\\ \frac{(\overline{d{h}_I}-\overline{d{h}_{3A}})}{\overline{d{h}_{3A}}}\times 100\%\end{array}\\ P{h}_{II}[\%]=\{\begin{array}{c}\frac{(\overline{d{h}_{II}}-\overline{d{h}_{1A}})}{\overline{d{h}_{1A}}}\times 100\%\\ \frac{(\overline{d{h}_{II}}-\overline{d{h}_{2A}})}{\overline{d{h}_{2A}}}\times 100\%\\ \frac{(\overline{d{h}_{II}}-\overline{d{h}_{3A}})}{\overline{d{h}_{3A}}}\times 100\%\end{array}\\ P{h}_{III}[\%]=\{\begin{array}{c}\frac{(\overline{d{h}_{III}}-\overline{d{h}_{1A}})}{\overline{d{h}_{1A}}}\times 100\%\\ \frac{(\overline{d{h}_{III}}-\overline{d{h}_{2A}})}{\overline{d{h}_{2A}}}\times 100\%\\ \frac{(\overline{d{h}_{III}}-\overline{d{h}_{3A}})}{\overline{d{h}_{3A}}}\times 100\%\end{array}\end{array}$$

(15)

where $(P{B}_I,P{B}_{II},P{B}_{III})$—the improvement of the positioning accuracy for the B-component for criteria I–III presented as a percentage, $(P{L}_I,P{L}_{II},P{L}_{III})$—the improvement of the positioning accuracy for the L-component for criteria I–III presented as a percentage, $(P{h}_I,P{h}_{II},P{h}_{III})$—the improvement of the positioning accuracy for the h-component for criteria I–III presented as a percentage, $(\overline{d{B}_{1A}},\overline{d{L}_{1A}},\overline{d{h}_{1A}},\overline{d{B}_{2A}},\overline{d{L}_{2A}},\overline{d{h}_{2A}},\overline{d{B}_{3A}},\overline{d{L}_{3A}},\overline{d{h}_{3A}})$—the average accuracy specified for individual baselines, $(\overline{d{B}_I},\overline{d{L}_I},\overline{d{h}_I})$—the mean positioning accuracy specified for weighting criterion I, $(\overline{d{B}_{II}},\overline{d{L}_{II}},\overline{d{h}_{II}})$—the mean positioning accuracy specified for weighting criterion II, $(\overline{d{B}_{III}},\overline{d{L}_{III}},\overline{d{h}_{III}})$—the mean positioning accuracy specified for weighting criterion III.

In the analysed example, the results of the parameters are presented as $(P{B}_I,P{B}_{II},P{B}_{III},P{L}_I,P{L}_{II},P{L}_{III},P{h}_I,P{h}_{II},P{h}_{III})$ in

Table 3. The results of reducing the position errors according to Equations (13)–(15).

Parameter	Percentage Value (%)	Criterion of Measurement Weights	Baseline
$P{B}_I$	83	I	$1A$
$P{B}_I$	89	I	$2A$
$P{B}_I$	86	I	$3A$
$P{B}_{II}$	91	II	$1A$
$P{B}_{II}$	94	II	$2A$
$P{B}_{II}$	92	II	$3A$
$P{B}_{III}$	55	III	$1A$
$P{B}_{III}$	71	III	$2A$
$P{B}_{III}$	62	III	$3A$
$P{L}_I$	86	I	$1A$
$P{L}_I$	88	I	$2A$
$P{L}_I$	87	I	$3A$
$P{L}_{II}$	92	II	$1A$
$P{L}_{II}$	94	II	$2A$
$P{L}_{II}$	93	II	$3A$
$P{L}_{III}$	62	III	$1A$
$P{L}_{III}$	67	III	$2A$
$P{L}_{III}$	64	III	$3A$
$P{h}_I$	86	I	$1A$
$P{h}_I$	87	I	$2A$
$P{h}_I$	87	I	$3A$
$P{h}_{II}$	92	II	$1A$
$P{h}_{II}$	93	II	$2A$
$P{h}_{II}$	93	II	$3A$
$P{h}_{III}$	63	III	$1A$
$P{h}_{III}$	67	III	$2A$
$P{h}_{III}$	65	III	$3A$

. On the basis of the results presented in Table 3, it can be concluded that the reduction in position errors using Equations (6)–(11) is quite significant. It is only worth emphasising that in the case of the B component, one can speak of a reduction at the level of 55–94%, depending on the weighting criterion and the baseline. For the L component, the improvement in positioning accuracy ranges from 62% to 94%. For the ellipsoidal height h, on the other hand, the improvement in positioning accuracy is between 63 and 93%. As for the individual weighting criteria, criterion II gives the best results in terms of improvement in positioning accuracy. In contrast, the positioning accuracy improvement results are the worst for weighting criterion III compared to the values for criteria I and II. The results in Table 3 show how important the measurement weighting scheme is for the differential DGPS technique.

In the analysis of the results, the parameter of the arithmetic mean for the resultant position errors from the individual criteria I–III was also determined, as recorded below:

$$\{\begin{array}{c}dBw=\frac{d{B}_I+d{B}_{II}+d{B}_{III}}{3}\\ dLw=\frac{d{L}_I+d{L}_{II}+d{L}_{III}}{3}\\ dhw=\frac{d{h}_I+d{h}_{II}+d{h}_{III}}{3}\end{array}$$

(16)

where

• $(dBw,dLw,dhw)$—the arithmetic mean of the resultant position errors for criteria I–III,
• $(d{B}_I,d{B}_{II},d{B}_{III})$—the resultant position errors for coordinate B for criteria I–III,
• $(d{L}_I,d{L}_{II},d{L}_{III})$—the resultant position errors for coordinate L for criteria I–III,
• $(d{h}_I,d{h}_{II},d{h}_{III})$—the resultant position errors for coordinate h for criteria I–III.

Figure 12 presents the results of the parameter calculations $(dBw,dLw,dhw)$. The results of the parameters $(dBw,dLw,dhw)$ produce the following values:

-

from −0.114 to −0.654 m for the parameter $dBw$,

-

from −0.197 to −0.539 m for the parameter $dLw$,

-

from 0.575 to −0.146 m for the parameter $dhw$.

For the analysed parameters $(dBw,dLw,dhw)$, the highest positioning accuracy is obviously for the L component, while the lowest is for the h coordinate.

6. Discussion

The scheme of the proposed mathematical model (Equations (6)–(16)) was verified and checked for the accuracy improvement of the DGPS method in an additional airborne experiment performed at the EPDE (European Poland Dęblin) airport in Dęblin. The flight was performed with a Cessna 172 aircraft and the duration of the experiment was from 09:48:00 a.m. to 11:05:44 a.m. according to GPST time. Additionally, in this case, three reference stations named VirA, VirB and VirC were placed along the flight route. The coordinates of the GNSS reference stations were as follows:

-

B = 51.816666667°, L = 21.866666666°, h = 150,000 m for VirA station;

-

B = 51.616666667°, L = 21.516666666°, h = 150,000 m for VirB station;

-

B = 51.46666666666°, L = 22.06666666667°, h = 150,000 m for VirC station.

Figure 13 presents the horizontal trajectory of the Cessna 172 aircraft with the localisation of the GNSS base stations.

RTKPOST software was also used in the calculations to obtain results from the DGPS and RTK–OTF techniques. In addition, the algorithm scheme (Equations (6)–(16)) was made in Scilab software and the calculation results are shown in a graphical and tabular form.

Figure 14, Figure 15 and Figure 16 show the values of the position errors $(d{B}_{1A},d{L}_{1A},d{h}_{1A},d{B}_{2A},d{L}_{2A},d{h}_{2A},d{B}_{3A},d{L}_{3A},d{h}_{3A})$, determined as the difference in the coordinates of the Cessna 172 (manufacturer: Cessna Aircraft Company, USA) aircraft between the DGPS and RTK–OTF solutions for the individual baseline. The determined position errors of the BLh coordinate between the DGPS and RTK–OTF solution are as follows:

For the B component:

-

from −6.302 to 6.075 m for the baseline $1A$;

-

from −6.262 to 6.170 m for the baseline $2A$;

-

from −6.276 to 6.253 m for the baseline $3A$.

For the L component:

-

from −4.714 to 0.547 m for the baseline $1A$;

-

from −4.746 to 0.542 m for the baseline $2A$;

-

from −4.749 to 0.538 m for the baseline $3A$.

For the h component:

-

from −12.987 to 7.130 m for the baseline $1A$;

-

from −13.175 to 7.208 m for the baseline $2A$;

-

from −13.352 to 6.941 m for the baseline $3A$.

Similarly, to the results from the flight experiment in Mielec, the dispersion of position errors for the L component was the smallest for the flight test carried out in Dęblin. On the other hand, the largest dispersion of position errors occurs for the h coordinate.

Figure 17 shows the results of the PDOP parameter for the flight experiment in Dęblin. It is worth emphasising that for the initial measurement epochs, the PDOP values are very high, since they amount to about 10. However, for the remaining measurement epochs, PDOP values are on average about 2.5, with the exception of single outliers reaching the level of 3.5. It should be added that the PDOP values from individual baselines are practically identical in the flight experiment in Dęblin; therefore, their correlation between each other is very close to one. It is important from the point of view of satellite-receiver system geometry for a single baseline. Comparing the PDOP results obtained in Mielec and Dęblin experiments, similar relations are visible. That is, large PDOP values are seen in the initial phase of the flight. On the other hand, for the remaining measurement epochs, the measurement conditions are favourable, with a PDOP value of about 2–2.5.

Table 4. The results of parameters $(rdB,rdL,rdh)$ in flight test in Dęblin.

Parameter	Value	Relation
$rdB$	1.029	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{2A}}}$
$rdB$	1.004	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{3A}}}$
$rdB$	0.975	$\frac{\overline{d{B}_{2A}}}{\overline{d{B}_{3A}}}$
$rdL$	1.001	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{2A}}}$
$rdL$	1.029	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{3A}}}$
$rdL$	1.029	$\frac{\overline{d{B}_{2A}}}{\overline{d{B}_{3A}}}$
$rdh$	0.976	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{2A}}}$
$rdh$	0.957	$\frac{\overline{d{B}_{1A}}}{\overline{d{B}_{3A}}}$
$rdh$	0.989	$\frac{\overline{d{B}_{2A}}}{\overline{d{B}_{3A}}}$

shows the results of determining the parameters of the relationship $(rdB,rdL,rdh)$. In this case, one can notice a regularity that basically all the parameter results are close to one. While in the flight test in Mielec, the parameter results $(rdB,rdL,rdh)$ were varied, in the flight experiment in Dęblin, this is no longer the case. Interestingly, the values of the $rdL$ correlations are always larger than one, on the other hand, the parameters $rdh$ are always smaller than one. This simply shows that the nature of the correlations of the $(rdB,rdL,rdh)$ parameters in the flight experiment in Dęblin is higher due to equal distribution of GNSS base stations and shorter distances between them than in the Mielec experiment.

Figure 18, Figure 19 and Figure 20 show the values of measurement weights for the flight experiment carried out in Dęblin. The calculations resulted in the following values of measurement weights:

For weighting criterion I:

-

from 0.025 to 0.624 for the weight $w_{1A}$,

-

from 0.025 to 0.245 for the weight $w_{2A}$,

-

from 0.023 to 0.628 for the weight $w_{3A}$.

For weighting criterion II:

-

from 0.010 to 0.186 for the weight $w_{1A}$,

-

from 0.010 to 0.127 for the weight $w_{2A}$,

-

from 0.009 to 0.186 for the weight $w_{3A}$.

For weighting criterion III:

-

from 0.111 to 0.250 for the weight $w_{1A}$, $w_{2A}$, $w_{3A}$.

It can be observed that the measurement weights determined as a function of the vector length have quite large values. In turn, the measurement weights as a function of the vector length error and the number of tracked GPS satellites are at a similar level, as in the Mielec flight experiment.

Figure 21, Figure 22 and Figure 23 present the results of the resultant position errors for the experiment in Dęblin, calculated according to Equations (7), (9) and (11). The resultant position errors for the flight experiment in Dęblin were as follows:

For the B coordinate:

-

from −2.268 to 0.824 m for weighting criterion I,

-

from −0.901 to 0.440 m for weighting criterion II,

-

from −2.355 to 3.083 m for weighting criterion III.

For the L coordinate:

-

from −1.556 to 0.072 m for weighting criterion I,

-

from −0.550 to 0.038 m for weighting criterion II,

-

from −2.842 to 0.271 m for weighting criterion III.

For the h coordinate:

-

from −3.388 to 0.942 m for weighting criterion I,

-

from −1.137 to 0.502 m for weighting criterion II,

-

from −6.586 to 2.950 m for weighting criterion III.

The results of the resultant position errors show that the smallest scatter occurs for the L component and the largest occurs for the h component. Moreover, the highest positioning accuracy is observed for the weighting criterion II. Therefore, the results obtained from the flight test in Dęblin only confirm the conclusions drawn for the test results determined in the flight experiment in Mielec (see the descriptions to Figure 9, Figure 10 and Figure 11).

Table 5. The results of reducing the position errors in flight test in Dęblin.

Parameter	Percentage Value (%)	Criterion of Measurement Weights	Baseline
$P{B}_I$	82	I	$1A$
$P{B}_I$	82	I	$2A$
$P{B}_I$	82	I	$3A$
$P{B}_{II}$	91	II	$1A$
$P{B}_{II}$	91	II	$2A$
$P{B}_{II}$	91	II	$3A$
$P{B}_{III}$	65	III	$1A$
$P{B}_{III}$	64	III	$2A$
$P{B}_{III}$	65	III	$3A$
$P{L}_I$	83	I	$1A$
$P{L}_I$	83	I	$2A$
$P{L}_I$	82	I	$3A$
$P{L}_{II}$	91	II	$1A$
$P{L}_{II}$	91	II	$2A$
$P{L}_{II}$	91	II	$3A$
$P{L}_{III}$	65	III	$1A$
$P{L}_{III}$	64	III	$2A$
$P{L}_{III}$	63	III	$3A$
$P{h}_I$	78	I	$1A$
$P{h}_I$	81	I	$2A$
$P{h}_I$	81	I	$3A$
$P{h}_{II}$	90	II	$1A$
$P{h}_{II}$	90	II	$2A$
$P{h}_{II}$	90	II	$3A$
$P{h}_{III}$	63	III	$1A$
$P{h}_{III}$	64	III	$2A$
$P{h}_{III}$	64	III	$3A$

presents the parameter results $(P{B}_I,P{B}_{II},P{B}_{III},P{L}_I,P{L}_{II},P{L}_{III},P{h}_I,P{h}_{II},P{h}_{III})$ for the Dęblin flight experiment. The results obtained only confirm the validity of the proposed research method. Moreover, the values of parameters $(P{B}_I,P{B}_{II},P{B}_{III},P{L}_I,P{L}_{II},P{L}_{III},P{h}_I,P{h}_{II},P{h}_{III})$ are at a similar level, as in the Mielec flight experiment. Based on Table 5, it can be seen that the positioning accuracy for weighting criterion II is the highest. On the other hand, for weighting criterion III, the worst results were obtained in terms of the improvement of positioning accuracy in relation to criteria I and II. For weighting criterion I, the improvement in accuracy ranges from 78% to 83%. For criterion II, the improvement in accuracy is 90–91%, while for criterion III, the improvement in DGPS positioning accuracy is 63–65%.

The final stage in the development of the test results for the flight experiment in Dęblin is the determination of the parameters $(dBw,dLw,dhw)$, as shown in Figure 24. The values of the parameters $(dBw,dLw,dhw)$ are as follows:

-

from −1.679 to 1.449 m for the parameter $dBw$,

-

from −1.270 to 0.127 m for the parameter $dLw$,

-

from −3.095 to 1.368 m for the parameter $dhw$.

As in the Mielec flight experiment and also in Dęblin experiment, the highest positioning accuracy is seen for the L component, while the lowest positioning accuracy is for the h coordinate.

The research problem shown in this paper concerned the accuracy of DGPS positioning in air navigation. The subject of accuracy research was analysed in scientific papers [6,13,32,33,34,39]. In these works, the accuracy of DGPS or DGLONASS positioning was also determined in relation to the differential RTK technique. The developed test methodology from Equations (6)–(16) can significantly improve the accuracy of DGPS positioning in air navigation, especially for a solution using several GNSS reference stations. Moreover, the developed algorithm shows that it is still possible to improve the accuracy of DGPS positioning in air navigation, especially when using different configuration of weighting functions. The submitted work is part of a research trend of studies on improving the accuracy of DGPS positioning, similar to articles [6,13,32,33,34,39].

A similar scope of research, as in the presented article, was conducted in works [50,51]. However, in these papers the research was carried out to improve the accuracy of DGPS positioning in static measurements for a network of GNSS referee stations. In contrast, our work shows an improvement of positioning accuracy in the kinematic mode, during flight experiments. Moreover, papers [50,51] present a mathematical model for the improvement of positioning accuracy using the weighting of measurements as a function of elevation angle to the satellite. On the other hand, our paper presents three weighting schemes, independent from each other, thanks to which it was possible to reduce the positioning errors. The approach presented in this paper and in papers [50,51] shows and emphasises how much of an important role in DGPS positioning model and weighting functions play, enabling us to improve the accuracy parameter.

7. Conclusions

This paper proposes a new research methodology for improving the accuracy of DGPS positioning using several GNSS reference stations. The paper proposes the use of a new weighting scheme for improving the determination of position errors in a differential DGPS technique. Namely, weighting functions in the form of the inverse of the vector length, the inverse of the square of the distance error and the inverse of the number of tracked GPS satellites for a single baseline were used in the calculations. On this basis, the equations of the linear combinations linking the weighting functions and the position errors determined for individual baselines were written. The weighting scheme for three independent baselines was used in the calculations. The proposed algorithm makes it possible to determine the resultant position errors for ellipsoidal BLh coordinates of the aircraft and significantly improve the accuracy of DGPS positioning.

The study used real observations and GPS navigation data from a Topcon HiperPro on-board receiver and the following three GNSS reference stations: VirA, VirB and VirC. Calculations were performed for the flight experiment in Mielec for the test flight carried out using a Seneca Piper PA34-200T aircraft. On the basis of the performed calculations, the following was found:

-

in the case of the B component, there is a 55–94% error reduction;

-

for the L component, the improvement in positioning accuracy ranges from 62% to 94%;

-

for the ellipsoidal height h, the improvement in positioning accuracy ranges from 63% to 93%.

In addition, it is shown that weighting criterion II enables an improvement in DGPS positioning accuracy of 91–94%. In contrast, for weighting criterion I, the accuracy improvement is in the range of 83–89%, while for weighting criterion III it is 55–71%.

As part of the research, DGPS data were also processed for an airborne experiment carried out with a Cessna 172 aircraft at Dęblin airport. In the case of this flight experiment, the proposed research method allowed for an accuracy improvement from 78 to 83% for weighting criterion I, 90–91% for weighting criterion II and 63–65% for weighting criterion III.

The research methodology shown in this paper proved to be effective for improving the accuracy of DGPS positioning. Therefore, the assumptions and purpose of the work were fully achieved. Thus, the obtained research results are satisfactory for the performance of a differential DGPS technique in air navigation. Moreover, two independent flight experiments carried out in Mielec and Dęblin confirmed the high efficiency of the proposed research method.

The mathematical model described in this paper was aimed at improving positioning accuracy in the kinematic mode for air navigation. The presented model has an application in positioning with GNSS satellite technology; hence, the presented work is focused on this problem. The research is based on the fact that in air navigation newer and newer navigation solutions are sought to improve the quality of positioning, including the accuracy parameter. The presented navigation solutions only emphasise that the model of weighting coefficients can be an interesting research method in the area of DGPS positioning in air navigation. The model of weighting coefficients may turn out to be an effective solution for improving the accuracy of DGPS positioning and, thus, for improving the quality of the determined coordinates of an aircraft.