CN108205151B - Low-cost GPS single-antenna attitude measurement method - Google Patents

Low-cost GPS single-antenna attitude measurement method Download PDF

Info

Publication number
CN108205151B
CN108205151B CN201810020784.2A CN201810020784A CN108205151B CN 108205151 B CN108205151 B CN 108205151B CN 201810020784 A CN201810020784 A CN 201810020784A CN 108205151 B CN108205151 B CN 108205151B
Authority
CN
China
Prior art keywords
receiver
data
satellite
difference
formula
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.)
Active
Application number
CN201810020784.2A
Other languages
Chinese (zh)
Other versions
CN108205151A (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.)
Chongqing University of Post and Telecommunications
Original Assignee
Chongqing University of Post and Telecommunications
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 Chongqing University of Post and Telecommunications filed Critical Chongqing University of Post and Telecommunications
Priority to CN201810020784.2A priority Critical patent/CN108205151B/en
Publication of CN108205151A publication Critical patent/CN108205151A/en
Application granted granted Critical
Publication of CN108205151B publication Critical patent/CN108205151B/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/53—Determining attitude
    • G01S19/54—Determining attitude using carrier phase measurements; using long or short baseline interferometry
    • G01S19/55—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)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention provides a low-cost GPS single-antenna attitude measurement method. According to the method, a double-difference carrier phase observation model is established by adopting a high-precision carrier phase as a main observation quantity, then an attitude measurement model is established by utilizing the relative positions of a single antenna at two different moments, and a single-antenna attitude measurement method is designed. Therefore, the attitude measurement cost is reduced, the defects that a multi-antenna structure is complex and easy to deform are overcome, and the application scene of GPS attitude measurement is effectively enlarged.

Description

