CN104714244A - Multi-system dynamic PPP resolving method based on robust self-adaption Kalman smoothing - Google Patents

Abstract

The invention discloses a multi-system dynamic PPP resolving method based on robust self-adaption Kalman smoothing. The method includes the steps that receiving machine outline coordinates and receiving machine clock bias of all systems are solved through selecting-weight-iteration pseudo-range single-point positioning, and accordingly all positioning error correction values are calculated according to an error correction model in combination with the satellite precise ephemeris and satellite precise clock bias; strict data quality control is conducted on observation data. Due to the fact that dynamic PPP accuracy is easily affected by undetected small cycle slips or the gross error and the like, an observation equation weight matrix is adjusted according to the observation value residual vectors, and the undetected small cycle slips or the gross error and other influence factors are removed; self-adaption factors are determined according to the state predictive information, and thus the influence on parameter estimation of the predictive information is controlled. By means of the method, when multi-system dynamic PPP is conducted through a single receiving machine, the feature that the number of multi-system satellites is increased greatly, on the basis that the stability of the satellite structure is guaranteed, the influence of the gross error is weaken effectively, the dynamic noise abnormity in dynamic positioning is improved, and finally the high-precision and high-stability multi-system dynamic PPP result is achieved.

Description

The dynamic PPP calculation method of a kind of multisystem based on robust adaptable Kalman filter

Technical field

The present invention relates to the dynamic PPP calculation method of a kind of multisystem, belong to GNSS dynamic precision One-Point Location field.

Background technology

Satellite positioning tech has become modern topmost navigator fix mode, no matter is in engineering survey, productive life or Military Application, all plays an important role.The ultimate principle of satnav is the distance between the satellite of measuring known location to receiver user, then adopts the method for distance resection according to the instantaneous position of multi-satellite, determines the position coordinates of tested point.According to the locator meams of satellite receiver devices, GNSS location technology can be divided into One-Point Location and Differential positioning.Traditional satellite One-Point Location technology refers to and utilizes pseudorange code to position, i.e. GNSS absolute fix, and this technology can adopt a receiver to complete real-time measurement.In actual applications, because pseudorange coded signal code length is larger, itself precision is lower, and by the satellite position that calculates in conjunction with broadcast ephemeris and the factor such as satellite clock correction precision is poor, pseudorange One-Point Location is difficult to reach higher positioning precision, be generally used for accuracy requirement lower, in the navigation application that requirement of real-time is higher.

Utilize carrier phase observation data to carry out carrier difference location, higher positioning precision can be reached, but need at least two receivers to carry out simultaneous observation, the positioning service of centimetre-sized could be realized by Double deference model.Along with the increase of receiver spacing, error correlativity reduces gradually, and arriving certain distance cannot weaken error effect.This location technology not only increases cost and the complexity of operation, and is also restricted in a lot of application scenario.The appearance of Static Precise Point Positioning (Precise Point Positioning, PPP) overcomes many defects of Differential positioning.PPP is proposed in 1997 by the people such as Zumberge of U.S. jet propulsion laboratory (JPL), and realizes on the data processing software GIPSY of their exploitation.PPP has the Static positioning accuracy of centimetre-sized and the dynamic locating accuracy of decimeter grade, can thoroughly break away from a large scale, long range measurements to the dependence of ground reference station, significantly improve operating efficiency, saved user cost.PPP is not only the study hotspot of field of locating technology in recent years, is also widely used in engineer applied fields such as various control survey, engineering survey, topographical surveyings.

In dynamic PPP location, the receiver coordinate of each epoch is in continuous change, and localizing environment is complicated.Utilize multisystem combined location, usable satellite number can be increased considerably, strengthen satellite geometry configuration, meet Kinematic Positioning demand further.For the object doing irregular movement, its accurate function model and probabilistic model are all difficult to obtain, and on the basis of traditional Kalman filter, adopt Robust filter to weaken the larger information of observation residual error, and utilize predicted state vector to determine adaptive factor, to obtain reliable Kinematic Positioning result.

Summary of the invention

Goal of the invention: in order to overcome the deficiencies in the prior art, the invention provides the dynamic PPP calculation method of a kind of multisystem based on robust adaptable Kalman filter, can effectively solve in current Kinematic Positioning, due to function model and probabilistic model inaccurate, and the phenomenon that satellite spatial structural instability causes positioning result poor to result.

Technical scheme: for achieving the above object, the technical solution used in the present invention is: the dynamic PPP calculation method of a kind of multisystem based on robust adaptable Kalman filter, first use three system broadcasts ephemeris computation co-ordinates of satellite and satellite clock correction separately, space-time datum unification is carried out to it; In conjunction with observation file in Pseudo-range Observations combination gained without ionosphere Pseudo-range Observations, carry out Iterated adjustment pseudorange One-Point Location, inverse reception epoch machine rough coordinates and each system receiver clock correction; Analyze from IGS etc. High Precision Satellite Ephemeris and the precise clock correction product that center website obtains each system, calculate each system-satellite precision coordinate and satellite precise clock correction by lagrange-interpolation; Carry out each Error Correction Model calculation of position errors corrected value subsequently, observation file, satellite precise coordinate and satellite precise clock correction also combines the quality control that observation file carries out observation data; Finally set up the dynamic PPP location model of multisystem, adopt robust adaptable Kalman filter to calculate residual values and the information of forecasting of observed reading information, determine weight of observation battle array and adaptive factor, realize the dynamic PPP location of high precision, high stability.

Specifically comprise the steps:

Step 1, uses three system broadcasts ephemeris computation co-ordinates of satellite and satellite clock correction separately, carries out space-time datum unification to it, obtains the co-ordinates of satellite after unified benchmark and satellite clock correction; According to P1 (C1) and P2 (C2) Pseudo-range Observations of observation file, carry out obtaining without ionospheric combination Pseudo-range Observations without ionospheric combination; According to L1 and the L2 observed reading of observation file, carry out obtaining without ionospheric combination carrier observations without ionospheric combination;

Step 2, what obtain according to step 1 carries out Iterated adjustment pseudorange One-Point Location without the co-ordinates of satellite after ionospheric combination Pseudo-range Observations and unified benchmark and satellite clock correction, obtains receiver rough coordinates and receiver clock-offsets;

Step 3, obtains three system-satellite precise ephemerises and satellite precise clock correction by network from analytic centres such as IGS, according to moment epoch of observation file, carries out satellite precise coordinate and satellite precise clock correction that Lagrange's interpolation obtains corresponding moment epoch;

Step 4, the satellite precise coordinate that the receiver rough coordinates that the receiver information, the step 2 that provide according to observation file obtain and receiver clock-offsets and step 3 obtain and precise clock correction, calculated every Correction of Errors value of Static Precise Point Positioning process by Error Correction Model, and carry out the quality control of observation data data in conjunction with observation data;

Step 5, according to step 1 obtain without ionospheric combination Pseudo-range Observations with without ionospheric combination carrier observations, the receiver rough coordinates that step 2 obtains and receiver clock-offsets, the every Correction of Errors value of calculating gained that the satellite precise coordinate that step 3 obtains and satellite precise clock correction and step 4 obtain, in conjunction with Static Precise Point Positioning model, build multisystem dynamic precision One-Point Location equation;

Step 6, adopt the multisystem dynamic precision One-Point Location equation that robust adaptive Kalman filter process of solution 5 obtains, according to the status information of the observation information obtained in filtering and prediction, the weights of adjustment observed reading and calculating adaptive factor, realize the Kinematic Positioning of high precision and high stability.

Obtain without ionospheric combination Pseudo-range Observations and the formula without ionospheric combination carrier observations without ionospheric combination in described step 1:

$\begin{array}{c}{P}_{\mathrm{IF}}=\frac{f_1^2\·P_1-f_2^2\·P_2}{f_1^2-f_2^2}\\ =\mathrm{\ρ}+\mathrm{cli}\·({\mathrm{dt}}_r-{\mathrm{dt}}^s)+\mathrm{dtrop}+{\mathrm{\ϵ}}_{P_{\mathrm{IF}}}\end{array}$

$\begin{array}{c}{L}_{\mathrm{IF}}=\frac{f_1^2\·L_1-f_2^2\·L_2}{f_1^2-f_2^2}\\ =\mathrm{\ρ}+\mathrm{cli}\·({\mathrm{dt}}_r-{\mathrm{dt}}^s)+\mathrm{dtrop}+\frac{f_1^2\·{\mathrm{\λ}}_1\·N_1-f_2^2\·{\mathrm{\λ}}_2\·N_2}{f_1^2-f_2^2}+{\mathrm{\ϵ}}_{L_{\mathrm{IF}}}\end{array}$

Wherein, P _iFindicate without ionospheric combination Pseudo-range Observations; f ₁, f ₂represent the carrier observations frequency of GPS/GLONASS/BDS system; P ₁, P ₂represent the Pseudo-range Observations of GPS/GLONASS/BDS system; If there is no P ₁/ P ₂during observed reading, C can be adopted ₁/ C ₂observed reading carries out using after DCB corrects; ρ is the distance of satellite and survey station; Cli represents the light velocity; Dt _rfor representing receiver clock-offsets; Dt ^srepresent satellite clock correction; d _tropit is tropospheric delay; L _iFindicate without ionospheric combination carrier observations; L ₁, L ₂represent the carrier phase observation data of GPS/GLONASS/BDS system; λ ₁, λ ₂represent the carrier phase wavelength of GPS/GLONASS/BDS system respectively, N ₁, N ₂represent respectively GPS/GLONASS/BDS system without ionosphere blur level; represent that pseudorange and carrier wave are without ionospheric combination observation noise respectively.

The satellite precise coordinate that described step 2 obtains and satellite precise clock correction:

$P_n\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{l}_i\left(x\right)\·y_i=\underset{i=0}{\overset{n}{\mathrm{\Σ}}}(y_i\·\underset{j=0,j\≠i}{\overset{n}{\mathrm{\Π}}}\frac{x-x_j}{x_i-x_j})$

In formula, P _nx () represents the satellite information in interpolation moment, i.e. the satellite orbit coordinate of each epoch and clock correction value after interpolation, and x is the interpolation moment, x _i(i=0,1 ..., n) be the node corresponding moment; y _i(i=0,1 ..., n) be satellite orbit coordinate corresponding to node time instance and clock correction value; l _i(x) (i=0,1 ..., n) be expressed as n interpolation basic differential polynomial.

The multisystem dynamic precision One-Point Location observation equation that described step 5 obtains is:

$P_{\mathrm{IF}}^g={\mathrm{\ρ}}^g+\mathrm{cli}\·{\mathrm{dt}}^g+d_{\mathrm{trop}}^g+{\mathrm{\ϵ}}_{P_{\mathrm{IF}}}^g$

${\mathrm{\Φ}}_{\mathrm{IF}}^g={\mathrm{\ρ}}^g+\mathrm{cli}\·{\mathrm{dt}}^g+d_{\mathrm{trop}}^g+N_{\mathrm{IF}}^g+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_{\mathrm{IF}}}^{\′g}$

$P_{\mathrm{IF}}^r={\mathrm{\ρ}}^r+\mathrm{cli}\·{\mathrm{dt}}^r+d_{\mathrm{trop}}^r+{\mathrm{\ϵ}}_{P_{\mathrm{IF}}}^r$

${\mathrm{\Φ}}_{\mathrm{IF}}^r={\mathrm{\ρ}}^r+\mathrm{cli}\·{\mathrm{dt}}^r+d_{\mathrm{trop}}^r+N_{\mathrm{IF}}^r+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_{\mathrm{IF}}}^r$

$P_{\mathrm{IF}}^c={\mathrm{\ρ}}^c+\mathrm{cli}\·{\mathrm{dt}}^c+d_{\mathrm{trop}}^g+{\mathrm{\ϵ}}_{P_{\mathrm{IF}}}^c$

${\mathrm{\Φ}}_{\mathrm{IF}}^c={\mathrm{\ρ}}^c+\mathrm{cli}\·{\mathrm{dt}}^c+d_{\mathrm{trop}}^g+N_{\mathrm{IF}}^c+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_{\mathrm{IF}}}^c$

In formula, subscript g, r, c represent GPS, GLONASS and BDS system respectively; Cli is the light velocity; P _iFbe after Correction of Errors without ionospheric combination Pseudo-range Observations; Φ _iFbe after Correction of Errors without ionospheric combination carrier observations; ρ is the distance (being tried to achieve by receiver rough coordinates and co-ordinates of satellite) of satellite and survey station; Dt is each system receiver clock correction; d _tropit is tropospheric delay; N _iFwithout ionospheric combination blur level; ε _iFmultipath Errors, observed reading noise and other residual errors.Location observation equation Satellite precise clock correction resolves gained positioning error corrected value with all the other and eliminates.

In described step 5 after the linearization of multisystem dynamic precision One-Point Location observation equation, can be expressed as:

V＝AX-L

In above formula, Observation Design matrix A can be expressed as:

$A=\left[\begin{array}{cccccccc}{a}_{x_1\mathrm{GPS}}& a_{y_1\mathrm{GPS}}& a_{z_1\mathrm{GPS}}& 1& 0& 0& M_{1\mathrm{GPS}}& ...\\ ...& ...& ...& ...& ...& ...& ...& ...\\ a_{x_{2n}\mathrm{GPS}}& a_{y_{2n}\mathrm{GPS}}& a_{z_{2n}\mathrm{GPS}}& 1& 0& 0& M_{2\mathrm{nGPS}}& ...\\ a_{x_1\mathrm{GLO}}& a_{y_1\mathrm{GLO}}& a_{z_1\mathrm{GLO}}& 0& 1& 0& M_{1\mathrm{GLO}}& ...\\ ...& ...& ...& ...& ...& ...& ...& ...\\ a_{x_{2j}\mathrm{GPS}}& a_{y_{2j}\mathrm{GPS}}& a_{z_{2j}\mathrm{GPS}}& 0& 1& 0& M_{2\mathrm{jGLO}}& ...\\ a_{x_1\mathrm{BDS}}& a_{y_1\mathrm{BDS}}& a_{z_1\mathrm{BDS}}& 0& 0& 1& M_{1\mathrm{BDS}}& ...\\ ...& ...& ...& ...& ...& ...& ...& ...\\ a_{x_{2k}\mathrm{BDS}}& a_{z_{2k}\mathrm{BDS}}& a_{z_{2k}\mathrm{BDS}}& 0& 0& 1& M_{2\mathrm{kBDS}}& ...\end{array}\right]$

State vector can be expressed as:

$\begin{array}{c}X=\\ \left[\begin{array}{cccccccccc}\mathrm{\δx}& \mathrm{\δy}& \mathrm{\δz}& \mathrm{cli}\·{\mathrm{dt}}^g& \mathrm{cli}\·{\mathrm{dt}}^r& \mathrm{cli}\·{\mathrm{dt}}^c& d_{\mathrm{trop}}& N_{\mathrm{IF}}^{g1}...N_{\mathrm{IF}}^{\mathrm{gn}}& N_{\mathrm{IF}}^{r1}...N_{\mathrm{IF}}^{\mathrm{rj}}& N_{\mathrm{IF}}^{c1}...N_{\mathrm{IF}}^{\mathrm{ck}}\end{array}\right]\end{array}$

In formula: V represents observed reading residual vector; A represents Observation Design matrix; X represents state vector; L represents observation vector; Subscript g, r, c represent GPS, GLONASS, BDS tri-system; N, j, k represent the GPS observed, GLONASS, BDS number of satellites; A battle array first three be classified as δ x, δ y after observation equation linearization, the coefficient of δ z; A battle array the 7th is classified as parameter d _tropcoefficient; Other after A battle array the 7th arranges are classified as the coefficient of blur level parameter, and coefficient is set to 1; δ x, δ y, δ z represents x respectively, and y, z defend in direction distance residual error; Cli represents the light velocity; Dt represents the receiver clock-offsets that three systems are corresponding; d _troprepresent tropospheric delay; N _iFindicate without ionosphere integer ambiguity.

Realize the dynamic positioning method of high precision and high stability in described step 6, comprise the following steps:

Step 6.1, according to the correlation parameter obtained in PPP observation equation, sets up dynamic Kalman filter model, and arranges the initial value of filtering parameter; The function model of Kalman filter is:

$\left\{\begin{array}{c}{X}_{k+1}={\mathrm{\Φ}}_k{X}_k+W_k\\ L_{k+1}=A_{k+1}{X}_{k+1}+V_{k+1}\end{array}\right.$

In formula: subscript k represents sequence number epoch; X _krepresent t _kstate vector in moment observation equation; X _k+1represent t _k+1state vector in moment observation equation; Φ _krepresent state-transition matrix, generally get unit matrix; W _krepresent state model input noise vector; L _k+1represent observation vector, be divided into pseudorange and carrier observations; A _k+1represent the matrix of coefficients in observation equation; V _k+1represent observed reading residual vector.

The statistical model of Kalman filter is:

$\left\{\begin{array}{c}E\left(W_k\right)=\stackrel{\→}{0},E\left(V_{k+1}\right)=\stackrel{\→}{0}\\ \mathrm{Cov}\left(W_k\right)=Q_k,\mathrm{Cov}\left(V_{k+1}\right)=R_{k+1},\mathrm{Cov}(W_k,V_{k+1})=0\end{array}\right.$

In formula: mathematical expectation is asked in E () expression; Covariance matrix is asked in Cov () expression; W _krepresent state model input noise vector; V _k+1represent observed reading residual vector; Q _krepresent the covariance matrix of state model input noise vector; R _k+1represent the covariance matrix of observed reading residual vector;

Step 6.2, arranges initial random model according to the experience weights ratio between elevation of satellite and system, the residual vector information of recycling observed reading, adopts Robust filter model to redefine weights:

${\stackrel{\&OverBar;}{p}}_i=\left\{\begin{array}{c}{p}_i,\left|{\tilde{v}}_i\right|\&le;k_0\\ p_i\frac{k_0}{\left|{\tilde{v}}_i\right|}{\left(\frac{k_1-\left|{\tilde{v}}_i\right|}{k_1-k_0}\right)}^2,k_0<\left|{\tilde{v}}_i\right|\&le;k_1\\ 0,\left|{\tilde{v}}_i\right|>k_1\end{array}\right.$

In formula: subscript i represents residual values number; p _irepresent priori weights, obtained by elevation angle weighting formula; represent the residual values of observed reading, i.e. observed reading residual vector V _kin element; k ₀, k ₁represent Robust filter threshold value, utilize the population distribution of residual values to determine; represent the weights after Robust filter adjustment;

According to the weights determined after Robust filter, then the state parameter robust solution that can obtain kth epoch is:

${\tilde{X}}_k={\left(A_k^T{\stackrel{\&OverBar;}{P}}_k{A}_k\right)}^{-1}{A}_k^T{\stackrel{\&OverBar;}{P}}_k{L}_k$

In formula: subscript k represents sequence number epoch; Subscript T representing matrix transposition; represent robust state vector; represent the power battle array after Robust filter adjustment; L _krepresent observation vector; A _krepresent the matrix of coefficients in observation equation.

Step 6.3, in adaptable Kalman filter, needs the Current observation information according to obtaining in filtering and predicted state information, determines that the value of adaptive factor is:

${\mathrm{\&alpha;}}_k=\left\{\begin{array}{c}1,\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|\&le;c_0\\ \frac{c_0}{\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|}{\left(\frac{c_1-\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|}{c_1-c_0}\right)}^2,c_0<\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|\&le;c_1\\ 0,\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|>c_1\end{array}\right.$

In formula:

$\mathrm{\&Delta;}{\tilde{X}}_k=\left|\right|{\tilde{X}}_k-{\stackrel{\&OverBar;}{X}}_k\left|\right|/\sqrt{\mathrm{tr}\left({\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}\right)}$

In formula: subscript k represents sequence number epoch; || || represent and ask norm; represent robust state vector; represent status prediction information vector; c ₀, c ₁represent the value threshold value of adaptive factor; represent the mark of state forecast vector covariance matrix; α _krepresent adaptive factor;

Step 6.4, according to the data that previous step obtains, upgrades observation information:

$\left\{\begin{array}{c}{K}_k=\frac{1}{{\mathrm{\&alpha;}}_k}{\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}{A}_k^T{(\frac{1}{{\mathrm{\&alpha;}}_k}{A}_k{\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}{A}_k^T+{\stackrel{\&OverBar;}{P}}_k^{-1})}^{-1}\\ {\hat{X}}_k={\stackrel{\&OverBar;}{X}}_k+K_k(L_k-A_k{\stackrel{\&OverBar;}{X}}_k)\\ {\mathrm{\&Sigma;}}_{{\hat{X}}_k}=(I-K_k{A}_k){\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}/{\mathrm{\&alpha;}}_k\end{array}\right.$

In formula: subscript T representing matrix transposition; K _krepresent gain matrix; α _krepresent adaptive factor; represent state forecast vector covariance matrix; A _krepresent the matrix of coefficients in observation equation; represent the power battle array after Robust filter adjustment; represent state estimation vector; represent state forecast vector; L _krepresent observation vector; represent state estimation vector covariance matrix; I representation unit matrix;

Step 6.5, utilizes the weights and adaptive factor determined in step 6.2 and step 6.3, in step 6.4, upgrades observation information, carries out filtering and resolves, thus realizes weakening rough error and dynamic noise to the impact of state parameter valuation.

Beneficial effect: the dynamic PPP calculation method of a kind of multisystem based on robust adaptable Kalman filter provided by the invention, on the basis that usable satellite quantity is enough sufficient, the residual information making full use of observed reading redefines probabilistic model, rejects the satellite that there is rough error in dynamic test; The information of the predicted state amount obtained according to Kalman filter, determines adaptive factor, reduces the impact of the larger information of forecasting of error.And adopt pseudorange One-Point Location simple epoch solution receiver general location and clock correction, can not transmit between epoch and accumulate positioning error, effectively prevent the phenomenon because local observational error causes positioning result to be dispersed.Compared with the dynamic PPP calculation method of routine, the method that this patent proposes can increase considerably usable satellite number, effectively weaken the impact of rough error, and improve the dynamic model noise abnormal conditions in Kinematic Positioning, finally reach the Kinematic Positioning result of high precision and high stability.

Accompanying drawing explanation

Fig. 1 is the multisystem dynamic PPP calculation method process flow diagram based on robust adaptable Kalman filter;

Fig. 2 is experiment Tianjin used CORS base station distribution plan;

Fig. 3 is base station 1:BC station GPS, GLONASS, BDS system usable satellite number;

Fig. 4 is base station 2:NH station GPS, GLONASS, BDS system usable satellite number;

Fig. 5 is the dynamic PPP positioning result of base station 1:BC station GPS, GLONASS, BDS tri-system;

Fig. 6 is the dynamic PPP positioning result of base station 2:NH station GPS, GLONASS, BDS tri-system;

Fig. 7 is dynamic test data: the dynamic PPP plane of the multisystem based on robust adaptable Kalman filter result;

Fig. 8 is dynamic test data: the dynamic PPP elevation of the multisystem based on robust adaptable Kalman filter result;

Fig. 9 is dynamic test data: based on the multisystem dynamic PPP elevation result of robust self-adaptation and general Kalman filter in motion process.

Embodiment

Below in conjunction with accompanying drawing, the present invention is further described.

The dynamic PPP calculation method of multisystem based on robust adaptable Kalman filter, as shown in Figure 1: first use three system broadcasts ephemeris computation co-ordinates of satellite and satellite clock correction separately, carries out space-time datum unification to it; In conjunction with observation file in Pseudo-range Observations combination gained without ionosphere Pseudo-range Observations, carry out Iterated adjustment pseudorange One-Point Location, inverse receiver epoch rough coordinates and each system receiver clock correction; Analyze from IGS etc. High Precision Satellite Ephemeris and the precise clock correction product that center website obtains each system, calculate each system-satellite precision coordinate and satellite precise clock correction by lagrange-interpolation; Carry out each Error Correction Model calculation of position errors corrected value subsequently, observation file, satellite precise coordinate and satellite precise clock correction also carries out the quality control of observation data in conjunction with observation data; Finally set up the dynamic PPP location model of multisystem, adopt robust adaptable Kalman filter to calculate residual values and the information of forecasting of observed reading information, determine weight of observation battle array and adaptive factor, realize the dynamic PPP location of multisystem of high precision, high stability.

Specifically comprise the following steps:

Step 1, uses three system broadcasts ephemeris computation co-ordinates of satellite and satellite clock correction separately, carries out space-time datum unification to it, obtains the co-ordinates of satellite after unified benchmark and satellite clock correction; According to P1 (C1) and P2 (C2) Pseudo-range Observations of observation file, carry out obtaining without ionospheric combination Pseudo-range Observations without ionospheric combination; According to L1 and the L2 observed reading of observation file, carry out obtaining without ionospheric combination carrier observations without ionospheric combination.

1), each navigation positioning system has respective time system, by the difference between each system time reference, by unified for each system time reference to GPST, thus eliminates the time deviation between each system.The time system information of GPS, GLONASS, BDS lists in table 1:

Table 1 GPS, GLONASS, BDS time system

In conjunction with the relation between GPST, GLONASST and BDT and UTC time, the time transforming relationship between three can be expressed as:

GPST＝GLONASST+1 ^s×n-19 ^s-3 ^h(1.1)

GPST＝BDT+14 ^s(1.2)

In formula, when GPST represents GPS; When GLONASST represents GLONASS; N represents jump second number; When BDST represents the Big Dipper.

Because the receiver clock-offsets value of each system also exists deviation, when carrying out Iterated adjustment pseudorange One-Point Location, receiver clock-offsets corresponding for each system need be resolved as solve for parameter; Resolve gained each system receiver clock correction substitute in follow-up Static Precise Point Positioning filtering equations as initial value participate in resolve.

For the conversion between GPS, GLONASS, BDS tri-coordinate systems, need logical zeroaxial translation, the conversion of the rotation of coordinate axis and the unit of coordinate axis scale realizes.Two arbitrary three dimensions rectangular coordinate system O-XYZ and O '-X ' Y ' Z ', when they exist three and above known point, can adopt boolean Sha seven parameter model to calculate and be converted to:

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}\mathrm{\&delta;x}\\ \mathrm{\&delta;y}\\ \mathrm{\&delta;z}\end{array}\right]+(1+m)\left[\begin{array}{ccc}1& {\mathrm{\&epsiv;}}_z& -{\mathrm{\&epsiv;}}_y\\ -{\mathrm{\&epsiv;}}_z& 1& {\mathrm{\&epsiv;}}_x\\ {\mathrm{\&epsiv;}}_y& -{\mathrm{\&epsiv;}}_x& 1\end{array}\right]\left[\begin{array}{c}{X}^{\&prime;}\\ Y^{\&prime;}\\ Z^{\&prime;}\end{array}\right]---\left(1.3\right)$

In formula, X, Y, Z and X ', Y ', Z ' represent three-dimensional coordinate under two three dimensions rectangular coordinate systems respectively; δ x, δ y, δ z represent the translation parameters between two coordinate origins respectively; ε _x, ε _y, ε _zdenotation coordination system rotates the rotation parameter produced; M represents dimensional variation parameter.

By long-term calculating and test, the WGS-84 coordinate system of the PZ-90 coordinate system that existing a lot of mechanism and scholar adopt GLONASS and GPS has carried out unified study on the transformation, and at present, the experience conversion formula that universally acknowledged precision is the highest is:

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{\mathrm{WGS}-84}=\left[\begin{array}{c}-0.47\\ -0.51\\ -1.56\end{array}\right]+(1+22\&times;{10}^{-9})\left[\begin{array}{ccc}1& -1.728\&times;{10}^{-6}& -0.017\&times;{10}^{-6}\\ 1.728\&times;{10}^{-6}& 1& 0.076\&times;{10}^{-6}\\ 0.017\&times;{10}^{-6}& -0.076\&times;{10}^{-6}& 1\end{array}\right]{\left[\begin{array}{c}U\\ V\\ W\end{array}\right]}_{\mathrm{PZ}-90}---\left(1.4\right)$

In formula, X, Y, Z represent the three-dimensional coordinate under WGS84 coordinate system; U, V, W represent the three-dimensional coordinate under PZ-90 coordinate system.

At present, between the CGCS2000 coordinate system that BDS is adopted and WGS-84, there is no clear and definite conversion formula, and because coordinate system parameters difference is between the two very little, therefore BDS and gps system are without the need to carrying out coordinate conversion operation.

2), utilize the satellite orbit coordinate in precise ephemeris and double frequency pseudorange observation value information, by " iono-free combination " model, eliminate the single order item of ionosphere delay error, calculate receiver general location and the receiver clock-offsets of meter accuracy.

$P_{\mathrm{IF}}=\frac{f_1^2{P}_1-f_2^2{P}_2}{f_1^2-f_2^2}---\left(1.6\right)$

In above formula, P _iFindicate without ionospheric combination Pseudo-range Observations; f ₁, f ₂represent the carrier observations frequency of GPS/GLONASS/BDS system; P ₁, P ₂represent the Pseudo-range Observations of GPS/GLONASS/BDS system; If there is no P ₁/ P ₂during observed reading, C can be adopted ₁/ C ₂observed reading carries out using after DCB corrects.In multisystem pseudorange One-Point Location equation, comprise 3 receiver location parameters, 3 receiver clock-offsets parameters.

Step 2, what obtain according to step 1 carries out Iterated adjustment pseudorange One-Point Location without the co-ordinates of satellite after ionospheric combination Pseudo-range Observations and unified benchmark and satellite clock correction, obtains receiver rough coordinates and receiver clock-offsets.

During due to Static Precise Point Positioning, only a receiver need be adopted to carry out data acquisition, cannot eliminate as Differential positioning or weaken fractional error.Therefore, need the accurately every Correction of Errors value of process, according to receiver rough coordinates and satellite exact position, error is divided three classes and processes: be relevant to satellite, to propagate with signal and to be correlated with and relevant with receiver.The conventional model of Static Precise Point Positioning is such as formula shown in (1.7) and (1.8):

$\begin{array}{c}{P}_{\mathrm{IF}}=\frac{f_1^2\&CenterDot;P_1-f_2^2\&CenterDot;P_2}{f_1^2-f_2^2}\\ =\mathrm{\&rho;}+\mathrm{cli}\&CenterDot;({\mathrm{dt}}_r-{\mathrm{dt}}^s)+\mathrm{dtrop}+{\mathrm{\&epsiv;}}_{P_{\mathrm{IF}}}\end{array}---\left(1.7\right)$

$\begin{array}{c}{L}_{\mathrm{IF}}=\frac{f_1^2\&CenterDot;L_1-f_2^2\&CenterDot;L_2}{f_1^2-f_2^2}\\ =\mathrm{\&rho;}+\mathrm{cli}\&CenterDot;({\mathrm{dt}}_r-{\mathrm{dt}}^s)+\mathrm{dtrop}+\frac{f_1^2\&CenterDot;{\mathrm{\&lambda;}}_1\&CenterDot;N_1-f_2^2\&CenterDot;{\mathrm{\&lambda;}}_2\&CenterDot;N_2}{f_1^2-f_2^2}+{\mathrm{\&epsiv;}}_{L_{\mathrm{IF}}}\end{array}---\left(1.8\right)$

In above formula, wherein, P _iFindicate without ionospheric combination Pseudo-range Observations; f ₁, f ₂represent the carrier observations frequency of GPS/GLONASS/BDS system; P ₁, P ₂represent the Pseudo-range Observations of GPS/GLONASS/BDS system; If there is no P ₁/ P ₂during observed reading, C can be adopted ₁/ C ₂observed reading carries out using after DCB corrects; ρ is the distance of satellite and survey station; Cli represents the light velocity; Dt _rfor representing receiver clock-offsets; Dt ^srepresent satellite clock correction; d _tropit is tropospheric delay; L _iFindicate without ionospheric combination carrier observations; L ₁, L ₂represent the carrier phase observation data of GPS/GLONASS/BDS system; N ₁, N ₂represent respectively GPS/GLONASS/BDS system without ionosphere blur level; λ ₁, λ _{, 2}represent the carrier phase wavelength of GPS/GLONASS/BDS system; represent respectively pseudorange and carrier wave without ionospheric combination observation noise.If observe n satellite, then observation equation number is 2n, and unknown number has, and (7+n) is individual, comprises 3 location parameters, 3 receiver clock-offsets, 1 troposphere wet stack emission error and n blur level parameters.

Step 3, obtains three system-satellite precise ephemerises and satellite precise clock correction by network from analytic centres such as IGS, according to the epoch time of observation file, carries out satellite precise coordinate and satellite precise clock correction that Lagrange's interpolation obtains the corresponding moment.

IGS website is downloaded precise ephemeris and precise clock correction file, utilizes Lagrange's interpolation to obtain the satellite information of high sampling rate, shown in (1.1):

$P_n\left(x\right)=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{l}_i\left(x\right)\&CenterDot;y_i=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}(y_i\&CenterDot;\underset{j=0,j\&NotEqual;i}{\overset{n}{\mathrm{\&Pi;}}}\frac{x-x_j}{x_i-x_j})---\left(1.5\right)$

In formula, represent the satellite information in interpolation moment, i.e. the satellite orbit coordinate of each epoch and clock correction value after interpolation, being the interpolation moment, is the node corresponding moment; The satellite orbit coordinate corresponding for node time instance and clock correction value; Be expressed as n interpolation basic differential polynomial; P _nx () represents the satellite information in interpolation moment, i.e. the satellite orbit coordinate of each epoch and clock correction value after interpolation.

When adopting Lagrange interpolation polynomial, should choose suitable exponent number, in addition, Lagrange's interpolation extrapolation accuracy is lower, therefore for single day 24 hours observation files, need carry out interpolation in conjunction with the front and back High Precision Satellite Ephemeris of two days and satellite precise clock correction file.

Step 4, the satellite precise coordinate that the receiver rough coordinates that the receiver information, the step 2 that provide according to observation file obtain and receiver clock-offsets and step 3 obtain and precise clock correction, calculated every Correction of Errors value of Static Precise Point Positioning process by Error Correction Model, and carry out the quality control of observation data data in conjunction with observation data.

Step 5, according to step 1 obtain without ionospheric combination Pseudo-range Observations with without ionospheric combination carrier observations, the receiver rough coordinates that step 2 obtains and receiver clock-offsets, the every Correction of Errors value of calculating gained that the satellite precise coordinate that step 3 obtains and satellite precise clock correction and step 4 obtain, in conjunction with Static Precise Point Positioning model, build multisystem dynamic precision One-Point Location equation.

After adopting High Precision Satellite Ephemeris and satellite precise clock correction product, multisystem dynamic precision One-Point Location observation equation can be expressed as:

$P_{\mathrm{IF}}^g={\mathrm{\&rho;}}^g+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^g+d_{\mathrm{trop}}^g+{\mathrm{\&epsiv;}}_{P_{\mathrm{IF}}}^g$

${\mathrm{\&Phi;}}_{\mathrm{IF}}^g={\mathrm{\&rho;}}^g+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^g+d_{\mathrm{trop}}^g+N_{\mathrm{IF}}^g+{\mathrm{\&epsiv;}}_{{\mathrm{\&Phi;}}_{\mathrm{IF}}}^{\&prime;g}$

$P_{\mathrm{IF}}^r={\mathrm{\&rho;}}^r+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^r+d_{\mathrm{trop}}^r+{\mathrm{\&epsiv;}}_{P_{\mathrm{IF}}}^r$

${\mathrm{\&Phi;}}_{\mathrm{IF}}^r={\mathrm{\&rho;}}^r+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^r+d_{\mathrm{trop}}^r+N_{\mathrm{IF}}^r+{\mathrm{\&epsiv;}}_{{\mathrm{\&Phi;}}_{\mathrm{IF}}}^r$

$P_{\mathrm{IF}}^c={\mathrm{\&rho;}}^c+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^c+d_{\mathrm{trop}}^g+{\mathrm{\&epsiv;}}_{P_{\mathrm{IF}}}^c$

${\mathrm{\&Phi;}}_{\mathrm{IF}}^c={\mathrm{\&rho;}}^c+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^c+d_{\mathrm{trop}}^g+N_{\mathrm{IF}}^c+{\mathrm{\&epsiv;}}_{{\mathrm{\&Phi;}}_{\mathrm{IF}}}^c$

In formula, subscript g, r, c represent GPS, GLONASS and BDS system respectively; Cli is the light velocity; P _iFbe after Correction of Errors without ionospheric combination Pseudo-range Observations; Φ _iFbe after Correction of Errors without ionospheric combination carrier observations; ρ is the distance of satellite and survey station; Dt is each system receiver clock correction; d _tropit is tropospheric delay; N _iFwithout ionospheric combination blur level; ε _iFmultipath Errors, observed reading noise and other residual errors.Wherein ρ is tried to achieve by receiver rough coordinates and co-ordinates of satellite; Satellite precise clock correction is eliminated in the lump with all the other positioning error corrected values.

In described step 6 after the linearization of multisystem dynamic precision One-Point Location observation equation, can be expressed as:

V＝AX-L

In above formula, Observation Design matrix A can be expressed as:

$A=\left[\begin{array}{cccccccc}{a}_{x_1\mathrm{GPS}}& a_{y_1\mathrm{GPS}}& a_{z_1\mathrm{GPS}}& 1& 0& 0& M_{1\mathrm{GPS}}& ...\\ ...& ...& ...& ...& ...& ...& ...& ...\\ a_{x_{2n}\mathrm{GPS}}& a_{y_{2n}\mathrm{GPS}}& a_{z_{2n}\mathrm{GPS}}& 1& 0& 0& M_{2\mathrm{nGPS}}& ...\\ a_{x_1\mathrm{GLO}}& a_{y_1\mathrm{GLO}}& a_{z_1\mathrm{GLO}}& 0& 1& 0& M_{1\mathrm{GLO}}& ...\\ ...& ...& ...& ...& ...& ...& ...& ...\\ a_{x_{2j}\mathrm{GPS}}& a_{y_{2j}\mathrm{GPS}}& a_{z_{2j}\mathrm{GPS}}& 0& 1& 0& M_{2\mathrm{jGLO}}& ...\\ a_{x_1\mathrm{BDS}}& a_{y_1\mathrm{BDS}}& a_{z_1\mathrm{BDS}}& 0& 0& 1& M_{1\mathrm{BDS}}& ...\\ ...& ...& ...& ...& ...& ...& ...& ...\\ a_{x_{2k}\mathrm{BDS}}& a_{z_{2k}\mathrm{BDS}}& a_{z_{2k}\mathrm{BDS}}& 0& 0& 1& M_{2\mathrm{kBDS}}& ...\end{array}\right]$

State vector can be expressed as:

$\begin{array}{c}X=\\ \left[\begin{array}{cccccccccc}\mathrm{\&delta;x}& \mathrm{\&delta;y}& \mathrm{\&delta;z}& \mathrm{cli}\&CenterDot;{\mathrm{dt}}^g& \mathrm{cli}\&CenterDot;{\mathrm{dt}}^r& \mathrm{cli}\&CenterDot;{\mathrm{dt}}^c& d_{\mathrm{trop}}& N_{\mathrm{IF}}^{g1}...N_{\mathrm{IF}}^{\mathrm{gn}}& N_{\mathrm{IF}}^{r1}...N_{\mathrm{IF}}^{\mathrm{rj}}& N_{\mathrm{IF}}^{c1}...N_{\mathrm{IF}}^{\mathrm{ck}}\end{array}\right]\end{array}$

In formula: V represents observed reading residual vector; A represents Observation Design matrix; X represents state vector; L represents observation vector; Subscript g, r, c represent GPS, GLONASS, BDS tri-system; N, j, k represent the GPS observed, GLONASS, BDS number of satellites; A battle array first three be classified as δ x, δ y after observation equation linearization, the coefficient of δ z; A battle array the 7th is classified as parameter d _tropcoefficient; Other after A battle array the 7th arranges are classified as the coefficient of blur level parameter, and coefficient is set to 1; δ x, δ y, δ z represents x respectively, and y, z defend in direction distance residual error; Cli represents the light velocity; Dt represents the receiver clock-offsets that three systems are corresponding; d _troprepresent tropospheric delay; N _iFindicate without ionosphere integer ambiguity.

Factor arrays for any satellite i is:

$A=\left[\begin{array}{ccccccc}{a}_i& b_i& c_i& 1& M_i& 1& 0\\ a_i& b_i& c_i& 1& M_i& 0& 0\end{array}\right]---\left(1.11\right)$

In formula: x _i, y _i, z _irepresent the coordinate of satellite i, X ^r, Y ^r, Z ^rrepresent the coordinate of receiver r, represent the distance of satellite to receiver; M _irepresent the troposphere wet stack emission projection function on zenith direction; The first row is carrier observations equation coefficients, and the second row is pseudorange observation equation coefficients.

