CN103969672B - A kind of multi-satellite system and strapdown inertial navigation system tight integration air navigation aid - Google Patents

A kind of multi-satellite system and strapdown inertial navigation system tight integration air navigation aid Download PDF

Info

Publication number
CN103969672B
CN103969672B CN201410204477.1A CN201410204477A CN103969672B CN 103969672 B CN103969672 B CN 103969672B CN 201410204477 A CN201410204477 A CN 201410204477A CN 103969672 B CN103969672 B CN 103969672B
Authority
CN
China
Prior art keywords
amp
satellite
centerdot
delta
system
Prior art date
Application number
CN201410204477.1A
Other languages
Chinese (zh)
Other versions
CN103969672A (en
Inventor
张涛
宫淑萍
徐晓苏
王立辉
李瑶
李佩娟
Original Assignee
东南大学
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 东南大学 filed Critical 东南大学
Priority to CN201410204477.1A priority Critical patent/CN103969672B/en
Publication of CN103969672A publication Critical patent/CN103969672A/en
Application granted granted Critical
Publication of CN103969672B publication Critical patent/CN103969672B/en

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39Determining 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/42Determining position
    • G01S19/48Determining position by combining or switching between position solutions derived from the satellite radio beacon positioning system and position solutions derived from a further system
    • G01S19/49Determining position by combining or switching between position solutions derived from the satellite radio beacon positioning system and position solutions derived from a further system whereby the further system is an inertial position system, e.g. loosely-coupled
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00Navigation; Navigational instruments not provided for in preceding groups G01C1/00-G01C19/00
    • G01C21/10Navigation; Navigational instruments not provided for in preceding groups G01C1/00-G01C19/00 by using measurements of speed or acceleration
    • G01C21/12Navigation; Navigational instruments not provided for in preceding groups G01C1/00-G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16Navigation; Navigational instruments not provided for in preceding 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/165Navigation; Navigational instruments not provided for in preceding 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

Abstract

The invention discloses a kind of multi-satellite system and strapdown inertial navigation system tight integration air navigation aid, first IMU output data strapdown is resolved and obtain SINS current pose, speed, positional information.Satellite ephemeris message is carried out ephemeris resolving and obtains each satellite position, speed, each satellite system space-time uniformity, carry out satellite and select star.Secondly, SINS position, velocity information that satellite position, velocity information and the strapdown resolving that utilization is chosen obtains are calculated SINS pseudorange, pseudorange rates.Again, utilize the raw measurement data that GNSS exports, i.e. survey code pseudorange and Doppler frequency shift, compare with SINS pseudorange, pseudorange rates, using difference as the observed quantity of wave filter, the optimal estimation of device after filtering, provide the compensation dosage of correction SINS, closed-loop corrected, obtain SINS attitude, speed, position optimal solution.The present invention solves the very few problem affecting tight integration navigation accuracy of number of satellite under the complex environments such as inertial navigation is blocked at high building, trees are covered with single combinations of satellites navigation.

Description

A kind of multi-satellite system and strapdown inertial navigation system tight integration air navigation aid

Technical field

The invention belongs to field of navigation technology, be specifically related to a kind of multi-satellite system and strapdown inertial navigation system tight integration Air navigation aid, is particularly well-suited to the combination solution in the case of number of satellite is very few under the complex environments such as high building blocks, trees are covered Calculate.

Background technology

Inertial navigation based on newton Classical Motion mechanics law, be a kind of utilize inertial measurement component measure carrier Angular movement and line kinematic parameter, on the basis of original state, carry out recursion to determine the autonomous of carrier attitude, speed and position Formula calculates air navigation aid.At present, inertial navigation technology has developed more ripe, but at navigation long-range, long-time and arm discharge Can't fully meet requirement when navigating Deng high accuracy, its major defect is that navigation positioning error accumulates in time.Navigation error Mainly determined by the dynamic characteristic of the precision being initially directed at, the error of inertial sensor and carrier movement track, the most solely After vertical work, error can dissipate quickly.The solution appearing as this problem of integrated navigation technology provides a kind of effective way, Satellite navigation can be real-time mensuration bearer rate and position, combine with inertial navigation system, by information fusion skill Art can further improve the precision of navigation system.At present, integrated navigation technology become airmanship development important directions it One.But in the case of the complex environment such as block, number is covered at high building under, number of satellite is very few, single satellite is used to strapdown Property navigation system integrated navigation calculation accuracy decline.

Summary of the invention

The very few shadow of number of satellite under complex environment such as block, trees are covered at high building for inertial navigation and single combinations of satellites navigation The problem ringing tight integration navigation accuracy, the present invention proposes a kind of method that multi-satellite system fits navigation with inertial navigation system.

A kind of multi-satellite system and strapdown inertial navigation system tight integration air navigation aid, comprise the following steps:

(1) to inertial measurement component (IMU) output data by strapdown resolve obtain inertial navigation current pose, speed, Position, carries out satellite data process to the ephemeris message of satellite systems receiver output, it is thus achieved that the speed of satellite, positional information;

(2) satellite velocities, position is utilized to be calculated by SINS pseudorange, pseudorange rates with SINS position, velocity information SINS pseudorange, pseudorange rates;

(3) the pseudorange message from satellite systems receiver output obtains survey code pseudorange and Doppler frequency shift, with SINS pseudorange, Pseudorange rates compares, and using difference as the observed quantity of wave filter, sets up tight integration state model and measurement model, through optimum Estimate and filtering, the amount of error correction of output system;Being corrected systematic error, output SINS attitude, speed, position are Excellent combination solves.

Described strapdown algorithm is angular movement and the line kinematic parameter utilizing IMU to measure carrier, at original state base Carry out recursion on plinth and determine that the autonomous type of the carrier attitude, speed and the position that are mounted with SINS infers navigation.

Described satellite data processes and includes that satellite selects star, satellite ephemeris to resolve and satellite system data merges.

It is by the position of SINS that described SINS pseudorange, pseudorange rates calculate process, the position calculation of satellite obtain corresponding to Pseudorange at SINS position;By the speed of SINS, the speed of satellite, it is thus achieved that SINS is relative to satellite relative motion i.e. SINS Pseudorange rate of change.

Described optimal estimation and filtering use Kalman filter.

Described tight integration state model and the method for building up of measurement model are as follows:

(1) multisystem tight integration state equation is set up

Strapdown error equation and satellite system are modeled and transfers state equation description to

X · I ( t ) X · GPS ( t ) X · GLO ( t ) X · Gal ( t ) X · BD ( t ) = F I ( t ) 0 0 0 0 0 F GPS ( t ) 0 0 0 0 0 F GLO ( t ) 0 0 0 0 0 F Gal ( t ) 0 0 0 0 0 F BD ( t ) X I ( t ) X GPS ( t ) X GLO ( t ) X Gal ( t ) X BD ( t ) + G I ( t ) 0 0 0 0 0 G GPS ( t ) 0 0 0 0 0 G GLO ( t ) 0 0 0 0 0 G Gal ( t ) 0 0 0 0 0 F BD ( t ) W I ( t ) W GPS ( t ) W GLO ( t ) W Gal ( t ) W BD ( t )

Wherein, X ∈ R23It is state variable, by 15 quantity of states of inertial navigation system (SINS), US Global location system System (GPS), GLONASS satellite system (GLONASS), GALILEO satellite system (Galileo), Beidou satellite system (BD) each 2 Individual quantity of state forms;

X = [ δV E δV N δV U φ E φ N φ U δL δλ δh ▿ bx ▿ by ▿ bz ϵ bx ϵ by ϵ bz δt GPS δt GLO δt Gal δt BD δt rGPS δt rGLO δt rGal δt rBD ]

In formula, δ VE,δVN,δVUIt is the velocity error northeastward on direction, three, sky,Three misalignments of strapdown Angle (platform stance angle), δ L, δ λ, δ h is that three site errors of strapdown are described by terrestrial coordinate system,It is to add table Three axial biased errors, εbxbybzIt is three axial drifts of gyro, δ tGPS、δtGLO、δtGal、δtBDBe respectively GPS, The distance of GLONASS, Galileo, BD Satellite clock errors equivalence, δ trGPS、δtrGLO、δtrGal、δtrBDBe respectively GPS, The range rate of GLONASS, Galileo, BD satellite clock frequency error equivalence, totally 23 dimension;

FIT () is inertial navigation system error state equation state matrix, and FGPS(t)、FGLO(t)、FGal(t)、FBD(t) difference It is the state matrix obtained after two quantity of states being correlated with GPS, GLONASS, Galileo, BD satellite model, GI(t) and GGPS (t)、GGLO(t)、GGal(t)、GBDT () is the input matrix of noise.δt*Represent that after modeling, it is actuated to white noise, δ tr*It is modeled as First-order Markov process, * represents each system (GPS, GLONASS, Galileo, BD)

F * = 0 1 0 - β tr * , G * = 1 0 0 1 , W * ( t ) = ω t * ω tr * ;

(2) multisystem tight integration measurement equation is set up

(a) pseudo range measurement establishing equation

δt*It is station clock error, represents station clock clock correction equivalent distances error, vρ*jRepresent other common error and exclusive error, Such as measurement error such as ionospheric error, tropospheric error and Multipath Errors;In tight integration emulates, vρ*jIt is considered one in vain The error of form of noise, the most multiple passage of clock correction error is unified to be arranged;In tight integration prototype experiment, then ionosphere must be prolonged Late etc. all multiple errors compensate, and the method for compensation is discussed below;

Taking satellite, being write as matrix form has

δρ GPS δρ GLO δρ Gal δ ρ BD = H GPS 1 GPS 0 0 0 H GLO 0 1 GLO 0 0 H Gal 0 0 1 Gal 0 H BD 0 0 0 1 BD δp δt GPS δt GLO δt Gal δt BD + υ ρGPS υ ρGLO υ ρGal υ ρBD

Wherein, δ p=[δ x δ y δ z]TFor User Status amount, H*For geometry observing matrix;Its dimension and selected system-satellite Number (n) is relevant is 2;H*=[e1 e2 e3]2*3, ekTie up for 2*1;

The measurement equation of pseudo range observed quantity

Zρ(t)=Hρ(t)X(t)+Vρ(t)

Wherein Hρ(t)=[0N*3 0N*3 Hρ1 0N*6 Hρ2 0N*3]N*23, N is 8, i.e. the sum of satellite;

Hρ1=(aji)N*3, H ρ 2 = H ρGPS H ρGLO H ρGal H ρBD

And

H ρGPS = 1 0 0 0 1 0 0 0 , H ρGLO = 1 0 0 0 1 0 0 0 , H ρGal = 1 0 0 0 1 0 0 0 , H ρBD = 1 0 0 0 1 0 0 0

aj1=(Rn+h)[-ej1sinLcosλ-ej2sinLsinλ]+[Rn*(1-e2)+h]ej3cosL

aj2=(Rn+h)[ej2cosLcosλ-ej1cosLsinλ]

aj3=ej1cosLcosλ+ej2cosLsinλ+ej3sinL

B () pseudorange rates measurement equation is set up

The measurement equation of pseudorange rates observed quantity

Z ρ · ( t ) = H ρ · ( t ) X ( t ) + V ρ · ( t )

In formula, H ρ · ( t ) = H ρ · 1 0 N * 3 H ρ · 2 0 N * 9 H ρ · 3

Wherein, H ρ · 1 = ( b ji ) N * 3 , H ρ · 2 = ( c ji ) N * 3 , H ρ · 3 = 1 GPS 1 GLO 1 Gal 1 BD

bj1=-e*j1sinλ+e*j2cosλ

bj2=-e*j2sinLcosλ-e*j1sinLsinλ+e*j3cosL

bj3=e*j1cosLcosλ+e*j2cosLsinλ+e*j3sinL

c j 1 = ( R n + h ) [ - e · * j 1 sin L cos λ - e · * j 2 sin L sin λ ] + [ R n * ( 1 - e 2 ) + h ] e · * j 3 cos L

c j 2 = ( R n + h ) [ e · * j 2 cos L cos λ - e · * j 1 cos L sin λ ]

c j 3 = e · * j 1 cos L cos λ + e · * j 2 cos L sin λ + e · * j 3 sin L

Combine (a), (b) measurement equation:

Z ( t ) = H ρ H ρ · X ( t ) + V ρ ( t ) V ρ · ( t )

Described satellite selects star i.e. to select suitable satellite to carry out navigation calculation, and detailed process is: first, it is determined that each satellite System-satellite number;If certain Satellite System satellites number is less than two, then this system is got rid of;Secondly, according to independent satellite With SINS integrated navigation, choose 4 star combinations;Combine with SINS according to double star base system, increase the synchronous error of satellite system Unknown number, at least needs 5 stars to carry out navigation calculation i.e. 3+2 pattern;6 are needed when combining with SINS according to three satellite systems Star carries out positioning calculation i.e. 2+2+2 pattern;8 stars are needed to carry out positioning calculation i.e. 2+ during according to four satellite system combinations 2+2+2 pattern;Finally, choosing on the premise of population of satellite mesh and each satellite system number meet according to maximum extremity multiaspect The selecting-star algorithm of body volumetric method or method of least square carries out selecting star.

Described satellite ephemeris computation refers to judge the form of the ephemeris message that satellite system exports, and defends according to difference The form of star system ephemeris message carries out ephemeris parsing;Satellite system ephemeris is divided into two kinds according to form: the broadcast ephemeris of class GPS The radio news program form of parameter format and class GLONASS;The packet parsing of the radio news program form of class GPS uses and defends Star position algorithm;The packet parsing of the radio news program form of class GLONASS uses derivation equation of motion integral operation method.

Described satellite system data fusion specially selects main satellite system, and other secondary satellite system is according to main satellite The reference time of system, coordinate system carry out the conversion of time, coordinate.

When described 3+2 pattern i.e. double star seat combines with SINS, a satellite system chooses 3 another system of satellite choosings Take 2 satellites;Described 2+2+2 pattern refers to that three satellite systems choose two satellites respectively;Described 2+2+2+2 pattern is four Individual satellite system chooses two satellites respectively.

Compared with prior art, present invention have the advantage that

(1) the single satellite system before extension and the work of inertial navigation integrated navigation, increase multiple satellite system to group Close the effect of navigation system, it is ensured that the precision of integrated navigation under the complex environment that number of satellite is few, improve integrated navigation system Reliability.

(2) present invention introduces 4 kinds of satellites in the whole world is and the research of SINS tight integration that in the long run, this grinds The achievement in research studied carefully can be transplanted on the combination research of any n global satellite system and strap-down inertial system easily, right The research of the application of satellite/inertia navigation system combination has great meaning.

Accompanying drawing explanation

Fig. 1 is principle of the invention block diagram;

Fig. 2 is that satellite data of the present invention processes content;

Fig. 3 is integrated navigation main program flow chart of the present invention;

Fig. 4 is ephemeris process of analysis figure of the present invention;

Fig. 5 is selecting-star algorithm flow chart of the present invention.

Detailed description of the invention

Below in conjunction with the accompanying drawings, it is further elucidated with the present invention.

As it is shown in figure 1, the present invention uses tight integration method to complete integrated navigation, in order to improve system under complex environment Precision, have employed the process of multi-satellite system data and multi-satellite system and the method for SINS integrated navigation, implements in invention Step is as follows:

(1) to inertial measurement component (IMU) output data by strapdown resolve obtain inertial navigation current pose, speed, Position.

SINS attitude angle calculates:

Using Quaternion Method to calculate attitude matrix, according to Euler's theorem, moving coordinate system is relative to the orientation etc. of reference frame Imitate and rotate an angle, θ in moving coordinate system around certain Equivalent Axis, if representing the unit vector in Equivalent Axis direction with u, then The orientation of moving coordinate system is determined by two parameters of u and θ completely.A quaternary number can be constructed with u and θ:

Q = cos θ 2 + u sin θ 2

Above formula derivation abbreviation can be obtained quaternion differential equation:

Q ( q · ) = 1 2 M * ( ω b ) Q ( q )

In formula M * ( ω b ) = 0 - ω nb bx - ω nb by - ω nb bz ω nb bx 0 ω nb bz - ω nb by ω nb by - ω nb bz 0 ω nb bx ω nb bz ω nb by - ω nb bx 0

Solve quaternion differential equation according to complete card approximatioss to obtain:

q ( t ) = { cos Δ θ 0 2 I + sin Δ θ 0 2 Δ θ 0 [ Δθ ] } q ( 0 )

In formula

Δ θ 0 = Δ θ x 2 + Δ θ y 2 + Δ θ z 2

[ Δθ ] = ∫ t 1 t 1 + h M * ( ω nb b ) dt = 0 - Δ θ x - Δ θ y - Δ θ z Δ θ x 0 Δ θ z - Δ θ y Δ θ y - Δ θ z 0 Δ θ x Δ θ z Δ θ y - Δ θ x 0

In formula

Δ θ i = ∫ t t + h ω nb bi dt , i = x , y , z .

The spin velocity making terrestrial coordinate system relative inertness coordinate system is ωie, (its value is 15.04088 °/h), L represents Local latitude, λ represents local longitude, then

ωie n: the spin velocity of terrestrial coordinate system relative inertness coordinate system vector in geographic coordinate system, for:

ω ie n = 0 ω ie cos L ω ie sin L T

ωie b: the spin velocity of terrestrial coordinate system relative inertness coordinate system vector in carrier coordinate system, for:

ω ie b = C n b ω ie n

Attitude matrix in formula, when carrier stationary, is determined by initial angle;When carrier rotates relative to geographic coordinate system, Attitude matrix and then change, quaternary number try to achieve (lower same) after immediately revising.

ωen n: geographical coordinate relatively spherical coordinate system rotational angular velocity vector in geographic coordinate system, for:

ω en n = - V N / R N V E / R E V E tan L / R E T

VE、VNIt is respectively east orientation and the north orientation speed of carrier movement;

RNFor the radius of curvature in reference ellipsoid meridian plane, RN=Re(1-2e+3esin2L);

REFor the radius of curvature in the plane normal of vertical meridian plane, RE=Re(1+esin2L);

Wherein ReMajor axis radius for reference ellipsoid;E is the ovality of ellipsoid.

Again because, L · = V N / R N , λ · = V E / ( R E cos L ) Then

ω en n = - L · λ · cos L λ · sin L T

ωen b: geographical coordinate relatively spherical coordinate system rotational angular velocity vector in carrier coordinate system, for:

ω en b = C n b ω en n

ωib b: gyro output angle speed, it is designated as

ω ib b = ω ib bx ω ib by ω ib bz T

ωnb b: carrier coordinate system, relative to vector in carrier coordinate system of the rotational angular velocity of geographic coordinate system, is designated as

ω nb b = ω nb bx ω nb by ω nb bz T

Then can obtain

ωnb bib bie ben b

After quaternary number is revised immediately, can be by first real-time update attitude matrix of quaternary number according to following formula

C n b = C 11 C 12 C 13 C 21 C 22 C 23 C 31 C 32 C 33 = q 0 2 + q 1 2 - q 2 2 - q 3 2 2 ( q 1 q 2 + q 0 q 3 ) 2 ( q 1 q 3 - q 0 q 2 ) 2 ( q 1 q 2 - q 0 q 3 ) q 0 2 - q 1 2 + q 2 2 - q 3 2 2 ( q 2 q 3 + q 0 q 1 ) 2 ( q 1 q 3 + q 0 q 2 ) 2 ( q 2 q 3 - q 0 q 1 ) q 0 2 - q 1 2 - q 2 2 + q 3 2

Real-time attitude angle can be extracted from attitude battle array

SINS speed calculation:

Ratio force vector in carrier coordinate system is fb, then geographic coordinate system has:

f n = C b n f b

Direction cosine matrix in formulaWhen carrier stationary, initial angle determine;When carrier is relative to geographic coordinate system During rotation, direction cosine matrixAnd then change, quaternary number try to achieve after immediately revising.

Had by inertial navigation fundamental equation:

V · n = f n - ( 2 ω ie n + ω en n ) × V n + g n

Write as component form to have:

V · E V · N V · U = f E f N f U + 0 ( λ · + 2 ω ie ) sin L - ( λ · + 2 ω ie ) cos L - ( λ · + 2 ω ie ) sin L 0 - L · ( λ · + 2 ω ie ) cos L L · 0 V E V N V U + 0 0 - g

In formula: fnThe projection fastened at navigation coordinate for carrier acceleration, fn=[fE fN fU]T;VnRepresent that hull is being led Velocity in boat coordinate system, Vn=[VE VN VU]T;gnFor gravity acceleration, gn=[0 0-g]T

Integration above formula, can try to achieve each velocity component V that carrier is fastened at navigation coordinateE、VN、VU

SINS position calculation:

The differential equation of longitude and latitude can be expressed as follows:

L · = V N R N + h λ · = V E ( R E + h ) cos L

In formula, h is height.

The more new formula of the longitude and latitude of integration above formula can obtain longitude and latitude:

L = ∫ L · dt + L ( 0 ) λ = ∫ λ · dt + λ ( 0 )

(2) the ephemeris message of satellite systems receiver output is carried out satellite data process, it is thus achieved that the speed of satellite, position Information;

Classify according to table 1 ephemeris format, it can be seen that satellite system ephemeris is divided into two kinds according to form: the broadcast of class GPS The radio news program form of ephemeris parameter form and class GLONASS.

Table 1 ephemeris format is classified

The parsing of the radio news program form ephemeris of class GPS:

A () makes

tk=t-toe

Wherein toeThe reference moment for satellite almanac data;tkIt it is normalized relative time.

B the mean angular velocity n of () satellite is

N=n0+Δn

WhereinIt it is the theoretical mean angular velocity of satellite;Δ n is the correction term of satellite mean angular velocity;μ= 3.986005×1014It it is geocentric gravitational constant.

C mean anomaly M in () observation moment is

M=M0+n*tk

Wherein M0For the mean anomaly with reference to the moment.

D () can the eccentric anomaly in calculating observation moment by mean anomaly

E=M+es*sinE

For above formula, eccentricity e of satellite orbitsThe most little (less than 0.02), according to solution by iterative method, must be complete Office's convergence, makes the initial value E of E0=M, typically can have Δ E/E < 10 after iteration ten times-10Precision.

E () calculates the earth's core of satellite to footpath r according to eccentric anomaly

R=a* (1-es*cosE)

WhereinSquare root for elliptical orbit major semiaxis.

F () calculates the true anomaly f of satellite according to eccentric anomaly

tan ( f / 2 ) = ( 1 + e s / 1 - e s * tan ( E / 2 ) )

Value during calculating when above formula is near E=π (180 °) can be very big, can use an easy method, i.e. Threshold value is set to E or tan (E/2).Or utilize following formula to calculate

sin f = 1 - e s 2 · sin E 1 - e s cos E cos f = cos E - e s 1 - e s cos E

G () calculates ascending node angular distance φ according to true anomaly and obtains

φ=f+ ω

Wherein ω is the argument of perigee of satellite orbit.

H () is revised satellite orbit parameter according to the perturbation correction term in ephemeris and ascending node angular distance and is calculated perturbation correction member

Ascending node angular distance correction value δ μ

δ μ=Cuc*cos(2φ)+Cus*sin(2φ)

The earth's core is to footpath correction value δ r

δ r=Crc*cos(2φ)+Crs*sin(2φ)

Orbit inclination angle correction value δ i

δ i=Cic*cos(2φ)+Cis*sin(2φ)

Wherein Cuc,CusMediation correction term for ascending node angular distance;Crc,CrsFor satellite the earth's core to the mediation correction term in footpath; Cic,CisMediation correction term for inclination of satellite orbit.

I () updates the orbit parameter of Current observation moment satellite

Ascending node angular distance φk φk=φ+δ μ

The earth's core is to footpath rk rk=r+ δ r

Orbit inclination angle ik ik=i0+idot*tk+δi

Wherein i0For with reference to moment satellite orbit plane relative to the inclination angle of earth equatorial plane;Idot is satellite orbit plane The rate of change at inclination angle.

J () calculates at normalization moment tkRight ascension of ascending node Ωk

Ω k = Ω e + Ω · * t k - ω ie * t

Wherein, ωieFor earth rotation angular speed, ωie=7.2921151467 × 10-5(rad/s);ΩeFor satellite orbit Right ascension of ascending node;Rate of change for right ascension of ascending node;

(k) coordinate of the satellite position

In the plane at satellite orbit place, setting up initial point is the earth's core, and Z axis orientation north is just, X-axis takes elliptic orbit major axis Coordinate system, obtaining coordinate of the satellite position is

P = r k * cos ( φ k ) r k * sin ( φ k ) 0

The similar process that l process that () satellite velocities calculates calculates with coordinate, it should be noted that should utilize star clock by mistake Difference derivative term carrys out calibration satellite speed with the derivative term of general-relativistic effect.

Satellite place coordinate system Satellite position coordinates and speed can be calculated according to above-mentioned steps.

The parsing of the radio news program form ephemeris of class GLONASS:

A () sets up equation of satellite motion

Assuming to there is an inertial coodinate system (0-xyz), zero is consistent with co-ordinates of satellite system, and coordinate axes points to and te Moment co-ordinates of satellite system coordinate axes points to identical.If any instant t, satellite position (X (t), Y (t), Z in body-fixed coordinate system (t)), in inertial system 0-xyz, corresponding position is (x (t), y (t), z (t)), then the relation between them is:

X ( t ) Y ( t ) Z ( t ) = cos ω ( t - t e ) sin ω ( t - t e ) 0 - sin ω ( t - t e ) cos ω ( t - t e ) 0 0 0 1 x ( t ) y ( t ) z ( t ) = R ( t ) x ( t ) y ( t ) z ( t )

Then its single order, second-order differential relational expression are:

X · ( t ) Y · ( t ) Z · ( t ) = R ( t ) x · ( t ) y · ( t ) z · ( t ) + R · ( t ) x ( t ) y ( t ) z ( t )

X · · ( t ) Y · · ( t ) Z · · ( t ) = R ( t ) x · · ( t ) y · · ( t ) z · · ( t ) + 2 R · ( t ) x · ( t ) y · ( t ) z · ( t ) + R · · ( t ) x ( t ) y ( t ) z ( t )

Known R ( t ) = cos ω ( t - t e ) sin ω ( t - t e ) 0 - sin ω ( t - t e ) cos ω ( t - t e ) 0 0 0 1 , Then

R · ( t ) = - ω sin ω ( t - t e ) ω cos ω ( t - t e ) 0 - ω cos ω ( t - t e ) ω sin ω ( t - t e ) 0 0 0 1 = ω 0 1 0 - 1 0 0 0 0 0 R ( t )

In like manner

R · · ( t ) = ω 2 1 0 0 0 1 0 0 0 0 R ( t )

To sum up can obtain:

X · · ( t ) Y · · ( t ) Z · · ( t ) = R ( t ) x · · ( t ) y · · ( t ) z · · ( t ) + 2 ω 0 1 0 - 1 0 0 0 0 0 X · ( t ) Y · ( t ) Z · ( t ) + ω 2 1 0 0 0 1 0 0 0 0 X ( t ) Y ( t ) Z ( t )

It is

X · · ( t ) Y · · ( t ) Z · · ( t ) = R ( t ) x · · ( t ) y · · ( t ) z · · ( t ) + ω 2 X ( t ) + 2 ω Y · ( t ) ω 2 Y ( t ) - 2 ω X · ( t ) 0

On the right of above formula, Section 1 is the expression in body-fixed coordinate system of the satellite motion acceleration.In theory, it is necessary to meter Calculate satellite by acceleration produced by whole active forces (earth center gravitation and whole perturbative force).Owing to the time of integration is shorter, Earth nonspherical gravitation perturbation (as long as typically considering the impact of C item), the Sun and the Moon Gravitational perturbation need only be calculated, just can reach Required precision.

In body-fixed coordinate system, if the quality of the earth is Me, GMe=3.9860044 × 1014m3/s2For Gravitational coefficient of the Earth, The acceleration that then earth center gravitation causes is:

a e = - GM e r 3 X ( t ) Y ( t ) Z ( t )

r = X 2 ( t ) + Y 2 ( t ) + Z 2 ( t )

The acceleration that the earth non-soccer star gravitation causes is:

a n = 1.5 GM e r 5 R e 2 C 2 0 X ( t ) ( r 2 - 5 Z ( t ) 2 ) / r 2 Y ( t ) ( r 2 - 5 Z ( t ) 2 ) / r 2 Z ( t ) ( 3 r 2 - 5 Z ( t ) 2 / r 2 )

Wherein, ReReference ellipsoid major radius (Re=6378136m),For the gravitational field coefficient of unnormalized, for normalization CoefficientTake advantage of 5 ( C 2 0 = 5 C - 2 0 = 1082625.7 × 10 - 9 ) .

The radio news program form ephemeris navigation message of class GLONASS provides teSolar-lunar perturbate with reference to the moment and cause Acceleration, the numerical value in body-fixed coordinate system isWithin a short period of time regards constant as.

To sum up, satellite equation of satellite motion in body-fixed coordinate system is:

X · · ( t ) Y · · ( t ) Z · · ( t ) = - GM e r 3 x · · ( t ) y · · ( t ) z · · ( t ) + 1.5 GM e r 5 R e 2 C 2 0 X ( t ) ( r 2 - 5 Z ( t ) 2 ) / r 2 Y ( t ) ( r 2 - 5 Z ( t ) 2 ) / r 2 Z ( t ) ( 3 r 2 - 5 Z ( t ) 2 ) / r 2 + X · · ( t e ) Y · · ( t e ) Z · · ( t e ) + ω 2 X ( t ) + 2 ω Y · ( t ) ω 2 Y ( t ) - 2 ω X · ( t ) 0

Satellite position during (b) calculating any time t and speed

According to the satellite equation of motion in body-fixed coordinate system, functional form can be able to and express tiMoment acceleration function is public Formula is as follows:

X · · i = f 1 ( t i , X i , Y i , Z i , X · i , Y · i , Z · i , X · · ( t e ) ) Y · · i = f 2 ( t i , X i , Y i , Z i , X · i , Y · i , Z · i , Y · · ( t e ) ) Z · · i = f 3 ( t i , X i , Y i , Z i , X · i , Y · i , Z · i , Z · · ( t e ) )

With t0Moment as original state, then tiThe integral equation of moment satellite position is:

r · ( t i ) = r · 0 + ∫ t 0 t i r · · dt

r ( t i ) = r ( t 0 ) + ∫ t 0 t i r · dt

r(t0),For the coordinate of the satellite position vector sum speed with reference to the moment, calculate t by twice Integral SolutioniMoment Satellite position, uses fixed step size quadravalence Long Beige-Ku Ta to be integrated track, and its integral formula is as follows:

r n + 1 = r n + 1 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 )

