GNSS time series analysis and modeling method considering variable amplitudes

Abstract

The present application relates to geophysics and geodetic data processing technical field, disclose the GNSS time series analysis and modeling method considering variable amplitude, the method comprises the following steps: S1, obtain the GNSS coordinate time series;S2, time series is modeled as the Gaussian process defined by mean function and covariance function;S3, construct the mean function describing long-term linear trend;S4, construct the product kernel by the periodic kernel and the non-periodic smooth kernel multiplication as the covariance function, to unify the periodicity and time-varying amplitude of the signal modeling;S5, the maximum likelihood estimation method is used to solve the hyperparameter in model;S6, the time series is decomposed by using the optimized model, and the time-varying amplitude periodic signal is obtained.The present application can integrally model the periodic signal and its time-varying amplitude by constructing the Gaussian process product kernel, so as to realize the accurate separation of signal component.

Description

Technical Field

[0001] This invention relates to the field of geophysics and geodesy data processing technology, specifically to a method for GNSS time series analysis and modeling that takes into account variable amplitude. Background Technology

[0002] GNSS technology, such as the Global Positioning System (GPS), has become a core tool for monitoring geodynamic phenomena such as crustal deformation, glacial isostatic adjustment, and sea-level changes. By conducting long-term monitoring of the coordinates of continuously operating GNSS stations worldwide, high-precision, high-temporal-resolution coordinate time series can be obtained. These time series contain rich geophysical information, primarily composed of a superposition of three parts: long-term linear trends (reflecting tectonic movements), periodic signals (mainly caused by changes in surface mass load), and noise. Therefore, accurately separating these signal components from the observational data is crucial for subsequent geophysical interpretation and scientific research.

[0003] In existing technologies, a common approach to processing such GNSS coordinate time series is to combine a deterministic function model with a random noise model. Specifically, the function model typically includes a linear or polynomial trend term describing the long-term stable motion of the station, and a set of harmonic terms (such as sine and cosine functions) describing periodic variations such as annual and semi-annual changes. The model parameters, including the linear rate, the amplitude and phase of the periodic terms, are generally estimated using the least squares method. This method is widely used due to its clear model structure and computational simplicity.

[0004] While existing deterministic function models can describe the main characteristics of GNSS time series to some extent, they still have some shortcomings. The core flaw of this method lies in its assumption that the amplitude of the periodic signal is constant. This assumption stems from the inherent limitations of the mathematical structure of the harmonic functions used, namely, that the amplitude parameter is a constant that does not change with time. However, in physical reality, the geophysical factors driving the periodic displacement of GNSS stations, such as seasonal hydrological loads, atmospheric loads, and non-tidal ocean loads, are not strictly the same year after year, but fluctuate interannually with factors such as climate change and extreme weather events. When a physical signal with inherently variable amplitude is forced to fit a mathematical model with constant amplitude, the amplitude variation that the model cannot capture will "leak" and remain in the fitted residual sequence. This not only leads to insufficient extraction of the physical characteristics of the periodic signal itself (such as the dynamic evolution of amplitude), but also contaminates the noise component, causing the noise level to be overestimated, and may introduce spurious periodicity into the residual, thus limiting the accuracy and reliability of subsequent geophysical research. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method for GNSS time series analysis and modeling that takes into account variable amplitude. This method solves the problem that existing technologies, which use function models with constant amplitude, cannot accurately model and extract periodic signals with time-varying amplitude characteristics from GNSS time series.

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] Considering the analysis and modeling methods for GNSS time series with variable amplitude, this method introduces a non-parametric model called Gaussian Process (GP). By designing specific mean and covariance functions, it achieves direct modeling and accurate decomposition of time-varying amplitude periodic signals in GNSS coordinate time series.

[0008] The technical solution of the present invention may include the following steps:

[0009] Includes the following steps:

[0010] S1. Obtain the time series of coordinates of the Global Navigation Satellite System;

[0011] S2. The acquired global satellite navigation system coordinate time series is modeled as a Gaussian process defined by the mean function and the covariance function;

[0012] S3. Construct a mean function to describe the long-term linear trend of the coordinate time series of the global satellite navigation system;

[0013] S4. Construct the covariance function. The covariance function takes the form of a product kernel. The product kernel is composed of a periodic kernel used to ensure the strict periodicity of the signal and an aperiodic smoothing kernel used to smooth the amplitude of the periodic signal.

[0014] S5. Integrate the constructed mean function and the constructed covariance function to form a complete Gaussian process model containing undetermined hyperparameters. Then, using the maximum likelihood estimation method, solve for the hyperparameters in the Gaussian process model based on the global satellite navigation system coordinate time series obtained in step S1 to obtain the optimal hyperparameters.

[0015] S6. Substitute the optimal hyperparameters obtained from the solution into the complete Gaussian process model to obtain an optimized Gaussian process model. Then, apply the optimized Gaussian process model to decompose the obtained global satellite navigation system coordinate time series to obtain the time-varying amplitude periodic signal.

[0016] Preferably, in step S1, obtaining the coordinate time series of the global satellite navigation system specifically involves obtaining the coordinate time series composed of the observation time t and the observation value vector Y corresponding to the observation time.

[0017] Preferably, in step S2, the global navigation satellite system coordinate time series is modeled as a Gaussian process, specifically as follows:

[0018] The generation process of the observation vector Y at any time t1 can be expressed as a linear superposition of a potential signal function f2(t1) as a Gaussian process and an independent and identically distributed observation noise ∈(t1), that is:

[0019] Y(t1) = f2(t1) + ∈(t1), where Y(t1) is the acquired coordinate time series; f2(t1) is the signal function at any time t1; and ∈(t1) is the observation noise at any time t1, assumed to be white noise following a zero-mean Gaussian distribution.

[0020] Preferably, in step S3, the mean function is constructed by setting the mean function to include only the long-term linear trend term, which is used to uniformly process the complex dynamic characteristics of the coordinate time series of the global satellite navigation system, including periodicity and its amplitude changes, by the constructed covariance function.

[0021] Preferably, in step S4, the non-periodic smoothing kernel is a Matern 3 / 2 kernel, and the periodic kernel is an exponential sine square kernel; wherein, the Matern 3 / 2 kernel is selected to introduce time dependence and model the smooth change of amplitude, and the exponential sine square kernel is selected to accurately capture the strictly repeating periodic features in the signal.

[0022] Preferably, in step S4, the complete expression for the covariance function is:

[0023]

[0024] In the formula, κ(t) i ,t j () represents any two observation times t i and t j The covariance between them; σ 2 is the signal variance, which controls the overall intensity or amplitude of the modulated periodic signal; l is the length scale parameter of the Matern 3 / 2 kernel, which determines the smoothness of the periodic signal amplitude change over time; p is the period length parameter of the periodic kernel; l p This is the smoothness parameter of the period kernel, which controls the smoothness of the waveform within the period; The white noise variance represents the level of random measurement error in the coordinate time series of a global navigation satellite system; δ ij Let t be the Kronecker function, when i = j. i With t j At the same moment, the Kronecker function δ ijThe value is 1, representing the variance of white noise. When i ≠ j, t is correctly added to the diagonal of the covariance matrix. i With t j For two different observation times, the Kronecker function δ ij The value is 0; |t i -t j | represents the time distance between two observation times; exp(·) is the natural exponential function; sin(·) is the sine function.

[0025] Preferably, in step S5, solving for the hyperparameters in the Gaussian process model includes the following steps:

[0026] First, we construct a marginal log-likelihood function for solving the hyperparameter φ in the Gaussian process model as the optimization objective function;

[0027] Then, a gradient-based numerical optimization algorithm is used to iteratively maximize the marginal log-likelihood function until convergence, thus obtaining the optimal hyperparameters.

[0028] Preferably, the expression for the marginal log-likelihood function is:

[0029]

[0030] In the formula, l(φ; Y) N ) represents the value of the marginal log-likelihood function; φ is the hyperparameter; Y N Let N be a column vector consisting of all observations; N is the total number of observation points; π is the mathematical constant pi; Σ N The N×N covariance matrix is ​​constructed from the covariance function based on all observation times; |Σ N | is the covariance matrix Σ N The determinant of; The covariance matrix Σ N The inverse matrix; Y is the column vector of observations N The transpose of .

[0031] Preferably, step S6, which involves decomposing the acquired global navigation satellite system coordinate time series to obtain a time-varying amplitude periodic signal, includes:

[0032] By applying the optimized Gaussian process model, the posterior mean of the signal function is calculated through posterior prediction. Based on this, the coordinate time series of the global satellite navigation system is accurately separated into three signal components with clear physical meanings: long-term linear trend, time-varying amplitude periodic signal, and residual noise component.

[0033] Preferably, the obtained time-varying amplitude periodic signal is used to reflect and quantify the dynamic evolution of the periodic signal amplitude in the time series of global satellite navigation system coordinates caused by seasonal hydrological load, glacier ablation environmental load changes, or slow-slip tectonic activity of faults.

[0034] The residual noise component is obtained by subtracting the posterior mean μ of the Gaussian process from the original observations. post The result (t) shows that this component represents the random fluctuations that the model failed to explain. The formula for this component is: R(t) = Y(t) - μ post (t), where R(t) is the value of the residual noise component extracted at observation time t; Y(t) is the original GNSS coordinate observation value obtained at observation time t; μ post f(t) represents the posterior mean of the Gaussian process at observation time t after all observation data have been given, and represents the optimal estimate of the signal function f2(t); t is the observation time.

[0035] This invention provides a method for GNSS time series analysis and modeling that takes into account variable amplitude. It has the following beneficial effects:

[0036] 1. This invention constructs a product kernel covariance function consisting of a periodic kernel and a non-periodic smooth kernel, enabling the model to directly describe the time-varying characteristics of the amplitude of a periodic signal. Compared to traditional fixed-amplitude models that attribute this variation to residuals, this invention can more realistically fit the non-stationary characteristics present in the observed data, thereby significantly reducing the residuals between the model and the observed values ​​and obtaining a more accurate overall fitting result.

[0037] 2. This invention simplifies the mean function to describe only the long-term linear trend, while entrusting all dynamic characteristics, such as periodicity and amplitude variation, to the covariance function of a Gaussian process. This specific model structure enables, after model optimization, the time-varying amplitude periodic signal, previously unidentifiable due to its mixing in the residuals, to be accurately separated from the original sequence as an independent, quantifiable component through posterior prediction.

[0038] 3. This invention can extract time-varying amplitude periodic signals as a definite signal component. This component is no longer a vague, unmodeled error, but a quantitative result that can be directly used for analysis. This makes it possible to compare and correlate the changes of this signal component with environmental factors such as seasonal hydrological loads and glacial mass migration, or geological tectonic activities such as slow slip, thereby transforming some signals that were originally considered noise into evidence with clear physical meaning that can be used for scientific research. Attached Figure Description

[0039] Figure 1 This is a schematic diagram of the system structure of the present invention;

[0040] Figure 2 This is a schematic diagram of the method flow of the present invention;

[0041] Figure 3 This is a comparison chart showing the signal fitting results of the three models of this invention on the simulated data;

[0042] Figure 4 This is a comparison chart showing the results of recovering the amplitude of analog signals using the three models of this invention;

[0043] Figure 5 This is the fitting result of the traditional model of this invention;

[0044] Figure 6 This is the fitting result of the GP residual model of the present invention;

[0045] Figure 7 This is the fitting result of the GP product kernel model of this invention;

[0046] Figure 8 This is a fitting residual plot of the traditional model of this invention;

[0047] Figure 9 This is a fitting residual plot of the GP residual model of the present invention;

[0048] Figure 10 This is a fitting residual plot of the GP product kernel model of this invention;

[0049] Figure 11 This is the amplitude sequence recovered by the GP residual model of this invention;

[0050] Figure 12 The amplitude sequence recovered by the GP product kernel model of this invention;

[0051] Figure 13 This is a time series diagram of the U-component coordinates of the LHAZ station in this invention. Detailed Implementation

[0052] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0053] Please see the appendix Figure 1 , Figure 1 This is a schematic diagram of a system structure according to an embodiment of the present invention. The present invention provides a system for time series analysis and modeling of a Global Navigation Satellite System (GNSS) considering variable amplitude, which may include:

[0054] The data input and preprocessing module 10 is used to obtain the original GNSS coordinate time series from external data sources and organize it into a standardized dataset containing observation times and corresponding observation value vectors.

[0055] The model building module 20 is used to construct a Gaussian process model defined by a mean function and a covariance function according to a preset technical solution. The mean function describes the long-term linear trend of the time series, while the covariance function, in the form of a product kernel, is used to uniformly model the periodicity of the signal and the time-varying characteristics of its amplitude.

[0056] The parameter optimization module 30, which is connected to the data input and preprocessing module 10 and the model building module 20, is used to receive the standardized dataset and model structure, and to use the maximum likelihood estimation method to execute the numerical optimization algorithm to solve the hyperparameters in the Gaussian process model.

[0057] The signal decomposition and output module 40, connected to the data input and preprocessing module 10 and the parameter optimization module 30, receives the standardized dataset and hyperparameters, applies the optimized Gaussian process model to decompose the time series, obtains the long-term linear trend, time-varying amplitude periodic signal, and residual noise component, and outputs the decomposition result. The residual noise component is obtained by subtracting the posterior mean μ of the Gaussian process from the original observations. post (t) is obtained.

[0058] Please refer to the appendix. Figure 2 , Figure 2 This is a schematic flowchart of a method according to an embodiment of the present invention. The present invention provides a method for time series analysis and modeling of global satellite navigation systems considering variable amplitude, which may include the following steps:

[0059] S1. Obtain the time series of coordinates of the Global Navigation Satellite System;

[0060] S2. The acquired global satellite navigation system coordinate time series is modeled as a Gaussian process defined by the mean function and the covariance function;

[0061] S3. Construct a mean function to describe the long-term linear trend of the coordinate time series of the global satellite navigation system;

[0062] S4. Construct the covariance function. The covariance function takes the form of a product kernel. The product kernel is composed of a periodic kernel used to ensure the strict periodicity of the signal and an aperiodic smoothing kernel used to smooth the amplitude of the periodic signal.

[0063] S5. Integrate the constructed mean function and the constructed covariance function to form a complete Gaussian process model containing undetermined hyperparameters, and use the maximum likelihood estimation method to solve for the hyperparameters in the Gaussian process model based on the global satellite navigation system coordinate time series obtained in step S1.

[0064] S6. Substitute the obtained hyperparameters into the complete Gaussian process model to obtain an optimized Gaussian process model. Then, apply the optimized Gaussian process model to decompose the obtained global satellite navigation system coordinate time series to obtain the time-varying amplitude periodic signal.

[0065] The specific implementation process of the method of the present invention will be described in detail below.

[0066] In step S1, the data input and preprocessing module 10 performs an operation to obtain a GNSS coordinate time series. This time series consists of a series of discrete observation times and their corresponding coordinate observation values.

[0067] In step S2, the model building module 20 performs an operation to establish a Gaussian process modeling framework. This framework represents the observed values ​​of the coordinate time series as a linear superposition of a potential signal function, which is a Gaussian process, and an independent and identically distributed observation noise, by subtracting the posterior mean μ of the Gaussian process from the original observed values. post (t) is obtained.

[0068] In step S3, the model building module 20 continues to perform operations to specifically construct the mean function and covariance function of the Gaussian process.

[0069] First, a mean function is constructed to describe the long-term linear trend of the coordinate time series. This mean function contains only one parameter describing the steady linear displacement rate. Next, a covariance function is constructed. This covariance function takes the form of a product kernel, consisting of multiplying a periodic kernel to ensure the strict periodicity of the signal with an aperiodic smoothing kernel to smooth the amplitude of the periodic signal. Through this structure, the covariance function can uniformly model the periodicity and the dynamic characteristics of its amplitude variation over time.

[0070] In step S4, the parameter optimization module 30 performs an operation to optimize and solve for the hyperparameters in the model. These hyperparameters include the linear rate of change parameter in the mean function, and parameters such as signal variance, length scale, period length, period smoothness, and white noise variance in the covariance function. This step uses the maximum likelihood estimation method, constructing the marginal log-likelihood function of the column vector composed of all observations, and applying a numerical optimization algorithm to determine the set of hyperparameters that maximize the signal function.