Low-cost GPS single-antenna attitude measurement method
Technical Field
The invention belongs to the technical field of navigation, and relates to a low-cost GPS single-antenna attitude measurement method.
Background
Attitude information of a carrier is an important parameter of modern navigation, and the parameters are beneficial to better understanding the motion state of the carrier. Attitude and attitude measurement systems based on navigation satellites have been widely used in various fields such as land, sea, aviation, aerospace, and the like. Such as ships and airplanes, need quick and accurate attitude and azimuth results, and whether the results influence the navigation route.
With the development of the carrier phase differential technology, the attitude measurement by using the GPS system becomes possible. The traditional inertial navigation system can accurately calculate the attitude angle of the carrier without being interfered by the outside and has higher concealment. However, the method has error accumulation along with the change of time, and usually an external device is required to correct the output attitude angle, so that the obtained attitude angle is not completely high-precision. By utilizing the advantage of the ranging signal precision of a satellite signal system and adopting a carrier phase with higher precision as an observed quantity for attitude information calculation, the attitude measurement is realized, and the requirements of high precision, no error accumulation, small volume, low cost and the like in modern navigation can be met.
The measurement of the attitude of a large ship carrier under a long baseline condition by installing a plurality of antennas is realized. The multi-antenna attitude measurement has the defects of complex structure, large volume, large installation difficulty, easy deformation and the like, and is probably not suitable for being applied to scenes with small requirements on the volume of the antenna, such as unmanned vehicles, small unmanned aerial vehicles and the like.
Therefore, the invention designs a low-cost GPS single-antenna attitude measurement method, which can effectively reduce the volume of the antenna, thus facilitating the installation and reducing the cost. Meanwhile, the application scene of attitude measurement can be enlarged, and the method has very important significance.
Disclosure of Invention
The invention aims to overcome the existing defects and provide a low-cost GPS single-antenna attitude measurement method, which takes the defects of overlarge volume, inconvenience in installation, easiness in deformation and the like of an antenna during multi-antenna measurement into consideration when the GPS is favorable for attitude measurement, effectively reduces the number of the antennas, selects a carrier phase as a main observed quantity and improves the attitude measurement precision.
In order to achieve the purpose, the invention provides the following technical scheme:
a low-cost GPS single-antenna attitude measurement method comprises the following steps:
the method comprises the following steps: a data acquisition platform is set up to record the original data of the GPS board card;
step two: preprocessing ephemeris data and observation data in the raw data acquired in the step 1;
step three: and calculating the positions of the satellites observed by the receiver by using the ephemeris data in the step 2.
Step four: and calculating the position of the receiver by using the observation data in the step 2.
Step five: and (3) analyzing the observation data at two adjacent moments respectively by using the observation data in the step (2), and establishing a double-difference carrier phase observation model.
Step six: and 5, deducing a model suitable for single antenna attitude measurement according to the double-difference carrier phase observation model established in the step 5.
Step seven: and fast integer ambiguity searching is carried out based on a least square algorithm, and integer ambiguity is solved to further improve the baseline vector estimation precision.
Step eight: and solving the attitude angle according to the transformation of the baseline vector and the coordinate system.
Further, in the first step, the specific method for setting up the data acquisition platform to record the original data of the GPS board card is as follows: the GPS board card is driven by the singlechip, and the board card output frequency is accurately controlled by the timer; acquiring original data broadcast by a GPS satellite in real time, and encapsulating 1 acquired epoch data each time into one frame of data; and outputting the data to the outside through the Bluetooth interface after the data packaging is finished.
Further, in the second step, the preprocessing of the data acquired by the GPS board card includes unit conversion of the original data according to a data manual.
Further, in step three, the sub-step of calculating the position of the satellite using the observation data includes:
step three (one) of calculating the average angular velocity n of the satellite0
Figure BDA0001543466610000021
In the formula, μ is an earth gravity constant in an earth-centered earth coordinate system, and a is a long axis number provided by ephemeris data.
Step three (two) calculating the time difference t from the observation epoch to the reference epochk
tk＝t-toe
Where t is the time at which the ephemeris data is received by the receiver, toeA reference time provided for ephemeris data.
Step three (three) of calculating the corrected average angular rate n
n＝n0+Δn
In the formula, n0The delta n is the difference between the average satellite motion rate provided by the ephemeris data and the calculated value.
Step three (four) of calculating observation epoch tkFlat near point angle Mk
Mk＝M0+n·tk
Where n is the calculated correction averageAngular rate, M0Is a reference time toeFlat approach angle of (1).
Step three (five) of calculating the angle E of the approximate pointk
Ek＝Mk-e·sinEk
In the formula, MkFor already calculated observation epoch tkMean anomaly angle, e is the eccentricity given by the ephemeris data.
Step three (six) of calculating the true proximal angle vk
Figure BDA0001543466610000031
Where e is the eccentricity given by the ephemeris data. EkIs the already calculated angle of approach point.
Step three (seven) of calculating latitude argument parameter phik
φk＝vk+ω
In the formula, vkIs the true proximal angle already calculated. Omega is a reference time toeThe orbital angular distance is provided by ephemeris data.
Step three (eight) calculating correction item
Figure BDA0001543466610000032
uk＝φk+δuk
rk＝A·(1-e·cosEk)+δrk
ik＝i0+IDOT·tk+δik
In the formula, δ uk、δrk、δikRespectively latitude argument phikCorrection term, radial correction term and orbital inclination correction term, ukTo corrected latitude argument, rkFor the corrected mirror image ikFor the corrected track inclination, IDOT is the track inclination change rate.
Step three (nine) calculating the coordinates of the satellite in the orbit plane
Figure BDA0001543466610000033
In the formula ukTo corrected latitude argument, rkIs a corrected mirror image.
Step three (ten) of calculating the position of the satellite in the geocentric geodetic coordinate system
Figure BDA0001543466610000034
In the formula, omega0Is a reference epoch toeThe diameter of the rising point of (A) is red,
Figure BDA0001543466610000035
is a reference epoch toeThe perturbation rate of change in declination.
Figure BDA0001543466610000036
In the formula, xk，ykFor the calculated coordinates of the satellite in the orbital plane, ikFor the calculated corrected track inclination, ΩkIs the ascending crossing declination of epoch.
Further, in step four, the specific method for calculating the receiver position by observing data is as follows: and setting a pseudo-range observation value as rho, the true distance from the receiver to the satellite as R, the clock error of the receiver as delta t and c as the light speed, and establishing a pseudo-range observation equation:
ρ＝R+c·δt
if a total of n satellites are observed at a certain time, n pseudo-range observation equations can be established. The above equation can be linearized, yielding the following set of equations:
Figure BDA0001543466610000041
in the formula, l, m, and n are unit vectors of the satellite with respect to the receiver.
Figure BDA0001543466610000042
During initial calculation, the pseudorange single-point positioning is adopted to solve the approximate coordinate, the geometric distance R from the receiver to the satellite is not accurate, and therefore the difference between the calculated result and the position coordinate of the real receiver is large. And obtaining an accurate positioning result through repeated iterative calculation. In the actual measurement process, an observation error model is required to be established to correct errors such as an ionosphere, a troposphere, satellite clock errors and the like in the GPS signal propagation process.
Further, in step five, a double-difference carrier phase observation model is established by using observation data at two different moments. The specific selection method comprises the following steps:
according to the carrier phase in the satellite observation data, firstly establishing a carrier phase observation model which can be expressed as:
Figure BDA0001543466610000043
in the formula (I), the compound is shown in the specification,
Figure BDA0001543466610000044
for the receiver tkThe carrier phase value actually measured at the position i of the moment;
Figure BDA0001543466610000045
is tkThe actual distance of the time satellite p from the receiver position i; deltai(tk) Is the receiver clock difference at position i;
Figure BDA0001543466610000046
is the integer ambiguity of the satellite p; deltap(tk) Is the clock error of satellite p;
Figure BDA0001543466610000047
is an ionospheric error;
Figure BDA0001543466610000048
is tropospheric error; c is the speed of light and λ is the wavelength of the carrier.
Let the receiver be at tkThe approximate position of the time is
Figure BDA0001543466610000049
The correction number is [ delta X ]i(tk)δYi(tk)δZi(tk)]Is represented by the formula
Figure BDA00015434666100000410
And (3) performing Taylor series expansion, and taking a first-order approximate expression as follows:
Figure BDA0001543466610000051
in the formula (I), the compound is shown in the specification,
Figure BDA0001543466610000053
[l m n]the unit vectors from satellite to receiver are respectively expressed as:
Figure BDA0001543466610000054
then, the carrier phase observation models of the receiver at the position i and the position j are subjected to subtraction to obtain single differential carrier phase observation models at different positions, linearization is performed according to the above formula, and the following formula can be obtained:
Figure BDA0001543466610000055
in the formula, deltaij(Δt)＝δi(tk)-δj(tm)，
Figure BDA0001543466610000057
Is the unit vector of the satellite p to the receiver j,
Figure BDA0001543466610000058
is the correction of the approximate coordinates of the receiver at position j, Δ t is the time difference between the receiver at position i and position j, tmThe time instant at which the receiver is at position j.
Finally, on the basis of the single-difference carrier phase observation model, performing single difference between the satellite p and the satellite q to obtain a double-difference carrier phase observation model which can be expressed as:
Figure BDA0001543466610000059
in the formula
Figure BDA00015434666100000510
The integer ambiguity after double differencing.
Further, in the sixth step, according to the double-difference carrier phase observation model, when 4 satellites are selected for resolving, 6 double-difference observation equations can be established, so that the single-antenna attitude measurement model is deduced. Selecting the satellite with the largest altitude angle as the reference satellite can obtain a simplified equation set:
Figure BDA0001543466610000061
the formula is abbreviated as:
Y＝A·X
in the formula: y is the observation data of the receiver; a is a design matrix formed by unit vectors from a receiver to a satellite, and X is a base line vector and a ambiguity solution to be solved.
According to the least square solution principle, X can be obtained:
X＝(AT·P·A)-1·AT·P·Y
neglecting the integer property of double difference integer ambiguity, making double difference integer ambiguity vector
Figure BDA0001543466610000062
Base line correction vector b ═ δ Xj δYj δZj]TThe result of the least squares solution can be expressed as:
Figure BDA0001543466610000063
Figure BDA0001543466610000064
commonly referred to as floating point solution, the variance and covariance matrix Q is:
Figure BDA0001543466610000065
in the formula (I), the compound is shown in the specification,
Figure BDA0001543466610000066
is an auto-covariance matrix of b,
Figure BDA0001543466610000067
is an auto-covariance matrix of a,
Figure BDA0001543466610000068
and
Figure BDA0001543466610000069
is the cross covariance matrix of a and b, and P is the weight matrix.
Further, in step seven, a specific method for performing fast ambiguity search based on least squares is as follows: obtaining a double-difference integer ambiguity vector by integer ambiguity resolution
Figure BDA00015434666100000610
And the integer characteristic of the ambiguity is utilized to further improve the baseline vector estimation precision. Derived from least squares search
Figure BDA00015434666100000611
The following can be obtained:
Figure BDA00015434666100000612
Figure BDA00015434666100000613
in the formula (I), the compound is shown in the specification,
Figure BDA00015434666100000614
is a fixed solution to the baseline vector modifier,
Figure BDA00015434666100000617
is composed of
Figure BDA00015434666100000615
The corresponding covariance matrix.
Further, in step eight, the specific method for solving the attitude angle according to the transformation between the baseline coordinate and the coordinate system is as follows: due to the fact that
Figure BDA00015434666100000616
For the correction of the baseline vector, let the coordinate of the receiver's position i in the geocentric geodetic coordinate system be [ Xm Ym Zm]TThen the coordinates of the receiver at position j relative to the coordinates of position i are:
Figure BDA0001543466610000071
in the formula (I), the compound is shown in the specification,
Figure BDA0001543466610000072
the coordinate of the receiver at position j is corrected by the amount already calculated.
Figure BDA0001543466610000073
The coordinates of the receiver at position i.
Assuming that the position i of the receiver is the origin of coordinates, the amount of coordinate correction and the baseline vector correction at the position j of the receiver
Figure BDA0001543466610000074
The equivalence, namely:
Figure BDA0001543466610000075
due to the base line
Figure BDA0001543466610000076
Is based on the geocentric geodetic coordinate system
Figure BDA0001543466610000077
Converting into a standing center coordinate system to obtain
Figure BDA0001543466610000078
Figure BDA0001543466610000079
Where λ is the longitude of the receiver at location j and φ is the latitude of the receiver at location j.
The coordinate vector in the station center coordinate system and the coordinate vector in the earth center geodetic coordinate system have the following conversion relationship:
Figure BDA00015434666100000710
in the formula (I), the compound is shown in the specification,
Figure BDA00015434666100000711
is a coordinate vector in the station-centric coordinate system,
Figure BDA00015434666100000712
is a coordinate vector in the geocentric geodetic coordinate system,
Figure BDA00015434666100000713
is a transformation matrix of a station center coordinate system and a geocentric geodetic coordinate system.
When the base line is connected withWhen the y-axis of the carrier is coincident, the course angle of the carrier can be measured
Figure BDA00015434666100000714
In this case, the normalized vector of the baseline in the carrier coordinate system can be expressed as [010 [ ]]TThe roll angle is set to 0 °.
Figure BDA00015434666100000715
The course angle can be obtained by the following formula:
Figure BDA0001543466610000081
the invention has the beneficial effects that: according to the method, a double-difference carrier phase observation model is established by adopting a high-precision carrier phase as a main observation quantity, and then an attitude measurement model based on a single antenna is designed. The method and the device reduce the cost of attitude measurement, overcome the defects of complex multi-antenna structure and easy deformation, and effectively expand the application scene of GPS attitude measurement.
Drawings
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be further described in detail with reference to the accompanying drawings, in which:
FIG. 1 is a schematic flow chart of a GPS single antenna attitude measurement method;
FIG. 2 is a schematic view of a carrier phase observation model;
FIG. 3 is a schematic diagram of a carrier phase single difference model;
FIG. 4 is a schematic diagram of a carrier phase double difference model;
Detailed Description
So that those skilled in the art can better understand the objects, aspects and advantages of the present invention, a full description of the invention, including the detailed description, can be had by referring to the accompanying drawings.
FIG. 1 is a schematic flow chart of a GPS single antenna attitude measurement method;
the method comprises the following steps: a singlechip is adopted to drive a GPS board card and accurately control the board card output frequency through a timer; acquiring original data broadcast by a GPS satellite in real time, and encapsulating 1 acquired epoch data each time into one frame of data; and outputting the data to the outside through the Bluetooth interface after the data packaging is finished.
Step two: the original data were unit transformed according to the data manual.
Step three: calculating the positions of the satellites according to the ephemeris data, and the substep of the step is as follows:
step three (one) of calculating the average angular velocity n of the satellite0
Figure BDA0001543466610000082
In the formula, μ is an earth gravity constant in an earth-centered earth-ground coordinate system, and a is a long axis number in the ephemeris parameters.
Step three (two) calculating the time difference t from the observation epoch to the reference epochk
tk＝t-toe
Where t is the receiver time at which the ephemeris data was received, toeIs the reference time provided by the ephemeris data.
Step three (three) of calculating the corrected average angular rate n
n＝n0+Δn
In the formula, n0For the satellite mean angular velocity, Δ n is the ephemeris parameter
Step three (four) of calculating observation epoch tkFlat near point angle Mk
Mk＝M0+n·tk
Wherein n is the corrected average angular rate, M0For reference time t given in ephemeris dataoeFlat approach angle of (1).
Step three (five) of calculating the angle E of the approximate pointk
Ek＝Mk-e·sinEk
In the formula, MkFor already calculated observation epoch tkMean anomaly angle, e eccentricity, given by ephemeris data.
Step three (six) of calculating the true proximal angle vk
Figure BDA0001543466610000091
Where e is the eccentricity given by the ephemeris data. EkFor calculated angle of approach point
Step three (seven) of calculating latitude argument parameter phik
φk＝vk+ω
Where ω is the reference time toeAngular distance of track. v. ofkIs the true proximal angle already calculated.
Step three (eight) calculating correction term delta uk，δrk，δik
Figure BDA0001543466610000092
uk＝φk+δuk
rk＝A·(1-e·cosEk)+δrk
ik＝i0+IDOT·tk+δik
In the formula, δ uk、δrk、δikRespectively latitude argument phikCorrection term, radial correction term and orbital inclination correction term, ukTo corrected latitude argument, rkFor the corrected mirror image ikFor the corrected track inclination, IDOT is the track inclination change rate.
Step three (nine) calculating the coordinates of the satellite in the orbit plane
Figure BDA0001543466610000093
In the formula ukTo corrected latitude argument, rkIs a corrected mirror image.
Step three (ten) of calculating the position of the satellite in the geocentric geodetic coordinate system, firstly calculating the ascent intersection declination diameter omega of the epochk
Figure BDA0001543466610000101
In the formula, omega0Is a reference epoch toeThe diameter of the rising point of (A) is red,
Figure BDA0001543466610000102
is a reference epoch toeThe perturbation rate of change in declination.
Figure BDA0001543466610000103
In the formula, xk，ykFor the calculated coordinates of the satellite in the orbital plane, ikFor the calculated corrected track inclination, ΩkIs the ascending crossing declination of epoch.
Step four: and setting a pseudo-range observation value as rho, a real distance from a receiver to a satellite as R, and a receiver clock error as deltat. Pseudorange observation equations may be established:
ρ＝R+c·δt
if n satellites are observed at a certain time, n pseudo-range observation equations can be established and linearized to obtain the following equation set:
Figure BDA0001543466610000104
in the formula, l, m, and n are unit vectors of the satellite with respect to the receiver. The formula can be solved as follows:
Figure BDA0001543466610000105
during initial calculation, the pseudorange single-point positioning is adopted to solve the approximate coordinate, the geometric distance R from the receiver to the satellite is not accurate, and therefore the difference between the calculated result and the position coordinate of the real receiver is large. And obtaining an accurate positioning result through repeated iterative calculation. In the actual measurement process, an observation error model is required to be established to correct errors such as an ionosphere, a troposphere, satellite clock errors and the like in the GPS signal propagation process.
Step five: firstly, a carrier phase observation model is established, and at t0To tkThe signal transmission between the satellite p and the receiver during this time is schematically illustrated in fig. 2. At t0The time receiver can only measure
Figure BDA0001543466610000106
I.e. the fractional part, the integer part of the signal propagation process
Figure BDA0001543466610000107
It could not be measured. From t0The receiver starts to count continuously when the time begins, and can measure t as long as the signal is not unlocked0To tkThe number of whole cycles of satellite signal propagation during this period
Figure BDA0001543466610000111
And fractional part
Figure BDA0001543466610000112
Thus, at t0The total variation of the carrier phase of any epoch can be expressed as:
Figure BDA0001543466610000113
in the formula (I), the compound is shown in the specification,
Figure BDA0001543466610000114
for receiver reference signal at tkThe phase of the time, i.e. the phase of the satellite signal at that time;
Figure BDA0001543466610000115
is that the receiver is at tkThe phase of the satellite signal received by the time receiver;
Figure BDA0001543466610000116
is that the receiver is from t0Time begins to tkThe carrier phase variation recorded at the moment is also the actual observed quantity output by the receiver;
Figure BDA0001543466610000117
is t0Time-of-day carrier integer ambiguity.
If the receiver observes the satellite p at two different positions i and j, two carrier phase observation equations with the above formula can be established, and a schematic diagram of a carrier phase single difference model is shown in fig. 3.
Since the two receiver locations are closely spaced, the transmission paths of the signals of satellite p to receiver location i and location j can be approximately considered to be identical. Therefore, most path propagation errors such as ionospheric and tropospheric errors can be eliminated by differencing. The following equation can be obtained after a single difference:
Figure BDA0001543466610000118
in the formula (I), the compound is shown in the specification,
Figure BDA0001543466610000119
representing the single-differenced value of the carrier phase observations.
Figure BDA00015434666100001110
Is the receiver at position i time tkThe amount of change in the carrier phase recorded at a time,
Figure BDA00015434666100001111
is the receiver at position j, tmThe amount of change in the carrier phase recorded at a time,
in the single difference carrier phase observation model, the receiver's coordinates at location i are assumed to be known and the receiver's rough coordinate at location j is assumed to be knownIs marked as
Figure BDA00015434666100001116
Is in the above formula
Figure BDA00015434666100001117
Taylor expansion is carried out to obtain a single difference carrier phase observation equation after linearization:
Figure BDA00015434666100001112
in the formula, deltaij(Δt)＝δi(tk)-δj(tm)，
Figure BDA00015434666100001114
Is the unit vector of the satellite p to the receiver j,
Figure BDA00015434666100001115
is the correction of the approximate coordinates of the receiver at position j, Δ t is the time difference between the receiver at position i and position j, tmThe time instant at which the receiver is at position j.
In the above equation, the single difference carrier phase model has removed most of the error, but the receiver clock difference still exists. On the basis of the single-difference model, the difference is performed again for different satellites to obtain a double-difference carrier phase model, and a schematic diagram of the carrier phase double-difference model is shown in fig. 4. Assuming that a satellite p is a reference satellite, performing single difference on the satellite p and a satellite q at a position i and a position j by a receiver respectively, firstly obtaining single-difference carrier phase observation models between different positions, and then performing single difference between different satellites to obtain a double-difference carrier phase observation model, wherein the following formula is as follows:
Figure BDA0001543466610000121
in the formula (I), the compound is shown in the specification,
Figure BDA0001543466610000122
the integer ambiguity after double differencing.
Step six: according to the double-difference carrier phase model principle, when 4 satellites are selected for resolving, 6 double-difference observation equations are established, and a simplified equation set can be obtained:
Figure BDA0001543466610000123
the formula is abbreviated as:
Y＝A·X
in the formula: y is the observation data of the receiver; a is a design matrix formed by unit vectors from a receiver to a satellite; and X is the baseline vector and ambiguity solution to be solved.
According to the least square solution principle, X can be obtained:
X＝(AT·P·A)-1·AT·P·Y
neglecting the integer property of double difference integer ambiguity, making double difference integer ambiguity vector
Figure BDA0001543466610000124
Base line correction vector b ═ δ Xj δYj δZj]TThe result of the least squares solution can be expressed as:
Figure BDA0001543466610000125
Figure BDA0001543466610000126
commonly referred to as floating point solution, the variance and covariance matrix Q is:
Figure BDA0001543466610000127
in the formula (I), the compound is shown in the specification,
Figure BDA0001543466610000128
the autocovariance matrix of b，
Figure BDA0001543466610000129
Is an auto-covariance matrix of a,
Figure BDA00015434666100001210
and
Figure BDA00015434666100001211
is the cross covariance matrix of a and b, and P is the weight matrix.
Step seven: the specific method for rapidly searching the ambiguity based on the least square algorithm comprises the following steps: obtaining a double-difference integer ambiguity vector by integer ambiguity resolution
Figure BDA0001543466610000131
And the integer characteristic of the ambiguity is utilized to further improve the baseline vector estimation precision. Derived from least squares search
Figure BDA0001543466610000132
The following can be obtained:
Figure BDA0001543466610000133
Figure BDA0001543466610000134
in the formula (I), the compound is shown in the specification,
Figure BDA0001543466610000135
is a fixed solution to the baseline vector modifier,
Figure BDA0001543466610000136
is composed of
Figure BDA0001543466610000137
The corresponding covariance matrix is then used as a basis,
Figure BDA0001543466610000138
is an auto-covariance matrix of a,
Figure BDA0001543466610000139
and
Figure BDA00015434666100001310
is the cross covariance matrix of a and b, and P is the weight matrix.
Step eight: the specific method for solving the attitude angle according to the transformation of the baseline coordinate and the coordinate system comprises the following steps: due to the fact that
Figure BDA00015434666100001311
For the correction of the baseline vector, let the coordinate of the receiver's position i in the geocentric geodetic coordinate system be [ Xm Ym Zm]TThen the coordinates of the receiver at position j relative to the coordinates of position i are:
Figure BDA00015434666100001312
in the formula (I), the compound is shown in the specification,
Figure BDA00015434666100001313
for the already calculated amount of coordinate correction of the receiver at position j,
Figure BDA00015434666100001314
the coordinates of the receiver at position i.
Assuming that the position i of the receiver is the origin of coordinates, the amount of coordinate correction and the baseline vector correction at the position j of the receiver
Figure BDA00015434666100001315
The equivalence, namely:
Figure BDA00015434666100001316
due to the base line
Figure BDA00015434666100001317
Is based on the geocentric geodetic coordinate system
Figure BDA00015434666100001318
Converting into a standing center coordinate system to obtain
Figure BDA00015434666100001319
Figure BDA00015434666100001320
Where λ is the longitude of the receiver at location j and φ is the latitude of the receiver at location j.
The coordinate vector in the station center coordinate system and the coordinate vector in the earth center geodetic coordinate system have the following conversion relationship:
Figure BDA0001543466610000141
in the formula (I), the compound is shown in the specification,
Figure BDA0001543466610000142
is a coordinate vector in the station-centric coordinate system,
Figure BDA0001543466610000143
is a coordinate vector in the geocentric geodetic coordinate system,
Figure BDA0001543466610000144
is a transformation matrix of a station center coordinate system and a geocentric geodetic coordinate system.
When the baseline coincides with the y-axis of the carrier, the heading angle of the carrier can be measured
Figure BDA0001543466610000145
In this case, the normalized vector of the baseline in the carrier coordinate system can be expressed as [010 [ ]]TThe roll angle is set to 0 °.
Figure BDA0001543466610000146
The course angle can be obtained by the following formula:
Figure BDA0001543466610000147