First acceleration is integrated, calculates a step-length rate integrating result, carry out speed further according to rate integrating result Degree integration, coordinates computed integral result, the speed then obtained according to the integral and calculating of this step, position quantity carry out next step Long integration, successively cycle calculations, so circulation is until tiMoment terminates.The number of times of integration and the step-length of acquirement and the time of integration Relevant computing formula:Finally obtain co-ordinates of satellite and the speed of t.

(3) satellite velocities, position is utilized to be calculated by SINS pseudorange, pseudorange rates with SINS position, velocity information SINS pseudorange, pseudorange rates;

SINS pseudorange ρIjCalculate:

Making carrier positions coordinate under ECEF coordinate system is (xI,yI,zI), strapdown be given;Jth satellite is at ginseng coordinate Position in system is (xsj,ysj,zsj), satellite ephemeris try to achieve.The carrier positions that is given of inertial navigation to the pseudorange of jth satellite For

ρ Ij = ( x I - x sj ) 2 + ( y I - y sj ) 2 + ( z I - z sj ) 2

SINS pseudorange ratesCalculate:

ρIjDerivativeFor:

ρ · Ij = d ( ( x I - x sj ) 2 + ( y I - y sj ) 2 + ( z I - z sj ) 2 ) / dt = ( x I - x sj ) ( x · I - x · sj ) + ( y I - y sj ) ( y · I - y · sj ) + ( z I - z sj ) ( z · I - z · sj ) ρ Ij = e j 1 ( x · I - x · sj ) + e j 2 ( y · I - y · sj ) + e j 3 ( z · I - z · sj )

Wherein,It is SINS speed under ECEF coordinate system, (xsj,ysj,zsj) it is that jth satellite is in reference Satellite velocities under co-ordinates of satellite system.And ( x I - x sj ) / ρ Ij = e j 1 ( y I - y sj ) / ρ Ij = e j 2 ( z I - z sj ) / ρ Ij = e j 3 .

(4) the pseudorange message from satellite systems receiver output obtains survey code pseudorange and Doppler frequency shift, with SINS pseudorange, Pseudorange rates compares, and using difference as the observed quantity of wave filter, sets up tight integration state model and measurement model.Through optimum Estimate and filtering, the amount of error correction of output system.Being corrected systematic error, output SINS attitude, speed, position are Excellent combination solves.

Strapdown inertial navigation system error equation and the foundation of model:

(a) velocity error equation

Velocity error equation general type is

Launch to obtain by formula

δ · V E = f N φ U - f U φ N + ( V N R E tan L - V U R E ) δ V E + ( 2 ω ie sin L + V E R E tan L ) δ V N - ( 2 ω ie cos L + V E R E ) δ V U + ( 2 ω ie ( V N cos L + V U sin L ) + V E V N R E sec 2 L ) δL + V E V U - V E V N tan L R E 2 δh + ▿ E δ · V N = f U φ E - f E φ U - 2 ( ω ie sin L + V E R E tan L ) δ V E - V U R N δ V N - V N R N δ V U - ( 2 ω ie cos L + V E R E sec 2 L ) V E δL + ( V N V U R N 2 + V E 2 tan L R E 2 ) δh + ▿ N δ · V U = f E φ N - f N φ E + 2 ( ω ie cos L + V E R E ) δV E + 2 V N R N δV N - 2 ω ie sin LV E δL - ( V N 2 R N 2 + V E 2 R E 2 ) δh + ▿ U

