CN105005099B - Atmospheric parameter calculation method based on strapdown inertial navigation and flight control system - Google Patents

Abstract

The invention discloses an atmospheric parameter calculation method based on a strapdown inertial navigation and flight control system. Based on definition of an attack angle and a sideslip angle, a mapping relation of true air speed projection to a machine system is established. Combined with a characteristic of stable high-altitude atmosphere flow speed and by utilization of inertial navigation information, a one-step prediction model of the atmosphere parameter is obtained. According to an aircraft pneumatic model, a function relation of a pneumatic parameter and the atmosphere parameter is established, a mathematical relation of an inertial navigation measuring valve and the atmosphere parameter is constructed, and the atmosphere parameter is solved through a Kalman filtering expansion method. The method provides a redundant means for measurement of the atmosphere parameter, solves measuring problems in hypersonic speed and high maneuver conditions and can provide an atmosphere parameter in real time in a flight envelope.

Description

A kind of atmospheric parameter calculation method based on inertial navigation and flight control system

Technical field

The present invention relates to a kind of atmospheric parameter calculation method, including the resolving of the angle of attack, yaw angle and true air speed, more particularly to A kind of atmospheric parameter calculation method based on inertial navigation and flight control system.

Background technology

Air data system is to complete the important airborne Aerial Electronic Equipment that atmospheric parameter is perceived, measures, resolved and exports, at present Mainly there are traditional air data system, the class of embedded air data system two.

Traditional air data system to stretch out the pitot of body as mark, and with reference to the other sensors (angle of attack/sideslip Angle/total temperature probe) direct measurement of stagnation pressure, static pressure, the angle of attack, yaw angle and stagnation temperature is realized, then calculated using atmosphere data Machine carries out the resolving and correction of correlation, completes the measurement of atmosphere data, and measuring principle is simple, development is earliest, technology maturation is stable, At home and abroad extensively apply on military secret and civil aircraft.But with the continuous development of modern technical aeronautics, traditional atmosphere data sensing Technology is gradually difficult to meet the flight demand of high performance aircraft, and such as under big angles-of-attack state, flow separation of being bullied affects, tradition Air data system is difficult to measure accurate pressure, and the device for stretching out will become the master for causing head vortex and medio-lateral instability Factor is wanted, causes its handling degradation；Under hypersonic flight state, prominent measurement apparatus are difficult in adapt to high temperature Environment, and the highly integrated Design of Aerodynamic Configuration of hypersonic aircraft is had a strong impact on, the measurement apparatus stretched out in addition are also difficult To meet Stealth demand.

Embedded air data system is a kind of by the pressure being embedded on aircraft front end (or wing) diverse location Sensor array obtains atmospheric parameter measuring the pressure distribution of aircraft surface by pressure distribution.The proposition of this technology With development, the General Promotion level of atmosphere data sensing technology.The device is not only convenient for stealthy, and efficiently solves and attacks greatly Atmosphere parameter measurement problem when angle, High Mach number flight, drastically increases the scope of application of air data system.But by pneumatic Time delay and pneumatic conduit frequency response characteristic are limited, and certainty of measurement declines when causing high motor-driven and high high-altitude flight, has a strong impact on The real-time of system.At the same time, during hypersonic flight, collector directly contact high temperature gas flow, easy initiating failure.

Strap-down inertial is a kind of self-aid navigation method that carrier navigation information is obtained according to Newton mechanics law, it Using the inertance element sensitive carrier movable information such as gyroscope, accelerometer, then computing is integrated by computer and is obtained The navigational parameters such as the attitude of carrier, speed and position.Flight control system is the operation according to pilot, to rudder face and is started Machine transmits control information, the system so as to implement control to the configuration of aircraft, flight attitude and kinematic parameter.

Both navigation modes respectively have feature：Strap-down inertial abundant information, navigation accuracy is higher；Flight control system Abundant information, accuracy is relatively low, but with atmospheric parameter close relation.Therefore, using strap-down inertial and flight control The information of system, estimates atmospheric parameter, can be load in full flight envelope on the premise of hardware device is not increased Body provides high-precision atmospheric parameter, with prominent using value.

The content of the invention

Goal of the invention：Atmospheric parameter asking of being difficult to measure in the case of hypersonic, high maneuver in order to overcome prior art Topic, provides in real time atmospheric parameter in flight envelope, there is provided a kind of atmospheric parameter measurement that need not additionally increase hardware device Redundant means；The method is applied to inertial navigation/flight control system, based on the angle of attack, the definition of yaw angle, sets up true air speed throwing Shadow, with reference to the speed stable feature of upper atmosphere flowing, using inertial navigation information atmospheric parameter is obtained to the mapping relations of body system One-step prediction model, according to flight vehicle aerodynamic model, set up the functional relation of aerodynamic parameter and atmospheric parameter, build inertial navigation Measuring value and the mathematical relationship of atmospheric parameter, are solved by the method for EKF to atmospheric parameter, it is to avoid Direct measurement to atmospheric parameter, solves the problems, such as the measurement of atmospheric parameter in the case of hypersonic, high maneuver.

Technical scheme：For achieving the above object, the technical solution used in the present invention is：

A kind of atmospheric parameter calculation method based on inertial navigation and flight control system, comprises the following steps：

Step 1, initial information are arranged：Quantity of state is angle of attack, sideslip angle beta and true air speed V_t, arrange angle of attack initial value be Carrier pitching angle theta, the initial value of sideslip angle beta is zero, and the angle of attack, sideslip angular unit are rad, true air speed V_tInitial value be bearer rate V_b, unit is ft/s；System noise variance matrix Q, measuring noise square difference battle array R and initial estimation error covariance matrix P are set；

Step 2, quantity of state rate of change are resolved：The mapping relations of body system are projected to according to true air speed, quantity of state is calculated and is existed The rate of change at current k moment, concrete steps include：

Step 201, the body system acceleration for being loaded into the output of k moment SINS, including heading accelerationFuselage right flank accelerationWith fuselage vertical accelerationUnit of acceleration is ft/s²(per square of second of foot)；

Step 202, the mapping relations that body system is projected to based on true air speed, carry out differential, calculate the rate of change of quantity of state：

In formula,WithFor the rate of change of k moment quantity of states；α_k、β_kAnd V_t,kFor the estimate of k moment quantity of states, The estimate of initial time quantity of state takes initial value；

Step 3, quantity of state one-step prediction：According to quantity of state rate of change, with reference to sampling step length T, a step is carried out to quantity of state Prediction：

In formula,WithFor the one-step prediction of k moment quantity of states；

Step 4, one-step prediction mean square error are resolved：Calculating state Matrix of shifting of a step F and discretization, with reference to the k moment Evaluated error covariance matrix P_k, resolve one-step prediction mean square error P_k+1,k, concrete steps include：

Step 401, the estimate according to k moment quantity of states and body system acceleration, calculate state Matrix of shifting of a step F, F In element be：

F_3,3=0；

Step 402, to F discretizations, obtain Φ_k+1,k=I+FT, I are unit matrix；

Step 403, with reference to the evaluated error covariance matrix P at k moment_k, resolve one-step prediction mean square error P_k+1,k：

In formula, P_kFor the evaluated error covariance matrix at k moment, the evaluated error covariance matrix of initial time is P；

Step 5, filtering gain are resolved：According to the information that SINS and flight control system are provided, set up pneumatic Parameter and the functional relation of atmospheric parameter, resolve measurement matrix, and concrete steps include：

Step 501, the body system angular speed of the SINS at loading subsequent time k+1 moment output, flight control The flight controlled quentity controlled variable of system and atmospheric density ρ, the unit of atmospheric density is slug/ft³(per cubic feet of slug)；

Step 502, according to aerodynamic model, with reference to the information that SINS and flight control system are provided, will be pneumatic Coefficient is rewritten as the function of atmospheric parameter, resolves measurement matrix H, and element is in H：

In formula, d_iFor the constant coefficient in function, by the function coefficients of Aerodynamic Coefficient, the speed of carrier, angular speed, attitude And rotary inertia is resolved and obtained, i=1～43；

Step 503, resolving filtering gain

Step 6, estimation mean square error are resolved：According to filtering gain K_k+1, measurement matrix H, one-step prediction mean square error P_k+1,k With measurement noise variance matrix R, resolve and estimate mean square error P_k+1：

Step 7, state estimation：The one-step prediction of bonding state amount, according to the aerodynamic model of carrier, resolves the k+1 moment Observed quantity estimate, so as to realize the estimation of the estimation of quantity of state, i.e. atmospheric parameter, concrete steps include：

Step 701, the body system angular acceleration for being loaded into the output of k+1 moment SINS, including heading angle adds SpeedFuselage right flank angular accelerationWith the vertical angular acceleration of fuselageAngular acceleration unit is rad/s²；It is loaded into k The flight controlled quentity controlled variable at+1 moment, body system acceleration

Step 702, by quantity of state the k moment one-step prediction, the flight controlled quentity controlled variable at k+1 moment, the body system at k+1 moment The body system angular acceleration at acceleration and k+1 moment is updated in aerodynamic model, resolves observed quantity estimate：

In formula, d₄₄For the constant in function, by the function coefficients of Aerodynamic Coefficient, the speed of carrier, angular speed, attitude and Rotary inertia is resolved and obtained；

Step 703, the observed quantity at note k+1 moment areDuring k+1 The observed quantity estimate at quarter isQuantity of state is pre- in a step at k moment Survey and beOne-step prediction of the quantity of state at the k+1 moment beThen There is state estimationAfter state estimation is obtained, return to step 2, after proceeding Continuous atmospheric parameter is resolved.

Beneficial effect：The atmospheric parameter calculation method based on inertial navigation and flight control system that the present invention is provided, with Prior art is compared, with following advantage：1st, existing binding kineticses equation carries out atmospheric parameter calculation method with winged control parameter In directly using flight control system provide power, moment information, true air speed is carried out in body system according to rigid body kinematics principle The resolving of lower three axle components, and then flow angle is resolved, but do not consider relation of power, the torque with atmospheric parameter itself, affect to calculate Method precision, for this problem, the present invention sets up the functional relation of aerodynamic parameter and atmospheric parameter according to carrier pneumatic model, Improve the accuracy of measurement equation, it is ensured that arithmetic accuracy；2nd, for existing air data system in motor-driven, the high ultrasound of height Speed, be difficult under super high altitude flight state to measure, easy failure, the problem of poor real, by inertial navigation and flight control system The mode of information fusion, while system autonomy is guaranteed, there is provided the Real-Time Atmospheric parameter in full flight envelope.

Description of the drawings

Fig. 1 is the algorithm flow chart of the atmospheric parameter calculation method of the present invention；

Fig. 2 is the simulation program structure figure of the atmospheric parameter calculation method of the present invention；

Fig. 3 is the result that the atmospheric parameter of the atmospheric parameter calculation method of the present invention is resolved, wherein：3 (a) is angle of attack contrast Curve, 3 (b) is yaw angle correlation curve, and 3 (c) is true air speed correlation curve.

Specific embodiment

The present invention is further described below in conjunction with the accompanying drawings.

A kind of atmospheric parameter calculation method based on inertial navigation and flight control system is illustrated in figure 1, including it is following Step：

Step 1, initial information are arranged：Quantity of state is angle of attack, sideslip angle beta and true air speed V_t, arrange angle of attack initial value be Carrier pitching angle theta, the initial value of sideslip angle beta is zero, and the angle of attack, sideslip angular unit are rad, true air speed V_tInitial value be bearer rate V_b, unit is ft/s；System noise variance matrix Q, measuring noise square difference battle array R and initial estimation error covariance matrix P are set；

Step 2, quantity of state rate of change are resolved：The mapping relations of body system are projected to according to true air speed, quantity of state is calculated and is existed The rate of change at current k moment, concrete steps include：

Step 201, the body system acceleration for being loaded into the output of k moment SINS, including heading acceleration Fuselage right flank accelerationWith fuselage vertical accelerationUnit of acceleration is ft/s²(per square of second of foot)；

Step 202, the mapping relations that body system is projected to based on true air speed, carry out differential, calculate the rate of change of quantity of state：

In formula,WithFor the rate of change of k moment quantity of states；α_k、β_kAnd V_t,kFor the estimate of k moment quantity of states, The estimate of initial time quantity of state takes initial value；

Step 3, quantity of state one-step prediction：According to quantity of state rate of change, with reference to sampling step length T, a step is carried out to quantity of state Prediction：

In formula,WithFor the one-step prediction of k moment quantity of states；

Step 4, one-step prediction mean square error are resolved：Calculating state Matrix of shifting of a step F and discretization, with reference to the k moment Evaluated error covariance matrix P_k, resolve one-step prediction mean square error P_k+1,k, concrete steps include：

Step 401, the estimate according to k moment quantity of states and body system acceleration, calculate state Matrix of shifting of a step F, F In element be：

F_3,3=0；

Step 402, to F discretizations, obtain Φ_k+1,k=I+FT, I are unit matrix；

Step 403, with reference to the evaluated error covariance matrix P at k moment_k, resolve one-step prediction mean square error P_k+1,k：

In formula, P_kFor the evaluated error covariance matrix at k moment, the evaluated error covariance matrix of initial time is P；

Step 5, filtering gain are resolved：According to the information that SINS and flight control system are provided, set up pneumatic Parameter and the functional relation of atmospheric parameter, resolve measurement matrix, and concrete steps include：

Step 501, the body system angular speed of the SINS at loading subsequent time k+1 moment output, flight control The flight controlled quentity controlled variable of system and atmospheric density ρ, the unit of atmospheric density is slug/ft³(per cubic feet of slug)；

Step 502, according to aerodynamic model, with reference to the information that SINS and flight control system are provided, will be pneumatic Coefficient is rewritten as the function of atmospheric parameter, resolves measurement matrix H, and element is in H：

In formula, d_iFor the constant coefficient in function, by the function coefficients of Aerodynamic Coefficient, the speed of carrier, angular speed, attitude And rotary inertia is resolved and obtained, i=1～43；

Step 503, resolving filtering gain

Step 6, estimation mean square error are resolved：According to filtering gain K_k+1, measurement matrix H, one-step prediction mean square error P_k+1,k With measurement noise variance matrix R, resolve and estimate mean square error P_k+1：

Step 7, state estimation：The one-step prediction of bonding state amount, according to the aerodynamic model of carrier, resolves the k+1 moment Observed quantity estimate, so as to realize the estimation of the estimation of quantity of state, i.e. atmospheric parameter, concrete steps include：

Step 701, the body system angular acceleration for being loaded into the output of k+1 moment SINS, including heading angle adds SpeedFuselage right flank angular accelerationWith the vertical angular acceleration of fuselageAngular acceleration unit is rad/s²；It is loaded into k The flight controlled quentity controlled variable at+1 moment, body system acceleration

Step 702, by quantity of state the k moment one-step prediction, the flight controlled quentity controlled variable at k+1 moment, the body system at k+1 moment The body system angular acceleration at acceleration and k+1 moment is updated in aerodynamic model, resolves observed quantity estimate：

In formula, d₄₄For the constant in function, by the function coefficients of Aerodynamic Coefficient, the speed of carrier, angular speed, attitude and Rotary inertia is resolved and obtained；

Step 703, the observed quantity at note k+1 moment areDuring k+1 The observed quantity estimate at quarter isQuantity of state is pre- in a step at k moment Survey and beOne-step prediction of the quantity of state at the k+1 moment be Then there is state estimationObtain state estimation after, return to step 2, continue into The follow-up atmospheric parameter of row is resolved.

