A GNSS Vertical Deformation Signal Extraction Method Based on Multi-Source Load Correction and CME Filtering

Abstract

This invention discloses a method for extracting GNSS vertical deformation signals based on multi-source load correction and CME filtering. It constructs a load correction model for common-mode error filtering and applies non-tidal atmospheric, oceanic, and hydrological load effects. Through a five-step collaborative process—constructing a standardized GNSS time-series dataset, performing multi-source load hierarchical correction, optimizing CME filtering, dynamically adapting a noise model, and extracting and estimating vertical velocity—the method quantifies the spatial heterogeneity of GNSS vertical time-series nonlinear signals, background noise, and vertical velocity. By combining the correction derived from GRACE with CME filtering, this invention enables the quantification of the spatial heterogeneity of GNSS vertical time-series nonlinear signals, background noise, and vertical velocity, accurately depicting tectonic signals within complex suture zones of glaciers and plateaus.

Description

Technical Field

[0001] This invention relates to the field of geodesy and geophysical monitoring technology, specifically to a method for extracting GNSS vertical deformation signals based on multi-source load correction and CME filtering. It is particularly applicable to geodynamic studies (such as mid-crustal flow and plate subduction response) in glacial plateau regions, including the Qinghai-Tibet Plateau, monitoring of cryosphere-hydrosphere-lithospheric interactions (such as elastic deformation caused by glacial melting and groundwater changes), and regional crustal stability assessment, providing key data support for the analysis of deep tectonic processes in the plateau and the assessment of geological hazard risks. Background Technology

[0002] One of the core objectives of Global Navigation Satellite Systems (GNSS) is to monitor long-term, stable tectonic movements of the Earth's crust (such as plateau uplift and plate compression). By analyzing long-term time series data from continuously operating GNSS observation stations, millimeter-level precision three-dimensional crustal displacement information can be obtained, providing crucial support for earthquake prediction, plate boundary activity analysis, glacier mass change, and groundwater storage change studies. In these applications, vertical deformation is particularly sensitive because it directly reflects changes in surface mass and the elastic response of the crust. However, compared to the horizontal component, GNSS vertical displacement is more susceptible to various non-tectonic factors, resulting in more complex signal noise and poorer data stability. Therefore, more refined signal extraction and error correction methods are needed. GNSS observation data from the Tibetan Plateau region contains rich tectonic deformation information, but due to strong interference from non-tectonic loads such as atmospheric, oceanic, and hydrological factors, as well as background noise (such as flicker noise and power-law noise), traditional GNSS data processing methods struggle to accurately extract true tectonic movement signals. Existing technologies suffer from the following problems:

[0003] First, nonlinear signals mask structural features.

[0004] GNSS receivers record the sum of all surface displacements. Among them, the non-tectonic nonlinear signals caused by geophysical loads have a much larger amplitude than long-term tectonic signals. Strong seasonal signals (annual or semi-annual cycles) can cause deviations in the long-term trend fitting of time series, resulting in inaccurate estimates of vertical velocity. Unmodeled nonlinear signals can significantly increase the uncertainty of GNSS velocity field estimation, causing the real tectonic deformation to be masked by strong noise background.

[0005] Second, the background noise has high structural complexity and high spatiotemporal correlation.

[0006] GNSS time series contain various noise components, including white noise, flicker noise, and power-law noise. For example, on the Tibetan Plateau, noise characteristics are closely related to local water-ice coupling, permafrost dynamics, and tectonic activity, exhibiting strong spatial variability. Existing noise models (such as WN+FN) are overly simplified and fail to effectively separate common-mode errors related to spatial location, leading to high uncertainty in velocity field estimation.

[0007] Third, the shortcomings of traditional load correction models and the neglect of spatial heterogeneity.

[0008] Taking the widely used GFZ hydrological load model (HYDL) as an example, it is mainly based on soil moisture and surface flow models, but it severely neglects the changes in the mass of glaciers and ice sheets, as well as changes in groundwater storage. In the Himalayas, where glaciers are widespread, and in areas with drastic groundwater changes, the incompleteness of this model results in a large amount of hydrological load deformation remaining uncorrected, persisting in the time series as "residual signals" and continuously interfering with the extraction of tectonic signals. Existing techniques often adopt a uniform correction strategy for the entire study area, ignoring the differences in the dominant noise sources in different regions.

[0009] Fourth, there is a lack of cross-scale signal separation technology.

[0010] Existing methods struggle to simultaneously handle seasonal (annual, semi-annual) and interannual signals, and fail to combine with wavelet spectral analysis to verify the interannual signal suppression effect. They also cannot quantify the noise reduction contribution of different correction schemes to the interannual cycle, resulting in insufficient long-term stability of GNSS time series and making it difficult to support detailed studies of geodynamics on the Tibetan Plateau.

[0011] In summary, existing technologies cannot fully address the interference problem between nonlinear signals and background noise in GNSS vertical time series data from glaciers and plateaus. There is an urgent need for a systematic processing method that can integrate multi-source load data, adapt to spatial heterogeneity of noise, and effectively reduce the uncertainty of velocity estimation, so as to accurately extract tectonic deformation signals and provide a reliable data foundation for geodynamic research on glaciers and plateaus. Summary of the Invention

[0012] The purpose of this invention is to provide a method for extracting GNSS vertical deformation signals based on multi-source load correction and CME filtering. A load correction model integrating GFZ and GRACEmascon data (CSR / JPL / GSFC) is constructed, followed by common-mode error (CME) filtering. Non-tidal atmospheric (NTAL), oceanic (NTOL), and hydrological (HYDL) load effects are applied. By combining the correction derived from GRACE with CME filtering, this invention can quantify the spatial heterogeneity of GNSS vertical time series nonlinear signals, background noise, and vertical velocity, accurately depicting tectonic signals within the complex suture zone of the glacier plateau, thus providing crucial data support for the geodynamic processes of the entire plateau.

[0013] To achieve the above-mentioned technical objectives, the technical solution adopted by the present invention is as follows:

[0014] A method for extracting GNSS vertical deformation signals based on multi-source load correction and CME filtering, the method comprising the following steps:

[0015] S1. Collect continuous GNSS observation data over a preset time span to ensure that the collected data contains complete seasonal and interannual signal cycles. After correcting outliers, construct a standardized GNSS vertical time series dataset without outliers.

[0016] S2. Calculate the non-tidal atmospheric load vertical displacement time series and the non-tidal ocean load vertical displacement time series for each GNSS station using the load Green's function. Superimpose the non-tidal atmospheric load vertical displacement time series and the non-tidal ocean load vertical displacement time series to obtain the total AO load displacement time series. Subtract the total AO load displacement time series from the standardized GNSS vertical time series to obtain the AO-corrected GNSS vertical time series.

[0017] S3. Based on the mass change time series collected by GRACE satellite, the equivalent water height reflecting the change of total land water storage is calculated, and the hydrological vertical displacement time series is obtained based on the load Green's function. The hydrological vertical displacement time series is superimposed on the AO-corrected GNSS vertical time series to obtain the multi-source load-layered corrected GNSS vertical time series.

[0018] S4. Select GNSS stations observed in the same period, perform EOF decomposition on the GNSS vertical time series with multi-source load hierarchical correction, calculate the CME time series based on the eigenvectors and eigenvalues ​​obtained from the decomposition, and then subtract the CME time series from the AOG-corrected GNSS vertical time series to obtain the AOG_CME-corrected GNSS vertical time series.

[0019] S5. Define multiple noise models according to the type of noise on the plateau. Fit the GNSS vertical time series after AOG_CME correction in step S4 using multiple noise models. Select the optimal noise model for each GNSS station based on the fitting structure to achieve dynamic adaptation.

[0020] S6 estimates the vertical velocity and its uncertainty for each GNSS station based on the optimal noise model.

[0021] Step S1 further includes:

[0022] Continuous GNSS observation data from the China Crustal Movement Observation Network and the Nevada Geodetic Laboratory were collected over a preset time span. GNSS stations with a time series length greater than a preset duration threshold were selected, and stations with a data missing rate greater than a missing rate threshold were removed to form a GNSS station network covering the entire Qinghai-Tibet Plateau.

[0023] The coordinate time series under the ITRF2014 framework was obtained using GipsyX software. Linear interpolation was used to replace gross errors. For step outliers caused by instrument changes or regional tectonic activities, step parameters were estimated and removed using known instrument replacement records or regional tectonic activity timestamps. For abrupt trend changes, piecewise linear fitting was used for repair. Finally, a standardized GNSS vertical time series without outliers was formed.

[0024] Step S2 further includes:

[0025] Three-hour resolution atmospheric pressure data with a time span consistent with the GNSS vertical time series were collected. Harmonic analysis was performed on 12 major atmospheric tidal components to remove tidal atmospheric effects and retain non-tidal components. Based on the elastic Earth load theory, the non-tidal atmospheric load model provided by GFZ was used to calculate the non-tidal atmospheric load vertical displacement time series for each GNSS station through the load Green's function. :

[0026] ;

[0027] in, At time t Non-tidal atmospheric pressure anomaly at the location For the load Green's function, The geocentric angle is used, and the integration range covers the entire globe;

[0028] Collect Max Using 3-hour resolution seafloor pressure data from the Max Planck Institute ocean model, and based on the GFZ non-tidal ocean loading model, the time series of vertical displacement of non-tidal ocean loading was calculated using the load Green's function. Among them, considering the boundary effect between the ocean and the land, a refined interpolation method is used to process the displacement data in the area near the coastline;

[0029] The total AO load displacement time series is obtained by superimposing the non-tidal atmospheric load vertical displacement time series and the non-tidal ocean load vertical displacement time series. The total AO load displacement time series is then subtracted from the standardized GNSS vertical time series to obtain the AO-corrected GNSS vertical time series. :

[0030] ;

[0031] in, To standardize GNSS time series, This is the AO-corrected GNSS vertical time series.

[0032] Step S3 further includes:

[0033] The mass change time series data collected by the GRACE satellite were acquired, and the acquired mass change time series data were preprocessed. The preprocessing process included: using a 30-day moving average to eliminate high-frequency noise in the mass change time series data, and using scale factor correction to eliminate leakage errors contained therein.

[0034] The preprocessed time series data of quality changes were converted into equivalent water height. :

[0035] ;

[0036] in, At time t The gravity anomaly at the location The density of water, It is the acceleration due to gravity;

[0037] Based on the load Green's function, the equivalent water height is converted into a hydrological vertical displacement time series using the following formula. :

[0038] ;

[0039] Hydrological vertical displacement time series Superimposed on the AO-corrected GNSS vertical time series, a multi-source load-layered corrected GNSS vertical time series is obtained:

[0040] .

[0041] Step S4 further includes:

[0042] Using a multi-source load-layered correction GNSS vertical time series, several GNSS stations observed simultaneously were selected for analysis. EOF decomposition yields spatial eigenvectors, temporal eigenvectors, and eigenvalues, as shown in the following formula:

[0043] ;

[0044] in, Let be an n×m GNSS station time series matrix, where n represents the number of GNSS stations and m represents the number of time points. Let n×n be the spatial eigenvector matrix. Let m be an n×m eigenvalue diagonal matrix. The time eigenvector matrix is ​​m×m;

[0045] Based on the variance contribution of the eigenvalues, the top k spatial modes are selected as the spatial modes of CME, and their products with the corresponding time eigenvectors and eigenvalues ​​are used as the CME time series. :

[0046] ;

[0047] in, These are the first k spatial feature vectors. For the corresponding eigenvalues, These are the feature vectors of the first k time periods;

[0048] The extracted CME time series is subtracted from the multi-source load-stratified GNSS vertical time series to obtain the AOG_CME-corrected GNSS vertical time series:

[0049] .

[0050] Step S5 further includes:

[0051] Five noise models are defined: the WN model containing only white noise with parameters equal to the standard deviation of white noise; the WN+FN model containing white noise and flicker noise with FN parameters equal to the spectral exponent and noise amplitude; the WN+FN+RWN model containing white noise, flicker noise, and red noise with RWN parameters equal to the spectral exponent and noise amplitude; the WN+GGM model containing white noise and generalized Gaussian-Markov noise with GGM parameters equal to the correlation time and noise amplitude; and the WN+PL model containing white noise and power-law noise with PL parameters equal to the spectral exponent and noise amplitude.

[0052] AOG_CME corrected GNSS time series for each GNSS station Five noise models were used for fitting, and the BIC_tp value of each noise model was calculated using the following formula:

[0053] ;

[0054] in, The amount of observation data from GNSS stations. The residual variance of the noise model. The number of parameters in the noise model;

[0055] The noise model with the smallest BIC_tp value is selected as the optimal noise model for the GNSS site.

[0056] Furthermore, in step S6, the velocity parameters are optimized using maximum likelihood estimation, as shown in the following formula:

[0057] ;

[0058] in, For a design matrix that includes a time trend item, This is the noise covariance matrix constructed based on the optimal noise model. GNSS vertical time series .

[0059] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0060] The GNSS vertical deformation signal extraction method based on multi-source load correction and CME filtering of this invention has significant innovation and practical value in terms of theoretical model, data processing flow and result reliability, mainly reflected in the following aspects:

[0061] First, multi-source load synergistic correction: For the first time, the GFZ model and GRACE data are integrated, achieving an RMS reduction rate of 35% in glacier-covered areas, which is significantly better than traditional methods. Compared with single-source correction methods, the integrated load correction of this invention can effectively eliminate the superposition effect of different load sources, making the GNSS vertical time series closer to the actual crustal movement.

[0062] Second, CME filtering region adaptation: The CME filtering algorithm based on EOF decomposition design is adapted to the regional characteristics of spatial correlation errors in glaciers and plateaus. After filtering, the velocity uncertainty is reduced by 26.9%, and the noise model complexity is significantly simplified.

[0063] Third, dynamic noise model selection: Based on the BIC_tp criterion, dynamic noise model adaptation is implemented to achieve "one model per station", improving the accuracy of velocity estimation and successfully identifying a weak tectonic uplift signal of +1.2 mm / yr; the optimal noise model for each station is determined, making velocity estimation and uncertainty assessment more consistent with the actual noise structure of each station. Compared with the traditional approach of using a unified noise model, this strategy significantly improves the model's adaptability and regional accuracy.

[0064] Fourth, standardized and portable process: The standardized process from data preprocessing to signal extraction can be transferred to other complex geological and climatic regions (such as the Andes Mountains and the Alps), and has broad application prospects; the modular structure of the method allows it to be integrated into existing GNSS data processing software systems (such as GAMIT / GLOBK, Bernese, PyGAMIT), and has good versatility and scalability. Attached Figure Description

[0065] Figure 1 This is a flowchart of the GNSS vertical deformation signal extraction method based on multi-source load correction and CME filtering according to the present invention.

[0066] Figure 2 The diagram illustrates the evaluation results of the optimal background noise model for GNSS time series, where (a) represents the original optimal background noise model for GNSS; (b) represents the optimal background noise model after removing GNSS AO; (c) represents the optimal background noise model after removing GNSS SAOH; (d) represents the optimal background noise model after removing AOG; and (e) represents the optimal background noise model for GNSS after removing AOG and CME. `raw`, `AO-removed`, `AOH-removed`, `AOG-removed`, and `AOG_CME-removed` represent the original image, the image after removing AO, the image after removing AOH, the image after removing AOG, and the image after removing both AOG and CME, respectively.

[0067] Figure 3 The diagram shows the GNSS velocity rate results; where (a)-(e) correspond to the original GNSS velocity rate, AO effect correction, AOH correction, AOG correction, and AOG and CME correction data, respectively.

[0068] Figure 4 The figure shows the uncertainty results; where (f)-(j) correspond to the original GNSS value, AO effect correction, AOH correction, AOG correction, AOG and CME correction data, respectively.

[0069] Figure 5The values ​​v and u represent the original vertical velocity field of the GNSS station under the WN+FN stochastic model after different surface load corrections (corresponding to (a)), velocity field mapping (corresponding to (b)), and uncertainty mapping (corresponding to (c)); the values ​​v and u represent the average velocity and uncertainty of all GNSS stations.

[0070] Figure 6 The vertical velocity and uncertainty of the AO-corrected vertical velocity field of the GNSS station under the WN+FN stochastic model are the vertical velocity and uncertainty after different surface load corrections (corresponding to (d)), velocity field mappings (corresponding to (e)), and uncertainty mappings (corresponding to (f)); the values ​​v and u represent the average velocity and uncertainty of all GNSS stations.

[0071] Figure 7 The vertical velocity and uncertainty of the AOH corrected vertical velocity field of the GNSS station under the WN+FN stochastic model are the vertical velocity and uncertainty after different surface load corrections (corresponding to (g)), velocity field mappings (corresponding to (h)), and uncertainty mappings (corresponding to (i)); the values ​​v and u represent the average velocity and uncertainty of all GNSS stations.

[0072] Figure 8 The vertical velocity and uncertainty of the AOG-corrected vertical velocity field of the GNSS station under the WN+FN stochastic model are the vertical velocity and uncertainty after different surface load corrections (corresponding to (j)), velocity field mappings (corresponding to (k)), and uncertainty mappings (corresponding to (l)); the values ​​v and u represent the average velocity and uncertainty of all GNSS stations.

[0073] Figure 9 The values ​​v and u represent the vertical velocity and uncertainty of the AOG+CME corrected vertical velocity field of the GNSS station under the WN+FN stochastic model after different surface load corrections (corresponding to (m)), velocity field mappings (corresponding to (n)), and uncertainty mappings (corresponding to (o)); the values ​​v and u represent the average velocity and uncertainty of all GNSS stations. Detailed Implementation

[0074] The embodiments of the present invention will be described in further detail below with reference to the accompanying drawings.

[0075] This invention discloses a method for extracting GNSS vertical deformation signals based on multi-source load correction and CME filtering. The method includes the following steps:

[0076] S1. Collect continuous GNSS observation data over a preset time span to ensure that the collected data contains complete seasonal and interannual signal cycles. After correcting outliers, construct a standardized GNSS vertical time series dataset without outliers.

[0077] S2. Calculate the non-tidal atmospheric load vertical displacement time series and the non-tidal ocean load vertical displacement time series for each GNSS station using the load Green's function. Superimpose the non-tidal atmospheric load vertical displacement time series and the non-tidal ocean load vertical displacement time series to obtain the total AO load displacement time series. Subtract the total AO load displacement time series from the standardized GNSS vertical time series to obtain the AO-corrected GNSS vertical time series.

[0078] S3. Based on the mass change time series collected by GRACE satellite, the equivalent water height reflecting the change of total land water storage is calculated, and the hydrological vertical displacement time series is obtained based on the load Green's function. The hydrological vertical displacement time series is superimposed on the AO-corrected GNSS vertical time series to obtain the multi-source load-layered corrected GNSS vertical time series.

[0079] S4. Select GNSS stations observed in the same period, perform EOF decomposition on the GNSS vertical time series with multi-source load hierarchical correction, calculate the CME time series based on the eigenvectors and eigenvalues ​​obtained from the decomposition, and then subtract the CME time series from the AOG-corrected GNSS vertical time series to obtain the AOG_CME-corrected GNSS vertical time series.

[0080] S5. Define multiple noise models according to the type of noise on the plateau. Fit the time series after AOG_CME correction in step S4 using multiple noise models. Select the optimal noise model for each GNSS station based on the fitting structure to achieve dynamic adaptation.

[0081] S6 estimates the vertical velocity and its uncertainty for each GNSS station based on the optimal noise model.

[0082] This invention employs a five-step collaborative process: "constructing a standardized GNSS time-series dataset - multi-source load hierarchical correction - CME filtering optimization - dynamic adaptation of the noise model - vertical velocity extraction and uncertainty estimation" (see [link to documentation]). Figure 1 This enables the quantization of the spatial heterogeneity of GNSS vertical time series nonlinear signals, background noise, and vertical velocity. The specific technical solution is as follows:

[0083] (I) Constructing a standardized GNSS time series dataset

[0084] (1) Establishing a GNSS station network

[0085] Continuous GNSS observation data were collected from the China Crustal Movement Observation Network (CMONOCI / II) and the Nevada Geodetic Laboratory (NGL), covering the period from March 1, 2002 to December 31, 2021. This ensured that the data contained complete seasonal and interannual signal cycles. GNSS stations with a time series length of ≥2.5 years were selected (to avoid misinterpretation of seasonal signals caused by short-cycle data). Stations with a data missing rate >10% were removed, ultimately forming a GNSS station network covering the entire Tibetan Plateau.

[0086] (2) High-precision data processing and outlier identification and repair

[0087] The coordinate time series under the ITRF2014 framework was obtained using GipsyX software. The TSAnalyzer tool was used to remove step anomalies caused by instrument variations and regional tectonic activities (such as earthquakes and fault activity), and long-term trends were removed through linear fitting (preliminary separation of tectonic signals). Specifically, outliers were replaced using linear interpolation; step anomalies were estimated and removed using known instrument replacement records or earthquake event timestamps; and abrupt trend changes were corrected using piecewise linear fitting after verification with regional tectonic activities, ultimately resulting in a standardized GNSS vertical time series without anomalies.

[0088] (ii) Multi-source load stratification correction: systematically eliminating nonlinear signals

[0089] (1) AO correction

[0090] AO correction aims to remove the influence of non-tidal atmospheric load (NTAL) and non-tidal ocean load (NTOL) on GNSS vertical time series. The specific implementation is as follows:

[0091] 1. NTAL displacement calculation

[0092] Atmospheric pressure data with a 3-hour resolution (spatial resolution 0.25°×0.25°) provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) were used, with the time span consistent with GNSS data. Harmonic analysis was performed on 12 major atmospheric tidal components (such as S1, S2, P1, etc.) to remove tidal atmospheric effects and retain non-tidal components. Based on the elastic earth load theory, the NTAL model provided by GFZ was used to calculate the NTAL vertical displacement time series for each GNSS station using the load Green's function, as shown in the following formula:

[0093] ;

[0094] in, At time t Non-tidal atmospheric pressure anomaly at the location For the load Green's function ( (The angle is the geocentric angle), and the integration range covers the entire globe.

[0095] 2. NTOL displacement calculation

[0096] Adopting Max Three-hour resolution seafloor pressure data (spatial resolution 1°×1°) from the Max Planck Institute Ocean Model (MPIOM). The NTOL model based on GFZ also calculates the NTOL vertical displacement time series using the load Green's function, considering the boundary effects between the ocean and land, and employing refined interpolation for areas near the coastline to ensure the accuracy of displacement calculations.

[0097] 3. AO correction implementation

[0098] The calculated NTAL and NTOL vertical displacement time series are superimposed to obtain the total AO load displacement time series. The total AO load displacement time series is then subtracted from the standardized GNSS vertical time series to obtain:

[0099] ;

[0100] in, To standardize GNSS time series, This is the GNSS vertical time series after AO correction. After correction, the RMS values ​​in the southern monsoon region (26°N-32°N) decreased by 10%-20%, and the reduction in the Nepal region (27°N-30°N, 82°E-88°E) exceeded 20%, effectively eliminating interference from high-frequency atmospheric pressure fluctuations.

[0101] (2) Hydrological load correction

[0102] To address the complexity of the hydrological load on the Qinghai-Tibet Plateau, this embodiment employs two correction schemes: AOH correction (AO+GFZ-HYDL) and AOG correction (AO+GRACE-HYDL) for comparison, as detailed below:

[0103] 1. Displacement calculation of GFZ-HYDL

[0104] The hydrological load model provided by GFZ was used, based on the 24-hour resolution hydrological variables (soil moisture content, snow depth, etc.) of the Land Surface Discharge Model (LSDM). Considering the mass redistribution of soil water and snow cover, the vertical displacement of HYDL was calculated through elastic load theory (a surface elastic load model based on the graph center (CF) framework). The parameter settings were optimized, focusing on the lake area (<34°N) within the plateau, to match the spatiotemporal variation characteristics of soil water around the lakes. The formula is as follows:

[0105] ;

[0106] in, It's the load size. It is the Earth's elastic modulus. It is the distance from the load point to the ground surface point, and A represents the surface area of ​​the ground surface affected by the load.

[0107] 2. Implementation of AOH Correction

[0108] The GFZ-HYDL displacement time series is superimposed on the AO-corrected GNSS vertical time series to obtain the AOH-corrected GNSS vertical time series:

[0109] ;

[0110] After correction, the RMS values ​​in the southern part of the plateau and the inland lake area decreased by 20%-35%.

[0111] 3. GRACE data preprocessing

[0112] GRACEmascon data (monthly resolution, 300 km spatial resolution) from the Center for Space Research (CSR), Jet Propulsion Laboratory (JPL), and Goddard Space Flight Center (GSFC) at the University of Texas at Austin, spanning from April 2002 to December 2021, were acquired. High-frequency noise in the GRACE data was eliminated using a 30-day moving average, and leakage errors were eliminated through scale factor correction (based on the global hydrological model GLDAS-2.1). For small-scale glaciers in the western Tibetan Plateau, the neighborhood averaging method was used to improve the spatial resolution of the data.

[0113] 4. Calculation of Equivalent Water Height (EWH):

[0114] The GRACEmascon (mass-dense) data is converted to EWH to reflect changes in total terrestrial water storage (TWS), as shown in the following formula:

[0115] ;

[0116] in, At time t The gravity anomaly at the location The density of water, This is the acceleration due to gravity.

[0117] Displacement modeling: Based on the load Green's function, EWH is converted into GRACE-HYDL vertical displacement, incorporating the combined effects of glacier mass loss, surface water (lakes, rivers), and groundwater changes. The formula is as follows:

[0118] .

[0119] 5. AOG Correction Implementation

[0120] The GRACE-HYDL displacement time series is superimposed on the AO-corrected GNSS vertical time series to obtain a multi-source load-layered corrected GNSS vertical time series:

[0121] ;

[0122] After AOG (AO combined with GRACE modeling of hydrological load) correction, the RMS value of the Himalayan glacier-covered area (25°N-28°N, 80°E-90°E) decreased by 25%-35%, an additional 10% reduction compared to AOH correction. This effectively captures the deformation caused by glacier mass loss. The AOG method has shown excellent performance in areas with complex hydrodynamics, such as areas affected by glacier-monsoon interactions.

[0123] (III) CME Filter Optimization: Eliminating Spatial Correlation Errors

[0124] Common-mode error (CME) is a key spatial correlation error affecting the accuracy of GNSS networks on the Tibetan Plateau. This step designs a CME filtering algorithm based on empirical orthogonal function (EOF) decomposition to accurately extract and eliminate CME. The specific process is as follows:

[0125] (1) CME extraction: Spatial correlation signal separation based on EOF decomposition

[0126] GNSS vertical time series corrected with AOG ( ), select GNSS stations (≥30, ensuring uniform spatial coverage) for observation during the same period, and conduct observations on them. Performing EOF decomposition yields eigenvectors (spatial modes) and eigenvalues ​​(variance contributions), as shown in the following formula:

[0127] ;

[0128] in, This is a time series matrix of GNSS stations (n ​​stations × m time points). Let n be the spatial eigenvector matrix (n×n). It is an eigenvalue diagonal matrix (n×m). It is the time eigenvector matrix (m×m).

[0129] Based on the variance contribution of the eigenvalues, the top k spatial modes (usually k=1-3, with a cumulative variance contribution ≥50%) are selected as the spatial modes of CME. The product of the corresponding time eigenvector and the eigenvalue is the CME time series, as shown in the following formula:

[0130] ;

[0131] in, These are the first k spatial feature vectors. For the corresponding eigenvalues, These are the first k time feature vectors. In the Qinghai-Tibet Plateau, CME mainly manifests as a uniform displacement change across the entire region (the first spatial mode variance contributes over 40%), and is related to regional atmospheric disturbances and satellite orbital deviations.

[0132] (2) Implementation of CME filtering

[0133] The extracted CME time series is subtracted from the AOG-corrected GNSS time series to obtain the AOG_CME-corrected time series:

[0134] .

[0135] (iv) Dynamic adaptation of noise model

[0136] To address the regional heterogeneity of noise types on the Tibetan Plateau, a dynamic noise model adaptation method is designed based on the optimized Bayesian information criterion (BIC_tp) to select the optimal noise model for each GNSS station. The specific process is as follows:

[0137] (1) Define five typical noise models to cover the main noise types on the Qinghai-Tibet Plateau:

[0138] WN model: contains only white noise, and the parameter is the standard deviation of white noise (σ_WN).

[0139] WN+FN model: white noise + flicker noise, where the FN parameter is the spectral index (α_FN≈1) and the noise amplitude (A_FN).

[0140] WN+FN+RWN model: white noise + flicker noise + red noise, RWN parameters are spectral index (α_RWN≈2) and noise amplitude (A_RWN).

[0141] WN+GGM model: white noise + generalized Gaussian-Markov noise, with GGM parameters being the correlation time (τ_GGM) and noise amplitude (A_GGM).

[0142] WN+PL model: white noise + power-law noise, where PL parameters are the spectral index (1<α_PL<2) and the noise amplitude (A_PL).

[0143] (2) Calculation of BIC_tp criterion

[0144] GNSS time series corrected by AOG_CME Five noise models were used for fitting, and the BIC_tp value of each model was calculated using the following formula:

[0145] ;

[0146] in, For the amount of observation data, For the model residual variance, The number of model parameters (e.g., WN+FN model) =3: σ_WN, α_FN, A_FN). The smaller the BIC_tp value, the better the model fits the data.

[0147] AOG correction effectively maintained the proportion of the WN+FN category at 61.1% and the proportion of the WN+GGM category at 27.5%, and completely eliminated the influence of the WN+FN+RWN category. Compared with the GFZ model, the removal of the RWN component highlights the superiority of GRACE in reducing noise caused by hydrological loads, demonstrating GRACE's advantage in comprehensively integrating terrestrial water storage data.

[0148] After applying AOG and CME corrections, the distribution of the optimal noise model changed significantly, as shown in Table 1. The proportion of WN+FN decreased from 61.1% to 51.7%, while the proportion of WN+GGM increased to 32.9%. In particular, the WN+FN+RWN category was completely removed, while the proportion of WN+PL increased to 15.4%. This change indicates that CME filtering is significantly effective in mitigating spatial correlation errors and significantly simplifies the complexity of stochastic models in the Tibetan Plateau region.

[0149] Table 1. Optimal noise model identified by BIC_tp before and after surface load correction.

[0150] Noise Model WN WN+FN WN+FN+RWN WN+GGM WN+PL Raw 0.0% 63.1% 1.3% 18.8% 16.8% AOremoved 0.0% 61.1% 1.3% 27.5% 10.1% AOH removed 0.0% 60.4% 1.3% 25.5% 12.8% AOG removed 0.0% 61.1% 0.0% 27.5% 11.4% AOG_CMEremoved 0.0% 51.7% 0.0% 32.9% 15.4%

[0151] (3) Model selection

[0152] The model with the smallest BIC_tp value was selected as the optimal noise model for this site. The results show:

[0153] In the southern monsoon region (26°N-32°N), the optimal model is WN+GGM (accounting for 32.9%). Due to the reduction of spatial correlation noise after CME filtering, GGM (reflecting mid-to-low frequency correlation noise) becomes dominant. In the northern permafrost region (34°N-40°N), the optimal model is WN+PL (accounting for 15.4%). PL noise matches the low-frequency characteristics of permafrost dynamics. In the southeastern seismically active region (90°E-98°E, 26°N-30°N), the optimal model is WN+FN (accounting for 51.7%). FN noise is correlated with high-frequency tectonic activity. Figure 2 This is a schematic diagram illustrating the evaluation results of the optimal background noise model for GNSS time series.

[0154] (v) Vertical velocity and structural signal extraction

[0155] (1) Vertical velocity extraction

[0156] Based on the optimal noise model, the vertical velocity (v) and its uncertainty of the GNSS station are estimated using Hector software. Considering the time correlation of noise, the velocity parameters are optimized through maximum likelihood estimation (MLE), as shown in the following formula:

