CN106996778A - Error parameter scaling method and device - Google Patents

Abstract

The invention discloses a kind of error parameter scaling method and device.This method includes：Set up comprising error parameter vector in the state equation and measurement equation of aircraft guidance system, the state equation and measurement equation, the error parameter vector is made up of multiple error parameters；Differentiate the observability of each error parameter；When there is observable error parameter, using preset duration as filtering cycle, using Kalman filtering or adaptive-filtering, the observable error parameter is calibrated.The present invention realize in real time, on-orbit calibration go out Guidance instrumentation error parameter purpose.

Description

Error parameter scaling method and device

Technical field

The present invention relates to aircraft navigation control technology, more particularly to a kind of error parameter scaling method and device.

Background technology

Inertia type instrument is the heart of SINS, and the size of its error will directly affect the essence of entering the orbit of spacecraft The size of degree and offset landings.Some spacecrafts for performing survey of deep space, after certain altitude is flown to, it is impossible to reuse group Navigation is closed to be modified；Or important strike weapon, behind the native country that flies out, in order to improve its antijamming capability, Neng Gouzheng Combinations of satellites navigation feature can also be closed by really performing strike task.Now, if it is possible to carry out an inertia system Guidance instrumentation Error coefficient on-orbit calibration, will be an important means of raising spacecraft navigation accuracy, and error coefficient mask work Core be method for parameter estimation research.

At present, the method comparative maturity of inertia system error parameter ground demarcation, but inertia system practical application During into aerial mission, due to being shaken by aircraft, space environment is changed, and the error parameter of ground demarcation can not meet winged Row device, spacecraft accurately determine the requirement of appearance positioning in space flight for a long time.

The content of the invention

Present invention solves the technical problem that being：Compared to prior art, there is provided a kind of error parameter scaling method and dress Put, realize the purpose that real-time, on-orbit calibration goes out the error parameter of Guidance instrumentation.

The above-mentioned purpose of the present invention is achieved by the following technical programs：

In a first aspect, the invention provides a kind of error parameter scaling method, including：

Set up in the state equation and measurement equation of aircraft guidance system, the state equation and measurement equation and include Error parameter vector, the error parameter vector is made up of multiple error parameters；

Differentiate the observability of each error parameter；

When there is observable error parameter, using preset duration as filtering cycle, Kalman filtering or adaptive is utilized Filtering, calibrates the observable error parameter.

Further, the error parameter includes the misaligned angle of the platform error, velocity error, site error, the scale of gyro Factor error, gyroscopic drift error, the drift error of the scale factor error of accelerometer and accelerometer；It is each described Error parameter is vector, and includes three durection components.

Further, the state equation is：

In formula (1), A represents state matrix；The transition matrix of launching inertial system is tied to for aircraft body；X represents to miss Poor parameter vector,φ represents the misaligned angle of the platform error, δ V Velocity error is represented, δ r represent site error, δ K_gRepresent the scale factor error of gyro, b₁Represent gyroscopic drift error, δ K_aTable Show the scale factor error of accelerometer,The drift error of accelerometer is represented, T represents transposition computing；η_gRepresent that gyro is surveyed The white noise of amount, ε_aFor the white noise of accelerometer measures；For error parameter vector X first derivative vector；The state Matrix A is：

In formula (2), The aircraft body system measured for gyro is relative to transmitting inertia The angular speed of system aircraft body system projection,Represent withFor the matrix of the elements in a main diagonal；f^bThe specific force measured for accelerometer aircraft body system projection,Represent SolveAntisymmetric matrix；I_3×3The unit matrix of expression 3 × 3；G Newtonian gravitational constant is represented, M represents earth quality, and x, y, z represents coordinate of the aircraft under launching inertial system, and r represents flight Distance of the device to launching inertial system origin；The transition matrix of launching inertial system is tied to for aircraft body；

The measurement equation is：

In formula (3), Z (t) represents guidance matrix of differences；φ^i″The attitude angle and inertial reference calculation posture measured for star sensor The difference at angle；Z_v(t) it is the GPS aircraft speeds measured and the difference of the speed of inertial reference calculation；Z_r(t) the aircraft position measured for GPS Put the difference with the position of inertial reference calculation；H_9×21Represent measurement matrix；X represents the error parameter vector；V_9×1Represent white noise to Amount；The measurement matrix H_9×21For：

In formula (4), I_3×3The unit matrix of expression 3 × 3；The white noise vector V_9×1For：

In formula (5),WithWhite noise, δ M are measured for the posture of star sensor_x、δM_yWith δ M_zFor GPS amount Degree of testing the speed white noise, δ x_G、δx_GWith δ z_GFor GPS adjustment location white noise；T represents transposition computing.

Further, the observability of each error parameter is differentiated, including：

Using the state matrix, transfer matrix is calculated, the calculation formula of the transfer matrix is：

In formula (6), Φ_k,k-1The expression state k-1 moment is to the transfer matrix at state k moment, I_21×21The list of expression 21 × 21 Bit matrix；T₁Represent preset duration；A_kRepresent the state matrix of k-th of preset duration；M is positive integer and m >=2；

According to the transfer matrix and the measurement matrix, singular value decomposition matrix is treated in calculating, described to treat singular value decomposition The calculation formula of matrix is：

In formula (7), Q represents to treat singular value decomposition matrix, H_9×21Represent measurement matrix, Φ_1,0The moment of expression state 0 is to shape The transfer matrix at the moment of state 1, Φ_m-1,m-2The transfer matrix at expression state m-2 moment to state m-1 moment, T represents that transposition is transported Calculate；

