zenith wet delay prediction method based on periodic space-time modeling and cascade residual correction

Abstract

The application discloses a zenith wet delay prediction method based on period space-time modeling and cascade residual correction. Firstly, the ZWD observation sequence of a GNSS reference station is used to extract the period change characteristics of ZWD through spectrum analysis, and a spatial variation model of the period parameter is established by using a spherical cap harmonic analysis function to obtain the ZWD period model value at any position and time. Then, the residual between the period model value and the actual ZWD is taken as a prediction object, a data set containing the ZWD residual, time information, spatial coordinates and meteorological parameters is constructed, and the CNN-BiLSTM-Attention network is inputted for first-stage residual prediction, and then the XGBoost model is used to correct the remaining residual in the second stage. Finally, the period model value and the two-stage residual prediction results are combined to obtain the ZWD prediction value at the target position and target time. The application can improve the ZWD modeling precision and space-time extrapolation capability while retaining the physical prior information, and can be used for regional high-precision GNSS troposphere delay modeling and correction.

Description

Technical Field

[0001] This invention belongs to the field of space geodesy and satellite navigation positioning technology, specifically involving a zenith wet delay prediction method based on periodic spatiotemporal modeling and cascaded residual correction. Background Technology

[0002] When GNSS signals propagate from a satellite to a ground receiver, they must pass through the Earth's neutral atmosphere; the resulting propagation delay is called tropospheric delay. At the zenith, the total tropospheric delay is approximately 2 to 2.5 meters. As the satellite's elevation angle decreases, this delay is amplified according to a mapping function, reaching over ten meters in low-elevation-angle observations. Tropospheric delay is one of the main error sources limiting convergence speed and elevation accuracy in current Precise Point Positioning (PPP) and Real-Time Kinematic (RTK) systems.

[0003] The tropospheric total delay (ZTD) is typically decomposed into two parts: the tropospheric hydrostatic delay (ZHD) and the tropospheric wet tropospheric delay (ZWD). The former, generated by dry gases, accounts for about 90% of the total, and its value is approximately proportional to surface air pressure, which can be corrected to the millimeter level using surface air pressure observations. The latter, generated by water vapor in the atmosphere, accounts for only about 10% of the total, with a value of approximately 0 to 40 centimeters. However, water vapor is extremely unevenly distributed in time and space and exhibits strong nonlinear variations from interannual and seasonal levels down to the hourly level, making it difficult to characterize with simple functions. Therefore, it has become the main source of error and technical challenge in tropospheric modeling.

[0004] Unlike ionospheric delay, tropospheric delay exhibits non-dispersive characteristics for navigation signals, meaning the delay is independent of signal frequency. Therefore, it cannot be eliminated through a linear combination of multi-frequency observations and must be corrected using a model. This characteristic dictates that improving the accuracy of the tropospheric wet delay model is essential for enhancing positioning accuracy. Furthermore, a definite conversion relationship exists between ZWD and Precipitable Water Vapor (PWV). A high-precision wet delay model can simultaneously provide observational constraints for ground-based water vapor inversion, short-term heavy precipitation nowcasting, and numerical weather prediction data assimilation, thus possessing dual value for both navigation and meteorology.

[0005] There are three main types of existing tropospheric wet delay correction techniques, each with its own significant shortcomings.

[0006] The first category consists of empirical models based on historical climate data, such as the Saastamoinen model, the Hopfield model, and the Global Pressure and Temperature (GPT) series of models. These models only require station coordinates and annual days for calculation, without the need for measured meteorological parameters, making them convenient to use. However, they are essentially climate mean fields, reflecting only the wet-delayed multi-year average seasonal variation and failing to capture actual weather processes. Their errors increase significantly when there are drastic changes in water vapor.

[0007] The second category is based on numerical weather prediction (NWP) models, such as methods that use reanalysis data or forecast field integration to obtain the delay. These methods have high accuracy, but high-precision reanalysis products have a release delay of several days, which cannot meet the needs of real-time applications; while the accuracy of real-time forecast fields is significantly lower than that of reanalysis products, and gridded products lack vertical representativeness in areas with complex terrain.

[0008] The third category is interpolation models based on regional benchmark station networks. These methods can achieve millimeter-level accuracy in densely populated areas, but their accuracy is highly dependent on the station network density. In sparsely populated areas such as oceans, deserts, and mountains, spatial interpolation errors increase rapidly. At the same time, most methods use uniform interpolation weights that are independent of location, failing to take into account the physical law that wet delay decays rapidly with altitude.

[0009] In recent years, with the development of machine learning methods, technical solutions have emerged that directly predict ZWD or ZTD using models such as Long Short-Term Memory (LSTM) networks, Random Forest (RF), and Gradient Boosting Decision Tree (GBDT). These purely data-driven solutions achieve high accuracy within the spatiotemporal range covered by the training samples, but they suffer from three common drawbacks: First, they lack large-scale physical constraints, leading to a significant drop in accuracy and insufficient robustness when extrapolating to stations not used in training or to unseen extreme weather events. Second, when using the raw wet delay as the learning objective, the model must simultaneously fit annual-cycle backgrounds and hourly-level disturbances of vastly different magnitudes, limiting its cross-scale fitting ability. Third, if historical and future information are not strictly distinguished when constructing lag features and rolling statistical features, label information leakage can easily occur, resulting in inflated offline evaluation metrics and insufficient actual forecasting capabilities.

[0010] Furthermore, in regional spatial field modeling, ordinary spherical harmonic functions are orthogonal basis function systems across the entire sphere. When used for modeling finite regions, they must be truncated to higher orders to characterize regional details, leading to parameter inflation and ill-conditioned design matrices. Spherical Cap Harmonic Analysis (SCHA) reconstructs a basis function system satisfying boundary conditions on a spherical cap domain, achieving better regional fitting capabilities and numerical conditions with fewer parameters. This method has previously been primarily applied to modeling regional geomagnetic fields, ionospheric electron content, and regional gravity fields. In recent years, its application to tropospheric ZTD spatial representation has also been reported. However, existing reports all directly use the delay quantity itself or its correction as the object of spherical cap harmonic expansion; no application to spherical cap harmonic coefficients in wet-delay time harmonic coefficient fields has been observed, nor has a technical solution combining it with exponential height scaling factors and two-level residual learning cascades been found.

[0011] In summary, how to fully utilize deep learning and ensemble learning to characterize high-frequency perturbations of wet latency while preserving physical priors and robustness to large-scale constraints, and at the same time strictly avoid information leakage and provide reliability criteria for spatial extrapolation, is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0012] To address the aforementioned technical problems, the present invention aims to provide a method for predicting tropospheric zenith wet delay based on periodic spatiotemporal modeling and cascaded residual correction. The core concept of this invention is: first, an interpretable physical prior model is used to remove the large-scale, low-frequency, and repeatable deterministic components from the wet delay; then, the learning task is reduced to residuals with smaller magnitudes and more stable statistical properties, which are then approximated step-by-step by a two-stage machine learning model. This simultaneously achieves the extrapolation robustness of the physical model and the high-frequency characterization capability of the data-driven model.

[0013] Technical solution The zenith wet delay prediction method based on periodic spatiotemporal modeling and cascaded residual correction includes the following steps: Step S1: Obtain GNSS tropospheric zenith total delay (ZTD) products aligned by station and time, and reanalysis meteorological data covering the same latitude and longitude range; read the station latitude, longitude, geodetic height, ZTD of each epoch, as well as surface air pressure, air temperature, and dew point temperature; interpolate the meteorological field to the station location; calculate the pressure-reduced zenith static delay (ZHD); subtract ZHD from ZTD to obtain ZWD; after quality control of ZWD and input fields, output the ZWD observation sequence and corresponding exogenous variables at the unified sampling time, as input for steps S2 and S4.

[0014] Step S2: Using the ZWD observation sequences of each station output in Step S1 as input, first determine the annual cycle and semi-annual cycle as the main cycles through spectrum analysis, and then establish a five-parameter harmonic regression model with a fixed cycle; calculate the basic angular frequency, and solve the constant term, annual cycle cosine term, annual cycle sine term, semi-annual cycle cosine term, and semi-annual cycle sine term for each station according to the least squares criterion, and output the five harmonic coefficients of each station.

[0015] Step S3: Using the five harmonic coefficients output from Step S2, along with the station coordinates and geodetic height, the center of the spherical cap is determined by the mean latitude and longitude of the stations involved in the fitting. The half-angle of the spherical cap is determined by the spherical angular distance from the center of the spherical cap to the farthest station. A SCHA (Spherical-Cap Harmonic Analysis) basis function that satisfies the Neumann boundary conditions is constructed. An exponential height scaling factor is superimposed, and the ridge regression regularization coefficient is selected through spatial block cross-validation. The baseline ZWD prediction value at any location and at any time within the modeling area is output.

[0016] Step S4: Subtract the ZWD observation value output in Step S1 from the baseline ZWD prediction value output in Step S3 to obtain the first-level residual objective; sort the exogenous variables and residual history sequences by station and time and apply causal lag, construct an input feature set that does not contain measured ZWD, ZTD, ZHD and their direct conversion quantities, where each time point consists of a 24-dimensional vector composed of 8 basic features, 9 water vapor pressure derived features and 7 residual history features, and output a fixed-length history sequence for use in Step S5.

[0017] Step S5: Input the 16 consecutive 6-hour time-series of each station output from Step S4, with each time-series consisting of 24 dimensions, into the first-level residual network composed of a one-dimensional convolutional neural network (CNN), a bidirectional long short-term memory network (BiLSTM), an additive attention convergence layer, and a fully connected output layer. After extracting local changes and bidirectional time-series information, output the first-level residual prediction value.

[0018] Step S6: Concatenate the input feature set output in step S4 with the first-level residual prediction value output in step S5 as the input to the gradient boosting decision tree implemented by XGBoost. Subtract the first-level residual prediction value from the first-level residual objective to obtain the second-level learning objective. Train and output the second-level residual prediction value.

[0019] Step S7: Add the baseline ZWD prediction value output in step S3, the first-level residual prediction value output in step S5, and the second-level residual prediction value output in step S6 to obtain the final ZWD prediction value. Based on the relationship between the point to be predicted and the convex hull of the set of stations to be fitted, and whether the spherical distance to the nearest fitted station exceeds the distance threshold, output the extrapolation reliability indicator of the prediction value.

[0020] Step S1 includes: Step S11: Obtain the tropospheric delay product in the standard exchange format (SINEX_TRO), parse its station information block and tropospheric solution block, and obtain the latitude, longitude, geodetic height and ZTD of each station for each epoch. For products that only provide geocentric rectangular coordinates, geodetic coordinates are calculated using an iterative method; for cases where multiple station codes exist for the same physical station, the first four digits of the station code are used to merge them into the same physical station. Step S12: For tropospheric products whose sampling time does not fall on a whole hour, take the sample closest to the target whole hour within the preset time tolerance range, and align its timestamp to the target whole hour so that data from different sources have a unified sampling time. Step S13: Clip the ERA5 reanalysis data according to the latitude and longitude range of the modeling area to obtain the surface analysis field including surface air pressure, two-meter air temperature, two-meter dew point temperature, and total water vapor in the whole layer, as well as the invariant field including the land-sea mask and surface potential. Divide the surface potential by the standard gravitational acceleration to obtain the reanalysis terrain height. Interpolate the above fields to the plane position of each reference station using bilinear interpolation. Step S14: Calculate the vapor pressure based on the two-meter dew point temperature, as follows: In the formula, Here, is the dew point temperature in degrees Celsius, and e is the vapor pressure in hectopascals. Step S15: Reduce the reanalysis surface pressure from the reanalysis topographic height to the station's geodetic height, as follows: In the formula, For further analysis of surface air pressure, h is the geodetic height of the station. To further analyze terrain elevation, Standard gravitational acceleration, Let be the gas constant for dry air, T be the air temperature, and q be the specific humidity. It is a false fever. To reduce the air pressure to the station altitude, virtual temperature is used instead of air temperature for the reduction, so that the air pressure reduction takes into account the influence of water vapor on air density. Step S16: Calculate ZHD from the reduced station pressure and derive the ZWD observation sequence for each station: In the formula, φ is the latitude of the measuring station. The geodetic height of the station is measured in kilometers. The ZTD calculated for GNSS; the resulting ZWD serves as the modeling target and evaluation benchmark for subsequent steps. Step S17: Apply quality control to the ZWD observation sequence, including: removing non-physical records with ZWD not greater than zero, removing records where the Precipitable Water Vapor (PWV) converted from ZWD exceeds the preset physical range, removing records with missing key fields, correcting the coordinates of stations found to have coordinate errors according to their original correct codes, and removing stations with sample sizes below the preset threshold during the evaluation period, thus obtaining the dataset used for modeling.

[0021] After this processing, the ZWD observation sequence output in step S1 serves as the time fitting target for step S2 and the evaluation benchmark for subsequent residual learning.

[0022] The exogenous variables include: station latitude, longitude, and geodetic height; air temperature converted from ERA5 2-meter air temperature; water vapor pressure calculated from ERA5 2-meter dew point temperature; surface air pressure corrected by station elevation; and sine and cosine time characteristics constructed from annual accumulated days.

[0023] In step S2, the fixed-period five-parameter harmonic regression model is as follows: In the formula, t represents the accumulated days of a year. The fundamental angular frequency, , , , , These are the five coefficients of the five-dimensional harmonic coefficients; The calculation formula is as follows: The annual cycle and the semi-annual cycle are determined as the main cycle components by the spectral analysis of the ZWD sequence; the five coefficients in equation (5) are solved according to the least squares criterion, and the missing epochs are not filled in the modeling target. Only the effective observation times are used to form the design matrix and complete the solution.

[0024] It should be noted that the harmonic analysis here uses periodic LS regression instead of Discrete Fourier Transform (DFT) inversion. The advantage of this method is that it is equally applicable to measured sequences with missing data and uneven sampling, and the five coefficients of the resulting five-dimensional harmonic coefficients are... , , , , It has a clear physical meaning, corresponding to the annual average level, the annual period cosine and sine terms, and the semi-annual period cosine and sine terms respectively; the number of rows of the design matrix G is equal to the number of effective observations of the station, the number of columns is 5, the observation vector y is the ZWD observation value at the same time, and the coefficient vector c is the five-dimensional harmonic coefficient.

[0025] The formula for solving the problem using the least squares method is as follows: Where y is the wet delay observation vector obtained in step S1, with dimension N; G is the vector composed of the values ​​of the five time basis functions shown in equation (6) at each observation time. Design matrix, where c is the vector of five-dimensional harmonic coefficients to be determined; when there are missing samples, only the valid times are used to construct G and y, and the missing values ​​are not filled into the target variable.

[0026] Step S3 includes: Step S31: After processing in step S2, the time variation of the wet delay has been compressed into five coefficients per station. Step S3 treats these five coefficients as spatial fields defined on the spherical cap domain and expands them using SCHA basis functions. Using the mean latitude and longitude of all participating stations as the center of the spherical cap, and the spherical angular distance from the center to the farthest station as the half-angle of the spherical cap, the geographical coordinates of each station are rotated into the co-latitude and azimuth of the spherical cap with the center as the pole. Rotate the geographic coordinates of the station or the point to be predicted into spherical coordinates: in, and Input the geographical latitude and longitude of the point. and The geographical latitude and longitude of the crown of the globe. The coaxial latitude of the spherical cap relative to the coronal center is the input point; Step S32: For each azimuth order m, using the Neumann boundary condition (i.e., the first derivative of the associated Legendre function along the colatitude direction is zero at the spherical cap boundary), solve for a set of non-integer degrees, and retain only the solution family where the difference between the degree and the order is even; as follows: in Let be the associated Legendre function of degree n and orientation order m. Given the spherical cap's half-angle, the boundary condition states that the first derivative of the function along the colatitude direction is... The value is zero; In one embodiment of the present invention, 0≤m≤M and 0≤n≤K are selected, and only the solution family where nm is an even number is retained.

[0027] Step S33: Using the Schmidt quasi-normalized associated Legendre function as the radial basis, expand the j-th harmonic coefficient into a spatially continuous field according to the following formula: in, and Let be the harmonic expansion coefficient of the spherical cap to be determined; The degree is a non-integer number determined by the boundary conditions; H is the reference height in meters; H is the exponential height scaling scale in meters; m is the azimuth order; k is the number index; M and K are the azimuth order and the upper limit of the number, respectively. Let be the Schmidt quasi-normalized associated Legendre function, and let represent the Schmidt quasi-normalization; in equation (9) is the associated Legendre function used for boundary conditions.

[0028] Step S34: Solve for the harmonic coefficients of the spherical cap using Ridge regression: In the formula, X is the design matrix composed of the basis functions described in step S33; Let be the spherical cap harmonic expansion coefficient vector of the spatial field of the j-th time harmonic coefficient; This is the observation vector of the harmonic coefficients of the j-th station obtained in step S2; The ridge regression regularization parameter corresponding to the j-th coefficient field is selected using spatial block cross-validation.

[0029] After the above spatial modeling is completed, the five coefficient spatial fields output in step S3 can take values ​​at any point to be predicted, and the baseline wet delay at any location and time can be reconstructed according to equation (11): in, The j-th time harmonic coefficient is the spatial field value at the spherical cap coordinates and the ground height. The coefficient is derived from step S3 for each station. Perform SCHA expansion and ridge regression; Let t be the geodetic height of the point to be predicted, and t be the annual day or equivalent time code of the target time. and Given the geographical latitude and longitude of the point to be predicted. and Given the coaxial latitude and azimuth of its spherical crown. The time basis function multiplied by the j-th coefficient in equation (5) is specifically: ,in .

[0030] Furthermore, the regularization coefficient in step S34 The selection process employs spatial block cross-validation: first, the stations participating in the fitting are sorted lexicographically by latitude and longitude; then, the sorted stations are divided into several non-overlapping spatial blocks. One spatial block is reserved as a validation block at a time, and the remaining blocks participate in the fitting. The block whose sum of the root mean square error (RMSE) of the coefficients of all validation blocks is selected. The spatial partitioning makes the stations within the same validation block geographically adjacent to each other, so that cross-validation measures the spatial extrapolation ability of the model rather than its random interpolation ability.

[0031] Step S4: Construction of Residual Target and Historical Input Information Based on measured wet delay Compared with baseline wet delay difference As a primary learning objective: Since the baseline model has already characterized the annual and semi-annual cycle variations in the wet delay, the residual weakens the main periodic components to some extent, and its variation amplitude is significantly reduced compared to the original wet delay. The stationarity of the time series is enhanced, thereby reducing the learning difficulty of the model and facilitating network training and convergence.

[0032] To avoid information leakage, step S4 imposes constraints on all derived features. The eight basic features at each time point are latitude, longitude, geodetic height, air temperature, station air pressure, water vapor pressure, and sine and cosine codes for annual accumulated days; the nine derived features of water vapor pressure are the first, second, and fourth order differences of water vapor pressure (three features), as well as the rolling mean and rolling standard deviation of time windows at times 4, 8, and 16 (six features); the seven residual history features are the first to third order lags of the first-level residuals (three features), the first to second order differences (two features), and the rolling mean and rolling standard deviation of the four time windows (two features), totaling 24 dimensions.

[0033] All records are first sorted by station and time. Lag terms, difference terms and rolling statistics are calculated by shifting the entire group of the same station by one time. When the time interval between adjacent records exceeds 6 hours, the continuous sequence is cut off, the historical window is reset for the new sequence, and the historical value of the previous sequence is not used. The measured ZWD, ZTD, ZHD and their direct conversion quantities are excluded from the output.

[0034] Step S5: First-level residual learning: CNN-BiLSTM-Attention residual prediction.

[0035] Step S5 takes a 24-dimensional historical feature sequence of 16 consecutive 6-hour time points from the same station as input, covering the 4 days prior to the target time; information after the current time, the measured ZWD of the target time, and future residuals are not included in the input. The sequence is first processed by a one-dimensional convolutional neural network (CNN) layer to extract local change patterns between adjacent time points, and then fed into a bidirectional long short-term memory (BiLSTM) network.

[0036] in, Let be the 24-dimensional input vector at time t. The number of time steps covered by the convolution kernel. and Here, BN represents batch normalization, and ReLU represents the modified linear unit activation function. The CNN outputs features, which are then used as input to the subsequent BiLSTM, and the CNN output sequences at each time step are processed. Inputting the data into a BiLSTM yields the forward and reverse hidden states. The forward hidden state... With reverse hidden state spliced ​​into formula (14) : Then, additive attention is used to weight and converge the hidden states at each time step according to equation (15), so that the model automatically emphasizes the historical moments that contribute more to the current residual; where, Rate attention For normalized weights, z is the context vector, and W, b, and v are learnable parameters. The converged z is passed through a fully connected layer to output the first-level residual prediction value of equation (16): in and These are the output layer parameters.

[0037] In one embodiment of the present invention, the CNN output has 48 channels and a kernel size of 3; the BiLSTM unidirectional hidden state dimension is 96. The fully connected hidden layer has a dimension of 64 and a Dropout probability of 0.15; the optimizer uses AdamW with an initial learning rate of The weight decays to The batch size is 1024, and the gradient norm pruning threshold is 5. The training set consists of data from January 1, 2020 to December 31, 2022; the fitting interval is from January 1, 2023 to September 30, 2023, and the early stopping interval is from October 1, 2023 to December 31, 2023, with no overlap between the two intervals; the maximum number of training epochs is 50, and early stopping occurs when the loss in the early stopping interval does not improve for 8 consecutive epochs.

[0038] Step S6: Second-level residual learning: Residual residual modeling based on XGBoost.

[0039] Step S6 receives the 24-dimensional current features from step S4 and the first-level residual prediction values ​​from step S5, constructs the second-level learning objective using equation (17), and fits equation (18) using GBDT tree-by-tree addition implemented by XGBoost.

[0040] In equation (17), The objective for the residuals learned by the second-level model is equal to the measured ZWD minus the baseline ZWD and the predicted residuals from the first level. Therefore, the input dimension of the second-level model is 24 + 1 = 25 dimensions, and the output is the predicted residuals from the second level.

[0041] In the formula, η is the learning rate. For the s-th regression tree, Let S be the set of candidate regression tree functions, and S be the final number of iteration trees. This represents the second-level XGBoost residual prediction value obtained by summing the results of S GBDT trees one by one with a learning rate η.

[0042] Its regularization objective function is: In equation (19), For single-sample loss function, It is the first The second-level "true value of residual" (measured ZWD) for each sample Baseline ZWD First-level CNN-BiLSTM-Attention residual prediction These are the corresponding second-level XGBoost residual predictions, and the two constitute the loss function l( , ); Let be the number of leaf nodes in the s-th tree. The leaf node weight vector, The penalty coefficient for the number of leaf nodes. is the L2 regularization coefficient for the leaf node weights.

[0043] In one embodiment of the present invention, parameters are determined by time continuity verification and early stopping in subsequent time periods: the maximum depth is 4, which is used to limit the complexity of a single tree; =0.02, used to control the step size of each round of addition update; the minimum weight sum of child nodes is 10, used to suppress splits due to insufficient sample size; row sampling is 0.90 and column sampling is 0.85, used to reduce inter-tree correlation; =0.5、 =10.0, used to constrain the number of leaf nodes and the weight of leaf nodes; the number of early stopping iterations is 250, used to stop iterations when the error no longer improves in subsequent time periods; the final number of trees is 122.

[0044] Step S7: Synthesis of ZWD prediction results and determination of extrapolation range.

[0045] Step S7: Given the geographical latitude of the point to be predicted. ,longitude , the earth is high Given the target time t, first calculate the baseline ZWD according to step S3; then construct the 24-dimensional historical features of the same station for 16 consecutive 6-hour periods up to t according to step S4, and load the CNN-BiLSTM-Attention model to obtain... Combine the current 24-dimensional features with The vectors are concatenated into a 25-dimensional vector and input into the XGBoost model to obtain the second-level residual prediction value of the XGBoost model for the target point after training. The three terms are added together according to Equation (20) to obtain the final ZWD prediction value. Then, the planar positions of the stations involved in the fitting are used to form a convex hull, and the spherical distance from the point to be predicted to the nearest fitted station is calculated according to Equation (21).

[0046] In the formula, R is the average radius of the Earth, and Δφ and Δλ are the difference in latitude and longitude, respectively. and The latitude of two points. Spherical distance. It is used to measure the spatial proximity between the point to be predicted and the nearest fitted station and compare it with a preset distance threshold; the 95th percentile of the nearest neighbor distance distribution of the training station is used as the distance threshold; when the point to be predicted is outside the convex hull or the distance to the nearest station exceeds the threshold, it is marked as an extrapolation point and output along with the predicted value, so that the weight of the predicted value can be adjusted according to reliability for downstream positioning or water vapor inversion.

[0047] Beneficial effects Compared with the prior art, the present invention has the following beneficial effects: 1. Significantly Improved Accuracy. This invention extracts the temporal background through periodic harmonic regression, describes the spatial background using SCHA, learns historical residuals using the first-level CNN-BiLSTM-Attention network, and further corrects the remaining errors using the second-level XGBoost, achieving a step-by-step improvement in accuracy. Experimental data shows that the root mean square error (RMSE) was 36.00 mm when using only SCHA as the baseline; after adding the first-level CNN-BiLSTM-Attention network, it decreased to 16.77 mm; and after adding the second-level XGBoost network, it further decreased to 16.49 mm, a reduction of approximately 54.2% relative to the baseline. The final coefficient of determination (COP)... The correlation coefficient (CC) was 0.906, the mean absolute error (MAE) was 12.08 mm.

[0048] 2. Spatial Continuity and Global Usability. The baseline model expands the station-by-station coefficients into a continuous spatial field over a spherical cap, freeing the model from limitations on existing station locations and allowing values ​​to be taken at any location within the modeling area. Experimental results show that the station-by-station time model only covers approximately 85.3% of the test samples, while the coverage reaches 100% after applying the spherical cap harmonic coefficients, with essentially the same accuracy, achieving a substantial expansion of the coverage area.

[0049] 3. Extrapolation robustness was tested using a double-blind method. This invention employs a dual hold-out validation method, both temporally and spatially: temporally, both the baseline and the two-stage residual model are fitted using data from 2020 to 2023, with 2024 data not participating in any training; spatially, a separate set of hold-out stations is established, whose observations are not included in the baseline fitting. Test results show that the RMSE of the spatially held-out stations is 16.92 mm, while that of the newly added stations that did not participate in the training is 15.99 mm, which is essentially equivalent to the 16.45 mm of the existing stations. This indicates that this invention does not rely on the memory of specific stations, but rather learns transferable physical and statistical laws.

[0050] 4. Strictly prevent information leakage. All lag, differential, and rolling statistical features are calculated by shifting the entire calculation one time step back within the same station group. All cross-station and cross-time discontinuities are cut off, and a list of prohibited features is provided. This design ensures that offline evaluation indicators are consistent with actual forecast capabilities, avoiding the problem of inflated indicators common in similar data-driven schemes.

[0051] 5. Consistent physical modeling targets and unified caliber. This invention does not directly target the ZTD given in publicly available products. Instead, it calculates the ZHD after pressure-height reduction using virtual temperature correction and then removes it, thus obtaining a well-defined ZWD observation sequence. Compared to schemes that directly model the ZTD, the evaluation index of this invention is not diluted by the static component, which accounts for about 90% of the total, and can truly reflect the ability to characterize water vapor, the most difficult component to model. At the same time, pressure-height reduction eliminates the systematic bias introduced by the inconsistency between the reanalysis topographic height and the actual height of the station, ensuring that the wet delay caliber of all stations in the network remains consistent.

[0052] 6. Strong interpretability. The five coefficients of the baseline model have clear physical meanings; the residual layer provides characteristic attribution through the Shapley Additive exPlanations (SHAP) method. The results show that the first-level residual output, the first lag of the wet delay residual, the first difference of water vapor pressure, and the annual accumulation periodic term are the main contributing factors, which are consistent with the physical understanding of the time inertia and seasonal rhythm of water vapor. The model behavior is verifiable.

[0053] 7. Decoupled structure, easy to deploy in engineering. The three-level architecture is connected by a fixed three-term summation contract, and each level can be trained, saved and replaced independently; the real-time prediction baseline is calculated analytically, and each of the two levels of residuals requires one forward inference, which reduces the amount of computation and facilitates business operation.

[0054] 8. Provide extrapolation reliability criteria. This invention simultaneously outputs the relationship between the convex hull and the interior / exterior of the predicted point and the distance to the nearest station, enabling downstream precision positioning or water vapor inversion applications to reasonably weight the predicted values, rather than using them equally without distinction. Attached Figure Description

[0055] Figure 1 This is a schematic diagram of the overall process of the modeling method described in this invention; Figure 2 This is a scatter plot of the density of observed and predicted ZWD values ​​for 2024 using the final model in this embodiment of the invention. Figure 3 This is a comparison chart of the cumulative contributions of the three-level residual architecture in this embodiment of the invention: the SCHA baseline, the first-level CNN-BiLSTM-Attention residual network, and the second-level XGBoost residual model. Detailed Implementation

[0056] The technical solution of the present invention will be further described below with reference to the accompanying drawings and specific embodiments. The embodiments described are only for explaining the present invention and are not intended to limit the scope of protection of the present invention.

[0057] Example The tropospheric products published by the GNSS reference station network provide ZTD (Zero Delay Time), not ZWD (Zero Delay Time) required by this invention. The two differ by about two meters in ZHD (Zero Delay Time). If ZTD is directly used as the modeling target, the static component will account for about 90% of the total, and the error of the wet component will be severely diluted in the statistical indicators. Therefore, this invention first constructs a physically consistent wet delay observation sequence as the modeling target and evaluation benchmark for subsequent levels.

[0058] This embodiment uses tropospheric products from the European Permanent GNSS Network (EPN) and meteorological parameters from the fifth-generation atmospheric reanalysis (ERA5) data as its data foundation, spanning from January 1, 2020 to December 31, 2024, with a sampling interval of 6 hours. Data from 2020 to 2023 were used for model fitting, while the entire year of 2024 data was used only for independent testing and did not participate in any training phase. The training set contains 360 stations and 1,696,695 records; the test set contains 385 stations and 514,526 records. After removing initial samples that could not fully construct historical time-series characteristics, the common samples participating in the unified evaluation amounted to 415,518 records, accounting for 80.76% of the test set. The station's latitude ranges from 30.6°N to 70.3°N, its longitude ranges from 15.2°W to 35.2°E, and its geodetic height ranges from 26.3 meters to 2508.4 meters. Figure 1 This is a schematic diagram of the overall process of the modeling method described in this invention.

[0059] The zenith wet delay prediction method based on periodic spatiotemporal modeling and cascaded residual correction includes the following steps: Step S1: Data Acquisition and Preprocessing The tropospheric products in this embodiment are taken from the publicly available product service of EPN and published in the SINEX_TRO format of the International GNSS Service (IGS). Weekly solutions are preferred, and daily solutions are used to supplement missing ones. After the product files are decompressed as needed, the latitude, longitude, and geodetic height of each station are read from its station information block, and the ZTD of each epoch is read from its tropospheric solution block. For earlier products that only provide Earth-Centered Earth-Fixed (ECEF) rectangular coordinates, geodetic coordinates are calculated using an iterative method. The hourly solution timestamps of this product are mostly marked at half-hours. Therefore, for the four target times of 0:00, 6:00, 12:00, and 18:00, the nearest sample is taken within a 40-minute tolerance and its timestamp is aligned to the target whole hour.

[0060] Meteorological parameters were obtained from the surface analysis field and invariant field of ERA5 atmospheric reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF). To ensure download feasibility, the data was cropped on the server according to the modeling region, with the cropping range from 28.5°N to 72.5°N and 17.5°W to 37.5°E. The parameters included surface air pressure, mean sea level pressure, 2-meter air temperature, 2-meter dew point temperature, and total atmospheric water vapor. The invariant field consisted of the land-sea mask and surface geopotential. Time sampling was consistent with the tropospheric product, using 0:00, 6:00, 12:00, and 18:00 daily. The reanalysis topographic height was obtained by dividing the surface geopotential by the standard gravitational acceleration of 9.80665 m / s². All fields were bilinearly interpolated to the plane position of each reference station.

[0061] The water vapor pressure is calculated from the 2-meter dew point temperature according to formula (1), and the surface air pressure is reduced from the reanalysis topographic height to the station geodetic height according to formula (2). Then, ZHD is calculated according to formula (3) and subtracted from ZTD to obtain the ZWD observation value.

[0062] In terms of quality control: records with ZWD not greater than zero were removed; wet delay was converted to PWV based on weighted average temperature, and records with PWV not in the 0 to 100 mm range were removed; the same physical stations were merged according to the first four digits of the station code, and duplicate records were removed; records with missing key fields such as air pressure, temperature, water vapor pressure, total water vapor, elevation difference, time period code, and various delay values ​​were removed; for stations with coordinate errors found during the merging process, the coordinates were corrected according to their original correct codes. In this embodiment, a total of 2 stations and 1057 records were corrected; stations with fewer than 30 samples in 2024 were removed. In this embodiment, a total of 2 stations and 55 records were removed.

[0063] After the above processing, the dataset from 2020 to 2023 contains 400 stations and 1,893,356 records, and the dataset from 2024 contains 385 stations and 514,526 records. The wet delay ranges from 0.0002 meters to 0.4068 meters, which is consistent with the physical magnitude of wet delay.

[0064] Step S2: Station-by-station harmonic fitting For the wet delay sequence of each station from 2020 to 2023, a five-parameter harmonic model is fitted according to Equations (4), (5) and (6) to obtain five coefficients: constant term, annual cycle cosine, annual cycle sine, semi-annual cycle cosine, and semi-annual cycle sine. Through this step, the time variation of each station is compressed into five scalars with physical meaning.

[0065] Step S3: SCHA Spatial Modeling The crown center was determined to be 49.023°N, 13.578°E based on the mean latitude and longitude of the stations involved in the fitting. The spherical angular distance from the crown center to the farthest station determined the crown half-angle to be 29.413°. The upper limit for the azimuth order was set to 4, and the upper limit for the degree was set to 8. A family of Neumann solutions with an even difference between the degree and order was used for selection, resulting in 19 sets of non-integer degrees, constituting 33 crown harmonic basis functions. The reference height was set to 172.969 meters, and the exponential height scaling scale was set to 3000 meters through a grid search. Regularization coefficients were selected using a 10-point logarithmic grid from 0 to 0.01, with five-fold spatial block cross-validation. All five coefficients were set to 0.01, corresponding to a design matrix condition number of 782.27, indicating good numerical conditions. The specific method for spatial partitioning is as follows: First, the stations are sorted lexicographically by latitude and longitude. Then, the sorting results are divided into five consecutive segments, ensuring that stations within the same segment are geographically close, thereby cross-validating the spatial extrapolation capability of the true metric. 330 stations participate in spatial fitting, while the remaining stations serve as a spatial holdout set, whose observations are not involved in baseline fitting.

[0066] In terms of fitting accuracy, the spatial reconstruction RMSE of the five harmonic coefficients were 2.92 mm, 3.10 mm, 1.48 mm, 1.17 mm and 1.56 mm, respectively, and the bias was less than 0.13 mm, indicating that the spherical cap harmonic basis adequately characterizes the spatial field of the coefficients.

[0067] Step S4: Construction of multi-scale temporal features The system constructs basic input information, vapor pressure time series information, and residual history information. The basic input information includes latitude, longitude, geodetic height, temperature, station air pressure, vapor pressure, and sine and cosine codes for the annual accumulated days, totaling 8 items. The vapor pressure time series information includes the 1st, 2nd, and 4th order differences of vapor pressure, as well as the rolling mean and rolling standard deviation for 4, 8, and 16 time windows, totaling 9 items. The residual history information includes the 1st to 3rd order lags of the residuals, the 1st to 2nd order differences, and the rolling mean and rolling standard deviation for 4 time windows, totaling 7 items. To avoid using information after the current time, all information involved in time series calculations is processed chronologically within the same station, and the corresponding series is lagged by one time point before calculation. When the time interval between two adjacent records from the same station exceeds 6 hours, they are divided into different continuous series, and lag, difference, and rolling statistics are not calculated across series. At the same time, a list of prohibited input variables is set to exclude variables such as measured ZWD, ZTD, and ZHD, which may directly or indirectly contain information about the prediction target, from the model input.

[0068] Step S5: Training the first-level residual network (CNN-BiLSTM-Attention) The first-level residual network was constructed, consisting of a one-dimensional convolutional neural network layer, a bidirectional long short-term memory network, an additive attention convergence layer, and a fully connected output layer. The input sequence length was 16 six-hour time points, covering the history of the four days prior to the target time point; the CNN output had 48 channels and a kernel size of 3; the BiLSTM unidirectional hidden state dimension was 96, with 1 layer; the fully connected hidden layer after attention convergence had a dimension of 64; and the Dropout probability was 0.15. The optimizer used was AdamW, with an initial learning rate of 0.0003, weight decay of 0.0001, a batch size of 1024, and a gradient norm pruning threshold of 5. The training set is from January 1, 2020 to December 31, 2022; the fitting interval is from January 1, 2023 to September 30, 2023, and the early stopping interval is from October 1, 2023 to December 31, 2023. The two intervals do not overlap; the maximum number of training rounds is 50. If the loss in the early stopping interval does not improve for 8 consecutive rounds, it is considered an early stopping; in this embodiment, the optimal validation loss is obtained in the 20th round.

[0069] Step S6: Training the second-level residual model (XGBoost) The second-level model input consists of all causal features from step S4 plus the first-level residual predictions from step S5. The maximum tree depth is 4, the learning rate is 0.02, the minimum sample weight sum of child nodes is 10, the row sampling ratio is 0.90, the column sampling ratio is 0.85, the first-order regularization coefficient γ is 0.5, the second-order regularization coefficient λ_reg is 10.0, and the number of early stopping rounds is 250. The fitting interval is from January 1, 2023 to September 30, 2023, and the early stopping interval is from October 1, 2023 to December 31, 2023; the two intervals do not overlap. In this embodiment, the final number of trees selected is 122.

[0070] Step S7: Prediction Result Synthesis and Extrapolation Judgment After training, the geographic latitude, longitude, geodetic height and target time of the point to be predicted are input into steps S3 to S6. 24-dimensional historical features of 16 consecutive 6-hour time periods are constructed according to step S4 to obtain the baseline prediction value, the first-level residual prediction value and the second-level residual prediction value respectively. The 24-dimensional current features and the first-level residual prediction value are concatenated into a 25-dimensional XGBoost input. The three are added according to formula (20) to obtain the final prediction result. The extrapolation state is determined according to formula (21) and the convex hull of 330 training stations. The nearest station distance threshold in this embodiment is 188.7 km.

[0071] Verification results For the ZWD observations of the independent test samples in 2024, the predicted values ​​are calculated using the method of this invention. The accuracy evaluation index is calculated according to formula (22), where N is the number of evaluation samples, y_i is the i-th measured value, y_hat_i is the corresponding predicted value, and y_bar is the measured mean.

[0072] The step-by-step accuracy contribution of the three-tier architecture is shown in Table 1. All results are from independent tests in 2024, and the unit is millimeters.

[0073] Table 1 shows that: First, the two-stage residual learning reduced the RMSE from 36.00 mm to 16.49 mm relative to the SCHA baseline, a reduction of approximately 54.2%. The first-stage network contributed the majority of the reduction, while the second-stage model further stabilized and reduced the RMSE by approximately 0.28 mm. Second, the accuracy of spatially reserved stations and completely new stations was basically equivalent to that of existing stations, and even better in some ranges, indicating that the model's accuracy did not come from memorizing specific stations. Third, the bias in each range was within 1 mm, and there was no systematic overestimation or underestimation. As a comparison, the persistence prediction RMSE directly extrapolated from the previous time-of-flight observation was 21.04 mm, and this invention reduced it by approximately 21.6% relative to that benchmark.

[0074] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit them. After reading this specification, those skilled in the art can make various modifications or variations to the technical solutions without departing from the concept of the present invention, and all such modifications or variations fall within the protection scope of the present invention.

Claims

1. A zenith wet delay prediction method based on periodic spatiotemporal modeling and cascaded residual correction, characterized in that, Includes the following steps: Step S1: Obtain the GNSS tropospheric zenith total delay (ZTD) product aligned by station and time, and reanalysis meteorological data covering the same latitude and longitude range; interpolate the meteorological field to the station location, calculate the pressure-reduced zenith static delay (ZHD), and then calculate the ZWD; output the ZWD observation sequence and corresponding exogenous variables at the unified sampling time, and provide them to steps S2 and S4. Step S2: Using the ZWD observation sequences of each station output in Step S1 as input, first determine the annual cycle and semi-annual cycle as the main cycles through spectrum analysis, and then establish a five-parameter harmonic regression model with a fixed cycle; solve and output the five harmonic coefficients of each station. Step S3: Using the five harmonic coefficients output from Step S2, the station coordinates, and the geodetic height as inputs, the center of the spherical cap is determined by the mean latitude and longitude of the stations involved in the fitting, and the half-angle of the spherical cap is determined by the spherical angular distance from the center of the spherical cap to the farthest station. The SCHA basis function that satisfies the Neumann boundary conditions is constructed, and the exponential height scaling factor is superimposed. The ridge regression regularization coefficient is selected through spatial block cross-validation, and the baseline ZWD prediction value at any location and at any time within the modeling area is output. Step S4: Subtract the ZWD observation value output in Step S1 from the baseline ZWD prediction value output in Step S3 to obtain the first-level residual objective; sort the exogenous variables and residual history sequences by station and time and apply causal lag, construct an input feature set that does not contain measured ZWD, ZTD, ZHD and their direct conversion quantities, where each time point consists of a 24-dimensional vector composed of 8 basic features, 9 water vapor pressure derived features and 7 residual history features, and output a fixed-length history sequence for use in Step S5; Step S5: Input the 16 consecutive 6-hour time-series of each station output from Step S4, with each time-series consisting of 24 dimensions, into the first-level residual network composed of a one-dimensional convolutional neural network (CNN), a bidirectional long short-term memory network (BiLSTM), an additive attention convergence layer, and a fully connected output layer. After extracting local changes and bidirectional temporal information, the first-level residual prediction value is output. Step S6: Concatenate the input feature set output in step S4 with the first-level residual prediction value output in step S5 as the input to the gradient boosting decision tree implemented by XGBoost. Subtract the first-level residual prediction value from the first-level residual objective to obtain the second-level learning objective. Train and output the second-level residual prediction value. Step S7: Add the baseline ZWD prediction value output in step S3, the first-level residual prediction value output in step S5, and the second-level residual prediction value output in step S6 to obtain the final ZWD prediction value. Based on the relationship between the point to be predicted and the convex hull of the set of stations to be fitted, and whether the spherical distance to the nearest fitted station exceeds the distance threshold, output the extrapolation reliability indicator of the prediction value.

2. The method according to claim 1, characterized in that, Step S1 is as follows: Step S11: Obtain the tropospheric delay product in standard exchange format, parse its station information block and tropospheric solution block to obtain the latitude, longitude, geodetic height and ZTD of each station for each epoch. Step S12: For tropospheric products whose sampling time does not fall on a whole hour, take the sample closest to the target whole hour within the preset time tolerance range, and align its timestamp to the target whole hour so that data from different sources have a unified sampling time. Step S13: Clip the ERA5 reanalysis data according to the latitude and longitude range of the modeling area to obtain the surface analysis field including surface air pressure, two-meter air temperature, two-meter dew point temperature, and total water vapor in the whole layer, as well as the invariant field including the land-sea mask and surface potential. Divide the surface potential by the standard gravitational acceleration to obtain the reanalysis terrain height. Interpolate the above fields to the plane position of each reference station using bilinear interpolation. Step S14: Calculate the vapor pressure based on the two-meter dew point temperature; Step S15: Reduce the reanalysis surface pressure from the reanalysis topographic height to the station geodetic height; Step S16: Calculate ZHD from the reduced station air pressure and derive the ZWD observation sequence for each station; Step S17: Apply quality control to the ZWD observation sequence, including: removing non-physical records with ZWD not greater than zero, removing records where the precipitable water obtained from ZWD conversion exceeds the preset physical range, removing records with missing key fields, correcting the coordinates of stations found to have coordinate errors according to their original correct codes, and removing stations with sample sizes below the preset threshold during the evaluation period, thus obtaining the dataset used for modeling.

3. The method according to claim 1, characterized in that, The fixed-period five-parameter harmonic regression model is as follows: In the formula, t represents the accumulated days of a year. The fundamental angular frequency, , , , , These are the five coefficients of the five-dimensional harmonic coefficients; The calculation formula is as follows: The annual and semi-annual cycles are determined as the main periodic components by spectral analysis of the ZWD sequence; the five coefficients are solved according to the least squares criterion, and the missing epochs are not filled in for modeling objectives. Only the effective observation times are used to construct the design matrix and complete the solution.

4. The method according to claim 1, characterized in that, Step S3 includes: Step S31: Using the mean latitude and longitude of all stations participating in the fitting as the center of the spherical cap, and the spherical angular distance from the center of the cap to the farthest station as the half-angle of the spherical cap, rotate the geographical coordinates of each station into the co-latitude and azimuth of the spherical cap with the center of the cap as the pole. Rotate the geographic coordinates of the station or the point to be predicted into spherical coordinates: in, and Input the geographical latitude and longitude of the point. and The geographical latitude and longitude of the crown of the globe. The coaxial latitude of the spherical cap relative to the coronal center is the input point; Step S32: For each azimuth order m, using the Neumann boundary condition (i.e., the first derivative of the associated Legendre function along the colatitude direction is zero at the spherical cap boundary), solve for a set of non-integer degrees, and retain only the solution family where the difference between the degree and the order is even; as follows: in Let be the associated Legendre function of degree n and orientation order m. Given the spherical cap's half-angle, the boundary condition states that the first derivative of the function along the colatitude direction is... The value is zero; Step S33: Using the Schmidt quasi-normalized associated Legendre function as the radial basis, expand the j-th harmonic coefficient into a spatially continuous field according to the following formula: in, and Let be the harmonic expansion coefficient of the spherical cap to be determined; The degree is a non-integer number determined by the boundary conditions; The reference height is in meters; H is the exponential height scaling scale, in meters; m is the azimuth order, k is the number index, and M and K are the azimuth order and the upper limit of the number, respectively. Let be the Schmidt quasi-normalized associated Legendre function, and denote the Schmidt quasi-normalization; Step S34: Solve for the harmonic coefficients of the spherical cap using Ridge Regression: In the formula, X is the design matrix composed of the basis functions described in step S33; Let be the spherical cap harmonic expansion coefficient vector of the spatial field of the j-th time harmonic coefficient; Let J be the harmonic coefficient observation vector of the j-th station; Let be the ridge regression regularization parameter corresponding to the j-th coefficient field.

5. The method according to claim 4, characterized in that, After spatial modeling is completed, the five coefficient spatial fields output in step S3 can take values ​​at any point to be predicted, and the baseline wet delay at any location and time can be reconstructed: in, The j-th time harmonic coefficient is the spatial field value at the spherical cap coordinates and the ground height. The coefficient is derived from step S3 for each station. Perform SCHA expansion and ridge regression; Let t be the geodetic height of the point to be predicted, and t be the annual day or equivalent time code of the target time. and Given the geographical latitude and longitude of the point to be predicted. and Given the coaxial latitude and azimuth of its spherical crown; The time basis function multiplied by the j-th coefficient is specifically... ,in .

6. The method according to claim 4, characterized in that, Regularization coefficient The selection process employs spatial block cross-validation: first, the stations participating in the fitting are sorted lexicographically by latitude and longitude; then, the sorted stations are divided into several non-overlapping spatial blocks. One spatial block is reserved as a validation block at a time, and the remaining blocks participate in the fitting. The block whose sum of the root mean square errors of the coefficients in each validation block is selected. .

7. The method according to claim 1, characterized in that, Step S4: Input the feature set. The eight basic features at each time point are latitude, longitude, geodetic height, air temperature, station air pressure, water vapor pressure, and sine and cosine codes for annual accumulated days; the nine derived features of water vapor pressure are the first, second, and fourth order differences of water vapor pressure, as well as the rolling mean and rolling standard deviation of time windows 4, 8, and 16; the seven residual history features are the first to third order lags and first to second order differences of the first-level residuals, as well as the rolling mean and rolling standard deviation of the four time windows, totaling 24 dimensions.

8. The method according to claim 1, characterized in that, Step S5 takes a 24-dimensional historical feature sequence of 16 consecutive 6-hour time points from the same station as input, covering the 4 days prior to the target time; information after the current time, the measured ZWD at the target time, and future residuals are not included in the input; the sequence is first processed by a one-dimensional convolutional neural network (CNN) layer to extract local change patterns between adjacent time points, and then fed into a bidirectional long short-term memory (BiLSTM) network. in, Let be the 24-dimensional input vector at time t. The number of time steps covered by the convolution kernel. and Here, BN represents batch normalization, and ReLU represents the modified linear unit activation function. The CNN outputs features, which are then used as input to the subsequent BiLSTM, and the CNN output sequences at each time step are processed. Inputting into a BiLSTM yields forward and backward hidden states; the forward hidden state is then converted into a reverse hidden state. With reverse hidden state spliced ​​as : Subsequently, additive attention is used to weight and converge the hidden states at each time step, so that the model automatically emphasizes the historical moments that contribute more to the current residual; in, Rate attention For normalized weights, z is the context vector, and W, b, and v are learnable parameters; the converged z is passed through a fully connected layer to output the first-level residual prediction value: in and These are the output layer parameters.

9. The method according to claim 1, characterized in that, Step S6 receives the 24-dimensional current features from step S4 and the first-level residual prediction values ​​from step S5, and constructs the second-level learning objective: in, The residual objective learned by the second-level model is equal to the measured ZWD minus the baseline ZWD and the predicted value of the first-level residual; therefore, the input dimension of the second-level model is 24+1=25 dimensions, and the output is the predicted value of the second-level residual. GBDT tree-by-tree additive fitting implemented using XGBoost: Where η is the learning rate. For the s-th regression tree, Let S be the set of candidate regression tree functions, and S be the final number of iteration trees. This represents the second-level XGBoost residual prediction value obtained by summing the results of S GBDT trees one by one with a learning rate η. Its regularization objective function is: in, For single-sample loss function, It is the first True values ​​of the second-order residuals for each sample These are the corresponding second-level XGBoost residual predictions, and the two constitute the loss function l( , ); Let be the number of leaf nodes in the s-th tree. The leaf node weight vector, The penalty coefficient for the number of leaf nodes. is the L2 regularization coefficient for the leaf node weights.

10. The method according to claim 1, characterized in that, Step S7: Given the geographical latitude of the point to be predicted. ,longitude , the earth is high Given the target time t, first calculate the baseline ZWD according to step S3; then construct the 24-dimensional historical features of the same station for 16 consecutive 6-hour periods up to t according to step S4, and load the CNN-BiLSTM-Attention model to obtain... Combine the current 24-dimensional features with The vectors are concatenated into a 25-dimensional vector and input into the XGBoost model to obtain the second-level residual prediction value of the XGBoost model for the target point after training. Adding the three terms together gives the final ZWD prediction value: Then, a convex hull is constructed using the planar positions of the stations involved in the fitting, and the spherical distance from the point to be predicted to the nearest fitted station is calculated: In the formula, R is the average radius of the Earth, and Δφ and Δλ are the difference in latitude and longitude, respectively. and Latitude of two points; spherical distance It is used to measure the spatial proximity between the point to be predicted and the nearest fitted station and compare it with a preset distance threshold. When the point to be predicted is outside the convex hull or the distance to the nearest station exceeds the threshold, it is marked as an extrapolation point and output along with the predicted value, so that the weight of the predicted value can be adjusted according to reliability for downstream positioning or water vapor inversion.

Images