CN102252677A

CN102252677A - Time series analysis-based variable proportion self-adaptive federal filtering method

Info

Publication number: CN102252677A
Application number: CN2011100958407A
Authority: CN
Inventors: 李宁; 袁克非; 袁赣南; 刘利强; 张勇刚
Original assignee: Harbin Engineering University
Current assignee: Harbin Engineering University
Priority date: 2011-04-18
Filing date: 2011-04-18
Publication date: 2011-11-23

Abstract

The invention discloses a time series analysis-based variable proportion self-adaptive federal filtering method, and the method is utilized for underwater multi-sensor integrated navigation systems. The time series analysis-based variable proportion self-adaptive federal filtering method is characterized by establishing system state equations and measurement equations according to error equations belonging to all navigation sensor systems, carrying out discretization processing of the established equations, creating discrete state-space models respectively corresponding to the all navigation sensor systems, acquiring information weight of the navigation sensor systems through an autoregressive model according to historical data of the navigation sensor systems, acquiring an information distribution ratio according to the obtained information weight and a law of information conservation, realizing global optimal estimates, and resetting a filtering value and an evaluated error covariance matrix by the global optimal estimates. The time series analysis-based variable proportion self-adaptive federal filtering method improves system navigation precision, system stability and fault tolerance, and can satisfy requirements which belong to underwater navigation devices and comprise high precision and high reliability requirements.

Description

Time sequence analysis-based variable-proportion self-adaptive federal filtering method

Technical Field

The invention belongs to the technical field of multi-sensor integrated navigation systems, relates to a filtering method of a multi-sensor integrated navigation system, and particularly relates to a variable proportion self-adaptive federal filtering method based on time series analysis.

Background

With the progress of navigation technology, the performance and application range of any single navigation device show obvious limitations, cannot meet the increasing precision requirement of a vehicle, cannot adapt to a complex application environment, and cannot completely meet the requirement of system reliability. Horst Ahlers states: "future belongs to multisensor. "multi-sensor integrated navigation system" is becoming a necessary trend. With the improvement of hardware conditions, the fusion method becomes an important factor that restricts the performance of the combined system. The Federal Kalman filtering method developed by N.A.Carlson is a fusion method successfully applied to engineering practice, and adopts a noise variance upper bound technology to eliminate the sub-filtering correlation so as to achieve global optimization.

However, in the existing federal filter structure, the information distribution proportion is fixed along with the determination of the structure, and the information distribution proportion cannot reflect the changes of the working state and the data state of the navigation sensor, accurately reflect the information weight distribution of the navigation sensor, accurately and truly reflect the condition of the system, and is not beneficial to the research on improving the precision, the reliability and the fault tolerance of the system.

Disclosure of Invention

The invention provides a variable proportion self-adaptive federal filtering method based on time sequence analysis, aiming at the problem that the information distribution proportion in the existing federal filtering method can not accurately reflect the information weight distribution of a navigation sensor along with the determination of the structure, so that the information distribution proportion can be adapted to the working condition of the navigation sensor and the running environment condition of a carrier, and the purposes of improving the precision, the reliability and the fault tolerance of a system are achieved.

A time sequence analysis-based variable proportion self-adaptive federal filtering method specifically comprises the following steps:

firstly, establishing a system state equation and a measurement equation according to an error equation of each navigation sensor system of the integrated navigation system;

discretizing the obtained system state equation and measurement equation, and establishing a discrete state space model corresponding to each navigation sensor system;

thirdly, acquiring an information weight of the navigation sensor system through an autoregressive model according to historical data of the navigation sensor system;

step four, obtaining an information distribution proportion according to the information weight and an information conservation law;

and fifthly, carrying out federal filtering according to the information distribution proportion to obtain global optimal estimation, and finally resetting the filtering value and the estimation error variance matrix by using the global optimal estimation.

The integrated navigation system in the first step is composed of an inertial navigation system INS, an electrostatic gyro monitor ESGM, a satellite positioning navigation system GPS and a Doppler velocity measurement system DVL according to the configuration condition of a navigation sensor of an underwater vehicle.

Respectively establishing an INS state equation, an INS/ESGM measurement equation, an INS/GPS measurement equation, an INS/DVL system state equation and a measurement equation according to each system error equation:

the INS state equation is as follows:

wherein the state quantity

δ L, δ λ, δ h represent changes in latitude L, longitude λ, altitude h, respectively, and δ V_E、δV_N、δV_URespectively represent east velocity V_EVelocity V in the north direction_NVelocity in the direction of the sky V_UChange of (phi)_E、φ_N、φ_URespectively representing pitch angle, roll angle, yaw angle, epsilon_bx、ε_by、ε_bzIs the component of gyro drift random constant error on three axes, epsilon_rx、ε_ry、ε_rzFor the components of the gyro drift markov process error in the three axes,