Claims (6)

1. a low-cost GPS single antenna attitude measurement method is characterized in that: comprises the following steps:
step one, a data acquisition platform is set up to record original data of a GPS board card;
step two, preprocessing ephemeris data and observation data in the raw data acquired in the step one;
step three, calculating the position of the satellite observed by the receiver by using the ephemeris data in the step two;
step four, calculating the position of the receiver by using the observation data in the step two;
analyzing observation data of two adjacent moments respectively by using the observation data in the step two, and establishing a double-difference carrier phase observation model;
step six, deducing a measurement model suitable for the single antenna attitude according to the double-difference carrier phase observation model established in the step five;
seventhly, fast ambiguity searching is carried out based on a least square algorithm, and ambiguity is solved to further improve baseline vector estimation accuracy;
step eight, solving an attitude angle according to the conversion of the baseline vector and the coordinate system;
in step five, according to the carrier phase in the satellite observation data, firstly, a carrier phase observation model is established, which can be expressed as:
Figure FDA0003522718280000011
in the formula:
Figure FDA0003522718280000012
for the receiver tkThe actual measured carrier phase value at the location i of the time instant,
Figure FDA0003522718280000013
is tkThe actual distance, delta, of the time satellite p from the receiver position ii(tk) For the receiver clock difference at position i,
Figure FDA0003522718280000014
is the integer ambiguity, delta, of the satellite pp(tk) For the clock offset of the satellite p,
Figure FDA0003522718280000015
in order to be an ionospheric error,
Figure FDA0003522718280000016
is tropospheric error, c is the speed of light, λ is the wavelength of the carrier;
then, an inter-station single-difference carrier phase observation model can be obtained by carrying out difference on the observed quantities of the receiver at two different moments, satellites observed by the receiver at two different moments are selected as reference quantities, and a double-difference carrier phase observation model can be obtained by carrying out single difference between the planets in the inter-station single-difference carrier phase observation model;
in the sixth step, according to the double-difference carrier phase observation model, when 4 satellites are selected for resolving, 6 double-difference observation equations can be established, so as to derive a single antenna attitude measurement model, as shown in the following formula:
Y＝A·X
in the formula: y is a known observed quantity, A is a design matrix formed by unit vectors from a receiver to a satellite, and X is a baseline vector to be solved and a ambiguity solution;
x can be obtained according to the least square solution principle:
X＝(AT·P·A)-1·AT·P·Y
in the formula, P is a weight matrix;
neglecting the integer property of double difference integer ambiguity, making double difference integer ambiguity vector
Figure FDA0003522718280000021
Base line correction vector b ═ δ Xj δYj δZj]TThe result of the least squares solution can be expressed as:
Figure FDA0003522718280000022
the following equation can be obtained after a single difference:
Figure FDA0003522718280000023
in the formula (I), the compound is shown in the specification,
Figure FDA0003522718280000024
representing a single-differenced value of the carrier phase observed quantity;
Figure FDA0003522718280000025
is the receiver at position i time tkThe amount of change in the carrier phase recorded at a time,
Figure FDA0003522718280000026
is the receiver at position j, tmThe amount of change in the carrier phase recorded at a time,
in the single-difference carrier-phase observation model, the coordinates of the receiver at position i are assumed to be known, and the approximate coordinates of the receiver at position j are assumed to be
Figure FDA00035227182800000212
Is in the above formula
Figure FDA00035227182800000211
Taylor expansion is carried out to obtain a single difference carrier phase observation party after linearizationThe process:
Figure FDA0003522718280000027
in the formula, deltaij(Δt)＝δi(tk)-δj(tm)，
Figure FDA0003522718280000028
Is the unit vector of the satellite p to the receiver j,
Figure FDA0003522718280000029
for receivers at positions j, tmThe correction of the approximate coordinates of the time, Δ t being the time difference between the receiver at position i and at position j, tmThe time instant of the receiver at position j;
step three: calculating the positions of the satellites according to the ephemeris data, and the substep of the step is as follows:
step three (one) of calculating the average angular velocity n of the satellite0
Figure FDA00035227182800000210
In the formula, mu is an earth gravity constant under an earth-centered earth-ground coordinate system, and A is a long semi-axis number in the forecast ephemeris parameter;
step three (two) calculating the time difference t from the observation epoch to the reference epochk
tk＝t-toe
Where t is the receiver time at which the ephemeris data was received, toeA reference time provided by ephemeris data;
step three (three) of calculating the corrected average angular rate n
n＝n0+Δn
In the formula, n0For the satellite mean angular velocity, Δ n is the ephemeris parameter
Step three (four) of calculating observation epoch tkNear and flatPoint angle Mk
Mk＝M0+n·tk
Wherein n is the corrected average angular rate, M0For reference time t given in ephemeris dataoeA flat proximal angle of;
step three (five) of calculating the angle E of the approximate pointk
Mk＝Ek-e·sinEk
In the formula, MkFor already calculated observation epoch tkMean anomaly, e is eccentricity, given by ephemeris data;
step three (six) of calculating the true proximal angle vk
Figure FDA0003522718280000031
Wherein e is eccentricity and is given by ephemeris data; ekFor calculated angle of approach point
Step three (seven) of calculating latitude argument parameter phik
φk＝vk+ω
Where ω is the reference time toeAngular distance of the track from the near place; v. ofkIs the true proximal angle already calculated;
step three (eight) calculating correction term delta uk，δrk，δik
uk＝φk+δuk
rk＝A·(1-e·cosEk)+δrk
ik＝i0+IDOT·tk+δik
In the formula, δ uk、δrk、δikRespectively latitude argument phikCorrection term, radial correction term and orbital inclination correction term, ukTo corrected latitude argument, rkFor the corrected mirror image ikFor the corrected track inclination angle, IDOT is the track inclination angle change rate;
step three (nine) calculating the coordinates of the satellite in the orbit plane
Figure FDA0003522718280000032
In the formula ukTo corrected latitude argument, rkIs a corrected mirror image;
step three (ten) of calculating the position of the satellite in the geocentric geodetic coordinate system, firstly calculating the ascent intersection declination diameter omega of the epochk
Figure FDA0003522718280000041
In the formula, omega0Is a reference epoch toeThe diameter of the rising point of (A) is red,
Figure FDA0003522718280000042
is a reference epoch toePerturbation rate of change of declination;
Figure FDA0003522718280000043
in the formula, xk，ykFor the calculated coordinates of the satellite in the orbital plane, ikFor the calculated corrected track inclination, ΩkIs the ascent diameter of the epoch ascending intersection point;
step seven: the specific method for rapidly searching the ambiguity based on the least square algorithm comprises the following steps: obtaining a double-difference integer ambiguity vector by integer ambiguity resolution
Figure FDA0003522718280000044
Derived from least squares search
Figure FDA0003522718280000045
The following can be obtained:
Figure FDA0003522718280000046
Figure FDA0003522718280000047
in the formula (I), the compound is shown in the specification,
Figure FDA0003522718280000048
is an auto-covariance matrix of b,
Figure FDA0003522718280000049
is an auto-covariance matrix of a,
Figure FDA00035227182800000410
is a fixed solution to the baseline vector modifier,
Figure FDA00035227182800000411
is composed of
Figure FDA00035227182800000412
The corresponding covariance matrix is then used as a basis,
Figure FDA00035227182800000413
is an auto-covariance matrix of a,
Figure FDA00035227182800000414
and
Figure FDA00035227182800000415
is the cross covariance matrix of a and b, and P is the weight matrix.
2. The low-cost GPS single-antenna attitude measurement method according to claim 1, characterized in that: in the first step, the singlechip drives the GPS board card and accurately controls the board card to output frequency through a timer; acquiring original data broadcast by a GPS satellite in real time, and encapsulating 1 acquired epoch data each time into one frame of data; and outputting the data to the outside through the Bluetooth interface after the data packaging is finished.
3. The low-cost GPS single-antenna attitude measurement method according to claim 2, characterized in that: in the second step, preprocessing the data acquired by the GPS board card comprises unit conversion of the original data according to a data manual.
4. The low-cost GPS single-antenna attitude measurement method according to claim 3, characterized in that: in step three, the positions of the satellites observed by the receiver are calculated according to the ephemeris data in the raw data.
5. The low-cost GPS single-antenna attitude measurement method according to claim 4, characterized in that: in the fourth step, a pseudo-range observation value is set as rho, the real distance from the receiver to the satellite is set as R, the receiver clock error is set as deltat, and a pseudo-range observation equation can be established:
ρ＝R+c·δt
setting a certain time to observe n satellites in total, then establishing n pseudo-range observation equations, and then solving an equation set; during initial calculation, pseudo-range point positioning is adopted to solve the approximate coordinate, the geometric distance R from the receiver to the satellite is not accurate, and therefore the difference between the calculated result and the position coordinate of the real receiver is large; and obtaining an accurate positioning result through repeated iterative calculation.
6. The low-cost GPS single-antenna attitude measurement method according to claim 1, characterized in that: in step eight, the change of the relative position of the receiver at two different moments can be deduced according to the baseline vector correction number solved in step seven, and meanwhile, the heading angle can be deduced by combining the conversion relation between the coordinate vector in the geocentric geodetic coordinate system and the coordinate vector in the standing-center coordinate system.
CN201810020784.2A 2018-01-10 2018-01-10 Low-cost GPS single-antenna attitude measurement method Active CN108205151B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810020784.2A CN108205151B (en) 2018-01-10 2018-01-10 Low-cost GPS single-antenna attitude measurement method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810020784.2A CN108205151B (en) 2018-01-10 2018-01-10 Low-cost GPS single-antenna attitude measurement method

Publications (2)

Publication Number Publication Date
CN108205151A CN108205151A (en) 2018-06-26
CN108205151B true CN108205151B (en) 2022-05-03

Family

ID=62606323

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810020784.2A Active CN108205151B (en) 2018-01-10 2018-01-10 Low-cost GPS single-antenna attitude measurement method

Country Status (1)

Country Link
CN (1) CN108205151B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110146050B (en) * 2019-02-18 2021-03-30 广东星舆科技有限公司 Communication base station antenna monitoring method
CN110133702B (en) * 2019-05-13 2022-12-27 桂林电子科技大学 Attitude measurement method and equipment based on orthogonal transformation
CN110646821A (en) * 2019-09-26 2020-01-03 北京交通大学 Train integrity detection method based on moving baseline calculation

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002040124A (en) * 2000-07-24 2002-02-06 Furuno Electric Co Ltd Measuring device for carrier phase relative position
CN102096084A (en) * 2010-12-09 2011-06-15 东南大学 Precise point positioning (PPP) method based on inter-satellite combination difference
CN106772478A (en) * 2016-11-11 2017-05-31 哈尔滨工程大学 The localization method of difference constraint between a kind of star based on epoch
CN107037464A (en) * 2017-05-24 2017-08-11 陈湘南 A kind of accident vehicle precision positioning method based on GNSS relative positionings
CN107356952A (en) * 2017-07-04 2017-11-17 北京航空航天大学 It is a kind of that the high-precision Relative Navigation based on GNSS is independently carried out using single-receiver

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5185610A (en) * 1990-08-20 1993-02-09 Texas Instruments Incorporated GPS system and method for deriving pointing or attitude from a single GPS receiver
JP2003232845A (en) * 2002-02-12 2003-08-22 Furuno Electric Co Ltd Detection device of azimuth and attitude of moving body
US7292185B2 (en) * 2005-10-04 2007-11-06 Csi Wireless Inc. Attitude determination exploiting geometry constraints
US8638257B2 (en) * 2009-10-15 2014-01-28 Novatel, Inc. Ultra short baseline GNSS receiver
CN105607106B (en) * 2015-12-18 2018-08-21 重庆邮电大学 A kind of low-cost and high-precision BD/MEMS fusions attitude measurement method
CN105424041A (en) * 2016-01-18 2016-03-23 重庆邮电大学 Pedestrian positioning algorithm based on BD/INS (Beidou/Inertial Navigation System) tight coupling

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002040124A (en) * 2000-07-24 2002-02-06 Furuno Electric Co Ltd Measuring device for carrier phase relative position
CN102096084A (en) * 2010-12-09 2011-06-15 东南大学 Precise point positioning (PPP) method based on inter-satellite combination difference
CN106772478A (en) * 2016-11-11 2017-05-31 哈尔滨工程大学 The localization method of difference constraint between a kind of star based on epoch
CN107037464A (en) * 2017-05-24 2017-08-11 陈湘南 A kind of accident vehicle precision positioning method based on GNSS relative positionings
CN107356952A (en) * 2017-07-04 2017-11-17 北京航空航天大学 It is a kind of that the high-precision Relative Navigation based on GNSS is independently carried out using single-receiver

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"卫星/惯性/视觉组合导航信息融合关键技术研究";董明;《万方》;20160914;第53-69页 *

