CN108873034A - A kind of implementation method of inertial navigation subcarrier ambiguity resolution - Google Patents

Abstract

A kind of implementation method of inertial navigation subcarrier ambiguity resolution disclosed by the invention, includes the following steps：First using Fast integer Ambiguity Resolution in the location information constraint auxiliary satellite navigation relative positioning of inertial reference calculation, dead reckoning, speed calculation are carried out using inertial navigation；Then again using its position as constraint solving integer ambiguity.The present invention for satellite navigation receiver carrier phase be accurately positioned in the Carrier Phase Ambiguity Resolution scheme deficiency that there are search spaces is big, the time is long, it is assisted using the inertial navigation information that Strapdown Inertial Navigation System provides, compressible search space, it is significantly reduced the complexity of solution, it realizes quickly approaching and solve for fuzziness, improves the speed and reliability of solution.

Description

A kind of implementation method of inertial navigation subcarrier ambiguity resolution

Technical field

The present invention relates to satellite navigation and positioning field, in particular to a kind of realization side of inertial navigation subcarrier ambiguity resolution Method.

Background technique

Currently, can be used for the method for carrier phase high-precision relative positioning ambiguity resolution mainly has：Least square search Adjustment of correlated observations method (LMABDA, Least-Squares AMBiguity Deeorrelation drops in method, least square Adjustment) etc..

One, least square search method

Least square search method is proposed by Hatch (1989,1990) earliest.The basic thought of this method is all double In poor fuzziness, only there are three be independent.As long as can determine that the integer ambiguity of three double differences, other double differences are fuzzy Degree can be with only determination.Least square search method can be divided into three steps:

1. determining the coordinate of unknown point and establishing Ambiguity Search Space.

The initial coordinate of unknown point can be used pseudorange double difference and is calculated using least square method.It is unknown acquiring After the initial coordinate of point, it can the precision of pseudo range difference solution is built as index (generally taking the three times standard deviation of each coordinate components) A three-dimensional coordinate search space is stood, is solved respectively with three basic double differences of eight apex coordinates in the space and selection Calculate corresponding fuzziness initial value.Then according to the fuzziness initial value being calculated on each vertex, these three double difference moulds are determined The respective max-int of paste degree parameterWith smallest positive integral valueThe fuzziness for needing to detect in this search space Combining sum K is：

2. least square is searched for.

Least square search for the step of be：

(1) one group of fuzziness (referred to as basic fuzziness group) to be detected is chosen from Ambiguity Search Space, utilizes phase Three double differences answered calculate dynamic point coordinate.

(2) integer ambiguity for calculating other double difference carrier phases using the dynamic point coordinate acquired is (referred to as remaining fuzzy Degree group).

(3) the double difference integer ambiguity according to obtained in (l), (2) is observed using the epoch all double difference carrier phases Value carries out least square resolving again, obtains the coordinate and corresponding residual vector V of dynamic point.

(4) variance of unit weight is calculated

V is observation residual vector in formula, and C is association's factor battle array of double difference observation, and n is double difference carrier phase observation data Number, u are the number of unknown number.

IfLess than a certain limit value, then by this group of fuzziness parameter andValue deposit destination file.Otherwise the group is obscured Parameter is spent to reject.

(5) (1) one (4) are repeated, until having detected all fuzziness combinations.

3. fixed integer ambiguity.

If certain epoch is only left one group of fuzziness parameter after carrying out step 2, then this group of fuzziness is correct fuzziness.It is no Then to saving in destination fileCarry out Ratio inspection

If Ratio is greater than a certain limit value (the general constant chosen greater than 2), then it is assumed thatMinimum corresponding fuzziness parameter Group is correct fuzziness, otherwise also needs to carry out least square to remaining fuzziness group using the data of next epoch to search Rope, until being greater than a certain limit value until being left unique one group or Ratio.

Two, adjustment of correlated observations method --- LAMBDA method drops in least square

Another carrier phase high-precision relative positioning Ambiguity Solution Methods are LAMBDA method, basic thought in the industry It is the improvement to least square adjustment.It is by the double difference observation equation after linearisation：

Y=Aa+Bb+ ε；

In formula, y is double difference carrier phase observation data vector, and a is double difference integer ambiguity vector, and b is unknown baseline component, A, B is the design matrix of fuzziness and baseline component, and ε is error vector.

The solution of above formula is to resolve least square problem

min||y-Aa-Bb||²,a∈Z^m,b∈Rⁿ；

Due to constraint a ∈ Z^m, above formula is integer least square problem, make traditional least square resolve method by To limit value.The resolving of integer least square will be carried out in two steps, the first step:With a ∈ R^mInstead of a ∈ Z^mWith ordinary least square solution The float-solution of calculation method resolving integer ambiguityTwo steps:Resolve least square problem：

It is solvedIt is the estimation of fuzziness integer least square.Above formula is resolved, first has to determine an integer vectors set, Then an integer vectors are therefrom selected, corresponding target function value is calculated, make the integer vectors of objective function minimum value It is the target to be searched for.This integer vectors set, that is, so-called search space must assure that comprising required solution, directly It connects as follows with objective function one ideal search space of above formula：

WhereinIt is variance matrix Q_aInverse matrix, it be one withCentered on polynary ellipsoid, its shape by Variance matrix Q_aIt determines, size can be by selecting χ²To change.Obviously, as long as can ensure that search space includes minimum of one net Lattice point (point that each coordinate components are all integer), it just includes required solution.In order to ensure this point, χ²Value cannot select It is too small, but can not select too big, in order to avoid include a large amount of unnecessary mesh points.

Although theoretically, search process recited above is all well and good, really can be used to resolve double difference fuzziness when It waits, efficiency is very low, and search is resulted even in when serious and is stopped.In order to avoid the generation of this situation, abandon highly relevant Double difference fuzziness, establish the lower fuzziness set of new degree of correlation.Integer Gauss transformation is used to a：

Wherein：Z^TIt is the fuzziness transformation matrix of n dimension, according to this transformation, double difference Ambiguity Search Space becomes：

In order to keep the integer characteristic of fuzziness, and original search space and transformed search space are had one by one Corresponding relationship, transformation matrix Z^TIt is necessary for integer battle array.After the transformation of fuzziness transformation matrix, fuzziness as a result will lead to Spheroid is become approximate spheres, constancy of volume by the extruding of search space.Correlation between fuzziness is reduced, and than original Double difference fuzziness it is more accurate.

That concludes LAMBDA method can be divided into following 4 steps：

1. seeking the float-solution of basic lineal vector and fuzziness using least square adjustment；

2. using integer Gauss transformation (transform), the correlation of integer ambiguity is reduced, improves Ambiguity Search Space；

3. search space after the conversion resolves fuzziness, pass through inverse transformation Z^-1It can be former the fuzziness conversion resolved Beginning fuzziness.Due to Z^-1It is only made of integer, therefore the integer characteristic of fuzziness can be kept.

4. seeking basic lineal vector using least square method with fixed solution of fuzzy degree.

In current existing scheme of the invention, more or less there are some problems in various methods.Least square is searched for Method, the selection of basic satellite be it is particularly important, resolve efficiency, same to reduce the quantity of fuzziness group substantially undetermined, improve When take into account and carry out positioning calculation using four satellites of basic satellite group and can have higher precision, to guarantee remaining double difference carrier wave phase The integer ambiguity of position can be resolved correctly, be elected to be four satellites of basic satellite group PDOP value answer it is moderate, can neither be too big It cannot be too small.In fact, it is very big for selecting different basic satellite computational efficiency difference.Thus how to select basic satellite is The matter of utmost importance that least square search method faces.If secondly there is week in the satellite in search process in basic satellite group Jump, and this cycle slip fails to detect, then most likely result in final search to fuzziness numerical value be wrong.Separately Outside, using least square search method, in addition to all fuzziness groups in first search epoch detection search space, subsequent epoch is equal It is only tested in previous epoch by the fuzziness that items are examined, computational efficiency is improved although doing so, if a certain A certain observation has biggish error in epoch, then refuses correct fuzziness group by a certain inspection item, thus Lead to search failure or acquires the integer ambiguity of mistake.

And for LAMBDA method, it is resolved using the observation of satellite navigation system merely, and prior information is had ignored Utilization, efficiency is lower, and often due to information content deficiency causes ambiguity search to fail.Their almost impossible single epoch solutions Calculate integer ambiguity.

Summary of the invention

The purpose of the present invention is to overcome the shortcomings of the existing technology and deficiency, provides a kind of inertial navigation subcarrier solution of fuzzy degree The implementation method of calculation.

The purpose of the present invention is realized by the following technical solution：

A kind of implementation method of inertial navigation subcarrier ambiguity resolution, includes the following steps：

First using Fast integer Ambiguity Resolution in the location information constraint auxiliary satellite navigation relative positioning of inertial reference calculation, make Dead reckoning, speed calculation are carried out with inertial navigation；

Then again using its position as constraint solving integer ambiguity.Mainly in GNSS short time losing lock and recapture It is effective that fuzziness is solved afterwards.

The implementation method of the inertial navigation subcarrier ambiguity resolution, specifically includes following steps：

The first step：Inertial alignment；The initial alignment refers to that each reference axis of determining inertial navigation system is sat relative to reference Mark mean to process；The alignment, process include that theta alignment (initial attitude for determining inertial navigation system) and initialization are led Speed that boat system defines, position；

Second step：Posture, speed and position are resolved using pure inertial guidance data；

Third step：Using pseudorange simultaneous error equation, Kalman filtering is carried out, obtains location error, and by location error It compensates to pure inertial navigation position；

4th step：Pure inertial navigation position is substituted into carrier phase observational equation, solves fuzziness, and integer ambiguity is returned again For carrier phase equation, accurate location error is obtained；

5th step：By accurate positional error compensation to pure inertial reference calculation value.

When solving integer ambiguity using carrier phase observational equation, the location information that pure inertial guidance data resolves has been substituted into, Substantially increase the calculation accuracy of integer ambiguity float-solution.

The integer ambiguity is directly rounded fixed fuzziness in solution process.During Carrier Phase Ambiguity Resolution, Since fuzziness float-solution precision is very high, it can be directly rounded fixed fuzziness, there is no need to scan for integer ambiguity.

Compared with prior art, the present invention having the following advantages that and beneficial effect：

1, the present invention exists for Carrier Phase Ambiguity Resolution scheme in the accurate positioning of satellite navigation receiver carrier phase and searches The deficiency that rope space is big, the time is long, is assisted, compressible search space using the inertial navigation information that Strapdown Inertial Navigation System provides, It is significantly reduced the complexity of solution, quickly approaching and solve for fuzziness is realized, improves the speed and reliability of solution.

2, integer ambiguity fixed rate of the invention is high, especially in the case where moonscope condition is bad, remains to lead to It crosses the reasonable scheme that resolves and fixes fuzziness.

3, the present invention can realize single epoch in the case where no other constraint conditions regardless of observation quality Integer ambiguity fix.

4, the present invention is without scanning for integer ambiguity, therefore efficiency of algorithm is high.

Detailed description of the invention

Fig. 1 is a kind of flow chart of the implementation method of inertial navigation subcarrier ambiguity resolution of the present invention.

Specific embodiment

Present invention will now be described in further detail with reference to the embodiments and the accompanying drawings, but embodiments of the present invention are unlimited In this.

Such as Fig. 1, a kind of implementation method of inertial navigation subcarrier ambiguity resolution is comprised the steps of：

The first step：Inertial alignment；

So-called initial alignment refers to determining the process that each reference axis of inertial navigation system is directed toward relative to reference frame.It is right Quasi- process includes theta alignment (initial attitude for determining inertial navigation system) and the speed that defines of initialization navigation system and position.I First carry out coarse alignment, then again using the result of coarse alignment as initial value, progress fine alignment.

Accelerometer and gyroscope on carrier measure each component and earth rate of specific force required for overcoming gravity Component, they use vector f respectively^bWithIt indicates.These vectors are related with gravitational vectors and earth angle rate vector respectively, working as Acceleration of gravity vector sum earth angle rate vector is denoted as g respectively in ground geographic coordinate systemⁿWithRecord system is b, navigation Referential is n, ω_ieFor rotational-angular velocity of the earth size, g is gravity acceleration magnitude, and B is latitude, then in navigation system,

It is knownWithDirection cosine matrix can be solved by reconstructing an equationSince rotational-angular velocity of the earth and acceleration of gravity are not parallel, we are accelerated using earth rate vector and gravity in this way The multiplication cross for spending vector is as follows to construct third equation：

Wherein, g^bFor the acceleration of gravity size in b system of carrier system.

Simultaneous equations (1), (2), obtain

Direction cosine matrix can be obtained：

Wherein,And f^bRespectively gyro and accelerometer output；

Substituting into above formula using the statistical average of (preset time period herein takes preceding 180s) in preset time period can To obtain coarse alignment result.

In order to improve alignment precision, reuses kalman filtering and carry out inertial navigation posture fine alignment.Using five state kalman Filter carries out fine alignment.Five states are respectively：Misalignment φ_x, φ_y, φ_z, indicate the attitude error to Barebone；δv_N,δv_E It is north orientation and east orientation speed error.Disturbance input is the measured value of gyro and Jia BiaoWith δ f^b。

The inertial system state error differential equation can be write as：

Wherein,Micro component, F for system mode are coefficient matrix, and δ x is system mode, G is coefficient matrix, u is used Guiding systems measured value.

In formula：

δ x=[φ_x φ_y φ_z δv_N δv_E]^T (6)

Wherein, φ_xFor the misalignment in the direction x, φ_yFor the misalignment in the direction y, φ_zFor the misalignment in the direction z, δ v_NFor north To velocity error, δ v_EFor east orientation speed error；For the measured value of gyro in the x direction,In y-direction for gyro Measured value,For the measured value of gyro in a z-direction,To add the measured value of table in the x direction,To add table in y Measured value on direction,To add the measured value of table in a z-direction；

In formula (8),For direction cosine matrix, expresses for convenience, matrix F is divided into four submatrixs, respectively F₁₁、F₁₂、F₂₁、F₂₂；

Wherein：

Wherein, ω_ieRotational-angular velocity of the earth, L are longitude, v_EFor east orientation speed, v_NFor north orientation speed, v_DIt is ground to speed, R_EFor chordwise curvature radius, h is height, R_NFor meridian circle radius, f_DIt is ground to acceleration, f_EFor east orientation acceleration, f_NFor north orientation Acceleration；

By earth radius R_EAnd R_NIt is considered as constant value R, height h is considered as 0.Only consider constant value gyroscopic drift and accelerometer constant value Zero bias, obtaining F matrix is：

Using horizontal velocity measurement error as observed quantity, then have：

Wherein Z is observation vector,For the horizontal velocity obtained by inertial reference calculation, v=[00]^TFor real standard speed, δ v_NFor north orientation speed error, δ v_EFor east orientation speed error, systematic observation equation is：

Z=H δ x+V_k,

Wherein δ x is system mode, V_kIndicate observation noise, H is observing matrix：

For system equation obtained as above and observational equation, carrying out discretization can be obtained：

δ x in formula_kFor k moment state, Φ_k,k-1For state-transition matrix, δ x_k-1For k-1 moment state, H_kIt is seen for the k moment Survey matrix, Γ_k-1For k-1 moment Tau matrix, W_k-1With W_kIt is all process noise, V_kTo measure noise, it is assumed to be irrelevant Zero mean Gaussian white noise sequence, and have：

Original state δ x₀With covariance matrix P₀It is known that then the more new formula of kalman filtering equations is：

Wherein, δ x_k-1For k-1 moment state, δ x_k-1To utilize k-1 status predication as a result, K_kFor kalman gain, P_k-1For The covariance matrix of k-1 state, Q_k-1For the state variable variance matrix at k-1 moment, P_k/k-1For δ x_k-1Variance matrix, R_kIt is seen for the k moment Survey noise variance matrix, Z_kFor k moment observation vector, P_kFor the covariance matrix of K state, I is unit vector.

Second step：Pure inertial reference calculation；

When do not defend lead information when, need pure inertial navigation individually resolve namely the time update.It mainly includes posture that time, which updates, It updates, speed updates and location updating.Coning compensation shakes compensation with stroke and uses Shuangzi sample backoff algorithm.

It is for strapdown inertia velocity differentials equation

Wherein,Expression of the speed in e system for carrier relative to the earth；For b system to the spin matrix of e system；f^b For specific force b system expression；g^eFor expression of the terrestrial gravitation in e system；For expression of the rotational-angular velocity of the earth in e system； ω_ieFor rotational-angular velocity of the earth.

The position differential equation is：

Wherein,For expression of the carrier positions in e system；

The posture differential equation：

Wherein σ is equivalent rotating vector,For the differential of equivalent rotating vector, σ is its mould, and ω is angular speed.Equivalent rotation The relationship for turning vector and quaternary number q is：

In order to avoid the Singular Value in calculating process, posture renewal uses Quaternion Algorithm, and uses Shuangzi sample circular cone Compensation.From t_k-1To t_kMoment, quaternary number renewal equation are：

Wherein,For t_k-1Moment b system to e system quaternary number,For the quaternary number at e system moment to t moment,For b It is quaternary number of the moment to t moment；

WhereinIt is the equivalent rotating vector from i system to b system, σ is its mould.Δθ₁、Δθ₂RespectivelyGyro output angle increment in period, T=t_k-t_k-1。

For the normalization for preventing truncated error influence quaternary number in calculating process, each quaternary number carries out four after the completion of updating First number standardization：

Wherein, q_iFor former quaternary number,For four cloud numbers after planning；

From quaternary number coordinates computed transition matrix

Wherein,For coordinate conversion matrix；

Speed, which updates to draw using Shuangzi sample, shakes backoff algorithm, in speed update cycle k-1 a to k, velocity differentials side The solution of journey is：

Wherein,For k moment carrier relative to the earth speed how the expression in e system,It is opposite for k-1 moment carrier In expression of the speed in e system of the earth,For the velocity compensation amount as caused by specific force,For brother's formula speed increment；

Wherein,For the rotational-angular velocity of the earth in e system,For bearer rate in e system,For the k-1 moment to k moment The coordinate conversion matrix of e system,For the coordinate conversion matrix of k-1 moment e system to b system,For in k-1 moment b system by than Velocity compensation amount caused by power, Δ v_kFor k moment bearer rate increment, Δ v_rotFor the rotation effect compensation term of speed, Δ v_scul For sculling compensation term；

Wherein (1) Δ θ, Δ θ (2), Δ v (1), Δ v (2) be respectively in a coning compensation period (i.e.) angle increase Amount and speed increment；Usual Ge Shi speed increment Δ v_g/corkItem variation is slower, and an order algorithm of Euler's integral can be used； Δθ_kFor k moment angle step, Δ v_kFor k moment speed increment, Δ v_rotFor the rotation effect compensation term of speed, Δ v_sculTo draw Paddle effect compensating item.

Location updating uses the method integrated to velocity vector.Since navigation calculation turnover rate is higher, ladder can be used Shape is integrated to improve integral accuracy.

Wherein,Expression of the position in e system for k moment carrier relative to the earth,It is opposite for k-1 moment carrier In expression of the position of the earth in e system, T is the period；

Third step：Using pseudorange simultaneous error equation, Kalman filtering is carried out, obtains location error, and by location error It compensates to pure inertial navigation position；

The measurement error as existing for inertial navigation system causes pure inertial reference calculation error occur.It can be filtered by kalman These errors are found out, pure inertial navigation is then added to and calculates in result, obtain more accurate result.

The error state vector of selection is：

X^e=[δ r^e δv^e δφ^e δω^b δf^b]^T,

Wherein, X^eFor error state vector in e system, δ r^eFor the location error in e system, δ v^eFor the velocity error in e system, δ φ^eFor misalignment, δ ω^bFor gyro output angle velocity error, δ f^bSpecific force error is exported for accelerometer.

Attitude error equations are：

Wherein,For the differential of attitude error, δ φ^eFor attitude error,The transition matrix of e system is tied up to for b,For Expression of the rotational-angular velocity of the earth error in b system,For earth rotation speedAntisymmetric matrix, expression formula is as follows：

Location error equation is：

Wherein,For the differential of location error, δ v^eFor the velocity error in e system；

Velocity error equation is：

Wherein,For the differential of velocity error,For b system to e system transition matrix,For b system to e system transition matrix Error, f^bFor the specific force of accelerometer output, δ f^bSpecific force error, δ g are exported for accelerometer^eIt is terrestrial gravitation error in e system In expression,For expression of the rotational-angular velocity of the earth in e system,For the velocity error in e system, matrix ψ^eFor misalignmentThe antisymmetric matrix constituted；

Wherein, δ φ^eFor misalignment,It is misalignment in x, y, the component on z coordinate axis；

Wherein,For specific force f^bIn the projection of e system.

Wherein, μ is earth universal gravitational constant, r=[x, y, z]^T。

The error equation of gyro and Jia Biao is：

Wherein,For expression of the rotational-angular velocity of the earth error in b system, B₀For the fixation zero bias of gyro, n is random Constant error.

The error for adding table also includes zero bias and random constant error, such as following formula：

Wherein, D₀For the fixation zero bias for adding table,For random constant error.The differential equation is：

In summary equation, error state vector X in available body-fixed coordinate system^eThe error state equation met：

Wherein,For the error state vector X of t moment^eDifferential,For the differential of location error,For speed mistake The differential of difference,For the differential of attitude error,For the differential of rotational-angular velocity of the earth error,For earth rotation angle speed Spend expression of the error in b systemError differential,Specific force error delta f is exported for accelerometer^bDifferential, F^eFor than Power f^bIn the projection of e system,For the velocity error in e system, δ φ^eFor attitude error,For b system to e system transition matrix, δ f^bFor Accelerometer exports specific force error,For earth rotation speedAntisymmetric matrix, δ r^eFor location error,Certainly for the earth Expression of the tarnsition velocity error in b system.

Write as matrix form：

Wherein, X^eFor the error state vector of t moment, δ r^eFor location error,For the velocity error in e system, δ φ^eFor Attitude error, δ ω^bFor the error of rotational-angular velocity of the earth, δ f^bSpecific force error is exported for accelerometer.

W (t)=[w_r w_v w_φ 0 0]^T,

Wherein, I_3×3It is the unit vector of 3*3,For earth rotation speedAntisymmetric matrix, F^eFor specific force f^bIn e system Projection,For b system to e system transition matrix, w_r、w_v、w_φThe respectively white Gaussian noise of position, speed and misalignment.

4th step：Pure inertial navigation position is substituted into carrier phase observational equation, solves fuzziness, and integer ambiguity is returned again For carrier phase equation, accurate location error is obtained；

To satellite navigation receiver r, the signal of satellite s, observing frequency i are received, observational equation is：

Wherein,For the pseudo range observed quantity of satellite s,For the carrier phase observed quantity of satellite s,It is satellite s to connecing The actual distance of receipts machine r, δ t_rFor satellite clock correction,For ionosphere delay,For tropospheric delay,For Multipath Errors, λ_iFor carrier wavelength,For carrier phase ambiguity, ε_pFor pseudorange noise, ε_φFor Carrier Phase Noise.

Behind ionosphere and the compensation of troposphere broadcast model, the initial position of receiver r is calculated by pseudorange One-Point Location x_r,0y_r,0z_r,0For：

[x_r,0 y_r,0 z_r,0]^T,

Geometric distance is between satellite s and the initial value of receiver r：

Wherein, x^s、y^s、z^sFor the coordinate of satellite s.

Its non-difference observation equation is (subscript on r and s is omitted in the case where not causing ambiguity),

Wherein,For the pseudo range observed quantity of satellite s,For the carrier phase observed quantity of satellite s,For satellite s with connect Geometric distance between the initial value of receipts machine r, δ t_rFor satellite clock correction,For the ionosphere delay of the calculation of initial value of receiver r,For the initial value tropospheric delay of receiver r,For range error,For ionosphere delay,For tropospheric delay,For Multipath Errors, λ_iFor carrier wavelength,For carrier phase ambiguity, ε_pFor pseudorange noise, ε_φIt makes an uproar for carrier phase Sound.

It can be obtained in the linearisation of receiver initial position：

Wherein, differential correctionalCarrier phase correctionδ X δ y δ z is location error, and linearisation coefficient is respectively：

For base station b, position is accurately known, and linearisation is carried out in its accurate location, location error [δ x_b δy_b δ z_b]^T=0, visual lines vector can be obtained according to above-mentioned steps；For movement station r, then the initial position of receiver is replaced The receiver location r that bar wall compensates is integrated and passed through for inertial navigation_ins,G：

Wherein, subscript G indicates that, using the position GPS as particle, subscript I is indicated using inertial navigation center as particle；It is arrived for b system E system transition matrix, Δ r^bFor coordinate of the bar wall in b system.Here bar wall compensation be by position by inertial navigation center compensation extremely GNSS antenna center.

Then the geometric distance of satellite s to receiver r are：

Wherein, x_ins,G、y_ins,G、z_ins,GThe receiver location that bar wall compensates, x are integrated and passed through for inertial navigation^s、y^s、z^s For the coordinate of satellite s.

Above-mentioned equation is changed accordingly.

Single poor observational equation is between standing：

Wherein,Single poor, l between pseudorange and the station of carrier phase_rb、m_rb、n_rbSingle poor linearisation system between station Number,

l_r m_r n_rAnd l_b m_b n_bThe linearisation coefficient of respectively receiver r and receiver b, δ x_rb、δy_rb、δz_rbFor baseline Vectorial coordinate,

Due to [δ x_b δy_b δz_b]^T=0, δ x_b δy_b δz_bFor the location error of receiver b,

[δx_rb δy_rb δz_rb]=[δ x_r δy_r δz_r]-[δx_b δy_b δz_b]=[δ x_r δy_r δz_r],

Its double difference pseudorange and carrier wave equation are：

Wherein,For the double difference for pseudorange and carrier phase,Coefficient, δ are linearized for double difference x_rb、δy_rb、δz_rbFor basic lineal vector coordinate,For ionosphere delay, tropospheric delay, Multipath Errors Double difference, ε_pFor pseudorange noise, ε_φFor Carrier Phase Noise；

WithSingle poor linearisation coefficient between satellite s and the station of satellite t.

When baseline length is shorter, ignore ionosphere delay double differenceTropospheric delay double differenceWhen baseline is longer, Ionosphere is eliminated by doing linear combination, or is estimated using ionosphere and troposphere as unknown number.N satellite is formed N-1 double-difference equation write as matrix form and be：

Wherein, Y is observation vector, δ P₁、δP₂、It is the pseudorange double difference and carrier wave double difference of 1 and 2 frequency ranges respectively.

Wherein, i is frequency, and n is number of satellites,For the double difference pseudorange difference of first satellite pair of i frequency range,For i Double difference carrier phase difference (the unit m), λ of first satellite pair of frequency range_iFor the wavelength of i frequency range, double difference pseudorange and carrier wave are seen Equation：

Wherein, I is unit battle array；Assuming that No. 1 satellite is reference star, satellite geometry matrix is：

Wherein,Be i satellite station between single poor linearisation coefficient, δ r is three-dimensional location error, δ I For the double difference ionospheric error of n-1 dimension, δ T is the double difference tropospheric error of n-1 dimension, N₁、N₂The double difference of respectively n-1 dimension is fuzzy Degree.The respectively observation noise battle array of pseudorange and carrier wave, V_PiIt is that the observation of i frequency range pseudorange is made an uproar Sound,It is the observation noise of i frequency range carrier phase.

Significantly, since bar wall compensation effect, herein δ r and before the δ r in formula^eIt is different.

Wherein, r_true,G、r_true,IRespectively indicate the true coordinate of receiver r using the position GPS as the coordinate of particle and Using inertial navigation center as the coordinate of particle, r_ins,G、r_ins,IRespectively indicate that the resolving of receiver r obtains using the position GPS as matter The coordinate of point and using inertial navigation center as the coordinate of particle, but due to simultaneously by r_true,IAnd r_ins,IThe compensation of bar wall is carried out, can be obtained It arrives：

δ r=r_true,G-r_ins,G=r_true,I-r_ins,I=δ r^e,

When solving integer ambiguity, there is following two mode：

1. using carrier phase observed quantity, fuzziness parameter is as state variable.

If it is short baseline case, ionospheric error δ I and tropospheric error δ T are available not as state variable State variable of the GNSS/INS tight integration under short baseline case is the error state variable of INS；

Wherein, x is error state variable, X^eFor the coordinate of receiver r, N₁、N₂For the fuzziness of 1,2 frequency ranges.

Observational equation is：

Wherein, Y is observation vector, and V is observation noise battle array, H and A matrix is same as above, λ_iThe wavelength of i frequency range, I be unit to Amount, X^eFor the coordinate of receiver r, δ r^eFor location error,For the velocity error in e system, δ φ^eFor attitude error, δ ω^bFor The error of rotational-angular velocity of the earth, δ f^bSpecific force error is exported for accelerometer.

If it is medium-long baselines situation, ionospheric error δ I and tropospheric error δ T are available as state variable State variable of the GNSS/INS tight integration under short baseline case be：

Wherein, x is error state variable, δ r^eFor location error,For the velocity error in e system, δ φ^eFor posture mistake Difference, δ ω^bFor the error of rotational-angular velocity of the earth, δ f^bSpecific force error is exported for accelerometer, δ I is ionosphere delay, and δ T is pair Tropospheric delay.

Observational equation is：

Wherein, μ₂For the ionosphere delay conversion coefficient of 2 frequency ranges.

2. not using carrier phase observed quantity, fuzziness parameter is not as state variable.

Corresponding to short baseline and Long baselines situation, above-mentioned state variable and observational equation become respectively：

Wherein survey noise battle array It is the observation noise of i frequency range pseudorange.

After Kalman filtering, location error δ r is obtained^eThe value of (ionosphere delay δ I, tropospheric delay δ T) simultaneously arrive by back substitution Carrier phase observational equation, i.e., by unknown number location error δ r (ionosphere delay δ I, troposphere in carrier phase observational equation Delay δ T) it is substituted for datum location error δ r^e(ionosphere delay δ I, tropospheric delay δ T), it is available：

Or (under the conditions of medium-long baselines)：

Wherein, ε_φFor Carrier Phase Noise.

Find out unknown numberAnd it is rounded and obtains integer ambiguity solution, then returned to fuzziness as datum generation State two formula, and by location error δ r^e(ionosphere delay δ I, tropospheric delay δ T) is resolved again as unknown number.It obtains Again the location error δ r resolved^eLater, final position vector solution is：

R=δ r+r_ins,G,

Wherein, r_ins,GRespectively indicate that the resolving of receiver r obtains using the position GPS as the coordinate of particle.

5th step：By accurate positional error compensation to pure inertial reference calculation value.

The above embodiment is a preferred embodiment of the present invention, but embodiments of the present invention are not by above-described embodiment Limitation, other any changes, modifications, substitutions, combinations, simplifications made without departing from the spirit and principles of the present invention, It should be equivalent substitute mode, be included within the scope of the present invention.

Claims (3)

1. a kind of implementation method of inertial navigation subcarrier ambiguity resolution, which is characterized in that include the following steps：

First using Fast integer Ambiguity Resolution in the location information constraint auxiliary satellite navigation relative positioning of inertial reference calculation, using used It leads and carries out dead reckoning, speed calculation；

Then again using its position as constraint solving integer ambiguity.

2. the implementation method of inertial navigation subcarrier ambiguity resolution according to claim 1, which is characterized in that specifically include with Lower step：

The first step：Inertial alignment；The initial alignment refers to each reference axis of determining inertial navigation system relative to reference frame The process of direction；The alignment, process include speed, the position that theta alignment is defined with initialization navigation system；

Second step：Posture, speed and position are resolved using pure inertial guidance data；

Third step：Using pseudorange simultaneous error equation, Kalman filtering is carried out, obtains location error, and by positional error compensation To pure inertial navigation position；

4th step：Pure inertial navigation position is substituted into carrier phase observational equation, solves fuzziness, and back substitution carries again by integer ambiguity Wave phase equation obtains accurate location error；

5th step：By accurate positional error compensation to pure inertial reference calculation value.

3. the implementation method of inertial navigation subcarrier ambiguity resolution according to claim 1 or claim 2, which is characterized in that described whole All fuzzinesses are directly rounded fixed fuzziness in solution process.