(b) attitude error equations

The attitude error equations of inertial navigation is

Formula is launched

(c) site error equation

By L · = V N R N , λ · = V E R E cos L Can obtain

δ L · = δV N R N δ λ · = δ V E R E sec L + V E R E sec L tan LδL δ h · = δ V U

Satellite system models:

The error state of satellite, in pseudorange, pseudorange rates combined system, generally takes two errors relevant with the time: one Individual be and clocking error equivalence range error δ tu, another is and the range error δ t of clock frequency error equivalenceru

A () GPS system models

The observed quantity of SINS/GPS tight integration have selected pseudorange and pseudorange rates, according to pseudorange, the mathematical model GPS of pseudorange rates Part needs to choose two relevant quantity of states, the clock correction equivalent distances error of clock of i.e. standing and clock frequency difference equivalent rate error.Its Model is described as

δ t · GPS = δ t rGPS + ω tGPS δ t · rGPS = - β trGPS δ t rGPS + ω trGPS

Above formula shows, in GPS, clock jitter is the differential equation of first order of white-noise excitation, and frequency deviation of clock is single order horse Markov process, βtrGPSFor the inverse correlation time.

(b) GLONASS system modelling

Mathematical model GLONASS part according to pseudorange, pseudorange rates needs to choose two relevant quantity of states, i.e. station clock Clock correction equivalent distances error and clock frequency difference equivalent rate error.Its model is described as