Also Published As

Publication number Publication date
CN108205151A (en) 2018-06-26

Similar Documents

Publication Publication Date Title
CN102230971B (en) GPS multi-antenna attitude determination method
CN106990424B (en) Double-antenna GPS attitude measurement method
US7292185B2 (en) Attitude determination exploiting geometry constraints
CN107193028B (en) Kalman relative positioning method based on GNSS
CN103675861B (en) Satellite autonomous orbit determination method based on satellite-borne GNSS multiple antennas
CN105425261B (en) Integrated navigation and localization method based on GPS/Beidou2/INS
CN101743453A (en) The post-mission high accuracy position and azimuth determining system
CN106255065A (en) Smart mobile phone and the seamless alignment system of mobile terminal indoor and outdoor and method thereof
CN108845345B (en) Double-antenna directional attitude measurement method based on GNSS speed measurement principle
CN108205151B (en) Low-cost GPS single-antenna attitude measurement method
CN111399020A (en) Directional attitude measurement system and method
CN112731502B (en) Unmanned aerial vehicle aerial refueling inertia-assisted Beidou three-frequency precise relative navigation method
CN110764127A (en) Relative orbit determination method for formation satellite easy for satellite-borne on-orbit real-time processing
CN115327588A (en) Network RTK-based high-precision positioning method for unmanned automatic operation special vehicle
Vana et al. Low-cost, dual-frequency PPP GNSS and MEMS-IMU integration performance in obstructed environments
CN115220078A (en) GNSS high-precision positioning method and navigation method based on carrier phase difference
CN117192578A (en) Shipborne measurement and control antenna shafting parameter calibration method for tracking unmanned aerial vehicle
CN115561796A (en) Real-time positioning method and system for power grid unmanned aerial vehicle routing inspection
CN113064195B (en) High-precision low-computation carrier attitude measurement method utilizing geometric characteristics of multiple antennas
CN116299599A (en) INS-assisted GNSS pseudo-range coarse difference detection method
CN115184977A (en) Integrated combined navigation device and navigation system
CN115267858A (en) Precise single-point positioning method assisted by regional navigation system
CN110646817A (en) Method for calculating positioning error and high-precision positioning method
Renga et al. Fusion of Multi-Antenna Carrier Phase Differential GPS and Inertial Measurements for Performance Evaluation of High Accuracy Integrated Aircraft Navigation Systems.
Kramlikh et al. Estimating the Inertial Characteristics of a Nanosatellite Using a Radio Compass Based on GNSS Technology

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