CN111142134B - Coordinate time series processing method and device - Google Patents
Coordinate time series processing method and device Download PDFInfo
- Publication number
- CN111142134B CN111142134B CN201911102015.8A CN201911102015A CN111142134B CN 111142134 B CN111142134 B CN 111142134B CN 201911102015 A CN201911102015 A CN 201911102015A CN 111142134 B CN111142134 B CN 111142134B
- Authority
- CN
- China
- Prior art keywords
- residual
- sequences
- principal component
- sequence
- new
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
- G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Image Analysis (AREA)
Abstract
The invention discloses a coordinate time series processing method and a device, wherein the method comprises the following steps: acquiring three residual error sequences corresponding to coordinate time sequences in three directions corresponding to the data from a target measuring station; calculating the main components of three residual error sequences; dividing the principal component into a first principal component and a second principal component, wherein the first principal component is used for representing a first part of white noise, and the second principal component is used for representing the colored noise and a second part of white noise; deleting the second principal component from the principal component, and reserving the first principal component to obtain an updated principal component; based on the updated principal components, three new residual sequences are calculated. By adopting the scheme of the invention, colored noise in the residual sequence can be eliminated, and three new residual sequences are obtained only through white noise, so that the error represented by the new residual sequences is smaller, and the accurate uncertainty of the GNSS reference station motion parameter estimation can be obtained.
Description
Technical Field
The invention relates to the field of GNSS data processing, in particular to a coordinate time series processing method and device of a GNSS observation station.
Background
GNSS reference station coordinate time series are widely used in geodetic and geodynamic studies, and contain white noise and colored noise (errors) in the GNSS station coordinate time series due to various factors including the geophysical environment and systematic errors associated with GNSS technology. Since noise in the GNSS coordinate time series affects the stability and reliability of the stations' coordinates, leading to misinterpretation of some geophysical phenomena, it is necessary to take measures to weaken the noise strength of the GNSS coordinate time series.
However, the current denoising method cannot remove the colored noise existing in the low frequency, and even if the intensity of the noise is weakened, the estimation of the noise amplitude, especially the residual colored noise, cannot be accurately obtained.
Disclosure of Invention
In view of the above, the main objective of the present invention is to provide a method and an apparatus for processing a coordinate time series, by which colored noise can be eliminated, and accurate motion parameters and uncertainty of the coordinate time series can be obtained.
In order to achieve the purpose, the technical scheme of the invention is realized as follows:
a method of coordinate time series processing, the method comprising:
acquiring data from a target survey station and preprocessing the data to obtain three residual error sequences corresponding to coordinate time sequences of the data in three directions; the residual error sequence represents an error between the actual value of the coordinate time sequence and the estimated value of the coordinate time sequence, and the error consists of white noise and colored noise;
calculating characteristic values and principal components of the three residual error sequences;
dividing the principal component into a first principal component and a second principal component based on the eigenvalue, wherein the first principal component is used for characterizing a first part of white noise, the second principal component is used for characterizing mixed noise of the colored noise and a second part of white noise, and the first part of white noise is different from the second part of white noise; deleting the second principal component from the principal component, and reserving the first principal component to obtain an updated principal component;
calculating three new residual sequences based on the updated principal component;
and replacing the original residual sequence with the three new residual sequences to obtain three new coordinate time sequences.
In the above scheme, the calculating the eigenvalues and principal components of the three residual sequences includes:
performing Fourier transform on the three residual error sequences to obtain three transformed residual error sequences;
calculating power spectrums of the three corresponding residual error sequences according to the three transformed residual error sequences;
carrying out standardization processing on the power spectrums of the three residual error sequences to obtain three standard residual error sequence power spectrums; the normalization processing is used for reducing the proportion of the power spectrum value of the residual sequence with the maximum power spectrum value in the sum of the power spectrum values of the three residual sequences from a first proportion to a second proportion to obtain an adjusted residual sequence; wherein the first ratio is greater than the second ratio;
and calculating eigenvalues and eigenvector matrixes of the three standard residual sequence power spectrums, and calculating principal components of the three standard residual sequence power spectrums based on the eigenvector matrixes.
In the above scheme, the calculating three new residual sequences based on the updated principal component includes:
calculating to obtain three new standard residual error sequence power spectrums according to the characteristic vector matrix and the updated principal component;
performing anti-standardization processing on the three new standard residual error sequence power spectrums to obtain three new anti-standardization residual error sequence power spectrums; wherein the de-normalization process is used to increase the power spectrum values of the adjusted residual sequence from the second scale to a third scale; wherein the second ratio is less than the third ratio;
carrying out inverse Fourier transform on the three new anti-standard residual sequence power spectrums to obtain three new residual sequences; wherein, the three new residual error sequences can represent the errors between the three coordinate time series actual values and the three coordinate time series estimated values.
In the foregoing solution, before performing inverse fourier transform on the three new anti-standard residual sequence power spectrums, the method further includes:
and calculating the phase corresponding to each frequency signal on each residual sequence power spectrum according to each residual sequence power spectrum obtained when the three residual sequences are subjected to Fourier transform.
In the above scheme, the performing inverse fourier transform on the three new anti-standard residual sequence power spectrums includes:
and performing inverse Fourier transform on the three new anti-standard residual sequences based on the phase corresponding to the frequency signal on the power spectrum of each of the three residual sequences.
In the above aspect, the apparatus includes:
the acquisition module is used for acquiring data from a target survey station and preprocessing the data to obtain three residual error sequences corresponding to coordinate time sequences of the data in three directions; the residual error sequence represents an error between the actual value of the coordinate time sequence and the estimated value of the coordinate time sequence, and the error consists of white noise and colored noise;
the principal component analysis module is used for calculating characteristic values and principal components of the three residual error sequences; dividing the principal component into a first principal component and a second principal component based on the eigenvalue, wherein the first principal component is used for characterizing a first part of white noise, the second principal component is used for characterizing mixed noise of the colored noise and a second part of white noise, and the first part of white noise is different from the second part of white noise; deleting the second principal component from the principal component, and reserving the first principal component to obtain an updated principal component; calculating three new residual sequences based on the updated principal component;
and the replacing module is used for replacing the original residual sequence with the three new residual sequences to obtain three new coordinate time sequences.
In the above solution, the apparatus further comprises: a power spectrum calculating module and a standardization module;
the calculation power spectrum module is used for performing Fourier transform on the three residual error sequences to obtain three transformed residual error sequences; calculating power spectrums of the three corresponding residual error sequences according to the three transformed residual error sequences;
the standardization module is used for standardizing the power spectrums of the three residual error sequences to obtain three standard residual error sequence power spectrums; the normalization processing is used for reducing the proportion of the power spectrum value of the residual sequence with the maximum power spectrum value in the sum of the power spectrum values of the three residual sequences from a first proportion to a second proportion to obtain an adjusted residual sequence; wherein the first ratio is greater than the second ratio;
the principal component analysis module is further used for solving eigenvalues and eigenvector matrixes of the three standard residual sequence power spectrums, and calculating principal components of the three standard residual sequence power spectrums based on the eigenvector matrixes.
In the above solution, the apparatus further comprises: an anti-normalization module;
the principal component analysis module is also used for calculating three new standard residual sequence power spectrums according to the characteristic vector matrix and the updated principal components;
the anti-standardization module is used for carrying out anti-standardization processing on the three new standard residual error sequences to obtain three new anti-standard residual error sequences; wherein the de-normalization process is used to increase the power spectrum values of the adjusted residual sequence from the second scale to a third scale; wherein the second ratio is less than the third ratio;
the power spectrum calculation module is also used for carrying out inverse Fourier transform on the three new anti-standard residual sequence power spectrums to obtain three new residual sequences; wherein, the three new residual error sequences can represent the errors between the three coordinate time series actual values and the three coordinate time series estimated values.
In the above scheme, the power spectrum calculating module is further configured to calculate a phase corresponding to each frequency signal in each residual sequence power spectrum according to each residual sequence power spectrum obtained when the fourier transform is performed on the three residual sequences.
In the foregoing solution, the power spectrum calculation module is further configured to perform inverse fourier transform on the three new anti-standard residual sequence power spectrums based on a phase corresponding to each frequency signal on each of the residual sequence power spectrums in the three residual sequences.
According to the coordinate time sequence processing method and device provided by the embodiment of the invention, three residual error sequences corresponding to coordinate time sequences in three directions corresponding to data are obtained from a target test station, and main components of the three residual error sequences are calculated; deleting a second principal component of mixed noise representing colored noise and partial white noise, and reserving the first principal component representing the white noise to obtain an updated principal component; based on the updated principal components, three new residual sequences are calculated. Therefore, colored noise in the residual sequence can be removed, three new residual sequences only containing white noise are obtained, the error represented by the new residual sequences is smaller, and the accurate uncertainty of the GNSS reference station motion parameter estimation can be obtained.
Drawings
FIG. 1 is a prior art flow chart of removing GNSS coordinate time series noise based on wave transform coefficient information entropy and modulo maximum;
FIG. 2 is a diagram illustrating the result of removing noise in a GNSS coordinate time series based on empirical mode decomposition provided in the prior art;
FIG. 3 is a flowchart of a coordinate time series processing method according to an embodiment of the present invention;
FIG. 4 is a flowchart of another method for processing a time series of coordinates according to an embodiment of the present invention;
FIG. 5 is a schematic diagram of a time series of coordinates in the direction N, E, U according to an embodiment of the present invention;
fig. 6 is a schematic diagram of a residual sequence in the N, E, U direction according to an embodiment of the present invention;
fig. 7 is a schematic diagram of a fourier transformed and logarithmized residual sequence in the N, E, U direction according to an embodiment of the present invention;
fig. 8 is a schematic diagram of a power spectrum of a new residual sequence in the N, E, U direction obtained after principal component analysis according to an embodiment of the present invention;
fig. 9 is a schematic diagram of power spectrums of three residual error sequences in the N, E, U direction after being subjected to principal component analysis processing according to an embodiment of the present invention;
fig. 10 is a block diagram of a coordinate time series processing apparatus according to an embodiment of the present invention.
Detailed Description
The following methods are mainly used for processing the GNSS survey station coordinate time series in the related art:
for a single station, some scholars propose removing noise in GNSS coordinate time series based on wavelet transformation. A common method is discrete wavelet transform, the method carries out multi-resolution decomposition on a GNSS coordinate time sequence based on Mallat algorithm, then carries out threshold quantization on the obtained high-frequency wavelet transform coefficient, and finally reconstructs the processed wavelet coefficient to obtain the de-noised coordinate time sequence. However, this method only considers white noise when setting the threshold, and has a good effect of attenuating high-frequency noise (white noise) in the coordinate time series, and the denoising effect is not ideal for low-frequency noise existing in the coordinate time series, such as flicker noise and random walk noise.
Denoising based on wavelet entropy and modulus maximum of wavelet transform coefficient. The method obtains the wavelet coefficient of the flicker noise according to the wavelet coefficient of the simulated white noise and the wavelet entropy of the wavelet coefficient of the coordinate time sequence, then reconstructs the wavelet coefficient of the flicker noise into the flicker noise and eliminates the flicker noise from the original coordinate time sequence. Finally, the remaining white noise is weakened by using a general discrete wavelet transform, and the specific processing flow is shown in fig. 1. The method needs to know the signal-to-noise ratio of white noise and the rest signal in advance to obtain the variance of the white noise, and can simulate a white noise sequence and obtain the wavelet transform coefficient thereof on the basis; in addition, when the method obtains the wavelet entropy of the colored noise through the wavelet entropy of the white noise, the influence of the wavelet coefficient modulus maximum of the periodic signal in the time sequence and the white noise existing in the low frequency are ignored. Therefore, the method is only suitable for short coordinate time series, otherwise, the periodic signal is likely to be influenced.
The noise in the GNSS coordinate time series is removed based on an Empirical Mode Decomposition (EMD) method. The method is similar to discrete wavelet transform, an original signal is decomposed into a series of Intrinsic Function Mode components (IMF) with frequencies from high to low, a residual error item is combined with adjacent low-frequency IMF components, and a nonlinear trend item of a sequence is judged through indexes such as IMF component variance contribution rate and the like. The low and medium frequency IMF components of the approach trend term are obvious periodic motion terms of the signal, and the high frequency part can be a signal high frequency periodic term which cannot be identified temporarily or noise information contained in the signal. By using the method, the nonlinear signal in the GNSS coordinate time series can be effectively extracted, and specific results can be seen in fig. 2. However, the power of the colored noise is concentrated in the low frequency, and the power of the white noise is maintained at the same level in the whole power spectrum, so the white noise and the colored noise existing in the low frequency cannot be removed by the EMD method.
And processing the GNSS survey station coordinate time sequence by using a singular spectrum analysis method. The method constructs a track matrix according to the observed time sequence, decomposes and reconstructs the track matrix, thereby extracting signals representing different components of the original time sequence, such as long-term trend signals, periodic signals, noise signals and the like, and analyzing the structure of the time sequence. Experimental results show that the sequence reconstructed by the singular spectrum analysis method is smoother than the original observed sequence, and has an obvious noise reduction effect. However, this method still cannot effectively separate the low frequency colored noise from the periodic signal.
In view of the foregoing problems, an embodiment of the present invention provides a coordinate time series processing method, as shown in fig. 3, and fig. 3 is a flowchart of a coordinate time series processing method provided in an embodiment of the present invention. The specific process is as follows:
s101, acquiring data from a target survey station and preprocessing the data to obtain three residual error sequences corresponding to coordinate time sequences of the data in three directions; and the residual error sequence represents the error between the actual value of the coordinate time sequence and the estimated value of the coordinate time sequence, and the error consists of white noise and colored noise.
The target survey station in the embodiment of the invention can be a GNSS survey station. A time series of coordinates in three directions in the time domain, here specifically N, E, U, which may represent three directions of three-dimensional coordinate axes, is obtained from the target station.
The coordinate time sequence of the single component of the GNSS single station can be written as the following motion model:
in the formula (1), x(1)+x(2)tiIn the form of a linear term, the term,
is a non-linear term; y (t)i) The coordinate time sequence of the GNSS single station and single component is obtained; t is tiIs time, i is 1, 2, 3 … m, m is epoch number; x is the number of(1)And x(2)Initial position and velocity, respectively; q-1 is the number of periodic signals contained in the time series; x is the number of(2k-1)And x(2k)Coefficients of a harmonic function describing the amplitude of the periodic signal; omegakThe/2 pi is the frequency of the signal; v. ofiIs an error, i.e., noise.
The matrix form of equation (1) is as follows:
Y＝BX+v (2)
in the above formula, Y is a time series of coordinates of three components N, E, U in this embodiment of the present invention. The matrix B is called a design matrix and can be obtained by formula (1), X is a parameter matrix to be estimated, and v is a residual sequence.
The specific form is as follows:
Y＝[y(t1),y(t2),y(t3),…,y(tm)]T，
and constructing a Y matrix and a B matrix according to the time corresponding to the GNSS observation station coordinate time sequence and the period of the nonlinear motion contained in the coordinate time sequence. In the case of white noise alone, the motion parameter X of the station is estimated by means of least squares, by equation (2), i.e. X ═ BTB)-1BTAnd Y. And then substituting the parameter estimation value X into the formula (1), and removing a linear term and a nonlinear term to obtain a residual sequence v ═ v [ v ]1,v2,v3,…vm]。
In the embodiment of the present invention, since the coordinate time series is a component in the three coordinate axis directions, m here is 3.
And S102, calculating characteristic values and principal components of the three residual sequences.
The embodiment of the invention calculates the characteristic value and the principal component through the obtained residual sequence, and specifically comprises the following steps:
and performing Fourier transform on the three residual error sequences to obtain three transformed residual error sequences.
The above residual sequence is substituted into formula (3),
py＝fft(v)(3)
in equation (3), fft is the fast fourier transform function in Matlab, py ═ py1,py2,py3,…,pym]Wherein py isiIs a plurality;
calculating power spectrums of the three corresponding residual error sequences according to the three transformed residual error sequences, specifically:
and (4) calculating the corresponding power spectrum of the three transformed residual sequences through a formula (4).
Py＝abs(py) (4)
Wherein abs is a function of the arithmetic square root of the sum of the squares of the real and imaginary parts of the complex number in Matlab, and Py ═ Py1,Py2,Py3,…,Pym]Is the power spectrum of the residual sequence.
Because the power curve of the residual sequence fluctuates greatly, the power and frequency of the three components are logarithmic with the base 10, so that the power spectrum curve changes slowly. Then, the residual power spectrum after the logarithm of the three components is formed into a matrix p ═ PyN,PyE,PyU]In which Py isN，PyE，PyULog power spectra of the residual sequence for north, east and vertical directions, respectively.
Carrying out standardization processing on the power spectrums of the three residual error sequences to obtain three standard residual error sequence power spectrums; the normalization processing is used for reducing the proportion of the power spectrum value of the residual sequence with the maximum power spectrum value in the sum of the power spectrum values of the three residual sequences from a first proportion to a second proportion to obtain an adjusted residual sequence; wherein the first ratio is greater than the second ratio.
In the examples of the present invention, since PyUIs far greater than PyNAnd PyEIf one proceeds directly from the covariance matrix of the three power spectra to perform principal component analysis, PyUWill obviously play a dominant role, while the other three power spectra are hard to be reflected in the principal component. In order to avoid that the extracted common noise information is dominated by noise information of a certain direction, the power spectra of the three residual sequences should be normalized. The z-score (zero-mean normalization) method is used herein without limitation. The specific formula is as follows:
wherein
Normalized power spectrum, p, for ith row and jth columnijIs PyNOr PyEOr PyUThe elements of row i and column j,
and σjAre each PyNOr PyEOr PyUMean and standard deviation of column j.
The power spectrum of each residual sequence after the standardization is transformed by using a formula (5), so that the value range of the power spectrum before transformation is far larger than PyNAnd PyEPy of (2)UThe ratio of the actual power spectrum values in the sum of all power spectrum values is changed from a first ratio to a second ratio, wherein the first ratio is much larger than the second ratio. It is understood that, prior to the change, PyUThe ratio of the power spectrum value of (2) to the total power spectrum value of the three components is the maximum and is a first ratio; after equation (5), the modified PyUCompared with the changed PyNAnd PyENot as large as before, so PyUThe proportion of power spectral values of (a) to the sum of all power spectral values over the three components is reduced by a little, the proportion at this time being recorded as a second proportion, the second proportion being smaller than the first proportion.
And calculating eigenvalues and eigenvector matrixes of the three standard residual sequence power spectrums, and calculating principal components of the three standard residual sequence power spectrums based on the eigenvector matrixes.
Calculating a variance-covariance matrix after power spectrum standardization of three residual error sequences in the embodiment of the invention, and solving a characteristic value lambda of the matrix [ lambda ═ lambda [ [ lambda ]1,λ2,λ3]And the eigenvector matrix V ═ V1,V2,V3]. Wherein the elements in λ are arranged in descending order, ViIs λiThe corresponding feature vector. The principal component was then calculated as follows:
a＝p*×V (6)
wherein p is*Is composed of
Corresponding normalized power spectrum, principal component a ═ a1，a2，a3]，aiIs the ith main component of m × 1.
S103, dividing the principal component into a first principal component and a second principal component based on the characteristic value, wherein the first principal component is used for representing the white noise, and the second principal component is used for representing mixed noise of the colored noise and partial white noise; and deleting the second principal component from the principal components, and reserving the first principal component to obtain the updated principal component.
In the embodiment of the present invention, the original three residual sequences are divided into 3 principal components by formula (6), and since colored noise occupies most of the principal components and the remaining components are white noise in the power spectrum of the residual sequences, the 3 principal components can be divided into the first principal component and the second principal component by at least three ways:
the first method is as follows: due to lambdaiAre elements arranged in descending order, so the last principal component can be directly selected as the first principal component, while the first three principal components are considered as a whole as the second principal component.
The second method comprises the following steps: and selecting k main components with the ratio of the k characteristic values in the sum of the characteristic values lower than a set threshold value, such as k main components with the ratio lower than 15%, as first main components, and the rest main components as second main components.
The first principal component corresponds to white noise, and the second principal component corresponds to mixed noise of colored noise and partial white noise. And removing the second principal component, reserving the first principal component, namely only reserving white noise, and taking the reserved first principal component as an updated principal component.
And S104, calculating three new residual sequence power spectrums based on the updated principal components.
And calculating to obtain three new standard residual error sequence power spectrums according to the characteristic vector matrix and the updated principal component.
For the updated first principal component, the power spectrum of the new residual sequence is calculated using equation (7).
p0＝[a2,a3]×[V2,V3]T (7)
Wherein, when the last of the three components is selected as the first main component, formula (7) is p0＝[a3]×[V2,V3]T，
Performing anti-standardization processing on the three new standard residual error sequence power spectrums to obtain three new anti-standardization residual error sequence power spectrums; wherein the de-normalization process is used to increase the power spectrum values of the adjusted residual sequence from the second scale to a third scale; wherein the second ratio is less than the third ratio.
Since the power spectrum has been normalized before the principal component analysis, p must be denormalized according to equation (8) after the new power spectrum of the standard residual sequence is obtained. Namely:
wherein the content of the first and second substances,
for the element power of ith row and j column of each component in three directions after the denormalization,
for the element values in the ith row and j column in the updated standard residual error sequence power spectrum,
and σjRespectively the mean value and the standard deviation of the jth column of the updated standard residual sequence power spectrum,
the finally obtained logarithmic power spectrums of the residual sequences in the north direction, the east direction and the vertical direction after the denormalization processing of the formula (8) are respectively obtained.
It can be understood that formula (8) is an inverse transformation of formula (5), and the ratio of the power spectrum value corresponding to the adjusted residual sequence to the total power spectrum value of the three new anti-standard residual sequences increases from the previous second ratio to the third ratio after the new anti-standard residual sequence power spectrum is subjected to the anti-standardization processing of formula (8).
Carrying out inverse Fourier transform on the three new anti-standard residual sequence power spectrums to obtain three new residual sequences; wherein, the three new residual error sequences can represent the errors between the three coordinate time series actual values and the three coordinate time series estimated values.
Before performing inverse fourier transform on the three new anti-standard residual sequence power spectrums, the method further includes:
and calculating the phase corresponding to each frequency signal on each residual sequence power spectrum according to each residual sequence power spectrum obtained when the three residual sequences are subjected to Fourier transform.
In order to more accurately restore the corresponding residual sequence according to the power spectrum, the phase of the signal corresponding to each frequency in the power spectrum is obtained first. According to the result of the formula (3), the phase of the signal at each frequency on the corresponding component is obtained by using the angle function in the Matlab tool, in the case of N components:
θN＝angle(pyN) (9)
py in formula (9)NFor the power spectra of the three components before principal component analysis, θN＝[θNf1,θNf2,θNf3,…,θNfm]Wherein thetaNfiFor frequency f of component NNi(i ═ 1, 2, … m) the corresponding phase of the signal.
The inverse fourier transform of the three new anti-standard residual sequences includes:
and performing inverse Fourier transform on the three new anti-standard residual sequence power spectrums based on the phase corresponding to each frequency signal on each residual sequence power spectrum in the three residual sequences.
Since the power spectrum curves of the three residual error sequences are changed slowly before the principal component analysis, the power spectrum and the frequency are logarithmic with a base 10, and an index of 10 is required to be taken for each power spectrum and frequency,
the phase and power corresponding to all frequencies calculated by the above equation (9) are subjected to inverse fourier transform in combination with equations (10) and (11). Taking the N component as an example, the method specifically includes:
in the above formula, i represents an imaginary number, PYN＝[PYN1，PYN2，PYN3，…PYNm]。
Representing the new residual sequence of the N component. ifft is the inverse fourier transform function in the Matlab tool.
And S105, replacing the original residual sequence with the three new residual sequences to obtain three new coordinate time sequences.
Adding the new residual sequence obtained after the processing in the steps into the linear term and the nonlinear term in the formula (1), replacing the previous residual sequence, and finally forming a new coordinate time sequence of the GNSS observation station
After the residual sequence is processed, the obtained new residual sequence only represents white noise, and at the moment, the accurate uncertainty of the GNSS reference station motion parameter estimation can be directly obtained by using a least square estimation method:
the complete process described above can also be implemented by the coordinate time series processing flow provided in fig. 4.
According to the coordinate time sequence processing method and device provided by the embodiment of the invention, three residual error sequences corresponding to coordinate time sequences in three directions corresponding to data are obtained from a target test station, and main components of the three residual error sequences are calculated; deleting a second principal component representing mixed noise of the colored noise and partial white noise, and reserving the first principal component representing the white noise to obtain an updated principal component; based on the updated principal components, three new residual sequences are calculated. Therefore, colored noise in the residual sequence can be removed, three new residual sequences only containing white noise are obtained, the error represented by the new residual sequences is smaller, and the accurate uncertainty of the GNSS reference station motion parameter estimation can be obtained.
An example of a coordinate time series processing method is provided in the embodiments of the present invention, as shown in fig. 5, fig. 5 is a schematic diagram of a time series of coordinates in N, E, U three directions.
FIG. 5 is a time series of coordinates for obtaining the AC62 station 2011-1-2016-12-N, E, U components from the https:// www.unavco.org/website. The apparent "jump" in the time series due to antenna replacement or earthquake is estimated and removed. Then, the gross error is removed by using an error criterion called triple, and finally a clean coordinate time sequence is obtained.
According to the motion model of the single component of the GNSS single station, namely formula (1) and the coordinate time series thereof, under the condition of only considering white noise and annual and semiannual motion of the survey station, the motion speed (linear motion), the annual and semiannual signal (nonlinear motion) amplitude of the three components of the target survey station are respectively estimated by using general least squares. And substituting the estimated values of the motion parameters and the coordinate time sequence into the formula (1) to obtain linear terms and nonlinear terms in the three-component time sequence, and further obtaining a residual sequence of the three components, wherein the residual sequence is shown in fig. 6, and the variation amplitude of each residual sequence is relatively large at the moment.
And (4) performing fast Fourier transform on the three residual sequences, and obtaining the phase of the signal corresponding to each frequency according to a formula (3). Meanwhile, the power spectrums of the three residual sequences are obtained according to the formula (4), then the logarithm with the base number of 10 is taken for the power spectrums, the result is shown in fig. 7, and then the power spectrums after the logarithm is taken for three parts form a matrix p ═ PyN,PyE,PyU]. In conjunction with fig. 6, it can be seen that the variation amplitude of each residual sequence after fourier transform and logarithm extraction in fig. 7 is significantly smaller than that of the residual sequence in fig. 6.
The power spectrum matrix is normalized according to equation (5). Then, the variance-covariance matrix of the matrix is calculated, and the eigenvalue lambda is obtained1,λ2,λ3]And the eigenvector matrix V ═ V1,V2,V3]. The principal component after power spectrum normalization is calculated according to equation (6). Due to lambda1＝1.9，λ2＝0.61，λ3Therefore, the 3 rd principal component is selected as the first principal component, i.e., the updated principal component, and is reconstructed according to the formula (7), and then is denormalized according to the formula (8), and the result is shown in fig. 8. It can be seen that the power spectra of the three residual sequences obtained by the reconstruction are stable, and the frequency is used as an independent variable and the power is used as the powerAnd (4) performing line fitting on the dependent variable to obtain the slopes of N, E, U components of 0.0217, 0.0243 and-0.0452 respectively, namely obtaining the spectral index of the new residual sequence of substantially 0, which also proves that the extracted new residual sequence is almost white noise.
The logarithmic power spectrum of the new standard residual sequence calculated according to the principal component analysis is calculated as an index of 10, inverse fourier transform is performed according to equations (9) to (11) using the phase obtained by equation (3) according to the obtained power spectrum before logarithm, and the obtained three residual sequences are new residual sequences from which colored noise has been removed, as shown in fig. 9. It can be seen that the fluctuation range of the three residual sequences is significantly reduced compared to fig. 6.
And adding the obtained new residual sequence with the previously obtained linear term and nonlinear term, namely replacing the original residual sequence in the original coordinate time sequence with the new residual sequence to obtain three new coordinate time sequences. Since the obtained new residual sequence contains almost no colored noise, the accurate motion parameters and uncertainty of the measuring station are obtained by using a general least square estimation method according to the formulas (12) - (14) under the condition of only considering white noise, and the result is shown in table 2. To analyze the changes in uncertainty in GNSS station motion parameters before and after the principal component processing, the motion parameters and their uncertainty of the AC62 station before the principal component processing were calculated using the vector software, taking into account the effects of colored noise, and the results before and after the processing were compared as described in table 1. It can be seen that before denoising by using the principal components, the uncertainty of the parameter estimation is larger; after colored noise is removed, the uncertainty of the parameter is significantly reduced and the estimated variation of the parameter is small. Therefore, the method can remove the influence of colored noise and improve the signal-to-noise ratio of the coordinate time sequence and the uncertainty of parameter estimation under the condition of less influence on the estimation of the motion parameters of the station.
Table 1 estimation and uncertainty of motion parameters for AC62 station before principal component processing
Table 2 principal component processed AC62 station motion parameter estimation and uncertainty
According to the coordinate time sequence processing method and device provided by the embodiment of the invention, three residual error sequences corresponding to coordinate time sequences in three directions corresponding to data are obtained from a target test station, and main components of the three residual error sequences are calculated; deleting a second principal component representing mixed noise of the colored noise and partial white noise, and reserving the first principal component representing the white noise to obtain an updated principal component; based on the updated principal components, three new residual sequences are calculated. Therefore, colored noise in the residual sequence can be removed, three new residual sequences can be obtained, the error represented by the new residual sequences is smaller, and the accurate uncertainty of the GNSS reference station motion parameter estimation can be obtained.
An embodiment of the present invention provides a coordinate time series processing apparatus, as shown in fig. 10, and as shown in fig. 10, the apparatus 10 includes: an acquisition module 11, a principal component analysis module 12 and a replacement module 13, in particular,
the acquisition module 11 is used for acquiring data from a target survey station and preprocessing the data to obtain three residual error sequences corresponding to coordinate time sequences of the data in three directions; and the residual error sequence represents the error between the actual value of the coordinate time sequence and the estimated value of the coordinate time sequence, and the error consists of white noise and colored noise.
And (3) constructing a matrix Y and a matrix B according to a formula (1), substituting a formula (2) by using a least square estimation method to calculate to obtain a motion parameter X, and obtaining a residual sequence according to the formula (1).
The principal component analysis module 12 is used for calculating characteristic values and principal components of the three residual error sequences; dividing the principal component into a first principal component and a second principal component based on the eigenvalue, wherein the first principal component is used for representing the white noise, and the second principal component is used for representing mixed noise of the colored noise and partial white noise; deleting the second principal component from the principal component, and reserving the first principal component to obtain an updated principal component; based on the updated principal components, three new residual sequences are calculated.
The method comprises the steps of carrying out Fourier transform on the residual sequence to obtain a power spectrum, then taking a logarithm with the base of 10 on the obtained power spectrum to enable the fluctuation of a power spectrum curve to be slow, and then carrying out standardization processing on the obtained power spectrum after taking the logarithm to avoid the phenomenon that the value range of the power spectrum of the residual sequence in a certain direction is too large to influence the functions of the power spectrums of components in other three directions in common noise.
And calculating a characteristic value and a characteristic vector matrix of the normalized power spectrum of the residual sequence, and calculating a corresponding principal component. The method comprises the steps of dividing the main components according to the characteristic values to obtain a first main component representing white noise and a second main component representing mixed noise of colored noise and white noise, removing the second main component, reserving the first main component to obtain an updated main component, and understanding that the updated main component represents the white noise. And calculating to obtain a new standard residual sequence according to the obtained updated principal component, and then performing anti-standardization processing and inverse Fourier transform on the new residual sequence. In order to restore the corresponding residual sequence more accurately according to the power spectrum, the phase of the signal corresponding to each frequency on the power spectrum of each residual sequence when the fourier transform is performed before is combined with the index of which the logarithmic power spectrum is used as 10 to perform the fourier inverse transform, and finally, a new residual sequence is obtained.
And the replacing module 13 is configured to replace the original residual sequence with the three new residual sequences to obtain three new coordinate time sequences.
And adding the obtained new residual sequence with the linear term and the nonlinear term in the formula (1) to obtain an updated coordinate time sequence. Because the new residual sequence in the current coordinate time sequence only corresponds to white noise and colored noise is reduced, the new coordinate time sequence can directly obtain the accurate uncertainty of the GNSS reference station motion parameter estimation by using a least square estimation method, and the uncertainty at the time is much smaller than the uncertainty before the principal component analysis, which can be specifically referred to the corresponding values in the above tables 1 and 2.
According to the coordinate time sequence processing method and device provided by the embodiment of the invention, three residual error sequences corresponding to coordinate time sequences in three directions corresponding to data are obtained from a target test station, and main components of the three residual error sequences are calculated; deleting a second principal component representing mixed noise of the colored noise and partial white noise, and reserving the first principal component representing the white noise to obtain an updated principal component; based on the updated principal components, three new residual sequences are calculated. Therefore, colored noise in the residual sequence can be removed, three new residual sequences only containing white noise are obtained, the error represented by the new residual sequences is smaller, and the accurate uncertainty of the GNSS reference station motion parameter estimation can be obtained.
The above description is only a preferred embodiment of the present invention, and is not intended to limit the scope of the present invention.
Claims (8)
1. A method for processing a time series of coordinates, the method comprising:
acquiring data from a target survey station and preprocessing the data to obtain three residual error sequences corresponding to coordinate time sequences of the data in three directions; the residual error sequence represents an error between the actual value of the coordinate time sequence and the estimated value of the coordinate time sequence, and the error consists of white noise and colored noise;
performing Fourier transform on the three residual error sequences to obtain three transformed residual error sequences;
calculating power spectrums of the three corresponding residual error sequences according to the three transformed residual error sequences;
carrying out standardization processing on the power spectrums of the three residual error sequences to obtain three standard residual error sequence power spectrums; the normalization processing is used for reducing the proportion of the power spectrum value of the residual sequence with the maximum power spectrum value in the sum of the power spectrum values of the three residual sequences from a first proportion to a second proportion to obtain an adjusted residual sequence; wherein the first ratio is greater than the second ratio;
calculating eigenvalues and eigenvector matrixes of the three standard residual sequence power spectrums, and calculating principal components of the three standard residual sequence power spectrums based on the eigenvector matrixes;
dividing the principal component into a first principal component and a second principal component based on the eigenvalue, wherein the first principal component is used for characterizing a first part of white noise, the second principal component is used for characterizing mixed noise of the colored noise and a second part of white noise, and the first part of white noise is different from the second part of white noise; deleting the second principal component from the principal component, and reserving the first principal component to obtain an updated principal component;
calculating three new residual sequences based on the updated principal component;
and replacing the original residual sequence with the three new residual sequences to obtain three new coordinate time sequences.
2. The method of claim 1, wherein computing three new residual sequences based on the updated principal component comprises:
calculating to obtain three new standard residual error sequence power spectrums according to the characteristic vector matrix and the updated principal component;
performing anti-standardization processing on the three new standard residual error sequence power spectrums to obtain three new anti-standardization residual error sequence power spectrums; wherein the de-normalization process is used to increase the power spectrum values of the adjusted residual sequence from the second scale to a third scale; wherein the second ratio is less than the third ratio;
carrying out inverse Fourier transform on the three new anti-standard residual sequence power spectrums to obtain three new residual sequences; wherein, the three new residual error sequences can represent the errors between the three coordinate time series actual values and the three coordinate time series estimated values.
3. The method of claim 2, wherein before performing the inverse fourier transform on the three new anti-standard residual sequence power spectra, further comprising:
and calculating the phase corresponding to each frequency signal on each residual sequence power spectrum according to each residual sequence power spectrum obtained when the three residual sequences are subjected to Fourier transform.
4. The method of claim 3, wherein said inverse Fourier transforming the three new anti-standard residual sequence power spectra comprises:
and performing inverse Fourier transform on the three new anti-standard residual sequence power spectrums based on the phase corresponding to the frequency signal on each residual sequence power spectrum in the three residual sequences.
5. A coordinate time series processing apparatus, characterized in that the apparatus comprises:
the acquisition module is used for acquiring data from a target survey station and preprocessing the data to obtain three residual error sequences corresponding to coordinate time sequences of the data in three directions; the residual error sequence represents an error between the actual value of the coordinate time sequence and the estimated value of the coordinate time sequence, and the error consists of white noise and colored noise;
the power spectrum calculation module is used for performing Fourier transform on the three residual error sequences to obtain three transformed residual error sequences; calculating power spectrums of the three corresponding residual error sequences according to the three transformed residual error sequences;
the standardization module is used for carrying out standardization processing on the power spectrums of the three residual error sequences to obtain three standard residual error sequence power spectrums; the normalization processing is used for reducing the proportion of the power spectrum value of the residual sequence with the maximum power spectrum value in the sum of the power spectrum values of the three residual sequences from a first proportion to a second proportion to obtain an adjusted residual sequence; wherein the first ratio is greater than the second ratio;
the principal component analysis module is used for solving eigenvalues and eigenvector matrixes of the three standard residual sequence power spectrums and calculating principal components of the three standard residual sequence power spectrums based on the eigenvector matrixes; dividing the principal component into a first principal component and a second principal component based on the eigenvalue, wherein the first principal component is used for characterizing a first part of white noise, the second principal component is used for characterizing mixed noise of the colored noise and a second part of white noise, and the first part of white noise is different from the second part of white noise; deleting the second principal component from the principal component, and reserving the first principal component to obtain an updated principal component; calculating three new residual sequences based on the updated principal component;
and the replacing module is used for replacing the original residual sequence with the three new residual sequences to obtain three new coordinate time sequences.
6. The apparatus of claim 5, further comprising: an anti-normalization module;
the principal component analysis module is also used for calculating three new standard residual sequence power spectrums according to the characteristic vector matrix and the updated principal components;
the anti-standardization module is used for carrying out anti-standardization processing on the three new standard residual error sequence power spectrums to obtain three new anti-standard residual error sequence power spectrums; wherein the de-normalization process is used to increase the power spectrum values of the adjusted residual sequence from the second scale to a third scale; wherein the second ratio is less than the third ratio;
the power spectrum calculation module is further configured to perform inverse fourier transform on the three new anti-standard residual sequence power spectrums to obtain three new residual sequences; wherein, the three new residual error sequences can represent the errors between the three coordinate time series actual values and the three coordinate time series estimated values.
7. The apparatus of claim 6,
the power spectrum calculating module is further configured to calculate a phase corresponding to each frequency signal on each residual sequence power spectrum according to each residual sequence power spectrum obtained when the three residual sequences are subjected to fourier transform.
8. The apparatus of claim 7,
the calculation power spectrum module is further configured to perform inverse fourier transform on the three new anti-standard residual sequence power spectrums based on a phase corresponding to the frequency signal on each of the three residual sequence power spectrums.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911102015.8A CN111142134B (en) | 2019-11-12 | 2019-11-12 | Coordinate time series processing method and device |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911102015.8A CN111142134B (en) | 2019-11-12 | 2019-11-12 | Coordinate time series processing method and device |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111142134A CN111142134A (en) | 2020-05-12 |
| CN111142134B true CN111142134B (en) | 2022-03-11 |
Family
ID=70517045
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911102015.8A Active CN111142134B (en) | 2019-11-12 | 2019-11-12 | Coordinate time series processing method and device |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111142134B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114844796A (en) * | 2022-04-29 | 2022-08-02 | 济南浪潮数据技术有限公司 | Method, device and medium for detecting abnormity of time-series KPI |
| CN117609710B (en) * | 2024-01-24 | 2024-04-12 | 中国电建集团西北勘测设计研究院有限公司 | Method and device for preventing normal jump of monitoring data from being removed |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104392414A (en) * | 2014-11-04 | 2015-03-04 | 河海大学 | Establishing method of regional CORS coordinate time series noise model |
| CN105572703A (en) * | 2015-12-17 | 2016-05-11 | 武汉大学 | GPS time sequence generalized common mode error extraction method |
| CN106814378A (en) * | 2017-01-17 | 2017-06-09 | 华东交通大学 | A kind of GNSS location time series cyclophysis method for digging |
| CN109188466A (en) * | 2018-09-29 | 2019-01-11 | 华东交通大学 | A kind of GNSS base station crust motion velocity field estimation method for taking nonlinear change into account |
| CN109709585A (en) * | 2018-12-04 | 2019-05-03 | 中铁第四勘察设计院集团有限公司 | The method for removing coloured noise in GPS coordinate time series |
-
2019
- 2019-11-12 CN CN201911102015.8A patent/CN111142134B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104392414A (en) * | 2014-11-04 | 2015-03-04 | 河海大学 | Establishing method of regional CORS coordinate time series noise model |
| CN105572703A (en) * | 2015-12-17 | 2016-05-11 | 武汉大学 | GPS time sequence generalized common mode error extraction method |
| CN106814378A (en) * | 2017-01-17 | 2017-06-09 | 华东交通大学 | A kind of GNSS location time series cyclophysis method for digging |
| CN109188466A (en) * | 2018-09-29 | 2019-01-11 | 华东交通大学 | A kind of GNSS base station crust motion velocity field estimation method for taking nonlinear change into account |
| CN109709585A (en) * | 2018-12-04 | 2019-05-03 | 中铁第四勘察设计院集团有限公司 | The method for removing coloured noise in GPS coordinate time series |
Non-Patent Citations (5)
| Title |
|---|
| Extracting the regional common-mode component of GPS station position time series from dense continuous network;Yunfeng Tian et.al.;《JOURNAL OF GEOPHYSICAL RESEARCH-SOLID EARTH》;20160228;第121卷(第2期);第1080—1096页 * |
| GPS坐标时间序列噪声估计及相关性分析;马俊等;《武汉大学学报•信息科学版》;20181031;第43卷(第10期);第1451—1457页 * |
| Spatio-temporal analysis of the land subsidence in the UK using Independent Component Analysis;Bin Liu et.al.;《2014 Third International Workshop on Earth Observation and Remote Sensing Applications (EORSA)》;20141020;第1—5页 * |
| 区域CORS站坐标时间序列特征分析;武曙光;《中国优秀硕士学位论文全文数据库 基础科学辑》;20170815(第8期);正文第7—15页 * |
| 测站坐标时间序列随机模型的解析表达;明锋等;《导航定位学报》;20150930;第3卷(第3期);第59—62页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111142134A (en) | 2020-05-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Pinnegar et al. | The S-transform with windows of arbitrary and varying shape | |
| CN110361778B (en) | Seismic data reconstruction method based on generation countermeasure network | |
| CN106845010B (en) | Low-frequency oscillation dominant mode identification method based on improved SVD noise reduction and Prony | |
| CN107480391B (en) | Near-fault non-stationary seismic oscillation simulation method based on data driving | |
| CN110069868B (en) | GNSS station nonlinear motion modeling method and device | |
| CN106441288A (en) | Adaptive wavelet denoising method for accelerometer | |
| CN113486574B (en) | Sound velocity profile completion method and device based on historical data and machine learning | |
| CN111142134B (en) | Coordinate time series processing method and device | |
| CN109709585B (en) | Method for removing colored noise in GPS coordinate time sequence | |
| CN105353408B (en) | A kind of Wigner higher-order spectrum seismic signal spectral factorization methods based on match tracing | |
| CN109214469B (en) | Multi-source signal separation method based on non-negative tensor decomposition | |
| CN110503060B (en) | Spectral signal denoising method and system | |
| CN113723171B (en) | Electroencephalogram signal denoising method based on residual error generation countermeasure network | |
| CN104143031B (en) | A kind of vegetation index time series data reconstructing method based on Multiscale Wavelet Decomposition | |
| CN112487604A (en) | Long-time nonlinear drift compensation method for output data of marine gravimeter | |
| CN111260131A (en) | Short-term traffic flow prediction method and device | |
| CN115758876A (en) | Method, system and computer equipment for forecasting accuracy of wind speed and wind direction | |
| Lombardi et al. | On-line bayesian estimation of signals in symmetric/spl alpha/-stable noise | |
| Hoffman | The effect of thinning and superobservations in a simple one-dimensional data analysis with mischaracterized error | |
| CN112766127B (en) | Lei Yundian charge positioning method based on complementary set modal decomposition and SG filtering | |
| CN113935246A (en) | Signal robust sparse time-frequency analysis method, terminal equipment and storage medium | |
| Cicone et al. | Jot: a variational signal decomposition into jump, oscillation and trend | |
| CN109725276B (en) | Optical fiber current transformer random error suppression method based on wavelet analysis | |
| CN111175608A (en) | Power distribution network harmonic responsibility quantitative division method based on accelerated independent component analysis | |
| CN113640891B (en) | Singular spectrum analysis-based transient electromagnetic detection data noise filtering method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |