A filtering method and system for establishing a nonlinear earth reference frame

Abstract

The present invention belongs to the field of geodesy and geodynamics technology, and discloses a filtering method and system for establishing a nonlinear earth reference frame. The present invention performs noise analysis on the frame point coordinate time series to obtain noise analysis result information, and introduces a first-order autoregressive model to model the frame point time-related noise during noise analysis; pre-processes the frame point coordinates to obtain unconstrained frame point coordinate observation values; performs parameter solution based on square root information filtering according to the noise analysis result information and the unconstrained frame point coordinate observation values, obtains parameter solution result information, and stores the correlation matrix as a temporary file; performs square root information smoothing processing according to the parameter solution result information and the temporary file, obtains the optimal estimate of the parameter to be estimated for each epoch, and completes the establishment of a nonlinear earth reference frame. The present invention can more accurately describe the nonlinear motion of the measuring station and improve the accuracy of the earth reference frame.

Description

Technical Field

[0001] The present invention belongs to the technical field of geodesy and geodynamics, and more specifically, relates to a filtering method and system for establishing a nonlinear earth reference frame. Background Art

[0002] The terrestrial reference frame (TRF) is the basis for describing the shape and changes of the earth and expressing geospatial information. It is also a key geospatial information infrastructure for expanding human activities and promoting social development. TRF consists of a series of reference stations fixed to the surface of the earth, which are specifically manifested as the high-precision coordinates and changes of the reference stations at a certain epoch under the reference system. Affected by geophysical effects or factors such as earthquakes, the positions of the reference stations show nonlinear changes, which leads to the instantaneous accuracy of the TRF station coordinates being only at the centimeter level, which cannot meet the needs of high-precision earth science research applications such as sea level change monitoring. How to improve the accuracy of establishing and maintaining the accuracy of the earth reference frame has become a key issue that needs to be urgently solved in geodesy and geodynamics.

[0003] At present, the existing technology usually adopts the method of establishing a nonlinear earth reference frame based on least squares or Kalman filtering. The nonlinear earth reference frame based on least squares expands the parameters of the traditional linear frame, uses harmonic functions to model the annual / semi-annual changes, or uses environmental load models to correct nonlinear changes. However, the ability of this method to capture nonlinear changes in the station is limited. The nonlinear earth reference frame based on filtering technology, that is, the Kalman filtering method is used to parameterize the frame with coordinate time series, and the expression form of this type of framework has undergone fundamental changes. In the context of millimeter-level earth science research, the nonlinear framework based on the filtering method has its natural advantages in the accurate expression of nonlinear motion of the station and geocentric motion, EOP prediction and computational efficiency. The accuracy of the state equation and random model is one of the main influencing factors for establishing a nonlinear earth reference frame based on filtering technology. Flicker noise is considered to be the best model to describe the time-correlated noise in the coordinate time series of GNSS stations. Since flicker noise is very difficult to implement in filtering, a series of problems need to be overcome. The existing filtering state equation usually uses a random walk model to approximate flicker noise. Therefore, due to the insufficient precision of the modeling of the time-correlated noise at the frame points, the existing technology has a technical problem of limited ability to capture nonlinear changes of the measuring stations when establishing the earth reference frame. Summary of the invention

[0004] The present invention solves the problem in the prior art that the nonlinear changes of frame points cannot be accurately described when establishing the earth reference frame, resulting in poor accuracy of the earth reference frame, by providing a filtering method and system for establishing a nonlinear earth reference frame.

[0005] The present invention provides a filtering method for establishing a nonlinear earth reference frame, comprising the following steps:

[0006] Perform noise analysis on the frame point coordinate time series to obtain noise analysis result information; in the noise analysis, a first-order autoregressive model is introduced to model the frame point time-related noise;

[0007] Preprocess the frame point coordinates to obtain unconstrained frame point coordinate observation values;

[0008] According to the noise analysis result information and the unconstrained frame point coordinate observation values, parameter solution is performed based on square root information filtering to obtain parameter solution result information, and a correlation matrix is ​​stored as a temporary file;

[0009] According to the parameter solution result information and the temporary file, square root information smoothing processing is performed to obtain the optimal estimate of the parameter to be estimated at each epoch, and the establishment of the nonlinear earth reference frame is completed.

[0010] Preferably, the noise analysis of the framework point coordinate time series includes: taking the framework point coordinates in the local coordinate system of the measuring station as input data, adopting a noise model combination of variable white noise and a first-order autoregressive model, and performing noise analysis using a maximum likelihood estimation method to obtain noise analysis results, wherein the noise analysis results include state transfer coefficients of the first-order autoregressive process of all measuring stations, the standard deviation of the driving white noise, and the coefficient factor of the variable white noise; converting the noise analysis results from the local coordinates of the measuring station to the geocentric Cartesian coordinates to obtain the noise analysis result information.

[0011] Preferably, before performing the noise analysis, the method further includes: performing a process of eliminating gross errors on the input data by least square fitting.

[0012] Preferably, the preprocessing of the framework point coordinates includes: eliminating a priori constraints attached to the framework point coordinates according to a constraint matrix.

[0013] Preferably, the parameter solution based on square root information filtering includes: measurement update and time update;

[0014] The measurement update includes: establishing an observation equation, parameterizing the unconstrained frame point coordinate observation value by epoch-by-epoch station coordinate decomposition, solving the measurement update parameters, and obtaining the information matrix after the measurement update; and using the current measurement updated information matrix as the prior information for the next epoch time update;

[0015] The time update includes: introducing a first-order autoregressive model to model the time-related noise of the frame point, using the state equation to predict the parameters to be estimated, and obtaining the information matrix after the time update; and using the information matrix after the current time update as the prior information for the next epoch measurement update;

[0016] The parameters to be estimated include the reference epoch coordinates of each epoch, the linear velocity of the frame point, the time-correlated noise parameter, the seasonal signal anniversary item and the seasonal signal semi-anniversary item.

[0017] Preferably, when establishing the observation equation, the unconstrained frame point coordinate observation value is parameterized as follows:

[0018]

[0019] In the formula, t epoch as the current epoch, for t The unconstrained frame point coordinate observations for the epoch, for t the reference epoch coordinates of the epoch, for t Epoch and reference epoch The time interval for t The linear velocity of the epoch frame point, for t The time-dependent noise parameters of the epoch, for t The annual term of the seasonal signal at the epoch, for t The semi-annual term of the seasonal signal of the epoch, They are t The translation, scaling and rotation parameters between the epoch combined reference frame and the known reference frame, for t The coordinates of the frame points of the epoch combined reference frame, is the coefficient factor of variable white noise, for t variable white noise term for the epoch;

[0020] During the measurement update process, the observation equation includes the earth orientation parameter observation equation, which is expressed as follows:

[0021]

[0022] In the formula, They are t Observed values ​​of the two components of the epoch polar motion parameter; They are t Observed values ​​of the polar motion rate of the two components of the epoch; for t Observations of changes in the Earth's rotation rate over the epoch; They are t Two components of the polar motion parameter in the epoch combination frame; They are t The polar motion rates of the two components in the epoch combination frame; for t Variations in the Earth's rotation rate in the epoch-combined frame; Both are rotation matrices The elements of correspond to the angles about the X and Y axes, respectively.

[0023] Preferably, during the time update process, the state equation used is as follows:

[0024]

[0025] In the formula, They are Epoch, the reference epoch coordinates of the epoch; , They are Epoch, The station linear velocity at the epoch; They are Epoch, Time-dependent noise parameters of the epoch; , and They are Epoch, Epoch and The annual term of the seasonal signal of the epoch; , and They are Epoch, Epoch and Semi-annual term of seasonal signal of epoch; is the state conversion coefficient of the time-dependent noise parameter; is the standard deviation of the coordinate parameter process noise, is the standard deviation of the linear velocity parameter process noise, for The driving white noise of the time-dependent noise parameter of the epoch, is the time interval between two consecutive epochs; and are the periods of the anniversary and semi-anniversary items respectively; τ is a set constant; e is a natural constant; π is the circular velocity constant.

[0026] Preferably, in the process of forecasting the estimated parameter, the current epoch time t Add 1 and judge Whether the coordinates or linear velocity of the frame point of the epoch have changed;

[0027] If the coordinates change, the standard deviation of the process noise of the epoch coordinate parameter is Adjust to the first optimal empirical value; if the linear velocity jumps, the standard deviation of the process noise of the linear velocity parameter of the epoch Adjust to the second optimal empirical value; after adjusting the above parameters, perform Forecast of the parameter to be estimated at the epoch;

[0028] If neither the coordinates nor the linear velocity jumps, the state equation and the noise analysis result information are directly used to perform Forecast of the parameter to be estimated at the epoch.

[0029] Preferably, when performing square root information smoothing, reverse epoch-by-epoch parameter smoothing is performed, and the optimal estimate of the parameter to be estimated in the previous epoch is calculated by the smoothing formula;

[0030] Determine whether the current epoch is the first epoch. If it is not the first epoch, decrement the epoch time by one and then recursively calculate again according to the smoothing formula until the current epoch is the first epoch;

[0031] The smoothing formula is expressed as:

[0032]

[0033] In the formula, for t The best estimate of the parameter to be estimated at the epoch, After time update t The autocorrelation matrix of the parameters to be estimated at the epoch, express The inverse matrix of After time update t The cross-correlation matrix of the parameters to be estimated between the epoch and the t+1 epoch; for The best estimate of the parameter to be estimated at the epoch, for The right side vector after the epoch time is updated;

[0034] The solution result information includes the parameter solution result of the last epoch, corresponding to the initial calculation in the smoothing calculation. ; The temporary files include , and .

[0035] In another aspect, the present invention provides a filtering system for establishing a nonlinear earth reference frame, comprising:

[0036] Noise analysis module, used to perform noise analysis on the frame point coordinate time series to obtain noise analysis result information;

[0037] A preprocessing module is used to preprocess the frame point coordinates to obtain unconstrained frame point coordinate observation values;

[0038] A parameter solution module, used for performing parameter solution based on square root information filtering according to the noise analysis result information and the unconstrained frame point coordinate observation values, obtaining parameter solution result information, and storing a correlation matrix as a temporary file;