is the component of the accelerometer error in three axes; process noise array W ═ W_bx w_by w_bz w_gx w_gy w_gz w_ax w_ay w_az]^T，w_bx、w_by、w_bzComponent on three axes of white noise, w, representing the gyro drift random constant_gx、w_gy、w_gzComponent of white noise on three axes, w, representing a gyro-drift Markov process_ax、w_ay、w_azA component of white noise on three axes representing accelerometer error; state transition arrayWherein

f_E、f_N、f_URespectively an east acceleration, a north acceleration, a sky acceleration, omega_ieThe rotational angular velocity of the earth; noise driving array

Wherein

Is a direction cosine matrix.

The INS/ESGM measurement equation is as follows: z₁(t)＝H₁(t)X₁(t)+V₁(t) wherein Z₁The difference value of the ESGM and the INS position and course information is an observed quantity V₁(t)＝[m₁，m₂，m₃，m₄]^TIs a zero mean value generated by the position error and the orientation error of the ESGMWhite gaussian noise of (1); h₁(t) is the measurement array, INS/ESGM State equation X₁(t)＝X(t)。

The INS/GPS measurement equation is as follows:

wherein the position quantity is measured

Measurement of speed v₂(t) white noise, INS/GPS equation of state X₂(t)＝X(t)；

The INS/DVL system state equation is as follows:

the INS/DVL measurement equation is as follows:

wherein, the noise of the Doppler velocity measurement system drives the array G_d(t)＝diag[1 1 0]，F_d(t) is the augmented state transition array of the Doppler velocimetry system, W_d(t) noise matrix of the augmentation Process of Doppler velocimetry System, X_d(t) is the Doppler velocimetry system augmented state quantity, X_d＝[δV_d δΔ δC]^T，δV_dThe measurement error is a speed deviation error in the Doppler velocity measurement system, delta is a drift angle error in the Doppler velocity measurement system, and delta C is a scale coefficient error in the Doppler velocity measurement system;

V₃(t)＝[m_VE m_VN]^Tto observe white noise, V_dE、V_dNRespectively represents the east speed and the north speed in the Doppler velocity measurement system, and K' is the heading angle of the ship.

In the second step, the discrete state space model corresponding to each navigation sensor system is:

where k represents the index value of the discrete time state, X_kTo representState quantity at time k, Z_kRepresents the observed quantity at time k, H_kA measuring matrix representing the time k, Γ_k/k-1Is a noise-driven array, W_k-1Is the system state noise at time k-1 and the process noise covariance matrix Q_k＝E[W_kW_k ^T]，V_kFor the system at time k, measure noise, and measure noise covariance matrix R_k＝E[V_kV_k ^T]D, state transition matrix phi_k/k-1The method is obtained by adopting a step-by-step cumulative discretization method:

where J is T/s, T is sampling time, s is step length divided in sampling time, I is unit array, F_i-1Is the state transition matrix at time i-1.

In the third step, the prediction error of each navigation sensor system is used as the information weight theta_i：

Wherein, y_i(τ) is the output sequence of the ith navigation sensor,for the predicted value of the output sequence of the ith navigation sensor, i is 0, 1, 2 and 3 respectively correspond to the inertial navigation system, the electrostatic gyro monitor, the satellite positioning navigation system and the DopplerA speed measuring system.

The information distribution proportion is in direct proportion to the information weight, and the weight sum of each navigation sensor system is

Wherein beta is₀Is the primary system weight, beta_iThe weight of each subsystem is obtained, n is the number of all subsystems, the main system is an inertial navigation system, the subsystems are an electrostatic gyro monitor, a satellite positioning navigation system and a Doppler velocity measurement system, and finally the weight of each navigation sensor system is obtained as follows:

the global optimal estimation in the fifth step is as follows:

wherein,

is a global estimate, P_gTo estimate the mean square error, P_i，kAn estimation error variance matrix P representing the ith navigation sensor system_kBy the ratio beta_iDividing the weighted estimation error variance matrix,

indicating the predicted state quantity of the ith navigation sensor system.

And resetting the filtering value and the estimation error variance matrix by using the global optimal solution.

Because the upper bound technology of variance is used when the system noise is distributed, all subsystems are not related, and the augmentation system is not coupled. Although the filtering results of the subsystems are suboptimal, the fusion is resynthesized so that the global estimate is optimal.

The method has the advantages and positive effects that:

(1) the selection of the information right is based on the source information of each navigation sensor system, so that errors caused by evaluating each navigation sensor system by a filter estimation value are avoided, and the structure of the integrated navigation system is more suitable for the real condition of the system; meanwhile, measurement forward data are used for obtaining the weight, so that the historical data of the navigation sensor system are fully utilized, and the utilization rate of information is improved.

(2) The method of the invention ensures that the navigation precision of the integrated navigation system is kept at a considerable level, the error is small and non-divergent relative to the INS error, and the defects of INS error accumulation, especially position error, are overcome.

(3) The method can keep the performance of the integrated navigation system stable under the condition of different navigation sensor configurations, so that the integrated navigation system has good filtering precision, good positioning capability and error correction capability, stronger anti-jamming capability and noise suppression capability, and the reliability and the usability of the integrated navigation system are improved.

