CN105372691A - Long baseline satellite formation GNSS relative positioning method based on ambiguity fixing - Google Patents

Abstract

A long baseline satellite formation GNSS relative positioning method based on ambiguity fixing is provided in order to improve the success rate of ambiguity fixing and the accuracy of relative positioning results. According to the technical scheme, the method comprises the following steps: first, collecting and pre-processing input data, and determining the absolute general orbit of a formation satellite; then, eliminating the geometric distance and clock error in differential observation data, estimating a single-difference phase ambiguity float solution and a single-difference ionosphere delay parameter, carrying out double-difference transform to get a double-difference wide-lane ambiguity float solution and a covariance matrix, and fixing the double-difference wide-lane integer ambiguity and the double-difference narrow-lane integer ambiguity; and finally, outputting the relative positioning result of ambiguity fixing. By adopting the method of the invention, the problem that ambiguity fixing strongly depends on a pseudo code with low observation precision due to equally-weighted pseudo code and phase processing in M-W combination in the traditional method is avoided, the success rate of long baseline satellite formation GNSS relative positioning ambiguity fixing and the accuracy of final relative positioning results are improved, calculation is stable, and the reliability of relative positioning results is improved.

Description

The Long baselines satellites formation GNSS relative positioning method that a kind of blur level is fixing

Technical field

The present invention relates to the relative positioning method of a kind of aerospace measurement field medium-long baselines satellites formation, specifically one have employed the Long baselines satellites formation relative positioning method of GPS (Global Position System) GNSS (GlobalNavigationSatelliteSystem) measurement means and blur level fixed policy.

Background technology

Spaceborne GNSS adopts star-star tracking mode, by installing GNSS receiver on the satellite of formation flight, adopting carrier phase difference GNSS technology, eliminating or slacken the impact of some common error, can provide millimetre-sized positioning precision.Spaceborne GNSS has the advantages such as round-the-clock, the covering of continuity, high precision, space-time is wide, utilizes the relative position of spaceborne GNSS high precision determination Satellite Formation Flying to be successfully applied to the tasks such as distributed radar interference, terrestrial gravitation field measurement.

The key of Long baselines satellites formation GNSS high precision relative positioning is that the high success rate of phase ambiguity is fixed, GNSS two difference phase ambiguity has integer characteristic, if blur level is accurately fixed, blur level parameter then in two poor observation model can be eliminated, now phase observations data can be converted into high precision relative distance, can realize high-precision relative positioning.In Long baselines GNSS blur level is fixing, major effect is from the ionosphere delay in GNSS signal communication process, major part in ionosphere delay is eliminated by difference processing, but ionosphere delay remaining after difference becomes greatly and gradually large along with the change of Satellite Formation Flying spacing.When analyzing GNSS blur level fixed effect, general with carrier phase wavelength X (about 20cm) for reference, when interstellar distance reach tens kilometers even up to a hundred kilometers time, ionosphere delay remaining after difference can reach several decimeters of magnitudes to rice, far more than λ/4, now it can not be left in the basket on the impact that Phase integer ambiguity is fixing.According to the literature, fix about Long baselines satellites formation GNSS relative positioning blur level at present, mainly contain following methods: a kind of method adopts Kuan Xianghezhai lane combined method, see " BerneseGPSSoftware:Version5.0 " that DachR etc. write in 2007 in Long baselines GPS relative positioning between Bernese software ground survey station in early days.First the method uses M-W (Melbourne-W ü bbena) to combine, eliminate the impact of geometric distance, ionosphere and clock correction, solve wide lane integer ambiguity, then construct only have electric eliminating absciss layer carrier phase to observe two differential mode types to solve narrow lane ambiguity.The shortcoming of the method is that in M-W combination, the combination of the power such as pseudo-code, phase place causes wide lane integer ambiguity to be fixedly strongly depend on the lower pseudo-code observation data of accuracy of observation, success ratio is limited, and once wide lane ambiguity cannot be fixed, narrow lane ambiguity can be caused also cannot to fix, reduce GNSS positioning precision.Another kind method is the sequential EKF filtering method that KroesR adopts in GRACE double star relative positioning, see the PhD dissertation " Preciserelativepositioningofformationflyingspacecraftusi ngGPS " of KroesR in 2006.The ionosphere delay random series that the method single order Gauss-Markov model is remaining after carrying out approximate representation difference, adopt sequential EKF filtering method, ionosphere delay is added in filter state parameter to be estimated as variable, while estimation relative orbital parameter, estimate ionosphere delay parameter, realize estimation and the calibration of ionosphere delay.The advantage of the method can differ from integer ambiguity by recurrence estimation pair online, shortcoming is if certain integer ambiguity solid error, near affecting, other integer ambiguity is fixing, positioning precision is caused to worsen, even can cause filtering divergence time serious, cause exporting relative positioning result accurately.

In sum, the precision of Long baselines satellites formation GNSS relative positioning method and reliability still need further improvement.

Summary of the invention

The technical problem to be solved in the present invention is: be fixedly strongly depend on for blur level in traditional Long baselines satellites formation GNSS relative positioning method the not high problem of positioning precision that pseudo-code causes, the Long baselines satellites formation GNSS relative positioning method that a kind of blur level is fixing is proposed, blur level is fixed by the pseudo-code after geometric distance, clock correction in elimination difference observation data, the combination of phase place differential weights value, pseudo-code in avoiding classic method M-W to combine, phase place by etc. power combination, improve blur level and be fixed into the precision of power and relative positioning result.The method adopts batch processing least square processing mode, and calculation stability, overcomes the shortcoming that Sequential filter disposal route is easily dispersed, and improves relative positioning reliability.

Technical scheme of the present invention comprises the following steps:

The first step, collects and pre-service input data.

1.1 collect observation data from satellite platform: GNSS observation data, satellite platform attitude data; Auxiliary data is collected: the change of GNSS satellite precise ephemeris, clock correction and antenna phase center PCV (PhaseCenterVariation) information, earth gravity field, earth rotation parameter, UTC (CoordinatedUniversalTime) time jump second data, JPL (JetPropulsionLaboratory) solar calendar, solar radiation flow, geomagnetic indices from website;

1.2 pairs of observation datas of collecting carry out pre-service: reject the outlier of GNSS observation data, and detecting phase cycle slip.Unruly-value rejecting and Cycle Slips Detection are shown in the monograph " MeasurementDataModelingandPrameterEstimation " (measurement data modeling and parameter estimation) Section 7.2 that ZhengmingWang etc. published in 2011.

Second step, determines the absolute outline track of Satellite Formation Flying.

Non-poor reduced-dynamic method (monograph " MeasurementDataModelingandPrameterEstimation " (measurement data modeling and parameter estimation) Section 7.3 see ZhengmingWang etc. published in 2011) is adopted to determine the absolute orbit position of Satellite Formation Flying A, B.Export the orbital mechanics model parameter y of Satellite Formation Flying A, B _a0, y _b0with the absolute orbit position r of each moment t _a(t), r _bt (), orbital position comprises x, y, z three cartesian component, and precision, at several cm, is mainly used as the outline orbital position of follow-up more high precision relative positioning.

3rd step, eliminates the geometric distance in difference observation data and clock correction.

3.1 single poor relative positionings.Method is: single poor be same GNSS satellite to two receivers observations, same kind observation data does difference, form difference observation data, adopt the poor observation equation of list (1) of t to eliminate GNSS satellite clock correction, retain the relative clock correction of receiver,

${\mathrm{\ΔP}}_1^j\left(t\right)=P_{1,A}^j\left(t\right)-P_{1,B}^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)+I_{AB}^j\left(t\right)+{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_1}^j\left(t\right)$

$\begin{array}{c}{\mathrm{\ΔP}}_2^j\left(t\right)=P_{2,A}^j\left(t\right)-P_{2,B}^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)+\frac{f_1^2}{f_2^2}\·I_{AB}^j\left(t\right)+{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_2}^j\left(t\right)\\ {\mathrm{\ΔL}}_1^j\left(t\right)=L_{1,A}^j\left(t\right)-L_{1,B}^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)-I_{AB}^j\left(t\right)+{\mathrm{\λ}}_1\·{\mathrm{\ΔN}}_1^{\left(k\right)}+{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_1}^j\left(t\right)\end{array}---\left(1\right)$