[0039] A smoothing processing and earth reference frame establishment module is used to perform square root information smoothing processing according to the parameter solution result information and the temporary file, obtain the optimal estimate of the parameter to be estimated at each epoch, and complete the establishment of a nonlinear earth reference frame;

[0040] The filtering system for establishing a nonlinear earth reference frame is used to execute the steps in the above-mentioned filtering method for establishing a nonlinear earth reference frame.

[0041] One or more technical solutions provided in the present invention have at least the following technical effects or advantages:

[0042] The present invention aims at the problem that the nonlinear changes of frame points cannot be accurately described when the earth reference frame is established based on the filtering method. The first-order autoregressive model (i.e., AR[1]) is used in the filtering algorithm to approximate the time-correlated noise of the GNSS station coordinates. That is, the time-correlated noise of the frame points is modeled by introducing the AR[1] model in the state equation, which can more accurately describe the nonlinear motion of the measuring station and ensure the accuracy of the reference frame establishment results. In addition, the present invention uses square root information filtering (SRIF) to establish the earth reference frame. Compared with the traditional Kalman filter, the square root information filter has the characteristics of higher numerical accuracy and strong stability, and can well overcome the problem of least squares time series analysis method ignoring the random model and the filter divergence problem that may occur in the Kalman filter time series analysis. The present invention provides important theoretical and technical support for the further development of establishing the earth reference frame based on filtering technology and the realization of the millimeter-level reference frame. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 A flowchart of a filtering method for establishing a nonlinear earth reference frame provided in Example 1 of the present invention. DETAILED DESCRIPTION

[0044] In order to better understand the above technical solution, the above technical solution will be described in detail below in conjunction with the accompanying drawings and specific implementation methods.

[0045] Embodiment 1:

[0046] Embodiment 1 provides a filtering method for establishing a nonlinear earth reference frame, comprising the following steps:

[0047] Perform noise analysis on the frame point coordinate time series to obtain noise analysis result information; in the noise analysis, a first-order autoregressive model is introduced to model the frame point time-related noise;

[0048] Preprocess the frame point coordinates to obtain unconstrained frame point coordinate observation values;

[0049] According to the noise analysis result information and the unconstrained frame point coordinate observation values, parameter solution is performed based on square root information filtering to obtain parameter solution result information, and a correlation matrix is ​​stored as a temporary file;

[0050] According to the parameter solution result information and the temporary file, square root information smoothing processing is performed to obtain the optimal estimate of the parameter to be estimated at each epoch, and the establishment of the nonlinear earth reference frame is completed.

[0051] The noise analysis of the framework point coordinate time series includes: taking the framework point coordinates in the local coordinate system of the measuring station as input data, adopting a noise model combination of variable white noise and a first-order autoregressive model, and performing noise analysis using a maximum likelihood estimation method to obtain a noise analysis result, wherein the noise analysis result includes the state transfer coefficient of the first-order autoregressive process of all measuring stations, the standard deviation of the driving white noise, and the coefficient factor of the variable white noise; converting the noise analysis result from the local coordinates of the measuring station to the geocentric Cartesian coordinates to obtain the noise analysis result information.

[0052] Before performing the noise analysis, the method may further include: performing a process of eliminating gross errors on the input data by least square fitting.

[0053] The preprocessing of the frame point coordinates includes: eliminating the prior constraints attached to the frame point coordinates according to the constraint matrix.

[0054] Wherein, the parameter solution based on square root information filtering includes: measurement update and time update. The measurement update includes: establishing observation equations, parameterizing the unconstrained frame point coordinate observations by epoch-by-epoch station coordinate decomposition, performing measurement update parameter solution, and obtaining the information matrix after measurement update; and using the current measurement updated information matrix as the prior information for the next epoch time update. The time update includes: introducing a first-order autoregressive model to model the frame point time-related noise, using the state equation to predict the parameters to be estimated, and obtaining the information matrix after time update; and using the current time updated information matrix as the prior information for the next epoch measurement update. Wherein, the parameters to be estimated include the reference epoch coordinates of each epoch, the linear velocity of the frame point, the time-related noise parameters, the seasonal signal anniversary term, and the seasonal signal semi-anniversary term.

[0055] Specifically, when establishing the observation equation, the unconstrained frame point coordinate observation value is parameterized as follows:

[0056]

[0057] In the formula, t epoch as the current epoch, for t The unconstrained frame point coordinate observations for the epoch, for t the reference epoch coordinates of the epoch, for t Epoch and reference epoch The time interval for t The linear velocity of the epoch frame point, for t The time-dependent noise parameters of the epoch, for t The annual term of the seasonal signal at the epoch, for t The semi-annual term of the seasonal signal of the epoch, They are t The translation, scaling and rotation parameters between the epoch combined reference frame and the known reference frame, for t The coordinates of the frame points of the epoch combined reference frame, is the coefficient factor of variable white noise, for t Variable white noise term for the epochs.

[0058] During the measurement update process, the observation equation includes the earth orientation parameter observation equation, which is expressed as follows:

[0059]

[0060] In the formula, They are t Observed values ​​of the two components of the epoch polar motion parameter; They are t Observed values ​​of the polar motion rate of the two components of the epoch; for t Observations of changes in the Earth's rotation rate over the epoch; They are t Two components of the polar motion parameter in the epoch combination frame; They are t The polar motion rates of the two components in the epoch combination frame; for t Variations in the Earth's rotation rate in the epoch-combined frame; Both are rotation matrices The elements of correspond to the angles about the X and Y axes, respectively.

