A denoising method for coordinate time series

Abstract

The application provides a coordinate time sequence denoising method, which comprises the following steps: firstly, obtaining an original coordinate time sequence and pre-processing the same; secondly, converting the pre-processed coordinate time sequence into a trajectory matrix, and then performing singular value decomposition on the coordinate time sequence based on the trajectory matrix to obtain decomposition components and calculate the MIC value between each decomposition component and time; finally, selecting the first P components with larger MIC values to reconstruct, taking the reconstructed coordinate time sequence as a denoising result, and thus realizing denoising of the original coordinate time sequence. The application determines the decomposition components for reconstruction through the MIC adaptive criterion, comprehensively considers the MIC values of the decomposition components, the reconstructed sequence and the residual sequence, guarantees that the reconstructed coordinate time sequence has a strong nonlinear relationship, and takes into account more decomposition components and less signal residues, thus realizing scientific and effective denoising.

Description

TECHNICAL FIELD

[0001] The present application relates to a coordinate time series denoising method, belonging to the coordinate time series denoising technical field. BACKGROUND

[0002] The coordinate time series contains a lot of important information related to plate tectonic movement, surface mass migration and other geodynamic processes. These information components are often time-dependent, which can provide rich data basis for geodynamic research. The time-dependent components in the coordinate time series mainly consist of the following parts: long-term changes in station coordinates caused by tectonic movement, seasonal changes in station coordinates caused by geophysical effects, post-seismic deformation caused by earthquakes and other time-dependent coordinate changes. However, these components that can reflect geophysical processes are often mixed with noise or even covered by noise. The noise in the coordinate time series mainly comes from observation errors and errors caused by imperfect data processing strategies. Therefore, effective denoising of the coordinate time series can achieve extraction of time-dependent components, thereby providing accurate data for geodynamic research. At the same time, the coordinate time series is also a key basic data for the establishment and maintenance of the Earth Reference Frame, and a scientific denoising method is helpful for accurate estimation of benchmark station coordinates and velocity and high-precision maintenance of the Earth Reference Frame.

[0003] There are many existing coordinate time series denoising methods. Some of them use local mean decomposition to decompose the coordinate time series data and remove the identified noise components, but this will result in the loss of some useful information. Some others use singular spectrum analysis to decompose and reconstruct the time series. As a data-driven non-parametric method, singular spectrum analysis can extract useful information from time series containing noise through decomposition and reconstruction of the time series without any prior information, thereby achieving denoising. However, the determination of the reconstructed components in singular spectrum analysis has a great influence on the denoising effect. The reconstructed components are currently determined according to the contribution rate, and the determination of the reconstructed components is often based on subjective judgment, lacking scientific adaptive criteria, which may result in loss of useful information and cannot guarantee the denoising effect. SUMMARY

[0004] The purpose of the present application is to provide a coordinate time series denoising method to solve the problem of loss of useful information due to lack of scientific adaptive criteria in the prior art coordinate time series data denoising process.

[0005] The present application provides a coordinate time series denoising method, which comprises the following steps:

[0006] 1) obtaining an original coordinate time series;

[0007] 2) convert the obtained coordinate time series into a corresponding trajectory matrix, perform singular value decomposition on the coordinate time series based on the converted trajectory matrix to obtain decomposition components;

[0008] 3) calculate the MIC value between each decomposition component and time to obtain the MIC value of each decomposition component, sort the decomposition components according to the MIC value of each decomposition component, select P decomposition components with larger MIC values to reconstruct, and generate a reconstructed sequence, which is a coordinate time series after noise reduction, to realize noise reduction on the original coordinate time series;

[0009] The selection of P is based on: the MIC value of the reconstructed sequence generated by the first P decomposition components is greater than a first set threshold, and the MIC value of the reconstructed sequence generated by the first P+1 decomposition components is less than the first set threshold, or the MIC value of the reconstructed sequence generated by the first P decomposition components is greater than the first set threshold and the MIC value of the residual sequence is less than a second set threshold, wherein the residual sequence is a sequence composed of the decomposition components remaining after selecting P decomposition components.

[0010] The present application realizes the noise reduction processing of the coordinate time series by using singular value decomposition and decomposition component reconstruction. When reconstructing the decomposition components, the MIC criterion is introduced, and the MIC values of the decomposition components, the reconstructed sequence and the residual sequence are considered comprehensively to ensure that the reconstructed coordinate time series has strong nonlinear relationship, while taking into account more decomposition components and less signal residue, and realizing scientific and effective noise reduction.

