CN112147663A - Satellite and inertia combined dynamic-alignment real-time precise relative positioning method - Google Patents

Satellite and inertia combined dynamic-alignment real-time precise relative positioning method Download PDF

Info

Publication number
CN112147663A
CN112147663A CN202011327522.4A CN202011327522A CN112147663A CN 112147663 A CN112147663 A CN 112147663A CN 202011327522 A CN202011327522 A CN 202011327522A CN 112147663 A CN112147663 A CN 112147663A
Authority
CN
China
Prior art keywords
satellite
epoch
relative positioning
dynamic
time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011327522.4A
Other languages
Chinese (zh)
Other versions
CN112147663B (en
Inventor
李青松
吴杰
董毅
王鼎杰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202011327522.4A priority Critical patent/CN112147663B/en
Publication of CN112147663A publication Critical patent/CN112147663A/en
Application granted granted Critical
Publication of CN112147663B publication Critical patent/CN112147663B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/45—Determining position by combining measurements of signals from the satellite radio beacon positioning system with a supplementary measurement
    • G01S19/47—Determining position by combining measurements of signals from the satellite radio beacon positioning system with a supplementary measurement the supplementary measurement being an inertial measurement, e.g. tightly coupled inertial
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/165—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/43—Determining position using carrier phase measurements, e.g. kinematic positioning; using long or short baseline interferometry
    • G01S19/44—Carrier phase ambiguity resolution; Floating ambiguity; LAMBDA [Least-squares AMBiguity Decorrelation Adjustment] method

Landscapes

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

Abstract

The invention provides a satellite and inertia combined dynamic-to-dynamic real-time precise relative positioning method, which obtains a relative positioning vector between a mobile station and a dynamic reference based on a carrier phase real-time dynamic differential relative positioning method, and obtains a measurement updating time interval delta based on a carrier phase GNSS/inertial navigation compact combined algorithmt TC Inner mobile station position increment, moving reference position increment and mobile station forecast time deltat p And (4) obtaining the forecasting time of the moving reference by a polynomial forecasting method based on a sliding window, and finally synthesizing and outputting a relative positioning result. The method not only can overcome the influence caused by position increment and GNSS original observation data broadcasting delay, but also can provide a precise relative positioning result with extremely high updating rate through low data broadcasting rate and sampling rate。

Description

Satellite and inertia combined dynamic-alignment real-time precise relative positioning method
Technical Field
The invention belongs to the technical field of satellite and inertia combined relative positioning navigation, and particularly relates to a satellite and inertia combined dynamic and dynamic precise relative positioning method with low data broadcasting rate, low sampling rate and extremely high updating rate.
Background
High update rate precision relative positioning is a technology required in many safety-related mobile-to-mobile applications, such as vehicle-to-vehicle cooperative safety applications, autonomous air refueling, and carrier-based aircraft landing. The currently common dynamic-to-dynamic relative positioning technique is a carrier-phase RTK (Real-Time Kinematic) differential relative positioning technique. However, the realization of high update rate output by using a carrier-phase RTK (Real-Time Kinematic) differential relative positioning technology mainly faces the problems of large data broadcasting rate, high requirement on the sampling rate of a receiver, difficulty in synchronizing satellite-guided observation data caused by broadcasting delay or communication data packet loss and the like.
To overcome the above disadvantages, Inertial Navigation System (INS) measurement information may be combined with Global Navigation Satellite System (GNSS) raw observation information to provide high update rate relative positioning results. There are currently documents reporting satellite/inertial combination Relative Positioning methods (see [1] Williamson, W. R.; Abdel-Hafez, M. F.; Rhee, I.; Song, E.; Wolfe, J. D.; Chichka, D. F.; Speyer, J. L., An Instrumentation System Applied to Formation flight. IEEE Transactions on Control Systems 2007, 15, (1), 75-85 [2] Alam, N.; Kelly, A. Dempster, A. G., An INS-air injection timing acquisition for Relative Positioning Systems, 14, J. H.; H. J. 1996; H. J. for alignment Systems, K. E. C. Adaptive GPS/INS integration for relative navigation. GPS Solutions 2016, 20, (1), 63-75.). However, these methods all need to broadcast the original inertial navigation measurement information, do not consider the problem of communication pressure and calculation burden caused by the data broadcast rate and the sampling rate, and do not consider the possibility of realizing a centimeter-level relative positioning solution with a very high update rate by using a low data broadcast rate and a low sampling rate.
Disclosure of Invention
Aiming at the defects in the prior art, the invention provides a satellite and inertia combined dynamic-dynamic real-time precise relative positioning method. The invention relates to a satellite and inertial navigation combined dynamic precise relative positioning method based on a carrier phase real-time kinematic (RTK) differential relative positioning method (DGNSS), a carrier phase GNSS/inertial navigation combined (TC-GNSS/INS) algorithm and a position increment polynomial forecasting method. The method can not only overcome the influence caused by position increment and GNSS original observation data broadcasting delay, but also provide a precise relative positioning result with extremely high updating rate through low data broadcasting rate and sampling rate.
In order to achieve the technical purpose, the technical scheme of the invention is as follows:
the real-time precise relative positioning method of the satellite and the inertia combined dynamic-alignment comprises the following steps:
the method comprises the steps that firstly, a mobile station receiver and a movable reference receiver synchronously sample to obtain original GNSS observation information, carrier phase cycle slip and pseudo range gross error in the original GNSS observation information obtained by respective sampling are respectively detected and eliminated, and carrier phase and pseudo range observation information of the mobile station receiver and the movable reference receiver after preprocessing are obtained.
Second, determining the carrier phase integer ambiguity
Figure 754053DEST_PATH_IMAGE001
Relative positioning is carried out by utilizing a carrier phase real-time dynamic differential relative positioning method to obtain a relative positioning vector between the mobile station and the moving referencedX RB,DGNSS 。
Third, calculating the measurement update time interval deltat TC Mobile station position increment in
Figure 835142DEST_PATH_IMAGE002
Moving reference position increment
Figure 529559DEST_PATH_IMAGE003
At the same time, the forecast time interval Delta of the mobile station is calculatedt p Increment of position within
Figure 823138DEST_PATH_IMAGE004
。
Fourthly, polynomial forecasting method based on sliding windowMethod for obtaining dynamic reference at forecast time interval Deltat p Increment of position within
Figure 535879DEST_PATH_IMAGE005
。
The fifth step is based ondX RB,DGNSS 、
Figure 725551DEST_PATH_IMAGE006
、
Figure 156533DEST_PATH_IMAGE007
、
Figure 4534DEST_PATH_IMAGE008
And
Figure 837361DEST_PATH_IMAGE009
and synthesizing and outputting the relative positioning result. In particular, the relative positioning results are synthesizeddX RB,DGNSS/INS Comprises the following steps:
Figure 994673DEST_PATH_IMAGE010
。
further, in the first step of the present invention, the fault information includes carrier phase cycle slip and pseudorange gross error. For carrier phase cycle slip, a GF method is adopted to be combined with inertial navigation auxiliary cycle slip detection to simultaneously detect and eliminate single-frequency or double-frequency carrier phase cycle slip of different satellites; and for pseudo range gross errors, rejecting the pseudo range gross errors after detection by adopting a single-station and double-station inertial navigation auxiliary gross error detection method.
Further, in the second step of the present invention, the carrier phase integer ambiguity is determined first
Figure 912950DEST_PATH_IMAGE011
Whether or not it is known. If carrier phase integer ambiguity
Figure 299063DEST_PATH_IMAGE012
If known, the carrier phase real-time dynamic differential relative positioning method is directly utilized to carry out relative positioning to obtain movementRelative positioning vector between station and moving referencedX RB,DGNSS . If carrier phase integer ambiguity
Figure 189659DEST_PATH_IMAGE013
If unknown, The carrier phase integer ambiguity resolution is solved using The least squares ambiguity reduction correlation adjustment method (LAMBDA, see Teunessen, P.J. G., The least-square ambiguity resolution adaptation: a method for The fast GPS inter-ambiguity resolution. Journal of geodety 1995, 70, (1), 65-82.) well known in The art based on The preprocessed carrier phase and pseudorange observations of The mobile station receiver and The mobile reference receiver
Figure 517872DEST_PATH_IMAGE014
Then, relative positioning is carried out by utilizing a carrier phase real-time dynamic differential relative positioning method to obtain a relative positioning vector between the mobile station and the moving referencedX RB,DGNSS 。
In particular, relative positioning vectors between a mobile station and a moving referencedX RB,DGNSS The acquisition method comprises the following steps:
firstly, the phase of the double-difference carrier between the stations is expressed as:
Figure 657867DEST_PATH_IMAGE015
wherein the DD operator expression is
Figure 362517DEST_PATH_IMAGE016
I.e. by
Figure 652160DEST_PATH_IMAGE017
，
Figure 151274DEST_PATH_IMAGE018
，
Figure 778565DEST_PATH_IMAGE019
，
Figure 21327DEST_PATH_IMAGE020
；
Figure 433985DEST_PATH_IMAGE021
、
Figure 104001DEST_PATH_IMAGE022
Respectively representing mobile stationsRObserved satelliteiAndjthe carrier phase observation of (a);
Figure 421850DEST_PATH_IMAGE023
、
Figure 202724DEST_PATH_IMAGE024
respectively representing moving referencesBObserved satelliteiAndjthe carrier phase observation of (a);
Figure 984735DEST_PATH_IMAGE025
、
Figure 576385DEST_PATH_IMAGE026
respectively representing mobile stationsRObserved satelliteiAndjthe distance between the star and the ground;
Figure 178267DEST_PATH_IMAGE027
、
Figure 762832DEST_PATH_IMAGE028
respectively representing moving referencesBObserved satelliteiAndjthe distance between the star and the ground;
Figure 133771DEST_PATH_IMAGE029
、
Figure 896322DEST_PATH_IMAGE030
respectively representing mobile stationsRObserved satelliteiAndjcarrier phase ambiguity of (a);
Figure 719921DEST_PATH_IMAGE031
、
Figure 311440DEST_PATH_IMAGE032
respectively representing moving referencesBObserved satelliteiAndjcarrier phase ambiguity of (a);
Figure 536884DEST_PATH_IMAGE033
、
Figure 985183DEST_PATH_IMAGE034
respectively representing mobile stationsRObserved satelliteiAndjthe measurement error of (2);
Figure 315320DEST_PATH_IMAGE035
、
Figure 241688DEST_PATH_IMAGE036
respectively representing moving referencesBObserved satelliteiAndjthe measurement error of (2).
Ignoring the inter-station sight vector change, the observation equation after linearization is expressed as:
Figure 321639DEST_PATH_IMAGE037
wherein
Figure 144102DEST_PATH_IMAGE038
And
Figure 942294DEST_PATH_IMAGE039
respectively representing mobile stationsRObserved satellite
Figure 157506DEST_PATH_IMAGE040
And
Figure 357543DEST_PATH_IMAGE041
the line-of-sight vector of (a),b RB representing the baseline vector to be found, i.e. the relative positioning vector between the mobile station and the moving reference.
Given at least 4 co-visitors, the relative position vector between the mobile station and the moving referenceb RB Obtained by a least squares method, and obtaining a relative position vector between the mobile station and the moving referenceb RB Is marked asdX RB,DGNSS . In which carrier phase ambiguity
Figure 944382DEST_PATH_IMAGE042
、
Figure 167553DEST_PATH_IMAGE043
Based on The double-difference carrier phase and The pseudo-range, The method of Lambda (see, e.g., Tennissen, P.J. G., The least-square algorithm calibration adaptation: a method for fast GPS integer estimation. Journal of geodety 1995, 70, (1), 65-82.) is used for calculation.
Further, the third step of the method is realized as follows:
firstly, determining a state vector, wherein the state vector comprises 9 navigation error states and 6 sensor deviation states, and is expressed as follows:
Figure 452035DEST_PATH_IMAGE044
wherein
Figure 506578DEST_PATH_IMAGE045
、
Figure 202002DEST_PATH_IMAGE046
、
Figure 240365DEST_PATH_IMAGE047
、
Figure 62959DEST_PATH_IMAGE048
And
Figure 972009DEST_PATH_IMAGE049
respectively representkThe position, the speed, the attitude error vector, the addition table and the gyro zero offset error of the epoch.
The use is carried out by a first stepConstructing an observation equation by the GNSS original observation information after the obstacle detection and the elimination, and then, obtaining the satelliteiIn thatkObserved information vector after epoch linearization
Figure 103913DEST_PATH_IMAGE050
Expressed as:
Figure 301676DEST_PATH_IMAGE051
whereinrA reference satellite number is indicated and,
Figure 177228DEST_PATH_IMAGE052
a time difference operator is represented by a time difference operator,
Figure 954167DEST_PATH_IMAGE053
after representing linearizationkSatellite to be detected of epochiAnd a reference satelliterThe single-differenced pseudoranges between the satellites,
Figure 256972DEST_PATH_IMAGE054
after linearizationkThe time difference term of the inter-satellite single difference carrier phase of the epoch.
Suppose the total number of satellites to be detected ism，kThe observation matrix of the epoch is:
Figure 4348DEST_PATH_IMAGE055
the time update and measurement update equations are expressed as:
Figure 621274DEST_PATH_IMAGE056
wherein
Figure 504917DEST_PATH_IMAGE057
To representk-1 epoch tokThe state transition matrix of the epoch is,W k-1to representk-1 epoch tokProcess noise moment of epochThe number of the arrays is determined,H k to representkA design matrix of the epoch is used,V k to representkThe measured noise matrix of the epoch is,X k-1to representk-1 epoch state vector.
Suppose thatPExpressing the variance matrix, performing Kalman filtering, to obtainkTo k+△t TC Epoch, measurement update time interval is Δt TC State vector delta of
Figure 729356DEST_PATH_IMAGE058
And its corresponding variance matrix
Figure 698449DEST_PATH_IMAGE059
The expression of (a) is as follows:
Figure 119066DEST_PATH_IMAGE060
Figure 857215DEST_PATH_IMAGE061
wherein
Figure 767402DEST_PATH_IMAGE062
Figure 427054DEST_PATH_IMAGE063
Figure 198832DEST_PATH_IMAGE064
Figure 57066DEST_PATH_IMAGE065
Figure 75838DEST_PATH_IMAGE066
Figure 19523DEST_PATH_IMAGE067
To representkEpoch tok+△t TC The state transition matrix of the epoch is,Ithe unit matrix is represented by a matrix of units,
Figure 578680DEST_PATH_IMAGE068
to representk+△t TC The matrix of the variance of the epoch is,Q k-1to representk-1 epoch tokThe process noise covariance matrix of the epoch,R k to representkThe measured noise covariance matrix of the epoch.
By using
Figure 494684DEST_PATH_IMAGE069
And
Figure 231827DEST_PATH_IMAGE070
is calculated at the mobile station and the moving reference respectively to obtain the mobile stationkTok+△t TC Epoch state vector delta
Figure 928387DEST_PATH_IMAGE071
Dynamic referencekTok+△t TC Epoch state vector delta
Figure 228919DEST_PATH_IMAGE072
And mobile stationk+△t TC Tok+△t TC +△t p Epoch state vector delta
Figure 796166DEST_PATH_IMAGE073
. The calculation expression is
Figure 219057DEST_PATH_IMAGE074
In the formula
Figure 613699DEST_PATH_IMAGE075
And
Figure 717921DEST_PATH_IMAGE076
indicating mobile station and moving reference, respectivelykEpoch tok+△t TC The state transition matrix of the epoch is,
Figure 405254DEST_PATH_IMAGE077
indicating a mobile stationk+△t TC To k+△t TC +△t p The state transition matrix of the epoch is,
Figure 999047DEST_PATH_IMAGE078
and
Figure 342303DEST_PATH_IMAGE079
indicating mobile station and moving reference, respectivelykTok+△t TC The process noise of the epoch is that of the epoch,
Figure 781375DEST_PATH_IMAGE080
indicating a mobile stationk+△t TC To k+△t TC +△t p The process noise of the epoch is that of the epoch,
Figure 339526DEST_PATH_IMAGE081
、
Figure 41903DEST_PATH_IMAGE082
and
Figure 669194DEST_PATH_IMAGE083
and the above
Figure 911956DEST_PATH_IMAGE084
In the same way, the first and second,Irepresenting an identity matrix. Based on the state vector increment calculated above, respectively extracting
Figure 839461DEST_PATH_IMAGE085
、
Figure 260209DEST_PATH_IMAGE086
And
Figure 312479DEST_PATH_IMAGE087
the first three-dimensional vector is the matrix measurement update time interval Δ to be obtainedt TC Mobile stationkTok+△t TC Epoch location increment
Figure 358932DEST_PATH_IMAGE088
Dynamic referencekTok+△t TC Epoch location increment
Figure 140943DEST_PATH_IMAGE089
And the mobile station in the forecast time interval deltat p In the interior of said container body,k+△t TC to k+△t TC +△t p Epoch location increment
Figure 185123DEST_PATH_IMAGE090
。
Further, in the sliding window-based polynomial forecasting method in the fourth step of the present invention, the one-dimensional direction polynomial model is defined as:
Figure 521426DEST_PATH_IMAGE091
wherein
Figure 122303DEST_PATH_IMAGE092
The variables to be modeled are represented by a graph,T i indicating the relative time with respect to the start of the sliding window,othe order of the order is represented,
Figure 493241DEST_PATH_IMAGE093
、
Figure 708322DEST_PATH_IMAGE094
...
Figure 63080DEST_PATH_IMAGE095
the coefficient of the polynomial is represented by,
Figure 185757DEST_PATH_IMAGE096
andT i the expression is:
Figure 880043DEST_PATH_IMAGE097
whereint real1,The true tail time representing the first position increment of the sliding window,t i real,indicating a sliding windowiThe true tail time of a position increment,
Figure 810566DEST_PATH_IMAGE098
to representt i real,The reference station corresponding to the moment is in the measurement update time interval deltat TC Increment of position within.
When the number of elements of the three-dimensional sliding window of the position vector is larger than that of the elements of the three-dimensional sliding window of the position vectoro+1, all polynomial coefficients can be solved using the least squares method. Suppose thatk+△t TC The relative time of the epoch with respect to the start of the sliding window ist c,real The moving reference is then at the forecast time interval Δt p Corresponding relative time ist c,real +△t p At this time, the dynamic reference position is increased
Figure 387040DEST_PATH_IMAGE099
The three-dimensional positions are all calculated by the following formula:
Figure 313408DEST_PATH_IMAGE100
in the formula
Figure 862201DEST_PATH_IMAGE101
、
Figure 481401DEST_PATH_IMAGE102
...
Figure 30325DEST_PATH_IMAGE103
、
Figure 963646DEST_PATH_IMAGE104
、
Figure 163684DEST_PATH_IMAGE105
...
Figure 688206DEST_PATH_IMAGE106
And
Figure 239273DEST_PATH_IMAGE107
、
Figure 523755DEST_PATH_IMAGE108
...
Figure 578299DEST_PATH_IMAGE109
respectively represent
Figure 476984DEST_PATH_IMAGE110
Is/are as followsx、yAndzthe polynomial coefficients of the three-dimensional directions,
Figure 249768DEST_PATH_IMAGE111
、
Figure 587209DEST_PATH_IMAGE112
and
Figure 246991DEST_PATH_IMAGE113
respectively represent
Figure 582158DEST_PATH_IMAGE114
Is/are as followsx、yAndzthree-dimensional directions at time intervals Δt p The internal dynamic reference position increment forecast value.
Compared with the prior art, the invention has the following advantages:
the method provided by the invention can output centimeter-level relative positions with extremely high update rate (100 Hz) equal to the sampling rate of inertial measurement information for dynamic-dynamic application.
The invention can provide precise relative positioning results with extremely high update rate through low data broadcasting rate and sampling rate.
The invention can overcome the influence of communication delay on real-time performance and has certain capability of tolerating data loss.
Drawings
FIG. 1 is a general block diagram of the present invention;
FIG. 2 is a flow chart of single station fault detection and rejection in accordance with the present invention;
FIG. 3 is a flow chart of the dual station fault detection and rejection of the present invention;
FIG. 4 is a schematic diagram of the combination of time and relative position at 103672.113s according to the present invention;
FIG. 5 is a diagram of the horizontal relative positioning error of a dynamic test of the present invention;
FIG. 6 is a vertical relative positioning error plot for a dynamic test of the present invention.
Detailed Description
In order to make the technical scheme and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
Example 1:
referring to fig. 1, a satellite and inertia combined dynamic-dynamic real-time precise relative positioning method includes the following steps:
the method comprises the steps that firstly, a mobile station receiver and a mobile reference receiver synchronously sample to obtain original GNSS observation information, and fault information gross errors of carrier phase cycle slip and pseudo-range observation information in the original GNSS observation information obtained by sampling are respectively detected and removed to obtain the carrier phase and pseudo-range observation information of the mobile station receiver and the mobile reference receiver after preprocessing.
The invention adopts a multi-method combined fault detection and elimination. GNSS observation data gross errors must be monitored and culled before being used for relative positioning. The fault information includes carrier phase cycle slip and pseudorange gross error. Carrier phase cycle slip and pseudorange gross errors are the most common faults.
For carrier phase cycle slip, the classical GF method is combined with inertial navigation assisted cycle slip detection to simultaneously detect and reject single or dual frequency cycle slips of different satellites. The GF test statistic and its variance can be expressed as
Figure 576659DEST_PATH_IMAGE115
（1）
WhereiniRepresenting the number of the satellite to be detected;f 1andf 2representing different frequency points of the signal and,
Figure 717790DEST_PATH_IMAGE116
and
Figure 481347DEST_PATH_IMAGE117
respectively representf 1And f 2the wavelength of the frequency point; delta t Representing a time difference operator;
Figure 537814DEST_PATH_IMAGE118
and
Figure 222873DEST_PATH_IMAGE119
respectively representing satellites
Figure 636537DEST_PATH_IMAGE120
Frequency pointf 1Andf 2is detected by the carrier-phase observation of (c),
Figure 520179DEST_PATH_IMAGE121
and
Figure 993886DEST_PATH_IMAGE122
respectively represent
Figure 979291DEST_PATH_IMAGE118
And
Figure 196645DEST_PATH_IMAGE119
error variance after time differentiation.
Figure 138057DEST_PATH_IMAGE118
And
Figure 782665DEST_PATH_IMAGE119
it may be in the form of non-differences, single differences or double differences. Non-difference means that the raw measurements are not combined. The single difference is divided into two forms of single difference between stars and single difference between stations. Single inter-satellite differences represent the differential combination of carrier phase measurements from different satellites at the same time. Single difference between stations means that carrier phase measurements received by different stations at the same time are differentially combined. The double differences represent that the single difference carrier phase measurement values between stations of different satellites at the same time are further differentially combined on the basis of the single difference between the stations. And (3) respectively substituting the carrier phase observed values in the four combination forms into formula (1) to calculate the test statistic and the variance thereof, and further obtaining a test threshold. Comparing the test statistic to a test threshold: if the absolute value of the test statistic is less than or equal to the test threshold, the test is passed without failure; if the absolute value of the test statistic is larger than the test threshold, the test fails and the fault occurs. This results in four forms of GF methods for cycle slip detection: non-differential GF method, inter-satellite single-difference GF method, inter-station single-difference GF method and inter-station double-difference GF method. The GF method has good performance for detecting small cycle slips, but cannot detect special dual-frequency cycle slips. In addition, this method also cannot distinguish between single frequency cycle slips and requires a dual frequency receiver.
The inertial navigation auxiliary cycle slip detection method can overcome the defects of the GF method, and the test statistic and the variance expression thereof of the method are
Figure 504633DEST_PATH_IMAGE123
（2）
WhereinrRepresents a reference satellite number;vrepresents a receiver number;crepresents the speed of light;frepresenting a frequency point;
Figure 463362DEST_PATH_IMAGE124
to representfThe wavelength of the frequency point;
Figure 72329DEST_PATH_IMAGE125
representing satellitesiAndrsingle difference satellite clock difference between the stars;
Figure 887838DEST_PATH_IMAGE126
presentation receivervSatellite to be detectediAnd a reference satelliterSingle difference between starsfA frequency point carrier phase observed value;
Figure 34786DEST_PATH_IMAGE127
to represent
Figure 593943DEST_PATH_IMAGE128
The measurement variance after time differentiation;
Figure 572263DEST_PATH_IMAGE129
presentation receivervMiddle satelliteiAndrsingle difference satellite-to-ground distance between the satellites;
Figure 309406DEST_PATH_IMAGE130
to represent
Figure 678071DEST_PATH_IMAGE131
Error variance after time differentiation.
Figure 40919DEST_PATH_IMAGE132
Can be expressed as:
Figure 873746DEST_PATH_IMAGE133
（3）
whereind i Representing satellitesiThe position of (a);d r representing reference satellitesrThe position of (a);
Figure 234320DEST_PATH_IMAGE134
indicating INS provided receivervThe position of (a). The method calculates a test threshold using the variance, and then compares the test statistic with the test threshold: if the absolute value of the test statistic is less than or equal to the test threshold, the test is passed without failure; if the absolute value of the test statistic is larger than the test threshold, the test fails and the fault occurs. For convenience of expression, the cycle slip detection method composed of the test statistic and the variance thereof is abbreviated as an inter-satellite single-difference inertial auxiliary cycle slip detection method.
When the observation information of the moving reference and the absolute position are received by the mobile station, the observation information of the two-station combination forms a test statistic and the variance thereof is
Figure 887018DEST_PATH_IMAGE135
（4）
Wherein
Figure 801360DEST_PATH_IMAGE136
Receiver with a plurality of receiversvAndwof a satelliteiAndrthe inter-station inter-satellite double-difference satellite-ground distance is expressed as follows:
Figure 488693DEST_PATH_IMAGE137
wherein
Figure 20169DEST_PATH_IMAGE138
Indicating reception via a data chainwIn the position of (a) in the first,
Figure 425742DEST_PATH_IMAGE139
presentation receivervAndwof a satelliteiAndrinter-station double difference between starsfThe observed value of the carrier phase of the frequency point,
Figure 864814DEST_PATH_IMAGE140
to represent
Figure 609916DEST_PATH_IMAGE141
The variance of the measurement after the time difference,
Figure 859763DEST_PATH_IMAGE142
to represent
Figure 752632DEST_PATH_IMAGE143
The error variance after the time difference is defined as described above. For convenience of expression, the cycle slip detection method based on the test statistic and the variance thereof is abbreviated as an interstation inter-satellite double-difference inertial auxiliary cycle slip detection method. Various types of cycle slips can be detected and eliminated by combining a GF method and an inertial navigation auxiliary method.
For pseudorange gross error, considering that the performance of the pure satellite navigation method is poor, the embodiment only adopts the inertial navigation assisted gross error detection method of the single station and the double station. The single station test statistic and its variance are
Figure 198657DEST_PATH_IMAGE144
（5）
WhereiniIndicating the number of the satellite to be detected,ra reference satellite number is indicated and,vrepresents a receiver number;cthe speed of light is indicated and is,fthe frequency points of the signals are represented,
Figure 595004DEST_PATH_IMAGE145
presentation receiver
Figure 530599DEST_PATH_IMAGE146
Middle satelliteiAnd rsingle difference between starsfThe observed value of the pseudo-range of the frequency point,
Figure 395917DEST_PATH_IMAGE147
representing satellitesiAndrthe clock error of the single-difference satellite between the satellites,
Figure 645633DEST_PATH_IMAGE148
representation receivervMiddle satelliteiAnd rthe single inter-satellite-to-ground distance, as defined above,
Figure 162065DEST_PATH_IMAGE149
presentation receivervMiddle satelliteiAnd rsingle difference between starsfIonospheric delay at the frequency points, cancellation using a double difference deionization layer combination,
Figure 2982DEST_PATH_IMAGE150
presentation receivervMiddle satelliteiAndrsingle difference between starsfTropospheric delay at the frequency points, compensated using the Saastamoinen model,to represent
Figure 392692DEST_PATH_IMAGE152
The error variance after the time difference, as defined above,
Figure 514363DEST_PATH_IMAGE153
to represent
Figure 791761DEST_PATH_IMAGE154
The measured variance after time differentiation. For convenience of expression, the cycle slip detection method based on the test statistic and the variance thereof is abbreviated as an inter-satellite single-difference inertial-assisted gross error detection method.
Under short baseline conditions, ionospheric and tropospheric delays can be assumed to be eliminated by double-difference combining, so that the two-station test statistic and its variance are
Figure 880940DEST_PATH_IMAGE155
（6）
In the formula
Figure 3616DEST_PATH_IMAGE156
Presentation watchDisplay receivervAndwof a satelliteiAnd rinter-station double difference between starsfThe observed value of the pseudo-range of the frequency point,
Figure 432324DEST_PATH_IMAGE157
presentation receivervAndwof a satelliteiAnd rthe inter-station inter-satellite double-difference inter-satellite distance between the stations is the same as the definition,
Figure 646003DEST_PATH_IMAGE158
to represent
Figure 222478DEST_PATH_IMAGE159
The variance of the error is determined by the error variance,
Figure 352108DEST_PATH_IMAGE160
to represent
Figure 697639DEST_PATH_IMAGE161
Error variance, as defined above. For convenience of expression, the cycle slip detection method based on the test statistic and the variance thereof is abbreviated as an interstation intersatellite double-difference inertial auxiliary gross error detection method.
All of the above test statistics obey a 0-mean global distribution with known variance in the original hypothesis. In an alternative assumption, the mean is equal to the gross error and the variance remains unchanged. Thus, the required false alarm rate is givenP FA The check threshold can be obtained as
Figure 316839DEST_PATH_IMAGE162
（7）
Wherein sigma T Represents the standard deviation of the test statistic Φ-1(x) Denotes Φ: (x) An inverse function, which is defined as
Figure 865763DEST_PATH_IMAGE163
（8）
All of the above methods may be combined to detect carrier phase cycle slip and pseudorange gross for single and dual station combined observation information.
In order to fully exploit the advantages of each method, the present embodiment designs a joint fault detection rule for single-station and dual-station, respectively, as shown in fig. 2 and 3.
The single-station joint fault detection steps in FIG. 2 are as follows:
firstly, the cycle slip of the carrier phase observation information is preliminarily detected by using a non-difference GF method.
Selecting the satellite with the largest elevation angle from the satellites without cycle slip after the non-difference GF detection as a reference satellite.
And thirdly, detecting the cycle slip of the carrier phase observation information by utilizing an inter-satellite single-difference inertial auxiliary cycle slip detection method based on the position information predicted by the inertial navigation system.
And fourthly, detecting the cycle slip of the carrier phase information by using a single difference GF method between the satellites.
And integrating the cycle slip detection results of the inter-satellite single-difference inertial auxiliary cycle slip detection method and the inter-satellite single-difference GF method. The two methods are judged whether the test passes: if yes, no cycle slip exists; if not, the satellite observation information is marked as invalid if the satellite observation information has cycle slip.
Sixthly, detecting the gross error of the pseudo-range observation information by using an inter-satellite single-difference inertial-assisted gross error detection method based on the position information forecasted by the inertial navigation system. Judging whether the test passes: if yes, no gross error exists; if not, the satellite observation information is marked as invalid if the satellite observation information has gross errors.
The double-station joint fault detection steps in fig. 3 are as follows:
firstly, preliminarily detecting cycle slip of carrier phase observation information by using an inter-station single difference GF method.
Selecting the satellite with the largest elevation angle from the satellites without cycle slip after the interstation single difference GF method detection as a reference satellite.
And thirdly, detecting the cycle slip of the carrier phase observation information by using an inter-station inter-satellite double-difference inertial auxiliary cycle slip detection method based on the position information predicted by the inertial navigation system.
And fourthly, detecting the cycle slip of the carrier phase information by using a GF method of the double difference between the satellites of the stations.
And integrating the cycle slip detection results of the inter-station inter-satellite double-difference inertial auxiliary cycle slip detection method and the inter-station inter-satellite double-difference GF method. The two methods are judged whether the test passes: if yes, no cycle slip exists; if not, the satellite observation information is marked as invalid if the satellite observation information has cycle slip.
And sixthly, detecting the gross error of the pseudo-range observation information by using an inter-station inter-satellite double-difference inertial auxiliary gross error detection method based on the position information predicted by the inertial navigation system. Judging whether the test passes: if yes, no gross error exists; if not, the satellite observation information is marked as invalid if the satellite observation information has gross errors.
Second, determining the carrier phase integer ambiguity
Figure 799084DEST_PATH_IMAGE164
Whether or not it is known. If not, based on the preprocessed carrier phase and pseudo-range observation information of the mobile station receiver and the mobile reference receiver, solving the carrier phase integer ambiguity by using a least squares ambiguity reduction correlation adjustment method (LAMBDA) known in the art
Figure 999121DEST_PATH_IMAGE164
Then, relative positioning is carried out by utilizing a carrier phase real-time dynamic differential relative positioning method to obtain a relative positioning vector between the mobile station and the moving referencedX RB,DGNSS (ii) a If yes, the relative positioning vector between the mobile station and the moving reference is obtained by directly utilizing the carrier phase real-time dynamic differential relative positioning method to carry out relative positioningdX RB,DGNSS 。
Because the ionospheric and tropospheric errors are highly correlated at the short baseline (< 10 km) between the mobile station and the mobile reference, other common errors such as receiver clock error, satellite clock error, ionospheric delay, tropospheric delay and the like can be eliminated by using the inter-station inter-satellite double-difference carrier phase.
The phase of the double-difference carrier wave between the stations is expressed as follows:
Figure 789223DEST_PATH_IMAGE165
（9）
wherein the DD operator expression is
Figure 74711DEST_PATH_IMAGE166
I.e. by
Figure 546143DEST_PATH_IMAGE167
，
Figure 351419DEST_PATH_IMAGE168
，
Figure 578001DEST_PATH_IMAGE169
，
Figure 350785DEST_PATH_IMAGE170
；φRepresenting an observation of the carrier phase,ρthe distance between the star and the ground is shown,λwhich is indicative of the wavelength of the signal,Nwhich is indicative of the carrier phase ambiguity,indicating measurement error, subscriptRAndBindicating mobile station and moving reference, respectively, superscriptiAndjindicating the number of satellites to be probed.
Ignoring the inter-station line-of-sight vector variation, the observation equation after linearization can be expressed as
Figure 360330DEST_PATH_IMAGE171
（10）
Wherein
Figure 534959DEST_PATH_IMAGE172
And
Figure 417596DEST_PATH_IMAGE173
respectively representing satellitesiAndjthe line-of-sight vector of (a),b RB representing the baseline vector to be found, i.e. the relative positioning vector between the mobile station and the moving reference. The carrier phase ambiguity can be computed based on The double-differenced carrier phase and pseudorange using The LAMBDA method known in The art (see Teunessen, P.J. G., The least-square ambiguity correction: a method for fast GPS inter-ambiguity, Journal of geodety 1995, 70, (1), 65-82.).
Given at least 4 co-visitors, the relative position vector between the mobile station and the moving referenceb RB Can be obtained by Least Squares (LSE) method, and the relative position vector between the mobile station and the moving reference is obtainedb RB Is marked asdX RB,DGNSS 。
Thirdly, calculating a measurement update time interval by utilizing a TC-GNSS/INS algorithm
Figure 677676DEST_PATH_IMAGE174
Mobile station position increment in
Figure 490911DEST_PATH_IMAGE175
And moving reference position increment
Figure 520047DEST_PATH_IMAGE176
. Meanwhile, the TC-GNSS/INS algorithm is used for calculating the forecast time interval of the mobile station
Figure 88431DEST_PATH_IMAGE177
Increment of position within
Figure 318031DEST_PATH_IMAGE178
。
A state vector is first determined.
To improve the accuracy of the position increment, carrier phase and pseudorange observations are used in the measurement equations. The state vector contains 9 navigation error states and 6 sensor bias states, which can be expressed as
Figure 731695DEST_PATH_IMAGE179
（11）
Wherein
Figure 818599DEST_PATH_IMAGE180
、
Figure 292306DEST_PATH_IMAGE181
、
Figure 261399DEST_PATH_IMAGE182
、
Figure 963907DEST_PATH_IMAGE183
And k respectively representkThe position, the speed, the attitude error vector, the addition table and the gyro zero offset error of the epoch.
The GNSS original observation information used after the first step of fault detection and elimination is used for constructing an observation equation, and then the satelliteiIn thatkObserved information vector after epoch linearization
Figure 233214DEST_PATH_IMAGE184
Can be expressed as:
Figure 877822DEST_PATH_IMAGE185
（12）
whereinrDenotes the reference satellite number, Δ t A time difference operator is represented by a time difference operator,
Figure 537474DEST_PATH_IMAGE186
after representing linearizationkSatellite to be detected of epochiAnd a reference satelliterThe single-differenced pseudoranges between the satellites,
Figure 292940DEST_PATH_IMAGE187
after linearizationkThe time difference term of the inter-satellite single difference carrier phase of the epoch.
For carrier phase, inter-epoch differences are used to remove the ambiguity parameters. For pseudoranges, ionospheric delays are cancelled using a dual-frequency deionization layer combination, and tropospheric delays are compensated using a Saastamoinen model.
Suppose the total number of satellites to be detected ism，kThe observation matrix of the epoch is:
Figure 901907DEST_PATH_IMAGE188
（13）
the time update and measurement update equations are expressed as:
Figure 717416DEST_PATH_IMAGE189
（14）
wherein
Figure 864364DEST_PATH_IMAGE190
To representk-1 epoch tokThe state transition matrix of the epoch is,W k-1to representk-1 epoch tokThe course of the epoch noise matrix is,H k to representkA design matrix of the epoch is used,V k to representkThe measured noise matrix of the epoch is,X k-1to representk-1 epoch state vector.
Suppose thatPRepresenting a variance matrix, a kalman filtering algorithm may be performed. FromkTo k+△t TC Epoch, measurement update time interval is Δt TC State vector delta of
Figure 423521DEST_PATH_IMAGE191
And its corresponding variance matrix
Figure 401842DEST_PATH_IMAGE192
The expression of (a) is as follows:
Figure 873405DEST_PATH_IMAGE193
（15）
Figure 507649DEST_PATH_IMAGE194
（16）
wherein
Figure 870497DEST_PATH_IMAGE195
（17）
Figure 703324DEST_PATH_IMAGE196
（18）
Figure 595057DEST_PATH_IMAGE197
（19）
Figure 716596DEST_PATH_IMAGE198
（20）
Figure 371218DEST_PATH_IMAGE199
（22）
In the formula
Figure 58552DEST_PATH_IMAGE200
To representkEpoch tok+△t TC The state transition matrix of the epoch is,Ithe unit matrix is represented by a matrix of units,
Figure 386765DEST_PATH_IMAGE201
to representk+△t TC The matrix of the variance of the epoch is,Q k-1to representk-1 epoch tokThe process noise covariance matrix of the epoch,R k to representkThe measured noise covariance matrix of the epoch.
By using
Figure 730022DEST_PATH_IMAGE202
And
Figure 700252DEST_PATH_IMAGE203
is calculated at the mobile station and the moving reference respectively to obtain the mobile stationkTok+△t TC Epoch state vector delta
Figure 992824DEST_PATH_IMAGE204
Dynamic referencekTok+△t TC Epoch state vector delta
Figure 695201DEST_PATH_IMAGE205
And mobile stationk+△t TC Tok+△t TC +△t p Epoch state vector delta
Figure 322491DEST_PATH_IMAGE206
. The calculation expression is as follows:
Figure 565253DEST_PATH_IMAGE207
（23）
in the formula
Figure 227179DEST_PATH_IMAGE208
And
Figure 100457DEST_PATH_IMAGE209
indicating mobile station and moving reference, respectivelykEpoch tok+△t TC The state transition matrix of the epoch is,
Figure 965776DEST_PATH_IMAGE210
indicating a mobile stationk+△t TC Tok+△t TC +△t p The state transition matrix of the epoch is,W R,K and W B,K indicating mobile station and moving reference, respectivelykTok+△t TC The process noise of the epoch is that of the epoch,
Figure 12229DEST_PATH_IMAGE211
indicating a mobile stationk+△t TC Tok+△t TC +△t p The process noise of the epoch is that of the epoch,
Figure 731924DEST_PATH_IMAGE212
、
Figure 572841DEST_PATH_IMAGE213
and
Figure 440303DEST_PATH_IMAGE214
and the above
Figure 775600DEST_PATH_IMAGE215
In the same way, the first and second,Irepresenting an identity matrix.
State vector increment based on the above calculation
Figure 349801DEST_PATH_IMAGE216
、
Figure 361619DEST_PATH_IMAGE217
And
Figure 450798DEST_PATH_IMAGE218
respectively extract
Figure 839054DEST_PATH_IMAGE219
、
Figure 267761DEST_PATH_IMAGE220
And
Figure 198284DEST_PATH_IMAGE221
the vector formed by the first three-dimensional structure is the vector which needs to be obtained
Matrix measurement update interval Δt TC Mobile stationkTok+△t TC Epoch location increment
Figure 40338DEST_PATH_IMAGE222
Dynamic referencekTok+△t TC Epoch location increment
Figure 966706DEST_PATH_IMAGE223
And the mobile station at the forecast time interval
Figure 249919DEST_PATH_IMAGE224
In the interior of said container body,k+△t TC tok+△t TC +△t p Epoch location increment
Figure 134699DEST_PATH_IMAGE225
。
Fourthly, obtaining the forecasting time delta of the dynamic reference in a polynomial forecasting method based on a sliding windowt p Increment of position within
Figure 683623DEST_PATH_IMAGE226
。
In order to obtain the synchronous relative position of the broadcasting gap of the moving reference, a polynomial forecasting method based on a sliding window is adopted.
The one-dimensional direction polynomial model in the polynomial forecasting method is defined as:
Figure 616944DEST_PATH_IMAGE227
wherein
Figure 816981DEST_PATH_IMAGE228
Representing a variable to be modeled;T i representing the relative time with respect to the start of the sliding window;orepresents the order;
Figure 607082DEST_PATH_IMAGE229
、
Figure 95833DEST_PATH_IMAGE230
...
Figure 629582DEST_PATH_IMAGE231
the representation represents polynomial coefficients.
Figure 434858DEST_PATH_IMAGE232
And
Figure 599123DEST_PATH_IMAGE233
the expression is:
Figure 371907DEST_PATH_IMAGE234
（24）
whereint real1,The true tail time representing the first position increment of the sliding window,t i real,indicating a sliding windowiThe true tail time of a position increment,
Figure 443769DEST_PATH_IMAGE235
to representt i real,The reference station corresponding to the moment is in the measurement update time interval deltat TC Increment of position within.
When the number of elements of the three-dimensional sliding window of the position vector is larger than that of the elements of the three-dimensional sliding window of the position vectoro+1, all polynomial coefficients can be solved using the least squares method. Suppose thatk+△t TC The relative time of the epoch with respect to the start of the sliding window ist c,real The moving reference is then at the forecast time interval Δt p Corresponding relative time ist c,real +△t p At this time, the dynamic reference position is increased
Figure 821660DEST_PATH_IMAGE236
The three-dimensional positions are all calculated by the following formula:
Figure 953564DEST_PATH_IMAGE237
（25）
in the formula
Figure 964377DEST_PATH_IMAGE238
、
Figure 777612DEST_PATH_IMAGE239
、
Figure 806748DEST_PATH_IMAGE240
、
Figure 109553DEST_PATH_IMAGE241
...
Figure 60192DEST_PATH_IMAGE242
And
Figure 739435DEST_PATH_IMAGE243
、
Figure 365020DEST_PATH_IMAGE244
...
Figure 838727DEST_PATH_IMAGE245
respectively represent
Figure 11082DEST_PATH_IMAGE246
Is/are as followsx、 yAndzthe polynomial coefficients of the three-dimensional directions,
Figure 228437DEST_PATH_IMAGE247
、
Figure 232165DEST_PATH_IMAGE248
and
Figure 814456DEST_PATH_IMAGE249
respectively represent
Figure 287157DEST_PATH_IMAGE250
Is/are as followsx、 yAndzthree-dimensional directional time interval deltat p The internal dynamic reference position increment forecast value.
The fifth step, synthesize the relative positioning resultdX RB,DGNSS/INS And outputting, wherein:
Figure 42624DEST_PATH_IMAGE251
（26）
example 2:
to test the method of example 1, this example performed a vehicle-to-vehicle test, where the mobile station and the moving reference were both vehicles and had many turns while moving on the square. A static reference station is installed at a known position near a test site to calculate the post-processing relative positions of a mobile station and a moving reference. The corresponding results are used to provide reference results for the position increments and the relative position. In the experiment, a GNSS/MEMS prototype system consisting of a sensor STIM300 MEMS and a ComNavOEM-K508 board was attached to the mobile station number two mobile reference. The GNSS receiver can provide observation information of 5 frequency points for real-time and post navigation (B1/B2/B3/L1/L2). Only the observed information of four frequency points is actually used in the subsequent analysis (B1/B3/L1/L2). The GNSS receiver maximum sampling rate is 10Hz and the MEMS sampling rate is 125Hz. digital International inc. Xtend-PKG 900 MHz RF station is used to broadcast and receive data packets. The GNSS raw observation information broadcasts differential data information in MSM 5 format in RTCM 3.2. To analyze the accuracy of the position increments and the combined synchronous relative position, post-processing of the mobile and static references is relied upon, and the relative position of the mobile and static references DGNSS is used as a reference. In the test process, the sampling rates of the mobile station and the GNSS receiver of the dynamic reference are set to be 1Hz, the broadcasting rate of the original observation data of the GNSS of the dynamic reference is 1Hz, and the incremental broadcasting rate of the position of the dynamic reference is 10 Hz.
Take the 103672.113s relative position calculation for the mobile station MEMS observations as an example. At this time, the positioning device is located in the mobile station GNSS observation data sampling gap and at the same time, the mobile reference GNSS original observation data and the position increment broadcasting gap are located, so that the current dynamic-to-dynamic relative positioning result needs to be obtained by combining the RTK DGNSS relative position result at the latest moment, the position increment obtained by the close combination extrapolation and the position increment obtained by polynomial prediction, and the specific combination principle is shown in fig. 4.
(1) First computing 103672.113sdX RB,DGNSS 、
Figure 104120DEST_PATH_IMAGE252
And
Figure 185209DEST_PATH_IMAGE253
time of day using the mobile station receiver sampled 103672.0 time GNSS raw observations, the mobile reference receiver sampled 103672.0 GNSS raw observations. Before GNSS observation information is used for positioning, carrier phase cycle slip and pseudo-range gross error are detected by adopting a multi-method joint fault detection and elimination method combining single station and double stations. Through comparison, all the test statistics exceed the test threshold and the satellite has no fault.
(2) At the moment, the ambiguity of the carrier phase whole cycle is known, and based on the received 103672.0s time moving reference and GNSS original observation information of the mobile station, 103672.0s time synchronization relative position is calculated by using an RTK DGNSS relative positioning methoddX RB,DGNSS . At the moment, the total number of Beidou satellites is 7, the total number of GPS satellites is 9, two reference satellites are removed, the double-frequency carrier double-difference observed values are 28 at the moment, and calculation can be carried out to obtain
Figure 128894DEST_PATH_IMAGE254
。
(3) The dynamic reference and mobile station position increment are calculated using the TC-GNSS/INS algorithm.
The dynamic reference is calculated by using TC-GNSS/INS algorithm
Figure 625735DEST_PATH_IMAGE255
The measurement update frequency is 1Hz, the extrapolated position increment is stored and broadcast at 10Hz, and the mobile station needs to use the received position increment at time 103672.1
Figure 89208DEST_PATH_IMAGE256
At this time
Figure 341198DEST_PATH_IMAGE257
. Meanwhile, the mobile station measurement updating frequency is 1Hz, and 103672.113s is very close to the whole second moment, so that the TC-GNSS/INS algorithm can be used for directly extrapolating from the 103672.0s moment to the current moment
Figure 772179DEST_PATH_IMAGE258
I.e. directly obtained at this time
Figure 72711DEST_PATH_IMAGE259
. Known from the calculation process
Figure 639958DEST_PATH_IMAGE260
。
(4) And predicting by using a moving reference position increment polynomial.
Because the incremental broadcasting interval is positioned at the moving reference position at the moment, a polynomial forecasting method is needed to be adopted for calculation
Figure 813582DEST_PATH_IMAGE261
. In the polynomial algorithm, the order is taken as the second order, and the number of sliding window epochs is taken as 10. Based on the sliding window data, the coefficients of the polynomial in each direction can be obtained by using a least square method as follows:
Figure 935121DEST_PATH_IMAGE262
the relative time of the current time relative to the first position increment of the sliding window is 1.013s, and three time polynomial modeling parameter values can be obtained according to the polynomial model as follows:
Figure 570502DEST_PATH_IMAGE263
and forecast time
Figure 523415DEST_PATH_IMAGE264
Thus, the position increment for forecasting in three directions can be:
Figure 54890DEST_PATH_IMAGE265
(5) and synthesizing relative positioning results at a very high updating rate.
Based on the component results of the foregoing calculations, the current relative position can be synthesized as:
Figure 460464DEST_PATH_IMAGE266
the above time examples intuitively show the process of realizing the relative positioning with extremely high update rate by the method provided by the invention.
Fig. 5 and 6 show the relative positioning error in the horizontal direction and the vertical direction, respectively. FIG. 5 shows a horizontal relative positioning error scatter plot, where all points are found to be within 0.1m, and most points are found to be within 0.05 m. FIG. 6 shows the variation of the vertical error with the test time, and it can be seen that the error in the vertical direction is mostly less than 0.1 m. It can be seen from the figure that the relative positioning error can be kept in the centimeter level under the condition of low data broadcasting rate and sampling rate.
TABLE 1 relative positioning accuracy at selected data broadcast and sample rates
Figure 912917DEST_PATH_IMAGE267
Table 1 shows the relative positioning accuracy statistics. It can be seen that both the positioning accuracy in each direction and the positioning accuracy in three dimensions are very high, further proving the high accuracy performance of the method of the invention.
TABLE 2 comparison of the results of the relative positioning of the method of the invention and the 10Hz pure RTK DGNSS
Figure 658019DEST_PATH_IMAGE268
Table 2 shows the results of comparing the method of the present invention with the 10Hz pure RTK DGNSS relative positioning method. Wherein, the common 10Hz pure RTK DGNSS relative positioning method is abbreviated as method 1, and the method of the invention is abbreviated as method 2. As can be seen from the table, the data broadcasting rate and the receiver sampling rate of the method of the present invention are greatly reduced by 76.24% and 90%, respectively. And the output update rate is greatly increased by 12.5 times.
TABLE 3 statistical accuracy under minimum delay conditions for incremental broadcast at different locations
Figure 422713DEST_PATH_IMAGE269
TABLE 4 statistical accuracy under minimum delay conditions for dissemination of different GNSS original observation information
Figure 50003DEST_PATH_IMAGE270
In order to verify the real-time performance of the method, the method is analyzed according to the simulation of position increment and GNSS original observation information time delay in the recorded data in the later analysis. Both the position increment and the GNSS raw observation information set a minimum time delay.
Tables 3 and 4 show the statistical accuracy of the position increment and the GNSS original observation information under different minimum time delay conditions. It is clear that the relative position accuracy deteriorates as the minimum delay of the simulation increases. For position increments, the isotropic accuracy can be maintained in centimeters when the minimum time delay is less than or equal to 0.4 s. This means that the method proposed in mode 3 can tolerate 4 consecutive position increment packet losses. For GNSS original observation information, the attenuation towards the vertical direction and the horizontal direction is not slow. When the minimum time delay is less than or equal to 2.5s, the accuracy of three directions can still be kept in centimeter level, which means that the proposed method can tolerate the loss of 2 consecutive GNSS observation data packets. Typically, the communication delay is typically on the order of milliseconds, and may sometimes reach 100 ms. Therefore, the method of the invention can effectively overcome the influence caused by time delay and has the capability of tolerating the loss of the data packet.
In summary, although the present invention has been described with reference to the preferred embodiments, it should be understood that various changes and modifications can be made by those skilled in the art without departing from the spirit and scope of the invention.

Claims (10)

1. The satellite and inertia combined dynamic-to-dynamic real-time precise relative positioning method is characterized by comprising the following steps of:
the method comprises the steps that firstly, a mobile station receiver and a dynamic reference receiver synchronously sample to obtain original GNSS observation information, carrier phase cycle slip and pseudo range gross error in the original GNSS observation information obtained by respective sampling are respectively detected and eliminated, and carrier phase and pseudo range observation information of the mobile station receiver and the dynamic reference receiver after preprocessing are obtained;
second, determining the carrier phase integer ambiguity
Figure 816296DEST_PATH_IMAGE001
Relative positioning is carried out by utilizing a carrier phase real-time dynamic differential relative positioning method to obtain a relative positioning vector between the mobile station and the moving referencedX RB,DGNSS ；
Third, calculating the measurement update time interval deltat TC Mobile station position increment in
Figure 264595DEST_PATH_IMAGE002
Moving reference position increment
Figure 575490DEST_PATH_IMAGE003
At the same time, the forecast time interval Delta of the mobile station is calculatedt p Increment of position within
Figure 236279DEST_PATH_IMAGE004
；
Fourthly, obtaining a dynamic reference forecasting time interval delta based on a polynomial forecasting method of a sliding windowt p Increment of position within
Figure 69892DEST_PATH_IMAGE005
；
The fifth step is based ondX RB,DGNSS 、、
Figure 221705DEST_PATH_IMAGE007
、
Figure 686184DEST_PATH_IMAGE008
And
Figure 620642DEST_PATH_IMAGE009
and synthesizing and outputting the relative positioning result.
2. The satellite and inertial combined dynamic-dynamic real-time precise relative positioning method according to claim 1, characterized in that: in the first step, the fault information includes carrier phase cycle slip and pseudorange gross error.
3. The satellite and inertial combined dynamic-dynamic real-time precise relative positioning method according to claim 2, characterized in that: in the first step, for carrier phase cycle slip, a GF method is adopted to be combined with inertial navigation auxiliary cycle slip detection to simultaneously detect and eliminate single-frequency or double-frequency carrier phase cycle slip of different satellites; and for pseudo range gross errors, rejecting the pseudo range gross errors after detection by adopting a single-station and double-station inertial navigation auxiliary gross error detection method.
4. A satellite and inertial combined dynamic real-time precise relative positioning method according to claim 1, 2 or 3, characterized in that: in the second step, the carrier phase integer ambiguity is first determined
Figure 879585DEST_PATH_IMAGE010
If the carrier phase integer ambiguity is known
Figure 915806DEST_PATH_IMAGE010
If known, the relative positioning vector between the mobile station and the moving reference is obtained by directly utilizing the carrier phase real-time dynamic differential relative positioning method to perform relative positioningdX RB,DGNSS (ii) a If carrier phase integer ambiguity
Figure 918397DEST_PATH_IMAGE010
Unknown, based on preprocessed carrier phase and pseudorange observations of the mobile station receiver and the mobile reference receiverSolving the carrier phase integer ambiguity by using least square ambiguity reduction correlation adjustment method, and then obtaining the relative positioning vector between the mobile station and the moving reference by using the relative positioning method of the real-time dynamic difference of the carrier phasedX RB,DGNSS 。
5. The satellite and inertial combined dynamic-dynamic real-time precise relative positioning method according to claim 4, characterized in that: relative positioning vector between mobile station and moving referencedX RB,DGNSS The obtaining method comprises the following steps:
the phase of the double-difference carrier wave between the stations is expressed as follows:
Figure 972940DEST_PATH_IMAGE011
wherein the DD operator expression is
Figure 668364DEST_PATH_IMAGE012
I.e. by
Figure 441148DEST_PATH_IMAGE013
，，
Figure 907212DEST_PATH_IMAGE015
，
Figure 507958DEST_PATH_IMAGE016
；
Figure 768038DEST_PATH_IMAGE017
、
Figure 378011DEST_PATH_IMAGE018
Respectively representing mobile stationsRObserved satelliteiAndjthe carrier phase observation of (a);
Figure 892300DEST_PATH_IMAGE019
、
Figure 929526DEST_PATH_IMAGE020
respectively representing moving referencesBObserved satelliteiAndjthe carrier phase observation of (a);
Figure 411323DEST_PATH_IMAGE021
、
Figure 559408DEST_PATH_IMAGE022
respectively representing mobile stationsRObserved satelliteiAndjthe distance between the star and the ground;
Figure 443050DEST_PATH_IMAGE023
、
Figure 651178DEST_PATH_IMAGE024
respectively representing moving referencesBObserved satelliteiAndjthe distance between the star and the ground;
Figure 368073DEST_PATH_IMAGE025
、
Figure 54270DEST_PATH_IMAGE026
respectively representing mobile stationsRObserved satelliteiAndjcarrier phase ambiguity of (a);
Figure 792418DEST_PATH_IMAGE027
、
Figure 171447DEST_PATH_IMAGE028
respectively representing moving referencesBObserved satelliteiAndjcarrier phase ambiguity of (a);
Figure 627836DEST_PATH_IMAGE029
、
Figure 134035DEST_PATH_IMAGE030
respectively representing mobile stationsRObserved satelliteiAndjthe measurement error of (2);
Figure 461111DEST_PATH_IMAGE031
、
Figure 11041DEST_PATH_IMAGE032
respectively representing moving referencesBObserved satelliteiAndjthe error in the measurement of (a) is,
Figure 954727DEST_PATH_IMAGE033
represents a signal wavelength;
ignoring the inter-station sight vector change, the observation equation after linearization is expressed as:
Figure 248305DEST_PATH_IMAGE034
wherein
Figure 961046DEST_PATH_IMAGE035
And
Figure 167030DEST_PATH_IMAGE036
respectively representing mobile stationsRObserved satellite
Figure 598012DEST_PATH_IMAGE037
And
Figure 429701DEST_PATH_IMAGE038
the line-of-sight vector of (a),b RB representing the baseline vector to be solvedQuantity, i.e. the relative positioning vector between the mobile station and the moving reference;
given at least 4 co-visitors, the relative position vector between the mobile station and the moving referenceb RB Obtained by a least squares method, and obtaining a relative position vector between the mobile station and the moving referenceb RB Is marked asdX RB,DGNSS 。
6. The satellite and inertial combined dynamic-dynamic real-time precise relative positioning method according to claim 5, characterized in that: the carrier phase ambiguity is calculated and obtained by a least square ambiguity reduction correlation adjustment method based on double-difference carrier phases and pseudo ranges.
7. The satellite and inertial combined dynamic-dynamic real-time precise relative positioning method according to claim 5, characterized in that: the third step is realized by the following steps:
determining a state vector, wherein the state vector comprises 9 navigation error states and 6 sensor deviation states, and is expressed as:
Figure 262528DEST_PATH_IMAGE039
wherein
Figure 888682DEST_PATH_IMAGE040
、
Figure 557691DEST_PATH_IMAGE041
、
Figure 193072DEST_PATH_IMAGE042
、
Figure 349247DEST_PATH_IMAGE043
And k respectively representkThe position, the speed, the attitude error vector, the adding table and the gyro zero offset error of the epoch;
constructing an observation equation by using GNSS original observation information after the first step of fault detection and elimination, and then obtaining the satelliteiIn thatkThe observation information vector after epoch linearization is represented as:
Figure 677460DEST_PATH_IMAGE044
whereinrIndicating the reference satellite number, indicating the time difference operator,
Figure 551875DEST_PATH_IMAGE045
after representing linearizationkSatellite to be detected of epochiAnd a reference satelliterThe single-differenced pseudoranges between the satellites,
Figure 256526DEST_PATH_IMAGE046
after linearizationkA time difference term of an inter-satellite single difference carrier phase of an epoch;
suppose the total number of satellites to be detected ism，kThe observation matrix of the epoch is:
Figure 376247DEST_PATH_IMAGE047
the time update and measurement update equations are expressed as:
Figure 609782DEST_PATH_IMAGE048
wherein
Figure 237072DEST_PATH_IMAGE049
To representk-1 epoch tokThe state transition matrix of the epoch is,W k-1to representk-1 epoch tokThe course of the epoch noise matrix is,H k to representkA design matrix of the epoch is used,V k to representkThe measured noise matrix of the epoch is,X k-1to representk-a state vector of 1 epoch;
suppose thatPExpressing the variance matrix, performing Kalman filtering, to obtainkTo k+△t TC Epoch, measurement update time interval is Δt TC State vector delta of
Figure 214256DEST_PATH_IMAGE050
And its corresponding variance matrix
Figure 610602DEST_PATH_IMAGE051
The expression of (a) is as follows:
Figure 765771DEST_PATH_IMAGE052
Figure 880357DEST_PATH_IMAGE053
wherein:
Figure 395652DEST_PATH_IMAGE054
Figure 646505DEST_PATH_IMAGE055
Figure 487422DEST_PATH_IMAGE056
Figure 823726DEST_PATH_IMAGE057
Figure 159023DEST_PATH_IMAGE058
Figure 529962DEST_PATH_IMAGE059
to representkEpoch tok+△t TC The state transition matrix of the epoch is,Ithe unit matrix is represented by a matrix of units,
Figure 276201DEST_PATH_IMAGE060
to representk+△t TC The matrix of the variance of the epoch is,Q k-1to representk-1 epoch tokThe process noise covariance matrix of the epoch,R k to representkA measured noise covariance matrix of epochs;
by using
Figure 834221DEST_PATH_IMAGE061
And
Figure 222477DEST_PATH_IMAGE062
is calculated at the mobile station and the moving reference respectively to obtain the mobile stationkTok+△t TC Epoch state vector delta
Figure 933075DEST_PATH_IMAGE063
Dynamic referencekTok+△t TC Epoch state vector delta
Figure 850216DEST_PATH_IMAGE064
And mobile stationk+△t TC Tok+△t TC +△t p Epoch state position vector delta
Figure 426690DEST_PATH_IMAGE065
Respectively extract
Figure 87479DEST_PATH_IMAGE066
、
Figure 167430DEST_PATH_IMAGE067
And
Figure 521051DEST_PATH_IMAGE068
the first three-dimensionally formed vector of (a) is the measurement update time interval Δ to be obtainedt TC Mobile stationkTok+△t TC Epoch location increment
Figure 67046DEST_PATH_IMAGE069
Moving referencekTok+△t TC Epoch location increment
Figure 531525DEST_PATH_IMAGE070
And the mobile station in the forecast time interval deltat p In the interior of said container body,k+△t TC tok+△t TC +△t p Epoch location increment
Figure 465983DEST_PATH_IMAGE071
。
8. The satellite and inertial combined dynamic-dynamic real-time precise relative positioning method according to claim 7, characterized in that: in the fourth step, the polynomial model in the sliding window-based polynomial forecasting method is defined as follows:
Figure 990505DEST_PATH_IMAGE072
wherein
Figure 10414DEST_PATH_IMAGE073
The variables to be modeled are represented by a graph,T i indicating the relative time with respect to the start of the sliding window,oto representThe order of the film is,
Figure 29317DEST_PATH_IMAGE074
the coefficient of the polynomial is represented by,
Figure 818281DEST_PATH_IMAGE075
andT i the expression is:
Figure 779284DEST_PATH_IMAGE076
whereint real1,The true tail time representing the first position increment of the sliding window,t i real,indicating a sliding windowiThe true tail time of a position increment,
Figure 286489DEST_PATH_IMAGE077
to representt i real,The reference station corresponding to the moment is in the measurement update time interval deltat TC Increment of position within.
9. The satellite and inertial combined dynamic-dynamic real-time precise relative positioning method according to claim 8, characterized in that: in the fourth step, when the number of elements of the three-dimensional sliding window of the position vector is more than that of the elementso+1, all polynomial coefficients are solved using the least squares method, assumingk+△t TC The relative time of the epoch with respect to the start of the sliding window ist c,real The moving reference is then at the forecast time interval Δt p Corresponding relative time ist c,real +△t p At this time, the dynamic reference position is increased
Figure 92771DEST_PATH_IMAGE078
The three-dimensional positions are all calculated by the following formula:
Figure 1821DEST_PATH_IMAGE079
in the formula
Figure 353299DEST_PATH_IMAGE080
、
Figure 613379DEST_PATH_IMAGE081
...
Figure 957772DEST_PATH_IMAGE082
、
Figure 721329DEST_PATH_IMAGE083
、
Figure 758555DEST_PATH_IMAGE084
...
Figure 240352DEST_PATH_IMAGE085
And
Figure 404748DEST_PATH_IMAGE086
、
Figure 22812DEST_PATH_IMAGE087
...
Figure 230939DEST_PATH_IMAGE088
respectively represent
Figure 200032DEST_PATH_IMAGE089
Is/are as followsx、 yAndzthe polynomial coefficients of the three-dimensional directions,
Figure 886228DEST_PATH_IMAGE090
、
Figure 643619DEST_PATH_IMAGE091
and
Figure 22647DEST_PATH_IMAGE092
respectively represent
Figure 213457DEST_PATH_IMAGE093
Is/are as followsx、 yAndzthree-dimensional directions at time intervals Δt p The internal dynamic reference position increment forecast value.
10. The satellite and inertial combined dynamic-dynamic real-time precise relative positioning method according to claim 1, characterized in that: in the fifth step, the relative positioning results are synthesizeddX RB,DGNSS/INS Comprises the following steps:
Figure 968924DEST_PATH_IMAGE094
。
CN202011327522.4A 2020-11-24 2020-11-24 Satellite and inertia combined dynamic-alignment real-time precise relative positioning method Active CN112147663B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011327522.4A CN112147663B (en) 2020-11-24 2020-11-24 Satellite and inertia combined dynamic-alignment real-time precise relative positioning method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011327522.4A CN112147663B (en) 2020-11-24 2020-11-24 Satellite and inertia combined dynamic-alignment real-time precise relative positioning method

Publications (2)

Publication Number Publication Date
CN112147663A true CN112147663A (en) 2020-12-29
CN112147663B CN112147663B (en) 2021-02-09

Family

ID=73887426

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011327522.4A Active CN112147663B (en) 2020-11-24 2020-11-24 Satellite and inertia combined dynamic-alignment real-time precise relative positioning method

Country Status (1)

Country Link
CN (1) CN112147663B (en)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113009537A (en) * 2021-02-18 2021-06-22 中国人民解放军国防科技大学 Inertial navigation auxiliary navigation relative positioning unit epoch partial ambiguity solving method
CN113203414A (en) * 2021-05-21 2021-08-03 北京交通大学 Train positioning method based on GPS + BDS PPP/IMU tight combination
CN113406670A (en) * 2021-06-17 2021-09-17 哈尔滨工程大学 Beidou non-motorized broadcast ephemeris fault real-time monitoring method
CN113703024A (en) * 2021-07-16 2021-11-26 交通运输部水运科学研究所 Beidou-based small ship dynamic management system and management method
CN115407371A (en) * 2022-09-02 2022-11-29 中国人民解放军国防科技大学 PPP-B2B-based real-time high-precision time transfer method and device
CN117168499A (en) * 2023-09-04 2023-12-05 武汉大学 High-frequency dynamic target reference position estimation method and computer readable medium

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104536026A (en) * 2015-01-08 2015-04-22 中国航空无线电电子研究所 Dynamic-to-dynamic real-time measurement system
CN107894232A (en) * 2017-09-29 2018-04-10 湖南航天机电设备与特种材料研究所 A kind of accurate method for locating speed measurement of GNSS/SINS integrated navigations and system
CN108873034A (en) * 2018-03-30 2018-11-23 广州海格通信集团股份有限公司 A kind of implementation method of inertial navigation subcarrier ambiguity resolution
CN109917436A (en) * 2019-04-28 2019-06-21 中国人民解放军国防科技大学 Satellite/inertia combined real-time precise relative motion datum positioning method
US20190324154A1 (en) * 2018-04-24 2019-10-24 Novatel Inc. Global navigation satellite system (gnss) antenna data link
WO2019216261A1 (en) * 2018-05-10 2019-11-14 ソニー株式会社 Information processing device, information processing method and program
CN111578935A (en) * 2020-05-08 2020-08-25 北京航空航天大学 Method for assisting GNSS ambiguity fixing by inertial navigation position increment
CN111965685A (en) * 2020-07-07 2020-11-20 北京自动化控制设备研究所 Low-orbit satellite/inertia combined navigation positioning method based on Doppler information

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104536026A (en) * 2015-01-08 2015-04-22 中国航空无线电电子研究所 Dynamic-to-dynamic real-time measurement system
CN107894232A (en) * 2017-09-29 2018-04-10 湖南航天机电设备与特种材料研究所 A kind of accurate method for locating speed measurement of GNSS/SINS integrated navigations and system
CN108873034A (en) * 2018-03-30 2018-11-23 广州海格通信集团股份有限公司 A kind of implementation method of inertial navigation subcarrier ambiguity resolution
US20190324154A1 (en) * 2018-04-24 2019-10-24 Novatel Inc. Global navigation satellite system (gnss) antenna data link
WO2019216261A1 (en) * 2018-05-10 2019-11-14 ソニー株式会社 Information processing device, information processing method and program
CN109917436A (en) * 2019-04-28 2019-06-21 中国人民解放军国防科技大学 Satellite/inertia combined real-time precise relative motion datum positioning method
CN111578935A (en) * 2020-05-08 2020-08-25 北京航空航天大学 Method for assisting GNSS ambiguity fixing by inertial navigation position increment
CN111965685A (en) * 2020-07-07 2020-11-20 北京自动化控制设备研究所 Low-orbit satellite/inertia combined navigation positioning method based on Doppler information

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
YI DONG 等: "Tightly Coupled GNSS/INS Integration with Robust Sequential Kalman Filter for Accurate Vehicular Navigation", 《SENSORS》 *
王凌轩 等: "INS辅助的实时动态GNSS单频周跳探测", 《武汉大学学报》 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113009537A (en) * 2021-02-18 2021-06-22 中国人民解放军国防科技大学 Inertial navigation auxiliary navigation relative positioning unit epoch partial ambiguity solving method
CN113009537B (en) * 2021-02-18 2023-10-31 中国人民解放军国防科技大学 Inertial navigation assisted defensive navigation relative positioning single epoch part ambiguity solving method
CN113203414A (en) * 2021-05-21 2021-08-03 北京交通大学 Train positioning method based on GPS + BDS PPP/IMU tight combination
CN113406670A (en) * 2021-06-17 2021-09-17 哈尔滨工程大学 Beidou non-motorized broadcast ephemeris fault real-time monitoring method
CN113703024A (en) * 2021-07-16 2021-11-26 交通运输部水运科学研究所 Beidou-based small ship dynamic management system and management method
CN115407371A (en) * 2022-09-02 2022-11-29 中国人民解放军国防科技大学 PPP-B2B-based real-time high-precision time transfer method and device
CN115407371B (en) * 2022-09-02 2023-08-15 中国人民解放军国防科技大学 PPP-B2B-based real-time high-precision time transmission method and device
CN117168499A (en) * 2023-09-04 2023-12-05 武汉大学 High-frequency dynamic target reference position estimation method and computer readable medium
CN117168499B (en) * 2023-09-04 2024-04-05 武汉大学 High-frequency dynamic target reference position estimation method and computer readable medium

Also Published As

Publication number Publication date
CN112147663B (en) 2021-02-09

Similar Documents

Publication Publication Date Title
CN112147663B (en) Satellite and inertia combined dynamic-alignment real-time precise relative positioning method
EP3805803A1 (en) Precise point position and real-time kinematic (ppp-rtk) positioning method and device
Xiaohong et al. Instantaneous re-initialization in real-time kinematic PPP with cycle slip fixing
AU2008260578B2 (en) Distance dependant error mitigation in real-time kinematic (RTK) positioning
CN111983654B (en) Method for constructing ionosphere phase scintillation factor in arctic region based on GNSS
US9494693B2 (en) Method, apparatus, and system for determining a position of an object having a global navigation satellite system receiver by processing undifferenced data like carrier-phase measurements and external products like ionosphere data
CN107193028B (en) Kalman relative positioning method based on GNSS
CN105758401A (en) Integrated navigation method and equipment based on multisource information fusion
CN111273687A (en) Multi-unmanned aerial vehicle collaborative relative navigation method based on GNSS observed quantity and inter-aircraft distance measurement
AU2009330687A1 (en) Navigation receiver and method for combined use of a standard RTK system and a global carrier-phase differential positioning system
CN111965673A (en) Time frequency transfer method of single-frequency precise single-point positioning algorithm based on multiple GNSS
CN111427068B (en) Method for monitoring integrity of ephemeris faults of type A satellites of dynamic-to-dynamic platform local augmentation
CN115267855B (en) Abnormal value detection method and differential positioning method in GNSS-INS tight combination
CN112526569B (en) Multi-epoch step-by-step ambiguity solving method for inertial navigation auxiliary navigation relative positioning
Li et al. Predicting atmospheric delays for rapid ambiguity resolution in precise point positioning
EP2092362B1 (en) Phase based measurement corrections
Dong et al. Low-latency, high-rate, high-precision relative positioning with moving base in real time
Yoder et al. Low-Cost Inertial Aiding for Deep-Urban Tightly Coupled Multi-Antenna Precise GNSS
US6704650B1 (en) Technique for accurate distance and velocity calculations using the global positioning system (GPS)
Chen et al. A new cycle slip detection and repair method for single-frequency GNSS data
CN112630811A (en) Real-time PPP-RTK combined positioning method
Gunning Safety critical bounds for precise positioning for aviation and autonomy
CN115267858A (en) Precise single-point positioning method assisted by regional navigation system
CN112113567B (en) Airborne comprehensive landing navigation method
KR102263393B1 (en) Method for generating corrections for rtk system using carrier phase observations of multiple receivers, and satellite navigation augmentation system using the method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant