CN111965669B - Separation method for observation pier thermal expansion signals in GNSS time sequence - Google Patents

Separation method for observation pier thermal expansion signals in GNSS time sequence Download PDF

Info

Publication number
CN111965669B
CN111965669B CN202010816037.7A CN202010816037A CN111965669B CN 111965669 B CN111965669 B CN 111965669B CN 202010816037 A CN202010816037 A CN 202010816037A CN 111965669 B CN111965669 B CN 111965669B
Authority
CN
China
Prior art keywords
time sequence
time
gnss
frequency
observation
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
Application number
CN202010816037.7A
Other languages
Chinese (zh)
Other versions
CN111965669A (en
Inventor
王锴华
肖玉钢
郭祚界
张辛
陈华
陈渠森
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Changjiang Spatial Information Technology Engineering Co ltd
Changjiang Institute of Survey Planning Design and Research Co Ltd
Original Assignee
Changjiang Spatial Information Technology Engineering Co ltd
Changjiang Institute of Survey Planning Design and Research Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Changjiang Spatial Information Technology Engineering Co ltd, Changjiang Institute of Survey Planning Design and Research Co Ltd filed Critical Changjiang Spatial Information Technology Engineering Co ltd
Priority to CN202010816037.7A priority Critical patent/CN111965669B/en
Publication of CN111965669A publication Critical patent/CN111965669A/en
Application granted granted Critical
Publication of CN111965669B publication Critical patent/CN111965669B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/21—Interference related issues ; Issues related to cross-correlation, spoofing or other methods of denial of service
    • 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/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/23—Testing, monitoring, correcting or calibrating of receiver elements
    • G01S19/235—Calibration of receiver components
    • 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/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain

Abstract

The invention discloses a method for separating thermal expansion signals of observation piers in a GNSS time sequence. The method comprises the following steps: selecting a comparison station near a GNSS reference station to be researched; step two: calculating to obtain a short baseline solution time sequence; step three: calculating the amplitude and phase of the low-frequency signal in the time sequence, effectively separating the low-frequency thermal expansion signal, and deducting the low-frequency thermal expansion signal from the sequence; step four: decomposing the time sequence into time-varying high-frequency signal components with different periods by adopting an SSAM method; step five: and determining the order of the significant signal component through omega correlation test, combining and reconstructing a time sequence, and effectively separating the thermal expansion signal of the high-frequency observation pier from other low-frequency signals and noise. The method has the advantage that a truer and more accurate observation pier thermal expansion signal can be extracted from the GNSS time sequence, and a more accurate and reliable GNSS reference station three-dimensional velocity field and uncertainty thereof can be calculated based on the true and accurate observation pier thermal expansion signal.

Description

Separation method for observation pier thermal expansion signals in GNSS time sequence
Technical Field
The invention relates to the field of nonlinear time series analysis, in particular to a method for separating thermal expansion signals of observation piers in a GNSS time series.
Background
The coordinate time sequence (time sequence for short) of the GNSS reference station which continuously runs can reflect the precise displacement of the earth surface where the reference station is located and the change of the earth surface along with the time, is one of important products of a GNSS reference station network, and is widely applied to the fields of coordinate reference frame establishment and maintenance, earth surface and sea level change monitoring in the global and regional range and the like. As an important component for fixedly connecting bedrock and a GNSS antenna, the GNSS observation pier can generate regular thermoelastic displacement due to periodic changes of ambient environment or surface temperature of the observation pier, and the regular thermoelastic displacement is represented as periodic oscillation in a time sequence, namely a Thermal Expansion (TEM) signal of the observation pier. The method has the advantages that TEM signals in the time sequence are effectively separated and extracted, the physical mechanism generated by part of nonlinear changes in the time sequence is beneficial to knowing, the accuracy and the reliability of the three-dimensional velocity field of the GNSS reference station for constructing the coordinate reference frame and the uncertainty estimation of the three-dimensional velocity field are improved, and the earth surface displacement information acquired based on the time sequence is more accurate and real.
TEM signals can be divided into low frequency (e.g., anniversary, half-week-year period) and high frequency portions (e.g., sunday, half-week-day period). The current research results show that in the regions with obvious seasonal temperature changes such as the middle hemisphere and the high latitude, the annual amplitude of the low-frequency TEM signal can maximally exceed 3 mm; in regions with large day-night temperature difference, such as inland regions and plateaus, the daily change of the coordinates of the reference station caused by the high-frequency TEM signal can even reach 5 mm. Currently, two problems still exist in extracting a TEM signal existing in a sequence by adopting an existing GNSS precision data processing model and a time sequence analysis method: firstly, nonlinear changes in a time sequence are from a plurality of sources, such as geophysical effects including system errors and environmental loads related to GNSS, local factors related to a survey station environment and the like, so that separation of millimeter-level TEM signals and other signals and noises with larger magnitudes is difficult; secondly, in a common time sequence product of the day solution, high-frequency TEM signals are aliased into false low-frequency signals due to low sampling rate, so that the separation of TEM signals with different significant periods is difficult.
Therefore, it is needed to develop a method for separating thermal expansion signals of observation piers in GNSS time series.
Disclosure of Invention
The invention aims to provide a method for separating thermal expansion signals of observation piers in a GNSS time sequence, which can extract the advantages of truer and more accurate thermal expansion signals of the observation piers from the GNSS time sequence, and can obtain a more accurate and reliable GNSS reference station three-dimensional velocity field and uncertainty estimation thereof based on the thermal expansion signals; the method is beneficial to revealing a nonlinear signal mechanism driven by temperature change in the time sequence, further weakening the influence of unmodeled geophysical effect on the nonlinear signal of the time sequence, acquiring a more accurate and reliable GNSS reference station three-dimensional velocity field and uncertainty estimation value thereof, and providing theoretical support for establishing an earth reference frame taking the nonlinear motion of the reference station into account.
In order to achieve the purpose, the technical scheme of the invention is as follows: a separation method for thermal expansion signals of observation piers in a GNSS time sequence is characterized by comprising the following steps: comprises the following steps of (a) carrying out,
the method comprises the following steps: selecting comparison stations which have consistent observation conditions and are fixed on bedrock within 150 meters near a GNSS reference station to be researched, and respectively erecting GNSS receivers and antennas of the same type on the two reference stations to simultaneously acquire data;
step two: forming a short baseline by two reference stations, selecting every 4 hours as a time interval based on GNSS double-difference observed values, and calculating by adopting a precise data processing model to obtain a baseline solution time sequence;
step three: the method comprises the steps of simultaneously calculating the amplitude and the phase of a low-frequency signal in a time sequence by adopting an MLE (maximum likelihood estimation) method, and effectively separating a thermal expansion signal of a low-frequency observation pier from the time sequence;
step four: selecting a proper window length according to the time sequence span, and decomposing the time sequence obtained in the step (3) into time-varying high-frequency signal components with different significant periods by adopting an SSAM method;
step five: and through omega correlation test, calculating the cross correlation coefficient among all components, determining the order of the components of the significant signal, combining and reconstructing a time sequence according to the order, and effectively separating the thermal expansion signal of the high-frequency observation pier from other low-frequency signals and noises.
In the above technical solution, in the third step, the MLE method specifically includes:
the following functional model is used to describe the position of the reference station at any time in the GNSS coordinate time series:
Figure BDA0002632720930000031
in the formula (1), tiAt the moment i, i is a positive integer; x (t)i) Is tiA reference station position at a time; x (t)0) Is a reference time position; vXIs a linear velocity;
Figure BDA0002632720930000032
respectively at jth frequency fjTotal amplitude and phase of the lower harmonic, generally fjTaking 1cycle/year (cpy) which represents the frequency of a regression year (period is about 365.2 days), wherein j is a positive integer; gkIndicating the amount of step due to various causes; t iskIs the moment when the step change occurs; h is a Havessit step function; before the step change, the value is 0, and after the step change, the value is 1; deltaPSD(ti) The displacement is the same shock and the displacement after the shock; epsilon is a random process; wherein, deltaPSD(ti) The model of (2) is as follows:
Figure BDA0002632720930000033
in the formula (2), nlAnd neRespectively the number of logarithm and exponential functions in the parameter model;
Figure BDA0002632720930000034
and
Figure BDA0002632720930000035
the amplitudes corresponding to the ith logarithm and the exponential function respectively;
Figure BDA0002632720930000036
and
Figure BDA0002632720930000037
respectively is the relaxation time after earthquake;
Figure BDA0002632720930000038
and
Figure BDA0002632720930000039
respectively corresponding to the ith logarithm and exponential function, and t is decimal year time; in the formula (1), the random process ε (t)i) From an amplitude of aw、bκA linear combination of white noise and colored noise description of (a);
writing equation (1) in matrix form:
Y＝Ax+r (3)
in the formula (3), Y is a GNSS coordinate time series observation value column vector, x is a parameter column vector to be estimated, A is a coefficient matrix of the parameter vector to be estimated, and r is a residual error;
wherein the X vector contains the parameter X (t) to be estimated0)、VX、
Figure BDA00026327209300000310
gk、aw、bκ(ii) a Least squares estimation of x-vectors
Figure BDA00026327209300000311
Expressed as:
Figure BDA00026327209300000312
in the formula (4), C is a covariance matrix, and the residual is defined as
Figure BDA00026327209300000313
Determining covariance matrix C by selecting different noise models, and calculating
Figure BDA0002632720930000041
The residual r is estimated and the MLE estimate is calculated as shown in (5):
Figure BDA0002632720930000042
and under the condition that the MLE estimation value is maximum, solving each parameter to be estimated and the uncertainty thereof under the assumption of the optimal noise model.
In the above technical solution, in the fourth step, the SSAM method specifically includes:
for a typical time series x with N observations, the symmetric hysteretic correlation Toeplitz matrix at a selected window M can be expressed as:
Figure BDA0002632720930000043
wherein the element cjIs defined as:
Figure BDA0002632720930000044
in formula (7): c. CjJ in (1) is arbitrarily taken from 0 to M-1;
estimating a matrix
Figure BDA0002632720930000045
Of the k-th order of eigenvalues λkAnd a feature vector EkAnd according to lambdakThe values are sorted in descending order of magnitude, where the k-th order component can be expressed as:
Figure BDA0002632720930000046
after recombining each component, the original time series x can be decomposed into different components xkIt can be expressed as:
Figure BDA0002632720930000047
if N is deleted in the time sequencelAn observed value, then a lag autocorrelation matrix
Figure BDA0002632720930000048
The elements in (1) can be defined as:
Figure BDA0002632720930000049
the corresponding k-th order component of ignoring missing data is:
Figure BDA0002632720930000051
definition f ∈ [0,1 ]]Fractional part specified for maximum missing observation at window size M, if NlIf f.M is smaller than
Figure BDA0002632720930000052
The reconstructed time series component can be expressed by equation (9) for one missing data.
In the above technical solution, in step five, the order of the reconstruction components is determined by calculating an ω correlation coefficient matrix between the reconstruction components, specifically, selecting the sum of the first r components with the largest autocorrelation coefficient and the smallest cross correlation coefficient; wherein, the method is used for evaluating two time sequences x with the length of N1And x2The ω correlation coefficient of inter-separability is defined as formula (12):
Figure BDA0002632720930000053
where ρ isω(x1,x2) For the correlation coefficient, the weight w can be calculated by the following equation (13):
Figure BDA0002632720930000054
the invention has the following advantages:
(1) the invention adopts the GNSS short baseline time sequence to separate the result more accurately: according to the method, the GNSS short baseline 4-hour solution time sequence without the observation pier measuring station is adopted, so that most of GNSS nonlinear signals and error sources irrelevant to the observation pier are eliminated, the obtained residual signals in the time sequence mainly comprise the thermal expansion displacement of the observation pier, the aliasing effect of high-frequency signals is weakened to the greatest extent, and the separation result based on the short baseline time sequence is more accurate;
(2) the low-frequency signal amplitude calculated by adopting the MLE method is more accurate: the invention adopts the most accurate MLE method acknowledged at present, can estimate linear trend term, periodic term, noise type and its characteristic isoparametric and its uncertainty at the same time, separate low-frequency TEM signal in the time sequence from observing the noise effectively, this method is regarded as a kind of unbiased estimation, the result is closer to the true value than the existing method;
(3) after the low-frequency signals in the time sequence are deducted, the high-frequency observation pier thermal expansion signals extracted by the SSAM method are closer to the real signal size: the time sequence is decomposed and reconstructed by a singular spectrum analysis method capable of processing discontinuous data, the time-varying high-frequency TEM signal is effectively separated from other low-frequency signals and noise, and in addition, the frequency mixing effect of the sunday and semisunday signals is avoided due to the sampling interval of 4 hours, so that the extracted high-frequency TEM signal result is more reliable.
The invention mainly aims at the problem of effective extraction of thermal expansion signals of observation piers in the current GNSS time sequence, provides a separation method of the thermal expansion signals of the observation piers in the GNSS time sequence, is beneficial to disclosing a temperature change driven nonlinear signal mechanism in the time sequence, further weakens the influence of unmodeled geophysical effect on the nonlinear signals of the time sequence, calculates a more accurate and reliable GNSS base station three-dimensional velocity field and uncertainty thereof, and provides theoretical support for establishing an earth reference frame taking account of the nonlinear motion of the base station.
Drawings
FIG. 1 is a vertical direction time series solution for a short baseline 4 hour period of GNSS in an embodiment of the present invention.
FIG. 2 shows the results of a baseline JOJO vertical time series decomposition of an embodiment of the present invention (only the first ten components are shown, 2015.5-2016.0).
Fig. 3 is a base-line WAWZ east time series ω correlation coefficient matrix according to an embodiment of the present invention.
FIG. 4 is a HEHEH baseline solution and temperature time series (2014/06/02-2014/07/04) after reconstruction in accordance with an embodiment of the present invention.
FIG. 5 is a time series difference between the short baselines Y2Y3, Y2YR (left panel) and WAWR, WAWZ (right panel) in accordance with an embodiment of the present invention.
Fig. 6 is a graph of the horizontal displacement profiles of baseline HEHE (upper panel, year 2015, 6 months to 12 months) and ZMZ2 (lower panel, year 2015, 6 months to 12 months) for different time periods in accordance with an embodiment of the present invention.
FIG. 7 is a process flow diagram of the present invention.
In fig. 1, the abscissa is time, the ordinate is vertical direction displacement (in mm), the black dots are time period solution results, the gray is uncertainty, and the short baselines are arranged in an increasing manner from top to bottom according to the difference in height of the observation piers.
In fig. 2, the abscissa is time and the ordinate is displacement (in meters), and the first ten components of the decomposed time series are arranged from top to bottom.
In fig. 3, the horizontal and vertical coordinates are the orders, the color depth corresponding to each square is the correlation coefficient, and the lighter the color, the larger the correlation number.
In fig. 4, the abscissa is time, and the ordinate is north direction, east direction, vertical direction displacement and ambient temperature in order from top to bottom; as can be seen in fig. 4, the reconstructed time series contains significant weekday and half-weekday cycle characteristics.
In fig. 5, the abscissa is time, and the ordinate is displacement in the north direction, the east direction, and the vertical direction in order from top to bottom; as can be seen in fig. 5, differences in the height of the observation mound of one meter height cause diurnal variations of 1mm in maximum amplitude, especially in summer where diurnal temperature changes greatly.
In fig. 6, the abscissa is the east-west displacement (east is positive), the ordinate is the north-south displacement (north is positive), solutions of 6 different time periods in a day are respectively represented by dots, triangles, diamonds, squares, pentagons and hexagons, and the lower right corner is a schematic diagram of the short-baseline observation pier.
Detailed Description
The embodiments of the present invention will be described in detail with reference to the accompanying drawings, which are not intended to limit the present invention, but are merely exemplary. While the advantages of the invention will be clear and readily understood by the description.
The principle of the invention is as follows: aiming at the problem of extraction of thermal expansion signals of observation piers in a GNSS reference station time sequence, the invention provides that a non-observation pier reference station is erected near the existing survey station to form a short base line, the short base line is obtained by adopting 4-hour observation time period, a short base line time sequence is obtained by calculation, and the influence of other nonlinear signal sources and high-frequency TEM signal aliasing effects on the time sequence is weakened or even eliminated to the maximum extent; calculating the signal amplitude and phase of the anniversary and semianniversary in the time sequence by adopting a Maximum Likelihood Estimation (MLE) method, and effectively extracting a low-frequency TEM signal in the time sequence; and decomposing and reconstructing the sequence according to a repetition period by adopting a Singular Spectrum Analysis (SSAM) method capable of processing discontinuous data, and effectively separating the high-frequency TEM signal from noise.
With reference to the accompanying drawings: a method for separating thermal expansion signals of observation piers in a GNSS time sequence comprises the following steps,
the method comprises the following steps: selecting a comparison station (without an observation pier) which has basically consistent observation conditions and is fixed on bedrock within 150 meters near a GNSS reference station (with an observation pier with a certain height) to be researched, and respectively erecting a GNSS receiver and antenna acquisition data of the same model at the two reference stations (namely the GNSS reference station to be researched and the GNSS receiver and antenna acquisition data erected at the comparison station are the same in model);
step two: two reference stations (namely a GNSS reference station to be researched and a contrast station fixed on bedrock) form a group of GNSS short baselines with obvious difference of observation piers, the resolving precision and the time resolution are comprehensively considered, and every 4 hours are selected to be a time period (the time period is too short, the whole cycle ambiguity is difficult to fix, and the resolution is poor; the time period is too long, such as 8 hours, a time period and 3 solutions in a day, the resolution of a time sequence is not enough to show the high-frequency signal characteristic to be researched, such as a 12-hour half-week day signal), and based on the GNSS double-difference observation value in each time period, universal precision data processing software and a model are adopted for calculating to obtain a short-baseline solution time sequence;
the time sequence obtained by the step can weaken or even eliminate most of system errors caused by the GNSS technology, such as receiver and antenna difference influence, high-order ionosphere errors, troposphere modeling errors and the like; meanwhile, the nonlinear displacement of the measuring station caused by the large-scale geophysical effect can also be counteracted in a short-baseline solution, such as environmental load and the like; in addition, the time resolution of one solution every 4 hours can effectively avoid the mixing effect of the TEM signals of the weekdays and the semiweekdays in the time sequence on the premise of ensuring the fixed success rate of the whole-cycle ambiguity and the resolving precision;
step three: simultaneously calculating the amplitude and phase of the low-frequency signal (such as anniversary and semianniversary) in the short-baseline solution time sequence obtained in the step two, the type of the random model and corresponding parameter information by adopting an MLE (Multi-level iterative optimization) method, and effectively separating the low-frequency signal (namely the low-frequency observation pier thermal expansion signal) from the short-baseline solution time sequence obtained in the step two; on the premise of assuming the random model type, the method can simultaneously estimate the characteristics of seasonal signals and random noise in the function model, select an optimal noise model according to a maximum likelihood criterion, extract low-frequency TEM signals in a time sequence, and estimate annual and semi-annual amplitude and uncertainty thereof more truly and accurately;
the time sequence comprises low-frequency signals, high-frequency signals and noise, and the separation of the method comprises the separation of the low-frequency signals and the high-frequency signals; the low-frequency signal can be described by trigonometric function harmonic waves (see formula 1), and the magnitude of the signal is determined after the amplitude and the phase are calculated by an MLE method; because the high-frequency signal may have a frequency mixing phenomenon, the method deducts the low-frequency signal from the original sequence, and then extracts the real high-frequency signal from the residual signal and noise by adopting an SSAM method;
step four: selecting a proper window length (the window size is related to the specifically adopted time sequence length and time resolution) according to the time sequence span, and decomposing the time sequence obtained in the step 3 into time-varying high-frequency signal components with different significant periods by adopting an SSAM method;
step five: through an omega correlation test, the cross correlation coefficient among all components is calculated, the high-frequency signal component order is determined, and the high-frequency signal component order is combined and reconstructed to reconstruct a time sequence, such as a cycle of a day of week and a cycle of a half day, so that the high-frequency observation pier thermal expansion signal is effectively separated from other low-frequency signals and noise (shown in figure 7).
Further, in step three, the MLE method is specifically as follows:
the MLE method can simultaneously estimate parameters such as linear trend term, periodic term, noise type and characteristics thereof and uncertainty information of the parameters, is the most accurate noise analysis method which is generally recognized at present, and has the central idea that: for the selected colored noise model or noise combination, selecting a group of undetermined parameters according to a maximum likelihood estimation criterion to enable the joint probability density of the coordinate time sequence residual error and the covariance thereof to be maximum; in general, the following functional model is used to describe the position of a reference station at any time in the GNSS coordinate time series:
Figure BDA0002632720930000091
in the formula (1), tiAt the moment i, i is a positive integer; x (t)i) Is tiA reference station position at a time; x (t)0) Is a reference time position; vXIs a linear velocity;
Figure BDA0002632720930000092
respectively at jth frequency fjTotal amplitude and phase of the lower harmonic, generally fjTaking 1cycle/year (cpy) which represents the frequency of a regression year (period is about 365.2 days), wherein j is a positive integer; gkIndicating the amount of step due to various causes; t iskIs the moment when the step change occurs; h is a Havessit step function; before the step change, the value is 0, and after the step change, the value is 1; deltaPSD(ti) The displacement is the same shock and the displacement after the shock; epsilon is a random process; n is a positive integer; n isgIs a positive integer;
wherein, deltaPSD(ti) The model description that can be given with reference to ITRF 2014:
Figure BDA0002632720930000101
in the formula (2), nlAnd neRespectively the number of logarithm and exponential functions in the parameter model;
Figure BDA0002632720930000102
and
Figure BDA0002632720930000103
the amplitudes corresponding to the ith logarithm and the exponential function respectively;
Figure BDA0002632720930000104
and
Figure BDA0002632720930000105
respectively is the relaxation time after earthquake;
Figure BDA0002632720930000106
and
Figure BDA0002632720930000107
respectively corresponding to the ith logarithm and exponential function; t is the decimal time of year; in the formula (1), the random process ε (t)i) From an amplitude of aw、bκA linear combination of white noise and colored noise description of (a);
writing equation (1) in matrix form:
Y＝Ax+r (3)
in the formula (3), Y is a GNSS coordinate time series observation value column vector; x is a parameter column vector to be estimated; a is a coefficient matrix of a parameter vector to be estimated; r is a residual error;
wherein the X vector includes X (t)0)、VX、
Figure BDA0002632720930000108
gk、aw、bκWaiting for estimating parameters; least squares estimation of x-vectors
Figure BDA0002632720930000109
The representation can be expressed as:
Figure BDA00026327209300001010
in the formula (4), the reaction mixture is,
Figure BDA00026327209300001011
a least squares estimate representing the x vector; t represents matrix transposition; y is a GNSS coordinate time series observation value column vector; a is a coefficient matrix of a parameter vector to be estimated;
where C is the covariance matrix and the residual is defined as
Figure BDA00026327209300001012
Determining covariance matrix C by selecting different noise models, and calculating
Figure BDA00026327209300001013
The residual r is estimated and the MLE estimate is calculated as shown in (5):
Figure BDA00026327209300001014
in the formula (5), the reaction mixture is,
Figure BDA00026327209300001015
a least squares estimate representing the x vector; c is a covariance matrix; r is a residual error; t represents matrix transposition; n represents the vector dimension of the observation value column of the GNSS coordinate time sequence, namely the length of the time sequence;
under the condition that the MLE estimated value is maximum, solving parameters to be estimated and uncertainty thereof under the assumption of the optimal noise model; the MLE method can process time sequences with non-uniform sampling and more missing data, and as an unbiased estimation, the MLE can simultaneously estimate parameters such as linear speed, periodic terms, step and the like and variables such as uncertainty, noise amplitude and the like, the advantages are obvious, the result is more stable and reliable, and especially the data with longer time span is processed.
Further, in step four, the SSAM method is specifically as follows:
for a typical time series x with N observations, the symmetric hysteretic correlation Toeplitz matrix at a selected window M can be expressed as:
Figure BDA0002632720930000111
wherein the element cjIs defined as:
Figure BDA0002632720930000112
in formula (7): c. CjJ in (2) is arbitrarily selected from 0 to M-1; when j is 0, cjCorresponds to c in the formula (6)0(ii) a M represents a selected window length; 1 ≦ M ≦ N; n represents the vector dimension of the observation value column of the GNSS coordinate time sequence, namely the length of the time sequence;
in the formula (6), c1I.e. cjJ in (1); i. k is a variable;
estimating a matrix
Figure BDA0002632720930000113
Of the k-th order of eigenvalues λkAnd a feature vector EkAnd according to lambdakThe values are sorted in descending order of magnitude, where the k-th order component can be expressed as:
Figure BDA0002632720930000114
in the formula (8), the reaction mixture is,
Figure BDA0002632720930000115
representation estimation matrix
Figure BDA0002632720930000116
The kth order component of (1); mIndicating a selected window length; 1 ≦ M ≦ N; n represents the vector dimension of the observation value column of the GNSS coordinate time sequence, namely the length of the time sequence;
Figure BDA0002632720930000117
is distinguished from EkAnother feature vector of (a); x is the number ofi+jA time series representing the time instant i + j;
after recombining each component, the original time series x can be decomposed into different components xkIt can be expressed as:
Figure BDA0002632720930000118
in the formula (9), the reaction mixture is,
Figure BDA0002632720930000119
time series components representing the original time series x decomposition; m represents a selected window length; 1 ≦ M ≦ N; n represents the vector dimension of the observation value column of the GNSS coordinate time sequence, namely the length of the time sequence;
Figure BDA0002632720930000121
representation differentiation from estimation matrix
Figure BDA0002632720930000122
Of the kth order component
Figure BDA0002632720930000123
Another k-th order component of (a);
Figure BDA0002632720930000124
is distinguished from EkAnother feature vector of (a);
if N is deleted in the time sequencelAn observed value, then a lag autocorrelation matrix
Figure BDA0002632720930000125
The elements in (1) can be defined as:
Figure BDA0002632720930000126
in the formula (10), the compound represented by the formula (10),
Figure BDA0002632720930000127
to be distinguished from cjAnother element of (1); n is a radical oflIs another time series length for distinguishing from N; n represents the vector dimension of the observation value column of the GNSS coordinate time sequence, namely the length of the time sequence;
Figure BDA0002632720930000128
a time series representing time l;
Figure BDA0002632720930000129
a time series representing time instants l + j; m represents a selected window length;
corresponding ignoring missing data (i.e. missing N)lObserved value of) is:
Figure BDA00026327209300001210
in formula (11), M represents a selected window length; 1 ≦ M ≦ N; n represents the vector dimension of the observation value column of the GNSS coordinate time sequence, namely the length of the time sequence; n is a radical oflIs another time series length for distinguishing from N;
Figure BDA00026327209300001211
distinguished from the estimation matrix
Figure BDA00026327209300001212
Of the kth order component
Figure BDA00026327209300001213
Another k-th order component of (a);
Figure BDA00026327209300001214
a time series representing the time of i + l;
Figure BDA00026327209300001215
is distinguished from EkAnother feature vector of (a);
definition f ∈ [0,1 ]]Fractional part specified for maximum missing observation at window size M, if NlIf f.M is smaller than
Figure BDA00026327209300001216
The reconstructed time series component can be expressed by equation (9) for one missing data.
Furthermore, in step five, after the original time series is decomposed into components with different significant periodic signals, the first several reconstructed components often contain most of the effective signals in the time series; the order of the reconstruction components can be determined by calculating an omega correlation coefficient matrix among the reconstruction components, specifically, the sum of the first r components with the maximum autocorrelation coefficient and the minimum cross correlation coefficient is selected as much as possible; wherein, the method can be used for evaluating two time sequences x with the length of N1And x2The ω correlation coefficient of inter-separability is defined as formula (12):
Figure BDA0002632720930000131
wherein in the formula (12), ρω(x1,x2) Is a correlation coefficient; w is used to distinguish labels; i is variable and is a positive integer; n represents the vector dimension of the observation value column of the GNSS coordinate time sequence, namely the length of the time sequence;
the weight w can be calculated by the following equation (13):
Figure BDA0002632720930000132
in the formula (13), i is more than 0 and less than M*；
M represents a selected window length;
M*is another selected window length to distinguish from M;
n represents the GNSS coordinate time series observation column vector dimension, i.e., the time series length.
Examples
In order to verify the effectiveness of the separation method of the thermal expansion signals of the observation piers in the GNSS time sequence, the method is applied to 9 groups of GNSS short baselines (the length of the baselines is not more than 150 meters) with obvious differences of the observation piers, and 2 groups of zero baselines and 1 group of short baselines without any difference of the observation piers are selected as a contrast. The technical scheme is explained, and the specific steps are as follows:
step 1, 9 groups of IGS reference stations with obvious observation pier differences and distances smaller than 150 meters are selected to form a GNSS short base line as an experimental group, 2 groups of zero base lines and 1 group of short base lines without any observation pier differences are selected as a comparison group, observation data of each reference station at a sampling rate of 30 seconds are obtained, the time span is from 1 month and 1 day of 2012 to 12 months and 31 days of 2016, and detailed information of each measurement station is shown in table 1. Wherein, the control group is marked as group A, and the experimental group is divided into B, C, D groups according to the difference (including height, material, structure difference and the like) of each short-baseline station to the observation pier. Meanwhile, the environmental temperature data recorded in the meteorological files of each station are obtained, and the detailed information of the height, material, structure and the like of the observation pier is determined, as shown in table 2.
TABLE 1 short Baseline reference station information for each experimental and control group
Figure BDA0002632720930000141
Figure BDA0002632720930000151
Table 2 detailed information of observation piers for each short baseline in control group a and experimental group B, C, D
Figure BDA0002632720930000152
Figure BDA0002632720930000161
Step 2, selecting a time interval of every 4 hours, and calculating a baseline solution time sequence by adopting GNSS precision data processing software GAMIT; the solving strategy and the model comprise: estimating the ambiguity of the whole cycle by using a final precise satellite orbit provided by IGS (geostationary orbit) without estimating troposphere delay parameters and adopting L1 and L2 observation values, wherein the cut-off angle of the antenna height is 15 degrees, and a 4-hour time interval solution is obtained by Kalman filtering; the time series preprocessing comprises the following steps: 1) removing abnormal values, and only preserving results of the sequence median values of +/-0.01 m and +/-0.015 m and uncertainty of <0.1m in the horizontal and vertical direction time sequences; 2) occasional errors, defined as results exceeding 4 times the magnitude of the root mean square error, are removed, where the root mean square error is replaced by a standard deviation approximation. In order not to affect the time-varying noise characteristics of the time series itself, filtering methods such as moving average are not used. For the obvious jump caused by the replacement of part of instruments, determining the jump size by checking a log file of the measuring station, and correcting a time sequence to weaken the influence of the jump size on seasonal signals and noise characteristics; small-magnitude jumps, for which the cause is difficult to determine, are not dealt with. The short baseline vertical time series are shown in fig. 1.
And 3, estimating noise characteristics under the assumption of different random models by adopting Hector software, determining an optimal noise model in each direction of the short baseline based on an MLE (maximum likelihood estimation) criterion, and calculating parameters to be estimated and uncertainty thereof in the formula (1), including amplitude and phase of low-frequency signals (such as anniversaries and semiweeks). Meanwhile, according to information such as temperature data and observation pier height, the thermal expansion model is used for calculating TEM displacement, annual and semi-annual amplitude and phase, the result calculated by the model is compared with the actual observation result and verified, the extraction effect of the low-frequency TEM signal is evaluated, and part of the comparison result is shown in Table 3.
TABLE 3 seasonal item parameter comparison of short baseline TEM model and GNSS observations for a set of 3D
Figure BDA0002632720930000171
And 4, after the low-frequency signal is deducted from the time sequence, decomposing the time sequence obtained in the step 3 into different time-varying signal components according to a significant period by adopting an SSAM method, wherein a partial baseline decomposition result is shown in figure 2.
And 5, taking the window size as 540 in order to decompose details in the original sequence as much as possible without influencing the calculation efficiency. By calculating the ω correlation coefficient matrix and the correlation test between the reconstructed components, the order of the maximum autocorrelation coefficient and the minimum cross correlation coefficient of each time series component is calculated to be 20, and a part of the baseline correlation matrix is shown in fig. 3. The sunday and semiweek day signals of the first 20 components were recombined to obtain a reconstructed time series including high frequency TEM signals of sunday, semiday, etc., and partial results are shown in fig. 4. Meanwhile, according to the difference of the structures of the short-baseline observation pillars in the horizontal direction, the difference between the extracted high-frequency TEM signal and the actual situation of the observation pillars in the vertical and horizontal directions is verified, and the difference is analyzed as shown in FIGS. 5 and 6.
The results of the above embodiment show that the method of the present invention separates low-frequency and high-frequency signals in the GNSS time series caused by the thermal expansion effect of the observation pier. The annual amplitude of the TEM signal in the time sequence calculated by adopting the MLE method can reach 1.86 +/-0.17 mm at most, the variation trend is basically consistent with the ambient temperature of the measuring station, and the random process of the time sequence with about 60 percent of the baseline direction can be described by band-pass noise or first-order Gaussian Markov noise. The comparison result of the TEM model considering the height of the auxiliary observation pier and the GNSS observation value shows that the coincidence degree of the TEM model and the GNSS observation value in the vertical direction is about 82%, which is obviously improved compared with the prior method, and the method provided by the invention is proved to be capable of effectively extracting the low-frequency TEM signal in the time sequence. In addition, high-frequency signals such as sunday, semisunday and the like in the time series are separated by using the SSAM method, and the results show that: in the vertical direction, the extracted signals are mainly consistent with the temperature change trend and have a linear positive correlation with the daily temperature change and the height difference of the observation pier; in the horizontal direction, the daily change of the extracted signal can reach 5mm, the change rule is the same as the diurnal periodic change of the surface temperature of the observation pier caused by direct sunlight, namely, the distribution in different time periods in one day has obvious directivity, the specific direction is consistent with the horizontal direction of the baseline observation pier, and particularly for a reference station with the observation pier in a horizontal direction asymmetric structure, the effectiveness of the method for separating the high-frequency TEM signal in the time sequence is proved.
The specific embodiments described herein are merely illustrative of the spirit of the invention. Various modifications or additions may be made to the described embodiments by those skilled in the art to other GNSS reference stations in other areas, or substituted in a similar manner, without departing from the spirit of the invention or exceeding the scope thereof as defined in the appended claims.

Claims (4)

1. A separation method for thermal expansion signals of observation piers in a GNSS time sequence is characterized by comprising the following steps: comprises the following steps of (a) carrying out,
the method comprises the following steps: selecting comparison stations which have consistent observation conditions and are fixed on bedrock within 150 meters near a GNSS reference station to be researched, and respectively erecting GNSS receivers and antennas of the same type on the two reference stations to acquire data;
step two: forming a short baseline by two reference stations, selecting every 4 hours as a time interval based on GNSS double-difference observed values, and calculating by adopting a precise data processing model to obtain a baseline solution time sequence;
step three: the method comprises the steps of simultaneously calculating the amplitude and the phase of a low-frequency signal in a time sequence by adopting an MLE (maximum likelihood estimation) method, and effectively separating a thermal expansion signal of a low-frequency observation pier from the time sequence;
step four: selecting a proper window length according to the time sequence span, and decomposing the time sequence obtained in the step three into time-varying high-frequency signal components with different significant periods by adopting an SSAM method;
step five: and through omega correlation test, calculating the cross correlation coefficient among all components, determining the order of the components of the significant signal, combining and reconstructing a time sequence according to the order, and effectively separating the thermal expansion signal of the high-frequency observation pier from other low-frequency signals and noises.
2. The method of claim 1, wherein the method comprises the following steps: in step three, the MLE method is specifically as follows:
the following functional model is used to describe the position of the reference station at any time in the GNSS coordinate time series:
Figure FDA0003137658500000011
in the formula (1), tiAt the moment i, i is a positive integer; x (t)i) Is tiA reference station position at a time; x (t)0) Is a reference time position; vXIs a linear velocity;
Figure FDA0003137658500000012
respectively at jth frequency fjTotal amplitude and phase of the lower harmonic, fj1cycle/year is taken to represent the frequency of the regression year, and j is a positive integer; gkIndicating the amount of step due to various causes; t iskIs the moment when the step change occurs; h is a Havessit step function; before the step change, the value is 0, and after the step change, the value is 1; deltaPSD(ti) The displacement is the same shock and the displacement after the shock; epsilon is a random process; wherein, deltaPSD(ti) The model of (2) is as follows:
Figure FDA0003137658500000021
in the formula (2), nlAnd neRespectively the number of logarithm and exponential functions in the parameter model;
Figure FDA0003137658500000022
and
Figure FDA0003137658500000023
the amplitudes corresponding to the ith logarithm and the exponential function respectively;
Figure FDA0003137658500000024
and
Figure FDA0003137658500000025
respectively is the relaxation time after earthquake;
Figure FDA0003137658500000026
and
Figure FDA0003137658500000027
respectively corresponding to the ith logarithm and exponential function, and t is decimal year time; in the formula (1), the random process ε (t)i) From an amplitude of aw、bκA linear combination of white noise and colored noise description of (a);
writing equation (1) in matrix form:
Y＝Ax+r (3)
in the formula (3), Y is a GNSS coordinate time series observation value column vector, x is a parameter column vector to be estimated, A is a coefficient matrix of the parameter vector to be estimated, and r is a residual error; wherein the X vector contains the parameter X (t) to be estimated0)、VX、
Figure FDA0003137658500000028
gk、aw、bκ(ii) a Least squares estimation of x-vectors
Figure FDA0003137658500000029
Expressed as:
Figure FDA00031376585000000210
in the formula (4), C is a covariance matrix, and the residual is defined as
Figure FDA00031376585000000211
Determining covariance matrix C by selecting different noise models, and calculating
Figure FDA00031376585000000212
The residual r is estimated and the MLE estimate is calculated as shown in (5):
Figure FDA00031376585000000213
in the formula (5), N represents a vector dimension of a GNSS coordinate time series observation value column, i.e., a time series length;
and under the condition that the MLE estimated value is maximum, selecting an optimal noise model according to a maximum likelihood criterion, and solving each parameter to be estimated and the uncertainty thereof under the assumption of the optimal noise model.
3. The method of claim 2, wherein the method comprises the steps of: in step four, the SSAM method is specifically as follows:
for a typical time series x with N observations, the symmetric hysteretic correlation Toeplitz matrix at a selected window M can be expressed as:
Figure FDA0003137658500000031
wherein the element cjIs defined as:
Figure FDA0003137658500000032
in formula (7): c. CjJ in (1) is arbitrarily taken from 0 to M-1; estimating a matrix
Figure FDA0003137658500000033
Of the k-th order of eigenvalues λkAnd a feature vector EkAnd according to lambdakThe values are sorted in descending order of magnitude, where the k-th order component can be expressed as:
Figure FDA0003137658500000034
after recombining each component, the original time series x can be decomposed into different components xkIt can be expressed as:
Figure FDA0003137658500000035
if N is deleted in the time sequencelAn observed value, then a lag autocorrelation matrix
Figure FDA0003137658500000036
The elements in (1) can be defined as:
Figure FDA0003137658500000037
the corresponding k-th order component of ignoring missing data is:
Figure FDA0003137658500000038
definition f ∈ [0,1 ]]Fractional part specified for maximum missing observation at window size M, if NlIf f.M is smaller than
Figure FDA0003137658500000039
The reconstructed time series component can be expressed by equation (9) for one missing data.
4. The method of claim 3, wherein the method comprises the following steps: in the fifth step, the order of the reconstruction components is determined by calculating an ω correlation coefficient matrix between the reconstruction components, specifically, selecting the sum of the first r components with the largest autocorrelation coefficient and the smallest cross correlation coefficient; therein is used forEvaluation of two time series x of length N1And x2The ω correlation coefficient of inter-separability is defined as formula (12):
Figure FDA0003137658500000041
where ρ isω(x1，x2) For the correlation coefficient, the weight w can be calculated by the following equation (13):
Figure FDA0003137658500000042
CN202010816037.7A 2020-08-14 2020-08-14 Separation method for observation pier thermal expansion signals in GNSS time sequence Active CN111965669B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010816037.7A CN111965669B (en) 2020-08-14 2020-08-14 Separation method for observation pier thermal expansion signals in GNSS time sequence

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010816037.7A CN111965669B (en) 2020-08-14 2020-08-14 Separation method for observation pier thermal expansion signals in GNSS time sequence