${\mathrm{\ΔL}}_2^j\left(t\right)=L_{2,A}^j\left(t\right)-L_{2,B}^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)-\frac{f_1^2}{f_2^2}\·I_{AB}^j\left(t\right)+{\mathrm{\λ}}_2\·{\mathrm{\ΔN}}_2^{\left(k\right)}+{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_2}^j\left(t\right)$

Wherein t represents the observation moment, and subscript j (1≤j≤M) represents a jth GNSS satellite, and M is the sum of satellite; Subscript 1,2 represents GNSS signal frequency f ₁, f ₂; Δ represents single poor; for in moment t frequency f ₁pseudo-code difference observation data; for in moment t frequency f ₂pseudo-code difference observation data; for in moment t frequency f ₁time-differenced phase observation data; for in moment t frequency f ₂time-differenced phase observation data; A, B represent two Satellite Formation Flyings; for Satellite Formation Flying A is in moment t frequency f ₁pseudo-code observation data; for Satellite Formation Flying A is in moment t frequency f ₂pseudo-code observation data; for Satellite Formation Flying A is in moment t frequency f ₁phase observations data; for Satellite Formation Flying A is in moment t frequency f ₂phase observations data; for Satellite Formation Flying B is in moment t frequency f ₁pseudo-code observation data; for Satellite Formation Flying B is in moment t frequency f ₂pseudo-code observation data; for Satellite Formation Flying B is in moment t frequency f ₁phase observations data; for Satellite Formation Flying B is in moment t frequency f ₂phase observations data; all from the first step through pretreated GNSS observation data; for the poor geometric distance of Satellite Formation Flying A, B list between moment t and GNSS satellite j; C represents the light velocity, δ t _aBfor the relative clock correction between Satellite Formation Flying A, B two GNSS receiver; for at the poor ionosphere delay of the list of moment t; for kth (1≤k≤n _b, n _bfor Continuous Tracking segmental arc sum) individual Continuous Tracking segmental arc f ₁the poor phase ambiguity of list of frequency; for a kth Continuous Tracking segmental arc f ₂the poor phase ambiguity of list of frequency; ε is observational error.

Iono-free combination is adopted to eliminate ionosphere delay obtain electric eliminating absciss layer list difference observation equation (2),

$\begin{array}{c}{\mathrm{\ΔP}}_{IF}^j\left(t\right)=\frac{f_1^2}{f_1^2-f_2^2}{\mathrm{\ΔP}}_1^j\left(t\right)-\frac{f_2^2}{f_1^2-f_2^2}{\mathrm{\ΔP}}_2^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)+{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_{IF}}^j\left(t\right)\\ {\mathrm{\ΔL}}_{IF}^j\left(t\right)=\frac{f_1^2}{f_1^2-f_2^2}{\mathrm{\ΔL}}_1^j\left(t\right)-\frac{f_2^2}{f_1^2-f_2^2}{\mathrm{\ΔL}}_2^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)+{\mathrm{\ΔA}}_{IF}^j+{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_{IF}}^j\left(t\right)\end{array}---\left(2\right)$

Wherein " IF " represents iono-free combination; list for moment t poor pseudo-code iono-free combination data; list for moment t poor phase place iono-free combination data; for the kth list that Continuous Tracking segmental arc is a corresponding poor phase place electric eliminating absciss layer blur level.Single poor geometric distance can be calculated by formula (3) below,

${\mathrm{\ρ}}_{AB}^j\left(t\right)=|r_B\left(t\right)-r_G^j\left(t\right)|-|r_A\left(t\right)-r_G^j\left(t\right)|---\left(3\right)$

Wherein || represent the length calculating trivector; r _a(t), r _bt () is for Satellite Formation Flying A, B are in the absolute orbit position of moment t; for GNSS satellite j is in the position of moment t, the interpolation of the precise ephemeris afterwards meter utilizing IGS to provide obtains.

Will at outline point place's linearization launches, and obtains formula (4),

${\mathrm{\ρ}}_{AB}^j\left(t\right)\≈{\mathrm{\ρ}}_{AB,0}^j\left(t\right)-e_B^j\left(t\right)\·{\mathrm{\Δr}}_B\left(t\right)---\left(4\right)$

represent in the direction of visual lines unit vector of moment t from satellite B to GNSS satellite j; for outline point, the absolute orbit position calculation of Satellite Formation Flying A, B that its initial value adopts second step to obtain obtains; Δ r _b(t) for satellite B is in the orbital position improvement of t,

Δr _B(t)＝r _AC(y _B0+Δy _B,t)-r _AC(y _B0,t)(5)

Wherein r _aCrepresent orbital mechanics integral function, r _aC(y _b0, t) represent track mechanical model parameter y _b0carry out satellite B that Adams-Cowell Multi-step Integration the obtains orbital position at moment t.Adams-Cowell Multistep Integrator (see " Astronomica Sinica " the 33rd volume the 4th periodical in 1992 carry " with once and Adams-Cowell method "), integration step gets 10 seconds; y _b0, Δ y _brepresent initial value and the improvement of the orbital mechanics model parameter vector of satellite B respectively, mainly comprise initial time orbital position and speed, solar light pressure coefficient, atmospherical drag coefficient and experience acceleration factor; y _b0adopt the orbital mechanics model parameter that second step exports; Δ y _bleast-squares estimation is adopted to obtain according to formula (6) below.

The electric eliminating absciss layer list in multiple moment difference observation equation (2) is organized into matrix form, obtains single poor observation corrected value vector z _sD,

z _SD＝H _SD·x _SD+e _SD(6)

Wherein " SD " represents single poor; Single poor observation corrected value vector [] ^trepresenting matrix transposition; Design matrix represent partial derivative; Error vector $e_{SD}={\[...,{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_{IF}}^j\left(t\right),{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_{IF}}^j\left(t\right),...\]}^T;$ Solve for parameter vector $x_{SD}={\[...,{\mathrm{\δt}}_{AB}\left(t\right),...,{\mathrm{\ΔA}}_{IF}^{\left(k\right)},...,{\mathrm{\Δy}}_B\]}^T.$ Now solve for parameter mainly comprises: the relative clock correction δ t of receiver _aB(t) (each observation moment 1), single poor phase place electric eliminating absciss layer blur level (each Continuous Tracking segmental arc 1) and orbital mechanics model parameter improve vectorial Δ y _b, wherein vectorial Δ y _bcomprise initial orbital position and speed parameter (3 location coordinates component and 3 speed coordinate components amount to 6 solve for parameters), solar light pressure coefficient (1 solve for parameter), atmospherical drag coefficient (1 solve for parameter) and experience acceleration parameter (every 15 minutes 1 group of solve for parameters).

3.2 adopt weighted least require method to estimate to obtain x _sD,

$x_{SD}={\[H_{SD}^T\·W_{SD}\·H_{SD}\]}^{-1}\·H_{SD}^T\·W_{SD}\·z_{SD}---\left(7\right)$

Wherein single poor weight of observation matrix σ _lrepresent phase observations precision, σ _prepresent pseudo-code accuracy of observation.

3.3 calculate single poor geometric distance.When orbital mechanics model parameter improves vectorial Δ y _bafter being estimated, satellite orbital position r _bt () is by orbital mechanics integration r _aC(y _b0+ Δ y _b, t) obtain, then by r _bt () substitutes into formula (3) and upgrades the single poor geometric distance of calculating

3.4 eliminate single poor geometric distance in difference observation data with clock correction item δ t _aB(t).The poor geometric distance of the list obtained 3.3 with the 3.2 clock correction item δ t obtained _aBt () deducts from the observation data formula (1), obtain the poor observation equation of list (8) after deduction geometric distance of multiple moment and clock correction,