[0157] ;

[0158] in, Design a matrix (including a time trend item). This is the noise covariance matrix (constructed based on the optimal noise model). for .

[0159] Figure 3 and Figure 4 The effects of surface loading deformation and CME on the GNSS velocity field and uncertainty are represented, respectively. The initial GNSS velocity field shows an average velocity of approximately 0.457 mm / yr, with significant spatial variability (uncertainty: ~0.573 mm / yr). After atmospheric loading (AO) correction, the average velocity increases to 0.555 mm / yr, and the uncertainty decreases to 0.567 mm / yr, with monsoon-driven high-frequency elastic deformation noise being partially suppressed. Further introduction of the Atmospheric and Hydrological Integrated Correction Model (AOH) increases the average velocity to 0.656 mm / yr, with the uncertainty slightly increasing to 0.573 mm / yr. After applying GRACE-based hydrological loading (AOG) correction, the average velocity decreases to 0.335 mm / yr, more closely resembling the actual tectonic signal, while the uncertainty remains at 0.569 mm / yr. These results demonstrate that the comprehensive constraint of GRACE data on Total Water Storage (TWS) effectively improves the accuracy and physical consistency of the velocity field. Compared to the GFZ model, GRACE performs better in addressing the water-ice equalization effect and eliminating anthropogenic signals. Finally, after introducing CME filtering, the overall uncertainty was reduced to 0.416 mm / yr, a decrease of approximately 26.9% compared to the GNSS results without AOG correction, validating the significant effect of spatial filtering in reducing regionally correlated noise and improving the reliability of the velocity field.

[0160] (2) Construction of signal extraction and analysis

[0161] Figures 5 to 9 These represent the vertical velocity and uncertainty of the GNSS station under the WN+FN stochastic model after different surface load corrections, velocity field mappings, and uncertainty mappings. Figure 5 Corresponding to the original vertical velocity field, Figure 6The corresponding AO has been corrected. Figure 7 The corresponding AOH has been corrected. Figure 8 The corresponding AOG has been corrected. Figure 9 The values ​​are AOG+CME corrected; v and u represent the average velocity and uncertainty of all GNSS stations. Results show a significant and continuous uplift in the southern part of the plateau (up to +2.3±0.4 mm / yr), a more stable central region, and slight subsidence in the north (-0.8±0.3 mm / yr”). After AO load correction, noise is significantly reduced, with an average decrease of approximately 23% in uncertainty at stations in the southern region and a reduction of approximately 13% in uplift rate.

[0162] The interpretation of the velocity field was significantly improved by joint correction using AOH and AOG. Correction based on GFZAOH effectively reduced measurement uncertainty; while AOG correction using GRACE constraints demonstrated superior performance in all regions, particularly in identifying compression signals masked by glacial hydrological influences in the southern inland Tibetan Plateau (26°N–33°N), revealing a coherent uplift of approximately +1.2 mm / yr. In the transition zone (26°N–33°N), the average velocity deviation observed between AOH and AOG was approximately 0.3 mm / yr, indicating that GRACE water storage constraints play a crucial role in extracting the true crustal movement signals in the cryosphere-tectonic interface region. Further introduction of AOG+CME filtering reduced the overall regional uncertainty by 25.6% (average σ = 0.41 mm / yr).

[0163] Although preferred embodiments of this application have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments as well as all changes and modifications falling within the scope of this application.

[0164] Obviously, those skilled in the art can make various modifications and variations to this application without departing from the spirit and scope of this application. Therefore, if such modifications and variations fall within the scope of the claims of this application and their equivalents, this application also intends to include such modifications and variations.

Claims

1. A method for GNSS vertical deformation signal extraction based on multi-source load correction and CME filtering, characterized in that, The method includes the following steps: S1. Collect continuous GNSS observation data over a preset time span to ensure that the collected data contains complete seasonal and interannual signal cycles. After correcting outliers, construct a standardized GNSS vertical time series dataset without outliers. S2. Calculate the non-tidal atmospheric load vertical displacement time series and the non-tidal ocean load vertical displacement time series for each GNSS station using the load Green's function. Superimpose the non-tidal atmospheric load vertical displacement time series and the non-tidal ocean load vertical displacement time series to obtain the total AO load displacement time series. Subtract the total AO load displacement time series from the standardized GNSS vertical time series to obtain the AO-corrected GNSS vertical time series. S3. Based on the mass change time series collected by GRACE satellite, the equivalent water height reflecting the change of total land water storage is calculated, and the hydrological vertical displacement time series is obtained based on the load Green's function. The hydrological vertical displacement time series is superimposed on the AO-corrected GNSS vertical time series to obtain the multi-source load-layered corrected GNSS vertical time series. S4. Select GNSS stations observed in the same period, perform EOF decomposition on the GNSS vertical time series with multi-source load hierarchical correction, calculate the CME time series based on the eigenvectors and eigenvalues ​​obtained from the decomposition, and then subtract the CME time series from the AOG-corrected GNSS vertical time series to obtain the AOG_CME-corrected GNSS vertical time series. S5. Define multiple noise models according to the type of noise on the plateau. Fit the time series after AOG_CME correction in step S4 using multiple noise models. Select the optimal noise model for each GNSS station based on the fitting structure to achieve dynamic adaptation. S6 estimates the vertical velocity and its uncertainty for each GNSS station based on the optimal noise model.

2. The GNSS vertical deformation signal extraction method based on multi-source load correction and CME filtering according to claim 1, characterized in that, Step S1 further includes: Continuous GNSS observation data from the China Crustal Movement Observation Network and the Nevada Geodetic Laboratory were collected over a preset time span. GNSS stations with a time series length greater than a preset duration threshold were selected, and stations with a data missing rate greater than a missing rate threshold were removed to form a GNSS station network covering the entire Qinghai-Tibet Plateau. The coordinate time series under the ITRF2014 framework was obtained using GipsyX software. Linear interpolation was used to replace gross errors. For step outliers caused by instrument changes or regional tectonic activities, step parameters were estimated and removed using known instrument replacement records or regional tectonic activity timestamps. For abrupt trend changes, piecewise linear fitting was used for repair. Finally, a standardized GNSS vertical time series without outliers was formed.

3. The GNSS vertical deformation signal extraction method based on multi-source load correction and CME filtering according to claim 1, characterized in that, Step S2 further includes: Three-hour resolution atmospheric pressure data with a time span consistent with the GNSS vertical time series were collected. Harmonic analysis was performed on 12 major atmospheric tidal components to remove tidal atmospheric effects and retain non-tidal components. Based on the elastic Earth load theory, the non-tidal atmospheric load model provided by GFZ was used to calculate the non-tidal atmospheric load vertical displacement time series for each GNSS station through the load Green's function. : ； in, At time t Non-tidal atmospheric pressure anomaly at the location For the load Green's function, The geocentric angle is used, and the integration range covers the entire globe; Collect Max Using 3-hour resolution seafloor pressure data from the Max Planck Institute ocean model, and based on the GFZ non-tidal ocean loading model, the time series of vertical displacement of non-tidal ocean loading was calculated using the load Green's function. Among them, considering the boundary effect between the ocean and the land, a refined interpolation method is used to process the displacement data in the area near the coastline; The total AO load displacement time series is obtained by superimposing the non-tidal atmospheric load vertical displacement time series and the non-tidal ocean load vertical displacement time series. The total AO load displacement time series is then subtracted from the standardized GNSS vertical time series to obtain the AO-corrected GNSS vertical time series. : ； in, To standardize GNSS time series, This is the AO-corrected GNSS vertical time series.

4. The GNSS vertical deformation signal extraction method based on multi-source load correction and CME filtering according to claim 1, characterized in that, Step S3 further includes: The mass change time series data collected by the GRACE satellite were acquired, and the acquired mass change time series data were preprocessed. The preprocessing process included: using a 30-day moving average to eliminate high-frequency noise in the mass change time series data, and using scale factor correction to eliminate leakage errors contained therein. The preprocessed time series data of quality changes were converted into equivalent water height. : ； in, At time t The gravity anomaly at the location The density of water, It is the acceleration due to gravity; Based on the load Green's function, the equivalent water height is converted into hydrological vertical displacement using the following formula. : ； Vertical displacement of hydrology Superimposed on the AO-corrected GNSS vertical time series, a multi-source load-layered corrected GNSS vertical time series is obtained: 。 5. The GNSS vertical deformation signal extraction method based on multi-source load correction and CME filtering according to claim 1, characterized in that, Step S4 further includes: GNSS vertical time series with multi-source load hierarchical correction Several GNSS stations observed during the same period were selected for analysis. EOF decomposition yields spatial eigenvectors, temporal eigenvectors, and eigenvalues, as shown in the following formula: ； in, Let be an n×m GNSS station time series matrix, where n represents the number of GNSS stations and m represents the number of time points. Let n×n be the spatial eigenvector matrix. Let m be an n×m eigenvalue diagonal matrix. The time eigenvector matrix is ​​m×m; Based on the variance contribution of the eigenvalues, the top k spatial modes are selected as the spatial modes of CME, and their products with the corresponding time eigenvectors and eigenvalues ​​are used as the CME time series. : ； in, These are the first k spatial feature vectors. For the corresponding eigenvalues, These are the feature vectors of the first k time periods; The extracted CME time series is subtracted from the multi-source load-stratified GNSS vertical time series to obtain the AOG_CME-corrected GNSS vertical time series: 。 6. The GNSS vertical deformation signal extraction method based on multi-source load correction and CME filtering according to claim 1, characterized in that, Step S5 further includes: Five noise models are defined: the WN model containing only white noise with parameters equal to the standard deviation of white noise; the WN+FN model containing white noise and flicker noise with FN parameters equal to the spectral exponent and noise amplitude; the WN+FN+RWN model containing white noise, flicker noise, and red noise with RWN parameters equal to the spectral exponent and noise amplitude; the WN+GGM model containing white noise and generalized Gaussian-Markov noise with GGM parameters equal to the correlation time and noise amplitude; and the WN+PL model containing white noise and power-law noise with PL parameters equal to the spectral exponent and noise amplitude. AOG_CME corrected GNSS vertical time series for each GNSS station Five noise models were used for fitting, and the BIC_tp value of each noise model was calculated using the following formula: ； in, The amount of observation data from GNSS stations. The residual variance of the noise model. The number of parameters in the noise model; The noise model with the smallest BIC_tp value is selected as the optimal noise model for the GNSS site.

7. The GNSS vertical deformation signal extraction method based on multi-source load correction and CME filtering according to claim 1, characterized in that, In step S6, the velocity parameters are optimized using maximum likelihood estimation, as shown in the following formula: ； in, For a design matrix that includes a time trend item, This is the noise covariance matrix constructed based on the optimal noise model. GNSS vertical time series .

Images