CN115327588A - Network RTK-based high-precision positioning method for unmanned automatic operation special vehicle - Google Patents

Network RTK-based high-precision positioning method for unmanned automatic operation special vehicle Download PDF

Info

Publication number
CN115327588A
CN115327588A CN202210820920.2A CN202210820920A CN115327588A CN 115327588 A CN115327588 A CN 115327588A CN 202210820920 A CN202210820920 A CN 202210820920A CN 115327588 A CN115327588 A CN 115327588A
Authority
CN
China
Prior art keywords
formula
ambiguity
difference
unmanned
special
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210820920.2A
Other languages
Chinese (zh)
Inventor
赵忠义
崔永志
王慎玉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Liaoning University of Technology
Original Assignee
Liaoning University of 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 Liaoning University of Technology filed Critical Liaoning University of Technology
Priority to CN202210820920.2A priority Critical patent/CN115327588A/en
Publication of CN115327588A publication Critical patent/CN115327588A/en
Pending legal-status Critical Current

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/40—Correcting position, velocity or attitude
    • G01S19/41—Differential correction, e.g. DGPS [differential GPS]
    • 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/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain
    • 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

Abstract

The invention discloses a network RTK-based high-precision positioning method for a special unmanned automatic operation vehicle, which comprises the following steps of: step 1: establishing a GNSS differential observation equation, and eliminating errors such as satellite clock error, atmospheric delay error and receiver clock error; step 2: estimating the integer ambiguity; and 3, step 3: performing decorrelation and search on the integer ambiguity; and 4, step 4: the ambiguity is checked; and 5: detecting and repairing cycle slip; and 6: the attitude information of the unmanned automatic operation special vehicle is sensed, meanwhile, the real-time motion measurement information is converted into a navigation coordinate system, and the position and the speed are updated. The invention can improve the positioning precision and reliability and greatly meet the requirements of the unmanned automatic operation special vehicle in practical application.

Description

Network RTK-based high-precision positioning method for unmanned automatic operation special vehicle
Technical Field
The invention relates to the technical field of positioning, in particular to a network RTK-based high-precision positioning method for a special vehicle for unmanned automatic operation.
Background
RTK, real-time-kinematic differential positioning technology (Real-time kinematic), is a technology for correcting a positioning result by using a GPS carrier phase observation value, so as to perform Real-time kinematic relative positioning. The basic operating principle of RTK: and a data communication link is established between the reference station and the nearby mobile user, the mobile station receives the carrier phase observation value of the reference station and the position information of the reference station through the communication link, then the length of the base station is solved by using a static relative measurement processing method, and finally the position information of the point to be measured is obtained. With the continuous development of the RTK technology, the RTK technology has been widely applied to the fields of transportation, engineering surveying and mapping, urban planning, agricultural machinery, geological disaster monitoring and the like, and gradually replaces the traditional navigation positioning method.
The modern society has higher and higher requirements on the positioning accuracy of navigation, and the single navigation technology has low positioning accuracy due to the defects of the single navigation technology and is difficult to adapt to complex working environments, so that the requirement of modern high-accuracy positioning cannot be met. Currently, the commonly used positioning methods mainly include: satellite positioning, inertial navigation positioning, visual mileage algorithm positioning, ultra wide band wireless positioning and the like.
The GPS system in the united states was the earliest and the most developed, and is commonly referred to as global navigation satellite system (gnss) together with the russian GLONASS navigation satellite system (GLONASS), the beidou navigation satellite system (BDS) in our country and the european union GALILEO system. Currently, the location services of GNSS have penetrated aspects of human social life, such as autopilot, transportation, engineering surveying and mapping, agricultural machinery, and geological disaster monitoring. Although the widely-used GNSS positioning system has the advantages of wide coverage area, short observation time, and capability of providing stable and high-precision positioning under the condition of good number of visible satellites, once the GNSS positioning system enters forest canyons, high-rise forests and indoor environments with complex and changeable structures, the GNSS positioning precision is finally low due to insufficient number of visible satellites, and even the positioning fails.
An Inertial Navigation System (INS) is composed of an Inertial Measurement Unit (IMU) and a computer, wherein the Inertial Measurement Unit integrates basic elements such as a gyroscope, an accelerometer, a geomagnetic sensor and the like, and is the core of a modern Inertial Navigation System. The inertial navigation measures the acceleration of a carrier by applying Newton's mechanical principle, and then obtains the position and the speed of the carrier by applying a differential calculation method in mathematics. The advantage of inertial navigation is that it can be protected from the external environment and does not rely on external navigation stations or radios during transmission. However, since the inertial element itself has a certain manufacturing error, the positioning deviation of the inertial navigation device will gradually accumulate over time. At present, due to the progress of the production process, the accuracy of the gyroscope is also increasing, and therefore, the platform type inertial navigation system is gradually replaced by the strapdown type inertial navigation system.
The visual mileage algorithm positioning is to use machine vision to obtain images of the surrounding environment, process the collected images through an algorithm, extract characteristic values, and finally form a visual mileage meter through a continuous iteration mode by adopting accumulation summation to obtain the position information of the carrier. Because the camera is installed on the vehicle body, when the unmanned vehicle moves at a relatively fast speed, the image shot by the camera becomes blurred, which causes the difficulty of the algorithm for extracting the characteristic value to be increased, even the algorithm to be invalid.
Ultra Wide Band (UWB) wireless positioning is to transmit data by sending and receiving extremely narrow pulses having nanosecond level, even lower than nanosecond level, thereby having bandwidth of GHz level, and mainly positioning by using time difference of arrival at a base station. If three base stations with known coordinate points receive signals and the distances between the transmitted pulse signals and the three base stations are different, the time points of the signals received by the three base stations are different. The ultra-wideband system has the advantages of strong penetrating power, good multipath resistance effect, high safety and the like, and is mainly applied to indoor positioning.
The unmanned special operation vehicle is mainly used for an unmanned operation platform for judging images, sound, air quality, target temperature and the like around an observation area, and if the unmanned special operation vehicle only depends on one navigation positioning, the requirement of the unmanned special operation vehicle on high precision is difficult to meet, and each navigation mode has certain defects. For the inertial navigation and satellite navigation which are most widely applied, the inertial navigation is not interfered by the external environment, but the error is accumulated continuously along with the time, the satellite navigation is greatly influenced by the environment, and once the inertial navigation and the satellite navigation enter a shielding environment, the satellite signal is easy to lose lock, so that the positioning loss is caused. Therefore, in practical engineering applications, especially for unmanned special work vehicles, a high-precision and high-reliability positioning method is needed.
Disclosure of Invention
The invention aims to provide a network RTK-based high-precision positioning method for an unmanned automatic operation special vehicle, which is used for solving the problem that the existing positioning method cannot meet the requirements of the unmanned automatic operation special vehicle on positioning precision and reliability in practical application.
In order to realize the purpose, the technical scheme of the invention is as follows:
a high-precision positioning method for a special unmanned automatic operation vehicle based on network RTK comprises the following steps:
step 1: establishing a GNSS differential observation equation, and eliminating errors such as satellite clock error, atmospheric delay error and receiver clock error;
step 2: estimating the integer ambiguity;
and step 3: performing decorrelation and search on the integer ambiguity;
and 4, step 4: checking the ambiguity;
and 5: detecting and repairing cycle slip;
step 6: the attitude information of the unmanned automatic operation special vehicle is sensed, meanwhile, the real-time motion measurement information is converted into a navigation coordinate system, and the position and the speed are updated.
Further, as a preferred technical scheme, the specific process of the step 1 is as follows:
step 1-1: the unified space-time reference comprises the unification of a time system and the unification of a coordinate system, and the unified relation is as follows:
GPST＝UTC+1s×n-19s (1)
BDT＝GPST-14s (2)
in the formula, when the GPST is a GPS, the GPS represents a special time system controlled by an atomic clock of a GPS main control station; UTC represents coordinated universal time; n represents an adjustment parameter, and a specific numerical value is issued by the international terrestrial rotation service bureau; BDT represents the time of big Dipper;
step 1-2: establishing a pseudo range and carrier phase observation equation:
Figure BDA0003742385530000031
Figure BDA0003742385530000032
in the formula, the superscript s represents the observation satellite, and the subscript r represents the receiver; p and
Figure BDA0003742385530000033
respectively are pseudo range and carrier phase observed quantity; rho is the real distance between the satellite signal transmitting time and the receiver; c is the speed of light in vacuum; dt is r And dt s Respectively a receiver clock error and a satellite clock error; t is tropospheric error of the signal in the process of propagation in the atmosphere; i is an ionospheric error of a signal in the process of propagation in the atmosphere; λ is carrier phase wavelength; n is the integer ambiguity in weeks; ε p and
Figure BDA0003742385530000034
other errors in pseudorange and carrier phase (including noise errors), respectively;
step 1-3: establishing a single difference observation model, carrying out inter-station difference between a reference station r and a mobile station q, and carrying out difference calculation on the same satellite s by the reference station receiver r and the mobile station receiver q, wherein the pseudo range and carrier phase observation equation of the reference station r is shown as (3) and (4), and the pseudo range and carrier phase observation equation of the mobile station q is as follows:
Figure BDA0003742385530000035
Figure BDA0003742385530000036
the difference between the formula (3) and the formula (5), and the difference between the formula (4) and the formula (6) can be the satellite clock difference dt q Thus, a single difference model is obtained:
Figure BDA0003742385530000041
Figure BDA0003742385530000042
in equations (7) and (8), Δ represents a single difference, and for the case of a short baseline, the difference between the ionospheric delay error and the tropospheric delay error between the reference station and the mobile station is small, so that the difference between the tropospheric delay and the ionospheric delay in equations (7) and (8) can be approximated to be 0, and it can be obtained:
Figure BDA0003742385530000043
Figure BDA0003742385530000044
step 1-4: establishing a double-difference observation model, and setting a reference satellite y, then the equation of the single difference of the reference station r and the mobile station q is similar to the equations (9) and (10),
Figure BDA0003742385530000045
Figure BDA0003742385530000046
by subtracting equations (9) and (11), (10) and (12), the receiver clock difference term is eliminated to obtain the double-difference observation equation:
Figure BDA0003742385530000047
Figure BDA0003742385530000048
in the formulae (13) and (14),
Figure BDA0003742385530000049
representing double difference symbols, wherein the integer property of double difference integer ambiguity is kept, and the double difference integer ambiguity can be effectively separated;
if N available satellites are provided, the number of observation equations for carrier phase and pseudorange is 2 (N-1), where the unknowns include 3 directional coordinates for baseline double-differenced pseudorange, and N-1 double-differenced ambiguities, so according to the principle of mathematically solving the equations: the equation number is greater than or equal to the unknown number. To solve for the floating solution of ambiguities, then: 2 (N-1) is more than or equal to 3+ (N-1), and N =4 is obtained by solving, so that 4 visible satellites are observed at the same time, and the minimum requirement for realizing differential positioning by a single system is met;
when 4 satellites a, b, c and d are used for resolving, pseudo-range double-difference observables are added, a 6-dimensional equation set can be established, and a satellite a with a large altitude angle is used as a main satellite, so that the following results are obtained:
Figure BDA0003742385530000051
where [ ln ] is the unit vector from the satellite to the receiver, [ δ X δ Y δ Z ] is the receiver coarse coordinate correction, [ Δ R is the double difference pseudorange observation, and equation (15) is written as a matrix:
Y＝HX (16)
in the formula (16), H is a design matrix formed by unit vectors from the receiver to the satellite, X is a base line vector to be solved and a ambiguity solution, and a float solution of ambiguity can be solved by using a least square method:
Figure BDA0003742385530000052
in the formula (17), the compound represented by the formula (I),
Figure BDA0003742385530000053
for floating-point solutions of ambiguities, let the baseline vector
Figure BDA0003742385530000054
Double difference integer ambiguity vector
Figure BDA0003742385530000055
The least squares result can be expressed as:
Figure BDA0003742385530000056
Figure BDA0003742385530000057
in equation (19), Q represents the covariance matrix of vector X,
Figure BDA0003742385530000058
a covariance matrix representing the baseline vector,
Figure BDA0003742385530000059
a covariance matrix representing the ambiguity vector,
Figure BDA00037423855300000510
and
Figure BDA00037423855300000511
a cross-covariance matrix representing the baseline vector and the ambiguity vector.
Further, as a preferred technical scheme, the specific process of step 2 is as follows:
step 2-1: solving floating point
Figure BDA00037423855300000512
Decomposing into the form of a baseline vector and an integer ambiguity vector to obtain:
Figure BDA00037423855300000513
in the equation (20), y represents a pseudo-range and carrier phase double-difference observed quantity,
Figure BDA00037423855300000514
is a baseline vector, M is a coefficient matrix of the baseline vector, N is an integer ambiguity coefficient matrix,
Figure BDA00037423855300000515
is a double difference integer ambiguity, and epsilon is an unmodeled error quantity;
step 2-2: converting the solved ambiguity problem into a least square problem:
Figure BDA00037423855300000516
in the formula (21), Q y Is the covariance matrix of y;
step 2-3: obtaining integer ambiguity float solution by equation (18) and equation (19)
Figure BDA0003742385530000061
Sum covariance matrix Q, then float solution with ambiguity vector
Figure BDA0003742385530000062
Sum covariance matrix
Figure BDA0003742385530000063
Integer solution for estimating ambiguities
Figure BDA0003742385530000064
Minimizing the sum of squared ambiguity residuals yields:
Figure BDA0003742385530000065
further, as a preferred technical solution, the specific process of step 3 is:
let transform matrix z satisfy:
Figure BDA0003742385530000066
then the corresponding covariance matrix is:
Figure BDA0003742385530000067
in the formula (24), Z is an integer matrix whose modulo is 1, all elements of the matrix are integers, and the modulo value is 1, | det (Z) | =1, and the correlation between the ambiguity and the target function Ω can be minimized by selecting an appropriate Z matrix and transforming;
Figure BDA0003742385530000068
obtaining an optimal solution of an integer matrix satisfying the formula (25), and setting a search space χ 2 Namely:
Figure BDA0003742385530000069
the search space is a multi-dimensional ellipsoid, shape and covariance matrix
Figure BDA00037423855300000610
In this regard, the magnitude is determined by the constant χ.
Further, as a preferred technical solution, the specific process of step 4 is:
setting the minimum value of the searched residual error as omega min Then the residual sub-minimum is Ω sec Then, it is possible to obtain:
Figure BDA00037423855300000611
in the formula (27), k is an empirical constant and cannot change with changes in model strength, measurement environment and the like, generally k =3 is used as a judgment condition for inspection, and ratio is greater than or equal to 3, and the optimal solution obtained by searching is considered to be a correct solution;
obtaining integer ambiguity fixing solution through ambiguity searching
Figure BDA00037423855300000612
By inverse transformation, a fixed solution of the original double-difference ambiguities can be obtained
Figure BDA00037423855300000613
Figure BDA00037423855300000614
Double-difference integer ambiguity vector obtained by integer ambiguity resolution
Figure BDA0003742385530000071
Further improving the accuracy of baseline vector estimation by integer property of ambiguity, obtained by LAMBDA search
Figure BDA0003742385530000072
It is possible to obtain:
Figure BDA0003742385530000073
Figure BDA0003742385530000074
in the above formula, the first and second carbon atoms are,
Figure BDA0003742385530000075
is a fixed solution to the baseline vector modifier,
Figure BDA0003742385530000076
is composed of
Figure BDA0003742385530000077
The covariance matrix of (2).
According to the above results, can obtain
Figure BDA0003742385530000078
Precise three-dimensional coordinates (X) of a rover q ,Y q ,Z q ) T Comprises the following steps:
(X q ,Y q ,Z q ) T ＝(X r0 ,Y r0 ,Z r0 ) T +(δX r ,δY r ,δZ r ) T (31)。
further, as a preferred technical solution, the specific process of step 6 is:
step 6-1: establishing an SINS coordinate system of the unmanned automatic operation special vehicle, wherein the SINS coordinate system comprises a carrier coordinate system, an inertia coordinate system and a geographical coordinate system;
step 6-2: calculating attitude parameters of the unmanned automatic operation special vehicle by adopting a quaternion method, wherein the attitude parameters comprise a course angle, a pitch angle and a roll angle of the unmanned automatic operation special vehicle, and updating the attitude of the strapdown inertial navigation;
step 6-3: updating the speed of the strapdown inertial navigation;
step 6-4: and updating the position of the strapdown inertial navigation.
Further, as a preferred technical scheme, the specific process of the step 6-2 is as follows:
the quaternion is a four-dimensional array, is composed of four elements, and comprises three virtual units and a real unit, and the mathematical expression of the quaternion is as follows:
Figure BDA0003742385530000079
in the formula (32), q 0 、q 1 、q 2 And q is 3 Is any real number, and is a real number,
Figure BDA00037423855300000710
and
Figure BDA00037423855300000711
are mutually perpendicular unit vectors; the matrix of quaternions is:
Figure BDA00037423855300000712
in the formula (33), θ is an angle,
Figure BDA0003742385530000081
the rotation matrix from which the quaternion representation can be derived is:
Figure BDA0003742385530000082
in the process of rotating from n to b, as three coordinate axes in two coordinate systems are mutually vertical in pairs,
Figure BDA0003742385530000083
and
Figure BDA0003742385530000084
two reciprocal coordinate system rotation processes are shown, which are related as follows:
Figure BDA0003742385530000085
in the formula (35), psi is a course angle of the unmanned special operating vehicle; theta is the pitch angle of the unmanned special operating vehicle; gamma is a roll angle of the unmanned special operation vehicle; the combined type (33) and the formula (34) are as follows:
Figure BDA0003742385530000086
γ＝arcsin(-C 13 )＝arcsin(-2(q 1 q 3 -q 0 q 2 )) (37)
Figure BDA0003742385530000087
observing the above formula can obtain the quaternion q 0 、q 1 、q 2 And q is 3 The values of the attitude angle psi, the pitch angle theta and the roll angle gamma of the unmanned special type operation vehicle can be calculated.
Further, as a preferred technical scheme, the specific process of the step 6-3 is
The velocity equation obtained from the specific force equation is:
Figure BDA0003742385530000088
in formula (39), f b Is the accelerometer output;
Figure BDA0003742385530000089
is the Coriolis acceleration;
Figure BDA00037423855300000810
for traction acceleration, speed updating can be completed through integral operation, and a specific expression is as follows:
Figure BDA00037423855300000811
in the formula (40), V m-1 The speed at the previous moment; Δ V smf Is a specific force integral increment; Δ V g/corm Velocity increments due to unwanted acceleration;
Figure BDA0003742385530000091
equations (39) to (41) are velocity updates for inertial navigation.
Further, as a preferred technical scheme, the specific process of the step 6-4 is as follows:
the longitude and latitude and the height of the unmanned vehicle are obtained through integral operation according to the speed of the unmanned automatic operation special vehicle under a navigation coordinate system, and the specific formula is as follows:
Figure BDA0003742385530000092
Figure BDA0003742385530000093
Figure BDA0003742385530000094
wherein L represents latitude, L m Represents the latitude at time m, L m-1 Representing the latitude at the m-1 moment; lambda [ alpha ] m Represents longitude at time m; lambda m-1 Represents longitude at time m-1; h is m Representing the m-time height, h, of the unmanned vehicle m-1 Representing the height of the unmanned vehicle at the m-1 moment; t is t m Represents time m, t m-1 Represents the m-1 time;
Figure BDA0003742385530000095
representing the speed n power of the x direction;
Figure BDA0003742385530000096
represents the speed n power in the y direction;
Figure BDA0003742385530000097
represents the speed n power in the z direction; r is n Represents the radius of the earth; h represents the satellite altitude.
Further, as a preferred technical solution, the method further includes an initial coarse alignment step, specifically: because the zero offset of the gyroscope is far greater than the rotation angular velocity of the earth, the gyroscope can hardly sense the rotation angular velocity of the earth, and therefore the magnitude of the initial course angle cannot be estimated; since accelerometers have a relatively high accuracy, the separation of gravitational acceleration sensed by accelerometers estimates the initial horizontal attitude angle by the accelerometer: a pitch angle theta and a roll angle gamma;
Figure BDA0003742385530000098
the three elements on the left side of the formula (45) are the output of the accelerometer under the carrier coordinate system; and three elements in brackets on the right side are projections of the gravity acceleration under the navigation coordinate system, under the static condition, the gravity acceleration in the general rough alignment only has a projection of-g in a Z axis, and the projections of other axes are all 0, namely g n ＝[0 0 -g] T (ii) a It is known that
Figure BDA0003742385530000099
Carry-over (2.50) gives:
Figure BDA00037423855300000910
Figure BDA00037423855300000911
Figure BDA0003742385530000101
g is the gravity acceleration, and the heading angle psi is an unknown quantity, which can be solved reversely:
Figure BDA0003742385530000102
Figure BDA0003742385530000103
compared with the prior art, the invention has the following beneficial effects:
(1) The invention utilizes the GPS double antenna and the MEMS inertial unit to carry out combined navigation positioning, eliminates main errors by establishing and deducing single difference and double difference observation equations, and combines the searching and the inspection of the ambiguity of the whole circle and the detection and the repair of the cycle slip, thereby effectively improving the positioning precision and the reliability and greatly meeting the requirements of the unmanned automatic operation special vehicle in practical application.
(2) The initial attitude of the strapdown inertial navigation system with certain accuracy is estimated by adopting initial coarse alignment, so that the initial attitude conversion rectangle is convenient to further determine, the error of the initial point is reduced as much as possible, and the condition that the subsequent accuracy of the navigation system is poorer and poorer is avoided.
Drawings
FIG. 1 is a schematic flow diagram of the process of the present invention;
FIG. 2 is a single difference observation model of the present invention;
FIG. 3 is a double-difference observation model of the present invention;
FIG. 4 is a schematic view of the integer ambiguity search process of the present invention;
FIG. 5 is a schematic diagram of the phase cycle slip of the present invention;
FIG. 6 is a schematic structural diagram of a strapdown inertial navigation system.
Detailed Description
The present invention will be described in further detail with reference to examples and drawings, but the present invention is not limited to these examples.
Examples
As shown in fig. 1, the method for high-precision positioning of a special vehicle for unmanned automated operation based on network RTK in this embodiment includes the following steps:
step 1: establishing a GNSS differential observation equation, and eliminating errors such as satellite clock error, atmospheric delay error and receiver clock error;
step 2: estimating the integer ambiguity;
and step 3: performing decorrelation and search on the integer ambiguity;
and 4, step 4: checking the ambiguity;
and 5: detecting and repairing cycle slip;
step 6: the attitude information of the unmanned automatic operation special vehicle is sensed, meanwhile, the real-time motion measurement information is converted into a navigation coordinate system, and the position and the speed are updated.
The specific method of this example is as follows:
firstly, the space-time reference is unified, including the unification of a time system and a coordinate system:
(1) Unification of time systems
The GPS establishes a special Time system controlled by an atomic clock of a GPS main control station, namely the Time of GPS (GPS Time, GPST), the Time of GPST and Coordinated Universal Time (UTC) is specified, the initial epoch is 00 min 00 s/6.1/6/1986 of UTC, the Time is continuously accumulated later, the difference between the GPST and the UTC is integral multiple of the second, and the relationship between the GPST and the UTC is as follows:
GPST＝UTC+1s×n-19s (1)
BDT＝GPST-14s (2)
in the formula, when the GPST is a GPS, the GPS represents a special time system controlled by an atomic clock of a GPS main control station; UTC represents coordinated universal time; n represents an adjustment parameter, and a specific numerical value is issued by the international terrestrial rotation service bureau; BDT represents the time of big Dipper; the initial epoch of the Beidou (BeiDou Time, BDT) is 00 minutes 00 seconds at 1 month and 1 day 2006, and the international unit System (SI) seconds are adopted as basic units to be continuously accumulated without leap seconds. The BDT is sourced to a coordinated universal Time maintained by the National Time Service Center (NTSC) of the chinese academy of sciences. GPST is 14 seconds faster than BDT because the two start epochs are different. For the unification of the coordinate system, the WGS-84 coordinate system is consistent with the CGCS2000 coordinate system, only the oblateness rate has slight difference among 4 basic constants, the difference of the oblateness rate does not influence the change of precision, and influences the geodetic latitude and the geodetic height to generate the difference of about 0.1mm, however, under the condition of the current positioning precision requirement, the difference can be almost ignored, so the WGS-84 coordinate system and the CGCS2000 coordinate system can be considered to be consistent.
In the GNSS data processing, a functional relation between the original observed quantity of the receiver and various factors influencing signal propagation is a GNSS observation equation, a satellite continuously broadcasts navigation signals in the GNSS positioning, and the receiver performs algorithm calculation on received pseudo-range, carrier phase, doppler and other information to realize the navigation positioning. The pseudo-range and carrier phase observations are two types of original observations, and the observation equation is as follows:
Figure BDA0003742385530000121
Figure BDA0003742385530000122
in the formula, the superscript s represents an observation satellite, and the subscript r represents a receiver; p and
Figure BDA0003742385530000123
respectively are pseudo range and carrier phase observed quantity; rho is the real distance between the satellite signal transmitting time and the receiver; c is the speed of light in vacuum; dt r And dt s Respectively a receiver clock error and a satellite clock error; t is tropospheric error of the signal in the process of propagation in the atmosphere; i is an ionized layer error of a signal in the process of propagation in the atmosphere; λ is carrier phase wavelength; n is the integer ambiguity in weeks; ε p and
Figure BDA0003742385530000124
other errors in pseudorange and carrier phase (including noise errors), respectively;
and establishing a single difference observation model. The reference station r and the mobile station q perform inter-station difference, the reference station receiver r and the mobile station receiver q perform difference calculation on the same satellite s, and a single difference observation model is shown in fig. 2. The pseudorange and carrier phase observation equations of the reference station r are expressed by the following equations (3) and (4), and the pseudorange and carrier phase observation equation of the mobile station q is expressed as follows:
Figure BDA0003742385530000125
Figure BDA0003742385530000126
the difference between the formula (3) and the formula (5), and the difference between the formula (4) and the formula (6) can be the satellite clock difference dt q Thus, a single difference model is obtained:
Figure BDA0003742385530000127
Figure BDA0003742385530000128
in equations (7) and (8), Δ represents a single difference, and for the case of a short baseline, the difference between the ionospheric delay error and the tropospheric delay error between the reference station and the mobile station is small, so that the difference between the tropospheric delay and the ionospheric delay in equations (7) and (8) can be approximated to be 0, and it can be obtained:
Figure BDA0003742385530000129
Figure BDA00037423855300001210
and establishing a double-difference observation model. As shown in fig. 3, double differencing is a method of inter-satellite differencing for another reference satellite, and both the atmospheric delay error and the receiver clock error can be eliminated by a double differencing observation model, and assuming that the reference satellite y, the equation for the single difference between the reference station r and the rover station q is similar to equations (9) and (10),
Figure BDA00037423855300001211
Figure BDA0003742385530000131
by subtracting equations (9) and (11), (10) and (12), the receiver clock difference term is eliminated to obtain the double-difference observation equation:
Figure BDA0003742385530000132
Figure BDA0003742385530000133
in the formulae (13) and (14),
Figure BDA0003742385530000134
representing double difference symbols, wherein the integer property of double difference integer ambiguity is kept, and the double difference integer ambiguity can be effectively separated;
with N available satellites, the number of carrier phase and pseudorange observation equations is 2 (N-1), where the unknowns include the 3 directional coordinates of the baseline double-differenced pseudorange, and also the N-1 double-differenced ambiguities, so according to the principle of mathematically solving the equations: the equation number is greater than or equal to the unknown number. To solve for the floating solution of ambiguities, then: 2 (N-1) is more than or equal to 3+ (N-1), and N =4 is obtained by solving, so that 4 visible satellites are observed at the same time, and the minimum requirement for realizing differential positioning by a single system is met;
when 4 satellites a, b, c and d are used for resolving, pseudo-range double-difference observed quantities are added, a 6-dimensional equation set can be established, and a satellite a with a large altitude angle is used as a main satellite, so that the following results are obtained:
Figure BDA0003742385530000135
wherein, [ l m n [ ]]Is a unit vector from satellite to receiver, [ delta X delta Y delta Z]For the receiver to be a coarse coordinate correction number,
Figure BDA0003742385530000136
for double-differenced pseudorange observations, write equation (15) in matrix form:
Y＝HX (16)
in the formula (16), H is a design matrix formed by unit vectors from a receiver to a satellite, X is a baseline vector to be solved and a ambiguity solution, and a floating solution of the ambiguity can be solved by using a least square method:
Figure BDA0003742385530000137
in the formula (17), the reaction mixture is,
Figure BDA0003742385530000138
for floating-point solutions of ambiguities, let the base-line vector
Figure BDA0003742385530000139
Double difference integer ambiguity vector
Figure BDA00037423855300001310
The least squares result can be expressed as:
Figure BDA0003742385530000141
Figure BDA0003742385530000142
in equation (19), Q represents the covariance matrix of vector X,
Figure BDA0003742385530000143
a covariance matrix representing the baseline vector,
Figure BDA0003742385530000144
a covariance matrix representing the ambiguity vector,
Figure BDA0003742385530000145
and
Figure BDA0003742385530000146
a cross-covariance matrix representing the baseline vector and the ambiguity vector.
The specific process of the integer ambiguity estimation in this embodiment is as follows:
solving floating point
Figure BDA0003742385530000147
Decomposing into the form of a baseline vector and an integer ambiguity vector to obtain:
Figure BDA0003742385530000148
in the equation (20), y represents a pseudo-range and carrier phase double-difference observed quantity,
Figure BDA0003742385530000149
is a baseline vector, M is a coefficient matrix of the baseline vector, N is an integer ambiguity coefficient matrix,
Figure BDA00037423855300001410
is double difference integer ambiguity, epsilon is error quantity not modeled;
converting the solution ambiguity problem into a least square problem:
Figure BDA00037423855300001411
in the formula (21), Q y Is the covariance matrix of y;
if the ambiguities are independent of each other, the covariance matrix of the ambiguities is a diagonal matrix, and the optimal solution can be obtained by a method of rounding nearby. However, in general, the ambiguity covariance matrix is not a diagonal matrix because of strong correlation between ambiguities, which results in a low success rate by rounding. Obtaining integer ambiguity floating solution by formula (18) and formula (19)
Figure BDA00037423855300001412
Sum covariance matrix Q, then float solution with ambiguity vector
Figure BDA00037423855300001413
Sum covariance matrix
Figure BDA00037423855300001414
Integer solution for estimating ambiguities
Figure BDA00037423855300001415
Minimizing the sum of squared ambiguity residuals yields:
Figure BDA00037423855300001416
the specific process of ambiguity reduction correlation and search in this embodiment is as follows:
as shown in fig. 4, the conversion is realized by Z-transform
Figure BDA00037423855300001417
And the matrix is changed into a diagonal matrix, the ambiguity search is simplified, the search time is shortened, and a transformation matrix z is set to satisfy the following conditions:
Figure BDA00037423855300001418
then the corresponding covariance matrix is:
Figure BDA00037423855300001419
in the formula (24), Z is an integer matrix whose modulo is 1, all elements of the matrix are integers, and the modulo value is 1, | det (Z) | =1, and the correlation between the ambiguity and the target function Ω can be minimized by selecting an appropriate Z matrix and transforming;
Figure BDA0003742385530000151
obtaining a whole satisfying the formula (25)Setting the optimal solution of the number matrix, and setting the search space x 2 Namely:
Figure BDA0003742385530000152
the search space is a multi-dimensional ellipsoid, the shape and the covariance matrix
Figure BDA0003742385530000153
The size of the relation is determined by a constant x, and if the value of x is too large, the value in a search space is too much, so that the search efficiency is reduced; conversely, if the value of χ is too small, the search space may not contain the complete ambiguity integer solution, so the value of χ is very important.
The ambiguity checking process of the embodiment includes:
setting the minimum value of the searched residual error to be omega min Then the residual sub-minimum is Ω sec Then, it is possible to obtain:
Figure BDA0003742385530000154
in the formula (27), k is an empirical constant and cannot change with changes in model strength, measurement environment and the like, generally k =3 is used as a judgment condition for inspection, and ratio is greater than or equal to 3, and the optimal solution obtained by searching is considered to be a correct solution;
obtaining integer ambiguity fixing solution through ambiguity searching
Figure BDA0003742385530000155
By inverse transformation, a fixed solution of the original double-difference ambiguity can be obtained
Figure BDA0003742385530000156
Figure BDA0003742385530000157
By integratingDouble-difference integer ambiguity vector obtained by ambiguity solution
Figure BDA0003742385530000158
Further improving the accuracy of baseline vector estimation by integer property of ambiguity, obtained by LAMBDA search
Figure BDA0003742385530000159
It is possible to obtain:
Figure BDA00037423855300001510
Figure BDA00037423855300001511
in the above formula, the first and second carbon atoms are,
Figure BDA00037423855300001512
is a fixed solution to the baseline vector correction,
Figure BDA00037423855300001513
is composed of
Figure BDA00037423855300001514
The covariance matrix of (2).
According to the above results, can obtain
Figure BDA00037423855300001515
Precise three-dimensional coordinates (X) of a rover q ,Y q ,Z q ) T Comprises the following steps:
(X q ,Y q ,Z q ) T ＝(X r0 ,Y r0 ,Z r0 ) T +(δX r ,δY r ,δZ r ) T (31)。
in this embodiment, in the process of performing continuous carrier phase measurement by the receiver, the environment of the receiver often causes signal blocking, and the receiver generates sharp jitter and signalThe reason that the signal-to-noise ratio is too low, the whole-Cycle counting is wrong, and the whole-Cycle jump occurs, that is, the Cycle Slip (Cycle Slip), as shown in fig. 5, the original phase value is a continuous curve function changing with time, when the Cycle Slip occurs, the curve changes significantly, and t is t 1 Observed value of carrier phase corresponding to epoch
Figure BDA0003742385530000161
And
Figure BDA0003742385530000162
subsequent phase observations still include missing cycle slip values.
The detection of cycle slip is also an important issue in GNSS data processing, since cycle slip can have a significant impact on GNSS data quality. Meanwhile, the problems related to cycle slip are difficult to solve due to the influence of factors such as the arbitrariness of cycle slip generation and noise. There are many commonly used methods for detecting and repairing cycle slip, but they can be generally classified into the following three categories:
(1) A method of using the observed value with time change, such as a polynomial fitting method;
(2) Using methods that combine different observations, such as ionospheric residuals;
(3) Methods for estimating residuals using observations, such as kalman filtering and inter-epoch differential observation parameter estimation, are provided.
The strapdown inertial navigation system SINS has the characteristics of small volume, high reliability, simple structure and the like, so the development is very rapid. The method mainly comprises three parts, namely position updating, speed updating and posture updating. Attitude updating is the core of the SINS, which not only senses attitude information of a carrier, but also converts real-time motion measurement information into a navigation coordinate system to realize position and speed updating. The specific calculation process is as follows: through initial alignment, information such as initial posture, position and speed of the carrier can be obtained and used as initial information of navigation calculation; and performing integral operation on the acceleration under the navigation coordinate system to obtain real-time speed and position information, then solving the updated direction cosine matrix or quaternion to finally obtain unmanned vehicle attitude information, wherein a schematic diagram of the strapdown inertial navigation system is shown in fig. 6.
In the navigation system of the special unmanned automatic operation vehicle, because the coordinate system referred by the sensor measurement value and the attitude observation value finally output by the navigation system is different, the coordinate system used by the SINS comprises:
(1) Carrier coordinate system (b): as the name implies, the carrier coordinate system is a coordinate system of the inertial navigation module mounted on the moving carrier, and is consistent with the moving state of the carrier, and is defined as follows:
origin O: the mass center of the unmanned special work vehicle;
an X axis: the device is parallel to the ground and points to the advancing direction of the unmanned special operating vehicle;
z axis: the center of mass of the unmanned special operation vehicle is vertically directed upwards;
y-axis: the right-hand rule was followed and perpendicular to the XOZ plane.
(2) Inertial coordinate system (i): belonging to a spatially-fixed coordinate system, in which the position of a spatial point is represented by (X, Y, Z); it is defined as follows:
origin O: the earth's geocentric;
an X axis: pointing to spring equinox;
z-axis: pointing to the north pole of the earth;
y-axis: following the right-hand rule and perpendicular to the XOZ plane, the plane XOY lies in the equatorial plane.
(3) Geographic coordinate system (g): selecting a northeast coordinate system;
origin O: the gravity center of the unmanned special operating vehicle;
an X axis: pointing in the east direction, namely E;
z-axis: pointing to the sky along the vertical line, namely U;
y-axis: pointing in the north direction, i.e. N.
In this embodiment, the attitude of the strapdown inertial navigation system includes a heading angle, a pitch angle, and a roll angle. The direction information and attitude matrix of the unmanned vehicle taking the northeast coordinate system as the reference system
Figure BDA0003742385530000171
Representing the rotational relationship of the coordinate system of the unmanned vehicle (carrier) to the coordinate system of the northeast sky. At present, three methods of an Euler angle method, a direction cosine matrix and a quaternion method are commonly used as attitude solutions, the quaternion method represents attitude information by 4 variables, compared with the other two methods, the method is very convenient in the aspect of real-time updating of the attitude information, and compared with an algorithm designed based on a Poisson equation, a quaternion algorithm can reduce the calculated amount by 30% on average, so that the quaternion method is selected as a main method for obtaining the attitude by an inertial sensor.
From a mathematical point of view, a quaternion is essentially a four-dimensional array consisting of four elements, including three imaginary units and one real unit. The mathematical expression is as follows:
Figure BDA0003742385530000172
in the formula (32), q 0 、q 1 、q 2 And q is 3 Is any real number, and is a real number,
Figure BDA0003742385530000173
and
Figure BDA0003742385530000174
are mutually perpendicular unit vectors; the matrix of quaternions is:
Figure BDA0003742385530000175
in the formula (33), θ is an angle,
Figure BDA0003742385530000176
the rotation matrix from which the quaternion representation can be derived is:
Figure BDA0003742385530000177
from n toIn the process of b, as three coordinate axes in two coordinate systems are mutually vertical in pairs,
Figure BDA0003742385530000181
and
Figure BDA0003742385530000182
two reciprocal coordinate system rotation processes are shown, which are related as follows:
Figure BDA0003742385530000183
in the formula (35), psi is a course angle of the unmanned special work vehicle; theta is the pitch angle of the unmanned special operation vehicle; gamma is a roll angle of the unmanned special operation vehicle; the joint type (33) and the formula (34) are as follows:
Figure BDA0003742385530000184
γ＝arcsin(-C 13 )＝arcsin(-2(q 1 q 3 -q 0 q 2 )) (37)
Figure BDA0003742385530000185
observing the above formula can obtain the quaternion q 0 、q 1 、q 2 And q is 3 The values of the attitude angle psi, the pitch angle theta and the roll angle gamma of the unmanned special operation vehicle can be calculated.
The specific process of the strapdown inertial navigation system of the embodiment for realizing speed update is as follows:
the velocity equation obtained from the specific force equation is:
Figure BDA0003742385530000186
in formula (39), f b Is the accelerometer output;
Figure BDA0003742385530000187
is the Coriolis acceleration;
Figure BDA0003742385530000188
for traction acceleration, speed updating can be completed through integral operation, and a specific expression is as follows:
Figure BDA0003742385530000189
in the formula (40), V m-1 The speed at the previous moment; Δ V smf Is a specific force integral increment; Δ V g/corm Velocity increments due to detrimental acceleration;
Figure BDA00037423855300001810
equations (39) to (41) are velocity updates for inertial navigation.
The specific process of implementing the position update by the strapdown inertial navigation system of the embodiment is as follows:
the longitude and latitude and the height of the unmanned vehicle are obtained through integral operation according to the speed of the unmanned automatic operation special vehicle under a navigation coordinate system, and the specific formula is as follows:
Figure BDA0003742385530000191
Figure BDA0003742385530000192
Figure BDA0003742385530000193
wherein L represents latitude, L m Represents the latitude at time m, L m-1 Represents the latitude at the m-1 moment; lambda m Represents longitude at time m; lambda [ alpha ] m-1 Represents longitude at time m-1; h is a total of m Representing the height of the unmanned vehicle at the m moment, h m-1 Representing the height of the unmanned vehicle at the m-1 moment; t is t m Represents time m, t m-1 Represents the m-1 time;
Figure BDA0003742385530000194
representing the speed of x direction to the power of n;
Figure BDA0003742385530000195
representing the speed n power in the y direction;
Figure BDA0003742385530000196
representing the velocity in the z direction to the power of n; r n Represents the radius of the earth; h represents the satellite altitude.
The quality of the positioning accuracy of the strapdown inertial navigation system depends on the effect of initial alignment, if the initial alignment is not performed, the larger the error of the starting point is, the worse the subsequent accuracy of the navigation system is, for this reason, the embodiment further includes an initial coarse alignment step, specifically: because the zero offset of the gyroscope is far greater than the rotation angular velocity of the earth, the gyroscope can hardly sense the rotation angular velocity of the earth, and therefore the magnitude of the initial course angle cannot be estimated; since accelerometers have a relatively high accuracy, the separation of gravitational acceleration that can be sensed by accelerometers allows the initial horizontal attitude angle to be estimated by the accelerometer: a pitch angle theta and a roll angle gamma;
Figure BDA0003742385530000197
the three elements on the left side of the formula (45) are the output of the accelerometer under the carrier coordinate system; and three elements in brackets on the right side are projections of the gravity acceleration under the navigation coordinate system, under the static condition, the gravity acceleration in the rough alignment generally has a projection of-g in the Z axis, and the projections of other axes are all 0, namely g n ＝[0 0 -g] T (ii) a It is known that
Figure BDA0003742385530000198
The belt (2.50) gives:
Figure BDA0003742385530000199
Figure BDA00037423855300001910
Figure BDA00037423855300001911
g is the gravity acceleration, and the heading angle psi is an unknown quantity, which can be solved reversely:
Figure BDA00037423855300001912
Figure BDA0003742385530000201
by adopting the method, the invention effectively improves the positioning accuracy and reliability by unifying the time system and the coordinate system of the BDS and the GPS, deducing the single-difference and double-difference observation equations, and combining the retrieval and the inspection of the ambiguity of the whole cycle and the detection and the repair of the cycle slip, and can greatly meet the requirements of the unmanned automatic operation special vehicle in practical application.
Although the invention has been described in detail with respect to the general description and the specific embodiments, it will be apparent to those skilled in the art that modifications and improvements may be made based on the invention. Accordingly, such modifications and improvements are intended to be within the scope of the invention as claimed.

Claims (10)

1. A high-precision positioning method for a special unmanned automatic operation vehicle based on network RTK is characterized by comprising the following steps:
step 1: establishing a GNSS differential observation equation, and eliminating errors such as satellite clock error, atmospheric delay error and receiver clock error;
and 2, step: estimating the integer ambiguity;
and step 3: performing decorrelation and search on the integer ambiguity;
and 4, step 4: checking the ambiguity;
and 5: detecting and repairing cycle slip;
step 6: the attitude information of the unmanned automatic operation special vehicle is sensed, meanwhile, the real-time motion measurement information is converted into a navigation coordinate system, and the position and the speed are updated.
2. The high-precision positioning method for the special unmanned automatic operation vehicle based on the network RTK as claimed in claim 1, characterized in that the specific process of step 1 is as follows:
step 1-1: the unified space-time reference comprises the unification of a time system and the unification of a coordinate system, wherein the unification of the time system is as follows:
GPST＝UTC+1s×n-19s (1)
BDT＝GPST-14s (2)
in the formula, when the GPST is a GPS, the GPS represents a special time system controlled by an atomic clock of a GPS main control station; UTC represents coordinated universal time; n represents an adjustment parameter, and a specific numerical value is issued by the international terrestrial rotation service bureau; BDT represents the time of big Dipper;
step 1-2: establishing a pseudo range and carrier phase observation equation:
Figure FDA0003742385520000011
Figure FDA0003742385520000012
in the formula, the superscript s represents the observation satellite and the subscript rRepresents a receiver; p and
Figure FDA0003742385520000013
respectively are pseudo range and carrier phase observed quantity; rho is the real distance between the satellite signal transmission time and the receiver; c is the speed of light in vacuum; dt r And dt s Respectively a receiver clock error and a satellite clock error; t is tropospheric error of the signal in the process of propagation in the atmosphere; i is an ionospheric error of a signal in the process of propagation in the atmosphere; λ is the carrier phase wavelength; n is the integer ambiguity in weeks; ε p and
Figure FDA0003742385520000021
other errors in pseudorange and carrier phase (including noise errors), respectively;
step 1-3: establishing a single difference observation model, carrying out inter-station difference between a reference station r and a mobile station q, and carrying out difference calculation on the same satellite s by the reference station receiver r and the mobile station receiver q, wherein the pseudo range and carrier phase observation equation of the reference station r is shown as (3) and (4), and the pseudo range and carrier phase observation equation of the mobile station q is as follows:
Figure FDA0003742385520000022
Figure FDA0003742385520000023
the difference between the formula (3) and the formula (5), and the difference between the formula (4) and the formula (6) can be the satellite clock difference dt q Thus, a single difference model is obtained:
Figure FDA0003742385520000024
Figure FDA0003742385520000025
in equations (7) and (8), Δ represents a single difference, and for the case of a short baseline, the ionospheric delay error and the tropospheric delay error of the reference station and the mobile station differ by a small amount, so that the difference between the tropospheric delay and the ionospheric delay in equations (7) and (8) can be approximated to be equal to 0, and it can be obtained:
Figure FDA0003742385520000026
Figure FDA0003742385520000027
step 1-4: establishing a double-difference observation model, and setting a reference satellite y, then the equation of the single difference of the reference station r and the mobile station q is similar to the equations (9) and (10),
Figure FDA0003742385520000028
Figure FDA0003742385520000029
by subtracting equations (9) and (11), (10) and (12), the receiver clock error term is eliminated to obtain a double-difference observation equation:
Figure FDA00037423855200000210
Figure FDA00037423855200000211
in formulas (13), (14),. DELTA.denotes a double difference sign, and double difference integer ambiguity retains the integer property and can be separated efficiently;
with N available satellites, the number of carrier phase and pseudorange observation equations is 2 (N-1), where the unknowns include the 3 directional coordinates of the baseline double-differenced pseudorange, and also the N-1 double-differenced ambiguities, so according to the principle of mathematically solving the equations: the equation number is more than or equal to the unknown number; to solve for the floating solution of ambiguities, then: 2 (N-1) is more than or equal to 3+ (N-1), and N =4 is obtained by solving, so that 4 visible satellites are observed at the same time, and the minimum requirement for realizing differential positioning by a single system is met;
when 4 satellites a, b, c and d are used for resolving, pseudo-range double-difference observed quantities are added, a 6-dimensional equation set can be established, and a satellite a with a large altitude angle is used as a main satellite, so that the following results are obtained:
Figure FDA0003742385520000031
where [ ln ] is the unit vector from the satellite to the receiver, [ δ X δ Y δ Z ] is the receiver coarse coordinate correction, [ Δ R is the double difference pseudorange observation, and equation (15) is written as a matrix:
Y＝HX (16)
in the formula (16), H is a design matrix formed by unit vectors from the receiver to the satellite, X is a base line vector to be solved and a ambiguity solution, and a float solution of ambiguity can be solved by using a least square method:
Figure FDA0003742385520000032
in the formula (17), the compound represented by the formula (I),
Figure FDA0003742385520000033
for floating-point solutions of ambiguities, let the base-line vector
Figure FDA0003742385520000034
Double difference integer ambiguity vector
Figure FDA0003742385520000035
Minimum sizeThe result of the second multiplication can be expressed as:
Figure FDA0003742385520000036
Figure FDA0003742385520000037
in equation (19), Q represents a covariance matrix of the vector X,
Figure FDA0003742385520000038
a covariance matrix representing the baseline vector,
Figure FDA0003742385520000039
a covariance matrix representing the ambiguity vector,
Figure FDA00037423855200000310
and
Figure FDA00037423855200000311
a cross-covariance matrix representing the baseline vector and the ambiguity vector.
3. The method for high-precision positioning of the special vehicle for unmanned automated operation based on network RTK according to claim 2, characterized in that the specific process of the step 2 is as follows:
step 2-1: solving floating point
Figure FDA00037423855200000312
Decomposing into the form of a baseline vector and an integer ambiguity vector to obtain:
Figure FDA00037423855200000313
in the formula (20), y is pseudo range and carrier phaseThe difference of the observed quantity is measured,
Figure FDA0003742385520000041
is a baseline vector, M is a coefficient matrix of the baseline vector, N is an integer ambiguity coefficient matrix,
Figure FDA0003742385520000042
is a double difference integer ambiguity, and epsilon is an unmodeled error quantity;
step 2-2: converting the solution ambiguity problem into a least square problem:
Figure FDA0003742385520000043
in the formula (21), Q y Is the covariance matrix of y;
step 2-3: obtaining integer ambiguity floating solution by formula (18) and formula (19)
Figure FDA0003742385520000044
Sum covariance matrix Q, then float solution with ambiguity vector
Figure FDA0003742385520000045
Sum covariance matrix
Figure FDA0003742385520000046
Integer solution for estimating ambiguities
Figure FDA0003742385520000047
Minimizing the sum of squared ambiguity residuals yields:
Figure FDA0003742385520000048
4. the high-precision positioning method for the special unmanned automated operation vehicle based on the network RTK as claimed in claim 3, wherein the specific process of the step 3 is as follows:
let transform matrix z satisfy:
Figure FDA0003742385520000049
then the corresponding covariance matrix is:
Figure FDA00037423855200000410
in the formula (24), Z is an integer matrix whose modulo is 1, all elements of the matrix are integers, and the modulo value is 1, | det (Z) | =1, and the correlation between the ambiguity and the target function Ω can be minimized by selecting an appropriate Z matrix and transforming;
Figure FDA00037423855200000411
obtaining an optimal solution of an integer matrix satisfying the formula (25), and setting a search space x 2 Namely:
Figure FDA00037423855200000412
the search space is a multi-dimensional ellipsoid, the shape and the covariance matrix
Figure FDA00037423855200000413
In this regard, the size is determined by the constant χ.
5. The high-precision positioning method for the special unmanned automated operation vehicle based on the network RTK as claimed in claim 4, wherein the specific process of the step 4 is as follows:
setting the minimum value of the searched residual error as omega min Then the residual sub-minimum is Ω sec Then, we can get:
Figure FDA00037423855200000414
in the formula (27), k is an empirical constant and cannot change with changes in model strength, measurement environment and the like, generally k =3 is used as a judgment condition for inspection, and ratio is greater than or equal to 3, and the optimal solution obtained by searching is considered to be a correct solution;
obtaining integer ambiguity fixing solution through ambiguity searching
Figure FDA0003742385520000051
By inverse transformation, a fixed solution of the original double-difference ambiguity can be obtained
Figure FDA0003742385520000052
Figure FDA0003742385520000053
Double-difference integer ambiguity vector obtained by integer ambiguity resolution
Figure FDA0003742385520000054
The integer characteristic of the ambiguity is used for further improving the estimation precision of the baseline vector, and the integer characteristic is obtained by LAMBDA search
Figure FDA0003742385520000055
It is possible to obtain:
Figure FDA0003742385520000056
Figure FDA0003742385520000057
in the above formula, the first and second carbon atoms are,
Figure FDA0003742385520000058
is a fixed solution to the baseline vector modifier,
Figure FDA0003742385520000059
is composed of
Figure FDA00037423855200000510
The covariance matrix of (2).
According to the above results, can obtain
Figure FDA00037423855200000511
Precise three-dimensional coordinates (X) of a rover q ,Y q ,Z q ) T Comprises the following steps:
(X q ,Y q ,Z q ) T ＝(X r0 ,Y r0 ,Z r0 ) T +(δX r ,δY r ,δZ r ) T (31)。
6. the method for high-precision positioning of the special unmanned automatic vehicle based on network RTK as claimed in claim 5, wherein the specific process of step 6 is:
step 6-1: establishing an SINS coordinate system of the unmanned automatic operation special vehicle, wherein the SINS coordinate system comprises a carrier coordinate system, an inertia coordinate system and a geographical coordinate system;
step 6-2: calculating attitude parameters of the unmanned automatic operation special vehicle by adopting a quaternion method, wherein the attitude parameters comprise a course angle, a pitch angle and a roll angle of the unmanned automatic operation special vehicle, and updating the attitude of the strapdown inertial navigation;
and 6-3: updating the speed of the strapdown inertial navigation;
step 6-4: and updating the position of the strapdown inertial navigation.
7. The high-precision positioning method for the special unmanned automated operation vehicle based on network RTK according to claim 6, characterized in that the specific process of the step 6-2 is as follows:
the quaternion is a four-dimensional array, is composed of four elements, and comprises three virtual units and a real unit, and the mathematical expression of the quaternion is as follows:
Figure FDA00037423855200000512
in the formula (32), q 0 、q 1 、q 2 And q is 3 Is any real number, and is a real number,
Figure FDA00037423855200000513
and
Figure FDA00037423855200000514
are mutually perpendicular unit vectors; the matrix of quaternions is:
Figure FDA0003742385520000061
in the formula (33), θ represents an angle,
Figure FDA0003742385520000062
the rotation matrix from which the quaternion representation can be derived is:
Figure FDA0003742385520000063
in the process of rotating from n to b, as three coordinate axes in two coordinate systems are mutually vertical in pairs,
Figure FDA0003742385520000064
and
Figure FDA0003742385520000065
show two reciprocityThe relationship between the coordinate system and the coordinate system is as follows:
Figure FDA0003742385520000066
in the formula (35), psi is a course angle of the unmanned special operating vehicle; theta is the pitch angle of the unmanned special operation vehicle; gamma is a roll angle of the unmanned special operation vehicle; the combined type (33) and the formula (34) are as follows:
Figure FDA0003742385520000067
γ＝arcsin(-C 13 )＝arcsin(-2(q 1 q 3 -q 0 q 2 )) (37)
Figure FDA0003742385520000068
observing the above formula can obtain the quaternion q 0 、q 1 、q 2 And q is 3 The values of the attitude angle psi, the pitch angle theta and the roll angle gamma of the unmanned special type operation vehicle can be calculated.
8. The method for high-precision positioning of special vehicle for unmanned automated operation based on network RTK according to claim 7, wherein the specific process of step 6-3 is
The velocity equation obtained from the specific force equation is:
Figure FDA0003742385520000069
in formula (39), f b Is the accelerometer output;
Figure FDA0003742385520000071
is the Coriolis acceleration;
Figure FDA0003742385520000072
for traction acceleration, speed updating can be completed through integral operation, and a specific expression is as follows:
Figure FDA0003742385520000073
in the formula (40), V m-1 The speed at the previous moment; Δ V smf Is the specific force integral increment; Δ V g/corm Velocity increments due to unwanted acceleration;
Figure FDA0003742385520000074
equations (39) to (41) are velocity updates for inertial navigation.
9. The high-precision positioning method for the special unmanned automated operation vehicle based on network RTK according to claim 8, characterized in that the specific process of the step 6-4 is as follows:
the longitude and latitude and the height of the unmanned vehicle are obtained through integral operation according to the speed of the unmanned automatic operation special vehicle under a navigation coordinate system, and the specific formula is as follows:
Figure FDA0003742385520000075
Figure FDA0003742385520000076
Figure FDA0003742385520000077
wherein L represents latitude, L m Represents the latitude at time m, L m-1 Represents the latitude at the m-1 moment; lambda [ alpha ] m Represents longitude at time m; lambda [ alpha ] m-1 Represents longitude at time m-1; h is a total of m Representing the height of the unmanned vehicle at the m moment, h m-1 Representing the height of the unmanned vehicle at the m-1 moment; t is t m Represents time m, t m-1 Represents the m-1 moment;
Figure FDA0003742385520000078
representing the speed of x direction to the power of n;
Figure FDA0003742385520000079
representing the speed n power in the y direction;
Figure FDA00037423855200000710
represents the speed n power in the z direction; r n Represents the radius of the earth; h represents the satellite altitude.
10. The method for high-precision positioning of the special unmanned automated operation vehicle based on the network RTK as claimed in claim 9, further comprising an initial coarse alignment step, specifically: because the zero offset of the gyroscope is far greater than the rotation angular velocity of the earth, the gyroscope can hardly sense the rotation angular velocity of the earth, and therefore the magnitude of the initial course angle cannot be estimated; since accelerometers have a relatively high accuracy, the separation of gravitational acceleration sensed by accelerometers estimates the initial horizontal attitude angle by the accelerometer: a pitch angle theta and a roll angle gamma;
Figure FDA0003742385520000081
the three elements on the left side of the formula (45) are the output of the accelerometer under the carrier coordinate system; and three elements in brackets on the right side are projections of the gravity acceleration under the navigation coordinate system, under the static condition, the gravity acceleration in the rough alignment generally has a projection of-g in the Z axis, and the projections of other axes are all 0, namely g n ＝[0 0 -g] T (ii) a It is known that
Figure FDA0003742385520000082
The belt (2.50) gives:
Figure FDA0003742385520000083
Figure FDA0003742385520000084
Figure FDA0003742385520000085
g is the gravity acceleration, and the heading angle psi is an unknown quantity, which can be solved reversely:
Figure FDA0003742385520000086
Figure FDA0003742385520000087
CN202210820920.2A 2022-07-12 2022-07-12 Network RTK-based high-precision positioning method for unmanned automatic operation special vehicle Pending CN115327588A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210820920.2A CN115327588A (en) 2022-07-12 2022-07-12 Network RTK-based high-precision positioning method for unmanned automatic operation special vehicle

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210820920.2A CN115327588A (en) 2022-07-12 2022-07-12 Network RTK-based high-precision positioning method for unmanned automatic operation special vehicle

Publications (1)

Publication Number Publication Date
CN115327588A true CN115327588A (en) 2022-11-11

Family

ID=83917942

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210820920.2A Pending CN115327588A (en) 2022-07-12 2022-07-12 Network RTK-based high-precision positioning method for unmanned automatic operation special vehicle

Country Status (1)

Country Link
CN (1) CN115327588A (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116030646A (en) * 2023-03-30 2023-04-28 南昌航天广信科技有限责任公司 Road condition information broadcasting method, system, computer and storage medium
CN116540278A (en) * 2023-07-06 2023-08-04 中国科学院空天信息创新研究院 Cloud edge end cooperative reference dynamic maintenance method and system
CN116699665A (en) * 2023-08-08 2023-09-05 山东科技大学 Unmanned ship positioning system and method suitable for offshore photovoltaic power plant environment
CN117250646A (en) * 2023-11-17 2023-12-19 毫厘智能科技(江苏)有限公司 Direction finding method and device based on chip, chip module and storage medium

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116030646A (en) * 2023-03-30 2023-04-28 南昌航天广信科技有限责任公司 Road condition information broadcasting method, system, computer and storage medium
CN116540278A (en) * 2023-07-06 2023-08-04 中国科学院空天信息创新研究院 Cloud edge end cooperative reference dynamic maintenance method and system
CN116540278B (en) * 2023-07-06 2023-09-08 中国科学院空天信息创新研究院 Cloud edge end cooperative reference dynamic maintenance method and system
CN116699665A (en) * 2023-08-08 2023-09-05 山东科技大学 Unmanned ship positioning system and method suitable for offshore photovoltaic power plant environment
CN117250646A (en) * 2023-11-17 2023-12-19 毫厘智能科技(江苏)有限公司 Direction finding method and device based on chip, chip module and storage medium
CN117250646B (en) * 2023-11-17 2024-02-02 毫厘智能科技(江苏)有限公司 Direction finding method and device based on chip, chip module and storage medium

Similar Documents

Publication Publication Date Title
CN101743453B (en) Post-mission high accuracy position and orientation system
Noureldin et al. Fundamentals of inertial navigation, satellite-based positioning and their integration
JP5989813B2 (en) Receiver positioning
CN115327588A (en) Network RTK-based high-precision positioning method for unmanned automatic operation special vehicle
US20090093959A1 (en) Real-time high accuracy position and orientation system
CN113203418B (en) GNSSINS visual fusion positioning method and system based on sequential Kalman filtering
CN111829512A (en) AUV navigation positioning method and system based on multi-sensor data fusion
CN111965685A (en) Low-orbit satellite/inertia combined navigation positioning method based on Doppler information
Yi et al. Tightly-coupled GPS/INS integration using unscented Kalman filter and particle filter
Chi et al. Enabling robust and accurate navigation for UAVs using real-time GNSS precise point positioning and IMU integration
Wen et al. 3D LiDAR aided GNSS real-time kinematic positioning
CN115683094A (en) Vehicle-mounted double-antenna tight coupling positioning method and system in complex environment
Du et al. Integration of PPP GPS and low cost IMU
CN108205151B (en) Low-cost GPS single-antenna attitude measurement method
Farkas et al. Small UAV’s position and attitude estimation using tightly coupled multi baseline multi constellation GNSS and inertial sensor fusion
Ellum et al. A mobile mapping system for the survey community
Hein et al. High‐Precision Kinematic GPS Differential Positioning and Integration of GPS with a Ring Laser Strapdown Inertial System
CN115220078A (en) GNSS high-precision positioning method and navigation method based on carrier phase difference
Consoli et al. A multi-antenna approach for UAV's attitude determination
Cannon et al. Development and testing of an integrated INS/GPS cross-linked system for sub-meter positioning of a CF-188 jet fighter
Kjørsvik et al. Tightly coupled precise point positioning and inertial navigation systems
Vana Continuous Urban Navigation with Next-Generation, Mass Market Navigation Sensors and Adaptive Filtering
CN117782080B (en) Real-time space-based navigation system and method based on PPP-B2B/INS
CN117031521B (en) Elastic fusion positioning method and system in indoor and outdoor seamless environment
O'Keefe et al. Comparing Multicarrier Ambiguity Resolution Methods for Geometry-Based GPS and Galileo Relative Positioning and Their Application to Low Earth Orbiting Satellite Attitude Determination.

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