Publications (2)

Publication Number Publication Date
CN111965669A CN111965669A (en) 2020-11-20
CN111965669B true CN111965669B (en) 2021-09-03

Family

ID=73365582

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010816037.7A Active CN111965669B (en) 2020-08-14 2020-08-14 Separation method for observation pier thermal expansion signals in GNSS time sequence

Country Status (1)

Country Link
CN (1) CN111965669B (en)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106597484A (en) * 2016-12-12 2017-04-26 武汉大学 Method for accurately quantifying influence of thermal expansion effect on GPS coordinate time series
KR20180013416A (en) * 2016-07-29 2018-02-07 한국해양과학기술원 Interference and signal quality performance analysis method of reference station selection for GNSS based correction
CN110398753A (en) * 2019-06-28 2019-11-01 武汉大学 GNSS survey station coordinate time sequence periodicity detection method and system

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102498415B (en) * 2009-09-19 2014-04-16 天宝导航有限公司 Gnss signal processing with rover ambiguity fixing

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20180013416A (en) * 2016-07-29 2018-02-07 한국해양과학기술원 Interference and signal quality performance analysis method of reference station selection for GNSS based correction
CN106597484A (en) * 2016-12-12 2017-04-26 武汉大学 Method for accurately quantifying influence of thermal expansion effect on GPS coordinate time series
CN110398753A (en) * 2019-06-28 2019-11-01 武汉大学 GNSS survey station coordinate time sequence periodicity detection method and system

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
GNSS坐标时间序列分析理论与方法及展望;姜卫平等;《武汉大学学报·信息科学版》;20181231;第43卷(第12期);2112-2123 *

Also Published As

Publication number Publication date
CN111965669A (en) 2020-11-20

Similar Documents

Publication Publication Date Title
Nilsson et al. Improved retrieval of land ice topography from CryoSat-2 data and its impact for volume-change estimation of the Greenland Ice Sheet
Williams et al. Error analysis of continuous GPS position time series
Reuter et al. A joint effort to deliver satellite retrieved atmospheric CO 2 concentrations for surface flux inversions: the ensemble median algorithm EMMA
CN110058236A (en) It is a kind of towards three-dimensional Ground Deformation estimation InSAR and GNSS determine Quan Fangfa
Zhang et al. Surface ice flow velocity and tide retrieval of the Amery ice shelf using precise point positioning
Van Dam et al. Predictions of crustal deformation and of geoid and sea-level variability caused by oceanic and atmospheric loading
Gong et al. Temporal filtering of InSAR data using statistical parameters from NWP models
CN110909449B (en) Multi-source data ionization layer region reporting method
He et al. Noise analysis for environmental loading effect on GPS position time series
Salstein et al. Uncertainties in atmospheric surface pressure fields from global analyses
Elyouncha et al. Joint retrieval of ocean surface wind and current vectors from satellite SAR data using a Bayesian inversion method
Xiang et al. Characterizing the seasonal hydrological loading over the asian continent using GPS, GRACE, and hydrological model
CN111965669B (en) Separation method for observation pier thermal expansion signals in GNSS time sequence
Yan et al. Correction of atmospheric delay error of airborne and spaceborne GNSS-R sea surface altimetry
Geiger et al. Coincident buoy‐and SAR‐derived surface fluxes in the western Weddell Sea during Ice Station Weddell 1992
Jakobsen et al. Analysis of local ionospheric time varying characteristics with singular value decomposition
Lindsay Ice deformation near SHEBA
Oktar et al. Research of behaviours of continuous GNSS stations by signal
CN111965670B (en) Method for quantifying aliasing of thermal expansion signals of GNSS time sequence high-frequency observation piers
Reid et al. A-chaim: Near-real-time data assimilation of the high latitude ionosphere with a particle filter
Olivieri et al. Spatial sea-level reconstruction in the Baltic Sea and in the pacific Ocean from tide gauges observations
Sherwood Climate signal mapping and an application to atmospheric tides
Hwang et al. Shallow‐water gravity anomalies from satellite altimetry: Case studies in the east china sea and Taiwan strait
Seaman et al. The impact of manually derived Southern Hemisphere sea level pressure data upon forecasts from a global model
Blanc et al. High‐resolution, high‐accuracy altimeter derived mean sea surface in the Norwegian Sea

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