${\mathrm{\ΔP}}_1^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·{\mathrm{\δt}}_{AB}\left(t\right)=I_{AB}^j\left(t\right)+{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_1}^j\left(t\right)$

$\begin{array}{c}{\mathrm{\ΔP}}_2^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·{\mathrm{\δt}}_{AB}\left(t\right)=\frac{f_1^2}{f_2^2}\·I_{AB}^j\left(t\right)+{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_2}^j\left(t\right)\\ {\mathrm{\ΔL}}_1^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·{\mathrm{\δt}}_{AB}\left(t\right)=-I_{AB}^j\left(t\right)+{\mathrm{\λ}}_1\·{\mathrm{\ΔN}}_1^{\left(k\right)}+{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_1}^j\left(t\right)\end{array}---\left(8\right)$

${\mathrm{\ΔL}}_2^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·{\mathrm{\δt}}_{AB}\left(t\right)=-\frac{f_1^2}{f_2^2}\·I_{AB}^j\left(t\right)+{\mathrm{\λ}}_2\·{\mathrm{\ΔN}}_2^{\left(k\right)}+{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_2}^j\left(t\right)$

Due to phase ambiguity in 3.2 single poor relative positionings estimation use only iono-free combination " floating point real number solution form ", be not fixed, therefore obtain formation relative orbit r _b(t)-r _at after the ratio of precision of (), " the fixed integer solution form " of 8.2 is much lower.However, r _b(t)-r _at () precision also reaches several mm magnitude, be far smaller than 1/4 of phase wave length λ, is used as to eliminate geometric distance, auxiliary wide lane ambiguity fixed precision is enough.

4th step, estimate sheet difference phase ambiguity floating-point solution and single poor ionosphere delay parameter.

The 4.1 ratio w that pseudo-code weights and phase place weights are set _p: w _l.Because phase observations precision is far above pseudo-code, type B error code weight w here _pwith phase place weight w _lratio w _p: w _l.Power (the w such as traditional M-W method pseudo-code and phase combination can only be _p: w _l=1:1), the feature of phase data accuracy of observation far above pseudo-code cannot be given full play of, and w of the present invention _p: w _lfreely can arrange according to the accuracy of observation of pseudo-code and phase data, usually get w _p: w _l=σ _l: σ _p, σ _lrepresent phase observations precision, σ _prepresent pseudo-code accuracy of observation, this is the main advantage of the present invention compared with classic method.

4.2 adopt weighted least require method estimate sheet difference phase ambiguity floating-point solution and single poor ionosphere delay parameter.In a kth Continuous Tracking segmental arc, the poor observation equation of list (8) of multiple moment deduction geometric distance and clock correction is organized into matrix form,

z _k＝H _k·x _k+e _k(9)

Wherein single poor observed reading corrects vector $z_k=\left[\begin{array}{c}\·\\ \·\\ \·\\ \mathrm{\Δ}{P}_1^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·\mathrm{\δ}{t}_{AB}\left(t\right)\\ \mathrm{\Δ}{P}_2^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·\mathrm{\δ}{t}_{AB}\left(t\right)\\ \mathrm{\Δ}{L}_1^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·\mathrm{\δ}{t}_{AB}\left(t\right)\\ \mathrm{\Δ}{L}_2^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·\mathrm{\δ}{t}_{AB}\left(t\right)\\ \·\\ \·\\ \·\end{array}\right];$ Error vector $e_k=\left[\begin{array}{c}\·\\ \·\\ \·\\ {\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_1}^j\left(t\right)\\ {\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_2}^j\left(t\right)\\ {\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_1}^j\left(t\right)\\ {\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_2}^j\left(t\right)\\ \·\\ \·\\ \·\end{array}\right];$ Design matrix $H_k=\∂z_k/\∂x_k;$ Solve for parameter vector $x_k={\[...,I_{AB}^j\left(t\right),...,{\mathrm{\ΔN}}_1^{\left(k\right)},{\mathrm{\ΔN}}_2^{\left(k\right)}\]}^T,$ Now solve for parameter mainly comprises: single poor ionosphere delay (each moment 1), 2 poor phase ambiguity parameters of list with for compensating and absorb the impact that satellites formation difference ionosphere delay under Long baselines situation brings.Obtain x _kweighted least square value for,

${\hat{x}}_k={\[H_k^T\·W_k\·H_k\]}^{-1}\·H_k^T\·W_k\·z_k---\left(10\right)$

Wherein weight of observation matrix diag () represents diagonal matrix.Output estimation obtains the poor phase ambiguity parameter of list of a kth Continuous Tracking segmental arc with single poor ionosphere delay parameter

5th step, adopts the conversion of two difference to obtain the wide lane ambiguity floating-point solution of two difference and covariance matrix.

5.1 kth (1≤k≤n that the 4th step is obtained _b) the poor phase ambiguity parameter of list of individual Continuous Tracking segmental arc subtract each other, obtain the corresponding poor wide lane ambiguity of list then n _bdan Chakuan lane phase ambiguity corresponding to individual Continuous Tracking segmental arc is respectively definition floating-point solution vector and covariance matrix the single poor wide lane ambiguity floating-point solution vector of structure because single differential mode type does not introduce correlativity, then single poor wide lane ambiguity is independent of one another, and variance is now $\mathrm{cov}\left({\hat{z}}_{W,SD}\right)=diag(...,{\mathrm{\σ}}_{W,SD}^2,...).$

5.2 according to the relation between two difference and single poor observation equation, structure two difference linear transformation operator T _dD, obtaining the wide lane ambiguity floating-point solution of two difference is,

${\hat{z}}_{W,DD}=T_{DD}\·{\hat{z}}_{W,SD}---\left(11\right)$

Wherein ${\hat{z}}_{W,DD}={\[...,\▿\mathrm{\Δ}{\hat{N}}_W^{\left(jk\right)},\▿\mathrm{\Δ}{\hat{N}}_W^{\left(jm\right)},...\]}^T,{\hat{z}}_{W,SD}={\[...,\mathrm{\Δ}{\hat{N}}_W^{\left(j\right)},\mathrm{\Δ}{\hat{N}}_W^{\left(k\right)},\mathrm{\Δ}{\hat{N}}_W^{\left(m\right)},...\]}^T,$

The wide lane ambiguity covariance matrix of two difference is calculated by formula (12)

$\mathrm{cov}\left({\hat{z}}_{W,DD}\right)=T_{DD}\·\mathrm{cov}\left({\hat{z}}_{W,SD}\right).T_{DD}^T---\left(12\right)$

6th step, fixing Shuan Chakuan lane integer ambiguity.

6.1 adopt integer least square method to fix Shuan Chakuan lane integer ambiguity.The wide lane ambiguity floating-point solution of two difference that 5th step is obtained and covariance matrix as input, adopt integer least square method to fix Shuan Chakuan lane integer ambiguity, namely search for integer space and obtain Shuan Chakuan lane integer ambiguity static solution $z_{W,DD}={\[...,\▿{\mathrm{\ΔN}}_W^{\left(jk\right)},\▿{\mathrm{\ΔN}}_W^{\left(jm\right)},...\]}^T,$ Formula (13) is made to reach minimum,

$\underset{z_{W,DD}\∈Z^n}{\min }\{{\[z_{W,DD}-{\hat{z}}_{W,DD}\]}^T\·{\mathrm{cov}}^{-1}\left({\hat{z}}_{W,DD}\right)\·\[z_{W,DD}-{\hat{z}}_{W,DD}\]\}---\left(13\right)$

Wherein cov ^-1represent the inverse of covariance matrix; Min{} represents and minimizes; Z ⁿrepresent integer space.Least square de-correlation (LAMBDA) is adopted to carry out searching for the integer solution obtaining Shuan Chakuan lane phase ambiguity, " TheLAMBDAmethodforintegerambiguityestimation:implementat ionaspects " (LAMBDA method estimates the realization of integer ambiguity) that LAMBDA method is write in 1996 see PauldeJonge etc.

6.2 inspection Shuan Chakuan lane integer ambiguities fix correctness.Calculate the blur level residual error amount R that the suboptimum integer ambiguity of LAMBDA method output is corresponding _sthe blur level residual error amount R corresponding with optimum integer ambiguity _bratio, if R _s/ R _b≤ k _s/B, Ze Shuanchakuan lane integer ambiguity static solution z _{w, DD}be used, turn the 7th step; If R _s/ R _b>k _s/B, then Shuan Chakuan lane integer ambiguity static solution z is thought _{w, DD}fixing incorrect, the wide lane ambiguity of two difference still adopts its floating-point solution even k _s/Bfor given threshold value (gets k _s/B=2.5) the 7th step, is turned.

7th step, fixing Shuan Chazhai lane integer ambiguity.

The 7.1 Shuan Chakuan lane integer ambiguity static solution z that the 6th step is obtained _{w, DD}deduct from double difference phase observation equation as known quantity, double difference observation comprises two difference pseudo-code and two difference phase place, is poor to the poor observation data of the list of different GNSS satellite again on the basis of single poor observation equation, the relative clock correction of cancellation receiver, namely

$\begin{array}{c}\▿{\mathrm{\ΔP}}_{IF}^{jk}\left(t\right)={\mathrm{\ρ}}_{AB}^{jk}\left(t\right)+{\mathrm{\ϵ}}_{\▿{\mathrm{\ΔP}}_{IF}}^{jk}\left(t\right)\\ \▿{\mathrm{\ΔL}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}\·\▿{\mathrm{\ΔN}}_W^{\left(jk\right)}={\mathrm{\ρ}}_{AB}^{jk}\left(t\right)+\frac{c}{2(f_1+f_2)}\·\▿{\mathrm{\ΔN}}_N^{\left(jk\right)}+{\mathrm{\ϵ}}_{\▿{\mathrm{\ΔL}}_{IF}}^{jk}\left(t\right)\end{array}---\left(14\right)$

Wherein represent two poor; Subscript j, k represent jth, a k GNSS satellite; represent the two difference pseudo-code electric eliminating absciss layer observed readings between GNSS satellite j, k, represent the two difference phase place electric eliminating absciss layer observed readings between GNSS satellite j, k, represent the wide lane ambiguity of two difference and the narrow lane ambiguity of two difference respectively; represent the two difference geometric distances between GNSS satellite j, k, ${\mathrm{\ρ}}_{AB}^{jk}\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)-{\mathrm{\ρ}}_{AB}^k\left(t\right).$

The double difference phase observation equation (14) in multiple moment is organized into matrix form, obtains two difference observation equation of multiple moment (15)

z _DD＝H _DD·x _DD+e _DD(15)

Wherein " DD " represents two poor; z _dDfor double difference observation corrects vector,

$z_{DD}={\[...,\▿{\mathrm{\ΔP}}_{IF}^{jk}\left(t\right)-P_{AB,0}^{jk}\left(t\right),\▿{\mathrm{\ΔL}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}\·\▿N_W^{\left(jk\right)}-{\mathrm{\ρ}}_{AB,0}^{jk}\left(t\right),\·\·\·\]}^T;$

Design matrix $H_{DD}=\∂z_{DD}/\∂x_{DD},$ Error vector $e_{DD}={\[...,{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_{IF}}^{jk}\left(t\right),{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_{IF}}^{jk}\left(t\right),...\]}^T;$ Solve for parameter vector $x_{DD}={\[...,\▿{\mathrm{\ΔN}}_N^{\left(jk\right)},...,{\mathrm{\Δy}}_B\]}^T.$

7.2 according to two difference observation equation of multiple moment (15), adopts the least square estimation method, estimates to obtain narrow lane phase ambiguity floating-point solution vector and covariance matrix

7.3 adopt integer least square method to fix Shuan Chazhai lane integer ambiguity.By narrow lane phase ambiguity floating-point solution and covariance matrix as input, adopt integer least square method to fix Shuan Chazhai lane integer ambiguity, search integer space obtains Shuan Chazhai lane integer ambiguity static solution $z_{N,DD}={\[...,\▿{\mathrm{\ΔN}}_N^{\left(jk\right)},...\]}^T,$ Formula (16) is made to reach minimum,

$\underset{z_{N,DD}\∈Z^n}{\min }\{{\[z_{N,DD}-{\hat{z}}_{N,DD}\]}^T\·{\mathrm{cov}}^{-1}\left({\hat{z}}_{N,DD}\right)\·\[z_{N,DD}-{\hat{z}}_{N,DD}\]\}---\left(16\right)$

7.4 inspection Shuan Chazhai lane integer ambiguities fix correctness.Calculate the blur level residual error amount R that suboptimum integer ambiguity is corresponding _sthe blur level residual error amount R corresponding with optimum integer ambiguity _bratio, if R _s/ R _b≤ k _s/B(get k _s/B=2.5), Ze Shuanchazhai lane integer ambiguity static solution z _{n, DD}be used, turn the 8th step; If R _s/ R _b>k _s/B, then Shuan Chazhai lane integer ambiguity static solution z is thought _{n, DD}fixing incorrect, the narrow lane ambiguity of two difference still adopts its floating-point solution even turn the 8th step.

8th step, exports the relative positioning result that blur level is fixing.

The 8.1 Shuan Chazhai lane integer ambiguity static solutions that the 7th step is obtained deduct from double difference phase observation equation (14) again as known quantity, obtain two differences observation equation (17) after deducting blur level.

$\begin{array}{c}\▿{\mathrm{\ΔP}}_{IF}^{jk}\left(t\right)={\mathrm{\ρ}}_{AB}^{jk}\left(t\right)+{\mathrm{\ϵ}}_{\▿{\mathrm{\ΔP}}_{IF}}^{jk}\left(t\right)\\ \▿{\mathrm{\ΔL}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}\·\▿{\mathrm{\ΔN}}_W^{\left(jk\right)}-\frac{c}{2(f_1+f_2)}\·\▿{\mathrm{\ΔN}}_N^{\left(jk\right)}={\mathrm{\ρ}}_{AB}^{jk}\left(t\right)+{\mathrm{\ϵ}}_{\▿{\mathrm{\ΔL}}_{IF}}^{jk}\left(t\right)\end{array}---\left(17\right)$

Two differences observation equation (17) after the deduction blur level in multiple moment are organized into matrix form, obtain two differences observation equation (18) after deduction blur level of multiple moment

$z_{DD}^{*}=H_{DD}^{*}\·{\mathrm{\Δy}}_B+e_{DD}---\left(18\right)$

Wherein double difference observation corrects vector,

$z_{DD}^{*}={\[...,\▿{\mathrm{\ΔP}}_{IF}^{jk}\left(t\right)-{\mathrm{\ρ}}_{AB,0}^{jk}\left(t\right),\▿{\mathrm{\ΔL}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}\·\▿{\mathrm{\ΔN}}_W^{\left(jk\right)}-\frac{c}{2(f_1-f_2)}\·\▿{\mathrm{\ΔN}}_N^{\left(jk\right)}-{\mathrm{\ρ}}_{AB,0}^{jk}\left(t\right),...\]}^T;$

Design matrix solve for parameter vector Δ y _bfor orbital mechanics model parameter improves vector.

8.2, according to two differences observation equation (18) after multiple moment deduction blur level, adopt the least square estimation method, estimate that obtaining orbital mechanics model parameter improves vectorial Δ y _b, satellite orbital position r _bt () is by orbital mechanics integration r _aC(y _b0+ Δ y _b, t) obtain, thus export final blur level fix after formation relative orbit r _b(t)-r _a(t).Compared with 3.1 single poor relative positionings, in 8.2, phase ambiguity is fixed with the form in Kuan Xianghezhai lane, therefore 8.2 obtain final blur level fix after formation relative orbit precision than 3.1 " floating-point solution form " much higher.

The present invention compared with prior art tool has the following advantages:

The present invention utilizes relative positioning floating-point solution to eliminate difference geometric distance and clock correction, pass through pseudo-code, the combination of phase place differential weights comes estimate sheet difference phase ambiguity and ionosphere delay parameter, and obtain the wide lane ambiguity floating-point solution of two difference and covariance matrix thereof through the conversion of two difference, and then utilize integer least square method to realize wide lane integer ambiguity to fix, avoid pseudo-code in classic method M-W combination, the power such as phase place processes the problem that the blur level caused fixedly is strongly depend on the lower pseudo-code of accuracy of observation, improve the precision that Long baselines satellites formation GNSS relative positioning blur level is fixed into power and final relative positioning result.

The present invention adopts batch processing least square processing mode, calculation stability, even if the solid error of local segmental arc blur level also can't cause the algorithm of whole time period not restrained, overcome the shortcoming that Sequential filter disposal route is easily dispersed, improve the reliability of relative positioning result.

Method of the present invention has versatility, is applicable to the high-precision GNSS relative positioning application of Long baselines satellites formation and other moving target.

Accompanying drawing explanation

Fig. 1 is Long baselines satellites formation GNSS relative positioning blur level fixed solution process flow diagram of the present invention;

Fig. 2 is the single poor ionosphere delay estimated result figure of GRACE of the present invention formation.

Fig. 3 is the present invention and classic method GRACE formation positioning precision comparison chart;

Embodiment

Below in conjunction with accompanying drawing, the present invention is described further.For GRACE satellites formation GNSS relative positioning, interstellar distance 140 km in this formation in January, 2006, belongs to typical Long baselines satellites formation, chooses totally 7 days measured datas in-orbit 1 to 7 February in 2006.As shown in Figure 1, the Long baselines satellites formation GNSS relative positioning method that a kind of blur level of the present invention is fixing comprises the steps:

The first step, collects and pre-service input data.

Collect observation data from satellite platform and collect auxiliary data (see table 1) from website.Pre-service is carried out to the GNSS observation data of 7 days GRACE satellites formations, mark outlier and cycle slip, and data pretreated every day are preserved with text form, export pretreated GNSS observation data file (* .edt).

Table 1 observation data and auxiliary positioning input data

Sequence number	Data type	Remarks
			1	GNSS observation data	Rinex (* .yyo), GNSS receiver, star passes up and down
2	Satellite platform attitude data	Hypercomplex number (* .att), star passes up and down
			3	GNSS satellite ephemeris	Rinex (* .sp3), website is downloaded
4	GNSS satellite clock correction	Rinex (* .clk), website is downloaded
			5	GNSS antenna PCV information	Igs08.atx, website is downloaded
6	Earth gravity field	GGM02C, website is downloaded
			7	Earth rotation parameter	IERS 2000A, website is downloaded
8	UTC time jump second	TAI-UTC, website is downloaded
			9	JPL solar calendar	DE405, website is downloaded
10	Solar radiation flow	Every day 1 record, website download
			11	Geomagnetic indices	Every 3 hours 1 records, website is downloaded

Second step, determines every absolute outline track of Satellite Formation Flying.

Adopt non-poor reduced-dynamic method, pretreated GNSS observation data file is processed, export outline orbital mechanics model parameter and the absolute orbit position of GRACEA and GRACEB satellite.GRACE for interstellar distance 140 km forms into columns, and in order to make relative positioning reach 1mm precision, now requires that absolute orbit Product Precision needs to be better than 10cm.Table 2 gives absolute orbit product and the science track comparison result in GRACE satellites formation 1 to 7 February in 2006, and visible 3d orbit precision on average reaches about 3.5cm, meets follow-up high precision relative positioning needs.

Table 2GRACE formation absolute orbit product and science track comparison result

Date	GRACE A/cm	GRACE B/cm
			2006-02-01	3.94	4.26
2006-02-02	3.67	4.01
			2006-02-03	3.41	3.15
2006-02-04	3.48	3.64
			2006-02-05	4.00	4.03
2006-02-06	3.49	3.07
			2006-02-07	2.71	3.06
Mean value	3.52	3.60

3rd step, eliminates the geometric distance in difference observation data and clock correction.

Carry out the poor relative positioning of list, improve orbital mechanics model parameter and the orbital position of the GRACEB satellite that second step obtains, obtain single poor relative positioning floating-point solution.Calculate single poor geometric distance, eliminate single poor geometric distance and clock correction in difference observation data.Table 4 gives the list poor relative positioning floating-point solution precision in GRACE satellites formation 1 to 7 February in 2006, its mean accuracy visible reaches about 4mm, be far smaller than 1/4 of phase wave length λ, be used as to eliminate geometric distance and assist wide lane ambiguity fixed precision enough.

The KBR of the single poor relative positioning floating-point solution of table 4GRACE formation checks precision

Date	KBR comparison standard deviation/mm
		2006-02-01	2.90
2006-02-02	3.70
		2006-02-03	4.38
2006-02-04	4.59
		2006-02-05	3.82
2006-02-06	4.23
		2006-02-07	4.17
Mean value	3.97

4th step, by pseudo-code, the combination of phase place differential weights value, adopts least square method estimate sheet difference phase ambiguity floating-point solution and single poor ionosphere delay parameter.

The pseudo-code weights of GRACE satellite and the ratio w of phase place weights are set _p: w _l=1:100, to the 3rd step obtain deducted geometric distance and clock correction after the poor pseudo-code of list and phase observations equation adopt least square method to solve, export single poor phase ambiguity floating-point solution and single poor ionosphere delay parameter estimation result.Fig. 2 gives the GRACE estimated result of forming into columns single poor ionosphere delay, visible difference ionosphere delay reaches about 3 meters, obviously be greater than 1/4 of phase wave length λ, therefore for satellites formation under Long baselines situation, ionosphere delay impact directly cannot be eliminated by difference, it can not be left in the basket to blur level fixed effect, must absorb by increasing parametric compensation.

5th step, adopts the conversion of two difference to obtain the wide lane ambiguity floating-point solution of two difference and covariance matrix.

The different frequency list difference phase ambiguity floating-point solution 4th step obtained is subtracted each other, and obtains single poor wide lane ambiguity floating-point solution, and according to two difference and single poor between relation, through the linear transformation of two difference, calculate the wide lane ambiguity floating-point solution of two difference and covariance matrix.

6th step, fixing Shuan Chakuan lane integer ambiguity.

The wide lane ambiguity floating-point solution of two difference 5th step obtained and covariance matrix thereof, as input, adopt integer least square method to fix Shuan Chakuan lane integer ambiguity.Table 5 gives GRACE satellites formation 1 to 7 February in 2006, and the wide lane ambiguity success ratio comparing result of the present invention and traditional M-W method, visible integer ambiguity of the present invention is fixed into power and improves 5% relative to traditional M-W method.

The GRACE of table 5 the present invention and classic method wide lane ambiguity success ratio of forming into columns contrasts

7th step, fixing Shuan Chazhai lane integer ambiguity.

The Shuan Chakuan lane integer ambiguity static solution 6th step obtained is being deducted from double difference phase observation equation as known quantity, adopt the least square estimation method to estimate to obtain narrow lane phase ambiguity floating-point solution vector and covariance matrix thereof, adopt integer least square method to fix Shuan Chazhai lane integer ambiguity.

8th step, exports the relative positioning result that blur level is fixing.

The Shuan Chazhai lane integer ambiguity static solution 7th step obtained is deducted from double difference phase observation equation as known quantity, adopt the least square estimation method, estimate obtain orbital mechanics model parameter improve vector, after orbit integration, export final blur level fix after formation relative orbit.Fig. 3 gives the GRACE formation relative positioning KBR comparison result of the present invention and classic method, visible the present invention is while raising integer ambiguity is fixed into power, also effectively improve the precision of relative positioning result, KBR comparison standard deviation is significantly reduced, average result is reduced to 0.78mm by the 0.91mm of classic method, positioning precision is significantly improved, and many days result of calculation is more stable.

Claims (4)

1. the Long baselines satellites formation GNSS relative positioning method that blur level is fixing, is characterized in that comprising the following steps:

The first step, collect and pre-service input data:

1.1 collect observation data from satellite platform: GNSS observation data, satellite platform attitude data; Auxiliary data is collected: GNSS satellite precise ephemeris, clock correction and antenna PCV information, earth gravity field, earth rotation parameter, UTC time jump second data, JPL solar calendar, solar radiation flow, geomagnetic indices from website;

1.2 pairs of observation datas of collecting carry out pre-service: reject the outlier of GNSS observation data, and detecting phase cycle slip;

Second step, determines the absolute orbit position of Satellite Formation Flying A, B, exports the orbital mechanics model parameter y of Satellite Formation Flying A, B _a0, y _b0with the absolute orbit position r of each moment t _a(t), r _b(t), orbital position comprises x, y, z three cartesian component;

3rd step, eliminate the geometric distance in difference observation data and clock correction:

3.1 single poor relative positionings, method is: single poor be same GNSS satellite to two receivers observations, same kind observation data does difference, form difference observation data, adopt the poor observation equation of list (1) of t to eliminate GNSS satellite clock correction, retain the relative clock correction of receiver

$\begin{array}{c}{\mathrm{\ΔP}}_1^j\left(t\right)=P_{1,A}^j\left(t\right)-P_{1,B}^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)+I_{AB}^j\left(t\right)+{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_1}^j\left(t\right)\\ {\mathrm{\ΔP}}_2^j\left(t\right)=P_{2,A}^j\left(t\right)-P_{2,B}^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)+\frac{f_1^2}{f_2^2}\·I_{AB}^j\left(t\right)+{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_2}^j\left(t\right)\\ {\mathrm{\ΔL}}_1^j\left(t\right)=L_{1,A}^j\left(t\right)-L_{1,B}^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)-I_{AB}^j\left(t\right)+{\mathrm{\λ}}_1\·{\mathrm{\ΔN}}_1^{\left(k\right)}+{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_1}^j\left(t\right)\\ {\mathrm{\ΔL}}_2^j\left(t\right)=L_{2,A}^j\left(t\right)-L_{2,B}^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)-\frac{f_1^2}{f_2^2}\·I_{AB}^j\left(t\right)+{\mathrm{\λ}}_2\·{\mathrm{\ΔN}}_2^{\left(k\right)}+{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_2}^j\left(t\right)\end{array}---\left(1\right)$

Wherein t represents the observation moment, and subscript j represents a jth GNSS satellite, and 1≤j≤M, M is the sum of satellite; Subscript 1,2 represents GNSS signal frequency f ₁, f ₂; Δ represents single poor; for in moment t frequency f ₁pseudo-code difference observation data; for in moment t frequency f ₂pseudo-code difference observation data; for in moment t frequency f ₁time-differenced phase observation data; for in moment t frequency f ₂phase data difference observation data; A, B represent two Satellite Formation Flyings; for Satellite Formation Flying A is in moment t frequency f ₁pseudo-code observation data; for Satellite Formation Flying A is in moment t frequency f ₂pseudo-code observation data; for Satellite Formation Flying A is in moment t frequency f ₁phase observations data; for Satellite Formation Flying A is in moment t frequency f ₂phase observations data; for Satellite Formation Flying B is in moment t frequency f ₁pseudo-code observation data; for Satellite Formation Flying B is in moment t frequency f ₂pseudo-code observation data; for Satellite Formation Flying B is in moment t frequency f ₁phase observations data; for Satellite Formation Flying B is in moment t frequency f ₂phase observations data; all from the first step through pretreated GNSS observation data; for the poor geometric distance of Satellite Formation Flying A, B list between moment t and GNSS satellite j; C represents the light velocity, δ t _aBfor the relative clock correction between Satellite Formation Flying A, B two GNSS receiver; for at the poor ionosphere delay of the list of moment t; for a kth Continuous Tracking segmental arc f ₁the poor phase ambiguity of list of frequency, 1≤k≤n _b, n _bfor Continuous Tracking segmental arc sum; for a kth Continuous Tracking segmental arc f ₂the poor phase ambiguity of list of frequency; ε is observational error;

Iono-free combination is adopted to eliminate ionosphere delay obtain electric eliminating absciss layer list difference observation equation (2),

$\begin{array}{c}{\mathrm{\ΔP}}_{IF}^j\left(t\right)=\frac{f_1^2}{f_1^2-f_2^2}{\mathrm{\ΔP}}_1^j\left(t\right)-\frac{f_2^2}{f_1^2-f_2^2}{\mathrm{\ΔP}}_2^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)+{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_{IF}}^j\left(t\right)\\ {\mathrm{\ΔL}}_{IF}^j\left(t\right)=\frac{f_1^2}{f_1^2-f_2^2}{\mathrm{\ΔL}}_1^j\left(t\right)-\frac{f_2^2}{f_1^2-f_2^2}{\mathrm{\ΔL}}_2^j\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)+c\·{\mathrm{\δt}}_{AB}\left(t\right)+{\mathrm{\ΔA}}_{IF}^j+{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_{IF}}^j\left(t\right)\end{array}---\left(2\right)$

Wherein " IF " represents iono-free combination; list for moment t poor pseudo-code iono-free combination data; list for moment t poor phase place iono-free combination data; for the kth list that Continuous Tracking segmental arc is a corresponding poor phase place electric eliminating absciss layer blur level; Single poor geometric distance calculated by formula (3),

${\mathrm{\ρ}}_{AB}^j\left(t\right)=|r_B\left(t\right)-r_G^j\left(t\right)|-|r_A\left(t\right)-r_G^j\left(t\right)|---\left(3\right)$

Wherein || represent the length calculating trivector; r _a(t), r _bt () is for Satellite Formation Flying A, B are in the absolute orbit position of moment t; for GNSS satellite j is in the position of moment t, the interpolation of the precise ephemeris afterwards meter utilizing IGS to provide obtains;

Will at outline point place's linearization launches, and obtains formula (4),

${\mathrm{\ρ}}_{AB}^j\left(t\right)\≈{\mathrm{\ρ}}_{AB,0}^j\left(t\right)-e_B^j\left(t\right)\·{\mathrm{\Δr}}_B\left(t\right)---\left(4\right)$

represent in the direction of visual lines unit vector of moment t from satellite B to GNSS satellite j; for outline point, the absolute orbit position calculation of Satellite Formation Flying A, B that its initial value adopts second step to obtain obtains; Δ r _b(t) for satellite B is in the orbital position improvement of t,

Δr _B(t)＝r _AC(y _B0+Δy _B,t)-r _AC(y _B0,t)(5)

Wherein r _aCrepresent orbital mechanics integral function, r _aC(y _b0, t) represent track mechanical model parameter y _b0carry out satellite B that Adams-Cowell Multi-step Integration the obtains orbital position at moment t; y _b0, Δ y _brepresent initial value and the improvement of the orbital mechanics model parameter vector of satellite B respectively, comprise initial time orbital position and speed, solar light pressure coefficient, atmospherical drag coefficient and experience acceleration factor; y _b0adopt the orbital mechanics model parameter that second step exports; Δ y _bleast-squares estimation is adopted to obtain according to formula (6);

The electric eliminating absciss layer list in multiple moment difference observation equation (2) is organized into matrix form, obtains single poor observation corrected value vector z _sD,

z _SD＝H _SD·x _SD+e _SD(6)

Wherein " SD " represents single poor; Single poor observation corrected value vector $z_{SD}={\[...,{\mathrm{\ΔP}}_{IF}^j\left(t\right)-{\mathrm{\ρ}}_{AB,0}^j\left(t\right),{\mathrm{\ΔL}}_{IF}^j\left(t\right)-{\mathrm{\ρ}}_{AB,0}^j\left(t\right),...\]}^T,$ [] ^trepresenting matrix transposition; Design matrix represent partial derivative; Error vector $e_{SD}={\[...,{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_{IF}}^j\left(t\right),{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_{IF}}^j\left(t\right),...\]}^T;$ Solve for parameter vector $x_{SD}={\[...,{\mathrm{\δt}}_{AB}\left(t\right),...,{\mathrm{\ΔA}}_{IF}^{\left(k\right)},...,{\mathrm{\Δy}}_B\]}^T;$ Now solve for parameter comprises: the relative clock correction δ t of receiver _aB(t), single poor phase place electric eliminating absciss layer blur level vectorial Δ y is improved with orbital mechanics model parameter _b;

3.2 adopt weighted least require method to estimate to obtain x _sD,

$x_{SD}={\[H_{SD}^T\·W_{SD}\·H_{SD}\]}^{-1}\·H_{SD}^T\·W_{SD}\·z_{SD}---\left(7\right)$

Wherein single poor weight of observation matrix $W_{SD}=diag(...,1/2{\mathrm{\σ}}_P^2,1/2{\mathrm{\σ}}_L^2,...),$ σ _lrepresent phase observations precision, σ _prepresent pseudo-code accuracy of observation;

3.3 calculate single poor geometric distance: satellite orbital position r _bt () is by Adams-Cowell Multi-step Integration r _aC(y _b0+ Δ y _b, t) obtain, then by r _bt () substitutes into formula (3) and upgrades the single poor geometric distance of calculating

3.4 eliminate single poor geometric distance in difference observation data with clock correction item δ t _aB(t): the poor geometric distance of the list obtained 3.3 with the 3.2 clock correction item δ t obtained _aBt () deducts from the observation data formula (1), obtain the poor observation equation of list (8) after deduction geometric distance of multiple moment and clock correction,

$\begin{array}{c}{\mathrm{\ΔP}}_1^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·{\mathrm{\δt}}_{AB}\left(t\right)=I_{AB}^j\left(t\right)+{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_1}^j\left(t\right)\\ {\mathrm{\ΔP}}_2^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·{\mathrm{\δt}}_{AB}\left(t\right)=\frac{f_1^2}{f_2^2}\·I_{AB}^j\left(t\right)+{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_2}^j\left(t\right)\\ {\mathrm{\ΔL}}_1^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·{\mathrm{\δt}}_{AB}\left(t\right)=-I_{AB}^j\left(t\right)+{\mathrm{\λ}}_1\·{\mathrm{\ΔN}}_1^{\left(k\right)}+{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_1}^j\left(t\right)\\ {\mathrm{\ΔL}}_2^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·{\mathrm{\δt}}_{AB}\left(t\right)=-\frac{f_1^2}{f_2^2}\·I_{AB}^j\left(t\right)+{\mathrm{\λ}}_2\·{\mathrm{\ΔN}}_2^{\left(k\right)}+{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_2}^j\left(t\right)\end{array}---\left(8\right)$

4th step, estimate sheet difference phase ambiguity floating-point solution and single poor ionosphere delay parameter:

4.1 arrange pseudo-code weight w _pwith phase place weight w _lratio, get w _p: w _l=σ _l: σ _p, σ _lrepresent phase observations precision, σ _prepresent pseudo-code accuracy of observation;

4.2 adopt weighted least require method estimate sheet difference phase ambiguity floating-point solution and single poor ionosphere delay parameter; In a kth Continuous Tracking segmental arc, the poor observation equation of list (8) of multiple moment deduction geometric distance and clock correction is organized into matrix form,

z _k＝H _k·x _k+e _k(9)

Wherein single poor observed reading corrects vector $z_k=\left[\begin{array}{c}\·\\ \·\\ \·\\ \mathrm{\Δ}{P}_1^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·\mathrm{\δ}{t}_{AB}\left(t\right)\\ \mathrm{\Δ}{P}_2^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·\mathrm{\δ}{t}_{AB}\left(t\right)\\ \mathrm{\Δ}{L}_1^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·\mathrm{\δ}{t}_{AB}\left(t\right)\\ \mathrm{\Δ}{L}_2^j\left(t\right)-{\mathrm{\ρ}}_{AB}^j\left(t\right)-c\·\mathrm{\δ}{t}_{AB}\left(t\right)\\ \·\\ \·\\ \·\end{array}\right];$ Error vector $e_k=\left[\begin{array}{c}\·\\ \·\\ \·\\ {\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_1}^j\left(t\right)\\ {\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_2}^j\left(t\right)\\ {\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_1}^j\left(t\right)\\ {\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_2}^j\left(t\right)\\ \·\\ \·\\ \·\end{array}\right];$ Design matrix solve for parameter vector now solve for parameter comprises: single poor ionosphere delay, 2 poor phase ambiguity parameters of list with obtain x _kweighted least square value for,

${\hat{x}}_k={\[H_k^T\·W_k\·H_k\]}^{-1}\·H_k^T\·W_k\·z_k---\left(10\right)$

Wherein weight of observation matrix $W_k=diag(...,w_P^2/w_L^2,w_P^2/w_L^2,1,1,...),$ Diag () represents diagonal matrix; Output estimation obtains the poor phase ambiguity parameter of list of a kth Continuous Tracking segmental arc with single poor ionosphere delay parameter

5th step, adopts the conversion of two difference to obtain the wide lane ambiguity floating-point solution of two difference and covariance matrix:

The poor phase ambiguity parameter of list of the 5.1 kth Continuous Tracking segmental arc that the 4th step is obtained subtract each other, obtain the corresponding poor wide lane ambiguity of list then n _bdan Chakuan lane phase ambiguity corresponding to individual Continuous Tracking segmental arc is respectively definition floating-point solution vector and covariance matrix the single poor wide lane ambiguity floating-point solution vector of structure variance is now $\mathrm{cov}\left({\hat{z}}_{W,SD}\right)=diag(...,{\mathrm{\σ}}_{W,SD}^2,...);$

5.2 according to the relation between two difference and single poor observation equation, structure two difference linear transformation operator T _dD, obtaining the wide lane ambiguity floating-point solution of two difference is,

${\hat{z}}_{W,DD}=T_{DD}\·{\hat{z}}_{W,SD}---\left(11\right)$

Wherein ${\hat{z}}_{W,DD}={\[...,\▿\mathrm{\Δ}{\hat{N}}_W^{\left(jk\right)},\▿\mathrm{\Δ}{\hat{N}}_W^{\left(jm\right)},...\]}^T,$ ${\hat{z}}_{W,SD}={\[...,\mathrm{\Δ}{\hat{N}}_W^{\left(j\right)},\mathrm{\Δ}{\hat{N}}_W^{\left(k\right)},\mathrm{\Δ}{\hat{N}}_W^{\left(m\right)},...\]}^T,$

The wide lane ambiguity covariance matrix of two difference is calculated by formula (12)

$\mathrm{cov}\left({\hat{z}}_{W,DD}\right)=T_{DD}\·\mathrm{cov}\left({\hat{z}}_{W,SD}\right)\·T_{DD}^T---\left(12\right)$

6th step, fixing Shuan Chakuan lane integer ambiguity:

6.1 adopt integer least square method to fix Shuan Chakuan lane integer ambiguity: the wide lane ambiguity floating-point solution of two difference the 5th step obtained and covariance matrix as input, adopt integer least square method to fix Shuan Chakuan lane integer ambiguity, namely search for integer space and obtain Shuan Chakuan lane integer ambiguity static solution formula (13) is made to reach minimum,

$\underset{z_{W,DD}\∈Z^n}{\min }\{{\[z_{W,DD}-{\hat{z}}_{W,DD}\]}^T\·{\mathrm{cov}}^{-1}\left({\hat{z}}_{W,DD}\right)\·\[z_{W,DD}-{\hat{z}}_{W,DD}\]\}---\left(13\right)$

Wherein cov ^-1represent the inverse of covariance matrix; Min{} represents and minimizes; Z ⁿrepresent integer space; Least square de-correlation LAMBDA is adopted to carry out searching for the integer solution obtaining Shuan Chakuan lane phase ambiguity;

6.2 inspection Shuan Chakuan lane integer ambiguities fix correctness: calculate the blur level residual error amount R that the suboptimum integer ambiguity of LAMBDA method output is corresponding _sthe blur level residual error amount R corresponding with optimum integer ambiguity _bratio, if R _s/ R _b≤ k _s/B, Ze Shuanchakuan lane integer ambiguity static solution z _{w, DD}be used, turn the 7th step; If R _s/ R _b>k _s/B, Ze Shuanchakuan lane integer ambiguity static solution z _{w, DD}fixing incorrect, the wide lane ambiguity of two difference still adopts its floating-point solution even k _s/Bfor given threshold value, turn the 7th step;

7th step, fixing Shuan Chazhai lane integer ambiguity:

The 7.1 Shuan Chakuan lane integer ambiguity static solution z that the 6th step is obtained _{w, DD}deduct from double difference phase observation equation as known quantity, double difference observation comprises two difference pseudo-code and two difference phase place, is poor to the poor observation data of the list of different GNSS satellite again on the basis of single poor observation equation, the relative clock correction of cancellation receiver, namely

$\▿{\mathrm{\ΔP}}_{IF}^{jk}\left(t\right)={\mathrm{\ρ}}_{AB}^{jk}\left(t\right)+{\mathrm{\ϵ}}_{\▿{\mathrm{\ΔP}}_{IF}}^{jk}\left(t\right)$

$\▿{\mathrm{\ΔL}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}\·\▿{\mathrm{\ΔN}}_W^{\left(jk\right)}={\mathrm{\ρ}}_{AB}^{jk}\left(t\right)+\frac{c}{2(f_1+f_2)}\·\▿{\mathrm{\ΔN}}_N^{\left(jk\right)}+{\mathrm{\ϵ}}_{\▿{\mathrm{\ΔL}}_{IF}}^{jk}\left(t\right)---\left(14\right)$

Wherein represent two poor; Subscript j, k represent jth, a k GNSS satellite; represent the two difference pseudo-code electric eliminating absciss layer observed readings between GNSS satellite j, k, represent the two difference phase place electric eliminating absciss layer observed readings between GNSS satellite j, k, $\▿{\mathrm{\ΔL}}_{IF}^{jk}\left(t\right)={\mathrm{\ΔL}}_{IF}^j\left(t\right)-{\mathrm{\ΔL}}_{IF}^k\left(t\right);$ represent the wide lane ambiguity of two difference and the narrow lane ambiguity of two difference respectively; represent the two difference geometric distances between GNSS satellite j, k, ${\mathrm{\ρ}}_{AB}^{jk}\left(t\right)={\mathrm{\ρ}}_{AB}^j\left(t\right)-{\mathrm{\ρ}}_{AB}^k\left(t\right);$

The double difference phase observation equation (14) in multiple moment is organized into matrix form, obtains two difference observation equation of multiple moment (15)

z _DD＝H _DD·x _DD+e _DD(15)

Wherein " DD " represents two poor; z _dDfor double difference observation corrects vector,

$z_{DD}={\[...,\▿{\mathrm{\ΔP}}_{IF}^{jk}\left(t\right)-P_{AB,0}^{jk}\left(t\right),\▿{\mathrm{\ΔL}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}\·\▿N_W^{\left(jk\right)}-{\mathrm{\ρ}}_{AB,0}^{jk}\left(t\right),\·\·\·\]}^T;$

Design matrix error vector $e_{DD}={\[...,{\mathrm{\ϵ}}_{{\mathrm{\ΔP}}_{IF}}^{jk}\left(t\right),{\mathrm{\ϵ}}_{{\mathrm{\ΔL}}_{IF}}^{jk}\left(t\right),...\]}^T;$ Solve for parameter vector $x_{DD}={\[...,\▿{\mathrm{\ΔN}}_N^{\left(jk\right)},...,{\mathrm{\Δy}}_B\]}^T;$

7.2 according to two difference observation equation of multiple moment (15), adopts the least square estimation method, estimates to obtain narrow lane phase ambiguity floating-point solution vector and covariance matrix

7.3 adopt integer least square method to fix Shuan Chazhai lane integer ambiguity; By narrow lane phase ambiguity floating-point solution and covariance matrix as input, adopt integer least square method to fix Shuan Chazhai lane integer ambiguity, search integer space obtains Shuan Chazhai lane integer ambiguity static solution formula (16) is made to reach minimum,

$\underset{z_{N,DD}\∈Z^n}{\min }\{{\[z_{N,DD}-{\hat{z}}_{N,DD}\]}^T\·{\mathrm{cov}}^{-1}\left({\hat{z}}_{N,DD}\right)\·\[z_{N,DD}-{\hat{z}}_{N,DD}\]\}---\left(16\right)$

7.4 inspection Shuan Chazhai lane integer ambiguities fix correctness: calculate the blur level residual error amount R that suboptimum integer ambiguity is corresponding _sthe blur level residual error amount R corresponding with optimum integer ambiguity _bratio, if R _s/ R _b≤ k _s/B, Ze Shuanchazhai lane integer ambiguity static solution z _{n, DD}be used, turn the 8th step; If R _s/ R _b>k _s/B, Ze Shuanchazhai lane integer ambiguity static solution z _{n, DD}fixing incorrect, the narrow lane ambiguity of two difference still adopts its floating-point solution even turn the 8th step;

8th step, exports the relative positioning result that blur level is fixing:

The 8.1 Shuan Chazhai lane integer ambiguity static solutions that the 7th step is obtained deduct from double difference phase observation equation (14) again as known quantity, obtain two differences observation equation (17) after deducting blur level:

$\▿{\mathrm{\ΔP}}_{IF}^{jk}\left(t\right)={\mathrm{\ρ}}_{AB}^{jk}\left(t\right)+{\mathrm{\ϵ}}_{\▿{\mathrm{\ΔP}}_{IF}}^{jk}\left(t\right)$

$\▿{\mathrm{\ΔL}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}\·\▿{\mathrm{\ΔN}}_W^{\left(jk\right)}-\frac{c}{2(f_1+f_2)}\·\▿{\mathrm{\ΔN}}_N^{\left(jk\right)}={\mathrm{\ρ}}_{AB}^{jk}\left(t\right)+{\mathrm{\ϵ}}_{\▿{\mathrm{\ΔL}}_{IF}}^{jk}\left(t\right)---\left(17\right)$

Two differences observation equation (17) after the deduction blur level in multiple moment are organized into matrix form, obtain two differences observation equation (18) after deduction blur level of multiple moment

$z_{DD}^{*}=H_{DD}^{*}\·{\mathrm{\Δy}}_B+e_{DD}---\left(18\right)$

Wherein double difference observation corrects vector,

$z_{DD}^{*}={\[...,\▿{\mathrm{\ΔP}}_{IF}^{jk}\left(t\right)-{\mathrm{\ρ}}_{AB,0}^{jk}\left(t\right),\▿{\mathrm{\ΔL}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}\·\▿{\mathrm{\ΔN}}_W^{\left(jk\right)}-\frac{c}{2(f_1-f_2)}\·\▿{\mathrm{\ΔN}}_N^{\left(jk\right)}-{\mathrm{\ρ}}_{AB,0}^{jk}\left(t\right),...\]}^T;$

Design matrix solve for parameter vector Δ y _bfor orbital mechanics model parameter improves vector;

8.2, according to two differences observation equation (18) after multiple moment deduction blur level, adopt the least square estimation method, estimate that obtaining orbital mechanics model parameter improves vectorial Δ y _b, satellite orbital position r _bt () is by orbital mechanics integration r _aC(y _b0+ Δ y _b, t) obtain, export final blur level fix after formation relative orbit r _b(t)-r _a(t).

2. the Long baselines satellites formation GNSS relative positioning method that a kind of blur level as claimed in claim 1 is fixing, it is characterized in that the method for the absolute orbit position of second step determination Satellite Formation Flying A, B is non-poor reduced-dynamic method, non-poor reduced-dynamic method is shown in the monograph measurement data modeling that ZhengmingWang etc. published in 2011 and parameter estimation Section 7.3.

3. the Long baselines satellites formation GNSS relative positioning method that a kind of blur level as claimed in claim 1 is fixing, is characterized in that the 3rd step Adams-Cowell Multistep Integrator step-length gets 10 seconds.

4. the Long baselines satellites formation GNSS relative positioning method that a kind of blur level as claimed in claim 1 is fixing, is characterized in that described k _s/B=2.5.