δ t · GLO = δ t rGLO + ω tGLO δ t · rGLO = - β trGLO δ t rGLO + ω trGLO

(c) Galileo system modelling

Mathematical model Galileo part according to pseudorange, pseudorange rates needs to choose two relevant quantity of states, i.e. station clock Clock correction equivalent distances error and clock frequency difference equivalent rate error.Its model is described as

δ t · Gal = δ t rGal + ω tGal δ t · rGal = - β trGal δ t rGal + ω trGal

(d) BD system modelling

Mathematical model BD part according to pseudorange, pseudorange rates needs to choose two relevant quantity of states, the clock correction of clock of i.e. standing Equivalent distances error and clock frequency difference equivalent rate error.Its model is described as

δ t · BD = δ t rBD + ω tBD δ t · rBD = - β trBD δ t rBD 1 + ω trBD

Multisystem tight integration state equation is set up:

Strapdown error equation and satellite system are modeled and transfers state equation description to

X · I ( t ) X · GPS ( t ) X · GLO ( t ) X · Gal ( t ) X · BD ( t ) = F I ( t ) 0 0 0 0 0 F GPS ( t ) 0 0 0 0 0 F GLO ( t ) 0 0 0 0 0 F Gal ( t ) 0 0 0 0 0 F BD ( t ) X I ( t ) X GPS ( t ) X GLO ( t ) X Gal ( t ) X BD ( t ) + G I ( t ) 0 0 0 0 0 G GPS ( t ) 0 0 0 0 0 G GLO ( t ) 0 0 0 0 0 G Gal ( t ) 0 0 0 0 0 F BD ( t ) W I ( t ) W GPS ( t ) W GLO ( t ) W Gal ( t ) W BD ( t )

Wherein, X ∈ R23It is state variable, by the quantity of state of 15 inertial navigations and the quantity of state of 2 GPS, 2 GLONASS Quantity of state, 2 Galileo, the quantity of state composition of 2 BD.

X = [ δV E δV N δV U φ E φ N φ U δL δλ δh ▿ bx ▿ by ▿ bz ϵ bx ϵ by ϵ bz δt GPS δt GLO δt Gal δt BD δt rGPS δt rGLO δt rGal δt rBD ] State Measure from the 1st to 21 dimension respectively:

δVE,δVN,δVUIt is the velocity error northeastward on direction, three, sky,Three misalignments of strapdown are (flat Platform attitude angle), δ L, δ λ, δ h is that three site errors of strapdown are described by terrestrial coordinate system,It is to add three axles of table To biased error, εbxbybzIt is three axial drifts of gyro, δ tGPS、δtGLO、δtGal、δtBDBe respectively GPS, The distance of GLONASS, Galileo, BD Satellite clock errors equivalence, δ trGPS、δtrGLO、δtrGal、δtrBDBe respectively GPS, The range rate of GLONASS, Galileo, BD satellite clock frequency error equivalence, totally 23 dimension.

In formula, FIT () is inertial navigation system error state equation state matrix, and FGPS(t)、FGLO(t)、FGal(t)、FBD(t) It is the state matrix obtained after two quantity of states being correlated with GPS, GLONASS, Galileo, BD satellite model respectively, GI(t) And GGPS(t)、GGLO(t)、GGal(t)、GBDT () is the input matrix of noise.δt*Represent that after modeling, it is actuated to white noise, δ tr* Being modeled as first-order Markov process, * represents each system (GPS, GLONASS, Galileo, BD).As follows

F * = 0 1 0 - β tr * , G * = 1 0 0 1 , W * ( t ) = ω t * ω tr * ;

Multisystem tight integration measurement equation is set up:

(a) pseudo range measurement establishing equation

δt*It is station clock error, represents station clock clock correction equivalent distances error, vρ*jRepresent other common error and exclusive error, Such as measurement error such as ionospheric error, tropospheric error and Multipath Errors.In tight integration emulates, vρ*jIt is considered one in vain The error of form of noise, the most multiple passage of clock correction error is unified to be arranged.In tight integration prototype experiment, then ionosphere must be prolonged Late etc. all multiple errors compensate, and the method for compensation is discussed below.

Taking satellite, being write as matrix form has

δρ GPS δρ GLO δρ Gal δ ρ BD = H GPS 1 GPS 0 0 0 H GLO 0 1 GLO 0 0 H Gal 0 0 1 Gal 0 H BD 0 0 0 1 BD δp δt GPS δt GLO δt Gal δt BD + υ ρGPS υ ρGLO υ ρGal υ ρBD

Wherein, δ p=[δ x δ y δ z]TFor User Status amount, H*For geometry observing matrix.Its dimension and selected system-satellite Number (n) is relevant is 2.H*=[e1 e2 e3]2*3, ekTie up for 2*1.

The measurement equation of pseudo range observed quantity

Zρ(t)=Hρ(t)X(t)+Vρ(t)

Wherein Hρ(t)=[0N*3 0N*3 Hρ1 0N*6 Hρ2 0N*3]N*23, N is 8, i.e. the sum of satellite.

Hρ1=(aji)N*3, H ρ 2 = H ρGPS H ρGLO H ρGal H ρBD

And

H ρGPS = 1 0 0 0 1 0 0 0 , H ρGLO = 1 0 0 0 1 0 0 0 , H ρGal = 1 0 0 0 1 0 0 0 , H ρBD = 1 0 0 0 1 0 0 0

aj1=(Rn+h)[-ej1sinLcosλ-ej2sinLsinλ]+[Rn*(1-e2)+h]ej3cosL

aj2=(Rn+h)[ej2cosLcosλ-ej1cosLsinλ]

aj3=ej1cosLcosλ+ej2cosLsinλ+ej3sinL

B () pseudorange rates measurement equation is set up

The measurement equation of pseudorange rates observed quantity

Z ρ · ( t ) = H ρ · ( t ) X ( t ) + V ρ · ( t )

In formula, H ρ · ( t ) = H ρ · 1 0 N * 3 H ρ · 2 0 N * 9 H ρ · 3

Wherein, H ρ · 1 = ( b ji ) N * 3 , H ρ · 2 = ( c ji ) N * 3 , H ρ · 3 = 1 GPS 1 GLO 1 Gal 1 BD

bj1=-e*j1sinλ+e*j2cosλ

bj2=-e*j2sinLcosλ-e*j1sinLsinλ+e*j3cosL

bj3=e*j1cosLcosλ+e*j2cosLsinλ+e*j3sinL

c j 1 = ( R n + h ) [ - e · * j 1 sin L cos λ - e · * j 2 sin L sin λ ] + [ R n * ( 1 - e 2 ) + h ] e · * j 3 cos L

c j 2 = ( R n + h ) [ e · * j 2 cos L cos λ - e · * j 1 cos L sin λ ]

c j 3 = e · * j 1 cos L cos λ + e · * j 2 cos L sin λ + e · * j 3 sin L

Combine (a), (b) measurement equation:

Z ( t ) = H ρ H ρ · X ( t ) + V ρ ( t ) V ρ · ( t )

(5) compensation terminates, and exports optimum combination result.Return continues to start to perform step (1).

Claims (9)

1. a multi-satellite system and strapdown inertial navigation system tight integration air navigation aid, it is characterised in that comprise the following steps:
(1) inertial measurement component (IMU) output data are resolved acquisition inertial navigation current pose, speed, position by strapdown, The ephemeris message of satellite systems receiver output is carried out satellite data process, it is thus achieved that the speed of satellite, positional information;
(2) satellite velocities, position is utilized to be calculated SINS puppet with SINS position, velocity information by SINS pseudorange, pseudorange rates Away from, pseudorange rates;
(3) the pseudorange message from satellite systems receiver output obtains survey code pseudorange and Doppler frequency shift, with SINS pseudorange, pseudorange Rate compares, and using difference as the observed quantity of wave filter, sets up tight integration state model and measurement model, through optimal estimation With filtering, the amount of error correction of output system;Systematic error is corrected, exports SINS attitude, speed, the optimal set of position Close and solve;Described tight integration state model and the method for building up of measurement model are as follows:
(3.1) multisystem tight integration state equation is set up
Strapdown error equation and satellite system are modeled and transfers state equation description to
X · I ( t ) X · G P S ( t ) X · G L O ( t ) X · G a l ( t ) X · B D ( t ) = F I ( t ) 0 0 0 0 0 F G P S ( t ) 0 0 0 0 0 F G L O ( t ) 0 0 0 0 0 F G a l ( t ) 0 0 0 0 0 F B D ( t ) X I ( t ) X G P S ( t ) X G L O ( t ) X G a l ( t ) X B D ( t ) + G I ( t ) 0 0 0 0 0 G G P S ( t ) 0 0 0 0 0 G G L O ( t ) 0 0 0 0 0 G G a l ( t ) 0 0 0 0 0 F B D ( t ) W I ( t ) W G P S ( t ) W G L O ( t ) W G a l ( t ) W B D ( t )
Wherein, X ∈ R23It is state variable, the 15 of SINS quantity of states, american global positioning system (GPS), GLONASS defends Star system (GLONASS), GALILEO satellite system (Galileo), Beidou satellite system (BD) each 2 quantity of states composition;
X = δV E δV N δV U φ E φ N φ U δ L δ λ δ h ▿ b x ▿ b y ▿ b z ϵ b x ϵ b y ϵ b z δt G P S δt G L O δt G a l δt B D δt r G P S δt r G L O δt r G a l δt r B D
In formula, δ VE,δVN,δVUIt is the velocity error northeastward on direction, three, sky,Three platform stances of strapdown Angle, δ L, δ λ, δ h is that three site errors of strapdown are described by terrestrial coordinate system,It is to add axial inclined of three, table Put error, εbxbybzIt is three axial drifts of gyro, δ tGPS、δtGLO、δtGal、δtBDBe respectively GPS, GLONASS, The distance of Galileo, BD Satellite clock errors equivalence, δ trGPS、δtrGLO、δtrGal、δtrBDBe respectively GPS, GLONASS, The range rate of Galileo, BD satellite clock frequency error equivalence, totally 23 dimension;FIT () is inertial navigation system error state equation State matrix, and FGPS(t)、FGLO(t)、FGal(t)、FBDT GPS, GLONASS, Galileo, BD satellite is correlated with by () respectively The state matrix obtained after two quantity of state modelings, GI(t) and GGPS(t)、GGLO(t)、GGal(t)、GBDT () is the input of noise Matrix;δt*Represent that after modeling, it is actuated to white noise,Be modeled as first-order Markov process, * represent each system GPS, GLONASS、Galileo、BD
F * = 0 1 0 - β t r * , G * = 1 0 0 1 , W * ( t ) = ω t * ω t r * ;
(3.2) multisystem tight integration measurement equation is set up
(a) pseudo range measurement establishing equation
δt*It is station clock error, represents station clock clock correction equivalent distances error, vρ*jRepresent other common error and exclusive error, such as electricity The measurement error such as absciss layer error, tropospheric error and Multipath Errors;In tight integration emulates, vρ*jIt is considered a white noise The error of form, the most multiple passage of clock correction error is unified to be arranged;In tight integration prototype experiment, then must be to ionosphere delay etc. All multiple errors compensate, and the method for compensation is discussed below;
Taking satellite, being write as matrix form has
δρ G P S δρ G L O δρ G a l δρ B D = H G P S 1 G P S 0 0 0 H G L O 0 1 G L O 0 0 H G a l 0 0 1 G a l 0 H B D 0 0 0 1 B D δ p δt G P S δt G L O δt G a l δt B D + υ ρ G P S υ ρ G L O υ ρ G a l υ ρ B D
Wherein, δ p=[δ x δ y δ z]TFor User Status amount, H*For geometry observing matrix;Its dimension and selected system-satellite number N () is relevant is 2;H*=[e1 e2 e3]2*3, ekTie up for 2*1;The measurement equation of pseudo range observed quantity
Zρ(t)=Hρ(t)X(t)+Vρ(t)
Wherein Hρ(t)=[0N*3 0N*3 Hρ1 0N*6 Hρ2 0N*3]N*23, N is 8, i.e. the sum of satellite;
Hρ=(aji)N*3,
And
H ρ G P S = 1 0 0 0 1 0 0 0 , H ρ G L O = 1 0 0 0 1 0 0 0 , H ρ G a l = 1 0 0 0 1 0 0 0 , H ρ B D = 1 0 0 0 1 0 0 0
aj1=(Rn+h)[-ej1sinLcosλ-ej2SinLsin λ]+[Rn*(1-e2)+h] ej3cosL
aj2=(Rn+ h) [ej2cosLcosλ-ej1cosLsinλ]
aj3=ej1cosLcosλ+ej2cosLsinλ+ej3sinL
B () pseudorange rates measurement equation is set up
The measurement equation of pseudorange rates observed quantity
Z ρ · ( t ) = H ρ · ( t ) X ( t ) + V ρ · ( t )
In formula,
Wherein,
bj1=-e*j1sinλ+e*j2cosλ
bj2=-e*j2sinLcosλ-e*j1sinLsinλ+e*j3cosL
bj3=e*j1cosLcosλ+e*j2cosLsinλ+e*j3sinL
e j 1 = ( R n + h ) [ - e · * j 1 sin L c o s λ - e · * j 2 sin L s i n λ ] + [ R n * ( 1 - e 2 ) + h ] e · * j 3 cos L
c j 2 = ( R n + h ) [ e · * j 2 cos L c o s λ - e · * j 1 cos L s i n λ ]
c j 3 = e · * j 1 cos L c o s λ + e · * j 2 cos L s i n λ + e · * j 3 sin L
Combine (a), (b) measurement equation:
Z ( t ) = H ρ H ρ · X ( t ) + V ρ ( t ) V ρ · ( t ) .
Multi-satellite system the most according to claim 1 and strapdown inertial navigation system tight integration air navigation aid, is characterized in that: Described strapdown resolves angular movement and the line kinematic parameter being to utilize IMU to measure carrier, carries out recursion on the basis of original state Determine that the autonomous type of the carrier attitude, speed and the position that are mounted with SINS infers navigation.
Multi-satellite system the most according to claim 1 and strapdown inertial navigation system tight integration air navigation aid, is characterized in that: Described satellite data processes and includes that satellite selects star, satellite ephemeris to resolve and satellite system data merges.
A kind of multi-satellite system the most according to claim 1 and strapdown inertial navigation system tight integration air navigation aid, it is special Levy and be: described SINS pseudorange, pseudorange rates calculating process are by the position of SINS, and the position calculation of satellite obtains corresponding to SINS Pseudorange at position;By the speed of SINS, the speed of satellite, it is thus achieved that SINS is relative to the satellite relative motion i.e. puppet of SINS Away from rate of change.
A kind of multi-satellite system the most according to claim 1 and strapdown inertial navigation system tight integration air navigation aid, it is special Levy and be: described optimal estimation and filtering use Kalman filter.
Multi-satellite system the most according to claim 3 and strapdown inertial navigation system tight integration air navigation aid, is characterized in that: Described satellite selects star i.e. to select suitable satellite to carry out navigation calculation, and detailed process is: first, it is determined that each Satellite System satellites Number;If certain Satellite System satellites number is less than two, then this system is got rid of;Secondly, according to independent satellite and SINS group Close navigation, choose 4 star combinations;Combine with SINS according to double star base system, increase the synchronous error unknown number of satellite system, extremely 5 stars are needed to carry out navigation calculation i.e. 3+2 pattern less;6 stars are needed to carry out when combining with SINS according to three satellite systems Positioning calculation i.e. 2+2+2 pattern;8 stars are needed to carry out positioning calculation i.e. 2+2+2+2 mould during according to four satellite system combinations Formula;Finally, choosing on the premise of population of satellite mesh and each satellite system number meet according to maximum extremity polyhedron volume The selecting-star algorithm of method or method of least square carries out selecting star.
A kind of multi-satellite system the most according to claim 3 and strapdown inertial navigation system tight integration air navigation aid, it is special Levy and be: described satellite ephemeris computation refers to judge the form of the ephemeris message that satellite system exports, according to different satellites The form of System almanac message carries out ephemeris parsing;Satellite system ephemeris is divided into two kinds according to form: the broadcast ephemeris ginseng of class GPS The radio news program form of number format and class GLONASS;The packet parsing of the radio news program form of class GPS uses satellite Position algorithm;The packet parsing of the radio news program form of class GLONASS uses derivation equation of motion integral operation method.
A kind of multi-satellite system the most according to claim 3 and strapdown inertial navigation system tight integration air navigation aid, it is special Levy and be: described satellite system data fusion specially selects main satellite system, and other secondary satellite system is according to main satellite system The reference time of system, coordinate system carry out the conversion of time, coordinate.
A kind of multi-satellite system the most according to claim 6 and strapdown inertial navigation system tight integration air navigation aid, it is special Levy and be: when described 3+2 pattern i.e. double star seat combines with SINS, a satellite system is chosen 3 another systems of satellite and chosen 2 Satellite;Described 2+2+2 pattern refers to that three satellite systems choose two satellites respectively;Described 2+2+2+2 pattern is four Satellite system chooses two satellites respectively.
CN201410204477.1A 2014-05-14 2014-05-14 A kind of multi-satellite system and strapdown inertial navigation system tight integration air navigation aid CN103969672B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201410204477.1A CN103969672B (en) 2014-05-14 2014-05-14 A kind of multi-satellite system and strapdown inertial navigation system tight integration air navigation aid

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201410204477.1A CN103969672B (en) 2014-05-14 2014-05-14 A kind of multi-satellite system and strapdown inertial navigation system tight integration air navigation aid