In order to evaluate the performance of true air speed calculation method proposed by the present invention, simulated program is devised, in matlab platforms On algorithm is verified, structure is as shown in Fig. 2 the simulated program is comprised the following steps：

(1) flight path is set, according to carrier pneumatic model, aircraft trim is carried out；

(2) flight simulation is carried out, flight controlled quentity controlled variable, inertial guidance data, atmospheric parameter with track matching is generated；

(3) to flight controlled quentity controlled variable, aerodynamic model, inertial guidance data Injection Error；

(4) carry out atmospheric parameter resolving using the flight controlled quentity controlled variable after Injection Error, aerodynamic model, inertial guidance data, with fly The actual value that row emulation is obtained is analyzed；

Experimental result is as shown in figure 3, show that atmospheric parameter calculation result of the invention is essentially coincided with actual value, the angle of attack is missed Difference average is 0.16268 °, and yaw angle error mean is 0.05505 °, and true air speed error mean is -4.82748m/s, it was demonstrated that The present invention resolves the correctness and validity of atmospheric parameter method.

The above is only the preferred embodiment of the present invention, it should be pointed out that：For the ordinary skill people of the art For member, under the premise without departing from the principles of the invention, some improvements and modifications can also be made, these improvements and modifications also should It is considered as protection scope of the present invention.

Claims (1)

1. a kind of atmospheric parameter calculation method based on inertial navigation and flight control system, it is characterised in that：Including following step Suddenly：

Step 1, initial information are arranged：Quantity of state is angle of attack, sideslip angle beta and true air speed V_t, arrange angle of attack initial value bow for carrier Elevation angle theta, the initial value of sideslip angle beta is zero, and the angle of attack, sideslip angular unit are rad, true air speed V_tInitial value be bearer rate V_b, unit For ft/s；System noise variance matrix Q, measurement noise variance matrix R and initial estimation error covariance matrix P are set；

Step 2, quantity of state rate of change are resolved：The mapping relations of body system are projected to according to true air speed, quantity of state is calculated in current k The rate of change at moment, concrete steps include：

Step 201, the body system acceleration for being loaded into the output of k moment SINS, including heading accelerationFuselage Right flank accelerationWith fuselage vertical accelerationUnit of acceleration is ft/s²；

Step 202, the mapping relations that body system is projected to based on true air speed, carry out differential, calculate the rate of change of quantity of state：

$\stackrel{\·}{\mathrm{\α}}=\frac{{\mathrm{cos\α}}_k{\stackrel{\·}{W}}_k-{\mathrm{sin\α}}_k{\stackrel{\·}{U}}_k}{{\mathrm{cos\β}}_k{V}_{t,k}}$

$\stackrel{\·}{\mathrm{\β}}=\frac{{\mathrm{cos\β}}_k{\stackrel{\·}{V}}_k-{\mathrm{sin\β}}_k{\mathrm{cos\α}}_k{\stackrel{\·}{U}}_k-{\mathrm{sin\β}}_k{\mathrm{sin\α}}_k{\stackrel{\·}{W}}_k}{V_{t,k}}$

${\stackrel{\·}{V}}_t={\mathrm{cos\α}}_k{\mathrm{cos\β}}_k{\stackrel{\·}{U}}_k+{\mathrm{sin\β}}_k{\stackrel{\·}{V}}_k+{\mathrm{cos\β}}_k{\mathrm{sin\α}}_k{\stackrel{\·}{W}}_k$

In formula,WithFor the rate of change of k moment quantity of states；α_k、β_kAnd V_t,kFor the estimate of k moment quantity of states, when initial The estimate for carving quantity of state takes initial value；

Step 3, quantity of state one-step prediction：According to quantity of state rate of change, with reference to sampling step length T, one-step prediction is carried out to quantity of state：

${\widetilde{\mathrm{\α}}}_{k+1}={\mathrm{\α}}_k+\stackrel{\·}{\mathrm{\α}}T$

${\widetilde{\mathrm{\β}}}_{k+1}={\mathrm{\β}}_k+\stackrel{\·}{\mathrm{\β}}T$

${\tilde{V}}_{t,k+1}=V_{t,k}+{\stackrel{\·}{V}}_{t}T$

In formula,WithFor the one-step prediction of k moment quantity of states；

Step 4, one-step prediction mean square error are resolved：Calculating state Matrix of shifting of a step F and discretization, with reference to the estimation at k moment Error covariance matrix P_k, resolve one-step prediction mean square error P_k+1,k, concrete steps include：

Step 401, the estimate according to k moment quantity of states and body system acceleration, calculate state Matrix of shifting of a step F, in F Element is：

$F_{1,1}=-\frac{{\mathrm{sin\α}}_k{\stackrel{\·}{W}}_k+{\mathrm{cos\α}}_k{\stackrel{\·}{U}}_k}{V_{t,k}{\mathrm{cos\β}}_k}$

$F_{1,2}=\frac{{\mathrm{cos\α}}_k{\stackrel{\·}{W}}_k-{\mathrm{sin\α}}_k{\stackrel{\·}{U}}_k}{V_{t,k}{\cos }^2{\mathrm{\β}}_k}$

$F_{1,3}=\frac{{\mathrm{sin\α}}_k{\stackrel{\·}{U}}_k-{\mathrm{cos\α}}_k{\stackrel{\·}{W}}_k}{V_{t,k}^2{\mathrm{cos\β}}_k}$

$F_{2,1}=\frac{{\mathrm{sin\α}}_k{\mathrm{sin\β}}_k{\stackrel{\·}{U}}_k-{\mathrm{sin\β}}_k{\mathrm{cos\α}}_k{\stackrel{\·}{W}}_k}{V_{t,k}}$

$F_{2,2}=-\frac{{\mathrm{cos\α}}_k{\mathrm{cos\β}}_k{\stackrel{\·}{U}}_k+{\mathrm{cos\β}}_k{\mathrm{sin\α}}_k{\stackrel{\·}{W}}_k+{\mathrm{sin\β}}_k{\stackrel{\·}{V}}_k}{V_{t,k}}$

$F_{2,3}=\frac{{\mathrm{cos\α}}_k{\mathrm{sin\β}}_k{\stackrel{\·}{U}}_k+{\mathrm{sin\β}}_k{\mathrm{sin\α}}_k{\stackrel{\·}{W}}_k-{\mathrm{cos\β}}_k{\stackrel{\·}{V}}_k}{V_{t,k}^2}$

$F_{3,1}={\mathrm{cos\α}}_k{\mathrm{cos\β}}_k{\stackrel{\·}{W}}_k-{\mathrm{sin\α}}_k{\mathrm{cos\β}}_k{\stackrel{\·}{U}}_k$

$F_{3,2}={\mathrm{cos\β}}_k{\stackrel{\·}{V}}_k-{\mathrm{cos\α}}_k{\mathrm{sin\β}}_k{\stackrel{\·}{U}}_k-{\mathrm{sin\β}}_k{\mathrm{sin\α}}_k{\stackrel{\·}{W}}_k$

F_3,3=0

Step 402, to F discretizations, obtain Φ_k+1,k=I+FT, I are unit matrix；

Step 403, with reference to the evaluated error covariance matrix P at k moment_k, resolve one-step prediction mean square error P_k+1,k：

$P_{k+1,k}={\mathrm{\Φ}}_{k+1,k}{P}_k{\mathrm{\Φ}}_{k+1,k}^T+Q$

In formula, P_kFor the evaluated error covariance matrix at k moment, the evaluated error covariance matrix of initial time is P；

Step 5, filtering gain are resolved：According to the information that SINS and flight control system are provided, aerodynamic parameter is set up With the functional relation of atmospheric parameter, measurement matrix is resolved, concrete steps include：

Step 501, body system angular speed, the flight control system of the output of the SINS at loading subsequent time k+1 moment Flight controlled quentity controlled variable and atmospheric density ρ, the unit of atmospheric density is slug/ft³；

Step 502, according to aerodynamic model, with reference to the information that SINS and flight control system are provided, by Aerodynamic Coefficient The function of atmospheric parameter is rewritten as, measurement matrix H is resolved, element is in H：

$H_{1,1}=d_1{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_3{\tilde{V}}_{t,k+1}^2+d_4{\tilde{V}}_{t,k+1}$

$H_{1,2}=d_1{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_2{\tilde{V}}_{t,k+1}^2$

$H_{1,3}=2d_1{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_2{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_3{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_4{\widetilde{\mathrm{\α}}}_{k+1}+d_5+2d_6{\tilde{V}}_{t,k+1}$

$H_{2,1}=d_8{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{10}{\tilde{V}}_{t,k+1}^2+d_{11}{\tilde{V}}_{t,k+1}$

$H_{2,2}=d_8{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_9{\tilde{V}}_{t,k+1}^2$

$H_{2,3}=2d_8{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_9{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_{10}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_{11}{\widetilde{\mathrm{\α}}}_{k+1}+d_{12}+2d_{13}{\tilde{V}}_{t,k+1}$

$H_{3,1}=d_{15}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{17}{\tilde{V}}_{t,k+1}^2+d_{18}{\tilde{V}}_{t,k+1}$

$H_{3,2}=d_{15}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{16}{\tilde{V}}_{t,k+1}^2$

$H_{3,3}=2d_{15}{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_{16}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_{17}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_{18}{\widetilde{\mathrm{\α}}}_{k+1}+d_{19}+2d_{20}{\tilde{V}}_{t,k+1}$

$H_{4,1}=d_{22}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{24}{\tilde{V}}_{t,k+1}^2+d_{25}{\tilde{V}}_{t,k+1}$

$H_{4,2}=d_{22}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{23}{\tilde{V}}_{t,k+1}^2$

$H_{4,3}=2d_{22}{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_{23}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_{24}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_{25}{\widetilde{\mathrm{\α}}}_{k+1}+d_{26}+2d_{27}{\tilde{V}}_{t,k+1}$

$H_{5,1}=2d_{29}{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{V}_{t,k+1}^2+2d_{30}{\widetilde{\mathrm{\α}}}_{k+1}{V}_{t,k+1}^2+d_{31}{\widetilde{\mathrm{\β}}}_{k+1}{V}_{t,k+1}^2+d_{33}{V}_{t,k+1}^2+d_{34}{V}_{t,k+1}$

$H_{5,2}=d_{29}{\widetilde{\mathrm{\α}}}_{k+1}^2{\tilde{V}}_{t,k+1}^2+d_{31}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{32}{\tilde{V}}_{t,k+1}^2$

$\begin{array}{c}{H}_{5,3}=2d_{29}{\widetilde{\mathrm{\α}}}_{k+1}^2{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_{30}{\widetilde{\mathrm{\α}}}_{k+1}^2{\tilde{V}}_{t,k+1}+2d_{31}{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_{32}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}\\ +2d_{33}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_{34}{\widetilde{\mathrm{\α}}}_{k+1}+d_{35}+2d_{36}{\tilde{V}}_{t,k+1}\end{array}$

$H_{6,1}=d_{38}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{40}{\tilde{V}}_{t,k+1}^2+d_{41}{\tilde{V}}_{t,k+1}$

$H_{6,2}=d_{38}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{39}{\tilde{V}}_{t,k+1}^2$

$H_{6,3}=2d_{38}{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_{39}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}+2d_{40}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_{41}{\widetilde{\mathrm{\α}}}_{k+1}+d_{42}+2d_{43}{\tilde{V}}_{t,k+1}$

In formula, d_iFor the constant coefficient in function, by the function coefficients of Aerodynamic Coefficient, the speed of carrier, angular speed, attitude and Rotary inertia is resolved and obtained, i=1～43；

Step 503, resolving filtering gain

Step 6, estimation mean square error are resolved：According to filtering gain K_k+1, measurement matrix H, one-step prediction mean square error P_k+1,kAnd survey Amount noise variance matrix R, resolves and estimates mean square error P_k+1：

$P_{k+1}=(I-K_{k+1}H)P_{k+1,k}{(I-K_{k+1}H)}^T+K_{k+1}{\mathrm{RK}}_{k+1}^T;$

Step 7, state estimation：The one-step prediction of bonding state amount, according to the aerodynamic model of carrier, resolves the observation at k+1 moment Amount estimate, so as to realize the estimation of the estimation of quantity of state, i.e. atmospheric parameter, concrete steps include：

Step 701, the body system angular acceleration for being loaded into the output of k+1 moment SINS, including heading angular accelerationFuselage right flank angular accelerationWith the vertical angular acceleration of fuselageAngular acceleration unit is rad/s²；When being loaded into k+1 The flight controlled quentity controlled variable at quarter, body system acceleration

Step 702, quantity of state is accelerated in the one-step prediction at k moment, the flight controlled quentity controlled variable at k+1 moment, the body system at k+1 moment The body system angular acceleration at degree and k+1 moment is updated in aerodynamic model, resolves observed quantity estimate：

${\widetilde{\stackrel{\·}{U}}}_{k+1}=d_1{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_2{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_3{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_4{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_5{\tilde{V}}_{t,k+1}+d_6{\tilde{V}}_{t,k+1}^2+d_7$

${\widetilde{\stackrel{\·}{V}}}_{k+1}=d_8{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_9{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{10}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{11}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_{12}{\tilde{V}}_{t,k+1}+d_{13}{\tilde{V}}_{t,k+1}^2+d_{14}$

${\widetilde{\stackrel{\·}{W}}}_{k+1}=d_{15}{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{16}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{17}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{18}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_{19}{\tilde{V}}_{t,k+1}+d_{20}{\tilde{V}}_{t,k+1}^2+d_{21}$

${\widetilde{\stackrel{\·}{P}}}_{k+1}=d_{22}{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{23}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{24}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{25}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_{26}{\tilde{V}}_{t,k+1}+d_{27}{\tilde{V}}_{t,k+1}^2+d_{28}$

$\begin{array}{c}{\widetilde{\stackrel{\·}{Q}}}_{k+1}=d_{29}{\widetilde{\mathrm{\α}}}_{k+1}^2{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{30}{\widetilde{\mathrm{\α}}}_{k+1}^2{\tilde{V}}_{t,k+1}^2+d_{31}{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{32}{\widetilde{\mathrm{\β}}}_{k+1}{V}_{t,k+1}^2\\ +d_{33}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{34}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_{35}{\tilde{V}}_{t,k+1}+d_{36}{\tilde{V}}_{t,k+1}^2+d_{37}\end{array}$

${\widetilde{\stackrel{\·}{R}}}_{k+1}=d_{38}{\widetilde{\mathrm{\α}}}_{k+1}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{39}{\widetilde{\mathrm{\β}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{40}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}^2+d_{41}{\widetilde{\mathrm{\α}}}_{k+1}{\tilde{V}}_{t,k+1}+d_{42}{\tilde{V}}_{t,k+1}+d_{43}{\tilde{V}}_{t,k+1}^2+d_{44}$

In formula, d₄₄For the constant in function, by the function coefficients of Aerodynamic Coefficient, the speed of carrier, angular speed, attitude and rotation Inertia is resolved and obtained；

Step 703, the observed quantity at note k+1 moment areThe k+1 moment Observed quantity estimate isOne-step prediction of the quantity of state at the k moment beOne-step prediction of the quantity of state at the k+1 moment beThen There is state estimationAfter state estimation is obtained, return to step 2, after proceeding Continuous atmospheric parameter is resolved.