[0071] In steps S5 and S6, the signal decomposition and output module 40 performs operations. This module substitutes the hyperparameters obtained in step S4 into the Gaussian process model to obtain an optimized model. Subsequently, the optimized model is applied to decompose the original time series, and the long-term linear trend component, the time-varying amplitude periodic signal component, and the residual noise component are calculated and separated through the posterior prediction function of the Gaussian process.

[0072] The residual noise component is obtained by subtracting the posterior mean μ of the Gaussian process from the original observations. post The component (t) is obtained. This component represents the random fluctuations that the model fails to explain, and its formula is: R(t) = Y(t) - μ. post (t), where R(t) is the value of the residual noise component extracted at observation time t; Y(t) is the original GNSS coordinate observation value obtained at observation time t; μ post f(t) represents the posterior mean of the Gaussian process at observation time t after all observation data have been given, and represents the optimal estimate of the signal function f2(t); t is the observation time.

[0073] Step S1 of the present invention is performed by the data input and preprocessing module 10, which is used to acquire and prepare the Global Navigation Satellite System (GNSS) coordinate time series required for subsequent analysis.

[0074] This step specifically includes acquiring raw observation data using receiver equipment from one or more GNSS continuous observation stations. This raw observation data undergoes a standardized, high-precision data processing procedure. In one specific embodiment, this processing procedure employs specialized software such as GAMIT / GLOBK, based on a double-difference positioning model, and performs unified calculations within the International Geodetic Reference Frame (e.g., ITRF2014) to obtain high-precision station coordinate time series.

[0075] The obtained coordinate time series typically contains three independent components: a north-south (N) component, an east-west (E) component, and a vertical (U) component. The method of this invention can be applied independently to each of these component time series. Each component time series consists of a series of discrete observation times and the corresponding coordinate observation values ​​for each time point.

[0076] To facilitate subsequent mathematical modeling, the data input and preprocessing module 10 organizes and outputs a standardized dataset from the selected time series component. This dataset comprises two core parts:

[0077] An observation time t = [t1, t2, ..., t N ] T It stores the specific time information for all N observation times.

[0078] An observation vector Y = [y1, y2, ..., y3] N ] T It stores N coordinate observation values ​​that correspond one-to-one with each time point in the observation time t.

[0079] This standardized dataset serves as the direct input for all subsequent model building, parameter optimization, and signal decomposition steps.

[0080] Step S2 of the present invention is performed by the model building module 20. This step is used to establish a Gaussian Process (GP) modeling framework, that is, to determine the form of the mathematical model used to describe the intrinsic structure of the GNSS coordinate time series.

[0081] This step mathematically represents the coordinate time series Y(t1) obtained in step S1 as a linear superposition of a signal function f2(t) and an observation noise ∈(t1). The relationship is as follows:

[0082] Y(t1) = f2(t1) + ∈(t1), where Y(t1) is the acquired coordinate time series; f2(t1) is the signal function at any time t1; and ∈(t1) is the observation noise at any time t1, which is usually assumed to be white noise following a zero-mean Gaussian distribution.

[0083] Preferably, the signal function f2(t) is modeled as a Gaussian process. A Gaussian process is defined as a set of random variables, where any finite number of random variables have a joint Gaussian distribution. This process is entirely composed of a mean function μ2(t) and a covariance function κ(t). i ,t j The kernel function (or kernel function) determines the kernel function.

[0084] Therefore, the signal function f2(t) follows a Gaussian process distribution, denoted as:

[0085]

[0086] In the formula, f2(t) is the signal function at observation time t, which represents the real signal without observation noise; ~ is a mathematical symbol indicating that the random variable on the left follows the probability distribution specified on the right. κ(t) is an abbreviation for Gaussian Process, representing a stochastic process completely defined by a mean function and a covariance function; μ2(t) is the prior mean function of the Gaussian process, used to describe the presupposition of the expected value of the signal function f2(t) before any data is observed; κ(t) i ,t j) is the covariance function (or kernel function) of a Gaussian process, used to calculate the variance at two observation times t. i and t j The signal function f2(t) i f2(t) and f2(t) j Covariance between ) and t; i ,t j These are two specific observation times; t is the observation time.

[0087] The mean function μ2(t1) defines the expected value of the signal function f2(t) at any time t1, i.e., E[f2(t)] = μ2(t1). It describes the average behavior or trend of the function as a whole.

[0088] Covariance function κ(t) i ,t j The definition is given for any two different times t. i and t j The signal function f2(t) i f2(t) and f2(t) j The covariance between ) is, i.e., Cov(f2(t) i ),f2(t j ))=κ(t i ,t j Cov(·,·) is the covariance operator used to calculate the covariance between two random variables as input; this function encodes prior information about the characteristics of the function, such as smoothness, periodicity, and amplitude variation.

[0089] The observation noise ∈(t) is assumed to be independent and identically distributed Gaussian white noise with a mean of 0 and a variance of . The Gaussian distribution, i.e.

[0090] By establishing the above framework, this step transforms the time series analysis problem from fitting a deterministic parameterized function to inferring posterior probabilities within a function space defined by a Gaussian process. This framework lays the foundation for the detailed design of the mean and covariance functions in subsequent steps S3 and S4, as well as for modeling time-varying amplitude signals.

[0091] In the process of building the Gaussian process model, model building module 20 specifically performs the construction of the mean function. This step aims to determine the functional form used to describe the overall average behavior of the GNSS coordinate time series.

[0092] Preferably, the mean function μ2(t) is constructed as a simplified form containing only a long-term linear trend term. The specific mathematical expression of this mean function is:

[0093] μ2(t) = vt;

[0094] In the formula, μ2(t) is the value of the mean function at the observation time t, which is used to describe the long-term linear trend of the time series; v is the linear rate of change parameter, which is used to describe the long-term, stable linear displacement rate caused by factors such as plate movement; and t is the observation time.

[0095] A key technical objective of this construction method is to separate the complex dynamic characteristics of GNSS time series, especially the time-varying nature of periodic signals and their amplitudes, from the modeling scope of the mean function.

[0096] By limiting the mean function to handle only linear trends, all nonlinear and dynamic signal features are uniformly described and modeled by the covariance function constructed in step S4. This division of responsibilities simplifies the overall structure of the model and lays the foundation for directly and explicitly handling time-varying amplitude problems through the covariance function. The parameter v, as an undetermined hyperparameter, will be optimized and solved together with other hyperparameters in step S5.

[0097] In the process of building the Gaussian process model, the model building module 20 specifically performs the covariance function construction operation. This step is used to mathematically describe the dynamic characteristics of all signals except for long-term linear trends, and is a key step in realizing the direct modeling of time-varying amplitude periodic signals.

[0098] Preferably, the covariance function K2 takes the form of a product kernel. This structure consists of multiplying a periodic kernel and a non-periodic smoothing kernel. The technical principle behind this design is that the periodic kernel describes the strictly repeating periodic characteristics of the signal, while the non-periodic smoothing kernel modulates the amplitude of the periodic signal, describing the smooth evolution of the amplitude itself over time. By multiplying the two, a unified and explicit model of a periodic signal whose amplitude is modulated by a time function can be achieved.

[0099] Specifically, the period kernel can be an exponential sine squared kernel. This kernel function sets the time interval for signal repetition (e.g., an annual or semi-annual period) through the period length parameter p, and the smoothness parameter l... p By controlling the smoothness of the waveform within a period, the periodicity of the signal can be accurately captured.

[0100] The non-periodic smoothing kernel can be the Matern 3 / 2 kernel. This kernel function controls the smoothness of the amplitude function's change over time through the length scale parameter l. The Matern 3 / 2 kernel is chosen based on its first-order mean-square differentiability, which is consistent with the physical reality of amplitude smoothing changes caused by factors such as environmental loads in geophysical processes.

[0101] Multiplying the periodic kernel by the aperiodic smoothing kernel and then adding the term representing the independent observation noise from step S2 yields the complete covariance function K2. This function is used to calculate the covariance function K2 for any two observation times t. i and t j Covariance κ(t) between i ,t j The specific mathematical expression is:

[0102]

[0103] In the formula, κ(t) i ,t j () represents any two observation times t i and t j The covariance between them; σ 2 is the signal variance, which controls the overall strength or amplitude of the modulated periodic signal; l is the length scale parameter of the Matern 3 / 2 kernel, which determines the smoothness of the periodic signal amplitude change over time; p is the period length parameter of the periodic kernel, used to set the time interval for signal repetition, such as an annual period or a semi-annual period; p This is the smoothness parameter of the period kernel, which controls the smoothness of the waveform within the period; δ represents the variance of white noise, indicating the level of random measurement error in the observed data; ij Let t be the Kronecker function, when i = j. i With t j At the same moment, the Kronecker function δ ij The value is 1, representing the variance of white noise. It is correctly added to the diagonal of the covariance matrix; when i ≠ j, t i With t j For different times, δ ij The value is 0; |t i -t j | represents the time distance between two observation times; exp(·) is the natural exponential function; sin(·) is the sine function.

[0104] The parameter set of the covariance function constructed in this way All of these are undetermined hyperparameters, which will be optimized and solved in step S5.

[0105] Step S4 of the present invention is performed by the parameter optimization module 30. The purpose of this step is to solve the Gaussian process model containing unknown parameters constructed in steps S3 and S4 to determine a set of hyperparameters.

[0106] This step first identifies all hyperparameters that need optimization. These hyperparameters are derived from the mean function constructed in step S3 and the covariance function constructed in step S4. In one specific implementation, the hyperparameters φ to be optimized include: the linear rate of change parameter v from the mean function, and the signal variance σ from the covariance function. 2 Length scale parameter l, period length parameter p, period smoothness parameter l p and white noise variance Therefore, the complete set of hyperparameters is

[0107] To solve for this set of hyperparameters, this invention employs the Maximum Likelihood Estimation (MLE) method.

[0108] Specifically, the parameter optimization module 30 uses the column vector Y composed of all observations. N The marginal log-likelihood function is used to determine the hyperparameter φ. This marginal log-likelihood function l(φ; Y) N The mathematical expression for ) is:

[0109]

[0110] In the formula, l(φ; Y) N ) represents the value of the marginal log-likelihood function; φ is the hyperparameter; Y N Let N be a column vector consisting of all observations; N is the total number of observation points; π is the mathematical constant pi; Σ N The N×N covariance matrix is ​​constructed from the covariance function based on all observation times; |Σ N | is the covariance matrix Σ N The determinant of; The covariance matrix Σ N The inverse matrix; Y is the column vector of observations N The transpose of .

[0111] The parameter optimization module 30 searches for the marginal log-likelihood function l(φ; Y) by executing a gradient-based numerical optimization algorithm, such as the L-BFGS-B algorithm or the conjugate gradient algorithm. N The hyperparameter φ that reaches its maximum value.

[0112] The final output of this step is a set of hyperparameters. These values ​​will be passed to the signal decomposition and output module 40 to construct a fully deterministic, optimized Gaussian process model for subsequent signal decomposition operations.