[0061] During the time update process, the state equation used is as follows:

[0062]

[0063] In the formula, They are Epoch, the reference epoch coordinates of the epoch; , They are Epoch, The station linear velocity at the epoch; They are Epoch, Time-dependent noise parameters of the epoch; , and They are Epoch, Epoch and The annual term of the seasonal signal of the epoch; , and They are Epoch, Epoch and Semi-annual term of seasonal signal of epoch; is the state conversion coefficient of the time-dependent noise parameter; is the standard deviation of the coordinate parameter process noise, is the standard deviation of the linear velocity parameter process noise, for The driving white noise of the time-dependent noise parameter of the epoch, is the time interval between two consecutive epochs; and are the periods of the anniversary and semi-anniversary items respectively; τ is a set constant; e is a natural constant; π is the circular velocity constant.

[0064] In the process of forecasting the estimated parameters, the current epoch time t Add 1 and judge Whether the coordinates or linear velocity of the frame point of the epoch have changed; if the coordinates have changed, the standard deviation of the process noise of the epoch coordinate parameters is Adjust to the first optimal empirical value; if the linear velocity jumps, the standard deviation of the process noise of the linear velocity parameter of the epoch Adjust to the second optimal empirical value; after adjusting the above parameters, perform If neither the coordinates nor the linear velocity change, the prediction of the parameters to be estimated at the epoch is performed directly according to the state equation and the noise analysis result information. Forecast of the parameter to be estimated at the epoch.

[0065] Specifically, when performing square root information smoothing, reverse epoch-by-epoch parameter smoothing is performed, and the optimal estimate of the parameter to be estimated in the previous epoch is calculated by the smoothing formula; it is determined whether the current epoch is the first epoch. If it is not the first epoch, the epoch time is reduced by one and then recursively calculated, and the calculation is performed again according to the smoothing formula until the current epoch is the first epoch.

[0066] The smoothing formula is expressed as:

[0067]

[0068] In the formula, for t The best estimate of the parameter to be estimated at the epoch, After time update t The autocorrelation matrix of the parameters to be estimated at the epoch, express The inverse matrix of After time update t The cross-correlation matrix of the parameters to be estimated between the epoch and the t+1 epoch; for The best estimate of the parameter to be estimated at the epoch, for The right side vector after the epoch time update.

[0069] The solution result information includes the parameter solution result of the last epoch, corresponding to the initial calculation in the smoothing calculation. ; The temporary files include , and .

[0070] Example 1 introduces a first-order autoregressive (AR[1]) model instead of a random walk to approximate the time-correlated noise in the frame point coordinate sequence, and models and modifies it in the random model according to the maximum likelihood estimation noise level, and further optimizes the filter state equation and its random model to improve the reference frame accuracy. On the one hand, it helps to further refine the velocity model of the existing GNSS station coordinate time series and improve the reliability of the linear velocity and its uncertainty estimation; on the other hand, compared with the traditional random walk filter state equation, the time-correlated noise model established by the AR[1] process significantly reduces the frame point coordinate residual.

[0071] The present invention is further described below. Figure 1 The present invention provides a filtering method for establishing a nonlinear earth reference frame, which specifically comprises the following steps:

[0072] Step 1: Noise analysis of frame point coordinate time series.

[0073] It includes taking the frame point coordinates in the local coordinate system (NEU) of the measuring station as input data to perform station coordinate sequence noise analysis, and converting the noise analysis results from NEU to geocentric Cartesian coordinates (XYZ).

[0074] The following is a detailed description of step 1.

[0075] The accuracy of the state equation and its random model is one of the main factors affecting the establishment of a reference frame based on filtering technology. In order to further optimize the filtering state equation and its random model, the present invention uses a first-order autoregression (AR[1]) instead of a random walk model to model the flicker noise in the GNSS station coordinate time series. The purpose of the station coordinate time series noise analysis is to determine the state transfer coefficient and process noise of the AR[1] model, and then optimize the filtering state equation.

[0076] First, extract the station coordinates and their covariance from the SINEX product (ftp: / / igs-rf.ign.fr / pub / repro3 / ) provided by the International GNSS service (IGS), and convert them from geocentric Cartesian coordinates (XYZ) to station local coordinates (NEU):

[0077] (1)

[0078] Where N, E and U are the coordinates of the frame point in the north, east and height components in the local coordinate system of the station, B and L are the geodetic latitude and longitude of the frame point, X, Y and Z are the coordinates of the frame point in the geocentric Cartesian coordinate system, , and The geocentric Cartesian coordinates specified when performing coordinate system conversion.

[0079] Then, the GNSS station coordinates in the local coordinate system of the station are used as input data, and the input data are processed by least squares fitting to eliminate gross errors, such as eliminating observations with residuals greater than 3 times the mean error. By specifying a noise model combination of variable white noise and first-order autoregression, the maximum likelihood estimation method is used to perform station coordinate time series noise analysis, thereby obtaining the state transfer coefficients and standard deviations of the driving white noise of the AR[1] process of all stations, as well as the coefficient factors of the variable white noise. Finally, the noise analysis results are converted from the local coordinates of the station to the geocentric Cartesian coordinates.

[0080] Step 2: Preprocessing of frame point coordinates.

[0081] According to the constraint matrix provided by the SINEX file, the prior constraints attached to the station coordinates are eliminated, and the unconstrained station coordinates in the XYZ coordinate system (i.e., unconstrained GNSS station coordinate observations) are obtained to establish the earth reference frame.

[0082] The following is a detailed description of step 2.

[0083] The purpose of frame point coordinate preprocessing is to eliminate the prior constraints attached to the station coordinates and obtain unconstrained XYZ station coordinates for establishing the earth reference frame.

[0084] A single geodetic technique cannot provide enough information to estimate the station coordinates and earth orientation parameters. Therefore, the least squares adjustment method is used to process the rank deficiency of the normal equation matrix when processing the GNSS raw data. In order to solve the rank deficiency of the normal equation matrix, IGS adds additional minimum constraints to the normal equation for parameter estimation when processing the raw data:

[0085] (2)

[0086] In the formula, is the normal equation matrix, is the additional constraint matrix, is the estimated value of the parameter to be estimated, is the prior value of the parameter to be estimated, is the right side vector of the normal equation, Estimate the parameters The covariance matrix of .

[0087] The SINEX format product provided by IGS contains all the parameters in equation (2), which can be used to recover the unconstrained normal equation and obtain the unconstrained frame point coordinate observation values:

[0088] (3)

[0089] In the formula, Estimate the parameters The inverse covariance matrix.

[0090] Step 3: Parameter solution based on square root information filter (SRIF).

[0091] Parameter solution using square root information filtering includes two parts: measurement update and time update.

[0092] The measurement update includes taking the frame point coordinates and covariance under XYZ as input data, parameterizing the frame point station coordinates of different components by epoch-by-epoch station coordinate decomposition, solving the measurement update parameters, and obtaining the information matrix after the measurement update; the information matrix after each parameter solution is completed (i.e., the information matrix after the current measurement update) is used as the prior information for the next epoch time update.

[0093] The time update includes introducing a first-order autoregressive model to model the time-related noise of the framework point, using the state equation to predict the evolution of the estimated parameters over time, and storing the matrix containing the parameter evolution (i.e., the information matrix after the time update) on a local computer for subsequent reverse smoothing processing. The information matrix after the current time update will be used as the prior information for the next epoch measurement update.

[0094] The present invention introduces a first-order autoregressive model (AR[1]) to model time-correlated noise. Compared with the traditional random walk model, it can effectively improve the modeling accuracy of time-correlated noise parameters.

[0095] The following is a detailed description of step 3.

[0096] The purpose of the measurement update is to establish the observation equation, parameterize the frame coordinate observations and solve the parameters to be estimated. The filter is initialized based on the coordinate observations of the first epoch station, and the frame point coordinates in the geocentric Cartesian coordinates are decomposed into:

[0097] (4)

[0098] In the formula, t epoch as the current epoch, t The unconstrained frame point coordinate observations for the epoch, for t the reference epoch coordinates of the epoch, fort The linear velocity of the epoch frame point, is the epoch t and the reference epoch The time interval for t The time-dependent noise parameters of the epoch, for t The annual term of the seasonal signal at the epoch, for t The semi-annual term of the seasonal signal of the epoch, They are t The translation, scaling and rotation parameters between the epoch combined reference frame and the known reference frame, for t The coordinates of the frame points of the epoch combined reference frame, for t The variable white noise term at the epoch, is a variable white noise coefficient factor. In the preferred embodiment, is the median of the variable white noise coefficient factor ( determined by step 1).

[0099] In addition, the Earth Orientation Parameter (EOP) observation equation is added to the measurement update observation equation:

[0100] (5)

[0101] In the formula, They are t Observed values ​​of the two components of the epoch polar motion parameter; They are t Observed values ​​of the polar motion rate of the two components of the epoch; for t Observations of changes in the Earth's rotation rate over the epoch; They are t Two components of the polar motion parameter in the epoch combination frame; They are t The polar motion rates of the two components in the epoch combination frame; for t Variations in the Earth's rotation rate in the epoch-combined frame; Both are rotation matrices The elements of correspond to the angles about the X and Y axes, respectively.

[0102] The initial information and observation equation of the parameter vector to be estimated can be expressed as:

[0103] (6)

[0104] In the formula, is the estimated value of the parameter to be estimated, is the prior value of the parameter to be estimated, is the correction number of the prior value of the parameter to be estimated, is the difference between the station coordinate observation and the Taylor expansion term, is the coefficient matrix of the observation equation, is the correction number of the station coordinate observation value.

[0105] and are the covariance matrices of the prior values ​​of the parameters to be estimated and the frame point coordinate observations. It should be noted that the variance of the frame point coordinate observations in the input data SINEX file is usually small, so it is necessary to Matrix multiplication by coefficient factors . and Square root decomposition of a matrix:

[0106] (7)

[0107] In the formula, for The matrix after square root decomposition, for The inverse matrix of for The transposed matrix of for The matrix after square root decomposition, for The inverse matrix of for The transposed matrix of .

