Abstract
In recent years, due to the development of precise point positioning (PPP), the inconsistency of map and satellite positioning results due to crustal deformation has come to the surface. In order to solve this problem, a parameter for correcting the crustal deformation is necessary, and if it is provided as a service, the validity period and the accuracy guarantee are also necessary. In this research, we considered the validity period of parameters and the estimation method of the coordinates used for calculating the amount of crustal deformation and actually constructed the crustal movement correction parameter generation system and accuracy evaluation system. The parameters generated by this method are called semi-dynamic reduction (SD/R) parameters. As a result of evaluating the accuracy of SD/R parameters, we achieved a target accuracy of 4 cm at 98.3% verification points in the horizontal direction. Therefore, the effectiveness of the proposed method was shown.
1 INTRODUCTION
Conventionally, when performing highly accurate positioning, a relative positioning representative of RTK (real-time kinematic)-GPS has been used. However, in recent years, methods called PPP (precise point positioning) and methods called PPP-RTK that further solve integer ambiguities have been developed. In PPP, a bias term called carrier phase ambiguity is included in the carrier phase observation amount used for positioning calculation. This bias term includes the uncalibrated phase delays (UPD) present in the satellite and receiver carrier phase, and by removing this UPD, the ambiguity can be obtained as an integer value. This technology is called PPP ambiguity resolution (PPP-AR). Two methods are mainly available for this1,2: integer recovery clock (IRC) method and fractional cycle bias (FCB) method. These PPP positioning methods are expected to be utilized in the future in fields such as agriculture and automatic driving as low-cost methods that do not require a reference station.
In Japan, service of the quasi-zenith satellite system “Michibiki” started on November 1, 2018 (JST), and centimeter-class positioning augmentation service (CLAS, MADOCA), which enables high-accuracy positioning, is provided. These services target PPP. Regarding positioning accuracy, for example, when CLAS is used, the horizontal direction RMS is 6 cm in static observation and 12 cm in kinematic observation.3
However, at the same time as this development, inconsistency between satellite positioning and maps have surfaced.4 In other words, since the independent positioning result is a four-dimensional observation value including time information, it is affected by crustal deformation. Therefore, due to the precision of positioning results, deviation from geospatial information became noticeable. Actually, the crustal deformation in Japan’s year may be larger than the positioning accuracy of the CLAS method in some cases. In Japan, the crustal deformation is monitored by the continuously operating reference station (CORS), and on the site5 of the Geospatial Information Authority of Japan (GSI), the amount of crustal deformation in various parts of Japan has been disclosed recently. In this site, you can get the 1-year fluctuation vector of crustal deformation in Japan as a text file or KML file. Values that can be acquired include a reference point code, a reference point name, reference point coordinates, and a fluctuation vector (east-west direction, north-south direction, and height direction). Figure 1 shows a histogram of the length of the fluctuation vector and the number of reference points using the fluctuation vector of the Tohoku region obtained from this site. Here, we use data from 2018/01/04 – 2018/01/18 to 2019/01/05 – 2019/01/19. In particular, in the Tohoku region of 2018, crustal movement amount of 6 cm or more per year is observed at about 77 reference points.
Histogram of the length of the fluctuation vector and the number of reference points
Therefore, although the positioning result was acquired with high accuracy, when plotted on a map due to the crustal deformation, the plotted point shifts. In particular, since satellite positioning results are often used together with geospatial information, it is natural to correct the amount of crustal deformation.
Regarding the correction of crustal deformation, in Japan, the GSI publishes correction parameters for crustal movement on April 1 (JST) every year. Corrections using these parameters are called semi-dynamic corrections and are used for correction of reference point survey to eliminate distortion due to crustal deformation.6 This correction is an important component of Japan’s semi-dynamic geodetic system and interpolates the coordinate difference between the current period and the original period by the Kriging method.7 In this paper, correction parameters used for semi-dynamic correction are called semi-dynamic (SD) parameters. However, since the SD parameter is generated based on the coordinate value of January 1 every year and is updated only once a year, crustal movement amount of up to 15 months may not be corrected in some cases. Therefore, when providing it as a service, it is necessary to shorten the update interval of parameters, to set a valid period, and to ensure accuracy and quality guarantee within their conditions.
In this study, we clarified the valid period of the crustal movement correction parameter, examined the method of estimating the coordinate value to be used for generating the parameter, and proposed the index for the accuracy guarantee. Based on this result, a crustal movement correction parameter generation system and an accuracy evaluation system were constructed, and the generated semi-dynamic reduction (SD/R) parameter and SD parameter were compared and evaluated.
2 CONSIDERATIONS
2.1 Estimation of the coordinate value of the reference station
The amount of crustal deformation is generally calculated by taking the difference between the coordinate values before and after the movement. To calculate the amount of crustal deformation across Japan, highly precise crustal deformation can be obtained by using the F3 solution of the CORS provided by the GSI. “F” in the F3 solution represents “Final.” “Final” means that the ephemeris used for analysis is IGS’s Final product.8 The F3 solution is obtained from the result of precision baseline analysis using the analysis software “Bernese GNSS Software,” from Bern University.9 However, the coordinate values of the F3 solution are compliant with ITRF2005. The fact that the F3 solution uses the old ITRF2005 is also mentioned by the GSI itself, which is analyzing and releasing, and it is peeling about 3 cm from ITRF2014. On the other hand, there is a movement to realize a new analysis method including improvement of this problem. This new method is supposed to make the coordinate values of the F3 solution conform to ITRF2014.10
Furthermore, in order to determine the valid period of the correction parameter, it is necessary to estimate the future coordinate value. However, since the F3 solution is analyzed using IGS’s Final products updated every 2 weeks, the coordinate values 2 weeks ago are analyzed and made public.
Therefore, it is necessary to estimate the F3 solution after about 3 weeks in order to generate the parameter of validity period 2 weeks. Figure 2 shows the relationship described above in chronological order.
Relationship between the interval in which F3 coordinates exist and the date on which coordinate values are to be estimated
There are many methods for estimating the coordinate values, such as linear regression analysis and time series analysis. For example, when the GSI creates SD parameters, it estimates the coordinate value on January 1 from the observation value of the previous year by extrapolation, using an estimation model M based on a straight line and a quadratic term of the Fourier series.11
The same model may be used, but in this case, we will estimate the coordinate value after about 3 weeks, so we investigated whether the above formula is really suitable for estimating the coordinate value. As a comparison target, the following linear regression model was used.
In the comparative evaluation of the estimation model, evaluation was performed with the mean value, the standard deviation, and the ratio of the reference stations with small error when the actual F3 solution was taken as the true value, with respect to the error of the latitude, the longitude, and the ellipsoidal height. The number of data points used for the estimation was six patterns of 100, 150, 200, 250, 300, and 365. Increasing the amount of data increases the possibility of including unsteady crustal deformation such as earthquakes, so the upper limit of the number of data is set to 365. The estimated date is three patterns of January 15, 2018, July 15, 2018, and December 15, 2018. In addition, the reference station where the F3 solution has few published coordinates and greatly affects the error is excluded (KYOUGOKU-A, SHIRIUCHI-A). The results of comparative evaluation of the estimation model are shown in Figures 3 through 11. Figure 4, 7, and 10 show the results of calculating which error is smaller between the GSI model and the linear regression model at each point of CORS. (Here, the error is the difference between the actual value and the estimated value.) In other words, a model with a higher ratio is a more effective model.
Average of error between F3 solution actual value and estimated value on January 15, 2018
Inverted ratio of error between F3 solution actual value and estimated value on January 15, 2018
Standard deviation of error between F3 solution actual value and estimated value on January 15, 2018
Average of error between F3 solution actual value and estimated value on July 1, 2018
Inverted ratio of error between F3 solution actual value and estimated value on July 1, 2018
Standard deviation of error between F3 solution actual value and estimated value on July 1, 2018
Average of error between F3 solution actual value and estimated value on December 15, 2018
Inverted ratio of error between F3 solution actual value and estimated value on December 15, 2018
Standard deviation of error between F3 solution actual value and estimated value on December 15, 2018
For the linear regression model, the error of the estimated coordinates is stable regardless of the amount of data. However, increasing the data does not significantly improve the error. On the other hand, for the GSI model, the error tends to decrease as data increases. In particular, using the data for 1 year, the GSI model has smaller errors in most cases. From the figure, the trend can be read, but it is difficult to identify the fine permutations. Therefore, scores are given in order of decreasing error index (absolute average of error + standard deviation) in Table 1.
Score of total errors
We ranked in descending order of score, with the result of capturing the features of Figures 3 to 10. The highest score was the linear regression model where the number of used data was 200. Therefore, the crustal movement correction parameter generation system was generated based on this model. The method in the generation system is presented in the next section.
2.2 Jump of coordinates of F3 solution
Since the F3 solution contains noise, estimation accuracy improvement can be expected by performing smoothing by moving average or filtering as preprocessing. For the analysis of the F3 solution, the coordinate value of “TSUKUBA 1” is used as a fixed point. In other words, the national GNSS network in Japan is based on only one station “TSUKUBA 1.” Therefore, if the coordinate value of “TSUKUBA 1” jumps, it will affect the F3 solution analyzed from the entire nationwide GNSS network. Moreover, the coordinate value of “TSUKUBA 1” is calculated in the network including IGS stations around Japan. However, when other IGS observation stations have confirmed that the quality of the observation data has deteriorated, the coordinate values of “TSUKUBA 1” are calculated excluding those stations. This causes a jump in the coordinate value of “TSUKUBA 1.” The GSI releases information about coordinate jumps a few days to a few weeks after the release of the F3 solution, but the F3 solution will not be reanalyzed.12
The coordinate values in the north-south direction and the east-west direction of the reference station “TSUKUBA 1” in 2018 are shown in the Figures 12 and 13.
Jump of north-south direction of F3 solution on TSUKUBA 1
Jump of east-west direction of F3 solution on TSUKUBA 1
The part marked with a dashed circle is marked based on the coordinate jump information published by the GSI. However, since the jump is obviously larger than the noise, it is expected that it will be better to remove it before smoothing.
3 METHODS
3.1 Estimation of theoretical variogram model
For the calculation of the crustal movement correction parameter, we used the ordinary Kriging method used in geostatistics.,13,14 In this method, weights are estimated statistically for surrounding observations and an estimate at any point is calculated by calculating the sum of the weights. SD parameters are estimated by the same method. In order to use the Kriging method, it is necessary to estimate a theoretical variogram which is statistically a spatial correlation. For the theoretical variogram estimation model, spherical model C3(r) is used, and parameters a and b of variogram function 
are estimated. Let r be the distance between the reference stations. In the method by Ichikawa and Hosoi,4 the parameter b of the scaling was omitted, but in this paper, we improved the scaling so that estimation can be done at the same time. In general, the spherical model is defined as C3(r) = 
,so if we let 
, we get the following.
Specifically, we derive the solution by estimating the parameters of this model by the method of least squares. The data to be fitted are empirical variograms calculated from the movement of the reference stations. The empirical variogram is classified by the distance between the reference stations, and the average value of the dissimilarity corresponds to the frequency. A representative value (basically an intermediate value) of each class is xi, and an average value of corresponding dissimilarities is yi. Calculate the following quantity in advance.
Here, the subscript of the sum of the above equation moves the index of the bin of the empirical variogram. Assuming η = (αγ − β2), the product of parameters b and c is given by
On the other hand, since the parameter c can be calculated as
the parameter b can also be calculated. Using the above equations, the theoretical variogram was calculated. The data uses the F3 solution on January 1, 2017, and March 31, 2018, and the width of the bin is 30 000 m. We calculated the spherical model in two steps. This is because even if all the data are used, the estimation result is affected by a pair of points that are too far apart. Therefore, in Step 1, a spherical model was estimated using all point pairs, and we calculated the range a, which is a criterion of the threshold between the points. In Step 2, we set the value close to range a as a threshold, exclude a pair of points with a distance greater than that value, and estimated the spherical model again. The actual calculation results were 897 km for latitude, 752 km for longitude, and 749 km for height; hence, we set the threshold for point distance to 700 km. Results are shown in Figure 14.
Comparison of empirical variogram and theoretical variogram
With a near field of about 500 km, it seems to be working well, including scaling, but there is room for improvement. Regarding latitude and ellipsoidal height, dissimilarities up to 165 km are large, and it is necessary to select another model to fit the empirical variogram.
If we can estimate the theoretical variogram model, we can estimate the amount of crustal movement at any point by calculating the ordinary Kriging method.
3.2 Generation of SD/R parameter
Simply calculating the amount of crustal movement by the Kriging method cannot accurately map PPP results to the map. In the first place, since the Kriging method is a method to estimate steady crustal movement, crustal movement due to a large earthquake cannot be accurately estimated in some cases. Since the SD parameter is a parameter including correction of such crustal movement, we propose a method of estimating and correcting the amount of crustal movement from January 1 of that year once it is corrected by SD parameters. Specifically, the correction amount is given to lattice points covering the whole of Japan called the regional mesh that the statistics bureau has opened for the SD parameter. Therefore, similarly to the SD parameters, the amount of crustal movement can also be obtained on the lattice points of the regional mesh and then combined with the SD parameters. Parameters generated by this method are called SD/R parameters.
4 EVALUATIONS
4.1 Comparison of SD/R parameter and SD parameter
In this section, comparison and evaluation of SD/R parameters and SD parameters are carried out, and the results are shown. In the evaluation method, the F3 solution is corrected with the crustal movement correction parameter to be evaluated. Next, the difference between the corrected F3 solution and the coordinate value of the Japanese Geodetic Datum 2011 (JGD2011) was evaluated as an evaluation index (or an error index). The JGD2011 is a Japanese geodetic datum and is defined by the Survey Act of Japan and the Survey Act Enforcement Order of Japan. In the JGD2011, surveying standards are defined by the Survey Act of Japan, and the origin of Japanese longitude and latitude is determined by very-long-baseline interferometry (VLBI) according to the enforcement order. The enforcement order also stipulates the Japanese standard origin and the Earth ellipsoid. Japan adopted the world geodetic datum from the Japanese Geodetic Datum 2000 (JGD2000), the previous geodetic datum.15-18 The JGD2011 was inherited from the concept of JGD2000 and was constructed by changing the geodetic origin defined in the Survey Act Enforcement Ordinance of Japan after the 2011 earthquake off the Pacific coast of Tohoku.
So that accuracy degradation due to the passage of SD parameters over time can be observed, we chose 3 days: April 8, 2017 (UTC), October 1, 2017 (UTC), and March 24, 2018 (UTC).
The F3 solution used also includes the reference station installed on the island. Regarding the generation of SD/R parameters, we estimate the amount of crustal movement using 200 points from 7 days ago. Evaluation results are shown in Figures 15 and 16.
Comparison of horizontal error between SD parameter and SD/R parameter
Comparison of height error between SD parameter and SD/R parameter
Figure 15 shows the histogram and the cumulative line graph on the horizontal direction error. There is no big difference on April 8, but we can confirm that the SD parameter has larger error with time. On the other hand, the SD/R parameter shows no large increase in the error amount even after a lapse of time, and the error is less than the target accuracy of 4 cm or less at most reference stations. The target of horizontal error of 4 cm or less is because we wanted to make it less than 10 cm together with CLAS stationary horizontal accuracy of 6 cm; there is no particular scientific basis. The same is true for the target accuracy of 8 cm in the height direction.
Figure 16 shows the histogram and the cumulative line graph on the height direction error. As for the result on October 1, it can be confirmed that the error is smaller in the SD/R parameter, but there is no particularly large difference.
Next, the ratio of the reference station at which the horizontal direction error becomes 4 cm or less is taken in Table 2, and the ratio of the reference station at which height error in the height direction is 8 cm or less is shown in Table 3. Although the error in the horizontal direction of the SD parameter increased from Table 2, the error of the SD/R parameter hardly changed, and the ratio of the reference station at 4 cm or less for all days was 98.3%. This shows the effectiveness of SD/R parameters in the horizontal direction. Regarding the error in the height direction, the ratio of the reference station at 8 cm or less was in any case more than 99%.
Percentage of reference stations where horizontal error is less than 4 cm by correction of SD/R parameter or SD parameter
Percentage of reference stations where height error is less than 8 cm by correction of SD/R parameter or SD parameter
4.2 Error index of SD/R parameter
When providing parameters as a service, it is sometimes necessary to create data according to a provided protocol. Given the usage scene of crustal movement correction, periodic network distribution can be mentioned. Therefore, it is better for the error index to be smaller in size. Furthermore, it is necessary to evaluate how much error exists for any service area, and it must be easy for users to use. Therefore, we think that it is preferable to calculate the error amount for each mesh of the regional mesh. Hence, we propose a method to estimate the amount of error on the regional mesh from the error amount at the reference station and an error index by the ordinary Kriging method of the point to block type. This method was realized as an accuracy evaluation system. We selected the data on March 31, 2018, which will actually increase the error and carried out the accuracy evaluation. The results are shown in Figure 17.
Error map (horizontal error)
5 CONCLUSIONS
When generating crustal movement correction parameters using the F3 solution provided by the GSI, we mentioned that the valid period is about 2 weeks at the shortest and it is necessary to estimate the coordinate value after 3 weeks for calculating the crustal movement amount. The coordinate values are shown to be effective by estimating with the linear regression model using the F3 solution for 200 days. In the estimation of the amount of crustal movement, we show the solution in the parameter estimation of the spherical model which is one of the theoretical variogram models and showed that the solution is effective. However, in the parameter generation method shown, the F3 solution is estimated with a linear model, but this model is not optimal for postseismic deformation after a big earthquake. For example, it is known that the postseismic deformation of a big thrust earthquake such as the 2011 earthquake off the Pacific coast of Tohoku is not linear but has a logarithmic + exponential trend.19 A method of generating SD/R parameters is shown, a crustal movement correction parameter generation system is actually constructed, and the generated SD/R parameters and SD parameters are compared and evaluated. Then, although the effect was about the same with respect to the error in the height direction, both correction parameters were highly effective. Regarding the error in the horizontal direction, the SD parameter deteriorates in accuracy, whereas the SD/R parameter continues to maintain high accuracy. In addition, the reference station within the error of 4 cm showed a high numerical value of 98.3% or more under any condition. Therefore, the validity of the SD/R parameter was fully shown. Finally, we proposed a lightweight error index and its method that can be evaluated with accuracy by combining the block to point type ordinary Kriging method and the regional mesh.
Since it is known that the weight of the Kriging method takes a negative value as a future task, we investigate the influence of that phenomenon on estimation. A number of papers have already been published regarding solutions to this problem. For example, Barnes and Johnson20 utilize the quadratic programming method to solve the condition that the weight is nonnegative and skillfully utilizes the block matrix for calculation. Szidarovszky et al21 deforms expressions without using the quadratic programming method and uses computational cost by using a search tree. If there is an influence on the accuracy of the crustal movement correction parameter, we would like to consider countermeasures.
HOW TO CITE THIS ARTICLE
Ichikawa Y, Hosoi M. Generation of crustal movement correction parameter considering valid period and accuracy guarantee. NAVIGATION-US. 2020;67:307-317. https://doi.org/10.1002/navi.358
ACKNOWLEDGMENTS
The authors are grateful to Dr Katsumi Nakane of AISAN TECHNOLOGY CO., LTD., who gave us the idea of SD/R. Thanks also to AISAN TECHNOLOGY CO., LTD., for giving us opportunities for research and development.
- Received May 12, 2019.
- Revision received November 15, 2019.
- Accepted January 17, 2020.
- © 2020 Institute of Navigation
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.