[0011] Further, the first set threshold in step 3) has a value range of [0.8, 1]; the second set threshold has a value range of [0, 0.3].

[0012] By the above method, the strong time correlation of the reconstructed sequence and the weak time correlation of the residual sequence can be ensured.

[0013] Further, the original coordinate time series data needs to be preprocessed before singular value decomposition in step 2), and the preprocessing includes detection and elimination of gross errors and jumps and interpolation of missing values.

[0014] By detecting and eliminating the gross errors and jumps of the original coordinate time series, a part of noise caused by measurement error can be reduced, the data quality can be improved, and the accuracy of the subsequent reconstructed sequence can be ensured; by interpolating the missing values, the uniformity of sampling can be ensured.

[0015] Further, the preprocessing further includes zero mean processing of the original coordinate time series.

[0016] In order to meet the requirement of data zero mean in the subsequent calculation process, the zero mean processing is carried out on the original coordinate time sequence, so that the processed data fluctuates up and down at the 0 value and is uniformly distributed, the linear trend item is eliminated, the coordinate time sequence is analyzed conveniently, and the subsequent operation is carried out smoothly.

[0017] Further, the trajectory matrix is an M-row and L-column matrix, wherein 1

[0018] Further, the calculation formula of the decomposition component is:

[0019]

[0020]

[0021] wherein, is the i th value of the k th reconstruction component, x i is the i th element in the preprocessed coordinate time sequence, D i is the trajectory matrix, E k is the eigenvector corresponding to the k th eigenvalue in the covariance matrix of the trajectory matrix, a k is the coefficient of the corresponding eigenvector E k is the j th element in E k is the i th element of a k . BRIEF DESCRIPTION OF DRAWINGS

[0022] Figure 1 is the flow chart of the coordinate time sequence denoising method of the application;

[0023] Fig. 2 (a) is a reconstructed sequence result graph after the simulation time sequence data is denoised by the method of the application;

[0024] Fig. 2 (b) is a reconstructed component result graph after the simulation time sequence data is denoised by the method of the application;

[0025] Figure 3 is the result graph after the measured data is denoised by the method of the application. DETAILED DESCRIPTION

[0026] The specific embodiment of the application will be further described below in combination with the drawings.

[0027] The application provides a coordinate time sequence denoising method, and the specific process is as follows: Figure 1 ​​The original coordinate time sequence is firstly acquired and preprocessed, the preprocessed coordinate time sequence is converted into a trajectory matrix, and then singular value decomposition is performed on the coordinate time sequence based on the trajectory matrix to obtain decomposition components, and finally the reconstruction components are determined according to the MIC criterion, and the reconstructed coordinate time sequence is generated as the denoising result. The reconstruction components are determined by the MIC criterion, and the MIC values of the reconstruction sequence and the residual sequence are comprehensively considered to ensure that the reconstructed coordinate time sequence has a strong nonlinear relationship, while more decomposition components and less signal residues are considered, and scientific and effective denoising is realized.

[0028] Step 1. Data acquisition and preprocessing

[0029] The original coordinate time sequence data is firstly acquired, mainly including coordinate time sequence data generated by GNSS, VLBI, SLR and DORIS technologies.

[0030] Taking GNSS as an example, due to the influence of ionospheric delay, multipath effect, satellite cycle slip, receiver error, environmental interference signal and other factors in the satellite signal receiving process, the original coordinate time sequence obtained will exist in the case of gross error, step and data loss. In order to ensure the uniformity of sampling, when there is data loss, the data at the missing time can be obtained by interpolation, for example, the missing value interpolation can be realized by using the nearest neighbor method, cubic polynomial interpolation, singular value iterative difference method; In order to avoid this problem as much as possible, the gross error and step in the data need to be detected and eliminated, for example, the gross error and step can be detected and eliminated by using least square residual method, IQR robust estimation method, etc. This can eliminate part of the noise in advance, which is more conducive to subsequent processing; In addition, since the window length set by singular value decomposition for different components is different, in order to simplify the components of the coordinate time sequence and improve the accuracy of decomposition, it needs to be zero mean, which can eliminate the linear trend item by linear fitting, so that the processed data fluctuates around 0 value and is uniformly distributed.

[0031] Step 2. Singular value decomposition is performed to solve the decomposition components

[0032] Suppose the preprocessed coordinate time sequence is X=(x1,x2,x3,…,x N , select a suitable window length M, convert X into an M row L column trajectory matrix D, as shown in formula (1), wherein M should satisfy 1

[0033]

[0034] Based on the trajectory matrix D, the singular value decomposition is performed on the preprocessed coordinate time series, and the decomposition components are obtained, and the specific calculation process is as follows:

[0035] As can be seen from formula (1), the row and column of D are sub-sequences of X, so it can be expanded in a set of standard orthogonal bases: as shown in formula (2):

[0036]

[0037] In the formula, A is a k coefficient matrix composed of a k corresponding to the eigenvector E k , E k is the eigenvector corresponding to the kth eigenvalue in the covariance matrix of the trajectory matrix, and the covariance matrix T x of the trajectory matrix D is as shown in formula (3):

[0038]

[0039] In the formula, Solving T x , the eigenvalue and eigenvector E k can be obtained.

[0040] At the same time, a k can also be regarded as the result of filtering X, and the ith element can be calculated by formula (4):

[0041]

[0042] In the formula, is the jth element in E k , and is the ith element of a k .

[0043] Through the above formulas (2)-(4), the ith value of the kth decomposition component is:

[0044]

[0045] Among them, is the ith value of the kth reconstruction component, and D i is the trajectory matrix.

[0046] Step 3. Data reconstruction

[0047] The number of components used for data reconstruction is determined by the MIC criterion. The MIC value between the M components obtained by decomposition and time is calculated respectively, and the specific calculation method is as follows:

[0048] The set of data points formed by epoch t and coordinate time series data X (t i ,x i The data is distributed in a two-dimensional space (R,S). The data space is divided using a grid of m x n. The frequency of the data point falling in the (r,s)th grid is used as an estimate of P(r,s), i.e.

[0049]

[0050] Then, the MIC value of the time series X is calculated using formula (7).

[0051]

[0052] The resolution of the grid is limited to m×n < B, where B is set to the power of 0.6 of the data volume.

[0053] After calculating the MIC values ​​of the M decomposed components, the decomposed components are arranged in descending order of MIC values. Then, the top P decomposed components with the largest MIC values ​​are selected and summed to obtain the reconstructed sequence. Where P is the order of reconstruction (0 < P ≤ M), the remaining decomposed components form a residual sequence; calculate the reconstructed sequence X. P MIC value MIC values ​​of residual sequences The MIC value of the reconstructed sequence MIC values ​​of residual sequences The calculation method is the same as the MIC value calculation method of the above-mentioned decomposed components, except that the input data is converted into a reconstructed sequence and a residual sequence.

[0054] The result calculated from the above formula and The number of P values ​​can be determined, and the specific selection criteria for P are as follows: Greater than the first set threshold Less than the first set threshold or Greater than the first set threshold and The value is less than the second set threshold; wherein, in order to ensure the strong temporal correlation of the reconstructed sequence and the weak temporal correlation of the residual sequence, the value range of the first set threshold is [0.8, 1], and the value range of the second set threshold is [0, 0.3]. As another implementation, the first threshold and the second threshold can be set according to the requirements of the strength of the correlation of the extracted components and the quality of the obtained coordinate time series data.

[0055] After the number of P is determined according to the selection basis of P, the first P decomposition components with larger MIC values are reconstructed according to the MIC values of the decomposition components, and the reconstructed sequence is the denoised sequence of the coordinate time sequence.

[0056] In order to further illustrate the denoising effect of the method of the application on the coordinate time sequence, the method of the application is verified by using simulation data and measured data. The simulation data is composed of a plurality of periodic terms (see Table 1) and white noise (White Noise, WH) and flicker noise (Flicker Noise, FL) of different magnitudes.

[0057] Table 1: Periodic terms and parameters in the simulation time sequence

[0058]

[0059]

[0060] Simulation data experiment:

[0061] In order to verify the performance of the method of the application under different noise sizes and components, two groups of experiments are set up. In the first group of experiments, the noise components are controlled to be unchanged (WH and FL are 40% and 60% respectively), and the simulation time sequences with signal-to-noise ratios of 10dB, 5dB, 0dB, -5dB and -10dB are denoised respectively. In the second group of experiments, the noise intensity is controlled to be unchanged (the signal-to-noise ratio is 0dB), and the simulation time sequences with different noise components (WH+FL: 100%+0, 80%+20%, 60%+40%, 40%+60%, 20%+80%, 0+100%) are denoised. The denoising results of the simulation time sequence with a signal-to-noise ratio of 5dB and a noise component of 40% WH and 60% FL are shown in Figs. Figure 2(a) and 2(b) It can be seen that the method of the application can well realize the denoising of the simulation time sequence. At the same time, the method of the application is compared with the traditional SSA method. Table 2 shows the denoising effects of different methods on the simulation time sequence under the same noise component and different signal-to-noise ratios, and Table 3 shows the denoising effects of different methods on the simulation time sequence under the same signal-to-noise ratio and different noise components.

[0062] Table 2:

[0063]

[0064] Table 3:

[0065]

[0066]

[0067] Simulation results show that, compared to the traditional SSA denoising method which automatically selects reconstructed components by setting the contribution rate of cutoff eigenvalues, the denoising method of this invention has significant advantages. Under different noise components and intensities, the correlation and signal-to-noise ratio between the denoised signal and the real time series are greatly improved. Furthermore, the method of this invention achieves comparable denoising performance to the traditional SSA method which manually selects reconstructed components, while also enabling adaptive denoising, avoiding the influence of subjective factors. Moreover, the introduction of the MIC criterion makes the denoising effect more scientific.

[0068] Experimental data:

[0069] In the experiment using measured data, the ZIMM station from the IGS station was selected as the experimental object to verify the noise reduction effect of the method of the present invention. The experimental results are as follows: Figure 3 As shown (the first and second thresholds are 0.999 and 0.01, respectively).

[0070] from Figure 3 As can be seen, the method of the present invention achieves denoising of the measured coordinate time series. After calculation, the MIC value of the denoised time series is 0.999, the MIC value of the residual series is 0.07, and the eigenvalue contribution rate of the reconstructed series can reach 83%. While ensuring that the reconstructed coordinate time series has a strong nonlinear relationship, it also takes into account more decomposition components and less signal residue, thus achieving scientific and effective denoising.

Claims

1. A method for denoising coordinate time series data, applied in GNSS, characterized in that, The method includes the following steps: 1) Receive GNSS satellite signals and obtain the raw coordinate time series; 2) Convert the acquired coordinate time series into the corresponding trajectory matrix, and perform singular value decomposition on the coordinate time series based on the converted trajectory matrix to obtain the decomposed components; 3) Calculate the MIC value between each decomposition component and time, obtain the MIC value of each decomposition component, and sort the decomposition components according to the size of the MIC value of each decomposition component. 4) Select the top P decomposition components with the largest MIC values ​​and sum them to generate the first reconstructed sequence. Then, select the top P+1 decomposition components with the largest MIC values ​​to form the second reconstructed sequence. Finally, select the remaining decomposition components after selecting the top P components to form the residual sequence. Calculate the MIC value of the first reconstructed sequence. MIC value of the second reconstructed sequence MIC values ​​of residual sequences When the value of P satisfies either condition one or condition two, the first reconstructed sequence is the denoised coordinate time series, thus achieving denoising of the original coordinate time series; otherwise, the value of P is incremented by 1, and step 4) is repeated; where condition one is: Greater than the first set threshold and Less than the first set threshold; condition two is: Greater than the first set threshold and It is less than the second set threshold.

