New GEONET analysis strategy at GSI: daily coordinates of over 1300 GNSS CORS in Japan throughout the last quarter century

The station coordinates derived from GNSS (Global Navigation Satellite System) with a conventional static method underpin the study of Earth and planetary science and surveying and mapping. For the continuous provision of high-quality coordinates, it is mandatory to cope with the increasing deviation from the global standard reference frame and the launch of modern GPS (Global Positioning System) satellites. To provide coordinates agreed with ITRF2014 (International Terrestrial Reference Frame 2014) at several millimeters for GEONET (GNSS Earth Observation Network System) in Japan, we developed a new analysis strategy named F5 and assessed the reprocessed station coordinates from 1996. The major updates in F5 can be summarized as incorporating global network processing and enhancements in troposphere modeling. As for the troposphere enhancements, a modern mapping function VMF1 (Vienna Mapping Functions 1) was employed and time intervals for troposphere estimates were shortened. Station coordinates in the global network showed a great agreement with ITRF2014 at several millimeters in the recent 20 years and comparable or slightly better performance with IGS (International GNSS Services) Analysis Centers. The RMS (root mean square) averaged over all GEONET stations indicated very high accuracy of 3.2 mm (horizontal) and 7.3 mm (vertical); the latter accounts for an improvement of roughly 10% from the previous strategy. Sensitivity tests about troposphere estimates revealed that the reduced RMS was completely due to the short time intervals, not the use of VMF1, which contributed to partly suppressing the spurious vertical annual deformation. These results confirm that F5 is sufficiently accurate for the requirements of individual applications and infer the capability of detecting smaller signals the previous strategy could not resolve.


Introduction
Since 1996, the Geospatial Information Authority of Japan (GSI) has operated a dense GNSS (Global Navigation Satellite System) CORS (Continuously Operating Reference Stations) network system, known as GEONET (GNSS Earth Observation Network System), which was initially established to understand the process of strain accumulation and release across Japan (e.g., Sagiya 2004;Tsuji and Hatanaka 2018). The observation data and their daily coordinates are publicly provided and widely used to maintain the basis of surveying and mapping. The current system consists of over 1300 CORS with an average spacing of 20 km (Fig. 1).
The daily coordinates of each station provide dense and continuous information about the complicated deformation field of Japan (e.g., Sagiya et al 2000;Hatanaka et al. 2003). Many challenges concerning geodesy and seismology have been resolved using the daily coordinates. These include understanding the crustal response of the loading effect by snow (Heki 2001) and heavy rainfall (Heki and Arief 2022;Zhan et al. 2021), identifying a rapidly deformed area by the secular plate motion, which is represented by Niigata-Kobe Tectonic Zone (Sagiya et al. 2000), comprehending the source process of the 2011 gigantic earthquake off the coast of Tohoku region (Ozawa et al. 2011;Suito 2017), and detecting slow slip events along subduction zones, such as Nankai trough (Kano and Kato 2020) and Ryukyu trench (Nishimura 2014). Regarding surveying and mapping, the daily coordinates have been used for semi-dynamic correction since 2010, which involves aligning positions at a current epoch obtained by GNSS to the map that is based on the survey at a reference epoch (Hiyama et al. 2010). Such geospatial information management helps support social and economic activities. Furthermore, recent applications have included use as the fundamental data for the Centimeter-Level Augmentation Service (Cabinet Office 2022), which precisely measures the positions of GNSS receivers in real time without the use of a reference station.
In other countries or regions, several efforts have been made to provide coordinates for a regionally densified CORS network. Each of these has its own inherent set of strengths and weaknesses. The location of 330 multi-GNSS stations in Europe, named EUREF (the Regional Reference Frame Sub-commission for Europe)  Bruyninx et al. 2019), is determined in the way that normal equations for the subnetworks are combined and solved with a minimum constraint condition imposed (Guidelines for the EPN Analysis Centres). In this approach, errors in the coordinates of stations considered in the minimum constraint do neither distort the network geometry nor significantly degrade the datum definition (Dach et al. 2015). However, the resultant reference frame is highly dependent on the choices of considered stations (Legrand and Bruyninx 2009). Meanwhile, some recent works highlighted the application of PPP (Precise Point Positioning) for regional network processing (e.g., Ebner and Featherstone 2008;Grinter and Janssen 2012). The PPP approach is superior to the differential one in that only a single receiver is necessary (Rizos et al. 2012), which enables to simplify the operation of routine analysis. However, long convergence times are necessary for the ambiguity float solution to ensure the sufficient accuracy individual utilization requires (Seepersad 2012;Grinter and Janssen 2012). Furthermore, it is mandatory to fully understand the implications of transforming between a global and a regional datum (Rizos et al. 2012).
Our previous analysis strategy "F3" needed to be updated owing to two reasons. The first was increasing deviation from the recent ITRF (International Terrestrial Reference Frame). F3 was based on the differential approach where a single reference station was fixed. There was a disadvantage that a common mode error appeared in the same direction and quantity at all of the GEONET stations due to the single fixing of mismodeled coordinates in the reference station. The coordinate at the reference station was determined under the regional network that covers East Asia and Oceania around Japan. The stations in the regional network have experienced either real or apparent shifts in their coordinates owing to earthquakes and equipment replacements since the end of IGS05 (Ferland 2006), i.e., the GNSS realization of ITRF2005 (Altamimi et al. 2007) where F3 was aligned. Namely, the post-seismic deformation associated with the 2011 Tohoku-oki earthquake even affected Northeastern China (Meng et al. 2019), and most of the antennas at the IGS (International GNSS Services) stations used for F3 have been replaced. As a result, the internal consistency within the network was no longer assured, and an increasing deviation from the recent ITRF had been observed particularly in the last decade.
The second reason was a forthcoming deterioration in accuracy. Bernese GNSS Software Version 5.0 (Dach et al. 2007) was used for F3 but it cannot process GPS (Global Positioning System) Block III satellites, which were thus excluded from the process of estimating the station coordinates. Currently, the absence of Block III has little influence on the derived coordinates; however, the forthcoming replacement with Block III/IIIF is expected to considerably degrade the accuracy in the future.
In this paper, we present details of our new analysis strategy of GEONET "F5" that overcome the problems found in F3 and achieve the mm-level repeatability and consistency with the ITRF. First, we provide an overview of our analysis strategies. Then, we describe the basic concept of the F5 processing scheme, which comprises two steps: reference station analysis (RSA) and GEONET station analysis (GSA). Next, we describe the specifications of each analysis. After that, we then assess the performance of F5 throughout the last quarter century in terms of station coordinates. Finally, we discuss the future prospects of our analysis strategies.

Overview of analysis strategies
In response to the deployment of GEONET in 1996, we started the routine analysis that continuously provides the coordinates of each station each day until the present. We previously updated the analysis strategy three times (Hatanaka et al. 2003;Geodetic Observation Center 2004;Nakagawa et al. 2009), and we reprocessed data from 1996 each time to generate station coordinates over the entire period in a consistent manner. In April 2021, we developed a new analysis strategy F5 and released the reprocessed products. Subsequent results are routinely appended by the automatic analysis. Anyone can access two types of products: F5 and R5. These depend on the final and rapid IGS products: satellite ephemeris and earth orientation parameter (EOP; Table 1). In addition, Q5, where IGS ultra-rapid products are applied, is available only for use within the government to rapidly grasp the crustal deformation field after an earthquake or volcanic eruption.

Basic concept of processing scheme
The series of our analysis strategies consists of two steps: the reference station analysis (RSA) and the GEONET station analysis (GSA). RSA offers the ITRF-aligned coordinates of three GEONET stations: "Tsukuba-1, " "Tsukuba-3, " and "TSKB" stations. GSA offers the In the RSA + GSA approach, errors in the coordinates of the reference station can propagate to all other GEONET stations as a common mode error. Nonetheless, we chose this approach in F5 because the two independent analyses simplify the troubleshooting of accuracy degradation due to the different noise characteristics in the coordinates. In addition, the regional network geometry is never distorted under the single fixing of the reference station.
To mitigate the common mode error and thus achieve excellent agreement with ITRF, we developed a new RSA scheme in which globally distributed IGS stations were used. Unlike the regional network of GSA, the global network of RSA has different noise characteristics due to longer baseline lengths and latitudes where stations locate. Therefore, a different set of parameters (e.g., threshold for quality checking, troposphere parameters) is adopted for each processing.

Reference station analysis (RSA)
The purpose of F5 RSA is to obtain the coordinates of the reference station for GSA that are aligned to ITRF2014. We use over 100 IGS stations that are globally distributed. For the reduced risk of abnormal solutions due to low-quality data, IGS stations with high-quality data are selected before the processing. The stations are selected according to the following steps. This procedure is applied once per year in the reprocessing term (1996 through 2020) and at recognizing a discrepancy of 1 cm (horizontal) or 2 cm (vertical) from ITRF2014 in the routine operation term (after 2021): 1. Choose IGb14 stations with an annual data acquisition rate above 70% (IGb14 is the revision of IGS14, which is the GNSS realization of ITRF2014; Rebischung and Schmid 2016; Rebischung 2020). 2. From the stations chosen in step 1, choose the IGS core sites that are used for IGb14. They have homogeneous spatial coverage on the Earth and thus are suitable for ensuring consistency with ITRF2014 (Rebischung and Schmid 2016). 3. From stations not chosen in step 2, exclude stations with a large discrepancy between the IGS daily combined solution and the IGb14 model. First, the root mean square difference (RMSD) between the combined solution and the model is calculated for each station over the year. Then, the frequency at each bin of RMSD is calculated and fitted by a gamma distribution. Finally, stations with RMSD above the 99 percentile are excluded. 4. From stations on each face of a geodesic dome as defined in Fig. 2, select one station that maximizes the distance from stations chosen by step 2, and add it to the list of stations chosen in step 2. The geodesic dome is defined in advance of the station selection and consists of 42 pentagonal and hexagonal faces. The number of faces is suitable to rapidly add the station without disrupting the spatial coverage during the iteration of steps 2-4. 5. Repeat steps 2-4 until 130 stations are chosen for the Bernese version 5.2 to be capable of processing in F5 RSA.
For the further stabilization of the reference station coordinates, additional five GEONET stations registered to IGS (AIRA, CCJ2, ISHI, STK2, TSK2, and their historical monuments) are included without coordinate constraints. Once a station nearby the reference station is removed due to missing observation, discontinuity is likely to arise in coordinates at the reference station. To prevent this, well-maintained stations are distributed nearby the reference station (Fig. 2).
The precise satellite orbit and EOP are given for the analysis. The IGS operational combined solutions are applied after February 14, 2015 whereas the repro2 solutions ) are applied before the day. The solutions before January 28, 2017 are transformed from the individual terrestrial reference frame into IGS14/IGb14 by applying the seven parameters of the Helmert transformation, which account for the origin shift, rotation, and scale shift between the pairs of reference frames, based on the method proposed by Hatanaka et al. (2003). The IGS antex (igs14_wwww. atx; Rebischung and Schmid 2016) is used for the satellite and receiver antenna phase center calibration. The a priori coordinates of each constraining station are given by the IGb14 model, which represents the station coordinates at an arbitrary epoch with the combination of the station coordinates and velocity at the reference epoch and the post-seismic deformation by major earthquakes Rebischung 2020). The reference frame is realized by a minimum constraint condition. Typically, there are two ways of implementation: nonet-translation and no-net-rotation, which respectively impose the constraints of translation and rotation so that the entire system becomes zero. We carried out RSA analysis with the different minimum constraints to examine which one is more suitable for stabilizing the station coordinates (Additional file 1: Text S1; Figure S1). Nonet-translation largely stabilized, as the result, and thus we adopted it as a minimum constraint condition. ZTD (zenith tropospheric delay) and atmospheric gradient parameters are estimated every 1 and 24 h respectively.
Global Mapping Function (Boehm et al. 2006b) was employed for a mapping function. The other models and conventions used for RSA are listed in Table 2. Complete information is given at ftp:// terras. gsi. go. jp/ data/ coord inates_ F5/ DOC/ GSI_ F5FIX. acn.

Processing flow in Bernese software
1. Preprocessing: Stations or satellites are excluded from the analysis if their quality indices retrieved by teqc (Estey and Meertens 1999), a common GNSS preprocessing tool, do not satisfy the criteria listed in Table 3. Then, cycle slips are detected and corrected by the "MAUPRP, " the sub-program of the Bernese software. 2. Baseline forming: The "SHORTEST" strategy of the Bernese software is used to create a set of the shortest baselines. Figure 2 shows an example of the station distribution and baseline network.

GEONET station analysis (GSA)
Similar to the previous strategy F3, we form three types of subnetworks called "backbone" (BB), "C, " and "A" clusters ( Fig. 3) considering the various observation periods and solve the normal equation for each subnetwork in a certain order. The number of GEONET stations has drastically increased from the end of the 1990s to the beginning of the 2000s. Some stations have been relocated during the recent 20 years and now there are around 1300 stations deployed. To cope with this and maintain the homogeneity of the solutions as much as possible over the entire period, the GEONET stations are divided into the following three types of clusters.
1. BB cluster: Approximately 20 stations are selected from the GEONET stations so that they are evenly distributed across the country. These stations should have been operated since GEONET was launched and should have particularly high-quality observing data. Then, we form a set of the shortest baseline network, namely the BB cluster (Fig. 3i), and estimate the coordinates of these stations by fixing the reference coordinate derived by RSA. The troposphere parameters are simultaneously estimated with the station coordinates. 2. C cluster: We select just under 1000 stations, which consists of the C cluster, mostly constructed in the 1990s, including stations eventually relocated. Roughly 40 stations within 1000 stations are directly linked with BB stations by the "SHORTEST" strategy of Bernese software. Then, a set of radial networks are built by connecting the rest of the 1,000 stations with selected 40 stations and BB stations (Fig. 3ii).
The station coordinates and troposphere parameters are estimated according to the fixed parameters of the BB cluster. 3. A cluster: The remaining 300 stations newly constructed after the 2000s are called the A cluster. We Fig. 3 Network combination procedure. The yellow, dark blue, and light blue circles denote the GEONET stations. The star denotes the reference station. The red lines indicate the baselines spanning each cluster, while the gray lines are those defined in the previous steps form a set of radial networks in the same way as forming the C cluster (Fig. 3iii) and estimate the station coordinates and troposphere parameters according to the fixed parameters of the center of networks.
Through these steps, we use VMF1 (Vienna Mapping Functions 1; Boehm et al. 2006a), which is newly supported by the Bernese Version 5.2. The use of NMF (Niell Mapping Function;Niell 1996), an empirical mapping function used in F3, induces a latitude-dependent spurious annual deformation across Japan (Munekane et al. 2008). This problem does not occur with the use of VMF1 owing to the direct calculation from the numerical weather prediction model with spatiotemporally high resolution (Munekane and Boehm 2010). Moreover, troposphere parameters are estimated with shorter time intervals in F5 compared with F3. In F3, the time intervals of ZTD and gradient parameters are 3 and 24 h respectively. In F5, the intervals are 1 and 3 h, respectively.
The precise satellite orbit and EOP used for GSA are the same as those used for RSA. We created the phase center variations and offsets (PCVs and PCOs) for most stations ourselves to consider the influence of the original antennas, radomes, and monuments. For typical combinations of the equipment, the absolute PCVs and PCOs aligned to ITRF2014 were calculated using data observed at a pair of stations ). The PCVs and PCOs for other combinations, which accounted for a small proportion of GEONET stations, were created via conversion of the PCVs and PCOs for F3. This procedure is as follows; the differences in the PCV and PCO between IGS05 and IGS14 were calculated for the reference combination (i.e., AOA Dorne-Margolin T antenna without any radomes) and were added to the PCV and PCO that are aligned to IGS05. The a priori coordinates for the processing of GSA were given constant values because they had little influence on the final solution under very loose constraints. The models and conventions are listed in Table 4. Complete information is available at ftp:// terras. gsi. go. jp/ data/ coord inates_ F5/ DOC/ GSI_ F5. acn.

Processing flow in Bernese software
The processing flow of the GSA procedure is as follows: 1. Preprocessing: same as RSA. 2. Baseline forming: create the BB, C, and A clusters as mentioned above. 3. Ambiguity resolution: Depending on the baseline length, different linear combinations are used. The Melbourne-Wübbena combination (Wübbena 1985;Melbourne 1985) and wide-lane combination are used for long (> 100 km) and short (< 300 km) baselines, respectively. After the wide-lane ambiguity resolution in the two ways, the narrow-lane ambiguity is resolved by using the ionosphere-free linear combination. 4. Parameter estimation: The station coordinates are estimated in conjunction with the troposphere parameters in the order of the BB, C, and A clusters.

Frame realization
The coordinates of the reference station must be consistent with ITRF2014 and stable because the subsequent GSA is processed by fixing the reference station. To Although Tsukuba-1 is usually selected as the reference station, we used TSKB, which is registered to IGS and thus exists in most ACs' and IGS combined solutions. IGS operational solution was adopted from February 14, 2015 whereas the repro2 solution was used before the day. The specifications of the repro2 and operational solutions are summarized in Table 5. Station coordinates aligned to the original reference frame were converted into IGS14/IGb14 by application of the Helmert transformation between corresponding ITRF to ITRF2014. The station coordinates of F3 were also employed for comparison after the coordinate transformation to ITRF2014. For quantitative assessments, statistics of mean error (ME), standard deviation (STD), and root mean square error (RMSE) were also evaluated (Additional file 1: Text S2 for the definition). The time series of F5 showed good consistency with the IGS combined solution at several millimeters in the 2000s or later (Figs. 4, 5). EMR, GRG, JPL, and MIT were more biased and scattered than F5. The large discrepancy of several centimeters found in the previous F3 solution was no longer observed in the new F5 solution (Fig. 5). Furthermore, the degree and number of outliers, which were significant in F3 for the recent several years, were largely decreased. Both ME and STD demonstrated a mm-level consistency and repeatability in F5 after the 2000s (Table 6a and b). They were slightly better than MIT and came close to COD. The RMSE for EW, NS, and UD components was 3.5, 4.0, and 7.2 mm in 2007-2009 and was 2.3, 2.0, and 4.6 mm in 2017-2019, respectively (Table 6c). These statistics were also comparable to each AC, except for ULR.
These results suggest that the coordinates estimated by F5 RSA are successfully aligned to ITRF2014 with mm-level accuracy. However, COD, ESA, GFZ, and GTZ showed excellent consistency and repeatability, which were superior to F5. Among them, COD and ESA incorporated GLONASS (Global Navigation Satellite System) in addition to GPS, which likely contributed to the improved accuracy. Moreover, at the end of the repro2 period, COD, GFZ, and GTZ used roughly twice as many stations as F5 where up to 130 stations were integrated. This could have enriched the robustness against contamination of low quality data; thus, contributing to excellent consistency and repeatability. In the 1990s, however, F5 was more biased and scattered than after the 2000s. This is probably caused by selective availability, which is an intentional degradation of GPS signals applied before May 2000, and inhomogeneous station distribution in the processing network.

Noise characteristics
The spectral features were obtained to understand the noise characteristics of RSA. Figure 6 plotted the power spectrum in 365 days shifting the time window of half a year from 1996 to 2020 (See Additional file 1: Text S2 for the detailed procedure). Each power spectrum was proportional to the inverse of frequency except for frequency-independent noise found at the highest frequency bands. These characteristics are consistent with past studies (Mao et al. 1999;Williams et al. 2004), which revealed that the noise hidden in station coordinates time series can be modeled by a combination of flicker noise and white noise when stations are globally distributed. The power spectrum after 2003 showed a relatively modest noise level than that before 2002; the noise level reduced to a quarter and a half for the horizontal and vertical components respectively. The possible reasons are the increased number of stations over 2003 and the discontinuation of selective availability. Figure 7 showed the F5 and ACs power spectrum stacked from 2003 through 2020, where the spectra in 365 days showed a modest noise level (See Additional file 1: Text S2 for the detailed procedure). F5 showed a comparable or slightly better noise level among EMR, JPL, and MIT, which only use GPS. COD and ESA, which incorporate GLONASS (Global Navigation Satellite System) in addition to GPS, exhibited spectra with lower noise levels than those of the F5 and GPS-dependent ACs.
F5 showed a spectral peak around 25.0-26.9 cpy (13.6-14.6 day) as well as most ACs spectrum (blue hatched   Fig. 7). Because the Mf and Mf′ tides (i.e., a pair of long-term tides with periods of 13.66 and 13.63 days respectively) have a particularly substantial tide potential, their mismodeling directly imposes a spectral peak around the fortnightly band . In addition to the direct error, errors in sub-daily EOP tidal variations with periods close to 12 and 24 h are absorbed into the resonant GPS orbit and daily EOP parameter estimates, which resulted in aliased signals around the fortnightly band as well as GPS draconitic period for 24 h sampled processing (Penna and Stewart 2003;Griffiths and Ray 2013 (Penna and Stewart 2003). As an exception, the JPL spectra did not show any peaks around this range, which can probably be attributed to the 30 h sampled processing ). COD and ESA exhibited a peak at 44.9-46.9 cpy (7.8-8.1 days) followed by weak signals at the second and third harmonics [90.8-92.8 cpy (3.9-4.0 days) and 134.7-138.7 cpy (2.6-2.7 days)] while those were absent in F5 and other ACs (pink hatched bands in Fig. 7). Because the GLONASS tracking stations have a strongly non-uniform distribution, the orbit error with 8-day period (i.e., the ground repeat period) appears to have caused the  Fig. 4, but at TSKB. The daily coordinates of F3 and their moving average were also plotted periodic error in the station coordinates . The peak was also absent from the GFZ spectra because GLONASS data were not used for the repro2, which accounted for most of the stacking period.

Performance of GSA
The accuracy in all GEONET stations in 2019 was assessed by RMS of the baseline vectors after removing the linear signal (See Additional file 1: Text S3 for the definition). F5 showed latitude-dependent RMS of 1.5-2.5 mm (horizontal) and 6-8 mm (vertical), which increased toward the south (Fig. 8a). The RMS averaged over the entire nation indicated 2.3 mm (EW) and 2.3 mm (NS), which correspond to 3.2 mm in the horizontal, and 7.3 mm in the vertical component. A similar spatial pattern was found in F3 (Fig. 8b), but F5 showed a clear nationwide improvement in the UD component (Fig. 8c). The averaged RMS in the UD component accounts for roughly 10% improvement (F5: 7.3 mm, F3: 7.9 mm), despite the deterioration in the southwestern islands of Japan.
The latitude-dependent distribution of RMS would be relevant to the zenith wet delay, which is the wet part of ZTD and is responsible for large spatiotemporal variation (Leick et al. 2015). In Japan, the spatiotemporal distribution of water vapor is highly variable in summer (Naito et al. 1998;Iwabuchi et al. 2000). As a matter of fact, the time series of the UD component at "Muroto-3" with respect to "Misumi" showed a large variance in summer (Fig. 9). Meanwhile, the variance in summer was largely mitigated in F5. In the next subsection, we will further discuss how and to what extent the different troposphere estimates impact the spatiotemporal variation of the positioning accuracy.

Sensitivity of troposphere estimates
As for troposphere estimates of F5 GSA, we newly employed VMF1 and shortened time intervals. To distinguish each contribution, we conducted two extra analyses: the same setting as F3 but VMF1 was applied (F3VMF1) and the same as F3VMF1 but the time intervals for F5 were applied (F3VMF1_SHORT). Figure 10 plotted the RMS in the UD component by the individual solution and their difference (F3VMF1 minus F3 and F3VMF1_SHORT minus F3VMF1). There was little difference in RMS between F3 and F3VMF1 (Fig. 10d). In contrast, the difference between F3VMF1 and F3VMF1_SHORT (Fig. 10e) showed a highly correlated spatial pattern with the difference between F5 and F3 (Fig. 8c). These results suggest that the reduced RMS in F5 was completely due to the shortening of time intervals, not the use of VMF1.
We obtained the implication of high-frequent troposphere estimates based on the precipitable water vapor (PWV) retrieved from ZTD. We used the method proposed by Ohtani and Naito (2000) for PWV retrieval. For its verification, surface specific humidity (SSH) provided by the Mesoscale Model in Japan Meteorological Agency (JMA; Ishida et al. 2022) was used. Figure 11 indicated the temporal evolution of PWV and SSH at the "Shiojiri" station on July 24, 2019. In F3VMF1_SHORT, PWV was maximum at 10 UTC (19JST) and then fell at the subsequent hour, which corresponds to the temporal changes reproduced by SSH (Fig. 11). This can be expressed as the thermal-induced local circulation; sea breeze driven by the differing heat capacities on land and sea helped concentrate water vapor in the inland area, and then water vapor was released through intensive rainfall, which was indeed observed in the neighboring AMeDAS (Automated Meteorological Data Acquisition System) station after PWV reached its peak. By contrast, such a rapid time evolution was not tracked by F3VMF1; a broad and moderate change was observed over the day. Based on these results, it is clear that the troposphere estimates with the short time intervals improve the accuracy of station coordinates via the realistic tracking of the temporal evolution of water vapor.
In the southwestern islands of Japan, significant outliers were rarely observed in F3VMF1_SHORT. On the day with the outlier, the temporal evolution of ZTD in F3VMF1_SHORT showed an instantaneous rise across the archipelago around 15UTC while a gradual transition was observed in F3VMF1 (Additional file 1: Figure Table 6 Station coordinate statistics at TSKB The colored bar corresponding to each statistics value is also drawn for ease of understanding S2). Based on the flat time series in SSH (Additional file 1: Figure S3), this sharp rise would not be relevant to a meteorological signal, but instability of the parameter estimation. Because stations in the BB cluster were combined at first and solved by fixing Tsukuba-1, spatially heterogeneous noises were radially redistributed across Japan. Therefore, the noise in the cluster was imposed on the furthest station, which could invoke the outlier in the southwestern islands.
We introduced an index of coordinates difference between August and February (CDAF), which is the coordinate averaged over August with respect to that in February. This index is capable of capturing the spurious annual deformation in NMF because it contains phase lags of 220 degrees in an annual cosine curve (Munekane et al. 2008), which corresponds to spurious subsidence in February and uplift in August. Figure 12 plotted of CDAF in the UD component by the individual solution and their difference (F3VMF1 minus F3 and F3VMF1_SHORT minus F3VMF1). The CDAF varied from station to station (Fig. 12a-c) because it involved not only spurious signals but also station-dependent crustal signals including co-/post-seismic deformation, SSE, and snow loading modulation. The ∆CDAF between F3 and F3VMF1 increased/decreased toward southwest/northeast from the reference station and was spatially correlated (Fig. 12d). Meanwhile, the ∆CDAF between F3VMF1 and F3VMF1_SHORT showed positive signatures in the entire domain and was spatially uncorrelated (Fig. 12e).
The past studies based on PPP analysis revealed that the spurious annual deformation reaches up to 3 mm half amplitude in the middle to high latitudes (Munekane et al. 2008), and VMF1 contributes to mitigating the spurious signals (Munekane and Boehm 2010). The decrease in ∆CDAF in the northern part of Japan was consistent with those studies but much less than 3 mm (6 mm in peak to peak). Instead, the ∆CDAF showed positive in the southwestern part of Japan (Fig. 12d). The time series at Gusukube certainly illustrated that the use of VMF1 mitigated the spurious annual deformation in F3 (Fig. 13). This discrepancy could be attributed to the nature of the network combination strategy in GSA. Under the single fixing of Tsukuba-1, spatially heterogeneous noises were radially redistributed within the whole network. As the result, the different noise characteristics between southwestern and northeastern Japan could be homogenized, which in turn mitigated the spurious movements in southwestern Japan. To confirm this hypothesis, a quantitative assessment of the error propagation is required.

Future prospects
For further improvement on the accuracy, several issues are still open to discussion. Using multi-GNSS is one of them. As we have shown previously, COD and ESA, which incorporated GLONASS in addition to GPS, highlighted the smaller noise level than other GPS-dependent ACs and F5. This implies that the increasing number of satellites is capable of satisfying growing demands on positioning accuracy. In fact, the increased number of ACs has integrated GLONASS and Galileo as well as GPS and has contributed to constructing ITRF2020, the successor of ITRF2014 and the latest ITRF (Rebischung 2021). Thus, processing with multi-GNSS constellations is important to improve the alignment to the ITRF successively updated in the future. The state-of-the-art modeling technique should be implemented in response Fig. 7 Power spectra for the station coordinates derived from the individual processing. The power spectra in 365 days were stacked over the stations and period from 2003 to 2021 shifting an interval of half a year. Except for GRG and GTZ, the repro2 solution is used up to 2014. ULR was eliminated because the noise level was too high to draw a consistent scale. The blue hatched area indicates the fortnightly band. The pink hatched area indicates the period of 8 days and its second and third harmonics to the update in IERS Conventions. For example, the IERS Conventions recommend applying higher-order ionosphere correction, modern ocean tidal loading models, and secular pole tide. Although the higher-order ionosphere correction was not taken into account in F5 RSA, it has diurnal, semiannual, and decadal signatures of several millimeters (Kedar et al. 2003). The implementation of minimum constraints condition is worth reconsidering as well. In COD, the geocenter motion parameter is introduced as an additional parameter and absorbs the detachment that ITRF in seasonal or shorter time scales is no longer the center of mass but the center of figure (Dong et al. 2003;Zajdel et al. 2019). However, ACs except for COD do not estimate the geocenter motion in repro2 and operational solutions. Further assessments are required for what type of parameters should be estimated in future strategies.
As a vital infrastructure for social and science communities, the daily coordinates in GEONET should have robustness as well as good accuracy. For this purpose, the use of QZSS (Quasi-Zenith Satellite System) is beneficial. As of 2022, four QZSS satellites are in commission and three more are planned to become available at the end of the 2023 fiscal year, which provides a preferable satellite constellation in the Asia-Oceania region. GEONET particularly embraces the benefit of the constellations due to

Conclusion
The daily station coordinates of over 1300 GEONET stations in Japan are routinely provided based on an analysis strategy. With the aim of alignment to ITRF2014 at several millimeters, we developed a new GNSS analysis strategy, named F5. Similar to the previous strategy F3, the RSA + GSA approach was employed. Namely, RSA provides the coordinate aligned to ITRF2014 at one of the GEONET stations, and GSA yields coordinates at all GEONET stations fixing the RSA-derived coordinate as a reference. To achieve excellent agreement with ITRF2014, we processed data of globally distributed IGS stations in RSA. The troposphere estimates were enhanced in terms of the application of the modern mapping function VMF1 and the shortening of time intervals.
Based on F5, we reprocessed the data from 1996 and assessed the obtained station coordinates. After the 2000s, the station coordinates in the global network were consistent with the IGS combined solution at several millimeters, which was comparable with IGS ACs' performance. In TSKB, which is the candidate reference station in GSA, the RMSE with respect to IGS combined solution was 3.5 mm (EW), 4.0 mm (NS), and 7.2 mm (UD) in 2007-2009 and was 2.3 mm, 2.0 mm, and 4.6 mm in 2017-2019, respectively. In the 1990s, however, the station coordinates were biased and scattered. This was probably due to selective availability and the poor spatial coverage of the IGS-registered stations. The power spectra stacked over the entire network showed a noise level as small as other GPS-dependent ACs and depicted a common peak at the fortnightly band.
The coordinates at all GEONET stations indicated the RMS of 3.2 mm and 7.3 mm for the horizontal and vertical components respectively. In particular, the vertical component was improved by roughly 10% from F3. The sensitivity tests about the troposphere estimates revealed that this improvement is completely owing to the shortening in time intervals not the application of VMF1. The troposphere estimates with the short time intervals ensured us to successfully model the change in water vapor accompanied by a local atmospheric circulation system. The sensitivity tests also showed that the use of VMF1 mitigated the spurious annual vertical deformation to some extent. However, the spatial distribution of the mitigation was different compared with the past studies probably due to the nature of the single fixing.
These results confirmed that F5 has sufficient accuracy to support the study of Earth and planetary science and maintain the national geodetic datum. In addition, the performance of F5 suggested the capability of detecting smaller signals that were not resolved by the previous strategy. As a vital geodetic infrastructure, further efforts are needed to provide excellent accuracy and robust positioning.
Additional file 1: supplementary texts and figures.