CN115373005A - High-precision product conversion method between satellite navigation signals - Google Patents

High-precision product conversion method between satellite navigation signals Download PDF

Info

Publication number
CN115373005A
CN115373005A CN202211049892.5A CN202211049892A CN115373005A CN 115373005 A CN115373005 A CN 115373005A CN 202211049892 A CN202211049892 A CN 202211049892A CN 115373005 A CN115373005 A CN 115373005A
Authority
CN
China
Prior art keywords
satellite
orbit
parameters
parameter
ambiguity
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.)
Pending
Application number
CN202211049892.5A
Other languages
Chinese (zh)
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.)
63921 Troops of PLA
Original Assignee
63921 Troops of PLA
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 63921 Troops of PLA filed Critical 63921 Troops of PLA
Priority to CN202211049892.5A priority Critical patent/CN115373005A/en
Publication of CN115373005A publication Critical patent/CN115373005A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39Determining 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
    • G01S19/42Determining position
    • G01S19/43Determining position using carrier phase measurements, e.g. kinematic positioning; using long or short baseline interferometry

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention relates to a high-precision product conversion method among satellite navigation signals, which comprises the following steps: performing orbit determination and real-time filtering resolving on the GNSS satellite by utilizing observation data of partial frequency points of the GNSS satellite by a wide area monitoring station to obtain a clock error parameter and resolving an inter-code deviation parameter and a phase deviation parameter; the method comprises the steps that a wide area monitoring station is utilized to conduct precise point positioning processing on observation data of partial frequency points of a GNSS satellite, and an ionosphere delay parameter and a troposphere delay parameter are calculated by a real-time filtering method; resolving pseudo range intra-frequency differences and phase intra-frequency differences of partial frequency points and all other frequency points by using observation data of all frequency points of the GNSS satellite by using a local area monitoring station; forecasting the pseudo-range frequency inner difference and the phase frequency inner difference to generate code deviation parameters and phase deviation parameters of all frequency points of the full arc section of the GNSS satellite; broadcasting satellite orbit parameters, clock error parameters, delay parameters and deviation parameters, and performing positioning calculation. The invention can provide high-precision space-time service for GNSS users with different frequency points all over the world.

Description

High-precision product conversion method between satellite navigation signals
Technical Field
The invention relates to the technical field of satellite navigation positioning, in particular to a high-precision product conversion method between satellite navigation signals.
Background
The satellite navigation system is an important space infrastructure which can provide meter-level positioning navigation time service for users. With social progress and technological development, the demand of people for navigation and location services is increased explosively, and a satellite navigation system is difficult to meet the requirements of users on real-time, accurate, ubiquitous, reliable and other service performances. The satellite navigation enhancement is to enhance and improve the performance of a satellite navigation system by technologies such as real-time monitoring of satellite navigation signals, error source correction and the like on the basis of the service of the satellite navigation system so as to improve the navigation and positioning performance of a user. Global uniformly distributed wide-area ground monitoring stations are established by scientific organizations at home and abroad, such as a Global Navigation Satellite System (GNSS) service organization and an international GNSS monitoring and evaluating organization, so that monitoring of GNSS Satellite downlink Navigation signals is realized, and original observation data, a data processing strategy and a high-precision Navigation Satellite product are freely disclosed to the scientific community and the public, so that high-precision modeling of parameters such as Satellite orbit parameters, clock error parameters, inter-code deviation and phase floating point deviation, ionospheric delay, tropospheric delay and the like is possible. However, it is often difficult for these ground stations to monitor GNSS downlink signals to cover all of the GNSS downlink navigation signals, and users cannot use high-precision products to achieve high-precision navigation positioning when using the ground stations because it is difficult to provide signals of high-precision products without monitoring.
Disclosure of Invention
In order to solve the problems in the prior art, an object of the present invention is to provide a method for converting high-precision products between satellite navigation signals, which can realize the conversion of high-precision products between GNSS signals and provide high-precision space-time services to GNSS users at different frequency points around the world, in consideration of the difference relationship between signals received by different ground stations.
In order to achieve the above object, the present invention provides a method for converting high precision products between satellite navigation signals, comprising:
s1, performing orbit determination on the GNSS satellite by utilizing observation data of partial frequency points of the GNSS satellite by using a wide area monitoring station;
s2, utilizing a wide area monitoring station to carry out real-time filtering resolving on observation data of partial frequency points of the GNSS satellite to obtain clock error parameters;
s3, decoding observation data of partial frequency points of the GNSS satellite by using a wide area monitoring station to calculate an inter-code deviation parameter and a phase deviation parameter;
s4, performing precise single-point positioning processing on observation data of partial frequency points of the GNSS satellite by using a wide area monitoring station, and solving ionospheric delay parameters and tropospheric delay parameters by using a real-time filtering method;
s5, resolving pseudo range frequency internal differences and phase frequency internal differences of partial frequency points and all other frequency points by using a local monitoring station to observation data of all frequency points of the GNSS satellite;
s6, forecasting the pseudo range frequency inner difference and the phase frequency inner difference to generate an intersymbol deviation parameter and a phase deviation parameter of all frequency points of the full arc section of the GNSS satellite;
and S7, broadcasting the satellite orbit parameters, the clock error parameters, the delay parameters and the deviation parameters, and performing positioning calculation.
Has the advantages that:
according to the scheme of the invention, the method considers the difference relationship among the signals received by different ground stations, can realize the conversion of high-precision products among GNSS signals, further realize the real-time calculation of the high-precision products of the unmonitored signals, provide the high-precision products for the users corresponding to the unmonitored signals, further improve the navigation and positioning performance of the users, and provide high-precision space-time services for the GNSS users with different frequency points in the world.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the embodiments will be briefly described below. It is obvious that the drawings in the following description are only some embodiments of the invention, and that for a person skilled in the art, other drawings can be derived from them without inventive effort.
Fig. 1 is a flow chart schematically illustrating a method for converting a high-precision product between satellite navigation signals according to an embodiment of the present invention.
Detailed Description
The description of the embodiments of this specification is intended to be taken in conjunction with the accompanying drawings, which are to be considered part of the complete specification. In the drawings, the shape or thickness of the embodiments may be exaggerated and simplified or conveniently indicated. Further, the components of the structures in the drawings are described separately, and it should be noted that the components not shown or described in the drawings are well known to those skilled in the art.
Any reference to directions and orientations to the description of the embodiments herein is merely for convenience of description and should not be construed as limiting the scope of the invention in any way. The following description of the preferred embodiments refers to combinations of features which may be present individually or in combination, and the invention is not particularly limited to the preferred embodiments. The scope of the invention is defined by the claims.
The embodiment of the invention discloses a method for converting high-precision products among satellite navigation signals, which comprehensively utilizes a wide area monitoring station (such as a global range monitoring station or a large range monitoring station) to calculate full-arc-band full-frequency-band high-precision space-time service products of GNSS satellites according to observation data of partial frequency points of the GNSS satellites and observation data of all frequency points of the GNSS satellites by a local area monitoring station (referring to a monitoring station with a smaller range, wherein the distribution scale of the monitoring station is far smaller than the distribution range of the wide area monitoring station), and provides high-precision space-time service for GNSS users with different frequency points in the world.
According to the conception of the invention, the method utilizes the observation data of a wide area monitoring station on partial frequency points of the GNSS satellite to calculate the shared parameters of different frequency points such as high-precision satellite orbit, clock error parameters, ionosphere delay, troposphere delay and the like, and phase deviation parameters and inter-code deviation parameters of partial frequency points of the GNSS satellite in real time according to the description of the following steps S1, S2, S3, S4 and the like, then utilizes a local area monitoring station to calculate the inter-code deviation parameter difference and the phase deviation parameter difference of other frequency points and partial frequency points of the GNSS satellite on the observation data of all frequency points of the GNSS satellite, finally generates the inter-code deviation parameters and the phase deviation parameters of all frequency points of the full arc section of the GNSS satellite, and broadcasts the product to provide high-precision space-time information service for users of all frequency points. The partial frequency points described herein refer to navigation signals of GNSS satellite partial frequency points that can be observed by a wide area monitoring station. All frequency points refer to all navigation signals broadcast by GNSS satellites. The other frequency points refer to navigation signals except for the part of the frequency points of the GNSS satellite which can be observed by the wide area monitoring station. Other frequency points can be observed by local monitoring stations within a small range.
Referring to fig. 1, the flow of the method specifically includes:
s1, performing orbit determination on the GNSS satellite by using observation data of partial frequency points of the GNSS satellite by using a wide area monitoring station.
Further, the observation data in step S1 is an ionosphere-free combination of pseudo-range and carrier phase of a partial frequency point of the GNSS satellite. In the step S1, the observation data is used for performing statistical dynamics orbit determination and orbit prediction on the GNSS satellite, the initial orbit, the dynamics parameters and the measurement parameters of the GNSS satellite are solved, and the precise prediction orbit of the GNSS satellite is obtained.
Specifically, (1) a differential equation of a t satellite state transition matrix is constructed by taking a satellite broadcast ephemeris as an initial orbit parameter for satellite orbit determination processing.
Assume that the stateful vector characterizes the position and velocity of the satellite at time t:
Figure BDA0003823378840000041
satisfy a first order differential equation
Figure BDA0003823378840000042
Then it is determined that,
Figure BDA0003823378840000043
assume a state transition matrix of
Figure BDA0003823378840000044
From the above equation,. Phi. (t, t) can be obtained 0 ) Differential equation of (a):
Figure BDA0003823378840000045
the initial conditions were:
Φ(t 0 ,t 0 )=I 3*3
(2) Calculating a differential equation of a satellite sensitive matrix at the moment t by taking a satellite broadcast ephemeris as an initial orbit parameter for satellite orbit determination processing
The sensitivity matrix refers to the dynamic parameter p of the satellite state vector pair i (i=1,…,n p ) The dependency of (2):
Figure BDA0003823378840000051
the differential equation of the sensitive matrix gives the partial derivative relation of the state vector to the parameters of the dynamic model, and the calculation process is similar to that of the dynamic model, including:
Figure BDA0003823378840000052
or:
Figure BDA0003823378840000053
due to t 0 The satellite state vector at the moment is independent of the kinetic parameter P, so the initial condition of the sensitivity matrix is
S(t 0 )=0
(3) Constructing a variational equation according to the differential equation of the t satellite state transition matrix and the sensitive matrix, and performing orbit integration to obtain an approximate value of the satellite orbit at the current moment and partial derivatives of the satellite orbit initial orbit parameter and perturbation acceleration parameter
Combining the differential equations of the state transition matrix and the sensitivity matrix together can yield the following form:
Figure BDA0003823378840000054
the above formula is a first-order initial value problem, and can be solved by a numerical integration method. If a direct integration method of a second-order differential equation is adopted for solving, the form of a variational equation needs to be converted, and phi and S are firstly decomposed into:
Figure BDA0003823378840000055
Figure BDA0003823378840000061
due to the fact that
Figure BDA0003823378840000062
The variational equation can be written as:
Figure BDA0003823378840000063
since acceleration is independent of speed, i.e.
Figure BDA0003823378840000064
Zero means the right side of the second order formal variational equation and
Figure BDA0003823378840000065
is irrelevant. The calculation of the satellite orbit may be done using a numerical scoring method (e.g., a longkutta integrator).
(4) Circularly processing the available wide area monitoring station pseudo-range phase ionosphere-free combined observation data, combining the satellite orbit approximate value and partial derivative obtained in the step (3), constructing a method equation, and calculating the satellite initial orbit parameter and dynamic parameter estimation result by using least square estimation
When the precise orbit determination of the navigation satellite is processed, besides the initial orbit and the dynamic parameters of the satellite, the measurement parameter q is also included i (i=1,…,n q ). The measurement parameters comprise troposphere delay scale factors, phase data ambiguity parameters, satellite and observation station clock errors, inter-satellite bidirectional distance measurement and equipment time delay of satellite-ground bidirectional clock error measurement results.
For the observation data L (t), it can be a non-linear function of the state vector and the measured parameters of the satellite, such as:
L(t)=G(y(t),q)
the state vector of the satellite can be obtained by numerical integration according to the initial orbit and the dynamic parameters of the satellite, and the above formula can be obtained:
L(t)=G(y(t 0 ),p,q)
reference initial values y for the initial orbit, kinetic and measurement parameters of the satellite * (t 0 )、p * And q is * For the above formula linearization, we can obtain:
Figure BDA0003823378840000066
wherein, Δ y (t) 0 ) Initial orbit parameter of satellite relative to y for precise orbit determination * (t 0 ) Δ p is the correction of the kinetic parameter with respect to p ×, Δ q is the correction of the measurement parameter with respect to q ×.
Order to
Figure BDA0003823378840000071
x=(Δy(t 0 ) Δp Δq) T
l(t)=L(t)-G(y * (t 0 ),p*,q*)
The above equation can be written as:
v(t)=B(t)·x-l(t)
it can be seen that x is the parameter vector to be estimated for precise orbit determination, B is the coefficient matrix of the observed data L (t), and L (t) is the constant matrix of the observed data L (t).
For all observations, the observation equation is:
Figure BDA0003823378840000072
wherein, P i As observation data 1 (t) i ) The weight of (c).
In addition, the initial orbit of the satellite, the dynamic parameter and the reference initial value y of the measurement parameter are considered * (t 0 )、p * And q, the weights are:
v 0 =x-x 0 ,P 0
x 0 is the initial value of the parameter to be estimated. The construction matrix is as follows:
Figure BDA0003823378840000073
Figure BDA0003823378840000074
by using the principle of least squares, the optimal estimation of the parameters to be estimated is as follows:
Figure BDA0003823378840000075
and S2, utilizing the wide area monitoring station to carry out real-time filtering calculation on observation data of partial frequency points of the GNSS satellite to obtain clock error parameters.
Further, the observation data in step S2 is an ionosphere-free combination of pseudo-range and carrier phase of a partial frequency point of the GNSS satellite. Specifically, step S2 includes: and S21, on the basis of the step S1, fixing the GNSS satellite precision forecasting orbit, using non-ionosphere combination of pseudo range of partial frequency points of the GNSS satellite and carrier phase of a wide area monitoring station as input, resolving GNSS satellite clock error parameters by adopting Kalman filtering epoch-by-epoch, and synchronously and concomitantly estimating ambiguity parameters of carrier phase data. And S22, extrapolating a GNSS satellite precise orbit to generate a forecast orbit, simultaneously carrying out gross error elimination, cycle slip detection and restoration on the basis of a non-differential deionization layer combined observation value, fixing the extrapolated forecast orbit, introducing common parameters such as wide area monitoring station coordinates and ERP (enterprise resource planning), adding various corrections such as troposphere delay, station tide correction, satellite antenna phase center deviation, relativistic period term and relativistic gravity delay, carrying out parameter estimation until a set residual error threshold is met, and exiting iteration to obtain a precise GNSS satellite clock error parameter.
Wherein, step S21 includes: and S211, preprocessing pseudo range and carrier phase data of a receiver of the current wide area monitoring station, wherein the preprocessing comprises two parts of gross error elimination and cycle slip detection. The method mainly uses a MW combination of a non-difference double-frequency pseudo-range phase and an ionosphere residual combination observed value to remove gross error and detect and repair cycle slip. Although this method can detect most cycle slips and gross errors, there will still be some small cycle slips which will be handled at a later stage by a strict quality control method.
And step S212, for the moment without ionosphere combination of the pseudo range of the receiver of the wide area monitoring station and the carrier phase data obtained in the step S211, carrying out Lagrange interpolation on the forecast orbit to obtain a GNSS satellite orbit at the satellite launching moment, and then calculating the satellite-ground distance by combining with the accurate coordinate of the ground monitoring receiver to obtain a prior residual sequence.
And S213, processing the observation data currently using all the wide area monitoring stations in the step S212 to obtain a coefficient matrix and a constant matrix, constructing a normal equation matrix, and estimating parameters of clock error, atmosphere and ambiguity of the GNSS satellite in real time by using Kalman filtering.
And S3, decoding the observation data of partial frequency points of the GNSS satellite by using the wide area monitoring station to calculate the inter-code deviation parameter and the phase deviation parameter.
Further, the process of calculating the inter-code deviation parameter in step S3 includes: calculating an ionospheric residual combination of pseudo-range data of a wide area monitoring station on partial frequency points of a GNSS satellite, and separating ionospheric delay from code deviation parameters in an ionospheric delay modeling manner to obtain high-precision estimation of the code deviation parameters;
specifically, the pseudo-range observation equation of the GNSS satellite is shown as follows:
Figure BDA0003823378840000091
wherein, P i To be a pseudo-range observation,
Figure BDA0003823378840000092
is the geometric distance, δ t, between receiver j and satellite k j Is the clock error, δ t, of receiver j k Clock error, Δ t, of satellite k cor For errors correctable by the model, d trop In order to delay the tropospheric delay,
Figure BDA0003823378840000093
for the channel delay, tau, at the ith frequency point of the satellite k i,j The channel time delay of the ith frequency point of the receiver j, the TEC is the total electron content on an inclined path, f i Carrier frequency of ith frequency point, epsilon i C is the speed of light for the measurement error including multipath;
according to the formula, the difference between different frequency points can obtain a dual-frequency non-geometric combination:
Figure BDA0003823378840000094
wherein the content of the first and second substances,
Figure BDA0003823378840000095
IFB j =τ 1,j2,j
according to the ionosphere map GIM issued by the IGS, the TEC in the above equation can be used as a known quantity, so as to obtain a combined value of the GNSS satellite and the receiver hardware delay:
Figure BDA0003823378840000096
according to a predetermined spherical harmonic model (low order or high order), the TEC is regarded as a known quantity including an ionosphere model parameter unknown; and establishing an equation set for hardware delay combination values obtained by all GNSS satellites and all receivers, wherein the equation set is rank-deficient one, and proper constraints are added to enable the equation set to be full-rank.
Further, the process of calculating the phase deviation parameter in step S3 includes: calculating a high-precision phase deviation parameter by taking the carrier phase data ambiguity parameter obtained in the step S2 as input and adopting a multi-station multi-satellite combined calculation mode;
specifically, ambiguity fixing of a system end is carried out together with orbit and clock error calculation, input data of orbit determination calculation comprise an initial orbit, an initial clock error, a plurality of globally distributed receiver observation data and accurate coordinates of a receiver, and non-differential ionosphere-free combination is adopted as basic observed quantity;
firstly, a satellite orbit and clock error floating solution is calculated, ambiguity fixing is not carried out in the step, and an ambiguity parameter result is in a real number form:
fitting the initial orbit to obtain the initial orbit number and the kinetic model parameters, and performing orbit integration to obtain a reference orbit, a state transition matrix and a sensitive matrix;
preprocessing observed data, eliminating obvious unreasonable gross errors, detecting and repairing cycle slip of phase observed data, performing antenna phase center correction, antenna phase winding correction, troposphere correction, relativity correction and tide correction on the observed data, and forming pseudo-range ionosphere-free combined observed quantity and phase ionosphere-free combined observed quantity by utilizing dual-frequency observed data;
establishing an observation equation according to a reference orbit, an initial clock error, a survey station coordinate and a non-ionosphere combined observed quantity, calculating an O-C value, further performing parameter estimation to obtain new orbit parameters and clock error parameters, and eliminating gross errors through 5-7 iterations to obtain a floating solution orbit and clock error result with better precision;
then, calculating a satellite orbit and clock error fixing solution, and estimating a satellite UPD parameter: integrating the floating-point solution orbit to obtain a reference orbit, a state transition matrix and a sensitive matrix with higher precision;
decomposing the non-ionized layer combined ambiguity obtained in the floating solution into a wide lane ambiguity and a narrow lane ambiguity, and selecting an independent base line to form a station-satellite double difference to enable the ambiguity to have a whole-cycle characteristic; because the wavelength of the ambiguity of the wide lane reaches 0.86m, the ambiguity of the wide lane is easy to fix, the ambiguity of the wide lane is obtained by calculation of MW combined observed quantity, in order to weaken pseudo-range noise and multi-path influence, after the ambiguity of the wide lane is subjected to inter-epoch smoothing, the ambiguity of the wide lane is locally fixed, and the ambiguity of the narrow lane is fixed by search through an LAMDA method;
constructing an observation equation according to the floating-point solution orbit clock error, the station coordinates and the ionosphere-free combined observed quantity, calculating an O-C value, performing parameter estimation by taking a double-difference integer ambiguity fixed result as a constraint condition to obtain new orbit parameters and clock error parameters, improving integer ambiguity fixed success rate through 3-5 iterations, and eliminating gross errors to obtain final fixed solution orbit and clock error;
finally, calculating UPD parameters of the GNSS satellite and the receiver according to the ambiguity fixing result, selecting a certain satellite or a certain receiver as a reference datum for solving the rank deficiency problem during resolving, and considering the UPD parameter as 0 or selecting the UPD mean value of the whole constellation as 0 as the datum; and converting the obtained UPD parameters of the wide and narrow lanes to obtain UPD parameter estimation values corresponding to two frequency points of the satellite.
And S4, fixing the GNSS satellite precision forecast orbit obtained in the step S1, the GNSS precision clock error obtained in the step S2 and the inter-code deviation parameter and the phase deviation parameter obtained in the step S3, performing precision point positioning processing on observation data of partial frequency points of the GNSS satellite by using a wide area monitoring station, and solving an ionosphere delay parameter and a troposphere delay parameter by adopting a real-time filtering method.
Further, the observation data in step S4 is a multifrequency non-differential combination of pseudo-range of the wide area monitoring station to partial frequency points of the GNSS satellite and carrier phase data.
And S4, after the resolving is finished, carrying out grid processing on the ionosphere delay parameter and the troposphere delay parameter in real time, and generating an ionosphere delay model background field, a troposphere delay model background field, grid-processed residual ionosphere inclined path delay and troposphere delay.
Specifically, by performing error correction on the original pseudo range and carrier observed quantity, and simultaneously taking the precise orbit and clock error as known values and linearizing the equation, the pseudo range and carrier observed equation of the non-combined PPP of any satellite navigation system can be simplified into
Figure BDA0003823378840000111
Figure BDA0003823378840000112
Figure BDA0003823378840000113
Figure BDA0003823378840000114
Figure BDA0003823378840000115
In the formula (I), the compound is shown in the specification,
Figure BDA0003823378840000116
and
Figure BDA0003823378840000117
to add various error-corrected pseudoranges and carrier observations, gamma i =f 1 2 /f i 2
Figure BDA0003823378840000118
The ionospheric delay of the L1 bin includes hardware Delay (DCB) at the receiver side and at the satellite side, wherein,
Figure BDA0003823378840000119
the pseudo-range hardware delay bias is respectively at the receiver end and the satellite end.
When the non-difference non-combination positioning is carried out, the processing mode is consistent with the processing mode of the combination of the deionization layer, the receiver clock error of other systems is expressed in the form of the receiver clock error plus the system time deviation, and the non-combination PPP pseudo range and the carrier equation of the satellite navigation system s can be further expressed in the form of
Figure BDA00038233788400001110
Figure BDA00038233788400001111
Figure BDA00038233788400001112
Figure BDA00038233788400001113
Figure BDA00038233788400001114
Figure BDA00038233788400001115
In the above formula, the system time deviation parameter is consistent with the ISB parameter form in the ionospheric elimination combination model, and meanwhile, the above formula uniformly summarizes the non-combination PPP model of any satellite navigation system, and is suitable for any satellite navigation system and multi-system combination PPP.
In summary, the non-combined PPP model and the estimated parameters are
Figure BDA0003823378840000121
Figure BDA0003823378840000122
Figure BDA0003823378840000123
Figure BDA0003823378840000124
Figure BDA0003823378840000125
At the moment, X is a parameter vector to be estimated in the equation, and parameters in the X vector are respectively a station coordinate correction number, a receiver clock error, a system time deviation containing hardware delay influence, zenith troposphere delay, ionosphere delay and ambiguity.
And S5, resolving pseudo range intra-frequency differences and phase intra-frequency differences of partial frequency points and all other frequency points by using observation data of all frequency points of the GNSS satellite through the local monitoring station.
Further, the observation data in step S5 is a multifrequency non-differential combination of pseudo ranges of all frequency points of the GNSS satellites and carrier phase data of the wide area monitoring station.
For example, the Beidou system B1 frequency point is taken as an example, and the assumption is made
Figure BDA0003823378840000126
And
Figure BDA0003823378840000127
respectively representing pseudo-range observation of the receiver r on the satellite I at the frequency point B1C and the frequency point B1I, and performing difference on the pseudo-range observation to obtain:
Figure BDA0003823378840000128
in the above formula, the first and second carbon atoms are,
Figure BDA0003823378840000129
is the intra-frequency difference parameter between the B1C frequency point of the satellite I and the B1I frequency point,
Figure BDA00038233788400001210
is the intra-frequency difference parameter between the B1 frequency point and the B1I frequency point of the receiver r.
And obtaining intra-frequency difference parameters of all satellites and the B1 frequency point and the B1I frequency point of the receiver by utilizing least square estimation through an observation simultaneous observation equation set of a plurality of monitoring receivers for a plurality of satellites.
An observation equation of the phase difference between the B1C frequency point and the B1I frequency point is as follows:
Figure BDA00038233788400001211
similarly, pseudo-range intra-frequency difference and phase intra-frequency difference parameters of all frequency points of the satellite i can be obtained by constructing an observation equation and utilizing least square estimation.
And S6, forecasting the pseudo range frequency internal difference and the phase frequency internal difference obtained in the step S5, and generating the code deviation parameters and the phase deviation parameters of all frequency points of the full arc section of the GNSS satellite.
Specifically, the original observation values of the Chinese area monitoring stations are subjected to difference to directly obtain pseudo range intra-frequency difference and phase intra-frequency difference, and then a time sequence forecasting method is used for obtaining the code deviation and the phase deviation value of the full arc section. A time series is a sequence that varies randomly over time. The current time series models mainly comprise 3 types: an AR (p) model (autoregressive model), an MA model (moving average model), and an ARMA model (autoregressive moving average model or hybrid model), wherein the AR (p) model is a well-developed model in the predictive model.
The AR (p) model may be defined as having a constant, normal, zero-mean timing { yt }, with
Figure BDA0003823378840000131
In the formula, y i (i = t, t-1, …, t-p) represents the value at time i;
Figure BDA0003823378840000132
is a coefficient; ε (k) is subject to NID (0, σ) ε 2 ) White noise of (2); p is the order.
Because the AR (p) model analyzes and forecasts the stabilized time sequence or the quasi-stabilized time sequence, and random fluctuation caused by random noise in the short-term forecast of the clock error time sequence possibly plays an important role (particularly a clock with poor stability). According to the characteristic that the satellite clock error accords with the characteristics of a quadratic polynomial in a short-term range, the modeling idea of the invention is as follows: a quadratic polynomial model is adopted to extract trend terms when the data are stabilized, and a difference AR (p) model is adopted to model random terms. Therefore, even for the clock with the clock difference trend not linearly changing, the AR (p) model can also model the residual error of the clock with the clock difference trend, so that the characteristic of improving the clock difference forecasting precision is achieved.
And S7, broadcasting the satellite orbit parameters, the clock error parameters, the delay parameters and the deviation parameters, and performing positioning calculation. Here, any communication format may be selected for distribution.
The sequence numbers of the above steps related to the method of the present invention do not mean the order of execution of the method, and the order of execution of the steps should be determined by their functions and inherent logic, and should not limit the implementation process of the embodiment of the present invention.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the present invention, and any modifications, equivalents, improvements and the like made within the spirit and principle of the present invention should be included in the scope of the present invention.