[0108] Multiply the above and the following equations by and Matrix, and let the matrix , , , , , then equation (6) can be expressed as:

[0109] (8)

[0110] Perform an orthogonal transformation on equation (8), that is, multiply both sides of the equation by the orthogonal matrix :

[0111] (9)

[0112] In the formula, is the information matrix after orthogonal transformation, is the right side vector after orthogonal transformation, is the parameter vector correction after orthogonal transformation, is the normalized observed posterior residual vector, is the residual vector correction number.

[0113] Define the objective function :

[0114] (10)

[0115] In the formula, is the matrix bi-norm operator.

[0116] To make Minimum, then the parameter solution in equation (10) needs to satisfy:

[0117] (11)

[0118] The above process completes the measurement update parameter solution:

[0119] (12)

[0120] Measure the updated information matrix It can be used as the prior information of the parameter vector when updating the next epoch time, where is the normalized observed posterior residual vector.

[0121] After completing the parameter solution through measurement update, determine the current t Is the epoch the last epoch? If t The epoch is not the last epoch, then the parameters to be estimated are predicted by performing a time update; if t If the epoch is the last epoch, the process proceeds to step 4 to perform square root information smoothing.

[0122] The purpose of time updating is to predict the parameters to be estimated in the next epoch based on geophysical principles. First, add 1 to the current epoch time t and determine Whether the frame point station coordinates or linear velocity of the epoch has changed. If the station coordinates have changed, the standard deviation of the process noise of the epoch reference epoch coordinate parameters is used. Adjust to the optimal empirical value of multiple test results (recorded as the first optimal empirical value), for example, the first optimal empirical value is 1m; if the linear velocity jumps, the standard deviation of the process noise of the linear velocity parameter of the epoch is Adjust to the optimal empirical value of multiple test results (recorded as the second optimal empirical value), for example, the second optimal empirical value is 0.3 m / yr; if the coordinate observation value and linear velocity of the epoch station do not jump, then use the state equation and the process noise determined in step 1 to calculate The prediction of the parameter to be estimated at the epoch. Subsequently, the information matrix after the prediction of the parameter to be estimated (i.e., the information matrix after the current time update) is used as the prior information for the next epoch measurement update, and the correlation matrix containing the evolution of the parameter to be estimated over time (i.e., the information matrix after all time updates) is stored locally for the square root information smoothing processing in step 4. In the above process, the state equation for the prediction of the parameter to be estimated is:

[0123] (13)

[0124] In the formula, They are Epoch, the reference epoch coordinates of the epoch; , They are Epoch, The linear velocity of the station (i.e., GNSS monitoring station) at the epoch; They are Time-dependent noise parameters of the epoch; , and They are Epoch, Epoch and The annual term of the seasonal signal of the epoch; , and They are Epoch, Epoch and Semi-annual term of seasonal signal of epoch; is the state transfer coefficient of the time-dependent noise parameter; is the standard deviation of the coordinate parameter process noise, is the standard deviation of the linear velocity parameter process noise, and is an infinitesimal quantity, for The driving white noise of the time-dependent noise parameters of the epochs, whose standard deviation is determined in the noise analysis of the frame point station coordinate time series in step 1. is the time interval between two consecutive epochs; and are the periods of the anniversary and semi-anniversary items, i.e., 1 year and 0.5 years respectively; τ is a set constant, usually a large constant; e It is a natural constant in mathematical symbols; π is the circumference constant.

[0125] In the formula, and They are and The parameter vector to be estimated for the epoch, is the state transfer matrix, is the process noise vector, and its covariance matrix ,in for The matrix after square root decomposition, for The inverse matrix of for The transposed matrix of .

[0126] in, From the AR[1] model = + To improve the stability of the Earth reference frame, Preferably the median of the frame points is used.

[0127] Multiply both ends of equation (13) by the square root decomposition Matrix, and then perform orthogonal transformation:

[0128] (14)

[0129] In the formula, for The information matrix after epoch measurement update, for The right-hand side vector after the epoch measurement update, for The matrix after square root decomposition, is the autocorrelation matrix of the estimated parameters after time update, is the cross-correlation matrix of the estimated parameters after time update, for The right-hand side vector after the epoch time is updated, for The information matrix after the epoch time is updated, for The updated vector of epoch times.

[0130] The right-hand side of equation (14) gives the smoothing algorithm for the epochs before and after the estimated parameter. , and The square root information is stored locally in the computer and used for smoothing the square root information in step 4. The following formula on the right side of equation (14) is the prediction of the estimated parameter updated at the SRIF time, which is passed as the prior information of the estimated parameter to the measurement update of the next epoch.

[0131] Step 4: Square root information smoothing (SRIS) processing.

[0132] It includes reading the matrix stored on the local computer, obtaining and outputting the optimal parameter estimates for each epoch, and then establishing the earth reference frame.

[0133] The following is a detailed description of step 4.

[0134] According to the basic principle of SRIF filter, in step 3 t The optimal estimate of the parameter to be estimated in the epoch is determined only by historical observations (including epoch t). In order to obtain the optimal estimate of the unknown parameter using all epoch observations, reverse epoch-by-epoch parameter smoothing is performed:

[0135] (15)

[0136] In the formula, , and is the temporary matrix stored locally in this calculation in step 3, is the autocorrelation matrix of the estimated parameter at epoch t after time update, express The inverse matrix of is the cross-correlation matrix of the estimated parameters between epoch t and epoch t+1 after time update, for The right-hand side vector after the epoch time is updated; is the optimal estimate of the parameter to be estimated at epoch t+1 in step 3, is the optimal estimate of the parameter to be estimated at epoch t.

[0137] When processing square root information smoothing (SRIS), first, the result is calculated based on the parameters of the last epoch in step 3 , and the temporary matrix stored for the epoch , , , calculate the optimal estimate of the parameter value to be estimated in the previous epoch by formula (15) , and output the optimal estimate; then, determine whether the current epoch is the first epoch. If it is not the first epoch, decrement the epoch time by one and then recursively calculate the parameters again according to formula (15) until the current epoch is the first epoch, and then end the program operation to complete the establishment of the earth reference frame.

[0138] Embodiment 2:

[0139] Embodiment 2 provides a system for establishing a nonlinear earth reference frame, comprising:

[0140] Noise analysis module, used to perform noise analysis on the frame point coordinate time series to obtain noise analysis result information;

[0141] A preprocessing module is used to preprocess the frame point coordinates to obtain unconstrained frame point coordinate observation values;

[0142] A parameter solution module, used for performing parameter solution based on square root information filtering according to the noise analysis result information and the unconstrained frame point coordinate observation values, obtaining parameter solution result information, and storing a correlation matrix as a temporary file;

[0143] A smoothing processing and earth reference frame establishment module is used to perform square root information smoothing processing according to the parameter solution result information and the temporary file, obtain the optimal estimate of the parameter to be estimated at each epoch, and complete the establishment of a nonlinear earth reference frame;

[0144] The filtering system for establishing a nonlinear earth reference frame provided in Example 2 is used to execute the steps in the filtering method for establishing a nonlinear earth reference frame as described in Example 1.

[0145] Since the functions of each module in the system provided in Example 2 correspond to the steps in the method provided in Example 1, Example 2 can be understood by referring to the description of Example 1, which will not be repeated here.

[0146] In addition, during specific implementation, the method proposed in the technical solution of the present invention can be automatically run by those skilled in the art using computer software technology. System devices for implementing the method, such as computer-readable storage media storing the corresponding computer program of the technical solution of the present invention and computer equipment running the corresponding computer program, should also be within the protection scope of the present invention.

[0147] In some possible embodiments, the present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps in a filtering method for establishing a nonlinear earth reference frame as described in Example 1.

[0148] In some possible embodiments, the present invention also provides a computer-readable storage medium having a computer program / instruction stored thereon, which, when executed by a processor, implements the steps in a filtering method for establishing a nonlinear earth reference frame as described in Example 1.

[0149] In some possible embodiments, the present invention also provides a computer program product, including a computer program / instruction, which, when executed by a processor, implements the steps in a filtering method for establishing a nonlinear earth reference frame as described in Example 1.

[0150] Finally, it should be noted that the above specific implementation methods are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to examples, those skilled in the art should understand that the technical solution of the present invention can be modified or replaced by equivalents without departing from the spirit and scope of the technical solution of the present invention, which should be included in the scope of the claims of the present invention.

Claims

1. A filtering method for establishing a nonlinear earth reference frame, characterized in that: The following steps are involved: Perform noise analysis on the frame point coordinate time series to obtain noise analysis result information; in the noise analysis, a first-order autoregressive model is introduced to model the frame point time-related noise; Preprocess the frame point coordinates to obtain unconstrained frame point coordinate observation values; According to the noise analysis result information and the unconstrained frame point coordinate observation values, parameter solution is performed based on square root information filtering to obtain parameter solution result information, and a correlation matrix is ​​stored as a temporary file; According to the parameter solution result information and the temporary file, square root information smoothing processing is performed to obtain the optimal estimate of the parameter to be estimated at each epoch, and the establishment of a nonlinear earth reference frame is completed; The parameter solution based on square root information filtering includes: measurement update and time update; The measurement update includes: establishing an observation equation, parameterizing the unconstrained frame point coordinate observation value by epoch-by-epoch station coordinate decomposition, solving the measurement update parameters, and obtaining the information matrix after the measurement update; and using the current measurement updated information matrix as the prior information for the next epoch time update; The time update includes: introducing a first-order autoregressive model to model the time-related noise of the frame point, using the state equation to predict the parameters to be estimated, and obtaining the information matrix after the time update; and using the information matrix after the current time update as the prior information for the next epoch measurement update; The parameters to be estimated include the reference epoch coordinates of each epoch, the linear velocity of the frame point, the time-correlated noise parameter, the seasonal signal anniversary term and the seasonal signal semi-anniversary term; When establishing the observation equation, the unconstrained frame point coordinate observation value is parameterized as follows: In the formula, t epoch as the current epoch, for t The unconstrained frame point coordinate observations for the epoch, for t the reference epoch coordinates of the epoch, for t Epoch and reference epoch The time interval for t The linear velocity of the epoch frame point, for t The time-dependent noise parameters of the epoch, for t The annual term of the seasonal signal at the epoch, for t The semi-annual term of the seasonal signal of the epoch, They are t The translation, scaling and rotation parameters between the epoch combined reference frame and the known reference frame, for t The coordinates of the frame points of the epoch combined reference frame, is the coefficient factor of variable white noise, for t variable white noise term for the epoch; During the measurement update process, the observation equation includes the earth orientation parameter observation equation, which is expressed as follows: In the formula, They are t Observed values ​​of the two components of the epoch polar motion parameter; They are t Observed values ​​of the polar motion rate of the two components of the epoch; for t Observations of changes in the Earth's rotation rate over the epoch; They are t Two components of the polar motion parameter in the epoch combination frame; They are t The polar motion rates of the two components in the epoch combination frame; for t Variations in the Earth's rotation rate in the epoch-combined frame; Both are rotation matrices The elements of correspond to the angles about the X and Y axes respectively; During the time update process, the state equation used is as follows: In the formula, They are Epoch, the reference epoch coordinates of the epoch; , They are Epoch, The station linear velocity at the epoch; They are Epoch, Time-dependent noise parameters of the epoch; , and They are Epoch, Epoch and The annual term of the seasonal signal of the epoch; , and They are Epoch, Epoch and Semi-annual term of seasonal signal of epoch; is the state conversion coefficient of the time-dependent noise parameter; is the standard deviation of the coordinate parameter process noise, is the standard deviation of the linear velocity parameter process noise, for The driving white noise of the time-dependent noise parameter of the epoch, is the time interval between two consecutive epochs; and are the periods of the anniversary and semi-anniversary items respectively; τ is a set constant; e is a natural constant; π is the circular velocity constant.

2. The filtering method for establishing a nonlinear earth reference frame according to claim 1, characterized in that: The noise analysis of the framework point coordinate time series includes: taking the framework point coordinates in the local coordinate system of the measuring station as input data, adopting a noise model combination of variable white noise and a first-order autoregressive model, and performing noise analysis using a maximum likelihood estimation method to obtain a noise analysis result, wherein the noise analysis result includes a state transfer coefficient of the first-order autoregressive process of all measuring stations, a standard deviation of the driving white noise, and a coefficient factor of the variable white noise; and converting the noise analysis result from the local coordinates of the measuring station to the geocentric Cartesian coordinates to obtain the noise analysis result information.

3. The filtering method for establishing a nonlinear earth reference frame according to claim 2, characterized in that: Before performing noise analysis, it also includes: eliminating gross errors from input data through least square fitting.

4. The filtering method for establishing a nonlinear earth reference frame according to claim 1, characterized in that: The preprocessing of the frame point coordinates includes: eliminating the prior constraints attached to the frame point coordinates according to the constraint matrix.

5. The filtering method for establishing a nonlinear earth reference frame according to claim 1, characterized in that: In the process of forecasting the estimated parameters, the current epoch time t Add 1 and judge Whether the coordinates or linear velocity of the frame point of the epoch have changed; If the coordinates change, the standard deviation of the process noise of the epoch coordinate parameter is Adjust to the first optimal empirical value; if the linear velocity jumps, the standard deviation of the process noise of the linear velocity parameter of the epoch Adjust to the second optimal empirical value; after adjusting the above parameters, perform Forecast of the parameters to be estimated at the epoch; If neither the coordinates nor the linear velocity jumps, the state equation and the noise analysis result information are directly used to perform Forecast of the parameter to be estimated at the epoch.

6. The filtering method for establishing a nonlinear earth reference frame according to claim 5, characterized in that: When performing square root information smoothing, reverse epoch-by-epoch parameter smoothing is performed, and the optimal estimate of the parameter to be estimated in the previous epoch is calculated by the smoothing formula; Determine whether the current epoch is the first epoch. If it is not the first epoch, decrement the epoch time by one and then recursively calculate again according to the smoothing formula until the current epoch is the first epoch; The smoothing formula is expressed as: In the formula, for t The best estimate of the parameter to be estimated at the epoch, After time update t The autocorrelation matrix of the parameters to be estimated at the epoch, express The inverse matrix of After time update t The cross-correlation matrix of the parameters to be estimated between the epoch and the t+1 epoch; for The best estimate of the parameter to be estimated at the epoch, for The right side vector after the epoch time is updated; The solution result information includes the parameter solution result of the last epoch, corresponding to the initial calculation in the smoothing calculation. ; The temporary files include , and .

7. A filtering system for establishing a nonlinear earth reference frame, characterized in that: include: Noise analysis module, used to perform noise analysis on the frame point coordinate time series to obtain noise analysis result information; A preprocessing module is used to preprocess the frame point coordinates to obtain unconstrained frame point coordinate observation values; A parameter solution module, used for performing parameter solution based on square root information filtering according to the noise analysis result information and the unconstrained frame point coordinate observation values, obtaining parameter solution result information, and storing a correlation matrix as a temporary file; A smoothing processing and earth reference frame establishment module is used to perform square root information smoothing processing according to the parameter solution result information and the temporary file, obtain the optimal estimate of the parameter to be estimated at each epoch, and complete the establishment of a nonlinear earth reference frame; The filtering system for establishing a nonlinear earth reference frame is used to execute the steps in the filtering method for establishing a nonlinear earth reference frame as described in any one of claims 1 to 6.

Images