Drawings

FIG. 1 is a schematic diagram of the Federal Filtering method of the present invention applied to a multi-sensor integrated navigation system;

FIG. 2 is a flow chart of the steps of the federated filtering method of the present invention;

FIG. 3 is a schematic representation of a federated filtering architecture employing the federated filtering method of the present invention;

FIG. 4 is a diagram illustrating a variation curve of an information distribution ratio;

FIG. 5 is a schematic diagram of a fusion system and INS system position error versus curve;

FIG. 6 is a schematic diagram of a velocity error versus INS system for a fusion system;

FIG. 7 is a schematic diagram of a position error versus curve for a fusion system applying the method of the present invention and a fusion system applying the classical federated filtering method.

Detailed Description

The present invention will be described in further detail with reference to the accompanying drawings and examples.

The multi-sensor combined Navigation System applied in the embodiment of the invention is an underwater multi-sensor combined Navigation System, and the combined Navigation System consists of an Inertial Navigation System (INS for short), an Electrostatic gyro Monitor (ESGM for short), a satellite Positioning Navigation System (GPS for short) and a Doppler Velocity measurement System (DVL for short).

The INS is composed of a gyroscope and an accelerometer, is an autonomous navigation system independent of external information, can provide various navigation parameters including speed, position, attitude and the like, has the advantages of interference resistance and all weather, and is suitable for serving as a reference navigation system of an underwater vehicle. The precision of the INS mainly depends on inertial devices, and navigation precision errors are increased along with the accumulation of time, so that the problem of error stability of long-time navigation needs to be considered. The ESGM is a dumbbell-type combination body formed by two-degree-of-freedom electrostatic gyroscopes and an indirect stable platform, and is used as an information source monitoring and compensating inertial navigation system with high precision and long-term stability. The ESGM utilizes the star navigation principle to monitor the INS parameters, and the INS provides an indirect stable platform and a solution parameter for the ESGM. The combined mode can greatly inhibit the divergence of navigation errors, prolong the readjustment period of the system and ensure the long-time high-precision cruising of the underwater vehicle. The GPS has the obvious advantages of real-time navigation, irrelevant positioning error and time and higher positioning and speed measurement precision. When the antenna floats out of the water surface or the detection buoy is thrown out, the accurate position of the current point is obtained, the navigation error is calibrated, and the system time is corrected. Based on Doppler effect, the DVL mostly adopts four-beam Doppler sonar with fixed transmitting direction, which not only can compensate attitude error of a carrier in drift angle measurement, but also can measure the ground speed of the carrier, thereby giving space absolute speed vector of the carrier.

As shown in fig. 1, navigation information measured by the integrated navigation system is preprocessed, mainly by performing outlier processing on data by a distance function method, and the like, then performing time and space calibration on the acquired navigation data, then fusing the calibrated effective data by a variable-proportion adaptive federal filtering method to obtain an optimal solution, and finally outputting the navigation data.

The invention discloses a time series analysis-based variable-proportion self-adaptive federal filtering method, which comprises the following steps as shown in figure 2.

Step one, establishing a system state equation and a measurement equation according to error equations of all navigation sensors of the integrated navigation system. In this embodiment, an INS state equation, an INS/ESGM measurement equation, an INS/GPS measurement equation, and an INS/DVL system state equation and measurement equation need to be established.

And A, establishing an INS state equation.

The inertial navigation system INS errors include three-channel position, velocity, attitude errors and inertial device errors. The three-channel position error is represented by latitude L, longitude lambda and altitude h, and the three-channel speed error is represented by east speed V_EVelocity V in the north direction_NVelocity V in the direction of the sky_UAnd characterizing, wherein attitude errors of the three channels are characterized by a pitch angle, a roll angle and a yaw angle, and error equations are respectively as follows:

<math> <mrow> <mi>δ</mi> <mover> <mi>λ</mi> <mo>·</mo> </mover> <mo>=</mo> <mfrac> <mrow> <mi>δ</mi> <msub> <mi>V</mi> <mi>E</mi> </msub> </mrow> <mrow> <msub> <mi>R</mi> <mi>N</mi> </msub> <mo>+</mo> <mi>h</mi> </mrow> </mfrac> <mi>sec</mi> <mi>L</mi> <mo>+</mo> <mfrac> <msub> <mi>V</mi> <mi>E</mi> </msub> <mrow> <msub> <mi>R</mi> <mi>N</mi> </msub> <mo>+</mo> <mi>h</mi> </mrow> </mfrac> <mi>sec</mi> <mi>LtgLδL</mi> </mrow> </math>

wherein,

latitude L, longitude lambda, altitude h, east speed V_EVelocity V in the north direction_NVelocity in the direction of the sky V_UThe rate of change of pitch angle, roll angle, yaw angle; δ L, δ λ, δ h represent changes in latitude L, longitude λ, altitude h, respectively, and δ V_E、δV_N、δV_URespectively represent east velocity V_EVelocity V in the north direction_NVelocity in the direction of the sky V_UA change in (c); f. of_E、f_N、f_URespectively an east acceleration, a north acceleration, a sky acceleration, omega_ieIs the rotational angular velocity of the earth, R_MRadius of curvature of meridian, R_NThe curvature radius of the unitary mortise ring.

The inertial device errors comprise installation errors, scale coefficient errors, random errors and the like, and the installation errors and the scale coefficient errors are different due to different conditions of device types, models, environments, processes and the like and have no uniformity, so that only the random errors are considered in the invention. The random error needs to take into account gyro drift and accelerometer error.

The gyro drift is: e ═ e-_b+ε_r+ε_g(ii) a Wherein epsilon_bDenotes a random constant, ∈_rRepresenting a first order Markov process,. epsilon_gWhite noise due to gyro drift. The triaxial gyro drift equation is:

wherein T is_gRepresenting the correlation time, w, of the three-axis gyro drift_gWhite noise representing a gyro drift markov process.

The accelerometer error can be considered as a first order markov process, and the error equation of the accelerometer is as follows:

wherein T is_aTime of correlation, w, representing accelerometer error_aWhite noise representing accelerometer error.

The INS state equation can be determined by an error equation of an inertial navigation system INS as follows:

wherein the state quantity

Including three channel position, velocity and attitude, gyro error and accelerometer error, epsilon_bx、ε_by、ε_bzIs the component of gyro drift random constant error on three axes, epsilon_rx、ε_ry、ε_rzFor the components of the gyro drift markov process error in the three axes,

the components of the accelerometer error are in three axes, namely the x-axis, the y-axis and the z-axis.

Process noise array W ═ W_bx w_by w_bz w_gx w_gy w_gz w_ax w_ay w_az]^T，w_bx、w_by、w_bzComponent on three axes of white noise, w, representing the gyro drift random constant, respectively_gx、w_gy、w_gzComponent of white noise on three axes, w, representing the gyro-drift Markov process_ax、w_ay、w_azThe components of white noise on the three axes representing the accelerometer error, respectively.

Noise driving array

Wherein

The direction cosine matrix is b, a ship body coordinate system is fixedly connected to a ship body, and n is a navigation coordinate system and is selected according to the working state of the navigation system.

State transition array

Wherein,

wherein T is_gx、T_gy、T_gzRespectively representing the components of the gyro drift with respect to time on three axes, T_ax、T_ay、T_azRespectively representing the components of the accelerometer error in three axes with respect to time;

and step B, establishing an INS/ESGM measuring equation.

The electrostatic gyro monitor ESGM monitors and compensates the inertial navigation system INS, but it can only provide position and heading information. INS/ESGM system stateThe equation of state uses INS equation of state such that X₁X. The difference value of the ESGM and the INS position and course information is the observed quantity Z₁：

Wherein L is_EL is divided intoThe latitude, lambda, of ESGM and INS are respectively expressed_Eλ represents longitude of ESGM and INS, h_EH represents the height of ESGM and INS respectively; k_EAnd K respectively represent heading angles of the ESGM and the INS. [ m ] of₁，m₂，m₃，m₄]^TIs zero-mean white gaussian noise generated by the ESGM position and orientation errors.

The corresponding INS/ESGM measurement equation is:

Z₁(t)＝H₁(t)X₁(t)+V₁(t)

wherein, V₁(t)＝[m₁，m₂，m₃，m₄]^TIs zero-mean white gaussian noise generated by the ESGM position and orientation errors. H₁(t) is the measurement array, INS/ESGM State equation X₁(t)＝X(t)。

And step C, determining an INS/GPS measurement equation.

The GPS errors include: SA error, ionosphere transmission error, troposphere transmission error, ephemeris error, and clock synchronization error, etc. In the combined mode of position and speed, the error of GPS is regarded as measuring noise, and the state is not expanded, and X is taken₂(t)＝X(t)。

Observed quantity Z₂The GPS and INS position difference and velocity difference. The INS/GPS system state equation uses the INS state equation: x₂(t) x (t), INS/GPS measurement equation:

Z_{2} = [\begin{matrix} Z_{p} (t) \\ Z_{V} (t) \end{matrix}] = H_{2} (t) X_{2} (t) + v_{2} (t)

wherein the position quantity is measured

Measurement of speed

Measuring array

v₂(t) represents a zero-mean Gaussian white noise matrix generated by the position error and the orientation error of the GPS. L is_GIndicating the latitude, λ, of the GPS_GIndicates the longitude, h, of the GPS_GIndicating the height, V, of the GPS_GEIndicating east velocity, V, in GPS_GNIndicating north velocity, V, in GPS_GUIndicating the speed of the sky in GPS.

And D, determining an INS/DVL system equation.