Treat that singular value decomposition matrix carries out singular value decomposition by described, to obtain the first singular value vector, the second singular value Vector and singular value matrix, the formula of the singular value decomposition is：

Q=U ∑s V^T (8)

In formula (8), Q represents to treat singular value decomposition matrix, and ∑ represents singular value matrix, and the elements in a main diagonal of ∑ is σ_i,i ∈[1,m]；U represents the first singular value vector, U=[u_i]=[u₁u₂...u_m]；V represents the second singular value vector, V=[v_i]= [v₁ v₂ ... v_m]；T represents transposition computing；

Utilize the σ_i、u_i、v_iAnd the guidance matrix of differences Z (t), calculate observability discriminant vector Y_21×1, institute The calculation formula for stating observability discriminant vector Y is：

In formula (9), T represents transposition computing；

By the observability discriminant vector Y_21×1Each element be compared with given threshold δ e；

As the observability discriminant vector Y_21×1L, the element at the column position of l ∈ [1,21] row the 1st is more than setting During threshold value δ e, l in decision errors parameter vector X, the error parameter Observable at the column position of l ∈ [1,21] row the 1st.

Second aspect, present invention also offers a kind of error parameter caliberating device, the device includes：

Build module, state equation and measurement equation for setting up aircraft guidance system, the state equation and amount Survey comprising error parameter vector in equation, the error parameter vector is made up of multiple error parameters；

Discrimination module, the observability for differentiating each error parameter；

Demarcating module, for when there is observable error parameter, using preset duration as filtering cycle, utilizing Kalman Filtering or adaptive-filtering, calibrate the observable error parameter.

Further, the error parameter includes the misaligned angle of the platform error, velocity error, site error, the scale of gyro Factor error, gyroscopic drift error, the drift error of the scale factor error of accelerometer and accelerometer；It is each described Error parameter is vector, and includes three durection components.

Further, the state equation is：

In formula (1), A represents state matrix；C_b ⁱThe transition matrix of launching inertial system is tied to for aircraft body；X represents to miss Poor parameter vector,φ represents the misaligned angle of the platform error, δ V Velocity error is represented, δ r represent site error, δ K_gRepresent the scale factor error of gyro, b₁Represent gyroscopic drift error, δ K_aTable Show the scale factor error of accelerometer,The drift error of accelerometer is represented, T represents transposition computing；η_gRepresent that gyro is surveyed The white noise of amount, ε_aFor the white noise of accelerometer measures；For error parameter vector X first derivative vector；The state square Battle array A be：

In formula (2), The aircraft body system measured for gyro is relative to launching inertial system Angular speed aircraft body system projection,Represent withFor the matrix of the elements in a main diagonal；f^bThe specific force measured for accelerometer aircraft body system projection,Represent SolveAntisymmetric matrix；I_3×3The unit matrix of expression 3 × 3；G Newtonian gravitational constant is represented, M represents earth quality, x,_y, z represents coordinate of the aircraft under launching inertial system, and r represents flight Distance of the device to launching inertial system origin；The transition matrix of launching inertial system is tied to for aircraft body；

The measurement equation is：

In formula (3), Z (t) represents guidance matrix of differences；φ^i″The attitude angle and inertial reference calculation posture measured for star sensor The difference at angle；Z_v(t) it is the GPS aircraft speeds measured and the difference of the speed of inertial reference calculation；Z_r(t) the aircraft position measured for GPS Put the difference with the position of inertial reference calculation；H_9×21Represent measurement matrix；X represents the error parameter vector；V_9×1Represent white noise to Amount；The measurement matrix H_9×21For：

In formula (4), I_3×3The unit matrix of expression 3 × 3；The white noise vector V_9×1For：

In formula (5),WithWhite noise, δ M are measured for the posture of star sensor_x、δM_yWith δ M_zFor GPS amount Degree of testing the speed white noise, δ x_G、δx_GWith δ z_GFor GPS adjustment location white noise；T represents transposition computing.

Further, the discrimination module includes：

First computing unit, for utilizing the state matrix, calculates transfer matrix, the calculation formula of the transfer matrix For：

In formula (6), Φ_k,k-1The expression state k-1 moment is to the transfer matrix at state k moment, I_21×21The list of expression 21 × 21 Bit matrix；T₁Represent preset duration；A_kRepresent the state matrix of k-th of preset duration；M is positive integer and m >=2；

Second computing unit, for according to the transfer matrix and the measurement matrix, calculating and treating singular value decomposition matrix, It is described treat singular value decomposition matrix calculation formula be：

In formula (7), Q represents to treat singular value decomposition matrix, H_9×21Represent measurement matrix, Φ_1,0The moment of expression state 0 is to shape The transfer matrix at the moment of state 1, Φ_m-1,m-2The transfer matrix at expression state m-2 moment to state m-1 moment, T represents that transposition is transported Calculate；

Resolving cell, for by it is described treat singular value decomposition matrix carry out singular value decomposition, with obtain the first singular value to Amount, the second singular value vector and singular value matrix, the formula of the singular value decomposition is：

Q=U ∑s V^T (8)

In formula (8), Q represents to treat singular value decomposition matrix, and ∑ represents singular value matrix, and the elements in a main diagonal of ∑ is σ_i,i ∈[1,m]；U represents the first singular value vector, U=[u_i]=[u₁ u₂ ... u_m]；V represents the second singular value vector, V=[v_i] =[v₁ v₂ ... v_m]；T represents transposition computing；

3rd computing unit, for utilizing the σ_i、u_i、v_iAnd the guidance matrix of differences Z (t), calculate Observable Property discriminant vector Y_21×1, the calculation formula of the observability discriminant vector Y is：

In formula (9), T represents transposition computing；

Comparing unit, for by the observability discriminant vector Y_21×1Each element compared with given threshold δ e Compared with；

Identifying unit, for as the observability discriminant vector Y_21×1L, at the column position of l ∈ [1,21] row the 1st When element is more than given threshold δ e, l in decision errors parameter vector X, the error at the column position of l ∈ [1,21] row the 1st Parameter Observable.

The present invention has the advantages that compared with prior art：

(1), the state equation and measurement equation of the invention by setting up aircraft guidance system；Differentiate each error parameter Observability；When there is observable error parameter, using Kalman filtering or adaptive-filtering, calibrate described considerable The error parameter of survey, can overcome the shortcomings of existing guidance instrument error isolation technics " world is inconsistent ", can be real-time, in-orbit The error coefficient of Guidance instrumentation is calibrated, algorithm is simple, is easy to engineering.

(2), the present invention can reduce inertia device Support expense, due to the ability with on-orbit calibration, so as to The number of times of inertia system ground periodic calibrating is enough reduced, is used manpower and material resources sparingly.

(3), the present invention can improve aircraft guidance precision, reduce the dependence to satellite navigation, extending space aircraft The scope and ability of execution task.

Brief description of the drawings

Fig. 1 is a kind of flow chart of error parameter scaling method in the embodiment of the present invention one；

Fig. 2 is a kind of structure chart of error parameter caliberating device in the embodiment of the present invention two.

Embodiment

The present invention is described in further detail with reference to the accompanying drawings and examples.It is understood that described herein Specific embodiment be used only for explaining the present invention, rather than limitation of the invention.It also should be noted that, for the ease of Description, part related to the present invention rather than entire infrastructure are illustrate only in accompanying drawing.

Embodiment one

Fig. 1 is a kind of flow chart of error parameter scaling method in the embodiment of the present invention one, and the present embodiment is applicable to The situation to the error parameter progress on-orbit calibration of Guidance instrumentation is needed, this method can be held by error parameter caliberating device OK, wherein the device can be realized by software and/or hardware.With reference to Fig. 1, the error parameter scaling method tool that the present embodiment is provided Body may include steps of：

In S110, the state equation and measurement equation for setting up aircraft guidance system, the state equation and measurement equation Comprising error parameter vector, the error parameter vector is made up of multiple error parameters.

Specifically, the error parameter include the misaligned angle of the platform error, velocity error, site error, the scale of gyro because Count the drift error of error, gyroscopic drift error, the scale factor error of accelerometer and accelerometer；Each mistake Poor parameter is vector, and includes three durection components.Each error parameter is vector, and includes three directions Component, i.e., each error parameter is vector, in respective coordinate system comprising three direction vectors.For example, gyroscopic drift error Three durection components are included in its coordinate system：X-axis is to component, y-axis to component and z-axis to component.In another example, accelerometer Scale factor error include three durection components in its coordinate system：X-axis is to component, y-axis to component and z-axis to component.

Specifically, the state equation is：

In formula (1), A represents state matrix；The transition matrix of launching inertial system is tied to for aircraft body；X represents to miss Poor parameter vector,φ represents the misaligned angle of the platform error, δ V represents velocity error, and δ r represent site error, δ K_gRepresent the scale factor error of gyro, b₁Represent gyroscopic drift error, δ K_a The scale factor error of accelerometer is represented,The drift error of accelerometer is represented, T represents transposition computing；η_gRepresent gyro The white noise of measurement, ε_aFor the white noise of accelerometer measures；For error parameter vector X first derivative vector；The shape State matrix A is：

In formula (2), The aircraft body system measured for gyro is relative to launching inertial system Angular speed aircraft body system projection,Represent withFor the matrix of the elements in a main diagonal；f^bThe specific force measured for accelerometer aircraft body system projection,Represent SolveAntisymmetric matrix；I_3×3The unit matrix of expression 3 × 3； G represents Newtonian gravitational constant, and M represents earth quality, x,_y, z represents coordinate of the aircraft under launching inertial system, and r represents flight Distance of the device to launching inertial system origin；The transition matrix of launching inertial system is tied to for aircraft body；

The measurement equation is：

In formula (3), Z (t) represents guidance matrix of differences；φⁱ" it is attitude angle and inertial reference calculation posture that star sensor is measured The difference at angle；Z_v(t) it is the GPS aircraft speeds measured and the difference of the speed of inertial reference calculation；Z_r(t) the aircraft position measured for GPS Put the difference with the position of inertial reference calculation；H_9×21Represent measurement matrix；X represents the error parameter vector；V_9×1Represent white noise to Amount；The measurement matrix H_9×21For：

In formula (4), I_3×3The unit matrix of expression 3 × 3；The white noise vector V_9×1For：

In formula (5),WithWhite noise, δ M are measured for the posture of star sensor_x、δM_yWith δ M_zFor GPS amount Degree of testing the speed white noise, δ x_G、δx_GWith δ z_GFor GPS adjustment location white noise；T represents transposition computing.

S120, the observability for differentiating each error parameter.

Optionally, the observability of each error parameter is differentiated, including：

Using the state matrix, transfer matrix is calculated, the calculation formula of the transfer matrix is：

In formula (6), Φ_k,k-1The expression state k-1 moment is to the transfer matrix at state k moment, I_21×21The list of expression 21 × 21 Bit matrix；T₁Represent preset duration；A_kRepresent the state matrix of k-th of preset duration；M is positive integer and m >=2.

Specifically, T₁Preset duration is represented, the preset duration is usually 200ms integral multiple, is 1s to the maximum.For example, If preset duration T₁=200ms, A_kRepresent the state matrix of k-th of 200ms period.

According to the transfer matrix and the measurement matrix, singular value decomposition matrix is treated in calculating, described to treat singular value decomposition The calculation formula of matrix is：

In formula (7), Q represents to treat singular value decomposition matrix, H_9×21Represent measurement matrix, Φ_1,0The moment of expression state 0 is to shape The transfer matrix at the moment of state 1, Φ_m-1,m-2The transfer matrix at expression state m-2 moment to state m-1 moment, T represents that transposition is transported Calculate.

Treat that singular value decomposition matrix carries out singular value decomposition by described, to obtain the first singular value vector, the second singular value Vector and singular value matrix, the formula of the singular value decomposition is：

Q=U ∑s V^T (8)

In formula (8), Q represents to treat singular value decomposition matrix, and ∑ represents singular value matrix, and the elements in a main diagonal of ∑ is σ_i,i ∈[1,m]；U represents the first singular value vector, U=[u_i]=[u₁ u₂ ... u_m]；V represents the second singular value vector, V=[v_i] =[v₁ v₂ ... v_m]；T represents transposition computing.

Utilize the σ_i、u_i、v_iAnd the guidance matrix of differences Z (t), calculate observability discriminant vector Y_21×1, institute The calculation formula for stating observability discriminant vector Y is：

In formula (9), T represents transposition computing.

By the observability discriminant vector Y_21×1Each element be compared with given threshold δ e.

Specifically, due to observability discriminant vector Y_21×1In each element scope between 0~1, thus, it is described Given threshold δ e are chosen as 0.1.In the present embodiment, the given threshold δ e=0.1 will the observability discriminant vector Y_21×1Each element be compared with given threshold 0.1.

As the observability discriminant vector Y_21×1L, the element at the column position of l ∈ [1,21] row the 1st is more than setting During threshold value δ e, l in decision errors parameter vector X, the error parameter Observable at the column position of l ∈ [1,21] row the 1st.

Specifically, in the present embodiment,X represents to miss Poor parameter vector, φ represents the misaligned angle of the platform error, and δ V represent velocity error, and δ r represent site error, δ K_gRepresent the mark of gyro Spend factor error, b₁Represent gyroscopic drift error, δ K_aThe scale factor error of accelerometer is represented,Represent accelerometer Drift error；And each error parameter includes three direction vectors in respective coordinate system.For example, when the observability differentiates Vectorial Y_21×1The column position of the 10th row the 1st at element when being more than given threshold δ e=0.1, the in decision errors parameter vector X The scale factor error δ K of the error parameter Observable at the column position of 10 row the 1st, i.e. gyro_gObservable.In another example, work as institute State observability discriminant vector Y_21×1The column position of the 19th row the 1st at element be more than given threshold δ e=0.1 when, decision errors The drift error of the error parameter Observable in parameter vector X at the column position of the 19th row the 1st, i.e. accelerometerIt is considerable Survey.

S130, when there is observable error parameter, using preset duration as filtering cycle, using Kalman filtering or from Adaptive filtering, calibrates the observable error parameter.

Specifically, in the present embodiment, when there is observable error parameter, with the preset duration in step S120 T₁For filtering cycle, the filtering cycle can be taken as 200ms integral multiple, be 1s to the maximum.Utilize Kalman filtering or adaptive Filtering, calibrates the observable error parameter.

The technical scheme of the present embodiment is by setting up the state equation and measurement equation of aircraft guidance system；Differentiate each The observability of error parameter；When there is observable error parameter, using Kalman filtering or adaptive-filtering, calibrate The observable error parameter, can overcome the shortcomings of existing guidance instrument error isolation technics " world is inconsistent ", can In real time, on-orbit calibration goes out the error coefficient of Guidance instrumentation, and algorithm is simple, is easy to engineering；Inertia device can be reduced and safeguard guarantor Barrier expense, due to the ability with on-orbit calibration, so as to reduce the number of times of inertia system ground periodic calibrating, saves manpower Material resources；Aircraft guidance precision can be improved, the dependence to satellite navigation is reduced, extending space aircraft performs the scope of task And ability.

Embodiment two

Fig. 2 is a kind of structure chart of error parameter caliberating device in the embodiment of the present invention two, and the present embodiment is applicable to Need the situation to the error parameter progress on-orbit calibration of Guidance instrumentation.With reference to Fig. 2, the error parameter demarcation that the present embodiment is provided Device specifically can be as follows：

Build module 210, state equation and measurement equation for setting up aircraft guidance system, the state equation and Comprising error parameter vector in measurement equation, the error parameter vector is made up of multiple error parameters；

Discrimination module 220, the observability for differentiating each error parameter；

Demarcating module 230, for when there is observable error parameter, using preset duration as filtering cycle, utilizing card Kalman Filtering or adaptive-filtering, calibrate the observable error parameter.

Optionally, the error parameter include the misaligned angle of the platform error, velocity error, site error, the scale of gyro because Count the drift error of error, gyroscopic drift error, the scale factor error of accelerometer and accelerometer；Each mistake Poor parameter is vector, and includes three durection components.

Optionally, the state equation is：

In formula (1), A represents state matrix；The transition matrix of launching inertial system is tied to for aircraft body；X represents to miss Poor parameter vector,φ represents the misaligned angle of the platform error, δ V Velocity error is represented, δ r represent site error, δ K_gRepresent the scale factor error of gyro, b₁Represent gyroscopic drift error, δ K_aTable Show the scale factor error of accelerometer,The drift error of accelerometer is represented, T represents transposition computing；η_gRepresent that gyro is surveyed The white noise of amount, ε_aFor the white noise of accelerometer measures；For error parameter vector X first derivative vector；The state square Battle array A be：

In formula (2), The aircraft body system measured for gyro is relative to launching inertial system Angular speed aircraft body system projection,Represent withFor the matrix of the elements in a main diagonal；f^bThe specific force measured for accelerometer aircraft body system projection,Represent SolveAntisymmetric matrix；I_3×3The unit matrix of expression 3 × 3； G represents Newtonian gravitational constant, and M represents earth quality, and x, y, z represents coordinate of the aircraft under launching inertial system, and r represents flight Distance of the device to launching inertial system origin；The transition matrix of launching inertial system is tied to for aircraft body；

The measurement equation is：

In formula (3), Z (t) represents guidance matrix of differences；φⁱ" it is attitude angle and inertial reference calculation posture that star sensor is measured The difference at angle；Z_v(t) it is the GPS aircraft speeds measured and the difference of the speed of inertial reference calculation；Z_r(t) the aircraft position measured for GPS Put the difference with the position of inertial reference calculation；H_9×21Represent measurement matrix；X represents the error parameter vector；V_9×1Represent white noise to Amount；The measurement matrix H_9×21For：

In formula (4), I_3×3The unit matrix of expression 3 × 3；The white noise vector V_9×1For：

In formula (5),WithWhite noise, δ M are measured for the posture of star sensor_x、δM_yWith δ M_zFor GPS amount Degree of testing the speed white noise, δ x_G、δx_GWith δ z_GFor GPS adjustment location white noise；T represents transposition computing.

Optionally, the discrimination module includes：

First computing unit, for utilizing the state matrix, calculates transfer matrix, the calculation formula of the transfer matrix For：

In formula (6), Φ_k,k-1The expression state k-1 moment is to the transfer matrix at state k moment, I_21×21The list of expression 21 × 21 Bit matrix；T₁Represent preset duration；A_kRepresent the state matrix of k-th of preset duration；M is positive integer and m >=2；

Second computing unit, for according to the transfer matrix and the measurement matrix, calculating and treating singular value decomposition matrix, It is described treat singular value decomposition matrix calculation formula be：

In formula (7), Q represents to treat singular value decomposition matrix, H_9×21Represent measurement matrix, Φ_1,0The moment of expression state 0 is to shape The transfer matrix at the moment of state 1, Φ_m-1,m-2The transfer matrix at expression state m-2 moment to state m-1 moment, T represents that transposition is transported Calculate；

Resolving cell, for by it is described treat singular value decomposition matrix carry out singular value decomposition, with obtain the first singular value to Amount, the second singular value vector and singular value matrix, the formula of the singular value decomposition is：

Q=U ∑s V^T (8)

In formula (8), Q represents to treat singular value decomposition matrix, and ∑ represents singular value matrix, and the elements in a main diagonal of ∑ is σ_i,i ∈[1,m]；U represents the first singular value vector, U=[u_i]=[u₁ u₂ ... u_m]；V represents the second singular value vector, V=[v_i] =[v₁ v₂ ... v_m]；T represents transposition computing；

3rd computing unit, for utilizing the σ_i、u_i、v_iAnd the guidance matrix of differences Z (t), calculate Observable Property discriminant vector Y_21×1, the calculation formula of the observability discriminant vector Y is：

In formula (9), T represents transposition computing；

Comparing unit, for by the observability discriminant vector Y_21×1Each element compared with given threshold δ e Compared with；

Identifying unit, for as the observability discriminant vector Y_21×1L, at the column position of l ∈ [1,21] row the 1st When element is more than given threshold δ e, l in decision errors parameter vector X, the error at the column position of l ∈ [1,21] row the 1st Parameter Observable.

The error parameter caliberating device that the present embodiment is provided, is demarcated with the error parameter that any embodiment of the present invention is provided Method belongs to same inventive concept, can perform the error parameter scaling method that any embodiment of the present invention is provided, possesses execution The corresponding functional module of error parameter scaling method and beneficial effect.Not ins and outs of detailed description in the present embodiment, can The error parameter scaling method provided referring to any embodiment of the present invention.

Note, above are only presently preferred embodiments of the present invention and institute's application technology principle.It will be appreciated by those skilled in the art that The invention is not restricted to specific embodiment described here, can carry out for a person skilled in the art it is various it is obvious change, Readjust and substitute without departing from protection scope of the present invention.Therefore, although the present invention is carried out by above example It is described in further detail, but the present invention is not limited only to above example, without departing from the inventive concept, also Other more equivalent embodiments can be included, and the scope of the present invention is determined by scope of the appended claims.

