New Methodology for Computing the Aircraft’s Position Based on the PPP Method in GPS and GLONASS Systems

: In the ﬁeld of air navigation, there is a constant pursuit for new navigation solutions for precise GNSS (Global Navigation Satellite System) positioning of aircraft. This study aims to present the results of research on the development of a new method for improving the performance of PPP (Precise Point Positioning) positioning in the GPS (Global Positioning System) and GLONASS (Globalnaja Nawigacionnaja Sputnikovaya Sistema) systems for air navigation. The research method is based on a linear combination of individual position solutions from the GPS and GLONASS systems. The paper shows a computational scheme based on the linear combination for geocentric XYZ coordinates of an aircraft. The algorithm of the new research method uses the weighted mean method to determine the resultant aircraft position. The research method was tested on GPS and GLONASS kinematic data from an airborne experiment carried out with a Seneca Piper PA34-200T aircraft at the Mielec airport. A dual-frequency dual-system GPS/GLONASS receiver was placed on-board the plane, which made it possible to record GNSS observations, which were then used to calculate the aircraft’s position in CSRS-PPP software. The calculated XYZ position coordinates from the CSRS-PPP software were then used in the weighted mean model’s developed optimization algorithm. The measurement weights are a function of the number of GPS and GLONASS satellites and the inverse of the mean error square. The obtained coordinates of aircraft from the research model were veriﬁed with the RTK-OTF solution. As a result of the research, the presented solution’s accuracy is better by 11–87% for the model with a weighting scheme as a function of the inverse of the mean error square. Moreover, using the XYZ position from the RTKLIB program, the research method’s accuracy increases from 45% to 82% for the model with a weighting scheme as a function of the inverse of the square of mean error. The developed method demonstrates high efﬁciency for improving the performance of GPS and GLONASS solutions for the PPP measurement technology in air navigation.


Introduction
To date, the trend of using GNSS (Global Navigation Satellite System) satellite technology in aviation has been exclusively based on single-frequency receivers, recording the publicly available L1-C/A code measurement. Consequently, ICAO (International Civil Aviation Organization) technical standards in the aspect of quality of GNSS positioning in aviation are referred to the L1 frequency [1]. Thus, for many years, only GNSS L1-C/A code observations were used in flight tests to determine the accuracy, reliability, continuity and availability of the aircraft navigation solution. However, the development of GNSS satellite receiver technology and the growth of the GNSS space segment made it possible to use observations on the other carrier frequencies as well. In practice, the use of GNSS dual-frequency observations in navigation calculations adds value to the aircraft position determination process.
For GNSS satellite observations, we distinguish both absolute and differential (relative) positioning methods in navigation calculations [2]. Among the absolute positioning methods, one of the precise measurement techniques is the dual-frequency PPP (Precise Point Positioning) method [3]. The PPP measurement technique uses code observations (P1 and P2) and phase observations (L1 and L2) from a dual-frequency GNSS satellite receiver on-board the aircraft [3]. Thus, during the flight operation, the aircraft crew should be equipped only with a dual-frequency GNSS receiver for precise computation of the aircraft flight path using the PPP method. A completely different situation occurs in the case of differential positioning, where the differential RTK (Real-Time Kinematic) technique with an OTF (On The Fly) solution for GNSS phase measurements is used as a precise measurement technique [4]. The differential RTK technique uses a GNSS satellite receiver on-board the aircraft and a GNSS satellite receiver installed at the reference station (base station) [5]. The differential RTK technique requires the use of a minimum of two GNSS satellite receivers, which, in turn, translates into increased costs and the need to purchase technical infrastructure. In order to reduce the financial costs and the number of personnel to operate the GNSS reference station infrastructure, the PPP measurement technique is increasingly used in flight tests.
The common feature of the absolute PPP method and the differential RTK technique is their use to determine an aircraft's reference position in flight tests. Consequently, both measurement methods are significant for the problem of determining the true trajectory of an aircraft determined from dual-frequency GNSS observations. Knowing the aircraft's true position is essential for both pilots and air traffic ground personnel, including controllers and navigators.

Related Papers
The problem of determining an aircraft's actual position using the PPP method has been presented and described in many research papers worldwide. Most research papers on PPP positioning in air navigation had described and presented position solutions using GPS (Global Positioning System) or GLONASS (Globalnaja Nawigacionnaja Sputnikovaya Sistema) navigation systems. That is, of course, related to the fact that only GPS and GLONASS systems are certified by ICAO for GNSS satellite navigation in aviation [1]. It can be concluded that GPS and GLONASS, as global GNSS navigation systems, provide continuous satellite positioning in aviation.
The literature review provides a multi-faceted approach to the application of the PPP measurement technique in air navigation. Reference [6] propose using a PPP solution to determine the position of an aircraft with respect to the differential solution from the AUS-POS program. The aircraft's position was determined in the CSRS-PPP (Canadian Spatial Reference System-Precise Point Positioning) program using the least-squares method for GPS observations. The maximum coordinate difference between the GAPS (GPS Analysis and Positioning Software) and AUSPOS (online GPS processing service) solution was up to ±0.5 m for the horizontal coordinates and ±1 m for the vertical coordinate. In turn, the paper [7] proposed the use of a PPP solution to determine the position of a helicopter with respect to a DD (Double Difference) differential solution in postprocessing. The helicopter position was determined in the P3 (software developed by the University of Calgary) software package using Kalman filtering for GPS observations. The maximum coordinate difference between the PPP solution and the differential DD solution reached ±0.5 m for all ENU (East-North-Up) components. A similar solution was proposed in work [8], in which helicopter positioning accuracy was determined using the P3 (software developed by the University of Calgary) software package with the PPP method. The initial positioning accuracy of the helicopter was about 2.7 m, while, in the final phase of the flight, it was higher than 0.1 m for BLh (B: Latitude, L: Longitude and h: ellipsoidal height) coordinates. In paper [9], P3 software was used to determine the positioning accuracy of an aircraft flying at high speed. The values of positioning errors of the aircraft between the PPP solution and the differential DD solution were up to ±0.2 m for the horizontal coordinates and ±1.1 m for the vertical coordinate. However, paper [10] used the PPP method to determine a helicopter's position with respect to the KGPS (Kinematic GPS) solution in the PNAV program. The helicopter position was determined using Kalman filtering for GPS observations. The maximum coordinate difference between the PPP solution and the KGPS PNAV solution was up to ±1 m for all ENU components. Furthermore, in work [11], the PPP method was used to determine the UAV (Unmanned Aerial Vehicle) position in GIPSY and RTKLIB software. The UAV reference position was determined using the DGPS (Differential GPS) technique for GPS phase observations. The paper presents different scenarios for UAV position determination for the PPP method: iterative method, filtering of GPS observations and the smoothing of code observations by GPS phase observations. The maximum coordinate difference between the PPP and DGPS solutions reached ±1 m for all ENU components. In paper [12], the PPP method's accuracy in PP-RTX software was investigated in more than 100 flight tests. The PP-RTX solution's typical accuracy with respect to the ASB (Applanix Smartbase) differential solution was up to ±0.2 m for all ENU components. In turn, the paper [13] proposed to test the PPP method in TriP, IT and GrafNav software in a flight experiment. This paper presented the results of the comparison of the aircraft coordinate determination between the TriP, IT, GrafNav software and the GPSurvey baseline solution. The typical positioning accuracy of the aircraft from the PPP method was up to ±0.7 m for all ENU components. On the other hand, the paper [14] presented the results of the positioning accuracy of the aircraft from the CSRS-PPP and GrafNav solutions concerning the DD differential solution. The average positioning accuracy from the PPP method was ±0.2 ÷ 0.4 m; however, numerous anomalies can be observed from the GrafNav solution when the accuracy decreases to the level of ±1.5 m. Interesting solutions for applying the PPP method for single-frequency GPS receivers were presented in papers [15,16]. The papers described the application of the method of smoothing code measurements with GPS phase measurements using IGS (International GNSS Service) products in air navigation. The typical GPS positioning accuracy for the PPP method in works [15,16] was about ±1 m for all ENU and BLh frame components. Another very interesting solution for the needs of aerial photogrammetry was the application of the PPP method in the process of digital aerotriangulation in the works [17,18]. The PPP method was used to determine the UAV position during the photogrammetric flight. The obtained PPP positioning accuracy was higher than ±0.3 m in [17] and higher than ±1.2 m in [19] in the ENU coordinate frame. The use of the PPP method in air navigation allowed for combining navigation solutions using different sensors, e.g., a GNSS receiver and an INS (Inertial Navigation System) system. The works [19][20][21][22] presented a combined solution of aircraft positioning using the PPP/INS method in relation to the differential GPS RTK solution. The obtained PPP/INS positioning accuracy was: higher than ±0.2 m in [19], higher than ±0.4 m in work [20] and higher than ±0.1 m in [21,22]. In paper [23], the IGS and CNES (Centre National d'Etudes Spatiales) Analysis Center products were applied in aerial triangulation and, also, to determine the coordinates of an airplane.
During in-flight tests in Poland, the PPP measurement technique has not been widely used as an aircraft flight reference position so far. An example of testing the vertical plane's actual flight trajectory was described in the paper [24]. In the previous flight experiments in Poland, the differential RTK technique in the OTF mode was used as the aircraft reference position. Such a solution was proposed, and in-flight tests described, among others, in research papers [25][26][27]. Therefore, in Polish scientific research, there is a lack of confirmation and verification of the PPP measurement method as an alternative measurement technique to determine the aircraft's real trajectory. The use of the PPP method in flight tests in Poland has been limited only to testing the precision of determining aircraft coordinates [3,28,29]. Moreover, in works [30,31], the PPP method was used to determine the aircraft's precise trajectory to test the accuracy of SPP (Single-Point Positioning) code positioning in air navigation in Polish aviation. It is worth mentioning that the PPP method was used in air navigation in Poland for positioning using the GLONASS navigation system. Thus, paper [32] presented a model for determining the In turn, paper [33] presented the results of GLONASS positioning in the PPP method in relation to the joint GPS/GLONASS solution in the RTKLIB program. Furthermore, paper [34] presented the results of PPP positioning in GLONASS from RTKLIB and CSRS-PPP solution.
The research problem of using the PPP method in air navigation, discussed in this paper, concerns both real-time positioning and post-processing mode. The presented publications  presented research works both in real-time and post-processing mode. Therefore, the PPP positioning problem in air navigation is dual, i.e., it concerns realtime applications and the post-processing mode. The use of both is crucial and has its properties. In real-time, this is a navigation domain and allows the user to read real-time positions using the PPP method. The pilot must have confidence and credibility in the computer software that performs the calculations, and above all, the software must be tested and validated for air navigation. On the other hand, post-processing for the PPP method is primarily a post-factum analysis of the determined plane coordinates, reduction of systematic errors, detection of outliers, etc. All this is reflected in the objective flight control laboratory. Therefore, in the case of the PPP method, application tools must be designed for real-time and post-processing calculations.
Based on the accumulated literature on PPP positioning in air navigation, it can be observed that research studies have used a solution from a single GNSS navigation system such as GPS or GLONASS. This raises the research problem of what optimal computational strategy to use to improve a single GPS or GLONASS PPP solution relative to a RTK GNSS solution. Thus, an optimization algorithm should be used to improve the performance of PPP positioning in GPS or GLONASS in air navigation. One way to solve this problem is to create a linear combination of the determined coordinates from the PPP GPS and PPP GLONASS methods. This solution will determine the relationship between individual PPP GPS and PPP GLONASS solutions in relation to the RTK GNSS solution. Moreover, this approach shows how to choose the optimal weighting strategy to integrate GPS/GLONASS solutions in the PPP method. The weighted model gives an additional one degree of freedom to solve the plane's accidental position for independent GPS and GLONASS data. Moreover, in such an approach, measurement scales may be selected further to increase the air navigation PPP method's positioning accuracy. The weighted average model is characterized by high efficiency and effectiveness in improving coordinates' performance with appropriate selection of measuring weights [35].
The presented paper aims to show a new computational strategy for PPP positioning in air navigation for a single GPS and GLONASS solution. For this purpose, the development of a new boundary condition combining single GPS and GLONASS solutions in the PPP method is proposed. This condition is based on the weighted mean model. The proposed algorithm was tested on GPS and GLONASS kinematic data from an airborne experiment carried out with a Seneca Piper PA34-200T aircraft at a civil airport in Mielec, Poland. The presented numerical solution was performed in Scilab based on the aircraft position results obtained from CSRS-PPP. The presented results are unique for the application of the PPP method in Polish aviation.

Research Method
This section shows the PPP positioning algorithm for a single GPS and GLONASS solution. Additionally, a new PPP positioning algorithm improving the GPS and GLONASS single solution is presented.

Positioning Model of the PPP Method for GPS and GLONASS Solutions
The basic observation equation of the PPP measurement technique can be written as follows [33,34,36]: for GPS or GLONASS: where: P 3 and L 3 -observation equations of the GPS or GLONASS linear Ionosphere-Free combination; P 3 = α·P 1 + β·P 2 ; L 3 = α·L 1 + β·L 2 ; (P 1 , P 2 )-GPS or GLONASS code observations; (L 1 , L 2 )-GPS or GLONASS phase observations; , GPS or GLONASS linear coefficient; (f 1 , f 2 )-frequencies of the GPS or GLONASS signals; d-geometric satellite-receiver distance in GPS or GLONASS system; takes into account the phase centre characteristics of the satellite antenna and the receiver [37], The unknown parameters in the PPP measurement technique in the GPS and GLONASS solutions are, respectively [38], residuals to the approximate coordinates of the receiver [δX, δY, δZ], the receiver's clock correction dtr, phase ambiguities B 3 (determined for each visible GPS and GLONASS satellite) and the wet component of the tropospheric delay ZWD. The above-mentioned unknown parameters are determined by the least-squares method in a sequential stochastic process in GPS and GLONASS systems, respectively, as follows [39]: where: Sx-vector of the unknown parameters, A-design matrix, C n -the process disturbance noise matrix, l-vector of difference between observations and modelled parameters.
The final aircraft coordinates from the PPP method from the GPS or GLONASS solutions are determined as recorded below [40]: where:

New Mathematical Scheme for Improving the GPS and GLONASS Solution in the PPP Method
The new concept of PPP positioning from the GPS and GLONASS solution was based on the resultant linear combination model in which weighting factors (a, b, c, d, e, f ) were used. In the new approach, the aircraft's resultant position is determined in geocentric XYZ coordinates based on a single PPP GPS and PPP GLONASS solution. The new algorithm's goal is to determine the best match between the PPP GPS and PPP GLONASS solution relative to the flight reference position calculated from the differential RTK technique. Moreover, the proposed algorithm determines the most optimal solution of the resultant position of the aircraft, i.e., it will be checked which measuring scales improve the performance of determining the resultant aircraft coordinates. The proposed calculation scheme of the resultant position of the aircraft is described below: where: X g = X r i (see Equation (3)), position of aircraft from GPS solution along X axis, Y g = Y r i (see Equation (3)), position of aircraft from GPS solution along Y axis, Z g = Z r i (see Equation (3)), position of aircraft from GPS solution along Z axis, X r = X r i (see Equation (3)), position of aircraft from GLONASS solution along X axis,  (3)), position of aircraft from GLONASS solution along Y axis, Z r = Z r i (see Equation (3)), position of aircraft from GLONASS solution along Z axis, (a, b, c, d, e, f )-determined linear coefficients, dimensionless values, (X k , Y k , Z k )-resultant position of aircraft. Equation (4) combines 2 independent GPS and GLONASS determinations of the aircraft position using the PPP method. In this way, the aircraft's resultant position is determined every given measurement epoch, typically in the time interval of 1 s. It should be noted that the critical parameters in Equation (4) are the values of linear coefficients. In the analyzed example, a weighting scheme based on the following weighting factors was proposed: where: ns g -number of GPS satellites, ns r -number of GLONASS satellites, Cl g -mean error of the determined coordinates from the PPP GPS solution, Cl r -mean error of the determined coordinates from the PPP GLONASS solution.
The calculation process of determining the resultant position of the aircraft according to Equation (4) uses weighting factors (a, b, c, d, e, f ). In the analyzed model, two approaches for determining the weighting factors (a, b, c, d, e, f ) were proposed. The first approach is based on a weighing scale as a function of the reciprocal of the number of GPS or GLONASS satellites. On the other hand, the second approach is based on the measurement weight as a function of the inverse of the square of the mean coordinate errors of the aircraft from the PPP GPS and PPP GLONASS solutions. On this basis, it is possible to optimize the calculation process of the aircraft's resultant position. It means that based on Equation(4), the resultant plane position for the PPP method will be calculated, and by the way, it will be shown which measurement weight is the best for this type of mathematical model, as in Equation (4). That is important, because we are looking for mathematical models in air navigation to improve the aircraft's position. All the more, determining the appropriate measurement weight here is crucial to improving the accuracy of PPP positioning in aeronautical navigation.
The algorithm for determining the values in the Cl g and Cl r parameters is described below [41]: where: (mx g , my g , mz g )-mean coordinate errors determined from the GPS variance-covariance matrix C x (see Equation (2) f rom PPP GPS solution (see Equation (2)), (mx r , my r , mz r )-mean coordinate errors determined from the GLONASS variancecovariance matrix C x (see Equation (2) f rom PPP GLON ASS solution (see Equation (2)).
In the next step, additional parameters for the resultant position of the aircraft were determined, i.e., -determination of the average measurement error, -determining the corrections, -determination of the mean error of the arithmetic mean. The scheme of the procedure for this stage was written in Equation (7): where: P g -matrix of weights for PPP GPS solution, P g = 1 ns g ∨ P g = 1 Cl 2 g (see Equation (5)), P r -matrix of weights for PPP GLONASS solution, P r = 1 ns r ∨ P r = 1 Cl 2 r (see Equation (5)), V g -all residuals between PPP GPS solution and resultant aircraft position, V g,x -residuals between PPP GPS solution and resultant aircraft position along X axis, V g,y -residuals between PPP GPS solution and resultant aircraft position along Y axis, V g,z -residuals between PPP GPS solution and resultant aircraft position along Z axis, V r -all residuals between PPP GLONASS solution and resultant aircraft position, V r,x -residuals between PPP GLONASS solution and resultant aircraft position along X axis, V r,y -residuals between PPP GLONASS solution and resultant aircraft position along Y axis, V r,z -residuals between PPP GLONASS solution and resultant aircraft position along Z axis, n-number of solutions, n = 2, (σ x , σ y , σ z )-mean errors of measurement along XYZ axis, (M x , M y , M z )-mean error of the arithmetic mean along XYZ axis.
Then, the calculations are verified in the form of the global chi-square test, in accordance with the algorithm [41]: where: x · (n − 1), σ 2 y · (n − 1) and σ 2 z · (n − 1) must be equal or less than the tabular value of χ 2 f ,1−α . If it did not happen, then the test of internal reliability of calculations is not fulfilled, and the algorithm with the confidence level of 0.95 does not meet the requirements.
The last element of the calculations of the presented research method is the external control of the proposed algorithm, i.e., determining the accuracy of the presented algorithm, as shown below [35]: where: (DX, DY, DZ)-position errors, and (X RTK , Y RTK , Z RTK )-reference trajectory of aircraft based on RTK-OTF solution [35].
Position errors are determined independently along the XYZ coordinate axes for the entire flight path for each measurement epoch. Accuracy determination will confirm whether the proposed weighing scheme is optimal for improving the PPP GPS and PPP GLONASS positioning performances. Furthermore, the accuracy analysis will also show which weighting scheme is better, whether using weighting in the function of the number of satellites or using the inverse square of the error of the mean of the determined coordinates.
Considering the entire algorithm of the research methodology, it can include elements of calculations carried out in real time and post-processing mode. Starting from the beginning of the algorithm, the classic PPP method, i.e., Equations (1)-(3), can be implemented both in real time and post-processing. The boundary condition is, of course, the use of a dual-frequency dual-system GPS/GLONASS receiver. Next, the new proposed algorithm, i.e., Equations (4)-(8), can also be performed in real time and post-processing mode. However, it should be noted that, in Equation (4), there are measurement weights from Equations (5) and (6). In this case, parallel programming must be performed for both weighting processes. Of course, this weighting condition should first be tested in postprocessing to determine the best measurement weight-fitting model for GPS + GLONASS data integration. Equation (9), in turn, determines the aircraft positioning accuracy for the PPP measurement technique. The algorithm allows for numerical analysis in postprocessing mode. However, it is possible to perform it in real time, but it involves sending RTK corrections from the reference station to the GNSS on-board receiver, which requires an appropriate internet connection, transmission technique and data transmission power. Therefore, it is a difficult task in real time, when the aircraft's position changes dynamically every 1 s.

Research Test
To analyze the application of the presented research method in aviation, GPS/GLONASS satellite data were used from a flight test performed with a Seneca Piper PA34-200T aircraft at the civil airport EPML (European Poland Mielec) in Mielec, Poland. The test flight took place in the morning and lasted from 12:47:50 to 14:02:24, according to GPST time. The horizontal trajectory of the test flight is shown in Figure 1. A Topcon HiperPro dual-frequency dual-system geodetic receiver was deployed on-board the Seneca Piper PA34-200T aircraft to record GPS/GLONASS code-phase observations on the L1/L2 frequency. The Topcon HiperPro receiver was located in the cockpit. The Topcon HiperPro receiver recorded GPS/GLONASS code-phase observations with an interval of 1 s. The flight test was performed in real time. The receiver was collecting real-time GPS/GLONASS data for the numerical analysis.
The stage of the first scientific research for this paper was carried out in geodetic software CSRS-PPP v.3 [42]. CSRS-PPP software used the PPP absolute positioning method to determine the coordinates of the aircraft. The position of the Seneca Piper PA34-200T aircraft was determined from the GPS and GLONASS solutions and expressed in geocentric coordinates, according to Equations (1)-(3). The IGS products, such as precision ephemeris, precision satellite clocks and satellite and receiver antenna phase center characteristics, were used in the calculations. In addition, the P1/P2 codes and L1/L2 phase observations in GPS and GLONASS from the Topcon HiperPro on-board receiver were used in the calculations. The calculations assumed a pseudo-range measurement error of 0.3 m in GPS and 0.6 m in GLONASS and a phase measurement error of 0.003 m in GPS and 0.006 m in GLONASS. The elevation mask was 10 degrees for the GPS and GLONASS observation cut-off. For the stochastic modeling of the determined parameters, the following settings were assumed [38,39]: the coordinates of the aircraft are modeled from white noise, -the receiver clock offset is modeled from white noise, -the phase ambiguity for each satellite is modeled as a constant value; -the ZWD is modeled as a random walk parameter.
Other model parameters, such as the multipath effect, tidal and geodynamic corrections and phase slip, were used in the calculations.
The second part of the research was carried out in Scilab v.6.0.0 [43] for Equations (4)(5)(6)(7)(8)(9). In the Scilab program, the Seneca Piper PA34-200T aircraft's position coordinates were imported from the GPS and GLONASS solutions, and then, the alignment model was set up according to Equation (4). The final calculations in Scilab were performed for mathematical Equations (4)- (9). Calculations in Scilab were performed for all readings of the aircraft position from the GPS and GLONASS solutions. The reference position of the aircraft was calculated using the RTK-OTF solution in Trimble Total Control v.2.700 software. A Topcon HiperPro dual-frequency dual-system geodetic receiver was deployed on-board the Seneca Piper PA34-200T aircraft to record GPS/GLONASS code-phase observations on the L1/L2 frequency. The Topcon HiperPro receiver was located in the cockpit. The Topcon HiperPro receiver recorded GPS/GLONASS code-phase observations with an interval of 1 s. The flight test was performed in real time. The receiver was collecting real-time GPS/GLONASS data for the numerical analysis.
The stage of the first scientific research for this paper was carried out in geodetic software CSRS-PPP v.3 [42]. CSRS-PPP software used the PPP absolute positioning method to determine the coordinates of the aircraft. The position of the Seneca Piper PA34-200T aircraft was determined from the GPS and GLONASS solutions and expressed in geocentric coordinates, according to Equations (1)-(3). The IGS products, such as precision ephemeris, precision satellite clocks and satellite and receiver antenna phase center characteristics, were used in the calculations. In addition, the P1/P2 codes and L1/L2 phase observations in GPS and GLONASS from the Topcon HiperPro on-board receiver were used in the calculations. The calculations assumed a pseudo-range measurement error of 0.3 m in GPS and 0.6 m in GLONASS and a phase measurement error of 0.003 m in GPS and 0.006 m in GLONASS. The elevation mask was 10 degrees for the GPS and GLONASS observation cut-off. For the stochastic modeling of the determined parameters, the following settings were assumed [38,39]: the coordinates of the aircraft are modeled from white noise, -the receiver clock offset is modeled from white noise, -the phase ambiguity for each satellite is modeled as a constant value; -the ZWD is modeled as a random walk parameter.
Other model parameters, such as the multipath effect, tidal and geodynamic corrections and phase slip, were used in the calculations.
The second part of the research was carried out in Scilab v.6.0.0 [43] for Equations (4)-(9). In the Scilab program, the Seneca Piper PA34-200T aircraft's position coordinates were imported from the GPS and GLONASS solutions, and then, the alignment model was set up according to Equation (4). The final calculations in Scilab were performed for mathematical Equations (4)- (9). For the purposes of scientific research in aviation, the state of the ionosphere was determined by the VTEC (Vertical TEC) parameter. The VTEC value was determined based on the GIM (Global Ionosphere Maps) ionosphere model from the CODE (Center for Orbit Determination in Europe) Analysis Center in Switzerland [44]. The VTEC parameter fluctuated during the flight tests from 16.8 TECU to 17.3 TECU.

Results
The presentation of the research results began with determining the mean reference position errors from the RTK-OTF solution. Figure 2 shows the results of the mean reference flight position errors. The mean errors were resolved as in Equation (6), and these values are the resultant mean errors for the XYZ coordinates. The RTK-OTF solution used precise GPS and GLONASS phase observations at the L1/L2 frequencies. Analyzing the obtained results, the resultant aeroplane position's mean errors were from 0.005 m to 0.063 m. Over 52% of the mean error results were less than 0.02 m. In turn, over 84% of the results of the mean error were less than 0.03 m. For the purposes of scientific research in aviation, the state of the ionosphere was determined by the VTEC (Vertical TEC) parameter. The VTEC value was determined based on the GIM (Global Ionosphere Maps) ionosphere model from the CODE (Center for Orbit Determination in Europe) Analysis Center in Switzerland [44]. The VTEC parameter fluctuated during the flight tests from 16.8 TECU to 17.3 TECU.

Results
The presentation of the research results began with determining the mean reference position errors from the RTK-OTF solution. Figure 2 shows the results of the mean reference flight position errors. The mean errors were resolved as in Equation (6), and these values are the resultant mean errors for the XYZ coordinates. The RTK-OTF solution used precise GPS and GLONASS phase observations at the L1/L2 frequencies. Analyzing the obtained results, the resultant aeroplane position's mean errors were from 0.005 m to 0.063 m. Over 52% of the mean error results were less than 0.02 m. In turn, over 84% of the results of the mean error were less than 0.03 m. In the next step, the measurement weights, i.e., the values of the coefficients  In the next step, the measurement weights, i.e., the values of the coefficients (a, b, c, d, e, f ) used in Equation (4), are shown. The weighting factors' values were calculated based on Equation(5). In Figure 3, the values of the coefficients (a, b, c, d, e, f ) for measuring the weights (a = c = e = 1 ns g = P g ; b = d = f = 1 ns r = P r ) are shown. For the solution of the PPP GPS position, the measurement weight was from 0.111 to 0.200. On the other hand, for the PPP GLONASS position solution, the measuring balance results were from 0.143 to 0.200. You can see that the measurement weights were not the same as the numbers of the GPS and GLONASS satellites, which also varied during the experiment.
Moreover, Figure 4 shows the results of the coefficient values (a, b, c, d, e, f ) for measuring weights (a = c = e = 1       The next test results concerned the distribution of corrections (V g,x , V g,y , V g,z , V r,x , V r,y , V r,z ) for the XYZ components. be seen that we deal with the results that best correspond to the expected values of the Z components. The V g,z corrections distribution is from +0.097 m to +0.201 m and V r,z corrections from −0.284 m to +0.026 m, respectively. m. The distribution of corrections along the Z-axis is very interesting. When we compare the results of corrections between Figures 5-7, we can see at once that the distributions in Figure 7 are an order of magnitude smaller than in Figures 5 and 6. That causes the scatter of the Z coordinate obtained results as the smallest around the expected value. Especially when we look at the distribution of corrections in Figure 7, it can be concluded that their dispersion is from −0.284 m to +0.201 m. It can be seen that we deal with the results that best correspond to the expected values of the Z components. The    , r y m. The distribution of corrections along the Z-axis is very interesting. When we compare the results of corrections between Figures 5-7, we can see at once that the distributions in Figure 7 are an order of magnitude smaller than in Figures 5 and 6. That causes the scatter of the Z coordinate obtained results as the smallest around the expected value. Especially when we look at the distribution of corrections in Figure 7, it can be concluded that their dispersion is from −0.284 m to +0.201 m. It can be seen that we deal with the results that best correspond to the expected values of the Z components. The g y r y V V along the Y-axis. Figure 6. Residuals (V g,y ; V r,y ) along the Y-axis.
Then, in Figures 8-10, the results of the calculated mean errors (σ x , σ y , σ z ) are presented separately for the resultant XYZ coordinates of the aircraft position. It should be added that the parameter values (σ x , σ y , σ z ) are shown for two weight models, i.e., P g = 1 ns g ∨ P g = 1 Cl 2 g and P r = 1 ns r . When observing Figures 8-10, it can be seen that, when we use the weight in the inverse of the mean error function, we get better results than when we weight in the joint function of the number of tracked GPS or GLONASS satellites. Therefore, the values σ x are from 0.246 m to 0.323 m when using measuring weights P g = 1 ns g and P r = 1 ns r . On the other hand, when we use measuring weights P g = 1 ) compared to weights (P g = 1 ns g ; P r = 1 ns r ) in the calculations. The mean errors σ y were from 0.571 m to 0.769 m when using weights P g = 1 ns g and P r = 1 ns r . On the other hand, when we used measuring weights P g = 1 by about 65% when we used weights (P g = 1 Cl 2 g ; P g = 1 Cl 2 r ) compared to the application of (P g = 1 ns g ; P r = 1 ns r ) in the calculations. The mean errors σ z had values from 0.001 m to 0.089 m when using weights P g = 1 ns g and P r = 1 ns r . When we used measuring weights P g = 1 Cl 2 g and P r = 1 Cl 2 r , then the σ z parameters ranged from 0.001 m to 0.029 m. In this case, it can also be said that the use of measuring weights (P g = 1 Cl 2 g ; P g = 1 Cl 2 r ) reduced the mean errors σ z by over 65% compared to the solution using (P g = 1 ns g ; P r = 1 ns r ).   The test results in Figures 11 and 12 show the (M x , M y , M z ) parameter values for the solution using measuring weights (P g = 1 ns g ; P r = 1 ns r ) and (P g = 1 ) improved by about 65% to 66% compared to the corresponding results for these quantities when weights (P g = 1 ns g ; P r = 1 ns r ) were used.       . Figure 11. Mean errors of the arithmetic means (M x , M y , M z ) from the adjustment scheme (P g = 1 ns g ; P r = 1 ns r ).    . Table 1 shows the results of the σ 2 0 · (n − 1) chi-square test's parameters and table values χ 2 f ,1−α at the 1 − α = 0.95 confidence level and for f = n − 1 degrees of freedom. The parameter results σ 2 0 · (n − 1) were smaller than the tabulated χ 2 f ,1−α values, so the global chi-square statistical test, which determines the internal reliability of the performed alignment, was fulfilled and implemented. In the analyzed case, the tabular value of the chisquare test was 3.841. Value σ 0 denotes parameters (σ x , σ y , σ z ), respectively. Moreover, the number of degrees of freedom f = n − 1 is equal to 1. The chi-square test in the analyzed example was carried out and checked for each measurement epoch. The chi-square test was fulfilled for all measurement epochs, making the obtained test results more reliable in the experiment.
The last stage of the research concerned determining the accuracy of the presented research method. In this, Equation (9) from the research methodology was used. The position errors DX were from +0.032 m to +0.527 m for weights (P g = 1 ns g ; P r = 1 ns r ). On the other hand, for weights (P g = 1 Cl 2 g ; P g = 1 Cl 2 r ), the position errors DX were from −0.244 m to +0.206 m. It could be observed that the positioning accuracy along the X-axis improved by about 87% in the case of using weights (P g = 1 Cl 2 g ; P g = 1 Cl 2 r ) in relation to the results of DX determined for the (P g = 1 ns g ; P r = 1 ns r ) weighing process. The DY position errors were from −1.051 m to −0.492 m for (P g = 1 ns g ; P r = 1 ns r ) weighing. On the other hand, for (P g = 1 ) in relation to the DY results determined for the (P g = 1 ns g ; P r = 1 ns r ). The DZ position errors were from −0.451 m to +0.183 m for weights (P g = 1 ns g ; P r = 1 ns r ). In turn, for (P g = 1 Cl 2 g ; P g = 1 Cl 2 r ), the DZ position errors were −0.396 m to +0.228 m. It is worth noting that the positioning accuracy along the Z-axis improved by more than 11% when weights (P g = 1 were used in relation to the results determined for weights (P g = 1 ns g ; P r = 1 ns r ) (Please see Figures 13-15).

Discussion
In the discussion of the results, the attention is mainly focused on determining the proposed test method's external reliability. For this purpose, the proposed algorithm was verified in terms of determining positioning accuracy. In the discussion, the aircraft po-

Discussion
In the discussion of the results, the attention is mainly focused on determining the proposed test method's external reliability. For this purpose, the proposed algorithm was verified in terms of determining positioning accuracy. In the discussion, the aircraft position's obtained results using the PPP method from the RTKLIB program [45,46] for the weighted average model used were analyzed. Such a check is crucial, because we obtain independent confirmation of the proposed mathematical algorithm's correctness (4)(5)(6)(7)(8)(9). In particular, in the proposed research method's external control, the key factor is Equation (9), thanks to which, the position errors can be determined [47]. Therefore, the weighted average model results will be compared to the reference coordinates determined by RTK-OTF. Figures 16-18 show the positioning accuracy between the weighted average model and the RTK technique as a function of the measurements used. The DX position errors were from +0.157 m to +1.668 m for weights (P g = 1 ns g ; P r = 1 ns r ). On the other hand, for (P g = 1        The DY position errors were from −0.745 m to +0.812 m for (P g = 1 ns g ; P r = 1 ns r ) weighing. In turn, for (P g = 1 for weights (P g = 1 ns g ; P r = 1 ns r ). The DZ position errors ranged from −2.021 m to −0.228 m for (P g = 1 ns g ; P r = 1 ns r ). On the other hand, for (P g = 1 Cl 2 g ; P g = 1 Cl 2 r ) weighing, the position errors were −1.604 m to +0.295 m. It is worth noting that the positioning accuracy along the Z-axis improved by over 59% when (P g = 1 Cl 2 g ; P g = 1 Cl 2 r ) weighting was used in relation to the DZ results determined for the (P g = 1 ns g ; P r = 1 ns r ) weighing process. The aircraft positioning results obtained during the research using the PPP method for the GPS and GLONASS solutions showed the high quality and accuracy of the moving object's determined coordinates. The presented research method significantly improved the accuracy of PPP GPS and PPP GLONASS positioning in relation to the RTK solution. The obtained research results showed that the applied computational scheme of the mathematical model fits well in the modern trend of scientific research on the application of the PPP method in air navigation, as in other works . That is important, because the developed computational algorithm improved the performance of position determination and, at the same time, the accuracy of PPP positioning, as in works . It should be mentioned that the research problem of improving the performance of the PPP method in air navigation is a current and extremely interesting topic; hence, there have been a large number of research papers  undertaken by scientists worldwide. The topic should be constantly developed in the aviation industry, especially in Poland, which needs modern solutions in implementing GNSS systems in aviation. The use of PPP measurement technology in the positioning of aircrafts can accelerate research on the determination of the quality of GNSS satellite positioning in aviation, especially the positioning accuracy parameter. The submitted research work shows that modern navigation solutions in Polish aviation are needed and require implementation.

Conclusions
The paper presented the results of research on the application of a new research method for improving the performance of computing aircraft coordinates from GPS and GLONASS solutions for the PPP measurement technique. The paper addressed the problem of improving GPS and GLONASS positioning accuracy in the PPP method by developing a new computational algorithm based on a linear combination of individual GNSS solutions, i.e., GPS and GLONASS. The algorithm used a weighted mean model between the XYZ geocentric coordinates from the GPS and GLONASS. The mathematical relationship was based on the determination of six linear coefficients (a, b, c, d, e, f ), which linked the XYZ position coordinate readings together. The whole computational scheme was solved by the weighted mean model, taking into account different measurement weights. The measurement weights were a function of the number of GPS and GLONASS satellites and, also, as the inverse of the square of the mean error. Moreover, the whole algorithm was tested on GPS and GLONASS kinematic data from an airborne experiment carried out with a Seneca Piper PA34-200T aircraft at the Mielec airport. A dual-frequency dual-system GPS/GLONASS receiver was placed on-board the aircraft, which made it possible to record the GNSS observations, which were then used to calculate the aircraft's position in CSRS-PPP software. The calculated XYZ position coordinates from the CSRS-PPP software were used in the developed optimization algorithm and implemented in the Scilab programming language. The accuracy of the presented research method was checked and verified with the results obtained from the RTK-OTF solution. Moreover, the proposed solution's accuracy was better by 11-87% for the model with the weighting scheme as a function of the inverse of the square of the mean error. The performance of the algorithm was also tested for results from another navigation software, i.e., RTKLIB. In this case, the algorithm's accuracy for the improvement of PPP positioning was better by 45-82% for the model with the weighting scheme as a function of the inverse of the square of the mean error. The presented method's effectiveness was relatively high and could be used to improve the performance of GPS and GLONASS in PPP positioning method for air navigation.