The Doppler log can provide earth speed and drift angle, and its measurement error mainly includes speed deviation error delta V_dThe scale coefficient error delta C and the drift angle error delta, the speed measurement error and the drift angle measurement error can be described by a first-order Markov process, the scale coefficient error is a random constant, and the error equation is as follows:

wherein, δ V_dDelta and delta C respectively represent the change rates of speed error, drift angle error and scale coefficient error, beta_d、β_ΔThe inverse of the correlation time constant, w, of the velocity error Markov process and the drift angle error Markov process, respectively_d、w_ΔWhite noise indicating a velocity error and a drift angle error. Thus, X is_d＝[δV_d δΔ δC]^TAnd as an expansion system state matrix, forming a system state equation of the INS/GPS together with an INS error equation:

wherein the matrix G_d＝diag[1 1 0]。F_dIs an extended state transfer array, W, of a Doppler velocity measurement system_dIs an amplification process noise array of the Doppler velocity measurement subsystem.

Observed quantity Z₃(t) is the speed difference between the DVL and the INS, and the observation equation of the DVL/INS is as follows:

Z_{3} (t) = [\begin{matrix} V_{E} - V_{dE} \\ V_{N} - V_{dN} \end{matrix}] = H_{3} (t) X_{3} (t) + V_{3} (t)

wherein the measuring array

V₃(t)＝[m_VE m_VN]^TFor white noise observation, K' is the course angle of the ship, V_dE、V_dNIndicating the eastbound and northbound velocities of the doppler measurements.

And step two, discretizing the state equation and the measurement equation to establish a discrete state space model.

The federal filtering method of the invention adjusts the variable proportion weight to make the action weight of the main filter and the sub-filter adapt to the working condition of the navigation equipment and the effective condition of the navigation data, and the main filter provides feedback information for the sub-filter and has good performance. The variable-scale adaptive federal filter structure applying the method of the invention is shown in fig. 3. The inertial navigation system INS supplies state quantity to a main filter, and the electrostatic gyro monitors ESGM and INS generate observed quantity Z₁Sending signals to the first sub-filter, and generating observation quantity Z by GPS and INS of satellite positioning and navigation system₂And sending a signal to the second sub-filter, wherein the Doppler velocity measurement systems DVL and INS generate an observed quantity Z3 and send a signal to the third sub-filter, and the signal filtered by the three sub-filters is sent to the main filter. Performing proportional weight calculation according to historical data of an inertial navigation system INS, an electrostatic gyro monitor ESGM, a satellite positioning navigation system GPS and a Doppler velocity measurement system DVL through an AutoRegressive (AR) model to obtain an information weight theta_iAnd then provided to the main filter as a basis for information distribution ratio calculationAfter one filtering period is completed, the main filter provides global estimated values to the INS and the three sub-filters

Estimating mean square error P_gAnd information distribution ratio beta_iTo reset. And obtaining the optimal estimation through continuous time updating and measurement updating.

The state equation and the measurement equation of the main filter and each subsystem are discretized, and the INS state equation, the INS/ESGM measurement equation, the INS/GPS measurement equation, the INS/DVL state equation and the measurement equation are discretized in the embodiment of the invention. The discrete state space model established after discretization of the corresponding state equation and the measurement equation is as follows:

where k represents the index value of the discrete time state, X_kRepresenting the state quantity at time k, Z_kRepresents the observed quantity at time k, H_kA measuring matrix representing the time k, Γ_k/k-1Is a noise-driven array, W_k-1Is system state noise and process noise covariance matrix Q_k＝E[W_kW_k ^T]，V_kNoise is measured for the system and a noise covariance matrix R is measured_k＝E[V_kV_k ^T]. State transition matrix phi_k/k-1The method is obtained by adopting a step-by-step cumulative discretization method:

And step three, obtaining an information weight value through an autoregressive model according to historical data measured by each navigation sensor system before.

The calculation of the information weight of the navigation sensor system is based on a plurality of historical data measured by the navigation sensor system before, and is obtained by the analysis of an autoregressive model (AR model). The following description will be given taking one of the navigation sensor systems as an example.

Let y (τ) be the output sequence of the navigation sensor system, and the AR model of y (τ) is expressed as:

wherein τ represents the sequence y (τ)) N denotes the order of the AR model, a_kIs that the mean value is zero and the variance is sigma²White Gaussian noise of (a)_kThe error of the model is represented by,

representing AR model coefficients of order N, the power spectral density S of y (tau) can be obtained_y(ω) is:

representing the AR model coefficients of order τ, e^-jωτDenotes the frequency response of the system, ω denotes the angular frequency and j denotes the imaginary unit. y (tau) autocorrelation function array R_N+1Comprises the following steps:

R_{N + 1} = [\begin{matrix} {\hat{r}}_{y} (0) & {\hat{r}}_{y} (1) & {\hat{r}}_{y} (2) & . . . & {\hat{r}}_{y} (N) \\ {\hat{r}}_{y} (1) & {\hat{r}}_{y} (0) & {\hat{r}}_{y} (2) & . . . & {\hat{r}}_{y} (N - 1) \\ {\hat{r}}_{y} (2) & {\hat{r}}_{y} (1) & {\hat{r}}_{y} (0) & . . . & {\hat{r}}_{y} (N - 2) \\ . & . & . & . & . \\ . & . & . & . & . \\ . & . & . & . & . \\ {\hat{r}}_{y} (N) & {\hat{r}}_{y} (N - 1) & {\hat{r}}_{y} (N - 2) & . . . & {\hat{r}}_{y} (0) \end{matrix}]

from the Yule-Waker equation, the autocorrelation function R of y (τ) is known_N+1Can calculate the model coefficient of AR model

Suppose that

As coefficients of AR model

I is equal to 1, 2, … N,

is the error power of the AR model of order N,

in order to be the predicted minimum error power,

according to the properties of the Toeplite matrix, determining the model coefficients by a Levinson-Durbin recursion algorithm:

the order is obtained by the Final Prediction Error (FPE) rule:

wherein xi is the length of the measurement data.

Thus, the predicted value of the sensor output is obtained:

since the AR model is an all-pole model built on a stationary solution space, the prediction error can reflect the smoothness of the sensor output. The prediction error of each sensor is taken as the information weight theta_iThe smaller the information weight, the smaller the prediction error, and the smoother the navigation sensor output. The information weight is as follows:

y_i(τ) is the output sequence of the ith navigation sensor,

the predicted value of the output sequence of the ith navigation sensor is that i is 0, 1, 2 and 3 respectively correspond to an inertial navigation system, an electrostatic gyro monitor, a satellite positioning navigation system and a Doppler velocity measurement system.

And step four, calculating the information distribution proportion according to the information weight and the information conservation law.

The information distribution proportion selection principle is as follows: the higher the accuracy of the subsystem is, the more reliable the output is, and the more the effect of the estimation value of the main system is reduced, the smaller the information distribution ratio is. The main system is a reference system, the embodiment of the invention is an inertial navigation system, and the subsystems are an electrostatic gyro monitor, a satellite positioning navigation system and a Doppler velocity measurement system. Therefore, the information distribution proportion is proportional to the information weight. Meanwhile, according to information conservation:wherein beta is₀Is the primary system weight, beta_iAnd n is the number of all subsystems. Therefore, the final weight of each navigation sensor system is obtained as follows:

and fifthly, carrying out federal filtering according to the information distribution proportion to obtain the global optimal estimation.

Firstly, information distribution is obtained:

X_i，k-1＝X_k-1

using the estimated error variance matrix P of the ith navigation sensor system_k-1By the ratio beta_iAssigned as the reset estimation error variance matrix P_i，k-1And the process noise covariance matrix Q of the ith navigation sensor system_k-1Distribution as reset process noise covariance matrix Q by proportion beta_i，k-1。X_i，k-1Indicating the state quantity of the i-th navigation sensor system, and k indicates the index value of the discrete-time state.

And secondly, updating time:

one-step prediction matrix for state of ith navigation sensor system, P_i，k/k-1Predicting mean square error matrix for the ith navigation sensor system in one step, phi_i，k/k-1For state-transfer arrays of the i-th navigation sensor system, Γ_i，k/k-1For the noise-driven array of the ith navigation sensor system,

is a state quantity X_i，k-1The predicted value of (2).

Then, observation updating is carried out:

K_{i, k} = P_{i, k / k - 1} H_{i, k}^{T} {(H_{i, k} P_{i, k / k - 1} H_{i, k}^{T} + R_{i, k})}^{- 1}

{\hat{X}}_{i, k} = {\hat{X}}_{i, k / k - 1} + K_{i, k} (Z_{i, k} - H_{i, k} {\hat{X}}_{i, k / k - 1})

P_{i, k} = (I - K_{i, k} H_{i, k}) P_{i, k / k - 1} {(I - K_{i, k} H_{i, k})}^{T} + K_{i, k} R_{i, k} K_{i, k}^{T}

K_i，kfilter gain for the ith navigation sensor system, H_i，kIs a measurement array of the ith navigation sensor system, R_i，kIs a system of the ith navigation sensor systemStatistical noise covariance matrix, Z_i，kIs the observed quantity of the ith navigation sensor system. k represents the index value of the discrete time state.

And finally obtaining a global optimization value:

is a global estimate, P_gTo estimate the mean square error, m is the number of sub-filters, and m is 3 in the embodiment of the present invention.

And resetting the filtering value and the estimation error variance matrix by using the global optimal solution. Because the upper bound technology of variance is used when the system noise is distributed, all subsystems are not related, and the augmentation system is not coupled. Although the filtering results of the subsystems are suboptimal, the fusion is resynthesized so that the global estimate is optimal.

And performing simulation experiments through the integrated navigation system in the mathematic analysis software matlab environment. The simulation experiment was 8 hours. The following parameters are set, and the variable-ratio federal filtering method provided by the invention is compared with the INS error result. INS constant drift of 0.1 degree/h, drift Markov process time constant T_g300s, adding table zero bias of 0.1mg, adding table Markov process time constant T_aIs 1000 s. ESGM white noise was 0.04 °/h. The GPS outputs add randomly distributed position and velocity errors. DVL error 0.4 m/s.

The proportional weight change curve is shown in fig. 4. From fig. 4, the proportional weight change under different working conditions of the sensor can be reflected, and when the INS plays a main role, the weight β of the INS₀If the ratio is lower than other ratio values, when the GPS or DVL is used for correction, the ratio weight value is higher than other systems, and the function of the correction system is highlighted. The system error curves are obtained as shown in fig. 5 and fig. 6. Wherein the red solid line represents the error of the underwater multi-sensor combined navigation system of the embodiment of the invention, and the black dotted line represents the INS system error. As shown in fig. 5 and 6, the absolute value of the position error peak of the combined system is 20m, and the absolute value of the velocity error peak is 0.8m/s, and does not diverge. And the INS position error is 500m, the speed error is 5m/s and is continuously increased. It can be seen that the navigation information precision obtained by the federal filtering method provided by the invention is superior to that of a single-sensor system, the system output is kept stable and does not diverge within a certain range, and the weak error accumulation (especially position error) of an inertial navigation system is overcomeAnd the reliability and stability of the system are improved. The position error variance calculated under the condition that the information weight of the navigation sensor is changed continuously, particularly under the condition of large jump of GPS and DVL correction switching is 8.92, and the speed error variance is 0.065, so that the performance of the combined navigation system using the method of the invention is relatively stable. As shown in fig. 7, compared with the integrated navigation system using the conventional federal filtering method, the dynamic performance and the position error amplitude of the integrated navigation system using the variable-ratio adaptive federal filtering method of the present invention have the advantages that the position error amplitude of the integrated navigation system using the method of the present invention changes little, and the method has obvious advantages, and can reflect the situation of the integrated navigation system more truly.

The invention obtains the information weight by analyzing a plurality of forward data of the navigation sensor and adopting time sequence analysis, and then calculates the filter weight, thereby changing the proportion of information distribution and achieving the purposes of improving the precision, reliability and fault tolerance of the system. Compared with other self-adaptive federal filtering, the system has the advantages that the information weight is determined according to the source information of the navigation sensor instead of the estimated value processed by Kalman filtering, so that the determination of the weight is closer to the real situation of the navigation sensor, and the situation of the system is reflected more truly.

Claims

1. A time series analysis-based variable proportion self-adaptive federal filtering method is characterized by comprising the following steps:

2. The variable-proportion self-adaptive federal filtering method based on time series analysis according to claim 1, wherein the integrated navigation system is an underwater multi-sensor integrated navigation system, and is composed of four navigation sensor systems: inertial navigation system INS, static gyro monitor ESGM, satellite positioning navigation system GPS and Doppler velocity measurement system DVL.

3. The time series analysis-based variable-scale adaptive federated filtering method according to claim 1 or 2, wherein the system state equation and the measurement equation of step one include: an INS state equation, an INS/ESGM measurement equation, an INS/GPS measurement equation, an INS/DVL system state equation and a measurement equation;

the INS state equation is as follows:

wherein the state quantity

is the component of the accelerometer error in three axes; process noise array W ═ W_bx w_by w_bz w_gx w_gy w_gz w_ax w_ay w_az]^T，w_bx、w_by、w_bzComponent on three axes of white noise, w, representing the gyro drift random constant_gx、w_gy、w_gzComponent of white noise on three axes, w, representing a gyro-drift Markov process_ax、w_ay、w_azA component of white noise on three axes representing accelerometer error; the three axes are an x axis, a y axis and a z axis; noise driving array

Wherein

Is a directional cosine matrix; state transition array

Wherein,

T_gx、T_gy、T_gzrespectively representing the components of the gyro drift with respect to time on three axes, T_ax、T_ay、T_azRespectively representing the components of the accelerometer error with respect to time in three axes,

wherein f is_EFor acceleration in east directionDegree f_NIs the north acceleration, f_UAcceleration in the direction of the sky, ω_ieIs the rotational angular velocity of the earth, R_MRadius of curvature of meridian, R_NThe curvature radius of the prime circle is;

the INS/ESGM measurement equation is as follows: z₁(t)＝H₁(t)X₁(t)+V₁(t)；

Wherein, V₁(t)＝[m₁，m₂，m₃，m₄]^TIs zero-mean white Gaussian noise, H, generated by the ESGM position and orientation errors₁(t) is a measurement array, the observed quantity Z₁(t) is the difference value of the ESGM, the INS and the course information; INS/ESGM State equation X₁(t)＝X(t)；

The INS/GPS measurement equation is as follows:

wherein the position quantity is measured

Measurement of speedMeasuring array

v₂(t) the observed quantity Z is a white Gaussian noise array of zero mean value generated by the position error and the azimuth error of the GPS₂(t) is the difference between the GPS and INS positions and velocities; l is_GIndicating the latitude, λ, of the GPS_GIndicates the longitude, h, of the GPS_GIndicating the height, V, of the GPS_GEIndicating east velocity, V, in GPS_GNIndicating north velocity, V, in GPS_GURepresenting the speed of the sky in GPS; INS/GPS equation of state X₂(t)＝X(t)；

The INS/DVL system state equation is as follows:

the INS/DVL measurement equation is as follows:

wherein the measuring array

4. The time series analysis-based variable-proportion adaptive federated filtering method according to claim 1 or 2, characterized in that the discrete state space model corresponding to each navigation sensor system in step two is:

where k represents the index value of the discrete time state, X_kRepresenting the state quantity at time k, Z_kRepresents the observed quantity at time k, H_kA measuring matrix representing the time k, Γ_k/k-1Is a noise-driven array, W_k-1Is the system state noise at time k-1 and the process noise covariance matrix Q_k＝E[W_kW_k ^T]，V_kFor the system at time k, measure noise, and measure noise covariance matrix R_k＝E[V_kV_k ^T]D, state transition matrix phi_k/k-1The method is obtained by adopting a step-by-step cumulative discretization method:

5. The time series analysis-based variable-proportion adaptive federated filtering method according to claim 1 or 2, wherein the information weight obtained by the autoregressive model in step three is specifically:

first, an autoregressive model of the navigation sensor system is determined:

where y (τ) is the output sequence of the navigation sensor system, τ represents the τ -th point of the sequence y (τ), N represents the order of the autoregressive model, a_kRepresenting the error of the model as mean zero and variance σ²Is highThe white noise of the white noise is generated,

representing the coefficients of an Nth order autoregressive model;

then, the autoregressive model coefficients were determined by the Levinson-Durbin recursion algorithm:

wherein,

for coefficients of an Nth order autoregressive modelI is equal to 1, 2, … N,

is the coefficient of an autoregressive model of order N-1Is determined by the estimated value of (c),

to predict the minimum error power of the autoregressive model of order N,to predict the minimum error power of the N-1 order autoregressive model,

as an intermediate parameter, the parameter is,

the 1 st element and the i +1 st element of the Nth row in the autocorrelation function matrix of the output sequence y (tau) of the navigation sensor respectively;

the order is obtained by the final prediction error rule:

wherein xi is the length of the measured data,

error for AR model of order NA differential power;

and further obtaining a predicted value of the output sequence of the navigation sensor system:

taking the prediction error of each navigation sensor system as an information weight theta_i：

y_i(τ) is the output sequence of the ith navigation sensor,

6. The time series analysis-based variable-proportion adaptive federal filtering method as claimed in claim 1 or 2, wherein the information distribution proportion in the step four is proportional to the information weight, and the sum of the weights of the navigation sensor systems is equal to

wherein, i is 0, 1, 2, 3 respectively correspond inertial navigation system, static top monitor, satellite positioning navigation system and doppler speed measuring system.

7. The time series analysis-based variable-proportion adaptive federated filtering method according to claim 1 or 2, characterized in that the global optimum estimation of step five is obtained specifically by the following process:

first, information allocation is obtained:

X_i，k-1＝X_k-1

wherein, P_i，k-1An estimation error variance matrix P representing the ith navigation sensor system_k-1By the ratio beta_iWeighted estimation error variance matrix, Q_i，k-1Covariance matrix Q representing process noise for the ith navigation sensor system_k-1By the ratio beta_iCovariance matrix of process noise after weight distribution, X_i，k-1An index value representing a state quantity of the i-th navigation sensor system, k representing a discrete-time state; the i is 0, 1, 2 and 3 respectively corresponding to an inertial navigation system, an electrostatic gyro monitor, a satellite positioning navigation system and a Doppler velocity measurement system;

and then time updating is carried out:

wherein,

is a state quantity X_i，k-1The predicted value of (2);

then, observation updating is carried out:

K_{i, k} = P_{i, k / k - 1} H_{i, k}^{T} {(H_{i, k} P_{i, k / k - 1} H_{i, k}^{T} + R_{i, k})}^{- 1}

{\hat{X}}_{i, k} = {\hat{X}}_{i, k / k - 1} + K_{i, k} (Z_{i, k} - H_{i, k} {\hat{X}}_{i, k / k - 1})

P_{i, k} = (I - K_{i, k} H_{i, k}) P_{i, k / k - 1} {(I - K_{i, k} H_{i, k})}^{T} + K_{i, k} R_{i, k} K_{i, k}^{T}

wherein, K_i，kFilter gain for the ith navigation sensor system, H_i，kIs a measurement array of the ith navigation sensor system, R_i，kIs the system measurement noise covariance matrix, Z, of the ith navigation sensor system_i，kIs the observed quantity of the ith navigation sensor system;

and finally obtaining global optimization estimation:

wherein,

is a global estimate, P_gTo estimate the mean square error, m is the number of filters.