2. The noise reduction method for coordinate time series according to claim 1, characterized in that, The first set threshold in step 3) has a value range of [0.8, 1]; the second set threshold has a value range of [0, 0.3].

3. The method for denoising coordinate time series according to claim 1, characterized in that, Before the singular value decomposition in step 2), the original coordinate time series data needs to be preprocessed. The preprocessing includes the detection and removal of gross errors and steps, as well as the interpolation of missing values.

4. The noise reduction method for coordinate time series according to claim 3, characterized in that, The preprocessing also includes zero-mean normalization of the original coordinate time series.

5. The method for denoising coordinate time series according to claim 3 or 4, characterized in that, The trajectory matrix is ​​an M-row, L-column matrix, where 1 < M < N / 2, L = N - M + 1, and N is the length of the preprocessed coordinate time series.

6. The method for denoising coordinate time series according to claim 1, 3, or 4, characterized in that, The formula for calculating the decomposed components is: in, x is the i-th value of the k-th reconstructed component. i Let D be the i-th element in the preprocessed coordinate time series. i Let E be the trajectory matrix. k Let a be the eigenvector corresponding to the k-th eigenvalue in the covariance matrix of the trajectory matrix. k For the corresponding feature vector E k coefficient, It is E k The j-th element in It is a k The i-th element.

Images