[0113] Step S5 of the present invention is executed by the signal decomposition and output module 40. This module receives the hyperparameters output by the parameter optimization module 30 and the raw observation data provided by the data input and preprocessing module 10. Its purpose is to apply the optimized Gaussian process model to accurately decompose the original time series and output the results.

[0114] The core of this step is utilizing the posterior prediction function of Gaussian process regression. Specifically, the signal decomposition and output module 40 uses the hyperparameters obtained in step S5 to construct a fully deterministic Gaussian process model. Subsequently, based on this model and the original observation data, the posterior distribution of the signal function f2(t) is calculated. This posterior distribution is itself also a Gaussian process, with the optimal estimate (posterior mean) μ. post f(t) is the posterior mean of the Gaussian process at observation time t after all observation data have been given, representing the optimal estimate of the signal function f2(t).

[0115] To obtain the posterior mean μ of the signal function post After (t), the signal decomposition and output module 40 performs the following operation to decompose the original time series into three independent signal components:

[0116] Extracting the long-term linear trend component: This component is directly given by the mean function μ2(t1) constructed in step S3 and optimized in step S5. Its value at any time t1 is given by v. opt ·t is calculated, where v opt This is the optimal linear rate of change parameter obtained from parameter optimization module 30. This component represents the overall linear evolution trend of the time series.

[0117] Extracting the time-varying amplitude periodic signal component: This component is obtained by extracting the posterior mean μ of the Gaussian process. post The result is obtained by subtracting the long-term linear trend component μ2(t) from (t). Since the posterior mean μ post f2(t) is the best estimate of the entire signal function f2(t), which includes the linear trend and all dynamic characteristics. Therefore, after subtracting the linear trend, the result is the pure time-varying amplitude periodic signal component.

[0118] The formula for calculating this component is: S periodic (t)=μ post (t)-μ2(t), where S periodic (t) represents the value of the time-varying amplitude periodic signal component extracted at observation time t; μ postf2(t) represents the posterior mean of the Gaussian process at observation time t, given all the observed data, and is the optimal estimate of the signal function f2(t); μ2(t) is the value of the prior mean function at observation time t, representing the long-term linear trend component of the time series; t is the observation time.

[0119] Finally, the signal decomposition and output module 40 outputs the three separated time series components (long-term linear trend, time-varying amplitude periodic signal, and residual noise component) as the final result of this invention. These decomposed components can be stored, visualized, or used for subsequent geophysical analysis and research.

[0120] To verify the performance of the GNSS time series analysis and modeling method considering variable amplitude provided in this invention (hereinafter referred to as the GP product kernel model), a set of simulated data experiments were conducted and compared with two existing technical methods. The two existing technical methods are: a traditional model based on the assumption of constant amplitude, and a GP residual model that models the residuals using a Gaussian process based on the traditional model.

[0121] Example 1:

[0122] Experimental Design and Objectives: The objective of this simulation experiment is to quantitatively evaluate the overall fitting accuracy and dynamic characteristic capture capability of the method of the present invention when processing periodic signals with time-varying amplitude characteristics.

[0123] To achieve this goal, a set of simulated observation data was first constructed. This simulated data was generated by superimposing a linear trend term, a time-varying amplitude periodic term, and a residual noise component. Specifically:

[0124] The linear trend term is set as a linear function with a fixed rate.

[0125] The time-varying amplitude period term is set as a sine function with a period of one year, and its amplitude is set to fluctuate and increase over time to simulate the dynamic changes in amplitude caused by factors such as environmental loads in real crustal movements.

[0126] The residual noise component is a linear superposition of white noise and flicker noise with fixed parameters.

[0127] The simulated data was processed using the method of this invention and two comparative methods.

[0128] For the traditional model, the least squares method is used to fit the parameters of a function model containing a linear trend term and a fixed-amplitude periodic term composed of sine and cosine harmonic terms. Specifically, for an annual periodic signal, the mathematical form of this periodic term can be expressed as A·sin(2πt) + B·cos(2πt), where t is time in years; A is the constant amplitude parameter of the sine harmonic term to be solved, whose value does not change with time; B is the constant amplitude parameter of the cosine harmonic term to be solved, whose value does not change with time; and π is pi.

[0129] For the GP residual model, the same function model as the traditional model is first used for fitting, and then a Gaussian process with a Matern 3 / 2 kernel as the covariance function is applied to model the fitting residual.

[0130] The GP product kernel model provided by this invention is directly processed using the methods in steps S2 to S5, that is, the mean function only contains a linear trend term, and the covariance function adopts the product form of the Matern 3 / 2 kernel and the period kernel.

[0131] For both the GP residual model and the GP product kernel model, their hyperparameters are optimized by maximizing the marginal likelihood.

[0132] Performance evaluation method: In order to quantitatively evaluate the performance of the three models, this embodiment uses two evaluation indicators.

[0133] The first metric is the Root Mean Square Error (RMSE), which quantifies the accuracy of the model's fit to the overall data. Its calculation formula is:

[0134]

[0135] In the formula, RMSE is the root mean square error, which is used to quantify the overall deviation between the model fit value and the true signal value; N1 is the total number of data points in the time series; i is the index of the data point, which ranges from 1 to N1. y represents the fitted or predicted value of the model at the i-th observation time; i Let be the true signal value at the i-th observation time. In the simulation experiment, this value is a pre-set signal value that does not contain residual noise components.

[0136] The second metric is the amplitude correlation coefficient, used to evaluate the model's ability to capture time-varying amplitude dynamics. This metric is calculated using the Pearson correlation coefficient, and its specific formula is as follows:

[0137]

[0138] In the formula, ρ(A) true A est ) represents the amplitude correlation coefficient, used to evaluate the degree of linear correlation between the amplitude sequence recovered by the model and the true amplitude sequence, and its value ranges from [-1, 1]; A true The true amplitude sequence is a pre-defined, baseline amplitude sequence that changes over time in a simulation experiment; A est To estimate the amplitude sequence, i.e., the sequence of amplitude changes over time recovered or estimated from simulated observation data by the evaluated model; Cov(A true A est (A) represents the true amplitude sequence. true With the estimated amplitude sequence A est Covariance between them; For the true amplitude sequence A true Standard deviation; To estimate the amplitude sequence A est The standard deviation.

[0139] Experimental Results and Analysis:

[0140] Please refer to the appendix. Figure 3 , Figure 3 This is a comparison chart of the results of fitting simulated data using three different models. In the chart, the real signal curve exhibits periodic fluctuations with time-varying amplitude. While the fitted curves of the traditional models can reflect the overall linear trend and periodicity, they deviate from the real signal curve in the details of amplitude changes over time. Compared to the traditional models, the fitted curve of the GP residual model shows improved local fit to the real signal curve. The fitted curve of the GP product kernel model provided by this invention has the highest fit to the real signal curve in both the overall trend and the details of amplitude changes.

[0141] Please refer to the appendix. Figure 4 , Figure 4 This is a comparison chart of the amplitude recovery results of three models for analog signals. In the chart, the true amplitude curve exhibits a pre-defined periodic change. The amplitude recovered by the traditional model is represented by a horizontal line, failing to reflect the dynamic changes of the true amplitude. The amplitude curve recovered by the GP residual model reflects the dynamic changes of the amplitude to some extent, but still differs in shape from the true amplitude curve. The amplitude curve recovered by the GP product kernel model provided by this invention has the highest degree of agreement with the true amplitude curve in both overall trend and detail.

[0142] Please refer to the table below, which lists the quantitative performance evaluation results of the three models in the simulation experiment.

[0143] Performance evaluation table for three models

[0144] Model RMSE ρ Traditional model 0.4013 0.0000 GP residual model 0.2995 0.7818 GP product kernel model 0.2072 0.8974

[0145] In terms of overall fitting accuracy, the RMSE value of the GP product kernel model of this invention is 0.2072, which is lower than that of the traditional model (0.4013) and the GP residual model (0.2995). Regarding the ability to capture the dynamic characteristics of time-varying periodic terms, the amplitude correlation coefficient of the GP product kernel model of this invention is 0.8974, which is higher than that of the traditional model (0.0000) and the GP residual model (0.7818), and closer to 1. These quantitative evaluation results demonstrate that the method of this invention possesses higher fitting accuracy and amplitude recovery accuracy when processing periodic signals with time-varying amplitudes.

[0146] Example 2:

[0147] To further verify the application effect of the technical solution provided by this invention in processing real Global Navigation Satellite System (GNSS) data, this invention selected the coordinate time series of LHAZ stations as experimental data. This embodiment compares and analyzes the method of this invention (hereinafter referred to as the GP product kernel model) with two existing technical methods (traditional model and GP residual model).

[0148] Please refer to the appendix. Figure 13 , Figure 13 This is a time series plot of the U-component coordinates of the LHAZ station according to an embodiment of the present invention. The time series data was collected from 2010 to 2022. As can be seen from the figure, the discrete observation points of this time series exhibit significant periodic fluctuations, and the amplitude of these fluctuations changes over time, showing significant and time-varying annual periodic signal characteristics caused by factors such as surface mass load.

[0149] Please refer to the appendix. Figure 5 , Figure 6 and Figure 7 This set of figures shows the fitting results of three models to the U component time series of the LHAZ station.

[0150] Please refer to the appendix. Figure 5 The figure shows the fitting result of the traditional model. Its fitting curve is a periodic curve with a constant amplitude, which fails to match the amplitude variation in the observed data.

[0151] Please refer to the appendix. Figure 6 The figure shows the fitting result of the GP residual model. Its fitting curve corrects the fitting result of the traditional model in local details, but there are still deviations from the overall shape of the observed data.

[0152] Please refer to the appendix. Figure 7The figure shows the fitting result of the GP product kernel model provided by this invention. Its fitting curve exhibits a high degree of morphological consistency with the discrete point distribution of the original observation data, especially at the peaks and troughs of the annual signal, where the model's fitting curve changes with the fluctuations in the observation data.

[0153] Please refer to the appendix. Figure 8 , Figure 9 and Figure 10 This set of figures compares the fitting residuals of the U component time series of the LHAZ station using three different models.

[0154] Please refer to the appendix. Figure 8 This figure shows the fitting residuals of a traditional model. The residual sequence exhibits a clear periodic structure that was not fully extracted by the model.

[0155] Please refer to the appendix. Figure 9 The figure shows the fitted residuals of the GP residual model. Compared with the traditional model, the amplitude of its residual sequence is reduced, but there is still some systematic fluctuation.

[0156] Please refer to the appendix. Figure 10 The figure shows the fitting residuals of the GP product kernel model of this invention. The amplitude of the residual sequence is the smallest among the three models, and its distribution is closer to random noise.

[0157] Please refer to the appendix. Figure 11 and Figure 12 The figures in this set are time series diagrams of the annual signal amplitude recovered by the GP residual model and the GP product kernel model of this invention.

[0158] Please refer to the appendix. Figure 11 The figure shows the amplitude sequence recovered by the GP residual model. The sequence exhibits some fluctuations, but its shape is rather irregular.

[0159] Please refer to the appendix. Figure 12 The figure shows the amplitude sequence recovered by the GP product kernel model of this invention. This sequence clearly demonstrates a dynamic process with a specific pattern of change. This result indicates that the covariance function construction method provided by this invention can directly and effectively extract the dynamic information of the amplitude evolution of periodic signals over time from the original GNSS time series.

[0160] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for analyzing and modeling GNSS time series with variable amplitude, characterized in that, Includes the following steps: S1. Obtain the time series of coordinates of the Global Navigation Satellite System; S2. The acquired global satellite navigation system coordinate time series is modeled as a Gaussian process defined by the mean function and the covariance function; S3. Construct a mean function to describe the long-term linear trend of the coordinate time series of the global satellite navigation system; S4. Construct the covariance function. The covariance function takes the form of a product kernel. The product kernel is composed of a periodic kernel used to ensure the periodicity of the signal and an aperiodic smoothing kernel used to smooth the amplitude of the periodic signal. S5. Integrate the constructed mean function and the constructed covariance function to form a complete Gaussian process model containing undetermined hyperparameters. Then, using the maximum likelihood estimation method, solve for the hyperparameters in the Gaussian process model based on the global satellite navigation system coordinate time series obtained in step S1 to obtain the optimal hyperparameters. S6. Substitute the optimal hyperparameters obtained from the solution into the complete Gaussian process model to obtain an optimized Gaussian process model. Then, apply the optimized Gaussian process model to decompose the obtained global satellite navigation system coordinate time series to obtain the time-varying amplitude periodic signal. In step S3, the specific mathematical expression of the mean function is: ； In the formula, The mean function at observation time The value of is used to describe the long-term linear trend of the time series; It is a linear rate of change parameter used to describe the long-term, stable linear displacement rate caused by factors such as plate tectonics; The observation time; In step S4, the non-periodic smoothing kernel is the Matern 3 / 2 kernel, and the periodic kernel is the exponential sine square kernel; the Matern 3 / 2 kernel is selected to introduce time dependence and model the smooth change of amplitude, and the exponential sine square kernel is selected to capture the repetitive periodic features in the signal. In step S4, the complete expression for the covariance function is: ； In the formula, For any two observation times and The covariance between them; The signal variance; The length scale parameter of the Matern3 / 2 core; The period length parameter is the periodicity of the periodic kernel; The smoothness parameter of the periodic kernel; The variance is the white noise variance. The Kronecker function; The time distance between two observation times; It is a natural exponential function; It is a sine function.

2. The method for GNSS time series analysis and modeling considering variable amplitude according to claim 1, characterized in that, In step S1, the time series of coordinates of the Global Navigation Satellite System is obtained, specifically: the time series of coordinates of the Global Navigation Satellite System is obtained from the observation time. The vector of observation values ​​corresponding to the observation time The coordinate time series that are formed together.

3. The method for GNSS time series analysis and modeling considering variable amplitude according to claim 1, characterized in that, In step S2, the time series coordinates of the global navigation satellite system are modeled as a Gaussian process, specifically as follows: The observation vector At any time The generation process can be expressed as a latent signal function that is a Gaussian process. With an independent and identically distributed observation noise The linear superposition, that is: In the formula, For obtaining the coordinate time series; For any time The signal function; For any time Observation noise.

4. The method for GNSS time series analysis and modeling considering variable amplitude according to claim 1, characterized in that, In step S3, the mean function is constructed. Specifically, the mean function is set to contain only a long-term linear trend term. This is used to uniformly process the complex dynamic characteristics of the coordinate time series of the global satellite navigation system, including periodicity and its amplitude changes.

5. The method for GNSS time series analysis and modeling considering variable amplitude according to claim 1, characterized in that, In step S5, solving for the hyperparameters in the Gaussian process model includes the following steps: First, we construct a hyperparameter model for solving the Gaussian process. The marginal log-likelihood function is used as the optimization objective function; Then, a gradient-based numerical optimization algorithm is used to iteratively maximize the marginal log-likelihood function until convergence, thus obtaining the optimal hyperparameters.

6. The method for GNSS time series analysis and modeling considering variable amplitude according to claim 5, characterized in that, The expression for the marginal log-likelihood function is: ； In the formula, The value of the marginal log-likelihood function; For hyperparameters; It is a column vector consisting of all the observations; This represents the total number of observation data points. Pi is a mathematical constant. The covariance function is constructed based on all observation times. × Covariance matrix; Covariance matrix The determinant; Covariance matrix The inverse matrix; Observation column vector The transpose of .

7. The method for GNSS time series analysis and modeling considering variable amplitude according to claim 1, characterized in that, Step S6, which involves decomposing the acquired global navigation satellite system coordinate time series to obtain a time-varying amplitude periodic signal, includes the following steps: By applying the optimized Gaussian process model, the posterior mean of the signal function is calculated through posterior prediction. Based on this, the coordinate time series of the global satellite navigation system is separated into three signal components with clear physical meanings: long-term linear trend, time-varying amplitude periodic signal, and residual noise component.

8. The method for GNSS time series analysis and modeling considering variable amplitude according to claim 7, characterized in that, The time-varying amplitude periodic signal is used to reflect and quantify the dynamic evolution of the periodic signal amplitude in the coordinate time series of the global satellite navigation system caused by changes in seasonal hydrological load, glacier ablation environmental load, or slow-slip tectonic activity.

Images