Step 7, adopt the multisystem dynamic precision One-Point Location equation that robust adaptive Kalman filter process of solution 6 obtains, according to the status information of the observation information obtained in filtering and prediction, the weights of adjustment observed reading and calculating adaptive factor, realize the Kinematic Positioning of high precision and high stability.

Usually adopt Kalman filter to resolve when Static Precise Point Positioning, reliable Kalman filter requires function model and probabilistic model accurately.In real life, general object is difficult to the motion guaranteeing to do rule, thus builds accurate function model very difficult; Probabilistic model is generally determined according to existing prior imformation, usually adopts the sine value of elevation angle to calculate, also there is certain deviation with actual conditions.So the status information of current observation information and prediction can be utilized to weaken fractional error in motion, namely adopt robust model and auto adapted filtering.

Realize the dynamic positioning method of high precision and high stability in described step 6, comprise the following steps:

Step 6.1, according to the correlation parameter obtained in PPP observation equation, sets up dynamic Kalman filter model, and arranges the initial value of filtering parameter; The function model of Kalman filter is:

$\left\{\begin{array}{c}{X}_{k+1}={\mathrm{\&Phi;}}_k{X}_k+W_k\\ L_{k+1}=A_{k+1}{X}_{k+1}+V_{k+1}\end{array}\right.$

In formula: subscript k represents sequence number epoch; X _krepresent t _kstate vector in moment observation equation; X _k+1represent t _k+1state vector in moment observation equation; Φ _krepresent state-transition matrix, generally get unit matrix; W _krepresent state model input noise vector, in the dynamic PPP of multisystem, the corresponding numerical value of coordinate gets empirical value 100; L _k+1represent observation vector, be divided into pseudorange and carrier observations; A _k+1represent the matrix of coefficients in observation equation; V _k+1represent observed reading residual vector.

The statistical model of Kalman filter is:

$\left\{\begin{array}{c}E\left(W_k\right)=\stackrel{\&RightArrow;}{0},E\left(V_{k+1}\right)=\stackrel{\&RightArrow;}{0}\\ \mathrm{Cov}\left(W_k\right)=Q_k,\mathrm{Cov}\left(V_{k+1}\right)=R_{k+1},\mathrm{Cov}(W_k,V_{k+1})=0\end{array}\right.$

In formula: mathematical expectation is asked in E () expression; Covariance matrix is asked in Cov () expression; W _krepresent state model input noise vector; V _k+1represent observed reading residual vector; Q _krepresent the covariance matrix of state model input noise vector; R _k+1represent the covariance matrix of observed reading residual vector;

Step 6.2, arranges initial random model according to the experience weights ratio between elevation of satellite and system, the residual vector information of recycling observed reading, adopts Robust filter model to redefine weights:

${\stackrel{\&OverBar;}{p}}_i=\left\{\begin{array}{c}{p}_i,\left|{\tilde{v}}_i\right|\&le;k_0\\ p_i\frac{k_0}{\left|{\tilde{v}}_i\right|}{\left(\frac{k_1-\left|{\tilde{v}}_i\right|}{k_1-k_0}\right)}^2,k_0<\left|{\tilde{v}}_i\right|\&le;k_1\\ 0,\left|{\tilde{v}}_i\right|>k_1\end{array}\right.$

In formula: subscript i represents residual values number; p _irepresent priori weights, obtained by elevation angle weighting formula; represent the residual values of observed reading, i.e. observed reading residual vector V _kin element; k ₀, k ₁represent Robust filter threshold value, utilize the population distribution of residual values to determine; represent the weights after Robust filter adjustment;

According to the weights determined after Robust filter, then the state parameter robust solution that can obtain kth epoch is:

${\tilde{X}}_k={\left(A_k^T{\stackrel{\&OverBar;}{P}}_k{A}_k\right)}^{-1}{A}_k^T{\stackrel{\&OverBar;}{P}}_k{L}_k$

In formula: subscript k represents sequence number epoch; Subscript T representing matrix transposition; represent robust state vector; represent the power battle array after Robust filter adjustment; L _krepresent observation vector; A _krepresent the matrix of coefficients in observation equation.

Step 6.3, in adaptable Kalman filter, needs the Current observation information according to obtaining in filtering and predicted state information, determines that the value of adaptive factor is:

${\mathrm{\&alpha;}}_k=\left\{\begin{array}{c}1,\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|\&le;c_0\\ \frac{c_0}{\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|}{\left(\frac{c_1-\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|}{c_1-c_0}\right)}^2,c_0<\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|\&le;c_1\\ 0,\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|>c_1\end{array}\right.$

In formula:

$\mathrm{\&Delta;}{\tilde{X}}_k=\left|\right|{\tilde{X}}_k-{\stackrel{\&OverBar;}{X}}_k\left|\right|/\sqrt{\mathrm{tr}\left({\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}\right)}$

In formula: subscript k represents sequence number epoch; || || represent and ask norm; represent robust state vector; represent status prediction information vector; c ₀, c ₁represent the value threshold value of adaptive factor; represent the mark of state forecast vector covariance matrix; α _krepresent adaptive factor;

Step 6.4, according to the data that previous step obtains, upgrades observation information:

$\left\{\begin{array}{c}{K}_k=\frac{1}{{\mathrm{\&alpha;}}_k}{\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}{A}_k^T{(\frac{1}{{\mathrm{\&alpha;}}_k}{A}_k{\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}{A}_k^T+{\stackrel{\&OverBar;}{P}}_k^{-1})}^{-1}\\ {\hat{X}}_k={\stackrel{\&OverBar;}{X}}_k+K_k(L_k-A_k{\stackrel{\&OverBar;}{X}}_k)\\ {\mathrm{\&Sigma;}}_{{\hat{X}}_k}=(I-K_k{A}_k){\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}/{\mathrm{\&alpha;}}_k\end{array}\right.$

In formula: subscript T representing matrix transposition; K _krepresent gain matrix; α _krepresent adaptive factor; represent state forecast vector covariance matrix; A _krepresent the matrix of coefficients in observation equation; represent the power battle array after Robust filter adjustment; represent state estimation vector; represent state forecast vector; L _krepresent observation vector; represent state estimation vector covariance matrix; I representation unit matrix;

Step 6.5, utilizes the weights and adaptive factor determined in step 6.2 and step 6.3, in step 6.4, upgrades observation information, carries out filtering and resolves, thus realizes weakening rough error and dynamic noise to the impact of state parameter valuation.

The above is only the preferred embodiment of the present invention; be noted that for those skilled in the art; under the premise without departing from the principles of the invention, can also make some improvements and modifications, these improvements and modifications also should be considered as protection scope of the present invention.

Claims (7)

1. based on the dynamic PPP calculation method of multisystem of robust adaptable Kalman filter, it is characterized in that: first use three system broadcasts ephemeris computation co-ordinates of satellite and satellite clock correction separately, space-time datum unification is carried out to it; In conjunction with observation file in Pseudo-range Observations combination gained without ionosphere Pseudo-range Observations, carry out Iterated adjustment pseudorange One-Point Location, inverse receiver epoch rough coordinates and each system receiver clock correction; Analyze from IGS etc. High Precision Satellite Ephemeris and the precise clock correction product that center website obtains each system, calculate each system-satellite precision coordinate and satellite precise clock correction by lagrange-interpolation; Subsequently according to each Error Correction Model calculation of position errors corrected value, and carry out the quality control of observation data in conjunction with observation data; Finally set up the dynamic PPP location model of multisystem, adopt robust adaptable Kalman filter to calculate residual values and the information of forecasting of observed reading information, determine weight of observation battle array and adaptive factor, realize the dynamic PPP location of high precision, high stability.

2. the dynamic PPP calculation method of the multisystem based on robust adaptable Kalman filter according to claim 1, is characterized in that, comprise the steps:

Step 1, uses three system broadcasts ephemeris computation co-ordinates of satellite and satellite clock correction separately, carries out space-time datum unification to it, obtains the co-ordinates of satellite after unified benchmark and satellite clock correction; According to P1 (C1) and P2 (C2) Pseudo-range Observations of observation file, carry out obtaining without ionospheric combination Pseudo-range Observations without ionospheric combination; According to L1 and the L2 observed reading of observation file, carry out obtaining without ionospheric combination carrier observations without ionospheric combination;

Step 2, what obtain according to step 1 carries out Iterated adjustment pseudorange One-Point Location without the co-ordinates of satellite after ionospheric combination Pseudo-range Observations and unified benchmark and satellite clock correction, obtains receiver rough coordinates and each system receiver clock correction;

Step 3, obtains three system-satellite precise ephemerises and satellite precise clock correction by network from analytic centres such as IGS, according to moment epoch of observation file, carries out satellite precise coordinate and satellite precise clock correction that Lagrange's interpolation obtains the corresponding moment;

Step 4, the satellite precise coordinate that the receiver rough coordinates that the receiver information, the step 2 that provide according to observation file obtain and receiver clock-offsets and step 3 obtain and precise clock correction, calculated every Correction of Errors value of Static Precise Point Positioning process by Error Correction Model, and carry out the quality control of observation data data in conjunction with observation data;

Step 5, according to step 1 obtain without ionospheric combination Pseudo-range Observations with without ionospheric combination carrier observations, the receiver rough coordinates that step 2 obtains and receiver clock-offsets, the every Correction of Errors value of calculating gained that the satellite precise coordinate that step 3 obtains and satellite precise clock correction and step 4 obtain, in conjunction with Static Precise Point Positioning model, build multisystem dynamic precision One-Point Location equation;

Step 6, adopt the multisystem dynamic precision One-Point Location equation that robust adaptive Kalman filter process of solution 5 obtains, according to the status information of the observation information obtained in filtering and prediction, the weights of adjustment observed reading and calculating adaptive factor, realize the Kinematic Positioning of high precision and high stability.

3. the dynamic PPP calculation method of the multisystem based on robust adaptable Kalman filter according to claim 2, is characterized in that: obtain without ionospheric combination Pseudo-range Observations and the formula without ionospheric combination carrier observations without ionospheric combination in described step 1:

$\begin{array}{c}{P}_{\mathrm{IF}}=\frac{{f_1}^2\&CenterDot;P_1-{f_2}^2\&CenterDot;P_2}{{f_1}^2-{f_2}^2}\\ =\mathrm{\&rho;}+\mathrm{cli}\&CenterDot;({\mathrm{dt}}_r-{\mathrm{dt}}^s)+\mathrm{dtrop}+{\mathrm{\&epsiv;}}_{P_{\mathrm{IF}}}\end{array}$

$\begin{array}{c}{L}_{\mathrm{IF}}=\frac{{f_1}^2\&CenterDot;L_1-{f_2}^2\&CenterDot;L_2}{{f_1}^2-{f_2}^2}\\ =\mathrm{\&rho;}+\mathrm{cli}\&CenterDot;({\mathrm{dt}}_r-{\mathrm{dt}}^s)+\mathrm{dtrop}+\frac{{f_1}^2\&CenterDot;{\mathrm{\&lambda;}}_1\&CenterDot;N_1-{f_2}^2\&CenterDot;{\mathrm{\&lambda;}}_2\&CenterDot;N_2}{{f_1}^2-{f_2}^2}+{\mathrm{\&epsiv;}}_{L_{\mathrm{IF}}}\end{array}$

Wherein, P _iFindicate without ionospheric combination Pseudo-range Observations; f ₁, f ₂represent the carrier observations frequency of GPS/GLONASS/BDS system; P ₁, P ₂represent the Pseudo-range Observations of GPS/GLONASS/BDS system; ρ is the distance of satellite and survey station; Cli represents the light velocity; Dt _rfor representing receiver clock-offsets; Dt ^srepresent satellite clock correction; d _tropit is tropospheric delay; L _iFindicate without ionospheric combination carrier observations; L ₁, L ₂represent the carrier phase observation data of GPS/GLONASS/BDS system; λ ₁, λ ₂represent the carrier phase wavelength of GPS/GLONASS/BDS system respectively, N ₁, N ₂represent respectively GPS/GLONASS/BDS system without ionosphere blur level; represent that pseudorange and carrier wave are without ionospheric combination observation noise respectively.

4. the dynamic PPP calculation method of the multisystem based on robust adaptable Kalman filter according to claim 3, is characterized in that: the satellite precise coordinate that described step 2 obtains and satellite precise clock correction:

$P_n\left(x\right)=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{l}_i\left(x\right)\&CenterDot;y_i=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}(y_i\&CenterDot;\underset{j=0,j\&NotEqual;i}{\overset{n}{\mathrm{\&Pi;}}}\frac{x-x_j}{x_i-x_j})$

In formula, P _nx () represents satellite orbit coordinate and the clock correction value of each epoch after interpolation, x is the interpolation moment, x _i(i=0,1 ..., n) be the node corresponding moment; y _i(i=0,1 ..., n) be satellite orbit coordinate corresponding to node time instance and clock correction value; l _i(x) (i=0,1 ..., n) be expressed as n interpolation basic differential polynomial.

5. the dynamic PPP calculation method of the multisystem based on robust adaptable Kalman filter according to claim 4, is characterized in that: the multisystem dynamic precision One-Point Location observation equation that described step 5 obtains is:

$P_{\mathrm{IF}}^g={\mathrm{\&rho;}}^g+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^g+d_{\mathrm{trop}}^g+{\mathrm{\&epsiv;}}_{P_{\mathrm{IF}}}^g$

${\mathrm{\&Phi;}}_{\mathrm{IF}}^g={\mathrm{\&rho;}}^g+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^g+d_{\mathrm{trop}}^g+N_{\mathrm{IF}}^g+{\mathrm{\&epsiv;}}_{{\mathrm{\&Phi;}}_{\mathrm{IF}}}^{\&prime;g}$

$P_{\mathrm{IF}}^r={\mathrm{\&rho;}}^r+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^r+d_{\mathrm{trop}}^r+{\mathrm{\&epsiv;}}_{P_{\mathrm{IF}}}^r$

${\mathrm{\&Phi;}}_{\mathrm{IF}}^r={\mathrm{\&rho;}}^r+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^r+d_{\mathrm{trop}}^r+N_{\mathrm{IF}}^r+{\mathrm{\&epsiv;}}_{{\mathrm{\&Phi;}}_{\mathrm{IF}}}^r$

$P_{\mathrm{IF}}^c={\mathrm{\&rho;}}^c+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^c+d_{\mathrm{trop}}^g+{\mathrm{\&epsiv;}}_{P_{\mathrm{IF}}}^c$

${\mathrm{\&Phi;}}_{\mathrm{IF}}^c={\mathrm{\&rho;}}^c+\mathrm{cli}\&CenterDot;{\mathrm{dt}}^c+d_{\mathrm{trop}}^g+N_{\mathrm{IF}}^c+{\mathrm{\&epsiv;}}_{{\mathrm{\&Phi;}}_{\mathrm{IF}}}^c$

In formula, subscript g, r, c represent GPS, GLONASS and BDS system respectively; Cli is the light velocity; P _iFbe after Correction of Errors without ionospheric combination Pseudo-range Observations; Φ _iFbe after Correction of Errors without ionospheric combination carrier observations; ρ is the distance of satellite and survey station; Dt is each system receiver clock correction; d _tropit is tropospheric delay; N _iFwithout ionospheric combination blur level; ε _iFmultipath Errors, observed reading noise and other residual errors; Location observation equation Satellite precise clock correction resolves gained positioning error corrected value with all the other and eliminates.

6. the dynamic PPP calculation method of the multisystem based on robust adaptable Kalman filter according to claim 5, is characterized in that: in described step 5 after the linearization of multisystem dynamic precision One-Point Location observation equation, can be expressed as:

V＝AX-L

In above formula, Observation Design matrix A can be expressed as:

$A=\left[\begin{array}{cccccccc}{a}_{x_1\mathrm{GPS}}& a_{y_1\mathrm{GPS}}& a_{z_1\mathrm{GPS}}& 1& 0& 0& M_{1\mathrm{GPS}}& \&CenterDot;\&CenterDot;\&CenterDot;\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ a_{x_{2n}\mathrm{GPS}}& a_{y_{2n}\mathrm{GPS}}& a_{z_{2n}\mathrm{GPS}}& 1& 0& 0& M_{2\mathrm{nGPS}}& \&CenterDot;\&CenterDot;\&CenterDot;\\ a_{x_1\mathrm{GLO}}& a_{y_1\mathrm{GLO}}& a_{z_1\mathrm{GLO}}& 0& 1& 0& M_{1\mathrm{GLO}}& \&CenterDot;\&CenterDot;\&CenterDot;\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ a_{x_{2j}\mathrm{GPS}}& a_{y_{2j}\mathrm{GPS}}& a_{z_{2j}\mathrm{GPS}}& 0& 1& 0& M_{2\mathrm{jGLO}}& \&CenterDot;\&CenterDot;\&CenterDot;\\ a_{x_1\mathrm{BDS}}& a_{y_1\mathrm{BDS}}& a_{z_1\mathrm{BDS}}& 0& 0& 1& M_{1\mathrm{BDS}}& \&CenterDot;\&CenterDot;\&CenterDot;\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ a_{x_{2k}\mathrm{BDS}}& a_{z_{2k}\mathrm{BDS}}& a_{z_{2k}\mathrm{BDS}}& 0& 0& 1& M_{2\mathrm{kBDS}}& \&CenterDot;\&CenterDot;\&CenterDot;\end{array}\right]$

State vector can be expressed as:

$\begin{array}{c}X=\\ \left[\begin{array}{cccccccccc}\mathrm{\&delta;x}& \mathrm{\&delta;y}& \mathrm{\&delta;z}& \mathrm{cli}\&CenterDot;{\mathrm{dt}}^g& \mathrm{cli}\&CenterDot;{\mathrm{dt}}^r& \mathrm{cli}\&CenterDot;{\mathrm{dt}}^c& d_{\mathrm{trop}}& N_{\mathrm{IF}}^{g1}\&CenterDot;\&CenterDot;\&CenterDot;N_{\mathrm{IF}}^{\mathrm{gn}}& N_{\mathrm{IF}}^{r1}\&CenterDot;\&CenterDot;\&CenterDot;N_{\mathrm{IF}}^{\mathrm{rj}}& N_{\mathrm{IF}}^{c1}\&CenterDot;\&CenterDot;\&CenterDot;N_{\mathrm{IF}}^{\mathrm{ck}}\end{array}\right]\end{array}$

In formula: V represents observed reading residual vector; A represents Observation Design matrix; X represents state vector; L represents observation vector; Subscript g, r, c represent GPS, GLONASS, BDS tri-system; N, j, k represent the GPS observed, GLONASS, BDS number of satellites; A battle array first three be classified as δ x, δ y after observation equation linearization, the coefficient of δ z; A battle array the 7th is classified as parameter d _tropcoefficient; Other after A battle array the 7th arranges are classified as the coefficient of blur level parameter, and coefficient is set to 1; δ x, δ y, δ z represents x respectively, and y, z defend in direction distance residual error; Cli represents the light velocity; Dt represents the receiver clock-offsets that three systems are corresponding; d _troprepresent tropospheric delay; N _iFindicate without ionosphere integer ambiguity.

7. the dynamic PPP calculation method of the multisystem based on robust adaptable Kalman filter according to claim 6, is characterized in that: the dynamic positioning method realizing high precision and high stability in described step 6, comprises the following steps:

Step 6.1, according to the correlation parameter obtained in PPP observation equation, sets up dynamic Kalman filter model, and arranges the initial value of filtering parameter; The function model of Kalman filter is:

$\left\{\begin{array}{c}{X}_{k+1}={\mathrm{\&Phi;}}_k{X}_k+W_k\\ L_{k+1}=A_{k+1}{X}_{k+1}+V_{k+1}\end{array}\right.$

In formula: subscript k represents sequence number epoch; X _krepresent t _kstate vector in moment observation equation; X _k+1represent t _k+1state vector in moment observation equation; Φ _krepresent state-transition matrix, generally get unit matrix; W _krepresent state model input noise vector; L _k+1represent observation vector, be divided into pseudorange and carrier observations; A _k+1represent the matrix of coefficients in observation equation; V _k+1represent observed reading residual vector.

The statistical model of Kalman filter is:

$\left\{\begin{array}{c}E\left(W_k\right)=\stackrel{\&RightArrow;}{0},E\left(V_{k+1}\right)=\stackrel{\&RightArrow;}{0}\\ \mathrm{Cov}\left(W_k\right)=Q_k,\mathrm{Cov}\left(V_{k+1}\right)=R_{k+1},\mathrm{Cov}(W_k,V_{k+1})=0\end{array}\right.$

In formula: mathematical expectation is asked in E () expression; Covariance matrix is asked in Cov () expression; W _krepresent state model input noise vector; V _k+1represent observed reading residual vector; Q _krepresent the covariance matrix of state model input noise vector; R _k+1represent the covariance matrix of observed reading residual vector;

Step 6.2, arranges initial random model according to the experience weights ratio between elevation of satellite and system, the residual vector information of recycling observed reading, adopts Robust filter model to redefine weights:

${\stackrel{\&OverBar;}{p}}_i=\left\{\begin{array}{c}{p}_i,\left|{\tilde{v}}_i\right|\&le;k_0\\ p_i\frac{k_0}{\left|{\tilde{v}}_i\right|}{\left(\frac{k_1-\left|{\tilde{v}}_i\right|}{k_1-k_0}\right)}^2,k_0<\left|{\tilde{v}}_i\right|\&le;k_1\\ 0,\left|{\tilde{v}}_i\right|>k_1\end{array}\right.$

In formula: subscript i represents residual values number; p _irepresent priori weights; represent the residual values of observed reading; k ₀, k ₁represent Robust filter threshold value; represent the weights after Robust filter adjustment;

According to the weights determined after Robust filter, then the state parameter robust solution that can obtain kth epoch is:

${\tilde{X}}_k={\left(A_k^T{\stackrel{\&OverBar;}{P}}_k{A}_k\right)}^{-1}{A}_k^T{\stackrel{\&OverBar;}{P}}_k{L}_k$

In formula: subscript k represents sequence number epoch; Subscript T representing matrix transposition; represent robust state vector; represent the power battle array after Robust filter adjustment; L _krepresent observation vector; A _krepresent the matrix of coefficients in observation equation.

Step 6.3, in adaptable Kalman filter, needs the Current observation information according to obtaining in filtering and predicted state information, determines that the value of adaptive factor is:

${\mathrm{\&alpha;}}_k=\left\{\begin{array}{c}1,\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|\&le;c_0\\ \frac{c_0}{\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|}{\left(\frac{c_1-\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|}{c_1-c_0}\right)}^2,c_0<\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|\&le;c_1\\ 0,\left|\mathrm{\&Delta;}{\tilde{X}}_k\right|>c_1\end{array}\right.$

In formula:

$\mathrm{\&Delta;}{\tilde{X}}_k=\left|\right|{\tilde{X}}_k-{\stackrel{\&OverBar;}{X}}_k\left|\right|/\sqrt{\mathrm{tr}\left({\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}\right)}$

In formula: subscript k represents sequence number epoch; || || represent and ask norm; represent robust state vector; represent status prediction information vector; c ₀, c ₁represent the value threshold value of adaptive factor; represent the mark of state forecast vector covariance matrix; α _krepresent adaptive factor;

Step 6.4, according to the data that previous step obtains, upgrades observation information:

$\left\{\begin{array}{c}{K}_k=\frac{1}{{\mathrm{\&alpha;}}_k}{\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}{A}_k^T{(\frac{1}{{\mathrm{\&alpha;}}_k}{A}_k{\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}{A}_k^T+{\stackrel{\&OverBar;}{P}}_k^{-1})}^{-1}\\ {\hat{X}}_k={\stackrel{\&OverBar;}{X}}_k+K_k(L_k-A_k{\stackrel{\&OverBar;}{X}}_k)\\ {\mathrm{\&Sigma;}}_{{\hat{X}}_k}=(I-K_k{A}_k){\mathrm{\&Sigma;}}_{{\stackrel{\&OverBar;}{X}}_k}/{\mathrm{\&alpha;}}_k\end{array}\right.$

In formula: subscript T representing matrix transposition; K _krepresent gain matrix; α _krepresent adaptive factor; represent state forecast vector covariance matrix; A _krepresent the matrix of coefficients in observation equation; represent the power battle array after Robust filter adjustment; represent state estimation vector; represent state forecast vector; L _krepresent observation vector; represent state estimation vector covariance matrix; I representation unit matrix;

Step 6.5, utilizes the weights and adaptive factor determined in step 6.2 and step 6.3, in step 6.4, upgrades observation information, carries out filtering and resolves, thus realizes weakening rough error and dynamic noise to the impact of state parameter valuation.