Abstract
GPS (global position system) monitoring has great theoretical and practical significance for understanding regional crustal deformation and disaster prevention and mitigation. Based on the GPS system differential model, this paper builds a seismic real-time data analysis model and discusses and analyzes the GNSS (global navigation satellite system) single-system single-difference model, single-system double-difference model, and dual-system double-difference model. When the model is processing the data of the GNSS dual-satellite positioning system, the proposed short-baseline ambiguity solution is studied to a certain extent and solves the problems that other common methods can only detect large cycle slips or require dual-frequency data. During the simulation process, through the study of the GPS coordinate time series, the stability of the site, the periodic change of the site, the movement law of the plate, and other information can be analyzed. The combined method is used to estimate the station noise, coordinate time series parameters, and the station velocity, and then the influence of the common mode error on the station noise, parameter estimation, and station velocity estimation is studied. The experimental results show that the average magnitude of white noise and colored noise decreases by 31.61% and 42.28%, respectively, which effectively reduces the uncertainty of parameter estimation and the uncertainty of station velocity estimation.
1. Introduction
GPS is widely used in the research of crustal deformation, geodynamics, seismology, and other fields due to its advantages of all-weather, high-precision, and real-time performance [1]. With the intensive construction of global permanent GPS continuous observation stations and the continuous improvement of GPS data processing accuracy, GPS can monitor global plate movement, regional crustal deformation [2], and active blocks and fault motions with higher precision, so that it can provide high-precision crustal surface displacement field changes with time [3]. With the development and progress of technology, the observation accuracy and data processing accuracy of GPS technology have been continuously improved [4–6], and its application in the field of geology is also changing with each passing day. These GPS continuous observation data distributed around the world are the basis of geoscience research [7].
GPS technology can provide a large amount of observational data on a global scale. Using these observational data, we can not only quantitatively study the movement of large regional plates but also study the movement and deformation of local plates and subplates [8]. Under the action of the earth’s internal and external forces, the surface of the earth’s crust will rise and fall, tilt, and dislocate. Crustal movements, such as faults, folds, and even earthquakes, are the process of continuous accumulation of crustal stress [9–11]. These crustal movements often lead to some geological disasters and various secondary disasters, such as land subsidence [12]. Various disasters caused by crustal movement have serious adverse effects on national economic development and people’s lives and properties [13]. Therefore, effective monitoring of crustal deformation and analysis of its movement laws have great theoretical and practical significance for understanding regional crustal deformation and disaster prevention and mitigation [14].
In order to study the high-precision data processing method of earthquake monitoring and conduct data analysis and research on the laws of global plate motion based on this, this paper constructs a real-time seismic data analysis model, which is based on 6 years of data from more than 700 GPS stations around the world. From the relative accuracy of the baseline, the repetition rate of the baseline, the numerical value of the baseline solution, the accuracy of the network adjustment solution, and the external coincidence accuracy of the coordinates, the subnetwork division by latitude and longitude and the subnetwork division by plate are compared and analyzed. The results are as shown: when , the value of m is large, which is not conducive to the execution efficiency of the program, so it is not listed in the table. The average relative accuracy of the baselines in the , , and directions of the solution results in solution is , , and , and the average solution time is 101.2 minutes; the average relative accuracy of the baselines in the , , and directions of the results is , , and , and the average solution time is 102.1 minutes; the calculation accuracy and solution efficiency are slightly better.
2. Related Work
In recent years, GPS technology has been more and more widely used in earth science, astronomy, and other fields, and the use of GPS technology to study the law of crustal movement has made great progress. In a word, GPS technology has made outstanding contributions to monitoring global and regional crustal movement with its advantages of high observation accuracy, low equipment cost, and all-weather large scale [15].
Gurbuz et al. [16] studied based on the original pseudorange and carrier phase observations, using the Kalman filtering technology and estimating ionospheric and tropospheric parameters at the same time, to achieve long baseline, dynamic high-precision relative positioning, and high-precision relative positioning for single-difference models. The positioning technology is studied, and the receiver-related information such as the relative hardware delay of the receiver is also obtained while the baseline component is obtained. Yabe et al. [17] demonstrated the equivalence of nondifferential and differential models, that is, using the same observations and estimating the same unknown parameters, the results are consistent, and a unified data processing model was proposed. Li et al. [18] analyzed the application fields of sparse fast Fourier transform and the discussion of common methods. Finally, based on the above content as the theoretical basis, this paper proposed a parallel code phase acquisition algorithm based on sparse fast Fourier transform. The algorithm adopts the idea of the regional remainder theorem.
Akhoondzadeh et al. [19] first proposed a fast acquisition algorithm for GPS satellite signal C/A code based on FFT (fast Fourier transform algorithm). For example, a fast acquisition method of pseudorandom noise based on piecewise-matched filter and FFT algorithm [20] averagely erases the signal sequence of the receiver and at the same time gives multiple local coding phases. The superposition of sequences improves the acquisition rate of satellite signals. The researchers propose a new FFT-based Galileo signal secondary code acquisition algorithm, which is based on a combination of serial and parallel searches, performing serial searching in the first stage and parallel searching in the second stage [21]. Scholars use circular correlation decomposition to improve the performance of FFT-based parallel code phase search to obtain GNSS signals, and propose to use the shift correlation of PRN (pseudorandom noise) codes to create a new local code, and use the standard FFT method to collect [22], shorten the acquisition time, overcome the ambiguity problem, and provide the speed of acquisition through key techniques such as frequency compensation, superposition processing, segmentation processing, and inversion position estimation [23].
3. Construction of Seismic Real-Time Data Analysis Model Based on GPS Monitoring
3.1. GPS Baseline Solution
The main purpose of using the GPS baseline to solve the carrier is to eliminate or weaken the influence of the ionospheric delay error, and the role of using the high-frequency carrier is to improve the navigation and positioning accuracy, because the high-frequency carrier can accurately measure the Doppler frequency shift and carrier phase . In the GPS navigation and positioning system , the purpose of the carrier is to transmit ranging code information and navigation messages. At the same time, in the carrier phase measurement, the carrier can also be used as a ranging signal .
Among them, the main purpose of the C/A code is to roughly measure the ranging and capture the ranging code of the precise code. The C/A code is only modulated on the carrier. Since the C/A code is not modulated on the carrier , it cannot accurately eliminate the ionospheric delay error . The P code modulates the carrier , respectively, which can completely eliminate the influence of the ionospheric delay error .
The ranging code is divided into fine code and coarse code. The rough code has a symbol rate of 5.11 MHz and a symbol length of 58.67 m; the fine code has a symbol rate of 0.511 MHz and a symbol length of 586.7 m. The ranging accuracy of the coarse code is 1/100 of the length of the symbol. For each GLONASS satellite, there is a C/A code with a length of 511 symbols to perform phase modulation on the carrier; GLONASS satellites use P codes similar to GPS satellite signals, but the P codes of GLONASS satellites are strictly confidential.
The deviation of the GLONASS satellite coordinate value in Figure 1 obtained by the integration method and the actual value increases with the increase of the integration time. Therefore, for the positioning with high-precision requirements, the integration length is preferably controlled within 30 minutes. It can be seen from the data in the table that when the value of the threshold coefficient gradually increases, the width of the fitting window tends to gradually decrease.
If other cycle slip detection methods can be combined to repair the cycle slips to within 5 cycles and then use the ionospheric residual method, all cycle slips can be detected. According to this idea, this paper uses the polynomial fitting method and the ionospheric residual method to jointly detect and repair the cycle slip. The specific operation is to first use the polynomial fitting method to detect the original carrier phase observations and repair the possible cycle slips of the original carrier phase observations to within 5 cycles and then use the ionospheric residual method for further processing to correctly detect and repair.
3.2. Seismic Network Distribution
In addition to the same advantages as the baseline solution of the seismic network, the overall solution mode also avoids the problem of inconsistent solutions of the same baseline in different time periods and the closed loop formed by the baselines in different time periods. The transform calculates the two-dimensional autocorrelation; finally, the one-dimensional covariance function can be obtained by using the rotation average of the two-dimensional autocorrelation function . Nonzero points in the data can be used to weight nonzero points in the data to regularize the resulting autocorrelation function when performing the power spectrum calculation.
Whether the receiver is single frequency or dual frequency, high precision or low precision, before solving the coordinate value of the monitoring station, it must first calculate the position of the observation satellites of the station and then observe the position of these satellites and the receiver according to the position of the station . The coordinate value of the monitoring station is calculated according to the relevant information such as quantity, and the calculation of the satellite position is calculated according to the relevant parameters in the satellite navigation message .
Assuming that GPS satellites and GLONASS satellites are synchronously observed at a certain time, and one GPS and GLONASS satellites are selected as the reference satellites of their respective systems, double-difference equations can be formed, which is better than that of satellites alone.
It eliminates the influence of the receiver clock and further reduces the influence of errors such as ionospheric delay and tropospheric delay, so that the ambiguity has the characteristics of an integer number of cycles and improves the relative positioning accuracy. The card reader of the marking machine itself couples the data and energy into the TAG (Thos A Grabbert) chip of the data acquisition system board through the antenna and can record the GPS data even when the data acquisition system board is not powered on. The DSP (digital signal processor) is used as the control part, and its control mechanism is that when the data in the TAG needs to be read, the DSP will be powered on and activated through the IO, which will reduce the power consumption of the entire board. A custom protocol stack is implemented in the FPGA (field programmable gate array). After the TAG is activated, the FPGA will read the time information in it and try to implement the acquisition algorithm in the FPGA based on this time.
In the network environment of Figure 2, in order to take advantage of distributed processing for baseline solution, the task itself must be decomposed, that is, parallel processing. The seismic wave data acquisition system includes the following parts: FPGA+DSP data processing mechanism, preamplifier filter circuit, A/D acquisition module, communication module, self-detection module, dip angle measurement module, temperature and humidity acquisition module, power transmission control module, RFID (radio frequency identification) transponder module, power supply detection module, power supply voltage conversion DC/DC module, and external GPS marking machine equipment. It can be seen that with the increase of sampling interval, the fluctuation range of ionospheric residuals tends to increase, which is also consistent with the actual situation.
The whole process of the baseline solution can be divided into three stages: data storage and transmission, data preprocessing, and baseline vector calculation. The data characteristics of different stages are different, and the distributed strategies adopted are also different. To play the best role of distributed computing, it is necessary to design different distributed computing methods according to the characteristics of different stages of baseline computing.
3.3. Real-Time Data Analysis
In the real-time data adjustment calculation of the GPS network, two methods of batch processing and sequential processing can be adopted. The batch method needs to prepare all the observation data of the processing period in advance and then process the data in batches. The sequential processing method means that with the arrival of the data of the first day, the data of the first day is processed first, and the adjustment results of the data of the first day are obtained. The adjustment results include the recorrection of the first day’s data and the correction of the second day, and these corrections are redistributed to the original approximate values.
It will also regard the data adjustment result of the first day as an approximate value and obtain a correction amount. Therefore, relatively speaking, the wavelet transform coefficient value of the signal must be larger than the wavelet coefficient value of the noise whose energy is dispersed and the amplitude is small, which is very prominent. If an appropriate threshold is selected, the noise can be effectively eliminated by thresholding the wavelet coefficients of the noise-containing signal .
GAMIT/GLOBK is a baseline processing software for positioning and orbit determination using GPS carrier phase observations. GLOBK uses the Kalman filter method for network adjustment. The monitoring project conducts a comprehensive monitoring of the tunnel area, tunnel low settlement, crack width, and tunnel displacement every three months, among which the GNSS technology is used for regional deformation monitoring. The equipment used is 4 sets of Trimble 5700 or R6, the monitoring station is observed for 7.8 hours, and the adoption rate is 15 seconds. The data processing is combined with the data of the continuous monitoring station. Figure 3 uses the corresponding software to process the GNSS data.
In this experiment, let the threshold coefficient start to take value, and the experiment obtains the minimum value of the fitting window width. The results are as shown: when , the value of is large, which is not conducive to the execution efficiency of the program, so it is not listed in the table. It can be seen from the data in the table that when the value of the threshold coefficient gradually increases, the width of the fitting window tends to gradually decrease, but the value of should not be too large, generally within the range of 2 to 9. From the perspective of program execution efficiency, when and , the program running time is 0.046850 seconds; when and , the program running time is 0.028068 seconds; when and , the program run time is 0.060400 seconds. Therefore, when and , the program execution efficiency is the highest. From the above discussion, it can be seen that in the case of only one day of data, different processing methods for different time periods not only can solve the Runge phenomenon caused by high interpolation order but also can significantly improve the overall interpolation accuracy. Although this interpolation processing method is a little troublesome, it brings about the improvement of the overall interpolation accuracy. In the case of only one day of data, this is a desirable method.
3.4. Error Factor Fitting
To use error factor filtering, it is necessary to first select the number of wavelet decomposition layers , then finely decompose the signal according to the frequency level, and shrink the high-frequency noise by selecting an appropriate threshold for each layer. The high-frequency coefficients from the first layer to the layer after thresholding are used to reconstruct the signal, and finally, the deformation signal of the structure is obtained. The multiresolution characteristics of wavelet and the adaptive noise canceller just complement each other.
Because the correlation of multipath errors for two consecutive days is only about 85% under good conditions and the random noise of the receiver is mainly high-frequency noise, which is irrelevant and inevitable, the adaptive noise canceller cannot suppress it. However, the structural deformation signal and multipath noise are mainly at low frequency , so it is considered to use the multiresolution property of wavelet to perform secondary filtering on the signal .
Since the GAMIT/GLOBK software limits the number of stations for a single day to 99, in the actual calculation, the number of points exceeds 65, and the calculation time is greatly increased. The data used in this paper has a total of about 700 stations. The rover R1 is placed at the opposite corner to the reference station, because the surrounding clearance is better, the multipath effect will not be generated, and the main error source is the system noise of the receiver. The mobile station R2 is set next to the warehouse, and the side of the warehouse is used to generate multipath signals. The results measured by R2 will be compared with R1. The three sampling intervals can detect the cycle slip value of 1 cycle, which shows that increasing the sampling interval has no obvious effect on the ability to detect cycle slips. All brackets are installed on the top. The base is set to keep the antenna level, the sampling frequency of the receiver is set to 1 Hz, and the sampling time is 90 min.
The measured data is analyzed by the joint denoising method proposed in Figure 4, and the filtering results are satisfactory. The detection of singularity is mainly through the wavelet multiscale decomposition of the signal, and the high-frequency part of each layer after the decomposition is observed. The protruding point of the waveform of the high-frequency part is the singularity point. In the short-baseline differential positioning, it can be considered that the tropospheric errors at both ends of the baseline are basically the same, and the influence of this error can be basically eliminated by the carrier phase difference.
The GPS data postprocessing software is used to eliminate the periodic components and common mode errors and other nonstructural deformation information in the continuous station displacement sequence, so that the horizontal movement speed of the continuous observation station in the GPS crustal monitoring is improved. At this time, the tropospheric delay error cannot be eliminated by differencing even under short-baseline conditions. When there is a sudden change in the signal, the coefficients after wavelet transformation have the maximum value of modulus. Therefore, the time when the signal mutation occurs can be determined by detecting the maximum point of the modulus. Using the adaptive noise cancellation filtering method for two consecutive days, the standard deviation error of the signal can be reduced from the original 13.48% to 4.26%, and the final error is reduced to 1.42% after the second wavelet denoising.
4. Application and Analysis of Seismic Real-Time Data Analysis Model Based on GPS Monitoring
4.1. GPS Data Sampling
The data from 001-003 days in 2019 of about 700 globally distributed GPS continuous operation base stations are selected for the experiment, and the above two subnetwork division schemes are used for data calculation. The repeatability, the accuracy of the three-dimensional coordinates, and the compliance with the weekly solution are compared. Note that this file does not list the number of estimated parameters, and we read the estimated parameters directly and finally count how many estimated parameters there are.
The first field of the parameter line is the parameter name, the second number is the parameter index, the third and fourth numbers are the effective start epoch of the parameter and end epoch, the fifth number is the parameter prior value, the sixth number is the parameter constraint value, and the seventh number is the group number corresponding to the parameter. The computer equipment used in the two subnetting schemes is DELL R710, a 2U rack-mounted workstation, and its basic configuration parameters are Intel Xeon 4-core processor, 8 GB DDR-3 1066 MHz memory. Table 1 uses GAMIT/GLOBK10.6 to perform data calculation on all subnets in their respective subnetting schemes simultaneously under the Centos 7 system.
The method of subnetting by latitude and longitude will have a uniform number of stations in the subnet, and the effective time to solve the entire GPS network (the time taken to solve the slowest subnet) will be lower. The reason for the uneven number of stations in each subnetwork divided according to plate distribution is that the subnetwork is limited by plate distribution, and some plates have fewer stations. The waveform with a sampling interval of 15 s fluctuates within 5 weeks, and the detection ability of the polynomial fitting method for the observation data with a sampling interval of 15 s is 5 weeks. According to the above discussion, it can be concluded that if the sampling interval of the receiver is set within 15 s, the cycle slip value of the original observation data can be repaired within 5 weeks by using the polynomial fitting method.
Figure 5 uses the built-in mk-net coordinate conversion tool to process the GAMIT/GLOBK data in the catalog after updating the tables catalog org file to extract the uric line data for coordinate transformation to obtain the required a priori coordinate file format and then perform st-filter processing. The results obtained from the org file have better accuracy. The wavelet coefficients at different scales are similar to the wavelet coefficients of the measured wind speed, with many peaks and basically occurring in the corresponding time period, reflecting the internal hierarchical structure of the turbulence. The wavelet coefficient probability density function comparison of wavelet simulation wind speed shows that the wind speed simulated by harmonic stacking is Gaussian, while the wind speed simulated by wavelet is non-Gaussian due to the existence of intermittent, which is consistent with the measured wind speed, and the inverse distance interpolation algorithm is used to eliminate the nonstructural noise in the GPS mobile observation regional stations, so that the horizontal movement speed in the mobile monitoring of GPS regional stations is not affected. The kurtosis between different scales shows an obvious decreasing trend, especially in the first four scales, the probability density distribution is already quite different from the Gaussian distribution, and the kurtosis of scale 2 reaches 12.4333. In addition, it can be seen from the skewness analysis that the data distribution of the measured wind speed and the simulated wind speed is basically symmetrical.
4.2. Real-Time Data Measurement
In this paper, the above two subnetting schemes are used to analyze the baseline repetition rate from four aspects, and perform statistical analysis from the aspects of maximum value, minimum value, and average value. The data sampling interval was set to 15 seconds, and the observation period was 3 hours. Figure 6 selects the phase observations of the first 400 epochs of the satellites PRN2 and PRN4 (PRN represents the pseudorandom number of GPS satellites) synchronously observed by the two stations as the research objects.
It can be clearly found that there is a big difference with the measured wind speed. Although the wavelet coefficient shows a certain intermittent in the low frequency stage, it is basically Gaussian white noise in the high-frequency stage, and there is no difference between different scales like the measured wind speed. From the wavelet coefficient probability density function and the kurtosis and skewness analysis results of the harmonic simulated wind speed, it can be seen that the harmonic superposition method is in good agreement with the Gaussian distribution, its kurtosis fluctuates around 3, and the skewness fluctuates around 0.
The maximum error of the solution results of the three components of the two network distribution schemes is less than 2 cm, the minimum value is below 2 mm, and the average value is between 4 mm and 5.2 mm. The calculation accuracy is high and meets the requirements of crustal movement research; in contrast, the solution accuracy of the subnetting scheme by latitude and longitude is slightly better. It can be seen that most of the differences in the adjustment results of the two subnetting schemes are less than 5 mm, and the two have a high degree of agreement, which also proves the reliability of the calculation results from the side.
4.3. Monitoring Coordinate Simulation
In order to facilitate the management of GNSS deformation monitoring coordinate result data and the corresponding operations, SQL Server 2008 is used as the database for data management. This module includes functions such as submission and storage, data addition, data modification, data deletion, Excel import, and data clearing. Since the power supply of most GPS stations in the earthquake area was interrupted 52 seconds after the earthquake, the power supply and data collection were resumed after about four hours, so there was no GPS data during the power outage. Due to the large deformation caused by the earthquake, the observation data before and after the earthquake are divided into two periods, namely the period before the earthquake and the period after the earthquake. In order to facilitate the analysis of GPS data, it is mainly considered that most of the stations in the earthquake area have relatively large movements during the earthquake, so relatively stable GPS stations are required as reference stations and data fusion.
Before data query, first set the name of the monitoring point to be queried and the time period to be queried. The monitoring point name in Figure 7 can be selected more than one. The time period of the query is divided into two types: one is to set the start time, the end time defaults to the present time, and we query all monitoring data of GPS1 and GPS3 monitoring points; the second is to set the query start and end time, respectively. At the same time, the nonstructural noise parameter information such as the annual term and the semiannual term in the displacement time series after removing the common mode error is fitted and finally removed from the displacement sequence. The implementation of the first two modules is to specify the name and period of monitoring points in advance and selectively modify and delete the data in the SQL (structured query language) database; data clearing is to clear the specified database of all monitoring data.
4.4. Example Application and Analysis
In view of the characteristics of many monitoring points, long monitoring period, and large amount of data in engineering buildings, the monitoring data query function has been developed, which is convenient for users to screen data according to the needs of the actual situation, improves the efficiency of data sorting, and sifts the screened data. The primary indicator of the first criterion is the magnitude of the baseline component error. If it is larger than your preset prior constraints on station coordinates and orbital parameters, then a quick look at the file or the autcIn-sum file will reveal a lot of data lost by it. If the data is randomly distributed, the prior weights are also correct and will be close to unity. In practice, using the default weighting method, a good solution can yield about 0.25. If it is greater than 0.5, it means that the cycle slip is not eliminated or is related to additional biased parameters or occurs, if the final batch solution meets both criteria, it is usually not necessary to carefully observe any other outputs.
The seismic wave data acquisition system uses the high-precision 24-bit ADC (analog to digital converter) in Table 2 as the core of data acquisition, which greatly improves the dynamic range of the acquired signal and meets the signal acquisition requirement of a dynamic range of 120 dB. The ADC is driven by the FPGA digital logic, and the data is transmitted to the DSP through the interface for data processing.
The FPGA frames the processed data and transmits it to the terminal or the previous node through the POE communication interface. In order to enable users to visually check whether the displacement of the monitoring point exceeds the specification tolerance or early warning value, a tolerance setting module is added on the basis of the function of drawing the displacement trend graph. The second is the self-defined tolerance. According to the needs of different engineering buildings and different user levels, the tolerance value can be set by yourself. The displacement trend graph can be zoomed in and out, saved, and checked for information such as the phase and displacement of the monitoring point. It can be seen that the displacement of the GPS4 monitoring point in the ninth phase of the -axis is 1.8 mm. In this paper, GAMIT10.6 is used for the baseline solution, the solution result is a single-day solution, and the number of calculations for the 18 subnetworks is about 20,000 per day.
Since the application of GPS analysis requires uniform sampling, for the uneven sampling in the coordinate time series of the reference site, interpolation processing is required to obtain a uniformly sampled coordinate time series. The interpolation method used in Figure 8 is singular spectrum iterative difference. Due to the large number of stations in this paper, 4 stations with obvious data missing are selected from all stations to show the interpolation effect. It can be seen that the four stations have obvious data missing before the interpolation. After the iterative interpolation, the coordinate time series becomes continuous and smooth, which is more in line with the data requirements of crustal movement research. Trimble 5700 or R6 was used, and the adoption rate was 15 s; continuous monitoring data was officially collected from July 2019, the Leica 902 dual-satellite receiver and antenna were used, and the adoption rate was 15 seconds. It is found that the maximum speed difference in the north-south direction is 0.6 mm, the minimum is 0 mm, and the root mean square of the speed difference is 0.26 mm; the maximum difference in the east-west direction is 0.5 mm, the minimum difference is 0 mm, and the speed difference RMS is 0.19 mm.
This paper uses 5% as the threshold, and data with a missing rate of more than 5% will be eliminated. After analysis, it was found that about 160 stations did not meet the requirements and were eliminated, and the eliminated stations accounted for 22.9% of all stations. The remaining number of stations is about 536. The difference between the filter value and the measured value of GPS1 and GPS4 is within ±0.4 mm, indicating that the Kalman filter can better reflect the actual change of GNSS deformation monitoring and can truly reflect the real-time dynamic data of the deformable body. The deformation monitoring data has a good denoising effect.
5. Conclusion
In this paper, through the analysis of actual data, the deformation trend of seismic information can be obtained by GNSS, which can meet the requirements of high precision and real-time seismic monitoring. In view of the characteristics of long deformation monitoring time and large amount of data, the SQL Server 2008 database is studied, and the GNSS monitoring data is efficiently and real-time managed, and auxiliary functions such as modification and Excel batch import are designed; in order to visualize the axial and point displacement trends of GNSS monitoring points, the effectiveness and real-time update of the controls are verified through actual development. It can be seen that the root mean square of the velocity difference in the vertical direction is 1.97 times that of the horizontal direction. At the same time, the single-day solution was calculated for 6 years of observation data from more than 700 GPS stations around the world, reaching millimeter-level accuracy. On the basis of the research on the subnetwork division scheme, GAMIT/GLOBK10.6 is used to calculate the baseline of each subnetwork, and the joint network adjustment is carried out for the whole network. It can be seen from the calculation results that the relative accuracy of the baseline reaches the order of ; the calculated value of the baseline is around 0.18; the average repetition rate of the baseline is 0.0028 m; the average value of the network adjustment is below 6 mm, and the maximum parts are around 4 mm. It can be seen that the solution results can meet the research needs of this paper.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The author declares no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
The work presented here was supported by the National Key Research and Development Program of China under grants 2021YFC3000705-06 and 2017YFC1500502-05.