Claims (10)

1. A high-precision product conversion method between satellite navigation signals comprises the following steps:
s1, performing orbit determination on the GNSS satellite by utilizing observation data of partial frequency points of the GNSS satellite by using a wide area monitoring station;
s2, utilizing a wide area monitoring station to carry out real-time filtering resolving on observation data of partial frequency points of the GNSS satellite to obtain clock error parameters;
s3, decoding observation data of partial frequency points of the GNSS satellite by using a wide area monitoring station to calculate an inter-code deviation parameter and a phase deviation parameter;
s4, performing precise point positioning processing on observation data of partial frequency points of the GNSS satellite by using a wide area monitoring station, and solving ionosphere delay parameters and troposphere delay parameters by adopting a real-time filtering method;
s5, resolving pseudo range frequency internal differences and phase frequency internal differences of partial frequency points and all other frequency points by using a local monitoring station to observation data of all frequency points of the GNSS satellite;
s6, forecasting the pseudo range frequency internal difference and the phase frequency internal difference to generate an inter-code deviation parameter and a phase deviation parameter of all frequency points of a full arc section of the GNSS satellite;
and S7, broadcasting the satellite orbit parameter, the clock error parameter, the delay parameter and the deviation parameter, and performing positioning calculation.
2. The method of claim 1, wherein the observed data in step S1 and step S2 is an ionosphere-free combination of pseudoranges and carrier phases for partial frequency points of GNSS satellites.
3. The method according to claim 2, wherein in step S1, the observation data is used to perform statistical dynamics orbit determination and orbit prediction on the GNSS satellite, and the initial orbit, the dynamics parameters and the measurement parameters of the GNSS satellite are solved to obtain the precision predicted orbit of the GNSS satellite.
4. The method of claim 3, wherein the measurement parameters include tropospheric delay scale factors, phase data ambiguity parameters, device time delays for satellite and monitoring station clock bias, inter-satellite two-way ranging, and two-way-satellite-to-ground clock bias measurements.
5. The method according to claim 3, wherein the step S2 comprises:
s21, fixing the GNSS satellite precision forecasting orbit, using non-ionosphere combination of pseudo-range of partial frequency points of the wide area monitoring station to the GNSS satellite and carrier phase as input, adopting Kalman filtering to calculate clock error parameters of the GNSS satellite one epoch by one epoch, and synchronously estimating ambiguity parameters of carrier phase data along with the clock error parameters;
s22, extrapolating a GNSS satellite precise orbit to generate a forecast orbit, simultaneously performing gross error elimination, cycle slip detection and restoration on the basis of a non-difference deionization layer combined observed value, fixing the extrapolated forecast orbit, introducing wide area monitoring station coordinates and ERP common parameters, adding troposphere delay, station tide correction, satellite antenna phase center deviation, relativistic period item and relativistic gravity delay, performing parameter estimation until a set residual error threshold value is met, and exiting iteration to obtain a precise GNSS satellite clock error parameter.
6. The method according to claim 5, wherein the step S21 comprises:
s211, preprocessing pseudo range and carrier phase data of a receiver of the current wide area monitoring station, wherein the preprocessing comprises removing gross error and detecting and repairing cycle slip through MW combination of non-difference double-frequency pseudo range phases and ionosphere residual error combination observed values;
s212, performing Lagrange interpolation on the forecast orbit to obtain a GNSS satellite orbit at the satellite emission moment at the moment of non-ionosphere combination of the pseudo range of the receiver of the wide area monitoring station and the carrier phase data obtained in the step S211, and then calculating the satellite-ground distance by combining with the accurate coordinate of the ground monitoring receiver to obtain a prior residual sequence;
s213, processing the current observation data of all the wide area monitoring stations in the step S212 to obtain a coefficient matrix and a constant matrix, constructing a normal equation matrix, and estimating parameters of clock error, atmosphere and ambiguity of the GNSS satellite in real time by using Kalman filtering.
7. The method according to claim 1, wherein the calculating process of the inter-code bias parameter in the step S3 comprises: calculating an ionospheric residual combination of pseudo-range data of a wide area monitoring station on partial frequency points of a GNSS satellite, and separating ionospheric delay from code deviation parameters by adopting an ionospheric delay modeling mode to obtain high-precision estimation of the code deviation parameters;
the pseudo-range observation equation of the GNSS satellite is shown as follows:
Figure FDA0003823378830000021
wherein, P i To be a pseudo-range observation,
Figure FDA0003823378830000022
is the geometric distance, δ t, between the receiver j and the satellite k j Is the clock difference, δ t, of receiver j k Clock error, Δ t, of satellite k cor For errors correctable by the model, d trop In order to delay the tropospheric delay,
Figure FDA0003823378830000023
for the channel delay, tau, at the ith frequency point of the satellite k i,j The channel time delay of the ith frequency point of the receiver j, the TEC is the total electron content on an inclined path, f i Carrier frequency of ith frequency point,ε i C is the speed of light for the measurement error including multipath;
according to the formula, the difference between different frequency points can obtain a double-frequency non-geometric combination:
Figure FDA0003823378830000031
wherein the content of the first and second substances,
Figure FDA0003823378830000032
IFB j =τ 1,j2,j
according to the ionosphere map GIM issued by the IGS, the TEC in the above equation can be used as a known quantity, so as to obtain a combined value of the GNSS satellite and the receiver hardware delay:
Figure FDA0003823378830000033
according to a predetermined spherical harmonic model (low order or high order), the TEC is considered as a known quantity containing ionospheric model parameter unknowns; and establishing an equation set for hardware delay combination values obtained by all GNSS satellites and all receivers, wherein the equation set is rank-deficient one, and proper constraints are added to enable the equation set to be full-rank.
8. The method according to claim 1, wherein the calculating process of the phase deviation parameter in the step S3 comprises: calculating a high-precision phase deviation parameter by taking the carrier phase data ambiguity parameter obtained in the step S2 as input and adopting a multi-station multi-satellite combined calculation mode;
the ambiguity fixing of a system end is carried out together with orbit and clock error calculation, input data of orbit determination calculation comprise an initial orbit, an initial clock error, observation data of a plurality of receivers distributed globally and accurate coordinates of the receivers, and a non-differential non-ionosphere combination is adopted as a basic observation quantity;
firstly, a satellite orbit and clock error floating solution is calculated, ambiguity fixing is not carried out in the step, and an ambiguity parameter result is in a real number form:
fitting the initial orbit to obtain the initial orbit root and the kinetic model parameters, and performing orbit integration to obtain a reference orbit, a state transition matrix and a sensitive matrix;
preprocessing observed data, eliminating obvious unreasonable gross errors, detecting and repairing cycle slip of phase observed data, performing antenna phase center correction, antenna phase winding correction, troposphere correction, relativity correction and tide correction on the observed data, and forming pseudo-range ionosphere-free combined observed quantity and phase ionosphere-free combined observed quantity by utilizing dual-frequency observed data;
establishing an observation equation according to a reference orbit, an initial clock error, a survey station coordinate and a non-ionosphere combined observed quantity, calculating an O-C value, further performing parameter estimation to obtain new orbit parameters and clock error parameters, and eliminating gross errors through 5-7 iterations to obtain a floating solution orbit and clock error result with better precision;
then, calculating a satellite orbit and clock error fixing solution, and estimating a satellite UPD parameter: integrating the floating-point solution orbit to obtain a reference orbit, a state transition matrix and a sensitive matrix with higher precision;
decomposing the non-ionized layer combined ambiguity obtained in the floating solution into a wide lane ambiguity and a narrow lane ambiguity, and selecting an independent base line to form a station-satellite double difference to enable the ambiguity to have a whole-cycle characteristic; because the wavelength of the ambiguity of the wide lane reaches 0.86m, the ambiguity of the wide lane is easy to fix, the ambiguity of the wide lane is obtained by calculation of MW combined observed quantity, in order to weaken pseudo-range noise and multi-path influence, after the ambiguity of the wide lane is subjected to inter-epoch smoothing, the ambiguity of the wide lane is locally fixed, and the ambiguity of the narrow lane is fixed by search through an LAMDA method;
constructing an observation equation according to the floating-point solution orbit clock error, the station coordinates and the ionosphere-free combined observed quantity, calculating an O-C value, performing parameter estimation by taking a double-difference integer ambiguity fixed result as a constraint condition to obtain new orbit parameters and clock error parameters, improving integer ambiguity fixed success rate through 3-5 iterations, and eliminating gross errors to obtain final fixed solution orbit and clock error;
finally, calculating UPD parameters of the GNSS satellite and the receiver according to the ambiguity fixing result, selecting a certain satellite or a certain receiver as a reference datum for solving the rank deficiency problem during resolving, and considering the UPD parameter as 0 or selecting the UPD mean value of the whole constellation as 0 as the datum; and converting the obtained UPD parameters of the wide and narrow lanes to obtain UPD parameter estimation values corresponding to two frequency points of the satellite.
9. The method according to claim 1, wherein the observation data in step S4 and step S5 is multi-frequency non-differential combinations of pseudo-range and carrier phase data of the wide area monitoring station for partial frequency points and all frequency points of the GNSS satellite, respectively.
10. The method according to claim 1, wherein after the step S4 of completing the calculation, the ionospheric delay parameter and the tropospheric delay parameter are subjected to gridding processing in real time to generate a background field of an ionospheric delay model and a background field of a tropospheric delay model, and residual ionospheric diagonal path delay and tropospheric delay of the gridding.
CN202211049892.5A 2022-08-30 2022-08-30 High-precision product conversion method between satellite navigation signals Pending CN115373005A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211049892.5A CN115373005A (en) 2022-08-30 2022-08-30 High-precision product conversion method between satellite navigation signals

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211049892.5A CN115373005A (en) 2022-08-30 2022-08-30 High-precision product conversion method between satellite navigation signals

Publications (1)

Publication Number Publication Date
CN115373005A true CN115373005A (en) 2022-11-22

Family

ID=84069791

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211049892.5A Pending CN115373005A (en) 2022-08-30 2022-08-30 High-precision product conversion method between satellite navigation signals

Country Status (1)

Country Link
CN (1) CN115373005A (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115877415A (en) * 2022-12-07 2023-03-31 中国科学院国家授时中心 Method and system for estimating clock difference of ultrahigh-precision space spacecraft
CN116243341A (en) * 2022-12-22 2023-06-09 国汽大有时空科技(安庆)有限公司 Nationwide integrated PPP-RTK service system construction method, device and system
CN116540280A (en) * 2023-07-07 2023-08-04 中国科学院空天信息创新研究院 Comprehensive processing method and system for state domain correction information of multi-frequency satellite navigation data
CN116540278A (en) * 2023-07-06 2023-08-04 中国科学院空天信息创新研究院 Cloud edge end cooperative reference dynamic maintenance method and system

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115877415A (en) * 2022-12-07 2023-03-31 中国科学院国家授时中心 Method and system for estimating clock difference of ultrahigh-precision space spacecraft
CN115877415B (en) * 2022-12-07 2023-08-18 中国科学院国家授时中心 Ultra-high precision space spacecraft clock error estimation method and system
CN116243341A (en) * 2022-12-22 2023-06-09 国汽大有时空科技(安庆)有限公司 Nationwide integrated PPP-RTK service system construction method, device and system
CN116243341B (en) * 2022-12-22 2023-12-05 国汽大有时空科技(安庆)有限公司 Nationwide integrated PPP-RTK service system construction method, device and system
CN116540278A (en) * 2023-07-06 2023-08-04 中国科学院空天信息创新研究院 Cloud edge end cooperative reference dynamic maintenance method and system
CN116540278B (en) * 2023-07-06 2023-09-08 中国科学院空天信息创新研究院 Cloud edge end cooperative reference dynamic maintenance method and system
CN116540280A (en) * 2023-07-07 2023-08-04 中国科学院空天信息创新研究院 Comprehensive processing method and system for state domain correction information of multi-frequency satellite navigation data
CN116540280B (en) * 2023-07-07 2023-09-15 中国科学院空天信息创新研究院 Comprehensive processing method and system for state domain correction information of multi-frequency satellite navigation data

Similar Documents

Publication Publication Date Title
CN110231037B (en) GNSS maneuvering satellite orbit determination method with additional clock error model constraint
CN109709591B (en) GNSS high-precision positioning method for intelligent terminal
EP3130943B1 (en) Navigation satellite system positioning involving the generation of tropospheric correction information
EP3035080B1 (en) Navigation satellite system positioning involving the generation of correction information
CN110275186B (en) LEO satellite enhanced GNSS ionosphere normalization and fusion modeling method
US9557419B2 (en) Methods for generating accuracy information on an ionosphere model for satellite navigation applications
CN115373005A (en) High-precision product conversion method between satellite navigation signals
CN107966722B (en) GNSS clock error resolving method
CN109143298B (en) Beidou and GPS observation value cycle slip detection and restoration method, equipment and storage equipment
CN111694030A (en) BDS local difference method and system based on grid virtual observation value
CN113358017B (en) Multi-station cooperative processing GNSS high-precision deformation monitoring method
Cheng et al. Statistical analysis and quality control for GPS fractional cycle bias and integer recovery clock estimation with raw and combined observation models
Chen et al. A geometry-free and ionosphere-free multipath mitigation method for BDS three-frequency ambiguity resolution
Landau et al. Trimble’s RTK and DGPS solutions in comparison with precise point positioning
CN111983641A (en) Method for generating Beidou satellite-based augmentation system integrity parameters in real time
CN114859389A (en) GNSS multi-system robust adaptive fusion RTK resolving method
Elsayed et al. Bounding of correlated double-differenced GNSS observation errors using NRTK for precise positioning of autonomous vehicles
CN117388883A (en) Beidou low-orbit PPP-RTK high-precision service method based on sparse foundation nodes
Wang et al. Comparison of three widely used multi‐GNSS real‐time single‐frequency precise point positioning models using the International GNSS Service real‐time service
Bisnath Relative Positioning and Real‐Time Kinematic (RTK)
Luo et al. Benefit of sparse reference network in BDS single point positioning with single-frequency measurements
CN115902968A (en) PPP terminal positioning method based on Beidou third GEO broadcast enhancement information
Li et al. A grid-based ionospheric weighted method for PPP-RTK with diverse network scales and ionospheric activity levels
Rovira-Garcia et al. A real-time world-wide ionospheric model for single and multi-frequency precise navigation
CN112528213B (en) Global ionosphere total electron content multilayer analysis method based on low earth orbit satellite

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