Publications (2)

Publication Number Publication Date
CN103969672A CN103969672A (en) 2014-08-06
CN103969672B true CN103969672B (en) 2016-11-02

Family

ID=51239400

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201410204477.1A CN103969672B (en) 2014-05-14 2014-05-14 A kind of multi-satellite system and strapdown inertial navigation system tight integration air navigation aid

Country Status (1)

Country Link
CN (1) CN103969672B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106767787A (en) * 2016-12-29 2017-05-31 北京时代民芯科技有限公司 A kind of close coupling GNSS/INS combined navigation devices

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104729530B (en) * 2014-10-13 2017-06-06 北京航空航天大学 A kind of Hardware In The Loop Simulation Method of inertial navigation/Big Dipper integrated navigation
CN104596514B (en) * 2015-01-12 2017-08-29 东南大学 The Real-time Noisy Reducer and method of accelerometer and gyroscope
CN105093249A (en) * 2015-08-12 2015-11-25 浙大正呈科技有限公司 Inertial navigation device
CN105182375B (en) * 2015-09-29 2017-08-25 中国电子科技集团公司第五十四研究所 The pseudolite systems receiver carrier wave tracing method aided in based on inertial navigation system
CN105424041A (en) * 2016-01-18 2016-03-23 重庆邮电大学 Pedestrian positioning algorithm based on BD/INS (Beidou/Inertial Navigation System) tight coupling
CN105737823B (en) * 2016-02-01 2018-09-21 东南大学 A kind of GPS/SINS/CNS Combinated navigation methods based on five rank CKF
CN105954783B (en) * 2016-04-26 2017-03-29 武汉大学 A kind of real-time tight integration of improvement GNSS/INS navigates the method for real-time performance
CN106595705B (en) * 2016-11-22 2019-12-20 北京航天自动控制研究所 Method for estimating initial reference deviation of inertia in flight based on GPS
CN109085619A (en) * 2017-06-14 2018-12-25 展讯通信(上海)有限公司 Localization method and device, storage medium, the receiver of multimode GNSS system
CN108469627A (en) * 2018-03-16 2018-08-31 中国电子科技集团公司第三十六研究所 Based on when frequency difference ground with frequency more stationary radiant sources localization method and system
CN108802805B (en) * 2018-06-07 2019-10-11 武汉大学 Six degree of freedom broadband integration earthquake monitoring system
CN108919370B (en) * 2018-07-25 2019-11-29 赛德雷特(珠海)航天科技有限公司 A kind of positioning device and method based on gravitation field measurement

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106767787A (en) * 2016-12-29 2017-05-31 北京时代民芯科技有限公司 A kind of close coupling GNSS/INS combined navigation devices

Also Published As

Publication number Publication date
CN103969672A (en) 2014-08-06

Similar Documents

Publication Publication Date Title
Grewal et al. Global navigation satellite systems, inertial navigation, and integration
Noureldin et al. Fundamentals of inertial navigation, satellite-based positioning and their integration
Hofmann-Wellenhof et al. Global positioning system: theory and practice
Jäggi et al. Phase center modeling for LEO GPS receiver antennas and its impact on precise orbit determination
Bisnath et al. Current state of precise point positioning and future prospects and limitations
Han et al. Integrated GPS/INS navigation system with dual-rate Kalman Filter
Rabbou et al. Tightly coupled integration of GPS precise point positioning and MEMS-based inertial systems
Hofmann-Wellenhof et al. GNSS–global navigation satellite systems: GPS, GLONASS, Galileo, and more
Haines et al. One-centimeter orbit determination for Jason-1: new GPS-based strategies
Hasan et al. A review of navigation systems (integration and algorithms)
US8232917B2 (en) Post-mission high accuracy position and orientation system
DE102011076602A1 (en) GNSS - Atmospheric estimation with a Federated-Ionospheric filter
CN101788296B (en) SINS/CNS deep integrated navigation system and realization method thereof
CN103575299B (en) Utilize dual-axis rotation inertial navigation system alignment and the error correcting method of External Observation information
Montenbruck et al. Precision real-time navigation of LEO satellites using global positioning system measurements
CN101609140B (en) Compatible navigation receiver positioning system and positioning method thereof
Schumacher Integration of a gps aided strapdown inertial navigation system for land vehicles
Wang et al. Quadratic extended Kalman filter approach for GPS/INS integration
EP2570823B1 (en) Method and apparatus for differential global positioning system (DGPS) - based real time attitude determination (RTAD)
CN101514900B (en) Method for initial alignment of a single-axis rotation strap-down inertial navigation system (SINS)
US20090093959A1 (en) Real-time high accuracy position and orientation system
CN103675861B (en) Satellite autonomous orbit determination method based on satellite-borne GNSS multiple antennas
CN103090866B (en) Method for restraining speed errors of single-shaft rotation optical fiber gyro strapdown inertial navigation system
CN106289246A (en) A kind of rods arm measure method based on position and orientation measurement system
Mulder et al. Non-linear aircraft flight path reconstruction review and new advances

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant