A regional troposphere delay modeling method based on GNSS observation data and fusion model

Abstract

This invention discloses a regional tropospheric delay modeling method based on GNSS observation data and a fusion model, comprising: collecting and preprocessing raw data from GNSS observation stations in the target area; dividing the data into modeling and independent test datasets; calculating the actual ZTD (Zero-Temperature Delay) for each observation station; constructing an empirical function model and solving for the model parameters as training labels for the subsequent initial prediction model; constructing and training an initial prediction model with the latitude, longitude, and elevation of the observation station as input and the prediction model parameters as the output target; constructing and training a residual improvement model with the initial predicted ZTD as input and the residual between the initial predicted ZTD and the actual ZTD as the output; combining the initial prediction model and the residual improvement model to obtain the final predicted ZTD for unknown observation stations; and verifying the accuracy and regional adaptability by comparing the error index with existing models. This method can better fit the nonlinear relationship between ZTD and spatiotemporal information, providing a more accurate solution for tropospheric delay modeling.

Description

Technical Field

[0001] This invention relates to the field of precise positioning and atmospheric parameter inversion technology for global satellite navigation systems, and particularly to a regional tropospheric delay modeling method based on GNSS observation data and fusion model. Background Technology

[0002] Tropospheric zenith delay is a core error source in GNSS signal propagation, and its accuracy directly determines the reliability of positioning results. In the elevation direction, a 1cm tropospheric delay error can lead to a positioning deviation of approximately 10cm, severely restricting the implementation of high-precision applications such as urban mapping benchmark establishment and large-scale engineering construction. With the deepening application of GNSS technology in meteorology, geology, and other fields, there is an urgent need to develop tropospheric delay modeling methods adapted to regional characteristics.

[0003] However, while many existing techniques exist for zenith delay modeling in the troposphere, none can meet the high-precision requirements of regions like China, which possess both complex topographical types and monsoon climate characteristics. Existing techniques can be broadly categorized into the following four types: 1. Measured meteorological parameter-driven models (such as the Saastamoinen model): These models rely on synchronous observation data such as surface temperature, air pressure, and humidity, and calculate the delay through physical formulas. However, in areas such as the Qinghai-Tibet Plateau and the Northwest Desert, the density of meteorological observation stations is less than 0.1 per 1000 km², and extreme weather can easily cause sensor failure, making the models unable to be effectively applied in about 60% of the country's land area.

[0004] 2. Global empirical models (such as GPT3 and CTrop models): These are built based on global meteorological statistics and do not fully consider regional topographic and climatic differences. In the low-latitude monsoon region of southeastern China, sudden changes in water vapor in summer lead to an RMSE as high as 45-50 mm, which is more than 70% different from regional models. In the Qinghai-Tibet Plateau region of China, due to defects in the elevation normalization algorithm, the cumulative error can reach more than 30 mm, which cannot meet the requirements for precise positioning.

[0005] 3. Regional empirical models (such as the SHAtrop model): optimized for East Asia, but still have significant limitations: during the spring and autumn monsoon transition period, ZTD predicts MAE of 4-5 cm; in high-altitude areas (>3000m), the interpolation accuracy distortion rate exceeds 20% because the influence of topographic relief on the vertical structure of the atmosphere is not taken into account; in low elevation angle (<15°) observation scenarios, the model bias is more than 3 times larger than the measured value.

[0006] 4. Single machine learning model (such as basic BP neural network): Although it can fit nonlinear relationships, it has inherent defects: it requires a very large amount of training data (millions of samples), and it is prone to overfitting in sparse data areas; it does not incorporate physical periodic characteristics, and its generalization ability to extreme weather events (such as typhoons and cold waves) is weak, with prediction error fluctuations reaching 20-30mm.

[0007] Furthermore, the China Crustal Movement Observation Network (CMONOC), also known as the Land State Network, is the world's largest regional GNSS observation network. It consists of 25 reference stations, 56 basic stations, and over a thousand regional stations, providing continuous observation data with millimeter-level accuracy. It covers the entire region of China from 20° to 50°N and 75° to 135°E. However, existing technologies lack suitable methods for adapting to this network. Therefore, how to integrate physical mechanisms and data-driven advantages based on this valuable data resource to construct a high-precision and highly adaptable regional tropospheric delay model has become a key issue that urgently needs to be addressed in this field. Summary of the Invention

[0008] To address the shortcomings of existing technologies in regional tropospheric delay modeling with complex landforms and monsoon climate characteristics, such as insufficient accuracy and poor spatiotemporal adaptability, this invention provides a regional tropospheric delay modeling method based on GNSS observation data and a fusion model. The method aims to achieve high-precision ZTD prediction under complex conditions through a three-level architecture of "empirical model fitting - neural network parameter prediction - residual correction optimization" based on the fusion model of empirical models and multi-stage machine learning.

[0009] To achieve the above-mentioned technical objectives, the present invention provides the following technical solution: A method for regional tropospheric delay modeling based on GNSS observation data and a fusion model, specifically including: S1. Obtain long-term raw observation data from GNSS continuously operating observation stations that uniformly cover the target area. Process the data to obtain the actual ZTD of each observation station and divide it into training dataset and test dataset. Select the observation data of the observation stations in the training dataset that meet the preset requirements for data integrity to form the prior dataset. S2. An empirical function model for fitting the tropospheric zenith delay (ZTD) of each observation station is constructed in the form of a periodic function. Based on the prior data set, n undetermined parameters of the empirical function model of the known observation station are calculated. S3. Using the longitude, latitude, and elevation of each observation station as input and the undetermined parameters of the empirical function as the output target, construct an initial ZTD prediction model consisting of n independent BP neural networks. Each BP neural network is responsible for predicting one undetermined parameter. Use the n undetermined parameters of the known empirical function model of the observation station as training labels to train the initial ZTD prediction model. S4. Based on the trained initial ZTD prediction model, obtain n undetermined parameters of the empirical function model of each observation station in the training dataset. Then, obtain the initial ZTD prediction of each observation station through the empirical function model of each observation station. Subtract the initial ZTD prediction from the actual ZTD to obtain the residual prediction value. Construct and train a residual improvement model with the spatiotemporal characteristics of the observation station, the initial ZTD prediction as input, and the residual prediction value as output. S5. For an unknown observation station in the current region, the initial predicted ZTD is obtained based on the empirical function model and the trained ZTD initial model. The residual sequence is obtained through the residual improvement model. The sum of the initial predicted ZTD and the residual predicted value is used as the final predicted ZTD of the unknown observation station. S6. Using error quantification evaluation metrics, the model is compared and verified with existing tropospheric delay models on independent test datasets to evaluate the overall model's accuracy and regional adaptability.

[0010] Further, in step S1, the acquisition of long-time-series raw observation data from GNSS continuously operating observation stations that uniformly cover the target area, and the processing to obtain the actual ZTD of each observation station, specifically involves: Based on the observation data of the land state network, in the target area with complex landforms and monsoon climate characteristics, GNSS observation data from N continuously operating observation stations are collected as raw observation data; cycle slip repair, gross error removal and spatiotemporal alignment are performed on the raw observation data of each observation station in sequence to obtain the actual ZTD of each observation station.

[0011] Further, in step S2, the empirical function model includes: a fundamental fixed component in the tropospheric zenith delay that does not change periodically with time, and a periodic fluctuation component in the tropospheric zenith delay that changes with different time scales, expressed by the formula: ; Where c is a constant term, representing the fundamental fixed component in the tropospheric delay that does not change periodically with time; , , These represent the periodic fluctuation components in the tropospheric zenith delay that vary with the semi-annual, annual, and daily cycles, respectively. , , These are the amplitude coefficients of each periodic term, representing the magnitude of the periodic fluctuations; , , These are the phase parameters for each periodic term, reflecting the time position of the fluctuation; doy is the yearly cumulative day. , , , , , c represents the parameters of the empirical function model to be solved.

[0012] More specifically, in step S2, the process of calculating the n undetermined parameters of the known observation station empirical function model based on the prior dataset is as follows: The nonlinear least squares method is used to fit the n undetermined parameters of the empirical function model. The sum of squared residuals between the actual ZTD of each observation station and the ZTD calculated by the model is used as the fitting objective function, and the undetermined parameters are optimized by the Levenberg-Marquardt iterative algorithm.

[0013] Furthermore, step S3 specifically includes: The latitude (B), longitude (L), and elevation (H) of a single observation station are used as inputs. Each BP neural network consists of an input layer, a hidden layer, and an output layer. The input layer has three nodes, corresponding to the latitude, longitude, and elevation components, while the output layer has one node, corresponding to the model parameters that each BP neural network is responsible for predicting. The number of nodes in the hidden layer is determined through five-fold cross-validation. The activation function between the input and hidden layers is the tansig function, and the activation function between the hidden and output layers is the purelin function. After normalizing the three components of the input features, the Levenberg-Marquardt algorithm is used for training. After each training iteration, the mean squared error (MSE) between the output model parameters and the training labels is calculated and recorded. Three training convergence conditions are set: the calculated MSE is less than a preset threshold, the maximum number of iterations is reached, and the prediction accuracy does not improve after multiple consecutive training iterations. Training is considered complete when any one of the convergence conditions is met.

[0014] Further, the residual correction model described in step S4 employs a backpropagation (BP) neural network, using the latitude, longitude, elevation, annual day, hour, and ZTD initial prediction values ​​of the observation station as six-dimensional input features, and outputting residual prediction values. The residual correction model has six nodes in the input layer, one node in the output layer, and the number of nodes in the hidden layer is determined through five-fold cross-validation. The activation function between the input layer and the hidden layer is the tansig function, and the activation function between the hidden layer and the output layer is the purelin function. The Levenberg-Marquardt algorithm is used for training. After each training iteration, the mean squared error (MSE) between the output residual prediction value and the training label is calculated and recorded, and three training convergence conditions are set: the calculated MSE is less than a preset threshold, the maximum number of iterations is reached, and the prediction accuracy does not improve after multiple consecutive training iterations. Training is considered complete when any one of these convergence conditions is met.

[0015] Furthermore, the error quantification evaluation indicators include: the root mean square error (RMSE) and the mean absolute error (MAE) between the final predicted ZTD and the actual ZTD of the independent test set observation stations.

[0016] Furthermore, this invention also discloses a regional tropospheric delay modeling system based on GNSS observation data and a fusion model, which specifically includes: The data management module is used to acquire raw GNSS observation data of the target area for preprocessing, divide the training dataset and the test dataset, and obtain the actual tropospheric zenith delay (ZTD) of each observation station. The empirical modeling and solution module is used to establish an empirical function model to fit the tropospheric zenith delay (ZTD) of each observation station, and to solve the model parameters of the empirical function model corresponding to the actual ZTD based on the training dataset, which are then used as training labels. The initial prediction module is used to construct an initial ZTD prediction model that takes the spatial coordinates of the observation station as input and solves the parameters of the empirical function model as the objective. After training based on the training labels and training dataset, it performs initial ZTD predictions for the observation stations in the current area. The residual correction module is used to build and train a residual improvement model that takes the spatiotemporal information of the observation station and the initial predicted ZTD as inputs and outputs the residual prediction value between the initial predicted ZTD and the actual ZTD. The online prediction module receives the spatial coordinates of the observation station to be observed, and outputs a high-precision final prediction ZTD by calling the trained neural network model and the residual improvement model. The model validation and evaluation module is used to call the test dataset to evaluate and compare the accuracy of the final ZTD prediction value.

[0017] By employing the above technical solution, the method proposed in this invention has at least the following beneficial effects: 1. Significantly superior accuracy compared to existing models: The overall RMSE of the test set is 27.20 mm and MAE is 19.67 mm, which are 7.5% and 8.3% lower than the CTrop model, respectively, 1.0% and 1.5% lower than the GPT3 model, and 1.5% and 1.8% lower than the SHAtrop model; In the low-latitude monsoon region (20°-30°N), the RMSE is 24.8% lower than CTrop, reaching the leading level in the region; 2. Excellent spatiotemporal adaptability: Addressing the differences in terrain and climate within the target area, in high-altitude regions (>2000m), MAE is reduced by 10.9% compared to CTrop, resolving the issue of insufficient elevation normalization; during the period of drastic water vapor changes in summer (June-August), RMSE is only 43.97mm, more stable than GPT3's 45.41mm, and the error fluctuation during the spring-autumn transition period is reduced by 12%; 3. No need to rely on real-time meteorological data: High-precision forecasts can be achieved solely through the spatiotemporal information of the stations, breaking through the dependence of traditional models on meteorological sensors. It can still be stably applied in remote areas, plateaus and mountains, and other areas with weak observation capabilities, with a wide range of applications and large coverage. 4. Outstanding generalization ability and robustness: The two-stage architecture of "empirical function model + residual correction model" combines the advantages of physical mechanisms and data-driven approaches, avoiding the overfitting problem of a single machine learning model. It reduces the prediction error of extreme weather events such as typhoons and rainstorms by more than 15%, and significantly improves the robustness of the model. 5. Wide range of applications: It can also be ported to GNSS precise point positioning (PPP), network RTK, atmospheric water vapor inversion and other fields. It can shorten the convergence time of PPP elevation direction by more than 15% and improve the water vapor inversion accuracy to the 2mm level. It is expected to provide strong technical support for the construction of urban surveying benchmarks and geological disaster monitoring. Attached Figure Description

[0018] The accompanying drawings, which are included to provide a further understanding of this application and form part of this application, illustrate exemplary embodiments and are used to explain this application, but do not constitute an undue limitation of this application. In the drawings: Figure 1 This is an overall flowchart of the method proposed in this invention; Figure 2 A graph showing the ZTD change over a single observation station over five years; Figure 3 This is a flowchart illustrating the initial prediction model for ZTD in the method proposed in this invention. Figure 4 This is a schematic diagram of the structure of each BP neural network in the initial prediction model of ZTD; Figure 5 This is a schematic diagram of the residual correction model in the method proposed in this invention; Figure 6 This is a schematic diagram of the three-level fusion model of "empirical model fitting - neural network parameter prediction - residual correction optimization" in this invention. Detailed Implementation

[0019] To better understand the purpose, technical solution, and advantages of this invention, it is now combined with... Figures 1-6 The present invention will be described in detail below with reference to specific embodiments. It should be noted that this detailed description should not be considered as a limitation of the present invention, but rather as a more detailed description of certain aspects, features, and embodiments of the present invention.

[0020] In this invention, unless otherwise stated, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art. While only preferred methods and materials have been described herein, any methods and materials similar to or equivalent to those described herein may be used in the implementation or testing of this invention.

[0021] Various modifications and variations can be made to the specific embodiments described in this specification without departing from the scope or spirit of the invention, as will be apparent to those skilled in the art. The specification and embodiments of this invention are merely exemplary.

[0022] This embodiment discloses a regional tropospheric delay modeling method based on GNSS observation data and a fusion model. Based on historical long-term GNSS observation data published by the Land State Network and a three-level fusion model with an architecture of "empirical model fitting - neural network parameter prediction - residual correction optimization," it achieves high-precision and high-reliability prediction of tropospheric zenith delay for target areas with complex landforms and monsoon climate characteristics, especially the Chinese region covered by the Land State Network. Compared with traditional single models, it can better handle nonlinear relationships and noise, and also has excellent generalization ability. Figure 1 As shown, the method proposed in this invention specifically includes the following steps: S1. Data Acquisition and Dataset Construction: Long-term raw observation data are acquired from continuously operating GNSS observation stations that uniformly cover the target area. In this embodiment, publicly available GNSS observation data from 2014 to 2019 across the country is acquired. This data covers 217 continuously operating observation stations, encompassing all landform types, including plains, plateaus, and coastal areas. Preprocessing, including cycle slip repair, gross error removal, and spatiotemporal alignment, is then performed to obtain the actual ZTD for each observation station. The observations from each station are then divided in a 6:4 ratio into a training dataset (including GNSS observation data from 155 observation stations from 2014 to 2018, totaling 1,883,076 samples) and a test dataset (including GNSS observation data from 62 observation stations only from 2019, totaling 136,958 samples). The observation data from observation stations in the training dataset that meet the preset data completeness requirements are selected to form the prior dataset. S2. Construction and parameter calculation of multi-period empirical model: An empirical function model for fitting the tropospheric zenith delay (ZTD) of each observation station is constructed in the form of a periodic function. Based on the prior data set, n undetermined parameters of the empirical function model of the known observation station are calculated. In this embodiment, GNSS observation data from observation stations with a data integrity of ≥80% (GNSS observation data from 139 observation stations selected from 155 observation stations in the training dataset) were first selected as the prior dataset. Tropospheric zenith delay (ZTD) variation maps for these observation stations were then plotted for each year between 2014 and 2018. The ZTD variation map for a single observation station is shown below. Figure 2 As shown, the tropospheric zenith delay exhibits the following variation pattern: The annual variation range was relatively consistent over the five years, with a relatively stable start to the year and small fluctuations. Significant fluctuations occurred between 150 and 250 cumulative days (ZTD), with several peaks significantly higher than surrounding values. This may be related to changes in meteorological conditions during this period, such as temperature, humidity, and air pressure, indicating the presence of factors influencing tropospheric zenith delay, such as local meteorological system activity. ZTD exhibits annual and semi-annual cycle characteristics, closely related to the annual cycle of climatic elements such as temperature and precipitation. For example, ZTD is generally higher during warm and humid seasons and relatively lower during cold and dry seasons. The diurnal cycle is also an important factor; the alternation of day and night causes periodic fluctuations in surface temperature and air pressure, thus affecting tropospheric zenith delay. For instance, daytime sunshine increases temperature and enhances water vapor activity, while the opposite occurs at night, indicating that the diurnal cycle plays a crucial role in characterizing the short-term dynamics of tropospheric delay.

[0023] Based on these patterns, this embodiment constructs an empirical model that includes annual, semi-annual, and daily cycles. It includes a fundamental fixed component in the tropospheric zenith delay that does not change periodically with time, and periodic fluctuation components in the tropospheric zenith delay that vary with different time scales. The formula is expressed as follows: ; Where c is a constant term, representing the basic fixed component in the tropospheric delay that does not change periodically with time, and is an overall offset. , , These represent the periodic fluctuation components in the tropospheric zenith delay that vary with the semi-annual, annual, and daily cycles, respectively. , , These are the amplitude coefficients of each periodic term, representing the magnitude of the periodic fluctuations; , , These are the phase parameters for each periodic term, reflecting the time position of the fluctuation; doy is the yearly cumulative day. , , , , , c represents the parameters of the empirical function model to be solved; this model helps to balance complexity and accuracy.

[0024] Furthermore, it should be noted that the year-day (doy) input to the empirical function model is not always an integer; it varies in decimal places depending on the hour for each date. Since the ZTD time resolution calculated based on GNSS data is per hour, the "hour" parameter also needs to be input. Figure 3 As shown, the empirical function model will directly convert it to the annual day (doy).

[0025] This method fits the tropospheric zenith delay (ZTD) of a single observation station using empirical function model parameters. It already reflects the periodic fluctuations of ZTD using periodic terms across three time scales. To more accurately describe the relationship between the spatiotemporal characteristics of the observation station (location, day of day) and ZTD, training labels are first constructed using data from known observation stations. In this embodiment, the least squares method is used to calculate seven model parameters from the empirical function models of 139 known observation stations with high data integrity, which are then used as training labels. These labels, along with the latitude, longitude, and elevation information of these observation stations, serve as the training file for the initial ZTD prediction model in step S3.

[0026] S3, ZTD initial prediction: such as Figure 3 As shown, using the longitude, latitude, and elevation of each observation station as input and the calculation of the undetermined parameters of the empirical function as the output target, an initial ZTD prediction model composed of n independent BP neural networks is constructed, with each BP neural network responsible for predicting one undetermined parameter; the n undetermined parameters of the known empirical function model of the observation station are used as training labels to train the initial ZTD prediction model. In a preferred embodiment, step S3 specifically includes: The inputs are the latitude B, longitude L, and elevation H of a single observation station; each BP neural network is as follows: Figure 4 As shown, it includes an input layer, a hidden layer, and an output layer; the input layer X has 3 nodes, corresponding to the latitude, longitude, and elevation components, and the output layer... There is one node, corresponding to the model parameters that each BP neural network is responsible for predicting; the number of nodes in the hidden layer (each node is shown as m1, m2, ..., m in the diagram) n The number of nodes in each hidden layer is determined by five-fold cross-validation (based on experience, the candidate range for the number of nodes in each hidden layer is 4-30); the activation function between the input layer and the hidden layer is the tansig function, and the activation function between the hidden layer and the output layer is the purelin function;

[0027] After normalizing the three components of the input features, the Levenberg-Marquardt algorithm is used for training. After each training iteration, the mean squared error (MSE) between the output model parameters and the training labels is calculated and recorded. Three training convergence conditions are set: the calculated MSE is less than a preset threshold, the maximum number of iterations is reached, and the prediction accuracy does not improve after multiple consecutive training iterations. Training is considered complete when any one of the convergence conditions is met.

[0028] It should be noted that in this embodiment, the lack of improvement in prediction accuracy after multiple training iterations means that the 139 sets of input data are divided into training set, validation set and test set in a ratio of 7:1.5:1.5 for each training session. When using the validation set for prediction, if the MSE value between the output model parameter value and the training label is difficult to improve after multiple training sessions (for example, although the MSE has not reached the preset threshold, it cannot be reduced further), then the training is considered to be complete.

[0029] In this embodiment, the 139 sets of data obtained in step S2 are divided into training, validation, and test sets in a ratio of 7:1.5:1.5 for each training iteration, resulting in training, validation, and test sample sizes of 97, 21, and 21, respectively. Since the selection of the number of hidden layer nodes is crucial to model performance during BP neural network training, too few nodes may lead to underfitting (failure to capture data features), while too many nodes may cause overfitting (high model complexity and poor generalization ability). Therefore, after each training iteration, the mean squared error (MSE) between the output model parameter values ​​and the training labels needs to be calculated and recorded. After multiple training iterations, combined with five-fold cross-validation, the optimal number of hidden layer nodes in the BP neural network corresponding to each model parameter is determined. The optimal number of nodes is shown in Table 1 below. Table 1. Optimal number of hidden layers for each parameter

[0030] In this embodiment, the optimal number of nodes is then brought back to each group of BP neural networks. The latitude, longitude, and elevation of a known observation station are input, and the Pearson correlation coefficient R between the model parameters (predicted values) obtained through each group of BP neural networks and the model parameters (actual values) obtained by the least squares method in step S2 is calculated. The results are shown in Table 2 below: Table 2 Correlation coefficients of each model parameter fit

[0031] As can be seen from Table 2, the prediction performance of each group of BP neural networks in the initial ZTD prediction model of this application varies for different model parameters: The correlation coefficients R for the model parameter c are very high, close to 1, indicating that the trained BP neural network has a good regression fit for c. The BP neural network can capture the complex relationship between the parameter c and the three-dimensional spatial position of the observation station, and thus accurately predict the parameter. Since c reflects the fixed and unchanging basic level of tropospheric zenith delay, it proves that the initial ZTD prediction model can accurately predict the long-term average level of tropospheric zenith delay. The correlation coefficient corresponding to the amplitude coefficient a1 of the semi-annual cycle term is 0.64, which is relatively low. This indicates that the prediction accuracy of the BP neural network for a1 is relatively poor. This may be because the semi-annual cycle term reflects the tropospheric zenith delay change, which is affected by a variety of complex factors. In addition to geographical location, it may also be related to short-term abnormal changes in atmospheric circulation, local climate events, etc. These factors are difficult to fully and accurately reflect by latitude, longitude and elevation alone. The correlation coefficient corresponding to the annual cycle amplitude coefficient a2 is relatively high. a2 corresponds to the annual cycle fluctuation amplitude of the tropospheric zenith delay. The annual cycle change is relatively stable because it is mainly affected by seasonal changes. Latitude, longitude, and elevation can reflect information such as the climate zone of the region to a certain extent and are related to seasonal changes. Therefore, the BP neural network is sufficient to make good use of these input information to learn the relationship between the annual cycle fluctuation amplitude and geographical location, thereby predicting the parameter a2 more accurately. The correlation coefficient corresponding to the daily periodic amplitude coefficient a3 is relatively moderate. a3 represents the daily periodic fluctuation amplitude of ZTD. The change in the daily periodic fluctuation amplitude is not only related to latitude and longitude and elevation, but also affected by diurnal temperature difference and local meteorological conditions (such as humidity and temperature changes in the morning and evening).

[0032] Regarding the phase parameters of the three periodic terms, the overall prediction performance of the BP neural network is acceptable: The phase parameter b1 of the semi-annual cycle determines the starting position or initial phase of the tropospheric delayed semi-annual cycle fluctuation, reflecting the "lead" or "lag" of the tropospheric delayed fluctuation on the time axis within six months. The poor prediction effect of b1 is because the tropospheric delayed fluctuation of the semi-annual cycle is affected by a variety of complex factors. In addition to geographical location (latitude, longitude, and elevation), it is also related to the semi-annual scale anomalies in atmospheric circulation and the non-uniformity of solar radiation, which makes it impossible to accurately fit the starting position of the semi-annual cycle fluctuation. The correlation coefficient corresponding to b2 is at a moderate level. b2 reflects the phase of the annual cycle fluctuation of ZTD. Similar to a2, the annual cycle change is mainly related to the seasonal change. Latitude, longitude and elevation can reflect the climate zone and seasonal characteristics of the region to a certain extent. However, the phase of the annual cycle fluctuation is also affected by a variety of factors such as the global climate system and the distribution of land and sea. The correlation coefficient for b3 reached 0.95, which indicates that although the phase of the diurnal cycle fluctuation is affected by a variety of high-frequency changing factors such as diurnal temperature range, humidity and solar radiation, it is relatively easy to predict.

[0033] Therefore, as can be seen from the above, relying solely on the initial ZTD prediction model to predict the ZTD of an unknown observation station based on its spatiotemporal characteristics without acquiring additional meteorological data will inevitably result in bias. Therefore, this application further designs a residual improvement model to correct the initial prediction results, as described in step S4.

[0034] S4. Residual Improvement and Optimization: Based on the trained initial ZTD prediction model, seven undetermined parameters of the empirical function model of each observation station in the training dataset are obtained. Then, the initial predicted ZTD of each observation station is obtained through the empirical function model of each observation station. The residual prediction value is obtained by subtracting the initial predicted ZTD from the actual ZTD. A residual improvement model is constructed and trained with the spatiotemporal characteristics of the observation station (i.e., latitude, longitude, elevation, annual day, hour) and the initial predicted ZTD as inputs and the residual prediction value as outputs. As a preferred embodiment, in this example, the residual improvement model is as follows: Figure 5 As shown, the residual correction model employs a backpropagation (BP) neural network, using the latitude, longitude, elevation, annual day, hour, and ZTD initial prediction values ​​of the observation station as six-dimensional input features, and outputting residual prediction values. The residual correction model has six nodes in the input layer, one node in the output layer, and the number of nodes in the hidden layer is determined by five-fold cross-validation. The activation function between the input and hidden layers is the tansig function, and the activation function between the hidden and output layers is the purelin function. The Levenberg-Marquardt algorithm is used for training. After each training iteration, the mean squared error (MSE) between the output residual prediction value and the training label is calculated and recorded, and three training convergence conditions are set: the calculated MSE is less than a preset threshold, the maximum number of iterations is reached, and the prediction accuracy does not improve after multiple consecutive training iterations. Training is considered complete when any one of these convergence conditions is met.

[0035] S5. Final Prediction: (e.g., ...) Figure 6 As shown, for an unknown observation station in the current region, based on the empirical function model and the trained ZTD initial model (in... Figure 6 The initial predicted ZTD (Zero Tolerance) is obtained by fusion model (represented as fusion model). The residual sequence is obtained by improving the model through residuals. The sum of the initial predicted ZTD and the residual predicted values ​​is the final predicted ZTD for the unknown observation station. The RMSE and MAE distributions of the fusion model before and after residual correction are shown. The correction does indeed optimize the prediction accuracy. The peak values ​​of RMSE and MAE are reduced after correction, which proves that residual correction effectively optimizes the overall prediction model.

[0036] S6. Performance Verification: Using error quantification evaluation metrics, the model is compared and verified with existing tropospheric delay models on independent test datasets to evaluate the overall model's accuracy and regional adaptability. The error quantification evaluation metrics include: root mean square error (RMSE) and mean absolute error (MAE) between the final predicted ZTD and the actual ZTD of the observation stations in the independent test set.

[0037] In this embodiment, to investigate the modeling and prediction effect of the proposed three-level fusion model of "empirical model fitting - neural network parameter prediction - residual correction optimization" on tropospheric zenith delay, three existing models (CTrop model, GPT3 model and SHAtrop model) were introduced to perform ZTD prediction based on data from independent test sets. The comparison results of the evaluation indicators are shown in Table 3 below: Table 3 Comparison of evaluation indicators for the four models

[0038] As shown in Table 3, the overall prediction accuracy of the model constructed by the method proposed in this invention is comprehensively superior to the three existing traditional models. Overall, the RMSE is reduced by 7.5%, 1.0%, and 1.5% compared to the CTrop, GPT3, and SHAtrop models, respectively, while the MAE is reduced by 8.3%, 1.5%, and 1.8% compared to these three models, respectively. Furthermore, to verify the adaptability of this model for tropospheric delay modeling in the Chinese region, the accuracy results are compared in both spatial and temporal dimensions, verifying the superiority of the proposed method under different spatiotemporal conditions.

[0039] The comparison of model prediction accuracy under different dimensions is shown in Table 4 below: Table 4 Comparison of prediction accuracy at different latitudes

[0040] As shown in Table 4, in the low-latitude range (20°-30°), the RMSE and MAE of the prediction model constructed by this method are lower than those of the CTrop and GPT3 models, showing a significant reduction compared to the CTrop model. RMSE is reduced by 24.8%, and MAE by 26.4%, significantly optimizing the high error problem of traditional empirical models in the low-latitude range. In the mid-latitude range (30°-40°), the RMSE and MAE of the prediction model constructed by this method are lower than those of the other three models. Compared to the CTrop, GPT3, and SHAtrop models, RMSE is reduced by 0.9%, 0.2%, and 0.8%, respectively, and MAE is reduced by 1.7%, 0.1%, and 0.9%, respectively, indicating a slight improvement in the mid-latitude range correction model. In the high-latitude range (40°-50°), the prediction accuracy of the prediction model constructed by this method is also the best. Compared with the CTrop model, GPT3 model and SHAtrop model, its RMSE is reduced by 4.7%, 4.7% and 9.6% respectively, and its MAE is reduced by 3.5%, 4.0% and 9.1% respectively. These results show that the corrected model outperforms other models in high-latitude regions (such as North China and Northeast China). By learning the potential patterns through neural networks and correcting through residuals, it can effectively correct the prediction errors of traditional models, demonstrating the superiority and universality of this method in different latitude regions.

[0041] This embodiment also compared the prediction accuracy at different elevations, as shown in Table 5 below: Table 5 Comparison of prediction accuracy at different elevations

[0042] As shown in Table 5, at observation stations with elevations below 1000m, the prediction model constructed by this method outperforms other models in terms of ZTD prediction accuracy. Compared to the CTrop, GPT3, and SHAtrop models, its RMSE is reduced by 2.2%, 1.4%, and 4.0%, respectively, and its MAE is reduced by 1.5%, 1.7%, and 3.6%, respectively. Within the 1000-2000m range, the prediction accuracy of the model constructed by this method is basically comparable to that of the GPT3 and SHAtrop models, but significantly higher than that of the CTrop model. With the elevation increase, the RMSE and MAE decreased by 13.6% and 14.5% respectively, reaching the best accuracy of other traditional models within this elevation range. However, at elevations above 2000m, the prediction model constructed by this method achieved similar accuracy to the SHAtrop model. While the RMSE of this model was slightly higher, its MAE was slightly lower, still showing a significant advantage over the other two models. Compared to the CTrop and GPT3 models, its RMSE decreased by 6.9% and 1.7% respectively, and its MAE decreased by 10.9% and 3.9% respectively. Looking at different elevations, the prediction performance of this method is comprehensively superior to the CTrop and GPT3 models, and superior to the SHAtrop model in low-altitude areas. In other elevation ranges, it is on par with the SHAtrop model, demonstrating the superiority and universality of this method across different elevation regions.

[0043] Furthermore, for regions with distinct seasonal climate characteristics (such as my country's humid summers, dry winters, and spring-autumn monsoon transitions), these factors lead to significant changes in atmospheric parameters (water vapor and temperature), directly affecting the spatiotemporal distribution of tropospheric zenith delay. Therefore, this embodiment divides the time series by season, specifically spring (March-May), summer (June-August), autumn (September-November), and winter (December-February), to explore the difference in prediction accuracy between this method and existing methods in different seasons, as shown in Table 6 below: Table 6 Comparison of prediction accuracy in different seasons

[0044] As shown in Table 6, in spring (March-May), the prediction accuracy of the model constructed by this method is comprehensively superior to other models. Compared with the CTrop, GPT3, and SHAtrop models, the RMSE is reduced by 9.3%, 2.6%, and 4.5%, respectively, and the MAE is reduced by 6.2%, 0.6%, and 2.6%, respectively. This indicates that the prediction effect of this method is better in spring, significantly superior to the CTrop and SHAtrop models, reflecting that the prediction model is more adaptable to the water vapor changes during the spring monsoon transition period in China. In summer, the prediction accuracy of the model constructed by this method is also comprehensively superior to other models. Compared with the CTrop, GPT3, and SHAtrop models, its RMSE is reduced by 4.1%, 3.2%, and 1.1%, respectively, and the MAE is reduced by [missing data]. The accuracy rates were 5.8%, 3.7%, and 3.5% lower, indicating that this method has an absolute advantage over other models in the summer season when water vapor changes are large and the delay fluctuations are severe. Traditional models are difficult to fit the nonlinear growth in summer. In autumn (September-November), the prediction accuracy of the prediction model constructed by this method is slightly lower than that of the GPT3 model and is basically equivalent to that of the SHAtrop model, but still higher than that of the CTrop model. Its RMSE and MAE are reduced by 3.0% and 4.9% respectively compared to the CTrop model. In winter (December-February), the prediction accuracy of the prediction model constructed by this method is better than that of the CTrop and GPT3 models and is comparable to that of the SHAtrop model. Compared with the former two, its RMSE is reduced by 13.3% and 2.3% respectively, and its MAE is reduced by 15.3% and 5.1% respectively. Looking at the comparison across the four seasons, this method maintains an absolute advantage over the CTrop model throughout the year. Compared to the GPT3 model, the corrected model shows a significant advantage in spring, summer, and winter, while its accuracy is slightly lower in autumn. Compared to the SHAtrop model, the corrected model is optimal in spring and summer, while the accuracy of the two models is roughly equivalent in autumn and winter. Therefore, overall, this method still demonstrates its superiority and universality across different seasons.

[0045] Furthermore, this embodiment discloses a regional tropospheric delay modeling system based on GNSS observation data and a fusion model, which specifically includes: The data management module is used to acquire raw GNSS observation data of the target area for preprocessing, divide the training dataset and the test dataset, and obtain the actual tropospheric zenith delay (ZTD) of each observation station. The empirical modeling and solution module is used to establish an empirical function model to fit the tropospheric zenith delay (ZTD) of each observation station, and to solve the model parameters of the empirical function model corresponding to the actual ZTD based on the training dataset, which are then used as training labels. The initial prediction module is used to construct an initial ZTD prediction model that takes the spatial coordinates of the observation station as input and solves the parameters of the empirical function model as the objective. After training based on the training labels and training dataset, it performs initial ZTD predictions for the observation stations in the current area. The residual correction module is used to build and train a residual improvement model that takes the spatiotemporal information of the observation station and the initial predicted ZTD as inputs and outputs the residual prediction value between the initial predicted ZTD and the actual ZTD. The online prediction module receives the spatial coordinates of the observation station to be observed, and outputs a high-precision final prediction ZTD by calling the trained neural network model and the residual improvement model. The model validation and evaluation module is used to call the test dataset to evaluate and compare the accuracy of the final ZTD prediction value.

[0046] In summary, this invention, based on the original fusion model composed of the empirical function model and the initial ZTD model, optimizes through residual correction. This not only preserves the empirical function model's ability to fit the potential mapping between the "spatial-temporal characteristics of the observation station and the tropospheric zenith delay," but also solves the prediction bias caused by limitations imposed by spatiotemporal characteristics. Compared to other existing models, this method maintains a leading level of prediction performance under different spatiotemporal conditions, overcoming the shortcomings of existing technologies in tropospheric delay modeling, such as insufficient accuracy and poor spatiotemporal adaptability. Therefore, this invention provides a fusion modeling method based on empirical models and multi-stage machine learning. Through a three-level architecture of "empirical model fitting—neural network parameter prediction—residual correction optimization," it achieves high-precision ZTD prediction for target areas with complex landforms and monsoon climate characteristics.

[0047] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the present invention. No reference numerals in the claims should be construed as limiting the scope of the claims.

[0048] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A method for regional tropospheric delay modeling based on GNSS observation data and a fusion model, characterized in that, Specifically, the following steps are included: S1. Obtain long-time raw observation data from GNSS continuously operating observation stations that uniformly cover the target area, process the data to obtain the actual ZTD of each observation station, and divide it into training dataset and test dataset. The prior dataset is composed of observation data from observation stations whose data integrity meets the preset requirements in the training dataset. S2. An empirical function model for fitting the tropospheric zenith delay (ZTD) of each observation station is constructed in the form of a periodic function. Based on the prior data set, n undetermined parameters of the empirical function model of the known observation station are calculated. S3. Using the longitude, latitude, and elevation of each observation station as input and the undetermined parameters of the empirical function as the output target, construct an initial ZTD prediction model consisting of n independent BP neural networks. Each BP neural network is responsible for predicting one undetermined parameter. Use the n undetermined parameters of the known empirical function model of the observation station as training labels to train the initial ZTD prediction model. S4. Based on the trained initial ZTD prediction model, obtain n undetermined parameters of the empirical function model of each observation station in the training dataset. Then, obtain the initial ZTD prediction of each observation station through the empirical function model of each observation station. Subtract the initial ZTD prediction from the actual ZTD to obtain the residual prediction value. Construct and train a residual improvement model with the spatiotemporal characteristics of the observation station, the initial ZTD prediction as input, and the residual prediction value as output. S5. For an unknown observation station in the current region, the initial predicted ZTD is obtained based on the empirical function model and the trained ZTD initial model. The residual sequence is obtained through the residual improvement model. The sum of the initial predicted ZTD and the residual predicted value is used as the final predicted ZTD of the unknown observation station. S6. Using error quantification evaluation metrics, the model is compared and verified with existing tropospheric delay models on independent test datasets to evaluate the overall model's accuracy and regional adaptability.

2. The regional tropospheric delay modeling method based on GNSS observation data and fusion model according to claim 1, characterized in that, In step S1, the process of acquiring long-time-series raw observation data from GNSS continuously operating observation stations that uniformly cover the target area and processing it to obtain the actual ZTD of each observation station specifically involves: Based on the observation data of the landform network, in the target area with complex landforms and monsoon climate characteristics, GNSS observation data from N continuously operating observation stations are collected as the raw observation data. The original observation data of each observation station were sequentially processed by cycle slip repair, gross error removal, and spatiotemporal alignment to obtain the actual ZTD of each observation station.

3. The regional tropospheric delay modeling method based on GNSS observation data and fusion model according to claim 1, characterized in that, In step S2, the empirical function model includes: a fundamental fixed component in the tropospheric zenith delay that does not change periodically with time, and a periodic fluctuation component in the tropospheric zenith delay that changes with different time scales, expressed by the formula: ； Where c is a constant term, representing the fundamental fixed component in the tropospheric delay that does not change periodically with time; , , These represent the periodic fluctuation components in the tropospheric zenith delay that vary with the semi-annual, annual, and daily cycles, respectively. , , These are the amplitude coefficients of each periodic term, representing the magnitude of the periodic fluctuations; , , These are the phase parameters for each periodic term, reflecting the time position of the fluctuation; doy is the yearly cumulative day. , , , , , c represents the parameters of the empirical function model to be solved.

4. The regional tropospheric delay modeling method based on GNSS observation data and fusion model according to claim 1, characterized in that, In step S2, the specific steps of calculating the n undetermined parameters of the known observation station empirical function model based on the prior dataset are as follows: The nonlinear least squares method is used to fit the n undetermined parameters of the empirical function model. The sum of squared residuals between the actual ZTD of each observation station and the ZTD calculated by the model is used as the fitting objective function, and the undetermined parameters are optimized by the Levenberg-Marquardt iterative algorithm.

5. The regional tropospheric delay modeling method based on GNSS observation data and fusion model according to claim 1, characterized in that, Step S3 specifically includes: The latitude (B), longitude (L), and elevation (H) of a single observation station are used as inputs. Each BP neural network consists of an input layer, a hidden layer, and an output layer. The input layer has three nodes, corresponding to the latitude, longitude, and elevation components, while the output layer has one node, corresponding to the model parameters that each BP neural network is responsible for predicting. The number of nodes in the hidden layer is determined through five-fold cross-validation. The activation function between the input and hidden layers is the tansig function, and the activation function between the hidden and output layers is the purelin function. After normalizing the three components of the input features, the Levenberg-Marquardt algorithm is used for training. After each training iteration, the mean squared error (MSE) between the output model parameters and the training labels is calculated and recorded. Three training convergence conditions are set: the calculated MSE is less than a preset threshold, the maximum number of iterations is reached, and the prediction accuracy does not improve after multiple consecutive training iterations. Training is considered complete when any one of the convergence conditions is met.

6. The regional tropospheric delay modeling method based on GNSS observation data and fusion model according to claim 1, characterized in that, The residual correction model described in step S4 employs a backpropagation (BP) neural network, using the latitude, longitude, elevation, annual day, hour, and ZTD initial prediction values ​​of the observation station as six-dimensional input features, and outputting residual prediction values. The residual correction model has six nodes in the input layer, one node in the output layer, and the number of nodes in the hidden layer is determined through five-fold cross-validation. The activation function between the input and hidden layers is the tansig function, and the activation function between the hidden and output layers is the purelin function. The Levenberg-Marquardt algorithm is used for training. After each training iteration, the mean squared error (MSE) between the output residual prediction value and the training label is calculated and recorded, and three training convergence conditions are set: the calculated MSE is less than a preset threshold, the maximum number of iterations is reached, and the prediction accuracy does not improve after multiple consecutive training iterations. Training is considered complete when any one of these convergence conditions is met.

7. The regional tropospheric delay modeling method based on GNSS observation data and fusion model according to claim 1, characterized in that, The error quantification evaluation indicators include: the root mean square error (RMSE) and the mean absolute error (MAE) between the final predicted ZTD and the actual ZTD of the independent test set observation stations.

8. A regional tropospheric delay modeling system based on GNSS observation data and a fusion model, characterized in that, Specifically, it includes: The data management module is used to acquire raw GNSS observation data of the target area for preprocessing, divide the training dataset and the test dataset, and obtain the actual tropospheric zenith delay (ZTD) of each observation station. The empirical modeling and solution module is used to establish an empirical function model to fit the tropospheric zenith delay (ZTD) of each observation station, and to solve the model parameters of the empirical function model corresponding to the actual ZTD based on the training dataset, which are then used as training labels. The initial prediction module is used to construct an initial ZTD prediction model that takes the spatial coordinates of the observation station as input and solves the parameters of the empirical function model as the objective. After training based on the training labels and training dataset, it performs initial ZTD predictions for the observation stations in the current area. The residual correction module is used to build and train a residual improvement model that takes the spatiotemporal information of the observation station and the initial predicted ZTD as inputs and outputs the residual prediction value between the initial predicted ZTD and the actual ZTD. The online prediction module receives the spatial coordinates of the observation station to be observed, and outputs a high-precision final prediction ZTD by calling the trained neural network model and the residual improvement model. The model validation and evaluation module is used to call the test dataset to evaluate and compare the accuracy of the final ZTD prediction value.

Images