Claims (8)

1. a kind of error parameter scaling method, it is characterised in that including：

Set up and error is included in the state equation and measurement equation of aircraft guidance system, the state equation and measurement equation Parameter vector, the error parameter vector is made up of multiple error parameters；

Differentiate the observability of each error parameter；

When there is observable error parameter, using preset duration as filtering cycle, using Kalman filtering or adaptive-filtering, Calibrate the observable error parameter.

2. according to the method described in claim 1, it is characterised in that the error parameter includes the misaligned angle of the platform error, speed Error, site error, the scale factor error of gyro, gyroscopic drift error, the scale factor error of accelerometer and acceleration Spend the drift error of meter；Each error parameter is vector, and includes three durection components.

3. according to the method described in claim 1, it is characterised in that the state equation is：

$\stackrel{\·}{X}=AX+\left[\begin{array}{c}\begin{array}{cc}{\left(C_b^i\right)}_{3\×3}& 0_{3\×3}\\ 0_{3\×3}& {\left(C_b^i\right)}_{3\×3}\end{array}\\ 0_{3\×6}\\ 0_{3\×6}\\ 0_{3\×6}\\ 0_{3\×6}\\ 0_{3\×6}\end{array}\right]\·\left[\begin{array}{c}{\mathrm{\η}}_g\\ {\mathrm{\ϵ}}_a\end{array}\right]---\left(1\right)$

In formula (1), A represents state matrix；The transition matrix of launching inertial system is tied to for aircraft body；X represents that error is joined Number vector,φ represents the misaligned angle of the platform error, and δ V are represented Velocity error, δ r represent site error, δ K_gRepresent the scale factor error of gyro, b₁Represent gyroscopic drift error, δ K_aRepresent to add The scale factor error of speedometer,The drift error of accelerometer is represented, T represents transposition computing；η_gRepresent gyro to measure White noise, ε_aFor the white noise of accelerometer measures；For error parameter vector X first derivative vector；The state matrix A For：

$A=\left[\begin{array}{c}\begin{array}{ccccccc}{0}_{3\×3}& 0_{3\×3}& 0_{3\×3}& {\left(F_1\right)}_{3\×3}& {\left(C_b^i\right)}_{3\×3}& 0_{3\×3}& 0_{3\×3}\\ {\left(F_2\right)}_{3\×3}& 0_{3\×3}& {\left(N^i\right)}_{3\×3}& 0_{3\×3}& 0_{3\×3}& {\left(F_3\right)}_{3\×3}& {\left(C_b^i\right)}_{3\×3}\\ 0_{3\×3}& I_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}\end{array}\\ 0_{3\×21}\\ 0_{3\×21}\\ 0_{3\×21}\\ 0_{3\×21}\end{array}\right]---\left(2\right)$

In formula (2), Angular speed of the aircraft body system relative to launching inertial system measured for gyro In the projection of aircraft body system,Represent withFor the matrix of the elements in a main diagonal； f^bThe specific force measured for accelerometer aircraft body system projection,Represent to solveAntisymmetric matrix； I_3×3The unit matrix of expression 3 × 3；G represents Newtonian gravitational constant, M Earth quality is represented, x, y, z represents coordinate of the aircraft under launching inertial system, and r represents aircraft to launching inertial system origin Distance；The transition matrix of launching inertial system is tied to for aircraft body；

The measurement equation is：

$Z\left(t\right)=\left[\begin{array}{c}{\mathrm{\φ}}^{i^{\′\′}}\\ Z_v\left(t\right)\\ Z_r\left(t\right)\end{array}\right]=H_{9\×21}X+V_{9\×1}---\left(3\right)$

In formula (3), Z (t) represents guidance matrix of differences；φ^i”For star sensor measure attitude angle and inertial reference calculation attitude angle it Difference；Z_v(t) it is the GPS aircraft speeds measured and the difference of the speed of inertial reference calculation；Z_r(t) for GPS measure position of aircraft and The difference of the position of inertial reference calculation；H_9×21Represent measurement matrix；X represents the error parameter vector；V_9×1Represent white noise vector； The measurement matrix H_9×21For：

$H_{9\×21}=\left[\begin{array}{ccccccc}{I}_{3\×3}& 0& 0& 0& 0& 0& 0\\ 0& I_{3\×3}& 0& 0& 0& 0& 0\\ 0& 0& I_{3\×3}& 0& 0& 0& 0\end{array}\right]---\left(4\right)$

In formula (4), I_3×3The unit matrix of expression 3 × 3；The white noise vector V_9×1For：

$V_{9\×1}={\left[\begin{array}{ccccccccc}-V_x^i& -V_y^i& -V_z^i& {\mathrm{\δM}}_x& {\mathrm{\δM}}_y& {\mathrm{\δM}}_z& {\mathrm{\δx}}_G& {\mathrm{\δy}}_G& {\mathrm{\δz}}_G\end{array}\right]}^T---\left(5\right)$

In formula (5),WithWhite noise, δ M are measured for the posture of star sensor_x、δM_yWith δ M_zFor GPS measurement speed Spend white noise, δ x_G、δx_GWith δ z_GFor GPS adjustment location white noise；T represents transposition computing.

4. method according to claim 3, it is characterised in that differentiate the observability of each error parameter, including：

Using the state matrix, transfer matrix is calculated, the calculation formula of the transfer matrix is：

${\mathrm{\Φ}}_{k,k-1}=I_{21\×21}+\underset{j=1}{\overset{4}{\Σ}}\frac{{A_k}^j{T_1}^j}{j!},k\∈\[1,...,m-1\]---\left(6\right)$

In formula (6), Φ_k,k-1The expression state k-1 moment is to the transfer matrix at state k moment, I_21×21The unit square of expression 21 × 21 Battle array；T₁Represent preset duration；A_kRepresent the state matrix of k-th of preset duration；M is positive integer and m >=2；

According to the transfer matrix and the measurement matrix, singular value decomposition matrix is treated in calculating, described to treat singular value decomposition matrix Calculation formula be：

$Q={\left[\begin{array}{cccc}{H}_{9\×21}^T& {\left(H_{9\×21}{\mathrm{\Φ}}_{1,0}\right)}^T& ...& {(H_{9\×2}{\mathrm{\Φ}}_{m-1,m-2}...{\mathrm{\Φ}}_{1,0})}^T\end{array}\right]}^T---\left(7\right)$

In formula (7), Q represents to treat singular value decomposition matrix, H_9×21Represent measurement matrix, Φ_1,0When the moment of expression state 0 is to state 1 The transfer matrix at quarter, Φ_m-1,m-2The transfer matrix at expression state m-2 moment to state m-1 moment, T represents transposition computing；

By it is described treat singular value decomposition matrix carry out singular value decomposition, with obtain the first singular value vector, the second singular value vector, And singular value matrix, the formula of the singular value decomposition is：

Q=U ∑s V^T (8)

In formula (8), Q represents to treat singular value decomposition matrix, and ∑ represents singular value matrix, and the elements in a main diagonal of ∑ is σ_i,i∈[1, m]；U represents the first singular value vector, U=[u_i]=[u₁ u₂ ... u_m]；V represents the second singular value vector, V=[v_i]=[v₁ v₂ ... v_m]；T represents transposition computing；

Utilize the σ_i、u_i、v_iAnd the guidance matrix of differences Z (t), calculate observability discriminant vector Y_21×1, it is described can Observation discriminant vector Y calculation formula is：

$Y_{21\×1}=\underset{i=1}{\overset{m}{\Σ}}\[\frac{u_i^{T}Z\left(t\right)}{{\mathrm{\σ}}_i}\]v_i---\left(9\right)$

In formula (9), T represents transposition computing；

By the observability discriminant vector Y_21×1Each element be compared with given threshold δ e；

As the observability discriminant vector Y_21×1L, the element at the column position of l ∈ [1,21] row the 1st is more than given threshold δ During e, l in decision errors parameter vector X, the error parameter Observable at the column position of l ∈ [1,21] row the 1st.

5. a kind of error parameter caliberating device, it is characterised in that including：

Build module, state equation and measurement equation for setting up aircraft guidance system, the state equation and measurement side Cheng Zhongjun includes error parameter vector, and the error parameter vector is made up of multiple error parameters；

Discrimination module, the observability for differentiating each error parameter；

Demarcating module, for when there is observable error parameter, using preset duration as filtering cycle, utilizing Kalman filtering Or adaptive-filtering, calibrate the observable error parameter.

6. device according to claim 5, it is characterised in that the error parameter includes the misaligned angle of the platform error, speed Error, site error, the scale factor error of gyro, gyroscopic drift error, the scale factor error of accelerometer and acceleration Spend the drift error of meter；Each error parameter is vector, and includes three durection components.

7. device according to claim 5, it is characterised in that the state equation is：

$\stackrel{\·}{X}=AX+\left[\begin{array}{c}\begin{array}{cc}{\left(C_b^i\right)}_{3\×3}& 0_{3\×3}\\ 0_{3\×3}& {\left(C_b^i\right)}_{3\×3}\end{array}\\ 0_{3\×6}\\ 0_{3\×6}\\ 0_{3\×6}\\ 0_{3\×6}\\ 0_{3\×6}\end{array}\right]\·\left[\begin{array}{c}{\mathrm{\η}}_g\\ {\mathrm{\ϵ}}_a\end{array}\right]---\left(1\right)$

In formula (1), A represents state matrix；The transition matrix of launching inertial system is tied to for aircraft body；X represents that error is joined Number vector,φ represents the misaligned angle of the platform error, and δ V are represented Velocity error, δ r represent site error, δ K_gRepresent the scale factor error of gyro, b₁Represent gyroscopic drift error, δ K_aRepresent to add The scale factor error of speedometer,The drift error of accelerometer is represented, T represents transposition computing；η_gRepresent gyro to measure White noise, ε_aFor the white noise of accelerometer measures；For error parameter vector X first derivative vector；The state matrix A For：

$A=\left[\begin{array}{c}\begin{array}{ccccccc}{0}_{3\×3}& 0_{3\×3}& 0_{3\×3}& {\left(F_1\right)}_{3\×3}& {\left(C_b^i\right)}_{3\×3}& 0_{3\×3}& 0_{3\×3}\\ {\left(F_2\right)}_{3\×3}& 0_{3\×3}& {\left(N^i\right)}_{3\×3}& 0_{3\×3}& 0_{3\×3}& {\left(F_3\right)}_{3\×3}& {\left(C_b^i\right)}_{3\×3}\\ 0_{3\×3}& I_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}\end{array}\\ 0_{3\×21}\\ 0_{3\×21}\\ 0_{3\×21}\\ 0_{3\×21}\end{array}\right]---\left(2\right)$

In formula (2), Angular speed of the aircraft body system relative to launching inertial system measured for gyro In the projection of aircraft body system,Represent withFor the matrix of the elements in a main diagonal； f^bThe specific force measured for accelerometer aircraft body system projection,Represent to solveAntisymmetric matrix； I_3×3The unit matrix of expression 3 × 3；G represents Newtonian gravitational constant, M Earth quality is represented, x, y, z represents coordinate of the aircraft under launching inertial system, and r represents aircraft to launching inertial system origin Distance；The transition matrix of launching inertial system is tied to for aircraft body；

The measurement equation is：

$Z\left(t\right)=\left[\begin{array}{c}{\mathrm{\φ}}^{i^{\′\′}}\\ Z_v\left(t\right)\\ Z_r\left(t\right)\end{array}\right]=H_{9\×21}X+V_{9\×1}---\left(3\right)$

In formula (3), Z (t) represents guidance matrix of differences；φ^i”For star sensor measure attitude angle and inertial reference calculation attitude angle it Difference；Z_v(t) it is the GPS aircraft speeds measured and the difference of the speed of inertial reference calculation；Z_r(t) for GPS measure position of aircraft and The difference of the position of inertial reference calculation；H_9×21Represent measurement matrix；X represents the error parameter vector；V_9×1Represent white noise vector； The measurement matrix H_9×21For：

$H_{9\×21}=\left[\begin{array}{ccccccc}{I}_{3\×3}& 0& 0& 0& 0& 0& 0\\ 0& I_{3\×3}& 0& 0& 0& 0& 0\\ 0& 0& I_{3\×3}& 0& 0& 0& 0\end{array}\right]---\left(4\right)$

In formula (4), I_3×3The unit matrix of expression 3 × 3；The white noise vector V_9×1For：

$V_{9\×1}={\left[\begin{array}{ccccccccc}-V_x^i& -V_y^i& -V_z^i& {\mathrm{\δM}}_x& {\mathrm{\δM}}_y& {\mathrm{\δM}}_z& {\mathrm{\δx}}_G& {\mathrm{\δy}}_G& {\mathrm{\δz}}_G\end{array}\right]}^T---\left(5\right)$

In formula (5),WithWhite noise, δ M are measured for the posture of star sensor_x、δM_yWith δ M_zFor GPS measurement speed Spend white noise, δ x_G、δx_GWith δ z_GFor GPS adjustment location white noise；T represents transposition computing.

8. device according to claim 7, it is characterised in that the discrimination module includes：

First computing unit, for utilizing the state matrix, calculates transfer matrix, the calculation formula of the transfer matrix is：

${\mathrm{\Φ}}_{k,k-1}=I_{21\×21}+\underset{j=1}{\overset{4}{\Σ}}\frac{{A_k}^j{T}_1^j}{j!},k\∈\[1,...,m-1\]---\left(6\right)$

In formula (6), Φ_k,k-1The expression state k-1 moment is to the transfer matrix at state k moment, I_21×21The unit square of expression 21 × 21 Battle array；T₁Represent preset duration；A_kRepresent the state matrix of k-th of preset duration；M is positive integer and m >=2；

Second computing unit, for according to the transfer matrix and the measurement matrix, singular value decomposition matrix to be treated in calculating, described The calculation formula for treating singular value decomposition matrix is：

$Q={\left[\begin{array}{cccc}{H}_{9\×21}^T& {\left(H_{9\×21}{\mathrm{\Φ}}_{1,0}\right)}^T& ...& {(H_{9\×2}{\mathrm{\Φ}}_{m-1,m-2}...{\mathrm{\Φ}}_{1,0})}^T\end{array}\right]}^T---\left(7\right)$

In formula (7), Q represents to treat singular value decomposition matrix, H_9×21Represent measurement matrix, Φ_1,0When the moment of expression state 0 is to state 1 The transfer matrix at quarter, Φ_m-1,m-2The transfer matrix at expression state m-2 moment to state m-1 moment, T represents transposition computing；

Resolving cell, for treating that singular value decomposition matrix carries out singular value decomposition by described, to obtain the first singular value vector, the Two singular value vectors and singular value matrix, the formula of the singular value decomposition is：

Q=U ∑s V^T (8)

In formula (8), Q represents to treat singular value decomposition matrix, and ∑ represents singular value matrix, and the elements in a main diagonal of ∑ is σ_i,i∈[1, m]；U represents the first singular value vector, U=[u_i]=[u₁ u₂ ... u_m]；V represents the second singular value vector, V=[v_i]=[v₁ v₂ ... v_m]；T represents transposition computing；

3rd computing unit, for utilizing the σ_i、u_i、v_iAnd the guidance matrix of differences Z (t), calculate observability and sentence Not vector Y_21×1, the calculation formula of the observability discriminant vector Y is：

$Y_{21\×1}=\underset{i=1}{\overset{m}{\Σ}}\[\frac{u_i^{T}Z\left(t\right)}{{\mathrm{\σ}}_i}\]v_i---\left(9\right)$

In formula (9), T represents transposition computing；

Comparing unit, for by the observability discriminant vector Y_21×1Each element be compared with given threshold δ e；

Identifying unit, for as the observability discriminant vector Y_21×1L, the element at the column position of l ∈ [1,21] row the 1st During more than given threshold δ e, l in decision errors parameter vector X, the error parameter at the column position of l ∈ [1,21] row the 1st Observable.