An Improved QZSS Satellite Clock Offsets Prediction Based on the Extreme Learning Machine Method

The Japanese Quasi-Zenith Satellite System (QZSS) has been developed as a GPS (Global Positioning System) complementary system to improve positioning accuracy in the Asia-Pacific region. However, the accuracies of ultra-rapid predicted clocks are not high enough for real-time applications, so it is still a challenge problem. This article focuses on the clock predictions of the new QZSS constellation. Based on the QZSS clock periodic characteristics, an improved clock prediction method combining the spectrum analysis model (SAM) and extreme learning machine (ELM) is proposed with abbreviation as iELM. The key parameters of iELM are selected carefully including the number of hidden layer nodes and activation function. Further, the input length of the iELM network is optimized and discussed thoroughly. For the purpose of assessing the performance of the proposed algorithm, the clock prediction accuracies are compared among GBU-P (the ultra-rapid predicted orbits/clocks provided by GFZ (Deutsches GeoForschungsZentrum)), SAM and iELM methods. It is demonstrated that iELM keeps in a high level of accuracy below 1.0ns with the predicted time increasing from 0 to 18h, and a little larger during the next six hours. And the QZO (Quasi-Zenith Orbit) satellites perform better than GEO (Geostationary Earth Orbit) satellite. Furthermore, precise point positioning (PPP) for both static and kinematic modes are experimentally studied for 13 IGS (International GNSS (Global Navigation Satellite System) Service) MGEX (Multi-GNSS Experiment) stations in the longitude range between 100°E and 180°E. In the static PPP mode, the iELM method is verified to be effective for the GPS/QZSS constellation as positioning accuracy is improved by 28.3%, 57.7% and 47.4% on average in the east, north and up component with respect to GPS/QZSS GBU-P results. Nearly 70.0% stations achieve sub-decimeter accuracy in the vertical component. As for the kinematic PPP, the iELM method based on GPS/QZSS observations performs much better than the others for shorter convergence time and better positioning accuracy. Compared with GPS/QZSS GBU-P, iELM takes at most a half time to get convergence, and the accuracy is enhanced by 27.6%, 23.7% and 13.9% on average in the east, north and up direction respectively.


I. INTRODUCTION
Global Navigation Satellite System (GNSS) signals with high elevation angles can help to reduce atmospheric errors, multipath errors, and other errors in positioning [1]. The GPS (Global Positioning System) satellites are evenly distributed The associate editor coordinating the review of this manuscript and approving it for publication was Woorham Bae . in six orbital planes, so the number of satellites that can be seen at the same time is not sufficient in some parts of the earth. So GLONASS (Global Navigation Satellite System) is designed with a bigger inclination of 64.8 • to cover more Russian territory with better signal quality. In Japan, QZSS (Quasi-Zenith Satellite System) has been developed as a GPS complementary system to increase the availability, reliability, integrity and accuracy of positioning in the VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Asia-Pacific region [2], especially in urban canyons and mountainous areas [3]. Now, QZSS is a four satellite constellation [4], and is planned to be extended to a seven satellite system by 2023 [5]. Except for J07, which is a GEO (Geostationary Earth Orbit) satellite, the other three QZSS satellites are in QZOs (Quasi-Zenith Orbits) with an inclination of around 43 • [6]. The QZSS satellites send GPS compatible signals on L1 (1575.42 MHz), L2 (1227.6 MHz), L5 (1176.45 MHz) and augmentation signals [7], i.e. the L1-SAIF (Submeter-class Augmentation with Integrity Function) signal and the LEX (L-band Experimental) signal [8]. Therefore, QZSS can be applied in an integrated way with GPS to obtain more accurate positioning results [9], [10]. Fig.1 gives the QZSS satellites visibility. The up subfigure shown in Fig.1 indicates the percentage of a 24h period with at least one QZSS satellite visible at an elevation angle of 70 • . It is noted that this percentage indeed gets 100% in Japan. The bottom subfigure in Fig.1 shows that there is an olive-shaped area in the Asia-Pacific region, indicating all the four QZSS satellites being visible simultaneously where a cut-off elevation angle of 7 • is considered. Further, the average number of visible GPS/QZSS satellites in the Asia-Pacific region is increased from 11 to 15 [11]. Here, the elevation angle of 7 • is an empirical value [12]. If an elevation angle is smaller than 5 • , the atmospheric errors and multi-path errors are very large, which may significantly degrade the quality of received signal. When a high elevation angle is selected, e.g. more than 15 • , the geometry of satellites is not so good, and the number of lost observations will increase. The higher accuracies of positioning and navigation are required in a high quality application [13]. Normally, the GNSS positioning performance is immediately decided by the accuracies of satellite orbit/clock products [14]. With the great efforts of some IGS (International GNSS Service) analysis centers, several types of precise orbit/clock products for multi-GNSS are provided to meet the accuracy demands for different applications [15]. Among these products, the ultra-rapid products (6h update) are generated and published for possible real-time applications by predicting orbits/clocks over 24h [16]. Fortunately, the ultra-rapid products for some new developed constellations, e.g., BDS (BeiDou Navigation Satellite System) and QZSS, are already available in some analysis centers, e.g. GFZ (Deutsches GeoForschungsZentrum) [17]. However, for a newly built constellation, at the very beginning stage of generating the ultra-rapid products, the accuracy of products should be verified by lots of positioning experiments. It is also well known that the ultra-rapid orbits have an internal consistency of better than 5.0cm, which can satisfy the requirements of realtime applications. However, the existing clock extrapolation methods would incur serious accuracy loss because of the discrete nature of satellite clock errors [18]. The ultra-rapid clock accuracy is not well enough, and it is the key problem of the current real-time precise point positioning (PPP) [19].
There are several traditional methods to predict satellite clock for both short-term and long-term [20], such as quadratic polynomial, grey model, spectrum analysis, and Kalman filter [21]. However, many disadvantages are exposed for these methods. For example, the errors of quadratic polynomial are accumulated continuously along with the prediction time [22]. For the Kalman filter, a better prediction can be obtained after a suitable stochastic model for the clock offset series is fixed. With the rapid development of computer and information technologies, the artificial intelligence becomes widely used in the prediction fields [23]. Some contributions focus on improving the clock prediction using artificial intelligence algorithms [24], for example, Back Propagation Neural Networks (BPNN) [25], Generalized Regression Neural Network (GRNN) [17], and Support Vector Machine (SVM) [26]. Based on these methods, feasible clock prediction results have been obtained for GPS and BDS with optimized parameters [27]. However, the computing complexity increases to some extent, which is important for real-time applications. To avoid large computation burden and empirical parameter setting, this article will propose a new and efficient prediction method based on the Extreme Learning Machine (ELM), and it will make the QZSS clock prediction to be more accurate.
It is also well known that the satellite clock offset is very complex, since its complicated intrinsic physical characteristic and the diverse space environment. The most distinct feature of clock is believed as its periodic characteristic, which should be studied firstly as prior information for the improved clock prediction. Commonly, a spectrum analysis with one or two periodic terms is applied to generate clock predictions for GPS and BDS [28]. For example, as for GPS [29], and BDS [30], the most pronounced period is 12h or 24h. Interestingly enough, some special periodic terms can also be detected from the large clock sample set with the aid of spectrum analysis, e.g. of about 1.3h for BDS-2 C11, which is caused by its separately hardware of the atomic oscillator. Similarly, satellite-specific periodic characteristic for QZSS satellite clock is a prerequisite for successful clock prediction.
After introduction for the QZSS constellation and clock prediction background in Section I, the statistic characteristics for the QZSS clock are studied in Section II, focusing on periodic terms determined by Fast Fourier Transform (FFT) method. Then, an improved clock prediction method employing ELM algorithm for the QZSS satellite is proposed, which is called as iELM for short in Section III. The performances of the iELM method are investigated to obtain clock prediction accuracy. Further, the static and kinematic PPP results are analyzed in Section IV. Finally, Section V concludes with a discussion.

II. PERIODIC CHARACTERISTIC ANALYSIS A. DATA COLLECTIONS
The QZSS constellation consists of three QZOs (J01, J02 and J03) and one GEO (J07), which are equipped with the same Rubidium Atomic Frequency Standard (RAFS) as GPS block IIF satellite [31]. There are several analysis centers providing the QZSS orbit/clock products, e.g. GFZ, CODE (Center for Orbit Determination in Europe), TUM (Technische University München), JAXA ( Japan Aerospace Exploration Agency), and WUH (Wuhan University) [31]. First QZSS archive products for Satellite Positioning, Navigation and Timing Service have been available since April 14, 2017 by Japan National Space Policy Secretariat Cabinet office. Studies show that the number and quality of products are different between each analysis center. For example, JAXA just provides products of J01, while GFZ contributes both final precise products and ultra-rapid ones including J01, J02, J03, as well as J07, of which precise clocks were available from September in 2018. Xie analyzed the number of clock samples among these analysis centers, and indicated that the integrity of the samples from 2017 to 2019 provided by GFZ is ensured [31]. Hence, GFZ orbit/clock products are the best choice in this article for the characteristics research. GBU products (the ultra-rapid observed orbits/clocks provided by GFZ) are collected from GFZ servers for nearly one year from DOY (Day of Year) 001 to 340 in 2019 to analyze and evaluate the QZSS satellite clocks characteristics and performances.

B. CLOCK OFFSET PREPROCESSING
There are lots of factors, such as space environment, clock switch, solution interrupt, and orbit maneuver, impacting to the performance of onboard atomic clocks and the quality of clock products, which lead to kinds of abnormal data, e.g., gross errors and outliers. Hence, before clock products analysis and prediction, preprocessing has to be carefully conducted to get rid of abnormal data and keep data continuous. The epoch-differenced clock series are obtained ahead as frequency series, then a robust but simple method, namely Median Absolute Deviation (MAD) is applied for blunders detection. The MAD can be expressed as follows, where f i is the epoch-differenced clock series with index i, m is the median of the epoch-differenced clock series. Then, the abnormal clock offset will be deleted when |f i | > (m+ n× MAD) is valid. n is an empirically integer usually setting as 3 in this processing. To avoid losing too much useful information, no more than 5.0% of total data are eliminated in the practice [17]. Preprocessing is employed for each QZSS satellite clock series using the GBU products with 300s sample interval for more than 340 days, and the percentages of data rejection for J01, J02, J03 and J07 are 0.3%, 0.6%, 1.1% and 2.8%, respectively. Data interpolating is not used here to avoid manmade influence on clock series, thus detected oultliers are got rid directly and the interrupts are dealt piecewise.

C. PERIODIC CHARACTERISTICS
Since the satellite orbits and clocks are coupled tightly, significant periodic characteristics are found in satellite clocks, which also reflect the performance of the satellite orbits. In order to identify the potential satellite-specific periodic terms, the following two steps are conducted. Firstly, daily fitting residuals are obtained to eliminate trend terms using quadratic polynomial method. Secondly, after relativity effects removed, the FFT method is used to obtain the main periods and amplitudes [17]. This task is conducted for the QZSS clocks using more than 300 days samples from DOY 001 to 340 in 2019. The relationships between amplitude and frequency are identified for each satellite given in Fig.2. It can be seen from Fig.2 that the pronounced two periodic terms are 24h (one cycle per revolution) and 12h (two cycles per revolution) for four satellites, which are related to the orbital operation and general relativistic effects. As for J02 and J07, the amplitude of 12h is larger than that of 24h. Hence, the main periodic term is 12h and the secondary period term is 24h. For J01 and J03, the amplitude of 8h is more obvious than that of 12h. Hence, their two main periodic terms are 24h and 8h. Further, the outstanding two periodic terms of each QZSS clock are listed in Table 1. Additionally, for J01 and J03, sharp pulses at each integer from 1 to 24 cycles per revolution are obviously seen from Fig.2, which indicate there are a series periods. However, the corresponding amplitudes are decreasing, with n increasing from 1 to 24. These obvious periods may be also caused by the hardware oscillation features of the atomic clock itself, and it is out of the research field of this article. Besides, since the clock periods are tightly related to the orbit solution, these periods in QZSS clocks can also be explained as some unknown orbit related errors are absorbed by satellite clocks when precise orbit/clock products are determined. Hence, the accuracy of clock prediction just using polynomial model cannot satisfy the requirements of real-time applications, the periodic terms should be considered as one of the most key influences on the realization of more precise clock prediction.

III. QZSS CLOCK PREDICTION A. SAM FOR CLOCK PREDICTION
Normally, a second-order polynomial is adopted as a basic model for GNSS clocks, especially for Rb atomic clocks [31]. Since the QZSS orbital errors are absorbed in the clock products, several periodic terms with pronounced amplitudes are detected, and it should be considered to build a more precise model [25]. Hence, a polynomial model is taken as the basic function, and one or two significant periodic terms are applied to model QZSS clocks, which is called spectrum analysis model (SAM for short) and expressed as follows, where a 0 is the clock bias correction coefficient, a 1 is the clock drift correction coefficient, and a 2 is clock drift rate correction coefficient. t 0 is the reference epoch of clock offset series. t i and clk(t i ) denote the i th epoch and its corresponding clock offset. l is the number of periodic terms. r and f r are the amplitude and frequency of periodic terms, respectively. φ r is its corresponding initial phase. ε i is the residual of the prediction model.

B. iELM METHOD FOR CLOCK PREDICTION
It is reported that the ultra-rapid products with the accuracy of about 2.0 -6.0ns is one of the most important products in decimeter level real-time positioning services [25]. However, evident nonlinear system signals are detected in ultra-rapid clock products. To satisfy the requirements of high precision application scenarios, such as landslide early warning and hydrographic surveying, an improved model is required to further enhance the clock prediction accuracy [17]. Thanks to the dramatic development of artificial intelligence, some new algorithms for GNSS clock prediction are proposed, and effective results for GPS and BDS clock prediction have been obtained. However, these artificial intelligence algorithms rely much on empirical parameters selection [26], and their learning speeds are much slower than the required, because many parameters have to be tuned iteratively. This is the noteworthy bottleneck in their applications for past decades.
To solve this problem, a fast learning algorithm called Extreme Learning Machine (ELM) is proposed by Huang [32]. The ELM has a lot of merits. Firstly, the connection weights between the input layer and the hidden layer can be set randomly, as well as the thresholds of the hidden layer. This is different from the conventional neural networks, which need to constantly adjust the weights and thresholds. Secondly, the connection weights between the hidden layer and the output layer are determined once by solving equations, rather than adjusted iteratively. Hence, the learning time of ELM can be reduced by half. Thirdly, the generalization performance of ELM is better than that of the conventional neural networks. Therefore, the ELM method is used to improve QZSS clock prediction in this article. As displayed in Fig.3, ELM is a feedforward neural network with three layers, including one input layer, one hidden layer and one output layer. Given N as arbitary distinct ∈ R m , ELM mathematical model with L hidden layer nodes can be written as follows, where is an activation function. w j = w 1j , · · · , w nj T is the weight vector connecting the j th hidden node and the input nodes. β j = β j1 , · · · , β jm T is the weight vector connecting the j th hidden node and the output nodes. b j is the threshold of the j th hidden node. L is the number of hidden layer nodes. The ELM network with L hidden layer nodes and activation function g (·) can approximate the N samples with zero error means that L o i − y i = 0, there exist β j , w j , and b j yield, Further, Equation (4) can be compactly written as, where H is the hidden layer output matrix with detail elements in Equation (6), and the j th column of H is the j th hidden node output with respect to inputs As soon as the number of hidden layer nodes L and the activation function g (·) are fixed, the process of ELM can be described simply in three steps. Firstly, the weight vector w j and the bias b j can be set up randomly, j = 1, · · · , L. Secondly, the output matrix H of the hidden layer can be calculated. Thirdly, the weight matrix β can be calculated by β = H † Y, where H † is the Moore-Penrose generalized inverse of the matrix H.
In order to better approach the remaining clock residuals after the removal of trend and periodic terms, an improved model (iELM for short) based on ELM and SAM methods is proposed. Since the ELM is a machine learning algorithm without explicit mathematical expression, the part of ELM is depicted as a function of f ELM . Hence, the improved model is proposed as, in which, parameters are the same as Equation (2). Further, the processing diagram of the improved model is given in Fig.4. The input data are clock offset series from GBU estimations. After preprocessing to get rid of outliers, clock characteristics are studied with the aid of spectrum analysis to obtain at least two significant periodic terms. The spectrum analysis residuals are considered as sample dataset for ELM. The commonly usage is a subset of sample data (70%) for training and the rest of sample data for testing to assess its generalization ability. Generally, in ELM training, w and b are generated randomly. As soon as nodes number is fixed in the hidden layer, β can be obtained with the help of theoretical calculation. This process is more simplified with faster training, compared with other neural networks [17].

C. OPTIMIZATION OF ELM PARAMETERS
It is well known that the parameters setting is crucial for a machine learning algorithm. Different hidden layer nodes number, activation function, and input length will result in different prediction. Since the ELM algorithm is a forward neural network with single hidden layer, so number of the hidden layer is known as one. According to the basic theory of ELM in Section III.B, the number of hidden layer nodes L and activation function are the key parameters in the ELM method. Experiments have proved that the prediction accuracy difference is no more than 10 −4 ns when changing L from 1 day (288 epochs) to 11 days (3168 epochs). However, the larger L is selected, the longer time is consuming for training. When L is smaller, the prediction variance is larger. So L is selected as 2 days (576 epochs) in this study. Meanwhile, experiments also carried out to investigate the impacts of activation functions, including radial basis, sine, and exponential. The results are consistency with the previous studies that ELM is always effective when the activation function is selected randomly [32]. Then, a normal used activation function, i.e., sigmoidal function is used in this study.
However, the input length has an evident influence on the clock prediction accuracy of the ELM method. Experiments indicate that the prediction accuracy is not always proportional to the input length, since the over fitting phenomenon in the training process of machine learning [33]. More data may introduce more noise, and the ELM network focuses much more on the complex details rather than the whole characteristic of the samples, resulting in decreased performances. Additionally, each atomic clock behaves unique characteristic, in order to optimize satellite-specific clock prediction, the input length is chosen carefully for each QZSS clock. The input length is varying from one day to eleven days, and the prediction time is the following day (24h). The statistic results are shown in Fig.5. The upper bound and the lower bound of the prediction accuracy is given as error bar and the average accuracy is presented as bar for each QZSS clock. It can be seen from Fig.5 that the best input length for training data is VOLUME 8, 2020

IV. EVALUATION OF THE PREDICTION MODEL
Generally, all of the satellite clock prediction accuracy can be evaluated by the deviation of prediction from its true value. However, the actual value of in-orbit clock is unknown. So a reference benchmark must be given. On the other hand, the predicted clock can be checked together with its related orbit by means of precise point positioning. In this article, both evaluations are carried out and the later one is more important, as it will put out a preliminary assessment of the positioning accuracy, which is the main goal of this research.

A. CLOCK COMPARISON
First of all, the GBU products are selected as the reference to study the performance of clock prediction. To demonstrate the effectiveness of the iELM model, the residuals between the predicted clocks and GBU clocks at each epoch with interval 300s are investigated and the results of DOY 282 in 2019 are shown as an example in Fig. 6. In the meantime, GBU-P (the ultra-rapid predicted orbits/clocks provided by GFZ) and SAM results are also plotted for comparison. In order to comply with IGS specification, the 24h extrapolated clocks with 300s sample interval are given in this article.
The reason why extrapolated to 24h is that IGS is a broad-cast service, a large amount of calculations are required to obtain the observed orbits/clocks, and the latency is 6h. Therefore, it is necessary to provide the predicted orbits/clocks for end-user applications in future hours. Additionally, in order to further improve the reliability of enduser applications, IGS extrapolates to 24h due to the possible interruption of communication [34].
As shown in Fig.6, the prediction errors gradually enlarge with the predicted time increasing when GBU-P and SAM are used. It is indicated that the clock prediction accuracy is of similar magnitude at the beginning 3h using both GBU-P and SAM, but the prediction errors are soon accumulated to a significant magnitude from 6h to 24h. It is evident that J01 behaves the worst among the four QZSS satellites.
The prediction residuals of J01 using GBU-P are more than 3.5ns when the prediction time is approaching 24h, while SAM obtains a little better accuracy as 2.5ns. Meanwhile, the other three QZSS satellites (J02, J03 and J07) predicted clocks using GBU-P have similar performances. The clock predictions are no more than 1.0ns when prediction time is within 18h. Comparably, as for the SAM method, there are two satellite clocks, i.e., J01 and J02, reaching to 1.0ns when the predicted time increases to 12h. The accuracy of SAM is up and down compared to GBU-P. The reason is GBU-P is optimized according to the satellite clock characteristics, while SAM employs two fixed periodic terms on a quadratic polynomial model. Comparatively, the performance is much better if the iELM algorithm is applied. The most important improvement is the prediction errors do not enlarge obviously from 0 to 18h. Moreover, the accuracies of predicted clocks keep in 1.0ns at 24h even for J01, which is the worst one using GBU-P and SAM. Fig.6 gives a straightforward and visualized example employing the extrapolated clocks on DOY 282 in 2019, which are the instantaneous predicted clocks using three methods, i.e. GBU-P, SAM and iELM. Although the similar performances can be found before 12h for GBU-P, SAM and iELM in Fig.6, it still can be seen that the clock residuals of the iELM method are smaller for shorter prediction times. In order to prove this result, statistical analysis is conducted, and the statistic result is given in Table 2.
The experiments on different 60 days from DOY 279 to 340 in 2019 are repeated, and then the Standard Deviations (STDs) are calculated at each epoch, and employed to describe the accuracies of the predicted clocks [25]. For example, 3h STD is calculated using 60 days clock residuals at epoch 36 (60min/5min*3 = 36), the specific formulas are as follows, where ST D j is the Standard Deviation at epoch j. ij is the difference between the predicted clock P ij and the reference clock R ij at epoch j on day i. j is the average value of ij with T days at epoch j. T is the number of days. Additionally, the ultra-rapid products are updated every 6h (3h at some analysis centers) [16], [35], and extrapolated to 24h [34], so the four typical extrapolated times, i.e., 3h, 6h, 12h and 24h, are selected. It can be seen from Table 2 that, in most cases, the iELM method performs better than GBU-P and SAM at 3h, 6h, 12h, and 24h. Since the predicted clock residuals are small within 3h, the improvements of the iELM method are accordingly not so large. However, the extrapolated clock errors are accumulated along with time, so the  iELM clocks are significantly better than GBU-P and SAM when the prediction time is larger than 18h.
Specifically, the accuracy of clock prediction is strongly related to predicted time, prediction algorithm and satellite type. GBU-P is suitable for short-term clock prediction (within 3h), which is no worse than the SAM method. But for the long-term, e.g., larger than 12h, GBU-P prediction performance is inferior to SAM. Generally, the SAM method can improve the accuracies of about 3.1%, 0.7%, 12.7%, 17.8% on average for 3h, 6h, 12h and 24h, compared with GBU-P. When the iELM method is employed, the accuracies of predicted clocks achieve 0.19ns, 0.22ns, 0.65ns, and 1.51ns for 3h, 6h, 12h, and 24h, respectively, which are much better than the results of GBU-P and SAM. The accordingly improved percentages are 33.2%, 49.0%, 27.4%, and 36.4% with respect to the SAM clocks.
Since the QZSS clock types are all of Rb atomic, the corresponding clock prediction errors can be further discussed based on satellite types, namely QZO and GEO. The clock prediction errors of QZOs (J01, J02 and J03) reach 0.14ns, 0.18ns, 0.31ns and 1.38ns at 3h, 6h, 12h and 24h, while clock prediction errors of GEO (J07) are 0.33ns, 0.36ns, 0.96ns, 1.90ns, respectively. It is obvious that the prediction performance of GEO (J07) is not so good as QZOs. The most possible reason is the inaccurate orbit/clock products of GEO lead to a poor prediction performance. It is well known that the orbit/clock products generated by analysis centers using a series of ground tracking measurements, and the observation geometry of GEO by ground tracking stations is worse than that of QZOs. Therefore, it is believed that the clock prediction accuracy will be further improved if more precise orbits/clocks are employed. Now, we examine the learning time of the ELM method, and compare it with the existing neural network methods. With hardware configurations of Intel Core i7-4500u (2.4 GHz) and 8.00 GB RAM (Random Access Memory), the average learning times of ELM, GRNN and BPNN are 0.28s, 1.68s and 1.76s respectively under the constraint that the residuals are less than 0.1ns. It can be seen that the learning times of GRNN and BPNN are much longer than twice of ELM learning time, because GRNN and BPNN both have two hidden layers (the theoretical structure for GRNN, an optimized structure for BPNN), and ELM has one hidden layer. Comparatively, SVM is not a neural network method, and its learning time (8.84s) is much longer than that of other methods.
Meanwhile, it should be noted that ultra-rapid observed orbits and clocks are calculated simultaneously [36]. With the same hardware configurations, the calculation time of ultrarapid observed orbits/clocks is about 60min (without using a cluster of processors), and the average learning time of the ELM methods is less than 0.3s. Therefore, the learning time of the ELM method is negligible compared with the calculation time of ultra-rapid observed orbits/clocks, which satisfies the requirements of the services.   same, all of the 13 stations are capable of tracking four QZSS satellites (see Fig.1).
The observations of tracking stations are collected from Crustal Dynamics Data Information System (CDDIS) servers from DOY 279 to 293 in 2019. For the purpose of a better PPP processing, the open-source software RTKLIB is further developed for this experiment with optimized parameters [37]. For example, the ISB (Inter-System Bias) estimation and observation weights between GPS and QZSS are adjusted carefully. In the multi-GNSS PPP processing, different receivers with different antennas may produce different hardware delays, but receiver manufacturers do not provide receiver hardware delay corrections [38]. Therefore, a random walk process is adopted for the ISB estimation in the GPS/QZSS PPP processing based on the RTKLIB manual [39]. Meanwhile, weights should be assigned to observations according to their relative quality. Firstly, it is known that code pseudo range and carrier phase range are two types of observations with different accuracies. Empirically, the weight of phase observation is 100 times than that of code observation [40]. Secondly, the lower elevation signals will be subject to more environmental interference. Therefore, the observation weight is also set based on the signal's elevation angle [41]. Here, the main processing strategies and parameters are listed in Table 3 for both static and kinematic PPP. The position residual and Root Mean Square (RMS) are calculated by comparing the PPP result with the station coordinate provided by IGS SINEX (Solution Independent Exchange Format) file in east, north, and up component, respectively.
In order to fully evaluate the performance of clock predictions, four strategies are designed and carried out for both the static and kinematic PPP. In static PPP, the station coordinates are estimated as constant values, which utilize the information from previous epochs in the estimation process. In kinematic PPP, the station coordinates are considered to be in continuous motion, and estimated as a white noise process without imposing any constraints between the epochs [42].

S.1: GPS PPP using GFZ predicted clock (GBU-P) S.2: GPS/QZSS PPP using GFZ predicted clock (GBU-P) S.3: GPS/QZSS PPP using SAM predicted clock (SAM) S.4: GPS/QZSS PPP using iELM predicted clock (iELM)
Generally, the 24h predicted clock product is employed in a PPP solution. It should be noted that the 24h predicted clock product is a clock series that contains 2880 clocks (corresponding to 2880 epochs) from 0 to 24h. For both static and kinematic PPP, the positioning estimation is processed with using corresponding clock at each epoch. Therefore, all clocks from 0 to 24h are used for one static or kinematic PPP solution, not using the predicted clock at 24h.
In order to ensure PPP results are affected only by clocks, the general models and parameters for the PPP processing are kept the same for the four strategies and listed in Table 3. For example, ionosphere delay is mostly eliminated with the help of ionosphere-free combination observation, while the tropospheric zenith delay is considered as a random walk parameter with optimal process noise. Further, the IGS tracking stations are set up with the strict environmental criteria, and the same observations are used for comparison, so the multi-path effects are almost the same. In addition, phase windup correction, tides correction recommended by IERS (International Earth Rotation and Reference Systems Service) including solid tide, ocean tide and pole tide, and antenna phase center correction from the newest IGS atx-file, are also modeled and corrected. Therefore, it is believed that the PPP improvements are contributed by the enhancement of the predicted clock.

2) STATIC PPP
The position RMS errors of static PPP solutions in east, north and up components for the 13 stations from DOY 279 to 293 in 2019 are presented in Fig.8, using the GBU-P and SAM methods, as well as the iELM method, respectively.
Firstly, the GBU-P orbit/clock products are employed, and the position RMS errors of static PPP solutions using GPS and GPS/QZSS observations are solved and compared. Since the QZSS constellation is a GPS complementary system in the Asia-Pacific region, the horizontal components of most stations are improved using GPS/QZSS observations. It is noted that some stations perform worse on the up component when QZSS is combined with GPS, e.g., MIZU, MRO1 and PERT. With the help of theoretical analysis, we find out that the average number of satellites used for positioning on each station increases from 11 to 14, and the related HDOP (Horizontal Dilution of Precision) and VDOP (Vertical Dilution of Precision) decrease correspondingly [43]. Therefore, the DOP analysis cannot help to interpret the positioning results at MIZU, MRO1 and PERT stations, and there may be other possible reasons. Firstly, due to some unknown interference, poor QZSS observations lead to poor results when GPS is combined with QZSS under a certain strategy. Secondly, the constraint on the tropospheric noise is too loose, which leads to systematic errors in the solution [44]. Thirdly, there may be undetected small cycle slips in the QZSS observations, and they lead to worse effects on the combined positioning [45]. All of these reasons require improved strategies, constraints and methods for specific stations, which will be studied as our future work.
Then, based on the GPS/QZSS observations, the results of GBU-P, SAM and iELM are compared. As soon as the SAM method is applied, a slight accuracy improvement is available (0.0719m) compared with GBU-P (0.0721m) in the east direction on average. And 12 stations are improved in the north direction, of which a half of stations are improved more than 50.0%. As for the up direction, there are four stations obtaining improvements, while the other nine stations get a similar accuracy with respect to GBU-P.
When the iELM method is employed, the positioning accuracy is further enhanced. Compared with SAM, there are 11 stations achieving improvements in the east component by the iELM method. The improvements reach 28.1% overall, and the best one is CUUT at the rate of 78.7%. Meanwhile, most of stations obtain improvements in the north component, at the percentage of 24.0% by iELM with respect to SAM. It should be noticed that near 40.0% stations achieve subcentimeter accuracy by iELM in the north component. While only one station (KIRI) has sub-centimeter accuracy in the north component when using GBU-P or SAM. As for the The overall positioning accuracies of daily static PPP are compared, and the results are presented in Table 4. From the details in Table 4, the accuracy of static PPP is significantly enhanced when using iELM method, especially in the north and up component. Currently, a static PPP solution usually uses clocks in one ultra-rapid file. It is noted that the sliding window mechanism is employed to generate the ultra-rapid files with updating every 6h (3h at some analysis centers) [16], [35]. Under perfect communication conditions, all the ultra-rapid files can be correctly received, and the dataset of static PPP processing can be changed as follows, the first 6h predicted clocks of four consecutive ultra-rapid files are spliced into a 24h extrapolated clock set for static PPP processing. In this way, a better positioning result can be obtained. Similarly, for some analysis centers that updates every 3h, we can further improve the positioning result by splicing the first 3h predicted clocks from eight consecutive ultra-rapid files into a 24h extrapolated clock set for static PPP processing.
Besides positioning accuracy, convergence time is also an important index. For the PPP processing, the convergence is usually defined as the position residual in each component is below a certain threshold value and maintains this accuracy during the subsequent time periods (20 epochs). According to practical engineering experiences, 2.0dm and 5.0dm are defined as the threshold values for the static and kinematic PPP, respectively [17]. Under this definition, the convergence time of static PPP is reduced by the iELM method from about 1.5h to less than 0.5h on average, except the improvement on the positioning accuracy.

3) KINEMATIC PPP
The kinematic PPP is employed to further verify the performance of the predicted clock. The details of the four strategies are described in Section IV.B Data Description. The observations of 13 stations in Fig.7 are collected from DOY 279 to 293 in 2019, which are used to carry out the kinematic PPP experiment. Compared with the IGS SINEX coordinates, the positioning residuals are solved for each station in the east, north and up component, respectively. One station is selected as an example for kinematic PPP experiment, namely MRO1, of which the position residuals  on each component are shown in Fig.9. Furthermore, the position RMS errors of kinematic PPP solutions are presented in Fig.10 for 13 stations. Fig.9 intuitively demonstrates that the iELM method provides a better accuracy and shorter convergence time. When GBU-P is used, the positioning residuals using GPS/QZSS observations have been decreased compared with GPS observations, especially in the vertical component. That is because the tracking geometry is much better with dropped PDOP (Position Dilution of Precision) when QZSS observations are involved. Based on the GPS/QZSS observations, positioning residuals using SAM are at similar magnitude with GBU-P, however, the convergence time of SAM is shorter than GBU-P in all three components. Compared with GBU-P and SAM, the iELM method performs the best in kinematic PPP for positioning accuracy and convergence time. As shown in Fig.9, the positioning residuals keep more stable with the accuracy less than 0.25m after convergence (within 2h) in horizontal components. In addition, although the up component has slightly larger absolute residuals than the horizontal ones, however, the improvement on the up component is still evident, and it is approximately at the rate of 26.0%.
The kinematic PPP positioning residuals RMSs in each component after convergence for the 13 stations are shown in Fig.10. When GBU-P is employed, the kinematic PPP positioning residuals using GPS/QZSS observations have been decreased by the percentage of 8.5%, 20.7%, and 22.4% in the east, north, and up component, compared with GPS observations. Based on GPS/QZSS observations, a similar level of positioning residuals is achieved by SAM with the accuracy changing no more than 2.0% in three components with respect to GBU-P. The further enhancement can be found by iELM as the improvement rates reach 27.5%, 23.5% and 15.1% on average, compared with SAM, in the east, north and up component, respectively. Statistically, results indicate that the kinematic PPP positioning residuals RMSs using iELM is less than 0.41m on average.
Moreover, the convergence times of kinematic PPP solutions are presented in Table 5. Generally, when the GBU-P method is employed, the convergence times of kinematic PPP solutions for most stations are more than 100min. Compared with GBU-P (GPS/QZSS), the convergence times of all stations are enhanced using SAM (GPS/QZSS) with the largest improvement of 46.2% at ALIC station. When the iELM model is used, the convergence time of kinematic PPP solutions decrease further, and most stations only need 60.0% of SAM convergence times.

V. CONCLUSION
GNSS applications are developing rapidly nowadays, especially for real-time scenarios. Although some institutes have already provided real-time GNSS orbit/clock streams, not all users can use them smoothly due to communication and equipment limitations. Therefore, the 6h updated ultra-rapid orbit/clock products are crucial to support real-time GNSS applications. However, it is a challenge task to obtain predicted clocks with a high accuracy, and the performance of ultra-rapid predicted clocks do not satisfy the requirements of some real-time applications. The reason is not only the characteristics of clock itself are very complex at the space environment, also as the clock is strongly related to the orbit solution. For the purpose of improving the quality of predicted clock, an improved model named iELM is proposed in this article.
For the purpose of obtaining key parameters of the iELM model, the periodic characteristics of QZSS clocks are studied using the FFT method. The results indicate that QZSS clocks have typical periodic terms corresponding with their orbit periods near 12h and 24h for QZOs and GEO. Further, the iELM model with satellite-specific periodic terms and coefficients are proposed. Since experiments indicate input length has great influence on the clock prediction accuracy, the best choices for the input length of each satellite are optimized. Then, the iELM predicted clocks are evaluated by comparing with GBU-P and SAM. When the iELM method is employed, the overall accuracy of predicted clocks is better than 1.0ns within 12h, which is significantly enhanced with respect to the GBU-P products. The accuracy of iELM predicted clock is considered well enough for some real-time applications.
Furthermore, the static and kinematic PPP are carried out for four strategies to investigate the behaviors of the iELM method. Considering the satellite visibility and geometry performance, totally 13 stations in the Japan and Asia-Pacific area distributed within the longitude range from 100 • E to 180 • E have selected for the PPP experiments. The only variable in different strategies is the predicted clock obtained by the GBU-P, SAM and iELM method respectively, while other parameters are fixed the same for comparison. In the static PPP, iELM enhances GBU-P as 28.3%, 57.7% and 47.4% in the east, north, and up component, based on the GPS/QZSS observations, which realizes sub-centimeter positioning in the horizontal components for near 40.0% stations and subdecimeter accuracy in the vertical component for most experimental stations. As for the kinematic PPP, the iELM method performs much better than GBU-P and SAM as a shorter convergence time and a higher positioning accuracy. The iELM method takes only a third to a half time for getting convergence. And the positioning accuracy is enhanced by iELM (GPS/QZSS) at the rate of 27.6%, 23.7% and 13.9% compared with GBU-P (GPS/QZSS) on average in the east, north and up component, respectively.